Image courtesy of Cormullion, code here.

This post is available as a Jupyter notebook here.

by Simon Byrne

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

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

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`

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

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

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

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

s, but the result is hampered by the slow convergence of the series.

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

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