Changeset 79020f in git

Timestamp:

Jan 6, 2012, 9:28:16 PM (12 years ago)

Author:

Oleksandr Motsak <motsak@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

40e88db63c509564c2b00b7142b39ef6b846f72e

Parents:

1cbb1f452324cf72b33ef5a1cb6256d97bcbcbd0

git-author:

Oleksandr Motsak <motsak@mathematik.uni-kl.de>2012-01-06 21:28:16+01:00

git-committer:

Oleksandr Motsak <motsak@mathematik.uni-kl.de>2012-01-06 21:59:37+01:00

Message:

fixed arithmetic overflow by rewriting nlModP

CHG: nlModP (longrat) was rewritten in order to simplify and generalize it
CHG/FIX: nlModP _works_ again as in the old Singular (uses n_Div)

Location:

libpolys

Files:

• coeffs/longrat.cc (modified) (1 diff)
• coeffs/longrat.h (modified) (1 diff)
• coeffs/modulop.cc (modified) (3 diffs)
• coeffs/modulop.h (modified) (1 diff)
• polys/ext_fields/algext.cc (modified) (1 diff)
• polys/ext_fields/transext.cc (modified) (2 diffs)

libpolys/coeffs/longrat.cc

@@ -1256 +1256 @@
 }
 
-int nlModP(number n, int p, const coeffs r)
-{
-  if (SR_HDL(n) & SR_INT)
-  {
-    long i=SR_TO_INT(n);
-    if (i<0L) return (((long)p)-((-i)%((long)p)));
-    return i%((long)p);
-  }
-  int iz=(int)mpz_fdiv_ui(n->z,(unsigned long)p);
-  if (n->s!=3)
-  {
-    int in=mpz_fdiv_ui(n->n,(unsigned long)p);
-    long  s;
-
-    long  u, v, u0, v0, u1, v1, u2, v2, q, r;
-
-     u1=1; v1=0;
-     u2=0; v2=1;
-     u = in; v = p;
-
-     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;
-     if (s < 0) s+=p;
-     u=(s*((long)iz)) % ((long)p); // BUG: integer overflow!!!
-     return (int)u;
-  }
-  return iz;
+// Map q \in QQ \to Zp
+// src = Q, dst = Zp (or an extension of Zp?)
+number nlModP(number q, const coeffs Q, const coeffs Zp)
+{
+  assume( getCoeffType(Q) == ID );
+
+  const int p = n_GetChar(Zp);
+  assume( p > 0 );
+
+  const long P = p;
+  assume( P > 0 );
+
+  // embedded long within q => only long numerator has to be converted
+  // to int (modulo char.)
+  if (SR_HDL(q) & SR_INT)
+  {
+    long i = SR_TO_INT(q);
+    if (i<0L)
+      return n_Init( static_cast<int>( P - ((-i)%P) ), Zp);
+    else
+      return n_Init( static_cast<int>( i % P ), Zp );
+  }
+
+  const unsigned long PP = p;
+
+  // numerator modulo char. should fit into int
+  number z = n_Init( static_cast<int>(mpz_fdiv_ui(q->z, PP)), Zp );  
+
+  // denominator != 1?
+  if (q->s!=3)
+  {
+    // denominator modulo char. should fit into int
+    number n = n_Init( static_cast<int>(mpz_fdiv_ui(q->n, PP)), Zp );
+
+    number res = n_Div( z, n, Zp );
+
+    n_Delete(&z, Zp);
+    n_Delete(&n, Zp);
+
+    return res;
+  }
+
+  return z;
 }
 

libpolys/coeffs/longrat.h

@@ -87 +87 @@
 const char *   nlRead (const char *s, number *a, const coeffs r);
 void     nlWrite(number &a, const coeffs r);
-int      nlModP(number n, int p, const coeffs r);
+
+/// Map q \in QQ \to Zp
+number nlModP(number q, const coeffs Q, const coeffs Zp);
+
 int      nlSize(number n, const coeffs r);
 number   nlGetDenom(number &n, const coeffs r);

libpolys/coeffs/modulop.cc

@@ -460 +460 @@
   //r->cfName = ndName;
   r->cfInpMult=ndInpMult;
-  r->cfInit_bigint= npMap0;
+  r->cfInit_bigint= nlModP; // npMap0;
 #ifdef NV_OPS
   if (c>NV_MAX_PRIME)
@@ -545 +545 @@
 #endif
 
-number npMap0(number from, const coeffs src, const coeffs dst_r)
-{
-  int nlModP(number n, int p, const coeffs r);
-
-  return npInit(nlModP(from, dst_r->npPrimeM, src),dst_r);
-}
-
 number npMapP(number from, const coeffs src, const coeffs dst_r)
 {
@@ -686 +679 @@
   if (nCoeff_is_Q(src))
   {
-    return npMap0;
+    return nlModP; // npMap0;
   }
   if ( nCoeff_is_Zp(src) )

libpolys/coeffs/modulop.h

@@ -55 +55 @@
 nMapFunc npSetMap(const coeffs src, const coeffs dst);
 number  npMapP(number from, const coeffs src, const coeffs r);
-number  npMap0(number from, const coeffs src, const coeffs r);
 /*-------specials for spolys, do NOT use otherwise--------------------------*/
 /* for npMultM, npSubM, npNegM, npEqualM : */

libpolys/polys/ext_fields/algext.cc

@@ -650 +650 @@
   if (n_IsZero(a, src)) return NULL;
   int p = rChar(dst->extRing);
-  int n = nlModP(a, p, src);
-  number q = n_Init(n, dst->extRing->cf);
-  poly result = p_One(dst->extRing);
-  p_SetCoeff(result, q, dst->extRing);
+
+  number q = nlModP(a, src, dst->extRing->cf);
+
+  poly result = p_NSet(q, dst->extRing);
+  
   return (number)result;
 }

libpolys/polys/ext_fields/transext.cc

@@ -1119 +1119 @@
   if (n_IsZero(a, src)) return NULL;
   int p = rChar(dst->extRing);
-  int n = nlModP(a, p, src);
-  number q = n_Init(n, dst->extRing->cf);
-  poly g;
+  number q = nlModP(a, src, dst->extRing->cf);
+
   if (n_IsZero(q, dst->extRing->cf))
   {
@@ -1127 +1126 @@
     return NULL;
   }
-  g = p_One(dst->extRing);
-  p_SetCoeff(g, q, dst->extRing);
+  
+  poly g = p_NSet(q, dst->extRing);
+
   fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
   NUM(f) = g; // DEN(f) = NULL; COM(f) = 0;