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

Generated by: LCOV version 1.13