Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - rootpol.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.2 lcov report (development 29115-f22e516b23) Lines: 1480 1543 95.9 %
Date: 2024-03-18 08:03:28 Functions: 116 119 97.5 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : /*******************************************************************/
      16             : /*                                                                 */
      17             : /*                ROOTS OF COMPLEX POLYNOMIALS                     */
      18             : /*  (original code contributed by Xavier Gourdon, INRIA RR 1852)   */
      19             : /*                                                                 */
      20             : /*******************************************************************/
      21             : #include "pari.h"
      22             : #include "paripriv.h"
      23             : 
      24             : #define DEBUGLEVEL DEBUGLEVEL_polroots
      25             : 
      26             : static const double pariINFINITY = 1./0.;
      27             : 
      28             : static long
      29     1148238 : isvalidcoeff(GEN x)
      30             : {
      31     1148238 :   switch (typ(x))
      32             :   {
      33     1125101 :     case t_INT: case t_REAL: case t_FRAC: return 1;
      34       23123 :     case t_COMPLEX: return isvalidcoeff(gel(x,1)) && isvalidcoeff(gel(x,2));
      35             :   }
      36          14 :   return 0;
      37             : }
      38             : 
      39             : static void
      40      247845 : checkvalidpol(GEN p, const char *f)
      41             : {
      42      247845 :   long i,n = lg(p);
      43     1349816 :   for (i=2; i<n; i++)
      44     1101978 :     if (!isvalidcoeff(gel(p,i))) pari_err_TYPE(f, gel(p,i));
      45      247838 : }
      46             : 
      47             : /********************************************************************/
      48             : /**                                                                **/
      49             : /**                   FAST ARITHMETIC over Z[i]                    **/
      50             : /**                                                                **/
      51             : /********************************************************************/
      52             : 
      53             : static GEN
      54    16189227 : ZX_to_ZiX(GEN Pr, GEN Pi)
      55             : {
      56    16189227 :   long i, lr = lg(Pr), li = lg(Pi), l = maxss(lr, li), m = minss(lr, li);
      57    16191642 :   GEN P = cgetg(l, t_POL);
      58    16196046 :   P[1] = Pr[1];
      59    65745331 :   for(i = 2; i < m; i++)
      60    49550270 :     gel(P,i) = signe(gel(Pi,i)) ? mkcomplex(gel(Pr,i), gel(Pi,i))
      61    49550270 :                                 : gel(Pr,i);
      62    22000157 :   for(     ; i < lr; i++)
      63     5805096 :     gel(P,i) = gel(Pr, i);
      64    16229995 :   for(     ; i < li; i++)
      65       34934 :     gel(P,i) = mkcomplex(gen_0, gel(Pi, i));
      66    16195061 :   return normalizepol_lg(P, l);
      67             : }
      68             : 
      69             : static GEN
      70    94677477 : ZiX_sqr(GEN P)
      71             : {
      72    94677477 :   pari_sp av = avma;
      73             :   GEN Pr2, Pi2, Qr, Qi;
      74    94677477 :   GEN Pr = real_i(P), Pi = imag_i(P);
      75    94673992 :   if (signe(Pi)==0) return gerepileupto(av, ZX_sqr(Pr));
      76    16251612 :   if (signe(Pr)==0) return gerepileupto(av, ZX_neg(ZX_sqr(Pi)));
      77    16197286 :   Pr2 = ZX_sqr(Pr); Pi2 = ZX_sqr(Pi);
      78    16188598 :   Qr = ZX_sub(Pr2, Pi2);
      79    16190025 :   if (degpol(Pr)==degpol(Pi))
      80    10531040 :     Qi = ZX_sub(ZX_sqr(ZX_add(Pr, Pi)), ZX_add(Pr2, Pi2));
      81             :   else
      82     5663303 :     Qi = ZX_shifti(ZX_mul(Pr, Pi), 1);
      83    16191661 :   return gerepilecopy(av, ZX_to_ZiX(Qr, Qi));
      84             : }
      85             : 
      86             : static GEN
      87    47328238 : graeffe(GEN p)
      88             : {
      89    47328238 :   pari_sp av = avma;
      90             :   GEN p0, p1, s0, s1;
      91    47328238 :   long n = degpol(p);
      92             : 
      93    47335336 :   if (!n) return RgX_copy(p);
      94    47335336 :   RgX_even_odd(p, &p0, &p1);
      95             :   /* p = p0(x^2) + x p1(x^2) */
      96    47340150 :   s0 = ZiX_sqr(p0);
      97    47348344 :   s1 = ZiX_sqr(p1);
      98    47347449 :   return gerepileupto(av, RgX_sub(s0, RgX_shift_shallow(s1,1)));
      99             : }
     100             : 
     101             : GEN
     102        5306 : ZX_graeffe(GEN p)
     103             : {
     104        5306 :   pari_sp av = avma;
     105             :   GEN p0, p1, s0, s1;
     106        5306 :   long n = degpol(p);
     107             : 
     108        5306 :   if (!n) return ZX_copy(p);
     109        5306 :   RgX_even_odd(p, &p0, &p1);
     110             :   /* p = p0(x^2) + x p1(x^2) */
     111        5306 :   s0 = ZX_sqr(p0);
     112        5306 :   s1 = ZX_sqr(p1);
     113        5306 :   return gerepileupto(av, ZX_sub(s0, RgX_shift_shallow(s1,1)));
     114             : }
     115             : GEN
     116          14 : polgraeffe(GEN p)
     117             : {
     118          14 :   pari_sp av = avma;
     119             :   GEN p0, p1, s0, s1;
     120          14 :   long n = degpol(p);
     121             : 
     122          14 :   if (typ(p) != t_POL) pari_err_TYPE("polgraeffe",p);
     123          14 :   n = degpol(p);
     124          14 :   if (!n) return gcopy(p);
     125          14 :   RgX_even_odd(p, &p0, &p1);
     126             :   /* p = p0(x^2) + x p1(x^2) */
     127          14 :   s0 = RgX_sqr(p0);
     128          14 :   s1 = RgX_sqr(p1);
     129          14 :   return gerepileupto(av, RgX_sub(s0, RgX_shift_shallow(s1,1)));
     130             : }
     131             : 
     132             : /********************************************************************/
     133             : /**                                                                **/
     134             : /**                       MODULUS OF ROOTS                         **/
     135             : /**                                                                **/
     136             : /********************************************************************/
     137             : 
     138             : /* Quick approximation to log2(|x|); first define y s.t. |y-x| < 2^-32 then
     139             :  * return y rounded to 2 ulp. In particular, if result < 2^21, absolute error
     140             :  * is bounded by 2^-31. If result > 2^21, it is correct to 2 ulp */
     141             : static double
     142   215843612 : mydbllog2i(GEN x)
     143             : {
     144             : #ifdef LONG_IS_64BIT
     145   186579935 :   const double W = 1/(4294967296. * 4294967296.); /* 2^-64 */
     146             : #else
     147    29263677 :   const double W = 1/4294967296.; /*2^-32*/
     148             : #endif
     149             :   GEN m;
     150   215843612 :   long lx = lgefint(x);
     151             :   double l;
     152   215843612 :   if (lx == 2) return -pariINFINITY;
     153   215053407 :   m = int_MSW(x);
     154   215053407 :   l = (double)(ulong)*m;
     155   215053407 :   if (lx == 3) return log2(l);
     156    67534308 :   l += ((double)(ulong)*int_precW(m)) * W;
     157             :   /* at least m = min(53,BIL) bits are correct in the mantissa, thus log2
     158             :    * is correct with error < log(1 + 2^-m) ~ 2^-m. Adding the correct
     159             :    * exponent BIL(lx-3) causes 1ulp further round-off error */
     160    67534308 :   return log2(l) + (double)(BITS_IN_LONG*(lx-3));
     161             : }
     162             : 
     163             : /* return log(|x|) or -pariINFINITY */
     164             : static double
     165     9199129 : mydbllogr(GEN x) {
     166     9199129 :   if (!signe(x)) return -pariINFINITY;
     167     9199129 :   return M_LN2*dbllog2r(x);
     168             : }
     169             : 
     170             : /* return log2(|x|) or -pariINFINITY */
     171             : static double
     172    52749834 : mydbllog2r(GEN x) {
     173    52749834 :   if (!signe(x)) return -pariINFINITY;
     174    52320980 :   return dbllog2r(x);
     175             : }
     176             : double
     177   290603258 : dbllog2(GEN z)
     178             : {
     179             :   double x, y;
     180   290603258 :   switch(typ(z))
     181             :   {
     182   215738275 :     case t_INT: return mydbllog2i(z);
     183       19719 :     case t_FRAC: return mydbllog2i(gel(z,1))-mydbllog2i(gel(z,2));
     184    48228210 :     case t_REAL: return mydbllog2r(z);
     185    26617054 :     default: /*t_COMPLEX*/
     186    26617054 :       x = dbllog2(gel(z,1));
     187    26699375 :       y = dbllog2(gel(z,2));
     188    26699209 :       if (x == -pariINFINITY) return y;
     189    26459841 :       if (y == -pariINFINITY) return x;
     190    26259336 :       if (fabs(x-y) > 10) return maxdd(x,y);
     191    25654760 :       return x + 0.5*log2(1 + exp2(2*(y-x)));
     192             :   }
     193             : }
     194             : static GEN /* beware overflow */
     195     6260252 : dblexp(double x) { return fabs(x) < 100.? dbltor(exp(x)): mpexp(dbltor(x)); }
     196             : 
     197             : /* find s such that  A_h <= 2^s <= 2 A_i  for one h and all i < n = deg(p),
     198             :  * with  A_i := (binom(n,i) lc(p) / p_i) ^ 1/(n-i), and  p = sum p_i X^i */
     199             : static long
     200    39135356 : findpower(GEN p)
     201             : {
     202    39135356 :   double x, L, mins = pariINFINITY;
     203    39135356 :   long n = degpol(p),i;
     204             : 
     205    39134926 :   L = dbllog2(gel(p,n+2)); /* log2(lc * binom(n,i)) */
     206   158042859 :   for (i=n-1; i>=0; i--)
     207             :   {
     208   118906731 :     L += log2((double)(i+1) / (double)(n-i));
     209   118906731 :     x = dbllog2(gel(p,i+2));
     210   118903796 :     if (x != -pariINFINITY)
     211             :     {
     212   118142069 :       double s = (L - x) / (double)(n-i);
     213   118142069 :       if (s < mins) mins = s;
     214             :     }
     215             :   }
     216    39136128 :   i = (long)ceil(mins);
     217    39136128 :   if (i - mins > 1 - 1e-12) i--;
     218    39136128 :   return i;
     219             : }
     220             : 
     221             : /* returns the exponent for logmodulus(), from the Newton diagram */
     222             : static long
     223     5275249 : newton_polygon(GEN p, long k)
     224             : {
     225     5275249 :   pari_sp av = avma;
     226     5275249 :   long n = degpol(p), i, j, h, l, *vertex = (long*)new_chunk(n+1);
     227     5275231 :   double *L = (double*)stack_malloc_align((n+1)*sizeof(double), sizeof(double));
     228             : 
     229             :   /* vertex[i] = 1 if i a vertex of convex hull, 0 otherwise */
     230    25400922 :   for (i=0; i<=n; i++) { L[i] = dbllog2(gel(p,2+i)); vertex[i] = 0; }
     231     5275201 :   vertex[0] = 1; /* sentinel */
     232    18829391 :   for (i=0; i < n; i=h)
     233             :   {
     234             :     double slope;
     235    13554190 :     h = i+1;
     236    13558907 :     while (L[i] == -pariINFINITY) { vertex[h] = 1; i = h; h = i+1; }
     237    13554190 :     slope = L[h] - L[i];
     238    36687858 :     for (j = i+2; j<=n; j++) if (L[j] != -pariINFINITY)
     239             :     {
     240    23127506 :       double pij = (L[j] - L[i])/(double)(j - i);
     241    23127506 :       if (slope < pij) { slope = pij; h = j; }
     242             :     }
     243    13554190 :     vertex[h] = 1;
     244             :   }
     245     5960092 :   h = k;   while (!vertex[h]) h++;
     246     5459519 :   l = k-1; while (!vertex[l]) l--;
     247     5275201 :   set_avma(av);
     248     5275233 :   return (long)floor((L[h]-L[l])/(double)(h-l) + 0.5);
     249             : }
     250             : 
     251             : /* change z into z*2^e, where z is real or complex of real */
     252             : static void
     253    34082815 : myshiftrc(GEN z, long e)
     254             : {
     255    34082815 :   if (typ(z)==t_COMPLEX)
     256             :   {
     257     6300385 :     if (signe(gel(z,1))) shiftr_inplace(gel(z,1), e);
     258     6300396 :     if (signe(gel(z,2))) shiftr_inplace(gel(z,2), e);
     259             :   }
     260             :   else
     261    27782430 :     if (signe(z)) shiftr_inplace(z, e);
     262    34082925 : }
     263             : 
     264             : /* return z*2^e, where z is integer or complex of integer (destroy z) */
     265             : static GEN
     266   130480526 : myshiftic(GEN z, long e)
     267             : {
     268   130480526 :   if (typ(z)==t_COMPLEX)
     269             :   {
     270    17643056 :     gel(z,1) = signe(gel(z,1))? mpshift(gel(z,1),e): gen_0;
     271    17642067 :     gel(z,2) = mpshift(gel(z,2),e);
     272    17639067 :     return z;
     273             :   }
     274   112837470 :   return signe(z)? mpshift(z,e): gen_0;
     275             : }
     276             : 
     277             : static GEN
     278     6738134 : RgX_gtofp_bit(GEN q, long bit) { return RgX_gtofp(q, nbits2prec(bit)); }
     279             : 
     280             : static GEN
     281   207550183 : mygprecrc(GEN x, long prec, long e)
     282             : {
     283             :   GEN y;
     284   207550183 :   switch(typ(x))
     285             :   {
     286   154662227 :     case t_REAL:
     287   154662227 :       if (!signe(x)) return real_0_bit(e);
     288   151517308 :       return realprec(x) == prec? x: rtor(x, prec);
     289    36243729 :     case t_COMPLEX:
     290    36243729 :       y = cgetg(3,t_COMPLEX);
     291    36243323 :       gel(y,1) = mygprecrc(gel(x,1),prec,e);
     292    36243277 :       gel(y,2) = mygprecrc(gel(x,2),prec,e);
     293    36242816 :       return y;
     294    16644227 :     default: return x;
     295             :   }
     296             : }
     297             : 
     298             : /* gprec behaves badly with the zero for polynomials.
     299             : The second parameter in mygprec is the precision in base 2 */
     300             : static GEN
     301    61051032 : mygprec(GEN x, long bit)
     302             : {
     303             :   long lx, i, e, prec;
     304             :   GEN y;
     305             : 
     306    61051032 :   if (bit < 0) bit = 0; /* should rarely happen */
     307    61051032 :   e = gexpo(x) - bit;
     308    61054155 :   prec = nbits2prec(bit);
     309    61059170 :   switch(typ(x))
     310             :   {
     311    44455625 :     case t_POL:
     312    44455625 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     313   161117211 :       for (i=2; i<lx; i++) gel(y,i) = mygprecrc(gel(x,i),prec,e);
     314    44455695 :       break;
     315             : 
     316    16603545 :     default: y = mygprecrc(x,prec,e);
     317             :   }
     318    61056077 :   return y;
     319             : }
     320             : 
     321             : /* normalize a polynomial p, that is change it with coefficients in Z[i],
     322             : after making product by 2^shift */
     323             : static GEN
     324    16829777 : pol_to_gaussint(GEN p, long shift)
     325             : {
     326    16829777 :   long i, l = lg(p);
     327    16829777 :   GEN q = cgetg(l, t_POL); q[1] = p[1];
     328    96365315 :   for (i=2; i<l; i++) gel(q,i) = gtrunc2n(gel(p,i), shift);
     329    16809497 :   return q;
     330             : }
     331             : 
     332             : /* returns a polynomial q in Z[i][x] keeping bit bits of p */
     333             : static GEN
     334    12578091 : eval_rel_pol(GEN p, long bit)
     335             : {
     336             :   long i;
     337    71112829 :   for (i = 2; i < lg(p); i++)
     338    58534683 :     if (gequal0(gel(p,i))) gel(p,i) = gen_0; /* bad behavior of gexpo */
     339    12578146 :   return pol_to_gaussint(p, bit-gexpo(p)+1);
     340             : }
     341             : 
     342             : /* returns p(R*x)/R^n (in R or R[i]), R = exp(lrho), bit bits of precision */
     343             : static GEN
     344     1537583 : homothetie(GEN p, double lrho, long bit)
     345             : {
     346             :   GEN q, r, t, iR;
     347     1537583 :   long n = degpol(p), i;
     348             : 
     349     1537583 :   iR = mygprec(dblexp(-lrho),bit);
     350     1537532 :   q = mygprec(p, bit);
     351     1537587 :   r = cgetg(n+3,t_POL); r[1] = p[1];
     352     1537583 :   t = iR; r[n+2] = q[n+2];
     353     6532254 :   for (i=n-1; i>0; i--)
     354             :   {
     355     4994714 :     gel(r,i+2) = gmul(t, gel(q,i+2));
     356     4994651 :     t = mulrr(t, iR);
     357             :   }
     358     1537540 :   gel(r,2) = gmul(t, gel(q,2)); return r;
     359             : }
     360             : 
     361             : /* change q in 2^(n*e) p(x*2^(-e)), n=deg(q)  [ ~as above with R = 2^-e ]*/
     362             : static void
     363     9526455 : homothetie2n(GEN p, long e)
     364             : {
     365     9526455 :   if (e)
     366             :   {
     367     7622758 :     long i,n = lg(p)-1;
     368    41705553 :     for (i=2; i<=n; i++) myshiftrc(gel(p,i), (n-i)*e);
     369             :   }
     370     9526621 : }
     371             : 
     372             : /* return 2^f * 2^(n*e) p(x*2^(-e)), n=deg(q) */
     373             : static void
     374    34873889 : homothetie_gauss(GEN p, long e, long f)
     375             : {
     376    34873889 :   if (e || f)
     377             :   {
     378    31169965 :     long i, n = lg(p)-1;
     379   161590673 :     for (i=2; i<=n; i++) gel(p,i) = myshiftic(gel(p,i), f+(n-i)*e);
     380             :   }
     381    34816809 : }
     382             : 
     383             : /* Lower bound on the modulus of the largest root z_0
     384             :  * k is set to an upper bound for #{z roots, |z-z_0| < eps} */
     385             : static double
     386    39134322 : lower_bound(GEN p, long *k, double eps)
     387             : {
     388    39134322 :   long n = degpol(p), i, j;
     389    39134998 :   pari_sp ltop = avma;
     390             :   GEN a, s, S, ilc;
     391             :   double r, R, rho;
     392             : 
     393    39134998 :   if (n < 4) { *k = n; return 0.; }
     394     8097406 :   S = cgetg(5,t_VEC);
     395     8098142 :   a = cgetg(5,t_VEC); ilc = gdiv(real_1(DEFAULTPREC), gel(p,n+2));
     396    40468171 :   for (i=1; i<=4; i++) gel(a,i) = gmul(ilc,gel(p,n+2-i));
     397             :   /* i = 1 split out from next loop for efficiency and initialization */
     398     8096856 :   s = gel(a,1);
     399     8096856 :   gel(S,1) = gneg(s); /* Newton sum S_i */
     400     8097541 :   rho = r = gtodouble(gabs(s,3));
     401     8098081 :   R = r / n;
     402    32385624 :   for (i=2; i<=4; i++)
     403             :   {
     404    24288013 :     s = gmulsg(i,gel(a,i));
     405    72782476 :     for (j=1; j<i; j++) s = gadd(s, gmul(gel(S,j),gel(a,i-j)));
     406    24272928 :     gel(S,i) = gneg(s); /* Newton sum S_i */
     407    24284285 :     r = gtodouble(gabs(s,3));
     408    24287543 :     if (r > 0.)
     409             :     {
     410    24241983 :       r = exp(log(r/n) / (double)i);
     411    24241983 :       if (r > R) R = r;
     412             :     }
     413             :   }
     414     8097611 :   if (R > 0. && eps < 1.2)
     415     8093885 :     *k = (long)floor((rho/R + n) / (1 + exp(-eps)*cos(eps)));
     416             :   else
     417        3726 :     *k = n;
     418     8097611 :   return gc_double(ltop, R);
     419             : }
     420             : 
     421             : /* return R such that exp(R - tau) <= rho_n(P) <= exp(R + tau)
     422             :  * P(0) != 0 and P non constant */
     423             : static double
     424     4251618 : logmax_modulus(GEN p, double tau)
     425             : {
     426             :   GEN r, q, aux, gunr;
     427     4251618 :   pari_sp av, ltop = avma;
     428     4251618 :   long i,k,n=degpol(p),nn,bit,M,e;
     429     4251605 :   double rho,eps, tau2 = (tau > 3.0)? 0.5: tau/6.;
     430             : 
     431     4251605 :   r = cgeti(BIGDEFAULTPREC);
     432     4251595 :   av = avma;
     433             : 
     434     4251595 :   eps = - 1/log(1.5*tau2); /* > 0 */
     435     4251595 :   bit = (long) ((double) n*log2(1./tau2)+3*log2((double) n))+1;
     436     4251595 :   gunr = real_1_bit(bit+2*n);
     437     4251568 :   aux = gdiv(gunr, gel(p,2+n));
     438     4251476 :   q = RgX_Rg_mul(p, aux); gel(q,2+n) = gunr;
     439     4251395 :   e = findpower(q);
     440     4251538 :   homothetie2n(q,e);
     441     4251601 :   affsi(e, r);
     442     4251603 :   q = pol_to_gaussint(q, bit);
     443     4251146 :   M = (long) (log2( log(4.*n) / (2*tau2) )) + 2;
     444     4251146 :   nn = n;
     445     4251146 :   for (i=0,e=0;;)
     446             :   { /* nn = deg(q) */
     447    39137022 :     rho = lower_bound(q, &k, eps);
     448    39134775 :     if (rho > exp2(-(double)e)) e = (long)-floor(log2(rho));
     449    39134775 :     affii(shifti(addis(r,e), 1), r);
     450    39090329 :     if (++i == M) break;
     451             : 
     452    34839537 :     bit = (long) ((double)k * log2(1./tau2) +
     453    34839537 :                      (double)(nn-k)*log2(1./eps) + 3*log2((double)nn)) + 1;
     454    34839537 :     homothetie_gauss(q, e, bit-(long)floor(dbllog2(gel(q,2+nn))+0.5));
     455    34859113 :     nn -= RgX_valrem(q, &q);
     456    34863026 :     q = gerepileupto(av, graeffe(q));
     457    34885272 :     tau2 *= 1.5; if (tau2 > 0.9) tau2 = 0.5;
     458    34885272 :     eps = -1/log(tau2); /* > 0 */
     459    34885272 :     e = findpower(q);
     460             :   }
     461     4250792 :   if (!signe(r)) return gc_double(ltop,0.);
     462     3857041 :   r = itor(r, DEFAULTPREC); shiftr_inplace(r, -M);
     463     3857557 :   return gc_double(ltop, -rtodbl(r) * M_LN2); /* -log(2) sum e_i 2^-i */
     464             : }
     465             : 
     466             : static GEN
     467       35234 : RgX_normalize1(GEN x)
     468             : {
     469       35234 :   long i, n = lg(x)-1;
     470             :   GEN y;
     471       35248 :   for (i = n; i > 1; i--)
     472       35241 :     if (!gequal0( gel(x,i) )) break;
     473       35234 :   if (i == n) return x;
     474          14 :   pari_warn(warner,"normalizing a polynomial with 0 leading term");
     475          14 :   if (i == 1) pari_err_ROOTS0("roots");
     476          14 :   y = cgetg(i+1, t_POL); y[1] = x[1];
     477          42 :   for (; i > 1; i--) gel(y,i) = gel(x,i);
     478          14 :   return y;
     479             : }
     480             : 
     481             : static GEN
     482       26978 : polrootsbound_i(GEN P, double TAU)
     483             : {
     484       26978 :   pari_sp av = avma;
     485             :   double d;
     486       26978 :   (void)RgX_valrem_inexact(P,&P);
     487       26978 :   P = RgX_normalize1(P);
     488       26978 :   switch(degpol(P))
     489             :   {
     490           7 :     case -1: pari_err_ROOTS0("roots");
     491         126 :     case 0:  set_avma(av); return gen_0;
     492             :   }
     493       26845 :   d = logmax_modulus(P, TAU) + TAU;
     494             :   /* not dblexp: result differs on ARM emulator */
     495       26845 :   return gerepileuptoleaf(av, mpexp(dbltor(d)));
     496             : }
     497             : GEN
     498       26985 : polrootsbound(GEN P, GEN tau)
     499             : {
     500       26985 :   if (typ(P) != t_POL) pari_err_TYPE("polrootsbound",P);
     501       26978 :   checkvalidpol(P, "polrootsbound");
     502       26978 :   return polrootsbound_i(P, tau? gtodouble(tau): 0.01);
     503             : }
     504             : 
     505             : /* log of modulus of the smallest root of p, with relative error tau */
     506             : static double
     507     1541547 : logmin_modulus(GEN p, double tau)
     508             : {
     509     1541547 :   pari_sp av = avma;
     510     1541547 :   if (gequal0(gel(p,2))) return -pariINFINITY;
     511     1541546 :   return gc_double(av, - logmax_modulus(RgX_recip_i(p),tau));
     512             : }
     513             : 
     514             : /* return the log of the k-th modulus (ascending order) of p, rel. error tau*/
     515             : static double
     516      581048 : logmodulus(GEN p, long k, double tau)
     517             : {
     518             :   GEN q;
     519      581048 :   long i, kk = k, imax, n = degpol(p), nn, bit, e;
     520      581048 :   pari_sp av, ltop=avma;
     521      581048 :   double r, tau2 = tau/6;
     522             : 
     523      581048 :   bit = (long)(n * (2. + log2(3.*n/tau2)));
     524      581048 :   av = avma;
     525      581048 :   q = gprec_w(p, nbits2prec(bit));
     526      581055 :   q = RgX_gtofp_bit(q, bit);
     527      581054 :   e = newton_polygon(q,k);
     528      581053 :   r = (double)e;
     529      581053 :   homothetie2n(q,e);
     530      581072 :   imax = (long)(log2(3./tau) + log2(log(4.*n)))+1;
     531     5275210 :   for (i=1; i<imax; i++)
     532             :   {
     533     4694160 :     q = eval_rel_pol(q,bit);
     534     4693818 :     kk -= RgX_valrem(q, &q);
     535     4693963 :     nn = degpol(q);
     536             : 
     537     4693960 :     q = gerepileupto(av, graeffe(q));
     538     4694196 :     e = newton_polygon(q,kk);
     539     4694191 :     r += e / exp2((double)i);
     540     4694191 :     q = RgX_gtofp_bit(q, bit);
     541     4694073 :     homothetie2n(q,e);
     542             : 
     543     4694138 :     tau2 *= 1.5; if (tau2 > 1.) tau2 = 1.;
     544     4694138 :     bit = 1 + (long)(nn*(2. + log2(3.*nn/tau2)));
     545             :   }
     546      581050 :   return gc_double(ltop, -r * M_LN2);
     547             : }
     548             : 
     549             : /* return the log of the k-th modulus r_k of p, rel. error tau, knowing that
     550             :  * rmin < r_k < rmax. This information helps because we may reduce precision
     551             :  * quicker */
     552             : static double
     553      581052 : logpre_modulus(GEN p, long k, double tau, double lrmin, double lrmax)
     554             : {
     555             :   GEN q;
     556      581052 :   long n = degpol(p), i, imax, imax2, bit;
     557      581052 :   pari_sp ltop = avma, av;
     558      581052 :   double lrho, aux, tau2 = tau/6.;
     559             : 
     560      581052 :   aux = (lrmax - lrmin) / 2. + 4*tau2;
     561      581052 :   imax = (long) log2(log((double)n)/ aux);
     562      581052 :   if (imax <= 0) return logmodulus(p,k,tau);
     563             : 
     564      572098 :   lrho  = (lrmin + lrmax) / 2;
     565      572098 :   av = avma;
     566      572098 :   bit = (long)(n*(2. + aux / M_LN2 - log2(tau2)));
     567      572098 :   q = homothetie(p, lrho, bit);
     568      572094 :   imax2 = (long)(log2(3./tau * log(4.*n))) + 1;
     569      572094 :   if (imax > imax2) imax = imax2;
     570             : 
     571     1524330 :   for (i=0; i<imax; i++)
     572             :   {
     573      952236 :     q = eval_rel_pol(q,bit);
     574      952232 :     q = gerepileupto(av, graeffe(q));
     575      952244 :     aux = 2*aux + 2*tau2;
     576      952244 :     tau2 *= 1.5;
     577      952244 :     bit = (long)(n*(2. + aux / M_LN2 - log2(1-exp(-tau2))));
     578      952244 :     q = RgX_gtofp_bit(q, bit);
     579             :   }
     580      572094 :   aux = exp2((double)imax);
     581      572094 :   return gc_double(ltop, lrho + logmodulus(q,k, aux*tau/3.) / aux);
     582             : }
     583             : 
     584             : static double
     585      858616 : ind_maxlog2(GEN q)
     586             : {
     587      858616 :   long i, k = -1;
     588      858616 :   double L = - pariINFINITY;
     589     2118884 :   for (i=0; i<=degpol(q); i++)
     590             :   {
     591     1260266 :     double d = dbllog2(gel(q,2+i));
     592     1260268 :     if (d > L) { L = d; k = i; }
     593             :   }
     594      858616 :   return k;
     595             : }
     596             : 
     597             : /* Returns k such that r_k e^(-tau) < R < r_{k+1} e^tau.
     598             :  * Assume that l <= k <= n-l */
     599             : static long
     600      965488 : dual_modulus(GEN p, double lrho, double tau, long l)
     601             : {
     602      965488 :   long i, imax, delta_k = 0, n = degpol(p), nn, v2, v, bit, ll = l;
     603      965487 :   double tau2 = tau * 7./8.;
     604      965487 :   pari_sp av = avma;
     605             :   GEN q;
     606             : 
     607      965487 :   bit = 6*n - 5*l + (long)(n*(-log2(tau2) + tau2 * 8./7.));
     608      965487 :   q = homothetie(p, lrho, bit);
     609      965473 :   imax = (long)(log(log(2.*n)/tau2)/log(7./4.)+1);
     610             : 
     611     7790609 :   for (i=0; i<imax; i++)
     612             :   {
     613     6931993 :     q = eval_rel_pol(q,bit); v2 = n - degpol(q);
     614     6930819 :     v = RgX_valrem(q, &q);
     615     6931396 :     ll -= maxss(v, v2); if (ll < 0) ll = 0;
     616             : 
     617     6931556 :     nn = degpol(q); delta_k += v;
     618     6931606 :     if (!nn) return delta_k;
     619             : 
     620     6824731 :     q = gerepileupto(av, graeffe(q));
     621     6825136 :     tau2 *= 7./4.;
     622     6825136 :     bit = 6*nn - 5*ll + (long)(nn*(-log2(tau2) + tau2 * 8./7.));
     623             :   }
     624      858616 :   return gc_long(av, delta_k + (long)ind_maxlog2(q));
     625             : }
     626             : 
     627             : /********************************************************************/
     628             : /**                                                                **/
     629             : /**              FACTORS THROUGH CIRCLE INTEGRATION                **/
     630             : /**                                                                **/
     631             : /********************************************************************/
     632             : /* l power of 2, W[step*j] = w_j; set f[j] = p(w_j)
     633             :  * if inv, w_j = exp(2IPi*j/l), else exp(-2IPi*j/l) */
     634             : 
     635             : static void
     636        7462 : fft2(GEN W, GEN p, GEN f, long step, long l)
     637             : {
     638             :   pari_sp av;
     639             :   long i, s1, l1, step2;
     640             : 
     641        7462 :   if (l == 2)
     642             :   {
     643        3766 :     gel(f,0) = gadd(gel(p,0), gel(p,step));
     644        3766 :     gel(f,1) = gsub(gel(p,0), gel(p,step)); return;
     645             :   }
     646        3696 :   av = avma;
     647        3696 :   l1 = l>>1; step2 = step<<1;
     648        3696 :   fft2(W,p,          f,   step2,l1);
     649        3696 :   fft2(W,p+step,     f+l1,step2,l1);
     650       32760 :   for (i = s1 = 0; i < l1; i++, s1 += step)
     651             :   {
     652       29064 :     GEN f0 = gel(f,i);
     653       29064 :     GEN f1 = gmul(gel(W,s1), gel(f,i+l1));
     654       29064 :     gel(f,i)    = gadd(f0, f1);
     655       29064 :     gel(f,i+l1) = gsub(f0, f1);
     656             :   }
     657        3696 :   gerepilecoeffs(av, f, l);
     658             : }
     659             : 
     660             : static void
     661    13731990 : fft(GEN W, GEN p, GEN f, long step, long l, long inv)
     662             : {
     663             :   pari_sp av;
     664             :   long i, s1, l1, l2, l3, step4;
     665             :   GEN f1, f2, f3, f02;
     666             : 
     667    13731990 :   if (l == 2)
     668             :   {
     669     6509091 :     gel(f,0) = gadd(gel(p,0), gel(p,step));
     670     6508854 :     gel(f,1) = gsub(gel(p,0), gel(p,step)); return;
     671             :   }
     672     7222899 :   av = avma;
     673     7222899 :   if (l == 4)
     674             :   {
     675             :     pari_sp av2;
     676     5067530 :     f1 = gadd(gel(p,0), gel(p,step<<1));
     677     5066946 :     f2 = gsub(gel(p,0), gel(p,step<<1));
     678     5067040 :     f3 = gadd(gel(p,step), gel(p,3*step));
     679     5066979 :     f02 = gsub(gel(p,step), gel(p,3*step));
     680     5066985 :     f02 = inv? mulcxI(f02): mulcxmI(f02);
     681     5067432 :     av2 = avma;
     682     5067432 :     gel(f,0) = gadd(f1, f3);
     683     5066800 :     gel(f,1) = gadd(f2, f02);
     684     5066960 :     gel(f,2) = gsub(f1, f3);
     685     5066856 :     gel(f,3) = gsub(f2, f02);
     686     5067083 :     gerepileallsp(av,av2,4,&gel(f,0),&gel(f,1),&gel(f,2),&gel(f,3));
     687     5067652 :     return;
     688             :   }
     689     2155369 :   l1 = l>>2; l2 = 2*l1; l3 = l1+l2; step4 = step<<2;
     690     2155369 :   fft(W,p,          f,   step4,l1,inv);
     691     2155749 :   fft(W,p+step,     f+l1,step4,l1,inv);
     692     2155781 :   fft(W,p+(step<<1),f+l2,step4,l1,inv);
     693     2155791 :   fft(W,p+3*step,   f+l3,step4,l1,inv);
     694     8022722 :   for (i = s1 = 0; i < l1; i++, s1 += step)
     695             :   {
     696     5866999 :     long s2 = s1 << 1, s3 = s1 + s2;
     697             :     GEN g02, g13, f13;
     698     5866999 :     f1 = gmul(gel(W,s1), gel(f,i+l1));
     699     5867185 :     f2 = gmul(gel(W,s2), gel(f,i+l2));
     700     5867147 :     f3 = gmul(gel(W,s3), gel(f,i+l3));
     701             : 
     702     5867224 :     f02 = gadd(gel(f,i),f2);
     703     5866575 :     g02 = gsub(gel(f,i),f2);
     704     5866654 :     f13 = gadd(f1,f3);
     705     5866577 :     g13 = gsub(f1,f3); g13 = inv? mulcxI(g13): mulcxmI(g13);
     706             : 
     707     5867074 :     gel(f,i)    = gadd(f02, f13);
     708     5866628 :     gel(f,i+l1) = gadd(g02, g13);
     709     5866700 :     gel(f,i+l2) = gsub(f02, f13);
     710     5866696 :     gel(f,i+l3) = gsub(g02, g13);
     711             :   }
     712     2155723 :   gerepilecoeffs(av, f, l);
     713             : }
     714             : 
     715             : #define code(t1,t2) ((t1 << 6) | t2)
     716             : 
     717             : static GEN
     718          98 : FFT_i(GEN W, GEN x)
     719             : {
     720          98 :   long i, l = lg(W), n = lg(x), tx = typ(x), tw, pa;
     721             :   GEN y, z, p, pol;
     722          98 :   if (l==1 || ((l-1) & (l-2))) pari_err_DIM("fft");
     723          84 :   tw = RgV_type(W, &p, &pol, &pa);
     724          84 :   if (tx == t_POL) { x++; n--; }
     725          49 :   else if (!is_vec_t(tx)) pari_err_TYPE("fft",x);
     726          84 :   if (n > l) pari_err_DIM("fft");
     727          84 :   if (n < l) {
     728           0 :     z = cgetg(l, t_VECSMALL); /* cf stackdummy */
     729           0 :     for (i = 1; i < n; i++) gel(z,i) = gel(x,i);
     730           0 :     for (     ; i < l; i++) gel(z,i) = gen_0;
     731             :   }
     732          84 :   else z = x;
     733          84 :   if (l == 2) return mkveccopy(gel(z,1));
     734          70 :   y = cgetg(l, t_VEC);
     735          70 :   if (tw==code(t_COMPLEX,t_INT) || tw==code(t_COMPLEX,t_REAL))
     736           0 :   {
     737           0 :     long inv = (l >= 5 && signe(imag_i(gel(W,1+(l>>2))))==1) ? 1 : 0;
     738           0 :     fft(W+1, z+1, y+1, 1, l-1, inv);
     739             :   } else
     740          70 :     fft2(W+1, z+1, y+1, 1, l-1);
     741          70 :   return y;
     742             : }
     743             : 
     744             : #undef code
     745             : 
     746             : GEN
     747          56 : FFT(GEN W, GEN x)
     748             : {
     749          56 :   if (!is_vec_t(typ(W))) pari_err_TYPE("fft",W);
     750          56 :   return FFT_i(W, x);
     751             : }
     752             : 
     753             : GEN
     754          56 : FFTinv(GEN W, GEN x)
     755             : {
     756          56 :   long l = lg(W), i;
     757             :   GEN w;
     758          56 :   if (!is_vec_t(typ(W))) pari_err_TYPE("fft",W);
     759          56 :   if (l==1 || ((l-1) & (l-2))) pari_err_DIM("fft");
     760          42 :   w = cgetg(l, t_VECSMALL); /* cf stackdummy */
     761          42 :   gel(w,1) = gel(W,1); /* w = gconj(W), faster */
     762        3773 :   for (i = 2; i < l; i++) gel(w, i) = gel(W, l-i+1);
     763          42 :   return FFT_i(w, x);
     764             : }
     765             : 
     766             : /* returns 1 if p has only real coefficients, 0 else */
     767             : static int
     768      916456 : isreal(GEN p)
     769             : {
     770             :   long i;
     771     4631123 :   for (i = lg(p)-1; i > 1; i--)
     772     3872514 :     if (typ(gel(p,i)) == t_COMPLEX) return 0;
     773      758609 :   return 1;
     774             : }
     775             : 
     776             : /* x non complex */
     777             : static GEN
     778      736089 : abs_update_r(GEN x, double *mu) {
     779      736089 :   GEN y = gtofp(x, DEFAULTPREC);
     780      736094 :   double ly = mydbllogr(y); if (ly < *mu) *mu = ly;
     781      736094 :   setabssign(y); return y;
     782             : }
     783             : /* return |x|, low accuracy. Set *mu = min(log(y), *mu) */
     784             : static GEN
     785     7726979 : abs_update(GEN x, double *mu) {
     786             :   GEN y, xr, yr;
     787             :   double ly;
     788     7726979 :   if (typ(x) != t_COMPLEX) return abs_update_r(x, mu);
     789     7003089 :   xr = gel(x,1);
     790     7003089 :   yr = gel(x,2);
     791     7003089 :   if (gequal0(xr)) return abs_update_r(yr,mu);
     792     7001276 :   if (gequal0(yr)) return abs_update_r(xr,mu);
     793             :   /* have to treat 0 specially: 0E-10 + 1e-20 = 0E-10 */
     794     6991143 :   xr = gtofp(xr, DEFAULTPREC);
     795     6992019 :   yr = gtofp(yr, DEFAULTPREC);
     796     6992071 :   y = sqrtr(addrr(sqrr(xr), sqrr(yr)));
     797     6991544 :   ly = mydbllogr(y); if (ly < *mu) *mu = ly;
     798     6991582 :   return y;
     799             : }
     800             : 
     801             : static void
     802      949445 : initdft(GEN *Omega, GEN *prim, long N, long Lmax, long bit)
     803             : {
     804      949445 :   long prec = nbits2prec(bit);
     805      949444 :   *Omega = grootsof1(Lmax, prec) + 1;
     806      949437 :   *prim = rootsof1u_cx(N, prec);
     807      949441 : }
     808             : 
     809             : static void
     810      470719 : parameters(GEN p, long *LMAX, double *mu, double *gamma,
     811             :            int polreal, double param, double param2)
     812             : {
     813             :   GEN q, pc, Omega, A, RU, prim, g, TWO;
     814      470719 :   long n = degpol(p), bit, NN, K, i, j, Lmax;
     815      470719 :   pari_sp av2, av = avma;
     816             : 
     817      470719 :   bit = gexpo(p) + (long)param2+8;
     818      656784 :   Lmax = 4; while (Lmax <= n) Lmax <<= 1;
     819      470725 :   NN = (long)(param*3.14)+1; if (NN < Lmax) NN = Lmax;
     820      470725 :   K = NN/Lmax; if (K & 1) K++;
     821      470725 :   NN = Lmax*K;
     822      470725 :   if (polreal) K = K/2+1;
     823             : 
     824      470725 :   initdft(&Omega, &prim, NN, Lmax, bit);
     825      470723 :   q = mygprec(p,bit) + 2;
     826      470725 :   A = cgetg(Lmax+1,t_VEC); A++;
     827      470725 :   pc= cgetg(Lmax+1,t_VEC); pc++;
     828     2830204 :   for (i=0; i <= n; i++) gel(pc,i)= gel(q,i);
     829      935700 :   for (   ; i<Lmax; i++) gel(pc,i) = gen_0;
     830             : 
     831      470724 :   *mu = pariINFINITY;
     832      470724 :   g = real_0_bit(-bit);
     833      470724 :   TWO = real2n(1, DEFAULTPREC);
     834      470729 :   av2 = avma;
     835      470729 :   RU = gen_1;
     836     1660857 :   for (i=0; i<K; i++)
     837             :   {
     838     1190134 :     if (i) {
     839      719412 :       GEN z = RU;
     840     3324319 :       for (j=1; j<n; j++)
     841             :       {
     842     2604901 :         gel(pc,j) = gmul(gel(q,j),z);
     843     2604885 :         z = gmul(z,RU); /* RU = prim^i, z=prim^(ij) */
     844             :       }
     845      719418 :       gel(pc,n) = gmul(gel(q,n),z);
     846             :     }
     847             : 
     848     1190134 :     fft(Omega,pc,A,1,Lmax,1);
     849     1190144 :     if (polreal && i>0 && i<K-1)
     850     1099795 :       for (j=0; j<Lmax; j++) g = addrr(g, divrr(TWO, abs_update(gel(A,j),mu)));
     851             :     else
     852     7817064 :       for (j=0; j<Lmax; j++) g = addrr(g, invr(abs_update(gel(A,j),mu)));
     853     1189824 :     RU = gmul(RU, prim);
     854     1190128 :     if (gc_needed(av,1))
     855             :     {
     856           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"parameters");
     857           0 :       gerepileall(av2,2, &g,&RU);
     858             :     }
     859             :   }
     860      470723 :   *gamma = mydbllog2r(divru(g,NN));
     861      470719 :   *LMAX = Lmax; set_avma(av);
     862      470719 : }
     863             : 
     864             : /* NN is a multiple of Lmax */
     865             : static void
     866      478720 : dft(GEN p, long k, long NN, long Lmax, long bit, GEN F, GEN H, long polreal)
     867             : {
     868             :   GEN Omega, q, qd, pc, pd, A, B, C, RU, aux, U, W, prim, prim2;
     869      478720 :   long n = degpol(p), i, j, K;
     870             :   pari_sp ltop;
     871             : 
     872      478720 :   initdft(&Omega, &prim, NN, Lmax, bit);
     873      478718 :   RU = cgetg(n+2,t_VEC) + 1;
     874             : 
     875      478718 :   K = NN/Lmax; if (polreal) K = K/2+1;
     876      478718 :   q = mygprec(p,bit);
     877      478717 :   qd = RgX_deriv(q);
     878             : 
     879      478713 :   A = cgetg(Lmax+1,t_VEC); A++;
     880      478714 :   B = cgetg(Lmax+1,t_VEC); B++;
     881      478715 :   C = cgetg(Lmax+1,t_VEC); C++;
     882      478715 :   pc = cgetg(Lmax+1,t_VEC); pc++;
     883      478715 :   pd = cgetg(Lmax+1,t_VEC); pd++;
     884      984320 :   gel(pc,0) = gel(q,2);  for (i=n+1; i<Lmax; i++) gel(pc,i) = gen_0;
     885     1463037 :   gel(pd,0) = gel(qd,2); for (i=n;   i<Lmax; i++) gel(pd,i) = gen_0;
     886             : 
     887      478716 :   ltop = avma;
     888      478716 :   W = cgetg(k+1,t_VEC);
     889      478716 :   U = cgetg(k+1,t_VEC);
     890     1149357 :   for (i=1; i<=k; i++) gel(W,i) = gel(U,i) = gen_0;
     891             : 
     892      478718 :   gel(RU,0) = gen_1;
     893      478718 :   prim2 = gen_1;
     894     1458654 :   for (i=0; i<K; i++)
     895             :   {
     896      979935 :     gel(RU,1) = prim2;
     897     4270146 :     for (j=1; j<n; j++) gel(RU,j+1) = gmul(gel(RU,j),prim2);
     898             :     /* RU[j] = prim^(ij)= prim2^j */
     899             : 
     900     4270133 :     for (j=1; j<n; j++) gel(pd,j) = gmul(gel(qd,j+2),gel(RU,j));
     901      979921 :     fft(Omega,pd,A,1,Lmax,1);
     902     5250055 :     for (j=1; j<=n; j++) gel(pc,j) = gmul(gel(q,j+2),gel(RU,j));
     903      979879 :     fft(Omega,pc,B,1,Lmax,1);
     904     7369996 :     for (j=0; j<Lmax; j++) gel(C,j) = ginv(gel(B,j));
     905     7369995 :     for (j=0; j<Lmax; j++) gel(B,j) = gmul(gel(A,j),gel(C,j));
     906      979850 :     fft(Omega,B,A,1,Lmax,1);
     907      979935 :     fft(Omega,C,B,1,Lmax,1);
     908             : 
     909      979936 :     if (polreal) /* p has real coefficients */
     910             :     {
     911      754171 :       if (i>0 && i<K-1)
     912             :       {
     913       99440 :         for (j=1; j<=k; j++)
     914             :         {
     915       83403 :           gel(W,j) = gadd(gel(W,j), gshift(mulreal(gel(A,j+1),gel(RU,j+1)),1));
     916       83403 :           gel(U,j) = gadd(gel(U,j), gshift(mulreal(gel(B,j),gel(RU,j)),1));
     917             :         }
     918             :       }
     919             :       else
     920             :       {
     921     1738502 :         for (j=1; j<=k; j++)
     922             :         {
     923     1000382 :           gel(W,j) = gadd(gel(W,j), mulreal(gel(A,j+1),gel(RU,j+1)));
     924     1000363 :           gel(U,j) = gadd(gel(U,j), mulreal(gel(B,j),gel(RU,j)));
     925             :         }
     926             :       }
     927             :     }
     928             :     else
     929             :     {
     930      591675 :       for (j=1; j<=k; j++)
     931             :       {
     932      365913 :         gel(W,j) = gadd(gel(W,j), gmul(gel(A,j+1),gel(RU,j+1)));
     933      365906 :         gel(U,j) = gadd(gel(U,j), gmul(gel(B,j),gel(RU,j)));
     934             :       }
     935             :     }
     936      979919 :     prim2 = gmul(prim2,prim);
     937      979933 :     gerepileall(ltop,3, &W,&U,&prim2);
     938             :   }
     939             : 
     940     1149349 :   for (i=1; i<=k; i++)
     941             :   {
     942      670638 :     aux=gel(W,i);
     943     1058312 :     for (j=1; j<i; j++) aux = gadd(aux, gmul(gel(W,i-j),gel(F,k+2-j)));
     944      670639 :     gel(F,k+2-i) = gdivgs(aux,-i*NN);
     945             :   }
     946     1149344 :   for (i=0; i<k; i++)
     947             :   {
     948      670632 :     aux=gel(U,k-i);
     949     1058305 :     for (j=1+i; j<k; j++) aux = gadd(aux,gmul(gel(F,2+j),gel(U,j-i)));
     950      670631 :     gel(H,i+2) = gdivgu(aux,NN);
     951             :   }
     952      478712 : }
     953             : 
     954             : #define NEWTON_MAX 10
     955             : static GEN
     956     2355235 : refine_H(GEN F, GEN G, GEN HH, long bit, long Sbit)
     957             : {
     958     2355235 :   GEN H = HH, D, aux;
     959     2355235 :   pari_sp ltop = avma;
     960             :   long error, i, bit1, bit2;
     961             : 
     962     2355235 :   D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F)); error = gexpo(D);
     963     2355199 :   bit2 = bit + Sbit;
     964     4307982 :   for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
     965             :   {
     966     1952785 :     if (gc_needed(ltop,1))
     967             :     {
     968           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"refine_H");
     969           0 :       gerepileall(ltop,2, &D,&H);
     970             :     }
     971     1952785 :     bit1 = -error + Sbit;
     972     1952785 :     aux = RgX_mul(mygprec(H,bit1), mygprec(D,bit1));
     973     1952777 :     aux = RgX_rem(mygprec(aux,bit1), mygprec(F,bit1));
     974             : 
     975     1952794 :     bit1 = -error*2 + Sbit; if (bit1 > bit2) bit1 = bit2;
     976     1952794 :     H = RgX_add(mygprec(H,bit1), aux);
     977     1952742 :     D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F));
     978     1952786 :     error = gexpo(D); if (error < -bit1) error = -bit1;
     979             :   }
     980     2355197 :   if (error > -bit/2) return NULL; /* FAIL */
     981     2354856 :   return gerepilecopy(ltop,H);
     982             : }
     983             : 
     984             : /* return 0 if fails, 1 else */
     985             : static long
     986      478713 : refine_F(GEN p, GEN *F, GEN *G, GEN H, long bit, double gamma)
     987             : {
     988      478713 :   GEN f0, FF, GG, r, HH = H;
     989      478713 :   long error, i, bit1 = 0, bit2, Sbit, Sbit2,  enh, normF, normG, n = degpol(p);
     990      478713 :   pari_sp av = avma;
     991             : 
     992      478713 :   FF = *F; GG = RgX_divrem(p, FF, &r);
     993      478718 :   error = gexpo(r); if (error <= -bit) error = 1-bit;
     994      478718 :   normF = gexpo(FF);
     995      478718 :   normG = gexpo(GG);
     996      478719 :   enh = gexpo(H); if (enh < 0) enh = 0;
     997      478719 :   Sbit = normF + 2*normG + enh + (long)(4.*log2((double)n)+gamma) + 1;
     998      478719 :   Sbit2 = enh + 2*(normF+normG) + (long)(2.*gamma+5.*log2((double)n)) + 1;
     999      478719 :   bit2 = bit + Sbit;
    1000     2833595 :   for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
    1001             :   {
    1002     2355224 :     if (bit1 == bit2 && i >= 2) { Sbit += n; Sbit2 += n; bit2 += n; }
    1003     2355224 :     if (gc_needed(av,1))
    1004             :     {
    1005           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"refine_F");
    1006           0 :       gerepileall(av,4, &FF,&GG,&r,&HH);
    1007             :     }
    1008             : 
    1009     2355224 :     bit1 = -error + Sbit2;
    1010     2355224 :     HH = refine_H(mygprec(FF,bit1), mygprec(GG,bit1), mygprec(HH,bit1),
    1011             :                   1-error, Sbit2);
    1012     2355238 :     if (!HH) return 0; /* FAIL */
    1013             : 
    1014     2354897 :     bit1 = -error + Sbit;
    1015     2354897 :     r = RgX_mul(mygprec(HH,bit1), mygprec(r,bit1));
    1016     2354859 :     f0 = RgX_rem(mygprec(r,bit1), mygprec(FF,bit1));
    1017             : 
    1018     2354869 :     bit1 = -2*error + Sbit; if (bit1 > bit2) bit1 = bit2;
    1019     2354869 :     FF = gadd(mygprec(FF,bit1),f0);
    1020             : 
    1021     2354859 :     bit1 = -3*error + Sbit; if (bit1 > bit2) bit1 = bit2;
    1022     2354859 :     GG = RgX_divrem(mygprec(p,bit1), mygprec(FF,bit1), &r);
    1023     2354871 :     error = gexpo(r); if (error < -bit1) error = -bit1;
    1024             :   }
    1025      478371 :   if (error>-bit) return 0; /* FAIL */
    1026      470717 :   *F = FF; *G = GG; return 1;
    1027             : }
    1028             : 
    1029             : /* returns F and G from the unit circle U such that |p-FG|<2^(-bit) |cd|,
    1030             : where cd is the leading coefficient of p */
    1031             : static void
    1032      470723 : split_fromU(GEN p, long k, double delta, long bit,
    1033             :             GEN *F, GEN *G, double param, double param2)
    1034             : {
    1035             :   GEN pp, FF, GG, H;
    1036      470723 :   long n = degpol(p), NN, bit2, Lmax;
    1037      470723 :   int polreal = isreal(p);
    1038             :   pari_sp ltop;
    1039             :   double mu, gamma;
    1040             : 
    1041      470722 :   pp = gdiv(p, gel(p,2+n));
    1042      470718 :   parameters(pp, &Lmax,&mu,&gamma, polreal,param,param2);
    1043             : 
    1044      470719 :   H  = cgetg(k+2,t_POL); H[1] = p[1];
    1045      470721 :   FF = cgetg(k+3,t_POL); FF[1]= p[1];
    1046      470724 :   gel(FF,k+2) = gen_1;
    1047             : 
    1048      470724 :   NN = (long)(0.5/delta); NN |= 1; if (NN < 2) NN = 2;
    1049      470724 :   NN *= Lmax; ltop = avma;
    1050             :   for(;;)
    1051             :   {
    1052      478719 :     bit2 = (long)(((double)NN*delta-mu)/M_LN2) + gexpo(pp) + 8;
    1053      478720 :     dft(pp, k, NN, Lmax, bit2, FF, H, polreal);
    1054      478713 :     if (refine_F(pp,&FF,&GG,H,bit,gamma)) break;
    1055        7995 :     NN <<= 1; set_avma(ltop);
    1056             :   }
    1057      470722 :   *G = gmul(GG,gel(p,2+n)); *F = FF;
    1058      470721 : }
    1059             : 
    1060             : static void
    1061      470723 : optimize_split(GEN p, long k, double delta, long bit,
    1062             :             GEN *F, GEN *G, double param, double param2)
    1063             : {
    1064      470723 :   long n = degpol(p);
    1065             :   GEN FF, GG;
    1066             : 
    1067      470723 :   if (k <= n/2)
    1068      363825 :     split_fromU(p,k,delta,bit,F,G,param,param2);
    1069             :   else
    1070             :   {
    1071      106898 :     split_fromU(RgX_recip_i(p),n-k,delta,bit,&FF,&GG,param,param2);
    1072      106897 :     *F = RgX_recip_i(GG);
    1073      106897 :     *G = RgX_recip_i(FF);
    1074             :   }
    1075      470721 : }
    1076             : 
    1077             : /********************************************************************/
    1078             : /**                                                                **/
    1079             : /**               SEARCH FOR SEPARATING CIRCLE                     **/
    1080             : /**                                                                **/
    1081             : /********************************************************************/
    1082             : 
    1083             : /* return p(2^e*x) *2^(-n*e) */
    1084             : static void
    1085           0 : scalepol2n(GEN p, long e)
    1086             : {
    1087           0 :   long i,n=lg(p)-1;
    1088           0 :   for (i=2; i<=n; i++) gel(p,i) = gmul2n(gel(p,i),(i-n)*e);
    1089           0 : }
    1090             : 
    1091             : /* returns p(x/R)*R^n; assume R is at the correct accuracy */
    1092             : static GEN
    1093     4086480 : scalepol(GEN p, GEN R, long bit)
    1094     4086480 : { return RgX_rescale(mygprec(p, bit), R); }
    1095             : 
    1096             : /* return (conj(a)X-1)^n * p[ (X-a) / (conj(a)X-1) ] */
    1097             : static GEN
    1098     1337321 : conformal_basecase(GEN p, GEN a)
    1099             : {
    1100             :   GEN z, r, ma, ca;
    1101     1337321 :   long i, n = degpol(p);
    1102             :   pari_sp av;
    1103             : 
    1104     1337319 :   if (n <= 0) return p;
    1105     1337319 :   ma = gneg(a); ca = conj_i(a);
    1106     1337323 :   av = avma;
    1107     1337323 :   z = deg1pol_shallow(ca, gen_m1, 0);
    1108     1337324 :   r = scalarpol_shallow(gel(p,2+n), 0);
    1109     3490241 :   for (i=n-1; ; i--)
    1110             :   {
    1111     3490241 :     r = RgX_addmulXn_shallow(r, gmul(ma,r), 1); /* r *= (X - a) */
    1112     3490209 :     r = gadd(r, gmul(z, gel(p,2+i)));
    1113     3490212 :     if (i == 0) return gerepileupto(av, r);
    1114     2152897 :     z = RgX_addmulXn_shallow(gmul(z,ca), gneg(z), 1); /* z *= conj(a)X - 1 */
    1115     2152922 :     if (gc_needed(av,2))
    1116             :     {
    1117           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"conformal_pol (%ld/%ld)",n-i, n);
    1118           0 :       gerepileall(av,2, &r,&z);
    1119             :     }
    1120             :   }
    1121             : }
    1122             : static GEN
    1123     1337441 : conformal_pol(GEN p, GEN a)
    1124             : {
    1125     1337441 :   pari_sp av = avma;
    1126     1337441 :   long d, nR, n = degpol(p), v;
    1127             :   GEN Q, R, S, T;
    1128     1337442 :   if (n < 35) return conformal_basecase(p, a);
    1129         119 :   d = (n+1) >> 1; v = varn(p);
    1130         119 :   Q = RgX_shift_shallow(p, -d);
    1131         119 :   R = RgXn_red_shallow(p, d);
    1132         119 :   Q = conformal_pol(Q, a);
    1133         119 :   R = conformal_pol(R, a);
    1134         119 :   S = gpowgs(deg1pol_shallow(gen_1, gneg(a), v), d);
    1135         119 :   T = RgX_recip_i(S);
    1136         119 :   if (typ(a) == t_COMPLEX) T = gconj(T);
    1137         119 :   if (odd(d)) T = RgX_neg(T);
    1138             :   /* S = (X - a)^d, T = (conj(a) X - 1)^d */
    1139         119 :   nR = n - degpol(R) - d; /* >= 0 */
    1140         119 :   if (nR) T = RgX_mul(T, gpowgs(deg1pol_shallow(gconj(a), gen_m1, v), nR));
    1141         119 :   return gerepileupto(av, RgX_add(RgX_mul(Q, S), RgX_mul(R, T)));
    1142             : }
    1143             : 
    1144             : static const double UNDEF = -100000.;
    1145             : 
    1146             : static double
    1147      470719 : logradius(double *radii, GEN p, long k, double aux, double *delta)
    1148             : {
    1149      470719 :   long i, n = degpol(p);
    1150             :   double lrho, lrmin, lrmax;
    1151      470719 :   if (k > 1)
    1152             :   {
    1153      273949 :     i = k-1; while (i>0 && radii[i] == UNDEF) i--;
    1154      200910 :     lrmin = logpre_modulus(p,k,aux, radii[i], radii[k]);
    1155             :   }
    1156             :   else /* k=1 */
    1157      269809 :     lrmin = logmin_modulus(p,aux);
    1158      470723 :   radii[k] = lrmin;
    1159             : 
    1160      470723 :   if (k+1<n)
    1161             :   {
    1162      568654 :     i = k+2; while (i<=n && radii[i] == UNDEF) i++;
    1163      380142 :     lrmax = logpre_modulus(p,k+1,aux, radii[k+1], radii[i]);
    1164             :   }
    1165             :   else /* k+1=n */
    1166       90581 :     lrmax = logmax_modulus(p,aux);
    1167      470718 :   radii[k+1] = lrmax;
    1168             : 
    1169      470718 :   lrho = radii[k];
    1170      791968 :   for (i=k-1; i>=1; i--)
    1171             :   {
    1172      321250 :     if (radii[i] == UNDEF || radii[i] > lrho)
    1173      235064 :       radii[i] = lrho;
    1174             :     else
    1175       86186 :       lrho = radii[i];
    1176             :   }
    1177      470718 :   lrho = radii[k+1];
    1178     1567491 :   for (i=k+1; i<=n; i++)
    1179             :   {
    1180     1096773 :     if (radii[i] == UNDEF || radii[i] < lrho)
    1181      545719 :       radii[i] = lrho;
    1182             :     else
    1183      551054 :       lrho = radii[i];
    1184             :   }
    1185      470718 :   *delta = (lrmax - lrmin) / 2;
    1186      470718 :   if (*delta > 1.) *delta = 1.;
    1187      470718 :   return (lrmin + lrmax) / 2;
    1188             : }
    1189             : 
    1190             : static void
    1191      470720 : update_radius(long n, double *radii, double lrho, double *par, double *par2)
    1192             : {
    1193      470720 :   double t, param = 0., param2 = 0.;
    1194             :   long i;
    1195     2359387 :   for (i=1; i<=n; i++)
    1196             :   {
    1197     1888699 :     radii[i] -= lrho;
    1198     1888699 :     t = fabs(rtodbl( invr(subsr(1, dblexp(radii[i]))) ));
    1199     1888667 :     param += t; if (t > 1.) param2 += log2(t);
    1200             :   }
    1201      470688 :   *par = param; *par2 = param2;
    1202      470688 : }
    1203             : 
    1204             : /* apply the conformal mapping then split from U */
    1205             : static void
    1206      445736 : conformal_mapping(double *radii, GEN ctr, GEN p, long k, long bit,
    1207             :                   double aux, GEN *F,GEN *G)
    1208             : {
    1209      445736 :   long bit2, n = degpol(p), i;
    1210      445736 :   pari_sp ltop = avma, av;
    1211             :   GEN q, FF, GG, a, R;
    1212             :   double lrho, delta, param, param2;
    1213             :   /* n * (2.*log2(2.732)+log2(1.5)) + 1 */
    1214      445736 :   bit2 = bit + (long)(n*3.4848775) + 1;
    1215      445736 :   a = sqrtr_abs( utor(3, 2*MEDDEFAULTPREC - 2) );
    1216      445738 :   a = divrs(a, -6);
    1217      445738 :   a = gmul(mygprec(a,bit2), mygprec(ctr,bit2)); /* a = -ctr/2sqrt(3) */
    1218             : 
    1219      445737 :   av = avma;
    1220      445737 :   q = conformal_pol(mygprec(p,bit2), a);
    1221     2190902 :   for (i=1; i<=n; i++)
    1222     1745170 :     if (radii[i] != UNDEF) /* update array radii */
    1223             :     {
    1224     1472040 :       pari_sp av2 = avma;
    1225     1472040 :       GEN t, r = dblexp(radii[i]), r2 = sqrr(r);
    1226             :       /* 2(r^2 - 1) / (r^2 - 3(r-1)) */
    1227     1471987 :       t = divrr(shiftr((subrs(r2,1)),1), subrr(r2, mulur(3,subrs(r,1))));
    1228     1472030 :       radii[i] = mydbllogr(addsr(1,t)) / 2;
    1229     1472016 :       set_avma(av2);
    1230             :     }
    1231      445732 :   lrho = logradius(radii, q,k,aux/10., &delta);
    1232      445734 :   update_radius(n, radii, lrho, &param, &param2);
    1233             : 
    1234      445733 :   bit2 += (long)(n * fabs(lrho)/M_LN2 + 1.);
    1235      445733 :   R = mygprec(dblexp(-lrho), bit2);
    1236      445737 :   q = scalepol(q,R,bit2);
    1237      445737 :   gerepileall(av,2, &q,&R);
    1238             : 
    1239      445737 :   optimize_split(q,k,delta,bit2,&FF,&GG,param,param2);
    1240      445736 :   bit2 += n; R = invr(R);
    1241      445733 :   FF = scalepol(FF,R,bit2);
    1242      445736 :   GG = scalepol(GG,R,bit2);
    1243             : 
    1244      445735 :   a = mygprec(a,bit2);
    1245      445736 :   FF = conformal_pol(FF,a);
    1246      445734 :   GG = conformal_pol(GG,a);
    1247             : 
    1248      445739 :   a = invr(subsr(1, gnorm(a)));
    1249      445736 :   FF = RgX_Rg_mul(FF, powru(a,k));
    1250      445735 :   GG = RgX_Rg_mul(GG, powru(a,n-k));
    1251             : 
    1252      445734 :   *F = mygprec(FF,bit+n);
    1253      445737 :   *G = mygprec(GG,bit+n); gerepileall(ltop,2, F,G);
    1254      445738 : }
    1255             : 
    1256             : /* split p, this time without scaling. returns in F and G two polynomials
    1257             :  * such that |p-FG|< 2^(-bit)|p| */
    1258             : static void
    1259      470722 : split_2(GEN p, long bit, GEN ctr, double thickness, GEN *F, GEN *G)
    1260             : {
    1261             :   GEN q, FF, GG, R;
    1262             :   double aux, delta, param, param2;
    1263      470722 :   long n = degpol(p), i, j, k, bit2;
    1264             :   double lrmin, lrmax, lrho, *radii;
    1265             : 
    1266      470722 :   radii = (double*) stack_malloc_align((n+1) * sizeof(double), sizeof(double));
    1267             : 
    1268     1418019 :   for (i=2; i<n; i++) radii[i] = UNDEF;
    1269      470722 :   aux = thickness/(double)(4 * n);
    1270      470722 :   lrmin = logmin_modulus(p, aux);
    1271      470720 :   lrmax = logmax_modulus(p, aux);
    1272      470719 :   radii[1] = lrmin;
    1273      470719 :   radii[n] = lrmax;
    1274      470719 :   i = 1; j = n;
    1275      470719 :   lrho = (lrmin + lrmax) / 2;
    1276      470719 :   k = dual_modulus(p, lrho, aux, 1);
    1277      470725 :   if (5*k < n || (n < 2*k && 5*k < 4*n))
    1278       74452 :     { lrmax = lrho; j=k+1; radii[j] = lrho; }
    1279             :   else
    1280      396273 :     { lrmin = lrho; i=k;   radii[i] = lrho; }
    1281      965494 :   while (j > i+1)
    1282             :   {
    1283      494771 :     if (i+j == n+1)
    1284      351329 :       lrho = (lrmin + lrmax) / 2;
    1285             :     else
    1286             :     {
    1287      143442 :       double kappa = 2. - log(1. + minss(i,n-j)) / log(1. + minss(j,n-i));
    1288      143442 :       if (i+j < n+1) lrho = lrmax * kappa + lrmin;
    1289      114497 :       else           lrho = lrmin * kappa + lrmax;
    1290      143442 :       lrho /= 1+kappa;
    1291             :     }
    1292      494771 :     aux = (lrmax - lrmin) / (4*(j-i));
    1293      494771 :     k = dual_modulus(p, lrho, aux, minss(i,n+1-j));
    1294      494769 :     if (k-i < j-k-1 || (k-i == j-k-1 && 2*k > n))
    1295      367802 :       { lrmax = lrho; j=k+1; radii[j] = lrho - aux; }
    1296             :     else
    1297      126967 :       { lrmin = lrho; i=k;   radii[i] = lrho + aux; }
    1298             :   }
    1299      470723 :   aux = lrmax - lrmin;
    1300             : 
    1301      470723 :   if (ctr)
    1302             :   {
    1303      445736 :     lrho = (lrmax + lrmin) / 2;
    1304     2190931 :     for (i=1; i<=n; i++)
    1305     1745195 :       if (radii[i] != UNDEF) radii[i] -= lrho;
    1306             : 
    1307      445736 :     bit2 = bit + (long)(n * fabs(lrho)/M_LN2 + 1.);
    1308      445736 :     R = mygprec(dblexp(-lrho), bit2);
    1309      445734 :     q = scalepol(p,R,bit2);
    1310      445737 :     conformal_mapping(radii, ctr, q, k, bit2, aux, &FF, &GG);
    1311             :   }
    1312             :   else
    1313             :   {
    1314       24987 :     lrho = logradius(radii, p, k, aux/10., &delta);
    1315       24986 :     update_radius(n, radii, lrho, &param, &param2);
    1316             : 
    1317       24985 :     bit2 = bit + (long)(n * fabs(lrho)/M_LN2 + 1.);
    1318       24985 :     R = mygprec(dblexp(-lrho), bit2);
    1319       24986 :     q = scalepol(p,R,bit2);
    1320       24986 :     optimize_split(q, k, delta, bit2, &FF, &GG, param, param2);
    1321             :   }
    1322      470723 :   bit  += n;
    1323      470723 :   bit2 += n; R = invr(mygprec(R,bit2));
    1324      470718 :   *F = mygprec(scalepol(FF,R,bit2), bit);
    1325      470723 :   *G = mygprec(scalepol(GG,R,bit2), bit);
    1326      470724 : }
    1327             : 
    1328             : /* procedure corresponding to steps 5,6,.. page 44 in RR n. 1852 */
    1329             : /* put in F and G two polynomial such that |p-FG|<2^(-bit)|p|
    1330             :  * where the maximum modulus of the roots of p is <=1.
    1331             :  * Assume sum of roots is 0. */
    1332             : static void
    1333      445736 : split_1(GEN p, long bit, GEN *F, GEN *G)
    1334             : {
    1335      445736 :   long i, imax, n = degpol(p), polreal = isreal(p), ep = gexpo(p), bit2 = bit+n;
    1336             :   GEN ctr, q, qq, FF, GG, v, gr, r, newq;
    1337             :   double lrmin, lrmax, lthick;
    1338      445736 :   const double LOG3 = 1.098613;
    1339             : 
    1340      445736 :   lrmax = logmax_modulus(p, 0.01);
    1341      445736 :   gr = mygprec(dblexp(-lrmax), bit2);
    1342      445737 :   q = scalepol(p,gr,bit2);
    1343             : 
    1344      445739 :   bit2 = bit + gexpo(q) - ep + (long)((double)n*2.*log2(3.)+1);
    1345      445738 :   v = cgetg(5,t_VEC);
    1346      445738 :   gel(v,1) = gen_2;
    1347      445738 :   gel(v,2) = gen_m2;
    1348      445738 :   gel(v,3) = mkcomplex(gen_0, gel(v,1));
    1349      445737 :   gel(v,4) = mkcomplex(gen_0, gel(v,2));
    1350      445737 :   q = mygprec(q,bit2); lthick = 0;
    1351      445738 :   newq = ctr = NULL; /* -Wall */
    1352      445738 :   imax = polreal? 3: 4;
    1353      807505 :   for (i=1; i<=imax; i++)
    1354             :   {
    1355      801009 :     qq = RgX_translate(q, gel(v,i));
    1356      801020 :     lrmin = logmin_modulus(qq,0.05);
    1357      801008 :     if (LOG3 > lrmin + lthick)
    1358             :     {
    1359      789203 :       double lquo = logmax_modulus(qq,0.05) - lrmin;
    1360      789206 :       if (lquo > lthick) { lthick = lquo; newq = qq; ctr = gel(v,i); }
    1361             :     }
    1362      801011 :     if (lthick > M_LN2) break;
    1363      411505 :     if (polreal && i==2 && lthick > LOG3 - M_LN2) break;
    1364             :   }
    1365      445740 :   bit2 = bit + gexpo(newq) - ep + (long)(n*LOG3/M_LN2 + 1);
    1366      445737 :   split_2(newq, bit2, ctr, lthick, &FF, &GG);
    1367      445738 :   r = gneg(mygprec(ctr,bit2));
    1368      445738 :   FF = RgX_translate(FF,r);
    1369      445737 :   GG = RgX_translate(GG,r);
    1370             : 
    1371      445739 :   gr = invr(gr); bit2 = bit - ep + gexpo(FF)+gexpo(GG);
    1372      445739 :   *F = scalepol(FF,gr,bit2);
    1373      445739 :   *G = scalepol(GG,gr,bit2);
    1374      445738 : }
    1375             : 
    1376             : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|,
    1377             : where the maximum modulus of the roots of p is < 0.5 */
    1378             : static int
    1379      446061 : split_0_2(GEN p, long bit, GEN *F, GEN *G)
    1380             : {
    1381             :   GEN q, b;
    1382      446061 :   long n = degpol(p), k, bit2, eq;
    1383      446061 :   double aux0 = dbllog2(gel(p,n+2)); /* != -oo */
    1384      446061 :   double aux1 = dbllog2(gel(p,n+1)), aux;
    1385             : 
    1386      446062 :   if (aux1 == -pariINFINITY) /* p1 = 0 */
    1387        9621 :     aux = 0;
    1388             :   else
    1389             :   {
    1390      436441 :     aux = aux1 - aux0; /* log2(p1/p0) */
    1391             :     /* beware double overflow */
    1392      436441 :     if (aux >= 0 && (aux > 1e4 || exp2(aux) > 2.5*n)) return 0;
    1393      436441 :     aux = (aux < -300)? 0.: n*log2(1 + exp2(aux)/(double)n);
    1394             :   }
    1395      446062 :   bit2 = bit+1 + (long)(log2((double)n) + aux);
    1396      446062 :   q = mygprec(p,bit2);
    1397      446063 :   if (aux1 == -pariINFINITY) b = NULL;
    1398             :   else
    1399             :   {
    1400      436442 :     b = gdivgs(gdiv(gel(q,n+1),gel(q,n+2)),-n);
    1401      436441 :     q = RgX_translate(q,b);
    1402             :   }
    1403      446062 :   gel(q,n+1) = gen_0; eq = gexpo(q);
    1404      446060 :   k = 0;
    1405      446614 :   while (k <= n/2 && (- gexpo(gel(q,k+2)) > bit2 + 2*(n-k) + eq
    1406      446486 :                       || gequal0(gel(q,k+2)))) k++;
    1407      446060 :   if (k > 0)
    1408             :   {
    1409         324 :     if (k > n/2) k = n/2;
    1410         324 :     bit2 += k<<1;
    1411         324 :     *F = pol_xn(k, 0);
    1412         324 :     *G = RgX_shift_shallow(q, -k);
    1413             :   }
    1414             :   else
    1415             :   {
    1416      445736 :     split_1(q,bit2,F,G);
    1417      445738 :     bit2 = bit + gexpo(*F) + gexpo(*G) - gexpo(p) + (long)aux+1;
    1418      445738 :     *F = mygprec(*F,bit2);
    1419             :   }
    1420      446062 :   *G = mygprec(*G,bit2);
    1421      446062 :   if (b)
    1422             :   {
    1423      436441 :     GEN mb = mygprec(gneg(b), bit2);
    1424      436441 :     *F = RgX_translate(*F, mb);
    1425      436442 :     *G = RgX_translate(*G, mb);
    1426             :   }
    1427      446062 :   return 1;
    1428             : }
    1429             : 
    1430             : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|.
    1431             :  * Assume max_modulus(p) < 2 */
    1432             : static void
    1433      446060 : split_0_1(GEN p, long bit, GEN *F, GEN *G)
    1434             : {
    1435             :   GEN FF, GG;
    1436             :   long n, bit2, normp;
    1437             : 
    1438      446060 :   if  (split_0_2(p,bit,F,G)) return;
    1439             : 
    1440           0 :   normp = gexpo(p);
    1441           0 :   scalepol2n(p,2); /* p := 4^(-n) p(4*x) */
    1442           0 :   n = degpol(p); bit2 = bit + 2*n + gexpo(p) - normp;
    1443           0 :   split_1(mygprec(p,bit2), bit2,&FF,&GG);
    1444           0 :   scalepol2n(FF,-2);
    1445           0 :   scalepol2n(GG,-2); bit2 = bit + gexpo(FF) + gexpo(GG) - normp;
    1446           0 :   *F = mygprec(FF,bit2);
    1447           0 :   *G = mygprec(GG,bit2);
    1448             : }
    1449             : 
    1450             : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P| */
    1451             : static void
    1452      471048 : split_0(GEN p, long bit, GEN *F, GEN *G)
    1453             : {
    1454      471048 :   const double LOG1_9 = 0.6418539;
    1455      471048 :   long n = degpol(p), k = 0;
    1456             :   GEN q;
    1457             : 
    1458      471048 :   while (gexpo(gel(p,k+2)) < -bit && k <= n/2) k++;
    1459      471048 :   if (k > 0)
    1460             :   {
    1461           0 :     if (k > n/2) k = n/2;
    1462           0 :     *F = pol_xn(k, 0);
    1463           0 :     *G = RgX_shift_shallow(p, -k);
    1464             :   }
    1465             :   else
    1466             :   {
    1467      471048 :     double lr = logmax_modulus(p, 0.05);
    1468      471042 :     if (lr < LOG1_9) split_0_1(p, bit, F, G);
    1469             :     else
    1470             :     {
    1471      416001 :       q = RgX_recip_i(p);
    1472      416006 :       lr = logmax_modulus(q,0.05);
    1473      416005 :       if (lr < LOG1_9)
    1474             :       {
    1475      391019 :         split_0_1(q, bit, F, G);
    1476      391022 :         *F = RgX_recip_i(*F);
    1477      391022 :         *G = RgX_recip_i(*G);
    1478             :       }
    1479             :       else
    1480       24986 :         split_2(p,bit,NULL, 1.2837,F,G);
    1481             :     }
    1482             :   }
    1483      471049 : }
    1484             : 
    1485             : /********************************************************************/
    1486             : /**                                                                **/
    1487             : /**                ERROR ESTIMATE FOR THE ROOTS                    **/
    1488             : /**                                                                **/
    1489             : /********************************************************************/
    1490             : 
    1491             : static GEN
    1492     1828552 : root_error(long n, long k, GEN roots_pol, long err, GEN shatzle)
    1493             : {
    1494     1828552 :   GEN rho, d, eps, epsbis, eps2, aux, rap = NULL;
    1495             :   long i, j;
    1496             : 
    1497     1828552 :   d = cgetg(n+1,t_VEC);
    1498    11732672 :   for (i=1; i<=n; i++)
    1499             :   {
    1500     9904526 :     if (i!=k)
    1501             :     {
    1502     8076149 :       aux = gsub(gel(roots_pol,i), gel(roots_pol,k));
    1503     8075803 :       gel(d,i) = gabs(mygprec(aux,31), DEFAULTPREC);
    1504             :     }
    1505             :   }
    1506     1828146 :   rho = gabs(mygprec(gel(roots_pol,k),31), DEFAULTPREC);
    1507     1828532 :   if (expo(rho) < 0) rho = real_1(DEFAULTPREC);
    1508     1828532 :   eps = mulrr(rho, shatzle);
    1509     1828473 :   aux = shiftr(powru(rho,n), err);
    1510             : 
    1511     5551902 :   for (j=1; j<=2 || (j<=5 && cmprr(rap, dbltor(1.2)) > 0); j++)
    1512             :   {
    1513     3723505 :     GEN prod = NULL; /* 1. */
    1514     3723505 :     long m = n;
    1515     3723505 :     epsbis = mulrr(eps, dbltor(1.25));
    1516    25709610 :     for (i=1; i<=n; i++)
    1517             :     {
    1518    21985852 :       if (i != k && cmprr(gel(d,i),epsbis) > 0)
    1519             :       {
    1520    18223043 :         GEN dif = subrr(gel(d,i),eps);
    1521    18220963 :         prod = prod? mulrr(prod, dif): dif;
    1522    18221971 :         m--;
    1523             :       }
    1524             :     }
    1525     3723758 :     eps2 = prod? divrr(aux, prod): aux;
    1526     3723552 :     if (m > 1) eps2 = sqrtnr(shiftr(eps2, 2*m-2), m);
    1527     3723552 :     rap = divrr(eps,eps2); eps = eps2;
    1528             :   }
    1529     1828348 :   return eps;
    1530             : }
    1531             : 
    1532             : /* round a complex or real number x to an absolute value of 2^(-bit) */
    1533             : static GEN
    1534     4173846 : mygprec_absolute(GEN x, long bit)
    1535             : {
    1536             :   long e;
    1537             :   GEN y;
    1538             : 
    1539     4173846 :   switch(typ(x))
    1540             :   {
    1541     2855060 :     case t_REAL:
    1542     2855060 :       e = expo(x) + bit;
    1543     2855060 :       return (e <= 0 || !signe(x))? real_0_bit(-bit): rtor(x, nbits2prec(e));
    1544     1189201 :     case t_COMPLEX:
    1545     1189201 :       if (gexpo(gel(x,2)) < -bit) return mygprec_absolute(gel(x,1),bit);
    1546     1156166 :       y = cgetg(3,t_COMPLEX);
    1547     1156169 :       gel(y,1) = mygprec_absolute(gel(x,1),bit);
    1548     1156081 :       gel(y,2) = mygprec_absolute(gel(x,2),bit);
    1549     1156170 :       return y;
    1550      129585 :     default: return x;
    1551             :   }
    1552             : }
    1553             : 
    1554             : static long
    1555      510651 : a_posteriori_errors(GEN p, GEN roots_pol, long err)
    1556             : {
    1557      510651 :   long i, n = degpol(p), e_max = -(long)EXPOBITS;
    1558             :   GEN sigma, shatzle;
    1559             : 
    1560      510651 :   err += (long)log2((double)n) + 1;
    1561      510651 :   if (err > -2) return 0;
    1562      510651 :   sigma = real2n(-err, LOWDEFAULTPREC);
    1563             :   /*  2 / ((s - 1)^(1/n) - 1) */
    1564      510650 :   shatzle = divur(2, subrs(sqrtnr(subrs(sigma,1),n), 1));
    1565     2339202 :   for (i=1; i<=n; i++)
    1566             :   {
    1567     1828542 :     pari_sp av = avma;
    1568     1828542 :     GEN x = root_error(n,i,roots_pol,err,shatzle);
    1569     1828345 :     long e = gexpo(x);
    1570     1828402 :     set_avma(av); if (e > e_max) e_max = e;
    1571     1828482 :     gel(roots_pol,i) = mygprec_absolute(gel(roots_pol,i), -e);
    1572             :   }
    1573      510660 :   return e_max;
    1574             : }
    1575             : 
    1576             : /********************************************************************/
    1577             : /**                                                                **/
    1578             : /**                           MAIN                                 **/
    1579             : /**                                                                **/
    1580             : /********************************************************************/
    1581             : static GEN
    1582     1536771 : append_clone(GEN r, GEN a) { a = gclone(a); vectrunc_append(r, a); return a; }
    1583             : 
    1584             : /* put roots in placeholder roots_pol so that |P - L_1...L_n| < 2^(-bit)|P|
    1585             :  * returns prod (x-roots_pol[i]) */
    1586             : static GEN
    1587     1452747 : split_complete(GEN p, long bit, GEN roots_pol)
    1588             : {
    1589     1452747 :   long n = degpol(p);
    1590             :   pari_sp ltop;
    1591             :   GEN p1, F, G, a, b, m1, m2;
    1592             : 
    1593     1452744 :   if (n == 1)
    1594             :   {
    1595      426611 :     a = gneg_i(gdiv(gel(p,2), gel(p,3)));
    1596      426620 :     (void)append_clone(roots_pol,a); return p;
    1597             :   }
    1598     1026133 :   ltop = avma;
    1599     1026133 :   if (n == 2)
    1600             :   {
    1601      555085 :     F = gsub(gsqr(gel(p,3)), gmul2n(gmul(gel(p,2),gel(p,4)), 2));
    1602      555078 :     F = gsqrt(F, nbits2prec(bit));
    1603      555077 :     p1 = ginv(gmul2n(gel(p,4),1));
    1604      555077 :     a = gneg_i(gmul(gadd(F,gel(p,3)), p1));
    1605      555083 :     b =        gmul(gsub(F,gel(p,3)), p1);
    1606      555080 :     a = append_clone(roots_pol,a);
    1607      555084 :     b = append_clone(roots_pol,b); set_avma(ltop);
    1608      555086 :     a = mygprec(a, 3*bit);
    1609      555086 :     b = mygprec(b, 3*bit);
    1610      555088 :     return gmul(gel(p,4), mkpoln(3, gen_1, gneg(gadd(a,b)), gmul(a,b)));
    1611             :   }
    1612      471048 :   split_0(p,bit,&F,&G);
    1613      471049 :   m1 = split_complete(F,bit,roots_pol);
    1614      471049 :   m2 = split_complete(G,bit,roots_pol);
    1615      471043 :   return gerepileupto(ltop, gmul(m1,m2));
    1616             : }
    1617             : 
    1618             : static GEN
    1619     4052217 : quicktofp(GEN x)
    1620             : {
    1621     4052217 :   const long prec = DEFAULTPREC;
    1622     4052217 :   switch(typ(x))
    1623             :   {
    1624     4031772 :     case t_INT: return itor(x, prec);
    1625        8809 :     case t_REAL: return rtor(x, prec);
    1626           0 :     case t_FRAC: return fractor(x, prec);
    1627       11636 :     case t_COMPLEX: {
    1628       11636 :       GEN a = gel(x,1), b = gel(x,2);
    1629             :       /* avoid problem with 0, e.g. x = 0 + I*1e-100. We don't want |x| = 0. */
    1630       11636 :       if (isintzero(a)) return cxcompotor(b, prec);
    1631       11594 :       if (isintzero(b)) return cxcompotor(a, prec);
    1632       11594 :       a = cxcompotor(a, prec);
    1633       11594 :       b = cxcompotor(b, prec); return sqrtr(addrr(sqrr(a), sqrr(b)));
    1634             :     }
    1635           0 :     default: pari_err_TYPE("quicktofp",x);
    1636             :       return NULL;/*LCOV_EXCL_LINE*/
    1637             :   }
    1638             : 
    1639             : }
    1640             : 
    1641             : /* bound log_2 |largest root of p| (Fujiwara's bound) */
    1642             : double
    1643     1097596 : fujiwara_bound(GEN p)
    1644             : {
    1645     1097596 :   pari_sp av = avma;
    1646     1097596 :   long i, n = degpol(p);
    1647             :   GEN cc;
    1648             :   double loglc, Lmax;
    1649             : 
    1650     1097595 :   if (n <= 0) pari_err_CONSTPOL("fujiwara_bound");
    1651     1097595 :   loglc = mydbllog2r( quicktofp(gel(p,n+2)) ); /* log_2 |lc(p)| */
    1652     1097564 :   cc = gel(p, 2);
    1653     1097564 :   if (gequal0(cc))
    1654      144274 :     Lmax = -pariINFINITY-1;
    1655             :   else
    1656      953295 :     Lmax = (mydbllog2r(quicktofp(cc)) - loglc - 1) / n;
    1657     3887606 :   for (i = 1; i < n; i++)
    1658             :   {
    1659     2790030 :     GEN y = gel(p,i+2);
    1660             :     double L;
    1661     2790030 :     if (gequal0(y)) continue;
    1662     2001416 :     L = (mydbllog2r(quicktofp(y)) - loglc) / (n-i);
    1663     2001421 :     if (L > Lmax) Lmax = L;
    1664             :   }
    1665     1097576 :   return gc_double(av, Lmax+1);
    1666             : }
    1667             : 
    1668             : /* Fujiwara's bound, real roots. Based on the following remark: if
    1669             :  *   p = x^n + sum a_i x^i and q = x^n + sum min(a_i,0)x^i
    1670             :  * then for all x >= 0, p(x) >= q(x). Thus any bound for the (positive) roots
    1671             :  * of q is a bound for the positive roots of p. */
    1672             : double
    1673      269190 : fujiwara_bound_real(GEN p, long sign)
    1674             : {
    1675      269190 :   pari_sp av = avma;
    1676             :   GEN x;
    1677      269190 :   long n = degpol(p), i, signodd, signeven;
    1678      269189 :   if (n <= 0) pari_err_CONSTPOL("fujiwara_bound");
    1679      269189 :   x = shallowcopy(p);
    1680      269189 :   if (gsigne(gel(x, n+2)) > 0)
    1681      269168 :   { signeven = 1; signodd = sign; }
    1682             :   else
    1683          21 :   { signeven = -1; signodd = -sign; }
    1684     1238712 :   for (i = 0; i < n; i++)
    1685             :   {
    1686      969521 :     if ((n - i) % 2)
    1687      536118 :     { if (gsigne(gel(x, i+2)) == signodd ) gel(x, i+2) = gen_0; }
    1688             :     else
    1689      433403 :     { if (gsigne(gel(x, i+2)) == signeven) gel(x, i+2) = gen_0; }
    1690             :   }
    1691      269191 :   return gc_double(av, fujiwara_bound(x));
    1692             : }
    1693             : 
    1694             : static GEN
    1695     2072241 : mygprecrc_special(GEN x, long prec, long e)
    1696             : {
    1697             :   GEN y;
    1698     2072241 :   switch(typ(x))
    1699             :   {
    1700       34223 :     case t_REAL:
    1701       34223 :       if (!signe(x)) return real_0_bit(minss(e, expo(x)));
    1702       33243 :       return (prec > realprec(x))? rtor(x, prec): x;
    1703       12399 :     case t_COMPLEX:
    1704       12399 :       y = cgetg(3,t_COMPLEX);
    1705       12399 :       gel(y,1) = mygprecrc_special(gel(x,1),prec,e);
    1706       12399 :       gel(y,2) = mygprecrc_special(gel(x,2),prec,e);
    1707       12399 :       return y;
    1708     2025619 :     default: return x;
    1709             :   }
    1710             : }
    1711             : 
    1712             : /* like mygprec but keep at least the same precision as before */
    1713             : static GEN
    1714      510659 : mygprec_special(GEN x, long bit)
    1715             : {
    1716             :   long lx, i, e, prec;
    1717             :   GEN y;
    1718             : 
    1719      510659 :   if (bit < 0) bit = 0; /* should not happen */
    1720      510659 :   e = gexpo(x) - bit;
    1721      510658 :   prec = nbits2prec(bit);
    1722      510656 :   switch(typ(x))
    1723             :   {
    1724      510656 :     case t_POL:
    1725      510656 :       y = cgetg_copy(x, &lx); y[1] = x[1];
    1726     2558099 :       for (i=2; i<lx; i++) gel(y,i) = mygprecrc_special(gel(x,i),prec,e);
    1727      510656 :       break;
    1728             : 
    1729           0 :     default: y = mygprecrc_special(x,prec,e);
    1730             :   }
    1731      510655 :   return y;
    1732             : }
    1733             : 
    1734             : static GEN
    1735      375334 : fix_roots1(GEN R)
    1736             : {
    1737      375334 :   long i, l = lg(R);
    1738      375334 :   GEN v = cgetg(l, t_VEC);
    1739     1670278 :   for (i=1; i < l; i++) { GEN r = gel(R,i); gel(v,i) = gcopy(r); gunclone(r); }
    1740      375334 :   return v;
    1741             : }
    1742             : static GEN
    1743      510658 : fix_roots(GEN R, long h, long bit)
    1744             : {
    1745             :   long i, j, c, n, prec;
    1746             :   GEN v, Z, gh;
    1747             : 
    1748      510658 :   if (h == 1) return fix_roots1(R);
    1749      135324 :   prec = nbits2prec(bit); Z = grootsof1(h, prec); gh = utoipos(h);
    1750      135324 :   n = lg(R)-1; v = cgetg(h*n + 1, t_VEC);
    1751      377174 :   for (c = i = 1; i <= n; i++)
    1752             :   {
    1753      241848 :     GEN s, r = gel(R,i);
    1754      241848 :     s = (h == 2)? gsqrt(r, prec): gsqrtn(r, gh, NULL, prec);
    1755      775481 :     for (j = 1; j <= h; j++) gel(v, c++) = gmul(s, gel(Z,j));
    1756      241835 :     gunclone(r);
    1757             :   }
    1758      135326 :   return v;
    1759             : }
    1760             : 
    1761             : static GEN
    1762      509845 : all_roots(GEN p, long bit)
    1763             : {
    1764      509845 :   long bit2, i, e, h, n = degpol(p), elc = gexpo(leading_coeff(p));
    1765      509845 :   GEN q, R, m, pd = RgX_deflate_max(p, &h);
    1766      509847 :   double fb = fujiwara_bound(pd);
    1767             :   pari_sp av;
    1768             : 
    1769      509848 :   if (fb < 0) fb = 0;
    1770      509848 :   bit2 = bit + maxss(gexpo(p), 0) + (long)ceil(log2(n / h) + 2 * fb);
    1771      510659 :   for (av = avma, i = 1, e = 0;; i++, set_avma(av))
    1772             :   {
    1773      510659 :     R = vectrunc_init(n+1);
    1774      510659 :     bit2 += e + (n << i);
    1775      510659 :     q = RgX_gtofp_bit(mygprec(pd,bit2), bit2);
    1776      510655 :     q[1] = evalsigne(1)|evalvarn(0);
    1777      510655 :     m = split_complete(q, bit2, R);
    1778      510658 :     R = fix_roots(R, h, bit2);
    1779      510659 :     q = mygprec_special(pd,bit2);
    1780      510655 :     q[1] = evalsigne(1)|evalvarn(0);
    1781      510655 :     e = gexpo(RgX_sub(q, m)) - elc + (long)log2((double)n) + 1;
    1782      510651 :     if (e < 0)
    1783             :     {
    1784      510651 :       if (e < -2*bit2) e = -2*bit2; /* avoid e = -oo */
    1785      510651 :       e = bit + a_posteriori_errors(p, R, e);
    1786      510660 :       if (e < 0) return R;
    1787             :     }
    1788         810 :     if (DEBUGLEVEL)
    1789           0 :       err_printf("all_roots: restarting, i = %ld, e = %ld\n", i,e);
    1790             :   }
    1791             : }
    1792             : 
    1793             : INLINE int
    1794      858319 : isexactscalar(GEN x) { long tx = typ(x); return is_rational_t(tx); }
    1795             : 
    1796             : static int
    1797      220853 : isexactpol(GEN p)
    1798             : {
    1799      220853 :   long i,n = degpol(p);
    1800     1070916 :   for (i=0; i<=n; i++)
    1801      858319 :     if (!isexactscalar(gel(p,i+2))) return 0;
    1802      212597 :   return 1;
    1803             : }
    1804             : 
    1805             : /* p(0) != 0 [for efficiency] */
    1806             : static GEN
    1807      212597 : solve_exact_pol(GEN p, long bit)
    1808             : {
    1809      212597 :   long i, j, k, m, n = degpol(p), iroots = 0;
    1810      212597 :   GEN ex, factors, v = zerovec(n);
    1811             : 
    1812      212597 :   factors = ZX_squff(Q_primpart(p), &ex);
    1813      425194 :   for (i=1; i<lg(factors); i++)
    1814             :   {
    1815      212597 :     GEN roots_fact = all_roots(gel(factors,i), bit);
    1816      212597 :     n = degpol(gel(factors,i));
    1817      212597 :     m = ex[i];
    1818      849419 :     for (j=1; j<=n; j++)
    1819     1273644 :       for (k=1; k<=m; k++) v[++iroots] = roots_fact[j];
    1820             :   }
    1821      212597 :   return v;
    1822             : }
    1823             : 
    1824             : /* return the roots of p with absolute error bit */
    1825             : static GEN
    1826      220853 : roots_com(GEN q, long bit)
    1827             : {
    1828             :   GEN L, p;
    1829      220853 :   long v = RgX_valrem_inexact(q, &p);
    1830      220853 :   int ex = isexactpol(p);
    1831      220853 :   if (!ex) p = RgX_normalize1(p);
    1832      220853 :   if (lg(p) == 3)
    1833           0 :     L = cgetg(1,t_VEC); /* constant polynomial */
    1834             :   else
    1835      220853 :     L = ex? solve_exact_pol(p,bit): all_roots(p,bit);
    1836      220853 :   if (v)
    1837             :   {
    1838        2786 :     GEN M, z, t = gel(q,2);
    1839             :     long i, x, y, l, n;
    1840             : 
    1841        2786 :     if (isrationalzero(t)) x = -bit;
    1842             :     else
    1843             :     {
    1844           7 :       n = gexpo(t);
    1845           7 :       x = n / v; l = degpol(q);
    1846          35 :       for (i = v; i <= l; i++)
    1847             :       {
    1848          28 :         t  = gel(q,i+2);
    1849          28 :         if (isrationalzero(t)) continue;
    1850          28 :         y = (n - gexpo(t)) / i;
    1851          28 :         if (y < x) x = y;
    1852             :       }
    1853             :     }
    1854        2786 :     z = real_0_bit(x); l = v + lg(L);
    1855        2786 :     M = cgetg(l, t_VEC); L -= v;
    1856        5635 :     for (i = 1; i <= v; i++) gel(M,i) = z;
    1857        8379 :     for (     ; i <  l; i++) gel(M,i) = gel(L,i);
    1858        2786 :     L = M;
    1859             :   }
    1860      220853 :   return L;
    1861             : }
    1862             : 
    1863             : static GEN
    1864     1174127 : tocomplex(GEN x, long l, long bit)
    1865             : {
    1866             :   GEN y;
    1867     1174127 :   if (typ(x) == t_COMPLEX)
    1868             :   {
    1869     1154825 :     if (signe(gel(x,1))) return mygprecrc(x, l, -bit);
    1870      136519 :     x = gel(x,2);
    1871      136519 :     y = cgetg(3,t_COMPLEX);
    1872      136522 :     gel(y,1) = real_0_bit(-bit);
    1873      136523 :     gel(y,2) = mygprecrc(x, l, -bit);
    1874             :   }
    1875             :   else
    1876             :   {
    1877       19302 :     y = cgetg(3,t_COMPLEX);
    1878       19302 :     gel(y,1) = mygprecrc(x, l, -bit);
    1879       19302 :     gel(y,2) = real_0_bit(-bit);
    1880             :   }
    1881      155825 :   return y;
    1882             : }
    1883             : 
    1884             : /* x,y are t_COMPLEX of t_REALs or t_REAL, compare wrt |Im x| - |Im y|,
    1885             :  * then Re x - Re y, up to 2^-e absolute error */
    1886             : static int
    1887     2152118 : cmp_complex_appr(void *E, GEN x, GEN y)
    1888             : {
    1889     2152118 :   long e = (long)E;
    1890             :   GEN z, xi, yi, xr, yr;
    1891             :   long sz, sxi, syi;
    1892     2152118 :   if (typ(x) == t_COMPLEX) { xr = gel(x,1); xi = gel(x,2); sxi = signe(xi); }
    1893      791426 :   else { xr = x; xi = NULL; sxi = 0; }
    1894     2152118 :   if (typ(y) == t_COMPLEX) { yr = gel(y,1); yi = gel(y,2); syi = signe(yi); }
    1895      525745 :   else { yr = y; yi = NULL; syi = 0; }
    1896             :   /* Compare absolute values of imaginary parts */
    1897     2152118 :   if (!sxi)
    1898             :   {
    1899      810598 :     if (syi && expo(yi) >= e) return -1;
    1900             :     /* |Im x| ~ |Im y| ~ 0 */
    1901             :   }
    1902     1341520 :   else if (!syi)
    1903             :   {
    1904       48113 :     if (sxi && expo(xi) >= e) return 1;
    1905             :     /* |Im x| ~ |Im y| ~ 0 */
    1906             :   }
    1907             :   else
    1908             :   {
    1909     1293407 :     z = addrr_sign(xi, 1, yi, -1); sz = signe(z);
    1910     1293374 :     if (sz && expo(z) >= e) return (int)sz;
    1911             :   }
    1912             :   /* |Im x| ~ |Im y|, sort according to real parts */
    1913     1282249 :   z = subrr(xr, yr); sz = signe(z);
    1914     1282267 :   if (sz && expo(z) >= e) return (int)sz;
    1915             :   /* Re x ~ Re y. Place negative imaginary part before positive */
    1916      572493 :   return (int) (sxi - syi);
    1917             : }
    1918             : 
    1919             : static GEN
    1920      509892 : clean_roots(GEN L, long l, long bit, long clean)
    1921             : {
    1922      509892 :   long i, n = lg(L), ex = 5 - bit;
    1923      509892 :   GEN res = cgetg(n,t_COL);
    1924     2338014 :   for (i=1; i<n; i++)
    1925             :   {
    1926     1828124 :     GEN c = gel(L,i);
    1927     1828124 :     if (clean && isrealappr(c,ex))
    1928             :     {
    1929      654001 :       if (typ(c) == t_COMPLEX) c = gel(c,1);
    1930      654001 :       c = mygprecrc(c, l, -bit);
    1931             :     }
    1932             :     else
    1933     1174122 :       c = tocomplex(c, l, bit);
    1934     1828122 :     gel(res,i) = c;
    1935             :   }
    1936      509890 :   gen_sort_inplace(res, (void*)ex, &cmp_complex_appr, NULL);
    1937      509888 :   return res;
    1938             : }
    1939             : 
    1940             : /* the vector of roots of p, with absolute error 2^(- prec2nbits(l)) */
    1941             : static GEN
    1942      220881 : roots_aux(GEN p, long l, long clean)
    1943             : {
    1944      220881 :   pari_sp av = avma;
    1945             :   long bit;
    1946             :   GEN L;
    1947             : 
    1948      220881 :   if (typ(p) != t_POL)
    1949             :   {
    1950          21 :     if (gequal0(p)) pari_err_ROOTS0("roots");
    1951          14 :     if (!isvalidcoeff(p)) pari_err_TYPE("roots",p);
    1952           7 :     return cgetg(1,t_COL); /* constant polynomial */
    1953             :   }
    1954      220860 :   if (!signe(p)) pari_err_ROOTS0("roots");
    1955      220860 :   checkvalidpol(p,"roots");
    1956      220853 :   if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
    1957      220853 :   if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
    1958      220853 :   bit = prec2nbits(l);
    1959      220853 :   L = roots_com(p, bit);
    1960      220853 :   return gerepilecopy(av, clean_roots(L, l, bit, clean));
    1961             : }
    1962             : GEN
    1963        7913 : roots(GEN p, long l) { return roots_aux(p,l, 0); }
    1964             : /* clean up roots. If root is real replace it by its real part */
    1965             : GEN
    1966      212968 : cleanroots(GEN p, long l) { return roots_aux(p,l, 1); }
    1967             : 
    1968             : /* private variant of conjvec. Allow non rational coefficients, shallow
    1969             :  * function. */
    1970             : GEN
    1971          84 : polmod_to_embed(GEN x, long prec)
    1972             : {
    1973          84 :   GEN v, T = gel(x,1), A = gel(x,2);
    1974             :   long i, l;
    1975          84 :   if (typ(A) != t_POL || varn(A) != varn(T))
    1976             :   {
    1977           7 :     checkvalidpol(T,"polmod_to_embed");
    1978           7 :     return const_col(degpol(T), A);
    1979             :   }
    1980          77 :   v = cleanroots(T,prec); l = lg(v);
    1981         231 :   for (i=1; i<l; i++) gel(v,i) = poleval(A,gel(v,i));
    1982          77 :   return v;
    1983             : }
    1984             : 
    1985             : GEN
    1986      289036 : QX_complex_roots(GEN p, long l)
    1987             : {
    1988      289036 :   pari_sp av = avma;
    1989             :   long bit, v;
    1990             :   GEN L;
    1991             : 
    1992      289036 :   if (!signe(p)) pari_err_ROOTS0("QX_complex_roots");
    1993      289036 :   if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
    1994      289036 :   if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
    1995      289036 :   bit = prec2nbits(l);
    1996      289035 :   v = RgX_valrem(p, &p);
    1997      289035 :   L = lg(p) > 3? all_roots(Q_primpart(p), bit): cgetg(1,t_COL);
    1998      289039 :   if (v) L = shallowconcat(const_vec(v, real_0_bit(-bit)), L);
    1999      289039 :   return gerepilecopy(av, clean_roots(L, l, bit, 1));
    2000             : }
    2001             : 
    2002             : /********************************************************************/
    2003             : /**                                                                **/
    2004             : /**                REAL ROOTS OF INTEGER POLYNOMIAL                **/
    2005             : /**                                                                **/
    2006             : /********************************************************************/
    2007             : 
    2008             : /* Count sign changes in the coefficients of (x+1)^deg(P)*P(1/(x+1)), P
    2009             :  * has no rational root. The inversion is implicit (we take coefficients
    2010             :  * backwards). */
    2011             : static long
    2012     1276912 : X2XP1(GEN P, GEN *Premapped)
    2013             : {
    2014     1276912 :   const pari_sp av = avma;
    2015     1276912 :   GEN v = shallowcopy(P);
    2016     1276921 :   long i, j, nb, s, dP = degpol(P), vlim = dP+2;
    2017             : 
    2018     8184937 :   for (j = 2; j < vlim; j++) gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
    2019     1276816 :   s = -signe(gel(v, vlim));
    2020     1276816 :   vlim--; nb = 0;
    2021     3805397 :   for (i = 1; i < dP; i++)
    2022             :   {
    2023     3390129 :     long s2 = -signe(gel(v, 2));
    2024     3390129 :     int flag = (s2 == s);
    2025    27709624 :     for (j = 2; j < vlim; j++)
    2026             :     {
    2027    24319471 :       gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
    2028    24319495 :       if (flag) flag = (s2 != signe(gel(v, j+1)));
    2029             :     }
    2030     3390153 :     if (s == signe(gel(v, vlim)))
    2031             :     {
    2032     1075598 :       if (++nb >= 2) return gc_long(av,2);
    2033      736149 :       s = -s;
    2034             :     }
    2035             :     /* if flag is set there will be no further sign changes */
    2036     3050704 :     if (flag && (!Premapped || !nb)) return gc_long(av, nb);
    2037     2528445 :     vlim--;
    2038     2528445 :     if (gc_needed(av, 3))
    2039             :     {
    2040           0 :       if (DEBUGMEM>1) pari_warn(warnmem, "X2XP1, i = %ld/%ld", i, dP-1);
    2041           0 :       if (!Premapped) setlg(v, vlim + 2);
    2042           0 :       v = gerepilecopy(av, v);
    2043             :     }
    2044             :   }
    2045      415268 :   if (vlim >= 2 && s == signe(gel(v, vlim))) nb++;
    2046      415268 :   if (Premapped && nb == 1) *Premapped = v; else set_avma(av);
    2047      415138 :   return nb;
    2048             : }
    2049             : 
    2050             : static long
    2051           0 : _intervalcmp(GEN x, GEN y)
    2052             : {
    2053           0 :   if (typ(x) == t_VEC) x = gel(x, 1);
    2054           0 :   if (typ(y) == t_VEC) y = gel(y, 1);
    2055           0 :   return gcmp(x, y);
    2056             : }
    2057             : 
    2058             : static GEN
    2059    11045079 : _gen_nored(void *E, GEN x) { (void)E; return x; }
    2060             : static GEN
    2061    24370256 : _mp_add(void *E, GEN x, GEN y) { (void)E; return mpadd(x, y); }
    2062             : static GEN
    2063           0 : _mp_sub(void *E, GEN x, GEN y) { (void)E; return mpsub(x, y); }
    2064             : static GEN
    2065     4316917 : _mp_mul(void *E, GEN x, GEN y) { (void)E; return mpmul(x, y); }
    2066             : static GEN
    2067     6228175 : _mp_sqr(void *E, GEN x) { (void)E; return mpsqr(x); }
    2068             : static GEN
    2069    14286973 : _gen_one(void *E) { (void)E; return gen_1; }
    2070             : static GEN
    2071      321027 : _gen_zero(void *E) { (void)E; return gen_0; }
    2072             : 
    2073             : static struct bb_algebra mp_algebra = { _gen_nored, _mp_add, _mp_sub,
    2074             :                          _mp_mul, _mp_sqr, _gen_one, _gen_zero };
    2075             : 
    2076             : static GEN
    2077    34332761 : _mp_cmul(void *E, GEN P, long a, GEN x) {(void)E; return mpmul(gel(P,a+2), x);}
    2078             : 
    2079             : /* Split the polynom P in two parts, whose coeffs have constant sign:
    2080             :  * P(X) = X^D*Pp + Pm. Also compute the two parts of the derivative of P,
    2081             :  * Pprimem = Pm', Pprimep = X*Pp'+ D*Pp => P' = X^(D-1)*Pprimep + Pprimem;
    2082             :  * Pprimep[i] = (i+D) Pp[i]. Return D */
    2083             : static long
    2084      164780 : split_pols(GEN P, GEN *pPp, GEN *pPm, GEN *pPprimep, GEN *pPprimem)
    2085             : {
    2086      164780 :   long i, D, dP = degpol(P), s0 = signe(gel(P,2));
    2087             :   GEN Pp, Pm, Pprimep, Pprimem;
    2088      505686 :   for(i=1; i <= dP; i++)
    2089      505689 :     if (signe(gel(P, i+2)) == -s0) break;
    2090      164780 :   D = i;
    2091      164780 :   Pm = cgetg(D + 2, t_POL);
    2092      164795 :   Pprimem = cgetg(D + 1, t_POL);
    2093      164792 :   Pp = cgetg(dP-D + 3, t_POL);
    2094      164791 :   Pprimep = cgetg(dP-D + 3, t_POL);
    2095      164798 :   Pm[1] = Pp[1] = Pprimem[1] = Pprimep[1] = P[1];
    2096      670486 :   for(i=0; i < D; i++)
    2097             :   {
    2098      505697 :     GEN c = gel(P, i+2);
    2099      505697 :     gel(Pm, i+2) = c;
    2100      505697 :     if (i) gel(Pprimem, i+1) = mului(i, c);
    2101             :   }
    2102      684887 :   for(; i <= dP; i++)
    2103             :   {
    2104      520105 :     GEN c = gel(P, i+2);
    2105      520105 :     gel(Pp, i+2-D) = c;
    2106      520105 :     gel(Pprimep, i+2-D) = mului(i, c);
    2107             :   }
    2108      164782 :   *pPm = normalizepol_lg(Pm, D+2);
    2109      164785 :   *pPprimem = normalizepol_lg(Pprimem, D+1);
    2110      164793 :   *pPp = normalizepol_lg(Pp, dP-D+3);
    2111      164796 :   *pPprimep = normalizepol_lg(Pprimep, dP-D+3);
    2112      164797 :   return dP - degpol(*pPp);
    2113             : }
    2114             : 
    2115             : static GEN
    2116     5166754 : bkeval_single_power(long d, GEN V)
    2117             : {
    2118     5166754 :   long mp = lg(V) - 2;
    2119     5166754 :   if (d > mp) return gmul(gpowgs(gel(V, mp+1), d/mp), gel(V, (d%mp)+1));
    2120     5166754 :   return gel(V, d+1);
    2121             : }
    2122             : 
    2123             : static GEN
    2124     5166397 : splitpoleval(GEN Pp, GEN Pm, GEN pows, long D, long bitprec)
    2125             : {
    2126     5166397 :   GEN vp = gen_bkeval_powers(Pp, degpol(Pp), pows, NULL, &mp_algebra, _mp_cmul);
    2127     5166563 :   GEN vm = gen_bkeval_powers(Pm, degpol(Pm), pows, NULL, &mp_algebra, _mp_cmul);
    2128     5166675 :   GEN xa = bkeval_single_power(D, pows);
    2129             :   GEN r;
    2130     5166740 :   if (!signe(vp)) return vm;
    2131     5166740 :   vp = gmul(vp, xa);
    2132     5165112 :   r = gadd(vp, vm);
    2133     5161280 :   if (gexpo(vp) - (signe(r)? gexpo(r): 0) > prec2nbits(realprec(vp)) - bitprec)
    2134      336897 :     return NULL;
    2135     4825930 :   return r;
    2136             : }
    2137             : 
    2138             : /* optimized Cauchy bound for P = X^D*Pp + Pm, D > deg(Pm) */
    2139             : static GEN
    2140      164797 : splitcauchy(GEN Pp, GEN Pm, long prec)
    2141             : {
    2142      164797 :   GEN S = gel(Pp,2), A = gel(Pm,2);
    2143      164797 :   long i, lPm = lg(Pm), lPp = lg(Pp);
    2144      502693 :   for (i=3; i < lPm; i++) { GEN c = gel(Pm,i); if (abscmpii(A, c) < 0) A = c; }
    2145      520146 :   for (i=3; i < lPp; i++) S = addii(S, gel(Pp, i));
    2146      164796 :   return subsr(1, rdivii(A, S, prec)); /* 1 + |Pm|_oo / |Pp|_1 */
    2147             : }
    2148             : 
    2149             : static GEN
    2150       14442 : ZX_deg1root(GEN P, long prec)
    2151             : {
    2152       14442 :   GEN a = gel(P,3), b = gel(P,2);
    2153       14442 :   if (is_pm1(a))
    2154             :   {
    2155       14442 :     b = itor(b, prec); if (signe(a) > 0) togglesign(b);
    2156       14442 :     return b;
    2157             :   }
    2158           0 :   return rdivii(negi(b), a, prec);
    2159             : }
    2160             : 
    2161             : /* Newton for polynom P, P(0)!=0, with unique sign change => one root in ]0,oo[
    2162             :  * P' has also at most one zero there */
    2163             : static GEN
    2164      164785 : polsolve(GEN P, long bitprec)
    2165             : {
    2166             :   pari_sp av;
    2167             :   GEN Pp, Pm, Pprimep, Pprimem, Pprime, Pprime2, ra, rb, rc, Pc;
    2168      164785 :   long dP = degpol(P), prec = nbits2prec(bitprec);
    2169             :   long expoold, iter, D, rt, s0, bitaddprec, cprec, PREC;
    2170             : 
    2171      164787 :   if (dP == 1) return ZX_deg1root(P, prec);
    2172      164787 :   Pprime = ZX_deriv(P);
    2173      164784 :   Pprime2 = ZX_deriv(Pprime);
    2174      164783 :   bitaddprec = 1 + 2*expu(dP); PREC = prec + nbits2prec(bitaddprec);
    2175      164781 :   D = split_pols(P, &Pp, &Pm, &Pprimep, &Pprimem); /* P = X^D*Pp + Pm */
    2176      164797 :   s0 = signe(gel(P, 2));
    2177      164797 :   rt = maxss(D, brent_kung_optpow(maxss(degpol(Pp), degpol(Pm)), 2, 1));
    2178      164797 :   rb = splitcauchy(Pp, Pm, DEFAULTPREC);
    2179      164789 :   for (cprec = DEFAULTPREC, expoold = LONG_MAX;;)
    2180           0 :   {
    2181      164789 :     GEN pows = gen_powers(rb, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
    2182      164794 :     Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec);
    2183      164788 :     if (!Pc) { cprec += EXTRAPREC64; rb = rtor(rb, cprec); continue; }
    2184      164788 :     if (signe(Pc) != s0) break;
    2185           0 :     shiftr_inplace(rb,1);
    2186             :   }
    2187      164788 :   for (iter = 0, ra = NULL;;)
    2188     1803643 :   {
    2189             :     GEN wdth;
    2190     1968431 :     iter++;
    2191     1968431 :     if (ra)
    2192      895835 :       rc = shiftr(addrr(ra, rb), -1);
    2193             :     else
    2194     1072596 :       rc = shiftr(rb, -1);
    2195             :     for(;;)
    2196           0 :     {
    2197     1968664 :       GEN pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
    2198     1968410 :       Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec+2);
    2199     1968198 :       if (Pc) break;
    2200           0 :       cprec += EXTRAPREC64;
    2201           0 :       rc = rtor(rc, cprec);
    2202             :     }
    2203     1968198 :     if (signe(Pc) == s0)
    2204      588905 :       ra = rc;
    2205             :     else
    2206     1379293 :       rb = rc;
    2207     1968198 :     if (!ra) continue;
    2208     1060364 :     wdth = subrr(rb, ra);
    2209     1060489 :     if (!(iter % 8))
    2210             :     {
    2211      165911 :       GEN m1 = poleval(Pprime, ra), M2;
    2212      165919 :       if (signe(m1) == s0) continue;
    2213      164806 :       M2 = poleval(Pprime2, rb);
    2214      164810 :       if (abscmprr(gmul(M2, wdth), shiftr(m1, 1)) > 0) continue;
    2215      161655 :       break;
    2216             :     }
    2217      894578 :     else if (gexpo(wdth) <= -bitprec)
    2218        3168 :       break;
    2219             :   }
    2220      164823 :   rc = rb; av = avma;
    2221     1353355 :   for(;; rc = gerepileuptoleaf(av, rc))
    2222     1353604 :   {
    2223             :     long exponew;
    2224     1518427 :     GEN Ppc, dist, rcold = rc;
    2225     1518427 :     GEN pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
    2226     1518086 :     Ppc = splitpoleval(Pprimep, Pprimem, pows, D-1, bitaddprec+4);
    2227     1517881 :     if (Ppc) Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec+4);
    2228     1517980 :     if (!Ppc || !Pc)
    2229             :     {
    2230      336905 :       if (cprec >= PREC)
    2231       40913 :         cprec += EXTRAPREC64;
    2232             :       else
    2233      295992 :         cprec = minss(2*cprec, PREC);
    2234      336912 :       rc = rtor(rc, cprec); continue; /* backtrack one step */
    2235             :     }
    2236     1181075 :     dist = typ(Ppc) == t_REAL? divrr(Pc, Ppc): divri(Pc, Ppc);
    2237     1181285 :     rc = subrr(rc, dist);
    2238     1180657 :     if (cmprr(ra, rc) > 0 || cmprr(rb, rc) < 0)
    2239             :     {
    2240           0 :       if (cprec >= PREC) break;
    2241           0 :       cprec = minss(2*cprec, PREC);
    2242           0 :       rc = rtor(rcold, cprec); continue; /* backtrack one step */
    2243             :     }
    2244     1181212 :     if (expoold == LONG_MAX) { expoold = expo(dist); continue; }
    2245      963778 :     exponew = expo(dist);
    2246      963778 :     if (exponew < -bitprec - 1)
    2247             :     {
    2248      230779 :       if (cprec >= PREC) break;
    2249       65987 :       cprec = minss(2*cprec, PREC);
    2250       65990 :       rc = rtor(rc, cprec); continue;
    2251             :     }
    2252      732999 :     if (exponew > expoold - 2)
    2253             :     {
    2254       52664 :       if (cprec >= PREC) break;
    2255       52664 :       expoold = LONG_MAX;
    2256       52664 :       cprec = minss(2*cprec, PREC);
    2257       52665 :       rc = rtor(rc, cprec); continue;
    2258             :     }
    2259      680335 :     expoold = exponew;
    2260             :   }
    2261      164792 :   return rtor(rc, prec);
    2262             : }
    2263             : 
    2264             : /* Return primpart(P(x / 2)) */
    2265             : static GEN
    2266      415677 : ZX_rescale2prim(GEN P)
    2267             : {
    2268      415677 :   long i, l = lg(P), v, n;
    2269             :   GEN Q;
    2270      415677 :   if (l==2) return pol_0(varn(P));
    2271      415677 :   Q = cgetg(l,t_POL); v = vali(gel(P,l-1));
    2272     1379067 :   for (i = l-2, n = 1; v > n && i >= 2; i--, n++)
    2273      963380 :     v = minss(v, vali(gel(P,i)) + n);
    2274      415687 :   gel(Q,l-1) = v? shifti(gel(P,l-1), -v): gel(P,l-1);
    2275     2680919 :   for (i = l-2, n = 1-v; i >= 2; i--, n++)
    2276     2265256 :     gel(Q,i) = shifti(gel(P,i), n);
    2277      415663 :   Q[1] = P[1]; return Q;
    2278             : }
    2279             : 
    2280             : /* assume Q0 has no rational root */
    2281             : static GEN
    2282      229466 : usp(GEN Q0, long flag, long bitprec)
    2283             : {
    2284      229466 :   const pari_sp av = avma;
    2285             :   GEN Qremapped, Q, c, Lc, Lk, sol;
    2286      229466 :   GEN *pQremapped = flag == 1? &Qremapped: NULL;
    2287      229466 :   const long prec = nbits2prec(bitprec), deg = degpol(Q0);
    2288      229464 :   long listsize = 64, nbr = 0, nb_todo, ind, indf, i, k, nb;
    2289             : 
    2290      229464 :   sol = zerocol(deg);
    2291      229468 :   Lc = zerovec(listsize);
    2292      229476 :   Lk = cgetg(listsize+1, t_VECSMALL);
    2293      229475 :   k = Lk[1] = 0;
    2294      229475 :   ind = 1; indf = 2;
    2295      229475 :   Q = Q0;
    2296      229475 :   c = gen_0;
    2297      229475 :   nb_todo = 1;
    2298     1506323 :   while (nb_todo)
    2299             :   {
    2300     1276858 :     GEN nc = gel(Lc, ind);
    2301             :     pari_sp av2;
    2302     1276858 :     if (Lk[ind] == k + 1)
    2303             :     {
    2304      415677 :       Q = Q0 = ZX_rescale2prim(Q0);
    2305      415681 :       c = gen_0;
    2306             :     }
    2307     1276862 :     if (!equalii(nc, c)) Q = ZX_translate(Q, subii(nc, c));
    2308     1276884 :     av2 = avma;
    2309     1276884 :     k = Lk[ind];
    2310     1276884 :     ind++;
    2311     1276884 :     c = nc;
    2312     1276884 :     nb_todo--;
    2313     1276884 :     nb = X2XP1(Q, pQremapped);
    2314             : 
    2315     1276778 :     if (nb == 1)
    2316             :     { /* exactly one root */
    2317      279766 :       GEN s = gen_0;
    2318      279766 :       if (flag == 0)
    2319             :       {
    2320           0 :         s = mkvec2(gmul2n(c,-k), gmul2n(addiu(c,1),-k));
    2321           0 :         s = gerepilecopy(av2, s);
    2322             :       }
    2323      279766 :       else if (flag == 1) /* Caveat: Qremapped is the reciprocal polynomial */
    2324             :       {
    2325      164786 :         s = polsolve(*pQremapped, bitprec+1);
    2326      164798 :         s = addir(c, divrr(s, addsr(1, s)));
    2327      164794 :         shiftr_inplace(s, -k);
    2328      164792 :         if (realprec(s) != prec) s = rtor(s, prec);
    2329      164794 :         s = gerepileupto(av2, s);
    2330             :       }
    2331      114980 :       else set_avma(av2);
    2332      279785 :       gel(sol, ++nbr) = s;
    2333             :     }
    2334      997012 :     else if (nb)
    2335             :     { /* unknown, add two nodes to refine */
    2336      523717 :       if (indf + 2 > listsize)
    2337             :       {
    2338        1055 :         if (ind>1)
    2339             :         {
    2340        3818 :           for (i = ind; i < indf; i++)
    2341             :           {
    2342        2763 :             gel(Lc, i-ind+1) = gel(Lc, i);
    2343        2763 :             Lk[i-ind+1] = Lk[i];
    2344             :           }
    2345        1055 :           indf -= ind-1;
    2346        1055 :           ind = 1;
    2347             :         }
    2348        1055 :         if (indf + 2 > listsize)
    2349             :         {
    2350           0 :           listsize *= 2;
    2351           0 :           Lc = vec_lengthen(Lc, listsize);
    2352           0 :           Lk = vecsmall_lengthen(Lk, listsize);
    2353             :         }
    2354       65812 :         for (i = indf; i <= listsize; i++) gel(Lc, i) = gen_0;
    2355             :       }
    2356      523717 :       gel(Lc, indf) = nc = shifti(c, 1);
    2357      523717 :       gel(Lc, indf + 1) = addiu(nc, 1);
    2358      523716 :       Lk[indf] = Lk[indf + 1] = k + 1;
    2359      523716 :       indf += 2;
    2360      523716 :       nb_todo += 2;
    2361             :     }
    2362     1276796 :     if (gc_needed(av, 2))
    2363             :     {
    2364           0 :       gerepileall(av, 6, &Q0, &Q, &c, &Lc, &Lk, &sol);
    2365           0 :       if (DEBUGMEM > 1) pari_warn(warnmem, "ZX_Uspensky", avma);
    2366             :     }
    2367             :   }
    2368      229465 :   setlg(sol, nbr+1);
    2369      229465 :   return gerepilecopy(av, sol);
    2370             : }
    2371             : 
    2372             : static GEN
    2373          14 : ZX_Uspensky_equal_yes(GEN a, long flag, long bit)
    2374             : {
    2375          14 :   if (flag == 2) return gen_1;
    2376           7 :   if (flag == 1 && typ(a) != t_REAL)
    2377             :   {
    2378           7 :     if (typ(a) == t_INT && !signe(a))
    2379           0 :       a = real_0_bit(bit);
    2380             :     else
    2381           7 :       a = gtofp(a, nbits2prec(bit));
    2382             :   }
    2383           7 :   return mkcol(a);
    2384             : }
    2385             : static GEN
    2386          21 : ZX_Uspensky_no(long flag)
    2387          21 : { return flag <= 1 ? cgetg(1, t_COL) : gen_0; }
    2388             : /* ZX_Uspensky(P, [a,a], flag) */
    2389             : static GEN
    2390          28 : ZX_Uspensky_equal(GEN P, GEN a, long flag, long bit)
    2391             : {
    2392          28 :   if (typ(a) != t_INFINITY && gequal0(poleval(P, a)))
    2393          14 :     return ZX_Uspensky_equal_yes(a, flag, bit);
    2394             :   else
    2395          14 :     return ZX_Uspensky_no(flag);
    2396             : }
    2397             : static int
    2398        3399 : sol_ok(GEN r, GEN a, GEN b) { return gcmp(a, r) <= 0 && gcmp(r, b) <= 0; }
    2399             : 
    2400             : /* P a ZX without real double roots; better if primitive and squarefree but
    2401             :  * caller should ensure that. If flag & 4 assume that P has no rational root
    2402             :  * (modest speedup) */
    2403             : GEN
    2404      182109 : ZX_Uspensky(GEN P, GEN ab, long flag, long bitprec)
    2405             : {
    2406      182109 :   pari_sp av = avma;
    2407             :   GEN a, b, res, sol;
    2408             :   double fb;
    2409             :   long l, nbz, deg;
    2410             : 
    2411      182109 :   if (ab)
    2412             :   {
    2413       80512 :     if (typ(ab) == t_VEC)
    2414             :     {
    2415       53317 :       if (lg(ab) != 3) pari_err_DIM("ZX_Uspensky");
    2416       53317 :       a = gel(ab, 1);
    2417       53317 :       b = gel(ab, 2);
    2418             :     }
    2419             :     else
    2420             :     {
    2421       27195 :       a = ab;
    2422       27195 :       b = mkoo();
    2423             :     }
    2424             :   }
    2425             :   else
    2426             :   {
    2427      101597 :     a = mkmoo();
    2428      101596 :     b = mkoo();
    2429             :   }
    2430      182108 :   if (flag & 4)
    2431             :   {
    2432      127461 :     if (gcmp(a, b) >= 0) { set_avma(av); return ZX_Uspensky_no(flag); }
    2433      127461 :     flag &= ~4;
    2434      127461 :     sol = cgetg(1, t_COL);
    2435             :   }
    2436             :   else
    2437             :   {
    2438       54647 :     switch (gcmp(a, b))
    2439             :     {
    2440           7 :       case 1: set_avma(av); return ZX_Uspensky_no(flag);
    2441          28 :       case 0: return gerepilecopy(av, ZX_Uspensky_equal(P, a, flag, bitprec));
    2442             :     }
    2443       54612 :     sol = nfrootsQ(P);
    2444             :   }
    2445      182073 :   nbz = 0; l = lg(sol);
    2446      182073 :   if (l > 1)
    2447             :   {
    2448             :     long i, j;
    2449        2706 :     P = RgX_div(P, roots_to_pol(sol, varn(P)));
    2450        2706 :     if (!RgV_is_ZV(sol)) P = Q_primpart(P);
    2451        6105 :     for (i = j = 1; i < l; i++)
    2452        3399 :       if (sol_ok(gel(sol,i), a, b)) gel(sol,j++) = gel(sol,i);
    2453        2706 :     setlg(sol, j);
    2454        2706 :     if (flag == 2) { nbz = j-1; sol = utoi(nbz); }
    2455        2538 :     else if (flag == 1) sol = RgC_gtofp(sol, nbits2prec(bitprec));
    2456             :   }
    2457      179367 :   else if (flag == 2) sol = gen_0;
    2458      182073 :   deg = degpol(P);
    2459      182073 :   if (deg == 0) return gerepilecopy(av, sol);
    2460      180112 :   if (typ(a) == t_INFINITY && typ(b) != t_INFINITY && gsigne(b))
    2461             :   {
    2462          28 :     fb = fujiwara_bound_real(P, -1);
    2463          28 :     if (fb <= -pariINFINITY) a = gen_0;
    2464          21 :     else if (fb < 0) a = gen_m1;
    2465          21 :     else a = negi(int2n((long)ceil(fb)));
    2466             :   }
    2467      180112 :   if (typ(b) == t_INFINITY && typ(a) != t_INFINITY && gsigne(a))
    2468             :   {
    2469          21 :     fb = fujiwara_bound_real(P, 1);
    2470          21 :     if (fb <= -pariINFINITY) b = gen_0;
    2471          21 :     else if (fb < 0) b = gen_1;
    2472           7 :     else b = int2n((long)ceil(fb));
    2473             :   }
    2474      180112 :   if (typ(a) != t_INFINITY && typ(b) != t_INFINITY)
    2475             :   {
    2476             :     GEN d, ad, bd, diff;
    2477             :     long i;
    2478             :     /* can occur if one of a,b was initially a t_INFINITY */
    2479       12491 :     if (gequal(a,b)) return gerepilecopy(av, sol);
    2480       12484 :     d = lcmii(Q_denom(a), Q_denom(b));
    2481       12484 :     if (is_pm1(d)) { d = NULL; ad = a; bd = b; }
    2482             :     else
    2483          21 :     { P = ZX_rescale(P, d); ad = gmul(a, d); bd = gmul(b, d); }
    2484       12484 :     diff = subii(bd, ad);
    2485       12484 :     P = ZX_affine(P, diff, ad);
    2486       12484 :     res = usp(P, flag, bitprec);
    2487       12484 :     if (flag <= 1)
    2488             :     {
    2489       34008 :       for (i = 1; i < lg(res); i++)
    2490             :       {
    2491       21818 :         GEN z = gmul(diff, gel(res, i));
    2492       21818 :         if (typ(z) == t_VEC)
    2493             :         {
    2494           0 :           gel(z, 1) = gadd(ad, gel(z, 1));
    2495           0 :           gel(z, 2) = gadd(ad, gel(z, 2));
    2496             :         }
    2497             :         else
    2498       21818 :           z = gadd(ad, z);
    2499       21818 :         if (d) z = gdiv(z, d);
    2500       21818 :         gel(res, i) = z;
    2501             :       }
    2502       12190 :       sol = shallowconcat(sol, res);
    2503             :     }
    2504             :     else
    2505         294 :       nbz += lg(res) - 1;
    2506             :   }
    2507      180105 :   if (typ(b) == t_INFINITY && (fb=fujiwara_bound_real(P, 1)) > -pariINFINITY)
    2508             :   {
    2509      131659 :     long bp = maxss((long)ceil(fb), 0);
    2510      131659 :     res = usp(ZX_unscale2n(P, bp), flag, bitprec);
    2511      131661 :     if (flag <= 1)
    2512       70847 :       sol = shallowconcat(sol, gmul2n(res, bp));
    2513             :     else
    2514       60814 :       nbz += lg(res)-1;
    2515             :   }
    2516      180105 :   if (typ(a) == t_INFINITY && (fb=fujiwara_bound_real(P,-1)) > -pariINFINITY)
    2517             :   {
    2518       85328 :     long i, bm = maxss((long)ceil(fb), 0);
    2519       85328 :     res = usp(ZX_unscale2n(ZX_z_unscale(P, -1), bm), flag, bitprec);
    2520       85330 :     if (flag <= 1)
    2521             :     {
    2522       74345 :       for (i = 1; i < lg(res); i++)
    2523             :       {
    2524       46561 :         GEN z = gneg(gmul2n(gel(res, i), bm));
    2525       46562 :         if (typ(z) == t_VEC) swap(gel(z, 1), gel(z, 2));
    2526       46562 :         gel(res, i) = z;
    2527             :       }
    2528       27784 :       sol = shallowconcat(res, sol);
    2529             :     }
    2530             :     else
    2531       57547 :       nbz += lg(res)-1;
    2532             :   }
    2533      180108 :   if (flag >= 2) return utoi(nbz);
    2534       83122 :   if (flag)
    2535       83122 :     sol = sort(sol);
    2536             :   else
    2537           0 :     sol = gen_sort(sol, (void *)_intervalcmp, cmp_nodata);
    2538       83122 :   return gerepileupto(av, sol);
    2539             : }
    2540             : 
    2541             : /* x a scalar */
    2542             : static GEN
    2543          42 : rootsdeg0(GEN x)
    2544             : {
    2545          42 :   if (!is_real_t(typ(x))) pari_err_TYPE("realroots",x);
    2546          35 :   if (gequal0(x)) pari_err_ROOTS0("realroots");
    2547          14 :   return cgetg(1,t_COL); /* constant polynomial */
    2548             : }
    2549             : static void
    2550      105596 : checkbound(GEN a)
    2551             : {
    2552      105596 :   switch(typ(a))
    2553             :   {
    2554      105596 :     case t_INT: case t_FRAC: case t_INFINITY: break;
    2555           0 :     default: pari_err_TYPE("polrealroots", a);
    2556             :   }
    2557      105596 : }
    2558             : static GEN
    2559       53610 : check_ab(GEN ab)
    2560             : {
    2561             :   GEN a, b;
    2562       53610 :   if (!ab) return NULL;
    2563       52798 :   if (typ(ab) != t_VEC || lg(ab) != 3) pari_err_TYPE("polrootsreal",ab);
    2564       52798 :   a = gel(ab,1); checkbound(a);
    2565       52798 :   b = gel(ab,2); checkbound(b);
    2566       52798 :   if (typ(a) == t_INFINITY && inf_get_sign(a) < 0 &&
    2567        1078 :       typ(b) == t_INFINITY && inf_get_sign(b) > 0) ab = NULL;
    2568       52798 :   return ab;
    2569             : }
    2570             : /* e^(1/h) assuming the h-th root is real, beware that sqrtnr assumes e >= 0 */
    2571             : static GEN
    2572       22214 : _sqrtnr(GEN e, long h)
    2573             : {
    2574             :   long s;
    2575             :   GEN r;
    2576       22214 :   if (h == 2) return sqrtr(e);
    2577          14 :   s = signe(e); setsigne(e, 1); /* e < 0 is possible, implies h is odd */
    2578          14 :   r = sqrtnr(e, h); if (s < 0) setsigne(r, -1);
    2579          14 :   return r;
    2580             : }
    2581             : GEN
    2582       50536 : realroots(GEN P, GEN ab, long prec)
    2583             : {
    2584       50536 :   pari_sp av = avma;
    2585       50536 :   GEN sol = NULL, fa, ex;
    2586             :   long i, j, v, l;
    2587             : 
    2588       50536 :   ab = check_ab(ab);
    2589       50536 :   if (typ(P) != t_POL) return rootsdeg0(P);
    2590       50515 :   if (!RgX_is_ZX(P)) P = RgX_rescale_to_int(P);
    2591       50515 :   switch(degpol(P))
    2592             :   {
    2593          14 :     case -1: return rootsdeg0(gen_0);
    2594           7 :     case 0: return rootsdeg0(gel(P,2));
    2595             :   }
    2596       50494 :   v = ZX_valrem(Q_primpart(P), &P);
    2597       50494 :   fa = ZX_squff(P, &ex); l = lg(fa); sol = cgetg(l + 1, t_VEC);
    2598      102529 :   for (i = 1; i < l; i++)
    2599             :   {
    2600       52035 :     GEN Pi = gel(fa, i), soli, soli2;
    2601             :     long n, h;
    2602       52035 :     if (ab) h = 1; else Pi = ZX_deflate_max(Pi, &h);
    2603       52035 :     soli = ZX_Uspensky(Pi, odd(h)? ab: gen_0, 1, prec2nbits(prec));
    2604       52035 :     n = lg(soli); soli2 = odd(h)? NULL: cgetg(n, t_COL);
    2605      118884 :     for (j = 1; j < n; j++)
    2606             :     {
    2607       66849 :       GEN r = gel(soli, j); /* != 0 */
    2608       66849 :       if (typ(r) != t_REAL) gel(soli, j) = r = gtofp(r, prec);
    2609       66849 :       if (h > 1)
    2610             :       {
    2611          77 :         gel(soli, j) = r = _sqrtnr(r, h);
    2612          77 :         if (soli2) gel(soli2, j) = negr(r);
    2613             :       }
    2614             :     }
    2615       52035 :     if (soli2) soli = shallowconcat(soli, soli2);
    2616       52035 :     if (ex[i] > 1) soli = shallowconcat1( const_vec(ex[i], soli) );
    2617       52035 :     gel(sol, i) = soli;
    2618             :   }
    2619       50494 :   if (v && (!ab || (gsigne(gel(ab,1)) <= 0 && gsigne(gel(ab,2)) >= 0)))
    2620          84 :     gel(sol, i++) = const_col(v, real_0(prec));
    2621       50494 :   setlg(sol, i); if (i == 1) { set_avma(av); return cgetg(1,t_COL); }
    2622       50480 :   return gerepileupto(av, sort(shallowconcat1(sol)));
    2623             : }
    2624             : GEN
    2625       47414 : ZX_realroots_irred(GEN P, long prec)
    2626             : {
    2627       47414 :   long dP = degpol(P), j, n, h;
    2628             :   GEN sol, sol2;
    2629             :   pari_sp av;
    2630       47413 :   if (dP == 1) retmkvec(ZX_deg1root(P, prec));
    2631       44322 :   av = avma; P = ZX_deflate_max(P, &h);
    2632       44322 :   if (h == dP)
    2633             :   {
    2634       11351 :     GEN r = _sqrtnr(ZX_deg1root(P, prec), h);
    2635       11352 :     return gerepilecopy(av, odd(h)? mkvec(r): mkvec2(negr(r), r));
    2636             :   }
    2637       32971 :   sol = ZX_Uspensky(P, odd(h)? NULL: gen_0, 1 | 4, prec2nbits(prec));
    2638       32971 :   n = lg(sol); sol2 = odd(h)? NULL: cgetg(n, t_COL);
    2639      131645 :   for (j = 1; j < n; j++)
    2640             :   {
    2641       98673 :     GEN r = gel(sol, j);
    2642       98673 :     if (typ(r) != t_REAL) gel(sol, j) = r = gtofp(r, prec);
    2643       98673 :     if (h > 1)
    2644             :     {
    2645       10786 :       gel(sol, j) = r = _sqrtnr(r, h);
    2646       10787 :       if (sol2) gel(sol2, j) = negr(r);
    2647             :     }
    2648             :   }
    2649       32972 :   if (sol2) sol = shallowconcat(sol, sol2);
    2650       32972 :   return gerepileupto(av, sort(sol));
    2651             : }
    2652             : 
    2653             : static long
    2654      116089 : ZX_sturm_i(GEN P, long flag)
    2655             : {
    2656             :   pari_sp av;
    2657      116089 :   long h, r, dP = degpol(P);
    2658      116089 :   if (dP == 1) return 1;
    2659      112957 :   av = avma; P = ZX_deflate_max(P, &h);
    2660      112957 :   if (h == dP)
    2661             :   { /* now deg P = 1 */
    2662       17521 :     if (odd(h))
    2663         651 :       r = 1;
    2664             :     else
    2665       16870 :       r = (signe(gel(P,2)) != signe(gel(P,3)))? 2: 0;
    2666       17521 :     return gc_long(av, r);
    2667             :   }
    2668       95436 :   if (odd(h))
    2669       73225 :     r = itou(ZX_Uspensky(P, NULL, flag, 0));
    2670             :   else
    2671       22211 :     r = 2*itou(ZX_Uspensky(P, gen_0, flag, 0));
    2672       95437 :   return gc_long(av,r);
    2673             : }
    2674             : /* P nonconstant, squarefree ZX */
    2675             : long
    2676        3074 : ZX_sturmpart(GEN P, GEN ab)
    2677             : {
    2678        3074 :   pari_sp av = avma;
    2679        3074 :   if (!check_ab(ab)) return ZX_sturm(P);
    2680        1667 :   return gc_long(av, itou(ZX_Uspensky(P, ab, 2, 0)));
    2681             : }
    2682             : /* P nonconstant, squarefree ZX */
    2683             : long
    2684        1414 : ZX_sturm(GEN P) { return ZX_sturm_i(P, 2); }
    2685             : /* P irreducible ZX */
    2686             : long
    2687      114675 : ZX_sturm_irred(GEN P) { return ZX_sturm_i(P, 2 + 4); }

Generated by: LCOV version 1.14