Happy π Day, 2024¶

By Christopher Phan, 2024-03-14 / 2024-W11-4

Happy π Day! It's March 14, i.e. 3/14 using US date conventions, or 03-14 using ISO 8601 (or RFC 3339, if you prefer) and omitting the year. Back on π Day 2021, I made a JavaScript demo of a Monte Carlo method to approximate π. To celebrate this wonderful geek holiday in 2024, I will present another silly way to calculate an approximation for this wonderful irrational number.

About this document¶

This document is a Jupyter notebook, which means that it includes cells with computer code and the output. For example:

In [1]:
# This is computer code

using Printf
@printf "The value of π is approximately %0.5f." pi

# Just below is the output

The value of π is approximately 3.14159.

A copy of this notebook is available on GitHub.

Julia¶

In this notebook, I am using the language Julia, because it is relatively easy to work with, has just-in-time compilation, and has built-in arbitrary precision numbers. (Also fun fact, the "Ju" in "Jupyter" stands for Julia. The other parts stand for Python and R.)

The arbitrary precision numbers will help. The default in Julia (and many other languages) is to store integers as having a fixed number of bits, which makes computation faster. For example, 42 is a signed 64-bit integer.

In [2]:
typeof(42)
Out[2]:
Int64
In [3]:
# calculation of 21!
a = 1
for k in range(1, 21)
    a = a * k
    @printf "%i! = %i\n" k a
end

1! = 1
2! = 2
3! = 6
4! = 24
5! = 120
6! = 720
7! = 5040
8! = 40320
9! = 362880
10! = 3628800
11! = 39916800
12! = 479001600
13! = 6227020800
14! = 87178291200
15! = 1307674368000
16! = 20922789888000
17! = 355687428096000
18! = 6402373705728000
19! = 121645100408832000
20! = 2432902008176640000
21! = -4249290049419214848

Note that $21!$ is shown to be negative! This is because of how 64-bit signed integers are represented internally. Basically, we ran out of room in the 64 bits to calculate $20! \times 21$. The solution is to use "big" integers, which can have as many bits as necessary. (The trade-off is that they are a bit slower.)

In [4]:
typeof(big(21))
Out[4]:
BigInt
In [5]:
# correct calculation of 21!
a = big(1)
for k in range(1, 21)
    a = a * k
    @printf "%i! = %i\n" k a
end

1! = 1
2! = 2
3! = 6
4! = 24
5! = 120
6! = 720
7! = 5040
8! = 40320
9! = 362880
10! = 3628800
11! = 39916800
12! = 479001600
13! = 6227020800
14! = 87178291200
15! = 1307674368000
16! = 20922789888000
17! = 355687428096000
18! = 6402373705728000
19! = 121645100408832000
20! = 2432902008176640000
21! = 51090942171709440000

Julia has a built-in factorial function:

In [6]:
factorial(big(21))
Out[6]:
51090942171709440000

As it happens, I don't have much experience writing code in Julia (I have mainly used Python, Rust, and Javascript), so this a good excuse to practice my Julia skills.

