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-bordeaux1.fr machine (x86_64 architecture), and agregate them in the final report:

The target is 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.8.0 lcov report (development 16393-29b9383) Lines: 527 563 93.6 %
Date: 2014-04-24 Functions: 38 39 97.4 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 276 348 79.3 %

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

Generated by: LCOV version 1.9