Coverage for /usr/lib/python3/dist-packages/mpmath/functions/zeta.py: 7%
829 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from __future__ import print_function
3from ..libmp.backend import xrange
4from .functions import defun, defun_wrapped, defun_static
6@defun
7def stieltjes(ctx, n, a=1):
8 n = ctx.convert(n)
9 a = ctx.convert(a)
10 if n < 0:
11 return ctx.bad_domain("Stieltjes constants defined for n >= 0")
12 if hasattr(ctx, "stieltjes_cache"):
13 stieltjes_cache = ctx.stieltjes_cache
14 else:
15 stieltjes_cache = ctx.stieltjes_cache = {}
16 if a == 1:
17 if n == 0:
18 return +ctx.euler
19 if n in stieltjes_cache:
20 prec, s = stieltjes_cache[n]
21 if prec >= ctx.prec:
22 return +s
23 mag = 1
24 def f(x):
25 xa = x/a
26 v = (xa-ctx.j)*ctx.ln(a-ctx.j*x)**n/(1+xa**2)/(ctx.exp(2*ctx.pi*x)-1)
27 return ctx._re(v) / mag
28 orig = ctx.prec
29 try:
30 # Normalize integrand by approx. magnitude to
31 # speed up quadrature (which uses absolute error)
32 if n > 50:
33 ctx.prec = 20
34 mag = ctx.quad(f, [0,ctx.inf], maxdegree=3)
35 ctx.prec = orig + 10 + int(n**0.5)
36 s = ctx.quad(f, [0,ctx.inf], maxdegree=20)
37 v = ctx.ln(a)**n/(2*a) - ctx.ln(a)**(n+1)/(n+1) + 2*s/a*mag
38 finally:
39 ctx.prec = orig
40 if a == 1 and ctx.isint(n):
41 stieltjes_cache[n] = (ctx.prec, v)
42 return +v
44@defun_wrapped
45def siegeltheta(ctx, t, derivative=0):
46 d = int(derivative)
47 if (t == ctx.inf or t == ctx.ninf):
48 if d < 2:
49 if t == ctx.ninf and d == 0:
50 return ctx.ninf
51 return ctx.inf
52 else:
53 return ctx.zero
54 if d == 0:
55 if ctx._im(t):
56 # XXX: cancellation occurs
57 a = ctx.loggamma(0.25+0.5j*t)
58 b = ctx.loggamma(0.25-0.5j*t)
59 return -ctx.ln(ctx.pi)/2*t - 0.5j*(a-b)
60 else:
61 if ctx.isinf(t):
62 return t
63 return ctx._im(ctx.loggamma(0.25+0.5j*t)) - ctx.ln(ctx.pi)/2*t
64 if d > 0:
65 a = (-0.5j)**(d-1)*ctx.polygamma(d-1, 0.25-0.5j*t)
66 b = (0.5j)**(d-1)*ctx.polygamma(d-1, 0.25+0.5j*t)
67 if ctx._im(t):
68 if d == 1:
69 return -0.5*ctx.log(ctx.pi)+0.25*(a+b)
70 else:
71 return 0.25*(a+b)
72 else:
73 if d == 1:
74 return ctx._re(-0.5*ctx.log(ctx.pi)+0.25*(a+b))
75 else:
76 return ctx._re(0.25*(a+b))
78@defun_wrapped
79def grampoint(ctx, n):
80 # asymptotic expansion, from
81 # http://mathworld.wolfram.com/GramPoint.html
82 g = 2*ctx.pi*ctx.exp(1+ctx.lambertw((8*n+1)/(8*ctx.e)))
83 return ctx.findroot(lambda t: ctx.siegeltheta(t)-ctx.pi*n, g)
86@defun_wrapped
87def siegelz(ctx, t, **kwargs):
88 d = int(kwargs.get("derivative", 0))
89 t = ctx.convert(t)
90 t1 = ctx._re(t)
91 t2 = ctx._im(t)
92 prec = ctx.prec
93 try:
94 if abs(t1) > 500*prec and t2**2 < t1:
95 v = ctx.rs_z(t, d)
96 if ctx._is_real_type(t):
97 return ctx._re(v)
98 return v
99 except NotImplementedError:
100 pass
101 ctx.prec += 21
102 e1 = ctx.expj(ctx.siegeltheta(t))
103 z = ctx.zeta(0.5+ctx.j*t)
104 if d == 0:
105 v = e1*z
106 ctx.prec=prec
107 if ctx._is_real_type(t):
108 return ctx._re(v)
109 return +v
110 z1 = ctx.zeta(0.5+ctx.j*t, derivative=1)
111 theta1 = ctx.siegeltheta(t, derivative=1)
112 if d == 1:
113 v = ctx.j*e1*(z1+z*theta1)
114 ctx.prec=prec
115 if ctx._is_real_type(t):
116 return ctx._re(v)
117 return +v
118 z2 = ctx.zeta(0.5+ctx.j*t, derivative=2)
119 theta2 = ctx.siegeltheta(t, derivative=2)
120 comb1 = theta1**2-ctx.j*theta2
121 if d == 2:
122 def terms():
123 return [2*z1*theta1, z2, z*comb1]
124 v = ctx.sum_accurately(terms, 1)
125 v = -e1*v
126 ctx.prec = prec
127 if ctx._is_real_type(t):
128 return ctx._re(v)
129 return +v
130 ctx.prec += 10
131 z3 = ctx.zeta(0.5+ctx.j*t, derivative=3)
132 theta3 = ctx.siegeltheta(t, derivative=3)
133 comb2 = theta1**3-3*ctx.j*theta1*theta2-theta3
134 if d == 3:
135 def terms():
136 return [3*theta1*z2, 3*z1*comb1, z3+z*comb2]
137 v = ctx.sum_accurately(terms, 1)
138 v = -ctx.j*e1*v
139 ctx.prec = prec
140 if ctx._is_real_type(t):
141 return ctx._re(v)
142 return +v
143 z4 = ctx.zeta(0.5+ctx.j*t, derivative=4)
144 theta4 = ctx.siegeltheta(t, derivative=4)
145 def terms():
146 return [theta1**4, -6*ctx.j*theta1**2*theta2, -3*theta2**2,
147 -4*theta1*theta3, ctx.j*theta4]
148 comb3 = ctx.sum_accurately(terms, 1)
149 if d == 4:
150 def terms():
151 return [6*theta1**2*z2, -6*ctx.j*z2*theta2, 4*theta1*z3,
152 4*z1*comb2, z4, z*comb3]
153 v = ctx.sum_accurately(terms, 1)
154 v = e1*v
155 ctx.prec = prec
156 if ctx._is_real_type(t):
157 return ctx._re(v)
158 return +v
159 if d > 4:
160 h = lambda x: ctx.siegelz(x, derivative=4)
161 return ctx.diff(h, t, n=d-4)
164_zeta_zeros = [
16514.134725142,21.022039639,25.010857580,30.424876126,32.935061588,
16637.586178159,40.918719012,43.327073281,48.005150881,49.773832478,
16752.970321478,56.446247697,59.347044003,60.831778525,65.112544048,
16867.079810529,69.546401711,72.067157674,75.704690699,77.144840069,
16979.337375020,82.910380854,84.735492981,87.425274613,88.809111208,
17092.491899271,94.651344041,95.870634228,98.831194218,101.317851006,
171103.725538040,105.446623052,107.168611184,111.029535543,111.874659177,
172114.320220915,116.226680321,118.790782866,121.370125002,122.946829294,
173124.256818554,127.516683880,129.578704200,131.087688531,133.497737203,
174134.756509753,138.116042055,139.736208952,141.123707404,143.111845808,
175146.000982487,147.422765343,150.053520421,150.925257612,153.024693811,
176156.112909294,157.597591818,158.849988171,161.188964138,163.030709687,
177165.537069188,167.184439978,169.094515416,169.911976479,173.411536520,
178174.754191523,176.441434298,178.377407776,179.916484020,182.207078484,
179184.874467848,185.598783678,187.228922584,189.416158656,192.026656361,
180193.079726604,195.265396680,196.876481841,198.015309676,201.264751944,
181202.493594514,204.189671803,205.394697202,207.906258888,209.576509717,
182211.690862595,213.347919360,214.547044783,216.169538508,219.067596349,
183220.714918839,221.430705555,224.007000255,224.983324670,227.421444280,
184229.337413306,231.250188700,231.987235253,233.693404179,236.524229666,
185]
187def _load_zeta_zeros(url):
188 import urllib
189 d = urllib.urlopen(url)
190 L = [float(x) for x in d.readlines()]
191 # Sanity check
192 assert round(L[0]) == 14
193 _zeta_zeros[:] = L
195@defun
196def oldzetazero(ctx, n, url='http://www.dtc.umn.edu/~odlyzko/zeta_tables/zeros1'):
197 n = int(n)
198 if n < 0:
199 return ctx.zetazero(-n).conjugate()
200 if n == 0:
201 raise ValueError("n must be nonzero")
202 if n > len(_zeta_zeros) and n <= 100000:
203 _load_zeta_zeros(url)
204 if n > len(_zeta_zeros):
205 raise NotImplementedError("n too large for zetazeros")
206 return ctx.mpc(0.5, ctx.findroot(ctx.siegelz, _zeta_zeros[n-1]))
208@defun_wrapped
209def riemannr(ctx, x):
210 if x == 0:
211 return ctx.zero
212 # Check if a simple asymptotic estimate is accurate enough
213 if abs(x) > 1000:
214 a = ctx.li(x)
215 b = 0.5*ctx.li(ctx.sqrt(x))
216 if abs(b) < abs(a)*ctx.eps:
217 return a
218 if abs(x) < 0.01:
219 # XXX
220 ctx.prec += int(-ctx.log(abs(x),2))
221 # Sum Gram's series
222 s = t = ctx.one
223 u = ctx.ln(x)
224 k = 1
225 while abs(t) > abs(s)*ctx.eps:
226 t = t * u / k
227 s += t / (k * ctx._zeta_int(k+1))
228 k += 1
229 return s
231@defun_static
232def primepi(ctx, x):
233 x = int(x)
234 if x < 2:
235 return 0
236 return len(ctx.list_primes(x))
238# TODO: fix the interface wrt contexts
239@defun_wrapped
240def primepi2(ctx, x):
241 x = int(x)
242 if x < 2:
243 return ctx._iv.zero
244 if x < 2657:
245 return ctx._iv.mpf(ctx.primepi(x))
246 mid = ctx.li(x)
247 # Schoenfeld's estimate for x >= 2657, assuming RH
248 err = ctx.sqrt(x,rounding='u')*ctx.ln(x,rounding='u')/8/ctx.pi(rounding='d')
249 a = ctx.floor((ctx._iv.mpf(mid)-err).a, rounding='d')
250 b = ctx.ceil((ctx._iv.mpf(mid)+err).b, rounding='u')
251 return ctx._iv.mpf([a,b])
253@defun_wrapped
254def primezeta(ctx, s):
255 if ctx.isnan(s):
256 return s
257 if ctx.re(s) <= 0:
258 raise ValueError("prime zeta function defined only for re(s) > 0")
259 if s == 1:
260 return ctx.inf
261 if s == 0.5:
262 return ctx.mpc(ctx.ninf, ctx.pi)
263 r = ctx.re(s)
264 if r > ctx.prec:
265 return 0.5**s
266 else:
267 wp = ctx.prec + int(r)
268 def terms():
269 orig = ctx.prec
270 # zeta ~ 1+eps; need to set precision
271 # to get logarithm accurately
272 k = 0
273 while 1:
274 k += 1
275 u = ctx.moebius(k)
276 if not u:
277 continue
278 ctx.prec = wp
279 t = u*ctx.ln(ctx.zeta(k*s))/k
280 if not t:
281 return
282 #print ctx.prec, ctx.nstr(t)
283 ctx.prec = orig
284 yield t
285 return ctx.sum_accurately(terms)
287# TODO: for bernpoly and eulerpoly, ensure that all exact zeros are covered
289@defun_wrapped
290def bernpoly(ctx, n, z):
291 # Slow implementation:
292 #return sum(ctx.binomial(n,k)*ctx.bernoulli(k)*z**(n-k) for k in xrange(0,n+1))
293 n = int(n)
294 if n < 0:
295 raise ValueError("Bernoulli polynomials only defined for n >= 0")
296 if z == 0 or (z == 1 and n > 1):
297 return ctx.bernoulli(n)
298 if z == 0.5:
299 return (ctx.ldexp(1,1-n)-1)*ctx.bernoulli(n)
300 if n <= 3:
301 if n == 0: return z ** 0
302 if n == 1: return z - 0.5
303 if n == 2: return (6*z*(z-1)+1)/6
304 if n == 3: return z*(z*(z-1.5)+0.5)
305 if ctx.isinf(z):
306 return z ** n
307 if ctx.isnan(z):
308 return z
309 if abs(z) > 2:
310 def terms():
311 t = ctx.one
312 yield t
313 r = ctx.one/z
314 k = 1
315 while k <= n:
316 t = t*(n+1-k)/k*r
317 if not (k > 2 and k & 1):
318 yield t*ctx.bernoulli(k)
319 k += 1
320 return ctx.sum_accurately(terms) * z**n
321 else:
322 def terms():
323 yield ctx.bernoulli(n)
324 t = ctx.one
325 k = 1
326 while k <= n:
327 t = t*(n+1-k)/k * z
328 m = n-k
329 if not (m > 2 and m & 1):
330 yield t*ctx.bernoulli(m)
331 k += 1
332 return ctx.sum_accurately(terms)
334@defun_wrapped
335def eulerpoly(ctx, n, z):
336 n = int(n)
337 if n < 0:
338 raise ValueError("Euler polynomials only defined for n >= 0")
339 if n <= 2:
340 if n == 0: return z ** 0
341 if n == 1: return z - 0.5
342 if n == 2: return z*(z-1)
343 if ctx.isinf(z):
344 return z**n
345 if ctx.isnan(z):
346 return z
347 m = n+1
348 if z == 0:
349 return -2*(ctx.ldexp(1,m)-1)*ctx.bernoulli(m)/m * z**0
350 if z == 1:
351 return 2*(ctx.ldexp(1,m)-1)*ctx.bernoulli(m)/m * z**0
352 if z == 0.5:
353 if n % 2:
354 return ctx.zero
355 # Use exact code for Euler numbers
356 if n < 100 or n*ctx.mag(0.46839865*n) < ctx.prec*0.25:
357 return ctx.ldexp(ctx._eulernum(n), -n)
358 # http://functions.wolfram.com/Polynomials/EulerE2/06/01/02/01/0002/
359 def terms():
360 t = ctx.one
361 k = 0
362 w = ctx.ldexp(1,n+2)
363 while 1:
364 v = n-k+1
365 if not (v > 2 and v & 1):
366 yield (2-w)*ctx.bernoulli(v)*t
367 k += 1
368 if k > n:
369 break
370 t = t*z*(n-k+2)/k
371 w *= 0.5
372 return ctx.sum_accurately(terms) / m
374@defun
375def eulernum(ctx, n, exact=False):
376 n = int(n)
377 if exact:
378 return int(ctx._eulernum(n))
379 if n < 100:
380 return ctx.mpf(ctx._eulernum(n))
381 if n % 2:
382 return ctx.zero
383 return ctx.ldexp(ctx.eulerpoly(n,0.5), n)
385# TODO: this should be implemented low-level
386def polylog_series(ctx, s, z):
387 tol = +ctx.eps
388 l = ctx.zero
389 k = 1
390 zk = z
391 while 1:
392 term = zk / k**s
393 l += term
394 if abs(term) < tol:
395 break
396 zk *= z
397 k += 1
398 return l
400def polylog_continuation(ctx, n, z):
401 if n < 0:
402 return z*0
403 twopij = 2j * ctx.pi
404 a = -twopij**n/ctx.fac(n) * ctx.bernpoly(n, ctx.ln(z)/twopij)
405 if ctx._is_real_type(z) and z < 0:
406 a = ctx._re(a)
407 if ctx._im(z) < 0 or (ctx._im(z) == 0 and ctx._re(z) >= 1):
408 a -= twopij*ctx.ln(z)**(n-1)/ctx.fac(n-1)
409 return a
411def polylog_unitcircle(ctx, n, z):
412 tol = +ctx.eps
413 if n > 1:
414 l = ctx.zero
415 logz = ctx.ln(z)
416 logmz = ctx.one
417 m = 0
418 while 1:
419 if (n-m) != 1:
420 term = ctx.zeta(n-m) * logmz / ctx.fac(m)
421 if term and abs(term) < tol:
422 break
423 l += term
424 logmz *= logz
425 m += 1
426 l += ctx.ln(z)**(n-1)/ctx.fac(n-1)*(ctx.harmonic(n-1)-ctx.ln(-ctx.ln(z)))
427 elif n < 1: # else
428 l = ctx.fac(-n)*(-ctx.ln(z))**(n-1)
429 logz = ctx.ln(z)
430 logkz = ctx.one
431 k = 0
432 while 1:
433 b = ctx.bernoulli(k-n+1)
434 if b:
435 term = b*logkz/(ctx.fac(k)*(k-n+1))
436 if abs(term) < tol:
437 break
438 l -= term
439 logkz *= logz
440 k += 1
441 else:
442 raise ValueError
443 if ctx._is_real_type(z) and z < 0:
444 l = ctx._re(l)
445 return l
447def polylog_general(ctx, s, z):
448 v = ctx.zero
449 u = ctx.ln(z)
450 if not abs(u) < 5: # theoretically |u| < 2*pi
451 j = ctx.j
452 v = 1-s
453 y = ctx.ln(-z)/(2*ctx.pi*j)
454 return ctx.gamma(v)*(j**v*ctx.zeta(v,0.5+y) + j**-v*ctx.zeta(v,0.5-y))/(2*ctx.pi)**v
455 t = 1
456 k = 0
457 while 1:
458 term = ctx.zeta(s-k) * t
459 if abs(term) < ctx.eps:
460 break
461 v += term
462 k += 1
463 t *= u
464 t /= k
465 return ctx.gamma(1-s)*(-u)**(s-1) + v
467@defun_wrapped
468def polylog(ctx, s, z):
469 s = ctx.convert(s)
470 z = ctx.convert(z)
471 if z == 1:
472 return ctx.zeta(s)
473 if z == -1:
474 return -ctx.altzeta(s)
475 if s == 0:
476 return z/(1-z)
477 if s == 1:
478 return -ctx.ln(1-z)
479 if s == -1:
480 return z/(1-z)**2
481 if abs(z) <= 0.75 or (not ctx.isint(s) and abs(z) < 0.9):
482 return polylog_series(ctx, s, z)
483 if abs(z) >= 1.4 and ctx.isint(s):
484 return (-1)**(s+1)*polylog_series(ctx, s, 1/z) + polylog_continuation(ctx, s, z)
485 if ctx.isint(s):
486 return polylog_unitcircle(ctx, int(s), z)
487 return polylog_general(ctx, s, z)
489@defun_wrapped
490def clsin(ctx, s, z, pi=False):
491 if ctx.isint(s) and s < 0 and int(s) % 2 == 1:
492 return z*0
493 if pi:
494 a = ctx.expjpi(z)
495 else:
496 a = ctx.expj(z)
497 if ctx._is_real_type(z) and ctx._is_real_type(s):
498 return ctx.im(ctx.polylog(s,a))
499 b = 1/a
500 return (-0.5j)*(ctx.polylog(s,a) - ctx.polylog(s,b))
502@defun_wrapped
503def clcos(ctx, s, z, pi=False):
504 if ctx.isint(s) and s < 0 and int(s) % 2 == 0:
505 return z*0
506 if pi:
507 a = ctx.expjpi(z)
508 else:
509 a = ctx.expj(z)
510 if ctx._is_real_type(z) and ctx._is_real_type(s):
511 return ctx.re(ctx.polylog(s,a))
512 b = 1/a
513 return 0.5*(ctx.polylog(s,a) + ctx.polylog(s,b))
515@defun
516def altzeta(ctx, s, **kwargs):
517 try:
518 return ctx._altzeta(s, **kwargs)
519 except NotImplementedError:
520 return ctx._altzeta_generic(s)
522@defun_wrapped
523def _altzeta_generic(ctx, s):
524 if s == 1:
525 return ctx.ln2 + 0*s
526 return -ctx.powm1(2, 1-s) * ctx.zeta(s)
528@defun
529def zeta(ctx, s, a=1, derivative=0, method=None, **kwargs):
530 d = int(derivative)
531 if a == 1 and not (d or method):
532 try:
533 return ctx._zeta(s, **kwargs)
534 except NotImplementedError:
535 pass
536 s = ctx.convert(s)
537 prec = ctx.prec
538 method = kwargs.get('method')
539 verbose = kwargs.get('verbose')
540 if (not s) and (not derivative):
541 return ctx.mpf(0.5) - ctx._convert_param(a)[0]
542 if a == 1 and method != 'euler-maclaurin':
543 im = abs(ctx._im(s))
544 re = abs(ctx._re(s))
545 #if (im < prec or method == 'borwein') and not derivative:
546 # try:
547 # if verbose:
548 # print "zeta: Attempting to use the Borwein algorithm"
549 # return ctx._zeta(s, **kwargs)
550 # except NotImplementedError:
551 # if verbose:
552 # print "zeta: Could not use the Borwein algorithm"
553 # pass
554 if abs(im) > 500*prec and 10*re < prec and derivative <= 4 or \
555 method == 'riemann-siegel':
556 try: # py2.4 compatible try block
557 try:
558 if verbose:
559 print("zeta: Attempting to use the Riemann-Siegel algorithm")
560 return ctx.rs_zeta(s, derivative, **kwargs)
561 except NotImplementedError:
562 if verbose:
563 print("zeta: Could not use the Riemann-Siegel algorithm")
564 pass
565 finally:
566 ctx.prec = prec
567 if s == 1:
568 return ctx.inf
569 abss = abs(s)
570 if abss == ctx.inf:
571 if ctx.re(s) == ctx.inf:
572 if d == 0:
573 return ctx.one
574 return ctx.zero
575 return s*0
576 elif ctx.isnan(abss):
577 return 1/s
578 if ctx.re(s) > 2*ctx.prec and a == 1 and not derivative:
579 return ctx.one + ctx.power(2, -s)
580 return +ctx._hurwitz(s, a, d, **kwargs)
582@defun
583def _hurwitz(ctx, s, a=1, d=0, **kwargs):
584 prec = ctx.prec
585 verbose = kwargs.get('verbose')
586 try:
587 extraprec = 10
588 ctx.prec += extraprec
589 # We strongly want to special-case rational a
590 a, atype = ctx._convert_param(a)
591 if ctx.re(s) < 0:
592 if verbose:
593 print("zeta: Attempting reflection formula")
594 try:
595 return _hurwitz_reflection(ctx, s, a, d, atype)
596 except NotImplementedError:
597 pass
598 if verbose:
599 print("zeta: Reflection formula failed")
600 if verbose:
601 print("zeta: Using the Euler-Maclaurin algorithm")
602 while 1:
603 ctx.prec = prec + extraprec
604 T1, T2 = _hurwitz_em(ctx, s, a, d, prec+10, verbose)
605 cancellation = ctx.mag(T1) - ctx.mag(T1+T2)
606 if verbose:
607 print("Term 1:", T1)
608 print("Term 2:", T2)
609 print("Cancellation:", cancellation, "bits")
610 if cancellation < extraprec:
611 return T1 + T2
612 else:
613 extraprec = max(2*extraprec, min(cancellation + 5, 100*prec))
614 if extraprec > kwargs.get('maxprec', 100*prec):
615 raise ctx.NoConvergence("zeta: too much cancellation")
616 finally:
617 ctx.prec = prec
619def _hurwitz_reflection(ctx, s, a, d, atype):
620 # TODO: implement for derivatives
621 if d != 0:
622 raise NotImplementedError
623 res = ctx.re(s)
624 negs = -s
625 # Integer reflection formula
626 if ctx.isnpint(s):
627 n = int(res)
628 if n <= 0:
629 return ctx.bernpoly(1-n, a) / (n-1)
630 if not (atype == 'Q' or atype == 'Z'):
631 raise NotImplementedError
632 t = 1-s
633 # We now require a to be standardized
634 v = 0
635 shift = 0
636 b = a
637 while ctx.re(b) > 1:
638 b -= 1
639 v -= b**negs
640 shift -= 1
641 while ctx.re(b) <= 0:
642 v += b**negs
643 b += 1
644 shift += 1
645 # Rational reflection formula
646 try:
647 p, q = a._mpq_
648 except:
649 assert a == int(a)
650 p = int(a)
651 q = 1
652 p += shift*q
653 assert 1 <= p <= q
654 g = ctx.fsum(ctx.cospi(t/2-2*k*b)*ctx._hurwitz(t,(k,q)) \
655 for k in range(1,q+1))
656 g *= 2*ctx.gamma(t)/(2*ctx.pi*q)**t
657 v += g
658 return v
660def _hurwitz_em(ctx, s, a, d, prec, verbose):
661 # May not be converted at this point
662 a = ctx.convert(a)
663 tol = -prec
664 # Estimate number of terms for Euler-Maclaurin summation; could be improved
665 M1 = 0
666 M2 = prec // 3
667 N = M2
668 lsum = 0
669 # This speeds up the recurrence for derivatives
670 if ctx.isint(s):
671 s = int(ctx._re(s))
672 s1 = s-1
673 while 1:
674 # Truncated L-series
675 l = ctx._zetasum(s, M1+a, M2-M1-1, [d])[0][0]
676 #if d:
677 # l = ctx.fsum((-ctx.ln(n+a))**d * (n+a)**negs for n in range(M1,M2))
678 #else:
679 # l = ctx.fsum((n+a)**negs for n in range(M1,M2))
680 lsum += l
681 M2a = M2+a
682 logM2a = ctx.ln(M2a)
683 logM2ad = logM2a**d
684 logs = [logM2ad]
685 logr = 1/logM2a
686 rM2a = 1/M2a
687 M2as = M2a**(-s)
688 if d:
689 tailsum = ctx.gammainc(d+1, s1*logM2a) / s1**(d+1)
690 else:
691 tailsum = 1/((s1)*(M2a)**s1)
692 tailsum += 0.5 * logM2ad * M2as
693 U = [1]
694 r = M2as
695 fact = 2
696 for j in range(1, N+1):
697 # TODO: the following could perhaps be tidied a bit
698 j2 = 2*j
699 if j == 1:
700 upds = [1]
701 else:
702 upds = [j2-2, j2-1]
703 for m in upds:
704 D = min(m,d+1)
705 if m <= d:
706 logs.append(logs[-1] * logr)
707 Un = [0]*(D+1)
708 for i in xrange(D): Un[i] = (1-m-s)*U[i]
709 for i in xrange(1,D+1): Un[i] += (d-(i-1))*U[i-1]
710 U = Un
711 r *= rM2a
712 t = ctx.fdot(U, logs) * r * ctx.bernoulli(j2)/(-fact)
713 tailsum += t
714 if ctx.mag(t) < tol:
715 return lsum, (-1)**d * tailsum
716 fact *= (j2+1)*(j2+2)
717 if verbose:
718 print("Sum range:", M1, M2, "term magnitude", ctx.mag(t), "tolerance", tol)
719 M1, M2 = M2, M2*2
720 if ctx.re(s) < 0:
721 N += N//2
725@defun
726def _zetasum(ctx, s, a, n, derivatives=[0], reflect=False):
727 """
728 Returns [xd0,xd1,...,xdr], [yd0,yd1,...ydr] where
730 xdk = D^k ( 1/a^s + 1/(a+1)^s + ... + 1/(a+n)^s )
731 ydk = D^k conj( 1/a^(1-s) + 1/(a+1)^(1-s) + ... + 1/(a+n)^(1-s) )
733 D^k = kth derivative with respect to s, k ranges over the given list of
734 derivatives (which should consist of either a single element
735 or a range 0,1,...r). If reflect=False, the ydks are not computed.
736 """
737 #print "zetasum", s, a, n
738 # don't use the fixed-point code if there are large exponentials
739 if abs(ctx.re(s)) < 0.5 * ctx.prec:
740 try:
741 return ctx._zetasum_fast(s, a, n, derivatives, reflect)
742 except NotImplementedError:
743 pass
744 negs = ctx.fneg(s, exact=True)
745 have_derivatives = derivatives != [0]
746 have_one_derivative = len(derivatives) == 1
747 if not reflect:
748 if not have_derivatives:
749 return [ctx.fsum((a+k)**negs for k in xrange(n+1))], []
750 if have_one_derivative:
751 d = derivatives[0]
752 x = ctx.fsum(ctx.ln(a+k)**d * (a+k)**negs for k in xrange(n+1))
753 return [(-1)**d * x], []
754 maxd = max(derivatives)
755 if not have_one_derivative:
756 derivatives = range(maxd+1)
757 xs = [ctx.zero for d in derivatives]
758 if reflect:
759 ys = [ctx.zero for d in derivatives]
760 else:
761 ys = []
762 for k in xrange(n+1):
763 w = a + k
764 xterm = w ** negs
765 if reflect:
766 yterm = ctx.conj(ctx.one / (w * xterm))
767 if have_derivatives:
768 logw = -ctx.ln(w)
769 if have_one_derivative:
770 logw = logw ** maxd
771 xs[0] += xterm * logw
772 if reflect:
773 ys[0] += yterm * logw
774 else:
775 t = ctx.one
776 for d in derivatives:
777 xs[d] += xterm * t
778 if reflect:
779 ys[d] += yterm * t
780 t *= logw
781 else:
782 xs[0] += xterm
783 if reflect:
784 ys[0] += yterm
785 return xs, ys
787@defun
788def dirichlet(ctx, s, chi=[1], derivative=0):
789 s = ctx.convert(s)
790 q = len(chi)
791 d = int(derivative)
792 if d > 2:
793 raise NotImplementedError("arbitrary order derivatives")
794 prec = ctx.prec
795 try:
796 ctx.prec += 10
797 if s == 1:
798 have_pole = True
799 for x in chi:
800 if x and x != 1:
801 have_pole = False
802 h = +ctx.eps
803 ctx.prec *= 2*(d+1)
804 s += h
805 if have_pole:
806 return +ctx.inf
807 z = ctx.zero
808 for p in range(1,q+1):
809 if chi[p%q]:
810 if d == 1:
811 z += chi[p%q] * (ctx.zeta(s, (p,q), 1) - \
812 ctx.zeta(s, (p,q))*ctx.log(q))
813 else:
814 z += chi[p%q] * ctx.zeta(s, (p,q))
815 z /= q**s
816 finally:
817 ctx.prec = prec
818 return +z
821def secondzeta_main_term(ctx, s, a, **kwargs):
822 tol = ctx.eps
823 f = lambda n: ctx.gammainc(0.5*s, a*gamm**2, regularized=True)*gamm**(-s)
824 totsum = term = ctx.zero
825 mg = ctx.inf
826 n = 0
827 while mg > tol:
828 totsum += term
829 n += 1
830 gamm = ctx.im(ctx.zetazero_memoized(n))
831 term = f(n)
832 mg = abs(term)
833 err = 0
834 if kwargs.get("error"):
835 sg = ctx.re(s)
836 err = 0.5*ctx.pi**(-1)*max(1,sg)*a**(sg-0.5)*ctx.log(gamm/(2*ctx.pi))*\
837 ctx.gammainc(-0.5, a*gamm**2)/abs(ctx.gamma(s/2))
838 err = abs(err)
839 return +totsum, err, n
841def secondzeta_prime_term(ctx, s, a, **kwargs):
842 tol = ctx.eps
843 f = lambda n: ctx.gammainc(0.5*(1-s),0.25*ctx.log(n)**2 * a**(-1))*\
844 ((0.5*ctx.log(n))**(s-1))*ctx.mangoldt(n)/ctx.sqrt(n)/\
845 (2*ctx.gamma(0.5*s)*ctx.sqrt(ctx.pi))
846 totsum = term = ctx.zero
847 mg = ctx.inf
848 n = 1
849 while mg > tol or n < 9:
850 totsum += term
851 n += 1
852 term = f(n)
853 if term == 0:
854 mg = ctx.inf
855 else:
856 mg = abs(term)
857 if kwargs.get("error"):
858 err = mg
859 return +totsum, err, n
861def secondzeta_exp_term(ctx, s, a):
862 if ctx.isint(s) and ctx.re(s) <= 0:
863 m = int(round(ctx.re(s)))
864 if not m & 1:
865 return ctx.mpf('-0.25')**(-m//2)
866 tol = ctx.eps
867 f = lambda n: (0.25*a)**n/((n+0.5*s)*ctx.fac(n))
868 totsum = ctx.zero
869 term = f(0)
870 mg = ctx.inf
871 n = 0
872 while mg > tol:
873 totsum += term
874 n += 1
875 term = f(n)
876 mg = abs(term)
877 v = a**(0.5*s)*totsum/ctx.gamma(0.5*s)
878 return v
880def secondzeta_singular_term(ctx, s, a, **kwargs):
881 factor = a**(0.5*(s-1))/(4*ctx.sqrt(ctx.pi)*ctx.gamma(0.5*s))
882 extraprec = ctx.mag(factor)
883 ctx.prec += extraprec
884 factor = a**(0.5*(s-1))/(4*ctx.sqrt(ctx.pi)*ctx.gamma(0.5*s))
885 tol = ctx.eps
886 f = lambda n: ctx.bernpoly(n,0.75)*(4*ctx.sqrt(a))**n*\
887 ctx.gamma(0.5*n)/((s+n-1)*ctx.fac(n))
888 totsum = ctx.zero
889 mg1 = ctx.inf
890 n = 1
891 term = f(n)
892 mg2 = abs(term)
893 while mg2 > tol and mg2 <= mg1:
894 totsum += term
895 n += 1
896 term = f(n)
897 totsum += term
898 n +=1
899 term = f(n)
900 mg1 = mg2
901 mg2 = abs(term)
902 totsum += term
903 pole = -2*(s-1)**(-2)+(ctx.euler+ctx.log(16*ctx.pi**2*a))*(s-1)**(-1)
904 st = factor*(pole+totsum)
905 err = 0
906 if kwargs.get("error"):
907 if not ((mg2 > tol) and (mg2 <= mg1)):
908 if mg2 <= tol:
909 err = ctx.mpf(10)**int(ctx.log(abs(factor*tol),10))
910 if mg2 > mg1:
911 err = ctx.mpf(10)**int(ctx.log(abs(factor*mg1),10))
912 err = max(err, ctx.eps*1.)
913 ctx.prec -= extraprec
914 return +st, err
916@defun
917def secondzeta(ctx, s, a = 0.015, **kwargs):
918 r"""
919 Evaluates the secondary zeta function `Z(s)`, defined for
920 `\mathrm{Re}(s)>1` by
922 .. math ::
924 Z(s) = \sum_{n=1}^{\infty} \frac{1}{\tau_n^s}
926 where `\frac12+i\tau_n` runs through the zeros of `\zeta(s)` with
927 imaginary part positive.
929 `Z(s)` extends to a meromorphic function on `\mathbb{C}` with a
930 double pole at `s=1` and simple poles at the points `-2n` for
931 `n=0`, 1, 2, ...
933 **Examples**
935 >>> from mpmath import *
936 >>> mp.pretty = True; mp.dps = 15
937 >>> secondzeta(2)
938 0.023104993115419
939 >>> xi = lambda s: 0.5*s*(s-1)*pi**(-0.5*s)*gamma(0.5*s)*zeta(s)
940 >>> Xi = lambda t: xi(0.5+t*j)
941 >>> chop(-0.5*diff(Xi,0,n=2)/Xi(0))
942 0.023104993115419
944 We may ask for an approximate error value::
946 >>> secondzeta(0.5+100j, error=True)
947 ((-0.216272011276718 - 0.844952708937228j), 2.22044604925031e-16)
949 The function has poles at the negative odd integers,
950 and dyadic rational values at the negative even integers::
952 >>> mp.dps = 30
953 >>> secondzeta(-8)
954 -0.67236328125
955 >>> secondzeta(-7)
956 +inf
958 **Implementation notes**
960 The function is computed as sum of four terms `Z(s)=A(s)-P(s)+E(s)-S(s)`
961 respectively main, prime, exponential and singular terms.
962 The main term `A(s)` is computed from the zeros of zeta.
963 The prime term depends on the von Mangoldt function.
964 The singular term is responsible for the poles of the function.
966 The four terms depends on a small parameter `a`. We may change the
967 value of `a`. Theoretically this has no effect on the sum of the four
968 terms, but in practice may be important.
970 A smaller value of the parameter `a` makes `A(s)` depend on
971 a smaller number of zeros of zeta, but `P(s)` uses more values of
972 von Mangoldt function.
974 We may also add a verbose option to obtain data about the
975 values of the four terms.
977 >>> mp.dps = 10
978 >>> secondzeta(0.5 + 40j, error=True, verbose=True)
979 main term = (-30190318549.138656312556 - 13964804384.624622876523j)
980 computed using 19 zeros of zeta
981 prime term = (132717176.89212754625045 + 188980555.17563978290601j)
982 computed using 9 values of the von Mangoldt function
983 exponential term = (542447428666.07179812536 + 362434922978.80192435203j)
984 singular term = (512124392939.98154322355 + 348281138038.65531023921j)
985 ((0.059471043 + 0.3463514534j), 1.455191523e-11)
987 >>> secondzeta(0.5 + 40j, a=0.04, error=True, verbose=True)
988 main term = (-151962888.19606243907725 - 217930683.90210294051982j)
989 computed using 9 zeros of zeta
990 prime term = (2476659342.3038722372461 + 28711581821.921627163136j)
991 computed using 37 values of the von Mangoldt function
992 exponential term = (178506047114.7838188264 + 819674143244.45677330576j)
993 singular term = (175877424884.22441310708 + 790744630738.28669174871j)
994 ((0.059471043 + 0.3463514534j), 1.455191523e-11)
996 Notice the great cancellation between the four terms. Changing `a`, the
997 four terms are very different numbers but the cancellation gives
998 the good value of Z(s).
1000 **References**
1002 A. Voros, Zeta functions for the Riemann zeros, Ann. Institute Fourier,
1003 53, (2003) 665--699.
1005 A. Voros, Zeta functions over Zeros of Zeta Functions, Lecture Notes
1006 of the Unione Matematica Italiana, Springer, 2009.
1007 """
1008 s = ctx.convert(s)
1009 a = ctx.convert(a)
1010 tol = ctx.eps
1011 if ctx.isint(s) and ctx.re(s) <= 1:
1012 if abs(s-1) < tol*1000:
1013 return ctx.inf
1014 m = int(round(ctx.re(s)))
1015 if m & 1:
1016 return ctx.inf
1017 else:
1018 return ((-1)**(-m//2)*\
1019 ctx.fraction(8-ctx.eulernum(-m,exact=True),2**(-m+3)))
1020 prec = ctx.prec
1021 try:
1022 t3 = secondzeta_exp_term(ctx, s, a)
1023 extraprec = max(ctx.mag(t3),0)
1024 ctx.prec += extraprec + 3
1025 t1, r1, gt = secondzeta_main_term(ctx,s,a,error='True', verbose='True')
1026 t2, r2, pt = secondzeta_prime_term(ctx,s,a,error='True', verbose='True')
1027 t4, r4 = secondzeta_singular_term(ctx,s,a,error='True')
1028 t3 = secondzeta_exp_term(ctx, s, a)
1029 err = r1+r2+r4
1030 t = t1-t2+t3-t4
1031 if kwargs.get("verbose"):
1032 print('main term =', t1)
1033 print(' computed using', gt, 'zeros of zeta')
1034 print('prime term =', t2)
1035 print(' computed using', pt, 'values of the von Mangoldt function')
1036 print('exponential term =', t3)
1037 print('singular term =', t4)
1038 finally:
1039 ctx.prec = prec
1040 if kwargs.get("error"):
1041 w = max(ctx.mag(abs(t)),0)
1042 err = max(err*2**w, ctx.eps*1.*2**w)
1043 return +t, err
1044 return +t
1047@defun_wrapped
1048def lerchphi(ctx, z, s, a):
1049 r"""
1050 Gives the Lerch transcendent, defined for `|z| < 1` and
1051 `\Re{a} > 0` by
1053 .. math ::
1055 \Phi(z,s,a) = \sum_{k=0}^{\infty} \frac{z^k}{(a+k)^s}
1057 and generally by the recurrence `\Phi(z,s,a) = z \Phi(z,s,a+1) + a^{-s}`
1058 along with the integral representation valid for `\Re{a} > 0`
1060 .. math ::
1062 \Phi(z,s,a) = \frac{1}{2 a^s} +
1063 \int_0^{\infty} \frac{z^t}{(a+t)^s} dt -
1064 2 \int_0^{\infty} \frac{\sin(t \log z - s
1065 \operatorname{arctan}(t/a)}{(a^2 + t^2)^{s/2}
1066 (e^{2 \pi t}-1)} dt.
1068 The Lerch transcendent generalizes the Hurwitz zeta function :func:`zeta`
1069 (`z = 1`) and the polylogarithm :func:`polylog` (`a = 1`).
1071 **Examples**
1073 Several evaluations in terms of simpler functions::
1075 >>> from mpmath import *
1076 >>> mp.dps = 25; mp.pretty = True
1077 >>> lerchphi(-1,2,0.5); 4*catalan
1078 3.663862376708876060218414
1079 3.663862376708876060218414
1080 >>> diff(lerchphi, (-1,-2,1), (0,1,0)); 7*zeta(3)/(4*pi**2)
1081 0.2131391994087528954617607
1082 0.2131391994087528954617607
1083 >>> lerchphi(-4,1,1); log(5)/4
1084 0.4023594781085250936501898
1085 0.4023594781085250936501898
1086 >>> lerchphi(-3+2j,1,0.5); 2*atanh(sqrt(-3+2j))/sqrt(-3+2j)
1087 (1.142423447120257137774002 + 0.2118232380980201350495795j)
1088 (1.142423447120257137774002 + 0.2118232380980201350495795j)
1090 Evaluation works for complex arguments and `|z| \ge 1`::
1092 >>> lerchphi(1+2j, 3-j, 4+2j)
1093 (0.002025009957009908600539469 + 0.003327897536813558807438089j)
1094 >>> lerchphi(-2,2,-2.5)
1095 -12.28676272353094275265944
1096 >>> lerchphi(10,10,10)
1097 (-4.462130727102185701817349e-11 - 1.575172198981096218823481e-12j)
1098 >>> lerchphi(10,10,-10.5)
1099 (112658784011940.5605789002 - 498113185.5756221777743631j)
1101 Some degenerate cases::
1103 >>> lerchphi(0,1,2)
1104 0.5
1105 >>> lerchphi(0,1,-2)
1106 -0.5
1108 Reduction to simpler functions::
1110 >>> lerchphi(1, 4.25+1j, 1)
1111 (1.044674457556746668033975 - 0.04674508654012658932271226j)
1112 >>> zeta(4.25+1j)
1113 (1.044674457556746668033975 - 0.04674508654012658932271226j)
1114 >>> lerchphi(1 - 0.5**10, 4.25+1j, 1)
1115 (1.044629338021507546737197 - 0.04667768813963388181708101j)
1116 >>> lerchphi(3, 4, 1)
1117 (1.249503297023366545192592 - 0.2314252413375664776474462j)
1118 >>> polylog(4, 3) / 3
1119 (1.249503297023366545192592 - 0.2314252413375664776474462j)
1120 >>> lerchphi(3, 4, 1 - 0.5**10)
1121 (1.253978063946663945672674 - 0.2316736622836535468765376j)
1123 **References**
1125 1. [DLMF]_ section 25.14
1127 """
1128 if z == 0:
1129 return a ** (-s)
1130 # Faster, but these cases are useful for testing right now
1131 if z == 1:
1132 return ctx.zeta(s, a)
1133 if a == 1:
1134 return ctx.polylog(s, z) / z
1135 if ctx.re(a) < 1:
1136 if ctx.isnpint(a):
1137 raise ValueError("Lerch transcendent complex infinity")
1138 m = int(ctx.ceil(1-ctx.re(a)))
1139 v = ctx.zero
1140 zpow = ctx.one
1141 for n in xrange(m):
1142 v += zpow / (a+n)**s
1143 zpow *= z
1144 return zpow * ctx.lerchphi(z,s, a+m) + v
1145 g = ctx.ln(z)
1146 v = 1/(2*a**s) + ctx.gammainc(1-s, -a*g) * (-g)**(s-1) / z**a
1147 h = s / 2
1148 r = 2*ctx.pi
1149 f = lambda t: ctx.sin(s*ctx.atan(t/a)-t*g) / \
1150 ((a**2+t**2)**h * ctx.expm1(r*t))
1151 v += 2*ctx.quad(f, [0, ctx.inf])
1152 if not ctx.im(z) and not ctx.im(s) and not ctx.im(a) and ctx.re(z) < 1:
1153 v = ctx.chop(v)
1154 return v