Gram points computed by iterative formula

(* Start *)
(* Gram points to arbitrary precision *)
c = 7/8;
x = N[Table[2*Pi*Exp[1]*Exp[ProductLog[(n - c)/Exp[1]]], {n, 1, 42}],
   40];
Do[n = N[Round[(x/(2*Pi)*Log[x/(2*Pi*Exp[1])] - 1/8 -
       RiemannSiegelTheta[x]/Pi) + Range[42], 10^-20], 20];
 x = N[Round[2*Pi*Exp[1]*Exp[ProductLog[(n - c)/Exp[1]]], 10^-20],
   20];, {i, 1, 20}]
n
x
RiemannSiegelTheta[x]/Pi
Round[x]


(* start *)
(* Gram points *)
Table[N[InverseFunction[RiemannSiegelTheta][Pi*(n - 1)], 30], {n, 0,
  12}]
Round[%]
(* end *)


(* start *)
(* complementary Gram points *)
Table[N[InverseFunction[RiemannSiegelTheta][Pi*(n - 1/2)], 30], {n, 0,
   12}]
Round[%]
(* end *)


(* Franca-Leclair points to arbitrary precision *)
c = 11/8;
x = N[Table[2*Pi*Exp[1]*Exp[ProductLog[(n - c)/Exp[1]]], {n, 1, 42}],
   40];
Do[n = N[Round[(x/(2*Pi)*Log[x/(2*Pi*Exp[1])] - 1/8 -
       RiemannSiegelTheta[x]/Pi) + Range[42], 10^-20], 20];
 x = N[Round[2*Pi*Exp[1]*Exp[ProductLog[(n - c)/Exp[1]]], 10^-20],
   20];, {i, 1, 12}]
n
x
RiemannSiegelTheta[x]/Pi
Round[x]
(* End *)

(* Start *)
(*Setting the variable start=1 gives the minus 1-th Gram \
point.*)Clear[n]; start = 1; end = 42;(*Interesting values of the \
variable c are:c=0 c=1/2 c=1/4 c=3/4*)c = 0; numberOfIterations = \
300; decimaldigits = 100; Monitor[
 x = Table[x = 1;
    Do[x = N[
        Round[2*Pi*
          E^(ProductLog[-((c - n +
                   RiemannSiegelTheta[x]/
                    Pi + (x*(-Log[x] + 1 + Log[2] + Log[Pi]))/(2*
                    Pi) + 2)/E)] + 1), 10^(-decimaldigits)],
        decimaldigits];, {i, 1, numberOfIterations}];
    x, {n, start, end}];, n]; Zeta[
 1/2 + I*x]; Round[x] (*_Mats Granvik_,Jun 17 2017*)
(* End *)

(* Start *)
(* Mathematica *)
(* Derivation of the Gram point approximation by Weisstein in \
Mathworld, https://oeis.org/A002505 *)
Clear[x, n, a, g];
Series[RiemannSiegelTheta[x], {x, Infinity, 12}]
a = Normal[Series[RiemannSiegelTheta[x], {x, Infinity, 0}]]
g = FullSimplify[(x /. Solve[a == n*Pi, x])[[1]]]
(* g above is equivalent to g[n_] = 2*E^(1 + ProductLog[(1 + \
8*n)/(8*E)])*Pi *)
n = Range[42] - 2;
N[g, 20]
Round[g]
(* End *)