Changeset 0b937c in git

Timestamp:

Jan 31, 2003, 10:09:46 AM (21 years ago)

Author:

Hans Schönemann <hannes@…>

Branches:

(u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')

Children:

6a874ee99f6ab05d5981f0dc589bc510d18d6f6c

Parents:

cd22d1a27db3332385586c32609b3ccce07f4c94

Message:

*hannes: NV_OPS, NV_MAX_PRIME


git-svn-id: file:///usr/local/Singular/svn/trunk@6454 2c84dea3-7e68-4137-9b89-c4e89433aadc

Location:

Singular

Files:

• modulop.cc (modified) (13 diffs)
• modulop.h (modified) (5 diffs)

Singular/modulop.cc

@@ -2 +2 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: modulop.cc,v 1.26 2001-08-27 14:47:12 Singular Exp $ */
+/* $Id: modulop.cc,v 1.27 2003-01-31 09:09:46 Singular Exp $ */
 /*
 * ABSTRACT: numbers modulo p (<=32003)
@@ -24 +24 @@
 int npMapPrime;
 
+#ifdef HAVE_DIV_MOD
+CARDINAL *npInvTable=NULL;
+#endif
+
+#if !defined(HAVE_DIV_MOD) || !defined(HAVE_MULT_MOD)
 CARDINAL *npExpTable=NULL;
 CARDINAL *npLogTable=NULL;
+#endif
 
 
@@ -34 +40 @@
 }
 
-unsigned long npMultMod(unsigned long a, unsigned long b)
-{
-  unsigned long c = a*b;
-  c = c % npPrimeM;
-  assume(c == (unsigned long) npMultM((number) a, (number) b));
-  return c;
-}
+//unsigned long npMultMod(unsigned long a, unsigned long b)
+//{
+//  unsigned long c = a*b;
+//  c = c % npPrimeM;
+//  assume(c == (unsigned long) npMultM((number) a, (number) b));
+//  return c;
+//}
 
 number npMult (number a,number b)
@@ -47 +53 @@
     return (number)0;
   else
-  {
     return npMultM(a,b);
-  }
 }
 
@@ -96 +100 @@
 }
 
+#ifdef HAVE_DIV_MOD
+#if 1 //ifdef HAVE_NTL // in ntl.a
+extern void XGCD(long& d, long& s, long& t, long a, long b);
+#else
+void XGCD(long& d, long& s, long& t, long a, long b)
+{
+   long  u, v, u0, v0, u1, v1, u2, v2, q, r;
+
+   long aneg = 0, bneg = 0;
+
+   if (a < 0) {
+      a = -a;
+      aneg = 1;
+   }
+
+   if (b < 0) {
+      b = -b;
+      bneg = 1;
+   }
+
+   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;
+   }
+
+   if (aneg)
+      u1 = -u1;
+
+   if (bneg)
+      v1 = -v1;
+
+   d = u;
+   s = u1;
+   t = v1;
+}
+#endif
+
+long InvMod(long a)
+{
+   long d, s, t;
+
+   XGCD(d, s, t, a, npPrimeM);
+   assume (d == 1);
+   if (s < 0)
+      return s + npPrimeM;
+   else
+      return s;
+}
+#endif
+
+inline number npInversM (number c)
+{
+#ifndef HAVE_DIV_MOD
+  return (number)npExpTable[npPminus1M - npLogTable[(int)c]];
+#else
+  CARDINAL inv=npInvTable[(int)c];
+  if (inv==0)
+  {
+    inv=InvMod((long)c);
+    npInvTable[(int)c]=inv;
+  }
+  return (number)inv;
+#endif
+}
+
 number npDiv (number a,number b)
 {
   if ((int)a==0)
     return (number)0;
+#ifndef HAVE_DIV_MOD
   else if ((int)b==0)
   {
@@ -112 +194 @@
     return (number)npExpTable[s];
   }
-}
-
+#else
+  number inv=npInversM(b);
+  return npMultM(a,inv);
+#endif
+}
 number  npInvers (number c)
 {
@@ -121 +206 @@
     return (number)0;
   }
-  return (number)npExpTable[npPminus1M - npLogTable[(int)c]];
+  return npInversM(c);
 }
 
@@ -195 +280 @@
     s = npEati(s, &n);
   }
-  *a = npDiv((number)z,(number)n);
+  if (n == 1)
+    *a = (number)z;
+  else 
+#ifdef NV_OPS
+    if (npPrimeM>NV_MAX_PRIME)
+      *a = nvDiv((number)z,(number)n);
+    else
+#endif    
+      *a = npDiv((number)z,(number)n);
   return s;
 }
@@ -216 +309 @@
 
 //  if (c==npPrimeM) return;
-//  if (npPrimeM > 1)
-//  {
-//    omFreeSize( (ADDRESS)npExpTable,npPrimeM*sizeof(CARDINAL) );
-//    omFreeSize( (ADDRESS)npLogTable,npPrimeM*sizeof(CARDINAL) );
-//  }
   if ((c>1) || (c<(-1)))
   {
@@ -226 +314 @@
     else     npPrimeM = -c;
     npPminus1M = npPrimeM - 1;
-//    npExpTable= (CARDINAL *)omAlloc( npPrimeM*sizeof(CARDINAL) );
-//     npLogTable= (CARDINAL *)omAlloc( npPrimeM*sizeof(CARDINAL) );
-//     omMarkAsStaticAddr(npExpTable);
-//     omMarkAsStaticAddr(npLogTable);
-//     memcpy(npExpTable,r->cf->npExpTable,npPrimeM*sizeof(CARDINAL));
-//     memcpy(npLogTable,r->cf->npLogTable,npPrimeM*sizeof(CARDINAL));
+#ifdef HAVE_DIV_MOD
+    npInvTable=r->cf->npInvTable;
+#endif
+#if !defined(HAVE_DIV_MOD) || !defined(HAVE_MULT_MOD)
     npExpTable=r->cf->npExpTable;
     npLogTable=r->cf->npLogTable;
     npGen = npExpTable[1];
+#endif
   }
   else
   {
     npPrimeM=0;
+#ifdef HAVE_DIV_MOD
+    npInvTable=NULL;
+#endif
+#if !defined(HAVE_DIV_MOD) || !defined(HAVE_MULT_MOD)
     npExpTable=NULL;
     npLogTable=NULL;
+#endif
   }
 }
@@ -253 +345 @@
     else     r->cf->npPrimeM = -c;
     r->cf->npPminus1M = r->cf->npPrimeM - 1;
-    r->cf->npExpTable= (CARDINAL *)omAlloc( r->cf->npPrimeM*sizeof(CARDINAL) );
-    r->cf->npLogTable= (CARDINAL *)omAlloc( r->cf->npPrimeM*sizeof(CARDINAL) );
-    r->cf->npExpTable[0] = 1;
-    r->cf->npLogTable[0] = 0;
-    if (r->cf->npPrimeM > 2)
+#ifdef NV_OPS
+    if (r->cf->npPrimeM <=NV_MAX_PRIME)
+#endif
     {
-      w = 1;
-      loop
+#if !defined(HAVE_DIV_MOD) || !defined(HAVE_MULT_MOD)
+      r->cf->npExpTable=(CARDINAL *)omAlloc( r->cf->npPrimeM*sizeof(CARDINAL) );
+      r->cf->npLogTable=(CARDINAL *)omAlloc( r->cf->npPrimeM*sizeof(CARDINAL) );
+      r->cf->npExpTable[0] = 1;
+      r->cf->npLogTable[0] = 0;
+      if (r->cf->npPrimeM > 2)
       {
-        r->cf->npLogTable[1] = 0;
-        w++;
-        i = 0;
+        w = 1;
         loop
         {
-          i++;
-          r->cf->npExpTable[i] = (int)(((long)w * (long)r->cf->npExpTable[i-1])
+          r->cf->npLogTable[1] = 0;
+          w++;
+          i = 0;
+          loop
+          {
+            i++;
+            r->cf->npExpTable[i] =(int)(((long)w * (long)r->cf->npExpTable[i-1])
                                  % r->cf->npPrimeM);
-          r->cf->npLogTable[r->cf->npExpTable[i]] = i;
-          if (/*(i == npPrimeM - 1 ) ||*/ (r->cf->npExpTable[i] == 1))
+            r->cf->npLogTable[r->cf->npExpTable[i]] = i;
+            if (/*(i == npPrimeM - 1 ) ||*/ (r->cf->npExpTable[i] == 1))
+              break;
+          }
+          if (i == r->cf->npPrimeM - 1)
             break;
         }
