PATCH: polcyclo_eval in roots of unity

[Jeroen Demeyer <jdemeyer@cage.ugent.be>]

Hello pari-dev,

This patch fixes the libpari function polcyclo_eval (the GP function 
polcyclo) such that evaluating in a root of unity now works. 
Additionally, I changed it such that it can compute n-th cyclotomic 
polynomials where n is a t_INT (instead of a long).

The following now all work correctly:

gp> polcyclo(20, I)
%8 = 5

gp> polcyclo(10, Mod(3,11))
%17 = Mod(6, 11)

gp> polcyclo(10, 2 + O(11))
%32 = O(11)

gp> z15 = Mod(t,polcyclo(15,t));
gp> polcyclo(15, z15)
%23 = 0
gp> polcyclo(30, z15)
%24 = Mod(2*t^7 - 2*t^5 - 2*t^3 + 2*t, t^8 - t^7 + t^5 - t^4 + t^3 - t + 1)

gp> polcyclo(30, -1.0)
%5 = 1.0000000000000000000000000000000000000

gp> polcyclo(2^128-1, Mod(x,x^10))
%13 = Mod(-x^7 - x^6 - x^5 + x^2 + x + 1, x^10)

However, there are still problems when working in a ring with 
zero-divisors, for example:

gp> polcyclo(10, Mod(3,4))
  ***   at top-level: polcyclo(10,Mod(3,4)
  ***                 ^--------------------
  *** polcyclo: impossible inverse modulo: Mod(2, 4).


Please test and let me know what you think.

Cheers,
Jeroen.

Index: src/basemath/arith2.c
===================================================================
--- src/basemath/arith2.c	(revision 11590)
+++ src/basemath/arith2.c	(working copy)
@@ -862,10 +862,10 @@
       GEN f=factor_Aurifeuille_prime(p,d[i]);
       B = Z_factor(f);
       A = merge_factor(A, B, (void*)&cmpii, cmp_nodata);
-      B = Z_factor(diviiexact(polcyclo_eval(d[i],p), f));
+      B = Z_factor(diviiexact(polcyclo_eval(utoipos(d[i]), p), f));
     }
     else
-      B = Z_factor(polcyclo_eval(d[i],p));
+      B = Z_factor(polcyclo_eval(utoipos(d[i]), p));
     A = merge_factor(A, B, (void*)&cmpii, cmp_nodata);
   }
   return gerepilecopy(av, A);
Index: src/basemath/bibli2.c
===================================================================
--- src/basemath/bibli2.c	(revision 11590)
+++ src/basemath/bibli2.c	(working copy)
@@ -389,56 +389,126 @@
   return gerepilecopy(av, RgX_inflate(T,q));
 }
 
-/* cyclotomic polynomial */
+/* n-th cyclotomic polynomial evaluated in x */
 GEN
