Implementing the Riemann-Siegel Formula to find zeros of the Riemann Zeta Function with Roots.jl

[matthewshawnkehoe]

I am interested in implementing a version of the Riemann-Siegel formula for calculating non-trivial zeros of the Riemann zeta function in Julia. I can "successfully" get "manually finding roots by stepping through an interval and looking for sign changes" to work. I can also "successfully" get find_zeros from the Roots.jl package to work (although, this is a "bad" approach for solving this problem and appears to auto-default to the Bisection method). This method might work better for this problem if one had access to other root-finding methods such as Brent's method.

However, I cannot seem to get roots.jl (from IntervalRootFinding) to work.The issue appears to be the way I pass through rtsArray = roots(f, a, b) on line 40 of the code below. Through looking at the example here, one can write

f(x) = sin(x) - 0.1*x^2 + 1

and then

roots(f, -10..10)

which returns

4-element Array{Root{Interval{Float64}},1}:
 Root([3.14959, 3.1496], :unique)
 Root([-4.42654, -4.42653], :unique)
 Root([-1.08205, -1.08204], :unique)
 Root([-3.10682, -3.10681], :unique)

I don't understand why I cannot do the same thing with RiemennZ (as is done in findRootsSlow and findRoots in my code below). The second version of my code is pasted below.

#=**
**
**    Riemann-Siegel Formula for roots of Zeta(s) on the critical line.
**
**************************************************************************
**    Axion004
**    10/29/2022
**
**    This program finds roots of Zeta(s) using the Riemann-Siegel
**    formula. The Riemann–Siegel theta function is approximated by
**    Stirling's approximation. It also uses an interpolation method to
**    locate zeros. The coefficients for R(t) are handled by the Taylor
**    Series approximation originally listed by Haselgrove in the 1960s.
**    It is necessary to approximate these coefficients in order to
**    increase computational speed.

Future Work

1. Implement a cleaner version of findRoots(). Investigate interval methods.
2. Use GPU and parallel processing capability. A lot of good ideas are discussed in
https://www.youtube.com/watch?v=8E0qdO_jRZ0.
3. Calculate all of the Riemann-Siegel coefficients explicitly.
4. Implement the Odlyzko–Schönhage algorithm.
**************************************************************************=#

using Roots
using IntervalArithmetic, IntervalRootFinding

"""
* Search for zeros of the Riemann-Siegel Z function in the interval [a,b] by
    the roots function in  IntervalRootFinding.jl. See
    https://github.com/JuliaIntervals/IntervalRootFinding.jl and
    https://juliaintervals.github.io/pages/packages/.
* @param a - the starting value of the interval.
* @param b - the ending value of the interval.
* @return an array of zeros of zeta(s) on the critical line.
"""
function findRootsIntervalArithmetic(a::Float64, b::Float64)
    f(x) = RiemennZ(x, 4);
    rtsArray = roots(f, a, b) # I don't know why this doesn't work.
    println(size(rtsArray))
    #for root in rtsArray
     #   println(root)
    #end
end

"""
* Search for zeros of the Riemann-Siegel Z function in the interval [a,b] by
    the find_zeros function in Roots.jl. See
    https://juliamath.github.io/Roots.jl/dev/roots/#Searching-for-all-zeros-in-an-interval.
* @param a - the starting value of the interval.
* @param b - the ending value of the interval.
* @return an array of zeros of zeta(s) on the critical line.
"""
function findRoots(a::Float64, b::Float64)
    f(x) = RiemennZ(x, 4);
    rtsArray = find_zeros(f, a, b) # This isn't accurate when b > 1000.
    # rtsArray = find_zeros(f, a, b, Roots.Brent()) # I don't know why this doesn't work.
    println(size(rtsArray))
    #for root in rtsArray
        #println(root)
    #end
end

