Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - rootpol.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.0 lcov report (development 23332-367b47754) Lines: 1512 1651 91.6 %
Date: 2018-12-10 05:41:52 Functions: 108 115 93.9 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /*******************************************************************/
      15             : /*                                                                 */
      16             : /*                ROOTS OF COMPLEX POLYNOMIALS                     */
      17             : /*  (original code contributed by Xavier Gourdon, INRIA RR 1852)   */
      18             : /*                                                                 */
      19             : /*******************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : static const double pariINFINITY = 1./0.;
      24             : 
      25             : static long
      26      101504 : isvalidcoeff(GEN x)
      27             : {
      28      101504 :   switch (typ(x))
      29             :   {
      30       86578 :     case t_INT: case t_REAL: case t_FRAC: return 1;
      31       14912 :     case t_COMPLEX: return isvalidcoeff(gel(x,1)) && isvalidcoeff(gel(x,2));
      32             :   }
      33          14 :   return 0;
      34             : }
      35             : 
      36             : static void
      37       14059 : checkvalidpol(GEN p, const char *f)
      38             : {
      39       14059 :   long i,n = lg(p);
      40       85718 :   for (i=2; i<n; i++)
      41       71666 :     if (!isvalidcoeff(gel(p,i))) pari_err_TYPE(f, gel(p,i));
      42       14052 : }
      43             : 
      44             : /********************************************************************/
      45             : /**                                                                **/
      46             : /**                   FAST ARITHMETIC over Z[i]                    **/
      47             : /**                                                                **/
      48             : /********************************************************************/
      49             : static THREAD long KARASQUARE_LIMIT, COOKSQUARE_LIMIT;
      50             : 
      51             : /* fast sum of x,y: t_INT or t_COMPLEX(t_INT) */
      52             : static GEN
      53    14649060 : addCC(GEN x, GEN y)
      54             : {
      55             :   GEN z;
      56             : 
      57    14649060 :   if (typ(x) == t_INT)
      58             :   {
      59    12489205 :     if (typ(y) == t_INT) return addii(x,y);
      60             :     /* ty == t_COMPLEX */
      61        2292 :     z = cgetg(3,t_COMPLEX);
      62        2292 :     gel(z,1) = addii(x, gel(y,1));
      63        2292 :     gel(z,2) = icopy(gel(y,2)); return z;
      64             :   }
      65             :   /* tx == t_COMPLEX */
      66     2159855 :   z = cgetg(3,t_COMPLEX);
      67     2159855 :   if (typ(y) == t_INT)
      68             :   {
      69       27046 :     gel(z,1) = addii(gel(x,1),y);
      70       27046 :     gel(z,2) = icopy(gel(x,2)); return z;
      71             :   }
      72             :   /* ty == t_COMPLEX */
      73     2132809 :   gel(z,1) = addii(gel(x,1),gel(y,1));
      74     2132809 :   gel(z,2) = addii(gel(x,2),gel(y,2)); return z;
      75             : }
      76             : /* fast product of x,y: t_INT or t_COMPLEX(t_INT) */
      77             : static GEN
      78    28697539 : mulCC(GEN x, GEN y)
      79             : {
      80             :   GEN z;
      81             : 
      82    28697539 :   if (typ(x) == t_INT)
      83             :   {
      84    23282029 :     if (typ(y) == t_INT) return mulii(x,y);
      85             :     /* ty == t_COMPLEX */
      86       18360 :     z = cgetg(3,t_COMPLEX);
      87       18360 :     gel(z,1) = mulii(x, gel(y,1));
      88       18360 :     gel(z,2) = mulii(x, gel(y,2)); return z;
      89             :   }
      90             :   /* tx == t_COMPLEX */
      91     5415510 :   z = cgetg(3,t_COMPLEX);
      92     5415510 :   if (typ(y) == t_INT)
      93             :   {
      94     1604129 :     gel(z,1) = mulii(gel(x,1),y);
      95     1604129 :     gel(z,2) = mulii(gel(x,2),y); return z;
      96             :   }
      97             :   /* ty == t_COMPLEX */
      98             :   {
      99     3811381 :     pari_sp av = avma, tetpil;
     100             :     GEN p1, p2;
     101             : 
     102     3811381 :     p1 = mulii(gel(x,1),gel(y,1));
     103     3811381 :     p2 = mulii(gel(x,2),gel(y,2));
     104     7622762 :     y = mulii(addii(gel(x,1),gel(x,2)),
     105     7622762 :               addii(gel(y,1),gel(y,2)));
     106     3811381 :     x = addii(p1,p2); tetpil = avma;
     107     3811381 :     gel(z,1) = subii(p1,p2);
     108     3811381 :     gel(z,2) = subii(y,x); gerepilecoeffssp(av,tetpil,z+1,2);
     109     3811381 :     return z;
     110             :   }
     111             : }
     112             : /* fast squaring x: t_INT or t_COMPLEX(t_INT) */
     113             : static GEN
     114    25766204 : sqrCC(GEN x)
     115             : {
     116             :   GEN z;
     117             : 
     118    25766204 :   if (typ(x) == t_INT) return sqri(x);
     119             :   /* tx == t_COMPLEX */
     120     5577223 :   z = cgetg(3,t_COMPLEX);
     121             :   {
     122     5577223 :     pari_sp av = avma, tetpil;
     123             :     GEN y, p1, p2;
     124             : 
     125     5577223 :     p1 = sqri(gel(x,1));
     126     5577223 :     p2 = sqri(gel(x,2));
     127     5577223 :     y = sqri(addii(gel(x,1),gel(x,2)));
     128     5577223 :     x = addii(p1,p2); tetpil = avma;
     129     5577223 :     gel(z,1) = subii(p1,p2);
     130     5577223 :     gel(z,2) = subii(y,x); gerepilecoeffssp(av,tetpil,z+1,2);
     131     5577223 :     return z;
     132             :   }
     133             : }
     134             : 
     135             : static void
     136     4543027 : set_karasquare_limit(long bit)
     137             : {
     138     4543027 :   if (bit<600)       { KARASQUARE_LIMIT=8; COOKSQUARE_LIMIT=400; }
     139        1204 :   else if (bit<2000) { KARASQUARE_LIMIT=4; COOKSQUARE_LIMIT=200; }
     140           0 :   else if (bit<3000) { KARASQUARE_LIMIT=4; COOKSQUARE_LIMIT=125; }
     141           0 :   else if (bit<5000) { KARASQUARE_LIMIT=2; COOKSQUARE_LIMIT= 75; }
     142           0 :   else               { KARASQUARE_LIMIT=1; COOKSQUARE_LIMIT= 50; }
     143     4543027 : }
     144             : 
     145             : /* assume lP > 0, lP = lgpol(P) */
     146             : static GEN
     147     9380040 : CX_square_spec(GEN P, long lP)
     148             : {
     149             :   GEN s, t;
     150     9380040 :   long i, j, l, nn, n = lP - 1;
     151             :   pari_sp av;
     152             : 
     153     9380040 :   nn = n<<1; s = cgetg(nn+3,t_POL); s[1] = evalsigne(1)|evalvarn(0);
     154     9380040 :   gel(s,2) = sqrCC(gel(P,0)); /* i = 0 */
     155    23428519 :   for (i=1; i<=n; i++)
     156             :   {
     157    14048479 :     av = avma; l = (i+1)>>1;
     158    14048479 :     t = mulCC(gel(P,0), gel(P,i)); /* j = 0 */
     159    14048479 :     for (j=1; j<l; j++) t = addCC(t, mulCC(gel(P,j), gel(P,i-j)));
     160    14048479 :     t = gmul2n(t,1);
     161    14048479 :     if ((i & 1) == 0) t = addCC(t, sqrCC(gel(P,i>>1)));
     162    14048479 :     gel(s,i+2) = gerepileupto(av, t);
     163             :   }
     164     9380040 :   gel(s,nn+2) = sqrCC(gel(P,n)); /* i = nn */
     165    16386164 :   for (   ; i<nn; i++)
     166             :   {
     167     7006124 :     av = avma; l = (i+1)>>1;
     168     7006124 :     t = mulCC(gel(P,i-n),gel(P,n)); /* j = i-n */
     169     7006124 :     for (j=i-n+1; j<l; j++) t = addCC(t, mulCC(gel(P,j),gel(P,i-j)));
     170     7006124 :     t = gmul2n(t,1);
     171     7006124 :     if ((i & 1) == 0) t = addCC(t, sqrCC(gel(P,i>>1)));
     172     7006124 :     gel(s,i+2) = gerepileupto(av, t);
     173             :   }
     174     9380040 :   return normalizepol_lg(s, nn+3);
     175             : }
     176             : /* nx = lgpol(x) */
     177             : static GEN
     178           0 : RgX_s_mulspec(GEN x, long nx, long s)
     179             : {
     180             :   GEN z, t;
     181             :   long i;
     182           0 :   if (!s || !nx) return pol_0(0);
     183           0 :   z = cgetg(nx+2, t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z + 2;
     184           0 :   for (i=0; i < nx; i++) gel(t,i) = gmulgs(gel(x,i), s);
     185           0 :   return z;
     186             : }
     187             : /* nx = lgpol(x), return x << s. Inefficient if s = 0... */
     188             : static GEN
     189           0 : RgX_shiftspec(GEN x, long nx, long s)
     190             : {
     191             :   GEN z, t;
     192             :   long i;
     193           0 :   if (!nx) return pol_0(0);
     194           0 :   z = cgetg(nx+2, t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z + 2;
     195           0 :   for (i=0; i < nx; i++) gel(t,i) = gmul2n(gel(x,i), s);
     196           0 :   return z;
     197             : }
     198             : 
     199             : /* spec function. nP = lgpol(P) */
     200             : static GEN
     201     9528482 : karasquare(GEN P, long nP)
     202             : {
     203             :   GEN Q, s0, s1, s2, a, t;
     204     9528482 :   long n0, n1, i, l, N, N0, N1, n = nP - 1; /* degree(P) */
     205             :   pari_sp av;
     206             : 
     207     9528482 :   if (n <= KARASQUARE_LIMIT) return nP? CX_square_spec(P, nP): pol_0(0);
     208      147476 :   av = avma;
     209      147476 :   n0 = (n>>1) + 1; n1 = nP - n0;
     210      147476 :   s0 = karasquare(P, n0); Q = P + n0;
     211      147476 :   s2 = karasquare(Q, n1);
     212      147476 :   s1 = RgX_addspec_shallow(P, Q, n0, n1);
     213      147476 :   s1 = RgX_sub(karasquare(s1+2, lgpol(s1)), RgX_add(s0,s2));
     214      147476 :   N = (n<<1) + 1;
     215      147476 :   a = cgetg(N + 2, t_POL); a[1] = evalsigne(1)|evalvarn(0);
     216      147476 :   t = a+2; l = lgpol(s0); s0 += 2; N0 = n0<<1;
     217      147476 :   for (i=0; i < l;  i++) gel(t,i) = gel(s0,i);
     218      147476 :   for (   ; i < N0; i++) gel(t,i) = gen_0;
     219      147476 :   t = a+2 + N0; l = lgpol(s2); s2 += 2; N1 = N - N0;
     220      147476 :   for (i=0; i < l;  i++) gel(t,i) = gel(s2,i);
     221      147476 :   for (   ; i < N1; i++) gel(t,i) = gen_0;
     222      147476 :   t = a+2 + n0; l = lgpol(s1); s1 += 2;
     223      147476 :   for (i=0; i < l; i++)  gel(t,i) = gadd(gel(t,i), gel(s1,i));
     224      147476 :   return gerepilecopy(av, normalizepol_lg(a, N+2));
     225             : }
     226             : /* spec function. nP = lgpol(P) */
     227             : static GEN
     228     9086054 : cook_square(GEN P, long nP)
     229             : {
     230             :   GEN Q, p0, p1, p2, p3, q, r, t, vp, vm;
     231     9086054 :   long n0, n3, i, j, n = nP - 1;
     232             :   pari_sp av;
     233             : 
     234     9086054 :   if (n <= COOKSQUARE_LIMIT) return  nP? karasquare(P, nP): pol_0(0);
     235           0 :   av = avma;
     236             : 
     237           0 :   n0 = (n+1) >> 2; n3 = n+1 - 3*n0;
     238           0 :   p0 = P;
     239           0 :   p1 = p0+n0;
     240           0 :   p2 = p1+n0;
     241           0 :   p3 = p2+n0; /* lgpol(p0,p1,p2) = n0, lgpol(p3) = n3 */
     242             : 
     243           0 :   q = cgetg(8,t_VEC) + 4;
     244           0 :   Q = cook_square(p0, n0);
     245           0 :   r = RgX_addspec_shallow(p0,p2, n0,n0);
     246           0 :   t = RgX_addspec_shallow(p1,p3, n0,n3);
     247           0 :   gel(q,-1) = RgX_sub(r,t);
     248           0 :   gel(q,1)  = RgX_add(r,t);
     249           0 :   r = RgX_addspec_shallow(p0,RgX_shiftspec(p2,n0, 2)+2, n0,n0);
     250           0 :   t = gmul2n(RgX_addspec_shallow(p1,RgX_shiftspec(p3,n3, 2)+2, n0,n3), 1);
     251           0 :   gel(q,-2) = RgX_sub(r,t);
     252           0 :   gel(q,2)  = RgX_add(r,t);
     253           0 :   r = RgX_addspec_shallow(p0,RgX_s_mulspec(p2,n0, 9)+2, n0,n0);
     254           0 :   t = gmulsg(3, RgX_addspec_shallow(p1,RgX_s_mulspec(p3,n3, 9)+2, n0,n3));
     255           0 :   gel(q,-3) = RgX_sub(r,t);
     256           0 :   gel(q,3)  = RgX_add(r,t);
     257             : 
     258           0 :   r = new_chunk(7);
     259           0 :   vp = cgetg(4,t_VEC);
     260           0 :   vm = cgetg(4,t_VEC);
     261           0 :   for (i=1; i<=3; i++)
     262             :   {
     263           0 :     GEN a = gel(q,i), b = gel(q,-i);
     264           0 :     a = cook_square(a+2, lgpol(a));
     265           0 :     b = cook_square(b+2, lgpol(b));
     266           0 :     gel(vp,i) = RgX_add(b, a);
     267           0 :     gel(vm,i) = RgX_sub(b, a);
     268             :   }
     269           0 :   gel(r,0) = Q;
     270           0 :   gel(r,1) = gdivgs(gsub(gsub(gmulgs(gel(vm,2),9),gel(vm,3)),
     271           0 :                      gmulgs(gel(vm,1),45)),
     272             :                 60);
     273           0 :   gel(r,2) = gdivgs(gadd(gadd(gmulgs(gel(vp,1),270),gmulgs(Q,-490)),
     274           0 :                      gadd(gmulgs(gel(vp,2),-27),gmulgs(gel(vp,3),2))),
     275             :                 360);
     276           0 :   gel(r,3) = gdivgs(gadd(gadd(gmulgs(gel(vm,1),13),gmulgs(gel(vm,2),-8)),
     277           0 :                     gel(vm,3)),
     278             :                 48);
     279           0 :   gel(r,4) = gdivgs(gadd(gadd(gmulgs(Q,56),gmulgs(gel(vp,1),-39)),
     280           0 :                      gsub(gmulgs(gel(vp,2),12),gel(vp,3))),
     281             :                 144);
     282           0 :   gel(r,5) = gdivgs(gsub(gadd(gmulgs(gel(vm,1),-5),gmulgs(gel(vm,2),4)),
     283           0 :                      gel(vm,3)),
     284             :                 240);
     285           0 :   gel(r,6) = gdivgs(gadd(gadd(gmulgs(Q,-20),gmulgs(gel(vp,1),15)),
     286           0 :                      gadd(gmulgs(gel(vp,2),-6),gel(vp,3))),
     287             :                 720);
     288           0 :   q = cgetg(2*n+3,t_POL); q[1] = evalsigne(1)|evalvarn(0);
     289           0 :   t = q+2;
     290           0 :   for (i=0; i<=2*n; i++) gel(t,i) = gen_0;
     291           0 :   for (i=0; i<=6; i++,t += n0)
     292             :   {
     293           0 :     GEN h = gel(r,i);
     294           0 :     long d = lgpol(h);
     295           0 :     h += 2;
     296           0 :     for (j=0; j<d; j++) gel(t,j) = gadd(gel(t,j), gel(h,j));
     297             :   }
     298           0 :   return gerepilecopy(av, normalizepol_lg(q, 2*n+3));
     299             : }
     300             : 
     301             : static GEN
     302     4543027 : graeffe(GEN p)
     303             : {
     304             :   GEN p0, p1, s0, s1;
     305     4543027 :   long n = degpol(p), n0, n1, i;
     306             : 
     307     4543027 :   if (!n) return gcopy(p);
     308     4543027 :   n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
     309     4543027 :   p0 = new_chunk(n0);
     310     4543027 :   p1 = new_chunk(n1);
     311    14719659 :   for (i=0; i<n1; i++)
     312             :   {
     313    10176632 :     p0[i] = p[2+(i<<1)];
     314    10176632 :     p1[i] = p[3+(i<<1)];
     315             :   }
     316     4543027 :   if (n1 != n0)
     317     2112429 :     p0[i] = p[2+(i<<1)];
     318     4543027 :   s0 = cook_square(p0, n0);
     319     4543027 :   s1 = cook_square(p1, n1);
     320     4543027 :   return RgX_sub(s0, RgX_shift_shallow(s1,1));
     321             : }
     322             : GEN
     323        7581 : ZX_graeffe(GEN p)
     324             : {
     325        7581 :   pari_sp av = avma;
     326             :   GEN p0, p1, s0, s1;
     327        7581 :   long n = degpol(p);
     328             : 
     329        7581 :   if (!n) return ZX_copy(p);
     330        7581 :   RgX_even_odd(p, &p0, &p1);
     331             :   /* p = p0(x^2) + x p1(x^2) */
     332        7581 :   s0 = ZX_sqr(p0);
     333        7581 :   s1 = ZX_sqr(p1);
     334        7581 :   return gerepileupto(av, ZX_sub(s0, RgX_shift_shallow(s1,1)));
     335             : }
     336             : GEN
     337          14 : polgraeffe(GEN p)
     338             : {
     339          14 :   pari_sp av = avma;
     340             :   GEN p0, p1, s0, s1;
     341          14 :   long n = degpol(p);
     342             : 
     343          14 :   if (typ(p) != t_POL) pari_err_TYPE("polgraeffe",p);
     344          14 :   n = degpol(p);
     345          14 :   if (!n) return gcopy(p);
     346          14 :   RgX_even_odd(p, &p0, &p1);
     347             :   /* p = p0(x^2) + x p1(x^2) */
     348          14 :   s0 = RgX_sqr(p0);
     349          14 :   s1 = RgX_sqr(p1);
     350          14 :   return gerepileupto(av, RgX_sub(s0, RgX_shift_shallow(s1,1)));
     351             : }
     352             : 
     353             : /********************************************************************/
     354             : /**                                                                **/
     355             : /**                       MODULUS OF ROOTS                         **/
     356             : /**                                                                **/
     357             : /********************************************************************/
     358             : 
     359             : /* Quick approximation to log2(|x|); first define y s.t. |y-x| < 2^-32 then
     360             :  * return y rounded to 2 ulp. In particular, if result < 2^21, absolute error
     361             :  * is bounded by 2^-31. If result > 2^21, it is correct to 2 ulp */
     362             : static double
     363    24260053 : mydbllog2i(GEN x)
     364             : {
     365             : #ifdef LONG_IS_64BIT
     366    20840340 :   const double W = 1/(4294967296. * 4294967296.); /* 2^-64 */
     367             : #else
     368     3419713 :   const double W = 1/4294967296.; /*2^-32*/
     369             : #endif
     370             :   GEN m;
     371    24260053 :   long lx = lgefint(x);
     372             :   double l;
     373    24260053 :   if (lx == 2) return -pariINFINITY;
     374    24092008 :   m = int_MSW(x);
     375    24092008 :   l = (double)(ulong)*m;
     376    24092008 :   if (lx == 3) return log2(l);
     377    10614060 :   l += ((double)(ulong)*int_precW(m)) * W;
     378             :   /* at least m = min(53,BIL) bits are correct in the mantissa, thus log2
     379             :    * is correct with error < log(1 + 2^-m) ~ 2^-m. Adding the correct
     380             :    * exponent BIL(lx-3) causes 1ulp further round-off error */
     381    10614060 :   return log2(l) + (double)(BITS_IN_LONG*(lx-3));
     382             : }
     383             : 
     384             : /* return log(|x|) or -pariINFINITY */
     385             : static double
     386     1282737 : mydbllogr(GEN x) {
     387     1282737 :   if (!signe(x)) return -pariINFINITY;
     388     1282737 :   return M_LN2*dbllog2r(x);
     389             : }
     390             : 
     391             : /* return log2(|x|) or -pariINFINITY */
     392             : static double
     393    14714273 : mydbllog2r(GEN x) {
     394    14714273 :   if (!signe(x)) return -pariINFINITY;
     395    14628596 :   return dbllog2r(x);
     396             : }
     397             : double
     398    42665807 : dbllog2(GEN z)
     399             : {
     400             :   double x, y;
     401    42665807 :   switch(typ(z))
     402             :   {
     403    24257267 :     case t_INT: return mydbllog2i(z);
     404        1393 :     case t_FRAC: return mydbllog2i(gel(z,1))-mydbllog2i(gel(z,2));
     405    13980689 :     case t_REAL: return mydbllog2r(z);
     406             :     default: /*t_COMPLEX*/
     407     4426458 :       x = dbllog2(gel(z,1));
     408     4426458 :       y = dbllog2(gel(z,2));
     409     4426458 :       if (x == -pariINFINITY) return y;
     410     4390080 :       if (y == -pariINFINITY) return x;
     411     4262552 :       if (fabs(x-y) > 10) return maxdd(x,y);
     412     4158079 :       return x + 0.5*log2(1 + exp2(2*(y-x)));
     413             :   }
     414             : }
     415             : static GEN /* beware overflow */
     416      657090 : dblexp(double x) { return fabs(x) < 100.? dbltor(exp(x)): mpexp(dbltor(x)); }
     417             : 
     418             : /* find s such that  A_h <= 2^s <= 2 A_i  for one h and all i < n = deg(p),
     419             :  * with  A_i := (binom(n,i) lc(p) / p_i) ^ 1/(n-i), and  p = sum p_i X^i */
     420             : static long
     421     3580641 : findpower(GEN p)
     422             : {
     423     3580641 :   double x, L, mins = pariINFINITY;
     424     3580641 :   long n = degpol(p),i;
     425             : 
     426     3580641 :   L = dbllog2(gel(p,n+2)); /* log2(lc * binom(n,i)) */
     427    17245438 :   for (i=n-1; i>=0; i--)
     428             :   {
     429    13664797 :     L += log2((double)(i+1) / (double)(n-i));
     430    13664797 :     x = dbllog2(gel(p,i+2));
     431    13664797 :     if (x != -pariINFINITY)
     432             :     {
     433    13578359 :       double s = (L - x) / (double)(n-i);
     434    13578359 :       if (s < mins) mins = s;
     435             :     }
     436             :   }
     437     3580641 :   i = (long)ceil(mins);
     438     3580641 :   if (i - mins > 1 - 1e-12) i--;
     439     3580641 :   return i;
     440             : }
     441             : 
     442             : /* returns the exponent for logmodulus(), from the Newton diagram */
     443             : static long
     444      572344 : newton_polygon(GEN p, long k)
     445             : {
     446      572344 :   pari_sp av = avma;
     447             :   double *logcoef, slope;
     448      572344 :   long n = degpol(p), i, j, h, l, *vertex;
     449             : 
     450      572344 :   logcoef = (double*)stack_malloc_align((n+1)*sizeof(double), sizeof(double));
     451      572344 :   vertex = (long*)new_chunk(n+1);
     452             : 
     453             :   /* vertex[i] = 1 if i a vertex of convex hull, 0 otherwise */
     454      572344 :   for (i=0; i<=n; i++) { logcoef[i] = dbllog2(gel(p,2+i)); vertex[i] = 0; }
     455      572344 :   vertex[0] = 1; /* sentinel */
     456     2353585 :   for (i=0; i < n; i=h)
     457             :   {
     458     1781241 :     slope = logcoef[i+1]-logcoef[i];
     459     8209356 :     for (j = h = i+1; j<=n; j++)
     460             :     {
     461     6428115 :       double pij = (logcoef[j]-logcoef[i])/(double)(j-i);
     462     6428115 :       if (slope < pij) { slope = pij; h = j; }
     463             :     }
     464     1781241 :     vertex[h] = 1;
     465             :   }
     466      572344 :   h = k;   while (!vertex[h]) h++;
     467      572344 :   l = k-1; while (!vertex[l]) l--;
     468      572344 :   set_avma(av);
     469      572344 :   return (long)floor((logcoef[h]-logcoef[l])/(double)(h-l) + 0.5);
     470             : }
     471             : 
     472             : /* change z into z*2^e, where z is real or complex of real */
     473             : static void
     474     3890916 : myshiftrc(GEN z, long e)
     475             : {
     476     3890916 :   if (typ(z)==t_COMPLEX)
     477             :   {
     478     1006563 :     if (signe(gel(z,1))) shiftr_inplace(gel(z,1), e);
     479     1006563 :     if (signe(gel(z,2))) shiftr_inplace(gel(z,2), e);
     480             :   }
     481             :   else
     482     2884353 :     if (signe(z)) shiftr_inplace(z, e);
     483     3890916 : }
     484             : 
     485             : /* return z*2^e, where z is integer or complex of integer (destroy z) */
     486             : static GEN
     487    14716441 : myshiftic(GEN z, long e)
     488             : {
     489    14716441 :   if (typ(z)==t_COMPLEX)
     490             :   {
     491     2946607 :     gel(z,1) = signe(gel(z,1))? mpshift(gel(z,1),e): gen_0;
     492     2946607 :     gel(z,2) = mpshift(gel(z,2),e);
     493     2946607 :     return z;
     494             :   }
     495    11769834 :   return signe(z)? mpshift(z,e): gen_0;
     496             : }
     497             : 
     498             : static GEN
     499      741892 : RgX_gtofp_bit(GEN q, long bit)
     500             : {
     501      741892 :   if (bit < 0) bit = 0;
     502      741892 :   return RgX_gtofp(q, nbits2prec(bit));
     503             : }
     504             : 
     505             : static GEN
     506    29224416 : mygprecrc(GEN x, long prec, long e)
     507             : {
     508             :   GEN y;
     509    29224416 :   switch(typ(x))
     510             :   {
     511    21659677 :     case t_REAL: return signe(x)? rtor(x, prec): real_0_bit(e);
     512             :     case t_COMPLEX:
     513     5873247 :       y = cgetg(3,t_COMPLEX);
     514     5873247 :       gel(y,1) = mygprecrc(gel(x,1),prec,e);
     515     5873247 :       gel(y,2) = mygprecrc(gel(x,2),prec,e);
     516     5873247 :       return y;
     517     1691492 :     default: return gcopy(x);
     518             :   }
     519             : }
     520             : 
     521             : /* gprec behaves badly with the zero for polynomials.
     522             : The second parameter in mygprec is the precision in base 2 */
     523             : static GEN
     524     7213678 : mygprec(GEN x, long bit)
     525             : {
     526             :   long lx, i, e, prec;
     527             :   GEN y;
     528             : 
     529     7213678 :   if (bit < 0) bit = 0; /* should rarely happen */
     530     7213678 :   e = gexpo(x) - bit;
     531     7213678 :   prec = nbits2prec(bit);
     532     7213678 :   switch(typ(x))
     533             :   {
     534             :     case t_POL:
     535     4532285 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     536     4532285 :       for (i=2; i<lx; i++) gel(y,i) = mygprecrc(gel(x,i),prec,e);
     537     4532285 :       break;
     538             : 
     539     2681393 :     default: y = mygprecrc(x,prec,e);
     540             :   }
     541     7213678 :   return y;
     542             : }
     543             : 
     544             : /* normalize a polynomial p, that is change it with coefficients in Z[i],
     545             : after making product by 2^shift */
     546             : static GEN
     547     1746482 : pol_to_gaussint(GEN p, long shift)
     548             : {
     549     1746482 :   long i, l = lg(p);
     550     1746482 :   GEN q = cgetg(l, t_POL); q[1] = p[1];
     551     1746482 :   for (i=2; i<l; i++) gel(q,i) = gtrunc2n(gel(p,i), shift);
     552     1746482 :   return q;
     553             : }
     554             : 
     555             : /* returns a polynomial q in Z[i][x] keeping bit bits of p */
     556             : static GEN
     557     1362282 : eval_rel_pol(GEN p, long bit)
     558             : {
     559             :   long i;
     560     9281273 :   for (i = 2; i < lg(p); i++)
     561     7918991 :     if (gequal0(gel(p,i))) gel(p,i) = gen_0; /* bad behavior of gexpo */
     562     1362282 :   return pol_to_gaussint(p, bit-gexpo(p)+1);
     563             : }
     564             : 
     565             : /* returns p(R*x)/R^n (in R or R[i]), R = exp(lrho), bit bits of precision */
     566             : static GEN
     567      154713 : homothetie(GEN p, double lrho, long bit)
     568             : {
     569             :   GEN q, r, t, iR;
     570      154713 :   long n = degpol(p), i;
     571             : 
     572      154713 :   iR = mygprec(dblexp(-lrho),bit);
     573      154713 :   q = mygprec(p, bit);
     574      154713 :   r = cgetg(n+3,t_POL); r[1] = p[1];
     575      154713 :   t = iR; r[n+2] = q[n+2];
     576      905830 :   for (i=n-1; i>0; i--)
     577             :   {
     578      751117 :     gel(r,i+2) = gmul(t, gel(q,i+2));
     579      751117 :     t = mulrr(t, iR);
     580             :   }
     581      154713 :   gel(r,2) = gmul(t, gel(q,2)); return r;
     582             : }
     583             : 
     584             : /* change q in 2^(n*e) p(x*2^(-e)), n=deg(q)  [ ~as above with R = 2^-e ]*/
     585             : static void
     586      956544 : homothetie2n(GEN p, long e)
     587             : {
     588      956544 :   if (e)
     589             :   {
     590      728967 :     long i,n = lg(p)-1;
     591      728967 :     for (i=2; i<=n; i++) myshiftrc(gel(p,i), (n-i)*e);
     592             :   }
     593      956544 : }
     594             : 
     595             : /* return 2^f * 2^(n*e) p(x*2^(-e)), n=deg(q) */
     596             : static void
     597     3196441 : homothetie_gauss(GEN p, long e, long f)
     598             : {
     599     3196441 :   if (e || f)
     600             :   {
     601     2904877 :     long i, n = lg(p)-1;
     602     2904877 :     for (i=2; i<=n; i++) gel(p,i) = myshiftic(gel(p,i), f+(n-i)*e);
     603             :   }
     604     3196441 : }
     605             : 
     606             : /* Lower bound on the modulus of the largest root z_0
     607             :  * k is set to an upper bound for #{z roots, |z-z_0| < eps} */
     608             : static double
     609     3580641 : lower_bound(GEN p, long *k, double eps)
     610             : {
     611     3580641 :   long n = degpol(p), i, j;
     612     3580641 :   pari_sp ltop = avma;
     613             :   GEN a, s, S, ilc;
     614             :   double r, R, rho;
     615             : 
     616     3580641 :   if (n < 4) { *k = n; return 0.; }
     617     1374025 :   S = cgetg(5,t_VEC);
     618     1374025 :   a = cgetg(5,t_VEC); ilc = gdiv(real_1(DEFAULTPREC), gel(p,n+2));
     619     1374025 :   for (i=1; i<=4; i++) gel(a,i) = gmul(ilc,gel(p,n+2-i));
     620             :   /* i = 1 split out from next loop for efficiency and initialization */
     621     1374025 :   s = gel(a,1);
     622     1374025 :   gel(S,1) = gneg(s); /* Newton sum S_i */
     623     1374025 :   rho = r = gtodouble(gabs(s,3));
     624     1374025 :   R = r / n;
     625     5496100 :   for (i=2; i<=4; i++)
     626             :   {
     627     4122075 :     s = gmulsg(i,gel(a,i));
     628     4122075 :     for (j=1; j<i; j++) s = gadd(s, gmul(gel(S,j),gel(a,i-j)));
     629     4122075 :     gel(S,i) = gneg(s); /* Newton sum S_i */
     630     4122075 :     r = gtodouble(gabs(s,3));
     631     4122075 :     if (r > 0.)
     632             :     {
     633     4113707 :       r = exp(log(r/n) / (double)i);
     634     4113707 :       if (r > R) R = r;
     635             :     }
     636             :   }
     637     1374025 :   if (R > 0. && eps < 1.2)
     638     1372785 :     *k = (long)floor((rho/R + n) / (1 + exp(-eps)*cos(eps)));
     639             :   else
     640        1240 :     *k = n;
     641     1374025 :   return gc_double(ltop, R);
     642             : }
     643             : 
     644             : /* return R such that exp(R - tau) <= rho_n(P) <= exp(R + tau)
     645             :  * P(0) != 0 and P non constant */
     646             : static double
     647      384200 : logmax_modulus(GEN p, double tau)
     648             : {
     649             :   GEN r, q, aux, gunr;
     650      384200 :   pari_sp av, ltop = avma;
     651      384200 :   long i,k,n=degpol(p),nn,bit,M,e;
     652      384200 :   double rho,eps, tau2 = (tau > 3.0)? 0.5: tau/6.;
     653             : 
     654      384200 :   r = cgeti(BIGDEFAULTPREC);
     655      384200 :   av = avma;
     656             : 
     657      384200 :   eps = - 1/log(1.5*tau2); /* > 0 */
     658      384200 :   bit = (long) ((double) n*log2(1./tau2)+3*log2((double) n))+1;
     659      384200 :   gunr = real_1_bit(bit+2*n);
     660      384200 :   aux = gdiv(gunr, gel(p,2+n));
     661      384200 :   q = RgX_Rg_mul(p, aux); gel(q,2+n) = gunr;
     662      384200 :   e = findpower(q);
     663      384200 :   homothetie2n(q,e);
     664      384200 :   affsi(e, r);
     665      384200 :   q = pol_to_gaussint(q, bit);
     666      384200 :   M = (long) (log2( log(4.*n) / (2*tau2) )) + 2;
     667      384200 :   nn = n;
     668      384200 :   for (i=0,e=0;;)
     669             :   { /* nn = deg(q) */
     670     6777082 :     rho = lower_bound(q, &k, eps);
     671     3580641 :     if (rho > exp2(-(double)e)) e = (long)-floor(log2(rho));
     672     3580641 :     affii(shifti(addis(r,e), 1), r);
     673     3580641 :     if (++i == M) break;
     674             : 
     675     9589323 :     bit = (long) ((double)k * log2(1./tau2) +
     676     6392882 :                      (double)(nn-k)*log2(1./eps) + 3*log2((double)nn)) + 1;
     677     3196441 :     homothetie_gauss(q, e, bit-(long)floor(dbllog2(gel(q,2+nn))+0.5));
     678     3196441 :     nn -= RgX_valrem(q, &q);
     679     3196441 :     set_karasquare_limit(gexpo(q));
     680     3196441 :     q = gerepileupto(av, graeffe(q));
     681     3196441 :     tau2 *= 1.5; if (tau2 > 0.9) tau2 = 0.5;
     682     3196441 :     eps = -1/log(tau2); /* > 0 */
     683     3196441 :     e = findpower(q);
     684             :   }
     685      384200 :   if (!signe(r)) return gc_double(ltop,0.);
     686      360986 :   r = itor(r, DEFAULTPREC); shiftr_inplace(r, -M);
     687      360986 :   return gc_double(ltop, -rtodbl(r) * M_LN2); /* -log(2) sum e_i 2^-i */
     688             : }
     689             : 
     690             : static GEN
     691       12959 : RgX_normalize1(GEN x)
     692             : {
     693       12959 :   long i, n = lg(x)-1;
     694             :   GEN y;
     695       12987 :   for (i = n; i > 1; i--)
     696       12980 :     if (!gequal0( gel(x,i) )) break;
     697       12959 :   if (i == n) return x;
     698          28 :   pari_warn(warner,"normalizing a polynomial with 0 leading term");
     699          28 :   if (i == 1) pari_err_ROOTS0("roots");
     700          28 :   y = cgetg(i+1, t_POL); y[1] = x[1];
     701          28 :   for (; i > 1; i--) gel(y,i) = gel(x,i);
     702          28 :   return y;
     703             : }
     704             : 
     705             : static GEN
     706        5699 : polrootsbound_i(GEN P, double TAU)
     707             : {
     708        5699 :   pari_sp av = avma;
     709             :   double d;
     710        5699 :   (void)RgX_valrem_inexact(P,&P);
     711        5699 :   P = RgX_normalize1(P);
     712        5699 :   switch(degpol(P))
     713             :   {
     714           7 :     case -1: pari_err_ROOTS0("roots");
     715          70 :     case 0:  set_avma(av); return gen_0;
     716             :   }
     717        5622 :   d = logmax_modulus(P, TAU) + TAU;
     718             :   /* not dblexp: result differs on ARM emulator */
     719        5622 :   return gerepileuptoleaf(av, mpexp(dbltor(d)));
     720             : }
     721             : GEN
     722        5706 : polrootsbound(GEN P, GEN tau)
     723             : {
     724        5706 :   if (typ(P) != t_POL) pari_err_TYPE("polrootsbound",P);
     725        5699 :   checkvalidpol(P, "polrootsbound");
     726        5699 :   return polrootsbound_i(P, tau? gtodouble(tau): 0.01);
     727             : }
     728             : 
     729             : /* log of modulus of the smallest root of p, with relative error tau */
     730             : static double
     731      135914 : logmin_modulus(GEN p, double tau)
     732             : {
     733      135914 :   pari_sp av = avma;
     734      135914 :   if (gequal0(gel(p,2))) return -pariINFINITY;
     735      135914 :   return gc_double(av, - logmax_modulus(RgX_recip_shallow(p),tau));
     736             : }
     737             : 
     738             : /* return the log of the k-th modulus (ascending order) of p, rel. error tau*/
     739             : static double
     740       62739 : logmodulus(GEN p, long k, double tau)
     741             : {
     742             :   GEN q;
     743       62739 :   long i, kk = k, imax, n = degpol(p), nn, bit, e;
     744       62739 :   pari_sp av, ltop=avma;
     745       62739 :   double r, tau2 = tau/6;
     746             : 
     747       62739 :   bit = (long)(n * (2. + log2(3.*n/tau2)));
     748       62739 :   av = avma;
     749       62739 :   q = gprec_w(p, nbits2prec(bit));
     750       62739 :   q = RgX_gtofp_bit(q, bit);
     751       62739 :   e = newton_polygon(q,k);
     752       62739 :   r = (double)e;
     753       62739 :   homothetie2n(q,e);
     754       62739 :   imax = (long)(log2(3./tau) + log2(log(4.*n)))+1;
     755      572344 :   for (i=1; i<imax; i++)
     756             :   {
     757      509605 :     q = eval_rel_pol(q,bit);
     758      509605 :     kk -= RgX_valrem(q, &q);
     759      509605 :     nn = degpol(q);
     760             : 
     761      509605 :     set_karasquare_limit(bit);
     762      509605 :     q = gerepileupto(av, graeffe(q));
     763      509605 :     e = newton_polygon(q,kk);
     764      509605 :     r += e / exp2((double)i);
     765      509605 :     q = RgX_gtofp_bit(q, bit);
     766      509605 :     homothetie2n(q,e);
     767             : 
     768      509605 :     tau2 *= 1.5; if (tau2 > 1.) tau2 = 1.;
     769      509605 :     bit = 1 + (long)(nn*(2. + log2(3.*nn/tau2)));
     770             :   }
     771       62739 :   return gc_double(ltop, -r * M_LN2);
     772             : }
     773             : 
     774             : /* return the log of the k-th modulus r_k of p, rel. error tau, knowing that
     775             :  * rmin < r_k < rmax. This information helps because we may reduce precision
     776             :  * quicker */
     777             : static double
     778       62739 : logpre_modulus(GEN p, long k, double tau, double lrmin, double lrmax)
     779             : {
     780             :   GEN q;
     781       62739 :   long n = degpol(p), i, imax, imax2, bit;
     782       62739 :   pari_sp ltop = avma, av;
     783       62739 :   double lrho, aux, tau2 = tau/6.;
     784             : 
     785       62739 :   aux = (lrmax - lrmin) / 2. + 4*tau2;
     786       62739 :   imax = (long) log2(log((double)n)/ aux);
     787       62739 :   if (imax <= 0) return logmodulus(p,k,tau);
     788             : 
     789       61527 :   lrho  = (lrmin + lrmax) / 2;
     790       61527 :   av = avma;
     791       61527 :   bit = (long)(n*(2. + aux / M_LN2 - log2(tau2)));
     792       61527 :   q = homothetie(p, lrho, bit);
     793       61527 :   imax2 = (long)(log2(3./tau * log(4.*n))) + 1;
     794       61527 :   if (imax > imax2) imax = imax2;
     795             : 
     796      187521 :   for (i=0; i<imax; i++)
     797             :   {
     798      125994 :     q = eval_rel_pol(q,bit);
     799      125994 :     set_karasquare_limit(bit);
     800      125994 :     q = gerepileupto(av, graeffe(q));
     801      125994 :     aux = 2*aux + 2*tau2;
     802      125994 :     tau2 *= 1.5;
     803      125994 :     bit = (long)(n*(2. + aux / M_LN2 - log2(1-exp(-tau2))));
     804      125994 :     q = RgX_gtofp_bit(q, bit);
     805             :   }
     806       61527 :   aux = exp2((double)imax);
     807       61527 :   return gc_double(ltop, lrho + logmodulus(q,k, aux*tau/3.) / aux);
     808             : }
     809             : 
     810             : static double
     811       77490 : ind_maxlog2(GEN q)
     812             : {
     813       77490 :   long i, k = -1;
     814       77490 :   double L = - pariINFINITY;
     815      208922 :   for (i=0; i<=degpol(q); i++)
     816             :   {
     817      131432 :     double d = dbllog2(gel(q,2+i));
     818      131432 :     if (d > L) { L = d; k = i; }
     819             :   }
     820       77490 :   return k;
     821             : }
     822             : 
     823             : /* Returns k such that r_k e^(-tau) < R < r_{k+1} e^tau.
     824             :  * Assume that l <= k <= n-l */
     825             : static long
     826       93186 : dual_modulus(GEN p, double lrho, double tau, long l)
     827             : {
     828       93186 :   long i, imax, delta_k = 0, n = degpol(p), nn, v2, v, bit, ll = l;
     829       93186 :   double tau2 = tau * 7./8.;
     830       93186 :   pari_sp av = avma;
     831             :   GEN q;
     832             : 
     833       93186 :   bit = 6*n - 5*l + (long)(n*(-log2(tau2) + tau2 * 8./7.));
     834       93186 :   q = homothetie(p, lrho, bit);
     835       93186 :   imax = (long)(log(log(2.*n)/tau2)/log(7./4.)+1);
     836             : 
     837      804173 :   for (i=0; i<imax; i++)
     838             :   {
     839      726683 :     q = eval_rel_pol(q,bit); v2 = n - degpol(q);
     840      726683 :     v = RgX_valrem(q, &q);
     841      726683 :     ll -= maxss(v, v2); if (ll < 0) ll = 0;
     842             : 
     843      726683 :     nn = degpol(q); delta_k += v;
     844      726683 :     if (!nn) return delta_k;
     845             : 
     846      710987 :     set_karasquare_limit(bit);
     847      710987 :     q = gerepileupto(av, graeffe(q));
     848      710987 :     tau2 *= 7./4.;
     849      710987 :     bit = 6*nn - 5*ll + (long)(nn*(-log2(tau2) + tau2 * 8./7.));
     850             :   }
     851       77490 :   return gc_long(av, delta_k + (long)ind_maxlog2(q));
     852             : }
     853             : 
     854             : /********************************************************************/
     855             : /**                                                                **/
     856             : /**              FACTORS THROUGH CIRCLE INTEGRATION                **/
     857             : /**                                                                **/
     858             : /********************************************************************/
     859             : /* l power of 2 */
     860             : static void
     861     2400093 : fft(GEN Omega, GEN p, GEN f, long step, long l)
     862             : {
     863             :   pari_sp ltop;
     864             :   long i, l1, l2, l3, rapi, step4;
     865             :   GEN f1, f2, f3, f02, f13, g02, g13, ff;
     866             : 
     867     2400093 :   if (l == 2)
     868             :   {
     869     1319932 :     gel(f,0) = gadd(gel(p,0),gel(p,step));
     870     1319932 :     gel(f,1) = gsub(gel(p,0),gel(p,step)); return;
     871             :   }
     872     1080161 :   if (l == 4)
     873             :   {
     874      610102 :     f1 = gadd(gel(p,0),   gel(p,step<<1));
     875      610102 :     f2 = gsub(gel(p,0),   gel(p,step<<1));
     876      610102 :     f3 = gadd(gel(p,step),gel(p,3*step));
     877      610102 :     f02= gsub(gel(p,step),gel(p,3*step));
     878      610102 :     f02 = mulcxI(f02);
     879      610102 :     gel(f,0) = gadd(f1, f3);
     880      610102 :     gel(f,1) = gadd(f2, f02);
     881      610102 :     gel(f,2) = gsub(f1, f3);
     882      610102 :     gel(f,3) = gsub(f2, f02); return;
     883             :   }
     884             : 
     885      470059 :   ltop = avma;
     886      470059 :   l1 = l>>2; l2 = 2*l1; l3 = l1+l2; step4 = step<<2;
     887      470059 :   fft(Omega,p,          f,   step4,l1);
     888      470059 :   fft(Omega,p+step,     f+l1,step4,l1);
     889      470059 :   fft(Omega,p+(step<<1),f+l2,step4,l1);
     890      470059 :   fft(Omega,p+3*step,   f+l3,step4,l1);
     891             : 
     892      470059 :   ff = cgetg(l+1,t_VEC);
     893     1890937 :   for (i=0; i<l1; i++)
     894             :   {
     895     1420878 :     rapi = step*i;
     896     1420878 :     f1 = gmul(gel(Omega,rapi),    gel(f,i+l1));
     897     1420878 :     f2 = gmul(gel(Omega,rapi<<1), gel(f,i+l2));
     898     1420878 :     f3 = gmul(gel(Omega,3*rapi),  gel(f,i+l3));
     899             : 
     900     1420878 :     f02 = gadd(gel(f,i),f2);
     901     1420878 :     g02 = gsub(gel(f,i),f2);
     902     1420878 :     f13 = gadd(f1,f3);
     903     1420878 :     g13 = mulcxI(gsub(f1,f3));
     904             : 
     905     1420878 :     gel(ff,i+1)    = gadd(f02, f13);
     906     1420878 :     gel(ff,i+l1+1) = gadd(g02, g13);
     907     1420878 :     gel(ff,i+l2+1) = gsub(f02, f13);
     908     1420878 :     gel(ff,i+l3+1) = gsub(g02, g13);
     909             :   }
     910      470059 :   ff = gerepilecopy(ltop,ff);
     911      470059 :   for (i=0; i<l; i++) f[i] = ff[i+1];
     912             : }
     913             : 
     914             : GEN
     915           0 : FFTinit(long k, long prec)
     916             : {
     917           0 :   if (k <= 0) pari_err_DOMAIN("FFTinit", "k", "<=", gen_0, stoi(k));
     918           0 :   return grootsof1(1L << k, prec);
     919             : }
     920             : 
     921             : GEN
     922           0 : FFT(GEN x, GEN Omega)
     923             : {
     924           0 :   long i, l = lg(Omega), n = lg(x);
     925             :   GEN y, z;
     926           0 :   if (!is_vec_t(typ(x))) pari_err_TYPE("FFT",x);
     927           0 :   if (typ(Omega) != t_VEC) pari_err_TYPE("FFT",Omega);
     928           0 :   if (n > l) pari_err_DIM("FFT");
     929             : 
     930           0 :   if (n < l) {
     931           0 :     z = cgetg(l, t_VECSMALL); /* cf stackdummy */
     932           0 :     for (i = 1; i < n; i++) z[i] = x[i];
     933           0 :     for (     ; i < l; i++) gel(z,i) = gen_0;
     934             :   }
     935           0 :   else z = x;
     936           0 :   y = cgetg(l, t_VEC);
     937           0 :   fft(Omega+1, z+1, y+1, 1, l-1);
     938           0 :   return y;
     939             : }
     940             : 
     941             : /* returns 1 if p has only real coefficients, 0 else */
     942             : static int
     943       83321 : isreal(GEN p)
     944             : {
     945             :   long i;
     946      462739 :   for (i = lg(p)-1; i > 1; i--)
     947      407789 :     if (typ(gel(p,i)) == t_COMPLEX) return 0;
     948       54950 :   return 1;
     949             : }
     950             : 
     951             : /* x non complex */
     952             : static GEN
     953       56292 : abs_update_r(GEN x, double *mu) {
     954       56292 :   GEN y = gtofp(x, DEFAULTPREC);
     955       56292 :   double ly = mydbllogr(y); if (ly < *mu) *mu = ly;
     956       56292 :   setabssign(y); return y;
     957             : }
     958             : /* return |x|, low accuracy. Set *mu = min(log(y), *mu) */
     959             : static GEN
     960     1138288 : abs_update(GEN x, double *mu) {
     961             :   GEN y, xr, yr;
     962             :   double ly;
     963     1138288 :   if (typ(x) != t_COMPLEX) return abs_update_r(x, mu);
     964     1085816 :   xr = gel(x,1);
     965     1085816 :   yr = gel(x,2);
     966     1085816 :   if (gequal0(xr)) return abs_update_r(yr,mu);
     967     1085662 :   if (gequal0(yr)) return abs_update_r(xr,mu);
     968             :   /* have to treat 0 specially: 0E-10 + 1e-20 = 0E-10 */
     969     1081996 :   xr = gtofp(xr, DEFAULTPREC);
     970     1081996 :   yr = gtofp(yr, DEFAULTPREC);
     971     1081996 :   y = sqrtr(addrr(sqrr(xr), sqrr(yr)));
     972     1081996 :   ly = mydbllogr(y); if (ly < *mu) *mu = ly;
     973     1081996 :   return y;
     974             : }
     975             : 
     976             : static void
     977       89173 : initdft(GEN *Omega, GEN *prim, long N, long Lmax, long bit)
     978             : {
     979       89173 :   long prec = nbits2prec(bit);
     980       89173 :   *Omega = grootsof1(Lmax, prec) + 1;
     981       89173 :   *prim = rootsof1u_cx(N, prec);
     982       89173 : }
     983             : 
     984             : static void
     985       43371 : parameters(GEN p, long *LMAX, double *mu, double *gamma,
     986             :            int polreal, double param, double param2)
     987             : {
     988             :   GEN q, pc, Omega, A, RU, prim, g, TWO;
     989       43371 :   long n = degpol(p), bit, NN, K, i, j, Lmax;
     990       43371 :   pari_sp av2, av = avma;
     991             : 
     992       43371 :   bit = gexpo(p) + (long)param2+8;
     993       43371 :   Lmax = 4; while (Lmax <= n) Lmax <<= 1;
     994       43371 :   NN = (long)(param*3.14)+1; if (NN < Lmax) NN = Lmax;
     995       43371 :   K = NN/Lmax; if (K & 1) K++;
     996       43371 :   NN = Lmax*K;
     997       43371 :   if (polreal) K = K/2+1;
     998             : 
     999       43371 :   initdft(&Omega, &prim, NN, Lmax, bit);
    1000       43371 :   q = mygprec(p,bit) + 2;
    1001       43371 :   A = cgetg(Lmax+1,t_VEC); A++;
    1002       43371 :   pc= cgetg(Lmax+1,t_VEC); pc++;
    1003       43371 :   for (i=0; i <= n; i++) gel(pc,i)= gel(q,i);
    1004       43371 :   for (   ; i<Lmax; i++) gel(pc,i) = gen_0;
    1005             : 
    1006       43371 :   *mu = pariINFINITY;
    1007       43371 :   g = real_0_bit(-bit);
    1008       43371 :   TWO = real2n(1, DEFAULTPREC);
    1009       43371 :   av2 = avma;
    1010       43371 :   RU = gen_1;
    1011      168152 :   for (i=0; i<K; i++)
    1012             :   {
    1013      124781 :     if (i) {
    1014       81410 :       GEN z = RU;
    1015      518409 :       for (j=1; j<n; j++)
    1016             :       {
    1017      436999 :         gel(pc,j) = gmul(gel(q,j),z);
    1018      436999 :         z = gmul(z,RU); /* RU = prim^i, z=prim^(ij) */
    1019             :       }
    1020       81410 :       gel(pc,n) = gmul(gel(q,n),z);
    1021             :     }
    1022             : 
    1023      124781 :     fft(Omega,pc,A,1,Lmax);
    1024      140696 :     if (polreal && i>0 && i<K-1)
    1025       15915 :       for (j=0; j<Lmax; j++) g = addrr(g, divrr(TWO, abs_update(gel(A,j),mu)));
    1026             :     else
    1027      108866 :       for (j=0; j<Lmax; j++) g = addrr(g, invr(abs_update(gel(A,j),mu)));
    1028      124781 :     RU = gmul(RU, prim);
    1029      124781 :     if (gc_needed(av,1))
    1030             :     {
    1031           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"parameters");
    1032           0 :       gerepileall(av2,2, &g,&RU);
    1033             :     }
    1034             :   }
    1035       43371 :   *gamma = mydbllog2r(divru(g,NN));
    1036       43371 :   *LMAX = Lmax; set_avma(av);
    1037       43371 : }
    1038             : 
    1039             : /* NN is a multiple of Lmax */
    1040             : static void
    1041       45802 : dft(GEN p, long k, long NN, long Lmax, long bit, GEN F, GEN H, long polreal)
    1042             : {
    1043             :   GEN Omega, q, qd, pc, pd, A, B, C, RU, aux, U, W, prim, prim2;
    1044       45802 :   long n = degpol(p), i, j, K;
    1045             :   pari_sp ltop;
    1046             : 
    1047       45802 :   initdft(&Omega, &prim, NN, Lmax, bit);
    1048       45802 :   RU = cgetg(n+2,t_VEC) + 1;
    1049             : 
    1050       45802 :   K = NN/Lmax; if (polreal) K = K/2+1;
    1051       45802 :   q = mygprec(p,bit);
    1052       45802 :   qd = RgX_deriv(q);
    1053             : 
    1054       45802 :   A = cgetg(Lmax+1,t_VEC); A++;
    1055       45802 :   B = cgetg(Lmax+1,t_VEC); B++;
    1056       45802 :   C = cgetg(Lmax+1,t_VEC); C++;
    1057       45802 :   pc = cgetg(Lmax+1,t_VEC); pc++;
    1058       45802 :   pd = cgetg(Lmax+1,t_VEC); pd++;
    1059       45802 :   pc[0] = q[2];  for (i=n+1; i<Lmax; i++) gel(pc,i) = gen_0;
    1060       45802 :   pd[0] = qd[2]; for (i=n;   i<Lmax; i++) gel(pd,i) = gen_0;
    1061             : 
    1062       45802 :   ltop = avma;
    1063       45802 :   W = cgetg(k+1,t_VEC);
    1064       45802 :   U = cgetg(k+1,t_VEC);
    1065       45802 :   for (i=1; i<=k; i++) gel(W,i) = gel(U,i) = gen_0;
    1066             : 
    1067       45802 :   gel(RU,0) = gen_1;
    1068       45802 :   prim2 = gen_1;
    1069      144571 :   for (i=0; i<K; i++)
    1070             :   {
    1071       98769 :     gel(RU,1) = prim2;
    1072       98769 :     for (j=1; j<n; j++) gel(RU,j+1) = gmul(gel(RU,j),prim2);
    1073             :     /* RU[j] = prim^(ij)= prim2^j */
    1074             : 
    1075       98769 :     for (j=1; j<n; j++) gel(pd,j) = gmul(gel(qd,j+2),gel(RU,j));
    1076       98769 :     fft(Omega,pd,A,1,Lmax);
    1077       98769 :     for (j=1; j<=n; j++) gel(pc,j) = gmul(gel(q,j+2),gel(RU,j));
    1078       98769 :     fft(Omega,pc,B,1,Lmax);
    1079       98769 :     for (j=0; j<Lmax; j++) gel(C,j) = ginv(gel(B,j));
    1080       98769 :     for (j=0; j<Lmax; j++) gel(B,j) = gmul(gel(A,j),gel(C,j));
    1081       98769 :     fft(Omega,B,A,1,Lmax);
    1082       98769 :     fft(Omega,C,B,1,Lmax);
    1083             : 
    1084       98769 :     if (polreal) /* p has real coefficients */
    1085             :     {
    1086       68856 :       if (i>0 && i<K-1)
    1087             :       {
    1088       35273 :         for (j=1; j<=k; j++)
    1089             :         {
    1090       29338 :           gel(W,j) = gadd(gel(W,j), gshift(mulreal(gel(A,j+1),gel(RU,j+1)),1));
    1091       29338 :           gel(U,j) = gadd(gel(U,j), gshift(mulreal(gel(B,j),gel(RU,j)),1));
    1092             :         }
    1093             :       }
    1094             :       else
    1095             :       {
    1096      176408 :         for (j=1; j<=k; j++)
    1097             :         {
    1098      119422 :           gel(W,j) = gadd(gel(W,j), mulreal(gel(A,j+1),gel(RU,j+1)));
    1099      119422 :           gel(U,j) = gadd(gel(U,j), mulreal(gel(B,j),gel(RU,j)));
    1100             :         }
    1101             :       }
    1102             :     }
    1103             :     else
    1104             :     {
    1105       95459 :       for (j=1; j<=k; j++)
    1106             :       {
    1107       59611 :         gel(W,j) = gadd(gel(W,j), gmul(gel(A,j+1),gel(RU,j+1)));
    1108       59611 :         gel(U,j) = gadd(gel(U,j), gmul(gel(B,j),gel(RU,j)));
    1109             :       }
    1110             :     }
    1111       98769 :     prim2 = gmul(prim2,prim);
    1112       98769 :     gerepileall(ltop,3, &W,&U,&prim2);
    1113             :   }
    1114             : 
    1115      132269 :   for (i=1; i<=k; i++)
    1116             :   {
    1117       86467 :     aux=gel(W,i);
    1118       86467 :     for (j=1; j<i; j++) aux = gadd(aux, gmul(gel(W,i-j),gel(F,k+2-j)));
    1119       86467 :     gel(F,k+2-i) = gdivgs(aux,-i*NN);
    1120             :   }
    1121      132269 :   for (i=0; i<k; i++)
    1122             :   {
    1123       86467 :     aux=gel(U,k-i);
    1124       86467 :     for (j=1+i; j<k; j++) aux = gadd(aux,gmul(gel(F,2+j),gel(U,j-i)));
    1125       86467 :     gel(H,i+2) = gdivgs(aux,NN);
    1126             :   }
    1127       45802 : }
    1128             : 
    1129             : #define NEWTON_MAX 10
    1130             : static GEN
    1131      244213 : refine_H(GEN F, GEN G, GEN HH, long bit, long Sbit)
    1132             : {
    1133      244213 :   GEN H = HH, D, aux;
    1134      244213 :   pari_sp ltop = avma;
    1135             :   long error, i, bit1, bit2;
    1136             : 
    1137      244213 :   D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F)); error = gexpo(D);
    1138      244213 :   bit2 = bit + Sbit;
    1139      457678 :   for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
    1140             :   {
    1141      213465 :     if (gc_needed(ltop,1))
    1142             :     {
    1143           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"refine_H");
    1144           0 :       gerepileall(ltop,2, &D,&H);
    1145             :     }
    1146      213465 :     bit1 = -error + Sbit;
    1147      213465 :     aux = RgX_mul(mygprec(H,bit1), mygprec(D,bit1));
    1148      213465 :     aux = RgX_rem(mygprec(aux,bit1), mygprec(F,bit1));
    1149             : 
    1150      213465 :     bit1 = -error*2 + Sbit; if (bit1 > bit2) bit1 = bit2;
    1151      213465 :     H = RgX_add(mygprec(H,bit1), aux);
    1152      213465 :     D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F));
    1153      213465 :     error = gexpo(D); if (error < -bit1) error = -bit1;
    1154             :   }
    1155      244213 :   if (error > -bit/2) return NULL; /* FAIL */
    1156      244058 :   return gerepilecopy(ltop,H);
    1157             : }
    1158             : 
    1159             : /* return 0 if fails, 1 else */
    1160             : static long
    1161       45802 : refine_F(GEN p, GEN *F, GEN *G, GEN H, long bit, double gamma)
    1162             : {
    1163       45802 :   GEN f0, FF, GG, r, HH = H;
    1164       45802 :   long error, i, bit1 = 0, bit2, Sbit, Sbit2,  enh, normF, normG, n = degpol(p);
    1165       45802 :   pari_sp av = avma;
    1166             : 
    1167       45802 :   FF = *F; GG = RgX_divrem(p, FF, &r);
    1168       45802 :   error = gexpo(r); if (error <= -bit) error = 1-bit;
    1169       45802 :   normF = gexpo(FF);
    1170       45802 :   normG = gexpo(GG);
    1171       45802 :   enh = gexpo(H); if (enh < 0) enh = 0;
    1172       45802 :   Sbit = normF + 2*normG + enh + (long)(4.*log2((double)n)+gamma) + 1;
    1173       45802 :   Sbit2 = enh + 2*(normF+normG) + (long)(2.*gamma+5.*log2((double)n)) + 1;
    1174       45802 :   bit2 = bit + Sbit;
    1175      289860 :   for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
    1176             :   {
    1177      244213 :     if (bit1 == bit2 && i >= 2) { Sbit += n; Sbit2 += n; bit2 += n; }
    1178      244213 :     if (gc_needed(av,1))
    1179             :     {
    1180           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"refine_F");
    1181           0 :       gerepileall(av,4, &FF,&GG,&r,&HH);
    1182             :     }
    1183             : 
    1184      244213 :     bit1 = -error + Sbit2;
    1185      244213 :     HH = refine_H(mygprec(FF,bit1), mygprec(GG,bit1), mygprec(HH,bit1),
    1186             :                   1-error, Sbit2);
    1187      244213 :     if (!HH) return 0; /* FAIL */
    1188             : 
    1189      244058 :     bit1 = -error + Sbit;
    1190      244058 :     r = RgX_mul(mygprec(HH,bit1), mygprec(r,bit1));
    1191      244058 :     f0 = RgX_rem(mygprec(r,bit1), mygprec(FF,bit1));
    1192             : 
    1193      244058 :     bit1 = -2*error + Sbit; if (bit1 > bit2) bit1 = bit2;
    1194      244058 :     FF = gadd(mygprec(FF,bit1),f0);
    1195             : 
    1196      244058 :     bit1 = -3*error + Sbit; if (bit1 > bit2) bit1 = bit2;
    1197      244058 :     GG = RgX_divrem(mygprec(p,bit1), mygprec(FF,bit1), &r);
    1198      244058 :     error = gexpo(r); if (error < -bit1) error = -bit1;
    1199             :   }
    1200       45647 :   if (error>-bit) return 0; /* FAIL */
    1201       43371 :   *F = FF; *G = GG; return 1;
    1202             : }
    1203             : 
    1204             : /* returns F and G from the unit circle U such that |p-FG|<2^(-bit) |cd|,
    1205             : where cd is the leading coefficient of p */
    1206             : static void
    1207       43371 : split_fromU(GEN p, long k, double delta, long bit,
    1208             :             GEN *F, GEN *G, double param, double param2)
    1209             : {
    1210             :   GEN pp, FF, GG, H;
    1211       43371 :   long n = degpol(p), NN, bit2, Lmax;
    1212       43371 :   int polreal = isreal(p);
    1213             :   pari_sp ltop;
    1214             :   double mu, gamma;
    1215             : 
    1216       43371 :   pp = gdiv(p, gel(p,2+n));
    1217       43371 :   parameters(pp, &Lmax,&mu,&gamma, polreal,param,param2);
    1218             : 
    1219       43371 :   H  = cgetg(k+2,t_POL); H[1] = p[1];
    1220       43371 :   FF = cgetg(k+3,t_POL); FF[1]= p[1];
    1221       43371 :   gel(FF,k+2) = gen_1;
    1222             : 
    1223       43371 :   NN = (long)(0.5/delta); NN |= 1; if (NN < 2) NN = 2;
    1224       43371 :   NN *= Lmax; ltop = avma;
    1225             :   for(;;)
    1226             :   {
    1227       48233 :     bit2 = (long)(((double)NN*delta-mu)/M_LN2) + gexpo(pp) + 8;
    1228       45802 :     dft(pp, k, NN, Lmax, bit2, FF, H, polreal);
    1229       45802 :     if (refine_F(pp,&FF,&GG,H,bit,gamma)) break;
    1230        2431 :     NN <<= 1; set_avma(ltop);
    1231             :   }
    1232       43371 :   *G = gmul(GG,gel(p,2+n)); *F = FF;
    1233       43371 : }
    1234             : 
    1235             : static void
    1236       43371 : optimize_split(GEN p, long k, double delta, long bit,
    1237             :             GEN *F, GEN *G, double param, double param2)
    1238             : {
    1239       43371 :   long n = degpol(p);
    1240             :   GEN FF, GG;
    1241             : 
    1242       43371 :   if (k <= n/2)
    1243       32121 :     split_fromU(p,k,delta,bit,F,G,param,param2);
    1244             :   else
    1245             :   {
    1246       11250 :     split_fromU(RgX_recip_shallow(p),n-k,delta,bit,&FF,&GG,param,param2);
    1247       11250 :     *F = RgX_recip_shallow(GG);
    1248       11250 :     *G = RgX_recip_shallow(FF);
    1249             :   }
    1250       43371 : }
    1251             : 
    1252             : /********************************************************************/
    1253             : /**                                                                **/
    1254             : /**               SEARCH FOR SEPARATING CIRCLE                     **/
    1255             : /**                                                                **/
    1256             : /********************************************************************/
    1257             : 
    1258             : /* return p(2^e*x) *2^(-n*e) */
    1259             : static void
    1260           0 : scalepol2n(GEN p, long e)
    1261             : {
    1262           0 :   long i,n=lg(p)-1;
    1263           0 :   for (i=2; i<=n; i++) gel(p,i) = gmul2n(gel(p,i),(i-n)*e);
    1264           0 : }
    1265             : 
    1266             : /* returns p(x/R)*R^n */
    1267             : static GEN
    1268      369813 : scalepol(GEN p, GEN R, long bit)
    1269             : {
    1270             :   GEN q,aux,gR;
    1271             :   long i;
    1272             : 
    1273      369813 :   aux = gR = mygprec(R,bit); q = mygprec(p,bit);
    1274     1652323 :   for (i=lg(p)-2; i>=2; i--)
    1275             :   {
    1276     1282510 :     gel(q,i) = gmul(aux,gel(q,i));
    1277     1282510 :     aux = gmul(aux,gR);
    1278             :   }
    1279      369813 :   return q;
    1280             : }
    1281             : 
    1282             : /* return (conj(a)X-1)^n * p[ (X-a) / (conj(a)X-1) ] */
    1283             : static GEN
    1284      119850 : conformal_pol(GEN p, GEN a)
    1285             : {
    1286      119850 :   GEN z, r, ma = gneg(a), ca = conj_i(a);
    1287      119850 :   long n = degpol(p), i;
    1288      119850 :   pari_sp av = avma;
    1289             : 
    1290      119850 :   z = mkpoln(2, ca, gen_m1);
    1291      119850 :   r = scalarpol(gel(p,2+n), 0);
    1292      406598 :   for (i=n-1; ; i--)
    1293             :   {
    1294      693346 :     r = RgX_addmulXn_shallow(r, gmul(ma,r), 1); /* r *= (X - a) */
    1295      406598 :     r = gadd(r, gmul(z, gel(p,2+i)));
    1296      526448 :     if (i == 0) return gerepileupto(av, r);
    1297      286748 :     z = RgX_addmulXn_shallow(gmul(z,ca), gneg(z), 1); /* z *= conj(a)X - 1 */
    1298      286748 :     if (gc_needed(av,2))
    1299             :     {
    1300           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"conformal_pol");
    1301           0 :       gerepileall(av,2, &r,&z);
    1302             :     }
    1303             :   }
    1304             : }
    1305             : 
    1306             : static const double UNDEF = -100000.;
    1307             : 
    1308             : static double
    1309       43371 : logradius(double *radii, GEN p, long k, double aux, double *delta)
    1310             : {
    1311       43371 :   long i, n = degpol(p);
    1312             :   double lrho, lrmin, lrmax;
    1313       43371 :   if (k > 1)
    1314             :   {
    1315       27356 :     i = k-1; while (i>0 && radii[i] == UNDEF) i--;
    1316       27356 :     lrmin = logpre_modulus(p,k,aux, radii[i], radii[k]);
    1317             :   }
    1318             :   else /* k=1 */
    1319       16015 :     lrmin = logmin_modulus(p,aux);
    1320       43371 :   radii[k] = lrmin;
    1321             : 
    1322       43371 :   if (k+1<n)
    1323             :   {
    1324       35383 :     i = k+2; while (i<=n && radii[i] == UNDEF) i++;
    1325       35383 :     lrmax = logpre_modulus(p,k+1,aux, radii[k+1], radii[i]);
    1326             :   }
    1327             :   else /* k+1=n */
    1328        7988 :     lrmax = logmax_modulus(p,aux);
    1329       43371 :   radii[k+1] = lrmax;
    1330             : 
    1331       43371 :   lrho = radii[k];
    1332       97558 :   for (i=k-1; i>=1; i--)
    1333             :   {
    1334       54187 :     if (radii[i] == UNDEF || radii[i] > lrho)
    1335       37815 :       radii[i] = lrho;
    1336             :     else
    1337       16372 :       lrho = radii[i];
    1338             :   }
    1339       43371 :   lrho = radii[k+1];
    1340      180470 :   for (i=k+1; i<=n; i++)
    1341             :   {
    1342      137099 :     if (radii[i] == UNDEF || radii[i] < lrho)
    1343       80074 :       radii[i] = lrho;
    1344             :     else
    1345       57025 :       lrho = radii[i];
    1346             :   }
    1347       43371 :   *delta = (lrmax - lrmin) / 2;
    1348       43371 :   if (*delta > 1.) *delta = 1.;
    1349       43371 :   return (lrmin + lrmax) / 2;
    1350             : }
    1351             : 
    1352             : static void
    1353       43371 : update_radius(long n, double *radii, double lrho, double *par, double *par2)
    1354             : {
    1355       43371 :   double t, param = 0., param2 = 0.;
    1356             :   long i;
    1357      278028 :   for (i=1; i<=n; i++)
    1358             :   {
    1359      234657 :     radii[i] -= lrho;
    1360      234657 :     t = fabs(rtodbl( invr(subsr(1, dblexp(radii[i]))) ));
    1361      234657 :     param += t; if (t > 1.) param2 += log2(t);
    1362             :   }
    1363       43371 :   *par = param; *par2 = param2;
    1364       43371 : }
    1365             : 
    1366             : /* apply the conformal mapping then split from U */
    1367             : static void
    1368       39950 : conformal_mapping(double *radii, GEN ctr, GEN p, long k, long bit,
    1369             :                   double aux, GEN *F,GEN *G)
    1370             : {
    1371       39950 :   long bit2, n = degpol(p), i;
    1372       39950 :   pari_sp ltop = avma, av;
    1373             :   GEN q, FF, GG, a, R;
    1374             :   double lrho, delta, param, param2;
    1375             :   /* n * (2.*log2(2.732)+log2(1.5)) + 1 */
    1376       39950 :   bit2 = bit + (long)(n*3.4848775) + 1;
    1377       39950 :   a = sqrtr_abs( utor(3, 2*MEDDEFAULTPREC - 2) );
    1378       39950 :   a = divrs(a, -6);
    1379       39950 :   a = gmul(mygprec(a,bit2), mygprec(ctr,bit2)); /* a = -ctr/2sqrt(3) */
    1380             : 
    1381       39950 :   av = avma;
    1382       39950 :   q = conformal_pol(mygprec(p,bit2), a);
    1383      243249 :   for (i=1; i<=n; i++)
    1384      203299 :     if (radii[i] != UNDEF) /* update array radii */
    1385             :     {
    1386      144449 :       pari_sp av2 = avma;
    1387      144449 :       GEN t, r = dblexp(radii[i]), r2 = sqrr(r);
    1388             :       /* 2(r^2 - 1) / (r^2 - 3(r-1)) */
    1389      144449 :       t = divrr(shiftr((subrs(r2,1)),1), subrr(r2, mulur(3,subrs(r,1))));
    1390      144449 :       radii[i] = mydbllogr(addsr(1,t)) / 2;
    1391      144449 :       set_avma(av2);
    1392             :     }
    1393       39950 :   lrho = logradius(radii, q,k,aux/10., &delta);
    1394       39950 :   update_radius(n, radii, lrho, &param, &param2);
    1395             : 
    1396       39950 :   bit2 += (long)(n * fabs(lrho)/M_LN2 + 1.);
    1397       39950 :   R = mygprec(dblexp(-lrho), bit2);
    1398       39950 :   q = scalepol(q,R,bit2);
    1399       39950 :   gerepileall(av,2, &q,&R);
    1400             : 
    1401       39950 :   optimize_split(q,k,delta,bit2,&FF,&GG,param,param2);
    1402       39950 :   bit2 += n; R = invr(R);
    1403       39950 :   FF = scalepol(FF,R,bit2);
    1404       39950 :   GG = scalepol(GG,R,bit2);
    1405             : 
    1406       39950 :   a = mygprec(a,bit2);
    1407       39950 :   FF = conformal_pol(FF,a);
    1408       39950 :   GG = conformal_pol(GG,a);
    1409             : 
    1410       39950 :   a = invr(subsr(1, gnorm(a)));
    1411       39950 :   FF = RgX_Rg_mul(FF, powru(a,k));
    1412       39950 :   GG = RgX_Rg_mul(GG, powru(a,n-k));
    1413             : 
    1414       39950 :   *F = mygprec(FF,bit+n);
    1415       39950 :   *G = mygprec(GG,bit+n); gerepileall(ltop,2, F,G);
    1416       39950 : }
    1417             : 
    1418             : /* split p, this time without scaling. returns in F and G two polynomials
    1419             :  * such that |p-FG|< 2^(-bit)|p| */
    1420             : static void
    1421       43371 : split_2(GEN p, long bit, GEN ctr, double thickness, GEN *F, GEN *G)
    1422             : {
    1423             :   GEN q, FF, GG, R;
    1424             :   double aux, delta, param, param2;
    1425       43371 :   long n = degpol(p), i, j, k, bit2;
    1426             :   double lrmin, lrmax, lrho, *radii;
    1427             : 
    1428       43371 :   radii = (double*) stack_malloc_align((n+1) * sizeof(double), sizeof(double));
    1429             : 
    1430       43371 :   for (i=2; i<n; i++) radii[i] = UNDEF;
    1431       43371 :   aux = thickness/(double)(4 * n);
    1432       43371 :   lrmin = logmin_modulus(p, aux);
    1433       43371 :   lrmax = logmax_modulus(p, aux);
    1434       43371 :   radii[1] = lrmin;
    1435       43371 :   radii[n] = lrmax;
    1436       43371 :   i = 1; j = n;
    1437       43371 :   lrho = (lrmin + lrmax) / 2;
    1438       43371 :   k = dual_modulus(p, lrho, aux, 1);
    1439       43371 :   if (5*k < n || (n < 2*k && 5*k < 4*n))
    1440        8942 :     { lrmax = lrho; j=k+1; radii[j] = lrho; }
    1441             :   else
    1442       34429 :     { lrmin = lrho; i=k;   radii[i] = lrho; }
    1443      136557 :   while (j > i+1)
    1444             :   {
    1445       49815 :     if (i+j == n+1)
    1446       22655 :       lrho = (lrmin + lrmax) / 2;
    1447             :     else
    1448             :     {
    1449       27160 :       double kappa = 2. - log(1. + minss(i,n-j)) / log(1. + minss(j,n-i));
    1450       27160 :       if (i+j < n+1) lrho = lrmax * kappa + lrmin;
    1451       20287 :       else           lrho = lrmin * kappa + lrmax;
    1452       27160 :       lrho /= 1+kappa;
    1453             :     }
    1454       49815 :     aux = (lrmax - lrmin) / (4*(j-i));
    1455       49815 :     k = dual_modulus(p, lrho, aux, minss(i,n+1-j));
    1456       49815 :     if (k-i < j-k-1 || (k-i == j-k-1 && 2*k > n))
    1457       34076 :       { lrmax = lrho; j=k+1; radii[j] = lrho - aux; }
    1458             :     else
    1459       15739 :       { lrmin = lrho; i=k;   radii[i] = lrho + aux; }
    1460             :   }
    1461       43371 :   aux = lrmax - lrmin;
    1462             : 
    1463       43371 :   if (ctr)
    1464             :   {
    1465       39950 :     lrho = (lrmax + lrmin) / 2;
    1466      243249 :     for (i=1; i<=n; i++)
    1467      203299 :       if (radii[i] != UNDEF) radii[i] -= lrho;
    1468             : 
    1469       39950 :     bit2 = bit + (long)(n * fabs(lrho)/M_LN2 + 1.);
    1470       39950 :     R = mygprec(dblexp(-lrho), bit2);
    1471       39950 :     q = scalepol(p,R,bit2);
    1472       39950 :     conformal_mapping(radii, ctr, q, k, bit2, aux, &FF, &GG);
    1473             :   }
    1474             :   else
    1475             :   {
    1476        3421 :     lrho = logradius(radii, p, k, aux/10., &delta);
    1477        3421 :     update_radius(n, radii, lrho, &param, &param2);
    1478             : 
    1479        3421 :     bit2 = bit + (long)(n * fabs(lrho)/M_LN2 + 1.);
    1480        3421 :     R = mygprec(dblexp(-lrho), bit2);
    1481        3421 :     q = scalepol(p,R,bit2);
    1482        3421 :     optimize_split(q, k, delta, bit2, &FF, &GG, param, param2);
    1483             :   }
    1484       43371 :   bit  += n;
    1485       43371 :   bit2 += n; R = invr(mygprec(R,bit2));
    1486       43371 :   *F = mygprec(scalepol(FF,R,bit2), bit);
    1487       43371 :   *G = mygprec(scalepol(GG,R,bit2), bit);
    1488       43371 : }
    1489             : 
    1490             : /* procedure corresponding to steps 5,6,.. page 44 in RR n. 1852 */
    1491             : /* put in F and G two polynomial such that |p-FG|<2^(-bit)|p|
    1492             :  * where the maximum modulus of the roots of p is <=1.
    1493             :  * Assume sum of roots is 0. */
    1494             : static void
    1495       39950 : split_1(GEN p, long bit, GEN *F, GEN *G)
    1496             : {
    1497       39950 :   long i, imax, n = degpol(p), polreal = isreal(p), ep = gexpo(p), bit2 = bit+n;
    1498             :   GEN ctr, q, qq, FF, GG, v, gr, r, newq;
    1499             :   double lrmin, lrmax, lthick;
    1500       39950 :   const double LOG3 = 1.098613;
    1501             : 
    1502       39950 :   lrmax = logmax_modulus(p, 0.01);
    1503       39950 :   gr = mygprec(dblexp(-lrmax), bit2);
    1504       39950 :   q = scalepol(p,gr,bit2);
    1505             : 
    1506       39950 :   bit2 = bit + gexpo(q) - ep + (long)((double)n*2.*log2(3.)+1);
    1507       39950 :   v = cgetg(5,t_VEC);
    1508       39950 :   gel(v,1) = gen_2;
    1509       39950 :   gel(v,2) = gen_m2;
    1510       39950 :   gel(v,3) = mkcomplex(gen_0, gel(v,1));
    1511       39950 :   gel(v,4) = mkcomplex(gen_0, gel(v,2));
    1512       39950 :   q = mygprec(q,bit2); lthick = 0;
    1513       39950 :   newq = ctr = NULL; /* -Wall */
    1514       39950 :   imax = polreal? 3: 4;
    1515       77970 :   for (i=1; i<=imax; i++)
    1516             :   {
    1517       76528 :     qq = RgX_translate(q, gel(v,i));
    1518       76528 :     lrmin = logmin_modulus(qq,0.05);
    1519       76528 :     if (LOG3 > lrmin + lthick)
    1520             :     {
    1521       74431 :       double lquo = logmax_modulus(qq,0.05) - lrmin;
    1522       74431 :       if (lquo > lthick) { lthick = lquo; newq = qq; ctr = gel(v,i); }
    1523             :     }
    1524       76528 :     if (lthick > M_LN2) break;
    1525       43389 :     if (polreal && i==2 && lthick > LOG3 - M_LN2) break;
    1526             :   }
    1527       39950 :   bit2 = bit + gexpo(newq) - ep + (long)(n*LOG3/M_LN2 + 1);
    1528       39950 :   split_2(newq, bit2, ctr, lthick, &FF, &GG);
    1529       39950 :   r = gneg(mygprec(ctr,bit2));
    1530       39950 :   FF = RgX_translate(FF,r);
    1531       39950 :   GG = RgX_translate(GG,r);
    1532             : 
    1533       39950 :   gr = invr(gr); bit2 = bit - ep + gexpo(FF)+gexpo(GG);
    1534       39950 :   *F = scalepol(FF,gr,bit2);
    1535       39950 :   *G = scalepol(GG,gr,bit2);
    1536       39950 : }
    1537             : 
    1538             : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|,
    1539             : where the maximum modulus of the roots of p is < 0.5 */
    1540             : static int
    1541       40085 : split_0_2(GEN p, long bit, GEN *F, GEN *G)
    1542             : {
    1543             :   GEN q, b;
    1544       40085 :   long n = degpol(p), k, bit2, eq;
    1545       40085 :   double aux0 = dbllog2(gel(p,n+2)); /* != -oo */
    1546       40085 :   double aux1 = dbllog2(gel(p,n+1)), aux;
    1547             : 
    1548       40085 :   if (aux1 == -pariINFINITY) /* p1 = 0 */
    1549        1041 :     aux = 0;
    1550             :   else
    1551             :   {
    1552       39044 :     aux = aux1 - aux0; /* log2(p1/p0) */
    1553             :     /* beware double overflow */
    1554       39044 :     if (aux >= 0 && (aux > 1e4 || exp2(aux) > 2.5*n)) return 0;
    1555       39044 :     aux = (aux < -300)? 0.: n*log2(1 + exp2(aux)/(double)n);
    1556             :   }
    1557       40085 :   bit2 = bit+1 + (long)(log2((double)n) + aux);
    1558       40085 :   q = mygprec(p,bit2);
    1559       40085 :   if (aux1 == -pariINFINITY) b = NULL;
    1560             :   else
    1561             :   {
    1562       39044 :     b = gdivgs(gdiv(gel(q,n+1),gel(q,n+2)),-n);
    1563       39044 :     q = RgX_translate(q,b);
    1564             :   }
    1565       40085 :   gel(q,n+1) = gen_0; eq = gexpo(q);
    1566       40085 :   k = 0;
    1567       80535 :   while (k <= n/2 && (- gexpo(gel(q,k+2)) > bit2 + 2*(n-k) + eq
    1568       40328 :                       || gequal0(gel(q,k+2)))) k++;
    1569       40085 :   if (k > 0)
    1570             :   {
    1571         135 :     if (k > n/2) k = n/2;
    1572         135 :     bit2 += k<<1;
    1573         135 :     *F = pol_xn(k, 0);
    1574         135 :     *G = RgX_shift_shallow(q, -k);
    1575             :   }
    1576             :   else
    1577             :   {
    1578       39950 :     split_1(q,bit2,F,G);
    1579       39950 :     bit2 = bit + gexpo(*F) + gexpo(*G) - gexpo(p) + (long)aux+1;
    1580       39950 :     *F = mygprec(*F,bit2);
    1581             :   }
    1582       40085 :   *G = mygprec(*G,bit2);
    1583       40085 :   if (b)
    1584             :   {
    1585       39044 :     GEN mb = mygprec(gneg(b), bit2);
    1586       39044 :     *F = RgX_translate(*F, mb);
    1587       39044 :     *G = RgX_translate(*G, mb);
    1588             :   }
    1589       40085 :   return 1;
    1590             : }
    1591             : 
    1592             : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|.
    1593             :  * Assume max_modulus(p) < 2 */
    1594             : static void
    1595       40085 : split_0_1(GEN p, long bit, GEN *F, GEN *G)
    1596             : {
    1597             :   GEN FF, GG;
    1598             :   long n, bit2, normp;
    1599             : 
    1600       40085 :   if  (split_0_2(p,bit,F,G)) return;
    1601             : 
    1602           0 :   normp = gexpo(p);
    1603           0 :   scalepol2n(p,2); /* p := 4^(-n) p(4*x) */
    1604           0 :   n = degpol(p); bit2 = bit + 2*n + gexpo(p) - normp;
    1605           0 :   split_1(mygprec(p,bit2), bit2,&FF,&GG);
    1606           0 :   scalepol2n(FF,-2);
    1607           0 :   scalepol2n(GG,-2); bit2 = bit + gexpo(FF) + gexpo(GG) - normp;
    1608           0 :   *F = mygprec(FF,bit2);
    1609           0 :   *G = mygprec(GG,bit2);
    1610             : }
    1611             : 
    1612             : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P| */
    1613             : static void
    1614       43506 : split_0(GEN p, long bit, GEN *F, GEN *G)
    1615             : {
    1616       43506 :   const double LOG1_9 = 0.6418539;
    1617       43506 :   long n = degpol(p), k = 0;
    1618             :   GEN q;
    1619             : 
    1620       43506 :   while (gexpo(gel(p,k+2)) < -bit && k <= n/2) k++;
    1621       43506 :   if (k > 0)
    1622             :   {
    1623           0 :     if (k > n/2) k = n/2;
    1624           0 :     *F = pol_xn(k, 0);
    1625           0 :     *G = RgX_shift_shallow(p, -k);
    1626             :   }
    1627             :   else
    1628             :   {
    1629       43506 :     double lr = logmax_modulus(p, 0.05);
    1630       43506 :     if (lr < LOG1_9) split_0_1(p, bit, F, G);
    1631             :     else
    1632             :     {
    1633       33418 :       q = RgX_recip_shallow(p);
    1634       33418 :       lr = logmax_modulus(q,0.05);
    1635       33418 :       if (lr < LOG1_9)
    1636             :       {
    1637       29997 :         split_0_1(q, bit, F, G);
    1638       29997 :         *F = RgX_recip_shallow(*F);
    1639       29997 :         *G = RgX_recip_shallow(*G);
    1640             :       }
    1641             :       else
    1642        3421 :         split_2(p,bit,NULL, 1.2837,F,G);
    1643             :     }
    1644             :   }
    1645       43506 : }
    1646             : 
    1647             : /********************************************************************/
    1648             : /**                                                                **/
    1649             : /**                ERROR ESTIMATE FOR THE ROOTS                    **/
    1650             : /**                                                                **/
    1651             : /********************************************************************/
    1652             : 
    1653             : static GEN
    1654      175281 : root_error(long n, long k, GEN roots_pol, long err, GEN shatzle)
    1655             : {
    1656      175281 :   GEN rho, d, eps, epsbis, eps2, aux, rap = NULL;
    1657             :   long i, j;
    1658             : 
    1659      175281 :   d = cgetg(n+1,t_VEC);
    1660     1860738 :   for (i=1; i<=n; i++)
    1661             :   {
    1662     1685457 :     if (i!=k)
    1663             :     {
    1664     1510176 :       aux = gsub(gel(roots_pol,i), gel(roots_pol,k));
    1665     1510176 :       gel(d,i) = gabs(mygprec(aux,31), DEFAULTPREC);
    1666             :     }
    1667             :   }
    1668      175281 :   rho = gabs(mygprec(gel(roots_pol,k),31), DEFAULTPREC);
    1669      175281 :   if (expo(rho) < 0) rho = real_1(DEFAULTPREC);
    1670      175281 :   eps = mulrr(rho, shatzle);
    1671      175281 :   aux = shiftr(powru(rho,n), err);
    1672             : 
    1673      543094 :   for (j=1; j<=2 || (j<=5 && cmprr(rap, dbltor(1.2)) > 0); j++)
    1674             :   {
    1675      367813 :     GEN prod = NULL; /* 1. */
    1676      367813 :     long m = n;
    1677      367813 :     epsbis = mulrr(eps, dbltor(1.25));
    1678     4423671 :     for (i=1; i<=n; i++)
    1679             :     {
    1680     4055858 :       if (i != k && cmprr(gel(d,i),epsbis) > 0)
    1681             :       {
    1682     3656871 :         GEN dif = subrr(gel(d,i),eps);
    1683     3656871 :         prod = prod? mulrr(prod, dif): dif;
    1684     3656871 :         m--;
    1685             :       }
    1686             :     }
    1687      367813 :     eps2 = prod? divrr(aux, prod): aux;
    1688      367813 :     if (m > 1) eps2 = sqrtnr(shiftr(eps2, 2*m-2), m);
    1689      367813 :     rap = divrr(eps,eps2); eps = eps2;
    1690             :   }
    1691      175281 :   return eps;
    1692             : }
    1693             : 
    1694             : /* round a complex or real number x to an absolute value of 2^(-bit) */
    1695             : static GEN
    1696      430421 : mygprec_absolute(GEN x, long bit)
    1697             : {
    1698             :   long e;
    1699             :   GEN y;
    1700             : 
    1701      430421 :   switch(typ(x))
    1702             :   {
    1703             :     case t_REAL:
    1704      280697 :       e = expo(x) + bit;
    1705      280697 :       return (e <= 0 || !signe(x))? real_0_bit(-bit): rtor(x, nbits2prec(e));
    1706             :     case t_COMPLEX:
    1707      129384 :       if (gexpo(gel(x,2)) < -bit) return mygprec_absolute(gel(x,1),bit);
    1708      125756 :       y = cgetg(3,t_COMPLEX);
    1709      125756 :       gel(y,1) = mygprec_absolute(gel(x,1),bit);
    1710      125756 :       gel(y,2) = mygprec_absolute(gel(x,2),bit);
    1711      125756 :       return y;
    1712       20340 :     default: return x;
    1713             :   }
    1714             : }
    1715             : 
    1716             : static long
    1717       43554 : a_posteriori_errors(GEN p, GEN roots_pol, long err)
    1718             : {
    1719       43554 :   long i, n = degpol(p), e_max = -(long)EXPOBITS;
    1720             :   GEN sigma, shatzle;
    1721             : 
    1722       43554 :   err += (long)log2((double)n) + 1;
    1723       43554 :   if (err > -2) return 0;
    1724       43554 :   sigma = real2n(-err, LOWDEFAULTPREC);
    1725             :   /*  2 / ((s - 1)^(1/n) - 1) */
    1726       43554 :   shatzle = divur(2, subrs(sqrtnr(subrs(sigma,1),n), 1));
    1727      218835 :   for (i=1; i<=n; i++)
    1728             :   {
    1729      175281 :     pari_sp av = avma;
    1730      175281 :     GEN x = root_error(n,i,roots_pol,err,shatzle);
    1731      175281 :     long e = gexpo(x);
    1732      175281 :     set_avma(av); if (e > e_max) e_max = e;
    1733      175281 :     gel(roots_pol,i) = mygprec_absolute(gel(roots_pol,i), -e);
    1734             :   }
    1735       43554 :   return e_max;
    1736             : }
    1737             : 
    1738             : /********************************************************************/
    1739             : /**                                                                **/
    1740             : /**                           MAIN                                 **/
    1741             : /**                                                                **/
    1742             : /********************************************************************/
    1743             : static GEN
    1744      140022 : append_clone(GEN r, GEN a) { a = gclone(a); vectrunc_append(r, a); return a; }
    1745             : 
    1746             : /* put roots in placeholder roots_pol so that |P - L_1...L_n| < 2^(-bit)|P|
    1747             :  * returns prod (x-roots_pol[i]) */
    1748             : static GEN
    1749      130566 : split_complete(GEN p, long bit, GEN roots_pol)
    1750             : {
    1751      130566 :   long n = degpol(p);
    1752             :   pari_sp ltop;
    1753             :   GEN p1, F, G, a, b, m1, m2;
    1754             : 
    1755      130566 :   if (n == 1)
    1756             :   {
    1757       34098 :     a = gneg_i(gdiv(gel(p,2), gel(p,3)));
    1758       34098 :     (void)append_clone(roots_pol,a); return p;
    1759             :   }
    1760       96468 :   ltop = avma;
    1761       96468 :   if (n == 2)
    1762             :   {
    1763       52962 :     F = gsub(gsqr(gel(p,3)), gmul2n(gmul(gel(p,2),gel(p,4)), 2));
    1764       52962 :     F = gsqrt(F, nbits2prec(bit));
    1765       52962 :     p1 = ginv(gmul2n(gel(p,4),1));
    1766       52962 :     a = gneg_i(gmul(gadd(F,gel(p,3)), p1));
    1767       52962 :     b =        gmul(gsub(F,gel(p,3)), p1);
    1768       52962 :     a = append_clone(roots_pol,a);
    1769       52962 :     b = append_clone(roots_pol,b); set_avma(ltop);
    1770       52962 :     a = mygprec(a, 3*bit);
    1771       52962 :     b = mygprec(b, 3*bit);
    1772       52962 :     return gmul(gel(p,4), mkpoln(3, gen_1, gneg(gadd(a,b)), gmul(a,b)));
    1773             :   }
    1774       43506 :   split_0(p,bit,&F,&G);
    1775       43506 :   m1 = split_complete(F,bit,roots_pol);
    1776       43506 :   m2 = split_complete(G,bit,roots_pol);
    1777       43506 :   return gerepileupto(ltop, gmul(m1,m2));
    1778             : }
    1779             : 
    1780             : static GEN
    1781      690213 : quicktofp(GEN x)
    1782             : {
    1783      690213 :   const long prec = DEFAULTPREC;
    1784      690213 :   switch(typ(x))
    1785             :   {
    1786      671897 :     case t_INT: return itor(x, prec);
    1787        7198 :     case t_REAL: return rtor(x, prec);
    1788           0 :     case t_FRAC: return fractor(x, prec);
    1789             :     case t_COMPLEX: {
    1790       11118 :       GEN a = gel(x,1), b = gel(x,2);
    1791             :       /* avoid problem with 0, e.g. x = 0 + I*1e-100. We don't want |x| = 0. */
    1792       11118 :       if (isintzero(a)) return cxcompotor(b, prec);
    1793       11076 :       if (isintzero(b)) return cxcompotor(a, prec);
    1794       11076 :       a = cxcompotor(a, prec);
    1795       11076 :       b = cxcompotor(b, prec); return sqrtr(addrr(sqrr(a), sqrr(b)));
    1796             :     }
    1797           0 :     default: pari_err_TYPE("quicktofp",x);
    1798             :       return NULL;/*LCOV_EXCL_LINE*/
    1799             :   }
    1800             : 
    1801             : }
    1802             : 
    1803             : /* bound log_2 |largest root of p| (Fujiwara's bound) */
    1804             : double
    1805      182823 : fujiwara_bound(GEN p)
    1806             : {
    1807      182823 :   pari_sp av = avma;
    1808      182823 :   long i, n = degpol(p);
    1809             :   GEN cc;
    1810             :   double loglc, Lmax;
    1811             : 
    1812      182823 :   if (n <= 0) pari_err_CONSTPOL("fujiwara_bound");
    1813      182823 :   loglc = mydbllog2r( quicktofp(gel(p,n+2)) ); /* log_2 |lc(p)| */
    1814      182823 :   cc = gel(p, 2);
    1815      182823 :   if (gequal0(cc))
    1816       21812 :     Lmax = -pariINFINITY-1;
    1817             :   else
    1818      161011 :     Lmax = (mydbllog2r(quicktofp(cc)) - loglc - 1) / n;
    1819      701329 :   for (i = 1; i < n; i++)
    1820             :   {
    1821      518506 :     GEN y = gel(p,i+2);
    1822             :     double L;
    1823      518506 :     if (gequal0(y)) continue;
    1824      346379 :     L = (mydbllog2r(quicktofp(y)) - loglc) / (n-i);
    1825      346379 :     if (L > Lmax) Lmax = L;
    1826             :   }
    1827      182823 :   return gc_double(av, Lmax+1);
    1828             : }
    1829             : 
    1830             : /* Fujiwara's bound, real roots. Based on the following remark: if
    1831             :  *   p = x^n + sum a_i x^i and q = x^n + sum min(a_i,0)x^i
    1832             :  * then for all x >= 0, p(x) >= q(x). Thus any bound for the (positive) roots
    1833             :  * of q is a bound for the positive roots of p. */
    1834             : double
    1835       68779 : fujiwara_bound_real(GEN p, long sign)
    1836             : {
    1837       68779 :   pari_sp av = avma;
    1838             :   GEN x;
    1839       68779 :   long n = degpol(p), i, signodd, signeven;
    1840       68779 :   if (n <= 0) pari_err_CONSTPOL("fujiwara_bound");
    1841       68779 :   x = shallowcopy(p);
    1842       68779 :   if (gsigne(gel(x, n+2)) > 0)
    1843       68779 :   { signeven = 1; signodd = sign; }
    1844             :   else
    1845           0 :   { signeven = -1; signodd = -sign; }
    1846      330316 :   for (i = 0; i < n; i++)
    1847             :   {
    1848      261537 :     if ((n - i) % 2)
    1849      133394 :     { if (gsigne(gel(x, i+2)) == signodd ) gel(x, i+2) = gen_0; }
    1850             :     else
    1851      128143 :     { if (gsigne(gel(x, i+2)) == signeven) gel(x, i+2) = gen_0; }
    1852             :   }
    1853       68779 :   return gc_double(av, fujiwara_bound(x));
    1854             : }
    1855             : 
    1856             : static GEN
    1857      242465 : mygprecrc_special(GEN x, long prec, long e)
    1858             : {
    1859             :   GEN y;
    1860      242465 :   switch(typ(x))
    1861             :   {
    1862             :     case t_REAL:
    1863       31449 :       if (!signe(x)) return real_0_bit(minss(e, expo(x)));
    1864       28964 :       return (prec > realprec(x))? rtor(x, prec): x;
    1865             :     case t_COMPLEX:
    1866       11815 :       y = cgetg(3,t_COMPLEX);
    1867       11815 :       gel(y,1) = mygprecrc_special(gel(x,1),prec,e);
    1868       11815 :       gel(y,2) = mygprecrc_special(gel(x,2),prec,e);
    1869       11815 :       return y;
    1870      199201 :     default: return x;
    1871             :   }
    1872             : }
    1873             : 
    1874             : /* like mygprec but keep at least the same precision as before */
    1875             : static GEN
    1876       43554 : mygprec_special(GEN x, long bit)
    1877             : {
    1878             :   long lx, i, e, prec;
    1879             :   GEN y;
    1880             : 
    1881       43554 :   if (bit < 0) bit = 0; /* should not happen */
    1882       43554 :   e = gexpo(x) - bit;
    1883       43554 :   prec = nbits2prec(bit);
    1884       43554 :   switch(typ(x))
    1885             :   {
    1886             :     case t_POL:
    1887       43554 :       y = cgetg_copy(x, &lx); y[1] = x[1];
    1888       43554 :       for (i=2; i<lx; i++) gel(y,i) = mygprecrc_special(gel(x,i),prec,e);
    1889       43554 :       break;
    1890             : 
    1891           0 :     default: y = mygprecrc_special(x,prec,e);
    1892             :   }
    1893       43554 :   return y;
    1894             : }
    1895             : 
    1896             : static GEN
    1897       29278 : fix_roots1(GEN r)
    1898             : {
    1899       29278 :   long i, l = lg(r);
    1900       29278 :   GEN allr = cgetg(l, t_VEC);
    1901      141855 :   for (i=1; i<l; i++)
    1902             :   {
    1903      112577 :     GEN t = gel(r,i);
    1904      112577 :     gel(allr,i) = gcopy(t); gunclone(t);
    1905             :   }
    1906       29278 :   return allr;
    1907             : }
    1908             : static GEN
    1909       43554 : fix_roots(GEN r, GEN *m, long h, long bit)
    1910             : {
    1911             :   long i, j, k, l, prec;
    1912             :   GEN allr, ro1;
    1913       43554 :   if (h == 1) return fix_roots1(r);
    1914       14276 :   prec = nbits2prec(bit);
    1915       14276 :   ro1 = grootsof1(h, prec) + 1;
    1916       14276 :   l = lg(r)-1;
    1917       14276 :   allr = cgetg(h*l+1, t_VEC);
    1918       41721 :   for (k=1,i=1; i<=l; i++)
    1919             :   {
    1920       27445 :     GEN p2, p1 = gel(r,i);
    1921       27445 :     p2 = (h == 2)? gsqrt(p1, prec): gsqrtn(p1, utoipos(h), NULL, prec);
    1922       27445 :     for (j=0; j<h; j++) gel(allr,k++) = gmul(p2, gel(ro1,j));
    1923       27445 :     gunclone(p1);
    1924             :   }
    1925       14276 :   *m = roots_to_pol(allr, 0);
    1926       14276 :   return allr;
    1927             : }
    1928             : 
    1929             : static GEN
    1930       42990 : all_roots(GEN p, long bit)
    1931             : {
    1932             :   GEN lc, pd, q, roots_pol, m;
    1933       42990 :   long bit0,  bit2, i, e, h, n = degpol(p);
    1934             :   double fb;
    1935             :   pari_sp av;
    1936             : 
    1937       42990 :   pd = RgX_deflate_max(p, &h); lc = leading_coeff(pd);
    1938       42990 :   fb = fujiwara_bound(pd);
    1939       42990 :   e = (fb < 0)? 0: (long)(2 * fb);
    1940       42990 :   bit0 = bit + gexpo(pd) - gexpo(lc) + (long)log2(n/h)+1+e;
    1941       42990 :   bit2 = bit0; e = 0;
    1942       43554 :   for (av=avma,i=1;; i++,set_avma(av))
    1943             :   {
    1944       44118 :     roots_pol = vectrunc_init(n+1);
    1945       43554 :     bit2 += e + (n << i);
    1946       43554 :     q = RgX_gtofp_bit(mygprec(pd,bit2), bit2);
    1947       43554 :     q[1] = evalsigne(1)|evalvarn(0);
    1948       43554 :     m = split_complete(q,bit2,roots_pol);
    1949       43554 :     roots_pol = fix_roots(roots_pol, &m, h, bit2);
    1950       43554 :     q = mygprec_special(p,bit2); lc = leading_coeff(q);
    1951       43554 :     q[1] = evalsigne(1)|evalvarn(0);
    1952       43554 :     if (h > 1) m = gmul(m,lc);
    1953             : 
    1954       43554 :     e = gexpo(gsub(q, m)) - gexpo(lc) + (long)log2((double)n) + 1;
    1955       43554 :     if (e < -2*bit2) e = -2*bit2; /* avoid e = -pariINFINITY */
    1956       43554 :     if (e < 0)
    1957             :     {
    1958       43554 :       e = bit + a_posteriori_errors(p,roots_pol,e);
    1959       86544 :       if (e < 0) return roots_pol;
    1960             :     }
    1961         564 :     if (DEBUGLEVEL > 7)
    1962           0 :       err_printf("all_roots: restarting, i = %ld, e = %ld\n", i,e);
    1963             :   }
    1964             : }
    1965             : 
    1966             : INLINE int
    1967       12199 : isexactscalar(GEN x) { long tx = typ(x); return is_rational_t(tx); }
    1968             : 
    1969             : static int
    1970        8346 : isexactpol(GEN p)
    1971             : {
    1972        8346 :   long i,n = degpol(p);
    1973       13285 :   for (i=0; i<=n; i++)
    1974       12199 :     if (!isexactscalar(gel(p,i+2))) return 0;
    1975        1086 :   return 1;
    1976             : }
    1977             : 
    1978             : static GEN
    1979        1086 : solve_exact_pol(GEN p, long bit)
    1980             : {
    1981        1086 :   long i, j, k, m, n = degpol(p), iroots = 0;
    1982        1086 :   GEN ex, factors, v = zerovec(n);
    1983             : 
    1984        1086 :   factors = ZX_squff(Q_primpart(p), &ex);
    1985        2172 :   for (i=1; i<lg(factors); i++)
    1986             :   {
    1987        1086 :     GEN roots_fact = all_roots(gel(factors,i), bit);
    1988        1086 :     n = degpol(gel(factors,i));
    1989        1086 :     m = ex[i];
    1990        4309 :     for (j=1; j<=n; j++)
    1991        3223 :       for (k=1; k<=m; k++) v[++iroots] = roots_fact[j];
    1992             :   }
    1993        1086 :   return v;
    1994             : }
    1995             : 
    1996             : /* return the roots of p with absolute error bit */
    1997             : static GEN
    1998        8346 : roots_com(GEN q, long bit)
    1999             : {
    2000             :   GEN L, p;
    2001        8346 :   long v = RgX_valrem_inexact(q, &p);
    2002        8346 :   int ex = isexactpol(p);
    2003        8346 :   if (!ex) p = RgX_normalize1(p);
    2004        8346 :   if (lg(p) == 3)
    2005           7 :     L = cgetg(1,t_VEC); /* constant polynomial */
    2006             :   else
    2007        8339 :     L = ex? solve_exact_pol(p,bit): all_roots(p,bit);
    2008        8346 :   if (v)
    2009             :   {
    2010         168 :     GEN M, z, t = gel(q,2);
    2011             :     long i, x, y, l, n;
    2012             : 
    2013         168 :     if (isrationalzero(t)) x = -bit;
    2014             :     else
    2015             :     {
    2016          14 :       n = gexpo(t);
    2017          14 :       x = n / v; l = degpol(q);
    2018          56 :       for (i = v; i <= l; i++)
    2019             :       {
    2020          42 :         t  = gel(q,i+2);
    2021          42 :         if (isrationalzero(t)) continue;
    2022          42 :         y = (n - gexpo(t)) / i;
    2023          42 :         if (y < x) x = y;
    2024             :       }
    2025             :     }
    2026         168 :     z = real_0_bit(x); l = v + lg(L);
    2027         168 :     M = cgetg(l, t_VEC); L -= v;
    2028         168 :     for (i = 1; i <= v; i++) gel(M,i) = z;
    2029         168 :     for (     ; i <  l; i++) gel(M,i) = gel(L,i);
    2030         168 :     L = M;
    2031             :   }
    2032        8346 :   return L;
    2033             : }
    2034             : 
    2035             : static GEN
    2036      133440 : tocomplex(GEN x, long l, long bit)
    2037             : {
    2038             :   GEN y;
    2039      133440 :   if (typ(x) == t_COMPLEX)
    2040             :   {
    2041      124181 :     if (signe(gel(x,1))) return mygprecrc(x, l, -bit);
    2042       20601 :     x = gel(x,2);
    2043       20601 :     y = cgetg(3,t_COMPLEX);
    2044       20601 :     gel(y,1) = real_0_bit(-bit);
    2045       20601 :     gel(y,2) = mygprecrc(x, l, -bit);
    2046             :   }
    2047             :   else
    2048             :   {
    2049        9259 :     y = cgetg(3,t_COMPLEX);
    2050        9259 :     gel(y,1) = mygprecrc(x, l, -bit);
    2051        9259 :     gel(y,2) = real_0_bit(-bit);
    2052             :   }
    2053       29860 :   return y;
    2054             : }
    2055             : 
    2056             : /* x,y are t_COMPLEX of t_REALs or t_REAL, compare lexicographically,
    2057             :  * up to 2^-e absolute error */
    2058             : static int
    2059      265732 : cmp_complex_appr(void *E, GEN x, GEN y)
    2060             : {
    2061      265732 :   long e = (long)E;
    2062             :   GEN z, xi, yi, xr, yr;
    2063             :   long sxi, syi;
    2064      265732 :   if (typ(x) == t_COMPLEX) { xr = gel(x,1); xi = gel(x,2); sxi = signe(xi); }
    2065       84225 :   else { xr = x; xi = NULL; sxi = 0; }
    2066      265732 :   if (typ(y) == t_COMPLEX) { yr = gel(y,1); yi = gel(y,2); syi = signe(yi); }
    2067       76047 :   else { yr = y; yi = NULL; syi = 0; }
    2068             :   /* Compare absolute values of imaginary parts */
    2069      265732 :   if (!sxi)
    2070             :   {
    2071       93722 :     if (syi && expo(yi) >= e) return -1;
    2072             :     /* |Im x| ~ |Im y| ~ 0 */
    2073             :   }
    2074      172010 :   else if (!syi)
    2075             :   {
    2076        4184 :     if (sxi && expo(xi) >= e) return 1;
    2077             :     /* |Im x| ~ |Im y| ~ 0 */
    2078             :   }
    2079             :   else
    2080             :   {
    2081             :     long sz;
    2082      167826 :     z = addrr_sign(xi, 1, yi, -1);
    2083      167826 :     sz = signe(z);
    2084      167826 :     if (sz && expo(z) >= e) return (int)sz;
    2085             :   }
    2086             :   /* |Im x| ~ |Im y|, sort according to real parts */
    2087      154306 :   z = subrr(xr, yr);
    2088      154306 :   if (expo(z) >= e) return (int)signe(z);
    2089             :   /* Re x ~ Re y. Place negative absolute value before positive */
    2090       52113 :   return (int) (sxi - syi);
    2091             : }
    2092             : 
    2093             : static GEN
    2094       43025 : clean_roots(GEN L, long l, long bit, long clean)
    2095             : {
    2096       43025 :   long i, n = lg(L), ex = 5 - bit;
    2097       43025 :   GEN res = cgetg(n,t_COL);
    2098      215526 :   for (i=1; i<n; i++)
    2099             :   {
    2100      172501 :     GEN c = gel(L,i);
    2101      172501 :     if (clean && isrealappr(c,ex))
    2102             :     {
    2103       39061 :       if (typ(c) == t_COMPLEX) c = gel(c,1);
    2104       39061 :       c = mygprecrc(c, l, -bit);
    2105             :     }
    2106             :     else
    2107      133440 :       c = tocomplex(c, l, bit);
    2108      172501 :     gel(res,i) = c;
    2109             :   }
    2110       43025 :   gen_sort_inplace(res, (void*)ex, &cmp_complex_appr, NULL);
    2111       43025 :   return res;
    2112             : }
    2113             : 
    2114             : /* the vector of roots of p, with absolute error 2^(- prec2nbits(l)) */
    2115             : static GEN
    2116        8374 : roots_aux(GEN p, long l, long clean)
    2117             : {
    2118        8374 :   pari_sp av = avma;
    2119             :   long bit;
    2120             :   GEN L;
    2121             : 
    2122        8374 :   if (typ(p) != t_POL)
    2123             :   {
    2124          21 :     if (gequal0(p)) pari_err_ROOTS0("roots");
    2125          14 :     if (!isvalidcoeff(p)) pari_err_TYPE("roots",p);
    2126           7 :     return cgetg(1,t_COL); /* constant polynomial */
    2127             :   }
    2128        8353 :   if (!signe(p)) pari_err_ROOTS0("roots");
    2129        8353 :   checkvalidpol(p,"roots");
    2130        8346 :   if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
    2131        8346 :   if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
    2132        8346 :   bit = prec2nbits(l);
    2133        8346 :   L = roots_com(p, bit);
    2134        8346 :   return gerepileupto(av, clean_roots(L, l, bit, clean));
    2135             : }
    2136             : GEN
    2137        8157 : roots(GEN p, long l) { return roots_aux(p,l, 0); }
    2138             : /* clean up roots. If root is real replace it by its real part */
    2139             : GEN
    2140         217 : cleanroots(GEN p, long l) { return roots_aux(p,l, 1); }
    2141             : 
    2142             : /* private variant of conjvec. Allow non rational coefficients, shallow
    2143             :  * function. */
    2144             : GEN
    2145          84 : polmod_to_embed(GEN x, long prec)
    2146             : {
    2147          84 :   GEN v, T = gel(x,1), A = gel(x,2);
    2148             :   long i, l;
    2149          84 :   if (typ(A) != t_POL || varn(A) != varn(T))
    2150             :   {
    2151           7 :     checkvalidpol(T,"polmod_to_embed");
    2152           7 :     return const_col(degpol(T), A);
    2153             :   }
    2154          77 :   v = cleanroots(T,prec); l = lg(v);
    2155          77 :   for (i=1; i<l; i++) gel(v,i) = poleval(A,gel(v,i));
    2156          77 :   return v;
    2157             : }
    2158             : 
    2159             : GEN
    2160       34679 : QX_complex_roots(GEN p, long l)
    2161             : {
    2162       34679 :   pari_sp av = avma;
    2163             :   long bit, v;
    2164             :   GEN L;
    2165             : 
    2166       34679 :   if (!signe(p)) pari_err_ROOTS0("QX_complex_roots");
    2167       34679 :   if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
    2168       34679 :   if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
    2169       34679 :   bit = prec2nbits(l);
    2170       34679 :   v = RgX_valrem(p, &p);
    2171       34679 :   L = lg(p) > 3? all_roots(Q_primpart(p), bit): cgetg(1,t_COL);
    2172       34679 :   if (v) L = shallowconcat(const_vec(v, real_0_bit(-bit)), L);
    2173       34679 :   return gerepileupto(av, clean_roots(L, l, bit, 1));
    2174             : }
    2175             : 
    2176             : /********************************************************************/
    2177             : /**                                                                **/
    2178             : /**                REAL ROOTS OF INTEGER POLYNOMIAL                **/
    2179             : /**                                                                **/
    2180             : /********************************************************************/
    2181             : 
    2182             : /* Count sign changes in the coefficients of (x+1)^deg(P)*P(1/(x+1))
    2183             :  * The inversion is implicit (we take coefficients backwards). Roots of P
    2184             :  * at 0 and 1 (mapped to oo and 0) are ignored here and must be dealt with
    2185             :  * by the caller. Set root1 if P(1) = 0 */
    2186             : static long
    2187      502557 : X2XP1(GEN P, long deg, int *root1, GEN *Premapped)
    2188             : {
    2189      502557 :   const pari_sp av = avma;
    2190      502557 :   GEN v = shallowcopy(P);
    2191             :   long i, j, vlim, nb, s;
    2192             : 
    2193      502557 :   for (i = 0, vlim = deg+2;;)
    2194             :   {
    2195      502711 :     for (j = 2; j < vlim; j++) gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
    2196      502634 :     s = -signe(gel(v, vlim));
    2197      502634 :     vlim--; i++; if (s) break;
    2198             :   }
    2199      502557 :   if (vlim == deg+1) *root1 = 0;
    2200             :   else
    2201             :   {
    2202          77 :     *root1 = 1;
    2203          77 :     if (Premapped) setlg(v, vlim + 2);
    2204             :   }
    2205             : 
    2206      502557 :   nb = 0;
    2207     1869721 :   for (; i < deg; i++)
    2208             :   {
    2209     1755022 :     long s2 = -signe(gel(v, 2));
    2210     1755022 :     int flag = (s2 == s);
    2211    21392611 :     for (j = 2; j < vlim; j++)
    2212             :     {
    2213    19637589 :       gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
    2214    19637589 :       if (flag) flag = (s2 != signe(gel(v, j+1)));
    2215             :     }
    2216     1755022 :     if (s == signe(gel(v, vlim)))
    2217             :     {
    2218      475411 :       if (++nb >= 2) return gc_long(av,2);
    2219      293968 :       s = -s;
    2220             :     }
    2221             :     /* if flag is set there will be no further sign changes */
    2222     1573579 :     if (flag && (!Premapped || !nb)) goto END;
    2223     1367164 :     vlim--;
    2224     1367164 :     if (gc_needed(av, 3))
    2225             :     {
    2226           0 :       if (DEBUGMEM>1) pari_warn(warnmem, "X2XP1, i = %ld/%ld", i, deg-1);
    2227           0 :       if (!Premapped) setlg(v, vlim + 2);
    2228           0 :       v = gerepilecopy(av, v);
    2229             :     }
    2230             :   }
    2231      114699 :   if (vlim >= 2 && s == signe(gel(v, vlim))) nb++;
    2232             : END:
    2233      321114 :   if (Premapped && nb == 1) *Premapped = v; else set_avma(av);
    2234      321114 :   return nb;
    2235             : }
    2236             : 
    2237             : static long
    2238           0 : _intervalcmp(GEN x, GEN y)
    2239             : {
    2240           0 :   if (typ(x) == t_VEC) x = gel(x, 1);
    2241           0 :   if (typ(y) == t_VEC) y = gel(y, 1);
    2242           0 :   return gcmp(x, y);
    2243             : }
    2244             : 
    2245             : static GEN
    2246     5248713 : _gen_nored(void *E, GEN x) { (void)E; return x; }
    2247             : static GEN
    2248    17602601 : _mp_add(void *E, GEN x, GEN y) { (void)E; return mpadd(x, y); }
    2249             : static GEN
    2250           0 : _mp_sub(void *E, GEN x, GEN y) { (void)E; return mpsub(x, y); }
    2251             : static GEN
    2252     3188839 : _mp_mul(void *E, GEN x, GEN y) { (void)E; return mpmul(x, y); }
    2253             : static GEN
    2254     4097747 : _mp_sqr(void *E, GEN x) { (void)E; return mpsqr(x); }
    2255             : static GEN
    2256     6735219 : _gen_one(void *E) { (void)E; return gen_1; }
    2257             : static GEN
    2258       34981 : _gen_zero(void *E) { (void)E; return gen_0; }
    2259             : 
    2260             : static struct bb_algebra mp_algebra = { _gen_nored, _mp_add, _mp_sub,
    2261             :                          _mp_mul, _mp_sqr, _gen_one, _gen_zero };
    2262             : 
    2263             : static GEN
    2264    22110750 : _mp_cmul(void *E, GEN P, long a, GEN x) {(void)E; return mpmul(gel(P,a+2), x);}
    2265             : 
    2266             : /* Split the polynom P in two parts, whose coeffs have constant sign:
    2267             :  * P(X) = X^D*Pp + Pm. Also compute the two parts of the derivative of P,
    2268             :  * Pprimem = Pm', Pprimep = X*Pp'+ D*Pp => P' = X^(D-1)*Pprimep + Pprimem;
    2269             :  * Pprimep[i] = (i+D) Pp[i]. Return D */
    2270             : static long
    2271       69684 : split_pols(GEN P, GEN *pPp, GEN *pPm, GEN *pPprimep, GEN *pPprimem)
    2272             : {
    2273       69684 :   long i, D, dP = degpol(P), s0 = signe(gel(P,2));
    2274             :   GEN Pp, Pm, Pprimep, Pprimem;
    2275      292098 :   for(i=1; i <= dP; i++)
    2276      292098 :     if (signe(gel(P, i+2)) == -s0) break;
    2277       69684 :   D = i;
    2278       69684 :   Pm = cgetg(D + 2, t_POL);
    2279       69684 :   Pprimem = cgetg(D + 1, t_POL);
    2280       69684 :   Pp = cgetg(dP-D + 3, t_POL);
    2281       69684 :   Pprimep = cgetg(dP-D + 3, t_POL);
    2282       69684 :   Pm[1] = Pp[1] = Pprimem[1] = Pprimep[1] = P[1];
    2283      361782 :   for(i=0; i < D; i++)
    2284             :   {
    2285      292098 :     GEN c = gel(P, i+2);
    2286      292098 :     gel(Pm, i+2) = c;
    2287      292098 :     if (i) gel(Pprimem, i+1) = mului(i, c);
    2288             :   }
    2289      393248 :   for(; i <= dP; i++)
    2290             :   {
    2291      323564 :     GEN c = gel(P, i+2);
    2292      323564 :     gel(Pp, i+2-D) = c;
    2293      323564 :     gel(Pprimep, i+2-D) = mului(i, c);
    2294             :   }
    2295       69684 :   *pPm = normalizepol_lg(Pm, D+2);
    2296       69684 :   *pPprimem = normalizepol_lg(Pprimem, D+1);
    2297       69684 :   *pPp = normalizepol_lg(Pp, dP-D+3);
    2298       69684 :   *pPprimep = normalizepol_lg(Pprimep, dP-D+3);
    2299       69684 :   return dP - degpol(*pPp);
    2300             : }
    2301             : 
    2302             : static GEN
    2303     2271565 : bkeval_single_power(long d, GEN V)
    2304             : {
    2305     2271565 :   long mp = lg(V) - 2;
    2306     2271565 :   if (d > mp) return gmul(gpowgs(gel(V, mp+1), d/mp), gel(V, (d%mp)+1));
    2307     2271565 :   return gel(V, d+1);
    2308             : }
    2309             : 
    2310             : static GEN
    2311     2271565 : splitpoleval(GEN Pp, GEN Pm, GEN pows, long D, long bitprec)
    2312             : {
    2313     2271565 :   GEN vp = gen_bkeval_powers(Pp, degpol(Pp), pows, NULL, &mp_algebra, _mp_cmul);
    2314     2271565 :   GEN vm = gen_bkeval_powers(Pm, degpol(Pm), pows, NULL, &mp_algebra, _mp_cmul);
    2315     2271565 :   GEN xa = bkeval_single_power(D, pows);
    2316             :   GEN r;
    2317     2271565 :   if (!signe(vp)) return vm;
    2318     2271565 :   vp = gmul(vp, xa);
    2319     2271565 :   r = gadd(vp, vm);
    2320     2271565 :   if (gexpo(vp) - (signe(r)? gexpo(r): 0) > prec2nbits(realprec(vp)) - bitprec)
    2321       52209 :     return NULL;
    2322     2219356 :   return r;
    2323             : }
    2324             : 
    2325             : /* optimized Cauchy bound for P = X^D*Pp + Pm, D > deg(Pm) */
    2326             : static GEN
    2327       69684 : splitcauchy(GEN Pp, GEN Pm, long prec)
    2328             : {
    2329       69684 :   GEN S = gel(Pp,2), A = gel(Pm,2);
    2330       69684 :   long i, lPm = lg(Pm), lPp = lg(Pp);
    2331       69684 :   for (i=3; i < lPm; i++) { GEN c = gel(Pm,i); if (abscmpii(A, c) < 0) A = c; }
    2332       69684 :   for (i=3; i < lPp; i++) S = addii(S, gel(Pp, i));
    2333       69684 :   return subsr(1, rdivii(A, S, prec)); /* 1 + |Pm|_oo / |Pp|_1 */
    2334             : }
    2335             : 
    2336             : /* Newton for polynom P, P(0)!=0, with unique sign change => one root in ]0,oo[
    2337             :  * P' has also at most one zero there */
    2338             : static GEN
    2339       69684 : polsolve(GEN P, long bitprec)
    2340             : {
    2341       69684 :   pari_sp av = avma, av2;
    2342             :   GEN Pp, Pm, Pprimep, Pprimem, Pprime, Pprime2, ra, rb, rc, Pc;
    2343       69684 :   long deg = degpol(P);
    2344       69684 :   long expoold = LONG_MAX, cprec = DEFAULTPREC, prec = nbits2prec(bitprec);
    2345             :   long iter, D, rt, s0, bitaddprec, addprec;
    2346       69684 :   if (deg == 1)
    2347           0 :     return gerepileuptoleaf(av, rdivii(negi(gel(P,2)), gel(P,3), prec));
    2348       69684 :   Pprime = ZX_deriv(P);
    2349       69684 :   Pprime2 = ZX_deriv(Pprime);
    2350       69684 :   bitaddprec = 1 + 2*expu(deg); addprec = nbits2prec(bitaddprec);
    2351       69684 :   D = split_pols(P, &Pp, &Pm, &Pprimep, &Pprimem); /* P = X^D*Pp + Pm */
    2352       69684 :   s0 = signe(gel(P, 2));
    2353       69684 :   rt = maxss(D, brent_kung_optpow(maxss(degpol(Pp), degpol(Pm)), 2, 1));
    2354       69684 :   rb = splitcauchy(Pp, Pm, DEFAULTPREC);
    2355             :   for(;;)
    2356           0 :   {
    2357       69684 :     GEN pows = gen_powers(rb, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
    2358       69684 :     Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec);
    2359       69684 :     if (!Pc) { cprec++; rb = rtor(rb, cprec); continue; }
    2360       69684 :     if (signe(Pc) != s0) break;
    2361           0 :     shiftr_inplace(rb,1);
    2362             :   }
    2363       69684 :   ra = NULL;
    2364       69684 :   iter = 0;
    2365             :   for(;;)
    2366     1117651 :   {
    2367             :     GEN wdth;
    2368     1187335 :     iter++;
    2369     1187335 :     if (ra)
    2370      345631 :       rc = shiftr(addrr(ra, rb), -1);
    2371             :     else
    2372      841704 :       rc = shiftr(rb, -1);
    2373             :     do
    2374           0 :     {
    2375     1187335 :       GEN pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
    2376     1187335 :       Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec+2);
    2377     1187335 :       if (Pc) break;
    2378           0 :       cprec++; rc = rtor(rc, cprec);
    2379             :     } while (1);
    2380     1187335 :     if (signe(Pc) == s0)
    2381      235272 :       ra = rc;
    2382             :     else
    2383      952063 :       rb = rc;
    2384     1187335 :     if (!ra) continue;
    2385      415315 :     wdth = subrr(rb, ra);
    2386      415315 :     if (!(iter % 8))
    2387             :     {
    2388       69964 :       GEN m1 = poleval(Pprime, ra), M2;
    2389       69964 :       if (signe(m1) == s0) continue;
    2390       69103 :       M2 = poleval(Pprime2, rb);
    2391       69103 :       if (abscmprr(gmul(M2, wdth), shiftr(m1, 1)) > 0) continue;
    2392       66548 :       break;
    2393             :     }
    2394      345351 :     else if (gexpo(wdth) <= -bitprec)
    2395        3136 :       break;
    2396             :   }
    2397       69684 :   rc = rb;
    2398       69684 :   av2 = avma;
    2399      437589 :   for(;; rc = gerepileuptoleaf(av2, rc))
    2400      437589 :   {
    2401             :     long exponew;
    2402      507273 :     GEN Ppc, dist, rcold = rc;
    2403      507273 :     GEN pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
    2404      507273 :     Ppc = splitpoleval(Pprimep, Pprimem, pows, D-1, bitaddprec+4);
    2405      507273 :     if (Ppc)
    2406      507273 :       Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec+4);
    2407      507273 :     if (!Ppc || !Pc)
    2408             :     {
    2409       52209 :       if (cprec >= prec+addprec)
    2410        3803 :         cprec += EXTRAPRECWORD;
    2411             :       else
    2412       48406 :         cprec = minss(2*cprec, prec+addprec);
    2413       52209 :       rc = rtor(rc, cprec); continue; /* backtrack one step */
    2414             :     }
    2415      455064 :     dist = typ(Ppc) == t_REAL? divrr(Pc, Ppc): divri(Pc, Ppc);
    2416      455064 :     rc = subrr(rc, dist);
    2417      455064 :     if (cmprr(ra, rc) > 0 || cmprr(rb, rc) < 0)
    2418             :     {
    2419           0 :       if (cprec >= prec+addprec) break;
    2420           0 :       cprec = minss(2*cprec, prec+addprec);
    2421           0 :       rc = rtor(rcold, cprec); continue; /* backtrack one step */
    2422             :     }
    2423      455064 :     if (expoold == LONG_MAX) { expoold = expo(dist); continue; }
    2424      381679 :     exponew = expo(dist);
    2425      381679 :     if (exponew < -bitprec - 1)
    2426             :     {
    2427      146372 :       if (cprec >= prec+addprec) break;
    2428       76688 :       cprec = minss(2*cprec, prec+addprec);
    2429       76688 :       rc = rtor(rc, cprec); continue;
    2430             :     }
    2431      235307 :     if (exponew > expoold - 2)
    2432             :     {
    2433        3701 :       if (cprec >= prec+addprec) break;
    2434        3701 :       expoold = LONG_MAX;
    2435        3701 :       cprec = minss(2*cprec, prec+addprec);
    2436        3701 :       rc = rtor(rc, cprec); continue;
    2437             :     }
    2438      231606 :     expoold = exponew;
    2439             :   }
    2440       69684 :   return gerepileuptoleaf(av, rtor(rc, prec));
    2441             : }
    2442             : 
    2443             : static GEN
    2444       68329 : usp(GEN Q0, long deg, long flag, long bitprec)
    2445             : {
    2446       68329 :   const pari_sp av = avma;
    2447       68329 :   GEN Qremapped, Q, c, Lc, Lk, sol = zerocol(deg);
    2448       68329 :   GEN *pQremapped = flag == 1? &Qremapped: NULL;
    2449       68329 :   const long prec = nbits2prec(bitprec);
    2450       68329 :   long listsize = 64, nbr = 0, nb_todo, ind, deg0, indf, i, k, nb;
    2451             : 
    2452       68329 :   deg0 = deg;
    2453       68329 :   Lc = zerovec(listsize);
    2454       68329 :   Lk = cgetg(listsize+1, t_VECSMALL);
    2455       68329 :   c = gen_0;
    2456       68329 :   k = Lk[1] = 0;
    2457       68329 :   ind = 1; indf = 2;
    2458       68329 :   Q = leafcopy(Q0);
    2459             : 
    2460       68329 :   nb_todo = 1;
    2461      639215 :   while (nb_todo)
    2462             :   {
    2463      502557 :     GEN nc = gel(Lc, ind);
    2464             :     pari_sp av2;
    2465             :     int root1;
    2466      502557 :     if (Lk[ind] == k + 1)
    2467             :     {
    2468      147186 :       deg0 = deg;
    2469      147186 :       setlg(Q0, deg + 3);
    2470      147186 :       Q0 = ZX_rescale2n(Q0, 1);
    2471      147186 :       Q = Q_primpart(Q0);
    2472      147186 :       c = gen_0;
    2473             :     }
    2474      502557 :     if (!equalii(nc, c)) Q = ZX_translate(Q, subii(nc, c));
    2475      502557 :     k = Lk[ind];
    2476      502557 :     c = nc;
    2477      502557 :     ind++;
    2478      502557 :     nb_todo--;
    2479             : 
    2480      502557 :     if (equalii(gel(Q, 2), gen_0))
    2481             :     { /* Q(0) = 0 */
    2482         168 :       GEN s = gmul2n(c, -k);
    2483             :       long j;
    2484         210 :       for (j = 1; j <= nbr; j++)
    2485         126 :         if (gequal(gel(sol, j), s)) break;
    2486         168 :       if (j > nbr) gel(sol, ++nbr) = s;
    2487             : 
    2488         168 :       deg0--;
    2489         168 :       for (j = 2; j <= deg0 + 2; j++) gel(Q, j) = gel(Q, j+1);
    2490         168 :       setlg(Q, j);
    2491             :     }
    2492             : 
    2493      502557 :     av2 = avma;
    2494      502557 :     nb = X2XP1(Q, deg0, &root1, pQremapped);
    2495             : 
    2496      502557 :     if      (nb == 0) set_avma(av2);/* no root in this open interval */
    2497      302407 :     else if (nb == 1) /* exactly one root */
    2498             :     {
    2499       85293 :       GEN s = gen_0;
    2500       85293 :       if (flag == 0)
    2501           0 :         s = mkvec2(gmul2n(c,-k), gmul2n(addiu(c,1),-k));
    2502       85293 :       else if (flag == 1) /* Caveat: Qremapped is the reciprocal polynomial */
    2503             :       {
    2504       69684 :         s = polsolve(*pQremapped, bitprec+1);
    2505       69684 :         s = addir(c, divrr(s, addsr(1, s)));
    2506       69684 :         shiftr_inplace(s, -k);
    2507       69684 :         if (realprec(s) != prec) s = rtor(s, prec);
    2508             :       }
    2509       85293 :       gel(sol, ++nbr) = gerepileupto(av2, s);
    2510             :     }
    2511             :     else
    2512             :     { /* unknown, add two nodes to refine */
    2513      217114 :       if (indf + 2 > listsize)
    2514             :       {
    2515         413 :         if (ind>1)
    2516             :         {
    2517        2940 :           for (i = ind; i < indf; i++)
    2518             :           {
    2519        2527 :             gel(Lc, i-ind+1) = gel(Lc, i);
    2520        2527 :             Lk[i-ind+1] = Lk[i];
    2521             :           }
    2522         413 :           indf -= ind-1;
    2523         413 :           ind = 1;
    2524             :         }
    2525         413 :         if (indf + 2 > listsize)
    2526             :         {
    2527           0 :           listsize *= 2;
    2528           0 :           Lc = vec_lengthen(Lc, listsize);
    2529           0 :           Lk = vecsmall_lengthen(Lk, listsize);
    2530             :         }
    2531         413 :         for (i = indf; i <= listsize; i++) gel(Lc, i) = gen_0;
    2532             :       }
    2533      217114 :       nc = shifti(c, 1);
    2534      217114 :       gel(Lc, indf) = nc;
    2535      217114 :       gel(Lc, indf + 1) = addiu(nc, 1);
    2536      217114 :       Lk[indf] = Lk[indf + 1] = k + 1;
    2537      217114 :       indf += 2;
    2538      217114 :       nb_todo += 2;
    2539             :     }
    2540      502557 :     if (root1)
    2541             :     { /* Q(1) = 0 */
    2542          77 :       GEN s = gmul2n(addiu(c,1), -k);
    2543             :       long j;
    2544         119 :       for (j = 1; j <= nbr; j++)
    2545          49 :         if (gequal(gel(sol, j), s)) break;
    2546          77 :       if (j > nbr) gel(sol, ++nbr) = s;
    2547             :     }
    2548             : 
    2549      502557 :     if (gc_needed(av, 2))
    2550             :     {
    2551           0 :       gerepileall(av, 6, &Q0, &Q, &c, &Lc, &Lk, &sol);
    2552           0 :       if (DEBUGMEM > 1) pari_warn(warnmem, "ZX_Uspensky", avma);
    2553             :     }
    2554             :   }
    2555             : 
    2556       68329 :   setlg(sol, nbr+1);
    2557       68329 :   return gerepilecopy(av, sol);
    2558             : }
    2559             : 
    2560             : static GEN
    2561        6118 : ZX_Uspensky_cst_pol(long nbz, long flag, long bitprec)
    2562             : {
    2563        6118 :   switch(flag)
    2564             :   {
    2565           0 :     case 0:  return zerocol(nbz);
    2566        1078 :     case 1:  retconst_col(nbz, real_0_bit(-bitprec));
    2567        5040 :     default: return utoi(nbz);
    2568             :   }
    2569             : }
    2570             : 
    2571             : /* ZX_Uspensky(P, [a,a], flag) */
    2572             : static GEN
    2573          28 : ZX_Uspensky_equal(GEN P, GEN a, long flag)
    2574             : {
    2575          28 :   if (typ(a) != t_INFINITY && gequal0(poleval(P, a)))
    2576          21 :     return flag <= 1 ? mkcol(a): gen_1;
    2577             :   else
    2578           7 :     return flag <= 1 ? cgetg(1, t_COL) : gen_0;
    2579             : }
    2580             : 
    2581             : GEN
    2582       73037 : ZX_Uspensky(GEN P, GEN ab, long flag, long bitprec)
    2583             : {
    2584       73037 :   pari_sp av = avma;
    2585       73037 :   GEN a, b, sol = NULL, Pcur;
    2586             :   double fb;
    2587             :   long nbz, deg;
    2588             : 
    2589       73037 :   deg = degpol(P);
    2590       73037 :   if (deg == 0) return flag <= 1 ? cgetg(1, t_COL) : gen_0;
    2591       73037 :   if (ab)
    2592             :   {
    2593       57049 :     if (typ(ab) == t_VEC)
    2594             :     {
    2595       47513 :       if (lg(ab) != 3) pari_err_DIM("ZX_Uspensky");
    2596       47513 :       a = gel(ab, 1);
    2597       47513 :       b = gel(ab, 2);
    2598             :     }
    2599             :     else
    2600             :     {
    2601        9536 :       a = ab;
    2602        9536 :       b = mkoo();
    2603             :     }
    2604             :   }
    2605             :   else
    2606             :   {
    2607       15988 :     a = mkmoo();
    2608       15988 :     b = mkoo();
    2609             :   }
    2610       73037 :   switch (gcmp(a, b))
    2611             :   {
    2612           7 :     case 1: set_avma(av); return flag <= 1 ? cgetg(1, t_COL) : gen_0;
    2613          21 :     case 0: return gerepilecopy(av, ZX_Uspensky_equal(P, a, flag));
    2614             :   }
    2615       73009 :   nbz = ZX_valrem(P, &Pcur);
    2616       73009 :   deg -= nbz;
    2617       73009 :   if (!nbz) Pcur = P;
    2618       73009 :   if (nbz && (gsigne(a) > 0 || gsigne(b) < 0)) nbz = 0;
    2619       73009 :   if (deg == 0) { set_avma(av); return ZX_Uspensky_cst_pol(nbz, flag, bitprec); }
    2620       71574 :   if (deg == 1)
    2621             :   {
    2622        9480 :     sol = gdiv(gneg(gel(Pcur, 2)), pollead(Pcur, -1));
    2623        9480 :     if (gcmp(a, sol) > 0 || gcmp(sol, b) > 0)
    2624             :     {
    2625        4683 :       set_avma(av);
    2626        4683 :       return ZX_Uspensky_cst_pol(nbz, flag, bitprec);
    2627             :     }
    2628        4797 :     if (flag >= 2) { set_avma(av); return utoi(nbz+1); }
    2629        2235 :     sol = gconcat(zerocol(nbz), mkcol(sol));
    2630        2235 :     if (flag == 1) sol = RgC_gtofp(sol, nbits2prec(bitprec));
    2631        2235 :     return gerepilecopy(av, sol);
    2632             :   }
    2633       62094 :   switch(flag)
    2634             :   {
    2635             :     case 0:
    2636           0 :       sol = zerocol(nbz);
    2637           0 :       break;
    2638             :     case 1:
    2639       48230 :       sol = const_col(nbz, real_0_bit(-bitprec));
    2640       48230 :       break;
    2641             :     /* case 2: nothing */
    2642             :   }
    2643             : 
    2644       62094 :   if (typ(a) == t_INFINITY && typ(b) != t_INFINITY && gsigne(b))
    2645             :   {
    2646          35 :     fb = fujiwara_bound_real(Pcur, -1);
    2647          35 :     if (fb > -pariINFINITY) a = negi(int2n((long)ceil(fb))); else a = gen_0;
    2648             :   }
    2649       62094 :   if (typ(b) == t_INFINITY && typ(a) != t_INFINITY && gsigne(a))
    2650             :   {
    2651          21 :     fb = fujiwara_bound_real(Pcur, 1);
    2652          21 :     if (fb > -pariINFINITY) b = int2n((long)ceil(fb)); else b = gen_0;
    2653             :   }
    2654             : 
    2655       62094 :   if (typ(a) != t_INFINITY && typ(b) != t_INFINITY)
    2656             :   {
    2657             :     GEN den, diff, unscaledres, co, Pdiv, ascaled;
    2658             :     pari_sp av1;
    2659             :     long i;
    2660        7287 :     if (gequal(a,b)) /* can occur if one of a,b was initially a t_INFINITY */
    2661           7 :       return gerepilecopy(av, ZX_Uspensky_equal(P, a, flag));
    2662        7280 :     den = lcmii(Q_denom(a), Q_denom(b));
    2663        7280 :     if (!is_pm1(den))
    2664             :     {
    2665          42 :       Pcur = ZX_rescale(Pcur, den);
    2666          42 :       ascaled = gmul(a, den);
    2667             :     }
    2668             :     else
    2669             :     {
    2670        7238 :       den = NULL;
    2671        7238 :       ascaled = a;
    2672             :     }
    2673        7280 :     diff = subii(den ? gmul(b,den) : b, ascaled);
    2674        7280 :     Pcur = ZX_unscale(ZX_translate(Pcur, ascaled), diff);
    2675        7280 :     av1 = avma;
    2676        7280 :     Pdiv = cgetg(deg+2, t_POL);
    2677        7280 :     Pdiv[1] = Pcur[1];
    2678        7280 :     co = gel(Pcur, deg+2);
    2679       56959 :     for (i = deg; --i >= 0; )
    2680             :     {
    2681       42399 :       gel(Pdiv, i+2) = co;
    2682       42399 :       co = addii(co, gel(Pcur, i+2));
    2683             :     }
    2684        7280 :     if (!signe(co))
    2685             :     {
    2686          77 :       Pcur = Pdiv;
    2687          77 :       deg--;
    2688          77 :       if (flag <= 1)
    2689          21 :         sol = gconcat(sol, b);
    2690             :       else
    2691          56 :         nbz++;
    2692             :     }
    2693             :     else
    2694        7203 :       set_avma(av1);
    2695        7280 :     unscaledres = usp(Pcur, deg, flag, bitprec);
    2696        7280 :     if (flag <= 1)
    2697             :     {
    2698       19719 :       for (i = 1; i < lg(unscaledres); i++)
    2699             :       {
    2700       12565 :         GEN z = gmul(diff, gel(unscaledres, i));
    2701       12565 :         if (typ(z) == t_VEC)
    2702             :         {
    2703           0 :           gel(z, 1) = gadd(ascaled, gel(z, 1));
    2704           0 :           gel(z, 2) = gadd(ascaled, gel(z, 2));
    2705             :         }
    2706             :         else
    2707       12565 :           z = gadd(ascaled, z);
    2708       12565 :         if (den) z = gdiv(z, den);
    2709       12565 :         gel(unscaledres, i) = z;
    2710             :       }
    2711        7154 :       sol = gconcat(sol, unscaledres);
    2712             :     }
    2713             :     else
    2714         126 :       nbz += lg(unscaledres) - 1;
    2715             :   }
    2716       62087 :   if (typ(b) == t_INFINITY && (fb=fujiwara_bound_real(Pcur, 1)) > -pariINFINITY)
    2717             :   {
    2718             :     GEN Pcurp, unscaledres;
    2719       53189 :     long bp = (long)ceil(fb);
    2720       53189 :     if (bp < 0) bp = 0;
    2721       53189 :     Pcurp = ZX_unscale2n(Pcur, bp);
    2722       53189 :     unscaledres = usp(Pcurp, deg, flag, bitprec);
    2723       53189 :     if (flag <= 1)
    2724       41069 :       sol = shallowconcat(sol, gmul2n(unscaledres, bp));
    2725             :     else
    2726       12120 :       nbz += lg(unscaledres)-1;
    2727             :   }
    2728       62087 :   if (typ(a) == t_INFINITY && (fb=fujiwara_bound_real(Pcur,-1)) > -pariINFINITY)
    2729             :   {
    2730             :     GEN Pcurm, unscaledres;
    2731        7860 :     long i, bm = (long)ceil(fb);
    2732        7860 :     if (bm < 0) bm = 0;
    2733        7860 :     Pcurm = ZX_unscale2n(ZX_z_unscale(Pcur, -1), bm);
    2734        7860 :     unscaledres = usp(Pcurm, deg, flag, bitprec);
    2735        7860 :     if (flag <= 1)
    2736             :     {
    2737        8554 :       for (i = 1; i < lg(unscaledres); i++)
    2738             :       {
    2739        5355 :         GEN z = gneg(gmul2n(gel(unscaledres, i), bm));
    2740        5355 :         if (typ(z) == t_VEC) swap(gel(z, 1), gel(z, 2));
    2741        5355 :         gel(unscaledres, i) = z;
    2742             :       }
    2743        3199 :       sol = shallowconcat(unscaledres, sol);
    2744             :     }
    2745             :     else
    2746        4661 :       nbz += lg(unscaledres)-1;
    2747             :   }
    2748       62087 :   if (flag >= 2) return utoi(nbz);
    2749       48230 :   if (flag)
    2750       48230 :     sol = sort(sol);
    2751             :   else
    2752           0 :     sol = gen_sort(sol, (void *)_intervalcmp, cmp_nodata);
    2753       48230 :   return gerepileupto(av, sol);
    2754             : }
    2755             : 
    2756             : /* x a scalar */
    2757             : static GEN
    2758          35 : rootsdeg0(GEN x)
    2759             : {
    2760          35 :   if (!is_real_t(typ(x))) pari_err_TYPE("realroots",x);
    2761          28 :   if (gequal0(x)) pari_err_ROOTS0("realroots");
    2762          14 :   return cgetg(1,t_COL); /* constant polynomial */
    2763             : }
    2764             : static void
    2765       79150 : checkbound(GEN a)
    2766             : {
    2767       79150 :   switch(typ(a))
    2768             :   {
    2769       79150 :     case t_INT: case t_FRAC: case t_INFINITY: break;
    2770           0 :     default: pari_err_TYPE("polrealroots", a);
    2771             :   }
    2772       79150 : }
    2773             : static GEN
    2774       47137 : check_ab(GEN ab)
    2775             : {
    2776             :   GEN a, b;
    2777       47137 :   if (!ab) return NULL;
    2778       39575 :   if (typ(ab) != t_VEC || lg(ab) != 3) pari_err_TYPE("polrootsreal",ab);
    2779       39575 :   a = gel(ab,1); checkbound(a);
    2780       39575 :   b = gel(ab,2); checkbound(b);
    2781       39848 :   if (typ(a) == t_INFINITY && inf_get_sign(a) < 0 &&
    2782         511 :       typ(b) == t_INFINITY && inf_get_sign(b) > 0) ab = NULL;
    2783       39575 :   return ab;
    2784             : }
    2785             : GEN
    2786       44690 : realroots(GEN P, GEN ab, long prec)
    2787             : {
    2788       44690 :   pari_sp av = avma;
    2789       44690 :   long nrr = 0;
    2790       44690 :   GEN sol = NULL, fa, ex;
    2791             :   long i, j, v;
    2792             : 
    2793       44690 :   ab = check_ab(ab);
    2794       44690 :   if (typ(P) != t_POL) return rootsdeg0(P);
    2795       44669 :   switch(degpol(P))
    2796             :   {
    2797           7 :     case -1: return rootsdeg0(gen_0);
    2798           7 :     case 0: return rootsdeg0(gel(P,2));
    2799             :   }
    2800       44655 :   if (!RgX_is_ZX(P)) P = RgX_rescale_to_int(P);
    2801       44655 :   P = Q_primpart(P);
    2802       44655 :   v = ZX_valrem(P,&P);
    2803       44655 :   if (v && (!ab || (gsigne(gel(ab,1)) <= 0 && gsigne(gel(ab,2)) >= 0)))
    2804        1351 :     sol = const_col(v, real_0(prec));
    2805       44655 :   fa = ZX_squff(P, &ex);
    2806       89058 :   for (i = 1; i < lg(fa); i++)
    2807             :   {
    2808       44403 :     GEN Pi = gel(fa, i), soli, soli2 = NULL;
    2809       44403 :     long n, nrri = 0, h;
    2810       44403 :     if (ab)
    2811       38241 :       h = 1;
    2812             :     else
    2813        6162 :       Pi = ZX_deflate_max(Pi, &h);
    2814       44403 :     soli = ZX_Uspensky(Pi, h%2 ? ab: gen_0, 1, prec2nbits(prec));
    2815       44403 :     n = lg(soli);
    2816       44403 :     if (!(h % 2)) soli2 = cgetg(n, t_COL);
    2817      102945 :     for (j = 1; j < n; j++)
    2818             :     {
    2819       58542 :       GEN elt = gel(soli, j); /* != 0 */
    2820       58542 :       if (typ(elt) != t_REAL)
    2821             :       {
    2822          84 :         nrri++; if (h > 1 && !(h % 2)) nrri++;
    2823          84 :         elt = gtofp(elt, prec);
    2824          84 :         gel(soli, j) = elt;
    2825             :       }
    2826       58542 :       if (h > 1)
    2827             :       {
    2828             :         GEN r;
    2829        3222 :         if (h == 2)
    2830        3208 :           r = sqrtr(elt);
    2831             :         else
    2832             :         {
    2833          14 :           if (signe(elt) < 0)
    2834           7 :             r = negr(sqrtnr(negr(elt), h));
    2835             :           else
    2836           7 :             r = sqrtnr(elt, h);
    2837             :         }
    2838        3222 :         gel(soli, j) = r;
    2839        3222 :         if (!(h % 2)) gel(soli2, j) = negr(r);
    2840             :       }
    2841             :     }
    2842       44403 :     if (!(h % 2)) soli = shallowconcat(soli, soli2);
    2843       44403 :     if (ex[i] > 1) soli = shallowconcat1( const_vec(ex[i], soli) );
    2844       44403 :     sol = sol? shallowconcat(sol, soli): soli;
    2845       44403 :     nrr += ex[i]*nrri;
    2846             :   }
    2847       44655 :   if (!sol) { set_avma(av); return cgetg(1,t_COL); }
    2848             : 
    2849       44641 :   if (DEBUGLEVEL > 4)
    2850             :   {
    2851           0 :     err_printf("Number of real roots: %d\n", lg(sol)-1);
    2852           0 :     err_printf(" -- of which 2-integral: %ld\n", nrr);
    2853             :   }
    2854       44641 :   return gerepileupto(av, sort(sol));
    2855             : }
    2856             : 
    2857             : /* P non-constant, squarefree ZX */
    2858             : long
    2859       19327 : ZX_sturm(GEN P)
    2860             : {
    2861       19327 :   pari_sp av = avma;
    2862             :   long h, r;
    2863       19327 :   P = ZX_deflate_max(P, &h);
    2864       19327 :   if (odd(h))
    2865       12194 :     r = itos(ZX_Uspensky(P, NULL, 2, 0));
    2866             :   else
    2867        7133 :     r = 2*itos(ZX_Uspensky(P, gen_0, 2, 0));
    2868       19327 :   return gc_long(av,r);
    2869             : }
    2870             : /* P non-constant, squarefree ZX */
    2871             : long
    2872        2447 : ZX_sturmpart(GEN P, GEN ab)
    2873             : {
    2874        2447 :   pari_sp av = avma;
    2875        2447 :   if (!check_ab(ab)) return ZX_sturm(P);
    2876        2153 :   return gc_long(av, itou(ZX_Uspensky(P, ab, 2, 0)));
    2877             : }

Generated by: LCOV version 1.13