Come traccia lo spettro zero di Riemann zeta con la trasformata di Fourier in Mathematica?

Question

Nel documento "The Riemann" di J. Brian Conrey in figura 6 è un grafico della trasformata di Fourier del termine di errore nel teorema dei numeri primi. Vedere la trama verso sinistra nell'immagine qui sotto:

In un post sul blog chiamato Primes out of Thin Air scritto da Chris King c'è un programma Matlab che traccia lo spettro. Guarda la trama a destra all'inizio del post. Una traduzione in Mathematica è possibile:

Mathematica:

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], yy+=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] 

Tuttavia, questo è quanto mi risulta la formulazione della matrice del seno trasformata di Fourier ed è quindi molto costoso da calcolare. NON consiglio di eseguirlo perché si è già bloccato il mio computer una volta.

C'è un modo in Mathematica che utilizza la trasformata di Fourier veloce, per tracciare lo spettro con picchi ai valori x uguale alla parte immaginaria di zeri zeri di Riemann?

Ho provato i comandi FourierDST e Fourier senza successo. Il problema sembra essere che la variabile yy nel codice sia inclusa in entrambi Sin[t*myspan/xres*yy] e (y[[Floor[Exp[yy]]]] - Exp[yy])/Exp[yy/2].

EDIT: 20.1.2012, ho cambiato la linea:

For[yy = 0, yy <= Log[x], 1/stpval++,

nelle seguenti:

For[yy = 0, yy/stpval <= Log[x], yy++,

EDIT: 2012/01/22, Da commento di Heike, modificata:

For[yy = 0, yy/stpval <= Log[x], yy++,

in:

For[yy=0, yy<=Log[x], yy+=1/stpval,

Answer 1

Che dire di questo? Ho riscritto il seno trasformare un po 'utilizzando l'identità Exp[a Log[x]]==x^a

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] 

che produce