songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #3614 from mglisse/Number_types-boost_to_interval-glisse 

Rewrite to_interval for boost.multiprecision types
This commit is contained in:
Laurent Rineau 2019-02-07 12:48:46 +00:00 committed by GitHub
parent b3fbbdfd83 811ab8ac43
commit 2f903894e5
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
2 changed files with 143 additions and 46 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

178
Number_types/include/CGAL/boost_mp.h
View File

@ -57,6 +57,9 @@
# include <CGAL/enable_warnings.h> # include <CGAL/enable_warnings.h>
#endif #endif
#ifdef CGAL_USE_MPFR
# include <mpfr.h>
#endif
// TODO: work on the coercions (end of the file) // TODO: work on the coercions (end of the file)
@ -136,6 +139,13 @@ struct AST_boost_mp <NT, boost::mpl::int_<boost::multiprecision::number_kind_rat
}; };
}; };
template <class Backend, boost::multiprecision::expression_template_option Eto>
struct Algebraic_structure_traits<boost::multiprecision::number<Backend, Eto> >
: AST_boost_mp <boost::multiprecision::number<Backend, Eto> > {};
template <class T1,class T2,class T3,class T4,class T5>
struct Algebraic_structure_traits<boost::multiprecision::detail::expression<T1,T2,T3,T4,T5> >
: Algebraic_structure_traits<typename boost::multiprecision::detail::expression<T1,T2,T3,T4,T5>::result_type > {};
// Real_embeddable_traits // Real_embeddable_traits
template <class NT> template <class NT>
@ -193,11 +203,22 @@ struct RET_boost_mp_base
: public CGAL::cpp98::unary_function< Type, std::pair< double, double > > { : public CGAL::cpp98::unary_function< Type, std::pair< double, double > > {
std::pair<double, double> std::pair<double, double>
operator()(const Type& x) const { operator()(const Type& x) const {
// See if https://github.com/boostorg/multiprecision/issues/108 suggests anything better
// assume the conversion is within 1 ulp // assume the conversion is within 1 ulp
// adding IA::smallest() doesn't work because inf-e=inf, even rounded down. // adding IA::smallest() doesn't work because inf-e=inf, even rounded down.
double d = x.template convert_to<double>(); double i = x.template convert_to<double>();
double s = i;
double inf = std::numeric_limits<double>::infinity(); double inf = std::numeric_limits<double>::infinity();
return std::pair<double, double> (nextafter (d, -inf), nextafter (d, inf)); int cmp = x.compare(i);
if (cmp > 0) {
s = nextafter(s, +inf);
CGAL_assertion(x.compare(s) < 0);
}
else if (cmp < 0) {
i = nextafter(i, -inf);
CGAL_assertion(x.compare(i) > 0);
}
return std::pair<double, double> (i, s);
} }
}; };
}; };
@ -211,26 +232,93 @@ struct RET_boost_mp <NT, boost::mpl::int_<boost::multiprecision::number_kind_int
template <class NT> template <class NT>
struct RET_boost_mp <NT, boost::mpl::int_<boost::multiprecision::number_kind_rational> > struct RET_boost_mp <NT, boost::mpl::int_<boost::multiprecision::number_kind_rational> >
: RET_boost_mp_base <NT> { : RET_boost_mp_base <NT> {};
#if BOOST_VERSION < 105700
typedef NT Type; #ifdef CGAL_USE_MPFR
// Because of these full specializations, things get instantiated more eagerly. Make them artificially partial if necessary.
template <>
struct RET_boost_mp <boost::multiprecision::mpz_int>
: RET_boost_mp_base <boost::multiprecision::mpz_int> {
typedef boost::multiprecision::mpz_int Type;
struct To_interval struct To_interval
: public CGAL::cpp98::unary_function< Type, std::pair< double, double > > { : public CGAL::cpp98::unary_function< Type, std::pair< double, double > > {
std::pair<double, double> std::pair<double, double>
operator()(const Type& x) const { operator()(const Type& x) const {
std::pair<double,double> p_num = CGAL::to_interval (numerator (x)); #if MPFR_VERSION_MAJOR >= 3
std::pair<double,double> p_den = CGAL::to_interval (denominator (x)); MPFR_DECL_INIT (y, 53); /* Assume IEEE-754 */
typedef Interval_nt<false> IA; int r = mpfr_set_z (y, x.backend().data(), MPFR_RNDA);
IA::Protector P; double i = mpfr_get_d (y, MPFR_RNDA); /* EXACT but can overflow */
// assume the conversion is within 1 ulp for integers, but the conversion from rational may be unsafe, see boost trac #10085 if (r == 0 && is_finite (i))
IA i_num (p_num.first, p_num.second); return std::pair<double, double>(i, i);
IA i_den (p_den.first, p_den.second); else
IA i = i_num / i_den; {
return std::pair<double, double>(i.inf(), i.sup()); double s = nextafter (i, 0);
if (i < 0)
return std::pair<double, double>(i, s);
else
return std::pair<double, double>(s, i);
}
#else
mpfr_t y;
mpfr_init2 (y, 53); /* Assume IEEE-754 */
mpfr_set_z (y, x.backend().data(), GMP_RNDD);
double i = mpfr_get_d (y, GMP_RNDD); /* EXACT but can overflow */
mpfr_set_z (y, x.backend().data(), GMP_RNDU);
double s = mpfr_get_d (y, GMP_RNDU); /* EXACT but can overflow */
mpfr_clear (y);
return std::pair<double, double>(i, s);
#endif
} }
}; };
#endif
}; };
template <>
struct RET_boost_mp <boost::multiprecision::mpq_rational>
: RET_boost_mp_base <boost::multiprecision::mpq_rational> {
typedef boost::multiprecision::mpq_rational Type;
struct To_interval
: public CGAL::cpp98::unary_function< Type, std::pair< double, double > > {
std::pair<double, double>
operator()(const Type& x) const {
# if MPFR_VERSION_MAJOR >= 3
mpfr_exp_t emin = mpfr_get_emin();
mpfr_set_emin(-1073);
MPFR_DECL_INIT (y, 53); /* Assume IEEE-754 */
int r = mpfr_set_q (y, x.backend().data(), MPFR_RNDA);
r = mpfr_subnormalize (y, r, MPFR_RNDA); /* Round subnormals */
double i = mpfr_get_d (y, MPFR_RNDA); /* EXACT but can overflow */
mpfr_set_emin(emin); /* Restore old value, users may care */
// With mpfr_set_emax(1024) we could drop the is_finite test
if (r == 0 && is_finite (i))
return std::pair<double, double>(i, i);
else
{
double s = nextafter (i, 0);
if (i < 0)
return std::pair<double, double>(i, s);
else
return std::pair<double, double>(s, i);
}
# else
mpfr_t y;
mpfr_init2 (y, 53); /* Assume IEEE-754 */
mpfr_set_q (y, x.backend().data(), GMP_RNDD);
double i = mpfr_get_d (y, GMP_RNDD); /* EXACT but can overflow */
mpfr_set_q (y, x.backend().data(), GMP_RNDU);
double s = mpfr_get_d (y, GMP_RNDU); /* EXACT but can overflow */
mpfr_clear (y);
return std::pair<double, double>(i, s);
# endif
}
};
};
#endif
template <class Backend, boost::multiprecision::expression_template_option Eto>
struct Real_embeddable_traits<boost::multiprecision::number<Backend, Eto> >
: RET_boost_mp <boost::multiprecision::number<Backend, Eto> > {};
template <class T1,class T2,class T3,class T4,class T5>
struct Real_embeddable_traits<boost::multiprecision::detail::expression<T1,T2,T3,T4,T5> >
: Real_embeddable_traits<typename boost::multiprecision::detail::expression<T1,T2,T3,T4,T5>::result_type > {};
// Modular_traits // Modular_traits
@ -262,6 +350,13 @@ struct MT_boost_mp <T, boost::mpl::int_<boost::multiprecision::number_kind_integ
}; };
}; };
template <class Backend, boost::multiprecision::expression_template_option Eto>
struct Modular_traits<boost::multiprecision::number<Backend, Eto> >
: MT_boost_mp <boost::multiprecision::number<Backend, Eto> > {};
template <class T1,class T2,class T3,class T4,class T5>
struct Modular_traits<boost::multiprecision::detail::expression<T1,T2,T3,T4,T5> >
: Modular_traits<typename boost::multiprecision::detail::expression<T1,T2,T3,T4,T5>::result_type > {};
// Split_double // Split_double
template <class NT, class = boost::mpl::int_<boost::multiprecision::number_category<NT>::value> > template <class NT, class = boost::mpl::int_<boost::multiprecision::number_category<NT>::value> >
@ -284,6 +379,13 @@ struct SD_boost_mp <NT, boost::mpl::int_<boost::multiprecision::number_kind_inte
} }
}; };
template <class Backend, boost::multiprecision::expression_template_option Eto>
struct Split_double<boost::multiprecision::number<Backend, Eto> >
: SD_boost_mp <boost::multiprecision::number<Backend, Eto> > {};
template <class T1,class T2,class T3,class T4,class T5>
struct Split_double<boost::multiprecision::detail::expression<T1,T2,T3,T4,T5> >
: Split_double<typename boost::multiprecision::detail::expression<T1,T2,T3,T4,T5>::result_type > {};
// Fraction_traits // Fraction_traits
@ -335,27 +437,6 @@ struct FT_boost_mp <NT, boost::mpl::int_<boost::multiprecision::number_kind_rati
}; };
}; };
template <class Backend, boost::multiprecision::expression_template_option Eto>
struct Algebraic_structure_traits<boost::multiprecision::number<Backend, Eto> >
: AST_boost_mp <boost::multiprecision::number<Backend, Eto> > {};
template <class T1,class T2,class T3,class T4,class T5>
struct Algebraic_structure_traits<boost::multiprecision::detail::expression<T1,T2,T3,T4,T5> >
: Algebraic_structure_traits<typename boost::multiprecision::detail::expression<T1,T2,T3,T4,T5>::result_type > {};
template <class Backend, boost::multiprecision::expression_template_option Eto>
struct Real_embeddable_traits<boost::multiprecision::number<Backend, Eto> >
: RET_boost_mp <boost::multiprecision::number<Backend, Eto> > {};
template <class T1,class T2,class T3,class T4,class T5>
struct Real_embeddable_traits<boost::multiprecision::detail::expression<T1,T2,T3,T4,T5> >
: Real_embeddable_traits<typename boost::multiprecision::detail::expression<T1,T2,T3,T4,T5>::result_type > {};
template <class Backend, boost::multiprecision::expression_template_option Eto>
struct Modular_traits<boost::multiprecision::number<Backend, Eto> >
: MT_boost_mp <boost::multiprecision::number<Backend, Eto> > {};
template <class T1,class T2,class T3,class T4,class T5>
struct Modular_traits<boost::multiprecision::detail::expression<T1,T2,T3,T4,T5> >
: Modular_traits<typename boost::multiprecision::detail::expression<T1,T2,T3,T4,T5>::result_type > {};
template <class Backend, boost::multiprecision::expression_template_option Eto> template <class Backend, boost::multiprecision::expression_template_option Eto>
struct Fraction_traits<boost::multiprecision::number<Backend, Eto> > struct Fraction_traits<boost::multiprecision::number<Backend, Eto> >
: FT_boost_mp <boost::multiprecision::number<Backend, Eto> > {}; : FT_boost_mp <boost::multiprecision::number<Backend, Eto> > {};
@ -363,13 +444,6 @@ template <class T1,class T2,class T3,class T4,class T5>
struct Fraction_traits<boost::multiprecision::detail::expression<T1,T2,T3,T4,T5> > struct Fraction_traits<boost::multiprecision::detail::expression<T1,T2,T3,T4,T5> >
: Fraction_traits<typename boost::multiprecision::detail::expression<T1,T2,T3,T4,T5>::result_type > {}; : Fraction_traits<typename boost::multiprecision::detail::expression<T1,T2,T3,T4,T5>::result_type > {};
template <class Backend, boost::multiprecision::expression_template_option Eto>
struct Split_double<boost::multiprecision::number<Backend, Eto> >
: SD_boost_mp <boost::multiprecision::number<Backend, Eto> > {};
template <class T1,class T2,class T3,class T4,class T5>
struct Split_double<boost::multiprecision::detail::expression<T1,T2,T3,T4,T5> >
: Split_double<typename boost::multiprecision::detail::expression<T1,T2,T3,T4,T5>::result_type > {};
// Coercions // Coercions
@ -440,7 +514,7 @@ typename boost::multiprecision::detail::expression<T1,T2,T3,T4,T5>::result_type,
boost::multiprecision::number<B, E> > boost::multiprecision::number<B, E> >
{ }; { };
// TODO: coercion with expressions, fix existing coercions // TODO: fix existing coercions
// (double -> rational is implicit only for 1.56+, see ticket #10082) // (double -> rational is implicit only for 1.56+, see ticket #10082)
// The real solution would be to avoid specializing Coercion_traits for all pairs of number types and let it auto-detect what works, so only broken types need an explicit specialization. // The real solution would be to avoid specializing Coercion_traits for all pairs of number types and let it auto-detect what works, so only broken types need an explicit specialization.
@ -458,7 +532,13 @@ struct Coercion_traits<boost::multiprecision::number<B1, E1>, int> { \
}; \ }; \
template <class B1, boost::multiprecision::expression_template_option E1> \ template <class B1, boost::multiprecision::expression_template_option E1> \
struct Coercion_traits<int, boost::multiprecision::number<B1, E1> > \ struct Coercion_traits<int, boost::multiprecision::number<B1, E1> > \
: Coercion_traits<boost::multiprecision::number<B1, E1>, int> {} : Coercion_traits<boost::multiprecision::number<B1, E1>, int> {}; \
template <class T1, class T2, class T3, class T4, class T5> \
struct Coercion_traits<boost::multiprecision::detail::expression<T1,T2,T3,T4,T5>, int> \
: Coercion_traits<typename boost::multiprecision::detail::expression<T1,T2,T3,T4,T5>::result_type, int>{}; \
template <class T1, class T2, class T3, class T4, class T5> \
struct Coercion_traits<int, boost::multiprecision::detail::expression<T1,T2,T3,T4,T5> > \
: Coercion_traits<typename boost::multiprecision::detail::expression<T1,T2,T3,T4,T5>::result_type, int>{}
CGAL_COERCE_INT(short); CGAL_COERCE_INT(short);
CGAL_COERCE_INT(int); CGAL_COERCE_INT(int);
@ -479,7 +559,13 @@ struct Coercion_traits<boost::multiprecision::number<B1, E1>, float> { \
}; \ }; \
template <class B1, boost::multiprecision::expression_template_option E1> \ template <class B1, boost::multiprecision::expression_template_option E1> \
struct Coercion_traits<float, boost::multiprecision::number<B1, E1> > \ struct Coercion_traits<float, boost::multiprecision::number<B1, E1> > \
: Coercion_traits<boost::multiprecision::number<B1, E1>, float> {} : Coercion_traits<boost::multiprecision::number<B1, E1>, float> {}; \
template <class T1, class T2, class T3, class T4, class T5> \
struct Coercion_traits<boost::multiprecision::detail::expression<T1,T2,T3,T4,T5>, float> \
: Coercion_traits<typename boost::multiprecision::detail::expression<T1,T2,T3,T4,T5>::result_type, float>{}; \
template <class T1, class T2, class T3, class T4, class T5> \
struct Coercion_traits<float, boost::multiprecision::detail::expression<T1,T2,T3,T4,T5> > \
: Coercion_traits<typename boost::multiprecision::detail::expression<T1,T2,T3,T4,T5>::result_type, float>{}
CGAL_COERCE_FLOAT(float); CGAL_COERCE_FLOAT(float);
CGAL_COERCE_FLOAT(double); CGAL_COERCE_FLOAT(double);

