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