Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - language - intnum.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 20068-a012b06) Lines: 1210 1242 97.4 %
Date: 2017-01-16 05:49:55 Functions: 98 99 99.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : #include "anal.h"
      17             : 
      18             : static const long EXTRAPREC =
      19             : #ifdef LONG_IS_64BIT
      20             :   1;
      21             : #else
      22             :   2;
      23             : #endif
      24             : 
      25             : /********************************************************************/
      26             : /**                NUMERICAL INTEGRATION (Romberg)                 **/
      27             : /********************************************************************/
      28             : typedef struct {
      29             :   void *E;
      30             :   GEN (*f)(void *E, GEN);
      31             : } invfun;
      32             : 
      33             : /* 1/x^2 f(1/x) */
      34             : static GEN
      35       12474 : _invf(void *E, GEN x)
      36             : {
      37       12474 :   invfun *S = (invfun*)E;
      38       12474 :   GEN y = ginv(x);
      39       12474 :   return gmul(S->f(S->E, y), gsqr(y));
      40             : }
      41             : 
      42             : /* h and s are arrays of the same length L > D. The h[i] are (decreasing)
      43             :  * step sizes, s[i] is the computed Riemann sum for step size h[i].
      44             :  * Interpolate the last D+1 values so that s ~ polynomial in h of degree D.
      45             :  * Guess that limit_{h->0} = s(0) */
      46             : static GEN
      47         105 : interp(GEN h, GEN s, long L, long bit, long D)
      48             : {
      49         105 :   pari_sp av = avma;
      50             :   long e1,e2;
      51         105 :   GEN dss, ss = polint_i(h + L-D,s + L-D, gen_0, D+1, &dss);
      52             : 
      53         105 :   e1 = gexpo(ss);
      54         105 :   e2 = gexpo(dss);
      55         105 :   if (DEBUGLEVEL>2)
      56             :   {
      57           0 :     err_printf("romb: iteration %ld, guess: %Ps\n", L,ss);
      58           0 :     err_printf("romb: relative error < 2^-%ld [target %ld bits]\n",e1-e2,bit);
      59             :   }
      60         105 :   if (e1-e2 <= bit && (L <= 10 || e1 >= -bit)) { avma = av; return NULL; }
      61          70 :   if (typ(ss) == t_COMPLEX && gequal0(gel(ss,2))) ss = gel(ss,1);
      62          70 :   return ss;
      63             : }
      64             : 
      65             : static GEN
      66           7 : qrom3(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, long bit)
      67             : {
      68           7 :   const long JMAX = 25, KLOC = 4;
      69             :   GEN ss,s,h,p1,p2,qlint,del,x,sum;
      70           7 :   long j, j1, it, sig, prec = nbits2prec(bit);
      71             : 
      72           7 :   a = gtofp(a,prec);
      73           7 :   b = gtofp(b,prec);
      74           7 :   qlint = subrr(b,a); sig = signe(qlint);
      75           7 :   if (!sig) return gen_0;
      76           7 :   if (sig < 0) { setabssign(qlint); swap(a,b); }
      77             : 
      78           7 :   s = new_chunk(JMAX+KLOC-1);
      79           7 :   h = new_chunk(JMAX+KLOC-1);
      80           7 :   gel(h,0) = real_1(prec);
      81             : 
      82           7 :   p1 = eval(E, a); if (p1 == a) p1 = rcopy(p1);
      83           7 :   p2 = eval(E, b);
      84           7 :   gel(s,0) = gmul2n(gmul(qlint,gadd(p1,p2)),-1);
      85          28 :   for (it=1,j=1; j<JMAX; j++, it<<=1) /* it = 2^(j-1) */
      86             :   {
      87             :     pari_sp av, av2;
      88          28 :     gel(h,j) = real2n(-2*j, prec); /* 2^(-2j) */
      89          28 :     av = avma; del = divru(qlint,it);
      90          28 :     x = addrr(a, shiftr(del,-1));
      91          28 :     av2 = avma;
      92         133 :     for (sum = gen_0, j1 = 1; j1 <= it; j1++, x = addrr(x,del))
      93             :     {
      94         105 :       sum = gadd(sum, eval(E, x));
      95         105 :       if ((j1 & 0x1ff) == 0) gerepileall(av2, 2, &sum,&x);
      96             :     }
      97          28 :     sum = gmul(sum,del);
      98          28 :     gel(s,j) = gerepileupto(av, gmul2n(gadd(gel(s,j-1), sum), -1));
      99          28 :     if (j >= KLOC && (ss = interp(h, s, j, bit-j-6, KLOC)))
     100           7 :       return gmulsg(sig,ss);
     101             :   }
     102           0 :   pari_err_IMPL("intnumromb recovery [too many iterations]");
     103             :   return NULL; /* LCOV_EXCL_LINE */
     104             : }
     105             : 
     106             : static GEN
     107          63 : qrom2(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, long bit)
     108             : {
     109          63 :   const long JMAX = 16, KLOC = 4;
     110             :   GEN ss,s,h,p1,qlint,del,ddel,x,sum;
     111          63 :   long j, j1, it, sig, prec = nbits2prec(bit);
     112             : 
     113          63 :   a = gtofp(a, prec);
     114          63 :   b = gtofp(b, prec);
     115          63 :   qlint = subrr(b,a); sig = signe(qlint);
     116          63 :   if (!sig)  return gen_0;
     117          63 :   if (sig < 0) { setabssign(qlint); swap(a,b); }
     118             : 
     119          63 :   s = new_chunk(JMAX+KLOC-1);
     120          63 :   h = new_chunk(JMAX+KLOC-1);
     121          63 :   gel(h,0) = real_1(prec);
     122             : 
     123          63 :   p1 = shiftr(addrr(a,b),-1);
     124          63 :   gel(s,0) = gmul(qlint, eval(E, p1));
     125         287 :   for (it=1, j=1; j<JMAX; j++, it*=3) /* it = 3^(j-1) */
     126             :   {
     127             :     pari_sp av, av2;
     128         287 :     gel(h,j) = divru(gel(h,j-1), 9); /* 3^(-2j) */
     129         287 :     av = avma; del = divru(qlint,3*it); ddel = shiftr(del,1);
     130         287 :     x = addrr(a, shiftr(del,-1));
     131         287 :     av2 = avma;
     132        7910 :     for (sum = gen_0, j1 = 1; j1 <= it; j1++)
     133             :     {
     134        7623 :       sum = gadd(sum, eval(E, x)); x = addrr(x,ddel);
     135        7623 :       sum = gadd(sum, eval(E, x)); x = addrr(x,del);
     136        7623 :       if ((j1 & 0x1ff) == 0) gerepileall(av2, 2, &sum,&x);
     137             :     }
     138         287 :     sum = gmul(sum,del); p1 = gdivgs(gel(s,j-1),3);
     139         287 :     gel(s,j) = gerepileupto(av, gadd(p1,sum));
     140         287 :     if (j >= KLOC && (ss = interp(h, s, j, bit-(3*j/2)+3, KLOC)))
     141          63 :       return gmulsg(sig, ss);
     142             :   }
     143           0 :   pari_err_IMPL("intnumromb recovery [too many iterations]");
     144             :   return NULL; /* LCOV_EXCL_LINE */
     145             : }
     146             : 
     147             : /* integrate after change of variables x --> 1/x */
     148             : static GEN
     149          28 : qromi(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, long bit)
     150             : {
     151          28 :   GEN A = ginv(b), B = ginv(a);
     152             :   invfun S;
     153          28 :   S.f = eval;
     154          28 :   S.E = E; return qrom2(&S, &_invf, A, B, bit);
     155             : }
     156             : 
     157             : /* a < b, assume b "small" (< 100 say) */
     158             : static GEN
     159          28 : rom_bsmall(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, long bit)
     160             : {
     161          28 :   if (gcmpgs(a,-100) >= 0) return qrom2(E,eval,a,b,bit);
     162           7 :   if (gcmpgs(b, -1) < 0)   return qromi(E,eval,a,b,bit); /* a<-100, b<-1 */
     163             :   /* a<-100, b>=-1, split at -1 */
     164           7 :   return gadd(qromi(E,eval,a,gen_m1,bit),
     165             :               qrom2(E,eval,gen_m1,b,bit));
     166             : }
     167             : 
     168             : static GEN
     169          35 : rombint(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, long bit)
     170             : {
     171          35 :   long l = gcmp(b,a);
     172             :   GEN z;
     173             : 
     174          35 :   if (!l) return gen_0;
     175          35 :   if (l < 0) swap(a,b);
     176          35 :   if (gcmpgs(b,100) >= 0)
     177             :   {
     178          14 :     if (gcmpgs(a,1) >= 0)
     179           7 :       z = qromi(E,eval,a,b,bit);
     180             :     else /* split at 1 */
     181           7 :       z = gadd(rom_bsmall(E,eval,a,gen_1,bit), qromi(E,eval,gen_1,b,bit));
     182             :   }
     183             :   else
     184          21 :     z = rom_bsmall(E,eval,a,b,bit);
     185          35 :   if (l < 0) z = gneg(z);
     186          35 :   return z;
     187             : }
     188             : 
     189             : GEN
     190          56 : intnumromb_bitprec(void *E, GEN (*f)(void *, GEN), GEN a,GEN b, long fl, long B)
     191             : {
     192          56 :   pari_sp av = avma;
     193             :   GEN z;
     194          56 :   switch(fl)
     195             :   {
     196           7 :     case 0: z = qrom3  (E, f, a, b, B); break;
     197          35 :     case 1: z = rombint(E, f, a, b, B); break;
     198           7 :     case 2: z = qromi  (E, f, a, b, B); break;
     199           7 :     case 3: z = qrom2  (E, f, a, b, B); break;
     200             :     default: pari_err_FLAG("intnumromb"); return NULL; /* LCOV_EXCL_LINE */
     201             :   }
     202          56 :   return gerepileupto(av, z);
     203             : }
     204             : GEN
     205           0 : intnumromb(void *E, GEN (*f)(void *, GEN), GEN a, GEN b, long flag, long prec)
     206           0 : { return intnumromb_bitprec(E,f,a,b,flag,prec2nbits(prec));}
     207             : GEN
     208          56 : intnumromb0_bitprec(GEN a, GEN b, GEN code, long flag, long bit)
     209          56 : { EXPR_WRAP(code, intnumromb_bitprec(EXPR_ARG, a, b, flag, bit)); }
     210             : 
     211             : /********************************************************************/
     212             : /**             NUMERICAL INTEGRATION (Gauss-Legendre)             **/
     213             : /********************************************************************/
     214             : GEN
     215          35 : intnumgaussinit(long n, long prec)
     216             : {
     217          35 :   pari_sp ltop = avma;
     218             :   GEN L, dp1, p1, p2, R, W;
     219          35 :   long prec0 = prec + EXTRAPRECWORD;
     220          35 :   long bitprec = prec2nbits(prec), i, d1;
     221          35 :   if (n <= 0) n = (long)(bitprec*0.2258);
     222          35 :   if (odd(n)) n++;
     223          35 :   if (n == 2) n = 4;
     224             :   /* n even >= 4, p1 is even */
     225          35 :   prec = nbits2prec(3*bitprec/2 + 32);
     226          35 :   L = pollegendre(n, 0); /* L_n = p1(x^2) */
     227          35 :   p1 = Q_remove_denom(RgX_deflate(L, 2), &dp1);
     228          35 :   d1 = vali(dp1);
     229          35 :   p2 = ZX_deriv(p1); /* L_n' = 2x p2(x^2) / 2^d1 */
     230          35 :   R = ZX_Uspensky(p1, gen_0, 1, 3*bitprec/2 + 32); /* positive roots of p1 */
     231          35 :   n >>= 1;
     232          35 :   W = cgetg(n+1, t_VEC);
     233         973 :   for (i = 1; i <= n; ++i)
     234             :   {
     235         938 :     GEN t, r2 = gel(R,i);
     236         938 :     if (typ(r2) != t_REAL) r2 = gtofp(r2, prec);
     237         938 :     gel(R,i) = sqrtr_abs(r2); /* positive root of L_n */
     238             :     /* 2 / (L'(r)^2(1-r^2)) =  2^(2d1 - 1) / (1-r2)r2 (p2(r2))^2 */
     239         938 :     t = mulrr(subrr(r2, sqrr(r2)), sqrr(poleval(p2, r2)));
     240         938 :     shiftr_inplace(t,1-2*d1);
     241         938 :     gel(W,i) = invr(t);
     242             :   }
     243          35 :   R = gprec_wtrunc(R,prec0);
     244          35 :   W = gprec_wtrunc(W,prec0);
     245          35 :   return gerepilecopy(ltop, mkvec2(R,W));
     246             : }
     247             : 
     248             : GEN
     249          56 : intnumgauss(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, GEN tab, long prec)
     250             : {
     251          56 :   pari_sp ltop = avma;
     252             :   GEN R, W, bma, bpa, S;
     253             :   long n, i;
     254          56 :   if (!tab)
     255           7 :     tab = intnumgaussinit(0,prec);
     256          49 :   else if (typ(tab) != t_INT)
     257             :   {
     258          35 :     if (typ(tab) != t_VEC || lg(tab) != 3)
     259           7 :       pari_err_TYPE("intnumgauss",tab);
     260             :   }
     261             :   else
     262          14 :     tab = intnumgaussinit(itos(tab),prec);
     263             : 
     264          49 :   R = gel(tab,1); n = lg(R)-1;
     265          49 :   W = gel(tab,2);
     266          49 :   a = gprec_w(a, prec+EXTRAPREC);
     267          49 :   b = gprec_w(b, prec+EXTRAPREC);
     268          49 :   bma = gmul2n(gsub(b,a), -1); /* (b-a)/2 */
     269          49 :   bpa = gadd(bma, a); /* (b+a)/2 */
     270          49 :   S = gen_0;
     271        1491 :   for (i = 1; i <= n; ++i)
     272             :   {
     273        1442 :     GEN r = gel(R,i);
     274        1442 :     GEN P = eval(E, gadd(bpa, gmul(bma, r)));
     275        1442 :     GEN M = eval(E, gsub(bpa, gmul(bma, r)));
     276        1442 :     S = gadd(S, gmul(gel(W,i), gadd(P,M)));
     277             :   }
     278          49 :   return gerepilecopy(ltop, gprec_wtrunc(gmul(bma,S), prec));
     279             : }
     280             : 
     281             : GEN
     282          56 : intnumgauss0(GEN a, GEN b, GEN code, GEN tab, long prec)
     283          56 : { EXPR_WRAP(code, intnumgauss(EXPR_ARG, a, b, tab, prec)); }
     284             : 
     285             : /********************************************************************/
     286             : /**                DOUBLE EXPONENTIAL INTEGRATION                  **/
     287             : /********************************************************************/
     288             : 
     289             : typedef struct _intdata {
     290             :   long eps;  /* bit accuracy of current precision */
     291             :   long l; /* table lengths */
     292             :   GEN tabx0; /* abscissa phi(0) for t = 0 */
     293             :   GEN tabw0; /* weight phi'(0) for t = 0 */
     294             :   GEN tabxp; /* table of abscissas phi(kh) for k > 0 */
     295             :   GEN tabwp; /* table of weights phi'(kh) for k > 0 */
     296             :   GEN tabxm; /* table of abscissas phi(kh) for k < 0, possibly empty */
     297             :   GEN tabwm; /* table of weights phi'(kh) for k < 0, possibly empty */
     298             :   GEN h; /* integration step */
     299             : } intdata;
     300             : 
     301             : static const long LGTAB = 8;
     302             : #define TABh(v) gel(v,1)
     303             : #define TABx0(v) gel(v,2)
     304             : #define TABw0(v) gel(v,3)
     305             : #define TABxp(v) gel(v,4)
     306             : #define TABwp(v) gel(v,5)
     307             : #define TABxm(v) gel(v,6)
     308             : #define TABwm(v) gel(v,7)
     309             : 
     310             : static int
     311       11088 : isinR(GEN z) { return is_real_t(typ(z)); }
     312             : static int
     313       10752 : isinC(GEN z)
     314       10752 : { return (typ(z) == t_COMPLEX)? isinR(gel(z,1)) && isinR(gel(z,2)): isinR(z); }
     315             : 
     316             : static int
     317        8820 : checktabsimp(GEN tab)
     318             : {
     319             :   long L, LN, LW;
     320        8820 :   if (!tab || typ(tab) != t_VEC) return 0;
     321        8820 :   if (lg(tab) != LGTAB) return 0;
     322        8820 :   if (typ(TABxp(tab)) != t_VEC) return 0;
     323        8820 :   if (typ(TABwp(tab)) != t_VEC) return 0;
     324        8820 :   if (typ(TABxm(tab)) != t_VEC) return 0;
     325        8820 :   if (typ(TABwm(tab)) != t_VEC) return 0;
     326        8820 :   L = lg(TABxp(tab)); if (lg(TABwp(tab)) != L) return 0;
     327        8820 :   LN = lg(TABxm(tab)); if (LN != 1 && LN != L) return 0;
     328        8820 :   LW = lg(TABwm(tab)); if (LW != 1 && LW != L) return 0;
     329        8820 :   return 1;
     330             : }
     331             : 
     332             : static int
     333         476 : checktabdoub(GEN tab)
     334             : {
     335             :   long L;
     336         476 :   if (typ(tab) != t_VEC) return 0;
     337         476 :   if (lg(tab) != LGTAB) return 0;
     338         476 :   L = lg(TABxp(tab));
     339         476 :   if (lg(TABwp(tab)) != L) return 0;
     340         476 :   if (lg(TABxm(tab)) != L) return 0;
     341         476 :   if (lg(TABwm(tab)) != L) return 0;
     342         476 :   return 1;
     343             : }
     344             : 
     345             : static int
     346        4445 : checktab(GEN tab)
     347             : {
     348        4445 :   if (typ(tab) != t_VEC) return 0;
     349        4445 :   if (lg(tab) != 3) return checktabsimp(tab);
     350          14 :   return checktabsimp(gel(tab,1))
     351           7 :       && checktabsimp(gel(tab,2));
     352             : }
     353             : 
     354             : /* the TUNE parameter is heuristic */
     355             : static void
     356         679 : intinit_start(intdata *D, long m, double TUNE, long prec)
     357             : {
     358         679 :   long l, n, bitprec = prec2nbits(prec);
     359         679 :   double d = bitprec*LOG10_2;
     360         679 :   GEN h, nh, pi = mppi(prec);
     361             : 
     362         679 :   n = (long)ceil(d*log(d) / TUNE); /* heuristic */
     363             :   /* nh ~ log(2npi/log(n)) */
     364         679 :   nh = logr_abs(divrr(mulur(2*n, pi), logr_abs(utor(n,prec))));
     365         679 :   h = divru(nh, n);
     366         679 :   if (m > 0) { h = gmul2n(h,-m); n <<= m; }
     367         679 :   D->h = h;
     368         679 :   D->eps = bitprec;
     369         679 :   D->l = l = n+1;
     370         679 :   D->tabxp = cgetg(l, t_VEC);
     371         679 :   D->tabwp = cgetg(l, t_VEC);
     372         679 :   D->tabxm = cgetg(l, t_VEC);
     373         679 :   D->tabwm = cgetg(l, t_VEC);
     374         679 : }
     375             : 
     376             : static GEN
     377         679 : intinit_end(intdata *D, long pnt, long mnt)
     378             : {
     379         679 :   GEN v = cgetg(LGTAB, t_VEC);
     380         679 :   if (pnt < 0) pari_err_DOMAIN("intnuminit","table length","<",gen_0,stoi(pnt));
     381         679 :   TABx0(v) = D->tabx0;
     382         679 :   TABw0(v) = D->tabw0;
     383         679 :   TABh(v) = D->h;
     384         679 :   TABxp(v) = D->tabxp; setlg(D->tabxp, pnt+1);
     385         679 :   TABwp(v) = D->tabwp; setlg(D->tabwp, pnt+1);
     386         679 :   TABxm(v) = D->tabxm; setlg(D->tabxm, mnt+1);
     387         679 :   TABwm(v) = D->tabwm; setlg(D->tabwm, mnt+1); return v;
     388             : }
     389             : 
     390             : /* divide by 2 in place */
     391             : static GEN
     392      257764 : divr2_ip(GEN x) { shiftr_inplace(x, -1); return x; }
     393             : 
     394             : /* phi(t)=tanh((Pi/2)sinh(t)): from -1 to 1, hence also from a to b compact
     395             :  * interval */
     396             : static GEN
     397         112 : inittanhsinh(long m, long prec)
     398             : {
     399         112 :   GEN et, ex, pi = mppi(prec);
     400         112 :   long k, nt = -1;
     401             :   intdata D;
     402             : 
     403         112 :   intinit_start(&D, m, 1.86, prec);
     404         112 :   D.tabx0 = real_0(prec);
     405         112 :   D.tabw0 = Pi2n(-1,prec);
     406         112 :   et = ex = mpexp(D.h);
     407       27160 :   for (k = 1; k < D.l; k++)
     408             :   {
     409             :     GEN xp, wp, ct, st, z;
     410             :     pari_sp av;
     411       27160 :     gel(D.tabxp,k) = cgetr(prec);
     412       27160 :     gel(D.tabwp,k) = cgetr(prec); av = avma;
     413       27160 :     ct = divr2_ip(addrr(et, invr(et))); /* ch(kh) */
     414       27160 :     st = subrr(et, ct); /* sh(kh) */
     415       27160 :     z = invr( addrs(mpexp(mulrr(pi, st)), 1) );
     416       27160 :     shiftr_inplace(z, 1);
     417       27160 :     xp = subsr(1, z);
     418       27160 :     wp = divr2_ip(mulrr(mulrr(pi,ct), mulrr(z, subsr(2, z))));
     419       27160 :     if (expo(wp) < -D.eps) { nt = k-1; break; }
     420       27160 :     affrr(xp, gel(D.tabxp,k));
     421       27160 :     if (absrnz_equal1(gel(D.tabxp,k))) { nt = k-1; break; }
     422       27048 :     affrr(wp, gel(D.tabwp,k)); et = gerepileuptoleaf(av, mulrr(et, ex));
     423             :   }
     424         112 :   return intinit_end(&D, nt, 0);
     425             : }
     426             : 
     427             : /* phi(t)=sinh(sinh(t)): from -oo to oo, slowly decreasing, at least
     428             :  * as 1/x^2. */
     429             : static GEN
     430          14 : initsinhsinh(long m, long prec)
     431             : {
     432             :   pari_sp av;
     433             :   GEN et, ct, st, ex;
     434          14 :   long k, nt = -1;
     435             :   intdata D;
     436             : 
     437          14 :   intinit_start(&D, m, 0.666, prec);
     438          14 :   D.tabx0 = real_0(prec);
     439          14 :   D.tabw0 = real_1(prec);
     440          14 :   et = ex = mpexp(D.h);
     441        8184 :   for (k = 1; k < D.l; k++)
     442             :   {
     443             :     GEN xp, wp, ext, exu;
     444        8184 :     gel(D.tabxp,k) = cgetr(prec);
     445        8184 :     gel(D.tabwp,k) = cgetr(prec); av = avma;
     446        8184 :     ct = divr2_ip(addrr(et, invr(et)));
     447        8184 :     st = subrr(et, ct);
     448        8184 :     ext = mpexp(st);
     449        8184 :     exu = invr(ext);
     450        8184 :     xp = divr2_ip(subrr(ext, exu));
     451        8184 :     wp = divr2_ip(mulrr(ct, addrr(ext, exu)));
     452        8184 :     if (expo(wp) - 2*expo(xp) < -D.eps) { nt = k-1; break; }
     453        8170 :     affrr(xp, gel(D.tabxp,k));
     454        8170 :     affrr(wp, gel(D.tabwp,k)); et = gerepileuptoleaf(av, mulrr(et, ex));
     455             :   }
     456          14 :   return intinit_end(&D, nt, 0);
     457             : }
     458             : 
     459             : /* phi(t)=2sinh(t): from -oo to oo, exponentially decreasing as exp(-x) */
     460             : static GEN
     461         126 : initsinh(long m, long prec)
     462             : {
     463             :   pari_sp av;
     464             :   GEN et, ex, eti, xp, wp;
     465         126 :   long k, nt = -1;
     466             :   intdata D;
     467             : 
     468         126 :   intinit_start(&D, m, 1.0, prec);
     469         126 :   D.tabx0 = real_0(prec);
     470         126 :   D.tabw0 = real2n(1, prec);
     471         126 :   et = ex = mpexp(D.h);
     472       38136 :   for (k = 1; k < D.l; k++)
     473             :   {
     474       38136 :     gel(D.tabxp,k) = cgetr(prec);
     475       38136 :     gel(D.tabwp,k) = cgetr(prec); av = avma;
     476       38136 :     eti = invr(et);
     477       38136 :     xp = subrr(et, eti);
     478       38136 :     wp = addrr(et, eti);
     479       38136 :     if (cmprs(xp, (long)(LOG2*(expo(wp)+D.eps) + 1)) > 0) { nt = k-1; break; }
     480       38010 :     affrr(xp, gel(D.tabxp,k));
     481       38010 :     affrr(wp, gel(D.tabwp,k)); et = gerepileuptoleaf(av, mulrr(et, ex));
     482             :   }
     483         126 :   return intinit_end(&D, nt, 0);
     484             : }
     485             : 
     486             : /* phi(t)=exp(2sinh(t)): from 0 to oo, slowly decreasing at least as 1/x^2 */
     487             : static GEN
     488         238 : initexpsinh(long m, long prec)
     489             : {
     490             :   GEN et, ex;
     491         238 :   long k, nt = -1;
     492             :   intdata D;
     493             : 
     494         238 :   intinit_start(&D, m, 1.05, prec);
     495         238 :   D.tabx0 = real_1(prec);
     496         238 :   D.tabw0 = real2n(1, prec);
     497         238 :   ex = mpexp(D.h);
     498         238 :   et = real_1(prec);
     499      120482 :   for (k = 1; k < D.l; k++)
     500             :   {
     501             :     GEN t, eti, xp;
     502      120482 :     et = mulrr(et, ex);
     503      120482 :     eti = invr(et); t = addrr(et, eti);
     504      120482 :     xp = mpexp(subrr(et, eti));
     505      120482 :     gel(D.tabxp,k) = xp;
     506      120482 :     gel(D.tabwp,k) = mulrr(xp, t);
     507      120482 :     gel(D.tabxm,k) = invr(xp);
     508      120482 :     gel(D.tabwm,k) = mulrr(gel(D.tabxm,k), t);
     509      120482 :     if (expo(gel(D.tabxm,k)) < -D.eps) { nt = k-1; break; }
     510             :   }
     511         238 :   return intinit_end(&D, nt, nt);
     512             : }
     513             : 
     514             : /* phi(t)=exp(t-exp(-t)) : from 0 to \infty, exponentially decreasing. */
     515             : static GEN
     516          84 : initexpexp(long m, long prec)
     517             : {
     518             :   pari_sp av;
     519             :   GEN et, ex;
     520          84 :   long k, nt = -1;
     521             :   intdata D;
     522             : 
     523          84 :   intinit_start(&D, m, 1.76, prec);
     524          84 :   D.tabx0 = mpexp(real_m1(prec));
     525          84 :   D.tabw0 = gmul2n(D.tabx0, 1);
     526          84 :   et = ex = mpexp(negr(D.h));
     527       36112 :   for (k = 1; k < D.l; k++)
     528             :   {
     529             :     GEN xp, xm, wp, wm, eti, kh;
     530       36112 :     gel(D.tabxp,k) = cgetr(prec);
     531       36112 :     gel(D.tabwp,k) = cgetr(prec);
     532       36112 :     gel(D.tabxm,k) = cgetr(prec);
     533       36112 :     gel(D.tabwm,k) = cgetr(prec); av = avma;
     534       36112 :     eti = invr(et); kh = mulur(k,D.h);
     535       36112 :     xp = mpexp(subrr(kh, et));
     536       36112 :     xm = mpexp(negr(addrr(kh, eti)));
     537       36112 :     wp = mulrr(xp, addsr(1, et));
     538       36112 :     wm = mulrr(xm, addsr(1, eti));
     539       36112 :     if (expo(xm) < -D.eps && cmprs(xp, (long)(LOG2*(expo(wp)+D.eps) + 1)) > 0) { nt = k-1; break; }
     540       36028 :     affrr(xp, gel(D.tabxp,k));
     541       36028 :     affrr(wp, gel(D.tabwp,k));
     542       36028 :     affrr(xm, gel(D.tabxm,k));
     543       36028 :     affrr(wm, gel(D.tabwm,k)); et = gerepileuptoleaf(av, mulrr(et, ex));
     544             :   }
     545          84 :   return intinit_end(&D, nt, nt);
     546             : }
     547             : 
     548             : /* phi(t)=(Pi/h)*t/(1-exp(-sinh(t))) from 0 to oo, sine oscillation */
     549             : static GEN
     550         105 : initnumsine(long m, long prec)
     551             : {
     552             :   pari_sp av;
     553         105 :   GEN invh, et, eti, ex, pi = mppi(prec);
     554         105 :   long exh, k, nt = -1;
     555             :   intdata D;
     556             : 
     557         105 :   intinit_start(&D, m, 0.666, prec);
     558         105 :   invh = invr(D.h);
     559         105 :   D.tabx0 = mulrr(pi, invh);
     560         105 :   D.tabw0 = gmul2n(D.tabx0,-1);
     561         105 :   exh = expo(invh); /*  expo(1/h) */
     562         105 :   et = ex = mpexp(D.h);
     563       89446 :   for (k = 1; k < D.l; k++)
     564             :   {
     565             :     GEN xp,xm, wp,wm, ct,st, extp,extp1,extp2, extm,extm1,extm2, kct, kpi;
     566       89446 :     gel(D.tabxp,k) = cgetr(prec);
     567       89446 :     gel(D.tabwp,k) = cgetr(prec);
     568       89446 :     gel(D.tabxm,k) = cgetr(prec);
     569       89446 :     gel(D.tabwm,k) = cgetr(prec); av = avma;
     570       89446 :     eti = invr(et); /* exp(-kh) */
     571       89446 :     ct = divr2_ip(addrr(et, eti)); /* ch(kh) */
     572       89446 :     st = divr2_ip(subrr(et, eti)); /* sh(kh) */
     573       89446 :     extp = mpexp(st);  extp1 = subsr(1, extp);
     574       89446 :     extp2 = invr(extp1); /* 1/(1-exp(sh(kh))) */
     575       89446 :     extm = invr(extp); extm1 = subsr(1, extm);
     576       89446 :     extm2 = invr(extm1);/* 1/(1-exp(sh(-kh))) */
     577       89446 :     kpi = mulur(k, pi);
     578       89446 :     kct = mulur(k, ct);
     579       89446 :     extm1 = mulrr(extm1, invh);
     580       89446 :     extp1 = mulrr(extp1, invh);
     581       89446 :     xp = mulrr(kpi, extm2); /* phi(kh) */
     582       89446 :     wp = mulrr(subrr(extm1, mulrr(kct, extm)), mulrr(pi, sqrr(extm2)));
     583       89446 :     xm = mulrr(negr(kpi), extp2); /* phi(-kh) */
     584       89446 :     wm = mulrr(addrr(extp1, mulrr(kct, extp)), mulrr(pi, sqrr(extp2)));
     585       89446 :     if (expo(wm) < -D.eps && expo(extm) + exh + expu(10 * k) < -D.eps) { nt = k-1; break; }
     586       89341 :     affrr(xp, gel(D.tabxp,k));
     587       89341 :     affrr(wp, gel(D.tabwp,k));
     588       89341 :     affrr(xm, gel(D.tabxm,k));
     589       89341 :     affrr(wm, gel(D.tabwm,k)); et = gerepileuptoleaf(av, mulrr(et, ex));
     590             :   }
     591         105 :   return intinit_end(&D, nt, nt);
     592             : }
     593             : 
     594             : /* End of initialization functions. These functions can be executed once
     595             :  * and for all for a given accuracy and type of integral ([a,b], [a,oo[ or
     596             :  * ]-oo,a], ]-oo,oo[) */
     597             : 
     598             : /* The numbers below can be changed, but NOT the ordering */
     599             : enum {
     600             :   f_REG    = 0, /* regular function */
     601             :   f_SING   = 1, /* algebraic singularity */
     602             :   f_YSLOW  = 2, /* +\infty, slowly decreasing, at least x^(-2)  */
     603             :   f_YVSLO  = 3, /* +\infty, very slowly decreasing, worse than x^(-2) */
     604             :   f_YFAST  = 4, /* +\infty, exponentially decreasing */
     605             :   f_YOSCS  = 5, /* +\infty, sine oscillating */
     606             :   f_YOSCC  = 6  /* +\infty, cosine oscillating */
     607             : };
     608             : /* is finite ? */
     609             : static int
     610         959 : is_fin_f(long c) { return c == f_REG || c == f_SING; }
     611             : /* is oscillatory ? */
     612             : static int
     613         140 : is_osc(long c) { long a = labs(c); return a == f_YOSCC|| a == f_YOSCS; }
     614             : 
     615             : /* All inner functions such as intn, etc... must be called with a
     616             :  * valid 'tab' table. The wrapper intnum provides a higher level interface */
     617             : 
     618             : /* compute \int_a^b f(t)dt with [a,b] compact and f nonsingular. */
     619             : static GEN
     620        3745 : intn(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, GEN tab)
     621             : {
     622             :   GEN tabx0, tabw0, tabxp, tabwp;
     623             :   GEN bpa, bma, bmb, S;
     624             :   long i;
     625        3745 :   pari_sp ltop = avma, av;
     626             : 
     627        3745 :   if (!checktabsimp(tab)) pari_err_TYPE("intnum",tab);
     628        3745 :   tabx0 = TABx0(tab); tabw0 = TABw0(tab);
     629        3745 :   tabxp = TABxp(tab); tabwp = TABwp(tab);
     630        3745 :   bpa = gmul2n(gadd(b, a), -1); /* (b+a)/2 */
     631        3745 :   bma = gsub(bpa, a); /* (b-a)/2 */
     632        3745 :   av = avma;
     633        3745 :   bmb = gmul(bma, tabx0); /* (b-a)/2 phi(0) */
     634             :   /* phi'(0) f( (b+a)/2 + (b-a)/2 * phi(0) ) */
     635        3745 :   S = gmul(tabw0, eval(E, gadd(bpa, bmb)));
     636      982639 :   for (i = lg(tabxp)-1; i > 0; i--)
     637             :   {
     638             :     GEN SP, SM;
     639      978894 :     bmb = gmul(bma, gel(tabxp,i));
     640      978894 :     SP = eval(E, gsub(bpa, bmb));
     641      978894 :     SM = eval(E, gadd(bpa, bmb));
     642      978894 :     S = gadd(S, gmul(gel(tabwp,i), gadd(SP, SM)));
     643      978894 :     if ((i & 0x7f) == 1) S = gerepileupto(av, S);
     644             :   }
     645        3745 :   return gerepileupto(ltop, gmul(S, gmul(bma, TABh(tab))));
     646             : }
     647             : 
     648             : /* compute \int_a^b f(t)dt with [a,b] compact, possible singularity with
     649             :  * exponent a[2] at lower extremity, b regular. Use tanh(sinh(t)). */
     650             : static GEN
     651          42 : intnsing(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, GEN tab, long prec)
     652             : {
     653             :   GEN tabx0, tabw0, tabxp, tabwp, ea, ba, S;
     654             :   long i;
     655          42 :   pari_sp ltop = avma, av;
     656             : 
     657          42 :   if (!checktabsimp(tab)) pari_err_TYPE("intnum",tab);
     658          42 :   tabx0 = TABx0(tab); tabw0 = TABw0(tab);
     659          42 :   tabxp = TABxp(tab); tabwp = TABwp(tab);
     660          42 :   ea = ginv(gaddsg(1, gel(a,2)));
     661          42 :   a = gel(a,1);
     662          42 :   ba = gdiv(gsub(b, a), gpow(gen_2, ea, prec));
     663          42 :   av = avma;
     664          42 :   S = gmul(gmul(tabw0, ba), eval(E, gadd(gmul(ba, addsr(1, tabx0)), a)));
     665       11046 :   for (i = lg(tabxp)-1; i > 0; i--)
     666             :   {
     667       11004 :     GEN p = addsr(1, gel(tabxp,i));
     668       11004 :     GEN m = subsr(1, gel(tabxp,i));
     669       11004 :     GEN bp = gmul(ba, gpow(p, ea, prec));
     670       11004 :     GEN bm = gmul(ba, gpow(m, ea, prec));
     671       11004 :     GEN SP = gmul(gdiv(bp, p), eval(E, gadd(bp, a)));
     672       11004 :     GEN SM = gmul(gdiv(bm, m), eval(E, gadd(bm, a)));
     673       11004 :     S = gadd(S, gmul(gel(tabwp,i), gadd(SP, SM)));
     674       11004 :     if ((i & 0x7f) == 1) S = gerepileupto(av, S);
     675             :   }
     676          42 :   return gerepileupto(ltop, gmul(gmul(S, TABh(tab)), ea));
     677             : }
     678             : 
     679      170450 : static GEN id(GEN x) { return x; }
     680             : 
     681             : /* compute  \int_a^oo f(t)dt if si>0 or \int_{-oo}^a f(t)dt if si<0$.
     682             :  * Use exp(2sinh(t)) for slowly decreasing functions, exp(1+t-exp(-t)) for
     683             :  * exponentially decreasing functions, and (pi/h)t/(1-exp(-sinh(t))) for
     684             :  * oscillating functions. */
     685             : static GEN
     686         476 : intninfpm(void *E, GEN (*eval)(void*, GEN), GEN a, long sb, GEN tab)
     687             : {
     688             :   GEN tabx0, tabw0, tabxp, tabwp, tabxm, tabwm;
     689             :   GEN S;
     690             :   long L, i;
     691         476 :   pari_sp av = avma;
     692             : 
     693         476 :   if (!checktabdoub(tab)) pari_err_TYPE("intnum",tab);
     694         476 :   tabx0 = TABx0(tab); tabw0 = TABw0(tab);
     695         476 :   tabxp = TABxp(tab); tabwp = TABwp(tab); L = lg(tabxp);
     696         476 :   tabxm = TABxm(tab); tabwm = TABwm(tab);
     697         476 :   if (gequal0(a))
     698             :   {
     699         238 :     GEN (*NEG)(GEN) = sb > 0? id: gneg;
     700         238 :     S = gmul(tabw0, eval(E, NEG(tabx0)));
     701      128814 :     for (i = 1; i < L; i++)
     702             :     {
     703      128576 :       GEN SP = eval(E, NEG(gel(tabxp,i)));
     704      128576 :       GEN SM = eval(E, NEG(gel(tabxm,i)));
     705      128576 :       S = gadd(S, gadd(gmul(gel(tabwp,i), SP), gmul(gel(tabwm,i), SM)));
     706      128576 :       if ((i & 0x7f) == 1) S = gerepileupto(av, S);
     707             :     }
     708             :   }
     709         238 :   else if (gexpo(a) <= 0 || is_osc(sb))
     710         105 :   { /* a small */
     711         105 :     GEN (*ADD)(GEN,GEN) = sb > 0? gadd: gsub;
     712         105 :     S = gmul(tabw0, eval(E, ADD(a, tabx0)));
     713       61460 :     for (i = 1; i < L; i++)
     714             :     {
     715       61355 :       GEN SP = eval(E, ADD(a, gel(tabxp,i)));
     716       61355 :       GEN SM = eval(E, ADD(a, gel(tabxm,i)));
     717       61355 :       S = gadd(S, gadd(gmul(gel(tabwp,i), SP), gmul(gel(tabwm,i), SM)));
     718       61355 :       if ((i & 0x7f) == 1) S = gerepileupto(av, S);
     719             :     }
     720             :   }
     721             :   else
     722             :   { /* a large, |a|*\int_sgn(a)^{oo} f(|a|*x)dx (sb > 0)*/
     723         133 :     GEN (*ADD)(long,GEN) = sb > 0? addsr: subsr;
     724         133 :     long sa = gsigne(a);
     725         133 :     GEN A = sa > 0? a: gneg(a);
     726         133 :     pari_sp av2 = avma;
     727         133 :     S = gmul(tabw0, eval(E, gmul(A, ADD(sa, tabx0))));
     728       98604 :     for (i = 1; i < L; i++)
     729             :     {
     730       98471 :       GEN SP = eval(E, gmul(A, ADD(sa, gel(tabxp,i))));
     731       98471 :       GEN SM = eval(E, gmul(A, ADD(sa, gel(tabxm,i))));
     732       98471 :       S = gadd(S, gadd(gmul(gel(tabwp,i), SP), gmul(gel(tabwm,i), SM)));
     733       98471 :       if ((i & 0x7f) == 1) S = gerepileupto(av2, S);
     734             :     }
     735         133 :     S = gmul(S,A);
     736             :   }
     737         476 :   return gerepileupto(av, gmul(S, TABh(tab)));
     738             : }
     739             : 
     740             : /* Compute  \int_{-oo}^oo f(t)dt
     741             :  * use sinh(sinh(t)) for slowly decreasing functions and sinh(t) for
     742             :  * exponentially decreasing functions.
     743             :  * HACK: in case TABwm(tab) contains something, assume function to be integrated
     744             :  * satisfies f(-x) = conj(f(x)).
     745             :  */
     746             : static GEN
     747         581 : intninfinf(void *E, GEN (*eval)(void*, GEN), GEN tab)
     748             : {
     749             :   GEN tabx0, tabw0, tabxp, tabwp, tabwm;
     750             :   GEN S;
     751             :   long L, i, spf;
     752         581 :   pari_sp ltop = avma;
     753             : 
     754         581 :   if (!checktabsimp(tab)) pari_err_TYPE("intnum",tab);
     755         581 :   tabx0 = TABx0(tab); tabw0 = TABw0(tab);
     756         581 :   tabxp = TABxp(tab); tabwp = TABwp(tab); L = lg(tabxp);
     757         581 :   tabwm = TABwm(tab);
     758         581 :   spf = (lg(tabwm) == lg(tabwp));
     759         581 :   S = gmul(tabw0, eval(E, tabx0));
     760         581 :   if (spf) S = gmul2n(real_i(S), -1);
     761      176099 :   for (i = L-1; i > 0; i--)
     762             :   {
     763      175518 :     GEN SP = eval(E, gel(tabxp,i));
     764      175518 :     if (spf)
     765      170044 :       S = gadd(S, real_i(gmul(gel(tabwp,i), SP)));
     766             :     else
     767             :     {
     768        5474 :       GEN SM = eval(E, negr(gel(tabxp,i)));
     769        5474 :       S = gadd(S, gmul(gel(tabwp,i), gadd(SP,SM)));
     770             :     }
     771      175518 :     if ((i & 0x7f) == 1) S = gerepileupto(ltop, S);
     772             :   }
     773         581 :   if (spf) S = gmul2n(S,1);
     774         581 :   return gerepileupto(ltop, gmul(S, TABh(tab)));
     775             : }
     776             : 
     777             : /* general num integration routine int_a^b f(t)dt, where a and b are as follows:
     778             :  - a scalar : the scalar, no singularity worse than logarithmic at a.
     779             :  - [a, e] : the scalar a, singularity exponent -1 < e <= 0.
     780             :  - +oo: slowly decreasing function (at least O(t^-2))
     781             :  - [[+oo], a], a nonnegative real : +oo, function behaving like exp(-a|t|)
     782             :  - [[+oo], e], e < -1 : +oo, function behaving like t^e
     783             :  - [[+oo], a*I], a > 0 real : +oo, function behaving like cos(at)
     784             :  - [[+oo], a*I], a < 0 real : +oo, function behaving like sin(at)
     785             :  and similarly at -oo */
     786             : static GEN
     787        1939 : f_getycplx(GEN a, long prec)
     788             : {
     789             :   long s;
     790             :   GEN tmp, a2R, a2I;
     791             : 
     792        1939 :   if (lg(a) == 2 || gequal0(gel(a,2))) return gen_1;
     793        1897 :   a2R = real_i(gel(a,2));
     794        1897 :   a2I = imag_i(gel(a,2));
     795        1897 :   s = gsigne(a2I); if (s < 0) a2I = gneg(a2I);
     796        1897 :   tmp = s ? ginv(a2I) : ginv(a2R);
     797        1897 :   if (gprecision(tmp) < prec) tmp = gprec_w(tmp, prec);
     798        1897 :   return tmp;
     799             : }
     800             : 
     801             : static void
     802          21 : err_code(GEN a, const char *name)
     803             : {
     804          21 :   char *s = stack_sprintf("intnum [incorrect %s]", name);
     805          21 :   pari_err_TYPE(s, a);
     806           0 : }
     807             : 
     808             : /* a = [[+/-oo], alpha]*/
     809             : static long
     810        2562 : code_aux(GEN a, const char *name)
     811             : {
     812        2562 :   GEN re, im, alpha = gel(a,2);
     813             :   long s;
     814        2562 :   if (!isinC(alpha)) err_code(a, name);
     815        2555 :   re = real_i(alpha);
     816        2555 :   im = imag_i(alpha);
     817        2555 :   s = gsigne(im);
     818        2555 :   if (s)
     819             :   {
     820         252 :     if(!gequal0(re))
     821           7 :       pari_warn(warner,"real(z)*imag(z)!=0 in endpoint code, real(z) ignored");
     822         252 :     return s > 0 ? f_YOSCC : f_YOSCS;
     823             :   }
     824        2303 :   if (gequal0(re) || gcmpgs(re, -2)<=0) return f_YSLOW;
     825        2065 :   if (gsigne(re) > 0) return f_YFAST;
     826         231 :   if (gcmpgs(re, -1) >= 0) err_code(a, name);
     827         231 :   return f_YVSLO;
     828             : }
     829             : 
     830             : static long
     831       10913 : transcode(GEN a, const char *name)
     832             : {
     833             :   GEN a1, a2;
     834       10913 :   switch(typ(a))
     835             :   {
     836        2667 :     case t_VEC: break;
     837             :     case t_INFINITY:
     838         140 :       return inf_get_sign(a) == 1 ? f_YSLOW: -f_YSLOW;
     839        8106 :     default: if (!isinC(a)) err_code(a,name);
     840        8099 :       return f_REG;
     841             :   }
     842        2667 :   switch(lg(a))
     843             :   {
     844          14 :     case 2: return gsigne(gel(a,1)) > 0 ? f_YSLOW : -f_YSLOW;
     845        2646 :     case 3: break;
     846           7 :     default: err_code(a,name);
     847             :   }
     848        2646 :   a1 = gel(a,1);
     849        2646 :   a2 = gel(a,2);
     850        2646 :   switch(typ(a1))
     851             :   {
     852             :     case t_VEC:
     853          14 :       if (lg(a1) != 2) err_code(a,name);
     854          14 :       return gsigne(gel(a1,1)) * code_aux(a, name);
     855             :     case t_INFINITY:
     856        2548 :       return inf_get_sign(a1) * code_aux(a, name);
     857             :     default:
     858          84 :       if (!isinC(a1) || !isinR(a2) || gcmpgs(a2, -1) <= 0) err_code(a,name);
     859          84 :       return gsigne(a2) < 0 ? f_SING : f_REG;
     860             :   }
     861             : }
     862             : 
     863             : /* computes the necessary tabs, knowing a, b and m */
     864             : static GEN
     865         343 : homtab(GEN tab, GEN k)
     866             : {
     867             :   GEN z;
     868         343 :   if (gequal0(k) || gequal(k, gen_1)) return tab;
     869         217 :   if (gsigne(k) < 0) k = gneg(k);
     870         217 :   z = cgetg(LGTAB, t_VEC);
     871         217 :   TABx0(z) = gmul(TABx0(tab), k);
     872         217 :   TABw0(z) = gmul(TABw0(tab), k);
     873         217 :   TABxp(z) = gmul(TABxp(tab), k);
     874         217 :   TABwp(z) = gmul(TABwp(tab), k);
     875         217 :   TABxm(z) = gmul(TABxm(tab), k);
     876         217 :   TABwm(z) = gmul(TABwm(tab), k);
     877         217 :   TABh(z) = rcopy(TABh(tab)); return z;
     878             : }
     879             : 
     880             : static GEN
     881         238 : expvec(GEN v, GEN ea, long prec)
     882             : {
     883         238 :   long lv = lg(v), i;
     884         238 :   GEN z = cgetg(lv, t_VEC);
     885         238 :   for (i = 1; i < lv; i++) gel(z,i) = gpow(gel(v,i),ea,prec);
     886         238 :   return z;
     887             : }
     888             : 
     889             : static GEN
     890      137561 : expscalpr(GEN vnew, GEN xold, GEN wold, GEN ea)
     891             : {
     892      137561 :   pari_sp av = avma;
     893      137561 :   return gerepileupto(av, gdiv(gmul(gmul(vnew, wold), ea), xold));
     894             : }
     895             : static GEN
     896         238 : expvecpr(GEN vnew, GEN xold, GEN wold, GEN ea)
     897             : {
     898         238 :   long lv = lg(vnew), i;
     899         238 :   GEN z = cgetg(lv, t_VEC);
     900      137680 :   for (i = 1; i < lv; i++)
     901      137442 :     gel(z,i) = expscalpr(gel(vnew,i), gel(xold,i), gel(wold,i), ea);
     902         238 :   return z;
     903             : }
     904             : 
     905             : /* here k < -1 */
     906             : static GEN
     907         119 : exptab(GEN tab, GEN k, long prec)
     908             : {
     909             :   GEN v, ea;
     910             : 
     911         119 :   if (gcmpgs(k, -2) <= 0) return tab;
     912         119 :   ea = ginv(gsubsg(-1, k));
     913         119 :   v = cgetg(LGTAB, t_VEC);
     914         119 :   TABx0(v) = gpow(TABx0(tab), ea, prec);
     915         119 :   TABw0(v) = expscalpr(TABx0(v), TABx0(tab), TABw0(tab), ea);
     916         119 :   TABxp(v) = expvec(TABxp(tab), ea, prec);
     917         119 :   TABwp(v) = expvecpr(TABxp(v), TABxp(tab), TABwp(tab), ea);
     918         119 :   TABxm(v) = expvec(TABxm(tab), ea, prec);
     919         119 :   TABwm(v) = expvecpr(TABxm(v), TABxm(tab), TABwm(tab), ea);
     920         119 :   TABh(v) = rcopy(TABh(tab));
     921         119 :   return v;
     922             : }
     923             : 
     924             : static GEN
     925         343 : init_fin(GEN b, long codeb, long m, long l, long prec)
     926             : {
     927         343 :   switch(labs(codeb))
     928             :   {
     929             :     case f_REG:
     930          70 :     case f_SING:  return inittanhsinh(m,l);
     931         112 :     case f_YSLOW: return initexpsinh(m,l);
     932          70 :     case f_YVSLO: return exptab(initexpsinh(m,l), gel(b,2), prec);
     933          42 :     case f_YFAST: return homtab(initexpexp(m,l), f_getycplx(b,l));
     934             :     /* f_YOSCS, f_YOSCC */
     935          49 :     default: return homtab(initnumsine(m,l),f_getycplx(b,l));
     936             :   }
     937             : }
     938             : 
     939             : static GEN
     940         623 : intnuminit_i(GEN a, GEN b, long m, long prec)
     941             : {
     942             :   long codea, codeb, l;
     943             :   GEN T, kma, kmb, tmp;
     944             : 
     945         623 :   if (m > 30) pari_err_OVERFLOW("intnuminit [m]");
     946         623 :   if (m < 0) pari_err_DOMAIN("intnuminit", "m", "<", gen_0, stoi(m));
     947         616 :   l = prec+EXTRAPREC;
     948         616 :   codea = transcode(a, "a");
     949         609 :   codeb = transcode(b, "b");
     950         595 :   if (labs(codea) > labs(codeb)) { swap(a, b); lswap(codea, codeb); }
     951         595 :   if (codea == f_REG)
     952             :   {
     953         322 :     T = init_fin(b, codeb, m,l,prec);
     954         322 :     switch(labs(codeb))
     955             :     {
     956          42 :       case f_YOSCS: if (gequal0(a)) break;
     957           7 :       case f_YOSCC: T = mkvec2(inittanhsinh(m,l), T);
     958             :     }
     959         322 :     return T;
     960             :   }
     961         273 :   if (codea == f_SING)
     962             :   {
     963          21 :     T = init_fin(b,codeb, m,l,prec);
     964          21 :     T = mkvec2(inittanhsinh(m,l), T);
     965          21 :     return T;
     966             :   }
     967             :   /* now a and b are infinite */
     968         252 :   if (codea * codeb > 0) return gen_0;
     969         238 :   kma = f_getycplx(a,l); codea = labs(codea);
     970         238 :   kmb = f_getycplx(b,l); codeb = labs(codeb);
     971         238 :   if (codea == f_YSLOW && codeb == f_YSLOW) return initsinhsinh(m, l);
     972         224 :   if (codea == f_YFAST && codeb == f_YFAST && gequal(kma, kmb))
     973         126 :     return homtab(initsinh(m,l), kmb);
     974          98 :   T = cgetg(3, t_VEC);
     975          98 :   switch (codea)
     976             :   {
     977             :     case f_YSLOW:
     978             :     case f_YVSLO:
     979          56 :       tmp = initexpsinh(m,l);
     980          56 :       gel(T,1) = codea == f_YSLOW? tmp: exptab(tmp, gel(a,2), prec);
     981          56 :       switch (codeb)
     982             :       {
     983          14 :         case f_YVSLO: gel(T,2) = exptab(tmp, gel(b,2), prec); return T;
     984          21 :         case f_YFAST: gel(T,2) = homtab(initexpexp(m,l), kmb); return T;
     985             :         case f_YOSCS:
     986          21 :         case f_YOSCC: gel(T,2) = homtab(initnumsine(m,l), kmb); return T;
     987             :       }
     988           0 :       break;
     989             :     case f_YFAST:
     990          21 :       tmp = initexpexp(m, l);
     991          21 :       gel(T,1) = homtab(tmp, kma);
     992          21 :       switch (codeb)
     993             :       {
     994           7 :         case f_YFAST: gel(T,2) = homtab(tmp, kmb); return T;
     995             :         case f_YOSCS:
     996          14 :         case f_YOSCC: gel(T,2) = homtab(initnumsine(m, l), kmb); return T;
     997             :       }
     998             :     case f_YOSCS: case f_YOSCC:
     999          21 :       tmp = initnumsine(m, l);
    1000          21 :       gel(T,1) = homtab(tmp,kma);
    1001          21 :       if (codea == f_YOSCC && codeb == f_YOSCC && !gequal(kma, kmb))
    1002          14 :         gel(T,2) = mkvec2(inittanhsinh(m,l), homtab(tmp,kmb));
    1003             :       else
    1004           7 :         gel(T,2) = homtab(tmp,kmb);
    1005          21 :       return T;
    1006             :   }
    1007             :   return gen_0; /* LCOV_EXCL_LINE */
    1008             : }
    1009             : GEN
    1010         497 : intnuminit(GEN a, GEN b, long m, long prec)
    1011             : {
    1012         497 :   pari_sp av = avma;
    1013         497 :   return gerepilecopy(av, intnuminit_i(a,b,m,prec));
    1014             : }
    1015             : 
    1016             : static GEN
    1017        4718 : intnuminit0(GEN a, GEN b, GEN tab, long prec)
    1018             : {
    1019             :   long m;
    1020        4718 :   if (!tab) m = 0;
    1021        4452 :   else if (typ(tab) != t_INT)
    1022             :   {
    1023        4445 :     if (!checktab(tab)) pari_err_TYPE("intnuminit0",tab);
    1024        4445 :     return tab;
    1025             :   }
    1026             :   else
    1027           7 :     m = itos(tab);
    1028         273 :   return intnuminit(a, b, m, prec);
    1029             : }
    1030             : 
    1031             : /* Assigns the values of the function weighted by w[k] at quadrature points x[k]
    1032             :  * [replacing the weights]. Return the index of the last non-zero coeff */
    1033             : static long
    1034         252 : weight(void *E, GEN (*eval)(void *, GEN), GEN x, GEN w)
    1035             : {
    1036         252 :   long k, l = lg(x);
    1037         252 :   for (k = 1; k < l; k++) gel(w,k) = gmul(gel(w,k), eval(E, gel(x,k)));
    1038         252 :   k--; while (k >= 1) if (!gequal0(gel(w,k--))) break;
    1039         252 :   return k;
    1040             : }
    1041             : /* compute the necessary tabs, weights multiplied by f(t) */
    1042             : static GEN
    1043         126 : intfuncinit_i(void *E, GEN (*eval)(void*, GEN), GEN tab)
    1044             : {
    1045         126 :   GEN tabxp = TABxp(tab), tabwp = TABwp(tab);
    1046         126 :   GEN tabxm = TABxm(tab), tabwm = TABwm(tab);
    1047         126 :   long L = weight(E, eval, tabxp, tabwp), L0 = lg(tabxp);
    1048             : 
    1049         126 :   TABw0(tab) = gmul(TABw0(tab), eval(E, TABx0(tab)));
    1050         126 :   if (lg(tabxm) > 1)
    1051           0 :     (void)weight(E, eval, tabxm, tabwm);
    1052             :   else
    1053             :   {
    1054         126 :     tabxm = gneg(tabxp);
    1055         126 :     tabwm = leafcopy(tabwp);
    1056         126 :     L = minss(L, weight(E, eval, tabxm, tabwm));
    1057         126 :     TABxm(tab) = tabxm;
    1058         126 :     TABwm(tab) = tabwm;
    1059             :   }
    1060         126 :   if (L < L0)
    1061             :   { /* catch up functions whose growth at oo was not adequately described */
    1062         126 :     setlg(tabxp, L+1);
    1063         126 :     setlg(tabwp, L+1);
    1064         126 :     if (lg(tabxm) > 1) { setlg(tabxm, L+1); setlg(tabwm, L+1); }
    1065             :   }
    1066         126 :   return tab;
    1067             : }
    1068             : 
    1069             : GEN
    1070         126 : intfuncinit(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, long m, long prec)
    1071             : {
    1072         126 :   pari_sp av = avma;
    1073         126 :   GEN tab = intnuminit_i(a, b, m, prec);
    1074             : 
    1075         126 :   if (lg(tab) == 3)
    1076           0 :     pari_err_IMPL("intfuncinit with hard endpoint behaviour");
    1077         252 :   if (is_fin_f(transcode(a,"intfuncinit")) ||
    1078         126 :       is_fin_f(transcode(b,"intfuncinit")))
    1079           0 :     pari_err_IMPL("intfuncinit with finite endpoints");
    1080         126 :   return gerepilecopy(av, intfuncinit_i(E, eval, tab));
    1081             : }
    1082             : 
    1083             : static GEN
    1084        4718 : intnum_i(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, GEN tab, long prec)
    1085             : {
    1086        4718 :   GEN S = gen_0, res1, res2, kma, kmb;
    1087        4718 :   long sb, sgns = 1, codea = transcode(a, "a"), codeb = transcode(b, "b");
    1088             : 
    1089        4718 :   if (codea == f_REG && typ(a) == t_VEC) a = gel(a,1);
    1090        4718 :   if (codeb == f_REG && typ(b) == t_VEC) b = gel(b,1);
    1091        4718 :   if (codea == f_REG && codeb == f_REG) return intn(E, eval, a, b, tab);
    1092         994 :   if (labs(codea) > labs(codeb)) { swap(a,b); lswap(codea,codeb); sgns = -1; }
    1093             :   /* now labs(codea) <= labs(codeb) */
    1094         994 :   if (codeb == f_SING)
    1095             :   {
    1096          21 :     if (codea == f_REG)
    1097           7 :       S = intnsing(E, eval, b, a, tab, prec), sgns = -sgns;
    1098             :     else
    1099             :     {
    1100          14 :       GEN c = gmul2n(gadd(gel(a,1), gel(b,1)), -1);
    1101          14 :       res1 = intnsing(E, eval, a, c, gel(tab,1), prec);
    1102          14 :       res2 = intnsing(E, eval, b, c, gel(tab,2), prec);
    1103          14 :       S = gsub(res1, res2);
    1104             :     }
    1105          21 :     return (sgns < 0) ? gneg(S) : S;
    1106             :   }
    1107             :   /* now b is infinite */
    1108         973 :   sb = codeb > 0 ? 1 : -1;
    1109         973 :   codeb = labs(codeb);
    1110         973 :   if (codea == f_REG && codeb != f_YOSCC
    1111         266 :       && (codeb != f_YOSCS || gequal0(a)))
    1112             :   {
    1113         266 :     S = intninfpm(E, eval, a, sb*codeb, tab);
    1114         266 :     return sgns*sb < 0 ? gneg(S) : S;
    1115             :   }
    1116         707 :   if (is_fin_f(codea))
    1117             :   { /* either codea == f_SING  or codea == f_REG and codeb = f_YOSCC
    1118             :      * or (codeb == f_YOSCS and !gequal0(a)) */
    1119             :     GEN c;
    1120          14 :     GEN pi2p = gmul(Pi2n(1,prec), f_getycplx(b, prec));
    1121          14 :     GEN pis2p = gmul2n(pi2p, -2);
    1122          14 :     c = real_i(codea == f_SING ? gel(a,1) : a);
    1123          14 :     switch(codeb)
    1124             :     {
    1125             :       case f_YOSCC: case f_YOSCS:
    1126           7 :         if (codeb == f_YOSCC) c = gadd(c, pis2p);
    1127           7 :         c = gdiv(c, pi2p);
    1128           7 :         if (sb > 0)
    1129           7 :           c = addsi(1, gceil(c));
    1130             :         else
    1131           0 :           c = subis(gfloor(c), 1);
    1132           7 :         c = gmul(pi2p, c);
    1133           7 :         if (codeb == f_YOSCC) c = gsub(c, pis2p);
    1134           7 :         break;
    1135           7 :       default: c = addsi(1, gceil(c));
    1136           7 :         break;
    1137             :     }
    1138          21 :     res1 = codea==f_SING? intnsing(E, eval, a, c, gel(tab,1), prec)
    1139          21 :                         : intn    (E, eval, a, c, gel(tab,1));
    1140          14 :     res2 = intninfpm(E, eval, c, sb*codeb,gel(tab,2));
    1141          14 :     if (sb < 0) res2 = gneg(res2);
    1142          14 :     res1 = gadd(res1, res2);
    1143          14 :     return sgns < 0 ? gneg(res1) : res1;
    1144             :   }
    1145             :   /* now a and b are infinite */
    1146         693 :   if (codea * sb > 0)
    1147             :   {
    1148          14 :     if (codea > 0) pari_warn(warner, "integral from oo to oo");
    1149          14 :     if (codea < 0) pari_warn(warner, "integral from -oo to -oo");
    1150          14 :     return gen_0;
    1151             :   }
    1152         679 :   if (sb < 0) sgns = -sgns;
    1153         679 :   codea = labs(codea);
    1154         679 :   kma = f_getycplx(a, prec);
    1155         679 :   kmb = f_getycplx(b, prec);
    1156         679 :   if ((codea == f_YSLOW && codeb == f_YSLOW)
    1157         672 :    || (codea == f_YFAST && codeb == f_YFAST && gequal(kma, kmb)))
    1158         581 :     S = intninfinf(E, eval, tab);
    1159             :   else
    1160             :   {
    1161          98 :     GEN pis2 = Pi2n(-1, prec);
    1162          98 :     GEN ca = (codea == f_YOSCC)? gmul(pis2, kma): gen_0;
    1163          98 :     GEN cb = (codeb == f_YOSCC)? gmul(pis2, kmb): gen_0;
    1164          98 :     GEN c = codea == f_YOSCC ? ca : cb;
    1165          98 :     GEN SP, SN = intninfpm(E, eval, c, -sb*codea, gel(tab,1)); /*signe(a)=-sb*/
    1166          98 :     if (codea != f_YOSCC)
    1167          84 :       SP = intninfpm(E, eval, cb, sb*codeb, gel(tab,2));
    1168             :     /* codea = codeb = f_YOSCC */
    1169          14 :     else if (gequal(kma, kmb))
    1170           0 :       SP = intninfpm(E, eval, cb, sb*codeb, gel(tab,2));
    1171             :     else
    1172             :     {
    1173          14 :       tab = gel(tab,2);
    1174          14 :       SP = intninfpm(E, eval, cb, sb*codeb, gel(tab,2));
    1175          14 :       SP = gadd(SP, intn(E, eval, ca, cb, gel(tab,1)));
    1176             :     }
    1177          98 :     S = gadd(SN, SP);
    1178             :   }
    1179         679 :   if (sgns < 0) S = gneg(S);
    1180         679 :   return S;
    1181             : }
    1182             : 
    1183             : GEN
    1184        4718 : intnum(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, GEN tab, long prec)
    1185             : {
    1186        4718 :   pari_sp ltop = avma;
    1187        4718 :   long l = prec+EXTRAPREC;
    1188             :   GEN S;
    1189             : 
    1190        4718 :   tab = intnuminit0(a, b, tab, prec);
    1191        4718 :   S = intnum_i(E, eval, gprec_w(a, l), gprec_w(b, l), tab, prec);
    1192        4718 :   return gerepilecopy(ltop, gprec_wtrunc(S, prec));
    1193             : }
    1194             : 
    1195             : typedef struct auxint_s {
    1196             :   GEN a, R, pi;
    1197             :   GEN (*f)(void*, GEN);
    1198             :   GEN (*w)(GEN, long);
    1199             :   long prec;
    1200             :   void *E;
    1201             : } auxint_t;
    1202             : 
    1203             : static GEN
    1204        3675 : auxcirc(void *E, GEN t)
    1205             : {
    1206        3675 :   auxint_t *D = (auxint_t*) E;
    1207             :   GEN s, c, z;
    1208        3675 :   mpsincos(mulrr(t, D->pi), &s, &c); z = mkcomplex(c,s);
    1209        3675 :   return gmul(z, D->f(D->E, gadd(D->a, gmul(D->R, z))));
    1210             : }
    1211             : 
    1212             : GEN
    1213           7 : intcirc(void *E, GEN (*eval)(void*, GEN), GEN a, GEN R, GEN tab, long prec)
    1214             : {
    1215             :   auxint_t D;
    1216             :   GEN z;
    1217             : 
    1218           7 :   D.a = a;
    1219           7 :   D.R = R;
    1220           7 :   D.pi = mppi(prec);
    1221           7 :   D.f = eval;
    1222           7 :   D.E = E;
    1223           7 :   z = intnum(&D, &auxcirc, real_m1(prec), real_1(prec), tab, prec);
    1224           7 :   return gmul2n(gmul(R, z), -1);
    1225             : }
    1226             : 
    1227             : GEN
    1228        4536 : intnum0(GEN a, GEN b, GEN code, GEN tab, long prec)
    1229        4536 : { EXPR_WRAP(code, intnum(EXPR_ARG, a, b, tab, prec)); }
    1230             : GEN
    1231           7 : intcirc0(GEN a, GEN R, GEN code, GEN tab, long prec)
    1232           7 : { EXPR_WRAP(code, intcirc(EXPR_ARG, a, R, tab, prec)); }
    1233             : GEN
    1234         126 : intfuncinit0(GEN a, GEN b, GEN code, long m, long prec)
    1235         126 : { EXPR_WRAP(code, intfuncinit(EXPR_ARG, a, b, m, prec)); }
    1236             : 
    1237             : #if 0
    1238             : /* Two variable integration */
    1239             : 
    1240             : typedef struct auxf_s {
    1241             :   GEN x;
    1242             :   GEN (*f)(void *, GEN, GEN);
    1243             :   void *E;
    1244             : } auxf_t;
    1245             : 
    1246             : typedef struct indi_s {
    1247             :   GEN (*c)(void*, GEN);
    1248             :   GEN (*d)(void*, GEN);
    1249             :   GEN (*f)(void *, GEN, GEN);
    1250             :   void *Ec;
    1251             :   void *Ed;
    1252             :   void *Ef;
    1253             :   GEN tabintern;
    1254             :   long prec;
    1255             : } indi_t;
    1256             : 
    1257             : static GEN
    1258             : auxf(GEN y, void *E)
    1259             : {
    1260             :   auxf_t *D = (auxf_t*) E;
    1261             :   return D->f(D->E, D->x, y);
    1262             : }
    1263             : 
    1264             : static GEN
    1265             : intnumdoubintern(GEN x, void *E)
    1266             : {
    1267             :   indi_t *D = (indi_t*) E;
    1268             :   GEN c = D->c(x, D->Ec), d = D->d(x, D->Ed);
    1269             :   auxf_t A;
    1270             : 
    1271             :   A.x = x;
    1272             :   A.f = D->f;
    1273             :   A.E = D->Ef;
    1274             :   return intnum(&A, &auxf, c, d, D->tabintern, D->prec);
    1275             : }
    1276             : 
    1277             : GEN
    1278             : intnumdoub(void *Ef, GEN (*evalf)(void *, GEN, GEN), void *Ec, GEN (*evalc)(void*, GEN), void *Ed, GEN (*evald)(void*, GEN), GEN a, GEN b, GEN tabext, GEN tabint, long prec)
    1279             : {
    1280             :   indi_t E;
    1281             : 
    1282             :   E.c = evalc;
    1283             :   E.d = evald;
    1284             :   E.f = evalf;
    1285             :   E.Ec = Ec;
    1286             :   E.Ed = Ed;
    1287             :   E.Ef = Ef;
    1288             :   E.prec = prec;
    1289             :   if (typ(tabint) == t_INT)
    1290             :   {
    1291             :     GEN C = evalc(a, Ec), D = evald(a, Ed);
    1292             :     if (typ(C) != t_VEC && typ(D) != t_VEC) { C = gen_0; D = gen_1; }
    1293             :     E.tabintern = intnuminit0(C, D, tabint, prec);
    1294             :   }
    1295             :   else E.tabintern = tabint;
    1296             :   return intnum(&E, &intnumdoubintern, a, b, tabext, prec);
    1297             : }
    1298             : 
    1299             : GEN
    1300             : intnumdoub0(GEN a, GEN b, int nc, int nd, int nf, GEN tabext, GEN tabint, long prec)
    1301             : {
    1302             :   GEN z;
    1303             :   push_lex(NULL);
    1304             :   push_lex(NULL);
    1305             :   z = intnumdoub(chf, &gp_eval2, chc, &gp_eval, chd, &gp_eval, a, b, tabext, tabint, prec);
    1306             :   pop_lex(1); pop_lex(1); return z;
    1307             : }
    1308             : #endif
    1309             : 
    1310             : 
    1311             : /* The quotient-difference algorithm. Given a vector M, convert the series
    1312             :  * S = \sum_{n >= 0} M[n+1]z^n into a continued fraction.
    1313             :  * Compute the c[n] such that
    1314             :  * S = c[1] / (1 + c[2]z / (1+c[3]z/(1+...c[lim]z))),
    1315             :  * Compute A[n] and B[n] such that
    1316             :  * S = M[1]/ (1+A[1]*z+B[1]*z^2 / (1+A[2]*z+B[2]*z^2/ (1+...1/(1+A[lim\2]*z)))),
    1317             :  * Assume lim <= #M.
    1318             :  * Does not work for certain M. */
    1319             : 
    1320             : /* Given a continued fraction CF output by the quodif program,
    1321             : convert it into an Euler continued fraction A(n), B(n), where
    1322             : $1/(1+c[2]z/(1+c[3]z/(1+..c[lim]z)))
    1323             : =1/(1+A[1]*z+B[1]*z^2/(1+A[2]*z+B[2]*z^2/(1+...1/(1+A[lim\2]*z)))). */
    1324             : static GEN
    1325        2359 : contfrac_Euler(GEN CF)
    1326             : {
    1327        2359 :   long lima, limb, i, lim = lg(CF)-1;
    1328             :   GEN A, B;
    1329        2359 :   lima = lim/2;
    1330        2359 :   limb = (lim - 1)/2;
    1331        2359 :   A = cgetg(lima+1, t_VEC);
    1332        2359 :   B = cgetg(limb+1, t_VEC);
    1333        2359 :   gel (A, 1) = gel(CF, 2);
    1334        2359 :   for (i=2; i <= lima; ++i) gel(A,i) = gadd(gel(CF, 2*i), gel(CF, 2*i-1));
    1335        2359 :   for (i=1; i <= limb; ++i) gel(B,i) = gneg(gmul(gel(CF, 2*i+1), gel(CF, 2*i)));
    1336        2359 :   return mkvec2(A, B);
    1337             : }
    1338             : 
    1339             : static GEN
    1340        2576 : contfracinit_i(GEN M, long lim)
    1341             : {
    1342             :   pari_sp av;
    1343             :   GEN e, q, c;
    1344             :   long lim2;
    1345             :   long j, k;
    1346        2576 :   e = zerovec(lim);
    1347        2576 :   c = zerovec(lim+1); gel(c, 1) = gel(M, 1);
    1348        2576 :   q = cgetg(lim+1, t_VEC);
    1349        2576 :   for (k = 1; k <= lim; ++k) gel(q, k) = gdiv(gel(M, k+1), gel(M, k));
    1350        2576 :   lim2 = lim/2; av = avma;
    1351      103452 :   for (j = 1; j <= lim2; ++j)
    1352             :   {
    1353      100876 :     long l = lim - 2*j;
    1354      100876 :     gel(c, 2*j) = gneg(gel(q, 1));
    1355     9781301 :     for (k = 0; k <= l; ++k)
    1356     9680425 :       gel(e, k+1) = gsub(gadd(gel(e, k+2), gel(q, k+2)), gel(q, k+1));
    1357     9680425 :     for (k = 0; k < l; ++k)
    1358     9579549 :       gel(q, k+1) = gdiv(gmul(gel(q, k+2), gel(e, k+2)), gel(e, k+1));
    1359      100876 :     gel(c, 2*j+1) = gneg(gel(e, 1));
    1360      100876 :     if (gc_needed(av, 3))
    1361             :     {
    1362          98 :       if (DEBUGMEM>1) pari_warn(warnmem,"contfracinit, %ld/%ld",j,lim2);
    1363          98 :       gerepileall(av, 3, &e, &c, &q);
    1364             :     }
    1365             :   }
    1366        2576 :   if (odd(lim)) gel(c, lim+1) = gneg(gel(q, 1));
    1367        2576 :   return c;
    1368             : }
    1369             : 
    1370             : GEN
    1371        2380 : contfracinit(GEN M, long lim)
    1372             : {
    1373        2380 :   pari_sp ltop = avma;
    1374             :   GEN c;
    1375        2380 :   switch(typ(M))
    1376             :   {
    1377             :     case t_RFRAC:
    1378           0 :       if (lim < 0) pari_err_TYPE("contfracinit",M);
    1379           0 :       M = gadd(M, zeroser(gvar(M), lim + 2)); /*fall through*/
    1380          28 :     case t_SER: M = gtovec(M); break;
    1381           7 :     case t_POL: M = gtovecrev(M); break;
    1382        2338 :     case t_VEC: case t_COL: break;
    1383           7 :     default: pari_err_TYPE("contfracinit", M);
    1384             :   }
    1385        2373 :   if (lim < 0)
    1386          28 :     lim = lg(M)-2;
    1387        2345 :   else if (lg(M)-1 <= lim)
    1388           0 :     pari_err_COMPONENT("contfracinit", "<", stoi(lg(M)-1), stoi(lim));
    1389        2373 :   if (lim < 0) retmkvec2(cgetg(1,t_VEC),cgetg(1,t_VEC));
    1390        2359 :   c = contfracinit_i(M, lim);
    1391        2359 :   return gerepilecopy(ltop, contfrac_Euler(c));
    1392             : }
    1393             : 
    1394             : /* Evaluate at 1/tinv the nlim first terms of the continued fraction output by
    1395             :  * contfracinit. */
    1396             : /* Not stack clean */
    1397             : GEN
    1398     1873029 : contfraceval_inv(GEN CF, GEN tinv, long nlim)
    1399             : {
    1400             :   pari_sp btop;
    1401             :   long j;
    1402     1873029 :   GEN S = gen_0, S1, S2, A, B;
    1403     1873029 :   if (typ(CF) != t_VEC || lg(CF) != 3) pari_err_TYPE("contfraceval", CF);
    1404     1873029 :   A = gel(CF, 1); if (typ(A) != t_VEC) pari_err_TYPE("contfraceval", CF);
    1405     1873029 :   B = gel(CF, 2); if (typ(B) != t_VEC) pari_err_TYPE("contfraceval", CF);
    1406     1873029 :   if (nlim < 0)
    1407          14 :     nlim = lg(A)-1;
    1408     1873015 :   else if (lg(A) <= nlim)
    1409           7 :     pari_err_COMPONENT("contfraceval", ">", stoi(lg(A)-1), stoi(nlim));
    1410     1873022 :   if (lg(B)+1 <= nlim)
    1411           0 :     pari_err_COMPONENT("contfraceval", ">", stoi(lg(B)), stoi(nlim));
    1412     1873022 :   btop = avma;
    1413     1873022 :   if (nlim <= 1) return lg(A)==1? gen_0: gdiv(tinv, gadd(gel(A, 1), tinv));
    1414     1861388 :   switch(nlim % 3)
    1415             :   {
    1416             :     case 2:
    1417      638880 :       S = gdiv(gel(B, nlim-1), gadd(gel(A, nlim), tinv));
    1418      638880 :       nlim--; break;
    1419             : 
    1420             :     case 0:
    1421      661709 :       S1 = gadd(gel(A, nlim), tinv);
    1422      661709 :       S2 = gadd(gmul(gadd(gel(A, nlim-1), tinv), S1), gel(B, nlim-1));
    1423      661709 :       S = gdiv(gmul(gel(B, nlim-2), S1), S2);
    1424      661709 :       nlim -= 2; break;
    1425             :   }
    1426             :   /* nlim = 1 (mod 3) */
    1427     8618111 :   for (j = nlim; j >= 4; j -= 3)
    1428             :   {
    1429             :     GEN S3;
    1430     6756723 :     S1 = gadd(gadd(gel(A, j), tinv), S);
    1431     6756723 :     S2 = gadd(gmul(gadd(gel(A, j-1), tinv), S1), gel(B, j-1));
    1432     6756723 :     S3 = gadd(gmul(gadd(gel(A, j-2), tinv), S2), gmul(gel(B, j-2), S1));
    1433     6756723 :     S = gdiv(gmul(gel(B, j-3), S2), S3);
    1434     6756723 :     if (gc_needed(btop, 3)) S = gerepilecopy(btop, S);
    1435             :   }
    1436     1861388 :   return gdiv(tinv, gadd(gadd(gel(A, 1), tinv), S));
    1437             : }
    1438             : 
    1439             : GEN
    1440          35 : contfraceval(GEN CF, GEN t, long nlim)
    1441             : {
    1442          35 :   pari_sp ltop = avma;
    1443          35 :   return gerepileupto(ltop, contfraceval_inv(CF, ginv(t), nlim));
    1444             : }
    1445             : 
    1446             : /* MONIEN SUMMATION */
    1447             : 
    1448             : /* basic Newton, find x ~ z such that Q(x) = 0 */
    1449             : static GEN
    1450        1561 : monrefine(GEN Q, GEN QP, GEN z, long prec)
    1451             : {
    1452        1561 :   pari_sp av = avma;
    1453        1561 :   GEN pr = poleval(Q, z);
    1454             :   for(;;)
    1455             :   {
    1456             :     GEN prnew;
    1457        6868 :     z = gsub(z, gdiv(pr, poleval(QP, z)));
    1458        6868 :     prnew = poleval(Q, z);
    1459        6868 :     if (gcmp(gabs(prnew, prec), gabs(pr, prec)) >= 0) break;
    1460        5307 :     pr = prnew;
    1461        5307 :   }
    1462        1561 :   z = gprec_w(z, 2*prec-2);
    1463        1561 :   z = gsub(z, gdiv(poleval(Q, z), poleval(QP, z)));
    1464        1561 :   return gerepileupto(av, z);
    1465             : }
    1466             : 
    1467             : static GEN
    1468         217 : RX_realroots(GEN x, long prec)
    1469         217 : { return realroots(gprec_wtrunc(x,prec), NULL, prec); }
    1470             : 
    1471             : /* (real) roots of Q, assuming QP = Q' and that half the roots are close to
    1472             :  * k+1, ..., k+m, m = deg(Q)/2-1. N.B. All roots are real and >= 1 */
    1473             : static GEN
    1474         119 : monroots(GEN Q, GEN QP, long k, long prec)
    1475             : {
    1476         119 :   long j, n = degpol(Q), m = n/2 - 1;
    1477         119 :   GEN v2, v1 = cgetg(m+1, t_VEC);
    1478         119 :   for (j = 1; j <= m; ++j) gel(v1, j) = monrefine(Q, QP, stoi(k+j), prec);
    1479         119 :   Q = gdivent(Q, roots_to_pol(v1, varn(Q)));
    1480         119 :   v2 = RX_realroots(Q, prec); settyp(v2, t_VEC);
    1481         119 :   return shallowconcat(v1, v2);
    1482             : }
    1483             : 
    1484             : static void
    1485         217 : Pade(GEN M, GEN *pP, GEN *pQ)
    1486             : {
    1487         217 :   pari_sp av = avma;
    1488         217 :   long n = lg(M)-2, i;
    1489         217 :   GEN v = contfracinit_i(M, n), P = pol_0(0), Q = pol_1(0);
    1490             :   /* evaluate continued fraction => Pade approximants */
    1491       12803 :   for (i = n-1; i >= 1; i--)
    1492             :   { /* S = P/Q: S -> v[i]*x / (1+S) */
    1493       12586 :     GEN R = RgX_shift_shallow(RgX_Rg_mul(Q,gel(v,i)), 1);
    1494       12586 :     Q = RgX_add(P,Q); P = R;
    1495       12586 :     if (gc_needed(av, 3))
    1496             :     {
    1497           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Pade, %ld/%ld",i,n-1);
    1498           0 :       gerepileall(av, 3, &P, &Q, &v);
    1499             :     }
    1500             :   }
    1501             :   /* S -> 1+S */
    1502         217 :   *pP = RgX_add(P,Q);
    1503         217 :   *pQ = Q;
    1504         217 : }
    1505             : 
    1506             : static GEN
    1507           7 : veczetaprime(GEN a, GEN b, long N, long prec)
    1508             : {
    1509           7 :   long newprec, fpr = prec2nbits(prec), pr = (long)ceil(fpr * 1.5);
    1510           7 :   long l = nbits2prec(pr), e = fpr / 2;
    1511             :   GEN eps, A, B;
    1512           7 :   newprec = nbits2prec(pr + BITS_IN_LONG);
    1513           7 :   a = gprec_w(a, newprec);
    1514           7 :   b = gprec_w(b, newprec);
    1515           7 :   eps = real2n(-e, l);
    1516           7 :   A = veczeta(a, gsub(b, eps), N, newprec);
    1517           7 :   B = veczeta(a, gadd(b, eps), N, newprec);
    1518           7 :   return gmul2n(RgV_sub(B, A), e-1);
    1519             : }
    1520             : 
    1521             : struct mon_w {
    1522             :   GEN w, a, b;
    1523             :   long n, j, prec;
    1524             : };
    1525             : 
    1526             : /* w(x) / x^(a*(j+k)+b), k >= 1; w a t_CLOSURE or t_INT [encodes log(x)^w] */
    1527             : static GEN
    1528       36508 : wrapmonw(void* E, GEN x)
    1529             : {
    1530       36508 :   struct mon_w *W = (struct mon_w*)E;
    1531       36508 :   long k, j = W->j, n = W->n, prec = W->prec, l = 2*n+4-j;
    1532      109524 :   GEN wx = typ(W->w) == t_CLOSURE? closure_callgen1prec(W->w, x, prec)
    1533       73016 :                                  : powgi(glog(x, prec), W->w);
    1534       36508 :   GEN v = cgetg(l, t_VEC);
    1535       36508 :   GEN xa = gpow(x, gneg(W->a), prec), w = gmul(wx, gpowgs(xa, j));
    1536       36508 :   w = gdiv(w, gpow(x,W->b,prec));
    1537       36508 :   for (k = 1; k < l; k++) { gel(v,k) = w; w = gmul(w, xa); }
    1538       36508 :   return v;
    1539             : }
    1540             : /* w(x) / x^(a*j+b) */
    1541             : static GEN
    1542       21387 : wrapmonw2(void* E, GEN x)
    1543             : {
    1544       21387 :   struct mon_w *W = (struct mon_w*)E;
    1545       21387 :   GEN wnx = closure_callgen1prec(W->w, x, W->prec);
    1546       21387 :   return gdiv(wnx, gpow(x, gadd(gmulgs(W->a, W->j), W->b), W->prec));
    1547             : }
    1548             : static GEN
    1549          35 : M_from_wrapmon(struct mon_w *S, GEN wfast, GEN n0)
    1550             : {
    1551          35 :   long j, N = 2*S->n+2;
    1552          35 :   GEN M = cgetg(N+1, t_VEC), faj = gsub(wfast, S->b);
    1553          42 :   for (j = 1; j <= N; j++)
    1554             :   {
    1555          42 :     faj = gsub(faj, S->a);
    1556          42 :     if (gcmpgs(faj, -2) <= 0)
    1557             :     {
    1558          28 :       S->j = j; setlg(M,j);
    1559          28 :       M = shallowconcat(M, sumnum((void*)S, wrapmonw, n0, NULL, S->prec));
    1560          28 :       break;
    1561             :     }
    1562          14 :     S->j = j;
    1563          14 :     gel(M,j) = sumnum((void*)S, wrapmonw2, mkvec2(n0,faj), NULL, S->prec);
    1564             :   }
    1565          28 :   return M;
    1566             : }
    1567             : 
    1568             : static void
    1569         168 : checkmonroots(GEN vr, long n)
    1570             : {
    1571         168 :   if (lg(vr) != n+1)
    1572           0 :     pari_err_IMPL("recovery when missing roots in sumnummonieninit");
    1573         168 : }
    1574             : 
    1575             : static GEN
    1576         175 : sumnummonieninit_i(GEN a, GEN b, GEN w, GEN n0, long prec)
    1577             : {
    1578         175 :   GEN c, M, P, Q, Qp, vr, vabs, vwt, ga = gadd(a, b);
    1579         175 :   double bit = 2*prec2nbits(prec) / gtodouble(ga), D = bit*LOG2;
    1580         175 :   double da = maxdd(1., gtodouble(a));
    1581         175 :   long n = (long)ceil(D / (da*(log(D)-1)));
    1582         175 :   long j, prec2, prec0 = prec + EXTRAPRECWORD;
    1583         175 :   double bit0 = ceil((2*n+1)*LOG2_10);
    1584         175 :   int neg = 1;
    1585             :   struct mon_w S;
    1586             : 
    1587         175 :   prec = nbits2prec(maxdd(2*da*bit, bit0));
    1588         175 :   prec2 = nbits2prec(maxdd(1.3*da*bit, bit0));
    1589         175 :   S.w = w;
    1590         175 :   S.a = a = gprec_w(a, 2*prec-2);
    1591         175 :   S.b = b = gprec_w(b, 2*prec-2);
    1592         175 :   S.n = n;
    1593         175 :   S.j = 1;
    1594         175 :   S.prec = prec;
    1595         175 :   if (typ(w) == t_INT)
    1596             :   { /* f(n) ~ \sum_{i > 0} f_i log(n)^k / n^(a*i + b); a > 0, a+b > 1 */
    1597         140 :     long k = itos(w);
    1598         140 :     if (k == 0)
    1599         133 :       M = veczeta(a, ga, 2*n+2, prec);
    1600           7 :     else if (k == 1)
    1601           7 :       M = veczetaprime(a, ga, 2*n+2, prec);
    1602             :     else
    1603           0 :       M = M_from_wrapmon(&S, gen_0, gen_1);
    1604         140 :     if (odd(k)) neg = 0;
    1605             :   }
    1606             :   else
    1607             :   {
    1608          35 :     GEN wfast = gen_0;
    1609          35 :     if (typ(w) == t_VEC) { wfast = gel(w,2); w = gel(w,1); }
    1610          35 :     M = M_from_wrapmon(&S, wfast, n0);
    1611             :   }
    1612             :   /* M[j] = sum(n >= n0, w(n) / n^(a*j+b) */
    1613         168 :   Pade(M, &P,&Q);
    1614         168 :   Qp = RgX_deriv(Q);
    1615         168 :   if (gequal1(a)) a = NULL;
    1616         168 :   if (!a && typ(w) == t_INT)
    1617             :   {
    1618         119 :     vabs = vr = monroots(Q, Qp, signe(w)? 1: 0, prec2);
    1619         119 :     checkmonroots(vr, n);
    1620         119 :     c = b;
    1621             :   }
    1622             :   else
    1623             :   {
    1624          49 :     vr = RX_realroots(Q, prec2); settyp(vr, t_VEC);
    1625          49 :     checkmonroots(vr, n);
    1626          49 :     if (!a) { vabs = vr; c = b; }
    1627             :     else
    1628             :     {
    1629          35 :       GEN ai = ginv(a);
    1630          35 :       vabs = cgetg(n+1, t_VEC);
    1631          35 :       for (j = 1; j <= n; j++) gel(vabs,j) = gpow(gel(vr,j), ai, prec2);
    1632          35 :       c = gdiv(b,a);
    1633             :     }
    1634             :   }
    1635             : 
    1636         168 :   vwt = cgetg(n+1, t_VEC);
    1637         168 :   c = gsubgs(c,1); if (gequal0(c)) c = NULL;
    1638        5012 :   for (j = 1; j <= n; j++)
    1639             :   {
    1640        4844 :     GEN r = gel(vr,j), t = gdiv(poleval(P,r), poleval(Qp,r));
    1641        4844 :     if (c) t = gmul(t, gpow(r, c, prec));
    1642        4844 :     gel(vwt,j) = neg? gneg(t): t;
    1643             :   }
    1644         168 :   if (typ(w) == t_INT && !equali1(n0))
    1645             :   {
    1646          28 :     GEN h = subiu(n0,1);
    1647          28 :     for (j = 1; j <= n; j++) gel(vabs,j) = gadd(gel(vabs,j), h);
    1648             :   }
    1649         168 :   return mkvec3(gprec_wtrunc(vabs,prec0), gprec_wtrunc(vwt,prec0), n0);
    1650             : }
    1651             : 
    1652             : GEN
    1653         105 : sumnummonieninit(GEN asymp, GEN w, GEN n0, long prec)
    1654             : {
    1655         105 :   pari_sp av = avma;
    1656         105 :   const char *fun = "sumnummonieninit";
    1657             :   GEN a, b;
    1658         105 :   if (!n0) n0 = gen_1; else if (typ(n0) != t_INT) pari_err_TYPE(fun, n0);
    1659         105 :   if (!asymp) a = b = gen_1;
    1660             :   else
    1661             :   {
    1662          77 :     if (typ(asymp) == t_VEC)
    1663             :     {
    1664          70 :       if (lg(asymp) != 3) pari_err_TYPE(fun, asymp);
    1665          70 :       a = gel(asymp,1);
    1666          70 :       b = gel(asymp,2);
    1667             :     }
    1668             :     else
    1669             :     {
    1670           7 :       a = gen_1;
    1671           7 :       b = asymp;
    1672             :     }
    1673          77 :     if (gsigne(a) <= 0) pari_err_DOMAIN(fun, "a", "<=", gen_0, a);
    1674          70 :     if (gcmpgs(gadd(a,b), 1) <= 0)
    1675           7 :       pari_err_DOMAIN(fun, "a+b", "<=", gen_1, mkvec2(a,b));
    1676             :   }
    1677          91 :   if (!w) w = gen_0;
    1678          42 :   else switch(typ(w))
    1679             :   {
    1680             :     case t_INT:
    1681           7 :       if (signe(w) < 0) pari_err_IMPL("log power < 0 in sumnummonieninit");
    1682          35 :     case t_CLOSURE: break;
    1683             :     case t_VEC:
    1684           7 :       if (lg(w) == 3 && typ(gel(w,1)) == t_CLOSURE) break;
    1685           0 :     default: pari_err_TYPE(fun, w);
    1686             :   }
    1687          91 :   return gerepilecopy(av, sumnummonieninit_i(a, b, w, n0, prec));
    1688             : }
    1689             : 
    1690             : GEN
    1691         175 : sumnummonien(void *E, GEN (*eval)(void*,GEN), GEN n0, GEN tab, long prec)
    1692             : {
    1693         175 :   pari_sp av = avma;
    1694             :   GEN vabs, vwt, S;
    1695             :   long l, i;
    1696         175 :   if (typ(n0) != t_INT) pari_err_TYPE("sumnummonien", n0);
    1697         175 :   if (!tab)
    1698          84 :     tab = sumnummonieninit_i(gen_1, gen_1, gen_0, n0, prec);
    1699             :   else
    1700             :   {
    1701          91 :     if (lg(tab) != 4 || typ(tab) != t_VEC) pari_err_TYPE("sumnummonien", tab);
    1702          91 :     if (!equalii(n0, gel(tab,3)))
    1703           7 :       pari_err(e_MISC, "incompatible initial value %Ps != %Ps", gel(tab,3),n0);
    1704             :   }
    1705         168 :   vabs= gel(tab,1); l = lg(vabs);
    1706         168 :   vwt = gel(tab,2);
    1707         168 :   if (typ(vabs) != t_VEC || typ(vwt) != t_VEC || lg(vwt) != l)
    1708           0 :     pari_err_TYPE("sumnummonien", tab);
    1709         168 :   S = gen_0;
    1710         168 :   for (i = 1; i < l; i++) S = gadd(S, gmul(gel(vwt,i), eval(E, gel(vabs,i))));
    1711         168 :   return gerepileupto(av, gprec_w(S, prec));
    1712             : }
    1713             : 
    1714             : static GEN
    1715         196 : get_oo(GEN fast) { return mkvec2(mkoo(), fast); }
    1716             : 
    1717             : GEN
    1718         119 : sumnuminit(GEN fast, long prec)
    1719             : {
    1720             :   pari_sp av;
    1721         119 :   GEN s, v, d, C, D, res = cgetg(6, t_VEC);
    1722         119 :   long bitprec = prec2nbits(prec), N, k, k2, m;
    1723             :   double w;
    1724             : 
    1725         119 :   d = mkfrac(gen_1, utoipos(4)); /* 1/4 */
    1726         119 :   gel(res, 1) = d;
    1727         119 :   av = avma;
    1728         119 :   w = gtodouble(glambertW(ginv(d), LOWDEFAULTPREC));
    1729         119 :   N = (long)ceil(LOG2*bitprec/(w*(1+w))+5);
    1730         119 :   k = (long)ceil(N*w); if (k&1) k--;
    1731             : 
    1732         119 :   prec += EXTRAPRECWORD;
    1733         119 :   s = RgX_to_ser(monomial(d,1,0), k+3);
    1734         119 :   s = gdiv(gasinh(s, prec), d); /* asinh(dx)/d */
    1735         119 :   s = gsub(ginv(gsubgs(gexp(s,prec), 1)), ginv(s));
    1736         119 :   k2 = k/2;
    1737         119 :   C = matpascal(k-1);
    1738         119 :   D = gpowers(ginv(gmul2n(d,1)), k-1);
    1739         119 :   v = cgetg(k2+1, t_VEC);
    1740        8331 :   for (m = 1; m <= k2; m++)
    1741             :   {
    1742        8212 :     pari_sp av = avma;
    1743        8212 :     GEN S = real_0(prec);
    1744             :     long j;
    1745      478289 :     for (j = m; j <= k2; j++)
    1746             :     { /* s[X^(2j-1)] * binomial(2*j-1, j-m) / (2d)^(2j-1) */
    1747      470077 :       GEN t = gmul(gmul(gel(s,2*j+1), gcoeff(C, 2*j,j-m+1)), gel(D, 2*j));
    1748      470077 :       S = odd(j)? gsub(S,t): gadd(S,t);
    1749             :     }
    1750        8212 :     if (odd(m)) S = gneg(S);
    1751        8212 :     gel(v,m) = gerepileupto(av, S);
    1752             :   }
    1753         119 :   v = RgC_gtofp(v,prec); settyp(v, t_VEC);
    1754         119 :   gel(res, 4) = gerepileupto(av, v);
    1755         119 :   gel(res, 2) = utoi(N);
    1756         119 :   gel(res, 3) = utoi(k);
    1757         119 :   if (!fast) fast = get_oo(gen_0);
    1758         119 :   gel(res, 5) = intnuminit(gel(res,2), fast, 0, prec);
    1759         119 :   return res;
    1760             : }
    1761             : 
    1762             : static int
    1763          28 : checksumtab(GEN T)
    1764             : {
    1765          28 :   if (typ(T) != t_VEC || lg(T) != 6) return 0;
    1766          21 :   return typ(gel(T,2))==t_INT && typ(gel(T,3))==t_INT && typ(gel(T,4))==t_VEC;
    1767             : }
    1768             : GEN
    1769         133 : sumnum(void *E, GEN (*eval)(void*, GEN), GEN a, GEN tab, long prec)
    1770             : {
    1771         133 :   pari_sp av = avma, av2;
    1772             :   GEN v, tabint, S, d, fast;
    1773             :   long as, N, k, m, prec2;
    1774         133 :   if (!a) { a = gen_1; fast = get_oo(gen_0); }
    1775         133 :   else switch(typ(a))
    1776             :   {
    1777             :   case t_VEC:
    1778          49 :     if (lg(a) != 3) pari_err_TYPE("sumnum", a);
    1779          49 :     fast = get_oo(gel(a,2));
    1780          49 :     a = gel(a,1); break;
    1781             :   default:
    1782          84 :     fast = get_oo(gen_0);
    1783             :   }
    1784         133 :   if (typ(a) != t_INT) pari_err_TYPE("sumnum", a);
    1785         133 :   if (!tab) tab = sumnuminit(fast, prec);
    1786          28 :   else if (!checksumtab(tab)) pari_err_TYPE("sumnum",tab);
    1787         126 :   as = itos(a);
    1788         126 :   d = gel(tab,1);
    1789         126 :   N = maxss(as, itos(gel(tab,2)));
    1790         126 :   k = itos(gel(tab,3));
    1791         126 :   v = gel(tab,4);
    1792         126 :   tabint = gel(tab,5);
    1793         126 :   prec2 = prec+EXTRAPRECWORD;
    1794         126 :   av2 = avma;
    1795         126 :   S = gmul(eval(E, stoi(N)), real2n(-1,prec2));
    1796       15557 :   for (m = as; m < N; m++)
    1797             :   {
    1798       15431 :     S = gadd(S, eval(E, stoi(m)));
    1799       15431 :     if (gc_needed(av, 2))
    1800             :     {
    1801           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"sumnum [1], %ld/%ld",m,N);
    1802           0 :       S = gerepileupto(av2, S);
    1803             :     }
    1804             :   }
    1805        9493 :   for (m = 1; m <= k/2; m++)
    1806             :   {
    1807        9367 :     GEN t = gmulsg(2*m-1, d);
    1808        9367 :     GEN s = gsub(eval(E, gsubsg(N,t)), eval(E, gaddsg(N,t)));
    1809        9367 :     S = gadd(S, gmul(gel(v,m), s));
    1810        9367 :     if (gc_needed(av2, 2))
    1811             :     {
    1812           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"sumnum [2], %ld/%ld",m,k/2);
    1813           0 :       S = gerepileupto(av2, S);
    1814             :     }
    1815             :   }
    1816         126 :   S = gadd(S, intnum(E, eval,stoi(N), fast, tabint, prec2));
    1817         126 :   return gerepilecopy(av, gprec_w(S, prec));
    1818             : }
    1819             : 
    1820             : GEN
    1821         175 : sumnummonien0(GEN a, GEN code, GEN tab, long prec)
    1822         175 : { EXPR_WRAP(code, sumnummonien(EXPR_ARG, a, tab, prec)); }
    1823             : GEN
    1824          91 : sumnum0(GEN a, GEN code, GEN tab, long prec)
    1825          91 : { EXPR_WRAP(code, sumnum(EXPR_ARG, a, tab, prec)); }
    1826             : 
    1827             : /* Abel-Plana summation */
    1828             : 
    1829             : static GEN
    1830          49 : intnumgauexpinit(long prec)
    1831             : {
    1832          49 :   pari_sp ltop = avma;
    1833             :   GEN V, N, E, P, Q, R, vabs, vwt;
    1834          49 :   long l, n, k, j, prec2, prec0 = prec + EXTRAPRECWORD, bit = prec2nbits(prec);
    1835             : 
    1836          49 :   n = (long)ceil(bit*0.226);
    1837          49 :   n |= 1; /* make n odd */
    1838          49 :   prec = nbits2prec(1.5*bit + 32);
    1839          49 :   prec2 = maxss(prec0, nbits2prec(1.15*bit + 32));
    1840          49 :   V = cgetg(n + 4, t_VEC);
    1841        3045 :   for (k = 1; k <= n + 3; ++k)
    1842        2996 :     gel(V, k) = gtofp(gdivgs(bernfrac(2*k), odd(k)? 2*k: -2*k), prec);
    1843          49 :   Pade(V, &P, &Q);
    1844          49 :   N = RgX_recip(gsub(P, Q));
    1845          49 :   E = RgX_recip(Q);
    1846          49 :   R = gdivgs(gdiv(N, RgX_deriv(E)), 2);
    1847          49 :   vabs = RX_realroots(E,prec2);
    1848          49 :   l = lg(vabs); settyp(vabs, t_VEC);
    1849          49 :   vwt = cgetg(l, t_VEC);
    1850        1498 :   for (j = 1; j < l; ++j)
    1851             :   {
    1852        1449 :     GEN a = gel(vabs,j);
    1853        1449 :     gel(vwt, j) = gprec_wtrunc(poleval(R, a), prec0);
    1854        1449 :     gel(vabs, j) = gprec_wtrunc(sqrtr_abs(a), prec0);
    1855             :   }
    1856          49 :   return gerepilecopy(ltop, mkvec2(vabs, vwt));
    1857             : }
    1858             : 
    1859             : /* Compute $\int_{-\infty}^\infty w(x)f(x)\,dx$, where $w(x)=x/(exp(2\pi x)-1)$
    1860             :  * for $x>0$ and $w(-x)=w(x)$. For Abel-Plana (sumnumap). */
    1861             : static GEN
    1862          49 : intnumgauexp(void *E, GEN (*eval)(void*,GEN), GEN gN, GEN tab, long prec)
    1863             : {
    1864          49 :   pari_sp av = avma;
    1865          49 :   GEN U = mkcomplex(gN, gen_0), V = mkcomplex(gel(U,1), gen_0), S = gen_0;
    1866          49 :   GEN vabs = gel(tab, 1), vwt = gel(tab, 2);
    1867          49 :   long l = lg(vabs), i;
    1868          49 :   if (lg(vwt) != l || typ(vabs) != t_VEC || typ(vwt) != t_VEC)
    1869           0 :     pari_err_TYPE("sumnumap", tab);
    1870        1498 :   for (i = 1; i < l; i++)
    1871             :   { /* I * (w_i/a_i) * (f(N + I*a_i) - f(N - I*a_i)) */
    1872        1449 :     GEN x = gel(vabs,i), w = gel(vwt,i), t;
    1873        1449 :     gel(U,2) = x;
    1874        1449 :     gel(V,2) = gneg(x);
    1875        1449 :     t = mulcxI(gsub(eval(E,U), eval(E,V)));
    1876        1449 :     if (typ(t) == t_COMPLEX && gequal0(gel(t,2))) t = gel(t,1);
    1877        1449 :     S = gadd(S, gmul(gdiv(w,x), t));
    1878             :   }
    1879          49 :   return gerepileupto(av, gprec_w(S, prec));
    1880             : }
    1881             : 
    1882             : GEN
    1883          49 : sumnumapinit(GEN fast, long prec)
    1884             : {
    1885          49 :   if (!fast) fast = mkoo();
    1886          49 :   retmkvec2(intnumgauexpinit(prec), intnuminit(gen_1, fast, 0, prec));
    1887             : }
    1888             : 
    1889             : typedef struct {
    1890             :   GEN (*f)(void *E, GEN);
    1891             :   void *E;
    1892             :   long N;
    1893             : } expfn;
    1894             : 
    1895             : /* f(Nx) */
    1896             : static GEN
    1897       31157 : _exfn(void *E, GEN x)
    1898             : {
    1899       31157 :   expfn *S = (expfn*)E;
    1900       31157 :   return S->f(S->E, gmulsg(S->N, x));
    1901             : }
    1902             : 
    1903             : GEN
    1904          56 : sumnumap(void *E, GEN (*eval)(void*,GEN), GEN a, GEN tab, long prec)
    1905             : {
    1906          56 :   pari_sp av = avma;
    1907             :   expfn T;
    1908             :   GEN S, fast, gN;
    1909             :   long as, m, N;
    1910          56 :   if (!a) { a = gen_1; fast = get_oo(gen_0); }
    1911          56 :   else switch(typ(a))
    1912             :   {
    1913             :     case t_VEC:
    1914          21 :       if (lg(a) != 3) pari_err_TYPE("sumnumap", a);
    1915          21 :       fast = get_oo(gel(a,2));
    1916          21 :       a = gel(a,1); break;
    1917             :     default:
    1918          35 :       fast = get_oo(gen_0);
    1919             :   }
    1920          56 :   if (typ(a) != t_INT) pari_err_TYPE("sumnumap", a);
    1921          56 :   if (!tab) tab = sumnumapinit(fast, prec);
    1922           7 :   else if (typ(tab) != t_VEC || lg(tab) != 3) pari_err_TYPE("sumnumap",tab);
    1923          49 :   as = itos(a);
    1924          49 :   T.N = N = maxss(as + 1, (long)ceil(prec2nbits(prec)*0.327));
    1925          49 :   T.E = E;
    1926          49 :   T.f = eval;
    1927          49 :   gN = stoi(N);
    1928          49 :   S = gtofp(gmul2n(eval(E, gN), -1), prec);
    1929          49 :   for (m = as; m < N; ++m) S = gadd(S, eval(E, stoi(m)));
    1930          49 :   S = gadd(S, gmulsg(N, intnum(&T, &_exfn, gen_1, fast, gel(tab, 2), prec)));
    1931          49 :   S = gadd(S, intnumgauexp(E, eval, gN, gel(tab, 1), prec));
    1932          49 :   return gerepileupto(av, S);
    1933             : }
    1934             : 
    1935             : GEN
    1936          56 : sumnumap0(GEN a, GEN code, GEN tab, long prec)
    1937          56 : { EXPR_WRAP(code, sumnumap(EXPR_ARG, a, tab, prec)); }
    1938             : 
    1939             : 
    1940             : /* max (1, |zeros|), P a t_POL or scalar */
    1941             : static GEN
    1942         189 : polmax(GEN P)
    1943             : {
    1944             :   GEN r;
    1945         189 :   if (typ(P) != t_POL || degpol(P) <= 0) return gen_1;
    1946         189 :   r = polrootsbound(P);
    1947         189 :   if (gcmpgs(r, 1) <= 0) return gen_1;
    1948         154 :   return r;
    1949             : }
    1950             : 
    1951             : /* max (1, |poles|), F a t_POL or t_RFRAC or scalar */
    1952             : static GEN
    1953         119 : ratpolemax(GEN F)
    1954             : {
    1955         119 :   if (typ(F) == t_POL) return gen_1;
    1956         119 :   return polmax(gel(F, 2));
    1957             : }
    1958             : /* max (1, |poles|, |zeros|), sets *p = max(1, |poles|)) */
    1959             : static GEN
    1960          35 : ratpolemax2(GEN F, GEN *p)
    1961             : {
    1962             :   GEN t;
    1963          35 :   if (typ(F) == t_POL) { if (p) *p = gen_1; return polmax(F); }
    1964          35 :   t = polmax(gel(F,2)); if (p) *p = t;
    1965          35 :   return gmax(t, polmax(gel(F,1)));
    1966             : }
    1967             : 
    1968             : static GEN
    1969       18074 : sercoeff(GEN x, long n)
    1970             : {
    1971       18074 :   long N = n - valp(x);
    1972       18074 :   return (N < 0)? gen_0: gel(x,N+2);
    1973             : }
    1974             : 
    1975             : /* Compute the integral from N to infinity of a rational function F, deg(F) < -1
    1976             :  * We must have N > 2 * max(1, |poles F|). */
    1977             : static GEN
    1978          56 : intnumainfrat(GEN F, long N, long prec)
    1979             : {
    1980          56 :   pari_sp av = avma;
    1981             :   GEN r, S, ser;
    1982          56 :   long B = prec2nbits(prec), v, k, lim;
    1983             : 
    1984          56 :   r = gdivsg(N, ratpolemax(F));
    1985          56 :   if (gcmpgs(r, 2) < 0) pari_err(e_MISC, "N too small in intnumainfrat");
    1986          56 :   lim = (long)ceil(B/dbllog2(r)) + 1;
    1987          56 :   ser = gmul(F, real_1(prec + EXTRAPREC));
    1988          56 :   ser = rfracrecip_to_ser_absolute(ser, lim + 2);
    1989          56 :   v = valp(ser);
    1990          56 :   S = gen_0;
    1991             :   /* goes down to 2, but coeffs are 0 in degree < v */
    1992        3388 :   for (k = lim+1; k >= v; k--) /* S <- (S + coeff(ser,k)/(k-1)) / N */
    1993        3332 :     S = gdivgs(gadd(S, gdivgs(sercoeff(ser,k), k-1)), N);
    1994          56 :   if (v-2) S = gdiv(S, powuu(N, v-2));
    1995          56 :   return gerepileupto(av, gprec_w(S, prec));
    1996             : }
    1997             : 
    1998             : /* sum_{n >= a} F(n) */
    1999             : GEN
    2000          77 : sumnumrat(GEN F, GEN ga, long prec)
    2001             : {
    2002          77 :   pari_sp av = avma;
    2003          77 :   long B = prec2nbits(prec), k, vx;
    2004             :   GEN nF, dF, S, S1, S2, intf, tmp;
    2005             :   long r, j, N, m, a;
    2006             : 
    2007          77 :   switch(typ(F))
    2008             :   {
    2009          70 :     case t_RFRAC: break;
    2010             :     case t_POL:
    2011             :     case t_INT:
    2012             :     case t_REAL:
    2013           7 :     case t_COMPLEX: if (gequal0(F)) return real_0(prec);
    2014           7 :     default: pari_err_TYPE("sumnumrat",F);
    2015             :   }
    2016          70 :   nF = gel(F,1);
    2017          70 :   dF = gel(F,2); vx = varn(dF);
    2018          70 :   switch(typ(ga))
    2019             :   {
    2020             :     case t_INT:
    2021          49 :       a = itos(ga);
    2022          49 :       if (a)
    2023             :       {
    2024          14 :         F = gsubst(F, vx, gaddgs(pol_x(vx), a));
    2025          14 :         nF = gel(F,1);
    2026          14 :         dF = gel(F,2);
    2027             :       }
    2028          49 :       break;
    2029             :     case t_INFINITY:
    2030          21 :       if (inf_get_sign(ga) == -1)
    2031             :       { /* F(x) + F(-x) */
    2032          14 :         a = 0;
    2033          14 :         F = gadd(F, gsubst(F, vx, RgX_neg(pol_x(vx))));
    2034          14 :         if (gequal0(F)) return real_0(prec);
    2035           7 :         break;
    2036             :       }
    2037             :     default:
    2038           7 :       pari_err_TYPE("sumnumrat",ga);
    2039             :       return NULL; /* LCOV_EXCL_LINE */
    2040             :   }
    2041          56 :   if (poldegree(F, -1) > -2) pari_err(e_MISC, "sum diverges in sumnumrat");
    2042          49 :   r = (long)ceil(gtodouble(ratpolemax(F)));
    2043          49 :   k = maxss(50, (long)ceil(0.241*B)); if (k&1L) k++;
    2044          49 :   tmp = gpow(gmul2n(Q_abs(bernfrac(k)), B), ginv(stoi(k)), LOWDEFAULTPREC);
    2045          49 :   N = maxss(2*r, r + 1 + itos(gceil(tmp)));
    2046          49 :   intf = intnumainfrat(F, N, prec);
    2047          49 :   S = gmul(real_1(prec), gsubst(F, vx, gaddgs(pol_x(vx), N)));
    2048          49 :   S = gtoser(S, vx, k + 2);
    2049          49 :   S1 = NULL;
    2050          49 :   if (typ(ga) == t_INFINITY)
    2051        1414 :     for (m = N; m >= -N; m--)
    2052             :     {
    2053        1407 :       GEN gm = stoi(m), d = poleval(dF, gm);
    2054        1407 :       if (gequal0(d)) continue;
    2055        1400 :       d = gdiv(poleval(nF, gm), d);
    2056        1400 :       if (labs(m) == N) d = gtofp(gmul2n(d,-1),prec);
    2057        1400 :       S1 = S1? gadd(S1, d): d;
    2058             :     }
    2059             :   else
    2060        4326 :     for (m = N; m >= 0; m--)
    2061             :     {
    2062        4284 :       GEN gm = utoi(m), d = poleval(dF, gm);
    2063        4284 :       if (gequal0(d)) continue;
    2064        4277 :       d = gdiv(poleval(nF, gm), d);
    2065        4277 :       S1 = S1? gadd(S1, d): gtofp(gmul2n(d,-1),prec);
    2066             :     }
    2067          49 :   S2 = gen_0;
    2068        2352 :   for (j = 2; j <= k; j += 2)
    2069        2303 :     S2 = gadd(S2, gdivgs(gmul(bernfrac(j), sercoeff(S, j-1)), j));
    2070          49 :   return gerepileupto(av, gadd(intf, gsub(S1, S2)));
    2071             : }
    2072             : 
    2073             : GEN
    2074          42 : sumaltrat(GEN F, GEN ga, long prec)
    2075             : {
    2076          42 :   pari_sp av = avma;
    2077             :   GEN G, res, X, X2;
    2078          42 :   long a = 0, vx;
    2079             : 
    2080          42 :   if (typ(F) != t_RFRAC) pari_err_TYPE("sumnumrat",F);
    2081          35 :   vx = varn(gel(F,2));
    2082          35 :   switch(typ(ga))
    2083             :   {
    2084             :     case t_INT:
    2085          21 :       a = itos(ga);
    2086          21 :       break;
    2087             :     case t_INFINITY:
    2088          14 :       if (inf_get_sign(ga) == -1)
    2089             :       {
    2090           7 :         F = gsub(F, gsubst(F, vx, gsubsg(-1, pol_x(vx))));
    2091           7 :         break;
    2092             :       }
    2093             :     default:
    2094           7 :       pari_err_TYPE("sumnumaltrat",ga);
    2095             :   }
    2096          28 :   X = pol_x(vx); X2 = gmul2n(X, 1);
    2097          28 :   if (a) F = gsubst(F, vx, gaddgs(X, a));
    2098          28 :   G = gsubst(F, vx, X2);
    2099          28 :   G = gsub(G, gsubst(F, vx, gaddgs(X2, 1)));
    2100          28 :   res = sumnumrat(G, gen_0, prec);
    2101          28 :   if (odd(a)) res = gneg(res);
    2102          28 :   return gerepileupto(av, res);
    2103             : }
    2104             : 
    2105             : /* prod_{n >= a} F(n) */
    2106             : GEN
    2107          14 : prodnumrat(GEN F, long a, long prec)
    2108             : {
    2109          14 :   pari_sp ltop = avma;
    2110          14 :   long B = prec2nbits(prec), k, vx;
    2111          14 :   GEN S, S1, S2, intf, G, tmp, F1 = gsubgs(F,1);
    2112             :   long r, j, N, m;
    2113             : 
    2114          14 :   if (poldegree(F1,-1) > -2) pari_err(e_MISC, "product diverges in prodnumrat");
    2115           7 :   vx = gvar(F);
    2116           7 :   if (a) F = gsubst(F, vx, gaddgs(pol_x(vx), a));
    2117           7 :   r = (long)ceil(gtodouble(ratpolemax2(F, NULL)));
    2118           7 :   k = maxss(50, (long)ceil(0.241*B)); if (k&1L) k++;
    2119           7 :   tmp = gpow(gmul2n(Q_abs(bernfrac(k)), B), ginv(stoi(k)), LOWDEFAULTPREC);
    2120           7 :   N = maxss(2*r, r + 1 + itos(gceil(tmp)));
    2121           7 :   G = gdiv(deriv(F, vx), F);
    2122           7 :   intf = intnumainfrat(gmul(pol_x(vx),G), N, prec);
    2123           7 :   intf = gneg(gadd(intf, gmulsg(N, glog(gsubst(F, vx, stoi(N)), prec))));
    2124           7 :   S = gmul(real_1(prec), gsubst(G, vx, gaddgs(pol_x(vx), N)));
    2125           7 :   S = gtoser(S, vx, k + 2);
    2126           7 :   S1 = gsqrt(gsubst(F, vx, stoi(N)), prec);
    2127           7 :   for (m = 0; m < N; ++m) S1 = gmul(S1, gsubst(F, vx, stoi(m)));
    2128           7 :   S2 = gen_0;
    2129         336 :   for (j = 2; j <= k; j += 2)
    2130         329 :     S2 = gadd(S2, gdivgs(gmul(bernfrac(j), sercoeff(S, j-2)), j*(j-1)));
    2131           7 :   return gerepileupto(ltop, gmul(S1, gexp(gsub(intf, S2), prec)));
    2132             : }
    2133             : 
    2134             : static GEN
    2135        3374 : sdmob(GEN ser, long n)
    2136             : {
    2137        3374 :   GEN D = divisorsu(n), S = gen_0;
    2138        3374 :   long i, l = lg(D);
    2139       19481 :   for (i = 1; i < l; ++i)
    2140             :   {
    2141       16107 :     long d = D[i], mob = moebiusu(d);
    2142       16107 :     if (mob) S = gadd(S, gdivgs(sercoeff(ser, n/d), mob*d));
    2143             :   }
    2144        3374 :   return S;
    2145             : }
    2146             : static GEN
    2147        3374 : logzetan(GEN s, long N, long prec)
    2148             : {
    2149             :   forprime_t T;
    2150        3374 :   GEN negs = gneg(s), Z = gzeta(s, prec);
    2151             :   long p;
    2152        3374 :   u_forprime_init(&T, 2, N);
    2153       40488 :   while ( (p = u_forprime_next(&T)) )
    2154       33740 :     Z = gmul(Z, gsubsg(1, gpow(utoipos(p), negs, prec)));
    2155        3374 :   return glog(Z, prec);
    2156             : }
    2157             : static GEN
    2158          42 : sumlogzeta(GEN ser, GEN s, long N, long vF, long lim, long prec)
    2159             : {
    2160          42 :   GEN z = gen_0;
    2161             :   long n;
    2162        3416 :   for (n = vF; n <= lim; ++n)
    2163        3374 :     z = gadd(z, gmul(logzetan(gmulsg(n, s), N, prec), sdmob(ser, n)));
    2164          42 :   return z;
    2165             : }
    2166             : 
    2167             : /* sum_{p prime, p >= a} F(p^s), F rational function */
    2168             : GEN
    2169          14 : sumeulerrat(GEN F, GEN s, long a, long prec)
    2170             : {
    2171          14 :   pari_sp av = avma;
    2172             :   forprime_t T;
    2173             :   GEN r, ser, res, rsg;
    2174             :   double rs;
    2175          14 :   long B = prec2nbits(prec), vx = gvar(F), vF, p, N, lim;
    2176             : 
    2177          14 :   if (!s) s = gen_1;
    2178          14 :   vF = -poldegree(F, -1);
    2179          14 :   rsg = real_i(s); rs = gtodouble(rsg);
    2180          14 :   if (vF <= 0 || vF*rs <= 1)
    2181           0 :     pari_err(e_MISC, "real(s) <= 1/v in sumeulerrat");
    2182          14 :   a = maxss(a, 2); r = ratpolemax(F);
    2183          14 :   if (rs <= dbllog2(r) / log2((double)a))
    2184           0 :       pari_err(e_MISC, "real(s) too small in sumeulerrat");
    2185          14 :   N = maxss(maxss(30, a), (long)ceil(2*gtodouble(r)));
    2186          14 :   lim = (long)ceil(B / dbllog2(gdiv(gpow(stoi(N), rsg, LOWDEFAULTPREC), r))) + 1;
    2187          14 :   ser = gmul(real_1(prec + EXTRAPREC), F);
    2188          14 :   ser = rfracrecip_to_ser_absolute(ser, lim + 2);
    2189          14 :   res = sumlogzeta(ser, s, N, vF, lim, prec);
    2190          14 :   u_forprime_init(&T, a, N);
    2191         168 :   while ( (p = u_forprime_next(&T)) )
    2192         140 :     res = gadd(res, gsubst(F, vx, gpow(utoipos(p), s, prec)));
    2193          14 :   return gerepileupto(av, gprec_w(res, prec));
    2194             : }
    2195             : 
    2196             : /* prod_{p prime, p >= a} F(p^s), F rational function */
    2197             : GEN
    2198          28 : prodeulerrat(GEN F, GEN s, long a, long prec)
    2199             : {
    2200          28 :   pari_sp ltop = avma;
    2201             :   forprime_t T;
    2202             :   GEN F1, r, r1, ser, res, rsg;
    2203             :   double rs;
    2204          28 :   long B = prec2nbits(prec), vx = gvar(F), vF, p, N, lim;
    2205             : 
    2206          28 :   if (!s) s = gen_1;
    2207          28 :   F1 = gsubgs(F, 1);
    2208          28 :   vF = -poldegree(F1, -1);
    2209          28 :   rsg = real_i(s); rs = gtodouble(rsg);
    2210          28 :   if (vF <= 0 || vF*rs <= 1)
    2211           0 :     pari_err(e_MISC, "real(s) <= 1/v in prodeulerrat");
    2212          28 :   r = ratpolemax2(F, &r1);
    2213          28 :   if (rs <= dbllog2(r1) / log2((double)a))
    2214           0 :     pari_err(e_MISC, "real(s) too small in prodeulerrat");
    2215          28 :   N = maxss(maxss(30, a), (long)ceil(2*gtodouble(r)));
    2216          28 :   lim = (long)ceil(B / dbllog2(gdiv(gpow(stoi(N), rsg, LOWDEFAULTPREC), r))) + 1;
    2217          28 :   ser = gmul(real_1(prec + EXTRAPREC), F1);
    2218          28 :   ser = gaddsg(1, rfracrecip_to_ser_absolute(ser, lim + 2));
    2219          28 :   ser = glog(ser, 0);
    2220          28 :   res = sumlogzeta(ser, s, N, vF, lim, prec);
    2221          28 :   res = gexp(res, prec);
    2222          28 :   u_forprime_init(&T, a, N);
    2223         336 :   while ( (p = u_forprime_next(&T)) )
    2224         280 :     res = gmul(res, gsubst(F, vx, gpow(utoipos(p), s, prec)));
    2225          28 :   return gerepileupto(ltop, gprec_w(res, prec));
    2226             : }

Generated by: LCOV version 1.11