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 92826 : QX_den_pval(GEN x, GEN p)
55 : {
56 92826 : long i, vmax = 0, l = lg(x);
57 334882 : for (i = 2; i < l; i++)
58 : {
59 242056 : GEN z = gel(x, i);
60 : long v;
61 242056 : if (typ(z)==t_FRAC && (v = Z_pval(gel(z, 2), p)) > vmax) vmax = v;
62 : }
63 92826 : return vmax;
64 : }
65 :
66 : static long
67 15807 : QXV_den_pval(GEN vT, GEN kT, GEN p)
68 : {
69 15807 : long k, vmax = 0, l = lg(kT);
70 75173 : for (k = 1; k < l; k++)
71 : {
72 59366 : long v = QX_den_pval(gel(vT, kT[k]), p);
73 59366 : if (v > vmax) vmax = v;
74 : }
75 15807 : 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 33460 : set_E(long pmodn, long n, long d, long f, long g)
84 : {
85 : long i, j;
86 33460 : GEN E = const_vecsmall(n-1, 0);
87 33460 : pari_sp av = avma;
88 33460 : GEN C = Fl_powers(g, f-1, n);
89 121856 : for (i = 1; i <= f; i++)
90 : {
91 88396 : ulong x = C[i];
92 1239140 : for (j = 1; j <= d; j++) { x = Fl_mul(x, pmodn, n); E[x] = i; }
93 : }
94 33460 : 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 98534 : list_el_n(ulong el0, ulong n, GEN d1, long n_el)
119 : {
120 98534 : GEN v = cgetg(n_el + 1, t_VECSMALL);
121 : forprime_t T;
122 : long i;
123 98534 : u_forprime_arith_init(&T, el0+n, ULONG_MAX, 1, n);
124 247044 : for (i = 1; i <= n_el; i++)
125 : {
126 : ulong el;
127 148510 : do el = u_forprime_next(&T); while (dvdiu(d1, el));
128 148510 : v[i] = el;
129 : }
130 98534 : return v;
131 : }
132 :
133 : /* return min el s.t. 2^63<el and el=1 (mod n). */
134 : static ulong
135 49267 : start_el_n(ulong n)
136 : {
137 49267 : ulong MAXHLONG = 1UL<<(BITS_IN_LONG-1), el = (MAXHLONG/n)*n + 1;
138 49267 : if ((el&1)==0) el += n; /* if el is even, then n is odd */
139 49267 : return el + (n << 1);
140 : }
141 :
142 : /* start probably catches d0*T_k(x). So small second is enough. */
143 : static ulong
144 49267 : get_n_el(GEN d0, ulong *psec)
145 : {
146 49267 : ulong start = ((lgefint(d0)-2)*BITS_IN_LONG)/(BITS_IN_LONG-1)+1, second = 1;
147 49267 : if (start>10) second++;
148 49267 : if (start>100) { start++; second++; }
149 49267 : if (start>500) { start++; second++; }
150 49267 : if (start>1000) { start++; second++; }
151 49267 : *psec = second; return start;
152 : }
153 :
154 : static long
155 77072 : FpX_degsub(GEN P, GEN Q, GEN p)
156 : {
157 77072 : pari_sp av = avma;
158 77072 : long d = degpol(FpX_sub(P, Q, p));
159 77072 : 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 100456 : gausspol_el(GEN H, ulong n, ulong d, ulong f, ulong g, ulong el)
170 : {
171 100456 : ulong j, k, z_n = rootsof1_Fl(n, el);
172 100456 : GEN vz_n, L = cgetg(1+f, t_VECSMALL), x = cgetg(1+f, t_VECSMALL), X;
173 :
174 100456 : vz_n = Fl_powers(z_n, n-1, el)+1;
175 366316 : for (k = 0; k < f; k++)
176 : {
177 265860 : ulong gk = Fl_powu(g, k, n), t = 0;
178 3725628 : for (j = 1; j <= d; j++)
179 3459768 : t = Fl_add(t, vz_n[Fl_mul(H[j], gk, n)], el);
180 265860 : L[1+k] = t;
181 265860 : x[1+(k+f-1)%f] = t;
182 : }
183 100456 : X = Flv_invVandermonde(L, 1, el);
184 100456 : return Flv_to_Flx(Flm_Flc_mul(X, x, el), 0);
185 : }
186 :
187 : static GEN
188 66920 : get_G(GEN H, GEN d0, GEN d1, GEN N, long k, ulong *pel, GEN *pM)
189 : {
190 66920 : long n = N[1], d = N[2], f = N[3], g = N[4], i;
191 66920 : GEN POL = cgetg(1+k, t_VEC), EL, G, M, x;
192 : pari_timer ti;
193 :
194 66920 : if (DEBUGLEVEL >= 6) timer_start(&ti);
195 66920 : EL = list_el_n(*pel, n, d1, k);
196 167376 : for (i = 1; i <= k; i++)
197 : {
198 100456 : ulong el = uel(EL,i);
199 100456 : x = gausspol_el(H, n, d, f, g, el);
200 100456 : gel(POL, i) = Flx_Fl_mul(x, umodiu(d0, el), el);
201 : }
202 66920 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "get_G : make data k=%ld",k);
203 66920 : G = nxV_chinese_center(POL, EL, &M);
204 66920 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "get_G : nxV_chinese_center k=%ld",k);
205 66920 : *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 49267 : get_d0_d1(GEN T, GEN P)
231 : {
232 49267 : long i, l = lg(P);
233 49267 : GEN x, y, dT = ZX_disc(T);
234 49267 : x = y = dT; setsigne(dT, 1);
235 114451 : for (i = 1; i < l; i++)
236 65184 : if (odd(Z_lvalrem(dT, P[i], &dT)))
237 : {
238 50739 : x = diviuexact(x, P[i]);
239 50739 : y = muliu(y, P[i]);
240 : }
241 49267 : return mkvec2(sqrti(x), sqrti(y)); /* x and y are square */
242 : }
243 :
244 : static GEN
245 33460 : gausspol(GEN T, GEN H, GEN N, ulong d, ulong f, ulong g)
246 : {
247 33460 : long n = N[1], el0 = N[2];
248 33460 : 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 33460 : d0 = gel(d0d1,1); /* d0*F is in Z[X] */
253 33460 : d1 = gel(d0d1,2); /* d1 has same prime factors as disc(T) */
254 33460 : Data = mkvecsmall4(n, d, f, g);
255 33460 : start = get_n_el(d0, &second);
256 33460 : el = start_el_n(n);
257 :
258 33460 : if (DEBUGLEVEL >= 6) timer_start(&ti);
259 33460 : if (DEBUGLEVEL == 2) err_printf("gausspol:start=(%ld,%ld)\n",start,second);
260 33460 : G1 = get_G(H, d0, d1, Data, start, &el, &M1);
261 33460 : for (n_el=second; n_el; n_el++)
262 : {
263 : GEN m, G2, M2;
264 33460 : G2 = get_G(H, d0, d1, Data, n_el, &el, &M2);
265 33460 : 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 33460 : F = RgX_Rg_div(G1, d0);
275 33460 : 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 33460 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "gausspol");
279 33460 : 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 48054 : mk_v_t_el(GEN vT, GEN Data, ulong el)
286 : {
287 48054 : pari_sp av = avma;
288 48054 : GEN H = gel(Data, 1), GH = gel(Data,2), i_t = gel(Data, 3), N=gel(Data, 6);
289 48054 : ulong n = N[1], d = N[2], mitk = N[5], f = N[3], i, k;
290 48054 : ulong z_n = rootsof1_Fl(n, el);
291 48054 : GEN vz_n = Fl_powers(z_n, n-1, el)+1;
292 48054 : GEN v_t_el = const_vecsmall(n-1, 0);
293 :
294 348257 : for (k = 1; k <= mitk; k++)
295 : {
296 300203 : if (k > 1 && !isintzero(gel(vT, k))) continue; /* k=1 is always handled */
297 2222451 : for (i=1; i<=f; i++)
298 : {
299 1922248 : ulong j, t = 0, x = Fl_mul(k, GH[i], n);
300 1922248 : long y = i_t[x]; /* x!=0, y!=0 */
301 1922248 : if (v_t_el[y]) continue;
302 5435404 : for (j = 1; j <= d; j++) t = Fl_add(t, vz_n[Fl_mul(x, H[j], n)], el);
303 262980 : v_t_el[y] = t;
304 : }
305 : }
306 48054 : return gerepileuptoleaf(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 31614 : get_vG(GEN vT, GEN Data, long n_el, ulong *pel, GEN *pM)
313 : {
314 31614 : GEN GH = gel(Data, 2), i_t = gel(Data, 3);
315 31614 : GEN d0 = gmael(Data, 4, 1), d1 = gmael(Data, 4, 2);
316 31614 : GEN kT = gel(Data, 5), N = gel(Data, 6);
317 31614 : long n = N[1], f = N[3], n_T = N[4], mitk = N[5];
318 31614 : GEN G = const_vec(mitk, gen_0), vPOL = cgetg(1+mitk, t_VEC);
319 : GEN EL, M, X, v_t_el;
320 31614 : GEN L = cgetg(1+f, t_VECSMALL), x = cgetg(1+f, t_VECSMALL);
321 : long i, j, k;
322 :
323 216870 : for (k=1; k<=mitk; k++) gel(vPOL, k) = cgetg(1+n_el, t_VEC);
324 31614 : EL = list_el_n(*pel, n, d1, n_el);
325 79668 : for (i=1; i<=n_el; i++)
326 : {
327 48054 : ulong el = uel(EL,i), d0model = umodiu(d0, el);
328 48054 : v_t_el = mk_v_t_el(vT, Data, el);
329 214716 : for (j = 1; j <= f; j++) L[j] = v_t_el[i_t[GH[j]]];
330 48054 : X = Flv_invVandermonde(L, 1, el);
331 229086 : for (k = 1; k <= n_T; k++)
332 : {
333 : GEN y;
334 181032 : long itk = kT[k];
335 181032 : if (!isintzero(gel(vT, itk))) continue;
336 700581 : for (j = 1; j <= f; j++) x[j] = v_t_el[i_t[Fl_mul(itk, GH[j], n)]];
337 133705 : y = Flv_to_Flx(Flm_Flc_mul(X, x, el), 0);
338 133705 : gmael(vPOL, itk, i) = Flx_Fl_mul(y, d0model, el);
339 : }
340 : }
341 150346 : for (k = 1; k <= n_T; k++)
342 : {
343 118732 : long itk = kT[k];
344 118732 : if (!isintzero(gel(vT, itk))) continue;
345 87224 : gel(G, itk) = nxV_chinese_center(gel(vPOL, itk), EL, &M);
346 : }
347 31614 : *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 15807 : get_i_t(long n, long p)
355 : {
356 : long a;
357 15807 : GEN v_t = const_vecsmall(n-1, 0);
358 15807 : GEN i_t = cgetg(n, t_VECSMALL); /* access i_t[a] with 1 <= a <= n-1 */
359 118319 : for (a = 1; a < n; a++)
360 : {
361 : long b;
362 854693 : while (a < n && v_t[a]) a++;
363 118319 : if (a==n) break;
364 102512 : b = a;
365 838886 : do { v_t[b] = 1; i_t[b] = a; b = Fl_mul(b, p, n); } while (b != a);
366 : }
367 15807 : 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 15807 : get_vT(GEN Data, int NEW)
375 : {
376 15807 : pari_sp av = avma;
377 15807 : GEN d0 = gmael(Data, 4, 1), kT = gel(Data, 5), N = gel(Data, 6);
378 15807 : ulong k, n = N[1], n_T = N[4], mitk = N[5];
379 15807 : GEN G1, M1, vT = const_vec(mitk, gen_0); /* vT[k]!=NULL ==> vT[k]=T_k */
380 15807 : ulong n_k = 0, el, n_el, start, second;
381 : pari_timer ti;
382 :
383 15807 : if (DEBUGLEVEL >= 6) timer_start(&ti);
384 15807 : if (!NEW) { gel(vT, 1) = pol_x(0); n_k++; }
385 15807 : start = get_n_el(d0, &second);
386 15807 : el = start_el_n(n);
387 15807 : if (DEBUGLEVEL == 2) err_printf("get_vT: start=(%ld,%ld)\n",start,second);
388 15807 : G1 = get_vG(vT, Data, start, &el, &M1);
389 15807 : for (n_el = second;; n_el++)
390 0 : {
391 : GEN G2, M2, m, m2;
392 15807 : G2 = get_vG(vT, Data, n_el, &el, &M2);
393 15807 : m = mulii(M1, M2); m2 = shifti(m,-1);
394 75173 : for (k = 1; k <= n_T; k++)
395 : {
396 59366 : long j = kT[k];
397 59366 : if (!isintzero(gel(vT,j))) continue;
398 43612 : if (FpX_degsub(gel(G1,j), gel(G2,j), M2) < 0) /* G1=G2 (mod M2) */
399 : {
400 43612 : gel(vT,j) = RgX_Rg_div(gel(G1,j), d0);
401 43612 : n_k++;
402 43612 : 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 15807 : if (n_k == n_T) break;
414 0 : M1 = m;
415 : }
416 15807 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "get_vT");
417 15807 : return gerepilecopy(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 15754 : get_kT(GEN i_t, long d)
424 15754 : { return vecsmall_uniq(vecsmall_shorten(i_t, d)); }
425 :
426 : static GEN
427 53 : get_kT_all(GEN GH, GEN i_t, long n, long d, long m)
428 : {
429 53 : long i, j, k = 0;
430 53 : GEN x = const_vecsmall(d*m, 0);
431 106 : for (i = 1; i <= m; i++)
432 755 : for (j = 1; j <= d; j++) x[++k] = i_t[Fl_mul(GH[i], j, n)];
433 53 : return vecsmall_uniq(x);
434 : }
435 :
436 : static GEN
437 53 : kT_from_kt_new(GEN gGH, GEN kt, GEN i_t, long n)
438 : {
439 53 : GEN EL = gel(gGH, 1);
440 53 : long i, k = 0, lEL = lg(EL), lkt = lg(kt);
441 53 : GEN x = cgetg(lEL+lkt, t_VECSMALL);
442 :
443 151 : for (i = 1; i < lEL; i++) x[++k] = i_t[EL[i]];
444 424 : for (i = 2; i < lkt; i++) if (n%kt[i]==0) x[++k] = kt[i];
445 53 : setlg(x, k+1); return vecsmall_uniq(x);
446 : }
447 :
448 : static GEN
449 53 : get_kTdiv(GEN kT, long n)
450 : {
451 53 : long i, k = 0, l = lg(kT);
452 53 : GEN div = cgetg(l, t_VECSMALL);
453 202 : for (i = 1; i < l; i++) if (n%kT[i]==0) div[++k] = kT[i];
454 53 : setlg(div, k+1);
455 53 : 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 833 : ZpX_roots_nonsep(GEN T, GEN R0, GEN p, GEN ps, GEN ps1)
462 : {
463 833 : long i, j, n = 0, lr = lg(R0);
464 833 : GEN v = cgetg(lr, t_VEC), R;
465 4375 : for (i = 1; i < lr; i++)
466 : {
467 3542 : GEN z, f = ZX_unscale_div(ZX_translate(T, gel(R0, i)), ps);
468 3542 : (void)ZX_pvalrem(f, p, &f);
469 3542 : gel(v, i) = z = FpX_roots(f, p);
470 3542 : n += lg(z)-1;
471 : }
472 833 : R = cgetg(n+1, t_VEC); n = 0;
473 4375 : for (i = 1; i < lr; i++)
474 : {
475 3542 : GEN z = gel(v, i);
476 3542 : long lz = lg(z);
477 8190 : for (j=1; j<lz; j++)
478 4648 : gel(R, ++n) = Fp_add(gel(R0, i), mulii(gel(z, j), ps), ps1);
479 : }
480 833 : return ZV_sort_uniq_shallow(R);
481 : }
482 : static GEN
483 49267 : ZpX_roots_all(GEN T, GEN p, long f, long *ptrs)
484 : {
485 49267 : pari_sp av = avma;
486 : pari_timer ti;
487 : GEN v, ps, ps1;
488 : long s;
489 :
490 49267 : if (DEBUGLEVEL >= 6) timer_start(&ti);
491 49267 : v = FpX_roots(T, p); /* FpX_roots returns sorted roots */
492 49267 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "FpX_roots, deg=%ld", degpol(T));
493 49267 : ps = NULL; ps1 = p;
494 50100 : for (s = 1; lg(v) != f+1; s++)
495 : {
496 833 : ps = ps1; ps1 = mulii(ps1, p); /* p^s, p^(s+1) */
497 833 : v = ZpX_roots_nonsep(T, v, p, ps, ps1);
498 833 : if (gc_needed(av, 1)) gerepileall(av, 3, &v, &ps, &ps1);
499 : }
500 49267 : *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 2513 : ZX_Zp_liftroots(GEN T, GEN v, GEN p, long r, long n)
508 : {
509 2513 : long i, l = lg(v);
510 2513 : GEN R = cgetg(l, t_VEC);
511 2513 : GEN pr = powiu(p, r), pn = powiu(p, n);
512 : pari_timer ti;
513 :
514 2513 : if (DEBUGLEVEL >= 6) timer_start(&ti);
515 9709 : for (i=1; i<l; i++)
516 : {
517 7196 : GEN x = gel(v, i), y, z;
518 7196 : GEN f = ZX_unscale_div(ZX_translate(T, x), pr);
519 7196 : (void)ZX_pvalrem(f, p, &f);
520 7196 : y = FpX_roots(f, p);
521 7196 : z = (n==r+1) ? y : ZX_Zp_root(f, gel(y, 1), p, n-r);
522 7196 : if (lg(y)!=2 || lg(z)!=2)
523 0 : pari_err_BUG("ZX_Zp_liftroots, roots are not separable");
524 7196 : gel(R, i) = Fp_add(x, mulii(gel(z, 1), pr), pn);
525 : }
526 2513 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "ZX_Zp_liftroots");
527 2513 : return R;
528 : }
529 :
530 : static GEN
531 33460 : set_R(GEN T, GEN F, GEN Rs, GEN p, long nf, long r, long s, long u)
532 : {
533 : long i;
534 33460 : GEN x, pr = powiu(p, r), prs = powiu(p, r+s), R = cgetg(1+nf, t_VEC), Rrs;
535 33460 : Rrs = r ? ZX_Zp_liftroots(T, Rs, p, s, r+s) : Rs;
536 33460 : x = gel(Rrs, 1);
537 121856 : for (i = 1; i <= nf; i++)
538 : {
539 88396 : x = FpX_eval(F, x, prs);
540 88396 : if (r) x = gel(Rrs, ZV_search(Rs, diviiexact(x, pr)));
541 88396 : gel(R, i) = x; /* R[i]=t_1^(g^i), t_1=Rrs[1], mod p^(r+s) */
542 : }
543 33460 : if (r+s < u) R = ZX_Zp_liftroots(T, R, p, r+s, u);
544 32116 : else if (r+s > u) R = FpV_red(R, powiu(p, u));
545 33460 : 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 34222 : set_e0_e1(ulong el, ulong e, GEN p)
575 : {
576 34222 : ulong s, d0, e0, e1, f0, n, phin, g, up = itou_or_0(p);
577 :
578 34222 : 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 34222 : ulong pmodel = up ? up%el : umodiu(p, el);
600 34222 : d0 = Fl_order(pmodel, el-1, el);
601 34222 : s = Z_lval(subiu(powiu(p, d0), 1), el);
602 34222 : if (e <= s) { e0 = e; e1 = 0; } else { e0 = s; e1 = e-s; }
603 34222 : phin = (el-1)*upowuu(el, e0-1);
604 : }
605 34222 : n = upowuu(el, e0); f0 = phin/d0;
606 34222 : g = (el==2) ? 1 : pgener_Zl(el);
607 34222 : 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 78511 : use_newton(long d, long f)
613 : {
614 78511 : if (2*d <= f) return 0;
615 73603 : else if (f <= d) return 1;
616 5935 : 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 14 : use_newton_general(long d, long f, long maxdeg)
629 : {
630 14 : if (maxdeg < 20) return 0;
631 14 : else if (f <= 50) return 1;
632 7 : else if (maxdeg < 30) return 0;
633 7 : else if (f <= 60) return 1;
634 7 : else if (maxdeg < 40) return 0;
635 7 : else if (f <= 70) return 1;
636 7 : else if (maxdeg < 50) return 0;
637 7 : else if (f <= 80) return 1;
638 7 : 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 7 : use_general(long d, long maxdeg)
648 : {
649 7 : if (d <= 50) return 1;
650 7 : else if (maxdeg <= 3*d) return 0;
651 7 : 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 70 : update_dfm(long *pd, long *pf, long *pm, long di, long fi)
661 : {
662 70 : long c = ugcd(*pd,di), d1 = *pd * di, f1 = *pf * fi;
663 70 : *pd = d1 / c; *pf = c * f1; *pm += d1 * d1 * f1;
664 70 : if (DEBUGLEVEL == 4) err_printf("(%ld,%ld), ",d1,f1);
665 70 : }
666 : /* assume ord(p mod f) > 1 */
667 : static ulong
668 94605 : set_action(GEN fn, GEN p, long d, long f)
669 : {
670 94605 : GEN EL = gel(fn, 1), E = gel(fn, 2);
671 94605 : long i, d0, f0, m0, m1, maxdeg, max, l = lg(EL);
672 94605 : ulong action = 0, up = itou_or_0(p);
673 94605 : GEN D = cgetg(l, t_VECSMALL), F = cgetg(l, t_VECSMALL);
674 :
675 94605 : d += 10*(lgefint(p)-3);
676 94605 : if (l == 2)
677 : { /* n is a prime power */
678 47670 : action |= (EL[1]==2 || !use_newton(d, f))? GENERAL: NEWTON_POWER;
679 47670 : return action;
680 : }
681 46935 : if (f <= 2) action |= NEWTON_GENERAL;
682 36540 : else if (d <= 10) action |= GENERAL;
683 5835 : else if (f <= 10) action |= NEWTON_GENERAL;
684 476 : else if (d <= 20) action |= GENERAL;
685 73 : else if (f <= 20) action |= NEWTON_GENERAL_NEW;
686 27 : else if (d <= 30) action |= GENERAL;
687 21 : else if (f <= 30) action |= NEWTON_GENERAL_NEW;
688 21 : else if (d <= 40) action |= GENERAL;
689 21 : else if (f <= 40) action |= NEWTON_GENERAL_NEW;
690 46935 : if (action) return action;
691 : /* can assume that d > 40 and f > 40 */
692 :
693 21 : maxdeg = max = 1;
694 91 : for (i = 1; i < l; i++)
695 : {
696 70 : long x, el = EL[i], e = E[i];
697 70 : long q = upowuu(el, e-1), ni = q * el, phini = ni - q;
698 70 : long di = Fl_order(umodiu(p, ni), phini, ni);
699 70 : D[i] = di; F[i] = phini / di;
700 70 : x = ugcd(max, di); max = max * (di / x); /* = lcm(d1,..di) */
701 70 : if (x > 1) maxdeg = max*x;
702 70 : if (DEBUGLEVEL == 4) err_printf("(%ld,%ld), ", D[i], F[i]);
703 : }
704 21 : if (maxdeg == 1) return action;
705 14 : if (up != 2)
706 : {
707 14 : if (use_newton_general(d, f, maxdeg))
708 : { /* does not decompose n */
709 7 : action |= (20 < d)? NEWTON_GENERAL_NEW: NEWTON_GENERAL;
710 7 : return action;
711 : }
712 7 : if (use_general(d, maxdeg)) action |= GENERAL;
713 : }
714 7 : if (l < 4) return action; /* n has two factors */
715 :
716 7 : d0 = f0 = 1; m0 = 0;
717 42 : for (i = 1; i < l; i++) update_dfm(&d0, &f0, &m0, D[i], F[i]);
718 7 : if (DEBUGLEVEL == 4) err_printf("(d0,f0)=(%ld,%ld)\n",d0,f0);
719 7 : d0 = f0 = 1; m1 = 0;
720 42 : for (i = l-1; i >= 1; i--) update_dfm(&d0, &f0, &m1, D[i], D[i]);
721 7 : if (DEBUGLEVEL == 4) err_printf("(d0,f0)=(%ld,%ld)\n",d0,f0);
722 7 : if (DEBUGLEVEL == 4) err_printf("(m0,m1)=(%lu,%lu) %ld\n",m0,m1,m0<=m1);
723 7 : if (m0 <= m1) action |= ASCENT;
724 7 : return action;
725 : }
726 :
727 : static GEN
728 10758 : FpX_Newton_perm(long d, GEN R, GEN v, GEN pu, GEN p)
729 : {
730 10758 : GEN S = cgetg(d+2, t_VEC);
731 : long k;
732 95189 : gel(S,1) = utoi(d); for (k = 1; k <= d; k++) gel(S, k+1) = gel(R, v[k]);
733 10758 : return FpX_red(FpX_fromNewton(RgV_to_RgX(S, 0), pu), p);
734 : }
735 : static GEN
736 81847 : Flx_Newton_perm(long d, GEN R, GEN v, ulong pu, ulong p)
737 : {
738 81847 : GEN S = cgetg(d+2, t_VEC);
739 : long k;
740 1177903 : S[1] = d; for (k = 1; k <= d; k++) gel(S, k+1) = gel(R, v[k]);
741 81847 : 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 6866 : FpX_pol_newton_general(GEN Data, GEN N2, GEN vT, GEN x)
748 : {
749 6866 : GEN i_t = gel(Data, 3), kT = gel(Data, 5), N = gel(Data, 6);
750 6866 : long k, d = N[2], n_T = N[4], mitk = N[5];
751 6866 : GEN p = gel(N2,1), pr = gel(N2,2), pu = gel(N2,3), pru = gel(N2,4);
752 6866 : GEN R = cgetg(1+mitk, t_VEC);
753 :
754 32579 : for (k = 1; k <= n_T; k++)
755 25713 : gel(R, kT[k]) = diviiexact(FpX_eval(gel(vT, kT[k]), x, pru), pr);
756 6866 : 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 2414 : FpX_factcyclo_newton_general(GEN Data)
765 : {
766 2414 : GEN GH = gel(Data, 2), fn = gel(Data, 3), p = gel(Data, 4);
767 2414 : long n = mael(Data, 5, 1), d = mael(Data, 5, 2), f = mael(Data, 5, 3);
768 2414 : long m = mael(Data, 5, 4), pmodn = umodiu(p, n);
769 2414 : 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 2414 : if (m != 1) m = f;
774 2414 : for (pu = p; cmpiu(pu,d) <= 0; u++) pu = mulii(pu, p); /* d<pu, pu=p^n */
775 :
776 2414 : H = Fl_powers(pmodn, d-1, n); /* H=<p> */
777 2414 : i_t = get_i_t(n, pmodn); /* i_t[1+i]=k ==> t_i=t_k */
778 2414 : kT = get_kT(i_t, d); n_T = lg(kT)-1; mitk = kT[n_T];
779 2414 : if (DEBUGLEVEL == 2) err_printf("kT=%Ps %ld elements\n",kT,n_T);
780 2414 : if (DEBUGLEVEL >= 6) timer_start(&ti);
781 2414 : T = galoissubcyclo(utoi(n), utoi(pmodn), 0, 0);
782 2414 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "galoissubcyclo");
783 2414 : d0d1 = get_d0_d1(T, gel(fn,1)); /* d0*T_k(x) is in Z[x] */
784 2414 : Data2 = mkvecn(6, H, GH, i_t, d0d1, kT, mkvecsmalln(5, n, d, f, n_T, mitk));
785 2414 : vT = get_vT(Data2, 0);
786 2414 : if (DEBUGLEVEL == 4) err_printf("vT=%Ps\n",vT);
787 2414 : r = QXV_den_pval(vT, kT, p);
788 2414 : Rs = ZpX_roots_all(T, p, f, &s);
789 2414 : if (DEBUGLEVEL >= 2) err_printf("(u,s,r)=(%ld,%ld,%ld)\n",u,s,r);
790 2414 : 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 2414 : R = (r+u==s) ? Rs : ZX_Zp_liftroots(T, Rs, p, s, r+u);
793 2414 : pr = powiu(p, r); pru = powiu(p, r+u); /* Usually, r=0, s=1, pr=1, pru=p */
794 11320 : for (k = 1; k<=n_T; k++)
795 : {
796 8906 : long itk = kT[k];
797 8906 : GEN z = r? RgX_Rg_mul(gel(vT, itk), pr): gel(vT, itk);
798 8906 : gel(vT, itk) = RgX_to_FpX(z, pru);
799 : }
800 2414 : Data3 = mkvec4(p, pr, pu, pru);
801 2414 : if (DEBUGLEVEL >= 6) timer_start(&ti);
802 9280 : for (i=1; i<=m; i++)
803 6866 : gel(R,i) = FpX_pol_newton_general(Data2, Data3, vT, gel(R,i));
804 2414 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "FpX_pol_newton_general");
805 2414 : 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 21024 : Fl_next_gen3(long x, long i, GEN v_t_p, ulong t, GEN Data)
871 : {
872 21024 : GEN vT = gel(Data, 1), gGH = gel(Data, 2), Rs = gel(Data, 3);
873 21024 : GEN Rrs = gel(Data, 4), i_t = gel(Data, 5);
874 21024 : long pu = mael(Data, 8, 2), pr = mael(Data, 9, 2), prs = mael(Data, 10, 2);
875 21024 : GEN EL = gel(gGH, 1), E = gel(gGH, 2), div = gel(Data, 12);
876 21024 : long n = mael(Data, 11, 1), r = mael(Data, 11, 2), mitk = mael(Data, 11, 5);
877 21024 : long j, k, l = lg(EL), ld = lg(div);
878 21024 : if (i >= l) return;
879 27936 : for (j = 0; j < E[i]; j++)
880 : {
881 20736 : if (j > 0)
882 : {
883 13536 : long itx = i_t[x];
884 13536 : t = Flx_eval(gel(vT, i_t[EL[i]]), t, prs); /* mod p^(r+s) */
885 13536 : if (r) t = Rrs[zv_search(Rs, t/pr)]; /* mod p^(r+s) */
886 13536 : if (itx <= mitk) v_t_p[itx] = t%pu; /* mod p^u */
887 54144 : for (k = 1; k < ld; k++)
888 : {
889 40608 : ulong y = Fl_mul(x, div[k], n), v;
890 40608 : long ity = i_t[y];
891 40608 : if (ity > mitk || v_t_p[ity]) continue;
892 2592 : v = Flx_eval(gel(vT, i_t[div[k]]), t, prs); /* mod p^(r+s) */
893 2592 : if (r) v /= pr; /* mod p^s */
894 2592 : v_t_p[ity] = v%pu;
895 : }
896 : }
897 20736 : Fl_next_gen3(x, i+1, v_t_p, t, Data);
898 20736 : 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 288 : 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 288 : GEN Rs = gel(Data, 3), Rrs = gel(Data, 4);
908 288 : ulong pu = mael(Data, 8, 2), pr = mael(Data, 9, 2), prs = mael(Data, 10, 2);
909 288 : GEN vT = gel(Data, 1), i_t = gel(Data, 5), div = gel(Data, 12);
910 288 : long k, r = mael(Data, 11, 2), mitk = mael(Data, 11, 5), ld = lg(div);
911 288 : GEN v_t_p = const_vecsmall(mitk, 0);
912 :
913 288 : v_t_p[1] = Rs[i] % pu; /* mod p^u, guaranteed u<=s */
914 288 : Fl_next_gen3(1, 1, v_t_p, Rrs[i], Data);
915 1152 : for (k = 1; k < ld; k++)
916 : {
917 864 : ulong itk = i_t[div[k]], x = Flx_eval(gel(vT, itk), Rrs[i], prs);
918 864 : if (r) x /= pr; /* mod p^s */
919 864 : v_t_p[itk] = x % pu;
920 : }
921 288 : return v_t_p;
922 : }
923 :
924 : /* Data = [GHgen, GH, fn, p, [n, d, f, m]] */
925 : static GEN
926 53 : newton_general_new_pre3(GEN Data)
927 : {
928 53 : GEN gGH = gel(Data, 1), GH = gel(Data, 2), fn = gel(Data, 3);
929 53 : GEN p = gel(Data, 4);
930 53 : long n = mael(Data, 5, 1), d = mael(Data, 5, 2), f = mael(Data, 5, 3);
931 53 : long pmodn = umodiu(p, n);
932 53 : 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 60 : for (pu = p; cmpiu(pu,d)<=0; u++) pu = mulii(pu, p); /* d<pu, pu=p^u */
938 :
939 53 : H = Fl_powers(pmodn, d-1, n); /* H=<p> */
940 53 : i_t = get_i_t(n, pmodn); /* i_t[1+i]=k ==> t_i=t_k */
941 53 : kt = get_kT_all(GH, i_t, n, d, 1); n_t = lg(kt)-1; mitk = kt[n_t];
942 53 : kT = kT_from_kt_new(gGH, kt, i_t, n); n_T = lg(kT)-1; miTk = kT[n_T];
943 53 : kTdiv = get_kTdiv(kT, n);
944 53 : if (DEBUGLEVEL == 2)
945 0 : err_printf("kt=%Ps %ld elements\nkT=%Ps %ld elements\n",kt,n_t,kT,n_T);
946 53 : if (DEBUGLEVEL == 2)
947 0 : err_printf("kTdiv=%Ps\n",kTdiv);
948 :
949 53 : if (DEBUGLEVEL >= 6) timer_start(&ti);
950 53 : T = galoissubcyclo(utoi(n), utoi(pmodn), 0, 0);
951 53 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "galoissubcyclo");
952 53 : d0d1 = get_d0_d1(T, gel(fn,1)); /* d0*T_k(x) is in Z[x] */
953 53 : Data2 = mkvecn(6, H, GH, i_t, d0d1, kT, mkvecsmalln(5, n, d, f, n_T, miTk));
954 53 : vT = get_vT(Data2, 1);
955 53 : if (DEBUGLEVEL == 4) err_printf("vT=%Ps\n",vT);
956 53 : r = QXV_den_pval(vT, kT, p);
957 53 : Rs = ZpX_roots_all(T, p, f, &s);
958 53 : if (DEBUGLEVEL >= 2) err_printf("(u,s,r)=(%ld,%ld,%ld)\n",u,s,r);
959 53 : 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 53 : Rrs = r ? ZX_Zp_liftroots(T, Rs, p, s, r+s) : Rs;
966 53 : pr = powiu(p, r); prs = powiu(p, r+s); /* Usually, r=0, s=1, pr=1, pru=p */
967 53 : 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 6 : ulong upr = itou(pr), uprs = itou(prs);
979 36 : for (k = 1; k <= n_T; k++)
980 : {
981 30 : long itk = kT[k];
982 30 : GEN z = r? RgX_muls(gel(vT, itk), upr): gel(vT, itk);
983 30 : gel(vT, itk) = RgX_to_Flx(z, uprs);
984 : }
985 6 : Rs = ZV_to_nv(Rs); Rrs = ZV_to_nv(Rrs);
986 : }
987 53 : 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 4483 : FpX_1st_lift_2(GEN z, GEN p)
1026 : {
1027 4483 : long i, r = lg(z);
1028 4483 : GEN x = cgetg(r, t_POL);
1029 4483 : x[1] = evalsigne(1) | evalvarn(0);
1030 4483 : if (odd(r))
1031 9962 : for (i = 2; i < r; i++) gel(x,i) = odd(i)? Fp_neg(gel(z,i), p): gel(z,i);
1032 : else
1033 16909 : for (i = 2; i < r; i++) gel(x,i) = odd(i)? gel(z,i): Fp_neg(gel(z,i), p);
1034 4483 : return x;
1035 : }
1036 :
1037 : static GEN
1038 3190 : FpX_1st_lift(GEN z, GEN p, ulong e, ulong el, GEN vP)
1039 : {
1040 3190 : GEN z2, z3, P = gel(vP, e);
1041 3190 : if (!gel(vP, e)) P = gel(vP, e) = FpX_polcyclo(e, p);
1042 3190 : z2 = RgX_inflate(z, el);
1043 3190 : z3 = FpX_normalize(FpX_gcd(P, z2, p), p);
1044 3190 : return FpX_div(z2, z3, p);
1045 : }
1046 :
1047 : static GEN
1048 11807 : FpX_lift(GEN z, GEN p, ulong e, ulong el, ulong r, ulong s, GEN vP)
1049 : {
1050 11807 : if (s == 0)
1051 : {
1052 7673 : z = (el==2) ? FpX_1st_lift_2(z, p) : FpX_1st_lift(z, p, e, el, vP);
1053 7673 : if (r >= 2) z = RgX_inflate(z, upowuu(el, r-1));
1054 : }
1055 : else
1056 4134 : z = RgX_inflate(z, upowuu(el, r-s));
1057 11807 : return z;
1058 : }
1059 :
1060 : /* e is the conductor of the splitting field of p in Q(zeta_n) */
1061 : static GEN
1062 10501 : FpX_conductor_lift(GEN z, GEN p, GEN fn, ulong e, GEN vP)
1063 : {
1064 10501 : GEN EL = gel(fn, 1), En = gel(fn, 2);
1065 10501 : long i, r = lg(EL), new_e = e;
1066 :
1067 32077 : for (i = 1; i < r; i++)
1068 : {
1069 21576 : long el = EL[i], en = En[i], ee = u_lval(e, el);
1070 21576 : if (ee < en)
1071 : {
1072 11807 : z = FpX_lift(z, p, new_e, el, en, ee, vP);
1073 11807 : new_e *= upowuu(el, en-ee);
1074 : }
1075 : }
1076 10501 : return z;
1077 : }
1078 :
1079 : /* R0 is mod p^u, d < p^u */
1080 : static GEN
1081 3547 : FpX_pol_newton(long j, GEN R0, GEN E, GEN D3, long d, long f, GEN p)
1082 : {
1083 3547 : long i, u = D3[3];
1084 3547 : GEN R = cgetg(1+f, t_VEC);
1085 15377 : for (i = 1; i <= f; i++) gel(R, i) = gel(R0, 1+(i+j)%f);
1086 3547 : 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 1896 : FpX_factcyclo_newton_pre(GEN Data, GEN N, GEN p, ulong m)
1092 : {
1093 1896 : GEN T = gel(Data, 1), F = gel(Data, 2), Rs = gel(Data, 3), D4 = gel(Data, 4);
1094 1896 : 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 1896 : ulong i, n = N[1], pmodn = umodiu(p, n);
1097 : pari_timer ti;
1098 :
1099 1896 : if (m!=1) m = nf;
1100 1896 : pols = cgetg(1+m, t_VEC);
1101 1896 : E = set_E(pmodn, n, d, nf, g);
1102 1896 : R = set_R(T, F, Rs, p, nf, r, s, u);
1103 1896 : Data3 = mkvecsmall3(r, s, u);
1104 1896 : if (DEBUGLEVEL >= 6) timer_start(&ti);
1105 5443 : for (i = 1; i <= m; i++) gel(pols, i) = FpX_pol_newton(i, R,E,Data3, d,nf,p);
1106 1896 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "FpX_pol_newton");
1107 1896 : 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 7 : FpX_factcyclo_lift(long n1, GEN v1, long n2, GEN v2, GEN p, long m)
1115 : {
1116 7 : pari_sp av = avma;
1117 7 : long d1 = degpol(gel(v1, 1)), lv1 = lg(v1), nf1 = eulerphiu(n1)/d1;
1118 7 : long d2 = degpol(gel(v2, 1)), lv2 = lg(v2), nf2 = eulerphiu(n2)/d2;
1119 7 : long e = ugcd(d1, d2), nf = nf1*nf2*e, need_factor = (e!=1);
1120 7 : long i, j, k = 0;
1121 7 : GEN z = NULL, v;
1122 : pari_timer ti;
1123 :
1124 7 : if (DEBUGLEVEL >= 6) timer_start(&ti);
1125 7 : if (m != 1) m = nf;
1126 7 : v = cgetg(1+m, t_VEC);
1127 7 : 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 91 : for (i = 1; i < lv1; i++)
1137 420 : for (j = 1; j < lv2; j++)
1138 : {
1139 336 : GEN x1 = gel(v1, i), x2 = gel(v2, j);
1140 336 : z = FpX_gcd(RgX_inflate(x1, n2), RgX_inflate(x2, n1), p);
1141 336 : z = FpX_normalize(z, p); /* needed */
1142 336 : 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 7 : 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 7 : return gerepilecopy(av, v);
1155 : }
1156 :
1157 : /* n is any integer prime to p; d>1 and f>1 */
1158 : static GEN
1159 1349 : FpX_factcyclo_gen(GEN GH, ulong n, GEN p, long m)
1160 : {
1161 1349 : GEN fn = factoru(n), fa = zm_to_ZM(fn);
1162 : GEN A, T, X, pd_n, v, pols;
1163 1349 : long i, j, pmodn = umodiu(p, n), phin = eulerphiu_fact(fn);
1164 1349 : long d = Fl_order(pmodn, phin, n), f = phin/d;
1165 : pari_timer ti;
1166 :
1167 1349 : if (m != 1) m = f;
1168 1349 : if (m != 1 && GH==NULL) /* FpX_factcyclo_fact is used */
1169 : {
1170 381 : GEN H = znstar_generate(n, mkvecsmall(pmodn));
1171 381 : GH = znstar_cosets(n, phin, H); /* representatives of G/H */
1172 : }
1173 :
1174 1349 : pols = cgetg(1+m, t_VEC);
1175 1349 : v = cgetg(1+d, t_VEC);
1176 1349 : pd_n = diviuexact(subis(powiu(p, d), 1), n); /* (p^d-1)/n */
1177 1349 : T = init_Fq(p, d, 1);
1178 1349 : A = pol_x(1); /* A is a generator of F_q, q=p^d */
1179 1349 : if (d == 1) A = FpX_rem(A, T, p);
1180 1349 : random_FpX(degpol(T), varn(T), p); /* skip 1st one */
1181 1349 : if (DEBUGLEVEL >= 6) timer_start(&ti);
1182 1791 : do X = FpXQ_pow(random_FpX(degpol(T), varn(T), p), pd_n, T, p);
1183 1791 : while(lg(X)<3 || equaliu(FpXQ_order(X, fa, T, p), n)==0); /* find zeta_n */
1184 1349 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "find X");
1185 :
1186 1349 : if (m == 1)
1187 : {
1188 2854 : for (j = 1; j <= d; j++)
1189 : {
1190 2183 : gel(v, j) = pol_x(0);
1191 2183 : gmael(v, j, 2) = FpX_neg(X, p);
1192 2183 : if (j < d) X = FpXQ_pow(X, p, T, p);
1193 : }
1194 671 : gel(pols, 1) = ZXX_evalx0(FpXQXV_prod(v, T, p));
1195 : }
1196 : else
1197 : {
1198 : GEN vX, vp;
1199 678 : if (DEBUGLEVEL >= 6) timer_start(&ti);
1200 678 : vX = FpXQ_powers(X, n, T, p)+1;
1201 678 : vp = Fl_powers(pmodn, d, n);
1202 7884 : for (i = 1; i <= m; i++)
1203 : {
1204 68694 : for (j = 1; j <= d; j++)
1205 : {
1206 61488 : gel(v, j) = pol_x(0);
1207 61488 : gmael(v, j, 2) = FpX_neg(gel(vX, Fl_mul(GH[i], vp[j], n)), p);
1208 : }
1209 7206 : gel(pols, i) = ZXX_evalx0(FpXQXV_prod(v, T, p));
1210 : }
1211 678 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "FpXQXV_prod");
1212 : }
1213 1349 : 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 33460 : FpX_factcyclo_newton_power(GEN N, GEN p, ulong m, int small)
1221 : {
1222 : GEN Rs, H, T, F, pr, prs, pu, Data;
1223 33460 : long n = N[1], el = N[2], phin = N[4], g = N[5];
1224 33460 : long pmodn = umodiu(p, n), pmodel = umodiu(p, el);
1225 33460 : long d = Fl_order(pmodel, el-1, el), nf = phin/d;
1226 33460 : long r, s = 0, u = 1;
1227 :
1228 33460 : if (m != 1) m = nf;
1229 35700 : for (pu = p; cmpiu(pu,d) <= 0; u++) pu = mulii(pu,p); /* d<p^u, pu=p^u */
1230 33460 : H = Fl_powers(pmodn, d, n);
1231 33460 : T = galoissubcyclo(utoi(n), utoi(pmodn), 0, 0);
1232 33460 : F = gausspol(T, H, N, d, nf, g);
1233 33460 : r = QX_den_pval(F, p);
1234 33460 : Rs = ZpX_roots_all(T, p, nf, &s);
1235 33460 : if (DEBUGLEVEL >= 2) err_printf("(u,s,r)=(%ld,%ld,%ld)\n",u,s,r);
1236 33460 : pr = powiu(p, r); prs = powiu(p, r+s); /* Usually, r=0, s=1, pr=1, prs=p */
1237 33460 : F = r ? RgX_to_FpX(RgX_Rg_mul(F, pr), prs) : RgX_to_FpX(F, prs);
1238 33460 : Data = mkvec4(T, F, Rs, mkvecsmalln(6, d, nf, g, r, s, u));
1239 33460 : if (small && lgefint(pu) == 3)
1240 31564 : return Flx_factcyclo_newton_pre(Data, N, uel(p,2), m);
1241 : else
1242 1896 : return FpX_factcyclo_newton_pre(Data, N, p, m);
1243 : }
1244 :
1245 : static GEN
1246 2802 : FpX_split(ulong n, GEN p, ulong m)
1247 : {
1248 : ulong i, j;
1249 2802 : GEN v, C, vx, z = rootsof1u_Fp(n, p);
1250 2802 : if (m == 1) return mkvec(deg1pol_shallow(gen_1, Fp_neg(z,p), 0));
1251 1401 : v = cgetg(m+1, t_VEC); C = coprimes_zv(n); vx = Fp_powers(z, n-1, p);
1252 18495 : for (i = j = 1; i <= n; i++)
1253 17094 : if (C[i]) gel(v, j++) = deg1pol_shallow(gen_1, Fp_neg(gel(vx,i+1), p), 0);
1254 1401 : 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 2658 : FpX_factcyclo_prime_power_i(long el, long e, GEN p, long m)
1261 : {
1262 2658 : GEN z = set_e0_e1(el, e, p), v;
1263 2658 : long n = z[1], e0 = z[2], e1 = z[3], phin = z[4], g = z[5];
1264 2658 : long d = z[6], f = z[7]; /* d and f for n=el^e0 */
1265 :
1266 2658 : if (f == 1) v = mkvec(FpX_polcyclo(n, p));
1267 2658 : else if (d == 1) v = FpX_split(n, p, (m==1)? 1: f);
1268 2658 : else if (el == 2) v = FpX_factcyclo_gen(NULL, n, p, m); /* d==2 in this case */
1269 2658 : else if (!use_newton(d,f)) v = FpX_factcyclo_gen(NULL, n, p, m);
1270 : else
1271 : {
1272 1896 : GEN N = mkvecsmall5(n, el, e0, phin, g);
1273 1896 : v = FpX_factcyclo_newton_power(N, p, m, 0);
1274 : }
1275 2658 : 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 2658 : return v;
1281 : }
1282 : static GEN
1283 14 : FpX_factcyclo_prime_power(long el, long e, GEN p, long m)
1284 : {
1285 14 : pari_sp av = avma;
1286 14 : return gerepilecopy(av, FpX_factcyclo_prime_power_i(el, e, p, m));
1287 : }
1288 :
1289 : static GEN
1290 7 : FpX_factcyclo_fact(GEN fn, GEN p, ulong m, long ascent)
1291 : {
1292 7 : GEN EL = gel(fn, 1), E = gel(fn, 2), v1 = NULL;
1293 7 : long i, l = lg(EL), n1 = 1;
1294 :
1295 21 : for (i = 1; i < l; i++)
1296 : {
1297 14 : long j = ascent? i: l-i, n2 = upowuu(EL[j], E[j]);
1298 14 : GEN v2 = FpX_factcyclo_prime_power(EL[j], E[j], p, m);
1299 14 : v1 = v1? FpX_factcyclo_lift(n1, v1, n2, v2, p, m): v2;
1300 14 : n1 *= n2;
1301 : }
1302 7 : return v1;
1303 : }
1304 :
1305 : static GEN
1306 15807 : Flv_FlvV_factorback(GEN g, GEN x, ulong q)
1307 34270 : { 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 15807 : get_GH_gen(long n, long pmodn)
1313 : {
1314 15807 : GEN G = znstar0(utoipos(n), 1), cycG = znstar_get_cyc(G);
1315 : GEN cycGH, gG, gGH, Ui, P;
1316 : ulong expG;
1317 15807 : P = hnfmodid(mkmat(Zideallog(G, utoi(pmodn))), cycG);
1318 15807 : cycGH = ZV_to_nv(ZM_snf_group(P, NULL, &Ui));
1319 15807 : expG = itou(cyc_get_expo(cycG));
1320 15807 : gG = ZV_to_Flv(znstar_get_gen(G), n);
1321 15807 : gGH = Flv_FlvV_factorback(gG, ZM_to_Flm(Ui, expG), n);
1322 15807 : 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 94605 : FpX_factcyclo_just_conductor_init(GEN *pData, ulong n, GEN p, ulong m)
1345 : {
1346 94605 : GEN fn = factoru(n), GH = NULL, GHgen = NULL;
1347 94605 : long phin = eulerphiu_fact(fn), pmodn = umodiu(p, n);
1348 94605 : long d = Fl_order(pmodn, phin, n), f = phin/d; /* d > 1 */
1349 94605 : ulong action = set_action(fn, p, d, f);
1350 94605 : if (action & ~NEWTON_POWER)
1351 : { /* needed for GENERAL* */
1352 60390 : GEN H = znstar_generate(n, mkvecsmall(pmodn));
1353 60390 : GH = znstar_cosets(n, phin, H); /* representatives of G/H */
1354 60390 : if (action & (NEWTON_GENERAL_NEW | NEWTON_GENERAL))
1355 15807 : GHgen = get_GH_gen(n, pmodn); /* gen and order of G/H */
1356 : }
1357 94605 : *pData = mkvec5(GHgen, GH, fn, p, mkvecsmall4(n, d, f, m));
1358 94605 : 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 94605 : return action;
1368 : }
1369 :
1370 : static GEN
1371 5698 : FpX_factcyclo_just_conductor(ulong n, GEN p, ulong m)
1372 : {
1373 : GEN Data, fn;
1374 5698 : ulong action = FpX_factcyclo_just_conductor_init(&Data, n, p, m);
1375 5698 : fn = gel(Data,3);
1376 5698 : if (action & GENERAL)
1377 587 : return FpX_factcyclo_gen(gel(Data,2), n, p, m);
1378 5111 : else if (action & NEWTON_POWER)
1379 2644 : return FpX_factcyclo_prime_power_i(ucoeff(fn,1,1), ucoeff(fn,1,2), p, m);
1380 2467 : else if (action & NEWTON_GENERAL)
1381 2414 : return FpX_factcyclo_newton_general(Data);
1382 53 : else if (action & NEWTON_GENERAL_NEW)
1383 46 : return FpX_factcyclo_newton_general_new3(Data);
1384 : else
1385 7 : return FpX_factcyclo_fact(fn, p, m, action & ASCENT);
1386 : }
1387 :
1388 : static GEN
1389 10656 : FpX_factcyclo_i(ulong n, GEN p, long fl)
1390 : {
1391 10656 : GEN fn = factoru(n), z;
1392 10656 : long phin = eulerphiu_fact(fn), pmodn = umodiu(p, n);
1393 10656 : ulong d = Fl_order(pmodn, phin, n), f = phin/d, fK;
1394 :
1395 10656 : if (DEBUGLEVEL >= 1) header(fn, n, d, f, p);
1396 10656 : if (f == 1) { z = FpX_polcyclo(n, p); return fl? z: mkvec(z); }
1397 8500 : else if (d == 1) /* p=1 (mod n), zeta_n in Z_p */
1398 774 : { z = FpX_split(n, p, fl? 1: f); return fl? gel(z,1): z; }
1399 7726 : fK = znstar_conductor(znstar_generate(n, mkvecsmall(pmodn)));
1400 7726 : if (fK != n && umodiu(p, fK) == 1)
1401 2028 : z = FpX_split(fK, p, fl? 1: f);
1402 : else
1403 5698 : z = FpX_factcyclo_just_conductor(fK, p, fl? 1: f);
1404 7726 : if (n > fK)
1405 : {
1406 4294 : GEN vP = const_vec(n, NULL);
1407 4294 : long i, l = fl? 2: lg(z);
1408 14795 : for (i = 1; i < l; i++)
1409 10501 : gel(z, i) = FpX_conductor_lift(gel(z, i), p, fn, fK, vP);
1410 : }
1411 7726 : return fl? gel(z,1): gen_sort(z,(void*)cmpii, &gen_cmp_RgX);
1412 : }
1413 :
1414 : GEN
1415 42 : FpX_factcyclo(ulong n, GEN p, ulong m)
1416 42 : { pari_sp av = avma; return gerepilecopy(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 24102 : Flx_pol_newton_general(GEN Data, GEN N2, GEN vT, ulong x)
1422 : {
1423 24102 : GEN i_t = gel(Data, 3), kT = gel(Data, 5), N = gel(Data, 6);
1424 24102 : long k, d = N[2], n_T = N[4], mitk = N[5];
1425 24102 : long p = N2[1], pr = N2[2], pu = N2[3], pru = N2[4];
1426 24102 : GEN R = cgetg(1+mitk, t_VECSMALL);
1427 :
1428 123717 : for (k = 1; k <= n_T; k++) uel(R,kT[k]) = Flx_eval(gel(vT, kT[k]), x, pru) / pr;
1429 24102 : 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 13340 : Flx_factcyclo_newton_general(GEN Data)
1438 : {
1439 13340 : GEN GH = gel(Data, 2), fn = gel(Data, 3), p = gel(Data, 4);
1440 13340 : ulong up = p[2], n = mael(Data, 5, 1), pmodn = up%n;
1441 13340 : long d = mael(Data, 5, 2), f = mael(Data, 5, 3), m = mael(Data, 5, 4);
1442 13340 : 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 13340 : if (m != 1) m = f;
1447 14593 : for (pu = p; cmpiu(pu,d) <= 0; u++) pu = muliu(pu, up); /* d<pu, pu=p^u */
1448 :
1449 13340 : H = Fl_powers(pmodn, d-1, n); /* H=<p> */
1450 13340 : i_t = get_i_t(n, pmodn); /* i_t[1+i]=k ==> t_i=t_k */
1451 13340 : kT = get_kT(i_t, d); n_T = lg(kT)-1; mitk = kT[n_T];
1452 13340 : if (DEBUGLEVEL == 2) err_printf("kT=%Ps %ld elements\n",kT,n_T);
1453 13340 : if (DEBUGLEVEL >= 6) timer_start(&ti);
1454 13340 : T = galoissubcyclo(utoi(n), utoi(pmodn), 0, 0);
1455 13340 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "galoissubcyclo");
1456 13340 : d0d1 = get_d0_d1(T, gel(fn,1)); /* d0*T_k(x) is in Z[x] */
1457 13340 : Data2 = mkvecn(6, H, GH, i_t, d0d1, kT, mkvecsmalln(5, n, d, f, n_T, mitk));
1458 13340 : vT = get_vT(Data2, 0);
1459 13340 : if (DEBUGLEVEL == 4) err_printf("vT=%Ps\n",vT);
1460 13340 : r = QXV_den_pval(vT, kT, p);
1461 13340 : Rs = ZpX_roots_all(T, p, f, &s);
1462 13340 : if (DEBUGLEVEL >= 2) err_printf("(u,s,r)=(%ld,%ld,%ld)\n",u,s,r);
1463 13340 : 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 13340 : R = (r+u==s) ? Rs : ZV_sort_shallow(ZX_Zp_liftroots(T, Rs, p, s, r+u));
1466 13340 : pr = powiu(p, r); pru = powiu(p, r+u); /* Usually, r=0, s=1, pr=1, pru=p */
1467 13340 : 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 13340 : ulong upr = itou(pr), upru = itou(pru), upu = itou(pu);
1483 63651 : for (k = 1; k <= n_T; k++)
1484 : {
1485 50311 : long itk = kT[k];
1486 50311 : GEN z = r? RgX_muls(gel(vT, itk), upr): gel(vT, itk);
1487 50311 : gel(vT, itk) = RgX_to_Flx(z, upru);
1488 : }
1489 13340 : Data3 = mkvecsmall4(up, upr, upu, upru);
1490 13340 : if (DEBUGLEVEL >= 6) timer_start(&ti);
1491 37442 : for (i = 1; i <= m; i++)
1492 24102 : gel(R,i) = Flx_pol_newton_general(Data2, Data3, vT, itou(gel(R,i)));
1493 13340 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "Flx_pol_newton_general");
1494 : }
1495 13340 : 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 336 : Flx_pol_newton_general_new3(GEN Data, long k)
1502 : {
1503 336 : GEN i_t = gel(Data,5), p = gel(Data,7), pu = gel(Data,8), prs = gel(Data,10);
1504 336 : 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 336 : : Fl_mk_v_t_p3(Data, k);
1507 336 : if (DEBUGLEVEL == 3) err_printf("v_t_p=%Ps\n",v_t_p);
1508 336 : 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 7 : Flx_factcyclo_newton_general_new3(GEN Data)
1514 : {
1515 7 : long i, f = mael(Data, 5, 3), m = mael(Data, 5, 4);
1516 : GEN Data2, pu, pols;
1517 : pari_timer ti;
1518 :
1519 7 : if (m != 1) m = f;
1520 7 : pols = cgetg(1+m, t_VEC);
1521 7 : 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 7 : if (DEBUGLEVEL >= 6) timer_start(&ti);
1525 7 : 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 343 : { for (i = 1; i <= m; i++)
1531 336 : gel(pols, i) = Flx_pol_newton_general_new3(Data2, i);
1532 7 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "Flx_pol_newton_general_new3"); }
1533 7 : return pols;
1534 : }
1535 :
1536 : /* return normalized z(-x) */
1537 : static GEN
1538 65972 : Flx_1st_lift_2(GEN z, ulong p)
1539 : {
1540 65972 : long i, r = lg(z);
1541 65972 : GEN x = cgetg(r, t_VECSMALL);
1542 65972 : x[1] = z[1];
1543 65972 : if (odd(r))
1544 194270 : for (i = 2; i<r; i++) uel(x,i) = odd(i)? Fl_neg(uel(z,i), p) : uel(z,i);
1545 : else
1546 231633 : for (i = 2; i<r; i++) uel(x,i) = odd(i)? uel(z,i): Fl_neg(uel(z,i), p);
1547 65972 : return x;
1548 : }
1549 :
1550 : /* el does not divides e.
1551 : * construct Phi_{e*el}(x) from Phi_e(x). */
1552 : static GEN
1553 42709 : Flx_1st_lift(GEN z, ulong p, ulong e, ulong el, GEN vP)
1554 : {
1555 42709 : GEN z2, z3, P = gel(vP, e);
1556 42709 : if (!gel(vP, e)) P = gel(vP, e) = Flx_polcyclo(e, p);
1557 42709 : z2 = Flx_inflate(z, el);
1558 42709 : z3 = Flx_normalize(Flx_gcd(P, z2, p), p);
1559 42709 : return Flx_div(z2, z3, p);
1560 : }
1561 :
1562 : static GEN
1563 158258 : Flx_lift(GEN z, ulong p, ulong e, ulong el, ulong r, ulong s, GEN vP)
1564 : {
1565 158258 : if (s == 0)
1566 : {
1567 108681 : z = (el==2) ? Flx_1st_lift_2(z, p) : Flx_1st_lift(z, p, e, el, vP);
1568 108681 : if (r >= 2) z = Flx_inflate(z, upowuu(el, r-1));
1569 : }
1570 : else
1571 49577 : z = Flx_inflate(z, upowuu(el, r-s));
1572 158258 : return z;
1573 : }
1574 :
1575 : /* e is the conductor of the splitting field of p in Q(zeta_n) */
1576 : static GEN
1577 140965 : Flx_conductor_lift(GEN z, ulong p, GEN fn, ulong e, GEN vP)
1578 : {
1579 140965 : GEN EL = gel(fn, 1), En = gel(fn, 2);
1580 140965 : long i, r = lg(EL), new_e = e;
1581 :
1582 432170 : for (i = 1; i < r; i++)
1583 : {
1584 291205 : long el = EL[i], en = En[i], ee = u_lval(e, el);
1585 291205 : if (ee < en)
1586 : {
1587 158258 : z = Flx_lift(z, p, new_e, el, en, ee, vP);
1588 158258 : new_e *= upowuu(el, en-ee);
1589 : }
1590 : }
1591 140965 : return z;
1592 : }
1593 :
1594 : /* R0 is mod p^u, d < p^u */
1595 : static GEN
1596 57409 : Flx_pol_newton(long j, GEN R0, GEN E, GEN D3, long d, long f, ulong p)
1597 : {
1598 57409 : ulong u = D3[3];
1599 57409 : GEN R = cgetg(f+1, t_VECSMALL);
1600 : long i;
1601 233473 : for (i = 1; i <= f; i++) R[i] = R0[1+(i+j)%f];
1602 57409 : 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 31564 : Flx_factcyclo_newton_pre(GEN Data, GEN N, ulong p, ulong m)
1608 : {
1609 31564 : GEN T = gel(Data, 1), F = gel(Data, 2), Rs = gel(Data, 3), D4 = gel(Data, 4);
1610 31564 : long d = D4[1], nf = D4[2], g = D4[3], r = D4[4], s = D4[5], u = D4[6];
1611 31564 : GEN pols, E, R, p0 = utoi(p), Data3;
1612 31564 : ulong i, n = N[1], pmodn = p%n;
1613 : pari_timer ti;
1614 :
1615 31564 : if (m != 1) m = nf;
1616 31564 : pols = cgetg(1+m, t_VEC);
1617 31564 : E = set_E(pmodn, n, d, nf, g);
1618 31564 : R = set_R(T, F, Rs, p0, nf, r, s, u);
1619 31564 : R = ZV_to_nv(R);
1620 31564 : Data3 = mkvecsmall3(r, s, u);
1621 31564 : if (DEBUGLEVEL >= 6) timer_start(&ti);
1622 88973 : for (i = 1; i <= m; i++) gel(pols, i) = Flx_pol_newton(i, R,E,Data3, d,nf,p);
1623 31564 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "Flx_pol_newton");
1624 31564 : 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 gerepilecopy(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 43996 : Flx_factcyclo_gen(GEN GH, ulong n, ulong p, ulong m)
1675 : {
1676 43996 : GEN fu = factoru(n), fa = zm_to_ZM(fu), p0 = utoi(p);
1677 : GEN T, X, pd_n, v, pols;
1678 43996 : ulong i, j, pmodn = p%n, phin = eulerphiu_fact(fu);
1679 43996 : ulong d = Fl_order(pmodn, phin, n), nf = phin/d;
1680 : pari_timer ti;
1681 :
1682 43996 : if (m != 1) m = nf;
1683 43996 : 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 43996 : pols = cgetg(1+m, t_VEC);
1690 43996 : v = cgetg(1+d, t_VEC);
1691 43996 : pd_n = diviuexact(subis(powiu(p0, d), 1), n); /* (p^d-1)/n */
1692 43996 : T = ZX_to_Flx(init_Fq(p0, d, evalvarn(1)), p);
1693 43996 : random_Flx(degpol(T), T[1], p); /* skip 1st one */
1694 43996 : if (DEBUGLEVEL >= 6) timer_start(&ti);
1695 82769 : do X = Flxq_pow(random_Flx(degpol(T), T[1], p), pd_n, T, p);
1696 82769 : while (lg(X)<3 || equaliu(Flxq_order(X, fa, T, p), n)==0); /* find zeta_n */
1697 43996 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "find X");
1698 :
1699 43996 : if (m == 1)
1700 : {
1701 102528 : for (j = 1; j <= d; j++)
1702 : {
1703 80481 : gel(v, j) = polx_FlxX(0, 1);
1704 80481 : gmael(v, j, 2) = Flx_neg(X, p);
1705 80481 : if (j < d) X = Flxq_powu(X, p, T, p);
1706 : }
1707 22047 : gel(pols, 1) = FlxX_to_Flx(FlxqXV_prod(v, T, p));
1708 : }
1709 : else
1710 : {
1711 : GEN vX, vp;
1712 21949 : if (DEBUGLEVEL >= 6) timer_start(&ti);
1713 21949 : vX = Flxq_powers(X, n, T, p)+1;
1714 21949 : vp = Fl_powers(pmodn, d, n);
1715 190844 : for (i = 1; i <= m; i++)
1716 : {
1717 757867 : for (j = 1; j <= d; j++)
1718 : {
1719 588972 : gel(v, j) = polx_FlxX(0, 1);
1720 588972 : gmael(v, j, 2) = Flx_neg(gel(vX, Fl_mul(GH[i], vp[j], n)), p);
1721 : }
1722 168895 : gel(pols, i) = FlxX_to_Flx(FlxqXV_prod(v, T, p));
1723 : }
1724 21949 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "FlxqXV_prod");
1725 : }
1726 43996 : return pols;
1727 : }
1728 :
1729 : static int
1730 1734612 : 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 34186 : Flx_split(ulong n, ulong p, ulong m)
1735 : {
1736 34186 : ulong i, j, z = rootsof1_Fl(n, p);
1737 : GEN v, C, vx;
1738 34186 : if (m == 1) return mkvec(mkvecsmall3(0, Fl_neg(z,p), 1));
1739 17016 : v = cgetg(m+1, t_VEC); C = coprimes_zv(n); vx = Fl_powers(z, n-1, p);
1740 223026 : for (i = j = 1; i <= n; i++)
1741 206010 : if (C[i]) gel(v, j++) = mkvecsmall3(0, Fl_neg(vx[i+1], p), 1);
1742 17016 : return gen_sort(v, (void*)cmpGuGu, &gen_cmp_RgX);
1743 : }
1744 :
1745 : /* d==1 or f==1 occurs */
1746 : static GEN
1747 31564 : Flx_factcyclo_prime_power_i(long el, long e, long p, long m)
1748 : {
1749 31564 : GEN p0 = utoipos(p), z = set_e0_e1(el, e, p0), v;
1750 31564 : long n = z[1], e0 = z[2], e1 = z[3], phin = z[4], g = z[5];
1751 31564 : long d = z[6], f = z[7]; /* d and f for n=el^e0 */
1752 :
1753 31564 : if (f == 1) v = mkvec(Flx_polcyclo(n, p));
1754 31564 : else if (d == 1) v = Flx_split(n, p, (m==1)?1:f);
1755 31564 : else if (el == 2) v = Flx_factcyclo_gen(NULL, n, p, m);/* d==2 in this case */
1756 31564 : else if (!use_newton(d, f)) v = Flx_factcyclo_gen(NULL, n, p, m);
1757 : else
1758 : {
1759 31564 : GEN N = mkvecsmall5(n, el, e0, phin, g);
1760 31564 : v = FpX_factcyclo_newton_power(N, p0, m, 1);
1761 31564 : 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 31564 : 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 31564 : 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 gerepilecopy(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 88907 : Flx_factcyclo_just_conductor(ulong n, ulong p, ulong m)
1803 : {
1804 : GEN Data, fn;
1805 88907 : ulong action = FpX_factcyclo_just_conductor_init(&Data, n, utoipos(p), m);
1806 88907 : fn = gel(Data,3);
1807 88907 : if (action & GENERAL)
1808 43996 : return Flx_factcyclo_gen(gel(Data,2), n, p, m);
1809 44911 : else if (action & NEWTON_POWER)
1810 31564 : return Flx_factcyclo_prime_power_i(ucoeff(fn,1,1), ucoeff(fn,1,2), p, m);
1811 13347 : else if (action & NEWTON_GENERAL)
1812 13340 : return Flx_factcyclo_newton_general(Data);
1813 7 : else if (action & NEWTON_GENERAL_NEW)
1814 7 : 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 156721 : Flx_factcyclo_i(ulong n, ulong p, ulong fl)
1821 : {
1822 156721 : GEN fn = factoru(n), z;
1823 156721 : ulong phin = eulerphiu_fact(fn), pmodn = p%n;
1824 156721 : ulong d = Fl_order(pmodn, phin, n), f = phin/d, fK;
1825 :
1826 156721 : if (DEBUGLEVEL >= 1) header(fn, n, d, f, utoi(p));
1827 156721 : if (f == 1) { z = Flx_polcyclo(n, p); return fl? z: mkvec(z); }
1828 123093 : if (d == 1) /* p=1 (mod n), zeta_n in Z_p */
1829 9488 : { z = Flx_split(n, p, fl? 1: f); return fl? gel(z,1): z; }
1830 113605 : fK = znstar_conductor(znstar_generate(n, mkvecsmall(pmodn)));
1831 113605 : if (fK != n && p % fK == 1)
1832 24698 : z = Flx_split(fK, p, fl? 1: f);
1833 : else
1834 88907 : z = Flx_factcyclo_just_conductor(fK, p, fl? 1: f);
1835 113605 : if (n > fK)
1836 : {
1837 58055 : GEN vP = const_vec(n, NULL);
1838 58055 : long i, l = fl? 2: lg(z);
1839 199020 : for (i = 1; i < l; i++)
1840 140965 : gel(z, i) = Flx_conductor_lift(gel(z, i), p, fn, fK, vP);
1841 : }
1842 113605 : return fl? gel(z,1): gen_sort(z,(void*)cmpGuGu, &gen_cmp_RgX);
1843 : }
1844 :
1845 : GEN
1846 308 : Flx_factcyclo(ulong n, ulong p, ulong m)
1847 308 : { pari_sp av = avma; return gerepilecopy(av, Flx_factcyclo_i(n, p, m)); }
1848 :
1849 : GEN
1850 167027 : factormodcyclo(long n, GEN p, long fl, long v)
1851 : {
1852 167027 : const char *fun = "factormodcyclo";
1853 167027 : pari_sp av = avma;
1854 : long i, l;
1855 : GEN z;
1856 167027 : if (v < 0) v = 0;
1857 167027 : if (fl < 0 || fl > 1) pari_err_FLAG(fun);
1858 167027 : if (n <= 0) pari_err_DOMAIN(fun, "n", "<=", gen_0, stoi(n));
1859 167027 : if (typ(p) != t_INT) pari_err_TYPE(fun, p);
1860 167027 : if (dvdui(n, p)) pari_err_COPRIME(fun, stoi(n), p);
1861 167027 : if (fl)
1862 : {
1863 83503 : if (lgefint(p) == 3)
1864 78203 : z = Flx_to_ZX(Flx_factcyclo_i(n, p[2], 1));
1865 : else
1866 5300 : z = FpX_factcyclo_i(n, p, 1);
1867 83503 : setvarn(z, v);
1868 83503 : return gerepileupto(av, FpX_to_mod(z, p));
1869 : }
1870 : else
1871 : {
1872 83524 : if (lgefint(p) == 3)
1873 78210 : z = FlxC_to_ZXC(Flx_factcyclo_i(n, p[2], 0));
1874 : else
1875 5314 : z = FpX_factcyclo_i(n, p, 0);
1876 465409 : l = lg(z); for (i = 1; i < l; i++) setvarn(gel(z, i), v);
1877 83524 : return gerepileupto(av, FpXC_to_mod(z, p));
1878 : }
1879 : }
|
