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.18.1 lcov report (development 30005-fc14bb602a) Lines: 1467 1530 95.9 %
Date: 2025-02-18 09:22:46 Functions: 116 119 97.5 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : /*******************************************************************/
      16             : /*                                                                 */
      17             : /*                ROOTS OF COMPLEX POLYNOMIALS                     */
      18             : /*  (original code contributed by Xavier Gourdon, INRIA RR 1852)   */
      19             : /*                                                                 */
      20             : /*******************************************************************/
      21             : #include "pari.h"
      22             : #include "paripriv.h"
      23             : 
      24             : #define DEBUGLEVEL DEBUGLEVEL_polroots
      25             : 
      26             : static const double pariINFINITY = 1./0.;
      27             : 
      28             : static long
      29     1227588 : isvalidcoeff(GEN x)
      30             : {
      31     1227588 :   switch (typ(x))
      32             :   {
      33     1203615 :     case t_INT: case t_REAL: case t_FRAC: return 1;
      34       23959 :     case t_COMPLEX: return isvalidcoeff(gel(x,1)) && isvalidcoeff(gel(x,2));
      35             :   }
      36          14 :   return 0;
      37             : }
      38             : 
      39             : static void
      40      266906 : checkvalidpol(GEN p, const char *f)
      41             : {
      42      266906 :   long i,n = lg(p);
      43     1446555 :   for (i=2; i<n; i++)
      44     1179656 :     if (!isvalidcoeff(gel(p,i))) pari_err_TYPE(f, gel(p,i));
      45      266899 : }
      46             : 
      47             : /********************************************************************/
      48             : /**                                                                **/
      49             : /**                   FAST ARITHMETIC over Z[i]                    **/
      50             : /**                                                                **/
      51             : /********************************************************************/
      52             : 
      53             : static GEN
      54    16457252 : ZX_to_ZiX(GEN Pr, GEN Pi)
      55             : {
      56    16457252 :   long i, lr = lg(Pr), li = lg(Pi), l = maxss(lr, li), m = minss(lr, li);
      57    16459701 :   GEN P = cgetg(l, t_POL);
      58    16464637 :   P[1] = Pr[1];
      59    66781312 :   for(i = 2; i < m; i++)
      60    50319526 :     gel(P,i) = signe(gel(Pi,i)) ? mkcomplex(gel(Pr,i), gel(Pi,i))
      61    50319526 :                                 : gel(Pr,i);
      62    22363985 :   for(     ; i < lr; i++)
      63     5902199 :     gel(P,i) = gel(Pr, i);
      64    16496741 :   for(     ; i < li; i++)
      65       34955 :     gel(P,i) = mkcomplex(gen_0, gel(Pi, i));
      66    16461786 :   return normalizepol_lg(P, l);
      67             : }
      68             : 
      69             : static GEN
      70    99019268 : ZiX_sqr(GEN P)
      71             : {
      72    99019268 :   pari_sp av = avma;
      73             :   GEN Pr2, Pi2, Qr, Qi;
      74    99019268 :   GEN Pr = real_i(P), Pi = imag_i(P);
      75    99012422 :   if (signe(Pi)==0) return gerepileupto(av, ZX_sqr(Pr));
      76    16521683 :   if (signe(Pr)==0) return gerepileupto(av, ZX_neg(ZX_sqr(Pi)));
      77    16464310 :   Pr2 = ZX_sqr(Pr); Pi2 = ZX_sqr(Pi);
      78    16456326 :   Qr = ZX_sub(Pr2, Pi2);
      79    16457373 :   if (degpol(Pr)==degpol(Pi))
      80    10703015 :     Qi = ZX_sub(ZX_sqr(ZX_add(Pr, Pi)), ZX_add(Pr2, Pi2));
      81             :   else
      82     5758598 :     Qi = ZX_shifti(ZX_mul(Pr, Pi), 1);
      83    16459397 :   return gerepilecopy(av, ZX_to_ZiX(Qr, Qi));
      84             : }
      85             : 
      86             : static GEN
      87    49497945 : graeffe(GEN p)
      88             : {
      89    49497945 :   pari_sp av = avma;
      90             :   GEN p0, p1, s0, s1;
      91    49497945 :   long n = degpol(p);
      92             : 
      93    49505120 :   if (!n) return RgX_copy(p);
      94    49505120 :   RgX_even_odd(p, &p0, &p1);
      95             :   /* p = p0(x^2) + x p1(x^2) */
      96    49512559 :   s0 = ZiX_sqr(p0);
      97    49518469 :   s1 = ZiX_sqr(p1);
      98    49517929 :   return gerepileupto(av, RgX_sub(s0, RgX_shift_shallow(s1,1)));
      99             : }
     100             : 
     101             : GEN
     102        5383 : ZX_graeffe(GEN p)
     103             : {
     104        5383 :   pari_sp av = avma;
     105             :   GEN p0, p1, s0, s1;
     106        5383 :   long n = degpol(p);
     107             : 
     108        5383 :   if (!n) return ZX_copy(p);
     109        5383 :   RgX_even_odd(p, &p0, &p1);
     110             :   /* p = p0(x^2) + x p1(x^2) */
     111        5383 :   s0 = ZX_sqr(p0);
     112        5383 :   s1 = ZX_sqr(p1);
     113        5383 :   return gerepileupto(av, ZX_sub(s0, RgX_shift_shallow(s1,1)));
     114             : }
     115             : GEN
     116          14 : polgraeffe(GEN p)
     117             : {
     118          14 :   pari_sp av = avma;
     119             :   GEN p0, p1, s0, s1;
     120          14 :   long n = degpol(p);
     121             : 
     122          14 :   if (typ(p) != t_POL) pari_err_TYPE("polgraeffe",p);
     123          14 :   n = degpol(p);
     124          14 :   if (!n) return gcopy(p);
     125          14 :   RgX_even_odd(p, &p0, &p1);
     126             :   /* p = p0(x^2) + x p1(x^2) */
     127          14 :   s0 = RgX_sqr(p0);
     128          14 :   s1 = RgX_sqr(p1);
     129          14 :   return gerepileupto(av, RgX_sub(s0, RgX_shift_shallow(s1,1)));
     130             : }
     131             : 
     132             : /********************************************************************/
     133             : /**                                                                **/
     134             : /**                       MODULUS OF ROOTS                         **/
     135             : /**                                                                **/
     136             : /********************************************************************/
     137             : 
     138             : /* Quick approximation to log2(|x|); first define y s.t. |y-x| < 2^-32 then
     139             :  * return y rounded to 2 ulp. In particular, if result < 2^21, absolute error
     140             :  * is bounded by 2^-31. If result > 2^21, it is correct to 2 ulp */
     141             : static double
     142   224334778 : mydbllog2i(GEN x)
     143             : {
     144             : #ifdef LONG_IS_64BIT
     145   192977732 :   const double W = 1/(4294967296. * 4294967296.); /* 2^-64 */
     146             : #else
     147    31357046 :   const double W = 1/4294967296.; /*2^-32*/
     148             : #endif
     149             :   GEN m;
     150   224334778 :   long lx = lgefint(x);
     151             :   double l;
     152   224334778 :   if (lx == 2) return -pariINFINITY;
     153   223517813 :   m = int_MSW(x);
     154   223517813 :   l = (double)(ulong)*m;
     155   223517813 :   if (lx == 3) return log2(l);
     156    70086218 :   l += ((double)(ulong)*int_precW(m)) * W;
     157             :   /* at least m = min(53,BIL) bits are correct in the mantissa, thus log2
     158             :    * is correct with error < log(1 + 2^-m) ~ 2^-m. Adding the correct
     159             :    * exponent BIL(lx-3) causes 1ulp further round-off error */
     160    70086218 :   return log2(l) + (double)(BITS_IN_LONG*(lx-3));
     161             : }
     162             : 
     163             : /* return log(|x|) or -pariINFINITY */
     164             : static double
     165     9540649 : mydbllogr(GEN x) {
     166     9540649 :   if (!signe(x)) return -pariINFINITY;
     167     9540649 :   return M_LN2*dbllog2r(x);
     168             : }
     169             : 
     170             : /* return log2(|x|) or -pariINFINITY */
     171             : static double
     172    55902260 : mydbllog2r(GEN x) {
     173    55902260 :   if (!signe(x)) return -pariINFINITY;
     174    55462310 :   return dbllog2r(x);
     175             : }
     176             : double
     177   300262029 : dbllog2(GEN z)
     178             : {
     179             :   double x, y;
     180   300262029 :   switch(typ(z))
     181             :   {
     182   224227300 :     case t_INT: return mydbllog2i(z);
     183       22363 :     case t_FRAC: return mydbllog2i(gel(z,1))-mydbllog2i(gel(z,2));
     184    48956683 :     case t_REAL: return mydbllog2r(z);
     185    27055683 :     default: /*t_COMPLEX*/
     186    27055683 :       x = dbllog2(gel(z,1));
     187    27136822 :       y = dbllog2(gel(z,2));
     188    27136681 :       if (x == -pariINFINITY) return y;
     189    26893344 :       if (y == -pariINFINITY) return x;
     190    26688170 :       if (fabs(x-y) > 10) return maxdd(x,y);
     191    26059812 :       return x + 0.5*log2(1 + exp2(2*(y-x)));
     192             :   }
     193             : }
     194             : static GEN /* beware overflow */
     195     6546655 : dblexp(double x) { return fabs(x) < 100.? dbltor(exp(x)): mpexp(dbltor(x)); }
     196             : 
     197             : /* find s such that  A_h <= 2^s <= 2 A_i  for one h and all i < n = deg(p),
     198             :  * with  A_i := (binom(n,i) lc(p) / p_i) ^ 1/(n-i), and  p = sum p_i X^i */
     199             : static long
     200    40955251 : findpower(GEN p)
     201             : {
     202    40955251 :   double x, L, mins = pariINFINITY;
     203    40955251 :   long n = degpol(p),i;
     204             : 
     205    40954708 :   L = dbllog2(gel(p,n+2)); /* log2(lc * binom(n,i)) */
     206   164659331 :   for (i=n-1; i>=0; i--)
     207             :   {
     208   123703226 :     L += log2((double)(i+1) / (double)(n-i));
     209   123703226 :     x = dbllog2(gel(p,i+2));
     210   123700540 :     if (x != -pariINFINITY)
     211             :     {
     212   122910018 :       double s = (L - x) / (double)(n-i);
     213   122910018 :       if (s < mins) mins = s;
     214             :     }
     215             :   }
     216    40956105 :   i = (long)ceil(mins);
     217    40956105 :   if (i - mins > 1 - 1e-12) i--;
     218    40956105 :   return i;
     219             : }
     220             : 
     221             : /* returns the exponent for logmodulus(), from the Newton diagram */
     222             : static long
     223     5500920 : newton_polygon(GEN p, long k)
     224             : {
     225     5500920 :   pari_sp av = avma;
     226     5500920 :   long n = degpol(p), i, j, h, l, *vertex = (long*)new_chunk(n+1);
     227     5500893 :   double *L = (double*)stack_malloc_align((n+1)*sizeof(double), sizeof(double));
     228             : 
     229             :   /* vertex[i] = 1 if i a vertex of convex hull, 0 otherwise */
     230    26426673 :   for (i=0; i<=n; i++) { L[i] = dbllog2(gel(p,2+i)); vertex[i] = 0; }
     231     5500895 :   vertex[0] = 1; /* sentinel */
     232    19584170 :   for (i=0; i < n; i=h)
     233             :   {
     234             :     double slope;
     235    14083275 :     h = i+1;
     236    14088149 :     while (L[i] == -pariINFINITY) { vertex[h] = 1; i = h; h = i+1; }
     237    14083275 :     slope = L[h] - L[i];
     238    37929551 :     for (j = i+2; j<=n; j++) if (L[j] != -pariINFINITY)
     239             :     {
     240    23840061 :       double pij = (L[j] - L[i])/(double)(j - i);
     241    23840061 :       if (slope < pij) { slope = pij; h = j; }
     242             :     }
     243    14083275 :     vertex[h] = 1;
     244             :   }
     245     6216450 :   h = k;   while (!vertex[h]) h++;
     246     5689371 :   l = k-1; while (!vertex[l]) l--;
     247     5500895 :   set_avma(av);
     248     5500922 :   return (long)floor((L[h]-L[l])/(double)(h-l) + 0.5);
     249             : }
     250             : 
     251             : /* change z into z*2^e, where z is real or complex of real */
     252             : static void
     253    35494757 : myshiftrc(GEN z, long e)
     254             : {
     255    35494757 :   if (typ(z)==t_COMPLEX)
     256             :   {
     257     6417353 :     if (signe(gel(z,1))) shiftr_inplace(gel(z,1), e);
     258     6417362 :     if (signe(gel(z,2))) shiftr_inplace(gel(z,2), e);
     259             :   }
     260             :   else
     261    29077404 :     if (signe(z)) shiftr_inplace(z, e);
     262    35494862 : }
     263             : 
     264             : /* return z*2^e, where z is integer or complex of integer (destroy z) */
     265             : static GEN
     266   135808929 : myshiftic(GEN z, long e)
     267             : {
     268   135808929 :   if (typ(z)==t_COMPLEX)
     269             :   {
     270    17901219 :     gel(z,1) = signe(gel(z,1))? mpshift(gel(z,1),e): gen_0;
     271    17900326 :     gel(z,2) = mpshift(gel(z,2),e);
     272    17897838 :     return z;
     273             :   }
     274   117907710 :   return signe(z)? mpshift(z,e): gen_0;
     275             : }
     276             : 
     277             : static GEN
     278     7017941 : RgX_gtofp_bit(GEN q, long bit) { return RgX_gtofp(q, nbits2prec(bit)); }
     279             : 
     280             : static GEN
     281   214822884 : mygprecrc(GEN x, long prec, long e)
     282             : {
     283             :   GEN y;
     284   214822884 :   switch(typ(x))
     285             :   {
     286   160319053 :     case t_REAL:
     287   160319053 :       if (!signe(x)) return real_0_bit(e);
     288   157033374 :       return realprec(x) == prec? x: rtor(x, prec);
     289    37145750 :     case t_COMPLEX:
     290    37145750 :       y = cgetg(3,t_COMPLEX);
     291    37145153 :       gel(y,1) = mygprecrc(gel(x,1),prec,e);
     292    37145183 :       gel(y,2) = mygprecrc(gel(x,2),prec,e);
     293    37145149 :       return y;
     294    17358081 :     default: return x;
     295             :   }
     296             : }
     297             : 
     298             : /* gprec behaves badly with the zero for polynomials.
     299             : The second parameter in mygprec is the precision in base 2 */
     300             : static GEN
     301    63655072 : mygprec(GEN x, long bit)
     302             : {
     303             :   long e, prec;
     304    63655072 :   if (bit < 0) bit = 0; /* should rarely happen */
     305    63655072 :   e = gexpo(x) - bit;
     306    63657010 :   prec = nbits2prec(bit);
     307    63662533 :   switch(typ(x))
     308             :   {
     309   167915051 :     case t_POL: pari_APPLY_pol_normalized(mygprecrc(gel(x,i),prec,e));
     310    17207244 :     default: return mygprecrc(x,prec,e);
     311             :   }
     312             : }
     313             : 
     314             : /* normalize a polynomial x, that is change it with coefficients in Z[i],
     315             : after making product by 2^shift */
     316             : static GEN
     317    17579131 : pol_to_gaussint(GEN x, long shift)
     318   100281296 : { pari_APPLY_pol_normalized(gtrunc2n(gel(x,i), shift)); }
     319             : 
     320             : /* returns a polynomial q in Z[i][x] keeping bit bits of p */
     321             : static GEN
     322    13128726 : eval_rel_pol(GEN p, long bit)
     323             : {
     324             :   long i;
     325    73951829 :   for (i = 2; i < lg(p); i++)
     326    60823106 :     if (gequal0(gel(p,i))) gel(p,i) = gen_0; /* bad behavior of gexpo */
     327    13128723 :   return pol_to_gaussint(p, bit-gexpo(p)+1);
     328             : }
     329             : 
     330             : /* returns p(R*x)/R^n (in R or R[i]), R = exp(lrho), bit bits of precision */
     331             : static GEN
     332     1608940 : homothetie(GEN p, double lrho, long bit)
     333             : {
     334             :   GEN q, r, t, iR;
     335     1608940 :   long n = degpol(p), i;
     336             : 
     337     1608939 :   iR = mygprec(dblexp(-lrho),bit);
     338     1608906 :   q = mygprec(p, bit);
     339     1608938 :   r = cgetg(n+3,t_POL); r[1] = p[1];
     340     1608933 :   t = iR; r[n+2] = q[n+2];
     341     6788233 :   for (i=n-1; i>0; i--)
     342             :   {
     343     5179326 :     gel(r,i+2) = gmul(t, gel(q,i+2));
     344     5179267 :     t = mulrr(t, iR);
     345             :   }
     346     1608907 :   gel(r,2) = gmul(t, gel(q,2)); return r;
     347             : }
     348             : 
     349             : /* change q in 2^(n*e) p(x*2^(-e)), n=deg(q)  [ ~as above with R = 2^-e ]*/
     350             : static void
     351     9950971 : homothetie2n(GEN p, long e)
     352             : {
     353     9950971 :   if (e)
     354             :   {
     355     7972007 :     long i,n = lg(p)-1;
     356    43466782 :     for (i=2; i<=n; i++) myshiftrc(gel(p,i), (n-i)*e);
     357             :   }
     358     9951111 : }
     359             : 
     360             : /* return 2^f * 2^(n*e) p(x*2^(-e)), n=deg(q) */
     361             : static void
     362    36496795 : homothetie_gauss(GEN p, long e, long f)
     363             : {
     364    36496795 :   if (e || f)
     365             :   {
     366    32603282 :     long i, n = lg(p)-1;
     367   168364171 :     for (i=2; i<=n; i++) gel(p,i) = myshiftic(gel(p,i), f+(n-i)*e);
     368             :   }
     369    36448702 : }
     370             : 
     371             : /* Lower bound on the modulus of the largest root z_0
     372             :  * k is set to an upper bound for #{z roots, |z-z_0| < eps} */
     373             : static double
     374    40954519 : lower_bound(GEN p, long *k, double eps)
     375             : {
     376    40954519 :   long n = degpol(p), i, j;
     377    40954776 :   pari_sp ltop = avma;
     378             :   GEN a, s, S, ilc;
     379             :   double r, R, rho;
     380             : 
     381    40954776 :   if (n < 4) { *k = n; return 0.; }
     382     8232223 :   S = cgetg(5,t_VEC);
     383     8232936 :   a = cgetg(5,t_VEC); ilc = gdiv(real_1(DEFAULTPREC), gel(p,n+2));
     384    41141894 :   for (i=1; i<=4; i++) gel(a,i) = gmul(ilc,gel(p,n+2-i));
     385             :   /* i = 1 split out from next loop for efficiency and initialization */
     386     8231739 :   s = gel(a,1);
     387     8231739 :   gel(S,1) = gneg(s); /* Newton sum S_i */
     388     8232447 :   rho = r = gtodouble(gabs(s,3));
     389     8232900 :   R = r / n;
     390    32925532 :   for (i=2; i<=4; i++)
     391             :   {
     392    24693011 :     s = gmulsg(i,gel(a,i));
     393    73992466 :     for (j=1; j<i; j++) s = gadd(s, gmul(gel(S,j),gel(a,i-j)));
     394    24676926 :     gel(S,i) = gneg(s); /* Newton sum S_i */
     395    24688386 :     r = gtodouble(gabs(s,3));
     396    24692632 :     if (r > 0.)
     397             :     {
     398    24646671 :       r = exp(log(r/n) / (double)i);
     399    24646671 :       if (r > R) R = r;
     400             :     }
     401             :   }
     402     8232521 :   if (R > 0. && eps < 1.2)
     403     8228794 :     *k = (long)floor((rho/R + n) / (1 + exp(-eps)*cos(eps)));
     404             :   else
     405        3727 :     *k = n;
     406     8232521 :   return gc_double(ltop, R);
     407             : }
     408             : 
     409             : /* return R such that exp(R - tau) <= rho_n(P) <= exp(R + tau)
     410             :  * P(0) != 0 and P non constant */
     411             : static double
     412     4450456 : logmax_modulus(GEN p, double tau)
     413             : {
     414             :   GEN r, q, aux, gunr;
     415     4450456 :   pari_sp av, ltop = avma;
     416     4450456 :   long i,k,n=degpol(p),nn,bit,M,e;
     417     4450446 :   double rho,eps, tau2 = (tau > 3.0)? 0.5: tau/6.;
     418             : 
     419     4450446 :   r = cgeti(BIGDEFAULTPREC);
     420     4450427 :   av = avma;
     421             : 
     422     4450427 :   eps = - 1/log(1.5*tau2); /* > 0 */
     423     4450427 :   bit = (long) ((double) n*log2(1./tau2)+3*log2((double) n))+1;
     424     4450427 :   gunr = real_1_bit(bit+2*n);
     425     4450393 :   aux = gdiv(gunr, gel(p,2+n));
     426     4450341 :   q = RgX_Rg_mul(p, aux); gel(q,2+n) = gunr;
     427     4450225 :   e = findpower(q);
     428     4450352 :   homothetie2n(q,e);
     429     4450402 :   affsi(e, r);
     430     4450400 :   q = pol_to_gaussint(q, bit);
     431     4450032 :   M = (long) (log2( log(4.*n) / (2*tau2) )) + 2;
     432     4450032 :   nn = n;
     433     4450032 :   for (i=0,e=0;;)
     434             :   { /* nn = deg(q) */
     435    40957035 :     rho = lower_bound(q, &k, eps);
     436    40954516 :     if (rho > exp2(-(double)e)) e = (long)-floor(log2(rho));
     437    40954516 :     affii(shifti(addis(r,e), 1), r);
     438    40913111 :     if (++i == M) break;
     439             : 
     440    36463412 :     bit = (long) ((double)k * log2(1./tau2) +
     441    36463412 :                      (double)(nn-k)*log2(1./eps) + 3*log2((double)nn)) + 1;
     442    36463412 :     homothetie_gauss(q, e, bit-(long)floor(dbllog2(gel(q,2+nn))+0.5));
     443    36479902 :     nn -= RgX_valrem(q, &q);
     444    36484036 :     q = gerepileupto(av, graeffe(q));
     445    36506453 :     tau2 *= 1.5; if (tau2 > 0.9) tau2 = 0.5;
     446    36506453 :     eps = -1/log(tau2); /* > 0 */
     447    36506453 :     e = findpower(q);
     448             :   }
     449     4449699 :   if (!signe(r)) return gc_double(ltop,0.);
     450     4027362 :   r = itor(r, DEFAULTPREC); shiftr_inplace(r, -M);
     451     4027812 :   return gc_double(ltop, -rtodbl(r) * M_LN2); /* -log(2) sum e_i 2^-i */
     452             : }
     453             : 
     454             : static GEN
     455       35853 : RgX_normalize1(GEN x)
     456             : {
     457       35853 :   long i, n = lg(x)-1;
     458             :   GEN y;
     459       35867 :   for (i = n; i > 1; i--)
     460       35860 :     if (!gequal0( gel(x,i) )) break;
     461       35853 :   if (i == n) return x;
     462          14 :   pari_warn(warner,"normalizing a polynomial with 0 leading term");
     463          14 :   if (i == 1) pari_err_ROOTS0("roots");
     464          14 :   y = cgetg(i+1, t_POL); y[1] = x[1];
     465          42 :   for (; i > 1; i--) gel(y,i) = gel(x,i);
     466          14 :   return y;
     467             : }
     468             : 
     469             : static GEN
     470       27259 : polrootsbound_i(GEN P, double TAU)
     471             : {
     472       27259 :   pari_sp av = avma;
     473             :   double d;
     474       27259 :   (void)RgX_valrem_inexact(P,&P);
     475       27259 :   P = RgX_normalize1(P);
     476       27259 :   switch(degpol(P))
     477             :   {
     478           7 :     case -1: pari_err_ROOTS0("roots");
     479         126 :     case 0:  set_avma(av); return gen_0;
     480             :   }
     481       27126 :   d = logmax_modulus(P, TAU) + TAU;
     482             :   /* not dblexp: result differs on ARM emulator */
     483       27128 :   return gerepileuptoleaf(av, mpexp(dbltor(d)));
     484             : }
     485             : GEN
     486       27266 : polrootsbound(GEN P, GEN tau)
     487             : {
     488       27266 :   if (typ(P) != t_POL) pari_err_TYPE("polrootsbound",P);
     489       27259 :   checkvalidpol(P, "polrootsbound");
     490       27258 :   return polrootsbound_i(P, tau? gtodouble(tau): 0.01);
     491             : }
     492             : 
     493             : /* log of modulus of the smallest root of p, with relative error tau */
     494             : static double
     495     1614748 : logmin_modulus(GEN p, double tau)
     496             : {
     497     1614748 :   pari_sp av = avma;
     498     1614748 :   if (gequal0(gel(p,2))) return -pariINFINITY;
     499     1614746 :   return gc_double(av, - logmax_modulus(RgX_recip_i(p),tau));
     500             : }
     501             : 
     502             : /* return the log of the k-th modulus (ascending order) of p, rel. error tau*/
     503             : static double
     504      606123 : logmodulus(GEN p, long k, double tau)
     505             : {
     506             :   GEN q;
     507      606123 :   long i, kk = k, imax, n = degpol(p), nn, bit, e;
     508      606123 :   pari_sp av, ltop=avma;
     509      606123 :   double r, tau2 = tau/6;
     510             : 
     511      606123 :   bit = (long)(n * (2. + log2(3.*n/tau2)));
     512      606123 :   av = avma;
     513      606123 :   q = gprec_w(p, nbits2prec(bit));
     514      606123 :   q = RgX_gtofp_bit(q, bit);
     515      606123 :   e = newton_polygon(q,k);
     516      606122 :   r = (double)e;
     517      606122 :   homothetie2n(q,e);
     518      606140 :   imax = (long)(log2(3./tau) + log2(log(4.*n)))+1;
     519     5500886 :   for (i=1; i<imax; i++)
     520             :   {
     521     4894767 :     q = eval_rel_pol(q,bit);
     522     4894453 :     kk -= RgX_valrem(q, &q);
     523     4894605 :     nn = degpol(q);
     524             : 
     525     4894607 :     q = gerepileupto(av, graeffe(q));
     526     4894796 :     e = newton_polygon(q,kk);
     527     4894808 :     r += e / exp2((double)i);
     528     4894808 :     q = RgX_gtofp_bit(q, bit);
     529     4894685 :     homothetie2n(q,e);
     530             : 
     531     4894746 :     tau2 *= 1.5; if (tau2 > 1.) tau2 = 1.;
     532     4894746 :     bit = 1 + (long)(nn*(2. + log2(3.*nn/tau2)));
     533             :   }
     534      606119 :   return gc_double(ltop, -r * M_LN2);
     535             : }
     536             : 
     537             : /* return the log of the k-th modulus r_k of p, rel. error tau, knowing that
     538             :  * rmin < r_k < rmax. This information helps because we may reduce precision
     539             :  * quicker */
     540             : static double
     541      606117 : logpre_modulus(GEN p, long k, double tau, double lrmin, double lrmax)
     542             : {
     543             :   GEN q;
     544      606117 :   long n = degpol(p), i, imax, imax2, bit;
     545      606116 :   pari_sp ltop = avma, av;
     546      606116 :   double lrho, aux, tau2 = tau/6.;
     547             : 
     548      606116 :   aux = (lrmax - lrmin) / 2. + 4*tau2;
     549      606116 :   imax = (long) log2(log((double)n)/ aux);
     550      606116 :   if (imax <= 0) return logmodulus(p,k,tau);
     551             : 
     552      596992 :   lrho  = (lrmin + lrmax) / 2;
     553      596992 :   av = avma;
     554      596992 :   bit = (long)(n*(2. + aux / M_LN2 - log2(tau2)));
     555      596992 :   q = homothetie(p, lrho, bit);
     556      596990 :   imax2 = (long)(log2(3./tau * log(4.*n))) + 1;
     557      596990 :   if (imax > imax2) imax = imax2;
     558             : 
     559     1583509 :   for (i=0; i<imax; i++)
     560             :   {
     561      986510 :     q = eval_rel_pol(q,bit);
     562      986503 :     q = gerepileupto(av, graeffe(q));
     563      986521 :     aux = 2*aux + 2*tau2;
     564      986521 :     tau2 *= 1.5;
     565      986521 :     bit = (long)(n*(2. + aux / M_LN2 - log2(1-exp(-tau2))));
     566      986521 :     q = RgX_gtofp_bit(q, bit);
     567             :   }
     568      596999 :   aux = exp2((double)imax);
     569      596999 :   return gc_double(ltop, lrho + logmodulus(q,k, aux*tau/3.) / aux);
     570             : }
     571             : 
     572             : static double
     573      902540 : ind_maxlog2(GEN q)
     574             : {
     575      902540 :   long i, k = -1;
     576      902540 :   double L = - pariINFINITY;
     577     2221885 :   for (i=0; i<=degpol(q); i++)
     578             :   {
     579     1319343 :     double d = dbllog2(gel(q,2+i));
     580     1319345 :     if (d > L) { L = d; k = i; }
     581             :   }
     582      902540 :   return k;
     583             : }
     584             : 
     585             : /* Returns k such that r_k e^(-tau) < R < r_{k+1} e^tau.
     586             :  * Assume that l <= k <= n-l */
     587             : static long
     588     1011951 : dual_modulus(GEN p, double lrho, double tau, long l)
     589             : {
     590     1011951 :   long i, imax, delta_k = 0, n = degpol(p), nn, v2, v, bit, ll = l;
     591     1011951 :   double tau2 = tau * 7./8.;
     592     1011951 :   pari_sp av = avma;
     593             :   GEN q;
     594             : 
     595     1011951 :   bit = 6*n - 5*l + (long)(n*(-log2(tau2) + tau2 * 8./7.));
     596     1011951 :   q = homothetie(p, lrho, bit);
     597     1011942 :   imax = (long)(log(log(2.*n)/tau2)/log(7./4.)+1);
     598             : 
     599     8150250 :   for (i=0; i<imax; i++)
     600             :   {
     601     7247710 :     q = eval_rel_pol(q,bit); v2 = n - degpol(q);
     602     7246572 :     v = RgX_valrem(q, &q);
     603     7247117 :     ll -= maxss(v, v2); if (ll < 0) ll = 0;
     604             : 
     605     7247309 :     nn = degpol(q); delta_k += v;
     606     7247343 :     if (!nn) return delta_k;
     607             : 
     608     7137931 :     q = gerepileupto(av, graeffe(q));
     609     7138308 :     tau2 *= 7./4.;
     610     7138308 :     bit = 6*nn - 5*ll + (long)(nn*(-log2(tau2) + tau2 * 8./7.));
     611             :   }
     612      902540 :   return gc_long(av, delta_k + (long)ind_maxlog2(q));
     613             : }
     614             : 
     615             : /********************************************************************/
     616             : /**                                                                **/
     617             : /**              FACTORS THROUGH CIRCLE INTEGRATION                **/
     618             : /**                                                                **/
     619             : /********************************************************************/
     620             : /* l power of 2, W[step*j] = w_j; set f[j] = p(w_j)
     621             :  * if inv, w_j = exp(2IPi*j/l), else exp(-2IPi*j/l) */
     622             : 
     623             : static void
     624        7462 : fft2(GEN W, GEN p, GEN f, long step, long l)
     625             : {
     626             :   pari_sp av;
     627             :   long i, s1, l1, step2;
     628             : 
     629        7462 :   if (l == 2)
     630             :   {
     631        3766 :     gel(f,0) = gadd(gel(p,0), gel(p,step));
     632        3766 :     gel(f,1) = gsub(gel(p,0), gel(p,step)); return;
     633             :   }
     634        3696 :   av = avma;
     635        3696 :   l1 = l>>1; step2 = step<<1;
     636        3696 :   fft2(W,p,          f,   step2,l1);
     637        3696 :   fft2(W,p+step,     f+l1,step2,l1);
     638       32760 :   for (i = s1 = 0; i < l1; i++, s1 += step)
     639             :   {
     640       29064 :     GEN f0 = gel(f,i);
     641       29064 :     GEN f1 = gmul(gel(W,s1), gel(f,i+l1));
     642       29064 :     gel(f,i)    = gadd(f0, f1);
     643       29064 :     gel(f,i+l1) = gsub(f0, f1);
     644             :   }
     645        3696 :   gerepilecoeffs(av, f, l);
     646             : }
     647             : 
     648             : static void
     649    14138300 : fft(GEN W, GEN p, GEN f, long step, long l, long inv)
     650             : {
     651             :   pari_sp av;
     652             :   long i, s1, l1, l2, l3, step4;
     653             :   GEN f1, f2, f3, f02;
     654             : 
     655    14138300 :   if (l == 2)
     656             :   {
     657     6622294 :     gel(f,0) = gadd(gel(p,0), gel(p,step));
     658     6622024 :     gel(f,1) = gsub(gel(p,0), gel(p,step)); return;
     659             :   }
     660     7516006 :   av = avma;
     661     7516006 :   if (l == 4)
     662             :   {
     663             :     pari_sp av2;
     664     5319744 :     f1 = gadd(gel(p,0), gel(p,step<<1));
     665     5319187 :     f2 = gsub(gel(p,0), gel(p,step<<1));
     666     5319179 :     f3 = gadd(gel(p,step), gel(p,3*step));
     667     5319176 :     f02 = gsub(gel(p,step), gel(p,3*step));
     668     5319126 :     f02 = inv? mulcxI(f02): mulcxmI(f02);
     669     5319663 :     av2 = avma;
     670     5319663 :     gel(f,0) = gadd(f1, f3);
     671     5318985 :     gel(f,1) = gadd(f2, f02);
     672     5319176 :     gel(f,2) = gsub(f1, f3);
     673     5319046 :     gel(f,3) = gsub(f2, f02);
     674     5319295 :     gerepileallsp(av,av2,4,&gel(f,0),&gel(f,1),&gel(f,2),&gel(f,3));
     675     5319889 :     return;
     676             :   }
     677     2196262 :   l1 = l>>2; l2 = 2*l1; l3 = l1+l2; step4 = step<<2;
     678     2196262 :   fft(W,p,          f,   step4,l1,inv);
     679     2196671 :   fft(W,p+step,     f+l1,step4,l1,inv);
     680     2196659 :   fft(W,p+(step<<1),f+l2,step4,l1,inv);
     681     2196678 :   fft(W,p+3*step,   f+l3,step4,l1,inv);
     682     8180390 :   for (i = s1 = 0; i < l1; i++, s1 += step)
     683             :   {
     684     5983786 :     long s2 = s1 << 1, s3 = s1 + s2;
     685             :     GEN g02, g13, f13;
     686     5983786 :     f1 = gmul(gel(W,s1), gel(f,i+l1));
     687     5984072 :     f2 = gmul(gel(W,s2), gel(f,i+l2));
     688     5983984 :     f3 = gmul(gel(W,s3), gel(f,i+l3));
     689             : 
     690     5984084 :     f02 = gadd(gel(f,i),f2);
     691     5983477 :     g02 = gsub(gel(f,i),f2);
     692     5983521 :     f13 = gadd(f1,f3);
     693     5983399 :     g13 = gsub(f1,f3); g13 = inv? mulcxI(g13): mulcxmI(g13);
     694             : 
     695     5983920 :     gel(f,i)    = gadd(f02, f13);
     696     5983476 :     gel(f,i+l1) = gadd(g02, g13);
     697     5983477 :     gel(f,i+l2) = gsub(f02, f13);
     698     5983477 :     gel(f,i+l3) = gsub(g02, g13);
     699             :   }
     700     2196604 :   gerepilecoeffs(av, f, l);
     701             : }
     702             : 
     703             : static GEN
     704          98 : FFT_i(GEN W, GEN x)
     705             : {
     706          98 :   long i, l = lg(W), n = lg(x), tx = typ(x), tw, pa;
     707             :   GEN y, z, p, pol;
     708          98 :   if (l==1 || ((l-1) & (l-2))) pari_err_DIM("fft");
     709          84 :   tw = RgV_type(W, &p, &pol, &pa);
     710          84 :   if (tx == t_POL) { x++; n--; }
     711          49 :   else if (!is_vec_t(tx)) pari_err_TYPE("fft",x);
     712          84 :   if (n > l) pari_err_DIM("fft");
     713          84 :   if (n < l) {
     714           0 :     z = cgetg(l, t_VECSMALL); /* cf stackdummy */
     715           0 :     for (i = 1; i < n; i++) gel(z,i) = gel(x,i);
     716           0 :     for (     ; i < l; i++) gel(z,i) = gen_0;
     717             :   }
     718          84 :   else z = x;
     719          84 :   if (l == 2) return mkveccopy(gel(z,1));
     720          70 :   y = cgetg(l, t_VEC);
     721          70 :   if (tw == RgX_type_code(t_COMPLEX,t_INT) ||
     722             :       tw == RgX_type_code(t_COMPLEX,t_REAL))
     723           0 :   {
     724           0 :     long inv = (l >= 5 && signe(imag_i(gel(W,1+(l>>2))))==1) ? 1 : 0;
     725           0 :     fft(W+1, z+1, y+1, 1, l-1, inv);
     726             :   } else
     727          70 :     fft2(W+1, z+1, y+1, 1, l-1);
     728          70 :   return y;
     729             : }
     730             : 
     731             : GEN
     732          56 : FFT(GEN W, GEN x)
     733             : {
     734          56 :   if (!is_vec_t(typ(W))) pari_err_TYPE("fft",W);
     735          56 :   return FFT_i(W, x);
     736             : }
     737             : 
     738             : GEN
     739          56 : FFTinv(GEN W, GEN x)
     740             : {
     741          56 :   long l = lg(W), i;
     742             :   GEN w;
     743          56 :   if (!is_vec_t(typ(W))) pari_err_TYPE("fft",W);
     744          56 :   if (l==1 || ((l-1) & (l-2))) pari_err_DIM("fft");
     745          42 :   w = cgetg(l, t_VECSMALL); /* cf stackdummy */
     746          42 :   gel(w,1) = gel(W,1); /* w = gconj(W), faster */
     747        3773 :   for (i = 2; i < l; i++) gel(w, i) = gel(W, l-i+1);
     748          42 :   return FFT_i(w, x);
     749             : }
     750             : 
     751             : /* returns 1 if p has only real coefficients, 0 else */
     752             : static int
     753      961637 : isreal(GEN p)
     754             : {
     755             :   long i;
     756     4857616 :   for (i = lg(p)-1; i > 1; i--)
     757     4057615 :     if (typ(gel(p,i)) == t_COMPLEX) return 0;
     758      800001 :   return 1;
     759             : }
     760             : 
     761             : /* x non complex */
     762             : static GEN
     763      778140 : abs_update_r(GEN x, double *mu) {
     764      778140 :   GEN y = gtofp(x, DEFAULTPREC);
     765      778145 :   double ly = mydbllogr(y); if (ly < *mu) *mu = ly;
     766      778145 :   setabssign(y); return y;
     767             : }
     768             : /* return |x|, low accuracy. Set *mu = min(log(y), *mu) */
     769             : static GEN
     770     7999573 : abs_update(GEN x, double *mu) {
     771             :   GEN y, xr, yr;
     772             :   double ly;
     773     7999573 :   if (typ(x) != t_COMPLEX) return abs_update_r(x, mu);
     774     7234898 :   xr = gel(x,1);
     775     7234898 :   yr = gel(x,2);
     776     7234898 :   if (gequal0(xr)) return abs_update_r(yr,mu);
     777     7232906 :   if (gequal0(yr)) return abs_update_r(xr,mu);
     778             :   /* have to treat 0 specially: 0E-10 + 1e-20 = 0E-10 */
     779     7221652 :   xr = gtofp(xr, DEFAULTPREC);
     780     7222501 :   yr = gtofp(yr, DEFAULTPREC);
     781     7222528 :   y = sqrtr(addrr(sqrr(xr), sqrr(yr)));
     782     7222134 :   ly = mydbllogr(y); if (ly < *mu) *mu = ly;
     783     7222212 :   return y;
     784             : }
     785             : 
     786             : static void
     787      995811 : initdft(GEN *Omega, GEN *prim, long N, long Lmax, long bit)
     788             : {
     789      995811 :   long prec = nbits2prec(bit);
     790      995809 :   *Omega = grootsof1(Lmax, prec) + 1;
     791      995808 :   *prim = rootsof1u_cx(N, prec);
     792      995813 : }
     793             : 
     794             : static void
     795      493731 : parameters(GEN p, long *LMAX, double *mu, double *gamma,
     796             :            int polreal, double param, double param2)
     797             : {
     798             :   GEN q, pc, Omega, A, RU, prim, g, TWO;
     799      493731 :   long n = degpol(p), bit, NN, K, i, j, Lmax;
     800      493731 :   pari_sp av2, av = avma;
     801             : 
     802      493731 :   bit = gexpo(p) + (long)param2+8;
     803      683159 :   Lmax = 4; while (Lmax <= n) Lmax <<= 1;
     804      493732 :   NN = (long)(param*3.14)+1; if (NN < Lmax) NN = Lmax;
     805      493732 :   K = NN/Lmax; if (K & 1) K++;
     806      493732 :   NN = Lmax*K;
     807      493732 :   if (polreal) K = K/2+1;
     808             : 
     809      493732 :   initdft(&Omega, &prim, NN, Lmax, bit);
     810      493733 :   q = mygprec(p,bit) + 2;
     811      493731 :   A = cgetg(Lmax+1,t_VEC); A++;
     812      493732 :   pc= cgetg(Lmax+1,t_VEC); pc++;
     813     2955089 :   for (i=0; i <= n; i++) gel(pc,i)= gel(q,i);
     814      967058 :   for (   ; i<Lmax; i++) gel(pc,i) = gen_0;
     815             : 
     816      493732 :   *mu = pariINFINITY;
     817      493732 :   g = real_0_bit(-bit);
     818      493732 :   TWO = real2n(1, DEFAULTPREC);
     819      493729 :   av2 = avma;
     820      493729 :   RU = gen_1;
     821     1737111 :   for (i=0; i<K; i++)
     822             :   {
     823     1243378 :     if (i) {
     824      749652 :       GEN z = RU;
     825     3440499 :       for (j=1; j<n; j++)
     826             :       {
     827     2690844 :         gel(pc,j) = gmul(gel(q,j),z);
     828     2690816 :         z = gmul(z,RU); /* RU = prim^i, z=prim^(ij) */
     829             :       }
     830      749655 :       gel(pc,n) = gmul(gel(q,n),z);
     831             :     }
     832             : 
     833     1243375 :     fft(Omega,pc,A,1,Lmax,1);
     834     1243392 :     if (polreal && i>0 && i<K-1)
     835     1139183 :       for (j=0; j<Lmax; j++) g = addrr(g, divrr(TWO, abs_update(gel(A,j),mu)));
     836             :     else
     837     8103580 :       for (j=0; j<Lmax; j++) g = addrr(g, invr(abs_update(gel(A,j),mu)));
     838     1243084 :     RU = gmul(RU, prim);
     839     1243382 :     if (gc_needed(av,1))
     840             :     {
     841           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"parameters");
     842           0 :       gerepileall(av2,2, &g,&RU);
     843             :     }
     844             :   }
     845      493733 :   *gamma = mydbllog2r(divru(g,NN));
     846      493730 :   *LMAX = Lmax; set_avma(av);
     847      493730 : }
     848             : 
     849             : /* NN is a multiple of Lmax */
     850             : static void
     851      502080 : dft(GEN p, long k, long NN, long Lmax, long bit, GEN F, GEN H, long polreal)
     852             : {
     853             :   GEN Omega, q, qd, pc, pd, A, B, C, RU, aux, U, W, prim, prim2;
     854      502080 :   long n = degpol(p), i, j, K;
     855             :   pari_sp ltop;
     856             : 
     857      502080 :   initdft(&Omega, &prim, NN, Lmax, bit);
     858      502080 :   RU = cgetg(n+2,t_VEC) + 1;
     859             : 
     860      502079 :   K = NN/Lmax; if (polreal) K = K/2+1;
     861      502079 :   q = mygprec(p,bit);
     862      502077 :   qd = RgX_deriv(q);
     863             : 
     864      502077 :   A = cgetg(Lmax+1,t_VEC); A++;
     865      502077 :   B = cgetg(Lmax+1,t_VEC); B++;
     866      502077 :   C = cgetg(Lmax+1,t_VEC); C++;
     867      502077 :   pc = cgetg(Lmax+1,t_VEC); pc++;
     868      502078 :   pd = cgetg(Lmax+1,t_VEC); pd++;
     869     1017406 :   gel(pc,0) = gel(q,2);  for (i=n+1; i<Lmax; i++) gel(pc,i) = gen_0;
     870     1519479 :   gel(pd,0) = gel(qd,2); for (i=n;   i<Lmax; i++) gel(pd,i) = gen_0;
     871             : 
     872      502080 :   ltop = avma;
     873      502080 :   W = cgetg(k+1,t_VEC);
     874      502081 :   U = cgetg(k+1,t_VEC);
     875     1200637 :   for (i=1; i<=k; i++) gel(W,i) = gel(U,i) = gen_0;
     876             : 
     877      502080 :   gel(RU,0) = gen_1;
     878      502080 :   prim2 = gen_1;
     879     1529381 :   for (i=0; i<K; i++)
     880             :   {
     881     1027303 :     gel(RU,1) = prim2;
     882     4441790 :     for (j=1; j<n; j++) gel(RU,j+1) = gmul(gel(RU,j),prim2);
     883             :     /* RU[j] = prim^(ij)= prim2^j */
     884             : 
     885     4441766 :     for (j=1; j<n; j++) gel(pd,j) = gmul(gel(qd,j+2),gel(RU,j));
     886     1027255 :     fft(Omega,pd,A,1,Lmax,1);
     887     5469029 :     for (j=1; j<=n; j++) gel(pc,j) = gmul(gel(q,j+2),gel(RU,j));
     888     1027229 :     fft(Omega,pc,B,1,Lmax,1);
     889     7658105 :     for (j=0; j<Lmax; j++) gel(C,j) = ginv(gel(B,j));
     890     7658034 :     for (j=0; j<Lmax; j++) gel(B,j) = gmul(gel(A,j),gel(C,j));
     891     1027182 :     fft(Omega,B,A,1,Lmax,1);
     892     1027304 :     fft(Omega,C,B,1,Lmax,1);
     893             : 
     894     1027302 :     if (polreal) /* p has real coefficients */
     895             :     {
     896      796073 :       if (i>0 && i<K-1)
     897             :       {
     898      102472 :         for (j=1; j<=k; j++)
     899             :         {
     900       85881 :           gel(W,j) = gadd(gel(W,j), gshift(mulreal(gel(A,j+1),gel(RU,j+1)),1));
     901       85881 :           gel(U,j) = gadd(gel(U,j), gshift(mulreal(gel(B,j),gel(RU,j)),1));
     902             :         }
     903             :       }
     904             :       else
     905             :       {
     906     1827860 :         for (j=1; j<=k; j++)
     907             :         {
     908     1048393 :           gel(W,j) = gadd(gel(W,j), mulreal(gel(A,j+1),gel(RU,j+1)));
     909     1048363 :           gel(U,j) = gadd(gel(U,j), mulreal(gel(B,j),gel(RU,j)));
     910             :         }
     911             :       }
     912             :     }
     913             :     else
     914             :     {
     915      604742 :       for (j=1; j<=k; j++)
     916             :       {
     917      373516 :         gel(W,j) = gadd(gel(W,j), gmul(gel(A,j+1),gel(RU,j+1)));
     918      373509 :         gel(U,j) = gadd(gel(U,j), gmul(gel(B,j),gel(RU,j)));
     919             :       }
     920             :     }
     921     1027284 :     prim2 = gmul(prim2,prim);
     922     1027291 :     gerepileall(ltop,3, &W,&U,&prim2);
     923             :   }
     924             : 
     925     1200629 :   for (i=1; i<=k; i++)
     926             :   {
     927      698556 :     aux=gel(W,i);
     928     1095504 :     for (j=1; j<i; j++) aux = gadd(aux, gmul(gel(W,i-j),gel(F,k+2-j)));
     929      698552 :     gel(F,k+2-i) = gdivgs(aux,-i*NN);
     930             :   }
     931     1200624 :   for (i=0; i<k; i++)
     932             :   {
     933      698552 :     aux=gel(U,k-i);
     934     1095503 :     for (j=1+i; j<k; j++) aux = gadd(aux,gmul(gel(F,2+j),gel(U,j-i)));
     935      698549 :     gel(H,i+2) = gdivgu(aux,NN);
     936             :   }
     937      502072 : }
     938             : 
     939             : #define NEWTON_MAX 10
     940             : static GEN
     941     2459693 : refine_H(GEN F, GEN G, GEN HH, long bit, long Sbit)
     942             : {
     943     2459693 :   GEN H = HH, D, aux;
     944     2459693 :   pari_sp ltop = avma;
     945             :   long error, i, bit1, bit2;
     946             : 
     947     2459693 :   D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F)); error = gexpo(D);
     948     2459664 :   bit2 = bit + Sbit;
     949     4495310 :   for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
     950             :   {
     951     2035643 :     if (gc_needed(ltop,1))
     952             :     {
     953           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"refine_H");
     954           0 :       gerepileall(ltop,2, &D,&H);
     955             :     }
     956     2035643 :     bit1 = -error + Sbit;
     957     2035643 :     aux = RgX_mul(mygprec(H,bit1), mygprec(D,bit1));
     958     2035654 :     aux = RgX_rem(mygprec(aux,bit1), mygprec(F,bit1));
     959             : 
     960     2035666 :     bit1 = -error*2 + Sbit; if (bit1 > bit2) bit1 = bit2;
     961     2035666 :     H = RgX_add(mygprec(H,bit1), aux);
     962     2035598 :     D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F));
     963     2035648 :     error = gexpo(D); if (error < -bit1) error = -bit1;
     964             :   }
     965     2459667 :   if (error > -bit/2) return NULL; /* FAIL */
     966     2459341 :   return gerepilecopy(ltop,H);
     967             : }
     968             : 
     969             : /* return 0 if fails, 1 else */
     970             : static long
     971      502075 : refine_F(GEN p, GEN *F, GEN *G, GEN H, long bit, double gamma)
     972             : {
     973      502075 :   GEN f0, FF, GG, r, HH = H;
     974      502075 :   long error, i, bit1 = 0, bit2, Sbit, Sbit2,  enh, normF, normG, n = degpol(p);
     975      502075 :   pari_sp av = avma;
     976             : 
     977      502075 :   FF = *F; GG = RgX_divrem(p, FF, &r);
     978      502079 :   error = gexpo(r); if (error <= -bit) error = 1-bit;
     979      502079 :   normF = gexpo(FF);
     980      502079 :   normG = gexpo(GG);
     981      502080 :   enh = gexpo(H); if (enh < 0) enh = 0;
     982      502080 :   Sbit = normF + 2*normG + enh + (long)(4.*log2((double)n)+gamma) + 1;
     983      502080 :   Sbit2 = enh + 2*(normF+normG) + (long)(2.*gamma+5.*log2((double)n)) + 1;
     984      502080 :   bit2 = bit + Sbit;
     985     2961438 :   for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
     986             :   {
     987     2459698 :     if (bit1 == bit2 && i >= 2) { Sbit += n; Sbit2 += n; bit2 += n; }
     988     2459698 :     if (gc_needed(av,1))
     989             :     {
     990           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"refine_F");
     991           0 :       gerepileall(av,4, &FF,&GG,&r,&HH);
     992             :     }
     993             : 
     994     2459698 :     bit1 = -error + Sbit2;
     995     2459698 :     HH = refine_H(mygprec(FF,bit1), mygprec(GG,bit1), mygprec(HH,bit1),
     996             :                   1-error, Sbit2);
     997     2459714 :     if (!HH) return 0; /* FAIL */
     998             : 
     999     2459388 :     bit1 = -error + Sbit;
    1000     2459388 :     r = RgX_mul(mygprec(HH,bit1), mygprec(r,bit1));
    1001     2459353 :     f0 = RgX_rem(mygprec(r,bit1), mygprec(FF,bit1));
    1002             : 
    1003     2459369 :     bit1 = -2*error + Sbit; if (bit1 > bit2) bit1 = bit2;
    1004     2459369 :     FF = gadd(mygprec(FF,bit1),f0);
    1005             : 
    1006     2459342 :     bit1 = -3*error + Sbit; if (bit1 > bit2) bit1 = bit2;
    1007     2459342 :     GG = RgX_divrem(mygprec(p,bit1), mygprec(FF,bit1), &r);
    1008     2459361 :     error = gexpo(r); if (error < -bit1) error = -bit1;
    1009             :   }
    1010      501740 :   if (error>-bit) return 0; /* FAIL */
    1011      493719 :   *F = FF; *G = GG; return 1;
    1012             : }
    1013             : 
    1014             : /* returns F and G from the unit circle U such that |p-FG|<2^(-bit) |cd|,
    1015             : where cd is the leading coefficient of p */
    1016             : static void
    1017      493732 : split_fromU(GEN p, long k, double delta, long bit,
    1018             :             GEN *F, GEN *G, double param, double param2)
    1019             : {
    1020             :   GEN pp, FF, GG, H;
    1021      493732 :   long n = degpol(p), NN, bit2, Lmax;
    1022      493732 :   int polreal = isreal(p);
    1023             :   pari_sp ltop;
    1024             :   double mu, gamma;
    1025             : 
    1026      493732 :   pp = gdiv(p, gel(p,2+n));
    1027      493731 :   parameters(pp, &Lmax,&mu,&gamma, polreal,param,param2);
    1028             : 
    1029      493730 :   H  = cgetg(k+2,t_POL); H[1] = p[1];
    1030      493731 :   FF = cgetg(k+3,t_POL); FF[1]= p[1];
    1031      493732 :   gel(FF,k+2) = gen_1;
    1032             : 
    1033      493732 :   NN = (long)(0.5/delta); NN |= 1; if (NN < 2) NN = 2;
    1034      493732 :   NN *= Lmax; ltop = avma;
    1035             :   for(;;)
    1036             :   {
    1037      502079 :     bit2 = (long)(((double)NN*delta-mu)/M_LN2) + gexpo(pp) + 8;
    1038      502080 :     dft(pp, k, NN, Lmax, bit2, FF, H, polreal);
    1039      502074 :     if (refine_F(pp,&FF,&GG,H,bit,gamma)) break;
    1040        8347 :     NN <<= 1; set_avma(ltop);
    1041             :   }
    1042      493732 :   *G = gmul(GG,gel(p,2+n)); *F = FF;
    1043      493731 : }
    1044             : 
    1045             : static void
    1046      493733 : optimize_split(GEN p, long k, double delta, long bit,
    1047             :             GEN *F, GEN *G, double param, double param2)
    1048             : {
    1049      493733 :   long n = degpol(p);
    1050             :   GEN FF, GG;
    1051             : 
    1052      493732 :   if (k <= n/2)
    1053      382825 :     split_fromU(p,k,delta,bit,F,G,param,param2);
    1054             :   else
    1055             :   {
    1056      110907 :     split_fromU(RgX_recip_i(p),n-k,delta,bit,&FF,&GG,param,param2);
    1057      110908 :     *F = RgX_recip_i(GG);
    1058      110908 :     *G = RgX_recip_i(FF);
    1059             :   }
    1060      493732 : }
    1061             : 
    1062             : /********************************************************************/
    1063             : /**                                                                **/
    1064             : /**               SEARCH FOR SEPARATING CIRCLE                     **/
    1065             : /**                                                                **/
    1066             : /********************************************************************/
    1067             : 
    1068             : /* return p(2^e*x) *2^(-n*e) */
    1069             : static void
    1070           0 : scalepol2n(GEN p, long e)
    1071             : {
    1072           0 :   long i,n=lg(p)-1;
    1073           0 :   for (i=2; i<=n; i++) gel(p,i) = gmul2n(gel(p,i),(i-n)*e);
    1074           0 : }
    1075             : 
    1076             : /* returns p(x/R)*R^n; assume R is at the correct accuracy */
    1077             : static GEN
    1078     4288511 : scalepol(GEN p, GEN R, long bit)
    1079     4288511 : { return RgX_rescale(mygprec(p, bit), R); }
    1080             : 
    1081             : /* return (conj(a)X-1)^n * p[ (X-a) / (conj(a)X-1) ] */
    1082             : static GEN
    1083     1403832 : conformal_basecase(GEN p, GEN a)
    1084             : {
    1085             :   GEN z, r, ma, ca;
    1086     1403832 :   long i, n = degpol(p);
    1087             :   pari_sp av;
    1088             : 
    1089     1403831 :   if (n <= 0) return p;
    1090     1403831 :   ma = gneg(a); ca = conj_i(a);
    1091     1403828 :   av = avma;
    1092     1403828 :   z = deg1pol_shallow(ca, gen_m1, 0);
    1093     1403830 :   r = scalarpol_shallow(gel(p,2+n), 0);
    1094     3641446 :   for (i=n-1; ; i--)
    1095             :   {
    1096     3641446 :     r = RgX_addmulXn_shallow(r, gmul(ma,r), 1); /* r *= (X - a) */
    1097     3641421 :     r = gadd(r, gmul(z, gel(p,2+i)));
    1098     3641419 :     if (i == 0) return gerepileupto(av, r);
    1099     2237599 :     z = RgX_addmulXn_shallow(gmul(z,ca), gneg(z), 1); /* z *= conj(a)X - 1 */
    1100     2237621 :     if (gc_needed(av,2))
    1101             :     {
    1102           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"conformal_pol (%ld/%ld)",n-i, n);
    1103           0 :       gerepileall(av,2, &r,&z);
    1104             :     }
    1105             :   }
    1106             : }
    1107             : static GEN
    1108     1403950 : conformal_pol(GEN p, GEN a)
    1109             : {
    1110     1403950 :   pari_sp av = avma;
    1111     1403950 :   long d, nR, n = degpol(p), v;
    1112             :   GEN Q, R, S, T;
    1113     1403951 :   if (n < 35) return conformal_basecase(p, a);
    1114         119 :   d = (n+1) >> 1; v = varn(p);
    1115         119 :   Q = RgX_shift_shallow(p, -d);
    1116         119 :   R = RgXn_red_shallow(p, d);
    1117         119 :   Q = conformal_pol(Q, a);
    1118         119 :   R = conformal_pol(R, a);
    1119         119 :   S = gpowgs(deg1pol_shallow(gen_1, gneg(a), v), d);
    1120         119 :   T = RgX_recip_i(S);
    1121         119 :   if (typ(a) == t_COMPLEX) T = gconj(T);
    1122         119 :   if (odd(d)) T = RgX_neg(T);
    1123             :   /* S = (X - a)^d, T = (conj(a) X - 1)^d */
    1124         119 :   nR = n - degpol(R) - d; /* >= 0 */
    1125         119 :   if (nR) T = RgX_mul(T, gpowgs(deg1pol_shallow(gconj(a), gen_m1, v), nR));
    1126         119 :   return gerepileupto(av, RgX_add(RgX_mul(Q, S), RgX_mul(R, T)));
    1127             : }
    1128             : 
    1129             : static const double UNDEF = -100000.;
    1130             : 
    1131             : static double
    1132      493725 : logradius(double *radii, GEN p, long k, double aux, double *delta)
    1133             : {
    1134      493725 :   long i, n = degpol(p);
    1135             :   double lrho, lrmin, lrmax;
    1136      493726 :   if (k > 1)
    1137             :   {
    1138      282079 :     i = k-1; while (i>0 && radii[i] == UNDEF) i--;
    1139      206798 :     lrmin = logpre_modulus(p,k,aux, radii[i], radii[k]);
    1140             :   }
    1141             :   else /* k=1 */
    1142      286928 :     lrmin = logmin_modulus(p,aux);
    1143      493728 :   radii[k] = lrmin;
    1144             : 
    1145      493728 :   if (k+1<n)
    1146             :   {
    1147      591054 :     i = k+2; while (i<=n && radii[i] == UNDEF) i++;
    1148      399319 :     lrmax = logpre_modulus(p,k+1,aux, radii[k+1], radii[i]);
    1149             :   }
    1150             :   else /* k+1=n */
    1151       94409 :     lrmax = logmax_modulus(p,aux);
    1152      493727 :   radii[k+1] = lrmax;
    1153             : 
    1154      493727 :   lrho = radii[k];
    1155      823490 :   for (i=k-1; i>=1; i--)
    1156             :   {
    1157      329763 :     if (radii[i] == UNDEF || radii[i] > lrho)
    1158      241998 :       radii[i] = lrho;
    1159             :     else
    1160       87765 :       lrho = radii[i];
    1161             :   }
    1162      493727 :   lrho = radii[k+1];
    1163     1637848 :   for (i=k+1; i<=n; i++)
    1164             :   {
    1165     1144121 :     if (radii[i] == UNDEF || radii[i] < lrho)
    1166      567642 :       radii[i] = lrho;
    1167             :     else
    1168      576479 :       lrho = radii[i];
    1169             :   }
    1170      493727 :   *delta = (lrmax - lrmin) / 2;
    1171      493727 :   if (*delta > 1.) *delta = 1.;
    1172      493727 :   return (lrmin + lrmax) / 2;
    1173             : }
    1174             : 
    1175             : static void
    1176      493727 : update_radius(long n, double *radii, double lrho, double *par, double *par2)
    1177             : {
    1178      493727 :   double t, param = 0., param2 = 0.;
    1179             :   long i;
    1180     2461252 :   for (i=1; i<=n; i++)
    1181             :   {
    1182     1967555 :     radii[i] -= lrho;
    1183     1967555 :     t = fabs(rtodbl( invr(subsr(1, dblexp(radii[i]))) ));
    1184     1967525 :     param += t; if (t > 1.) param2 += log2(t);
    1185             :   }
    1186      493697 :   *par = param; *par2 = param2;
    1187      493697 : }
    1188             : 
    1189             : /* apply the conformal mapping then split from U */
    1190             : static void
    1191      467906 : conformal_mapping(double *radii, GEN ctr, GEN p, long k, long bit,
    1192             :                   double aux, GEN *F,GEN *G)
    1193             : {
    1194      467906 :   long bit2, n = degpol(p), i;
    1195      467906 :   pari_sp ltop = avma, av;
    1196             :   GEN q, FF, GG, a, R;
    1197             :   double lrho, delta, param, param2;
    1198             :   /* n * (2.*log2(2.732)+log2(1.5)) + 1 */
    1199      467906 :   bit2 = bit + (long)(n*3.4848775) + 1;
    1200      467906 :   a = sqrtr_abs( utor(3, precdbl(MEDDEFAULTPREC)) );
    1201      467906 :   a = divrs(a, -6);
    1202      467906 :   a = gmul(mygprec(a,bit2), mygprec(ctr,bit2)); /* a = -ctr/2sqrt(3) */
    1203             : 
    1204      467906 :   av = avma;
    1205      467906 :   q = conformal_pol(mygprec(p,bit2), a);
    1206     2288670 :   for (i=1; i<=n; i++)
    1207     1820771 :     if (radii[i] != UNDEF) /* update array radii */
    1208             :     {
    1209     1540878 :       pari_sp av2 = avma;
    1210     1540878 :       GEN t, r = dblexp(radii[i]), r2 = sqrr(r);
    1211             :       /* 2(r^2 - 1) / (r^2 - 3(r-1)) */
    1212     1540815 :       t = divrr(shiftr((subrs(r2,1)),1), subrr(r2, mulur(3,subrs(r,1))));
    1213     1540873 :       radii[i] = mydbllogr(addsr(1,t)) / 2;
    1214     1540853 :       set_avma(av2);
    1215             :     }
    1216      467899 :   lrho = logradius(radii, q,k,aux/10., &delta);
    1217      467901 :   update_radius(n, radii, lrho, &param, &param2);
    1218             : 
    1219      467901 :   bit2 += (long)(n * fabs(lrho)/M_LN2 + 1.);
    1220      467901 :   R = mygprec(dblexp(-lrho), bit2);
    1221      467900 :   q = scalepol(q,R,bit2);
    1222      467906 :   gerepileall(av,2, &q,&R);
    1223             : 
    1224      467907 :   optimize_split(q,k,delta,bit2,&FF,&GG,param,param2);
    1225      467906 :   bit2 += n; R = invr(R);
    1226      467902 :   FF = scalepol(FF,R,bit2);
    1227      467903 :   GG = scalepol(GG,R,bit2);
    1228             : 
    1229      467905 :   a = mygprec(a,bit2);
    1230      467903 :   FF = conformal_pol(FF,a);
    1231      467905 :   GG = conformal_pol(GG,a);
    1232             : 
    1233      467907 :   a = invr(subsr(1, gnorm(a)));
    1234      467899 :   FF = RgX_Rg_mul(FF, powru(a,k));
    1235      467904 :   GG = RgX_Rg_mul(GG, powru(a,n-k));
    1236             : 
    1237      467905 :   *F = mygprec(FF,bit+n);
    1238      467905 :   *G = mygprec(GG,bit+n); gerepileall(ltop,2, F,G);
    1239      467907 : }
    1240             : 
    1241             : /* split p, this time without scaling. returns in F and G two polynomials
    1242             :  * such that |p-FG|< 2^(-bit)|p| */
    1243             : static void
    1244      493733 : split_2(GEN p, long bit, GEN ctr, double thickness, GEN *F, GEN *G)
    1245             : {
    1246             :   GEN q, FF, GG, R;
    1247             :   double aux, delta, param, param2;
    1248      493733 :   long n = degpol(p), i, j, k, bit2;
    1249             :   double lrmin, lrmax, lrho, *radii;
    1250             : 
    1251      493734 :   radii = (double*) stack_malloc_align((n+1) * sizeof(double), sizeof(double));
    1252             : 
    1253     1473896 :   for (i=2; i<n; i++) radii[i] = UNDEF;
    1254      493733 :   aux = thickness/(double)(4 * n);
    1255      493733 :   lrmin = logmin_modulus(p, aux);
    1256      493729 :   lrmax = logmax_modulus(p, aux);
    1257      493728 :   radii[1] = lrmin;
    1258      493728 :   radii[n] = lrmax;
    1259      493728 :   i = 1; j = n;
    1260      493728 :   lrho = (lrmin + lrmax) / 2;
    1261      493728 :   k = dual_modulus(p, lrho, aux, 1);
    1262      493733 :   if (5*k < n || (n < 2*k && 5*k < 4*n))
    1263       77988 :     { lrmax = lrho; j=k+1; radii[j] = lrho; }
    1264             :   else
    1265      415745 :     { lrmin = lrho; i=k;   radii[i] = lrho; }
    1266     1011957 :   while (j > i+1)
    1267             :   {
    1268      518223 :     if (i+j == n+1)
    1269      372383 :       lrho = (lrmin + lrmax) / 2;
    1270             :     else
    1271             :     {
    1272      145840 :       double kappa = 2. - log(1. + minss(i,n-j)) / log(1. + minss(j,n-i));
    1273      145840 :       if (i+j < n+1) lrho = lrmax * kappa + lrmin;
    1274      116270 :       else           lrho = lrmin * kappa + lrmax;
    1275      145840 :       lrho /= 1+kappa;
    1276             :     }
    1277      518223 :     aux = (lrmax - lrmin) / (4*(j-i));
    1278      518223 :     k = dual_modulus(p, lrho, aux, minss(i,n+1-j));
    1279      518224 :     if (k-i < j-k-1 || (k-i == j-k-1 && 2*k > n))
    1280      386911 :       { lrmax = lrho; j=k+1; radii[j] = lrho - aux; }
    1281             :     else
    1282      131313 :       { lrmin = lrho; i=k;   radii[i] = lrho + aux; }
    1283             :   }
    1284      493734 :   aux = lrmax - lrmin;
    1285             : 
    1286      493734 :   if (ctr)
    1287             :   {
    1288      467907 :     lrho = (lrmax + lrmin) / 2;
    1289     2288713 :     for (i=1; i<=n; i++)
    1290     1820806 :       if (radii[i] != UNDEF) radii[i] -= lrho;
    1291             : 
    1292      467907 :     bit2 = bit + (long)(n * fabs(lrho)/M_LN2 + 1.);
    1293      467907 :     R = mygprec(dblexp(-lrho), bit2);
    1294      467906 :     q = scalepol(p,R,bit2);
    1295      467906 :     conformal_mapping(radii, ctr, q, k, bit2, aux, &FF, &GG);
    1296             :   }
    1297             :   else
    1298             :   {
    1299       25827 :     lrho = logradius(radii, p, k, aux/10., &delta);
    1300       25827 :     update_radius(n, radii, lrho, &param, &param2);
    1301             : 
    1302       25827 :     bit2 = bit + (long)(n * fabs(lrho)/M_LN2 + 1.);
    1303       25827 :     R = mygprec(dblexp(-lrho), bit2);
    1304       25827 :     q = scalepol(p,R,bit2);
    1305       25827 :     optimize_split(q, k, delta, bit2, &FF, &GG, param, param2);
    1306             :   }
    1307      493733 :   bit  += n;
    1308      493733 :   bit2 += n; R = invr(mygprec(R,bit2));
    1309      493731 :   *F = mygprec(scalepol(FF,R,bit2), bit);
    1310      493732 :   *G = mygprec(scalepol(GG,R,bit2), bit);
    1311      493733 : }
    1312             : 
    1313             : /* procedure corresponding to steps 5,6,.. page 44 in RR n. 1852 */
    1314             : /* put in F and G two polynomial such that |p-FG|<2^(-bit)|p|
    1315             :  * where the maximum modulus of the roots of p is <=1.
    1316             :  * Assume sum of roots is 0. */
    1317             : static void
    1318      467905 : split_1(GEN p, long bit, GEN *F, GEN *G)
    1319             : {
    1320      467905 :   long i, imax, n = degpol(p), polreal = isreal(p), ep = gexpo(p), bit2 = bit+n;
    1321             :   GEN ctr, q, qq, FF, GG, v, gr, r, newq;
    1322             :   double lrmin, lrmax, lthick;
    1323      467906 :   const double LOG3 = 1.098613;
    1324             : 
    1325      467906 :   lrmax = logmax_modulus(p, 0.01);
    1326      467906 :   gr = mygprec(dblexp(-lrmax), bit2);
    1327      467906 :   q = scalepol(p,gr,bit2);
    1328             : 
    1329      467905 :   bit2 = bit + gexpo(q) - ep + (long)((double)n*2.*log2(3.)+1);
    1330      467905 :   v = cgetg(5,t_VEC);
    1331      467905 :   gel(v,1) = gen_2;
    1332      467905 :   gel(v,2) = gen_m2;
    1333      467905 :   gel(v,3) = mkcomplex(gen_0, gel(v,1));
    1334      467905 :   gel(v,4) = mkcomplex(gen_0, gel(v,2));
    1335      467906 :   q = mygprec(q,bit2); lthick = 0;
    1336      467907 :   newq = ctr = NULL; /* -Wall */
    1337      467907 :   imax = polreal? 3: 4;
    1338      840643 :   for (i=1; i<=imax; i++)
    1339             :   {
    1340      834089 :     qq = RgX_translate(q, gel(v,i));
    1341      834096 :     lrmin = logmin_modulus(qq,0.05);
    1342      834082 :     if (LOG3 > lrmin + lthick)
    1343             :     {
    1344      821928 :       double lquo = logmax_modulus(qq,0.05) - lrmin;
    1345      821932 :       if (lquo > lthick) { lthick = lquo; newq = qq; ctr = gel(v,i); }
    1346             :     }
    1347      834086 :     if (lthick > M_LN2) break;
    1348      423841 :     if (polreal && i==2 && lthick > LOG3 - M_LN2) break;
    1349             :   }
    1350      467904 :   bit2 = bit + gexpo(newq) - ep + (long)(n*LOG3/M_LN2 + 1);
    1351      467906 :   split_2(newq, bit2, ctr, lthick, &FF, &GG);
    1352      467906 :   r = gneg(mygprec(ctr,bit2));
    1353      467905 :   FF = RgX_translate(FF,r);
    1354      467906 :   GG = RgX_translate(GG,r);
    1355             : 
    1356      467907 :   gr = invr(gr); bit2 = bit - ep + gexpo(FF)+gexpo(GG);
    1357      467907 :   *F = scalepol(FF,gr,bit2);
    1358      467906 :   *G = scalepol(GG,gr,bit2);
    1359      467906 : }
    1360             : 
    1361             : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|,
    1362             : where the maximum modulus of the roots of p is < 0.5 */
    1363             : static int
    1364      468228 : split_0_2(GEN p, long bit, GEN *F, GEN *G)
    1365             : {
    1366             :   GEN q, b;
    1367      468228 :   long n = degpol(p), k, bit2, eq;
    1368      468228 :   double aux0 = dbllog2(gel(p,n+2)); /* != -oo */
    1369      468229 :   double aux1 = dbllog2(gel(p,n+1)), aux;
    1370             : 
    1371      468229 :   if (aux1 == -pariINFINITY) /* p1 = 0 */
    1372        9886 :     aux = 0;
    1373             :   else
    1374             :   {
    1375      458343 :     aux = aux1 - aux0; /* log2(p1/p0) */
    1376             :     /* beware double overflow */
    1377      458343 :     if (aux >= 0 && (aux > 1e4 || exp2(aux) > 2.5*n)) return 0;
    1378      458343 :     aux = (aux < -300)? 0.: n*log2(1 + exp2(aux)/(double)n);
    1379             :   }
    1380      468229 :   bit2 = bit+1 + (long)(log2((double)n) + aux);
    1381      468229 :   q = mygprec(p,bit2);
    1382      468229 :   if (aux1 == -pariINFINITY) b = NULL;
    1383             :   else
    1384             :   {
    1385      458343 :     b = gdivgs(gdiv(gel(q,n+1),gel(q,n+2)),-n);
    1386      458344 :     q = RgX_translate(q,b);
    1387             :   }
    1388      468230 :   gel(q,n+1) = gen_0; eq = gexpo(q);
    1389      468229 :   k = 0;
    1390      468783 :   while (k <= n/2 && (- gexpo(gel(q,k+2)) > bit2 + 2*(n-k) + eq
    1391      468655 :                       || gequal0(gel(q,k+2)))) k++;
    1392      468229 :   if (k > 0)
    1393             :   {
    1394         324 :     if (k > n/2) k = n/2;
    1395         324 :     bit2 += k<<1;
    1396         324 :     *F = pol_xn(k, 0);
    1397         324 :     *G = RgX_shift_shallow(q, -k);
    1398             :   }
    1399             :   else
    1400             :   {
    1401      467905 :     split_1(q,bit2,F,G);
    1402      467905 :     bit2 = bit + gexpo(*F) + gexpo(*G) - gexpo(p) + (long)aux+1;
    1403      467907 :     *F = mygprec(*F,bit2);
    1404             :   }
    1405      468230 :   *G = mygprec(*G,bit2);
    1406      468229 :   if (b)
    1407             :   {
    1408      458343 :     GEN mb = mygprec(gneg(b), bit2);
    1409      458343 :     *F = RgX_translate(*F, mb);
    1410      458344 :     *G = RgX_translate(*G, mb);
    1411             :   }
    1412      468230 :   return 1;
    1413             : }
    1414             : 
    1415             : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|.
    1416             :  * Assume max_modulus(p) < 2 */
    1417             : static void
    1418      468228 : split_0_1(GEN p, long bit, GEN *F, GEN *G)
    1419             : {
    1420             :   GEN FF, GG;
    1421             :   long n, bit2, normp;
    1422             : 
    1423      468228 :   if  (split_0_2(p,bit,F,G)) return;
    1424             : 
    1425           0 :   normp = gexpo(p);
    1426           0 :   scalepol2n(p,2); /* p := 4^(-n) p(4*x) */
    1427           0 :   n = degpol(p); bit2 = bit + 2*n + gexpo(p) - normp;
    1428           0 :   split_1(mygprec(p,bit2), bit2,&FF,&GG);
    1429           0 :   scalepol2n(FF,-2);
    1430           0 :   scalepol2n(GG,-2); bit2 = bit + gexpo(FF) + gexpo(GG) - normp;
    1431           0 :   *F = mygprec(FF,bit2);
    1432           0 :   *G = mygprec(GG,bit2);
    1433             : }
    1434             : 
    1435             : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P| */
    1436             : static void
    1437      494056 : split_0(GEN p, long bit, GEN *F, GEN *G)
    1438             : {
    1439      494056 :   const double LOG1_9 = 0.6418539;
    1440      494056 :   long n = degpol(p), k = 0;
    1441             :   GEN q;
    1442             : 
    1443      494056 :   while (gexpo(gel(p,k+2)) < -bit && k <= n/2) k++;
    1444      494056 :   if (k > 0)
    1445             :   {
    1446           0 :     if (k > n/2) k = n/2;
    1447           0 :     *F = pol_xn(k, 0);
    1448           0 :     *G = RgX_shift_shallow(p, -k);
    1449             :   }
    1450             :   else
    1451             :   {
    1452      494056 :     double lr = logmax_modulus(p, 0.05);
    1453      494054 :     if (lr < LOG1_9) split_0_1(p, bit, F, G);
    1454             :     else
    1455             :     {
    1456      436616 :       q = RgX_recip_i(p);
    1457      436617 :       lr = logmax_modulus(q,0.05);
    1458      436617 :       if (lr < LOG1_9)
    1459             :       {
    1460      410790 :         split_0_1(q, bit, F, G);
    1461      410791 :         *F = RgX_recip_i(*F);
    1462      410791 :         *G = RgX_recip_i(*G);
    1463             :       }
    1464             :       else
    1465       25827 :         split_2(p,bit,NULL, 1.2837,F,G);
    1466             :     }
    1467             :   }
    1468      494057 : }
    1469             : 
    1470             : /********************************************************************/
    1471             : /**                                                                **/
    1472             : /**                ERROR ESTIMATE FOR THE ROOTS                    **/
    1473             : /**                                                                **/
    1474             : /********************************************************************/
    1475             : 
    1476             : static GEN
    1477     1897661 : root_error(long n, long k, GEN roots_pol, long err, GEN shatzle)
    1478             : {
    1479     1897661 :   GEN rho, d, eps, epsbis, eps2, aux, rap = NULL;
    1480             :   long i, j;
    1481             : 
    1482     1897661 :   d = cgetg(n+1,t_VEC);
    1483    12088141 :   for (i=1; i<=n; i++)
    1484             :   {
    1485    10190837 :     if (i!=k)
    1486             :     {
    1487     8293353 :       aux = gsub(gel(roots_pol,i), gel(roots_pol,k));
    1488     8293012 :       gel(d,i) = gabs(mygprec(aux,31), DEFAULTPREC);
    1489             :     }
    1490             :   }
    1491     1897304 :   rho = gabs(mygprec(gel(roots_pol,k),31), DEFAULTPREC);
    1492     1897647 :   if (expo(rho) < 0) rho = real_1(DEFAULTPREC);
    1493     1897647 :   eps = mulrr(rho, shatzle);
    1494     1897586 :   aux = shiftr(powru(rho,n), err);
    1495             : 
    1496     5759484 :   for (j=1; j<=2 || (j<=5 && cmprr(rap, dbltor(1.2)) > 0); j++)
    1497             :   {
    1498     3861990 :     GEN prod = NULL; /* 1. */
    1499     3861990 :     long m = n;
    1500     3861990 :     epsbis = mulrr(eps, dbltor(1.25));
    1501    26426160 :     for (i=1; i<=n; i++)
    1502             :     {
    1503    22563906 :       if (i != k && cmprr(gel(d,i),epsbis) > 0)
    1504             :       {
    1505    18662414 :         GEN dif = subrr(gel(d,i),eps);
    1506    18660191 :         prod = prod? mulrr(prod, dif): dif;
    1507    18661635 :         m--;
    1508             :       }
    1509             :     }
    1510     3862254 :     eps2 = prod? divrr(aux, prod): aux;
    1511     3862024 :     if (m > 1) eps2 = sqrtnr(shiftr(eps2, 2*m-2), m);
    1512     3862024 :     rap = divrr(eps,eps2); eps = eps2;
    1513             :   }
    1514     1897426 :   return eps;
    1515             : }
    1516             : 
    1517             : /* round a complex or real number x to an absolute value of 2^(-bit) */
    1518             : static GEN
    1519     4297139 : mygprec_absolute(GEN x, long bit)
    1520             : {
    1521             :   long e;
    1522             :   GEN y;
    1523             : 
    1524     4297139 :   switch(typ(x))
    1525             :   {
    1526     2949962 :     case t_REAL:
    1527     2949962 :       e = expo(x) + bit;
    1528     2949962 :       return (e <= 0 || !signe(x))? real_0_bit(-bit): rtor(x, nbits2prec(e));
    1529     1217087 :     case t_COMPLEX:
    1530     1217087 :       if (gexpo(gel(x,2)) < -bit) return mygprec_absolute(gel(x,1),bit);
    1531     1182497 :       y = cgetg(3,t_COMPLEX);
    1532     1182503 :       gel(y,1) = mygprec_absolute(gel(x,1),bit);
    1533     1182511 :       gel(y,2) = mygprec_absolute(gel(x,2),bit);
    1534     1182511 :       return y;
    1535      130090 :     default: return x;
    1536             :   }
    1537             : }
    1538             : 
    1539             : static long
    1540      530511 : a_posteriori_errors(GEN p, GEN roots_pol, long err)
    1541             : {
    1542      530511 :   long i, n = degpol(p), e_max = -(long)EXPOBITS;
    1543             :   GEN sigma, shatzle;
    1544             : 
    1545      530511 :   err += (long)log2((double)n) + 1;
    1546      530511 :   if (err > -2) return 0;
    1547      530511 :   sigma = real2n(-err, LOWDEFAULTPREC);
    1548             :   /*  2 / ((s - 1)^(1/n) - 1) */
    1549      530511 :   shatzle = divur(2, subrs(sqrtnr(subrs(sigma,1),n), 1));
    1550     2428170 :   for (i=1; i<=n; i++)
    1551             :   {
    1552     1897659 :     pari_sp av = avma;
    1553     1897659 :     GEN x = root_error(n,i,roots_pol,err,shatzle);
    1554     1897421 :     long e = gexpo(x);
    1555     1897474 :     set_avma(av); if (e > e_max) e_max = e;
    1556     1897565 :     gel(roots_pol,i) = mygprec_absolute(gel(roots_pol,i), -e);
    1557             :   }
    1558      530511 :   return e_max;
    1559             : }
    1560             : 
    1561             : /********************************************************************/
    1562             : /**                                                                **/
    1563             : /**                           MAIN                                 **/
    1564             : /**                                                                **/
    1565             : /********************************************************************/
    1566             : static GEN
    1567     1602259 : append_clone(GEN r, GEN a) { a = gclone(a); vectrunc_append(r, a); return a; }
    1568             : 
    1569             : /* put roots in placeholder roots_pol so that |P - L_1...L_n| < 2^(-bit)|P|
    1570             :  * returns prod (x-roots_pol[i]) */
    1571             : static GEN
    1572     1518616 : split_complete(GEN p, long bit, GEN roots_pol)
    1573             : {
    1574     1518616 :   long n = degpol(p);
    1575             :   pari_sp ltop;
    1576             :   GEN p1, F, G, a, b, m1, m2;
    1577             : 
    1578     1518616 :   if (n == 1)
    1579             :   {
    1580      446854 :     a = gneg_i(gdiv(gel(p,2), gel(p,3)));
    1581      446862 :     (void)append_clone(roots_pol,a); return p;
    1582             :   }
    1583     1071762 :   ltop = avma;
    1584     1071762 :   if (n == 2)
    1585             :   {
    1586      577706 :     F = gsub(gsqr(gel(p,3)), gmul2n(gmul(gel(p,2),gel(p,4)), 2));
    1587      577697 :     F = gsqrt(F, nbits2prec(bit));
    1588      577701 :     p1 = ginv(gmul2n(gel(p,4),1));
    1589      577698 :     a = gneg_i(gmul(gadd(F,gel(p,3)), p1));
    1590      577707 :     b =        gmul(gsub(F,gel(p,3)), p1);
    1591      577706 :     a = append_clone(roots_pol,a);
    1592      577708 :     b = append_clone(roots_pol,b); set_avma(ltop);
    1593      577708 :     a = mygprec(a, 3*bit);
    1594      577706 :     b = mygprec(b, 3*bit);
    1595      577707 :     return gmul(gel(p,4), mkpoln(3, gen_1, gneg(gadd(a,b)), gmul(a,b)));
    1596             :   }
    1597      494056 :   split_0(p,bit,&F,&G);
    1598      494057 :   m1 = split_complete(F,bit,roots_pol);
    1599      494056 :   m2 = split_complete(G,bit,roots_pol);
    1600      494057 :   return gerepileupto(ltop, gmul(m1,m2));
    1601             : }
    1602             : 
    1603             : static GEN
    1604     6453310 : quicktofp(GEN x)
    1605             : {
    1606     6453310 :   const long prec = DEFAULTPREC;
    1607     6453310 :   switch(typ(x))
    1608             :   {
    1609     6431854 :     case t_INT: return itor(x, prec);
    1610        9064 :     case t_REAL: return rtor(x, prec);
    1611           0 :     case t_FRAC: return fractor(x, prec);
    1612       12395 :     case t_COMPLEX: {
    1613       12395 :       GEN a = gel(x,1), b = gel(x,2);
    1614             :       /* avoid problem with 0, e.g. x = 0 + I*1e-100. We don't want |x| = 0. */
    1615       12395 :       if (isintzero(a)) return cxcompotor(b, prec);
    1616       12353 :       if (isintzero(b)) return cxcompotor(a, prec);
    1617       12353 :       a = cxcompotor(a, prec);
    1618       12353 :       b = cxcompotor(b, prec); return sqrtr(addrr(sqrr(a), sqrr(b)));
    1619             :     }
    1620           0 :     default: pari_err_TYPE("quicktofp",x);
    1621             :       return NULL;/*LCOV_EXCL_LINE*/
    1622             :   }
    1623             : 
    1624             : }
    1625             : 
    1626             : /* bound log_2 |largest root of p| (Fujiwara's bound) */
    1627             : double
    1628     2008800 : fujiwara_bound(GEN p)
    1629             : {
    1630     2008800 :   pari_sp av = avma;
    1631     2008800 :   long i, n = degpol(p);
    1632             :   GEN cc;
    1633             :   double loglc, Lmax;
    1634             : 
    1635     2008799 :   if (n <= 0) pari_err_CONSTPOL("fujiwara_bound");
    1636     2008799 :   loglc = mydbllog2r( quicktofp(gel(p,n+2)) ); /* log_2 |lc(p)| */
    1637     2008759 :   cc = gel(p, 2);
    1638     2008759 :   if (gequal0(cc))
    1639      614888 :     Lmax = -pariINFINITY-1;
    1640             :   else
    1641     1393890 :     Lmax = (mydbllog2r(quicktofp(cc)) - loglc - 1) / n;
    1642     6814833 :   for (i = 1; i < n; i++)
    1643             :   {
    1644     4806088 :     GEN y = gel(p,i+2);
    1645             :     double L;
    1646     4806088 :     if (gequal0(y)) continue;
    1647     3050711 :     L = (mydbllog2r(quicktofp(y)) - loglc) / (n-i);
    1648     3050689 :     if (L > Lmax) Lmax = L;
    1649             :   }
    1650     2008745 :   return gc_double(av, Lmax+1);
    1651             : }
    1652             : 
    1653             : /* Fujiwara's bound, real roots. Based on the following remark: if
    1654             :  *   p = x^n + sum a_i x^i and q = x^n + sum min(a_i,0)x^i
    1655             :  * then for all x >= 0, p(x) >= q(x). Thus any bound for the (positive) roots
    1656             :  * of q is a bound for the positive roots of p. */
    1657             : double
    1658     1151914 : fujiwara_bound_real(GEN p, long sign)
    1659             : {
    1660     1151914 :   pari_sp av = avma;
    1661             :   GEN x;
    1662     1151914 :   long n = degpol(p), i, signodd, signeven;
    1663     1151912 :   if (n <= 0) pari_err_CONSTPOL("fujiwara_bound");
    1664     1151912 :   x = shallowcopy(p);
    1665     1151910 :   if (gsigne(gel(x, n+2)) > 0)
    1666     1151889 :   { signeven = 1; signodd = sign; }
    1667             :   else
    1668          21 :   { signeven = -1; signodd = -sign; }
    1669     4953023 :   for (i = 0; i < n; i++)
    1670             :   {
    1671     3801113 :     if ((n - i) % 2)
    1672     2205079 :     { if (gsigne(gel(x, i+2)) == signodd ) gel(x, i+2) = gen_0; }
    1673             :     else
    1674     1596034 :     { if (gsigne(gel(x, i+2)) == signeven) gel(x, i+2) = gen_0; }
    1675             :   }
    1676     1151910 :   return gc_double(av, fujiwara_bound(x));
    1677             : }
    1678             : 
    1679             : static GEN
    1680     2160142 : mygprecrc_special(GEN x, long prec, long e)
    1681             : {
    1682             :   GEN y;
    1683     2160142 :   switch(typ(x))
    1684             :   {
    1685       37263 :     case t_REAL:
    1686       37263 :       if (!signe(x)) return real_0_bit(minss(e, expo(x)));
    1687       36115 :       return (prec > realprec(x))? rtor(x, prec): x;
    1688       13678 :     case t_COMPLEX:
    1689       13678 :       y = cgetg(3,t_COMPLEX);
    1690       13678 :       gel(y,1) = mygprecrc_special(gel(x,1),prec,e);
    1691       13678 :       gel(y,2) = mygprecrc_special(gel(x,2),prec,e);
    1692       13678 :       return y;
    1693     2109201 :     default: return x;
    1694             :   }
    1695             : }
    1696             : 
    1697             : /* like mygprec but keep at least the same precision as before */
    1698             : static GEN
    1699      530517 : mygprec_special(GEN x, long bit)
    1700             : {
    1701      530517 :   long e = gexpo(x) - bit, prec = nbits2prec(bit);
    1702      530512 :   switch(typ(x))
    1703             :   {
    1704     2663298 :     case t_POL: pari_APPLY_pol_normalized(mygprecrc_special(gel(x,i),prec,e));
    1705           0 :     default: return mygprecrc_special(x,prec,e);
    1706             :   }
    1707             : }
    1708             : 
    1709             : static GEN
    1710      394028 : fix_roots1(GEN R)
    1711             : {
    1712      394028 :   long i, l = lg(R);
    1713      394028 :   GEN v = cgetg(l, t_VEC);
    1714     1751249 :   for (i=1; i < l; i++) { GEN r = gel(R,i); gel(v,i) = gcopy(r); gunclone(r); }
    1715      394034 :   return v;
    1716             : }
    1717             : static GEN
    1718      530507 : fix_roots(GEN R, long h, long bit)
    1719             : {
    1720             :   long i, j, c, n, prec;
    1721             :   GEN v, Z, gh;
    1722             : 
    1723      530507 :   if (h == 1) return fix_roots1(R);
    1724      136479 :   prec = nbits2prec(bit); Z = grootsof1(h, prec); gh = utoipos(h);
    1725      136480 :   n = lg(R)-1; v = cgetg(h*n + 1, t_VEC);
    1726      381529 :   for (c = i = 1; i <= n; i++)
    1727             :   {
    1728      245046 :     GEN s, r = gel(R,i);
    1729      245046 :     s = (h == 2)? gsqrt(r, prec): gsqrtn(r, gh, NULL, prec);
    1730      785489 :     for (j = 1; j <= h; j++) gel(v, c++) = gmul(s, gel(Z,j));
    1731      245025 :     gunclone(r);
    1732             :   }
    1733      136483 :   return v;
    1734             : }
    1735             : 
    1736             : static GEN
    1737      529465 : all_roots(GEN p, long bit)
    1738             : {
    1739      529465 :   long bit2, i, e, h, n = degpol(p), elc = gexpo(leading_coeff(p));
    1740      529464 :   GEN q, R, m, pd = RgX_deflate_max(p, &h);
    1741      529465 :   double fb = fujiwara_bound(pd);
    1742             :   pari_sp av;
    1743             : 
    1744      529460 :   if (fb < 0) fb = 0;
    1745      529460 :   bit2 = bit + maxss(gexpo(p), 0) + (long)ceil(log2(n / h) + 2 * fb);
    1746      530507 :   for (av = avma, i = 1, e = 0;; i++, set_avma(av))
    1747             :   {
    1748      530507 :     R = vectrunc_init(n+1);
    1749      530508 :     bit2 += e + (n << i);
    1750      530508 :     q = RgX_gtofp_bit(mygprec(pd,bit2), bit2);
    1751      530509 :     q[1] = evalsigne(1)|evalvarn(0);
    1752      530509 :     m = split_complete(q, bit2, R);
    1753      530507 :     R = fix_roots(R, h, bit2);
    1754      530517 :     q = mygprec_special(pd,bit2);
    1755      530514 :     q[1] = evalsigne(1)|evalvarn(0);
    1756      530514 :     e = gexpo(RgX_sub(q, m)) - elc + (long)log2((double)n) + 1;
    1757      530511 :     if (e < 0)
    1758             :     {
    1759      530511 :       if (e < -2*bit2) e = -2*bit2; /* avoid e = -oo */
    1760      530511 :       e = bit + a_posteriori_errors(p, R, e);
    1761      530511 :       if (e < 0) return R;
    1762             :     }
    1763        1044 :     if (DEBUGLEVEL)
    1764           0 :       err_printf("all_roots: restarting, i = %ld, e = %ld\n", i,e);
    1765             :   }
    1766             : }
    1767             : 
    1768             : INLINE int
    1769      931280 : isexactscalar(GEN x) { long tx = typ(x); return is_rational_t(tx); }
    1770             : 
    1771             : static int
    1772      239634 : isexactpol(GEN p)
    1773             : {
    1774      239634 :   long i,n = degpol(p);
    1775     1162320 :   for (i=0; i<=n; i++)
    1776      931280 :     if (!isexactscalar(gel(p,i+2))) return 0;
    1777      231040 :   return 1;
    1778             : }
    1779             : 
    1780             : /* p(0) != 0 [for efficiency] */
    1781             : static GEN
    1782      231040 : solve_exact_pol(GEN p, long bit)
    1783             : {
    1784      231040 :   long i, j, k, m, n = degpol(p), iroots = 0;
    1785      231040 :   GEN ex, factors, v = zerovec(n);
    1786             : 
    1787      231040 :   factors = ZX_squff(Q_primpart(p), &ex);
    1788      462080 :   for (i=1; i<lg(factors); i++)
    1789             :   {
    1790      231040 :     GEN roots_fact = all_roots(gel(factors,i), bit);
    1791      231040 :     n = degpol(gel(factors,i));
    1792      231040 :     m = ex[i];
    1793      922042 :     for (j=1; j<=n; j++)
    1794     1382004 :       for (k=1; k<=m; k++) v[++iroots] = roots_fact[j];
    1795             :   }
    1796      231040 :   return v;
    1797             : }
    1798             : 
    1799             : /* return the roots of p with absolute error bit */
    1800             : static GEN
    1801      239634 : roots_com(GEN q, long bit)
    1802             : {
    1803             :   GEN L, p;
    1804      239634 :   long v = RgX_valrem_inexact(q, &p);
    1805      239634 :   int ex = isexactpol(p);
    1806      239634 :   if (!ex) p = RgX_normalize1(p);
    1807      239634 :   if (lg(p) == 3)
    1808           0 :     L = cgetg(1,t_VEC); /* constant polynomial */
    1809             :   else
    1810      239634 :     L = ex? solve_exact_pol(p,bit): all_roots(p,bit);
    1811      239634 :   if (v)
    1812             :   {
    1813        3935 :     GEN M, z, t = gel(q,2);
    1814             :     long i, x, y, l, n;
    1815             : 
    1816        3935 :     if (isrationalzero(t)) x = -bit;
    1817             :     else
    1818             :     {
    1819           7 :       n = gexpo(t);
    1820           7 :       x = n / v; l = degpol(q);
    1821          35 :       for (i = v; i <= l; i++)
    1822             :       {
    1823          28 :         t  = gel(q,i+2);
    1824          28 :         if (isrationalzero(t)) continue;
    1825          28 :         y = (n - gexpo(t)) / i;
    1826          28 :         if (y < x) x = y;
    1827             :       }
    1828             :     }
    1829        3935 :     z = real_0_bit(x); l = v + lg(L);
    1830        3935 :     M = cgetg(l, t_VEC); L -= v;
    1831        7933 :     for (i = 1; i <= v; i++) gel(M,i) = z;
    1832       11826 :     for (     ; i <  l; i++) gel(M,i) = gel(L,i);
    1833        3935 :     L = M;
    1834             :   }
    1835      239634 :   return L;
    1836             : }
    1837             : 
    1838             : static GEN
    1839     1200150 : tocomplex(GEN x, long l, long bit)
    1840             : {
    1841             :   GEN y;
    1842     1200150 :   if (typ(x) == t_COMPLEX)
    1843             :   {
    1844     1180743 :     if (signe(gel(x,1))) return mygprecrc(x, l, -bit);
    1845      137756 :     x = gel(x,2);
    1846      137756 :     y = cgetg(3,t_COMPLEX);
    1847      137764 :     gel(y,1) = real_0_bit(-bit);
    1848      137762 :     gel(y,2) = mygprecrc(x, l, -bit);
    1849             :   }
    1850             :   else
    1851             :   {
    1852       19407 :     y = cgetg(3,t_COMPLEX);
    1853       19407 :     gel(y,1) = mygprecrc(x, l, -bit);
    1854       19407 :     gel(y,2) = real_0_bit(-bit);
    1855             :   }
    1856      157167 :   return y;
    1857             : }
    1858             : 
    1859             : /* x,y are t_COMPLEX of t_REALs or t_REAL, compare wrt |Im x| - |Im y|,
    1860             :  * then Re x - Re y, up to 2^-e absolute error */
    1861             : static int
    1862     2227047 : cmp_complex_appr(void *E, GEN x, GEN y)
    1863             : {
    1864     2227047 :   long e = (long)E;
    1865             :   GEN z, xi, yi, xr, yr;
    1866             :   long sz, sxi, syi;
    1867     2227047 :   if (typ(x) == t_COMPLEX) { xr = gel(x,1); xi = gel(x,2); sxi = signe(xi); }
    1868      836744 :   else { xr = x; xi = NULL; sxi = 0; }
    1869     2227047 :   if (typ(y) == t_COMPLEX) { yr = gel(y,1); yi = gel(y,2); syi = signe(yi); }
    1870      558400 :   else { yr = y; yi = NULL; syi = 0; }
    1871             :   /* Compare absolute values of imaginary parts */
    1872     2227047 :   if (!sxi)
    1873             :   {
    1874      856117 :     if (syi && expo(yi) >= e) return -1;
    1875             :     /* |Im x| ~ |Im y| ~ 0 */
    1876             :   }
    1877     1370930 :   else if (!syi)
    1878             :   {
    1879       50191 :     if (sxi && expo(xi) >= e) return 1;
    1880             :     /* |Im x| ~ |Im y| ~ 0 */
    1881             :   }
    1882             :   else
    1883             :   {
    1884     1320739 :     z = addrr_sign(xi, 1, yi, -1); sz = signe(z);
    1885     1320705 :     if (sz && expo(z) >= e) return (int)sz;
    1886             :   }
    1887             :   /* |Im x| ~ |Im y|, sort according to real parts */
    1888     1329086 :   z = subrr(xr, yr); sz = signe(z);
    1889     1329082 :   if (sz && expo(z) >= e) return (int)sz;
    1890             :   /* Re x ~ Re y. Place negative imaginary part before positive */
    1891      585309 :   return (int) (sxi - syi);
    1892             : }
    1893             : 
    1894             : static GEN
    1895      529509 : clean_roots(GEN L, long l, long bit, long clean)
    1896             : {
    1897      529509 :   long i, n = lg(L), ex = 5 - bit;
    1898      529509 :   GEN res = cgetg(n,t_COL);
    1899     2427140 :   for (i=1; i<n; i++)
    1900             :   {
    1901     1897631 :     GEN c = gel(L,i);
    1902     1897631 :     if (clean && isrealappr(c,ex))
    1903             :     {
    1904      697490 :       if (typ(c) == t_COMPLEX) c = gel(c,1);
    1905      697490 :       c = mygprecrc(c, l, -bit);
    1906             :     }
    1907             :     else
    1908     1200147 :       c = tocomplex(c, l, bit);
    1909     1897631 :     gel(res,i) = c;
    1910             :   }
    1911      529509 :   gen_sort_inplace(res, (void*)ex, &cmp_complex_appr, NULL);
    1912      529504 :   return res;
    1913             : }
    1914             : 
    1915             : /* the vector of roots of p, with absolute error 2^(- prec2nbits(l)) */
    1916             : static GEN
    1917      239662 : roots_aux(GEN p, long l, long clean)
    1918             : {
    1919      239662 :   pari_sp av = avma;
    1920             :   long bit;
    1921             :   GEN L;
    1922             : 
    1923      239662 :   if (typ(p) != t_POL)
    1924             :   {
    1925          21 :     if (gequal0(p)) pari_err_ROOTS0("roots");
    1926          14 :     if (!isvalidcoeff(p)) pari_err_TYPE("roots",p);
    1927           7 :     return cgetg(1,t_COL); /* constant polynomial */
    1928             :   }
    1929      239641 :   if (!signe(p)) pari_err_ROOTS0("roots");
    1930      239641 :   checkvalidpol(p,"roots");
    1931      239634 :   if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
    1932      239634 :   if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
    1933      239634 :   bit = prec2nbits(l);
    1934      239634 :   L = roots_com(p, bit);
    1935      239634 :   return gerepilecopy(av, clean_roots(L, l, bit, clean));
    1936             : }
    1937             : GEN
    1938        8018 : roots(GEN p, long l) { return roots_aux(p,l, 0); }
    1939             : /* clean up roots. If root is real replace it by its real part */
    1940             : GEN
    1941      231644 : cleanroots(GEN p, long l) { return roots_aux(p,l, 1); }
    1942             : 
    1943             : /* private variant of conjvec. Allow non rational coefficients, shallow
    1944             :  * function. */
    1945             : GEN
    1946          84 : polmod_to_embed(GEN x, long prec)
    1947             : {
    1948          84 :   GEN v, T = gel(x,1), A = gel(x,2);
    1949             :   long i, l;
    1950          84 :   if (typ(A) != t_POL || varn(A) != varn(T))
    1951             :   {
    1952           7 :     checkvalidpol(T,"polmod_to_embed");
    1953           7 :     return const_col(degpol(T), A);
    1954             :   }
    1955          77 :   v = cleanroots(T,prec); l = lg(v);
    1956         231 :   for (i=1; i<l; i++) gel(v,i) = poleval(A,gel(v,i));
    1957          77 :   return v;
    1958             : }
    1959             : 
    1960             : GEN
    1961      289876 : QX_complex_roots(GEN p, long l)
    1962             : {
    1963      289876 :   pari_sp av = avma;
    1964             :   long bit, v;
    1965             :   GEN L;
    1966             : 
    1967      289876 :   if (!signe(p)) pari_err_ROOTS0("QX_complex_roots");
    1968      289876 :   if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
    1969      289876 :   if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
    1970      289876 :   bit = prec2nbits(l);
    1971      289876 :   v = RgX_valrem(p, &p);
    1972      289877 :   L = lg(p) > 3? all_roots(Q_primpart(p), bit): cgetg(1,t_COL);
    1973      289875 :   if (v) L = shallowconcat(const_vec(v, real_0_bit(-bit)), L);
    1974      289875 :   return gerepilecopy(av, clean_roots(L, l, bit, 1));
    1975             : }
    1976             : 
    1977             : /********************************************************************/
    1978             : /**                                                                **/
    1979             : /**                REAL ROOTS OF INTEGER POLYNOMIAL                **/
    1980             : /**                                                                **/
    1981             : /********************************************************************/
    1982             : 
    1983             : /* Count sign changes in the coefficients of (x+1)^deg(P)*P(1/(x+1)), P
    1984             :  * has no rational root. The inversion is implicit (we take coefficients
    1985             :  * backwards). */
    1986             : static long
    1987     5208444 : X2XP1(GEN P, GEN *Premapped)
    1988             : {
    1989     5208444 :   const pari_sp av = avma;
    1990     5208444 :   GEN v = shallowcopy(P);
    1991     5208531 :   long i, j, nb, s, dP = degpol(P), vlim = dP+2;
    1992             : 
    1993    31967192 :   for (j = 2; j < vlim; j++) gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
    1994     5208118 :   s = -signe(gel(v, vlim));
    1995     5208118 :   vlim--; nb = 0;
    1996    15093680 :   for (i = 1; i < dP; i++)
    1997             :   {
    1998    13080909 :     long s2 = -signe(gel(v, 2));
    1999    13080909 :     int flag = (s2 == s);
    2000    87463086 :     for (j = 2; j < vlim; j++)
    2001             :     {
    2002    74382078 :       gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
    2003    74382177 :       if (flag) flag = (s2 != signe(gel(v, j+1)));
    2004             :     }
    2005    13081008 :     if (s == signe(gel(v, vlim)))
    2006             :     {
    2007     4518040 :       if (++nb >= 2) return gc_long(av,2);
    2008     3322356 :       s = -s;
    2009             :     }
    2010             :     /* if flag is set there will be no further sign changes */
    2011    11885324 :     if (flag && (!Premapped || !nb)) return gc_long(av, nb);
    2012     9885130 :     vlim--;
    2013     9885130 :     if (gc_needed(av, 3))
    2014             :     {
    2015           0 :       if (DEBUGMEM>1) pari_warn(warnmem, "X2XP1, i = %ld/%ld", i, dP-1);
    2016           0 :       if (!Premapped) setlg(v, vlim + 2);
    2017           0 :       v = gerepilecopy(av, v);
    2018             :     }
    2019             :   }
    2020     2012771 :   if (vlim >= 2 && s == signe(gel(v, vlim))) nb++;
    2021     2012771 :   if (Premapped && nb == 1) *Premapped = v; else set_avma(av);
    2022     2012475 :   return nb;
    2023             : }
    2024             : 
    2025             : static long
    2026           0 : _intervalcmp(GEN x, GEN y)
    2027             : {
    2028           0 :   if (typ(x) == t_VEC) x = gel(x, 1);
    2029           0 :   if (typ(y) == t_VEC) y = gel(y, 1);
    2030           0 :   return gcmp(x, y);
    2031             : }
    2032             : 
    2033             : static GEN
    2034    11120853 : _gen_nored(void *E, GEN x) { (void)E; return x; }
    2035             : static GEN
    2036    24576412 : _mp_add(void *E, GEN x, GEN y) { (void)E; return mpadd(x, y); }
    2037             : static GEN
    2038           0 : _mp_sub(void *E, GEN x, GEN y) { (void)E; return mpsub(x, y); }
    2039             : static GEN
    2040     4360522 : _mp_mul(void *E, GEN x, GEN y) { (void)E; return mpmul(x, y); }
    2041             : static GEN
    2042     6273900 : _mp_sqr(void *E, GEN x) { (void)E; return mpsqr(x); }
    2043             : static GEN
    2044    14372396 : _gen_one(void *E) { (void)E; return gen_1; }
    2045             : static GEN
    2046      323472 : _gen_zero(void *E) { (void)E; return gen_0; }
    2047             : 
    2048             : static struct bb_algebra mp_algebra = { _gen_nored, _mp_add, _mp_sub,
    2049             :                          _mp_mul, _mp_sqr, _gen_one, _gen_zero };
    2050             : 
    2051             : static GEN
    2052    34598447 : _mp_cmul(void *E, GEN P, long a, GEN x) {(void)E; return mpmul(gel(P,a+2), x);}
    2053             : 
    2054             : /* Split the polynom P in two parts, whose coeffs have constant sign:
    2055             :  * P(X) = X^D*Pp + Pm. Also compute the two parts of the derivative of P,
    2056             :  * Pprimem = Pm', Pprimep = X*Pp'+ D*Pp => P' = X^(D-1)*Pprimep + Pprimem;
    2057             :  * Pprimep[i] = (i+D) Pp[i]. Return D */
    2058             : static long
    2059      165671 : split_pols(GEN P, GEN *pPp, GEN *pPm, GEN *pPprimep, GEN *pPprimem)
    2060             : {
    2061      165671 :   long i, D, dP = degpol(P), s0 = signe(gel(P,2));
    2062             :   GEN Pp, Pm, Pprimep, Pprimem;
    2063      509959 :   for(i=1; i <= dP; i++)
    2064      509963 :     if (signe(gel(P, i+2)) == -s0) break;
    2065      165670 :   D = i;
    2066      165670 :   Pm = cgetg(D + 2, t_POL);
    2067      165679 :   Pprimem = cgetg(D + 1, t_POL);
    2068      165678 :   Pp = cgetg(dP-D + 3, t_POL);
    2069      165682 :   Pprimep = cgetg(dP-D + 3, t_POL);
    2070      165691 :   Pm[1] = Pp[1] = Pprimem[1] = Pprimep[1] = P[1];
    2071      675655 :   for(i=0; i < D; i++)
    2072             :   {
    2073      509976 :     GEN c = gel(P, i+2);
    2074      509976 :     gel(Pm, i+2) = c;
    2075      509976 :     if (i) gel(Pprimem, i+1) = mului(i, c);
    2076             :   }
    2077      689722 :   for(; i <= dP; i++)
    2078             :   {
    2079      524057 :     GEN c = gel(P, i+2);
    2080      524057 :     gel(Pp, i+2-D) = c;
    2081      524057 :     gel(Pprimep, i+2-D) = mului(i, c);
    2082             :   }
    2083      165665 :   *pPm = normalizepol_lg(Pm, D+2);
    2084      165672 :   *pPprimem = normalizepol_lg(Pprimem, D+1);
    2085      165681 :   *pPp = normalizepol_lg(Pp, dP-D+3);
    2086      165681 :   *pPprimep = normalizepol_lg(Pprimep, dP-D+3);
    2087      165681 :   return dP - degpol(*pPp);
    2088             : }
    2089             : 
    2090             : static GEN
    2091     5195907 : bkeval_single_power(long d, GEN V)
    2092             : {
    2093     5195907 :   long mp = lg(V) - 2;
    2094     5195907 :   if (d > mp) return gmul(gpowgs(gel(V, mp+1), d/mp), gel(V, (d%mp)+1));
    2095     5195907 :   return gel(V, d+1);
    2096             : }
    2097             : 
    2098             : static GEN
    2099     5195889 : splitpoleval(GEN Pp, GEN Pm, GEN pows, long D, long bitprec)
    2100             : {
    2101     5195889 :   GEN vp = gen_bkeval_powers(Pp, degpol(Pp), pows, NULL, &mp_algebra, _mp_cmul);
    2102     5195772 :   GEN vm = gen_bkeval_powers(Pm, degpol(Pm), pows, NULL, &mp_algebra, _mp_cmul);
    2103     5195779 :   GEN xa = bkeval_single_power(D, pows);
    2104             :   GEN r;
    2105     5195934 :   if (!signe(vp)) return vm;
    2106     5195934 :   vp = gmul(vp, xa);
    2107     5194772 :   r = gadd(vp, vm);
    2108     5191335 :   if (gexpo(vp) - (signe(r)? gexpo(r): 0) > prec2nbits(realprec(vp)) - bitprec)
    2109      339935 :     return NULL;
    2110     4852517 :   return r;
    2111             : }
    2112             : 
    2113             : /* optimized Cauchy bound for P = X^D*Pp + Pm, D > deg(Pm) */
    2114             : static GEN
    2115      165682 : splitcauchy(GEN Pp, GEN Pm, long prec)
    2116             : {
    2117      165682 :   GEN S = gel(Pp,2), A = gel(Pm,2);
    2118      165682 :   long i, lPm = lg(Pm), lPp = lg(Pp);
    2119      506900 :   for (i=3; i < lPm; i++) { GEN c = gel(Pm,i); if (abscmpii(A, c) < 0) A = c; }
    2120      524082 :   for (i=3; i < lPp; i++) S = addii(S, gel(Pp, i));
    2121      165678 :   return subsr(1, rdivii(A, S, prec)); /* 1 + |Pm|_oo / |Pp|_1 */
    2122             : }
    2123             : 
    2124             : static GEN
    2125       15325 : ZX_deg1root(GEN P, long prec)
    2126             : {
    2127       15325 :   GEN a = gel(P,3), b = gel(P,2);
    2128       15325 :   if (is_pm1(a))
    2129             :   {
    2130       15325 :     b = itor(b, prec); if (signe(a) > 0) togglesign(b);
    2131       15325 :     return b;
    2132             :   }
    2133           0 :   return rdivii(negi(b), a, prec);
    2134             : }
    2135             : 
    2136             : /* Newton for polynom P, P(0)!=0, with unique sign change => one root in ]0,oo[
    2137             :  * P' has also at most one zero there */
    2138             : static GEN
    2139      165675 : polsolve(GEN P, long bitprec)
    2140             : {
    2141             :   pari_sp av;
    2142             :   GEN Pp, Pm, Pprimep, Pprimem, Pprime, Pprime2, ra, rb, rc, Pc;
    2143      165675 :   long dP = degpol(P), prec = nbits2prec(bitprec);
    2144             :   long expoold, iter, D, rt, s0, bitaddprec, cprec, PREC;
    2145             : 
    2146      165676 :   if (dP == 1) return ZX_deg1root(P, prec);
    2147      165676 :   Pprime = ZX_deriv(P);
    2148      165673 :   Pprime2 = ZX_deriv(Pprime);
    2149      165674 :   bitaddprec = 1 + 2*expu(dP); PREC = prec + nbits2prec(bitaddprec);
    2150      165671 :   D = split_pols(P, &Pp, &Pm, &Pprimep, &Pprimem); /* P = X^D*Pp + Pm */
    2151      165681 :   s0 = signe(gel(P, 2));
    2152      165681 :   rt = maxss(D, brent_kung_optpow(maxss(degpol(Pp), degpol(Pm)), 2, 1));
    2153      165681 :   rb = splitcauchy(Pp, Pm, DEFAULTPREC);
    2154      165673 :   for (cprec = DEFAULTPREC, expoold = LONG_MAX;;)
    2155           0 :   {
    2156      165673 :     GEN pows = gen_powers(rb, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
    2157      165684 :     Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec);
    2158      165677 :     if (!Pc) { cprec += EXTRAPREC64; rb = rtor(rb, cprec); continue; }
    2159      165677 :     if (signe(Pc) != s0) break;
    2160           0 :     shiftr_inplace(rb,1);
    2161             :   }
    2162      165677 :   for (iter = 0, ra = NULL;;)
    2163     1810570 :   {
    2164             :     GEN wdth;
    2165     1976247 :     iter++;
    2166     1976247 :     if (ra)
    2167      901493 :       rc = shiftr(addrr(ra, rb), -1);
    2168             :     else
    2169     1074754 :       rc = shiftr(rb, -1);
    2170             :     for(;;)
    2171           0 :     {
    2172     1976461 :       GEN pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
    2173     1976209 :       Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec+2);
    2174     1976043 :       if (Pc) break;
    2175           0 :       cprec += EXTRAPREC64;
    2176           0 :       rc = rtor(rc, cprec);
    2177             :     }
    2178     1976043 :     if (signe(Pc) == s0)
    2179      592880 :       ra = rc;
    2180             :     else
    2181     1383163 :       rb = rc;
    2182     1976043 :     if (!ra) continue;
    2183     1066938 :     wdth = subrr(rb, ra);
    2184     1067069 :     if (!(iter % 8))
    2185             :     {
    2186      166859 :       GEN m1 = poleval(Pprime, ra), M2;
    2187      166859 :       if (signe(m1) == s0) continue;
    2188      165704 :       M2 = poleval(Pprime2, rb);
    2189      165699 :       if (abscmprr(gmul(M2, wdth), shiftr(m1, 1)) > 0) continue;
    2190      162544 :       break;
    2191             :     }
    2192      900210 :     else if (gexpo(wdth) <= -bitprec)
    2193        3169 :       break;
    2194             :   }
    2195      165713 :   rc = rb; av = avma;
    2196     1362612 :   for(;; rc = gerepileuptoleaf(av, rc))
    2197     1362815 :   {
    2198             :     long exponew;
    2199     1528528 :     GEN Ppc, dist, rcold = rc;
    2200     1528528 :     GEN pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
    2201     1528258 :     Ppc = splitpoleval(Pprimep, Pprimem, pows, D-1, bitaddprec+4);
    2202     1528017 :     if (Ppc) Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec+4);
    2203     1528131 :     if (!Ppc || !Pc)
    2204             :     {
    2205      339943 :       if (cprec >= PREC)
    2206       44179 :         cprec += EXTRAPREC64;
    2207             :       else
    2208      295764 :         cprec = minss(2*cprec, PREC);
    2209      339947 :       rc = rtor(rc, cprec); continue; /* backtrack one step */
    2210             :     }
    2211     1188188 :     dist = typ(Ppc) == t_REAL? divrr(Pc, Ppc): divri(Pc, Ppc);
    2212     1188377 :     rc = subrr(rc, dist);
    2213     1187847 :     if (cmprr(ra, rc) > 0 || cmprr(rb, rc) < 0)
    2214             :     {
    2215           0 :       if (cprec >= PREC) break;
    2216           0 :       cprec = minss(2*cprec, PREC);
    2217           0 :       rc = rtor(rcold, cprec); continue; /* backtrack one step */
    2218             :     }
    2219     1188322 :     if (expoold == LONG_MAX) { expoold = expo(dist); continue; }
    2220      969326 :     exponew = expo(dist);
    2221      969326 :     if (exponew < -bitprec - 1)
    2222             :     {
    2223      231232 :       if (cprec >= PREC) break;
    2224       65552 :       cprec = minss(2*cprec, PREC);
    2225       65554 :       rc = rtor(rc, cprec); continue;
    2226             :     }
    2227      738094 :     if (exponew > expoold - 2)
    2228             :     {
    2229       53328 :       if (cprec >= PREC) break;
    2230       53328 :       expoold = LONG_MAX;
    2231       53328 :       cprec = minss(2*cprec, PREC);
    2232       53329 :       rc = rtor(rc, cprec); continue;
    2233             :     }
    2234      684766 :     expoold = exponew;
    2235             :   }
    2236      165680 :   return rtor(rc, prec);
    2237             : }
    2238             : 
    2239             : /* Return primpart(P(x / 2)) */
    2240             : static GEN
    2241     1937569 : ZX_rescale2prim(GEN P)
    2242             : {
    2243     1937569 :   long i, l = lg(P), v, n;
    2244             :   GEN Q;
    2245     1937569 :   if (l==2) return pol_0(varn(P));
    2246     1937569 :   Q = cgetg(l,t_POL); v = vali(gel(P,l-1));
    2247     9716770 :   for (i = l-2, n = 1; v > n && i >= 2; i--, n++)
    2248     7779159 :     v = minss(v, vali(gel(P,i)) + n);
    2249     1937611 :   gel(Q,l-1) = v? shifti(gel(P,l-1), -v): gel(P,l-1);
    2250    11644345 :   for (i = l-2, n = 1-v; i >= 2; i--, n++)
    2251     9706797 :     gel(Q,i) = shifti(gel(P,i), n);
    2252     1937548 :   Q[1] = P[1]; return Q;
    2253             : }
    2254             : 
    2255             : /* assume Q0 has no rational root */
    2256             : static GEN
    2257      940287 : usp(GEN Q0, long flag, long bitprec)
    2258             : {
    2259      940287 :   const pari_sp av = avma;
    2260             :   GEN Qremapped, Q, c, Lc, Lk, sol;
    2261      940287 :   GEN *pQremapped = flag == 1? &Qremapped: NULL;
    2262      940287 :   const long prec = nbits2prec(bitprec), deg = degpol(Q0);
    2263      940288 :   long listsize = 64, nbr = 0, nb_todo, ind, indf, i, k, nb;
    2264             : 
    2265      940288 :   sol = zerocol(deg);
    2266      940316 :   Lc = zerovec(listsize);
    2267      940337 :   Lk = cgetg(listsize+1, t_VECSMALL);
    2268      940335 :   k = Lk[1] = 0;
    2269      940335 :   ind = 1; indf = 2;
    2270      940335 :   Q = Q0;
    2271      940335 :   c = gen_0;
    2272      940335 :   nb_todo = 1;
    2273     6148569 :   while (nb_todo)
    2274             :   {
    2275     5208257 :     GEN nc = gel(Lc, ind);
    2276             :     pari_sp av2;
    2277     5208257 :     if (Lk[ind] == k + 1)
    2278             :     {
    2279     1937569 :       Q = Q0 = ZX_rescale2prim(Q0);
    2280     1937562 :       c = gen_0;
    2281             :     }
    2282     5208250 :     if (!equalii(nc, c)) Q = ZX_translate(Q, subii(nc, c));
    2283     5208313 :     av2 = avma;
    2284     5208313 :     k = Lk[ind];
    2285     5208313 :     ind++;
    2286     5208313 :     c = nc;
    2287     5208313 :     nb_todo--;
    2288     5208313 :     nb = X2XP1(Q, pQremapped);
    2289             : 
    2290     5208098 :     if (nb == 1)
    2291             :     { /* exactly one root */
    2292     1602112 :       GEN s = gen_0;
    2293     1602112 :       if (flag == 0)
    2294             :       {
    2295           0 :         s = mkvec2(gmul2n(c,-k), gmul2n(addiu(c,1),-k));
    2296           0 :         s = gerepilecopy(av2, s);
    2297             :       }
    2298     1602112 :       else if (flag == 1) /* Caveat: Qremapped is the reciprocal polynomial */
    2299             :       {
    2300      165674 :         s = polsolve(*pQremapped, bitprec+1);
    2301      165685 :         s = addir(c, divrr(s, addsr(1, s)));
    2302      165687 :         shiftr_inplace(s, -k);
    2303      165687 :         if (realprec(s) != prec) s = rtor(s, prec);
    2304      165689 :         s = gerepileupto(av2, s);
    2305             :       }
    2306     1436438 :       else set_avma(av2);
    2307     1602143 :       gel(sol, ++nbr) = s;
    2308             :     }
    2309     3605986 :     else if (nb)
    2310             :     { /* unknown, add two nodes to refine */
    2311     2134089 :       if (indf + 2 > listsize)
    2312             :       {
    2313        1788 :         if (ind>1)
    2314             :         {
    2315        5297 :           for (i = ind; i < indf; i++)
    2316             :           {
    2317        3509 :             gel(Lc, i-ind+1) = gel(Lc, i);
    2318        3509 :             Lk[i-ind+1] = Lk[i];
    2319             :           }
    2320        1788 :           indf -= ind-1;
    2321        1788 :           ind = 1;
    2322             :         }
    2323        1788 :         if (indf + 2 > listsize)
    2324             :         {
    2325           0 :           listsize *= 2;
    2326           0 :           Lc = vec_lengthen(Lc, listsize);
    2327           0 :           Lk = vecsmall_lengthen(Lk, listsize);
    2328             :         }
    2329      112711 :         for (i = indf; i <= listsize; i++) gel(Lc, i) = gen_0;
    2330             :       }
    2331     2134089 :       gel(Lc, indf) = nc = shifti(c, 1);
    2332     2134084 :       gel(Lc, indf + 1) = addiu(nc, 1);
    2333     2134087 :       Lk[indf] = Lk[indf + 1] = k + 1;
    2334     2134087 :       indf += 2;
    2335     2134087 :       nb_todo += 2;
    2336             :     }
    2337     5208127 :     if (gc_needed(av, 2))
    2338             :     {
    2339           0 :       gerepileall(av, 6, &Q0, &Q, &c, &Lc, &Lk, &sol);
    2340           0 :       if (DEBUGMEM > 1) pari_warn(warnmem, "ZX_Uspensky", avma);
    2341             :     }
    2342             :   }
    2343      940312 :   setlg(sol, nbr+1);
    2344      940315 :   return gerepilecopy(av, sol);
    2345             : }
    2346             : 
    2347             : static GEN
    2348          14 : ZX_Uspensky_equal_yes(GEN a, long flag, long bit)
    2349             : {
    2350          14 :   if (flag == 2) return gen_1;
    2351           7 :   if (flag == 1 && typ(a) != t_REAL)
    2352             :   {
    2353           7 :     if (typ(a) == t_INT && !signe(a))
    2354           0 :       a = real_0_bit(bit);
    2355             :     else
    2356           7 :       a = gtofp(a, nbits2prec(bit));
    2357             :   }
    2358           7 :   return mkcol(a);
    2359             : }
    2360             : static GEN
    2361          21 : ZX_Uspensky_no(long flag)
    2362          21 : { return flag <= 1 ? cgetg(1, t_COL) : gen_0; }
    2363             : /* ZX_Uspensky(P, [a,a], flag) */
    2364             : static GEN
    2365          28 : ZX_Uspensky_equal(GEN P, GEN a, long flag, long bit)
    2366             : {
    2367          28 :   if (typ(a) != t_INFINITY && gequal0(poleval(P, a)))
    2368          14 :     return ZX_Uspensky_equal_yes(a, flag, bit);
    2369             :   else
    2370          14 :     return ZX_Uspensky_no(flag);
    2371             : }
    2372             : static int
    2373        3350 : sol_ok(GEN r, GEN a, GEN b) { return gcmp(a, r) <= 0 && gcmp(r, b) <= 0; }
    2374             : 
    2375             : /* P a ZX without real double roots; better if primitive and squarefree but
    2376             :  * caller should ensure that. If flag & 4 assume that P has no rational root
    2377             :  * (modest speedup) */
    2378             : GEN
    2379     1061714 : ZX_Uspensky(GEN P, GEN ab, long flag, long bitprec)
    2380             : {
    2381     1061714 :   pari_sp av = avma;
    2382             :   GEN a, b, res, sol;
    2383             :   double fb;
    2384             :   long l, nbz, deg;
    2385             : 
    2386     1061714 :   if (ab)
    2387             :   {
    2388      956935 :     if (typ(ab) == t_VEC)
    2389             :     {
    2390      929361 :       if (lg(ab) != 3) pari_err_DIM("ZX_Uspensky");
    2391      929361 :       a = gel(ab, 1);
    2392      929361 :       b = gel(ab, 2);
    2393             :     }
    2394             :     else
    2395             :     {
    2396       27574 :       a = ab;
    2397       27574 :       b = mkoo();
    2398             :     }
    2399             :   }
    2400             :   else
    2401             :   {
    2402      104779 :     a = mkmoo();
    2403      104782 :     b = mkoo();
    2404             :   }
    2405     1061717 :   if (flag & 4)
    2406             :   {
    2407      128805 :     if (gcmp(a, b) >= 0) { set_avma(av); return ZX_Uspensky_no(flag); }
    2408      128805 :     flag &= ~4;
    2409      128805 :     sol = cgetg(1, t_COL);
    2410             :   }
    2411             :   else
    2412             :   {
    2413      932912 :     switch (gcmp(a, b))
    2414             :     {
    2415           7 :       case 1: set_avma(av); return ZX_Uspensky_no(flag);
    2416          28 :       case 0: return gerepilecopy(av, ZX_Uspensky_equal(P, a, flag, bitprec));
    2417             :     }
    2418      932876 :     sol = nfrootsQ(P);
    2419             :   }
    2420     1061683 :   nbz = 0; l = lg(sol);
    2421     1061683 :   if (l > 1)
    2422             :   {
    2423             :     long i, j;
    2424        2699 :     P = RgX_div(P, roots_to_pol(sol, varn(P)));
    2425        2699 :     if (!RgV_is_ZV(sol)) P = Q_primpart(P);
    2426        6049 :     for (i = j = 1; i < l; i++)
    2427        3350 :       if (sol_ok(gel(sol,i), a, b)) gel(sol,j++) = gel(sol,i);
    2428        2699 :     setlg(sol, j);
    2429        2699 :     if (flag == 2) { nbz = j-1; sol = utoi(nbz); }
    2430        2552 :     else if (flag == 1) sol = RgC_gtofp(sol, nbits2prec(bitprec));
    2431             :   }
    2432     1058984 :   else if (flag == 2) sol = gen_0;
    2433     1061683 :   deg = degpol(P);
    2434     1061682 :   if (deg == 0) return gerepilecopy(av, sol);
    2435     1059742 :   if (typ(a) == t_INFINITY && typ(b) != t_INFINITY && gsigne(b))
    2436             :   {
    2437          28 :     fb = fujiwara_bound_real(P, -1);
    2438          28 :     if (fb <= -pariINFINITY) a = gen_0;
    2439          21 :     else if (fb < 0) a = gen_m1;
    2440          21 :     else a = negi(int2n((long)ceil(fb)));
    2441             :   }
    2442     1059742 :   if (typ(b) == t_INFINITY && typ(a) != t_INFINITY && gsigne(a))
    2443             :   {
    2444          21 :     fb = fujiwara_bound_real(P, 1);
    2445          21 :     if (fb <= -pariINFINITY) b = gen_0;
    2446          21 :     else if (fb < 0) b = gen_1;
    2447           7 :     else b = int2n((long)ceil(fb));
    2448             :   }
    2449     1059742 :   if (typ(a) != t_INFINITY && typ(b) != t_INFINITY)
    2450             :   {
    2451             :     GEN d, ad, bd, diff;
    2452             :     long i;
    2453             :     /* can occur if one of a,b was initially a t_INFINITY */
    2454       12582 :     if (gequal(a,b)) return gerepilecopy(av, sol);
    2455       12575 :     d = lcmii(Q_denom(a), Q_denom(b));
    2456       12575 :     if (is_pm1(d)) { d = NULL; ad = a; bd = b; }
    2457             :     else
    2458          14 :     { P = ZX_rescale(P, d); ad = gmul(a, d); bd = gmul(b, d); }
    2459       12575 :     diff = subii(bd, ad);
    2460       12575 :     P = ZX_affine(P, diff, ad);
    2461       12575 :     res = usp(P, flag, bitprec);
    2462       12575 :     if (flag <= 1)
    2463             :     {
    2464       34176 :       for (i = 1; i < lg(res); i++)
    2465             :       {
    2466       21916 :         GEN z = gmul(diff, gel(res, i));
    2467       21916 :         if (typ(z) == t_VEC)
    2468             :         {
    2469           0 :           gel(z, 1) = gadd(ad, gel(z, 1));
    2470           0 :           gel(z, 2) = gadd(ad, gel(z, 2));
    2471             :         }
    2472             :         else
    2473       21916 :           z = gadd(ad, z);
    2474       21916 :         if (d) z = gdiv(z, d);
    2475       21916 :         gel(res, i) = z;
    2476             :       }
    2477       12260 :       sol = shallowconcat(sol, res);
    2478             :     }
    2479             :     else
    2480         315 :       nbz += lg(res) - 1;
    2481             :   }
    2482     1059735 :   if (typ(b) == t_INFINITY && (fb=fujiwara_bound_real(P, 1)) > -pariINFINITY)
    2483             :   {
    2484      839330 :     long bp = maxss((long)ceil(fb), 0);
    2485      839330 :     res = usp(ZX_unscale2n(P, bp), flag, bitprec);
    2486      839338 :     if (flag <= 1)
    2487       70991 :       sol = shallowconcat(sol, gmul2n(res, bp));
    2488             :     else
    2489      768347 :       nbz += lg(res)-1;
    2490             :   }
    2491     1059740 :   if (typ(a) == t_INFINITY && (fb=fujiwara_bound_real(P,-1)) > -pariINFINITY)
    2492             :   {
    2493       88424 :     long i, bm = maxss((long)ceil(fb), 0);
    2494       88424 :     res = usp(ZX_unscale2n(ZX_z_unscale(P, -1), bm), flag, bitprec);
    2495       88423 :     if (flag <= 1)
    2496             :     {
    2497       74926 :       for (i = 1; i < lg(res); i++)
    2498             :       {
    2499       46943 :         GEN z = gneg(gmul2n(gel(res, i), bm));
    2500       46943 :         if (typ(z) == t_VEC) swap(gel(z, 1), gel(z, 2));
    2501       46943 :         gel(res, i) = z;
    2502             :       }
    2503       27983 :       sol = shallowconcat(res, sol);
    2504             :     }
    2505             :     else
    2506       60440 :       nbz += lg(res)-1;
    2507             :   }
    2508     1059738 :   if (flag >= 2) return utoi(nbz);
    2509       83334 :   if (flag)
    2510       83334 :     sol = sort(sol);
    2511             :   else
    2512           0 :     sol = gen_sort(sol, (void *)_intervalcmp, cmp_nodata);
    2513       83335 :   return gerepileupto(av, sol);
    2514             : }
    2515             : 
    2516             : /* x a scalar */
    2517             : static GEN
    2518          42 : rootsdeg0(GEN x)
    2519             : {
    2520          42 :   if (!is_real_t(typ(x))) pari_err_TYPE("realroots",x);
    2521          35 :   if (gequal0(x)) pari_err_ROOTS0("realroots");
    2522          14 :   return cgetg(1,t_COL); /* constant polynomial */
    2523             : }
    2524             : static void
    2525     1856430 : checkbound(GEN a)
    2526             : {
    2527     1856430 :   switch(typ(a))
    2528             :   {
    2529     1856430 :     case t_INT: case t_FRAC: case t_INFINITY: break;
    2530           0 :     default: pari_err_TYPE("polrealroots", a);
    2531             :   }
    2532     1856430 : }
    2533             : static GEN
    2534      929651 : check_ab(GEN ab)
    2535             : {
    2536             :   GEN a, b;
    2537      929651 :   if (!ab) return NULL;
    2538      928216 :   if (typ(ab) != t_VEC || lg(ab) != 3) pari_err_TYPE("polrootsreal",ab);
    2539      928216 :   a = gel(ab,1); checkbound(a);
    2540      928216 :   b = gel(ab,2); checkbound(b);
    2541      928216 :   if (typ(a) == t_INFINITY && inf_get_sign(a) < 0 &&
    2542         448 :       typ(b) == t_INFINITY && inf_get_sign(b) > 0) ab = NULL;
    2543      928216 :   return ab;
    2544             : }
    2545             : /* e^(1/h) assuming the h-th root is real, beware that sqrtnr assumes e >= 0 */
    2546             : static GEN
    2547       22762 : _sqrtnr(GEN e, long h)
    2548             : {
    2549             :   long s;
    2550             :   GEN r;
    2551       22762 :   if (h == 2) return sqrtr(e);
    2552          14 :   s = signe(e); setsigne(e, 1); /* e < 0 is possible, implies h is odd */
    2553          14 :   r = sqrtnr(e, h); if (s < 0) setsigne(r, -1);
    2554          14 :   return r;
    2555             : }
    2556             : GEN
    2557       50613 : realroots(GEN P, GEN ab, long prec)
    2558             : {
    2559       50613 :   pari_sp av = avma;
    2560       50613 :   GEN sol = NULL, fa, ex;
    2561             :   long i, j, v, l;
    2562             : 
    2563       50613 :   ab = check_ab(ab);
    2564       50613 :   if (typ(P) != t_POL) return rootsdeg0(P);
    2565       50592 :   if (!RgX_is_ZX(P)) P = RgX_rescale_to_int(P);
    2566       50592 :   switch(degpol(P))
    2567             :   {
    2568          14 :     case -1: return rootsdeg0(gen_0);
    2569           7 :     case 0: return rootsdeg0(gel(P,2));
    2570             :   }
    2571       50571 :   v = ZX_valrem(Q_primpart(P), &P);
    2572       50571 :   fa = ZX_squff(P, &ex); l = lg(fa); sol = cgetg(l + 1, t_VEC);
    2573      102683 :   for (i = 1; i < l; i++)
    2574             :   {
    2575       52112 :     GEN Pi = gel(fa, i), soli, soli2;
    2576             :     long n, h;
    2577       52112 :     if (ab) h = 1; else Pi = ZX_deflate_max(Pi, &h);
    2578       52112 :     soli = ZX_Uspensky(Pi, odd(h)? ab: gen_0, 1, prec2nbits(prec));
    2579       52112 :     n = lg(soli); soli2 = odd(h)? NULL: cgetg(n, t_COL);
    2580      119073 :     for (j = 1; j < n; j++)
    2581             :     {
    2582       66961 :       GEN r = gel(soli, j); /* != 0 */
    2583       66961 :       if (typ(r) != t_REAL) gel(soli, j) = r = gtofp(r, prec);
    2584       66961 :       if (h > 1)
    2585             :       {
    2586          77 :         gel(soli, j) = r = _sqrtnr(r, h);
    2587          77 :         if (soli2) gel(soli2, j) = negr(r);
    2588             :       }
    2589             :     }
    2590       52112 :     if (soli2) soli = shallowconcat(soli, soli2);
    2591       52112 :     if (ex[i] > 1) soli = shallowconcat1( const_vec(ex[i], soli) );
    2592       52112 :     gel(sol, i) = soli;
    2593             :   }
    2594       50571 :   if (v && (!ab || (gsigne(gel(ab,1)) <= 0 && gsigne(gel(ab,2)) >= 0)))
    2595          84 :     gel(sol, i++) = const_col(v, real_0(prec));
    2596       50571 :   setlg(sol, i); if (i == 1) { set_avma(av); return cgetg(1,t_COL); }
    2597       50557 :   return gerepileupto(av, sort(shallowconcat1(sol)));
    2598             : }
    2599             : GEN
    2600       48431 : ZX_realroots_irred(GEN P, long prec)
    2601             : {
    2602       48431 :   long dP = degpol(P), j, n, h;
    2603             :   GEN sol, sol2;
    2604             :   pari_sp av;
    2605       48431 :   if (dP == 1) retmkvec(ZX_deg1root(P, prec));
    2606       45137 :   av = avma; P = ZX_deflate_max(P, &h);
    2607       45138 :   if (h == dP)
    2608             :   {
    2609       12031 :     GEN r = _sqrtnr(ZX_deg1root(P, prec), h);
    2610       12031 :     return gerepilecopy(av, odd(h)? mkvec(r): mkvec2(negr(r), r));
    2611             :   }
    2612       33107 :   sol = ZX_Uspensky(P, odd(h)? NULL: gen_0, 1 | 4, prec2nbits(prec));
    2613       33107 :   n = lg(sol); sol2 = odd(h)? NULL: cgetg(n, t_COL);
    2614      132562 :   for (j = 1; j < n; j++)
    2615             :   {
    2616       99455 :     GEN r = gel(sol, j);
    2617       99455 :     if (typ(r) != t_REAL) gel(sol, j) = r = gtofp(r, prec);
    2618       99455 :     if (h > 1)
    2619             :     {
    2620       10654 :       gel(sol, j) = r = _sqrtnr(r, h);
    2621       10654 :       if (sol2) gel(sol2, j) = negr(r);
    2622             :     }
    2623             :   }
    2624       33107 :   if (sol2) sol = shallowconcat(sol, sol2);
    2625       33107 :   return gerepileupto(av, sort(sol));
    2626             : }
    2627             : 
    2628             : static long
    2629      120244 : ZX_sturm_i(GEN P, long flag)
    2630             : {
    2631             :   pari_sp av;
    2632      120244 :   long h, r, dP = degpol(P);
    2633      120244 :   if (dP == 1) return 1;
    2634      116916 :   av = avma; P = ZX_deflate_max(P, &h);
    2635      116916 :   if (h == dP)
    2636             :   { /* now deg P = 1 */
    2637       18052 :     if (odd(h))
    2638         665 :       r = 1;
    2639             :     else
    2640       17387 :       r = (signe(gel(P,2)) != signe(gel(P,3)))? 2: 0;
    2641       18052 :     return gc_long(av, r);
    2642             :   }
    2643       98864 :   if (odd(h))
    2644       76204 :     r = itou(ZX_Uspensky(P, NULL, flag, 0));
    2645             :   else
    2646       22660 :     r = 2*itou(ZX_Uspensky(P, gen_0, flag, 0));
    2647       98863 :   return gc_long(av,r);
    2648             : }
    2649             : /* P nonconstant, squarefree ZX */
    2650             : long
    2651      879038 : ZX_sturmpart(GEN P, GEN ab)
    2652             : {
    2653      879038 :   pari_sp av = avma;
    2654      879038 :   if (!check_ab(ab)) return ZX_sturm(P);
    2655      877634 :   return gc_long(av, itou(ZX_Uspensky(P, ab, 2, 0)));
    2656             : }
    2657             : /* P nonconstant, squarefree ZX */
    2658             : long
    2659        3740 : ZX_sturm(GEN P) { return ZX_sturm_i(P, 2); }
    2660             : /* P irreducible ZX */
    2661             : long
    2662      116504 : ZX_sturm_irred(GEN P) { return ZX_sturm_i(P, 2 + 4); }

Generated by: LCOV version 1.16