-        if (i == r->cf->npPrimeM - 1)
-          break;
       }
-    }
-    else
-    {
-      r->cf->npExpTable[1] = 1;
-      r->cf->npLogTable[1] = 0;
+      else
+      {
+        r->cf->npExpTable[1] = 1;
+        r->cf->npLogTable[1] = 0;
+      }
+#endif
+#ifdef HAVE_DIV_MOD
+      r->cf->npInvTable=(CARDINAL*)omAlloc0( r->cf->npPrimeM*sizeof(CARDINAL) );
+#endif
     }
   }
@@ -348 +450 @@
 #if defined(LDEBUG)
   res->debug=123456;
+#endif
+#ifndef NL_OLDCOPY
+  res->ref=1;
 #endif
   dest = &(res->z);
@@ -412 +517 @@
   return NULL;      /* default */
 }
+
+// -----------------------------------------------------------
+//  operation for very large primes (32003< p < 2^31-1)
+// ----------------------------------------------------------
+#ifdef NV_OPS
+
+number nvMult (number a,number b)
+{
+  if (((int)a == 0) || ((int)b == 0))
+    return (number)0;
+  else
+    return nvMultM(a,b);
+}
+
+long nvInvMod(long a)
+{
+   long  s, t;
+
+   long  u, v, u0, v0, u1, v1, u2, v2, q, r;
+
+   u1=1; v1=0;
+   u2=0; v2=1;
+   u = a; v = npPrimeM;
+
+   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;
+   }
+
+   s = u1;
+   //t = v1;
+   if (s < 0)
+      return s + npPrimeM;
+   else
+      return s;
+}
+
+inline number nvInversM (number c)
+{
+  long inv=InvMod((long)c);
+  return (number)inv;
+}
+
+number nvDiv (number a,number b)
+{
+  if ((int)a==0)
+    return (number)0;
+  else if ((int)b==0)
+  {
+    WerrorS("div by 0");
+    return (number)0;
+  }
+  else
+  {
+    number inv=nvInversM(b);
+    return nvMultM(a,inv);
+  }
+}
+number  nvInvers (number c)
+{
+  if ((int)c==0)
+  {
+    WerrorS("1/0");
+    return (number)0;
+  }
+  return nvInversM(c);
+}
+#endif

