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 4.141592653589793

Float32 ( 1 ) + pi 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 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:

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

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 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 ⁡ ( π 6 ) = 1 3 . \tan\left( \frac{\pi}{6} \right) = \frac{1}{\sqrt{3}}. tan ( 6 π ​ ) = 3 ​ 1 ​ .

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

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

series1 = 6 atan( 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 / 3 1/\sqrt{3} 1/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 = 6 atan( Taylor1( BigFloat , 99 ) ) 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:

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

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

π 4 = 4 arctan ⁡ ( 1 5 ) − arctan ⁡ ( 1 239 ) . \frac{\pi}{4} = 4 \arctan\left(\frac{1}{5}\right) - \arctan\left(\frac{1}{239}\right). 4 π ​ = 4 arctan ( 5 1 ​ ) − arctan ( 2 3 9 1 ​ ) .

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 : = ∑ n = 1 ∞ 1 n 2 , S := \sum_{n=1}^\infty \frac{1}{n^2}, S : = n = 1 ∑ ∞ ​ n 2 1 ​ ,

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

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

We will evalute S N S_N SN​ numerically, and use the following analytical bound for T N T_N TN​:

1 N + 1 ≤ T N ≤ 1 N \frac{1}{N+1} \le T_N \le \frac{1}{N} N + 1 1 ​ ≤ T N ​ ≤ N 1 ​

.

This is obtained by approximating the sum in T N T_N TN​ using integrals from below and above:

∫ x = N + 1 ∞ 1 x 2 d x ≤ T N ≤ ∫ x = N ∞ 1 x 2 d x . \int_{x=N+1}^\infty \frac{1}{x^2} dx \le T_N \le \int_{x=N}^\infty \frac{1}{x^2} dx. ∫ x = N + 1 ∞ ​ x 2 1 ​ d x ≤ T N ​ ≤ ∫ x = N ∞ ​ x 2 1 ​ d x .

S N S_N SN​ 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 S_N SN​, 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 T_N TN​ to obtain the bounds on S S 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) 6

N = 10000 S = forward_sum(N, Interval) S += 1 /(N+ 1 ) .. 1 /N π_interval = √( 6 S) [ 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) Interval( 0 ) + 1 / 3 ^ 2 Interval( 0.1111111111111111 , 0.11111111111111112 )

This is an interval containing the true real number 1 / 9 1/9 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 r r is A ( r ) = π r 2 A(r) = \pi r^2 A(r)=πr2. We will calculate the area of one quadrant of a circle of radius r = 2 r=2 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 r=2 r=2 is given by x 2 + y 2 = 2 2 = 4 x^2 + y^2 = 2^2 = 4 x2+y2=22=4, so

π = 1 4 A ( 2 ) = ∫ x = 0 2 y ( x ) d x = ∫ x = 0 2 4 − x 2 . \pi = \frac{1}{4} A(2) = \int_{x=0}^2 y(x) \, dx = \int_{x=0}^2 \sqrt{4 - x^2}. π = 4 1 ​ A ( 2 ) = ∫ x = 0 2 ​ y ( x ) d x = ∫ x = 0 2 ​ 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 x 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 f 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: