Line data Source code
1 : /* Copyright (C) 2000 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; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : #define DEBUGLEVEL DEBUGLEVEL_subcyclo
19 :
20 : /*************************************************************************/
21 : /** **/
22 : /** Routines for handling subgroups of (Z/nZ)^* **/
23 : /** without requiring discrete logarithms. **/
24 : /** **/
25 : /*************************************************************************/
26 : /* Subgroups are [gen,ord,bits] where
27 : * gen is a vecsmall of generators
28 : * ord is theirs relative orders
29 : * bits is a bit vector of the elements, of length(n). */
30 :
31 : /*The algorithm is similar to testpermutation*/
32 : static void
33 1157798 : znstar_partial_coset_func(long n, GEN H, void (*func)(void *data,long c)
34 : , void *data, long d, long c)
35 : {
36 : GEN gen, ord, cache;
37 : long i, j, card;
38 :
39 1157798 : if (!d) { (*func)(data,c); return; }
40 :
41 938769 : cache = const_vecsmall(d,c);
42 938769 : (*func)(data,c); /* AFTER cache: may contain gerepileupto statement */
43 938769 : gen = gel(H,1);
44 938769 : ord = gel(H,2);
45 990927 : card = ord[1]; for (i = 2; i <= d; i++) card *= ord[i];
46 70745045 : for(i=1; i<card; i++)
47 : {
48 69806276 : long k, m = i;
49 70297700 : for(j=1; j<d && m%ord[j]==0 ;j++) m /= ord[j];
50 69806276 : cache[j] = Fl_mul(cache[j],gen[j],n);
51 70297700 : for (k=1; k<j; k++) cache[k] = cache[j];
52 69806276 : (*func)(data, cache[j]);
53 : }
54 : }
55 :
56 : static void
57 237180 : znstar_coset_func(long n, GEN H, void (*func)(void *data,long c)
58 : , void *data, long c)
59 237180 : { znstar_partial_coset_func(n, H, func,data, lg(gel(H,1))-1, c); }
60 :
61 : /* Add the element of the bitvec of the coset c modulo the subgroup of H
62 : * generated by the first d generators to the bitvec bits.*/
63 :
64 : static void
65 920618 : znstar_partial_coset_bits_inplace(long n, GEN H, GEN bits, long d, long c)
66 : {
67 920618 : pari_sp av = avma;
68 920618 : znstar_partial_coset_func(n,H, (void (*)(void *,long)) &F2v_set,
69 : (void *) bits, d, c);
70 920618 : set_avma(av);
71 920618 : }
72 :
73 : static void
74 456208 : znstar_coset_bits_inplace(long n, GEN H, GEN bits, long c)
75 456208 : { znstar_partial_coset_bits_inplace(n, H, bits, lg(gel(H,1))-1, c); }
76 :
77 : static GEN
78 464410 : znstar_partial_coset_bits(long n, GEN H, long d, long c)
79 : {
80 464410 : GEN bits = zero_F2v(n);
81 464410 : znstar_partial_coset_bits_inplace(n,H,bits,d,c);
82 464410 : return bits;
83 : }
84 :
85 : /* Compute the bitvec of the elements of the subgroup of H generated by the
86 : * first d generators.*/
87 : static GEN
88 464410 : znstar_partial_bits(long n, GEN H, long d)
89 464410 : { return znstar_partial_coset_bits(n, H, d, 1); }
90 :
91 : /* Compute the bitvec of the elements of H. */
92 : GEN
93 0 : znstar_bits(long n, GEN H)
94 0 : { return znstar_partial_bits(n,H,lg(gel(H,1))-1); }
95 :
96 : /* Compute the subgroup of (Z/nZ)^* generated by the elements of
97 : * the vecsmall V */
98 : GEN
99 218669 : znstar_generate(long n, GEN V)
100 : {
101 218669 : pari_sp av = avma;
102 218669 : GEN gen = cgetg(lg(V),t_VECSMALL);
103 218669 : GEN ord = cgetg(lg(V),t_VECSMALL), res = mkvec2(gen,ord);
104 218669 : GEN bits = znstar_partial_bits(n,NULL,0);
105 218669 : long i, r = 0;
106 1586638 : for(i=1; i<lg(V); i++)
107 : {
108 1367969 : ulong v = uel(V,i), g = v;
109 1367969 : long o = 0;
110 32793015 : while (!F2v_coeff(bits, (long)g)) { g = Fl_mul(g, v, (ulong)n); o++; }
111 1367969 : if (!o) continue;
112 245741 : r++;
113 245741 : gen[r] = v;
114 245741 : ord[r] = o+1;
115 245741 : cgiv(bits); bits = znstar_partial_bits(n,res,r);
116 : }
117 218669 : setlg(gen,r+1);
118 218669 : setlg(ord,r+1); return gerepilecopy(av, mkvec3(gen,ord,bits));
119 : }
120 :
121 : static ulong
122 140741 : znstar_order(GEN H) { return zv_prod(gel(H,2)); }
123 :
124 : /* Return the lists of element of H.
125 : * This can be implemented with znstar_coset_func instead. */
126 : GEN
127 3240 : znstar_elts(long n, GEN H)
128 : {
129 3240 : long card = znstar_order(H);
130 3240 : GEN gen = gel(H,1), ord = gel(H,2);
131 3240 : GEN sg = cgetg(1 + card, t_VECSMALL);
132 : long k, j, l;
133 3240 : sg[1] = 1;
134 4968 : for (j = 1, l = 1; j < lg(gen); j++)
135 : {
136 1728 : long c = l * (ord[j]-1);
137 3924 : for (k = 1; k <= c; k++) sg[++l] = Fl_mul(sg[k], gen[j], n);
138 : }
139 3240 : vecsmall_sort(sg); return sg;
140 : }
141 :
142 : /* Take a znstar H and n dividing the modulus of H.
143 : * Output H reduced to modulus n */
144 : GEN
145 156 : znstar_reduce_modulus(GEN H, long n)
146 : {
147 156 : pari_sp ltop=avma;
148 156 : GEN gen=cgetg(lgcols(H),t_VECSMALL);
149 : long i;
150 546 : for(i=1; i < lg(gen); i++)
151 390 : gen[i] = mael(H,1,i)%n;
152 156 : return gerepileupto(ltop, znstar_generate(n,gen));
153 : }
154 :
155 : /* Compute conductor of H, bits = H[3] */
156 : long
157 147630 : znstar_conductor_bits(GEN bits)
158 : {
159 147630 : pari_sp av = avma;
160 147630 : long i, f = 1, cnd0 = bits[1];
161 147630 : GEN F = factoru(cnd0), P = gel(F,1), E = gel(F,2);
162 397058 : for (i = lg(P)-1; i > 0; i--)
163 : {
164 249428 : long p = P[i], e = E[i], cnd = cnd0;
165 281528 : for ( ; e >= 2; e--)
166 : {
167 71930 : long q = cnd / p;
168 71930 : if (!F2v_coeff(bits, 1 + q)) break;
169 32100 : cnd = q;
170 : }
171 249428 : if (e == 1)
172 : {
173 209598 : if (p == 2) e = 0;
174 : else
175 : {
176 184974 : long h, g = pgener_Fl(p), q = cnd / p;
177 184974 : h = Fl_mul(g-1, Fl_inv(q % p, p), p); /* 1+h*q = g (mod p) */
178 184974 : if (F2v_coeff(bits, 1 + h*q)) e = 0;
179 : }
180 : }
181 249428 : if (e) f *= upowuu(p, e);
182 : }
183 147630 : return gc_long(av,f);
184 : }
185 : long
186 147066 : znstar_conductor(GEN H) { return znstar_conductor_bits(gel(H,3)); }
187 :
188 : /* Compute the orbits of a subgroups of Z/nZ given by a generator
189 : * or a set of generators given as a vector.
190 : */
191 : GEN
192 94787 : znstar_cosets(long n, long phi_n, GEN H)
193 : {
194 94787 : long k, c = 0, card = znstar_order(H), index = phi_n/card;
195 94787 : GEN cosets = cgetg(index+1,t_VECSMALL);
196 94787 : pari_sp ltop = avma;
197 94787 : GEN bits = zero_F2v(n);
198 550959 : for (k = 1; k <= index; k++)
199 : {
200 1177062 : for (c++ ; F2v_coeff(bits,c) || ugcd(c,n)!=1; c++);
201 456172 : cosets[k]=c;
202 456172 : znstar_coset_bits_inplace(n, H, bits, c);
203 : }
204 94787 : set_avma(ltop); return cosets;
205 : }
206 :
207 : static GEN
208 6 : znstar_quotient(long n, long phi_n, GEN H, GEN R, ulong l)
209 : {
210 6 : long i, j, k, c = 0, card = znstar_order(H), index = phi_n/card;
211 6 : GEN cosets = cgetg(index+1,t_VECSMALL), mult = cgetg(index+1, t_VEC);
212 6 : GEN bits = zero_F2v(n), vbits= zero_F2m_copy(n,index);
213 42 : for (k = 1; k <= index; k++)
214 : {
215 102 : for (c++ ; F2v_coeff(bits,c) || ugcd(c,n)!=1; c++);
216 36 : cosets[k]=c;
217 36 : znstar_coset_bits_inplace(n, H, gel(vbits,k), c);
218 36 : F2v_add_inplace(bits, gel(vbits,k));
219 : }
220 :
221 42 : for (k = 1; k <= index; k++)
222 : {
223 36 : GEN v = cgetg(index+1, t_VECSMALL);
224 252 : for (j = 1; j <= index; j++)
225 : {
226 216 : long s = Fl_mul(cosets[k],cosets[j],n);
227 756 : for (i = 1; i <= index; i++)
228 756 : if (F2v_coeff(gel(vbits,i),s)) break;
229 216 : v[j] = R[i];
230 : }
231 36 : gel(mult,k) = v;
232 : }
233 6 : return Flv_Flm_polint(R, mult, l, 0);
234 : }
235 :
236 : /* return n s.t. x^n in H */
237 : static ulong
238 1164228 : order_H_x(GEN H, ulong x, GEN D, ulong *xn)
239 : {
240 1164228 : ulong i, l = lg(D), f = H[1], y = x; /* = x^D[1] */
241 1164228 : *xn = x; if (F2v_coeff(H, y)) return D[1];
242 22368108 : for (i = 2; i < l; i++)
243 : { /* TODO: could cache the x^delta[i] incrementally */
244 22368108 : y = Fl_mul(y, Fl_powu(x, D[i]-D[i-1], f), f);
245 22368108 : if (F2v_coeff(H, y))
246 : {
247 1164228 : if (xn) *xn = y;
248 1164228 : return D[i];
249 : }
250 : }
251 0 : pari_err_BUG("znsubgroupgenerators [order_H_x]");
252 : return 0;/*LCOV_EXCL_LINE*/
253 : }
254 : /* If H2 = (0), return 0. Else find an x in H2 of maximal order o in H2/H1;
255 : * if flag is set, make sure that x has also order o in (Z/fZ)^* */
256 : static ulong
257 60 : max_order_ele(GEN H1, GEN H2, GEN D, ulong *po, long flag)
258 : {
259 60 : ulong x, h, n = F2v_hamming(H2), f = H1[1], O = 0, X = 0, XO = 0;
260 60 : if (!n) return 0;
261 34869072 : for (x = 2; x < f; x++) if (F2v_coeff(H2, x))
262 : {
263 1164228 : ulong xo, o = order_H_x(H1, x, D, &xo);
264 1164228 : if (o > O) { O = o; X = x; XO = xo; if (o == n) break; }
265 : }
266 : /* X of maximal order O in H2/H1 */
267 48 : *po = O; if (!flag || XO == 1) return X;
268 : /* find h in H1 s.t. (Xh)^o=1 */
269 18 : x = Fl_inv(XO, f);
270 1592238 : for (h = 1;; h++) /* stops for h < f */
271 1592238 : if (F2v_coeff(H1, h) && x == Fl_powu(h, O, f)) break;
272 18 : return Fl_mul(X, h, f);
273 : }
274 : /* x not in H. Replace H in place by the subgroup generated by H and x,
275 : * which has order o > 1 in G/H */
276 : static void
277 48 : enlarge_H(GEN H, ulong x, ulong o)
278 : {
279 48 : pari_sp av = avma;
280 48 : ulong i, j, l = lg(H), f = H[1];
281 48 : GEN H1 = vecsmall_copy(H), y = Fl_powers(x, o-1, f);
282 34869120 : for (i = 1; i < f; i++)
283 34869072 : if (F2v_coeff(H, i))
284 650160 : for (j = 2; j <= o; j++) F2v_set(H1, Fl_mul(i, y[j], f));
285 635704 : for (i = 2; i < l; i++) H[i] = H1[i];
286 48 : set_avma(av);
287 48 : }
288 : /* H F2v subgroup of (Z/fZ)^*: a in H <=> H[a] = 1
289 : * Return generators. If flag is set, they generate the cyclic components */
290 : GEN
291 18 : znsubgroupgenerators(GEN H, long flag)
292 : {
293 18 : pari_sp av = avma;
294 : ulong f, g, o;
295 18 : GEN H1, H2, z = const_vecsmall(0,0), D;
296 18 : switch(typ(H))
297 : {
298 12 : case t_VECSMALL: H = Flv_to_F2v(H); break;
299 6 : case t_VEC: if (RgV_is_ZV(H)) { H = ZV_to_F2v(H); break; }
300 6 : default: pari_err_TYPE("znsubgroupgenerators", H);
301 : }
302 12 : f = H[1]; z = cgetg(1, t_VECSMALL);
303 12 : D = divisorsu(F2v_hamming(H));
304 12 : H1 = zero_F2v(f); F2v_set(H1, 1);
305 12 : H2 = H;
306 60 : while ((g = max_order_ele(H1, H2, D, &o, flag)))
307 : {
308 48 : z = vecsmall_append(z, g); enlarge_H(H1, g, o);
309 48 : F2v_negimply_inplace(H2, H1); /* H2 <- H2 - H1 */
310 : }
311 12 : return gerepileupto(av, zv_to_ZV(z));
312 : }
313 :
314 : /*************************************************************************/
315 : /** **/
316 : /** znstar/HNF interface **/
317 : /** **/
318 : /*************************************************************************/
319 :
320 : static long
321 2538 : mod_to_small(GEN x)
322 2538 : { return itos(typ(x) == t_INTMOD ? gel(x,2): x); }
323 :
324 : static GEN
325 1872 : vecmod_to_vecsmall(GEN x)
326 4410 : { pari_APPLY_long(mod_to_small(gel(x,i))) }
327 :
328 : /* Convert a true znstar output by znstar to a `small znstar' */
329 : GEN
330 1872 : znstar_small(GEN zn)
331 : {
332 1872 : GEN Z = cgetg(4,t_VEC);
333 1872 : gel(Z,1) = icopy(gmael3(zn,3,1,1));
334 1872 : gel(Z,2) = vec_to_vecsmall(gel(zn,2));
335 1872 : gel(Z,3) = vecmod_to_vecsmall(gel(zn,3)); return Z;
336 : }
337 :
338 : /* Compute generators for the subgroup of (Z/nZ)* given in HNF. */
339 : GEN
340 3600 : znstar_hnf_generators(GEN Z, GEN M)
341 : {
342 3600 : long j, h, l = lg(M);
343 3600 : GEN gen = cgetg(l, t_VECSMALL);
344 3600 : pari_sp ltop = avma;
345 3600 : GEN zgen = gel(Z,3);
346 3600 : ulong n = itou(gel(Z,1));
347 7866 : for (j = 1; j < l; j++)
348 : {
349 4266 : GEN Mj = gel(M,j);
350 4266 : gen[j] = 1;
351 10560 : for (h = 1; h < l; h++)
352 : {
353 6294 : ulong u = itou(gel(Mj,h));
354 6294 : if (!u) continue;
355 4632 : gen[j] = Fl_mul(uel(gen,j), Fl_powu(uel(zgen,h), u, n), n);
356 : }
357 : }
358 3600 : set_avma(ltop); return gen;
359 : }
360 :
361 : GEN
362 3240 : znstar_hnf(GEN Z, GEN M)
363 3240 : { return znstar_generate(itos(gel(Z,1)),znstar_hnf_generators(Z,M)); }
364 :
365 : GEN
366 3240 : znstar_hnf_elts(GEN Z, GEN H)
367 : {
368 3240 : pari_sp ltop = avma;
369 3240 : GEN G = znstar_hnf(Z,H);
370 3240 : long n = itos(gel(Z,1));
371 3240 : return gerepileupto(ltop, znstar_elts(n,G));
372 : }
373 :
374 : /*************************************************************************/
375 : /** **/
376 : /** polsubcyclo **/
377 : /** **/
378 : /*************************************************************************/
379 :
380 : static GEN
381 42774 : gscycloconductor(GEN g, long n, long flag)
382 : {
383 42774 : if (flag==2) retmkvec2(gcopy(g), stoi(n));
384 42768 : return g;
385 : }
386 :
387 : static long
388 1146804 : lift_check_modulus(GEN H, long n)
389 : {
390 : long h;
391 1146804 : switch(typ(H))
392 : {
393 90 : case t_INTMOD:
394 90 : if (!equalsi(n, gel(H,1)))
395 6 : pari_err_MODULUS("galoissubcyclo", stoi(n), gel(H,1));
396 84 : H = gel(H,2); /* fall through */
397 1146798 : case t_INT:
398 1146798 : h = smodis(H,n);
399 1146798 : if (ugcd(h,n) != 1) pari_err_COPRIME("galoissubcyclo", H,stoi(n));
400 1146792 : return h ? h: 1;
401 : }
402 0 : pari_err_TYPE("galoissubcyclo [subgroup]", H);
403 : return 0;/*LCOV_EXCL_LINE*/
404 : }
405 :
406 : /* Compute z^ex using the baby-step/giant-step table powz
407 : * with only one multiply.
408 : * In the modular case, the result is not reduced. */
409 : static GEN
410 4862464 : polsubcyclo_powz(GEN powz, long ex)
411 : {
412 4862464 : long m = lg(gel(powz,1))-1, q = ex/m, r = ex%m; /*ex=m*q+r*/
413 4862464 : GEN g = gmael(powz,1,r+1), G = gmael(powz,2,q+1);
414 4862464 : return (lg(powz)==4)? mulreal(g,G): gmul(g,G);
415 : }
416 :
417 : static GEN
418 69710 : polsubcyclo_complex_bound(pari_sp av, GEN V, long prec)
419 : {
420 69710 : GEN pol = real_i(roots_to_pol(V,0));
421 69710 : return gerepileuptoint(av, ceil_safe(gsupnorm(pol,prec)));
422 : }
423 :
424 : /* Newton sums mod le. if le==NULL, works with complex instead */
425 : static GEN
426 54004 : polsubcyclo_cyclic(long n, long d, long m ,long z, long g, GEN powz, GEN le)
427 : {
428 54004 : GEN V = cgetg(d+1,t_VEC);
429 54004 : ulong base = 1;
430 : long i,k;
431 : pari_timer ti;
432 54004 : if (DEBUGLEVEL >= 6) timer_start(&ti);
433 262424 : for (i=1; i<=d; i++, base = Fl_mul(base,z,n))
434 : {
435 208420 : pari_sp av = avma;
436 208420 : long ex = base;
437 208420 : GEN s = gen_0;
438 971884 : for (k=0; k<m; k++, ex = Fl_mul(ex,g,n))
439 : {
440 763464 : s = gadd(s, polsubcyclo_powz(powz,ex));
441 763464 : if ((k&0xff)==0) s = gerepileupto(av,s);
442 : }
443 208420 : if (le) s = modii(s, le);
444 208420 : gel(V,i) = gerepileupto(av, s);
445 : }
446 54004 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "polsubcyclo_cyclic");
447 54004 : return V;
448 : }
449 :
450 : struct _subcyclo_orbits_u
451 : {
452 : GEN bab, gig;
453 : ulong l;
454 : ulong s;
455 : long m;
456 : };
457 :
458 : static void
459 0 : _Fl_subcyclo_orbits(void *E, long k)
460 : {
461 0 : struct _subcyclo_orbits_u *D = (struct _subcyclo_orbits_u *) E;
462 0 : ulong l = D->l;
463 0 : long q = k/D->m, r = k%D->m; /*k=m*q+r*/
464 0 : ulong g = uel(D->bab,r+1), G = uel(D->gig,q+1);
465 0 : D->s = Fl_add(D->s, Fl_mul(g, G, l), l);
466 0 : }
467 :
468 : /* Newton sums mod le. if le==NULL, works with complex instead */
469 : static GEN
470 0 : Fl_polsubcyclo_orbits(long n, GEN H, GEN O, ulong z, ulong l)
471 : {
472 0 : long i, d = lg(O);
473 0 : GEN V = cgetg(d,t_VECSMALL);
474 : struct _subcyclo_orbits_u D;
475 0 : long m = 1+usqrt(n);
476 0 : D.l = l;
477 0 : D.m = m;
478 0 : D.bab = Fl_powers(z, m, l);
479 0 : D.gig = Fl_powers(uel(D.bab,m+1), m-1, l);
480 0 : for(i=1; i<d; i++)
481 : {
482 0 : D.s = 0;
483 0 : znstar_coset_func(n, H, _Fl_subcyclo_orbits, (void *) &D, O[i]);
484 0 : uel(V,i) = D.s;
485 : }
486 0 : return V;
487 : }
488 :
489 : struct _subcyclo_orbits_s
490 : {
491 : GEN powz;
492 : GEN *s;
493 : ulong count;
494 : pari_sp ltop;
495 : };
496 :
497 : static void
498 4099000 : _subcyclo_orbits(struct _subcyclo_orbits_s *data, long k)
499 : {
500 4099000 : GEN powz = data->powz;
501 4099000 : GEN *s = data->s;
502 :
503 4099000 : if (!data->count) data->ltop = avma;
504 4099000 : *s = gadd(*s, polsubcyclo_powz(powz,k));
505 4099000 : data->count++;
506 4099000 : if ((data->count & 0xffUL) == 0) *s = gerepileupto(data->ltop, *s);
507 4099000 : }
508 :
509 : /* Newton sums mod le. if le==NULL, works with complex instead */
510 : static GEN
511 85416 : polsubcyclo_orbits(long n, GEN H, GEN O, GEN powz, GEN le)
512 : {
513 85416 : long i, d = lg(O);
514 85416 : GEN V = cgetg(d,t_VEC);
515 : struct _subcyclo_orbits_s data;
516 85416 : long lle = le?lg(le)*2+1: 2*lg(gmael(powz,1,2))+3;/*dvmdii uses lx+ly space*/
517 85416 : data.powz = powz;
518 322596 : for(i=1; i<d; i++)
519 : {
520 237180 : GEN s = gen_0;
521 237180 : pari_sp av = avma;
522 237180 : (void)new_chunk(lle);
523 237180 : data.count = 0;
524 237180 : data.s = &s;
525 237180 : znstar_coset_func(n, H, (void (*)(void *,long)) _subcyclo_orbits,
526 237180 : (void *) &data, O[i]);
527 237180 : set_avma(av); /* HACK */
528 237180 : gel(V,i) = le? modii(s,le): gcopy(s);
529 : }
530 85416 : return V;
531 : }
532 :
533 : static GEN
534 70460 : polsubcyclo_start(long n, long d, long o, long e, GEN borne, long *ptr_val,long *ptr_l)
535 : {
536 : pari_sp av;
537 : GEN le, z, gl;
538 : long i, l, val;
539 70460 : l = e*n+1;
540 267758 : while(!uisprime(l)) { l += n; e++; }
541 70460 : if (DEBUGLEVEL >= 4) err_printf("Subcyclo: prime l=%ld\n",l);
542 70460 : gl = utoipos(l); av = avma;
543 70460 : if (!borne)
544 : { /* Use vecmax(Vec((x+o)^d)) = max{binomial(d,i)*o^i ;1<=i<=d} */
545 0 : i = d-(1+d)/(1+o);
546 0 : borne = mulii(binomial(utoipos(d),i),powuu(o,i));
547 : }
548 70460 : if (DEBUGLEVEL >= 4) err_printf("Subcyclo: bound=2^%ld\n",expi(borne));
549 70460 : val = logint(shifti(borne,2), gl) + 1;
550 70460 : set_avma(av);
551 70460 : if (DEBUGLEVEL >= 4) err_printf("Subcyclo: val=%ld\n",val);
552 70460 : le = powiu(gl,val);
553 70460 : z = utoipos( Fl_powu(pgener_Fl(l), e, l) );
554 70460 : z = Zp_sqrtnlift(gen_1,utoipos(n),z,gl,val);
555 70460 : *ptr_val = val;
556 70460 : *ptr_l = l;
557 70460 : return gmodulo(z,le);
558 : }
559 :
560 : /*Fill in the powz table:
561 : * powz[1]: baby-step
562 : * powz[2]: giant-step
563 : * powz[3] exists only if the field is real (value is ignored). */
564 : static GEN
565 69710 : polsubcyclo_complex_roots(long n, long real, long prec)
566 : {
567 69710 : long i, m = (long)(1+sqrt((double) n));
568 69710 : GEN bab, gig, powz = cgetg(real?4:3, t_VEC);
569 :
570 69710 : bab = cgetg(m+1,t_VEC);
571 69710 : gel(bab,1) = gen_1;
572 69710 : gel(bab,2) = rootsof1u_cx(n, prec); /* = e_n(1) */
573 322818 : for (i=3; i<=m; i++) gel(bab,i) = gmul(gel(bab,2),gel(bab,i-1));
574 69710 : gig = cgetg(m+1,t_VEC);
575 69710 : gel(gig,1) = gen_1;
576 69710 : gel(gig,2) = gmul(gel(bab,2),gel(bab,m));;
577 322818 : for (i=3; i<=m; i++) gel(gig,i) = gmul(gel(gig,2),gel(gig,i-1));
578 69710 : gel(powz,1) = bab;
579 69710 : gel(powz,2) = gig;
580 69710 : if (real) gel(powz,3) = gen_0;
581 69710 : return powz;
582 : }
583 :
584 : static GEN
585 575926 : muliimod_sz(GEN x, GEN y, GEN l, long siz)
586 : {
587 575926 : pari_sp av = avma;
588 : GEN p1;
589 575926 : (void)new_chunk(siz); /* HACK */
590 575926 : p1 = mulii(x,y);
591 575926 : set_avma(av); return modii(p1,l);
592 : }
593 :
594 : static GEN
595 69710 : polsubcyclo_roots(long n, GEN zl)
596 : {
597 69710 : GEN le = gel(zl,1), z = gel(zl,2);
598 69710 : long i, lle = lg(le)*3; /*Assume dvmdii use lx+ly space*/
599 69710 : long m = (long)(1+sqrt((double) n));
600 69710 : GEN bab, gig, powz = cgetg(3,t_VEC);
601 : pari_timer ti;
602 69710 : if (DEBUGLEVEL >= 6) timer_start(&ti);
603 69710 : bab = cgetg(m+1,t_VEC);
604 69710 : gel(bab,1) = gen_1;
605 69710 : gel(bab,2) = icopy(z);
606 322818 : for (i=3; i<=m; i++) gel(bab,i) = muliimod_sz(z,gel(bab,i-1),le,lle);
607 69710 : gig = cgetg(m+1,t_VEC);
608 69710 : gel(gig,1) = gen_1;
609 69710 : gel(gig,2) = muliimod_sz(z,gel(bab,m),le,lle);;
610 322818 : for (i=3; i<=m; i++) gel(gig,i) = muliimod_sz(gel(gig,2),gel(gig,i-1),le,lle);
611 69710 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "polsubcyclo_roots");
612 69710 : gel(powz,1) = bab;
613 69710 : gel(powz,2) = gig; return powz;
614 : }
615 :
616 : GEN
617 384 : galoiscyclo(long n, long v)
618 : {
619 384 : ulong av = avma;
620 : GEN grp, G, z, le, L, elts;
621 : long val, l, i, j, k;
622 384 : GEN zn = znstar(stoi(n));
623 384 : long card = itos(gel(zn,1));
624 384 : GEN gen = vec_to_vecsmall(lift_shallow(gel(zn,3)));
625 384 : GEN ord = vec_to_vecsmall(gel(zn,2));
626 384 : GEN T = polcyclo(n,v);
627 384 : long d = degpol(T);
628 384 : GEN borneabs = powuu(2,d);
629 384 : z = polsubcyclo_start(n,card/2,2,2*usqrt(d),borneabs,&val,&l);
630 384 : le = gel(z,1); z = gel(z,2);
631 384 : L = cgetg(1+card,t_VEC);
632 384 : gel(L,1) = z;
633 774 : for (j = 1, i = 1; j < lg(gen); j++)
634 : {
635 390 : long c = i * (ord[j]-1);
636 1860 : for (k = 1; k <= c; k++) gel(L,++i) = Fp_powu(gel(L,k), gen[j], le);
637 : }
638 384 : G = abelian_group(ord);
639 384 : elts = group_elts(G, card); /*not stack clean*/
640 384 : grp = cgetg(9, t_VEC);
641 384 : gel(grp,1) = T;
642 384 : gel(grp,2) = mkvec3(stoi(l), stoi(val), icopy(le));
643 384 : gel(grp,3) = L;
644 384 : gel(grp,4) = FpV_invVandermonde(L, NULL, le);
645 384 : gel(grp,5) = gen_1;
646 384 : gel(grp,6) = elts;
647 384 : gel(grp,7) = gel(G,1);
648 384 : gel(grp,8) = gel(G,2);
649 384 : return gerepilecopy(av, grp);
650 : }
651 :
652 : /* Convert a bnrinit(Q,n) to an abelian group similar to znstar(n), with
653 : * t_INTMOD generators; set cx = 0 if the class field is real and to 1
654 : * otherwise */
655 : static GEN
656 48 : bnr_to_abgrp(GEN bnr, long *cx)
657 : {
658 48 : GEN gen, F, v, bid, G, Ui = NULL;
659 : long l, i;
660 48 : checkbnr(bnr);
661 48 : bid = bnr_get_bid(bnr);
662 48 : G = bnr_get_clgp(bnr);
663 48 : if (lg(G) == 4)
664 24 : gen = abgrp_get_gen(G);
665 : else
666 : {
667 24 : Ui = gmael(bnr,4,3);
668 24 : if (ZM_isidentity(Ui)) Ui = NULL;
669 24 : gen = bid_get_gen(bid);
670 : }
671 48 : F = bid_get_ideal(bid);
672 48 : if (lg(F) != 2)
673 6 : pari_err_DOMAIN("bnr_to_abgrp", "bnr", "!=", strtoGENstr("Q"), bnr);
674 : /* F is the finite part of the conductor, cx is the infinite part*/
675 42 : F = gcoeff(F, 1, 1);
676 42 : *cx = signe(gel(bid_get_arch(bid), 1));
677 42 : l = lg(gen); v = cgetg(l, t_VEC);
678 174 : for (i = 1; i < l; ++i)
679 : {
680 132 : GEN x = gel(gen,i);
681 132 : if (typ(x) == t_COL) x = gel(x,1);
682 132 : gel(v,i) = gmodulo(absi_shallow(x), F);
683 : }
684 42 : if (Ui)
685 : { /* from bid.gen to bnr.gen (maybe one less) */
686 18 : GEN w = v;
687 18 : l = lg(Ui); v = cgetg(l, t_VEC);
688 72 : for (i = 1; i < l; i++) gel(v,i) = factorback2(w, gel(Ui, i));
689 : }
690 42 : return mkvec3(bnr_get_no(bnr), bnr_get_cyc(bnr), v);
691 : }
692 :
693 : static long
694 43146 : _itos(const char *fun, GEN n)
695 : {
696 43146 : if (is_bigint(n))
697 6 : pari_err_IMPL(stack_sprintf("conductor f > %ld in %s", LONG_MAX, fun));
698 43140 : return itos(n);
699 : }
700 : long
701 43158 : subcyclo_nH(const char *fun, GEN N, GEN *psg)
702 : {
703 43158 : GEN V, Z = NULL, H = *psg;
704 43158 : long i, l, n = 0, complex = 1;
705 43158 : switch(typ(N))
706 : {
707 42780 : case t_INT:
708 42780 : n = _itos(fun, N);
709 42774 : if (n < 1) pari_err_DOMAIN(fun, "degree", "<=", gen_0, N);
710 42768 : break;
711 378 : case t_VEC:
712 378 : if (lg(N)==7)
713 48 : N = bnr_to_abgrp(N,&complex);
714 330 : else if (checkznstar_i(N))
715 12 : N = mkvec3(znstar_get_no(N), znstar_get_cyc(N),
716 : gmodulo(znstar_get_gen(N), znstar_get_N(N)));
717 372 : if (lg(N)==4)
718 : { /* abgrp */
719 372 : GEN gen = abgrp_get_gen(N), z;
720 372 : if (typ(gen)!=t_VEC) pari_err_TYPE(fun,gen);
721 372 : Z = N;
722 372 : if (lg(gen) == 1) { n = 1; break; }
723 372 : z = gel(gen,1);
724 372 : if (typ(z) == t_INTMOD) { n = _itos(fun, gel(z,1)); break; }
725 : }
726 : default: /*fall through*/
727 6 : pari_err_TYPE(fun,N);
728 : return 0;/*LCOV_EXCL_LINE*/
729 : }
730 43134 : if (!H) H = gen_1;
731 43134 : switch(typ(H))
732 : {
733 42612 : case t_INTMOD: case t_INT:
734 42612 : V = mkvecsmall( lift_check_modulus(H,n) );
735 42606 : break;
736 12 : case t_VECSMALL:
737 12 : l = lg(H); V = leafcopy(H);
738 36 : for (i = 1; i < l; i++) { V[i] %= n; if (V[i] < 0) V[i] += n; }
739 12 : break;
740 132 : case t_VEC: case t_COL:
741 132 : l = lg(H); V = cgetg(l,t_VECSMALL);
742 1104318 : for(i=1; i < l; i++) V[i] = lift_check_modulus(gel(H,i),n);
743 126 : break;
744 378 : case t_MAT:
745 378 : l = lg(H);
746 378 : if (l == 1 || l != lgcols(H))
747 6 : pari_err_TYPE(stack_strcat(fun," [H not in HNF]"),H);
748 372 : if (!Z) pari_err_TYPE(stack_strcat(fun," [N not a bnrinit or znstar]"),H);
749 366 : if (lg(gel(Z,2)) != l) pari_err_DIM(fun);
750 360 : V = znstar_hnf_generators(znstar_small(Z),H);
751 360 : break;
752 0 : default:
753 0 : pari_err_TYPE(fun,H);
754 : return 0;/*LCOV_EXCL_LINE*/
755 : }
756 43104 : if (!complex) V = vecsmall_append(V, n-1); /*add complex conjugation*/
757 43104 : *psg = V; return n;
758 : }
759 :
760 : GEN
761 42828 : galoissubcyclo(GEN N, GEN sg, long flag, long v)
762 : {
763 42828 : pari_sp ltop = avma, av;
764 : GEN H, B, zl, L, T, le, powz, O;
765 : long i, card, phi_n, val,l, n, cnd, complex;
766 : pari_timer ti;
767 :
768 42828 : if (flag<0 || flag>3) pari_err_FLAG("galoissubcyclo");
769 42828 : if (v < 0) v = 0;
770 42828 : n = subcyclo_nH("galoissubcyclo", N, &sg);
771 42786 : if (n==1)
772 : {
773 18 : set_avma(ltop); if (flag == 1) return gen_1;
774 12 : return gscycloconductor(deg1pol_shallow(gen_1, gen_m1, v), 1, flag);
775 : }
776 42768 : H = znstar_generate(n, sg);
777 42768 : if (DEBUGLEVEL >= 6)
778 : {
779 0 : err_printf("Subcyclo: elements:");
780 0 : for (i=1;i<n;i++)
781 0 : if (F2v_coeff(gel(H,3),i)) err_printf(" %ld",i);
782 0 : err_printf("\n");
783 : }
784 : /* field is real iff z -> conj(z) = z^-1 = z^(n-1) is in H */
785 42768 : complex = !F2v_coeff(gel(H,3),n-1);
786 42768 : if (DEBUGLEVEL >= 6) err_printf("Subcyclo: complex=%ld\n",complex);
787 42768 : if (DEBUGLEVEL >= 1) timer_start(&ti);
788 42768 : cnd = znstar_conductor(H);
789 42768 : if (DEBUGLEVEL >= 1) timer_printf(&ti, "znstar_conductor");
790 42768 : if (flag == 1) return gc_stoi(ltop, cnd);
791 42768 : if (cnd == 1)
792 : {
793 60 : set_avma(ltop); if (flag == 1) return gen_1;
794 60 : return gscycloconductor(deg1pol_shallow(gen_1,gen_m1,v),1,flag);
795 : }
796 42708 : if (n != cnd)
797 : {
798 150 : H = znstar_reduce_modulus(H, cnd);
799 150 : n = cnd;
800 : }
801 42708 : card = znstar_order(H);
802 42708 : phi_n = eulerphiu(n);
803 42708 : if (card == phi_n)
804 : {
805 0 : set_avma(ltop);
806 0 : return gscycloconductor(polcyclo(n,v),n,flag);
807 : }
808 42708 : O = znstar_cosets(n, phi_n, H);
809 42708 : if (DEBUGLEVEL >= 1) timer_printf(&ti, "znstar_cosets");
810 42708 : if (DEBUGLEVEL >= 6) err_printf("Subcyclo: orbits=%Ps\n",O);
811 42708 : if (DEBUGLEVEL >= 4)
812 0 : err_printf("Subcyclo: %ld orbits with %ld elements each\n",phi_n/card,card);
813 42708 : av = avma;
814 42708 : powz = polsubcyclo_complex_roots(n,!complex,LOWDEFAULTPREC);
815 42708 : L = polsubcyclo_orbits(n,H,O,powz,NULL);
816 42708 : B = polsubcyclo_complex_bound(av,L,LOWDEFAULTPREC);
817 42708 : zl = polsubcyclo_start(n,phi_n/card,card,1,B,&val,&l);
818 42708 : powz = polsubcyclo_roots(n,zl);
819 42708 : le = gel(zl,1);
820 42708 : L = polsubcyclo_orbits(n,H,O,powz,le);
821 42708 : if (DEBUGLEVEL >= 6) timer_start(&ti);
822 42708 : T = FpV_roots_to_pol(L,le,v);
823 42708 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "roots_to_pol");
824 42708 : T = FpX_center(T,le,shifti(le,-1));
825 42708 : if (flag==3)
826 : {
827 : GEN L2, aut;
828 6 : if (Flx_is_squarefree(ZX_to_Flx(T, l),l))
829 6 : L2 = ZV_to_Flv(L,l);
830 : else
831 : {
832 : ulong z;
833 0 : do l+=n; while (!uisprime(l) || !Flx_is_squarefree(ZX_to_Flx(T, l), l));
834 0 : z = rootsof1_Fl(n,l);
835 0 : L2 = Fl_polsubcyclo_orbits(n,H,O,z,l);
836 0 : if (DEBUGLEVEL >= 4)
837 0 : err_printf("galoissubcyclo: switching to unramified prime %lu\n",l);
838 : }
839 6 : aut = znstar_quotient(n, phi_n, H, L2, l);
840 6 : return gerepileupto(ltop, galoisinitfromaut(T, aut, l));
841 : }
842 42702 : return gerepileupto(ltop, gscycloconductor(T,n,flag));
843 : }
844 :
845 : /* Z = znstar(n) cyclic. n = 1,2,4,p^a or 2p^a,
846 : * and d | phi(n) = 1,1,2,(p-1)p^(a-1) */
847 : static GEN
848 29000 : polsubcyclo_g(long n, long d, GEN Z, long v)
849 : {
850 29000 : pari_sp ltop = avma;
851 : long o, p, r, g, gd, l , val;
852 : GEN zl, L, T, le, B, powz;
853 : pari_timer ti;
854 29000 : if (d==1) return deg1pol_shallow(gen_1,gen_m1,v); /* get rid of n=1,2 */
855 29000 : if ((n & 3) == 2) n >>= 1;
856 : /* n = 4 or p^a, p odd */
857 29000 : o = itos(gel(Z,1));
858 29000 : g = itos(gmael3(Z,3,1,2));
859 29000 : p = n / ugcd(n,o); /* p^a / gcd(p^a,phi(p^a)) = p*/
860 29000 : r = ugcd(d,n); /* = p^(v_p(d)) < n */
861 29000 : n = r*p; /* n is now the conductor */
862 29000 : o = n-r; /* = phi(n) */
863 29000 : if (o == d) return polcyclo(n,v);
864 27002 : o /= d;
865 27002 : gd = Fl_powu(g%n, d, n);
866 : /*FIXME: If degree is small, the computation of B is a waste of time*/
867 27002 : powz = polsubcyclo_complex_roots(n,(o&1)==0,LOWDEFAULTPREC);
868 27002 : L = polsubcyclo_cyclic(n,d,o,g,gd,powz,NULL);
869 27002 : B = polsubcyclo_complex_bound(ltop,L,LOWDEFAULTPREC);
870 27002 : zl = polsubcyclo_start(n,d,o,1,B,&val,&l);
871 27002 : le = gel(zl,1);
872 27002 : powz = polsubcyclo_roots(n,zl);
873 27002 : L = polsubcyclo_cyclic(n,d,o,g,gd,powz,le);
874 27002 : if (DEBUGLEVEL >= 6) timer_start(&ti);
875 27002 : T = FpV_roots_to_pol(L,le,v);
876 27002 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "roots_to_pol");
877 27002 : return gerepileupto(ltop, FpX_center(T,le,shifti(le,-1)));
878 : }
879 :
880 : GEN
881 29048 : polsubcyclo(long n, long d, long v)
882 : {
883 29048 : pari_sp ltop = avma;
884 : GEN L, Z;
885 29048 : if (v<0) v = 0;
886 29048 : if (d<=0) pari_err_DOMAIN("polsubcyclo","d","<=",gen_0,stoi(d));
887 29042 : if (n<=0) pari_err_DOMAIN("polsubcyclo","n","<=",gen_0,stoi(n));
888 29036 : Z = znstar(stoi(n));
889 29036 : if (!dvdis(gel(Z,1), d)) { set_avma(ltop); return cgetg(1, t_VEC); }
890 29030 : if (lg(gel(Z,2)) == 2)
891 : { /* faster but Z must be cyclic */
892 29000 : set_avma(ltop);
893 29000 : return polsubcyclo_g(n, d, Z, v);
894 : }
895 30 : L = subgrouplist(gel(Z,2), mkvec(stoi(d)));
896 30 : if (lg(L) == 2)
897 6 : return gerepileupto(ltop, galoissubcyclo(Z, gel(L,1), 0, v));
898 : else
899 : {
900 24 : GEN V = cgetg(lg(L),t_VEC);
901 : long i;
902 312 : for (i=1; i< lg(V); i++) gel(V,i) = galoissubcyclo(Z, gel(L,i), 0, v);
903 24 : return gerepileupto(ltop, V);
904 : }
905 : }
906 :
907 : struct aurifeuille_t {
908 : GEN z, le;
909 : ulong l;
910 : long e;
911 : };
912 :
913 : /* Let z a primitive n-th root of 1, n > 1, A an integer such that
914 : * Aurifeuillian factorization of Phi_n(A) exists ( z.A is a square in Q(z) ).
915 : * Let G(p) the Gauss sum mod p prime:
916 : * sum_x (x|p) z^(xn/p) for p odd, i - 1 for p = 2 [ i := z^(n/4) ]
917 : * We have N(-1) = Nz = 1 (n != 1,2), and
918 : * G^2 = (-1|p) p for p odd, G^2 = -2i for p = 2
919 : * In particular, for odd A, (-1|A) A = g^2 is a square. If A = prod p^{e_p},
920 : * sigma_j(g) = \prod_p (sigma_j G(p)))^e_p = \prod_p (j|p)^e_p g = (j|A) g
921 : * n odd : z^2 is a primitive root, A = g^2
922 : * Phi_n(A) = N(A - z^2) = N(g - z) N(g + z)
923 : *
924 : * n = 2 (4) : -z^2 is a primitive root, -A = g^2
925 : * Phi_n(A) = N(A - (-z^2)) = N(g^2 - z^2) [ N(-1) = 1 ]
926 : * = N(g - z) N(g + z)
927 : *
928 : * n = 4 (8) : i z^2 primitive root, -Ai = g^2
929 : * Phi_n(A) = N(A - i z^2) = N(-Ai - z^2) = N(g - z) N(g + z)
930 : * sigma_j(g) / g = (j|A) if j = 1 (4)
931 : * (-j|A)i if j = 3 (4)
932 : * */
933 : /* factor Phi_n(A), Astar: A* = squarefree kernel of A, P = odd prime divisors
934 : * of n */
935 : static GEN
936 366 : factor_Aurifeuille_aux(GEN A, long Astar, long n, GEN P,
937 : struct aurifeuille_t *S)
938 : {
939 : pari_sp av;
940 366 : GEN f, a, b, s, powers, z = S->z, le = S->le;
941 366 : long j, k, maxjump, lastj, e = S->e;
942 366 : ulong l = S->l;
943 : char *invertible;
944 :
945 366 : if ((n & 7) == 4)
946 : { /* A^* even */
947 300 : GEN i = Fp_powu(z, n>>2, le), z2 = Fp_sqr(z, le);
948 :
949 300 : invertible = stack_malloc(n); /* even indices unused */
950 1380 : for (j = 1; j < n; j+=2) invertible[j] = 1;
951 348 : for (k = 1; k < lg(P); k++)
952 : {
953 48 : long p = P[k];
954 144 : for (j = p; j < n; j += 2*p) invertible[j] = 0;
955 : }
956 300 : lastj = 1; maxjump = 2;
957 1080 : for (j= 3; j < n; j+=2)
958 780 : if (invertible[j]) {
959 684 : long jump = j - lastj;
960 684 : if (jump > maxjump) maxjump = jump;
961 684 : lastj = j;
962 : }
963 300 : powers = cgetg(maxjump+1, t_VEC); /* powers[k] = z^k, odd indices unused */
964 300 : gel(powers,2) = z2;
965 348 : for (k = 4; k <= maxjump; k+=2)
966 96 : gel(powers,k) = odd(k>>1)? Fp_mul(gel(powers, k-2), z2, le)
967 48 : : Fp_sqr(gel(powers, k>>1), le);
968 :
969 300 : if (Astar == 2)
970 : { /* important special case (includes A=2), split for efficiency */
971 282 : if (!equalis(A, 2))
972 : {
973 12 : GEN f = sqrti(shifti(A,-1)), mf = Fp_neg(f,le), fi = Fp_mul(f,i,le);
974 12 : a = Fp_add(mf, fi, le);
975 12 : b = Fp_sub(mf, fi, le);
976 : }
977 : else
978 : {
979 270 : a = subiu(i,1);
980 270 : b = subsi(-1,i);
981 : }
982 282 : av = avma;
983 282 : s = z; f = subii(a, s); lastj = 1;
984 780 : for (j = 3, k = 0; j < n; j+=2)
985 498 : if (invertible[j])
986 : {
987 438 : s = Fp_mul(gel(powers, j-lastj), s, le); /* z^j */
988 438 : lastj = j;
989 438 : f = Fp_mul(f, subii((j & 3) == 1? a: b, s), le);
990 438 : if (++k == 0x1ff) { gerepileall(av, 2, &s, &f); k = 0; }
991 : }
992 : }
993 : else
994 : {
995 18 : GEN ma, mb, B = Fp_mul(A, i, le), gl = utoipos(l);
996 : long t;
997 18 : Astar >>= 1;
998 18 : t = Astar & 3; if (Astar < 0) t = 4-t; /* t = 1 or 3 */
999 18 : if (t == 1) B = Fp_neg(B, le);
1000 18 : a = Zp_sqrtlift(B, Fp_sqrt(B, gl), gl, e);
1001 18 : b = Fp_mul(a, i, le);
1002 18 : ma = Fp_neg(a, le);
1003 18 : mb = Fp_neg(b, le);
1004 18 : av = avma;
1005 18 : s = z; f = subii(a, s); lastj = 1;
1006 300 : for (j = 3, k = 0; j<n; j+=2)
1007 282 : if (invertible[j])
1008 : {
1009 : GEN t;
1010 246 : if ((j & 3) == 1) t = (kross(j, Astar) < 0)? ma: a;
1011 132 : else t = (kross(j, Astar) < 0)? mb: b;
1012 246 : s = Fp_mul(gel(powers, j-lastj), s, le); /* z^j */
1013 246 : lastj = j;
1014 246 : f = Fp_mul(f, subii(t, s), le);
1015 246 : if (++k == 0x1ff) { gerepileall(av, 2, &s, &f); k = 0; }
1016 : }
1017 : }
1018 : }
1019 : else /* A^* odd */
1020 : {
1021 : ulong g;
1022 66 : if ((n & 3) == 2)
1023 : { /* A^* = 3 (mod 4) */
1024 0 : A = negi(A); Astar = -Astar;
1025 0 : z = Fp_neg(z, le);
1026 0 : n >>= 1;
1027 : }
1028 : /* A^* = 1 (mod 4) */
1029 66 : g = Fl_sqrt(umodiu(A,l), l);
1030 66 : a = Zp_sqrtlift(A, utoipos(g), utoipos(l), e);
1031 66 : b = negi(a);
1032 :
1033 66 : invertible = stack_malloc(n);
1034 1146 : for (j = 1; j < n; j++) invertible[j] = 1;
1035 168 : for (k = 1; k < lg(P); k++)
1036 : {
1037 102 : long p = P[k];
1038 414 : for (j = p; j < n; j += p) invertible[j] = 0;
1039 : }
1040 66 : lastj = 2; maxjump = 1;
1041 1014 : for (j= 3; j < n; j++)
1042 948 : if (invertible[j]) {
1043 636 : long jump = j - lastj;
1044 636 : if (jump > maxjump) maxjump = jump;
1045 636 : lastj = j;
1046 : }
1047 66 : powers = cgetg(maxjump+1, t_VEC); /* powers[k] = z^k */
1048 66 : gel(powers,1) = z;
1049 138 : for (k = 2; k <= maxjump; k++)
1050 144 : gel(powers,k) = odd(k)? Fp_mul(gel(powers, k-1), z, le)
1051 72 : : Fp_sqr(gel(powers, k>>1), le);
1052 66 : av = avma;
1053 66 : s = z; f = subii(a, s); lastj = 1;
1054 1080 : for(j = 2, k = 0; j < n; j++)
1055 1014 : if (invertible[j])
1056 : {
1057 702 : s = Fp_mul(gel(powers, j-lastj), s, le);
1058 702 : lastj = j;
1059 702 : f = Fp_mul(f, subii(kross(j,Astar)==1? a: b, s), le);
1060 702 : if (++k == 0x1ff) { gerepileall(av, 2, &s, &f); k = 0; }
1061 : }
1062 : }
1063 366 : return f;
1064 : }
1065 :
1066 : /* d != 2 mod 4; fd = factoru(odd(d)? d: d / 4) */
1067 : static void
1068 366 : Aurifeuille_init(GEN a, long d, GEN fd, struct aurifeuille_t *S)
1069 : {
1070 366 : GEN bound, zl, sqrta = sqrtr_abs(itor(a, LOWDEFAULTPREC));
1071 366 : ulong phi = eulerphiu_fact(fd);
1072 366 : if (!odd(d)) phi <<= 1; /* eulerphi(d) */
1073 366 : bound = ceil_safe(powru(addrs(sqrta,1), phi));
1074 366 : zl = polsubcyclo_start(d, 0, 0, 1, bound, &(S->e), (long*)&(S->l));
1075 366 : S->le = gel(zl,1);
1076 366 : S->z = gel(zl,2);
1077 366 : }
1078 :
1079 : GEN
1080 312 : factor_Aurifeuille_prime(GEN p, long d)
1081 : {
1082 312 : pari_sp av = avma;
1083 : struct aurifeuille_t S;
1084 : GEN fd;
1085 : long pp;
1086 312 : if ((d & 3) == 2) { d >>= 1; p = negi(p); }
1087 312 : fd = factoru(odd(d)? d: d>>2);
1088 312 : pp = itos(p);
1089 312 : Aurifeuille_init(p, d, fd, &S);
1090 312 : return gerepileuptoint(av, factor_Aurifeuille_aux(p, pp, d, gel(fd,1), &S));
1091 : }
1092 :
1093 : /* an algebraic factor of Phi_d(a), a != 0 */
1094 : GEN
1095 54 : factor_Aurifeuille(GEN a, long d)
1096 : {
1097 54 : pari_sp av = avma;
1098 : GEN fd, P, A;
1099 54 : long i, lP, va = vali(a), sa, astar, D;
1100 : struct aurifeuille_t S;
1101 :
1102 54 : if (d <= 0)
1103 0 : pari_err_DOMAIN("factor_Aurifeuille", "degre", "<=",gen_0,stoi(d));
1104 54 : if ((d & 3) == 2) { d >>= 1; a = negi(a); }
1105 54 : if ((va & 1) == (d & 1)) { set_avma(av); return gen_1; }
1106 54 : sa = signe(a);
1107 54 : if (odd(d))
1108 : {
1109 : long a4;
1110 24 : if (d == 1)
1111 : {
1112 0 : if (!Z_issquareall(a, &A)) return gen_1;
1113 0 : return gerepileuptoint(av, addiu(A,1));
1114 : }
1115 24 : A = va? shifti(a, -va): a;
1116 24 : a4 = mod4(A); if (sa < 0) a4 = 4 - a4;
1117 24 : if (a4 != 1) { set_avma(av); return gen_1; }
1118 : }
1119 30 : else if ((d & 7) == 4)
1120 30 : A = shifti(a, -va);
1121 : else
1122 : {
1123 0 : set_avma(av); return gen_1;
1124 : }
1125 : /* v_2(d) = 0 or 2. Kill 2 from factorization (minor efficiency gain) */
1126 54 : fd = factoru(odd(d)? d: d>>2); P = gel(fd,1); lP = lg(P);
1127 54 : astar = sa;
1128 54 : if (odd(va)) astar <<= 1;
1129 126 : for (i = 1; i < lP; i++)
1130 72 : if (odd( (Z_lvalrem(A, P[i], &A)) ) ) astar *= P[i];
1131 54 : if (sa < 0)
1132 : { /* negate in place if possible */
1133 12 : if (A == a) A = icopy(A);
1134 12 : setabssign(A);
1135 : }
1136 54 : if (!Z_issquare(A)) { set_avma(av); return gen_1; }
1137 :
1138 54 : D = odd(d)? 1: 4;
1139 126 : for (i = 1; i < lP; i++) D *= P[i];
1140 54 : if (D != d) { a = powiu(a, d/D); d = D; }
1141 :
1142 54 : Aurifeuille_init(a, d, fd, &S);
1143 54 : return gerepileuptoint(av, factor_Aurifeuille_aux(a, astar, d, P, &S));
1144 : }
|
