(*Mathematica 8.0.1*)
(*Program 1*)
(*Series expansion of the Riemann zeta function in the critical strip \
with imaginary part greater than 10 and arbitrary real part.The \
series expansion is done in two steps:First as a Taylor polynomial to \
length firstNN=10 and then of the reciprocal of that Taylor \
polynomial to length secondNN=20.*)(*Mathematica 8.0.1*)
(*start*)
Clear[cons, ww, div, real, secondNN, firstNN, zz, z, t, b, n, k, nn, 
  x, g1, g2];
cons = 10;
ww = 500;
div = 10;
real = 0;
(*secondNN must be much greater than firstNN*)
firstNN = 10;
secondNN = 20;
Monitor[TableForm[zz = Table[Clear[t, b, n, k, nn, x];
     z = N[cons + w/div, 20];
     polynomial = 
      Normal[Series[Zeta[(real + x + z*N[I, 20])], {x, 0, firstNN}]];
     digits = 20;
     b = With[{nn = secondNN}, 
       CoefficientList[Series[1/polynomial, {x, 0, nn}], x]];
     nn = Length[b];
     x = z*I + N[b[[nn - 1]]/b[[nn]], digits], {w, 0, ww}]];, w]
"Plot of the real part:"
g1 = ListLinePlot[Flatten[Re[zz]], DataRange -> {cons, cons + ww/div},
   PlotRange -> {-2, 2}]
"Plot of the imaginary part:"
g2 = ListLinePlot[Flatten[Im[zz]], DataRange -> {cons, cons + ww/div}]
zz
(*end*)


(*Program 2*)
(*The limiting ratio of the form (n-1)*a (n-1)/a(n) applied to the \
first column of the matrix inverse of the binomial matrix multiplied \
elementwise with the coefficients in the Taylor polynomial of the \
Riemann zeta function,expanded at an arbitrary point in the critical \
strip (with imaginary part greater than 10),appended with trailing \
zeros.The trailing zeros...0,0,0,0,0,0,0,... in the vector "a" give \
better convergence to the Riemann zeta zeros*)
(*start*)
Clear[start, end, realPart, zz, d, digits, a, nn, A, n, k, b, z, list,
   x];
start = 10;
end = 60;
realPart = 0;
(*realPart=1/2;*)
(*realPart=1;*)
(*d must be much greater than zz*)
Monitor[list = Table[zz = 10;
   d = 20;
   digits = 30;
   a = Flatten[{CoefficientList[
        Normal[Series[
          Zeta[realPart + x + z*N[I, digits]], {x, 0, zz}]], x]*
       Range[0, zz]!, Range[d]*0}];
   nn = Length[a];
   A = Table[
     Table[If[n >= k, Binomial[n - 1, k - 1]*a[[n - k + 1]], 0], {k, 
       1, Length[a]}], {n, 1, Length[a]}];
   Quiet[b = Inverse[A][[All, 1]]];
   z*I + N[(nn - 1)*b[[nn - 1]]/b[[nn]], digits], {z, start, end, 
    1/10}], z*10]
"Plot of the real part:"
ListLinePlot[Re[list], PlotRange -> {-1, 3}, DataRange -> {start, end}]
"Plot of the imaginary part:"
ListLinePlot[Im[list], DataRange -> {start, end}]
(*end*)


(*Program 3*)
(*Newton Raphson iteration*)
(*start*)
Clear[f, g, t, min, max];
f[t_] = Zeta[N[1/2] + I*t]/Zeta'[1/2 + I*t];
g[t_] = Im[(1/2 + I*t) - f[t]];
min = 10;
max = 60;
"Plot of the imaginary part"
Plot[g[g[g[g[t]]]], {t, min, max}]
(*Plot[g[g[g[g[g[g[g[t]]]]]]],{t,min,max}]*)
(*end*)

(*Program 4*)
(*An approximation of the logarithmic derivative of the Riemann zeta \
function giving faster Newton Raphson iteration*)
(*The approximation of the logarithmic derivative is based on the \
formula:Limit[Zeta[s]*Zeta[c]/Zeta[s+c+-1]-Zeta[c],c->1]^(-1)=-Zeta[s]\
/Zeta'[s]*)
(*start*)
Clear[f, g, t, s, c, min, max];
c = 1 + 1/1000;
f[t_] = -1/((Zeta[N[1/2] + I*t]*Zeta[c])/Zeta[1/2 + I*t + c - 1] - 
     Zeta[c]);
g[t_] = Im[(1/2 + I*t) - f[t]];
min = 10;
max = 60;
"Plot of imaginary part"
Plot[g[g[g[g[t]]]], {t, min, max}]
(*Plot[g[g[g[g[g[g[g[t]]]]]]],{t,min,max}]*)
(*end*)