# Gram points computed by iterative formula

Jun 18th, 2017
212
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. (* Start *)
2. (* Gram points to arbitrary precision *)
3. c = 7/8;
4. x = N[Table[2*Pi*Exp[1]*Exp[ProductLog[(n - c)/Exp[1]]], {n, 1, 42}],
5. 40];
6. Do[n = N[Round[(x/(2*Pi)*Log[x/(2*Pi*Exp[1])] - 1/8 -
7. RiemannSiegelTheta[x]/Pi) + Range[42], 10^-20], 20];
8. x = N[Round[2*Pi*Exp[1]*Exp[ProductLog[(n - c)/Exp[1]]], 10^-20],
9. 20];, {i, 1, 20}]
10. n
11. x
12. RiemannSiegelTheta[x]/Pi
13. Round[x]
14.
15.
16. (* start *)
17. (* Gram points *)
18. Table[N[InverseFunction[RiemannSiegelTheta][Pi*(n - 1)], 30], {n, 0,
19. 12}]
20. Round[%]
21. (* end *)
22.
23.
24. (* start *)
25. (* complementary Gram points *)
26. Table[N[InverseFunction[RiemannSiegelTheta][Pi*(n - 1/2)], 30], {n, 0,
27. 12}]
28. Round[%]
29. (* end *)
30.
31.
32. (* Franca-Leclair points to arbitrary precision *)
33. c = 11/8;
34. x = N[Table[2*Pi*Exp[1]*Exp[ProductLog[(n - c)/Exp[1]]], {n, 1, 42}],
35. 40];
36. Do[n = N[Round[(x/(2*Pi)*Log[x/(2*Pi*Exp[1])] - 1/8 -
37. RiemannSiegelTheta[x]/Pi) + Range[42], 10^-20], 20];
38. x = N[Round[2*Pi*Exp[1]*Exp[ProductLog[(n - c)/Exp[1]]], 10^-20],
39. 20];, {i, 1, 12}]
40. n
41. x
42. RiemannSiegelTheta[x]/Pi
43. Round[x]
44. (* End *)
45.
46. (* Start *)
47. (*Setting the variable start=1 gives the minus 1-th Gram \
48. point.*)Clear[n]; start = 1; end = 42;(*Interesting values of the \
49. variable c are:c=0 c=1/2 c=1/4 c=3/4*)c = 0; numberOfIterations = \
50. 300; decimaldigits = 100; Monitor[
51. x = Table[x = 1;
52. Do[x = N[
53. Round[2*Pi*
54. E^(ProductLog[-((c - n +
55. RiemannSiegelTheta[x]/
56. Pi + (x*(-Log[x] + 1 + Log[2] + Log[Pi]))/(2*
57. Pi) + 2)/E)] + 1), 10^(-decimaldigits)],
58. decimaldigits];, {i, 1, numberOfIterations}];
59. x, {n, start, end}];, n]; Zeta[
60. 1/2 + I*x]; Round[x] (*_Mats Granvik_,Jun 17 2017*)
61. (* End *)
62.
63. (* Start *)
64. (* Mathematica *)
65. (* Derivation of the Gram point approximation by Weisstein in \
66. Mathworld, https://oeis.org/A002505 *)
67. Clear[x, n, a, g];
68. Series[RiemannSiegelTheta[x], {x, Infinity, 12}]
69. a = Normal[Series[RiemannSiegelTheta[x], {x, Infinity, 0}]]
70. g = FullSimplify[(x /. Solve[a == n*Pi, x])[[1]]]
71. (* g above is equivalent to g[n_] = 2*E^(1 + ProductLog[(1 + \
72. 8*n)/(8*E)])*Pi *)
73. n = Range[42] - 2;
74. N[g, 20]
75. Round[g]
76. (* End *)
RAW Paste Data