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 1349303 : 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 1349303 : if (!d) { (*func)(data,c); return; }
40 :
41 1093889 : cache = const_vecsmall(d,c);
42 1093889 : (*func)(data,c); /* AFTER cache: may contain gerepileupto statement */
43 1093889 : gen = gel(H,1);
44 1093889 : ord = gel(H,2);
45 1154740 : card = ord[1]; for (i = 2; i <= d; i++) card *= ord[i];
46 82530811 : for(i=1; i<card; i++)
47 : {
48 81436922 : long k, m = i;
49 82010250 : for(j=1; j<d && m%ord[j]==0 ;j++) m /= ord[j];
50 81436922 : cache[j] = Fl_mul(cache[j],gen[j],n);
51 82010250 : for (k=1; k<j; k++) cache[k] = cache[j];
52 81436922 : (*func)(data, cache[j]);
53 : }
54 : }
55 :
56 : static void
57 275778 : znstar_coset_func(long n, GEN H, void (*func)(void *data,long c)
58 : , void *data, long c)
59 275778 : { 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 1073525 : znstar_partial_coset_bits_inplace(long n, GEN H, GEN bits, long d, long c)
66 : {
67 1073525 : pari_sp av = avma;
68 1073525 : znstar_partial_coset_func(n,H, (void (*)(void *,long)) &F2v_set,
69 : (void *) bits, d, c);
70 1073525 : set_avma(av);
71 1073525 : }
72 :
73 : static void
74 531855 : znstar_coset_bits_inplace(long n, GEN H, GEN bits, long c)
75 531855 : { znstar_partial_coset_bits_inplace(n, H, bits, lg(gel(H,1))-1, c); }
76 :
77 : static GEN
78 541670 : znstar_partial_coset_bits(long n, GEN H, long d, long c)
79 : {
80 541670 : GEN bits = zero_F2v(n);
81 541670 : znstar_partial_coset_bits_inplace(n,H,bits,d,c);
82 541670 : 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 541670 : znstar_partial_bits(long n, GEN H, long d)
89 541670 : { 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 255036 : znstar_generate(long n, GEN V)
100 : {
101 255036 : pari_sp av = avma;
102 255036 : GEN gen = cgetg(lg(V),t_VECSMALL);
103 255036 : GEN ord = cgetg(lg(V),t_VECSMALL), res = mkvec2(gen,ord);
104 255036 : GEN bits = znstar_partial_bits(n,NULL,0);
105 255036 : long i, r = 0;
106 1850922 : for(i=1; i<lg(V); i++)
107 : {
108 1595886 : ulong v = uel(V,i), g = v;
109 1595886 : long o = 0;
110 38258254 : while (!F2v_coeff(bits, (long)g)) { g = Fl_mul(g, v, (ulong)n); o++; }
111 1595886 : if (!o) continue;
112 286634 : r++;
113 286634 : gen[r] = v;
114 286634 : ord[r] = o+1;
115 286634 : cgiv(bits); bits = znstar_partial_bits(n,res,r);
116 : }
117 255036 : setlg(gen,r+1);
118 255036 : setlg(ord,r+1); return gerepilecopy(av, mkvec3(gen,ord,bits));
119 : }
120 :
121 : static ulong
122 164044 : 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 3780 : znstar_elts(long n, GEN H)
128 : {
129 3780 : long card = znstar_order(H);
130 3780 : GEN gen = gel(H,1), ord = gel(H,2);
131 3780 : GEN sg = cgetg(1 + card, t_VECSMALL);
132 : long k, j, l;
133 3780 : sg[1] = 1;
134 5796 : for (j = 1, l = 1; j < lg(gen); j++)
135 : {
136 2016 : long c = l * (ord[j]-1);
137 4578 : for (k = 1; k <= c; k++) sg[++l] = Fl_mul(sg[k], gen[j], n);
138 : }
139 3780 : 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 182 : znstar_reduce_modulus(GEN H, long n)
146 : {
147 182 : pari_sp ltop=avma;
148 182 : GEN gen=cgetg(lgcols(H),t_VECSMALL);
149 : long i;
150 637 : for(i=1; i < lg(gen); i++)
151 455 : gen[i] = mael(H,1,i)%n;
152 182 : return gerepileupto(ltop, znstar_generate(n,gen));
153 : }
154 :
155 : /* Compute conductor of H, bits = H[3] */
156 : long
157 172145 : znstar_conductor_bits(GEN bits)
158 : {
159 172145 : pari_sp av = avma;
160 172145 : long i, f = 1, cnd0 = bits[1];
161 172145 : GEN F = factoru(cnd0), P = gel(F,1), E = gel(F,2);
162 462988 : for (i = lg(P)-1; i > 0; i--)
163 : {
164 290843 : long p = P[i], e = E[i], cnd = cnd0;
165 328293 : for ( ; e >= 2; e--)
166 : {
167 83879 : long q = cnd / p;
168 83879 : if (!F2v_coeff(bits, 1 + q)) break;
169 37450 : cnd = q;
170 : }
171 290843 : if (e == 1)
172 : {
173 244414 : if (p == 2) e = 0;
174 : else
175 : {
176 215693 : long h, g = pgener_Fl(p), q = cnd / p;
177 215693 : h = Fl_mul(g-1, Fl_inv(q % p, p), p); /* 1+h*q = g (mod p) */
178 215693 : if (F2v_coeff(bits, 1 + h*q)) e = 0;
179 : }
180 : }
181 290843 : if (e) f *= upowuu(p, e);
182 : }
183 172145 : return gc_long(av,f);
184 : }
185 : long
186 171487 : 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 110514 : znstar_cosets(long n, long phi_n, GEN H)
193 : {
194 110514 : long k, c = 0, card = znstar_order(H), index = phi_n/card;
195 110514 : GEN cosets = cgetg(index+1,t_VECSMALL);
196 110514 : pari_sp ltop = avma;
197 110514 : GEN bits = zero_F2v(n);
198 642327 : for (k = 1; k <= index; k++)
199 : {
200 1371942 : for (c++ ; F2v_coeff(bits,c) || ugcd(c,n)!=1; c++);
201 531813 : cosets[k]=c;
202 531813 : znstar_coset_bits_inplace(n, H, bits, c);
203 : }
204 110514 : set_avma(ltop); return cosets;
205 : }
206 :
207 : static GEN
208 7 : znstar_quotient(long n, long phi_n, GEN H, GEN R, ulong l)
209 : {
210 7 : long i, j, k, c = 0, card = znstar_order(H), index = phi_n/card;
211 7 : GEN cosets = cgetg(index+1,t_VECSMALL), mult = cgetg(index+1, t_VEC);
212 7 : GEN bits = zero_F2v(n), vbits= zero_F2m_copy(n,index);
213 49 : for (k = 1; k <= index; k++)
214 : {
215 119 : for (c++ ; F2v_coeff(bits,c) || ugcd(c,n)!=1; c++);
216 42 : cosets[k]=c;
217 42 : znstar_coset_bits_inplace(n, H, gel(vbits,k), c);
218 42 : F2v_add_inplace(bits, gel(vbits,k));
219 : }
220 :
221 49 : for (k = 1; k <= index; k++)
222 : {
223 42 : GEN v = cgetg(index+1, t_VECSMALL);
224 294 : for (j = 1; j <= index; j++)
225 : {
226 252 : long s = Fl_mul(cosets[k],cosets[j],n);
227 882 : for (i = 1; i <= index; i++)
228 882 : if (F2v_coeff(gel(vbits,i),s)) break;
229 252 : v[j] = R[i];
230 : }
231 42 : gel(mult,k) = v;
232 : }
233 7 : return Flv_Flm_polint(R, mult, l, 0);
234 : }
235 :
236 : /* return n s.t. x^n in H */
237 : static ulong
238 1358266 : order_H_x(GEN H, ulong x, GEN D, ulong *xn)
239 : {
240 1358266 : ulong i, l = lg(D), f = H[1], y = x; /* = x^D[1] */
241 1358266 : *xn = x; if (F2v_coeff(H, y)) return D[1];
242 26096126 : for (i = 2; i < l; i++)
243 : { /* TODO: could cache the x^delta[i] incrementally */
244 26096126 : y = Fl_mul(y, Fl_powu(x, D[i]-D[i-1], f), f);
245 26096126 : if (F2v_coeff(H, y))
246 : {
247 1358266 : if (xn) *xn = y;
248 1358266 : 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 70 : max_order_ele(GEN H1, GEN H2, GEN D, ulong *po, long flag)
258 : {
259 70 : ulong x, h, n = F2v_hamming(H2), f = H1[1], O = 0, X = 0, XO = 0;
260 70 : if (!n) return 0;
261 40680584 : for (x = 2; x < f; x++) if (F2v_coeff(H2, x))
262 : {
263 1358266 : ulong xo, o = order_H_x(H1, x, D, &xo);
264 1358266 : if (o > O) { O = o; X = x; XO = xo; if (o == n) break; }
265 : }
266 : /* X of maximal order O in H2/H1 */
267 56 : *po = O; if (!flag || XO == 1) return X;
268 : /* find h in H1 s.t. (Xh)^o=1 */
269 21 : x = Fl_inv(XO, f);
270 1857611 : for (h = 1;; h++) /* stops for h < f */
271 1857611 : if (F2v_coeff(H1, h) && x == Fl_powu(h, O, f)) break;
272 21 : 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 56 : enlarge_H(GEN H, ulong x, ulong o)
278 : {
279 56 : pari_sp av = avma;
280 56 : ulong i, j, l = lg(H), f = H[1];
281 56 : GEN H1 = vecsmall_copy(H), y = Fl_powers(x, o-1, f);
282 40680640 : for (i = 1; i < f; i++)
283 40680584 : if (F2v_coeff(H, i))
284 758520 : for (j = 2; j <= o; j++) F2v_set(H1, Fl_mul(i, y[j], f));
285 726520 : for (i = 2; i < l; i++) H[i] = H1[i];
286 56 : set_avma(av);
287 56 : }
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 21 : znsubgroupgenerators(GEN H, long flag)
292 : {
293 21 : pari_sp av = avma;
294 : ulong f, g, o;
295 21 : GEN H1, H2, z = const_vecsmall(0,0), D;
296 21 : switch(typ(H))
297 : {
298 14 : case t_VECSMALL: H = Flv_to_F2v(H); break;
299 7 : case t_VEC: if (RgV_is_ZV(H)) { H = ZV_to_F2v(H); break; }
300 7 : default: pari_err_TYPE("znsubgroupgenerators", H);
301 : }
302 14 : f = H[1]; z = cgetg(1, t_VECSMALL);
303 14 : D = divisorsu(F2v_hamming(H));
304 14 : H1 = zero_F2v(f); F2v_set(H1, 1);
305 14 : H2 = H;
306 70 : while ((g = max_order_ele(H1, H2, D, &o, flag)))
307 : {
308 56 : z = vecsmall_append(z, g); enlarge_H(H1, g, o);
309 56 : F2v_negimply_inplace(H2, H1); /* H2 <- H2 - H1 */
310 : }
311 14 : return gerepileupto(av, zv_to_ZV(z));
312 : }
313 :
314 : /*************************************************************************/
315 : /** **/
316 : /** znstar/HNF interface **/
317 : /** **/
318 : /*************************************************************************/
319 :
320 : static long
321 2961 : mod_to_small(GEN x)
322 2961 : { return itos(typ(x) == t_INTMOD ? gel(x,2): x); }
323 :
324 : static GEN
325 2184 : vecmod_to_vecsmall(GEN x)
326 5145 : { 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 2184 : znstar_small(GEN zn)
331 : {
332 2184 : GEN Z = cgetg(4,t_VEC);
333 2184 : gel(Z,1) = icopy(gmael3(zn,3,1,1));
334 2184 : gel(Z,2) = vec_to_vecsmall(gel(zn,2));
335 2184 : 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 4200 : znstar_hnf_generators(GEN Z, GEN M)
341 : {
342 4200 : long j, h, l = lg(M);
343 4200 : GEN gen = cgetg(l, t_VECSMALL);
344 4200 : pari_sp ltop = avma;
345 4200 : GEN zgen = gel(Z,3);
346 4200 : ulong n = itou(gel(Z,1));
347 9177 : for (j = 1; j < l; j++)
348 : {
349 4977 : GEN Mj = gel(M,j);
350 4977 : gen[j] = 1;
351 12320 : for (h = 1; h < l; h++)
352 : {
353 7343 : ulong u = itou(gel(Mj,h));
354 7343 : if (!u) continue;
355 5404 : gen[j] = Fl_mul(uel(gen,j), Fl_powu(uel(zgen,h), u, n), n);
356 : }
357 : }
358 4200 : set_avma(ltop); return gen;
359 : }
360 :
361 : GEN
362 3780 : znstar_hnf(GEN Z, GEN M)
363 3780 : { return znstar_generate(itos(gel(Z,1)),znstar_hnf_generators(Z,M)); }
364 :
365 : GEN
366 3780 : znstar_hnf_elts(GEN Z, GEN H)
367 : {
368 3780 : pari_sp ltop = avma;
369 3780 : GEN G = znstar_hnf(Z,H);
370 3780 : long n = itos(gel(Z,1));
371 3780 : return gerepileupto(ltop, znstar_elts(n,G));
372 : }
373 :
374 : /*************************************************************************/
375 : /** **/
376 : /** polsubcyclo **/
377 : /** **/
378 : /*************************************************************************/
379 :
380 : static GEN
381 49813 : gscycloconductor(GEN g, long n, long flag)
382 : {
383 49813 : if (flag==2) retmkvec2(gcopy(g), stoi(n));
384 49806 : return g;
385 : }
386 :
387 : static long
388 1337848 : lift_check_modulus(GEN H, long n)
389 : {
390 : long h;
391 1337848 : switch(typ(H))
392 : {
393 105 : case t_INTMOD:
394 105 : if (!equalsi(n, gel(H,1)))
395 7 : pari_err_MODULUS("galoissubcyclo", stoi(n), gel(H,1));
396 98 : H = gel(H,2); /* fall through */
397 1337841 : case t_INT:
398 1337841 : h = smodis(H,n);
399 1337841 : if (ugcd(h,n) != 1) pari_err_COPRIME("galoissubcyclo", H,stoi(n));
400 1337834 : 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 5648028 : polsubcyclo_powz(GEN powz, long ex)
411 : {
412 5648028 : long m = lg(gel(powz,1))-1, q = ex/m, r = ex%m; /*ex=m*q+r*/
413 5648028 : GEN g = gmael(powz,1,r+1), G = gmael(powz,2,q+1);
414 5648028 : return (lg(powz)==4)? mulreal(g,G): gmul(g,G);
415 : }
416 :
417 : static GEN
418 81042 : polsubcyclo_complex_bound(pari_sp av, GEN V, long prec)
419 : {
420 81042 : GEN pol = real_i(roots_to_pol(V,0));
421 81042 : 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 62598 : polsubcyclo_cyclic(long n, long d, long m ,long z, long g, GEN powz, GEN le)
427 : {
428 62598 : GEN V = cgetg(d+1,t_VEC);
429 62598 : ulong base = 1;
430 : long i,k;
431 : pari_timer ti;
432 62598 : if (DEBUGLEVEL >= 6) timer_start(&ti);
433 303578 : for (i=1; i<=d; i++, base = Fl_mul(base,z,n))
434 : {
435 240980 : pari_sp av = avma;
436 240980 : long ex = base;
437 240980 : GEN s = gen_0;
438 1110272 : for (k=0; k<m; k++, ex = Fl_mul(ex,g,n))
439 : {
440 869292 : s = gadd(s, polsubcyclo_powz(powz,ex));
441 869292 : if ((k&0xff)==0) s = gerepileupto(av,s);
442 : }
443 240980 : if (le) s = modii(s, le);
444 240980 : gel(V,i) = gerepileupto(av, s);
445 : }
446 62598 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "polsubcyclo_cyclic");
447 62598 : 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 4778736 : _subcyclo_orbits(struct _subcyclo_orbits_s *data, long k)
499 : {
500 4778736 : GEN powz = data->powz;
501 4778736 : GEN *s = data->s;
502 :
503 4778736 : if (!data->count) data->ltop = avma;
504 4778736 : *s = gadd(*s, polsubcyclo_powz(powz,k));
505 4778736 : data->count++;
506 4778736 : if ((data->count & 0xffUL) == 0) *s = gerepileupto(data->ltop, *s);
507 4778736 : }
508 :
509 : /* Newton sums mod le. if le==NULL, works with complex instead */
510 : static GEN
511 99486 : polsubcyclo_orbits(long n, GEN H, GEN O, GEN powz, GEN le)
512 : {
513 99486 : long i, d = lg(O);
514 99486 : GEN V = cgetg(d,t_VEC);
515 : struct _subcyclo_orbits_s data;
516 99486 : long lle = le?lg(le)*2+1: 2*lg(gmael(powz,1,2))+3;/*dvmdii uses lx+ly space*/
517 99486 : data.powz = powz;
518 375264 : for(i=1; i<d; i++)
519 : {
520 275778 : GEN s = gen_0;
521 275778 : pari_sp av = avma;
522 275778 : (void)new_chunk(lle);
523 275778 : data.count = 0;
524 275778 : data.s = &s;
525 275778 : znstar_coset_func(n, H, (void (*)(void *,long)) _subcyclo_orbits,
526 275778 : (void *) &data, O[i]);
527 275778 : set_avma(av); /* HACK */
528 275778 : gel(V,i) = le? modii(s,le): gcopy(s);
529 : }
530 99486 : return V;
531 : }
532 :
533 : static GEN
534 81903 : 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 81903 : l = e*n+1;
540 311390 : while(!uisprime(l)) { l += n; e++; }
541 81903 : if (DEBUGLEVEL >= 4) err_printf("Subcyclo: prime l=%ld\n",l);
542 81903 : gl = utoipos(l); av = avma;
543 81903 : 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 81903 : if (DEBUGLEVEL >= 4) err_printf("Subcyclo: bound=2^%ld\n",expi(borne));
549 81903 : val = logint(shifti(borne,2), gl) + 1;
550 81903 : set_avma(av);
551 81903 : if (DEBUGLEVEL >= 4) err_printf("Subcyclo: val=%ld\n",val);
552 81903 : le = powiu(gl,val);
553 81903 : z = utoipos( Fl_powu(pgener_Fl(l), e, l) );
554 81903 : z = Zp_sqrtnlift(gen_1,utoipos(n),z,gl,val);
555 81903 : *ptr_val = val;
556 81903 : *ptr_l = l;
557 81903 : 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 81042 : polsubcyclo_complex_roots(long n, long real, long prec)
566 : {
567 81042 : long i, m = (long)(1+sqrt((double) n));
568 81042 : GEN bab, gig, powz = cgetg(real?4:3, t_VEC);
569 :
570 81042 : bab = cgetg(m+1,t_VEC);
571 81042 : gel(bab,1) = gen_1;
572 81042 : gel(bab,2) = rootsof1u_cx(n, prec); /* = e_n(1) */
573 374982 : for (i=3; i<=m; i++) gel(bab,i) = gmul(gel(bab,2),gel(bab,i-1));
574 81042 : gig = cgetg(m+1,t_VEC);
575 81042 : gel(gig,1) = gen_1;
576 81042 : gel(gig,2) = gmul(gel(bab,2),gel(bab,m));;
577 374982 : for (i=3; i<=m; i++) gel(gig,i) = gmul(gel(gig,2),gel(gig,i-1));
578 81042 : gel(powz,1) = bab;
579 81042 : gel(powz,2) = gig;
580 81042 : if (real) gel(powz,3) = gen_0;
581 81042 : return powz;
582 : }
583 :
584 : static GEN
585 668922 : muliimod_sz(GEN x, GEN y, GEN l, long siz)
586 : {
587 668922 : pari_sp av = avma;
588 : GEN p1;
589 668922 : (void)new_chunk(siz); /* HACK */
590 668922 : p1 = mulii(x,y);
591 668922 : set_avma(av); return modii(p1,l);
592 : }
593 :
594 : static GEN
595 81042 : polsubcyclo_roots(long n, GEN zl)
596 : {
597 81042 : GEN le = gel(zl,1), z = gel(zl,2);
598 81042 : long i, lle = lg(le)*3; /*Assume dvmdii use lx+ly space*/
599 81042 : long m = (long)(1+sqrt((double) n));
600 81042 : GEN bab, gig, powz = cgetg(3,t_VEC);
601 : pari_timer ti;
602 81042 : if (DEBUGLEVEL >= 6) timer_start(&ti);
603 81042 : bab = cgetg(m+1,t_VEC);
604 81042 : gel(bab,1) = gen_1;
605 81042 : gel(bab,2) = icopy(z);
606 374982 : for (i=3; i<=m; i++) gel(bab,i) = muliimod_sz(z,gel(bab,i-1),le,lle);
607 81042 : gig = cgetg(m+1,t_VEC);
608 81042 : gel(gig,1) = gen_1;
609 81042 : gel(gig,2) = muliimod_sz(z,gel(bab,m),le,lle);;
610 374982 : for (i=3; i<=m; i++) gel(gig,i) = muliimod_sz(gel(gig,2),gel(gig,i-1),le,lle);
611 81042 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "polsubcyclo_roots");
612 81042 : gel(powz,1) = bab;
613 81042 : gel(powz,2) = gig; return powz;
614 : }
615 :
616 : GEN
617 434 : galoiscyclo(long n, long v)
618 : {
619 434 : ulong av = avma;
620 : GEN grp, G, z, le, L, elts;
621 : long val, l, i, j, k;
622 434 : GEN zn = znstar(stoi(n));
623 434 : long card = itos(gel(zn,1));
624 434 : GEN gen = vec_to_vecsmall(lift_shallow(gel(zn,3)));
625 434 : GEN ord = vec_to_vecsmall(gel(zn,2));
626 434 : GEN T = polcyclo(n,v);
627 434 : long d = degpol(T);
628 434 : GEN borneabs = powuu(2,d);
629 434 : z = polsubcyclo_start(n,card/2,2,2*usqrt(d),borneabs,&val,&l);
630 434 : le = gel(z,1); z = gel(z,2);
631 434 : L = cgetg(1+card,t_VEC);
632 434 : gel(L,1) = z;
633 882 : for (j = 1, i = 1; j < lg(gen); j++)
634 : {
635 448 : long c = i * (ord[j]-1);
636 2156 : for (k = 1; k <= c; k++) gel(L,++i) = Fp_powu(gel(L,k), gen[j], le);
637 : }
638 434 : G = abelian_group(ord);
639 434 : elts = group_elts(G, card); /*not stack clean*/
640 434 : grp = cgetg(9, t_VEC);
641 434 : gel(grp,1) = T;
642 434 : gel(grp,2) = mkvec3(stoi(l), stoi(val), icopy(le));
643 434 : gel(grp,3) = L;
644 434 : gel(grp,4) = FpV_invVandermonde(L, NULL, le);
645 434 : gel(grp,5) = gen_1;
646 434 : gel(grp,6) = elts;
647 434 : gel(grp,7) = gel(G,1);
648 434 : gel(grp,8) = gel(G,2);
649 434 : 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 56 : bnr_to_abgrp(GEN bnr, long *cx)
657 : {
658 56 : GEN gen, F, v, bid, G, Ui = NULL;
659 : long l, i;
660 56 : checkbnr(bnr);
661 56 : bid = bnr_get_bid(bnr);
662 56 : G = bnr_get_clgp(bnr);
663 56 : if (lg(G) == 4)
664 28 : gen = abgrp_get_gen(G);
665 : else
666 : {
667 28 : Ui = gmael(bnr,4,3);
668 28 : if (ZM_isidentity(Ui)) Ui = NULL;
669 28 : gen = bid_get_gen(bid);
670 : }
671 56 : F = bid_get_ideal(bid);
672 56 : if (lg(F) != 2)
673 7 : 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 49 : F = gcoeff(F, 1, 1);
676 49 : *cx = signe(gel(bid_get_arch(bid), 1));
677 49 : l = lg(gen); v = cgetg(l, t_VEC);
678 203 : for (i = 1; i < l; ++i)
679 : {
680 154 : GEN x = gel(gen,i);
681 154 : if (typ(x) == t_COL) x = gel(x,1);
682 154 : gel(v,i) = gmodulo(absi_shallow(x), F);
683 : }
684 49 : if (Ui)
685 : { /* from bid.gen to bnr.gen (maybe one less) */
686 21 : GEN w = v;
687 21 : l = lg(Ui); v = cgetg(l, t_VEC);
688 84 : for (i = 1; i < l; i++) gel(v,i) = factorback2(w, gel(Ui, i));
689 : }
690 49 : return mkvec3(bnr_get_no(bnr), bnr_get_cyc(bnr), v);
691 : }
692 :
693 : static long
694 50247 : _itos(const char *fun, GEN n)
695 : {
696 50247 : if (is_bigint(n))
697 7 : pari_err_IMPL(stack_sprintf("conductor f > %ld in %s", LONG_MAX, fun));
698 50240 : return itos(n);
699 : }
700 : long
701 50261 : subcyclo_nH(const char *fun, GEN N, GEN *psg)
702 : {
703 50261 : GEN V, Z = NULL, H = *psg;
704 50261 : long i, l, n = 0, complex = 1;
705 50261 : switch(typ(N))
706 : {
707 49820 : case t_INT:
708 49820 : n = _itos(fun, N);
709 49813 : if (n < 1) pari_err_DOMAIN(fun, "degree", "<=", gen_0, N);
710 49806 : break;
711 441 : case t_VEC:
712 441 : if (lg(N)==7)
713 56 : N = bnr_to_abgrp(N,&complex);
714 385 : else if (checkznstar_i(N))
715 14 : N = mkvec3(znstar_get_no(N), znstar_get_cyc(N),
716 : gmodulo(znstar_get_gen(N), znstar_get_N(N)));
717 434 : if (lg(N)==4)
718 : { /* abgrp */
719 434 : GEN gen = abgrp_get_gen(N), z;
720 434 : if (typ(gen)!=t_VEC) pari_err_TYPE(fun,gen);
721 434 : Z = N;
722 434 : if (lg(gen) == 1) { n = 1; break; }
723 434 : z = gel(gen,1);
724 434 : if (typ(z) == t_INTMOD) { n = _itos(fun, gel(z,1)); break; }
725 : }
726 : default: /*fall through*/
727 7 : pari_err_TYPE(fun,N);
728 : return 0;/*LCOV_EXCL_LINE*/
729 : }
730 50233 : if (!H) H = gen_1;
731 50233 : switch(typ(H))
732 : {
733 49624 : case t_INTMOD: case t_INT:
734 49624 : V = mkvecsmall( lift_check_modulus(H,n) );
735 49617 : break;
736 14 : case t_VECSMALL:
737 14 : l = lg(H); V = leafcopy(H);
738 42 : for (i = 1; i < l; i++) { V[i] %= n; if (V[i] < 0) V[i] += n; }
739 14 : break;
740 154 : case t_VEC: case t_COL:
741 154 : l = lg(H); V = cgetg(l,t_VECSMALL);
742 1288371 : for(i=1; i < l; i++) V[i] = lift_check_modulus(gel(H,i),n);
743 147 : break;
744 441 : case t_MAT:
745 441 : l = lg(H);
746 441 : if (l == 1 || l != lgcols(H))
747 7 : pari_err_TYPE(stack_strcat(fun," [H not in HNF]"),H);
748 434 : if (!Z) pari_err_TYPE(stack_strcat(fun," [N not a bnrinit or znstar]"),H);
749 427 : if (lg(gel(Z,2)) != l) pari_err_DIM(fun);
750 420 : V = znstar_hnf_generators(znstar_small(Z),H);
751 420 : break;
752 0 : default:
753 0 : pari_err_TYPE(fun,H);
754 : return 0;/*LCOV_EXCL_LINE*/
755 : }
756 50198 : if (!complex) V = vecsmall_append(V, n-1); /*add complex conjugation*/
757 50198 : *psg = V; return n;
758 : }
759 :
760 : GEN
761 49876 : galoissubcyclo(GEN N, GEN sg, long flag, long v)
762 : {
763 49876 : 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 49876 : if (flag<0 || flag>3) pari_err_FLAG("galoissubcyclo");
769 49876 : if (v < 0) v = 0;
770 49876 : n = subcyclo_nH("galoissubcyclo", N, &sg);
771 49827 : if (n==1)
772 : {
773 21 : set_avma(ltop); if (flag == 1) return gen_1;
774 14 : return gscycloconductor(deg1pol_shallow(gen_1, gen_m1, v), 1, flag);
775 : }
776 49806 : H = znstar_generate(n, sg);
777 49806 : 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 49806 : complex = !F2v_coeff(gel(H,3),n-1);
786 49806 : if (DEBUGLEVEL >= 6) err_printf("Subcyclo: complex=%ld\n",complex);
787 49806 : if (DEBUGLEVEL >= 1) timer_start(&ti);
788 49806 : cnd = znstar_conductor(H);
789 49806 : if (DEBUGLEVEL >= 1) timer_printf(&ti, "znstar_conductor");
790 49806 : if (flag == 1) return gc_stoi(ltop, cnd);
791 49806 : if (cnd == 1)
792 : {
793 63 : set_avma(ltop); if (flag == 1) return gen_1;
794 63 : return gscycloconductor(deg1pol_shallow(gen_1,gen_m1,v),1,flag);
795 : }
796 49743 : if (n != cnd)
797 : {
798 175 : H = znstar_reduce_modulus(H, cnd);
799 175 : n = cnd;
800 : }
801 49743 : card = znstar_order(H);
802 49743 : phi_n = eulerphiu(n);
803 49743 : if (card == phi_n)
804 : {
805 0 : set_avma(ltop);
806 0 : return gscycloconductor(polcyclo(n,v),n,flag);
807 : }
808 49743 : O = znstar_cosets(n, phi_n, H);
809 49743 : if (DEBUGLEVEL >= 1) timer_printf(&ti, "znstar_cosets");
810 49743 : if (DEBUGLEVEL >= 6) err_printf("Subcyclo: orbits=%Ps\n",O);
811 49743 : if (DEBUGLEVEL >= 4)
812 0 : err_printf("Subcyclo: %ld orbits with %ld elements each\n",phi_n/card,card);
813 49743 : av = avma;
814 49743 : powz = polsubcyclo_complex_roots(n,!complex,LOWDEFAULTPREC);
815 49743 : L = polsubcyclo_orbits(n,H,O,powz,NULL);
816 49743 : B = polsubcyclo_complex_bound(av,L,LOWDEFAULTPREC);
817 49743 : zl = polsubcyclo_start(n,phi_n/card,card,1,B,&val,&l);
818 49743 : powz = polsubcyclo_roots(n,zl);
819 49743 : le = gel(zl,1);
820 49743 : L = polsubcyclo_orbits(n,H,O,powz,le);
821 49743 : if (DEBUGLEVEL >= 6) timer_start(&ti);
822 49743 : T = FpV_roots_to_pol(L,le,v);
823 49743 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "roots_to_pol");
824 49743 : T = FpX_center(T,le,shifti(le,-1));
825 49743 : if (flag==3)
826 : {
827 : GEN L2, aut;
828 7 : if (Flx_is_squarefree(ZX_to_Flx(T, l),l))
829 7 : 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 7 : aut = znstar_quotient(n, phi_n, H, L2, l);
840 7 : return gerepileupto(ltop, galoisinitfromaut(T, aut, l));
841 : }
842 49736 : 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 33527 : polsubcyclo_g(long n, long d, GEN Z, long v)
849 : {
850 33527 : 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 33527 : if (d==1) return deg1pol_shallow(gen_1,gen_m1,v); /* get rid of n=1,2 */
855 33527 : if ((n & 3) == 2) n >>= 1;
856 : /* n = 4 or p^a, p odd */
857 33527 : o = itos(gel(Z,1));
858 33527 : g = itos(gmael3(Z,3,1,2));
859 33527 : p = n / ugcd(n,o); /* p^a / gcd(p^a,phi(p^a)) = p*/
860 33527 : r = ugcd(d,n); /* = p^(v_p(d)) < n */
861 33527 : n = r*p; /* n is now the conductor */
862 33527 : o = n-r; /* = phi(n) */
863 33527 : if (o == d) return polcyclo(n,v);
864 31299 : o /= d;
865 31299 : gd = Fl_powu(g%n, d, n);
866 : /*FIXME: If degree is small, the computation of B is a waste of time*/
867 31299 : powz = polsubcyclo_complex_roots(n,(o&1)==0,LOWDEFAULTPREC);
868 31299 : L = polsubcyclo_cyclic(n,d,o,g,gd,powz,NULL);
869 31299 : B = polsubcyclo_complex_bound(ltop,L,LOWDEFAULTPREC);
870 31299 : zl = polsubcyclo_start(n,d,o,1,B,&val,&l);
871 31299 : le = gel(zl,1);
872 31299 : powz = polsubcyclo_roots(n,zl);
873 31299 : L = polsubcyclo_cyclic(n,d,o,g,gd,powz,le);
874 31299 : if (DEBUGLEVEL >= 6) timer_start(&ti);
875 31299 : T = FpV_roots_to_pol(L,le,v);
876 31299 : if (DEBUGLEVEL >= 6) timer_printf(&ti, "roots_to_pol");
877 31299 : return gerepileupto(ltop, FpX_center(T,le,shifti(le,-1)));
878 : }
879 :
880 : GEN
881 33583 : polsubcyclo(long n, long d, long v)
882 : {
883 33583 : pari_sp ltop = avma;
884 : GEN L, Z;
885 33583 : if (v<0) v = 0;
886 33583 : if (d<=0) pari_err_DOMAIN("polsubcyclo","d","<=",gen_0,stoi(d));
887 33576 : if (n<=0) pari_err_DOMAIN("polsubcyclo","n","<=",gen_0,stoi(n));
888 33569 : Z = znstar(stoi(n));
889 33569 : if (!dvdis(gel(Z,1), d)) { set_avma(ltop); return cgetg(1, t_VEC); }
890 33562 : if (lg(gel(Z,2)) == 2)
891 : { /* faster but Z must be cyclic */
892 33527 : set_avma(ltop);
893 33527 : return polsubcyclo_g(n, d, Z, v);
894 : }
895 35 : L = subgrouplist(gel(Z,2), mkvec(stoi(d)));
896 35 : if (lg(L) == 2)
897 7 : return gerepileupto(ltop, galoissubcyclo(Z, gel(L,1), 0, v));
898 : else
899 : {
900 28 : GEN V = cgetg(lg(L),t_VEC);
901 : long i;
902 364 : for (i=1; i< lg(V); i++) gel(V,i) = galoissubcyclo(Z, gel(L,i), 0, v);
903 28 : 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 427 : factor_Aurifeuille_aux(GEN A, long Astar, long n, GEN P,
937 : struct aurifeuille_t *S)
938 : {
939 : pari_sp av;
940 427 : GEN f, a, b, s, powers, z = S->z, le = S->le;
941 427 : long j, k, maxjump, lastj, e = S->e;
942 427 : ulong l = S->l;
943 : char *invertible;
944 :
945 427 : if ((n & 7) == 4)
946 : { /* A^* even */
947 350 : GEN i = Fp_powu(z, n>>2, le), z2 = Fp_sqr(z, le);
948 :
949 350 : invertible = stack_malloc(n); /* even indices unused */
950 1610 : for (j = 1; j < n; j+=2) invertible[j] = 1;
951 406 : for (k = 1; k < lg(P); k++)
952 : {
953 56 : long p = P[k];
954 168 : for (j = p; j < n; j += 2*p) invertible[j] = 0;
955 : }
956 350 : lastj = 1; maxjump = 2;
957 1260 : for (j= 3; j < n; j+=2)
958 910 : if (invertible[j]) {
959 798 : long jump = j - lastj;
960 798 : if (jump > maxjump) maxjump = jump;
961 798 : lastj = j;
962 : }
963 350 : powers = cgetg(maxjump+1, t_VEC); /* powers[k] = z^k, odd indices unused */
964 350 : gel(powers,2) = z2;
965 406 : for (k = 4; k <= maxjump; k+=2)
966 112 : gel(powers,k) = odd(k>>1)? Fp_mul(gel(powers, k-2), z2, le)
967 56 : : Fp_sqr(gel(powers, k>>1), le);
968 :
969 350 : if (Astar == 2)
970 : { /* important special case (includes A=2), split for efficiency */
971 329 : if (!equalis(A, 2))
972 : {
973 14 : GEN f = sqrti(shifti(A,-1)), mf = Fp_neg(f,le), fi = Fp_mul(f,i,le);
974 14 : a = Fp_add(mf, fi, le);
975 14 : b = Fp_sub(mf, fi, le);
976 : }
977 : else
978 : {
979 315 : a = subiu(i,1);
980 315 : b = subsi(-1,i);
981 : }
982 329 : av = avma;
983 329 : s = z; f = subii(a, s); lastj = 1;
984 910 : for (j = 3, k = 0; j < n; j+=2)
985 581 : if (invertible[j])
986 : {
987 511 : s = Fp_mul(gel(powers, j-lastj), s, le); /* z^j */
988 511 : lastj = j;
989 511 : f = Fp_mul(f, subii((j & 3) == 1? a: b, s), le);
990 511 : if (++k == 0x1ff) { gerepileall(av, 2, &s, &f); k = 0; }
991 : }
992 : }
993 : else
994 : {
995 21 : GEN ma, mb, B = Fp_mul(A, i, le), gl = utoipos(l);
996 : long t;
997 21 : Astar >>= 1;
998 21 : t = Astar & 3; if (Astar < 0) t = 4-t; /* t = 1 or 3 */
999 21 : if (t == 1) B = Fp_neg(B, le);
1000 21 : a = Zp_sqrtlift(B, Fp_sqrt(B, gl), gl, e);
1001 21 : b = Fp_mul(a, i, le);
1002 21 : ma = Fp_neg(a, le);
1003 21 : mb = Fp_neg(b, le);
1004 21 : av = avma;
1005 21 : s = z; f = subii(a, s); lastj = 1;
1006 350 : for (j = 3, k = 0; j<n; j+=2)
1007 329 : if (invertible[j])
1008 : {
1009 : GEN t;
1010 287 : if ((j & 3) == 1) t = (kross(j, Astar) < 0)? ma: a;
1011 154 : else t = (kross(j, Astar) < 0)? mb: b;
1012 287 : s = Fp_mul(gel(powers, j-lastj), s, le); /* z^j */
1013 287 : lastj = j;
1014 287 : f = Fp_mul(f, subii(t, s), le);
1015 287 : if (++k == 0x1ff) { gerepileall(av, 2, &s, &f); k = 0; }
1016 : }
1017 : }
1018 : }
1019 : else /* A^* odd */
1020 : {
1021 : ulong g;
1022 77 : 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 77 : g = Fl_sqrt(umodiu(A,l), l);
1030 77 : a = Zp_sqrtlift(A, utoipos(g), utoipos(l), e);
1031 77 : b = negi(a);
1032 :
1033 77 : invertible = stack_malloc(n);
1034 1337 : for (j = 1; j < n; j++) invertible[j] = 1;
1035 196 : for (k = 1; k < lg(P); k++)
1036 : {
1037 119 : long p = P[k];
1038 483 : for (j = p; j < n; j += p) invertible[j] = 0;
1039 : }
1040 77 : lastj = 2; maxjump = 1;
1041 1183 : for (j= 3; j < n; j++)
1042 1106 : if (invertible[j]) {
1043 742 : long jump = j - lastj;
1044 742 : if (jump > maxjump) maxjump = jump;
1045 742 : lastj = j;
1046 : }
1047 77 : powers = cgetg(maxjump+1, t_VEC); /* powers[k] = z^k */
1048 77 : gel(powers,1) = z;
1049 161 : for (k = 2; k <= maxjump; k++)
1050 168 : gel(powers,k) = odd(k)? Fp_mul(gel(powers, k-1), z, le)
1051 84 : : Fp_sqr(gel(powers, k>>1), le);
1052 77 : av = avma;
1053 77 : s = z; f = subii(a, s); lastj = 1;
1054 1260 : for(j = 2, k = 0; j < n; j++)
1055 1183 : if (invertible[j])
1056 : {
1057 819 : s = Fp_mul(gel(powers, j-lastj), s, le);
1058 819 : lastj = j;
1059 819 : f = Fp_mul(f, subii(kross(j,Astar)==1? a: b, s), le);
1060 819 : if (++k == 0x1ff) { gerepileall(av, 2, &s, &f); k = 0; }
1061 : }
1062 : }
1063 427 : return f;
1064 : }
1065 :
1066 : /* d != 2 mod 4; fd = factoru(odd(d)? d: d / 4) */
1067 : static void
1068 427 : Aurifeuille_init(GEN a, long d, GEN fd, struct aurifeuille_t *S)
1069 : {
1070 427 : GEN bound, zl, sqrta = sqrtr_abs(itor(a, LOWDEFAULTPREC));
1071 427 : ulong phi = eulerphiu_fact(fd);
1072 427 : if (!odd(d)) phi <<= 1; /* eulerphi(d) */
1073 427 : bound = ceil_safe(powru(addrs(sqrta,1), phi));
1074 427 : zl = polsubcyclo_start(d, 0, 0, 1, bound, &(S->e), (long*)&(S->l));
1075 427 : S->le = gel(zl,1);
1076 427 : S->z = gel(zl,2);
1077 427 : }
1078 :
1079 : GEN
1080 364 : factor_Aurifeuille_prime(GEN p, long d)
1081 : {
1082 364 : pari_sp av = avma;
1083 : struct aurifeuille_t S;
1084 : GEN fd;
1085 : long pp;
1086 364 : if ((d & 3) == 2) { d >>= 1; p = negi(p); }
1087 364 : fd = factoru(odd(d)? d: d>>2);
1088 364 : pp = itos(p);
1089 364 : Aurifeuille_init(p, d, fd, &S);
1090 364 : 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 63 : factor_Aurifeuille(GEN a, long d)
1096 : {
1097 63 : pari_sp av = avma;
1098 : GEN fd, P, A;
1099 63 : long i, lP, va = vali(a), sa, astar, D;
1100 : struct aurifeuille_t S;
1101 :
1102 63 : if (d <= 0)
1103 0 : pari_err_DOMAIN("factor_Aurifeuille", "degre", "<=",gen_0,stoi(d));
1104 63 : if ((d & 3) == 2) { d >>= 1; a = negi(a); }
1105 63 : if ((va & 1) == (d & 1)) { set_avma(av); return gen_1; }
1106 63 : sa = signe(a);
1107 63 : if (odd(d))
1108 : {
1109 : long a4;
1110 28 : if (d == 1)
1111 : {
1112 0 : if (!Z_issquareall(a, &A)) return gen_1;
1113 0 : return gerepileuptoint(av, addiu(A,1));
1114 : }
1115 28 : A = va? shifti(a, -va): a;
1116 28 : a4 = mod4(A); if (sa < 0) a4 = 4 - a4;
1117 28 : if (a4 != 1) { set_avma(av); return gen_1; }
1118 : }
1119 35 : else if ((d & 7) == 4)
1120 35 : 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 63 : fd = factoru(odd(d)? d: d>>2); P = gel(fd,1); lP = lg(P);
1127 63 : astar = sa;
1128 63 : if (odd(va)) astar <<= 1;
1129 147 : for (i = 1; i < lP; i++)
1130 84 : if (odd( (Z_lvalrem(A, P[i], &A)) ) ) astar *= P[i];
1131 63 : if (sa < 0)
1132 : { /* negate in place if possible */
1133 14 : if (A == a) A = icopy(A);
1134 14 : setabssign(A);
1135 : }
1136 63 : if (!Z_issquare(A)) { set_avma(av); return gen_1; }
1137 :
1138 63 : D = odd(d)? 1: 4;
1139 147 : for (i = 1; i < lP; i++) D *= P[i];
1140 63 : if (D != d) { a = powiu(a, d/D); d = D; }
1141 :
1142 63 : Aurifeuille_init(a, d, fd, &S);
1143 63 : return gerepileuptoint(av, factor_Aurifeuille_aux(a, astar, d, P, &S));
1144 : }
|
