Changeset c9b1d66 in git

Timestamp:

Dec 16, 2010, 2:22:47 PM (12 years ago)

Author:

Martin Lee <martinlee84@…>

Branches:

(u'spielwiese', '91fdef05f09f54b8d58d92a472e9c4a43aa4656f')

Children:

8af49d09040b6ced3fa697fae0b4338a908944a5

Parents:

d08ed84a213b4110fe1dbe0d38276dce7695a481

Message:

added division with remainder via Newton inversion


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

Location:

factory

Files:

• facHensel.cc (modified) (1 diff)
• facHensel.h (modified) (2 diffs)

factory/facHensel.cc

@@ -395 +395 @@
     buf2= prodMod (tmp2, M);
     return mulMod (buf1, buf2, M);
+  }
+}
+
+
+CanonicalForm reverse (const CanonicalForm& F, int d)
+{
+  if (d == 0)
+    return F;
+  CanonicalForm A= F;
+  Variable y= F.mvar();
+  Variable x= Variable (1);
+  A= swapvar (A, x, y);
+  CanonicalForm result= 0;
+  CFIterator i= A;
+  while (d - i.exp() < 0)
+    i++;
+
+  for (; i.hasTerms() && (d - i.exp() >= 0); i++)
+    result += swapvar (i.coeff(),x,y)*power (x, d - i.exp());
+  return result;
+}
+
+CanonicalForm
+newtonInverse (const CanonicalForm& F, const int n, const CanonicalForm& M)
+{
+  int l= ilog2(n);
+
+  CanonicalForm g= mod (F, M)[0] [0];
+
+  ASSERT (!g.isZero(), "expected a unit");
+
+  Variable alpha;
+
+  if (!g.isOne())
+    g = 1/g;
+  Variable x= Variable (1);
+  CanonicalForm result;
+  int exp= 0;
+  if (n & 1)
+  {
+    result= g;
+    exp= 1;
+  }
+  CanonicalForm h;
+
+  for (int i= 1; i <= l; i++)
+  {
+    h= mulMod2 (g, mod (F, power (x, (1 << i))), M);
+    h= mod (h, power (x, (1 << i)) - 1);
+    h= div (h, power (x, (1 << (i - 1))));
+    h= mod (h, M);
+    g -= power (x, (1 << (i - 1)))*
+         mod (mulMod2 (g, h, M), power (x, (1 << (i - 1))));
+
+    if (n & (1 << i))
+    {
+      if (exp)
+      {
+        h= mulMod2 (result, mod (F, power (x, exp + (1 << i))), M);
+        h= mod (h, power (x, exp + (1 << i)) - 1);
+        h= div (h, power (x, exp));
+        h= mod (h, M);
+        result -= power(x, exp)*mod (mulMod2 (g, h, M),
+                                       power (x, (1 << i)));
+        exp += (1 << i);
+      }
+      else
+      {
+        exp= (1 << i);
+        result= g;
+      }
+    }
+  }
+
+  return result;
+}
+
+CanonicalForm
+newtonDiv (const CanonicalForm& F, const CanonicalForm& G, const CanonicalForm&
+           M)
+{
+  ASSERT (getCharacteristic() > 0, "positive characteristic expected");
+  ASSERT (CFFactory::gettype() != GaloisFieldDomain, "no GF expected");
+
+  CanonicalForm A= mod (F, M);
+  CanonicalForm B= mod (G, M);
+
+  Variable x= Variable (1);
+  int degA= degree (A, x);
+  int degB= degree (B, x);
+  int m= degA - degB;
+  if (m < 0)
+    return 0;
+
+  CanonicalForm Q;
+  if (degB <= 1)
+  {
+    CanonicalForm R;
+    divrem2 (A, B, Q, R, M);
+  }
+  else
+  {
+    CanonicalForm R= reverse (A, degA);
+    CanonicalForm revB= reverse (B, degB);
+    revB= newtonInverse (revB, m + 1, M);
+    Q= mulMod2 (R, revB, M);
+    Q= mod (Q, power (x, m + 1));
+    Q= reverse (Q, m);
+  }
+
+  return Q;
+}
+
+void
+newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
+              CanonicalForm& R, const CanonicalForm& M)
+{
+  CanonicalForm A= mod (F, M);
+  CanonicalForm B= mod (G, M);
+  Variable x= Variable (1);
+  int degA= degree (A, x);
+  int degB= degree (B, x);
+  int m= degA - degB;
+
+  if (m < 0)
+  {
+    R= A;
+    Q= 0;
+    return;
+  }
+
+  if (degB <= 1)
+  {
+     divrem2 (A, B, Q, R, M);
+  }
+  else
+  {
+    R= reverse (A, degA);
+
+    CanonicalForm revB= reverse (B, degB);
+    revB= newtonInverse (revB, m + 1, M);
+    Q= mulMod2 (R, revB, M);
+
+    Q= mod (Q, power (x, m + 1));
+    Q= reverse (Q, m);
+
+    R= A - mulMod2 (Q, B, M);
   }
 }

factory/facHensel.h

@@ -29 +29 @@
 #include "algext.h"
 
-CanonicalForm mul (const CanonicalForm& F, const CanonicalForm& G);
-CanonicalForm mulMod3 (const CanonicalForm& F, const CanonicalForm& G, const CFList& MOD);
 /// multiplication of univariate polys over a finite field using NTL, if we are
 /// in GF factory's default multiplication is used.
@@ -139 +137 @@
                                    ///< Variables of level higher than 1
             );
+
+
+/// division with remainder of @a F by
+/// @a G wrt Variable (1) modulo @a M using Newton inversion
+///
+/// @return @a Q returns the dividend, @a R returns the remainder.
+/// @sa divrem2(), newtonDiv()
+void
+newtonDivrem (const CanonicalForm& F, ///< [in] bivariate, compressed polynomial
+              const CanonicalForm& G, ///< [in] bivariate, compressed polynomial
+                                      ///< which is monic in Variable (1)
+              CanonicalForm& Q,       ///< [in,out] dividend
+              CanonicalForm& R,       ///< [in,out] remainder, degree (R, 1) <
+                                      ///< degree (G, 1)
+              const CanonicalForm& M  ///< [in] power of Variable (2)
+             );
+
+/// division of @a F by
+/// @a G wrt Variable (1) modulo @a M using Newton inversion
+///
+/// @return @a newtonDiv returns the dividend
+/// @sa divrem2(), newtonDivrem()
+CanonicalForm
+newtonDiv (const CanonicalForm& F, ///< [in] bivariate, compressed polynomial
+           const CanonicalForm& G, ///< [in] bivariate, compressed polynomial
+                                   ///< which is monic in Variable (1)
+           const CanonicalForm& M  ///< [in] power of Variable (2)
+          );
 
 /// reduce @a F modulo elements in @a M.