zetam001.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 > 0
    l = [0.5, -1.0 / (a + 1.0)]
    for i in range(1, n):
        l.append(0.5 * ((i + 0.0) ** a + (i + 1.0) ** a))
        l.append(-df(i, a + 1.0) / (a + 1.0))
    bs = [1.0 / 6, -1.0 / 30, 1.0 / 42, -1.0 / 30, 5.0 / 66, -691.0 / 2730, 7.0 / 6]
    t0 = 0.5 * a
    for i in range(len(bs)):
        a0 = a - 1 - 2 * i
        l.append(-bs[i] * t0 * (n ** a0))
        t0 *= a0 * (a0 - 1.0) / ((i * 2 + 3) * (i * 2 + 4))

    return math.fsum(l)

zetam_vs = {
    2.5 : 0.00851692877785033054,
    2.0 : 0,
    0.5 : -0.207886224977354566,
    0.0 : -0.5,
    -0.5 : -1.4603545088095868,
}



a = 0.5
n = 9

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

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