Changeset 66

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

Location:

include/fflas-ffpack

Files:

• fflas.h (modified) (2 diffs)
• fflas_fgemm.inl (modified) (1 diff)
• ffpack_frobenius.inl (modified) (1 diff)
• modular-balanced.h (modified) (16 diffs)
• modular-int.h (modified) (3 diffs)
• modular-positive.h (modified) (4 diffs)
• modular-randiter.h (modified) (3 diffs)
• nonzero-randiter.h (added)

include/fflas-ffpack/fflas.h

@@ -26 +26 @@
 #include "linbox/field/modular-double.h"
 #include "linbox/field/modular-float.h"
+#include "linbox/field/modular-balanced-double.h"
+#include "linbox/field/modular-balanced-float.h"
 namespace LinBox {
 #else
@@ -35 +37 @@
         
 #ifndef __LINBOX_STRASSEN_OPTIMIZATION
-#define WINOTHRESHOLD 2000
+#define WINOTHRESHOLD 1000
 #else
 #define WINOTHRESHOLD __LINBOX_WINOTHRESHOLD

include/fflas-ffpack/fflas_fgemm.inl

@@ -95 +95 @@
                                        Bdd, dldb, betad, Cd, n, kmax,base );
 
-                        typename Field::Element tmp;
                         MatD2MatF (F, C, ldc, Cd, n, m, n);
                         MatF2MatD (F, Cd, n, C, ldc, m, n);

include/fflas-ffpack/ffpack_frobenius.inl

@@ -43 +43 @@
         // Picking a random noc x N block vector U^T
         typename Field::RandIter g (F);
+        typename Field::NonZeroRandIter nzg (F,g);
         for (size_t i = 0; i < noc; ++i)
                 for (size_t j = 0; j < N; ++j)
                         g.random( *(K + i*ldk +j) );
         for (size_t i = 0; i < noc; ++i)
-                g.nonzerorandom (*(K + i*ldk +i));
+                nzg.random (*(K + i*ldk +i));
         
         // Computing the bloc Krylov matrix [U AU .. A^(c-1) U]^T

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()

include/fflas-ffpack/modular-int.h

@@ -18 +18 @@
 #include <sys/time.h>
 #include "fflas-ffpack/modular-randiter.h"
+#include "fflas-ffpack/nonzero-randiter.h"
 template< class Element >
 class Modular;
@@ -26 +27 @@
         typedef unsigned long long uint64;
         typedef unsigned long FieldInt;
-
+        
 protected:
         
-                int modulus;
-
-                double modulusinv;
-                
-        
-                int _two64;
+        int modulus;
+        double modulusinv;
+        int _two64;
                 
         public:        
                 
-                typedef int Element;
-                typedef ModularRandIter<int> RandIter;
-
-                const bool balanced;
-
-                
-                Modular (int value)  : modulus(value), balanced(false) {
-                        modulusinv = 1 / ((double) value); 
-                        _two64 = (int) ((uint64) (-1) % (uint64) value);
-                        _two64 += 1;
-                        if (_two64 >= value) _two64 -= value;
-                }
-
-                Modular(const Modular<int>& mf) : modulus(mf.modulus),modulusinv(mf.modulusinv),_two64(mf._two64), balanced(false){}
-
-                const Modular &operator=(const Modular<int> &F) {
-                        modulus = F.modulus;
-                        modulusinv = F.modulusinv;
-                        _two64 = F._two64;
-                        return *this;
-                }
-
-        
-                FieldInt &cardinality (FieldInt &c) const{ 
-                        return c = modulus;
-                }
-
-                FieldInt &characteristic (FieldInt &c) const {
-                        return c = modulus; 
-                }
-
-                FieldInt &convert (FieldInt &x, const Element &y) const { 
-                        return x = y;
-                }
+        typedef int Element;
+        typedef ModularRandIter<int> RandIter;
+        typedef NonzeroRandIter<Modular<int>,RandIter> NonZeroRandIter;
+        
+        const bool balanced;
+        
+        
+        Modular (int value)  : modulus(value), balanced(false) {
+                modulusinv = 1 / ((double) value); 
+                _two64 = (int) ((uint64) (-1) % (uint64) value);
+                _two64 += 1;
+                if (_two64 >= value) _two64 -= value;
+        }
+        
+        Modular(const Modular<int>& mf) : modulus(mf.modulus),modulusinv(mf.modulusinv),_two64(mf._two64), balanced(false){}
+        
+        const Modular &operator=(const Modular<int> &F) {
+                modulus = F.modulus;
+                modulusinv = F.modulusinv;
+                _two64 = F._two64;
+                return *this;
+        }
+
+        
+        FieldInt &cardinality (FieldInt &c) const{ 
+                return c = modulus;
+        }
+
+        FieldInt &characteristic (FieldInt &c) const {
+                return c = modulus; 
+        }
+
+        FieldInt &convert (FieldInt &x, const Element &y) const { 
+                return x = y;
+        }
         Element &convert (Element &x, const Element &y) const { 
                 return x = y;
         }
 
