Some fun with Π in Julia

Image courtesy of Cormullion, code here.
This post is available as a Jupyter notebook here.

π in Julia
π via inline assembly instructions
π using a Taylor series expansions
1. Madhava's formula
2. Machin's approach
Finding guaranteed bounds on π
1. Using standard floating-point arithmetic
Summing a series using interval arithmetic
Calculating an area

π in Julia

Like most technical languages, Julia provides a variable constant for π. However Julia's handling is a bit special.

pi

π = 3.1415926535897...

It can also be accessed via the unicode symbol (you can get it at the REPL or in a notebook via the TeX completion \pi followed by a tab)

π

π = 3.1415926535897...

You'll notice that it doesn't print like an ordinary floating point number: that's because it isn't one.

typeof(pi)

Irrational{:π}

π and a few other irrational constants are instead stored as special Irrational values, rather than being rounded to Float64. These act like ordinary numeric values, except that they can are converted automatically to any floating point type without any intermediate rounding:

1 + pi # integers are promoted to Float64 by default

4.141592653589793

Float32(1) + pi # Float32

4.141593f0

This is particularly useful for use with arbitrary-precision BigFloats, as π can be evaluated to full precision (rather than be truncated to Float64 and converted back).

BigFloat(1) + pi # 256 bits by default

4.141592653589793238462643383279502884197169399375105820974944592307816406286198

If π were stored as a Float64, we would instead get

BigFloat(1) + Float64(pi)

4.141592653589793115997963468544185161590576171875000000000000000000000000000000

In fact BigFloat (which uses the MPFR library) will compute π on demand to the current precision, which is set via setprecision. This provides an easy way to get its digits:

# to 1024 bits
setprecision(BigFloat, 1024) do
    BigFloat(pi)
end

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724586997

The last few printed digits may be incorrect due to the conversion from the internal binary format of BigFloat to the decimal representation used for printing. This is just a presentation issue, however – the internal binary representation is correctly rounded to the last bit.

Another neat property of Irrationals is that inequalities are correct:

Float64(pi) < pi < nextfloat(Float64(pi))

true

π via inline assembly instructions

by Simon Byrne

Julia provides a very low-level llvmcall interface, which allows the user to directly write LLVM intermediate representation, including the use of inline assembly. The following snippet calls the fldpi instruction ("floating point load pi") which loads the constant π onto the floating point register stack (this works only on x86 and x86_64 architectures)

function asm_pi()
    Base.llvmcall(
    """ %pi = call double asm "fldpi", "={st}"()
        ret double %pi""",
    Float64, Tuple{})
end

asm_pi (generic function with 1 method)

asm_pi()

3.141592653589793

We can look at the actual resulting code that is generated:

@code_native asm_pi()

    	.section	__TEXT,__text,regular,pure_instructions
    Filename: In[10]
    	pushq	%rbp
    	movq	%rsp, %rbp
    Source line: 2
    	fldpi
    	fstpl	-8(%rbp)
    	movsd	-8(%rbp), %xmm0         ## xmm0 = mem[0],zero
    	popq	%rbp
    	retq

If you're wondering what the rest of these instructions are doing:

the pushq and movq adds to the call stack frame.
fldpi pushes π to the x87 floating point register stack

x87 is the older legacy floating point instruction set dating back to the original Intel 8087 coprocessor.

fstpl and movsd moves the value to the SSE floating point register xmm0

Julia, like most modern software, uses the newer SSE instruction set for its floating point operations. This also allows us to take advantage of things like SIMD operations.

popq and retq pops the call stack frame.

π using a Taylor series expansions

by Luis Benet, Instituto de Ciencias Físicas, Universidad Nacional Autónoma de México (UNAM))

This will demonstrate how to evaluate π using various Taylor series expansions via the TaylorSeries.jl package.

using TaylorSeries

Madhava's formula

One of the standard trigonmetric identities is

