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 17242-6994acc) Lines: 850 947 89.8 %
Date: 2014-12-22 Functions: 87 92 94.6 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 396 624 63.5 %

           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                 :            : /**                                                                **/
      19                 :            : /**                NUMERICAL INTEGRATION (Romberg)                 **/
      20                 :            : /**                                                                **/
      21                 :            : /********************************************************************/
      22                 :            : typedef struct {
      23                 :            :   void *E;
      24                 :            :   GEN (*f)(void *E, GEN);
      25                 :            : } invfun;
      26                 :            : 
      27                 :            : /* 1/x^2 f(1/x) */
      28                 :            : static GEN
      29                 :     413343 : _invf(void *E, GEN x)
      30                 :            : {
      31                 :     413343 :   invfun *S = (invfun*)E;
      32                 :     413343 :   GEN y = ginv(x);
      33                 :     413343 :   return gmul(S->f(S->E, y), gsqr(y));
      34                 :            : }
      35                 :            : 
      36                 :            : static GEN
      37                 :        203 : interp(GEN h, GEN s, long j, long lim, long KLOC)
      38                 :            : {
      39                 :        203 :   pari_sp av = avma;
      40                 :            :   long e1,e2;
      41                 :        203 :   GEN dss, ss = polint_i(h+j-KLOC,s+j-KLOC,gen_0,KLOC+1,&dss);
      42                 :            : 
      43                 :        203 :   e1 = gexpo(ss);
      44                 :        203 :   e2 = gexpo(dss);
      45 [ +  + ][ +  + ]:        203 :   if (e1-e2 <= lim && (j <= 10 || e1 >= -lim)) { avma = av; return NULL; }
                 [ +  - ]
      46 [ -  + ][ #  # ]:         28 :   if (typ(ss) == t_COMPLEX && gequal0(gel(ss,2))) ss = gel(ss,1);
      47                 :        203 :   return ss;
      48                 :            : }
      49                 :            : 
      50                 :            : static GEN
      51                 :          7 : qrom3(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, long prec)
      52                 :            : {
      53                 :          7 :   const long JMAX = 25, KLOC = 4;
      54                 :            :   GEN ss,s,h,p1,p2,qlint,del,x,sum;
      55                 :            :   long j, j1, it, sig;
      56                 :            : 
      57                 :          7 :   a = gtofp(a,prec);
      58                 :          7 :   b = gtofp(b,prec);
      59                 :          7 :   qlint = subrr(b,a); sig = signe(qlint);
      60         [ -  + ]:          7 :   if (!sig)  return gen_0;
      61         [ -  + ]:          7 :   if (sig < 0) { setabssign(qlint); swap(a,b); }
      62                 :            : 
      63                 :          7 :   s = new_chunk(JMAX+KLOC-1);
      64                 :          7 :   h = new_chunk(JMAX+KLOC-1);
      65                 :          7 :   gel(h,0) = real_1(prec);
      66                 :            : 
      67         [ -  + ]:          7 :   p1 = eval(E, a); if (p1 == a) p1 = rcopy(p1);
      68                 :          7 :   p2 = eval(E, b);
      69                 :          7 :   gel(s,0) = gmul2n(gmul(qlint,gadd(p1,p2)),-1);
      70         [ +  - ]:         91 :   for (it=1,j=1; j<JMAX; j++, it<<=1) /* it = 2^(j-1) */
      71                 :            :   {
      72                 :            :     pari_sp av, av2;
      73                 :         91 :     gel(h,j) = real2n(-2*j, prec); /* 2^(-2j) */
      74                 :         91 :     av = avma; del = divru(qlint,it);
      75                 :         91 :     x = addrr(a, shiftr(del,-1));
      76                 :         91 :     av2 = avma;
      77         [ +  + ]:      57428 :     for (sum = gen_0, j1 = 1; j1 <= it; j1++, x = addrr(x,del))
      78                 :            :     {
      79                 :      57337 :       sum = gadd(sum, eval(E, x));
      80         [ +  + ]:      57337 :       if ((j1 & 0x1ff) == 0) gerepileall(av2, 2, &sum,&x);
      81                 :            :     }
      82                 :         91 :     sum = gmul(sum,del);
      83                 :         91 :     gel(s,j) = gerepileupto(av, gmul2n(gadd(gel(s,j-1), sum), -1));
      84         [ -  + ]:         91 :     if (DEBUGLEVEL>3) err_printf("qrom3: iteration %ld: %Ps\n", j,s[j]);
      85                 :            : 
      86 [ +  + ][ +  + ]:         91 :     if (j >= KLOC && (ss = interp(h, s, j, prec2nbits(prec)-j-6, KLOC)))
      87                 :          7 :       return gmulsg(sig,ss);
      88                 :            :   }
      89                 :          7 :   return NULL;
      90                 :            : }
      91                 :            : 
      92                 :            : static GEN
      93                 :         21 : qrom2(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, long prec)
      94                 :            : {
      95                 :         21 :   const long JMAX = 16, KLOC = 4;
      96                 :            :   GEN ss,s,h,p1,qlint,del,ddel,x,sum;
      97                 :            :   long j, j1, it, sig;
      98                 :            : 
      99                 :         21 :   a = gtofp(a, prec);
     100                 :         21 :   b = gtofp(b, prec);
     101                 :         21 :   qlint = subrr(b,a); sig = signe(qlint);
     102         [ -  + ]:         21 :   if (!sig)  return gen_0;
     103         [ -  + ]:         21 :   if (sig < 0) { setabssign(qlint); swap(a,b); }
     104                 :            : 
     105                 :         21 :   s = new_chunk(JMAX+KLOC-1);
     106                 :         21 :   h = new_chunk(JMAX+KLOC-1);
     107                 :         21 :   gel(h,0) = real_1(prec);
     108                 :            : 
     109                 :         21 :   p1 = shiftr(addrr(a,b),-1);
     110                 :         21 :   gel(s,0) = gmul(qlint, eval(E, p1));
     111         [ +  - ]:        196 :   for (it=1, j=1; j<JMAX; j++, it*=3) /* it = 3^(j-1) */
     112                 :            :   {
     113                 :            :     pari_sp av, av2;
     114                 :        196 :     gel(h,j) = divru(gel(h,j-1), 9); /* 3^(-2j) */
     115                 :        196 :     av = avma; del = divru(qlint,3*it); ddel = shiftr(del,1);
     116                 :        196 :     x = addrr(a, shiftr(del,-1));
     117                 :        196 :     av2 = avma;
     118         [ +  + ]:     344638 :     for (sum = gen_0, j1 = 1; j1 <= it; j1++)
     119                 :            :     {
     120                 :     344442 :       sum = gadd(sum, eval(E, x)); x = addrr(x,ddel);
     121                 :     344442 :       sum = gadd(sum, eval(E, x)); x = addrr(x,del);
     122         [ +  + ]:     344442 :       if ((j1 & 0x1ff) == 0) gerepileall(av2, 2, &sum,&x);
     123                 :            :     }
     124                 :        196 :     sum = gmul(sum,del); p1 = gdivgs(gel(s,j-1),3);
     125                 :        196 :     gel(s,j) = gerepileupto(av, gadd(p1,sum));
     126         [ -  + ]:        196 :     if (DEBUGLEVEL>3) err_printf("qrom2: iteration %ld: %Ps\n", j,s[j]);
     127                 :            : 
     128 [ +  + ][ +  + ]:        196 :     if (j >= KLOC && (ss = interp(h, s, j, prec2nbits(prec)-(3*j/2)-6, KLOC)))
     129                 :         21 :       return gmulsg(sig, ss);
     130                 :            :   }
     131                 :         21 :   return NULL;
     132                 :            : }
     133                 :            : 
     134                 :            : /* integrate after change of variables x --> 1/x */
     135                 :            : static GEN
     136                 :          7 : qromi(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, long prec)
     137                 :            : {
     138                 :          7 :   GEN A = ginv(b), B = ginv(a);
     139                 :            :   invfun S;
     140                 :          7 :   S.f = eval;
     141                 :          7 :   S.E = E; return qrom2(&S, &_invf, A, B, prec);
     142                 :            : }
     143                 :            : 
     144                 :            : /* a < b, assume b "small" (< 100 say) */
     145                 :            : static GEN
     146                 :          7 : rom_bsmall(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, long prec)
     147                 :            : {
     148         [ +  - ]:          7 :   if (gcmpgs(a,-100) >= 0) return qrom2(E,eval,a,b,prec);
     149 [ #  # ][ #  # ]:          0 :   if (b == gen_1 || gcmpgs(b, -1) >= 0) /* a < -100, b >= -1 */
     150                 :          0 :     return gadd(qromi(E,eval,a,gen_m1,prec), /* split at -1 */
     151                 :            :                 qrom2(E,eval,gen_m1,b,prec));
     152                 :            :   /* a < -100, b < -1 */
     153                 :          7 :   return qromi(E,eval,a,b,prec);
     154                 :            : }
     155                 :            : 
     156                 :            : static GEN
     157                 :          7 : rombint(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, long prec)
     158                 :            : {
     159                 :          7 :   long l = gcmp(b,a);
     160                 :            :   GEN z;
     161                 :            : 
     162         [ -  + ]:          7 :   if (!l) return gen_0;
     163         [ -  + ]:          7 :   if (l < 0) swap(a,b);
     164         [ -  + ]:          7 :   if (gcmpgs(b,100) >= 0)
     165                 :            :   {
     166         [ #  # ]:          0 :     if (gcmpgs(a,1) >= 0)
     167                 :          0 :       z = qromi(E,eval,a,b,prec);
     168                 :            :     else /* split at 1 */
     169                 :          0 :       z = gadd(rom_bsmall(E,eval,a,gen_1,prec), qromi(E,eval,gen_1,b,prec));
     170                 :            :   }
     171                 :            :   else
     172                 :          7 :     z = rom_bsmall(E,eval,a,b,prec);
     173         [ -  + ]:          7 :   if (l < 0) z = gneg(z);
     174                 :          7 :   return z;
     175                 :            : }
     176                 :            : 
     177                 :            : /********************************************************************/
     178                 :            : /**                                                                **/
     179                 :            : /**                DOUBLE EXPONENTIAL INTEGRATION                  **/
     180                 :            : /**                                                                **/
     181                 :            : /********************************************************************/
     182                 :            : 
     183                 :            : /* The init functions have the following purposes:
     184                 :            : * 1) They fill the value tabx0 = phi(0) and arrays of abcissas
     185                 :            : *   tabxp[] = phi(k/2^m) (positive x) and also of tabxm[] = phi(-k/2^m)
     186                 :            : *   (negative x) unless the phi function is odd, in which case this is useless.
     187                 :            : * 2) They fill the corresponding arrays of weights tabw0 = phi'(0) and
     188                 :            : *   tabwp[] = phi'(k/2^m) (and possibly also of tabwm[] = phi'(-k/2^m)).
     189                 :            : * 3) They set eps to the desired accuracy (depending on the GP default).
     190                 :            : * 4) They compute nt which says that the weights tabwp[k] and tabwm[k] are
     191                 :            : *   negligible with respect to eps if k > nt. In particular the tabxx[] arrays
     192                 :            : *   are indexed from 1 to nt+1. */
     193                 :            : 
     194                 :            : typedef struct _intdata {
     195                 :            :   long m;    /* integration step h = 1/2^m */
     196                 :            :   long eps;  /* bit accuracy of current precision */
     197                 :            :   GEN tabx0; /* abcissa phi(0) for t = 0 */
     198                 :            :   GEN tabw0; /* weight phi'(0) for t = 0 */
     199                 :            :   GEN tabxp; /* table of abcissas phi(kh) for k > 0 */
     200                 :            :   GEN tabwp; /* table of weights phi'(kh) for k > 0 */
     201                 :            :   GEN tabxm; /* table of abcissas phi(kh) for k < 0 */
     202                 :            :   GEN tabwm; /* table of weights phi'(kh) for k < 0 */
     203                 :            : } intdata;
     204                 :            : 
     205                 :            : #define TABm(v)  gel(v,1)
     206                 :            : #define TABx0(v) gel(v,2)
     207                 :            : #define TABw0(v) gel(v,3)
     208                 :            : #define TABxp(v) gel(v,4)
     209                 :            : #define TABwp(v) gel(v,5)
     210                 :            : #define TABxm(v) gel(v,6)
     211                 :            : #define TABwm(v) gel(v,7)
     212                 :            : 
     213                 :            : static int
     214                 :      19740 : isinR(GEN z)
     215                 :            : {
     216                 :      19740 :   long tz = typ(z);
     217 [ +  + ][ +  + ]:      19740 :   return (tz == t_INT || tz == t_REAL || tz == t_FRAC);
                 [ +  - ]
     218                 :            : }
     219                 :            : 
     220                 :            : static int
     221                 :      19257 : isinC(GEN z)
     222                 :            : {
     223 [ +  + ][ +  - ]:      19257 :   return (typ(z) == t_COMPLEX)? isinR(gel(z,1)) && isinR(gel(z,2)):
                 [ +  - ]
     224                 :            :                                 isinR(z);
     225                 :            : }
     226                 :            : 
     227                 :            : static int
     228                 :      18809 : checktabsimp(GEN tab)
     229                 :            : {
     230                 :            :   long L, LN, LW;
     231 [ +  - ][ -  + ]:      18809 :   if (!tab || typ(tab) != t_VEC) return 0;
     232         [ -  + ]:      18809 :   if (lg(tab) != 8) return 0;
     233         [ -  + ]:      18809 :   if (typ(TABm(tab))!= t_INT) return 0;
     234         [ -  + ]:      18809 :   if (typ(TABxp(tab)) != t_VEC) return 0;
     235         [ -  + ]:      18809 :   if (typ(TABwp(tab)) != t_VEC) return 0;
     236         [ -  + ]:      18809 :   if (typ(TABxm(tab)) != t_VEC) return 0;
     237         [ -  + ]:      18809 :   if (typ(TABwm(tab)) != t_VEC) return 0;
     238         [ -  + ]:      18809 :   L = lg(TABxp(tab)); if (lg(TABwp(tab)) != L) return 0;
     239 [ +  + ][ -  + ]:      18809 :   LN = lg(TABxm(tab)); if (LN != 1 && LN != L) return 0;
     240 [ +  + ][ -  + ]:      18809 :   LW = lg(TABwm(tab)); if (LW != 1 && LW != L) return 0;
     241                 :      18809 :   return 1;
     242                 :            : }
     243                 :            : 
     244                 :            : static int
     245                 :        259 : checktabdoub(GEN tab)
     246                 :            : {
     247                 :            :   long L;
     248         [ -  + ]:        259 :   if (typ(tab) != t_VEC) return 0;
     249         [ -  + ]:        259 :   if (lg(tab) != 8) return 0;
     250         [ -  + ]:        259 :   if (typ(TABm(tab)) != t_INT) return 0;
     251                 :        259 :   L = lg(TABxp(tab));
     252         [ -  + ]:        259 :   if (lg(TABwp(tab)) != L) return 0;
     253         [ -  + ]:        259 :   if (lg(TABxm(tab)) != L) return 0;
     254         [ -  + ]:        259 :   if (lg(TABwm(tab)) != L) return 0;
     255                 :        259 :   return 1;
     256                 :            : }
     257                 :            : 
     258                 :            : static int
     259                 :       9380 : checktab(GEN tab)
     260                 :            : {
     261         [ -  + ]:       9380 :   if (typ(tab) != t_VEC) return 0;
     262         [ +  - ]:       9380 :   if (lg(tab) != 3) return checktabsimp(tab);
     263                 :       9380 :   return checktabsimp(gel(tab,1))
     264 [ #  # ][ #  # ]:          0 :       && checktabsimp(gel(tab,2));
     265                 :            : }
     266                 :            : 
     267                 :            : static long
     268                 :        287 : findmforinit(long m, long prec)
     269                 :            : {
     270                 :            :   long p, r;
     271                 :            : 
     272         [ +  + ]:        287 :   if (m <= 0)
     273                 :            :   {
     274                 :        266 :     p = (long)prec2nbits_mul(prec, 0.3);
     275                 :        266 :     m = 2; r = 4;
     276         [ +  + ]:       1631 :     while (r < p) { m++; r <<= 1; }
     277                 :            :   }
     278                 :        287 :   return m;
     279                 :            : }
     280                 :            : 
     281                 :            : long
     282                 :         14 : intnumstep(long prec) { return findmforinit(0, prec); }
     283                 :            : 
     284                 :            : static void
     285                 :        273 : intinit_start(intdata *D, long m0, long flext, long prec)
     286                 :            : {
     287                 :        273 :   long m = findmforinit(m0, prec), lim = 20L<<m;
     288         [ -  + ]:        273 :   if (flext > 0) lim = lim << (2*flext);
     289                 :        273 :   D->m = m;
     290                 :        273 :   D->eps = prec2nbits(prec);
     291                 :        273 :   D->tabxp = cgetg(lim+1, t_VEC);
     292                 :        273 :   D->tabwp = cgetg(lim+1, t_VEC);
     293                 :        273 :   D->tabxm = cgetg(lim+1, t_VEC);
     294                 :        273 :   D->tabwm = cgetg(lim+1, t_VEC);
     295                 :        273 : }
     296                 :            : 
     297                 :            : static GEN
     298                 :        273 : intinit_end(intdata *D, long pnt, long mnt)
     299                 :            : {
     300                 :        273 :   GEN v = cgetg(8, t_VEC);
     301         [ -  + ]:        273 :   if (pnt < 0) pari_err_DOMAIN("intnuminit","table length","<",gen_0,stoi(pnt));
     302                 :        273 :   gel(v,1) = stoi(D->m);
     303                 :        273 :   TABx0(v) = D->tabx0;
     304                 :        273 :   TABw0(v) = D->tabw0;
     305                 :        273 :   TABxp(v) = D->tabxp; setlg(D->tabxp, pnt+1);
     306                 :        273 :   TABwp(v) = D->tabwp; setlg(D->tabwp, pnt+1);
     307                 :        273 :   TABxm(v) = D->tabxm; setlg(D->tabxm, mnt+1);
     308                 :        273 :   TABwm(v) = D->tabwm; setlg(D->tabwm, mnt+1); return v;
     309                 :            : }
     310                 :            : 
     311                 :            : static const long EXTRAPREC =
     312                 :            : #ifdef LONG_IS_64BIT
     313                 :            :   1;
     314                 :            : #else
     315                 :            :   2;
     316                 :            : #endif
     317                 :            : 
     318                 :            : /* divide by 2 in place */
     319                 :            : static GEN
     320                 :     549710 : divr2_ip(GEN x) { shiftr_inplace(x, -1); return x; }
     321                 :            : 
     322                 :            : /* phi(t)=tanh((3/2)sinh(t)) : from -1 to 1, hence also from a to b compact
     323                 :            :  * interval. */
     324                 :            : static GEN
     325                 :         63 : inittanhsinh(long m, long prec)
     326                 :            : {
     327                 :         63 :   pari_sp av, ltop = avma;
     328                 :            :   GEN h, et, ct, st, ext, ex, xp, wp;
     329                 :         63 :   long k, nt = -1, lim;
     330                 :         63 :   intdata D; intinit_start(&D, m, 0, prec);
     331                 :            : 
     332                 :         63 :   lim = lg(D.tabxp) - 1;
     333                 :         63 :   D.tabx0 = real_0(prec);
     334                 :         63 :   D.tabw0 = divr2_ip(stor(3, prec));
     335                 :         63 :   h = real2n(-D.m, prec);
     336                 :         63 :   et = ex = mpexp(h);
     337         [ +  - ]:      33537 :   for (k = 1; k <= lim; k++)
     338                 :            :   {
     339                 :      33537 :     gel(D.tabxp,k) = cgetr(prec+EXTRAPREC);
     340                 :      33537 :     gel(D.tabwp,k) = cgetr(prec+EXTRAPREC); av = avma;
     341                 :      33537 :     ct = divr2_ip(addrr(et, invr(et)));
     342                 :      33537 :     st = subrr(et, ct);
     343                 :      33537 :     ext = invr( addrs(mpexp(mulur(3, st)), 1) );
     344                 :      33537 :     shiftr_inplace(ext, 1);
     345                 :      33537 :     xp = subsr(1, ext);
     346                 :      33537 :     wp = divr2_ip(mulur(3, mulrr(ct, mulrr(ext, addsr(1, xp)))));
     347         [ +  + ]:      33537 :     if (expo(wp) < -D.eps) { nt = k-1; break; }
     348                 :      33474 :     affrr(xp, gel(D.tabxp,k));
     349                 :      33474 :     affrr(wp, gel(D.tabwp,k)); et = gerepileuptoleaf(av, mulrr(et, ex));
     350                 :            :   }
     351                 :         63 :   return gerepilecopy(ltop, intinit_end(&D, nt, 0));
     352                 :            : }
     353                 :            : 
     354                 :            : /* phi(t)=sinh(sinh(t)) : from -\infty to \infty, slowly decreasing, at least
     355                 :            :  * as 1/x^2. */
     356                 :            : static GEN
     357                 :         21 : initsinhsinh(long m, long prec)
     358                 :            : {
     359                 :         21 :   pari_sp av, ltop = avma;
     360                 :            :   GEN h, et, ct, st, ext, exu, ex, xp, wp;
     361                 :         21 :   long k, nt = -1, lim;
     362                 :         21 :   intdata D; intinit_start(&D, m, 0, prec);
     363                 :            : 
     364                 :         21 :   lim = lg(D.tabxp) - 1;
     365                 :         21 :   D.tabx0 = real_0(prec);
     366                 :         21 :   D.tabw0 = real_1(prec);
     367                 :         21 :   h = real2n(-D.m, prec);
     368                 :         21 :   et = ex = mpexp(h);
     369         [ +  - ]:      78813 :   for (k = 1; k <= lim; k++)
     370                 :            :   {
     371                 :      78813 :     gel(D.tabxp,k) = cgetr(prec+EXTRAPREC);
     372                 :      78813 :     gel(D.tabwp,k) = cgetr(prec+EXTRAPREC); av = avma;
     373                 :      78813 :     ct = divr2_ip(addrr(et, invr(et)));
     374                 :      78813 :     st = subrr(et, ct);
     375                 :      78813 :     ext = mpexp(st);
     376                 :      78813 :     exu = invr(ext);
     377                 :      78813 :     xp = divr2_ip(subrr(ext, exu));
     378                 :      78813 :     wp = divr2_ip(mulrr(ct, addrr(ext, exu)));
     379         [ +  + ]:      78813 :     if (expo(wp) - 2*expo(xp) < -D.eps) { nt = k-1; break; }
     380                 :      78792 :     affrr(xp, gel(D.tabxp,k));
     381                 :      78792 :     affrr(wp, gel(D.tabwp,k)); et = gerepileuptoleaf(av, mulrr(et, ex));
     382                 :            :   }
     383                 :         21 :   return gerepilecopy(ltop, intinit_end(&D, nt, 0));
     384                 :            : }
     385                 :            : 
     386                 :            : /* phi(t)=2sinh(t) : from -\infty to \infty, exponentially decreasing as
     387                 :            :  * exp(-x). */
     388                 :            : static GEN
     389                 :         14 : initsinh(long m, long prec)
     390                 :            : {
     391                 :         14 :   pari_sp av, ltop = avma;
     392                 :            :   GEN h, et, ex, eti, xp, wp;
     393                 :         14 :   long k, nt = -1, lim;
     394                 :         14 :   intdata D; intinit_start(&D, m, 0, prec);
     395                 :            : 
     396                 :         14 :   lim = lg(D.tabxp) - 1;
     397                 :         14 :   D.tabx0 = real_0(prec);
     398                 :         14 :   D.tabw0 = real2n(1, prec);
     399                 :         14 :   h = real2n(-D.m, prec);
     400                 :         14 :   et = ex = mpexp(h);
     401         [ +  - ]:      10052 :   for (k = 1; k <= lim; k++)
     402                 :            :   {
     403                 :      10052 :     gel(D.tabxp,k) = cgetr(prec+EXTRAPREC);
     404                 :      10052 :     gel(D.tabwp,k) = cgetr(prec+EXTRAPREC); av = avma;
     405                 :      10052 :     eti = invr(et);
     406                 :      10052 :     xp = subrr(et, eti);
     407                 :      10052 :     wp = addrr(et, eti);
     408         [ +  + ]:      10052 :     if (cmprs(xp, (long)(LOG2*(expo(wp)+D.eps) + 1)) > 0) { nt = k-1; break; }
     409                 :      10038 :     affrr(xp, gel(D.tabxp,k));
     410                 :      10038 :     affrr(wp, gel(D.tabwp,k)); et = gerepileuptoleaf(av, mulrr(et, ex));
     411                 :            :   }
     412                 :         14 :   return gerepilecopy(ltop, intinit_end(&D, nt, 0));
     413                 :            : }
     414                 :            : 
     415                 :            : /* phi(t)=exp(2sinh(t)) : from 0 to \infty, slowly decreasing at least as
     416                 :            :  * 1/x^2. */
     417                 :            : static GEN
     418                 :         21 : initexpsinh(long m, long prec)
     419                 :            : {
     420                 :         21 :   pari_sp ltop = avma;
     421                 :            :   GEN h, et, eti, ex, xp;
     422                 :         21 :   long k, nt = -1, lim;
     423                 :         21 :   intdata D; intinit_start(&D, m, 0, prec);
     424                 :            : 
     425                 :         21 :   lim = lg(D.tabxp) - 1;
     426                 :         21 :   D.tabx0 = real_1(prec);
     427                 :         21 :   D.tabw0 = real2n(1, prec);
     428                 :         21 :   h = real2n(-D.m, prec);
     429                 :         21 :   ex = mpexp(h);
     430                 :         21 :   et = real_1(prec);
     431         [ +  - ]:      30772 :   for (k = 1; k <= lim; k++)
     432                 :            :   {
     433                 :            :     GEN t;
     434                 :      30772 :     et = mulrr(et, ex);
     435                 :      30772 :     eti = invr(et); t = addrr(et, eti);
     436                 :      30772 :     xp = mpexp(subrr(et, eti));
     437                 :      30772 :     gel(D.tabxp,k) = xp;
     438                 :      30772 :     gel(D.tabwp,k) = mulrr(xp, t);
     439                 :      30772 :     gel(D.tabxm,k) = invr(xp);
     440                 :      30772 :     gel(D.tabwm,k) = mulrr(gel(D.tabxm,k), t);
     441         [ +  + ]:      30772 :     if (expo(gel(D.tabxm,k)) < -D.eps) { nt = k-1; break; }
     442                 :            :   }
     443                 :         21 :   return gerepilecopy(ltop, intinit_end(&D, nt, nt));
     444                 :            : }
     445                 :            : 
     446                 :            : /* phi(t)=exp(t-exp(-t)) : from 0 to \infty, exponentially decreasing. */
     447                 :            : static GEN
     448                 :         56 : initexpexp(long m, long prec)
     449                 :            : {
     450                 :         56 :   pari_sp av, ltop = avma;
     451                 :            :   GEN kh, h, et, eti, ex, xp, xm, wp, wm;
     452                 :         56 :   long k, nt = -1, lim;
     453                 :         56 :   intdata D; intinit_start(&D, m, 0, prec);
     454                 :            : 
     455                 :         56 :   lim = lg(D.tabxp) - 1;
     456                 :         56 :   D.tabx0 = mpexp(real_m1(prec));
     457                 :         56 :   D.tabw0 = gmul2n(D.tabx0, 1);
     458                 :         56 :   h = real2n(-D.m, prec);
     459                 :         56 :   et = ex = mpexp(negr(h));
     460         [ +  - ]:      40264 :   for (k = 1; k <= lim; k++)
     461                 :            :   {
     462                 :      40264 :     gel(D.tabxp,k) = cgetr(prec+EXTRAPREC);
     463                 :      40264 :     gel(D.tabwp,k) = cgetr(prec+EXTRAPREC);
     464                 :      40264 :     gel(D.tabxm,k) = cgetr(prec+EXTRAPREC);
     465                 :      40264 :     gel(D.tabwm,k) = cgetr(prec+EXTRAPREC); av = avma;
     466                 :      40264 :     eti = invr(et); kh = mulur(k,h);
     467                 :      40264 :     xp = mpexp(subrr(kh, et));
     468                 :      40264 :     xm = mpexp(negr(addrr(kh, eti)));
     469                 :      40264 :     wp = mulrr(xp, addsr(1, et));
     470                 :      40264 :     wm = mulrr(xm, addsr(1, eti));
     471 [ +  + ][ +  + ]:      40264 :     if (expo(xm) < -D.eps && cmprs(xp, (long)(LOG2*(expo(wp)+D.eps) + 1)) > 0) { nt = k-1; break; }
     472                 :      40208 :     affrr(xp, gel(D.tabxp,k));
     473                 :      40208 :     affrr(wp, gel(D.tabwp,k));
     474                 :      40208 :     affrr(xm, gel(D.tabxm,k));
     475                 :      40208 :     affrr(wm, gel(D.tabwm,k)); et = gerepileuptoleaf(av, mulrr(et, ex));
     476                 :            :   }
     477                 :         56 :   return gerepilecopy(ltop, intinit_end(&D, nt, nt));
     478                 :            : }
     479                 :            : 
     480                 :            : /* phi(t)=(Pi/h)t/(1-exp(-sinh(t))) : from 0 to \infty, sine oscillation. */
     481                 :            : static GEN
     482                 :         84 : initnumsine(long m, long prec)
     483                 :            : {
     484                 :         84 :   pari_sp av, ltop = avma;
     485                 :            :   GEN h, et, eti, ex, st, ct, extp, extm, extp1, extm1, extp2, extm2, kpi, kct;
     486                 :         84 :   GEN xp, xm, wp, wm, pi = mppi(prec);
     487                 :         84 :   long k, nt = -1, lim;
     488                 :         84 :   intdata D; intinit_start(&D, m, 0, prec);
     489                 :            : 
     490                 :         84 :   lim = lg(D.tabxp) - 1;
     491                 :         84 :   D.tabx0 = gmul2n(pi, D.m);
     492                 :         84 :   D.tabw0 = gmul2n(pi, D.m - 1);
     493                 :         84 :   h = real2n(-D.m, prec);
     494                 :         84 :   et = ex = mpexp(h);
     495         [ +  - ]:     123067 :   for (k = 1; k <= lim; k++)
     496                 :            :   {
     497                 :     123067 :     gel(D.tabxp,k) = cgetr(prec+EXTRAPREC);
     498                 :     123067 :     gel(D.tabwp,k) = cgetr(prec+EXTRAPREC);
     499                 :     123067 :     gel(D.tabxm,k) = cgetr(prec+EXTRAPREC);
     500                 :     123067 :     gel(D.tabwm,k) = cgetr(prec+EXTRAPREC); av = avma;
     501                 :     123067 :     eti = invr(et); /* exp(-kh) */
     502                 :     123067 :     ct = divr2_ip(addrr(et, eti));
     503                 :     123067 :     st = divr2_ip(subrr(et, eti));
     504                 :     123067 :     extp = mpexp(st);  extp1 = subsr(1, extp); extp2 = invr(extp1);
     505                 :     123067 :     extm = invr(extp); extm1 = subsr(1, extm); extm2 = invr(extm1);
     506                 :     123067 :     kpi = mulur(k, pi);
     507                 :     123067 :     kct = mulur(k, ct);
     508                 :     123067 :     shiftr_inplace(extm1, D.m);
     509                 :     123067 :     shiftr_inplace(extp1, D.m);
     510                 :     123067 :     xp = mulrr(kpi, extm2);
     511                 :     123067 :     wp = mulrr(subrr(extm1, mulrr(kct, extm)), mulrr(pi, sqrr(extm2)));
     512                 :     123067 :     xm = mulrr(negr(kpi), extp2);
     513                 :     123067 :     wm = mulrr(addrr(extp1, mulrr(kct, extp)), mulrr(pi, sqrr(extp2)));
     514 [ +  + ][ +  + ]:     123067 :     if (expo(wm) < -D.eps && expo(extm) + D.m + expu(10 * k) < -D.eps) { nt = k-1; break; }
     515                 :     122983 :     affrr(xp, gel(D.tabxp,k));
     516                 :     122983 :     affrr(wp, gel(D.tabwp,k));
     517                 :     122983 :     affrr(xm, gel(D.tabxm,k));
     518                 :     122983 :     affrr(wm, gel(D.tabwm,k)); et = gerepileuptoleaf(av, mulrr(et, ex));
     519                 :            :   }
     520                 :         84 :   return gerepilecopy(ltop, intinit_end(&D, nt, nt));
     521                 :            : }
     522                 :            : 
     523                 :            : static GEN
     524                 :         84 : suminit_start(GEN sig)
     525                 :            : {
     526                 :            :   GEN sig2;
     527                 :            : 
     528         [ +  + ]:         84 :   if (typ(sig) == t_VEC)
     529                 :            :   {
     530         [ -  + ]:         28 :     if (lg(sig) != 3) pari_err_TYPE("sumnum",sig);
     531                 :         28 :     sig2 = gel(sig,2);
     532                 :         28 :     sig  = gel(sig,1);
     533         [ -  + ]:         28 :     if (!isinR(sig2)) pari_err_TYPE("sumnum",sig2);
     534         [ +  + ]:         28 :     if (gsigne(sig2) > 0) sig2 = mulcxmI(sig2);
     535                 :            :   }
     536                 :         56 :   else sig2 = gen_0;
     537         [ -  + ]:         84 :   if (!isinR(sig)) pari_err_TYPE("sumnum",sig);
     538                 :         84 :   return mkvec2(mkvec(gen_1), sig2);
     539                 :            : }
     540                 :            : 
     541                 :            : /* phi(t) depending on sig[2] as in intnum, with weights phi'(t)tanh(Pi*phi(t))
     542                 :            :  * (sgn >= 0) or phi'(t)/cosh(Pi*phi(t)) (otherwise), for use in sumnumall.
     543                 :            :  * integrations are done from 0 to +infty (flii is set to 0), except if slowly
     544                 :            :    decreasing, from -infty to +infty (flii is set to 1). */
     545                 :            : GEN
     546                 :         35 : sumnuminit(GEN sig, long m, long sgn, long prec)
     547                 :            : {
     548                 :         35 :   pari_sp ltop = avma;
     549                 :         35 :   GEN b, t, tab, tabxp, tabwp, tabxm, tabwm, pi = mppi(prec);
     550                 :            :   long L, k, eps, flii;
     551                 :            : 
     552                 :         35 :   b = suminit_start(sig);
     553                 :         35 :   flii = gequal0(gel(b,2));
     554         [ +  + ]:         35 :   if (flii)
     555                 :         21 :     tab = intnuminit(mkvec(gen_m1), mkvec(gen_1), m, prec);
     556                 :            :   else
     557                 :         14 :     tab = intnuminit(gen_0, b, m, prec);
     558                 :         35 :   eps = prec2nbits(prec);
     559                 :         35 :   t = gmul(pi, TABx0(tab));
     560         [ +  + ]:         35 :   if (sgn < 0) TABw0(tab) = gdiv(TABw0(tab), gcosh(t, prec));
     561                 :         28 :   else         TABw0(tab) = gmul(TABw0(tab), gtanh(t, prec));
     562                 :         35 :   tabxp = TABxp(tab); L = lg(tabxp);
     563                 :         35 :   tabwp = TABwp(tab);
     564                 :         35 :   tabxm = TABxm(tab);
     565                 :         35 :   tabwm = TABwm(tab);
     566         [ +  + ]:     155197 :   for (k = 1; k < L; k++)
     567                 :            :   {
     568         [ +  + ]:     155162 :     if (cmprs(gel(tabxp,k), eps) < 0)
     569                 :            :     {
     570                 :      34055 :       t = mulrr(pi, gel(tabxp,k));
     571                 :      43834 :       gel(tabwp,k) = (sgn < 0)? divrr(gel(tabwp,k), gcosh(t, prec))
     572         [ +  + ]:      43834 :                               : mulrr(gel(tabwp,k), gtanh(t, prec));
     573                 :            :     }
     574                 :            :     else
     575         [ +  + ]:     121107 :       if (sgn < 0) gel(tabwp,k) = real_0_bit(-eps);
     576         [ +  + ]:     155162 :     if (!flii)
     577                 :            :     {
     578                 :      76370 :       t = mulrr(pi, gel(tabxm,k));
     579                 :      76370 :       gel(tabwm,k) = (sgn < 0)? divrr(gel(tabwm,k), gcosh(t, prec))
     580         [ -  + ]:      76370 :                               : mulrr(gel(tabwm,k), gtanh(t, prec));
     581                 :            :     }
     582                 :            :   }
     583                 :         35 :   return gerepilecopy(ltop, tab);
     584                 :            : }
     585                 :            : 
     586                 :            : /* End of initialization functions. These functions can be executed once
     587                 :            :  * and for all for a given accuracy, type of integral ([a,b], [a,\infty[ or
     588                 :            :  * ]-\infty,a], ]-\infty,\infty[) and of integrand in the noncompact case
     589                 :            :  * (slowly decreasing, exponentially decreasing, oscillating with a fixed
     590                 :            :  * oscillating factor such as sin(x)). */
     591                 :            : 
     592                 :            : /* In the following integration functions the parameters are as follows:
     593                 :            : * 1) The parameter denoted by m is the most crucial and difficult to
     594                 :            : * determine in advance: h = 1/2^m is the integration step size. Usually
     595                 :            : * m = floor(log(D)/log(2)), where D is the number of decimal digits of accuracy
     596                 :            : * is plenty for very regulat functions, for instance m = 6 for 100D, and m = 9
     597                 :            : * for 1000D, but values of m 1 or 2 less are often sufficient, while for
     598                 :            : * singular functions, 1 or 2 more may be necessary. The best test is to take 2
     599                 :            : * or 3 consecutive values of m and look. Note that the number of function
     600                 :            : * evaluations, hence the time doubles when m increases by 1. */
     601                 :            : 
     602                 :            : /* All inner functions such as intn, etc... must be called with a
     603                 :            :  * valid 'tab' table. The wrapper intnum provides a higher level interface */
     604                 :            : 
     605                 :            : /* compute $\int_a^b f(t)dt$ with [a,b] compact and f nonsingular. */
     606                 :            : static GEN
     607                 :       9373 : intn(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, GEN tab)
     608                 :            : {
     609                 :            :   GEN tabx0, tabw0, tabxp, tabwp;
     610                 :            :   GEN bpa, bma, bmb, S, SP, SM;
     611                 :            :   long m, k, L, i;
     612                 :       9373 :   pari_sp ltop = avma, av;
     613                 :            : 
     614         [ -  + ]:       9373 :   if (!checktabsimp(tab)) pari_err_TYPE("intnum",tab);
     615         [ -  + ]:       9373 :   if (!isinC(a)) pari_err_TYPE("intnum",a);
     616         [ -  + ]:       9373 :   if (!isinC(b)) pari_err_TYPE("intnum",b);
     617                 :       9373 :   m = itos(TABm(tab));
     618                 :       9373 :   tabx0 = TABx0(tab); tabw0 = TABw0(tab);
     619                 :       9373 :   tabxp = TABxp(tab); tabwp = TABwp(tab); L = lg(tabxp);
     620                 :       9373 :   bpa = gmul2n(gadd(b, a), -1); /* (b+a)/2 */
     621                 :       9373 :   bma = gsub(bpa, a); /* (b-a)/2 */
     622                 :       9373 :   bmb = gmul(bma, tabx0); /* (b-a)/2 phi(0) */
     623                 :       9373 :   av = avma;
     624                 :            :   /* phi'(0) f( (b+a)/2 + (b-a)/2 * phi(0) ) */
     625                 :       9373 :   S = gmul(tabw0, eval(E, gadd(bpa, bmb)));
     626         [ +  + ]:      74963 :   for (k = 1; k <= m; k++)
     627                 :            :   {
     628                 :      65590 :     long pas = 1L<<(m-k);
     629         [ +  + ]:   12393073 :     for (i = pas; i < L; i += pas)
     630 [ +  + ][ +  + ]:   12327483 :       if (i & pas || k == 1)
     631                 :            :       {
     632                 :    6224624 :         bmb = gmul(bma, gel(tabxp,i));
     633                 :    6224624 :         SP = eval(E, gsub(bpa, bmb));
     634                 :    6224624 :         SM = eval(E, gadd(bpa, bmb));
     635                 :    6224624 :         S = gadd(S, gmul(gel(tabwp,i), gadd(SP, SM)));
     636         [ +  + ]:    6224624 :         if ((i & 0x7f) == 1) S = gerepileupto(av, S);
     637                 :            :       }
     638                 :            :   }
     639                 :       9373 :   return gerepileupto(ltop, gmul(S, gmul2n(bma, -m)));
     640                 :            : }
     641                 :            : 
     642                 :            : /* compute $\int_{a[1]}^{b} f(t)dt$ with [a,b] compact, possible
     643                 :            :  *  singularity with exponent a[2] at lower extremity, b regular.
     644                 :            :  *  Use tanh(sinh(t)). */
     645                 :            : static GEN
     646                 :          7 : intnsing(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, GEN tab, long prec)
     647                 :            : {
     648                 :            :   GEN tabx0, tabw0, tabxp, tabwp, ea, ba, bm, bp, S, tra, SP, SM;
     649                 :            :   long m, k, L, i;
     650                 :          7 :   pari_sp ltop = avma, av;
     651                 :            : 
     652         [ -  + ]:          7 :   if (!checktabsimp(tab)) pari_err_TYPE("intnum",tab);
     653                 :          7 :   m = itos(TABm(tab));
     654                 :          7 :   tabx0 = TABx0(tab); tabw0 = TABw0(tab);
     655                 :          7 :   tabxp = TABxp(tab); tabwp = TABwp(tab); L = lg(tabxp);
     656                 :          7 :   tra = gel(a,1);
     657                 :          7 :   ea = ginv(gaddsg(1, gel(a,2)));
     658                 :          7 :   ba = gdiv(gsub(b, tra), gpow(gen_2, ea, prec));
     659                 :          7 :   av = avma;
     660                 :          7 :   S = gmul(gmul(tabw0, ba), eval(E, gadd(gmul(ba, gaddsg(1, tabx0)), tra)));
     661         [ +  + ]:         56 :   for (k = 1; k <= m; k++)
     662                 :            :   {
     663                 :         49 :     long pas = 1L<<(m-k);
     664         [ +  + ]:       9268 :     for (i = pas; i < L; i += pas)
     665 [ +  + ][ +  + ]:       9219 :       if (i & pas || k == 1) /* i = odd multiple of pas = 2^(m-k) */
     666                 :            :       {
     667                 :       4655 :         GEN p = addsr(1, gel(tabxp,i));
     668                 :       4655 :         GEN m = subsr(1, gel(tabxp,i));
     669                 :       4655 :         bp = gmul(ba, gpow(p, ea, prec));
     670                 :       4655 :         bm = gmul(ba, gpow(m, ea, prec));
     671                 :       4655 :         SP = gmul(gdiv(bp, p), eval(E, gadd(bp, tra)));
     672                 :       4655 :         SM = gmul(gdiv(bm, m), eval(E, gadd(bm, tra)));
     673                 :       4655 :         S = gadd(S, gmul(gel(tabwp,i), gadd(SP, SM)));
     674         [ +  + ]:       4655 :         if ((i & 0x7f) == 1) S = gerepileupto(av, S);
     675                 :            :       }
     676                 :            :   }
     677                 :          7 :   return gerepileupto(ltop, gmul(gmul2n(S, -m), ea));
     678                 :            : }
     679                 :            : 
     680                 :            : /* compute  $\int_a^\infty f(t)dt$ if $si=1$ or $\int_{-\infty}^a f(t)dt$
     681                 :            :    if $si=-1$. Use exp(2sinh(t)) for slowly decreasing functions,
     682                 :            :    exp(1+t-exp(-t)) for exponentially decreasing functions, and
     683                 :            :    (pi/h)t/(1-exp(-sinh(t))) for oscillating functions. */
     684                 :            : 
     685                 :            : static GEN
     686                 :        259 : intninfpm(void *E, GEN (*eval)(void*, GEN), GEN a, long si, GEN tab)
     687                 :            : {
     688                 :            :   GEN tabx0, tabw0, tabxp, tabwp, tabxm, tabwm;
     689                 :            :   GEN S, SP, SM;
     690                 :        259 :   long m, L, k, h = 0, pas, i;
     691                 :        259 :   pari_sp ltop = avma, av;
     692                 :            : 
     693         [ -  + ]:        259 :   if (!checktabdoub(tab)) pari_err_TYPE("intnum",tab);
     694                 :        259 :   m = itos(TABm(tab));
     695                 :        259 :   tabx0 = TABx0(tab); tabw0 = TABw0(tab);
     696                 :        259 :   tabxp = TABxp(tab); tabwp = TABwp(tab); L = lg(tabxp);
     697                 :        259 :   tabxm = TABxm(tab); tabwm = TABwm(tab);
     698         [ +  + ]:        259 :   if (si < 0) { tabxp = gneg(tabxp); tabxm = gneg(tabxm); }
     699                 :        259 :   av = avma;
     700                 :        259 :   S = gmul(tabw0, eval(E, gadd(a, gmulsg(si, tabx0))));
     701         [ +  + ]:       2128 :   for (k = 1; k <= m; k++)
     702                 :            :   {
     703                 :       1869 :     h++; pas = 1L<<(m-k);
     704         [ +  + ]:     575946 :     for (i = pas; i < L; i += pas)
     705 [ +  + ][ +  + ]:     574077 :       if (i & pas || k == 1)
     706                 :            :       {
     707                 :     288939 :         SP = eval(E, gadd(a, gel(tabxp,i)));
     708                 :     288939 :         SM = eval(E, gadd(a, gel(tabxm,i)));
     709                 :     288939 :         S = gadd(S, gadd(gmul(gel(tabwp,i), SP), gmul(gel(tabwm,i), SM)));
     710         [ +  + ]:     288939 :         if ((i & 0x7f) == 1) S = gerepileupto(av, S);
     711                 :            :       }
     712                 :            :   }
     713                 :        259 :   return gerepileupto(ltop, gmul2n(S, -h));
     714                 :            : }
     715                 :            : 
     716                 :            : /* compute  $\int_{-\infty}^\infty f(t)dt$
     717                 :            :  * use sinh(sinh(t)) for slowly decreasing functions and sinh(t) for
     718                 :            :  * exponentially decreasing functions.
     719                 :            :  * HACK: in case TABwm(tab) contains something, assume function to be integrated
     720                 :            :  * satisfies f(-x) = conj(f(x)).
     721                 :            :  * Usually flag < 0, but flag > 0 is used in sumnumall. */
     722                 :            : static GEN
     723                 :         49 : intninfinfintern(void *E, GEN (*eval)(void*, GEN), GEN tab, long flag)
     724                 :            : {
     725                 :            :   GEN tabx0, tabw0, tabxp, tabwp, tabwm;
     726                 :            :   GEN S, SP, SM;
     727                 :            :   long m, L, k, i, spf;
     728                 :         49 :   pari_sp ltop = avma;
     729                 :            : 
     730         [ -  + ]:         49 :   if (!checktabsimp(tab)) pari_err_TYPE("intnum",tab);
     731                 :         49 :   m = itos(TABm(tab));
     732                 :         49 :   tabx0 = TABx0(tab); tabw0 = TABw0(tab);
     733                 :         49 :   tabxp = TABxp(tab); tabwp = TABwp(tab); L = lg(tabxp);
     734                 :         49 :   tabwm = TABwm(tab);
     735                 :         49 :   spf = (lg(tabwm) == lg(tabwp));
     736         [ +  + ]:         49 :   S = flag > 0 ? gen_0 : gmul(tabw0, eval(E, tabx0));
     737         [ +  + ]:         49 :   if (spf) S = gmul2n(real_i(S), -1);
     738         [ +  + ]:        462 :   for (k = 1; k <= m; k++)
     739                 :            :   {
     740                 :        413 :     long pas = 1L<<(m-k);
     741         [ +  + ]:     282331 :     for (i = pas; i < L; i += pas)
     742 [ +  + ][ +  + ]:     281918 :       if (i & pas || k == 1)
     743                 :            :       {
     744                 :     141351 :         SP = eval(E, gel(tabxp,i));
     745         [ +  + ]:     141351 :         if (spf) S = gadd(S, real_i(gmul(gel(tabwp,i), SP)));
     746                 :            :         else
     747                 :            :         {
     748                 :      31283 :           SM = eval(E, negr(gel(tabxp,i)));
     749         [ +  + ]:      31283 :           if (flag > 0) SM = gneg(SM);
     750                 :      31283 :           S = gadd(S, gmul(gel(tabwp,i), gadd(SP, SM)));
     751                 :            :         }
     752         [ +  + ]:     141351 :         if ((i & 0x7f) == 1) S = gerepileupto(ltop, S);
     753                 :            :       }
     754                 :            :   }
     755         [ +  + ]:         49 :   if (spf) m--;
     756                 :         49 :   return gerepileupto(ltop, gmul2n(S, -m));
     757                 :            : }
     758                 :            : 
     759                 :            : static GEN
     760                 :         42 : intninfinf(void *E, GEN (*eval)(void*, GEN), GEN tab)
     761                 :            : {
     762                 :         42 :   return intninfinfintern(E, eval, tab, -1);
     763                 :            : }
     764                 :            : 
     765                 :            : /* general num integration routine int_a^b f(t)dt, where a and b are as follows:
     766                 :            :  (1) a scalar : the scalar, no singularity worse than logarithmic at a.
     767                 :            :  (2) [a, e] : the scalar a, singularity exponent -1 < e <= 0.
     768                 :            :  (3) [1], [-1] : +\infty, -\infty, slowly decreasing function.
     769                 :            :  (4) [[+-1], a], a nonnegative real : +-\infty, function behaving like
     770                 :            :       exp(-a|t|) at +-\infty.
     771                 :            :  (5) [[+-1], e], e < -1 : +-\infty, function behaving like t^e
     772                 :            :       at +-\infty.
     773                 :            :  (5) [[+-1], a*I], a real : +-\infty, function behaving like cos(at) if a>0
     774                 :            :      and like sin(at) if a < 0 at +-\infty.
     775                 :            : */
     776                 :            : 
     777                 :            : /* FIXME: The numbers below can be changed, but NOT the ordering */
     778                 :            : enum {
     779                 :            :   f_REG    = 0, /* regular function */
     780                 :            :   f_SING   = 1, /* algebraic singularity */
     781                 :            :   f_YSLOW  = 2, /* +\infty, slowly decreasing */
     782                 :            :   f_YVSLO  = 3, /* +\infty, very slowly decreasing */
     783                 :            :   f_YFAST  = 4, /* +\infty, exponentially decreasing */
     784                 :            :   f_YOSCS  = 5, /* +\infty, sine oscillating */
     785                 :            :   f_YOSCC  = 6  /* +\infty, cosine oscillating */
     786                 :            : };
     787                 :            : /* is c finite */
     788                 :            : static int
     789 [ +  + ][ +  + ]:        560 : is_fin_f(long c) { return c == f_REG || c == f_SING; }
     790                 :            : /* oscillating case: valid for +oo (c > 0) or -oo (c < 0) */
     791                 :            : static int
     792 [ +  - ][ -  + ]:      10066 : is_osc_f(long c) { c = labs(c); return c == f_YOSCS || c == f_YOSCC; }
     793                 :            : 
     794                 :            : static GEN
     795                 :        406 : f_getycplx(GEN a, long prec)
     796                 :            : {
     797                 :            :   long s;
     798                 :            :   GEN tmp, a2R, a2I;
     799                 :            : 
     800 [ +  + ][ -  + ]:        406 :   if (lg(a) == 2 || gequal0(gel(a,2))) return gen_1;
     801                 :        350 :   a2R = real_i(gel(a,2));
     802                 :        350 :   a2I = imag_i(gel(a,2));
     803         [ +  + ]:        350 :   s = gsigne(a2I); if (s < 0) a2I = gneg(a2I);
     804         [ +  + ]:        350 :   tmp = s ? ginv(a2I) : ginv(a2R);
     805         [ +  - ]:        350 :   if (gprecision(tmp) < prec) tmp = gprec_w(tmp, prec);
     806                 :        406 :   return tmp;
     807                 :            : }
     808                 :            : 
     809                 :            : static void
     810                 :          0 : err_code(GEN a, const char *name)
     811                 :            : {
     812                 :          0 :   char *s = stack_sprintf("intnum [incorrect %s]", name);
     813                 :          0 :   pari_err_TYPE(s, a);
     814                 :          0 : }
     815                 :            : 
     816                 :            : /* a = [[+/-1], alpha]*/
     817                 :            : static long
     818                 :        497 : code_aux(GEN a, const char *name)
     819                 :            : {
     820                 :        497 :   GEN re, im, alpha = gel(a,2);
     821                 :            :   long s;
     822         [ -  + ]:        497 :   if (!isinC(alpha)) err_code(a, name);
     823                 :        497 :   re = real_i(alpha);
     824                 :        497 :   im = imag_i(alpha);
     825                 :        497 :   s = gsigne(im);
     826         [ +  + ]:        497 :   if (s)
     827                 :            :   {
     828         [ -  + ]:        273 :     if(!gequal0(re))
     829                 :          0 :       pari_warn(warner,"real(z)*imag(z)!=0 in endpoint code, real(z) ignored");
     830         [ +  + ]:        273 :     return s > 0 ? f_YOSCC : f_YOSCS;
     831                 :            :   }
     832 [ +  - ][ -  + ]:        224 :   if (gequal0(re) || gcmpgs(re, -2)<=0) return f_YSLOW;
     833         [ +  + ]:        224 :   if (gsigne(re) > 0) return f_YFAST;
     834         [ -  + ]:         14 :   if (gcmpgs(re, -1) >= 0) err_code(a, name);
     835                 :        497 :   return f_YVSLO;
     836                 :            : }
     837                 :            : 
     838                 :            : static long
     839                 :      19768 : transcode(GEN a, const char *name)
     840                 :            : {
     841                 :            :   GEN a1, a2;
     842      [ +  +  + ]:      19768 :   switch(typ(a))
     843                 :            :   {
     844                 :        581 :     case t_VEC: break;
     845         [ +  - ]:         21 :     case t_INFINITY: return inf_get_sign(a) == 1 ? f_YSLOW: -f_YSLOW;
     846                 :      19166 :     default: return f_REG;
     847                 :            :   }
     848         [ -  + ]:        581 :   if (typ(a) != t_VEC) return f_REG;
     849      [ +  +  - ]:        581 :   switch(lg(a))
     850                 :            :   {
     851         [ +  + ]:         70 :     case 2: return gsigne(gel(a,1)) > 0 ? f_YSLOW : -f_YSLOW;
     852                 :        511 :     case 3: break;
     853                 :          0 :     default: err_code(a,name);
     854                 :            :   }
     855                 :        511 :   a1 = gel(a,1);
     856                 :        511 :   a2 = gel(a,2);
     857      [ +  +  + ]:        511 :   switch(typ(a1))
     858                 :            :   {
     859                 :            :     case t_VEC:
     860         [ -  + ]:        350 :       if (lg(a1) != 2) err_code(a,name);
     861                 :        350 :       return gsigne(gel(a1,1)) * code_aux(a, name);
     862                 :            :     case t_INFINITY:
     863                 :        147 :       return inf_get_sign(a1) * code_aux(a, name);
     864                 :            :     default:
     865 [ +  - ][ +  - ]:         14 :       if (!isinC(a1) || !isinR(a2) || gcmpgs(a2, -1) <= 0) err_code(a,name);
                 [ -  + ]
     866                 :      19768 :       return gsigne(a2) < 0 ? f_SING : f_REG;
     867                 :            :   }
     868                 :            : }
     869                 :            : 
     870                 :            : /* computes the necessary tabs, knowing a, b and m */
     871                 :            : static GEN
     872                 :        182 : homtab(GEN tab, GEN k)
     873                 :            : {
     874                 :            :   GEN z;
     875 [ +  - ][ +  + ]:        182 :   if (gequal0(k) || gequal(k, gen_1)) return tab;
     876         [ -  + ]:         21 :   if (gsigne(k) < 0) k = gneg(k);
     877                 :         21 :   z = cgetg(8, t_VEC);
     878                 :         21 :   TABm(z)  = icopy(TABm(tab));
     879                 :         21 :   TABx0(z) = gmul(TABx0(tab), k);
     880                 :         21 :   TABw0(z) = gmul(TABw0(tab), k);
     881                 :         21 :   TABxp(z) = gmul(TABxp(tab), k);
     882                 :         21 :   TABwp(z) = gmul(TABwp(tab), k);
     883                 :         21 :   TABxm(z) = gmul(TABxm(tab), k);
     884                 :        182 :   TABwm(z) = gmul(TABwm(tab), k); return z;
     885                 :            : }
     886                 :            : 
     887                 :            : static GEN
     888                 :         14 : expvec(GEN v, GEN ea, long prec)
     889                 :            : {
     890                 :         14 :   long lv = lg(v), i;
     891                 :         14 :   GEN z = cgetg(lv, t_VEC);
     892         [ +  + ]:      47502 :   for (i = 1; i < lv; i++) gel(z,i) = gpow(gel(v,i),ea,prec);
     893                 :         14 :   return z;
     894                 :            : }
     895                 :            : 
     896                 :            : static GEN
     897                 :      47495 : expscalpr(GEN vnew, GEN xold, GEN wold, GEN ea)
     898                 :            : {
     899                 :      47495 :   pari_sp av = avma;
     900                 :      47495 :   return gerepileupto(av, gdiv(gmul(gmul(vnew, wold), ea), xold));
     901                 :            : }
     902                 :            : static GEN
     903                 :         14 : expvecpr(GEN vnew, GEN xold, GEN wold, GEN ea)
     904                 :            : {
     905                 :         14 :   long lv = lg(vnew), i;
     906                 :         14 :   GEN z = cgetg(lv, t_VEC);
     907         [ +  + ]:      47502 :   for (i = 1; i < lv; i++)
     908                 :      47488 :     gel(z,i) = expscalpr(gel(vnew,i), gel(xold,i), gel(wold,i), ea);
     909                 :         14 :   return z;
     910                 :            : }
     911                 :            : 
     912                 :            : /* here k < -1 */
     913                 :            : static GEN
     914                 :          7 : exptab(GEN tab, GEN k, long prec)
     915                 :            : {
     916                 :            :   GEN v, ea;
     917                 :            : 
     918         [ -  + ]:          7 :   if (gcmpgs(k, -2) <= 0) return tab;
     919                 :          7 :   ea = ginv(gsubsg(-1, k));
     920                 :          7 :   v = cgetg(8, t_VEC);
     921                 :          7 :   TABm(v) = icopy(TABm(tab));
     922                 :          7 :   TABx0(v) = gpow(TABx0(tab), ea, prec);
     923                 :          7 :   TABw0(v) = expscalpr(TABx0(v), TABx0(tab), TABw0(tab), ea);
     924                 :          7 :   TABxp(v) = expvec(TABxp(tab), ea, prec);
     925                 :          7 :   TABwp(v) = expvecpr(TABxp(v), TABxp(tab), TABwp(tab), ea);
     926                 :          7 :   TABxm(v) = expvec(TABxm(tab), ea, prec);
     927                 :          7 :   TABwm(v) = expvecpr(TABxm(v), TABxm(tab), TABwm(tab), ea);
     928                 :          7 :   return v;
     929                 :            : }
     930                 :            : 
     931                 :            : GEN
     932                 :        252 : intnuminit(GEN a, GEN b, long m, long prec)
     933                 :            : {
     934                 :            :   long codea, codeb, l;
     935                 :            :   GEN T, U, km, kma, kmb, tmp;
     936                 :            : 
     937         [ -  + ]:        252 :   if (m > 30) pari_err_OVERFLOW("intnuminit [m]");
     938                 :        252 :   l = prec+EXTRAPREC;
     939                 :        252 :   codea = transcode(a, "a");
     940                 :        252 :   codeb = transcode(b, "b");
     941 [ +  + ][ +  + ]:        252 :   if (is_fin_f(codea) && is_fin_f(codeb)) return inittanhsinh(m, l);
     942         [ -  + ]:        196 :   if (labs(codea) > labs(codeb)) { swap(a, b); lswap(codea, codeb); }
     943         [ +  + ]:        196 :   if (codea == f_REG)
     944                 :            :   {
     945                 :        133 :     km = f_getycplx(b, l);
     946   [ +  +  +  +  :        133 :     switch(labs(codeb))
                   +  - ]
     947                 :            :     {
     948                 :         14 :       case f_YSLOW: return initexpsinh(m, l);
     949                 :          7 :       case f_YVSLO: return exptab(initexpsinh(m, l), gel(b,2), prec);
     950                 :         56 :       case f_YFAST: return homtab(initexpexp(m, l), km);
     951                 :            :       case f_YOSCS:
     952 [ +  - ][ +  - ]:         49 :         if (typ(a) == t_VEC || gequal0(a)) return homtab(initnumsine(m, l), km);
     953                 :            :             /* fall through */
     954                 :            :       case f_YOSCC:
     955                 :          7 :         T = cgetg(3, t_VEC);
     956                 :          7 :         gel(T,1) = inittanhsinh(m, l);
     957                 :          7 :         gel(T,2) = homtab(initnumsine(m, l), km);
     958                 :          7 :         return T;
     959                 :            :     }
     960                 :            :   }
     961         [ -  + ]:         63 :   if (codea == f_SING)
     962                 :            :   {
     963                 :          0 :     km = f_getycplx(b, l);
     964                 :          0 :     T = cgetg(3, t_VEC);
     965                 :          0 :     gel(T,1) = inittanhsinh(m, l);
     966   [ #  #  #  #  :          0 :     switch(labs(codeb))
                      # ]
     967                 :            :     {
     968                 :          0 :       case f_YSLOW: gel(T,2) = initexpsinh(m, l); break;
     969                 :          0 :       case f_YVSLO: gel(T,2) = exptab(initexpsinh(m, l), gel(b,2), prec); break;
     970                 :          0 :       case f_YFAST: gel(T,2) = homtab(initexpexp(m, l), km); break;
     971                 :            :       case f_YOSCS: case f_YOSCC:
     972                 :          0 :         gel(T,2) = homtab(initnumsine(m, l), km); break;
     973                 :            :     }
     974                 :          0 :     return T;
     975                 :            :   }
     976         [ -  + ]:         63 :   if (codea * codeb > 0) return gen_0;
     977                 :         63 :   kma = f_getycplx(a, l);
     978                 :         63 :   kmb = f_getycplx(b, l);
     979                 :         63 :   codea = labs(codea);
     980                 :         63 :   codeb = labs(codeb);
     981 [ +  + ][ +  - ]:         63 :   if (codea == f_YSLOW && codeb == f_YSLOW) return initsinhsinh(m, l);
     982 [ +  + ][ +  - ]:         42 :   if (codea == f_YFAST && codeb == f_YFAST && gequal(kma, kmb))
                 [ +  - ]
     983                 :         14 :     return homtab(initsinh(m, l), kmb);
     984                 :         28 :   T = cgetg(3, t_VEC);
     985   [ -  -  -  +  :         28 :   switch (codea)
                      - ]
     986                 :            :   {
     987                 :          0 :     case f_YSLOW: gel(T,1) = initexpsinh(m, l);
     988   [ #  #  #  # ]:          0 :       switch (codeb)
     989                 :            :       {
     990                 :          0 :         case f_YVSLO: gel(T,2) = exptab(gel(T,1), gel(b,2), prec); return T;
     991                 :          0 :         case f_YFAST: gel(T,2) = homtab(initexpexp(m, l), kmb); return T;
     992                 :            :         case f_YOSCS: case f_YOSCC:
     993                 :          0 :           gel(T,2) = homtab(initnumsine(m, l), kmb); return T;
     994                 :            :       }
     995                 :            :     case f_YVSLO:
     996                 :          0 :       tmp = initexpsinh(m, l);
     997                 :          0 :       gel(T,1) = exptab(tmp, gel(a,2), prec);
     998   [ #  #  #  # ]:          0 :       switch (codeb)
     999                 :            :       {
    1000                 :          0 :         case f_YVSLO: gel(T,2) = exptab(tmp, gel(b,2), prec); return T;
    1001                 :          0 :         case f_YFAST: gel(T,2) = homtab(initexpexp(m, l), kmb); return T;
    1002                 :            :         case f_YOSCS:
    1003                 :          0 :         case f_YOSCC: gel(T,2) = homtab(initnumsine(m, l), kmb); return T;
    1004                 :            :       }
    1005                 :            :     case f_YFAST:
    1006                 :          0 :       tmp = initexpexp(m, l);
    1007                 :          0 :       gel(T,1) = homtab(tmp, kma);
    1008      [ #  #  # ]:          0 :       switch (codeb)
    1009                 :            :       {
    1010                 :          0 :         case f_YFAST: gel(T,2) = homtab(tmp, kmb); return T;
    1011                 :            :         case f_YOSCS:
    1012                 :          0 :         case f_YOSCC: gel(T,2) = homtab(initnumsine(m, l), kmb); return T;
    1013                 :            :       }
    1014                 :         28 :     case f_YOSCS: case f_YOSCC: tmp = initnumsine(m, l);
    1015                 :         28 :       gel(T,1) = homtab(tmp, kma);
    1016 [ +  - ][ +  - ]:         28 :       if (codea == f_YOSCC && codeb == f_YOSCC && !gequal(kma, kmb))
                 [ -  + ]
    1017                 :            :       {
    1018                 :          0 :         U = cgetg(3, t_VEC);
    1019                 :          0 :         gel(U,1) = inittanhsinh(m, l);
    1020                 :          0 :         gel(U,2) = homtab(tmp, kmb);
    1021                 :          0 :         gel(T,2) = U;
    1022                 :            :       }
    1023                 :         28 :       else gel(T,2) = homtab(tmp, kmb);
    1024                 :         28 :       return T;
    1025                 :            :   }
    1026                 :        252 :   return gen_0; /* not reached */
    1027                 :            : }
    1028                 :            : 
    1029                 :            : GEN
    1030                 :       9541 : intnuminit0(GEN a, GEN b, GEN tab, long prec)
    1031                 :            : {
    1032                 :            :   long m;
    1033         [ +  + ]:       9541 :   if (!tab) m = 0;
    1034         [ +  + ]:       9352 :   else if (typ(tab) != t_INT)
    1035                 :            :   {
    1036         [ -  + ]:       9345 :     if (!checktab(tab)) pari_err_TYPE("intnuminit0",tab);
    1037                 :       9345 :     return tab;
    1038                 :            :   }
    1039                 :            :   else
    1040                 :          7 :     m = itos(tab);
    1041                 :       9541 :   return intnuminit(a, b, m, prec);
    1042                 :            : }
    1043                 :            : GEN
    1044                 :         49 : sumnuminit0(GEN a, GEN tab, long sgn, long prec)
    1045                 :            : {
    1046                 :            :   long m;
    1047         [ +  + ]:         49 :   if (!tab) m = 0;
    1048         [ +  + ]:         42 :   else if (typ(tab) != t_INT)
    1049                 :            :   {
    1050         [ -  + ]:         35 :     if (!checktab(tab)) pari_err_TYPE("sumnuminit0",tab);
    1051                 :         35 :     return tab;
    1052                 :            :   }
    1053                 :            :   else
    1054                 :          7 :     m = itos(tab);
    1055                 :         49 :   return sumnuminit(a, m, sgn, prec);
    1056                 :            : }
    1057                 :            : 
    1058                 :            : /* User-defined change of variable phi(t) = f(t), where t goes from -oo to +oo,
    1059                 :            :  * and a and b are as in intnuminit. If [a,b] compact, assume phi(t) odd. */
    1060                 :            : static int
    1061                 :      10038 : condfin(long code, GEN xw, long eps, long m, long k)
    1062                 :            : {
    1063                 :            :   GEN x, w;
    1064                 :      10038 :   eps -= 8; /* for safety. Lose 8 bits, but took 1 whole word extra. */
    1065                 :      10038 :   x = gel(xw,1);
    1066                 :      10038 :   w = gel(xw,2);
    1067   [ +  +  +  -  :      10038 :   switch(labs(code))
                      - ]
    1068                 :            :   {
    1069                 :            :     case f_REG: case f_SING:
    1070                 :         14 :       return gexpo(w) < -eps;
    1071                 :            :     case f_YSLOW: case f_YVSLO:
    1072                 :       5012 :       return gexpo(w) - 2*gexpo(x) < -eps;
    1073                 :            :     case f_YFAST:
    1074                 :       5012 :       return cmprs(x, (long)(LOG2 * (gexpo(w) + eps) + 1)) > 0;
    1075                 :            :     case f_YOSCS: case f_YOSCC:
    1076                 :          0 :       return gexpo(x) + m + expu(10 * k) < - eps;
    1077                 :      10038 :     default: return 0;
    1078                 :            :   }
    1079                 :            : }
    1080                 :            : 
    1081                 :            : /* eps = 2^(-k). Return f'(a) ~ (f(a+eps) - f(a-eps)) / 2eps*/
    1082                 :            : static GEN
    1083                 :      20062 : myderiv_num(void *E, GEN (*eval)(void*, GEN), GEN a, GEN eps, long k, long prec)
    1084                 :            : {
    1085                 :      20062 :   GEN d = gmul2n(gsub(eval(E, gadd(a,eps)), eval(E, gsub(a,eps))), k-1);
    1086                 :      20062 :   return gprec_w(d, prec);
    1087                 :            : }
    1088                 :            : /* zp = z to a higher accuracy (enough to evaluate numerical derivative) */
    1089                 :            : static GEN
    1090                 :      20062 : ffprime(void *E, GEN (*eval)(void*, GEN), GEN z, GEN zp, GEN eps, long h, long precl)
    1091                 :            : {
    1092                 :      20062 :   GEN f = eval(E, z), fp = myderiv_num(E, eval, zp, eps, h, precl);
    1093                 :      20062 :   return mkvec2(f, fp);
    1094                 :            : }
    1095                 :            : /* v = [f(z), f'(z)]. Return h := z/(1-f(z)), h + h^2 zf'(z) */
    1096                 :            : static GEN
    1097                 :          0 : ffmodify(GEN v, GEN z)
    1098                 :            : {
    1099                 :          0 :   GEN f = gel(v,1), fp = gel(v,2), h = ginv(gsubsg(1, f));
    1100                 :          0 :   return mkvec2(gmul(z, h), gadd(h, gmul(gsqr(h), gmul(z,fp))));
    1101                 :            : }
    1102                 :            : GEN
    1103                 :         14 : intnuminitgen(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, long m,
    1104                 :            :               long flext, long prec)
    1105                 :            : {
    1106                 :            :   enum {
    1107                 :            :     f_COMP, /* [a,b] */
    1108                 :            :     f_SEMI, /* [a,+-oo[, no oscillation */
    1109                 :            :     f_OSC1, /* [a,+-oo[, oscillation */
    1110                 :            :     f_INF , /* ]-oo,+oo[, no oscillation */
    1111                 :            :     f_OSC2  /* ]-oo,+oo[, oscillation */
    1112                 :            :   };
    1113                 :         14 :   pari_sp ltop = avma;
    1114                 :            :   GEN hnpr, eps, v;
    1115                 :         14 :   long k, h, newprec, lim, precl = prec+EXTRAPREC;
    1116                 :         14 :   long flag, codea = transcode(a, "a"), codeb = transcode(b, "b");
    1117                 :            :   int NOT_OSC, NOT_ODD;
    1118                 :         14 :   intdata D; intinit_start(&D, m, flext, precl);
    1119                 :            : 
    1120                 :         14 :   flag = f_SEMI;
    1121 [ +  - ][ -  + ]:         14 :   if (is_osc_f(codea) || is_osc_f(codeb)) flag = f_OSC1;
    1122 [ +  - ][ -  + ]:         14 :   if (is_fin_f(codea) && is_fin_f(codeb)) flag = f_COMP;
    1123 [ -  + ][ #  # ]:         14 :   else if (!is_fin_f(codea) && !is_fin_f(codeb))
    1124                 :            :   {
    1125         [ #  # ]:          0 :     if (codea * codeb > 0) return gen_0;
    1126         [ #  # ]:          0 :     if (codea != -codeb)
    1127                 :          0 :       pari_err_TYPE("intnuminitgen [infinities of different types]",
    1128                 :            :                     mkvec2(a,b));
    1129         [ #  # ]:          0 :     flag = (flag == f_SEMI) ? f_INF : f_OSC2;
    1130                 :            :   }
    1131 [ +  - ][ -  + ]:         14 :   NOT_OSC = (flag == f_COMP || flag == f_SEMI || flag == f_INF);
                 [ #  # ]
    1132 [ -  + ][ #  # ]:         14 :   NOT_ODD = (flag == f_SEMI || flag == f_OSC1);
    1133                 :            : 
    1134                 :         14 :   newprec = (3*precl - 1)>>1;
    1135                 :         14 :   h = prec2nbits(precl)/2;
    1136                 :         14 :   eps = real2n(-h, newprec);
    1137                 :            : 
    1138 [ -  + ][ #  # ]:         14 :   if (NOT_OSC || !gequal1(eval(E, gen_0)))
    1139                 :         14 :   {
    1140                 :         14 :     GEN a0 = real_0(precl), a0n = real_0(newprec), xw, xwmod;
    1141                 :         14 :     xw = ffprime(E, eval, a0, a0n, eps, h, precl);
    1142         [ -  + ]:         14 :     xwmod = NOT_OSC? xw: ffmodify(xw, a0);
    1143                 :         14 :     D.tabx0 = gel(xwmod,1);
    1144                 :         14 :     D.tabw0 = gel(xwmod,2);
    1145                 :            :   }
    1146                 :            :   else
    1147                 :            :   {
    1148                 :          0 :     GEN xw = gdiv(pol_x(0), gsubsg(1, eval(E, gadd(pol_x(0), zeroser(0, 4)))));
    1149                 :          0 :     D.tabx0 = gprec_w(polcoeff0(xw, 0, 0), precl);
    1150                 :          0 :     D.tabw0 = gprec_w(polcoeff0(xw, 1, 0), precl);
    1151                 :            :   }
    1152                 :            :   /* precl <= newprec */
    1153                 :         14 :   hnpr = real2n(-D.m, newprec);
    1154                 :         14 :   lim = lg(D.tabxp) - 1;
    1155         [ +  - ]:      10024 :   for (k = 1; k <= lim; k++)
    1156                 :            :   {
    1157                 :      10024 :     GEN akn = mulur(k, hnpr), ak = rtor(akn, precl), xw, xwmod;
    1158                 :            :     int finb;
    1159                 :      10024 :     xw = ffprime(E, eval, ak, akn, eps, h, precl);
    1160         [ -  + ]:      10024 :     xwmod = NOT_OSC? xw: ffmodify(xw, ak);
    1161                 :      10024 :     gel(D.tabxp,k) = gel(xwmod,1);
    1162                 :      10024 :     gel(D.tabwp,k) = gel(xwmod,2);
    1163         [ -  + ]:      10024 :     finb = condfin(codeb, is_osc_f(codeb)? xw: xwmod, D.eps, D.m, k);
    1164         [ +  - ]:      10024 :     if (NOT_ODD)
    1165                 :            :     {
    1166                 :      10024 :       ak = negr(ak); akn = negr(akn);
    1167                 :      10024 :       xw = ffprime(E, eval, ak, akn, eps, h, precl);
    1168         [ -  + ]:      10024 :       xwmod = NOT_OSC? xw: ffmodify(xw, ak);
    1169                 :      10024 :       gel(D.tabxm,k) = gel(xwmod,1);
    1170                 :      10024 :       gel(D.tabwm,k) = gel(xwmod,2);
    1171 [ +  + ][ -  + ]:      10024 :       if (finb && condfin(codea, is_osc_f(codeb)? xw: xwmod, D.eps, D.m, k))
                 [ +  - ]
    1172                 :         14 :         break;
    1173                 :            :     }
    1174         [ #  # ]:          0 :     else if (finb) break;
    1175                 :            :   }
    1176         [ +  - ]:         14 :   v = intinit_end(&D, k-1, NOT_ODD? k-1: 0);
    1177         [ -  + ]:         14 :   if (!NOT_OSC)
    1178                 :            :   {
    1179                 :          0 :     GEN C = Pi2n(D.m, precl);
    1180                 :          0 :     TABx0(v) = gmul(TABx0(v), C);
    1181                 :          0 :     TABw0(v) = gmul(TABw0(v), C);
    1182                 :          0 :     TABxp(v) = RgV_Rg_mul(TABxp(v), C);
    1183                 :          0 :     TABwp(v) = RgV_Rg_mul(TABwp(v), C);
    1184                 :          0 :     TABxm(v) = RgV_Rg_mul(TABxm(v), C);
    1185                 :          0 :     TABwm(v) = RgV_Rg_mul(TABwm(v), C);
    1186                 :            :   }
    1187                 :         14 :   return gerepilecopy(ltop, v);
    1188                 :            : }
    1189                 :            : 
    1190                 :            : /* Assigns the values of the function weighted by w[k] at quadrature points x[k]
    1191                 :            :  * [replacing the weights]. Return the index of the last non-zero coeff */
    1192                 :            : static long
    1193                 :          7 : weight(void *E, GEN (*eval)(void *, GEN), GEN x, GEN w)
    1194                 :            : {
    1195                 :          7 :   long k, l = lg(x);
    1196         [ +  + ]:       5026 :   for (k = 1; k < l; k++) gel(w,k) = gmul(gel(w,k), eval(E, gel(x,k)));
    1197 [ +  - ][ +  - ]:          7 :   k--; while (k >= 1) if (!gequal0(gel(w,k--))) break;
    1198                 :          7 :   return k;
    1199                 :            : }
    1200                 :            : /* compute the necessary tabs, weights multiplied by f(t).
    1201                 :            :  * If flag set, assumes that f(-t) = conj(f(t)). */
    1202                 :            : static GEN
    1203                 :          7 : intfuncinitintern(void *E, GEN (*eval)(void*, GEN), GEN tab, long flag)
    1204                 :            : {
    1205                 :          7 :   GEN tabxp = TABxp(tab), tabwp = TABwp(tab);
    1206                 :          7 :   GEN tabxm = TABxm(tab), tabwm = TABwm(tab);
    1207                 :          7 :   long L = weight(E, eval, tabxp, tabwp), L0 = lg(tabxp);
    1208                 :            : 
    1209                 :          7 :   TABw0(tab) = gmul(TABw0(tab), eval(E, TABx0(tab)));
    1210         [ -  + ]:          7 :   if (lg(tabxm) > 1) (void)weight(E, eval, tabxm, tabwm);
    1211                 :            :   else
    1212                 :            :   {
    1213                 :          7 :     tabxm = gneg(tabxp);
    1214         [ +  - ]:          7 :     if (flag) tabwm = gconj(tabwp);
    1215                 :            :     else
    1216                 :            :     {
    1217                 :            :       long L2;
    1218                 :          0 :       tabwm = leafcopy(tabwp);
    1219                 :          0 :       L2 = weight(E, eval, tabxm, tabwm);
    1220         [ #  # ]:          0 :       if (L > L2) L = L2;
    1221                 :            :     }
    1222                 :          7 :     TABxm(tab) = tabxm;
    1223                 :          7 :     TABwm(tab) = tabwm;
    1224                 :            :   }
    1225         [ +  - ]:          7 :   if (L < L0)
    1226                 :            :   { /* catch up functions whose growth at oo was not adequately described */
    1227                 :          7 :     setlg(tabxp, L+1);
    1228                 :          7 :     setlg(tabwp, L+1);
    1229         [ +  - ]:          7 :     if (lg(tabxm) > 1) { setlg(tabxm, L+1); setlg(tabwm, L+1); }
    1230                 :            :   }
    1231                 :          7 :   return tab;
    1232                 :            : }
    1233                 :            : 
    1234                 :            : GEN
    1235                 :          7 : intfuncinit(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, long m, long flag, long prec)
    1236                 :            : {
    1237                 :          7 :   pari_sp ltop = avma;
    1238                 :          7 :   GEN T, tab = intnuminit(a, b, m, prec);
    1239                 :            : 
    1240         [ +  - ]:          7 :   if (lg(tab) != 3) T = intfuncinitintern(E, eval, tab, flag);
    1241                 :            :   else
    1242                 :            :   {
    1243                 :          0 :     T = cgetg(3, t_VEC);
    1244                 :          0 :     gel(T,1) = intfuncinitintern(E, eval, gel(tab,1), flag);
    1245                 :          0 :     gel(T,2) = intfuncinitintern(E, eval, gel(tab,2), flag);
    1246                 :            :   }
    1247                 :          7 :   return gerepilecopy(ltop, T);
    1248                 :            : }
    1249                 :            : 
    1250                 :            : static GEN
    1251                 :       9590 : intnum_i(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, GEN tab, long prec)
    1252                 :            : {
    1253                 :       9590 :   GEN tmp, S = gen_0, res1, res2, tm, pi2, pi2p, pis2, pis2p, kma, kmb;
    1254                 :            :   GEN SP, SN;
    1255                 :       9590 :   long tmpi, sgns = 1, codea = transcode(a, "a"), codeb = transcode(b, "b");
    1256                 :            : 
    1257 [ +  + ][ -  + ]:       9590 :   if (codea == f_REG && typ(a) == t_VEC) a = gel(a,1);
    1258 [ +  + ][ -  + ]:       9590 :   if (codeb == f_REG && typ(b) == t_VEC) b = gel(b,1);
    1259 [ +  + ][ +  + ]:       9590 :   if (codea == f_REG && codeb == f_REG) return intn(E, eval, a, b, tab);
    1260         [ +  + ]:        224 :   if (labs(codea) > labs(codeb)) { swap(a, b); lswap(codea, codeb); sgns = -1; }
    1261                 :            :   /* now labs(codea) <= labs(codeb) */
    1262         [ +  + ]:        224 :   if (codeb == f_SING)
    1263                 :            :   {
    1264         [ +  - ]:          7 :     if (codea == f_REG)
    1265                 :          7 :       S = intnsing(E, eval, b, a, tab, prec), sgns = -sgns;
    1266                 :            :     else
    1267                 :            :     {
    1268                 :          0 :       tmp = gmul2n(gadd(gel(a,1), gel(b,1)), -1);
    1269                 :          0 :       res1 = intnsing(E, eval, a, tmp, tab, prec);
    1270                 :          0 :       res2 = intnsing(E, eval, b, tmp, tab, prec);
    1271                 :          0 :       S = gsub(res1, res2);
    1272                 :            :     }
    1273         [ -  + ]:          7 :     return (sgns < 0) ? gneg(S) : S;
    1274                 :            :   }
    1275                 :            :   /* now b is infinite */
    1276         [ +  - ]:        217 :   tmpi = codeb > 0 ? 1 : -1;
    1277 [ +  + ][ +  + ]:        217 :   if (codea == f_REG && labs(codeb) != f_YOSCC
    1278 [ +  + ][ +  - ]:        140 :       && (labs(codeb) != f_YOSCS || gequal0(a)))
    1279                 :            :   {
    1280                 :        140 :     S = intninfpm(E, eval, a, tmpi, tab);
    1281         [ -  + ]:        140 :     return sgns*tmpi < 0 ? gneg(S) : S;
    1282                 :            :   }
    1283                 :         77 :   pi2 = Pi2n(1, prec); pis2 = Pi2n(-1, prec);
    1284         [ +  + ]:         77 :   if (is_fin_f(codea))
    1285                 :            :   { /* either codea == f_SING  or codea == f_REG and codeb = f_YOSCC
    1286                 :            :      * or (codeb == f_YOSCS and !gequal0(a)) */
    1287                 :          7 :     pi2p = gmul(pi2, f_getycplx(b, prec));
    1288                 :          7 :     pis2p = gmul2n(pi2p, -2);
    1289         [ -  + ]:          7 :     tm = real_i(codea == f_SING ? gel(a,1) : a);
    1290         [ +  - ]:          7 :     if (labs(codeb) == f_YOSCC) tm = gadd(tm, pis2p);
    1291                 :          7 :     tm = gdiv(tm, pi2p);
    1292         [ +  - ]:          7 :     if (tmpi > 0)
    1293                 :          7 :       tm = addsi(1, gceil(tm));
    1294                 :            :     else
    1295                 :          0 :       tm = subis(gfloor(tm), 1);
    1296                 :          7 :     tm = gmul(pi2p, tm);
    1297         [ +  - ]:          7 :     if (labs(codeb) == f_YOSCC) tm = gsub(tm, pis2p);
    1298                 :          7 :     res1 = codea==f_SING? intnsing(E, eval, a,  tm,  gel(tab,1), prec)
    1299         [ -  + ]:          7 :                         : intn    (E, eval, a,  tm,  gel(tab,1));
    1300                 :          7 :     res2 = intninfpm(E, eval, tm, tmpi,gel(tab,2));
    1301         [ -  + ]:          7 :     if (tmpi < 0) res2 = gneg(res2);
    1302                 :          7 :     res1 = gadd(res1, res2);
    1303         [ -  + ]:          7 :     return sgns < 0 ? gneg(res1) : res1;
    1304                 :            :   }
    1305                 :            :   /* now a and b are infinite */
    1306         [ -  + ]:         70 :   if (codea * codeb > 0)
    1307                 :            :   {
    1308                 :          0 :     pari_warn(warner, "integral from infty to infty or from -infty to -infty");
    1309                 :          0 :     return gen_0;
    1310                 :            :   }
    1311         [ -  + ]:         70 :   if (codea > 0) { lswap(codea, codeb); swap(a, b); sgns = -sgns; }
    1312                 :            :   /* now codea < 0 < codeb */
    1313                 :         70 :   codea = -codea;
    1314                 :         70 :   kma = f_getycplx(a, prec);
    1315                 :         70 :   kmb = f_getycplx(b, prec);
    1316 [ -  + ][ #  # ]:         70 :   if ((codea == f_YSLOW && codeb == f_YSLOW)
    1317 [ +  + ][ +  - ]:         70 :    || (codea == f_YFAST && codeb == f_YFAST && gequal(kma, kmb)))
                 [ +  - ]
    1318                 :         14 :     S = intninfinf(E, eval, tab);
    1319                 :            :   else
    1320                 :            :   {
    1321         [ +  + ]:         56 :     GEN coupea = (codea == f_YOSCC)? gmul(pis2, kma): gen_0;
    1322         [ +  + ]:         56 :     GEN coupeb = (codeb == f_YOSCC)? gmul(pis2, kmb): gen_0;
    1323         [ +  + ]:         56 :     GEN coupe = codea == f_YOSCC ? coupea : coupeb;
    1324                 :         56 :     SN = intninfpm(E, eval, coupe, -1, gel(tab,1));
    1325         [ +  + ]:         56 :     if (codea != f_YOSCC)
    1326                 :         28 :       SP = intninfpm(E, eval, coupeb, 1, gel(tab,2));
    1327                 :            :     else
    1328                 :            :     {
    1329         [ -  + ]:         28 :       if (codeb != f_YOSCC) pari_err_BUG("code error in intnum");
    1330         [ +  - ]:         28 :       if (gequal(kma, kmb))
    1331                 :         28 :         SP = intninfpm(E, eval, coupeb, 1, gel(tab,2));
    1332                 :            :       else
    1333                 :            :       {
    1334                 :          0 :         tab = gel(tab,2);
    1335                 :          0 :         SP = intninfpm(E, eval, coupeb, 1, gel(tab,2));
    1336                 :          0 :         SP = gadd(SP, intn(E, eval, coupea, coupeb, gel(tab,1)));
    1337                 :            :       }
    1338                 :            :     }
    1339                 :         56 :     S = gadd(SN, SP);
    1340                 :            :   }
    1341         [ -  + ]:         70 :   if (sgns < 0) S = gneg(S);
    1342                 :       9590 :   return S;
    1343                 :            : }
    1344                 :            : 
    1345                 :            : GEN
    1346                 :       9513 : intnum(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, GEN tab, long prec)
    1347                 :            : {
    1348                 :       9513 :   pari_sp ltop = avma;
    1349                 :       9513 :   long l = prec+EXTRAPREC;
    1350                 :            :   GEN S;
    1351                 :            : 
    1352                 :       9513 :   tab = intnuminit0(a, b, tab, prec);
    1353                 :       9513 :   S = intnum_i(E, eval, gprec_w(a, l), gprec_w(b, l), tab, l);
    1354                 :       9513 :   return gerepilecopy(ltop, gprec_wtrunc(S, prec));
    1355                 :            : }
    1356                 :            : 
    1357                 :            : typedef struct auxint_s {
    1358                 :            :   GEN a, R, pi;
    1359                 :            :   GEN (*f)(void*, GEN);
    1360                 :            :   GEN (*w)(GEN, long);
    1361                 :            :   long prec;
    1362                 :            :   void *E;
    1363                 :            : } auxint_t;
    1364                 :            : 
    1365                 :            : static GEN
    1366                 :       9317 : auxcirc(void *E, GEN t)
    1367                 :            : {
    1368                 :       9317 :   auxint_t *D = (auxint_t*) E;
    1369                 :            :   GEN s, c, z;
    1370                 :       9317 :   mpsincos(mulrr(t, D->pi), &s, &c); z = mkcomplex(c,s);
    1371                 :       9317 :   return gmul(z, D->f(D->E, gadd(D->a, gmul(D->R, z))));
    1372                 :            : }
    1373                 :            : 
    1374                 :            : GEN
    1375                 :          7 : intcirc(void *E, GEN (*eval)(void*, GEN), GEN a, GEN R, GEN tab, long prec)
    1376                 :            : {
    1377                 :            :   auxint_t D;
    1378                 :            :   GEN z;
    1379                 :            : 
    1380                 :          7 :   D.a = a;
    1381                 :          7 :   D.R = R;
    1382                 :          7 :   D.pi = mppi(prec);
    1383                 :          7 :   D.f = eval;
    1384                 :          7 :   D.E = E;
    1385                 :          7 :   z = intnum(&D, &auxcirc, real_m1(prec), real_1(prec), tab, prec);
    1386                 :          7 :   return gmul2n(gmul(R, z), -1);
    1387                 :            : }
    1388                 :            : 
    1389                 :            : static GEN
    1390                 :         42 : gettmpP(GEN x) { return mkvec2(mkvec(gen_1), x); }
    1391                 :            : 
    1392                 :            : static GEN
    1393                 :         70 : gettmpN(GEN tmpP) { return mkvec2(gneg(gel(tmpP,1)), gel(tmpP,2)); }
    1394                 :            : 
    1395                 :            : /* w(Rt) f(a+it) */
    1396                 :            : static GEN
    1397                 :     236789 : auxinv(void *E, GEN t)
    1398                 :            : {
    1399                 :     236789 :   auxint_t *D = (auxint_t*) E;
    1400                 :     236789 :   GEN tmp = D->w(gmul(D->R, t), D->prec);
    1401                 :     236789 :   return gmul(tmp, D->f(D->E, gadd(D->a, mulcxI(t))));
    1402                 :            : }
    1403                 :            : static GEN
    1404                 :         35 : intinvintern(void *E, GEN (*eval)(void*, GEN), GEN sig, GEN x, GEN tab, long prec)
    1405                 :            : {
    1406                 :            :   auxint_t D;
    1407                 :            :   GEN z, zR, zI, tmpP, tmpN;
    1408                 :            : 
    1409 [ +  - ][ +  - ]:         35 :   if (lg(sig) != 3 || !isinR(gel(sig,1)) || !isinR(gel(sig,2)))
                 [ -  + ]
    1410                 :          0 :     pari_err_TYPE("integral transform",sig);
    1411         [ -  + ]:         35 :   if (gsigne(gel(sig,2)) < 0)
    1412                 :          0 :     pari_err_OVERFLOW("integral transform [exponential increase]");
    1413                 :         35 :   D.a = gel(sig,1);
    1414                 :         35 :   D.prec = prec;
    1415                 :         35 :   D.f = eval;
    1416                 :         35 :   D.E = E;
    1417         [ +  + ]:         35 :   if (gequal0(gel(sig,2)))
    1418                 :            :   {
    1419                 :         28 :     D.R = x;
    1420                 :         28 :     tmpP = gettmpP(mulcxI(gabs(x, prec)));
    1421                 :         28 :     tmpN = gettmpN(tmpP);
    1422                 :         28 :     tab = intnuminit0(tmpN, tmpP, tab, prec);
    1423                 :         28 :     D.w = gcos;
    1424                 :         28 :     zR = intnum_i(&D, &auxinv, tmpN, tmpP, tab, prec);
    1425                 :         28 :     gel(tmpP,2) = gneg(gel(tmpP,2));
    1426                 :         28 :     D.w = gsin;
    1427                 :         28 :     zI = intnum_i(&D, &auxinv, gettmpN(tmpP), tmpP, tab, prec);
    1428                 :         28 :     z = gadd(zR, mulcxI(zI));
    1429                 :            :   }
    1430                 :            :   else
    1431                 :            :   {
    1432                 :          7 :     D.R = mulcxI(x);
    1433                 :          7 :     tmpP = gettmpP(gel(sig,2));
    1434                 :          7 :     D.w = gexp;
    1435                 :          7 :     z = intnum(&D, &auxinv, gettmpN(tmpP), tmpP, tab, prec);
    1436                 :            :   }
    1437                 :         35 :   return gdiv(gmul(gexp(gmul(gel(sig,1), x), prec), z), Pi2n(1, prec));
    1438                 :            : }
    1439                 :            : 
    1440                 :            : /* If sig = [sigR, e]: if e = 0, slowly decreasing, if e > 0, exponentially
    1441                 :            :  * decreasing like exp(-e*t). If sig is real, identical to [sig, 1]. */
    1442                 :            : GEN
    1443                 :          7 : intmellininv(void *E, GEN (*eval)(void*, GEN), GEN sig, GEN x, GEN tab, long prec)
    1444                 :            : {
    1445         [ +  - ]:          7 :   if (typ(sig) != t_VEC) sig = mkvec2(sig, gen_1);
    1446                 :          7 :   return intinvintern(E, eval, sig, gneg(glog(x, prec)), tab, prec);
    1447                 :            : }
    1448                 :            : 
    1449                 :            : /* If sig = [sigR, e]: if e = 0, slowly decreasing, if e > 0, exponentially
    1450                 :            :  * decreasing like exp(-e*t). If sig is real, identical to [sig, 0]. */
    1451                 :            : GEN
    1452                 :         28 : intlaplaceinv(void *E, GEN (*eval)(void*, GEN), GEN sig, GEN x, GEN tab, long prec)
    1453                 :            : {
    1454         [ +  - ]:         28 :   if (typ(sig) != t_VEC) sig = mkvec2(sig, gen_0);
    1455                 :         28 :   return intinvintern(E, eval, sig, x, tab, prec);
    1456                 :            : }
    1457                 :            : 
    1458                 :            : /* assume tab computed with additional weights f(sig + I*T) */
    1459                 :            : typedef struct auxmel_s {
    1460                 :            :   GEN L;
    1461                 :            :   long prec;
    1462                 :            : } auxmel_t;
    1463                 :            : 
    1464                 :            : static GEN
    1465                 :       5019 : auxmelshort(void *E, GEN t)
    1466                 :            : {
    1467                 :       5019 :   auxmel_t *D = (auxmel_t*) E;
    1468                 :       5019 :   return gexp(gmul(D->L, t), D->prec);
    1469                 :            : }
    1470                 :            : 
    1471                 :            : GEN
    1472                 :          7 : intmellininvshort(GEN sig, GEN x, GEN tab, long prec)
    1473                 :            : {
    1474                 :            :   auxmel_t D;
    1475                 :          7 :   GEN z, tmpP, LX = gneg(glog(x, prec));
    1476                 :            : 
    1477         [ +  - ]:          7 :   if (typ(sig) != t_VEC) sig = mkvec2(sig, gen_1);
    1478 [ +  - ][ +  - ]:          7 :   if (lg(sig) != 3 || !isinR(gel(sig,1)) || !isinR(gel(sig,2)))
                 [ -  + ]
    1479                 :          0 :     pari_err_TYPE("intmellininvshort",sig);
    1480         [ -  + ]:          7 :   if (gsigne(gel(sig,2)) <= 0)
    1481                 :          0 :     pari_err_OVERFLOW("intinvmellinshort [need exponential decrease]");
    1482                 :          7 :   D.L = mulcxI(LX);
    1483                 :          7 :   D.prec = prec;
    1484                 :          7 :   tmpP = gettmpP(gel(sig,2));
    1485                 :          7 :   z = intnum_i(&D, &auxmelshort, gettmpN(tmpP), tmpP, tab, prec);
    1486                 :          7 :   return gdiv(gmul(gexp(gmul(gel(sig,1), LX), prec), z), Pi2n(1, prec));
    1487                 :            : }
    1488                 :            : 
    1489                 :            : /* a as in intnum. flag = 0 for sin, flag = 1 for cos. */
    1490                 :            : static GEN
    1491                 :         56 : mytra(GEN a, GEN x, long flag, const char *name)
    1492                 :            : {
    1493                 :            :   GEN b, xa;
    1494                 :         56 :   long s, codea = transcode(a, name);
    1495                 :            : 
    1496   [ +  -  -  - ]:         56 :   switch (labs(codea))
    1497                 :            :   {
    1498                 :         56 :     case f_REG: case f_SING: case f_YFAST: return a;
    1499                 :            :     case f_YSLOW: case f_YVSLO:
    1500                 :          0 :       xa = real_i(x); s = gsigne(xa);
    1501         [ #  # ]:          0 :       if (!s) pari_err_DOMAIN("Fourier transform","Re(x)","=",gen_0,x);
    1502         [ #  # ]:          0 :       if (s < 0) xa = gneg(xa);
    1503                 :          0 :       b = cgetg(3, t_VEC);
    1504         [ #  # ]:          0 :       gel(b,1) = mkvec( codea > 0 ? gen_1 : gen_m1 );
    1505         [ #  # ]:          0 :       gel(b,2) = (flag? mulcxI(xa): mulcxmI(xa));
    1506                 :          0 :       return b;
    1507                 :            :     case f_YOSCS: case f_YOSCC:
    1508                 :          0 :       pari_err_IMPL("Fourier transform of oscillating functions");
    1509                 :            :   }
    1510                 :         56 :   return NULL;
    1511                 :            : }
    1512                 :            : 
    1513                 :            : /* w(ta) f(t) */
    1514                 :            : static GEN
    1515                 :      40236 : auxfour(void *E, GEN t)
    1516                 :            : {
    1517                 :      40236 :   auxint_t *D = (auxint_t*) E;
    1518                 :      40236 :   return gmul(D->w(gmul(t, D->a), D->prec), D->f(D->E, t));
    1519                 :            : }
    1520                 :            : 
    1521                 :            : GEN
    1522                 :         14 : intfouriersin(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, GEN x, GEN tab, long prec)
    1523                 :            : {
    1524                 :            :   auxint_t D;
    1525                 :            :   GEN tmp;
    1526                 :            : 
    1527         [ -  + ]:         14 :   if (gequal0(x)) return gcopy(x);
    1528                 :         14 :   tmp = gmul(x, Pi2n(1, prec));
    1529                 :         14 :   D.a = tmp;
    1530                 :         14 :   D.R = NULL;
    1531                 :         14 :   D.prec = prec;
    1532                 :         14 :   D.f = eval;
    1533                 :         14 :   D.E = E;
    1534                 :         14 :   a = mytra(a,tmp,0,"a");
    1535                 :         14 :   b = mytra(b,tmp,0,"b");
    1536                 :         14 :   D.w = gsin;
    1537                 :         14 :   return intnum(&D, &auxfour, a, b, tab, prec);
    1538                 :            : }
    1539                 :            : 
    1540                 :            : GEN
    1541                 :         14 : intfouriercos(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, GEN x, GEN tab, long prec)
    1542                 :            : {
    1543                 :            :   auxint_t D;
    1544                 :            :   GEN tmp;
    1545                 :            : 
    1546         [ -  + ]:         14 :   if (gequal0(x)) return intnum(E, eval, a, b, tab, prec);
    1547                 :         14 :   tmp = gmul(x, Pi2n(1, prec));
    1548                 :         14 :   D.a = tmp;
    1549                 :         14 :   D.R = NULL;
    1550                 :         14 :   D.prec = prec;
    1551                 :         14 :   D.f = eval;
    1552                 :         14 :   D.E = E;
    1553                 :         14 :   a = mytra(a,tmp,1,"a");
    1554                 :         14 :   b = mytra(b,tmp,1,"b");
    1555                 :         14 :   D.w = gcos;
    1556                 :         14 :   return intnum(&D, &auxfour, a, b, tab, prec);
    1557                 :            : }
    1558                 :            : 
    1559                 :            : GEN
    1560                 :          7 : intfourierexp(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, GEN x, GEN tab,
    1561                 :            :               long prec)
    1562                 :            : {
    1563                 :          7 :   pari_sp ltop = avma;
    1564                 :          7 :   GEN R = intfouriercos(E, eval, a, b, x, tab, prec);
    1565                 :          7 :   GEN I = intfouriersin(E, eval, a, b, x, tab, prec);
    1566                 :          7 :   return gerepileupto(ltop, gadd(R, mulcxmI(I)));
    1567                 :            : }
    1568                 :            : 
    1569                 :            : GEN
    1570                 :         28 : intnumromb(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, long flag, long prec)
    1571                 :            : {
    1572                 :         28 :   pari_sp av = avma;
    1573                 :            :   GEN z;
    1574   [ +  +  +  +  :         28 :   switch(flag)
                      - ]
    1575                 :            :   {
    1576                 :          7 :     case 0: z = qrom3  (E, eval, a, b, prec); break;
    1577                 :          7 :     case 1: z = rombint(E, eval, a, b, prec); break;
    1578                 :          7 :     case 2: z = qromi  (E, eval, a, b, prec); break;
    1579                 :          7 :     case 3: z = qrom2  (E, eval, a, b, prec); break;
    1580                 :          0 :     default: pari_err_FLAG("intnumromb"); return NULL; /* not reached */
    1581                 :            :   }
    1582         [ -  + ]:         28 :   if (!z) pari_err_IMPL("intnumromb recovery [too many iterations]");
    1583                 :         28 :   return gerepileupto(av, z);
    1584                 :            : }
    1585                 :            : 
    1586                 :            : GEN
    1587                 :         28 : intnumromb0(GEN a, GEN b, GEN code, long flag, long prec)
    1588                 :         28 : { EXPR_WRAP(code, intnumromb(EXPR_ARG, a, b, flag, prec)); }
    1589                 :            : GEN
    1590                 :       9471 : intnum0(GEN a, GEN b, GEN code, GEN tab, long prec)
    1591                 :       9471 : { EXPR_WRAP(code, intnum(EXPR_ARG, a, b, tab, prec)); }
    1592                 :            : GEN
    1593                 :          7 : intcirc0(GEN a, GEN R, GEN code, GEN tab, long prec)
    1594                 :          7 : { EXPR_WRAP(code, intcirc(EXPR_ARG, a, R, tab, prec)); }
    1595                 :            : GEN
    1596                 :          7 : intmellininv0(GEN sig, GEN x, GEN code, GEN tab, long prec)
    1597                 :          7 : { EXPR_WRAP(code, intmellininv(EXPR_ARG, sig, x, tab, prec)); }
    1598                 :            : GEN
    1599                 :         28 : intlaplaceinv0(GEN sig, GEN x, GEN code, GEN tab, long prec)
    1600                 :         28 : { EXPR_WRAP(code, intlaplaceinv(EXPR_ARG, sig, x, tab, prec)); }
    1601                 :            : GEN
    1602                 :          7 : intfourcos0(GEN a, GEN b, GEN x, GEN code, GEN tab, long prec)
    1603                 :          7 : { EXPR_WRAP(code, intfouriercos(EXPR_ARG, a, b, x, tab, prec)); }
    1604                 :            : GEN
    1605                 :          7 : intfoursin0(GEN a, GEN b, GEN x, GEN code, GEN tab, long prec)
    1606                 :          7 : { EXPR_WRAP(code, intfouriersin(EXPR_ARG, a, b, x, tab, prec)); }
    1607                 :            : GEN
    1608                 :          7 : intfourexp0(GEN a, GEN b, GEN x, GEN code, GEN tab, long prec)
    1609                 :          7 : { EXPR_WRAP(code, intfourierexp(EXPR_ARG, a, b, x, tab, prec)); }
    1610                 :            : 
    1611                 :            : GEN
    1612                 :         14 : intnuminitgen0(GEN a, GEN b, GEN code, long m, long flag, long prec)
    1613                 :         14 : { EXPR_WRAP(code, intnuminitgen(EXPR_ARG, a, b, m, flag, prec)); }
    1614                 :            : 
    1615                 :            : /* m and flag reversed on purpose */
    1616                 :            : GEN
    1617                 :          7 : intfuncinit0(GEN a, GEN b, GEN code, long flag, long m, long prec)
    1618                 :          7 : { EXPR_WRAP(code, intfuncinit(EXPR_ARG, a, b, m, flag? 1: 0, prec)); }
    1619                 :            : 
    1620                 :            : #if 0
    1621                 :            : /* Two variable integration */
    1622                 :            : 
    1623                 :            : typedef struct auxf_s {
    1624                 :            :   GEN x;
    1625                 :            :   GEN (*f)(void *, GEN, GEN);
    1626                 :            :   void *E;
    1627                 :            : } auxf_t;
    1628                 :            : 
    1629                 :            : typedef struct indi_s {
    1630                 :            :   GEN (*c)(void*, GEN);
    1631                 :            :   GEN (*d)(void*, GEN);
    1632                 :            :   GEN (*f)(void *, GEN, GEN);
    1633                 :            :   void *Ec;
    1634                 :            :   void *Ed;
    1635                 :            :   void *Ef;
    1636                 :            :   GEN tabintern;
    1637                 :            :   long prec;
    1638                 :            : } indi_t;
    1639                 :            : 
    1640                 :            : static GEN
    1641                 :            : auxf(GEN y, void *E)
    1642                 :            : {
    1643                 :            :   auxf_t *D = (auxf_t*) E;
    1644                 :            :   return D->f(D->E, D->x, y);
    1645                 :            : }
    1646                 :            : 
    1647                 :            : static GEN
    1648                 :            : intnumdoubintern(GEN x, void *E)
    1649                 :            : {
    1650                 :            :   indi_t *D = (indi_t*) E;
    1651                 :            :   GEN c = D->c(x, D->Ec), d = D->d(x, D->Ed);
    1652                 :            :   auxf_t A;
    1653                 :            : 
    1654                 :            :   A.x = x;
    1655                 :            :   A.f = D->f;
    1656                 :            :   A.E = D->Ef;
    1657                 :            :   return intnum(&A, &auxf, c, d, D->tabintern, D->prec);
    1658                 :            : }
    1659                 :            : 
    1660                 :            : GEN
    1661                 :            : 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)
    1662                 :            : {
    1663                 :            :   indi_t E;
    1664                 :            : 
    1665                 :            :   E.c = evalc;
    1666                 :            :   E.d = evald;
    1667                 :            :   E.f = evalf;
    1668                 :            :   E.Ec = Ec;
    1669                 :            :   E.Ed = Ed;
    1670                 :            :   E.Ef = Ef;
    1671                 :            :   E.prec = prec;
    1672                 :            :   if (typ(tabint) == t_INT)
    1673                 :            :   {
    1674                 :            :     GEN C = evalc(a, Ec), D = evald(a, Ed);
    1675                 :            :     if (typ(C) != t_VEC && typ(D) != t_VEC) { C = gen_0; D = gen_1; }
    1676                 :            :     E.tabintern = intnuminit0(C, D, tabint, prec);
    1677                 :            :   }
    1678                 :            :   else E.tabintern = tabint;
    1679                 :            :   return intnum(&E, &intnumdoubintern, a, b, tabext, prec);
    1680                 :            : }
    1681                 :            : 
    1682                 :            : GEN
    1683                 :            : intnumdoub0(GEN a, GEN b, int nc, int nd, int nf, GEN tabext, GEN tabint, long prec)
    1684                 :            : {
    1685                 :            :   GEN z;
    1686                 :            :   push_lex(NULL);
    1687                 :            :   push_lex(NULL);
    1688                 :            :   z = intnumdoub(chf, &gp_eval2, chc, &gp_eval, chd, &gp_eval, a, b, tabext, tabint, prec);
    1689                 :            :   pop_lex(1); pop_lex(1); return z;
    1690                 :            : }
    1691                 :            : #endif
    1692                 :            : 
    1693                 :            : /* Numerical summation routine assuming f holomorphic for Re(s) >= sig.
    1694                 :            :  * Computes sum_{n>=a} f(n)  if sgn >= 0,
    1695                 :            :  *          sum_{n>=a} (-1)^n f(n) otherwise,  where a is real.
    1696                 :            :  * Variant of Abel-Plana. */
    1697                 :            : 
    1698                 :            : static GEN
    1699                 :      52528 : auxsum(void *E, GEN t)
    1700                 :            : {
    1701                 :      52528 :   auxint_t *D = (auxint_t*) E;
    1702                 :      52528 :   GEN z = mkcomplex(D->a, t);
    1703                 :      52528 :   return D->f(D->E, z);
    1704                 :            : }
    1705                 :            : /* assume that conj(f(z)) = f(conj(z)) */
    1706                 :            : static GEN
    1707                 :     257838 : auxsumintern1(void *E, GEN t, long sgn)
    1708                 :            : {
    1709                 :     257838 :   auxint_t *D = (auxint_t*) E;
    1710                 :     257838 :   GEN z = mkcomplex(D->a, t), u = D->f(D->E, z);
    1711         [ +  + ]:     257838 :   return sgn > 0 ? imag_i(u): real_i(u);
    1712                 :            : }
    1713                 :            : /* no assumption */
    1714                 :            : static GEN
    1715                 :          0 : auxsumintern(void *E, GEN t, long sgn)
    1716                 :            : {
    1717                 :          0 :   auxint_t *D = (auxint_t*) E;
    1718                 :          0 :   GEN u,v, z = mkcomplex(D->a, t);
    1719                 :          0 :   u = D->f(D->E, z); gel(z,2) = gneg(t);
    1720         [ #  # ]:          0 :   v = D->f(D->E, z); return sgn > 0 ? gsub(u, v) : gadd(u, v);
    1721                 :            : }
    1722                 :            : static GEN
    1723                 :          0 : auxsum0(void *E, GEN t) { return auxsumintern(E, t, 1); }
    1724                 :            : static GEN
    1725                 :     231567 : auxsum1(void *E, GEN t) { return auxsumintern1(E, t, 1); }
    1726                 :            : static GEN
    1727                 :          0 : auxsumalt0(void *E, GEN t) { return auxsumintern(E, t, -1); }
    1728                 :            : static GEN
    1729                 :      26271 : auxsumalt1(void *E, GEN t) { return auxsumintern1(E, t, -1); }
    1730                 :            : 
    1731                 :            : static GEN
    1732                 :         49 : sumnumall(void *E, GEN (*eval)(void*, GEN), GEN a, GEN sig, GEN tab, long flag, long sgn, long prec)
    1733                 :            : {
    1734                 :            :   GEN SI, S, nsig, b, signew;
    1735                 :         49 :   long si = 1, flii;
    1736                 :         49 :   pari_sp ltop = avma;
    1737                 :            :   auxint_t D;
    1738                 :            : 
    1739                 :         49 :   b = suminit_start(sig);
    1740                 :         49 :   flii = gequal0(gel(b,2));
    1741         [ -  + ]:         49 :   if (!is_scalar_t(typ(a))) pari_err_TYPE("sumnum",a);
    1742                 :         49 :   tab = sumnuminit0(sig, tab, sgn, prec);
    1743                 :            : 
    1744         [ +  + ]:         49 :   signew = (typ(sig) == t_VEC) ? gel(sig,1) : sig;
    1745                 :         49 :   a = gceil(a); nsig = gmax(subis(a, 1), gceil(gsub(signew, ghalf)));
    1746         [ +  + ]:         49 :   if (sgn < 0) {
    1747         [ -  + ]:          7 :     if (mpodd(nsig)) nsig = addsi(1, nsig);
    1748         [ +  - ]:          7 :     si = mpodd(a) ? -1 : 1;
    1749                 :            :   }
    1750                 :         49 :   SI = real_0(prec);
    1751         [ +  + ]:        147 :   while (cmpii(a, nsig) <= 0)
    1752                 :            :   {
    1753         [ +  + ]:         98 :     SI = (si < 0) ? gsub(SI, eval(E, a)) : gadd(SI, eval(E, a));
    1754         [ +  + ]:         98 :     a = addsi(1, a); if (sgn < 0) si = -si;
    1755                 :            :   }
    1756                 :         49 :   D.a = gadd(nsig, ghalf);
    1757                 :         49 :   D.R = gen_0;
    1758                 :         49 :   D.f = eval;
    1759                 :         49 :   D.E = E;
    1760                 :         49 :   D.prec = prec;
    1761         [ +  + ]:         49 :   if (!flii)
    1762 [ +  - ][ +  - ]:         14 :     S = intnum_i(&D, sgn > 0? (flag ? &auxsum1 : &auxsum0)
    1763         [ #  # ]:          0 :                             : (flag ? &auxsumalt1 : &auxsumalt0),
    1764                 :            :                       gen_0, b, tab, prec);
    1765                 :            :   else
    1766                 :            :   {
    1767         [ +  + ]:         35 :     if (flag)
    1768                 :            :     {
    1769                 :         28 :       GEN emp = leafcopy(tab); TABwm(emp) = TABwp(emp);
    1770         [ +  + ]:         28 :       S = gmul2n(intninfinf(&D, sgn > 0? &auxsum1: &auxsumalt1, emp),-1);
    1771                 :            :     }
    1772                 :            :     else
    1773                 :          7 :       S = intninfinfintern(&D, &auxsum, tab, sgn);
    1774                 :            :   }
    1775         [ +  + ]:         49 :   if (flag) S = gneg(S);
    1776                 :            :   else
    1777                 :            :   {
    1778                 :          7 :     S = gmul2n(S, -1);
    1779         [ -  + ]:          7 :     S = (sgn < 0) ? gneg(S): mulcxI(S);
    1780                 :            :   }
    1781                 :         49 :   return gerepileupto(ltop, gadd(SI, S));
    1782                 :            : }
    1783                 :            : GEN
    1784                 :         42 : sumnum(void *E, GEN (*f)(void *, GEN), GEN a,GEN sig,GEN tab,long flag,long prec)
    1785                 :         42 : { return sumnumall(E,f,a,sig,tab,flag,1,prec); }
    1786                 :            : GEN
    1787                 :          7 : sumnumalt(void *E, GEN (*f)(void *, GEN),GEN a,GEN s,GEN tab,long flag,long prec)
    1788                 :          7 : { return sumnumall(E,f,a,s,tab,flag,-1,prec); }
    1789                 :            : 
    1790                 :            : GEN
    1791                 :         42 : sumnum0(GEN a, GEN sig, GEN code, GEN tab, long flag, long prec)
    1792                 :         42 : { EXPR_WRAP(code, sumnum(EXPR_ARG, a, sig, tab, flag, prec)); }
    1793                 :            : GEN
    1794                 :          7 : sumnumalt0(GEN a, GEN sig, GEN code, GEN tab, long flag, long prec)
    1795                 :          7 : { EXPR_WRAP(code, sumnumalt(EXPR_ARG, a, sig, tab, flag, prec)); }

Generated by: LCOV version 1.9