Monty Carlo Pi Approximation

Topic: Maths
Date: 03/06/2017

What is PI?

To put it simply, PI is the ratio of a circles circumference to its diameter. We can find the irrational constant with the following calculation. Assuming we have a circle with a diameter $D$, we then measure the circumference $C$ of the circle. $\pi$ can then be calculated by $C/D$.

$$ \begin{aligned} let \quad D &= 5, \\ C &= 15.72 \\ \therefore \quad \pi &\approx \frac{C}{D} = \frac{15.72}{5} = 3.144 \end{aligned} $$

Its not particularly accurate, but neither was the measurement I got for the circumference. Now we know how to calculate PI by 2 measurable quantities, it doesn't take much to see why the circumference of a circle is found by the formula $C = 2 \pi r = \pi D$. So if we had a circle with a diameter of 1, it's circumference would be equal to $\pi$.

As a side note, we can find the equation for the area of a circle by integrating the circumference equation (integration gives the area under a curve).

$$ \begin{aligned} C &= \pi D = 2\pi r\\ A &= \int 2\pi r \, dr = \pi r^2 \end{aligned} $$

Approximating PI

While the above method will give us a value for $\pi$, it's not exactly easy/possible to calculate via a computer, i.e., how do you measure the circumference of a circle draw on a screen automatically? The answer is... you don't. There are many other ways of calculating the decimal places of $\pi$. The Gauss-Legendre algorithm is a particularly interesting example of how to calculate $\pi$ (more can be found here). It's an iterative algorithm, whereby each pass through doubles the precision for the value of $\pi$. Since I write code in Erlang for a living, I have also provided an example solution for the algorithm here. It is worth noting, that running the algorithm past 3 iterations will just return the same number. This is because of the limit to floating point numbers in Erlang.

Monty Carlo Method

For this post however, I am going to use a pretty basic Monty Carlo method. The idea is to place a unit circle inside a 1x1 square. We can then get a formula for $\pi$ based off the ratio of the circle's area and the square's area. $A_c$ is the area of the circle, $A_s$ is the area of the square.

$$ \begin{aligned} A_c &= \pi r^2 \\ A_s &= (2r)^2 = 4r^2 \\ \frac{A_c}{A_s} &= \frac{\pi r^2}{4r^2} = \frac{\pi}{4} \end{aligned} $$

Given the above, we can see that $\pi$ can be calculated as $4 \times \frac{A_c}{A_s}$. The simulation below uses random sampling to produce a value for this ratio thus letting us approximate $\pi$ (albeit very badly).



$\pi \approx$ 4 * Inside / All

Accuracy = | Math.PI - PI Approximation |


Gauss-Legendre algorithm | Other algorithms | Monty Carlo Pi Approximation | Gauss-Legendre (Erlang Implentation)

Concrete Mathematics: Chapter 2: Part 3

Topic: Maths
Date: 07/10/2016

Manipulation of Sums

As the book covers the distributive and associative laws very well, this section will only focus on the commutative law. The book states that the sum of integers in the finite set $K$ can be transformed with a permutation $p(k).$ The law looks as follows

$$ \sum \limits_{k \in K} a_k = \sum \limits_{p(k) \in K} a_{p(k)}. $$

Put simply, this means that the sum remains the same as long as all the elements in the set get summed, just in a different permutation. A simple example of this is $p(k) = -k$, which is the same set just in reverse. The book gives an example of the set $K = \{ -1, 0, 1 \}$, so the permutation $p(k) = -k$ gives

$$ a_{-1} + a_0 + a_1 = a_1 + a_0 + a_{-1} $$

While this example isn't too taxing, the following is a little harder to understand. This next example is how Gauss's trick from Chapter 1 can be found using the different manipulation laws. The first step in the example is the most confusing as it is the one which features the commutative law. We start with

$$ S = \sum \limits_{0 \le k \le n}(a + bk). $$

The book says that we now replace $k$ with $n - k.$ Sadly it does little to explain why this is the case. Well, we can show the permutation more clearly by creating a new variable called $k_1$ and saying $k_1 = n - k$, which we can rearrange for $k$ and plug it back into the sum.

