[python3] Calc Riemann Zeta function

complex.py

# complex arith for programming with other languages
# - required functions: exp(f), log(f), sin(f), cos(f), atan2(f), pow(f1, f2)
import math


# [equality for complex]
def ceq(a, b):
    return a.real == b.real and a.imag == b.imag

# [add, sub, mul for complex]
def cadd(a, b):
    a, b = complex(a), complex(b)
    return (a.real + b.real) + (a.imag + b.imag) * 1j

def csub(a, b):
    a, b = complex(a), complex(b)
    return (a.real - b.real) + (a.imag - b.imag) * 1j

def cmul(a, b):
    a, b = complex(a), complex(b)
    real = a.real * b.real - a.imag * b.imag
    imag = a.real * b.imag + a.imag * b.real
    return real + imag * 1j

print([cadd(1+1j, 2+2j), (1+1j) + (2+2j)])
print([cmul(1+1j, 2+2j), (1+1j) * (2+2j)])


# [div for complex]
def conj(c):
    c = complex(c)
    return c.real - c.imag * 1j
def norm(c):
    c = complex(c)
    return c.real * c.real + c.imag * c.imag
def cinv(c):
    c = complex(c)
    # 1/(a+i*b) = (a-i*b)/[(a+i*b)*(a-i*b)]
    #           = (a-i*b)/(a*a + b*b)
    base = norm(c)
    return (c.real / base) - (c.imag / base) * 1j

def cdiv(a, b):
    #return cmul(a, cinv(b))
    base = norm(b)
    c = cmul(a, conj(b))
    return (c.real / base) + (c.imag / base) * 1j

print([cdiv(1+2j, 2+1j), (1+2j) / (2+1j)])


# [pow for complex]
def c2p(c):
    c = complex(c)
    r = math.pow(norm(c), 0.5)
    theta = math.atan2(c.imag, c.real)
    return r, theta # "r" and "theta" also called as "abs" and "arg"
def p2c(r, theta):
    return r * math.cos(theta) + r * math.sin(theta) * 1j

def cpow(a, b):
    r, t = c2p(a)
    b = complex(b)
    c, d = b.real, b.imag
    # a^b = (r*exp(i*t))^(c+i*d)
    #     = r^(c+i*d) * (exp(i*t))^(c+i*d)
    #     = r^c * r^(i*d) * exp(-d*t + i*c*t)
    #     = r^c * exp(log(r) * i*d) * exp(-d*t + i*c*t)
    #     = r^c * exp(log(r) * i*d -d*t + i*c*t)
    #     = r^c * exp(-d*t) * exp(i*(d*log(r) + c*t))
    #     = r^c * exp(-d*t) * [cos(d*log(r) + c*t) + i*sin(d*log(r) + c*t)]
    R = math.pow(r, c) * math.exp(-d * t)
    T = d * math.log(r) + c * t
    return p2c(R, T)

print([cpow(1+2j, 2+1j), (1+2j) ** (2+1j)])
print([cpow(2, 3), (2+0j) ** (3+0j)])


# [exp, log for complex]
def cexp(c):
    #return cpow(math.e, c)
    # exp(a+i*b) = exp(a) * exp(i*b)
    #            = exp(a) * [cos(b) + i*sin(b)]
    c = complex(c)
    return p2c(math.exp(c.real), c.imag)
    
def clog(c, n=0):
    r, t = c2p(c)
    # log(r*exp(i*t)) = log(r) + log(exp(i*t))
    #                 = log(r) + i*(t + 2*PI*n)
    # (n = ..., -2, -1, 0, +1, +2, ...)
    return math.log(r) + (t + 2 * n * math.pi) * 1j
 
print(cexp(clog(1+1j)))
print(clog(cexp(1+1j)))

# [Riemann sphere projection (single infinity)]
def cproj(c):
    inf = float("inf")
    if c.real == inf or c.real == -inf or c.imag == inf or c.imag == -inf:
        return complex("inf")
    return c

zeta-plot.py

# plotting zeta function on complex plane


# [zeta function impl]
from itertools import count, islice

def binom(n, k):
    v = 1
    for i in range(k):
        v *= (n - i) / (i + 1)
    return v


def zeta(s, t=100):
    if s == 1: return complex("inf")
    term = (1 / 2 ** (n + 1) * sum((-1) ** k * binom(n, k) * (k + 1) ** -s 
                                   for k in range(n + 1)) for n in count(0))
    return sum(islice(term, t)) / (1 - 2 ** (1 - s))


# [correct data with numpy]
import numpy

