The Mobius function Blog

Fourier transform of PNT error term plot

Mathematica 8 code from Heike’s answer at stackoverflow:

Clear[f]
scale = 1000000;
f = ConstantArray[0, scale];
f[[1]] = N@MangoldtLambda[1];
Monitor[Do[
f[[i]] = N@MangoldtLambda[i] + f[[i - 1]], {i, 2, scale}], i]
xres = .002;
xlist = Exp[Range[0, Log[scale], xres]];
tmax = 60;
tres = .015;
Monitor[errList =
Table[(xlist^(-1/2 + I t).(f[[Floor[xlist]]] - xlist)), {t,
Range[0, 60, tres]}];, t]
ListLinePlot[Im[errList]/Length[xlist], DataRange -> {0, 60},
PlotRange -> {-.09, .02}, Frame -> True, Axes -> False]

Clear[f]
scale = 1000000;
f = ConstantArray[0, scale];
f[[1]] = N@MoebiusMu[1];
Monitor[Do[f[[i]] = N@MoebiusMu[i] + f[[i - 1]], {i, 2, scale}], i]
xres = .002;
xlist = Exp[Range[0, Log[scale], xres]];
tmax = 60;
tres = .015;
Monitor[errList =
Table[(xlist^(-1/2 + I t).(f[[Floor[xlist]]] - 0*xlist)), {t,
Range[0, 60, tres]}];, t]
ListLinePlot[Im[errList]/Length[xlist], DataRange -> {0, 60},
PlotRange -> {-.02, .10}, Frame -> True, Axes -> False]

Clear[n, d];
Clear[f]
scale = 1000000;
f = ConstantArray[0, scale];
f[[1]] = Sum[d*MoebiusMu@d, {d, Divisors[1]}];
Monitor[Do[
f[[i]] = Sum[d*MoebiusMu@d, {d, Divisors[i]}] + f[[i - 1]], {i, 2,
scale}], i]
xres = .002;
xlist = Exp[Range[0, Log[scale], xres]];
tmax = 60;
tres = .015;
Monitor[errList =
Table[(xlist^(-1/2 + I t).(f[[Floor[xlist]]] - 0*xlist)), {t,
Range[0, 60, tres]}];, t]
ListLinePlot[Im[errList]/Length[xlist]/100000, DataRange -> {0, 60},
PlotRange -> {-.05, .05}, Frame -> True, Axes -> False]

Excel cell formula, European dot-comma style:

=IF(OR(ROW()=1;COLUMN()=1);1;IF(ROW()=COLUMN();-(SUM(INDIRECT(ADDRESS(1;COLUMN())&":"&ADDRESS(COLUMN()-1;COLUMN())))+SUM(INDIRECT(ADDRESS(1;COLUMN())&":"&ADDRESS(COLUMN()-1;COLUMN())))/(COLUMN()-1));IF(ROW()>COLUMN();INDIRECT(ADDRESS(ROW()-COLUMN();COLUMN()));INDIRECT(ADDRESS(COLUMN();COLUMN()-ROW())))))

Excel cell formula, American dot-comma style:
=IF(OR(ROW()=1,COLUMN()=1),1,IF(ROW()=COLUMN(),-(SUM(INDIRECT(ADDRESS(1,COLUMN())&":"&ADDRESS(COLUMN()-1,COLUMN())))+SUM(INDIRECT(ADDRESS(1,COLUMN())&":"&ADDRESS(COLUMN()-1,COLUMN())))/(COLUMN()-1)),IF(ROW()>COLUMN(),INDIRECT(ADDRESS(ROW()-COLUMN(),COLUMN())),INDIRECT(ADDRESS(COLUMN(),COLUMN()-ROW())))))

Which outputs in the main diagonal:
1, -2, -3, 0, -5, 6, -7, 0, 0, 10, -11,…

I have written this code just to remember it. DO NOT run it! It eats up all of your computers memory.

Mathematica 8

