Как рассчитать спектр нулевого нуля Римана с преобразованием Фурье в Mathematica?

В статье "Гипотеза Римана" Дж. Брайана Конри на рисунке 6 приведен график преобразования Фурье члена ошибки в теореме о простых числах. Смотрите рисунок слева на изображении ниже:

В сообщении в блоге, которое называется Primes из Thin Air, написанном Chris King, есть программа Matlab, которая отображает спектр. См. График справа в начале сообщения. Возможно перевод в Mathematica:

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, использующем быстрое преобразование Фурье, для построения спектра с шипами при значениях х, равных мнимой части дзета-нулей Римана?

Я пробовал команды FourierDST и Fourier без успеха. Проблема заключается в том, что переменная yy в коде включена как в Sin[t*myspan/xres*yy], так и (y[[Floor[Exp[yy]]]] - Exp[yy])/Exp[yy/2].

EDIT: 20.1.2012, я изменил строку:

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

в следующее:

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

EDIT: 22.1.2012, Из комментария Heike, изменено:

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

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

Ответ 1