$$ \sum \limits_{0 \le k \le n}(a + bk) = \sum \limits_{0 \le n - k_1 \le n}(a + b(n - k_1)) $$

Now it pays to look at the bounds of the sum, $0 \le n - k_1 \le n$ is technically the same as $0 \le k_1 \le n$ (We know this to be true as $n - k_1$ is just the reverse of $k_1$). So we can substitute this back into the summation.

$$ S = \sum \limits_{0 \le k_1 \le n}(a + b(n - k_1)) $$

Now as $k_1$ is just a simple variable, we can replace its letter with $k$. This gives us

$$ S = \sum \limits_{0 \le k \le n}(a + b(n - k)) $$

Just to prove that these are the same, we can expand out the sum to show the individual components for each summation.

$$ \begin{aligned} \sum \limits_{0 \le k \le n}(a + bk) &= (a + b \cdot 0) + (a + b \cdot 1) + \ldots + (a + b \cdot n) \\ \sum \limits_{0 \le k \le n}(a + b(n - k)) &= (a + b \cdot (n - 0)) + (a + b \cdot (n - 1)) + \ldots + (a + b \cdot (n - n)) \end{aligned} $$

You can see that these are essentially the same just with the later summing down. The book now walks us through the other laws with the same example, but those are not too hard to understand, so I won't be covering them here.

Links with Integration

The book suggests we can use the concept of integrals to solve sums. This is first introduced by showing an example in method 4 of chapter 2.5. Unfortunately, it doesn't quite work out as it leaves a small error term. Fortunately, the book's next section shows the real connection between integration and summation.

It is worth noting that the book uses different notation than some people might be used to, e.g $D(f(x)) = \frac{d}{dx} f(x).$ The equation for a derivative that the book uses is more commonly known as

$$ \frac{d}{dx} f(x) = \lim_{\Delta x \to 0} \frac{f(x + \Delta x) - f(x)}{\Delta x}. $$

Now the book introduces finite calculus. This is the same as infinite calculus (integration/differentiation) but restricted to integers. This means that $\Delta x = 1$ is the lowest value we can get to with a limit approaching 0 using integers. Therefore the derivative for finite calculus is

$$ \Delta f(x) = f(x + 1) - f(x). $$

This isn't super useful yet, as it just turns equations into sums. Take $f(x) = x^2$ for example

$$ \begin{aligned} \Delta f(x) &= (x + 1)^2 - x^2 \\ &= (x + 1)(x + 1) - x^2 \\ &= x^2 + 2x + 1 - x^2 \\ &= 2x + 1. \end{aligned} $$

Equation 2.48 tells us that this means

$$ \sum \limits_{k = 0}^{x - 1} 2k + 1 = x^2 $$

Finite 'integration' to remove summation is corollary to this delta operation the book has defined. Unfortunately, it doesn't correlate 1:1. Instead we need to use falling powers to remove the summation cleanly.

Falling Powers

Falling powers/factorials are a handy dandy way to keep the existing rules of infinite calculus but for finite calculus instead. A falling power is explained by equation 2.43, but I prefer to have examples:

$$ \begin{aligned} x^{\underline{3}} &= x(x - 1)(x - 2) \\ x^{\underline{2}} &= x(x - 1) \\ x^{\underline{1}} &= x \\ x^{\underline{0}} &= 1 \\ x^{\underline{-1}} &= \frac{1}{x + 1} \\ x^{\underline{-2}} &= \frac{1}{(x + 1)(x + 2)} \\ x^{\underline{-3}} &= \frac{1}{(x + 1)(x + 2)(x + 3)}. \end{aligned} $$

Falling powers can be used with the rule

$$ \Delta (x^{\underline{m}}) = mx^{\underline{m - 1}} $$

So we can solve sums like

$$ \begin{aligned} \sum \limits_{0 \le k < n} k^{\underline{m}} &= \frac{k^{\underline{m + 1}}}{m + 1} = \frac{n^{\underline{m + 1}}}{m + 1} \\ \sum \limits_{0 \le k < n} k &= \frac{n^{\underline{2}}}{2} = \frac{n(n - 1)}{2} \end{aligned} $$

