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