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

Generated by: LCOV version 1.11