11
Number_types/test/Number_types/utilities.cpp
View File

@ -10,6 +10,7 @@
#include <CGAL/Lazy_exact_nt.h> #include <CGAL/Lazy_exact_nt.h>
#include <CGAL/Interval_nt.h> #include <CGAL/Interval_nt.h>
#include <CGAL/Sqrt_extension.h> #include <CGAL/Sqrt_extension.h>
#include <CGAL/boost_mp.h>
#ifdef CGAL_USE_GMP #ifdef CGAL_USE_GMP
#include <CGAL/Gmpz.h> #include <CGAL/Gmpz.h>
@ -82,6 +83,12 @@ int main()
TESTIT(CGAL::Lazy_exact_nt<QMPF>, "Lazy_exact_nt<Quotient<MP_Float> >") TESTIT(CGAL::Lazy_exact_nt<QMPF>, "Lazy_exact_nt<Quotient<MP_Float> >")
TESTIT(CGAL::Interval_nt<>, "Interval_nt<>") TESTIT(CGAL::Interval_nt<>, "Interval_nt<>")
// Boost.Multiprecision
#ifdef CGAL_USE_BOOST_MP
TESTIT(boost::multiprecision::cpp_int, "cpp_int")
TESTIT(boost::multiprecision::cpp_rational, "cpp_rational")
#endif
// GMP based NTs // GMP based NTs
#ifdef CGAL_USE_GMP #ifdef CGAL_USE_GMP
TESTIT(CGAL::Gmpz, "Gmpz") TESTIT(CGAL::Gmpz, "Gmpz")
@ -90,6 +97,10 @@ int main()
TESTIT(CGAL::Mpzf, "Mpzf") TESTIT(CGAL::Mpzf, "Mpzf")
# endif # endif
TESTIT(CGAL::Gmpq, "Gmpq") TESTIT(CGAL::Gmpq, "Gmpq")
# ifdef CGAL_USE_BOOST_MP
TESTIT(boost::multiprecision::mpz_int, "mpz_int")
TESTIT(boost::multiprecision::mpq_rational, "mpq_rational")
# endif
#endif // CGAL_USE_GMP #endif // CGAL_USE_GMP
#ifdef CGAL_USE_GMPXX #ifdef CGAL_USE_GMPXX
TESTIT(mpz_class, "mpz_class") TESTIT(mpz_class, "mpz_class")