-                double & convert (double &x, const Element &y) const { 
-                        return x = (double) y;
-                }
+        double & convert (double &x, const Element &y) const { 
+                return x = (double) y;
+        }
 
         float & convert (float &x, const Element &y) const { 
@@ -84 +83 @@
         }
                 
-                std::ostream &write (std::ostream &os) const {
-                        return os << "int mod " << modulus;
+        std::ostream &write (std::ostream &os) const {
+                return os << "int mod " << modulus;
+        }
+                
+        std::istream &read (std::istream &is) {
+                is >> modulus; 
+                modulusinv = 1 /((double) modulus );
+                _two64 = (int) ((uint64) (-1) % (uint64) modulus);
+                _two64 += 1;
+                if (_two64 >= modulus) _two64 -= modulus;
+                        
+                return is;
+        }
+                
+        std::ostream &write (std::ostream &os, const Element &x) const {
+                return os << x;
+        }
+
+        std::istream &read (std::istream &is, Element &x) const {
+                FieldInt tmp;
+                is >> tmp;
+                init(x,tmp); 
+                return is;
+        }
+                
+
+        template<class Element1>
+        Element &init (Element & x, const Element1 &y) const {
+                x = y % modulus;
+                if (x < 0) x += modulus;
+                return x;
+        }
+
+        Element &init (Element &x, const double &y) const  { 
+                double z = fmod(y, (double)modulus);
+                if (z < 0) z += (double)modulus;
+                z += 0.5;
+                return x = static_cast<long>(z); //rounds towards 0
+        }
+        Element &init (Element &x, const float  &y) const  { 
+                float z = fmod(y, (float)modulus);
+                if (z < 0) z += (float)modulus;
+                z += 0.5;
+                return x = static_cast<long>(z); //rounds towards 0
+        }
+
+        Element &init (Element &x, const FieldInt &y) const  {
+                x = y % modulus;
+                if (x < 0) x += modulus;
+                return x;
+        }
+
+        inline Element& init(Element& x, int y =0) const {
+                x = y % modulus;
+                if ( x < 0 ) x += modulus;
+                return x;
+        }
+
+        inline Element& init(Element& x, long y) const {
+                x = y % modulus;
+                if ( x < 0 ) x += modulus; 
+                return x;
+        }
+                
+        inline Element& assign(Element& x, const Element& y) const {
+                return x = y;
+        }
+                                                                        
+                
+        inline bool areEqual (const Element &x, const Element &y) const {
+                return x == y;
+        }
+
+        inline  bool isZero (const Element &x) const {
+                return x == 0; 
+        }
+                
+        inline bool isOne (const Element &x) const {
+                return x == 1; 
+        }
+
+        inline Element &add (Element &x, const Element &y, const Element &z) const {
+                x = y + z;
+                if ( x >= modulus ) x -= modulus;
+                return x;
+        }
+ 
+        inline Element &sub (Element &x, const Element &y, const Element &z) const {
+                x = y - z;
+                if (x < 0) x += modulus;
+                return x;
+        }
+                
+        inline Element &mul (Element &x, const Element &y, const Element &z) const {
+                int q;
+
+                q  = (int) ((((double) y)*((double) z)) * modulusinv);  // q could be off by (+/-) 1
+                x = (int) (y*z - q*modulus);
+                        
+                        
+                if (x >= modulus)
+                        x -= modulus;
+                else if (x < 0)
+                        x += modulus;
+
+                return x;
+        }
+ 
+        inline Element &div (Element &x, const Element &y, const Element &z) const {
+                Element temp;
+                inv (temp, z);
+                return mul (x, y, temp);
+        }
+ 
+        inline Element &neg (Element &x, const Element &y) const {
+                if(y == 0) return x=0;
+                else return x = modulus-y;
+        }
+ 
+        inline Element &inv (Element &x, const Element &y) const {
+                int d, t;                       
+                XGCD(d, x, t, y, modulus);
+                if (x < 0)
+                        x += modulus;
+                return x;               
+                                                              
+        }
+
+        inline Element &axpy (Element &r, 
+                              const Element &a, 
+                              const Element &x, 
+                              const Element &y) const {
+                int q;
+                        
+                q  = (int) (((((double) a) * ((double) x)) + (double)y) * modulusinv);  // q could be off by (+/-) 1
+                r = (int) (a * x + y - q*modulus);
+                        
+                        
+                if (r >= modulus)
+                        r -= modulus;
+                else if (r < 0)
+                        r += modulus;
+
+                return r;
+
+        }
+
+        inline Element &addin (Element &x, const Element &y) const {
+                x += y;
+                if (  x >= modulus ) x -= modulus;
+                return x;
+        }
+ 
+        inline Element &subin (Element &x, const Element &y) const {
+                x -= y;
+                if (x < 0) x += modulus;
+                return x;
+        }
+ 
+        inline Element &mulin (Element &x, const Element &y) const {
+                return mul(x,x,y);
+        }
+ 
+        inline Element &divin (Element &x, const Element &y) const {
+                return div(x,x,y);
+        }
+ 
+        inline Element &negin (Element &x) const {
+                if (x == 0) return x; 
+                else return x = modulus - x; 
+        }
+ 
+        inline Element &invin (Element &x) const {
+                return inv (x, x);
+        }
+
+        inline Element &axpyin (Element &r, const Element &a, const Element &x) const {
+                int q;
+                        
+                q  = (int) (((((double) a) * ((double) x)) + (double) r) * modulusinv);  // q could be off by (+/-) 1
+                r = (int) (a * x + r - q*modulus);
+                        
+                        
+                if (r >= modulus)
+                        r -= modulus;
+                else if (r < 0)
+                        r += modulus;
+
+                return r;
+        }
+
+private:
+
+        static void XGCD(int& d, int& s, int& t, int a, int b) {
+                int  u, v, u0, v0, u1, v1, u2, v2, q, r;
+                        
+                int aneg = 0, bneg = 0;
+                        
+                if (a < 0) {
+                        a = -a;
+                        aneg = 1;
                 }
-                
-                std::istream &read (std::istream &is) {
-                        is >> modulus; 
-                        modulusinv = 1 /((double) modulus );
-                        _two64 = (int) ((uint64) (-1) % (uint64) modulus);
-                        _two64 += 1;
-                        if (_two64 >= modulus) _two64 -= modulus;
-                        
-                        return is;
+                        
+                if (b < 0) {
+                        b = -b;
+                        bneg = 1;
                 }
-                
-                std::ostream &write (std::ostream &os, const Element &x) const {
-                        return os << x;
+                        
+                u1 = 1; v1 = 0;
+                u2 = 0; v2 = 1;
+                u = a; v = b;
+                        
+                while (v != 0) {
+                        q = u / v;
+                        r = u % v;
+                        u = v;
+                        v = r;
+                        u0 = u2;
+                        v0 = v2;
+                        u2 =  u1 - q*u2;
+                        v2 = v1- q*v2;
+                        u1 = u0;
+                        v1 = v0;
                 }
-
-                std::istream &read (std::istream &is, Element &x) const {
-                        FieldInt tmp;
-                        is >> tmp;
-                        init(x,tmp); 
-                        return is;
-                }
-                
-
-                template<class Element1>
-                Element &init (Element & x, const Element1 &y) const {
-                        x = y % modulus;
-                        if (x < 0) x += modulus;
-                        return x;
-                }
-
-                Element &init (Element &x, const double &y) const  { 
-                  double z = fmod(y, (double)modulus);
-                  if (z < 0) z += (double)modulus;
-                  z += 0.5;
-                  return x = static_cast<long>(z); //rounds towards 0
-                }
-                Element &init (Element &x, const float  &y) const  { 
-                        float z = fmod(y, (float)modulus);
-                  if (z < 0) z += (float)modulus;
-                  z += 0.5;
-                  return x = static_cast<long>(z); //rounds towards 0