-polcyclo_eval(long n, GEN x)
+polcyclo_eval(GEN n, GEN x)
 {
-  pari_sp av= avma;
-  GEN P, md, xd, yn, yd;
-  long l, s, i, j, q, mu, tx;
+  pari_sp av = avma;
+  GEN P, s;
+  GEN xd, md, ypos, yneg;
+  long l, i, j, tx, root_of_unity = 0;
+  GEN X = NULL;  /* Assign NULL to prevent compiler warning */
 
-  if (!x) return polcyclo(n, 0);
-  tx = typ(x);
-  if (gcmpX(x)) return polcyclo(n, varn(x));
-  if (n <= 0) pari_err(talker, "argument must be positive in polcyclo");
-  if (n == 1) return gsubgs(x, 1);
-  /* n >= 2 */
-  P = gel(factoru(n), 1); l = lg(P)-1;
-  s = P[1]; for (i = 2; i <= l; i++) s *= P[i];
-  q = n/s;
-  if (tx == t_INT && is_pm1(x))
+  if (typ(n) != t_INT) pari_err(typeer, "polcyclo_eval");
+
+  if (!x) return polcyclo(itos(n), 0);
+  if (gcmpX(x)) return polcyclo(itos(n), varn(x));  /* x is a variable */
+
+  /* Special case n <= 1, so in what follows
+   * we may assume that n >= 2. */
+  if (signe(n) <= 0) pari_err(talker, "argument must be positive in polcyclo_eval");
+  if (cmpiu(n,1) == 0) return gsubgs(x, 1);  /* n == 1: return x - 1 */
+
+  /* Factor n and let s be the largest squarefree divisor of n */
+  P = gel(Z_factor(n), 1); l = lg(P)-1;
+  s = gel(P,1); for (i = 2; i <= l; i++) s = mulii(s, gel(P,i));
+
+  /* If n is not equal to s then n is not squarefree.  By replacing x
+   * with x^(n/s), we may assume that n is squarefree in what follows. */
+  if (cmpii(n,s) != 0)
   {
-    avma = av;
-    if (signe(x) > 0 || !odd(q)) return l == 1? utoipos(P[1]): gen_1;
-    /* return Phi_s(-1) */
-    if (n == 2) return gen_0;
-    if (!odd(n) && l == 2) return utoipos(P[2]);
-    return gen_1;
+    x = powgi(x, divii(n,s));
+    n = s;
   }
-  if (q != 1) { x = gpowgs(x, q); n = s; } /* replace n by squarefree part */
-  if (tx == t_POL || tx == t_MAT || lg(x) > n)
-    return gerepileupto(av, poleval(polcyclo(n,0), x));
 
-  xd = cgetg((1<<l) + 1, t_VEC); /* the x^d, where d | n */
+  /* Handle the case x == 1 */
+  if (gcmp1(x))
+  {
+    /* n is prime => return n.
+     * However, we multiply with x to keep the type.
+     * For example if x is Mod(1,11) we want to return Mod(n, 11). */
+    if (l == 1) return gerepileupto(av, gmul(x, n));
+
+    /* n is composite => return 1. */
+    avma = av; return x;
+  }
+
+  /* Heuristics: would it be more efficient to compute the n-th
+   * cyclotomic polynomial and then substituting the GEN x? */
+  tx = typ(x);
+  if (tx == t_POL || tx == t_MAT || cmpui(lg(x), n) > 0)
+    return gerepileupto(av, poleval(polcyclo(itos(n), 0), x));
+
+  xd = cgetg((1<<l) + 1, t_VEC);      /* the x^d, where d | n */
   md = cgetg((1<<l) + 1, t_VECSMALL); /* the mu(d), where d | n */
+
+  /* We compute the n-th cyclotomic polynomial using the formula
+   * Prod_{d|n} (x^d-1)^mu(n/d).
+   * Note that n is squarefree, therefore mu(n/d) is either 1 or -1.
+   * We store the factors with mu(d) == 1 in ypos, the factors with
+   * mu(d) == -1 in yneg.  Then at the end we return ypos/yneg if
+   * mu(n) == 1 and yneg/ypos if mu(n) == -1. */
+  /* Initialize (handle the divisor d == 1) */
   gel(xd, 1) = x;
   md[1] = 1;
-  mu = odd(l)? -1: 1; /* mu(n) */
-  if (mu == 1) { yd = gen_1; yn = gsubgs(x,1); }
-  else         { yn = gen_1; yd = gsubgs(x,1); }
-  for (i = 1; i <= l; i++) /* compute Prod_{d|n} (x^d-1)^mu(n/d) */
+  ypos = gsubgs(x,1);
+  yneg = gen_1;
+
+  for (i = 1; i <= l; i++)
   {
-    long ti = 1<<(i-1), p = P[i];
+    long ti = 1<<(i-1);
+    GEN p = gel(P,i);
     for (j = 1; j <= ti; j++) {
-      GEN X = gpowgs(gel(xd,j), p);
-      gel(xd,ti+j) = X;
+      X = powgi(gel(xd,j), p);
+      gel(xd, ti+j) = X;
       md[ti+j] = -md[j];
-      if (mu == md[ti+j]) yn = gmul(yn, gsubgs(X,1));
-      else                yd = gmul(yd, gsubgs(X,1));
+      /* Handle the case where X == 1, i.e. when x^d = 1.
+       * We let the variable root_of_unity denote the smallest index
+       * ti+j such that X == 1.  If no such index exists, then
+       * root_of_unity remains 0.  Furthermore, we do not multiply
+       * with X - 1 if X == 1, we handle these factors at the end. */
+      if (gcmp1(X))
+      {
+        if (root_of_unity == 0) root_of_unity = ti+j;
+      }
+      else
+      {
+      	if (md[ti+j] == 1) ypos = gmul(ypos, gsubgs(X, 1));
+      	else               yneg = gmul(yneg, gsubgs(X, 1));
+      }
     }
   }
-  return gerepileupto(av, gdiv(yn,yd));
+
+  /* Compute result in ypos */
+  if (odd(l))
+    ypos = gdiv(yneg, ypos);  /* mu(n) == -1 */
+  else
+    ypos = gdiv(ypos, yneg);  /* mu(n) == 1 */
+
+  if (root_of_unity)
+  {
+    /* x is a root of unity.  If root_of_unity == (1<<l), then x was a
+     * primitive n-th root of unity, so the result should be zero.
+     * We return x^n - 1 to preserve the type (instead of returning gen_0).
+     * Note that X is equal to x^n, which should be (inexactly) 1. */
+    if (root_of_unity == (1<<l)) return gerepileupto(av, gsubgs(X, 1));
+
+    /* x is a primitive d-th root of unity, where d|n and d<n.
+     * We first multiply ypos by X = x^n = 1 to preserve the type. */
+    ypos = gmul(ypos, X);
+
+    /* Check bits of root_of_unity.  If root_of_unity = (1<<l) - (1<<(i-1))
+     * for some i < l, then n/d == P[i] and we need to multiply with P[i];
+     * otherwise n/d is composite and nothing more needs to be done. */
+    for (i = 1; i <= l; i++)
+    {
+      if (root_of_unity + (1<<(i-1)) == (1<<l)) {ypos = gmul(ypos, gel(P,i)); break;}
+    }
+  }
+
+  return gerepileupto(av, ypos);
 }