Singular/modulop.h

@@ -4 +4 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: modulop.h,v 1.16 2001-10-09 16:36:09 Singular Exp $ */
+/* $Id: modulop.h,v 1.17 2003-01-31 09:09:46 Singular Exp $ */
 /*
 * ABSTRACT: numbers modulo p (<=32003)
@@ -10 +10 @@
 #include "structs.h"
 
+// defines are in struct.h
 // define if a*b is with mod instead of tables
-#define HAVE_MULT_MOD
+//#define HAVE_MULT_MOD 
+// define if a/b is with mod instead of tables
+//#define HAVE_DIV_MOD 
+// define if an if should be used
+//#define HAVE_GENERIC_ADD
+
+// enable large primes (32003 < p < 2^31-)
+#define NV_OPS
+#define NV_MAX_PRIME 32003
 
 extern int npPrimeM;
@@ -46 +55 @@
 /*-------specials for spolys, do NOT use otherwise--------------------------*/
 /* for npMultM, npSubM, npNegM, npEqualM : */
+#ifdef HAVE_DIV_MOD
+extern CARDINAL *npInvTable;
+#else
+#ifndef HAVE_MULT_MOD
 extern int npPminus1M;
 extern CARDINAL *npExpTable;
 extern CARDINAL *npLogTable;
+#endif
+#endif
 
+#if 0
+inline number npMultM(number a, number b)
+// return (a*b)%n
+{
+   double ab;
+   long q, res;
+
+   ab = ((double) ((int)a)) * ((double) ((int)b));
+   q  = (long) (ab/((double) npPrimeM));  // q could be off by (+/-) 1
+   res = (long) (ab - ((double) q)*((double) npPrimeM));
+   res += (res >> 31) & npPrimeM;
+   res -= npPrimeM;
+   res += (res >> 31) & npPrimeM;
+   return (number)res;
+}
+#endif
 #ifdef HAVE_MULT_MOD
-inline number npMultM(number a, number b)
+static inline number npMultM(number a, number b)
 {
-  return (number)
+  return (number) 
     ((((unsigned long) a)*((unsigned long) b)) % ((unsigned long) npPrimeM));
 }
 #else
-inline number npMultM(number a, number b)
+static inline number npMultM(number a, number b)
 {
   int x = npLogTable[(int)a]+npLogTable[(int)b];
@@ -84 +115 @@
 }
 #endif
-
-inline number npAddM(number a, number b)
+#ifdef HAVE_GENERIC_ADD
+static inline number npAddM(number a, number b)
 {
   int r = (int)a + (int)b;
   return (number)(r >= npPrimeM ? r - npPrimeM : r);
 }
+static inline number npSubM(number a, number b)
+{
+  return (number)((int)a<(int)b ?
+                       npPrimeM-(int)b+(int)a : (int)a-(int)b);
+}
+#else
+static inline number npAddM(number a, number b)
+{
+   int res = (int)a + (int)b;
+   res -= npPrimeM;
+   res += (res >> 31) & npPrimeM;
+   return (number)res;
+}
+static inline number npSubM(number a, number b)
+{
+   int res = (int)a - (int)b;
+   res += (res >> 31) & npPrimeM;
+   return (number)res;
+}
+#endif
 
-inline BOOLEAN npIsZeroM (number  a)
+static inline BOOLEAN npIsZeroM (number  a)
 {
   return 0 == (int)a;
@@ -103 +154 @@
 */
 
-#define npSubM(a,b)    (number)((int)a<(int)b ?\
-                       npPrimeM-(int)b+(int)a : (int)a-(int)b)
-
 #define npNegM(A)      (number)(npPrimeM-(int)(A))
 #define npEqualM(A,B)  ((int)A==(int)B)
-#define npIsZeroM(a)   (0 == (int)a)
+
+
+#ifdef NV_OPS
+static inline number nvMultM(number a, number b)
+{
+  return (number) 
+    ((((unsigned long long) a)*((unsigned long long) b)) % ((unsigned long long) npPrimeM));
+}
+number  nvMult        (number a, number b);
+number  nvDiv         (number a, number b);
+number  nvInvers      (number c);
 #endif
 
+#endif