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 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 - polarit3.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 19352-1b11b25) Lines: 1368 1593 85.9 %
Date: 2016-08-25 06:11:27 Functions: 128 144 88.9 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000-2005  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             : /***********************************************************************/
      15             : /**                                                                   **/
      16             : /**               ARITHMETIC OPERATIONS ON POLYNOMIALS                **/
      17             : /**                         (third part)                              **/
      18             : /**                                                                   **/
      19             : /***********************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : /************************************************************************
      24             :  **                                                                    **
      25             :  **                      Ring membership                               **
      26             :  **                                                                    **
      27             :  ************************************************************************/
      28             : struct charact {
      29             :   GEN q;
      30             :   int isprime;
      31             : };
      32             : static void
      33        3290 : char_update_prime(struct charact *S, GEN p)
      34             : {
      35        3290 :   if (!S->isprime) { S->isprime = 1; S->q = p; }
      36        3290 :   if (!equalii(p, S->q)) pari_err_MODULUS("characteristic", S->q, p);
      37        3283 : }
      38             : static void
      39        3423 : char_update_int(struct charact *S, GEN n)
      40             : {
      41        3423 :   if (S->isprime)
      42             :   {
      43        3423 :     if (dvdii(n, S->q)) return;
      44           0 :     pari_err_MODULUS("characteristic", S->q, n);
      45             :   }
      46        3423 :   S->q = gcdii(S->q, n);
      47             : }
      48             : static void
      49      559965 : charact(struct charact *S, GEN x)
      50             : {
      51      559965 :   const long tx = typ(x);
      52             :   long i, l;
      53      559965 :   switch(tx)
      54             :   {
      55        3157 :     case t_INTMOD:char_update_int(S, gel(x,1)); break;
      56        3227 :     case t_FFELT: char_update_prime(S, gel(x,4)); break;
      57             :     case t_COMPLEX: case t_QUAD:
      58             :     case t_POLMOD: case t_POL: case t_SER: case t_RFRAC:
      59             :     case t_VEC: case t_COL: case t_MAT:
      60        6454 :       l = lg(x);
      61        6454 :       for (i=lontyp[tx]; i < l; i++) charact(S,gel(x,i));
      62        6447 :       break;
      63             :     case t_LIST:
      64           7 :       x = list_data(x);
      65           7 :       if (x) charact(S, x);
      66           7 :       break;
      67             :   }
      68      559951 : }
      69             : static void
      70       30296 : charact_res(struct charact *S, GEN x)
      71             : {
      72       30296 :   const long tx = typ(x);
      73             :   long i, l;
      74       30296 :   switch(tx)
      75             :   {
      76         266 :     case t_INTMOD:char_update_int(S, gel(x,1)); break;
      77           0 :     case t_FFELT: char_update_prime(S, gel(x,4)); break;
      78          63 :     case t_PADIC: char_update_prime(S, gel(x,2)); break;
      79             :     case t_COMPLEX: case t_QUAD:
      80             :     case t_POLMOD: case t_POL: case t_SER: case t_RFRAC:
      81             :     case t_VEC: case t_COL: case t_MAT:
      82        9317 :       l = lg(x);
      83        9317 :       for (i=lontyp[tx]; i < l; i++) charact_res(S,gel(x,i));
      84        9317 :       break;
      85             :     case t_LIST:
      86           0 :       x = list_data(x);
      87           0 :       if (x) charact_res(S, x);
      88           0 :       break;
      89             :   }
      90       30296 : }
      91             : GEN
      92        5061 : characteristic(GEN x)
      93             : {
      94             :   struct charact S;
      95        5061 :   S.q = gen_0; S.isprime = 0;
      96        5061 :   charact(&S, x); return S.q;
      97             : }
      98             : GEN
      99        2324 : residual_characteristic(GEN x)
     100             : {
     101             :   struct charact S;
     102        2324 :   S.q = gen_0; S.isprime = 0;
     103        2324 :   charact_res(&S, x); return S.q;
     104             : }
     105             : 
     106             : int
     107    16536253 : Rg_is_Fp(GEN x, GEN *pp)
     108             : {
     109             :   GEN mod;
     110    16536253 :   switch(typ(x))
     111             :   {
     112             :   case t_INTMOD:
     113     3849629 :     mod = gel(x,1);
     114     3849629 :     if (!*pp) *pp = mod;
     115     3797353 :     else if (mod != *pp && !equalii(mod, *pp))
     116             :     {
     117           0 :       if (DEBUGLEVEL) pari_warn(warner,"different moduli in Rg_is_Fp");
     118           0 :       return 0;
     119             :     }
     120     3849629 :     return 1;
     121             :   case t_INT:
     122    10368169 :     return 1;
     123     2318455 :   default: return 0;
     124             :   }
     125             : }
     126             : 
     127             : int
     128     2349608 : RgX_is_FpX(GEN x, GEN *pp)
     129             : {
     130     2349608 :   long i, lx = lg(x);
     131     8975794 :   for (i=2; i<lx; i++)
     132     8334430 :     if (!Rg_is_Fp(gel(x, i), pp))
     133     1708244 :       return 0;
     134      641364 :   return 1;
     135             : }
     136             : 
     137             : int
     138     1729929 : RgV_is_FpV(GEN x, GEN *pp)
     139             : {
     140     1729929 :   long i, lx = lg(x);
     141     9299386 :   for (i=1; i<lx; i++)
     142     8179661 :     if (!Rg_is_Fp(gel(x,i), pp)) return 0;
     143     1119725 :   return 1;
     144             : }
     145             : 
     146             : int
     147      737691 : RgM_is_FpM(GEN x, GEN *pp)
     148             : {
     149      737691 :   long i, lx = lg(x);
     150     1854980 :   for (i=1; i<lx; i++)
     151     1727416 :     if (!RgV_is_FpV(gel(x, i), pp)) return 0;
     152      127564 :   return 1;
     153             : }
     154             : 
     155             : int
     156       95634 : Rg_is_FpXQ(GEN x, GEN *pT, GEN *pp)
     157             : {
     158             :   GEN pol, mod, p;
     159       95634 :   switch(typ(x))
     160             :   {
     161             :   case t_INTMOD:
     162       22092 :     return Rg_is_Fp(x, pp);
     163             :   case t_INT:
     164       21630 :     return 1;
     165             :   case t_POL:
     166        9240 :     return RgX_is_FpX(x, pp);
     167             :   case t_FFELT:
     168       10010 :     mod = FF_1(x); p = FF_p_i(x);
     169       10010 :     if (!*pp) *pp = p;
     170       10010 :     if (!*pT) *pT = mod;
     171       10010 :     if ((p != *pp && !equalii(p, *pp)) || (mod != *pT && !gequal(mod, *pT)))
     172             :     {
     173           0 :       if (DEBUGLEVEL) pari_warn(warner,"different moduli in Rg_is_FpXQ");
     174           0 :       return 0;
     175             :     }
     176       10010 :     return 1;
     177             :   case t_POLMOD:
     178        5397 :     mod = gel(x,1); pol = gel(x, 2);
     179        5397 :     if (!RgX_is_FpX(mod, pp)) return 0;
     180        5397 :     if (typ(pol)==t_POL)
     181             :     {
     182        5327 :       if (!RgX_is_FpX(pol, pp)) return 0;
     183             :     }
     184          70 :     else if (!Rg_is_Fp(pol, pp)) return 0;
     185        5376 :     if (!*pT) *pT = mod;
     186        4718 :     else if (mod != *pT && !gequal(mod, *pT))
     187             :     {
     188           0 :       if (DEBUGLEVEL) pari_warn(warner,"different moduli in Rg_is_FpXQ");
     189           0 :       return 0;
     190             :     }
     191        5376 :     return 1;
     192             : 
     193       27265 :   default: return 0;
     194             :   }
     195             : }
     196             : 
     197             : int
     198       36232 : RgX_is_FpXQX(GEN x, GEN *pT, GEN *pp)
     199             : {
     200       36232 :   long i, lx = lg(x);
     201      102550 :   for (i = 2; i < lx; i++)
     202       95081 :     if (!Rg_is_FpXQ(gel(x,i), pT, pp)) return 0;
     203        7469 :   return 1;
     204             : }
     205             : 
     206             : /************************************************************************
     207             :  **                                                                    **
     208             :  **                      Ring conversion                               **
     209             :  **                                                                    **
     210             :  ************************************************************************/
     211             : 
     212             : /* p > 0 a t_INT, return lift(x * Mod(1,p)).
     213             :  * If x is an INTMOD, assume modulus is a multiple of p. */
     214             : GEN
     215    26640211 : Rg_to_Fp(GEN x, GEN p)
     216             : {
     217    26640211 :   if (lgefint(p) == 3) return utoi(Rg_to_Fl(x, uel(p,2)));
     218      386702 :   switch(typ(x))
     219             :   {
     220       25731 :     case t_INT: return modii(x, p);
     221             :     case t_FRAC: {
     222          60 :       pari_sp av = avma;
     223          60 :       GEN z = modii(gel(x,1), p);
     224          60 :       if (z == gen_0) return gen_0;
     225          60 :       return gerepileuptoint(av, remii(mulii(z, Fp_inv(gel(x,2), p)), p));
     226             :     }
     227           0 :     case t_PADIC: return padic_to_Fp(x, p);
     228             :     case t_INTMOD: {
     229      360911 :       GEN q = gel(x,1), a = gel(x,2);
     230      360911 :       if (equalii(q, p)) return icopy(a);
     231           7 :       if (!dvdii(q,p)) pari_err_MODULUS("Rg_to_Fp", q, p);
     232           0 :       return remii(a, p);
     233             :     }
     234           0 :     default: pari_err_TYPE("Rg_to_Fp",x);
     235           0 :       return NULL; /* not reached */
     236             :   }
     237             : }
     238             : /* If x is a POLMOD, assume modulus is a multiple of T. */
     239             : GEN
     240       46264 : Rg_to_FpXQ(GEN x, GEN T, GEN p)
     241             : {
     242       46264 :   long ta, tx = typ(x), v = get_FpX_var(T);
     243             :   GEN a, b;
     244       46264 :   if (is_const_t(tx))
     245             :   {
     246       41777 :     if (tx == t_FFELT)
     247             :     {
     248       21812 :       GEN z = FF_to_FpXQ(x);
     249       21812 :       setvarn(z, v);
     250       21812 :       return z;
     251             :     }
     252       19965 :     return scalar_ZX(Rg_to_Fp(x, p), v);
     253             :   }
     254        4487 :   switch(tx)
     255             :   {
     256             :     case t_POLMOD:
     257        3724 :       b = gel(x,1);
     258        3724 :       a = gel(x,2); ta = typ(a);
     259        3724 :       if (is_const_t(ta)) return scalar_ZX(Rg_to_Fp(a, p), v);
     260        3682 :       b = RgX_to_FpX(b, p); if (varn(b) != v) break;
     261        3682 :       a = RgX_to_FpX(a, p); if (ZX_equal(b,get_FpX_mod(T))) return a;
     262           0 :       if (signe(FpX_rem(b,T,p))==0) return FpX_rem(a, T, p);
     263           0 :       break;
     264             :     case t_POL:
     265         763 :       if (varn(x) != v) break;
     266         763 :       return FpX_rem(RgX_to_FpX(x,p), T, p);
     267             :     case t_RFRAC:
     268           0 :       a = Rg_to_FpXQ(gel(x,1), T,p);
     269           0 :       b = Rg_to_FpXQ(gel(x,2), T,p);
     270           0 :       return FpXQ_div(a,b, T,p);
     271             :   }
     272           0 :   pari_err_TYPE("Rg_to_FpXQ",x);
     273           0 :   return NULL; /* not reached */
     274             : }
     275             : GEN
     276      140588 : RgX_to_FpX(GEN x, GEN p)
     277             : {
     278             :   long i, l;
     279      140588 :   GEN z = cgetg_copy(x, &l); z[1] = x[1];
     280      140588 :   for (i = 2; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
     281      140588 :   return FpX_renormalize(z, l);
     282             : }
     283             : 
     284             : GEN
     285        1099 : RgV_to_FpV(GEN x, GEN p)
     286             : {
     287        1099 :   long i, l = lg(x);
     288        1099 :   GEN z = cgetg(l, t_VEC);
     289        1099 :   for (i = 1; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
     290        1099 :   return z;
     291             : }
     292             : 
     293             : GEN
     294      607977 : RgC_to_FpC(GEN x, GEN p)
     295             : {
     296      607977 :   long i, l = lg(x);
     297      607977 :   GEN z = cgetg(l, t_COL);
     298      607977 :   for (i = 1; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
     299      607977 :   return z;
     300             : }
     301             : 
     302             : GEN
     303       44663 : RgM_to_FpM(GEN x, GEN p)
     304             : {
     305       44663 :   long i, l = lg(x);
     306       44663 :   GEN z = cgetg(l, t_MAT);
     307       44663 :   for (i = 1; i < l; i++) gel(z,i) = RgC_to_FpC(gel(x,i), p);
     308       44663 :   return z;
     309             : }
     310             : GEN
     311       12615 : RgV_to_Flv(GEN x, ulong p)
     312             : {
     313       12615 :   long l = lg(x), i;
     314       12615 :   GEN a = cgetg(l, t_VECSMALL);
     315       12615 :   for (i=1; i<l; i++) a[i] = Rg_to_Fl(gel(x,i), p);
     316       12615 :   return a;
     317             : }
     318             : GEN
     319        1768 : RgM_to_Flm(GEN x, ulong p)
     320             : {
     321             :   long l, i;
     322        1768 :   GEN a = cgetg_copy(x, &l);
     323        1768 :   for (i=1; i<l; i++) gel(a,i) = RgV_to_Flv(gel(x,i), p);
     324        1768 :   return a;
     325             : }
     326             : 
     327             : GEN
     328         336 : RgX_to_FpXQX(GEN x, GEN T, GEN p)
     329             : {
     330         336 :   long i, l = lg(x);
     331         336 :   GEN z = cgetg(l, t_POL); z[1] = x[1];
     332         336 :   for (i = 2; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T,p);
     333         336 :   return FpXQX_renormalize(z, l);
     334             : }
     335             : GEN
     336         665 : RgX_to_FqX(GEN x, GEN T, GEN p)
     337             : {
     338         665 :   long i, l = lg(x);
     339         665 :   GEN z = cgetg(l, t_POL); z[1] = x[1];
     340         665 :   if (T)
     341       11263 :     for (i = 2; i < l; i++)
     342       10640 :       gel(z,i) = simplify_shallow(Rg_to_FpXQ(gel(x,i), T, p));
     343             :   else
     344        1554 :     for (i = 2; i < l; i++)
     345        1512 :       gel(z,i) = Rg_to_Fp(gel(x,i), p);
     346         665 :   return FpXQX_renormalize(z, l);
     347             : }
     348             : 
     349             : /* lg(V) > 1 */
     350             : GEN
     351      848638 : FpXV_FpC_mul(GEN V, GEN W, GEN p)
     352             : {
     353      848638 :   pari_sp av = avma;
     354      848638 :   long i, l = lg(V);
     355      848638 :   GEN z = ZX_Z_mul(gel(V,1),gel(W,1));
     356     4173540 :   for(i=2; i<l; i++)
     357             :   {
     358     3324902 :     z = ZX_add(z, ZX_Z_mul(gel(V,i),gel(W,i)));
     359     3324902 :     if ((i & 7) == 0) z = gerepileupto(av, z);
     360             :   }
     361      848638 :   return gerepileupto(av, FpX_red(z,p));
     362             : }
     363             : 
     364             : GEN
     365        2275 : FqX_Fq_add(GEN y, GEN x, GEN T, GEN p)
     366             : {
     367        2275 :   long i, lz = lg(y);
     368             :   GEN z;
     369        2275 :   if (!T) return FpX_Fp_add(y, x, p);
     370        2275 :   if (lz == 2) return scalarpol(x, varn(y));
     371        2275 :   z = cgetg(lz,t_POL); z[1] = y[1];
     372        2275 :   gel(z,2) = Fq_add(gel(y,2),x, T, p);
     373        2275 :   if (lz == 3) z = FpXX_renormalize(z,lz);
     374             :   else
     375        1239 :     for(i=3;i<lz;i++) gel(z,i) = gcopy(gel(y,i));
     376        2275 :   return z;
     377             : }
     378             : 
     379             : GEN
     380       66600 : FqX_Fq_mul_to_monic(GEN P, GEN U, GEN T, GEN p)
     381             : {
     382             :   long i, lP;
     383       66600 :   GEN res = cgetg_copy(P, &lP); res[1] = P[1];
     384       66600 :   for(i=2; i<lP-1; i++) gel(res,i) = Fq_mul(U,gel(P,i), T,p);
     385       66600 :   gel(res,lP-1) = gen_1; return res;
     386             : }
     387             : 
     388             : GEN
     389        7659 : FpXQX_normalize(GEN z, GEN T, GEN p)
     390             : {
     391             :   GEN lc;
     392        7659 :   if (lg(z) == 2) return z;
     393        7659 :   lc = leading_coeff(z);
     394        7659 :   if (typ(lc) == t_POL)
     395             :   {
     396        3916 :     if (lg(lc) > 3) /* non-constant */
     397        3901 :       return FqX_Fq_mul_to_monic(z, Fq_inv(lc,T,p), T,p);
     398             :     /* constant */
     399          15 :     lc = gel(lc,2);
     400          15 :     z = shallowcopy(z);
     401          15 :     gel(z, lg(z)-1) = lc;
     402             :   }
     403             :   /* lc a t_INT */
     404        3758 :   if (equali1(lc)) return z;
     405         609 :   return FqX_Fq_mul_to_monic(z, Fp_inv(lc,p), T,p);
     406             : }
     407             : 
     408             : GEN
     409      119658 : FqX_eval(GEN x, GEN y, GEN T, GEN p)
     410             : {
     411             :   pari_sp av;
     412             :   GEN p1, r;
     413      119658 :   long j, i=lg(x)-1;
     414      119658 :   if (i<=2)
     415       24136 :     return (i==2)? Fq_red(gel(x,2), T, p): gen_0;
     416       95522 :   av=avma; p1=gel(x,i);
     417             :   /* specific attention to sparse polynomials (see poleval)*/
     418             :   /*You've guessed it! It's a copy-paste(tm)*/
     419      283087 :   for (i--; i>=2; i=j-1)
     420             :   {
     421      187880 :     for (j=i; !signe(gel(x,j)); j--)
     422         315 :       if (j==2)
     423             :       {
     424         182 :         if (i!=j) y = Fq_pow(y,utoipos(i-j+1), T, p);
     425         182 :         return gerepileupto(av, Fq_mul(p1,y, T, p));
     426             :       }
     427      187565 :     r = (i==j)? y: Fq_pow(y, utoipos(i-j+1), T, p);
     428      187565 :     p1 = Fq_add(Fq_mul(p1,r,T,p), gel(x,j), T, p);
     429             :   }
     430       95340 :   return gerepileupto(av, p1);
     431             : }
     432             : 
     433             : GEN
     434       30380 : FqXY_evalx(GEN Q, GEN x, GEN T, GEN p)
     435             : {
     436       30380 :   long i, lb = lg(Q);
     437             :   GEN z;
     438       30380 :   if (!T) return FpXY_evalx(Q, x, p);
     439       20720 :   z = cgetg(lb, t_POL); z[1] = Q[1];
     440      115969 :   for (i=2; i<lb; i++)
     441             :   {
     442       95249 :     GEN q = gel(Q,i);
     443       95249 :     gel(z,i) = typ(q) == t_INT? modii(q,p): FqX_eval(q, x, T, p);
     444             :   }
     445       20720 :   return FpXQX_renormalize(z, lb);
     446             : }
     447             : 
     448             : /* Q an FpXY, evaluate at (X,Y) = (x,y) */
     449             : GEN
     450       12733 : FqXY_eval(GEN Q, GEN y, GEN x, GEN T, GEN p)
     451             : {
     452       12733 :   pari_sp av = avma;
     453       12733 :   if (!T) return FpXY_eval(Q, y, x, p);
     454         336 :   return gerepileupto(av, FqX_eval(FqXY_evalx(Q, x, T, p), y, T, p));
     455             : }
     456             : 
     457             : /* a X^d */
     458             : GEN
     459     2416798 : monomial(GEN a, long d, long v)
     460             : {
     461     2416798 :   long i, lP = d+3;
     462             :   GEN P;
     463     2416798 :   if (d < 0) {
     464          56 :     if (isrationalzero(a)) return pol_0(v);
     465          56 :     P = cgetg(3, t_RFRAC);
     466          56 :     gel(P,1) = a;
     467          56 :     gel(P,2) = monomial(gen_1, -d, v);
     468             :   } else {
     469     2416742 :     P = cgetg(lP, t_POL);
     470     2416742 :     if (gequal0(a))
     471             :     {
     472         266 :       if (isexactzero(a)) return scalarpol_shallow(a,v);
     473           0 :       P[1] = evalsigne(0) | evalvarn(v);
     474             :     }
     475             :     else
     476     2416476 :       P[1] = evalsigne(1) | evalvarn(v);
     477     2416476 :     lP--; gel(P,lP) = a;
     478     2416476 :     for (i=2; i<lP; i++) gel(P,i) = gen_0;
     479             :   }
     480     2416532 :   return P;
     481             : }
     482             : GEN
     483     7307226 : monomialcopy(GEN a, long d, long v)
     484             : {
     485     7307226 :   long i, lP = d+3;
     486             :   GEN P;
     487     7307226 :   if (d < 0) {
     488           7 :     if (isrationalzero(a)) return pol_0(v);
     489           7 :     P = cgetg(3, t_RFRAC);
     490           7 :     gel(P,1) = gcopy(a);
     491           7 :     gel(P,2) = monomial(gen_1, -d, v);
     492             :   } else {
     493     7307219 :     P = cgetg(lP, t_POL);
     494     7307219 :     if (gequal0(a))
     495             :     {
     496           7 :       if (isexactzero(a)) return scalarpol(a,v);
     497           0 :       P[1] = evalsigne(0) | evalvarn(v);
     498             :     }
     499             :     else
     500     7307212 :       P[1] = evalsigne(1) | evalvarn(v);
     501     7307212 :     lP--; gel(P,lP) = gcopy(a);
     502     7307212 :     for (i=2; i<lP; i++) gel(P,i) = gen_0;
     503             :   }
     504     7307219 :   return P;
     505             : }
     506             : GEN
     507       19390 : pol_x_powers(long N, long v)
     508             : {
     509       19390 :   GEN L = cgetg(N+1,t_VEC);
     510             :   long i;
     511       19390 :   for (i=1; i<=N; i++) gel(L,i) = monomial(gen_1, i-1, v);
     512       19390 :   return L;
     513             : }
     514             : 
     515             : GEN
     516           0 : FqXQ_powers(GEN x, long l, GEN S, GEN T, GEN p)
     517             : {
     518           0 :   return T ? FpXQXQ_powers(x, l, S, T, p): FpXQ_powers(x, l, S, p);
     519             : }
     520             : 
     521             : GEN
     522           0 : FqXQ_matrix_pow(GEN y, long n, long m, GEN S, GEN T, GEN p)
     523             : {
     524           0 :   return T ? FpXQXQ_matrix_pow(y, n, m, S, T, p): FpXQ_matrix_pow(y, n, m, S, p);
     525             : }
     526             : 
     527             : /*******************************************************************/
     528             : /*                                                                 */
     529             : /*                             Fq                                  */
     530             : /*                                                                 */
     531             : /*******************************************************************/
     532             : 
     533             : GEN
     534     6966160 : Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p)
     535             : {
     536             :   (void)T;
     537     6966160 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     538             :   {
     539     2695468 :     case 0: return Fp_add(x,y,p);
     540      204218 :     case 1: return FpX_Fp_add(x,y,p);
     541      337828 :     case 2: return FpX_Fp_add(y,x,p);
     542     3728646 :     case 3: return FpX_add(x,y,p);
     543             :   }
     544           0 :   return NULL;
     545             : }
     546             : 
     547             : GEN
     548     4733499 : Fq_sub(GEN x, GEN y, GEN T/*unused*/, GEN p)
     549             : {
     550             :   (void)T;
     551     4733499 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     552             :   {
     553      160812 :     case 0: return Fp_sub(x,y,p);
     554        1231 :     case 1: return FpX_Fp_sub(x,y,p);
     555        8627 :     case 2: return Fp_FpX_sub(x,y,p);
     556     4562829 :     case 3: return FpX_sub(x,y,p);
     557             :   }
     558           0 :   return NULL;
     559             : }
     560             : 
     561             : GEN
     562      529167 : Fq_neg(GEN x, GEN T/*unused*/, GEN p)
     563             : {
     564             :   (void)T;
     565      529167 :   return (typ(x)==t_POL)? FpX_neg(x,p): Fp_neg(x,p);
     566             : }
     567             : 
     568             : GEN
     569       11505 : Fq_halve(GEN x, GEN T/*unused*/, GEN p)
     570             : {
     571             :   (void)T;
     572       11505 :   return (typ(x)==t_POL)? FpX_halve(x,p): Fp_halve(x,p);
     573             : }
     574             : 
     575             : /* If T==NULL do not reduce*/
     576             : GEN
     577    42811568 : Fq_mul(GEN x, GEN y, GEN T, GEN p)
     578             : {
     579    42811568 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     580             :   {
     581     2702880 :     case 0: return Fp_mul(x,y,p);
     582      197035 :     case 1: return FpX_Fp_mul(x,y,p);
     583      122937 :     case 2: return FpX_Fp_mul(y,x,p);
     584    39788716 :     case 3: if (T) return FpXQ_mul(x,y,T,p);
     585     3060049 :             else return FpX_mul(x,y,p);
     586             :   }
     587           0 :   return NULL;
     588             : }
     589             : 
     590             : /* If T==NULL do not reduce*/
     591             : GEN
     592      854223 : Fq_mulu(GEN x, ulong y, /*unused*/GEN T, GEN p)
     593             : {
     594             :   (void) T;
     595      854223 :   return typ(x)==t_POL ? FpX_Fp_mul(x,utoi(y),p): Fp_mulu(x, y, p);
     596             : }
     597             : 
     598             : /* y t_INT */
     599             : GEN
     600       51749 : Fq_Fp_mul(GEN x, GEN y, GEN T/*unused*/, GEN p)
     601             : {
     602             :   (void)T;
     603      103498 :   return (typ(x) == t_POL)? FpX_Fp_mul(x,y,p)
     604       51749 :                           : Fp_mul(x,y,p);
     605             : }
     606             : /* If T==NULL do not reduce*/
     607             : GEN
     608      264871 : Fq_sqr(GEN x, GEN T, GEN p)
     609             : {
     610      264871 :   if (typ(x) == t_POL)
     611             :   {
     612       12416 :     if (T) return FpXQ_sqr(x,T,p);
     613           0 :     else return FpX_sqr(x,p);
     614             :   }
     615             :   else
     616      252455 :     return Fp_sqr(x,p);
     617             : }
     618             : 
     619             : GEN
     620           0 : Fq_neg_inv(GEN x, GEN T, GEN p)
     621             : {
     622           0 :   if (typ(x) == t_INT) return Fp_inv(Fp_neg(x,p),p);
     623           0 :   return FpXQ_inv(FpX_neg(x,p),T,p);
     624             : }
     625             : 
     626             : GEN
     627           0 : Fq_invsafe(GEN x, GEN pol, GEN p)
     628             : {
     629           0 :   if (typ(x) == t_INT) return Fp_invsafe(x,p);
     630           0 :   return FpXQ_invsafe(x,pol,p);
     631             : }
     632             : 
     633             : GEN
     634       22793 : Fq_inv(GEN x, GEN pol, GEN p)
     635             : {
     636       22793 :   if (typ(x) == t_INT) return Fp_inv(x,p);
     637       18110 :   return FpXQ_inv(x,pol,p);
     638             : }
     639             : 
     640             : GEN
     641      481607 : Fq_div(GEN x, GEN y, GEN pol, GEN p)
     642             : {
     643      481607 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     644             :   {
     645      453670 :     case 0: return Fp_div(x,y,p);
     646       23002 :     case 1: return FpX_Fp_mul(x,Fp_inv(y,p),p);
     647         147 :     case 2: return FpX_Fp_mul(FpXQ_inv(y,pol,p),x,p);
     648        4788 :     case 3: return FpXQ_div(x,y,pol,p);
     649             :   }
     650           0 :   return NULL;
     651             : }
     652             : 
     653             : GEN
     654       11718 : Fq_pow(GEN x, GEN n, GEN pol, GEN p)
     655             : {
     656       11718 :   if (typ(x) == t_INT) return Fp_pow(x,n,p);
     657        8806 :   return FpXQ_pow(x,n,pol,p);
     658             : }
     659             : 
     660             : GEN
     661       10647 : Fq_powu(GEN x, ulong n, GEN pol, GEN p)
     662             : {
     663       10647 :   if (typ(x) == t_INT) return Fp_powu(x,n,p);
     664         322 :   return FpXQ_powu(x,n,pol,p);
     665             : }
     666             : 
     667             : GEN
     668      698772 : Fq_sqrt(GEN x, GEN T, GEN p)
     669             : {
     670      698772 :   if (typ(x) == t_INT)
     671             :   {
     672      692370 :     if (!T || odd(get_FpX_degree(T))) return Fp_sqrt(x,p);
     673           0 :     x = scalarpol_shallow(x, get_FpX_var(T));
     674             :   }
     675        6402 :   return FpXQ_sqrt(x,T,p);
     676             : }
     677             : GEN
     678       60452 : Fq_sqrtn(GEN x, GEN n, GEN T, GEN p, GEN *zeta)
     679             : {
     680       60452 :   if (typ(x) == t_INT)
     681             :   {
     682             :     long d;
     683       59647 :     if (!T) return Fp_sqrtn(x,n,p,zeta);
     684         182 :     d = get_FpX_degree(T);
     685         182 :     if (ugcd(umodiu(n,d),d) == 1)
     686             :     {
     687          14 :       if (!zeta)
     688           7 :         return Fp_sqrtn(x,n,p,NULL);
     689             :       else
     690             :       {
     691             :         /* gcd(n,p-1)=gcd(n,p^d-1) <=> same number of solutions if Fp and F_{p^d} */
     692           7 :         if (equalii(gcdii(subiu(p,1),n), gcdii(subiu(Fp_powu(p,d,n), 1), n)))
     693           7 :           return Fp_sqrtn(x,n,p,zeta);
     694             :       }
     695             :     }
     696         168 :     x = scalarpol(x, get_FpX_var(T)); /* left on stack */
     697             :   }
     698         973 :   return FpXQ_sqrtn(x,n,T,p,zeta);
     699             : }
     700             : 
     701             : struct _Fq_field
     702             : {
     703             :   GEN T, p;
     704             : };
     705             : 
     706             : static GEN
     707       49696 : _Fq_red(void *E, GEN x)
     708       49696 : { struct _Fq_field *s = (struct _Fq_field *)E;
     709       49696 :   return Fq_red(x, s->T, s->p);
     710             : }
     711             : 
     712             : static GEN
     713      613027 : _Fq_add(void *E, GEN x, GEN y)
     714             : {
     715             :   (void) E;
     716      613027 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     717             :   {
     718           0 :     case 0: return addii(x,y);
     719           0 :     case 1: return ZX_Z_add(x,y);
     720          90 :     case 2: return ZX_Z_add(y,x);
     721      612937 :     default: return ZX_add(x,y);
     722             :   }
     723             : }
     724             : 
     725             : static GEN
     726        1287 : _Fq_neg(void *E, GEN x) { (void) E; return typ(x)==t_POL?ZX_neg(x):negi(x); }
     727             : 
     728             : static GEN
     729      637241 : _Fq_mul(void *E, GEN x, GEN y)
     730             : {
     731             :   (void) E;
     732      637241 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     733             :   {
     734           0 :     case 0: return mulii(x,y);
     735          42 :     case 1: return ZX_Z_mul(x,y);
     736          28 :     case 2: return ZX_Z_mul(y,x);
     737      637171 :     default: return ZX_mul(x,y);
     738             :   }
     739             : }
     740             : 
     741             : static GEN
     742        1210 : _Fq_inv(void *E, GEN x)
     743        1210 : { struct _Fq_field *s = (struct _Fq_field *)E;
     744        1210 :   return Fq_inv(x,s->T,s->p);
     745             : }
     746             : 
     747             : static int
     748       25182 : _Fq_equal0(GEN x) { return signe(x)==0; }
     749             : 
     750             : static GEN
     751        1534 : _Fq_s(void *E, long x) { (void) E; return stoi(x); }
     752             : 
     753             : static const struct bb_field Fq_field={_Fq_red,_Fq_add,_Fq_mul,_Fq_neg,
     754             :                                        _Fq_inv,_Fq_equal0,_Fq_s};
     755             : 
     756         427 : const struct bb_field *get_Fq_field(void **E, GEN T, GEN p)
     757             : {
     758         427 :   GEN z = new_chunk(sizeof(struct _Fq_field));
     759         427 :   struct _Fq_field *e = (struct _Fq_field *) z;
     760         427 :   e->T = T; e->p  = p; *E = (void*)e;
     761         427 :   return &Fq_field;
     762             : }
     763             : 
     764             : /*******************************************************************/
     765             : /*                                                                 */
     766             : /*                             Fq[X]                               */
     767             : /*                                                                 */
     768             : /*******************************************************************/
     769             : /* P(X + c) */
     770             : GEN
     771           0 : FpX_translate(GEN P, GEN c, GEN p)
     772             : {
     773           0 :   pari_sp av = avma;
     774             :   GEN Q, *R;
     775             :   long i, k, n;
     776             : 
     777           0 :   if (!signe(P) || !signe(c)) return ZX_copy(P);
     778           0 :   Q = leafcopy(P);
     779           0 :   R = (GEN*)(Q+2); n = degpol(P);
     780           0 :   for (i=1; i<=n; i++)
     781             :   {
     782           0 :     for (k=n-i; k<n; k++)
     783           0 :       R[k] = Fp_add(R[k], Fp_mul(c, R[k+1], p), p);
     784             : 
     785           0 :     if (gc_needed(av,2))
     786             :     {
     787           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FpX_translate, i = %ld/%ld", i,n);
     788           0 :       Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
     789             :     }
     790             :   }
     791           0 :   return gerepilecopy(av, FpX_renormalize(Q, lg(Q)));
     792             : }
     793             : /* P(X + c), c an Fq */
     794             : GEN
     795       40348 : FqX_translate(GEN P, GEN c, GEN T, GEN p)
     796             : {
     797       40348 :   pari_sp av = avma;
     798             :   GEN Q, *R;
     799             :   long i, k, n;
     800             : 
     801             :   /* signe works for t_(INT|POL) */
     802       40348 :   if (!signe(P) || !signe(c)) return RgX_copy(P);
     803       40348 :   Q = leafcopy(P);
     804       40348 :   R = (GEN*)(Q+2); n = degpol(P);
     805      182987 :   for (i=1; i<=n; i++)
     806             :   {
     807      787136 :     for (k=n-i; k<n; k++)
     808      644497 :       R[k] = Fq_add(R[k], Fq_mul(c, R[k+1], T, p), T, p);
     809             : 
     810      142639 :     if (gc_needed(av,2))
     811             :     {
     812           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FqX_translate, i = %ld/%ld", i,n);
     813           0 :       Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
     814             :     }
     815             :   }
     816       40348 :   return gerepilecopy(av, FpXQX_renormalize(Q, lg(Q)));
     817             : }
     818             : 
     819             : GEN
     820         588 : FqV_roots_to_pol(GEN V, GEN T, GEN p, long v)
     821             : {
     822         588 :   pari_sp ltop = avma;
     823             :   long k;
     824             :   GEN W;
     825         588 :   if (lgefint(p) == 3)
     826             :   {
     827         549 :     ulong pp = p[2];
     828         549 :     GEN Tl = ZX_to_Flx(T, pp);
     829         549 :     GEN Vl = FqV_to_FlxV(V, T, p);
     830         549 :     Tl = FlxqV_roots_to_pol(Vl, Tl, pp, v);
     831         549 :     return gerepileupto(ltop, FlxX_to_ZXX(Tl));
     832             :   }
     833          39 :   W = cgetg(lg(V),t_VEC);
     834         255 :   for(k=1; k < lg(V); k++)
     835         216 :     gel(W,k) = deg1pol_shallow(gen_1,Fq_neg(gel(V,k),T,p),v);
     836          39 :   return gerepileupto(ltop, FpXQXV_prod(W, T, p));
     837             : }
     838             : 
     839             : GEN
     840       10605 : FqV_red(GEN z, GEN T, GEN p)
     841             : {
     842       10605 :   long i, l = lg(z);
     843       10605 :   GEN res = cgetg(l, typ(z));
     844       10605 :   for(i=1;i<l;i++) gel(res,i) = Fq_red(gel(z,i),T,p);
     845       10605 :   return res;
     846             : }
     847             : 
     848             : GEN
     849           0 : FqC_add(GEN x, GEN y, GEN T, GEN p)
     850             : {
     851           0 :   long i, lx = lg(x);
     852             :   GEN z;
     853           0 :   if (!T) return FpC_add(x, y, p);
     854           0 :   z = cgetg(lx, t_COL);
     855           0 :   for (i = 1; i < lx; i++) gel(z, i) = Fq_add(gel(x, i), gel(y, i), T, p);
     856           0 :   return z;
     857             : }
     858             : 
     859             : GEN
     860           0 : FqC_sub(GEN x, GEN y, GEN T, GEN p)
     861             : {
     862           0 :   long i, lx = lg(x);
     863             :   GEN z;
     864           0 :   if (!T) return FpC_sub(x, y, p);
     865           0 :   z = cgetg(lx, t_COL);
     866           0 :   for (i = 1; i < lx; i++) gel(z, i) = Fq_sub(gel(x, i), gel(y, i), T, p);
     867           0 :   return z;
     868             : }
     869             : 
     870             : GEN
     871           0 : FqC_Fq_mul(GEN x, GEN y, GEN T, GEN p)
     872             : {
     873           0 :   long i, l = lg(x);
     874             :   GEN z;
     875           0 :   if (!T) return FpC_Fp_mul(x, y, p);
     876           0 :   z = cgetg(l, t_COL);
     877           0 :   for (i=1;i<l;i++) gel(z,i) = Fq_mul(gel(x,i),y,T,p);
     878           0 :   return z;
     879             : }
     880             : 
     881             : GEN
     882         549 : FqV_to_FlxV(GEN v, GEN T, GEN pp)
     883             : {
     884         549 :   long j, N = lg(v);
     885         549 :   long vT = evalvarn(get_FpX_var(T));
     886         549 :   ulong p = pp[2];
     887         549 :   GEN y = cgetg(N, t_VEC);
     888        2748 :   for (j=1; j<N; j++)
     889        4398 :     gel(y,j) = (typ(gel(v,j))==t_INT?  Z_to_Flx(gel(v,j), p, vT)
     890        2199 :                                     : ZX_to_Flx(gel(v,j), p));
     891         549 :   return y;
     892             : }
     893             : 
     894             : GEN
     895        1932 : FqC_to_FlxC(GEN v, GEN T, GEN pp)
     896             : {
     897        1932 :   long j, N = lg(v);
     898        1932 :   long vT = evalvarn(get_FpX_var(T));
     899        1932 :   ulong p = pp[2];
     900        1932 :   GEN y = cgetg(N, t_COL);
     901       37688 :   for (j=1; j<N; j++)
     902       82227 :     gel(y,j) = (typ(gel(v,j))==t_INT?  Z_to_Flx(gel(v,j), p, vT)
     903       46471 :                                     : ZX_to_Flx(gel(v,j), p));
     904        1932 :   return y;
     905             : }
     906             : 
     907             : GEN
     908         392 : FqM_to_FlxM(GEN x, GEN T, GEN pp)
     909             : {
     910         392 :   long j, n = lg(x);
     911         392 :   GEN y = cgetg(n,t_MAT);
     912         392 :   if (n == 1) return y;
     913        2324 :   for (j=1; j<n; j++)
     914        1932 :     gel(y,j) = FqC_to_FlxC(gel(x,j), T, pp);
     915         392 :   return y;
     916             : }
     917             : 
     918             : 
     919             : /*******************************************************************/
     920             : /*  Isomorphisms between finite fields                             */
     921             : /*******************************************************************/
     922             : 
     923             : /* compute the reciprocical isomorphism of S mod T,p, i.e. V such that
     924             :    V(S)=X  mod T,p*/
     925             : 
     926             : GEN
     927          28 : Flxq_ffisom_inv(GEN S,GEN T, ulong p)
     928             : {
     929          28 :   pari_sp ltop = avma;
     930          28 :   long n = degpol(T);
     931          28 :   GEN M = Flxq_matrix_pow(S,n,n,T,p);
     932          28 :   GEN V = Flm_Flc_invimage(M, vecsmall_ei(n, 2), p);
     933          28 :   return gerepileupto(ltop, Flv_to_Flx(V, T[1]));
     934             : }
     935             : 
     936             : GEN
     937        1253 : FpXQ_ffisom_inv(GEN S,GEN T, GEN p)
     938             : {
     939        1253 :   pari_sp ltop = avma;
     940        1253 :   long n = degpol(T);
     941        1253 :   GEN V, M = FpXQ_matrix_pow(S,n,n,T,p);
     942        1253 :   V = FpM_FpC_invimage(M, col_ei(n, 2), p);
     943        1253 :   return gerepilecopy(ltop, RgV_to_RgX(V, varn(T)));
     944             : }
     945             : 
     946             : /* Let M the matrix of the x^p Frobenius automorphism.
     947             :  * Compute x^(p^i) for i=0..r */
     948             : static GEN
     949          28 : FpM_Frobenius(GEN M, long r, GEN p, long v)
     950             : {
     951          28 :   GEN W, V = cgetg(r+2,t_VEC);
     952             :   long i;
     953          28 :   gel(V,1) = pol_x(v); if (!r) return V;
     954          28 :   gel(V,2) = RgV_to_RgX(gel(M,2),v);
     955          28 :   W = gel(M,2);
     956         280 :   for (i = 3; i <= r+1; ++i)
     957             :   {
     958         252 :     W = FpM_FpC_mul(M,W,p);
     959         252 :     gel(V,i) = RgV_to_RgX(W,v);
     960             :   }
     961          28 :   return V;
     962             : }
     963             : 
     964             : /* Let M the matrix of the x^p Frobenius automorphism.
     965             :  * Compute x^(p^i) for i=0..r */
     966             : static GEN
     967        2226 : Flm_Frobenius(GEN M, long r, ulong p, long v)
     968             : {
     969        2226 :   GEN W, V = cgetg(r+2,t_VEC);
     970             :   long i;
     971        2226 :   gel(V,1) = polx_Flx(v); if (!r) return V;
     972        2226 :   gel(V,2) = Flv_to_Flx(gel(M,2),v);
     973        2226 :   W = gel(M,2);
     974        7462 :   for (i = 3; i <= r+1; ++i)
     975             :   {
     976        5236 :     W = Flm_Flc_mul(M,W,p);
     977        5236 :     gel(V,i) = Flv_to_Flx(W,v);
     978             :   }
     979        2226 :   return V;
     980             : }
     981             : 
     982             : /* Let P a polynomial != 0 and M the matrix of the x^p Frobenius automorphism in
     983             :  * FFp[X]/T. Compute P(M)
     984             :  * V=FpM_Frobenius(M, p, degpol(P), v)
     985             :  * not stack clean
     986             :  */
     987             : 
     988             : static GEN
     989          28 : FpXQV_FpX_Frobenius(GEN V, GEN P, GEN T, GEN p)
     990             : {
     991             :   pari_sp btop;
     992             :   long i;
     993          28 :   long l = get_FpX_degree(T);
     994          28 :   long v = get_FpX_var(T);
     995             :   GEN M,W,Mi;
     996             :   GEN *gptr[2];
     997          28 :   long lV=lg(V);
     998          28 :   GEN  PV=RgX_to_RgC(P, lgpol(P));
     999          28 :   M=cgetg(l+1,t_VEC);
    1000          28 :   gel(M,1) = scalar_ZX_shallow(FpX_eval(P,gen_1,p),v);
    1001          28 :   gel(M,2) = FpXV_FpC_mul(V,PV,p);
    1002          28 :   btop=avma;
    1003          28 :   gptr[0]=&Mi;
    1004          28 :   gptr[1]=&W;
    1005          28 :   W = leafcopy(V);
    1006         280 :   for(i=3;i<=l;i++)
    1007             :   {
    1008             :     long j;
    1009             :     pari_sp bbot;
    1010         252 :     GEN W2=cgetg(lV,t_VEC);
    1011        4816 :     for(j=1;j<lV;j++)
    1012        4564 :       gel(W2,j) = FpXQ_mul(gel(W,j),gel(V,j),T,p);
    1013         252 :     bbot=avma;
    1014         252 :     Mi=FpXV_FpC_mul(W2,PV,p);
    1015         252 :     W=gcopy(W2);
    1016         252 :     gerepilemanysp(btop,bbot,gptr,2);
    1017         252 :     btop=(pari_sp)W;
    1018         252 :     gel(M,i) = Mi;
    1019             :   }
    1020          28 :   return RgXV_to_RgM(M,l);
    1021             : }
    1022             : 
    1023             : static GEN
    1024        2226 : FlxqV_Flx_Frobenius(GEN V, GEN P, GEN T, ulong p)
    1025             : {
    1026             :   pari_sp btop;
    1027             :   long i;
    1028        2226 :   long l = get_Flx_degree(T);
    1029        2226 :   long v = get_Flx_var(T);
    1030             :   GEN M,W,Mi;
    1031        2226 :   GEN PV=Flx_to_Flv(P, lgpol(P));
    1032             :   GEN *gptr[2];
    1033        2226 :   long lV=lg(V);
    1034        2226 :   M=cgetg(l+1,t_VEC);
    1035        2226 :   gel(M,1) = Fl_to_Flx(Flx_eval(P,1,p),v);
    1036        2226 :   gel(M,2) = FlxV_Flc_mul(V,PV,p);
    1037        2226 :   btop=avma;
    1038        2226 :   gptr[0]=&Mi;
    1039        2226 :   gptr[1]=&W;
    1040        2226 :   W=gcopy(V);
    1041       12810 :   for(i=3;i<=l;i++)
    1042             :   {
    1043             :     long j;
    1044             :     pari_sp bbot;
    1045       10584 :     GEN W2=cgetg(lV,t_VEC);
    1046      196070 :     for(j=1;j<lV;j++)
    1047      185486 :       gel(W2,j) = Flxq_mul(gel(W,j),gel(V,j),T,p);
    1048       10584 :     bbot=avma;
    1049       10584 :     Mi=FlxV_Flc_mul(W2,PV,p);
    1050       10584 :     W=gcopy(W2);
    1051       10584 :     gerepilemanysp(btop,bbot,gptr,2);
    1052       10584 :     btop=(pari_sp)W;
    1053       10584 :     gel(M,i) = Mi;
    1054             :   }
    1055        2226 :   return FlxV_to_Flm(M,l);
    1056             : }
    1057             : 
    1058             : /* Let M the matrix of the Frobenius automorphism of Fp[X]/(T).
    1059             :  * Compute M^d
    1060             :  * TODO: use left-right binary (tricky!)
    1061             :  */
    1062             : GEN
    1063           0 : Flm_Frobenius_pow(GEN M, long d, GEN T, ulong p)
    1064             : {
    1065           0 :   pari_sp ltop=avma;
    1066           0 :   long i,l=degpol(T);
    1067           0 :   GEN R, W = gel(M,2);
    1068           0 :   for (i = 2; i <= d; ++i) W = Flm_Flc_mul(M,W,p);
    1069           0 :   R=Flxq_matrix_pow(Flv_to_Flx(W,T[2]),l,l,T,p);
    1070           0 :   return gerepileupto(ltop,R);
    1071             : }
    1072             : 
    1073             : GEN
    1074           0 : FpM_Frobenius_pow(GEN M, long d, GEN T, GEN p)
    1075             : {
    1076           0 :   pari_sp ltop=avma;
    1077           0 :   long i,l=degpol(T);
    1078           0 :   GEN R, W = gel(M,2);
    1079           0 :   for (i = 2; i <= d; ++i) W = FpM_FpC_mul(M,W,p);
    1080           0 :   R=FpXQ_matrix_pow(RgV_to_RgX(W,varn(T)),l,l,T,p);
    1081           0 :   return gerepilecopy(ltop,R);
    1082             : }
    1083             : 
    1084             : /* Essentially we want to compute
    1085             :  * FqM_ker(MA-pol_x(v),U,l)
    1086             :  * To avoid use of matrix in Fq we procede as follows:
    1087             :  * We compute FpM_ker(U(MA),l) and then we recover
    1088             :  * the eigen value by Galois action, see formula.
    1089             :  */
    1090             : 
    1091             : static GEN
    1092        2226 : Flx_intersect_ker(GEN P, GEN MA, GEN U, ulong p)
    1093             : {
    1094        2226 :   pari_sp ltop = avma;
    1095        2226 :   long i, vp = P[1], vu = U[1], r = degpol(U);
    1096             :   GEN A, R;
    1097             :   ulong ib0;
    1098             :   pari_timer T;
    1099             :   GEN M, V;
    1100        2226 :   if (DEBUGLEVEL>=4) timer_start(&T);
    1101        2226 :   V = Flm_Frobenius(MA, r, p, U[1]);
    1102        2226 :   if (DEBUGLEVEL>=4) timer_printf(&T,"pol[Frobenius]");
    1103        2226 :   M = FlxqV_Flx_Frobenius(V, U, P, p);
    1104        2226 :   if (p==2)
    1105        1386 :     A = F2m_to_Flm(F2m_ker(Flm_to_F2m(M)));
    1106             :   else
    1107         840 :     A = Flm_ker(M,p);
    1108        2226 :   if (DEBUGLEVEL>=4) timer_printf(&T,"matrix polcyclo");
    1109        2226 :   if (lg(A)!=r+1) pari_err_IRREDPOL("FpX_ffintersect", Flx_to_ZX(P));
    1110        2226 :   A = gerepileupto(ltop,A);
    1111             :   /*The formula is
    1112             :    * a_{r-1} = -\phi(a_0)/b_0
    1113             :    * a_{i-1} = \phi(a_i)+b_ia_{r-1}  i=r-1 to 1
    1114             :    * Where a_0=A[1] and b_i=U[i+2] */
    1115        2226 :   ib0 = Fl_inv(Fl_neg(U[2], p), p);
    1116        2226 :   R = cgetg(r+1,t_MAT);
    1117        2226 :   gel(R,1) = gel(A,1);
    1118        2226 :   gel(R,r) = Flm_Flc_mul(MA, Flv_Fl_mul(gel(A,1),ib0, p), p);
    1119        5236 :   for(i=r-1; i>1; i--)
    1120             :   {
    1121        3010 :     gel(R,i) = Flm_Flc_mul(MA,gel(R,i+1),p);
    1122        3010 :     Flv_add_inplace(gel(R,i), Flv_Fl_mul(gel(R,r), U[i+2], p), p);
    1123             :   }
    1124        2226 :   return gerepileupto(ltop, Flm_to_FlxX(Flm_transpose(R),vp,vu));
    1125             : }
    1126             : 
    1127             : static GEN
    1128        1442 : FpX_intersect_ker(GEN P, GEN MA, GEN U, GEN l)
    1129             : {
    1130        1442 :   pari_sp ltop = avma;
    1131        1442 :   long i, vp = varn(P), vu = varn(U), r = degpol(U);
    1132             :   GEN V, A, R, ib0;
    1133             :   pari_timer T;
    1134        1442 :   if (lgefint(l)==3)
    1135             :   {
    1136        1414 :     ulong p = l[2];
    1137        1414 :     GEN res = Flx_intersect_ker(ZX_to_Flx(P,p), ZM_to_Flm(MA,p), ZX_to_Flx(U,p), p);
    1138        1414 :     return gerepileupto(ltop, FlxX_to_ZXX(res));
    1139             :   }
    1140          28 :   if (DEBUGLEVEL>=4) timer_start(&T);
    1141          28 :   V = FpM_Frobenius(MA,r,l,vu);
    1142          28 :   if (DEBUGLEVEL>=4) timer_printf(&T,"pol[Frobenius]");
    1143          28 :   A = FpM_ker(FpXQV_FpX_Frobenius(V, U, P, l), l);
    1144          28 :   if (DEBUGLEVEL>=4) timer_printf(&T,"matrix polcyclo");
    1145          28 :   if (lg(A)!=r+1) pari_err_IRREDPOL("FpX_ffintersect", P);
    1146          28 :   A = gerepileupto(ltop,A);
    1147             :   /*The formula is
    1148             :    * a_{r-1} = -\phi(a_0)/b_0
    1149             :    * a_{i-1} = \phi(a_i)+b_ia_{r-1}  i=r-1 to 1
    1150             :    * Where a_0=A[1] and b_i=U[i+2] */
    1151          28 :   ib0 = Fp_inv(negi(gel(U,2)),l);
    1152          28 :   R = cgetg(r+1,t_MAT);
    1153          28 :   gel(R,1) = gel(A,1);
    1154          28 :   gel(R,r) = FpM_FpC_mul(MA, FpC_Fp_mul(gel(A,1),ib0,l), l);
    1155         252 :   for(i=r-1;i>1;i--)
    1156         672 :     gel(R,i) = FpC_add(FpM_FpC_mul(MA,gel(R,i+1),l),
    1157         448 :         FpC_Fp_mul(gel(R,r), gel(U,i+2), l),l);
    1158          28 :   return gerepilecopy(ltop,RgM_to_RgXX(shallowtrans(R),vp,vu));
    1159             : }
    1160             : 
    1161             : /* n must divide both the degree of P and Q.  Compute SP and SQ such
    1162             :   that the subfield of FF_l[X]/(P) generated by SP and the subfield of
    1163             :   FF_l[X]/(Q) generated by SQ are isomorphic of degree n.  P and Q do
    1164             :   not need to be of the same variable.  if MA (resp. MB) is not NULL,
    1165             :   must be the matrix of the Frobenius map in FF_l[X]/(P) (resp.
    1166             :   FF_l[X]/(Q) ).  */
    1167             : /* Note on the implementation choice:
    1168             :  * We assume the prime p is very large
    1169             :  * so we handle Frobenius as matrices.
    1170             :  */
    1171             : 
    1172             : void
    1173        3332 : Flx_ffintersect(GEN P, GEN Q, long n, ulong l,GEN *SP, GEN *SQ, GEN MA, GEN MB)
    1174             : {
    1175        3332 :   pari_sp ltop = avma;
    1176        3332 :   long vp = P[1], vq = Q[1], np = degpol(P), nq = degpol(Q), e;
    1177             :   ulong pg;
    1178             :   GEN A, B, Ap, Bp;
    1179        3332 :   if (np<=0) pari_err_IRREDPOL("FpX_ffintersect", P);
    1180        3332 :   if (nq<=0) pari_err_IRREDPOL("FpX_ffintersect", Q);
    1181        3332 :   if (n<=0 || np%n || nq%n)
    1182           0 :     pari_err_TYPE("FpX_ffintersect [bad degrees]",stoi(n));
    1183        3332 :   e = u_lvalrem(n, l, &pg);
    1184        3332 :   if(!MA) MA = Flx_matFrobenius(P,l);
    1185        3332 :   if(!MB) MB = Flx_matFrobenius(Q,l);
    1186        3332 :   A = Ap = pol0_Flx(vp);
    1187        3332 :   B = Bp = pol0_Flx(vq);
    1188        3332 :   if (pg > 1)
    1189             :   {
    1190             :     pari_timer T;
    1191        2128 :     GEN ipg = utoipos(pg);
    1192        2128 :     if (l%pg == 1)
    1193             :     /* No need to use relative extension, so don't. (Well, now we don't
    1194             :      * in the other case either, but this special case is more efficient) */
    1195             :     {
    1196             :       GEN L;
    1197             :       ulong z, An, Bn;
    1198        1722 :       z = Fl_neg(rootsof1_Fl(pg, l), l);
    1199        1722 :       if (DEBUGLEVEL>=4) timer_start(&T);
    1200        1722 :       A = Flm_ker(Flm_Fl_add(MA, z, l),l);
    1201        1722 :       if (lg(A)!=2) pari_err_IRREDPOL("FpX_ffintersect",P);
    1202        1722 :       A = Flv_to_Flx(gel(A,1),vp);
    1203             : 
    1204        1722 :       B = Flm_ker(Flm_Fl_add(MB, z, l),l);
    1205        1722 :       if (lg(B)!=2) pari_err_IRREDPOL("FpX_ffintersect",Q);
    1206        1722 :       B = Flv_to_Flx(gel(B,1),vq);
    1207             : 
    1208        1722 :       if (DEBUGLEVEL>=4) timer_printf(&T, "FpM_ker");
    1209        1722 :       An = Flxq_powu(A,pg,P,l)[2];
    1210        1722 :       Bn = Flxq_powu(B,pg,Q,l)[2];
    1211        1722 :       if (!Bn) pari_err_IRREDPOL("FpX_ffintersect", mkvec2(P,Q));
    1212        1722 :       z = Fl_div(An,Bn,l);
    1213        1722 :       L = Fp_sqrtn(utoi(z),ipg,utoipos(l),NULL);
    1214        1722 :       if (!L) pari_err_IRREDPOL("FpX_ffintersect", mkvec2(P,Q));
    1215        1722 :       if (DEBUGLEVEL>=4) timer_printf(&T, "Fp_sqrtn");
    1216        1722 :       B = Flx_Fl_mul(B,itou(L),l);
    1217             :     }
    1218             :     else
    1219             :     {
    1220             :       GEN L, An, Bn, z, U;
    1221         406 :       U = gmael(Flx_factor(ZX_to_Flx(polcyclo(pg, fetch_var()),l),l),1,1);
    1222         406 :       A = Flx_intersect_ker(P, MA, U, l);
    1223         406 :       B = Flx_intersect_ker(Q, MB, U, l);
    1224         406 :       if (DEBUGLEVEL>=4) timer_start(&T);
    1225         406 :       An = gel(FlxYqq_pow(A,ipg,P,U,l),2);
    1226         406 :       Bn = gel(FlxYqq_pow(B,ipg,Q,U,l),2);
    1227         406 :       if (DEBUGLEVEL>=4) timer_printf(&T,"pows [P,Q]");
    1228         406 :       z = Flxq_div(An,Bn,U,l);
    1229         406 :       L = Flxq_sqrtn(z,ipg,U,l,NULL);
    1230         406 :       if (!L) pari_err_IRREDPOL("FpX_ffintersect", mkvec2(P,Q));
    1231         406 :       if (DEBUGLEVEL>=4) timer_printf(&T,"FpXQ_sqrtn");
    1232         406 :       B = FlxqX_Flxq_mul(B,L,U,l);
    1233         406 :       A = FlxY_evalx(A,0,l);
    1234         406 :       B = FlxY_evalx(B,0,l);
    1235         406 :       (void)delete_var();
    1236             :     }
    1237             :   }
    1238        3332 :   if (e)
    1239             :   {
    1240             :     GEN VP, VQ, Ay, By;
    1241        1253 :     ulong lmun = l-1;
    1242             :     long j;
    1243        1253 :     MA = Flm_Fl_add(MA,lmun,l);
    1244        1253 :     MB = Flm_Fl_add(MB,lmun,l);
    1245        1253 :     Ay = pol1_Flx(vp);
    1246        1253 :     By = pol1_Flx(vq);
    1247        1253 :     VP = vecsmall_ei(np, 1);
    1248        1253 :     VQ = np == nq? VP: vecsmall_ei(nq, 1); /* save memory */
    1249        2681 :     for(j=0;j<e;j++)
    1250             :     {
    1251        1428 :       if (j)
    1252             :       {
    1253         175 :         Ay = Flxq_mul(Ay,Flxq_powu(Ap,lmun,P,l),P,l);
    1254         175 :         VP = Flx_to_Flv(Ay,np);
    1255             :       }
    1256        1428 :       Ap = Flm_Flc_invimage(MA,VP,l);
    1257        1428 :       Ap = Flv_to_Flx(Ap,vp);
    1258             : 
    1259        1428 :       if (j)
    1260             :       {
    1261         175 :         By = Flxq_mul(By,Flxq_powu(Bp,lmun,Q,l),Q,l);
    1262         175 :         VQ = Flx_to_Flv(By,nq);
    1263             :       }
    1264        1428 :       Bp = Flm_Flc_invimage(MB,VQ,l);
    1265        1428 :       Bp = Flv_to_Flx(Bp,vq);
    1266             :     }
    1267             :   }
    1268        3332 :   *SP = Flx_add(A,Ap,l);
    1269        3332 :   *SQ = Flx_add(B,Bp,l);
    1270        3332 :   gerepileall(ltop,2,SP,SQ);
    1271        3332 : }
    1272             : 
    1273             : /* Let l be a prime number, P, Q in ZZ[X].  P and Q are both
    1274             :  * irreducible modulo l and degree(P) divides degree(Q).  Output a
    1275             :  * monomorphism between FF_l[X]/(P) and FF_l[X]/(Q) as a polynomial R such
    1276             :  * that Q | P(R) mod l.  If P and Q have the same degree, it is of course an
    1277             :  * isomorphism.  */
    1278             : GEN
    1279          28 : Flx_ffisom(GEN P,GEN Q,ulong l)
    1280             : {
    1281          28 :   pari_sp av = avma;
    1282             :   GEN SP, SQ, R;
    1283          28 :   Flx_ffintersect(P,Q,degpol(P),l,&SP,&SQ,NULL,NULL);
    1284          28 :   R = Flxq_ffisom_inv(SP,P,l);
    1285          28 :   return gerepileupto(av, Flx_Flxq_eval(R,SQ,Q,l));
    1286             : }
    1287             : 
    1288             : void
    1289        1148 : FpX_ffintersect(GEN P, GEN Q, long n, GEN l, GEN *SP, GEN *SQ, GEN MA, GEN MB)
    1290             : {
    1291        1148 :   pari_sp ltop = avma;
    1292             :   long vp, vq, np, nq, e;
    1293             :   ulong pg;
    1294             :   GEN A, B, Ap, Bp;
    1295        1148 :   vp = varn(P); np = degpol(P);
    1296        1148 :   vq = varn(Q); nq = degpol(Q);
    1297        1148 :   if (np<=0) pari_err_IRREDPOL("FpX_ffintersect", P);
    1298        1148 :   if (nq<=0) pari_err_IRREDPOL("FpX_ffintersect", Q);
    1299        1148 :   if (n<=0 || np%n || nq%n)
    1300           0 :     pari_err_TYPE("FpX_ffintersect [bad degrees]",stoi(n));
    1301        1148 :   e = u_pvalrem(n, l, &pg);
    1302        1148 :   if(!MA) MA = FpX_matFrobenius(P, l);
    1303        1148 :   if(!MB) MB = FpX_matFrobenius(Q, l);
    1304        1148 :   A = Ap = pol_0(vp);
    1305        1148 :   B = Bp = pol_0(vq);
    1306        1148 :   if (pg > 1)
    1307             :   {
    1308        1022 :     GEN ipg = utoipos(pg);
    1309             :     pari_timer T;
    1310        1022 :     if (umodiu(l,pg) == 1)
    1311             :     /* No need to use relative extension, so don't. (Well, now we don't
    1312             :      * in the other case either, but this special case is more efficient) */
    1313             :     {
    1314             :       GEN L, An, Bn, z;
    1315         301 :       z = negi( rootsof1u_Fp(pg, l) );
    1316         301 :       if (DEBUGLEVEL>=4) timer_start(&T);
    1317         301 :       A = FpM_ker(RgM_Rg_add_shallow(MA, z),l);
    1318         301 :       if (lg(A)!=2) pari_err_IRREDPOL("FpX_ffintersect",P);
    1319         301 :       A = RgV_to_RgX(gel(A,1),vp);
    1320             : 
    1321         301 :       B = FpM_ker(RgM_Rg_add_shallow(MB, z),l);
    1322         301 :       if (lg(B)!=2) pari_err_IRREDPOL("FpX_ffintersect",Q);
    1323         301 :       B = RgV_to_RgX(gel(B,1),vq);
    1324             : 
    1325         301 :       if (DEBUGLEVEL>=4) timer_printf(&T, "FpM_ker");
    1326         301 :       An = gel(FpXQ_pow(A,ipg,P,l),2);
    1327         301 :       Bn = gel(FpXQ_pow(B,ipg,Q,l),2);
    1328         301 :       if (!signe(Bn)) pari_err_IRREDPOL("FpX_ffintersect", mkvec2(P,Q));
    1329         301 :       z = Fp_div(An,Bn,l);
    1330         301 :       L = Fp_sqrtn(z,ipg,l,NULL);
    1331         301 :       if (!L) pari_err_IRREDPOL("FpX_ffintersect", mkvec2(P,Q));
    1332         301 :       if (DEBUGLEVEL>=4) timer_printf(&T, "Fp_sqrtn");
    1333         301 :       B = FpX_Fp_mul(B,L,l);
    1334             :     }
    1335             :     else
    1336             :     {
    1337             :       GEN L, An, Bn, z, U;
    1338         721 :       U = gmael(FpX_factor(polcyclo(pg,fetch_var()),l),1,1);
    1339         721 :       A = FpX_intersect_ker(P, MA, U, l);
    1340         721 :       B = FpX_intersect_ker(Q, MB, U, l);
    1341         721 :       if (DEBUGLEVEL>=4) timer_start(&T);
    1342         721 :       An = gel(FpXYQQ_pow(A,ipg,P,U,l),2);
    1343         721 :       Bn = gel(FpXYQQ_pow(B,ipg,Q,U,l),2);
    1344         721 :       if (DEBUGLEVEL>=4) timer_printf(&T,"pows [P,Q]");
    1345         721 :       if (!signe(Bn)) pari_err_IRREDPOL("FpX_ffintersect", mkvec2(P,Q));
    1346         721 :       z = Fq_div(An,Bn,U,l);
    1347         721 :       L = Fq_sqrtn(z,ipg,U,l,NULL);
    1348         721 :       if (!L) pari_err_IRREDPOL("FpX_ffintersect", mkvec2(P,Q));
    1349         721 :       if (DEBUGLEVEL>=4) timer_printf(&T,"FpXQ_sqrtn");
    1350         721 :       B = FqX_Fq_mul(B,L,U,l);
    1351         721 :       A = FpXY_evalx(A,gen_0,l);
    1352         721 :       B = FpXY_evalx(B,gen_0,l);
    1353         721 :       (void)delete_var();
    1354             :     }
    1355             :   }
    1356        1148 :   if (e)
    1357             :   {
    1358         133 :     GEN VP, VQ, Ay, By, lmun = addis(l,-1);
    1359             :     long j;
    1360         133 :     MA = RgM_Rg_add_shallow(MA,gen_m1);
    1361         133 :     MB = RgM_Rg_add_shallow(MB,gen_m1);
    1362         133 :     Ay = pol_1(vp);
    1363         133 :     By = pol_1(vq);
    1364         133 :     VP = col_ei(np, 1);
    1365         133 :     VQ = np == nq? VP: col_ei(nq, 1); /* save memory */
    1366         343 :     for(j=0;j<e;j++)
    1367             :     {
    1368         210 :       if (j)
    1369             :       {
    1370          77 :         Ay = FpXQ_mul(Ay,FpXQ_pow(Ap,lmun,P,l),P,l);
    1371          77 :         VP = RgX_to_RgC(Ay,np);
    1372             :       }
    1373         210 :       Ap = FpM_FpC_invimage(MA,VP,l);
    1374         210 :       Ap = RgV_to_RgX(Ap,vp);
    1375             : 
    1376         210 :       if (j)
    1377             :       {
    1378          77 :         By = FpXQ_mul(By,FpXQ_pow(Bp,lmun,Q,l),Q,l);
    1379          77 :         VQ = RgX_to_RgC(By,nq);
    1380             :       }
    1381         210 :       Bp = FpM_FpC_invimage(MB,VQ,l);
    1382         210 :       Bp = RgV_to_RgX(Bp,vq);
    1383             :     }
    1384             :   }
    1385        1148 :   *SP = FpX_add(A,Ap,l);
    1386        1148 :   *SQ = FpX_add(B,Bp,l);
    1387        1148 :   gerepileall(ltop,2,SP,SQ);
    1388        1148 : }
    1389             : /* Let l be a prime number, P, Q in ZZ[X].  P and Q are both
    1390             :  * irreducible modulo l and degree(P) divides degree(Q).  Output a
    1391             :  * monomorphism between FF_l[X]/(P) and FF_l[X]/(Q) as a polynomial R such
    1392             :  * that Q | P(R) mod l.  If P and Q have the same degree, it is of course an
    1393             :  * isomorphism.  */
    1394             : GEN
    1395        1127 : FpX_ffisom(GEN P,GEN Q,GEN l)
    1396             : {
    1397        1127 :   pari_sp av = avma;
    1398             :   GEN SP, SQ, R;
    1399        1127 :   FpX_ffintersect(P,Q,degpol(P),l,&SP,&SQ,NULL,NULL);
    1400        1127 :   R = FpXQ_ffisom_inv(SP,P,l);
    1401        1127 :   return gerepileupto(av, FpX_FpXQ_eval(R,SQ,Q,l));
    1402             : }
    1403             : 
    1404             : /* Let l be a prime number, P a ZX irreducible modulo l, MP the matrix of the
    1405             :  * Frobenius automorphism of F_l[X]/(P).
    1406             :  * Factor P over the subfield of F_l[X]/(P) of index d. */
    1407             : static GEN
    1408          21 : FpX_factorgalois(GEN P, GEN l, long d, long w, GEN MP)
    1409             : {
    1410          21 :   pari_sp ltop = avma;
    1411             :   GEN R, V, Tl, z, M;
    1412          21 :   long k, n = degpol(P), m = n/d;
    1413          21 :   long v = varn(P);
    1414             : 
    1415             :   /* x - y */
    1416          21 :   if (m == 1) return deg1pol_shallow(gen_1, deg1pol_shallow(subis(l,1), gen_0, w), v);
    1417           0 :   M = FpM_Frobenius_pow(MP,d,P,l);
    1418             : 
    1419           0 :   Tl = leafcopy(P); setvarn(Tl,w);
    1420           0 :   V = cgetg(m+1,t_VEC);
    1421           0 :   gel(V,1) = pol_x(w);
    1422           0 :   z = gel(M,2);
    1423           0 :   gel(V,2) = RgV_to_RgX(z,w);
    1424           0 :   for(k=3;k<=m;k++)
    1425             :   {
    1426           0 :     z = FpM_FpC_mul(M,z,l);
    1427           0 :     gel(V,k) = RgV_to_RgX(z,w);
    1428             :   }
    1429           0 :   R = FqV_roots_to_pol(V,Tl,l,v);
    1430           0 :   return gerepileupto(ltop,R);
    1431             : }
    1432             : /* same: P is an Flx, MP an Flm */
    1433             : static GEN
    1434        3304 : Flx_factorgalois(GEN P, ulong l, long d, long w, GEN MP)
    1435             : {
    1436        3304 :   pari_sp ltop = avma;
    1437             :   GEN R, V, Tl, z, M;
    1438        3304 :   long k, n = degpol(P), m = n/d;
    1439        3304 :   long v = P[1];
    1440             : 
    1441        3304 :   if (m == 1) {
    1442        3304 :     R = polx_Flx(v);
    1443        3304 :     gel(R,2) = z = polx_Flx(w); z[3] = l - 1; /* - y */
    1444        3304 :     gel(R,3) = pol1_Flx(w);
    1445        3304 :     return R; /* x - y */
    1446             :   }
    1447           0 :   M = Flm_Frobenius_pow(MP,d,P,l);
    1448             : 
    1449           0 :   Tl = leafcopy(P); setvarn(Tl,w);
    1450           0 :   V = cgetg(m+1,t_VEC);
    1451           0 :   gel(V,1) = polx_Flx(w);
    1452           0 :   z = gel(M,2);
    1453           0 :   gel(V,2) = Flv_to_Flx(z,w);
    1454           0 :   for(k=3;k<=m;k++)
    1455             :   {
    1456           0 :     z = Flm_Flc_mul(M,z,l);
    1457           0 :     gel(V,k) = Flv_to_Flx(z,w);
    1458             :   }
    1459           0 :   R = FlxqV_roots_to_pol(V,Tl,l,v);
    1460           0 :   return gerepileupto(ltop,R);
    1461             : }
    1462             : 
    1463             : GEN
    1464       24801 : Flx_factorff_irred(GEN P, GEN Q, ulong p)
    1465             : {
    1466       24801 :   pari_sp ltop = avma, av;
    1467             :   GEN SP, SQ, MP, MQ, M, FP, FQ, E, V, IR, res;
    1468       24801 :   long np = degpol(P), nq = degpol(Q), d = cgcd(np,nq);
    1469       24801 :   long i, vp = P[1], vq = Q[1];
    1470       24801 :   if (d==1) retmkcol(Flx_to_FlxX(P, vq));
    1471        3304 :   FQ = Flx_matFrobenius(Q,p);
    1472        3304 :   av = avma;
    1473        3304 :   FP = Flx_matFrobenius(P,p);
    1474        3304 :   Flx_ffintersect(P,Q,d,p,&SP,&SQ, FP, FQ);
    1475        3304 :   E = Flx_factorgalois(P,p,d,vq, FP);
    1476        3304 :   E = FlxX_to_Flm(E,np);
    1477        3304 :   MP= Flxq_matrix_pow(SP,np,d,P,p);
    1478        3304 :   IR= gel(Flm_indexrank(MP,p),1);
    1479        3304 :   E = rowpermute(E, IR);
    1480        3304 :   M = rowpermute(MP,IR);
    1481        3304 :   M = Flm_inv(M,p);
    1482        3304 :   MQ= Flxq_matrix_pow(SQ,nq,d,Q,p);
    1483        3304 :   M = Flm_mul(MQ,M,p);
    1484        3304 :   M = Flm_mul(M,E,p);
    1485        3304 :   M = gerepileupto(av,M);
    1486        3304 :   V = cgetg(d+1,t_VEC);
    1487        3304 :   gel(V,1) = M;
    1488       12145 :   for(i=2;i<=d;i++)
    1489        8841 :     gel(V,i) = Flm_mul(FQ,gel(V,i-1),p);
    1490        3304 :   res = cgetg(d+1,t_COL);
    1491       15449 :   for(i=1;i<=d;i++)
    1492       12145 :     gel(res,i) = Flm_to_FlxX(gel(V,i),vp,vq);
    1493        3304 :   return gerepileupto(ltop,res);
    1494             : }
    1495             : 
    1496             : /* P,Q irreducible over F_p. Factor P over FF_p[X] / Q  [factors are ordered as
    1497             :  * a Frobenius cycle] */
    1498             : GEN
    1499         973 : FpX_factorff_irred(GEN P, GEN Q, GEN p)
    1500             : {
    1501         973 :   pari_sp ltop = avma, av;
    1502             :   GEN res;
    1503         973 :   long np = degpol(P), nq = degpol(Q), d = cgcd(np,nq);
    1504         973 :   if (d==1) return mkcolcopy(P);
    1505             : 
    1506         700 :   if (lgefint(p)==3)
    1507             :   {
    1508         679 :     ulong pp = p[2];
    1509         679 :     GEN F = Flx_factorff_irred(ZX_to_Flx(P,pp), ZX_to_Flx(Q,pp), pp);
    1510         679 :     long i, lF = lg(F);
    1511         679 :     res = cgetg(lF, t_COL);
    1512        6748 :     for(i=1; i<lF; i++)
    1513        6069 :       gel(res,i) = FlxX_to_ZXX(gel(F,i));
    1514             :   }
    1515             :   else
    1516             :   {
    1517             :     GEN SP, SQ, MP, MQ, M, FP, FQ, E, V, IR;
    1518          21 :     long i, vp = varn(P), vq = varn(Q);
    1519          21 :     FQ = FpX_matFrobenius(Q,p);
    1520          21 :     av = avma;
    1521          21 :     FP = FpX_matFrobenius(P,p);
    1522          21 :     FpX_ffintersect(P,Q,d,p,&SP,&SQ,FP,FQ);
    1523             : 
    1524          21 :     E = FpX_factorgalois(P,p,d,vq,FP);
    1525          21 :     E = RgXX_to_RgM(E,np);
    1526          21 :     MP= FpXQ_matrix_pow(SP,np,d,P,p);
    1527          21 :     IR= gel(FpM_indexrank(MP,p),1);
    1528          21 :     E = rowpermute(E, IR);
    1529          21 :     M = rowpermute(MP,IR);
    1530          21 :     M = FpM_inv(M,p);
    1531          21 :     MQ= FpXQ_matrix_pow(SQ,nq,d,Q,p);
    1532          21 :     M = FpM_mul(MQ,M,p);
    1533          21 :     M = FpM_mul(M,E,p);
    1534          21 :     M = gerepileupto(av,M);
    1535          21 :     V = cgetg(d+1,t_VEC);
    1536          21 :     gel(V,1) = M;
    1537         287 :     for(i=2;i<=d;i++)
    1538         266 :       gel(V,i) = FpM_mul(FQ,gel(V,i-1),p);
    1539          21 :     res = cgetg(d+1,t_COL);
    1540         308 :     for(i=1;i<=d;i++)
    1541         287 :       gel(res,i) = RgM_to_RgXX(gel(V,i),vp,vq);
    1542             :   }
    1543         700 :   return gerepilecopy(ltop,res);
    1544             : }
    1545             : 
    1546             : /*******************************************************************/
    1547             : /*                                                                 */
    1548             : /*                          MODULAR GCD                            */
    1549             : /*                                                                 */
    1550             : /*******************************************************************/
    1551             : /* return z = a mod q, b mod p (p,q) = 1. qinv = 1/q mod p */
    1552             : static GEN
    1553    18183445 : Fl_chinese_coprime(GEN a, ulong b, GEN q, ulong p, ulong qinv, GEN pq)
    1554             : {
    1555    18183445 :   ulong d, amod = umodiu(a, p);
    1556    18183445 :   pari_sp av = avma;
    1557             :   GEN ax;
    1558             : 
    1559    18183445 :   if (b == amod) return NULL;
    1560     2384581 :   d = (b > amod)? b - amod: p - (amod - b); /* (b - a) mod p */
    1561     2384581 :   (void)new_chunk(lgefint(pq)<<1); /* HACK */
    1562     2384581 :   ax = mului(Fl_mul(d,qinv,p), q); /* d mod p, 0 mod q */
    1563     2384581 :   avma = av; return addii(a, ax); /* in ]-q, pq[ assuming a in -]-q,q[ */
    1564             : }
    1565             : GEN
    1566       17325 : Z_init_CRT(ulong Hp, ulong p) { return stoi(Fl_center(Hp, p, p>>1)); }
    1567             : GEN
    1568     3114337 : ZX_init_CRT(GEN Hp, ulong p, long v)
    1569             : {
    1570     3114337 :   long i, l = lg(Hp), lim = (long)(p>>1);
    1571     3114337 :   GEN H = cgetg(l, t_POL);
    1572     3114337 :   H[1] = evalsigne(1) | evalvarn(v);
    1573    11019296 :   for (i=2; i<l; i++)
    1574     7904959 :     gel(H,i) = stoi(Fl_center(Hp[i], p, lim));
    1575     3114337 :   return H;
    1576             : }
    1577             : 
    1578             : /* assume lg(Hp) > 1 */
    1579             : GEN
    1580      105877 : ZM_init_CRT(GEN Hp, ulong p)
    1581             : {
    1582      105877 :   long i,j, m = lgcols(Hp), l = lg(Hp), lim = (long)(p>>1);
    1583      105877 :   GEN c,cp,H = cgetg(l, t_MAT);
    1584      871690 :   for (j=1; j<l; j++)
    1585             :   {
    1586      765813 :     cp = gel(Hp,j);
    1587      765813 :     c = cgetg(m, t_COL);
    1588      765813 :     gel(H,j) = c;
    1589      765813 :     for (i=1; i<m; i++) gel(c,i) = stoi(Fl_center(cp[i],p, lim));
    1590             :   }
    1591      105877 :   return H;
    1592             : }
    1593             : 
    1594             : int
    1595       58488 : Z_incremental_CRT(GEN *H, ulong Hp, GEN *ptq, ulong p)
    1596             : {
    1597       58488 :   GEN h, q = *ptq, qp = muliu(q,p), lim = shifti(qp,-1);
    1598       58488 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1599       58488 :   int stable = 1;
    1600       58488 :   h = Fl_chinese_coprime(*H,Hp,q,p,qinv,qp);
    1601       58488 :   if (h)
    1602             :   {
    1603       14573 :     if (cmpii(h,lim) > 0) h = subii(h,qp);
    1604       14573 :     *H = h; stable = 0;
    1605             :   }
    1606       58488 :   *ptq = qp; return stable;
    1607             : }
    1608             : 
    1609             : static int
    1610     3126563 : ZX_incremental_CRT_raw(GEN *ptH, GEN Hp, GEN q, GEN qp, ulong p)
    1611             : {
    1612     3126563 :   GEN H = *ptH, h, lim = shifti(qp,-1);
    1613     3126563 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1614     3126563 :   long i, l = lg(H), lp = lg(Hp);
    1615     3126563 :   int stable = 1;
    1616             : 
    1617     3126563 :   if (l < lp)
    1618             :   { /* degree increases */
    1619           0 :     GEN x = cgetg(lp, t_POL);
    1620           0 :     for (i=1; i<l; i++)  x[i] = H[i];
    1621           0 :     for (   ; i<lp; i++) gel(x,i) = gen_0;
    1622           0 :     *ptH = H = x;
    1623           0 :     stable = 0;
    1624     3126563 :   } else if (l > lp)
    1625             :   { /* degree decreases */
    1626           0 :     GEN x = cgetg(l, t_VECSMALL);
    1627           0 :     for (i=1; i<lp; i++)  x[i] = Hp[i];
    1628           0 :     for (   ; i<l; i++) x[i] = 0;
    1629           0 :     Hp = x; lp = l;
    1630             :   }
    1631    11397217 :   for (i=2; i<lp; i++)
    1632             :   {
    1633     8270654 :     h = Fl_chinese_coprime(gel(H,i),Hp[i],q,p,qinv,qp);
    1634     8270654 :     if (h)
    1635             :     {
    1636      215885 :       if (cmpii(h,lim) > 0) h = subii(h,qp);
    1637      215885 :       gel(H,i) = h; stable = 0;
    1638             :     }
    1639             :   }
    1640     3126563 :   return stable;
    1641             : }
    1642             : 
    1643             : int
    1644     3117200 : ZX_incremental_CRT(GEN *ptH, GEN Hp, GEN *ptq, ulong p)
    1645             : {
    1646     3117200 :   GEN q = *ptq, qp = muliu(q,p);
    1647     3117200 :   int stable = ZX_incremental_CRT_raw(ptH, Hp, q, qp, p);
    1648     3117200 :   *ptq = qp; return stable;
    1649             : }
    1650             : 
    1651             : int
    1652      111483 : ZM_incremental_CRT(GEN *pH, GEN Hp, GEN *ptq, ulong p)
    1653             : {
    1654      111483 :   GEN h, H = *pH, q = *ptq, qp = muliu(q, p), lim = shifti(qp,-1);
    1655      111483 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1656      111483 :   long i,j, l = lg(H), m = lgcols(H);
    1657      111483 :   int stable = 1;
    1658      983449 :   for (j=1; j<l; j++)
    1659    10726269 :     for (i=1; i<m; i++)
    1660             :     {
    1661     9854303 :       h = Fl_chinese_coprime(gcoeff(H,i,j), coeff(Hp,i,j),q,p,qinv,qp);
    1662     9854303 :       if (h)
    1663             :       {
    1664     2154123 :         if (cmpii(h,lim) > 0) h = subii(h,qp);
    1665     2154123 :         gcoeff(H,i,j) = h; stable = 0;
    1666             :       }
    1667             :     }
    1668      111483 :   *ptq = qp; return stable;
    1669             : }
    1670             : 
    1671             : /* record the degrees of Euclidean remainders (make them as large as
    1672             :  * possible : smaller values correspond to a degenerate sequence) */
    1673             : static void
    1674        1561 : Flx_resultant_set_dglist(GEN a, GEN b, GEN dglist, ulong p)
    1675             : {
    1676             :   long da,db,dc, ind;
    1677        1561 :   pari_sp av = avma;
    1678             : 
    1679        1561 :   if (lgpol(a)==0 || lgpol(b)==0) return;
    1680        1561 :   da = degpol(a);
    1681        1561 :   db = degpol(b);
    1682        1561 :   if (db > da)
    1683           0 :   { swapspec(a,b, da,db); }
    1684        1561 :   else if (!da) return;
    1685        1561 :   ind = 0;
    1686        9814 :   while (db)
    1687             :   {
    1688        6692 :     GEN c = Flx_rem(a,b, p);
    1689        6692 :     a = b; b = c; dc = degpol(c);
    1690        6692 :     if (dc < 0) break;
    1691             : 
    1692        6692 :     ind++;
    1693        6692 :     if (dc > dglist[ind]) dglist[ind] = dc;
    1694        6692 :     if (gc_needed(av,2))
    1695             :     {
    1696           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant_all");
    1697           0 :       gerepileall(av, 2, &a,&b);
    1698             :     }
    1699        6692 :     db = dc; /* = degpol(b) */
    1700             :   }
    1701        1561 :   if (ind+1 > lg(dglist)) setlg(dglist,ind+1);
    1702        1561 :   avma = av; return;
    1703             : }
    1704             : /* assuming the PRS finishes on a degree 1 polynomial C0 + C1X, with
    1705             :  * "generic" degree sequence as given by dglist, set *Ci and return
    1706             :  * resultant(a,b). Modular version of Collins's subresultant */
    1707             : static ulong
    1708        6996 : Flx_resultant_all(GEN a, GEN b, long *C0, long *C1, GEN dglist, ulong p)
    1709             : {
    1710             :   long da,db,dc, ind;
    1711        6996 :   ulong lb, res, g = 1UL, h = 1UL, ca = 1UL, cb = 1UL;
    1712        6996 :   int s = 1;
    1713        6996 :   pari_sp av = avma;
    1714             : 
    1715        6996 :   *C0 = 1; *C1 = 0;
    1716        6996 :   if (lgpol(a)==0 || lgpol(b)==0) return 0;
    1717        6996 :   da = degpol(a);
    1718        6996 :   db = degpol(b);
    1719        6996 :   if (db > da)
    1720             :   {
    1721           0 :     swapspec(a,b, da,db);
    1722           0 :     if (both_odd(da,db)) s = -s;
    1723             :   }
    1724        6996 :   else if (!da) return 1; /* = a[2] ^ db, since 0 <= db <= da = 0 */
    1725        6996 :   ind = 0;
    1726       40727 :   while (db)
    1727             :   { /* sub-resultant algo., applied to ca * a and cb * b, ca,cb scalars,
    1728             :      * da = deg a, db = deg b */
    1729       27099 :     GEN c = Flx_rem(a,b, p);
    1730       27099 :     long delta = da - db;
    1731             : 
    1732       27099 :     if (both_odd(da,db)) s = -s;
    1733       27099 :     lb = Fl_mul(b[db+2], cb, p);
    1734       27099 :     a = b; b = c; dc = degpol(c);
    1735       27099 :     ind++;
    1736       27099 :     if (dc != dglist[ind]) { avma = av; return 0; } /* degenerates */
    1737       26735 :     if (g == h)
    1738             :     { /* frequent */
    1739       24621 :       ulong cc = Fl_mul(ca, Fl_powu(Fl_div(lb,g,p), delta+1, p), p);
    1740       24621 :       ca = cb;
    1741       24621 :       cb = cc;
    1742             :     }
    1743             :     else
    1744             :     {
    1745        2114 :       ulong cc = Fl_mul(ca, Fl_powu(lb, delta+1, p), p);
    1746        2114 :       ulong ghdelta = Fl_mul(g, Fl_powu(h, delta, p), p);
    1747        2114 :       ca = cb;
    1748        2114 :       cb = Fl_div(cc, ghdelta, p);
    1749             :     }
    1750       26735 :     da = db; /* = degpol(a) */
    1751       26735 :     db = dc; /* = degpol(b) */
    1752             : 
    1753       26735 :     g = lb;
    1754       26735 :     if (delta == 1)
    1755       17482 :       h = g; /* frequent */
    1756             :     else
    1757        9253 :       h = Fl_mul(h, Fl_powu(Fl_div(g,h,p), delta, p), p);
    1758             : 
    1759       26735 :     if (gc_needed(av,2))
    1760             :     {
    1761           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant_all");
    1762           0 :       gerepileall(av, 2, &a,&b);
    1763             :     }
    1764             :   }
    1765        6632 :   if (da > 1) return 0; /* Failure */
    1766             :   /* last non-constant polynomial has degree 1 */
    1767        6632 :   *C0 = Fl_mul(ca, a[2], p);
    1768        6632 :   *C1 = Fl_mul(ca, a[3], p);
    1769        6632 :   res = Fl_mul(cb, b[2], p);
    1770        6632 :   if (s == -1) res = p - res;
    1771        6632 :   avma = av; return res;
    1772             : }
    1773             : 
    1774             : /* Q a vector of polynomials representing B in Fp[X][Y], evaluate at X = x,
    1775             :  * Return 0 in case of degree drop. */
    1776             : static GEN
    1777        8557 : FlxY_evalx_drop(GEN Q, ulong x, ulong p)
    1778             : {
    1779             :   GEN z;
    1780        8557 :   long i, lb = lg(Q);
    1781        8557 :   ulong leadz = Flx_eval(leading_coeff(Q), x, p);
    1782        8557 :   long vs=mael(Q,2,1);
    1783        8557 :   if (!leadz) return zero_Flx(vs);
    1784             : 
    1785        8557 :   z = cgetg(lb, t_VECSMALL); z[1] = vs;
    1786        8557 :   for (i=2; i<lb-1; i++) z[i] = Flx_eval(gel(Q,i), x, p);
    1787        8557 :   z[i] = leadz; return z;
    1788             : }
    1789             : 
    1790             : GEN
    1791       17836 : FpXY_Fq_evaly(GEN Q, GEN y, GEN T, GEN p, long vx)
    1792             : {
    1793       17836 :   pari_sp av = avma;
    1794       17836 :   long i, lb = lg(Q);
    1795             :   GEN z;
    1796       17836 :   if (!T) return FpXY_evaly(Q, y, p, vx);
    1797        1148 :   if (lb == 2) return pol_0(vx);
    1798        1148 :   z = gel(Q, lb-1);
    1799        1148 :   if (lb == 3 || !signe(y)) return typ(z)==t_INT? scalar_ZX(z, vx): ZX_copy(z);
    1800             : 
    1801        1148 :   if (typ(z) == t_INT) z = scalar_ZX_shallow(z, vx);
    1802       28084 :   for (i=lb-2; i>=2; i--)
    1803             :   {
    1804       26936 :     GEN c = gel(Q,i);
    1805       26936 :     z = FqX_Fq_mul(z, y, T, p);
    1806       26936 :     z = typ(c) == t_INT? FqX_Fq_add(z,c,T,p): FqX_add(z,c,T,p);
    1807             :   }
    1808        1148 :   return gerepileupto(av, z);
    1809             : }
    1810             : 
    1811             : static GEN
    1812        5950 : ZX_norml1(GEN x)
    1813             : {
    1814        5950 :   long i, l = lg(x);
    1815             :   GEN s;
    1816             : 
    1817        5950 :   if (l == 2) return gen_0;
    1818        5950 :   s = gel(x, l-1); /* != 0 */
    1819       31794 :   for (i = l-2; i > 1; i--) {
    1820       25844 :     GEN xi = gel(x,i);
    1821       25844 :     if (!signe(x)) continue;
    1822       25844 :     s = addii_sign(s,1, xi,1);
    1823             :   }
    1824        5950 :   return s;
    1825             : }
    1826             : 
    1827             : /* Interpolate at roots of 1 and use Hadamard bound for univariate resultant:
    1828             :  *   bound = N_2(A)^degpol B N_2(B)^degpol(A),  where
    1829             :  *     N_2(A) = sqrt(sum (N_1(Ai))^2)
    1830             :  * Return e such that Res(A, B) < 2^e */
    1831             : ulong
    1832      582501 : ZX_ZXY_ResBound(GEN A, GEN B, GEN dB)
    1833             : {
    1834      582501 :   pari_sp av = avma;
    1835      582501 :   GEN a = gen_0, b = gen_0;
    1836      582501 :   long i , lA = lg(A), lB = lg(B);
    1837             :   double loga, logb;
    1838      582501 :   for (i=2; i<lA; i++) a = addii(a, sqri(gel(A,i)));
    1839     2079668 :   for (i=2; i<lB; i++)
    1840             :   {
    1841     1497167 :     GEN t = gel(B,i);
    1842     1497167 :     if (typ(t) == t_POL) t = ZX_norml1(t);
    1843     1497167 :     b = addii(b, sqri(t));
    1844             :   }
    1845      582501 :   loga = dbllog2(a);
    1846      582501 :   logb = dbllog2(b); if (dB) logb -= 2 * dbllog2(dB);
    1847      582501 :   i = (long)((degpol(B) * loga + degpol(A) * logb) / 2);
    1848      582501 :   avma = av; return (i <= 0)? 1: 1 + (ulong)i;
    1849             : }
    1850             : 
    1851             : /* return Res(a(Y), b(n,Y)) over Fp. la = leading_coeff(a) [for efficiency] */
    1852             : static ulong
    1853      322364 : Flx_FlxY_eval_resultant(GEN a, GEN b, ulong n, ulong p, ulong la)
    1854             : {
    1855      322364 :   GEN ev = FlxY_evalx(b, n, p);
    1856      322364 :   long drop = lg(b) - lg(ev);
    1857      322364 :   ulong r = Flx_resultant(a, ev, p);
    1858      322364 :   if (drop && la != 1) r = Fl_mul(r, Fl_powu(la, drop,p),p);
    1859      322364 :   return r;
    1860             : }
    1861             : static GEN
    1862           4 : FpX_FpXY_eval_resultant(GEN a, GEN b, GEN n, GEN p, GEN la, long db, long vX)
    1863             : {
    1864           4 :   GEN ev = FpXY_evaly(b, n, p, vX);
    1865           4 :   long drop = db-degpol(ev);
    1866           4 :   GEN r = FpX_resultant(a, ev, p);
    1867           4 :   if (drop && !gequal1(la)) r = Fp_mul(r, Fp_powu(la, drop,p),p);
    1868           4 :   return r;
    1869             : }
    1870             : 
    1871             : /* assume dres := deg(Res_X(a,b), Y) <= deg(a,X) * deg(b,Y) < p */
    1872             : /* Return a Fly */
    1873             : static GEN
    1874        7573 : Flx_FlxY_resultant_polint(GEN a, GEN b, ulong p, ulong dres, long sx)
    1875             : {
    1876        7573 :   ulong i, n, la = Flx_lead(a);
    1877        7573 :   GEN  x = cgetg(dres+2, t_VECSMALL);
    1878        7573 :   GEN  y = cgetg(dres+2, t_VECSMALL);
    1879             :  /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
    1880             :   * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
    1881      166132 :   for (i=0,n = 1; i < dres; n++)
    1882             :   {
    1883      158559 :     x[++i] = n;   y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
    1884      158559 :     x[++i] = p-n; y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
    1885             :   }
    1886        7573 :   if (i == dres)
    1887             :   {
    1888        5246 :     x[++i] = 0;   y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
    1889             :   }
    1890        7573 :   return Flv_polint(x,y, p, sx);
    1891             : }
    1892             : 
    1893             : static GEN
    1894        4679 : FlxX_pseudorem(GEN x, GEN y, ulong p)
    1895             : {
    1896        4679 :   long vx = varn(x), dx, dy, dz, i, lx, dp;
    1897        4679 :   pari_sp av = avma, av2;
    1898             : 
    1899        4679 :   if (!signe(y)) pari_err_INV("FlxX_pseudorem",y);
    1900        4679 :   (void)new_chunk(2);
    1901        4680 :   dx=degpol(x); x = RgX_recip_shallow(x)+2;
    1902        4678 :   dy=degpol(y); y = RgX_recip_shallow(y)+2; dz=dx-dy; dp = dz+1;
    1903        4677 :   av2 = avma;
    1904             :   for (;;)
    1905             :   {
    1906       33553 :     gel(x,0) = Flx_neg(gel(x,0), p); dp--;
    1907      132062 :     for (i=1; i<=dy; i++)
    1908      194816 :       gel(x,i) = Flx_add( Flx_mul(gel(y,0), gel(x,i), p),
    1909       97408 :                               Flx_mul(gel(x,0), gel(y,i), p), p );
    1910      552478 :     for (   ; i<=dx; i++)
    1911      518928 :       gel(x,i) = Flx_mul(gel(y,0), gel(x,i), p);
    1912       35453 :     do { x++; dx--; } while (dx >= 0 && lg(gel(x,0))==2);
    1913       33550 :     if (dx < dy) break;
    1914       28873 :     if (gc_needed(av2,1))
    1915             :     {
    1916           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FlxX_pseudorem dx = %ld >= %ld",dx,dy);
    1917           0 :       gerepilecoeffs(av2,x,dx+1);
    1918             :     }
    1919       28876 :   }
    1920        4677 :   if (dx < 0) return zero_Flx(0);
    1921        4677 :   lx = dx+3; x -= 2;
    1922        4677 :   x[0]=evaltyp(t_POL) | evallg(lx);
    1923        4677 :   x[1]=evalsigne(1) | evalvarn(vx);
    1924        4677 :   x = RgX_recip_shallow(x);
    1925        4682 :   if (dp)
    1926             :   { /* multiply by y[0]^dp   [beware dummy vars from FpX_FpXY_resultant] */
    1927         970 :     GEN t = Flx_powu(gel(y,0), dp, p);
    1928        3879 :     for (i=2; i<lx; i++)
    1929        2909 :       gel(x,i) = Flx_mul(gel(x,i), t, p);
    1930             :   }
    1931        4682 :   return gerepilecopy(av, x);
    1932             : }
    1933             : 
    1934             : /* return a Flx */
    1935             : GEN
    1936        1516 : FlxX_resultant(GEN u, GEN v, ulong p, long sx)
    1937             : {
    1938        1516 :   pari_sp av = avma, av2;
    1939             :   long degq,dx,dy,du,dv,dr,signh;
    1940             :   GEN z,g,h,r,p1;
    1941             : 
    1942        1516 :   dx=degpol(u); dy=degpol(v); signh=1;
    1943        1517 :   if (dx < dy)
    1944             :   {
    1945          28 :     swap(u,v); lswap(dx,dy);
    1946          28 :     if (both_odd(dx, dy)) signh = -signh;
    1947             :   }
    1948        1517 :   if (dy < 0) return zero_Flx(sx);
    1949        1517 :   if (dy==0) return gerepileupto(av, Flx_powu(gel(v,2),dx,p));
    1950             : 
    1951        1517 :   g = h = pol1_Flx(sx); av2 = avma;
    1952             :   for(;;)
    1953             :   {
    1954        4679 :     r = FlxX_pseudorem(u,v,p); dr = lg(r);
    1955        4676 :     if (dr == 2) { avma = av; return zero_Flx(sx); }
    1956        4676 :     du = degpol(u); dv = degpol(v); degq = du-dv;
    1957        4679 :     u = v; p1 = g; g = leading_coeff(u);
    1958        4678 :     switch(degq)
    1959             :     {
    1960           0 :       case 0: break;
    1961             :       case 1:
    1962        3449 :         p1 = Flx_mul(h,p1, p); h = g; break;
    1963             :       default:
    1964        1229 :         p1 = Flx_mul(Flx_powu(h,degq,p), p1, p);
    1965        1229 :         h = Flx_div(Flx_powu(g,degq,p), Flx_powu(h,degq-1,p), p);
    1966             :     }
    1967        4678 :     if (both_odd(du,dv)) signh = -signh;
    1968        4677 :     v = FlxY_Flx_div(r, p1, p);
    1969        4678 :     if (dr==3) break;
    1970        3162 :     if (gc_needed(av2,1))
    1971             :     {
    1972           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"resultant_all, dr = %ld",dr);
    1973           0 :       gerepileall(av2,4, &u, &v, &g, &h);
    1974             :     }
    1975        3163 :   }
    1976        1516 :   z = gel(v,2);
    1977        1516 :   if (dv > 1) z = Flx_div(Flx_powu(z,dv,p), Flx_powu(h,dv-1,p), p);
    1978        1516 :   if (signh < 0) z = Flx_neg(z,p);
    1979        1516 :   return gerepileupto(av, z);
    1980             : }
    1981             : 
    1982             : /* Warning:
    1983             :  * This function switches between valid and invalid variable ordering*/
    1984             : 
    1985             : static GEN
    1986        1621 : FlxY_to_FlyX(GEN b, long sv)
    1987             : {
    1988        1621 :   long i, n=-1;
    1989        1621 :   long sw = b[1]&VARNBITS;
    1990        1621 :   for(i=2;i<lg(b);i++) n = maxss(n,lgpol(gel(b,i)));
    1991        1621 :   return Flm_to_FlxX(Flm_transpose(FlxX_to_Flm(b,n)),sv,sw);
    1992             : }
    1993             : 
    1994             : /* Return a Fly*/
    1995             : GEN
    1996        1618 : Flx_FlxY_resultant(GEN a, GEN b, ulong pp)
    1997             : {
    1998        1618 :   pari_sp ltop=avma;
    1999        1618 :   long dres = degpol(a)*degpol(b);
    2000        1620 :   long sx=a[1], sy=b[1]&VARNBITS;
    2001             :   GEN z;
    2002        1620 :   b = FlxY_to_FlyX(b,sx);
    2003        1617 :   if ((ulong)dres >= pp)
    2004        1513 :     z = FlxX_resultant(Fly_to_FlxY(a, sy), b, pp, sx);
    2005             :   else
    2006         104 :     z = Flx_FlxY_resultant_polint(a, b, pp, (ulong)dres, sy);
    2007        1620 :   return gerepileupto(ltop,z);
    2008             : }
    2009             : 
    2010             : /* return a t_POL (in variable v >= 0) whose coeffs are the coeffs of b,
    2011             :  * in variable v. This is an incorrect PARI object if initially varn(b) << v.
    2012             :  * We could return a vector of coeffs, but it is convenient to have degpol()
    2013             :  * and friends available. Even in that case, it will behave nicely with all
    2014             :  * functions treating a polynomial as a vector of coeffs (eg poleval).
    2015             :  * FOR INTERNAL USE! */
    2016             : GEN
    2017        5362 : swap_vars(GEN b0, long v)
    2018             : {
    2019        5362 :   long i, n = RgX_degree(b0, v);
    2020             :   GEN b, x;
    2021        5362 :   if (n < 0) return pol_0(v);
    2022        5362 :   b = cgetg(n+3, t_POL); x = b + 2;
    2023        5362 :   b[1] = evalsigne(1) | evalvarn(v);
    2024        5362 :   for (i=0; i<=n; i++) gel(x,i) = polcoeff_i(b0, i, v);
    2025        5362 :   return b;
    2026             : }
    2027             : 
    2028             : /* assume varn(b) << varn(a) */
    2029             : /* return a FpY*/
    2030             : GEN
    2031        1587 : FpX_FpXY_resultant(GEN a, GEN b, GEN p)
    2032             : {
    2033        1587 :   long i,n,dres, db, vY = varn(b), vX = varn(a);
    2034             :   GEN la,x,y;
    2035             : 
    2036        1587 :   if (lgefint(p) == 3)
    2037             :   {
    2038        1586 :     ulong pp = uel(p,2);
    2039        1586 :     b = ZXX_to_FlxX(b, pp, vX);
    2040        1592 :     a = ZX_to_Flx(a, pp);
    2041        1592 :     x = Flx_FlxY_resultant(a, b, pp);
    2042        1594 :     return Flx_to_ZX(x);
    2043             :   }
    2044           1 :   db = RgXY_degreex(b);
    2045           1 :   dres = degpol(a)*degpol(b);
    2046           1 :   la = leading_coeff(a);
    2047           1 :   x = cgetg(dres+2, t_VEC);
    2048           1 :   y = cgetg(dres+2, t_VEC);
    2049             :  /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
    2050             :   * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
    2051           3 :   for (i=0,n = 1; i < dres; n++)
    2052             :   {
    2053           2 :     gel(x,++i) = utoipos(n);
    2054           2 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b,gel(x,i),p,la,db,vY);
    2055           2 :     gel(x,++i) = subis(p,n);
    2056           2 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b,gel(x,i),p,la,db,vY);
    2057             :   }
    2058           1 :   if (i == dres)
    2059             :   {
    2060           0 :     gel(x,++i) = gen_0;
    2061           0 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b, gel(x,i), p,la,db,vY);
    2062             :   }
    2063           1 :   return FpV_polint(x,y, p, vY);
    2064             : }
    2065             : 
    2066             : GEN
    2067         427 : FpX_direct_compositum(GEN a, GEN b, GEN p)
    2068             : {
    2069         427 :   GEN x = deg1pol_shallow(gen_1, pol_x(varn(a)), fetch_var_higher()); /* x+y */
    2070         427 :   x = FpX_FpXY_resultant(a, poleval(b,x),p);
    2071         427 :   (void)delete_var(); return x;
    2072             : }
    2073             : 
    2074             : /* 0, 1, -1, 2, -2, ... */
    2075             : #define next_lambda(a) (a>0 ? -a : 1-a)
    2076             : GEN
    2077           0 : FpX_compositum(GEN a, GEN b, GEN p)
    2078             : {
    2079           0 :   long k, v = fetch_var_higher();
    2080           0 :   for (k = 1;; k = next_lambda(k))
    2081             :   {
    2082           0 :     GEN x = deg1pol_shallow(gen_1, gmulsg(k, pol_x(v)), 0); /* x + k y */
    2083           0 :     GEN C = FpX_FpXY_resultant(a, poleval(b,x),p);
    2084           0 :     if (FpX_is_squarefree(C, p)) { (void)delete_var(); return C; }
    2085           0 :   }
    2086             : }
    2087             : 
    2088             : /* Assume A in Z[Y], B in Q[Y][X], and Res_Y(A, B) in Z[X].
    2089             :  * If lambda = NULL, return Res_Y(A,B).
    2090             :  * Otherwise, find a small lambda (start from *lambda, use the sequence above)
    2091             :  * such that R(X) = Res_Y(A(Y), B(X + lambda Y)) is squarefree, reset *lambda
    2092             :  * to the chosen value and return R
    2093             :  *
    2094             :  * If LERS is non-NULL, set it to the Last non-constant polynomial in the
    2095             :  * Euclidean Remainder Sequence */
    2096             : GEN
    2097        4172 : ZX_ZXY_resultant_all(GEN A, GEN B0, long *plambda, GEN *LERS)
    2098             : {
    2099        4172 :   int checksqfree = plambda? 1: 0, stable;
    2100        4172 :   long lambda = plambda? *plambda: 0, cnt = 0;
    2101             :   ulong bound, dp;
    2102        4172 :   pari_sp av = avma, av2 = 0;
    2103        4172 :   long i,n, lb, degA = degpol(A), dres = degA*degpol(B0);
    2104        4172 :   long v = fetch_var_higher();
    2105        4172 :   long vX = varn(B0), vY = varn(A); /* assume vY has lower priority */
    2106        4172 :   long sX = evalvarn(vX);
    2107             :   GEN x, y, dglist, dB, B, q, a, b, ev, H, H0, H1, Hp, H0p, H1p, C0, C1;
    2108             :   forprime_t S;
    2109             : 
    2110        4172 :   dglist = Hp = H0p = H1p = C0 = C1 = NULL; /* gcc -Wall */
    2111        4172 :   if (LERS)
    2112             :   {
    2113        1092 :     if (!checksqfree)
    2114           0 :       pari_err_BUG("ZX_ZXY_resultant_all [LERS != NULL needs lambda]");
    2115        1092 :     C0 = cgetg(dres+2, t_VECSMALL);
    2116        1092 :     C1 = cgetg(dres+2, t_VECSMALL);
    2117        1092 :     dglist = cgetg(dres+1, t_VECSMALL);
    2118             :   }
    2119        4172 :   x = cgetg(dres+2, t_VECSMALL);
    2120        4172 :   y = cgetg(dres+2, t_VECSMALL);
    2121        4172 :   B0 = Q_remove_denom(B0, &dB);
    2122        4172 :   if (!dB) B0 = leafcopy(B0);
    2123        4172 :   A = leafcopy(A);
    2124        4172 :   B = B0;
    2125        4172 :   setvarn(A,v);
    2126             :   /* make sure p large enough */
    2127             : INIT:
    2128             :   /* always except the first time */
    2129        5327 :   if (av2) { avma = av2; lambda = next_lambda(lambda); }
    2130        5327 :   if (lambda) B = RgX_translate(B0, monomial(stoi(lambda), 1, vY));
    2131        5327 :   B = swap_vars(B, vY); setvarn(B,v);
    2132             :   /* B0(lambda v + x, v) */
    2133        5327 :   if (DEBUGLEVEL>4 && checksqfree) err_printf("Trying lambda = %ld\n", lambda);
    2134        5327 :   av2 = avma;
    2135             : 
    2136        5327 :   if (degA <= 3)
    2137             :   { /* sub-resultant faster for small degrees */
    2138        3332 :     if (LERS)
    2139             :     { /* implies checksqfree */
    2140        1533 :       H = resultant_all(A,B,&q);
    2141        1533 :       if (typ(q) != t_POL || degpol(q)!=1) goto INIT;
    2142         945 :       H0 = gel(q,2);
    2143         945 :       if (typ(H0) == t_POL) setvarn(H0,vX); else H0 = scalarpol(H0,vX);
    2144         945 :       H1 = gel(q,3);
    2145         945 :       if (typ(H1) == t_POL) setvarn(H1,vX); else H1 = scalarpol(H1,vX);
    2146             :     }
    2147             :     else
    2148        1799 :       H = resultant(A,B);
    2149        2744 :     if (checksqfree && !ZX_is_squarefree(H)) goto INIT;
    2150        2569 :     goto END;
    2151             :   }
    2152             : 
    2153        1995 :   H = H0 = H1 = NULL;
    2154        1995 :   lb = lg(B);
    2155        1995 :   bound = ZX_ZXY_ResBound(A, B, dB);
    2156        1995 :   if (DEBUGLEVEL>4) err_printf("bound for resultant coeffs: 2^%ld\n",bound);
    2157        1995 :   dp = 1;
    2158        1995 :   init_modular_big(&S);
    2159             :   while (1)
    2160             :   {
    2161        7866 :     ulong p = u_forprime_next(&S);
    2162        7866 :     if (dB) { dp = umodiu(dB, p); if (!dp) continue; }
    2163             : 
    2164        7866 :     a = ZX_to_Flx(A, p);
    2165        7866 :     b = ZXX_to_FlxX(B, p, varn(A));
    2166        7866 :     if (LERS)
    2167             :     {
    2168             :       GEN Hi;
    2169         397 :       if (degpol(a) < degA || lg(b) < lb) continue; /* p | lc(A)lc(B) */
    2170         397 :       if (checksqfree)
    2171             :       { /* find degree list for generic Euclidean Remainder Sequence */
    2172         196 :         long goal = minss(degpol(a), degpol(b)); /* longest possible */
    2173         196 :         for (n=1; n <= goal; n++) dglist[n] = 0;
    2174         196 :         setlg(dglist, 1);
    2175        1645 :         for (n=0; n <= dres; n++)
    2176             :         {
    2177        1561 :           ev = FlxY_evalx_drop(b, n, p);
    2178        1561 :           Flx_resultant_set_dglist(a, ev, dglist, p);
    2179        1561 :           if (lg(dglist)-1 == goal) break;
    2180             :         }
    2181             :         /* last pol in ERS has degree > 1 ? */
    2182         196 :         goal = lg(dglist)-1;
    2183         196 :         if (degpol(B) == 1) { if (!goal) goto INIT; }
    2184             :         else
    2185             :         {
    2186         189 :           if (goal <= 1) goto INIT;
    2187         182 :           if (dglist[goal] != 0 || dglist[goal-1] != 1) goto INIT;
    2188             :         }
    2189         189 :         if (DEBUGLEVEL>4)
    2190           0 :           err_printf("Degree list for ERS (trials: %ld) = %Ps\n",n+1,dglist);
    2191             :       }
    2192             : 
    2193        7386 :       for (i=0,n = 0; i <= dres; n++)
    2194             :       {
    2195        6996 :         ev = FlxY_evalx_drop(b, n, p);
    2196        6996 :         x[++i] = n; y[i] = Flx_resultant_all(a, ev, C0+i, C1+i, dglist, p);
    2197        6996 :         if (!C1[i]) i--; /* C1(i) = 0. No way to recover C0(i) */
    2198             :       }
    2199         390 :       Hi = Flv_Flm_polint(x, mkvec3(y,C0,C1), p, 0);
    2200         390 :       Hp = gel(Hi,1); H0p = gel(Hi,2); H1p = gel(Hi,3);
    2201             :     }
    2202             :     else
    2203             :     {
    2204        7469 :       long dropa = degA - degpol(a), dropb = lb - lg(b);
    2205        7469 :       Hp = Flx_FlxY_resultant_polint(a, b, p, (ulong)dres, sX);
    2206        7469 :       if (dropa && dropb)
    2207           0 :         Hp = zero_Flx(sX);
    2208             :       else {
    2209        7469 :         if (dropa)
    2210             :         { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
    2211           0 :           GEN c = gel(b,lb-1); /* lc(B) */
    2212           0 :           if (!odd(lb)) c = Flx_neg(c, p); /* deg B = lb - 3 */
    2213           0 :           if (!Flx_equal1(c)) {
    2214           0 :             c = Flx_powu(c, dropa, p);
    2215           0 :             if (!Flx_equal1(c)) Hp = Flx_mul(Hp, c, p);
    2216             :           }
    2217             :         }
    2218        7469 :         else if (dropb)
    2219             :         { /* multiply by lc(A)^(deg B - deg b) */
    2220           0 :           ulong c = a[degA+2]; /* lc(A) */
    2221           0 :           c = Fl_powu(c, dropb, p);
    2222           0 :           if (c != 1) Hp = Flx_Fl_mul(Hp, c, p);
    2223             :         }
    2224             :       }
    2225             :     }
    2226        7859 :     if (!H && degpol(Hp) != dres) continue;
    2227        7859 :     if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
    2228        7859 :     if (checksqfree) {
    2229        1050 :       if (!Flx_is_squarefree(Hp, p)) goto INIT;
    2230         665 :       if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
    2231         665 :       checksqfree = 0;
    2232             :     }
    2233             : 
    2234        7474 :     if (!H)
    2235             :     { /* initialize */
    2236        1603 :       q = utoipos(p); stable = 0;
    2237        1603 :       H = ZX_init_CRT(Hp, p,vX);
    2238        1603 :       if (LERS) {
    2239         189 :         H0= ZX_init_CRT(H0p, p,vX);
    2240         189 :         H1= ZX_init_CRT(H1p, p,vX);
    2241             :       }
    2242             :     }
    2243             :     else
    2244             :     {
    2245        5871 :       if (LERS) {
    2246         201 :         GEN qp = muliu(q,p);
    2247         402 :         stable  = ZX_incremental_CRT_raw(&H, Hp, q,qp, p)
    2248         201 :                 & ZX_incremental_CRT_raw(&H0,H0p, q,qp, p)
    2249         201 :                 & ZX_incremental_CRT_raw(&H1,H1p, q,qp, p);
    2250         201 :         q = qp;
    2251             :       }
    2252             :       else
    2253        5670 :         stable = ZX_incremental_CRT(&H, Hp, &q, p);
    2254             :     }
    2255             :     /* could make it probabilistic for H ? [e.g if stable twice, etc]
    2256             :      * Probabilistic anyway for H0, H1 */
    2257        7474 :     if (DEBUGLEVEL>5 && (stable ||  ++cnt==100))
    2258           0 :     { cnt=0; err_printf("%ld%%%s ",100*expi(q)/bound,stable?"s":""); }
    2259        7474 :     if (stable && (ulong)expi(q) >= bound) break; /* DONE */
    2260        5871 :     if (gc_needed(av,2))
    2261             :     {
    2262           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_rnfequation");
    2263           0 :       gerepileall(av2, LERS? 4: 2, &H, &q, &H0, &H1);
    2264             :     }
    2265        5871 :   }
    2266             : END:
    2267        4172 :   if (DEBUGLEVEL>5) err_printf(" done\n");
    2268        4172 :   setvarn(H, vX); (void)delete_var();
    2269        4172 :   if (plambda) *plambda = lambda;
    2270        4172 :   if (LERS)
    2271             :   {
    2272        1092 :     *LERS = mkvec2(H0,H1);
    2273        1092 :     gerepileall(av, 2, &H, LERS);
    2274        1092 :     return H;
    2275             :   }
    2276        3080 :   return gerepilecopy(av, H);
    2277             : }
    2278             : 
    2279             : GEN
    2280        2478 : ZX_ZXY_rnfequation(GEN A, GEN B, long *lambda)
    2281             : {
    2282        2478 :   return ZX_ZXY_resultant_all(A, B, lambda, NULL);
    2283             : }
    2284             : 
    2285             : /* If lambda = NULL, return caract(Mod(A, T)), T,A in Z[X].
    2286             :  * Otherwise find a small lambda such that caract (Mod(A + lambda X, T)) is
    2287             :  * squarefree */
    2288             : GEN
    2289        2247 : ZXQ_charpoly_sqf(GEN A, GEN T, long *lambda, long v)
    2290             : {
    2291        2247 :   pari_sp av = avma;
    2292             :   GEN R, a;
    2293             :   long dA;
    2294             :   int delvar;
    2295             : 
    2296        2247 :   if (v < 0) v = 0;
    2297        2247 :   switch (typ(A))
    2298             :   {
    2299        2247 :     case t_POL: dA = degpol(A); if (dA > 0) break;
    2300           0 :       A = constant_coeff(A);
    2301             :     default:
    2302           0 :       if (lambda) { A = scalar_ZX_shallow(A,varn(T)); dA = 0; break;}
    2303           0 :       return gerepileupto(av, gpowgs(gsub(pol_x(v), A), degpol(T)));
    2304             :   }
    2305        2247 :   delvar = 0;
    2306        2247 :   if (varn(T) == 0)
    2307             :   {
    2308        1799 :     long v0 = fetch_var(); delvar = 1;
    2309        1799 :     T = leafcopy(T); setvarn(T,v0);
    2310        1799 :     A = leafcopy(A); setvarn(A,v0);
    2311             :   }
    2312        2247 :   R = ZX_ZXY_rnfequation(T, deg1pol_shallow(gen_1, gneg_i(A), 0), lambda);
    2313        2247 :   if (delvar) (void)delete_var();
    2314        2247 :   setvarn(R, v); a = leading_coeff(T);
    2315        2247 :   if (!gequal1(a)) R = gdiv(R, powiu(a, dA));
    2316        2247 :   return gerepileupto(av, R);
    2317             : }
    2318             : 
    2319             : 
    2320             : /* charpoly(Mod(A,T)), A may be in Q[X], but assume T and result are integral */
    2321             : GEN
    2322       56875 : ZXQ_charpoly(GEN A, GEN T, long v)
    2323             : {
    2324       56875 :   return (degpol(T) < 16) ? RgXQ_charpoly(A,T,v): ZXQ_charpoly_sqf(A,T, NULL, v);
    2325             : }
    2326             : 
    2327             : GEN
    2328       46844 : QXQ_charpoly(GEN A, GEN T, long v)
    2329             : {
    2330       46844 :   pari_sp av = avma;
    2331       46844 :   GEN den, B = Q_remove_denom(A, &den);
    2332       46844 :   GEN P = ZXQ_charpoly(B, T, v);
    2333       46844 :   return gerepilecopy(av, den ? RgX_rescale(P, ginv(den)): P);
    2334             : }
    2335             : 
    2336             : static GEN
    2337     1170229 : trivial_case(GEN A, GEN B)
    2338             : {
    2339             :   long d;
    2340     1170229 :   if (typ(A) == t_INT) return powiu(A, degpol(B));
    2341     1167764 :   d = degpol(A);
    2342     1167764 :   if (d == 0) return trivial_case(gel(A,2),B);
    2343     1166531 :   if (d < 0) return gen_0;
    2344     1166516 :   return NULL;
    2345             : }
    2346             : 
    2347             : static long
    2348      582669 : get_nbprimes(ulong bound, ulong *pt_start)
    2349             : {
    2350             : #ifdef LONG_IS_64BIT
    2351      498774 :   ulong pstart = 4611686018427388039UL;
    2352             : #else
    2353       83895 :   ulong pstart = 1073741827UL;
    2354             : #endif
    2355      582669 :   *pt_start = pstart;
    2356      582669 :   return (bound/expu(pstart))+1;
    2357             : }
    2358             : 
    2359             : static ulong
    2360      990531 : ZX_resultant_prime(GEN a, GEN b, GEN dB, long degA, long degB, ulong p)
    2361             : {
    2362      990531 :   pari_sp av = avma;
    2363             :   ulong H;
    2364             :   long dropa, dropb;
    2365      990531 :   ulong dp = dB ? umodiu(dB, p): 1;
    2366      990570 :   if (!b) b = Flx_deriv(a, p);
    2367      990595 :   dropa = degA - degpol(a);
    2368      990603 :   dropb = degB - degpol(b);
    2369      990603 :   if (dropa && dropb) /* p | lc(A), p | lc(B) */
    2370           0 :   { avma = av; return 0; }
    2371      990603 :   H = Flx_resultant(a, b, p);
    2372      990595 :   if (dropa)
    2373             :   { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
    2374           0 :     ulong c = b[degB+2]; /* lc(B) */
    2375           0 :     if (odd(degB)) c = p - c;
    2376           0 :     c = Fl_powu(c, dropa, p);
    2377           0 :     if (c != 1) H = Fl_mul(H, c, p);
    2378             :   }
    2379      990595 :   else if (dropb)
    2380             :   { /* multiply by lc(A)^(deg B - deg b) */
    2381           0 :     ulong c = a[degA+2]; /* lc(A) */
    2382           0 :     c = Fl_powu(c, dropb, p);
    2383           0 :     if (c != 1) H = Fl_mul(H, c, p);
    2384             :   }
    2385      990579 :   if (dp != 1) H = Fl_mul(H, Fl_powu(Fl_inv(dp,p), degA, p), p);
    2386      990579 :   avma = av; return H;
    2387             : }
    2388             : 
    2389             : /* If B=NULL, assume B=A' */
    2390             : static GEN
    2391      874538 : ZX_resultant_slice(GEN A, GEN B, GEN dB, GEN P, GEN *mod)
    2392             : {
    2393      874538 :   pari_sp av = avma;
    2394      874538 :   long degA, degB, i, n = lg(P)-1;
    2395             :   GEN H, T;
    2396             : 
    2397      874538 :   degA = degpol(A);
    2398      874553 :   degB = B ? degpol(B): degA - 1;
    2399      874583 :   if (n == 1)
    2400             :   {
    2401      827711 :     ulong Hp, p = uel(P,1);
    2402             :     GEN a, b;
    2403      827711 :     a = ZX_to_Flx(A, p), b = B ? ZX_to_Flx(B, p): NULL;
    2404      827689 :     Hp = ZX_resultant_prime(a, b, dB, degA, degB, p);
    2405      827714 :     avma = av;
    2406      827714 :     *mod = utoi(p); return utoi(Hp);
    2407             :   }
    2408       46872 :   T = ZV_producttree(P);
    2409       46871 :   A = ZX_nv_mod_tree(A, P, T);
    2410       46874 :   if (B) B = ZX_nv_mod_tree(B, P, T);
    2411       46874 :   H = cgetg(n+1, t_VECSMALL);
    2412      209741 :   for(i=1; i <= n; i++)
    2413             :   {
    2414      162864 :     ulong p = P[i];
    2415      162864 :     GEN a = gel(A, i), b = B ? gel(B, i): NULL;
    2416      162864 :     H[i] = ZX_resultant_prime(a, b, dB, degA, degB, p);
    2417             :   }
    2418       46877 :   H = ZV_chinese_tree(H, P, T, mod);
    2419       46878 :   gerepileall(av, 2, &H, mod);
    2420       46875 :   return H;
    2421             : }
    2422             : 
    2423             : GEN
    2424      390152 : ZX_resultant_worker(GEN P, GEN A, GEN B, GEN dB)
    2425             : {
    2426      390152 :   GEN V = cgetg(3, t_VEC);
    2427      390206 :   if (isintzero(B)) B = NULL;
    2428      390181 :   if (isintzero(dB)) dB = NULL;
    2429      390221 :   gel(V,1) = ZX_resultant_slice(A,B,dB,P,&gel(V,2));
    2430      390158 :   return V;
    2431             : }
    2432             : 
    2433             : static GEN
    2434      874658 : primelist_disc(ulong *p, long n, GEN dB)
    2435             : {
    2436      874658 :   GEN P = cgetg(n+1, t_VECSMALL);
    2437             :   long i;
    2438     1865316 :   for (i=1; i <= n; i++, *p = unextprime(*p+1))
    2439             :   {
    2440      990658 :     if (dB && umodiu(dB, *p)==0) { i--; continue; }
    2441      990658 :     P[i] = *p;
    2442             :   }
    2443      874658 :   return P;
    2444             : }
    2445             : 
    2446             : /* Res(A, B/dB), assuming the A,B in Z[X] and result is integer */
    2447             : /* if B=NULL, take B = A' */
    2448             : GEN
    2449      585149 : ZX_resultant_all(GEN A, GEN B, GEN dB, ulong bound)
    2450             : {
    2451             :   ulong p;
    2452      585149 :   pari_sp av = avma;
    2453             :   long n, m;
    2454             :   GEN  H, P, mod;
    2455      585149 :   int is_disc = !B;
    2456      585149 :   if (is_disc) B = ZX_deriv(A);
    2457             : 
    2458      585149 :   if ((H = trivial_case(A,B)) || (H = trivial_case(B,A))) return H;
    2459      582669 :   if (!bound) bound = ZX_ZXY_ResBound(A, B, dB);
    2460      582669 :   n = get_nbprimes(bound+1, &p);/* +1 to account for sign */
    2461      582669 :   if (is_disc)
    2462       25356 :     B = NULL;
    2463      582669 :   m = minss(degpol(A)+(B ? degpol(B): 0), n);
    2464      582669 :   if (m == 1)
    2465             :   {
    2466      470503 :     GEN P = primelist_disc(&p, n, dB);
    2467      470503 :     H = ZX_resultant_slice(A, B, dB, P, &mod);
    2468             :   }
    2469             :   else
    2470             :   {
    2471      112166 :     long i, s = n/m, r = n - m*s, di = 0;
    2472      112166 :     GEN worker = strtoclosure("_ZX_resultant_worker", 3, A, B?B:gen_0, dB?dB:gen_0);
    2473             :     struct pari_mt pt;
    2474             :     long pending;
    2475      112166 :     if (DEBUGLEVEL > 4)
    2476           0 :       err_printf("ZX_resultant: bound 2^%ld, nb primes: %ld\n",bound, n);
    2477      112166 :     H = cgetg(m+1+!!r, t_VEC); P = cgetg(m+1+!!r, t_VEC);
    2478      112166 :     mt_queue_start(&pt, worker);
    2479      544458 :     for (i=1; i<=m || pending; i++)
    2480             :     {
    2481             :       GEN done;
    2482      432292 :       mt_queue_submit(&pt, i, i<=m ? mkvec(primelist_disc(&p, s, dB)): NULL);
    2483      432292 :       done = mt_queue_get(&pt, NULL, &pending);
    2484      432292 :       if (done)
    2485             :       {
    2486      390252 :         di++;
    2487      390252 :         gel(H, di) = gel(done,1);
    2488      390252 :         gel(P, di) = gel(done,2);
    2489      390252 :         if (DEBUGLEVEL>5) err_printf("%ld%% ",100*di/m);
    2490             :       }
    2491             :     }
    2492      112166 :     mt_queue_end(&pt);
    2493      112166 :     if (r)
    2494             :     {
    2495       13903 :       GEN Pr = primelist_disc(&p, r, dB);
    2496       13903 :       gel(H, m+1) = ZX_resultant_slice(A, B, dB, Pr, &gel(P, m+1));
    2497             :     }
    2498      112166 :     H = ZV_chinese(H, P, &mod);
    2499      112166 :     if (DEBUGLEVEL>5) err_printf("done\n");
    2500             :   }
    2501      582669 :   H = Fp_center(H, mod, shifti(mod,-1));
    2502      582669 :   return gerepileuptoint(av, H);
    2503             : }
    2504             : 
    2505             : /* A0 and B0 in Q[X] */
    2506             : GEN
    2507       10468 : QX_resultant(GEN A0, GEN B0)
    2508             : {
    2509             :   GEN s, a, b, A, B;
    2510       10468 :   pari_sp av = avma;
    2511             : 
    2512       10468 :   A = Q_primitive_part(A0, &a);
    2513       10468 :   B = Q_primitive_part(B0, &b);
    2514       10468 :   s = ZX_resultant(A, B);
    2515       10468 :   if (!signe(s)) { avma = av; return gen_0; }
    2516       10468 :   if (a) s = gmul(s, gpowgs(a,degpol(B)));
    2517       10468 :   if (b) s = gmul(s, gpowgs(b,degpol(A)));
    2518       10468 :   return gerepileupto(av, s);
    2519             : }
    2520             : 
    2521             : GEN
    2522       23729 : ZX_resultant(GEN A, GEN B) { return ZX_resultant_all(A,B,NULL,0); }
    2523             : 
    2524             : GEN
    2525           0 : QXQ_intnorm(GEN A, GEN B)
    2526             : {
    2527             :   GEN c, n, R, lB;
    2528           0 :   long dA = degpol(A), dB = degpol(B);
    2529           0 :   pari_sp av = avma;
    2530           0 :   if (dA < 0) return gen_0;
    2531           0 :   A = Q_primitive_part(A, &c);
    2532           0 :   if (!c || typ(c) == t_INT) {
    2533           0 :     n = c;
    2534           0 :     R = ZX_resultant(B, A);
    2535             :   } else {
    2536           0 :     n = gel(c,1);
    2537           0 :     R = ZX_resultant_all(B, A, gel(c,2), 0);
    2538             :   }
    2539           0 :   if (n && !equali1(n)) R = mulii(R, powiu(n, dB));
    2540           0 :   lB = leading_coeff(B);
    2541           0 :   if (!equali1(lB)) R = diviiexact(R, powiu(lB, dA));
    2542           0 :   return gerepileuptoint(av, R);
    2543             : }
    2544             : 
    2545             : GEN
    2546           0 : QXQ_norm(GEN A, GEN B)
    2547             : {
    2548             :   GEN c, R, lB;
    2549           0 :   long dA = degpol(A), dB = degpol(B);
    2550           0 :   pari_sp av = avma;
    2551           0 :   if (dA < 0) return gen_0;
    2552           0 :   A = Q_primitive_part(A, &c);
    2553           0 :   R = ZX_resultant(B, A);
    2554           0 :   if (c) R = gmul(R, gpowgs(c, dB));
    2555           0 :   lB = leading_coeff(B);
    2556           0 :   if (!equali1(lB)) R = gdiv(R, gpowgs(lB, dA));
    2557           0 :   return gerepileupto(av, R);
    2558             : }
    2559             : 
    2560             : /* assume x has integral coefficients */
    2561             : GEN
    2562       26490 : ZX_disc_all(GEN x, ulong bound)
    2563             : {
    2564       26490 :   pari_sp av = avma;
    2565             :   GEN l, R;
    2566       26490 :   long s, d = degpol(x);
    2567       26490 :   if (d <= 1) return d ? gen_1: gen_0;
    2568       25356 :   s = (d & 2) ? -1: 1;
    2569       25356 :   l = leading_coeff(x);
    2570       25356 :   R = ZX_resultant_all(x, NULL, NULL, bound);
    2571       25356 :   if (is_pm1(l))
    2572       22570 :   { if (signe(l) < 0) s = -s; }
    2573             :   else
    2574        2786 :     R = diviiexact(R,l);
    2575       25356 :   if (s == -1) togglesign_safe(&R);
    2576       25356 :   return gerepileuptoint(av,R);
    2577             : }
    2578       24320 : GEN ZX_disc(GEN x) { return ZX_disc_all(x,0); }
    2579             : 
    2580             : GEN
    2581           0 : QX_disc(GEN x)
    2582             : {
    2583           0 :   pari_sp av = avma;
    2584           0 :   GEN c, d = ZX_disc( Q_primitive_part(x, &c) );
    2585           0 :   if (c) d = gmul(d, gpowgs(c, 2*degpol(x) - 2));
    2586           0 :   return gerepileupto(av, d);
    2587             : }
    2588             : 
    2589             : /* lift(1 / Mod(A,B)). B a ZX, A a scalar or a QX */
    2590             : GEN
    2591       10145 : QXQ_inv(GEN A, GEN B)
    2592             : {
    2593             :   GEN D, cU, q, U, V;
    2594             :   ulong p;
    2595       10145 :   pari_sp av2, av = avma;
    2596             :   forprime_t S;
    2597             :   pari_timer ti;
    2598       10145 :   if (is_scalar_t(typ(A))) return scalarpol(ginv(A), varn(B));
    2599             :   /* A a QX, B a ZX */
    2600       10145 :   if (degpol(A) < 15) return RgXQ_inv(A,B);
    2601         147 :   A = Q_primitive_part(A, &D);
    2602             :   /* A, B in Z[X] */
    2603         147 :   init_modular_small(&S);
    2604         147 :   if (DEBUGLEVEL>5) timer_start(&ti);
    2605         147 :   av2 = avma; U = NULL;
    2606        4674 :   while ((p = u_forprime_next(&S)))
    2607             :   {
    2608             :     GEN a, b, qp, Up, Vp;
    2609             :     int stable;
    2610             : 
    2611        4527 :     a = ZX_to_Flx(A, p);
    2612        4527 :     b = ZX_to_Flx(B, p);
    2613             :     /* if p | Res(A/G, B/G), discard */
    2614        4674 :     if (!Flx_extresultant(b,a,p, &Vp,&Up)) continue;
    2615             : 
    2616        4527 :     if (!U)
    2617             :     { /* First time */
    2618         147 :       U = ZX_init_CRT(Up,p,varn(A));
    2619         147 :       V = ZX_init_CRT(Vp,p,varn(A));
    2620         147 :       q = utoipos(p); continue;
    2621             :     }
    2622        4380 :     if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_inv: mod %ld (bound 2^%ld)", p,expi(q));
    2623        4380 :     qp = muliu(q,p);
    2624        8760 :     stable = ZX_incremental_CRT_raw(&U, Up, q,qp, p)
    2625        4380 :            & ZX_incremental_CRT_raw(&V, Vp, q,qp, p);
    2626        4380 :     if (stable)
    2627             :     { /* all stable: check divisibility */
    2628         147 :       GEN res = ZX_add(ZX_mul(A,U), ZX_mul(B,V));
    2629         147 :       if (degpol(res) == 0) {
    2630         147 :         res = gel(res,2);
    2631         147 :         D = D? gmul(D, res): res;
    2632         294 :         break;
    2633             :       } /* DONE */
    2634           0 :       if (DEBUGLEVEL) err_printf("QXQ_inv: char 0 check failed");
    2635             :     }
    2636        4233 :     q = qp;
    2637        4233 :     if (gc_needed(av,1))
    2638             :     {
    2639           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"QXQ_inv");
    2640           0 :       gerepileall(av2, 3, &q,&U,&V);
    2641             :     }
    2642             :   }
    2643         147 :   if (!p) pari_err_OVERFLOW("QXQ_inv [ran out of primes]");
    2644         147 :   cU = ZX_content(U);
    2645         147 :   if (!is_pm1(cU)) { U = Q_div_to_int(U, cU); D = gdiv(D, cU); }
    2646         147 :   return gerepileupto(av, RgX_Rg_div(U, D));
    2647             : }
    2648             : 
    2649             : /************************************************************************
    2650             :  *                                                                      *
    2651             :  *                   IRREDUCIBLE POLYNOMIAL / Fp                        *
    2652             :  *                                                                      *
    2653             :  ************************************************************************/
    2654             : 
    2655             : /* irreducible (unitary) polynomial of degree n over Fp */
    2656             : GEN
    2657           0 : ffinit_rand(GEN p,long n)
    2658             : {
    2659             :   for(;;) {
    2660           0 :     pari_sp av = avma;
    2661           0 :     GEN pol = ZX_add(monomial(gen_1, n, 0), random_FpX(n-1,0, p));
    2662           0 :     if (FpX_is_irred(pol, p)) return pol;
    2663           0 :     avma = av;
    2664           0 :   }
    2665             : }
    2666             : 
    2667             : /* return an extension of degree 2^l of F_2, assume l > 0
    2668             :  * Not stack clean. */
    2669             : static GEN
    2670         278 : f2init(long l)
    2671             : {
    2672             :   GEN Q, T, S;
    2673             :   long i, v;
    2674             : 
    2675         278 :   if (l == 1) return polcyclo(3, 0);
    2676         250 :   v = fetch_var_higher();
    2677         252 :   S = mkpoln(4, gen_1,gen_1,gen_0,gen_0); /* y(y^2 + y) */
    2678         250 :   Q = mkpoln(3, gen_1,gen_1, S); /* x^2 + x + y(y^2+y) */
    2679         251 :   setvarn(Q, v);
    2680             : 
    2681             :   /* x^4+x+1, irred over F_2, minimal polynomial of a root of Q */
    2682         251 :   T = mkpoln(5, gen_1,gen_0,gen_0,gen_1,gen_1);
    2683         252 :   setvarn(T, v);
    2684             :   /* Q = x^2 + x + a(y) irred. over K = F2[y] / (T(y))
    2685             :    * ==> x^2 + x + a(y) b irred. over K for any root b of Q
    2686             :    * ==> x^2 + x + (b^2+b)b */
    2687         252 :   for (i=2; i<l; i++) T = FpX_FpXY_resultant(T, Q, gen_2); /* minpoly(b) / F2*/
    2688         253 :   (void)delete_var(); setvarn(T,0); return T;
    2689             : }
    2690             : 
    2691             : /* return an extension of degree p^l of F_p, assume l > 0
    2692             :  * Not stack clean. */
    2693             : GEN
    2694           0 : ffinit_Artin_Shreier(GEN ip, long l)
    2695             : {
    2696           0 :   long i, v, p = itos(ip);
    2697           0 :   GEN T, Q, xp = monomial(gen_1,p,0); /* x^p */
    2698           0 :   T = ZX_sub(xp, deg1pol_shallow(gen_1,gen_1,0)); /* x^p - x - 1 */
    2699           0 :   if (l == 1) return T;
    2700             : 
    2701           0 :   v = fetch_var_higher();
    2702           0 :   setvarn(xp, v);
    2703           0 :   Q = ZX_sub(monomial(gen_1,2*p-1,0), monomial(gen_1,p,0));
    2704           0 :   Q = gsub(xp, deg1pol_shallow(gen_1, Q, v)); /* x^p - x - (y^(2p-1)-y^p) */
    2705           0 :   for (i = 2; i <= l; ++i) T = FpX_FpXY_resultant(T, Q, ip);
    2706           0 :   (void)delete_var(); setvarn(T,0); return T;
    2707             : }
    2708             : 
    2709             : /* check if polsubcyclo(n,l,0) is irreducible modulo p */
    2710             : static long
    2711       11432 : fpinit_check(GEN p, long n, long l)
    2712             : {
    2713             :   ulong q;
    2714       11432 :   if (!uisprime(n)) return 0;
    2715        5726 :   q = umodiu(p,n); if (!q) return 0;
    2716        5159 :   return cgcd((n-1)/Fl_order(q, n-1, n), l) == 1;
    2717             : }
    2718             : 
    2719             : /* let k=2 if p%4==1, and k=4 else and assume k*p does not divide l.
    2720             :  * Return an irreducible polynomial of degree l over F_p.
    2721             :  * Variant of Adleman and Lenstra "Finding irreducible polynomials over
    2722             :  * finite fields", ACM, 1986 (5) 350--355.
    2723             :  * Not stack clean */
    2724             : static GEN
    2725        2933 : fpinit(GEN p, long l)
    2726             : {
    2727        2933 :   ulong n = 1+l;
    2728        2933 :   while (!fpinit_check(p,n,l)) n += l;
    2729        2933 :   if (DEBUGLEVEL>=4) err_printf("FFInit: using polsubcyclo(%ld, %ld)\n",n,l);
    2730        2933 :   return FpX_red(polsubcyclo(n,l,0),p);
    2731             : }
    2732             : 
    2733             : static GEN
    2734        2691 : ffinit_fact(GEN p, long n)
    2735             : {
    2736        2691 :   GEN P, F = gel(factoru_pow(n),3);
    2737             :   long i;
    2738        2691 :   if (!odd(n) && absequaliu(p, 2))
    2739         279 :     P = f2init(vals(n)); /* if n is even, F[1] = 2^vals(n)*/
    2740             :   else
    2741        2410 :     P = fpinit(p, F[1]);
    2742        3118 :   for (i = 2; i < lg(F); ++i)
    2743         427 :     P = FpX_direct_compositum(fpinit(p, F[i]), P, p);
    2744        2691 :   return P;
    2745             : }
    2746             : 
    2747             : static GEN
    2748          96 : ffinit_nofact(GEN p, long n)
    2749             : {
    2750          96 :   GEN P, Q = NULL;
    2751          96 :   if (lgefint(p)==3)
    2752             :   {
    2753           0 :     ulong pp = p[2], q;
    2754           0 :     long v = u_lvalrem(n,pp,&q);
    2755           0 :     if (v>0)
    2756             :     {
    2757           0 :       Q = (pp == 2)? f2init(v): fpinit(p,n/q);
    2758           0 :       n = q;
    2759             :     }
    2760             :   }
    2761             :   /* n coprime to p */
    2762          96 :   if (n==1) P = Q;
    2763             :   else
    2764             :   {
    2765          96 :     P = fpinit(p, n);
    2766          96 :     if (Q) P = FpX_direct_compositum(P, Q, p);
    2767             :   }
    2768          96 :   return P;
    2769             : }
    2770             : 
    2771             : static GEN
    2772        3522 : init_Fq_i(GEN p, long n, long v)
    2773             : {
    2774             :   GEN P;
    2775        3522 :   if (n <= 0) pari_err_DOMAIN("ffinit", "degree", "<=", gen_0, stoi(n));
    2776        3522 :   if (typ(p) != t_INT) pari_err_TYPE("ffinit",p);
    2777        3522 :   if (signe(p) <= 0) pari_err_PRIME("ffinit",p);
    2778        3522 :   if (v < 0) v = 0;
    2779        3522 :   if (n == 1) return pol_x(v);
    2780        3347 :   if (fpinit_check(p, n+1, n)) return polcyclo(n+1, v);
    2781        2787 :   if (lgefint(p)-2 <= expu(n))
    2782        2691 :     P = ffinit_fact(p,n);
    2783             :   else
    2784          96 :     P = ffinit_nofact(p,n);
    2785        2787 :   setvarn(P, v); return P;
    2786             : }
    2787             : GEN
    2788        3395 : init_Fq(GEN p, long n, long v)
    2789             : {
    2790        3395 :   pari_sp av = avma;
    2791        3395 :   return gerepileupto(av, init_Fq_i(p, n, v));
    2792             : }
    2793             : GEN
    2794         126 : ffinit(GEN p, long n, long v)
    2795             : {
    2796         126 :   pari_sp av = avma;
    2797         126 :   return gerepileupto(av, FpX_to_mod(init_Fq_i(p, n, v), p));
    2798             : }
    2799             : 
    2800             : GEN
    2801        3178 : ffnbirred(GEN p, long n)
    2802             : {
    2803        3178 :   pari_sp av = avma;
    2804             :   long j;
    2805        3178 :   GEN s = gen_0, dk, pd;
    2806        3178 :   dk = divisorsu(n);
    2807       10535 :   for (j = 1; j < lg(dk); ++j)
    2808             :   {
    2809        7357 :     long d = dk[j];
    2810        7357 :     long m = moebiusu(d);
    2811        7357 :     if (!m) continue;
    2812        6797 :     pd = powiu(p, n/d);
    2813        6797 :     s = m>0 ? addii(s, pd): subii(s,pd);
    2814             :   }
    2815        3178 :   return gerepileuptoint(av, divis(s, n));
    2816             : }
    2817             : 
    2818             : GEN
    2819         427 : ffsumnbirred(GEN p, long n)
    2820             : {
    2821         427 :   pari_sp av = avma;
    2822             :   long i,j;
    2823         427 :   GEN v,q, t = gen_0;
    2824         427 :   v = cgetg(n+1,t_VECSMALL); v[1] = 1;
    2825         427 :   q = cgetg(n+1,t_VEC); gel(q,1) = p;
    2826        1470 :   for(i=2; i<=n; i++)
    2827             :   {
    2828        1043 :     v[i] = moebiusu(i);
    2829        1043 :     gel(q,i) = mulii(gel(q,i-1), p);
    2830             :   }
    2831        1897 :   for(i=1; i<=n; i++)
    2832             :   {
    2833        1470 :     GEN s = gen_0;
    2834        1470 :     GEN dk = divisorsu(i);
    2835        4445 :     for (j = 1; j < lg(dk); ++j)
    2836             :     {
    2837        2975 :       long d = dk[j], m = v[d];
    2838        2975 :       if (!m) continue;
    2839        2709 :       s = m>0 ? addii(s, gel(q, i/d)): subii(s, gel(q, i/d));
    2840             :     }
    2841        1470 :     t = addii(t, divis(s, i));
    2842             :   }
    2843         427 :   return gerepileuptoint(av, t);
    2844             : }
    2845             : 
    2846             : GEN
    2847         140 : ffnbirred0(GEN p, long n, long flag)
    2848             : {
    2849         140 :   if (typ(p) != t_INT) pari_err_TYPE("ffnbirred", p);
    2850         140 :   if (n <= 0) pari_err_DOMAIN("ffnbirred", "degree", "<=", gen_0, stoi(n));
    2851         140 :   switch(flag)
    2852             :   {
    2853             :     case 0:
    2854          70 :       return ffnbirred(p, n);
    2855             :     case 1:
    2856          70 :       return ffsumnbirred(p, n);
    2857             :     default:
    2858           0 :       pari_err_FLAG("ffnbirred");
    2859             :   }
    2860           0 :   return NULL; /* NOT REACHED */
    2861             : }

Generated by: LCOV version 1.11