Euler Maclaurin formula for Riemann zeta function

(*start*)
Clear[n, k, s, r, q];
Print["choose a complex number to calculate for"]
s = 1/2 - I*100
Print["Increase these number for better precision"]
k = 1000
q = 100
Print["The Euler Maclaurin formula compared to Mathematica 8 inbuilt \
formula"]
N[Sum[1/n^s, {n, 1, k}] + k^(1 - s)/(s - 1) - (k^(-s))/2 +
  Sum[BernoulliB[2*r]/((2*r)!)*Product[s + i, {i, 0, 2*r - 2}]*
    k^(-s - 2*r + 1), {r, 1, q - 1}], 20]
N[Zeta[s], 20]
%% - %
Print["The difference above is zero"]
(*end*)

(*start 28.1.2018*)
Clear[n, k, s, r, q];
Print["choose a complex number to calculate for"]
s = 1/2 - I*100
Print["Increase these number for better precision"]
k = 100
q = 100
Print["The Euler Maclaurin formula compared to Mathematica 8 inbuilt \
formula"]
s = 1/2 + I*10;
N[Sum[1/n^s, {n, 1, k}] + k^(1 - s)/(s - 1) +
  Sum[BernoulliB[r]/(r!)*Product[s + i, {i, 0, r - 2}]*
    k^(-s - r + 1), {r, 1, 2*(q - 1)}], 20]
N[Zeta[s], 20]
%% - %
Print["The difference above is zero"]
(*end*)
(* In agreement with *)
(*start*)
nn = 12;
TableForm[
 A = Table[
   Table[If[n >= k, BernoulliB[n - k], 0], {k, 1, nn}], {n, 1, nn}]]
TableForm[
  K = Inverse[
    Table[Table[
      If[n >= k, Binomial[n, k]/(n - k + 1), 0], {k, 1, nn}], {n, 1,
      nn}]]];
TableForm[
 B = Table[
   Table[If[n >= k, K[[n, k]]/Binomial[n, k], 0], {k, 1, nn}], {n, 1,
    nn}]]
TableForm[B - A]
(*end*)