Progress Report

March 31, 2002

Riemann audio

I have made much progress on this project. My principal reference for the implementation of the zeta has been Glen Pugh's master thesis.

Like the reference, I am implementing the Riemann-Siegel formula, and putting quite a bit of care that my implementation is right. Thus far, I have uncovered two bugs in Pugh's thesis (I have corresponded with him about them.)

One of the bugs I discovered when I was trying to double-check the coefficents of the remainder term of the Riemann-Siegel formula by re-deriving them according to the method he describes. I got different results than he did: it turned out his description of the derivation is slightly incorrect. (His coefficients however, agree with other published sources, and are probably correct. Unfortunately, I am still unable to independently derive them.)

The other bug was a small typo in an equation.

I am putting more care into my implementation because I am not root-counting (verifying that all the zeros lie on the critical line), like he was in his thesis. An implementation of zeta might be wrong (slightly off by a few decimal digits) and still work for root-counting, so long as roots aren't made to appear or disappear. I want to have the correct value of zeta (up to the sensitivity of human hearing.)

I have chosen not to implement the Odlyzko-Schonhage (Schoenhage) FFT-based method for the zeta. It is so complicated that I cannot understand it.

I have not followed up on a suggestion of Professor Edelman to ask Odlyzko for his code. His code (like Pugh's) is intended for root-counting, and I imagine I would have a hard time modifying it for my purposes because of its complexity.

Currently, my progress is as follows:

• Implemented the complex log gamma function, based on the (real-valued) implementation in Numerical Recipes. Still have not verified (to my satisfaction) that the implentation is correct. This implementation turned out to be easier than expected because of the <complex> class built-in to the C++ Standard Template Library.
• Derived the Chebyshev coefficients for the remainder term of the Riemann-Siegel formula. (More specifically, the Chebyshev coefficients for each of the "C_i(p)" functions in the remainder term.) I am using Chebyshev expansions as opposed to power series (in contrast to Pugh's implementation) because a Chebyshev expansion spreads the error equally on the domain of each of the C_i(p) functions. (The domain is [0..1].) I am slightly worried that the 0/0 singularities at 0.25 and 0.75 of the Psi function has caused the Chebyshev expansions to be messed up.

I expect to have completed within a week an implementation of zeta with which I am satisfied.

Mandelbrot

I investigated with the freeware program FRACTINT the region I proposed to magnify. It turns out that 17-18 decimal places (of the long double type) is insufficient to see anything interesting.

Therefore, in order to do this "right", I will need to implement the algorithms with arbitrary precision arithmetic. This project turned out to be harder than I expected. It is currently low priority and will wait until the Riemann audio project is near completion.