mirror of https://github.com/CGAL/cgal
interval_support, is now ported from EXACUS2CGAL
This commit is contained in:
Ralf Schindlbeck
parent
7eee16103a
commit
edd70f970a
|
|
@ -2917,6 +2917,7 @@ Nef_S2/test/Nef_S2/CMakeLists.txt -text
|
|||
Number_types/doc_tex/NumberTypeSupport_ref/CORE_BigFloat.tex -text
|
||||
Number_types/doc_tex/NumberTypeSupport_ref/fundamental_types.tex -text
|
||||
Number_types/doc_tex/NumberTypeSupport_ref/open.tex -text
|
||||
Number_types/include/CGAL/Float_traits.h -text
|
||||
Number_types/include/CGAL/Number_types/core_interval_support.h -text
|
||||
Number_types/include/CGAL/Sqrt_extension/Algebraic_extension_traits.h -text
|
||||
Number_types/include/CGAL/Sqrt_extension/Algebraic_structure_traits.h -text
|
||||
|
|
@ -2927,6 +2928,9 @@ Number_types/include/CGAL/Sqrt_extension/Real_embeddable_traits.h -text
|
|||
Number_types/include/CGAL/Sqrt_extension/Scalar_factor_traits.h -text
|
||||
Number_types/include/CGAL/Sqrt_extension/Sqrt_extension_type.h -text
|
||||
Number_types/include/CGAL/Sqrt_extension/io.h -text
|
||||
Number_types/include/CGAL/core_interval_support.h -text
|
||||
Number_types/include/CGAL/interval_support.h -text
|
||||
Number_types/include/CGAL/leda_interval_support.h -text
|
||||
Number_types/test/Number_types/CMakeLists.txt -text
|
||||
Optimisation_doc/doc_tex/Bounding_volumes/annulus.eps -text svneol=unset#application/postscript
|
||||
Optimisation_doc/doc_tex/Bounding_volumes/annulus.gif -text svneol=unset#image/gif
|
||||
|
|
|
|||
|
|
@ -0,0 +1,477 @@
|
|||
// TODO: Add licence
|
||||
//
|
||||
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
|
||||
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
|
||||
//
|
||||
// $URL:$
|
||||
// $Id: 4515 rschindl$
|
||||
//
|
||||
//
|
||||
// Author(s) : Ralf Schindlbeck <rschindl@mpi-inf.mpg.de>
|
||||
//
|
||||
// ============================================================================
|
||||
|
||||
// TODO: The comments are all original EXACUS comments and aren't adapted. So
|
||||
// they may be wrong now.
|
||||
|
||||
|
||||
/*! \file CGAL/Float_traits.h
|
||||
\brief Defines class CGAL::Float_traits.
|
||||
|
||||
Provides dependent types and function objects for all the functions
|
||||
beyond operators of the different float types.
|
||||
*/
|
||||
#ifndef CGAL_FLOAT_TRAITS_H
|
||||
#define CGAL_FLOAT_TRAITS_H
|
||||
|
||||
#include <CGAL/basic.h>
|
||||
|
||||
// #include <NiX/enums.h>
|
||||
|
||||
#include <stack>
|
||||
|
||||
// namespace NiX {
|
||||
CGAL_BEGIN_NAMESPACE
|
||||
|
||||
//compensation for the #include <NiX/enums.h>
|
||||
|
||||
//Enumerator indicating the desired rounding mode for an operation
|
||||
//and the global rounding mode, if exists.
|
||||
enum Rounding_mode {
|
||||
TO_NEAREST = 0, //!< round to nearest number
|
||||
TO_ZERO = 1, //!< round towards zero
|
||||
TO_P_INF = 2, //!< round up (towards positive infinity)
|
||||
TO_N_INF = 3, //!< round down (towards negative infinity)
|
||||
TO_INF = 4 //!< round away from zero
|
||||
};
|
||||
//Enumerator indicating a special value represented by a floating point value.
|
||||
enum Special_value {
|
||||
SV_NONE = 0, //!< indicates non-specialness of a value
|
||||
SV_N_ZERO = -1, //!< zero with a negative sign
|
||||
SV_P_ZERO = 1, //!< zero with a positive sign
|
||||
SV_N_INF = -2, //!< negative infinity
|
||||
SV_P_INF = 2, //!< positive infinity
|
||||
SV_NAN = 3 //!< "Not Any Number (not even infinity)": error indicator
|
||||
};
|
||||
|
||||
typedef long Precision_type;
|
||||
|
||||
// The tags for Float type corresponding to the float type concepts
|
||||
// ================================================================
|
||||
|
||||
//! corresponds to the \c FloatNumber concept.
|
||||
struct Float_number_tag {};
|
||||
|
||||
//! corresponds to the \c ExactFloatNumber concept.
|
||||
struct Exact_float_number_tag {};
|
||||
|
||||
//! corresponds to the \c FixedPointNumber concept.
|
||||
struct Fixed_point_number_tag {};
|
||||
|
||||
// The float traits kernel class
|
||||
// =========================================================================
|
||||
|
||||
// This class provides implementation dependent access to the precision and
|
||||
// rounding mode settings along with information about the float type, the
|
||||
// corresponding integer type and the float concept.
|
||||
|
||||
template< class FT_ >
|
||||
class Float_traits_kernel {
|
||||
public:
|
||||
typedef FT_ FT;
|
||||
typedef CGAL::Null_tag IT;
|
||||
typedef CGAL::Null_tag Float_type;
|
||||
|
||||
typedef CGAL::Null_tag Precision_type;
|
||||
typedef CGAL::Null_functor Get_precision;
|
||||
typedef CGAL::Null_functor Set_precision;
|
||||
typedef CGAL::Null_functor Get_rounding_mode;
|
||||
typedef CGAL::Null_functor Set_rounding_mode;
|
||||
};
|
||||
|
||||
// The float traits template
|
||||
// =========================
|
||||
|
||||
// The general float traits is only declared and not defined.
|
||||
// That leads to an immediate error message at the first attempt to use
|
||||
// a non-supported float type.
|
||||
// Documentation for the concept of being a valid specialization
|
||||
// of Float_traits is provided in NumeriX/manual/AlgebraicKernel/AK_Float.dox
|
||||
|
||||
template< class FT_ >
|
||||
class Float_traits;
|
||||
|
||||
|
||||
// The float traits base class
|
||||
// =========================================================================
|
||||
|
||||
|
||||
//! The template specialization that can be used for types that are not any
|
||||
//! of the float type concepts. All functors are set to \c CGAL::Null_functor.
|
||||
//! See also \link NiX_Float_traits_functors concept Float_traits \endlink .
|
||||
//! \ingroup NiX_Float_traits_bases
|
||||
//
|
||||
template< class FT_kernel_ >
|
||||
class Float_traits_base {
|
||||
private:
|
||||
typedef FT_kernel_ FT_kernel;
|
||||
|
||||
public:
|
||||
typedef CGAL::Null_functor Construct;
|
||||
typedef CGAL::Null_functor Get_mantissa;
|
||||
typedef CGAL::Null_functor Get_exponent;
|
||||
typedef CGAL::Null_functor Get_mantissa_and_exponent;
|
||||
typedef CGAL::Null_functor Round_relative;
|
||||
typedef CGAL::Null_functor Round_absolute;
|
||||
|
||||
typedef CGAL::Null_functor Add;
|
||||
typedef CGAL::Null_functor Sub;
|
||||
typedef CGAL::Null_functor Mul;
|
||||
typedef CGAL::Null_functor Div;
|
||||
typedef CGAL::Null_functor Sqrt;
|
||||
typedef CGAL::Null_functor Mul_by_pow_of_2;
|
||||
typedef CGAL::Null_functor Div_by_pow_of_2;
|
||||
|
||||
typedef CGAL::Null_functor Floor;
|
||||
typedef CGAL::Null_functor Ceil;
|
||||
|
||||
typedef CGAL::Null_functor Construct_special_value;
|
||||
|
||||
typedef CGAL::Null_functor Is_special_value;
|
||||
typedef CGAL::Null_functor Get_special_value;
|
||||
typedef CGAL::Null_functor Is_zero;
|
||||
typedef CGAL::Null_functor Is_inf;
|
||||
|
||||
typedef CGAL::Null_functor Set_exponent;
|
||||
|
||||
class Rounding_mode_controller {
|
||||
public:
|
||||
|
||||
inline Rounding_mode_controller( const Rounding_mode& mode ) {
|
||||
// store current rounding mode
|
||||
rounding_modes.push( get_rounding_mode() );
|
||||
|
||||
#ifndef NDEBUG
|
||||
// store scope id
|
||||
scope_id = rounding_modes.size();
|
||||
#endif
|
||||
set_rounding_mode( mode );
|
||||
};
|
||||
|
||||
inline ~Rounding_mode_controller() {
|
||||
// restore old rounding mode
|
||||
set_rounding_mode( rounding_modes.top() );
|
||||
|
||||
rounding_modes.pop();
|
||||
};
|
||||
|
||||
inline void set_rounding_mode(
|
||||
const Rounding_mode& mode ) const {
|
||||
#ifndef NDEBUG
|
||||
if( scope_id != rounding_modes.size() )
|
||||
CGAL_error_msg( "Rounding mode controller out of scope" );
|
||||
#endif
|
||||
typename FT_kernel::Set_rounding_mode set_mode;
|
||||
set_mode( mode );
|
||||
};
|
||||
|
||||
inline static Rounding_mode get_rounding_mode() {
|
||||
typename FT_kernel::Get_rounding_mode get_mode;
|
||||
return get_mode();
|
||||
};
|
||||
|
||||
private:
|
||||
static std::stack<Rounding_mode> rounding_modes;
|
||||
#ifndef NDEBUG
|
||||
unsigned scope_id;
|
||||
#endif
|
||||
};
|
||||
|
||||
class Precision_controller {
|
||||
public:
|
||||
typedef typename FT_kernel::Precision_type PT;
|
||||
|
||||
Precision_controller( const PT& precision ) {
|
||||
// store current precision
|
||||
precisions.push( get_precision() );
|
||||
|
||||
#ifndef NDEBUG
|
||||
// store scope id
|
||||
scope_id = precisions.size();
|
||||
#endif
|
||||
set_precision( precision );
|
||||
};
|
||||
|
||||
inline ~Precision_controller() {
|
||||
// restore previous precision
|
||||
set_precision( precisions.top() );
|
||||
|
||||
precisions.pop();
|
||||
};
|
||||
|
||||
inline void set_precision( const PT& precision ) const {
|
||||
#ifndef NDEBUG
|
||||
if( scope_id != precisions.size() )
|
||||
CGAL_error_msg( "Precision controller out of scope" );
|
||||
#endif
|
||||
typename FT_kernel::Set_precision set_prec;
|
||||
set_prec( precision );
|
||||
};
|
||||
|
||||
inline static PT get_precision() {
|
||||
typename FT_kernel::Get_precision get_prec;
|
||||
return get_prec();
|
||||
};
|
||||
|
||||
private:
|
||||
static std::stack<PT> precisions;
|
||||
#ifndef NDEBUG
|
||||
unsigned scope_id;
|
||||
#endif
|
||||
};
|
||||
};
|
||||
|
||||
template< class FT_kernel >
|
||||
std::stack< typename FT_kernel::Precision_type >
|
||||
CGAL::Float_traits_base< FT_kernel >::Precision_controller::precisions =
|
||||
std::stack< typename FT_kernel::Precision_type >();
|
||||
|
||||
template< class FT_kernel >
|
||||
std::stack<Rounding_mode>
|
||||
CGAL::Float_traits_base< FT_kernel >::Rounding_mode_controller::rounding_modes =
|
||||
std::stack<Rounding_mode>();
|
||||
|
||||
|
||||
// namespace Intern {
|
||||
namespace CGALi {
|
||||
// IT_float_traits provide functionality for arbitary integer values
|
||||
// such as Division (without and with rounding) and Shifting
|
||||
// which are not explizit definied for integer value, but needed
|
||||
// for the float_type
|
||||
template< class IT_ >
|
||||
class IT_float_traits {
|
||||
public:
|
||||
typedef IT_ IT;
|
||||
|
||||
//! returns the length (in bits) of an integer value $v$.
|
||||
typedef CGAL::Null_functor Length;
|
||||
|
||||
//! calculates for an integer value $e$, $2^e$
|
||||
class Pow_2 {
|
||||
public:
|
||||
IT operator() ( const IT& e ) const {
|
||||
if( e<0 ) return 0;
|
||||
if( e==0 ) return 1;
|
||||
IT a = 1;
|
||||
IT p = 2;
|
||||
IT n = e;
|
||||
while( n > 0 ) {
|
||||
if( n%2 == 1 )
|
||||
a = p*a;
|
||||
n = n/2;
|
||||
p = p*p;
|
||||
}
|
||||
return a;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
class Ilog_2 {
|
||||
public:
|
||||
IT operator()( const IT& a) const {
|
||||
IT log=0;
|
||||
IT pow = 1;
|
||||
while(pow < a)
|
||||
{
|
||||
pow = pow * 2;
|
||||
log = log + 1;
|
||||
}
|
||||
return log;
|
||||
}
|
||||
};
|
||||
|
||||
class Floor_Ilog_2 {
|
||||
public:
|
||||
IT operator()( const IT& a) const {
|
||||
IT log=0;
|
||||
IT pow = 1;
|
||||
IT floor_log;
|
||||
while(pow <= a)
|
||||
{
|
||||
pow = pow * 2;
|
||||
floor_log = log;
|
||||
log = log + 1;
|
||||
}
|
||||
return floor_log;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
//! count the number of zeros at the end in the binary reprasentation of an integer
|
||||
class Zeros {
|
||||
public:
|
||||
IT operator() ( const IT& value ) const {
|
||||
if ( value == IT(0) ) return 0;
|
||||
IT zeros = -1;
|
||||
IT d = value;
|
||||
IT f = 1;
|
||||
do {
|
||||
zeros++;
|
||||
d /= IT(2);
|
||||
f *= IT(2);
|
||||
} while ( d * f == value ) ;
|
||||
|
||||
return zeros;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
//! Division of two integer values $a/b$
|
||||
//! return an integer div which satisfy the equation
|
||||
//! a = div * b + r (r < b )
|
||||
//! Division by zero is not definied.
|
||||
class Integer_Div {
|
||||
public:
|
||||
void help ( const IT&a, const IT& b, IT& sol, IT& sum ) {
|
||||
if ( a<b ){
|
||||
sol = 0;
|
||||
sum = 0;
|
||||
return;
|
||||
}
|
||||
|
||||
if ( a==b ) {
|
||||
sol = 1;
|
||||
sum = b;
|
||||
return ;
|
||||
}
|
||||
|
||||
sol = IT( 1 );
|
||||
sum = b;
|
||||
while ( (sum+sum) <= a ) {
|
||||
sum += sum;
|
||||
sol += sol;
|
||||
}
|
||||
|
||||
return ;
|
||||
}
|
||||
|
||||
IT operator() ( const IT& a, const IT& b ) {
|
||||
typename CGAL::Real_embeddable_traits<IT>::Abs abs;
|
||||
|
||||
CGAL_precondition( b != 0 );
|
||||
|
||||
IT sol = 0; // solution of division
|
||||
IT sum_rem = abs( a );
|
||||
IT abs_b = abs( b );
|
||||
do {
|
||||
IT tmpsol, tmpsum;
|
||||
|
||||
help( sum_rem, abs_b, tmpsol, tmpsum );
|
||||
sol += tmpsol;
|
||||
sum_rem -= tmpsum;
|
||||
} while ( sum_rem >= abs_b );
|
||||
// false => terminate
|
||||
|
||||
|
||||
if ( ( a < 0 && b > 0 )
|
||||
|| ( a > 0 && b < 0 ) ) {
|
||||
sol = (IT)0 - sol;
|
||||
}
|
||||
|
||||
return sol;
|
||||
}
|
||||
}; // class Integer_Div
|
||||
|
||||
class Mod {
|
||||
public:
|
||||
IT operator() ( const IT& a, const IT&b ) {
|
||||
typename IT_float_traits<IT>::Integer_Div int_div;
|
||||
IT divisor = int_div( a, b );
|
||||
return (a-(b*divisor));
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
|
||||
//! shifts an integer value by shift to the left
|
||||
//! (multiplication by $2^shift$)
|
||||
class Shift_left_by {
|
||||
public:
|
||||
IT operator() ( const IT& value, const IT& shift ) const {
|
||||
typename IT_float_traits<IT>::Pow_2 pow_2;
|
||||
return value* pow_2( shift );
|
||||
}
|
||||
}; // class Shift_left_by
|
||||
|
||||
//! shifts an integer value by shift to the right
|
||||
//! (division by $2^shift$)
|
||||
class Shift_right_by {
|
||||
public:
|
||||
IT operator() ( const IT& value, const IT& shift ) const {
|
||||
typename IT_float_traits<IT>::Pow_2 pow_2;
|
||||
typename IT_float_traits<IT>::Integer_Div Integer_div;
|
||||
return Integer_div( value, pow_2( shift ) );
|
||||
}
|
||||
}; // class Shift_right_by
|
||||
|
||||
}; // class IT_float_traits
|
||||
|
||||
|
||||
|
||||
} // namespace CGALi
|
||||
|
||||
|
||||
//! For the exact binary number types of EXACUS,
|
||||
//! this base class provides implementations of
|
||||
//! \c NiX::NT_traits::Floor_log2_abs and \c NiX::NT_traits::Ceil_log2_abs
|
||||
//! based on \c NiX::FloatTraits::Get_mantissa and
|
||||
//! \c NiX::FloatTraits::Get_exponent together with the respective
|
||||
//! \c NiX::NT_traits functor for the mantissa.
|
||||
//! \ingroup NiX_NT_traits_bases
|
||||
template <class NT>
|
||||
class NT_traits_log2_float_traits_base {
|
||||
//TODO
|
||||
//porting Floor_log2_abs and Ceil_log2_abs from EXACUS2CGAL
|
||||
public:
|
||||
//! <tt>NiX::NT_traits::Floor_log2_abs()(NT x)</tt> implemented as
|
||||
//! <tt>NiX::NT_traits<Get_mantissa::result_type>::Floor_log2_abs()(Get_mantissa()(x)) + Get_exponent()(x)</tt>
|
||||
// class Floor_log2_abs {
|
||||
// public:
|
||||
// //! type for the \c AdaptableUnaryFunction concept.
|
||||
// typedef typename Float_traits<NT>::Get_exponent::result_type
|
||||
// result_type;
|
||||
// //! type for the \c AdaptableUnaryFunction concept.
|
||||
// typedef NT argument_type;
|
||||
// //! the function call.
|
||||
// result_type operator() (argument_type x) {
|
||||
// typedef typename Float_traits<NT>::Get_mantissa Get_mantissa;
|
||||
// typedef typename Float_traits<NT>::Get_exponent Get_exponent;
|
||||
// typename NT_traits<typename Get_mantissa::result_type>
|
||||
// ::Floor_log2_abs log;
|
||||
// return log(Get_mantissa()(x)) + Get_exponent()(x);
|
||||
// }
|
||||
// };
|
||||
//! <tt>NiX::NT_traits::Ceil_log2_abs()(NT x)</tt> implemented as
|
||||
//! <tt>NiX::NT_traits<Get_mantissa::result_type>::Ceil_log2_abs()(Get_mantissa()(x)) + Get_exponent()(x)</tt>
|
||||
// class Ceil_log2_abs {
|
||||
// public:
|
||||
// //! type for the \c AdaptableUnaryFunction concept.
|
||||
// typedef typename Float_traits<NT>::Get_exponent::result_type
|
||||
// result_type;
|
||||
// //! type for the \c AdaptableUnaryFunction concept.
|
||||
// typedef NT argument_type;
|
||||
// //! the function call.
|
||||
// result_type operator() (argument_type x) {
|
||||
// typedef typename Float_traits<NT>::Get_mantissa Get_mantissa;
|
||||
// typedef typename Float_traits<NT>::Get_exponent Get_exponent;
|
||||
// typename NT_traits<typename Get_mantissa::result_type>
|
||||
// ::Ceil_log2_abs log;
|
||||
// return log(Get_mantissa()(x)) + Get_exponent()(x);
|
||||
// }
|
||||
// };
|
||||
};
|
||||
|
||||
// } // namespace NiX
|
||||
CGAL_END_NAMESPACE
|
||||
#endif // CGAL_FT_TRAITS_H
|
||||
// EOF
|
||||
|
|
@ -0,0 +1,290 @@
|
|||
// TODO: Add licence
|
||||
//
|
||||
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
|
||||
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
|
||||
//
|
||||
// $URL:$
|
||||
// $Id: 4515 rschindl$
|
||||
//
|
||||
//
|
||||
// Author(s) : Ralf Schindlbeck <rschindl@mpi-inf.mpg.de>
|
||||
//
|
||||
// ============================================================================
|
||||
|
||||
// TODO: The comments are all original EXACUS comments and aren't adapted. So
|
||||
// they may be wrong now.
|
||||
|
||||
|
||||
|
||||
/*! \file CGAL/core_interval_support.h
|
||||
This is experimental
|
||||
use CORE::BigFloat as interval type.
|
||||
*/
|
||||
|
||||
|
||||
#ifndef CGAL_CORE_INTERVAL_SUPPORT_H
|
||||
#define CGAL_CORE_INTERVAL_SUPPORT_H
|
||||
|
||||
#include <CGAL/basic.h>
|
||||
|
||||
#ifndef CGAL_USE_CORE
|
||||
#warning This header file needs CORE installed in order to work properly.
|
||||
#else // CGAL_USE_CORE
|
||||
|
||||
//#include <NiX/CORE/BigFloat.h>
|
||||
//#include <NiX/CORE/Expr.h>
|
||||
//#include <NiX/CORE/BigInt.h>
|
||||
//#include <NiX/CORE/BigRat.h>
|
||||
|
||||
#include <CGAL/CORE/BigFloat.h>
|
||||
|
||||
// namespace NiX{
|
||||
CGAL_BEGIN_NAMESPACE
|
||||
|
||||
template<typename BFI> long get_significant_bits(BFI bfi);
|
||||
|
||||
CORE::BigFloat
|
||||
inline
|
||||
round(const CORE::BigFloat& x, long rel_prec = CORE::defRelPrec.toLong() ){
|
||||
CGAL_postcondition(rel_prec >= 0);
|
||||
// since there is not rel prec defined if in_zero(x)
|
||||
if (x.isZeroIn()) return x;
|
||||
if (CGAL::get_significant_bits(x) <= rel_prec) return x;
|
||||
|
||||
typedef CORE::BigFloat BF;
|
||||
typedef CORE::BigFloat BFI;
|
||||
typedef CORE::BigInt Integer;
|
||||
BF xr;
|
||||
|
||||
CORE::BigInt m = x.m();
|
||||
long err = x.err();
|
||||
long exp = x.exp();
|
||||
|
||||
long shift = ::CORE::bitLength(m) - rel_prec - 1;
|
||||
if( shift > 0 ){
|
||||
Integer new_m = m >> shift ;
|
||||
if(err == 0){
|
||||
xr = BF(new_m,1,0)*BF::exp2(exp*14+shift);
|
||||
}else{
|
||||
xr = BF(new_m,2,0)*BF::exp2(exp*14+shift);
|
||||
}
|
||||
}else{
|
||||
// noting to do
|
||||
xr = x;
|
||||
}
|
||||
|
||||
CGAL_postcondition(CGAL::abs(CGAL::get_significant_bits(xr) - rel_prec) <= 1);
|
||||
CGAL_postcondition(BF(xr.m()-xr.err(),0,xr.exp()) <= BF(x.m()-x.err(),0,x.exp()));
|
||||
CGAL_postcondition(BF(xr.m()+xr.err(),0,xr.exp()) >= BF(x.m()+x.err(),0,x.exp()));
|
||||
return xr;
|
||||
}
|
||||
|
||||
template<typename BFI> class Bigfloat_interval_traits;
|
||||
|
||||
template<> class Bigfloat_interval_traits<CORE::BigFloat>
|
||||
|
||||
{
|
||||
public:
|
||||
typedef CORE::BigFloat NT;
|
||||
typedef CORE::BigFloat BF;
|
||||
|
||||
typedef Bigfloat_interval_traits<NT> Self;
|
||||
|
||||
class Get_significant_bits {
|
||||
public:
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef NT argument_type;
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef long result_type;
|
||||
|
||||
long operator()( NT x) const {
|
||||
//std::cout << ":::::::::::::::::::: "<< std::endl;
|
||||
if(x.err() == 0 ) {
|
||||
//std::cout << "Exact: \n" << x.m() << std::endl;
|
||||
return ::CORE::bitLength(x.m());
|
||||
}
|
||||
else {
|
||||
//std::cout << "With error:\n" << x.m() << std::endl << x.err() << std::endl;
|
||||
//std::cout << "bitlength m: " << ::CORE::bitLength(x.m()) << "\nbitlength e: " << ::CORE::bitLength(x.err()) << std::endl;
|
||||
return ::CORE::bitLength(x.m()) - ::CORE::bitLength(x.err());
|
||||
}
|
||||
|
||||
}
|
||||
};
|
||||
|
||||
class Set_precision {
|
||||
public:
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef long argument_type;
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef long result_type;
|
||||
|
||||
long operator() ( long prec ) const {
|
||||
long result = ::CORE::defRelPrec.toLong();
|
||||
::CORE::defRelPrec = prec;
|
||||
::CORE::defBFdivRelPrec = prec;
|
||||
return result;
|
||||
}
|
||||
};
|
||||
|
||||
class Get_precision {
|
||||
public:
|
||||
// type for the \c AdaptableGenerator concept.
|
||||
typedef long result_type;
|
||||
|
||||
long operator() () const {
|
||||
return ::CORE::defRelPrec.toLong();
|
||||
}
|
||||
};
|
||||
|
||||
class Upper {
|
||||
public:
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef NT argument_type;
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef BF result_type;
|
||||
|
||||
BF operator() ( NT x ) const {
|
||||
CORE::BigFloat result = ::CORE::BigFloat(x.m()+x.err(),0,x.exp());
|
||||
CGAL_postcondition(result >= x);
|
||||
return result;
|
||||
}
|
||||
};
|
||||
|
||||
class Lower {
|
||||
public:
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef NT argument_type;
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef BF result_type;
|
||||
|
||||
BF operator() ( NT x ) const {
|
||||
CORE::BigFloat result = ::CORE::BigFloat(x.m()-x.err(),0,x.exp());
|
||||
CGAL_postcondition(result <= x);
|
||||
return result;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
class In_zero {
|
||||
public:
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef NT argument_type;
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef bool result_type;
|
||||
|
||||
bool operator() ( NT x ) const {
|
||||
return x.isZeroIn();
|
||||
}
|
||||
};
|
||||
|
||||
class Overlap {
|
||||
public:
|
||||
// type for the \c AdaptableBinaryFunction concept.
|
||||
typedef NT first_argument_type;
|
||||
// type for the \c AdaptableBinaryFunction concept.
|
||||
typedef NT second_argument_type;
|
||||
// type for the \c AdaptableBinaryFunction concept.
|
||||
typedef bool result_type;
|
||||
|
||||
bool operator() ( NT x, NT y ) const {
|
||||
Self::In_zero in_zero;
|
||||
bool result = in_zero(x-y);
|
||||
return result;
|
||||
}
|
||||
};
|
||||
|
||||
class Hull {
|
||||
public:
|
||||
// type for the \c AdaptableBinaryFunction concept.
|
||||
typedef NT first_argument_type;
|
||||
// type for the \c AdaptableBinaryFunction concept.
|
||||
typedef NT second_argument_type;
|
||||
// type for the \c AdaptableBinaryFunction concept.
|
||||
typedef NT result_type;
|
||||
|
||||
NT operator() ( NT x, NT y ) const {
|
||||
CORE::BigFloat result;
|
||||
|
||||
// Unfortunately, CORE::centerize(x,y) has bugs.
|
||||
if ((x.m() == y.m()) && (x.err() == y.err()) &&
|
||||
(x.exp() == y.exp())) {
|
||||
|
||||
//rep(x) == rep(y)
|
||||
return x;
|
||||
}
|
||||
Self::Lower lower_functor;
|
||||
Self::Upper upper_functor;
|
||||
CORE::BigFloat lower = std::min(lower_functor(x),
|
||||
lower_functor(y));
|
||||
CORE::BigFloat upper = std::max(upper_functor(x),
|
||||
upper_functor(y));
|
||||
CORE::BigFloat mid = (lower + upper)/2;
|
||||
CORE::BigFloat err = CGAL::round((upper - lower)/2,0);
|
||||
err = CORE::BigFloat(0,err.m().longValue()+err.err(),err.exp());
|
||||
result = mid + err;
|
||||
|
||||
CGAL_postcondition(lower_functor(result) <=
|
||||
std::min(lower_functor(x),
|
||||
lower_functor(y)));
|
||||
CGAL_postcondition(upper_functor(result) >=
|
||||
std::max(upper_functor(x),
|
||||
upper_functor(y)));
|
||||
return result ;
|
||||
}
|
||||
};
|
||||
|
||||
class Singleton {
|
||||
public:
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef NT argument_type;
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef bool result_type;
|
||||
|
||||
bool operator() ( NT x ) const {
|
||||
return (x.err() == 0);
|
||||
}
|
||||
};
|
||||
|
||||
class Width {
|
||||
public:
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef NT argument_type;
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef BF result_type;
|
||||
|
||||
BF operator() ( NT x ) const {
|
||||
unsigned long err = 2*x.err();
|
||||
return BF(CORE::BigInt(err),0,x.exp());
|
||||
}
|
||||
};
|
||||
|
||||
class Convert_to_bfi {
|
||||
public:
|
||||
|
||||
typedef NT result_type;
|
||||
|
||||
|
||||
NT operator() (const ::CORE::Expr& x){
|
||||
return round(CORE::BigFloat(x));
|
||||
}
|
||||
|
||||
NT operator() (const ::CORE::BigInt& x){
|
||||
return round(CORE::BigFloat(x));
|
||||
}
|
||||
|
||||
NT operator() (const ::CORE::BigRat& x){
|
||||
return round(CORE::BigFloat(x));
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
};
|
||||
|
||||
|
||||
|
||||
// } // namespace NiX
|
||||
CGAL_END_NAMESPACE
|
||||
|
||||
#endif // CGAL_USE_CORE
|
||||
#endif // CGAL_CORE_INTERVAL_SUPPORT_H
|
||||
|
|
@ -0,0 +1,380 @@
|
|||
// TODO: Add licence
|
||||
//
|
||||
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
|
||||
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
|
||||
//
|
||||
// $URL:$
|
||||
// $Id: 4515 rschindl$
|
||||
//
|
||||
//
|
||||
// Author(s) : Ralf Schindlbeck <rschindl@mpi-inf.mpg.de>
|
||||
//
|
||||
// ============================================================================
|
||||
|
||||
// TODO: The comments are all original EXACUS comments and aren't adapted. So
|
||||
// they may be wrong now.
|
||||
|
||||
|
||||
/*! \file CGAL/interval_support.h
|
||||
This is experimental
|
||||
*/
|
||||
|
||||
#ifndef CGAL_INTERVAL_SUPPORT_H
|
||||
#define CGAL_INTERVAL_SUPPORT_H
|
||||
|
||||
#include <CGAL/basic.h>
|
||||
#include <CGAL/Arithmetic_kernel.h>
|
||||
|
||||
// #include <NiX/Float_traits.h>
|
||||
#include <CGAL/Float_traits.h>
|
||||
|
||||
///////// covnersion tools:
|
||||
// namespace NiX{
|
||||
CGAL_BEGIN_NAMESPACE
|
||||
|
||||
template <class NT> class Get_arithmetic_traits;
|
||||
|
||||
template <typename A,typename B,typename C,typename D,typename E>
|
||||
class Algebraic_real;
|
||||
|
||||
// } // namespace NiX
|
||||
// namespace CGAL {
|
||||
template <typename A,typename B> class Sqrt_extension;
|
||||
// } // namespace CGAL
|
||||
// namespace NiX {
|
||||
// Bigfloat_interval_functions
|
||||
|
||||
template<typename BigfloatInterval> class Bigfloat_interval_traits;
|
||||
|
||||
template<typename BFI> long get_significant_bits(BFI bfi) {
|
||||
typename Bigfloat_interval_traits<BFI>::Get_significant_bits
|
||||
get_significant_bits;
|
||||
return get_significant_bits(bfi);
|
||||
}
|
||||
|
||||
/*
|
||||
template<typename BFI> long set_precision(BFI bfi,long prec) {
|
||||
typename Float_traits<typename Bigfloat_interval_traits<BFI>::BF>
|
||||
::Set_precision set_precision;
|
||||
return set_precision(prec);
|
||||
}
|
||||
|
||||
template<typename BFI> long get_precision(BFI bfi) {
|
||||
typename Float_traits<typename Bigfloat_interval_traits<BFI>::BF>
|
||||
::Get_precision get_precision;
|
||||
return get_precision();
|
||||
}
|
||||
*/
|
||||
|
||||
// ONLY FOR TESTING THIS IS THE WRONG FILE!!!
|
||||
template<typename BF> long set_precision(BF bfi,long prec) {
|
||||
typename Float_traits<BF>::Set_precision set_precision;
|
||||
return set_precision(prec);
|
||||
}
|
||||
|
||||
template<typename BF> long get_precision(BF bfi) {
|
||||
typename Float_traits<BF>::Get_precision get_precision;
|
||||
return get_precision();
|
||||
}
|
||||
|
||||
|
||||
template<typename BFI>
|
||||
typename Bigfloat_interval_traits<BFI>::BF upper(BFI bfi) {
|
||||
typename Bigfloat_interval_traits<BFI>::Upper upper;
|
||||
return upper(bfi);
|
||||
}
|
||||
|
||||
template<typename BFI>
|
||||
typename Bigfloat_interval_traits<BFI>::BF lower(BFI bfi) {
|
||||
typename Bigfloat_interval_traits<BFI>::Lower lower;
|
||||
return lower(bfi);
|
||||
}
|
||||
|
||||
template<typename BFI> BFI hull(BFI bfi1, BFI bfi2) {
|
||||
typename Bigfloat_interval_traits<BFI>::Hull hull;
|
||||
return hull(bfi1, bfi2);
|
||||
}
|
||||
|
||||
template<typename BFI> bool in_zero(BFI bfi) {
|
||||
typename Bigfloat_interval_traits<BFI>::In_zero in_zero;
|
||||
return in_zero(bfi);
|
||||
}
|
||||
|
||||
template<typename BFI> bool overlap(BFI bfi1, BFI bfi2) {
|
||||
typename Bigfloat_interval_traits<BFI>::Overlap overlap;
|
||||
return overlap(bfi1, bfi2);
|
||||
}
|
||||
|
||||
template<typename BFI> typename Bigfloat_interval_traits<BFI>::BF
|
||||
width(BFI bfi) {
|
||||
|
||||
typename Bigfloat_interval_traits<BFI>::Width width;
|
||||
return width(bfi);
|
||||
}
|
||||
|
||||
template<typename BFI> bool singleton(BFI bfi) {
|
||||
typename Bigfloat_interval_traits<BFI>::Singleton singleton;
|
||||
return singleton(bfi);
|
||||
}
|
||||
|
||||
template<typename NT> typename CGALi::Get_arithmetic_kernel<NT>::Bigfloat_interval
|
||||
convert_to_bfi(NT x) {
|
||||
|
||||
typedef typename
|
||||
CGALi::Get_arithmetic_kernel<NT>::Bigfloat_interval BFI;
|
||||
typename Bigfloat_interval_traits<BFI>::CGALi::Get_arithmetic_kernel
|
||||
get_arithmetic_traits;
|
||||
return get_arithmetic_traits(x);
|
||||
}
|
||||
|
||||
#if 0
|
||||
//TODO
|
||||
//porting Specialisation for double-intervals from EXACUS2CGAL
|
||||
|
||||
// (Parital) Specialisation for double-intervals
|
||||
|
||||
#include <NiX/Interval.h>
|
||||
|
||||
template<>
|
||||
class Bigfloat_interval_traits<Interval>
|
||||
{
|
||||
public:
|
||||
typedef Interval NT;
|
||||
|
||||
typedef double BF;
|
||||
|
||||
/* ??
|
||||
class Get_significant_bits {
|
||||
public:
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef NT argument_type;
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef long result_type;
|
||||
|
||||
long operator()( NT x) const {}
|
||||
};
|
||||
*/
|
||||
|
||||
class Upper {
|
||||
public:
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef NT argument_type;
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef BF result_type;
|
||||
|
||||
BF operator() ( NT a ) const {
|
||||
return a.upper();
|
||||
}
|
||||
};
|
||||
|
||||
class Lower {
|
||||
public:
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef NT argument_type;
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef BF result_type;
|
||||
|
||||
BF operator() ( NT a ) const {
|
||||
return a.lower();
|
||||
}
|
||||
};
|
||||
|
||||
class In_zero {
|
||||
public:
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef NT argument_type;
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef bool result_type;
|
||||
|
||||
bool operator() ( NT x ) const {
|
||||
return ::boost::numeric::in_zero(x);
|
||||
}
|
||||
};
|
||||
|
||||
class Overlap {
|
||||
public:
|
||||
// type for the \c AdaptableBinaryFunction concept.
|
||||
typedef NT first_argument_type;
|
||||
// type for the \c AdaptableBinaryFunction concept.
|
||||
typedef NT second_argument_type;
|
||||
// type for the \c AdaptableBinaryFunction concept.
|
||||
typedef bool result_type;
|
||||
|
||||
bool operator() ( NT x, NT y ) const {
|
||||
return ::boost::numeric::overlap(x,y);
|
||||
}
|
||||
};
|
||||
|
||||
class Hull {
|
||||
public:
|
||||
// type for the \c AdaptableBinaryFunction concept.
|
||||
typedef NT first_argument_type;
|
||||
// type for the \c AdaptableBinaryFunction concept.
|
||||
typedef NT second_argument_type;
|
||||
// type for the \c AdaptableBinaryFunction concept.
|
||||
typedef NT result_type;
|
||||
|
||||
NT operator() ( NT x, NT y ) const {
|
||||
return ::boost::numeric::hull(x,y);
|
||||
}
|
||||
};
|
||||
|
||||
class Singleton {
|
||||
public:
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef NT argument_type;
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef bool result_type;
|
||||
|
||||
bool operator() ( NT a ) const {
|
||||
return ::boost::numeric::singleton(a);
|
||||
}
|
||||
};
|
||||
|
||||
class Width {
|
||||
public:
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef NT argument_type;
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef BF result_type;
|
||||
|
||||
BF operator() ( NT a ) const {
|
||||
return ::boost::numeric::width(a);
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
class Convert_to_bfi {
|
||||
public:
|
||||
|
||||
typedef NT result_type;
|
||||
|
||||
template<typename NTX>
|
||||
NT operator() ( NTX x) {
|
||||
return to_interval( x );
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
};
|
||||
#endif
|
||||
|
||||
|
||||
#if 0
|
||||
template <class COEFF, class HP >
|
||||
leda_bigfloat_interval
|
||||
inline
|
||||
convert_to_bfi(const Algebraic_real< COEFF, ::leda::real, ::leda::rational, HP >& x) {
|
||||
|
||||
//std::cout << "convert_to_bfi(const Algebraic_real<COEFF, FWS, RAT, HP >& x)"<<std::endl;
|
||||
typedef ::leda::real FWS;
|
||||
typedef ::leda::rational RAT;
|
||||
typedef Algebraic_real< COEFF, FWS, RAT, HP > ALG;
|
||||
|
||||
// if zero done
|
||||
if(NiX::sign(x) == CGAL::ZERO)
|
||||
return (leda_bigfloat_interval(0));
|
||||
CGAL_postcondition(NiX::sign(x.low()) == NiX::sign(x.high()));
|
||||
|
||||
|
||||
typedef ::leda::bigfloat BF;
|
||||
long final_prec = BF::set_precision(BF::get_precision()+4);
|
||||
|
||||
::leda::bigfloat rerror(2,-final_prec);
|
||||
|
||||
leda_bigfloat_interval bfi(convert_to_bfi(x.low()).lower(), convert_to_bfi(x.high()).upper());
|
||||
|
||||
//while( NiX::width(bfi) >= rerror){
|
||||
//while( NiX::width(bfi) >= NiX::median(NiX::abs(bfi))*rerror){
|
||||
//std::cout << "rel error: "<< (NiX::width(bfi) / NiX::abs(bfi)).upper() <<std::endl;
|
||||
while( (NiX::width(bfi) / NiX::abs(bfi)).upper() > rerror ){
|
||||
|
||||
//std::cout <<" diff " << NiX::width(bfi) - NiX::abs(bfi).upper()*rerror << std::endl;
|
||||
x.refine();
|
||||
bfi = leda_bigfloat_interval(convert_to_bfi(x.low()).lower(), convert_to_bfi(x.high()).upper());
|
||||
//std::cout <<"bfi: "<< bfi << std::endl;
|
||||
//std::cout << "final_prec: "<< final_prec << std::endl;
|
||||
//std::cout << "rerror : "<< rerror << std::endl;
|
||||
//std::cout << "rel error: "<< (NiX::width(bfi) / NiX::abs(bfi)).upper() <<std::endl;
|
||||
}
|
||||
//std::cout << "convert_to_bfi(const Algebraic_real<COEFF, FWS, RAT, HP >& x) end"<<std::endl;
|
||||
BF::set_precision(final_prec);
|
||||
|
||||
#ifndef NDEBUG
|
||||
::leda::rational lower_num(bfi.lower().get_significant());
|
||||
::leda::rational lower_denom(NiX::ipower(RAT(2),NiX::abs(bfi.lower().get_exponent().to_long())));
|
||||
::leda::rational lower;
|
||||
if(bfi.lower().get_exponent().to_long() < 0 )
|
||||
lower = lower_num / lower_denom ;
|
||||
else
|
||||
lower = lower_num * lower_denom;
|
||||
CGAL_postcondition( x.compare(lower) == CGAL::LARGER );
|
||||
|
||||
::leda::rational upper_num(bfi.upper().get_significant());
|
||||
::leda::rational upper_denom(NiX::ipower(RAT(2),NiX::abs(bfi.upper().get_exponent().to_long())));
|
||||
::leda::rational upper;
|
||||
if(bfi.upper().get_exponent().to_long() < 0 )
|
||||
upper = upper_num / upper_denom ;
|
||||
else
|
||||
upper = upper_num * upper_denom;
|
||||
CGAL_postcondition( x.compare(upper) == CGAL::SMALLER );
|
||||
#endif
|
||||
|
||||
return bfi;
|
||||
}
|
||||
#endif
|
||||
|
||||
template <class NTX>
|
||||
typename CGALi::Get_arithmetic_kernel<NTX>::Arithmetic_kernel::Bigfloat_interval
|
||||
convert_to_bfi(const NTX& x) {
|
||||
typedef typename CGALi::Get_arithmetic_kernel<NTX>::Arithmetic_kernel AT;
|
||||
typedef typename AT::Bigfloat_interval BFI;
|
||||
typename Bigfloat_interval_traits<BFI>::Convert_to_bfi convert_to_bfi;
|
||||
return convert_to_bfi(x);
|
||||
}
|
||||
|
||||
template <typename NT, typename ROOT>
|
||||
typename CGALi::Get_arithmetic_kernel<NT>::Arithmetic_kernel::Bigfloat_interval
|
||||
convert_to_bfi(const CGAL::Sqrt_extension<NT,ROOT>& x) {
|
||||
typedef typename CGALi::Get_arithmetic_kernel<NT>::Arithmetic_kernel AT;
|
||||
typedef typename AT::Bigfloat_interval BFI;
|
||||
if(x.is_extended()){
|
||||
BFI a0(convert_to_bfi(x.a0()));
|
||||
BFI a1(convert_to_bfi(x.a1()));
|
||||
BFI root(CGAL::sqrt(convert_to_bfi(x.root())));
|
||||
return a0+a1*root;
|
||||
}else{
|
||||
return convert_to_bfi(x.a0());
|
||||
}
|
||||
}
|
||||
|
||||
template <class COEFF, class FWS, class RAT, class T1, class T2>
|
||||
typename CGALi::Get_arithmetic_kernel<COEFF>::Arithmetic_kernel::Bigfloat_interval
|
||||
inline
|
||||
convert_to_bfi(const Algebraic_real< COEFF, FWS, RAT,T1,T2>& x) {
|
||||
typedef typename CGALi::Get_arithmetic_kernel<COEFF>::Arithmetic_kernel AT;
|
||||
typedef typename AT::Bigfloat_interval BFI;
|
||||
typedef typename AT::Bigfloat BF;
|
||||
BFI result = x.convert_to_bfi();
|
||||
|
||||
#ifndef NDEBUG
|
||||
CGAL::set_precision(BF(),CGAL::get_precision(BF())*2);
|
||||
CGAL_postcondition(CGAL::lower(result) <= CGAL::lower(CGAL::convert_to_bfi(x.low() )));
|
||||
CGAL_postcondition(CGAL::upper(result) >= CGAL::upper(CGAL::convert_to_bfi(x.high())));
|
||||
CGAL::set_precision(BF(),CGAL::get_precision(BF())/2);
|
||||
#endif
|
||||
return result;
|
||||
};
|
||||
|
||||
// } // namespace NiX
|
||||
CGAL_END_NAMESPACE
|
||||
|
||||
#ifdef CGAL_USE_LEDA
|
||||
#include <CGAL/leda_interval_support.h>
|
||||
#endif // CGAL_USE_LEDA
|
||||
|
||||
#ifdef CGAL_USE_CORE
|
||||
#include <CGAL/core_interval_support.h>
|
||||
#endif // CGAL_USE_CORE
|
||||
|
||||
#endif // CGAL_INTERVAL_SUPPORT_H
|
||||
|
|
@ -0,0 +1,585 @@
|
|||
// TODO: Add licence
|
||||
//
|
||||
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
|
||||
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
|
||||
//
|
||||
// $URL:$
|
||||
// $Id: 4515 rschindl$
|
||||
//
|
||||
//
|
||||
// Author(s) : Ralf Schindlbeck <rschindl@mpi-inf.mpg.de>
|
||||
//
|
||||
// ============================================================================
|
||||
|
||||
// TODO: The comments are all original EXACUS comments and aren't adapted. So
|
||||
// they may be wrong now.
|
||||
|
||||
|
||||
/*! \file CGAL/leda_interval_support.h
|
||||
*/
|
||||
|
||||
#ifndef CGAL_LEDA_INTERVAL_SUPPORT_H
|
||||
#define CGAL_LEDA_INTERVAL_SUPPORT_H
|
||||
|
||||
#include <CGAL/basic.h>
|
||||
#include <CGAL/leda_bigfloat.h>
|
||||
#include <boost/numeric/interval.hpp>
|
||||
|
||||
#include <CGAL/leda_integer.h>
|
||||
#include <CGAL/leda_rational.h>
|
||||
#include <CGAL/leda_real.h>
|
||||
|
||||
// namespace NiX {
|
||||
CGAL_BEGIN_NAMESPACE
|
||||
|
||||
// namespace Intern {
|
||||
namespace CGALi {
|
||||
|
||||
|
||||
struct Rounding_for_leda_bigfloat {
|
||||
private: typedef leda::bigfloat T;
|
||||
public:
|
||||
Rounding_for_leda_bigfloat(){};
|
||||
~Rounding_for_leda_bigfloat(){};
|
||||
|
||||
T conv_down(const T& a){
|
||||
return round(a,leda::bigfloat::get_precision(),leda::TO_N_INF);
|
||||
};
|
||||
T conv_up (const T& a){
|
||||
return round(a,leda::bigfloat::get_precision(),leda::TO_P_INF);
|
||||
};
|
||||
// mathematical operations
|
||||
T add_down(const T& a, const T& b){
|
||||
return add(a,b,leda::bigfloat::get_precision(),leda::TO_N_INF);
|
||||
};
|
||||
T add_up (const T& a, const T& b){
|
||||
return add(a,b,leda::bigfloat::get_precision(),leda::TO_P_INF);
|
||||
};
|
||||
T sub_down(const T& a, const T& b){
|
||||
return sub(a, b, leda::bigfloat::get_precision(),leda::TO_N_INF);
|
||||
};
|
||||
T sub_up (const T& a, const T& b){
|
||||
return sub(a, b, leda::bigfloat::get_precision(),leda::TO_P_INF);
|
||||
};
|
||||
T mul_down(const T& a, const T& b){
|
||||
return mul(a, b, leda::bigfloat::get_precision(),leda::TO_N_INF);
|
||||
};
|
||||
T mul_up (const T& a, const T& b){
|
||||
return mul(a, b, leda::bigfloat::get_precision(),leda::TO_P_INF);
|
||||
};
|
||||
T div_down(const T& a, const T& b){
|
||||
return div(a, b, leda::bigfloat::get_precision(),leda::TO_N_INF);
|
||||
};
|
||||
T div_up (const T& a, const T& b){
|
||||
return div(a, b, leda::bigfloat::get_precision(),leda::TO_P_INF);
|
||||
};
|
||||
T sqrt_down(const T& a){
|
||||
return sqrt(a, leda::bigfloat::get_precision(),leda::TO_N_INF);
|
||||
};
|
||||
T sqrt_up (const T& a){
|
||||
return sqrt(a, leda::bigfloat::get_precision(),leda::TO_P_INF);
|
||||
};
|
||||
|
||||
T median(const T& a, const T& b){ return (a+b)/2; };
|
||||
T int_down(const T& a) { return T(floor(a));};
|
||||
T int_up (const T& a) { return T(ceil(a)); };
|
||||
|
||||
/*
|
||||
T exp_down(T);
|
||||
T exp_up (T);
|
||||
T log_down(T);
|
||||
T log_up (T);
|
||||
T cos_down(T);
|
||||
T cos_up (T);
|
||||
T tan_down(T);
|
||||
T tan_up (T);
|
||||
T asin_down(T); // [-1;1]
|
||||
T asin_up (T); // [-1;1]
|
||||
T acos_down(T); // [-1;1]
|
||||
T acos_up (T); // [-1;1]
|
||||
T atan_down(T); // [-?;+?]
|
||||
T atan_up (T); // [-?;+?]
|
||||
T sinh_down(T); // [-?;+?]
|
||||
T sinh_up (T); // [-?;+?]
|
||||
T cosh_down(T); // [-?;+?]
|
||||
T cosh_up (T); // [-?;+?]
|
||||
T tanh_down(T); // [-?;+?]
|
||||
T tanh_up (T); // [-?;+?]
|
||||
T asinh_down(T); // [-?;+?]
|
||||
T asinh_up (T); // [-?;+?]
|
||||
T acosh_down(T); // [1;+?]
|
||||
T acosh_up (T); // [1;+?]
|
||||
T atanh_down(T); // [-1;1]
|
||||
T atanh_up (T); // [-1;1]
|
||||
*/
|
||||
|
||||
// unprotected rounding class
|
||||
//typedef ... unprotected_rounding;
|
||||
};
|
||||
|
||||
class Checking_for_leda_bigfloat {
|
||||
|
||||
typedef leda::bigfloat T;
|
||||
|
||||
public:
|
||||
|
||||
static T pos_inf() {
|
||||
T b = T(5) / T(0);
|
||||
CGAL_assertion(leda::ispInf(b));
|
||||
return b;
|
||||
}
|
||||
|
||||
static T neg_inf() {
|
||||
T b = T(-5) / T(0);
|
||||
CGAL_assertion(leda::isnInf(b));
|
||||
return b;
|
||||
}
|
||||
|
||||
static T nan() {
|
||||
T b = T(0)*pos_inf();
|
||||
CGAL_assertion(leda::isNaN(b));
|
||||
return b;
|
||||
}
|
||||
|
||||
static bool is_nan(const T& b) {
|
||||
return leda::isNaN(b);
|
||||
}
|
||||
|
||||
static T empty_lower() {
|
||||
return T(1);
|
||||
}
|
||||
|
||||
static T empty_upper() {
|
||||
return T(0);
|
||||
}
|
||||
|
||||
static bool is_empty(const T& a, const T& b) {
|
||||
return a==T(1) && b == T(0);
|
||||
}
|
||||
};
|
||||
|
||||
} // namespace CGALi
|
||||
// } // namespace NiX
|
||||
CGAL_END_NAMESPACE
|
||||
|
||||
namespace boost {
|
||||
namespace numeric {
|
||||
|
||||
|
||||
inline
|
||||
std::ostream& operator <<
|
||||
(std::ostream& os, const boost::numeric::interval<leda::bigfloat>& x)
|
||||
{
|
||||
os << "["
|
||||
<< x.lower().get_significant() << "*2^" << x.lower().get_exponent()
|
||||
<< " , "
|
||||
<< x.upper().get_significant() << "*2^" << x.upper().get_exponent()
|
||||
<< "]";
|
||||
return os;
|
||||
}
|
||||
|
||||
|
||||
}//namespace numeric
|
||||
}//namespace boost
|
||||
|
||||
// namespace NiX {
|
||||
CGAL_BEGIN_NAMESPACE
|
||||
|
||||
typedef boost::numeric::interval
|
||||
< leda::bigfloat,
|
||||
boost::numeric::interval_lib::policies
|
||||
< CGALi::Rounding_for_leda_bigfloat,
|
||||
CGALi::Checking_for_leda_bigfloat > >
|
||||
leda_bigfloat_interval;
|
||||
|
||||
|
||||
|
||||
|
||||
/*
|
||||
CGAL::Sign inline sign(const leda::bigfloat& x){
|
||||
if (x < 0 ) return CGAL::NEGATIVE;
|
||||
if (x > 0 ) return CGAL::POSITIVE;
|
||||
return CGAL::ZERO;
|
||||
}
|
||||
*/
|
||||
|
||||
leda_bigfloat_interval inline ipower(const leda_bigfloat_interval& x, int i ){
|
||||
return ::boost::numeric::pow(x,i);
|
||||
}
|
||||
|
||||
::leda::bigfloat inline median(const leda_bigfloat_interval& x){
|
||||
return ::boost::numeric::median(x);
|
||||
}
|
||||
|
||||
// forward
|
||||
template<typename BFI> class Bigfloat_interval_traits;
|
||||
template<>
|
||||
class Bigfloat_interval_traits<leda_bigfloat_interval>
|
||||
|
||||
{
|
||||
public:
|
||||
typedef leda_bigfloat_interval NT;
|
||||
|
||||
typedef leda::bigfloat BF;
|
||||
|
||||
class Get_significant_bits {
|
||||
public:
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef NT argument_type;
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef long result_type;
|
||||
|
||||
long operator()( NT x) const {
|
||||
leda::bigfloat lower = x.lower();
|
||||
leda::bigfloat upper = x.upper();
|
||||
|
||||
leda::integer lower_m = lower.get_significant();
|
||||
leda::integer upper_m = upper.get_significant();
|
||||
|
||||
leda::integer lower_exp = lower.get_exponent();
|
||||
leda::integer upper_exp = upper.get_exponent();
|
||||
|
||||
long shift = (upper_exp - lower_exp).to_long();
|
||||
if(shift >= 0 ) upper_m = (upper_m << shift);
|
||||
else lower_m = (lower_m << -shift);
|
||||
|
||||
//CGAL_postcondition(lower_m.length() == upper_m.length());
|
||||
|
||||
leda::integer err = lower_m-upper_m;
|
||||
|
||||
return std::max(lower_m.length()-err.length(),0);
|
||||
|
||||
}
|
||||
};
|
||||
|
||||
class Upper {
|
||||
public:
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef NT argument_type;
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef BF result_type;
|
||||
|
||||
BF operator() ( NT a ) const {
|
||||
return a.upper();
|
||||
}
|
||||
};
|
||||
|
||||
class Lower {
|
||||
public:
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef NT argument_type;
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef BF result_type;
|
||||
|
||||
BF operator() ( NT a ) const {
|
||||
return a.lower();
|
||||
}
|
||||
};
|
||||
|
||||
class Set_precision {
|
||||
public:
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef long argument_type;
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef long result_type;
|
||||
|
||||
long operator() ( long prec ) const {
|
||||
return BF::set_precision(prec);
|
||||
}
|
||||
};
|
||||
|
||||
class Get_precision {
|
||||
public:
|
||||
// type for the \c AdaptableGenerator concept.
|
||||
typedef long result_type;
|
||||
|
||||
long operator() () const {
|
||||
return BF::get_precision();
|
||||
}
|
||||
};
|
||||
|
||||
class In_zero {
|
||||
public:
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef NT argument_type;
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef bool result_type;
|
||||
|
||||
bool operator() ( NT x ) const {
|
||||
return ::boost::numeric::in_zero(x);
|
||||
}
|
||||
};
|
||||
|
||||
class Overlap {
|
||||
public:
|
||||
// type for the \c AdaptableBinaryFunction concept.
|
||||
typedef NT first_argument_type;
|
||||
// type for the \c AdaptableBinaryFunction concept.
|
||||
typedef NT second_argument_type;
|
||||
// type for the \c AdaptableBinaryFunction concept.
|
||||
typedef bool result_type;
|
||||
|
||||
bool operator() ( NT x, NT y ) const {
|
||||
return ::boost::numeric::overlap(x,y);
|
||||
}
|
||||
};
|
||||
|
||||
class Hull {
|
||||
public:
|
||||
// type for the \c AdaptableBinaryFunction concept.
|
||||
typedef NT first_argument_type;
|
||||
// type for the \c AdaptableBinaryFunction concept.
|
||||
typedef NT second_argument_type;
|
||||
// type for the \c AdaptableBinaryFunction concept.
|
||||
typedef NT result_type;
|
||||
|
||||
NT operator() ( NT x, NT y ) const {
|
||||
return ::boost::numeric::hull(x,y);
|
||||
}
|
||||
};
|
||||
|
||||
class Singleton {
|
||||
public:
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef NT argument_type;
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef bool result_type;
|
||||
|
||||
bool operator() ( NT a ) const {
|
||||
return ::boost::numeric::singleton(a);
|
||||
}
|
||||
};
|
||||
|
||||
class Width {
|
||||
public:
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef NT argument_type;
|
||||
// type for the \c AdaptableUnaryFunction concept.
|
||||
typedef BF result_type;
|
||||
|
||||
BF operator() ( NT a ) const {
|
||||
return ::boost::numeric::width(a);
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
class Convert_to_bfi {
|
||||
public:
|
||||
|
||||
typedef NT result_type;
|
||||
|
||||
|
||||
NT operator() ( const leda::real& x ) {
|
||||
long current_prec = ::leda::bigfloat::get_precision();
|
||||
//x.improve_approximation_to(current_prec);
|
||||
x.guarantee_relative_error(current_prec);
|
||||
|
||||
leda::bigfloat bnum = x.to_bigfloat();
|
||||
leda::bigfloat berr = x.get_bigfloat_error();
|
||||
|
||||
leda::bigfloat low
|
||||
= leda::sub(bnum,berr,current_prec,LEDA::TO_N_INF);
|
||||
leda::bigfloat upp
|
||||
= leda::add(bnum,berr,current_prec,LEDA::TO_P_INF);
|
||||
leda_bigfloat_interval bfi(low,upp) ;
|
||||
|
||||
// std::cout <<"x: "<< x << std::endl;
|
||||
// std::cout <<"bfi.lower(): "<< bfi.lower() << std::endl;
|
||||
// std::cout <<"bfi.upper(): "<< bfi.upper() << std::endl;
|
||||
|
||||
CGAL_postcondition( bfi.lower() <= x );
|
||||
CGAL_postcondition( bfi.upper() >= x );
|
||||
|
||||
return bfi;
|
||||
}
|
||||
|
||||
|
||||
NT operator() (const ::leda::integer& x) {
|
||||
long current_prec = ::leda::bigfloat::get_precision();
|
||||
leda_bigfloat_interval bfi;
|
||||
long length = x.length();
|
||||
|
||||
if(length > current_prec) {
|
||||
::leda::integer significant
|
||||
= CGAL::abs(x) >> (length - current_prec);
|
||||
::leda::bigfloat lower,upper;
|
||||
if(x > 0){
|
||||
lower = ::leda::bigfloat(significant,length - current_prec);
|
||||
upper = ::leda::bigfloat(significant+1,length - current_prec);
|
||||
}else{
|
||||
lower
|
||||
= -::leda::bigfloat(significant+1,length - current_prec);
|
||||
upper
|
||||
= -::leda::bigfloat(significant,length - current_prec);
|
||||
}
|
||||
bfi = leda_bigfloat_interval(lower,upper);
|
||||
}else{
|
||||
::leda::bigfloat bf(x);
|
||||
bfi = leda_bigfloat_interval(bf,bf);
|
||||
}
|
||||
CGAL_postcondition( bfi.lower() <= x );
|
||||
CGAL_postcondition( bfi.upper() >= x );
|
||||
return bfi;
|
||||
}
|
||||
|
||||
|
||||
NT operator() (const ::leda::rational& x) {
|
||||
long old_prec = ::leda::bigfloat::get_precision();
|
||||
::leda::bigfloat::set_precision(old_prec*2);
|
||||
Bigfloat_interval_traits<NT>::Convert_to_bfi convert_to_bfi;
|
||||
leda_bigfloat_interval num = convert_to_bfi(x.numerator());
|
||||
leda_bigfloat_interval den = convert_to_bfi(x.denominator());
|
||||
::leda::bigfloat::set_precision(old_prec);
|
||||
leda_bigfloat_interval bfi = num/den;
|
||||
CGAL_postcondition( bfi.lower() <= x );
|
||||
CGAL_postcondition( bfi.upper() >= x );
|
||||
return bfi;
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
};
|
||||
|
||||
template <> class Algebraic_structure_traits< leda_bigfloat_interval >
|
||||
: public Algebraic_structure_traits_base< leda_bigfloat_interval,
|
||||
Field_with_sqrt_tag > {
|
||||
public:
|
||||
typedef Tag_false Is_exact;
|
||||
typedef Tag_true Is_numerical_sensitive;
|
||||
|
||||
class Sqrt
|
||||
: public Unary_function< Type, Type > {
|
||||
public:
|
||||
Type operator()( const Type& x ) const {
|
||||
return ::boost::numeric::sqrt(x);
|
||||
}
|
||||
};
|
||||
};
|
||||
|
||||
template <> class Real_embeddable_traits< leda_bigfloat_interval >
|
||||
: public Real_embeddable_traits_base< leda_bigfloat_interval > {
|
||||
public:
|
||||
|
||||
class Abs
|
||||
: public Unary_function< Type, Type > {
|
||||
public:
|
||||
Type operator()( const Type& x ) const {
|
||||
return ::boost::numeric::abs(x);
|
||||
}
|
||||
};
|
||||
|
||||
class To_double
|
||||
: public Unary_function< Type, double > {
|
||||
public:
|
||||
double operator()( const Type& x ) const {
|
||||
return CGAL::to_double(::boost::numeric::median(x));
|
||||
}
|
||||
};
|
||||
|
||||
class To_interval
|
||||
: public Unary_function< Type, std::pair< double, double > > {
|
||||
public:
|
||||
std::pair<double, double> operator()( const Type& x ) const {
|
||||
// TODO
|
||||
// return ::boost::numeric::hull(CGAL::to_interval(x.lower()),CGAL::to_interval(x.upper()));
|
||||
|
||||
return std::pair< double, double >( CGAL::to_double( x.lower() ),
|
||||
CGAL::to_double( x.upper() ) );
|
||||
}
|
||||
};
|
||||
};
|
||||
|
||||
|
||||
|
||||
// template <>
|
||||
// class NT_traits<leda_bigfloat_interval>
|
||||
// : public NT_traits_base<leda_bigfloat_interval,
|
||||
// Field_with_sqrt_tag>,
|
||||
// public NT_traits_comparable_base<leda_bigfloat_interval>
|
||||
//
|
||||
// {
|
||||
// public:
|
||||
// typedef leda_bigfloat_interval NT;
|
||||
//
|
||||
// class Abs {
|
||||
// public:
|
||||
// // type for the \c AdaptableUnaryFunction concept.
|
||||
// typedef NT argument_type;
|
||||
// // type for the \c AdaptableUnaryFunction concept.
|
||||
// typedef NT result_type;
|
||||
// NT operator()( NT a) const {
|
||||
// return ::boost::numeric::abs(a);
|
||||
// };
|
||||
// };
|
||||
//
|
||||
// class Sqrt {
|
||||
// public:
|
||||
// // type for the \c AdaptableUnaryFunction concept.
|
||||
// typedef NT argument_type;
|
||||
// // type for the \c AdaptableUnaryFunction concept.
|
||||
// typedef NT result_type;
|
||||
// NT operator()( NT a) const {
|
||||
// return ::boost::numeric::sqrt(a);
|
||||
// };
|
||||
// };
|
||||
//
|
||||
// class To_double {
|
||||
// public:
|
||||
// // the result type.
|
||||
// typedef double result_type;
|
||||
// // the argument type.
|
||||
// typedef NT argument_type;
|
||||
//
|
||||
// double operator()(const NT& a) const {
|
||||
// return NiX::to_double(::boost::numeric::median(a));
|
||||
// }
|
||||
// };
|
||||
// class To_Interval {
|
||||
// public:
|
||||
// // the result type.
|
||||
// typedef Interval result_type;
|
||||
// // the argument type.
|
||||
// typedef NT argument_type;
|
||||
//
|
||||
// Interval operator()(const NT& a) const {
|
||||
// return ::boost::numeric::hull(NiX::to_Interval(a.lower()),
|
||||
// NiX::to_Interval(a.upper()));
|
||||
// }
|
||||
// };
|
||||
|
||||
//TODO
|
||||
//porting Floor_log2_abs and Ceil_log2_abs from EXACUS2CGAL
|
||||
// class Floor_log2_abs {
|
||||
// public:
|
||||
// typedef NT argument_type;
|
||||
// typedef long result_type;
|
||||
// result_type operator() (argument_type x) {
|
||||
// CGAL_precondition(! ::boost::numeric::in_zero(x));
|
||||
// return floor_log2_abs(::boost::numeric::abs(x).lower());
|
||||
// }
|
||||
// };
|
||||
// class Ceil_log2_abs {
|
||||
// public:
|
||||
// typedef NT argument_type;
|
||||
// typedef long result_type;
|
||||
// result_type operator() (argument_type x) {
|
||||
// CGAL_precondition(!(::boost::numeric::in_zero(x) &&
|
||||
// ::boost::numeric::singleton(x)));
|
||||
// return ceil_log2_abs(::boost::numeric::abs(x).upper());
|
||||
// }
|
||||
// };
|
||||
// };
|
||||
|
||||
::leda::bigfloat inline relative_error(const leda_bigfloat_interval& x){
|
||||
if(in_zero(x)){
|
||||
return CGAL::abs(x).upper();
|
||||
}else{
|
||||
return (width(x) / CGAL::abs(x)).upper();
|
||||
}
|
||||
}
|
||||
|
||||
// } // namespace NiX
|
||||
CGAL_END_NAMESPACE
|
||||
|
||||
|
||||
#endif // CGAL_LEDA_INTERVAL_SUPPORT_H
|
||||
Loading…
Reference in New Issue