Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - language - intnum.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 17837-9309d7c) Lines: 972 1021 95.2 %
Date: 2015-05-22 Functions: 77 79 97.5 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 479 634 75.6 %

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

Generated by: LCOV version 1.9