Great! That's the meat of this chapter, anything else isn't too hard to pick up after understanding this.

Understanding Limits

Here I'm just going to briefly mention the limits shown in the infinite sums section of the chapter. The first issue I had was with a limit between two bounds (I forgot to subtract the equation evaluated at 0). This was with the summation

$$ \sum \limits_{k \ge 0} \frac{1}{(k + 1)(k + 2)}. $$

Next the summation is solved by the following steps

$$ \begin{aligned} \sum \limits_{k \ge 0} \frac{1}{(k + 1)(k + 2)} &= \sum \limits_{k \ge 0} k^{\underline{-2}} \\ &= \lim_{n \to \infty} \sum \limits^{n - 1}_{k = 0}k^{\underline{-2}} \\ &= \lim_{n \to \infty} \frac{k^{\underline{-1}}}{-1} \bigg|^n_0 \end{aligned} $$

Now we need to evaluate the limits at $n$ and $0$

$$ \begin{aligned} \lim_{n \to \infty} \frac{n^{\underline{-1}}}{-1} - \lim_{n \to \infty} \frac{0^{\underline{-1}}}{-1} &= \lim_{n \to \infty} -\frac{1}{k + 1} - \lim_{n \to \infty}\left(- \frac{1}{0 + 1} \right) \\ &= 0 - (-1) \\ &= 1 \end{aligned} $$

Which is the same answer the book gives. So I'm happy.


Difference Quotient (Clip 5) | Commutative law

Concrete Mathematics: Chapter 2: Part 2

Topic: Maths
Date: 29/09/2016

Summation Factor

As mentioned in the previous blog, I need to go over a few things that were mentioned but never explained. The first is the summation factor. This is what is used in order to reduce a recursive equation down to a more simple version. So we can reduce a recursion from

$$ a_nT_n = b_nT_{n-1} + c_n $$

to the form

$$ S_n = S_{n-1} + s_nc_n $$

by multiplying by a summation factor and then substituting in $S_n = s_na_nT_n.$ This greatly reduces the recurrence so it is much easier to convert it into a summation with sigma notation. The above would then translate to the following summation:

$$ S_n = s_0a_0T_0 + \sum \limits^n_{k=1}s_kc_k = s_1a_1T_0 + \sum \limits^n_{k=1}s_kc_k $$

It is important to note that there is the term $s_0a_0T_0$ as this is used when the starting terms are non-zero or the summation factor is undefined. The book makes subtle use of this term with the quicksort example. I will talk a bit about that later on.

Now we know how to use the summation factor, we need to learn how to get the summation factor. In the first part for this chapter I just multiplied the tower of Hanoi recursion by $2^{-n}$ without explaining where that was from. Well, it turns out that was found using the equation

$$ s_n = \frac{a_{n-1}a_{n-2} \ldots a_1}{b_nb_{n-1} \ldots b_2} $$

Okay. Let's try that and see if we can get $2^{-n}.$

First we note that $a_n=1$, $b_n=2$ and $c_n=1$. This means that the formula is the following

$$ \begin{aligned} s_n &= \frac{1 \cdot 1 \cdot \ldots \cdot 1}{2 \cdot 2 \cdot \ldots \cdot 2} &(\text{Substitute in the values}) \\ s_n &= \frac{1}{2^{n-1}} &(\text{There are}\ n-1\ \text{terms of}\ b_n) \\ s_n &= 2^{-(n-1)} &(\text{Rearrange}) \end{aligned} $$

That's weird, we ended up with a different summation factor than the one that's in the book. Not to worry, the book tells us that it's ok to use any multiple of the summation factor and it should still work. So let's carry on.

$$ \begin{aligned} T_0 &= 0; \\ T_n &= 2T_{n-1} + 1, &\text{for}\ n>0 \end{aligned} $$

Multiply by summation factor, $2^{-(n-1)}$

$$ \begin{aligned} T_0/2^{n-1} &= 0 \\ T_n/2^{n-1} &= 2T_{n-1}/2^{n-1} + 1/2^{n-1} \\ &= T_{n-1}/2^{n-2} + 2^{-(n-1)} \\ \end{aligned} $$

