Select Page
This post originally appeared on the Julialang.org blog.

# Some fun with π in Julia

This post is available as a Jupyter notebook here

## π in Julia

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

## π via inline assembly instructions

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

1. the `pushq` and `movq` adds to the call stack frame.
2. `fldpi` pushes π to the x87 floating point register stack
3. `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.
4. `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
``````

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