Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - factcyclo.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.2 lcov report (development 29115-f22e516b23) Lines: 915 1077 85.0 %
Date: 2024-03-28 08:06:56 Functions: 77 85 90.6 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000, 2012  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /* written by Takashi Fukuda
      15             :  *  2019.10.27 : Flx_factcyclo_gen, FpX_factcyclo_gen
      16             :  *  2019.10.28 : Flx_factcyclo_lift, FpX_factcyclo_lift
      17             :  *  2019.11.3  : Flx_factcyclo_newton, FpX_factcyclo_newton
      18             :  *  2019.11.12 : gausspol for prime
      19             :  *  2019.11.13 : gausspol for prime power
      20             :  *  2019.11.14 : ZpX_roots_nonsep with ZX_Zp_root
      21             :  *  2019.11.15 : test ZpX_roots_nonsep with polrootspadic
      22             :  *  2019.11.16 : accept variable number
      23             :  *  2019.11.17 : gen_ascent
      24             :  *  2019.11.20 : ZpX_roots_nonsep with FpX_roots
      25             :  *  2021.7.19  : Flx_factcyclo_newton_general
      26             :  *  2021.7.22  : Flx_conductor_lift */
      27             : 
      28             : #include "pari.h"
      29             : #include "paripriv.h"
      30             : 
      31             : #define DEBUGLEVEL DEBUGLEVEL_factcyclo
      32             : 
      33             : #define GENERAL            1
      34             : #define NEWTON_POWER       2
      35             : #define NEWTON_GENERAL     4
      36             : #define NEWTON_GENERAL_NEW 8
      37             : #define ASCENT            16
      38             : 
      39             : #define Flx_polcyclo(n, p) ZX_to_Flx(polcyclo(n, 0), p)
      40             : #define FpX_polcyclo(n, p) FpX_red(polcyclo(n, 0), p)
      41             : 
      42             : /* 0 <= z[i] <= ULONG_MAX */
      43             : static GEN
      44           0 : ZX_to_nx(GEN z)
      45             : {
      46           0 :   long i, r = lg(z);
      47           0 :   GEN x = cgetg(r, t_VECSMALL);
      48           0 :   for (i = 2; i < r; i++) x[i] = itou(gel(z, i));
      49           0 :   return x;
      50             : }
      51             : 
      52             : static long
      53       92826 : QX_den_pval(GEN x, GEN p)
      54             : {
      55       92826 :   long i, vmax = 0, l = lg(x);
      56      334882 :   for (i = 2; i < l; i++)
      57             :   {
      58      242056 :     GEN z = gel(x, i);
      59             :     long v;
      60      242056 :     if (typ(z)==t_FRAC && (v = Z_pval(gel(z, 2), p)) > vmax) vmax = v;
      61             :   }
      62       92826 :   return vmax;
      63             : }
      64             : 
      65             : static long
      66       15807 : QXV_den_pval(GEN vT, GEN kT, GEN p)
      67             : {
      68       15807 :   long k, vmax = 0, l = lg(kT);
      69       75173 :   for (k = 1; k < l; k++)
      70             :   {
      71       59366 :     long v = QX_den_pval(gel(vT, kT[k]), p);
      72       59366 :     if (v > vmax) vmax = v;
      73             :   }
      74       15807 :   return vmax;
      75             : }
      76             : 
      77             : /* n=el^e, p^d=1 (mod n), d is in [1,el-1], phi(n)=d*f.
      78             :  *  E[i] : 1 <= i <= n-1
      79             :  *  E[g^i*p^j mod n] = i+1  (0 <= i <= f-1, 0 <= j <= d-1)
      80             :  *  We use E[i] (1 <= i <= d). Namely i < el. */
      81             : static GEN
      82       33460 : set_E(long pmodn, long n, long d, long f, long g)
      83             : {
      84             :   long i, j;
      85       33460 :   GEN E = const_vecsmall(n-1, 0);
      86       33460 :   pari_sp av = avma;
      87       33460 :   GEN C = Fl_powers(g, f-1, n);
      88      121856 :   for (i = 1; i <= f; i++)
      89             :   {
      90       88396 :     ulong x = C[i];
      91     1239140 :     for (j = 1; j <= d; j++) { x = Fl_mul(x, pmodn, n); E[x] = i; }
      92             :   }
      93       33460 :   return gc_const(av, E);
      94             : }
      95             : 
      96             : /* x1, x2 of the same degree; gcd(p1,p2) = 1, m = p1*p2, m2 = m >> 1*/
      97             : static GEN
      98           0 : ZX_chinese_center(GEN x1, GEN p1, GEN x2, GEN p2, GEN m, GEN m2)
      99             : {
     100           0 :   long i, l = lg(x1);
     101           0 :   GEN x = cgetg(l, t_POL);
     102             :   GEN y1, y2, q1, q2;
     103           0 :   (void)bezout(p1, p2, &q1, &q2);
     104           0 :   y1 = Fp_mul(p2, q2, m);
     105           0 :   y2 = Fp_mul(p1, q1, m);
     106           0 :   for (i = 2; i < l; i++)
     107             :   {
     108           0 :     GEN y = Fp_add(mulii(gel(x1, i), y1), mulii(gel(x2, i), y2), m);
     109           0 :     if (cmpii(y, m2) > 0) y = subii(y, m);
     110           0 :     gel(x, i) = y;
     111             :   }
     112           0 :   x[1] = x1[1]; return x;
     113             : }
     114             : 
     115             : /* find n_el primes el such that el=1 (mod n) and el does not divide d(T) */
     116             : static GEN
     117       98534 : list_el_n(ulong el0, ulong n, GEN d1, long n_el)
     118             : {
     119       98534 :   GEN v = cgetg(n_el + 1, t_VECSMALL);
     120             :   forprime_t T;
     121             :   long i;
     122       98534 :   u_forprime_arith_init(&T, el0+n, ULONG_MAX, 1, n);
     123      247044 :   for (i = 1; i <= n_el; i++)
     124             :   {
     125             :     ulong el;
     126      148510 :     do el = u_forprime_next(&T); while (dvdiu(d1, el));
     127      148510 :     v[i] = el;
     128             :   }
     129       98534 :   return v;
     130             : }
     131             : 
     132             : /* return min el s.t. 2^63<el and el=1 (mod n). */
     133             : static ulong
     134       49267 : start_el_n(ulong n)
     135             : {
     136       49267 :   ulong MAXHLONG = 1UL<<(BITS_IN_LONG-1), el = (MAXHLONG/n)*n + 1;
     137       49267 :   if ((el&1)==0) el += n; /* if el is even, then n is odd */
     138       49267 :   return el + (n << 1);
     139             : }
     140             : 
     141             : /* start probably catches d0*T_k(x). So small second is enough. */
     142             : static ulong
     143       49267 : get_n_el(GEN d0, ulong *psec)
     144             : {
     145       49267 :   ulong start = ((lgefint(d0)-2)*BITS_IN_LONG)/(BITS_IN_LONG-1)+1, second = 1;
     146       49267 :   if (start>10) second++;
     147       49267 :   if (start>100)  { start++; second++; }
     148       49267 :   if (start>500)  { start++; second++; }
     149       49267 :   if (start>1000) { start++; second++; }
     150       49267 :   *psec = second; return start;
     151             : }
     152             : 
     153             : static long
     154       77072 : FpX_degsub(GEN P, GEN Q, GEN p)
     155             : {
     156       77072 :   pari_sp av = avma;
     157       77072 :   long d = degpol(FpX_sub(P, Q, p));
     158       77072 :   return gc_long(av, d);
     159             : }
     160             : 
     161             : /* n=odd prime power, F=Q(zeta_n), G=G(F/Q)=(Z/nZ)^*, H=<p>, K <--> H,
     162             :  * t_k=Tr_{F/K}(zeta_n^k), T=minimal pol. of t_1 over Q.
     163             :  * g is a given generator of G(K/Q).
     164             :  * There exists a unique G(x) in Q[x] s.t. t_g=G(t_1).
     165             :  * return G(x) mod el assuming that el does not divide d(T), in which case
     166             :  * T(x) is separable over F_el and so Vandermonde mod el is regular. */
     167             : static GEN
     168      100456 : gausspol_el(GEN H, ulong n, ulong d, ulong f, ulong g, ulong el)
     169             : {
     170      100456 :   ulong j, k, z_n = rootsof1_Fl(n, el);
     171      100456 :   GEN vz_n, L = cgetg(1+f, t_VECSMALL), x = cgetg(1+f, t_VECSMALL), X;
     172             : 
     173      100456 :   vz_n = Fl_powers(z_n, n-1, el)+1;
     174      366316 :   for (k = 0; k < f; k++)
     175             :   {
     176      265860 :     ulong gk = Fl_powu(g, k, n), t = 0;
     177     3725628 :     for (j = 1; j <= d; j++)
     178     3459768 :       t = Fl_add(t, vz_n[Fl_mul(H[j], gk, n)], el);
     179      265860 :     L[1+k] = t;
     180      265860 :     x[1+(k+f-1)%f] = t;
     181             :   }
     182      100456 :   X = Flv_invVandermonde(L, 1, el);
     183      100456 :   return Flv_to_Flx(Flm_Flc_mul(X, x, el), 0);
     184             : }
     185             : 
     186             : static GEN
     187       66920 : get_G(GEN H, GEN d0, GEN d1, GEN N, long k, ulong *pel, GEN *pM)
     188             : {
     189       66920 :   long n = N[1], d = N[2], f = N[3], g = N[4], i;
     190       66920 :   GEN POL = cgetg(1+k, t_VEC), EL, G, M, x;
     191             :   pari_timer ti;
     192             : 
     193       66920 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
     194       66920 :   EL = list_el_n(*pel, n, d1, k);
     195      167376 :   for (i = 1; i <= k; i++)
     196             :   {
     197      100456 :     ulong el = uel(EL,i);
     198      100456 :     x = gausspol_el(H, n, d, f, g, el);
     199      100456 :     gel(POL, i) = Flx_Fl_mul(x, umodiu(d0, el), el);
     200             :   }
     201       66920 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "get_G : make data k=%ld",k);
     202       66920 :   G = nxV_chinese_center(POL, EL, &M);
     203       66920 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "get_G : nxV_chinese_center k=%ld",k);
     204       66920 :   *pel = EL[k]; *pM = M; return G;
     205             : }
     206             : 
     207             : static long
     208           0 : Q_size(GEN z)
     209             : {
     210           0 :   if (typ(z)==t_INT) return lgefint(z) - 2;
     211           0 :   return maxss(lgefint(gel(z,1)), lgefint(gel(z,2))) - 2;
     212             : }
     213             : /* return max log_a(|x[i]|), a=2^(BITS_IN_LONG-1) */
     214             : static long
     215           0 : ZX_size(GEN x)
     216             : {
     217           0 :   long i, l = lg(x), max = 0;
     218           0 :   for (i = 2; i < l; i++)
     219             :   {
     220           0 :     long y = lgefint(gel(x,i)) - 2;
     221           0 :     if (y > max) max = y;
     222             :   }
     223           0 :   return max;
     224             : }
     225             : 
     226             : /* d0 is a multiple of (O_K:Z[t_1]). i.e. d0*T_k(x) is in Z[x].
     227             :  * d1 has the same prime factors as d(T); d0 d1 = d(T)^2 */
     228             : static GEN
     229       49267 : get_d0_d1(GEN T, GEN P)
     230             : {
     231       49267 :   long i, l = lg(P);
     232       49267 :   GEN x, y, dT = ZX_disc(T);
     233       49267 :   x = y = dT; setsigne(dT, 1);
     234      114451 :   for (i = 1; i < l; i++)
     235       65184 :     if (odd(Z_lvalrem(dT, P[i], &dT)))
     236             :     {
     237       50739 :       x = diviuexact(x, P[i]);
     238       50739 :       y = muliu(y, P[i]);
     239             :     }
     240       49267 :   return mkvec2(sqrti(x), sqrti(y));  /* x and y are square */
     241             : }
     242             : 
     243             : static GEN
     244       33460 : gausspol(GEN T, GEN H, GEN N, ulong d, ulong f, ulong g)
     245             : {
     246       33460 :   long n = N[1], el0 = N[2];
     247       33460 :   GEN F, G1, M1, d0, d1, Data, d0d1 = get_d0_d1(T, mkvecsmall(el0));
     248             :   ulong el, n_el, start, second;
     249             :   pari_timer ti;
     250             : 
     251       33460 :   d0 = gel(d0d1,1); /* d0*F is in Z[X] */
     252       33460 :   d1 = gel(d0d1,2); /* d1 has same prime factors as disc(T) */
     253       33460 :   Data = mkvecsmall4(n, d, f, g);
     254       33460 :   start = get_n_el(d0, &second);
     255       33460 :   el = start_el_n(n);
     256             : 
     257       33460 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
     258       33460 :   if (DEBUGLEVEL == 2) err_printf("gausspol:start=(%ld,%ld)\n",start,second);
     259       33460 :   G1 = get_G(H, d0, d1, Data, start, &el, &M1);
     260       33460 :   for (n_el=second; n_el; n_el++)
     261             :   {
     262             :     GEN m, G2, M2;
     263       33460 :     G2 = get_G(H, d0, d1, Data, n_el, &el, &M2);
     264       33460 :     if (FpX_degsub(G1, G2, M2) < 0) break;  /* G1 = G2 (mod M2) */
     265           0 :     if (DEBUGLEVEL == 2)
     266           0 :       err_printf("G1:%ld, G2:%ld\n",ZX_size(G1), ZX_size(G2));
     267           0 :     if (DEBUGLEVEL >= 6) timer_start(&ti);
     268           0 :     m = mulii(M1, M2);
     269           0 :     G2 = ZX_chinese_center(G1, M1, G2, M2, m, shifti(m,-1));
     270           0 :     if (DEBUGLEVEL >= 6) timer_printf(&ti, "ZX_chinese_center");
     271           0 :     G1 = G2; M1 = m;
     272             :   }
     273       33460 :   F = RgX_Rg_div(G1, d0);
     274       33460 :   if (DEBUGLEVEL == 2)
     275           0 :     err_printf("G1:%ld, d0:%ld, M1:%ld, F:%ld\n",
     276             :                ZX_size(G1), Q_size(d0), Q_size(M1), ZX_size(F));
     277       33460 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "gausspol");
     278       33460 :   return F;
     279             : }
     280             : 
     281             : /* Data = [H, GH, i_t, d0d1, kT, [n, d, f, n_T, mitk]]
     282             :  * v_t_el[k]=t_k mod el, 1<= k <= n-1 */
     283             : static GEN
     284       48054 : mk_v_t_el(GEN vT, GEN Data, ulong el)
     285             : {
     286       48054 :   pari_sp av = avma;
     287       48054 :   GEN H = gel(Data, 1), GH = gel(Data,2), i_t = gel(Data, 3), N=gel(Data, 6);
     288       48054 :   ulong n = N[1],  d = N[2], mitk = N[5], f = N[3], i, k;
     289       48054 :   ulong z_n = rootsof1_Fl(n, el);
     290       48054 :   GEN vz_n = Fl_powers(z_n, n-1, el)+1;
     291       48054 :   GEN v_t_el = const_vecsmall(n-1, 0);
     292             : 
     293      348257 :   for (k = 1; k <= mitk; k++)
     294             :   {
     295      300203 :     if (k > 1 && !isintzero(gel(vT, k))) continue; /* k=1 is always handled */
     296     2222451 :     for (i=1; i<=f; i++)
     297             :     {
     298     1922248 :       ulong j, t = 0, x = Fl_mul(k, GH[i], n);
     299     1922248 :       long y = i_t[x]; /* x!=0, y!=0 */
     300     1922248 :       if (v_t_el[y]) continue;
     301     5435404 :       for (j = 1; j <= d; j++) t = Fl_add(t, vz_n[Fl_mul(x, H[j], n)], el);
     302      262980 :       v_t_el[y] = t;
     303             :     }
     304             :   }
     305       48054 :   return gerepileuptoleaf(av, v_t_el);
     306             : }
     307             : 
     308             : /* G=[[G_1,...,G_d],M,el]
     309             :  * Data = [H, GH, i_t, d0d1, kT, [n, d, f, n_T, mitk]] */
     310             : static GEN
     311       31614 : get_vG(GEN vT, GEN Data, long n_el, ulong *pel, GEN *pM)
     312             : {
     313       31614 :   GEN GH = gel(Data, 2), i_t = gel(Data, 3);
     314       31614 :   GEN d0 = gmael(Data, 4, 1), d1 = gmael(Data, 4, 2);
     315       31614 :   GEN kT = gel(Data, 5), N = gel(Data, 6);
     316       31614 :   long n = N[1], f = N[3], n_T = N[4], mitk = N[5];
     317       31614 :   GEN G = const_vec(mitk, gen_0), vPOL = cgetg(1+mitk, t_VEC);
     318             :   GEN EL, M, X, v_t_el;
     319       31614 :   GEN L = cgetg(1+f, t_VECSMALL), x = cgetg(1+f, t_VECSMALL);
     320             :   long i, j, k;
     321             : 
     322      216870 :   for (k=1; k<=mitk; k++) gel(vPOL, k) = cgetg(1+n_el, t_VEC);
     323       31614 :   EL = list_el_n(*pel, n, d1, n_el);
     324       79668 :   for (i=1; i<=n_el; i++)
     325             :   {
     326       48054 :     ulong el = uel(EL,i), d0model = umodiu(d0, el);
     327       48054 :     v_t_el = mk_v_t_el(vT, Data, el);
     328      214716 :     for (j = 1; j <= f; j++) L[j] = v_t_el[i_t[GH[j]]];
     329       48054 :     X = Flv_invVandermonde(L, 1, el);
     330      229086 :     for (k = 1; k <= n_T; k++)
     331             :     {
     332             :       GEN y;
     333      181032 :       long itk = kT[k];
     334      181032 :       if (!isintzero(gel(vT, itk))) continue;
     335      700581 :       for (j = 1; j <= f; j++) x[j] = v_t_el[i_t[Fl_mul(itk, GH[j], n)]];
     336      133705 :       y = Flv_to_Flx(Flm_Flc_mul(X, x, el), 0);
     337      133705 :       gmael(vPOL, itk, i) = Flx_Fl_mul(y, d0model, el);
     338             :     }
     339             :   }
     340      150346 :   for (k = 1; k <= n_T; k++)
     341             :   {
     342      118732 :     long itk = kT[k];
     343      118732 :     if (!isintzero(gel(vT, itk))) continue;
     344       87224 :     gel(G, itk) = nxV_chinese_center(gel(vPOL, itk), EL, &M);
     345             :   }
     346       31614 :   *pel = EL[n_el]; *pM = M; return G;
     347             : }
     348             : 
     349             : /* F=Q(zeta_n), H=<p> in (Z/nZ)^*, K<-->H, t_k=Tr_{F/K}(zeta_n^k).
     350             :  * i_t[i]=k ==> iH=kH, i.e. t_i=t_k. We use t_k instead of t_i:
     351             :  * the number of k << the number of i. */
     352             : static GEN
     353       15807 : get_i_t(long n, long p)
     354             : {
     355             :   long a;
     356       15807 :   GEN v_t = const_vecsmall(n-1, 0);
     357       15807 :   GEN i_t = cgetg(n, t_VECSMALL);  /* access i_t[a] with 1 <= a <= n-1 */
     358      118319 :   for (a = 1; a < n; a++)
     359             :   {
     360             :     long b;
     361      854693 :     while (a < n && v_t[a]) a++;
     362      118319 :     if (a==n) break;
     363      102512 :     b = a;
     364      838886 :     do { v_t[b] = 1; i_t[b] = a; b = Fl_mul(b, p, n); } while (b != a);
     365             :   }
     366       15807 :   return i_t;
     367             : }
     368             : 
     369             : /* get T_k(x) 1<= k <= d. d0*T_k(x) is in Z[x].
     370             :  * T_0(x)=T_n(x)=f.
     371             :  * Data = [H, GH, i_t, d0d1, kT, [n, d, f, n_T, mitk]] */
     372             : static GEN
     373       15807 : get_vT(GEN Data, int NEW)
     374             : {
     375       15807 :   pari_sp av = avma;
     376       15807 :   GEN d0 = gmael(Data, 4, 1), kT = gel(Data, 5), N = gel(Data, 6);
     377       15807 :   ulong k, n = N[1], n_T = N[4], mitk = N[5];
     378       15807 :   GEN G1, M1, vT = const_vec(mitk, gen_0); /* vT[k]!=NULL ==> vT[k]=T_k */
     379       15807 :   ulong n_k = 0, el, n_el, start, second;
     380             :   pari_timer ti;
     381             : 
     382       15807 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
     383       15807 :   if (!NEW) { gel(vT, 1) = pol_x(0); n_k++; }
     384       15807 :   start = get_n_el(d0, &second);
     385       15807 :   el = start_el_n(n);
     386       15807 :   if (DEBUGLEVEL == 2) err_printf("get_vT: start=(%ld,%ld)\n",start,second);
     387       15807 :   G1 = get_vG(vT, Data, start, &el, &M1);
     388       15807 :   for (n_el = second;; n_el++)
     389           0 :   {
     390             :     GEN G2, M2, m, m2;
     391       15807 :     G2 = get_vG(vT, Data, n_el, &el, &M2);
     392       15807 :     m = mulii(M1, M2); m2 = shifti(m,-1);
     393       75173 :     for (k = 1; k <= n_T; k++)
     394             :     {
     395       59366 :       long j = kT[k];
     396       59366 :       if (!isintzero(gel(vT,j))) continue;
     397       43612 :       if (FpX_degsub(gel(G1,j), gel(G2,j), M2) < 0) /* G1=G2 (mod M2) */
     398             :       {
     399       43612 :         gel(vT,j) = RgX_Rg_div(gel(G1,j), d0);
     400       43612 :         n_k++;
     401       43612 :         if (DEBUGLEVEL == 2)
     402           0 :           err_printf("G1:%ld, d0:%ld, M1:%ld, vT[%ld]:%ld words\n",
     403           0 :             ZX_size(gel(G1,j)), Q_size(d0), Q_size(M1), j, ZX_size(gel(vT,j)));
     404             :       }
     405             :       else
     406             :       {
     407           0 :         if (DEBUGLEVEL == 2)
     408           0 :           err_printf("G1:%ld, G2:%ld\n", ZX_size(gel(G1,j)),ZX_size(gel(G2,j)));
     409           0 :         gel(G1, j) = ZX_chinese_center(gel(G1,j),M1, gel(G2,j),M2, m,m2);
     410             :       }
     411             :     }
     412       15807 :     if (n_k == n_T) break;
     413           0 :     M1 = m;
     414             :   }
     415       15807 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "get_vT");
     416       15807 :   return gerepilecopy(av, vT);
     417             : }
     418             : 
     419             : /* return sorted kT={i_t[k] | 1<=k<=d}
     420             :  * {T_k(x) | k in kT} are all the different T_k(x) (1<=k<=d) */
     421             : static GEN
     422       15754 : get_kT(GEN i_t, long d)
     423       15754 : { return vecsmall_uniq(vecsmall_shorten(i_t, d)); }
     424             : 
     425             : static GEN
     426          53 : get_kT_all(GEN GH, GEN i_t, long n, long d, long m)
     427             : {
     428          53 :   long i, j, k = 0;
     429          53 :   GEN x = const_vecsmall(d*m, 0);
     430         106 :   for (i = 1; i <= m; i++)
     431         755 :     for (j = 1; j <= d; j++) x[++k] = i_t[Fl_mul(GH[i], j, n)];
     432          53 :   return vecsmall_uniq(x);
     433             : }
     434             : 
     435             : static GEN
     436          53 : kT_from_kt_new(GEN gGH, GEN kt, GEN i_t, long n)
     437             : {
     438          53 :   GEN EL = gel(gGH, 1);
     439          53 :   long i, k = 0, lEL = lg(EL), lkt = lg(kt);
     440          53 :   GEN x = cgetg(lEL+lkt, t_VECSMALL);
     441             : 
     442         151 :   for (i = 1; i < lEL; i++) x[++k] = i_t[EL[i]];
     443         424 :   for (i = 2; i < lkt; i++) if (n%kt[i]==0) x[++k] = kt[i];
     444          53 :   setlg(x, k+1); return vecsmall_uniq(x);
     445             : }
     446             : 
     447             : static GEN
     448          53 : get_kTdiv(GEN kT, long n)
     449             : {
     450          53 :   long i, k = 0, l = lg(kT);
     451          53 :   GEN div = cgetg(l, t_VECSMALL);
     452         202 :   for (i = 1; i < l; i++) if (n%kT[i]==0) div[++k] = kT[i];
     453          53 :   setlg(div, k+1);
     454          53 :   return div;
     455             : }
     456             : 
     457             : /* T is separable over Zp but not separable over Fp.
     458             :  * receive all roots mod p^s and return all roots mod p^(s+1) */
     459             : static GEN
     460         833 : ZpX_roots_nonsep(GEN T, GEN R0, GEN p, GEN ps, GEN ps1)
     461             : {
     462         833 :   long i, j, n = 0, lr = lg(R0);
     463         833 :   GEN v = cgetg(lr, t_VEC), R;
     464        4375 :   for (i = 1; i < lr; i++)
     465             :   {
     466        3542 :     GEN z, f = ZX_unscale_div(ZX_translate(T, gel(R0, i)), ps);
     467        3542 :     (void)ZX_pvalrem(f, p, &f);
     468        3542 :     gel(v, i) = z = FpX_roots(f, p);
     469        3542 :     n += lg(z)-1;
     470             :   }
     471         833 :   R = cgetg(n+1, t_VEC); n = 0;
     472        4375 :   for (i = 1; i < lr; i++)
     473             :   {
     474        3542 :     GEN z = gel(v, i);
     475        3542 :     long lz = lg(z);
     476        8190 :     for (j=1; j<lz; j++)
     477        4648 :       gel(R, ++n) = Fp_add(gel(R0, i), mulii(gel(z, j), ps), ps1);
     478             :   }
     479         833 :   return ZV_sort_uniq_shallow(R);
     480             : }
     481             : static GEN
     482       49267 : ZpX_roots_all(GEN T, GEN p, long f, long *ptrs)
     483             : {
     484       49267 :   pari_sp av = avma;
     485             :   pari_timer ti;
     486             :   GEN v, ps, ps1;
     487             :   long s;
     488             : 
     489       49267 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
     490       49267 :   v = FpX_roots(T, p); /* FpX_roots returns sorted roots */
     491       49267 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "FpX_roots, deg=%ld", degpol(T));
     492       49267 :   ps = NULL; ps1 = p;
     493       50100 :   for (s = 1; lg(v) != f+1; s++)
     494             :   {
     495         833 :     ps = ps1; ps1 = mulii(ps1, p); /* p^s, p^(s+1) */
     496         833 :     v = ZpX_roots_nonsep(T, v, p, ps, ps1);
     497         833 :     if (gc_needed(av, 1)) gerepileall(av, 3, &v, &ps, &ps1);
     498             :   }
     499       49267 :   *ptrs = s; return v;
     500             : }
     501             : /* x : roots of T in Zp. r < n.
     502             :  * receive vec of x mod p^r and return vec of x mod p^n.
     503             :  * under the assumtion lg(v)-1==degpol(T), x is uniquely determined by
     504             :  * x mod p^r. Namely, x mod p^n is unique. */
     505             : static GEN
     506        2513 : ZX_Zp_liftroots(GEN T, GEN v, GEN p, long r, long n)
     507             : {
     508        2513 :   long i, l = lg(v);
     509        2513 :   GEN R = cgetg(l, t_VEC);
     510        2513 :   GEN pr = powiu(p, r), pn = powiu(p, n);
     511             :   pari_timer ti;
     512             : 
     513        2513 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
     514        9709 :   for (i=1; i<l; i++)
     515             :   {
     516        7196 :     GEN x = gel(v, i), y, z;
     517        7196 :     GEN f = ZX_unscale_div(ZX_translate(T, x), pr);
     518        7196 :     (void)ZX_pvalrem(f, p, &f);
     519        7196 :     y = FpX_roots(f, p);
     520        7196 :     z = (n==r+1) ? y : ZX_Zp_root(f, gel(y, 1), p, n-r);
     521        7196 :     if (lg(y)!=2 || lg(z)!=2)
     522           0 :       pari_err_BUG("ZX_Zp_liftroots, roots are not separable");
     523        7196 :     gel(R, i) = Fp_add(x, mulii(gel(z, 1), pr), pn);
     524             :   }
     525        2513 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "ZX_Zp_liftroots");
     526        2513 :   return R;
     527             : }
     528             : 
     529             : static GEN
     530       33460 : set_R(GEN T, GEN F, GEN Rs, GEN p, long nf, long r, long s, long u)
     531             : {
     532             :   long i;
     533       33460 :   GEN x, pr = powiu(p, r), prs = powiu(p, r+s), R = cgetg(1+nf, t_VEC), Rrs;
     534       33460 :   Rrs = r ? ZX_Zp_liftroots(T, Rs, p, s, r+s) : Rs;
     535       33460 :   x = gel(Rrs, 1);
     536      121856 :   for (i = 1; i <= nf; i++)
     537             :   {
     538       88396 :     x = FpX_eval(F, x, prs);
     539       88396 :     if (r) x = gel(Rrs, ZV_search(Rs, diviiexact(x, pr)));
     540       88396 :     gel(R, i) = x;  /* R[i]=t_1^(g^i), t_1=Rrs[1], mod p^(r+s) */
     541             :   }
     542       33460 :   if (r+s < u) R = ZX_Zp_liftroots(T, R, p, r+s, u);
     543       32116 :   else if (r+s > u) R = FpV_red(R, powiu(p, u));
     544       33460 :   return R;  /* mod p^u, accuracy for pol_newton */
     545             : }
     546             : 
     547             : /* Preparation for factoring polcyclo(el^e) mod p
     548             :  * f_0(x) : irred factor of polcyclo(el^e0) mod p
     549             :  * f_1(x)=f_0(x^(el^e1)) : irred factor of polcyclo(el^e) mod p
     550             :  *
     551             :  * p^d0 = 1 (mod el^e0), p^d = 1 (mod el^e)
     552             :  *
     553             :  * case el=2, 2^s || (p-1), s>=2
     554             :  * d=1 (1 <= e <= s), d=2^(e-s) (s < e)
     555             :  * e0=e, e1=0 if e <= s
     556             :  * e0=s, e1=e-s if s < e
     557             :  * d0=1
     558             :  *
     559             :  * case el=2, 2^s || (p+1), s>=2
     560             :  * d=1 (e == 1), d=2 (2 <= e <= s), d=2^(e-s) (s < e)
     561             :  * e0=e, e1=0 if e <= s+1
     562             :  * e0=s+1, e1=e-s-1 if s+1 < e
     563             :  * d0=1 if e==1,  d0=2 if e>1
     564             :  *
     565             :  * case el>2, el^s || (p^d0-1)
     566             :  * d=d0 (1 <= e <= s), d=d0*el^(e-s) (s < e)
     567             :  * e0=e, e1=0 if e <= s
     568             :  * e0=s, e1=e-s if s < e
     569             :  * d0 is min d s.t. p^d=1 (mod el)
     570             :  *
     571             :  * We do not need d. So d is not returned. */
     572             : static GEN
     573       34222 : set_e0_e1(ulong el, ulong e, GEN p)
     574             : {
     575       34222 :   ulong s, d0, e0, e1, f0, n, phin, g, up = itou_or_0(p);
     576             : 
     577       34222 :   if (el == 2)
     578             :   {
     579           0 :     ulong pmod4 = up ? up&3 : mod4(p);
     580           0 :     if (pmod4 == 3)  /* p+1 = 0 (mod 4) */
     581             :     {
     582           0 :       s = up ? vals(up+1) : vali(addiu(p, 1));
     583           0 :       if (e <= s+1) { e0 = e; e1 = 0;}
     584           0 :       else { e0 = s+1; e1= e-s-1;}
     585           0 :       d0 = e == 1? 1: 2;
     586             :     }
     587             :     else  /* p-1 = 0 (mod 4) */
     588             :     {
     589           0 :       s = up ? vals(up-1) : vali(subiu(p, 1));
     590           0 :       if (e <= s) { e0 = e; e1 = 0; }
     591           0 :       else { e0 = s; e1 = e-s; }
     592           0 :       d0 = 1;
     593             :     }
     594           0 :     phin = 1L<<(e0-1);
     595             :   }
     596             :   else  /* el is odd */
     597             :   {
     598       34222 :     ulong pmodel = up ? up%el : umodiu(p, el);
     599       34222 :     d0 = Fl_order(pmodel, el-1, el);
     600       34222 :     s = Z_lval(subiu(powiu(p, d0), 1), el);
     601       34222 :     if (e <= s) { e0 = e; e1 = 0; } else { e0 = s; e1 = e-s; }
     602       34222 :     phin = (el-1)*upowuu(el, e0-1);
     603             :   }
     604       34222 :   n = upowuu(el, e0); f0 = phin/d0;
     605       34222 :   g = (el==2) ? 1 : pgener_Zl(el);
     606       34222 :   return mkvecsmalln(7, n, e0, e1, phin, g, d0, f0);
     607             : }
     608             : 
     609             : /* return 1 if newton is fast, return 0 if gen is fast */
     610             : static int
     611       78511 : use_newton(long d, long f)
     612             : {
     613       78511 :   if (2*d <= f) return 0;
     614       73603 :   else if (f <= d) return 1;
     615        5935 :   else if (d <= 50) return 0;
     616           0 :   else if (f <= 60) return 1;
     617           0 :   else if (d <= 90) return 0;
     618           0 :   else if (f <= 150) return 1;
     619           0 :   else if (d <= 150) return 0;
     620           0 :   else if (f <= 200) return 1;
     621           0 :   else if (200*d <= f*f) return 0;
     622           0 :   else return 1;
     623             : }
     624             : 
     625             : /* return 1 if newton_general is fast, return 0 otherwise. Assume f > 40 */
     626             : static int
     627          14 : use_newton_general(long d, long f, long maxdeg)
     628             : {
     629          14 :   if (maxdeg < 20) return 0;
     630          14 :   else if (f <= 50) return 1;
     631           7 :   else if (maxdeg < 30) return 0;
     632           7 :   else if (f <= 60) return 1;
     633           7 :   else if (maxdeg < 40) return 0;
     634           7 :   else if (f <= 70) return 1;
     635           7 :   else if (maxdeg < 50) return 0;
     636           7 :   else if (f <= 80) return 1;
     637           7 :   else if (d < 200) return 0;
     638           0 :   else if (f <= 100) return 1;
     639           0 :   else if (d < 300) return 0;
     640           0 :   else if (f <= 120) return 1;
     641           0 :   else if (6*maxdeg < f*f) return 0;
     642           0 :   else return 1;
     643             : }
     644             : 
     645             : static int
     646           7 : use_general(long d, long maxdeg)
     647             : {
     648           7 :   if (d <= 50) return 1;
     649           7 :   else if (maxdeg <= 3*d) return 0;
     650           7 :   else if (d <= 200) return 1;
     651           0 :   else if (maxdeg <= 10*d) return 0;
     652           0 :   else if (d <= 500) return 1;
     653           0 :   else if (maxdeg <= 20*d) return 0;
     654           0 :   else if (d <= 1000) return 1;
     655           0 :   else return 0;
     656             : }
     657             : 
     658             : static void
     659          70 : update_dfm(long *pd, long *pf, long *pm, long di, long fi)
     660             : {
     661          70 :   long c = ugcd(*pd,di), d1 = *pd * di, f1 = *pf * fi;
     662          70 :   *pd = d1 / c; *pf = c * f1; *pm += d1 * d1 * f1;
     663          70 :   if (DEBUGLEVEL == 4) err_printf("(%ld,%ld), ",d1,f1);
     664          70 : }
     665             : /* assume ord(p mod f) > 1 */
     666             : static ulong
     667       94605 : set_action(GEN fn, GEN p, long d, long f)
     668             : {
     669       94605 :   GEN EL = gel(fn, 1), E = gel(fn, 2);
     670       94605 :   long i, d0, f0, m0, m1, maxdeg, max, l = lg(EL);
     671       94605 :   ulong action = 0, up = itou_or_0(p);
     672       94605 :   GEN D = cgetg(l, t_VECSMALL), F = cgetg(l, t_VECSMALL);
     673             : 
     674       94605 :   d += 10*(lgefint(p)-3);
     675       94605 :   if (l == 2)
     676             :   { /* n is a prime power */
     677       47670 :     action |= (EL[1]==2 || !use_newton(d, f))? GENERAL: NEWTON_POWER;
     678       47670 :     return action;
     679             :   }
     680       46935 :   if (f <= 2) action |= NEWTON_GENERAL;
     681       36540 :   else if (d <= 10) action |= GENERAL;
     682        5835 :   else if (f <= 10) action |= NEWTON_GENERAL;
     683         476 :   else if (d <= 20) action |= GENERAL;
     684          73 :   else if (f <= 20) action |= NEWTON_GENERAL_NEW;
     685          27 :   else if (d <= 30) action |= GENERAL;
     686          21 :   else if (f <= 30) action |= NEWTON_GENERAL_NEW;
     687          21 :   else if (d <= 40) action |= GENERAL;
     688          21 :   else if (f <= 40) action |= NEWTON_GENERAL_NEW;
     689       46935 :   if (action) return action;
     690             :   /* can assume that d > 40 and f > 40 */
     691             : 
     692          21 :   maxdeg = max = 1;
     693          91 :   for (i = 1; i < l; i++)
     694             :   {
     695          70 :     long x, el = EL[i], e = E[i];
     696          70 :     long q = upowuu(el, e-1), ni = q * el, phini = ni - q;
     697          70 :     long di = Fl_order(umodiu(p, ni), phini, ni);
     698          70 :     D[i] = di; F[i] = phini / di;
     699          70 :     x = ugcd(max, di); max = max * (di / x); /* = lcm(d1,..di) */
     700          70 :     if (x > 1) maxdeg = max*x;
     701          70 :     if (DEBUGLEVEL == 4) err_printf("(%ld,%ld), ", D[i], F[i]);
     702             :   }
     703          21 :   if (maxdeg == 1) return action;
     704          14 :   if (up != 2)
     705             :   {
     706          14 :     if (use_newton_general(d, f, maxdeg))
     707             :     { /* does not decompose n */
     708           7 :       action |= (20 < d)? NEWTON_GENERAL_NEW: NEWTON_GENERAL;
     709           7 :       return action;
     710             :     }
     711           7 :     if (use_general(d, maxdeg)) action |= GENERAL;
     712             :   }
     713           7 :   if (l < 4) return action; /* n has two factors */
     714             : 
     715           7 :   d0 = f0 = 1; m0 = 0;
     716          42 :   for (i = 1; i < l; i++) update_dfm(&d0, &f0, &m0, D[i], F[i]);
     717           7 :   if (DEBUGLEVEL == 4) err_printf("(d0,f0)=(%ld,%ld)\n",d0,f0);
     718           7 :   d0 = f0 = 1; m1 = 0;
     719          42 :   for (i = l-1; i >= 1; i--) update_dfm(&d0, &f0, &m1, D[i], D[i]);
     720           7 :   if (DEBUGLEVEL == 4) err_printf("(d0,f0)=(%ld,%ld)\n",d0,f0);
     721           7 :   if (DEBUGLEVEL == 4) err_printf("(m0,m1)=(%lu,%lu) %ld\n",m0,m1,m0<=m1);
     722           7 :   if (m0 <= m1) action |= ASCENT;
     723           7 :   return action;
     724             : }
     725             : 
     726             : static GEN
     727       10758 : FpX_Newton_perm(long d, GEN R, GEN v, GEN pu, GEN p)
     728             : {
     729       10758 :   GEN S = cgetg(d+2, t_VEC);
     730             :   long k;
     731       95189 :   gel(S,1) = utoi(d); for (k = 1; k <= d; k++) gel(S, k+1) = gel(R, v[k]);
     732       10758 :   return FpX_red(FpX_fromNewton(RgV_to_RgX(S, 0), pu), p);
     733             : }
     734             : static GEN
     735       81847 : Flx_Newton_perm(long d, GEN R, GEN v, ulong pu, ulong p)
     736             : {
     737       81847 :   GEN S = cgetg(d+2, t_VEC);
     738             :   long k;
     739     1177903 :   S[1] = d; for (k = 1; k <= d; k++) gel(S, k+1) = gel(R, v[k]);
     740       81847 :   return Flx_red(Flx_fromNewton(Flv_to_Flx(S, 0), pu), p);
     741             : }
     742             : 
     743             : /*  Data = [H, GH, i_t, d0, kT, [n, d, f, n_T, mitk]]
     744             :     N2 = [p, pr, pu, pru] */
     745             : static GEN
     746        6866 : FpX_pol_newton_general(GEN Data, GEN N2, GEN vT, GEN x)
     747             : {
     748        6866 :   GEN i_t = gel(Data, 3), kT = gel(Data, 5), N = gel(Data, 6);
     749        6866 :   long k, d = N[2], n_T = N[4], mitk = N[5];
     750        6866 :   GEN p = gel(N2,1), pr = gel(N2,2), pu = gel(N2,3), pru = gel(N2,4);
     751        6866 :   GEN R = cgetg(1+mitk, t_VEC);
     752             : 
     753       32579 :   for (k = 1; k <= n_T; k++)
     754       25713 :     gel(R, kT[k]) = diviiexact(FpX_eval(gel(vT, kT[k]), x, pru), pr);
     755        6866 :   return FpX_Newton_perm(d, R, i_t, pu, p);
     756             : }
     757             : 
     758             : /* n is any integer prime to p, but must be equal to the conductor
     759             :  * of the splitting field of p in Q(zeta_n).
     760             :  * GH=G/H={g_1, ... ,g_f}
     761             :  * Data = [GHgen, GH, fn, p, [n, d, f, m]] */
     762             : static GEN
     763        2414 : FpX_factcyclo_newton_general(GEN Data)
     764             : {
     765        2414 :   GEN GH = gel(Data, 2), fn = gel(Data, 3), p = gel(Data, 4);
     766        2414 :   long n = mael(Data, 5, 1), d = mael(Data, 5, 2), f = mael(Data, 5, 3);
     767        2414 :   long m = mael(Data, 5, 4), pmodn = umodiu(p, n);
     768        2414 :   long i, k, n_T, mitk, r, s = 0, u = 1;
     769             :   GEN vT, kT, H, i_t, T, d0d1, Data2, Data3, R, Rs, pr, pu, pru;
     770             :   pari_timer ti;
     771             : 
     772        2414 :   if (m != 1) m = f;
     773        2414 :   for (pu = p; cmpiu(pu,d) <= 0; u++) pu = mulii(pu, p);  /* d<pu, pu=p^n */
     774             : 
     775        2414 :   H = Fl_powers(pmodn, d-1, n); /* H=<p> */
     776        2414 :   i_t = get_i_t(n, pmodn); /* i_t[1+i]=k ==> t_i=t_k */
     777        2414 :   kT = get_kT(i_t, d); n_T = lg(kT)-1; mitk = kT[n_T];
     778        2414 :   if (DEBUGLEVEL == 2) err_printf("kT=%Ps %ld elements\n",kT,n_T);
     779        2414 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
     780        2414 :   T = galoissubcyclo(utoi(n), utoi(pmodn), 0, 0);
     781        2414 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "galoissubcyclo");
     782        2414 :   d0d1 = get_d0_d1(T, gel(fn,1)); /* d0*T_k(x) is in Z[x] */
     783        2414 :   Data2 = mkvecn(6, H, GH, i_t, d0d1, kT, mkvecsmalln(5, n, d, f, n_T, mitk));
     784        2414 :   vT = get_vT(Data2, 0);
     785        2414 :   if (DEBUGLEVEL == 4) err_printf("vT=%Ps\n",vT);
     786        2414 :   r = QXV_den_pval(vT, kT, p);
     787        2414 :   Rs = ZpX_roots_all(T, p, f, &s);
     788        2414 :   if (DEBUGLEVEL >= 2) err_printf("(u,s,r)=(%ld,%ld,%ld)\n",u,s,r);
     789        2414 :   if (r+u < s) pari_err_BUG("FpX_factcyclo_newton_general (T is not separable mod p^(r+u))");
     790             :   /* R and vT are mod p^(r+u) */
     791        2414 :   R = (r+u==s) ? Rs : ZX_Zp_liftroots(T, Rs, p, s, r+u);
     792        2414 :   pr = powiu(p, r); pru = powiu(p, r+u); /* Usually, r=0, s=1, pr=1, pru=p */
     793       11320 :   for (k = 1; k<=n_T; k++)
     794             :   {
     795        8906 :     long itk = kT[k];
     796        8906 :     GEN z = r? RgX_Rg_mul(gel(vT, itk), pr): gel(vT, itk);
     797        8906 :     gel(vT, itk) = RgX_to_FpX(z, pru);
     798             :   }
     799        2414 :   Data3 = mkvec4(p, pr, pu, pru);
     800        2414 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
     801        9280 :   for (i=1; i<=m; i++)
     802        6866 :     gel(R,i) = FpX_pol_newton_general(Data2, Data3, vT, gel(R,i));
     803        2414 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "FpX_pol_newton_general");
     804        2414 :   return R;
     805             : }
     806             : 
     807             : /* Data = [vT, gGH, Rs, Rrs, i_t, kt, p, pu, pr, prs,
     808             :           [n, r, s, n_t,mitk], div] */
     809             : static void
     810       10937 : Fp_next_gen3(long x, long i, GEN v_t_p, GEN t, GEN Data)
     811             : {
     812       10937 :   GEN vT = gel(Data, 1), gGH = gel(Data, 2), Rs = gel(Data, 3);
     813       10937 :   GEN Rrs = gel(Data, 4), i_t = gel(Data, 5);
     814       10937 :   GEN pu = gel(Data, 8), pr = gel(Data, 9), prs = gel(Data, 10);
     815       10937 :   GEN EL = gel(gGH, 1), E = gel(gGH, 2), div = gel(Data, 12);
     816       10937 :   long n = mael(Data, 11, 1), r = mael(Data, 11, 2), mitk = mael(Data, 11, 5);
     817       10937 :   long j, k, l = lg(EL), ld = lg(div);
     818       10937 :   if (i >= l) return;
     819       14099 :   for (j = 0; j < E[i]; j++)
     820             :   {
     821       10544 :     if (j > 0)
     822             :     {
     823        6989 :       long itx = i_t[x];
     824        6989 :       t = FpX_eval(gel(vT, i_t[EL[i]]), t, prs); /* mod p^(r+s) */
     825        6989 :       if (r) t = gel(Rrs, ZV_search(Rs, diviiexact(t, pr))); /* mod p^(r+s) */
     826        6989 :       if (itx <= mitk) gel(v_t_p, itx) = Fp_red(t, pu); /* mod p^u */
     827       16722 :       for (k = 1; k<ld; k++)
     828             :       {
     829        9733 :         ulong y = Fl_mul(x, div[k], n);
     830        9733 :         long ity = i_t[y];
     831             :         GEN v;
     832        9733 :         if (ity > mitk || !isintzero(gel(v_t_p, ity))) continue;
     833         653 :         v = FpX_eval(gel(vT, i_t[div[k]]), t, prs); /* mod p^(r+s) */
     834         653 :         if (r) v = diviiexact(v, pr); /* mod p^s */
     835         653 :         gel(v_t_p, ity) = Fp_red(v, pu);
     836             :       }
     837             :     }
     838       10544 :     Fp_next_gen3(x, i+1, v_t_p, t, Data);
     839       10544 :     x = Fl_mul(x, EL[i], n);
     840             :   }
     841             : }
     842             : 
     843             : /* Data = [vT, gGH, Rs, Rrs, i_t, kt, p, pu, pr, prs,
     844             :           [n, r, s, n_t, mitk], div] */
     845             : static GEN
     846         393 : Fp_mk_v_t_p3(GEN Data, long i)
     847             : { /* v_t_p[k] != gen_0 => v_t_p[k] = t_k mod p^u */
     848         393 :   GEN Rs = gel(Data, 3), Rrs = gel(Data, 4);
     849         393 :   GEN pu = gel(Data, 8), pr = gel(Data, 9), prs = gel(Data, 10);
     850         393 :   GEN vT = gel(Data, 1), i_t = gel(Data, 5), div = gel(Data, 12);
     851         393 :   long k, r = mael(Data, 11, 2), mitk = mael(Data, 11, 5), ld = lg(div);
     852         393 :   GEN v_t_p = const_vec(mitk, gen_0);
     853             : 
     854         393 :   gel(v_t_p, 1) = Fp_red(gel(Rs, i), pu); /* mod p^u, guaranteed u<=s */
     855         393 :   Fp_next_gen3(1, 1, v_t_p, gel(Rrs, i), Data);
     856         758 :   for (k = 1; k<ld; k++)
     857             :   {
     858         365 :     ulong itk = i_t[div[k]];
     859         365 :     GEN x = FpX_eval(gel(vT, itk), gel(Rrs, i), prs);
     860         365 :     if (r) x = diviiexact(x, pr); /* mod p^s */
     861         365 :     gel(v_t_p, itk) = Fp_red(x, pu);
     862             :   }
     863         393 :   return v_t_p;
     864             : }
     865             : 
     866             : /* Data = [vT, gGH, Rs, Rrs, i_t, kt, p, pu, pr, prs,
     867             :           [n, r, s, n_t,mitk], div] */
     868             : static void
     869       21024 : Fl_next_gen3(long x, long i, GEN v_t_p, ulong t, GEN Data)
     870             : {
     871       21024 :   GEN vT = gel(Data, 1), gGH = gel(Data, 2), Rs = gel(Data, 3);
     872       21024 :   GEN Rrs = gel(Data, 4), i_t = gel(Data, 5);
     873       21024 :   long pu = mael(Data, 8, 2), pr = mael(Data, 9, 2), prs = mael(Data, 10, 2);
     874       21024 :   GEN EL = gel(gGH, 1), E = gel(gGH, 2), div = gel(Data, 12);
     875       21024 :   long n = mael(Data, 11, 1), r = mael(Data, 11, 2), mitk = mael(Data, 11, 5);
     876       21024 :   long j, k, l = lg(EL), ld = lg(div);
     877       21024 :   if (i >= l) return;
     878       27936 :   for (j = 0; j < E[i]; j++)
     879             :   {
     880       20736 :     if (j > 0)
     881             :     {
     882       13536 :       long itx = i_t[x];
     883       13536 :       t = Flx_eval(gel(vT, i_t[EL[i]]), t, prs); /* mod p^(r+s) */
     884       13536 :       if (r) t = Rrs[zv_search(Rs, t/pr)]; /* mod p^(r+s) */
     885       13536 :       if (itx <= mitk) v_t_p[itx] = t%pu; /* mod p^u */
     886       54144 :       for (k = 1; k < ld; k++)
     887             :       {
     888       40608 :         ulong y = Fl_mul(x, div[k], n), v;
     889       40608 :         long ity = i_t[y];
     890       40608 :         if (ity > mitk || v_t_p[ity]) continue;
     891        2592 :         v = Flx_eval(gel(vT, i_t[div[k]]), t, prs); /* mod p^(r+s) */
     892        2592 :         if (r) v /= pr; /* mod p^s */
     893        2592 :         v_t_p[ity] = v%pu;
     894             :       }
     895             :     }
     896       20736 :     Fl_next_gen3(x, i+1, v_t_p, t, Data);
     897       20736 :     x = Fl_mul(x, EL[i], n);
     898             :   }
     899             : }
     900             : 
     901             : /* Data = [vT, gGH, Rs, Rrs, i_t, kt, p, pu, pr, prs,
     902             :           [n, r, s, n_t, mitk], div] */
     903             : static GEN
     904         288 : Fl_mk_v_t_p3(GEN Data, long i)
     905             : { /* v_t_p[k] != 0 => v_t_p[k] = t_k mod p^u */
     906         288 :   GEN Rs = gel(Data, 3), Rrs = gel(Data, 4);
     907         288 :   ulong pu = mael(Data, 8, 2), pr = mael(Data, 9, 2), prs = mael(Data, 10, 2);
     908         288 :   GEN vT = gel(Data, 1), i_t = gel(Data, 5), div = gel(Data, 12);
     909         288 :   long k, r = mael(Data, 11, 2), mitk = mael(Data, 11, 5), ld = lg(div);
     910         288 :   GEN v_t_p = const_vecsmall(mitk, 0);
     911             : 
     912         288 :   v_t_p[1] = Rs[i] % pu; /* mod p^u, guaranteed u<=s */
     913         288 :   Fl_next_gen3(1, 1, v_t_p, Rrs[i], Data);
     914        1152 :   for (k = 1; k < ld; k++)
     915             :   {
     916         864 :     ulong itk = i_t[div[k]], x = Flx_eval(gel(vT, itk), Rrs[i], prs);
     917         864 :     if (r) x /= pr; /* mod p^s */
     918         864 :     v_t_p[itk] = x % pu;
     919             :   }
     920         288 :   return v_t_p;
     921             : }
     922             : 
     923             : /* Data = [GHgen, GH, fn, p, [n, d, f, m]] */
     924             : static GEN
     925          53 : newton_general_new_pre3(GEN Data)
     926             : {
     927          53 :   GEN gGH = gel(Data, 1), GH = gel(Data, 2), fn = gel(Data, 3);
     928          53 :   GEN p = gel(Data, 4);
     929          53 :   long n = mael(Data, 5, 1), d = mael(Data, 5, 2), f = mael(Data, 5, 3);
     930          53 :   long pmodn = umodiu(p, n);
     931          53 :   long k, n_t, n_T, mitk, miTk, r, s = 0, u = 1;
     932             :   GEN vT, kt, kT, H, i_t, T, d0d1, Data2, Rs, Rrs, kTdiv;
     933             :   GEN pr, pu, prs;
     934             :   pari_timer ti;
     935             : 
     936          60 :   for (pu = p; cmpiu(pu,d)<=0; u++) pu = mulii(pu, p);  /* d<pu, pu=p^u */
     937             : 
     938          53 :   H = Fl_powers(pmodn, d-1, n); /* H=<p> */
     939          53 :   i_t = get_i_t(n, pmodn); /* i_t[1+i]=k ==> t_i=t_k */
     940          53 :   kt = get_kT_all(GH, i_t, n, d, 1); n_t = lg(kt)-1; mitk = kt[n_t];
     941          53 :   kT = kT_from_kt_new(gGH, kt, i_t, n); n_T = lg(kT)-1; miTk = kT[n_T];
     942          53 :   kTdiv = get_kTdiv(kT, n);
     943          53 :   if (DEBUGLEVEL == 2)
     944           0 :     err_printf("kt=%Ps %ld elements\nkT=%Ps %ld elements\n",kt,n_t,kT,n_T);
     945          53 :   if (DEBUGLEVEL == 2)
     946           0 :     err_printf("kTdiv=%Ps\n",kTdiv);
     947             : 
     948          53 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
     949          53 :   T = galoissubcyclo(utoi(n), utoi(pmodn), 0, 0);
     950          53 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "galoissubcyclo");
     951          53 :   d0d1 = get_d0_d1(T, gel(fn,1)); /* d0*T_k(x) is in Z[x] */
     952          53 :   Data2 = mkvecn(6, H, GH, i_t, d0d1, kT, mkvecsmalln(5, n, d, f, n_T, miTk));
     953          53 :   vT = get_vT(Data2, 1);
     954          53 :   if (DEBUGLEVEL == 4) err_printf("vT=%Ps\n",vT);
     955          53 :   r = QXV_den_pval(vT, kT, p);
     956          53 :   Rs = ZpX_roots_all(T, p, f, &s);
     957          53 :   if (DEBUGLEVEL >= 2) err_printf("(u,s,r)=(%ld,%ld,%ld)\n",u,s,r);
     958          53 :   if (s < u)
     959             :   {
     960           0 :     Rs = ZV_sort_shallow(ZX_Zp_liftroots(T, Rs, p, s, u));
     961           0 :     s = u;
     962             :   }
     963             :   /* Rs : mod p^s, Rrs : mod p^(r+s), vT : mod p^(r+s) */
     964          53 :   Rrs = r ? ZX_Zp_liftroots(T, Rs, p, s, r+s) : Rs;
     965          53 :   pr = powiu(p, r); prs = powiu(p, r+s); /* Usually, r=0, s=1, pr=1, pru=p */
     966          53 :   if (lgefint(prs)>3)  /* ULONG_MAX < p^(r+s), usually occur */
     967             :   {
     968         166 :     for (k = 1; k <= n_T; k++)
     969             :     {
     970         119 :       long itk = kT[k];
     971         119 :       GEN z = r? RgX_Rg_mul(gel(vT, itk), pr): gel(vT, itk);
     972         119 :       gel(vT, itk) = RgX_to_FpX(z, prs);
     973             :     }
     974             :   }
     975             :   else  /* p^(r+s) < ULONG_MAX, frequently occur */
     976             :   {
     977           6 :     ulong upr = itou(pr), uprs = itou(prs);
     978          36 :     for (k = 1; k <= n_T; k++)
     979             :     {
     980          30 :       long itk = kT[k];
     981          30 :       GEN z = r? RgX_muls(gel(vT, itk), upr): gel(vT, itk);
     982          30 :       gel(vT, itk) = RgX_to_Flx(z, uprs);
     983             :     }
     984           6 :     Rs = ZV_to_nv(Rs); Rrs = ZV_to_nv(Rrs);
     985             :   }
     986          53 :   return mkvecn(12, vT, gGH, Rs, Rrs, i_t, kt, p, pu, pr, prs,
     987             :                 mkvecsmalln(6, n, r, s, n_t, mitk, d), kTdiv);
     988             : }
     989             : 
     990             : /* Data=[vT, gGH, Rs, Rrs, i_t, kt, p, pu, pr, prs,
     991             :  *      [n, r, s, n_t, mitk, d], div] */
     992             : static GEN
     993         345 : FpX_pol_newton_general_new3(GEN Data, long k)
     994             : {
     995         345 :   GEN i_t = gel(Data, 5), p = gel(Data, 7), pu = gel(Data, 8);
     996         345 :   long d = mael(Data, 11, 6);
     997         345 :   GEN v_t_p = Fp_mk_v_t_p3(Data, k);
     998         345 :   if (DEBUGLEVEL == 3) err_printf("v_t_p=%Ps\n",v_t_p);
     999         345 :   return FpX_Newton_perm(d, v_t_p, i_t, pu, p);
    1000             : }
    1001             : 
    1002             : /* Data = [GHgen, GH, fn, p, [n, d, f, m]] */
    1003             : static GEN
    1004          46 : FpX_factcyclo_newton_general_new3(GEN Data)
    1005             : {
    1006          46 :   long i, f = mael(Data, 5, 3), m = mael(Data, 5, 4);
    1007             :   GEN Data2, pols;
    1008             :   pari_timer ti;
    1009             : 
    1010          46 :   if (m != 1) m = f;
    1011          46 :   pols = cgetg(1+m, t_VEC);
    1012          46 :   Data2 = newton_general_new_pre3(Data);
    1013             :   /* Data2 = [vT, gGH, Rs, Rrs, i_t, kt, p, pu, pr, prs,
    1014             :               [n, r, s, n_t, mitk, d], div] */
    1015          46 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
    1016         391 :   for (i = 1; i <= m; i++)
    1017         345 :     gel(pols, i) = FpX_pol_newton_general_new3(Data2, i);
    1018          46 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "FpX_pol_newton_general_new3");
    1019          46 :   return pols;
    1020             : }
    1021             : 
    1022             : /* return normalized z(-x) */
    1023             : static GEN
    1024        4483 : FpX_1st_lift_2(GEN z, GEN p)
    1025             : {
    1026        4483 :   long i, r = lg(z);
    1027        4483 :   GEN x = cgetg(r, t_POL);
    1028        4483 :   x[1] = evalsigne(1) | evalvarn(0);
    1029        4483 :   if (odd(r))
    1030        9962 :     for (i = 2; i < r; i++) gel(x,i) = odd(i)? Fp_neg(gel(z,i), p): gel(z,i);
    1031             :   else
    1032       16909 :     for (i = 2; i < r; i++) gel(x,i) = odd(i)? gel(z,i): Fp_neg(gel(z,i), p);
    1033        4483 :   return x;
    1034             : }
    1035             : 
    1036             : static GEN
    1037        3190 : FpX_1st_lift(GEN z, GEN p, ulong e, ulong el, GEN vP)
    1038             : {
    1039        3190 :    GEN z2, z3, P = gel(vP, e);
    1040        3190 :    if (!gel(vP, e)) P = gel(vP, e) = FpX_polcyclo(e, p);
    1041        3190 :    z2 = RgX_inflate(z, el);
    1042        3190 :    z3 = FpX_normalize(FpX_gcd(P, z2, p), p);
    1043        3190 :    return FpX_div(z2, z3, p);
    1044             : }
    1045             : 
    1046             : static GEN
    1047       11807 : FpX_lift(GEN z, GEN p, ulong e, ulong el, ulong r, ulong s, GEN vP)
    1048             : {
    1049       11807 :   if (s == 0)
    1050             :   {
    1051        7673 :     z = (el==2) ? FpX_1st_lift_2(z, p) : FpX_1st_lift(z, p, e, el, vP);
    1052        7673 :     if (r >= 2) z = RgX_inflate(z, upowuu(el, r-1));
    1053             :   }
    1054             :   else
    1055        4134 :     z = RgX_inflate(z, upowuu(el, r-s));
    1056       11807 :   return z;
    1057             : }
    1058             : 
    1059             : /* e is the conductor of the splitting field of p in Q(zeta_n) */
    1060             : static GEN
    1061       10501 : FpX_conductor_lift(GEN z, GEN p, GEN fn, ulong e, GEN vP)
    1062             : {
    1063       10501 :   GEN EL = gel(fn, 1), En = gel(fn, 2);
    1064       10501 :   long i, r = lg(EL), new_e = e;
    1065             : 
    1066       32077 :   for (i = 1; i < r; i++)
    1067             :   {
    1068       21576 :     long el = EL[i], en = En[i], ee = u_lval(e, el);
    1069       21576 :     if (ee < en)
    1070             :     {
    1071       11807 :       z = FpX_lift(z, p, new_e, el, en, ee, vP);
    1072       11807 :       new_e *= upowuu(el, en-ee);
    1073             :     }
    1074             :   }
    1075       10501 :   return z;
    1076             : }
    1077             : 
    1078             : /* R0 is mod p^u, d < p^u */
    1079             : static GEN
    1080        3547 : FpX_pol_newton(long j, GEN R0, GEN E, GEN D3, long d, long f, GEN p)
    1081             : {
    1082        3547 :   long i, u = D3[3];
    1083        3547 :   GEN R = cgetg(1+f, t_VEC);
    1084       15377 :   for (i = 1; i <= f; i++) gel(R, i) = gel(R0, 1+(i+j)%f);
    1085        3547 :   return FpX_Newton_perm(d, R, E, powiu(p, u), p);
    1086             : }
    1087             : 
    1088             : /* Data = [T, F, Rs, [d, nf, g, r, s, u]], nf>1 */
    1089             : static GEN
    1090        1896 : FpX_factcyclo_newton_pre(GEN Data, GEN N, GEN p, ulong m)
    1091             : {
    1092        1896 :   GEN T = gel(Data, 1), F = gel(Data, 2), Rs = gel(Data, 3), D4 = gel(Data, 4);
    1093        1896 :   long d = D4[1], nf = D4[2], g = D4[3], r = D4[4], s = D4[5], u = D4[6];
    1094             :   GEN pols, E, R, Data3;
    1095        1896 :   ulong i, n = N[1], pmodn = umodiu(p, n);
    1096             :   pari_timer ti;
    1097             : 
    1098        1896 :   if (m!=1) m = nf;
    1099        1896 :   pols = cgetg(1+m, t_VEC);
    1100        1896 :   E = set_E(pmodn, n, d, nf, g);
    1101        1896 :   R = set_R(T, F, Rs, p, nf, r, s, u);
    1102        1896 :   Data3 = mkvecsmall3(r, s, u);
    1103        1896 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
    1104        5443 :   for (i = 1; i <= m; i++) gel(pols, i) = FpX_pol_newton(i, R,E,Data3, d,nf,p);
    1105        1896 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "FpX_pol_newton");
    1106        1896 :   return pols;
    1107             : }
    1108             : 
    1109             : /* gcd(n1,n2)=1, e = gcd(d1,d2).
    1110             :  * phi(n1) = d1*nf1, phi(n2) = d2*nf2.
    1111             :  * d = lcm(d1,d2) = d1*d2/e, nf = phi(n1)*phi(n2)/d = nf1*nf2*e. */
    1112             : static GEN
    1113           7 : FpX_factcyclo_lift(long n1, GEN v1, long n2, GEN v2, GEN p, long m)
    1114             : {
    1115           7 :   pari_sp av = avma;
    1116           7 :   long d1 = degpol(gel(v1, 1)), lv1 = lg(v1), nf1 = eulerphiu(n1)/d1;
    1117           7 :   long d2 = degpol(gel(v2, 1)), lv2 = lg(v2), nf2 = eulerphiu(n2)/d2;
    1118           7 :   long e = ugcd(d1, d2), nf = nf1*nf2*e, need_factor = (e!=1);
    1119           7 :   long i, j, k = 0;
    1120           7 :   GEN z = NULL, v;
    1121             :   pari_timer ti;
    1122             : 
    1123           7 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
    1124           7 :   if (m != 1) m = nf;
    1125           7 :   v = cgetg(1+m, t_VEC);
    1126           7 :   if (m == 1)
    1127             :   {
    1128           0 :     GEN x1 = gel(v1, 1), x2 = gel(v2, 1);
    1129           0 :     z = FpX_gcd(RgX_inflate(x1, n2), RgX_inflate(x2, n1), p);
    1130           0 :     z = FpX_normalize(z, p);  /* FpX_gcd sometimes returns non-monic */
    1131           0 :     gel(v, 1) = (!need_factor)? z : gmael(FpX_factor(z, p), 1, 1);
    1132             :   }
    1133             :   else
    1134             :   {
    1135          91 :     for (i = 1; i < lv1; i++)
    1136         420 :       for (j = 1; j < lv2; j++)
    1137             :       {
    1138         336 :         GEN x1 = gel(v1, i), x2 = gel(v2, j);
    1139         336 :         z = FpX_gcd(RgX_inflate(x1, n2), RgX_inflate(x2, n1), p);
    1140         336 :         z = FpX_normalize(z, p);  /* needed */
    1141         336 :         if (!need_factor) gel(v, ++k) = z;
    1142             :         else
    1143             :         {
    1144           0 :           GEN z1 = gel(FpX_factor(z, p), 1);
    1145           0 :           long i, l = lg(z1);
    1146           0 :           for (i = 1; i < l; i++) gel(v, ++k) = gel(z1, i);
    1147             :         }
    1148             :       }
    1149             :   }
    1150           7 :   if (DEBUGLEVEL >= 6)
    1151           0 :     timer_printf(&ti, "FpX_factcyclo_lift (%ld,%ld)*(%ld,%ld)-->(%ld,%ld)-->(%ld,%ld)",
    1152           0 :         d1, nf1, d2, nf2, degpol(z), nf1*nf2, d1*d2/e, nf);
    1153           7 :   return gerepilecopy(av, v);
    1154             : }
    1155             : 
    1156             : /* n is any integer prime to p; d>1 and f>1 */
    1157             : static GEN
    1158        1349 : FpX_factcyclo_gen(GEN GH, ulong n, GEN p, long m)
    1159             : {
    1160        1349 :   GEN fn = factoru(n), fa = zm_to_ZM(fn);
    1161             :   GEN A, T, X, pd_n, v, pols;
    1162        1349 :   long i, j, pmodn = umodiu(p, n), phin = eulerphiu_fact(fn);
    1163        1349 :   long d = Fl_order(pmodn, phin, n), f = phin/d;
    1164             :   pari_timer ti;
    1165             : 
    1166        1349 :   if (m != 1) m = f;
    1167        1349 :   if (m != 1 && GH==NULL) /* FpX_factcyclo_fact is used */
    1168             :   {
    1169         381 :     GEN H = znstar_generate(n, mkvecsmall(pmodn));
    1170         381 :     GH = znstar_cosets(n, phin, H); /* representatives of G/H */
    1171             :   }
    1172             : 
    1173        1349 :   pols = cgetg(1+m, t_VEC);
    1174        1349 :   v = cgetg(1+d, t_VEC);
    1175        1349 :   pd_n = diviuexact(subis(powiu(p, d), 1), n); /* (p^d-1)/n */
    1176        1349 :   T = init_Fq(p, d, 1);
    1177        1349 :   A = pol_x(1);  /* A is a generator of F_q, q=p^d */
    1178        1349 :   if (d == 1) A = FpX_rem(A, T, p);
    1179        1349 :   random_FpX(degpol(T), varn(T), p); /* skip 1st one */
    1180        1349 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
    1181        1791 :   do X = FpXQ_pow(random_FpX(degpol(T), varn(T), p), pd_n, T, p);
    1182        1791 :   while(lg(X)<3 || equaliu(FpXQ_order(X, fa, T, p), n)==0); /* find zeta_n */
    1183        1349 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "find X");
    1184             : 
    1185        1349 :   if (m == 1)
    1186             :   {
    1187        2854 :     for (j = 1; j <= d; j++)
    1188             :     {
    1189        2183 :       gel(v, j) = pol_x(0);
    1190        2183 :       gmael(v, j, 2) = FpX_neg(X, p);
    1191        2183 :       if (j < d) X = FpXQ_pow(X, p, T, p);
    1192             :     }
    1193         671 :     gel(pols, 1) = ZXX_evalx0(FpXQXV_prod(v, T, p));
    1194             :   }
    1195             :   else
    1196             :   {
    1197             :     GEN vX, vp;
    1198         678 :     if (DEBUGLEVEL >= 6) timer_start(&ti);
    1199         678 :     vX = FpXQ_powers(X, n, T, p)+1;
    1200         678 :     vp = Fl_powers(pmodn, d, n);
    1201        7884 :     for (i = 1; i <= m; i++)
    1202             :     {
    1203       68694 :       for (j = 1; j <= d; j++)
    1204             :       {
    1205       61488 :         gel(v, j) = pol_x(0);
    1206       61488 :         gmael(v, j, 2) = FpX_neg(gel(vX, Fl_mul(GH[i], vp[j], n)), p);
    1207             :       }
    1208        7206 :       gel(pols, i) = ZXX_evalx0(FpXQXV_prod(v, T, p));
    1209             :     }
    1210         678 :     if (DEBUGLEVEL >= 6) timer_printf(&ti, "FpXQXV_prod");
    1211             :   }
    1212        1349 :   return pols;
    1213             : }
    1214             : 
    1215             : static GEN Flx_factcyclo_newton_pre(GEN Data, GEN N, ulong p, ulong m);
    1216             : /* n is an odd prime power prime to p and equal to the conductor
    1217             :  * of the splitting field of p in Q(zeta_n). d>1 and nf>1 */
    1218             : static GEN
    1219       33460 : FpX_factcyclo_newton_power(GEN N, GEN p, ulong m, int small)
    1220             : {
    1221             :   GEN Rs, H, T, F, pr, prs, pu, Data;
    1222       33460 :   long n = N[1], el = N[2], phin = N[4], g = N[5];
    1223       33460 :   long pmodn = umodiu(p, n), pmodel = umodiu(p, el);
    1224       33460 :   long d = Fl_order(pmodel, el-1, el), nf = phin/d;
    1225       33460 :   long r, s = 0, u = 1;
    1226             : 
    1227       33460 :   if (m != 1) m = nf;
    1228       35700 :   for (pu = p; cmpiu(pu,d) <= 0; u++) pu = mulii(pu,p);  /* d<p^u, pu=p^u */
    1229       33460 :   H = Fl_powers(pmodn, d, n);
    1230       33460 :   T = galoissubcyclo(utoi(n), utoi(pmodn), 0, 0);
    1231       33460 :   F = gausspol(T, H, N, d, nf, g);
    1232       33460 :   r = QX_den_pval(F, p);
    1233       33460 :   Rs = ZpX_roots_all(T, p, nf, &s);
    1234       33460 :   if (DEBUGLEVEL >= 2) err_printf("(u,s,r)=(%ld,%ld,%ld)\n",u,s,r);
    1235       33460 :   pr = powiu(p, r); prs = powiu(p, r+s); /* Usually, r=0, s=1, pr=1, prs=p */
    1236       33460 :   F = r ? RgX_to_FpX(RgX_Rg_mul(F, pr), prs) : RgX_to_FpX(F, prs);
    1237       33460 :   Data = mkvec4(T, F, Rs, mkvecsmalln(6, d, nf, g, r, s, u));
    1238       33460 :   if (small && lgefint(pu) == 3)
    1239       31564 :     return Flx_factcyclo_newton_pre(Data, N, uel(p,2), m);
    1240             :   else
    1241        1896 :     return FpX_factcyclo_newton_pre(Data, N, p, m);
    1242             : }
    1243             : 
    1244             : static GEN
    1245        2802 : FpX_split(ulong n, GEN p, ulong m)
    1246             : {
    1247             :   ulong i, j;
    1248        2802 :   GEN v, C, vx, z = rootsof1u_Fp(n, p);
    1249        2802 :   if (m == 1) return mkvec(deg1pol_shallow(gen_1, Fp_neg(z,p), 0));
    1250        1401 :   v = cgetg(m+1, t_VEC); C = coprimes_zv(n); vx = Fp_powers(z, n-1, p);
    1251       18495 :   for (i = j = 1; i <= n; i++)
    1252       17094 :     if (C[i]) gel(v, j++) = deg1pol_shallow(gen_1, Fp_neg(gel(vx,i+1), p), 0);
    1253        1401 :   return gen_sort(v, (void*)cmpii, &gen_cmp_RgX);
    1254             : }
    1255             : 
    1256             : /* n is a prime power prime to p. n is not necessarily equal to the
    1257             :  * conductor of the splitting field of p in Q(zeta_n). */
    1258             : static GEN
    1259        2658 : FpX_factcyclo_prime_power_i(long el, long e, GEN p, long m)
    1260             : {
    1261        2658 :   GEN z = set_e0_e1(el, e, p), v;
    1262        2658 :   long n = z[1], e0 = z[2], e1 = z[3], phin = z[4], g = z[5];
    1263        2658 :   long d = z[6], f = z[7]; /* d and f for n=el^e0 */
    1264             : 
    1265        2658 :   if (f == 1) v = mkvec(FpX_polcyclo(n, p));
    1266        2658 :   else if (d == 1) v = FpX_split(n, p, (m==1)? 1: f);
    1267        2658 :   else if (el == 2) v = FpX_factcyclo_gen(NULL, n, p, m); /* d==2 in this case */
    1268        2658 :   else if (!use_newton(d,f)) v = FpX_factcyclo_gen(NULL, n, p, m);
    1269             :   else
    1270             :   {
    1271        1896 :     GEN N = mkvecsmall5(n, el, e0, phin, g);
    1272        1896 :     v = FpX_factcyclo_newton_power(N, p, m, 0);
    1273             :   }
    1274        2658 :   if (e1)
    1275             :   {
    1276           0 :     long i, l = lg(v), ele1 = upowuu(el, e1);
    1277           0 :     for (i = 1; i < l; i++) gel(v, i) = RgX_inflate(gel(v, i), ele1);
    1278             :   }
    1279        2658 :   return v;
    1280             : }
    1281             : static GEN
    1282          14 : FpX_factcyclo_prime_power(long el, long e, GEN p, long m)
    1283             : {
    1284          14 :   pari_sp av = avma;
    1285          14 :   return gerepilecopy(av, FpX_factcyclo_prime_power_i(el, e, p, m));
    1286             : }
    1287             : 
    1288             : static GEN
    1289           7 : FpX_factcyclo_fact(GEN fn, GEN p, ulong m, long ascent)
    1290             : {
    1291           7 :   GEN EL = gel(fn, 1), E = gel(fn, 2), v1 = NULL;
    1292           7 :   long i, l = lg(EL), n1 = 1;
    1293             : 
    1294          21 :   for (i = 1; i < l; i++)
    1295             :   {
    1296          14 :     long j = ascent? i: l-i, n2 = upowuu(EL[j], E[j]);
    1297          14 :     GEN v2 = FpX_factcyclo_prime_power(EL[j], E[j], p, m);
    1298          14 :     v1 = v1? FpX_factcyclo_lift(n1, v1, n2, v2, p, m): v2;
    1299          14 :     n1 *= n2;
    1300             :   }
    1301           7 :   return v1;
    1302             : }
    1303             : 
    1304             : static GEN
    1305       15807 : Flv_FlvV_factorback(GEN g, GEN x, ulong q)
    1306       34270 : { pari_APPLY_ulong(Flv_factorback(g, gel(x,i), q)) }
    1307             : 
    1308             : /* return the structure and the generators of G/H. G=(Z/nZ)^, H=<p>.
    1309             :  * For efficiency assume p mod n != 1 (trivial otherwise) */
    1310             : static GEN
    1311       15807 : get_GH_gen(long n, long pmodn)
    1312             : {
    1313       15807 :   GEN G = znstar0(utoipos(n), 1), cycG = znstar_get_cyc(G);
    1314             :   GEN cycGH, gG, gGH, Ui, P;
    1315             :   ulong expG;
    1316       15807 :   P = hnfmodid(mkmat(Zideallog(G, utoi(pmodn))), cycG);
    1317       15807 :   cycGH = ZV_to_nv(ZM_snf_group(P, NULL, &Ui));
    1318       15807 :   expG = itou(cyc_get_expo(cycG));
    1319       15807 :   gG = ZV_to_Flv(znstar_get_gen(G), n);
    1320       15807 :   gGH = Flv_FlvV_factorback(gG, ZM_to_Flm(Ui, expG), n);
    1321       15807 :   return mkvec2(gGH, cycGH);
    1322             : }
    1323             : 
    1324             : /* 1st output */
    1325             : static void
    1326           0 : header(GEN fn, long n, long d, long f, GEN p)
    1327             : {
    1328           0 :   GEN EL = gel(fn, 1), E = gel(fn, 2);
    1329           0 :   long i, l = lg(EL)-1;
    1330           0 :   err_printf("n=%lu=", n);
    1331           0 :   for (i = 1; i <= l; i++)
    1332             :   {
    1333           0 :     long el = EL[i], e = E[i];
    1334           0 :     err_printf("%ld", el);
    1335           0 :     if (e > 1) err_printf("^%ld", e);
    1336           0 :     if (i < l) err_printf("*");
    1337             :   }
    1338           0 :   err_printf(", p=%Ps, phi(%lu)=%lu*%lu\n", p, n, d, f);
    1339           0 :   err_printf("(n,d,f) : (%ld,%ld,%ld) --> ",n,d,f);
    1340           0 : }
    1341             : 
    1342             : static ulong
    1343       94605 : FpX_factcyclo_just_conductor_init(GEN *pData, ulong n, GEN p, ulong m)
    1344             : {
    1345       94605 :   GEN fn = factoru(n), GH = NULL, GHgen = NULL;
    1346       94605 :   long phin = eulerphiu_fact(fn), pmodn = umodiu(p, n);
    1347       94605 :   long d = Fl_order(pmodn, phin, n), f = phin/d; /* d > 1 */
    1348       94605 :   ulong action = set_action(fn, p, d, f);
    1349       94605 :   if (action & ~NEWTON_POWER)
    1350             :   { /* needed for GENERAL* */
    1351       60390 :     GEN H = znstar_generate(n, mkvecsmall(pmodn));
    1352       60390 :     GH = znstar_cosets(n, phin, H); /* representatives of G/H */
    1353       60390 :     if (action & (NEWTON_GENERAL_NEW | NEWTON_GENERAL))
    1354       15807 :       GHgen = get_GH_gen(n, pmodn);  /* gen and order of G/H */
    1355             :   }
    1356       94605 :   *pData = mkvec5(GHgen, GH, fn, p, mkvecsmall4(n, d, f, m));
    1357       94605 :   if (DEBUGLEVEL >= 1)
    1358             :   {
    1359           0 :     err_printf("(%ld,%ld,%ld)  action=%ld\n", n, d, f, action);
    1360           0 :     if (GHgen)
    1361             :     {
    1362           0 :       GEN cycGH = gel(GHgen,2), gGH = gel(GHgen,1);
    1363           0 :       err_printf("G(K/Q)=%Ps gen=%Ps\n", zv_to_ZV(cycGH), zv_to_ZV(gGH));
    1364             :     }
    1365             :   }
    1366       94605 :   return action;
    1367             : }
    1368             : 
    1369             : static GEN
    1370        5698 : FpX_factcyclo_just_conductor(ulong n, GEN p, ulong m)
    1371             : {
    1372             :   GEN Data, fn;
    1373        5698 :   ulong action = FpX_factcyclo_just_conductor_init(&Data, n, p, m);
    1374        5698 :   fn = gel(Data,3);
    1375        5698 :   if (action & GENERAL)
    1376         587 :     return FpX_factcyclo_gen(gel(Data,2), n, p, m);
    1377        5111 :   else if (action & NEWTON_POWER)
    1378        2644 :     return FpX_factcyclo_prime_power_i(ucoeff(fn,1,1), ucoeff(fn,1,2), p, m);
    1379        2467 :   else if (action & NEWTON_GENERAL)
    1380        2414 :     return FpX_factcyclo_newton_general(Data);
    1381          53 :   else if (action & NEWTON_GENERAL_NEW)
    1382          46 :     return FpX_factcyclo_newton_general_new3(Data);
    1383             :   else
    1384           7 :     return FpX_factcyclo_fact(fn, p, m, action & ASCENT);
    1385             : }
    1386             : 
    1387             : static GEN
    1388       10656 : FpX_factcyclo_i(ulong n, GEN p, long fl)
    1389             : {
    1390       10656 :   GEN fn = factoru(n), z;
    1391       10656 :   long phin = eulerphiu_fact(fn), pmodn = umodiu(p, n);
    1392       10656 :   ulong d = Fl_order(pmodn, phin, n), f = phin/d, fK;
    1393             : 
    1394       10656 :   if (DEBUGLEVEL >= 1) header(fn, n, d, f, p);
    1395       10656 :   if (f == 1) { z = FpX_polcyclo(n, p); return fl? z: mkvec(z); }
    1396        8500 :   else if (d == 1) /* p=1 (mod n), zeta_n in Z_p */
    1397         774 :   { z = FpX_split(n, p, fl? 1: f); return fl? gel(z,1): z; }
    1398        7726 :   fK = znstar_conductor(znstar_generate(n, mkvecsmall(pmodn)));
    1399        7726 :   if (fK != n && umodiu(p, fK) == 1)
    1400        2028 :     z = FpX_split(fK, p, fl? 1: f);
    1401             :   else
    1402        5698 :     z = FpX_factcyclo_just_conductor(fK, p, fl? 1: f);
    1403        7726 :   if (n > fK)
    1404             :   {
    1405        4294 :     GEN vP = const_vec(n, NULL);
    1406        4294 :     long i, l = fl? 2: lg(z);
    1407       14795 :     for (i = 1; i < l; i++)
    1408       10501 :       gel(z, i) = FpX_conductor_lift(gel(z, i), p, fn, fK, vP);
    1409             :   }
    1410        7726 :   return fl? gel(z,1): gen_sort(z,(void*)cmpii, &gen_cmp_RgX);
    1411             : }
    1412             : 
    1413             : GEN
    1414          42 : FpX_factcyclo(ulong n, GEN p, ulong m)
    1415          42 : { pari_sp av = avma; return gerepilecopy(av, FpX_factcyclo_i(n, p, m)); }
    1416             : 
    1417             : /*  Data = [H, GH, i_t, d0, kT, [n, d, f, n_T, mitk]]
    1418             :  *  N2 = [p, pr, pu, pru] */
    1419             : static GEN
    1420       24102 : Flx_pol_newton_general(GEN Data, GEN N2, GEN vT, ulong x)
    1421             : {
    1422       24102 :   GEN i_t = gel(Data, 3), kT = gel(Data, 5), N = gel(Data, 6);
    1423       24102 :   long k, d = N[2], n_T = N[4], mitk = N[5];
    1424       24102 :   long p = N2[1], pr = N2[2], pu = N2[3], pru = N2[4];
    1425       24102 :   GEN R = cgetg(1+mitk, t_VECSMALL);
    1426             : 
    1427      123717 :   for (k = 1; k <= n_T; k++) uel(R,kT[k]) = Flx_eval(gel(vT, kT[k]), x, pru) / pr;
    1428       24102 :   return Flx_Newton_perm(d, R, i_t, pu, p);
    1429             : }
    1430             : 
    1431             : /* n is any integer prime to p, but must be equal to the conductor
    1432             :  * of the splitting field of p in Q(zeta_n).
    1433             :  * GH=G/H={g_1, ... ,g_f}
    1434             :  * Data = [GHgen, GH, fn, p, [n, d, f, m]] */
    1435             : static GEN
    1436       13340 : Flx_factcyclo_newton_general(GEN Data)
    1437             : {
    1438       13340 :   GEN GH = gel(Data, 2), fn = gel(Data, 3), p = gel(Data, 4);
    1439       13340 :   ulong up = p[2], n = mael(Data, 5, 1), pmodn = up%n;
    1440       13340 :   long d = mael(Data, 5, 2), f = mael(Data, 5, 3), m = mael(Data, 5, 4);
    1441       13340 :   long i, k, n_T, mitk, r, s = 0, u = 1;
    1442             :   GEN vT, kT, H, i_t, T, d0d1, Data2, Data3, R, Rs, pr, pu, pru;
    1443             :   pari_timer ti;
    1444             : 
    1445       13340 :   if (m != 1) m = f;
    1446       14593 :   for (pu = p; cmpiu(pu,d) <= 0; u++) pu = muliu(pu, up);  /* d<pu, pu=p^u */
    1447             : 
    1448       13340 :   H = Fl_powers(pmodn, d-1, n); /* H=<p> */
    1449       13340 :   i_t = get_i_t(n, pmodn); /* i_t[1+i]=k ==> t_i=t_k */
    1450       13340 :   kT = get_kT(i_t, d); n_T = lg(kT)-1; mitk = kT[n_T];
    1451       13340 :   if (DEBUGLEVEL == 2) err_printf("kT=%Ps %ld elements\n",kT,n_T);
    1452       13340 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
    1453       13340 :   T = galoissubcyclo(utoi(n), utoi(pmodn), 0, 0);
    1454       13340 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "galoissubcyclo");
    1455       13340 :   d0d1 = get_d0_d1(T, gel(fn,1)); /* d0*T_k(x) is in Z[x] */
    1456       13340 :   Data2 = mkvecn(6, H, GH, i_t, d0d1, kT, mkvecsmalln(5, n, d, f, n_T, mitk));
    1457       13340 :   vT = get_vT(Data2, 0);
    1458       13340 :   if (DEBUGLEVEL == 4) err_printf("vT=%Ps\n",vT);
    1459       13340 :   r = QXV_den_pval(vT, kT, p);
    1460       13340 :   Rs = ZpX_roots_all(T, p, f, &s);
    1461       13340 :   if (DEBUGLEVEL >= 2) err_printf("(u,s,r)=(%ld,%ld,%ld)\n",u,s,r);
    1462       13340 :   if (r+u < s) pari_err_BUG("Flx_factcyclo_newton_general, T is not separable mod p^(r+u)");
    1463             :   /* R and vT are mod p^(r+u) */
    1464       13340 :   R = (r+u==s) ? Rs : ZV_sort_shallow(ZX_Zp_liftroots(T, Rs, p, s, r+u));
    1465       13340 :   pr = powiu(p, r); pru = powiu(p, r+u); /* Usually, r=0, s=1, pr=1, pru=p */
    1466       13340 :   if (lgefint(pru) > 3)  /* ULONG_MAX < p^(r+u), probably won't occur */
    1467             :   {
    1468           0 :     for (k = 1; k <= n_T; k++)
    1469             :     {
    1470           0 :       long itk = kT[k];
    1471           0 :       GEN z = r? RgX_Rg_mul(gel(vT, itk), pr): gel(vT, itk);
    1472           0 :       gel(vT, itk) = RgX_to_FpX(z, pru);
    1473             :     }
    1474           0 :     Data3 = mkvec4(p, pr, pu, pru);
    1475           0 :     for (i = 1; i <= m; i++)
    1476           0 :       gel(R,i) = ZX_to_nx(FpX_pol_newton_general(Data2, Data3, vT, gel(R,i)));
    1477           0 :     if (DEBUGLEVEL >= 6) timer_printf(&ti, "FpX_pol_newton_general");
    1478             :   }
    1479             :   else
    1480             :   {
    1481       13340 :     ulong upr = itou(pr), upru = itou(pru), upu = itou(pu);
    1482       63651 :     for (k = 1; k <= n_T; k++)
    1483             :     {
    1484       50311 :       long itk = kT[k];
    1485       50311 :       GEN z = r? RgX_muls(gel(vT, itk), upr): gel(vT, itk);
    1486       50311 :       gel(vT, itk) = RgX_to_Flx(z, upru);
    1487             :     }
    1488       13340 :     Data3 = mkvecsmall4(up, upr, upu, upru);
    1489       13340 :     if (DEBUGLEVEL >= 6) timer_start(&ti);
    1490       37442 :     for (i = 1; i <= m; i++)
    1491       24102 :       gel(R,i) = Flx_pol_newton_general(Data2, Data3, vT, itou(gel(R,i)));
    1492       13340 :     if (DEBUGLEVEL >= 6) timer_printf(&ti, "Flx_pol_newton_general");
    1493             :   }
    1494       13340 :   return R;
    1495             : }
    1496             : 
    1497             : /*  Data=[vT, gGH, Rs, Rrs, i_t, kt, p, pu, pr, prs,
    1498             :  *       [n, r, s, n_t, mitk, d], div] */
    1499             : static GEN
    1500         336 : Flx_pol_newton_general_new3(GEN Data, long k)
    1501             : {
    1502         336 :   GEN i_t = gel(Data,5), p = gel(Data,7), pu = gel(Data,8), prs = gel(Data,10);
    1503         336 :   long d = mael(Data, 11, 6);
    1504          48 :   GEN v_t_p = (lgefint(prs)>3)? ZV_to_nv(Fp_mk_v_t_p3(Data, k))
    1505         336 :                               : Fl_mk_v_t_p3(Data, k);
    1506         336 :   if (DEBUGLEVEL == 3) err_printf("v_t_p=%Ps\n",v_t_p);
    1507         336 :   return Flx_Newton_perm(d, v_t_p, i_t, pu[2], p[2]);
    1508             : }
    1509             : 
    1510             : /* Data = [GHgen, GH, fn, p, [n, d, f, m]] */
    1511             : static GEN
    1512           7 : Flx_factcyclo_newton_general_new3(GEN Data)
    1513             : {
    1514           7 :   long i, f = mael(Data, 5, 3), m = mael(Data, 5, 4);
    1515             :   GEN Data2, pu, pols;
    1516             :   pari_timer ti;
    1517             : 
    1518           7 :   if (m != 1) m = f;
    1519           7 :   pols = cgetg(1+m, t_VEC);
    1520           7 :   Data2 = newton_general_new_pre3(Data); pu = gel(Data2, 8);
    1521             :   /* Data2 = [vT, gGH, Rs, Rrs, i_t, kt, p, pu, pr, prs,
    1522             :               [n, r, s, n_t, mitk, d], div] */
    1523           7 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
    1524           7 :   if (lgefint(pu) > 3)  /* ULONG_MAX < p^u, probably won't occur */
    1525           0 :   { for (i = 1; i <= m; i++)
    1526           0 :       gel(pols, i) = ZX_to_nx(FpX_pol_newton_general_new3(Data2, i));
    1527           0 :     if (DEBUGLEVEL >= 6) timer_printf(&ti, "FpX_pol_newton_general_new3"); }
    1528             :   else
    1529         343 :   { for (i = 1; i <= m; i++)
    1530         336 :       gel(pols, i) = Flx_pol_newton_general_new3(Data2, i);
    1531           7 :     if (DEBUGLEVEL >= 6) timer_printf(&ti, "Flx_pol_newton_general_new3"); }
    1532           7 :   return pols;
    1533             : }
    1534             : 
    1535             : /* return normalized z(-x) */
    1536             : static GEN
    1537       65972 : Flx_1st_lift_2(GEN z, ulong p)
    1538             : {
    1539       65972 :   long i, r = lg(z);
    1540       65972 :   GEN x = cgetg(r, t_VECSMALL);
    1541       65972 :   if (odd(r))
    1542      194270 :     for (i = 2; i<r; i++) uel(x,i) = odd(i)? Fl_neg(uel(z,i), p) : uel(z,i);
    1543             :   else
    1544      231633 :     for (i = 2; i<r; i++) uel(x,i) = odd(i)? uel(z,i): Fl_neg(uel(z,i), p);
    1545       65972 :   return x;
    1546             : }
    1547             : 
    1548             : /* el does not divides e.
    1549             :  * construct Phi_{e*el}(x) from Phi_e(x). */
    1550             : static GEN
    1551       42709 : Flx_1st_lift(GEN z, ulong p, ulong e, ulong el, GEN vP)
    1552             : {
    1553       42709 :    GEN z2, z3, P = gel(vP, e);
    1554       42709 :    if (!gel(vP, e)) P = gel(vP, e) = Flx_polcyclo(e, p);
    1555       42709 :    z2 = Flx_inflate(z, el);
    1556       42709 :    z3 = Flx_normalize(Flx_gcd(P, z2, p), p);
    1557       42709 :    return Flx_div(z2, z3, p);
    1558             : }
    1559             : 
    1560             : static GEN
    1561      158258 : Flx_lift(GEN z, ulong p, ulong e, ulong el, ulong r, ulong s, GEN vP)
    1562             : {
    1563      158258 :   if (s == 0)
    1564             :   {
    1565      108681 :     z = (el==2) ? Flx_1st_lift_2(z, p) : Flx_1st_lift(z, p, e, el, vP);
    1566      108681 :     if (r >= 2) z = Flx_inflate(z, upowuu(el, r-1));
    1567             :   }
    1568             :   else
    1569       49577 :     z = Flx_inflate(z, upowuu(el, r-s));
    1570      158258 :   return z;
    1571             : }
    1572             : 
    1573             : /* e is the conductor of the splitting field of p in Q(zeta_n) */
    1574             : static GEN
    1575      140965 : Flx_conductor_lift(GEN z, ulong p, GEN fn, ulong e, GEN vP)
    1576             : {
    1577      140965 :   GEN EL = gel(fn, 1), En = gel(fn, 2);
    1578      140965 :   long i, r = lg(EL), new_e = e;
    1579             : 
    1580      432170 :   for (i = 1; i < r; i++)
    1581             :   {
    1582      291205 :     long el = EL[i], en = En[i], ee = u_lval(e, el);
    1583      291205 :     if (ee < en)
    1584             :     {
    1585      158258 :       z = Flx_lift(z, p, new_e, el, en, ee, vP);
    1586      158258 :       new_e *= upowuu(el, en-ee);
    1587             :     }
    1588             :   }
    1589      140965 :   return z;
    1590             : }
    1591             : 
    1592             : /* R0 is mod p^u, d < p^u */
    1593             : static GEN
    1594       57409 : Flx_pol_newton(long j, GEN R0, GEN E, GEN D3, long d, long f, ulong p)
    1595             : {
    1596       57409 :   ulong u = D3[3];
    1597       57409 :   GEN R = cgetg(f+1, t_VECSMALL);
    1598             :   long i;
    1599      233473 :   for (i = 1; i <= f; i++) R[i] = R0[1+(i+j)%f];
    1600       57409 :   return Flx_Newton_perm(d, R, E, upowuu(p,u), p);
    1601             : }
    1602             : 
    1603             : /* Data = [T, F, Rs, [d, nf, g, r, s, u]], nf>1 */
    1604             : static GEN
    1605       31564 : Flx_factcyclo_newton_pre(GEN Data, GEN N, ulong p, ulong m)
    1606             : {
    1607       31564 :   GEN T = gel(Data, 1), F = gel(Data, 2), Rs = gel(Data, 3), D4 = gel(Data, 4);
    1608       31564 :   long d = D4[1], nf = D4[2], g = D4[3], r = D4[4], s = D4[5], u = D4[6];
    1609       31564 :   GEN pols, E, R, p0 = utoi(p), Data3;
    1610       31564 :   ulong i, n = N[1], pmodn = p%n;
    1611             :   pari_timer ti;
    1612             : 
    1613       31564 :   if (m != 1) m = nf;
    1614       31564 :   pols = cgetg(1+m, t_VEC);
    1615       31564 :   E = set_E(pmodn, n, d, nf, g);
    1616       31564 :   R = set_R(T, F, Rs, p0, nf, r, s, u);
    1617       31564 :   R = ZV_to_nv(R);
    1618       31564 :   Data3 = mkvecsmall3(r, s, u);
    1619       31564 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
    1620       88973 :   for (i = 1; i <= m; i++) gel(pols, i) = Flx_pol_newton(i, R,E,Data3, d,nf,p);
    1621       31564 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "Flx_pol_newton");
    1622       31564 :   return pols;
    1623             : }
    1624             : 
    1625             : /* gcd(n1,n2)=1, e = gcd(d1,d2).
    1626             :  * phi(n1) = d1*nf1, phi(n2) = d2*nf2.
    1627             :  * d = lcm(d1,d2) = d1*d2/e, nf = phi(n1)*phi(n2)/d = nf1*nf2*e. */
    1628             : static GEN
    1629           0 : Flx_factcyclo_lift(long n1, GEN v1, long n2, GEN v2, long p, long m)
    1630             : {
    1631           0 :   pari_sp av = avma;
    1632           0 :   long d1 = degpol(gel(v1, 1)), lv1 = lg(v1), nf1 = eulerphiu(n1)/d1;
    1633           0 :   long d2 = degpol(gel(v2, 1)), lv2 = lg(v2), nf2 = eulerphiu(n2)/d2;
    1634           0 :   long e = ugcd(d1, d2), nf = nf1*nf2*e, need_factor = (e!=1);
    1635           0 :   long i, j, k = 0;
    1636           0 :   GEN v, z = NULL;
    1637             :   pari_timer ti;
    1638             : 
    1639           0 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
    1640           0 :   if (m != 1) m = nf;
    1641           0 :   v = cgetg(1+m, t_VEC);
    1642           0 :   if (m == 1)
    1643             :   {
    1644           0 :     GEN x1 = gel(v1, 1), x2 = gel(v2, 1);
    1645           0 :     z = Flx_gcd(Flx_inflate(x1, n2), Flx_inflate(x2, n1), p);
    1646           0 :     z = Flx_normalize(z, p);  /* Flx_gcd sometimes returns non-monic */
    1647           0 :     gel(v, 1) = (!need_factor)? z : gmael(Flx_factor(z, p), 1, 1);
    1648             :   }
    1649             :   else
    1650           0 :     for (i = 1; i < lv1; i++)
    1651           0 :       for (j = 1; j < lv2; j++)
    1652             :       {
    1653           0 :         GEN x1 = gel(v1, i), x2 = gel(v2, j);
    1654           0 :         z = Flx_gcd(Flx_inflate(x1, n2), Flx_inflate(x2, n1), p);
    1655           0 :         z = Flx_normalize(z, p);  /* needed */
    1656           0 :         if (!need_factor) gel(v, ++k) = z;
    1657             :         else
    1658             :         {
    1659           0 :           GEN z1 = gel(Flx_factor(z, p), 1);
    1660           0 :           long i, l = lg(z1);
    1661           0 :           for (i = 1; i < l; i++) gel(v, ++k) = gel(z1, i);
    1662             :         }
    1663             :       }
    1664           0 :   if (DEBUGLEVEL >= 6)
    1665           0 :     timer_printf(&ti, "Flx_factcyclo_lift (%ld,%ld)*(%ld,%ld)-->(%ld,%ld)-->(%ld,%ld)",
    1666           0 :         d1, nf1, d2, nf2, degpol(z), nf1*nf2, d1*d2/e, nf);
    1667           0 :   return gerepilecopy(av, v);
    1668             : }
    1669             : 
    1670             : /* factor polcyclo(n) mod p based on an idea of Bill Allombert; d>1 and nf>1 */
    1671             : static GEN
    1672       43996 : Flx_factcyclo_gen(GEN GH, ulong n, ulong p, ulong m)
    1673             : {
    1674       43996 :   GEN fu = factoru(n), fa = zm_to_ZM(fu), p0 = utoi(p);
    1675             :   GEN T, X, pd_n, v, pols;
    1676       43996 :   ulong i, j, pmodn = p%n, phin = eulerphiu_fact(fu);
    1677       43996 :   ulong d = Fl_order(pmodn, phin, n), nf = phin/d;
    1678             :   pari_timer ti;
    1679             : 
    1680       43996 :   if (m != 1) m = nf;
    1681       43996 :   if (m != 1 && !GH) /* Flx_factcyclo_fact is used */
    1682             :   {
    1683           0 :     GEN H = znstar_generate(n, mkvecsmall(pmodn));
    1684           0 :     GH = znstar_cosets(n, phin, H); /* representatives of G/H */
    1685             :   }
    1686             : 
    1687       43996 :   pols = cgetg(1+m, t_VEC);
    1688       43996 :   v = cgetg(1+d, t_VEC);
    1689       43996 :   pd_n = diviuexact(subis(powiu(p0, d), 1), n); /* (p^d-1)/n */
    1690       43996 :   T = ZX_to_Flx(init_Fq(p0, d, evalvarn(1)), p);
    1691       43996 :   random_Flx(degpol(T), T[1], p);  /* skip 1st one */
    1692       43996 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
    1693       82769 :   do X = Flxq_pow(random_Flx(degpol(T), T[1], p), pd_n, T, p);
    1694       82769 :   while (lg(X)<3 || equaliu(Flxq_order(X, fa, T, p), n)==0); /* find zeta_n */
    1695       43996 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "find X");
    1696             : 
    1697       43996 :   if (m == 1)
    1698             :   {
    1699      102528 :     for (j = 1; j <= d; j++)
    1700             :     {
    1701       80481 :       gel(v, j) = polx_FlxX(0, 1);
    1702       80481 :       gmael(v, j, 2) = Flx_neg(X, p);
    1703       80481 :       if (j < d) X = Flxq_powu(X, p, T, p);
    1704             :     }
    1705       22047 :     gel(pols, 1) = FlxX_to_Flx(FlxqXV_prod(v, T, p));
    1706             :   }
    1707             :   else
    1708             :   {
    1709             :     GEN vX, vp;
    1710       21949 :     if (DEBUGLEVEL >= 6) timer_start(&ti);
    1711       21949 :     vX = Flxq_powers(X, n, T, p)+1;
    1712       21949 :     vp = Fl_powers(pmodn, d, n);
    1713      190844 :     for (i = 1; i <= m; i++)
    1714             :     {
    1715      757867 :       for (j = 1; j <= d; j++)
    1716             :       {
    1717      588972 :         gel(v, j) = polx_FlxX(0, 1);
    1718      588972 :         gmael(v, j, 2) = Flx_neg(gel(vX, Fl_mul(GH[i], vp[j], n)), p);
    1719             :       }
    1720      168895 :       gel(pols, i) = FlxX_to_Flx(FlxqXV_prod(v, T, p));
    1721             :     }
    1722       21949 :     if (DEBUGLEVEL >= 6) timer_printf(&ti, "FlxqXV_prod");
    1723             :   }
    1724       43996 :   return pols;
    1725             : }
    1726             : 
    1727             : static int
    1728     1734612 : cmpGuGu(GEN a, GEN b) { return (ulong)a < (ulong)b? -1: (a == b? 0: 1); }
    1729             : 
    1730             : /* p=1 (mod n). If m!=1, then m=phi(n) */
    1731             : static GEN
    1732       34186 : Flx_split(ulong n, ulong p, ulong m)
    1733             : {
    1734       34186 :   ulong i, j, z = rootsof1_Fl(n, p);
    1735             :   GEN v, C, vx;
    1736       34186 :   if (m == 1) return mkvec(mkvecsmall3(0, Fl_neg(z,p), 1));
    1737       17016 :   v = cgetg(m+1, t_VEC); C = coprimes_zv(n); vx = Fl_powers(z, n-1, p);
    1738      223026 :   for (i = j = 1; i <= n; i++)
    1739      206010 :     if (C[i]) gel(v, j++) = mkvecsmall3(0, Fl_neg(vx[i+1], p), 1);
    1740       17016 :   return gen_sort(v, (void*)cmpGuGu, &gen_cmp_RgX);
    1741             : }
    1742             : 
    1743             : /* d==1 or f==1 occurs */
    1744             : static GEN
    1745       31564 : Flx_factcyclo_prime_power_i(long el, long e, long p, long m)
    1746             : {
    1747       31564 :   GEN p0 = utoipos(p), z = set_e0_e1(el, e, p0), v;
    1748       31564 :   long n = z[1], e0 = z[2], e1 = z[3], phin = z[4], g = z[5];
    1749       31564 :   long d = z[6], f = z[7]; /* d and f for n=el^e0 */
    1750             : 
    1751       31564 :   if (f == 1) v = mkvec(Flx_polcyclo(n, p));
    1752       31564 :   else if (d == 1) v = Flx_split(n, p, (m==1)?1:f);
    1753       31564 :   else if (el == 2) v = Flx_factcyclo_gen(NULL, n, p, m);/* d==2 in this case */
    1754       31564 :   else if (!use_newton(d, f)) v = Flx_factcyclo_gen(NULL, n, p, m);
    1755             :   else
    1756             :   {
    1757       31564 :     GEN N = mkvecsmall5(n, el, e0, phin, g);
    1758       31564 :     v = FpX_factcyclo_newton_power(N, p0, m, 1);
    1759       31564 :     if (typ(gel(v,1)) == t_POL)
    1760             :     { /* ZX: convert to Flx */
    1761           0 :       long i, l = lg(v);
    1762           0 :       for (i = 1; i < l; i++) gel(v,i) = ZX_to_nx(gel(v,i));
    1763             :     }
    1764             :   }
    1765       31564 :   if (e1)
    1766             :   {
    1767           0 :     long i, l = lg(v), ele1 = upowuu(el, e1);
    1768           0 :     for (i = 1; i < l; i++) gel(v, i) = Flx_inflate(gel(v, i), ele1);
    1769             :   }
    1770       31564 :   return v;
    1771             : }
    1772             : static GEN
    1773           0 : Flx_factcyclo_prime_power(long el, long e, long p, long m)
    1774             : {
    1775           0 :   pari_sp av = avma;
    1776           0 :   return gerepilecopy(av, Flx_factcyclo_prime_power_i(el, e, p, m));
    1777             : }
    1778             : 
    1779             : static GEN
    1780           0 : Flx_factcyclo_fact(GEN fn, ulong p, ulong m, long ascent)
    1781             : {
    1782           0 :   GEN EL = gel(fn, 1), E = gel(fn, 2), v1, v2;
    1783           0 :   long l = lg(EL), i, j, n1, n2;
    1784             : 
    1785           0 :   j = ascent? 1: l-1;
    1786           0 :   n1 = upowuu(EL[j], E[j]);
    1787           0 :   v1 = Flx_factcyclo_prime_power(EL[j], E[j], p, m);
    1788           0 :   for (i = 2; i < l; i++)
    1789             :   {
    1790           0 :     j = ascent? i: l-i;
    1791           0 :     n2 = upowuu(EL[j], E[j]);
    1792           0 :     v2 = Flx_factcyclo_prime_power(EL[j], E[j], p, m);
    1793           0 :     v1 = Flx_factcyclo_lift(n1, v1, n2, v2, p, m);
    1794           0 :     n1 *= n2;
    1795             :   }
    1796           0 :   return v1;
    1797             : }
    1798             : 
    1799             : static GEN
    1800       88907 : Flx_factcyclo_just_conductor(ulong n, ulong p, ulong m)
    1801             : {
    1802             :   GEN Data, fn;
    1803       88907 :   ulong action = FpX_factcyclo_just_conductor_init(&Data, n, utoipos(p), m);
    1804       88907 :   fn = gel(Data,3);
    1805       88907 :   if (action & GENERAL)
    1806       43996 :     return Flx_factcyclo_gen(gel(Data,2), n, p, m);
    1807       44911 :   else if (action & NEWTON_POWER)
    1808       31564 :     return Flx_factcyclo_prime_power_i(ucoeff(fn,1,1), ucoeff(fn,1,2), p, m);
    1809       13347 :   else if (action & NEWTON_GENERAL)
    1810       13340 :     return Flx_factcyclo_newton_general(Data);
    1811           7 :   else if (action & NEWTON_GENERAL_NEW)
    1812           7 :     return Flx_factcyclo_newton_general_new3(Data);
    1813             :   else
    1814           0 :     return Flx_factcyclo_fact(fn, p, m, action & ASCENT);
    1815             : }
    1816             : 
    1817             : static GEN
    1818      156721 : Flx_factcyclo_i(ulong n, ulong p, ulong fl)
    1819             : {
    1820      156721 :   GEN fn = factoru(n), z;
    1821      156721 :   ulong phin = eulerphiu_fact(fn), pmodn = p%n;
    1822      156721 :   ulong d = Fl_order(pmodn, phin, n), f = phin/d, fK;
    1823             : 
    1824      156721 :   if (DEBUGLEVEL >= 1) header(fn, n, d, f, utoi(p));
    1825      156721 :   if (f == 1) { z = Flx_polcyclo(n, p); return fl? z: mkvec(z); }
    1826      123093 :   if (d == 1) /* p=1 (mod n), zeta_n in Z_p */
    1827        9488 :   { z = Flx_split(n, p, fl? 1: f); return fl? gel(z,1): z; }
    1828      113605 :   fK = znstar_conductor(znstar_generate(n, mkvecsmall(pmodn)));
    1829      113605 :   if (fK != n && p % fK == 1)
    1830       24698 :     z = Flx_split(fK, p, fl? 1: f);
    1831             :   else
    1832       88907 :     z = Flx_factcyclo_just_conductor(fK, p, fl? 1: f);
    1833      113605 :   if (n > fK)
    1834             :   {
    1835       58055 :     GEN vP = const_vec(n, NULL);
    1836       58055 :     long i, l = fl? 2: lg(z);
    1837      199020 :     for (i = 1; i < l; i++)
    1838      140965 :       gel(z, i) = Flx_conductor_lift(gel(z, i), p, fn, fK, vP);
    1839             :   }
    1840      113605 :   return fl? gel(z,1): gen_sort(z,(void*)cmpGuGu, &gen_cmp_RgX);
    1841             : }
    1842             : 
    1843             : GEN
    1844         308 : Flx_factcyclo(ulong n, ulong p, ulong m)
    1845         308 : { pari_sp av = avma; return gerepilecopy(av, Flx_factcyclo_i(n, p, m)); }
    1846             : 
    1847             : GEN
    1848      167027 : factormodcyclo(long n, GEN p, long fl, long v)
    1849             : {
    1850      167027 :   const char *fun = "factormodcyclo";
    1851      167027 :   pari_sp av = avma;
    1852             :   long i, l;
    1853             :   GEN z;
    1854      167027 :   if (v < 0) v = 0;
    1855      167027 :   if (fl < 0 || fl > 1) pari_err_FLAG(fun);
    1856      167027 :   if (n <= 0) pari_err_DOMAIN(fun, "n", "<=", gen_0, stoi(n));
    1857      167027 :   if (typ(p) != t_INT) pari_err_TYPE(fun, p);
    1858      167027 :   if (dvdui(n, p)) pari_err_COPRIME(fun, stoi(n), p);
    1859      167027 :   if (fl)
    1860             :   {
    1861       83503 :     if (lgefint(p) == 3)
    1862       78203 :       z = Flx_to_ZX(Flx_factcyclo_i(n, p[2], 1));
    1863             :     else
    1864        5300 :       z = FpX_factcyclo_i(n, p, 1);
    1865       83503 :     setvarn(z, v);
    1866       83503 :     return gerepileupto(av, FpX_to_mod(z, p));
    1867             :   }
    1868             :   else
    1869             :   {
    1870       83524 :     if (lgefint(p) == 3)
    1871       78210 :       z = FlxC_to_ZXC(Flx_factcyclo_i(n, p[2], 0));
    1872             :     else
    1873        5314 :       z = FpX_factcyclo_i(n, p, 0);
    1874      465409 :     l = lg(z); for (i = 1; i < l; i++) setvarn(gel(z, i), v);
    1875       83524 :     return gerepileupto(av, FpXC_to_mod(z, p));
    1876             :   }
    1877             : }

Generated by: LCOV version 1.14