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