Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - subcyclo.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.2 lcov report (development 29115-f22e516b23) Lines: 608 656 92.7 %
Date: 2024-03-19 08:06:49 Functions: 46 49 93.9 %
Legend: Lines: hit not hit

          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     1349177 : 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     1349177 :   if (!d) { (*func)(data,c); return; }
      40             : 
      41     1093812 :   cache = const_vecsmall(d,c);
      42     1093812 :   (*func)(data,c);  /* AFTER cache: may contain gerepileupto statement */
      43     1093812 :   gen = gel(H,1);
      44     1093812 :   ord = gel(H,2);
      45     1154663 :   card = ord[1]; for (i = 2; i <= d; i++) card *= ord[i];
      46    82530657 :   for(i=1; i<card; i++)
      47             :   {
      48    81436845 :     long k, m = i;
      49    82010173 :     for(j=1; j<d && m%ord[j]==0 ;j++) m /= ord[j];
      50    81436845 :     cache[j] = Fl_mul(cache[j],gen[j],n);
      51    82010173 :     for (k=1; k<j; k++) cache[k] = cache[j];
      52    81436845 :     (*func)(data, cache[j]);
      53             :   }
      54             : }
      55             : 
      56             : static void
      57      275736 : znstar_coset_func(long n, GEN H, void (*func)(void *data,long c)
      58             :     , void *data, long c)
      59      275736 : { 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     1073441 : znstar_partial_coset_bits_inplace(long n, GEN H, GEN bits, long d, long c)
      66             : {
      67     1073441 :   pari_sp av = avma;
      68     1073441 :   znstar_partial_coset_func(n,H, (void (*)(void *,long)) &F2v_set,
      69             :       (void *) bits, d, c);
      70     1073441 :   set_avma(av);
      71     1073441 : }
      72             : 
      73             : static void
      74      531834 : znstar_coset_bits_inplace(long n, GEN H, GEN bits, long c)
      75      531834 : { znstar_partial_coset_bits_inplace(n, H, bits, lg(gel(H,1))-1, c); }
      76             : 
      77             : static GEN
      78      541607 : znstar_partial_coset_bits(long n, GEN H, long d, long c)
      79             : {
      80      541607 :   GEN bits = zero_F2v(n);
      81      541607 :   znstar_partial_coset_bits_inplace(n,H,bits,d,c);
      82      541607 :   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      541607 : znstar_partial_bits(long n, GEN H, long d)
      89      541607 : { 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      254987 : znstar_generate(long n, GEN V)
     100             : {
     101      254987 :   pari_sp av = avma;
     102      254987 :   GEN gen = cgetg(lg(V),t_VECSMALL);
     103      254987 :   GEN ord = cgetg(lg(V),t_VECSMALL), res = mkvec2(gen,ord);
     104      254987 :   GEN bits = znstar_partial_bits(n,NULL,0);
     105      254987 :   long i, r = 0;
     106     1850824 :   for(i=1; i<lg(V); i++)
     107             :   {
     108     1595837 :     ulong v = uel(V,i), g = v;
     109     1595837 :     long o = 0;
     110    38258191 :     while (!F2v_coeff(bits, (long)g)) { g = Fl_mul(g, v, (ulong)n); o++; }
     111     1595837 :     if (!o) continue;
     112      286620 :     r++;
     113      286620 :     gen[r] = v;
     114      286620 :     ord[r] = o+1;
     115      286620 :     cgiv(bits); bits = znstar_partial_bits(n,res,r);
     116             :   }
     117      254987 :   setlg(gen,r+1);
     118      254987 :   setlg(ord,r+1); return gerepilecopy(av, mkvec3(gen,ord,bits));
     119             : }
     120             : 
     121             : static ulong
     122      164016 : 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        3766 : znstar_elts(long n, GEN H)
     128             : {
     129        3766 :   long card = znstar_order(H);
     130        3766 :   GEN gen = gel(H,1), ord = gel(H,2);
     131        3766 :   GEN sg = cgetg(1 + card, t_VECSMALL);
     132             :   long k, j, l;
     133        3766 :   sg[1] = 1;
     134        5775 :   for (j = 1, l = 1; j < lg(gen); j++)
     135             :   {
     136        2009 :     long c = l * (ord[j]-1);
     137        4564 :     for (k = 1; k <= c; k++) sg[++l] = Fl_mul(sg[k], gen[j], n);
     138             :   }
     139        3766 :   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      172103 : znstar_conductor_bits(GEN bits)
     158             : {
     159      172103 :   pari_sp av = avma;
     160      172103 :   long i, f = 1, cnd0 = bits[1];
     161      172103 :   GEN F = factoru(cnd0), P = gel(F,1), E = gel(F,2);
     162      462897 :   for (i = lg(P)-1; i > 0; i--)
     163             :   {
     164      290794 :     long p = P[i], e = E[i], cnd = cnd0;
     165      328237 :     for (  ; e >= 2; e--)
     166             :     {
     167       83872 :       long q = cnd / p;
     168       83872 :       if (!F2v_coeff(bits, 1 + q)) break;
     169       37443 :       cnd = q;
     170             :     }
     171      290794 :     if (e == 1)
     172             :     {
     173      244365 :       if (p == 2) e = 0;
     174             :       else
     175             :       {
     176      215679 :         long h, g = pgener_Fl(p), q = cnd / p;
     177      215679 :         h = Fl_mul(g-1, Fl_inv(q % p, p), p); /* 1+h*q = g (mod p) */
     178      215679 :         if (F2v_coeff(bits, 1 + h*q)) e = 0;
     179             :       }
     180             :     }
     181      290794 :     if (e) f *= upowuu(p, e);
     182             :   }
     183      172103 :   return gc_long(av,f);
     184             : }
     185             : long
     186      171452 : 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      110507 : znstar_cosets(long n, long phi_n, GEN H)
     193             : {
     194      110507 :   long k, c = 0, card = znstar_order(H), index = phi_n/card;
     195      110507 :   GEN cosets = cgetg(index+1,t_VECSMALL);
     196      110507 :   pari_sp ltop = avma;
     197      110507 :   GEN bits = zero_F2v(n);
     198      642299 :   for (k = 1; k <= index; k++)
     199             :   {
     200     1371921 :     for (c++ ; F2v_coeff(bits,c) || ugcd(c,n)!=1; c++);
     201      531792 :     cosets[k]=c;
     202      531792 :     znstar_coset_bits_inplace(n, H, bits, c);
     203             :   }
     204      110507 :   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        2954 : mod_to_small(GEN x)
     322        2954 : { return itos(typ(x) == t_INTMOD ? gel(x,2): x); }
     323             : 
     324             : static GEN
     325        2177 : vecmod_to_vecsmall(GEN x)
     326        5131 : { 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        2177 : znstar_small(GEN zn)
     331             : {
     332        2177 :   GEN Z = cgetg(4,t_VEC);
     333        2177 :   gel(Z,1) = icopy(gmael3(zn,3,1,1));
     334        2177 :   gel(Z,2) = vec_to_vecsmall(gel(zn,2));
     335        2177 :   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        4186 : znstar_hnf_generators(GEN Z, GEN M)
     341             : {
     342        4186 :   long j, h, l = lg(M);
     343        4186 :   GEN gen = cgetg(l, t_VECSMALL);
     344        4186 :   pari_sp ltop = avma;
     345        4186 :   GEN zgen = gel(Z,3);
     346        4186 :   ulong n = itou(gel(Z,1));
     347        9149 :   for (j = 1; j < l; j++)
     348             :   {
     349        4963 :     GEN Mj = gel(M,j);
     350        4963 :     gen[j] = 1;
     351       12292 :     for (h = 1; h < l; h++)
     352             :     {
     353        7329 :       ulong u = itou(gel(Mj,h));
     354        7329 :       if (!u) continue;
     355        5390 :       gen[j] = Fl_mul(uel(gen,j), Fl_powu(uel(zgen,h), u, n), n);
     356             :     }
     357             :   }
     358        4186 :   set_avma(ltop); return gen;
     359             : }
     360             : 
     361             : GEN
     362        3766 : znstar_hnf(GEN Z, GEN M)
     363        3766 : { return znstar_generate(itos(gel(Z,1)),znstar_hnf_generators(Z,M)); }
     364             : 
     365             : GEN
     366        3766 : znstar_hnf_elts(GEN Z, GEN H)
     367             : {
     368        3766 :   pari_sp ltop = avma;
     369        3766 :   GEN G = znstar_hnf(Z,H);
     370        3766 :   long n = itos(gel(Z,1));
     371        3766 :   return gerepileupto(ltop, znstar_elts(n,G));
     372             : }
     373             : 
     374             : /*************************************************************************/
     375             : /**                                                                     **/
     376             : /**                     polsubcyclo                                     **/
     377             : /**                                                                     **/
     378             : /*************************************************************************/
     379             : 
     380             : static GEN
     381       49778 : gscycloconductor(GEN g, long n, long flag)
     382             : {
     383       49778 :   if (flag==2) retmkvec2(gcopy(g), stoi(n));
     384       49771 :   return g;
     385             : }
     386             : 
     387             : static long
     388     1337813 : lift_check_modulus(GEN H, long n)
     389             : {
     390             :   long h;
     391     1337813 :   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     1337806 :     case t_INT:
     398     1337806 :       h = smodis(H,n);
     399     1337806 :       if (ugcd(h,n) != 1) pari_err_COPRIME("galoissubcyclo", H,stoi(n));
     400     1337799 :       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     5647160 : polsubcyclo_powz(GEN powz, long ex)
     411             : {
     412     5647160 :   long m = lg(gel(powz,1))-1, q = ex/m, r = ex%m; /*ex=m*q+r*/
     413     5647160 :   GEN g = gmael(powz,1,r+1), G = gmael(powz,2,q+1);
     414     5647160 :   return (lg(powz)==4)? mulreal(g,G): gmul(g,G);
     415             : }
     416             : 
     417             : static GEN
     418       80986 : polsubcyclo_complex_bound(pari_sp av, GEN V, long prec)
     419             : {
     420       80986 :   GEN pol = real_i(roots_to_pol(V,0));
     421       80986 :   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       62500 : polsubcyclo_cyclic(long n, long d, long m ,long z, long g, GEN powz, GEN le)
     427             : {
     428       62500 :   GEN V = cgetg(d+1,t_VEC);
     429       62500 :   ulong base = 1;
     430             :   long i,k;
     431             :   pari_timer ti;
     432       62500 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
     433      303228 :   for (i=1; i<=d; i++, base = Fl_mul(base,z,n))
     434             :   {
     435      240728 :     pari_sp av = avma;
     436      240728 :     long ex = base;
     437      240728 :     GEN s = gen_0;
     438     1109236 :     for (k=0; k<m; k++, ex = Fl_mul(ex,g,n))
     439             :     {
     440      868508 :       s = gadd(s, polsubcyclo_powz(powz,ex));
     441      868508 :       if ((k&0xff)==0) s = gerepileupto(av,s);
     442             :     }
     443      240728 :     if (le) s = modii(s, le);
     444      240728 :     gel(V,i) = gerepileupto(av, s);
     445             :   }
     446       62500 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "polsubcyclo_cyclic");
     447       62500 :   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     4778652 : _subcyclo_orbits(struct _subcyclo_orbits_s *data, long k)
     499             : {
     500     4778652 :   GEN powz = data->powz;
     501     4778652 :   GEN *s = data->s;
     502             : 
     503     4778652 :   if (!data->count) data->ltop = avma;
     504     4778652 :   *s = gadd(*s, polsubcyclo_powz(powz,k));
     505     4778652 :   data->count++;
     506     4778652 :   if ((data->count & 0xffUL) == 0) *s = gerepileupto(data->ltop, *s);
     507     4778652 : }
     508             : 
     509             : /* Newton sums mod le. if le==NULL, works with complex instead */
     510             : static GEN
     511       99472 : polsubcyclo_orbits(long n, GEN H, GEN O, GEN powz, GEN le)
     512             : {
     513       99472 :   long i, d = lg(O);
     514       99472 :   GEN V = cgetg(d,t_VEC);
     515             :   struct _subcyclo_orbits_s data;
     516       99472 :   long lle = le?lg(le)*2+1: 2*lg(gmael(powz,1,2))+3;/*dvmdii uses lx+ly space*/
     517       99472 :   data.powz = powz;
     518      375208 :   for(i=1; i<d; i++)
     519             :   {
     520      275736 :     GEN s = gen_0;
     521      275736 :     pari_sp av = avma;
     522      275736 :     (void)new_chunk(lle);
     523      275736 :     data.count = 0;
     524      275736 :     data.s     = &s;
     525      275736 :     znstar_coset_func(n, H, (void (*)(void *,long)) _subcyclo_orbits,
     526      275736 :       (void *) &data, O[i]);
     527      275736 :     set_avma(av); /* HACK */
     528      275736 :     gel(V,i) = le? modii(s,le): gcopy(s);
     529             :   }
     530       99472 :   return V;
     531             : }
     532             : 
     533             : static GEN
     534       81819 : 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       81819 :   l = e*n+1;
     540      311166 :   while(!uisprime(l)) { l += n; e++; }
     541       81819 :   if (DEBUGLEVEL >= 4) err_printf("Subcyclo: prime l=%ld\n",l);
     542       81819 :   gl = utoipos(l); av = avma;
     543       81819 :   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       81819 :   if (DEBUGLEVEL >= 4) err_printf("Subcyclo: bound=2^%ld\n",expi(borne));
     549       81819 :   val = logint(shifti(borne,2), gl) + 1;
     550       81819 :   set_avma(av);
     551       81819 :   if (DEBUGLEVEL >= 4) err_printf("Subcyclo: val=%ld\n",val);
     552       81819 :   le = powiu(gl,val);
     553       81819 :   z = utoipos( Fl_powu(pgener_Fl(l), e, l) );
     554       81819 :   z = Zp_sqrtnlift(gen_1,utoipos(n),z,gl,val);
     555       81819 :   *ptr_val = val;
     556       81819 :   *ptr_l = l;
     557       81819 :   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       80986 : polsubcyclo_complex_roots(long n, long real, long prec)
     566             : {
     567       80986 :   long i, m = (long)(1+sqrt((double) n));
     568       80986 :   GEN bab, gig, powz = cgetg(real?4:3, t_VEC);
     569             : 
     570       80986 :   bab = cgetg(m+1,t_VEC);
     571       80986 :   gel(bab,1) = gen_1;
     572       80986 :   gel(bab,2) = rootsof1u_cx(n, prec); /* = e_n(1) */
     573      374849 :   for (i=3; i<=m; i++) gel(bab,i) = gmul(gel(bab,2),gel(bab,i-1));
     574       80986 :   gig = cgetg(m+1,t_VEC);
     575       80986 :   gel(gig,1) = gen_1;
     576       80986 :   gel(gig,2) = gmul(gel(bab,2),gel(bab,m));;
     577      374849 :   for (i=3; i<=m; i++) gel(gig,i) = gmul(gel(gig,2),gel(gig,i-1));
     578       80986 :   gel(powz,1) = bab;
     579       80986 :   gel(powz,2) = gig;
     580       80986 :   if (real) gel(powz,3) = gen_0;
     581       80986 :   return powz;
     582             : }
     583             : 
     584             : static GEN
     585      668712 : muliimod_sz(GEN x, GEN y, GEN l, long siz)
     586             : {
     587      668712 :   pari_sp av = avma;
     588             :   GEN p1;
     589      668712 :   (void)new_chunk(siz); /* HACK */
     590      668712 :   p1 = mulii(x,y);
     591      668712 :   set_avma(av); return modii(p1,l);
     592             : }
     593             : 
     594             : static GEN
     595       80986 : polsubcyclo_roots(long n, GEN zl)
     596             : {
     597       80986 :   GEN le = gel(zl,1), z = gel(zl,2);
     598       80986 :   long i, lle = lg(le)*3; /*Assume dvmdii use lx+ly space*/
     599       80986 :   long m = (long)(1+sqrt((double) n));
     600       80986 :   GEN bab, gig, powz = cgetg(3,t_VEC);
     601             :   pari_timer ti;
     602       80986 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
     603       80986 :   bab = cgetg(m+1,t_VEC);
     604       80986 :   gel(bab,1) = gen_1;
     605       80986 :   gel(bab,2) = icopy(z);
     606      374849 :   for (i=3; i<=m; i++) gel(bab,i) = muliimod_sz(z,gel(bab,i-1),le,lle);
     607       80986 :   gig = cgetg(m+1,t_VEC);
     608       80986 :   gel(gig,1) = gen_1;
     609       80986 :   gel(gig,2) = muliimod_sz(z,gel(bab,m),le,lle);;
     610      374849 :   for (i=3; i<=m; i++) gel(gig,i) = muliimod_sz(gel(gig,2),gel(gig,i-1),le,lle);
     611       80986 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "polsubcyclo_roots");
     612       80986 :   gel(powz,1) = bab;
     613       80986 :   gel(powz,2) = gig; return powz;
     614             : }
     615             : 
     616             : GEN
     617         406 : galoiscyclo(long n, long v)
     618             : {
     619         406 :   ulong av = avma;
     620             :   GEN grp, G, z, le, L, elts;
     621             :   long val, l, i, j, k;
     622         406 :   GEN zn = znstar(stoi(n));
     623         406 :   long card = itos(gel(zn,1));
     624         406 :   GEN gen = vec_to_vecsmall(lift_shallow(gel(zn,3)));
     625         406 :   GEN ord = vec_to_vecsmall(gel(zn,2));
     626         406 :   GEN T = polcyclo(n,v);
     627         406 :   long d = degpol(T);
     628         406 :   GEN borneabs = powuu(2,d);
     629         406 :   z = polsubcyclo_start(n,card/2,2,2*usqrt(d),borneabs,&val,&l);
     630         406 :   le = gel(z,1); z = gel(z,2);
     631         406 :   L = cgetg(1+card,t_VEC);
     632         406 :   gel(L,1) = z;
     633         854 :   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         406 :   G = abelian_group(ord);
     639         406 :   elts = group_elts(G, card); /*not stack clean*/
     640         406 :   grp = cgetg(9, t_VEC);
     641         406 :   gel(grp,1) = T;
     642         406 :   gel(grp,2) = mkvec3(stoi(l), stoi(val), icopy(le));
     643         406 :   gel(grp,3) = L;
     644         406 :   gel(grp,4) = FpV_invVandermonde(L,  NULL, le);
     645         406 :   gel(grp,5) = gen_1;
     646         406 :   gel(grp,6) = elts;
     647         406 :   gel(grp,7) = gel(G,1);
     648         406 :   gel(grp,8) = gel(G,2);
     649         406 :   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       50205 : _itos(const char *fun, GEN n)
     695             : {
     696       50205 :   if (is_bigint(n))
     697           0 :     pari_err_IMPL(stack_sprintf("conductor f > %ld in %s", LONG_MAX, fun));
     698       50205 :   return itos(n);
     699             : }
     700             : long
     701       50219 : subcyclo_nH(const char *fun, GEN N, GEN *psg)
     702             : {
     703       50219 :   GEN V, Z = NULL, H = *psg;
     704       50219 :   long i, l, n = 0, complex = 1;
     705       50219 :   switch(typ(N))
     706             :   {
     707       49778 :     case t_INT:
     708       49778 :       n = _itos(fun, N);
     709       49778 :       if (n < 1) pari_err_DOMAIN(fun, "degree", "<=", gen_0, N);
     710       49771 :       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       50198 :   if (!H) H = gen_1;
     731       50198 :   switch(typ(H))
     732             :   {
     733       49589 :      case t_INTMOD: case t_INT:
     734       49589 :       V = mkvecsmall( lift_check_modulus(H,n) );
     735       49582 :       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       50163 :   if (!complex) V = vecsmall_append(V, n-1); /*add complex conjugation*/
     757       50163 :   *psg = V; return n;
     758             : }
     759             : 
     760             : GEN
     761       49841 : galoissubcyclo(GEN N, GEN sg, long flag, long v)
     762             : {
     763       49841 :   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       49841 :   if (flag<0 || flag>3) pari_err_FLAG("galoissubcyclo");
     769       49841 :   if (v < 0) v = 0;
     770       49841 :   n = subcyclo_nH("galoissubcyclo", N, &sg);
     771       49792 :   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       49771 :   H = znstar_generate(n, sg);
     777       49771 :   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       49771 :   complex = !F2v_coeff(gel(H,3),n-1);
     786       49771 :   if (DEBUGLEVEL >= 6) err_printf("Subcyclo: complex=%ld\n",complex);
     787       49771 :   if (DEBUGLEVEL >= 1) timer_start(&ti);
     788       49771 :   cnd = znstar_conductor(H);
     789       49771 :   if (DEBUGLEVEL >= 1) timer_printf(&ti, "znstar_conductor");
     790       49771 :   if (flag == 1)  return gc_stoi(ltop, cnd);
     791       49771 :   if (cnd == 1)
     792             :   {
     793          35 :     set_avma(ltop); if (flag == 1) return gen_1;
     794          35 :     return gscycloconductor(deg1pol_shallow(gen_1,gen_m1,v),1,flag);
     795             :   }
     796       49736 :   if (n != cnd)
     797             :   {
     798         175 :     H = znstar_reduce_modulus(H, cnd);
     799         175 :     n = cnd;
     800             :   }
     801       49736 :   card = znstar_order(H);
     802       49736 :   phi_n = eulerphiu(n);
     803       49736 :   if (card == phi_n)
     804             :   {
     805           0 :     set_avma(ltop);
     806           0 :     return gscycloconductor(polcyclo(n,v),n,flag);
     807             :   }
     808       49736 :   O = znstar_cosets(n, phi_n, H);
     809       49736 :   if (DEBUGLEVEL >= 1) timer_printf(&ti, "znstar_cosets");
     810       49736 :   if (DEBUGLEVEL >= 6) err_printf("Subcyclo: orbits=%Ps\n",O);
     811       49736 :   if (DEBUGLEVEL >= 4)
     812           0 :     err_printf("Subcyclo: %ld orbits with %ld elements each\n",phi_n/card,card);
     813       49736 :   av = avma;
     814       49736 :   powz = polsubcyclo_complex_roots(n,!complex,LOWDEFAULTPREC);
     815       49736 :   L = polsubcyclo_orbits(n,H,O,powz,NULL);
     816       49736 :   B = polsubcyclo_complex_bound(av,L,LOWDEFAULTPREC);
     817       49736 :   zl = polsubcyclo_start(n,phi_n/card,card,1,B,&val,&l);
     818       49736 :   powz = polsubcyclo_roots(n,zl);
     819       49736 :   le = gel(zl,1);
     820       49736 :   L = polsubcyclo_orbits(n,H,O,powz,le);
     821       49736 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
     822       49736 :   T = FpV_roots_to_pol(L,le,v);
     823       49736 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "roots_to_pol");
     824       49736 :   T = FpX_center(T,le,shifti(le,-1));
     825       49736 :   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       49729 :   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       33478 : polsubcyclo_g(long n, long d, GEN Z, long v)
     849             : {
     850       33478 :   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       33478 :   if (d==1) return deg1pol_shallow(gen_1,gen_m1,v); /* get rid of n=1,2 */
     855       33478 :   if ((n & 3) == 2) n >>= 1;
     856             :   /* n = 4 or p^a, p odd */
     857       33478 :   o = itos(gel(Z,1));
     858       33478 :   g = itos(gmael3(Z,3,1,2));
     859       33478 :   p = n / ugcd(n,o); /* p^a / gcd(p^a,phi(p^a)) = p*/
     860       33478 :   r = ugcd(d,n); /* = p^(v_p(d)) < n */
     861       33478 :   n = r*p; /* n is now the conductor */
     862       33478 :   o = n-r; /* = phi(n) */
     863       33478 :   if (o == d) return polcyclo(n,v);
     864       31250 :   o /= d;
     865       31250 :   gd = Fl_powu(g%n, d, n);
     866             :   /*FIXME: If degree is small, the computation of B is a waste of time*/
     867       31250 :   powz = polsubcyclo_complex_roots(n,(o&1)==0,LOWDEFAULTPREC);
     868       31250 :   L = polsubcyclo_cyclic(n,d,o,g,gd,powz,NULL);
     869       31250 :   B = polsubcyclo_complex_bound(ltop,L,LOWDEFAULTPREC);
     870       31250 :   zl = polsubcyclo_start(n,d,o,1,B,&val,&l);
     871       31250 :   le = gel(zl,1);
     872       31250 :   powz = polsubcyclo_roots(n,zl);
     873       31250 :   L = polsubcyclo_cyclic(n,d,o,g,gd,powz,le);
     874       31250 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
     875       31250 :   T = FpV_roots_to_pol(L,le,v);
     876       31250 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "roots_to_pol");
     877       31250 :   return gerepileupto(ltop, FpX_center(T,le,shifti(le,-1)));
     878             : }
     879             : 
     880             : GEN
     881       33534 : polsubcyclo(long n, long d, long v)
     882             : {
     883       33534 :   pari_sp ltop = avma;
     884             :   GEN L, Z;
     885       33534 :   if (v<0) v = 0;
     886       33534 :   if (d<=0) pari_err_DOMAIN("polsubcyclo","d","<=",gen_0,stoi(d));
     887       33527 :   if (n<=0) pari_err_DOMAIN("polsubcyclo","n","<=",gen_0,stoi(n));
     888       33520 :   Z = znstar(stoi(n));
     889       33520 :   if (!dvdis(gel(Z,1), d)) { set_avma(ltop); return cgetg(1, t_VEC); }
     890       33513 :   if (lg(gel(Z,2)) == 2)
     891             :   { /* faster but Z must be cyclic */
     892       33478 :     set_avma(ltop);
     893       33478 :     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             : }

Generated by: LCOV version 1.14