Line data Source code
1 : /* Copyright (C) 2022 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : #define DEBUGLEVEL DEBUGLEVEL_trans
19 :
20 : /********************************************************/
21 : /* Hurwitz zeta function */
22 : /********************************************************/
23 : struct Qp_zetahurwitz_t { GEN B, _1, s1; };
24 : static void
25 301 : Qp_zetahurwitz_init(struct Qp_zetahurwitz_t *S, long prec, GEN s)
26 : {
27 301 : GEN B, C = gen_1, s1 = gsubgs(s, 1), p = padic_p(s);
28 301 : long j, J = ((equaliu(p,2)? (prec >> 1): prec) + 2) >> 1;
29 301 : if (gequal0(s1)) s1 = NULL;
30 301 : B = bernvec(J);
31 3346 : for (j = 1; j <= J; j++)
32 : {
33 3045 : GEN t = (j == 1 && !s1)? s: gmul(gaddgs(s, 2*j-3), gaddgs(s, 2*j-2));
34 3045 : C = gdivgunextu(gmul(C, t), 2*j-1);
35 3045 : gel(B, j+1) = gmul(gel(B, j+1), C); /* B_{2j} * binomial(1-s, 2j) */
36 : }
37 301 : S->_1 = cvtop(gen_1, p, prec);
38 301 : S->s1 = s1;
39 301 : S->B = B;
40 301 : }
41 :
42 : /* v_p(x) < (p==2)?-1: 0; s1 = s-1 or NULL (if s=1) */
43 : static GEN
44 1218 : Qp_zetahurwitz_0(struct Qp_zetahurwitz_t *S, GEN x)
45 : {
46 1218 : GEN z, x2, x2j, s1 = S->s1;
47 1218 : long j, J = lg(S->B) - 2;
48 :
49 1218 : x = cvtop2(ginv(x), S->_1); z = gmul2n(x, -1);
50 1218 : z = s1? gmul(s1, z): gadd(Qp_log(x), z);
51 1218 : x2j = x2 = gsqr(x); z = gaddgs(z, 1);
52 1218 : for (j = 1;; j++)
53 : {
54 21525 : z = gadd(z, gmul(gel(S->B, j + 1), x2j));
55 21525 : if (j == J) break;
56 20307 : x2j = gmul(x2, x2j);
57 : }
58 1218 : if (s1) z = gmul(gdiv(z, s1), Qp_exp(gmul(s1, Qp_log(x))));
59 1218 : return z;
60 : }
61 : /* private (absolute) padicprec */
62 : static long
63 546 : pprec(GEN x) { return maxss(valp(x) + precp(x), 1); }
64 :
65 : /* L_p(s, (D, .)); assume s != 1 if D = 1 */
66 : static GEN
67 14 : Qp_zeta_i(GEN s, long D)
68 : {
69 14 : pari_sp av = avma;
70 14 : GEN z, va, gp = padic_p(s);
71 14 : ulong a, p = itou(gp), m;
72 14 : long prec = pprec(s);
73 : struct Qp_zetahurwitz_t S;
74 :
75 14 : if (D <= 0) pari_err_DOMAIN("p-adic L-function", "D", "<=", gen_0, stoi(D));
76 14 : if (!uposisfundamental(D))
77 0 : pari_err_TYPE("p-adic L-function [D not fundamental]", stoi(D));
78 14 : Qp_zetahurwitz_init(&S, prec, s);
79 14 : m = ulcm(D, p == 2? 4: p); va = coprimes_zv(m);
80 42 : for (a = 1, z = gen_0; a <= (m >> 1); a++)
81 28 : if (va[a])
82 : {
83 21 : GEN h = Qp_zetahurwitz_0(&S, uutoQ(a, m));
84 21 : if (D != 1 && kross(D, a) < 0) h = gneg(h);
85 21 : z = gadd(z, h);
86 : }
87 14 : z = gdivgs(gmul2n(z, 1), m);
88 14 : if (D != 1) z = gmul(z, Qp_exp(gmul(gsubsg(1, s), Qp_log(cvstop2(m, s)))));
89 14 : return gc_upto(av, z);
90 : }
91 : GEN
92 14 : Qp_zeta(GEN s) { return Qp_zeta_i(s, 1); }
93 :
94 : /* s a t_PADIC; gc_upto-safe. Could be exported */
95 : static GEN
96 287 : Qp_zetahurwitz_ii(GEN s, GEN x, long k)
97 : {
98 287 : GEN gp = padic_p(s);
99 287 : long p = itou(gp), prec = pprec(s);
100 : struct Qp_zetahurwitz_t S;
101 287 : Qp_zetahurwitz_init(&S, prec, s);
102 287 : if (typ(x) != t_PADIC) x = gadd(x, zeropadic_shallow(gp, prec));
103 287 : if (valp(x) >= ((p==2)? -1: 0))
104 : {
105 280 : GEN z = gen_0;
106 280 : long j, M = (p==2)? 4: p;
107 1869 : for (j = 0; j < M; j++)
108 : {
109 1589 : GEN y = gaddsg(j, x);
110 1589 : if (valp(y) <= 0)
111 : {
112 1190 : GEN tmp = Qp_zetahurwitz_0(&S, gdivgu(y, M));
113 1190 : if (k) tmp = gmul(tmp, gpowgs(teich(y), k));
114 1190 : z = gadd(z, tmp);
115 : }
116 : }
117 280 : return gdivgu(z, M);
118 : }
119 7 : if (valp(s) < 0) pari_err_DOMAIN("Qp_zetahurwitz", "v(s)", "<", gen_0, s);
120 7 : return Qp_zetahurwitz_0(&S, x);
121 : }
122 :
123 : /* x or s must be p-adic */
124 : static GEN
125 287 : Qp_zetahurwitz_i(GEN s, GEN x, long k)
126 : {
127 287 : if (typ(x) == t_PADIC)
128 : {
129 245 : pari_sp av = avma;
130 245 : GEN p = padic_p(x);
131 245 : long e = pprec(x);
132 245 : e += sdivsi(e, subis(p, 1));
133 245 : s = gadd(s, zeropadic_shallow(p, e));
134 245 : return gc_upto(av, Qp_zetahurwitz_ii(s, x, k));
135 : }
136 42 : return Qp_zetahurwitz_ii(s, x, k);
137 : }
138 :
139 : GEN
140 238 : Qp_zetahurwitz(GEN s, GEN x, long k)
141 : {
142 238 : pari_sp av = avma;
143 238 : return gc_upto(av, Qp_zetahurwitz_i(s, x, k));
144 : }
145 :
146 : static void
147 8596 : binsplit(GEN *pP, GEN *pR, GEN aN2, GEN isqaN, GEN s, long j, long k, long prec)
148 : {
149 8596 : if (j + 1 == k)
150 : {
151 4347 : long j2 = j << 1;
152 : GEN P;
153 4347 : if (!j) P = gdiv(s, aN2);
154 : else
155 : {
156 4249 : P = gmul(gaddgs(s, j2-1), gaddgs(s, j2));
157 4249 : P = gdivgunextu(gmul(P, isqaN), j2+1);
158 : }
159 4347 : if (pP) *pP = P;
160 4347 : if (pR) *pR = gmul(bernreal(j2+2, prec), P);
161 : }
162 : else
163 : {
164 : GEN P1, R1, P2, R2;
165 4249 : binsplit(&P1,pR? &R1: NULL, aN2, isqaN, s, j, (j+k) >> 1, prec);
166 4249 : binsplit(pP? &P2: NULL, pR? &R2: NULL, aN2, isqaN, s, (j+k) >> 1, k, prec);
167 4249 : if (pP) *pP = gmul(P1,P2);
168 4249 : if (pR) *pR = gadd(R1, gmul(P1, R2));
169 : }
170 8596 : }
171 :
172 : /* a0 + a1 x + O(x^e), e >= 0 */
173 : static GEN
174 77 : deg1ser_shallow(GEN a1, GEN a0, long v, long e)
175 77 : { return RgX_to_ser(deg1pol_shallow(a1, a0, v), e+2); }
176 :
177 : static long
178 539 : hurwitz_cutoff(GEN s, long bit)
179 : {
180 588 : return typ(s) == t_COMPLEX &&
181 49 : fabs(gtodouble(gel(s,2))) > 5.37 * pow(bit, 1.4) / mt_nbthreads();
182 : }
183 :
184 : /* New zetahurwitz, from Fredrik Johansson. */
185 : GEN
186 672 : zetahurwitz(GEN s, GEN x, long der, long bitprec)
187 : {
188 672 : pari_sp av = avma, av2;
189 672 : GEN a, ra, ra0, Nx, S1, S2, S3, N2, rx, sch = NULL, s0 = s, x0 = x, y;
190 672 : long j, k, m, N, prec0 = nbits2prec(bitprec), prec = prec0, fli = 0;
191 : pari_timer T;
192 :
193 672 : if (typ(s) == t_PADIC || typ(x) == t_PADIC)
194 49 : return gc_upto(av, Qp_zetahurwitz_i(s, x, der));
195 623 : if (der < 0) pari_err_DOMAIN("zetahurwitz", "der", "<", gen_0, stoi(der));
196 623 : if (der)
197 : {
198 : GEN z;
199 21 : if (!is_scalar_t(typ(s)))
200 : {
201 7 : z = deriv(zetahurwitz(s, x, der - 1, bitprec), -1);
202 7 : z = gdiv(z, deriv(s, -1));
203 : }
204 : else
205 : {
206 14 : if (gequal1(s)) pari_err_DOMAIN("zetahurwitz", "s", "=", gen_1, s0);
207 14 : s = deg1ser_shallow(gen_1, s, 0, der+2);
208 14 : z = zetahurwitz(s, x, 0, bitprec + der * log2(der));
209 14 : z = gmul(mpfact(der), polcoef_i(z, der, -1));
210 : }
211 21 : return gc_upto(av,z);
212 : }
213 602 : switch(typ(x))
214 : {
215 441 : case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: break;
216 161 : default:
217 161 : if (!(y = toser_i(x))) pari_err_TYPE("zetahurwitz", x);
218 154 : x = y; x0 = polcoef_i(x, 0, -1); break;
219 : }
220 595 : rx = real_i(x0);
221 595 : if (typ(x) != t_SER && typ(rx) == t_INT && signe(rx) <= 0
222 84 : && gequal0(imag_i(x0)))
223 0 : pari_err_DOMAIN("zetahurwitz","x", "=",
224 : strtoGENstr("nonpositive integer"), x0);
225 595 : rx = grndtoi(rx, NULL);
226 595 : if (typ(rx) != t_INT) pari_err_TYPE("zetahurwitz", x);
227 595 : switch (typ(s))
228 : {
229 : long v, pr;
230 525 : case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX:
231 525 : if (!der && hurwitz_cutoff(s, bitprec))
232 7 : return zetahurwitzlarge(s, x, prec);
233 518 : break;
234 70 : default:
235 70 : if (!(y = toser_i(s))) pari_err_TYPE("zetahurwitz", s);
236 70 : if (valser(y) < 0) pari_err_DOMAIN("zetahurwitz", "val(s)", "<", gen_0, s);
237 70 : s0 = polcoef_i(y, 0, -1);
238 70 : switch(typ(s0))
239 : {
240 63 : case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: break;
241 0 : case t_PADIC: pari_err_IMPL("zetahurwitz(t_SER of t_PADIC)");
242 7 : default: pari_err_TYPE("zetahurwitz", s0);
243 : }
244 63 : sch = gequal0(s0)? y: serchop0(y);
245 63 : v = valser(sch);
246 63 : pr = (lg(y) + v + 1) / v;
247 63 : if (gequal1(s0)) pr += v;
248 63 : s = deg1ser_shallow(gen_1, s0, 0, pr);
249 : }
250 581 : a = gneg(s0); ra = real_i(a); ra0 = ground(ra);
251 581 : if (gequal1(s0) && (!sch || gequal0(sch)))
252 14 : pari_err_DOMAIN("zetahurwitz", "s", "=", gen_1, s0);
253 567 : fli = (gsigne(ra0) >= 0 && gexpo(gsub(a, ra0)) < 17 - bitprec);
254 567 : if (!sch && fli)
255 : { /* a ~ non negative integer */
256 14 : k = itos(gceil(ra)) + 1;
257 14 : if (odd(k)) k++;
258 14 : N = 1;
259 : }
260 : else
261 : {
262 553 : GEN C, ix = imag_i(x0);
263 553 : double c = (typ(s) == t_INT)? 1: 20 * log((double)bitprec);
264 553 : double rs = gtodouble(ra) + 1;
265 : long k0, bit;
266 553 : if (fli) a = gadd(a, ghalf); /* hack */
267 553 : if (rs > 0)
268 : {
269 49 : bitprec += (long)ceil(rs * expu(bitprec));
270 49 : prec = nbits2prec(bitprec);
271 49 : x = gprec_w(x, prec);
272 49 : s = gprec_w(s, prec);
273 49 : if (sch) sch = gprec_w(sch, prec);
274 : }
275 553 : k = bitprec * M_LN2 / (1 + dbllambertW0(M_PI / c));
276 553 : k0 = itos(gceil(gadd(ra, ghalf))) + 1;
277 553 : k = maxss(k0, k);
278 553 : if (odd(k)) k++;
279 : /* R_k < 2 |binom(a,k+1) B_{k+2}/(k+2)| * |N + x - 1|^(ra - k - 1)*/
280 553 : C = binomial(a, k+1); C = polcoef_i(C, 0, -1);
281 553 : C = gmul(C, gdivgu(bernfrac(k+2), k+2));
282 553 : bit = bitprec;
283 : /* if a < 0 and |x| >> 1, |zeta(s,x)| ~ |x|^ra is small: compensate for
284 : * correct relative accuracy */
285 553 : if (rs < 0)
286 : {
287 497 : double dx = dbllog2(x0);
288 497 : if (dx > 0) bit -= dx * gtodouble(ra);
289 : }
290 : /* + 1 i from the factor 2 in RHS above */
291 553 : C = gmul2n(gabs(C,LOWDEFAULTPREC), bit + 1);
292 553 : C = gpow(C, ginv(gsubsg(k+1, ra)), LOWDEFAULTPREC);
293 : /* need |N + x - 1|^2 > C^2 to have R_k < 2^(-bit) */
294 553 : if (!gequal0(ix))
295 : {
296 133 : GEN tmp = gsub(gsqr(C), gsqr(ix));
297 133 : C = (gsigne(tmp) <= 0)? NULL: gsqrt(tmp, LOWDEFAULTPREC);
298 : }
299 553 : N = 1;
300 553 : if (C)
301 : { /* now need |N + re(x) - 1| > C */
302 490 : C = gadd(C, gsubsg(1, rx));
303 490 : if (!is_real_t(typ(C))) pari_err_TYPE("zetahurwitz",s);
304 490 : if (gsigne(C) > 0) N = itos(gceil(C));
305 : }
306 553 : if (N == 1 && signe(a) > 0)
307 : { /* May reduce k if 2Pix > a */
308 : /* Need 2 |x^(-K) (B_K/K) binom(a, K-1)| < 2^-bit |x|^-rs |zeta(s,x)|
309 : * with K = k+2; N = 1; |zeta(s,x)| ~ |x|^(rs-1);
310 : * if a > 0, (B_K/K) binom(a, K-1) < 2 |a / 2Pi|^K */
311 0 : double dx = dbllog2(x0), d = 1 + dx + log2(M_PI) - dbllog2(s0);
312 0 : if (d > 0)
313 : { /* d ~ log2 |2Pi x / a| */
314 0 : long K = (long)ceil((bitprec + 1 + dx) / d);
315 0 : K = maxss(k0, K);
316 0 : if (odd(K)) K++;
317 0 : if (K < k) k = K;
318 : }
319 : }
320 : }
321 567 : if (gsigne(rx) < 0) N = maxss(N, 1 - itos(rx));
322 567 : a = gneg(s);
323 567 : if (DEBUGLEVEL>2) timer_start(&T);
324 567 : incrprec(prec);
325 567 : Nx = gaddsg(N - 1, x);
326 567 : Nx = typ(Nx) == t_SER? RgX_gtofp(Nx, prec): gtofp(Nx, prec);
327 567 : S1 = S3 = gpow(Nx, a, prec);
328 567 : av2 = avma;
329 567 : if (gequal1(x)) S1 = dirpowerssum(N, a, 0, prec);
330 : else
331 6790 : for (m = N - 2; m >= 0; m--)
332 : {
333 6279 : S1 = gadd(S1, gpow(gaddsg(m,x), a, prec));
334 6279 : if ((m & 0xff) == 0) S1 = gc_upto(av2, S1);
335 : }
336 567 : if (DEBUGLEVEL>2) timer_printf(&T,"sum from 0 to N - 1");
337 567 : constbern(k >> 1);
338 567 : N2 = ginv(gsqr(Nx));
339 567 : if (typ(s0) == t_INT)
340 : {
341 469 : S2 = divru(bernreal(k, prec), k);
342 10920 : for (j = k - 2; j >= 2; j -= 2)
343 : {
344 10451 : GEN t = gsubgs(a, j), u = gmul(t, gaddgs(t, 1));
345 10451 : u = gmul(gdivgunextu(u, j), gmul(S2, N2));
346 10451 : S2 = gadd(divru(bernreal(j, prec), j), u);
347 : }
348 469 : S2 = gmul(S2, gdiv(a, Nx));
349 : }
350 : else
351 : {
352 98 : binsplit(NULL,&S2, gmul2n(Nx,1), N2, s, 0, k >> 1, prec);
353 98 : S2 = gneg(S2);
354 : }
355 567 : S2 = gadd(ghalf, S2);
356 567 : if (DEBUGLEVEL>2) timer_printf(&T,"Bernoulli sum");
357 567 : S2 = gmul(S3, gadd(gdiv(Nx, gaddsg(1, a)), S2));
358 567 : S1 = gprec_wtrunc(gsub(S1, S2), prec0);
359 567 : if (sch) return gc_upto(av, gsubst(S1, 0, sch));
360 511 : return gc_GEN(av, S1);
361 : }
362 :
363 : /* New Lerch, inspired by Fredrik Johansson. */
364 :
365 : GEN
366 153919 : lerch_worker(GEN t, GEN E)
367 : {
368 153919 : GEN n, d, T, s = gel(E,1), a = gmul(gel(E,2), t), z = gel(E,3);
369 153919 : long p = itos(gel(E,4)), prec = labs(p);
370 153919 : d = gadd(gexp(t, prec), z);
371 153919 : T = p > 0? t: gneg(t);
372 153919 : if (typ(s) == t_INT)
373 58576 : n = gmul(gpow(T, s, prec), gexp(a, prec));
374 : else /* save one exp */
375 95343 : n = gexp(gadd(gmul(s, glog(T, prec)), a), prec);
376 153919 : return gdiv(n, d);
377 : }
378 :
379 : /* tab already computed with N = #tab[1] even */
380 : static GEN
381 441 : parintnumgauss(GEN f, GEN a, GEN b, GEN tab, long prec)
382 : {
383 441 : GEN R = gel(tab, 1), W = gel(tab, 2), bma, bpa, S = gen_0, VP, VM, V;
384 441 : long n = lg(R) - 1, i, prec2 = prec + EXTRAPREC64;
385 441 : a = gprec_wensure(a, prec2);
386 441 : b = gprec_wensure(b, prec2);
387 441 : VP = cgetg(n + 1, t_VEC); bma = gmul2n(gsub(b, a), -1);
388 441 : VM = cgetg(n + 1, t_VEC); bpa = gadd(bma, a);
389 28203 : for (i = 1; i <= n; i++)
390 : {
391 27762 : GEN h = gmul(bma, gel(R, i));
392 27762 : gel(VP, i) = gadd(bpa, h);
393 27762 : gel(VM, i) = gsub(bpa, h);
394 : }
395 441 : V = gadd(parapply(f, VP), parapply(f, VM));
396 28203 : for (i = 1; i <= n; i++)
397 : {
398 27762 : S = gadd(S, gmul(gel(W, i), gel(V, i)));
399 27762 : S = gprec_wensure(S, prec2);
400 : }
401 441 : return gprec_wtrunc(gmul(bma, S), prec);
402 : }
403 :
404 : /* Assume tab computed and a >= 0 */
405 : static GEN
406 119 : parintnum(GEN f, GEN a, GEN tab)
407 : {
408 : pari_sp av;
409 119 : GEN tabx0 = gel(tab, 2), tabw0 = gel(tab, 3), tabxm = gel(tab, 6);
410 119 : GEN tabxp = gel(tab, 4), tabwp = gel(tab, 5), tabwm = gel(tab, 7);
411 119 : GEN VP = tabxp, VM = tabxm, x0 = tabx0, S;
412 119 : long prec = gprecision(tabw0), L = lg(tabxp), i, fla = 0;
413 119 : if (!gequal0(a))
414 : {
415 91 : if (gexpo(a) <= 0)
416 : {
417 63 : x0 = gadd(a, x0);
418 30357 : for (i = 1; i < L; i++)
419 : {
420 30294 : gel(VP, i) = gadd(a, gel(VP, i));
421 30294 : gel(VM, i) = gadd(a, gel(VM, i));
422 : }
423 : }
424 : else
425 : {
426 28 : x0 = gmul(a, gaddsg(1, x0)); fla = 1;
427 5404 : for (i = 1; i < L; i++)
428 : {
429 5376 : gel(VP, i) = gmul(a, gaddsg(1, gel(VP, i)));
430 5376 : gel(VM, i) = gmul(a, gaddsg(1, gel(VM, i)));
431 : }
432 : }
433 : }
434 119 : VP = parapply(f, VP);
435 119 : VM = parapply(f, VM); av = avma;
436 119 : S = gmul(tabw0, closure_callgen1(f, x0));
437 49257 : for (i = 1; i < L; i++)
438 : {
439 49138 : S = gadd(S, gadd(gmul(gel(tabwp, i), gel(VP, i)),
440 49138 : gmul(gel(tabwm, i), gel(VM, i))));
441 49138 : if ((i & 0x7f) == 1) S = gc_upto(av, S);
442 49138 : S = gprec_wensure(S, prec);
443 : }
444 119 : if (fla) S = gmul(S, a);
445 119 : return gmul(S, gel(tab, 1));
446 : }
447 :
448 : static GEN
449 84 : refine(GEN A)
450 : {
451 84 : long n = lg(A) - 1, i;
452 84 : GEN B = cgetg(2 * n, t_VEC);
453 231 : for (i = 1; i < n; i++)
454 : {
455 147 : gel(B, 2 * i - 1) = gel(A, i);
456 147 : gel(B, 2 * i) = gmul2n(gadd(gel(A, i), gel(A, i + 1)), -1);
457 : }
458 84 : gel(B, 2 * n - 1) = gel(A, n); return B;
459 : }
460 :
461 : /* Here L = [a1, a2, a3,...] integration vertices. Refine by splitting
462 : * intervals. */
463 : static GEN
464 84 : parintnumgaussadapt(GEN f, GEN L, GEN tab, long bit)
465 : {
466 84 : GEN Rold = gen_0, Rnew;
467 84 : long i, ct = 0, prec = nbits2prec(bit);
468 168 : while (ct <= 5)
469 : {
470 168 : Rnew = gen_0;
471 609 : for (i = 1; i < lg(L) - 1; i++)
472 441 : Rnew = gadd(Rnew, parintnumgauss(f, gel(L, i), gel(L, i + 1), tab, prec));
473 168 : if (ct && gexpo(gsub(Rnew, Rold)) - gexpo(Rnew) < 10 - bit) return Rnew;
474 84 : ct++; Rold = Rnew; L = refine(L);
475 : }
476 0 : if (DEBUGLEVEL) err_printf("intnumgaussadapt: possible accuracy loss");
477 0 : return Rnew; /* Possible accuracy loss */
478 : }
479 :
480 : /* Here b = [oo, r], so refine by increasing integration step m */
481 : static GEN
482 42 : parintnumadapt(GEN f, GEN a, GEN b, GEN tab, long bit)
483 : {
484 42 : GEN Rold = gen_0, Rnew;
485 42 : long m = 0, prec = nbits2prec(bit);
486 42 : if (!tab) tab = intnuminit(gen_0, b, 0, prec);
487 119 : while (m <= 5)
488 : {
489 119 : Rnew = parintnum(f, a, tab);
490 119 : if (m && gexpo(gsub(Rnew, Rold)) - gexpo(Rnew) < 10 - bit) return Rnew;
491 77 : m++; Rold = Rnew; tab = intnuminit(gen_0, b, m, prec);
492 : }
493 0 : if (DEBUGLEVEL) err_printf("intnumadapt: possible accuracy loss");
494 0 : return Rnew; /* Possible accuracy loss */
495 : }
496 :
497 : static int
498 210 : iscplx(GEN z) { long t = typ(z); return is_real_t(t) || t == t_COMPLEX; }
499 :
500 : static GEN
501 14 : lerch_easy(GEN z, GEN s, GEN a, long B)
502 : {
503 14 : long n, prec = nbits2prec(B + 32);
504 14 : GEN zn, ms = gneg(s), S = gpow(a, ms, prec);
505 14 : zn = z = gtofp(z, prec);
506 3808 : for (n = 1;; n++, zn = gmul(zn, z))
507 : {
508 3808 : S = gadd(S, gmul(zn, gpow(gaddgs(a, n), ms, prec)));
509 3808 : if (gexpo(zn) <= - B - 5) return S;
510 : }
511 : }
512 :
513 : static GEN
514 112 : _lerchphi(GEN z, GEN s, GEN a, long prec)
515 : {
516 112 : GEN res = NULL, L, LT, J, rs, mleft, left, right, top, w, Linf, tabg;
517 : GEN E, f, fm;
518 112 : long B = prec2nbits(prec), MB = 3 - B, NB, prec2;
519 : entree *ep;
520 :
521 112 : if (gexpo(z) < MB) return gpow(a, gneg(s), prec);
522 112 : if (gexpo(gsubgs(z, 1)) < MB) return zetahurwitz(s, a, 0, B); /* z ~ 1 */
523 112 : if (gexpo(gaddgs(z, 1)) < MB) /* z ~ -1 */
524 : {
525 7 : GEN tmp = gsub(zetahurwitz(s, gmul2n(a, -1), 0, B),
526 : zetahurwitz(s, gmul2n(gaddgs(a, 1), -1), 0, B));
527 7 : return gmul(gpow(gen_2, gneg(s), prec), tmp);
528 : }
529 105 : if (gcmpgs(gmulsg(10, gabs(z, prec)), 9) <= 0) /* |z| <= 9/10 */
530 14 : return lerch_easy(z, s, a, B);
531 91 : if (gcmpgs(real_i(a), 2) < 0)
532 49 : return gadd(gpow(a, gneg(s), prec),
533 : gmul(z, _lerchphi(z, s, gaddgs(a, 1), prec)));
534 42 : NB = (long)ceil(B + M_PI * fabs(gtodouble(imag_i(s))));
535 42 : prec2 = nbits2prec(NB);
536 42 : z = gprec_w(z, prec2); /* |z| > 9/10 */
537 42 : s = gprec_w(s, prec2);
538 42 : a = gprec_w(a, prec2); /* Re(a) >= 2 */
539 42 : rs = ground(real_i(s)); L = glog(z, prec2); /* Re(L) > -0.11 */
540 42 : ep = is_entry("_lerch_worker");
541 42 : E = mkvec4(gsubgs(s, 1), gsubsg(1, a), gneg(z), stoi(prec2));
542 42 : f = snm_closure(ep, mkvec(E));
543 42 : E = shallowcopy(E); gel(E,4) = stoi(-prec2);
544 42 : fm = snm_closure(ep, mkvec(E));
545 42 : Linf = mkvec2(mkoo(), real_i(a));
546 42 : if (gexpo(gsub(s, rs)) < MB && gcmpgs(rs, 1) >= 0)
547 : { /* s ~ positive integer */
548 14 : if (gcmp(gabs(imag_i(L), prec2), sstoQ(1, 4)) < 0 && gsigne(real_i(L)) >= 0)
549 7 : { /* Re(L) >= 0, |Im(L)| < 1/4 */
550 7 : GEN t = gsigne(imag_i(z)) > 0 ? gen_m1: gen_1;
551 7 : GEN LT1 = gaddgs(gabs(L, prec2), 1);
552 7 : LT = mkvec4(gen_0, mkcomplex(gen_0, t), mkcomplex(LT1, t), LT1);
553 7 : tabg = intnumgaussinit(2*(NB >> 2) + 60, prec2);
554 7 : J = parintnumgaussadapt(f, LT, tabg, NB);
555 7 : J = gadd(J, parintnumadapt(f, LT1, Linf, NULL, NB));
556 : }
557 7 : else J = parintnumadapt(f, gen_0, Linf, NULL, NB);
558 14 : return gdiv(J, ggamma(s, prec2));
559 : }
560 28 : tabg = intnumgaussinit(2*(NB >> 2) + 60, prec2);
561 28 : if (gcmp(gabs(imag_i(L), prec2), ghalf) > 0) /* |Im(L)| > 1/2 */
562 14 : left = right = top = gmin(gmul2n(gabs(imag_i(L), prec2), -1), gen_1);
563 : else
564 : {
565 14 : res = gdiv(gpow(gneg(L), s, prec2), gmul(L, gpow(z, a, prec2)));
566 14 : left = gaddgs(gmax(gen_0, gneg(real_i(L))), 1);
567 14 : top = gaddgs(gabs(imag_i(L), prec2), 1);
568 14 : right = gaddgs(gabs(L, prec2), 1);
569 : }
570 28 : w = expIPiC(gsubgs(s, 1), prec2);
571 28 : mleft = gneg(left);
572 28 : if (gexpo(imag_i(z)) < MB && gexpo(imag_i(a)) < MB && gexpo(imag_i(s)) < MB
573 7 : && gcmpgs(real_i(z), 1) < 0)
574 : { /* (z, s, a) real, z < 1 */
575 7 : LT = mkvec3(right, mkcomplex(right, top), mkcomplex(mleft, top));
576 7 : J = imag_i(gdiv(parintnumgaussadapt(f, LT, tabg, NB), w));
577 7 : LT = mkvec2(mkcomplex(mleft, top), mleft);
578 7 : J = gmul2n(gadd(J, imag_i(parintnumgaussadapt(fm, LT, tabg, NB))), 1);
579 7 : J = mulcxI(J);
580 : }
581 : else
582 : {
583 21 : GEN mtop = gneg(top);
584 21 : LT = mkvec3(right, mkcomplex(right, top), mkcomplex(mleft, top));
585 21 : J = gdiv(parintnumgaussadapt(f, LT, tabg, NB), w);
586 21 : LT = mkvec2(mkcomplex(mleft, top), mkcomplex(mleft, mtop));
587 21 : J = gadd(J, parintnumgaussadapt(fm, LT, tabg, NB));
588 21 : LT = mkvec3(mkcomplex(mleft, mtop), mkcomplex(right, mtop), right);
589 21 : J = gadd(J, gmul(parintnumgaussadapt(f, LT, tabg, NB), w));
590 : }
591 28 : J = gadd(J, gmul(gsub(w, ginv(w)), parintnumadapt(f, right, Linf, NULL, NB)));
592 28 : J = gdiv(J, PiI2(prec2)); if (res) J = gadd(J, res);
593 28 : return gneg(gmul(ggamma(gsubsg(1, s), prec2), J));
594 : }
595 : /* lerchphi(z,-k,a)=
596 : * -1/(z-1)*sum(q=0,k,(z/(z-1))^q*sum(j=0,q,(-1)^j*(j+a)^k*binomial(q,j)))
597 : * zetahurwitz(-k,a)=-B(k+1,a)/(k+1) */
598 : GEN
599 56 : lerchphi(GEN z, GEN s, GEN a, long prec)
600 : {
601 56 : pari_sp av = avma;
602 56 : if (!iscplx(z)) pari_err_TYPE("lerchphi", z);
603 56 : if (!iscplx(s)) pari_err_TYPE("lerchphi", s);
604 56 : if (!iscplx(a)) pari_err_TYPE("lerchphi", a);
605 56 : return gc_upto(av, _lerchphi(z, s, a, prec));
606 : }
607 :
608 : GEN
609 14 : lerchzeta(GEN s, GEN a, GEN lam, long prec)
610 : {
611 14 : pari_sp av = avma;
612 14 : GEN z = gexp(gmul(PiI2(prec), lam), prec);
613 14 : if (!iscplx(z)) pari_err_TYPE("lerchzeta", z);
614 14 : if (!iscplx(s)) pari_err_TYPE("lerchzeta", s);
615 14 : if (!iscplx(a)) pari_err_TYPE("lerchzeta", a);
616 14 : if (hurwitz_cutoff(s, prec)) return lerchzetalarge(s, a, lam, prec);
617 7 : return gc_upto(av, _lerchphi(z, s, a, prec));
618 : }
|