+
 /********************************************************************/
 /**                                                                **/
 /**                  HILBERT & PASCAL MATRICES                     **/
Index: src/functions/polynomials/polcyclo
===================================================================
--- src/functions/polynomials/polcyclo	(revision 11590)
+++ src/functions/polynomials/polcyclo	(working copy)
@@ -1,20 +1,20 @@
 Function: polcyclo
 Section: polynomials
 C-Name: polcyclo_eval
-Prototype: LDG
+Prototype: GDG
 Help: polcyclo(n,{a = 'x}): n-th cyclotomic polynomial evaluated at a.
 Description:
-  (small,?var):gen     polcyclo($1,$2)
-  (small,gen):gen      polcyclo_eval($1,$2)
+  (gen,?var):gen     polcyclo($1,$2)
+  (gen,gen):gen      polcyclo_eval($1,$2)
 Doc: $n$-th cyclotomic polynomial, evaluated at $a$ (\kbd{'x} by default). The
  integer $n$ must be positive.
 
  Algorithm used: reduce to the case where $n$ is squarefree; to compute the
  cyclotomic polynomial, use $\Phi_{np}(x)=\Phi_n(x^p)/\Phi(x)$; to compute
  it evaluated, use $\Phi_n(x) = \prod_{d\mid n} (x^d-1)^{\mu(n/d)}$. In the
- evaluated case, the algorithm can deal with all rational values $a$;
- otherwise it assumes that $a^d - 1$ is invertible for all $d\mid n$. If this
- is not the case, use \kbd{subst(polcyclo(n),x,a)}.
+ evaluated case, the algorithm assumes the following: for all $d\mid n$,
+ either $a^d = 1$ or $a^d - 1$ is invertible.
+ If this is not the case, use \kbd{subst(polcyclo(n),x,a)}.
 
 Variant: The variant \fun{GEN}{polcyclo}{long n, long v} returns the $n$-th
  cyclotomic polynomial in variable $v$.
Index: src/headers/paridecl.h
===================================================================
--- src/headers/paridecl.h	(revision 11590)
+++ src/headers/paridecl.h	(working copy)
@@ -1074,7 +1074,7 @@
 GEN     convol(GEN x, GEN y);
 int     gen_cmp_RgX(void *data, GEN x, GEN y);
 GEN     polcyclo(long n, long v);
-GEN     polcyclo_eval(long n, GEN x);
+GEN     polcyclo_eval(GEN n, GEN x);
 GEN     dirdiv(GEN x, GEN y);
 GEN     dirmul(GEN x, GEN y);
 GEN     gen_indexsort(GEN x, void *E, int (*cmp)(void*,GEN,GEN));