Advertisement
am1x

zetam001.py

Aug 14th, 2022
829
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 0.85 KB | Source Code | 0 0
  1. import math
  2.  
  3.  
  4. #calculates (x+1)^a - x^a
  5. def df(x, a):
  6.     lm = 0.5 * a * math.log(x * (x + 1.0))
  7.     ld = 0.5 * a * math.log1p(1.0 / x)
  8.     return 2 * math.exp(lm) * math.sinh(ld)
  9.  
  10.  
  11. def zetam(a, n):
  12.     a += 0.0
  13.     assert a + 1 > 0
  14.     l = [0.5, -1.0 / (a + 1.0)]
  15.     for i in range(1, n):
  16.         l.append(0.5 * ((i + 0.0) ** a + (i + 1.0) ** a))
  17.         l.append(-df(i, a + 1.0) / (a + 1.0))
  18.     bs = [1.0 / 6, -1.0 / 30, 1.0 / 42, -1.0 / 30, 5.0 / 66, -691.0 / 2730, 7.0 / 6]
  19.     t0 = 0.5 * a
  20.     for i in range(len(bs)):
  21.         a0 = a - 1 - 2 * i
  22.         l.append(-bs[i] * t0 * (n ** a0))
  23.         t0 *= a0 * (a0 - 1.0) / ((i * 2 + 3) * (i * 2 + 4))
  24.  
  25.     return math.fsum(l)
  26.  
  27. zetam_vs = {
  28.     2.5 : 0.00851692877785033054,
  29.     2.0 : 0,
  30.     0.5 : -0.207886224977354566,
  31.     0.0 : -0.5,
  32.     -0.5 : -1.4603545088095868,
  33. }
  34.  
  35.  
  36.  
  37. a = 0.5
  38. n = 9
  39.  
  40. print("a = ", a, ", n = ", n)
  41. r0 = zetam(a, n)
  42. rd = zetam_vs.get(a)
  43.  
  44. print("{:.30f}".format(r0), r0 - rd)
  45.  
  46.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement