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.12.1 lcov report (development 24988-2584e74448) Lines: 1466 1531 95.8 %
Date: 2020-01-26 05:57:03 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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /*******************************************************************/
      15             : /*                                                                 */
      16             : /*                ROOTS OF COMPLEX POLYNOMIALS                     */
      17             : /*  (original code contributed by Xavier Gourdon, INRIA RR 1852)   */
      18             : /*                                                                 */
      19             : /*******************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : static const double pariINFINITY = 1./0.;
      24             : 
      25             : static long
      26      117620 : isvalidcoeff(GEN x)
      27             : {
      28      117620 :   switch (typ(x))
      29             :   {
      30      102666 :     case t_INT: case t_REAL: case t_FRAC: return 1;
      31       14940 :     case t_COMPLEX: return isvalidcoeff(gel(x,1)) && isvalidcoeff(gel(x,2));
      32             :   }
      33          14 :   return 0;
      34             : }
      35             : 
      36             : static void
      37       16002 : checkvalidpol(GEN p, const char *f)
      38             : {
      39       16002 :   long i,n = lg(p);
      40      103721 :   for (i=2; i<n; i++)
      41       87726 :     if (!isvalidcoeff(gel(p,i))) pari_err_TYPE(f, gel(p,i));
      42       15995 : }
      43             : 
      44             : /********************************************************************/
      45             : /**                                                                **/
      46             : /**                   FAST ARITHMETIC over Z[i]                    **/
      47             : /**                                                                **/
      48             : /********************************************************************/
      49             : 
      50             : static GEN
      51     3418983 : ZX_to_ZiX(GEN Pr, GEN Pi)
      52             : {
      53     3418983 :   long i, lr = lg(Pr), li = lg(Pi), l = maxss(lr, li), m = minss(lr, li);
      54     3418983 :   GEN P = cgetg(l, t_POL);
      55    13633476 :   for(i = 2; i < m; i++)
      56    30590533 :     gel(P,i) = signe(gel(Pi,i)) ? mkcomplex(gel(Pr,i), gel(Pi,i))
      57    20376040 :                                 : gel(Pr,i);
      58     4643309 :   for(     ; i < lr; i++)
      59     1224326 :     gel(P,i) = gel(Pr, i);
      60     3427103 :   for(     ; i < li; i++)
      61        8120 :     gel(P,i) = mkcomplex(gen_0, gel(Pi, i));
      62     3418983 :   return P;
      63             : }
      64             : 
      65             : 
      66             : static GEN
      67    11578592 : ZiX_sqr(GEN P)
      68             : {
      69    11578592 :   pari_sp av = avma;
      70             :   GEN Pr2, Pi2, Qr, Qi;
      71    11578592 :   GEN Pr = real_i(P), Pi = imag_i(P);
      72    11578592 :   if (signe(Pi)==0) return gerepileupto(av, ZX_sqr(Pr));
      73     3427010 :   if (signe(Pr)==0) return gerepileupto(av, ZX_neg(ZX_sqr(Pi)));
      74     3418983 :   Pr2 = ZX_sqr(Pr); Pi2 = ZX_sqr(Pi);
      75     3418983 :   Qr = ZX_sub(Pr2, Pi2);
      76     3418983 :   if (degpol(Pr)==degpol(Pi))
      77     2227264 :     Qi = ZX_sub(ZX_sqr(ZX_add(Pr, Pi)), ZX_add(Pr2, Pi2));
      78             :   else
      79     1191719 :     Qi = ZX_shifti(ZX_mul(Pr, Pi), 1);
      80     3418983 :   return gerepilecopy(av, ZX_to_ZiX(Qr, Qi));
      81             : }
      82             : 
      83             : static GEN
      84     5789296 : graeffe(GEN p)
      85             : {
      86     5789296 :   pari_sp av = avma;
      87             :   GEN p0, p1, s0, s1;
      88     5789296 :   long n = degpol(p);
      89             : 
      90     5789296 :   if (!n) return RgX_copy(p);
      91     5789296 :   RgX_even_odd(p, &p0, &p1);
      92             :   /* p = p0(x^2) + x p1(x^2) */
      93     5789296 :   s0 = ZiX_sqr(p0);
      94     5789296 :   s1 = ZiX_sqr(p1);
      95     5789296 :   return gerepileupto(av, RgX_sub(s0, RgX_shift_shallow(s1,1)));
      96             : }
      97             : 
      98             : GEN
      99        7840 : ZX_graeffe(GEN p)
     100             : {
     101        7840 :   pari_sp av = avma;
     102             :   GEN p0, p1, s0, s1;
     103        7840 :   long n = degpol(p);
     104             : 
     105        7840 :   if (!n) return ZX_copy(p);
     106        7840 :   RgX_even_odd(p, &p0, &p1);
     107             :   /* p = p0(x^2) + x p1(x^2) */
     108        7840 :   s0 = ZX_sqr(p0);
     109        7840 :   s1 = ZX_sqr(p1);
     110        7840 :   return gerepileupto(av, ZX_sub(s0, RgX_shift_shallow(s1,1)));
     111             : }
     112             : GEN
     113          14 : polgraeffe(GEN p)
     114             : {
     115          14 :   pari_sp av = avma;
     116             :   GEN p0, p1, s0, s1;
     117          14 :   long n = degpol(p);
     118             : 
     119          14 :   if (typ(p) != t_POL) pari_err_TYPE("polgraeffe",p);
     120          14 :   n = degpol(p);
     121          14 :   if (!n) return gcopy(p);
     122          14 :   RgX_even_odd(p, &p0, &p1);
     123             :   /* p = p0(x^2) + x p1(x^2) */
     124          14 :   s0 = RgX_sqr(p0);
     125          14 :   s1 = RgX_sqr(p1);
     126          14 :   return gerepileupto(av, RgX_sub(s0, RgX_shift_shallow(s1,1)));
     127             : }
     128             : 
     129             : /********************************************************************/
     130             : /**                                                                **/
     131             : /**                       MODULUS OF ROOTS                         **/
     132             : /**                                                                **/
     133             : /********************************************************************/
     134             : 
     135             : /* Quick approximation to log2(|x|); first define y s.t. |y-x| < 2^-32 then
     136             :  * return y rounded to 2 ulp. In particular, if result < 2^21, absolute error
     137             :  * is bounded by 2^-31. If result > 2^21, it is correct to 2 ulp */
     138             : static double
     139    31540293 : mydbllog2i(GEN x)
     140             : {
     141             : #ifdef LONG_IS_64BIT
     142    27195548 :   const double W = 1/(4294967296. * 4294967296.); /* 2^-64 */
     143             : #else
     144     4344745 :   const double W = 1/4294967296.; /*2^-32*/
     145             : #endif
     146             :   GEN m;
     147    31540293 :   long lx = lgefint(x);
     148             :   double l;
     149    31540293 :   if (lx == 2) return -pariINFINITY;
     150    31417828 :   m = int_MSW(x);
     151    31417828 :   l = (double)(ulong)*m;
     152    31417828 :   if (lx == 3) return log2(l);
     153    13910610 :   l += ((double)(ulong)*int_precW(m)) * W;
     154             :   /* at least m = min(53,BIL) bits are correct in the mantissa, thus log2
     155             :    * is correct with error < log(1 + 2^-m) ~ 2^-m. Adding the correct
     156             :    * exponent BIL(lx-3) causes 1ulp further round-off error */
     157    13910610 :   return log2(l) + (double)(BITS_IN_LONG*(lx-3));
     158             : }
     159             : 
     160             : /* return log(|x|) or -pariINFINITY */
     161             : static double
     162     1675165 : mydbllogr(GEN x) {
     163     1675165 :   if (!signe(x)) return -pariINFINITY;
     164     1675165 :   return M_LN2*dbllog2r(x);
     165             : }
     166             : 
     167             : /* return log2(|x|) or -pariINFINITY */
     168             : static double
     169    14399787 : mydbllog2r(GEN x) {
     170    14399787 :   if (!signe(x)) return -pariINFINITY;
     171    14302529 :   return dbllog2r(x);
     172             : }
     173             : double
     174    50642181 : dbllog2(GEN z)
     175             : {
     176             :   double x, y;
     177    50642181 :   switch(typ(z))
     178             :   {
     179    31537171 :     case t_INT: return mydbllog2i(z);
     180        1561 :     case t_FRAC: return mydbllog2i(gel(z,1))-mydbllog2i(gel(z,2));
     181    13574733 :     case t_REAL: return mydbllog2r(z);
     182             :     default: /*t_COMPLEX*/
     183     5528716 :       x = dbllog2(gel(z,1));
     184     5528716 :       y = dbllog2(gel(z,2));
     185     5528716 :       if (x == -pariINFINITY) return y;
     186     5482172 :       if (y == -pariINFINITY) return x;
     187     5418912 :       if (fabs(x-y) > 10) return maxdd(x,y);
     188     5300882 :       return x + 0.5*log2(1 + exp2(2*(y-x)));
     189             :   }
     190             : }
     191             : static GEN /* beware overflow */
     192      843934 : dblexp(double x) { return fabs(x) < 100.? dbltor(exp(x)): mpexp(dbltor(x)); }
     193             : 
     194             : /* find s such that  A_h <= 2^s <= 2 A_i  for one h and all i < n = deg(p),
     195             :  * with  A_i := (binom(n,i) lc(p) / p_i) ^ 1/(n-i), and  p = sum p_i X^i */
     196             : static long
     197     4537321 : findpower(GEN p)
     198             : {
     199     4537321 :   double x, L, mins = pariINFINITY;
     200     4537321 :   long n = degpol(p),i;
     201             : 
     202     4537321 :   L = dbllog2(gel(p,n+2)); /* log2(lc * binom(n,i)) */
     203    22185528 :   for (i=n-1; i>=0; i--)
     204             :   {
     205    17648207 :     L += log2((double)(i+1) / (double)(n-i));
     206    17648207 :     x = dbllog2(gel(p,i+2));
     207    17648207 :     if (x != -pariINFINITY)
     208             :     {
     209    17542442 :       double s = (L - x) / (double)(n-i);
     210    17542442 :       if (s < mins) mins = s;
     211             :     }
     212             :   }
     213     4537321 :   i = (long)ceil(mins);
     214     4537321 :   if (i - mins > 1 - 1e-12) i--;
     215     4537321 :   return i;
     216             : }
     217             : 
     218             : /* returns the exponent for logmodulus(), from the Newton diagram */
     219             : static long
     220      744782 : newton_polygon(GEN p, long k)
     221             : {
     222      744782 :   pari_sp av = avma;
     223             :   double *logcoef, slope;
     224      744782 :   long n = degpol(p), i, j, h, l, *vertex;
     225             : 
     226      744782 :   logcoef = (double*)stack_malloc_align((n+1)*sizeof(double), sizeof(double));
     227      744782 :   vertex = (long*)new_chunk(n+1);
     228             : 
     229             :   /* vertex[i] = 1 if i a vertex of convex hull, 0 otherwise */
     230      744782 :   for (i=0; i<=n; i++) { logcoef[i] = dbllog2(gel(p,2+i)); vertex[i] = 0; }
     231      744782 :   vertex[0] = 1; /* sentinel */
     232     3063869 :   for (i=0; i < n; i=h)
     233             :   {
     234     2319087 :     slope = logcoef[i+1]-logcoef[i];
     235    10749319 :     for (j = h = i+1; j<=n; j++)
     236             :     {
     237     8430232 :       double pij = (logcoef[j]-logcoef[i])/(double)(j-i);
     238     8430232 :       if (slope < pij) { slope = pij; h = j; }
     239             :     }
     240     2319087 :     vertex[h] = 1;
     241             :   }
     242      744782 :   h = k;   while (!vertex[h]) h++;
     243      744782 :   l = k-1; while (!vertex[l]) l--;
     244      744782 :   set_avma(av);
     245      744782 :   return (long)floor((logcoef[h]-logcoef[l])/(double)(h-l) + 0.5);
     246             : }
     247             : 
     248             : /* change z into z*2^e, where z is real or complex of real */
     249             : static void
     250     5015698 : myshiftrc(GEN z, long e)
     251             : {
     252     5015698 :   if (typ(z)==t_COMPLEX)
     253             :   {
     254     1288721 :     if (signe(gel(z,1))) shiftr_inplace(gel(z,1), e);
     255     1288721 :     if (signe(gel(z,2))) shiftr_inplace(gel(z,2), e);
     256             :   }
     257             :   else
     258     3726977 :     if (signe(z)) shiftr_inplace(z, e);
     259     5015698 : }
     260             : 
     261             : /* return z*2^e, where z is integer or complex of integer (destroy z) */
     262             : static GEN
     263    18975427 : myshiftic(GEN z, long e)
     264             : {
     265    18975427 :   if (typ(z)==t_COMPLEX)
     266             :   {
     267     3620747 :     gel(z,1) = signe(gel(z,1))? mpshift(gel(z,1),e): gen_0;
     268     3620747 :     gel(z,2) = mpshift(gel(z,2),e);
     269     3620747 :     return z;
     270             :   }
     271    15354680 :   return signe(z)? mpshift(z,e): gen_0;
     272             : }
     273             : 
     274             : static GEN
     275      957525 : RgX_gtofp_bit(GEN q, long bit) { return RgX_gtofp(q, nbits2prec(bit)); }
     276             : 
     277             : static GEN
     278    36566784 : mygprecrc(GEN x, long prec, long e)
     279             : {
     280             :   GEN y;
     281    36566784 :   switch(typ(x))
     282             :   {
     283             :     case t_REAL:
     284    27016161 :       if (!signe(x)) return real_0_bit(e);
     285    26340507 :       return realprec(x) == prec? x: rtor(x, prec);
     286             :     case t_COMPLEX:
     287     7457469 :       y = cgetg(3,t_COMPLEX);
     288     7457469 :       gel(y,1) = mygprecrc(gel(x,1),prec,e);
     289     7457469 :       gel(y,2) = mygprecrc(gel(x,2),prec,e);
     290     7457469 :       return y;
     291     2093154 :     default: return x;
     292             :   }
     293             : }
     294             : 
     295             : /* gprec behaves badly with the zero for polynomials.
     296             : The second parameter in mygprec is the precision in base 2 */
     297             : static GEN
     298     8494969 : mygprec(GEN x, long bit)
     299             : {
     300             :   long lx, i, e, prec;
     301             :   GEN y;
     302             : 
     303     8494969 :   if (bit < 0) bit = 0; /* should rarely happen */
     304     8494969 :   e = gexpo(x) - bit;
     305     8494969 :   prec = nbits2prec(bit);
     306     8494969 :   switch(typ(x))
     307             :   {
     308             :     case t_POL:
     309     5631329 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     310     5631329 :       for (i=2; i<lx; i++) gel(y,i) = mygprecrc(gel(x,i),prec,e);
     311     5631329 :       break;
     312             : 
     313     2863640 :     default: y = mygprecrc(x,prec,e);
     314             :   }
     315     8494969 :   return y;
     316             : }
     317             : 
     318             : /* normalize a polynomial p, that is change it with coefficients in Z[i],
     319             : after making product by 2^shift */
     320             : static GEN
     321     2245742 : pol_to_gaussint(GEN p, long shift)
     322             : {
     323     2245742 :   long i, l = lg(p);
     324     2245742 :   GEN q = cgetg(l, t_POL); q[1] = p[1];
     325     2245742 :   for (i=2; i<l; i++) gel(q,i) = gtrunc2n(gel(p,i), shift);
     326     2245742 :   return q;
     327             : }
     328             : 
     329             : /* returns a polynomial q in Z[i][x] keeping bit bits of p */
     330             : static GEN
     331     1758937 : eval_rel_pol(GEN p, long bit)
     332             : {
     333             :   long i;
     334    12108456 :   for (i = 2; i < lg(p); i++)
     335    10349519 :     if (gequal0(gel(p,i))) gel(p,i) = gen_0; /* bad behavior of gexpo */
     336     1758937 :   return pol_to_gaussint(p, bit-gexpo(p)+1);
     337             : }
     338             : 
     339             : /* returns p(R*x)/R^n (in R or R[i]), R = exp(lrho), bit bits of precision */
     340             : static GEN
     341      198933 : homothetie(GEN p, double lrho, long bit)
     342             : {
     343             :   GEN q, r, t, iR;
     344      198933 :   long n = degpol(p), i;
     345             : 
     346      198933 :   iR = mygprec(dblexp(-lrho),bit);
     347      198933 :   q = mygprec(p, bit);
     348      198933 :   r = cgetg(n+3,t_POL); r[1] = p[1];
     349      198933 :   t = iR; r[n+2] = q[n+2];
     350     1183836 :   for (i=n-1; i>0; i--)
     351             :   {
     352      984903 :     gel(r,i+2) = gmul(t, gel(q,i+2));
     353      984903 :     t = mulrr(t, iR);
     354             :   }
     355      198933 :   gel(r,2) = gmul(t, gel(q,2)); return r;
     356             : }
     357             : 
     358             : /* change q in 2^(n*e) p(x*2^(-e)), n=deg(q)  [ ~as above with R = 2^-e ]*/
     359             : static void
     360     1231587 : homothetie2n(GEN p, long e)
     361             : {
     362     1231587 :   if (e)
     363             :   {
     364      934301 :     long i,n = lg(p)-1;
     365      934301 :     for (i=2; i<=n; i++) myshiftrc(gel(p,i), (n-i)*e);
     366             :   }
     367     1231587 : }
     368             : 
     369             : /* return 2^f * 2^(n*e) p(x*2^(-e)), n=deg(q) */
     370             : static void
     371     4050516 : homothetie_gauss(GEN p, long e, long f)
     372             : {
     373     4050516 :   if (e || f)
     374             :   {
     375     3684814 :     long i, n = lg(p)-1;
     376     3684814 :     for (i=2; i<=n; i++) gel(p,i) = myshiftic(gel(p,i), f+(n-i)*e);
     377             :   }
     378     4050516 : }
     379             : 
     380             : /* Lower bound on the modulus of the largest root z_0
     381             :  * k is set to an upper bound for #{z roots, |z-z_0| < eps} */
     382             : static double
     383     4537321 : lower_bound(GEN p, long *k, double eps)
     384             : {
     385     4537321 :   long n = degpol(p), i, j;
     386     4537321 :   pari_sp ltop = avma;
     387             :   GEN a, s, S, ilc;
     388             :   double r, R, rho;
     389             : 
     390     4537321 :   if (n < 4) { *k = n; return 0.; }
     391     1837724 :   S = cgetg(5,t_VEC);
     392     1837724 :   a = cgetg(5,t_VEC); ilc = gdiv(real_1(DEFAULTPREC), gel(p,n+2));
     393     1837724 :   for (i=1; i<=4; i++) gel(a,i) = gmul(ilc,gel(p,n+2-i));
     394             :   /* i = 1 split out from next loop for efficiency and initialization */
     395     1837724 :   s = gel(a,1);
     396     1837724 :   gel(S,1) = gneg(s); /* Newton sum S_i */
     397     1837724 :   rho = r = gtodouble(gabs(s,3));
     398     1837724 :   R = r / n;
     399     7350896 :   for (i=2; i<=4; i++)
     400             :   {
     401     5513172 :     s = gmulsg(i,gel(a,i));
     402     5513172 :     for (j=1; j<i; j++) s = gadd(s, gmul(gel(S,j),gel(a,i-j)));
     403     5513172 :     gel(S,i) = gneg(s); /* Newton sum S_i */
     404     5513172 :     r = gtodouble(gabs(s,3));
     405     5513172 :     if (r > 0.)
     406             :     {
     407     5502318 :       r = exp(log(r/n) / (double)i);
     408     5502318 :       if (r > R) R = r;
     409             :     }
     410             :   }
     411     1837724 :   if (R > 0. && eps < 1.2)
     412     1836309 :     *k = (long)floor((rho/R + n) / (1 + exp(-eps)*cos(eps)));
     413             :   else
     414        1415 :     *k = n;
     415     1837724 :   return gc_double(ltop, R);
     416             : }
     417             : 
     418             : /* return R such that exp(R - tau) <= rho_n(P) <= exp(R + tau)
     419             :  * P(0) != 0 and P non constant */
     420             : static double
     421      486805 : logmax_modulus(GEN p, double tau)
     422             : {
     423             :   GEN r, q, aux, gunr;
     424      486805 :   pari_sp av, ltop = avma;
     425      486805 :   long i,k,n=degpol(p),nn,bit,M,e;
     426      486805 :   double rho,eps, tau2 = (tau > 3.0)? 0.5: tau/6.;
     427             : 
     428      486805 :   r = cgeti(BIGDEFAULTPREC);
     429      486805 :   av = avma;
     430             : 
     431      486805 :   eps = - 1/log(1.5*tau2); /* > 0 */
     432      486805 :   bit = (long) ((double) n*log2(1./tau2)+3*log2((double) n))+1;
     433      486805 :   gunr = real_1_bit(bit+2*n);
     434      486805 :   aux = gdiv(gunr, gel(p,2+n));
     435      486805 :   q = RgX_Rg_mul(p, aux); gel(q,2+n) = gunr;
     436      486805 :   e = findpower(q);
     437      486805 :   homothetie2n(q,e);
     438      486805 :   affsi(e, r);
     439      486805 :   q = pol_to_gaussint(q, bit);
     440      486805 :   M = (long) (log2( log(4.*n) / (2*tau2) )) + 2;
     441      486805 :   nn = n;
     442      486805 :   for (i=0,e=0;;)
     443             :   { /* nn = deg(q) */
     444     8587837 :     rho = lower_bound(q, &k, eps);
     445     4537321 :     if (rho > exp2(-(double)e)) e = (long)-floor(log2(rho));
     446     4537321 :     affii(shifti(addis(r,e), 1), r);
     447     4537321 :     if (++i == M) break;
     448             : 
     449    12151548 :     bit = (long) ((double)k * log2(1./tau2) +
     450     8101032 :                      (double)(nn-k)*log2(1./eps) + 3*log2((double)nn)) + 1;
     451     4050516 :     homothetie_gauss(q, e, bit-(long)floor(dbllog2(gel(q,2+nn))+0.5));
     452     4050516 :     nn -= RgX_valrem(q, &q);
     453     4050516 :     q = gerepileupto(av, graeffe(q));
     454     4050516 :     tau2 *= 1.5; if (tau2 > 0.9) tau2 = 0.5;
     455     4050516 :     eps = -1/log(tau2); /* > 0 */
     456     4050516 :     e = findpower(q);
     457             :   }
     458      486805 :   if (!signe(r)) return gc_double(ltop,0.);
     459      458703 :   r = itor(r, DEFAULTPREC); shiftr_inplace(r, -M);
     460      458703 :   return gc_double(ltop, -rtodbl(r) * M_LN2); /* -log(2) sum e_i 2^-i */
     461             : }
     462             : 
     463             : static GEN
     464       13773 : RgX_normalize1(GEN x)
     465             : {
     466       13773 :   long i, n = lg(x)-1;
     467             :   GEN y;
     468       13787 :   for (i = n; i > 1; i--)
     469       13780 :     if (!gequal0( gel(x,i) )) break;
     470       13773 :   if (i == n) return x;
     471          14 :   pari_warn(warner,"normalizing a polynomial with 0 leading term");
     472          14 :   if (i == 1) pari_err_ROOTS0("roots");
     473          14 :   y = cgetg(i+1, t_POL); y[1] = x[1];
     474          14 :   for (; i > 1; i--) gel(y,i) = gel(x,i);
     475          14 :   return y;
     476             : }
     477             : 
     478             : static GEN
     479        6223 : polrootsbound_i(GEN P, double TAU)
     480             : {
     481        6223 :   pari_sp av = avma;
     482             :   double d;
     483        6223 :   (void)RgX_valrem_inexact(P,&P);
     484        6223 :   P = RgX_normalize1(P);
     485        6223 :   switch(degpol(P))
     486             :   {
     487           7 :     case -1: pari_err_ROOTS0("roots");
     488          56 :     case 0:  set_avma(av); return gen_0;
     489             :   }
     490        6160 :   d = logmax_modulus(P, TAU) + TAU;
     491             :   /* not dblexp: result differs on ARM emulator */
     492        6160 :   return gerepileuptoleaf(av, mpexp(dbltor(d)));
     493             : }
     494             : GEN
     495        6230 : polrootsbound(GEN P, GEN tau)
     496             : {
     497        6230 :   if (typ(P) != t_POL) pari_err_TYPE("polrootsbound",P);
     498        6223 :   checkvalidpol(P, "polrootsbound");
     499        6223 :   return polrootsbound_i(P, tau? gtodouble(tau): 0.01);
     500             : }
     501             : 
     502             : /* log of modulus of the smallest root of p, with relative error tau */
     503             : static double
     504      172962 : logmin_modulus(GEN p, double tau)
     505             : {
     506      172962 :   pari_sp av = avma;
     507      172962 :   if (gequal0(gel(p,2))) return -pariINFINITY;
     508      172962 :   return gc_double(av, - logmax_modulus(RgX_recip_shallow(p),tau));
     509             : }
     510             : 
     511             : /* return the log of the k-th modulus (ascending order) of p, rel. error tau*/
     512             : static double
     513       81454 : logmodulus(GEN p, long k, double tau)
     514             : {
     515             :   GEN q;
     516       81454 :   long i, kk = k, imax, n = degpol(p), nn, bit, e;
     517       81454 :   pari_sp av, ltop=avma;
     518       81454 :   double r, tau2 = tau/6;
     519             : 
     520       81454 :   bit = (long)(n * (2. + log2(3.*n/tau2)));
     521       81454 :   av = avma;
     522       81454 :   q = gprec_w(p, nbits2prec(bit));
     523       81454 :   q = RgX_gtofp_bit(q, bit);
     524       81454 :   e = newton_polygon(q,k);
     525       81454 :   r = (double)e;
     526       81454 :   homothetie2n(q,e);
     527       81454 :   imax = (long)(log2(3./tau) + log2(log(4.*n)))+1;
     528      744782 :   for (i=1; i<imax; i++)
     529             :   {
     530      663328 :     q = eval_rel_pol(q,bit);
     531      663328 :     kk -= RgX_valrem(q, &q);
     532      663328 :     nn = degpol(q);
     533             : 
     534      663328 :     q = gerepileupto(av, graeffe(q));
     535      663328 :     e = newton_polygon(q,kk);
     536      663328 :     r += e / exp2((double)i);
     537      663328 :     q = RgX_gtofp_bit(q, bit);
     538      663328 :     homothetie2n(q,e);
     539             : 
     540      663328 :     tau2 *= 1.5; if (tau2 > 1.) tau2 = 1.;
     541      663328 :     bit = 1 + (long)(nn*(2. + log2(3.*nn/tau2)));
     542             :   }
     543       81454 :   return gc_double(ltop, -r * M_LN2);
     544             : }
     545             : 
     546             : /* return the log of the k-th modulus r_k of p, rel. error tau, knowing that
     547             :  * rmin < r_k < rmax. This information helps because we may reduce precision
     548             :  * quicker */
     549             : static double
     550       81454 : logpre_modulus(GEN p, long k, double tau, double lrmin, double lrmax)
     551             : {
     552             :   GEN q;
     553       81454 :   long n = degpol(p), i, imax, imax2, bit;
     554       81454 :   pari_sp ltop = avma, av;
     555       81454 :   double lrho, aux, tau2 = tau/6.;
     556             : 
     557       81454 :   aux = (lrmax - lrmin) / 2. + 4*tau2;
     558       81454 :   imax = (long) log2(log((double)n)/ aux);
     559       81454 :   if (imax <= 0) return logmodulus(p,k,tau);
     560             : 
     561       79962 :   lrho  = (lrmin + lrmax) / 2;
     562       79962 :   av = avma;
     563       79962 :   bit = (long)(n*(2. + aux / M_LN2 - log2(tau2)));
     564       79962 :   q = homothetie(p, lrho, bit);
     565       79962 :   imax2 = (long)(log2(3./tau * log(4.*n))) + 1;
     566       79962 :   if (imax > imax2) imax = imax2;
     567             : 
     568      243530 :   for (i=0; i<imax; i++)
     569             :   {
     570      163568 :     q = eval_rel_pol(q,bit);
     571      163568 :     q = gerepileupto(av, graeffe(q));
     572      163568 :     aux = 2*aux + 2*tau2;
     573      163568 :     tau2 *= 1.5;
     574      163568 :     bit = (long)(n*(2. + aux / M_LN2 - log2(1-exp(-tau2))));
     575      163568 :     q = RgX_gtofp_bit(q, bit);
     576             :   }
     577       79962 :   aux = exp2((double)imax);
     578       79962 :   return gc_double(ltop, lrho + logmodulus(q,k, aux*tau/3.) / aux);
     579             : }
     580             : 
     581             : static double
     582       98814 : ind_maxlog2(GEN q)
     583             : {
     584       98814 :   long i, k = -1;
     585       98814 :   double L = - pariINFINITY;
     586      267057 :   for (i=0; i<=degpol(q); i++)
     587             :   {
     588      168243 :     double d = dbllog2(gel(q,2+i));
     589      168243 :     if (d > L) { L = d; k = i; }
     590             :   }
     591       98814 :   return k;
     592             : }
     593             : 
     594             : /* Returns k such that r_k e^(-tau) < R < r_{k+1} e^tau.
     595             :  * Assume that l <= k <= n-l */
     596             : static long
     597      118971 : dual_modulus(GEN p, double lrho, double tau, long l)
     598             : {
     599      118971 :   long i, imax, delta_k = 0, n = degpol(p), nn, v2, v, bit, ll = l;
     600      118971 :   double tau2 = tau * 7./8.;
     601      118971 :   pari_sp av = avma;
     602             :   GEN q;
     603             : 
     604      118971 :   bit = 6*n - 5*l + (long)(n*(-log2(tau2) + tau2 * 8./7.));
     605      118971 :   q = homothetie(p, lrho, bit);
     606      118971 :   imax = (long)(log(log(2.*n)/tau2)/log(7./4.)+1);
     607             : 
     608     1030855 :   for (i=0; i<imax; i++)
     609             :   {
     610      932041 :     q = eval_rel_pol(q,bit); v2 = n - degpol(q);
     611      932041 :     v = RgX_valrem(q, &q);
     612      932041 :     ll -= maxss(v, v2); if (ll < 0) ll = 0;
     613             : 
     614      932041 :     nn = degpol(q); delta_k += v;
     615      932041 :     if (!nn) return delta_k;
     616             : 
     617      911884 :     q = gerepileupto(av, graeffe(q));
     618      911884 :     tau2 *= 7./4.;
     619      911884 :     bit = 6*nn - 5*ll + (long)(nn*(-log2(tau2) + tau2 * 8./7.));
     620             :   }
     621       98814 :   return gc_long(av, delta_k + (long)ind_maxlog2(q));
     622             : }
     623             : 
     624             : /********************************************************************/
     625             : /**                                                                **/
     626             : /**              FACTORS THROUGH CIRCLE INTEGRATION                **/
     627             : /**                                                                **/
     628             : /********************************************************************/
     629             : /* l power of 2, W[step*j] = w_j; set f[j] = p(w_j)
     630             :  * if inv, w_j = exp(2IPi*j/l), else exp(-2IPi*j/l) */
     631             : 
     632             : static void
     633        7462 : fft2(GEN W, GEN p, GEN f, long step, long l)
     634             : {
     635             :   pari_sp av;
     636             :   long i, s1, l1, step2;
     637             : 
     638        7462 :   if (l == 2)
     639             :   {
     640        3766 :     gel(f,0) = gadd(gel(p,0), gel(p,step));
     641        3766 :     gel(f,1) = gsub(gel(p,0), gel(p,step)); return;
     642             :   }
     643        3696 :   av = avma;
     644        3696 :   l1 = l>>1; step2 = step<<1;
     645        3696 :   fft2(W,p,          f,   step2,l1);
     646        3696 :   fft2(W,p+step,     f+l1,step2,l1);
     647       32760 :   for (i = s1 = 0; i < l1; i++, s1 += step)
     648             :   {
     649       29064 :     GEN f0 = gel(f,i);
     650       29064 :     GEN f1 = gmul(gel(W,s1), gel(f,i+l1));
     651       29064 :     gel(f,i)    = gadd(f0, f1);
     652       29064 :     gel(f,i+l1) = gsub(f0, f1);
     653             :   }
     654        3696 :   gerepilecoeffs(av, f, l);
     655             : }
     656             : 
     657             : static void
     658     3158389 : fft(GEN W, GEN p, GEN f, long step, long l, long inv)
     659             : {
     660             :   pari_sp av;
     661             :   long i, s1, l1, l2, l3, step4;
     662             :   GEN f1, f2, f3, f02;
     663             : 
     664     3158389 :   if (l == 2)
     665             :   {
     666     1746808 :     gel(f,0) = gadd(gel(p,0), gel(p,step));
     667     1746808 :     gel(f,1) = gsub(gel(p,0), gel(p,step)); return;
     668             :   }
     669     1411581 :   av = avma;
     670     1411581 :   if (l == 4)
     671             :   {
     672             :     pari_sp av2;
     673      787894 :     f1 = gadd(gel(p,0), gel(p,step<<1));
     674      787894 :     f2 = gsub(gel(p,0), gel(p,step<<1));
     675      787894 :     f3 = gadd(gel(p,step), gel(p,3*step));
     676      787894 :     f02 = gsub(gel(p,step), gel(p,3*step));
     677      787894 :     f02 = inv? mulcxI(f02): mulcxmI(f02);
     678      787894 :     av2 = avma;
     679      787894 :     gel(f,0) = gadd(f1, f3);
     680      787894 :     gel(f,1) = gadd(f2, f02);
     681      787894 :     gel(f,2) = gsub(f1, f3);
     682      787894 :     gel(f,3) = gsub(f2, f02);
     683      787894 :     gerepileallsp(av,av2,4,&gel(f,0),&gel(f,1),&gel(f,2),&gel(f,3));
     684      787894 :     return;
     685             :   }
     686      623687 :   l1 = l>>2; l2 = 2*l1; l3 = l1+l2; step4 = step<<2;
     687      623687 :   fft(W,p,          f,   step4,l1,inv);
     688      623687 :   fft(W,p+step,     f+l1,step4,l1,inv);
     689      623687 :   fft(W,p+(step<<1),f+l2,step4,l1,inv);
     690      623687 :   fft(W,p+3*step,   f+l3,step4,l1,inv);
     691     2506451 :   for (i = s1 = 0; i < l1; i++, s1 += step)
     692             :   {
     693     1882764 :     long s2 = s1 << 1, s3 = s1 + s2;
     694             :     GEN g02, g13, f13;
     695     1882764 :     f1 = gmul(gel(W,s1), gel(f,i+l1));
     696     1882764 :     f2 = gmul(gel(W,s2), gel(f,i+l2));
     697     1882764 :     f3 = gmul(gel(W,s3), gel(f,i+l3));
     698             : 
     699     1882764 :     f02 = gadd(gel(f,i),f2);
     700     1882764 :     g02 = gsub(gel(f,i),f2);
     701     1882764 :     f13 = gadd(f1,f3);
     702     1882764 :     g13 = gsub(f1,f3); g13 = inv? mulcxI(g13): mulcxmI(g13);
     703             : 
     704     1882764 :     gel(f,i)    = gadd(f02, f13);
     705     1882764 :     gel(f,i+l1) = gadd(g02, g13);
     706     1882764 :     gel(f,i+l2) = gsub(f02, f13);
     707     1882764 :     gel(f,i+l3) = gsub(g02, g13);
     708             :   }
     709      623687 :   gerepilecoeffs(av, f, l);
     710             : }
     711             : 
     712             : #define code(t1,t2) ((t1 << 6) | t2)
     713             : 
     714             : static GEN
     715          84 : FFT_i(GEN W, GEN x)
     716             : {
     717          84 :   long i, l = lg(W), n = lg(x), tx = typ(x), tw, pa;
     718             :   GEN y, z, p, pol;
     719          84 :   if ((l-1) & (l-2)) pari_err_DIM("fft");
     720          70 :   tw = RgV_type(W, &p, &pol, &pa);
     721          70 :   if (tx == t_POL) { x++; n--; }
     722          35 :   else if (!is_vec_t(tx)) pari_err_TYPE("fft",x);
     723          70 :   if (n > l) pari_err_DIM("fft");
     724             : 
     725          70 :   if (n < l) {
     726           0 :     z = cgetg(l, t_VECSMALL); /* cf stackdummy */
     727           0 :     for (i = 1; i < n; i++) gel(z,i) = gel(x,i);
     728           0 :     for (     ; i < l; i++) gel(z,i) = gen_0;
     729             :   }
     730          70 :   else z = x;
     731          70 :   y = cgetg(l, t_VEC);
     732          70 :   if (tw==code(t_COMPLEX,t_INT) || tw==code(t_COMPLEX,t_REAL))
     733           0 :   {
     734           0 :     long inv = (l >= 5 && signe(imag_i(gel(W,1+(l>>2))))==1) ? 1 : 0;
     735           0 :     fft(W+1, z+1, y+1, 1, l-1, inv);
     736             :   } else
     737          70 :     fft2(W+1, z+1, y+1, 1, l-1);
     738          70 :   return y;
     739             : }
     740             : 
     741             : #undef code
     742             : 
     743             : GEN
     744          42 : FFT(GEN W, GEN x)
     745             : {
     746          42 :   if (!is_vec_t(typ(W))) pari_err_TYPE("fft",W);
     747          42 :   return FFT_i(W, x);
     748             : }
     749             : 
     750             : GEN
     751          42 : FFTinv(GEN W, GEN x)
     752             : {
     753          42 :   long l = lg(W), i;
     754             :   GEN w;
     755          42 :   if (!is_vec_t(typ(W))) pari_err_TYPE("fft",W);
     756          42 :   w = cgetg(l, t_VECSMALL); /* cf stackdummy */
     757          42 :   gel(w,1) = gel(W,1); /* w = gconj(W), faster */
     758          42 :   for (i = 2; i < l; i++) gel(w, i) = gel(W, l-i+1);
     759          42 :   return FFT_i(w, x);
     760             : }
     761             : 
     762             : /* returns 1 if p has only real coefficients, 0 else */
     763             : static int
     764      105829 : isreal(GEN p)
     765             : {
     766             :   long i;
     767      600297 :   for (i = lg(p)-1; i > 1; i--)
     768      529663 :     if (typ(gel(p,i)) == t_COMPLEX) return 0;
     769       70634 :   return 1;
     770             : }
     771             : 
     772             : /* x non complex */
     773             : static GEN
     774       70332 : abs_update_r(GEN x, double *mu) {
     775       70332 :   GEN y = gtofp(x, DEFAULTPREC);
     776       70332 :   double ly = mydbllogr(y); if (ly < *mu) *mu = ly;
     777       70332 :   setabssign(y); return y;
     778             : }
     779             : /* return |x|, low accuracy. Set *mu = min(log(y), *mu) */
     780             : static GEN
     781     1489704 : abs_update(GEN x, double *mu) {
     782             :   GEN y, xr, yr;
     783             :   double ly;
     784     1489704 :   if (typ(x) != t_COMPLEX) return abs_update_r(x, mu);
     785     1423310 :   xr = gel(x,1);
     786     1423310 :   yr = gel(x,2);
     787     1423310 :   if (gequal0(xr)) return abs_update_r(yr,mu);
     788     1423128 :   if (gequal0(yr)) return abs_update_r(xr,mu);
     789             :   /* have to treat 0 specially: 0E-10 + 1e-20 = 0E-10 */
     790     1419372 :   xr = gtofp(xr, DEFAULTPREC);
     791     1419372 :   yr = gtofp(yr, DEFAULTPREC);
     792     1419372 :   y = sqrtr(addrr(sqrr(xr), sqrr(yr)));
     793     1419372 :   ly = mydbllogr(y); if (ly < *mu) *mu = ly;
     794     1419372 :   return y;
     795             : }
     796             : 
     797             : static void
     798      113340 : initdft(GEN *Omega, GEN *prim, long N, long Lmax, long bit)
     799             : {
     800      113340 :   long prec = nbits2prec(bit);
     801      113340 :   *Omega = grootsof1(Lmax, prec) + 1;
     802      113340 :   *prim = rootsof1u_cx(N, prec);
     803      113340 : }
     804             : 
     805             : static void
     806       55076 : parameters(GEN p, long *LMAX, double *mu, double *gamma,
     807             :            int polreal, double param, double param2)
     808             : {
     809             :   GEN q, pc, Omega, A, RU, prim, g, TWO;
     810       55076 :   long n = degpol(p), bit, NN, K, i, j, Lmax;
     811       55076 :   pari_sp av2, av = avma;
     812             : 
     813       55076 :   bit = gexpo(p) + (long)param2+8;
     814       55076 :   Lmax = 4; while (Lmax <= n) Lmax <<= 1;
     815       55076 :   NN = (long)(param*3.14)+1; if (NN < Lmax) NN = Lmax;
     816       55076 :   K = NN/Lmax; if (K & 1) K++;
     817       55076 :   NN = Lmax*K;
     818       55076 :   if (polreal) K = K/2+1;
     819             : 
     820       55076 :   initdft(&Omega, &prim, NN, Lmax, bit);
     821       55076 :   q = mygprec(p,bit) + 2;
     822       55076 :   A = cgetg(Lmax+1,t_VEC); A++;
     823       55076 :   pc= cgetg(Lmax+1,t_VEC); pc++;
     824       55076 :   for (i=0; i <= n; i++) gel(pc,i)= gel(q,i);
     825       55076 :   for (   ; i<Lmax; i++) gel(pc,i) = gen_0;
     826             : 
     827       55076 :   *mu = pariINFINITY;
     828       55076 :   g = real_0_bit(-bit);
     829       55076 :   TWO = real2n(1, DEFAULTPREC);
     830       55076 :   av2 = avma;
     831       55076 :   RU = gen_1;
     832      213557 :   for (i=0; i<K; i++)
     833             :   {
     834      158481 :     if (i) {
     835      103405 :       GEN z = RU;
     836      676321 :       for (j=1; j<n; j++)
     837             :       {
     838      572916 :         gel(pc,j) = gmul(gel(q,j),z);
     839      572916 :         z = gmul(z,RU); /* RU = prim^i, z=prim^(ij) */
     840             :       }
     841      103405 :       gel(pc,n) = gmul(gel(q,n),z);
     842             :     }
     843             : 
     844      158481 :     fft(Omega,pc,A,1,Lmax,1);
     845      179806 :     if (polreal && i>0 && i<K-1)
     846       21325 :       for (j=0; j<Lmax; j++) g = addrr(g, divrr(TWO, abs_update(gel(A,j),mu)));
     847             :     else
     848      137156 :       for (j=0; j<Lmax; j++) g = addrr(g, invr(abs_update(gel(A,j),mu)));
     849      158481 :     RU = gmul(RU, prim);
     850      158481 :     if (gc_needed(av,1))
     851             :     {
     852           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"parameters");
     853           0 :       gerepileall(av2,2, &g,&RU);
     854             :     }
     855             :   }
     856       55076 :   *gamma = mydbllog2r(divru(g,NN));
     857       55076 :   *LMAX = Lmax; set_avma(av);
     858       55076 : }
     859             : 
     860             : /* NN is a multiple of Lmax */
     861             : static void
     862       58264 : dft(GEN p, long k, long NN, long Lmax, long bit, GEN F, GEN H, long polreal)
     863             : {
     864             :   GEN Omega, q, qd, pc, pd, A, B, C, RU, aux, U, W, prim, prim2;
     865       58264 :   long n = degpol(p), i, j, K;
     866             :   pari_sp ltop;
     867             : 
     868       58264 :   initdft(&Omega, &prim, NN, Lmax, bit);
     869       58264 :   RU = cgetg(n+2,t_VEC) + 1;
     870             : 
     871       58264 :   K = NN/Lmax; if (polreal) K = K/2+1;
     872       58264 :   q = mygprec(p,bit);
     873       58264 :   qd = RgX_deriv(q);
     874             : 
     875       58264 :   A = cgetg(Lmax+1,t_VEC); A++;
     876       58264 :   B = cgetg(Lmax+1,t_VEC); B++;
     877       58264 :   C = cgetg(Lmax+1,t_VEC); C++;
     878       58264 :   pc = cgetg(Lmax+1,t_VEC); pc++;
     879       58264 :   pd = cgetg(Lmax+1,t_VEC); pd++;
     880       58264 :   gel(pc,0) = gel(q,2);  for (i=n+1; i<Lmax; i++) gel(pc,i) = gen_0;
     881       58264 :   gel(pd,0) = gel(qd,2); for (i=n;   i<Lmax; i++) gel(pd,i) = gen_0;
     882             : 
     883       58264 :   ltop = avma;
     884       58264 :   W = cgetg(k+1,t_VEC);
     885       58264 :   U = cgetg(k+1,t_VEC);
     886       58264 :   for (i=1; i<=k; i++) gel(W,i) = gel(U,i) = gen_0;
     887             : 
     888       58264 :   gel(RU,0) = gen_1;
     889       58264 :   prim2 = gen_1;
     890      184554 :   for (i=0; i<K; i++)
     891             :   {
     892      126290 :     gel(RU,1) = prim2;
     893      126290 :     for (j=1; j<n; j++) gel(RU,j+1) = gmul(gel(RU,j),prim2);
     894             :     /* RU[j] = prim^(ij)= prim2^j */
     895             : 
     896      126290 :     for (j=1; j<n; j++) gel(pd,j) = gmul(gel(qd,j+2),gel(RU,j));
     897      126290 :     fft(Omega,pd,A,1,Lmax,1);
     898      126290 :     for (j=1; j<=n; j++) gel(pc,j) = gmul(gel(q,j+2),gel(RU,j));
     899      126290 :     fft(Omega,pc,B,1,Lmax,1);
     900      126290 :     for (j=0; j<Lmax; j++) gel(C,j) = ginv(gel(B,j));
     901      126290 :     for (j=0; j<Lmax; j++) gel(B,j) = gmul(gel(A,j),gel(C,j));
     902      126290 :     fft(Omega,B,A,1,Lmax,1);
     903      126290 :     fft(Omega,C,B,1,Lmax,1);
     904             : 
     905      126290 :     if (polreal) /* p has real coefficients */
     906             :     {
     907       89346 :       if (i>0 && i<K-1)
     908             :       {
     909       50531 :         for (j=1; j<=k; j++)
     910             :         {
     911       42087 :           gel(W,j) = gadd(gel(W,j), gshift(mulreal(gel(A,j+1),gel(RU,j+1)),1));
     912       42087 :           gel(U,j) = gadd(gel(U,j), gshift(mulreal(gel(B,j),gel(RU,j)),1));
     913             :         }
     914             :       }
     915             :       else
     916             :       {
     917      229458 :         for (j=1; j<=k; j++)
     918             :         {
     919      157000 :           gel(W,j) = gadd(gel(W,j), mulreal(gel(A,j+1),gel(RU,j+1)));
     920      157000 :           gel(U,j) = gadd(gel(U,j), mulreal(gel(B,j),gel(RU,j)));
     921             :         }
     922             :       }
     923             :     }
     924             :     else
     925             :     {
     926      123172 :       for (j=1; j<=k; j++)
     927             :       {
     928       77784 :         gel(W,j) = gadd(gel(W,j), gmul(gel(A,j+1),gel(RU,j+1)));
     929       77784 :         gel(U,j) = gadd(gel(U,j), gmul(gel(B,j),gel(RU,j)));
     930             :       }
     931             :     }
     932      126290 :     prim2 = gmul(prim2,prim);
     933      126290 :     gerepileall(ltop,3, &W,&U,&prim2);
     934             :   }
     935             : 
     936      172333 :   for (i=1; i<=k; i++)
     937             :   {
     938      114069 :     aux=gel(W,i);
     939      114069 :     for (j=1; j<i; j++) aux = gadd(aux, gmul(gel(W,i-j),gel(F,k+2-j)));
     940      114069 :     gel(F,k+2-i) = gdivgs(aux,-i*NN);
     941             :   }
     942      172333 :   for (i=0; i<k; i++)
     943             :   {
     944      114069 :     aux=gel(U,k-i);
     945      114069 :     for (j=1+i; j<k; j++) aux = gadd(aux,gmul(gel(F,2+j),gel(U,j-i)));
     946      114069 :     gel(H,i+2) = gdivgs(aux,NN);
     947             :   }
     948       58264 : }
     949             : 
     950             : #define NEWTON_MAX 10
     951             : static GEN
     952      301843 : refine_H(GEN F, GEN G, GEN HH, long bit, long Sbit)
     953             : {
     954      301843 :   GEN H = HH, D, aux;
     955      301843 :   pari_sp ltop = avma;
     956             :   long error, i, bit1, bit2;
     957             : 
     958      301843 :   D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F)); error = gexpo(D);
     959      301843 :   bit2 = bit + Sbit;
     960      565233 :   for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
     961             :   {
     962      263390 :     if (gc_needed(ltop,1))
     963             :     {
     964           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"refine_H");
     965           0 :       gerepileall(ltop,2, &D,&H);
     966             :     }
     967      263390 :     bit1 = -error + Sbit;
     968      263390 :     aux = RgX_mul(mygprec(H,bit1), mygprec(D,bit1));
     969      263390 :     aux = RgX_rem(mygprec(aux,bit1), mygprec(F,bit1));
     970             : 
     971      263390 :     bit1 = -error*2 + Sbit; if (bit1 > bit2) bit1 = bit2;
     972      263390 :     H = RgX_add(mygprec(H,bit1), aux);
     973      263390 :     D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F));
     974      263390 :     error = gexpo(D); if (error < -bit1) error = -bit1;
     975             :   }
     976      301843 :   if (error > -bit/2) return NULL; /* FAIL */
     977      301709 :   return gerepilecopy(ltop,H);
     978             : }
     979             : 
     980             : /* return 0 if fails, 1 else */
     981             : static long
     982       58264 : refine_F(GEN p, GEN *F, GEN *G, GEN H, long bit, double gamma)
     983             : {
     984       58264 :   GEN f0, FF, GG, r, HH = H;
     985       58264 :   long error, i, bit1 = 0, bit2, Sbit, Sbit2,  enh, normF, normG, n = degpol(p);
     986       58264 :   pari_sp av = avma;
     987             : 
     988       58264 :   FF = *F; GG = RgX_divrem(p, FF, &r);
     989       58264 :   error = gexpo(r); if (error <= -bit) error = 1-bit;
     990       58264 :   normF = gexpo(FF);
     991       58264 :   normG = gexpo(GG);
     992       58264 :   enh = gexpo(H); if (enh < 0) enh = 0;
     993       58264 :   Sbit = normF + 2*normG + enh + (long)(4.*log2((double)n)+gamma) + 1;
     994       58264 :   Sbit2 = enh + 2*(normF+normG) + (long)(2.*gamma+5.*log2((double)n)) + 1;
     995       58264 :   bit2 = bit + Sbit;
     996      359973 :   for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
     997             :   {
     998      301843 :     if (bit1 == bit2 && i >= 2) { Sbit += n; Sbit2 += n; bit2 += n; }
     999      301843 :     if (gc_needed(av,1))
    1000             :     {
    1001           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"refine_F");
    1002           0 :       gerepileall(av,4, &FF,&GG,&r,&HH);
    1003             :     }
    1004             : 
    1005      301843 :     bit1 = -error + Sbit2;
    1006      301843 :     HH = refine_H(mygprec(FF,bit1), mygprec(GG,bit1), mygprec(HH,bit1),
    1007             :                   1-error, Sbit2);
    1008      301843 :     if (!HH) return 0; /* FAIL */
    1009             : 
    1010      301709 :     bit1 = -error + Sbit;
    1011      301709 :     r = RgX_mul(mygprec(HH,bit1), mygprec(r,bit1));
    1012      301709 :     f0 = RgX_rem(mygprec(r,bit1), mygprec(FF,bit1));
    1013             : 
    1014      301709 :     bit1 = -2*error + Sbit; if (bit1 > bit2) bit1 = bit2;
    1015      301709 :     FF = gadd(mygprec(FF,bit1),f0);
    1016             : 
    1017      301709 :     bit1 = -3*error + Sbit; if (bit1 > bit2) bit1 = bit2;
    1018      301709 :     GG = RgX_divrem(mygprec(p,bit1), mygprec(FF,bit1), &r);
    1019      301709 :     error = gexpo(r); if (error < -bit1) error = -bit1;
    1020             :   }
    1021       58130 :   if (error>-bit) return 0; /* FAIL */
    1022       55076 :   *F = FF; *G = GG; return 1;
    1023             : }
    1024             : 
    1025             : /* returns F and G from the unit circle U such that |p-FG|<2^(-bit) |cd|,
    1026             : where cd is the leading coefficient of p */
    1027             : static void
    1028       55076 : split_fromU(GEN p, long k, double delta, long bit,
    1029             :             GEN *F, GEN *G, double param, double param2)
    1030             : {
    1031             :   GEN pp, FF, GG, H;
    1032       55076 :   long n = degpol(p), NN, bit2, Lmax;
    1033       55076 :   int polreal = isreal(p);
    1034             :   pari_sp ltop;
    1035             :   double mu, gamma;
    1036             : 
    1037       55076 :   pp = gdiv(p, gel(p,2+n));
    1038       55076 :   parameters(pp, &Lmax,&mu,&gamma, polreal,param,param2);
    1039             : 
    1040       55076 :   H  = cgetg(k+2,t_POL); H[1] = p[1];
    1041       55076 :   FF = cgetg(k+3,t_POL); FF[1]= p[1];
    1042       55076 :   gel(FF,k+2) = gen_1;
    1043             : 
    1044       55076 :   NN = (long)(0.5/delta); NN |= 1; if (NN < 2) NN = 2;
    1045       55076 :   NN *= Lmax; ltop = avma;
    1046             :   for(;;)
    1047             :   {
    1048       61452 :     bit2 = (long)(((double)NN*delta-mu)/M_LN2) + gexpo(pp) + 8;
    1049       58264 :     dft(pp, k, NN, Lmax, bit2, FF, H, polreal);
    1050       58264 :     if (refine_F(pp,&FF,&GG,H,bit,gamma)) break;
    1051        3188 :     NN <<= 1; set_avma(ltop);
    1052             :   }
    1053       55076 :   *G = gmul(GG,gel(p,2+n)); *F = FF;
    1054       55076 : }
    1055             : 
    1056             : static void
    1057       55076 : optimize_split(GEN p, long k, double delta, long bit,
    1058             :             GEN *F, GEN *G, double param, double param2)
    1059             : {
    1060       55076 :   long n = degpol(p);
    1061             :   GEN FF, GG;
    1062             : 
    1063       55076 :   if (k <= n/2)
    1064       41490 :     split_fromU(p,k,delta,bit,F,G,param,param2);
    1065             :   else
    1066             :   {
    1067       13586 :     split_fromU(RgX_recip_shallow(p),n-k,delta,bit,&FF,&GG,param,param2);
    1068       13586 :     *F = RgX_recip_shallow(GG);
    1069       13586 :     *G = RgX_recip_shallow(FF);
    1070             :   }
    1071       55076 : }
    1072             : 
    1073             : /********************************************************************/
    1074             : /**                                                                **/
    1075             : /**               SEARCH FOR SEPARATING CIRCLE                     **/
    1076             : /**                                                                **/
    1077             : /********************************************************************/
    1078             : 
    1079             : /* return p(2^e*x) *2^(-n*e) */
    1080             : static void
    1081           0 : scalepol2n(GEN p, long e)
    1082             : {
    1083           0 :   long i,n=lg(p)-1;
    1084           0 :   for (i=2; i<=n; i++) gel(p,i) = gmul2n(gel(p,i),(i-n)*e);
    1085           0 : }
    1086             : 
    1087             : /* returns p(x/R)*R^n; assume R is at the correct accuracy */
    1088             : static GEN
    1089      469746 : scalepol(GEN p, GEN R, long bit)
    1090      469746 : { return RgX_rescale(mygprec(p, bit), R); }
    1091             : 
    1092             : /* return (conj(a)X-1)^n * p[ (X-a) / (conj(a)X-1) ] */
    1093             : static GEN
    1094      152322 : conformal_basecase(GEN p, GEN a)
    1095             : {
    1096             :   GEN z, r, ma, ca;
    1097      152322 :   long i, n = degpol(p);
    1098             :   pari_sp av;
    1099             : 
    1100      152322 :   if (n <= 0) return p;
    1101      152322 :   ma = gneg(a); ca = conj_i(a);
    1102      152322 :   av = avma;
    1103      152322 :   z = deg1pol_shallow(ca, gen_m1, 0);
    1104      152322 :   r = scalarpol_shallow(gel(p,2+n), 0);
    1105      526135 :   for (i=n-1; ; i--)
    1106             :   {
    1107      899948 :     r = RgX_addmulXn_shallow(r, gmul(ma,r), 1); /* r *= (X - a) */
    1108      526135 :     r = gadd(r, gmul(z, gel(p,2+i)));
    1109      526135 :     if (i == 0) return gerepileupto(av, r);
    1110      373813 :     z = RgX_addmulXn_shallow(gmul(z,ca), gneg(z), 1); /* z *= conj(a)X - 1 */
    1111      373813 :     if (gc_needed(av,2))
    1112             :     {
    1113           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"conformal_pol (%ld/%ld)",n-i, n);
    1114           0 :       gerepileall(av,2, &r,&z);
    1115             :     }
    1116             :   }
    1117             : }
    1118             : static GEN
    1119      152385 : conformal_pol(GEN p, GEN a)
    1120             : {
    1121      152385 :   pari_sp av = avma;
    1122      152385 :   long d, nR, n = degpol(p), v;
    1123             :   GEN Q, R, S, T;
    1124      152385 :   if (n < 35) return conformal_basecase(p, a);
    1125          63 :   d = (n+1) >> 1; v = varn(p);
    1126          63 :   Q = RgX_shift_shallow(p, -d);
    1127          63 :   R = RgXn_red_shallow(p, d);
    1128          63 :   Q = conformal_pol(Q, a);
    1129          63 :   R = conformal_pol(R, a);
    1130          63 :   S = gpowgs(deg1pol_shallow(gen_1, gneg(a), v), d);
    1131          63 :   T = RgX_recip_shallow(S);
    1132          63 :   if (typ(a) == t_COMPLEX) T = gconj(T);
    1133          63 :   if (odd(d)) T = RgX_neg(T);
    1134             :   /* S = (X - a)^d, T = (conj(a) X - 1)^d */
    1135          63 :   nR = n - degpol(R) - d; /* >= 0 */
    1136          63 :   if (nR) T = RgX_mul(T, gpowgs(deg1pol_shallow(gconj(a), gen_m1, v), nR));
    1137          63 :   return gerepileupto(av, RgX_add(RgX_mul(Q, S), RgX_mul(R, T)));
    1138             : }
    1139             : 
    1140             : static const double UNDEF = -100000.;
    1141             : 
    1142             : static double
    1143       55076 : logradius(double *radii, GEN p, long k, double aux, double *delta)
    1144             : {
    1145       55076 :   long i, n = degpol(p);
    1146             :   double lrho, lrmin, lrmax;
    1147       55076 :   if (k > 1)
    1148             :   {
    1149       35823 :     i = k-1; while (i>0 && radii[i] == UNDEF) i--;
    1150       35823 :     lrmin = logpre_modulus(p,k,aux, radii[i], radii[k]);
    1151             :   }
    1152             :   else /* k=1 */
    1153       19253 :     lrmin = logmin_modulus(p,aux);
    1154       55076 :   radii[k] = lrmin;
    1155             : 
    1156       55076 :   if (k+1<n)
    1157             :   {
    1158       45631 :     i = k+2; while (i<=n && radii[i] == UNDEF) i++;
    1159       45631 :     lrmax = logpre_modulus(p,k+1,aux, radii[k+1], radii[i]);
    1160             :   }
    1161             :   else /* k+1=n */
    1162        9445 :     lrmax = logmax_modulus(p,aux);
    1163       55076 :   radii[k+1] = lrmax;
    1164             : 
    1165       55076 :   lrho = radii[k];
    1166      127028 :   for (i=k-1; i>=1; i--)
    1167             :   {
    1168       71952 :     if (radii[i] == UNDEF || radii[i] > lrho)
    1169       50155 :       radii[i] = lrho;
    1170             :     else
    1171       21797 :       lrho = radii[i];
    1172             :   }
    1173       55076 :   lrho = radii[k+1];
    1174      231006 :   for (i=k+1; i<=n; i++)
    1175             :   {
    1176      175930 :     if (radii[i] == UNDEF || radii[i] < lrho)
    1177      102924 :       radii[i] = lrho;
    1178             :     else
    1179       73006 :       lrho = radii[i];
    1180             :   }
    1181       55076 :   *delta = (lrmax - lrmin) / 2;
    1182       55076 :   if (*delta > 1.) *delta = 1.;
    1183       55076 :   return (lrmin + lrmax) / 2;
    1184             : }
    1185             : 
    1186             : static void
    1187       55076 : update_radius(long n, double *radii, double lrho, double *par, double *par2)
    1188             : {
    1189       55076 :   double t, param = 0., param2 = 0.;
    1190             :   long i;
    1191      358034 :   for (i=1; i<=n; i++)
    1192             :   {
    1193      302958 :     radii[i] -= lrho;
    1194      302958 :     t = fabs(rtodbl( invr(subsr(1, dblexp(radii[i]))) ));
    1195      302958 :     param += t; if (t > 1.) param2 += log2(t);
    1196             :   }
    1197       55076 :   *par = param; *par2 = param2;
    1198       55076 : }
    1199             : 
    1200             : /* apply the conformal mapping then split from U */
    1201             : static void
    1202       50753 : conformal_mapping(double *radii, GEN ctr, GEN p, long k, long bit,
    1203             :                   double aux, GEN *F,GEN *G)
    1204             : {
    1205       50753 :   long bit2, n = degpol(p), i;
    1206       50753 :   pari_sp ltop = avma, av;
    1207             :   GEN q, FF, GG, a, R;
    1208             :   double lrho, delta, param, param2;
    1209             :   /* n * (2.*log2(2.732)+log2(1.5)) + 1 */
    1210       50753 :   bit2 = bit + (long)(n*3.4848775) + 1;
    1211       50753 :   a = sqrtr_abs( utor(3, 2*MEDDEFAULTPREC - 2) );
    1212       50753 :   a = divrs(a, -6);
    1213       50753 :   a = gmul(mygprec(a,bit2), mygprec(ctr,bit2)); /* a = -ctr/2sqrt(3) */
    1214             : 
    1215       50753 :   av = avma;
    1216       50753 :   q = conformal_pol(mygprec(p,bit2), a);
    1217      313852 :   for (i=1; i<=n; i++)
    1218      263099 :     if (radii[i] != UNDEF) /* update array radii */
    1219             :     {
    1220      185461 :       pari_sp av2 = avma;
    1221      185461 :       GEN t, r = dblexp(radii[i]), r2 = sqrr(r);
    1222             :       /* 2(r^2 - 1) / (r^2 - 3(r-1)) */
    1223      185461 :       t = divrr(shiftr((subrs(r2,1)),1), subrr(r2, mulur(3,subrs(r,1))));
    1224      185461 :       radii[i] = mydbllogr(addsr(1,t)) / 2;
    1225      185461 :       set_avma(av2);
    1226             :     }
    1227       50753 :   lrho = logradius(radii, q,k,aux/10., &delta);
    1228       50753 :   update_radius(n, radii, lrho, &param, &param2);
    1229             : 
    1230       50753 :   bit2 += (long)(n * fabs(lrho)/M_LN2 + 1.);
    1231       50753 :   R = mygprec(dblexp(-lrho), bit2);
    1232       50753 :   q = scalepol(q,R,bit2);
    1233       50753 :   gerepileall(av,2, &q,&R);
    1234             : 
    1235       50753 :   optimize_split(q,k,delta,bit2,&FF,&GG,param,param2);
    1236       50753 :   bit2 += n; R = invr(R);
    1237       50753 :   FF = scalepol(FF,R,bit2);
    1238       50753 :   GG = scalepol(GG,R,bit2);
    1239             : 
    1240       50753 :   a = mygprec(a,bit2);
    1241       50753 :   FF = conformal_pol(FF,a);
    1242       50753 :   GG = conformal_pol(GG,a);
    1243             : 
    1244       50753 :   a = invr(subsr(1, gnorm(a)));
    1245       50753 :   FF = RgX_Rg_mul(FF, powru(a,k));
    1246       50753 :   GG = RgX_Rg_mul(GG, powru(a,n-k));
    1247             : 
    1248       50753 :   *F = mygprec(FF,bit+n);
    1249       50753 :   *G = mygprec(GG,bit+n); gerepileall(ltop,2, F,G);
    1250       50753 : }
    1251             : 
    1252             : /* split p, this time without scaling. returns in F and G two polynomials
    1253             :  * such that |p-FG|< 2^(-bit)|p| */
    1254             : static void
    1255       55076 : split_2(GEN p, long bit, GEN ctr, double thickness, GEN *F, GEN *G)
    1256             : {
    1257             :   GEN q, FF, GG, R;
    1258             :   double aux, delta, param, param2;
    1259       55076 :   long n = degpol(p), i, j, k, bit2;
    1260             :   double lrmin, lrmax, lrho, *radii;
    1261             : 
    1262       55076 :   radii = (double*) stack_malloc_align((n+1) * sizeof(double), sizeof(double));
    1263             : 
    1264       55076 :   for (i=2; i<n; i++) radii[i] = UNDEF;
    1265       55076 :   aux = thickness/(double)(4 * n);
    1266       55076 :   lrmin = logmin_modulus(p, aux);
    1267       55076 :   lrmax = logmax_modulus(p, aux);
    1268       55076 :   radii[1] = lrmin;
    1269       55076 :   radii[n] = lrmax;
    1270       55076 :   i = 1; j = n;
    1271       55076 :   lrho = (lrmin + lrmax) / 2;
    1272       55076 :   k = dual_modulus(p, lrho, aux, 1);
    1273       55076 :   if (5*k < n || (n < 2*k && 5*k < 4*n))
    1274       11084 :     { lrmax = lrho; j=k+1; radii[j] = lrho; }
    1275             :   else
    1276       43992 :     { lrmin = lrho; i=k;   radii[i] = lrho; }
    1277      174047 :   while (j > i+1)
    1278             :   {
    1279       63895 :     if (i+j == n+1)
    1280       27395 :       lrho = (lrmin + lrmax) / 2;
    1281             :     else
    1282             :     {
    1283       36500 :       double kappa = 2. - log(1. + minss(i,n-j)) / log(1. + minss(j,n-i));
    1284       36500 :       if (i+j < n+1) lrho = lrmax * kappa + lrmin;
    1285       27428 :       else           lrho = lrmin * kappa + lrmax;
    1286       36500 :       lrho /= 1+kappa;
    1287             :     }
    1288       63895 :     aux = (lrmax - lrmin) / (4*(j-i));
    1289       63895 :     k = dual_modulus(p, lrho, aux, minss(i,n+1-j));
    1290       63895 :     if (k-i < j-k-1 || (k-i == j-k-1 && 2*k > n))
    1291       44205 :       { lrmax = lrho; j=k+1; radii[j] = lrho - aux; }
    1292             :     else
    1293       19690 :       { lrmin = lrho; i=k;   radii[i] = lrho + aux; }
    1294             :   }
    1295       55076 :   aux = lrmax - lrmin;
    1296             : 
    1297       55076 :   if (ctr)
    1298             :   {
    1299       50753 :     lrho = (lrmax + lrmin) / 2;
    1300      313852 :     for (i=1; i<=n; i++)
    1301      263099 :       if (radii[i] != UNDEF) radii[i] -= lrho;
    1302             : 
    1303       50753 :     bit2 = bit + (long)(n * fabs(lrho)/M_LN2 + 1.);
    1304       50753 :     R = mygprec(dblexp(-lrho), bit2);
    1305       50753 :     q = scalepol(p,R,bit2);
    1306       50753 :     conformal_mapping(radii, ctr, q, k, bit2, aux, &FF, &GG);
    1307             :   }
    1308             :   else
    1309             :   {
    1310        4323 :     lrho = logradius(radii, p, k, aux/10., &delta);
    1311        4323 :     update_radius(n, radii, lrho, &param, &param2);
    1312             : 
    1313        4323 :     bit2 = bit + (long)(n * fabs(lrho)/M_LN2 + 1.);
    1314        4323 :     R = mygprec(dblexp(-lrho), bit2);
    1315        4323 :     q = scalepol(p,R,bit2);
    1316        4323 :     optimize_split(q, k, delta, bit2, &FF, &GG, param, param2);
    1317             :   }
    1318       55076 :   bit  += n;
    1319       55076 :   bit2 += n; R = invr(mygprec(R,bit2));
    1320       55076 :   *F = mygprec(scalepol(FF,R,bit2), bit);
    1321       55076 :   *G = mygprec(scalepol(GG,R,bit2), bit);
    1322       55076 : }
    1323             : 
    1324             : /* procedure corresponding to steps 5,6,.. page 44 in RR n. 1852 */
    1325             : /* put in F and G two polynomial such that |p-FG|<2^(-bit)|p|
    1326             :  * where the maximum modulus of the roots of p is <=1.
    1327             :  * Assume sum of roots is 0. */
    1328             : static void
    1329       50753 : split_1(GEN p, long bit, GEN *F, GEN *G)
    1330             : {
    1331       50753 :   long i, imax, n = degpol(p), polreal = isreal(p), ep = gexpo(p), bit2 = bit+n;
    1332             :   GEN ctr, q, qq, FF, GG, v, gr, r, newq;
    1333             :   double lrmin, lrmax, lthick;
    1334       50753 :   const double LOG3 = 1.098613;
    1335             : 
    1336       50753 :   lrmax = logmax_modulus(p, 0.01);
    1337       50753 :   gr = mygprec(dblexp(-lrmax), bit2);
    1338       50753 :   q = scalepol(p,gr,bit2);
    1339             : 
    1340       50753 :   bit2 = bit + gexpo(q) - ep + (long)((double)n*2.*log2(3.)+1);
    1341       50753 :   v = cgetg(5,t_VEC);
    1342       50753 :   gel(v,1) = gen_2;
    1343       50753 :   gel(v,2) = gen_m2;
    1344       50753 :   gel(v,3) = mkcomplex(gen_0, gel(v,1));
    1345       50753 :   gel(v,4) = mkcomplex(gen_0, gel(v,2));
    1346       50753 :   q = mygprec(q,bit2); lthick = 0;
    1347       50753 :   newq = ctr = NULL; /* -Wall */
    1348       50753 :   imax = polreal? 3: 4;
    1349      100351 :   for (i=1; i<=imax; i++)
    1350             :   {
    1351       98633 :     qq = RgX_translate(q, gel(v,i));
    1352       98633 :     lrmin = logmin_modulus(qq,0.05);
    1353       98633 :     if (LOG3 > lrmin + lthick)
    1354             :     {
    1355       95784 :       double lquo = logmax_modulus(qq,0.05) - lrmin;
    1356       95784 :       if (lquo > lthick) { lthick = lquo; newq = qq; ctr = gel(v,i); }
    1357             :     }
    1358       98633 :     if (lthick > M_LN2) break;
    1359       56660 :     if (polreal && i==2 && lthick > LOG3 - M_LN2) break;
    1360             :   }
    1361       50753 :   bit2 = bit + gexpo(newq) - ep + (long)(n*LOG3/M_LN2 + 1);
    1362       50753 :   split_2(newq, bit2, ctr, lthick, &FF, &GG);
    1363       50753 :   r = gneg(mygprec(ctr,bit2));
    1364       50753 :   FF = RgX_translate(FF,r);
    1365       50753 :   GG = RgX_translate(GG,r);
    1366             : 
    1367       50753 :   gr = invr(gr); bit2 = bit - ep + gexpo(FF)+gexpo(GG);
    1368       50753 :   *F = scalepol(FF,gr,bit2);
    1369       50753 :   *G = scalepol(GG,gr,bit2);
    1370       50753 : }
    1371             : 
    1372             : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|,
    1373             : where the maximum modulus of the roots of p is < 0.5 */
    1374             : static int
    1375       50888 : split_0_2(GEN p, long bit, GEN *F, GEN *G)
    1376             : {
    1377             :   GEN q, b;
    1378       50888 :   long n = degpol(p), k, bit2, eq;
    1379       50888 :   double aux0 = dbllog2(gel(p,n+2)); /* != -oo */
    1380       50888 :   double aux1 = dbllog2(gel(p,n+1)), aux;
    1381             : 
    1382       50888 :   if (aux1 == -pariINFINITY) /* p1 = 0 */
    1383        1297 :     aux = 0;
    1384             :   else
    1385             :   {
    1386       49591 :     aux = aux1 - aux0; /* log2(p1/p0) */
    1387             :     /* beware double overflow */
    1388       49591 :     if (aux >= 0 && (aux > 1e4 || exp2(aux) > 2.5*n)) return 0;
    1389       49591 :     aux = (aux < -300)? 0.: n*log2(1 + exp2(aux)/(double)n);
    1390             :   }
    1391       50888 :   bit2 = bit+1 + (long)(log2((double)n) + aux);
    1392       50888 :   q = mygprec(p,bit2);
    1393       50888 :   if (aux1 == -pariINFINITY) b = NULL;
    1394             :   else
    1395             :   {
    1396       49591 :     b = gdivgs(gdiv(gel(q,n+1),gel(q,n+2)),-n);
    1397       49591 :     q = RgX_translate(q,b);
    1398             :   }
    1399       50888 :   gel(q,n+1) = gen_0; eq = gexpo(q);
    1400       50888 :   k = 0;
    1401      102141 :   while (k <= n/2 && (- gexpo(gel(q,k+2)) > bit2 + 2*(n-k) + eq
    1402       51131 :                       || gequal0(gel(q,k+2)))) k++;
    1403       50888 :   if (k > 0)
    1404             :   {
    1405         135 :     if (k > n/2) k = n/2;
    1406         135 :     bit2 += k<<1;
    1407         135 :     *F = pol_xn(k, 0);
    1408         135 :     *G = RgX_shift_shallow(q, -k);
    1409             :   }
    1410             :   else
    1411             :   {
    1412       50753 :     split_1(q,bit2,F,G);
    1413       50753 :     bit2 = bit + gexpo(*F) + gexpo(*G) - gexpo(p) + (long)aux+1;
    1414       50753 :     *F = mygprec(*F,bit2);
    1415             :   }
    1416       50888 :   *G = mygprec(*G,bit2);
    1417       50888 :   if (b)
    1418             :   {
    1419       49591 :     GEN mb = mygprec(gneg(b), bit2);
    1420       49591 :     *F = RgX_translate(*F, mb);
    1421       49591 :     *G = RgX_translate(*G, mb);
    1422             :   }
    1423       50888 :   return 1;
    1424             : }
    1425             : 
    1426             : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|.
    1427             :  * Assume max_modulus(p) < 2 */
    1428             : static void
    1429       50888 : split_0_1(GEN p, long bit, GEN *F, GEN *G)
    1430             : {
    1431             :   GEN FF, GG;
    1432             :   long n, bit2, normp;
    1433             : 
    1434       50888 :   if  (split_0_2(p,bit,F,G)) return;
    1435             : 
    1436           0 :   normp = gexpo(p);
    1437           0 :   scalepol2n(p,2); /* p := 4^(-n) p(4*x) */
    1438           0 :   n = degpol(p); bit2 = bit + 2*n + gexpo(p) - normp;
    1439           0 :   split_1(mygprec(p,bit2), bit2,&FF,&GG);
    1440           0 :   scalepol2n(FF,-2);
    1441           0 :   scalepol2n(GG,-2); bit2 = bit + gexpo(FF) + gexpo(GG) - normp;
    1442           0 :   *F = mygprec(FF,bit2);
    1443           0 :   *G = mygprec(GG,bit2);
    1444             : }
    1445             : 
    1446             : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P| */
    1447             : static void
    1448       55211 : split_0(GEN p, long bit, GEN *F, GEN *G)
    1449             : {
    1450       55211 :   const double LOG1_9 = 0.6418539;
    1451       55211 :   long n = degpol(p), k = 0;
    1452             :   GEN q;
    1453             : 
    1454       55211 :   while (gexpo(gel(p,k+2)) < -bit && k <= n/2) k++;
    1455       55211 :   if (k > 0)
    1456             :   {
    1457           0 :     if (k > n/2) k = n/2;
    1458           0 :     *F = pol_xn(k, 0);
    1459           0 :     *G = RgX_shift_shallow(p, -k);
    1460             :   }
    1461             :   else
    1462             :   {
    1463       55211 :     double lr = logmax_modulus(p, 0.05);
    1464       55211 :     if (lr < LOG1_9) split_0_1(p, bit, F, G);
    1465             :     else
    1466             :     {
    1467       41414 :       q = RgX_recip_shallow(p);
    1468       41414 :       lr = logmax_modulus(q,0.05);
    1469       41414 :       if (lr < LOG1_9)
    1470             :       {
    1471       37091 :         split_0_1(q, bit, F, G);
    1472       37091 :         *F = RgX_recip_shallow(*F);
    1473       37091 :         *G = RgX_recip_shallow(*G);
    1474             :       }
    1475             :       else
    1476        4323 :         split_2(p,bit,NULL, 1.2837,F,G);
    1477             :     }
    1478             :   }
    1479       55211 : }
    1480             : 
    1481             : /********************************************************************/
    1482             : /**                                                                **/
    1483             : /**                ERROR ESTIMATE FOR THE ROOTS                    **/
    1484             : /**                                                                **/
    1485             : /********************************************************************/
    1486             : 
    1487             : static GEN
    1488      209190 : root_error(long n, long k, GEN roots_pol, long err, GEN shatzle)
    1489             : {
    1490      209190 :   GEN rho, d, eps, epsbis, eps2, aux, rap = NULL;
    1491             :   long i, j;
    1492             : 
    1493      209190 :   d = cgetg(n+1,t_VEC);
    1494     2279840 :   for (i=1; i<=n; i++)
    1495             :   {
    1496     2070650 :     if (i!=k)
    1497             :     {
    1498     1861460 :       aux = gsub(gel(roots_pol,i), gel(roots_pol,k));
    1499     1861460 :       gel(d,i) = gabs(mygprec(aux,31), DEFAULTPREC);
    1500             :     }
    1501             :   }
    1502      209190 :   rho = gabs(mygprec(gel(roots_pol,k),31), DEFAULTPREC);
    1503      209190 :   if (expo(rho) < 0) rho = real_1(DEFAULTPREC);
    1504      209190 :   eps = mulrr(rho, shatzle);
    1505      209190 :   aux = shiftr(powru(rho,n), err);
    1506             : 
    1507      648708 :   for (j=1; j<=2 || (j<=5 && cmprr(rap, dbltor(1.2)) > 0); j++)
    1508             :   {
    1509      439518 :     GEN prod = NULL; /* 1. */
    1510      439518 :     long m = n;
    1511      439518 :     epsbis = mulrr(eps, dbltor(1.25));
    1512     5383982 :     for (i=1; i<=n; i++)
    1513             :     {
    1514     4944464 :       if (i != k && cmprr(gel(d,i),epsbis) > 0)
    1515             :       {
    1516     4474441 :         GEN dif = subrr(gel(d,i),eps);
    1517     4474441 :         prod = prod? mulrr(prod, dif): dif;
    1518     4474441 :         m--;
    1519             :       }
    1520             :     }
    1521      439518 :     eps2 = prod? divrr(aux, prod): aux;
    1522      439518 :     if (m > 1) eps2 = sqrtnr(shiftr(eps2, 2*m-2), m);
    1523      439518 :     rap = divrr(eps,eps2); eps = eps2;
    1524             :   }
    1525      209190 :   return eps;
    1526             : }
    1527             : 
    1528             : /* round a complex or real number x to an absolute value of 2^(-bit) */
    1529             : static GEN
    1530      518555 : mygprec_absolute(GEN x, long bit)
    1531             : {
    1532             :   long e;
    1533             :   GEN y;
    1534             : 
    1535      518555 :   switch(typ(x))
    1536             :   {
    1537             :     case t_REAL:
    1538      340019 :       e = expo(x) + bit;
    1539      340019 :       return (e <= 0 || !signe(x))? real_0_bit(-bit): rtor(x, nbits2prec(e));
    1540             :     case t_COMPLEX:
    1541      156660 :       if (gexpo(gel(x,2)) < -bit) return mygprec_absolute(gel(x,1),bit);
    1542      152705 :       y = cgetg(3,t_COMPLEX);
    1543      152705 :       gel(y,1) = mygprec_absolute(gel(x,1),bit);
    1544      152705 :       gel(y,2) = mygprec_absolute(gel(x,2),bit);
    1545      152705 :       return y;
    1546       21876 :     default: return x;
    1547             :   }
    1548             : }
    1549             : 
    1550             : static long
    1551       49175 : a_posteriori_errors(GEN p, GEN roots_pol, long err)
    1552             : {
    1553       49175 :   long i, n = degpol(p), e_max = -(long)EXPOBITS;
    1554             :   GEN sigma, shatzle;
    1555             : 
    1556       49175 :   err += (long)log2((double)n) + 1;
    1557       49175 :   if (err > -2) return 0;
    1558       49175 :   sigma = real2n(-err, LOWDEFAULTPREC);
    1559             :   /*  2 / ((s - 1)^(1/n) - 1) */
    1560       49175 :   shatzle = divur(2, subrs(sqrtnr(subrs(sigma,1),n), 1));
    1561      258365 :   for (i=1; i<=n; i++)
    1562             :   {
    1563      209190 :     pari_sp av = avma;
    1564      209190 :     GEN x = root_error(n,i,roots_pol,err,shatzle);
    1565      209190 :     long e = gexpo(x);
    1566      209190 :     set_avma(av); if (e > e_max) e_max = e;
    1567      209190 :     gel(roots_pol,i) = mygprec_absolute(gel(roots_pol,i), -e);
    1568             :   }
    1569       49175 :   return e_max;
    1570             : }
    1571             : 
    1572             : /********************************************************************/
    1573             : /**                                                                **/
    1574             : /**                           MAIN                                 **/
    1575             : /**                                                                **/
    1576             : /********************************************************************/
    1577             : static GEN
    1578      169284 : append_clone(GEN r, GEN a) { a = gclone(a); vectrunc_append(r, a); return a; }
    1579             : 
    1580             : /* put roots in placeholder roots_pol so that |P - L_1...L_n| < 2^(-bit)|P|
    1581             :  * returns prod (x-roots_pol[i]) */
    1582             : static GEN
    1583      159597 : split_complete(GEN p, long bit, GEN roots_pol)
    1584             : {
    1585      159597 :   long n = degpol(p);
    1586             :   pari_sp ltop;
    1587             :   GEN p1, F, G, a, b, m1, m2;
    1588             : 
    1589      159597 :   if (n == 1)
    1590             :   {
    1591       39488 :     a = gneg_i(gdiv(gel(p,2), gel(p,3)));
    1592       39488 :     (void)append_clone(roots_pol,a); return p;
    1593             :   }
    1594      120109 :   ltop = avma;
    1595      120109 :   if (n == 2)
    1596             :   {
    1597       64898 :     F = gsub(gsqr(gel(p,3)), gmul2n(gmul(gel(p,2),gel(p,4)), 2));
    1598       64898 :     F = gsqrt(F, nbits2prec(bit));
    1599       64898 :     p1 = ginv(gmul2n(gel(p,4),1));
    1600       64898 :     a = gneg_i(gmul(gadd(F,gel(p,3)), p1));
    1601       64898 :     b =        gmul(gsub(F,gel(p,3)), p1);
    1602       64898 :     a = append_clone(roots_pol,a);
    1603       64898 :     b = append_clone(roots_pol,b); set_avma(ltop);
    1604       64898 :     a = mygprec(a, 3*bit);
    1605       64898 :     b = mygprec(b, 3*bit);
    1606       64898 :     return gmul(gel(p,4), mkpoln(3, gen_1, gneg(gadd(a,b)), gmul(a,b)));
    1607             :   }
    1608       55211 :   split_0(p,bit,&F,&G);
    1609       55211 :   m1 = split_complete(F,bit,roots_pol);
    1610       55211 :   m2 = split_complete(G,bit,roots_pol);
    1611       55211 :   return gerepileupto(ltop, gmul(m1,m2));
    1612             : }
    1613             : 
    1614             : static GEN
    1615      770002 : quicktofp(GEN x)
    1616             : {
    1617      770002 :   const long prec = DEFAULTPREC;
    1618      770002 :   switch(typ(x))
    1619             :   {
    1620      750801 :     case t_INT: return itor(x, prec);
    1621        8034 :     case t_REAL: return rtor(x, prec);
    1622           0 :     case t_FRAC: return fractor(x, prec);
    1623             :     case t_COMPLEX: {
    1624       11167 :       GEN a = gel(x,1), b = gel(x,2);
    1625             :       /* avoid problem with 0, e.g. x = 0 + I*1e-100. We don't want |x| = 0. */
    1626       11167 :       if (isintzero(a)) return cxcompotor(b, prec);
    1627       11125 :       if (isintzero(b)) return cxcompotor(a, prec);
    1628       11125 :       a = cxcompotor(a, prec);
    1629       11125 :       b = cxcompotor(b, prec); return sqrtr(addrr(sqrr(a), sqrr(b)));
    1630             :     }
    1631           0 :     default: pari_err_TYPE("quicktofp",x);
    1632             :       return NULL;/*LCOV_EXCL_LINE*/
    1633             :   }
    1634             : 
    1635             : }
    1636             : 
    1637             : /* bound log_2 |largest root of p| (Fujiwara's bound) */
    1638             : double
    1639      198571 : fujiwara_bound(GEN p)
    1640             : {
    1641      198571 :   pari_sp av = avma;
    1642      198571 :   long i, n = degpol(p);
    1643             :   GEN cc;
    1644             :   double loglc, Lmax;
    1645             : 
    1646      198571 :   if (n <= 0) pari_err_CONSTPOL("fujiwara_bound");
    1647      198571 :   loglc = mydbllog2r( quicktofp(gel(p,n+2)) ); /* log_2 |lc(p)| */
    1648      198571 :   cc = gel(p, 2);
    1649      198571 :   if (gequal0(cc))
    1650       24996 :     Lmax = -pariINFINITY-1;
    1651             :   else
    1652      173575 :     Lmax = (mydbllog2r(quicktofp(cc)) - loglc - 1) / n;
    1653      785776 :   for (i = 1; i < n; i++)
    1654             :   {
    1655      587205 :     GEN y = gel(p,i+2);
    1656             :     double L;
    1657      587205 :     if (gequal0(y)) continue;
    1658      397856 :     L = (mydbllog2r(quicktofp(y)) - loglc) / (n-i);
    1659      397856 :     if (L > Lmax) Lmax = L;
    1660             :   }
    1661      198571 :   return gc_double(av, Lmax+1);
    1662             : }
    1663             : 
    1664             : /* Fujiwara's bound, real roots. Based on the following remark: if
    1665             :  *   p = x^n + sum a_i x^i and q = x^n + sum min(a_i,0)x^i
    1666             :  * then for all x >= 0, p(x) >= q(x). Thus any bound for the (positive) roots
    1667             :  * of q is a bound for the positive roots of p. */
    1668             : double
    1669       72077 : fujiwara_bound_real(GEN p, long sign)
    1670             : {
    1671       72077 :   pari_sp av = avma;
    1672             :   GEN x;
    1673       72077 :   long n = degpol(p), i, signodd, signeven;
    1674       72077 :   if (n <= 0) pari_err_CONSTPOL("fujiwara_bound");
    1675       72077 :   x = shallowcopy(p);
    1676       72077 :   if (gsigne(gel(x, n+2)) > 0)
    1677       72077 :   { signeven = 1; signodd = sign; }
    1678             :   else
    1679           0 :   { signeven = -1; signodd = -sign; }
    1680      353569 :   for (i = 0; i < n; i++)
    1681             :   {
    1682      281492 :     if ((n - i) % 2)
    1683      143591 :     { if (gsigne(gel(x, i+2)) == signodd ) gel(x, i+2) = gen_0; }
    1684             :     else
    1685      137901 :     { if (gsigne(gel(x, i+2)) == signeven) gel(x, i+2) = gen_0; }
    1686             :   }
    1687       72077 :   return gc_double(av, fujiwara_bound(x));
    1688             : }
    1689             : 
    1690             : static GEN
    1691      242095 : mygprecrc_special(GEN x, long prec, long e)
    1692             : {
    1693             :   GEN y;
    1694      242095 :   switch(typ(x))
    1695             :   {
    1696             :     case t_REAL:
    1697       32245 :       if (!signe(x)) return real_0_bit(minss(e, expo(x)));
    1698       29852 :       return (prec > realprec(x))? rtor(x, prec): x;
    1699             :     case t_COMPLEX:
    1700       11818 :       y = cgetg(3,t_COMPLEX);
    1701       11818 :       gel(y,1) = mygprecrc_special(gel(x,1),prec,e);
    1702       11818 :       gel(y,2) = mygprecrc_special(gel(x,2),prec,e);
    1703       11818 :       return y;
    1704      198032 :     default: return x;
    1705             :   }
    1706             : }
    1707             : 
    1708             : /* like mygprec but keep at least the same precision as before */
    1709             : static GEN
    1710       49175 : mygprec_special(GEN x, long bit)
    1711             : {
    1712             :   long lx, i, e, prec;
    1713             :   GEN y;
    1714             : 
    1715       49175 :   if (bit < 0) bit = 0; /* should not happen */
    1716       49175 :   e = gexpo(x) - bit;
    1717       49175 :   prec = nbits2prec(bit);
    1718       49175 :   switch(typ(x))
    1719             :   {
    1720             :     case t_POL:
    1721       49175 :       y = cgetg_copy(x, &lx); y[1] = x[1];
    1722       49175 :       for (i=2; i<lx; i++) gel(y,i) = mygprecrc_special(gel(x,i),prec,e);
    1723       49175 :       break;
    1724             : 
    1725           0 :     default: y = mygprecrc_special(x,prec,e);
    1726             :   }
    1727       49175 :   return y;
    1728             : }
    1729             : 
    1730             : static GEN
    1731       33692 : fix_roots1(GEN r)
    1732             : {
    1733       33692 :   long i, l = lg(r);
    1734       33692 :   GEN allr = cgetg(l, t_VEC);
    1735      172215 :   for (i=1; i<l; i++)
    1736             :   {
    1737      138523 :     GEN t = gel(r,i);
    1738      138523 :     gel(allr,i) = gcopy(t);
    1739      138523 :     gunclone(t);
    1740             :   }
    1741       33692 :   return allr;
    1742             : }
    1743             : static GEN
    1744       49175 : fix_roots(GEN r, long h, long bit)
    1745             : {
    1746             :   long i, j, k, l, prec;
    1747             :   GEN allr, ro1;
    1748       49175 :   if (h == 1) return fix_roots1(r);
    1749       15483 :   prec = nbits2prec(bit);
    1750       15483 :   ro1 = grootsof1(h, prec) + 1;
    1751       15483 :   l = lg(r)-1;
    1752       15483 :   allr = cgetg(h*l+1, t_VEC);
    1753       46244 :   for (k=1,i=1; i<=l; i++)
    1754             :   {
    1755       30761 :     GEN p2, p1 = gel(r,i);
    1756       30761 :     p2 = (h == 2)? gsqrt(p1, prec): gsqrtn(p1, utoipos(h), NULL, prec);
    1757       30761 :     for (j=0; j<h; j++) gel(allr,k++) = gmul(p2, gel(ro1,j));
    1758       30761 :     gunclone(p1);
    1759             :   }
    1760       15483 :   return allr;
    1761             : }
    1762             : 
    1763             : static GEN
    1764       48742 : all_roots(GEN p, long bit)
    1765             : {
    1766             :   GEN pd, q, roots_pol, m;
    1767       48742 :   long bit2, i, e, h, elc, n = degpol(p);
    1768             :   double fb;
    1769             :   pari_sp av;
    1770             : 
    1771       48742 :   pd = RgX_deflate_max(p, &h); elc = gexpo(leading_coeff(pd));
    1772       48742 :   fb = fujiwara_bound(pd);
    1773       48742 :   e = (fb < 0)? 0: (long)(2 * fb);
    1774       48742 :   bit2 = bit + gexpo(pd) + (long)log2(n/h)+1+e;
    1775       48742 :   e = 0;
    1776       49175 :   for (av=avma,i=1;; i++,set_avma(av))
    1777             :   {
    1778       49608 :     roots_pol = vectrunc_init(n+1);
    1779       49175 :     bit2 += e + (n << i);
    1780       49175 :     q = RgX_gtofp_bit(mygprec(pd,bit2), bit2);
    1781       49175 :     q[1] = evalsigne(1)|evalvarn(0);
    1782       49175 :     m = split_complete(q,bit2,roots_pol);
    1783       49175 :     roots_pol = fix_roots(roots_pol, h, bit2);
    1784       49175 :     q = mygprec_special(pd,bit2);
    1785       49175 :     q[1] = evalsigne(1)|evalvarn(0);
    1786       49175 :     e = gexpo(RgX_sub(q, m)) - elc + (long)log2((double)n) + 1;
    1787       49175 :     if (e < -2*bit2) e = -2*bit2; /* avoid e = -oo */
    1788       49175 :     if (e < 0)
    1789             :     {
    1790       49175 :       e = bit + a_posteriori_errors(p,roots_pol,e);
    1791       97917 :       if (e < 0) return roots_pol;
    1792             :     }
    1793         433 :     if (DEBUGLEVEL > 7)
    1794           0 :       err_printf("all_roots: restarting, i = %ld, e = %ld\n", i,e);
    1795             :   }
    1796             : }
    1797             : 
    1798             : INLINE int
    1799       20254 : isexactscalar(GEN x) { long tx = typ(x); return is_rational_t(tx); }
    1800             : 
    1801             : static int
    1802        9765 : isexactpol(GEN p)
    1803             : {
    1804        9765 :   long i,n = degpol(p);
    1805       22469 :   for (i=0; i<=n; i++)
    1806       20254 :     if (!isexactscalar(gel(p,i+2))) return 0;
    1807        2215 :   return 1;
    1808             : }
    1809             : 
    1810             : static GEN
    1811        2215 : solve_exact_pol(GEN p, long bit)
    1812             : {
    1813        2215 :   long i, j, k, m, n = degpol(p), iroots = 0;
    1814        2215 :   GEN ex, factors, v = zerovec(n);
    1815             : 
    1816        2215 :   factors = ZX_squff(Q_primpart(p), &ex);
    1817        4430 :   for (i=1; i<lg(factors); i++)
    1818             :   {
    1819        2215 :     GEN roots_fact = all_roots(gel(factors,i), bit);
    1820        2215 :     n = degpol(gel(factors,i));
    1821        2215 :     m = ex[i];
    1822       12046 :     for (j=1; j<=n; j++)
    1823        9831 :       for (k=1; k<=m; k++) v[++iroots] = roots_fact[j];
    1824             :   }
    1825        2215 :   return v;
    1826             : }
    1827             : 
    1828             : /* return the roots of p with absolute error bit */
    1829             : static GEN
    1830        9765 : roots_com(GEN q, long bit)
    1831             : {
    1832             :   GEN L, p;
    1833        9765 :   long v = RgX_valrem_inexact(q, &p);
    1834        9765 :   int ex = isexactpol(p);
    1835        9765 :   if (!ex) p = RgX_normalize1(p);
    1836        9765 :   if (lg(p) == 3)
    1837           0 :     L = cgetg(1,t_VEC); /* constant polynomial */
    1838             :   else
    1839        9765 :     L = ex? solve_exact_pol(p,bit): all_roots(p,bit);
    1840        9765 :   if (v)
    1841             :   {
    1842         182 :     GEN M, z, t = gel(q,2);
    1843             :     long i, x, y, l, n;
    1844             : 
    1845         182 :     if (isrationalzero(t)) x = -bit;
    1846             :     else
    1847             :     {
    1848           7 :       n = gexpo(t);
    1849           7 :       x = n / v; l = degpol(q);
    1850          35 :       for (i = v; i <= l; i++)
    1851             :       {
    1852          28 :         t  = gel(q,i+2);
    1853          28 :         if (isrationalzero(t)) continue;
    1854          28 :         y = (n - gexpo(t)) / i;
    1855          28 :         if (y < x) x = y;
    1856             :       }
    1857             :     }
    1858         182 :     z = real_0_bit(x); l = v + lg(L);
    1859         182 :     M = cgetg(l, t_VEC); L -= v;
    1860         182 :     for (i = 1; i <= v; i++) gel(M,i) = z;
    1861         182 :     for (     ; i <  l; i++) gel(M,i) = gel(L,i);
    1862         182 :     L = M;
    1863             :   }
    1864        9765 :   return L;
    1865             : }
    1866             : 
    1867             : static GEN
    1868      163202 : tocomplex(GEN x, long l, long bit)
    1869             : {
    1870             :   GEN y;
    1871      163202 :   if (typ(x) == t_COMPLEX)
    1872             :   {
    1873      151507 :     if (signe(gel(x,1))) return mygprecrc(x, l, -bit);
    1874       22531 :     x = gel(x,2);
    1875       22531 :     y = cgetg(3,t_COMPLEX);
    1876       22531 :     gel(y,1) = real_0_bit(-bit);
    1877       22531 :     gel(y,2) = mygprecrc(x, l, -bit);
    1878             :   }
    1879             :   else
    1880             :   {
    1881       11695 :     y = cgetg(3,t_COMPLEX);
    1882       11695 :     gel(y,1) = mygprecrc(x, l, -bit);
    1883       11695 :     gel(y,2) = real_0_bit(-bit);
    1884             :   }
    1885       34226 :   return y;
    1886             : }
    1887             : 
    1888             : /* x,y are t_COMPLEX of t_REALs or t_REAL, compare wrt |Im x| - |Im y|,
    1889             :  * then Re x - Re y, up to 2^-e absolute error */
    1890             : static int
    1891      332039 : cmp_complex_appr(void *E, GEN x, GEN y)
    1892             : {
    1893      332039 :   long e = (long)E;
    1894             :   GEN z, xi, yi, xr, yr;
    1895             :   long sz, sxi, syi;
    1896      332039 :   if (typ(x) == t_COMPLEX) { xr = gel(x,1); xi = gel(x,2); sxi = signe(xi); }
    1897       94498 :   else { xr = x; xi = NULL; sxi = 0; }
    1898      332039 :   if (typ(y) == t_COMPLEX) { yr = gel(y,1); yi = gel(y,2); syi = signe(yi); }
    1899       84997 :   else { yr = y; yi = NULL; syi = 0; }
    1900             :   /* Compare absolute values of imaginary parts */
    1901      332039 :   if (!sxi)
    1902             :   {
    1903      108886 :     if (syi && expo(yi) >= e) return -1;
    1904             :     /* |Im x| ~ |Im y| ~ 0 */
    1905             :   }
    1906      223153 :   else if (!syi)
    1907             :   {
    1908        5073 :     if (sxi && expo(xi) >= e) return 1;
    1909             :     /* |Im x| ~ |Im y| ~ 0 */
    1910             :   }
    1911             :   else
    1912             :   {
    1913      218080 :     z = addrr_sign(xi, 1, yi, -1); sz = signe(z);
    1914      218080 :     if (sz && expo(z) >= e) return (int)sz;
    1915             :   }
    1916             :   /* |Im x| ~ |Im y|, sort according to real parts */
    1917      186984 :   z = subrr(xr, yr); sz = signe(z);
    1918      186984 :   if (sz && expo(z) >= e) return (int)sz;
    1919             :   /* Re x ~ Re y. Place negative imaginary part before positive */
    1920       70841 :   return (int) (sxi - syi);
    1921             : }
    1922             : 
    1923             : static GEN
    1924       48819 : clean_roots(GEN L, long l, long bit, long clean)
    1925             : {
    1926       48819 :   long i, n = lg(L), ex = 5 - bit;
    1927       48819 :   GEN res = cgetg(n,t_COL);
    1928      256145 :   for (i=1; i<n; i++)
    1929             :   {
    1930      207326 :     GEN c = gel(L,i);
    1931      207326 :     if (clean && isrealappr(c,ex))
    1932             :     {
    1933       44124 :       if (typ(c) == t_COMPLEX) c = gel(c,1);
    1934       44124 :       c = mygprecrc(c, l, -bit);
    1935             :     }
    1936             :     else
    1937      163202 :       c = tocomplex(c, l, bit);
    1938      207326 :     gel(res,i) = c;
    1939             :   }
    1940       48819 :   gen_sort_inplace(res, (void*)ex, &cmp_complex_appr, NULL);
    1941       48819 :   return res;
    1942             : }
    1943             : 
    1944             : /* the vector of roots of p, with absolute error 2^(- prec2nbits(l)) */
    1945             : static GEN
    1946        9793 : roots_aux(GEN p, long l, long clean)
    1947             : {
    1948        9793 :   pari_sp av = avma;
    1949             :   long bit;
    1950             :   GEN L;
    1951             : 
    1952        9793 :   if (typ(p) != t_POL)
    1953             :   {
    1954          21 :     if (gequal0(p)) pari_err_ROOTS0("roots");
    1955          14 :     if (!isvalidcoeff(p)) pari_err_TYPE("roots",p);
    1956           7 :     return cgetg(1,t_COL); /* constant polynomial */
    1957             :   }
    1958        9772 :   if (!signe(p)) pari_err_ROOTS0("roots");
    1959        9772 :   checkvalidpol(p,"roots");
    1960        9765 :   if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
    1961        9765 :   if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
    1962        9765 :   bit = prec2nbits(l);
    1963        9765 :   L = roots_com(p, bit);
    1964        9765 :   return gerepilecopy(av, clean_roots(L, l, bit, clean));
    1965             : }
    1966             : GEN
    1967        9556 : roots(GEN p, long l) { return roots_aux(p,l, 0); }
    1968             : /* clean up roots. If root is real replace it by its real part */
    1969             : GEN
    1970         237 : cleanroots(GEN p, long l) { return roots_aux(p,l, 1); }
    1971             : 
    1972             : /* private variant of conjvec. Allow non rational coefficients, shallow
    1973             :  * function. */
    1974             : GEN
    1975          84 : polmod_to_embed(GEN x, long prec)
    1976             : {
    1977          84 :   GEN v, T = gel(x,1), A = gel(x,2);
    1978             :   long i, l;
    1979          84 :   if (typ(A) != t_POL || varn(A) != varn(T))
    1980             :   {
    1981           7 :     checkvalidpol(T,"polmod_to_embed");
    1982           7 :     return const_col(degpol(T), A);
    1983             :   }
    1984          77 :   v = cleanroots(T,prec); l = lg(v);
    1985          77 :   for (i=1; i<l; i++) gel(v,i) = poleval(A,gel(v,i));
    1986          77 :   return v;
    1987             : }
    1988             : 
    1989             : GEN
    1990       39054 : QX_complex_roots(GEN p, long l)
    1991             : {
    1992       39054 :   pari_sp av = avma;
    1993             :   long bit, v;
    1994             :   GEN L;
    1995             : 
    1996       39054 :   if (!signe(p)) pari_err_ROOTS0("QX_complex_roots");
    1997       39054 :   if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
    1998       39054 :   if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
    1999       39054 :   bit = prec2nbits(l);
    2000       39054 :   v = RgX_valrem(p, &p);
    2001       39054 :   L = lg(p) > 3? all_roots(Q_primpart(p), bit): cgetg(1,t_COL);
    2002       39054 :   if (v) L = shallowconcat(const_vec(v, real_0_bit(-bit)), L);
    2003       39054 :   return gerepilecopy(av, clean_roots(L, l, bit, 1));
    2004             : }
    2005             : 
    2006             : /********************************************************************/
    2007             : /**                                                                **/
    2008             : /**                REAL ROOTS OF INTEGER POLYNOMIAL                **/
    2009             : /**                                                                **/
    2010             : /********************************************************************/
    2011             : 
    2012             : /* Count sign changes in the coefficients of (x+1)^deg(P)*P(1/(x+1)), P
    2013             :  * has no rational root. The inversion is implicit (we take coefficients
    2014             :  * backwards). */
    2015             : static long
    2016      514223 : X2XP1(GEN P, GEN *Premapped)
    2017             : {
    2018      514223 :   const pari_sp av = avma;
    2019      514223 :   GEN v = shallowcopy(P);
    2020      514223 :   long i, j, nb, s, dP = degpol(P), vlim = dP+2;
    2021             : 
    2022      514223 :   for (j = 2; j < vlim; j++) gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
    2023      514223 :   s = -signe(gel(v, vlim));
    2024      514223 :   vlim--; nb = 0;
    2025     1890680 :   for (i = 1; i < dP; i++)
    2026             :   {
    2027     1774804 :     long s2 = -signe(gel(v, 2));
    2028     1774804 :     int flag = (s2 == s);
    2029    21011398 :     for (j = 2; j < vlim; j++)
    2030             :     {
    2031    19236594 :       gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
    2032    19236594 :       if (flag) flag = (s2 != signe(gel(v, j+1)));
    2033             :     }
    2034     1774804 :     if (s == signe(gel(v, vlim)))
    2035             :     {
    2036      485027 :       if (++nb >= 2) return gc_long(av,2);
    2037      299868 :       s = -s;
    2038             :     }
    2039             :     /* if flag is set there will be no further sign changes */
    2040     1589645 :     if (flag && (!Premapped || !nb)) return gc_long(av, nb);
    2041     1376457 :     vlim--;
    2042     1376457 :     if (gc_needed(av, 3))
    2043             :     {
    2044           0 :       if (DEBUGMEM>1) pari_warn(warnmem, "X2XP1, i = %ld/%ld", i, dP-1);
    2045           0 :       if (!Premapped) setlg(v, vlim + 2);
    2046           0 :       v = gerepilecopy(av, v);
    2047             :     }
    2048             :   }
    2049      115876 :   if (vlim >= 2 && s == signe(gel(v, vlim))) nb++;
    2050      115876 :   if (Premapped && nb == 1) *Premapped = v; else set_avma(av);
    2051      115876 :   return nb;
    2052             : }
    2053             : 
    2054             : static long
    2055           0 : _intervalcmp(GEN x, GEN y)
    2056             : {
    2057           0 :   if (typ(x) == t_VEC) x = gel(x, 1);
    2058           0 :   if (typ(y) == t_VEC) y = gel(y, 1);
    2059           0 :   return gcmp(x, y);
    2060             : }
    2061             : 
    2062             : static GEN
    2063     5225567 : _gen_nored(void *E, GEN x) { (void)E; return x; }
    2064             : static GEN
    2065    16812157 : _mp_add(void *E, GEN x, GEN y) { (void)E; return mpadd(x, y); }
    2066             : static GEN
    2067           0 : _mp_sub(void *E, GEN x, GEN y) { (void)E; return mpsub(x, y); }
    2068             : static GEN
    2069     3053027 : _mp_mul(void *E, GEN x, GEN y) { (void)E; return mpmul(x, y); }
    2070             : static GEN
    2071     3976589 : _mp_sqr(void *E, GEN x) { (void)E; return mpsqr(x); }
    2072             : static GEN
    2073     6717447 : _gen_one(void *E) { (void)E; return gen_1; }
    2074             : static GEN
    2075       36437 : _gen_zero(void *E) { (void)E; return gen_0; }
    2076             : 
    2077             : static struct bb_algebra mp_algebra = { _gen_nored, _mp_add, _mp_sub,
    2078             :                          _mp_mul, _mp_sqr, _gen_one, _gen_zero };
    2079             : 
    2080             : static GEN
    2081    21316036 : _mp_cmul(void *E, GEN P, long a, GEN x) {(void)E; return mpmul(gel(P,a+2), x);}
    2082             : 
    2083             : /* Split the polynom P in two parts, whose coeffs have constant sign:
    2084             :  * P(X) = X^D*Pp + Pm. Also compute the two parts of the derivative of P,
    2085             :  * Pprimem = Pm', Pprimep = X*Pp'+ D*Pp => P' = X^(D-1)*Pprimep + Pprimem;
    2086             :  * Pprimep[i] = (i+D) Pp[i]. Return D */
    2087             : static long
    2088       69531 : split_pols(GEN P, GEN *pPp, GEN *pPm, GEN *pPprimep, GEN *pPprimem)
    2089             : {
    2090       69531 :   long i, D, dP = degpol(P), s0 = signe(gel(P,2));
    2091             :   GEN Pp, Pm, Pprimep, Pprimem;
    2092      282790 :   for(i=1; i <= dP; i++)
    2093      282790 :     if (signe(gel(P, i+2)) == -s0) break;
    2094       69531 :   D = i;
    2095       69531 :   Pm = cgetg(D + 2, t_POL);
    2096       69531 :   Pprimem = cgetg(D + 1, t_POL);
    2097       69531 :   Pp = cgetg(dP-D + 3, t_POL);
    2098       69531 :   Pprimep = cgetg(dP-D + 3, t_POL);
    2099       69531 :   Pm[1] = Pp[1] = Pprimem[1] = Pprimep[1] = P[1];
    2100      352321 :   for(i=0; i < D; i++)
    2101             :   {
    2102      282790 :     GEN c = gel(P, i+2);
    2103      282790 :     gel(Pm, i+2) = c;
    2104      282790 :     if (i) gel(Pprimem, i+1) = mului(i, c);
    2105             :   }
    2106      380880 :   for(; i <= dP; i++)
    2107             :   {
    2108      311349 :     GEN c = gel(P, i+2);
    2109      311349 :     gel(Pp, i+2-D) = c;
    2110      311349 :     gel(Pprimep, i+2-D) = mului(i, c);
    2111             :   }
    2112       69531 :   *pPm = normalizepol_lg(Pm, D+2);
    2113       69531 :   *pPprimem = normalizepol_lg(Pprimem, D+1);
    2114       69531 :   *pPp = normalizepol_lg(Pp, dP-D+3);
    2115       69531 :   *pPprimep = normalizepol_lg(Pprimep, dP-D+3);
    2116       69531 :   return dP - degpol(*pPp);
    2117             : }
    2118             : 
    2119             : static GEN
    2120     2270158 : bkeval_single_power(long d, GEN V)
    2121             : {
    2122     2270158 :   long mp = lg(V) - 2;
    2123     2270158 :   if (d > mp) return gmul(gpowgs(gel(V, mp+1), d/mp), gel(V, (d%mp)+1));
    2124     2270158 :   return gel(V, d+1);
    2125             : }
    2126             : 
    2127             : static GEN
    2128     2270158 : splitpoleval(GEN Pp, GEN Pm, GEN pows, long D, long bitprec)
    2129             : {
    2130     2270158 :   GEN vp = gen_bkeval_powers(Pp, degpol(Pp), pows, NULL, &mp_algebra, _mp_cmul);
    2131     2270158 :   GEN vm = gen_bkeval_powers(Pm, degpol(Pm), pows, NULL, &mp_algebra, _mp_cmul);
    2132     2270158 :   GEN xa = bkeval_single_power(D, pows);
    2133             :   GEN r;
    2134     2270158 :   if (!signe(vp)) return vm;
    2135     2270158 :   vp = gmul(vp, xa);
    2136     2270158 :   r = gadd(vp, vm);
    2137     2270158 :   if (gexpo(vp) - (signe(r)? gexpo(r): 0) > prec2nbits(realprec(vp)) - bitprec)
    2138       51919 :     return NULL;
    2139     2218239 :   return r;
    2140             : }
    2141             : 
    2142             : /* optimized Cauchy bound for P = X^D*Pp + Pm, D > deg(Pm) */
    2143             : static GEN
    2144       69531 : splitcauchy(GEN Pp, GEN Pm, long prec)
    2145             : {
    2146       69531 :   GEN S = gel(Pp,2), A = gel(Pm,2);
    2147       69531 :   long i, lPm = lg(Pm), lPp = lg(Pp);
    2148       69531 :   for (i=3; i < lPm; i++) { GEN c = gel(Pm,i); if (abscmpii(A, c) < 0) A = c; }
    2149       69531 :   for (i=3; i < lPp; i++) S = addii(S, gel(Pp, i));
    2150       69531 :   return subsr(1, rdivii(A, S, prec)); /* 1 + |Pm|_oo / |Pp|_1 */
    2151             : }
    2152             : 
    2153             : static GEN
    2154        3746 : ZX_deg1root(GEN P, long prec)
    2155             : {
    2156        3746 :   GEN a = gel(P,3), b = gel(P,2);
    2157        3746 :   if (is_pm1(a))
    2158             :   {
    2159        3746 :     b = itor(b, prec); if (signe(a) > 0) togglesign(b);
    2160        3746 :     return b;
    2161             :   }
    2162           0 :   return rdivii(negi(b), a, prec);
    2163             : }
    2164             : 
    2165             : /* Newton for polynom P, P(0)!=0, with unique sign change => one root in ]0,oo[
    2166             :  * P' has also at most one zero there */
    2167             : static GEN
    2168       69531 : polsolve(GEN P, long bitprec)
    2169             : {
    2170             :   pari_sp av;
    2171             :   GEN Pp, Pm, Pprimep, Pprimem, Pprime, Pprime2, ra, rb, rc, Pc;
    2172       69531 :   long dP = degpol(P), prec = nbits2prec(bitprec);
    2173             :   long expoold, iter, D, rt, s0, bitaddprec, cprec, PREC;
    2174             : 
    2175       69531 :   if (dP == 1) return ZX_deg1root(P, prec);
    2176       69531 :   Pprime = ZX_deriv(P);
    2177       69531 :   Pprime2 = ZX_deriv(Pprime);
    2178       69531 :   bitaddprec = 1 + 2*expu(dP); PREC = prec + nbits2prec(bitaddprec);
    2179       69531 :   D = split_pols(P, &Pp, &Pm, &Pprimep, &Pprimem); /* P = X^D*Pp + Pm */
    2180       69531 :   s0 = signe(gel(P, 2));
    2181       69531 :   rt = maxss(D, brent_kung_optpow(maxss(degpol(Pp), degpol(Pm)), 2, 1));
    2182       69531 :   rb = splitcauchy(Pp, Pm, DEFAULTPREC);
    2183       69531 :   for (cprec = DEFAULTPREC, expoold = LONG_MAX;;)
    2184           0 :   {
    2185       69531 :     GEN pows = gen_powers(rb, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
    2186       69531 :     Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec);
    2187       69531 :     if (!Pc) { cprec += EXTRAPRECWORD; rb = rtor(rb, cprec); continue; }
    2188       69531 :     if (signe(Pc) != s0) break;
    2189           0 :     shiftr_inplace(rb,1);
    2190             :   }
    2191       69531 :   for (iter = 0, ra = NULL;;)
    2192     1117724 :   {
    2193             :     GEN wdth;
    2194     1187255 :     iter++;
    2195     1187255 :     if (ra)
    2196      344609 :       rc = shiftr(addrr(ra, rb), -1);
    2197             :     else
    2198      842646 :       rc = shiftr(rb, -1);
    2199             :     for(;;)
    2200           0 :     {
    2201     1187255 :       GEN pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
    2202     1187255 :       Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec+2);
    2203     1187255 :       if (Pc) break;
    2204           0 :       cprec += EXTRAPRECWORD;
    2205           0 :       rc = rtor(rc, cprec);
    2206             :     }
    2207     1187255 :     if (signe(Pc) == s0)
    2208      234778 :       ra = rc;
    2209             :     else
    2210      952477 :       rb = rc;
    2211     1187255 :     if (!ra) continue;
    2212      414140 :     wdth = subrr(rb, ra);
    2213      414140 :     if (!(iter % 8))
    2214             :     {
    2215       69797 :       GEN m1 = poleval(Pprime, ra), M2;
    2216       69797 :       if (signe(m1) == s0) continue;
    2217       68894 :       M2 = poleval(Pprime2, rb);
    2218       68894 :       if (abscmprr(gmul(M2, wdth), shiftr(m1, 1)) > 0) continue;
    2219       66395 :       break;
    2220             :     }
    2221      344343 :     else if (gexpo(wdth) <= -bitprec)
    2222        3136 :       break;
    2223             :   }
    2224       69531 :   rc = rb; av = avma;
    2225      437155 :   for(;; rc = gerepileuptoleaf(av, rc))
    2226      437155 :   {
    2227             :     long exponew;
    2228      506686 :     GEN Ppc, dist, rcold = rc;
    2229      506686 :     GEN pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
    2230      506686 :     Ppc = splitpoleval(Pprimep, Pprimem, pows, D-1, bitaddprec+4);
    2231      506686 :     if (Ppc) Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec+4);
    2232      506686 :     if (!Ppc || !Pc)
    2233             :     {
    2234       51919 :       if (cprec >= PREC)
    2235        3783 :         cprec += EXTRAPRECWORD;
    2236             :       else
    2237       48136 :         cprec = minss(2*cprec, PREC);
    2238       51919 :       rc = rtor(rc, cprec); continue; /* backtrack one step */
    2239             :     }
    2240      454767 :     dist = typ(Ppc) == t_REAL? divrr(Pc, Ppc): divri(Pc, Ppc);
    2241      454767 :     rc = subrr(rc, dist);
    2242      454767 :     if (cmprr(ra, rc) > 0 || cmprr(rb, rc) < 0)
    2243             :     {
    2244           0 :       if (cprec >= PREC) break;
    2245           0 :       cprec = minss(2*cprec, PREC);
    2246           0 :       rc = rtor(rcold, cprec); continue; /* backtrack one step */
    2247             :     }
    2248      454767 :     if (expoold == LONG_MAX) { expoold = expo(dist); continue; }
    2249      381339 :     exponew = expo(dist);
    2250      381339 :     if (exponew < -bitprec - 1)
    2251             :     {
    2252      146057 :       if (cprec >= PREC) break;
    2253       76526 :       cprec = minss(2*cprec, PREC);
    2254       76526 :       rc = rtor(rc, cprec); continue;
    2255             :     }
    2256      235282 :     if (exponew > expoold - 2)
    2257             :     {
    2258        3897 :       if (cprec >= PREC) break;
    2259        3897 :       expoold = LONG_MAX;
    2260        3897 :       cprec = minss(2*cprec, PREC);
    2261        3897 :       rc = rtor(rc, cprec); continue;
    2262             :     }
    2263      231385 :     expoold = exponew;
    2264             :   }
    2265       69531 :   return rtor(rc, prec);
    2266             : }
    2267             : 
    2268             : /* Return primpart(P(x / 2)) */
    2269             : static GEN
    2270      153934 : ZX_rescale2prim(GEN P)
    2271             : {
    2272      153934 :   long i, l = lg(P), v, n;
    2273             :   GEN Q;
    2274      153934 :   if (l==2) return pol_0(varn(P));
    2275      153934 :   Q = cgetg(l,t_POL); v = vali(gel(P,l-1));
    2276      504780 :   for (i = l-2, n = 1; v > n && i >= 2; i--, n++)
    2277      350846 :     v = minss(v, vali(gel(P,i)) + n);
    2278      153934 :   gel(Q,l-1) = v? shifti(gel(P,l-1), -v): gel(P,l-1);
    2279     1239938 :   for (i = l-2, n = 1-v; i >= 2; i--, n++)
    2280     1086004 :     gel(Q,i) = shifti(gel(P,i), n);
    2281      153934 :   Q[1] = P[1]; return Q;
    2282             : }
    2283             : 
    2284             : /* assume Q0 has no rational root */
    2285             : static GEN
    2286       70369 : usp(GEN Q0, long flag, long bitprec)
    2287             : {
    2288       70369 :   const pari_sp av = avma;
    2289             :   GEN Qremapped, Q, c, Lc, Lk, sol;
    2290       70369 :   GEN *pQremapped = flag == 1? &Qremapped: NULL;
    2291       70369 :   const long prec = nbits2prec(bitprec), deg = degpol(Q0);
    2292       70369 :   long listsize = 64, nbr = 0, nb_todo, ind, indf, i, k, nb;
    2293             : 
    2294       70369 :   sol = zerocol(deg);
    2295       70369 :   Lc = zerovec(listsize);
    2296       70369 :   Lk = cgetg(listsize+1, t_VECSMALL);
    2297       70369 :   k = Lk[1] = 0;
    2298       70369 :   ind = 1; indf = 2;
    2299       70369 :   Q = Q0;
    2300       70369 :   c = gen_0;
    2301       70369 :   nb_todo = 1;
    2302      654961 :   while (nb_todo)
    2303             :   {
    2304      514223 :     GEN nc = gel(Lc, ind);
    2305             :     pari_sp av2;
    2306      514223 :     if (Lk[ind] == k + 1)
    2307             :     {
    2308      153934 :       Q = Q0 = ZX_rescale2prim(Q0);
    2309      153934 :       c = gen_0;
    2310             :     }
    2311      514223 :     if (!equalii(nc, c)) Q = ZX_translate(Q, subii(nc, c));
    2312      514223 :     av2 = avma;
    2313      514223 :     k = Lk[ind];
    2314      514223 :     ind++;
    2315      514223 :     c = nc;
    2316      514223 :     nb_todo--;
    2317      514223 :     nb = X2XP1(Q, pQremapped);
    2318             : 
    2319      514223 :     if (nb == 1)
    2320             :     { /* exactly one root */
    2321       86639 :       GEN s = gen_0;
    2322       86639 :       if (flag == 0)
    2323           0 :         s = mkvec2(gmul2n(c,-k), gmul2n(addiu(c,1),-k));
    2324       86639 :       else if (flag == 1) /* Caveat: Qremapped is the reciprocal polynomial */
    2325             :       {
    2326       69531 :         s = polsolve(*pQremapped, bitprec+1);
    2327       69531 :         s = addir(c, divrr(s, addsr(1, s)));
    2328       69531 :         shiftr_inplace(s, -k);
    2329       69531 :         if (realprec(s) != prec) s = rtor(s, prec);
    2330             :       }
    2331       86639 :       gel(sol, ++nbr) = gerepileupto(av2, s);
    2332             :     }
    2333      427584 :     else if (nb)
    2334             :     { /* unknown, add two nodes to refine */
    2335      221927 :       if (indf + 2 > listsize)
    2336             :       {
    2337         429 :         if (ind>1)
    2338             :         {
    2339        2377 :           for (i = ind; i < indf; i++)
    2340             :           {
    2341        1948 :             gel(Lc, i-ind+1) = gel(Lc, i);
    2342        1948 :             Lk[i-ind+1] = Lk[i];
    2343             :           }
    2344         429 :           indf -= ind-1;
    2345         429 :           ind = 1;
    2346             :         }
    2347         429 :         if (indf + 2 > listsize)
    2348             :         {
    2349           0 :           listsize *= 2;
    2350           0 :           Lc = vec_lengthen(Lc, listsize);
    2351           0 :           Lk = vecsmall_lengthen(Lk, listsize);
    2352             :         }
    2353         429 :         for (i = indf; i <= listsize; i++) gel(Lc, i) = gen_0;
    2354             :       }
    2355      221927 :       gel(Lc, indf) = nc = shifti(c, 1);
    2356      221927 :       gel(Lc, indf + 1) = addiu(nc, 1);
    2357      221927 :       Lk[indf] = Lk[indf + 1] = k + 1;
    2358      221927 :       indf += 2;
    2359      221927 :       nb_todo += 2;
    2360             :     }
    2361      514223 :     if (gc_needed(av, 2))
    2362             :     {
    2363           0 :       gerepileall(av, 6, &Q0, &Q, &c, &Lc, &Lk, &sol);
    2364           0 :       if (DEBUGMEM > 1) pari_warn(warnmem, "ZX_Uspensky", avma);
    2365             :     }
    2366             :   }
    2367       70369 :   setlg(sol, nbr+1);
    2368       70369 :   return gerepilecopy(av, sol);
    2369             : }
    2370             : 
    2371             : static GEN
    2372          14 : ZX_Uspensky_equal_yes(GEN a, long flag, long bit)
    2373             : {
    2374          14 :   if (flag == 2) return gen_1;
    2375           7 :   if (flag == 1 && typ(a) != t_REAL)
    2376             :   {
    2377           7 :     if (typ(a) == t_INT && !signe(a))
    2378           0 :       a = real_0_bit(bit);
    2379             :     else
    2380           7 :       a = gtofp(a, nbits2prec(bit));
    2381             :   }
    2382           7 :   return mkcol(a);
    2383             : }
    2384             : static GEN
    2385          14 : ZX_Uspensky_no(long flag)
    2386          14 : { return flag <= 1 ? cgetg(1, t_COL) : gen_0; }
    2387             : /* ZX_Uspensky(P, [a,a], flag) */
    2388             : static GEN
    2389          21 : ZX_Uspensky_equal(GEN P, GEN a, long flag, long bit)
    2390             : {
    2391          21 :   if (typ(a) != t_INFINITY && gequal0(poleval(P, a)))
    2392          14 :     return ZX_Uspensky_equal_yes(a, flag, bit);
    2393             :   else
    2394           7 :     return ZX_Uspensky_no(flag);
    2395             : }
    2396             : static int
    2397        2484 : sol_ok(GEN r, GEN a, GEN b) { return gcmp(a, r) <= 0 && gcmp(r, b) <= 0; }
    2398             : 
    2399             : /* P a ZX without real double roots; better if primitive and squarefree but
    2400             :  * caller should ensure that. If flag & 4 assume that P has no rational root
    2401             :  * (modest speedup) */
    2402             : GEN
    2403       65299 : ZX_Uspensky(GEN P, GEN ab, long flag, long bitprec)
    2404             : {
    2405       65299 :   pari_sp av = avma;
    2406             :   GEN a, b, res, sol;
    2407             :   double fb;
    2408             :   long l, nbz, deg;
    2409             : 
    2410       65299 :   if (ab)
    2411             :   {
    2412       49758 :     if (typ(ab) == t_VEC)
    2413             :     {
    2414       47525 :       if (lg(ab) != 3) pari_err_DIM("ZX_Uspensky");
    2415       47525 :       a = gel(ab, 1);
    2416       47525 :       b = gel(ab, 2);
    2417             :     }
    2418             :     else
    2419             :     {
    2420        2233 :       a = ab;
    2421        2233 :       b = mkoo();
    2422             :     }
    2423             :   }
    2424             :   else
    2425             :   {
    2426       15541 :     a = mkmoo();
    2427       15541 :     b = mkoo();
    2428             :   }
    2429       65299 :   if (flag & 4)
    2430             :   {
    2431       16913 :     if (gcmp(a, b) >= 0) { set_avma(av); return ZX_Uspensky_no(flag); }
    2432       16913 :     flag &= ~4;
    2433       16913 :     sol = cgetg(1, t_COL);
    2434             :   }
    2435             :   else
    2436             :   {
    2437       48386 :     switch (gcmp(a, b))
    2438             :     {
    2439           7 :       case 1: set_avma(av); return ZX_Uspensky_no(flag);
    2440          21 :       case 0: return gerepilecopy(av, ZX_Uspensky_equal(P, a, flag, bitprec));
    2441             :     }
    2442       48358 :     sol = nfrootsQ(P);
    2443             :   }
    2444       65271 :   nbz = 0; l = lg(sol);
    2445       65271 :   if (l > 1)
    2446             :   {
    2447             :     long i, j;
    2448        2015 :     P = RgX_div(P, roots_to_pol(sol, varn(P)));
    2449        2015 :     if (!RgV_is_ZV(sol)) P = Q_primpart(P);
    2450        4499 :     for (i = j = 1; i < l; i++)
    2451        2484 :       if (sol_ok(gel(sol,i), a, b)) gel(sol,j++) = gel(sol,i);
    2452        2015 :     setlg(sol, j);
    2453        2015 :     if (flag == 2) { nbz = j-1; sol = utoi(nbz); }
    2454        1854 :     else if (flag == 1) sol = RgC_gtofp(sol, nbits2prec(bitprec));
    2455             :   }
    2456       65271 :   deg = degpol(P);
    2457       65271 :   if (deg == 0) return gerepilecopy(av, sol);
    2458       63795 :   if (typ(a) == t_INFINITY && typ(b) != t_INFINITY && gsigne(b))
    2459             :   {
    2460           7 :     fb = fujiwara_bound_real(P, -1);
    2461           7 :     if (fb <= -pariINFINITY) a = gen_0;
    2462           0 :     else if (fb < 0) a = gen_m1;
    2463           0 :     else a = negi(int2n((long)ceil(fb)));
    2464             :   }
    2465       63795 :   if (typ(b) == t_INFINITY && typ(a) != t_INFINITY && gsigne(a))
    2466             :   {
    2467          21 :     fb = fujiwara_bound_real(P, 1);
    2468          21 :     if (fb <= -pariINFINITY) b = gen_0;
    2469          21 :     else if (fb < 0) b = gen_1;
    2470           7 :     else b = int2n((long)ceil(fb));
    2471             :   }
    2472       63795 :   if (typ(a) != t_INFINITY && typ(b) != t_INFINITY)
    2473             :   {
    2474             :     GEN d, ad, bd, diff;
    2475             :     long i;
    2476             :     /* can occur if one of a,b was initially a t_INFINITY */
    2477        7210 :     if (gequal(a,b)) return gerepilecopy(av, sol);
    2478        7210 :     d = lcmii(Q_denom(a), Q_denom(b));
    2479        7210 :     if (is_pm1(d)) { d = NULL; ad = a; bd = b; }
    2480             :     else
    2481          21 :     { P = ZX_rescale(P, d); ad = gmul(a, d); bd = gmul(b, d); }
    2482        7210 :     diff = subii(bd, ad);
    2483        7210 :     P = ZX_unscale(ZX_translate(P, ad), diff);
    2484        7210 :     res = usp(P, flag, bitprec);
    2485        7210 :     if (flag <= 1)
    2486             :     {
    2487       19397 :       for (i = 1; i < lg(res); i++)
    2488             :       {
    2489       12243 :         GEN z = gmul(diff, gel(res, i));
    2490       12243 :         if (typ(z) == t_VEC)
    2491             :         {
    2492           0 :           gel(z, 1) = gadd(ad, gel(z, 1));
    2493           0 :           gel(z, 2) = gadd(ad, gel(z, 2));
    2494             :         }
    2495             :         else
    2496       12243 :           z = gadd(ad, z);
    2497       12243 :         if (d) z = gdiv(z, d);
    2498       12243 :         gel(res, i) = z;
    2499             :       }
    2500        7154 :       sol = shallowconcat(sol, res);
    2501             :     }
    2502             :     else
    2503          56 :       nbz += lg(res) - 1;
    2504             :   }
    2505       63795 :   if (typ(b) == t_INFINITY && (fb=fujiwara_bound_real(P, 1)) > -pariINFINITY)
    2506             :   {
    2507       54098 :     long bp = maxss((long)ceil(fb), 0);
    2508       54098 :     res = usp(ZX_unscale2n(P, bp), flag, bitprec);
    2509       54098 :     if (flag <= 1)
    2510       41063 :       sol = shallowconcat(sol, gmul2n(res, bp));
    2511             :     else
    2512       13035 :       nbz += lg(res)-1;
    2513             :   }
    2514       63795 :   if (typ(a) == t_INFINITY && (fb=fujiwara_bound_real(P,-1)) > -pariINFINITY)
    2515             :   {
    2516        9061 :     long i, bm = maxss((long)ceil(fb), 0);
    2517        9061 :     res = usp(ZX_unscale2n(ZX_z_unscale(P, -1), bm), flag, bitprec);
    2518        9061 :     if (flag <= 1)
    2519             :     {
    2520        9129 :       for (i = 1; i < lg(res); i++)
    2521             :       {
    2522        5824 :         GEN z = gneg(gmul2n(gel(res, i), bm));
    2523        5824 :         if (typ(z) == t_VEC) swap(gel(z, 1), gel(z, 2));
    2524        5824 :         gel(res, i) = z;
    2525             :       }
    2526        3305 :       sol = shallowconcat(res, sol);
    2527             :     }
    2528             :     else
    2529        5756 :       nbz += lg(res)-1;
    2530             :   }
    2531       63795 :   if (flag >= 2) return utoi(nbz);
    2532       48252 :   if (flag)
    2533       48252 :     sol = sort(sol);
    2534             :   else
    2535           0 :     sol = gen_sort(sol, (void *)_intervalcmp, cmp_nodata);
    2536       48252 :   return gerepileupto(av, sol);
    2537             : }
    2538             : 
    2539             : /* x a scalar */
    2540             : static GEN
    2541          35 : rootsdeg0(GEN x)
    2542             : {
    2543          35 :   if (!is_real_t(typ(x))) pari_err_TYPE("realroots",x);
    2544          28 :   if (gequal0(x)) pari_err_ROOTS0("realroots");
    2545          14 :   return cgetg(1,t_COL); /* constant polynomial */
    2546             : }
    2547             : static void
    2548       80002 : checkbound(GEN a)
    2549             : {
    2550       80002 :   switch(typ(a))
    2551             :   {
    2552       80002 :     case t_INT: case t_FRAC: case t_INFINITY: break;
    2553           0 :     default: pari_err_TYPE("polrealroots", a);
    2554             :   }
    2555       80002 : }
    2556             : static GEN
    2557       40519 : check_ab(GEN ab)
    2558             : {
    2559             :   GEN a, b;
    2560       40519 :   if (!ab) return NULL;
    2561       40001 :   if (typ(ab) != t_VEC || lg(ab) != 3) pari_err_TYPE("polrootsreal",ab);
    2562       40001 :   a = gel(ab,1); checkbound(a);
    2563       40001 :   b = gel(ab,2); checkbound(b);
    2564       40771 :   if (typ(a) == t_INFINITY && inf_get_sign(a) < 0 &&
    2565        1505 :       typ(b) == t_INFINITY && inf_get_sign(b) > 0) ab = NULL;
    2566       40001 :   return ab;
    2567             : }
    2568             : /* e^(1/h) assuming the h-th root is real, beware that sqrtnr assumes e >= 0 */
    2569             : static GEN
    2570        3599 : _sqrtnr(GEN e, long h)
    2571             : {
    2572             :   long s;
    2573             :   GEN r;
    2574        3599 :   if (h == 2) return sqrtr(e);
    2575          14 :   s = signe(e); setsigne(e, 1); /* e < 0 is possible, implies h is odd */
    2576          14 :   r = sqrtnr(e, h); if (s < 0) setsigne(r, -1);
    2577          14 :   return r;
    2578             : }
    2579             : GEN
    2580       37660 : realroots(GEN P, GEN ab, long prec)
    2581             : {
    2582       37660 :   pari_sp av = avma;
    2583       37660 :   GEN sol = NULL, fa, ex;
    2584             :   long i, j, v, l;
    2585             : 
    2586       37660 :   ab = check_ab(ab);
    2587       37660 :   if (typ(P) != t_POL) return rootsdeg0(P);
    2588       37639 :   switch(degpol(P))
    2589             :   {
    2590           7 :     case -1: return rootsdeg0(gen_0);
    2591           7 :     case 0: return rootsdeg0(gel(P,2));
    2592             :   }
    2593       37625 :   if (!RgX_is_ZX(P)) P = RgX_rescale_to_int(P);
    2594       37625 :   v = ZX_valrem(Q_primpart(P), &P);
    2595       37625 :   fa = ZX_squff(P, &ex); l = lg(fa); sol = cgetg(l + 1, t_VEC);
    2596       76369 :   for (i = 1; i < l; i++)
    2597             :   {
    2598       38744 :     GEN Pi = gel(fa, i), soli, soli2;
    2599             :     long n, h;
    2600       38744 :     if (ab) h = 1; else Pi = ZX_deflate_max(Pi, &h);
    2601       38744 :     soli = ZX_Uspensky(Pi, odd(h)? ab: gen_0, 1, prec2nbits(prec));
    2602       38744 :     n = lg(soli); soli2 = odd(h)? NULL: cgetg(n, t_COL);
    2603       82543 :     for (j = 1; j < n; j++)
    2604             :     {
    2605       43799 :       GEN r = gel(soli, j); /* != 0 */
    2606       43799 :       if (typ(r) != t_REAL) gel(soli, j) = r = gtofp(r, prec);
    2607       43799 :       if (h > 1)
    2608             :       {
    2609          77 :         gel(soli, j) = r = _sqrtnr(r, h);
    2610          77 :         if (soli2) gel(soli2, j) = negr(r);
    2611             :       }
    2612             :     }
    2613       38744 :     if (soli2) soli = shallowconcat(soli, soli2);
    2614       38744 :     if (ex[i] > 1) soli = shallowconcat1( const_vec(ex[i], soli) );
    2615       38744 :     gel(sol, i) = soli;
    2616             :   }
    2617       37625 :   if (v && (!ab || (gsigne(gel(ab,1)) <= 0 && gsigne(gel(ab,2)) >= 0)))
    2618          42 :     gel(sol, i++) = const_col(v, real_0(prec));
    2619       37625 :   setlg(sol, i); if (i == 1) { set_avma(av); return cgetg(1,t_COL); }
    2620       37611 :   return gerepileupto(av, sort(shallowconcat1(sol)));
    2621             : }
    2622             : GEN
    2623        7506 : ZX_realroots_irred(GEN P, long prec)
    2624             : {
    2625        7506 :   long dP = degpol(P), j, n, h;
    2626             :   GEN sol, sol2;
    2627             :   pari_sp av;
    2628        7506 :   if (dP == 1) retmkvec(ZX_deg1root(P, prec));
    2629        5952 :   av = avma; P = ZX_deflate_max(P, &h);
    2630        5952 :   if (h == dP)
    2631             :   {
    2632        2192 :     GEN r = _sqrtnr(ZX_deg1root(P, prec), h);
    2633        2192 :     return gerepilecopy(av, odd(h)? mkvec(r): mkvec2(negr(r), r));
    2634             :   }
    2635        3760 :   sol = ZX_Uspensky(P, odd(h)? NULL: gen_0, 1 | 4, prec2nbits(prec));
    2636        3760 :   n = lg(sol); sol2 = odd(h)? NULL: cgetg(n, t_COL);
    2637       17606 :   for (j = 1; j < n; j++)
    2638             :   {
    2639       13846 :     GEN r = gel(sol, j);
    2640       13846 :     if (typ(r) != t_REAL) gel(sol, j) = r = gtofp(r, prec);
    2641       13846 :     if (h > 1)
    2642             :     {
    2643        1330 :       gel(sol, j) = r = _sqrtnr(r, h);
    2644        1330 :       if (sol2) gel(sol2, j) = negr(r);
    2645             :     }
    2646             :   }
    2647        3760 :   if (sol2) sol = shallowconcat(sol, sol2);
    2648        3760 :   return gerepileupto(av, sort(sol));
    2649             : }
    2650             : 
    2651             : static long
    2652       21770 : ZX_sturm_i(GEN P, long flag)
    2653             : {
    2654             :   pari_sp av;
    2655       21770 :   long h, r, dP = degpol(P);
    2656       21770 :   if (dP == 1) return 1;
    2657       20027 :   av = avma; P = ZX_deflate_max(P, &h);
    2658       20027 :   if (h == dP)
    2659             :   { /* now deg P = 1 */
    2660        6384 :     if (odd(h))
    2661         385 :       r = 1;
    2662             :     else
    2663        5999 :       r = (signe(gel(P,2)) != signe(gel(P,3)))? 2: 0;
    2664        6384 :     return gc_long(av, r);
    2665             :   }
    2666       13643 :   if (odd(h))
    2667       11746 :     r = itou(ZX_Uspensky(P, NULL, flag, 0));
    2668             :   else
    2669        1897 :     r = 2*itou(ZX_Uspensky(P, gen_0, flag, 0));
    2670       13643 :   return gc_long(av,r);
    2671             : }
    2672             : /* P non-constant, squarefree ZX */
    2673             : long
    2674        2859 : ZX_sturmpart(GEN P, GEN ab)
    2675             : {
    2676        2859 :   pari_sp av = avma;
    2677        2859 :   if (!check_ab(ab)) return ZX_sturm(P);
    2678        2012 :   return gc_long(av, itou(ZX_Uspensky(P, ab, 2, 0)));
    2679             : }
    2680             : /* P non-constant, squarefree ZX */
    2681             : long
    2682         847 : ZX_sturm(GEN P) { return ZX_sturm_i(P, 2); }
    2683             : /* P irreducible ZX */
    2684             : long
    2685       20923 : ZX_sturm_irred(GEN P) { return ZX_sturm_i(P, 2 + 4); }

Generated by: LCOV version 1.13