use boost::ordered_field_operators1

modular arithmetic needs ieee double precision

Modular_arithmetic/include/CGAL/Modular_arithmetic/Modular_type.h

@@ -27,14 +27,19 @@
 
 CGAL_BEGIN_NAMESPACE
 
+namespace CGALi{
+
+struct Modular_arithmetic_needs_ieee_double_precision{
+    Modular_arithmetic_needs_ieee_double_precision(){
+        CGAL::force_ieee_double_precision();
+    }
+};
+}
+
 class Modular;
     
 Modular operator + (const Modular&);
 Modular operator - (const Modular&);
-Modular operator + (const Modular&, const Modular&);
-Modular operator - (const Modular&, const Modular&);
-Modular operator * (const Modular&, const Modular&);
-Modular operator / (const Modular&, const Modular&);
 
 std::ostream& operator << (std::ostream& os, const Modular& p);
 std::istream& operator >> (std::istream& is, Modular& p);
@@ -50,7 +55,9 @@ std::istream& operator >> (std::istream& is, Modular& p);
  * 
  * \see Modular_traits
  */
-class Modular{
+class Modular:
+    boost::ordered_field_operators1< Modular,
+    boost::ordered_field_operators2< Modular, int > >{
     
 public:
     typedef Modular Self;
@@ -58,6 +65,8 @@ public:
 
 private:
     static const double  CST_CUT; 
+    static const CGALi::Modular_arithmetic_needs_ieee_double_precision 
+    modular_arithmetic_needs_ieee_double_precision;
 private:
    
     static int prime_int;
@@ -201,28 +210,40 @@ public:
         x() = MOD_add(x(),p2.x()); 
         return (*this); 
     }
-
     Self& operator -= (const Self& p2){ 
         x() = MOD_add(x(),MOD_negate(p2.x())); 
         return (*this); 
     }
-
     Self& operator *= (const Self& p2){ 
         x() = MOD_mul(x(),p2.x()); 
         return (*this); 
     }
-
     Self& operator /= (const Self& p2) { 
         x() = MOD_div(x(),p2.x()); 
         return (*this); 
     }
+    // 
+    Self& operator += (int p2) { 
+        x() = MOD_add(x(),Modular(p2).x()); 
+        return (*this); 
+    }
+    Self& operator -= (int p2){ 
+        x() = MOD_add(x(),Modular(-p2).x()); 
+        return (*this); 
+    }
+
+    Self& operator *= (int p2){ 
+        x() = MOD_mul(x(),Modular(p2).x()); 
+        return (*this); 
+    }
+
+    Self& operator /= (int p2) { 
+        x() = MOD_div(x(),Modular(p2).x()); 
+        return (*this); 
+    }
   
     friend Self operator + (const Self&);
     friend Self operator - (const Self&);                
-    friend Self operator + (const Self&, const Self&);  
-    friend Self operator - (const Self&, const Self&);  
-    friend Self operator * (const Self&, const Self&);    
-    friend Self operator / (const Self& p1, const Self& p2);
 };
 
 inline Modular operator + (const Modular& p1)
@@ -235,68 +256,17 @@ inline Modular operator - (const Modular& p1){
     return r; 
 }
 
-inline Modular operator + (const Modular& p1,const Modular& p2) { 
-    typedef Modular MOD; 
-    Modular r; 
-    r.x() = MOD::MOD_add(p1.x(),p2.x()); 
-    return r; 
-} 
-
-inline Modular operator - (const Modular& p1, const Modular& p2) { 
-    return p1+(-p2);  
-}
-
-inline Modular operator * (const Modular& p1, const Modular& p2) { 
-    typedef Modular MOD;
-    Modular r;
-    r.x() = MOD::MOD_mul(p1.x(),p2.x()); 
-    return r;
-}
-
-inline Modular operator / (const Modular& p1, const Modular& p2) { 
-    typedef Modular MOD; 
-    Modular r;
-    r.x() = MOD::MOD_div(p1.x(),p2.x()); 
-    return r;
-}
-
 inline bool operator == (const Modular& p1, const Modular& p2)
 { return ( p1.x()==p2.x() ); }   
+inline bool operator == (const Modular& p1, int p2)
+{ return ( p1 == Modular(p2) ); }   
 
-inline bool operator != (const Modular& p1, const Modular& p2)
-{ return ( p1.x()!=p2.x() ); }
 
