[-*- Mode: lispcc -*-]
[-----grammar
(gz start (prog))
(gz prog (global-star
::pr("[[global-star('','\n','\n')]]")))
(gz global (literal)(quotbrak))
(gz literal (id )
(astring
::pr((::c "my_astring->print_string();")))
)
(gz quotbrak (f id f modify-star j morestuff j
::pr("<[[id]] [[modify-star('',' ','')]]>[[morestuff]]"))
( f id f modify-star j :x j
::pr("<[[id]] [[modify-star('',' ','')]]/>"))
(f :p morestuff j
::pr("([[morestuff]])"))
(f :quote morestuff j
::pr("\x22[[morestuff]]\x22"))
(f :d id morestuff j
::pr("[[morestuff]]
\n"))
(f :s id morestuff j
::pr("[[morestuff]]"))
(f :c global-star j
::pr("[[global-star('','','')]]"))
)
(gz modify (f id qastring j ::pr ("[[id]]=[[qastring]]")))
(gz qastring (astring ::pr("\x22[[astring]]\x22")))
(gz morestuff (global-star
::pr("[[global-star('',' ','')]]")))
]
[""]
""
[(setq ss-verify-run nil)]
(html()
(head ()
(title () Tabulating values of the Riemann-Siegel Z function along the critical line)
[""]
(meta ((http-equiv "Content-Type")
(content "text/html; charset=iso-8859-1")
) :x)
(link ((rel "contents")(href "#TOC")):x)
(link ((rel "chapter")(href "#Introduction")(title "Introduction")):x)
(link ((rel "chapter")(href "#Evaluating")(title "Riemann-Siegel formula")):x)
(link ((rel "chapter")(href "#Theoretical")(title "Verifying the remainder terms")):x)
(link ((rel "chapter")(href "#Chebyshev")(title "Chebyshev polynomials for the remainder terms")):x)
(link ((rel "chapter")(href "#Implementation")(title "Implementation details")):x)
(link ((rel "chapter")(href "#Parallel")(title "Parallel Execution")):x)
(link ((rel "chapter")(href "#Tabulations")(title "Tabulations")):x)
(link ((rel "chapter")(href "#Analyses")(title "Analyses")):x)
(link ((rel "chapter")(href "#Source")(title "Source Code")):x)
(link ((rel "chapter")(href "#References")(title "References")):x)
(link ((rel "author") (href "http://web.mit.edu/kenta/www/")
(title "Ken Takusagawa")):x)
(link ((rel "start") (href "#")):x)
(link ((rel "parent")(href "./")):x)
(style ((type "text/css"))
.par { margin: 1.33em 0}
.center { text-align: center}
table {margin-left: auto;
margin-right:auto}
.indentFirst { margin: 1.33em 0 0 3em}
.indent { margin: 0 0 0 3em}
.indentMore { margin: 0 0 0 5em }
)
)
[body { background: silver}]
[
c0.001.ppm.gif
c1.001.ppm.gif
c2.001.ppm.gif
c3.001.ppm.gif
c4.001.ppm.gif
c5.001.ppm.gif
four1.gif
four2.gif
fourhi.gif
vcss.gif
]
[
asympt.png
brec.png
calligraphic-I.png
cform.png
gue1.png (width "648")(height "477")
gue2.png
product001.png
siegel1.png
siegel2.png
thetae.png
zetaZ.png
]
(body ()
(h1() Tabulating values of the Riemann-Siegel Z function along the
critical line)
(h2() Ken Takusagawa)
[(div() there once was a (a ((href "http://www.mit.edu/"))cat) named mike.)]
(div((id "Abstract"))
(h2() Abstract)
This project makes three main contributions:
(ul()
(li() Tabulation
of 2,354,000,000 values of the Riemann-Siegel Z
function in the domain 2.7 "<" (em()t) "<" 29,143,636.6, and
tabulation of
(a ((href "other-files/sample_31801772196387.flac"))5,000,000
values)
of the Riemann-Siegel Z function
in neighborhood of the trillionth (:p (:c 10(sup()12)th)) zero.)
(li() Tabulation of
(a ((href "other-files/47-coefficients.txt" ))
47 terms of the remainder series) of the
Riemann-Siegel formula.)
(li() Tabulation of the
(a ((href "other-files/coff-all.h"))
coefficients of Chebyshev polynomials) for
the first 6 terms of the remainder series of the Riemann-Siegel formula.)
))
[>>> samnumber (x)
ans = 267653395646.999
octave:10>
octave:10> samnumber(x+ 5000000)
ans = 267653436311.838
x = 31801772196387
]
(div((id "TOC"))
(ol()
(li() (a ((href "#Introduction"))Introduction))
(li() (a ((href "#Evaluating")) Riemann-Siegel formula))
(li() (a ((href "#Theoretical")) Verifying the remainder terms))
(li() (a ((href "#Chebyshev")) Chebyshev polynomials for the remainder terms))
(li() (a ((href "#Implementation")) Implementation details))
(li() (a ((href "#Parallel")) Parallel Execution))
(li() (a ((href "#Tabulations")) Tabulations))
(li() (a ((href "#Analyses")) Analyses))
(li() (a ((href "#Source"))Source Code ))
(li() (a ((href "#References"))References ))
[(li() (a ((href "#")) ))]
))
(div((id "Introduction"))
(h2()Introduction)
(:d par
Riemann's zeta function is defined as the analytic continuation of
(:d center
(img ((src "other-files/zetadef2m.png")
(alt "\zeta(s)=\sum_{n=1}^{\infty}\frac{1}{s^n}")
(width "159")(height "70")
):x))
to the entire complex plane.
(:p The series itself only converges for (:c Re
(:p s)) ">"
1.) The (em()Riemann conjecture) states that all non-trivial
zeros of the zeta function occur along the line (:c Re(:p s)) = 0.5.
This line is known as the (em()critical line.)
)
(:d par
The values (:p the real and imaginary parts) of the zeta
function evaluated along the critical line oscillate around zero, as the
plot below illustrates.)
(table((border "1"))
(tr()(td()(img ((src "other-files/zetaplot.png")
(alt "plot of zeta")
(width "366")(height "327")
):x))))
(:d par
I was inspired from the visual appearance of these plots to try to
produce a (em()sound file) of the zeta function along the critical
line. A similar attempt was made by
(a ((href
"http://www.maths.ex.ac.uk/~mwatkins/zeta/munafo-zetasound.htm"))
Robert Munafo).)
(:d par
Therefore, the (em()real) motivation for this project is
(em()artistic). I just wanted to hear what it sounded like,
(:p and to take advantage of a computational resource that I would
not likely have again any time in the near future).
I did not want discover any deep mathematical
revelations, although the tabulations may be useful to other
mathematicians for that endeavor.
There is no attempt to solve any particularly
difficult or interesting parallel computer programming problem;
the parallelization implemented was of the
(:quote embarrassingly) variety and relatively trivial.)
(:d par
To borrow a quote from George Mallory, I performed this
calculation of the Z
function
(:quote because it's there.)
)
(div ()
Useful links
(ul()
(li()(a ((href "http://www.math.ubc.ca/~pugh/RiemannZeta/RiemannZetaLong.html"))Glen Pugh))
(li()(a ((href "http://www.dtc.umn.edu/~odlyzko/"))Andrew Odlyzko))
))
)
(div ((id "Evaluating"))
(h2() Riemann-Siegel formula)
(:d par
The Riemann-Siegel formula is an approximation formula for the
Zeta function specifically on the critical line. Evaluation of
the formula
requires (:c O(:p (:c (em()t)(sup()1/2)))) operations to evaluate a single value
of (:c Zeta(:p 0.5 + (em()t) *i)). This is in contrast to the faster Odlyzko -
"Schönhage" method, which requires (:c O(:p (:c (em()t)(sup() "ε"))))
operations (:p where "ε" can be made arbitrarily small at
the cost of a higher space requirement).
However, the Odlyzko -
"Schönhage" is very complicated, so I chose not to implement
it.
)
(:d par The Riemann-Siegel formula
(:p excerpted from
(a ((href "other-files/masters.thesis.pdf"))(:p Pugh
1998)))
for calculating (:c Zeta (:p 0.5 + (em()t) * i)) for (em()t) ">" 0 is:)
(:d center
(img ((src "other-files/siegel1.png")
(alt "Riemann-Siegel formula, part 1")
(longdesc "other-files/masters.thesis.pdf")
(width "811")(height "313")
):x))
(:d center
(img ((src "other-files/siegel2.png")
(alt "Riemann-Siegel formula, part 2")
(width "762")(height "465")
):x))
(:d par Note that the calligraphic I
(:p (img ((src
"other-files/calligraphic-I.png")(alt "I")
(width "14")(height "15")):x))
function signifies (:quote take the imaginary part).
The capital Pi (:p "Π") function is (:c "Π"(:p s)) = (:c Gamma
(:p s+1)). [The C(sub() n) functions are functions of (em() p).]
The superscripts above the Psi (:p "Ψ") functions signify
repeated differentiation with respect to (em()p).)
(:d par The (:c C(sub() n)) functions are undefined at (em()p)= 0.25 and
(em()p)= 0.75. The values at those points are the limiting values
(:p from either side).)
(:d par And finally, (:c Zeta (:p 0.5 + (em()t) * i)) is calculated as:)
(:d center(img ((src "other-files/zetaZ.png")
(alt "Riemann-Siegel formula, part 3")
(width "361")(height "109")
):x))
(:d par Note that the sum in the (:c Z(:p (em()t))) function is the
only computationally expensive part of the calculation. After
(:c Z(:p t)) is known, (:c Zeta (:p 0.5 + (em()t) * i)) can be calculated easily.
Therefore, the remainder of this project deals with tabulating
values of the Z function, and not the Zeta function.
(:p The tabulated values of Z require half as much space to store
because
they are not complex.) Also, the zeros of the Z function
coincide with the zeros of the Zeta function along the critical line.
)
)
(div((id "Theoretical"))
(h2() Verifying the remainder terms)
(:d par
The first part of this project was to verify (:p and extend) the
remainder terms (:p i.e., the (:c C(sub() n)) functions) of the
Riemann-Siegel formula. Both references to the formula I found on
the web, namely (:p Pugh 1998) at
(a ((href
"http://mathworld.wolfram.com/Riemann-SiegelFunctions.html"))
Mathworld) both referenced Edwards's 1974 book for the remainder
terms. It was conceivable that an arithmetic or typographical
error could have been propagated through the documents.)
(:d par Note that I am verifying the low-level arithmetic and
algebra of the derivation of the remainder terms, and not the
high-level complex analysis upon which the Riemann-Siegel formula
is based.)
(:d par My verification followed the derivation described in
(:p Pugh 1998)
sections 3.8-3.9. Note that there is one typographical error and
one confusing point in his derivation,
(a ((href "other-files/pugh-email.txt"))
see my e-mail conversation with him.) The formulae below have
the typographical error corrected.
)
(:d par The first step is to calculate the sequence
(:c b(sub()n))
from the following recurrence:)
(:d center
(img ((src "other-files/brec.png")
(alt "recurrence for b")(width "250")(height "74")
):x))
(:d par where (:c b(sub()-2)) = 0, (:c b(sub()-1)) = 0, and
(:c b(sub()0)) = 1.)
(:d par The second step is to calculate the sequence
c(sub()n) from the following formula:)
(:d center
(img ((src "other-files/cform.png")
(alt "formula for c")
(width "264")(height "78")
):x))
(:d par My Haskell program
(a ((href
"other-files/zetaR.el.hs"))zetaR.el.hs) calculates these
symbolic expressions and produces output which can be fed into
Matlab. (:p The programming language (a ((href
"http://www.haskell.org/"))
Haskell,)
with its built in support for arbitrary
precision arithmetic on rationals and (:quote lazy) infinite
lists (:p used to implement infinite series) was quite
useful for this part of the project.) I calculated the
first 201 terms (:p up to (em()n) = 200) of each sequence.
)
(:d par The remainder of the calculation was done with Matlab's
symbolic toolbox, and is documented in the Matlab diary file
(a ((href
"other-files/47-coefficients.txt"))47-coefficients.txt).)
(:d par The following product is formed:)
(:d center
(img ((src "other-files/product001.png")
(alt
"e^{-i(\frac{1}{48t}+\frac{7}{5760t^3}+\cdots)}\sum_{n=0}^{K}b_nc_n")
(width "331")(height "86")
):x))
(:d par into which we substitute
(em()t) = (:c 2 "π"/"ω"(sup()2)), and the
notation (:quote (:c e(sup()X))) denotes the Taylor expansion of the
exponential function about zero.
)
(:d par The series inside the parenthesis are the negative order
terms of the series expansion of the Theta function:)
(:d center
(img ((src "other-files/thetae.png")
(alt "Formula Theta")
(width "320")(height "69")
):x))
(:d center
(img ((src "other-files/asympt.png")
(alt "Asypmtotic Series for Theta")(width "405")(height "64")
):x))
(:d par The series follows from the Stirling approximation of the
Gamma function. Coaxing this series out of Matlab was somewhat
involved, the line in the
(a ((href
"other-files/47-coefficients.txt"))diary) which produces 150
terms of the series is (code()
"gg=sym(maple('asympt',rr,x,150))"))
(:d par Finally, after forming the product, we can collect by
powers of "ω",
(:p The Matlab command was
(code() "maple('coeff',infisubs*psum,omega,i)")),
yielding the n-th remainder term as the
coefficient to (:c "ω"(sup()n)). I ran this calculation until
Maple broke: For the 47th term, it died with an error message
complaining about (:quote integer too large in context).)
(:d par In the
(a ((href "other-files/47-coefficients.txt"))diary) file, the
notation (code()Pn), e.g. (code()P7), denotes the (em()n)-th
derivative of the Psi (:p "Ψ") function.
The first few terms beyond those given in Pugh's thesis
are:)
(ul()
(li()
(:c C(sub()5)) =
(code()
-1/978447237120/pi^10*P15 - 7/849346560/pi^8*P11 - 5/3072/pi^4*P3
- 901/82575360 /pi^6*P7 ))
(li()(:c C(sub()6)) =
(code()
1/563585608581120/pi^12*P18 + 18889/237817036800/pi^8*P10 +
17/652298158080/pi^10*P14 + 367/7864320/pi^6*P6 + 5/2048/pi^4*P2
))
(li()(:c C(sub()7)) =
(code()
-5/2048/pi^4*P1 - 407/2621440/pi^6*P5 - 6649/11890851840/pi^8*P9 -
2131/5707608883200/pi^10*P13 - 1/15655155793920/pi^12*P17 -
1/378729528966512640/pi^14*P21))
(li()(:c C(sub()8)) =
(code()
1/290864278246281707520/pi^16*P24 +
11153/8766887244595200/pi^12*P16 + 26405/8455716864/pi^8*P8 +
427/1048576/pi^6*P4 + 41/32768/pi^4*P0 +
88651/22830435532800/pi^10*P12 + 23/180347394745958400/pi^14*P20
))
)
(:d par
More terms are given at the end of the
(a ((href "other-files/47-coefficients.txt"))diary file).
(:c C(sub()5)) agrees with a formula given at
(a ((href
"http://mathworld.wolfram.com/Riemann-SiegelFunctions.html"))
Mathworld). (:c C(sub()6)) through (:c C(sub()46)) have not been
published before to the best of my knowledge.)
(:d par This project would have benefitted from a parallelized
symbolic computation package.
The timings of the symbolic
calculations are included in the
(a ((href
"other-files/47-coefficients.txt"))diary) file.
Some of the commands listed in the
(a ((href
"other-files/47-coefficients.txt"))diary)
took longer than 10 minutes. The most time-consuming calculation
(:p not shown in diary file) was the dot product
(:c "Σ"b(sub()n)c(sub()n),) which took over 30 minutes.)
(:d par The directory
(a ((href "other-files/remainder_terms/"))remainder_terms/) contains Matlab
scripts used in the computation.)
[(p() Extra terms calculated, there is a need for parallel maple)]
)
(div((id "Chebyshev"))
(h2() Chebyshev polynomials for the remainder terms)
(:d par It would be very unwise to implement the remainder terms exactly as
the expressions which result from the symbolic computations above.
The high derivatives of the Psi function yield extremely long
expressions, which result in numerical instability, especially near
the singular points (em()p)= 0.25 and
(em()p)= 0.75.
)
(:d par The wise thing to do is to symbolically compute some sort of series
expansion of the remainder terms. Pugh, in his master's thesis,
chose a Taylor expansion (:p around (em()p) = 0.5). This project
does slightly better by computing a Chebyshev polynomial for each
of the remainder terms.)
(:d par (:p Although, since very little of the computation time is
spent on evaluating the remainder terms, the computation
time saved by using a Chebyshev polynomial instead of a
Taylor polynomial is negligible.))
(:d par The algorithm described in (em()Numerical Recipes) was used
to generate the Chebyshev coefficients. The calculations were
carried out with 200 digits of precision using the (code()vpa)
command in Matlab. These computations were extremely time
consuming as well, taking on the order of 30 minutes per remainder
term. Fortunately, the task is embarrassingly parallel, and was
carried out on 6 computers, one for each term calculated. With
these terms calculated here, hopefully no one will have to
calculate them again.)
(:d par The Chebyshev coefficients for the first 6 remainder terms
(:p namely (:c C(sub()0)) through (:c C(sub()5))) are given in
(a ((href "other-files/coff-all.h"))coff-all.h). There are 101
coefficients for each term. For my computation, I truncated the
polynomial to
20 terms.)
(:d par Plots of the first six remainder terms look as follows:)
(table((border "1"))
(tr()(td()(img ((src "other-files/c0.001.ppm.gif")
(alt "C0")(width "586")(height "500")
):x))))
(table((border "1"))
(tr()(td()(img ((src "other-files/c1.001.ppm.gif")
(alt "C1")(width "618")(height "496")
):x))))
(table((border "1"))
(tr()(td()(img ((src "other-files/c2.001.ppm.gif")
(alt "C2")(width "595")(height "496")
):x))))
(table((border "1"))
(tr()(td()(img ((src "other-files/c3.001.ppm.gif")
(alt "C3")(width "590")(height "496")
):x))))
(table((border "1"))
(tr()(td()(img ((src "other-files/c4.001.ppm.gif")
(alt "C4")(width "583")(height "491")
):x))))
(table((border "1"))
(tr()(td()
(img ((src "other-files/c5.001.ppm.gif")
(alt "C5")(width "590")(height "496")
):x))))
(:d par The directory
(a ((href "other-files/cheby/")) cheby/) contains Matlab
scripts used in this computation.)
)
(div((id "Implementation"))
(h2() Implementation details)
(:d par
We now discuss a few of the implementation details of the ideas
presented above.)
(h3() Verifying correctness)
(:d par
For low values of (em()t), the implementation was verified by
comparing with Matlab's (code()zeta) function. Matlab's
(code()zeta) function uses an algorithm which requires
(:c O(:p (em()t))) operations per point, so can only be used for low
values of t.)
(:d par
For high values of (em()t), I evaluated Z at
the (a ((href "other-files/zeros3.gz"))10,000 zeros)
in the neighborhood of the trillionth zero
previously computed by A. Odlyzko. If my program is implemented
correctly, all of these evaluations should come out zero, which
they did.
The maximum deviation from zero was 1.4235e-05, and the
root-mean-square (:p RMS) deviation over all 10,000 zeros was
9.7010e-07.)
(h3((id "scaling")) Scaling horizontally)
(:d par
The density of the zeros of the Z function
(:p and consequently the Zeta function) increases as (:c O
(:p log (em()t)).) This (:quote sort of) follows from the
(:quote (:c (em()t)/2) log (:p (:c(em()t)/(:p 2 pi)))) leading term of the
series expansion of the Theta function above.
)
(:d par
Recall that the purpose of this project is to produce a sound file.
The effect of the increasing density of the sound file is that the
frequency of the sound will increase logarithmically.
Unfortunately, this will make the sound file annoying or painful
to listen, so we must counteract this effect.)
(:d par
The countermeasure I implemented was to scale the (em()t) axis by
1/log(em()t). Specifically:)
(ul()
(li() Let N be the sample number, where N= 0,1,2,...)
(li() Let (em()y) = 3+(:p N/4) )
(li() Let (em()t) = (em()y)/log (em()y) )
(li() Then, the N-th sample is the value of (:c Z(:p (em()t))).))
(:d par
These parameters were chosen so that the sound file is not too
annoying to listen, when played at 44.1 kHz.
The frequency does slowly change over the
course of the sound file. This is because (em()x) / log(em()x) is
not an exact inverse of (em()x) * log(em()x).
)
(h3() Scaling Vertically)
(:d par The sound file stores samples as 16-bit signed integers, so
we must scale the value returned by the Z function into the range
-32768 ... 32767.)
(:d par The maximum absolute value of the Z function encountered in
the first 2,354,000,000 samples was 74.2677, which occurred at
(em()t) = 15457423.7117588, or sample number 1207244046.
Therefore, the scaling constant was chosen to be 256.
Note that the reciprocal of the scaling constant can be
exactly represented in binary, so the value of Z can be
reconstructed exactly (:p up to the limit of roundoff error).)
(:d par The maximum absolute value of the Z function encountered in
the neighborhood of the trillionth zero was -134.2658, which occurred at
(em()t) = 267653412338.13435665, or sample number 31801774248668.
Therefore, the scaling constant was chosen to be 128.)
(:d par It is interesting to examine these peaks with
(a ((href
"http://www.math.ubc.ca/~pugh/RiemannZeta/RiemannZetaLong.html"))
Glen Pugh's Zeta plotter).)
(:d par Unfortunately, these peaks were discovered by a two pass
algorithm: First compute all the values and determine the min and
max, then compute all the values again scaling by the right
constant.)
(:d par Extreme values of Z occur quite rarely, so the sound file is
rather quiet as a result.)
(h3()Key optimization)
(:d par Almost all of the computation time of my program is spent in
an inner loop which computes the sum of the Riemann-Siegel
formula. This inner loop has two key optimizations. The relevant
lines of code are:)
(:d indentFirst (code()"for (int n = 1; (n <= N); ++n) {"))
(:d indentMore
(code ()
"sumZ += (cosl ((theta_t - (t * lookup[n]._logl))) / lookup[n]._sqrtl)"))
(:d indent (code()"}"))
(ol()
(li() The values of (code()log) and (code()sqrt) are computed via lookup
table. )
(li() The lookup table is stored as an array of records
(:p as opposed to parallel arrays), so that entries for
(code()_logl) and (code()_sqrtl) will end up on the same cache
line.)
)
(:d par
These two optimizations more than doubled the speed of my program.)
)
(div((id "Parallel"))
(h2()Parallel Execution)
(:d par Multiple evaluations of the Z function is an embarrassingly
parallel task. My program, written with MPI, has a parent process which
deals out small chunks of 1,000,000 samples of the
domain to each of the child processes. As each child process
finishes, it is assigned a new chunk. The total CPU time consumed
in calculating the first 2,354,000,000 samples was approximately 100
hours, which was spread out over 6 to 9 processors.
(:p The entire computation was not done in a single run.)
)
(:d par A mildly clever innovation I implemented was to prevent a
multi-gigabyte temporary file from forming as the calculation
progressed. As each chunk is completed, it is immediately
losslessly compressed into the
(a ((href "http://flac.sourceforge.net"))FLAC) audio format, which
saves storage by a factor of 6. The compression was performed by
a Perl script which periodically polled the file system
(:p every 10 minutes) looking for completed chunks.)
(:d par The CPU time consumed for tabulating 5,000,000 values in
the neighborhood of the trillionth zero was about 35 hours, spread
out over 6 processors. The chunk size was 10,000.
Clearly high values take longer to calculate.)
)
(div((id "Tabulations"))
(h2()Tabulations)
(:d par This project produced two tabulations of the
Riemann-Siegel Z function
(:p from which the values of the Riemann Zeta function along the
critical line (:c Re(:p s)) = 0.5 can easily calculated) which may be
useful to curious mathematicians, or entertaining to curious
listeners. These tabulations are stored as audio files. For
completeness, the issues of horizontal and vertical scaling
mentioned in the previous section will
be repeated here. The tabulations produced are:)
(ol()
(li()
A tabulation from sample number N=0 to
sample number N=2,353,999,999. The domain and range of the Z
values stored in the file are:
(ul()
(li() Let N be the sample number, where N= 0,1,2,...)
(li() Let (em()y) = 3+(:p N/4) )
(li () Let (em()t) = (em()y)/log (em()y) )
(li() Then, the N-th sample has the value of
256*(:c Z(:p (em()t))), rounded to the nearest 16-bit integer))
(:d par Therefore, the range of evaluated values is approximately
2.7 "<" (em()t) "<" 29,143,636.6.)
)
(li()
A (a ((href "other-files/sample_31801772196387.flac"))
downloadable tabulation
)from sample number N=31801772196387 to
N=31801772196387+4999999. The domain and range of the Z
values stored in the file are:
(ul()
(li() Let N be the sample number, where N=
31801772196387, 31801772196388, 31801772196389, ...)
(li() Let (em()y) = 3+(:p N/4) )
(li() Let (em()t) = (em()y)/log (em()y) )
(li() Then, the N-th sample has the value of
128*(:c Z(:p (em()t))), rounded to the nearest 16-bit integer))
(:d par Therefore, the range of evaluated values is approximately
2.67653395646999e+11 "<" (em()t) "<"2.67653436311838e+11 . )
))
(:d par The tabulations are losslessly compressed into the
(a ((href "http://flac.sourceforge.net"))FLAC) audio format. The
FLAC command line to uncompress the files to stdout is
(code() "flac -d -c -fr -fl filename.flac"),
(:p You may wish to change the (code() "-fl") flag to
(code()"-fb") if you have a big-endian machine. Intel machines
are little-endian.) Note that the first tabulation would produce
a 4.4 gigabyte uncompressed file, which might break your operating
system's limit of 2 GB files. In such a case, you should pipe the
output directly to whatever analysis program you are running.)
(h3()Obtaining the long tabulation)
(:d par If you are interested in obtaining the first long tabulation of
2.354 billion values, please contact me at my permanent e-mail address
(code ()kenta at cs dot stanford dot edu). The tabulation is
approximately 800MB compressed. The easiest way to transfer the
file would be for you to set up an anonymous FTP account where I
can dump the file.)
(h3()Audio files)
(:d par In the limited space of my account are contained
(a ((href "other-files/riemann_low.mp3")) riemann_low.mp3) which
is the initial part of the first tabulation, and
(a ((href "other-files/riemann_hi.mp3")) riemann_hi.mp3) which is
the entirety of the second tabulation. It is interesting to note
that these files sound radically different. Over the course of a
trillion zeros, the sound of the Z function slowly changes. The
sound of the Z function does capture the sound of the Zeta
function: My
(a ((href "other-files/riemann_low.mp3"))
first file of the Z function) sounds quite similar to
(a ((href
"http://www.maths.ex.ac.uk/~mwatkins/zeta/munafo-zetasound.htm"))
Robert Munafo's calculation of the Zeta function.)
)
(:d par I also have available as a 350MB MP3 file
the 2,354,000,000 sample tabulation. Here I reveal the that the
number 2,354,000,000 was not chosen at random: The audio file,
played at 44.1 kHz, plays for 14.8 hours, or 0.618 day. The ratio
0.618 is approximately the golden ratio, which is aesthetically
the (:quote most irrational) number. That
is, if the sound file is played continuously repeatedly for many
days, but the listener is only awake or present for certain hours
of the day, the ratio 0.618 maximizes the (:quote coverage) heard
on different days.
)
(:d par I have not yet listened to the file in its entirety.
Unfortunately, the ambient noise
produced by my computer (:p by the fan and the hard drives)
masks the sound a lot, so it is recommended that the audio file
(:p 350 MBmb) be written or downloaded to a portable MP3 device
(:p such as a MP3 CD player) and
played. Please contact me at the address given above to obtain
the file.)
)
[(div((id "Practical"))
(h2 () Practical)
(ul()
(li() how to scale vertically?)
(li() how to scale horizontally?)
(li() parallelism)
(li() architecture of how streams come together)
(li() checking the log gamma)
(li() checking around the trillionth zero, "31801772196387")
(li() comparison of mp3 encoding)
(li() critical optimizations, and timings)
)
(div()(h3()Things produced)
(ul()
(li() samples 0 .. 2,354,000,000)
(li() rms=??)
(li() in the area of the (:c 10(sup()12)th) zero)
)
))]
(div((id "Analyses"))
(h2() Analyses)
(:d par The following are some simple analyses of the Z
tabulations. They only represent examples of things which could
be done
by curious mathematicians seeking to investigate the Riemann Zeta
function empirically.
For brevity, the tabulation from 0 to 2,354,000,000
will be referred to as the LOW dataset, and the tabulation in the
region of the trillionth zero will be referred to as the HIGH
dataset.)
(table((border "1") )
(tr()(th()) (th()LOW)(th()HIGH))
(tr() (td()Samples)(td() 2354000000)(td()5000000) )
(tr()(td()Min) (td() -68.13)(td()-134.27))
(tr()(td()Max) (td() 74.27)(td()109.58))
(tr()(td()Mean) (td() -1.1331e-05) (td()0.0017698) )
(tr()(td()Root Mean Square value)(td()3.9442)(td()5.0539))
(tr()(td()Approx. num zeros)(td() 66548266)(td()158297))
(tr()(td()Largest spacing between zeros)
(td() 159 samples)(td()114 samples))
(tr()(td()Mean spacing between zeros)
(td()35.373
samples)
(td()31.586 samples))
)
(:d par The number of zeros is approximate because of the 16 bit
precision of the data. The function may creep just barely above
(:p or below) zero and not cause a sign change.)
(h3() Distribution of the spacings between zeros)
(:d par The following plot gives the probability density of the
spacings of the zeros for the two data sets. The mean spacing has
been normalized to unity for each data set (:p separately). Note
that the sample data was already scaled horizontally
by a factor of log(em()t)
(:p that is, the density of zeros). The
LOW data set is given in red as
(:quote(code() line 1 *)), and the HIGH data
set is given in blue as (:quote(code() line 2 +)) .)
(table((border "1"))
(tr()(td()
(img ((src "other-files/gue2.png")
(alt "distribution of zeros")
(width "648")(height "477")
):x))))
(:d par The LOW data set has a long tail which extends off the
right hand edge of the graph.
The image above is quite similar to a computation performed by
Odlyzko in the region of the (:c 10(sup()16)th) zero
(:p excerpted from Odlyzko 2001).)
(table((border "1"))
(tr()(td()
(img ((src "other-files/odlyzko_gue.png")
(alt "distribution of zeros by Odlyzko")
(width "747")(height "703")
):x))))
(h3()Fourier Transform)
(:d par The following two plots are the power spectrum
(:p on a
logarithmic
scale)
of the LOW
data set of (:c 2(sup()24)) samples between sample number 168442784 and
sample number 185219999. The peak frequency can be clearly heard
in the audio file.)
(table((border "1"))
(tr()(td()(img ((src "other-files/four1.gif")
(alt "power spectrum")(width "735")(height "612")
):x))))
(:d par The second plot is a zoom-up of the peak
frequency, showing an interesting crown-like structure.)
(table((border "1"))
(tr()(td()
(img ((src "other-files/four2.gif")
(alt "power spectrum zoom")(width "740")(height "612")
):x))))
(:d par This last plot is the power spectrum of the HIGH data set.)
(table((border "1"))
(tr()(td()
(img ((src "other-files/fourhi.gif")(width "590")(height "496")
(alt "power spectrum")
):x))))
(:d par The Fourier transforms were performed with
(a ((href "http://www.fftw.org/"))FFTW). The size 16,777,216
transform took 17.5 seconds to compute.)
)
(div((id "Source"))
(h2()Source Code)
(:d par The C++ code to tabulate values of the Riemann-Siegel Z function is
contained in
(a ((href "other-files/riemann.src.tar.gz"))riemann.src.tar.gz).
A hodgepodge of other files used in this project are available in
(a ((href "other-files/"))this directory).
)
)
[largest space 159
total zeros 66548266
numsamples 2354000000 min -19252 max 19359
summ -6828225 sumsq 2399938674803437]
(div ((id "References"))
(h2()References)
(p() Edwards, H. M. Riemann's Zeta Function. New York: Academic
Press, 1974.)
(p() Fast algorithms for multiple evaluations of the Riemann zeta
function, A. M. Odlyzko and A. Schoenhage, Trans. Am. Math. Soc.,
309 (:p 1988), pp. 797-809.)
(p() The 10^22-nd zero of the Riemann zeta function, A. M. Odlyzko.
Dynamical, Spectral, and Arithmetic Zeta Functions, M. van
Frankenhuysen and M. L. Lapidus, eds., Amer. Math. Soc.,
Contemporary Math. series, no. 290, 2001, pp. 139-144.)
(p() (:quote The Riemann-Siegel formula and large scale computations of
the Riemann Zeta function), Glendon Pugh, Master's thesis,
University of New Brunswick, 1998)
)
(hr():x)
(p() (a ((href "http://validator.w3.org/check/referer"
["http://validator.w3.org/check?uri=http://www.mit.edu/~kenta/six/parallel/2-Final-Report.html"]))
(img ((src "other-files/valid-xhtml10.png")
(alt "Valid XHTML 1.0!")
(height "31")
(width "88")):x))
(a ((href "http://jigsaw.w3.org/css-validator/check/referer"))
(img ((src "other-files/vcss.gif")(width "88")(height "31")
(alt "Valid CSS!")) :x))
This page is best viewed with a standards-compliant browser such as
(a ((href "http://www.mozilla.org/"))Mozilla)
(:p with which you may wish to turn on the
(strong()Site Navigation Bar) in
the View Menu). The internal anchors which are used in the
(a ((href "#TOC")) Table of Contents)
are known not to work for
older versions of Netscape, which do not support the ID attribute
on the DIV tag.
)
))