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 : }
|