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 : /********************************************************************/
16 : /** **/
17 : /** TRANSCENDENTAL FONCTIONS **/
18 : /** (part 3) **/
19 : /** **/
20 : /********************************************************************/
21 : #include "pari.h"
22 : #include "paripriv.h"
23 :
24 : #define DEBUGLEVEL DEBUGLEVEL_trans
25 :
26 : #define HALF_E 1.3591409 /* exp(1) / 2 */
27 :
28 : /***********************************************************************/
29 : /** **/
30 : /** BESSEL FUNCTIONS **/
31 : /** **/
32 : /***********************************************************************/
33 :
34 : static GEN
35 900443 : _abs(GEN x)
36 900443 : { return gabs(gtofp(x,LOWDEFAULTPREC), LOWDEFAULTPREC); }
37 : /* can we use asymptotic expansion ? */
38 : static int
39 408964 : bessel_asymp(GEN n, GEN z, long bit)
40 : {
41 : GEN Z, N;
42 408964 : long t = typ(n);
43 408964 : if (!is_real_t(t) && t != t_COMPLEX) return 0;
44 408782 : Z = _abs(z); N = gaddgs(_abs(n), 1);
45 408782 : return gcmpgs(gdiv(Z, gsqr(N)), (bit+10)/2) >= 0; }
46 :
47 : /* Region I: 0 < Arg z <= Pi, II: -Pi < Arg z <= 0 */
48 : static int
49 518 : regI(GEN z)
50 : {
51 518 : long s = gsigne(imag_i(z));
52 518 : return (s > 0 || (s == 0 && gsigne(real_i(z)) < 0)) ? 1 : 2;
53 : }
54 : /* Region 1: Re(z) >= 0, 2: Re(z) < 0, Im(z) >= 0, 3: Re(z) < 0, Im(z) < 0 */
55 : static int
56 81899 : regJ(GEN z)
57 : {
58 81899 : if (gsigne(real_i(z)) >= 0) return 1;
59 336 : return gsigne(imag_i(z)) >= 0 ? 2 : 3;
60 : }
61 :
62 : /* Hankel's expansions:
63 : * a_k(n) = \prod_{0 <= j < k} (4n^2 - (2j+1)^2)
64 : * C(k)[n,z] = a_k(n) / (k! (8 z)^k)
65 : * A(z) = exp(-z) sum_{k >= 0} C(k)
66 : * A(-z) = exp(z) sum_{k >= 0} (-1)^k C(k)
67 : * J_n(z) ~ [1] (A(z/i) / r + A(-z/i) r) / sqrt(2Pi z)
68 : * [2] (A(z/i) r^3 + A(-z/i) r) / sqrt(2Pi z)
69 : * [3] (A(z/i) / r + A(-z/i) / r^3) / sqrt(2Pi z)
70 : * Y_n(z) ~ [1] i(A(z/i) / r + A(-z/i) r) / sqrt(2Pi z)
71 : * [2] i(A(z/i) (r^3-2/r) + A(-z/i) r) / sqrt(2Pi z)
72 : * [3] i(-A(z/i)/r + A(-z/i)(2r-1/r^3)) / sqrt(2Pi z)
73 : * K_n(z) ~ A(z) Pi / sqrt(2 Pi z)
74 : * I_n(z) ~ [I] (A(-z) + r^2 A(z)) / sqrt(2 Pi z)
75 : * [II](A(-z) + r^(-2) A(z)) / sqrt(2 Pi z) */
76 :
77 : /* set [A(z), A(-z), exp((2*nu+1)*I*Pi/4)] */
78 : static void
79 82879 : hankel_ABr(GEN *pA, GEN *pB, GEN *pr, GEN n, GEN z, long bit)
80 : {
81 82879 : GEN E, P, C, Q = gen_0, zi = ginv(gmul2n(z, 3));
82 82879 : GEN K = gaddgs(_abs(n), 1), n2 = gmul2n(gsqr(n),2);
83 82879 : long prec = nbits2prec(bit), B = bit + 4, m;
84 :
85 82879 : P = C = real_1_bit(bit);
86 82879 : for (m = 1;; m += 2)
87 : {
88 5475740 : C = gmul(C, gdivgu(gmul(gsub(n2, sqru(2*m - 1)), zi), m));
89 5475740 : Q = gadd(Q, C);
90 5475740 : C = gmul(C, gdivgu(gmul(gsub(n2, sqru(2*m + 1)), zi), m + 1));
91 5475740 : P = gadd(P, C);
92 5475740 : if (gexpo(C) < -B && gcmpgs(K, m) <= 0) break;
93 : }
94 82879 : E = gexp(z, prec);
95 82879 : *pA = gdiv(gadd(P, Q), E);
96 82879 : *pB = gmul(gsub(P, Q), E);
97 82879 : *pr = gexp(mulcxI(gmul(gaddgs(gmul2n(n,1), 1), Pi2n(-2, prec))), prec);
98 82879 : }
99 :
100 : /* sqrt(2*Pi*z) */
101 : static GEN
102 82879 : sqz(GEN z, long bit)
103 : {
104 82879 : long prec = nbits2prec(bit);
105 82879 : return gsqrt(gmul(Pi2n(1, prec), z), prec);
106 : }
107 :
108 : static GEN
109 462 : besskasymp(GEN nu, GEN z, long bit)
110 : {
111 : GEN A, B, r;
112 462 : long prec = nbits2prec(bit);
113 462 : hankel_ABr(&A,&B,&r, nu, z, bit);
114 462 : return gdiv(gmul(A, mppi(prec)), sqz(z, bit));
115 : }
116 :
117 : static GEN
118 518 : bessiasymp(GEN nu, GEN z, long bit)
119 : {
120 : GEN A, B, r, R, r2;
121 518 : hankel_ABr(&A,&B,&r, nu, z, bit);
122 518 : r2 = gsqr(r);
123 518 : R = regI(z) == 1 ? gmul(A, r2) : gdiv(A, r2);
124 518 : return gdiv(gadd(B, R), sqz(z, bit));
125 : }
126 :
127 : static GEN
128 81433 : bessjasymp(GEN nu, GEN z, long bit)
129 : {
130 : GEN A, B, r, R;
131 81433 : long reg = regJ(z);
132 81433 : hankel_ABr(&A,&B,&r, nu, mulcxmI(z), bit);
133 81433 : if (reg == 1) R = gadd(gdiv(A, r), gmul(B, r));
134 168 : else if (reg == 2) R = gadd(gmul(A, gpowgs(r, 3)), gmul(B, r));
135 56 : else R = gadd(gdiv(A, r), gdiv(B, gpowgs(r, 3)));
136 81433 : return gdiv(R, sqz(z, bit));
137 : }
138 :
139 : static GEN
140 466 : bessyasymp(GEN nu, GEN z, long bit)
141 : {
142 : GEN A, B, r, R;
143 466 : long reg = regJ(z);
144 466 : hankel_ABr(&A,&B,&r, nu, mulcxmI(z), bit);
145 466 : if (reg == 1) R = gsub(gmul(B, r), gdiv(A, r));
146 168 : else if (reg == 2)
147 112 : R = gadd(gmul(A, gsub(gpowgs(r, 3), gmul2n(ginv(r), 1))), gmul(B, r));
148 : else
149 56 : R = gsub(gmul(B, gsub(gmul2n(r, 1), ginv(gpowgs(r, 3)))), gdiv(A, r));
150 466 : return gdiv(mulcxI(R), sqz(z, bit));
151 : }
152 :
153 : /* n! sum_{0 <= k <= m} x^k / (k!*(k+n)!) */
154 : static GEN
155 314931 : _jbessel(GEN n, GEN x, long m)
156 : {
157 314931 : pari_sp av = avma;
158 314931 : GEN s = gen_1;
159 : long k;
160 :
161 109563028 : for (k = m; k >= 1; k--)
162 : {
163 109248097 : s = gaddsg(1, gdiv(gmul(x,s), gmulgu(gaddgs(n, k), k)));
164 109248097 : if (gc_needed(av,1))
165 : {
166 0 : if (DEBUGMEM>1) pari_warn(warnmem,"besselj");
167 0 : s = gc_upto(av, s);
168 : }
169 : }
170 314931 : return s;
171 : }
172 :
173 : /* max(2, L * approximate solution to x log x = B) */
174 : static long
175 325504 : bessel_get_lim(double B, double L)
176 325504 : { return maxss(2, L * exp(dbllambertW0(B))); }
177 :
178 : static GEN
179 42 : vjbesselh(void* E, GEN z, long prec){return jbesselh((GEN)E,z,prec);}
180 : static GEN
181 126 : vjbessel(void* E, GEN z, long prec) {return jbessel((GEN)E,z,prec);}
182 : static GEN
183 42 : vibessel(void* E, GEN z, long prec) {return ibessel((GEN)E,z,prec);}
184 : static GEN
185 126 : vnbessel(void* E, GEN z, long prec) {return ybessel((GEN)E,z,prec);}
186 : static GEN
187 42 : vkbessel(void* E, GEN z, long prec) {return kbessel((GEN)E,z,prec);}
188 :
189 : /* if J != 0 BesselJ, else BesselI. */
190 : static GEN
191 397197 : jbesselintern(GEN n, GEN z, long J, long prec)
192 : {
193 397197 : const char *f = J? "besselj": "besseli";
194 : long i, ki;
195 397197 : pari_sp av = avma;
196 : GEN y;
197 :
198 397197 : switch(typ(z))
199 : {
200 396791 : case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
201 : {
202 396791 : int flz0 = gequal0(z);
203 : long lim, k, precnew, bit;
204 : GEN p1, p2;
205 : double az, L;
206 :
207 396791 : i = precision(z); if (i) prec = i;
208 396791 : if (flz0 && gequal0(n)) return real_1(prec);
209 396791 : bit = prec2nbits(prec);
210 396791 : if (bessel_asymp(n, z, bit))
211 : {
212 81951 : GEN R = J? bessjasymp(n, z, bit): bessiasymp(n, z, bit);
213 81951 : if (typ(R) == t_COMPLEX && isexactzero(imag_i(n))
214 81265 : && gsigne(real_i(z)) > 0
215 81111 : && isexactzero(imag_i(z))) R = gcopy(gel(R,1));
216 81951 : return gc_upto(av, R);
217 : }
218 314840 : p2 = gpow(gmul2n(z,-1),n,prec);
219 314812 : p2 = gdiv(p2, ggamma(gaddgs(n,1),prec));
220 314812 : if (flz0) return gc_upto(av, p2);
221 314812 : az = dblmodulus(z); L = HALF_E * az;
222 314812 : precnew = prec;
223 314812 : if (az >= 1.0) precnew += 1 + nbits2extraprec((long)(az/M_LN2));
224 314812 : if (issmall(n,&ki)) {
225 313881 : k = labs(ki);
226 313881 : n = utoi(k);
227 : } else {
228 847 : i = precision(n);
229 847 : if (i && i < precnew) n = gtofp(n,precnew);
230 : }
231 314728 : z = gtofp(z,precnew);
232 314728 : lim = bessel_get_lim(prec2nbits_mul(prec,M_LN2/2) / L, L);
233 314728 : z = gmul2n(gsqr(z),-2); if (J) z = gneg(z);
234 314728 : p1 = gprec_wtrunc(_jbessel(n,z,lim), prec);
235 314728 : return gc_upto(av, gmul(p2,p1));
236 : }
237 :
238 14 : case t_PADIC: pari_err_IMPL(stack_strcat("p-adic ",f));
239 392 : default:
240 : {
241 : long v, k, m;
242 392 : if (!(y = toser_i(z))) break;
243 238 : if (issmall(n,&ki)) n = utoi(labs(ki));
244 210 : y = gmul2n(gsqr(y),-2); if (J) y = gneg(y);
245 210 : v = valser(y);
246 210 : if (v < 0) pari_err_DOMAIN(f, "valuation", "<", gen_0, z);
247 203 : if (v == 0) pari_err_IMPL(stack_strcat(f, " around a!=0"));
248 203 : m = lg(y) - 2;
249 203 : k = m - (v >> 1);
250 203 : if (k <= 0) { set_avma(av); return scalarser(gen_1, varn(z), v); }
251 203 : setlg(y, k+2); return gc_upto(av, _jbessel(n, y, m));
252 : }
253 : }
254 154 : return trans_evalgen(f, (void*)n, J? vjbessel: vibessel, z, prec);
255 : }
256 : GEN
257 384972 : jbessel(GEN n, GEN z, long prec) { return jbesselintern(n,z,1,prec); }
258 : GEN
259 896 : ibessel(GEN n, GEN z, long prec) { return jbesselintern(n,z,0,prec); }
260 :
261 : /* k > 0 */
262 : static GEN
263 119 : _jbesselh(long k, GEN z, long prec)
264 : {
265 119 : GEN s, c, p0, p1, zinv = ginv(z);
266 : long i;
267 :
268 119 : gsincos(z,&s,&c,prec);
269 119 : p1 = gmul(zinv,s);
270 119 : p0 = p1; p1 = gmul(zinv,gsub(p0,c));
271 1134 : for (i = 2; i <= k; i++)
272 : {
273 1015 : GEN p2 = gsub(gmul(gmulsg(2*i-1,zinv), p1), p0);
274 1015 : p0 = p1; p1 = p2;
275 : }
276 119 : return p1;
277 : }
278 :
279 : /* J_{n+1/2}(z) */
280 : GEN
281 315 : jbesselh(GEN n, GEN z, long prec)
282 : {
283 : long k, i;
284 : pari_sp av;
285 : GEN y;
286 :
287 315 : if (typ(n)!=t_INT) pari_err_TYPE("jbesselh",n);
288 203 : k = itos(n);
289 203 : if (k < 0) return jbessel(gadd(ghalf,n), z, prec);
290 :
291 203 : switch(typ(z))
292 : {
293 133 : case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
294 : {
295 : long pr;
296 : GEN p1;
297 133 : if (gequal0(z))
298 : {
299 7 : av = avma;
300 7 : p1 = gmul(gsqrt(gdiv(z,mppi(prec)),prec),gpowgs(z,k));
301 7 : p1 = gdiv(p1, mulu_interval(k+1, 2*k+1)); /* x k! / (2k+1)! */
302 7 : return gc_upto(av, gmul2n(p1,2*k));
303 : }
304 126 : if ( (pr = precision(z)) ) prec = pr;
305 126 : if (bessel_asymp(n, z, prec2nbits(prec)))
306 7 : return jbessel(gadd(ghalf,n), z, prec);
307 119 : y = cgetc(prec); av = avma;
308 119 : p1 = gsqrt(gdiv(z, Pi2n(-1,prec)), prec);
309 119 : if (!k)
310 21 : p1 = gmul(p1, gsinc(z, prec));
311 : else
312 : {
313 98 : long bits = BITS_IN_LONG + 2*k * (log2(k) - dbllog2(z));
314 98 : if (bits > 0)
315 : {
316 98 : prec += nbits2extraprec(bits);
317 98 : if (pr) z = gtofp(z, prec);
318 : }
319 98 : p1 = gmul(p1, _jbesselh(k,z,prec));
320 : }
321 119 : set_avma(av); return affc_fixlg(p1, y);
322 : }
323 0 : case t_PADIC: pari_err_IMPL("p-adic jbesselh function");
324 70 : default:
325 : {
326 : long t, v;
327 70 : av = avma; if (!(y = toser_i(z))) break;
328 35 : if (gequal0(y)) return gc_upto(av, gpowgs(y,k));
329 35 : v = valser(y);
330 35 : if (v < 0) pari_err_DOMAIN("besseljh","valuation", "<", gen_0, z);
331 28 : t = lg(y)-2;
332 28 : if (v) y = sertoser(y, t + (2*k+1)*v);
333 28 : if (!k)
334 7 : y = gsinc(y,prec);
335 : else
336 : {
337 21 : GEN T, a = _jbesselh(k, y, prec);
338 21 : if (v) y = sertoser(y, t + k*v); /* lower precision */
339 21 : y = gdiv(a, gpowgs(y, k));
340 21 : T = cgetg(k+1, t_VECSMALL);
341 168 : for (i = 1; i <= k; i++) T[i] = 2*i+1;
342 21 : y = gmul(y, zv_prod_Z(T));
343 : }
344 28 : return gc_upto(av, y);
345 : }
346 : }
347 35 : return trans_evalgen("besseljh",(void*)n, vjbesselh, z, prec);
348 : }
349 :
350 : static GEN
351 0 : kbessel2(GEN nu, GEN x, long prec)
352 : {
353 0 : pari_sp av = avma;
354 0 : GEN p1, a, x2 = gshift(x,1);
355 :
356 0 : a = gtofp(gaddgs(gshift(nu,1), 1), prec);
357 0 : p1 = hyperu(gshift(a,-1), a, x2, prec);
358 0 : p1 = gmul(gmul(p1, gpow(x2,nu,prec)), sqrtr(mppi(prec)));
359 0 : return gc_upto(av, gmul(p1, gexp(gneg(x),prec)));
360 : }
361 :
362 : /* special case of hyperu */
363 : static GEN
364 14 : kbessel1(GEN nu, GEN gx, long prec)
365 : {
366 : GEN x, y, zf, r, u, pi, nu2;
367 14 : long bit, k, k2, n2, n, l = (typ(gx)==t_REAL)? realprec(gx): prec;
368 : pari_sp av;
369 :
370 14 : if (typ(nu)==t_COMPLEX) return kbessel2(nu, gx, l);
371 14 : y = cgetr(l); av = avma;
372 14 : x = gtofp(gx, l);
373 14 : nu = gtofp(nu,l); nu2 = sqrr(nu);
374 14 : shiftr_inplace(nu2,2); togglesign(nu2); /* nu2 = -4nu^2 */
375 14 : n = (long) (prec2nbits_mul(l,M_LN2) + M_PI*fabs(rtodbl(nu))) / 2;
376 14 : bit = prec2nbits(l) - 1;
377 14 : l += EXTRAPREC64;
378 14 : pi = mppi(l); n2 = n<<1; r = gmul2n(x,1);
379 14 : if (cmprs(x, n) < 0)
380 : {
381 14 : pari_sp av2 = avma;
382 14 : GEN q, v, c, s = real_1(l), t = real_0(l);
383 1246 : for (k = n2, k2 = 2*n2-1; k > 0; k--, k2 -= 2)
384 : {
385 1232 : GEN ak = divri(addri(nu2, sqru(k2)), mulss(n2<<2, -k));
386 1232 : s = addsr(1, mulrr(ak,s));
387 1232 : t = addsr(k2,mulrr(ak,t));
388 1232 : if (gc_needed(av2,3)) (void)gc_all(av2, 2, &s,&t);
389 : }
390 14 : shiftr_inplace(t, -1);
391 14 : q = utor(n2, l);
392 14 : zf = sqrtr(divru(pi,n2));
393 14 : u = gprec_wensure(mulrr(zf, s), l);
394 14 : v = gprec_wensure(divrs(addrr(mulrr(t,zf),mulrr(u,nu)),-n2), l);
395 : for(;;)
396 301 : {
397 315 : GEN p1, e, f, d = real_1(l);
398 : pari_sp av3;
399 315 : c = divur(5,q); if (expo(c) >= -1) c = real2n(-1,l);
400 315 : p1 = subsr(1, divrr(r,q)); if (cmprr(c,p1)>0) c = p1;
401 315 : togglesign(c); av3 = avma;
402 315 : e = u;
403 315 : f = v;
404 315 : for (k = 1;; k++)
405 35230 : {
406 35545 : GEN w = addrr(gmul2n(mulur(2*k-1,u), -1), mulrr(subrs(q,k),v));
407 35545 : w = addrr(w, mulrr(nu, subrr(u,gmul2n(v,1))));
408 35545 : u = divru(mulrr(q,v), k);
409 35545 : v = divru(w,k);
410 35545 : d = mulrr(d,c);
411 35545 : e = addrr(e, mulrr(d,u));
412 35545 : f = addrr(f, p1 = mulrr(d,v));
413 35545 : if (expo(p1) - expo(f) <= 1-prec2nbits(realprec(p1))) break;
414 35230 : if (gc_needed(av3,3)) (void)gc_all(av3,5,&u,&v,&d,&e,&f);
415 : }
416 315 : u = e;
417 315 : v = f;
418 315 : q = mulrr(q, addrs(c,1));
419 315 : if (expo(r) - expo(subrr(q,r)) >= bit) break;
420 301 : (void)gc_all(av2, 3, &u,&v,&q);
421 : }
422 14 : u = mulrr(u, gpow(divru(x,n),nu,prec));
423 : }
424 : else
425 : {
426 0 : GEN s, zz = ginv(gmul2n(r,2));
427 0 : pari_sp av2 = avma;
428 0 : s = real_1(l);
429 0 : for (k = n2, k2 = 2*n2-1; k > 0; k--, k2 -= 2)
430 : {
431 0 : GEN ak = divru(mulrr(addri(nu2, sqru(k2)), zz), k);
432 0 : s = subsr(1, mulrr(ak,s));
433 0 : if (gc_needed(av2,3)) s = gc_leaf(av2, s);
434 : }
435 0 : zf = sqrtr(divrr(pi,r));
436 0 : u = mulrr(s, zf);
437 : }
438 14 : affrr(mulrr(u, mpexp(negr(x))), y);
439 14 : set_avma(av); return y;
440 : }
441 :
442 : /* sum_{k=0}^m x^k (H(k)+H(k+n)) / (k! (k+n)!)
443 : * + sum_{k=0}^{n-1} (-x)^(k-n) (n-k-1)!/k! */
444 : static GEN
445 10860 : _kbessel(long n, GEN x, long m, long prec)
446 : {
447 : GEN p1, p2, s, H;
448 10860 : long k, M = m + n, exact = (M <= prec2nbits(prec));
449 : pari_sp av;
450 :
451 10860 : H = cgetg(M+2,t_VEC); gel(H,1) = gen_0;
452 10860 : if (exact)
453 : {
454 10853 : gel(H,2) = s = gen_1;
455 479307 : for (k=2; k<=M; k++) gel(H,k+1) = s = gdivgu(gaddsg(1,gmulsg(k,s)),k);
456 : }
457 : else
458 : {
459 7 : gel(H,2) = s = real_1(prec);
460 2877 : for (k=2; k<=M; k++) gel(H,k+1) = s = divru(addsr(1,mulur(k,s)),k);
461 : }
462 10860 : s = gadd(gel(H,m+1), gel(H,M+1)); av = avma;
463 430240 : for (k = m; k > 0; k--)
464 : {
465 419380 : s = gadd(gadd(gel(H,k),gel(H,k+n)), gdiv(gmul(x,s),mulss(k,k+n)));
466 419380 : if (gc_needed(av,1))
467 : {
468 0 : if (DEBUGMEM>1) pari_warn(warnmem,"_kbessel");
469 0 : s = gc_upto(av, s);
470 : }
471 : }
472 10860 : p1 = exact? mpfact(n): mpfactr(n,prec);
473 10860 : s = gdiv(s,p1);
474 10860 : if (n)
475 : {
476 8330 : x = gneg(ginv(x));
477 8330 : p2 = gmulsg(n, gdiv(x,p1));
478 8330 : s = gadd(s,p2);
479 62804 : for (k=n-1; k>0; k--)
480 : {
481 54474 : p2 = gmul(p2, gmul(mulss(k,n-k),x));
482 54474 : s = gadd(s,p2);
483 : }
484 : }
485 10860 : return s;
486 : }
487 :
488 : /* N = 1: Bessel N, else Bessel K */
489 : static GEN
490 12432 : kbesselintern(GEN n, GEN z, long N, long prec)
491 : {
492 12432 : const char *f = N? "besseln": "besselk";
493 : long i, k, ki, lim, precnew, fl2, ex, bit;
494 12432 : pari_sp av = avma;
495 : GEN p1, p2, y, p3, pp, pm, s, c;
496 : double az;
497 :
498 12432 : switch(typ(z))
499 : {
500 12047 : case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
501 12047 : if (gequal0(z)) pari_err_DOMAIN(f, "argument", "=", gen_0, z);
502 12047 : i = precision(z); if (i) prec = i;
503 12047 : i = precision(n); if (i && prec > i) prec = i;
504 12047 : bit = prec2nbits(prec);
505 12047 : if (bessel_asymp(n, z, bit))
506 : {
507 928 : GEN R = N? bessyasymp(n, z, bit): besskasymp(n, z, bit);
508 928 : if (typ(R) == t_COMPLEX && isexactzero(imag_i(n))
509 221 : && gsigne(real_i(z)) > 0
510 81 : && isexactzero(imag_i(z))) R = gcopy(gel(R,1));
511 928 : return gc_upto(av, R);
512 : }
513 : /* heuristic threshold */
514 11119 : if (!N && !gequal0(n) && gexpo(z) > bit/16 + gexpo(n))
515 14 : return kbessel1(n,z,prec);
516 11098 : az = dblmodulus(z); precnew = prec;
517 11098 : if (az >= 1) precnew += 1 + nbits2extraprec((long)((N?az:2*az)/M_LN2));
518 11098 : z = gtofp(z, precnew);
519 11098 : if (issmall(n,&ki))
520 : {
521 10776 : GEN z2 = gmul2n(z, -1), Z;
522 10776 : double B, L = HALF_E * az;
523 10776 : k = labs(ki);
524 10776 : B = prec2nbits_mul(prec,M_LN2/2) / L;
525 10776 : if (!N) B += 0.367879; /* exp(-1) */
526 10776 : lim = bessel_get_lim(B, L);
527 10776 : Z = gsqr(z2); if (N) Z = gneg(Z);
528 10776 : p1 = gmul(gpowgs(z2,k), _kbessel(k, Z, lim, precnew));
529 10776 : p2 = gadd(mpeuler(precnew), glog(z2,precnew));
530 10776 : p3 = jbesselintern(stoi(k),z,N,precnew);
531 10776 : p2 = gsub(gmul2n(p1,-1),gmul(p2,p3));
532 10776 : p2 = gprec_wtrunc(p2, prec);
533 10776 : if (N)
534 : {
535 8361 : p2 = gdiv(p2, Pi2n(-1,prec));
536 8361 : if (ki >= 0 || !odd(k)) p2 = gneg(p2);
537 : } else
538 2415 : if (odd(k)) p2 = gneg(p2);
539 10776 : return gc_GEN(av, p2);
540 : }
541 :
542 259 : n = gtofp(n, precnew);
543 259 : gsincos(gmul(n,mppi(precnew)), &s,&c,precnew);
544 259 : ex = gexpo(s);
545 259 : if (ex < 0) precnew += nbits2extraprec(N? -ex: -2*ex);
546 259 : if (i && i < precnew) {
547 84 : n = gtofp(n,precnew);
548 84 : z = gtofp(z,precnew);
549 84 : gsincos(gmul(n,mppi(precnew)), &s,&c,precnew);
550 : }
551 :
552 259 : pp = jbesselintern(n, z,N,precnew);
553 259 : pm = jbesselintern(gneg(n),z,N,precnew);
554 259 : if (N)
555 189 : p1 = gsub(gmul(c,pp),pm);
556 : else
557 70 : p1 = gmul(gsub(pm,pp), Pi2n(-1,precnew));
558 259 : p1 = gdiv(p1, s);
559 259 : return gc_GEN(av, gprec_wtrunc(p1,prec));
560 :
561 14 : case t_PADIC: pari_err_IMPL(stack_strcat("p-adic ",f));
562 371 : default:
563 371 : if (!(y = toser_i(z))) break;
564 217 : if (issmall(n,&ki))
565 : {
566 105 : long v, mv, k = labs(ki), m = lg(y)-2;
567 105 : y = gmul2n(gsqr(y),-2); if (N) y = gneg(y);
568 105 : v = valser(y);
569 105 : if (v < 0) pari_err_DOMAIN(f, "valuation", "<", gen_0, z);
570 91 : if (v == 0) pari_err_IMPL(stack_strcat(f, " around a!=0"));
571 91 : mv = m - (v >> 1);
572 91 : if (mv <= 0) { set_avma(av); return scalarser(gen_1, varn(z), v); }
573 84 : setlg(y, mv+2); return gc_GEN(av, _kbessel(k, y, m, prec));
574 : }
575 98 : if (!issmall(gmul2n(n,1),&ki))
576 70 : pari_err_DOMAIN(f, "2n mod Z", "!=", gen_0, n);
577 28 : k = labs(ki); n = gmul2n(stoi(k),-1);
578 28 : fl2 = (k&3)==1;
579 28 : pm = jbesselintern(gneg(n), y, N, prec);
580 28 : if (N) p1 = pm;
581 : else
582 : {
583 7 : pp = jbesselintern(n, y, N, prec);
584 7 : p2 = gpowgs(y,-k); if (fl2 == 0) p2 = gneg(p2);
585 7 : p3 = gmul2n(diviiexact(mpfact(k + 1),mpfact((k + 1) >> 1)),-(k + 1));
586 7 : p3 = gdivgu(gmul2n(gsqr(p3),1),k);
587 7 : p2 = gmul(p2,p3);
588 7 : p1 = gsub(pp,gmul(p2,pm));
589 : }
590 28 : return gc_upto(av, fl2? gneg(p1): gcopy(p1));
591 : }
592 154 : return trans_evalgen(f, (void*)n, N? vnbessel: vkbessel, z, prec);
593 : }
594 :
595 : GEN
596 3115 : kbessel(GEN n, GEN z, long prec) { return kbesselintern(n,z,0,prec); }
597 : GEN
598 9317 : ybessel(GEN n, GEN z, long prec) { return kbesselintern(n,z,1,prec); }
599 : /* J + iN */
600 : GEN
601 224 : hbessel1(GEN n, GEN z, long prec)
602 : {
603 224 : pari_sp av = avma;
604 224 : GEN J = jbessel(n,z,prec);
605 196 : GEN Y = ybessel(n,z,prec);
606 182 : return gc_upto(av, gadd(J, mulcxI(Y)));
607 : }
608 : /* J - iN */
609 : GEN
610 224 : hbessel2(GEN n, GEN z, long prec)
611 : {
612 224 : pari_sp av = avma;
613 224 : GEN J = jbessel(n,z,prec);
614 196 : GEN Y = ybessel(n,z,prec);
615 182 : return gc_upto(av, gadd(J, mulcxmI(Y)));
616 : }
617 :
618 : static GEN
619 1008 : besselrefine(GEN z, GEN nu, GEN (*B)(GEN,GEN,long), long bit)
620 : {
621 1008 : GEN z0 = gprec_w(z, DEFAULTPREC), nu1 = gaddgs(nu, 1), t;
622 1008 : long e, n, c, j, prec = DEFAULTPREC;
623 :
624 1008 : t = gdiv(B(nu1, z0, prec), B(nu, z0, prec));
625 1008 : t = gadd(z0, gdiv(gsub(gsqr(z0), gsqr(nu)), gsub(gdiv(nu, z0), t)));
626 1008 : e = gexpo(t) - 2 * gexpo(z0) - 1; if (e < 0) e = 0;
627 1008 : n = expu(bit + 32 - e);
628 1008 : c = 1 + e + ((bit - e) >> n);
629 8064 : for (j = 1; j <= n; j++)
630 : {
631 7056 : c = 2 * c - e;
632 7056 : prec = nbits2prec(c); z = gprec_w(z, prec);
633 7056 : t = gdiv(B(nu1, z, prec), B(nu, z, prec));
634 7056 : z = gsub(z, ginv(gsub(gdiv(nu, z), t)));
635 : }
636 1008 : return gprec_w(z, nbits2prec(bit));
637 : }
638 :
639 : /* solve tan(fi) - fi = y, y >= 0; Temme's method */
640 : static double
641 700 : fi(double y)
642 : {
643 : double p, pp, r;
644 700 : if (y == 0) return 0;
645 700 : if (y > 100000) return M_PI/2;
646 700 : if (y < 1)
647 : {
648 455 : p = pow(3*y, 1.0/3); pp = p * p;
649 455 : p = p * (1 + pp * (-210 * pp + (27 - 2*pp)) / 1575);
650 : }
651 : else
652 : {
653 245 : p = 1 / (y + M_PI/2); pp = p * p;
654 245 : p = M_PI/2 - p*(1 + pp*(2310 + pp*(3003 + pp*(4818 + pp*(8591 + pp*16328)))) / 3465);
655 : }
656 700 : pp = (y + p) * (y + p); r = (p - atan(p + y)) / pp;
657 700 : return p - (1 + pp) * r * (1 + r / (p + y));
658 : }
659 :
660 : static GEN
661 1022 : besselzero(GEN nu, long n, GEN (*B)(GEN,GEN,long), long bit)
662 : {
663 1022 : pari_sp av = avma;
664 1022 : long prec = nbits2prec(bit);
665 1022 : int J = B == jbessel;
666 : GEN z;
667 1022 : if (n <= 0) pari_err_DOMAIN("besselzero", "n", "<=", gen_0, stoi(n));
668 1008 : if (n > LONG_MAX / 4) pari_err_OVERFLOW("besselzero");
669 1008 : if (is_real_t(typ(nu)) && gsigne(nu) >= 0)
670 1008 : { /* Temme */
671 1008 : double x, c, b, a = gtodouble(nu), t = J? 0.25: 0.75;
672 1008 : if (n >= 3*a - 8)
673 : {
674 308 : double aa = a*a, mu = 4*aa, mu2 = mu*mu, p, p0, p1, q1;
675 308 : p = 7 * mu - 31; p0 = mu-1;
676 308 : if (1 + p == p) /* p large */
677 0 : p1 = q1 = 0;
678 : else
679 : {
680 308 : p1 = 4 * (253 * mu2 - 3722 * mu + 17869) / (15 * p);
681 308 : q1 = 1.6 * (83 * mu2 - 982 * mu + 3779) / p;
682 : }
683 308 : b = (n + a/2 - t) * M_PI;
684 308 : c = 1 / (64 * b * b);
685 308 : x = b - p0 * (1 - p1 * c) / (8 * b * (1 - q1 * c));
686 : }
687 : else
688 : {
689 700 : double u, v, w, xx, bb = a >= 3? pow(a, -2./3): 1;
690 700 : if (n == 1)
691 336 : x = J? -2.33811: -1.17371;
692 : else
693 : {
694 364 : double pp1 = 5./48, qq1 = -5./36, y = 3./8 * M_PI;
695 364 : x = 4 * y * (n - t); v = 1 / (x*x);
696 364 : x = - pow(x, 2.0/3) * (1 + v * (pp1 + qq1 * v));
697 : }
698 700 : u = x * bb; v = fi(2.0/3 * pow(-u, 1.5));
699 700 : w = 1 / cos(v); xx = 1 - w*w; c = sqrt(u/xx);
700 700 : x = w * (a + c / (48*a*u) * (-5/u-c * (-10/xx + 6)));
701 : }
702 1008 : z = dbltor(x);
703 : }
704 : else
705 : { /* generic, hope for the best */
706 0 : long a = 4 * n - (J? 1: 3);
707 : GEN b, m;
708 0 : b = gmul(mppi(prec), gmul2n(gaddgs(gmul2n(nu,1), a), -2));
709 0 : m = gmul2n(gsqr(nu),2);
710 0 : z = gsub(b, gdiv(gsubgs(m, 1), gmul2n(b, 3)));
711 : }
712 1008 : return gc_GEN(av, besselrefine(z, nu, B, bit));
713 : }
714 : GEN
715 511 : besseljzero(GEN nu, long k, long b) { return besselzero(nu, k, jbessel, b); }
716 : GEN
717 511 : besselyzero(GEN nu, long k, long b) { return besselzero(nu, k, ybessel, b); }
718 :
719 : /***********************************************************************/
720 : /** INCOMPLETE GAMMA FUNCTION **/
721 : /***********************************************************************/
722 : /* mx ~ |x|, b = bit accuracy */
723 : static int
724 16765 : gamma_use_asymp(GEN x, long b)
725 : {
726 : long e;
727 16765 : if (is_real_t(typ(x)))
728 : {
729 12733 : pari_sp av = avma;
730 12733 : return gc_int(av, gcmpgs(R_abs_shallow(x), 3*b / 4) >= 0);
731 : }
732 4032 : e = gexpo(x); return e >= b || dblmodulus(x) >= 3*b / 4;
733 : }
734 : /* x a t_REAL */
735 : static GEN
736 28 : eint1r_asymp(GEN x, GEN expx, long prec)
737 : {
738 28 : pari_sp av = avma, av2;
739 : GEN S, q, z, ix;
740 28 : long oldeq = LONG_MAX, esx = -prec2nbits(prec), j;
741 :
742 28 : if (realprec(x) < prec + EXTRAPREC64) x = rtor(x, prec+EXTRAPREC64);
743 28 : ix = invr(x); q = z = negr(ix);
744 28 : av2 = avma; S = addrs(q, 1);
745 28 : for (j = 2;; j++)
746 1211 : {
747 1239 : long eq = expo(q); if (eq < esx) break;
748 1211 : if ((j & 3) == 0)
749 : { /* guard against divergence */
750 294 : if (eq > oldeq) return gc_NULL(av); /* regressing, abort */
751 294 : oldeq = eq;
752 : }
753 1211 : q = mulrr(q, mulru(z, j)); S = addrr(S, q);
754 1211 : if (gc_needed(av2, 1)) (void)gc_all(av2, 2, &S, &q);
755 : }
756 28 : if (DEBUGLEVEL > 2) err_printf("eint1: using asymp\n");
757 28 : S = expx? divrr(S, expx): mulrr(S, mpexp(negr(x)));
758 28 : return gc_leaf(av, mulrr(S, ix));
759 : }
760 : /* cf incgam_asymp(0, x); z = -1/x
761 : * exp(-x)/x * (1 + z + 2! z^2 + ...) */
762 : static GEN
763 105 : eint1_asymp(GEN x, GEN expx, long prec)
764 : {
765 105 : pari_sp av = avma, av2;
766 : GEN S, q, z, ix;
767 105 : long oldeq = LONG_MAX, esx = -prec2nbits(prec), j;
768 :
769 105 : if (typ(x) != t_REAL) x = gtofp(x, prec+EXTRAPREC64);
770 105 : if (typ(x) == t_REAL) return eint1r_asymp(x, expx, prec);
771 105 : ix = ginv(x); q = z = gneg_i(ix);
772 105 : av2 = avma; S = gaddgs(q, 1);
773 105 : for (j = 2;; j++)
774 5824 : {
775 5929 : long eq = gexpo(q); if (eq < esx) break;
776 5824 : if ((j & 3) == 0)
777 : { /* guard against divergence */
778 1442 : if (eq > oldeq) return gc_NULL(av); /* regressing, abort */
779 1442 : oldeq = eq;
780 : }
781 5824 : q = gmul(q, gmulgu(z, j)); S = gadd(S, q);
782 5824 : if (gc_needed(av2, 1)) (void)gc_all(av2, 2, &S, &q);
783 : }
784 105 : if (DEBUGLEVEL > 2) err_printf("eint1: using asymp\n");
785 105 : S = expx? gdiv(S, expx): gmul(S, gexp(gneg_i(x), prec));
786 105 : return gc_upto(av, gmul(S, ix));
787 : }
788 :
789 : /* eint1(x) = incgam(0, x); typ(x) = t_REAL, x > 0 */
790 : static GEN
791 6524 : eint1p(GEN x, GEN expx)
792 : {
793 : pari_sp av;
794 6524 : long prec = realprec(x), bit = prec2nbits(prec), i;
795 : double mx;
796 : GEN z, S, t, H, run;
797 :
798 6524 : if (gamma_use_asymp(x, bit)
799 28 : && (z = eint1r_asymp(x, expx, prec))) return z;
800 6496 : mx = rtodbl(x);
801 6496 : if (mx > 1)
802 3591 : prec += nbits2extraprec((mx+log(mx))/M_LN2 + 10);
803 : else
804 2905 : prec += EXTRAPREC64;
805 6496 : bit = prec2nbits(prec);
806 6496 : run = real_1(prec); x = rtor(x, prec);
807 6496 : av = avma; S = z = t = H = run;
808 618178 : for (i = 2; expo(S) - expo(t) <= bit; i++)
809 : {
810 611682 : H = addrr(H, divru(run,i)); /* H = sum_{k<=i} 1/k */
811 611682 : z = divru(mulrr(x,z), i); /* z = x^(i-1)/i! */
812 611682 : t = mulrr(z, H); S = addrr(S, t);
813 611682 : if ((i & 0x1ff) == 0) (void)gc_all(av, 4, &z,&t,&S,&H);
814 : }
815 6496 : return subrr(mulrr(x, divrr(S,expx? expx: mpexp(x))),
816 : addrr(mplog(x), mpeuler(prec)));
817 : }
818 : /* eint1(x) = incgam(0, x); typ(x) = t_REAL, x < 0
819 : * rewritten from code contributed by Manfred Radimersky */
820 : static GEN
821 140 : eint1m(GEN x, GEN expx)
822 : {
823 140 : GEN p1, q, S, y, z = cgetg(3, t_COMPLEX);
824 140 : long l = realprec(x), n = prec2nbits(l), j;
825 140 : pari_sp av = avma;
826 :
827 140 : y = rtor(x, l + EXTRAPREC64); setsigne(y,1); /* |x| */
828 140 : if (gamma_use_asymp(y, n))
829 : { /* ~eint1_asymp: asymptotic expansion */
830 14 : p1 = q = invr(y); S = addrs(q, 1);
831 560 : for (j = 2; expo(q) >= -n; j++) {
832 546 : q = mulrr(q, mulru(p1, j));
833 546 : S = addrr(S, q);
834 : }
835 14 : y = mulrr(p1, expx? divrr(S, expx): mulrr(S, mpexp(y)));
836 : }
837 : else
838 : {
839 126 : p1 = q = S = y;
840 24248 : for (j = 2; expo(q) - expo(S) >= -n; j++) {
841 24122 : p1 = mulrr(y, divru(p1, j)); /* (-x)^j/j! */
842 24122 : q = divru(p1, j);
843 24122 : S = addrr(S, q);
844 : }
845 126 : y = addrr(S, addrr(logr_abs(x), mpeuler(l)));
846 : }
847 140 : y = gc_leaf(av, y); togglesign(y);
848 140 : gel(z, 1) = y;
849 140 : y = mppi(l); setsigne(y, -1);
850 140 : gel(z, 2) = y; return z;
851 : }
852 :
853 : /* real(z*log(z)-z), z = x+iy */
854 : static double
855 8372 : mygamma(double x, double y)
856 : {
857 8372 : if (x == 0.) return -(M_PI/2)*fabs(y);
858 8372 : return (x/2)*log(x*x+y*y)-x-y*atan(y/x);
859 : }
860 :
861 : /* x^s exp(-x) */
862 : static GEN
863 10843 : expmx_xs(GEN s, GEN x, GEN logx, long prec)
864 : {
865 : GEN z;
866 10843 : long ts = typ(s);
867 10843 : if (ts == t_INT || (ts == t_FRAC && absequaliu(gel(s,2), 2)))
868 5264 : z = gmul(gexp(gneg(x), prec), gpow(x, s, prec));
869 : else
870 5579 : z = gexp(gsub(gmul(s, logx? logx: glog(x,prec+EXTRAPREC64)), x), prec);
871 10843 : return z;
872 : }
873 :
874 : /* Not yet: doesn't work at low accuracy
875 : #define INCGAM_CF
876 : */
877 :
878 : #ifdef INCGAM_CF
879 : /* Is s very close to a nonpositive integer ? */
880 : static int
881 : isgammapole(GEN s, long bitprec)
882 : {
883 : pari_sp av = avma;
884 : GEN t = imag_i(s);
885 : long e, b = bitprec - 10;
886 :
887 : if (gexpo(t) > - b) return 0;
888 : s = real_i(s);
889 : if (gsigne(s) > 0 && gexpo(s) > -b) return 0;
890 : (void)grndtoi(s, &e); return gc_bool(av, e < -b);
891 : }
892 :
893 : /* incgam using the continued fraction. x a t_REAL or t_COMPLEX, mx ~ |x|.
894 : * Assume precision(s), precision(x) >= prec */
895 : static GEN
896 : incgam_cf(GEN s, GEN x, double mx, long prec)
897 : {
898 : GEN ms, y, S;
899 : long n, i, j, LS, bitprec = prec2nbits(prec);
900 : double rs, is, m;
901 :
902 : if (typ(s) == t_COMPLEX)
903 : {
904 : rs = gtodouble(gel(s,1));
905 : is = gtodouble(gel(s,2));
906 : }
907 : else
908 : {
909 : rs = gtodouble(s);
910 : is = 0.;
911 : }
912 : if (isgammapole(s, bitprec)) LS = 0;
913 : else
914 : {
915 : double bit, LGS = mygamma(rs,is);
916 : LS = LGS <= 0 ? 0: ceil(LGS);
917 : bit = (LGS - (rs-1)*log(mx) + mx)/M_LN2;
918 : if (bit > 0)
919 : {
920 : prec += nbits2extraprec((long)bit);
921 : x = gtofp(x, prec);
922 : if (isinexactreal(s)) s = gtofp(s, prec);
923 : }
924 : }
925 : /* |ln(2*gamma(s)*sin(s*Pi))| <= ln(2) + |lngamma(s)| + |Im(s)*Pi|*/
926 : m = bitprec*M_LN2 + LS + M_LN2 + fabs(is)*M_PI + mx;
927 : if (rs < 1) m += (1 - rs)*log(mx);
928 : m /= 4;
929 : n = (long)(1 + m*m/mx);
930 : y = expmx_xs(gsubgs(s,1), x, NULL, prec);
931 : if (rs >= 0 && bitprec >= 512)
932 : {
933 : GEN A = cgetg(n+1, t_VEC), B = cgetg(n+1, t_VEC);
934 : ms = gsubsg(1, s);
935 : for (j = 1; j <= n; ++j)
936 : {
937 : gel(A,j) = ms;
938 : gel(B,j) = gmulsg(j, gsubgs(s,j));
939 : ms = gaddgs(ms, 2);
940 : }
941 : S = contfraceval_inv(mkvec2(A,B), x, -1);
942 : }
943 : else
944 : {
945 : GEN x_s = gsub(x, s);
946 : pari_sp av2 = avma;
947 : S = gdiv(gsubgs(s,n), gaddgs(x_s,n<<1));
948 : for (i=n-1; i >= 1; i--)
949 : {
950 : S = gdiv(gsubgs(s,i), gadd(gaddgs(x_s,i<<1),gmulsg(i,S)));
951 : if (gc_needed(av2,3))
952 : {
953 : if(DEBUGMEM>1) pari_warn(warnmem,"incgam_cf");
954 : S = gc_upto(av2, S);
955 : }
956 : }
957 : S = gaddgs(S,1);
958 : }
959 : return gmul(y, S);
960 : }
961 : #endif
962 :
963 : static double
964 6419 : findextraincgam(GEN s, GEN x)
965 : {
966 6419 : double sig = gtodouble(real_i(s)), t = gtodouble(imag_i(s));
967 6419 : double xr = gtodouble(real_i(x)), xi = gtodouble(imag_i(x));
968 6419 : double exd = 0., Nx = xr*xr + xi*xi, D = Nx - t*t;
969 : long n;
970 :
971 6419 : if (xr < 0)
972 : {
973 833 : long ex = gexpo(x);
974 833 : if (ex > 0 && ex > gexpo(s)) exd = sqrt(Nx)*log(Nx)/2; /* |x| log |x| */
975 : }
976 6419 : if (D <= 0.) return exd;
977 4977 : n = (long)(sqrt(D)-sig);
978 4977 : if (n <= 0) return exd;
979 1841 : return maxdd(exd, (n*log(Nx)/2 - mygamma(sig+n, t) + mygamma(sig, t)) / M_LN2);
980 : }
981 :
982 : /* use exp(-x) * (x^s/s) * sum_{k >= 0} x^k / prod(i=1, k, s+i) */
983 : static GEN
984 6426 : incgamc_i(GEN s, GEN x, long *ptexd, long prec)
985 : {
986 : GEN S, t, y;
987 : long l, n, i, exd;
988 6426 : pari_sp av = avma, av2;
989 :
990 6426 : if (gequal0(x))
991 : {
992 7 : if (ptexd) *ptexd = 0.;
993 7 : return gtofp(x, prec);
994 : }
995 6419 : l = precision(x);
996 6419 : if (!l) l = prec;
997 6419 : n = -prec2nbits(l)-1;
998 6419 : exd = (long)findextraincgam(s, x);
999 6419 : if (ptexd) *ptexd = exd;
1000 6419 : if (exd > 0)
1001 : {
1002 1666 : long p = l + nbits2extraprec(exd);
1003 1666 : x = gtofp(x, p);
1004 1666 : if (isinexactreal(s)) s = gtofp(s, p);
1005 : }
1006 4753 : else x = gtofp(x, l+EXTRAPREC64);
1007 6419 : av2 = avma;
1008 6419 : S = gdiv(x, gaddsg(1,s));
1009 6419 : t = gaddsg(1, S);
1010 770875 : for (i=2; gexpo(S) >= n; i++)
1011 : {
1012 764456 : S = gdiv(gmul(x,S), gaddsg(i,s)); /* x^i / ((s+1)...(s+i)) */
1013 764456 : t = gadd(S,t);
1014 764456 : if (gc_needed(av2,3))
1015 : {
1016 0 : if(DEBUGMEM>1) pari_warn(warnmem,"incgamc");
1017 0 : (void)gc_all(av2, 2, &S, &t);
1018 : }
1019 : }
1020 6419 : y = expmx_xs(s, x, NULL, prec);
1021 6419 : return gc_upto(av, gmul(gdiv(y,s), t));
1022 : }
1023 :
1024 : GEN
1025 2226 : incgamc(GEN s, GEN x, long prec)
1026 2226 : { return incgamc_i(s, x, NULL, prec); }
1027 :
1028 : /* incgamma using asymptotic expansion:
1029 : * exp(-x)x^(s-1)(1 + (s-1)/x + (s-1)(s-2)/x^2 + ...) */
1030 : static GEN
1031 2716 : incgam_asymp(GEN s, GEN x, long prec)
1032 : {
1033 2716 : pari_sp av = avma, av2;
1034 : GEN S, q, cox, invx;
1035 2716 : long oldeq = LONG_MAX, eq, esx, j;
1036 2716 : int flint = (typ(s) == t_INT && signe(s) > 0);
1037 :
1038 2716 : x = gtofp(x,prec+EXTRAPREC64);
1039 2716 : invx = ginv(x);
1040 2716 : esx = -prec2nbits(prec);
1041 2716 : av2 = avma;
1042 2716 : q = gmul(gsubgs(s,1), invx);
1043 2716 : S = gaddgs(q, 1);
1044 2716 : for (j = 2;; j++)
1045 : {
1046 123879 : eq = gexpo(q); if (eq < esx) break;
1047 121282 : if (!flint && (j & 3) == 0)
1048 : { /* guard against divergence */
1049 15778 : if (eq > oldeq) return gc_NULL(av); /* regressing, abort */
1050 15659 : oldeq = eq;
1051 : }
1052 121163 : q = gmul(q, gmul(gsubgs(s,j), invx));
1053 121163 : S = gadd(S, q);
1054 121163 : if (gc_needed(av2, 1)) (void)gc_all(av2, 2, &S, &q);
1055 : }
1056 2597 : if (DEBUGLEVEL > 2) err_printf("incgam: using asymp\n");
1057 2597 : cox = expmx_xs(gsubgs(s,1), x, NULL, prec);
1058 2597 : return gc_upto(av, gmul(cox, S));
1059 : }
1060 :
1061 : /* gasx = incgam(s-n,x). Compute incgam(s,x)
1062 : * = (s-1)(s-2)...(s-n)gasx + exp(-x)x^(s-1) *
1063 : * (1 + (s-1)/x + ... + (s-1)(s-2)...(s-n+1)/x^(n-1)) */
1064 : static GEN
1065 546 : incgam_asymp_partial(GEN s, GEN x, GEN gasx, long n, long prec)
1066 : {
1067 : pari_sp av;
1068 546 : GEN S, q, cox, invx, s1 = gsubgs(s, 1), sprod;
1069 : long j;
1070 546 : cox = expmx_xs(s1, x, NULL, prec);
1071 546 : if (n == 1) return gadd(cox, gmul(s1, gasx));
1072 546 : invx = ginv(x);
1073 546 : av = avma;
1074 546 : q = gmul(s1, invx);
1075 546 : S = gaddgs(q, 1);
1076 52164 : for (j = 2; j < n; j++)
1077 : {
1078 51618 : q = gmul(q, gmul(gsubgs(s, j), invx));
1079 51618 : S = gadd(S, q);
1080 51618 : if (gc_needed(av, 2)) (void)gc_all(av, 2, &S, &q);
1081 : }
1082 546 : sprod = gmul(gmul(q, gpowgs(x, n-1)), gsubgs(s, n));
1083 546 : return gadd(gmul(cox, S), gmul(sprod, gasx));
1084 : }
1085 :
1086 : /* Assume s != 0; called when Re(s) <= 1/2 */
1087 : static GEN
1088 2401 : incgamspec(GEN s, GEN x, GEN g, long prec)
1089 : {
1090 2401 : GEN q, S, cox = gen_0, P, sk, S1, S2, S3, F3, logx, mx;
1091 2401 : long n, esk, E, k = itos(ground(gneg(real_i(s)))); /* >= 0 */
1092 :
1093 2401 : if (k && gexpo(x) > 0)
1094 : {
1095 245 : GEN xk = gdivgu(x, k);
1096 245 : long bitprec = prec2nbits(prec);
1097 245 : double d = (gexpo(xk) > bitprec)? bitprec*M_LN2: log(dblmodulus(xk));
1098 245 : d = k * (d + 1) / M_LN2;
1099 245 : if (d > 0) prec += nbits2extraprec((long)d);
1100 245 : if (isinexactreal(s)) s = gtofp(s, prec);
1101 : }
1102 2401 : x = gtofp(x, maxss(precision(x), prec) + EXTRAPREC64);
1103 2401 : sk = gaddgs(s, k); /* |Re(sk)| <= 1/2 */
1104 2401 : logx = glog(x, prec);
1105 2401 : mx = gneg(x);
1106 2401 : if (k == 0) { S = gen_0; P = gen_1; }
1107 : else
1108 : {
1109 : long j;
1110 854 : q = ginv(s); S = q; P = s;
1111 16926 : for (j = 1; j < k; j++)
1112 : {
1113 16072 : GEN sj = gaddgs(s, j);
1114 16072 : q = gmul(q, gdiv(x, sj));
1115 16072 : S = gadd(S, q);
1116 16072 : P = gmul(P, sj);
1117 : }
1118 854 : cox = expmx_xs(s, x, logx, prec); /* x^s exp(-x) */
1119 854 : S = gmul(S, gneg(cox));
1120 : }
1121 2401 : if (k && gequal0(sk))
1122 175 : return gadd(S, gdiv(eint1(x, prec), P));
1123 2226 : esk = gexpo(sk);
1124 2226 : if (esk > -7)
1125 : {
1126 1015 : GEN a, b, PG = gmul(sk, P);
1127 1015 : if (g) g = gmul(g, PG);
1128 1015 : a = incgam0(gaddgs(sk,1), x, g, prec);
1129 1015 : if (k == 0) cox = expmx_xs(s, x, logx, prec);
1130 1015 : b = gmul(gpowgs(x, k), cox);
1131 1015 : return gadd(S, gdiv(gsub(a, b), PG));
1132 : }
1133 1211 : E = prec2nbits(prec) + 1;
1134 1211 : if (gexpo(x) > 0)
1135 : {
1136 420 : long X = (long)(dblmodulus(x)/M_LN2);
1137 420 : prec += 2*nbits2extraprec(X);
1138 420 : x = gtofp(x, prec); mx = gneg(x);
1139 420 : logx = glog(x, prec); sk = gtofp(sk, prec);
1140 420 : E += X;
1141 : }
1142 1211 : if (isinexactreal(sk)) sk = gtofp(sk, prec+EXTRAPREC64);
1143 : /* |sk| < 2^-7 is small, guard against cancellation */
1144 1211 : F3 = gexpm1(gmul(sk, logx), prec);
1145 : /* ( gamma(1+sk) - exp(sk log(x))) ) / sk */
1146 1211 : S1 = gdiv(gsub(ggamma1m1(sk, prec+EXTRAPREC64), F3), sk);
1147 1211 : q = x; S3 = gdiv(x, gaddsg(1,sk));
1148 255523 : for (n = 2; gexpo(q) - gexpo(S3) > -E; n++)
1149 : {
1150 254312 : q = gmul(q, gdivgu(mx, n));
1151 254312 : S3 = gadd(S3, gdiv(q, gaddsg(n, sk)));
1152 : }
1153 1211 : S2 = gadd(gadd(S1, S3), gmul(F3, S3));
1154 1211 : return gadd(S, gdiv(S2, P));
1155 : }
1156 :
1157 : /* return |x| */
1158 : double
1159 14431332 : dblmodulus(GEN x)
1160 : {
1161 14431332 : if (typ(x) == t_COMPLEX)
1162 : {
1163 1807267 : double a = gtodouble(gel(x,1));
1164 1807267 : double b = gtodouble(gel(x,2));
1165 1807267 : return sqrt(a*a + b*b);
1166 : }
1167 : else
1168 12624065 : return fabs(gtodouble(x));
1169 : }
1170 :
1171 : /* Driver routine. If g != NULL, assume that g=gamma(s,prec). */
1172 : GEN
1173 11564 : incgam0(GEN s, GEN x, GEN g, long prec)
1174 : {
1175 : pari_sp av;
1176 : long E, l;
1177 : GEN z, rs, is;
1178 :
1179 11564 : if (gequal0(x)) return g? gcopy(g): ggamma(s,prec);
1180 11564 : if (gequal0(s)) return eint1(x, prec);
1181 9744 : l = precision(s); if (!l) l = prec;
1182 9744 : E = prec2nbits(l);
1183 9744 : if (gamma_use_asymp(x, E) ||
1184 8323 : (typ(s) == t_INT && signe(s) > 0 && gexpo(x) >= expi(s)))
1185 2716 : if ((z = incgam_asymp(s, x, l))) return z;
1186 7147 : av = avma; E++;
1187 7147 : rs = real_i(s);
1188 7147 : is = imag_i(s);
1189 : #ifdef INCGAM_CF
1190 : /* Can one use continued fraction ? */
1191 : if (gequal0(is) && gequal0(imag_i(x)) && gsigne(x) > 0)
1192 : {
1193 : double sd = gtodouble(rs), LB, UB;
1194 : double xd = gtodouble(real_i(x));
1195 : if (sd > 0) {
1196 : LB = 15 + 0.1205*E;
1197 : UB = 5 + 0.752*E;
1198 : } else {
1199 : LB = -6 + 0.1205*E;
1200 : UB = 5 + 0.752*E + fabs(sd)/54.;
1201 : }
1202 : if (xd >= LB && xd <= UB)
1203 : {
1204 : if (DEBUGLEVEL > 2) err_printf("incgam: using continued fraction\n");
1205 : return gc_upto(av, incgam_cf(s, x, xd, prec));
1206 : }
1207 : }
1208 : #endif
1209 7147 : if (gsigne(rs) > 0 && gexpo(rs) >= -1)
1210 : { /* use complementary incomplete gamma */
1211 4746 : long n, egs, exd, precg, es = gexpo(s);
1212 4746 : if (es < 0) {
1213 602 : l += nbits2extraprec(-es) + 1;
1214 602 : x = gtofp(x, l);
1215 602 : if (isinexactreal(s)) s = gtofp(s, l);
1216 : }
1217 4746 : n = itos(gceil(rs));
1218 4746 : if (n > 100)
1219 : {
1220 : GEN gasx;
1221 546 : n -= 100;
1222 546 : if (es > 0)
1223 : {
1224 546 : es = mygamma(gtodouble(rs) - n, gtodouble(is)) / M_LN2;
1225 546 : if (es > 0)
1226 : {
1227 546 : l += nbits2extraprec(es);
1228 546 : x = gtofp(x, l);
1229 546 : if (isinexactreal(s)) s = gtofp(s, l);
1230 : }
1231 : }
1232 546 : gasx = incgam0(gsubgs(s, n), x, NULL, prec);
1233 546 : return gc_upto(av, incgam_asymp_partial(s, x, gasx, n, prec));
1234 : }
1235 4200 : if (DEBUGLEVEL > 2) err_printf("incgam: using power series 1\n");
1236 : /* egs ~ expo(gamma(s)) */
1237 4200 : precg = g? precision(g): 0;
1238 4200 : egs = g? gexpo(g): (long)(mygamma(gtodouble(rs), gtodouble(is)) / M_LN2);
1239 4200 : if (egs > 0) {
1240 1946 : l += nbits2extraprec(egs) + 1;
1241 1946 : x = gtofp(x, l);
1242 1946 : if (isinexactreal(s)) s = gtofp(s, l);
1243 1946 : if (precg < l) g = NULL;
1244 : }
1245 4200 : z = incgamc_i(s, x, &exd, l);
1246 4200 : if (exd > 0)
1247 : {
1248 896 : l += nbits2extraprec(exd);
1249 896 : if (isinexactreal(s)) s = gtofp(s, l);
1250 896 : if (precg < l) g = NULL;
1251 : }
1252 : else
1253 : { /* gamma(s) negligible ? Compute to lower accuracy */
1254 3304 : long e = gexpo(z) - egs;
1255 3304 : if (e > 3)
1256 : {
1257 420 : E -= e;
1258 420 : if (E <= 0) g = gen_0; else if (!g) g = ggamma(s, nbits2prec(E));
1259 : }
1260 : }
1261 : /* worry about possible cancellation */
1262 4200 : if (!g) g = ggamma(s, maxss(l,precision(z)));
1263 4200 : return gc_upto(av, gsub(g,z));
1264 : }
1265 2401 : if (DEBUGLEVEL > 2) err_printf("incgam: using power series 2\n");
1266 2401 : return gc_upto(av, incgamspec(s, x, g, l));
1267 : }
1268 :
1269 : GEN
1270 1106 : incgam(GEN s, GEN x, long prec) { return incgam0(s, x, NULL, prec); }
1271 :
1272 : /* x a t_REAL */
1273 : GEN
1274 2940 : mpeint1(GEN x, GEN expx)
1275 : {
1276 2940 : long s = signe(x);
1277 : pari_sp av;
1278 : GEN z;
1279 2940 : if (!s) pari_err_DOMAIN("eint1", "x","=",gen_0, x);
1280 2933 : if (s < 0) return eint1m(x, expx);
1281 2793 : z = cgetr(realprec(x));
1282 2793 : av = avma; affrr(eint1p(x, expx), z);
1283 2793 : set_avma(av); return z;
1284 : }
1285 :
1286 : static GEN
1287 357 : cxeint1(GEN x, long prec)
1288 : {
1289 357 : pari_sp av = avma, av2;
1290 : GEN q, S, run, z, H;
1291 357 : long n, E = prec2nbits(prec);
1292 :
1293 357 : if (gamma_use_asymp(x, E) && (z = eint1_asymp(x, NULL, prec))) return z;
1294 252 : E++;
1295 252 : if (gexpo(x) > 0)
1296 : { /* take cancellation into account, log2(\sum |x|^n / n!) = |x| / log(2) */
1297 42 : double dbx = dblmodulus(x);
1298 42 : long X = (long)((dbx + log(dbx))/M_LN2 + 10);
1299 42 : prec += nbits2extraprec(X);
1300 42 : x = gtofp(x, prec); E += X;
1301 : }
1302 252 : if (DEBUGLEVEL > 2) err_printf("eint1: using power series\n");
1303 252 : run = real_1(prec);
1304 252 : av2 = avma;
1305 252 : S = z = q = H = run;
1306 48384 : for (n = 2; gexpo(q) - gexpo(S) >= -E; n++)
1307 : {
1308 48132 : H = addrr(H, divru(run, n)); /* H = sum_{k<=n} 1/k */
1309 48132 : z = gdivgu(gmul(x,z), n); /* z = x^(n-1)/n! */
1310 48132 : q = gmul(z, H); S = gadd(S, q);
1311 48132 : if ((n & 0x1ff) == 0) (void)gc_all(av2, 4, &z, &q, &S, &H);
1312 : }
1313 252 : S = gmul(gmul(x, S), gexp(gneg_i(x), prec));
1314 252 : return gc_upto(av, gsub(S, gadd(glog(x, prec), mpeuler(prec))));
1315 : }
1316 :
1317 : GEN
1318 3297 : eint1(GEN x, long prec)
1319 : {
1320 3297 : switch(typ(x))
1321 : {
1322 357 : case t_COMPLEX: return cxeint1(x, prec);
1323 2541 : case t_REAL: break;
1324 399 : default: x = gtofp(x, prec);
1325 : }
1326 2940 : return mpeint1(x,NULL);
1327 : }
1328 :
1329 : GEN
1330 49 : veceint1(GEN C, GEN nmax, long prec)
1331 : {
1332 49 : if (!nmax) return eint1(C,prec);
1333 7 : if (typ(nmax) != t_INT) pari_err_TYPE("veceint1",nmax);
1334 7 : if (typ(C) != t_REAL) {
1335 7 : C = gtofp(C, prec);
1336 7 : if (typ(C) != t_REAL) pari_err_TYPE("veceint1",C);
1337 : }
1338 7 : if (signe(C) <= 0) pari_err_DOMAIN("veceint1", "argument", "<=", gen_0,C);
1339 7 : return mpveceint1(C, NULL, itos(nmax));
1340 : }
1341 :
1342 : /* j > 0, a t_REAL. Return sum_{m >= 0} a^m / j(j+1)...(j+m)).
1343 : * Stop when expo(summand) < E; note that s(j-1) = (a s(j) + 1) / (j-1). */
1344 : static GEN
1345 231 : mp_sum_j(GEN a, long j, long E, long prec)
1346 : {
1347 231 : pari_sp av = avma;
1348 231 : GEN q = divru(real_1(prec), j), s = q;
1349 : long m;
1350 4290 : for (m = 0;; m++)
1351 : {
1352 4290 : if (expo(q) < E) break;
1353 4059 : q = mulrr(q, divru(a, m+j));
1354 4059 : s = addrr(s, q);
1355 : }
1356 231 : return gc_leaf(av, s);
1357 : }
1358 : /* Return the s_a(j), j <= J */
1359 : static GEN
1360 231 : sum_jall(GEN a, long J, long prec)
1361 : {
1362 231 : GEN s = cgetg(J+1, t_VEC);
1363 231 : long j, E = -prec2nbits(prec) - 5;
1364 231 : gel(s, J) = mp_sum_j(a, J, E, prec);
1365 9624 : for (j = J-1; j; j--)
1366 9393 : gel(s,j) = divru(addrs(mulrr(a, gel(s,j+1)), 1), j);
1367 231 : return s;
1368 : }
1369 :
1370 : /* T a dense t_POL with t_REAL coeffs. Return T(n) [faster than poleval] */
1371 : static GEN
1372 364903 : rX_s_eval(GEN T, long n)
1373 : {
1374 364903 : long i = lg(T)-1;
1375 364903 : GEN c = gel(T,i);
1376 7536888 : for (i--; i>=2; i--) c = gadd(mulrs(c,n),gel(T,i));
1377 364903 : return c;
1378 : }
1379 :
1380 : /* C>0 t_REAL, eC = exp(C). Return eint1(n*C) for 1<=n<=N. Absolute accuracy */
1381 : GEN
1382 238 : mpveceint1(GEN C, GEN eC, long N)
1383 : {
1384 238 : const long prec = realprec(C);
1385 238 : long Nmin = 15; /* >= 1. E.g. between 10 and 30, but little effect */
1386 238 : GEN en, v, w = cgetg(N+1, t_VEC);
1387 : pari_sp av0;
1388 : double DL;
1389 : long n, j, jmax, jmin;
1390 238 : if (!N) return w;
1391 368641 : for (n = 1; n <= N; n++) gel(w,n) = cgetr(prec);
1392 238 : av0 = avma;
1393 238 : if (N < Nmin) Nmin = N;
1394 238 : if (!eC) eC = mpexp(C);
1395 238 : en = eC; affrr(eint1p(C, en), gel(w,1));
1396 3500 : for (n = 2; n <= Nmin; n++)
1397 : {
1398 : pari_sp av2;
1399 3262 : en = mulrr(en,eC); /* exp(n C) */
1400 3262 : av2 = avma;
1401 3262 : affrr(eint1p(mulru(C,n), en), gel(w,n));
1402 3262 : set_avma(av2);
1403 : }
1404 238 : if (Nmin == N) { set_avma(av0); return w; }
1405 :
1406 231 : DL = prec2nbits_mul(prec, M_LN2) + 5;
1407 231 : jmin = ceil(DL/log((double)N)) + 1;
1408 231 : jmax = ceil(DL/log((double)Nmin)) + 1;
1409 231 : v = sum_jall(C, jmax, prec);
1410 231 : en = powrs(eC, -N); /* exp(-N C) */
1411 231 : affrr(eint1p(mulru(C,N), invr(en)), gel(w,N));
1412 6041 : for (j = jmin, n = N-1; j <= jmax; j++)
1413 : {
1414 5810 : long limN = maxss((long)ceil(exp(DL/j)), Nmin);
1415 : GEN polsh;
1416 5810 : setlg(v, j+1);
1417 5810 : polsh = RgV_to_RgX_reverse(v, 0);
1418 370713 : for (; n >= limN; n--)
1419 : {
1420 364903 : pari_sp av2 = avma;
1421 364903 : GEN S = divri(mulrr(en, rX_s_eval(polsh, -n)), powuu(n,j));
1422 : /* w[n+1] - exp(-n C) * polsh(-n) / (-n)^j */
1423 364903 : GEN c = odd(j)? addrr(gel(w,n+1), S) : subrr(gel(w,n+1), S);
1424 364903 : affrr(c, gel(w,n)); set_avma(av2);
1425 364903 : en = mulrr(en,eC); /* exp(-n C) */
1426 : }
1427 : }
1428 231 : set_avma(av0); return w;
1429 : }
1430 :
1431 : /* erfc via numerical integration : assume real(x)>=1 */
1432 : static GEN
1433 14 : cxerfc_r1(GEN x, long prec)
1434 : {
1435 : GEN h, h2, eh2, denom, res, lambda;
1436 : long u, v;
1437 14 : const double D = prec2nbits_mul(prec, M_LN2);
1438 14 : const long npoints = (long)ceil(D/M_PI)+1;
1439 14 : pari_sp av = avma;
1440 : {
1441 14 : double t = exp(-2*M_PI*M_PI/D); /* ~exp(-2*h^2) */
1442 14 : v = 30; /* bits that fit in both long and double mantissa */
1443 14 : u = (long)floor(t*(1L<<v));
1444 : /* define exp(-2*h^2) to be u*2^(-v) */
1445 : }
1446 14 : incrprec(prec);
1447 14 : x = gtofp(x,prec);
1448 14 : eh2 = sqrtr_abs(rtor(shiftr(dbltor(u),-v),prec));
1449 14 : h2 = negr(logr_abs(eh2));
1450 14 : h = sqrtr_abs(h2);
1451 14 : lambda = gdiv(x,h);
1452 14 : denom = gsqr(lambda);
1453 : { /* res = h/x + 2*x*h*sum(k=1,npoints,exp(-(k*h)^2)/(lambda^2+k^2)); */
1454 : GEN Uk; /* = exp(-(kh)^2) */
1455 14 : GEN Vk = eh2;/* = exp(-(2k+1)h^2) */
1456 14 : pari_sp av2 = avma;
1457 : long k;
1458 : /* k = 0 moved out for efficiency */
1459 14 : denom = gaddsg(1,denom);
1460 14 : Uk = Vk;
1461 14 : Vk = mulur(u,Vk); shiftr_inplace(Vk, -v);
1462 14 : res = gdiv(Uk, denom);
1463 420 : for (k = 1; k < npoints; k++)
1464 : {
1465 406 : if ((k & 255) == 0) (void)gc_all(av2,4,&denom,&Uk,&Vk,&res);
1466 406 : denom = gaddsg(2*k+1,denom);
1467 406 : Uk = mpmul(Uk,Vk);
1468 406 : Vk = mulur(u,Vk); shiftr_inplace(Vk, -v);
1469 406 : res = gadd(res, gdiv(Uk, denom));
1470 : }
1471 : }
1472 14 : res = gmul(res, gshift(lambda,1));
1473 : /* 0 term : */
1474 14 : res = gadd(res, ginv(lambda));
1475 14 : res = gmul(res, gdiv(gexp(gneg(gsqr(x)), prec), mppi(prec)));
1476 14 : if (rtodbl(real_i(x)) < sqrt(D))
1477 : {
1478 14 : GEN t = gmul(divrr(Pi2n(1,prec),h), x);
1479 14 : res = gsub(res, gdivsg(2, cxexpm1(t, prec)));
1480 : }
1481 14 : return gc_upto(av,res);
1482 : }
1483 :
1484 : static GEN
1485 7 : sererfc(GEN x, long prec)
1486 : {
1487 7 : GEN u, z = invr(sqrtr_abs(Pi2n(-2,prec)));
1488 7 : setsigne(z, -1); /* -2/sqrt(Pi) */
1489 7 : z = gmul(z, integser(gmul(derivser(x), gexp(gneg(gsqr(x)), prec))));
1490 7 : u = polcoef_i(x, 0, varn(x));
1491 7 : if (!gequal0(u)) z = gadd(z, gerfc(u, prec));
1492 7 : return z;
1493 : }
1494 :
1495 : GEN
1496 70 : gerfc(GEN x, long prec)
1497 : {
1498 : GEN z, xr, xi, res;
1499 : long s;
1500 : pari_sp av;
1501 :
1502 70 : switch(typ(x))
1503 : {
1504 63 : case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX:
1505 63 : break;
1506 7 : default:
1507 7 : av = avma;
1508 7 : if ((z = toser_i(x))) return gc_upto(av, sererfc(z,prec));
1509 0 : return trans_eval("erfc",gerfc,x,prec);
1510 : }
1511 : /* x a complex scalar */
1512 63 : x = trans_fix_arg(&prec,&x,&xr,&xi, &av,&res);
1513 63 : s = signe(xr);
1514 63 : if (s > 0 || (s == 0 && signe(xi) >= 0)) {
1515 98 : if (cmprs(xr, 1) > 0) /* use numerical integration */
1516 14 : z = cxerfc_r1(x, prec);
1517 : else
1518 : { /* erfc(x) = incgam(1/2,x^2)/sqrt(Pi) */
1519 35 : GEN sqrtpi = sqrtr(mppi(prec));
1520 35 : z = incgam0(ghalf, gsqr(x), sqrtpi, prec);
1521 35 : z = gdiv(z, sqrtpi);
1522 : }
1523 : }
1524 : else
1525 : { /* erfc(-x)=2-erfc(x) */
1526 : /* FIXME could decrease prec
1527 : long size = nbits2extraprec((imag(x)^2-real(x)^2)/log(2));
1528 : prec = size > 0 ? prec : prec + size;
1529 : */
1530 : /* NOT gsubsg(2, ...) : would create a result of
1531 : * huge accuracy if re(x)>>1, rounded to 2 by subsequent affc_fixlg... */
1532 14 : z = gsub(real2n(1,prec+EXTRAPREC64), gerfc(gneg(x), prec));
1533 : }
1534 63 : set_avma(av); return affc_fixlg(z, res);
1535 : }
1536 :
1537 : /***********************************************************************/
1538 : /** **/
1539 : /** RIEMANN ZETA FUNCTION **/
1540 : /** **/
1541 : /***********************************************************************/
1542 : static const double log2PI = 1.83787706641;
1543 :
1544 : static double
1545 4585 : get_xinf(double beta)
1546 : {
1547 4585 : const double maxbeta = 0.06415003; /* 3^(-2.5) */
1548 : double x0, y0, x1;
1549 :
1550 4585 : if (beta < maxbeta) return beta + pow(3*beta, 1.0/3.0);
1551 4585 : x0 = beta + M_PI/2.0;
1552 : for(;;)
1553 : {
1554 7539 : y0 = x0*x0;
1555 7539 : x1 = (beta+atan(x0)) * (1+y0) / y0 - 1/x0;
1556 7539 : if (0.99*x0 < x1) return x1;
1557 2954 : x0 = x1;
1558 : }
1559 : }
1560 : /* optimize for zeta( s + it, prec ), assume |s-1| > 0.1
1561 : * (if gexpo(u = s-1) < -5, we use the functional equation s->1-s) */
1562 : static int
1563 18494 : optim_zeta(GEN S, long prec, long *pp, long *pn)
1564 : {
1565 : double s, t, alpha, beta, n, B;
1566 : long p;
1567 18494 : if (typ(S) == t_REAL) {
1568 10129 : t = 0.;
1569 10129 : s = rtodbl(S);
1570 : } else {
1571 8365 : t = fabs( rtodbl(gel(S,2)) );
1572 8365 : if (t > 2500) return 0; /* lfunlarge */
1573 8316 : s = rtodbl(gel(S,1));
1574 : }
1575 :
1576 18445 : B = prec2nbits_mul(prec, M_LN2);
1577 18445 : if (s > 0 && !t) /* positive real input */
1578 : {
1579 10080 : beta = B + 0.61 + s*(log2PI - log(s));
1580 10080 : if (beta > 0)
1581 : {
1582 5110 : p = (long)ceil(beta / 2.0);
1583 5110 : n = fabs(s + 2*p-1)/(2*M_PI);
1584 : }
1585 : else
1586 : {
1587 4970 : p = 0;
1588 4970 : n = exp((B - M_LN2) / s);
1589 : }
1590 : }
1591 8365 : else if (s <= 0 || t < 0.01) /* s < 0 may occur if s ~ 0 */
1592 3773 : { /* TODO: the crude bounds below are generally valid. Optimize ? */
1593 3773 : double l,l2, la = 1.; /* heuristic */
1594 3773 : double rlog, ilog; dblclog(s-1,t, &rlog,&ilog);
1595 3773 : l2 = (s - 0.5)*rlog - t*ilog; /* = Re( (S - 1/2) log (S-1) ) */
1596 3773 : l = (B - l2 + s*log2PI) / (2. * (1.+ log((double)la)));
1597 3773 : l2 = dblcabs(s, t)/2;
1598 3773 : if (l < l2) l = l2;
1599 3773 : p = (long) ceil(l); if (p < 2) p = 2;
1600 3773 : n = 1 + dblcabs(p+s/2.-.25, t/2) * la / M_PI;
1601 : }
1602 : else
1603 : {
1604 4592 : double sn = dblcabs(s, t), L = log(sn/s);
1605 4592 : alpha = B - 0.39 + L + s*(log2PI - log(sn));
1606 4592 : beta = (alpha+s)/t - atan(s/t);
1607 4592 : p = 0;
1608 4592 : if (beta > 0)
1609 : {
1610 4585 : beta = 1.0 - s + t * get_xinf(beta);
1611 4585 : if (beta > 0) p = (long)ceil(beta / 2.0);
1612 : }
1613 : else
1614 7 : if (s < 1.0) p = 1;
1615 4592 : n = p? dblcabs(s + 2*p-1, t) / (2*M_PI) : exp((B-M_LN2+L) / s);
1616 : }
1617 18445 : *pp = p;
1618 18445 : *pn = (long)ceil(n);
1619 18445 : if (*pp < 0 || *pn < 0) pari_err_OVERFLOW("zeta");
1620 18445 : return 1;
1621 : }
1622 :
1623 : /* zeta(a*j+b), j=0..N-1, b>1, using sumalt. Johansonn's thesis, Algo 4.7.1 */
1624 : static GEN
1625 711 : veczetas(long a, long b, long N, long prec)
1626 : {
1627 711 : const long n = ceil(2 + prec2nbits_mul(prec, M_LN2/1.7627));
1628 711 : pari_sp av = avma;
1629 711 : GEN c, d, z = zerovec(N);
1630 : long j, k;
1631 712 : c = d = int2n(2*n-1);
1632 88875 : for (k = n; k > 1; k--)
1633 : {
1634 88163 : GEN u, t = divii(d, powuu(k,b));
1635 88162 : if (!odd(k)) t = negi(t);
1636 88205 : gel(z,1) = addii(gel(z,1), t);
1637 88233 : u = powuu(k,a);
1638 4113974 : for (j = 1; j < N; j++)
1639 : {
1640 4069691 : t = divii(t,u); if (!signe(t)) break;
1641 4021665 : gel(z,j+1) = addii(gel(z,j+1), t);
1642 : }
1643 89364 : c = muluui(k,2*k-1,c);
1644 88218 : c = diviuuexact(c, 2*(n-k+1),n+k-1);
1645 88186 : d = addii(d,c);
1646 88162 : if (gc_needed(av,3))
1647 : {
1648 7 : if(DEBUGMEM>1) pari_warn(warnmem,"zetaBorwein, k = %ld", k);
1649 7 : (void)gc_all(av, 3, &c,&d,&z);
1650 : }
1651 : }
1652 : /* k = 1 */
1653 53672 : for (j = 1; j <= N; j++) gel(z,j) = addii(gel(z,j), d);
1654 713 : d = addiu(d, 1);
1655 53702 : for (j = 0, k = b - 1; j < N; j++, k += a)
1656 52990 : gel(z,j+1) = rdivii(shifti(gel(z,j+1), k), subii(shifti(d,k), d), prec);
1657 712 : return z;
1658 : }
1659 : /* zeta(a*j+b), j=0..N-1, b > 1, a*(N-1) + b > 1, using sumalt.
1660 : * a <= 0 is allowed (including the silly a = 0) */
1661 : GEN
1662 224 : veczeta(GEN a, GEN b, long N, long prec)
1663 : {
1664 224 : pari_sp av = avma;
1665 : long n, j, k;
1666 : GEN L, c, d, z;
1667 224 : if (typ(a) == t_INT && typ(b) == t_INT)
1668 126 : return gc_GEN(av, veczetas(itos(a), itos(b), N, prec));
1669 98 : z = zerovec(N);
1670 98 : n = ceil(2 + prec2nbits_mul(prec, M_LN2/1.7627));
1671 98 : c = d = int2n(2*n-1);
1672 14197 : for (k = n; k; k--)
1673 : {
1674 : GEN u, t;
1675 14099 : L = logr_abs(utor(k, prec)); /* log(k) */
1676 14099 : t = gdiv(d, gexp(gmul(b, L), prec)); /* d / k^b */
1677 14099 : if (!odd(k)) t = gneg(t);
1678 14099 : gel(z,1) = gadd(gel(z,1), t);
1679 14099 : u = gexp(gmul(a, L), prec);
1680 898001 : for (j = 1; j < N; j++)
1681 : {
1682 890325 : t = gdiv(t,u); if (gexpo(t) < 0) break;
1683 883902 : gel(z,j+1) = gadd(gel(z,j+1), t);
1684 : }
1685 14099 : c = muluui(k,2*k-1,c);
1686 14099 : c = diviuuexact(c, 2*(n-k+1),n+k-1);
1687 14099 : d = addii(d,c);
1688 14099 : if (gc_needed(av,3))
1689 : {
1690 0 : if(DEBUGMEM>1) pari_warn(warnmem,"veczeta, k = %ld", k);
1691 0 : (void)gc_all(av, 3, &c,&d,&z);
1692 : }
1693 : }
1694 98 : L = mplog2(prec);
1695 5824 : for (j = 0; j < N; j++)
1696 : {
1697 5726 : GEN u = gsubgs(gadd(b, gmulgu(a,j)), 1);
1698 5726 : GEN w = gexp(gmul(u, L), prec); /* 2^u */
1699 5726 : gel(z,j+1) = gdiv(gmul(gel(z,j+1), w), gmul(d,gsubgs(w,1)));
1700 : }
1701 98 : return gc_GEN(av, z);
1702 : }
1703 :
1704 : GEN
1705 29496 : constzeta(long n, long prec)
1706 : {
1707 29496 : GEN o = zetazone, z;
1708 29496 : long l = o? lg(o): 0;
1709 : pari_sp av;
1710 29496 : if (l > n)
1711 : {
1712 29076 : long p = realprec(gel(o,1));
1713 29076 : if (p >= prec) return o;
1714 : }
1715 585 : n = maxss(n, l + 15);
1716 585 : av = avma; z = veczetas(1, 2, n-1, prec);
1717 586 : zetazone = gclone(vec_prepend(z, mpeuler(prec)));
1718 586 : set_avma(av); guncloneNULL(o); return zetazone;
1719 : }
1720 :
1721 : /* zeta(s) using sumalt, case h=0,N=1. Assume s > 1 */
1722 : static GEN
1723 598 : zetaBorwein(long s, long prec)
1724 : {
1725 598 : pari_sp av = avma;
1726 598 : const long n = ceil(2 + prec2nbits_mul(prec, M_LN2/1.7627));
1727 : long k;
1728 598 : GEN c, d, z = gen_0;
1729 598 : c = d = int2n(2*n-1);
1730 135672 : for (k = n; k; k--)
1731 : {
1732 135074 : GEN t = divii(d, powuu(k,s));
1733 135074 : z = odd(k)? addii(z,t): subii(z,t);
1734 135074 : c = muluui(k,2*k-1,c);
1735 135074 : c = diviuuexact(c, 2*(n-k+1),n+k-1);
1736 135074 : d = addii(d,c);
1737 135074 : if (gc_needed(av,3))
1738 : {
1739 0 : if(DEBUGMEM>1) pari_warn(warnmem,"zetaBorwein, k = %ld", k);
1740 0 : (void)gc_all(av, 3, &c,&d,&z);
1741 : }
1742 : }
1743 598 : return rdivii(shifti(z, s-1), subii(shifti(d,s-1), d), prec);
1744 : }
1745 :
1746 : /* assume k != 1 */
1747 : GEN
1748 4760 : szeta(long k, long prec)
1749 : {
1750 4760 : pari_sp av = avma;
1751 : GEN z;
1752 :
1753 4760 : if (!k) { z = real2n(-1, prec); setsigne(z,-1); return z; }
1754 4753 : if (k < 0)
1755 : {
1756 0 : if (!odd(k)) return gen_0;
1757 : /* the one value such that k < 0 and 1 - k < 0, due to overflow */
1758 0 : if ((ulong)k == (HIGHBIT | 1))
1759 0 : pari_err_OVERFLOW("zeta [large negative argument]");
1760 0 : k = 1-k;
1761 0 : z = bernreal(k, prec); togglesign(z);
1762 0 : return gc_leaf(av, divru(z, k));
1763 : }
1764 : /* k > 1 */
1765 4753 : if (k > prec2nbits(prec)+1) return real_1(prec);
1766 4746 : if (zetazone && realprec(gel(zetazone,1)) >= prec && lg(zetazone) > k)
1767 140 : return rtor(gel(zetazone, k), prec);
1768 4606 : if (!odd(k))
1769 : {
1770 : GEN B;
1771 2835 : if (!bernzone) constbern(0);
1772 2835 : if (k < lg(bernzone))
1773 2380 : B = gel(bernzone, k>>1);
1774 : else
1775 : {
1776 455 : if (bernbitprec(k) > prec2nbits(prec))
1777 0 : return gc_upto(av, invr(inv_szeta_euler(k, prec)));
1778 455 : B = bernfrac(k);
1779 : }
1780 : /* B = B_k */
1781 2835 : z = gmul(powru(Pi2n(1, prec + EXTRAPREC64), k), B);
1782 2835 : z = k < 410? rtor(divri(z, mpfact(k)), prec)
1783 2835 : : divrr(z, mpfactr(k,prec));
1784 2835 : setsigne(z, 1); shiftr_inplace(z, -1);
1785 : }
1786 : else
1787 : {
1788 1771 : double p = prec2nbits_mul(prec,0.393); /* bit / log_2(3+sqrt(8)) */
1789 1771 : p = log2(p * log(p));
1790 3542 : z = (p * k > prec2nbits(prec))? invr(inv_szeta_euler(k, prec))
1791 1771 : : zetaBorwein(k, prec);
1792 : }
1793 4606 : return gc_leaf(av, z);
1794 : }
1795 :
1796 : /* Ensure |1-s| >= 1/32 and (|s| <= 1/32 or real(s) >= 1/2) */
1797 : static int
1798 18508 : zeta_funeq(GEN *ps)
1799 : {
1800 18508 : GEN s = *ps, u;
1801 18508 : if (typ(s) == t_REAL)
1802 : {
1803 10136 : u = subsr(1, s);
1804 10136 : if (expo(u) >= -5
1805 10108 : && ((signe(s) > 0 && expo(s) >= -1) || expo(s) <= -5)) return 0;
1806 : }
1807 : else
1808 : {
1809 8372 : GEN sig = gel(s,1);
1810 8372 : if (fabs(rtodbl(gel(s,2))) > 2500) return 0; /* lfunlarge */
1811 8323 : u = gsubsg(1, s);
1812 8323 : if (gexpo(u) >= -5
1813 8316 : && ((signe(sig) > 0 && expo(sig) >= -1) || gexpo(s) <= -5)) return 0;
1814 : }
1815 1526 : *ps = u; return 1;
1816 : }
1817 : /* s0 a t_INT, t_REAL or t_COMPLEX.
1818 : * If a t_INT, assume it's not a trivial case (i.e we have s0 > 1, odd) */
1819 : static GEN
1820 18515 : czeta(GEN s0, long prec)
1821 : {
1822 18515 : GEN ms, s, u, y, res, tes, sig, tau, invn2, ns, Ns, funeq_factor = NULL;
1823 : long i, nn, lim, lim2;
1824 18515 : pari_sp av0 = avma, av, av2;
1825 : pari_timer T;
1826 :
1827 18515 : if (DEBUGLEVEL>2) timer_start(&T);
1828 18515 : s = trans_fix_arg(&prec,&s0,&sig,&tau,&av,&res);
1829 18515 : if (typ(s0) == t_INT) return gc_upto(av0, gzeta(s0, prec));
1830 18508 : if (zeta_funeq(&s)) /* s -> 1-s */
1831 : { /* Gamma(s) (2Pi)^-s 2 cos(Pi s/2) [new s] */
1832 1526 : GEN t = gmul(ggamma(s,prec), pow2Pis(gsubgs(s0,1), prec));
1833 1526 : sig = real_i(s);
1834 1526 : funeq_factor = gmul2n(gmul(t, gsin(gmul(Pi2n(-1,prec),s0), prec)), 1);
1835 : }
1836 18508 : if (gcmpgs(sig, prec2nbits(prec) + 1) > 0) { /* zeta(s) = 1 */
1837 14 : if (!funeq_factor) { set_avma(av0); return real_1(prec); }
1838 7 : return gc_upto(av0, funeq_factor);
1839 : }
1840 18494 : if (!optim_zeta(s, prec, &lim, &nn))
1841 : {
1842 49 : long bit = prec2nbits(prec);
1843 49 : y = lfun(lfuninit(gen_1, cgetg(1,t_VEC), 0, bit), s, bit);
1844 49 : if (funeq_factor) y = gmul(y, funeq_factor);
1845 49 : set_avma(av); return affc_fixlg(y,res);
1846 : }
1847 18445 : if (DEBUGLEVEL>2) err_printf("lim, nn: [%ld, %ld]\n", lim, nn);
1848 18445 : ms = gneg(s);
1849 18445 : if (umuluu_le(nn, prec, 10000000))
1850 : {
1851 18445 : incrprec(prec); /* one extra word of precision */
1852 18445 : Ns = vecpowug(nn, ms, prec);
1853 18445 : ns = gel(Ns,nn); setlg(Ns, nn);
1854 18445 : y = gadd(gmul2n(ns, -1), RgV_sum(Ns));
1855 : }
1856 : else
1857 : {
1858 0 : Ns = dirpowerssum(nn, ms, 0, prec);
1859 0 : incrprec(prec); /* one extra word of precision */
1860 0 : ns = gpow(utor(nn, prec), ms, prec);
1861 0 : y = gsub(Ns, gmul2n(ns, -1));
1862 : }
1863 18445 : if (DEBUGLEVEL>2) timer_printf(&T,"sum from 1 to N");
1864 18445 : constbern(lim);
1865 18445 : if (DEBUGLEVEL>2) timer_start(&T);
1866 18445 : invn2 = divri(real_1(prec), sqru(nn)); lim2 = lim<<1;
1867 18445 : tes = bernfrac(lim2);
1868 : {
1869 : GEN s1, s2, s3, s4, s5;
1870 18445 : s2 = gmul(s, gsubgs(s,1));
1871 18445 : s3 = gmul2n(invn2,3);
1872 18445 : av2 = avma;
1873 18445 : s1 = gsubgs(gmul2n(s,1), 1);
1874 18445 : s4 = gmul(invn2, gmul2n(gaddsg(4*lim-2,s1),1));
1875 18445 : s5 = gmul(invn2, gadd(s2, gmulsg(lim2, gaddgs(s1, lim2))));
1876 1584104 : for (i = lim2-2; i>=2; i -= 2)
1877 : {
1878 1565659 : s5 = gsub(s5, s4);
1879 1565659 : s4 = gsub(s4, s3);
1880 1565659 : tes = gadd(bernfrac(i), gdivgunextu(gmul(s5,tes), i+1));
1881 1565659 : if (gc_needed(av2,3))
1882 : {
1883 0 : if(DEBUGMEM>1) pari_warn(warnmem,"czeta i = %ld", i);
1884 0 : (void)gc_all(av2,3, &tes,&s5,&s4);
1885 : }
1886 : }
1887 18445 : u = gmul(gmul(tes,invn2), gmul2n(s2, -1));
1888 18445 : tes = gmulsg(nn, gaddsg(1, u));
1889 : }
1890 18445 : if (DEBUGLEVEL>2) timer_printf(&T,"Bernoulli sum");
1891 : /* y += tes n^(-s) / (s-1) */
1892 18445 : y = gadd(y, gmul(tes, gdiv(ns, gsubgs(s,1))));
1893 18445 : if (funeq_factor) y = gmul(y, funeq_factor);
1894 18445 : set_avma(av); return affc_fixlg(y,res);
1895 : }
1896 : /* v a t_VEC/t_COL; is v[i] = a + b (i-1) for some a,b ? */
1897 : int
1898 42 : RgV_is_arithprog(GEN v, GEN *a, GEN *b)
1899 : {
1900 42 : pari_sp av = avma, av2;
1901 42 : long i, n = lg(v)-1;
1902 42 : if (n == 0) { *a = *b = gen_0; return 1; }
1903 42 : *a = gel(v,1);
1904 42 : if (n == 1) { * b = gen_0; return 1; }
1905 42 : *b = gsub(gel(v,2), *a); av2 = avma;
1906 77 : for (i = 2; i < n; i++)
1907 35 : if (!gequal(*b, gsub(gel(v,i+1), gel(v,i)))) return gc_int(av,0);
1908 42 : return gc_int(av2,1);
1909 : }
1910 :
1911 : GEN
1912 23583 : gzeta(GEN x, long prec)
1913 : {
1914 23583 : pari_sp av = avma;
1915 : GEN y;
1916 23583 : if (gequal1(x)) pari_err_DOMAIN("zeta", "argument", "=", gen_1, x);
1917 23548 : switch(typ(x))
1918 : {
1919 4774 : case t_INT:
1920 4774 : if (is_bigint(x))
1921 : {
1922 21 : if (signe(x) > 0) return real_1(prec);
1923 14 : if (mod2(x) == 0) return real_0(prec);
1924 7 : pari_err_OVERFLOW("zeta [large negative argument]");
1925 : }
1926 4753 : return szeta(itos(x),prec);
1927 18515 : case t_REAL: case t_COMPLEX: return czeta(x,prec);
1928 14 : case t_PADIC: return Qp_zeta(x);
1929 49 : case t_VEC: case t_COL:
1930 : {
1931 : GEN a, b;
1932 49 : long n = lg(x) - 1;
1933 49 : if (n > 1 && RgV_is_arithprog(x, &b, &a))
1934 : {
1935 42 : if (!is_real_t(typ(a)) || !is_real_t(typ(b))
1936 35 : || gcmpgs(gel(x,1), 1) <= 0
1937 42 : || gcmpgs(gel(x,n), 1) <= 0) { set_avma(av); break; }
1938 21 : a = veczeta(a, b, n, prec);
1939 21 : settyp(a, typ(x)); return a;
1940 : }
1941 : }
1942 : default:
1943 203 : if (!(y = toser_i(x))) break;
1944 35 : if (gequal1(y))
1945 7 : pari_err_DOMAIN("zeta", "argument", "=", gen_1, y);
1946 28 : return gc_upto(av, lfun(gen_1,y,prec2nbits(prec)));
1947 : }
1948 189 : return trans_eval("zeta",gzeta,x,prec);
1949 : }
1950 :
1951 : /***********************************************************************/
1952 : /** **/
1953 : /** FONCTIONS POLYLOGARITHME **/
1954 : /** **/
1955 : /***********************************************************************/
1956 :
1957 : /* smallish k such that bernbitprec(K) > bit + Kdz, K = 2k+4 */
1958 : static long
1959 21 : get_k(double dz, long bit)
1960 : {
1961 : long a, b;
1962 21 : for (b = 128;; b <<= 1)
1963 21 : if (bernbitprec(b) > bit + b*dz) break;
1964 21 : if (b == 128) return 128;
1965 0 : a = b >> 1;
1966 0 : while (b - a > 64)
1967 : {
1968 0 : long c = (a+b) >> 1;
1969 0 : if (bernbitprec(c) > bit + c*dz) b = c; else a = c;
1970 : }
1971 0 : return b >> 1;
1972 : }
1973 :
1974 : /* m >= 2. Validity domain |log x| < 2*Pi, contains log |x| < 5.44,
1975 : * Li_m(x = e^z) = sum_{n >= 0} zeta(m-n) z^n / n!
1976 : * with zeta(1) := H_{m-1} - log(-z) */
1977 : static GEN
1978 21 : cxpolylog(long m, GEN x, long prec)
1979 : {
1980 : long li, n, k, ksmall, real;
1981 : GEN vz, z, Z, h, q, s, S;
1982 : pari_sp av;
1983 : double dz;
1984 : pari_timer T;
1985 :
1986 21 : if (gequal1(x)) return szeta(m,prec);
1987 : /* x real <= 1 ==> Li_m(x) real */
1988 21 : real = (typ(x) == t_REAL && (expo(x) < 0 || signe(x) <= 0));
1989 :
1990 21 : vz = constzeta(m, prec);
1991 21 : z = glog(x,prec);
1992 : /* n = 0 */
1993 21 : q = gen_1; s = gel(vz, m);
1994 28 : for (n=1; n < m-1; n++)
1995 : {
1996 7 : q = gdivgu(gmul(q,z),n);
1997 7 : s = gadd(s, gmul(gel(vz,m-n), real? real_i(q): q));
1998 : }
1999 : /* n = m-1 */
2000 21 : q = gdivgu(gmul(q,z),n); /* multiply by "zeta(1)" */
2001 21 : h = gmul(q, gsub(harmonic(m-1), glog(gneg_i(z),prec)));
2002 21 : s = gadd(s, real? real_i(h): h);
2003 : /* n = m */
2004 21 : q = gdivgu(gmul(q,z),m);
2005 21 : s = gadd(s, gdivgs(real? real_i(q): q, -2)); /* zeta(0) = -1/2 */
2006 : /* n = m+1 */
2007 21 : q = gdivgu(gmul(q,z),m+1); /* = z^(m+1) / (m+1)! */
2008 21 : s = gadd(s, gdivgs(real? real_i(q): q, -12)); /* zeta(-1) = -1/12 */
2009 :
2010 21 : li = -(prec2nbits(prec)+1);
2011 21 : if (DEBUGLEVEL) timer_start(&T);
2012 21 : dz = dbllog2(z) - log2PI; /* ~ log2(|z|/2Pi) */
2013 : /* sum_{k >= 1} zeta(-1-2k) * z^(2k+m+1) / (2k+m+1)!
2014 : * = 2 z^(m-1) sum_{k >= 1} zeta(2k+2) * Z^(k+1) / (2k+2)..(2k+1+m)), where
2015 : * Z = -(z/2Pi)^2. Stop at 2k = (li - (m-1)*Lz - m) / dz, Lz = log2 |z| */
2016 : /* We cut the sum in two: small values of k first */
2017 21 : Z = gsqr(z); av = avma;
2018 21 : ksmall = get_k(dz, prec2nbits(prec));
2019 21 : constbern(ksmall);
2020 469 : for(k = 1; k < ksmall; k++)
2021 : {
2022 469 : GEN t = q = gdivgunextu(gmul(q,Z), 2*k+m); /* z^(2k+m+1)/(2k+m+1)! */
2023 469 : if (real) t = real_i(t);
2024 469 : t = gmul(t, gdivgu(bernfrac(2*k+2), 2*k+2)); /* - t * zeta(1-(2k+2)) */
2025 469 : s = gsub(s, t);
2026 469 : if (gexpo(t) < li) return s;
2027 : /* large values ? */
2028 448 : if ((k & 0x1ff) == 0) (void)gc_all(av, 2, &s, &q);
2029 : }
2030 0 : if (DEBUGLEVEL>2) timer_printf(&T, "polylog: small k <= %ld", k);
2031 0 : Z = gneg(gsqr(gdiv(z, Pi2n(1,prec))));
2032 0 : q = gmul(gpowgs(z, m-1), gpowgs(Z, k+1)); /* Z^(k+1) * z^(m-1) */
2033 0 : S = gen_0; av = avma;
2034 0 : for(;; k++)
2035 0 : {
2036 0 : GEN t = q;
2037 : long b;
2038 0 : if (real) t = real_i(t);
2039 0 : b = prec + gexpo(t) / BITS_IN_LONG; /* decrease accuracy */
2040 0 : if (b == 2) break;
2041 : /* t * zeta(2k+2) / (2k+2)..(2k+1+m) */
2042 0 : t = gdiv(t, mulri(inv_szeta_euler(2*k+2, b),
2043 0 : mulu_interval(2*k+2, 2*k+1+m)));
2044 0 : S = gadd(S, t); if (gexpo(t) < li) break;
2045 0 : q = gmul(q, Z);
2046 0 : if ((k & 0x1ff) == 0) (void)gc_all(av, 2, &S, &q);
2047 : }
2048 0 : if (DEBUGLEVEL>2) timer_printf(&T, "polylog: large k <= %ld", k);
2049 0 : return gadd(s, gmul2n(S,1));
2050 : }
2051 :
2052 : static GEN
2053 42 : Li1(GEN x, long prec) { return gneg(glog(gsubsg(1, x), prec)); }
2054 :
2055 : static GEN
2056 203 : polylog(long m, GEN x, long prec)
2057 : {
2058 : long l, e, i, G, sx;
2059 : pari_sp av, av1;
2060 : GEN X, Xn, z, p1, p2, y, res;
2061 :
2062 203 : if (m < 0) pari_err_DOMAIN("polylog", "index", "<", gen_0, stoi(m));
2063 203 : if (!m) return mkfrac(gen_m1,gen_2);
2064 203 : if (gequal0(x)) return gcopy(x);
2065 203 : if (m==1) { av = avma; return gc_upto(av, Li1(x, prec)); }
2066 :
2067 168 : l = precision(x);
2068 168 : if (!l) l = prec; else prec = l;
2069 168 : res = cgetc(l); av = avma;
2070 168 : x = gtofp(x, l+EXTRAPREC64);
2071 168 : e = gexpo(gnorm(x));
2072 168 : if (!e || e == -1) {
2073 21 : y = cxpolylog(m,x,prec);
2074 21 : set_avma(av); return affc_fixlg(y, res);
2075 : }
2076 147 : X = (e > 0)? ginv(x): x;
2077 147 : G = -prec2nbits(l);
2078 147 : av1 = avma;
2079 147 : y = Xn = X;
2080 147 : for (i=2; ; i++)
2081 : {
2082 68159 : Xn = gmul(X,Xn); p2 = gdiv(Xn,powuu(i,m));
2083 68159 : y = gadd(y,p2);
2084 68159 : if (gexpo(p2) <= G) break;
2085 :
2086 68012 : if (gc_needed(av1,1))
2087 : {
2088 0 : if(DEBUGMEM>1) pari_warn(warnmem,"polylog");
2089 0 : (void)gc_all(av1,2, &y, &Xn);
2090 : }
2091 : }
2092 147 : if (e < 0) { set_avma(av); return affc_fixlg(y, res); }
2093 :
2094 28 : sx = gsigne(imag_i(x));
2095 28 : if (!sx)
2096 : {
2097 28 : if (m&1) sx = gsigne(gsub(gen_1, real_i(x)));
2098 21 : else sx = - gsigne(real_i(x));
2099 : }
2100 28 : z = divri(mppi(l), mpfact(m-1)); setsigne(z, sx);
2101 28 : z = mkcomplex(gen_0, z);
2102 :
2103 28 : if (m == 2)
2104 : { /* same formula as below, written more efficiently */
2105 21 : y = gneg_i(y);
2106 21 : if (typ(x) == t_REAL && signe(x) < 0)
2107 7 : p1 = logr_abs(x);
2108 : else
2109 14 : p1 = gsub(glog(x,l), z);
2110 21 : p1 = gmul2n(gsqr(p1), -1); /* = (log(-x))^2 / 2 */
2111 :
2112 21 : p1 = gadd(p1, divru(sqrr(mppi(l)), 6));
2113 21 : p1 = gneg_i(p1);
2114 : }
2115 : else
2116 : {
2117 7 : GEN logx = glog(x,l), logx2 = gsqr(logx), vz = constzeta(m, l);
2118 7 : p1 = mkfrac(gen_m1,gen_2);
2119 14 : for (i = m-2; i >= 0; i -= 2)
2120 7 : p1 = gadd(gel(vz, m-i), gmul(p1, gdivgunextu(logx2, i+1)));
2121 7 : if (m&1) p1 = gmul(logx,p1); else y = gneg_i(y);
2122 7 : p1 = gadd(gmul2n(p1,1), gmul(z,gpowgs(logx,m-1)));
2123 7 : if (typ(x) == t_REAL && signe(x) < 0) p1 = real_i(p1);
2124 : }
2125 28 : y = gadd(y,p1);
2126 28 : set_avma(av); return affc_fixlg(y, res);
2127 : }
2128 : static GEN
2129 119 : RIpolylog(long m, GEN x, long real, long prec)
2130 : {
2131 119 : GEN y = polylog(m, x, prec);
2132 119 : return real? real_i(y): imag_i(y);
2133 : }
2134 : GEN
2135 21 : dilog(GEN x, long prec) { return gpolylog(2, x, prec); }
2136 :
2137 : /* x a floating point number, t_REAL or t_COMPLEX of t_REAL */
2138 : static GEN
2139 42 : logabs(GEN x)
2140 : {
2141 : GEN y;
2142 42 : if (typ(x) == t_COMPLEX)
2143 : {
2144 7 : y = logr_abs( cxnorm(x) );
2145 7 : shiftr_inplace(y, -1);
2146 : } else
2147 35 : y = logr_abs(x);
2148 42 : return y;
2149 : }
2150 :
2151 : static GEN
2152 21 : polylogD(long m, GEN x, long flag, long prec)
2153 : {
2154 21 : long fl = 0, k, l, m2;
2155 : pari_sp av;
2156 : GEN p1, p2, y;
2157 :
2158 21 : if (gequal0(x)) return gcopy(x);
2159 21 : m2 = m&1;
2160 21 : if (gequal1(x) && m>=2) return m2? szeta(m,prec): gen_0;
2161 21 : av = avma; l = precision(x);
2162 21 : if (!l) { l = prec; x = gtofp(x,l); }
2163 21 : p1 = logabs(x);
2164 21 : if (signe(p1) > 0) { x = ginv(x); fl = !m2; } else setabssign(p1);
2165 : /* |x| <= 1, p1 = - log|x| >= 0 */
2166 21 : p2 = gen_1;
2167 21 : y = RIpolylog(m, x, m2, l);
2168 84 : for (k = 1; k < m; k++)
2169 : {
2170 63 : GEN t = RIpolylog(m-k, x, m2, l);
2171 63 : p2 = gdivgu(gmul(p2,p1), k); /* (-log|x|)^k / k! */
2172 63 : y = gadd(y, gmul(p2, t));
2173 : }
2174 21 : if (m2)
2175 : {
2176 14 : p1 = flag? gdivgs(p1, -2*m): gdivgs(logabs(gsubsg(1,x)), m);
2177 14 : y = gadd(y, gmul(p2, p1));
2178 : }
2179 21 : if (fl) y = gneg(y);
2180 21 : return gc_upto(av, y);
2181 : }
2182 :
2183 : static GEN
2184 14 : polylogP(long m, GEN x, long prec)
2185 : {
2186 14 : long fl = 0, k, l, m2;
2187 : pari_sp av;
2188 : GEN p1,y;
2189 :
2190 14 : if (gequal0(x)) return gcopy(x);
2191 14 : m2 = m&1;
2192 14 : if (gequal1(x) && m>=2) return m2? szeta(m,prec): gen_0;
2193 14 : av = avma; l = precision(x);
2194 14 : if (!l) { l = prec; x = gtofp(x,l); }
2195 14 : p1 = logabs(x);
2196 14 : if (signe(p1) > 0) { x = ginv(x); fl = !m2; setsigne(p1, -1); }
2197 : /* |x| <= 1 */
2198 14 : y = RIpolylog(m, x, m2, l);
2199 14 : if (m==1)
2200 : {
2201 7 : shiftr_inplace(p1, -1); /* log |x| / 2 */
2202 7 : y = gadd(y, p1);
2203 : }
2204 : else
2205 : { /* m >= 2, \sum_{0 <= k <= m} 2^k B_k/k! (log |x|)^k Li_{m-k}(x),
2206 : with Li_0(x) := -1/2 */
2207 7 : GEN u, t = RIpolylog(m-1, x, m2, l);
2208 7 : u = gneg_i(p1); /* u = 2 B1 log |x| */
2209 7 : y = gadd(y, gmul(u, t));
2210 7 : if (m > 2)
2211 : {
2212 : GEN p2;
2213 7 : shiftr_inplace(p1, 1); /* 2log|x| <= 0 */
2214 7 : constbern(m>>1);
2215 7 : p1 = sqrr(p1);
2216 7 : p2 = shiftr(p1,-1);
2217 21 : for (k = 2; k < m; k += 2)
2218 : {
2219 14 : if (k > 2) p2 = gdivgunextu(gmul(p2,p1),k-1); /* 2^k/k! log^k |x|*/
2220 14 : t = RIpolylog(m-k, x, m2, l);
2221 14 : u = gmul(p2, bernfrac(k));
2222 14 : y = gadd(y, gmul(u, t));
2223 : }
2224 : }
2225 : }
2226 14 : if (fl) y = gneg(y);
2227 14 : return gc_upto(av, y);
2228 : }
2229 :
2230 : static GEN
2231 175 : gpolylog_i(void *E, GEN x, long prec)
2232 : {
2233 175 : pari_sp av = avma;
2234 175 : long i, n, v, m = (long)E;
2235 : GEN a, y;
2236 :
2237 175 : if (m <= 0)
2238 : {
2239 28 : a = gmul(x, poleval(eulerianpol(-m, 0), x));
2240 28 : return gc_upto(av, gdiv(a, gpowgs(gsubsg(1, x), 1-m)));
2241 : }
2242 147 : switch(typ(x))
2243 : {
2244 84 : case t_REAL: case t_COMPLEX: return polylog(m,x,prec);
2245 7 : case t_INTMOD: case t_PADIC: pari_err_IMPL( "padic polylogarithm");
2246 56 : default:
2247 56 : av = avma; if (!(y = toser_i(x))) break;
2248 21 : if (!m) { set_avma(av); return mkfrac(gen_m1,gen_2); }
2249 21 : if (m==1) return gc_upto(av, Li1(y, prec));
2250 21 : if (gequal0(y)) return gc_GEN(av, y);
2251 21 : v = valser(y);
2252 21 : if (v < 0) pari_err_DOMAIN("polylog","valuation", "<", gen_0, x);
2253 14 : if (v > 0) {
2254 7 : n = (lg(y)-3 + v) / v;
2255 7 : a = zeroser(varn(y), lg(y)-2);
2256 35 : for (i=n; i>=1; i--)
2257 28 : a = gmul(y, gadd(a, powis(utoipos(i),-m)));
2258 : } else { /* v == 0 */
2259 7 : long vy = varn(y);
2260 7 : GEN a0 = polcoef_i(y, 0, -1), t = gdiv(derivser(y), y);
2261 7 : a = Li1(y, prec);
2262 14 : for (i=2; i<=m; i++)
2263 7 : a = gadd(gpolylog(i, a0, prec), integ(gmul(t, a), vy));
2264 : }
2265 14 : return gc_upto(av, a);
2266 : }
2267 35 : return trans_evalgen("polylog", E, gpolylog_i, x, prec);
2268 : }
2269 : GEN
2270 133 : gpolylog(long m, GEN x, long prec) { return gpolylog_i((void*)m, x, prec); }
2271 :
2272 : GEN
2273 147 : polylog0(long m, GEN x, long flag, long prec)
2274 : {
2275 147 : switch(flag)
2276 : {
2277 105 : case 0: return gpolylog(m,x,prec);
2278 14 : case 1: return polylogD(m,x,0,prec);
2279 7 : case 2: return polylogD(m,x,1,prec);
2280 14 : case 3: return polylogP(m,x,prec);
2281 7 : default: pari_err_FLAG("polylog");
2282 : }
2283 : return NULL; /* LCOV_EXCL_LINE */
2284 : }
|