"""
    * Extremely inefficient way of calculating roots manually.
"""
function findRootsSlow(a::Float64, b::Float64)
    f(x) = RiemennZ(x, 4);
    step_size = 0.01
    while a < b
        a += step_size
        # This is a very inefficient way of looking for roots
        if sign(f(a)) != sign(f(a + step_size))
            x = fzero(f, [a, a + step_size])
            # Print the roots to the console
            println(x)
            #io = open("roots.txt","a")
            #    println(io,x)
            #close(io)
        end
    end
end

"""
    * Riemann-Siegel theta function using the approximation by the
      Stirling series.
    * @param t - the value of t inside the Z(t) function.
    * @return Stirling's approximation for theta(t).
"""
function theta(t::Float64)
    return (t/2.0 * log(t/(2.0*pi)) - t/2.0 - pi/8.0
            + 1.0/(48.0*t) + 7.0/(5760*t^3))
end

"""
    * Computes the floor of the absolute value term passed in as t.
    * @param t - the value of t inside the Z(t) function.
    * @return floor of the absolute value of t.
"""
function fAbs(t::Float64)
    return floor(abs(t))
end


"""
    * Riemann-Siegel Z(t) function implemented per the Riemenn Siegel
         formula. See http://mathworld.wolfram.com/Riemann-SiegelFormula.html
         for details.
    * @param t - the value of t inside the Z(t) function.
         @param r - referenced for calculating the remainder terms by the
         Taylor series approximations.
    * @return the approximate value of Z(t) through the Riemann-Siegel
         formula
"""
function RiemennZ(t::Float64, r::Int)
    twopi = pi * 2.0;
    val = sqrt(t/twopi)
    n = fAbs(val)
    sum = 0.0

    for i in 1:n
          sum += (cos(theta(t) - t * log(i))) / sqrt(i)
    end
    sum = 2.0 * sum

    frac = val - n
    k = 0
    R = 0.0

    # Calculate every remainder term through the Taylor Series coefficients.
    # The coefficients are defined in the Riemann-Siegel C function.
    while k <= r
        R = R + C(k, 2.0*frac-1.0) * (t / twopi)^(k * (-0.5))
        k += 1
    end

    remainder = (-1)^(n-1) * (t / twopi)^(-0.25) * R
    return sum + remainder
end


