Calculating rigorous bounds on $\pi$ with interval arithmetic

Department of Physics, Faculty of Sciences, National University of Mexico (UNAM)

Following the book Validated Numerics (Princeton, 2011) by Warwick Tucker, we find rigorous (i.e., guaranteed, or validated) bounds on $\pi$ using interval arithmetic, via the ValidatedNumerics.jl package.



In [7]:

    
using ValidatedNumerics

Calculating $\pi$ via a simple sum

There are many ways to calculate $\pi$. For illustrative purposes, we will use the following sum

$$ S := \sum_{n=1}^\infty \frac{1}{n^2}.$$

It is known that the exact value is $S = \frac{\pi^2}{6}$. Thus, if we can calculate rigorous bounds on $S$, then we can find rigorous bounds on $\pi$.

The idea is to split $S$ up into two parts, $S = S_N + T_N$, with $ S_N := \sum_{n=1}^N \frac{1}{n^2}$ and $T_N := S - S_N = \sum_{n=N+1}^\infty n^{-2}$.

By bounding $T_N$ using integrals from below and above, we can see that $\frac{1}{N+1} \le T_N \le \frac{1}{N}$. Rigorous bounds on $S_N$ are found by calculating it using interval arithmetic.

$S_N$ may be calculated by summing either forwards or backwards. A naive (non-interval) version could be the following:



In [13]:

    
function forward_sum_naive(N)
    S_N = 0.0

    for i in 1:N
        S_N += 1./(i^2)
    end

    S_N
end









    Out[13]:





forward_sum_naive (generic function with 1 method)



In [14]:

    
S = forward_sum_naive(10000)
err = abs(S - pi^2/6.)  # error
S, err









    Out[14]:





(1.6448340718480652,9.999500016122376e-5)

Interval arithmetic

To find rigorous bounds for $S_N$, we use interval arithmetic: each term is enclosed in an interval that is guaranteed to contain the true real value. A first idea is simply to wrap each term in the @interval macro, which converts its arguments into containing intervals:



In [15]:

    
function forward_S_N(N)
    S_N = @interval(0.0)

    for i in 1:N
        S_N += @interval( 1./(i^2) )
    end

    S_N
end









    Out[15]:





forward_S_N (generic function with 1 method)



In [16]:

    
N = 10^5
@time rigorous_approx_S_N = forward_S_N(N)









    



elapsed time: 0.787551828 seconds (236886656 bytes allocated, 24.83% gc time)






    Out[16]:





[1.6449240668871583, 1.644924066909359]

We incorporate the respective bound on $T_N$ to obtain the bounds on $S$, and hence on $\pi$. We can also optimize the code by creating the interval @interval(1) only once:



In [28]:

    
function forward_sum(N)
    S_N = @interval(0.0)
    interval_one = @interval(1.)

    for i in 1:N
        S_N += interval_one / (i^2)
    end
    
    T_N = interval_one / @interval(N, N+1)
    S = S_N + T_N

    sqrt(6*S)
end

N = 10^6
@time S = forward_sum(N)
S, diam(S)









    



elapsed time: 4.652939682 seconds (2080127080 bytes allocated, 27.31% gc time)






    Out[28]:





([3.1415926534833463, 3.1415926536963346],2.1298829366855898e-10)

We can ask for the midpoint--radius representation, which shows that the calculated bounds are correct to around 10 decimal places:



In [24]:

    
midpoint_radius(S)









    Out[24]:





(3.141592653589824,5.834666083615048e-11)

We may check that the true value of $\pi$ is indeed contained in the interval:



In [25]:

    
big(pi) ∈ S









    Out[25]:





true

We may repeat the calculation, but now summing in the opposite direction. Due to the way that floating-point arithmetic works, this gives a more accurate answer.



In [32]:

    
function reverse_sum(N)
    S_N = @interval(0.0)
    interval_one = @interval(1.)

    for i in N:-1:1
        S_N += interval_one / (i^2)
    end
    
    T_N = interval_one / @interval(N, N+1)
    S = S_N + T_N

    sqrt(6*S)
end

N = 10^6
@time S = reverse_sum(N)
S, diam(S)









    



elapsed time: 5.638430487 seconds (2080153808 bytes allocated, 28.92% gc time)






    Out[32]:





([3.1415926535893144, 3.141592653590272],9.57456336436735e-13)

Note that the sqrt function is guaranteed (by the IEEE 754 standard) to give correctly-rounded results, so the resulting bounds on $\pi$ are guaranteed to be correct.

Note also that due to the way that the ValidatedNumerics package works, we can make the code simpler, at the expense of some performance. Only two explicit calls to the @interval macro are now required:



In [31]:

    
function reverse_sum2(N)
    S_N = @interval(0.0)

    for i in N:-1:1
        S_N += 1 / (i^2)
    end
    
    T_N = 1 / @interval(N, N+1)
    S = S_N + T_N

    sqrt(6*S)
end

N = 10^6
@time S = reverse_sum2(N)
S, diam(S)









    



elapsed time: 6.699738655 seconds (2264155172 bytes allocated, 27.45% gc time)






    Out[31]:





([3.1415926535893144, 3.141592653590272],9.57456336436735e-13)