x = numpy.arange(-3, 3, 0.1)
y = numpy.arange(-30, 30, 1)
X, Y = numpy.meshgrid(x, y)

@numpy.vectorize
def z(real, imag):
    r = zeta(real + imag * 1j).real
    return min(max(-10, r), 10)

Z = z(X, Y)
#print(Z)


# [plot with matplotlib]
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import pyplot, cm

fig = pyplot.figure()
ax = fig.gca(projection="3d")
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm, 
                       linewidth=0, antialiased=False)
ax.set_zlim(-10, 10)


pyplot.show()

zeta-real-plot.png

zeta.py

# Programming for Calc Riemann Zeta function

from itertools import count, islice


# Basic definition of zeta func
# zeta(s) = sum(1/n^s)
# (Re(s) > 1)
def zeta1(s, t=10000):
    term = (1 / (n ** s) for n in count(1))
    return sum(islice(term, t))


print(zeta1(2)) # => 1.6449...
print((zeta1(2) * 6) ** 0.5) # => PI = 3.1415...
#print(zeta1(0.5+14.134725142j)) # invalid


# formula (20) in http://mathworld.wolfram.com/RiemannZetaFunction.html
#
#   sum((-1)^n/n^s) + sum(1/n^2) = 2 * sum(1/n^s, n=2,4,6,...)
#                                = 2 * sum(1/(2*n)^s)
#                                = 2 * 2^-s * sum(1/n^s)
#
#   sum((-1)^n/n^s) + zeta(s) = 2^(1-s) * zeta(s)
#
#   zeta(s) = 1/(1 - 2^(1-s)) * sum((-1^(n-1)) / n^s)
#   (Re(s) > 0 and s != 1)
def zeta2(s, t=10000):
    if s == 1: return float("inf")
    #term = ((-1)**(n - 1) / (n ** s) for n in count(1))
    #return sum(islice(term, t)) / (1 - 2**(1 - s))
    term = ((-1) ** n * n ** -s for n in count(1))
    return sum(islice(term, t)) / (2 ** (1 - s) -  1)


print(zeta2(2))
print((zeta2(2) * 6) ** 0.5)
print(abs(zeta2(0.5+14.134725142j))) # => 0
print(abs(zeta2(0.5-14.134725142j))) # => 0
#print(zeta2(0)) # invalid


# (utility) binomial coefficient
def binom(n, k):
    v = 1
    for i in range(k):
        v *= (n - i) / (i + 1)
    return v


# formula (21) in http://mathworld.wolfram.com/RiemannZetaFunction.html
# Global zeta function by Knopp and Hasse (s != 1)
def zeta3(s, t=100):
    if s == 1: return float("inf")
    term = (1 / 2 ** (n + 1) * sum((-1) ** k * binom(n, k) * (k + 1) ** -s 
                                   for k in range(n + 1)) for n in count(0))
    return sum(islice(term, t)) / (1 - 2 ** (1 - s))


print(zeta3(2))
print((zeta3(2) * 6) ** 0.5)
print(abs(zeta3(0.5+14.134725142j))) # => 0
print(abs(zeta3(0.5-14.134725142j))) # => 0
print(zeta3(1)) # => inf
print(zeta3(0)) # => -1/2
print(zeta3(-1)) # => -1/12 = 0.08333...
print(zeta3(-2)) # => 0

zeta.py.note.txt

# [Note on relationship of binomial and zeta function]

# [1. binomial definition inside]
#
# binom(a, b) as bn(a, b) = a!/((a-b)!*b!)
#
# bn(a, 0) = 1
# bn(a, 1) = (a - 0) / 1
# bn(a, 2) = (a-0)*(a-1) / (2*1)
# bn(a, 3) = (a-0)*(a-1)*(a-2) / (3*2*1)
# ...

# [2. case study: varying a side]
#
# bn(1,0)=1, bn(1,1)=1, bn(1, b >= 2)=0
# bn(2,0)=1, bin(2,1)=2, bn(2,2)=1, bn(2, b >= 3)=0
# bn(3,0)=1, bin(3,1)=3, bn(3,2)=3, bn(3,3)=1, bn(3, b >= 4)=0
# ...

