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-bordeaux1.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 17909-2b849a1) Lines: 1541 1695 90.9 %
Date: 2015-07-05 Functions: 110 116 94.8 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 882 1061 83.1 %

           Branch data     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                 :   11518741 : addCC(GEN x, GEN y)
      35                 :            : {
      36                 :            :   GEN z;
      37                 :            : 
      38         [ +  + ]:   11518741 :   if (typ(x) == t_INT)
      39                 :            :   {
      40         [ +  + ]:   10442146 :     if (typ(y) == t_INT) return addii(x,y);
      41                 :            :     /* ty == t_COMPLEX */
      42                 :       1036 :     z = cgetg(3,t_COMPLEX);
      43                 :       1036 :     gel(z,1) = addii(x, gel(y,1));
      44                 :       1036 :     gel(z,2) = icopy(gel(y,2)); return z;
      45                 :            :   }
      46                 :            :   /* tx == t_COMPLEX */
      47                 :    1076595 :   z = cgetg(3,t_COMPLEX);
      48         [ +  + ]:    1076595 :   if (typ(y) == t_INT)
      49                 :            :   {
      50                 :      14667 :     gel(z,1) = addii(gel(x,1),y);
      51                 :      14667 :     gel(z,2) = icopy(gel(x,2)); return z;
      52                 :            :   }
      53                 :            :   /* ty == t_COMPLEX */
      54                 :    1061928 :   gel(z,1) = addii(gel(x,1),gel(y,1));
      55                 :   11518741 :   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                 :   21134436 : mulCC(GEN x, GEN y)
      60                 :            : {
      61                 :            :   GEN z;
      62                 :            : 
      63         [ +  + ]:   21134436 :   if (typ(x) == t_INT)
      64                 :            :   {
      65         [ +  + ]:   18543293 :     if (typ(y) == t_INT) return mulii(x,y);
      66                 :            :     /* ty == t_COMPLEX */
      67                 :       3744 :     z = cgetg(3,t_COMPLEX);
      68                 :       3744 :     gel(z,1) = mulii(x, gel(y,1));
      69                 :       3744 :     gel(z,2) = mulii(x, gel(y,2)); return z;
      70                 :            :   }
      71                 :            :   /* tx == t_COMPLEX */
      72                 :    2591143 :   z = cgetg(3,t_COMPLEX);
      73         [ +  + ]:    2591143 :   if (typ(y) == t_INT)
      74                 :            :   {
      75                 :     747304 :     gel(z,1) = mulii(gel(x,1),y);
      76                 :     747304 :     gel(z,2) = mulii(gel(x,2),y); return z;
      77                 :            :   }
      78                 :            :   /* ty == t_COMPLEX */
      79                 :            :   {
      80                 :    1843839 :     pari_sp av = avma, tetpil;
      81                 :            :     GEN p1, p2;
      82                 :            : 
      83                 :    1843839 :     p1 = mulii(gel(x,1),gel(y,1));
      84                 :    1843839 :     p2 = mulii(gel(x,2),gel(y,2));
      85                 :    1843839 :     y = mulii(addii(gel(x,1),gel(x,2)),
      86                 :    3687678 :               addii(gel(y,1),gel(y,2)));
      87                 :    1843839 :     x = addii(p1,p2); tetpil = avma;
      88                 :    1843839 :     gel(z,1) = subii(p1,p2);
      89                 :    1843839 :     gel(z,2) = subii(y,x); gerepilecoeffssp(av,tetpil,z+1,2);
      90                 :   21134436 :     return z;
      91                 :            :   }
      92                 :            : }
      93                 :            : /* fast squaring x: t_INT or t_COMPLEX(t_INT) */
      94                 :            : static GEN
      95                 :   16377680 : sqrCC(GEN x)
      96                 :            : {
      97                 :            :   GEN z;
      98                 :            : 
      99         [ +  + ]:   16377680 :   if (typ(x) == t_INT) return sqri(x);
     100                 :            :   /* tx == t_COMPLEX */
     101                 :    2506896 :   z = cgetg(3,t_COMPLEX);
     102                 :            :   {
     103                 :    2506896 :     pari_sp av = avma, tetpil;
     104                 :            :     GEN y, p1, p2;
     105                 :            : 
     106                 :    2506896 :     p1 = sqri(gel(x,1));
     107                 :    2506896 :     p2 = sqri(gel(x,2));
     108                 :    2506896 :     y = sqri(addii(gel(x,1),gel(x,2)));
     109                 :    2506896 :     x = addii(p1,p2); tetpil = avma;
     110                 :    2506896 :     gel(z,1) = subii(p1,p2);
     111                 :    2506896 :     gel(z,2) = subii(y,x); gerepilecoeffssp(av,tetpil,z+1,2);
     112                 :   16377680 :     return z;
     113                 :            :   }
     114                 :            : }
     115                 :            : 
     116                 :            : static void
     117                 :    2571917 : set_karasquare_limit(long bit)
     118                 :            : {
     119         [ +  + ]:    2571917 :   if (bit<600)       { KARASQUARE_LIMIT=8; COOKSQUARE_LIMIT=400; }
     120         [ +  - ]:       2184 :   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                 :    2571917 : }
     125                 :            : 
     126                 :            : /* assume lP > 0, lP = lgpol(P) */
     127                 :            : static GEN
     128                 :    5506002 : CX_square_spec(GEN P, long lP)
     129                 :            : {
     130                 :            :   GEN s, t;
     131                 :    5506002 :   long i, j, l, nn, n = lP - 1;
     132                 :            :   pari_sp av;
     133                 :            : 
     134                 :    5506002 :   nn = n<<1; s = cgetg(nn+3,t_POL); s[1] = evalsigne(1)|evalvarn(0);
     135                 :    5506002 :   gel(s,2) = sqrCC(gel(P,0)); /* i = 0 */
     136         [ +  + ]:   15121697 :   for (i=1; i<=n; i++)
     137                 :            :   {
     138                 :    9615695 :     av = avma; l = (i+1)>>1;
     139                 :    9615695 :     t = mulCC(gel(P,0), gel(P,i)); /* j = 0 */
     140         [ +  + ]:   13643482 :     for (j=1; j<l; j++) t = addCC(t, mulCC(gel(P,j), gel(P,i-j)));
     141                 :    9615695 :     t = gmul2n(t,1);
     142         [ +  + ]:    9615695 :     if ((i & 1) == 0) t = addCC(t, sqrCC(gel(P,i>>1)));
     143                 :    9615695 :     gel(s,i+2) = gerepileupto(av, t);
     144                 :            :   }
     145                 :    5506002 :   gel(s,nn+2) = sqrCC(gel(P,n)); /* i = nn */
     146         [ +  + ]:   10871678 :   for (   ; i<nn; i++)
     147                 :            :   {
     148                 :    5365676 :     av = avma; l = (i+1)>>1;
     149                 :    5365676 :     t = mulCC(gel(P,i-n),gel(P,n)); /* j = i-n */
     150         [ +  + ]:    7490954 :     for (j=i-n+1; j<l; j++) t = addCC(t, mulCC(gel(P,j),gel(P,i-j)));
     151                 :    5365676 :     t = gmul2n(t,1);
     152         [ +  + ]:    5365676 :     if ((i & 1) == 0) t = addCC(t, sqrCC(gel(P,i>>1)));
     153                 :    5365676 :     gel(s,i+2) = gerepileupto(av, t);
     154                 :            :   }
     155                 :    5506002 :   return normalizepol_lg(s, nn+3);
     156                 :            : }
     157                 :            : /* not stack clean */
     158                 :            : static GEN
     159                 :     181658 : RgX_addspec(GEN x, long nx, GEN y, long ny)
     160                 :            : {
     161                 :            :   GEN z, t;
     162                 :            :   long i;
     163         [ +  + ]:     181658 :   if (nx == ny) {
     164                 :      90230 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
     165         [ +  + ]:     694236 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
     166                 :      90230 :     return normalizepol_lg(z, nx+2);
     167                 :            :   }
     168         [ +  - ]:      91428 :   if (ny < nx) {
     169                 :      91428 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
     170         [ +  + ]:     640065 :     for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
     171         [ +  + ]:     182856 :     for (   ; i < nx; i++) gel(t,i) = gel(x,i);
     172                 :      91428 :     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                 :     181658 :     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                 :    5688808 : karasquare(GEN P, long nP)
     206                 :            : {
     207                 :            :   GEN Q, s0, s1, s2, a, t;
     208                 :    5688808 :   long n0, n1, i, l, N, N0, N1, n = nP - 1; /* degree(P) */
     209                 :            :   pari_sp av;
     210                 :            : 
     211 [ +  + ][ +  + ]:    5688808 :   if (n <= KARASQUARE_LIMIT) return nP? CX_square_spec(P, nP): pol_0(0);
     212                 :     181658 :   av = avma;
     213                 :     181658 :   n0 = (n>>1) + 1; n1 = nP - n0;
     214                 :     181658 :   s0 = karasquare(P, n0); Q = P + n0;
     215                 :     181658 :   s2 = karasquare(Q, n1);
     216                 :     181658 :   s1 = RgX_addspec(P, n0, Q, n1);
     217                 :     181658 :   s1 = RgX_sub(karasquare(s1+2, lgpol(s1)), RgX_add(s0,s2));
     218                 :     181658 :   N = (n<<1) + 1;
     219                 :     181658 :   a = cgetg(N + 2, t_POL); a[1] = evalsigne(1)|evalvarn(0);
     220                 :     181658 :   t = a+2; l = lgpol(s0); s0 += 2; N0 = n0<<1;
     221         [ +  + ]:    2440304 :   for (i=0; i < l;  i++) gel(t,i) = gel(s0,i);
     222         [ +  + ]:     411154 :   for (   ; i < N0; i++) gel(t,i) = gen_0;
     223                 :     181658 :   t = a+2 + N0; l = lgpol(s2); s2 += 2; N1 = N - N0;
     224         [ +  + ]:    2170510 :   for (i=0; i < l;  i++) gel(t,i) = gel(s2,i);
     225         [ +  + ]:     316434 :   for (   ; i < N1; i++) gel(t,i) = gen_0;
     226                 :     181658 :   t = a+2 + n0; l = lgpol(s1); s1 += 2;
     227         [ +  + ]:    2292597 :   for (i=0; i < l; i++)  gel(t,i) = gadd(gel(t,i), gel(s1,i));
     228                 :    5688808 :   return gerepilecopy(av, normalizepol_lg(a, N+2));
     229                 :            : }
     230                 :            : /* spec function. nP = lgpol(P) */
     231                 :            : static GEN
     232                 :    5143834 : cook_square(GEN P, long nP)
     233                 :            : {
     234                 :            :   GEN Q, p0, p1, p2, p3, q, r, t, vp, vm;
     235                 :    5143834 :   long n0, n3, i, j, n = nP - 1;
     236                 :            :   pari_sp av;
     237                 :            : 
     238 [ +  - ][ +  - ]:    5143834 :   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                 :    5143834 :   return gerepilecopy(av, normalizepol_lg(q, 2*n+3));
     303                 :            : }
     304                 :            : 
     305                 :            : static GEN
     306                 :    2571917 : graeffe(GEN p)
     307                 :            : {
     308                 :            :   GEN p0, p1, s0, s1;
     309                 :    2571917 :   long n = degpol(p), n0, n1, i;
     310                 :            : 
     311         [ -  + ]:    2571917 :   if (!n) return gcopy(p);
     312                 :    2571917 :   n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
     313                 :    2571917 :   p0 = new_chunk(n0);
     314                 :    2571917 :   p1 = new_chunk(n1);
     315         [ +  + ]:    8868543 :   for (i=0; i<n1; i++)
     316                 :            :   {
     317                 :    6296626 :     p0[i] = p[2+(i<<1)];
     318                 :    6296626 :     p1[i] = p[3+(i<<1)];
     319                 :            :   }
     320         [ +  + ]:    2571917 :   if (n1 != n0)
     321                 :    1299942 :     p0[i] = p[2+(i<<1)];
     322                 :    2571917 :   s0 = cook_square(p0, n0);
     323                 :    2571917 :   s1 = cook_square(p1, n1);
     324                 :    2571917 :   return RgX_sub(s0, RgX_shift_shallow(s1,1));
     325                 :            : }
     326                 :            : GEN
     327                 :       4655 : ZX_graeffe(GEN p)
     328                 :            : {
     329                 :       4655 :   pari_sp av = avma;
     330                 :            :   GEN p0, p1, s0, s1;
     331                 :       4655 :   long n = degpol(p);
     332                 :            : 
     333         [ -  + ]:       4655 :   if (!n) return ZX_copy(p);
     334                 :       4655 :   RgX_even_odd(p, &p0, &p1);
     335                 :            :   /* p = p0(x^2) + x p1(x^2) */
     336                 :       4655 :   s0 = ZX_sqr(p0);
     337                 :       4655 :   s1 = ZX_sqr(p1);
     338                 :       4655 :   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                 :   14456780 : mydbllog2i(GEN x)
     368                 :            : {
     369                 :            : #ifdef LONG_IS_64BIT
     370                 :   12395644 :   const double W = 1/(4294967296. * 4294967296.); /* 2^-64 */
     371                 :            : #else
     372                 :    2061136 :   const double W = 1/4294967296.; /*2^-32*/
     373                 :            : #endif
     374                 :            :   GEN m;
     375                 :   14456780 :   long lx = lgefint(x);
     376                 :            :   double l;
     377         [ +  + ]:   14456780 :   if (lx == 2) return -pariINFINITY;
     378                 :   14357718 :   m = int_MSW(x);
     379                 :   14357718 :   l = (double)(ulong)*m;
     380         [ +  + ]:   14357718 :   if (lx == 3) return log2(l);
     381                 :    6816063 :   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                 :   14456780 :   return log2(l) + (double)(BITS_IN_LONG*(lx-3));
     386                 :            : }
     387                 :            : 
     388                 :            : /* return log(|x|) or -pariINFINITY */
     389                 :            : static double
     390                 :     868355 : mydbllogr(GEN x) {
     391         [ -  + ]:     868355 :   if (!signe(x)) return -pariINFINITY;
     392                 :     868355 :   return LOG2*dbllog2r(x);
     393                 :            : }
     394                 :            : 
     395                 :            : /* return log2(|x|) or -pariINFINITY */
     396                 :            : static double
     397                 :    7128573 : mydbllog2r(GEN x) {
     398         [ +  + ]:    7128573 :   if (!signe(x)) return -pariINFINITY;
     399                 :    7128573 :   return dbllog2r(x);
     400                 :            : }
     401                 :            : static double
     402         [ +  + ]:    3940156 : dbllog2mp(GEN x) { return typ(x) == t_INT? mydbllog2i(x): mydbllog2r(x); }
     403                 :            : double
     404                 :   19236902 : dbllog2(GEN z)
     405                 :            : {
     406                 :            :   double x, y;
     407   [ +  +  +  + ]:   19236902 :   switch(typ(z))
     408                 :            :   {
     409                 :   11143084 :     case t_INT: return mydbllog2i(z);
     410                 :          7 :     case t_FRAC: return mydbllog2i(gel(z,1))-mydbllog2i(gel(z,2));
     411                 :    6123733 :     case t_REAL: return mydbllog2r(z);
     412                 :            :     default: /*t_COMPLEX*/
     413                 :    1970078 :       x = dbllog2mp(gel(z,1));
     414                 :    1970078 :       y = dbllog2mp(gel(z,2));
     415         [ +  + ]:    1970078 :       if (fabs(x-y) > 10) return maxdd(x,y);
     416                 :   19236902 :       return x + 0.5*log2(1 + exp2(2*(y-x)));
     417                 :            :   }
     418                 :            : }
     419                 :            : static GEN /* beware overflow */
     420         [ +  + ]:     400197 : 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                 :    1970063 : findpower(GEN p)
     426                 :            : {
     427                 :    1970063 :   double x, L, mins = pariINFINITY;
     428                 :    1970063 :   long n = degpol(p),i;
     429                 :            : 
     430                 :    1970063 :   L = dbllog2(gel(p,n+2)); /* log2(lc * binom(n,i)) */
     431         [ +  + ]:   10302077 :   for (i=n-1; i>=0; i--)
     432                 :            :   {
     433                 :    8332014 :     L += log2((double)(i+1) / (double)(n-i));
     434                 :    8332014 :     x = dbllog2(gel(p,i+2));
     435         [ +  + ]:    8332014 :     if (x != -pariINFINITY)
     436                 :            :     {
     437                 :    8271710 :       double s = (L - x) / (double)(n-i);
     438         [ +  + ]:    8271710 :       if (s < mins) mins = s;
     439                 :            :     }
     440                 :            :   }
     441                 :    1970063 :   i = (long)ceil(mins);
     442         [ +  + ]:    1970063 :   if (i - mins > 1 - 1e-12) i--;
     443                 :    1970063 :   return i;
     444                 :            : }
     445                 :            : 
     446                 :            : /* returns the exponent for logmodulus(), from the Newton diagram */
     447                 :            : static long
     448                 :     343962 : newton_polygon(GEN p, long k)
     449                 :            : {
     450                 :     343962 :   pari_sp av = avma;
     451                 :            :   double *logcoef, slope;
     452                 :     343962 :   long n = degpol(p), i, j, h, l, *vertex;
     453                 :            : 
     454                 :     343962 :   init_dalloc();
     455                 :     343962 :   logcoef = (double*)stack_malloc((n+1)*sizeof(double));
     456                 :     343962 :   vertex = (long*)new_chunk(n+1);
     457                 :            : 
     458                 :            :   /* vertex[i] = 1 if i a vertex of convex hull, 0 otherwise */
     459         [ +  + ]:    1897081 :   for (i=0; i<=n; i++) { logcoef[i] = dbllog2(gel(p,2+i)); vertex[i] = 0; }
     460                 :     343962 :   vertex[0] = 1; /* sentinel */
     461         [ +  + ]:    1461372 :   for (i=0; i < n; i=h)
     462                 :            :   {
     463                 :    1117410 :     slope = logcoef[i+1]-logcoef[i];
     464         [ +  + ]:    5385762 :     for (j = h = i+1; j<=n; j++)
     465                 :            :     {
     466                 :    4268352 :       double pij = (logcoef[j]-logcoef[i])/(double)(j-i);
     467         [ +  + ]:    4268352 :       if (slope < pij) { slope = pij; h = j; }
     468                 :            :     }
     469                 :    1117410 :     vertex[h] = 1;
     470                 :            :   }
     471         [ +  + ]:     370489 :   h = k;   while (!vertex[h]) h++;
     472         [ +  + ]:     361793 :   l = k-1; while (!vertex[l]) l--;
     473                 :     343962 :   avma = av;
     474                 :     343962 :   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                 :    2395339 : myshiftrc(GEN z, long e)
     480                 :            : {
     481         [ +  + ]:    2395339 :   if (typ(z)==t_COMPLEX)
     482                 :            :   {
     483         [ +  + ]:     464970 :     if (signe(gel(z,1))) shiftr_inplace(gel(z,1), e);
     484         [ +  + ]:     464970 :     if (signe(gel(z,2))) shiftr_inplace(gel(z,2), e);
     485                 :            :   }
     486                 :            :   else
     487         [ +  + ]:    1930369 :     if (signe(z)) shiftr_inplace(z, e);
     488                 :    2395339 : }
     489                 :            : 
     490                 :            : /* return z*2^e, where z is integer or complex of integer (destroy z) */
     491                 :            : static GEN
     492                 :    8904240 : myshiftic(GEN z, long e)
     493                 :            : {
     494         [ +  + ]:    8904240 :   if (typ(z)==t_COMPLEX)
     495                 :            :   {
     496         [ +  + ]:    1300021 :     gel(z,1) = signe(gel(z,1))? mpshift(gel(z,1),e): gen_0;
     497                 :    1300021 :     gel(z,2) = mpshift(gel(z,2),e);
     498                 :    1300021 :     return z;
     499                 :            :   }
     500         [ +  + ]:    8904240 :   return signe(z)? mpshift(z,e): gen_0;
     501                 :            : }
     502                 :            : 
     503                 :            : /* as real_1 with precision in bits, not in words */
     504                 :            : static GEN
     505                 :     403401 : myreal_1(long bit)
     506                 :            : {
     507         [ -  + ]:     403401 :   if (bit < 0) bit = 0;
     508                 :     403401 :   return real_1(nbits2prec(bit));
     509                 :            : }
     510                 :            : static GEN
     511                 :     434989 : RgX_gtofp_bit(GEN q, long bit)
     512                 :            : {
     513         [ -  + ]:     434989 :   if (bit < 0) bit = 0;
     514                 :     434989 :   return RgX_gtofp(q, nbits2prec(bit));
     515                 :            : }
     516                 :            : 
     517                 :            : static GEN
     518                 :   16187154 : mygprecrc(GEN x, long prec, long e)
     519                 :            : {
     520                 :            :   GEN y;
     521      [ +  +  + ]:   16187154 :   switch(typ(x))
     522                 :            :   {
     523         [ +  + ]:   12657078 :     case t_REAL: return signe(x)? rtor(x, prec): real_0_bit(e);
     524                 :            :     case t_COMPLEX:
     525                 :    2673647 :       y = cgetg(3,t_COMPLEX);
     526                 :    2673647 :       gel(y,1) = mygprecrc(gel(x,1),prec,e);
     527                 :    2673647 :       gel(y,2) = mygprecrc(gel(x,2),prec,e);
     528                 :    2673647 :       return y;
     529                 :   16187154 :     default: return gcopy(x);
     530                 :            :   }
     531                 :            : }
     532                 :            : 
     533                 :            : /* gprec behaves badly with the zero for polynomials.
     534                 :            : The second parameter in mygprec is the precision in base 2 */
     535                 :            : static GEN
     536                 :    4239605 : mygprec(GEN x, long bit)
     537                 :            : {
     538                 :            :   long lx, i, e, prec;
     539                 :            :   GEN y;
     540                 :            : 
     541         [ -  + ]:    4239605 :   if (bit < 0) bit = 0; /* should rarely happen */
     542                 :    4239605 :   e = gexpo(x) - bit;
     543                 :    4239605 :   prec = nbits2prec(bit);
     544         [ +  + ]:    4239605 :   switch(typ(x))
     545                 :            :   {
     546                 :            :     case t_POL:
     547                 :    2521311 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     548         [ +  + ]:   11565821 :       for (i=2; i<lx; i++) gel(y,i) = mygprecrc(gel(x,i),prec,e);
     549                 :    2521311 :       break;
     550                 :            : 
     551                 :    1718294 :     default: y = mygprecrc(x,prec,e);
     552                 :            :   }
     553                 :    4239605 :   return y;
     554                 :            : }
     555                 :            : 
     556                 :            : /* normalize a polynomial p, that is change it with coefficients in Z[i],
     557                 :            : after making product by 2^shift */
     558                 :            : static GEN
     559                 :    1032504 : pol_to_gaussint(GEN p, long shift)
     560                 :            : {
     561                 :    1032504 :   long i, l = lg(p);
     562                 :    1032504 :   GEN q = cgetg(l, t_POL); q[1] = p[1];
     563         [ +  + ]:    7734967 :   for (i=2; i<l; i++) gel(q,i) = gtrunc2n(gel(p,i), shift);
     564                 :    1032504 :   return q;
     565                 :            : }
     566                 :            : 
     567                 :            : /* returns a polynomial q in Z[i][x] keeping bit bits of p */
     568                 :            : static GEN
     569                 :     822097 : eval_rel_pol(GEN p, long bit)
     570                 :            : {
     571                 :            :   long i;
     572         [ +  + ]:    6078198 :   for (i = 2; i < lg(p); i++)
     573         [ +  + ]:    5256101 :     if (gequal0(gel(p,i))) gel(p,i) = gen_0; /* bad behaviour of gexpo */
     574                 :     822097 :   return pol_to_gaussint(p, bit-gexpo(p)+1);
     575                 :            : }
     576                 :            : 
     577                 :            : /* returns p(R*x)/R^n (in R or R[i]), R = exp(lrho), bit bits of precision */
     578                 :            : static GEN
     579                 :      91102 : homothetie(GEN p, double lrho, long bit)
     580                 :            : {
     581                 :            :   GEN q, r, t, iR;
     582                 :      91102 :   long n = degpol(p), i;
     583                 :            : 
     584                 :      91102 :   iR = mygprec(dblexp(-lrho),bit);
     585                 :      91102 :   q = mygprec(p, bit);
     586                 :      91102 :   r = cgetg(n+3,t_POL); r[1] = p[1];
     587                 :      91102 :   t = iR; r[n+2] = q[n+2];
     588         [ +  + ]:     612018 :   for (i=n-1; i>0; i--)
     589                 :            :   {
     590                 :     520916 :     gel(r,i+2) = gmul(t, gel(q,i+2));
     591                 :     520916 :     t = mulrr(t, iR);
     592                 :            :   }
     593                 :      91102 :   gel(r,2) = gmul(t, gel(q,2)); return r;
     594                 :            : }
     595                 :            : 
     596                 :            : /* change q in 2^(n*e) p(x*2^(-e)), n=deg(q)  [ ~as above with R = 2^-e ]*/
     597                 :            : static void
     598                 :     554369 : homothetie2n(GEN p, long e)
     599                 :            : {
     600         [ +  + ]:     554369 :   if (e)
     601                 :            :   {
     602                 :     417800 :     long i,n = lg(p)-1;
     603         [ +  + ]:    2813139 :     for (i=2; i<=n; i++) myshiftrc(gel(p,i), (n-i)*e);
     604                 :            :   }
     605                 :     554369 : }
     606                 :            : 
     607                 :            : /* return 2^f * 2^(n*e) p(x*2^(-e)), n=deg(q) */
     608                 :            : static void
     609                 :    1759656 : homothetie_gauss(GEN p, long e, long f)
     610                 :            : {
     611 [ +  + ][ +  + ]:    1759656 :   if (e || f)
     612                 :            :   {
     613                 :    1611878 :     long i, n = lg(p)-1;
     614         [ +  + ]:   10516118 :     for (i=2; i<=n; i++) gel(p,i) = myshiftic(gel(p,i), f+(n-i)*e);
     615                 :            :   }
     616                 :    1759656 : }
     617                 :            : 
     618                 :            : /* Lower bound on the modulus of the largest root z_0
     619                 :            :  * k is set to an upper bound for #{z roots, |z-z_0| < eps} */
     620                 :            : static double
     621                 :    1970063 : lower_bound(GEN p, long *k, double eps)
     622                 :            : {
     623                 :    1970063 :   long n = degpol(p), i, j;
     624                 :    1970063 :   pari_sp ltop = avma;
     625                 :            :   GEN a, s, S, ilc;
     626                 :            :   double r, R, rho;
     627                 :            : 
     628         [ +  + ]:    1970063 :   if (n < 4) { *k = n; return 0.; }
     629                 :     918438 :   S = cgetg(5,t_VEC);
     630                 :     918438 :   a = cgetg(5,t_VEC); ilc = gdiv(real_1(DEFAULTPREC), gel(p,n+2));
     631         [ +  + ]:    4592190 :   for (i=1; i<=4; i++) gel(a,i) = gmul(ilc,gel(p,n+2-i));
     632                 :            :   /* i = 1 split out from next loop for efficiency and initialization */
     633                 :     918438 :   s = gel(a,1);
     634                 :     918438 :   gel(S,1) = gneg(s); /* Newton sum S_i */
     635                 :     918438 :   rho = r = gtodouble(gabs(s,3));
     636                 :     918438 :   R = r / n;
     637         [ +  + ]:    3673752 :   for (i=2; i<=4; i++)
     638                 :            :   {
     639                 :    2755314 :     s = gmulsg(i,gel(a,i));
     640         [ +  + ]:    8265942 :     for (j=1; j<i; j++) s = gadd(s, gmul(gel(S,j),gel(a,i-j)));
     641                 :    2755314 :     gel(S,i) = gneg(s); /* Newton sum S_i */
     642                 :    2755314 :     r = gtodouble(gabs(s,3));
     643         [ +  + ]:    2755314 :     if (r > 0.)
     644                 :            :     {
     645                 :    2749335 :       r = exp(log(r/n) / (double)i);
     646         [ +  + ]:    2749335 :       if (r > R) R = r;
     647                 :            :     }
     648                 :            :   }
     649 [ +  + ][ +  + ]:     918438 :   if (R > 0. && eps < 1.2)
     650                 :     917668 :     *k = (long)floor((rho/R + n) / (1 + exp(-eps)*cos(eps)));
     651                 :            :   else
     652                 :        770 :     *k = n;
     653                 :    1970063 :   avma = ltop; return R;
     654                 :            : }
     655                 :            : 
     656                 :            : /* log of modulus of the largest root of p with relative error tau. Assume
     657                 :            :  * P(0) != 0 and P non constant */
     658                 :            : static double
     659                 :     210407 : logmax_modulus(GEN p, double tau)
     660                 :            : {
     661                 :            :   GEN r,q,aux,gunr;
     662                 :     210407 :   pari_sp av, ltop = avma;
     663                 :     210407 :   long i,k,n=degpol(p),nn,bit,M,e;
     664         [ -  + ]:     210407 :   double rho,eps, tau2 = (tau > 3.0)? 0.5: tau/6.;
     665                 :            : 
     666                 :     210407 :   r = cgeti(BIGDEFAULTPREC);
     667                 :     210407 :   av = avma;
     668                 :            : 
     669                 :     210407 :   eps = - 1/log(1.5*tau2); /* > 0 */
     670                 :     210407 :   bit = (long) ((double) n*log2(1./tau2)+3*log2((double) n))+1;
     671                 :     210407 :   gunr = myreal_1(bit+2*n);
     672                 :     210407 :   aux = gdiv(gunr, gel(p,2+n));
     673                 :     210407 :   q = RgX_Rg_mul(p, aux); gel(q,2+n) = gunr;
     674                 :     210407 :   e = findpower(q);
     675                 :     210407 :   homothetie2n(q,e);
     676                 :     210407 :   affsi(e, r);
     677                 :     210407 :   q = pol_to_gaussint(q, bit);
     678                 :     210407 :   M = (long) (log2( log(4.*n) / (2*tau2) )) + 2;
     679                 :     210407 :   nn = n;
     680                 :     210407 :   for (i=0,e=0;;)
     681                 :            :   { /* nn = deg(q) */
     682                 :    1970063 :     rho = lower_bound(q, &k, eps);
     683         [ +  + ]:    1970063 :     if (rho > exp2(-(double)e)) e = (long)-floor(log2(rho));
     684                 :    1970063 :     affii(shifti(addis(r,e), 1), r);
     685         [ +  + ]:    1970063 :     if (++i == M) break;
     686                 :            : 
     687                 :    5278968 :     bit = (long) ((double)k * log2(1./tau2) +
     688                 :    3519312 :                      (double)(nn-k)*log2(1./eps) + 3*log2((double)nn)) + 1;
     689                 :    1759656 :     homothetie_gauss(q, e, bit-(long)floor(dbllog2(gel(q,2+nn))+0.5));
     690                 :    1759656 :     nn -= RgX_valrem(q, &q);
     691                 :    1759656 :     set_karasquare_limit(gexpo(q));
     692                 :    1759656 :     q = gerepileupto(av, graeffe(q));
     693         [ +  + ]:    1759656 :     tau2 *= 1.5; if (tau2 > 0.9) tau2 = 0.5;
     694                 :    1759656 :     eps = -1/log(tau2); /* > 0 */
     695                 :    1759656 :     e = findpower(q);
     696                 :    1759656 :   }
     697         [ +  + ]:     210407 :   if (!signe(r)) { avma = ltop; return 0.; }
     698                 :     195578 :   r = itor(r, DEFAULTPREC); shiftr_inplace(r, -M);
     699                 :     210407 :   avma = ltop; return -rtodbl(r) * LOG2; /* -log(2) sum e_i 2^-i */
     700                 :            : }
     701                 :            : /* assume P non constant */
     702                 :            : GEN
     703                 :       1610 : logmax_modulus_bound(GEN P)
     704                 :            : {
     705                 :       1610 :   (void)RgX_valrem_inexact(P,&P);
     706         [ +  + ]:       1610 :   return lg(P)==3? gen_0: dblexp(logmax_modulus(P, 0.01) + 0.01);
     707                 :            : }
     708                 :            : 
     709                 :            : /* log of modulus of the smallest root of p, with relative error tau */
     710                 :            : static double
     711                 :      73498 : logmin_modulus(GEN p, double tau)
     712                 :            : {
     713                 :      73498 :   pari_sp av = avma;
     714                 :            :   double r;
     715                 :            : 
     716         [ -  + ]:      73498 :   if (gequal0(gel(p,2))) return -pariINFINITY;
     717                 :      73498 :   r = - logmax_modulus(RgX_recip_shallow(p),tau);
     718                 :      73498 :   avma = av; return r;
     719                 :            : }
     720                 :            : 
     721                 :            : /* return the log of the k-th modulus (ascending order) of p, rel. error tau*/
     722                 :            : static double
     723                 :      37694 : logmodulus(GEN p, long k, double tau)
     724                 :            : {
     725                 :            :   GEN q;
     726                 :      37694 :   long i, kk = k, imax, n = degpol(p), nn, bit, e;
     727                 :      37694 :   pari_sp av, ltop=avma;
     728                 :      37694 :   double r, tau2 = tau/6;
     729                 :            : 
     730                 :      37694 :   bit = (long)(n * (2. + log2(3.*n/tau2)));
     731                 :      37694 :   av = avma;
     732                 :      37694 :   q = gprec_w(p, nbits2prec(bit));
     733                 :      37694 :   q = RgX_gtofp_bit(q, bit);
     734                 :      37694 :   e = newton_polygon(q,k);
     735                 :      37694 :   r = (double)e;
     736                 :      37694 :   homothetie2n(q,e);
     737                 :      37694 :   imax = (long)(log2(3./tau) + log2(log(4.*n)))+1;
     738         [ +  + ]:     343962 :   for (i=1; i<imax; i++)
     739                 :            :   {
     740                 :     306268 :     q = eval_rel_pol(q,bit);
     741                 :     306268 :     kk -= RgX_valrem(q, &q);
     742                 :     306268 :     nn = degpol(q);
     743                 :            : 
     744                 :     306268 :     set_karasquare_limit(bit);
     745                 :     306268 :     q = gerepileupto(av, graeffe(q));
     746                 :     306268 :     e = newton_polygon(q,kk);
     747                 :     306268 :     r += e / exp2((double)i);
     748                 :     306268 :     q = RgX_gtofp_bit(q, bit);
     749                 :     306268 :     homothetie2n(q,e);
     750                 :            : 
     751         [ -  + ]:     306268 :     tau2 *= 1.5; if (tau2 > 1.) tau2 = 1.;
     752                 :     306268 :     bit = 1 + (long)(nn*(2. + log2(3.*nn/tau2)));
     753                 :            :   }
     754                 :      37694 :   avma = ltop; return -r * LOG2;
     755                 :            : }
     756                 :            : 
     757                 :            : /* return the log of the k-th modulus r_k of p, rel. error tau, knowing that
     758                 :            :  * rmin < r_k < rmax. This information helps because we may reduce precision
     759                 :            :  * quicker */
     760                 :            : static double
     761                 :      37694 : logpre_modulus(GEN p, long k, double tau, double lrmin, double lrmax)
     762                 :            : {
     763                 :            :   GEN q;
     764                 :      37694 :   long n = degpol(p), i, imax, imax2, bit;
     765                 :      37694 :   pari_sp ltop = avma, av;
     766                 :      37694 :   double lrho, aux, tau2 = tau/6.;
     767                 :            : 
     768                 :      37694 :   aux = (lrmax - lrmin) / 2. + 4*tau2;
     769                 :      37694 :   imax = (long) log2(log((double)n)/ aux);
     770         [ +  + ]:      37694 :   if (imax <= 0) return logmodulus(p,k,tau);
     771                 :            : 
     772                 :      36827 :   lrho  = (lrmin + lrmax) / 2;
     773                 :      36827 :   av = avma;
     774                 :      36827 :   bit = (long)(n*(2. + aux / LOG2 - log2(tau2)));
     775                 :      36827 :   q = homothetie(p, lrho, bit);
     776                 :      36827 :   imax2 = (long)(log2(3./tau * log(4.*n))) + 1;
     777         [ -  + ]:      36827 :   if (imax > imax2) imax = imax2;
     778                 :            : 
     779         [ +  + ]:     116118 :   for (i=0; i<imax; i++)
     780                 :            :   {
     781                 :      79291 :     q = eval_rel_pol(q,bit);
     782                 :      79291 :     set_karasquare_limit(bit);
     783                 :      79291 :     q = gerepileupto(av, graeffe(q));
     784                 :      79291 :     aux = 2*aux + 2*tau2;
     785                 :      79291 :     tau2 *= 1.5;
     786                 :      79291 :     bit = (long)(n*(2. + aux / LOG2 - log2(1-exp(-tau2))));
     787                 :      79291 :     q = RgX_gtofp_bit(q, bit);
     788                 :            :   }
     789                 :      36827 :   aux = exp2((double)imax);
     790                 :      36827 :   aux = logmodulus(q,k, aux*tau/3.) / aux;
     791                 :      37694 :   avma = ltop; return lrho + aux;
     792                 :            : }
     793                 :            : 
     794                 :            : static double
     795                 :      44439 : ind_maxlog2(GEN q)
     796                 :            : {
     797                 :      44439 :   long i, k = -1;
     798                 :      44439 :   double L = - pariINFINITY;
     799         [ +  + ]:     124463 :   for (i=0; i<=degpol(q); i++)
     800                 :            :   {
     801                 :      80024 :     double d = dbllog2(gel(q,2+i));
     802         [ +  + ]:      80024 :     if (d > L) { L = d; k = i; }
     803                 :            :   }
     804                 :      44439 :   return k;
     805                 :            : }
     806                 :            : 
     807                 :            : /* Returns k such that r_k e^(-tau) < R < r_{k+1} e^tau.
     808                 :            :  * Assume that l <= k <= n-l */
     809                 :            : static long
     810                 :      54275 : dual_modulus(GEN p, double lrho, double tau, long l)
     811                 :            : {
     812                 :      54275 :   long i, imax, delta_k = 0, n = degpol(p), nn, v2, v, bit, ll = l;
     813                 :      54275 :   double tau2 = tau * 7./8.;
     814                 :      54275 :   pari_sp av = avma;
     815                 :            :   GEN q;
     816                 :            : 
     817                 :      54275 :   bit = 6*n - 5*l + (long)(n*(-log2(tau2) + tau2 * 8./7.));
     818                 :      54275 :   q = homothetie(p, lrho, bit);
     819                 :      54275 :   imax = (long)(log(log(2.*n)/tau2)/log(7./4.)+1);
     820                 :            : 
     821         [ +  + ]:     480977 :   for (i=0; i<imax; i++)
     822                 :            :   {
     823                 :     436538 :     q = eval_rel_pol(q,bit); v2 = n - degpol(q);
     824                 :     436538 :     v = RgX_valrem(q, &q);
     825         [ +  + ]:     436538 :     ll -= maxss(v, v2); if (ll < 0) ll = 0;
     826                 :            : 
     827                 :     436538 :     nn = degpol(q); delta_k += v;
     828         [ +  + ]:     436538 :     if (!nn) return delta_k;
     829                 :            : 
     830                 :     426702 :     set_karasquare_limit(bit);
     831                 :     426702 :     q = gerepileupto(av, graeffe(q));
     832                 :     426702 :     tau2 *= 7./4.;
     833                 :     426702 :     bit = 6*nn - 5*ll + (long)(nn*(-log2(tau2) + tau2 * 8./7.));
     834                 :            :   }
     835                 :      54275 :   avma = av; return delta_k + (long)ind_maxlog2(q);
     836                 :            : }
     837                 :            : 
     838                 :            : /********************************************************************/
     839                 :            : /**                                                                **/
     840                 :            : /**              FACTORS THROUGH CIRCLE INTEGRATION                **/
     841                 :            : /**                                                                **/
     842                 :            : /********************************************************************/
     843                 :            : /* l power of 2 */
     844                 :            : static void
     845                 :    1957422 : fft(GEN Omega, GEN p, GEN f, long step, long l)
     846                 :            : {
     847                 :            :   pari_sp ltop;
     848                 :            :   long i, l1, l2, l3, rapi, step4;
     849                 :            :   GEN f1, f2, f3, f02, f13, g02, g13, ff;
     850                 :            : 
     851         [ +  + ]:    1957422 :   if (l == 2)
     852                 :            :   {
     853                 :    1136688 :     gel(f,0) = gadd(gel(p,0),gel(p,step));
     854                 :    1136688 :     gel(f,1) = gsub(gel(p,0),gel(p,step)); return;
     855                 :            :   }
     856         [ +  + ]:     820734 :   if (l == 4)
     857                 :            :   {
     858                 :     408991 :     f1 = gadd(gel(p,0),   gel(p,step<<1));
     859                 :     408991 :     f2 = gsub(gel(p,0),   gel(p,step<<1));
     860                 :     408991 :     f3 = gadd(gel(p,step),gel(p,3*step));
     861                 :     408991 :     f02= gsub(gel(p,step),gel(p,3*step));
     862                 :     408991 :     f02 = mulcxI(f02);
     863                 :     408991 :     gel(f,0) = gadd(f1, f3);
     864                 :     408991 :     gel(f,1) = gadd(f2, f02);
     865                 :     408991 :     gel(f,2) = gsub(f1, f3);
     866                 :     408991 :     gel(f,3) = gsub(f2, f02); return;
     867                 :            :   }
     868                 :            : 
     869                 :     411743 :   ltop = avma;
     870                 :     411743 :   l1 = l>>2; l2 = 2*l1; l3 = l1+l2; step4 = step<<2;
     871                 :     411743 :   fft(Omega,p,          f,   step4,l1);
     872                 :     411743 :   fft(Omega,p+step,     f+l1,step4,l1);
     873                 :     411743 :   fft(Omega,p+(step<<1),f+l2,step4,l1);
     874                 :     411743 :   fft(Omega,p+3*step,   f+l3,step4,l1);
     875                 :            : 
     876                 :     411743 :   ff = cgetg(l+1,t_VEC);
     877         [ +  + ]:    1773047 :   for (i=0; i<l1; i++)
     878                 :            :   {
     879                 :    1361304 :     rapi = step*i;
     880                 :    1361304 :     f1 = gmul(gel(Omega,rapi),    gel(f,i+l1));
     881                 :    1361304 :     f2 = gmul(gel(Omega,rapi<<1), gel(f,i+l2));
     882                 :    1361304 :     f3 = gmul(gel(Omega,3*rapi),  gel(f,i+l3));
     883                 :            : 
     884                 :    1361304 :     f02 = gadd(gel(f,i),f2);
     885                 :    1361304 :     g02 = gsub(gel(f,i),f2);
     886                 :    1361304 :     f13 = gadd(f1,f3);
     887                 :    1361304 :     g13 = mulcxI(gsub(f1,f3));
     888                 :            : 
     889                 :    1361304 :     gel(ff,i+1)    = gadd(f02, f13);
     890                 :    1361304 :     gel(ff,i+l1+1) = gadd(g02, g13);
     891                 :    1361304 :     gel(ff,i+l2+1) = gsub(f02, f13);
     892                 :    1361304 :     gel(ff,i+l3+1) = gsub(g02, g13);
     893                 :            :   }
     894                 :     411743 :   ff = gerepilecopy(ltop,ff);
     895         [ +  + ]:    7402638 :   for (i=0; i<l; i++) f[i] = ff[i+1];
     896                 :            : }
     897                 :            : 
     898                 :            : /* e(1/N) */
     899                 :            : static GEN
     900                 :     107910 : RUgen(long N, long bit)
     901                 :            : {
     902         [ +  + ]:     107910 :   if (N == 2) return real_m1(nbits2prec(bit));
     903         [ +  + ]:     104148 :   if (N == 4) return gen_I();
     904                 :     107910 :   return expIr(divru(Pi2n(1, nbits2prec(bit)), N));
     905                 :            : }
     906                 :            : 
     907                 :            : /* N=2^k. returns a vector RU which contains exp(2*i*k*Pi/N), k=0..N-1 */
     908                 :            : static GEN
     909                 :      51359 : initRU(long N, long bit)
     910                 :            : {
     911                 :      51359 :   GEN *RU, z = RUgen(N, bit);
     912                 :      51359 :   long i, N2 = (N>>1), N4 = (N>>2), N8 = (N>>3);
     913                 :            : 
     914                 :      51359 :   RU = (GEN*)cgetg(N+1,t_VEC); RU++;
     915                 :            : 
     916                 :      51359 :   RU[0] = myreal_1(bit);
     917                 :      51359 :   RU[1] = z;
     918         [ +  + ]:      74312 :   for (i=1; i<N8; i++)
     919                 :            :   {
     920                 :      22953 :     GEN t = RU[i];
     921                 :      22953 :     RU[i+1] = gmul(z, t);
     922                 :      22953 :     RU[N4-i] = mkcomplex(gel(t,2), gel(t,1));
     923                 :            :   }
     924         [ +  + ]:     186891 :   for (i=0; i<N4; i++) RU[i+N4] = mulcxI(RU[i]);
     925         [ +  + ]:     322423 :   for (i=0; i<N2; i++) RU[i+N2] = gneg(RU[i]);
     926                 :      51359 :   return (GEN)RU;
     927                 :            : }
     928                 :            : 
     929                 :            : /* as above, N arbitrary */
     930                 :            : static GEN
     931                 :       5192 : initRUgen(long N, long bit)
     932                 :            : {
     933                 :       5192 :   GEN *RU = (GEN*)cgetg(N+1,t_VEC), z = RUgen(N,bit);
     934                 :       5192 :   long i, k = (N+3)>>1;
     935                 :       5192 :   RU[0] = gen_1;
     936                 :       5192 :   RU[1] = z;
     937         [ +  + ]:       7878 :   for (i=2; i<k; i++) RU[i] = gmul(z, RU[i-1]);
     938         [ +  + ]:       7247 :   for (   ; i<N; i++) RU[i] = gconj(RU[N-i]);
     939                 :       5192 :   return (GEN)RU;
     940                 :            : }
     941                 :            : 
     942                 :            : GEN
     943                 :          0 : FFTinit(long k, long prec)
     944                 :            : {
     945         [ #  # ]:          0 :   if (k <= 0) pari_err_DOMAIN("FFTinit", "k", "<=", gen_0, stoi(k));
     946                 :          0 :   return initRU(1L << k, prec2nbits(prec)) - 1;
     947                 :            : }
     948                 :            : 
     949                 :            : GEN
     950                 :          0 : FFT(GEN x, GEN Omega)
     951                 :            : {
     952                 :          0 :   long i, l = lg(Omega), n = lg(x);
     953                 :            :   GEN y, z;
     954         [ #  # ]:          0 :   if (!is_vec_t(typ(x))) pari_err_TYPE("FFT",x);
     955         [ #  # ]:          0 :   if (typ(Omega) != t_VEC) pari_err_TYPE("FFT",Omega);
     956         [ #  # ]:          0 :   if (n > l) pari_err_DIM("FFT");
     957                 :            : 
     958         [ #  # ]:          0 :   if (n < l) {
     959                 :          0 :     z = cgetg(l, t_VECSMALL); /* cf stackdummy */
     960         [ #  # ]:          0 :     for (i = 1; i < n; i++) z[i] = x[i];
     961         [ #  # ]:          0 :     for (     ; i < l; i++) gel(z,i) = gen_0;
     962                 :            :   }
     963                 :          0 :   else z = x;
     964                 :          0 :   y = cgetg(l, t_VEC);
     965                 :          0 :   fft(Omega+1, z+1, y+1, 1, l-1);
     966                 :          0 :   return y;
     967                 :            : }
     968                 :            : 
     969                 :            : /* returns 1 if p has only real coefficients, 0 else */
     970                 :            : static int
     971                 :      47188 : isreal(GEN p)
     972                 :            : {
     973                 :            :   long i;
     974         [ +  + ]:     316448 :   for (i = lg(p)-1; i > 1; i--)
     975         [ +  + ]:     280225 :     if (typ(gel(p,i)) == t_COMPLEX) return 0;
     976                 :      47188 :   return 1;
     977                 :            : }
     978                 :            : 
     979                 :            : /* x non complex */
     980                 :            : static GEN
     981                 :      35786 : abs_update_r(GEN x, double *mu) {
     982                 :      35786 :   GEN y = gtofp(x, DEFAULTPREC);
     983         [ +  + ]:      35786 :   double ly = mydbllogr(y); if (ly < *mu) *mu = ly;
     984                 :      35786 :   setabssign(y); return y;
     985                 :            : }
     986                 :            : /* return |x|, low accuracy. Set *mu = min(log(y), *mu) */
     987                 :            : static GEN
     988                 :     783100 : abs_update(GEN x, double *mu) {
     989                 :            :   GEN y, xr, yr;
     990                 :            :   double ly;
     991         [ +  + ]:     783100 :   if (typ(x) != t_COMPLEX) return abs_update_r(x, mu);
     992                 :     748308 :   xr = gel(x,1);
     993                 :     748308 :   yr = gel(x,2);
     994         [ +  + ]:     748308 :   if (gequal0(xr)) return abs_update_r(yr,mu);
     995         [ +  + ]:     748182 :   if (gequal0(yr)) return abs_update_r(xr,mu);
     996                 :            :   /* have to treat 0 specially: 0E-10 + 1e-20 = 0E-10 */
     997                 :     747314 :   xr = gtofp(xr, DEFAULTPREC);
     998                 :     747314 :   yr = gtofp(yr, DEFAULTPREC);
     999                 :     747314 :   y = sqrtr(addrr(sqrr(xr), sqrr(yr)));
    1000         [ +  + ]:     747314 :   ly = mydbllogr(y); if (ly < *mu) *mu = ly;
    1001                 :     783100 :   return y;
    1002                 :            : }
    1003                 :            : 
    1004                 :            : static void
    1005                 :      24619 : parameters(GEN p, long *LMAX, double *mu, double *gamma,
    1006                 :            :            int polreal, double param, double param2)
    1007                 :            : {
    1008                 :            :   GEN q, pc, Omega, A, RU, prim, g, ONE,TWO;
    1009                 :      24619 :   long n = degpol(p), bit, NN, K, i, j, Lmax;
    1010                 :      24619 :   pari_sp av2, av = avma;
    1011                 :            : 
    1012                 :      24619 :   bit = gexpo(p) + (long)param2+8;
    1013         [ +  + ]:      49809 :   Lmax = 4; while (Lmax <= n) Lmax <<= 1;
    1014         [ +  + ]:      24619 :   NN = (long)(param*3.14)+1; if (NN < Lmax) NN = Lmax;
    1015         [ +  + ]:      24619 :   K = NN/Lmax; if (K & 1) K++;
    1016                 :      24619 :   NN = Lmax*K;
    1017         [ +  + ]:      24619 :   if (polreal) K = K/2+1;
    1018                 :            : 
    1019                 :      24619 :   Omega = initRU(Lmax,bit);
    1020                 :      24619 :   prim = RUgen(NN, bit);
    1021                 :            : 
    1022                 :      24619 :   q = mygprec(p,bit) + 2;
    1023                 :      24619 :   A = cgetg(Lmax+1,t_VEC); A++;
    1024                 :      24619 :   pc= cgetg(Lmax+1,t_VEC); pc++;
    1025         [ +  + ]:     201725 :   for (i=0; i <= n; i++) gel(pc,i)= gel(q,i);
    1026         [ +  + ]:      91569 :   for (   ; i<Lmax; i++) gel(pc,i) = gen_0;
    1027                 :            : 
    1028                 :      24619 :   *mu = pariINFINITY;
    1029                 :      24619 :   g = real_0_bit(-bit);
    1030                 :      24619 :   ONE = real_1(DEFAULTPREC);
    1031                 :      24619 :   TWO = real2n(1, DEFAULTPREC);
    1032                 :      24619 :   av2 = avma;
    1033                 :      24619 :   RU = myreal_1(bit);
    1034         [ +  + ]:      95285 :   for (i=0; i<K; i++)
    1035                 :            :   {
    1036         [ +  + ]:      70666 :     if (i) {
    1037                 :      46047 :       GEN z = RU;
    1038         [ +  + ]:     355571 :       for (j=1; j<n; j++)
    1039                 :            :       {
    1040                 :     309524 :         gel(pc,j) = gmul(gel(q,j),z);
    1041                 :     309524 :         z = gmul(z,RU); /* RU = prim^i, z=prim^(ij) */
    1042                 :            :       }
    1043                 :      46047 :       gel(pc,n) = gmul(gel(q,n),z);
    1044                 :            :     }
    1045                 :            : 
    1046                 :      70666 :     fft(Omega,pc,A,1,Lmax);
    1047 [ +  + ][ +  + ]:      82952 :     if (polreal && i>0 && i<K-1)
                 [ +  + ]
    1048         [ +  + ]:     230970 :       for (j=0; j<Lmax; j++) g = addrr(g, divrr(TWO, abs_update(gel(A,j),mu)));
    1049                 :            :     else
    1050         [ +  + ]:     622796 :       for (j=0; j<Lmax; j++) g = addrr(g, divrr(ONE, abs_update(gel(A,j),mu)));
    1051                 :      70666 :     RU = gmul(RU, prim);
    1052         [ -  + ]:      70666 :     if (gc_needed(av,1))
    1053                 :            :     {
    1054         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"parameters");
    1055                 :          0 :       gerepileall(av2,2, &g,&RU);
    1056                 :            :     }
    1057                 :            :   }
    1058                 :      24619 :   *gamma = mydbllog2r(divru(g,NN));
    1059                 :      24619 :   *LMAX = Lmax; avma = av;
    1060                 :      24619 : }
    1061                 :            : 
    1062                 :            : /* NN is a multiple of Lmax */
    1063                 :            : static void
    1064                 :      26740 : dft(GEN p, long k, long NN, long Lmax, long bit, GEN F, GEN H, long polreal)
    1065                 :            : {
    1066                 :            :   GEN Omega, q, qd, pc, pd, A, B, C, RU, aux, U, W, prim, prim2;
    1067                 :      26740 :   long n = degpol(p), i, j, K;
    1068                 :            :   pari_sp ltop;
    1069                 :            : 
    1070                 :      26740 :   Omega = initRU(Lmax,bit);
    1071                 :      26740 :   prim = RUgen(NN, bit);
    1072                 :      26740 :   RU = cgetg(n+2,t_VEC); RU++;
    1073                 :            : 
    1074         [ +  + ]:      26740 :   K = NN/Lmax; if (polreal) K = K/2+1;
    1075                 :      26740 :   q = mygprec(p,bit);
    1076                 :      26740 :   qd = RgX_deriv(q);
    1077                 :            : 
    1078                 :      26740 :   A = cgetg(Lmax+1,t_VEC); A++;
    1079                 :      26740 :   B = cgetg(Lmax+1,t_VEC); B++;
    1080                 :      26740 :   C = cgetg(Lmax+1,t_VEC); C++;
    1081                 :      26740 :   pc = cgetg(Lmax+1,t_VEC); pc++;
    1082                 :      26740 :   pd = cgetg(Lmax+1,t_VEC); pd++;
    1083         [ +  + ]:     109500 :   pc[0] = q[2];  for (i=n+1; i<Lmax; i++) gel(pc,i) = gen_0;
    1084         [ +  + ]:     136240 :   pd[0] = qd[2]; for (i=n;   i<Lmax; i++) gel(pd,i) = gen_0;
    1085                 :            : 
    1086                 :      26740 :   ltop = avma;
    1087                 :      26740 :   W = cgetg(k+1,t_VEC);
    1088                 :      26740 :   U = cgetg(k+1,t_VEC);
    1089         [ +  + ]:      82873 :   for (i=1; i<=k; i++) gel(W,i) = gel(U,i) = gen_0;
    1090                 :            : 
    1091                 :      26740 :   gel(RU,0) = gen_1;
    1092                 :      26740 :   prim2 = myreal_1(bit);
    1093         [ +  + ]:      86686 :   for (i=0; i<K; i++)
    1094                 :            :   {
    1095                 :      59946 :     gel(RU,1) = prim2;
    1096         [ +  + ]:     496705 :     for (j=1; j<n; j++) gel(RU,j+1) = gmul(gel(RU,j),prim2);
    1097                 :            :     /* RU[j] = prim^(ij)= prim2^j */
    1098                 :            : 
    1099         [ +  + ]:     496705 :     for (j=1; j<n; j++) gel(pd,j) = gmul(gel(qd,j+2),gel(RU,j));
    1100                 :      59946 :     fft(Omega,pd,A,1,Lmax);
    1101         [ +  + ]:     556651 :     for (j=1; j<=n; j++) gel(pc,j) = gmul(gel(q,j+2),gel(RU,j));
    1102                 :      59946 :     fft(Omega,pc,B,1,Lmax);
    1103         [ +  + ]:     841506 :     for (j=0; j<Lmax; j++) gel(C,j) = ginv(gel(B,j));
    1104         [ +  + ]:     841506 :     for (j=0; j<Lmax; j++) gel(B,j) = gmul(gel(A,j),gel(C,j));
    1105                 :      59946 :     fft(Omega,B,A,1,Lmax);
    1106                 :      59946 :     fft(Omega,C,B,1,Lmax);
    1107                 :            : 
    1108         [ +  + ]:      59946 :     if (polreal) /* p has real coefficients */
    1109                 :            :     {
    1110 [ +  + ][ +  + ]:      50370 :       if (i>0 && i<K-1)
    1111                 :            :       {
    1112         [ +  + ]:      37556 :         for (j=1; j<=k; j++)
    1113                 :            :         {
    1114                 :      31811 :           gel(W,j) = gadd(gel(W,j), gshift(real_i(gmul(gel(A,j+1),gel(RU,j+1))),1));
    1115                 :      31811 :           gel(U,j) = gadd(gel(U,j), gshift(real_i(gmul(gel(B,j),gel(RU,j))),1));
    1116                 :            :         }
    1117                 :            :       }
    1118                 :            :       else
    1119                 :            :       {
    1120         [ +  + ]:     132513 :         for (j=1; j<=k; j++)
    1121                 :            :         {
    1122                 :      87888 :           gel(W,j) = gadd(gel(W,j), real_i(gmul(gel(A,j+1),gel(RU,j+1))));
    1123                 :      87888 :           gel(U,j) = gadd(gel(U,j), real_i(gmul(gel(B,j),gel(RU,j))));
    1124                 :            :         }
    1125                 :            :       }
    1126                 :            :     }
    1127                 :            :     else
    1128                 :            :     {
    1129         [ +  + ]:      43808 :       for (j=1; j<=k; j++)
    1130                 :            :       {
    1131                 :      28487 :         gel(W,j) = gadd(gel(W,j), gmul(gel(A,j+1),gel(RU,j+1)));
    1132                 :      28487 :         gel(U,j) = gadd(gel(U,j), gmul(gel(B,j),gel(RU,j)));
    1133                 :            :       }
    1134                 :            :     }
    1135                 :      59946 :     prim2 = gmul(prim2,prim);
    1136                 :      59946 :     gerepileall(ltop,3, &W,&U,&prim2);
    1137                 :            :   }
    1138                 :            : 
    1139         [ +  + ]:      82873 :   for (i=1; i<=k; i++)
    1140                 :            :   {
    1141                 :      56133 :     aux=gel(W,i);
    1142         [ +  + ]:     138835 :     for (j=1; j<i; j++) aux = gadd(aux, gmul(gel(W,i-j),gel(F,k+2-j)));
    1143                 :      56133 :     gel(F,k+2-i) = gdivgs(aux,-i*NN);
    1144                 :            :   }
    1145         [ +  + ]:      82873 :   for (i=0; i<k; i++)
    1146                 :            :   {
    1147                 :      56133 :     aux=gel(U,k-i);
    1148         [ +  + ]:     138835 :     for (j=1+i; j<k; j++) aux = gadd(aux,gmul(gel(F,2+j),gel(U,j-i)));
    1149                 :      56133 :     gel(H,i+2) = gdivgs(aux,NN);
    1150                 :            :   }
    1151                 :      26740 : }
    1152                 :            : 
    1153                 :            : #define NEWTON_MAX 10
    1154                 :            : static GEN
    1155                 :     135692 : refine_H(GEN F, GEN G, GEN HH, long bit, long Sbit)
    1156                 :            : {
    1157                 :     135692 :   GEN H = HH, D, aux;
    1158                 :     135692 :   pari_sp ltop = avma;
    1159                 :            :   long error, i, bit1, bit2;
    1160                 :            : 
    1161                 :     135692 :   D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F)); error = gexpo(D);
    1162                 :     135692 :   bit2 = bit + Sbit;
    1163 [ +  + ][ +  + ]:     254501 :   for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
                 [ +  + ]
    1164                 :            :   {
    1165         [ -  + ]:     118809 :     if (gc_needed(ltop,1))
    1166                 :            :     {
    1167         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"refine_H");
    1168                 :          0 :       gerepileall(ltop,2, &D,&H);
    1169                 :            :     }
    1170                 :     118809 :     bit1 = -error + Sbit;
    1171                 :     118809 :     aux = RgX_mul(mygprec(H,bit1), mygprec(D,bit1));
    1172                 :     118809 :     aux = RgX_rem(mygprec(aux,bit1), mygprec(F,bit1));
    1173                 :            : 
    1174         [ +  + ]:     118809 :     bit1 = -error*2 + Sbit; if (bit1 > bit2) bit1 = bit2;
    1175                 :     118809 :     H = RgX_add(mygprec(H,bit1), aux);
    1176                 :     118809 :     D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F));
    1177         [ +  + ]:     118809 :     error = gexpo(D); if (error < -bit1) error = -bit1;
    1178                 :            :   }
    1179         [ +  + ]:     135692 :   if (error > -bit/2) return NULL; /* FAIL */
    1180                 :     135692 :   return gerepilecopy(ltop,H);
    1181                 :            : }
    1182                 :            : 
    1183                 :            : /* return 0 if fails, 1 else */
    1184                 :            : static long
    1185                 :      26740 : refine_F(GEN p, GEN *F, GEN *G, GEN H, long bit, double gamma)
    1186                 :            : {
    1187                 :      26740 :   GEN f0, FF, GG, r, HH = H;
    1188                 :      26740 :   long error, i, bit1 = 0, bit2, Sbit, Sbit2,  enh, normF, normG, n = degpol(p);
    1189                 :      26740 :   pari_sp av = avma;
    1190                 :            : 
    1191                 :      26740 :   FF = *F; GG = RgX_divrem(p, FF, &r);
    1192         [ -  + ]:      26740 :   error = gexpo(r); if (error <= -bit) error = 1-bit;
    1193                 :      26740 :   normF = gexpo(FF);
    1194                 :      26740 :   normG = gexpo(GG);
    1195         [ +  + ]:      26740 :   enh = gexpo(H); if (enh < 0) enh = 0;
    1196                 :      26740 :   Sbit = normF + 2*normG + enh + (long)(4.*log2((double)n)+gamma) + 1;
    1197                 :      26740 :   Sbit2 = enh + 2*(normF+normG) + (long)(2.*gamma+5.*log2((double)n)) + 1;
    1198                 :      26740 :   bit2 = bit + Sbit;
    1199 [ +  + ][ +  + ]:     162383 :   for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
                 [ +  + ]
    1200                 :            :   {
    1201 [ +  + ][ +  + ]:     135692 :     if (bit1 == bit2 && i >= 2) { Sbit += n; Sbit2 += n; bit2 += n; }
    1202         [ -  + ]:     135692 :     if (gc_needed(av,1))
    1203                 :            :     {
    1204         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"refine_F");
    1205                 :          0 :       gerepileall(av,4, &FF,&GG,&r,&HH);
    1206                 :            :     }
    1207                 :            : 
    1208                 :     135692 :     bit1 = -error + Sbit2;
    1209                 :     135692 :     HH = refine_H(mygprec(FF,bit1), mygprec(GG,bit1), mygprec(HH,bit1),
    1210                 :            :                   1-error, Sbit2);
    1211         [ +  + ]:     135692 :     if (!HH) return 0; /* FAIL */
    1212                 :            : 
    1213                 :     135643 :     bit1 = -error + Sbit;
    1214                 :     135643 :     r = RgX_mul(mygprec(HH,bit1), mygprec(r,bit1));
    1215                 :     135643 :     f0 = RgX_rem(mygprec(r,bit1), mygprec(FF,bit1));
    1216                 :            : 
    1217         [ +  + ]:     135643 :     bit1 = -2*error + Sbit; if (bit1 > bit2) bit1 = bit2;
    1218                 :     135643 :     FF = gadd(mygprec(FF,bit1),f0);
    1219                 :            : 
    1220         [ +  + ]:     135643 :     bit1 = -3*error + Sbit; if (bit1 > bit2) bit1 = bit2;
    1221                 :     135643 :     GG = RgX_divrem(mygprec(p,bit1), mygprec(FF,bit1), &r);
    1222         [ +  + ]:     135643 :     error = gexpo(r); if (error < -bit1) error = -bit1;
    1223                 :            :   }
    1224         [ +  + ]:      26691 :   if (error>-bit) return 0; /* FAIL */
    1225                 :      26740 :   *F = FF; *G = GG; return 1;
    1226                 :            : }
    1227                 :            : 
    1228                 :            : /* returns F and G from the unit circle U such that |p-FG|<2^(-bit) |cd|,
    1229                 :            : where cd is the leading coefficient of p */
    1230                 :            : static void
    1231                 :      24619 : split_fromU(GEN p, long k, double delta, long bit,
    1232                 :            :             GEN *F, GEN *G, double param, double param2)
    1233                 :            : {
    1234                 :            :   GEN pp, FF, GG, H;
    1235                 :      24619 :   long n = degpol(p), NN, bit2, Lmax;
    1236                 :      24619 :   int polreal = isreal(p);
    1237                 :            :   pari_sp ltop;
    1238                 :            :   double mu, gamma;
    1239                 :            : 
    1240                 :      24619 :   pp = gdiv(p, gel(p,2+n));
    1241                 :      24619 :   parameters(pp, &Lmax,&mu,&gamma, polreal,param,param2);
    1242                 :            : 
    1243                 :      24619 :   H  = cgetg(k+2,t_POL); H[1] = p[1];
    1244                 :      24619 :   FF = cgetg(k+3,t_POL); FF[1]= p[1];
    1245                 :      24619 :   gel(FF,k+2) = gen_1;
    1246                 :            : 
    1247         [ +  + ]:      24619 :   NN = (long)(0.5/delta); NN |= 1; if (NN < 2) NN = 2;
    1248                 :      24619 :   NN *= Lmax; ltop = avma;
    1249                 :            :   for(;;)
    1250                 :            :   {
    1251                 :      26740 :     bit2 = (long)(((double)NN*delta-mu)/LOG2) + gexpo(pp) + 8;
    1252                 :      26740 :     dft(pp, k, NN, Lmax, bit2, FF, H, polreal);
    1253         [ +  + ]:      26740 :     if (refine_F(pp,&FF,&GG,H,bit,gamma)) break;
    1254                 :       2121 :     NN <<= 1; avma = ltop;
    1255                 :       2121 :   }
    1256                 :      24619 :   *G = gmul(GG,gel(p,2+n)); *F = FF;
    1257                 :      24619 : }
    1258                 :            : 
    1259                 :            : static void
    1260                 :      24619 : optimize_split(GEN p, long k, double delta, long bit,
    1261                 :            :             GEN *F, GEN *G, double param, double param2)
    1262                 :            : {
    1263                 :      24619 :   long n = degpol(p);
    1264                 :            :   GEN FF, GG;
    1265                 :            : 
    1266         [ +  + ]:      24619 :   if (k <= n/2)
    1267                 :      18861 :     split_fromU(p,k,delta,bit,F,G,param,param2);
    1268                 :            :   else
    1269                 :            :   {
    1270                 :       5758 :     split_fromU(RgX_recip_shallow(p),n-k,delta,bit,&FF,&GG,param,param2);
    1271                 :       5758 :     *F = RgX_recip_shallow(GG);
    1272                 :       5758 :     *G = RgX_recip_shallow(FF);
    1273                 :            :   }
    1274                 :      24619 : }
    1275                 :            : 
    1276                 :            : /********************************************************************/
    1277                 :            : /**                                                                **/
    1278                 :            : /**               SEARCH FOR SEPARATING CIRCLE                     **/
    1279                 :            : /**                                                                **/
    1280                 :            : /********************************************************************/
    1281                 :            : 
    1282                 :            : /* return p(2^e*x) *2^(-n*e) */
    1283                 :            : static void
    1284                 :          0 : scalepol2n(GEN p, long e)
    1285                 :            : {
    1286                 :          0 :   long i,n=lg(p)-1;
    1287         [ #  # ]:          0 :   for (i=2; i<=n; i++) gel(p,i) = gmul2n(gel(p,i),(i-n)*e);
    1288                 :          0 : }
    1289                 :            : 
    1290                 :            : /* returns p(x/R)*R^n */
    1291                 :            : static GEN
    1292                 :     209271 : scalepol(GEN p, GEN R, long bit)
    1293                 :            : {
    1294                 :            :   GEN q,aux,gR;
    1295                 :            :   long i;
    1296                 :            : 
    1297                 :     209271 :   aux = gR = mygprec(R,bit); q = mygprec(p,bit);
    1298         [ +  + ]:    1037129 :   for (i=lg(p)-2; i>=2; i--)
    1299                 :            :   {
    1300                 :     827858 :     gel(q,i) = gmul(aux,gel(q,i));
    1301                 :     827858 :     aux = gmul(aux,gR);
    1302                 :            :   }
    1303                 :     209271 :   return q;
    1304                 :            : }
    1305                 :            : 
    1306                 :            : /* return (conj(a)X-1)^n * p[ (X-a) / (conj(a)X-1) ] */
    1307                 :            : static GEN
    1308                 :      67707 : conformal_pol(GEN p, GEN a, long bit)
    1309                 :            : {
    1310                 :      67707 :   GEN z, r, ma = gneg(a), ca = gconj(a);
    1311                 :      67707 :   long n = degpol(p), i;
    1312                 :      67707 :   pari_sp av = avma;
    1313                 :            : 
    1314                 :      67707 :   z = mkpoln(2, ca, negr(myreal_1(bit)));
    1315                 :      67707 :   r = scalarpol(gel(p,2+n), 0);
    1316                 :      67707 :   for (i=n-1; ; i--)
    1317                 :            :   {
    1318                 :     261442 :     r = addmulXn(r, gmul(ma,r), 1); /* r *= (X - a) */
    1319                 :     261442 :     r = gadd(r, gmul(z, gel(p,2+i)));
    1320         [ +  + ]:     261442 :     if (i == 0) return gerepileupto(av, r);
    1321                 :     193735 :     z = addmulXn(gmul(z,ca), gneg(z), 1); /* z *= conj(a)X - 1 */
    1322         [ -  + ]:     193735 :     if (gc_needed(av,2))
    1323                 :            :     {
    1324         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"conformal_pol");
    1325                 :          0 :       gerepileall(av,2, &r,&z);
    1326                 :            :     }
    1327                 :     193735 :   }
    1328                 :            : }
    1329                 :            : 
    1330                 :            : static const double UNDEF = -100000.;
    1331                 :            : 
    1332                 :            : static double
    1333                 :      24619 : logradius(double *radii, GEN p, long k, double aux, double *delta)
    1334                 :            : {
    1335                 :      24619 :   long i, n = degpol(p);
    1336                 :            :   double lrho, lrmin, lrmax;
    1337         [ +  + ]:      24619 :   if (k > 1)
    1338                 :            :   {
    1339 [ +  - ][ +  + ]:      25994 :     i = k-1; while (i>0 && radii[i] == UNDEF) i--;
    1340                 :      16778 :     lrmin = logpre_modulus(p,k,aux, radii[i], radii[k]);
    1341                 :            :   }
    1342                 :            :   else /* k=1 */
    1343                 :       7841 :     lrmin = logmin_modulus(p,aux);
    1344                 :      24619 :   radii[k] = lrmin;
    1345                 :            : 
    1346         [ +  + ]:      24619 :   if (k+1<n)
    1347                 :            :   {
    1348 [ +  - ][ +  + ]:      47972 :     i = k+2; while (i<=n && radii[i] == UNDEF) i++;
    1349                 :      20916 :     lrmax = logpre_modulus(p,k+1,aux, radii[k+1], radii[i]);
    1350                 :            :   }
    1351                 :            :   else /* k+1=n */
    1352                 :       3703 :     lrmax = logmax_modulus(p,aux);
    1353                 :      24619 :   radii[k+1] = lrmax;
    1354                 :            : 
    1355                 :      24619 :   lrho = radii[k];
    1356         [ +  + ]:      60067 :   for (i=k-1; i>=1; i--)
    1357                 :            :   {
    1358 [ +  + ][ +  + ]:      35448 :     if (radii[i] == UNDEF || radii[i] > lrho)
    1359                 :      24600 :       radii[i] = lrho;
    1360                 :            :     else
    1361                 :      10848 :       lrho = radii[i];
    1362                 :            :   }
    1363                 :      24619 :   lrho = radii[k+1];
    1364         [ +  + ]:     117039 :   for (i=k+1; i<=n; i++)
    1365                 :            :   {
    1366 [ +  + ][ +  + ]:      92420 :     if (radii[i] == UNDEF || radii[i] < lrho)
    1367                 :      58142 :       radii[i] = lrho;
    1368                 :            :     else
    1369                 :      34278 :       lrho = radii[i];
    1370                 :            :   }
    1371                 :      24619 :   *delta = (lrmax - lrmin) / 2;
    1372         [ +  + ]:      24619 :   if (*delta > 1.) *delta = 1.;
    1373                 :      24619 :   return (lrmin + lrmax) / 2;
    1374                 :            : }
    1375                 :            : 
    1376                 :            : static void
    1377                 :      24619 : update_radius(long n, double *radii, double lrho, double *par, double *par2)
    1378                 :            : {
    1379                 :      24619 :   double t, param = 0., param2 = 0.;
    1380                 :            :   long i;
    1381         [ +  + ]:     177106 :   for (i=1; i<=n; i++)
    1382                 :            :   {
    1383                 :     152487 :     radii[i] -= lrho;
    1384                 :     152487 :     t = fabs(rtodbl( invr(subsr(1, dblexp(radii[i]))) ));
    1385         [ +  + ]:     152487 :     param += t; if (t > 1.) param2 += log2(t);
    1386                 :            :   }
    1387                 :      24619 :   *par = param; *par2 = param2;
    1388                 :      24619 : }
    1389                 :            : 
    1390                 :            : /* apply the conformal mapping then split from U */
    1391                 :            : static void
    1392                 :      22569 : conformal_mapping(double *radii, GEN ctr, GEN p, long k, long bit,
    1393                 :            :                   double aux, GEN *F,GEN *G)
    1394                 :            : {
    1395                 :      22569 :   long bit2, n = degpol(p), i;
    1396                 :      22569 :   pari_sp ltop = avma, av;
    1397                 :            :   GEN q, FF, GG, a, R;
    1398                 :            :   double lrho, delta, param, param2;
    1399                 :            :   /* n * (2.*log2(2.732)+log2(1.5)) + 1 */
    1400                 :      22569 :   bit2 = bit + (long)(n*3.4848775) + 1;
    1401                 :      22569 :   a = sqrtr_abs( stor(3, 2*MEDDEFAULTPREC - 2) );
    1402                 :      22569 :   a = divrs(a, -6);
    1403                 :      22569 :   a = gmul(mygprec(a,bit2), mygprec(ctr,bit2)); /* a = -ctr/2sqrt(3) */
    1404                 :            : 
    1405                 :      22569 :   av = avma;
    1406                 :      22569 :   q = conformal_pol(mygprec(p,bit2), a, bit2);
    1407         [ +  + ]:     153290 :   for (i=1; i<=n; i++)
    1408         [ +  + ]:     130721 :     if (radii[i] != UNDEF) /* update array radii */
    1409                 :            :     {
    1410                 :      85255 :       pari_sp av2 = avma;
    1411                 :      85255 :       GEN t, r = dblexp(radii[i]), r2 = sqrr(r);
    1412                 :            :       /* 2(r^2 - 1) / (r^2 - 3(r-1)) */
    1413                 :      85255 :       t = divrr(shiftr((subrs(r2,1)),1), subrr(r2, mulur(3,subrs(r,1))));
    1414                 :      85255 :       radii[i] = mydbllogr(addsr(1,t)) / 2;
    1415                 :      85255 :       avma = av2;
    1416                 :            :     }
    1417                 :      22569 :   lrho = logradius(radii, q,k,aux/10., &delta);
    1418                 :      22569 :   update_radius(n, radii, lrho, &param, &param2);
    1419                 :            : 
    1420                 :      22569 :   bit2 += (long)(n * fabs(lrho)/LOG2 + 1.);
    1421                 :      22569 :   R = mygprec(dblexp(-lrho), bit2);
    1422                 :      22569 :   q = scalepol(q,R,bit2);
    1423                 :      22569 :   gerepileall(av,2, &q,&R);
    1424                 :            : 
    1425                 :      22569 :   optimize_split(q,k,delta,bit2,&FF,&GG,param,param2);
    1426                 :      22569 :   bit2 += n; R = invr(R);
    1427                 :      22569 :   FF = scalepol(FF,R,bit2);
    1428                 :      22569 :   GG = scalepol(GG,R,bit2);
    1429                 :            : 
    1430                 :      22569 :   a = mygprec(a,bit2);
    1431                 :      22569 :   FF = conformal_pol(FF,a,bit2);
    1432                 :      22569 :   GG = conformal_pol(GG,a,bit2);
    1433                 :            : 
    1434                 :      22569 :   a = invr(subsr(1, gnorm(a)));
    1435                 :      22569 :   FF = RgX_Rg_mul(FF, powru(a,k));
    1436                 :      22569 :   GG = RgX_Rg_mul(GG, powru(a,n-k));
    1437                 :            : 
    1438                 :      22569 :   *F = mygprec(FF,bit+n);
    1439                 :      22569 :   *G = mygprec(GG,bit+n); gerepileall(ltop,2, F,G);
    1440                 :      22569 : }
    1441                 :            : 
    1442                 :            : /* split p, this time without scaling. returns in F and G two polynomials
    1443                 :            :  * such that |p-FG|< 2^(-bit)|p| */
    1444                 :            : static void
    1445                 :      24619 : split_2(GEN p, long bit, GEN ctr, double thickness, GEN *F, GEN *G)
    1446                 :            : {
    1447                 :            :   GEN q, FF, GG, R;
    1448                 :            :   double aux, delta, param, param2;
    1449                 :      24619 :   long n = degpol(p), i, j, k, bit2;
    1450                 :            :   double lrmin, lrmax, lrho, *radii;
    1451                 :            : 
    1452                 :      24619 :   init_dalloc();
    1453                 :      24619 :   radii = (double*) stack_malloc((n+1) * sizeof(double));
    1454                 :            : 
    1455         [ +  + ]:     127868 :   for (i=2; i<n; i++) radii[i] = UNDEF;
    1456                 :      24619 :   aux = thickness/(double)(4 * n);
    1457                 :      24619 :   lrmin = logmin_modulus(p, aux);
    1458                 :      24619 :   lrmax = logmax_modulus(p, aux);
    1459                 :      24619 :   radii[1] = lrmin;
    1460                 :      24619 :   radii[n] = lrmax;
    1461                 :      24619 :   i = 1; j = n;
    1462                 :      24619 :   lrho = (lrmin + lrmax) / 2;
    1463                 :      24619 :   k = dual_modulus(p, lrho, aux, 1);
    1464 [ +  + ][ +  + ]:      24619 :   if (5*k < n || (n < 2*k && 5*k < 4*n))
                 [ +  + ]
    1465                 :       5093 :     { lrmax = lrho; j=k+1; radii[j] = lrho; }
    1466                 :            :   else
    1467                 :      19526 :     { lrmin = lrho; i=k;   radii[i] = lrho; }
    1468         [ +  + ]:      54275 :   while (j > i+1)
    1469                 :            :   {
    1470         [ +  + ]:      29656 :     if (i+j == n+1)
    1471                 :      10419 :       lrho = (lrmin + lrmax) / 2;
    1472                 :            :     else
    1473                 :            :     {
    1474                 :      19237 :       double kappa = 2. - log(1. + minss(i,n-j)) / log(1. + minss(j,n-i));
    1475         [ +  + ]:      19237 :       if (i+j < n+1) lrho = lrmax * kappa + lrmin;
    1476                 :      13643 :       else           lrho = lrmin * kappa + lrmax;
    1477                 :      19237 :       lrho /= 1+kappa;
    1478                 :            :     }
    1479                 :      29656 :     aux = (lrmax - lrmin) / (4*(j-i));
    1480                 :      29656 :     k = dual_modulus(p, lrho, aux, minss(i,n+1-j));
    1481 [ +  + ][ +  + ]:      29656 :     if (k-i < j-k-1 || (k-i == j-k-1 && 2*k > n))
                 [ +  + ]
    1482                 :      20319 :       { lrmax = lrho; j=k+1; radii[j] = lrho - aux; }
    1483                 :            :     else
    1484                 :       9337 :       { lrmin = lrho; i=k;   radii[i] = lrho + aux; }
    1485                 :            :   }
    1486                 :      24619 :   aux = lrmax - lrmin;
    1487                 :            : 
    1488         [ +  + ]:      24619 :   if (ctr)
    1489                 :            :   {
    1490                 :      22569 :     lrho = (lrmax + lrmin) / 2;
    1491         [ +  + ]:     153290 :     for (i=1; i<=n; i++)
    1492         [ +  + ]:     130721 :       if (radii[i] != UNDEF) radii[i] -= lrho;
    1493                 :            : 
    1494                 :      22569 :     bit2 = bit + (long)(n * fabs(lrho)/LOG2 + 1.);
    1495                 :      22569 :     R = mygprec(dblexp(-lrho), bit2);
    1496                 :      22569 :     q = scalepol(p,R,bit2);
    1497                 :      22569 :     conformal_mapping(radii, ctr, q, k, bit2, aux, &FF, &GG);
    1498                 :            :   }
    1499                 :            :   else
    1500                 :            :   {
    1501                 :       2050 :     lrho = logradius(radii, p, k, aux/10., &delta);
    1502                 :       2050 :     update_radius(n, radii, lrho, &param, &param2);
    1503                 :            : 
    1504                 :       2050 :     bit2 = bit + (long)(n * fabs(lrho)/LOG2 + 1.);
    1505                 :       2050 :     R = mygprec(dblexp(-lrho), bit2);
    1506                 :       2050 :     q = scalepol(p,R,bit2);
    1507                 :       2050 :     optimize_split(q, k, delta, bit2, &FF, &GG, param, param2);
    1508                 :            :   }
    1509                 :      24619 :   bit  += n;
    1510                 :      24619 :   bit2 += n; R = invr(mygprec(R,bit2));
    1511                 :      24619 :   *F = mygprec(scalepol(FF,R,bit2), bit);
    1512                 :      24619 :   *G = mygprec(scalepol(GG,R,bit2), bit);
    1513                 :      24619 : }
    1514                 :            : 
    1515                 :            : /* procedure corresponding to steps 5,6,.. page 44 in RR n. 1852 */
    1516                 :            : /* put in F and G two polynomial such that |p-FG|<2^(-bit)|p|
    1517                 :            :  * where the maximum modulus of the roots of p is <=1.
    1518                 :            :  * Assume sum of roots is 0. */
    1519                 :            : static void
    1520                 :      22569 : split_1(GEN p, long bit, GEN *F, GEN *G)
    1521                 :            : {
    1522                 :      22569 :   long i, imax, n = degpol(p), polreal = isreal(p), ep = gexpo(p), bit2 = bit+n;
    1523                 :            :   GEN TWO, ctr, q, qq, FF, GG, v, gr, r, newq;
    1524                 :            :   double lrmin, lrmax, lthick;
    1525                 :      22569 :   const double LOG3 = 1.098613;
    1526                 :            : 
    1527                 :      22569 :   lrmax = logmax_modulus(p, 0.01);
    1528                 :      22569 :   gr = mygprec(dblexp(-lrmax), bit2);
    1529                 :      22569 :   q = scalepol(p,gr,bit2);
    1530                 :            : 
    1531                 :      22569 :   bit2 = bit + gexpo(q) - ep + (long)((double)n*2.*log2(3.)+1);
    1532                 :      22569 :   TWO = myreal_1(bit2); setexpo(TWO,1);
    1533                 :      22569 :   v = cgetg(5,t_VEC);
    1534                 :      22569 :   gel(v,1) = TWO;
    1535                 :      22569 :   gel(v,2) = negr(TWO);
    1536                 :      22569 :   gel(v,3) = mkcomplex(gen_0, gel(v,1));
    1537                 :      22569 :   gel(v,4) = mkcomplex(gen_0, gel(v,2));
    1538                 :      22569 :   q = mygprec(q,bit2); lthick = 0;
    1539                 :      22569 :   newq = ctr = NULL; /* -Wall */
    1540         [ +  + ]:      22569 :   imax = polreal? 3: 4;
    1541         [ +  + ]:      41549 :   for (i=1; i<=imax; i++)
    1542                 :            :   {
    1543                 :      41038 :     qq = RgX_translate(q, gel(v,i));
    1544                 :      41038 :     lrmin = logmin_modulus(qq,0.05);
    1545         [ +  + ]:      41038 :     if (LOG3 > lrmin + lthick)
    1546                 :            :     {
    1547                 :      40035 :       double lquo = logmax_modulus(qq,0.05) - lrmin;
    1548         [ +  + ]:      40035 :       if (lquo > lthick) { lthick = lquo; newq = qq; ctr = gel(v,i); }
    1549                 :            :     }
    1550         [ +  + ]:      41038 :     if (lthick > LOG2) break;
    1551 [ +  + ][ +  + ]:      22381 :     if (polreal && i==2 && lthick > LOG3 - LOG2) break;
                 [ +  + ]
    1552                 :            :   }
    1553                 :      22569 :   bit2 = bit + gexpo(newq) - ep + (long)(n*LOG3/LOG2 + 1);
    1554                 :      22569 :   split_2(newq, bit2, ctr, lthick, &FF, &GG);
    1555                 :      22569 :   r = gneg(mygprec(ctr,bit2));
    1556                 :      22569 :   FF = RgX_translate(FF,r);
    1557                 :      22569 :   GG = RgX_translate(GG,r);
    1558                 :            : 
    1559                 :      22569 :   gr = invr(gr); bit2 = bit - ep + gexpo(FF)+gexpo(GG);
    1560                 :      22569 :   *F = scalepol(FF,gr,bit2);
    1561                 :      22569 :   *G = scalepol(GG,gr,bit2);
    1562                 :      22569 : }
    1563                 :            : 
    1564                 :            : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|,
    1565                 :            : where the maximum modulus of the roots of p is < 0.5 */
    1566                 :            : static int
    1567                 :      22569 : split_0_2(GEN p, long bit, GEN *F, GEN *G)
    1568                 :            : {
    1569                 :            :   GEN q, b, FF, GG;
    1570                 :      22569 :   long n = degpol(p), k, bit2, eq;
    1571                 :      22569 :   double aux = dbllog2(gel(p,n+1)) - dbllog2(gel(p,n+2));
    1572                 :            : 
    1573                 :            :   /* beware double overflow */
    1574 [ +  + ][ +  - ]:      22569 :   if (aux >= 0 && (aux > 1e4 || exp2(aux) > 2.5*n)) return 0;
                 [ -  + ]
    1575                 :            : 
    1576         [ +  + ]:      22569 :   aux = (aux < -300)? 0.: n*log2(1 + exp2(aux)/(double)n);
    1577                 :      22569 :   bit2 = bit+1 + (long)(log2((double)n) + aux);
    1578                 :            : 
    1579                 :      22569 :   q = mygprec(p,bit2);
    1580                 :      22569 :   b = gdivgs(gdiv(gel(q,n+1),gel(q,n+2)),-n);
    1581                 :      22569 :   q = RgX_translate(q,b); gel(q,n+1) = gen_0; eq = gexpo(q);
    1582                 :      22569 :   k = 0;
    1583 [ +  - ][ -  + ]:      22569 :   while (k <= n/2 && (- gexpo(gel(q,k+2)) > bit2 + 2*(n-k) + eq
    1584         [ -  + ]:      22569 :                       || gequal0(gel(q,k+2)))) k++;
    1585         [ -  + ]:      22569 :   if (k > 0)
    1586                 :            :   {
    1587         [ #  # ]:          0 :     if (k > n/2) k = n/2;
    1588                 :          0 :     bit2 += k<<1;
    1589                 :          0 :     FF = monomial(myreal_1(bit2), k, 0);
    1590                 :          0 :     GG = RgX_shift_shallow(q, -k);
    1591                 :            :   }
    1592                 :            :   else
    1593                 :            :   {
    1594                 :      22569 :     split_1(q,bit2,&FF,&GG);
    1595                 :      22569 :     bit2 = bit + gexpo(FF) + gexpo(GG) - gexpo(p) + (long)aux+1;
    1596                 :      22569 :     FF = mygprec(FF,bit2);
    1597                 :            :   }
    1598                 :      22569 :   GG = mygprec(GG,bit2); b = mygprec(gneg(b),bit2);
    1599                 :      22569 :   *F = RgX_translate(FF, b);
    1600                 :      22569 :   *G = RgX_translate(GG, b); return 1;
    1601                 :            : }
    1602                 :            : 
    1603                 :            : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|.
    1604                 :            :  * Assume max_modulus(p) < 2 */
    1605                 :            : static void
    1606                 :      22569 : split_0_1(GEN p, long bit, GEN *F, GEN *G)
    1607                 :            : {
    1608                 :            :   GEN FF, GG;
    1609                 :            :   long n, bit2, normp;
    1610                 :            : 
    1611         [ -  + ]:      45138 :   if  (split_0_2(p,bit,F,G)) return;
    1612                 :            : 
    1613                 :          0 :   normp = gexpo(p);
    1614                 :          0 :   scalepol2n(p,2); /* p := 4^(-n) p(4*x) */
    1615                 :          0 :   n = degpol(p); bit2 = bit + 2*n + gexpo(p) - normp;
    1616                 :          0 :   split_1(mygprec(p,bit2), bit2,&FF,&GG);
    1617                 :          0 :   scalepol2n(FF,-2);
    1618                 :          0 :   scalepol2n(GG,-2); bit2 = bit + gexpo(FF) + gexpo(GG) - normp;
    1619                 :          0 :   *F = mygprec(FF,bit2);
    1620                 :      22569 :   *G = mygprec(GG,bit2);
    1621                 :            : }
    1622                 :            : 
    1623                 :            : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P| */
    1624                 :            : static void
    1625                 :      24619 : split_0(GEN p, long bit, GEN *F, GEN *G)
    1626                 :            : {
    1627                 :      24619 :   const double LOG1_9 = 0.6418539;
    1628                 :      24619 :   long n = degpol(p), k = 0;
    1629                 :            :   GEN q;
    1630                 :            : 
    1631 [ -  + ][ #  # ]:      24619 :   while (gexpo(gel(p,k+2)) < -bit && k <= n/2) k++;
    1632         [ -  + ]:      24619 :   if (k > 0)
    1633                 :            :   {
    1634         [ #  # ]:          0 :     if (k > n/2) k = n/2;
    1635                 :          0 :     *F = monomial(myreal_1(bit), k, 0);
    1636                 :          0 :     *G = RgX_shift_shallow(p, -k);
    1637                 :            :   }
    1638                 :            :   else
    1639                 :            :   {
    1640                 :      24619 :     double lr = logmax_modulus(p, 0.05);
    1641         [ +  + ]:      24619 :     if (lr < LOG1_9) split_0_1(p, bit, F, G);
    1642                 :            :     else
    1643                 :            :     {
    1644                 :      19768 :       q = RgX_recip_shallow(p);
    1645                 :      19768 :       lr = logmax_modulus(q,0.05);
    1646         [ +  + ]:      19768 :       if (lr < LOG1_9)
    1647                 :            :       {
    1648                 :      17718 :         split_0_1(q, bit, F, G);
    1649                 :      17718 :         *F = RgX_recip_shallow(*F);
    1650                 :      17718 :         *G = RgX_recip_shallow(*G);
    1651                 :            :       }
    1652                 :            :       else
    1653                 :       2050 :         split_2(p,bit,NULL, 1.2837,F,G);
    1654                 :            :     }
    1655                 :            :   }
    1656                 :      24619 : }
    1657                 :            : 
    1658                 :            : /********************************************************************/
    1659                 :            : /**                                                                **/
    1660                 :            : /**                ERROR ESTIMATE FOR THE ROOTS                    **/
    1661                 :            : /**                                                                **/
    1662                 :            : /********************************************************************/
    1663                 :            : 
    1664                 :            : static GEN
    1665                 :      78346 : root_error(long n, long k, GEN roots_pol, long err, GEN shatzle)
    1666                 :            : {
    1667                 :      78346 :   GEN rho, d, eps, epsbis, eps2, aux, rap = NULL;
    1668                 :            :   long i, j;
    1669                 :            : 
    1670                 :      78346 :   d = cgetg(n+1,t_VEC);
    1671         [ +  + ]:    1245294 :   for (i=1; i<=n; i++)
    1672                 :            :   {
    1673         [ +  + ]:    1166948 :     if (i!=k)
    1674                 :            :     {
    1675                 :    1088602 :       aux = gsub(gel(roots_pol,i), gel(roots_pol,k));
    1676                 :    1088602 :       gel(d,i) = gabs(mygprec(aux,31), DEFAULTPREC);
    1677                 :            :     }
    1678                 :            :   }
    1679                 :      78346 :   rho = gabs(mygprec(gel(roots_pol,k),31), DEFAULTPREC);
    1680         [ +  + ]:      78346 :   if (expo(rho) < 0) rho = real_1(DEFAULTPREC);
    1681                 :      78346 :   eps = mulrr(rho, shatzle);
    1682                 :      78346 :   aux = shiftr(powru(rho,n), err);
    1683                 :            : 
    1684 [ +  + ][ +  - ]:     247549 :   for (j=1; j<=2 || (j<=5 && cmprr(rap, dbltor(1.2)) > 0); j++)
                 [ +  + ]
    1685                 :            :   {
    1686                 :     169203 :     GEN prod = NULL; /* 1. */
    1687                 :     169203 :     long m = n;
    1688                 :     169203 :     epsbis = mulrr(eps, dbltor(1.25));
    1689         [ +  + ]:    3067733 :     for (i=1; i<=n; i++)
    1690                 :            :     {
    1691 [ +  + ][ +  + ]:    2898530 :       if (i != k && cmprr(gel(d,i),epsbis) > 0)
    1692                 :            :       {
    1693                 :    2704655 :         GEN dif = subrr(gel(d,i),eps);
    1694         [ +  + ]:    2704655 :         prod = prod? mulrr(prod, dif): dif;
    1695                 :    2704655 :         m--;
    1696                 :            :       }
    1697                 :            :     }
    1698         [ +  + ]:     169203 :     eps2 = prod? divrr(aux, prod): aux;
    1699         [ +  + ]:     169203 :     if (m > 1) eps2 = sqrtnr(shiftr(eps2, 2*m-2), m);
    1700                 :     169203 :     rap = divrr(eps,eps2); eps = eps2;
    1701                 :            :   }
    1702                 :      78346 :   return eps;
    1703                 :            : }
    1704                 :            : 
    1705                 :            : /* round a complex or real number x to an absolute value of 2^(-bit) */
    1706                 :            : static GEN
    1707                 :     181956 : mygprec_absolute(GEN x, long bit)
    1708                 :            : {
    1709                 :            :   long e;
    1710                 :            :   GEN y;
    1711                 :            : 
    1712      [ +  +  + ]:     181956 :   switch(typ(x))
    1713                 :            :   {
    1714                 :            :     case t_REAL:
    1715                 :     123744 :       e = expo(x) + bit;
    1716 [ +  + ][ -  + ]:     123744 :       return (e <= 0 || !signe(x))? real_0_bit(-bit): rtor(x, nbits2prec(e));
    1717                 :            :     case t_COMPLEX:
    1718         [ +  + ]:      52512 :       if (gexpo(gel(x,2)) < -bit) return mygprec_absolute(gel(x,1),bit);
    1719                 :      51098 :       y = cgetg(3,t_COMPLEX);
    1720                 :      51098 :       gel(y,1) = mygprec_absolute(gel(x,1),bit);
    1721                 :      51098 :       gel(y,2) = mygprec_absolute(gel(x,2),bit);
    1722                 :      51098 :       return y;
    1723                 :     181956 :     default: return x;
    1724                 :            :   }
    1725                 :            : }
    1726                 :            : 
    1727                 :            : static long
    1728                 :      11736 : a_posteriori_errors(GEN p, GEN roots_pol, long err)
    1729                 :            : {
    1730                 :      11736 :   long i, n = degpol(p), e_max = -(long)EXPOBITS;
    1731                 :            :   GEN sigma, shatzle;
    1732                 :            : 
    1733                 :      11736 :   err += (long)log2((double)n) + 1;
    1734         [ -  + ]:      11736 :   if (err > -2) return 0;
    1735                 :      11736 :   sigma = real2n(-err, LOWDEFAULTPREC);
    1736                 :            :   /*  2 / ((s - 1)^(1/n) - 1) */
    1737                 :      11736 :   shatzle = divur(2, subrs(sqrtnr(subrs(sigma,1),n), 1));
    1738         [ +  + ]:      90082 :   for (i=1; i<=n; i++)
    1739                 :            :   {
    1740                 :      78346 :     pari_sp av = avma;
    1741                 :      78346 :     GEN x = root_error(n,i,roots_pol,err,shatzle);
    1742                 :      78346 :     long e = gexpo(x);
    1743         [ +  + ]:      78346 :     avma = av; if (e > e_max) e_max = e;
    1744                 :      78346 :     gel(roots_pol,i) = mygprec_absolute(gel(roots_pol,i), -e);
    1745                 :            :   }
    1746                 :      11736 :   return e_max;
    1747                 :            : }
    1748                 :            : 
    1749                 :            : /********************************************************************/
    1750                 :            : /**                                                                **/
    1751                 :            : /**                           MAIN                                 **/
    1752                 :            : /**                                                                **/
    1753                 :            : /********************************************************************/
    1754                 :            : static GEN
    1755                 :      58231 : append_clone(GEN r, GEN a) { a = gclone(a); vectrunc_append(r, a); return a; }
    1756                 :            : 
    1757                 :            : /* put roots in placeholder roots_pol so that |P - L_1...L_n| < 2^(-bit)|P|
    1758                 :            :  * returns prod (x-roots_pol[i]) */
    1759                 :            : static GEN
    1760                 :      60974 : split_complete(GEN p, long bit, GEN roots_pol)
    1761                 :            : {
    1762                 :      60974 :   long n = degpol(p);
    1763                 :            :   pari_sp ltop;
    1764                 :            :   GEN p1, F, G, a, b, m1, m2;
    1765                 :            : 
    1766         [ +  + ]:      60974 :   if (n == 1)
    1767                 :            :   {
    1768                 :      14479 :     a = gneg_i(gdiv(gel(p,2), gel(p,3)));
    1769                 :      14479 :     (void)append_clone(roots_pol,a); return p;
    1770                 :            :   }
    1771                 :      46495 :   ltop = avma;
    1772         [ +  + ]:      46495 :   if (n == 2)
    1773                 :            :   {
    1774                 :      21876 :     F = gsub(gsqr(gel(p,3)), gmul2n(gmul(gel(p,2),gel(p,4)), 2));
    1775                 :      21876 :     F = gsqrt(F, nbits2prec(bit));
    1776                 :      21876 :     p1 = ginv(gmul2n(gel(p,4),1));
    1777                 :      21876 :     a = gneg_i(gmul(gadd(F,gel(p,3)), p1));
    1778                 :      21876 :     b =        gmul(gsub(F,gel(p,3)), p1);
    1779                 :      21876 :     a = append_clone(roots_pol,a);
    1780                 :      21876 :     b = append_clone(roots_pol,b); avma = ltop;
    1781                 :      21876 :     a = mygprec(a, 3*bit);
    1782                 :      21876 :     b = mygprec(b, 3*bit);
    1783                 :      21876 :     return gmul(gel(p,4), mkpoln(3, gen_1, gneg(gadd(a,b)), gmul(a,b)));
    1784                 :            :   }
    1785                 :      24619 :   split_0(p,bit,&F,&G);
    1786                 :      24619 :   m1 = split_complete(F,bit,roots_pol);
    1787                 :      24619 :   m2 = split_complete(G,bit,roots_pol);
    1788                 :      60974 :   return gerepileupto(ltop, gmul(m1,m2));
    1789                 :            : }
    1790                 :            : 
    1791                 :            : static GEN
    1792                 :     353747 : quicktofp(GEN x)
    1793                 :            : {
    1794                 :     353747 :   const long prec = DEFAULTPREC;
    1795   [ +  +  -  +  :     353747 :   switch(typ(x))
                      - ]
    1796                 :            :   {
    1797                 :     350534 :     case t_INT: return itor(x, prec);
    1798                 :       3192 :     case t_REAL: return rtor(x, prec);
    1799                 :          0 :     case t_FRAC: return fractor(x, prec);
    1800                 :            :     case t_COMPLEX: {
    1801                 :         21 :       GEN a = gel(x,1), b = gel(x,2);
    1802                 :            :       /* avoid problem with 0, e.g. x = 0 + I*1e-100. We don't want |x| = 0. */
    1803         [ -  + ]:         21 :       if (isintzero(a)) return cxcompotor(b, prec);
    1804         [ -  + ]:         21 :       if (isintzero(b)) return cxcompotor(a, prec);
    1805                 :         21 :       a = cxcompotor(a, prec);
    1806                 :         21 :       b = cxcompotor(b, prec); return sqrtr(addrr(sqrr(a), sqrr(b)));
    1807                 :            :     }
    1808                 :          0 :     default: pari_err_TYPE("quicktofp",x);
    1809                 :     353747 :       return NULL;/*not reached*/
    1810                 :            :   }
    1811                 :            : 
    1812                 :            : }
    1813                 :            : 
    1814                 :            : /* bound log_2 |largest root of p| (Fujiwara's bound) */
    1815                 :            : double
    1816                 :      71958 : fujiwara_bound(GEN p)
    1817                 :            : {
    1818                 :      71958 :   pari_sp av = avma;
    1819                 :      71958 :   long i, n = degpol(p);
    1820                 :            :   GEN cc;
    1821                 :            :   double loglc, Lmax;
    1822                 :            : 
    1823         [ -  + ]:      71958 :   if (n <= 0) pari_err_CONSTPOL("fujiwara_bound");
    1824                 :      71958 :   loglc = mydbllog2r( quicktofp(gel(p,n+2)) ); /* log_2 |lc(p)| */
    1825                 :      71958 :   cc = gel(p, 2);
    1826         [ +  + ]:      71958 :   if (gequal0(cc))
    1827                 :       5759 :     Lmax = -pariINFINITY-1;
    1828                 :            :   else
    1829                 :      66199 :     Lmax = (mydbllog2r(quicktofp(cc)) - loglc - 1) / n;
    1830         [ +  + ]:     338163 :   for (i = 1; i < n; i++)
    1831                 :            :   {
    1832                 :     266205 :     GEN y = gel(p,i+2);
    1833                 :            :     double L;
    1834         [ +  + ]:     266205 :     if (gequal0(y)) continue;
    1835                 :     215590 :     L = (mydbllog2r(quicktofp(y)) - loglc) / (n-i);
    1836         [ +  + ]:     215590 :     if (L > Lmax) Lmax = L;
    1837                 :            :   }
    1838                 :      71958 :   avma = av; return Lmax + 1;
    1839                 :            : }
    1840                 :            : 
    1841                 :            : /* Fujiwara's bound, real roots. Based on the following remark: if
    1842                 :            :  *   p = x^n + sum a_i x^i and q = x^n + sum min(a_i,0)x^i
    1843                 :            :  * then for all x >= 0, p(x) >= q(x). Thus any bound for the (positive) roots
    1844                 :            :  * of q is a bound for the positive roots of p. */
    1845                 :            : double
    1846                 :       8188 : fujiwara_bound_real(GEN p, long sign)
    1847                 :            : {
    1848                 :       8188 :   pari_sp av = avma;
    1849                 :            :   GEN x;
    1850                 :       8188 :   long n = degpol(p), i, signodd, signeven;
    1851                 :            :   double fb;
    1852         [ -  + ]:       8188 :   if (n <= 0) pari_err_CONSTPOL("fujiwara_bound");
    1853                 :       8188 :   x = shallowcopy(p);
    1854         [ +  - ]:       8188 :   if (gsigne(gel(x, n+2)) > 0)
    1855                 :            :   {
    1856                 :       8188 :     signeven = 1;
    1857                 :       8188 :     signodd = sign;
    1858                 :            :   }
    1859                 :            :   else
    1860                 :            :   {
    1861                 :          0 :     signeven = -1;
    1862                 :          0 :     signodd = -sign;
    1863                 :            :   }
    1864         [ +  + ]:      50215 :   for (i = 0; i < n; i++)
    1865                 :            :   {
    1866         [ +  + ]:      42027 :     if ((n - i) % 2)
    1867                 :            :     {
    1868         [ +  + ]:      22541 :       if (gsigne(gel(x, i+2)) == signodd ) gel(x, i+2) = gen_0;
    1869                 :            :     }
    1870                 :            :     else
    1871                 :            :     {
    1872         [ +  + ]:      19486 :       if (gsigne(gel(x, i+2)) == signeven) gel(x, i+2) = gen_0;
    1873                 :            :     }
    1874                 :            :   }
    1875                 :       8188 :   fb = fujiwara_bound(x);
    1876                 :       8188 :   avma = av; return fb;
    1877                 :            : }
    1878                 :            : 
    1879                 :            : static GEN
    1880                 :      90124 : mygprecrc_special(GEN x, long prec, long e)
    1881                 :            : {
    1882                 :            :   GEN y;
    1883      [ +  +  + ]:      90124 :   switch(typ(x))
    1884                 :            :   {
    1885                 :            :     case t_REAL:
    1886         [ -  + ]:       3225 :       if (!signe(x)) return real_0_bit(minss(e, expo(x)));
    1887         [ +  - ]:       3225 :       return (prec > realprec(x))? rtor(x, prec): x;
    1888                 :            :     case t_COMPLEX:
    1889                 :         21 :       y = cgetg(3,t_COMPLEX);
    1890                 :         21 :       gel(y,1) = mygprecrc_special(gel(x,1),prec,e);
    1891                 :         21 :       gel(y,2) = mygprecrc_special(gel(x,2),prec,e);
    1892                 :         21 :       return y;
    1893                 :      90124 :     default: return x;
    1894                 :            :   }
    1895                 :            : }
    1896                 :            : 
    1897                 :            : /* like mygprec but keep at least the same precision as before */
    1898                 :            : static GEN
    1899                 :      11736 : mygprec_special(GEN x, long bit)
    1900                 :            : {
    1901                 :            :   long lx, i, e, prec;
    1902                 :            :   GEN y;
    1903                 :            : 
    1904         [ -  + ]:      11736 :   if (bit < 0) bit = 0; /* should not happen */
    1905                 :      11736 :   e = gexpo(x) - bit;
    1906                 :      11736 :   prec = nbits2prec(bit);
    1907         [ +  - ]:      11736 :   switch(typ(x))
    1908                 :            :   {
    1909                 :            :     case t_POL:
    1910                 :      11736 :       y = cgetg_copy(x, &lx); y[1] = x[1];
    1911         [ +  + ]:     101818 :       for (i=2; i<lx; i++) gel(y,i) = mygprecrc_special(gel(x,i),prec,e);
    1912                 :      11736 :       break;
    1913                 :            : 
    1914                 :          0 :     default: y = mygprecrc_special(x,prec,e);
    1915                 :            :   }
    1916                 :      11736 :   return y;
    1917                 :            : }
    1918                 :            : 
    1919                 :            : static GEN
    1920                 :       6544 : fix_roots1(GEN r)
    1921                 :            : {
    1922                 :       6544 :   long i, l = lg(r);
    1923                 :       6544 :   GEN allr = cgetg(l, t_VEC);
    1924         [ +  + ]:      51297 :   for (i=1; i<l; i++)
    1925                 :            :   {
    1926                 :      44753 :     GEN t = gel(r,i);
    1927                 :      44753 :     gel(allr,i) = gcopy(t); gunclone(t);
    1928                 :            :   }
    1929                 :       6544 :   return allr;
    1930                 :            : }
    1931                 :            : static GEN
    1932                 :      11736 : fix_roots(GEN r, GEN *m, long h, long bit)
    1933                 :            : {
    1934                 :            :   long i, j, k, l, prec;
    1935                 :            :   GEN allr, ro1;
    1936         [ +  + ]:      11736 :   if (h == 1) return fix_roots1(r);
    1937                 :       5192 :   ro1 = initRUgen(h, bit);
    1938                 :       5192 :   prec = nbits2prec(bit);
    1939                 :       5192 :   l = lg(r)-1;
    1940                 :       5192 :   allr = cgetg(h*l+1, t_VEC);
    1941         [ +  + ]:      18670 :   for (k=1,i=1; i<=l; i++)
    1942                 :            :   {
    1943                 :      13478 :     GEN p2, p1 = gel(r,i);
    1944         [ +  + ]:      13478 :     p2 = (h == 2)? gsqrt(p1, prec): gsqrtn(p1, utoipos(h), NULL, prec);
    1945         [ +  + ]:      47071 :     for (j=0; j<h; j++) gel(allr,k++) = gmul(p2, gel(ro1,j));
    1946                 :      13478 :     gunclone(p1);
    1947                 :            :   }
    1948                 :       5192 :   *m = roots_to_pol(allr, 0);
    1949                 :      11736 :   return allr;
    1950                 :            : }
    1951                 :            : 
    1952                 :            : static GEN
    1953                 :        357 : RgX_normalize1(GEN x)
    1954                 :            : {
    1955                 :        357 :   long i, n = lg(x)-1;
    1956                 :            :   GEN y;
    1957         [ +  - ]:        364 :   for (i = n; i > 1; i--)
    1958         [ +  + ]:        364 :     if (!gequal0( gel(x,i) )) break;
    1959         [ +  + ]:        357 :   if (i == n) return x;
    1960                 :          7 :   pari_warn(warner,"normalizing a polynomial with 0 leading term");
    1961         [ -  + ]:          7 :   if (i == 1) pari_err_ROOTS0("roots");
    1962                 :          7 :   y = cgetg(i+1, t_POL); y[1] = x[1];
    1963         [ +  + ]:         14 :   for (; i > 1; i--) gel(y,i) = gel(x,i);
    1964                 :        357 :   return y;
    1965                 :            : }
    1966                 :            : 
    1967                 :            : static GEN
    1968                 :      11570 : all_roots(GEN p, long bit)
    1969                 :            : {
    1970                 :            :   GEN lc, pd, q, roots_pol, m;
    1971                 :      11570 :   long bit0,  bit2, i, e, h, n = degpol(p);
    1972                 :            :   pari_sp av;
    1973                 :            : 
    1974                 :      11570 :   pd = RgX_deflate_max(p, &h); lc = leading_term(pd);
    1975         [ +  + ]:      11570 :   e = (long)(2 * fujiwara_bound(pd)); if (e < 0) e = 0;
    1976                 :      11570 :   bit0 = bit + gexpo(pd) - gexpo(lc) + (long)log2(n/h)+1+e;
    1977                 :      11570 :   bit2 = bit0; e = 0;
    1978                 :      11570 :   for (av=avma,i=1;; i++,avma=av)
    1979                 :            :   {
    1980                 :      11736 :     roots_pol = vectrunc_init(n+1);
    1981                 :      11736 :     bit2 += e + (n << i);
    1982                 :      11736 :     q = RgX_gtofp_bit(mygprec(pd,bit2), bit2);
    1983                 :      11736 :     q[1] = evalsigne(1)|evalvarn(0);
    1984                 :      11736 :     m = split_complete(q,bit2,roots_pol);
    1985                 :      11736 :     roots_pol = fix_roots(roots_pol, &m, h, bit2);
    1986                 :      11736 :     q = mygprec_special(p,bit2); lc = leading_term(q);
    1987                 :      11736 :     q[1] = evalsigne(1)|evalvarn(0);
    1988         [ +  + ]:      11736 :     if (h > 1) m = gmul(m,lc);
    1989                 :            : 
    1990                 :      11736 :     e = gexpo(gsub(q, m)) - gexpo(lc) + (long)log2((double)n) + 1;
    1991         [ -  + ]:      11736 :     if (e < -2*bit2) e = -2*bit2; /* avoid e = -pariINFINITY */
    1992         [ +  - ]:      11736 :     if (e < 0)
    1993                 :            :     {
    1994                 :      11736 :       e = bit + a_posteriori_errors(p,roots_pol,e);
    1995         [ +  + ]:      11736 :       if (e < 0) return roots_pol;
    1996                 :            :     }
    1997         [ -  + ]:        166 :     if (DEBUGLEVEL > 7)
    1998                 :          0 :       err_printf("all_roots: restarting, i = %ld, e = %ld\n", i,e);
    1999                 :        166 :   }
    2000                 :            : }
    2001                 :            : 
    2002                 :            : INLINE int
    2003                 :       2370 : isexactscalar(GEN x) { long tx = typ(x); return is_rational_t(tx); }
    2004                 :            : 
    2005                 :            : static int
    2006                 :        848 : isexactpol(GEN p)
    2007                 :            : {
    2008                 :        848 :   long i,n = degpol(p);
    2009         [ +  + ]:       2861 :   for (i=0; i<=n; i++)
    2010         [ +  + ]:       2370 :     if (!isexactscalar(gel(p,i+2))) return 0;
    2011                 :        848 :   return 1;
    2012                 :            : }
    2013                 :            : 
    2014                 :            : static long
    2015                 :       5674 : isvalidcoeff(GEN x)
    2016                 :            : {
    2017      [ +  +  + ]:       5674 :   switch (typ(x))
    2018                 :            :   {
    2019                 :       5618 :     case t_INT: case t_REAL: case t_FRAC: return 1;
    2020 [ +  - ][ +  - ]:         42 :     case t_COMPLEX: return isvalidcoeff(gel(x,1)) && isvalidcoeff(gel(x,2));
    2021                 :            :   }
    2022                 :       5674 :   return 0;
    2023                 :            : }
    2024                 :            : 
    2025                 :            : static void
    2026                 :        855 : checkvalidpol(GEN p)
    2027                 :            : {
    2028                 :        855 :   long i,n = lg(p);
    2029         [ +  + ]:       6424 :   for (i=2; i<n; i++)
    2030         [ +  + ]:       5576 :     if (!isvalidcoeff(gel(p,i))) pari_err_TYPE("roots", gel(p,i));
    2031                 :        848 : }
    2032                 :            : 
    2033                 :            : static GEN
    2034                 :        491 : solve_exact_pol(GEN p, long bit)
    2035                 :            : {
    2036                 :        491 :   long i, j, k, m, n = degpol(p), iroots = 0;
    2037                 :        491 :   GEN ex, factors, v = zerovec(n);
    2038                 :            : 
    2039                 :        491 :   factors = ZX_squff(Q_primpart(p), &ex);
    2040         [ +  + ]:        982 :   for (i=1; i<lg(factors); i++)
    2041                 :            :   {
    2042                 :        491 :     GEN roots_fact = all_roots(gel(factors,i), bit);
    2043                 :        491 :     n = degpol(gel(factors,i));
    2044                 :        491 :     m = ex[i];
    2045         [ +  + ]:       1971 :     for (j=1; j<=n; j++)
    2046         [ +  + ]:       2960 :       for (k=1; k<=m; k++) v[++iroots] = roots_fact[j];
    2047                 :            :   }
    2048                 :        491 :   return v;
    2049                 :            : }
    2050                 :            : 
    2051                 :            : /* return the roots of p with absolute error bit */
    2052                 :            : static GEN
    2053                 :        848 : roots_com(GEN q, long bit)
    2054                 :            : {
    2055                 :            :   GEN L, p;
    2056                 :        848 :   long v = RgX_valrem_inexact(q, &p);
    2057                 :        848 :   int ex = isexactpol(p);
    2058         [ +  + ]:        848 :   if (!ex) p = RgX_normalize1(p);
    2059         [ +  + ]:        848 :   if (lg(p) == 3)
    2060                 :          7 :     L = cgetg(1,t_VEC); /* constant polynomial */
    2061                 :            :   else
    2062         [ +  + ]:        841 :     L = ex? solve_exact_pol(p,bit): all_roots(p,bit);
    2063         [ +  + ]:        848 :   if (v)
    2064                 :            :   {
    2065                 :         84 :     GEN M, z, t = gel(q,2);
    2066                 :            :     long i, x, y, l, n;
    2067                 :            : 
    2068         [ +  + ]:         84 :     if (isrationalzero(t)) x = -bit;
    2069                 :            :     else
    2070                 :            :     {
    2071                 :         14 :       n = gexpo(t);
    2072                 :         14 :       x = n / v; l = degpol(q);
    2073         [ +  + ]:         56 :       for (i = v; i <= l; i++)
    2074                 :            :       {
    2075                 :         42 :         t  = gel(q,i+2);
    2076         [ -  + ]:         42 :         if (isrationalzero(t)) continue;
    2077                 :         42 :         y = (n - gexpo(t)) / i;
    2078         [ +  + ]:         42 :         if (y < x) x = y;
    2079                 :            :       }
    2080                 :            :     }
    2081                 :         84 :     z = real_0_bit(x); l = v + lg(L);
    2082                 :         84 :     M = cgetg(l, t_VEC); L -= v;
    2083         [ +  + ]:        168 :     for (i = 1; i <= v; i++) gel(M,i) = z;
    2084         [ +  + ]:        245 :     for (     ; i <  l; i++) gel(M,i) = gel(L,i);
    2085                 :         84 :     L = M;
    2086                 :            :   }
    2087                 :        848 :   return L;
    2088                 :            : }
    2089                 :            : 
    2090                 :            : static GEN
    2091                 :      53995 : tocomplex(GEN x, long l, long bit)
    2092                 :            : {
    2093                 :            :   GEN y;
    2094         [ +  + ]:      53995 :   if (typ(x) == t_COMPLEX)
    2095                 :            :   {
    2096         [ +  + ]:      50478 :     if (signe(gel(x,1))) return mygprecrc(x, l, -bit);
    2097                 :       5344 :     x = gel(x,2);
    2098                 :       5344 :     y = cgetg(3,t_COMPLEX);
    2099                 :       5344 :     gel(y,1) = real_0_bit(-bit);
    2100                 :       5344 :     gel(y,2) = mygprecrc(x, l, -bit);
    2101                 :            :   }
    2102                 :            :   else
    2103                 :            :   {
    2104                 :       3517 :     y = cgetg(3,t_COMPLEX);
    2105                 :       3517 :     gel(y,1) = mygprecrc(x, l, -bit);
    2106                 :       3517 :     gel(y,2) = real_0_bit(-bit);
    2107                 :            :   }
    2108                 :      53995 :   return y;
    2109                 :            : }
    2110                 :            : 
    2111                 :            : /* x,y are t_COMPLEX of t_REALs or t_REAL, compare lexicographically,
    2112                 :            :  * up to 2^-e absolute error */
    2113                 :            : static int
    2114                 :     157083 : cmp_complex_appr(void *E, GEN x, GEN y)
    2115                 :            : {
    2116                 :     157083 :   long e = (long)E;
    2117                 :            :   GEN z, xi, yi, xr, yr;
    2118                 :            :   long sxi, syi;
    2119         [ +  + ]:     157083 :   if (typ(x) == t_COMPLEX) { xr = gel(x,1); xi = gel(x,2); sxi = signe(xi); }
    2120                 :      58789 :   else { xr = x; xi = NULL; sxi = 0; }
    2121         [ +  + ]:     157083 :   if (typ(y) == t_COMPLEX) { yr = gel(y,1); yi = gel(y,2); syi = signe(yi); }
    2122                 :      52132 :   else { yr = y; yi = NULL; syi = 0; }
    2123                 :            :   /* Compare absolute values of imaginary parts */
    2124         [ +  + ]:     157083 :   if (!sxi)
    2125                 :            :   {
    2126 [ +  + ][ +  - ]:      65484 :     if (syi && expo(yi) >= e) return -1;
    2127                 :            :     /* |Im x| ~ |Im y| ~ 0 */
    2128                 :            :   }
    2129         [ +  + ]:      91599 :   else if (!syi)
    2130                 :            :   {
    2131 [ +  - ][ +  - ]:       2545 :     if (sxi && expo(xi) >= e) return 1;
    2132                 :            :     /* |Im x| ~ |Im y| ~ 0 */
    2133                 :            :   }
    2134                 :            :   else
    2135                 :            :   {
    2136                 :            :     long sz;
    2137                 :      89054 :     z = addrr_sign(xi, 1, yi, -1);
    2138                 :      89054 :     sz = signe(z);
    2139 [ +  + ][ +  - ]:      89054 :     if (sz && expo(z) >= e) return (int)sz;
    2140                 :            :   }
    2141                 :            :   /* |Im x| ~ |Im y|, sort according to real parts */
    2142                 :      93624 :   z = subrr(xr, yr);
    2143         [ +  + ]:      93624 :   if (expo(z) >= e) return (int)signe(z);
    2144                 :            :   /* Re x ~ Re y. Place negative absolute value before positive */
    2145                 :     157083 :   return (int) (sxi - syi);
    2146                 :            : }
    2147                 :            : 
    2148                 :            : static GEN
    2149                 :      11577 : clean_roots(GEN L, long l, long bit, long clean)
    2150                 :            : {
    2151                 :      11577 :   long i, n = lg(L), ex = 5 - bit;
    2152                 :      11577 :   GEN res = cgetg(n,t_COL);
    2153         [ +  + ]:      88633 :   for (i=1; i<n; i++)
    2154                 :            :   {
    2155                 :      77056 :     GEN c = gel(L,i);
    2156 [ +  + ][ +  + ]:      77056 :     if (clean && isrealappr(c,ex))
    2157                 :            :     {
    2158         [ -  + ]:      23061 :       if (typ(c) == t_COMPLEX) c = gel(c,1);
    2159                 :      23061 :       c = mygprecrc(c, l, -bit);
    2160                 :            :     }
    2161                 :            :     else
    2162                 :      53995 :       c = tocomplex(c, l, bit);
    2163                 :      77056 :     gel(res,i) = c;
    2164                 :            :   }
    2165                 :      11577 :   gen_sort_inplace(res, (void*)ex, &cmp_complex_appr, NULL);
    2166                 :      11577 :   return res;
    2167                 :            : }
    2168                 :            : 
    2169                 :            : /* the vector of roots of p, with absolute error 2^(- prec2nbits(l)) */
    2170                 :            : static GEN
    2171                 :        876 : roots_aux(GEN p, long l, long clean)
    2172                 :            : {
    2173                 :        876 :   pari_sp av = avma;
    2174                 :            :   long bit;
    2175                 :            :   GEN L;
    2176                 :            : 
    2177         [ +  + ]:        876 :   if (typ(p) != t_POL)
    2178                 :            :   {
    2179         [ +  + ]:         21 :     if (gequal0(p)) pari_err_ROOTS0("roots");
    2180         [ +  + ]:         14 :     if (!isvalidcoeff(p)) pari_err_TYPE("roots",p);
    2181                 :          7 :     return cgetg(1,t_COL); /* constant polynomial */
    2182                 :            :   }
    2183         [ -  + ]:        855 :   if (!signe(p)) pari_err_ROOTS0("roots");
    2184                 :        855 :   checkvalidpol(p);
    2185         [ -  + ]:        848 :   if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
    2186         [ -  + ]:        848 :   if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
    2187                 :        848 :   bit = prec2nbits(l);
    2188                 :        848 :   L = roots_com(p, bit);
    2189                 :        855 :   return gerepileupto(av, clean_roots(L, l, bit, clean));
    2190                 :            : }
    2191                 :            : GEN
    2192                 :        729 : roots(GEN p, long l) { return roots_aux(p,l, 0); }
    2193                 :            : /* clean up roots. If root is real replace it by its real part */
    2194                 :            : GEN
    2195                 :        147 : cleanroots(GEN p, long l) { return roots_aux(p,l, 1); }
    2196                 :            : 
    2197                 :            : /* private variant of conjvec. Allow non rational coefficients, shallow
    2198                 :            :  * function. */
    2199                 :            : GEN
    2200                 :         77 : polmod_to_embed(GEN x, long prec)
    2201                 :            : {
    2202                 :         77 :   GEN v, T = gel(x,1), A = gel(x,2);
    2203                 :            :   long i, l;
    2204 [ +  - ][ -  + ]:         77 :   if (typ(A) != t_POL || varn(A) != varn(T))
    2205                 :            :   {
    2206                 :          0 :     checkvalidpol(T);
    2207                 :          0 :     return const_col(degpol(T), A);
    2208                 :            :   }
    2209                 :         77 :   v = cleanroots(T,prec); l = lg(v);
    2210         [ +  + ]:        231 :   for (i=1; i<l; i++) gel(v,i) = poleval(A,gel(v,i));
    2211                 :         77 :   return v;
    2212                 :            : }
    2213                 :            : 
    2214                 :            : GEN
    2215                 :      10729 : QX_complex_roots(GEN p, long l)
    2216                 :            : {
    2217                 :      10729 :   pari_sp av = avma;
    2218                 :            :   long bit;
    2219                 :            :   GEN L;
    2220                 :            : 
    2221         [ -  + ]:      10729 :   if (!signe(p)) pari_err_ROOTS0("QX_complex_roots");
    2222         [ -  + ]:      10729 :   if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
    2223         [ -  + ]:      10729 :   if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
    2224                 :      10729 :   bit = prec2nbits(l);
    2225                 :      10729 :   L = all_roots(Q_primpart(p), bit);
    2226                 :      10729 :   return gerepileupto(av, clean_roots(L, l, bit, 1));
    2227                 :            : }
    2228                 :            : 
    2229                 :            : /********************************************************************/
    2230                 :            : /**                                                                **/
    2231                 :            : /**                REAL ROOTS OF INTEGER POLYNOMIAL                **/
    2232                 :            : /**                                                                **/
    2233                 :            : /********************************************************************/
    2234                 :            : 
    2235                 :            : /* Count sign changes in the coefficients of (x+1)^deg(P)*P(1/(x+1))
    2236                 :            :  * The inversion is implicit (we take coefficients backwards) */
    2237                 :            : static long
    2238                 :      57409 : X2XP1(GEN P, long deg, GEN *Premapped)
    2239                 :            : {
    2240                 :      57409 :   pari_sp av = avma;
    2241                 :      57409 :   GEN v = shallowcopy(P);
    2242                 :            :   long i, j, vlim, nb, s, s2;
    2243                 :            :   char flag;
    2244                 :            : 
    2245                 :      57409 :   vlim = deg+2;
    2246                 :      57409 :   nb = 0;
    2247                 :      57409 :   i = 0;
    2248                 :            :   do
    2249                 :            :   {
    2250         [ +  + ]:     724275 :     for (j = 2; j < vlim; j++) gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
    2251                 :      57458 :     s = -signe(gel(v, vlim));
    2252                 :      57458 :     vlim--;
    2253                 :      57458 :     i++;
    2254                 :            :   }
    2255         [ +  + ]:      57458 :   while (!s);
    2256 [ +  + ][ +  + ]:      57409 :   if (vlim != deg + 1 && Premapped) setlg(v, vlim + 2);
    2257                 :            : 
    2258         [ +  + ]:     374759 :   for (; i <= deg - 1; i++)
    2259                 :            :   {
    2260                 :     353530 :     s2 = -signe(gel(v, 2));
    2261                 :     353530 :     flag = (s2 == s);
    2262         [ +  + ]:    4801906 :     for (j = 2; j < vlim; j++)
    2263                 :            :     {
    2264                 :    4448376 :       gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
    2265         [ +  + ]:    4448376 :       if (flag) flag = (s2 != signe(gel(v, j+1)));
    2266                 :            :     }
    2267         [ +  + ]:     353530 :     if (s == signe(gel(v, vlim)))
    2268                 :            :     {
    2269                 :      56029 :       nb++;
    2270         [ +  + ]:      56029 :       if (nb >= 2) { avma = av; return 2; }
    2271                 :      36264 :       s = -s;
    2272                 :            :     }
    2273         [ -  + ]:     333765 :     if (gc_needed(av, 3))
    2274                 :            :     {
    2275         [ #  # ]:          0 :       if (!Premapped) setlg(v, vlim);
    2276                 :          0 :       v = gerepileupto(av, v);
    2277         [ #  # ]:          0 :       if (DEBUGMEM > 1) pari_warn(warnmem, "X2XP1");
    2278                 :            :     }
    2279 [ +  + ][ +  + ]:     333765 :     if (flag && !Premapped) goto END;
    2280                 :     317350 :     vlim--;
    2281                 :            :   }
    2282         [ +  + ]:      21229 :   if (s == signe(gel(v, vlim))) nb++;
    2283                 :            : END:
    2284 [ +  + ][ +  + ]:      37644 :   if (Premapped && nb == 1)
    2285                 :       6807 :     *Premapped = gerepileupto(av, v);
    2286                 :            :   else
    2287                 :      30837 :     avma = av;
    2288                 :      57409 :   return nb;
    2289                 :            : }
    2290                 :            : 
    2291                 :            : static long
    2292                 :          0 : _intervalcmp(GEN x, GEN y)
    2293                 :            : {
    2294         [ #  # ]:          0 :   if (typ(x) == t_VEC) x = gel(x, 1);
    2295         [ #  # ]:          0 :   if (typ(y) == t_VEC) y = gel(y, 1);
    2296                 :          0 :   return gcmp(x, y);
    2297                 :            : }
    2298                 :            : 
    2299                 :            : static GEN
    2300                 :     626403 : _gen_nored(void *E, GEN x) { (void)E; return x; }
    2301                 :            : static GEN
    2302                 :    3422428 : _mp_add(void *E, GEN x, GEN y) { (void)E; return mpadd(x, y); }
    2303                 :            : static GEN
    2304                 :     726678 : _mp_mul(void *E, GEN x, GEN y) { (void)E; return mpmul(x, y); }
    2305                 :            : static GEN
    2306                 :     708995 : _mp_sqr(void *E, GEN x) { (void)E; return mpsqr(x); }
    2307                 :            : static GEN
    2308                 :     719750 : _gen_one(void *E) { (void)E; return gen_1; }
    2309                 :            : static GEN
    2310                 :       8688 : _gen_zero(void *E) { (void)E; return gen_0; }
    2311                 :            : 
    2312                 :            : static struct bb_algebra mp_algebra = { _gen_nored,_mp_add,_mp_mul,_mp_sqr,_gen_one,_gen_zero };
    2313                 :            : 
    2314                 :            : static GEN
    2315                 :    3885314 : _mp_cmul(void *E, GEN P, long a, GEN x) {(void)E; return mpmul(gel(P,a+2), x);}
    2316                 :            : 
    2317                 :            : /* Split the polynom P in two parts, each with unique sign.
    2318                 :            :  * Moreover compute the two parts of the derivative of P. */
    2319                 :            : static long
    2320                 :       6807 : split_polynoms(GEN P, long deg, long s0, GEN *Pp, GEN *Pm, GEN *Pprimep, GEN *Pprimem)
    2321                 :            : {
    2322                 :       6807 :   long i, degneg, v = evalvarn(varn(P));
    2323 [ +  + ][ +  - ]:      51440 :   for(i=1; i <= deg; i++) if (signe(gel(P, i+2)) != s0) break;
    2324                 :       6807 :   degneg = i;
    2325                 :       6807 :   *Pm = cgetg(degneg + 2, t_POL);
    2326                 :       6807 :   (*Pm)[1] = v;
    2327                 :       6807 :   *Pprimem = cgetg(degneg + 1, t_POL);
    2328                 :       6807 :   (*Pprimem)[1] = v;
    2329         [ +  + ]:      58247 :   for(i=0; i < degneg; i++)
    2330                 :            :   {
    2331                 :      51440 :     GEN elt = gel(P, i+2);
    2332                 :      51440 :     gel(*Pm, i+2) = elt;
    2333         [ +  + ]:      51440 :     if (i) gel(*Pprimem, i+1) = mului(i, elt);
    2334                 :            :   }
    2335                 :       6807 :   *Pp = cgetg(deg - degneg + 3, t_POL);
    2336                 :       6807 :   (*Pp)[1] = v;
    2337                 :       6807 :   *Pprimep = cgetg(deg - degneg + 3, t_POL);
    2338                 :       6807 :   (*Pprimep)[1] = v;
    2339         [ +  + ]:      58039 :   for(i=degneg; i <= deg; i++)
    2340                 :            :   {
    2341                 :      51232 :     GEN elt = gel(P, i+2);
    2342                 :      51232 :     gel(*Pp, i+2-degneg) = elt;
    2343                 :      51232 :     gel(*Pprimep, i+2-degneg) = mului(i, elt);
    2344                 :            :   }
    2345                 :       6807 :   *Pm = normalizepol_lg(*Pm, degneg+2);
    2346                 :       6807 :   *Pprimem = normalizepol_lg(*Pprimem, degneg+1);
    2347                 :       6807 :   *Pp = normalizepol_lg(*Pp, deg-degneg+3);
    2348                 :       6807 :   *Pprimep = normalizepol_lg(*Pprimep, deg-degneg+3);
    2349                 :       6807 :   return degpol(*Pm) + 1;
    2350                 :            : }
    2351                 :            : 
    2352                 :            : static GEN
    2353                 :     235787 : bkeval_single_power(long d, GEN V)
    2354                 :            : {
    2355                 :     235787 :   long mp = lg(V) - 2;
    2356         [ -  + ]:     235787 :   if (d > mp) return gmul(gpowgs(gel(V, mp+1), d/mp), gel(V, (d%mp)+1));
    2357                 :     235787 :   return gel(V, d+1);
    2358                 :            : }
    2359                 :            : 
    2360                 :            : static GEN
    2361                 :     235787 : splitpoleval(GEN Pp, GEN Pm, GEN pows, long deg, long degneg, long bitprec)
    2362                 :            : {
    2363                 :     235787 :   GEN vp = gen_bkeval_powers(Pp, deg-degneg, pows, NULL, &mp_algebra, _mp_cmul);
    2364                 :     235787 :   GEN vm = gen_bkeval_powers(Pm, degneg-1, pows, NULL, &mp_algebra, _mp_cmul);
    2365                 :     235787 :   GEN xa = bkeval_single_power(degneg, pows);
    2366                 :            :   GEN r;
    2367                 :     235787 :   vp = gmul(vp, xa);
    2368                 :     235787 :   r = gadd(vp, vm);
    2369 [ +  + ][ +  + ]:     235787 :   if (expo(vp) - (signe(r) ? expo(r) : 0) > prec2nbits(realprec(vp)) - bitprec)
    2370                 :      17396 :     return NULL;
    2371                 :     235787 :   return r;
    2372                 :            : }
    2373                 :            : 
    2374                 :            : /* Newton for polynom P. One solution between 0 and infinity, P' has also at
    2375                 :            :  * most one zero. P is almost certainly without zero coefficients. */
    2376                 :            : static GEN
    2377                 :       6807 : polsolve(GEN P, long bitprec)
    2378                 :            : {
    2379                 :       6807 :   pari_sp av = avma;
    2380                 :            :   GEN Pp, Pm, Pprimep, Pprimem, pows;
    2381                 :            :   GEN Pprime, Pprime2;
    2382                 :       6807 :   GEN ra, rb, rc, rcold = NULL, Pc, Ppc;
    2383                 :       6807 :   long deg = degpol(P), degneg, rt;
    2384                 :       6807 :   long s0, cprec = DEFAULTPREC, prec = nbits2prec(bitprec);
    2385                 :            :   long bitaddprec, addprec;
    2386                 :       6807 :   long expoold = LONG_MAX, iter;
    2387         [ -  + ]:       6807 :   if (deg == 1)
    2388                 :            :   {
    2389                 :          0 :     ra = negr(divrr(itor(gel(P, 2), bitprec), itor(gel(P, 3), bitprec)));
    2390                 :          0 :     return gerepileuptoleaf(av, ra);
    2391                 :            :   }
    2392                 :       6807 :   Pprime = ZX_deriv(P); Pprime2 = ZX_deriv(Pprime);
    2393                 :       6807 :   bitaddprec = 1 + 2*expu(deg); addprec = nbits2prec(bitaddprec);
    2394                 :       6807 :   s0 = signe(gel(P, 2));
    2395                 :       6807 :   degneg = split_polynoms(P, deg, s0, &Pp, &Pm, &Pprimep, &Pprimem);
    2396                 :       6807 :   rt = maxss(degneg, brent_kung_optpow(maxss(deg-degneg, degneg-1), 2, 1));
    2397                 :            :   {
    2398                 :            :     /* Optimized Cauchy bound */
    2399                 :       6807 :     GEN summin = gen_0, absmaj;
    2400                 :            :     long i;
    2401                 :            : 
    2402                 :       6807 :     absmaj = gel(P, 2);
    2403         [ +  + ]:      51440 :     for(i=1; i < degneg; i++)
    2404                 :            :     {
    2405                 :      44633 :       GEN elt = gel(P, i+2);
    2406         [ +  + ]:      44633 :       if (absi_cmp(absmaj, elt) < 0) absmaj = elt;
    2407                 :            :     }
    2408                 :       6807 :     summin = gel(P, i+2);
    2409         [ +  + ]:      51232 :     while (++i <= deg) summin = addii(summin, gel(P, i+2));
    2410                 :       6807 :     rb = subsr(1, rdivii(absmaj, summin, cprec));
    2411                 :            :     do
    2412                 :            :     {
    2413                 :       6807 :       pows = gen_powers(rb, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
    2414                 :       6807 :       Pc = splitpoleval(Pp, Pm, pows, deg, degneg, bitaddprec);
    2415         [ -  + ]:       6807 :       if (!Pc) { cprec++; rb = rtor(rb, cprec); continue; }
    2416         [ +  - ]:       6807 :       if (signe(Pc) != s0) break;
    2417                 :          0 :       shiftr_inplace(rb,1);
    2418                 :            :     }
    2419                 :          0 :     while (1);
    2420                 :            :   }
    2421                 :       6807 :   ra = NULL;
    2422                 :       6807 :   iter = 0;
    2423                 :            :   while (1)
    2424                 :            :   {
    2425                 :            :     GEN wdth;
    2426                 :      66272 :     iter++;
    2427         [ +  + ]:      66272 :     if (ra)
    2428                 :      43323 :       rc = shiftr(addrr(ra, rb), -1);
    2429                 :            :     else
    2430                 :      22949 :       rc = shiftr(rb, -1);
    2431                 :            :     do
    2432                 :            :     {
    2433                 :      66272 :       pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
    2434                 :      66272 :       Pc = splitpoleval(Pp, Pm, pows, deg, degneg, bitaddprec+2);
    2435         [ +  - ]:      66272 :       if (Pc) break;
    2436                 :          0 :       cprec++; rc = rtor(rc, cprec);
    2437                 :          0 :     } while (1);
    2438         [ +  + ]:      66272 :     if (signe(Pc) == s0)
    2439                 :      27784 :       ra = rc;
    2440                 :            :     else
    2441                 :      38488 :       rb = rc;
    2442         [ +  + ]:      66272 :     if (!ra) continue;
    2443                 :      50130 :     wdth = subrr(rb, ra);
    2444         [ +  + ]:      50130 :     if (!(iter % 8))
    2445                 :            :     {
    2446                 :       7451 :       GEN m1 = poleval(Pprime, ra), M2;
    2447         [ +  + ]:       7451 :       if (signe(m1) == s0) continue;
    2448                 :       7192 :       M2 = poleval(Pprime2, rb);
    2449         [ +  + ]:       7192 :       if (absr_cmp(gmul(M2, wdth), shiftr(m1, 1)) > 0) continue;
    2450                 :       6807 :       break;
    2451                 :            :     }
    2452         [ -  + ]:      42679 :     else if (gexpo(wdth) <= -bitprec)
    2453                 :          0 :       break;
    2454                 :      59465 :   }
    2455                 :       6807 :   rc = rb;
    2456                 :       6807 :   iter = 0;
    2457                 :            :   while (1)
    2458                 :            :   {
    2459                 :            :     long exponew;
    2460                 :            :     GEN dist;
    2461                 :      81354 :     iter++;
    2462                 :      81354 :     rcold = rc;
    2463                 :      81354 :     pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
    2464                 :      81354 :     Ppc = splitpoleval(Pprimep, Pprimem, pows, deg-1, degneg-1, bitaddprec+4);
    2465         [ +  - ]:      81354 :     if (Ppc)
    2466                 :      81354 :       Pc = splitpoleval(Pp, Pm, pows, deg, degneg, bitaddprec+4);
    2467 [ +  - ][ +  + ]:      81354 :     if (!Ppc || !Pc)
    2468                 :            :     {
    2469         [ +  + ]:      17396 :       if (cprec >= prec+addprec)
    2470                 :       2081 :         cprec++;
    2471                 :            :       else
    2472                 :      15315 :         cprec = minss(2*cprec, prec+addprec);
    2473                 :      17396 :       rc = rtor(rc, cprec); /* Backtrack one step */
    2474                 :      17396 :       continue;
    2475                 :            :     }
    2476                 :      63958 :     dist = divrr(Pc, Ppc);
    2477                 :      63958 :     rc = subrr(rc, dist);
    2478 [ +  - ][ -  + ]:      63958 :     if (cmprr(ra, rc) > 0 || cmprr(rb, rc) < 0)
    2479                 :            :     {
    2480         [ #  # ]:          0 :       if (cprec >= prec+addprec) break;
    2481                 :          0 :       cprec = minss(2*cprec, prec+addprec);
    2482                 :          0 :       rc = rtor(rcold, cprec); /* Backtrack one step */
    2483                 :          0 :       continue;
    2484                 :            :     }
    2485         [ +  + ]:      63958 :     if (expoold == LONG_MAX) { expoold = expo(dist); continue; }
    2486                 :      55090 :     exponew = expo(dist);
    2487         [ +  + ]:      55090 :     if (exponew < -bitprec - 1)
    2488                 :            :     {
    2489         [ +  + ]:       7757 :       if (cprec >= prec+addprec) break;
    2490                 :        950 :       cprec = minss(2*cprec, prec+addprec);
    2491                 :        950 :       rc = rtor(rc, cprec);
    2492                 :        950 :       continue;
    2493                 :            :     }
    2494         [ +  + ]:      47333 :     if (exponew > expoold - 2)
    2495                 :            :     {
    2496         [ -  + ]:       2061 :       if (cprec >= prec+addprec) break;
    2497                 :       2061 :       cprec = minss(2*cprec, prec+addprec);
    2498                 :       2061 :       rc = rtor(rc, cprec);
    2499                 :       2061 :       expoold = LONG_MAX;
    2500                 :       2061 :       continue;
    2501                 :            :     }
    2502                 :      45272 :     expoold = exponew;
    2503                 :      74547 :   }
    2504                 :       6807 :   return gerepileuptoleaf(av, rtor(rc, prec));
    2505                 :            : }
    2506                 :            : 
    2507                 :            : static GEN
    2508                 :       6851 : usp(GEN Q0, long deg, long *nb_donep, long flag, long bitprec)
    2509                 :            : {
    2510                 :            :   pari_sp av;
    2511                 :            :   GEN Q, sol;
    2512                 :       6851 :   long nb_todo, nbr = 0, ind, deg0, indf, i, k, nb, j;
    2513                 :       6851 :   long listsize = 64, nb_done = 0;
    2514                 :            :   GEN c, Lc, Lk;
    2515                 :            : 
    2516                 :       6851 :   av = avma;
    2517                 :            : 
    2518                 :       6851 :   sol = const_col(deg, gen_0);
    2519                 :       6851 :   deg0 = deg;
    2520                 :       6851 :   Lc = const_vec(listsize, gen_0);
    2521                 :       6851 :   Lk = cgetg(listsize+1, t_VECSMALL);
    2522                 :       6851 :   c = gen_0;
    2523                 :       6851 :   k = Lk[1] = 0;
    2524                 :       6851 :   ind = 1; indf = 2;
    2525                 :       6851 :   Q = gcopy(Q0);
    2526                 :            : 
    2527                 :       6851 :   nb_todo = 1;
    2528         [ +  + ]:      64260 :   while (nb_todo)
    2529                 :            :   {
    2530                 :      57409 :     GEN nc = gel(Lc, ind), Qremapped;
    2531         [ +  + ]:      57409 :     if (Lk[ind] == k + 1)
    2532                 :            :     {
    2533                 :      22178 :       deg0 = deg;
    2534                 :      22178 :       setlg(Q0, deg + 3);
    2535                 :      22178 :       Q0 = ZX_rescale2n(Q0, 1);
    2536                 :      22178 :       Q = Q_primpart(Q0);
    2537                 :      22178 :       c = gen_0;
    2538                 :            :     }
    2539                 :            : 
    2540         [ +  + ]:      57409 :     if (!equalii(nc, c)) Q = ZX_translate(Q, subii(nc, c));
    2541                 :            : 
    2542                 :      57409 :     k = Lk[ind];
    2543                 :      57409 :     c = nc;
    2544                 :      57409 :     ind++;
    2545                 :      57409 :     nb_todo--;
    2546                 :            : 
    2547         [ +  + ]:      57409 :     if (equalii(gel(Q, 2), gen_0))
    2548                 :            :     {
    2549                 :        126 :       GEN newsol = gmul2n(c, -k);
    2550         [ +  + ]:        168 :       for (j = 1; j <= nbr; j++)
    2551         [ +  + ]:         77 :         if (gequal(gel(sol, j), newsol)) break;
    2552         [ +  + ]:        126 :       if (j > nbr) gel(sol, ++nbr) = newsol;
    2553                 :            : 
    2554                 :        126 :       deg0--;
    2555         [ +  + ]:        504 :       for (j = 2; j <= deg0 + 2; j++) gel(Q, j) = gel(Q, j+1);
    2556                 :        126 :       setlg(Q, j);
    2557                 :            :     }
    2558                 :            : 
    2559         [ +  + ]:      57409 :     nb = X2XP1(Q, deg0, flag == 1 ? &Qremapped : NULL);
    2560                 :      57409 :     nb_done++;
    2561                 :            : 
    2562      [ +  +  + ]:      57409 :     switch (nb)
    2563                 :            :     {
    2564                 :            :       case 0:
    2565                 :      19359 :         break;
    2566                 :            :       case 1:
    2567      [ -  +  + ]:      12771 :         switch(flag)
    2568                 :            :         {
    2569                 :            :         case 0:
    2570                 :            :           {
    2571                 :            :             GEN low, hi;
    2572                 :          0 :             low = gmul2n(c, -k);
    2573                 :          0 :             hi  = gmul2n(addiu(c,1), -k);
    2574                 :          0 :             gel(sol, ++nbr) = mkvec2(low, hi);
    2575                 :            :           }
    2576                 :          0 :           break;
    2577                 :            :         case 1:
    2578                 :            :           { /* Caveat emptor: Qremapped is the reciprocal polynomial */
    2579                 :       6807 :             GEN sr = polsolve(Qremapped, bitprec+1);
    2580                 :       6807 :             sr = divrr(sr, addsr(1, sr));
    2581                 :       6807 :             sr = gmul2n(addir(c, sr), -k);
    2582                 :       6807 :             gel(sol, ++nbr) = rtor(sr, nbits2prec(bitprec));
    2583                 :            :           }
    2584                 :       6807 :           break;
    2585                 :            :         default:
    2586                 :       5964 :           gel(sol, ++nbr) = gen_0;
    2587                 :            :         }
    2588                 :      12771 :       break;
    2589                 :            : 
    2590                 :            :       default:
    2591         [ +  + ]:      25279 :         if (indf + 2 > listsize)
    2592                 :            :         {
    2593         [ +  - ]:        168 :           if (ind>1)
    2594                 :            :           {
    2595         [ +  + ]:        833 :             for (i = ind; i < indf; i++)
    2596                 :            :             {
    2597                 :        665 :               gel(Lc, i-ind+1) = gel(Lc, i);
    2598                 :        665 :               Lk[i-ind+1] = Lk[i];
    2599                 :            :             }
    2600                 :        168 :           indf -= ind-1; ind = 1;
    2601                 :            :           }
    2602         [ -  + ]:        168 :           if (indf + 2 > listsize)
    2603                 :            :           {
    2604                 :          0 :             listsize *= 2;
    2605                 :          0 :             Lc = vec_lengthen(Lc, listsize);
    2606                 :          0 :             Lk = vecsmall_lengthen(Lk, listsize);
    2607                 :            :           }
    2608         [ +  + ]:      10255 :           for (i = indf; i <= listsize; i++) gel(Lc, i) = gen_0;
    2609                 :            :         }
    2610                 :      25279 :         nc = shifti(c, 1);
    2611                 :      25279 :         gel(Lc, indf) = nc;
    2612                 :      25279 :         gel(Lc, indf + 1) = addis(nc, 1);
    2613                 :      25279 :         Lk[indf] = Lk[indf + 1] = k + 1;
    2614                 :      25279 :         indf += 2;
    2615                 :      25279 :         nb_todo += 2;
    2616                 :            :     }
    2617                 :            : 
    2618         [ -  + ]:      57409 :     if (gc_needed(av, 2))
    2619                 :            :     {
    2620                 :          0 :       gerepileall(av, 6, &Q0, &Q, &c, &Lc, &Lk, &sol);
    2621         [ #  # ]:      57409 :       if (DEBUGMEM > 1) pari_warn(warnmem, "ZX_uspensky", avma);
    2622                 :            :     }
    2623                 :            :   }
    2624                 :            : 
    2625                 :       6851 :   setlg(sol, nbr+1);
    2626                 :       6851 :   *nb_donep += nb_done;
    2627                 :       6851 :   return gerepilecopy(av, sol);
    2628                 :            : }
    2629                 :            : 
    2630                 :            : static GEN
    2631                 :       1554 : ZX_uspensky_cst_pol(long nbz, long flag, long bitprec)
    2632                 :            : {
    2633      [ -  -  + ]:       1554 :   switch(flag)
    2634                 :            :   {
    2635                 :          0 :     case 0:  return zerocol(nbz);
    2636         [ #  # ]:          0 :     case 1:  retconst_col(nbz, real_0_bit(bitprec));
    2637                 :       1554 :     default: return utoi(nbz);
    2638                 :            :   }
    2639                 :            : }
    2640                 :            : 
    2641                 :            : GEN
    2642                 :       8513 : ZX_uspensky(GEN P, GEN ab, long flag, long bitprec)
    2643                 :            : {
    2644                 :       8513 :   pari_sp av = avma;
    2645                 :       8513 :   GEN a, b, sol = NULL, Pcur;
    2646                 :            :   double fb;
    2647                 :       8513 :   long nbz, deg, nb_done = 0;
    2648                 :            : 
    2649                 :       8513 :   deg = degpol(P);
    2650 [ +  + ][ +  - ]:       8513 :   if (deg == 0) return flag <= 1 ? cgetg(1, t_COL) : gen_0;
    2651         [ +  + ]:       7946 :   if (ab)
    2652                 :            :   {
    2653         [ +  + ]:       3299 :     if (typ(ab) == t_VEC)
    2654                 :            :     {
    2655         [ -  + ]:        119 :       if (lg(ab) != 3) pari_err_DIM("ZX_uspensky");
    2656                 :        119 :       a = gel(ab, 1);
    2657                 :        119 :       b = gel(ab, 2);
    2658                 :            :     }
    2659                 :            :     else
    2660                 :            :     {
    2661                 :       3180 :       a = ab;
    2662                 :       3180 :       b = mkoo();
    2663                 :            :     }
    2664                 :            :   }
    2665                 :            :   else
    2666                 :            :   {
    2667                 :       4647 :     a = mkmoo();
    2668                 :       4647 :     b = mkoo();
    2669                 :            :   }
    2670      [ +  +  + ]:       7946 :   switch (gcmp(a, b))
    2671                 :            :   {
    2672         [ -  + ]:          7 :     case 1: avma = av; return flag <= 1 ? cgetg(1, t_COL) : gen_0;
    2673                 :            :     case 0:
    2674 [ +  - ][ +  + ]:         21 :       if (typ(a) != t_INFINITY && gequal0(poleval(P, a)))
    2675         [ +  + ]:         14 :       { avma = av; return flag <= 1 ? mkcolcopy(a): gen_1; }
    2676                 :            :       else
    2677         [ +  - ]:          7 :       { avma = av; return flag <= 1 ? cgetg(1, t_COL) : gen_0; }
    2678                 :            :   }
    2679                 :       7918 :   nbz = ZX_valrem(P, &Pcur);
    2680                 :       7918 :   deg -= nbz;
    2681         [ +  + ]:       7918 :   if (!nbz) Pcur = P;
    2682 [ +  + ][ +  - ]:       7918 :   if (nbz && (gsigne(a) > 0 || gsigne(b) < 0)) nbz = 0;
                 [ -  + ]
    2683         [ +  + ]:       7918 :   if (deg == 0) { avma = av; return ZX_uspensky_cst_pol(nbz, flag, bitprec); }
    2684         [ +  + ]:       7351 :   if (deg == 1)
    2685                 :            :   {
    2686                 :       2739 :     sol = gdiv(gneg(gel(Pcur, 2)), pollead(Pcur, -1));
    2687 [ +  + ][ -  + ]:       2739 :     if (gcmp(a, sol) > 0 || gcmp(sol, b) > 0)
    2688                 :            :     {
    2689                 :        987 :       avma = av;
    2690                 :        987 :       return ZX_uspensky_cst_pol(nbz, flag, bitprec);
    2691                 :            :     }
    2692         [ +  + ]:       1752 :     if (flag >= 2) { avma = av; return utoi(nbz+1); }
    2693                 :        737 :     sol = concat(zerocol(nbz), mkcol(sol));
    2694         [ +  - ]:        737 :     if (flag == 1) sol = RgC_gtofp(sol, nbits2prec(bitprec));
    2695                 :        737 :     return gerepilecopy(av, sol);
    2696                 :            :   }
    2697      [ -  +  + ]:       4612 :   switch(flag)
    2698                 :            :   {
    2699                 :            :     case 0:
    2700                 :          0 :       sol = zerocol(nbz);
    2701                 :          0 :       break;
    2702                 :            :     case 1:
    2703                 :       1364 :       sol = const_col(nbz, real_0_bit(bitprec));
    2704                 :       1364 :       break;
    2705                 :            :     /* case 2: nothing */
    2706                 :            :   }
    2707                 :            : 
    2708 [ +  + ][ +  + ]:       4612 :   if (typ(a) == t_INFINITY && typ(b) != t_INFINITY && gsigne(b))
                 [ +  - ]
    2709                 :            :   {
    2710                 :         28 :     fb = fujiwara_bound_real(Pcur, -1);
    2711         [ +  + ]:         28 :     if (fb > -pariINFINITY)
    2712                 :          7 :       a = negi(int2n((long)ceil(fb)));
    2713                 :            :     else
    2714                 :         21 :       a = gen_0;
    2715                 :            :   }
    2716 [ +  + ][ +  + ]:       4612 :   if (typ(b) == t_INFINITY && typ(a) != t_INFINITY && gsigne(a))
                 [ +  + ]
    2717                 :            :   {
    2718                 :         14 :     fb = fujiwara_bound_real(Pcur, 1);
    2719         [ +  - ]:         14 :     if (fb > -pariINFINITY)
    2720                 :         14 :       b = int2n((long)ceil(fb));
    2721                 :            :     else
    2722                 :          0 :       b = gen_0;
    2723                 :            :   }
    2724                 :            : 
    2725 [ +  + ][ +  + ]:       4612 :   if (typ(a) != t_INFINITY && typ(b) != t_INFINITY)
    2726                 :            :   {
    2727                 :            :     pari_sp av1;
    2728                 :         91 :     GEN den = lcmii(denom(a), denom(b)), diff, unscaledres, co, Pdiv;
    2729                 :            :     GEN ascaled;
    2730                 :            :     long i;
    2731         [ +  + ]:         91 :     if (!is_pm1(den))
    2732                 :            :     {
    2733                 :         28 :       Pcur = ZX_rescale(Pcur, den);
    2734                 :         28 :       ascaled = gmul(a, den);
    2735                 :            :     }
    2736                 :            :     else
    2737                 :            :     {
    2738                 :         63 :       den = NULL;
    2739                 :         63 :       ascaled = a;
    2740                 :            :     }
    2741         [ +  + ]:         91 :     diff = subii(den ? gmul(b,den) : b, ascaled);
    2742                 :         91 :     Pcur = ZX_unscale(ZX_translate(Pcur, ascaled), diff);
    2743                 :         91 :     av1 = avma;
    2744                 :         91 :     Pdiv = cgetg(deg+2, t_POL);
    2745                 :         91 :     Pdiv[1] = Pcur[1];
    2746                 :         91 :     co = gel(Pcur, deg+2);
    2747         [ +  + ]:        364 :     for (i = deg; --i >= 0; )
    2748                 :            :     {
    2749                 :        273 :       gel(Pdiv, i+2) = co;
    2750                 :        273 :       co = addii(co, gel(Pcur, i+2));
    2751                 :            :     }
    2752         [ +  + ]:         91 :     if (!signe(co))
    2753                 :            :     {
    2754                 :         49 :       Pcur = Pdiv;
    2755                 :         49 :       deg--;
    2756         [ +  + ]:         49 :       if (flag <= 1)
    2757                 :         14 :         sol = concat(sol, b);
    2758                 :            :       else
    2759                 :         35 :         nbz++;
    2760                 :            :     }
    2761                 :            :     else
    2762                 :         42 :       avma = av1;
    2763                 :         91 :     unscaledres = usp(Pcur, deg, &nb_done, flag, bitprec);
    2764         [ +  + ]:         91 :     if (flag <= 1)
    2765                 :            :     {
    2766         [ +  + ]:         63 :       for (i = 1; i < lg(unscaledres); i++)
    2767                 :            :       {
    2768                 :         35 :         GEN z = gmul(diff, gel(unscaledres, i));
    2769         [ -  + ]:         35 :         if (typ(z) == t_VEC)
    2770                 :            :         {
    2771                 :          0 :           gel(z, 1) = gadd(ascaled, gel(z, 1));
    2772                 :          0 :           gel(z, 2) = gadd(ascaled, gel(z, 2));
    2773                 :            :         }
    2774                 :            :         else
    2775                 :         35 :           z = gadd(ascaled, z);
    2776         [ +  + ]:         35 :         if (den) z = gdiv(z, den);
    2777                 :         35 :         gel(unscaledres, i) = z;
    2778                 :            :       }
    2779                 :         28 :       sol = concat(sol, unscaledres);
    2780                 :            :     }
    2781                 :            :     else
    2782                 :         63 :       nbz += lg(unscaledres) - 1;
    2783                 :            :   }
    2784 [ +  + ][ +  + ]:       4612 :   if (typ(b) == t_INFINITY && (fb=fujiwara_bound_real(Pcur, 1)) > -pariINFINITY)
    2785                 :            :   {
    2786                 :            :     GEN Pcurp, unscaledres;
    2787                 :       3660 :     long bp = (long)ceil(fb);
    2788         [ +  + ]:       3660 :     if (bp < 0) bp = 0;
    2789                 :       3660 :     Pcurp = ZX_unscale2n(Pcur, bp);
    2790                 :       3660 :     unscaledres = usp(Pcurp, deg, &nb_done, flag, bitprec);
    2791         [ +  + ]:       3660 :     if (flag <= 1)
    2792                 :       1329 :       sol = concat(sol, gmul2n(unscaledres, bp));
    2793                 :            :     else
    2794                 :       2331 :       nbz += lg(unscaledres) - 1;
    2795                 :            :   }
    2796 [ +  + ][ +  + ]:       4612 :   if (typ(a) == t_INFINITY && (fb=fujiwara_bound_real(Pcur,-1)) > -pariINFINITY)
    2797                 :            :   {
    2798                 :            :     GEN Pcurm, unscaledres;
    2799                 :       3100 :     long i, bm = (long)ceil(fb);
    2800         [ +  + ]:       3100 :     if (bm < 0) bm = 0;
    2801                 :       3100 :     Pcurm = ZX_unscale(Pcur, gen_m1);
    2802                 :       3100 :     Pcurm = ZX_unscale2n(Pcurm, bm);
    2803                 :       3100 :     unscaledres = usp(Pcurm,deg,&nb_done,flag,bitprec);
    2804         [ +  + ]:       3100 :     if (flag <= 1)
    2805                 :            :     {
    2806         [ +  + ]:       3168 :       for (i = 1; i < lg(unscaledres); i++)
    2807                 :            :       {
    2808                 :       2133 :         GEN z = gneg(gmul2n(gel(unscaledres, i), bm));
    2809         [ -  + ]:       2133 :         if (typ(z) == t_VEC) swap(gel(z, 1), gel(z, 2));
    2810                 :       2133 :         gel(unscaledres, i) = z;
    2811                 :            :       }
    2812                 :       1035 :       sol = concat(unscaledres, sol);
    2813                 :            :     }
    2814                 :            :     else
    2815                 :       2065 :       nbz += lg(unscaledres) - 1;
    2816                 :            :   }
    2817                 :            : 
    2818         [ -  + ]:       4612 :   if (DEBUGLEVEL > 4)
    2819                 :          0 :     err_printf("ZX_uspensky: Number of visited nodes: %d\n", nb_done);
    2820                 :            : 
    2821         [ +  + ]:       4612 :   if (flag >= 2) return utoi(nbz);
    2822         [ +  - ]:       1364 :   if (flag)
    2823                 :       1364 :     sol = sort(sol);
    2824                 :            :   else
    2825                 :          0 :     sol = gen_sort(sol, (void *)_intervalcmp, cmp_nodata);
    2826                 :       8513 :   return gerepileupto(av, sol);
    2827                 :            : }
    2828                 :            : 
    2829                 :            : /* x a scalar */
    2830                 :            : static GEN
    2831                 :         35 : rootsdeg0(GEN x)
    2832                 :            : {
    2833         [ +  + ]:         35 :   if (!is_rational_t(typ(x))) pari_err_TYPE("realroots",x);
    2834         [ +  + ]:         28 :   if (gequal0(x)) pari_err_ROOTS0("realroots");
    2835                 :         14 :   return cgetg(1,t_COL); /* constant polynomial */
    2836                 :            : }
    2837                 :            : static void
    2838                 :       1190 : checkbound(GEN a)
    2839                 :            : {
    2840         [ +  - ]:       1190 :   switch(typ(a))
    2841                 :            :   {
    2842                 :       1190 :     case t_INT: case t_FRAC: case t_INFINITY: break;
    2843                 :          0 :     default: pari_err_TYPE("polrealroots", a);
    2844                 :            :   }
    2845                 :       1190 : }
    2846                 :            : static GEN
    2847                 :       3221 : check_ab(GEN ab)
    2848                 :            : {
    2849                 :            :   GEN a, b;
    2850         [ +  + ]:       3221 :   if (!ab) return NULL;
    2851 [ +  - ][ -  + ]:        595 :   if (typ(ab) != t_VEC || lg(ab) != 3) pari_err_TYPE("polrootsreal",ab);
    2852                 :        595 :   a = gel(ab,1); checkbound(a);
    2853                 :        595 :   b = gel(ab,2); checkbound(b);
    2854 [ +  + ][ +  - ]:        595 :   if (typ(a) == t_INFINITY && inf_get_sign(a) < 0 &&
                 [ +  + ]
    2855         [ +  - ]:        476 :       typ(b) == t_INFINITY && inf_get_sign(b) > 0) ab = NULL;
    2856                 :       3221 :   return ab;
    2857                 :            : }
    2858                 :            : GEN
    2859                 :       2675 : realroots(GEN P, GEN ab, long prec)
    2860                 :            : {
    2861                 :       2675 :   pari_sp av = avma;
    2862                 :       2675 :   long nrr = 0;
    2863                 :       2675 :   GEN sol = NULL, fa, ex;
    2864                 :            :   long i, j, k;
    2865                 :            : 
    2866                 :       2675 :   ab = check_ab(ab);
    2867         [ +  + ]:       2675 :   if (typ(P) != t_POL) return rootsdeg0(P);
    2868      [ +  +  + ]:       2654 :   switch(degpol(P))
    2869                 :            :   {
    2870                 :          7 :     case -1: return rootsdeg0(gen_0);
    2871                 :          7 :     case 0: return rootsdeg0(gel(P,2));
    2872                 :            :   }
    2873                 :       2640 :   P = Q_primpart(P);
    2874         [ -  + ]:       2640 :   if (!RgX_is_ZX(P)) pari_err_TYPE("realroots",P);
    2875                 :       2640 :   fa = ZX_squff(P, &ex);
    2876         [ +  + ]:       5287 :   for (i = 1; i < lg(fa); i++)
    2877                 :            :   {
    2878                 :       2647 :     GEN Pi = gel(fa, i), soli, soli2 = NULL;
    2879                 :       2647 :     long n, nrri = 0, h, nbz;
    2880         [ +  + ]:       2647 :     if (ab)
    2881                 :         42 :       h = 1;
    2882                 :            :     else
    2883                 :       2605 :       Pi = RgX_deflate_max(Pi, &h);
    2884         [ +  + ]:       2647 :     if (!signe(gel(Pi, 2)))
    2885                 :            :     {
    2886                 :        567 :       Pi = RgX_shift_shallow(Pi, -1);
    2887                 :        567 :       nbz = 1;
    2888                 :            :     }
    2889                 :            :     else
    2890                 :       2080 :       nbz = 0;
    2891         [ +  + ]:       2647 :     soli = ZX_uspensky(Pi, h%2 ? ab: gen_0, 1, prec2nbits(prec));
    2892                 :       2647 :     n = lg(soli);
    2893         [ +  + ]:       2647 :     if (!(h % 2)) soli2 = cgetg(n, t_COL);
    2894         [ +  + ]:       9316 :     for (j = 1; j < n; j++)
    2895                 :            :     {
    2896                 :       6669 :       GEN elt = gel(soli, j);
    2897         [ +  + ]:       6669 :       if (typ(elt) != t_REAL)
    2898                 :            :       {
    2899                 :         63 :         nrri++;
    2900                 :         63 :         elt = gtofp(elt, prec);
    2901                 :         63 :         gel(soli, j) = elt;
    2902                 :            :       }
    2903         [ +  + ]:       6669 :       if (h > 1)
    2904                 :            :       { /* note: elt != 0 because we are square free */
    2905                 :            :         GEN r;
    2906         [ +  + ]:       1262 :         if (h == 2)
    2907                 :       1248 :           r = sqrtr(elt);
    2908                 :            :         else
    2909                 :            :         {
    2910         [ +  + ]:         14 :           if (signe(elt) < 0)
    2911                 :          7 :             r = negr(sqrtnr(negr(elt), h));
    2912                 :            :           else
    2913                 :          7 :             r = sqrtnr(elt, h);
    2914                 :            :         }
    2915                 :       1262 :         gel(soli, j) = r;
    2916         [ +  + ]:       1262 :         if (!(h % 2)) gel(soli2, j) = negr(r);
    2917                 :            :       }
    2918                 :            :     }
    2919         [ +  + ]:       2647 :     if (!(h % 2)) soli = shallowconcat(soli, soli2);
    2920         [ +  + ]:       2647 :     if (nbz) soli = shallowconcat(soli, real_0(prec));
    2921         [ +  + ]:       5294 :     for (k = 1; k <= ex[i]; k++)
    2922         [ +  + ]:       2647 :       sol = sol ? shallowconcat(sol, soli) : soli;
    2923                 :       2647 :     nrr += ex[i]*nrri;
    2924                 :            :   }
    2925                 :            : 
    2926         [ -  + ]:       2640 :   if (DEBUGLEVEL > 4)
    2927                 :            :   {
    2928                 :          0 :     err_printf("Number of real roots: %d\n", lg(sol)-1);
    2929                 :          0 :     err_printf(" -- of which 2-integral: %ld\n", nrr);
    2930                 :            :   }
    2931                 :            : 
    2932                 :       2654 :   return gerepileupto(av, sort(sol));
    2933                 :            : }
    2934                 :            : 
    2935                 :            : /* P non-constant, squarefree ZX */
    2936                 :            : long
    2937                 :       5754 : ZX_sturm(GEN P)
    2938                 :            : {
    2939                 :       5754 :   pari_sp av = avma;
    2940                 :            :   long h, r;
    2941                 :       5754 :   P = RgX_deflate_max(P, &h);
    2942         [ +  + ]:       5754 :   if (odd(h))
    2943                 :       3367 :     r = itos(ZX_uspensky(P, NULL, 2, 0));
    2944                 :            :   else
    2945                 :       2387 :     r = 2*itos(ZX_uspensky(P, gen_0, 2, 0));
    2946                 :       5754 :   avma = av; return r;
    2947                 :            : }
    2948                 :            : /* P non-constant, squarefree ZX */
    2949                 :            : long
    2950                 :        546 : ZX_sturmpart(GEN P, GEN ab)
    2951                 :            : {
    2952                 :        546 :   pari_sp av = avma;
    2953                 :            :   long r;
    2954         [ +  + ]:        546 :   if (!check_ab(ab)) return ZX_sturm(P);
    2955                 :         77 :   r = itos(ZX_uspensky(P, ab, 2, 0));
    2956                 :        546 :   avma = av; return r;
    2957                 :            : }

Generated by: LCOV version 1.9