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