Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - language - intnum.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 19588-1c9967d) Lines: 952 1007 94.5 %
Date: 2016-09-24 05:54:29 Functions: 78 81 96.3 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.11