Cosine (Taylor's version)¶

For this exercise, I will use the Taylor series of cosine centered at zero. Note that $$\frac{d^n}{dx^n} \cos x = \begin{cases} (-1)^{(n + 1)/2} \sin x &\text{ if $n$ is odd,}\cr (-1)^{n/2} \cos x &\text{ if $n$ is even.}\end{cases}$$ At $x = 0$, we have $$\left. \frac{d^n}{dx^n} \cos x \right|_{x = 0} = \begin{cases} 0 &\text{ if $n$ is odd,}\cr (-1)^{n/2} &\text{ if $n$ is even.}\end{cases}$$ Therefore, "Taylor's version" of cosine is $$\cos x = \sum_{n = 0}^\infty \left. \frac{d^n}{dx^n} \cos x \right |_{x = 0} \frac{x^n}{n!} = \sum_{n = 0}^\infty (-1)^n \frac{x^{2n}}{(2n)!}.$$

(I really wish I was was clever enough to have come up with this "Taylor's version" joke, but alas, it's just something I saw on social media.)

We'll set $$T_n(x) = \frac{x^{2n}}{(2n)!}$$ and $$C_n(t) = \sum_{k = 0}^{n-1} (-1)^k T_k(x),$$ which is a (Taylor polynomial) approximation of $\cos x$.

We can implement this in Julia:

In [7]:
# Note: functions automatically return last value

function partial_cos_term(n, x) # T_n(x)
    x^(2 * n)//factorial(2 * n)
end 

# NOTE: Julia's range(0, N) returns values 0 through N inclusive
# (unlike Python)

function partial_cos(n, x) # C_n(x)
    sum([(-1)^(k) * partial_cos_term(k, x) for k in range(0, n - 1)])    
end
Out[7]:
partial_cos (generic function with 1 method)

The // operator in Julia represents exact (rational) division, e.g.

In [8]:
1//2 + 2//3
Out[8]:
7//6

Squeezing $\cos x$¶

(The material in this section was inspired by an argument that appears on p. 3 of Real and Complex Analysis, 3rd Edition by Walter Rudin (McGraw Hill 1987), often affectionaltely known as "Papa Rudin". Fun fact: Rudin's house is architectually signifcant enough to have its own article in Wikipedia.)

Note that if $0 < x < 2k + 1$, then $$x^{2k + 2} < (2k + 1)(2k + 2) x^{2k}.$$

So, $$T_{k + 1}(x) = \frac{x^{2k + 2}}{(2k + 2)!} < \frac{(2k + 1)(2k + 2)x^{2k}}{(2k + 2)!} = T_k(x).$$

This gives us some inequalities for $C_n(x)$: For any $n \geq 1$ and $0 < x \leq 2n + 1$, we have $$C_{2n + 2}(x) = C_{2n}(x) + T_{2n}(x) - T_{2n + 1}(x) > C_{2n}(t)$$ and $$C_{2n + 3}(x) = C_{2n + 1}(x) - T_{2n + 1}(x) + T_{2n + 2}(x) < C_{2n + 1}(x).$$

Therefore, for every $0 < x \leq 2$, we have $$C_2(x) < C_4(x) < C_6(x) < \dots < \cos t < \dots < C_7(x) < C_5(x) < C_3(x).$$

We can illustrate this with some plots:

In [9]:
using Plots
using LaTeXStrings

x = [2 * big(x) // 100 for x in range(1, 100)]
cosx = cos.(float.(x)) # func. applies func elementwise
y = [
    log10.(float.(partial_cos.(big(3), x)) - cosx),    
    log10.(float.(partial_cos.(big(5), x)) - cosx),
    log10.(float.(partial_cos.(big(7), x)) - cosx),
    log10.(float.(partial_cos.(big(9), x)) - cosx),
]
plot(x, y, 
    title="Plot of \$\\log_{10}(C_n(x) - \\cos(x))\$ for odd \$n\$",
    label=["n = 3" "n = 5" "n = 7" "n = 9"]
)
Out[9]:
In [10]:
y = [
    log10.(cosx - float.(partial_cos.(big(2), x))),
    log10.(cosx - float.(partial_cos.(big(4), x))),    
    log10.(cosx - float.(partial_cos.(big(6), x))),
    log10.(cosx - float.(partial_cos.(big(8), x))),  
]
plot(x, y, 
    title="Plot of \$\\log_{10}(\\cos(x) - C_n(x))\$ for even \$n\$",
    label=["n = 2" "n = 4" "n = 6" "n = 8"]
)
Out[10]:

Furthermore, remember that for any real $x$, $\lim_{n \rightarrow \infty} C_n(x) = \cos x$. Therefore, for any real $0 < x \leq 2$:

  • If $\cos x > 0$, then there exists an integer $N > 0$ such that $C_{2n}(x) > 0$ for $n \geq N$.
  • If $\cos x < 0$, then there exists an integer $N > 0$ such that $C_{2n + 1}(x) < 0$ for $n \geq N$.

Approximating $\pi/2$ means finding approximate solutions to $\cos x = 0$ where $1 \leq x \leq 2$¶

What does this have to do with approximating $\pi$? Well, one definition of $\pi$ (e.g. in the aforementioned Papa Rudin) is that $\pi$ is the smallest positive real number such that $\cos(\pi/2) = 0$.

Note that $\cos 0 = 1$, and that when $0 < x \leq 1$, we have $$0 < 1 - \frac{x^2}{2} = C_2(x) < \cos x,$$ and thus, $\cos x > 0$ when $0 \leq x \leq 1$.

Therefore, $\pi/2 > 1$.

Now, $-\sin x < 0$ when $0 < x < \pi$, so that means $\cos$ is strictly decreasing on the interval $[1, 2]$. As $$\cos 2 < C_3(2) = 1 - \frac{4}{2} + \frac{16}{24} = -\frac{1}{3} < 0,$$ we know that $\cos$ has a unique zero on $[1, 2]$, i.e. $1 \leq \pi/2 \leq 2$.

A way to tighten bounds¶

So far, we've narrowed down the value of $\pi/2$ to the interval $[1, 2]$. How do we make a bound better.

For $1 \leq a < b \leq 2$ and $n \geq 1$, let $P(n, a, b)$ be the statement

$C_{2n}(a) > 0$ and $C_{2n + 1}(b) < 0$.

Note that:

  • We showed above that $P(1, 1, 2)$ is true.
  • $P(n, a, b)$ implies $\cos a > 0$ and $\cos b < 0$, i.e. that $a < \pi/2 < b$.

Now, assume that $P(n, a, b)$ is true for some rational $1 \leq a < b \leq 2$ and $n \geq 1$. Let $c = (a + b)/2$. Exactly one of the following conditions is true:

  1. $C_{2n + 3}(c) < 0$. Then $P(n + 1, a, c)$ is true.
  2. $C_{2n + 3}(c) > 0$ and $C_{2n + 2}(c) > 0$. Then $P(n + 1, c, b)$ is true.
  3. $C_{2n + 3}(c) > 0$ and $C_{2n + 2}(c) < 0$. In this case, $P(n + 1, a, b)$ is true.

In cases (1) and (2), we've sucessfully cut in half the possible values of $\pi/2$. Case (3) seems like let-down, but remember that for $n$ sufficently large, either $C_{2n} > 0$ or $C_{2n + 1} < 0$. So by incrementing $n$, we will eventually land in (1) or (2).

By repeatedly iterating this process, starting with $n = 1$, $a = 1$ and $b = 2$, we can find $\pi/2$ to any desired acccuracy.

Here is the Julia implementation of this procedure:

In [11]:
function narrow(n, a, b)
    c = (a + b) // 2
    if partial_cos(2 * big(n) + 3, big(c)) < 0
        (a, c)
    elseif partial_cos(2 * big(n) + 1, big(c)) > 0
        (c, b)
    else
        (a, b)
    end
end
Out[11]:
narrow (generic function with 1 method)

Working through the first few steps¶

Let's work through the first few steps. We start with $P(1, 1, 2)$.

In [12]:
narrow(1, 1, 2)
Out[12]:
(3//2, 2)

So, we know that $P(2, 3/2, 2)$ is true. We've narrowed the value of $\pi/2$ to $[3/2, 2]$. We repeat:

In [13]:
narrow(2, 3//2, 2)
Out[13]:
(3//2, 7//4)
In [14]:
narrow(3, 3//2, 7//4)
Out[14]:
(3//2, 13//8)
In [15]:
narrow(4, 3//2, 13//8)
Out[15]:
(25//16, 13//8)
In [16]:
narrow(5, 25//16, 13//8)
Out[16]:
(25//16, 51//32)
In [17]:
narrow(6, 25//16, 51//32)
Out[17]:
(25//16, 101//64)

Approximating $\pi$¶

Now it's time to automate this.

One last thing before seeing the conclusion: when we determine that $P(n, a, b)$ is true, then we are approximating $\pi/2 \approx (a + b)/2$. But this is the same thing as approximating $\pi \approx a + b$.

I will use some apply narrow repeatedly to show some $P(80, u, v)$. I will use some Julia code to make a table showing, for each $P(n, a, b)$ obtained:

  1. $n$
  2. The numerator of the rational approximation $a + b$ of $\pi$
  3. The denominator of the rational approximation $a + b$ of $\pi$
  4. A decimal representation of $a + b$
  5. $\log_{10}(b - a)$, which gives you an idea how many decimal places of precision we have achieved.
In [18]:
largest_k = big(80)

# create an array of strings
table = Array{String, 2}(undef, largest_k + 1, 5)

# create header row
# Everything is 1-indexed in Julia. 
# It messes me up sometimes!
table[1, 1] = "n"
table[1, 2] = "numerator"
table[1, 3] = "denominator"
table[1, 4] = "decimal"
table[1, 5] = "log_10(b-a)"

for k in range(1, 5)
    table[2, k] = "----"
end

# fill in the numbers in the table
lower = big(1)
upper = big(2)
k = big(2)
while k < largest_k + 1
    pi_midpt = lower + upper # approximation for pi
    pi_flt = float(pi_midpt) # as a float
    lg_err = log10(float(upper - lower))
    # OK, so because (lower, upper) is an interval for pi/2
    # (2 * lower, 2 * upper) is an interval for pi
    # Hence |lower + upper - pi| < upper - lower.
    # Let err = upper - lower.
    # Now, to guarantee |lower + upper - pi| to be precise
    # to N digits, we need:
    #    err < 0.5 * 10^(-N)
    #  2*err < 10^(-N)
    # log10(2) + log10(err) < -N
    # N < -log10(2) - log10(err) 
    # 
    # So, we'll display the larger of 
    # floor(-log10(2) - log10(err)) or 2 digits
    # past the decimal point.
    num_digits = max(Int(floor(-log10(2) - lg_err)), 2)

    # populate table row
    table[k + 1, 1] = @sprintf("%2i", k)
    table[k + 1, 2] = @sprintf("%i", numerator(pi_midpt))
    table[k + 1, 3] =  @sprintf("%i", denominator(pi_midpt))
    table[k + 1, 4] = @sprintf("%.*f", num_digits, pi_flt)
    table[k + 1, 5] = @sprintf("%.5f", lg_err)
    k += 1
    if k <= largest_k # otherwise the loop is ending
        (lower, upper) = narrow(k, lower, upper)
    end
end

# determine the maximum width of the text in each column
max_widths = [
    maximum([length(table[k, j]) 
    for k in range(1, largest_k + 1)]) for j in range(1, 5)
]

# the second row of the table will have entries like 
# "-------:" or ":-------" (markdown syntax
# for right and left justification, respectively)
for k in range(1, 3)
    # String concatenation is represented by * in Julia, not +
    # This makes WAY more sense; ask your friendly local algebraist why.
    table[2, k] = "-"^(max_widths[k] - 1) * ":"
end
for k in range(4, 5)
    table[2, k] = ":" * "-"^(max_widths[k] - 1)
end
for k in range(1, largest_k + 1)
    # We use some printf spacing goodness to make everything
    # look really nice
    for j in range(1, 3)
        @printf "| %*s " max_widths[j] table[k, j] 
    end
    for j in range(4, 5)
        @printf "| %-*s " max_widths[j] table[k, j]
    end
    print("|\n")
end

|    n |                numerator |              denominator | decimal                   | log_10(b-a) |
| ---: | -----------------------: | -----------------------: | :------------------------ | :---------- |
|    2 |                        3 |                        1 | 3.00                      | 0.00000     |
|    3 |                        7 |                        2 | 3.50                      | -0.30103    |
|    4 |                       13 |                        4 | 3.25                      | -0.60206    |
|    5 |                       25 |                        8 | 3.12                      | -0.90309    |
|    6 |                       51 |                       16 | 3.19                      | -1.20412    |
|    7 |                      101 |                       32 | 3.16                      | -1.50515    |
|    8 |                      201 |                       64 | 3.14                      | -1.80618    |
|    9 |                      403 |                      128 | 3.15                      | -2.10721    |
|   10 |                      805 |                      256 | 3.14                      | -2.40824    |
|   11 |                     1609 |                      512 | 3.14                      | -2.70927    |
|   12 |                     3217 |                     1024 | 3.14                      | -3.01030    |
|   13 |                     6433 |                     2048 | 3.141                     | -3.31133    |
|   14 |                    12867 |                     4096 | 3.141                     | -3.61236    |
|   15 |                    25735 |                     8192 | 3.141                     | -3.91339    |
|   16 |                    51471 |                    16384 | 3.142                     | -4.21442    |
|   17 |                   102943 |                    32768 | 3.1416                    | -4.51545    |
|   18 |                   205887 |                    65536 | 3.1416                    | -4.81648    |
|   19 |                   411775 |                   131072 | 3.1416                    | -5.11751    |
|   20 |                   823549 |                   262144 | 3.14159                   | -5.41854    |
|   21 |                  1647099 |                   524288 | 3.14159                   | -5.71957    |
|   22 |                  3294199 |                  1048576 | 3.14159                   | -6.02060    |
|   23 |                  6588397 |                  2097152 | 3.141593                  | -6.32163    |
|   24 |                 13176795 |                  4194304 | 3.141593                  | -6.62266    |
|   25 |                 26353589 |                  8388608 | 3.141593                  | -6.92369    |
|   26 |                 52707179 |                 16777216 | 3.141593                  | -7.22472    |
|   27 |                105414357 |                 33554432 | 3.1415927                 | -7.52575    |
|   28 |                210828715 |                 67108864 | 3.1415927                 | -7.82678    |
|   29 |                421657429 |                134217728 | 3.1415927                 | -8.12781    |
|   30 |                843314857 |                268435456 | 3.14159266                | -8.42884    |
|   31 |               1686629713 |                536870912 | 3.14159265                | -8.72987    |
|   32 |               3373259427 |               1073741824 | 3.14159265                | -9.03090    |
|   33 |               6746518853 |               2147483648 | 3.141592654               | -9.33193    |
|   34 |              13493037705 |               4294967296 | 3.141592654               | -9.63296    |
|   35 |              26986075409 |               8589934592 | 3.141592654               | -9.93399    |
|   36 |              53972150819 |              17179869184 | 3.141592654               | -10.23502   |
|   37 |             107944301637 |              34359738368 | 3.1415926536              | -10.53605   |
|   38 |             215888603273 |              68719476736 | 3.1415926536              | -10.83708   |
|   39 |             431777206545 |             137438953472 | 3.1415926536              | -11.13811   |
|   40 |             863554413089 |             274877906944 | 3.14159265359             | -11.43914   |
|   41 |            1727108826179 |             549755813888 | 3.14159265359             | -11.74017   |
|   42 |            3454217652357 |            1099511627776 | 3.14159265359             | -12.04120   |
|   43 |            6908435304715 |            2199023255552 | 3.141592653590            | -12.34223   |
|   44 |           13816870609431 |            4398046511104 | 3.141592653590            | -12.64326   |
|   45 |           27633741218861 |            8796093022208 | 3.141592653590            | -12.94429   |
|   46 |           55267482437723 |           17592186044416 | 3.141592653590            | -13.24532   |
|   47 |          110534964875445 |           35184372088832 | 3.1415926535898           | -13.54635   |
|   48 |          221069929750889 |           70368744177664 | 3.1415926535898           | -13.84738   |
|   49 |          442139859501777 |          140737488355328 | 3.1415926535898           | -14.14841   |
|   50 |          884279719003555 |          281474976710656 | 3.14159265358979          | -14.44944   |
|   51 |         1768559438007111 |          562949953421312 | 3.14159265358979          | -14.75047   |
|   52 |         3537118876014221 |         1125899906842624 | 3.14159265358979          | -15.05150   |
|   53 |         7074237752028441 |         2251799813685248 | 3.141592653589794         | -15.35253   |
|   54 |        14148475504056881 |         4503599627370496 | 3.141592653589793         | -15.65356   |
|   55 |        28296951008113761 |         9007199254740992 | 3.141592653589793         | -15.95459   |
|   56 |        56593902016227523 |        18014398509481984 | 3.141592653589793         | -16.25562   |
|   57 |       113187804032455045 |        36028797018963968 | 3.1415926535897933        | -16.55665   |
|   58 |       226375608064910089 |        72057594037927936 | 3.1415926535897932        | -16.85768   |
|   59 |       452751216129820177 |       144115188075855872 | 3.1415926535897932        | -17.15871   |
|   60 |       905502432259640355 |       288230376151711744 | 3.14159265358979324       | -17.45974   |
|   61 |      1811004864519280711 |       576460752303423488 | 3.14159265358979324       | -17.76077   |
|   62 |      3622009729038561421 |      1152921504606846976 | 3.14159265358979324       | -18.06180   |
|   63 |      7244019458077122843 |      2305843009213693952 | 3.141592653589793239      | -18.36283   |
|   64 |     14488038916154245685 |      4611686018427387904 | 3.141592653589793239      | -18.66386   |
|   65 |     28976077832308491369 |      9223372036854775808 | 3.141592653589793238      | -18.96489   |
|   66 |     57952155664616982739 |     18446744073709551616 | 3.141592653589793238      | -19.26592   |
|   67 |    115904311329233965479 |     36893488147419103232 | 3.1415926535897932385     | -19.56695   |
|   68 |    231808622658467930957 |     73786976294838206464 | 3.1415926535897932385     | -19.86798   |
|   69 |    463617245316935861913 |    147573952589676412928 | 3.1415926535897932385     | -20.16901   |
|   70 |    927234490633871723825 |    295147905179352825856 | 3.14159265358979323846    | -20.47004   |
|   71 |   1854468981267743447651 |    590295810358705651712 | 3.14159265358979323846    | -20.77107   |
|   72 |   3708937962535486895301 |   1180591620717411303424 | 3.14159265358979323846    | -21.07210   |
|   73 |   7417875925070973790601 |   2361183241434822606848 | 3.141592653589793238462   | -21.37313   |
|   74 |  14835751850141947581203 |   4722366482869645213696 | 3.141592653589793238463   | -21.67416   |
|   75 |  29671503700283895162407 |   9444732965739290427392 | 3.141592653589793238463   | -21.97519   |
|   76 |  59343007400567790324813 |  18889465931478580854784 | 3.141592653589793238463   | -22.27622   |
|   77 | 118686014801135580649625 |  37778931862957161709568 | 3.1415926535897932384626  | -22.57725   |
|   78 | 237372029602271161299249 |  75557863725914323419136 | 3.1415926535897932384626  | -22.87828   |
|   79 | 474744059204542322598499 | 151115727451828646838272 | 3.1415926535897932384626  | -23.17931   |
|   80 | 949488118409084645196999 | 302231454903657293676544 | 3.14159265358979323846265 | -23.48034   |

The approximation of $\pi$ we found is $\pi \approx 3.14159265358979323846265$.

Thank you¶

I hope you have enjoyed this π Day journey. I hope you are able to celebrate π Day with a yummy treat as well.