Merge pull request #1572 from luis4a0/Algebraic_kernel_d-add_assertions-Luis

Check for the MPFR version in RS-based AK

Algebraic_kernel_d/include/CGAL/RS/algebraic_1.h

@@ -192,12 +192,13 @@ boost::totally_ordered<Algebraic_1<Polynomial_,
 #else
 #define CGAL_RS_DBL_PREC        53
 #endif
-        double to_double(){
+        // This function is const because left and right are mutable.
+        double to_double()const{
                 typedef Real_embeddable_traits<Bound>                   RT;
                 typedef typename RT::To_double                          TD;
-                Refiner()(get_pol(),get_left(),get_right(),CGAL_RS_DBL_PREC);
+                Refiner()(pol,left,right,CGAL_RS_DBL_PREC);
-                CGAL_assertion(TD()(get_left())==TD()(get_right()));
+                CGAL_assertion(TD()(left)==TD()(right));
-                return TD()(get_left());
+                return TD()(left);
         }
         std::pair<double,double> to_interval()const{
                 typedef Real_embeddable_traits<Bound>                   RT;
@@ -266,7 +267,7 @@ public INTERN_RET::Real_embeddable_traits_base<
 
         class To_double:public std::unary_function<Type,double>{
                 public:
-                double operator()(Type a)const{return a.to_double();}
+                double operator()(const Type &a)const{return a.to_double();}
         };
 
         class To_interval:

Algebraic_kernel_d/include/CGAL/RS/algebraic_z_1.h

@@ -184,7 +184,7 @@ boost::totally_ordered<Algebraic_z_1<Polynomial_,
 #else
 #define CGAL_RS_DBL_PREC        53
 #endif
-        double to_double(){
+        double to_double()const{
                 typedef Real_embeddable_traits<Bound>                   RT;
                 typedef typename RT::To_double                          TD;
                 ZRefiner()(get_zpol(),
@@ -258,7 +258,7 @@ public INTERN_RET::Real_embeddable_traits_base<
 
         class To_double:public std::unary_function<Type,double>{
                 public:
-                double operator()(Type a)const{return a.to_double();}
+                double operator()(const Type &a)const{return a.to_double();}
         };
 
         class To_interval:

Algebraic_kernel_d/include/CGAL/RS/bisection_refiner_1.h

@@ -47,7 +47,7 @@ operator()(const Polynomial_&,Bound_&,Bound_&,int){
 
 // This works with any type of polynomial, but only for Gmpfr bounds.
 // TODO: Beyond writing generically, optimize this function. This would
-// remove the need for the next function, which essentially the same.
+// remove the need for the next function, which is essentially the same.
 template<>
 void
 Bisection_refiner_1<Polynomial<Gmpz>,Gmpfr>::
@@ -58,8 +58,6 @@ operator()(const Polynomial<Gmpz> &pol,Gmpfr &left,Gmpfr &right,int prec){
         typedef CGAL::RS_AK1::Signat_1<Polynomial,Gmpfr>
                                                         Signat;
         CGAL_precondition(left<=right);
-        // TODO: add precondition to check whether the interval is a point
-        // or the evaluations on its endpoints have different signs
         //std::cout<<"refining ["<<left<<","<<right<<"]"<<std::endl;
 
         CGAL_RS_GMPFR_MAKE_UNIQUE(left,temp_left);
@@ -72,6 +70,7 @@ operator()(const Polynomial<Gmpz> &pol,Gmpfr &left,Gmpfr &right,int prec){
         mpfr_t center;
 
         sl=signof(left);
+        CGAL_precondition(sl!=signof(right)||(left==right&&sl==ZERO));
         if(sl==ZERO)
                 return;
         pl=left.get_precision();
@@ -120,8 +119,6 @@ operator()(const Polynomial<Gmpq> &pol,Gmpfr &left,Gmpfr &right,int prec){
         typedef CGAL::RS_AK1::Signat_1<Polynomial,Gmpfr>
                                                         Signat;
         CGAL_precondition(left<=right);
-        // TODO: add precondition to check whether the interval is a point
-        // or the evaluations on its endpoints have different signs
         //std::cout<<"refining ["<<left<<","<<right<<"]"<<std::endl;
 
         CGAL_RS_GMPFR_MAKE_UNIQUE(left,temp_left);
@@ -134,6 +131,7 @@ operator()(const Polynomial<Gmpq> &pol,Gmpfr &left,Gmpfr &right,int prec){
         mpfr_t center;
 
         sl=signof(left);
+        CGAL_precondition(sl!=signof(right)||(left==right&&sl==ZERO));
         if(sl==ZERO)
                 return;
         pl=left.get_precision();

Algebraic_kernel_d/include/CGAL/RS/rs23_k_isolator_1.h

@@ -80,6 +80,7 @@ RS23_k_isolator_1(const CGAL::Polynomial<CGAL::Gmpz> &p):_polynomial(p){
         for(int j=0;j<numsols;++j){
                 Gmpfr left(intervals[j].inf());
                 Gmpfr right(intervals[j].sup());
+                CGAL_assertion(left<=right);
                 KRefiner()(p,left,right,53);
                 _real_roots.push_back(intervals[j]);
         }

Algebraic_kernel_d/include/CGAL/RS/rs2_calls.h

@@ -33,6 +33,19 @@
 #define CGALRS_PTR(a)   void *a
 #endif
 
+// RS3 does not work with MPFR 3.1.3 to 3.1.5. In case RS3 is enabled and
+// the version of MPFR is one of those buggy versions, abort the compilation
+// and instruct the user to update MPFR or don't use RS3.
+#ifdef CGAL_USE_RS3
+#include <boost/static_assert.hpp>
+BOOST_STATIC_ASSERT_MSG(
+        MPFR_VERSION_MAJOR!=3 ||
+        MPFR_VERSION_MINOR!=1 ||
+        MPFR_VERSION_PATCHLEVEL<3 || MPFR_VERSION_PATCHLEVEL>5,
+        "RS3 does not work with MPFR versions 3.1.3 to 3.1.5. "#
+        "Please update MPFR or disable RS3.");
+#endif // CGAL_USE_RS3
+
 namespace CGAL{
 namespace RS2{
 
@@ -120,6 +133,7 @@ struct RS2_calls{
                         // Construct Gmpfr's with pointers to endpoints.
                         Gmpfr left(&(root_pointer->left),root_prec);
                         Gmpfr right(&(root_pointer->right),root_prec);
+                        CGAL_assertion(left<=right);
                         // Copy them, to have the data out of RS memory.
                         *x++=Gmpfi(left,right,root_prec+1);
                         ident_node=rs_export_list_vect_ibfr_nextnode

Algebraic_kernel_d/include/CGAL/RS/rs3_k_refiner_1.h

@@ -22,10 +22,16 @@
 #define CGAL_RS_RS3_K_REFINER_1_H
 
 #include <CGAL/Polynomial_traits_d.h>
+#include "polynomial_converter_1.h"
 #include "rs2_calls.h"
 #include <rs3_fncts.h>
 #include "Gmpfr_make_unique.h"
 
+// If we want assertions, we need to evaluate.
+#ifndef CGAL_NO_PRECONDITIONS
+#include "signat_1.h"
+#endif
+
 namespace CGAL{
 namespace RS3{
 
@@ -51,8 +57,15 @@ operator()
         typedef Polynomial_traits_d<Polynomial>         Ptraits;
         typedef Ptraits::Degree                         Degree;
         CGAL_precondition(left<=right);
-        // TODO: add precondition to check whether the interval is a point
+#ifndef CGAL_NO_PRECONDITIONS
-        // or the evaluations on its endpoints have different signs
+        typedef Ptraits::Make_square_free               Sfpart;
+        typedef CGAL::RS_AK1::Signat_1<Polynomial,Gmpfr>
+                                                        Signat;
+        Polynomial sfpp=Sfpart()(pol);
+        Signat signof(sfpp);
+        CGAL::Sign sl=signof(left);
+        CGAL_precondition(sl!=signof(right)||(left==right&&sl==ZERO));
+#endif
         //std::cout<<"refining ["<<left<<","<<right<<"]"<<std::endl;
         int deg=Degree()(pol);
         mpz_t* coefficients=(mpz_t*)malloc((deg+1)*sizeof(mpz_t));
@@ -98,12 +111,21 @@ operator()
         typedef Polynomial_traits_d<ZPolynomial>        ZPtraits;
         typedef ZPtraits::Degree                        ZDegree;
         CGAL_precondition(left<=right);
-        // TODO: add precondition to check whether the interval is a point
+#ifndef CGAL_NO_PRECONDITIONS
-        // or the evaluations on its endpoints have different signs
+        typedef ZPtraits::Make_square_free              ZSfpart;
+        typedef CGAL::RS_AK1::Signat_1<ZPolynomial,Gmpfr>
+                                                        Signat;
+#endif
         //std::cout<<"refining ["<<left<<","<<right<<"]"<<std::endl;
         Polynomial<Gmpz> zpol=CGAL::RS_AK1::Polynomial_converter_1<
                                         CGAL::Polynomial<Gmpq>,
                                         CGAL::Polynomial<Gmpz> >()(qpol);
+#ifndef CGAL_NO_PRECONDITIONS
+        ZPolynomial zsfpp=ZSfpart()(zpol);
+        Signat signof(zsfpp);
+        CGAL::Sign sl=signof(left);
+        CGAL_precondition(sl!=signof(right)||(left==right&&sl==ZERO));
+#endif
         int deg=ZDegree()(zpol);
         mpz_t* coefficients=(mpz_t*)malloc((deg+1)*sizeof(mpz_t));
         __mpfi_struct interval;

Algebraic_kernel_d/include/CGAL/RS/rs3_refiner_1.h

@@ -26,6 +26,11 @@
 #include <rs3_fncts.h>
 #include "Gmpfr_make_unique.h"
 
+// If we want assertions, we need to evaluate.
+#ifndef CGAL_NO_PRECONDITIONS
+#include "signat_1.h"
+#endif
+
 namespace CGAL{
 namespace RS3{
 
@@ -51,8 +56,15 @@ operator()
         typedef Polynomial_traits_d<Polynomial>         Ptraits;
         typedef Ptraits::Degree                         Degree;
         CGAL_precondition(left<=right);
-        // TODO: add precondition to check whether the interval is a point
+#ifndef CGAL_NO_PRECONDITIONS
-        // or the evaluations on its endpoints have different signs
+        typedef Ptraits::Make_square_free               Sfpart;
+        typedef CGAL::RS_AK1::Signat_1<Polynomial,Gmpfr>
+                                                        Signat;
+        Polynomial sfpp=Sfpart()(pol);
+        Signat signof(sfpp);
+        CGAL::Sign sl=signof(left);
+        CGAL_precondition(sl!=signof(right)||(left==right&&sl==ZERO));
+#endif
         //std::cout<<"refining ["<<left<<","<<right<<"]"<<std::endl;
         int deg=Degree()(pol);
         mpz_t* coefficients=(mpz_t*)malloc((deg+1)*sizeof(mpz_t));
@@ -104,13 +116,22 @@ operator()
         typedef Polynomial_traits_d<ZPolynomial>        ZPtraits;
         typedef ZPtraits::Degree                        ZDegree;
         CGAL_precondition(left<=right);
-        // TODO: add precondition to check whether the interval is a point
+#ifndef CGAL_NO_PRECONDITIONS
-        // or the evaluations on its endpoints have different signs
+        typedef ZPtraits::Make_square_free              ZSfpart;
+        typedef CGAL::RS_AK1::Signat_1<ZPolynomial,Gmpfr>
+                                                        Signat;
+#endif
         //std::cout<<"refining ["<<left<<","<<right<<"]"<<std::endl;
         // Construct a Gmpz polynomial from the original Gmpq polynomial.
         Polynomial<Gmpz> zpol=CGAL::RS_AK1::Polynomial_converter_1<
                                         CGAL::Polynomial<Gmpq>,
                                         CGAL::Polynomial<Gmpz> >()(qpol);
+#ifndef CGAL_NO_PRECONDITIONS
+        ZPolynomial zsfpp=ZSfpart()(zpol);
+        Signat signof(zsfpp);
+        CGAL::Sign sl=signof(left);
+        CGAL_precondition(sl!=signof(right)||(left==right&&sl==ZERO));
+#endif
         int deg=ZDegree()(zpol);
         mpz_t* coefficients=(mpz_t*)malloc((deg+1)*sizeof(mpz_t));
         __mpfi_struct interval;