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