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)