# Some fun with π in Julia

*This post is available as a Jupyter notebook here*

## π 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 `BigFloat`

s, 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 `Irrational`

s is that inequalities are correct:

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

```
true
```

## π via inline assembly instructions

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 (“**f**loating point **l**oa**d** **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

*(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

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 `BigFloat`

s:

```
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*:

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

### Machin’s approach

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

begin{equation} frac{pi}{4} = 4 arctanleft(frac{1}{5}right) – arctanleft(frac{1}{239}right). end{equation}

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 π

*(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

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$:

.

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

$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 principal, we could attain arbitrarily good precision with higher-precision `BigFloat`

s, 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)
```

It’s 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="")
```

0.5 1.0 1.5 0.0 0.5 1.0 1.5

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

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
```

0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 y61

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
```