scale = 10^6;
start = 1;
fin = 50;
its = 490;
xres = 600;
y = N[Accumulate[Table[MangoldtLambda[i], {i, 1, scale}]], 10];
x = scale;
a = 1;
myspan = 800;
xres = 4000;
xx = N[Range[a, myspan, (myspan - a)/(xres - 1)]];
stpval = 10^4;
F = Range[1, xres]*0;

For[t = 1, t <= xres, t++,
For[yy = 0, yy <= Log[x], 1/stpval++,
F[[t]] =
F[[t]] +
Sin[t*myspan/xres*yy]*(y[[Floor[Exp[yy]]]] - Exp[yy])/Exp[yy/2];
]
]
F = F/Log[x];
ListLinePlot[F]

Mathematica 8

Clear[a, b, c, d, e, f];
a = 2^N[1/2 + I*Im[ZetaZero[1]]/2, 50]*Cos[2*Pi*1/2];
b = 2^N[ZetaZero[2], 50]*Cos[2*Pi*2/2];
c = a + b;
d = 2^N[ZetaZero[1], 50]*Cos[2*Pi*1/2];
e = 2^N[ZetaZero[2], 50]*Cos[2*Pi*2/2];
f = d + e;
g = Re[c + f]
N[Im[ZetaZero[1]], 50]

Which outputs:
-0.134729308890067286283770094746871131437286485045
14.134725141734693790457251983562470270784257115699

Calculation in Wolfram Alpha

Mathematica 8

Clear[nn, t, n, k, a, b];
nn = 15;
t[n_, 1] = 1;
t[1, k_] = 1;
t[n_, k_] :=
t[n, k] =
If[n >= k, -Sum[t[n - i, k], {i, 1, k - 1}], -Sum[
t[k - i, n], {i, 1, n - 1}]];
MatrixForm[Table[Table[t[n, k], {k, 1, nn}], {n, 1, nn}]];
MatrixForm[
Table[Table[t[n, k]*Exp[I*2*Pi*k/n], {k, 1, n}], {n, 1, nn}]]
MatrixForm[
Chop[Table[
N[Total[Table[t[n, k]*Exp[I*2*Pi*k/n], {k, 1, n}]]], {n, 1, nn}]]]

Consider the infinite matrix starting:

Notice the minus sign in front of the matrix. I was too lazy to change the sign in front of every element in the matrix. List the largest eigenvalues from each matrix:

{-1., 1.41421, 2.65544, 3.43931, 4.77106, 5.24392, 6.84437, 7.15539, \
7.47476, 7.57341, 10.9223, 11.096, 12.9021, 13.0453, 13.259, 13.4055, \
16.9724, 17.0824, 18.9443, 19.0552, 19.2282, 19.307, 22.9972, \
23.0759, 23.1576, 23.2173, 23.2976, 23.3972, 29.0103, 29.0407, \
30.963, 31.0104, 31.1008, 31.1505, 31.268, 31.34, 37.0284, 37.0658, \
37.1289, 37.174, 41.029, 41.0503, 42.9921, 43.0326, 43.0807, 43.1149, \
46.996, 47.0293, 47.0619, 47.1025, 47.1582, 47.2011, 53.0192, \
53.0497, 53.1076, 53.1419, 53.1893, 53.2117, 59.0477, 59.0681, \
61.0248, 61.0474, 61.0812, 61.1071}

As rounded to nearest integer the values are:
{-1, 1, 3, 3, 5, 5, 7, 7, 7, 8, 11, 11, 13, 13, 13, 13, 17, 17, 19, \
19, 19, 19, 23, 23, 23, 23, 23, 23, 29, 29, 31, 31, 31, 31, 31, 31, \
37, 37, 37, 37, 41, 41, 43, 43, 43, 43, 47, 47, 47, 47, 47, 47, 53, \
53, 53, 53, 53, 53, 59, 59, 61, 61, 61, 61}

Which is close to the previous prime number.

Plotting them we get:

A useful lower bound for the largest eigenvalue would say something about the distribution of prime numbers.

Or as a ListPlot:

Mathematica 8 code for plot of largest eigenvalues:

