Changeset 66 for include/fflas-ffpack/modular-balanced.h

Timestamp:

06/05/08 14:58:48 (6 months ago)

Author:

pernet

Message:

Update the randiter:
* add the non-zero randiter class (following LinBox? structure)
* modify finite fields accordingly
* Update implementation of finite fields to fix bugs, improve efficiency and readability
* minor tweeks

Files:

• include/fflas-ffpack/modular-balanced.h (modified) (16 diffs)

include/fflas-ffpack/modular-balanced.h

@@ -17 +17 @@
 #include <math.h>
 #include "fflas-ffpack/modular-randiter.h"
+#include "fflas-ffpack/nonzero-randiter.h"
 
 template <class Element>
@@ -26 +27 @@
 protected:
         double modulus;
-        double inv_modulus;
         double half_mod;
-
+        unsigned long lmodulus;
        
 public:        
@@ -34 +34 @@
         typedef unsigned long FieldInt;
         typedef ModularBalancedRandIter<double> RandIter;
+        typedef NonzeroRandIter<ModularBalanced<double>, ModularBalancedRandIter<double> >  NonZeroRandIter;
 
         const bool balanced;
 
-        ModularBalanced<double> (int p)  : modulus((double)p), inv_modulus(1./(double)p), half_mod( (p-1.)/2), balanced(true) {}
-
-        ModularBalanced<double>(const ModularBalanced<double>& mf) : modulus(mf.modulus), inv_modulus(mf.inv_modulus), half_mod(mf.half_mod), balanced(true){}
+        ModularBalanced<double> (int p) : modulus((double)p),
+                half_mod( (p-1.)/2),
+                balanced(true) {}
+
+        ModularBalanced<double>(const ModularBalanced<double>& mf)
+        : modulus(mf.modulus), half_mod(mf.half_mod), balanced(true){}
         
         const ModularBalanced<double> &operator=(const ModularBalanced<double> &F) {
                         modulus = F.modulus;
-                        inv_modulus = F.inv_modulus;
                         half_mod = F.half_mod;
                         return *this;
@@ -70 +73 @@
                 
         std::ostream &write (std::ostream &os) const {
-                return os << "int mod " << (int)modulus;
+                return os << "double mod " << (int)modulus;
         }
         
@@ -96 +99 @@
         }
 
-        inline Element& init(Element& x, double y =0) const {             
+        inline Element& init(Element& x, const double y =0) const {               
                 double tmp;
                 
@@ -131 +134 @@
                 if ( x > half_mod ) return x -= modulus;
                 if ( x < -half_mod ) return x += modulus; 
-                else return x;
+                return x;
         }
  
@@ -137 +140 @@
                 x = y - z;
                 if (x > half_mod) return x -= modulus;
-                else if (x < -half_mod) return x += modulus;
-                else return x;
+                if (x < -half_mod) return x += modulus;
+                return x;
         }
                 
         inline Element &mul (Element &x, const Element &y, const Element &z) const {            
                 double tmp= y*z;
-                //x= tmp - floor(tmp*inv_modulus)*modulus;
                 return init (x, tmp);
         }
@@ -187 +189 @@
                               const Element &y) const {
                 double tmp= a*x+y;
-                //return r = tmp- floor(tmp*inv_modulus)*modulus; 
                 return init( r, tmp);
-
         }
 
@@ -196 +196 @@
                 if ( x > half_mod ) return x -= modulus;
                 if ( x < -half_mod ) return x += modulus; 
-                else return x;
+                return x;
         }
  
@@ -203 +203 @@
                 if ( x > half_mod ) return x -= modulus;
                 if ( x < -half_mod ) return x += modulus; 
-                else return x;
+                return x;
         }
  
@@ -224 +224 @@
                 
         inline Element &axpyin (Element &r, const Element &a, const Element &x) const {
-                double tmp=r+a*x;
-                //return r= tmp- floor(tmp*inv_modulus)*modulus; 
-                return init( r, tmp );
+                r += a * x;
+                return init( r, r);
         }
 
@@ -238 +237 @@
 protected:
         float modulus;
-        float inv_modulus;
         float half_mod;
 
@@ -246 +244 @@
         typedef unsigned long FieldInt;
         typedef ModularBalancedRandIter<float> RandIter;
+        typedef NonzeroRandIter<ModularBalanced<float>, RandIter> NonZeroRandIter;
 
         const bool balanced;
 
-        ModularBalanced<float> (int p)  : modulus((float)p), inv_modulus(1./(float)p), half_mod( (p-1.)/2), balanced(true) {}
-
-        ModularBalanced<float>(const ModularBalanced<float>& mf) : modulus(mf.modulus), inv_modulus(mf.inv_modulus), half_mod(mf.half_mod), balanced(true){}
+        ModularBalanced<float> (int p)  : modulus((float)p), half_mod( (p-1.)/2), balanced(true) {}
+
+        ModularBalanced<float>(const ModularBalanced<float>& mf) : modulus(mf.modulus), half_mod(mf.half_mod), balanced(true){}
         
         const ModularBalanced<float> &operator=(const ModularBalanced<float> &F) {
                         modulus = F.modulus;
-                        inv_modulus = F.inv_modulus;
                         half_mod = F.half_mod;
                         return *this;
@@ -361 +359 @@
                 
         inline Element &mul (Element &x, const Element &y, const Element &z) const {            
-                float tmp= y*z;
-                //x= tmp - floor(tmp*inv_modulus)*modulus;
-                return init (x, tmp);
+                x = y*z;
+                return init (x, x);
         }
  
         inline Element &div (Element &x, const Element &y, const Element &z) const {
-                Element temp;
-                inv (temp, z);
-                return mul (x, y, temp);
+                inv (x, z);
+                return mulin (x, y);
         }
  
@@ -405 +401 @@
                               const Element &x, 
                               const Element &y) const {
-                float tmp= a*x+y;
-                //return r = tmp- floor(tmp*inv_modulus)*modulus; 
-                return init( r, tmp);
+                r = a * x + y;
+                return init( r, r);
 
         }
@@ -443 +438 @@
                 
         inline Element &axpyin (Element &r, const Element &a, const Element &x) const {
-                float tmp=r+a*x;
-                //return r= tmp- floor(tmp*inv_modulus)*modulus; 
-                return init( r, tmp );
+                r += a * x;
+                return init( r, r);
         }
         static inline float getMaxModulus()