(*------------------------------------------------------------------------------------*)
(* Mathematica notebook to compute zeta(s) using Theorem 1.1 in .... .                  *) 
(* The code is intended for illustration and debugging mainly, so                       *)
(* it is not necessarily efficient.                                                     *)
(*------------------------------------------------------------------------------------*)

Clear[K,v,m,s0,s,t,σ,f,g,gDirect,gEM,gDer,emCor,emCorr,BlockVal,BlockValEM,BlockValDirect];

m = 4; (* number of derivative to use for g_K(z) *)
(* differentiate f_s0(z) w.r.t. z, then set z <- 0, treat s0 symbolically *)
f[s_,z_]=Exp[s z]/(1+z)^s; 
Table[fDer1[j]=Simplify[D[f[s0,z]/j!,{z,j}]/.z->0],{j,0,m}]; 
(* set s, then set s0 <- s in the derivatives of f_s(z) *)
σ=0.5`30;t =1000000.`30;s = σ+I t;
Table[fDer[j]=fDer1[j]/.s0->s,{j,0,m}]; 
(* set u0, v0, and M, which is the original length of the main sum *)
u0 = 6*Max[6,Ceiling[(Abs[s]+3)^0.5],Ceiling[σ]];
v0=10*(m+1)*u0; 
M = Max[v0,20Ceiling[Abs[s]+3]];   

(* compute the j-th derivative of g_K(z) w.r.t. z, use EM summation, useful for small z *) 
gEM[K_,z_,j_]:=If[j==0,NSum[Exp[k z],{k,0,K-1},Method->Integrate],NSum[k^j Exp[k z],{k,0,K-1},Method->Integrate]];
(* compute g_K(z) using formula for geometric sum *)
g[K_,z_]=(Exp[K z]-1)/(Exp[z]-1); 
(* compute derivatives of g_K(z) w.r.t. z by differentiating geometric sum formula *)
(* Note. fullsimplify seems to take too long if m is large, it can be replaced by simplify which is much faster *)
Table[gDer[K_,z_,j]=FullSimplify[D[g[K,z],{z,j}]],{j,0,m}];  

(* compute \sum_{k=0}^{K-1} (1+k/v)^{-s} in terms of derivatives of g_K(z) *)
BlockVal[z_,K_]:=Sum[fDer[j]*gDer[K,z,j]/v^j,{j,0,m}];
(* compute \sum_{k=0}^{K-1} (1+k/v)^{-s} in terms of derivatives of g_K(z) *)
BlockValEM[z_,K_]:=Sum[fDer[j]*gDirect[K,z,j]/v^j,{j,0,m}]; 

(* compute the j-th derivative of g_K(z) w.r.t. z use direct summation, used for debugging *) 
gDirect[K_,z_,j_]:=If[j==0,Sum[ Exp[k z],{k,0,K-1}],Sum[k^j  Exp[k z],{k,0,K-1}]];
(* compute \sum_{k=0}^{K-1} (1+k/v)^{-s} by direct summation, used for debugging *)
BlockValDirect[v_,K_]:=Sum[(1+k/v)^(-s),{k,0,K-1}];

(* use L1 correction terms in the Euler-Macaluarin formula for zeta *)
L1 = 6; Table[emCor[ell]=BernoulliB[2ell]/(2ell)!* M^(-s)*Product [(s+l)/M,{l,0,2ell-2}],{ell,1,L1}];

(* compute the bound for T_{M,m}(s,u0,v0) *)
R = -1;BBound = 0; v = v0; 
While[v<M-1,K = Ceiling[v/u0];BBound += Min[K,Abs[Csc[t/(2.v)]]]/v^σ;R+=1;v+=K;];
epsm=0.56*Exp[0.78(m+1)]/(m+1)^((m+1)/2)*Abs[s]^((m+1)/2)/u0^(m+1); 
BBound*=epsm;

(* main for loop to compute the main sum \sum_{n=1}^{M-1} n^(-s) *)
initSum=0.;remSum=0.;v = v0;
For[n=1,n< v0,n++,initSum += n^(-s)];
For[r=0,r≤R,r++,K =Min[ Ceiling[v/u0],M-v];
targ=Mod[t/v,2.Pi]; If[targ>Pi,targ-=2.Pi];
addVal= If[Abs[targ]<10(m+1)/(K-1), BlockValEM[-(σ/v+I  targ),K],BlockVal[-(σ/v+I targ),K]];
remSum+=addVal*v^(-s); v+=K];
bpts = M^(-s)/2+M^(1-s)/(s-1); (* boundary term *)
emCorr=Sum[emCor[ell],{ell,1,L1}]; (* Euler-Maclaurin correction terms *)
zetaVal = initSum+remSum+bpts+emCorr
Abs[Zeta[s]-zetaVal]