"""
    * C terms for the Riemann-Siegel formula. See
         Haselgrove (1960): Tables of the Riemann Zeta Function,
         doi:10.2307/3614148  for details.
         Calculates the Taylor Series coefficients for C0, C1, C2, C3,
         and C4.
    * @param n - the number of coefficient terms to use.
    * @param z - referenced per the Taylor series calculations.
    * @return the Taylor series approximation of the remainder terms.
"""
function C(n::Int, z::Float64)
    if n==0
        return (.38268343236508977173 * z^0.0
            +.43724046807752044936 * z^2.0
            +.13237657548034352332 * z^4.0
            -.01360502604767418865 * z^6.0
            -.01356762197010358089 * z^8.0
            -.00162372532314446528 * z^10.0
            +.00029705353733379691 * z^12.0
            +.00007943300879521470 * z^14.0
            +.00000046556124614505 * z^16.0
            -.00000143272516309551 * z^18.0
            -.00000010354847112313 * z^20.0
            +.00000001235792708386 * z^22.0
            +.00000000178810838580 * z^24.0
            -.00000000003391414390 * z^26.0
            -.00000000001632663390 * z^28.0
            -.00000000000037851093 * z^30.0
            +.00000000000009327423 * z^32.0
            +.00000000000000522184 * z^34.0
            -.00000000000000033507 * z^36.0
            -.00000000000000003412 * z^38.0
            +.00000000000000000058 * z^40.0
            +.00000000000000000015 * z^42.0)
    elseif n==1
        return (-.02682510262837534703 * z^1.0
            +.01378477342635185305 * z^3.0
            +.03849125048223508223 * z^5.0
            +.00987106629906207647 * z^7.0
            -.00331075976085840433 * z^9.0
            -.00146478085779541508 * z^11.0
            -.00001320794062487696 * z^13.0
            +.00005922748701847141 * z^15.0
            +.00000598024258537345 * z^17.0
            -.00000096413224561698 * z^19.0
            -.00000018334733722714 * z^21.0
            +.00000000446708756272 * z^23.0
            +.00000000270963508218 * z^25.0
            +.00000000007785288654 * z^27.0
            -.00000000002343762601 * z^29.0
            -.00000000000158301728 * z^31.0
            +.00000000000012119942 * z^33.0
            +.00000000000001458378 * z^35.0
            -.00000000000000028786 * z^37.0
            -.00000000000000008663 * z^39.0
            -.00000000000000000084 * z^41.0
            +.00000000000000000036 * z^43.0
            +.00000000000000000001 * z^45.0)
    elseif n==2
        return(+.00518854283029316849 * z^0.0
            +.00030946583880634746 * z^2.0
            -.01133594107822937338 * z^4.0
            +.00223304574195814477 * z^6.0
            +.00519663740886233021 * z^8.0
            +.00034399144076208337 * z^10.0
            -.00059106484274705828 * z^12.0
            -.00010229972547935857 * z^14.0
            +.00002088839221699276 * z^16.0
            +.00000592766549309654 * z^18.0
            -.00000016423838362436 * z^20.0
            -.00000015161199700941 * z^22.0
            -.00000000590780369821 * z^24.0
            +.00000000209115148595 * z^26.0
            +.00000000017815649583 * z^28.0
            -.00000000001616407246 * z^30.0
            -.00000000000238069625 * z^32.0
            +.00000000000005398265 * z^34.0
            +.00000000000001975014 * z^36.0
            +.00000000000000023333 * z^38.0
            -.00000000000000011188 * z^40.0
            -.00000000000000000416 * z^42.0
            +.00000000000000000044 * z^44.0
            +.00000000000000000003 * z^46.0)
    elseif n==3
        return (-.00133971609071945690 * z^1.0
            +.00374421513637939370 * z^3.0
            -.00133031789193214681 * z^5.0
            -.00226546607654717871 * z^7.0
            +.00095484999985067304 * z^9.0
            +.00060100384589636039 * z^11.0
            -.00010128858286776622 * z^13.0
            -.00006865733449299826 * z^15.0
            +.00000059853667915386 * z^17.0
            +.00000333165985123995 * z^19.0
            +.00000021919289102435 * z^21.0
            -.00000007890884245681 * z^23.0
            -.00000000941468508130 * z^25.0
            +.00000000095701162109 * z^27.0
            +.00000000018763137453 * z^29.0
            -.00000000000443783768 * z^31.0
            -.00000000000224267385 * z^33.0
            -.00000000000003627687 * z^35.0
            +.00000000000001763981 * z^37.0
            +.00000000000000079608 * z^39.0
            -.00000000000000009420 * z^41.0
            -.00000000000000000713 * z^43.0
            +.00000000000000000033 * z^45.0
            +.00000000000000000004 * z^47.0)
    else
        return (+.00046483389361763382 * z^0.0
            -.00100566073653404708 * z^2.0
            +.00024044856573725793 * z^4.0
            +.00102830861497023219 * z^6.0
            -.00076578610717556442 * z^8.0
            -.00020365286803084818 * z^10.0
            +.00023212290491068728 * z^12.0
            +.00003260214424386520 * z^14.0
            -.00002557906251794953 * z^16.0
            -.00000410746443891574 * z^18.0
            +.00000117811136403713 * z^20.0
            +.00000024456561422485 * z^22.0
            -.00000002391582476734 * z^24.0
            -.00000000750521420704 * z^26.0
            +.00000000013312279416 * z^28.0
            +.00000000013440626754 * z^30.0
            +.00000000000351377004 * z^32.0
            -.00000000000151915445 * z^34.0
            -.00000000000008915418 * z^36.0
            +.00000000000001119589 * z^38.0
            +.00000000000000105160 * z^40.0
            -.00000000000000005179 * z^42.0
            -.00000000000000000807 * z^44.0
            +.00000000000000000011 * z^46.0
            +.00000000000000000004 * z^48.0)
    end
end