\tan\left( \frac{\pi}{6} \right) = \frac{1}{\sqrt{3}}.

Therefore, by taking the Taylor expansion of $6 \arctan(x)$ around 0 we may obtain the value of $\pi$ , by evaluating it at $1/\sqrt{3}$ , a value which is within the radius of convergence.

We obtain the Taylor series of order 37th, using BigFloats:

series1 = 6atan( Taylor1(BigFloat, 37) )
convert(Taylor1{Rational{BigInt}},series1)

6//1 t - 2//1 t³ + 6//5 t⁵ - 6//7 t⁷ + 2//3 t⁹ - 6//11 t¹¹ + 6//13 t¹³ - 2//5 t¹⁵ + 6//17 t¹⁷ - 6//19 t¹⁹ + 2//7 t²¹ - 6//23 t²³ + 6//25 t²⁵ - 2//9 t²⁷ + 6//29 t²⁹ - 6//31 t³¹ + 2//11 t³³ - 6//35 t³⁵ + 6//37 t³⁷ + 𝒪(t³⁸)

Note that the series above has only odd powers, so we will be using in this case 18 coefficients.

Evaluating that expression in $1/\sqrt{3}$ we get

pi_approx1 = evaluate(series1, 1/sqrt(big(3)))

3.141592653647826046431202390582141253830948237428790668441592864548346569098516

Then, the 37th order Taylor expansion yields a value which differs from $\pi$ in:

abs(pi - pi_approx1)

5.803280796855900730263836963377883805368484746664827224053016281231814650118929e-11

To obtain more accurate results, we may simply increase the order of the expansion:

series2 = 6atan( Taylor1(BigFloat,99) ) # 49 coefficients of the series
pi_approx2 = evaluate(series2, 1/sqrt(BigInt(3)))

3.141592653589793238462643347272152237127662423839333289949470742535834074912581

abs(pi - pi_approx2)

3.600735064706950697553577253102547384977198233137361734413175534929622111373249e-26

This formulation is one of the Madhava or Gregory–Leibniz series:

\pi = 6 \sum_{n=0}^{\infty} (-1)^n \frac{(1/\sqrt{3})^{2n+1}}{2n+1}.

Machin's approach

Following the same idea, John Machin derived an algorithm which converges much faster, using the identity

\frac{\pi}{4} = 4 \arctan\left(\frac{1}{5}\right) - \arctan\left(\frac{1}{239}\right).

Following what we did above, using again a 37th Taylor expansion:

ser = atan( Taylor1(BigFloat, 37) )
pi_approx3 = 4*( 4*evaluate(ser, 1/big(5)) - evaluate(ser, 1/big(239)) )
3.141592653589793238462643383496777424642594661632063407072684671069773618535135

abs(pi - pi_approx3)
2.17274540445425262256957586097740078761957212248936631045983596428448951876822e-28

Finding guaranteed bounds on π

by David P. Sanders, Department of Physics, Faculty of Sciences, National University of Mexico (UNAM)

Using standard floating-point arithmetic

We will calculate guaranteed (i.e., validated, or mathematically rigorous) bounds on $\pi$ using just floating-point arithmetic. This requires "directed rounding", i.e. the ability to control in which direction floating-point operations are rounded.

This is based on the book Validated Numerics (Princeton, 2011) by Warwick Tucker.

Consider the infinite series

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

whose exact value is known to be $S = \frac{\pi^2}{6}$ . Thus, if finding guaranteed bounds on $S$ will give guaranteed bounds on $\pi$ .

The idea is to split $S$ up into two parts, $S = S_N + T_N$ , where $S_N := \sum_{n=1}^N \frac{1}{n^2}$ contains the first $N$ terms, and $T_N := S - S_N = \sum_{n=N+1}^\infty \frac{1}{n^2}$ contains the rest (an infinite number of terms).

We will evalute $S_N$ numerically, and use the following analytical bound for $T_N$ :

