Peter Bruin on Wed, 07 Jun 2017 19:43:54 +0200


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Speed up RgX_mul, RgX_sqr, RgX_divrem over Z, Fp, FpXQ


Hello,

The following resultant computation (with output of degree 9480)
currently takes about 10 minutes:

g = x^11+t^13+3*x^7+(2*(t+1)^8+2*t)^7*x^6+(2*t^8+t^6+2*t)^45*x^3+(t+1)^5*t^4*(2*x+x)+t^12+2;
f = Mod(g^2+x^4*(t+2)^4-t^12, 5);
r = polresultant(f, f', 'x);

Claus Fieker pointed out that this only takes about 0.3 seconds in Bill
Hart's implementation of the subresultant algorithm in Julia.

It turns out that the PARI implementation can also be significantly sped
up by dispatching to more specialised functions in RgX_mul, RgX_sqr and
RgX_divrem.  After applying the two attached patches, the time for the
above computation drops to about 0.7 seconds.

A possible further speed-up over small prime finite fields would be to
implement a function Flx_to_mod_inplace that combines the calls to
Flx_to_ZX_inplace and FpX_to_mod.

Thanks,

Peter


>From c83deaa76e23b9951745f17766538b695c4366aa Mon Sep 17 00:00:00 2001
From: Peter Bruin <P.J.Bruin@math.leidenuniv.nl>
Date: Wed, 7 Jun 2017 02:27:22 +0200
Subject: [PATCH 1/2] check for ZX, FpX, FpXQ in RgX_mul and RgX_sqr

---
 src/basemath/RgX.c  | 55 +++++++++++++++++++++++++++++++++++++++++++++++++++--
 src/basemath/gen1.c | 17 +----------------
 src/test/32/ell     |  4 ++--
 3 files changed, 56 insertions(+), 20 deletions(-)

diff --git a/src/basemath/RgX.c b/src/basemath/RgX.c
index 3c02342..4c16aee 100644
--- a/src/basemath/RgX.c
+++ b/src/basemath/RgX.c
@@ -1514,14 +1514,65 @@ RgX_mul_normalized(GEN A, long a, GEN B, long b)
 GEN
 RgX_mul(GEN x, GEN y)
 {
-  GEN z = RgX_mulspec(y+2, x+2, lgpol(y), lgpol(x));
+  GEN p, pol, xx, z;
+  pari_sp av;
+  if (RgX_is_ZX(x) && RgX_is_ZX(y)) return ZX_mul(x, y);
+  av = avma;
+  p = NULL;
+  if (RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p))
+  {
+    if (lgefint(p) == 3)
+    {
+      ulong pp = uel(p, 2);
+      z = Flx_to_ZX_inplace(Flx_mul(RgX_to_Flx(x, pp),
+                                    RgX_to_Flx(y, pp), pp));
+    }
+    else
+      z = FpX_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
+    return gerepileupto(av, FpX_to_mod(z, p));
+  }
+  p = NULL; pol = NULL; xx = x;
+  if (ff_poltype(&xx, &p, &pol) && ff_poltype(&y, &p, &pol))
+  {
+    z = ZX_mul(xx, y);
+    if (p) z = FpX_to_mod(z, p);
+    if (pol) z = Kronecker_to_mod(z, pol);
+    return gerepileupto(av, z);
+  }
+  avma = av;
+  z = RgX_mulspec(y+2, x+2, lgpol(y), lgpol(x));
   setvarn(z,varn(x)); return z;
 }
 
 GEN
 RgX_sqr(GEN x)
 {
-  GEN z = RgX_sqrspec(x+2, lgpol(x));
+  GEN p, pol, z;
+  pari_sp av;
+  if (RgX_is_ZX(x)) return ZX_sqr(x);
+  av = avma;
+  p = NULL;
+  if (RgX_is_FpX(x, &p))
+  {
+    if (lgefint(p) == 3)
+    {
+      ulong pp = uel(p, 2);
+      z = Flx_to_ZX_inplace(Flx_sqr(RgX_to_Flx(x, pp), pp));
+    }
+    else
+      z = FpX_sqr(RgX_to_FpX(x, p), p);
+    return gerepileupto(av, FpX_to_mod(z, p));
+  }
+  p = NULL; pol = NULL;
+  if (ff_poltype(&x, &p, &pol))
+  {
+    z = ZX_sqr(x);
+    if (p) z = FpX_to_mod(z, p);
+    if (pol) z = Kronecker_to_mod(z, pol);
+    return gerepileupto(av, z);
+  }
+  avma = av;
+  z = RgX_sqrspec(x+2, lgpol(x));
   setvarn(z,varn(x)); return z;
 }
 
diff --git a/src/basemath/gen1.c b/src/basemath/gen1.c
index c2223a9..ff288b0 100644
--- a/src/basemath/gen1.c
+++ b/src/basemath/gen1.c
@@ -1880,7 +1880,6 @@ gmul(GEN x, GEN y)
         if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
         else                     return RgX_Rg_mul(y, x);
       }
