¿Cómo se traza el espectro cero zeta de Riemann con la transformada de Fourier en Mathematica?

Question

En el documento "The Riemann Hypothesis" de J. Brian Conrey en la figura 6, hay un gráfico de la transformada de Fourier del término de error en el teorema del número primo. Ver la trama a la izquierda en la siguiente imagen:

En un blog llamado Primes out of Thin Air escrito por Chris King, hay un programa de Matlab que traza el espectro. Ver el diagrama a la derecha al comienzo de la publicación. Una traducción en Mathematica es posible:

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] 

Sin embargo, esto es lo que tengo entendido que la formulación de matriz de la transformada de Fourier de seno y por lo tanto es muy costoso de calcular. NO recomiendo ejecutarlo porque ya se bloqueó mi computadora una vez.

¿Existe alguna forma en Mathematica que utilice la Transformada rápida de Fourier para trazar el espectro con picos en valores x iguales a la parte imaginaria de los zeta zeros de Riemann?

He intentado los comandos FourierDST y Fourier sin éxito. El problema parece ser que la variable yy en el código está incluida tanto en Sin[t*myspan/xres*yy] como en (y[[Floor[Exp[yy]]]] - Exp[yy])/Exp[yy/2].

EDIT: 20.1.2012, me cambió la línea:

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

en lo siguiente:

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

EDIT: 01/22/2012, Desde el comentario de Heike, modificado:

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

en:

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

Answer 1

¿Y esto? He reescrito el seno transformar ligeramente usando la identidad 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] 

que produce