\frac{1}{N+1} \le T_N \le \frac{1}{N}

This is obtained by approximating the sum in $T_N$ using integrals from below and above:

\int_{x=N+1}^\infty \frac{1}{x^2} dx \le T_N \le \int_{x=N}^\infty \frac{1}{x^2} dx.

$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 =
    setrounding(Float64, RoundDown) do
        forward_sum(N)
    end

upperbound_S_N =
    setrounding(Float64, RoundUp) do
        forward_sum(N)
    end

(lowerbound_S_N, upperbound_S_N)

(1.6449330667377557,1.644933066959796)

We incorporate the respective bound on $T_N$ to obtain the bounds on $S$ , and hence on $\pi$ :

N = 10^6

lower_π =
    setrounding(Float64, RoundDown) do
        lower_bound = forward_sum(N) + 1/(N+1)
        sqrt(6 * lower_bound)
    end

upper_π =
    setrounding(Float64, RoundUp) do
        upper_bound = forward_sum(N) + 1/N
        sqrt(6 * upper_bound)
    end

(lower_π, upper_π, lowerbound_S_N)

(3.1415926534833463,3.1415926536963346,1.6449330667377557)

upper_π - lower_π

2.1298829366855898e-10

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

lower_π < pi < upper_π

true

Summing in the opposite direction turns out to give a more accurate answer:

N = 10^6

lower_π =
    setrounding(Float64, RoundDown) do
        lower_bound = reverse_sum(N) + 1/(N+1)
        sqrt(6 * lower_bound)
    end

upper_π =
    setrounding(Float64, RoundUp) do
        upper_bound = reverse_sum(N) + 1/N
        sqrt(6 * upper_bound)
    end

(lower_π, upper_π)
(3.1415926535893144,3.141592653590272)

upper_π - lower_π

9.57456336436735e-13

lower_π < pi < upper_π
true

In principle, we could attain arbitrarily good precision with higher-precision BigFloats, but the result is hampered by the slow convergence of the series.

Summing a series using interval arithmetic

We repeat the calculation using interval arithmetic, provided by the ValidatedNumerics.jl package.

using ValidatedNumerics

setdisplay(:standard)  # abbreviated display of intervals
6

N = 10000
S = forward_sum(N, Interval)
S += 1/(N+1) .. 1/N  # interval bound on the remainder of the series
π_interval = √(6S)

[3.14159, 3.1416]

Here we used an abbreviated display for the interval. Let's see the whole thing:

setdisplay(:full)
π_interval

Interval(3.1415926488148807, 3.141592658365341)

Its diameter (width) is

diam(π_interval)
9.550460422502738e-9

Thus, the result is correct to approximately 8 decimals.

In this calculation, we used the fact that arithmetic operations of intervals with numbers automatically promote the numbers to an interval:

setdisplay(:full)  # full interval display
Interval(0) + 1/3^2

Interval(0.1111111111111111, 0.11111111111111112)

