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