Let $S_n = T_n/2^{n-1}$

$$ \begin{aligned} S_n &= S_{n-1} + 2^{-(n-1)} \\ &= \sum \limits^n_{k=1}2^{-(k-1)} \end{aligned} $$

Convert back to $T_n$

$$ \begin{aligned} T_n/2^{n-1} &= \sum \limits^n_{k=1}2^{-(k-1)} \\ T_n &= 2^{n-1} \sum \limits^n_{k=1}2^{-(k-1)} \end{aligned} $$

Now we can try some values to verify that this will indeed work. First $n=3,$

$$ \begin{aligned} T_n &= 2^2(2^{-0} \cdot 2^{-1} \cdot 2^{-2}) \\ &= 4\left(1 \cdot \frac{1}{2} \cdot \frac{1}{4}\right) \\ &= 7 \end{aligned} $$

Great! Checking other $n$ shows that this is a summation solution for the tower of Hanoi. However it is not as elegant as the solution the book gave when it used $2^{-n}$ as the summation factor. So make sure to play with the summation factor when converting recurrences to sums.

Problem with quicksort example

So far this book has stumped me in a few places, but I can generally work through the problem until the answer clicks. However, the one I have spent the most time on was a problem with the quicksort example.

While working through their examples, I couldn't see where they gained the extra term for the final solution. Fortunately, I found the solution on math.stackexchange. I turns out that the summation factor is undefined for the early terms and so we need to add them to the summation which doesn't include the terms.

$$ \begin{aligned} S_n &= S_2 + \sum_{k=3}^ns_kc_k \\ &= s_2a_2C_2 + \sum_{k=3}^ns_kc_k\\ &= \frac13 \cdot 2 \cdot 3 + \sum_{k=3}^n \frac{4k}{k(k+1)} \\ &= 2 + 4 \sum_{k=3}^n \frac1{k+1} \end{aligned} $$

Then convert back to $C_n$

$$ \begin{aligned} C_n &= \frac{S_n}{s_na_n}\\ &= \frac{n+1}2 \left(2 + 4 \sum_{k=3}^n \frac1{k+1} \right) \\ &= (n+1) \left(1 + 2 \sum_{k=1}^n \frac1{k+1}-2 \sum_{k=1}^2 \frac1{k+1} \right) \\ &= (n+1) \left(1 + 2 \sum_{k=1}^n \frac1{k+1} - \frac53 \right)\\ &= (n+1) \left(2 \sum_{k=1}^n \frac1{k+1} - \frac23 \right) \\ &= 2(n+1) \sum_{k=1}^n \frac1{k+1} - \frac23(n+1) \end{aligned} $$

Many thanks to Brian M. Scott for answering the question posed by Cyril on StackExchange. I'd still be trying to figure it out otherwise.


Summation factor problem | Quicksort problem

Concrete Mathematics: Chapter 2: Part 1

Topic: Maths
Date: 07/09/2016

Sigma Notation

Chapter 2 is about understanding the sigma notation for sums. The general delimited form looks as follows:


However, the book declares that the delimited form is not as explicit when compared to a more generalized form and so should only be used for a finished equation rather than when working with them. A more generalized form can look like the following:

$$\sum\limits_{1 \le k \le n}{a_k}$$

When comparing the two there is a big difference in the notation. This is where you can clearly see that the generalized form can be more useful. Below is the same summation (squaring the odd numbers from 1 to 100) but using the two different notations.

$$ \begin{aligned} \text{Delimited} &= \sum_{k = 0}^{49}(2k + 1)^2 \\ \text{Generalized} &= \sum_{ \scriptstyle1 \le k \le 100 \atop\scriptstyle k \text{ odd}} k^2 \end{aligned} $$

Iverson Bracket

This chapter also mentions the Iverson bracket, which is essentially a mathematical if statement where it returns 1 for true and 0 for false. The iverson bracket looks as follows:

$$ [p \text{ prime}] = \begin{cases} 0, & \text{if}\ p \text{ is a prime number} \\ 1, & \text{if}\ p \text{ is not a prime number} \end{cases} $$

