zetam002.py

import math


#calculates (x+1)^a - x^a
def df(x, a):
    lm = 0.5 * a * math.log(x * (x + 1.0))
    ld = 0.5 * a * math.log1p(1.0 / x)
    return 2 * math.exp(lm) * math.sinh(ld)


def zetam(a, n):
    a += 0.0
    assert a != -1

    l = []
    for i in range(n):
        l.append(df(2 * i + 1.0, a))
    l.append(-0.5 * (2 * n + 1.0) ** a)

    q = 20
    l1 = [0.25 * ((-1) ** i) * df(2 * n + i + 1.0, a) for i in range(q)]
    for i in range(q):
        l.append(l1[0])
        for j in range(q - i - 1):
            l1[j] = 0.5 * (l1[j] + l1[j + 1])

    return math.fsum(l) / df(1.0, a + 1.0)

zetam_vs = {
    2.5 : 0.00851692877785033054,
    2.0 : 0,
    0.5 : -0.207886224977354566,
    0.0 : -0.5,
    -0.5 : -1.4603545088095868,
    -2.0 : math.pi ** 2 / 6
}


a = -2.0
n = 6
q = 20

print("a =", a, ", n =", n, ", q =", q)
r0 = zetam(a, n)
rd = zetam_vs.get(a, 0.0)

print("{:.30f}".format(r0), r0 - rd)