wolfram-mathematica - 如何在 Mathematica 中用傅立叶变换绘制黎曼 zeta 零谱？

Question

在图 6 中 J. Brian Conrey 的论文“The Riemann Hypothesis”中，有素数定理中误差项的傅立叶变换图。请参阅下图中左侧的图:



在一篇名为 Primes out of Thin Air 的博客文章中由 Chris King 编写，有一个绘制频谱的 Matlab 程序。请参阅帖子开头右侧的情节。可以翻译成 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]

然而，据我所知，这是傅立叶正弦变换的矩阵公式，因此计算成本非常高。我不建议运行它，因为它已经使我的计算机崩溃了一次。

在 Mathematica 中是否有一种方法可以利用快速傅立叶变换来绘制 x 值处的峰值等于黎曼 zeta 零点的虚部的频谱？

我试过命令 FourierDST和 Fourier没有成功。问题似乎是变量yy在代码中同时包含 Sin[t*myspan/xres*yy]和 (y[[Floor[Exp[yy]]]] - Exp[yy])/Exp[yy/2] .

编辑:2012 年 1 月 1 日，我更改了行:
For[yy = 0, yy <= Log[x], 1/stpval++,
进入以下内容:
For[yy = 0, yy/stpval <= Log[x], yy++,
编辑:22.1.2012，来自 Heike 的评论，更改为:
For[yy = 0, yy/stpval <= Log[x], yy++,
进入:
For[yy=0, yy<=Log[x], yy+=1/stpval,

Answer 1

这个怎么办？我使用身份 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]

产生