-      if (RgX_is_ZX(x) && RgX_is_ZX(y)) return ZX_mul(x,y);
       return RgX_mul(x, y);
 
     case t_SER: {
@@ -2245,21 +2244,7 @@ gsqr(GEN x)
 
     case t_FFELT: return FF_sqr(x);
 
-    case t_POL:
-    {
-      GEN a = x, p = NULL, pol = NULL;
-      av = avma;
-      if (RgX_is_ZX(x)) return ZX_sqr(x);
-      if (ff_poltype(&x,&p,&pol))
-      {
-        z = ZX_sqr(x);
-        if (p) z = FpX_to_mod(z,p);
-        if (pol) z = Kronecker_to_mod(z,pol);
-        z = gerepileupto(av, z);
-      }
-      else { avma = av; z = RgX_sqr(a); }
-      return z;
-    }
+    case t_POL: return RgX_sqr(x);
 
     case t_SER:
       lx = lg(x);
diff --git a/src/test/32/ell b/src/test/32/ell
index 8954671..26c5492 100644
--- a/src/test/32/ell
+++ b/src/test/32/ell
@@ -140,8 +140,8 @@ x^6 + 884*z^4*x^4 + (50*z^4 + 809*z^3 + 859*z^2 + 150*z + 50)*x^2 + (5*z^4 +
 *z, z^2 + 5)*x^6 + Mod(-150, z^2 + 5)*x^4 + Mod(60*z, z^2 + 5)*x^2 + Mod(25,
  z^2 + 5)]
 (E)->ellxn(E,3):Mod(4, 1009):[Mod(1, 1009)*x^9 + Mod(961, 1009)*x^7 + Mod(48
-0, 1009)*x^5 + Mod(286, 1009)*x^3 + Mod(286, 1009)*x, 9*x^8 + Mod(144, 1009)
-*x^6 + Mod(480, 1009)*x^4 + Mod(241, 1009)*x^2 + Mod(256, 1009)]
+0, 1009)*x^5 + Mod(286, 1009)*x^3 + Mod(286, 1009)*x, Mod(9, 1009)*x^8 + Mod
+(144, 1009)*x^6 + Mod(480, 1009)*x^4 + Mod(241, 1009)*x^2 + Mod(256, 1009)]
 (E)->ellxn(E,3):z:[Mod(1, 1009)*x^9 + 997*z*x^7 + 30*z^2*x^5 + 36*z^3*x^3 + 
 9*z^4*x, 9*x^8 + 36*z*x^6 + 30*z^2*x^4 + 997*z^3*x^2 + z^4]
 (E)->ellmul(E,[0,0],0):1:[0]
-- 
2.10.2

>From c33013f45b88da612d17ffbb56d6a58d41695a07 Mon Sep 17 00:00:00 2001
From: Peter Bruin <P.J.Bruin@math.leidenuniv.nl>
Date: Wed, 7 Jun 2017 02:32:51 +0200
Subject: [PATCH 2/2] check for FpX in RgX_divrem

---
 src/basemath/RgX.c | 33 ++++++++++++++++++++++++++++-----
 1 file changed, 28 insertions(+), 5 deletions(-)

diff --git a/src/basemath/RgX.c b/src/basemath/RgX.c
index 4c16aee..b989fee 100644
--- a/src/basemath/RgX.c
+++ b/src/basemath/RgX.c
@@ -1636,17 +1636,16 @@ RgX_div_by_X_x(GEN a, GEN x, GEN *r)
   return z;
 }
 /* Polynomial division x / y:
- *   if z = ONLY_REM  return remainder, otherwise return quotient
- *   if z != NULL set *z to remainder
- *   *z is the last object on stack (and thus can be disposed of with cgiv
- *   instead of gerepile) */
+ *   if pr = ONLY_REM return remainder, otherwise return quotient
+ *   if pr = ONLY_DIVIDES return quotient if division is exact, else NULL
+ *   if pr != NULL set *pr to remainder, as the last object on stack */
 /* assume, typ(x) = typ(y) = t_POL, same variable */
 GEN
 RgX_divrem(GEN x, GEN y, GEN *pr)
 {
   pari_sp avy, av, av1;
   long dx,dy,dz,i,j,sx,lr;
-  GEN z,p1,p2,rem,y_lead,mod;
+  GEN z,p1,p2,rem,y_lead,mod,p;
   GEN (*f)(GEN,GEN);
 
   if (!signe(y)) pari_err_INV("RgX_divrem",y);
@@ -1682,6 +1681,30 @@ RgX_divrem(GEN x, GEN y, GEN *pr)
 
   /* x,y in R[X], y non constant */
   av = avma;
+  p = NULL;
+  if (RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p) && p)
+  {
+    if (lgefint(p) == 3)
+    {
+      ulong pp = uel(p, 2);
+      z = Flx_divrem(RgX_to_Flx(x, pp), RgX_to_Flx(y, pp), pp, pr);
+      if (!z) { avma = av; return NULL; }
+      z = Flx_to_ZX_inplace(z);
+      if (pr && pr != ONLY_REM && pr != ONLY_DIVIDES)
+	*pr = Flx_to_ZX_inplace(*pr);
+    }
+    else
+    {
+      z = FpX_divrem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, pr);
+      if (!z) { avma = av; return NULL; }
+    }
+    z = FpX_to_mod(z, p);
+    if (!pr || pr == ONLY_REM || pr == ONLY_DIVIDES)
+      return gerepileupto(av, z);
+    *pr = FpX_to_mod(*pr, p);
+    gerepileall(av, 2, pr, &z);
+    return z;
+  }
   switch(typ(y_lead))
   {
     case t_REAL:
-- 
2.10.2