From edd70f970a2c348b6d39cc7575a6fadb484c19d0 Mon Sep 17 00:00:00 2001 From: Ralf Schindlbeck Date: Thu, 14 Feb 2008 09:33:09 +0000 Subject: [PATCH] interval_support, is now ported from EXACUS2CGAL --- .gitattributes | 4 + Number_types/include/CGAL/Float_traits.h | 477 ++++++++++++++ .../include/CGAL/core_interval_support.h | 290 +++++++++ Number_types/include/CGAL/interval_support.h | 380 ++++++++++++ .../include/CGAL/leda_interval_support.h | 585 ++++++++++++++++++ 5 files changed, 1736 insertions(+) create mode 100644 Number_types/include/CGAL/Float_traits.h create mode 100644 Number_types/include/CGAL/core_interval_support.h create mode 100644 Number_types/include/CGAL/interval_support.h create mode 100644 Number_types/include/CGAL/leda_interval_support.h diff --git a/.gitattributes b/.gitattributes index 35d2d994b1f..afd6814839b 100644 --- a/.gitattributes +++ b/.gitattributes @@ -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 diff --git a/Number_types/include/CGAL/Float_traits.h b/Number_types/include/CGAL/Float_traits.h new file mode 100644 index 00000000000..c72b8ed99fb --- /dev/null +++ b/Number_types/include/CGAL/Float_traits.h @@ -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 +// +// ============================================================================ + +// 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 + +// #include + +#include + +// namespace NiX { +CGAL_BEGIN_NAMESPACE + +//compensation for the #include + +//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_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 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 +CGAL::Float_traits_base< FT_kernel >::Rounding_mode_controller::rounding_modes = + std::stack(); + + +// 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::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::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::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::Pow_2 pow_2; + typename IT_float_traits::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_traits_log2_float_traits_base { +//TODO +//porting Floor_log2_abs and Ceil_log2_abs from EXACUS2CGAL +public: + //! NiX::NT_traits::Floor_log2_abs()(NT x) implemented as + //! NiX::NT_traits::Floor_log2_abs()(Get_mantissa()(x)) + Get_exponent()(x) +// class Floor_log2_abs { +// public: +// //! type for the \c AdaptableUnaryFunction concept. +// typedef typename Float_traits::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::Get_mantissa Get_mantissa; +// typedef typename Float_traits::Get_exponent Get_exponent; +// typename NT_traits +// ::Floor_log2_abs log; +// return log(Get_mantissa()(x)) + Get_exponent()(x); +// } +// }; + //! NiX::NT_traits::Ceil_log2_abs()(NT x) implemented as + //! NiX::NT_traits::Ceil_log2_abs()(Get_mantissa()(x)) + Get_exponent()(x) +// class Ceil_log2_abs { +// public: +// //! type for the \c AdaptableUnaryFunction concept. +// typedef typename Float_traits::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::Get_mantissa Get_mantissa; +// typedef typename Float_traits::Get_exponent Get_exponent; +// typename NT_traits +// ::Ceil_log2_abs log; +// return log(Get_mantissa()(x)) + Get_exponent()(x); +// } +// }; +}; + +// } // namespace NiX +CGAL_END_NAMESPACE +#endif // CGAL_FT_TRAITS_H +// EOF diff --git a/Number_types/include/CGAL/core_interval_support.h b/Number_types/include/CGAL/core_interval_support.h new file mode 100644 index 00000000000..e29182356fe --- /dev/null +++ b/Number_types/include/CGAL/core_interval_support.h @@ -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 +// +// ============================================================================ + +// 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 + +#ifndef CGAL_USE_CORE +#warning This header file needs CORE installed in order to work properly. +#else // CGAL_USE_CORE + +//#include +//#include +//#include +//#include + +#include + +// namespace NiX{ +CGAL_BEGIN_NAMESPACE + +template 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 class Bigfloat_interval_traits; + +template<> class Bigfloat_interval_traits + + { + public: + typedef CORE::BigFloat NT; + typedef CORE::BigFloat BF; + + typedef Bigfloat_interval_traits 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 diff --git a/Number_types/include/CGAL/interval_support.h b/Number_types/include/CGAL/interval_support.h new file mode 100644 index 00000000000..eabd86ea1a3 --- /dev/null +++ b/Number_types/include/CGAL/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 +// +// ============================================================================ + +// 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 +#include + +// #include +#include + +///////// covnersion tools: +// namespace NiX{ +CGAL_BEGIN_NAMESPACE + +template class Get_arithmetic_traits; + +template +class Algebraic_real; + +// } // namespace NiX +// namespace CGAL { +template class Sqrt_extension; +// } // namespace CGAL +// namespace NiX { +// Bigfloat_interval_functions + +template class Bigfloat_interval_traits; + +template long get_significant_bits(BFI bfi) { + typename Bigfloat_interval_traits::Get_significant_bits + get_significant_bits; + return get_significant_bits(bfi); +} + +/* + template long set_precision(BFI bfi,long prec) { + typename Float_traits::BF> + ::Set_precision set_precision; + return set_precision(prec); +} + +template long get_precision(BFI bfi) { + typename Float_traits::BF> + ::Get_precision get_precision; + return get_precision(); +} +*/ + +// ONLY FOR TESTING THIS IS THE WRONG FILE!!! + template long set_precision(BF bfi,long prec) { + typename Float_traits::Set_precision set_precision; + return set_precision(prec); +} + +template long get_precision(BF bfi) { + typename Float_traits::Get_precision get_precision; + return get_precision(); +} + + +template + typename Bigfloat_interval_traits::BF upper(BFI bfi) { + typename Bigfloat_interval_traits::Upper upper; + return upper(bfi); +} + +template + typename Bigfloat_interval_traits::BF lower(BFI bfi) { + typename Bigfloat_interval_traits::Lower lower; + return lower(bfi); +} + +template BFI hull(BFI bfi1, BFI bfi2) { + typename Bigfloat_interval_traits::Hull hull; + return hull(bfi1, bfi2); +} + +template bool in_zero(BFI bfi) { + typename Bigfloat_interval_traits::In_zero in_zero; + return in_zero(bfi); +} + +template bool overlap(BFI bfi1, BFI bfi2) { + typename Bigfloat_interval_traits::Overlap overlap; + return overlap(bfi1, bfi2); +} + +template typename Bigfloat_interval_traits::BF + width(BFI bfi) { + + typename Bigfloat_interval_traits::Width width; + return width(bfi); +} + +template bool singleton(BFI bfi) { + typename Bigfloat_interval_traits::Singleton singleton; + return singleton(bfi); +} + +template typename CGALi::Get_arithmetic_kernel::Bigfloat_interval + convert_to_bfi(NT x) { + + typedef typename + CGALi::Get_arithmetic_kernel::Bigfloat_interval BFI; + typename Bigfloat_interval_traits::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 + +template<> + class Bigfloat_interval_traits + { + 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 + NT operator() ( NTX x) { + return to_interval( x ); + } + + }; + + }; +#endif + + +#if 0 +template +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& x)"< 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() < 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() <& x) end"< +typename CGALi::Get_arithmetic_kernel::Arithmetic_kernel::Bigfloat_interval +convert_to_bfi(const NTX& x) { + typedef typename CGALi::Get_arithmetic_kernel::Arithmetic_kernel AT; + typedef typename AT::Bigfloat_interval BFI; + typename Bigfloat_interval_traits::Convert_to_bfi convert_to_bfi; + return convert_to_bfi(x); +} + +template +typename CGALi::Get_arithmetic_kernel::Arithmetic_kernel::Bigfloat_interval +convert_to_bfi(const CGAL::Sqrt_extension& x) { + typedef typename CGALi::Get_arithmetic_kernel::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 +typename CGALi::Get_arithmetic_kernel::Arithmetic_kernel::Bigfloat_interval +inline +convert_to_bfi(const Algebraic_real< COEFF, FWS, RAT,T1,T2>& x) { + typedef typename CGALi::Get_arithmetic_kernel::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 +#endif // CGAL_USE_LEDA + +#ifdef CGAL_USE_CORE +#include +#endif // CGAL_USE_CORE + +#endif // CGAL_INTERVAL_SUPPORT_H diff --git a/Number_types/include/CGAL/leda_interval_support.h b/Number_types/include/CGAL/leda_interval_support.h new file mode 100644 index 00000000000..e67a6d5e7e2 --- /dev/null +++ b/Number_types/include/CGAL/leda_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 +// +// ============================================================================ + +// 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 +#include +#include + +#include +#include +#include + +// 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& 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 class Bigfloat_interval_traits; + template<> + class Bigfloat_interval_traits + + { + 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::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 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 +// : public NT_traits_base, +// public NT_traits_comparable_base +// +// { +// 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