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 488848 : _abs(GEN x)
36 488848 : { return gabs(gtofp(x,LOWDEFAULTPREC), LOWDEFAULTPREC); }
37 : /* can we use asymptotic expansion ? */
38 : static int
39 405158 : bessel_asymp(GEN z, long bit)
40 405158 : { return gcmpgs(_abs(z), (bit+10)/2) >= 0; }
41 :
42 : /* Region I: 0 < Arg z <= Pi, II: -Pi < Arg z <= 0 */
43 : static int
44 532 : regI(GEN z)
45 : {
46 532 : long s = gsigne(imag_i(z));
47 532 : return (s > 0 || (s == 0 && gsigne(real_i(z)) < 0)) ? 1 : 2;
48 : }
49 : /* Region 1: Re(z) >= 0, 2: Re(z) < 0, Im(z) >= 0, 3: Re(z) < 0, Im(z) < 0 */
50 : static int
51 82696 : regJ(GEN z)
52 : {
53 82696 : if (gsigne(real_i(z)) >= 0) return 1;
54 336 : return gsigne(imag_i(z)) >= 0 ? 2 : 3;
55 : }
56 :
57 : /* Hankel's expansions:
58 : * a_k(n) = \prod_{0 <= j < k} (4n^2 - (2j+1)^2)
59 : * C(k)[n,z] = a_k(n) / (k! (8 z)^k)
60 : * A(z) = exp(-z) sum_{k >= 0} C(k)
61 : * A(-z) = exp(z) sum_{k >= 0} (-1)^k C(k)
62 : * J_n(z) ~ [1] (A(z/i) / r + A(-z/i) r) / sqrt(2Pi z)
63 : * [2] (A(z/i) r^3 + A(-z/i) r) / sqrt(2Pi z)
64 : * [3] (A(z/i) / r + A(-z/i) / r^3) / sqrt(2Pi z)
65 : * Y_n(z) ~ [1] i(A(z/i) / r + A(-z/i) r) / sqrt(2Pi z)
66 : * [2] i(A(z/i) (r^3-2/r) + A(-z/i) r) / sqrt(2Pi z)
67 : * [3] i(-A(z/i)/r + A(-z/i)(2r-1/r^3)) / sqrt(2Pi z)
68 : * K_n(z) ~ A(z) Pi / sqrt(2 Pi z)
69 : * I_n(z) ~ [I] (A(-z) + r^2 A(z)) / sqrt(2 Pi z)
70 : * [II](A(-z) + r^(-2) A(z)) / sqrt(2 Pi z) */
71 :
72 : /* set [A(z), A(-z), exp((2*nu+1)*I*Pi/4)] */
73 : static void
74 83690 : hankel_ABr(GEN *pA, GEN *pB, GEN *pr, GEN n, GEN z, long bit)
75 : {
76 83690 : GEN E, P, C, Q = gen_0, zi = ginv(gmul2n(z, 3));
77 83690 : GEN K = gaddgs(_abs(n), 1), n2 = gmul2n(gsqr(n),2);
78 83690 : long prec = nbits2prec(bit), B = bit + 4, m;
79 :
80 83690 : P = C = real_1_bit(bit);
81 83690 : for (m = 1;; m += 2)
82 : {
83 5585602 : C = gmul(C, gdivgu(gmul(gsub(n2, sqru(2*m - 1)), zi), m));
84 5585602 : Q = gadd(Q, C);
85 5585602 : C = gmul(C, gdivgu(gmul(gsub(n2, sqru(2*m + 1)), zi), m + 1));
86 5585602 : P = gadd(P, C);
87 5585602 : if (gexpo(C) < -B && gcmpgs(K, m) <= 0) break;
88 : }
89 83690 : E = gexp(z, prec);
90 83690 : *pA = gdiv(gadd(P, Q), E);
91 83690 : *pB = gmul(gsub(P, Q), E);
92 83690 : *pr = gexp(mulcxI(gmul(gaddgs(gmul2n(n,1), 1), Pi2n(-2, prec))), prec);
93 83690 : }
94 :
95 : /* sqrt(2*Pi*z) */
96 : static GEN
97 83690 : sqz(GEN z, long bit)
98 : {
99 83690 : long prec = nbits2prec(bit);
100 83690 : return gsqrt(gmul(Pi2n(1, prec), z), prec);
101 : }
102 :
103 : static GEN
104 462 : besskasymp(GEN nu, GEN z, long bit)
105 : {
106 : GEN A, B, r;
107 462 : long prec = nbits2prec(bit);
108 462 : hankel_ABr(&A,&B,&r, nu, z, bit);
109 462 : return gdiv(gmul(A, mppi(prec)), sqz(z, bit));
110 : }
111 :
112 : static GEN
113 532 : bessiasymp(GEN nu, GEN z, long bit)
114 : {
115 : GEN A, B, r, R, r2;
116 532 : hankel_ABr(&A,&B,&r, nu, z, bit);
117 532 : r2 = gsqr(r);
118 532 : R = regI(z) == 1 ? gmul(A, r2) : gdiv(A, r2);
119 532 : return gdiv(gadd(B, R), sqz(z, bit));
120 : }
121 :
122 : static GEN
123 82226 : bessjasymp(GEN nu, GEN z, long bit)
124 : {
125 : GEN A, B, r, R;
126 82226 : long reg = regJ(z);
127 82226 : hankel_ABr(&A,&B,&r, nu, mulcxmI(z), bit);
128 82226 : if (reg == 1) R = gadd(gdiv(A, r), gmul(B, r));
129 168 : else if (reg == 2) R = gadd(gmul(A, gpowgs(r, 3)), gmul(B, r));
130 56 : else R = gadd(gdiv(A, r), gdiv(B, gpowgs(r, 3)));
131 82226 : return gdiv(R, sqz(z, bit));
132 : }
133 :
134 : static GEN
135 470 : bessyasymp(GEN nu, GEN z, long bit)
136 : {
137 : GEN A, B, r, R;
138 470 : long reg = regJ(z);
139 470 : hankel_ABr(&A,&B,&r, nu, mulcxmI(z), bit);
140 470 : if (reg == 1) R = gsub(gmul(B, r), gdiv(A, r));
141 168 : else if (reg == 2)
142 112 : R = gadd(gmul(A, gsub(gpowgs(r, 3), gmul2n(ginv(r), 1))), gmul(B, r));
143 : else
144 56 : R = gsub(gmul(B, gsub(gmul2n(r, 1), ginv(gpowgs(r, 3)))), gdiv(A, r));
145 470 : return gdiv(mulcxI(R), sqz(z, bit));
146 : }
147 :
148 : /* n! sum_{0 <= k <= m} x^k / (k!*(k+n)!) */
149 : static GEN
150 318172 : _jbessel(GEN n, GEN x, long m)
151 : {
152 318172 : pari_sp av = avma;
153 318172 : GEN s = gen_1;
154 : long k;
155 :
156 118000260 : for (k = m; k >= 1; k--)
157 : {
158 117682088 : s = gaddsg(1, gdiv(gmul(x,s), gmulgu(gaddgs(n, k), k)));
159 117682088 : if (gc_needed(av,1))
160 : {
161 0 : if (DEBUGMEM>1) pari_warn(warnmem,"besselj");
162 0 : s = gerepileupto(av, s);
163 : }
164 : }
165 318172 : return s;
166 : }
167 :
168 : /* max(2, L * approximate solution to x log x = B) */
169 : static long
170 320887 : bessel_get_lim(double B, double L)
171 320887 : { return maxss(2, L * exp(dbllambertW0(B))); }
172 :
173 42 : static GEN vjbesselh(void* E, GEN z, long prec){return jbesselh((GEN)E,z,prec);}
174 126 : static GEN vjbessel(void* E, GEN z, long prec) {return jbessel((GEN)E,z,prec);}
175 42 : static GEN vibessel(void* E, GEN z, long prec) {return ibessel((GEN)E,z,prec);}
176 126 : static GEN vnbessel(void* E, GEN z, long prec) {return nbessel((GEN)E,z,prec);}
177 42 : static GEN vkbessel(void* E, GEN z, long prec) {return kbessel((GEN)E,z,prec);}
178 :
179 : /* if J != 0 BesselJ, else BesselI. */
180 : static GEN
181 401245 : jbesselintern(GEN n, GEN z, long J, long prec)
182 : {
183 401245 : const char *f = J? "besselj": "besseli";
184 : long i, ki;
185 401245 : pari_sp av = avma;
186 : GEN y;
187 :
188 401245 : switch(typ(z))
189 : {
190 400839 : case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
191 : {
192 400839 : int flz0 = gequal0(z);
193 : long lim, k, precnew, bit;
194 : GEN p1, p2;
195 : double az, L;
196 :
197 400839 : i = precision(z); if (i) prec = i;
198 400839 : if (flz0 && gequal0(n)) return real_1(prec);
199 400839 : bit = prec2nbits(prec);
200 400839 : if (bessel_asymp(z, bit))
201 : {
202 82758 : GEN R = J? bessjasymp(n, z, bit): bessiasymp(n, z, bit);
203 82758 : if (typ(R) == t_COMPLEX && isexactzero(imag_i(n))
204 82072 : && gsigne(real_i(z)) > 0
205 81918 : && isexactzero(imag_i(z))) R = gcopy(gel(R,1));
206 82758 : return gerepileupto(av, R);
207 : }
208 318081 : p2 = gpow(gmul2n(z,-1),n,prec);
209 318053 : p2 = gdiv(p2, ggamma(gaddgs(n,1),prec));
210 318053 : if (flz0) return gerepileupto(av, p2);
211 318053 : az = dblmodulus(z); L = HALF_E * az;
212 318053 : precnew = prec;
213 318053 : if (az >= 1.0) precnew += 1 + nbits2extraprec((long)(az/M_LN2));
214 318053 : if (issmall(n,&ki)) {
215 317136 : k = labs(ki);
216 317136 : n = utoi(k);
217 : } else {
218 833 : i = precision(n);
219 833 : if (i && i < precnew) n = gtofp(n,precnew);
220 : }
221 317969 : z = gtofp(z,precnew);
222 317969 : lim = bessel_get_lim(prec2nbits_mul(prec,M_LN2/2) / L, L);
223 317969 : z = gmul2n(gsqr(z),-2); if (J) z = gneg(z);
224 317969 : p1 = gprec_wtrunc(_jbessel(n,z,lim), prec);
225 317969 : return gerepileupto(av, gmul(p2,p1));
226 : }
227 :
228 14 : case t_PADIC: pari_err_IMPL(stack_strcat("p-adic ",f));
229 392 : default:
230 : {
231 : long v, k, m;
232 392 : if (!(y = toser_i(z))) break;
233 238 : if (issmall(n,&ki)) n = utoi(labs(ki));
234 210 : y = gmul2n(gsqr(y),-2); if (J) y = gneg(y);
235 210 : v = valp(y);
236 210 : if (v < 0) pari_err_DOMAIN(f, "valuation", "<", gen_0, z);
237 203 : if (v == 0) pari_err_IMPL(stack_strcat(f, " around a!=0"));
238 203 : m = lg(y) - 2;
239 203 : k = m - (v >> 1);
240 203 : if (k <= 0) { set_avma(av); return scalarser(gen_1, varn(z), v); }
241 203 : setlg(y, k+2); return gerepileupto(av, _jbessel(n, y, m));
242 : }
243 : }
244 154 : return trans_evalgen(f, (void*)n, J? vjbessel: vibessel, z, prec);
245 : }
246 : GEN
247 396885 : jbessel(GEN n, GEN z, long prec) { return jbesselintern(n,z,1,prec); }
248 : GEN
249 889 : ibessel(GEN n, GEN z, long prec) { return jbesselintern(n,z,0,prec); }
250 :
251 : /* k > 0 */
252 : static GEN
253 119 : _jbesselh(long k, GEN z, long prec)
254 : {
255 119 : GEN s, c, p0, p1, zinv = ginv(z);
256 : long i;
257 :
258 119 : gsincos(z,&s,&c,prec);
259 119 : p1 = gmul(zinv,s);
260 119 : p0 = p1; p1 = gmul(zinv,gsub(p0,c));
261 1134 : for (i = 2; i <= k; i++)
262 : {
263 1015 : GEN p2 = gsub(gmul(gmulsg(2*i-1,zinv), p1), p0);
264 1015 : p0 = p1; p1 = p2;
265 : }
266 119 : return p1;
267 : }
268 :
269 : /* J_{n+1/2}(z) */
270 : GEN
271 315 : jbesselh(GEN n, GEN z, long prec)
272 : {
273 : long k, i;
274 : pari_sp av;
275 : GEN y;
276 :
277 315 : if (typ(n)!=t_INT) pari_err_TYPE("jbesselh",n);
278 203 : k = itos(n);
279 203 : if (k < 0) return jbessel(gadd(ghalf,n), z, prec);
280 :
281 203 : switch(typ(z))
282 : {
283 133 : case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
284 : {
285 : long pr;
286 : GEN p1;
287 133 : if (gequal0(z))
288 : {
289 7 : av = avma;
290 7 : p1 = gmul(gsqrt(gdiv(z,mppi(prec)),prec),gpowgs(z,k));
291 7 : p1 = gdiv(p1, mulu_interval(k+1, 2*k+1)); /* x k! / (2k+1)! */
292 7 : return gerepileupto(av, gmul2n(p1,2*k));
293 : }
294 126 : if ( (pr = precision(z)) ) prec = pr;
295 126 : if (bessel_asymp(z, prec2nbits(prec)))
296 7 : return jbessel(gadd(ghalf,n), z, prec);
297 119 : y = cgetc(prec); av = avma;
298 119 : p1 = gsqrt(gdiv(z, Pi2n(-1,prec)), prec);
299 119 : if (!k)
300 21 : p1 = gmul(p1, gsinc(z, prec));
301 : else
302 : {
303 98 : long bits = BITS_IN_LONG + 2*k * (log2(k) - dbllog2(z));
304 98 : if (bits > 0)
305 : {
306 98 : prec += nbits2extraprec(bits);
307 98 : if (pr) z = gtofp(z, prec);
308 : }
309 98 : p1 = gmul(p1, _jbesselh(k,z,prec));
310 : }
311 119 : set_avma(av); return affc_fixlg(p1, y);
312 : }
313 0 : case t_PADIC: pari_err_IMPL("p-adic jbesselh function");
314 70 : default:
315 : {
316 : long t, v;
317 70 : av = avma; if (!(y = toser_i(z))) break;
318 35 : if (gequal0(y)) return gerepileupto(av, gpowgs(y,k));
319 35 : v = valp(y);
320 35 : if (v < 0) pari_err_DOMAIN("besseljh","valuation", "<", gen_0, z);
321 28 : t = lg(y)-2;
322 28 : if (v) y = sertoser(y, t + (2*k+1)*v);
323 28 : if (!k)
324 7 : y = gsinc(y,prec);
325 : else
326 : {
327 21 : GEN T, a = _jbesselh(k, y, prec);
328 21 : if (v) y = sertoser(y, t + k*v); /* lower precision */
329 21 : y = gdiv(a, gpowgs(y, k));
330 21 : T = cgetg(k+1, t_VECSMALL);
331 168 : for (i = 1; i <= k; i++) T[i] = 2*i+1;
332 21 : y = gmul(y, zv_prod_Z(T));
333 : }
334 28 : return gerepileupto(av, y);
335 : }
336 : }
337 35 : return trans_evalgen("besseljh",(void*)n, vjbesselh, z, prec);
338 : }
339 :
340 : static GEN
341 0 : kbessel2(GEN nu, GEN x, long prec)
342 : {
343 0 : pari_sp av = avma;
344 0 : GEN p1, a, x2 = gshift(x,1);
345 :
346 0 : a = gtofp(gaddgs(gshift(nu,1), 1), prec);
347 0 : p1 = hyperu(gshift(a,-1), a, x2, prec);
348 0 : p1 = gmul(gmul(p1, gpow(x2,nu,prec)), sqrtr(mppi(prec)));
349 0 : return gerepileupto(av, gmul(p1, gexp(gneg(x),prec)));
350 : }
351 :
352 : /* special case of hyperu */
353 : static GEN
354 14 : kbessel1(GEN nu, GEN gx, long prec)
355 : {
356 : GEN x, y, zf, r, u, pi, nu2;
357 14 : long bit, k, k2, n2, n, l = (typ(gx)==t_REAL)? realprec(gx): prec;
358 : pari_sp av;
359 :
360 14 : if (typ(nu)==t_COMPLEX) return kbessel2(nu, gx, l);
361 14 : y = cgetr(l); av = avma;
362 14 : x = gtofp(gx, l);
363 14 : nu = gtofp(nu,l); nu2 = sqrr(nu);
364 14 : shiftr_inplace(nu2,2); togglesign(nu2); /* nu2 = -4nu^2 */
365 14 : n = (long) (prec2nbits_mul(l,M_LN2) + M_PI*fabs(rtodbl(nu))) / 2;
366 14 : bit = prec2nbits(l) - 1;
367 14 : l += EXTRAPRECWORD;
368 14 : pi = mppi(l); n2 = n<<1; r = gmul2n(x,1);
369 14 : if (cmprs(x, n) < 0)
370 : {
371 14 : pari_sp av2 = avma;
372 14 : GEN q, v, c, s = real_1(l), t = real_0(l);
373 1246 : for (k = n2, k2 = 2*n2-1; k > 0; k--, k2 -= 2)
374 : {
375 1232 : GEN ak = divri(addri(nu2, sqru(k2)), mulss(n2<<2, -k));
376 1232 : s = addsr(1, mulrr(ak,s));
377 1232 : t = addsr(k2,mulrr(ak,t));
378 1232 : if (gc_needed(av2,3)) gerepileall(av2, 2, &s,&t);
379 : }
380 14 : shiftr_inplace(t, -1);
381 14 : q = utor(n2, l);
382 14 : zf = sqrtr(divru(pi,n2));
383 14 : u = gprec_wensure(mulrr(zf, s), l);
384 14 : v = gprec_wensure(divrs(addrr(mulrr(t,zf),mulrr(u,nu)),-n2), l);
385 : for(;;)
386 301 : {
387 315 : GEN p1, e, f, d = real_1(l);
388 : pari_sp av3;
389 315 : c = divur(5,q); if (expo(c) >= -1) c = real2n(-1,l);
390 315 : p1 = subsr(1, divrr(r,q)); if (cmprr(c,p1)>0) c = p1;
391 315 : togglesign(c); av3 = avma;
392 315 : e = u;
393 315 : f = v;
394 315 : for (k = 1;; k++)
395 34475 : {
396 34790 : GEN w = addrr(gmul2n(mulur(2*k-1,u), -1), mulrr(subrs(q,k),v));
397 34790 : w = addrr(w, mulrr(nu, subrr(u,gmul2n(v,1))));
398 34790 : u = divru(mulrr(q,v), k);
399 34790 : v = divru(w,k);
400 34790 : d = mulrr(d,c);
401 34790 : e = addrr(e, mulrr(d,u));
402 34790 : f = addrr(f, p1 = mulrr(d,v));
403 34790 : if (expo(p1) - expo(f) <= 1-prec2nbits(realprec(p1))) break;
404 34475 : if (gc_needed(av3,3)) gerepileall(av3,5,&u,&v,&d,&e,&f);
405 : }
406 315 : u = e;
407 315 : v = f;
408 315 : q = mulrr(q, addrs(c,1));
409 315 : if (expo(r) - expo(subrr(q,r)) >= bit) break;
410 301 : gerepileall(av2, 3, &u,&v,&q);
411 : }
412 14 : u = mulrr(u, gpow(divru(x,n),nu,prec));
413 : }
414 : else
415 : {
416 0 : GEN s, zz = ginv(gmul2n(r,2));
417 0 : pari_sp av2 = avma;
418 0 : s = real_1(l);
419 0 : for (k = n2, k2 = 2*n2-1; k > 0; k--, k2 -= 2)
420 : {
421 0 : GEN ak = divru(mulrr(addri(nu2, sqru(k2)), zz), k);
422 0 : s = subsr(1, mulrr(ak,s));
423 0 : if (gc_needed(av2,3)) s = gerepileuptoleaf(av2, s);
424 : }
425 0 : zf = sqrtr(divrr(pi,r));
426 0 : u = mulrr(s, zf);
427 : }
428 14 : affrr(mulrr(u, mpexp(negr(x))), y);
429 14 : set_avma(av); return y;
430 : }
431 :
432 : /* sum_{k=0}^m x^k (H(k)+H(k+n)) / (k! (k+n)!)
433 : * + sum_{k=0}^{n-1} (-x)^(k-n) (n-k-1)!/k! */
434 : static GEN
435 3002 : _kbessel(long n, GEN x, long m, long prec)
436 : {
437 : GEN p1, p2, s, H;
438 3002 : long k, M = m + n, exact = (M <= bit_accuracy(prec));
439 : pari_sp av;
440 :
441 3002 : H = cgetg(M+2,t_VEC); gel(H,1) = gen_0;
442 3002 : if (exact)
443 : {
444 3002 : gel(H,2) = s = gen_1;
445 107856 : for (k=2; k<=M; k++) gel(H,k+1) = s = gdivgu(gaddsg(1,gmulsg(k,s)),k);
446 : }
447 : else
448 : {
449 0 : gel(H,2) = s = real_1(prec);
450 0 : for (k=2; k<=M; k++) gel(H,k+1) = s = divru(addsr(1,mulur(k,s)),k);
451 : }
452 3002 : s = gadd(gel(H,m+1), gel(H,M+1)); av = avma;
453 109378 : for (k = m; k > 0; k--)
454 : {
455 106376 : s = gadd(gadd(gel(H,k),gel(H,k+n)), gdiv(gmul(x,s),mulss(k,k+n)));
456 106376 : if (gc_needed(av,1))
457 : {
458 0 : if (DEBUGMEM>1) pari_warn(warnmem,"_kbessel");
459 0 : s = gerepileupto(av, s);
460 : }
461 : }
462 3002 : p1 = exact? mpfact(n): mpfactr(n,prec);
463 3002 : s = gdiv(s,p1);
464 3002 : if (n)
465 : {
466 472 : x = gneg(ginv(x));
467 472 : p2 = gmulsg(n, gdiv(x,p1));
468 472 : s = gadd(s,p2);
469 1480 : for (k=n-1; k>0; k--)
470 : {
471 1008 : p2 = gmul(p2, gmul(mulss(k,n-k),x));
472 1008 : s = gadd(s,p2);
473 : }
474 : }
475 3002 : return s;
476 : }
477 :
478 : /* N = 1: Bessel N, else Bessel K */
479 : static GEN
480 4578 : kbesselintern(GEN n, GEN z, long N, long prec)
481 : {
482 4578 : const char *f = N? "besseln": "besselk";
483 : long i, k, ki, lim, precnew, fl2, ex, bit;
484 4578 : pari_sp av = avma;
485 : GEN p1, p2, y, p3, pp, pm, s, c;
486 : double az;
487 :
488 4578 : switch(typ(z))
489 : {
490 4193 : case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
491 4193 : if (gequal0(z)) pari_err_DOMAIN(f, "argument", "=", gen_0, z);
492 4193 : i = precision(z); if (i) prec = i;
493 4193 : i = precision(n); if (i && prec > i) prec = i;
494 4193 : bit = prec2nbits(prec);
495 4193 : if (bessel_asymp(z, bit))
496 : {
497 932 : GEN R = N? bessyasymp(n, z, bit): besskasymp(n, z, bit);
498 932 : if (typ(R) == t_COMPLEX && isexactzero(imag_i(n))
499 225 : && gsigne(real_i(z)) > 0
500 85 : && isexactzero(imag_i(z))) R = gcopy(gel(R,1));
501 932 : return gerepileupto(av, R);
502 : }
503 : /* heuristic threshold */
504 3261 : if (!N && !gequal0(n) && gexpo(z) > bit/16 + gexpo(n))
505 14 : return kbessel1(n,z,prec);
506 3240 : az = dblmodulus(z); precnew = prec;
507 3240 : if (az >= 1) precnew += 1 + nbits2extraprec((long)((N?az:2*az)/M_LN2));
508 3240 : z = gtofp(z, precnew);
509 3240 : if (issmall(n,&ki))
510 : {
511 2918 : GEN z2 = gmul2n(z, -1), Z;
512 2918 : double B, L = HALF_E * az;
513 2918 : k = labs(ki);
514 2918 : B = prec2nbits_mul(prec,M_LN2/2) / L;
515 2918 : if (!N) B += 0.367879; /* exp(-1) */
516 2918 : lim = bessel_get_lim(B, L);
517 2918 : Z = gsqr(z2); if (N) Z = gneg(Z);
518 2918 : p1 = gmul(gpowgs(z2,k), _kbessel(k, Z, lim, precnew));
519 2918 : p2 = gadd(mpeuler(precnew), glog(z2,precnew));
520 2918 : p3 = jbesselintern(stoi(k),z,N,precnew);
521 2918 : p2 = gsub(gmul2n(p1,-1),gmul(p2,p3));
522 2918 : p2 = gprec_wtrunc(p2, prec);
523 2918 : if (N)
524 : {
525 510 : p2 = gdiv(p2, Pi2n(-1,prec));
526 510 : if (ki >= 0 || !odd(k)) p2 = gneg(p2);
527 : } else
528 2408 : if (odd(k)) p2 = gneg(p2);
529 2918 : return gerepilecopy(av, p2);
530 : }
531 :
532 259 : n = gtofp(n, precnew);
533 259 : gsincos(gmul(n,mppi(precnew)), &s,&c,precnew);
534 259 : ex = gexpo(s);
535 259 : if (ex < 0) precnew += nbits2extraprec(N? -ex: -2*ex);
536 259 : if (i && i < precnew) {
537 84 : n = gtofp(n,precnew);
538 84 : z = gtofp(z,precnew);
539 84 : gsincos(gmul(n,mppi(precnew)), &s,&c,precnew);
540 : }
541 :
542 259 : pp = jbesselintern(n, z,N,precnew);
543 259 : pm = jbesselintern(gneg(n),z,N,precnew);
544 259 : if (N)
545 189 : p1 = gsub(gmul(c,pp),pm);
546 : else
547 70 : p1 = gmul(gsub(pm,pp), Pi2n(-1,precnew));
548 259 : p1 = gdiv(p1, s);
549 259 : return gerepilecopy(av, gprec_wtrunc(p1,prec));
550 :
551 14 : case t_PADIC: pari_err_IMPL(stack_strcat("p-adic ",f));
552 371 : default:
553 371 : if (!(y = toser_i(z))) break;
554 217 : if (issmall(n,&ki))
555 : {
556 105 : long v, mv, k = labs(ki), m = lg(y)-2;
557 105 : y = gmul2n(gsqr(y),-2); if (N) y = gneg(y);
558 105 : v = valp(y);
559 105 : if (v < 0) pari_err_DOMAIN(f, "valuation", "<", gen_0, z);
560 91 : if (v == 0) pari_err_IMPL(stack_strcat(f, " around a!=0"));
561 91 : mv = m - (v >> 1);
562 91 : if (mv <= 0) { set_avma(av); return scalarser(gen_1, varn(z), v); }
563 84 : setlg(y, mv+2); return gerepilecopy(av, _kbessel(k, y, m, prec));
564 : }
565 98 : if (!issmall(gmul2n(n,1),&ki))
566 70 : pari_err_DOMAIN(f, "2n mod Z", "!=", gen_0, n);
567 28 : k = labs(ki); n = gmul2n(stoi(k),-1);
568 28 : fl2 = (k&3)==1;
569 28 : pm = jbesselintern(gneg(n), y, N, prec);
570 28 : if (N) p1 = pm;
571 : else
572 : {
573 7 : pp = jbesselintern(n, y, N, prec);
574 7 : p2 = gpowgs(y,-k); if (fl2 == 0) p2 = gneg(p2);
575 7 : p3 = gmul2n(diviiexact(mpfact(k + 1),mpfact((k + 1) >> 1)),-(k + 1));
576 7 : p3 = gdivgu(gmul2n(gsqr(p3),1),k);
577 7 : p2 = gmul(p2,p3);
578 7 : p1 = gsub(pp,gmul(p2,pm));
579 : }
580 28 : return gerepileupto(av, fl2? gneg(p1): gcopy(p1));
581 : }
582 154 : return trans_evalgen(f, (void*)n, N? vnbessel: vkbessel, z, prec);
583 : }
584 :
585 : GEN
586 3108 : kbessel(GEN n, GEN z, long prec) { return kbesselintern(n,z,0,prec); }
587 : GEN
588 1470 : ybessel(GEN n, GEN z, long prec) { return kbesselintern(n,z,1,prec); }
589 : /* J + iN */
590 : GEN
591 224 : hbessel1(GEN n, GEN z, long prec)
592 : {
593 224 : pari_sp av = avma;
594 224 : GEN J = jbessel(n,z,prec);
595 196 : GEN Y = ybessel(n,z,prec);
596 182 : return gerepileupto(av, gadd(J, mulcxI(Y)));
597 : }
598 : /* J - iN */
599 : GEN
600 224 : hbessel2(GEN n, GEN z, long prec)
601 : {
602 224 : pari_sp av = avma;
603 224 : GEN J = jbessel(n,z,prec);
604 196 : GEN Y = ybessel(n,z,prec);
605 182 : return gerepileupto(av, gadd(J, mulcxmI(Y)));
606 : }
607 :
608 : static GEN
609 28 : besselrefine(GEN z, GEN nu, GEN (*B)(GEN,GEN,long), long bit)
610 : {
611 28 : GEN z0 = gprec_w(z, DEFAULTPREC), nu1 = gaddgs(nu, 1), t;
612 28 : long e, n, c, j, prec = DEFAULTPREC;
613 :
614 28 : t = gdiv(B(nu1, z0, prec), B(nu, z0, prec));
615 28 : t = gadd(z0, gdiv(gsub(gsqr(z0), gsqr(nu)), gsub(gdiv(nu, z0), t)));
616 28 : e = gexpo(t) - 2 * gexpo(z0) - 1; if (e < 0) e = 0;
617 28 : n = expu(bit + 32 - e);
618 28 : c = 1 + e + ((bit - e) >> n);
619 224 : for (j = 1; j <= n; j++)
620 : {
621 196 : c = 2 * c - e;
622 196 : prec = nbits2prec(c); z = gprec_w(z, prec);
623 196 : t = gdiv(B(nu1, z, prec), B(nu, z, prec));
624 196 : z = gsub(z, ginv(gsub(gdiv(nu, z), t)));
625 : }
626 28 : return gprec_w(z, nbits2prec(bit));
627 : }
628 :
629 : static GEN
630 42 : besselzero(GEN nu, long n, GEN (*B)(GEN,GEN,long), long bit)
631 : {
632 42 : pari_sp av = avma;
633 42 : long prec = nbits2prec(bit), a;
634 : GEN z;
635 42 : if (n <= 0) pari_err_DOMAIN("besselzero", "n", "<=", gen_0, stoi(n));
636 28 : if (n > LONG_MAX / 4) pari_err_OVERFLOW("besselzero");
637 28 : a = 4 * n - (B == jbessel? 1: 3);
638 28 : z = gmul(mppi(prec), gmul2n(gaddgs(gmul2n(nu,1), a), -2));
639 28 : return gerepilecopy(av, besselrefine(z, nu, B, bit));
640 : }
641 : GEN
642 21 : besseljzero(GEN nu, long k, long b) { return besselzero(nu, k, jbessel, b); }
643 : GEN
644 21 : besselyzero(GEN nu, long k, long b) { return besselzero(nu, k, ybessel, b); }
645 :
646 : /***********************************************************************/
647 : /** INCOMPLETE GAMMA FUNCTION **/
648 : /***********************************************************************/
649 : /* mx ~ |x|, b = bit accuracy */
650 : static int
651 16422 : gamma_use_asymp(GEN x, long b)
652 : {
653 : long e;
654 16422 : if (is_real_t(typ(x)))
655 : {
656 12390 : pari_sp av = avma;
657 12390 : return gc_int(av, gcmpgs(R_abs_shallow(x), 3*b / 4) >= 0);
658 : }
659 4032 : e = gexpo(x); return e >= b || dblmodulus(x) >= 3*b / 4;
660 : }
661 : /* x a t_REAL */
662 : static GEN
663 28 : eint1r_asymp(GEN x, GEN expx, long prec)
664 : {
665 28 : pari_sp av = avma, av2;
666 : GEN S, q, z, ix;
667 28 : long oldeq = LONG_MAX, esx = -prec2nbits(prec), j;
668 :
669 28 : if (realprec(x) < prec + EXTRAPREC64) x = rtor(x, prec+EXTRAPREC64);
670 28 : ix = invr(x); q = z = negr(ix);
671 28 : av2 = avma; S = addrs(q, 1);
672 28 : for (j = 2;; j++)
673 1211 : {
674 1239 : long eq = expo(q); if (eq < esx) break;
675 1211 : if ((j & 3) == 0)
676 : { /* guard against divergence */
677 294 : if (eq > oldeq) return gc_NULL(av); /* regressing, abort */
678 294 : oldeq = eq;
679 : }
680 1211 : q = mulrr(q, mulru(z, j)); S = addrr(S, q);
681 1211 : if (gc_needed(av2, 1)) gerepileall(av2, 2, &S, &q);
682 : }
683 28 : if (DEBUGLEVEL > 2) err_printf("eint1: using asymp\n");
684 28 : S = expx? divrr(S, expx): mulrr(S, mpexp(negr(x)));
685 28 : return gerepileuptoleaf(av, mulrr(S, ix));
686 : }
687 : /* cf incgam_asymp(0, x); z = -1/x
688 : * exp(-x)/x * (1 + z + 2! z^2 + ...) */
689 : static GEN
690 105 : eint1_asymp(GEN x, GEN expx, long prec)
691 : {
692 105 : pari_sp av = avma, av2;
693 : GEN S, q, z, ix;
694 105 : long oldeq = LONG_MAX, esx = -prec2nbits(prec), j;
695 :
696 105 : if (typ(x) != t_REAL) x = gtofp(x, prec+EXTRAPREC64);
697 105 : if (typ(x) == t_REAL) return eint1r_asymp(x, expx, prec);
698 105 : ix = ginv(x); q = z = gneg_i(ix);
699 105 : av2 = avma; S = gaddgs(q, 1);
700 105 : for (j = 2;; j++)
701 5824 : {
702 5929 : long eq = gexpo(q); if (eq < esx) break;
703 5824 : if ((j & 3) == 0)
704 : { /* guard against divergence */
705 1442 : if (eq > oldeq) return gc_NULL(av); /* regressing, abort */
706 1442 : oldeq = eq;
707 : }
708 5824 : q = gmul(q, gmulgu(z, j)); S = gadd(S, q);
709 5824 : if (gc_needed(av2, 1)) gerepileall(av2, 2, &S, &q);
710 : }
711 105 : if (DEBUGLEVEL > 2) err_printf("eint1: using asymp\n");
712 105 : S = expx? gdiv(S, expx): gmul(S, gexp(gneg_i(x), prec));
713 105 : return gerepileupto(av, gmul(S, ix));
714 : }
715 :
716 : /* eint1(x) = incgam(0, x); typ(x) = t_REAL, x > 0 */
717 : static GEN
718 6188 : eint1p(GEN x, GEN expx)
719 : {
720 : pari_sp av;
721 6188 : long l = realprec(x), bit = prec2nbits(l), prec, i;
722 : double mx;
723 : GEN z, S, t, H, run;
724 :
725 6188 : if (gamma_use_asymp(x, bit)
726 28 : && (z = eint1r_asymp(x, expx, l))) return z;
727 6160 : mx = rtodbl(x);
728 6160 : prec = l + nbits2extraprec((mx+log(mx))/M_LN2 + 10);
729 6160 : bit = prec2nbits(prec);
730 6160 : run = real_1(prec); x = rtor(x, prec);
731 6160 : av = avma; S = z = t = H = run;
732 623624 : for (i = 2; expo(S) - expo(t) <= bit; i++)
733 : {
734 617464 : H = addrr(H, divru(run,i)); /* H = sum_{k<=i} 1/k */
735 617464 : z = divru(mulrr(x,z), i); /* z = x^(i-1)/i! */
736 617464 : t = mulrr(z, H); S = addrr(S, t);
737 617464 : if ((i & 0x1ff) == 0) gerepileall(av, 4, &z,&t,&S,&H);
738 : }
739 6160 : return subrr(mulrr(x, divrr(S,expx? expx: mpexp(x))),
740 : addrr(mplog(x), mpeuler(prec)));
741 : }
742 : /* eint1(x) = incgam(0, x); typ(x) = t_REAL, x < 0
743 : * rewritten from code contributed by Manfred Radimersky */
744 : static GEN
745 140 : eint1m(GEN x, GEN expx)
746 : {
747 140 : GEN p1, q, S, y, z = cgetg(3, t_COMPLEX);
748 140 : long l = realprec(x), n = prec2nbits(l), j;
749 140 : pari_sp av = avma;
750 :
751 140 : y = rtor(x, l + EXTRAPREC64); setsigne(y,1); /* |x| */
752 140 : if (gamma_use_asymp(y, n))
753 : { /* ~eint1_asymp: asymptotic expansion */
754 14 : p1 = q = invr(y); S = addrs(q, 1);
755 560 : for (j = 2; expo(q) >= -n; j++) {
756 546 : q = mulrr(q, mulru(p1, j));
757 546 : S = addrr(S, q);
758 : }
759 14 : y = mulrr(p1, expx? divrr(S, expx): mulrr(S, mpexp(y)));
760 : }
761 : else
762 : {
763 126 : p1 = q = S = y;
764 24248 : for (j = 2; expo(q) - expo(S) >= -n; j++) {
765 24122 : p1 = mulrr(y, divru(p1, j)); /* (-x)^j/j! */
766 24122 : q = divru(p1, j);
767 24122 : S = addrr(S, q);
768 : }
769 126 : y = addrr(S, addrr(logr_abs(x), mpeuler(l)));
770 : }
771 140 : y = gerepileuptoleaf(av, y); togglesign(y);
772 140 : gel(z, 1) = y;
773 140 : y = mppi(l); setsigne(y, -1);
774 140 : gel(z, 2) = y; return z;
775 : }
776 :
777 : /* real(z*log(z)-z), z = x+iy */
778 : static double
779 7714 : mygamma(double x, double y)
780 : {
781 7714 : if (x == 0.) return -(M_PI/2)*fabs(y);
782 7714 : return (x/2)*log(x*x+y*y)-x-y*atan(y/x);
783 : }
784 :
785 : /* x^s exp(-x) */
786 : static GEN
787 10066 : expmx_xs(GEN s, GEN x, GEN logx, long prec)
788 : {
789 : GEN z;
790 10066 : long ts = typ(s);
791 10066 : if (ts == t_INT || (ts == t_FRAC && absequaliu(gel(s,2), 2)))
792 5257 : z = gmul(gexp(gneg(x), prec), gpow(x, s, prec));
793 : else
794 4809 : z = gexp(gsub(gmul(s, logx? logx: glog(x,prec+EXTRAPREC64)), x), prec);
795 10066 : return z;
796 : }
797 :
798 : /* Not yet: doesn't work at low accuracy
799 : #define INCGAM_CF
800 : */
801 :
802 : #ifdef INCGAM_CF
803 : /* Is s very close to a nonpositive integer ? */
804 : static int
805 : isgammapole(GEN s, long bitprec)
806 : {
807 : pari_sp av = avma;
808 : GEN t = imag_i(s);
809 : long e, b = bitprec - 10;
810 :
811 : if (gexpo(t) > - b) return 0;
812 : s = real_i(s);
813 : if (gsigne(s) > 0 && gexpo(s) > -b) return 0;
814 : (void)grndtoi(s, &e); return gc_bool(av, e < -b);
815 : }
816 :
817 : /* incgam using the continued fraction. x a t_REAL or t_COMPLEX, mx ~ |x|.
818 : * Assume precision(s), precision(x) >= prec */
819 : static GEN
820 : incgam_cf(GEN s, GEN x, double mx, long prec)
821 : {
822 : GEN ms, y, S;
823 : long n, i, j, LS, bitprec = prec2nbits(prec);
824 : double rs, is, m;
825 :
826 : if (typ(s) == t_COMPLEX)
827 : {
828 : rs = gtodouble(gel(s,1));
829 : is = gtodouble(gel(s,2));
830 : }
831 : else
832 : {
833 : rs = gtodouble(s);
834 : is = 0.;
835 : }
836 : if (isgammapole(s, bitprec)) LS = 0;
837 : else
838 : {
839 : double bit, LGS = mygamma(rs,is);
840 : LS = LGS <= 0 ? 0: ceil(LGS);
841 : bit = (LGS - (rs-1)*log(mx) + mx)/M_LN2;
842 : if (bit > 0)
843 : {
844 : prec += nbits2extraprec((long)bit);
845 : x = gtofp(x, prec);
846 : if (isinexactreal(s)) s = gtofp(s, prec);
847 : }
848 : }
849 : /* |ln(2*gamma(s)*sin(s*Pi))| <= ln(2) + |lngamma(s)| + |Im(s)*Pi|*/
850 : m = bitprec*M_LN2 + LS + M_LN2 + fabs(is)*M_PI + mx;
851 : if (rs < 1) m += (1 - rs)*log(mx);
852 : m /= 4;
853 : n = (long)(1 + m*m/mx);
854 : y = expmx_xs(gsubgs(s,1), x, NULL, prec);
855 : if (rs >= 0 && bitprec >= 512)
856 : {
857 : GEN A = cgetg(n+1, t_VEC), B = cgetg(n+1, t_VEC);
858 : ms = gsubsg(1, s);
859 : for (j = 1; j <= n; ++j)
860 : {
861 : gel(A,j) = ms;
862 : gel(B,j) = gmulsg(j, gsubgs(s,j));
863 : ms = gaddgs(ms, 2);
864 : }
865 : S = contfraceval_inv(mkvec2(A,B), x, -1);
866 : }
867 : else
868 : {
869 : GEN x_s = gsub(x, s);
870 : pari_sp av2 = avma;
871 : S = gdiv(gsubgs(s,n), gaddgs(x_s,n<<1));
872 : for (i=n-1; i >= 1; i--)
873 : {
874 : S = gdiv(gsubgs(s,i), gadd(gaddgs(x_s,i<<1),gmulsg(i,S)));
875 : if (gc_needed(av2,3))
876 : {
877 : if(DEBUGMEM>1) pari_warn(warnmem,"incgam_cf");
878 : S = gerepileupto(av2, S);
879 : }
880 : }
881 : S = gaddgs(S,1);
882 : }
883 : return gmul(y, S);
884 : }
885 : #endif
886 :
887 : static double
888 5642 : findextraincgam(GEN s, GEN x)
889 : {
890 5642 : double sig = gtodouble(real_i(s)), t = gtodouble(imag_i(s));
891 5642 : double xr = gtodouble(real_i(x)), xi = gtodouble(imag_i(x));
892 5642 : double exd = 0., Nx = xr*xr + xi*xi, D = Nx - t*t;
893 : long n;
894 :
895 5642 : if (xr < 0)
896 : {
897 833 : long ex = gexpo(x);
898 833 : if (ex > 0 && ex > gexpo(s)) exd = sqrt(Nx)*log(Nx)/2; /* |x| log |x| */
899 : }
900 5642 : if (D <= 0.) return exd;
901 4200 : n = (long)(sqrt(D)-sig);
902 4200 : if (n <= 0) return exd;
903 1512 : return maxdd(exd, (n*log(Nx)/2 - mygamma(sig+n, t) + mygamma(sig, t)) / M_LN2);
904 : }
905 :
906 : /* use exp(-x) * (x^s/s) * sum_{k >= 0} x^k / prod(i=1, k, s+i) */
907 : static GEN
908 5649 : incgamc_i(GEN s, GEN x, long *ptexd, long prec)
909 : {
910 : GEN S, t, y;
911 : long l, n, i, exd;
912 5649 : pari_sp av = avma, av2;
913 :
914 5649 : if (gequal0(x))
915 : {
916 7 : if (ptexd) *ptexd = 0.;
917 7 : return gtofp(x, prec);
918 : }
919 5642 : l = precision(x);
920 5642 : if (!l) l = prec;
921 5642 : n = -prec2nbits(l)-1;
922 5642 : exd = (long)findextraincgam(s, x);
923 5642 : if (ptexd) *ptexd = exd;
924 5642 : if (exd > 0)
925 : {
926 1393 : long p = l + nbits2extraprec(exd);
927 1393 : x = gtofp(x, p);
928 1393 : if (isinexactreal(s)) s = gtofp(s, p);
929 : }
930 4249 : else x = gtofp(x, l+EXTRAPREC64);
931 5642 : av2 = avma;
932 5642 : S = gdiv(x, gaddsg(1,s));
933 5642 : t = gaddsg(1, S);
934 753599 : for (i=2; gexpo(S) >= n; i++)
935 : {
936 747957 : S = gdiv(gmul(x,S), gaddsg(i,s)); /* x^i / ((s+1)...(s+i)) */
937 747957 : t = gadd(S,t);
938 747957 : if (gc_needed(av2,3))
939 : {
940 0 : if(DEBUGMEM>1) pari_warn(warnmem,"incgamc");
941 0 : gerepileall(av2, 2, &S, &t);
942 : }
943 : }
944 5642 : y = expmx_xs(s, x, NULL, prec);
945 5642 : return gerepileupto(av, gmul(gdiv(y,s), t));
946 : }
947 :
948 : GEN
949 1456 : incgamc(GEN s, GEN x, long prec)
950 1456 : { return incgamc_i(s, x, NULL, prec); }
951 :
952 : /* incgamma using asymptotic expansion:
953 : * exp(-x)x^(s-1)(1 + (s-1)/x + (s-1)(s-2)/x^2 + ...) */
954 : static GEN
955 2716 : incgam_asymp(GEN s, GEN x, long prec)
956 : {
957 2716 : pari_sp av = avma, av2;
958 : GEN S, q, cox, invx;
959 2716 : long oldeq = LONG_MAX, eq, esx, j;
960 2716 : int flint = (typ(s) == t_INT && signe(s) > 0);
961 :
962 2716 : x = gtofp(x,prec+EXTRAPREC64);
963 2716 : invx = ginv(x);
964 2716 : esx = -prec2nbits(prec);
965 2716 : av2 = avma;
966 2716 : q = gmul(gsubgs(s,1), invx);
967 2716 : S = gaddgs(q, 1);
968 2716 : for (j = 2;; j++)
969 : {
970 123879 : eq = gexpo(q); if (eq < esx) break;
971 121282 : if (!flint && (j & 3) == 0)
972 : { /* guard against divergence */
973 15778 : if (eq > oldeq) return gc_NULL(av); /* regressing, abort */
974 15659 : oldeq = eq;
975 : }
976 121163 : q = gmul(q, gmul(gsubgs(s,j), invx));
977 121163 : S = gadd(S, q);
978 121163 : if (gc_needed(av2, 1)) gerepileall(av2, 2, &S, &q);
979 : }
980 2597 : if (DEBUGLEVEL > 2) err_printf("incgam: using asymp\n");
981 2597 : cox = expmx_xs(gsubgs(s,1), x, NULL, prec);
982 2597 : return gerepileupto(av, gmul(cox, S));
983 : }
984 :
985 : /* gasx = incgam(s-n,x). Compute incgam(s,x)
986 : * = (s-1)(s-2)...(s-n)gasx + exp(-x)x^(s-1) *
987 : * (1 + (s-1)/x + ... + (s-1)(s-2)...(s-n+1)/x^(n-1)) */
988 : static GEN
989 546 : incgam_asymp_partial(GEN s, GEN x, GEN gasx, long n, long prec)
990 : {
991 : pari_sp av;
992 546 : GEN S, q, cox, invx, s1 = gsubgs(s, 1), sprod;
993 : long j;
994 546 : cox = expmx_xs(s1, x, NULL, prec);
995 546 : if (n == 1) return gadd(cox, gmul(s1, gasx));
996 546 : invx = ginv(x);
997 546 : av = avma;
998 546 : q = gmul(s1, invx);
999 546 : S = gaddgs(q, 1);
1000 52164 : for (j = 2; j < n; j++)
1001 : {
1002 51618 : q = gmul(q, gmul(gsubgs(s, j), invx));
1003 51618 : S = gadd(S, q);
1004 51618 : if (gc_needed(av, 2)) gerepileall(av, 2, &S, &q);
1005 : }
1006 546 : sprod = gmul(gmul(q, gpowgs(x, n-1)), gsubgs(s, n));
1007 546 : return gadd(gmul(cox, S), gmul(sprod, gasx));
1008 : }
1009 :
1010 : /* Assume s != 0; called when Re(s) <= 1/2 */
1011 : static GEN
1012 2401 : incgamspec(GEN s, GEN x, GEN g, long prec)
1013 : {
1014 2401 : GEN q, S, cox = gen_0, P, sk, S1, S2, S3, F3, logx, mx;
1015 2401 : long n, esk, E, k = itos(ground(gneg(real_i(s)))); /* >= 0 */
1016 :
1017 2401 : if (k && gexpo(x) > 0)
1018 : {
1019 245 : GEN xk = gdivgu(x, k);
1020 245 : long bitprec = prec2nbits(prec);
1021 245 : double d = (gexpo(xk) > bitprec)? bitprec*M_LN2: log(dblmodulus(xk));
1022 245 : d = k * (d + 1) / M_LN2;
1023 245 : if (d > 0) prec += nbits2extraprec((long)d);
1024 245 : if (isinexactreal(s)) s = gtofp(s, prec);
1025 : }
1026 2401 : x = gtofp(x, maxss(precision(x), prec) + EXTRAPREC64);
1027 2401 : sk = gaddgs(s, k); /* |Re(sk)| <= 1/2 */
1028 2401 : logx = glog(x, prec);
1029 2401 : mx = gneg(x);
1030 2401 : if (k == 0) { S = gen_0; P = gen_1; }
1031 : else
1032 : {
1033 : long j;
1034 854 : q = ginv(s); S = q; P = s;
1035 16926 : for (j = 1; j < k; j++)
1036 : {
1037 16072 : GEN sj = gaddgs(s, j);
1038 16072 : q = gmul(q, gdiv(x, sj));
1039 16072 : S = gadd(S, q);
1040 16072 : P = gmul(P, sj);
1041 : }
1042 854 : cox = expmx_xs(s, x, logx, prec); /* x^s exp(-x) */
1043 854 : S = gmul(S, gneg(cox));
1044 : }
1045 2401 : if (k && gequal0(sk))
1046 175 : return gadd(S, gdiv(eint1(x, prec), P));
1047 2226 : esk = gexpo(sk);
1048 2226 : if (esk > -7)
1049 : {
1050 1015 : GEN a, b, PG = gmul(sk, P);
1051 1015 : if (g) g = gmul(g, PG);
1052 1015 : a = incgam0(gaddgs(sk,1), x, g, prec);
1053 1015 : if (k == 0) cox = expmx_xs(s, x, logx, prec);
1054 1015 : b = gmul(gpowgs(x, k), cox);
1055 1015 : return gadd(S, gdiv(gsub(a, b), PG));
1056 : }
1057 1211 : E = prec2nbits(prec) + 1;
1058 1211 : if (gexpo(x) > 0)
1059 : {
1060 420 : long X = (long)(dblmodulus(x)/M_LN2);
1061 420 : prec += 2*nbits2extraprec(X);
1062 420 : x = gtofp(x, prec); mx = gneg(x);
1063 420 : logx = glog(x, prec); sk = gtofp(sk, prec);
1064 420 : E += X;
1065 : }
1066 1211 : if (isinexactreal(sk)) sk = gtofp(sk, prec+EXTRAPREC64);
1067 : /* |sk| < 2^-7 is small, guard against cancellation */
1068 1211 : F3 = gexpm1(gmul(sk, logx), prec);
1069 : /* ( gamma(1+sk) - exp(sk log(x))) ) / sk */
1070 1211 : S1 = gdiv(gsub(ggamma1m1(sk, prec+EXTRAPREC64), F3), sk);
1071 1211 : q = x; S3 = gdiv(x, gaddsg(1,sk));
1072 255523 : for (n = 2; gexpo(q) - gexpo(S3) > -E; n++)
1073 : {
1074 254312 : q = gmul(q, gdivgu(mx, n));
1075 254312 : S3 = gadd(S3, gdiv(q, gaddsg(n, sk)));
1076 : }
1077 1211 : S2 = gadd(gadd(S1, S3), gmul(F3, S3));
1078 1211 : return gadd(S, gdiv(S2, P));
1079 : }
1080 :
1081 : /* return |x| */
1082 : double
1083 9419000 : dblmodulus(GEN x)
1084 : {
1085 9419000 : if (typ(x) == t_COMPLEX)
1086 : {
1087 1505049 : double a = gtodouble(gel(x,1));
1088 1505049 : double b = gtodouble(gel(x,2));
1089 1505049 : return sqrt(a*a + b*b);
1090 : }
1091 : else
1092 7913951 : return fabs(gtodouble(x));
1093 : }
1094 :
1095 : /* Driver routine. If g != NULL, assume that g=gamma(s,prec). */
1096 : GEN
1097 11543 : incgam0(GEN s, GEN x, GEN g, long prec)
1098 : {
1099 : pari_sp av;
1100 : long E, l;
1101 : GEN z, rs, is;
1102 :
1103 11543 : if (gequal0(x)) return g? gcopy(g): ggamma(s,prec);
1104 11543 : if (gequal0(s)) return eint1(x, prec);
1105 9737 : l = precision(s); if (!l) l = prec;
1106 9737 : E = prec2nbits(l);
1107 9737 : if (gamma_use_asymp(x, E) ||
1108 8316 : (typ(s) == t_INT && signe(s) > 0 && gexpo(x) >= expi(s)))
1109 2716 : if ((z = incgam_asymp(s, x, l))) return z;
1110 7140 : av = avma; E++;
1111 7140 : rs = real_i(s);
1112 7140 : is = imag_i(s);
1113 : #ifdef INCGAM_CF
1114 : /* Can one use continued fraction ? */
1115 : if (gequal0(is) && gequal0(imag_i(x)) && gsigne(x) > 0)
1116 : {
1117 : double sd = gtodouble(rs), LB, UB;
1118 : double xd = gtodouble(real_i(x));
1119 : if (sd > 0) {
1120 : LB = 15 + 0.1205*E;
1121 : UB = 5 + 0.752*E;
1122 : } else {
1123 : LB = -6 + 0.1205*E;
1124 : UB = 5 + 0.752*E + fabs(sd)/54.;
1125 : }
1126 : if (xd >= LB && xd <= UB)
1127 : {
1128 : if (DEBUGLEVEL > 2) err_printf("incgam: using continued fraction\n");
1129 : return gerepileupto(av, incgam_cf(s, x, xd, prec));
1130 : }
1131 : }
1132 : #endif
1133 7140 : if (gsigne(rs) > 0 && gexpo(rs) >= -1)
1134 : { /* use complementary incomplete gamma */
1135 4739 : long n, egs, exd, precg, es = gexpo(s);
1136 4739 : if (es < 0) {
1137 595 : l += nbits2extraprec(-es) + 1;
1138 595 : x = gtofp(x, l);
1139 595 : if (isinexactreal(s)) s = gtofp(s, l);
1140 : }
1141 4739 : n = itos(gceil(rs));
1142 4739 : if (n > 100)
1143 : {
1144 : GEN gasx;
1145 546 : n -= 100;
1146 546 : if (es > 0)
1147 : {
1148 546 : es = mygamma(gtodouble(rs) - n, gtodouble(is)) / M_LN2;
1149 546 : if (es > 0)
1150 : {
1151 546 : l += nbits2extraprec(es);
1152 546 : x = gtofp(x, l);
1153 546 : if (isinexactreal(s)) s = gtofp(s, l);
1154 : }
1155 : }
1156 546 : gasx = incgam0(gsubgs(s, n), x, NULL, prec);
1157 546 : return gerepileupto(av, incgam_asymp_partial(s, x, gasx, n, prec));
1158 : }
1159 4193 : if (DEBUGLEVEL > 2) err_printf("incgam: using power series 1\n");
1160 : /* egs ~ expo(gamma(s)) */
1161 4193 : precg = g? precision(g): 0;
1162 4193 : egs = g? gexpo(g): (long)(mygamma(gtodouble(rs), gtodouble(is)) / M_LN2);
1163 4193 : if (egs > 0) {
1164 1946 : l += nbits2extraprec(egs) + 1;
1165 1946 : x = gtofp(x, l);
1166 1946 : if (isinexactreal(s)) s = gtofp(s, l);
1167 1946 : if (precg < l) g = NULL;
1168 : }
1169 4193 : z = incgamc_i(s, x, &exd, l);
1170 4193 : if (exd > 0)
1171 : {
1172 896 : l += nbits2extraprec(exd);
1173 896 : if (isinexactreal(s)) s = gtofp(s, l);
1174 896 : if (precg < l) g = NULL;
1175 : }
1176 : else
1177 : { /* gamma(s) negligible ? Compute to lower accuracy */
1178 3297 : long e = gexpo(z) - egs;
1179 3297 : if (e > 3)
1180 : {
1181 420 : E -= e;
1182 420 : if (E <= 0) g = gen_0; else if (!g) g = ggamma(s, nbits2prec(E));
1183 : }
1184 : }
1185 : /* worry about possible cancellation */
1186 4193 : if (!g) g = ggamma(s, maxss(l,precision(z)));
1187 4193 : return gerepileupto(av, gsub(g,z));
1188 : }
1189 2401 : if (DEBUGLEVEL > 2) err_printf("incgam: using power series 2\n");
1190 2401 : return gerepileupto(av, incgamspec(s, x, g, l));
1191 : }
1192 :
1193 : GEN
1194 1106 : incgam(GEN s, GEN x, long prec) { return incgam0(s, x, NULL, prec); }
1195 :
1196 : /* x a t_REAL */
1197 : GEN
1198 2156 : mpeint1(GEN x, GEN expx)
1199 : {
1200 2156 : long s = signe(x);
1201 : pari_sp av;
1202 : GEN z;
1203 2156 : if (!s) pari_err_DOMAIN("eint1", "x","=",gen_0, x);
1204 2149 : if (s < 0) return eint1m(x, expx);
1205 2009 : z = cgetr(realprec(x));
1206 2009 : av = avma; affrr(eint1p(x, expx), z);
1207 2009 : set_avma(av); return z;
1208 : }
1209 :
1210 : static GEN
1211 357 : cxeint1(GEN x, long prec)
1212 : {
1213 357 : pari_sp av = avma, av2;
1214 : GEN q, S, run, z, H;
1215 357 : long n, E = prec2nbits(prec);
1216 :
1217 357 : if (gamma_use_asymp(x, E) && (z = eint1_asymp(x, NULL, prec))) return z;
1218 252 : E++;
1219 252 : if (gexpo(x) > 0)
1220 : { /* take cancellation into account, log2(\sum |x|^n / n!) = |x| / log(2) */
1221 42 : double dbx = dblmodulus(x);
1222 42 : long X = (long)((dbx + log(dbx))/M_LN2 + 10);
1223 42 : prec += nbits2extraprec(X);
1224 42 : x = gtofp(x, prec); E += X;
1225 : }
1226 252 : if (DEBUGLEVEL > 2) err_printf("eint1: using power series\n");
1227 252 : run = real_1(prec);
1228 252 : av2 = avma;
1229 252 : S = z = q = H = run;
1230 48384 : for (n = 2; gexpo(q) - gexpo(S) >= -E; n++)
1231 : {
1232 48132 : H = addrr(H, divru(run, n)); /* H = sum_{k<=n} 1/k */
1233 48132 : z = gdivgu(gmul(x,z), n); /* z = x^(n-1)/n! */
1234 48132 : q = gmul(z, H); S = gadd(S, q);
1235 48132 : if ((n & 0x1ff) == 0) gerepileall(av2, 4, &z, &q, &S, &H);
1236 : }
1237 252 : S = gmul(gmul(x, S), gexp(gneg_i(x), prec));
1238 252 : return gerepileupto(av, gsub(S, gadd(glog(x, prec), mpeuler(prec))));
1239 : }
1240 :
1241 : GEN
1242 2513 : eint1(GEN x, long prec)
1243 : {
1244 2513 : switch(typ(x))
1245 : {
1246 357 : case t_COMPLEX: return cxeint1(x, prec);
1247 1771 : case t_REAL: break;
1248 385 : default: x = gtofp(x, prec);
1249 : }
1250 2156 : return mpeint1(x,NULL);
1251 : }
1252 :
1253 : GEN
1254 49 : veceint1(GEN C, GEN nmax, long prec)
1255 : {
1256 49 : if (!nmax) return eint1(C,prec);
1257 7 : if (typ(nmax) != t_INT) pari_err_TYPE("veceint1",nmax);
1258 7 : if (typ(C) != t_REAL) {
1259 7 : C = gtofp(C, prec);
1260 7 : if (typ(C) != t_REAL) pari_err_TYPE("veceint1",C);
1261 : }
1262 7 : if (signe(C) <= 0) pari_err_DOMAIN("veceint1", "argument", "<=", gen_0,C);
1263 7 : return mpveceint1(C, NULL, itos(nmax));
1264 : }
1265 :
1266 : /* j > 0, a t_REAL. Return sum_{m >= 0} a^m / j(j+1)...(j+m)).
1267 : * Stop when expo(summand) < E; note that s(j-1) = (a s(j) + 1) / (j-1). */
1268 : static GEN
1269 259 : mp_sum_j(GEN a, long j, long E, long prec)
1270 : {
1271 259 : pari_sp av = avma;
1272 259 : GEN q = divru(real_1(prec), j), s = q;
1273 : long m;
1274 259 : for (m = 0;; m++)
1275 : {
1276 4871 : if (expo(q) < E) break;
1277 4612 : q = mulrr(q, divru(a, m+j));
1278 4612 : s = addrr(s, q);
1279 : }
1280 259 : return gerepileuptoleaf(av, s);
1281 : }
1282 : /* Return the s_a(j), j <= J */
1283 : static GEN
1284 259 : sum_jall(GEN a, long J, long prec)
1285 : {
1286 259 : GEN s = cgetg(J+1, t_VEC);
1287 259 : long j, E = -prec2nbits(prec) - 5;
1288 259 : gel(s, J) = mp_sum_j(a, J, E, prec);
1289 11103 : for (j = J-1; j; j--)
1290 10844 : gel(s,j) = divru(addrs(mulrr(a, gel(s,j+1)), 1), j);
1291 259 : return s;
1292 : }
1293 :
1294 : /* T a dense t_POL with t_REAL coeffs. Return T(n) [faster than poleval] */
1295 : static GEN
1296 481831 : rX_s_eval(GEN T, long n)
1297 : {
1298 481831 : long i = lg(T)-1;
1299 481831 : GEN c = gel(T,i);
1300 10743228 : for (i--; i>=2; i--) c = gadd(mulrs(c,n),gel(T,i));
1301 481831 : return c;
1302 : }
1303 :
1304 : /* C>0 t_REAL, eC = exp(C). Return eint1(n*C) for 1<=n<=N. Absolute accuracy */
1305 : GEN
1306 266 : mpveceint1(GEN C, GEN eC, long N)
1307 : {
1308 266 : const long prec = realprec(C);
1309 266 : long Nmin = 15; /* >= 1. E.g. between 10 and 30, but little effect */
1310 266 : GEN en, v, w = cgetg(N+1, t_VEC);
1311 : pari_sp av0;
1312 : double DL;
1313 : long n, j, jmax, jmin;
1314 266 : if (!N) return w;
1315 486017 : for (n = 1; n <= N; n++) gel(w,n) = cgetr(prec);
1316 266 : av0 = avma;
1317 266 : if (N < Nmin) Nmin = N;
1318 266 : if (!eC) eC = mpexp(C);
1319 266 : en = eC; affrr(eint1p(C, en), gel(w,1));
1320 3920 : for (n = 2; n <= Nmin; n++)
1321 : {
1322 : pari_sp av2;
1323 3654 : en = mulrr(en,eC); /* exp(n C) */
1324 3654 : av2 = avma;
1325 3654 : affrr(eint1p(mulru(C,n), en), gel(w,n));
1326 3654 : set_avma(av2);
1327 : }
1328 266 : if (Nmin == N) { set_avma(av0); return w; }
1329 :
1330 259 : DL = prec2nbits_mul(prec, M_LN2) + 5;
1331 259 : jmin = ceil(DL/log((double)N)) + 1;
1332 259 : jmax = ceil(DL/log((double)Nmin)) + 1;
1333 259 : v = sum_jall(C, jmax, prec);
1334 259 : en = powrs(eC, -N); /* exp(-N C) */
1335 259 : affrr(eint1p(mulru(C,N), invr(en)), gel(w,N));
1336 7038 : for (j = jmin, n = N-1; j <= jmax; j++)
1337 : {
1338 6779 : long limN = maxss((long)ceil(exp(DL/j)), Nmin);
1339 : GEN polsh;
1340 6779 : setlg(v, j+1);
1341 6779 : polsh = RgV_to_RgX_reverse(v, 0);
1342 488610 : for (; n >= limN; n--)
1343 : {
1344 481831 : pari_sp av2 = avma;
1345 481831 : GEN S = divri(mulrr(en, rX_s_eval(polsh, -n)), powuu(n,j));
1346 : /* w[n+1] - exp(-n C) * polsh(-n) / (-n)^j */
1347 481831 : GEN c = odd(j)? addrr(gel(w,n+1), S) : subrr(gel(w,n+1), S);
1348 481831 : affrr(c, gel(w,n)); set_avma(av2);
1349 481831 : en = mulrr(en,eC); /* exp(-n C) */
1350 : }
1351 : }
1352 259 : set_avma(av0); return w;
1353 : }
1354 :
1355 : /* erfc via numerical integration : assume real(x)>=1 */
1356 : static GEN
1357 14 : cxerfc_r1(GEN x, long prec)
1358 : {
1359 : GEN h, h2, eh2, denom, res, lambda;
1360 : long u, v;
1361 14 : const double D = prec2nbits_mul(prec, M_LN2);
1362 14 : const long npoints = (long)ceil(D/M_PI)+1;
1363 14 : pari_sp av = avma;
1364 : {
1365 14 : double t = exp(-2*M_PI*M_PI/D); /* ~exp(-2*h^2) */
1366 14 : v = 30; /* bits that fit in both long and double mantissa */
1367 14 : u = (long)floor(t*(1L<<v));
1368 : /* define exp(-2*h^2) to be u*2^(-v) */
1369 : }
1370 14 : incrprec(prec);
1371 14 : x = gtofp(x,prec);
1372 14 : eh2 = sqrtr_abs(rtor(shiftr(dbltor(u),-v),prec));
1373 14 : h2 = negr(logr_abs(eh2));
1374 14 : h = sqrtr_abs(h2);
1375 14 : lambda = gdiv(x,h);
1376 14 : denom = gsqr(lambda);
1377 : { /* res = h/x + 2*x*h*sum(k=1,npoints,exp(-(k*h)^2)/(lambda^2+k^2)); */
1378 : GEN Uk; /* = exp(-(kh)^2) */
1379 14 : GEN Vk = eh2;/* = exp(-(2k+1)h^2) */
1380 14 : pari_sp av2 = avma;
1381 : long k;
1382 : /* k = 0 moved out for efficiency */
1383 14 : denom = gaddsg(1,denom);
1384 14 : Uk = Vk;
1385 14 : Vk = mulur(u,Vk); shiftr_inplace(Vk, -v);
1386 14 : res = gdiv(Uk, denom);
1387 420 : for (k = 1; k < npoints; k++)
1388 : {
1389 406 : if ((k & 255) == 0) gerepileall(av2,4,&denom,&Uk,&Vk,&res);
1390 406 : denom = gaddsg(2*k+1,denom);
1391 406 : Uk = mpmul(Uk,Vk);
1392 406 : Vk = mulur(u,Vk); shiftr_inplace(Vk, -v);
1393 406 : res = gadd(res, gdiv(Uk, denom));
1394 : }
1395 : }
1396 14 : res = gmul(res, gshift(lambda,1));
1397 : /* 0 term : */
1398 14 : res = gadd(res, ginv(lambda));
1399 14 : res = gmul(res, gdiv(gexp(gneg(gsqr(x)), prec), mppi(prec)));
1400 14 : if (rtodbl(real_i(x)) < sqrt(D))
1401 : {
1402 14 : GEN t = gmul(divrr(Pi2n(1,prec),h), x);
1403 14 : res = gsub(res, gdivsg(2, cxexpm1(t, prec)));
1404 : }
1405 14 : return gerepileupto(av,res);
1406 : }
1407 :
1408 : GEN
1409 56 : gerfc(GEN x, long prec)
1410 : {
1411 : GEN z, xr, xi, res;
1412 : long s;
1413 : pari_sp av;
1414 :
1415 56 : x = trans_fix_arg(&prec,&x,&xr,&xi, &av,&res);
1416 56 : s = signe(xr);
1417 56 : if (s > 0 || (s == 0 && signe(xi) >= 0)) {
1418 84 : if (cmprs(xr, 1) > 0) /* use numerical integration */
1419 14 : z = cxerfc_r1(x, prec);
1420 : else
1421 : { /* erfc(x) = incgam(1/2,x^2)/sqrt(Pi) */
1422 28 : GEN sqrtpi = sqrtr(mppi(prec));
1423 28 : z = incgam0(ghalf, gsqr(x), sqrtpi, prec);
1424 28 : z = gdiv(z, sqrtpi);
1425 : }
1426 : }
1427 : else
1428 : { /* erfc(-x)=2-erfc(x) */
1429 : /* FIXME could decrease prec
1430 : long size = nbits2extraprec((imag(x)^2-real(x)^2)/log(2));
1431 : prec = size > 0 ? prec : prec + size;
1432 : */
1433 : /* NOT gsubsg(2, ...) : would create a result of
1434 : * huge accuracy if re(x)>>1, rounded to 2 by subsequent affc_fixlg... */
1435 14 : z = gsub(real2n(1,prec+EXTRAPREC64), gerfc(gneg(x), prec));
1436 : }
1437 56 : set_avma(av); return affc_fixlg(z, res);
1438 : }
1439 :
1440 : /***********************************************************************/
1441 : /** **/
1442 : /** RIEMANN ZETA FUNCTION **/
1443 : /** **/
1444 : /***********************************************************************/
1445 : static const double log2PI = 1.83787706641;
1446 :
1447 : static double
1448 4578 : get_xinf(double beta)
1449 : {
1450 4578 : const double maxbeta = 0.06415003; /* 3^(-2.5) */
1451 : double x0, y0, x1;
1452 :
1453 4578 : if (beta < maxbeta) return beta + pow(3*beta, 1.0/3.0);
1454 4578 : x0 = beta + M_PI/2.0;
1455 : for(;;)
1456 : {
1457 10486 : y0 = x0*x0;
1458 7532 : x1 = (beta+atan(x0)) * (1+y0) / y0 - 1/x0;
1459 7532 : if (0.99*x0 < x1) return x1;
1460 2954 : x0 = x1;
1461 : }
1462 : }
1463 : /* optimize for zeta( s + it, prec ), assume |s-1| > 0.1
1464 : * (if gexpo(u = s-1) < -5, we use the functional equation s->1-s) */
1465 : static void
1466 8575 : optim_zeta(GEN S, long prec, long *pp, long *pn)
1467 : {
1468 : double s, t, alpha, beta, n, B;
1469 : long p;
1470 8575 : if (typ(S) == t_REAL) {
1471 266 : s = rtodbl(S);
1472 266 : t = 0.;
1473 : } else {
1474 8309 : s = rtodbl(gel(S,1));
1475 8309 : t = fabs( rtodbl(gel(S,2)) );
1476 : }
1477 :
1478 8575 : B = prec2nbits_mul(prec, M_LN2);
1479 8575 : if (s > 0 && !t) /* positive real input */
1480 : {
1481 238 : beta = B + 0.61 + s*(log2PI - log(s));
1482 476 : if (beta > 0)
1483 : {
1484 231 : p = (long)ceil(beta / 2.0);
1485 231 : n = fabs(s + 2*p-1)/(2*M_PI);
1486 : }
1487 : else
1488 : {
1489 7 : p = 0;
1490 7 : n = exp((B - M_LN2) / s);
1491 : }
1492 : }
1493 8337 : else if (s <= 0 || t < 0.01) /* s < 0 may occur if s ~ 0 */
1494 3752 : { /* TODO: the crude bounds below are generally valid. Optimize ? */
1495 3752 : double l,l2, la = 1.; /* heuristic */
1496 3752 : double rlog, ilog; dcxlog(s-1,t, &rlog,&ilog);
1497 3752 : l2 = (s - 0.5)*rlog - t*ilog; /* = Re( (S - 1/2) log (S-1) ) */
1498 3752 : l = (B - l2 + s*log2PI) / (2. * (1.+ log((double)la)));
1499 3752 : l2 = dabs(s, t)/2;
1500 3752 : if (l < l2) l = l2;
1501 3752 : p = (long) ceil(l); if (p < 2) p = 2;
1502 3752 : n = 1 + dabs(p+s/2.-.25, t/2) * la / M_PI;
1503 : }
1504 : else
1505 : {
1506 4585 : double sn = dabs(s, t), L = log(sn/s);
1507 4585 : alpha = B - 0.39 + L + s*(log2PI - log(sn));
1508 4585 : beta = (alpha+s)/t - atan(s/t);
1509 4585 : p = 0;
1510 4585 : if (beta > 0)
1511 : {
1512 4578 : beta = 1.0 - s + t * get_xinf(beta);
1513 4578 : if (beta > 0) p = (long)ceil(beta / 2.0);
1514 : }
1515 : else
1516 7 : if (s < 1.0) p = 1;
1517 4585 : n = p? dabs(s + 2*p-1, t) / (2*M_PI) : exp((B-M_LN2+L) / s);
1518 : }
1519 8575 : *pp = p;
1520 8575 : *pn = (long)ceil(n);
1521 8575 : if (*pp < 0 || *pn < 0) pari_err_OVERFLOW("zeta");
1522 8575 : }
1523 :
1524 : /* zeta(a*j+b), j=0..N-1, b>1, using sumalt. Johansonn's thesis, Algo 4.7.1 */
1525 : static GEN
1526 723 : veczetas(long a, long b, long N, long prec)
1527 : {
1528 723 : const long n = ceil(2 + prec2nbits_mul(prec, M_LN2/1.7627));
1529 723 : pari_sp av = avma;
1530 723 : GEN c, d, z = zerovec(N);
1531 : long j, k;
1532 723 : c = d = int2n(2*n-1);
1533 85914 : for (k = n; k > 1; k--)
1534 : {
1535 85191 : GEN u, t = divii(d, powuu(k,b));
1536 85193 : if (!odd(k)) t = negi(t);
1537 85347 : gel(z,1) = addii(gel(z,1), t);
1538 85105 : u = powuu(k,a);
1539 3917942 : for (j = 1; j < N; j++)
1540 : {
1541 3877840 : t = divii(t,u); if (!signe(t)) break;
1542 3835536 : gel(z,j+1) = addii(gel(z,j+1), t);
1543 : }
1544 85085 : c = muluui(k,2*k-1,c);
1545 85102 : c = diviuuexact(c, 2*(n-k+1),n+k-1);
1546 85146 : d = addii(d,c);
1547 85130 : if (gc_needed(av,3))
1548 : {
1549 9 : if(DEBUGMEM>1) pari_warn(warnmem,"zetaBorwein, k = %ld", k);
1550 9 : gerepileall(av, 3, &c,&d,&z);
1551 : }
1552 : }
1553 : /* k = 1 */
1554 52436 : for (j = 1; j <= N; j++) gel(z,j) = addii(gel(z,j), d);
1555 692 : d = addiu(d, 1);
1556 52687 : for (j = 0, k = b - 1; j < N; j++, k += a)
1557 51964 : gel(z,j+1) = rdivii(shifti(gel(z,j+1), k), subii(shifti(d,k), d), prec);
1558 723 : return z;
1559 : }
1560 : /* zeta(a*j+b), j=0..N-1, b>1, using sumalt */
1561 : GEN
1562 217 : veczeta(GEN a, GEN b, long N, long prec)
1563 : {
1564 217 : pari_sp av = avma;
1565 : long n, j, k;
1566 217 : GEN L, c, d, z = zerovec(N);
1567 217 : if (typ(a) == t_INT && typ(b) == t_INT)
1568 126 : return gerepilecopy(av, veczetas(itos(a), itos(b), N, prec));
1569 91 : n = ceil(2 + prec2nbits_mul(prec, M_LN2/1.7627));
1570 91 : c = d = int2n(2*n-1);
1571 13819 : for (k = n; k; k--)
1572 : {
1573 : GEN u, t;
1574 13728 : L = logr_abs(utor(k, prec)); /* log(k) */
1575 13728 : t = gdiv(d, gexp(gmul(b, L), prec)); /* d / k^b */
1576 13728 : if (!odd(k)) t = gneg(t);
1577 13728 : gel(z,1) = gadd(gel(z,1), t);
1578 13728 : u = gexp(gmul(a, L), prec);
1579 896888 : for (j = 1; j < N; j++)
1580 : {
1581 889583 : t = gdiv(t,u); if (gexpo(t) < 0) break;
1582 883160 : gel(z,j+1) = gadd(gel(z,j+1), t);
1583 : }
1584 13728 : c = muluui(k,2*k-1,c);
1585 13728 : c = diviuuexact(c, 2*(n-k+1),n+k-1);
1586 13728 : d = addii(d,c);
1587 13728 : if (gc_needed(av,3))
1588 : {
1589 0 : if(DEBUGMEM>1) pari_warn(warnmem,"veczeta, k = %ld", k);
1590 0 : gerepileall(av, 3, &c,&d,&z);
1591 : }
1592 : }
1593 91 : L = mplog2(prec);
1594 5796 : for (j = 0; j < N; j++)
1595 : {
1596 5705 : GEN u = gsubgs(gadd(b, gmulgu(a,j)), 1);
1597 5705 : GEN w = gexp(gmul(u, L), prec); /* 2^u */
1598 5705 : gel(z,j+1) = gdiv(gmul(gel(z,j+1), w), gmul(d,gsubgs(w,1)));
1599 : }
1600 91 : return gerepilecopy(av, z);
1601 : }
1602 :
1603 : GEN
1604 26310 : constzeta(long n, long prec)
1605 : {
1606 26310 : GEN o = zetazone, z;
1607 26310 : long l = o? lg(o): 0;
1608 : pari_sp av;
1609 26310 : if (l > n)
1610 : {
1611 25864 : long p = realprec(gel(o,1));
1612 25864 : if (p >= prec) return o;
1613 : }
1614 597 : n = maxss(n, l + 15);
1615 597 : av = avma; z = veczetas(1, 2, n-1, prec);
1616 597 : zetazone = gclone(vec_prepend(z, mpeuler(prec)));
1617 597 : set_avma(av); guncloneNULL(o); return zetazone;
1618 : }
1619 :
1620 : /* zeta(s) using sumalt, case h=0,N=1. Assume s > 1 */
1621 : static GEN
1622 598 : zetaBorwein(long s, long prec)
1623 : {
1624 598 : pari_sp av = avma;
1625 598 : const long n = ceil(2 + prec2nbits_mul(prec, M_LN2/1.7627));
1626 : long k;
1627 598 : GEN c, d, z = gen_0;
1628 598 : c = d = int2n(2*n-1);
1629 135672 : for (k = n; k; k--)
1630 : {
1631 135074 : GEN t = divii(d, powuu(k,s));
1632 135074 : z = odd(k)? addii(z,t): subii(z,t);
1633 135074 : c = muluui(k,2*k-1,c);
1634 135074 : c = diviuuexact(c, 2*(n-k+1),n+k-1);
1635 135074 : d = addii(d,c);
1636 135074 : if (gc_needed(av,3))
1637 : {
1638 0 : if(DEBUGMEM>1) pari_warn(warnmem,"zetaBorwein, k = %ld", k);
1639 0 : gerepileall(av, 3, &c,&d,&z);
1640 : }
1641 : }
1642 598 : return rdivii(shifti(z, s-1), subii(shifti(d,s-1), d), prec);
1643 : }
1644 :
1645 : /* assume k != 1 */
1646 : GEN
1647 4739 : szeta(long k, long prec)
1648 : {
1649 4739 : pari_sp av = avma;
1650 : GEN z;
1651 :
1652 4739 : if (!k) { z = real2n(-1, prec); setsigne(z,-1); return z; }
1653 4732 : if (k < 0)
1654 : {
1655 252 : if (!odd(k)) return gen_0;
1656 : /* the one value such that k < 0 and 1 - k < 0, due to overflow */
1657 252 : if ((ulong)k == (HIGHBIT | 1))
1658 0 : pari_err_OVERFLOW("zeta [large negative argument]");
1659 252 : k = 1-k;
1660 252 : z = bernreal(k, prec); togglesign(z);
1661 252 : return gerepileuptoleaf(av, divru(z, k));
1662 : }
1663 : /* k > 1 */
1664 4480 : if (k > prec2nbits(prec)+1) return real_1(prec);
1665 4473 : if (zetazone && realprec(gel(zetazone,1)) >= prec && lg(zetazone) > k)
1666 119 : return rtor(gel(zetazone, k), prec);
1667 4354 : if (!odd(k))
1668 : {
1669 : GEN B;
1670 2583 : if (!bernzone) constbern(0);
1671 2583 : if (k < lg(bernzone))
1672 2128 : B = gel(bernzone, k>>1);
1673 : else
1674 : {
1675 455 : if (bernbitprec(k) > prec2nbits(prec))
1676 0 : return gerepileupto(av, invr(inv_szeta_euler(k, prec)));
1677 455 : B = bernfrac(k);
1678 : }
1679 : /* B = B_k */
1680 2583 : z = gmul(powru(Pi2n(1, prec + EXTRAPRECWORD), k), B);
1681 2583 : z = divrr(z, mpfactr(k,prec));
1682 2583 : setsigne(z, 1); shiftr_inplace(z, -1);
1683 : }
1684 : else
1685 : {
1686 1771 : double p = prec2nbits_mul(prec,0.393); /* bit / log_2(3+sqrt(8)) */
1687 1771 : p = log2(p * log(p));
1688 3542 : z = (p * k > prec2nbits(prec))? invr(inv_szeta_euler(k, prec))
1689 1771 : : zetaBorwein(k, prec);
1690 : }
1691 4354 : return gerepileuptoleaf(av, z);
1692 : }
1693 :
1694 : /* Ensure |1-s| >= 1/32 and (|s| <= 1/32 or real(s) >= 1/2) */
1695 : static int
1696 8589 : zeta_funeq(GEN *ps)
1697 : {
1698 8589 : GEN s = *ps, u;
1699 8589 : if (typ(s) == t_REAL)
1700 : {
1701 273 : u = subsr(1, s);
1702 273 : if (expo(u) >= -5
1703 266 : && ((signe(s) > 0 && expo(s) >= -1) || expo(s) <= -5)) return 0;
1704 : }
1705 : else
1706 : {
1707 8316 : GEN sig = gel(s,1);
1708 8316 : u = gsubsg(1, s);
1709 8316 : if (gexpo(u) >= -5
1710 8309 : && ((signe(sig) > 0 && expo(sig) >= -1) || gexpo(s) <= -5)) return 0;
1711 : }
1712 1505 : *ps = u; return 1;
1713 : }
1714 : /* s0 a t_INT, t_REAL or t_COMPLEX.
1715 : * If a t_INT, assume it's not a trivial case (i.e we have s0 > 1, odd) */
1716 : static GEN
1717 8596 : czeta(GEN s0, long prec)
1718 : {
1719 8596 : GEN ms, s, u, y, res, tes, sig, tau, invn2, ns, Ns, funeq_factor = NULL;
1720 : long i, nn, lim, lim2;
1721 8596 : pari_sp av0 = avma, av, av2;
1722 : pari_timer T;
1723 :
1724 8596 : if (DEBUGLEVEL>2) timer_start(&T);
1725 8596 : s = trans_fix_arg(&prec,&s0,&sig,&tau,&av,&res);
1726 8596 : if (typ(s0) == t_INT) return gerepileupto(av0, gzeta(s0, prec));
1727 8589 : if (zeta_funeq(&s)) /* s -> 1-s */
1728 : { /* Gamma(s) (2Pi)^-s 2 cos(Pi s/2) [new s] */
1729 1505 : GEN t = gmul(ggamma(s,prec), pow2Pis(gsubgs(s0,1), prec));
1730 1505 : sig = real_i(s);
1731 1505 : funeq_factor = gmul2n(gmul(t, gsin(gmul(Pi2n(-1,prec),s0), prec)), 1);
1732 : }
1733 8589 : if (gcmpgs(sig, prec2nbits(prec) + 1) > 0) { /* zeta(s) = 1 */
1734 14 : if (!funeq_factor) { set_avma(av0); return real_1(prec); }
1735 7 : return gerepileupto(av0, funeq_factor);
1736 : }
1737 8575 : optim_zeta(s, prec, &lim, &nn);
1738 8575 : if (DEBUGLEVEL>2) err_printf("lim, nn: [%ld, %ld]\n", lim, nn);
1739 :
1740 8575 : ms = gneg(s);
1741 8575 : if (umuluu_le(nn, prec, 10000000))
1742 : {
1743 8575 : incrprec(prec); /* one extra word of precision */
1744 8575 : Ns = vecpowug(nn, ms, prec);
1745 8575 : ns = gel(Ns,nn); setlg(Ns, nn);
1746 8575 : y = gadd(gmul2n(ns, -1), RgV_sum(Ns));
1747 : }
1748 : else
1749 : {
1750 0 : Ns = dirpowerssum(nn, ms, prec);
1751 0 : incrprec(prec); /* one extra word of precision */
1752 0 : ns = gpow(utor(nn, prec), ms, prec);
1753 0 : y = gsub(Ns, gmul2n(ns, -1));
1754 : }
1755 8575 : if (DEBUGLEVEL>2) timer_printf(&T,"sum from 1 to N");
1756 8575 : constbern(lim);
1757 8575 : if (DEBUGLEVEL>2) timer_start(&T);
1758 8575 : invn2 = divri(real_1(prec), sqru(nn)); lim2 = lim<<1;
1759 8575 : tes = bernfrac(lim2);
1760 : {
1761 : GEN s1, s2, s3, s4, s5;
1762 8575 : s2 = gmul(s, gsubgs(s,1));
1763 8575 : s3 = gmul2n(invn2,3);
1764 8575 : av2 = avma;
1765 8575 : s1 = gsubgs(gmul2n(s,1), 1);
1766 8575 : s4 = gmul(invn2, gmul2n(gaddsg(4*lim-2,s1),1));
1767 8575 : s5 = gmul(invn2, gadd(s2, gmulsg(lim2, gaddgs(s1, lim2))));
1768 805263 : for (i = lim2-2; i>=2; i -= 2)
1769 : {
1770 796688 : s5 = gsub(s5, s4);
1771 796688 : s4 = gsub(s4, s3);
1772 796688 : tes = gadd(bernfrac(i), gdivgunextu(gmul(s5,tes), i+1));
1773 796688 : if (gc_needed(av2,3))
1774 : {
1775 0 : if(DEBUGMEM>1) pari_warn(warnmem,"czeta i = %ld", i);
1776 0 : gerepileall(av2,3, &tes,&s5,&s4);
1777 : }
1778 : }
1779 8575 : u = gmul(gmul(tes,invn2), gmul2n(s2, -1));
1780 8575 : tes = gmulsg(nn, gaddsg(1, u));
1781 : }
1782 8575 : if (DEBUGLEVEL>2) timer_printf(&T,"Bernoulli sum");
1783 : /* y += tes n^(-s) / (s-1) */
1784 8575 : y = gadd(y, gmul(tes, gdiv(ns, gsubgs(s,1))));
1785 8575 : if (funeq_factor) y = gmul(y, funeq_factor);
1786 8575 : set_avma(av); return affc_fixlg(y,res);
1787 : }
1788 : /* v a t_VEC/t_COL */
1789 : int
1790 14 : RgV_is_arithprog(GEN v, GEN *a, GEN *b)
1791 : {
1792 14 : pari_sp av = avma, av2;
1793 14 : long i, n = lg(v)-1;
1794 14 : if (n == 0) { *a = *b = gen_0; return 1; }
1795 14 : *a = gel(v,1);
1796 14 : if (n == 1) { * b = gen_0; return 1; }
1797 14 : *b = gsub(gel(v,2), *a); av2 = avma;
1798 28 : for (i = 2; i < n; i++)
1799 14 : if (!gequal(*b, gsub(gel(v,i+1), gel(v,i)))) return gc_int(av,0);
1800 14 : return gc_int(av2,1);
1801 : }
1802 :
1803 : GEN
1804 13349 : gzeta(GEN x, long prec)
1805 : {
1806 13349 : pari_sp av = avma;
1807 : GEN y;
1808 13349 : if (gequal1(x)) pari_err_DOMAIN("zeta", "argument", "=", gen_1, x);
1809 13328 : switch(typ(x))
1810 : {
1811 4501 : case t_INT:
1812 4501 : if (is_bigint(x))
1813 : {
1814 21 : if (signe(x) > 0) return real_1(prec);
1815 14 : if (mod2(x) == 0) return real_0(prec);
1816 7 : pari_err_OVERFLOW("zeta [large negative argument]");
1817 : }
1818 4480 : return szeta(itos(x),prec);
1819 8596 : case t_REAL: case t_COMPLEX: return czeta(x,prec);
1820 14 : case t_PADIC: return Qp_zeta(x);
1821 21 : case t_VEC: case t_COL:
1822 : {
1823 : GEN a, b;
1824 21 : long n = lg(x) - 1;
1825 21 : if (n > 1 && RgV_is_arithprog(x, &a, &b))
1826 : {
1827 14 : if (!is_real_t(typ(a)) || !is_real_t(typ(b)) || gcmpgs(a, 1) <= 0)
1828 0 : { set_avma(av); break; }
1829 14 : a = veczeta(b, a, n, prec);
1830 14 : settyp(a, typ(x)); return a;
1831 : }
1832 : }
1833 : default:
1834 203 : if (!(y = toser_i(x))) break;
1835 35 : if (gequal1(y))
1836 7 : pari_err_DOMAIN("zeta", "argument", "=", gen_1, y);
1837 28 : return gerepileupto(av, lfun(gen_1,y,prec2nbits(prec)));
1838 : }
1839 168 : return trans_eval("zeta",gzeta,x,prec);
1840 : }
1841 :
1842 : /***********************************************************************/
1843 : /** **/
1844 : /** FONCTIONS POLYLOGARITHME **/
1845 : /** **/
1846 : /***********************************************************************/
1847 :
1848 : /* smallish k such that bernbitprec(K) > bit + Kdz, K = 2k+4 */
1849 : static long
1850 21 : get_k(double dz, long bit)
1851 : {
1852 : long a, b;
1853 21 : for (b = 128;; b <<= 1)
1854 21 : if (bernbitprec(b) > bit + b*dz) break;
1855 21 : if (b == 128) return 128;
1856 0 : a = b >> 1;
1857 0 : while (b - a > 64)
1858 : {
1859 0 : long c = (a+b) >> 1;
1860 0 : if (bernbitprec(c) > bit + c*dz) b = c; else a = c;
1861 : }
1862 0 : return b >> 1;
1863 : }
1864 :
1865 : /* m >= 2. Validity domain |log x| < 2*Pi, contains log |x| < 5.44,
1866 : * Li_m(x = e^z) = sum_{n >= 0} zeta(m-n) z^n / n!
1867 : * with zeta(1) := H_{m-1} - log(-z) */
1868 : static GEN
1869 21 : cxpolylog(long m, GEN x, long prec)
1870 : {
1871 : long li, n, k, ksmall, real;
1872 : GEN vz, z, Z, h, q, s, S;
1873 : pari_sp av;
1874 : double dz;
1875 : pari_timer T;
1876 :
1877 21 : if (gequal1(x)) return szeta(m,prec);
1878 : /* x real <= 1 ==> Li_m(x) real */
1879 21 : real = (typ(x) == t_REAL && (expo(x) < 0 || signe(x) <= 0));
1880 :
1881 21 : vz = constzeta(m, prec);
1882 21 : z = glog(x,prec);
1883 : /* n = 0 */
1884 21 : q = gen_1; s = gel(vz, m);
1885 28 : for (n=1; n < m-1; n++)
1886 : {
1887 7 : q = gdivgu(gmul(q,z),n);
1888 7 : s = gadd(s, gmul(gel(vz,m-n), real? real_i(q): q));
1889 : }
1890 : /* n = m-1 */
1891 21 : q = gdivgu(gmul(q,z),n); /* multiply by "zeta(1)" */
1892 21 : h = gmul(q, gsub(harmonic(m-1), glog(gneg_i(z),prec)));
1893 21 : s = gadd(s, real? real_i(h): h);
1894 : /* n = m */
1895 21 : q = gdivgu(gmul(q,z),m);
1896 21 : s = gadd(s, gdivgs(real? real_i(q): q, -2)); /* zeta(0) = -1/2 */
1897 : /* n = m+1 */
1898 21 : q = gdivgu(gmul(q,z),m+1); /* = z^(m+1) / (m+1)! */
1899 21 : s = gadd(s, gdivgs(real? real_i(q): q, -12)); /* zeta(-1) = -1/12 */
1900 :
1901 21 : li = -(prec2nbits(prec)+1);
1902 21 : if (DEBUGLEVEL) timer_start(&T);
1903 21 : dz = dbllog2(z) - log2PI; /* ~ log2(|z|/2Pi) */
1904 : /* sum_{k >= 1} zeta(-1-2k) * z^(2k+m+1) / (2k+m+1)!
1905 : * = 2 z^(m-1) sum_{k >= 1} zeta(2k+2) * Z^(k+1) / (2k+2)..(2k+1+m)), where
1906 : * Z = -(z/2Pi)^2. Stop at 2k = (li - (m-1)*Lz - m) / dz, Lz = log2 |z| */
1907 : /* We cut the sum in two: small values of k first */
1908 21 : Z = gsqr(z); av = avma;
1909 21 : ksmall = get_k(dz, prec2nbits(prec));
1910 21 : constbern(ksmall);
1911 469 : for(k = 1; k < ksmall; k++)
1912 : {
1913 469 : GEN t = q = gdivgunextu(gmul(q,Z), 2*k+m); /* z^(2k+m+1)/(2k+m+1)! */
1914 469 : if (real) t = real_i(t);
1915 469 : t = gmul(t, gdivgu(bernfrac(2*k+2), 2*k+2)); /* - t * zeta(1-(2k+2)) */
1916 469 : s = gsub(s, t);
1917 469 : if (gexpo(t) < li) return s;
1918 : /* large values ? */
1919 448 : if ((k & 0x1ff) == 0) gerepileall(av, 2, &s, &q);
1920 : }
1921 0 : if (DEBUGLEVEL>2) timer_printf(&T, "polylog: small k <= %ld", k);
1922 0 : Z = gneg(gsqr(gdiv(z, Pi2n(1,prec))));
1923 0 : q = gmul(gpowgs(z, m-1), gpowgs(Z, k+1)); /* Z^(k+1) * z^(m-1) */
1924 0 : S = gen_0; av = avma;
1925 0 : for(;; k++)
1926 0 : {
1927 0 : GEN t = q;
1928 : long b;
1929 0 : if (real) t = real_i(t);
1930 0 : b = prec + gexpo(t) / BITS_IN_LONG; /* decrease accuracy */
1931 0 : if (b == 2) break;
1932 : /* t * zeta(2k+2) / (2k+2)..(2k+1+m) */
1933 0 : t = gdiv(t, mulri(inv_szeta_euler(2*k+2, b),
1934 0 : mulu_interval(2*k+2, 2*k+1+m)));
1935 0 : S = gadd(S, t); if (gexpo(t) < li) break;
1936 0 : q = gmul(q, Z);
1937 0 : if ((k & 0x1ff) == 0) gerepileall(av, 2, &S, &q);
1938 : }
1939 0 : if (DEBUGLEVEL>2) timer_printf(&T, "polylog: large k <= %ld", k);
1940 0 : return gadd(s, gmul2n(S,1));
1941 : }
1942 :
1943 : static GEN
1944 42 : Li1(GEN x, long prec) { return gneg(glog(gsubsg(1, x), prec)); }
1945 :
1946 : static GEN
1947 203 : polylog(long m, GEN x, long prec)
1948 : {
1949 : long l, e, i, G, sx;
1950 : pari_sp av, av1;
1951 : GEN X, Xn, z, p1, p2, y, res;
1952 :
1953 203 : if (m < 0) pari_err_DOMAIN("polylog", "index", "<", gen_0, stoi(m));
1954 203 : if (!m) return mkfrac(gen_m1,gen_2);
1955 203 : if (gequal0(x)) return gcopy(x);
1956 203 : if (m==1) { av = avma; return gerepileupto(av, Li1(x, prec)); }
1957 :
1958 168 : l = precision(x);
1959 168 : if (!l) l = prec; else prec = l;
1960 168 : res = cgetc(l); av = avma;
1961 168 : x = gtofp(x, l+EXTRAPRECWORD);
1962 168 : e = gexpo(gnorm(x));
1963 168 : if (!e || e == -1) {
1964 21 : y = cxpolylog(m,x,prec);
1965 21 : set_avma(av); return affc_fixlg(y, res);
1966 : }
1967 147 : X = (e > 0)? ginv(x): x;
1968 147 : G = -prec2nbits(l);
1969 147 : av1 = avma;
1970 147 : y = Xn = X;
1971 147 : for (i=2; ; i++)
1972 : {
1973 68159 : Xn = gmul(X,Xn); p2 = gdiv(Xn,powuu(i,m));
1974 68159 : y = gadd(y,p2);
1975 68159 : if (gexpo(p2) <= G) break;
1976 :
1977 68012 : if (gc_needed(av1,1))
1978 : {
1979 0 : if(DEBUGMEM>1) pari_warn(warnmem,"polylog");
1980 0 : gerepileall(av1,2, &y, &Xn);
1981 : }
1982 : }
1983 147 : if (e < 0) { set_avma(av); return affc_fixlg(y, res); }
1984 :
1985 28 : sx = gsigne(imag_i(x));
1986 28 : if (!sx)
1987 : {
1988 28 : if (m&1) sx = gsigne(gsub(gen_1, real_i(x)));
1989 21 : else sx = - gsigne(real_i(x));
1990 : }
1991 28 : z = divri(mppi(l), mpfact(m-1)); setsigne(z, sx);
1992 28 : z = mkcomplex(gen_0, z);
1993 :
1994 28 : if (m == 2)
1995 : { /* same formula as below, written more efficiently */
1996 21 : y = gneg_i(y);
1997 21 : if (typ(x) == t_REAL && signe(x) < 0)
1998 7 : p1 = logr_abs(x);
1999 : else
2000 14 : p1 = gsub(glog(x,l), z);
2001 21 : p1 = gmul2n(gsqr(p1), -1); /* = (log(-x))^2 / 2 */
2002 :
2003 21 : p1 = gadd(p1, divru(sqrr(mppi(l)), 6));
2004 21 : p1 = gneg_i(p1);
2005 : }
2006 : else
2007 : {
2008 7 : GEN logx = glog(x,l), logx2 = gsqr(logx), vz = constzeta(m, l);
2009 7 : p1 = mkfrac(gen_m1,gen_2);
2010 14 : for (i = m-2; i >= 0; i -= 2)
2011 7 : p1 = gadd(gel(vz, m-i), gmul(p1, gdivgunextu(logx2, i+1)));
2012 7 : if (m&1) p1 = gmul(logx,p1); else y = gneg_i(y);
2013 7 : p1 = gadd(gmul2n(p1,1), gmul(z,gpowgs(logx,m-1)));
2014 7 : if (typ(x) == t_REAL && signe(x) < 0) p1 = real_i(p1);
2015 : }
2016 28 : y = gadd(y,p1);
2017 28 : set_avma(av); return affc_fixlg(y, res);
2018 : }
2019 : static GEN
2020 119 : RIpolylog(long m, GEN x, long real, long prec)
2021 : {
2022 119 : GEN y = polylog(m, x, prec);
2023 119 : return real? real_i(y): imag_i(y);
2024 : }
2025 : GEN
2026 21 : dilog(GEN x, long prec) { return gpolylog(2, x, prec); }
2027 :
2028 : /* x a floating point number, t_REAL or t_COMPLEX of t_REAL */
2029 : static GEN
2030 42 : logabs(GEN x)
2031 : {
2032 : GEN y;
2033 42 : if (typ(x) == t_COMPLEX)
2034 : {
2035 7 : y = logr_abs( cxnorm(x) );
2036 7 : shiftr_inplace(y, -1);
2037 : } else
2038 35 : y = logr_abs(x);
2039 42 : return y;
2040 : }
2041 :
2042 : static GEN
2043 21 : polylogD(long m, GEN x, long flag, long prec)
2044 : {
2045 21 : long fl = 0, k, l, m2;
2046 : pari_sp av;
2047 : GEN p1, p2, y;
2048 :
2049 21 : if (gequal0(x)) return gcopy(x);
2050 21 : m2 = m&1;
2051 21 : if (gequal1(x) && m>=2) return m2? szeta(m,prec): gen_0;
2052 21 : av = avma; l = precision(x);
2053 21 : if (!l) { l = prec; x = gtofp(x,l); }
2054 21 : p1 = logabs(x);
2055 21 : if (signe(p1) > 0) { x = ginv(x); fl = !m2; } else setabssign(p1);
2056 : /* |x| <= 1, p1 = - log|x| >= 0 */
2057 21 : p2 = gen_1;
2058 21 : y = RIpolylog(m, x, m2, l);
2059 84 : for (k = 1; k < m; k++)
2060 : {
2061 63 : GEN t = RIpolylog(m-k, x, m2, l);
2062 63 : p2 = gdivgu(gmul(p2,p1), k); /* (-log|x|)^k / k! */
2063 63 : y = gadd(y, gmul(p2, t));
2064 : }
2065 21 : if (m2)
2066 : {
2067 14 : p1 = flag? gdivgs(p1, -2*m): gdivgs(logabs(gsubsg(1,x)), m);
2068 14 : y = gadd(y, gmul(p2, p1));
2069 : }
2070 21 : if (fl) y = gneg(y);
2071 21 : return gerepileupto(av, y);
2072 : }
2073 :
2074 : static GEN
2075 14 : polylogP(long m, GEN x, long prec)
2076 : {
2077 14 : long fl = 0, k, l, m2;
2078 : pari_sp av;
2079 : GEN p1,y;
2080 :
2081 14 : if (gequal0(x)) return gcopy(x);
2082 14 : m2 = m&1;
2083 14 : if (gequal1(x) && m>=2) return m2? szeta(m,prec): gen_0;
2084 14 : av = avma; l = precision(x);
2085 14 : if (!l) { l = prec; x = gtofp(x,l); }
2086 14 : p1 = logabs(x);
2087 14 : if (signe(p1) > 0) { x = ginv(x); fl = !m2; setsigne(p1, -1); }
2088 : /* |x| <= 1 */
2089 14 : y = RIpolylog(m, x, m2, l);
2090 14 : if (m==1)
2091 : {
2092 7 : shiftr_inplace(p1, -1); /* log |x| / 2 */
2093 7 : y = gadd(y, p1);
2094 : }
2095 : else
2096 : { /* m >= 2, \sum_{0 <= k <= m} 2^k B_k/k! (log |x|)^k Li_{m-k}(x),
2097 : with Li_0(x) := -1/2 */
2098 7 : GEN u, t = RIpolylog(m-1, x, m2, l);
2099 7 : u = gneg_i(p1); /* u = 2 B1 log |x| */
2100 7 : y = gadd(y, gmul(u, t));
2101 7 : if (m > 2)
2102 : {
2103 : GEN p2;
2104 7 : shiftr_inplace(p1, 1); /* 2log|x| <= 0 */
2105 7 : constbern(m>>1);
2106 7 : p1 = sqrr(p1);
2107 7 : p2 = shiftr(p1,-1);
2108 21 : for (k = 2; k < m; k += 2)
2109 : {
2110 14 : if (k > 2) p2 = gdivgunextu(gmul(p2,p1),k-1); /* 2^k/k! log^k |x|*/
2111 14 : t = RIpolylog(m-k, x, m2, l);
2112 14 : u = gmul(p2, bernfrac(k));
2113 14 : y = gadd(y, gmul(u, t));
2114 : }
2115 : }
2116 : }
2117 14 : if (fl) y = gneg(y);
2118 14 : return gerepileupto(av, y);
2119 : }
2120 :
2121 : static GEN
2122 175 : gpolylog_i(void *E, GEN x, long prec)
2123 : {
2124 175 : pari_sp av = avma;
2125 175 : long i, n, v, m = (long)E;
2126 : GEN a, y;
2127 :
2128 175 : if (m <= 0)
2129 : {
2130 28 : a = gmul(x, poleval(eulerianpol(-m, 0), x));
2131 28 : return gerepileupto(av, gdiv(a, gpowgs(gsubsg(1, x), 1-m)));
2132 : }
2133 147 : switch(typ(x))
2134 : {
2135 84 : case t_REAL: case t_COMPLEX: return polylog(m,x,prec);
2136 7 : case t_INTMOD: case t_PADIC: pari_err_IMPL( "padic polylogarithm");
2137 56 : default:
2138 56 : av = avma; if (!(y = toser_i(x))) break;
2139 21 : if (!m) { set_avma(av); return mkfrac(gen_m1,gen_2); }
2140 21 : if (m==1) return gerepileupto(av, Li1(y, prec));
2141 21 : if (gequal0(y)) return gerepilecopy(av, y);
2142 21 : v = valp(y);
2143 21 : if (v < 0) pari_err_DOMAIN("polylog","valuation", "<", gen_0, x);
2144 14 : if (v > 0) {
2145 7 : n = (lg(y)-3 + v) / v;
2146 7 : a = zeroser(varn(y), lg(y)-2);
2147 35 : for (i=n; i>=1; i--)
2148 28 : a = gmul(y, gadd(a, powis(utoipos(i),-m)));
2149 : } else { /* v == 0 */
2150 7 : long vy = varn(y);
2151 7 : GEN a0 = polcoef_i(y, 0, -1), t = gdiv(derivser(y), y);
2152 7 : a = Li1(y, prec);
2153 14 : for (i=2; i<=m; i++)
2154 7 : a = gadd(gpolylog(i, a0, prec), integ(gmul(t, a), vy));
2155 : }
2156 14 : return gerepileupto(av, a);
2157 : }
2158 35 : return trans_evalgen("polylog", E, gpolylog_i, x, prec);
2159 : }
2160 : GEN
2161 133 : gpolylog(long m, GEN x, long prec) { return gpolylog_i((void*)m, x, prec); }
2162 :
2163 : GEN
2164 147 : polylog0(long m, GEN x, long flag, long prec)
2165 : {
2166 147 : switch(flag)
2167 : {
2168 105 : case 0: return gpolylog(m,x,prec);
2169 14 : case 1: return polylogD(m,x,0,prec);
2170 7 : case 2: return polylogD(m,x,1,prec);
2171 14 : case 3: return polylogP(m,x,prec);
2172 7 : default: pari_err_FLAG("polylog");
2173 : }
2174 : return NULL; /* LCOV_EXCL_LINE */
2175 : }
2176 :
2177 : GEN
2178 22411 : upper_to_cx(GEN x, long *prec)
2179 : {
2180 22411 : long tx = typ(x), l;
2181 22411 : if (tx == t_QUAD) { x = quadtofp(x, *prec); tx = typ(x); }
2182 22411 : switch(tx)
2183 : {
2184 22390 : case t_COMPLEX:
2185 22390 : if (gsigne(gel(x,2)) > 0) break; /*fall through*/
2186 : case t_REAL: case t_INT: case t_FRAC:
2187 14 : pari_err_DOMAIN("modular function", "Im(argument)", "<=", gen_0, x);
2188 7 : default:
2189 7 : pari_err_TYPE("modular function", x);
2190 : }
2191 22390 : l = precision(x); if (l) *prec = l;
2192 22390 : return x;
2193 : }
2194 :
2195 : /* sqrt(3)/2 */
2196 : static GEN
2197 2189 : sqrt32(long prec) { GEN z = sqrtr_abs(utor(3,prec)); setexpo(z, -1); return z; }
2198 : /* exp(i x), x = k pi/12 */
2199 : static GEN
2200 3532 : e12(ulong k, long prec)
2201 : {
2202 : int s, sPi, sPiov2;
2203 : GEN z, t;
2204 3532 : k %= 24;
2205 3532 : if (!k) return gen_1;
2206 3532 : if (k == 12) return gen_m1;
2207 3532 : if (k >12) { s = 1; k = 24 - k; } else s = 0; /* x -> 2pi - x */
2208 3532 : if (k > 6) { sPi = 1; k = 12 - k; } else sPi = 0; /* x -> pi - x */
2209 3532 : if (k > 3) { sPiov2 = 1; k = 6 - k; } else sPiov2 = 0; /* x -> pi/2 - x */
2210 3532 : z = cgetg(3, t_COMPLEX);
2211 3532 : switch(k)
2212 : {
2213 630 : case 0: gel(z,1) = icopy(gen_1); gel(z,2) = gen_0; break;
2214 782 : case 1: t = gmul2n(addrs(sqrt32(prec), 1), -1);
2215 782 : gel(z,1) = sqrtr(t);
2216 782 : gel(z,2) = gmul2n(invr(gel(z,1)), -2); break;
2217 :
2218 1407 : case 2: gel(z,1) = sqrt32(prec);
2219 1407 : gel(z,2) = real2n(-1, prec); break;
2220 :
2221 713 : case 3: gel(z,1) = sqrtr_abs(real2n(-1,prec));
2222 713 : gel(z,2) = rcopy(gel(z,1)); break;
2223 : }
2224 3532 : if (sPiov2) swap(gel(z,1), gel(z,2));
2225 3532 : if (sPi) togglesign(gel(z,1));
2226 3532 : if (s) togglesign(gel(z,2));
2227 3532 : return z;
2228 : }
2229 : /* z a t_FRAC */
2230 : static GEN
2231 12548 : expIPifrac(GEN z, long prec)
2232 : {
2233 12548 : GEN n = gel(z,1), d = gel(z,2);
2234 12548 : ulong r, q = uabsdivui_rem(12, d, &r);
2235 12548 : if (!r) return e12(q * umodiu(n, 24), prec); /* 12 | d */
2236 9016 : n = centermodii(n, shifti(d,1), d);
2237 9016 : return expIr(divri(mulri(mppi(prec), n), d));
2238 : }
2239 : /* exp(i Pi z), z a t_INT or t_FRAC */
2240 : static GEN
2241 10997 : expIPiQ(GEN z, long prec)
2242 : {
2243 10997 : if (typ(z) == t_INT) return mpodd(z)? gen_m1: gen_1;
2244 1958 : return expIPifrac(z, prec);
2245 : }
2246 :
2247 : /* convert power of 2 t_REAL to rational */
2248 : static GEN
2249 1719 : real2nQ(GEN x)
2250 : {
2251 1719 : long e = expo(x);
2252 : GEN z;
2253 1719 : if (e < 0)
2254 426 : z = mkfrac(signe(x) < 0? gen_m1: gen_1, int2n(-e));
2255 : else
2256 : {
2257 1293 : z = int2n(e);
2258 1293 : if (signe(x) < 0) togglesign_safe(&z);
2259 : }
2260 1719 : return z;
2261 : }
2262 : /* x a real number */
2263 : GEN
2264 180775 : expIPiR(GEN x, long prec)
2265 : {
2266 180775 : if (typ(x) == t_REAL && absrnz_equal2n(x)) x = real2nQ(x);
2267 180775 : switch(typ(x))
2268 : {
2269 3036 : case t_INT: return mpodd(x)? gen_m1: gen_1;
2270 1941 : case t_FRAC: return expIPifrac(x, prec);
2271 : }
2272 175798 : return expIr(mulrr(mppi(prec), x));
2273 : }
2274 : /* z a t_COMPLEX */
2275 : GEN
2276 232740 : expIPiC(GEN z, long prec)
2277 : {
2278 : GEN pi, r, x, y;
2279 232740 : if (typ(z) != t_COMPLEX) return expIPiR(z, prec);
2280 52252 : x = gel(z,1);
2281 52252 : y = gel(z,2); if (gequal0(y)) return expIPiR(x, prec);
2282 51965 : pi = mppi(prec);
2283 51965 : r = gmul(pi, y); togglesign(r); r = mpexp(r); /* exp(-pi y) */
2284 51965 : if (typ(x) == t_REAL && absrnz_equal2n(x)) x = real2nQ(x);
2285 51965 : switch(typ(x))
2286 : {
2287 4550 : case t_INT: if (mpodd(x)) togglesign(r);
2288 4550 : return r;
2289 8649 : case t_FRAC: return gmul(r, expIPifrac(x, prec));
2290 : }
2291 38766 : return gmul(r, expIr(mulrr(pi, x)));
2292 : }
2293 : /* exp(I x y), more efficient for x in R, y pure imaginary */
2294 : GEN
2295 147 : expIxy(GEN x, GEN y, long prec) { return gexp(gmul(x, mulcxI(y)), prec); }
2296 :
2297 : static GEN
2298 13549 : qq(GEN x, long prec)
2299 : {
2300 13549 : long tx = typ(x);
2301 : GEN y;
2302 :
2303 13549 : if (is_scalar_t(tx))
2304 : {
2305 13507 : if (tx == t_PADIC) return x;
2306 13493 : x = upper_to_cx(x, &prec);
2307 13479 : return cxtoreal(expIPiC(gmul2n(x,1), prec)); /* e(x) */
2308 : }
2309 42 : if (! ( y = toser_i(x)) ) pari_err_TYPE("modular function", x);
2310 42 : return y;
2311 : }
2312 :
2313 : /* return (y * X^d) + x. Assume d > 0, x != 0, valp(x) = 0 */
2314 : static GEN
2315 21 : ser_addmulXn(GEN y, GEN x, long d)
2316 : {
2317 21 : long i, lx, ly, l = valp(y) + d; /* > 0 */
2318 : GEN z;
2319 :
2320 21 : lx = lg(x);
2321 21 : ly = lg(y) + l; if (lx < ly) ly = lx;
2322 21 : if (l > lx-2) return gcopy(x);
2323 21 : z = cgetg(ly,t_SER);
2324 77 : for (i=2; i<=l+1; i++) gel(z,i) = gel(x,i);
2325 70 : for ( ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-l));
2326 21 : z[1] = x[1]; return z;
2327 : }
2328 :
2329 : /* q a t_POL s.t. q(0) != 0, v > 0, Q = x^v*q; return \prod_i (1-Q^i) */
2330 : static GEN
2331 28 : RgXn_eta(GEN q, long v, long lim)
2332 : {
2333 28 : pari_sp av = avma;
2334 : GEN qn, ps, y;
2335 : ulong vps, vqn, n;
2336 :
2337 28 : if (!degpol(q) && isint1(gel(q,2))) return eta_ZXn(v, lim+v);
2338 7 : y = qn = ps = pol_1(0);
2339 7 : vps = vqn = 0;
2340 7 : for(n = 0;; n++)
2341 7 : { /* qn = q^n, ps = (-1)^n q^(n(3n+1)/2),
2342 : * vps, vqn valuation of ps, qn HERE */
2343 14 : pari_sp av2 = avma;
2344 14 : ulong vt = vps + 2*vqn + v; /* valuation of t at END of loop body */
2345 : long k1, k2;
2346 : GEN t;
2347 14 : vqn += v; vps = vt + vqn; /* valuation of qn, ps at END of body */
2348 14 : k1 = lim + v - vt + 1;
2349 14 : k2 = k1 - vqn; /* = lim + v - vps + 1 */
2350 14 : if (k1 <= 0) break;
2351 14 : t = RgX_mul(q, RgX_sqr(qn));
2352 14 : t = RgXn_red_shallow(t, k1);
2353 14 : t = RgX_mul(ps,t);
2354 14 : t = RgXn_red_shallow(t, k1);
2355 14 : t = RgX_neg(t); /* t = (-1)^(n+1) q^(n(3n+1)/2 + 2n+1) */
2356 14 : t = gerepileupto(av2, t);
2357 14 : y = RgX_addmulXn_shallow(t, y, vt);
2358 14 : if (k2 <= 0) break;
2359 :
2360 7 : qn = RgX_mul(qn,q);
2361 7 : ps = RgX_mul(t,qn);
2362 7 : ps = RgXn_red_shallow(ps, k2);
2363 7 : y = RgX_addmulXn_shallow(ps, y, vps);
2364 :
2365 7 : if (gc_needed(av,1))
2366 : {
2367 0 : if(DEBUGMEM>1) pari_warn(warnmem,"eta, n = %ld", n);
2368 0 : gerepileall(av, 3, &y, &qn, &ps);
2369 : }
2370 : }
2371 7 : return y;
2372 : }
2373 :
2374 : static GEN
2375 16455 : inteta(GEN q)
2376 : {
2377 16455 : long tx = typ(q);
2378 : GEN ps, qn, y;
2379 :
2380 16455 : y = gen_1; qn = gen_1; ps = gen_1;
2381 16455 : if (tx==t_PADIC)
2382 : {
2383 28 : if (valp(q) <= 0) pari_err_DOMAIN("eta", "v_p(q)", "<=",gen_0,q);
2384 : for(;;)
2385 56 : {
2386 77 : GEN t = gneg_i(gmul(ps,gmul(q,gsqr(qn))));
2387 77 : y = gadd(y,t); qn = gmul(qn,q); ps = gmul(t,qn);
2388 77 : t = y;
2389 77 : y = gadd(y,ps); if (gequal(t,y)) return y;
2390 : }
2391 : }
2392 :
2393 16427 : if (tx == t_SER)
2394 : {
2395 : ulong vps, vqn;
2396 42 : long l = lg(q), v, n;
2397 : pari_sp av;
2398 :
2399 42 : v = valp(q); /* handle valuation separately to avoid overflow */
2400 42 : if (v <= 0) pari_err_DOMAIN("eta", "v_p(q)", "<=",gen_0,q);
2401 35 : y = ser2pol_i(q, l); /* t_SER inefficient when input has low degree */
2402 35 : n = degpol(y);
2403 35 : if (n <= (l>>2))
2404 : {
2405 28 : GEN z = RgXn_eta(y, v, l-2);
2406 28 : setvarn(z, varn(y)); return RgX_to_ser(z, l+v);
2407 : }
2408 7 : q = leafcopy(q); av = avma;
2409 7 : setvalp(q, 0);
2410 7 : y = scalarser(gen_1, varn(q), l+v);
2411 7 : vps = vqn = 0;
2412 7 : for(n = 0;; n++)
2413 7 : { /* qn = q^n, ps = (-1)^n q^(n(3n+1)/2) */
2414 14 : ulong vt = vps + 2*vqn + v;
2415 : long k;
2416 : GEN t;
2417 14 : t = gneg_i(gmul(ps,gmul(q,gsqr(qn))));
2418 : /* t = (-1)^(n+1) q^(n(3n+1)/2 + 2n+1) */
2419 14 : y = ser_addmulXn(t, y, vt);
2420 14 : vqn += v; vps = vt + vqn;
2421 14 : k = l+v - vps; if (k <= 2) return y;
2422 :
2423 7 : qn = gmul(qn,q); ps = gmul(t,qn);
2424 7 : y = ser_addmulXn(ps, y, vps);
2425 7 : setlg(q, k);
2426 7 : setlg(qn, k);
2427 7 : setlg(ps, k);
2428 7 : if (gc_needed(av,3))
2429 : {
2430 0 : if(DEBUGMEM>1) pari_warn(warnmem,"eta");
2431 0 : gerepileall(av, 3, &y, &qn, &ps);
2432 : }
2433 : }
2434 : }
2435 : {
2436 16385 : long l = -prec2nbits(precision(q));
2437 16385 : pari_sp av = avma;
2438 :
2439 : for(;;)
2440 53573 : {
2441 69958 : GEN t = gneg_i(gmul(ps,gmul(q,gsqr(qn))));
2442 : /* qn = q^n
2443 : * ps = (-1)^n q^(n(3n+1)/2)
2444 : * t = (-1)^(n+1) q^(n(3n+1)/2 + 2n+1) */
2445 69958 : y = gadd(y,t); qn = gmul(qn,q); ps = gmul(t,qn);
2446 69958 : y = gadd(y,ps);
2447 69958 : if (gexpo(ps)-gexpo(y) < l) return y;
2448 53573 : if (gc_needed(av,3))
2449 : {
2450 0 : if(DEBUGMEM>1) pari_warn(warnmem,"eta");
2451 0 : gerepileall(av, 3, &y, &qn, &ps);
2452 : }
2453 : }
2454 : }
2455 : }
2456 :
2457 : GEN
2458 77 : eta(GEN x, long prec)
2459 : {
2460 77 : pari_sp av = avma;
2461 77 : GEN z = inteta( qq(x,prec) );
2462 49 : if (typ(z) == t_SER) return gerepilecopy(av, z);
2463 14 : return gerepileupto(av, z);
2464 : }
2465 :
2466 : /* s(h,k) = sum(n = 1, k-1, (n/k)*(frac(h*n/k) - 1/2))
2467 : * Knuth's algorithm. h integer, k integer > 0, (h,k) = 1 */
2468 : GEN
2469 6258 : sumdedekind_coprime(GEN h, GEN k)
2470 : {
2471 6258 : pari_sp av = avma;
2472 : GEN s2, s1, p, pp;
2473 : long s;
2474 6258 : if (lgefint(k) == 3 && uel(k,2) <= (2*(ulong)LONG_MAX) / 3)
2475 : {
2476 6251 : ulong kk = k[2], hh = umodiu(h, kk);
2477 : long s1, s2;
2478 : GEN v;
2479 6251 : if (signe(k) < 0) { k = negi(k); hh = Fl_neg(hh, kk); }
2480 6251 : v = u_sumdedekind_coprime(hh, kk);
2481 6251 : s1 = v[1]; s2 = v[2];
2482 6251 : return gerepileupto(av, gdiv(addis(mulis(k,s1), s2), muluu(12, kk)));
2483 : }
2484 7 : s = 1;
2485 7 : s1 = gen_0; p = gen_1; pp = gen_0;
2486 7 : s2 = h = modii(h, k);
2487 35 : while (signe(h)) {
2488 28 : GEN r, nexth, a = dvmdii(k, h, &nexth);
2489 28 : if (is_pm1(h)) s2 = s == 1? addii(s2, p): subii(s2, p);
2490 28 : s1 = s == 1? addii(s1, a): subii(s1, a);
2491 28 : s = -s;
2492 28 : k = h; h = nexth;
2493 28 : r = addii(mulii(a,p), pp); pp = p; p = r;
2494 : }
2495 : /* at this point p = original k */
2496 7 : if (s == -1) s1 = subiu(s1, 3);
2497 7 : return gerepileupto(av, gdiv(addii(mulii(p,s1), s2), muliu(p,12)));
2498 : }
2499 : /* as above, for ulong arguments.
2500 : * k integer > 0, 0 <= h < k, (h,k) = 1. Returns [s1,s2] such that
2501 : * s(h,k) = (s2 + k s1) / (12k). Requires max(h + k/2, k) < LONG_MAX
2502 : * to avoid overflow, in particular k <= LONG_MAX * 2/3 is fine */
2503 : GEN
2504 6251 : u_sumdedekind_coprime(long h, long k)
2505 : {
2506 6251 : long s = 1, s1 = 0, s2 = h, p = 1, pp = 0;
2507 11424 : while (h) {
2508 5173 : long r, nexth = k % h, a = k / h; /* a >= 1, a >= 2 if h = 1 */
2509 5173 : if (h == 1) s2 += p * s; /* occurs exactly once, last step */
2510 5173 : s1 += a * s;
2511 5173 : s = -s;
2512 5173 : k = h; h = nexth;
2513 5173 : r = a*p + pp; pp = p; p = r; /* p >= pp >= 0 */
2514 : }
2515 : /* in the above computation, p increases from 1 to original k,
2516 : * -k/2 <= s2 <= h + k/2, and |s1| <= k */
2517 6251 : if (s < 0) s1 -= 3; /* |s1| <= k+3 ? */
2518 : /* But in fact, |s2 + p s1| <= k^2 + 1/2 - 3k; if (s < 0), we have
2519 : * |s2| <= k/2 and it follows that |s1| < k here as well */
2520 : /* p = k; s(h,k) = (s2 + p s1)/12p. */
2521 6251 : return mkvecsmall2(s1, s2);
2522 : }
2523 : GEN
2524 28 : sumdedekind(GEN h, GEN k)
2525 : {
2526 28 : pari_sp av = avma;
2527 : GEN d;
2528 28 : if (typ(h) != t_INT) pari_err_TYPE("sumdedekind",h);
2529 28 : if (typ(k) != t_INT) pari_err_TYPE("sumdedekind",k);
2530 28 : d = gcdii(h,k);
2531 28 : if (!is_pm1(d))
2532 : {
2533 7 : h = diviiexact(h, d);
2534 7 : k = diviiexact(k, d);
2535 : }
2536 28 : return gerepileupto(av, sumdedekind_coprime(h,k));
2537 : }
2538 :
2539 : /* eta(x); assume Im x >> 0 (e.g. x in SL2's standard fundamental domain) */
2540 : static GEN
2541 17297 : eta_reduced(GEN x, long prec)
2542 : {
2543 17297 : GEN z = expIPiC(gdivgu(x, 12), prec); /* e(x/24) */
2544 17297 : if (24 * gexpo(z) >= -prec2nbits(prec))
2545 16364 : z = gmul(z, inteta( gpowgs(z,24) ));
2546 17297 : return z;
2547 : }
2548 :
2549 : /* x = U.z (flag = 1), or x = U^(-1).z (flag = 0)
2550 : * Return [s,t] such that eta(z) = eta(x) * sqrt(s) * exp(I Pi t) */
2551 : static GEN
2552 17311 : eta_correction(GEN x, GEN U, long flag)
2553 : {
2554 : GEN a,b,c,d, s,t;
2555 : long sc;
2556 17311 : a = gcoeff(U,1,1);
2557 17311 : b = gcoeff(U,1,2);
2558 17311 : c = gcoeff(U,2,1);
2559 17311 : d = gcoeff(U,2,2);
2560 : /* replace U by U^(-1) */
2561 17311 : if (flag) {
2562 8939 : swap(a,d);
2563 8939 : togglesign_safe(&b);
2564 8939 : togglesign_safe(&c);
2565 : }
2566 17311 : sc = signe(c);
2567 17311 : if (!sc) {
2568 11081 : if (signe(d) < 0) togglesign_safe(&b);
2569 11081 : s = gen_1;
2570 11081 : t = uutoQ(umodiu(b, 24), 12);
2571 : } else {
2572 6230 : if (sc < 0) {
2573 1722 : togglesign_safe(&a);
2574 1722 : togglesign_safe(&b);
2575 1722 : togglesign_safe(&c);
2576 1722 : togglesign_safe(&d);
2577 : } /* now c > 0 */
2578 6230 : s = mulcxmI(gadd(gmul(c,x), d));
2579 6230 : t = gadd(gdiv(addii(a,d),muliu(c,12)), sumdedekind_coprime(negi(d),c));
2580 : /* correction : exp(I Pi (((a+d)/12c) + s(-d,c)) ) sqrt(-i(cx+d)) */
2581 : }
2582 17311 : return mkvec2(s, t);
2583 : }
2584 :
2585 : /* returns the true value of eta(x) for Im(x) > 0, using reduction to
2586 : * standard fundamental domain */
2587 : GEN
2588 8869 : trueeta(GEN x, long prec)
2589 : {
2590 8869 : pari_sp av = avma;
2591 : GEN U, st, s, t;
2592 :
2593 8869 : if (!is_scalar_t(typ(x))) pari_err_TYPE("trueeta",x);
2594 8869 : x = upper_to_cx(x, &prec);
2595 8869 : x = cxredsl2(x, &U);
2596 8869 : st = eta_correction(x, U, 1);
2597 8869 : x = eta_reduced(x, prec);
2598 8869 : s = gel(st, 1);
2599 8869 : t = gel(st, 2);
2600 8869 : x = gmul(x, expIPiQ(t, prec));
2601 8869 : if (s != gen_1) x = gmul(x, gsqrt(s, prec));
2602 8869 : return gerepileupto(av, x);
2603 : }
2604 :
2605 : GEN
2606 77 : eta0(GEN x, long flag,long prec)
2607 77 : { return flag? trueeta(x,prec): eta(x,prec); }
2608 :
2609 : /* eta(q) = 1 + \sum_{n>0} (-1)^n * (q^(n(3n-1)/2) + q^(n(3n+1)/2)) */
2610 : static GEN
2611 7 : ser_eta(long prec)
2612 : {
2613 7 : GEN e = cgetg(prec+2, t_SER), ed = e+2;
2614 : long n, j;
2615 7 : e[1] = evalsigne(1)|_evalvalp(0)|evalvarn(0);
2616 7 : gel(ed,0) = gen_1;
2617 483 : for (n = 1; n < prec; n++) gel(ed,n) = gen_0;
2618 49 : for (n = 1, j = 0; n < prec; n++)
2619 : {
2620 : GEN s;
2621 49 : j += 3*n-2; /* = n*(3*n-1) / 2 */;
2622 49 : if (j >= prec) break;
2623 42 : s = odd(n)? gen_m1: gen_1;
2624 42 : gel(ed, j) = s;
2625 42 : if (j+n >= prec) break;
2626 42 : gel(ed, j+n) = s;
2627 : }
2628 7 : return e;
2629 : }
2630 :
2631 : static GEN
2632 476 : coeffEu(GEN fa)
2633 : {
2634 476 : pari_sp av = avma;
2635 476 : return gerepileuptoint(av, mului(65520, usumdivk_fact(fa,11)));
2636 : }
2637 : /* E12 = 1 + q*E/691 */
2638 : static GEN
2639 7 : ser_E(long prec)
2640 : {
2641 7 : GEN e = cgetg(prec+2, t_SER), ed = e+2;
2642 7 : GEN F = vecfactoru_i(2, prec); /* F[n] = factoru(n+1) */
2643 : long n;
2644 7 : e[1] = evalsigne(1)|_evalvalp(0)|evalvarn(0);
2645 7 : gel(ed,0) = utoipos(65520);
2646 483 : for (n = 1; n < prec; n++) gel(ed,n) = coeffEu(gel(F,n));
2647 7 : return e;
2648 : }
2649 : /* j = E12/Delta + 432000/691, E12 = 1 + q*E/691 */
2650 : static GEN
2651 7 : ser_j2(long prec, long v)
2652 : {
2653 7 : pari_sp av = avma;
2654 7 : GEN iD = gpowgs(ginv(ser_eta(prec)), 24); /* q/Delta */
2655 7 : GEN J = gmul(ser_E(prec), iD);
2656 7 : setvalp(iD,-1); /* now 1/Delta */
2657 7 : J = gadd(gdivgu(J, 691), iD);
2658 7 : J = gerepileupto(av, J);
2659 7 : if (prec > 1) gel(J,3) = utoipos(744);
2660 7 : setvarn(J,v); return J;
2661 : }
2662 :
2663 : /* j(q) = \sum_{n >= -1} c(n)q^n,
2664 : * \sum_{n = -1}^{N-1} c(n) (-10n \sigma_3(N-n) + 21 \sigma_5(N-n))
2665 : * = c(N) (N+1)/24 */
2666 : static GEN
2667 14 : ser_j(long prec, long v)
2668 : {
2669 : GEN j, J, S3, S5, F;
2670 : long i, n;
2671 14 : if (prec > 64) return ser_j2(prec, v);
2672 7 : S3 = cgetg(prec+1, t_VEC);
2673 7 : S5 = cgetg(prec+1,t_VEC);
2674 7 : F = vecfactoru_i(1, prec);
2675 35 : for (n = 1; n <= prec; n++)
2676 : {
2677 28 : GEN fa = gel(F,n);
2678 28 : gel(S3,n) = mului(10, usumdivk_fact(fa,3));
2679 28 : gel(S5,n) = mului(21, usumdivk_fact(fa,5));
2680 : }
2681 7 : J = cgetg(prec+2, t_SER),
2682 7 : J[1] = evalvarn(v)|evalsigne(1)|evalvalp(-1);
2683 7 : j = J+3;
2684 7 : gel(j,-1) = gen_1;
2685 7 : gel(j,0) = utoipos(744);
2686 7 : gel(j,1) = utoipos(196884);
2687 21 : for (n = 2; n < prec; n++)
2688 : {
2689 14 : pari_sp av = avma;
2690 14 : GEN c, s3 = gel(S3,n+1), s5 = gel(S5,n+1);
2691 14 : c = addii(s3, s5);
2692 49 : for (i = 0; i < n; i++)
2693 : {
2694 35 : s3 = gel(S3,n-i); s5 = gel(S5,n-i);
2695 35 : c = addii(c, mulii(gel(j,i), subii(s5, mului(i,s3))));
2696 : }
2697 14 : gel(j,n) = gerepileuptoint(av, diviuexact(muliu(c,24), n+1));
2698 : }
2699 7 : return J;
2700 : }
2701 :
2702 : GEN
2703 42 : jell(GEN x, long prec)
2704 : {
2705 42 : long tx = typ(x);
2706 42 : pari_sp av = avma;
2707 : GEN q, h, U;
2708 :
2709 42 : if (!is_scalar_t(tx))
2710 : {
2711 : long v;
2712 21 : if (gequalX(x)) return ser_j(precdl, varn(x));
2713 21 : q = toser_i(x); if (!q) pari_err_TYPE("ellj",x);
2714 14 : v = fetch_var_higher();
2715 14 : h = ser_j(lg(q)-2, v);
2716 14 : h = gsubst(h, v, q);
2717 14 : delete_var(); return gerepileupto(av, h);
2718 : }
2719 21 : if (tx == t_PADIC)
2720 : {
2721 7 : GEN p2, p1 = gdiv(inteta(gsqr(x)), inteta(x));
2722 7 : p1 = gmul2n(gsqr(p1),1);
2723 7 : p1 = gmul(x,gpowgs(p1,12));
2724 7 : p2 = gaddsg(768,gadd(gsqr(p1),gdivsg(4096,p1)));
2725 7 : p1 = gmulsg(48,p1);
2726 7 : return gerepileupto(av, gadd(p2,p1));
2727 : }
2728 : /* Let h = Delta(2x) / Delta(x), then j(x) = (1 + 256h)^3 / h */
2729 14 : x = upper_to_cx(x, &prec);
2730 7 : x = cxredsl2(x, &U); /* forget about Ua : j has weight 0 */
2731 : { /* cf eta_reduced, raised to power 24
2732 : * Compute
2733 : * t = (inteta(q(2x)) / inteta(q(x))) ^ 24;
2734 : * then
2735 : * h = t * (q(2x) / q(x) = t * q(x);
2736 : * but inteta(q) costly and useless if expo(q) << 1 => inteta(q) = 1.
2737 : * log_2 ( exp(-2Pi Im tau) ) < -prec2nbits(prec)
2738 : * <=> Im tau > prec2nbits(prec) * log(2) / 2Pi */
2739 7 : long C = (long)prec2nbits_mul(prec, M_LN2/(2*M_PI));
2740 7 : q = expIPiC(gmul2n(x,1), prec); /* e(x) */
2741 7 : if (gcmpgs(gel(x,2), C) > 0) /* eta(q(x)) = 1 : no need to compute q(2x) */
2742 0 : h = q;
2743 : else
2744 : {
2745 7 : GEN t = gdiv(inteta(gsqr(q)), inteta(q));
2746 7 : h = gmul(q, gpowgs(t, 24));
2747 : }
2748 : }
2749 : /* real_1 important ! gaddgs(, 1) could increase the accuracy ! */
2750 7 : return gerepileupto(av, gdiv(gpowgs(gadd(gmul2n(h,8), real_1(prec)), 3), h));
2751 : }
2752 :
2753 : static GEN
2754 8372 : to_form(GEN a, GEN w, GEN C)
2755 8372 : { return mkvec3(a, w, diviiexact(C, a)); }
2756 : static GEN
2757 8372 : form_to_quad(GEN f, GEN sqrtD)
2758 : {
2759 8372 : long a = itos(gel(f,1)), a2 = a << 1;
2760 8372 : GEN b = gel(f,2);
2761 8372 : return mkcomplex(gdivgs(b, -a2), gdivgs(sqrtD, a2));
2762 : }
2763 : static GEN
2764 8372 : eta_form(GEN f, GEN sqrtD, GEN *s_t, long prec)
2765 : {
2766 8372 : GEN U, t = form_to_quad(redimagsl2(f, &U), sqrtD);
2767 8372 : *s_t = eta_correction(t, U, 0);
2768 8372 : return eta_reduced(t, prec);
2769 : }
2770 :
2771 : /* eta(t/p)eta(t/q) / (eta(t)eta(t/pq)), t = (-w + sqrt(D)) / 2a */
2772 : GEN
2773 2093 : double_eta_quotient(GEN a, GEN w, GEN D, long p, long q, GEN pq, GEN sqrtD)
2774 : {
2775 2093 : GEN C = shifti(subii(sqri(w), D), -2);
2776 : GEN d, t, z, zp, zq, zpq, s_t, s_tp, s_tpq, s, sp, spq;
2777 2093 : long prec = realprec(sqrtD);
2778 :
2779 2093 : z = eta_form(to_form(a, w, C), sqrtD, &s_t, prec);
2780 2093 : s = gel(s_t, 1);
2781 2093 : zp = eta_form(to_form(mului(p, a), w, C), sqrtD, &s_tp, prec);
2782 2093 : sp = gel(s_tp, 1);
2783 2093 : zpq = eta_form(to_form(mulii(pq, a), w, C), sqrtD, &s_tpq, prec);
2784 2093 : spq = gel(s_tpq, 1);
2785 2093 : if (p == q) {
2786 0 : z = gdiv(gsqr(zp), gmul(z, zpq));
2787 0 : t = gsub(gmul2n(gel(s_tp,2), 1),
2788 0 : gadd(gel(s_t,2), gel(s_tpq,2)));
2789 0 : if (sp != gen_1) z = gmul(z, sp);
2790 : } else {
2791 : GEN s_tq, sq;
2792 2093 : zq = eta_form(to_form(mului(q, a), w, C), sqrtD, &s_tq, prec);
2793 2093 : sq = gel(s_tq, 1);
2794 2093 : z = gdiv(gmul(zp, zq), gmul(z, zpq));
2795 2093 : t = gsub(gadd(gel(s_tp,2), gel(s_tq,2)),
2796 2093 : gadd(gel(s_t,2), gel(s_tpq,2)));
2797 2093 : if (sp != gen_1) z = gmul(z, gsqrt(sp, prec));
2798 2093 : if (sq != gen_1) z = gmul(z, gsqrt(sq, prec));
2799 : }
2800 2093 : d = NULL;
2801 2093 : if (s != gen_1) d = gsqrt(s, prec);
2802 2093 : if (spq != gen_1) {
2803 2065 : GEN x = gsqrt(spq, prec);
2804 2065 : d = d? gmul(d, x): x;
2805 : }
2806 2093 : if (d) z = gdiv(z, d);
2807 2093 : return gmul(z, expIPiQ(t, prec));
2808 : }
2809 :
2810 : typedef struct { GEN u; long v, t; } cxanalyze_t;
2811 :
2812 : /* Check whether a t_COMPLEX, t_REAL or t_INT z != 0 can be written as
2813 : * z = u * 2^(v/2) * exp(I Pi/4 t), u > 0, v = 0,1 and -3 <= t <= 4.
2814 : * Allow z t_INT/t_REAL to simplify handling of eta_correction() output */
2815 : static int
2816 70 : cxanalyze(cxanalyze_t *T, GEN z)
2817 : {
2818 : GEN a, b;
2819 : long ta, tb;
2820 :
2821 70 : T->v = 0;
2822 70 : if (is_intreal_t(typ(z)))
2823 : {
2824 63 : T->u = mpabs_shallow(z);
2825 63 : T->t = signe(z) < 0? 4: 0;
2826 63 : return 1;
2827 : }
2828 7 : a = gel(z,1); ta = typ(a);
2829 7 : b = gel(z,2); tb = typ(b);
2830 :
2831 7 : T->t = 0;
2832 7 : if (ta == t_INT && !signe(a))
2833 : {
2834 0 : T->u = R_abs_shallow(b);
2835 0 : T->t = gsigne(b) < 0? -2: 2;
2836 0 : return 1;
2837 : }
2838 7 : if (tb == t_INT && !signe(b))
2839 : {
2840 0 : T->u = R_abs_shallow(a);
2841 0 : T->t = gsigne(a) < 0? 4: 0;
2842 0 : return 1;
2843 : }
2844 7 : if (ta != tb || ta == t_REAL) { T->u = z; return 0; }
2845 : /* a,b both non zero, both t_INT or t_FRAC */
2846 7 : if (ta == t_INT)
2847 : {
2848 7 : if (!absequalii(a, b)) return 0;
2849 7 : T->u = absi_shallow(a);
2850 7 : T->v = 1;
2851 7 : if (signe(a) == signe(b))
2852 0 : { T->t = signe(a) < 0? -3: 1; }
2853 : else
2854 7 : { T->t = signe(a) < 0? 3: -1; }
2855 : }
2856 : else
2857 : {
2858 0 : if (!absequalii(gel(a,2), gel(b,2)) || !absequalii(gel(a,1),gel(b,1)))
2859 0 : return 0;
2860 0 : T->u = absfrac_shallow(a);
2861 0 : T->v = 1;
2862 0 : a = gel(a,1);
2863 0 : b = gel(b,1);
2864 0 : if (signe(a) == signe(b))
2865 0 : { T->t = signe(a) < 0? -3: 1; }
2866 : else
2867 0 : { T->t = signe(a) < 0? 3: -1; }
2868 : }
2869 7 : return 1;
2870 : }
2871 :
2872 : /* z * sqrt(st_b) / sqrt(st_a) exp(I Pi (t + t0)). Assume that
2873 : * sqrt2 = gsqrt(gen_2, prec) or NULL */
2874 : static GEN
2875 35 : apply_eta_correction(GEN z, GEN st_a, GEN st_b, GEN t0, GEN sqrt2, long prec)
2876 : {
2877 35 : GEN t, s_a = gel(st_a, 1), s_b = gel(st_b, 1);
2878 : cxanalyze_t Ta, Tb;
2879 : int ca, cb;
2880 :
2881 35 : t = gsub(gel(st_b,2), gel(st_a,2));
2882 35 : if (t0 != gen_0) t = gadd(t, t0);
2883 35 : ca = cxanalyze(&Ta, s_a);
2884 35 : cb = cxanalyze(&Tb, s_b);
2885 35 : if (ca || cb)
2886 35 : { /* compute sqrt(s_b) / sqrt(s_a) in a more efficient way:
2887 : * sb = ub sqrt(2)^vb exp(i Pi/4 tb) */
2888 35 : GEN u = gdiv(Tb.u,Ta.u);
2889 35 : switch(Tb.v - Ta.v)
2890 : {
2891 0 : case -1: u = gmul2n(u,-1); /* fall through: write 1/sqrt2 = sqrt2/2 */
2892 7 : case 1: u = gmul(u, sqrt2? sqrt2: sqrtr_abs(real2n(1, prec)));
2893 : }
2894 35 : if (!isint1(u)) z = gmul(z, gsqrt(u, prec));
2895 35 : t = gadd(t, gmul2n(stoi(Tb.t - Ta.t), -3));
2896 : }
2897 : else
2898 : {
2899 0 : z = gmul(z, gsqrt(s_b, prec));
2900 0 : z = gdiv(z, gsqrt(s_a, prec));
2901 : }
2902 35 : return gmul(z, expIPiQ(t, prec));
2903 : }
2904 :
2905 : /* sqrt(2) eta(2x) / eta(x) */
2906 : GEN
2907 7 : weberf2(GEN x, long prec)
2908 : {
2909 7 : pari_sp av = avma;
2910 : GEN z, sqrt2, a,b, Ua,Ub, st_a,st_b;
2911 :
2912 7 : x = upper_to_cx(x, &prec);
2913 7 : a = cxredsl2(x, &Ua);
2914 7 : b = cxredsl2(gmul2n(x,1), &Ub);
2915 7 : if (gequal(a,b)) /* not infrequent */
2916 0 : z = gen_1;
2917 : else
2918 7 : z = gdiv(eta_reduced(b,prec), eta_reduced(a,prec));
2919 7 : st_a = eta_correction(a, Ua, 1);
2920 7 : st_b = eta_correction(b, Ub, 1);
2921 7 : sqrt2 = sqrtr_abs(real2n(1, prec));
2922 7 : z = apply_eta_correction(z, st_a, st_b, gen_0, sqrt2, prec);
2923 7 : return gerepileupto(av, gmul(z, sqrt2));
2924 : }
2925 :
2926 : /* eta(x/2) / eta(x) */
2927 : GEN
2928 14 : weberf1(GEN x, long prec)
2929 : {
2930 14 : pari_sp av = avma;
2931 : GEN z, a,b, Ua,Ub, st_a,st_b;
2932 :
2933 14 : x = upper_to_cx(x, &prec);
2934 14 : a = cxredsl2(x, &Ua);
2935 14 : b = cxredsl2(gmul2n(x,-1), &Ub);
2936 14 : if (gequal(a,b)) /* not infrequent */
2937 0 : z = gen_1;
2938 : else
2939 14 : z = gdiv(eta_reduced(b,prec), eta_reduced(a,prec));
2940 14 : st_a = eta_correction(a, Ua, 1);
2941 14 : st_b = eta_correction(b, Ub, 1);
2942 14 : z = apply_eta_correction(z, st_a, st_b, gen_0, NULL, prec);
2943 14 : return gerepileupto(av, z);
2944 : }
2945 : /* exp(-I*Pi/24) * eta((x+1)/2) / eta(x) */
2946 : GEN
2947 14 : weberf(GEN x, long prec)
2948 : {
2949 14 : pari_sp av = avma;
2950 : GEN z, t0, a,b, Ua,Ub, st_a,st_b;
2951 14 : x = upper_to_cx(x, &prec);
2952 14 : a = cxredsl2(x, &Ua);
2953 14 : b = cxredsl2(gmul2n(gaddgs(x,1),-1), &Ub);
2954 14 : if (gequal(a,b)) /* not infrequent */
2955 7 : z = gen_1;
2956 : else
2957 7 : z = gdiv(eta_reduced(b,prec), eta_reduced(a,prec));
2958 14 : st_a = eta_correction(a, Ua, 1);
2959 14 : st_b = eta_correction(b, Ub, 1);
2960 14 : t0 = mkfrac(gen_m1, utoipos(24));
2961 14 : z = apply_eta_correction(z, st_a, st_b, t0, NULL, prec);
2962 14 : if (typ(z) == t_COMPLEX && isexactzero(real_i(x)))
2963 0 : z = gerepilecopy(av, gel(z,1));
2964 : else
2965 14 : z = gerepileupto(av, z);
2966 14 : return z;
2967 : }
2968 : GEN
2969 35 : weber0(GEN x, long flag,long prec)
2970 : {
2971 35 : switch(flag)
2972 : {
2973 14 : case 0: return weberf(x,prec);
2974 14 : case 1: return weberf1(x,prec);
2975 7 : case 2: return weberf2(x,prec);
2976 0 : default: pari_err_FLAG("weber");
2977 : }
2978 : return NULL; /* LCOV_EXCL_LINE */
2979 : }
2980 :
2981 : /* check |q| < 1 */
2982 : static GEN
2983 21 : check_unit_disc(const char *fun, GEN q, long prec)
2984 : {
2985 21 : GEN Q = gtofp(q, prec), Qlow;
2986 21 : Qlow = (prec > LOWDEFAULTPREC)? gtofp(Q,LOWDEFAULTPREC): Q;
2987 21 : if (gcmp(gnorm(Qlow), gen_1) >= 0)
2988 0 : pari_err_DOMAIN(fun, "abs(q)", ">=", gen_1, q);
2989 21 : return Q;
2990 : }
2991 :
2992 : GEN
2993 14 : theta(GEN q, GEN z, long prec)
2994 : {
2995 : long l, n;
2996 14 : pari_sp av = avma, av2;
2997 : GEN s, c, snz, cnz, s2z, c2z, ps, qn, y, zy, ps2, k, zold;
2998 :
2999 14 : l = precision(q);
3000 14 : n = precision(z); if (n && n < l) l = n;
3001 14 : if (l) prec = l;
3002 14 : z = gtofp(z, prec);
3003 14 : q = check_unit_disc("theta", q, prec);
3004 14 : zold = NULL; /* gcc -Wall */
3005 14 : zy = imag_i(z);
3006 14 : if (gequal0(zy)) k = gen_0;
3007 : else
3008 : {
3009 7 : GEN lq = glog(q,prec); k = roundr(divrr(zy, real_i(lq)));
3010 7 : if (signe(k)) { zold = z; z = gadd(z, mulcxmI(gmul(lq,k))); }
3011 : }
3012 14 : qn = gen_1;
3013 14 : ps2 = gsqr(q);
3014 14 : ps = gneg_i(ps2);
3015 14 : gsincos(z, &s, &c, prec);
3016 14 : s2z = gmul2n(gmul(s,c),1); /* sin 2z */
3017 14 : c2z = gsubgs(gmul2n(gsqr(c),1), 1); /* cos 2z */
3018 14 : snz = s;
3019 14 : cnz = c; y = s;
3020 14 : av2 = avma;
3021 14 : for (n = 3;; n += 2)
3022 259 : {
3023 : long e;
3024 273 : s = gadd(gmul(snz, c2z), gmul(cnz,s2z));
3025 273 : qn = gmul(qn,ps);
3026 273 : y = gadd(y, gmul(qn, s));
3027 273 : e = gexpo(s); if (e < 0) e = 0;
3028 273 : if (gexpo(qn) + e < -prec2nbits(prec)) break;
3029 :
3030 259 : ps = gmul(ps,ps2);
3031 259 : c = gsub(gmul(cnz, c2z), gmul(snz,s2z));
3032 259 : snz = s; /* sin nz */
3033 259 : cnz = c; /* cos nz */
3034 259 : if (gc_needed(av,2))
3035 : {
3036 0 : if (DEBUGMEM>1) pari_warn(warnmem,"theta (n = %ld)", n);
3037 0 : gerepileall(av2, 5, &snz, &cnz, &ps, &qn, &y);
3038 : }
3039 : }
3040 14 : if (signe(k))
3041 : {
3042 7 : y = gmul(y, gmul(powgi(q,sqri(k)),
3043 : gexp(gmul(mulcxI(zold),shifti(k,1)), prec)));
3044 7 : if (mod2(k)) y = gneg_i(y);
3045 : }
3046 14 : return gerepileupto(av, gmul(y, gmul2n(gsqrt(gsqrt(q,prec),prec),1)));
3047 : }
3048 :
3049 : GEN
3050 7 : thetanullk(GEN q, long k, long prec)
3051 : {
3052 : long l, n;
3053 7 : pari_sp av = avma;
3054 : GEN p1, ps, qn, y, ps2;
3055 :
3056 7 : if (k < 0)
3057 0 : pari_err_DOMAIN("thetanullk", "k", "<", gen_0, stoi(k));
3058 7 : l = precision(q);
3059 7 : if (l) prec = l;
3060 7 : q = check_unit_disc("thetanullk", q, prec);
3061 :
3062 7 : if (!odd(k)) { set_avma(av); return gen_0; }
3063 7 : qn = gen_1;
3064 7 : ps2 = gsqr(q);
3065 7 : ps = gneg_i(ps2);
3066 7 : y = gen_1;
3067 7 : for (n = 3;; n += 2)
3068 280 : {
3069 : GEN t;
3070 287 : qn = gmul(qn,ps);
3071 287 : ps = gmul(ps,ps2);
3072 287 : t = gmul(qn, powuu(n, k)); y = gadd(y, t);
3073 287 : if (gexpo(t) < -prec2nbits(prec)) break;
3074 : }
3075 7 : p1 = gmul2n(gsqrt(gsqrt(q,prec),prec),1);
3076 7 : if (k&2) y = gneg_i(y);
3077 7 : return gerepileupto(av, gmul(p1, y));
3078 : }
3079 :
3080 : /* q2 = q^2 */
3081 : static GEN
3082 13472 : vecthetanullk_loop(GEN q2, long k, long prec)
3083 : {
3084 13472 : GEN ps, qn = gen_1, y = const_vec(k, gen_1);
3085 13472 : pari_sp av = avma;
3086 13472 : const long bit = prec2nbits(prec);
3087 : long i, n;
3088 :
3089 13472 : if (gexpo(q2) < -2*bit) return y;
3090 13472 : ps = gneg_i(q2);
3091 13472 : for (n = 3;; n += 2)
3092 95595 : {
3093 109067 : GEN t = NULL/*-Wall*/, P = utoipos(n), N2 = sqru(n);
3094 109067 : qn = gmul(qn,ps);
3095 109067 : ps = gmul(ps,q2);
3096 327201 : for (i = 1; i <= k; i++)
3097 : {
3098 218134 : t = gmul(qn, P); gel(y,i) = gadd(gel(y,i), t);
3099 218134 : P = mulii(P, N2);
3100 : }
3101 109067 : if (gexpo(t) < -bit) return y;
3102 95595 : if (gc_needed(av,2))
3103 : {
3104 0 : if (DEBUGMEM>1) pari_warn(warnmem,"vecthetanullk_loop, n = %ld",n);
3105 0 : gerepileall(av, 3, &qn, &ps, &y);
3106 : }
3107 : }
3108 : }
3109 : /* [d^i theta/dz^i(q, 0), i = 1, 3, .., 2*k - 1] */
3110 : GEN
3111 0 : vecthetanullk(GEN q, long k, long prec)
3112 : {
3113 0 : long i, l = precision(q);
3114 0 : pari_sp av = avma;
3115 : GEN p1, y;
3116 :
3117 0 : if (l) prec = l;
3118 0 : q = check_unit_disc("vecthetanullk", q, prec);
3119 0 : y = vecthetanullk_loop(gsqr(q), k, prec);
3120 0 : p1 = gmul2n(gsqrt(gsqrt(q,prec),prec),1);
3121 0 : for (i = 2; i <= k; i += 2) gel(y,i) = gneg_i(gel(y,i));
3122 0 : return gerepileupto(av, gmul(p1, y));
3123 : }
3124 :
3125 : /* [d^i theta/dz^i(q, 0), i = 1, 3, .., 2*k - 1], q = exp(2iPi tau) */
3126 : GEN
3127 0 : vecthetanullk_tau(GEN tau, long k, long prec)
3128 : {
3129 0 : long i, l = precision(tau);
3130 0 : pari_sp av = avma;
3131 : GEN q4, y;
3132 :
3133 0 : if (l) prec = l;
3134 0 : if (typ(tau) != t_COMPLEX || gsigne(gel(tau,2)) <= 0)
3135 0 : pari_err_DOMAIN("vecthetanullk_tau", "imag(tau)", "<=", gen_0, tau);
3136 0 : q4 = expIPiC(gmul2n(tau,-1), prec); /* q^(1/4) */
3137 0 : y = vecthetanullk_loop(gpowgs(q4,8), k, prec);
3138 0 : for (i = 2; i <= k; i += 2) gel(y,i) = gneg_i(gel(y,i));
3139 0 : return gerepileupto(av, gmul(gmul2n(q4,1), y));
3140 : }
3141 :
3142 : /* Return E_k(tau). Slow if tau is not in standard fundamental domain */
3143 : GEN
3144 13741 : cxEk(GEN tau, long k, long prec)
3145 : {
3146 : pari_sp av;
3147 : GEN q, y, qn;
3148 13741 : long n, b, l = precision(tau);
3149 :
3150 13741 : if (l) prec = l;
3151 13741 : b = bit_accuracy(prec);
3152 : /* sum n^(k-1) x^n <= x(1 + (k!-1)x) / (1-x)^k (cf Eulerian polynomials)
3153 : * S = \sum_{n > 0} n^(k-1) |q^n/(1-q^n)| <= x(1+(k!-1)x) / (1-x)^(k+1),
3154 : * where x = |q| = exp(-2Pi Im(tau)) < 1. Neglegt 2/zeta(1-k) * S if
3155 : * (2Pi)^k/(k-1)! x < 2^(-b-1) and k! x < 1. Use log2((2Pi)^k/(k-1)!) < 10 */
3156 13741 : if (gcmpgs(imag_i(tau), (M_LN2 / (2*M_PI)) * (b+1+10)) > 0)
3157 17 : return real_1(prec);
3158 13724 : if (k == 2)
3159 : { /* -theta^(3)(tau/2) / theta^(1)(tau/2). Assume that Im tau > 0 */
3160 13472 : y = vecthetanullk_loop(qq(tau,prec), 2, prec);
3161 13472 : return gdiv(gel(y,2), gel(y,1));
3162 : }
3163 252 : q = cxtoreal(expIPiC(gneg(gmul2n(tau,1)), prec));
3164 252 : av = avma; y = gen_0; qn = q;
3165 252 : for(n = 1;; n++)
3166 2605 : { /* compute y := sum_{n>0} n^(k-1) / (q^n-1) */
3167 2857 : GEN t = gdiv(powuu(n,k-1), gsubgs(qn,1));
3168 2857 : if (gequal0(t) || gexpo(t) <= - prec2nbits(prec) - 5) break;
3169 2605 : y = gadd(y, t);
3170 2605 : qn = gmul(q,qn);
3171 2605 : if (gc_needed(av,2))
3172 : {
3173 0 : if(DEBUGMEM>1) pari_warn(warnmem,"elleisnum");
3174 0 : gerepileall(av, 2, &y,&qn);
3175 : }
3176 : }
3177 252 : return gadd(gen_1, gmul(y, gdiv(gen_2, szeta(1-k, prec))));
3178 : }
|