# Main
function main()
    println("Zeroes on the critical line Zeta(1/2 + it).")
    a = 1.0
    b = 100.0
    #findRootsSlow(a, b)
    findRoots(a, b)
    #findRootsIntervalArithmetic(a, b)
end

main()

Uncommenting findRootsIntervalArithmetic(a, b)on line 288 returns

ERROR: MethodError: no method matching roots(::var"#f#72", ::Float64, ::Float64)
Closest candidates are:
  roots(::Function, ::Function, ::Any) at ~/.julia/packages/IntervalRootFinding/OdQlQ/src/roots.jl:126
  roots(::Function, ::Function, ::Any, ::Type{C}) where C<:Contractor at ~/.julia/packages/IntervalRootFinding/OdQlQ/src/roots.jl:126
  roots(::Function, ::Function, ::Any, ::Type{C}, ::Float64) where C<:Contractor at ~/.julia/packages/IntervalRootFinding/OdQlQ/src/roots.jl:139
  ...
Stacktrace:
 [1] findRootsIntervalArithmetic(a::Float64, b::Float64)
   @ Main ~/RiemannSiegel/RiemannSiegel.jl:40
 [2] main()
   @ Main ~/RiemannSiegel/RiemannSiegel.jl:288
 [3] top-level scope
   @ ~/RiemannSiegel/RiemannSiegel.jl:291

So, I am likely doing something wrong on line 40 when I write rtsArray = roots(f, a, b). I think that I have incorrectly called roots and this is expected to error out.

[matthewshawnkehoe]

Another "bad" approach would be to subdivide the native find_zeros function with smaller intervals.

function findRoots(a::Float64, b::Float64, intervalSize::Float64)
    # This function is inaccurate when the interval size is large. 
    # So, subdivide the interval into smaller intervals.
    startPoint = a
    endPoint = (b % intervalSize) + intervalSize
    numIntervals = b ÷ intervalSize
    rtsArray = zeros(0)
    @inbounds for i=1:numIntervals
        f(x) = RiemennZ(x, 4);
        tempArray = find_zeros(f, startPoint, endPoint)
        startPoint = endPoint
        endPoint = endPoint + intervalSize
        append!(rtsArray, tempArray)
    end
    for root in rtsArray
        println(root)
    end
    println(size(rtsArray))
end 

The test sample of

a= 1.0
b = 74920.9 
intervalSize = 1.0
findRoots(a, b, intervalSize)

correctly returns the first 100,000 zeros on the critical line. Although, this probably isn't the way you want to approach this computation.

[lbenet]

Thanks for opening this issue.

First, let me note that roots in IntervalRootFinding expects an interval as the second argument (and not two floats), as you may notice in the example you quoted and in the error reported. This means that you have to change the line

rtsArray = roots(f, a, b)

in findRootsIntervalArithmetic, to

rtsArray = roots(f, a .. b)

Once you do that, this will still trigger an error, but that is related to the evaluation of RiemennZ using intervals, which currently does not work:

julia> RiemennZ(1..2, 4)
ERROR: MethodError: no method matching RiemennZ(::Interval{Float64}, ::Int64)

This is related to the fact that you have only defined a method for RiemennZ whose first entry is a Float64 and not an interval, which is needed by IntervalRootFinding. So you have to be sure that such an evaluation works, which may not be trivial.

I hope this helps to continue on this.

[Kolaru]

I was about to answer exactly that, but Luis beat me to it :)

I have tried your example a bit further, and you will face another issue a bit further down the line. Name IntervalRootFinding will find no roots, because

julia> RiemennZ(1..100, 4)
∅

This is because there seems to be a bug in IntervalArithmetic i.e. currently

julia> (-1)^(-1..1)
∅

where it should return (0..1) I believe.

Even so, you will face other problem because now everything is replaced by intervals. In particular, n in RiemannZ is no longer an integer but an interval. So what the loop for i in 1:n does is rather... unpredictable[1].

I hope this helps you go forward. The bottomline is simply that if you want tu use intervals, you have to remember they are propagated to result in guaranteed computations.

[1] Because it should probably error to be fair.