-// left hand side
+inline bool operator < (const Modular& p1, const Modular& p2)
-inline bool operator == (int num, const Modular& p) 
+{ return ( p1.x() < p2.x() ); }   
-{ return ( Modular(num) == p );}
+inline bool operator < (const Modular& p1, int p2)
-inline bool operator != (int num, const Modular& p) 
+{ return ( p1.x() < Modular(p2).x() ); }   
-{ return ( Modular(num) != p );}
 
-// right hand side
-inline bool operator == (const Modular& p, int num) 
-{ return ( Modular(num) == p );}
-inline bool operator != (const Modular& p, int num) 
-{ return ( Modular(num) != p );} 
-
-// left hand side
-inline Modular operator + (int num, const Modular& p2)
-{ return (Modular(num) + p2); }
-inline Modular operator - (int num, const Modular& p2)
-{ return (Modular(num) - p2); }
-inline Modular operator * (int num, const Modular& p2)
-{ return (Modular(num) * p2); }
-inline Modular operator / (int num, const Modular& p2)
-{ return (Modular(num)/p2); }
-
-// right hand side
-inline Modular operator + (const Modular& p1, int num)
-{ return (p1 + Modular(num)); }
-inline Modular operator - (const Modular& p1, int num)
-{ return (p1 - Modular(num)); }
-inline Modular operator * (const Modular& p1, int num)
-{ return (p1 * Modular(num)); }
-inline Modular operator / (const Modular& p1, int num)
-{ return (p1 / Modular(num)); }
 
 // I/O 
 inline std::ostream& operator << (std::ostream& os, const Modular& p) {   
@@ -305,7 +275,6 @@ inline std::ostream& operator << (std::ostream& os, const Modular& p) {
     return os;
 }
 
-
 inline std::istream& operator >> (std::istream& is, Modular& p) {
     typedef Modular MOD;
     char ch;

Modular_arithmetic/src/CGAL/Modular_type.cpp

@@ -27,10 +27,13 @@
 #include <CGAL/Modular_arithmetic/Modular_type.h>
 
 namespace CGAL{
-    int Modular::prime_int = 67111067;
-    double Modular::prime =67111067.0;
-    double Modular::prime_inv =1/67111067.0;
 
-    const double Modular::CST_CUT = std::ldexp( 3., 51 );
+int Modular::prime_int = 67111067;
+double Modular::prime =67111067.0;
+double Modular::prime_inv =1/67111067.0;
+const double Modular::CST_CUT = std::ldexp( 3., 51 );
 
+const CGALi::Modular_arithmetic_needs_ieee_double_precision 
+Modular::modular_arithmetic_needs_ieee_double_precision 
+= CGALi::Modular_arithmetic_needs_ieee_double_precision();
 }

Modular_arithmetic/test/Modular_arithmetic/Modular.cpp

@@ -9,6 +9,8 @@
 #include <CGAL/Modular.h>
 
 #include <CGAL/Test/_test_algebraic_structure.h>
+#include <CGAL/Arithmetic_kernel.h>
+#include <CGAL/number_utils.h>
 
 int main()
 {   
@@ -76,6 +78,39 @@ int main()
 
     //cout << x << endl;
 
+    assert(7 == NT::set_current_prime(old_prime));
+
+    typedef CGAL::Arithmetic_kernel AK;
+    typedef AK::Integer Integer; 
+    
+    Integer       int_x(7);
+    Integer       prime(NT::get_current_prime());
+    CGAL::Modular mod_x(7);
+    for(int i = 0; i < 10000; i++){
+        assert(mod_x == CGAL::modular_image(int_x));
+        int_x *= int_x; int_x = CGAL::mod(int_x, prime);
+        mod_x *= mod_x; 
+        assert(mod_x == CGAL::modular_image(int_x));
+        int_x += int_x; int_x = CGAL::mod(int_x, prime);
+        mod_x += mod_x; 
+        assert(mod_x == CGAL::modular_image(int_x));
+
+        assert(mod_x == CGAL::modular_image(int_x));
+        int_x *= int_x; int_x = CGAL::mod(int_x, prime);
+        mod_x *= mod_x; 
+
+        assert(mod_x == CGAL::modular_image(int_x));
+        int_x -= int_x; int_x = CGAL::mod(int_x, prime);
+        mod_x -= mod_x; 
+    }
+    //  CGAL::force_ieee_double_precision();
+    {
+        CGAL::Modular::set_current_prime(67111043);
+        CGAL::Modular x(-33546401);
+        CGAL::Modular y(23950928);
+        CGAL::Modular q = CGAL::integral_division(x,y);
+        assert(x == q*y);
+    }
 }