# [3. apply negative a to bn(a,b) definition]
#
# bn(-1, 0) = 1
# bn(-1, 1) = (-1-0) / 1 = -1
# bn(-1, 2) = (-1-0)*(-1-1) / (2*1) = 1
# bn(-1, 3) = (-1-0)*(-1-1)*(-1-2) / (3*2*1) = -1
# ...
# bn(-1, k) = -1^k
#
# bn(-2, 0) = 1
# bn(-2, 1) = (-2-0) / 1 = -2
# bn(-2, 2) = (-2-0)*(-2-1) / (2*1) = 3
# bn(-2, 3) = (-2-0)*(-2-1)*(-2-2) / (3*2*1) = -4
# ...
# bn(-2, k) = (-1^k)*k
#
# bn(-a, k) = -a*(-a-1)*(-a-2)*...*(-a-(k-1)) / k!
#           = (-1^k) * (a+k-1)*(a+k-2)*...*(a-0) / k! ... (numerator count = k)
#           = (-1^k) * ((a+k-1)!/(a-1)!) / k!
#           = (-1^k) * (a+k-1)! / ((a-1)!*k!)
#           = (-1^k) * bn(a+k-1, k)
#
# bn(-a, 1) = -a
# bn(-a, 2) = a*(a+1)/2
# bn(-a, 3) = -a*(a+1)*(a+2)/3!
# ...

# [4. polynomial with binonial]
#
# (1+x)^1 = 1+x = bn(1,0)*1 + bn(1,1)*x
# (1+x)^2 = 1+2x+x^2 = bn(2,0)*1 + bn(2,1)*x + bn(2,2)*x^2
# (1+x)^3 = bn(3,0)*1 + bn(3,1)*x + bn(3,2)*x^2 + bn(3,3)*x^3
# ...
# (1+x)^k = sum(n=0..inf, bn(k, n) * x^n) 
#                  ... polynomial with binomial (sum up to infinity, not k)

# [5. negative polynomials]
#
# (1+x)^-1 = sum(n=0..inf, bn(-1, n) * x^n) = 1 - x + x^2 - x^3 + ...
# (1+x)^-2 = sum(n=0..inf, bn(-2, n) * x^n) = 1 - 2*x + 3*x^2 - 4*x^3 + ...
 
# [6. apply x=1 to the negative polynomials]
#
# (1+1)^-1 = 1 - 1 + 1 - 1 + ... 
#      1/2 =
# 
# (1+1)^-2 = 1 - 2 + 3 - 4 + ...
#      1/4 =

# [7. convert alternating sum to positive sums]
#
#  (1 - 1 + 1 - 1 + ... ) = (1 + 1 + 1 + 1 + ...) + (0 - 2 + 0 - 2 ...)
#                                               ... (focus only even terms)
#                         = (1 + 1 + 1 + 1 + ...) + (- 2 - 2 - 2 - 2 - ...)
#                         = (1 + 1 + 1 + 1 + ...) - 2*(1 + 1 + 1 + 1 +  ...)
#                     1/2 = - (1 + 1 + 1 + 1 + ...)
#                    -1/2 = 1 + 1 + 1 + 1 + ...  => zeta(0)
# 
# (1 - 2 + 3 - 4 + ...) = (1 + 2 + 3 + 4 + ...) + (0 - 4 + 0 - 8 + ...)
#                       = (1 + 2 + 3 + 4 + ...) + (- 4 - 8 - 12 - 16 -  ...)
#                       = (1 + 2 + 3 + 4 + ...) - 4*(1 + 2 + 3 + 4 + ...)
#                   1/4 = -3*(1 + 2 + 3 + 4 + ...)
#                 -1/12 = 1 + 2 + 3 + 4 + ...  => zeta(-1)

# [8. zeta function transform with binomial]
#
# zeta(s) = sum(n=1..inf, n^-s)
#         = 1 + 2^-s + sum(n=3..inf, n^-s)
#         = 1 + 2^-s + sum(n=2..inf, (n+1)^-s)
#         = 1 + 2^-s + sum(n=2..inf, n^-s * (1+n^-1)^-s)
#                                ... apply "polynomial with binomial"
#         = 1 + 2^-s + sum(n=2..inf, n^-s * sum(k=0..inf, bn(-s,k)*n^-k))
#                                ... swap sum n and k
#         = 1 + 2^-s + sum(k=0..inf, bn(-s, k)*sum(n=2..inf, n^-s * n^-k))
#         = 1 + 2^-s + sum(k=0..inf, bn(-s, k)*sum(n=2..inf, n^-(s+k)))
#                                ... apply zeta(s) = 1 + sum(n=2..inf, n^-s)
#         = 1 + 2^-s + sum(k=0..inf, bn(-s, k) * (zeta(s+k) - 1))
#                                ... extract k=0
#         = 1 + 2^-s + b(-s,0)*(zeta(s)-1) + 
#                                 sum(k=1..inf, bn(-s, k) * (zeta(s+k) - 1))
#         = 1 + 2^-s + zeta(s) - 1 + sum(k=1..inf, bn(-s, k) * (zeta(s+k) - 1))
#
#       0 = 2^-s + sum(k=1..inf, bn(-s, k) * (zeta(s+k) - 1))
#                                 ... extract sum
#       0 = 2^-s + bn(-s,1)*(zeta(s+1)-1) + bn(-s,2)*(zeta(s+2)-1) + ...
#       0 = 2^-s - s*(zeta(s+1)-1) + s*(s+1)/2*(zeta(s+2)-1) + ...
#                                 ... s => s-1
#       0 = 2^(1-s) - (s-1)*(zeta(s)-1) + (s-1)*s/2*(zeta(s+1)-1) - 
#                                     (s-1)*s*(s+1)/3!*(zeta(s+2)-1) + ...
#       0 = 2^(1-s)/(s-1) - zeta(s)+1 + s/2*(zeta(s+1)-1) -
#                                     s*(s+1)/3!*(zeta(s+2)-1) + ...
#
# zeta(s) = 2^(1-s)/(s-1) + 1 + 
#      sum(k=1..inf, -1^(k+1) * (s+k-1)!/((s-1)!k! * (k+1)) * (zeta(s+k)-1))
#
# zeta(s) = 2^(1-s)/(s-1) + 1 + sum(k=1..inf, -bn(-s,k)*(zeta(s+k)-1) / (k+1))

