Rigorous approximation of π\pi

Download the notebook for this tutorial.

Introduction

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 IntervalArithmetic.jl package. First, let's import the package

using IntervalArithmetic

Calculating π\pi via a simple sum

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

S=n=11n2. S = \sum_{n=1}^\infty \frac{1}{n^2}.

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

The idea is to split SS up into two parts, S=SN+TNS = S_N + T_N, with SN=n=1N1n2 S_N = \sum_{n=1}^N \frac{1}{n^2} and TN=SSN=n=N+1n2T_N = S - S_N = \sum_{n=N+1}^\infty n^{-2}.

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

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

function forward_sum_naive(N)
    S_N = 0.0

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

    S_N
end

S = forward_sum_naive(10000)
err = abs(S - pi^2/6.)  # error
@show S, err
(S, err) = (1.6448340718480652, 9.999500016122376e-5)

Interval method

To find rigorous bounds for SNS_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:

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

N = 10^5
@time rigorous_approx_S_N = forward_S_N(N)
  2.710170 seconds (16.48 M allocations: 705.208 MiB, 5.64% gc time)
[1.64492, 1.64493]

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

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)
@show S, diam(S)
  0.657669 seconds (5.45 M allocations: 178.428 MiB, 7.03% gc time)
(S, diam(S)) = ([3.14159, 3.1416], 2.1298829366855898e-10)

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

midpoint_radius(S)
(3.14159265358984, 1.0649436887888442e-10)

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

big(π) ∈ S
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.

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)
@show S, diam(S)
  0.169420 seconds (3.03 M allocations: 77.610 MiB, 11.94% gc time)
(S, diam(S)) = ([3.14159, 3.1416], 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 IntervalArithmetic 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:

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)
@show S, diam(S)
  0.137809 seconds (2.06 M allocations: 48.472 MiB, 15.39% gc time)
(S, diam(S)) = ([3.14159, 3.1416], 9.57456336436735e-13)