Introduction to the Bernoulli function
Peter Luschny

A companion to arXiv:2009.06743v2 [math.HO]

Numbering of formulas as in version 2.

Browsable at BernoulliFunctionNotebook.

Source at BernoulliFunctionSource.

Needs Wolfram Language kernel for Jupyter. See instructions at GitHub.

Preliminaries

(* This is used only for some fancy plot legends. *)
(* You will have to do this only once: *)
(* ResourceFunction["MaTeXInstall"][] *)
(* Thereafter you have to execute only: *)
<<MaTeX`
(* MaTeX["x^2"] *)
(* But you not really need these fancy plot legends. *)
(* Comment them out to use this notebbok without MaTex. *)
MaTeX@HoldForm[Integrate[Sin[x], {x, 0, Infinity}]]
(* http://szhorvat.net/pelican/latex-typesetting-in-mathematica.html *)

(* Always execute first. *)

(* 0^0 = 1 *)
Unprotect[Power]; Power[0, 0] = 1; 

(* tau = 2*Pi*I *)
tau := 2*Pi*I;
\[Tau][s_] := tau^(-s) + (-tau)^(-s); 

(* Macro to get a better format *)
Poly[a_] := PolynomialForm[Collect[Expand[Simplify[FunctionExpand[a]]], x], 
            TraditionalOrder -> True];

(* Macro to display a sequence table *)
SeqGrid[T_] := With[{y = Length[T]},
               Grid[{Table[i, {i, 0, y-1}], T}, Frame -> All]];
SeqGrid1[T_] := With[{y = Length[T]},
               Grid[{Table[i, {i, y}], T}, Frame -> All]];             
SeqGrid2[T_] := With[{y = Length[T] - 1},
               Grid[{Table[2 i, {i, 0, y}], T}, Frame -> All]];

Stieljes constants and zeta function

Formula 1

Laurent series of the Riemann zeta function

Series[Zeta[s] - 1/(s - 1), {s, 1, 5}] // TeXForm

$$\gamma -(s-1) \gamma _1+\frac{1}{2} (s-1)^2 \gamma _2-\frac{1}{6} (s-1)^3 \gamma _3+\frac{1}{24} (s-1)^4 \gamma _4-\frac{1}{120} (s-1)^5 \gamma _5+O\left((s-1)^6\right)$$

Formula 2

Stieltjes constants, defined via integral

gamma[s0_, prec_] :=  (* parameter 'prec' for setting the precision *)
Module[{s, z}, {s} = SetPrecision[{s0}, prec $MachinePrecision];
    (-4 Pi /(s + 1)) NIntegrate[Log[1/2 + I z]^(s + 1)/(Exp[-Pi z] + Exp[Pi z])^2,
    {z, 0, Infinity},  WorkingPrecision -> prec $MachinePrecision ]];

Table[gamma[n, 2], {n, 0, 5}]

Table[Re[gamma[n, 2]], {n, 0, 5}]  // MatrixForm  
Table[StieltjesGamma[k], {k, 0, 5}] // N // SeqGrid

Bernoulli constants and Bernoulli function

Formula 3

Bernoulli constants, defined via integral

HoldForm[2 Pi Integrate[Log[1/2 + I z]^s / (Exp[-Pi z] + Exp[Pi z])^2, 
     {z, -Infinity, Infinity}]] // TeXForm

$$2 \pi \int_{-\infty }^{\infty } \frac{\log ^s\left(\frac{1}{2}+i z\right)}{(\exp (-\pi z)+\exp (\pi z))^2} \, dz$$

beta[s0_, prec_] :=  (* parameter 'prec' for setting the precision *)
Module[{s, z}, {s} = SetPrecision[{s0}, prec $MachinePrecision];
    4 Pi NIntegrate[Log[1/2 + I z]^s / (Exp[-Pi z] + Exp[Pi z])^2,
    {z, 0, Infinity},  WorkingPrecision -> prec $MachinePrecision ]];

Table[beta[n, 2], {n, 0, 5}] // MatrixForm

Table[Re[beta[n, 2]], {n, 0, 5}] // MatrixForm // TeXForm

$$\left( \begin{array}{c} 1.0000000000000000000000000000000 \\ -0.57721566490153286060651209008240 \\ 0.14563169096735344972117275174980 \\ 0.029071089578616955453591158105638 \\ -0.008215337681213383464640186171014 \\ -0.011626850327336500287340850887630 \\ \end{array} \right)$$

Memoization of the real part of the Bernoulli constants with single precision. (For quick checking and plotting, for higher precision use the two parameter form.)

beta[n_] := beta[n] = Re[beta[n, 1]];

Plot[beta[s], {s, -0.6, 4}, WorkingPrecision -> 30,
Epilog -> {PointSize[0.01], Red, Point[Table[{k, beta[k]}, {k, 0, 4}]]},  
PlotLegends -> Placed["Bernoulli constants", {Scaled[{0.9, 0.9}], {0.9, 0.9}}]]
(* Export["Fig2BernoulliConstants.eps", %] *)

Formula 4

Bernoulli function, definition.

Bbeta[s_] := Sum[beta[j] s^j / j!, {j, 0, 50}];

Table[Bbeta[n], {n, 0, 10, 2}] // N // SeqGrid2

Table[BernoulliB[n, 1], {n, 0, 10, 2}] // N // SeqGrid2

Formula 5

Bernoulli function using the Riemann zeta representation.

Limit[-n Zeta[1 - n], n -> 0]
B[s_] := -s Zeta[1 - s];
B[0] := 1;

Clear[s]; B[s] // FullSimplify // TeXForm

1

$$-s \zeta (1-s)$$

From now onwards we use Mathematica's efficient implementation of the Zeta function for the Bernoulli function.

Plot[B[s], {s, 0, 12}, PlotRange -> {-1/3, 1},
Epilog -> {PointSize[0.01], Red, Point[Table[{k,BernoulliB[k,1]}, {k, 0, 12}]]}]

Bernoulli function on the critical line and the zeta zeros which are the same as the Bernoulli zeros.

ReImPlot[B[1/2 + I t], {t, 20, 52},
Epilog -> {PointSize[0.01], Red, Point[Table[{Im[ZetaZero[k]], 0}, {k, 11}]]}]
(* Export["Fig10BernoullisZetaZeros.eps", %] *)

Table[Im[ZetaZero[n]], {n, 1, 11}] // N // SeqGrid1

Bernoulli function rises up at Riemann's critical line.

ComplexPlot3D[B[s], {s, -60 - 40 I, 1/2 + 40 I}, 
WorkingPrecision -> 40,
AxesLabel -> {"Re(s)", "Im(s)"},
ImageSize -> Large,
ViewPoint -> {2.6, -0.9, 0.6} ]
(* Export["Fig27BernoulliTsunami.pdf", %] *)