Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - language - intnum.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.1 lcov report (development 24988-2584e74448) Lines: 1451 1487 97.6 %
Date: 2020-01-26 05:57:03 Functions: 121 122 99.2 %
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             : 
      17             : static const long EXTRAPREC =
      18             : #ifdef LONG_IS_64BIT
      19             :   1;
      20             : #else
      21             :   2;
      22             : #endif
      23             : 
      24             : static GEN
      25             : intlin(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, GEN tab, long prec);
      26             : 
      27             : /********************************************************************/
      28             : /**                NUMERICAL INTEGRATION (Romberg)                 **/
      29             : /********************************************************************/
      30             : typedef struct {
      31             :   void *E;
      32             :   GEN (*f)(void *E, GEN);
      33             : } invfun;
      34             : 
      35             : /* 1/x^2 f(1/x) */
      36             : static GEN
      37       12474 : _invf(void *E, GEN x)
      38             : {
      39       12474 :   invfun *S = (invfun*)E;
      40       12474 :   GEN y = ginv(x);
      41       12474 :   return gmul(S->f(S->E, y), gsqr(y));
      42             : }
      43             : 
      44             : /* h and s are arrays of the same length L > D. The h[i] are (decreasing)
      45             :  * step sizes, s[i] is the computed Riemann sum for step size h[i].
      46             :  * Interpolate the last D+1 values so that s ~ polynomial in h of degree D.
      47             :  * Guess that limit_{h->0} = s(0) */
      48             : static GEN
      49         105 : interp(GEN h, GEN s, long L, long bit, long D)
      50             : {
      51         105 :   pari_sp av = avma;
      52             :   long e1,e2;
      53         105 :   GEN ss = polintspec(h + L-D, s + L-D, gen_0, D+1, &e2);
      54             : 
      55         105 :   e1 = gexpo(ss);
      56         105 :   if (DEBUGLEVEL>2)
      57             :   {
      58           0 :     err_printf("romb: iteration %ld, guess: %Ps\n", L,ss);
      59           0 :     err_printf("romb: relative error < 2^-%ld [target %ld bits]\n",e1-e2,bit);
      60             :   }
      61         105 :   if (e1-e2 <= bit && (L <= 10 || e1 >= -bit)) return gc_NULL(av);
      62          70 :   return cxtoreal(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             : /* P_N(z) / P'_N(z) if flag = 0, else N! P_N(z) */
     215             : static GEN
     216       10430 : Legendreeval(long N, GEN z, GEN z2, long flag)
     217             : {
     218       10430 :   GEN u0 = z, u1 = subrs(mulur(3, z2), 1), u2;
     219             :   long n;
     220      871017 :   for (n = 2; n < N; n++)
     221             :   {
     222      860587 :     u2 = subrr(mulrr(mulur(2*n+1, z), u1), mulir(sqru(n), u0));
     223      860587 :     u0 = u1; u1 = u2;
     224             :   }
     225       10430 :   if (flag) return u1;
     226        9387 :   return divrr(mulrr(subrs(z2, 1), u1),
     227             :                mulur(N, subrr(mulrr(z, u1), mulur(N, u0))));
     228             : }
     229             : 
     230             : /* Roots of Legendre Polynomials. */
     231             : static GEN
     232        1043 : Legendreroot(long N, double dz, long bit)
     233             : {
     234        1043 :   GEN Z = cgetr(nbits2prec(bit)), z = dbltor(dz), z2;
     235        1043 :   pari_sp av = avma;
     236        1043 :   long pr, j, e = - dblexpo(1 - dz*dz), n = 1 + expu(bit + 32 - e);
     237             : 
     238        1043 :   pr = 1 + e + ((bit - e) >> n);
     239       10430 :   for (j = 1; j <= n; j++)
     240             :   {
     241        9387 :     pr = 2 * pr - e;
     242        9387 :     z = rtor(z, nbits2prec(pr));
     243        9387 :     z2 = sqrr(z);
     244        9387 :     z = subrr(z, Legendreeval(N, z, z2, 0));
     245             :   }
     246        1043 :   affrr(z, Z); set_avma(av); return Z;
     247             : }
     248             : GEN
     249          56 : intnumgaussinit(long N, long prec)
     250             : {
     251          56 :   pari_sp av = avma;
     252             :   long N2, j, k, l, bit;
     253             :   GEN V, W, F;
     254             : 
     255          56 :   prec += EXTRAPREC;
     256          56 :   bit = prec2nbits(prec);
     257          56 :   if (N <= 0)
     258             :   {
     259          14 :     N = (long)(bit * 0.2258);
     260          14 :     if (odd(N)) N++;
     261             :   }
     262          56 :   if (N == 1) retmkvec2(mkvec(gen_0), mkvec(gen_2));
     263          49 :   if (N == 2)
     264             :   {
     265           7 :     V = mkvec(divru(sqrtr(utor(3,prec)), 3));
     266           7 :     W = mkvec(gen_1); return gerepilecopy(av, mkvec2(V, W));
     267             :   }
     268          42 :   N2 = N >> 1; l = (N+3)>> 1;
     269          42 :   V = cgetg(l, t_VEC);
     270          42 :   W = cgetg(l, t_VEC); F = sqrr(mpfactr(N-1, prec));
     271          42 :   if (!odd(N)) k = 1;
     272             :   else
     273             :   {
     274           7 :     GEN c = sqrr(divrr(sqrr(mpfactr(N2, prec)), F));
     275           7 :     shiftr_inplace(c, 2*(N-1));
     276           7 :     gel(V, 1) = gen_0;
     277           7 :     gel(W, 1) = c; k = 2;
     278             :   }
     279        1085 :   for (j = 4*N2-1; j >= 3; k++, j -= 4)
     280             :   {
     281        1043 :     GEN w, z2, z = Legendreroot(N, cos(M_PI * j / (4*N+2)), bit);
     282        1043 :     pari_sp av = avma;
     283        1043 :     gel(V, k) = z; z2 = sqrr(z);
     284        1043 :     w = divrr(subsr(1, z2), sqrr(Legendreeval(N-1, z, z2, 1)));
     285        1043 :     gel(W, k) = gerepileuptoleaf(av, w);
     286             :   }
     287          42 :   W = RgV_Rg_mul(W, divri(shiftr(F, 1), sqru(N)));
     288          42 :   return gerepilecopy(av, mkvec2(V, W));
     289             : }
     290             : 
     291             : GEN
     292          77 : intnumgauss(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, GEN tab, long prec)
     293             : {
     294          77 :   pari_sp ltop = avma;
     295             :   GEN R, W, bma, bpa, S;
     296          77 :   long n, i, prec2 = prec + EXTRAPREC;
     297          77 :   if (!tab)
     298           7 :     tab = intnumgaussinit(0,prec);
     299          70 :   else if (typ(tab) != t_INT)
     300             :   {
     301          35 :     if (typ(tab) != t_VEC || lg(tab) != 3
     302          28 :         || typ(gel(tab,1)) != t_VEC
     303          28 :         || typ(gel(tab,2)) != t_VEC
     304          28 :         || lg(gel(tab,1)) != lg(gel(tab,2)))
     305           7 :       pari_err_TYPE("intnumgauss",tab);
     306             :   }
     307             :   else
     308          35 :     tab = intnumgaussinit(itos(tab),prec);
     309             : 
     310          70 :   R = gel(tab,1); n = lg(R)-1;
     311          70 :   W = gel(tab,2);
     312          70 :   a = gprec_wensure(a, prec2);
     313          70 :   b = gprec_wensure(b, prec2);
     314          70 :   bma = gmul2n(gsub(b,a), -1); /* (b-a)/2 */
     315          70 :   bpa = gadd(bma, a); /* (b+a)/2 */
     316          70 :   if (!signe(gel(R,1)))
     317             :   { /* R[1] = 0, use middle node only once */
     318          14 :     S = gmul(gel(W,1), eval(E, bpa));
     319          14 :     i = 2;
     320             :   }
     321             :   else
     322             :   {
     323          56 :     S = gen_0;
     324          56 :     i = 1;
     325             :   }
     326        1722 :   for (; i <= n; ++i)
     327             :   {
     328        1652 :     GEN h = gmul(bma, gel(R,i)); /* != 0 */
     329        1652 :     GEN P = eval(E, gadd(bpa, h));
     330        1652 :     GEN M = eval(E, gsub(bpa, h));
     331        1652 :     S = gadd(S, gmul(gel(W,i), gadd(P,M)));
     332        1652 :     S = gprec_wensure(S, prec2);
     333             :   }
     334          70 :   return gerepilecopy(ltop, gprec_wtrunc(gmul(bma,S), prec));
     335             : }
     336             : 
     337             : GEN
     338          77 : intnumgauss0(GEN a, GEN b, GEN code, GEN tab, long prec)
     339          77 : { EXPR_WRAP(code, intnumgauss(EXPR_ARG, a, b, tab, prec)); }
     340             : 
     341             : /********************************************************************/
     342             : /**                DOUBLE EXPONENTIAL INTEGRATION                  **/
     343             : /********************************************************************/
     344             : 
     345             : typedef struct _intdata {
     346             :   long eps;  /* bit accuracy of current precision */
     347             :   long l; /* table lengths */
     348             :   GEN tabx0; /* abscissa phi(0) for t = 0 */
     349             :   GEN tabw0; /* weight phi'(0) for t = 0 */
     350             :   GEN tabxp; /* table of abscissas phi(kh) for k > 0 */
     351             :   GEN tabwp; /* table of weights phi'(kh) for k > 0 */
     352             :   GEN tabxm; /* table of abscissas phi(kh) for k < 0, possibly empty */
     353             :   GEN tabwm; /* table of weights phi'(kh) for k < 0, possibly empty */
     354             :   GEN h; /* integration step */
     355             : } intdata;
     356             : 
     357             : static const long LGTAB = 8;
     358             : #define TABh(v) gel(v,1)
     359             : #define TABx0(v) gel(v,2)
     360             : #define TABw0(v) gel(v,3)
     361             : #define TABxp(v) gel(v,4)
     362             : #define TABwp(v) gel(v,5)
     363             : #define TABxm(v) gel(v,6)
     364             : #define TABwm(v) gel(v,7)
     365             : 
     366             : static int
     367       25765 : isinR(GEN z) { return is_real_t(typ(z)); }
     368             : static int
     369       23581 : isinC(GEN z)
     370       23581 : { return (typ(z) == t_COMPLEX)? isinR(gel(z,1)) && isinR(gel(z,2)): isinR(z); }
     371             : 
     372             : static int
     373        9666 : checktabsimp(GEN tab)
     374             : {
     375             :   long L, LN, LW;
     376        9666 :   if (!tab || typ(tab) != t_VEC) return 0;
     377        9666 :   if (lg(tab) != LGTAB) return 0;
     378        9666 :   if (typ(TABxp(tab)) != t_VEC) return 0;
     379        9666 :   if (typ(TABwp(tab)) != t_VEC) return 0;
     380        9666 :   if (typ(TABxm(tab)) != t_VEC) return 0;
     381        9666 :   if (typ(TABwm(tab)) != t_VEC) return 0;
     382        9666 :   L = lg(TABxp(tab)); if (lg(TABwp(tab)) != L) return 0;
     383        9666 :   LN = lg(TABxm(tab)); if (LN != 1 && LN != L) return 0;
     384        9666 :   LW = lg(TABwm(tab)); if (LW != 1 && LW != L) return 0;
     385        9666 :   return 1;
     386             : }
     387             : 
     388             : static int
     389         539 : checktabdoub(GEN tab)
     390             : {
     391             :   long L;
     392         539 :   if (typ(tab) != t_VEC) return 0;
     393         539 :   if (lg(tab) != LGTAB) return 0;
     394         539 :   L = lg(TABxp(tab));
     395         539 :   if (lg(TABwp(tab)) != L) return 0;
     396         539 :   if (lg(TABxm(tab)) != L) return 0;
     397         539 :   if (lg(TABwm(tab)) != L) return 0;
     398         539 :   return 1;
     399             : }
     400             : 
     401             : static int
     402        4672 : checktab(GEN tab)
     403             : {
     404        4672 :   if (typ(tab) != t_VEC) return 0;
     405        4672 :   if (lg(tab) != 3) return checktabsimp(tab);
     406           7 :   return checktabsimp(gel(tab,1))
     407           7 :       && checktabsimp(gel(tab,2));
     408             : }
     409             : 
     410             : /* the TUNE parameter is heuristic */
     411             : static void
     412        1169 : intinit_start(intdata *D, long m, double TUNE, long prec)
     413             : {
     414        1169 :   long l, n, bitprec = prec2nbits(prec);
     415        1169 :   double d = bitprec*LOG10_2;
     416        1169 :   GEN h, nh, pi = mppi(prec);
     417             : 
     418        1169 :   n = (long)ceil(d*log(d) / TUNE); /* heuristic */
     419             :   /* nh ~ log(2npi/log(n)) */
     420        1169 :   nh = logr_abs(divrr(mulur(2*n, pi), logr_abs(utor(n,prec))));
     421        1169 :   h = divru(nh, n);
     422        1169 :   if (m > 0) { h = gmul2n(h,-m); n <<= m; }
     423        1169 :   D->h = h;
     424        1169 :   D->eps = bitprec;
     425        1169 :   D->l = l = n+1;
     426        1169 :   D->tabxp = cgetg(l, t_VEC);
     427        1169 :   D->tabwp = cgetg(l, t_VEC);
     428        1169 :   D->tabxm = cgetg(l, t_VEC);
     429        1169 :   D->tabwm = cgetg(l, t_VEC);
     430        1169 : }
     431             : 
     432             : static GEN
     433        1169 : intinit_end(intdata *D, long pnt, long mnt)
     434             : {
     435        1169 :   GEN v = cgetg(LGTAB, t_VEC);
     436        1169 :   if (pnt < 0) pari_err_DOMAIN("intnuminit","table length","<",gen_0,stoi(pnt));
     437        1169 :   TABx0(v) = D->tabx0;
     438        1169 :   TABw0(v) = D->tabw0;
     439        1169 :   TABh(v) = D->h;
     440        1169 :   TABxp(v) = D->tabxp; setlg(D->tabxp, pnt+1);
     441        1169 :   TABwp(v) = D->tabwp; setlg(D->tabwp, pnt+1);
     442        1169 :   TABxm(v) = D->tabxm; setlg(D->tabxm, mnt+1);
     443        1169 :   TABwm(v) = D->tabwm; setlg(D->tabwm, mnt+1); return v;
     444             : }
     445             : 
     446             : /* divide by 2 in place */
     447             : static GEN
     448      463970 : divr2_ip(GEN x) { shiftr_inplace(x, -1); return x; }
     449             : 
     450             : /* phi(t)=tanh((Pi/2)sinh(t)): from -1 to 1, hence also from a to b compact
     451             :  * interval */
     452             : static GEN
     453         532 : inittanhsinh(long m, long prec)
     454             : {
     455         532 :   GEN e, ei, ek, eik, pi = mppi(prec);
     456         532 :   long k, nt = -1;
     457             :   intdata D;
     458             : 
     459         532 :   intinit_start(&D, m, 1.86, prec);
     460         532 :   D.tabx0 = real_0(prec);
     461         532 :   D.tabw0 = Pi2n(-1,prec);
     462         532 :   e = mpexp(D.h); ek = mulrr(pi, e);
     463         532 :   ei = invr(e); eik = mulrr(pi, ei);
     464      129164 :   for (k = 1; k < D.l; k++)
     465             :   {
     466             :     GEN xp, wp, ct, st, z;
     467             :     pari_sp av;
     468      129164 :     gel(D.tabxp,k) = cgetr(prec);
     469      129164 :     gel(D.tabwp,k) = cgetr(prec); av = avma;
     470      129164 :     ct = divr2_ip(addrr(ek, eik)); /* Pi ch(kh) */
     471      129164 :     st = subrr(ek, ct); /* Pi sh(kh) */
     472      129164 :     z = invr( addrs(mpexp(st), 1) );
     473      129164 :     shiftr_inplace(z, 1); if (expo(z) < -D.eps) { nt = k-1; break; }
     474      128632 :     xp = subsr(1, z);
     475      128632 :     wp = divr2_ip(mulrr(ct, subsr(1, sqrr(xp))));
     476      128632 :     affrr(xp, gel(D.tabxp,k)); mulrrz(ek, e, ek);
     477      128632 :     affrr(wp, gel(D.tabwp,k)); mulrrz(eik, ei, eik); set_avma(av);
     478             :   }
     479         532 :   return intinit_end(&D, nt, 0);
     480             : }
     481             : 
     482             : /* phi(t)=sinh(sinh(t)): from -oo to oo, slowly decreasing, at least
     483             :  * as 1/x^2. */
     484             : static GEN
     485          14 : initsinhsinh(long m, long prec)
     486             : {
     487             :   pari_sp av;
     488             :   GEN et, ct, st, ex;
     489          14 :   long k, nt = -1;
     490             :   intdata D;
     491             : 
     492          14 :   intinit_start(&D, m, 0.666, prec);
     493          14 :   D.tabx0 = real_0(prec);
     494          14 :   D.tabw0 = real_1(prec);
     495          14 :   et = ex = mpexp(D.h);
     496        8184 :   for (k = 1; k < D.l; k++)
     497             :   {
     498             :     GEN xp, wp, ext, exu;
     499        8184 :     gel(D.tabxp,k) = cgetr(prec);
     500        8184 :     gel(D.tabwp,k) = cgetr(prec); av = avma;
     501        8184 :     ct = divr2_ip(addrr(et, invr(et)));
     502        8184 :     st = subrr(et, ct);
     503        8184 :     ext = mpexp(st);
     504        8184 :     exu = invr(ext);
     505        8184 :     xp = divr2_ip(subrr(ext, exu));
     506        8184 :     wp = divr2_ip(mulrr(ct, addrr(ext, exu)));
     507        8184 :     if (expo(wp) - 2*expo(xp) < -D.eps) { nt = k-1; break; }
     508        8170 :     affrr(xp, gel(D.tabxp,k));
     509        8170 :     affrr(wp, gel(D.tabwp,k)); et = gerepileuptoleaf(av, mulrr(et, ex));
     510             :   }
     511          14 :   return intinit_end(&D, nt, 0);
     512             : }
     513             : 
     514             : /* phi(t)=2sinh(t): from -oo to oo, exponentially decreasing as exp(-x) */
     515             : static GEN
     516         126 : initsinh(long m, long prec)
     517             : {
     518             :   pari_sp av;
     519             :   GEN et, ex, eti, xp, wp;
     520         126 :   long k, nt = -1;
     521             :   intdata D;
     522             : 
     523         126 :   intinit_start(&D, m, 1.0, prec);
     524         126 :   D.tabx0 = real_0(prec);
     525         126 :   D.tabw0 = real2n(1, prec);
     526         126 :   et = ex = mpexp(D.h);
     527       38136 :   for (k = 1; k < D.l; k++)
     528             :   {
     529       38136 :     gel(D.tabxp,k) = cgetr(prec);
     530       38136 :     gel(D.tabwp,k) = cgetr(prec); av = avma;
     531       38136 :     eti = invr(et);
     532       38136 :     xp = subrr(et, eti);
     533       38136 :     wp = addrr(et, eti);
     534       38136 :     if (cmprs(xp, (long)(M_LN2*(expo(wp)+D.eps) + 1)) > 0) { nt = k-1; break; }
     535       38010 :     affrr(xp, gel(D.tabxp,k));
     536       38010 :     affrr(wp, gel(D.tabwp,k)); et = gerepileuptoleaf(av, mulrr(et, ex));
     537             :   }
     538         126 :   return intinit_end(&D, nt, 0);
     539             : }
     540             : 
     541             : /* phi(t)=exp(2sinh(t)): from 0 to oo, slowly decreasing at least as 1/x^2 */
     542             : static GEN
     543         245 : initexpsinh(long m, long prec)
     544             : {
     545             :   GEN et, ex;
     546         245 :   long k, nt = -1;
     547             :   intdata D;
     548             : 
     549         245 :   intinit_start(&D, m, 1.05, prec);
     550         245 :   D.tabx0 = real_1(prec);
     551         245 :   D.tabw0 = real2n(1, prec);
     552         245 :   ex = mpexp(D.h);
     553         245 :   et = real_1(prec);
     554      113177 :   for (k = 1; k < D.l; k++)
     555             :   {
     556             :     GEN t, eti, xp;
     557      113177 :     et = mulrr(et, ex);
     558      113177 :     eti = invr(et); t = addrr(et, eti);
     559      113177 :     xp = mpexp(subrr(et, eti));
     560      113177 :     gel(D.tabxp,k) = xp;
     561      113177 :     gel(D.tabwp,k) = mulrr(xp, t);
     562      113177 :     gel(D.tabxm,k) = invr(xp);
     563      113177 :     gel(D.tabwm,k) = mulrr(gel(D.tabxm,k), t);
     564      113177 :     if (expo(gel(D.tabxm,k)) < -D.eps) { nt = k-1; break; }
     565             :   }
     566         245 :   return intinit_end(&D, nt, nt);
     567             : }
     568             : 
     569             : /* phi(t)=exp(t-exp(-t)) : from 0 to +oo, exponentially decreasing. */
     570             : static GEN
     571         140 : initexpexp(long m, long prec)
     572             : {
     573             :   pari_sp av;
     574             :   GEN et, ex;
     575         140 :   long k, nt = -1;
     576             :   intdata D;
     577             : 
     578         140 :   intinit_start(&D, m, 1.76, prec);
     579         140 :   D.tabx0 = mpexp(real_m1(prec));
     580         140 :   D.tabw0 = gmul2n(D.tabx0, 1);
     581         140 :   et = ex = mpexp(negr(D.h));
     582       44408 :   for (k = 1; k < D.l; k++)
     583             :   {
     584             :     GEN xp, xm, wp, wm, eti, kh;
     585       44408 :     gel(D.tabxp,k) = cgetr(prec);
     586       44408 :     gel(D.tabwp,k) = cgetr(prec);
     587       44408 :     gel(D.tabxm,k) = cgetr(prec);
     588       44408 :     gel(D.tabwm,k) = cgetr(prec); av = avma;
     589       44408 :     eti = invr(et); kh = mulur(k,D.h);
     590       44408 :     xp = mpexp(subrr(kh, et));
     591       44408 :     xm = mpexp(negr(addrr(kh, eti)));
     592       44408 :     wp = mulrr(xp, addsr(1, et));
     593       44408 :     if (expo(xm) < -D.eps && cmprs(xp, (long)(M_LN2*(expo(wp)+D.eps) + 1)) > 0) { nt = k-1; break; }
     594       44268 :     wm = mulrr(xm, addsr(1, eti));
     595       44268 :     affrr(xp, gel(D.tabxp,k));
     596       44268 :     affrr(wp, gel(D.tabwp,k));
     597       44268 :     affrr(xm, gel(D.tabxm,k));
     598       44268 :     affrr(wm, gel(D.tabwm,k)); et = gerepileuptoleaf(av, mulrr(et, ex));
     599             :   }
     600         140 :   return intinit_end(&D, nt, nt);
     601             : }
     602             : 
     603             : /* phi(t)=(Pi/h)*t/(1-exp(-sinh(t))) from 0 to oo, sine oscillation */
     604             : static GEN
     605         112 : initnumsine(long m, long prec)
     606             : {
     607             :   pari_sp av;
     608         112 :   GEN invh, et, eti, ex, pi = mppi(prec);
     609         112 :   long exh, k, nt = -1;
     610             :   intdata D;
     611             : 
     612         112 :   intinit_start(&D, m, 0.666, prec);
     613         112 :   invh = invr(D.h);
     614         112 :   D.tabx0 = mulrr(pi, invh);
     615         112 :   D.tabw0 = gmul2n(D.tabx0,-1);
     616         112 :   exh = expo(invh); /*  expo(1/h) */
     617         112 :   et = ex = mpexp(D.h);
     618       90811 :   for (k = 1; k < D.l; k++)
     619             :   {
     620             :     GEN xp,xm, wp,wm, ct,st, extp,extp1,extp2, extm,extm1,extm2, kct, kpi;
     621       90811 :     gel(D.tabxp,k) = cgetr(prec);
     622       90811 :     gel(D.tabwp,k) = cgetr(prec);
     623       90811 :     gel(D.tabxm,k) = cgetr(prec);
     624       90811 :     gel(D.tabwm,k) = cgetr(prec); av = avma;
     625       90811 :     eti = invr(et); /* exp(-kh) */
     626       90811 :     ct = divr2_ip(addrr(et, eti)); /* ch(kh) */
     627       90811 :     st = divr2_ip(subrr(et, eti)); /* sh(kh) */
     628       90811 :     extp = mpexp(st);  extp1 = subsr(1, extp);
     629       90811 :     extp2 = invr(extp1); /* 1/(1-exp(sh(kh))) */
     630       90811 :     extm = invr(extp); extm1 = subsr(1, extm);
     631       90811 :     extm2 = invr(extm1);/* 1/(1-exp(sh(-kh))) */
     632       90811 :     kpi = mulur(k, pi);
     633       90811 :     kct = mulur(k, ct);
     634       90811 :     extm1 = mulrr(extm1, invh);
     635       90811 :     extp1 = mulrr(extp1, invh);
     636       90811 :     xp = mulrr(kpi, extm2); /* phi(kh) */
     637       90811 :     wp = mulrr(subrr(extm1, mulrr(kct, extm)), mulrr(pi, sqrr(extm2)));
     638       90811 :     xm = mulrr(negr(kpi), extp2); /* phi(-kh) */
     639       90811 :     wm = mulrr(addrr(extp1, mulrr(kct, extp)), mulrr(pi, sqrr(extp2)));
     640       90811 :     if (expo(wm) < -D.eps && expo(extm) + exh + expu(10 * k) < -D.eps) { nt = k-1; break; }
     641       90699 :     affrr(xp, gel(D.tabxp,k));
     642       90699 :     affrr(wp, gel(D.tabwp,k));
     643       90699 :     affrr(xm, gel(D.tabxm,k));
     644       90699 :     affrr(wm, gel(D.tabwm,k)); et = gerepileuptoleaf(av, mulrr(et, ex));
     645             :   }
     646         112 :   return intinit_end(&D, nt, nt);
     647             : }
     648             : 
     649             : /* End of initialization functions. These functions can be executed once
     650             :  * and for all for a given accuracy and type of integral ([a,b], [a,oo[ or
     651             :  * ]-oo,a], ]-oo,oo[) */
     652             : 
     653             : /* The numbers below can be changed, but NOT the ordering */
     654             : enum {
     655             :   f_REG     = 0, /* regular function */
     656             :   f_SER     = 1, /* power series */
     657             :   f_SINGSER = 2, /* algebraic singularity, power series endpoint */
     658             :   f_SING    = 3, /* algebraic singularity */
     659             :   f_YSLOW   = 4, /* oo, slowly decreasing, at least x^(-2)  */
     660             :   f_YVSLO   = 5, /* oo, very slowly decreasing, worse than x^(-2) */
     661             :   f_YFAST   = 6, /* oo, exponentially decreasing */
     662             :   f_YOSCS   = 7, /* oo, sine oscillating */
     663             :   f_YOSCC   = 8  /* oo, cosine oscillating */
     664             : };
     665             : /* is finite ? */
     666             : static int
     667         973 : is_fin_f(long c) { return c == f_REG || c == f_SER || c == f_SING; }
     668             : /* is oscillatory ? */
     669             : static int
     670         140 : is_osc(long c) { long a = labs(c); return a == f_YOSCC|| a == f_YOSCS; }
     671             : 
     672             : /* All inner functions such as intn, etc... must be called with a
     673             :  * valid 'tab' table. The wrapper intnum provides a higher level interface */
     674             : 
     675             : /* compute \int_a^b f(t)dt with [a,b] compact and f nonsingular. */
     676             : static GEN
     677        4049 : intn(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, GEN tab)
     678             : {
     679             :   GEN tabx0, tabw0, tabxp, tabwp;
     680             :   GEN bpa, bma, bmb, S;
     681             :   long i, prec;
     682        4049 :   pari_sp ltop = avma, av;
     683             : 
     684        4049 :   if (!checktabsimp(tab)) pari_err_TYPE("intnum",tab);
     685        4049 :   tabx0 = TABx0(tab); tabw0 = TABw0(tab); prec = gprecision(tabw0);
     686        4049 :   tabxp = TABxp(tab); tabwp = TABwp(tab);
     687        4049 :   bpa = gmul2n(gadd(b, a), -1); /* (b+a)/2 */
     688        4049 :   bma = gsub(bpa, a); /* (b-a)/2 */
     689        4049 :   av = avma;
     690        4049 :   bmb = gmul(bma, tabx0); /* (b-a)/2 phi(0) */
     691             :   /* phi'(0) f( (b+a)/2 + (b-a)/2 * phi(0) ) */
     692        4049 :   S = gmul(tabw0, eval(E, gadd(bpa, bmb)));
     693     1033587 :   for (i = lg(tabxp)-1; i > 0; i--)
     694             :   {
     695             :     GEN SP, SM;
     696     1029538 :     bmb = gmul(bma, gel(tabxp,i));
     697     1029538 :     SP = eval(E, gsub(bpa, bmb));
     698     1029538 :     SM = eval(E, gadd(bpa, bmb));
     699     1029538 :     S = gadd(S, gmul(gel(tabwp,i), gadd(SP, SM)));
     700     1029538 :     if ((i & 0x7f) == 1) S = gerepileupto(av, S);
     701     1029538 :     S = gprec_wensure(S, prec);
     702             :   }
     703        4049 :   return gerepileupto(ltop, gmul(S, gmul(bma, TABh(tab))));
     704             : }
     705             : 
     706             : /* compute \int_a^b f(t)dt with [a,b] compact, possible singularity with
     707             :  * exponent a[2] at lower extremity, b regular. Use tanh(sinh(t)). */
     708             : static GEN
     709         357 : intnsing(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, GEN tab)
     710             : {
     711             :   GEN tabx0, tabw0, tabxp, tabwp, ea, ba, S;
     712             :   long i, prec;
     713         357 :   pari_sp ltop = avma, av;
     714             : 
     715         357 :   if (!checktabsimp(tab)) pari_err_TYPE("intnum",tab);
     716         357 :   tabx0 = TABx0(tab); tabw0 = TABw0(tab); prec = gprecision(tabw0);
     717         357 :   tabxp = TABxp(tab); tabwp = TABwp(tab);
     718         357 :   ea = ginv(gaddsg(1, gel(a,2)));
     719         357 :   a = gel(a,1);
     720         357 :   ba = gdiv(gsub(b, a), gpow(gen_2, ea, prec));
     721         357 :   av = avma;
     722         357 :   S = gmul(gmul(tabw0, ba), eval(E, gadd(gmul(ba, addsr(1, tabx0)), a)));
     723       88669 :   for (i = lg(tabxp)-1; i > 0; i--)
     724             :   {
     725       88312 :     GEN p = addsr(1, gel(tabxp,i));
     726       88312 :     GEN m = subsr(1, gel(tabxp,i));
     727       88312 :     GEN bp = gmul(ba, gpow(p, ea, prec));
     728       88312 :     GEN bm = gmul(ba, gpow(m, ea, prec));
     729       88312 :     GEN SP = gmul(gdiv(bp, p), eval(E, gadd(bp, a)));
     730       88312 :     GEN SM = gmul(gdiv(bm, m), eval(E, gadd(bm, a)));
     731       88312 :     S = gadd(S, gmul(gel(tabwp,i), gadd(SP, SM)));
     732       88312 :     if ((i & 0x7f) == 1) S = gerepileupto(av, S);
     733       88312 :     S = gprec_wensure(S, prec);
     734             :   }
     735         357 :   return gerepileupto(ltop, gmul(gmul(S, TABh(tab)), ea));
     736             : }
     737             : 
     738      187922 : static GEN id(GEN x) { return x; }
     739             : 
     740             : /* compute  \int_a^oo f(t)dt if si>0 or \int_{-oo}^a f(t)dt if si<0$.
     741             :  * Use exp(2sinh(t)) for slowly decreasing functions, exp(1+t-exp(-t)) for
     742             :  * exponentially decreasing functions, and (pi/h)t/(1-exp(-sinh(t))) for
     743             :  * oscillating functions. */
     744             : static GEN
     745         539 : intninfpm(void *E, GEN (*eval)(void*, GEN), GEN a, long sb, GEN tab)
     746             : {
     747             :   GEN tabx0, tabw0, tabxp, tabwp, tabxm, tabwm;
     748             :   GEN S;
     749             :   long L, i, prec;
     750         539 :   pari_sp av = avma;
     751             : 
     752         539 :   if (!checktabdoub(tab)) pari_err_TYPE("intnum",tab);
     753         539 :   tabx0 = TABx0(tab); tabw0 = TABw0(tab); prec = gprecision(tabw0);
     754         539 :   tabxp = TABxp(tab); tabwp = TABwp(tab); L = lg(tabxp);
     755         539 :   tabxm = TABxm(tab); tabwm = TABwm(tab);
     756         539 :   if (gequal0(a))
     757             :   {
     758         294 :     GEN (*NEG)(GEN) = sb > 0? id: gneg;
     759         294 :     S = gmul(tabw0, eval(E, NEG(tabx0)));
     760      137578 :     for (i = 1; i < L; i++)
     761             :     {
     762      137284 :       GEN SP = eval(E, NEG(gel(tabxp,i)));
     763      137284 :       GEN SM = eval(E, NEG(gel(tabxm,i)));
     764      137284 :       S = gadd(S, gadd(gmul(gel(tabwp,i), SP), gmul(gel(tabwm,i), SM)));
     765      137284 :       if ((i & 0x7f) == 1) S = gerepileupto(av, S);
     766      137284 :       S = gprec_wensure(S, prec);
     767             :     }
     768             :   }
     769         245 :   else if (gexpo(a) <= 0 || is_osc(sb))
     770         112 :   { /* a small */
     771         112 :     GEN (*ADD)(GEN,GEN) = sb > 0? gadd: gsub;
     772         112 :     S = gmul(tabw0, eval(E, ADD(a, tabx0)));
     773       62846 :     for (i = 1; i < L; i++)
     774             :     {
     775       62734 :       GEN SP = eval(E, ADD(a, gel(tabxp,i)));
     776       62734 :       GEN SM = eval(E, ADD(a, gel(tabxm,i)));
     777       62734 :       S = gadd(S, gadd(gmul(gel(tabwp,i), SP), gmul(gel(tabwm,i), SM)));
     778       62734 :       if ((i & 0x7f) == 1) S = gerepileupto(av, S);
     779       62734 :       S = gprec_wensure(S, prec);
     780             :     }
     781             :   }
     782             :   else
     783             :   { /* a large, |a|*\int_sgn(a)^{oo} f(|a|*x)dx (sb > 0)*/
     784         133 :     GEN (*ADD)(long,GEN) = sb > 0? addsr: subsr;
     785         133 :     long sa = gsigne(a);
     786         133 :     GEN A = sa > 0? a: gneg(a);
     787         133 :     pari_sp av2 = avma;
     788         133 :     S = gmul(tabw0, eval(E, gmul(A, ADD(sa, tabx0))));
     789       88705 :     for (i = 1; i < L; i++)
     790             :     {
     791       88572 :       GEN SP = eval(E, gmul(A, ADD(sa, gel(tabxp,i))));
     792       88572 :       GEN SM = eval(E, gmul(A, ADD(sa, gel(tabxm,i))));
     793       88572 :       S = gadd(S, gadd(gmul(gel(tabwp,i), SP), gmul(gel(tabwm,i), SM)));
     794       88572 :       if ((i & 0x7f) == 1) S = gerepileupto(av2, S);
     795       88572 :       S = gprec_wensure(S, prec);
     796             :     }
     797         133 :     S = gmul(S,A);
     798             :   }
     799         539 :   return gerepileupto(av, gmul(S, TABh(tab)));
     800             : }
     801             : 
     802             : /* Compute  \int_{-oo}^oo f(t)dt
     803             :  * use sinh(sinh(t)) for slowly decreasing functions and sinh(t) for
     804             :  * exponentially decreasing functions.
     805             :  * HACK: in case TABwm(tab) contains something, assume function to be integrated
     806             :  * satisfies f(-x) = conj(f(x)).
     807             :  */
     808             : static GEN
     809         581 : intninfinf(void *E, GEN (*eval)(void*, GEN), GEN tab)
     810             : {
     811             :   GEN tabx0, tabw0, tabxp, tabwp, tabwm;
     812             :   GEN S;
     813             :   long L, i, prec, spf;
     814         581 :   pari_sp ltop = avma;
     815             : 
     816         581 :   if (!checktabsimp(tab)) pari_err_TYPE("intnum",tab);
     817         581 :   tabx0 = TABx0(tab); tabw0 = TABw0(tab); prec = gprecision(tabw0);
     818         581 :   tabxp = TABxp(tab); tabwp = TABwp(tab); L = lg(tabxp);
     819         581 :   tabwm = TABwm(tab);
     820         581 :   spf = (lg(tabwm) == lg(tabwp));
     821         581 :   S = gmul(tabw0, eval(E, tabx0));
     822         581 :   if (spf) S = gmul2n(real_i(S), -1);
     823      176099 :   for (i = L-1; i > 0; i--)
     824             :   {
     825      175518 :     GEN SP = eval(E, gel(tabxp,i));
     826      175518 :     if (spf)
     827      170044 :       S = gadd(S, real_i(gmul(gel(tabwp,i), SP)));
     828             :     else
     829             :     {
     830        5474 :       GEN SM = eval(E, negr(gel(tabxp,i)));
     831        5474 :       S = gadd(S, gmul(gel(tabwp,i), gadd(SP,SM)));
     832             :     }
     833      175518 :     if ((i & 0x7f) == 1) S = gerepileupto(ltop, S);
     834      175518 :     S = gprec_wensure(S, prec);
     835             :   }
     836         581 :   if (spf) S = gmul2n(S,1);
     837         581 :   return gerepileupto(ltop, gmul(S, TABh(tab)));
     838             : }
     839             : 
     840             : /* general num integration routine int_a^b f(t)dt, where a and b are as follows:
     841             :  - a scalar : the scalar, no singularity worse than logarithmic at a.
     842             :  - [a, e] : the scalar a, singularity exponent -1 < e <= 0.
     843             :  - +oo: slowly decreasing function (at least O(t^-2))
     844             :  - [[+oo], a], a nonnegative real : +oo, function behaving like exp(-a|t|)
     845             :  - [[+oo], e], e < -1 : +oo, function behaving like t^e
     846             :  - [[+oo], a*I], a > 0 real : +oo, function behaving like cos(at)
     847             :  - [[+oo], a*I], a < 0 real : +oo, function behaving like sin(at)
     848             :  and similarly at -oo */
     849             : static GEN
     850        2002 : f_getycplx(GEN a, long prec)
     851             : {
     852             :   GEN a2R, a2I;
     853             :   long s;
     854             : 
     855        2002 :   if (lg(a) == 2 || gequal0(gel(a,2))) return gen_1;
     856        1960 :   a2R = real_i(gel(a,2));
     857        1960 :   a2I = imag_i(gel(a,2));
     858        1960 :   s = gsigne(a2I); if (s < 0) a2I = gneg(a2I);
     859        1960 :   return ginv(gprec_wensure(s ? a2I : a2R, prec));
     860             : }
     861             : 
     862             : static void
     863          14 : err_code(GEN a, const char *name)
     864             : {
     865          14 :   char *s = stack_sprintf("intnum [incorrect %s]", name);
     866          14 :   pari_err_TYPE(s, a);
     867           0 : }
     868             : 
     869             : /* a = [[+/-oo], alpha]*/
     870             : static long
     871        4333 : code_aux(GEN a, const char *name)
     872             : {
     873        4333 :   GEN re, im, alpha = gel(a,2);
     874             :   long s;
     875        4333 :   if (!isinC(alpha)) err_code(a, name);
     876        4333 :   re = real_i(alpha);
     877        4333 :   im = imag_i(alpha);
     878        4333 :   s = gsigne(im);
     879        4333 :   if (s)
     880             :   {
     881         385 :     if (!gequal0(re)) err_code(a, name);
     882         378 :     return s > 0 ? f_YOSCC : f_YOSCS;
     883             :   }
     884        3948 :   if (gequal0(re) || gcmpgs(re, -2)<=0) return f_YSLOW;
     885        3584 :   if (gsigne(re) > 0) return f_YFAST;
     886         343 :   if (gcmpgs(re, -1) >= 0) err_code(a, name);
     887         343 :   return f_YVSLO;
     888             : }
     889             : 
     890             : static long
     891       24127 : transcode(GEN a, const char *name)
     892             : {
     893             :   GEN a1, a2;
     894       24127 :   switch(typ(a))
     895             :   {
     896        5859 :     case t_VEC: break;
     897             :     case t_INFINITY:
     898         217 :       return inf_get_sign(a) == 1 ? f_YSLOW: -f_YSLOW;
     899             :     case t_SER: case t_POL: case t_RFRAC:
     900         259 :       return f_SER;
     901       17792 :     default: if (!isinC(a)) err_code(a,name);
     902       17792 :       return f_REG;
     903             :   }
     904        5859 :   switch(lg(a))
     905             :   {
     906          21 :     case 2: return gsigne(gel(a,1)) > 0 ? f_YSLOW : -f_YSLOW;
     907        5831 :     case 3: break;
     908           7 :     default: err_code(a,name);
     909             :   }
     910        5831 :   a1 = gel(a,1);
     911        5831 :   a2 = gel(a,2);
     912        5831 :   switch(typ(a1))
     913             :   {
     914             :     case t_VEC:
     915          21 :       if (lg(a1) != 2) err_code(a,name);
     916          21 :       return gsigne(gel(a1,1)) * code_aux(a, name);
     917             :     case t_INFINITY:
     918        4312 :       return inf_get_sign(a1) * code_aux(a, name);
     919             :     case t_SER: case t_POL: case t_RFRAC:
     920          42 :       if (!isinR(a2)) err_code(a,name);
     921          42 :       if (gcmpgs(a2, -1) <= 0)
     922           0 :         pari_err_IMPL("intnum with diverging non constant limit");
     923          42 :       return gsigne(a2) < 0 ? f_SINGSER : f_SER;
     924             :     default:
     925        1456 :       if (!isinC(a1) || !isinR(a2) || gcmpgs(a2, -1) <= 0) err_code(a,name);
     926        1456 :       return gsigne(a2) < 0 ? f_SING : f_REG;
     927             :   }
     928             : }
     929             : 
     930             : /* computes the necessary tabs, knowing a, b and m */
     931             : static GEN
     932         413 : homtab(GEN tab, GEN k)
     933             : {
     934             :   GEN z;
     935         413 :   if (gequal0(k) || gequal(k, gen_1)) return tab;
     936         217 :   if (gsigne(k) < 0) k = gneg(k);
     937         217 :   z = cgetg(LGTAB, t_VEC);
     938         217 :   TABx0(z) = gmul(TABx0(tab), k);
     939         217 :   TABw0(z) = gmul(TABw0(tab), k);
     940         217 :   TABxp(z) = gmul(TABxp(tab), k);
     941         217 :   TABwp(z) = gmul(TABwp(tab), k);
     942         217 :   TABxm(z) = gmul(TABxm(tab), k);
     943         217 :   TABwm(z) = gmul(TABwm(tab), k);
     944         217 :   TABh(z) = rcopy(TABh(tab)); return z;
     945             : }
     946             : 
     947             : static GEN
     948         238 : expvec(GEN v, GEN ea, long prec)
     949             : {
     950         238 :   long lv = lg(v), i;
     951         238 :   GEN z = cgetg(lv, t_VEC);
     952         238 :   for (i = 1; i < lv; i++) gel(z,i) = gpow(gel(v,i),ea,prec);
     953         238 :   return z;
     954             : }
     955             : 
     956             : static GEN
     957      128643 : expscalpr(GEN vnew, GEN xold, GEN wold, GEN ea)
     958             : {
     959      128643 :   pari_sp av = avma;
     960      128643 :   return gerepileupto(av, gdiv(gmul(gmul(vnew, wold), ea), xold));
     961             : }
     962             : static GEN
     963         238 : expvecpr(GEN vnew, GEN xold, GEN wold, GEN ea)
     964             : {
     965         238 :   long lv = lg(vnew), i;
     966         238 :   GEN z = cgetg(lv, t_VEC);
     967      128762 :   for (i = 1; i < lv; i++)
     968      128524 :     gel(z,i) = expscalpr(gel(vnew,i), gel(xold,i), gel(wold,i), ea);
     969         238 :   return z;
     970             : }
     971             : 
     972             : /* here k < -1 */
     973             : static GEN
     974         119 : exptab(GEN tab, GEN k, long prec)
     975             : {
     976             :   GEN v, ea;
     977             : 
     978         119 :   if (gcmpgs(k, -2) <= 0) return tab;
     979         119 :   ea = ginv(gsubsg(-1, k));
     980         119 :   v = cgetg(LGTAB, t_VEC);
     981         119 :   TABx0(v) = gpow(TABx0(tab), ea, prec);
     982         119 :   TABw0(v) = expscalpr(TABx0(v), TABx0(tab), TABw0(tab), ea);
     983         119 :   TABxp(v) = expvec(TABxp(tab), ea, prec);
     984         119 :   TABwp(v) = expvecpr(TABxp(v), TABxp(tab), TABwp(tab), ea);
     985         119 :   TABxm(v) = expvec(TABxm(tab), ea, prec);
     986         119 :   TABwm(v) = expvecpr(TABxm(v), TABxm(tab), TABwm(tab), ea);
     987         119 :   TABh(v) = rcopy(TABh(tab));
     988         119 :   return v;
     989             : }
     990             : 
     991             : static GEN
     992         833 : init_fin(GEN b, long codeb, long m, long l, long prec)
     993             : {
     994         833 :   switch(labs(codeb))
     995             :   {
     996             :     case f_REG:
     997         497 :     case f_SING:  return inittanhsinh(m,l);
     998         119 :     case f_YSLOW: return initexpsinh(m,l);
     999          70 :     case f_YVSLO: return exptab(initexpsinh(m,l), gel(b,2), prec);
    1000          98 :     case f_YFAST: return homtab(initexpexp(m,l), f_getycplx(b,l));
    1001             :     /* f_YOSCS, f_YOSCC */
    1002          49 :     default: return homtab(initnumsine(m,l),f_getycplx(b,l));
    1003             :   }
    1004             : }
    1005             : 
    1006             : static GEN
    1007        1120 : intnuminit_i(GEN a, GEN b, long m, long prec)
    1008             : {
    1009             :   long codea, codeb, l;
    1010             :   GEN T, kma, kmb, tmp;
    1011             : 
    1012        1120 :   if (m > 30) pari_err_OVERFLOW("intnuminit [m]");
    1013        1120 :   if (m < 0) pari_err_DOMAIN("intnuminit", "m", "<", gen_0, stoi(m));
    1014        1113 :   l = prec+EXTRAPREC;
    1015        1113 :   codea = transcode(a, "a"); if (codea == f_SER) codea = f_REG;
    1016        1099 :   codeb = transcode(b, "b"); if (codeb == f_SER) codeb = f_REG;
    1017        1099 :   if (codea == f_SINGSER || codeb == f_SINGSER)
    1018           7 :     pari_err_IMPL("intnuminit with singularity at non constant limit");
    1019        1092 :   if (labs(codea) > labs(codeb)) { swap(a, b); lswap(codea, codeb); }
    1020        1092 :   if (codea == f_REG)
    1021             :   {
    1022         742 :     T = init_fin(b, codeb, m,l,prec);
    1023         742 :     switch(labs(codeb))
    1024             :     {
    1025          42 :       case f_YOSCS: if (gequal0(a)) break;
    1026           7 :       case f_YOSCC: T = mkvec2(inittanhsinh(m,l), T);
    1027             :     }
    1028         742 :     return T;
    1029             :   }
    1030         350 :   if (codea == f_SING)
    1031             :   {
    1032          91 :     T = init_fin(b,codeb, m,l,prec);
    1033          91 :     T = mkvec2(labs(codeb) == f_SING? T: inittanhsinh(m,l), T);
    1034          91 :     return T;
    1035             :   }
    1036             :   /* now a and b are infinite */
    1037         259 :   if (codea * codeb > 0) return gen_0;
    1038         245 :   kma = f_getycplx(a,l); codea = labs(codea);
    1039         245 :   kmb = f_getycplx(b,l); codeb = labs(codeb);
    1040         245 :   if (codea == f_YSLOW && codeb == f_YSLOW) return initsinhsinh(m, l);
    1041         231 :   if (codea == f_YFAST && codeb == f_YFAST && gequal(kma, kmb))
    1042         126 :     return homtab(initsinh(m,l), kmb);
    1043         105 :   T = cgetg(3, t_VEC);
    1044         105 :   switch (codea)
    1045             :   {
    1046             :     case f_YSLOW:
    1047             :     case f_YVSLO:
    1048          56 :       tmp = initexpsinh(m,l);
    1049          56 :       gel(T,1) = codea == f_YSLOW? tmp: exptab(tmp, gel(a,2), prec);
    1050          56 :       switch (codeb)
    1051             :       {
    1052          14 :         case f_YVSLO: gel(T,2) = exptab(tmp, gel(b,2), prec); return T;
    1053          21 :         case f_YFAST: gel(T,2) = homtab(initexpexp(m,l), kmb); return T;
    1054             :         /* YOSC[CS] */
    1055          21 :         default: gel(T,2) = homtab(initnumsine(m,l), kmb); return T;
    1056             :       }
    1057             :       break;
    1058             :     case f_YFAST:
    1059          21 :       tmp = initexpexp(m, l);
    1060          21 :       gel(T,1) = homtab(tmp, kma);
    1061          21 :       switch (codeb)
    1062             :       {
    1063           7 :         case f_YFAST: gel(T,2) = homtab(tmp, kmb); return T;
    1064             :         /* YOSC[CS] */
    1065          14 :         default: gel(T,2) = homtab(initnumsine(m, l), kmb); return T;
    1066             :       }
    1067             :     default: /* YOSC[CS] */
    1068          28 :       tmp = initnumsine(m, l);
    1069          28 :       gel(T,1) = homtab(tmp,kma);
    1070          28 :       if (codea == f_YOSCC && codeb == f_YOSCC && !gequal(kma, kmb))
    1071          14 :         gel(T,2) = mkvec2(inittanhsinh(m,l), homtab(tmp,kmb));
    1072             :       else
    1073          14 :         gel(T,2) = homtab(tmp,kmb);
    1074          28 :       return T;
    1075             :   }
    1076             : }
    1077             : GEN
    1078         980 : intnuminit(GEN a, GEN b, long m, long prec)
    1079             : {
    1080         980 :   pari_sp av = avma;
    1081         980 :   return gerepilecopy(av, intnuminit_i(a,b,m,prec));
    1082             : }
    1083             : 
    1084             : static GEN
    1085        5400 : intnuminit0(GEN a, GEN b, GEN tab, long prec)
    1086             : {
    1087             :   long m;
    1088        5400 :   if (!tab) m = 0;
    1089        4686 :   else if (typ(tab) != t_INT)
    1090             :   {
    1091        4672 :     if (!checktab(tab)) pari_err_TYPE("intnuminit0",tab);
    1092        4672 :     return tab;
    1093             :   }
    1094             :   else
    1095          14 :     m = itos(tab);
    1096         728 :   return intnuminit(a, b, m, prec);
    1097             : }
    1098             : 
    1099             : /* Assigns the values of the function weighted by w[k] at quadrature points x[k]
    1100             :  * [replacing the weights]. Return the index of the last non-zero coeff */
    1101             : static long
    1102         252 : weight(void *E, GEN (*eval)(void *, GEN), GEN x, GEN w)
    1103             : {
    1104         252 :   long k, l = lg(x);
    1105         252 :   for (k = 1; k < l; k++) gel(w,k) = gmul(gel(w,k), eval(E, gel(x,k)));
    1106         252 :   k--; while (k >= 1) if (!gequal0(gel(w,k--))) break;
    1107         252 :   return k;
    1108             : }
    1109             : /* compute the necessary tabs, weights multiplied by f(t) */
    1110             : static GEN
    1111         126 : intfuncinit_i(void *E, GEN (*eval)(void*, GEN), GEN tab)
    1112             : {
    1113         126 :   GEN tabxp = TABxp(tab), tabwp = TABwp(tab);
    1114         126 :   GEN tabxm = TABxm(tab), tabwm = TABwm(tab);
    1115         126 :   long L, L0 = lg(tabxp);
    1116             : 
    1117         126 :   TABw0(tab) = gmul(TABw0(tab), eval(E, TABx0(tab)));
    1118         126 :   if (lg(tabxm) == 1)
    1119             :   {
    1120         126 :     TABxm(tab) = tabxm = gneg(tabxp);
    1121         126 :     TABwm(tab) = tabwm = leafcopy(tabwp);
    1122             :   }
    1123             :   /* update wp and wm in place */
    1124         126 :   L = minss(weight(E, eval, tabxp, tabwp), weight(E, eval, tabxm, tabwm));
    1125         126 :   if (L < L0)
    1126             :   { /* catch up functions whose growth at oo was not adequately described */
    1127         126 :     setlg(tabxp, L+1);
    1128         126 :     setlg(tabwp, L+1);
    1129         126 :     if (lg(tabxm) > 1) { setlg(tabxm, L+1); setlg(tabwm, L+1); }
    1130             :   }
    1131         126 :   return tab;
    1132             : }
    1133             : 
    1134             : GEN
    1135         140 : intfuncinit(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, long m, long prec)
    1136             : {
    1137         140 :   pari_sp av = avma;
    1138         140 :   GEN tab = intnuminit_i(a, b, m, prec);
    1139             : 
    1140         140 :   if (lg(tab) == 3)
    1141           7 :     pari_err_IMPL("intfuncinit with hard endpoint behavior");
    1142         259 :   if (is_fin_f(transcode(a,"intfuncinit")) ||
    1143         126 :       is_fin_f(transcode(b,"intfuncinit")))
    1144           7 :     pari_err_IMPL("intfuncinit with finite endpoints");
    1145         126 :   return gerepilecopy(av, intfuncinit_i(E, eval, tab));
    1146             : }
    1147             : 
    1148             : static GEN
    1149        5400 : intnum_i(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, GEN tab, long prec)
    1150             : {
    1151        5400 :   GEN S = gen_0, kma, kmb;
    1152        5400 :   long sb, sgns = 1, codea = transcode(a, "a"), codeb = transcode(b, "b");
    1153             : 
    1154        5400 :   if (codea == f_REG && typ(a) == t_VEC) a = gel(a,1);
    1155        5400 :   if (codeb == f_REG && typ(b) == t_VEC) b = gel(b,1);
    1156        5400 :   if (codea == f_REG && codeb == f_REG) return intn(E, eval, a, b, tab);
    1157        1372 :   if (codea == f_SER || codeb == f_SER) return intlin(E, eval, a, b, tab, prec);
    1158        1302 :   if (labs(codea) > labs(codeb)) { swap(a,b); lswap(codea,codeb); sgns = -1; }
    1159             :   /* now labs(codea) <= labs(codeb) */
    1160        1302 :   if (codeb == f_SING)
    1161             :   {
    1162         266 :     if (codea == f_REG)
    1163         189 :       S = intnsing(E, eval, b, a, tab), sgns = -sgns;
    1164             :     else
    1165             :     {
    1166          77 :       GEN c = gmul2n(gadd(gel(a,1), gel(b,1)), -1);
    1167          77 :       S = gsub(intnsing(E, eval, a, c, gel(tab,1)),
    1168          77 :                intnsing(E, eval, b, c, gel(tab,2)));
    1169             :     }
    1170         266 :     return (sgns < 0) ? gneg(S) : S;
    1171             :   }
    1172             :   /* now b is infinite */
    1173        1036 :   sb = codeb > 0 ? 1 : -1;
    1174        1036 :   codeb = labs(codeb);
    1175        1036 :   if (codea == f_REG && codeb != f_YOSCC
    1176         322 :       && (codeb != f_YOSCS || gequal0(a)))
    1177             :   {
    1178         322 :     S = intninfpm(E, eval, a, sb*codeb, tab);
    1179         322 :     return sgns*sb < 0 ? gneg(S) : S;
    1180             :   }
    1181         714 :   if (is_fin_f(codea))
    1182             :   { /* either codea == f_SING  or codea == f_REG and codeb = f_YOSCC
    1183             :      * or (codeb == f_YOSCS and !gequal0(a)) */
    1184          21 :     GEN S2, c = real_i(codea == f_SING? gel(a,1): a);
    1185          21 :     switch(codeb)
    1186             :     {
    1187             :       case f_YOSCC: case f_YOSCS:
    1188             :       {
    1189           7 :         GEN pi2p = gmul(Pi2n(1,prec), f_getycplx(b, prec));
    1190           7 :         GEN pis2p = gmul2n(pi2p, -2);
    1191           7 :         if (codeb == f_YOSCC) c = gadd(c, pis2p);
    1192           7 :         c = gdiv(c, pi2p);
    1193           7 :         c = sb > 0? addiu(gceil(c), 1): subiu(gfloor(c), 1);
    1194           7 :         c = gmul(pi2p, c);
    1195           7 :         if (codeb == f_YOSCC) c = gsub(c, pis2p);
    1196           7 :         break;
    1197             :       }
    1198             :       default:
    1199          14 :         c = sb > 0? addiu(gceil(c), 1): subiu(gfloor(c), 1);
    1200          14 :         break;
    1201             :     }
    1202          35 :     S = codea==f_SING? intnsing(E, eval, a, c, gel(tab,1))
    1203          35 :                      : intn    (E, eval, a, c, gel(tab,1));
    1204          21 :     S2 = intninfpm(E, eval, c, sb*codeb, gel(tab,2));
    1205          21 :     if (sb < 0) S2 = gneg(S2);
    1206          21 :     S = gadd(S, S2);
    1207          21 :     return sgns < 0 ? gneg(S) : S;
    1208             :   }
    1209             :   /* now a and b are infinite */
    1210         693 :   if (codea * sb > 0)
    1211             :   {
    1212          14 :     if (codea > 0) pari_warn(warner, "integral from oo to oo");
    1213          14 :     if (codea < 0) pari_warn(warner, "integral from -oo to -oo");
    1214          14 :     return gen_0;
    1215             :   }
    1216         679 :   if (sb < 0) sgns = -sgns;
    1217         679 :   codea = labs(codea);
    1218         679 :   kma = f_getycplx(a, prec);
    1219         679 :   kmb = f_getycplx(b, prec);
    1220         679 :   if ((codea == f_YSLOW && codeb == f_YSLOW)
    1221         672 :    || (codea == f_YFAST && codeb == f_YFAST && gequal(kma, kmb)))
    1222         581 :     S = intninfinf(E, eval, tab);
    1223             :   else
    1224             :   {
    1225          98 :     GEN pis2 = Pi2n(-1, prec);
    1226          98 :     GEN ca = (codea == f_YOSCC)? gmul(pis2, kma): gen_0;
    1227          98 :     GEN cb = (codeb == f_YOSCC)? gmul(pis2, kmb): gen_0;
    1228          98 :     GEN c = codea == f_YOSCC ? ca : cb; /*signe(a)=-sb*/
    1229          98 :     GEN SP, SN = intninfpm(E, eval, c, -sb*codea, gel(tab,1));
    1230          98 :     if (codea != f_YOSCC)
    1231          84 :       SP = intninfpm(E, eval, cb, sb*codeb, gel(tab,2));
    1232             :     /* codea = codeb = f_YOSCC */
    1233          14 :     else if (gequal(kma, kmb))
    1234           0 :       SP = intninfpm(E, eval, cb, sb*codeb, gel(tab,2));
    1235             :     else
    1236             :     {
    1237          14 :       tab = gel(tab,2);
    1238          14 :       SP = intninfpm(E, eval, cb, sb*codeb, gel(tab,2));
    1239          14 :       SP = gadd(SP, intn(E, eval, ca, cb, gel(tab,1)));
    1240             :     }
    1241          98 :     S = gadd(SN, SP);
    1242             :   }
    1243         679 :   if (sgns < 0) S = gneg(S);
    1244         679 :   return S;
    1245             : }
    1246             : 
    1247             : GEN
    1248        5428 : intnum(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, GEN tab, long prec)
    1249             : {
    1250        5428 :   pari_sp ltop = avma;
    1251        5428 :   long l = prec+EXTRAPREC;
    1252        5428 :   GEN na = NULL, nb = NULL, S;
    1253             : 
    1254        5428 :   if (transcode(a,"a") == f_SINGSER) {
    1255          21 :     long v = gvar(gel(a,1));
    1256          21 :     if (v != NO_VARIABLE) {
    1257          21 :       na = cgetg(3,t_VEC);
    1258          21 :       gel(na,1) = polcoef(gel(a,1),0,v);
    1259          21 :       gel(na,2) = gel(a,2);
    1260             :     }
    1261          21 :     a = gel(a,1);
    1262             :   }
    1263        5428 :   if (transcode(b,"b") == f_SINGSER) {
    1264          14 :     long v = gvar(gel(b,1));
    1265          14 :     if (v != NO_VARIABLE) {
    1266          14 :       nb = cgetg(3,t_VEC);
    1267          14 :       gel(nb,1) = polcoef(gel(b,1),0,v);
    1268          14 :       gel(nb,2) = gel(b,2);
    1269             :     }
    1270          14 :     b = gel(b,1);
    1271             :   }
    1272        5428 :   if (na || nb) {
    1273          28 :     if (tab && typ(tab) != t_INT)
    1274           7 :       pari_err_IMPL("non integer tab argument");
    1275          21 :     S = intnum(E, eval, na ? na : a, nb ? nb : b, tab, prec);
    1276          21 :     if (na) S = gsub(S, intnum(E, eval, na, a, tab, prec));
    1277          21 :     if (nb) S = gsub(S, intnum(E, eval, b, nb, tab, prec));
    1278          21 :     return gerepilecopy(ltop, S);
    1279             :   }
    1280        5400 :   tab = intnuminit0(a, b, tab, prec);
    1281        5400 :   S = intnum_i(E, eval, gprec_wensure(a, l), gprec_wensure(b, l), tab, prec);
    1282        5400 :   return gerepilecopy(ltop, gprec_wtrunc(S, prec));
    1283             : }
    1284             : 
    1285             : typedef struct auxint_s {
    1286             :   GEN a, R, mult;
    1287             :   GEN (*f)(void*, GEN);
    1288             :   GEN (*w)(GEN, long);
    1289             :   long prec;
    1290             :   void *E;
    1291             : } auxint_t;
    1292             : 
    1293             : static GEN
    1294        3675 : auxcirc(void *E, GEN t)
    1295             : {
    1296        3675 :   auxint_t *D = (auxint_t*) E;
    1297             :   GEN s, c, z;
    1298        3675 :   mpsincos(mulrr(t, D->mult), &s, &c); z = mkcomplex(c,s);
    1299        3675 :   return gmul(z, D->f(D->E, gadd(D->a, gmul(D->R, z))));
    1300             : }
    1301             : 
    1302             : GEN
    1303           7 : intcirc(void *E, GEN (*eval)(void*, GEN), GEN a, GEN R, GEN tab, long prec)
    1304             : {
    1305             :   auxint_t D;
    1306             :   GEN z;
    1307             : 
    1308           7 :   D.a = a;
    1309           7 :   D.R = R;
    1310           7 :   D.mult = mppi(prec);
    1311           7 :   D.f = eval;
    1312           7 :   D.E = E;
    1313           7 :   z = intnum(&D, &auxcirc, real_m1(prec), real_1(prec), tab, prec);
    1314           7 :   return gmul2n(gmul(R, z), -1);
    1315             : }
    1316             : 
    1317             : static GEN
    1318       36750 : auxlin(void *E, GEN t)
    1319             : {
    1320       36750 :   auxint_t *D = (auxint_t*) E;
    1321       36750 :   return D->f(D->E, gadd(D->a, gmul(D->mult, t)));
    1322             : }
    1323             : 
    1324             : static GEN
    1325          70 : intlin(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, GEN tab, long prec)
    1326             : {
    1327             :   auxint_t D;
    1328             :   GEN z;
    1329             : 
    1330          70 :   if (typ(a) == t_VEC) a = gel(a,1);
    1331          70 :   if (typ(b) == t_VEC) b = gel(b,1);
    1332          70 :   z = toser_i(a); if (z) a = z;
    1333          70 :   z = toser_i(b); if (z) b = z;
    1334          70 :   D.a = a;
    1335          70 :   D.mult = gsub(b,a);
    1336          70 :   D.f = eval;
    1337          70 :   D.E = E;
    1338          70 :   z = intnum(&D, &auxlin, real_0(prec), real_1(prec), tab, prec);
    1339          70 :   return gmul(D.mult, z);
    1340             : }
    1341             : 
    1342             : GEN
    1343        4627 : intnum0(GEN a, GEN b, GEN code, GEN tab, long prec)
    1344        4627 : { EXPR_WRAP(code, intnum(EXPR_ARG, a, b, tab, prec)); }
    1345             : GEN
    1346           7 : intcirc0(GEN a, GEN R, GEN code, GEN tab, long prec)
    1347           7 : { EXPR_WRAP(code, intcirc(EXPR_ARG, a, R, tab, prec)); }
    1348             : GEN
    1349         140 : intfuncinit0(GEN a, GEN b, GEN code, long m, long prec)
    1350         140 : { EXPR_WRAP(code, intfuncinit(EXPR_ARG, a, b, m, prec)); }
    1351             : 
    1352             : #if 0
    1353             : /* Two variable integration */
    1354             : 
    1355             : typedef struct auxf_s {
    1356             :   GEN x;
    1357             :   GEN (*f)(void *, GEN, GEN);
    1358             :   void *E;
    1359             : } auxf_t;
    1360             : 
    1361             : typedef struct indi_s {
    1362             :   GEN (*c)(void*, GEN);
    1363             :   GEN (*d)(void*, GEN);
    1364             :   GEN (*f)(void *, GEN, GEN);
    1365             :   void *Ec;
    1366             :   void *Ed;
    1367             :   void *Ef;
    1368             :   GEN tabintern;
    1369             :   long prec;
    1370             : } indi_t;
    1371             : 
    1372             : static GEN
    1373             : auxf(GEN y, void *E)
    1374             : {
    1375             :   auxf_t *D = (auxf_t*) E;
    1376             :   return D->f(D->E, D->x, y);
    1377             : }
    1378             : 
    1379             : static GEN
    1380             : intnumdoubintern(GEN x, void *E)
    1381             : {
    1382             :   indi_t *D = (indi_t*) E;
    1383             :   GEN c = D->c(x, D->Ec), d = D->d(x, D->Ed);
    1384             :   auxf_t A;
    1385             : 
    1386             :   A.x = x;
    1387             :   A.f = D->f;
    1388             :   A.E = D->Ef;
    1389             :   return intnum(&A, &auxf, c, d, D->tabintern, D->prec);
    1390             : }
    1391             : 
    1392             : GEN
    1393             : 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)
    1394             : {
    1395             :   indi_t E;
    1396             : 
    1397             :   E.c = evalc;
    1398             :   E.d = evald;
    1399             :   E.f = evalf;
    1400             :   E.Ec = Ec;
    1401             :   E.Ed = Ed;
    1402             :   E.Ef = Ef;
    1403             :   E.prec = prec;
    1404             :   if (typ(tabint) == t_INT)
    1405             :   {
    1406             :     GEN C = evalc(a, Ec), D = evald(a, Ed);
    1407             :     if (typ(C) != t_VEC && typ(D) != t_VEC) { C = gen_0; D = gen_1; }
    1408             :     E.tabintern = intnuminit0(C, D, tabint, prec);
    1409             :   }
    1410             :   else E.tabintern = tabint;
    1411             :   return intnum(&E, &intnumdoubintern, a, b, tabext, prec);
    1412             : }
    1413             : 
    1414             : GEN
    1415             : intnumdoub0(GEN a, GEN b, int nc, int nd, int nf, GEN tabext, GEN tabint, long prec)
    1416             : {
    1417             :   GEN z;
    1418             :   push_lex(NULL);
    1419             :   push_lex(NULL);
    1420             :   z = intnumdoub(chf, &gp_eval2, chc, &gp_eval, chd, &gp_eval, a, b, tabext, tabint, prec);
    1421             :   pop_lex(1); pop_lex(1); return z;
    1422             : }
    1423             : #endif
    1424             : 
    1425             : 
    1426             : /* The quotient-difference algorithm. Given a vector M, convert the series
    1427             :  * S = \sum_{n >= 0} M[n+1]z^n into a continued fraction.
    1428             :  * Compute the c[n] such that
    1429             :  * S = c[1] / (1 + c[2]z / (1+c[3]z/(1+...c[lim]z))),
    1430             :  * Compute A[n] and B[n] such that
    1431             :  * 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)))),
    1432             :  * Assume lim <= #M.
    1433             :  * Does not work for certain M. */
    1434             : 
    1435             : /* Given a continued fraction CF output by the quodif program,
    1436             : convert it into an Euler continued fraction A(n), B(n), where
    1437             : $1/(1+c[2]z/(1+c[3]z/(1+..c[lim]z)))
    1438             : =1/(1+A[1]*z+B[1]*z^2/(1+A[2]*z+B[2]*z^2/(1+...1/(1+A[lim\2]*z)))). */
    1439             : static GEN
    1440        3731 : contfrac_Euler(GEN CF)
    1441             : {
    1442        3731 :   long lima, limb, i, lim = lg(CF)-1;
    1443             :   GEN A, B;
    1444        3731 :   lima = lim/2;
    1445        3731 :   limb = (lim - 1)/2;
    1446        3731 :   A = cgetg(lima+1, t_VEC);
    1447        3731 :   B = cgetg(limb+1, t_VEC);
    1448        3731 :   gel (A, 1) = gel(CF, 2);
    1449        3731 :   for (i=2; i <= lima; ++i) gel(A,i) = gadd(gel(CF, 2*i), gel(CF, 2*i-1));
    1450        3731 :   for (i=1; i <= limb; ++i) gel(B,i) = gneg(gmul(gel(CF, 2*i+1), gel(CF, 2*i)));
    1451        3731 :   return mkvec2(A, B);
    1452             : }
    1453             : 
    1454             : static GEN
    1455        4004 : contfracinit_i(GEN M, long lim)
    1456             : {
    1457             :   pari_sp av;
    1458             :   GEN e, q, c;
    1459             :   long lim2, j, k;
    1460        4004 :   e = zerovec(lim);
    1461        4004 :   c = zerovec(lim+1); gel(c, 1) = gel(M, 1);
    1462        4004 :   q = cgetg(lim+1, t_VEC);
    1463        4004 :   for (k = 1; k <= lim; ++k) gel(q, k) = gdiv(gel(M, k+1), gel(M, k));
    1464        4004 :   lim2 = lim/2; av = avma;
    1465      137885 :   for (j = 1; j <= lim2; ++j)
    1466             :   {
    1467      133881 :     long l = lim - 2*j;
    1468      133881 :     gel(c, 2*j) = gneg(gel(q, 1));
    1469    11305306 :     for (k = 0; k <= l; ++k)
    1470    11171425 :       gel(e, k+1) = gsub(gadd(gel(e, k+2), gel(q, k+2)), gel(q, k+1));
    1471    11171425 :     for (k = 0; k < l; ++k)
    1472    11037544 :       gel(q, k+1) = gdiv(gmul(gel(q, k+2), gel(e, k+2)), gel(e, k+1));
    1473      133881 :     gel(c, 2*j+1) = gneg(gel(e, 1));
    1474      133881 :     if (gc_needed(av, 3))
    1475             :     {
    1476          98 :       if (DEBUGMEM>1) pari_warn(warnmem,"contfracinit, %ld/%ld",j,lim2);
    1477          98 :       gerepileall(av, 3, &e, &c, &q);
    1478             :     }
    1479             :   }
    1480        4004 :   if (odd(lim)) gel(c, lim+1) = gneg(gel(q, 1));
    1481        4004 :   return c;
    1482             : }
    1483             : 
    1484             : GEN
    1485        3752 : contfracinit(GEN M, long lim)
    1486             : {
    1487        3752 :   pari_sp ltop = avma;
    1488             :   GEN c;
    1489        3752 :   switch(typ(M))
    1490             :   {
    1491             :     case t_RFRAC:
    1492           7 :       if (lim < 0) pari_err_TYPE("contfracinit",M);
    1493           7 :       M = gtoser(M, varn(gel(M,2)), lim+3); /*fall through*/
    1494          35 :     case t_SER: M = gtovec(M); break;
    1495           7 :     case t_POL: M = RgX_to_RgC(M, degpol(M)+1); break;
    1496        3703 :     case t_VEC: case t_COL: break;
    1497           7 :     default: pari_err_TYPE("contfracinit", M);
    1498             :   }
    1499        3745 :   if (lim < 0)
    1500             :   {
    1501          28 :     lim = lg(M)-2;
    1502          28 :     if (lim < 0) retmkvec2(cgetg(1,t_VEC),cgetg(1,t_VEC));
    1503             :   }
    1504        3717 :   else if (lg(M)-1 <= lim)
    1505           0 :     pari_err_COMPONENT("contfracinit", "<", stoi(lg(M)-1), stoi(lim));
    1506        3731 :   c = contfracinit_i(M, lim);
    1507        3731 :   return gerepilecopy(ltop, contfrac_Euler(c));
    1508             : }
    1509             : 
    1510             : /* Evaluate at 1/tinv the nlim first terms of the continued fraction output by
    1511             :  * contfracinit. */
    1512             : /* Not stack clean */
    1513             : GEN
    1514     3000225 : contfraceval_inv(GEN CF, GEN tinv, long nlim)
    1515             : {
    1516             :   pari_sp btop;
    1517             :   long j;
    1518     3000225 :   GEN S = gen_0, S1, S2, A, B;
    1519     3000225 :   if (typ(CF) != t_VEC || lg(CF) != 3) pari_err_TYPE("contfraceval", CF);
    1520     3000243 :   A = gel(CF, 1); if (typ(A) != t_VEC) pari_err_TYPE("contfraceval", CF);
    1521     3000243 :   B = gel(CF, 2); if (typ(B) != t_VEC) pari_err_TYPE("contfraceval", CF);
    1522     3000243 :   if (nlim < 0)
    1523          14 :     nlim = lg(A)-1;
    1524     3000229 :   else if (lg(A) <= nlim)
    1525           7 :     pari_err_COMPONENT("contfraceval", ">", stoi(lg(A)-1), stoi(nlim));
    1526     3002268 :   if (lg(B)+1 <= nlim)
    1527           0 :     pari_err_COMPONENT("contfraceval", ">", stoi(lg(B)), stoi(nlim));
    1528     3002268 :   btop = avma;
    1529     3002268 :   if (nlim <= 1) return lg(A)==1? gen_0: gdiv(tinv, gadd(gel(A, 1), tinv));
    1530     2911515 :   switch(nlim % 3)
    1531             :   {
    1532             :     case 2:
    1533     1045175 :       S = gdiv(gel(B, nlim-1), gadd(gel(A, nlim), tinv));
    1534     1045193 :       nlim--; break;
    1535             : 
    1536             :     case 0:
    1537      994643 :       S1 = gadd(gel(A, nlim), tinv);
    1538      994455 :       S2 = gadd(gmul(gadd(gel(A, nlim-1), tinv), S1), gel(B, nlim-1));
    1539      994375 :       S = gdiv(gmul(gel(B, nlim-2), S1), S2);
    1540      994695 :       nlim -= 2; break;
    1541             :   }
    1542             :   /* nlim = 1 (mod 3) */
    1543    13540274 :   for (j = nlim; j >= 4; j -= 3)
    1544             :   {
    1545             :     GEN S3;
    1546    10630467 :     S1 = gadd(gadd(gel(A, j), tinv), S);
    1547    10591776 :     S2 = gadd(gmul(gadd(gel(A, j-1), tinv), S1), gel(B, j-1));
    1548    10587947 :     S3 = gadd(gmul(gadd(gel(A, j-2), tinv), S2), gmul(gel(B, j-2), S1));
    1549    10588912 :     S = gdiv(gmul(gel(B, j-3), S2), S3);
    1550    10628600 :     if (gc_needed(btop, 3)) S = gerepilecopy(btop, S);
    1551             :   }
    1552     2909807 :   return gdiv(tinv, gadd(gadd(gel(A, 1), tinv), S));
    1553             : }
    1554             : 
    1555             : GEN
    1556          35 : contfraceval(GEN CF, GEN t, long nlim)
    1557             : {
    1558          35 :   pari_sp ltop = avma;
    1559          35 :   return gerepileupto(ltop, contfraceval_inv(CF, ginv(t), nlim));
    1560             : }
    1561             : 
    1562             : /* MONIEN SUMMATION */
    1563             : 
    1564             : /* basic Newton, find x ~ z such that Q(x) = 0 */
    1565             : static GEN
    1566        2352 : monrefine(GEN Q, GEN QP, GEN z, long prec)
    1567             : {
    1568        2352 :   pari_sp av = avma;
    1569        2352 :   GEN pr = poleval(Q, z);
    1570             :   for(;;)
    1571        8736 :   {
    1572             :     GEN prnew;
    1573       11088 :     z = gsub(z, gdiv(pr, poleval(QP, z)));
    1574       11088 :     prnew = poleval(Q, z);
    1575       11088 :     if (gcmp(gabs(prnew, prec), gabs(pr, prec)) >= 0) break;
    1576        8736 :     pr = prnew;
    1577             :   }
    1578        2352 :   z = gprec_wensure(z, 2*prec-2);
    1579        2352 :   z = gsub(z, gdiv(poleval(Q, z), poleval(QP, z)));
    1580        2352 :   return gerepileupto(av, z);
    1581             : }
    1582             : 
    1583             : static GEN
    1584         273 : RX_realroots(GEN x, long prec)
    1585         273 : { return realroots(gprec_wtrunc(x,prec), NULL, prec); }
    1586             : 
    1587             : /* (real) roots of Q, assuming QP = Q' and that half the roots are close to
    1588             :  * k+1, ..., k+m, m = deg(Q)/2-1. N.B. All roots are real and >= 1 */
    1589             : static GEN
    1590         175 : monroots(GEN Q, GEN QP, long k, long prec)
    1591             : {
    1592         175 :   long j, n = degpol(Q), m = n/2 - 1;
    1593         175 :   GEN v2, v1 = cgetg(m+1, t_VEC);
    1594         175 :   for (j = 1; j <= m; ++j) gel(v1, j) = monrefine(Q, QP, stoi(k+j), prec);
    1595         175 :   Q = gdivent(Q, roots_to_pol(v1, varn(Q)));
    1596         175 :   v2 = RX_realroots(Q, prec); settyp(v2, t_VEC);
    1597         175 :   return shallowconcat(v1, v2);
    1598             : }
    1599             : 
    1600             : static void
    1601         273 : Pade(GEN M, GEN *pP, GEN *pQ)
    1602             : {
    1603         273 :   pari_sp av = avma;
    1604         273 :   long n = lg(M)-2, i;
    1605         273 :   GEN v = contfracinit_i(M, n), P = pol_0(0), Q = pol_1(0);
    1606             :   /* evaluate continued fraction => Pade approximants */
    1607       16289 :   for (i = n-1; i >= 1; i--)
    1608             :   { /* S = P/Q: S -> v[i]*x / (1+S) */
    1609       16016 :     GEN R = RgX_shift_shallow(RgX_Rg_mul(Q,gel(v,i)), 1);
    1610       16016 :     Q = RgX_add(P,Q); P = R;
    1611       16016 :     if (gc_needed(av, 3))
    1612             :     {
    1613           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Pade, %ld/%ld",i,n-1);
    1614           0 :       gerepileall(av, 3, &P, &Q, &v);
    1615             :     }
    1616             :   }
    1617             :   /* S -> 1+S */
    1618         273 :   *pP = RgX_add(P,Q);
    1619         273 :   *pQ = Q;
    1620         273 : }
    1621             : 
    1622             : static GEN
    1623           7 : veczetaprime(GEN a, GEN b, long N, long prec)
    1624             : {
    1625           7 :   long newprec, fpr = prec2nbits(prec), pr = (long)ceil(fpr * 1.5);
    1626           7 :   long l = nbits2prec(pr), e = fpr / 2;
    1627             :   GEN eps, A, B;
    1628           7 :   newprec = nbits2prec(pr + BITS_IN_LONG);
    1629           7 :   a = gprec_wensure(a, newprec);
    1630           7 :   b = gprec_wensure(b, newprec);
    1631           7 :   eps = real2n(-e, l);
    1632           7 :   A = veczeta(a, gsub(b, eps), N, newprec);
    1633           7 :   B = veczeta(a, gadd(b, eps), N, newprec);
    1634           7 :   return gmul2n(RgV_sub(B, A), e-1);
    1635             : }
    1636             : 
    1637             : struct mon_w {
    1638             :   GEN w, a, b;
    1639             :   long n, j, prec;
    1640             : };
    1641             : 
    1642             : /* w(x) / x^(a*(j+k)+b), k >= 1; w a t_CLOSURE or t_INT [encodes log(x)^w] */
    1643             : static GEN
    1644       34384 : wrapmonw(void* E, GEN x)
    1645             : {
    1646       34384 :   struct mon_w *W = (struct mon_w*)E;
    1647       34384 :   long k, j = W->j, n = W->n, prec = W->prec, l = 2*n+4-j;
    1648      103152 :   GEN wx = typ(W->w) == t_CLOSURE? closure_callgen1prec(W->w, x, prec)
    1649       68768 :                                  : powgi(glog(x, prec), W->w);
    1650       34384 :   GEN v = cgetg(l, t_VEC);
    1651       34384 :   GEN xa = gpow(x, gneg(W->a), prec), w = gmul(wx, gpowgs(xa, j));
    1652       34384 :   w = gdiv(w, gpow(x,W->b,prec));
    1653       34384 :   for (k = 1; k < l; k++) { gel(v,k) = w; w = gmul(w, xa); }
    1654       34384 :   return v;
    1655             : }
    1656             : /* w(x) / x^(a*j+b) */
    1657             : static GEN
    1658       18819 : wrapmonw2(void* E, GEN x)
    1659             : {
    1660       18819 :   struct mon_w *W = (struct mon_w*)E;
    1661       18819 :   GEN wnx = closure_callgen1prec(W->w, x, W->prec);
    1662       18819 :   return gdiv(wnx, gpow(x, gadd(gmulgs(W->a, W->j), W->b), W->prec));
    1663             : }
    1664             : static GEN
    1665          35 : M_from_wrapmon(struct mon_w *S, GEN wfast, GEN n0)
    1666             : {
    1667          35 :   long j, N = 2*S->n+2;
    1668          35 :   GEN M = cgetg(N+1, t_VEC), faj = gsub(wfast, S->b);
    1669          42 :   for (j = 1; j <= N; j++)
    1670             :   {
    1671          42 :     faj = gsub(faj, S->a);
    1672          42 :     if (gcmpgs(faj, -2) <= 0)
    1673             :     {
    1674          28 :       S->j = j; setlg(M,j);
    1675          28 :       M = shallowconcat(M, sumnum((void*)S, wrapmonw, n0, NULL, S->prec));
    1676          28 :       break;
    1677             :     }
    1678          14 :     S->j = j;
    1679          14 :     gel(M,j) = sumnum((void*)S, wrapmonw2, mkvec2(n0,faj), NULL, S->prec);
    1680             :   }
    1681          28 :   return M;
    1682             : }
    1683             : 
    1684             : static void
    1685         224 : checkmonroots(GEN vr, long n)
    1686             : {
    1687         224 :   if (lg(vr) != n+1)
    1688           0 :     pari_err_IMPL("recovery when missing roots in sumnummonieninit");
    1689         224 : }
    1690             : 
    1691             : static GEN
    1692         231 : sumnummonieninit_i(GEN a, GEN b, GEN w, GEN n0, long prec)
    1693             : {
    1694         231 :   GEN c, M, P, Q, Qp, vr, vabs, vwt, ga = gadd(a, b);
    1695         231 :   double bit = 2*prec2nbits(prec) / gtodouble(ga), D = bit*M_LN2;
    1696         231 :   double da = maxdd(1., gtodouble(a));
    1697         231 :   long n = (long)ceil(D / (da*(log(D)-1)));
    1698         231 :   long j, prec2, prec0 = prec + EXTRAPREC;
    1699         231 :   double bit0 = ceil((2*n+1)*LOG2_10);
    1700         231 :   int neg = 1;
    1701             :   struct mon_w S;
    1702             : 
    1703             :   /* 2.05 is heuristic; with 2.03, sumnummonien(n=1,1/n^2) loses
    1704             :    * 19 decimals at \p1500 */
    1705         231 :   prec = nbits2prec(maxdd(2.05*da*bit, bit0));
    1706         231 :   prec2 = nbits2prec(maxdd(1.3*da*bit, bit0));
    1707         231 :   S.w = w;
    1708         231 :   S.a = a = gprec_wensure(a, 2*prec-2);
    1709         231 :   S.b = b = gprec_wensure(b, 2*prec-2);
    1710         231 :   S.n = n;
    1711         231 :   S.j = 1;
    1712         231 :   S.prec = prec;
    1713         231 :   if (typ(w) == t_INT)
    1714             :   { /* f(n) ~ \sum_{i > 0} f_i log(n)^k / n^(a*i + b); a > 0, a+b > 1 */
    1715         196 :     long k = itos(w);
    1716         196 :     if (k == 0)
    1717         189 :       M = veczeta(a, ga, 2*n+2, prec);
    1718           7 :     else if (k == 1)
    1719           7 :       M = veczetaprime(a, ga, 2*n+2, prec);
    1720             :     else
    1721           0 :       M = M_from_wrapmon(&S, gen_0, gen_1);
    1722         196 :     if (odd(k)) neg = 0;
    1723             :   }
    1724             :   else
    1725             :   {
    1726          35 :     GEN wfast = gen_0;
    1727          35 :     if (typ(w) == t_VEC) { wfast = gel(w,2); w = gel(w,1); }
    1728          35 :     M = M_from_wrapmon(&S, wfast, n0);
    1729             :   }
    1730             :   /* M[j] = sum(n >= n0, w(n) / n^(a*j+b) */
    1731         224 :   Pade(M, &P,&Q);
    1732         224 :   Qp = RgX_deriv(Q);
    1733         224 :   if (gequal1(a)) a = NULL;
    1734         224 :   if (!a && typ(w) == t_INT)
    1735             :   {
    1736         175 :     vabs = vr = monroots(Q, Qp, signe(w)? 1: 0, prec2);
    1737         175 :     checkmonroots(vr, n);
    1738         175 :     c = b;
    1739             :   }
    1740             :   else
    1741             :   {
    1742          49 :     vr = RX_realroots(Q, prec2); settyp(vr, t_VEC);
    1743          49 :     checkmonroots(vr, n);
    1744          49 :     if (!a) { vabs = vr; c = b; }
    1745             :     else
    1746             :     {
    1747          35 :       GEN ai = ginv(a);
    1748          35 :       vabs = cgetg(n+1, t_VEC);
    1749          35 :       for (j = 1; j <= n; j++) gel(vabs,j) = gpow(gel(vr,j), ai, prec2);
    1750          35 :       c = gdiv(b,a);
    1751             :     }
    1752             :   }
    1753             : 
    1754         224 :   vwt = cgetg(n+1, t_VEC);
    1755         224 :   c = gsubgs(c,1); if (gequal0(c)) c = NULL;
    1756        6783 :   for (j = 1; j <= n; j++)
    1757             :   {
    1758        6559 :     GEN r = gel(vr,j), t = gdiv(poleval(P,r), poleval(Qp,r));
    1759        6559 :     if (c) t = gmul(t, gpow(r, c, prec));
    1760        6559 :     gel(vwt,j) = neg? gneg(t): t;
    1761             :   }
    1762         224 :   if (typ(w) == t_INT && !equali1(n0))
    1763             :   {
    1764          84 :     GEN h = subiu(n0,1);
    1765          84 :     for (j = 1; j <= n; j++) gel(vabs,j) = gadd(gel(vabs,j), h);
    1766             :   }
    1767         224 :   return mkvec3(gprec_wtrunc(vabs,prec0), gprec_wtrunc(vwt,prec0), n0);
    1768             : }
    1769             : 
    1770             : GEN
    1771         168 : sumnummonieninit(GEN asymp, GEN w, GEN n0, long prec)
    1772             : {
    1773         168 :   pari_sp av = avma;
    1774         168 :   const char *fun = "sumnummonieninit";
    1775             :   GEN a, b;
    1776         168 :   if (!n0) n0 = gen_1; else if (typ(n0) != t_INT) pari_err_TYPE(fun, n0);
    1777         168 :   if (!asymp) a = b = gen_1;
    1778             :   else
    1779             :   {
    1780         140 :     if (typ(asymp) == t_VEC)
    1781             :     {
    1782          70 :       if (lg(asymp) != 3) pari_err_TYPE(fun, asymp);
    1783          70 :       a = gel(asymp,1);
    1784          70 :       b = gel(asymp,2);
    1785             :     }
    1786             :     else
    1787             :     {
    1788          70 :       a = gen_1;
    1789          70 :       b = asymp;
    1790             :     }
    1791         140 :     if (gsigne(a) <= 0) pari_err_DOMAIN(fun, "a", "<=", gen_0, a);
    1792         133 :     if (!isinR(b)) pari_err_TYPE(fun, b);
    1793         126 :     if (gcmpgs(gadd(a,b), 1) <= 0)
    1794           7 :       pari_err_DOMAIN(fun, "a+b", "<=", gen_1, mkvec2(a,b));
    1795             :   }
    1796         147 :   if (!w) w = gen_0;
    1797          42 :   else switch(typ(w))
    1798             :   {
    1799             :     case t_INT:
    1800           7 :       if (signe(w) < 0) pari_err_IMPL("log power < 0 in sumnummonieninit");
    1801          35 :     case t_CLOSURE: break;
    1802             :     case t_VEC:
    1803           7 :       if (lg(w) == 3 && typ(gel(w,1)) == t_CLOSURE) break;
    1804           0 :     default: pari_err_TYPE(fun, w);
    1805             :   }
    1806         147 :   return gerepilecopy(av, sumnummonieninit_i(a, b, w, n0, prec));
    1807             : }
    1808             : 
    1809             : GEN
    1810         231 : sumnummonien(void *E, GEN (*eval)(void*,GEN), GEN n0, GEN tab, long prec)
    1811             : {
    1812         231 :   pari_sp av = avma;
    1813             :   GEN vabs, vwt, S;
    1814             :   long l, i;
    1815         231 :   if (typ(n0) != t_INT) pari_err_TYPE("sumnummonien", n0);
    1816         231 :   if (!tab)
    1817          84 :     tab = sumnummonieninit_i(gen_1, gen_1, gen_0, n0, prec);
    1818             :   else
    1819             :   {
    1820         147 :     if (lg(tab) != 4 || typ(tab) != t_VEC) pari_err_TYPE("sumnummonien", tab);
    1821         147 :     if (!equalii(n0, gel(tab,3)))
    1822           7 :       pari_err(e_MISC, "incompatible initial value %Ps != %Ps", gel(tab,3),n0);
    1823             :   }
    1824         224 :   vabs= gel(tab,1); l = lg(vabs);
    1825         224 :   vwt = gel(tab,2);
    1826         224 :   if (typ(vabs) != t_VEC || typ(vwt) != t_VEC || lg(vwt) != l)
    1827           0 :     pari_err_TYPE("sumnummonien", tab);
    1828         224 :   S = gen_0;
    1829        6783 :   for (i = 1; i < l; i++)
    1830             :   {
    1831        6559 :     S = gadd(S, gmul(gel(vwt,i), eval(E, gel(vabs,i))));
    1832        6559 :     S = gprec_wensure(S, prec);
    1833             :   }
    1834         224 :   return gerepilecopy(av, gprec_wtrunc(S, prec));
    1835             : }
    1836             : 
    1837             : static GEN
    1838         196 : get_oo(GEN fast) { return mkvec2(mkoo(), fast); }
    1839             : 
    1840             : GEN
    1841         119 : sumnuminit(GEN fast, long prec)
    1842             : {
    1843             :   pari_sp av;
    1844         119 :   GEN s, v, d, C, res = cgetg(6, t_VEC);
    1845         119 :   long bitprec = prec2nbits(prec), N, k, k2, m;
    1846             :   double w;
    1847             : 
    1848         119 :   gel(res, 1) = d = mkfrac(gen_1, utoipos(4)); /* 1/4 */
    1849         119 :   av = avma;
    1850         119 :   w = gtodouble(glambertW(ginv(d), LOWDEFAULTPREC));
    1851         119 :   N = (long)ceil(M_LN2*bitprec/(w*(1+w))+5);
    1852         119 :   k = (long)ceil(N*w); if (k&1) k--;
    1853             : 
    1854         119 :   prec += EXTRAPREC;
    1855         119 :   k2 = k/2;
    1856         119 :   s = RgX_to_ser(monomial(real_1(prec),1,0), k+3);
    1857         119 :   s = gmul2n(gasinh(s, prec), 2); /* asinh(x)/d, d = 1/4 */
    1858         119 :   gel(s, 2) = utoipos(4);
    1859         119 :   s = gsub(ser_inv(gexpm1(s,prec)), ser_inv(s));
    1860         119 :   C = matpascal(k-1);
    1861         119 :   v = cgetg(k2+1, t_VEC);
    1862        8449 :   for (m = 1; m <= k2; m++)
    1863             :   {
    1864        8330 :     pari_sp av = avma;
    1865        8330 :     GEN S = real_0(prec);
    1866             :     long j;
    1867      484169 :     for (j = m; j <= k2; j++)
    1868             :     { /* s[X^(2j-1)] * binomial(2*j-1, j-m) */
    1869      475839 :       GEN t = gmul(gel(s,2*j+1), gcoeff(C, 2*j,j-m+1));
    1870      475839 :       t = gmul2n(t, 1-2*j);
    1871      475839 :       S = odd(j)? gsub(S,t): gadd(S,t);
    1872             :     }
    1873        8330 :     if (odd(m)) S = gneg(S);
    1874        8330 :     gel(v,m) = gerepileupto(av, S);
    1875             :   }
    1876         119 :   v = RgC_gtofp(v,prec); settyp(v, t_VEC);
    1877         119 :   gel(res, 4) = gerepileupto(av, v);
    1878         119 :   gel(res, 2) = utoi(N);
    1879         119 :   gel(res, 3) = utoi(k);
    1880         119 :   if (!fast) fast = get_oo(gen_0);
    1881         119 :   gel(res, 5) = intnuminit(gel(res,2), fast, 0, prec - EXTRAPREC);
    1882         119 :   return res;
    1883             : }
    1884             : 
    1885             : static int
    1886          28 : checksumtab(GEN T)
    1887             : {
    1888          28 :   if (typ(T) != t_VEC || lg(T) != 6) return 0;
    1889          21 :   return typ(gel(T,2))==t_INT && typ(gel(T,3))==t_INT && typ(gel(T,4))==t_VEC;
    1890             : }
    1891             : GEN
    1892         133 : sumnum(void *E, GEN (*eval)(void*, GEN), GEN a, GEN tab, long prec)
    1893             : {
    1894         133 :   pari_sp av = avma, av2;
    1895             :   GEN v, tabint, S, d, fast;
    1896             :   long as, N, k, m, prec2;
    1897         133 :   if (!a) { a = gen_1; fast = get_oo(gen_0); }
    1898         133 :   else switch(typ(a))
    1899             :   {
    1900             :   case t_VEC:
    1901          49 :     if (lg(a) != 3) pari_err_TYPE("sumnum", a);
    1902          49 :     fast = get_oo(gel(a,2));
    1903          49 :     a = gel(a,1); break;
    1904             :   default:
    1905          84 :     fast = get_oo(gen_0);
    1906             :   }
    1907         133 :   if (typ(a) != t_INT) pari_err_TYPE("sumnum", a);
    1908         133 :   if (!tab) tab = sumnuminit(fast, prec);
    1909          28 :   else if (!checksumtab(tab)) pari_err_TYPE("sumnum",tab);
    1910         126 :   as = itos(a);
    1911         126 :   d = gel(tab,1);
    1912         126 :   N = maxss(as, itos(gel(tab,2)));
    1913         126 :   k = itos(gel(tab,3));
    1914         126 :   v = gel(tab,4);
    1915         126 :   tabint = gel(tab,5);
    1916         126 :   prec2 = prec+EXTRAPREC;
    1917         126 :   av2 = avma;
    1918         126 :   S = gmul(eval(E, stoi(N)), real2n(-1,prec2));
    1919       15765 :   for (m = as; m < N; m++)
    1920             :   {
    1921       15639 :     S = gadd(S, eval(E, stoi(m)));
    1922       15639 :     if (gc_needed(av, 2))
    1923             :     {
    1924           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"sumnum [1], %ld/%ld",m,N);
    1925           0 :       S = gerepileupto(av2, S);
    1926             :     }
    1927       15639 :     S = gprec_wensure(S, prec2);
    1928             :   }
    1929        9611 :   for (m = 1; m <= k/2; m++)
    1930             :   {
    1931        9485 :     GEN t = gmulsg(2*m-1, d);
    1932        9485 :     GEN s = gsub(eval(E, gsubsg(N,t)), eval(E, gaddsg(N,t)));
    1933        9485 :     S = gadd(S, gmul(gel(v,m), s));
    1934        9485 :     if (gc_needed(av2, 2))
    1935             :     {
    1936           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"sumnum [2], %ld/%ld",m,k/2);
    1937           0 :       S = gerepileupto(av2, S);
    1938             :     }
    1939        9485 :     S = gprec_wensure(S, prec2);
    1940             :   }
    1941         126 :   S = gadd(S, intnum(E, eval,stoi(N), fast, tabint, prec2));
    1942         126 :   return gerepilecopy(av, gprec_wtrunc(S, prec));
    1943             : }
    1944             : 
    1945             : GEN
    1946         175 : sumnummonien0(GEN a, GEN code, GEN tab, long prec)
    1947         175 : { EXPR_WRAP(code, sumnummonien(EXPR_ARG, a, tab, prec)); }
    1948             : GEN
    1949          91 : sumnum0(GEN a, GEN code, GEN tab, long prec)
    1950          91 : { EXPR_WRAP(code, sumnum(EXPR_ARG, a, tab, prec)); }
    1951             : 
    1952             : /* Abel-Plana summation */
    1953             : 
    1954             : static GEN
    1955          49 : intnumgauexpinit(long prec)
    1956             : {
    1957          49 :   pari_sp ltop = avma;
    1958             :   GEN V, N, E, P, Q, R, vabs, vwt;
    1959          49 :   long l, n, k, j, prec2, prec0 = prec + EXTRAPREC, bit = prec2nbits(prec);
    1960             : 
    1961          49 :   n = (long)ceil(bit*0.226);
    1962          49 :   n |= 1; /* make n odd */
    1963          49 :   prec = nbits2prec(1.5*bit + 32);
    1964          49 :   prec2 = maxss(prec0, nbits2prec(1.15*bit + 32));
    1965          49 :   constbern(n+3);
    1966          49 :   V = cgetg(n + 4, t_VEC);
    1967        3045 :   for (k = 1; k <= n + 3; ++k)
    1968        2996 :     gel(V, k) = gtofp(gdivgs(bernfrac(2*k), odd(k)? 2*k: -2*k), prec);
    1969          49 :   Pade(V, &P, &Q);
    1970          49 :   N = RgX_recip(gsub(P, Q));
    1971          49 :   E = RgX_recip(Q);
    1972          49 :   R = gdivgs(gdiv(N, RgX_deriv(E)), 2);
    1973          49 :   vabs = RX_realroots(E,prec2);
    1974          49 :   l = lg(vabs); settyp(vabs, t_VEC);
    1975          49 :   vwt = cgetg(l, t_VEC);
    1976        1498 :   for (j = 1; j < l; ++j)
    1977             :   {
    1978        1449 :     GEN a = gel(vabs,j);
    1979        1449 :     gel(vwt, j) = gprec_wtrunc(poleval(R, a), prec0);
    1980        1449 :     gel(vabs, j) = gprec_wtrunc(sqrtr_abs(a), prec0);
    1981             :   }
    1982          49 :   return gerepilecopy(ltop, mkvec2(vabs, vwt));
    1983             : }
    1984             : 
    1985             : /* Compute \int_{-oo}^oo w(x)f(x) dx, where w(x)=x/(exp(2pi x)-1)
    1986             :  * for x>0 and w(-x)=w(x). For Abel-Plana (sumnumap). */
    1987             : static GEN
    1988          49 : intnumgauexp(void *E, GEN (*eval)(void*,GEN), GEN gN, GEN tab, long prec)
    1989             : {
    1990          49 :   pari_sp av = avma;
    1991          49 :   GEN U = mkcomplex(gN, NULL), V = mkcomplex(gN, NULL), S = gen_0;
    1992          49 :   GEN vabs = gel(tab, 1), vwt = gel(tab, 2);
    1993          49 :   long l = lg(vabs), i;
    1994          49 :   if (lg(vwt) != l || typ(vabs) != t_VEC || typ(vwt) != t_VEC)
    1995           0 :     pari_err_TYPE("sumnumap", tab);
    1996        1498 :   for (i = 1; i < l; i++)
    1997             :   { /* I * (w_i/a_i) * (f(N + I*a_i) - f(N - I*a_i)) */
    1998        1449 :     GEN x = gel(vabs,i), w = gel(vwt,i), t;
    1999        1449 :     gel(U,2) = x;
    2000        1449 :     gel(V,2) = gneg(x);
    2001        1449 :     t = mulcxI(gsub(eval(E,U), eval(E,V)));
    2002        1449 :     S = gadd(S, gmul(gdiv(w,x), cxtoreal(t)));
    2003        1449 :     S = gprec_wensure(S, prec);
    2004             :   }
    2005          49 :   return gerepilecopy(av, gprec_wtrunc(S, prec));
    2006             : }
    2007             : 
    2008             : GEN
    2009          49 : sumnumapinit(GEN fast, long prec)
    2010             : {
    2011          49 :   if (!fast) fast = mkoo();
    2012          49 :   retmkvec2(intnumgauexpinit(prec), intnuminit(gen_1, fast, 0, prec));
    2013             : }
    2014             : 
    2015             : typedef struct {
    2016             :   GEN (*f)(void *E, GEN);
    2017             :   void *E;
    2018             :   long N;
    2019             : } expfn;
    2020             : 
    2021             : /* f(Nx) */
    2022             : static GEN
    2023       31157 : _exfn(void *E, GEN x)
    2024             : {
    2025       31157 :   expfn *S = (expfn*)E;
    2026       31157 :   return S->f(S->E, gmulsg(S->N, x));
    2027             : }
    2028             : 
    2029             : GEN
    2030          56 : sumnumap(void *E, GEN (*eval)(void*,GEN), GEN a, GEN tab, long prec)
    2031             : {
    2032          56 :   pari_sp av = avma;
    2033             :   expfn T;
    2034             :   GEN S, fast, gN;
    2035             :   long as, m, N;
    2036          56 :   if (!a) { a = gen_1; fast = get_oo(gen_0); }
    2037          56 :   else switch(typ(a))
    2038             :   {
    2039             :     case t_VEC:
    2040          21 :       if (lg(a) != 3) pari_err_TYPE("sumnumap", a);
    2041          21 :       fast = get_oo(gel(a,2));
    2042          21 :       a = gel(a,1); break;
    2043             :     default:
    2044          35 :       fast = get_oo(gen_0);
    2045             :   }
    2046          56 :   if (typ(a) != t_INT) pari_err_TYPE("sumnumap", a);
    2047          56 :   if (!tab) tab = sumnumapinit(fast, prec);
    2048          14 :   else if (typ(tab) != t_VEC || lg(tab) != 3) pari_err_TYPE("sumnumap",tab);
    2049          49 :   as = itos(a);
    2050          49 :   T.N = N = maxss(as + 1, (long)ceil(prec2nbits(prec)*0.327));
    2051          49 :   T.E = E;
    2052          49 :   T.f = eval;
    2053          49 :   gN = stoi(N);
    2054          49 :   S = gtofp(gmul2n(eval(E, gN), -1), prec);
    2055        4109 :   for (m = as; m < N; ++m)
    2056             :   {
    2057        4060 :     S = gadd(S, eval(E, stoi(m)));
    2058        4060 :     S = gprec_wensure(S, prec);
    2059             :   }
    2060          49 :   S = gadd(S, gmulsg(N, intnum(&T, &_exfn, gen_1, fast, gel(tab, 2), prec)));
    2061          49 :   S = gadd(S, intnumgauexp(E, eval, gN, gel(tab, 1), prec));
    2062          49 :   return gerepileupto(av, S);
    2063             : }
    2064             : 
    2065             : GEN
    2066          56 : sumnumap0(GEN a, GEN code, GEN tab, long prec)
    2067          56 : { EXPR_WRAP(code, sumnumap(EXPR_ARG, a, tab, prec)); }
    2068             : 
    2069             : 
    2070             : /* max (1, |zeros|), P a t_POL or scalar */
    2071             : static double
    2072         133 : polmax(GEN P)
    2073             : {
    2074         133 :   pari_sp av = avma;
    2075             :   double r;
    2076         133 :   if (typ(P) != t_POL || degpol(P) <= 0) return 1.0;
    2077         133 :   r = gtodouble(polrootsbound(P, NULL));
    2078         133 :   return gc_double(av, maxdd(r, 1.0));
    2079             : }
    2080             : 
    2081             : /* max (1, |poles|), F a t_POL or t_RFRAC or scalar */
    2082             : static double
    2083          21 : ratpolemax(GEN F)
    2084             : {
    2085          21 :   if (typ(F) == t_POL) return 1.0;
    2086          21 :   return polmax(gel(F,2));
    2087             : }
    2088             : /* max (1, |poles|, |zeros|)) */
    2089             : static double
    2090          42 : ratpolemax2(GEN F)
    2091             : {
    2092          42 :   if (typ(F) == t_POL) return polmax(F);
    2093          42 :   return maxdd(polmax(gel(F,1)), polmax(gel(F,2)));
    2094             : }
    2095             : 
    2096             : static GEN
    2097       22099 : sercoeff(GEN x, long n)
    2098             : {
    2099       22099 :   long N = n - valp(x);
    2100       22099 :   return (N < 0)? gen_0: gel(x,N+2);
    2101             : }
    2102             : 
    2103             : /* Compute the integral from N to infinity of a rational function F, deg(F) < -1
    2104             :  * We must have N > 2 * r, r = max(1, |poles F|). */
    2105             : static GEN
    2106          28 : intnumainfrat(GEN F, long N, double r, long prec)
    2107             : {
    2108          28 :   long B = prec2nbits(prec), v, k, lim;
    2109             :   GEN S, ser;
    2110          28 :   pari_sp av = avma;
    2111             : 
    2112          28 :   lim = (long)ceil(B/log2(N/r));
    2113          28 :   ser = gmul(F, real_1(prec + EXTRAPREC));
    2114          28 :   ser = rfracrecip_to_ser_absolute(ser, lim+2);
    2115          28 :   v = valp(ser);
    2116          28 :   S = gdivgs(sercoeff(ser,lim+1), lim*N);
    2117             :   /* goes down to 2, but coeffs are 0 in degree < v */
    2118        1673 :   for (k = lim; k >= v; k--) /* S <- (S + coeff(ser,k)/(k-1)) / N */
    2119        1645 :     S = gdivgs(gadd(S, gdivgs(sercoeff(ser,k), k-1)), N);
    2120          28 :   if (v-2) S = gdiv(S, powuu(N, v-2));
    2121          28 :   return gerepilecopy(av, gprec_wtrunc(S, prec));
    2122             : }
    2123             : 
    2124             : static GEN
    2125          28 : rfrac_eval0(GEN R)
    2126             : {
    2127          28 :   GEN N, n, D = gel(R,2), d = constant_coeff(D);
    2128          28 :   if (gcmp0(d)) return NULL;
    2129          21 :   N = gel(R,1);
    2130          21 :   n = typ(N)==t_POL? constant_coeff(N): N;
    2131          21 :   return gdiv(n, d);
    2132             : }
    2133             : static GEN
    2134        2093 : rfrac_eval(GEN R, GEN a)
    2135             : {
    2136        2093 :   GEN D = gel(R,2), d = poleval(D,a);
    2137        2093 :   return gcmp0(d)? NULL: gdiv(poleval(gel(R,1),a), d);
    2138             : }
    2139             : /* R = \sum_i vR[i], eval at a omitting poles */
    2140             : static GEN
    2141        2093 : RFRAC_eval(GEN R, GEN vR, GEN a)
    2142             : {
    2143        2093 :   GEN S = rfrac_eval(R,a);
    2144        2093 :   if (!S && vR)
    2145             :   {
    2146           0 :     long i, l = lg(vR);
    2147           0 :     for (i = 1; i < l; i++)
    2148             :     {
    2149           0 :       GEN z = rfrac_eval(gel(vR,i), a);
    2150           0 :       if (z) S = S? gadd(S,z): z;
    2151             :     }
    2152             :   }
    2153        2093 :   return S;
    2154             : }
    2155             : static GEN
    2156        2093 : add_RFRAC_eval(GEN S, GEN R, GEN vR, GEN a)
    2157             : {
    2158        2093 :   GEN z = RFRAC_eval(R, vR, a);
    2159        2093 :   return z? gadd(S, z): S;
    2160             : }
    2161             : static GEN
    2162          21 : add_sumrfrac(GEN S, GEN R, GEN vR, long b)
    2163             : {
    2164             :   long m;
    2165          21 :   for (m = b; m >= 1; m--) S = add_RFRAC_eval(S,R,vR,utoipos(m));
    2166          21 :   return S;
    2167             : }
    2168             : static void
    2169          28 : get_kN(long r, long B, long *pk, long *pN)
    2170             : {
    2171          28 :   long k = maxss(50, (long)ceil(0.241*B));
    2172             :   GEN z;
    2173          28 :   if (k&1L) k++;
    2174          28 :   *pk = k; constbern(k >> 1);
    2175          28 :   z = sqrtnr_abs(gmul2n(gtofp(bernfrac(k), LOWDEFAULTPREC), B), k);
    2176          28 :   *pN = maxss(2*r, r + 1 + itos(gceil(z)));
    2177          28 : }
    2178             : /* F a t_RFRAC, F0 = F(0) or NULL [pole], vF a vector of t_RFRAC summing to F
    2179             :  * or NULL [F atomic] */
    2180             : static GEN
    2181          28 : sumnumrat_i(GEN F, GEN F0, GEN vF, long prec)
    2182             : {
    2183          28 :   long B = prec2nbits(prec), vx, j, k, N;
    2184             :   GEN S, S1, S2, intf, _1;
    2185             :   double r;
    2186          28 :   if (poldegree(F, -1) > -2) pari_err(e_MISC, "sum diverges in sumnumrat");
    2187          21 :   vx = varn(gel(F,2));
    2188          21 :   r = ratpolemax(F);
    2189          21 :   get_kN((long)ceil(r), B, &k,&N);
    2190          21 :   intf = intnumainfrat(F, N, r, prec);
    2191             :   /* N > ratpolemax(F) is not a pole */
    2192          21 :   _1 = real_1(prec);
    2193          21 :   S1 = gmul2n(gmul(_1, gsubst(F, vx, utoipos(N))), -1);
    2194          21 :   S1 = add_sumrfrac(S1, F, vF, N-1);
    2195          21 :   if (F0) S1 = gadd(S1, F0);
    2196          21 :   S = gmul(_1, gsubst(F, vx, gaddgs(pol_x(vx), N)));
    2197          21 :   S = rfrac_to_ser(S, k + 2);
    2198          21 :   S2 = gen_0;
    2199        1008 :   for (j = 2; j <= k; j += 2)
    2200         987 :     S2 = gadd(S2, gmul(gdivgs(bernfrac(j),j), sercoeff(S, j-1)));
    2201          21 :   return gadd(intf, gsub(S1, S2));
    2202             : }
    2203             : /* sum_{n >= a} F(n) */
    2204             : GEN
    2205          56 : sumnumrat(GEN F, GEN a, long prec)
    2206             : {
    2207          56 :   pari_sp av = avma;
    2208             :   long vx;
    2209             :   GEN vF, F0;
    2210             : 
    2211          56 :   switch(typ(F))
    2212             :   {
    2213          42 :     case t_RFRAC: break;
    2214             :     case t_INT: case t_REAL: case t_COMPLEX: case t_POL:
    2215          14 :       if (gequal0(F)) return real_0(prec);
    2216           7 :     default: pari_err_TYPE("sumnumrat",F);
    2217             :   }
    2218          42 :   vx = varn(gel(F,2));
    2219          42 :   switch(typ(a))
    2220             :   {
    2221             :     case t_INT:
    2222          21 :       if (signe(a)) F = gsubst(F, vx, deg1pol_shallow(gen_1,a,vx));
    2223          21 :       F0 = rfrac_eval0(F);
    2224          21 :       vF = NULL;
    2225          21 :       break;
    2226             :     case t_INFINITY:
    2227          21 :       if (inf_get_sign(a) == -1)
    2228             :       { /* F(x) + F(-x). Could divide degree by 2, as G(x^2): pb with poles */
    2229          14 :         GEN F2 = gsubst(F, vx, RgX_neg(pol_x(vx)));
    2230          14 :         vF = mkvec2(F,F2);
    2231          14 :         F = gadd(F, F2);
    2232          14 :         if (gequal0(F)) { set_avma(av); return real_0(prec); }
    2233           7 :         F0 = rfrac_eval0(gel(vF,1));
    2234           7 :         break;
    2235             :       }
    2236             :     default:
    2237           7 :       pari_err_TYPE("sumnumrat",a);
    2238             :       return NULL; /* LCOV_EXCL_LINE */
    2239             :   }
    2240          28 :   return gerepileupto(av, sumnumrat_i(F, F0, vF, prec));
    2241             : }
    2242             : /* deg ((a / b) - 1), assuming b a t_POL of positive degree in main variable  */
    2243             : static long
    2244          56 : rfracm1_degree(GEN a, GEN b)
    2245             : {
    2246             :   long da, db;
    2247          56 :   if (typ(a) != t_POL || varn(a) != varn(b)) return 0;
    2248          56 :   da = degpol(a);
    2249          56 :   db = degpol(b); if (da != db) return maxss(da - db, 0);
    2250          56 :   return degpol(RgX_sub(a,b)) - db;
    2251             : }
    2252             : 
    2253             : /* prod_{n >= a} F(n) */
    2254             : GEN
    2255          28 : prodnumrat(GEN F, long a, long prec)
    2256             : {
    2257          28 :   pari_sp ltop = avma;
    2258          28 :   long B = prec2nbits(prec), j, k, m, N, vx;
    2259             :   GEN S, S1, S2, intf, G;
    2260             :   double r;
    2261             : 
    2262          28 :   switch(typ(F))
    2263             :   {
    2264          14 :     case t_RFRAC: break;
    2265             :     case t_INT: case t_REAL: case t_COMPLEX: case t_POL:
    2266          14 :       if (gequal1(F)) return real_1(prec);
    2267           7 :     default: pari_err_TYPE("prodnumrat",F);
    2268             :   }
    2269          14 :   if (rfracm1_degree(gel(F,1), gel(F,2)) > -2)
    2270           7 :     pari_err(e_MISC, "product diverges in prodnumrat");
    2271           7 :   vx = varn(gel(F,2));
    2272           7 :   if (a) F = gsubst(F, vx, gaddgs(pol_x(vx), a));
    2273           7 :   r = ratpolemax2(F);
    2274           7 :   get_kN((long)ceil(r), B, &k,&N);
    2275           7 :   G = gdiv(deriv(F, vx), F);
    2276           7 :   intf = intnumainfrat(gmul(pol_x(vx),G), N, r, prec);
    2277           7 :   intf = gneg(gadd(intf, gmulsg(N, glog(gsubst(F, vx, stoi(N)), prec))));
    2278           7 :   S = gmul(real_1(prec), gsubst(G, vx, gaddgs(pol_x(vx), N)));
    2279           7 :   S = rfrac_to_ser(S, k + 2);
    2280           7 :   S1 = gsqrt(gsubst(F, vx, utoipos(N)), prec);
    2281           7 :   for (m = 0; m < N; m++) S1 = gmul(S1, gsubst(F, vx, utoi(m)));
    2282           7 :   S2 = gen_0;
    2283         336 :   for (j = 2; j <= k; j += 2)
    2284         329 :     S2 = gadd(S2, gmul(gdivgs(bernfrac(j),j*(j-1)), sercoeff(S, j-2)));
    2285           7 :   return gerepileupto(ltop, gmul(S1, gexp(gsub(intf, S2), prec)));
    2286             : }
    2287             : 
    2288             : /* fan = factoru(n); sum_{d | n} mu(d)/d * s[n/d] */
    2289             : static GEN
    2290        4956 : sdmob(GEN s, long n, GEN fan)
    2291             : {
    2292        4956 :   GEN D = divisorsu_moebius(gel(fan,1)), S = sercoeff(s, n); /* d = 1 */
    2293        4956 :   long i, l = lg(D);
    2294       19110 :   for (i = 2; i < l; i++)
    2295       14154 :     S = gadd(S, gdivgs(sercoeff(s, n/labs(D[i])), D[i]));
    2296        4956 :   return S;
    2297             : }
    2298             : /* log (zeta(s) * prod_i (1 - P[i]^-s) */
    2299             : static GEN
    2300        4151 : logzetan(GEN s, GEN P, long prec)
    2301             : {
    2302        4151 :   GEN Z = gzeta(s, prec);
    2303        4151 :   long i, l = lg(P);
    2304        4151 :   for (i = 1; i < l; i++) Z = gsub(Z, gdiv(Z, gpow(gel(P,i), s, prec)));
    2305        4151 :   return glog(Z, prec);
    2306             : }
    2307             : static GEN
    2308          56 : sumlogzeta(GEN ser, GEN s, GEN P, double rs, double lN, long vF, long lim,
    2309             :            long prec)
    2310             : {
    2311          56 :   GEN z = gen_0, v = vecfactoru_i(vF,lim);
    2312          56 :   pari_sp av = avma;
    2313             :   long i, n;
    2314          56 :   if (typ(s) == t_INT) constbern((itos(s) * lim + 1) >> 1);
    2315        5012 :   for (n = lim, i = lg(v)-1; n >= vF; n--, i--)
    2316             :   {
    2317        4956 :     GEN t = sdmob(ser, n, gel(v,i));
    2318        4956 :     if (!gequal0(t))
    2319             :     { /* (n Re(s) - 1) log2(N) bits cancel in logzetan */
    2320        4151 :       long prec2 = prec + nbits2extraprec((n*rs-1) * lN);
    2321        4151 :       GEN L = logzetan(gmulsg(n,gprec_wensure(s,prec2)), P, prec2);
    2322        4151 :       z = gerepileupto(av, gadd(z, gmul(L, t)));
    2323        4151 :       z = gprec_wensure(z, prec);
    2324             :     }
    2325             :   }
    2326          56 :   return gprec_wtrunc(z, prec);
    2327             : }
    2328             : 
    2329             : static GEN
    2330         595 : rfrac_evalfp(GEN F, GEN x, long prec)
    2331             : {
    2332         595 :   GEN N = gel(F,1), D = gel(F,2), a, b = poleval(D, x);
    2333         595 :   a = (typ(N) == t_POL && varn(N) == varn(D))? poleval(N, x): N;
    2334        1120 :   if (typ(a) != t_INT || typ(b) != t_INT ||
    2335         980 :       (lgefint(a) <= prec && lgefint(b) <= prec)) return gdiv(a, b);
    2336          70 :   return rdivii(a, b, prec + EXTRAPRECWORD);
    2337             : }
    2338             : 
    2339             : /* { F(p^s): p in P, p >= a }, F t_RFRAC */
    2340             : static GEN
    2341          56 : vFps(GEN P, long a, GEN F, GEN s, long prec)
    2342             : {
    2343          56 :   long i, j, l = lg(P);
    2344          56 :   GEN v = cgetg(l, t_VEC);
    2345         658 :   for (i = j = 1; i < l; i++)
    2346             :   {
    2347         602 :     GEN p = gel(P,i); if (cmpiu(p, a) < 0) continue;
    2348         595 :     gel(v, j++) = rfrac_evalfp(F, gpow(p, s, prec), prec);
    2349             :   }
    2350          56 :   setlg(v, j); return v;
    2351             : }
    2352             : 
    2353             : static void
    2354          98 : euler_set_Fs(GEN *F, GEN *s)
    2355             : {
    2356          98 :   if (!*s) *s = gen_1;
    2357          98 :   if (typ(*F) == t_RFRAC)
    2358             :   {
    2359             :     long m;
    2360          70 :     *F = rfrac_deflate_max(*F, &m);
    2361          70 :     if (m != 1) *s = gmulgs(*s, m);
    2362             :   }
    2363          98 : }
    2364             : /* sum_{p prime, p >= a} F(p^s), F rational function */
    2365             : GEN
    2366          42 : sumeulerrat(GEN F, GEN s, long a, long prec)
    2367             : {
    2368          42 :   pari_sp av = avma;
    2369             :   GEN ser, z, P;
    2370             :   double r, rs, RS, lN;
    2371          42 :   long B = prec2nbits(prec), prec2 = prec + EXTRAPREC, vF, N, lim;
    2372             : 
    2373          42 :   euler_set_Fs(&F, &s);
    2374          42 :   switch(typ(F))
    2375             :   {
    2376          28 :     case t_RFRAC: break;
    2377             :     case t_INT: case t_REAL: case t_COMPLEX: case t_POL:
    2378          14 :       if (gequal0(F)) return real_0(prec);
    2379           7 :     default: pari_err_TYPE("sumeulerrat",F);
    2380             :   }
    2381             :   /* F t_RFRAC */
    2382          28 :   if (a < 2) a = 2;
    2383          28 :   vF = -poldegree(F, -1);
    2384          28 :   rs = gtodouble(real_i(s));
    2385          28 :   r = polmax(gel(F,2));
    2386          28 :   N = maxss(30, a); lN = log2((double)N);
    2387          28 :   RS = maxdd(1./vF, log2(r) / lN);
    2388          28 :   if (rs <= RS)
    2389           7 :     pari_err_DOMAIN("sumeulerrat", "real(s)", "<=",  dbltor(RS), dbltor(rs));
    2390          21 :   lim = (long)ceil(B / (rs*lN - log2(r)));
    2391          21 :   ser = gmul(real_1(prec2), F);
    2392          21 :   ser = rfracrecip_to_ser_absolute(ser, lim+1);
    2393          21 :   P = primes_interval(gen_2, utoipos(N));
    2394          21 :   z = sumlogzeta(ser, s, P, rs, lN, vF, lim, prec);
    2395          21 :   z = gadd(z, vecsum(vFps(P, a, F, s, prec)));
    2396          21 :   return gerepilecopy(av, gprec_wtrunc(z, prec));
    2397             : }
    2398             : 
    2399             : /* F = N/D; return F'/F. Assume D a t_POL */
    2400             : static GEN
    2401          35 : rfrac_logderiv(GEN N, GEN D)
    2402             : {
    2403          35 :   if (typ(N) != t_POL || varn(N) != varn(D)) return gdiv(gneg(RgX_deriv(D)), D);
    2404          35 :   if (!degpol(D)) return gdiv(RgX_deriv(N), N);
    2405          14 :   return gdiv(RgX_sub(RgX_mul(RgX_deriv(N), D), RgX_mul(RgX_deriv(D), N)),
    2406             :               RgX_mul(N, D));
    2407             : }
    2408             : 
    2409             : /* prod_{p prime, p >= a} F(p^s), F rational function */
    2410             : GEN
    2411          56 : prodeulerrat(GEN F, GEN s, long a, long prec)
    2412             : {
    2413          56 :   pari_sp ltop = avma;
    2414             :   GEN DF, NF, ser, P, z;
    2415             :   double r, rs, RS, lN;
    2416          56 :   long B = prec2nbits(prec), prec2 = prec + EXTRAPREC, vF, N, lim;
    2417             : 
    2418          56 :   euler_set_Fs(&F, &s);
    2419          56 :   switch(typ(F))
    2420             :   {
    2421          42 :     case t_RFRAC: break;
    2422             :     case t_INT: case t_REAL: case t_COMPLEX: case t_POL:
    2423          14 :       if (gequal1(F)) return real_1(prec);
    2424           7 :     default: pari_err_TYPE("prodeulerrat",F);
    2425             :   } /* F t_RFRAC */
    2426          42 :   NF = gel(F, 1);
    2427          42 :   DF = gel(F, 2);
    2428          42 :   rs = gtodouble(real_i(s));
    2429          42 :   vF = - rfracm1_degree(NF, DF);
    2430          42 :   if (rs * vF <= 1) pari_err(e_MISC, "product diverges in prodeulerrat");
    2431          35 :   r = ratpolemax2(F);
    2432          35 :   N = maxss(maxss(30, a), (long)ceil(2*r)); lN = log2((double)N);
    2433          35 :   RS = maxdd(1./vF, log2(r) / lN);
    2434          35 :   if (rs <= RS)
    2435           0 :     pari_err_DOMAIN("prodeulerrat", "real(s)", "<=",  dbltor(RS), dbltor(rs));
    2436          35 :   lim = (long)ceil(B / (rs*lN - log2(r)));
    2437          35 :   (void)rfracrecip(&NF, &DF); /* returned value is 0 */
    2438          35 :   if (!RgX_is_ZX(DF) || !is_pm1(gel(DF,2))
    2439          35 :       || lim * log2(r) > 4 * B) NF = gmul(NF, real_1(prec2));
    2440          35 :   ser = integser(rfrac_to_ser(rfrac_logderiv(NF,DF), lim+3));
    2441             :   /* ser = log f, f = F(1/x) + O(x^(lim+1)) */
    2442          35 :   P = primes_interval(gen_2, utoipos(N));
    2443          35 :   z = gexp(sumlogzeta(ser, s, P, rs, lN, vF, lim, prec), prec);
    2444          35 :   z = gmul(z, vecprod(vFps(P, a, F, s, prec)));
    2445          35 :   return gerepilecopy(ltop, gprec_wtrunc(z, prec));
    2446             : }
    2447             : 
    2448             : /* Compute $\sum_{n\ge a}c(n)$ using Lagrange extrapolation.
    2449             : Assume that the $N$-th remainder term of the series has a
    2450             : regular asymptotic expansion in integral powers of $1/N$. */
    2451             : static GEN
    2452          35 : sumnumlagrange1init(GEN c1, long flag, long prec)
    2453             : {
    2454          35 :   pari_sp av = avma;
    2455             :   GEN V, W, T;
    2456             :   double c1d;
    2457          35 :   long B = prec2nbits(prec), prec2;
    2458             :   ulong n, N;
    2459          35 :   c1d = c1 ? gtodouble(c1) : 0.332;
    2460          35 :   N = (ulong)ceil(c1d*B); if ((N&1L) == 0) N++;
    2461          35 :   prec2 = nbits2prec(B+(long)ceil(1.8444*N) + 16);
    2462          35 :   W = vecbinome(N);
    2463          35 :   T = vecpowuu(N, N);
    2464          35 :   V = cgetg(N+1, t_VEC); gel(V,N) = gel(T,N);
    2465        3773 :   for (n = N-1; n >= 1; n--)
    2466             :   {
    2467        3738 :     pari_sp av = avma;
    2468        3738 :     GEN t = mulii(gel(W, n+1), gel(T,n));
    2469        3738 :     if (!odd(n)) togglesign_safe(&t);
    2470        3738 :     if (flag) t = addii(gel(V, n+1), t);
    2471        3738 :     gel(V, n) = gerepileuptoint(av, t);
    2472             :   }
    2473          35 :   V = gdiv(RgV_gtofp(V, prec2), mpfact(N));
    2474          35 :   return gerepilecopy(av, mkvec4(gen_1, stoi(prec2), gen_1, V));
    2475             : }
    2476             : 
    2477             : static GEN
    2478           7 : sumnumlagrange2init(GEN c1, long flag, long prec)
    2479             : {
    2480           7 :   pari_sp av = avma;
    2481             :   GEN V, W, T, told;
    2482           7 :   double c1d = c1 ? gtodouble(c1) : 0.228;
    2483           7 :   long B = prec2nbits(prec), prec2;
    2484             :   ulong n, N;
    2485             : 
    2486           7 :   N = (ulong)ceil(c1d*B); if ((N&1L) == 0) N++;
    2487           7 :   prec2 = nbits2prec(B+(long)ceil(1.18696*N) + 16);
    2488           7 :   W = vecbinome(2*N);
    2489           7 :   T = vecpowuu(N, 2*N);
    2490           7 :   V = cgetg(N+1, t_VEC); gel(V, N) = told = gel(T,N);
    2491         623 :   for (n = N-1; n >= 1; n--)
    2492             :   {
    2493         616 :     GEN tnew = mulii(gel(W, N-n+1), gel(T,n));
    2494         616 :     if (!odd(n)) togglesign_safe(&tnew);
    2495         616 :     told = addii(told, tnew);
    2496         616 :     if (flag) told = addii(gel(V, n+1), told);
    2497         616 :     gel(V, n) = told; told = tnew;
    2498             :   }
    2499           7 :   V = gdiv(RgV_gtofp(V, prec2), mpfact(2*N));
    2500           7 :   return gerepilecopy(av, mkvec4(gen_2, stoi(prec2), gen_1, V));
    2501             : }
    2502             : 
    2503             : static GEN
    2504     1844164 : _mpmul(GEN x, GEN y)
    2505             : {
    2506     1844164 :   if (!x) return y;
    2507     1839306 :   return y? mpmul(x, y): x;
    2508             : }
    2509             : /* Used only for al = 2, 1, 1/2, 1/3, 1/4. */
    2510             : static GEN
    2511          49 : sumnumlagrangeinit_i(GEN al, GEN c1, long flag, long prec)
    2512             : {
    2513          49 :   pari_sp av = avma;
    2514             :   GEN V, W;
    2515          49 :   double c1d = 0.0, c2;
    2516          49 :   long B = prec2nbits(prec), B1, prec2, dal;
    2517             :   ulong j, n, N;
    2518             : 
    2519          49 :   if (typ(al) == t_INT)
    2520             :   {
    2521          28 :     switch(itos_or_0(al))
    2522             :     {
    2523          21 :       case 1: return sumnumlagrange1init(c1, flag, prec);
    2524           7 :       case 2: return sumnumlagrange2init(c1, flag, prec);
    2525           0 :       default: pari_err_IMPL("sumnumlagrange for this alpha");
    2526             :     }
    2527             :   }
    2528          21 :   if (typ(al) != t_FRAC) pari_err_TYPE("sumnumlagrangeinit",al);
    2529          21 :   dal = itos_or_0(gel(al,2));
    2530          21 :   if (dal > 4 || !equali1(gel(al,1)))
    2531           7 :     pari_err_IMPL("sumnumlagrange for this alpha");
    2532          14 :   switch(dal)
    2533             :   {
    2534           7 :     case 2: c2 = 2.6441; c1d = 0.62; break;
    2535           7 :     case 3: c2 = 3.1578; c1d = 1.18; break;
    2536           0 :     case 4: c2 = 3.5364; c1d = 3.00; break;
    2537             :     default: return NULL; /* LCOV_EXCL_LINE */
    2538             :   }
    2539          14 :   if (c1)
    2540             :   {
    2541           0 :     c1d = gtodouble(c1);
    2542           0 :     if (c1d <= 0)
    2543           0 :       pari_err_DOMAIN("sumnumlagrangeinit", "c1", "<=", gen_0, c1);
    2544             :   }
    2545          14 :   N = (ulong)ceil(c1d*B); if ((N&1L) == 0) N++;
    2546          14 :   B1 = B + (long)ceil(c2*N) + 16;
    2547          14 :   prec2 = nbits2prec(B1);
    2548          14 :   V = vecpowug(N, al, prec2);
    2549          14 :   W = cgetg(N+1, t_VEC);
    2550        4872 :   for (n = 1; n <= N; ++n)
    2551             :   {
    2552        4858 :     pari_sp av2 = avma;
    2553        4858 :     GEN t = NULL, vn = gel(V, n);
    2554     1853880 :     for (j = 1; j <= N; j++)
    2555     1849022 :       if (j != n) t = _mpmul(t, mpsub(vn, gel(V, j)));
    2556        4858 :     gel(W, n) = gerepileuptoleaf(av2, mpdiv(gpowgs(vn, N-1), t));
    2557             :   }
    2558          14 :   if (flag)
    2559          14 :     for (n = N-1; n >= 1; n--) gel(W, n) = gadd(gel(W, n+1), gel(W, n));
    2560          14 :   return gerepilecopy(av, mkvec4(al, stoi(prec2), gen_1, W));
    2561             : }
    2562             : 
    2563             : GEN
    2564          63 : sumnumlagrangeinit(GEN al, GEN c1, long prec)
    2565             : {
    2566          63 :   pari_sp ltop = avma;
    2567             :   GEN V, W, S, be;
    2568             :   long n, prec2, fl, N;
    2569             : 
    2570          63 :   if (!al) return sumnumlagrange1init(c1, 1, prec);
    2571          49 :   if (typ(al) != t_VEC) al = mkvec2(gen_1, al);
    2572          35 :   else if (lg(al) != 3) pari_err_TYPE("sumnumlagrangeinit",al);
    2573          49 :   be = gel(al,2);
    2574          49 :   al = gel(al,1);
    2575          49 :   if (gequal0(be)) return sumnumlagrangeinit_i(al, c1, 1, prec);
    2576          14 :   V = sumnumlagrangeinit_i(al, c1, 0, prec);
    2577          14 :   switch(typ(be))
    2578             :   {
    2579           0 :     case t_CLOSURE: fl = 1; break;
    2580           7 :     case t_INT: case t_FRAC: case t_REAL: fl = 0; break;
    2581           7 :     default: pari_err_TYPE("sumnumlagrangeinit", be);
    2582             :              return NULL; /* LCOV_EXCL_LINE */
    2583             :   }
    2584           7 :   prec2 = itos(gel(V, 2));
    2585           7 :   W = gel(V, 4);
    2586           7 :   N = lg(W) - 1;
    2587           7 :   S = gen_0; V = cgetg(N+1, t_VEC);
    2588         910 :   for (n = N; n >= 1; n--)
    2589             :   {
    2590         903 :     GEN tmp, gn = stoi(n);
    2591         903 :     tmp = fl ? closure_callgen1prec(be, gn, prec2) : gpow(gn, gneg(be), prec2);
    2592         903 :     tmp = gdiv(gel(W, n), tmp);
    2593         903 :     S = gadd(S, tmp);
    2594         903 :     gel(V, n) = (n == N)? tmp: gadd(gel(V, n+1), tmp);
    2595             :   }
    2596           7 :   return gerepilecopy(ltop, mkvec4(al, stoi(prec2), S, V));
    2597             : }
    2598             : 
    2599             : /* - sum_{n=1}^{as-1} f(n) */
    2600             : static GEN
    2601          14 : sumaux(void *E, GEN (*eval)(void*,GEN,long), long as, long prec)
    2602             : {
    2603          14 :   GEN S = gen_0;
    2604             :   long n;
    2605          14 :   if (as > 1)
    2606             :   {
    2607          14 :     for (n = 1; n < as; ++n)
    2608             :     {
    2609           7 :       S = gadd(S, eval(E, stoi(n), prec));
    2610           7 :       S = gprec_wensure(S, prec);
    2611             :     }
    2612           7 :     S = gneg(S);
    2613             :   }
    2614             :   else
    2615           7 :     for (n = as; n <= 0; ++n)
    2616             :     {
    2617           0 :       S = gadd(S, eval(E, stoi(n), prec));
    2618           0 :       S = gprec_wensure(S, prec);
    2619             :     }
    2620          14 :   return S;
    2621             : }
    2622             : 
    2623             : GEN
    2624          84 : sumnumlagrange(void *E, GEN (*eval)(void*,GEN,long), GEN a, GEN tab, long prec)
    2625             : {
    2626          84 :   pari_sp av = avma;
    2627             :   GEN s, S, al, V;
    2628             :   long as, prec2;
    2629             :   ulong n, l;
    2630             : 
    2631          84 :   if (typ(a) != t_INT) pari_err_TYPE("sumnumlagrange", a);
    2632          84 :   if (!tab) tab = sumnumlagrangeinit(NULL, tab, prec);
    2633          70 :   else if (lg(tab) != 5 || typ(gel(tab,2)) != t_INT || typ(gel(tab,4)) != t_VEC)
    2634           0 :     pari_err_TYPE("sumnumlagrange", tab);
    2635             : 
    2636          84 :   as = itos(a);
    2637          84 :   al = gel(tab, 1);
    2638          84 :   prec2 = itos(gel(tab, 2));
    2639          84 :   S = gel(tab, 3);
    2640          84 :   V = gel(tab, 4);
    2641          84 :   l = lg(V);
    2642          84 :   if (gequal(al, gen_2))
    2643             :   {
    2644          14 :     s = sumaux(E, eval, as, prec2);
    2645          14 :     as = 1;
    2646             :   }
    2647             :   else
    2648          70 :     s = gen_0;
    2649       16464 :   for (n = 1; n < l; n++)
    2650             :   {
    2651       16380 :     s = gadd(s, gmul(gel(V, n), eval(E, stoi(n+as-1), prec2)));
    2652       16380 :     s = gprec_wensure(s, prec);
    2653             :   }
    2654          84 :   if (!gequal1(S)) s = gdiv(s,S);
    2655          84 :   return gerepilecopy(av, gprec_wtrunc(s, prec));
    2656             : }
    2657             : 
    2658             : GEN
    2659          84 : sumnumlagrange0(GEN a, GEN code, GEN tab, long prec)
    2660          84 : { EXPR_WRAP(code, sumnumlagrange(EXPR_ARGPREC, a, tab, prec)); }

Generated by: LCOV version 1.13