# [9. expand the domain of zeta(s)]
#
# If s > 0, s+k > 1 from k=1..inf.
# When zeta(s > 0), there are zeta(s+k > 1) of the rhs.
# So it can apply other zeta(s > 1) funcs: e.g. sum(n=1..inf, n^-s)
#
# -  zeta(s > 0) is generated by zeta(s > 1)
#
# In the same way, zeta(s > -1) is generated by zeta(s > 0), 
#                  zeta(s > -2) is generated by zeta(s > -1), ...
#
# -  zeta(s) is generated from zeta(s > 1)

# [generating zeta(s)]
 
# [10. Prepare when s->1]
#
# (s-1)*zeta(s) = 2^(1-s) + 
#                   (s-1)*{1 + sum(k=1..inf, -bn(-s,k)*(zeta(s+k)-1) / (k+1))}
#
# lim(s->1, (s-1)*zeta(s)) = 2^(1-1) + 
#                   (1-1)*{1 + sum(k=1..inf, -bn(-1,k)*(zeta(1+k)-1)/(k+1))}
#                          = 1 + 0 *
#                              {1 + sum(k=1..inf, -1^(k+1)*(zeta(1+k)-1)/(k+1)}
#                          = 1
#
# lim(s->1, (s-1)*zeta(s)) = 1
#
# (It is for tackling to extract zeta(1) terms especially as 0*zeta(1) case)

# [11. calc zeta(n) from function series]
#
# Apply term: 0*zeta(1) => 1 
#             0*zeta(s) => 0  (s != 1)
# 
# zeta(s) = 2^(1-s)/(s-1) + 1 + 
#      sum(k=1..inf, -1^(k+1) * (s+k-1)!/((s-1)!k! * (k+1)) * (zeta(s+k)-1))
#          = 2^(1-s)/(s-1) + 1 + 
#                   s/2*(zeta(s+1)-1) - s*(s+1)/3!*(zeta(s+2)-1) + 
#                   s*(s+1)*(s+2)/4!*(zeta(s+3)-1) - ...
#
# zeta(0) = 2^(1-0)/(0-1) + 1 + 1/2*(0*zeta(1)-0) - 1/3!(0*zeta(2)-0) + 
#                                                   1*2/4!*(0*zeta(3)-0) - ...
#         = -2 + 1 + 1/2 - 0 + 0 - ...
#         = -1/2
# 
# zeta(-1) = 2^(1+1)/(-1-1) + 1 + -1/2*(zeta(0)-1) - -1/3!*(0*zeta(1)-0) + 
#                                                   -1*1/4!*(0*zeta(2)-0) - ...
#          = -2 + 1 + -1/2*(-1/2-1) + 1/6*(1-0) + 0 - ...
#          = -24/12 + 12/12 + 9/12 + 2/12 + 0 - ...
#          = -1/12
#
# zeta(-2) = 2^(1+2)/(-2-1) + 1 + -2/2*(zeta(-1)-1) - -2*-1/3!*(zeta(0)-1) +
#                                             -2*-1/4!*(0*zeta(1)-0) - ...
#          = -8/3 + 1 + -1*(-1/12-1) - 1/3*(-1/2-1) + 1/12*(1-0) - 0 + ....
#          = -32/12 + 12/12 + 13/12 + 6/12 + 1/12 - 0 + ...
#          = 0