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