Làm thế nào cốt truyện Riemann zeta không phổ với biến đổi Fourier trong Mathematica?

Question

Trong bài báo "Giả thuyết Riemann" của J. Brian Conrey trong hình 6 có một âm mưu của biến đổi Fourier của thuật ngữ lỗi trong định lý số nguyên tố. Xem cốt truyện sang bên trái trong hình dưới đây:

Trong một bài viết trên blog được gọi là Primes out of Thin Air được viết bởi Chris vua có một chương trình Matlab rằng lô quang phổ. Xem cốt truyện ở bên phải ở đầu bài đăng. Một bản dịch sang Mathematica có thể:

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] 

Tuy nhiên, đây là như tôi hiểu nó xây dựng ma trận của sin biến đổi Fourier và do đó nó là rất tốn kém để tính toán. Tôi KHÔNG khuyên bạn nên chạy nó vì nó đã bị hỏng máy tính của tôi một lần.

Có cách nào trong Mathematica sử dụng chuyển đổi Fourier nhanh, để vẽ phổ với các gai tại các giá trị x bằng phần ảo của số không Riemann zeta không?

Tôi đã thử các lệnh FourierDST và Fourier mà không thành công. Vấn đề có vẻ là biến số yy trong mã được bao gồm trong cả hai Sin[t*myspan/xres*yy] và (y[[Floor[Exp[yy]]]] - Exp[yy])/Exp[yy/2].

EDIT: 2012/01/20, tôi đã thay đổi dòng:

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

như sau:

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

EDIT: 2012/01/22, Từ bình luận Heike của, thay đổi:

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

thành:

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

Answer 1

gì về điều này? Tôi đã viết lại sin chuyển hơi bằng cách sử dụng danh tính 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] 

trong đó sản xuất