With this we can write more expressive delimited sums. Take this for example, a sum for finding the reciprocal of primes smaller than or equal to N.

$$ \sum \limits_{p}[p \text{ prime}][p \le N]/p $$

Tower of Hanoi

With these new tools, the book takes a brief look back at the tower of Hanoi example. As this is the original tower of Hanoi puzzle the solution is different from the previous blog in this series, so ill recap. The recursion is as follows:

$$ \begin{aligned} T_o &= 0;\\ T_n &= 2T_{n-1}, \quad \text{for }n>0 \end{aligned} $$

Now the book mysteriously tell us to divide both sides by $2^n$ (actually multiplying by $2^{-n}$), this is so we can remove the $2$ which is multiplying the $T_{n-1}$. This puts out equation in a form which can easily be converted to a summation. I'll talk a little more about how we find out what number to multiply by in the next blog. Anyway lets do the division.

$$ \begin{aligned} T_0/2^0 &= 0;\\ T_n/2^n &= 2T_{n-1}/2^n + 1/2^n, \quad \text{for }n>0 \\ \textrm{Remember that}\ 2^n/2 &= 2^{n-1} \\ \therefore\ T_n/2^n &= T_{n-1}/2^{n-1} + 1/2^n, \quad \text{for }n>0 \end{aligned} $$

Now we can simplify the equation a little bit by setting $S_n=T_n/2^n$, which gives us the following recursive solutions.

$$ \begin{aligned} S_0 &= 0;\\ S_n &= S_{n-1} + 1/2^n, & \text{for }n>0 \\ &= S_{n-1} + 2^{-n}, & \text{for }n>0 \end{aligned} $$

This gives us the summation below.

$$ S_n = \sum \limits_{k = 1}^{n} 2^{-k} $$

This doesn't exactly look like $2^n-1$. Fortunately the book gives us a small amount of closure by telling us what this series is equal to.

$$\sum \limits_{k=1}^{n}2^{-k} = (\frac{1}{2})^{-1} + (\frac{1}{2})^{-2} + \ldots + (\frac{1}{2})^{-n} = 1 - (\frac{1}{2})^n.$$

Now we can show that the previous sum is equal to $2^n - 1.$

$$ \begin{aligned} S_n &= \sum \limits_{k = 1}^{n} 2^{-k} \\ &= 1 - (\frac{1}{2})^n \\ \end{aligned} $$

Now solve for $T_n$:

$$ \begin{aligned} T_n &= 2^nS_n \\ &= 2^n (1-(\frac{1}{2})^n) \\ &= 2^n (1-2^{-n}) \\ &= 2^n - 2^n \cdot 2^{-n} \\ &= 2^n - 1 \end{aligned} $$

There were a few leaps of faith during this blog post, where we just followed what the book told us to do. These parts will be addressed in later blogs for the same chapter.

Concrete Mathematics: Chapter 1

Topic: Maths
Date: 04/09/2016

This series of blog posts is dedicated to the Concrete Mathematics book by Graham, Knuth and Patashnik. I should note that I'm a beginner when it comes to pretty much all of the contents of this book. I have an A level in maths and a passion for the subject, but I'm by no means good. The book is complex, so these blogs will document my struggles working through it.

Fortunately I didn't struggle too much with the first chapter. Therefore these notes are to just jog my memory of what I learnt.

Recurrent Problems: Tower of hanoi

The tower of hanoi puzzle is the first example that the book works through. So for a little bit of extra information, I'll work through the same problem but the pieces have to move through peg B on their way to peg C. This is the same puzzle as is posed by warmup exercise 2.

The original challenge is to move all the pieces from one peg to another to produce the same order but on a different peg (peg C in this case). You can't place a larger piece on top of a smaller piece, and the idea is to solve it in as fewer moves as possible. However in our version we must make sure to move each piece via peg B.

The starting position for the board is as follows:

tower of hanoi example at the start

Because of our new rules, the only legal move we can make right now is to move the top piece of the tower to peg B. Next we have to move the same piece to peg C, in order to get the next piece of the tower on peg B. If we continue on we should eventually reach the end position, which looks as follows:

tower of hanoi example completed