This is an interval containing the true real number $1/9$ (written 1//9 in Julia):

1//9 ∈ convert(Interval{Float64}, 1/3^2)

true

Finally, we can check that the true value of $\pi$ is indeed inside our interval:

pi ∈ π_interval

true

Calculating an area

Although the calculation above is simple, the derivation of the series itself is not. In this section, we will use a more natural way to calculate $\pi$ , namely that the area of a circle of radius $r$ is $A(r) = \pi r^2$ . We will calculate the area of one quadrant of a circle of radius $r=2$ , which is equal to $\pi$ :

using Plots; gr();

f(x) = √(4 - x^2)

f (generic function with 1 method)

plot(f, 0, 2, aspect_ratio=:equal, fill=(0, :orange), alpha=0.2, label="")

The circle of radius $r=2$ is given by $x^2 + y^2 = 2^2 = 4$ , so

\pi = \frac{1}{4} A(2) = \int_{x=0}^2 y(x) \, dx = \int_{x=0}^2 \sqrt{4 - x^2}.

In calculus, we learn that we can approximate integrals using Riemann sums. Interval arithmetic allows us to make these Riemann sums rigorous in a very simple way, as follows.

We split up the $x$ axis into intervals, for example of equal width:

function make_intervals(N=10)
    xs = linspace(0, 2, N+1)
    return [xs[i]..xs[i+1] for i in 1:length(xs)-1]
end

intervals = make_intervals()

10-element Array{ValidatedNumerics.Interval{Float64},1}:
 Interval(0.0, 0.2)                               
 Interval(0.19999999999999998, 0.4)               
 Interval(0.39999999999999997, 0.6000000000000001)
 Interval(0.6, 0.8)                               
 Interval(0.7999999999999999, 1.0)                
 Interval(1.0, 1.2000000000000002)                
 Interval(1.2, 1.4000000000000001)                
 Interval(1.4, 1.6)                               
 Interval(1.5999999999999999, 1.8)                
 Interval(1.7999999999999998, 2.0)

Given one of those intervals, we evaluate the function of interest:

II = intervals[1]

Interval(0.0, 0.2)

f(II)

Interval(1.9899748742132397, 2.0)

The result is an interval that is guaranteed to contain the true range of the function $f$ over that interval. So the lower and upper bounds of the intervals may be used as lower and upper bounds of the height of the box in a Riemann integral:

intervals = make_intervals(30)

p = plot(aspect_ratio=:equal)
for X in intervals
    Y = f(X)

    plot!(IntervalBox(X, Interval(0, Y.lo)), c=:blue, label="", alpha=0.1)
    plot!(IntervalBox(X, Interval(Y.lo, Y.hi)), c=:red, label="", alpha=0.1)
end

plot!(f, 0, 2)

p

Now we just sum up the areas:

N = 20
intervals = make_intervals(N)

width = 2/N
width * sum(√(4 - X^2) for X in intervals)

Interval(3.0284648797549782, 3.2284648797549846)

As we increase the number of sub-intervals, the approximation gets better and better:

setdisplay(:standard, sigfigs=5)

println("N \t area interval \t \t diameter")
for N in 50:50:1000
    intervals = make_intervals(N)
    area = (2/N) * sum(√(4 - X^2) for X in intervals)

    println("$N \t $area \t $(diam(area))")
end

N 	 area interval 	 	 diameter
50 	 [3.0982, 3.1783] 	 0.0800000000000165
100 	 [3.1204, 3.1605] 	 0.040000000000032454
150 	 [3.1276, 3.1543] 	 0.02666666666670814
200 	 [3.1311, 3.1512] 	 0.02000000000006308
250 	 [3.1332, 3.1493] 	 0.016000000000075065
300 	 [3.1346, 3.1481] 	 0.013333333333415354
350 	 [3.1356, 3.1472] 	 0.011428571428676815
400 	 [3.1364, 3.1465] 	 0.010000000000123688
450 	 [3.137, 3.146] 	 0.008888888889027502
500 	 [3.1374, 3.1455] 	 0.008000000000148333
550 	 [3.1378, 3.1452] 	 0.007272727272884527
600 	 [3.1381, 3.1449] 	 0.006666666666829357
650 	 [3.1384, 3.1446] 	 0.006153846154013376
700 	 [3.1386, 3.1444] 	 0.0057142857144931725
750 	 [3.1388, 3.1443] 	 0.005333333333562784
800 	 [3.139, 3.1441] 	 0.005000000000246363
850 	 [3.1391, 3.1439] 	 0.004705882353203794
900 	 [3.1393, 3.1438] 	 0.004444444444719142
950 	 [3.1394, 3.1437] 	 0.004210526316076102
1000 	 [3.1395, 3.1436] 	 0.004000000000294435

Some fun with Π in Julia

14 March 2017 | Simon Byrne, Luis Benet and David Sanders