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 floating-point arithmetic and directed rounding.
Consider the 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}$. We will evalute $S_N$ numerically and find an analytical bound for $T_N$.
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}$.
$S_N$ may be calculated easily by summing either forwards or backwards:
function forward_sum(N, T=Float64)
total = zero(T)
for i in 1:N
total += one(T) / (i^2)
end
total
end
function reverse_sum(N, T=Float64)
total = zero(T)
for i in N:-1:1
total += one(T) / (i^2)
end
total
end
reverse_sum (generic function with 2 methods)
To find rigorous bounds for $S_N$, we use "directed rounding", that is, we round downwards for the lower bound and upwards for the upper bound:
N = 10^6
lowerbound_S_N =
with_rounding(Float64, RoundDown) do
forward_sum(N)
end
upperbound_S_N =
with_rounding(Float64, RoundUp) do
forward_sum(N)
end
1.644933066959796
We incorporate the respective bound on $T_N$ to obtain the bounds on $S$, and hence on $\pi$:
N = 10^6
lower_pi =
with_rounding(Float64, RoundDown) do
lower_bound = forward_sum(N) + 1/(N+1)
sqrt(6 * lower_bound)
end
upper_pi =
with_rounding(Float64, RoundUp) do
upper_bound = forward_sum(N) + 1/N
sqrt(6 * upper_bound)
end
lower_pi, upper_pi
(3.1415926534833463,3.1415926536963346)
We may check that the true value of $\pi$ is indeed contained in the interval:
lower_pi < big(pi) < upper_pi
true
We may repeat the calculation summing in the opposite direction for a more accurate answer:
N = 10^6
lower_pi =
with_rounding(Float64, RoundDown) do
lower_bound = reverse_sum(N) + 1/(N+1)
sqrt(6 * lower_bound)
end
upper_pi =
with_rounding(Float64, RoundUp) do
upper_bound = reverse_sum(N) + 1/N
sqrt(6 * upper_bound)
end
lower_pi, upper_pi, lower_pi < big(pi) < upper_pi
(3.1415926535893144,3.141592653590272,true)
BigFloat
¶It is not clear (to me) if the standard sqrt
operation on Float64
s respects the rounding mode. We may use BigFloat
s to guarantee a correctly-rounded sqrt
function, although the calculation is much slower:
set_bigfloat_precision(53) # same precision as `Float64`
N = 10^6
lower_pi =
with_rounding(BigFloat, RoundDown) do
lower_bound = reverse_sum(N, BigFloat) + one(BigFloat)/(N+1)
sqrt(6 * lower_bound)
end
upper_pi =
with_rounding(BigFloat, RoundUp) do
upper_bound = reverse_sum(N, BigFloat) + one(BigFloat)/N
sqrt(6 * upper_bound)
end
lower_pi, upper_pi, lower_pi < big(pi) < upper_pi
(3.1415926535893144e+00 with 53 bits of precision,3.1415926535902718e+00 with 53 bits of precision,true)
(Note that the result seems to be the same as the result with Float64
, so perhaps the sqrt
function does respect the rounding mode.)
In principal, we could attain arbitrarily good precision with higher-precision BigFloat
s, but the result is hampered by the slow convergence of the series.