Clear[a, t, n, k, i, j]
t[n_, 1] = -1;
t[1, k_] = -1;
t[n_, k_] :=
t[n, k] =
If[n >= k, -Sum[t[n - i, k], {i, 1, k - 1}], -Sum[
t[k - i, n], {i, 1, n - 1}]];
nn = 64;
a = Range[1, nn]*0;
b = Range[1, nn]*0;
Do[m = Table[Table[t[n, k], {k, 1, j}], {n, 1, j}];
a[[j]] = Eigenvalues[m], {j, 1, nn}]

Round[Table[Max[N[a[[i]]]], {i, 1, nn}], 0.00001]
ListLinePlot[Table[Max[N[a[[i]]]], {i, 1, nn}]]
ListLinePlot[
Flatten[{0, Differences[Table[Max[N[a[[i]]]], {i, 1, nn}]]}],
Ticks -> {Range[nn]}]

Excel cell formulas for the matrix:

European:
=if(or(row()=1; column()=1); 1; if(row()>=column(); -sum(indirect(address(row()-column()+1; column(); 4)&":"&address(row()-1; column(); 4); 4)); -sum(indirect(address(column()-row()+1; row(); 4)&":"&address(column()-1; row(); 4); 4))))

American:
=if(or(row()=1, column()=1), 1, if(row()>=column(), -sum(indirect(address(row()-column()+1, column(), 4)&":"&address(row()-1, column(), 4), 4)), -sum(indirect(address(column()-row()+1, row(), 4)&":"&address(column()-1, row(), 4), 4))))

Edit 20.3.2012: Added a picture of the eigenvalues as radicals:

The first of these plots gives a rough picture of the von Mangoldt function.
I first learned it from J. Brian Conreys paper:

The Riemann Hypothesis, J.Brian Conrey

and later understood it from Edwin Chen’s blog:

Prime Numbers and the Riemann Zeta Function, Edwin Chen

But I have simplified it to let “x” be raised to Riemann zeta zero without the other terms
as in Chen’s blog post


Clear[n, k, A, a, B, b, kk, ii, i]
kk = 50;(*number of Zeta zeros used*)
ii = 10;(*reciprocal of x-axis spacing*)
n = Table[i, {i, 0, 42, 1/ii}];
A = Table[Re[n^N[ZetaZero[k], 12]], {k, 1, kk}];
a = Table[{-i/ii, Plus @@ A[[All, i]]}, {i, 1, Length[n]}];
a[[ii + 1]] = {-(ii + 1)/ii, 0};
ListLinePlot[-a]


Check out Heike’s answer on stack overflow:
Riemann zeta zero spectrum

Mathematica 8:
Clear[a, t, n, k, i, j]
t[n_, 1] = 1;
t[1, k_] = 1;
t[n_, k_] :=
t[n, k] =
If[n >= k, -Sum[t[n - i, k], {i, 1, k - 1}], -Sum[
t[k - i, n], {i, 1, n - 1}]];
nn = 64;
a = Range[1, nn]*0;
b = Range[1, nn]*0;
Do[
m = Table[Table[t[n, k], {k, 1, j}], {n, 1, j}];
a[[j]] = Eigenvalues[m],
{j, 1, nn}]

Round[Table[-Min[N[a[[i]]]], {i, 1, nn}], 0.00001]
(*ListLinePlot[-Table[Min[N[a[[i]]]],{i,1,nn}]]*)
ListLinePlot[
Flatten[{0, -Differences[Table[Min[N[a[[i]]]], {i, 1, nn}]]}],
Ticks -> {Range[nn]}]
ListLinePlot[
Flatten[{0, Differences[Table[Max[Abs[N[a[[i]]]]], {i, 1, nn}]]}],
Ticks -> {Range[nn]}]

2.6173553781947777777935357228301908029227443780483

N[Pi^(Im[ZetaZero[2]]/Im[ZetaZero[3]]), 50]

N[Pi^(Im[ZetaZero[2]]/Im[ZetaZero[3]]), 50]

Wolfram Alpha