"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "mpmath/functions/zeta.py" between
mpmath-1.0.0.tar.gz and mpmath-1.1.0.tar.gz

About: mpmath is a Python library for arbitrary-precision floating-point arithmetic.

zeta.py  (mpmath-1.0.0):zeta.py  (mpmath-1.1.0)
skipping to change at line 447 skipping to change at line 447
else: else:
raise ValueError raise ValueError
if ctx._is_real_type(z) and z < 0: if ctx._is_real_type(z) and z < 0:
l = ctx._re(l) l = ctx._re(l)
return l return l
def polylog_general(ctx, s, z): def polylog_general(ctx, s, z):
v = ctx.zero v = ctx.zero
u = ctx.ln(z) u = ctx.ln(z)
if not abs(u) < 5: # theoretically |u| < 2*pi if not abs(u) < 5: # theoretically |u| < 2*pi
raise NotImplementedError("polylog for arbitrary s and z") j = ctx.j
v = 1-s
y = ctx.ln(-z)/(2*ctx.pi*j)
return ctx.gamma(v)*(j**v*ctx.zeta(v,0.5+y) + j**-v*ctx.zeta(v,0.5-y))/(
2*ctx.pi)**v
t = 1 t = 1
k = 0 k = 0
while 1: while 1:
term = ctx.zeta(s-k) * t term = ctx.zeta(s-k) * t
if abs(term) < ctx.eps: if abs(term) < ctx.eps:
break break
v += term v += term
k += 1 k += 1
t *= u t *= u
t /= k t /= k
skipping to change at line 482 skipping to change at line 485
if s == -1: if s == -1:
return z/(1-z)**2 return z/(1-z)**2
if abs(z) <= 0.75 or (not ctx.isint(s) and abs(z) < 0.9): if abs(z) <= 0.75 or (not ctx.isint(s) and abs(z) < 0.9):
return polylog_series(ctx, s, z) return polylog_series(ctx, s, z)
if abs(z) >= 1.4 and ctx.isint(s): if abs(z) >= 1.4 and ctx.isint(s):
return (-1)**(s+1)*polylog_series(ctx, s, 1/z) + polylog_continuation(ct x, s, z) return (-1)**(s+1)*polylog_series(ctx, s, 1/z) + polylog_continuation(ct x, s, z)
if ctx.isint(s): if ctx.isint(s):
return polylog_unitcircle(ctx, int(s), z) return polylog_unitcircle(ctx, int(s), z)
return polylog_general(ctx, s, z) return polylog_general(ctx, s, z)
#raise NotImplementedError("polylog for arbitrary s and z")
# This could perhaps be used in some cases
#from quadrature import quad
#return quad(lambda t: t**(s-1)/(exp(t)/z-1),[0,inf])/gamma(s)
@defun_wrapped @defun_wrapped
def clsin(ctx, s, z, pi=False): def clsin(ctx, s, z, pi=False):
if ctx.isint(s) and s < 0 and int(s) % 2 == 1: if ctx.isint(s) and s < 0 and int(s) % 2 == 1:
return z*0 return z*0
if pi: if pi:
a = ctx.expjpi(z) a = ctx.expjpi(z)
else: else:
a = ctx.expj(z) a = ctx.expj(z)
if ctx._is_real_type(z) and ctx._is_real_type(s): if ctx._is_real_type(z) and ctx._is_real_type(s):
return ctx.im(ctx.polylog(s,a)) return ctx.im(ctx.polylog(s,a))
skipping to change at line 628 skipping to change at line 626
# TODO: implement for derivatives # TODO: implement for derivatives
if d != 0: if d != 0:
raise NotImplementedError raise NotImplementedError
res = ctx.re(s) res = ctx.re(s)
negs = -s negs = -s
# Integer reflection formula # Integer reflection formula
if ctx.isnpint(s): if ctx.isnpint(s):
n = int(res) n = int(res)
if n <= 0: if n <= 0:
return ctx.bernpoly(1-n, a) / (n-1) return ctx.bernpoly(1-n, a) / (n-1)
if not (atype == 'Q' or atype == 'Z'):
raise NotImplementedError
t = 1-s t = 1-s
# We now require a to be standardized # We now require a to be standardized
v = 0 v = 0
shift = 0 shift = 0
b = a b = a
while ctx.re(b) > 1: while ctx.re(b) > 1:
b -= 1 b -= 1
v -= b**negs v -= b**negs
shift -= 1 shift -= 1
while ctx.re(b) <= 0: while ctx.re(b) <= 0:
v += b**negs v += b**negs
b += 1 b += 1
shift += 1 shift += 1
# Rational reflection formula # Rational reflection formula
if atype == 'Q' or atype == 'Z': try:
try: p, q = a._mpq_
p, q = a._mpq_ except:
except: assert a == int(a)
assert a == int(a) p = int(a)
p = int(a) q = 1
q = 1 p += shift*q
p += shift*q assert 1 <= p <= q
assert 1 <= p <= q g = ctx.fsum(ctx.cospi(t/2-2*k*b)*ctx._hurwitz(t,(k,q)) \
g = ctx.fsum(ctx.cospi(t/2-2*k*b)*ctx._hurwitz(t,(k,q)) \ for k in range(1,q+1))
for k in range(1,q+1)) g *= 2*ctx.gamma(t)/(2*ctx.pi*q)**t
g *= 2*ctx.gamma(t)/(2*ctx.pi*q)**t v += g
v += g return v
return v
# General reflection formula
# Note: clcos/clsin can raise NotImplementedError
else:
C1, C2 = ctx.cospi_sinpi(0.5*t)
# Clausen functions; could maybe use polylog directly
if C1: C1 *= ctx.clcos(t, 2*a, pi=True)
if C2: C2 *= ctx.clsin(t, 2*a, pi=True)
v += 2*ctx.gamma(t)/(2*ctx.pi)**t*(C1+C2)
return v
def _hurwitz_em(ctx, s, a, d, prec, verbose): def _hurwitz_em(ctx, s, a, d, prec, verbose):
# May not be converted at this point # May not be converted at this point
a = ctx.convert(a) a = ctx.convert(a)
tol = -prec tol = -prec
# Estimate number of terms for Euler-Maclaurin summation; could be improved # Estimate number of terms for Euler-Maclaurin summation; could be improved
M1 = 0 M1 = 0
M2 = prec // 3 M2 = prec // 3
N = M2 N = M2
lsum = 0 lsum = 0
skipping to change at line 693 skipping to change at line 683
# l = ctx.fsum((-ctx.ln(n+a))**d * (n+a)**negs for n in range(M1,M2)) # l = ctx.fsum((-ctx.ln(n+a))**d * (n+a)**negs for n in range(M1,M2))
#else: #else:
# l = ctx.fsum((n+a)**negs for n in range(M1,M2)) # l = ctx.fsum((n+a)**negs for n in range(M1,M2))
lsum += l lsum += l
M2a = M2+a M2a = M2+a
logM2a = ctx.ln(M2a) logM2a = ctx.ln(M2a)
logM2ad = logM2a**d logM2ad = logM2a**d
logs = [logM2ad] logs = [logM2ad]
logr = 1/logM2a logr = 1/logM2a
rM2a = 1/M2a rM2a = 1/M2a
M2as = rM2a**s M2as = M2a**(-s)
if d: if d:
tailsum = ctx.gammainc(d+1, s1*logM2a) / s1**(d+1) tailsum = ctx.gammainc(d+1, s1*logM2a) / s1**(d+1)
else: else:
tailsum = 1/((s1)*(M2a)**s1) tailsum = 1/((s1)*(M2a)**s1)
tailsum += 0.5 * logM2ad * M2as tailsum += 0.5 * logM2ad * M2as
U = [1] U = [1]
r = M2as r = M2as
fact = 2 fact = 2
for j in range(1, N+1): for j in range(1, N+1):
# TODO: the following could perhaps be tidied a bit # TODO: the following could perhaps be tidied a bit
 End of changes. 5 change blocks. 
30 lines changed or deleted 21 lines changed or added

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)