If we have counted correctly and played optimally, this should have only taken 26 moves for 3 game pieces. Cool. While this is great and all, the purpose of the book is to teach us how to solve recurrent problems. So the real question is, how do we find out the minimum number of turns for $n$ pieces. It's best to start by defining some terms.

$$\begin{aligned} n &= \textrm{Number of pieces} \\ T_n &= \textrm{Number of moves to solve the puzzle} \end{aligned}$$

Now we have some terms we can look at the solutions for a different $n$. The very simplest is $n=0$, here we have no pieces on peg A to move. This means that the puzzle is essentially already solved, and it took 0 turns. Next we can look at $n=1$. With this we are essentially moving along 1 peg each turn until we end up on peg C. So the first move is to peg B, then finally to C. So when $n=1$, $T_n$ is equal to 2. Below is a table showing a few different values for $n$ as well as the corresponding $T_n$.

Number of pieces, $n$ Number of moves, $T_n$
0 0
1 2
2 8
3 26

Eagle eyed people will notice that this is the same as $3^n-1$, crafty people will put the sequence into the on-line encyclopedia of integer sequences (OEIS) and discover the same.

Sadly just noticing the pattern isn't exactly proof that the sequence will continue with our expectations. Instead we need to prove it. This chapter teaches us how to do just that via an induction proof.

An induction proof is where we prove a basis case (the smallest value) such as $n_0=0$, then we prove for values of $n$ which are greater than $n_0$, making sure that the values between $n$ and $n_0$ have also been proved. This is the inductive step. So if $P(0)$ is true and it implies $P(1)$, likewise that $P(1)$ implies $P(2)$. We can assume that $P(n)$ implies $P(n+1)$ and the proof holds true. This is what we'll use later on.

Noticing the recurrence

As this chapter is about recurrent problems, we should find the recurrence in the puzzle. Let's think about the problem when $n$ is 1 and 2. More specifically the end state when $n=1$, and whether or not this state occurs within the puzzle when $n=2$. If you look hard enough you will find that this occurs on the 2nd and 8th move in the $n_2$ puzzle.

recursion in hanoi solution

Great! But how what about in between? Well, it helps to notice that it wouldn't make a difference if the puzzle was mirrored. This means that it will still take $T_{n-1}$ turns even if the pieces are going in the other direction. An example of this in action is for turns 2 & 3 when $n=2$.

recursion in hanoi solution 2

On turn 2 we can see that we have already used $T_{n-1}$ turns, and the next step moves the largest piece to peg B. This adds 1 turn. Now we can think about $T_{n-1}$ but mirrored! So we complete yet another journey back to peg A with $n-1$ pieces, so we can move our largest piece to the finish line, and then the final $T_{n-1}$ to complete the tower.

All in all this turns out to be:

$$\begin{aligned} T_0 &= 0 \\ T_n &= T_{n-1} + 1 + T_{n-1} + 1 + T_{n-1}, \quad \text{for } n > 0. \\ &= 3T_{n-1} + 2 \end{aligned}$$

We can now cancel this down by simply adding 1 to both sides of the equation:

$$\begin{aligned} T_0 + 1 &= 1 \\ T_n + 1 &= 3T_{n-1} + 3, \quad \text{for } n > 0. \\ T_n + 1 &= 3(T_{n-1} + 1) \\ \end{aligned}$$

Let $U_n = T_n + 1:$

$$\begin{aligned} U_0 &= 1\\ U_n &= 3U_{n-1}, \quad \text{for } n > 0. \end{aligned}$$

Now we can see the recurrence as $U_n = 3^n$ (this is more explicit when shown like so, $U_n = 3 \cdot 3^{n-1}$), so it follows that $T_n = 3^n - 1.$ Now all that is left to do is to prove the inductive step. This is done by 'plugging in' the non-recursive into the recursive solution. Thus proving that $n-1$ implies $n$.

$$\begin{aligned} T_n &= 3T_{n-1} + 2 \\ &= 3(3^{n-1} - 1) + 2 \\ &= 3^n - 3 + 2 \\ &= 3^n - 1 \end{aligned}$$

et voila, or should I say Q.E.D.