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

Generated by: LCOV version 1.13