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

Generated by: LCOV version 1.14