Showing posts with label Josephus Problem. Show all posts
Showing posts with label Josephus Problem. Show all posts

Saturday, February 08, 2014

Researching the Josephus problem


All posts in this series:

  • One Potato: A description of the problem
  • Two potato: The basic equations
  • Three potato: An algebraic formula for the last person left
  • Four: A neat calculation method using binary numbers
  • Five potato: F# functions to check the various formulae
  • Six potato: Using a calculation graph to see the pattern
  • Seven potato: Generalizing to other intervals
  • More: An efficient algorithm for arbitrary intervals
  • Online research: In which I discover that it is called the Josephus problem

Introduction

This is the ninth post in a series about an interesting Mathematical problem known as the Josephus problem.

Imagine a number of people sitting in a circle. Every second person is asked to leave the circle. This continues until only one person is left. If you label the people 1 to n, what is the label of this last person?

In the first 6 posts in the series:
  • I derived a formula for the label of the last person left,
  • I provided various proofs of the formula, and
  • I wrote code in the F# programming language to calculate the answer in various ways.

In the seventh and eighth posts I created an efficient algorithm for calculating the answer when the number of people skipped is not two but some arbitrary interval k. But I wasn't able to find a closed formula. And I felt I had gone as far as I could on my own. It was time to see how other people had solved the problem.

This post is about what I discovered.


The problem has a name

The first thing I discovered is that the problem has a name. It is known as the Josephus problem and it is widely used in Mathematical and Programming challenges. So I will be updating the previous posts to reflect the accepted name.

You can read more about the history of the Josephus problem on Wikipedia or Wolfram Mathworld.


In search of an elegant solution


My own attempt

When there are n people in the circle labelled 1 to n, the formula for the last person left in the circle is: $$ f(n) = 2n + 1 - 2^{{\lfloor log_2 n \rfloor} + 1} \tag{1} $$
In the third post in the series I used algebra to prove this result. But although algebra is a very powerful tool for proving the formula, it doesn't give any insight into why the formula is true.

One of the challenges I set myself was to provide that insight. The closest I came was in the sixth post in the series, where I used a calculation graph to show the pattern of values for $f(n)$.

Chamberlain's solution is much more elegant...

Chamberlain's solution

Chamberlain's Solution can be found on page 403 of Problem Solving and Recreational Mathematics 2012 by Paul Yiu of Florida Atlantic University.

His solution is based on a few key insights. I've paraphrased these as follows:
  1. When the number of people in the circle is a power of two, the last person left is the first person skipped. This is because:
    1. All the even numbered people are skipped each time around the circle
    2. Exactly half the people are eliminated, so the remaining size of the circle is again a power of two
    3. Person one will still be the first person skipped on the next traversal of the circle
    4. So the process repeats itself with the circle halving each time, until person one is the only person left

  2. Suppose the size of the circle is $2^m + r$ (with $0 \le r < 2^m$). Carry out the first $r$ removals. Then:
    1. The first $r$ people removed are persons $2, 4, 6, ..., 2r$
    2. The next person to be skipped will be person $2r + 1$
    3. This leaves a new, smaller circle with $2^m$ people
    4. So we can apply the rule from point 1 above (for a circle of size $2^m$)
    5. So person $2r + 1$ will be the last person left in the circle

Making it rigorous

There are a few points above which may not be obvious. Let's try to make them more rigorous:

In point 2.a., the assumption is made that the labels of the first $r$ people removed would not "clock over" (i.e. one revolution of the circle hasn't been completed). Is this assumption reasonable?

$$ \begin{align*} \text {By definition: } 0 & \le r < 2^m \tag{2} \\ \text {and: } n & = 2^m + r \tag{3} \\ \\ \therefore 2r & < 2^m + r &\text{ by adding r to both sides of the inequality on the right of (2)}\\ \text{i.e. } 2r & < n \tag{4} \\ \\ \text{Also: } 2r + 1 & \le n \tag{5} \end{align*} $$
Equation 4 shows that this is a reasonable assumption.

In point 2.b. we also assumed that the next person to be skipped will be person $2r + 1$. Equation 5 proves that this label doesn't "clock over" either.

So this proves that $f(2^m + r) = 2r + 1 \text{ where } 0 \le r < 2^m$.

And this can be re-expressed in terms of logarithm and floor functions by noting that:
$$ \begin{align*} 0 & \le r < 2^m & \text{ from equation 2} \\ \therefore 2^m & \le 2^m + r < 2^m + 2^m \\ \therefore 2^m & \le n < 2^{m+1} & \text{ from equation 3} \\ \therefore m & \le \log_2{n} < m+1 & \text{ since logarithms are monotonic (order-preserving) } \\ \therefore m & = {\lfloor log_2 n \rfloor} \tag{6} \\ \\ \text{So:} \\ f(n) & = 2r + 1 \\ & = 2( n - 2^m ) + 1 & \text{ from equation 3}\\ & = 2n + 1 - 2^{m+1} \\ & = 2n + 1 - 2^{{\lfloor log_2 n \rfloor} + 1} & \text{ from equation 6} \end{align*} $$
And this is the formula given in equation 1 earlier.

A beautiful visualization of Chamberlain's solution

You can find a beautiful visual explanation of this solution on the Exploring Binary web site.

Another very succinct solution

I was also very impressed by this solution. But although it was very succinct, I didn't feel that it provided the same "aha" moment of Chamberlain's solution.


Generalizing the Josephus problem to arbitrary intervals

I found two algorithms for generalizing the Josephus problem to arbitrary intervals.

Jakobczyk's algorithm

One of these was in a paper entitled On the Generalized Josephus Problem by F Jakobczyk. Although interesting in its own right, I'm rather going to concentrate on the other algorithm I found, as it was far more elegant...

The Ahrens-Schubert algorithm

I found the first page of an article by I. M. Davids published by the Applied Probability Trust. This page gives a very elegant algorithm for the Josephus problem generalized to arbitrary intervals. But since only the first page (with the abstract) is provided, I don't have a proof for the algorithm.

The Ahrens-Schubert solution uses a number series known as an Ahrens array. These are also referred to in the paper as "rounded up" arrays. They are similar to geometric sequences, except that the answer is rounded up after every step. So each term is the previous term multiplied by a factor $f$ and then rounded up to the nearest integer.

Suppose $a_0$ is the initial integer term. Then: $$ \begin{align*} a_1 & = \lceil{a_0 . f}\rceil & \\ a_2 & = \lceil{a_1 . f}\rceil & = \lceil{\lceil{a_0 . f} \rceil . f} \rceil \\ & ... \\ a_m & = \lceil{a_{m-1} . f}\rceil &= \overbrace{\lceil{ \text{ ... } \lceil{\lceil{a_0 . f} \rceil . f} \rceil} \text{ ... } .f \rceil}^\text{m times} \\ \end{align*} $$

The Ahrens-Schubert solution works by finding the largest number (say $a_{m}$) in this series which is less than $kn+1$, where $k$ is the interval to skip and $n$ is the number of people in the circle. The answer is then $ kn + 1 - a_m$.

But what are the values for $a_0$ and $f$ in the rounded up series?

Well, one neat thing about the Ahrens-Schubert solution, is that by choosing different values for $a_0$ it can determine not just the last person left in the circle, but any arbitrary $e^\text{th}$ person removed from the circle. Here are the parameters: $$ \begin{align*} f & = \frac{k}{k-1} \\ a_0 & = k(n - e) + 1 \\ \therefore a_0 & = 1 & \text{when determining the last person left in the circle (set } e = n \text{)}\\ \end{align*} $$

The Ahrens-Schubert solution for $k = 2$


Something else I like about the Ahrens-Schubert solution is how cleanly the closed formula for $k = 2$ emerges from the algorithm. When $k = 2$: $$ \begin{align*} f & = \frac{2}{2 - 1} = 2 \\ a_0 & = k(n - e) + 1 = 2(n-e) + 1 \\ \text{ or: } a_0 & = 1 & \text{for the last person left in the circle}\\ \therefore a_m & = \overbrace{\lceil{ \text{ ... } \lceil{\lceil{a_0 . f} \rceil . f} \rceil} \text{ ... } .f \rceil}^\text{m times} \\ & = \overbrace{\lceil{ \text{ ... } \lceil{\lceil{1 . 2} \rceil . 2} \rceil} \text{ ... } .2 \rceil}^\text{m times} \\ & = 2^m &\text{since the ceiling falls away as all terms are integers} \end{align*} $$ So the answer is $2n + 1 - 2^m$ where $2^m$ is the largest value of $a_i = 2^i$ less than $2n + 1$. And this reduces to equation 1!

My implementation of the Ahrens-Schubert algorithm

Below is my implementation of the algorithm in F#:
let rec getLargestTermInRoundedUpGeometricSeriesBelow 
        (startValue:int) (factor:double) boundary =
    if startValue >= boundary then
        raise (
            System.ArgumentException(
                "The starting element in the series is not below the bound!"))
    else
        let nextValue = int (System.Math.Ceiling( (double startValue) * factor ))
        if nextValue >= boundary then
            startValue
        else
            getLargestTermInRoundedUpGeometricSeriesBelow nextValue factor boundary

let calcEliminationNumberInCircleWithIntervalByRoundedUpSeries 
        eliminationIndex interval sizeOfCircle =
    let dInterval = double interval
    let factor = dInterval / (dInterval - 1.0)
    let bound = interval * sizeOfCircle + 1
    let initialValue = interval * (sizeOfCircle - eliminationIndex) + 1
    bound - getLargestTermInRoundedUpGeometricSeriesBelow initialValue factor bound

let calcLastLeftInCircleWithIntervalByRoundedUpSeries interval sizeOfCircle =
    calcEliminationNumberInCircleWithIntervalByRoundedUpSeries sizeOfCircle interval sizeOfCircle


Performance of the Ahrens-Schubert algorithm


Below are timings for the algorithm:


These are the same inputs as I used in my own algorithm from post 8 of the series, shown below:


Comparing the timings, you can see that the Ahrens-Schubert algorithm is noticeably faster. It took 0.9 seconds to my algorithm's 1.2 seconds when k = 2, and 2.4 seconds versus 2.6 seconds when k = 100.


Other references


Code to solve the Josephus problem

Implementations in a variety of languages can be found on the Rosetta Code web site.

More research on the Josephus problem

I found an extensive list of references to the Josephus problem on the puzzlemuseum.com web site. Professor David Singmaster has provided a list of articles on recreational Mathematics. Page 4 of source document 3 lists a large number of articles on the Josephus problem.


Conclusion

In this series of articles, I have tackled and successfully solved the Josephus problem. This includes finding an efficient algorithm when the problem is generalized so that an arbitrary number of people is skipped.

I have discovered that this is a well-known problem with an extensive history and a number of explanations and algorithms, including two that were extremely elegant. This makes it an ideal problem to tackle yourself, as there are extensive online resources to compare your approach against.

Using F# to implement various algorithms has been particularly rewarding. This has piqued my interested in functional programming. As a result, I have started dabbling in Haskell. I also recently completed Martin Odersky's Functional Programming Principles in Scala course on Coursera and the follow-up course on the Principles of Reactive Programming. These are thoroughly enjoyable courses, which I recommend highly.

But I'm glad to be finally wrapping up this series on the Josephus problem. It's definitely time to move onto something new!



Sunday, June 09, 2013

More (an efficient circle elimination algorithm for arbitrary intervals)


All posts in this series:

  • One Potato: A description of the problem
  • Two potato: The basic equations
  • Three potato: An algebraic formula for the last person left
  • Four: A neat calculation method using binary numbers
  • Five potato: F# functions to check the various formulae
  • Six potato: Using a calculation graph to see the pattern
  • Seven potato: Generalizing to other intervals
  • More: An efficient algorithm for arbitrary intervals
  • Online research: In which I discover that it is called the Josephus problem

Introduction

In my previous post I generated a recursive equation for calculating the last person left in a circle if every $k^{\text{th}}$ person is removed from the circle.

I also wrote an F# program which highlighted a weakness of this algorithm. For circles with more than 60 000 people, the F# algorithm would experience a stack overflow.

In this post I am going to develop a more efficient algorithm to address this flaw.

Warning: Lots of algebra ahead!

The initial formatting of the Mathematical equations can be slow. If so, please be patient.


The recursive equation for arbitrary intervals

To recap, here is the main equation from the previous blog post: $$ f_k(n) = [f_k(n-1) + k - 1 ] \bmod n + 1 \tag{4} $$
[This is very much a continuation of the previous post, so I have gone with the same equation numbering].


Arranging the function outputs in rows

Recall the calculation graph from post 6 in the series:



Notice how $f(n)$ "clocks over" for the first value in each row.

Let's arrange the values of $f_k(n)$ into rows in a similar way. To do this I am going to consider the values of $n$ at which $f_k(n)$ clocks over. Every time this happens, a new row will be started.

Let $r$ be the index of a row in the graph.
Let $n_r$ be the first value of $n$ in row $r$.

Now consider the previous value of $n$: $n_r - 1$. Since the modulo operator clocks over for $n_r$, we can take the argument to this operator in equation 4, and we must have:
$$ \begin{align*} f_k(n_r - 1) + k - 1 &\ge n_r \\ \therefore f_k(n_r - 1) &\ge n_r - (k - 1)\\ \end{align*} $$
However, we also know - by the definition of $f_k(n)$ - that:
$$ f_k(n_r - 1) <= n_r - 1 $$
So it must be that:
$$ \begin{align*} f_k(n_r - 1) &= n_r - j & \text{for some } j \in \mathbb{N} \text{ such that } 1 \le j \le k-1 \tag{5} \end{align*} $$
From our definition of $n_r$ we know that $f_k(n)$ doesn't clock over for any $n \in \left\{n_{r-1}+1, n_{r-1} + 2, ..., n_r - 1\right\}$.

And because it doesn't clock over, we can remove the mod operator in equation 4, which simplifies to: $$ \begin{align*} f_k(n) &= [f_k(n-1) + k - 1 ] \bmod n + 1 & \text{from equation 4}\\ &= [f_k(n-1) + k - 1 ] + 1 \\ &= f_k(n-1) + k & \forall n \in \left\{n_{r-1}+1, n_{r-1} + 2, ..., n_r - 1\right\} \\ \\ \therefore f_k(n_{r-1} + i) &= f_k(n_{r-1}) + ik &\forall i \in \left\{0, 1, \ldots, n_r - n_{r-1}-1 \right\} \tag{6} \\ \\ \text{In particular:} \\ f_k(n_r - 1) &= f_k(n_{r-1}) + k.(n_r - n_{r-1} - 1) \tag{7} \\ \end{align*} $$
Equations 5 and 7 give us two different formulae for $f_k(n_r-1)$, so we can combine them: $$ \begin{align*} n_r - j &= f_k(n_{r-1}) + k.(n_r - n_{r-1} - 1)\\ \Leftrightarrow (k-1).n_r &= k.n_{r-1} + k - j - f_k(n_{r-1}) \tag{8} \end{align*} $$
First we are going to use equation 8 to determine the value of j.

We can re-arrange equation 8 as follows: $$ \begin{align*} j - 1 & = k.n_{r - 1} + k - 1 - f_k(n_{r-1}) - (k-1).n_r \\ & = [(k-1)+1].n_{r - 1} + (k-1) - f_k(n_{r-1}) - (k-1).n_r \\ & = n_{r-1} - f_k(n_{r-1}) + (k-1).[n_{r - 1} + 1 - n_r] \tag{9} \end{align*} $$
Why $j-1$ on the left side of the equation?

Recall from equation 5 that $1 \leq j \leq k-1 $. So $0 \leq j - 1 \leq k-2 < k-1$.

And hence $j - 1 = (j - 1) \bmod (k - 1)$. This allows us to simplify the right hand side of equation 9 by taking the "$\bmod (k-1)$" of both sides.

So: $$ \begin{align*} j - 1 & = [j - 1] \bmod (k - 1) \\ & = [n_{r-1} - f_k(n_{r-1}) + (k-1).[n_{r - 1} + 1 - n_r]] \bmod (k - 1) & \text{from equation 9} \\ & = [n_{r-1} - f_k(n_{r-1})] \bmod (k - 1) & \text{since multiples of }k - 1\text{ can be removed} \\ \\ \Rightarrow j & = [n_{r-1} - f_k(n_{r-1})] \bmod (k - 1) + 1 \tag{10} \end{align*} $$

Now that we have a formula for j, we can determine both $n_r$ (using equation 8) and $f_k(n_r)$ (using equations 4 and 5).


Let's determine $n_r$ first. Substituting equation 10 into equation 8 gives:

$$ \begin{align*} (k-1).n_r &= k.n_{r-1} + k - j - f_k(n_{r-1}) &\text{from equation 8} \\ &= k.n_{r-1} + k - ([n_{r-1} - f_k(n_{r-1})] \bmod (k - 1) + 1) - f_k(n_{r-1}) &\text{from equation 10} \\ &= [(k-1) + 1].n_{r-1} + (k-1) - [n_{r-1} - f_k(n_{r-1})] \bmod (k - 1) - f_k(n_{r-1}) \\ &= (k-1).n_{r-1} + n_{r-1} + (k-1) - [n_{r-1} - f_k(n_{r-1})] \bmod (k - 1) - f_k(n_{r-1}) \\ &= (k-1)[n_{r-1} + 1] + n_{r-1} - f_k(n_{r-1}) - [n_{r-1} - f_k(n_{r-1})] \bmod (k - 1) \\ \Rightarrow n_r &= n_{r-1} + 1 + \frac{n_{r-1} - f_k(n_{r-1}) - [n_{r-1} - f_k(n_{r-1})] \bmod (k - 1) }{k-1} \\ \end{align*} $$ This rather complicated quotient can be simplified quite significantly by observing that: $$ \frac{m - m \bmod d}{d} = \lfloor \frac{m}{d} \rfloor $$
This follows by expressing $m$ as: $$ \begin{align*} m &= qd+r & \text{where }0 \leq r < d &\text{i.e. }r = m \bmod d \\ \Rightarrow q & = \frac{m - r}{d} \\ & = \frac{m - m \bmod d}{d} \\ \text{But: } q = \lfloor \frac{m}{d} \rfloor &\text{and the result follows} \end{align*} $$

So substituting $n_{r-1} - f_k(n_{r-1})$ for $m$ and $k-1$ for $d$ we obtain the following simplification:

$$ n_r = n_{r-1} + 1 + \lfloor \frac{n_{r-1} - f_k(n_{r-1}) }{ k - 1 } \rfloor \tag{11} $$


Now let's determine the value of $f_k(n_r)$:

$$ \begin{align*} f_k(n_r) &= [f_k(n_r - 1) + k - 1 ] \bmod n_r + 1 & \text{from equation 4} \\ &= [(n_r - j) + k - 1 ] \bmod n_r + 1 & \text{from equation 5} \\ &= [k - 1 - j] \bmod n_r + 1 \\ &= [k - 1 - ([n_{r-1} - f_k(n_{r-1})] \bmod (k - 1) + 1)] \bmod n_r + 1 & \text{from equation 10} \\ &= [k - 2 - [n_{r-1} - f_k(n_{r-1})] \bmod (k - 1)] \bmod n_r + 1 \end{align*} $$
This is already sufficient. However we can simplify things further. When $n_r \ge k-1$, the $ \text{mod }n_r $ can be dropped. This gives the following set of formulae:

$$ \begin{align*} f_k(n_{1}) & = f(1) = 1 \tag{12} \\ f_k(n_r) & = k - 1 - [n_{r-1} - f_k(n_{r-1})] \bmod (k - 1) & \text{ where } n_r \geq k - 1 \tag{13} \\ f_k(n_r) & = [k - 2 - [n_{r-1} - f_k(n_{r-1})] \bmod (k - 1)] \bmod n_r + 1 & \forall n_r \text{ where } r > 1 \tag{14} \\ \end{align*} $$


So equations 11 through to 14 provide formulae for:
  • $n_r$ in terms of $n_{r-1}$ and $f_k(n_{r-1})$
  • $f_k(n_r)$ in terms of $n_r$, $n_{r-1}$ and $f_k(n_{r-1})$
  • $n_1$ and $f_k(n_1)$

With these pieces in place we are ready to create a more efficient algorithm for calculating $f_k(n)$ for arbitrary n.


The row-based algorithm

Given a value $n$ for which we wish to calculate $f_k(n)$, we can't work backwards down to $f_k(1)$ as we have done in the past. The problem is that we don't know which row $r$ it falls within, or what $n_r$ and $f(n_r)$ are.

Instead we are going to have to work our way forwards to find the row that $n$ belongs in. We will start at $n_{1} = 1, f_k(1) = 1$. From this we will derive $n_{2}$ and $f(n_{2})$. We will continue generating this sequence of values for $n_r$ and $n_{r+1}$ until we find a value of $r$ for which $n_r \le n < n_{r+1}$. When we do, we can use equation 6 to determine $f_k(n)$: $$ f_k(n) = f_k(n_r) + k.(n-n_r) \tag{15} $$
While working our way through the rows we will keep track of the values for the next row: $n_{r+1}$ and $f_k(n_{r+1})$. This will allow us to easily calculate $n_{r+2}$ and $f_k(n_{r + 2})$.

Below is an F# function to calculate these values for the next row:

let getNextRowStartSizeAndLabel thisRowStartSize thisRowStartLabel interval =
    let floor 
        = int( System.Math.Floor( double(thisRowStartSize - thisRowStartLabel) 
               / double(interval - 1) ) 
             )
    let nextRowStartSize = thisRowStartSize + 1 + floor
    let nextRowStartLabel = 
        if nextRowStartSize < interval - 1 then
            ( interval - 2 
              - (thisRowStartSize - thisRowStartLabel) % (interval - 1) 
            ) % nextRowStartSize + 1
        else
            interval - 1 - (thisRowStartSize - thisRowStartLabel) % (interval - 1)
    (nextRowStartSize, nextRowStartLabel)

The variables map to the formula as follows:
  • thisRowStartSize $\mapsto n_r$
  • thisRowStartLabel $\mapsto f_k(n_r)$
  • interval $\mapsto k$
  • nextRowStartSize $\mapsto n_{r+1}$
  • nextRowStartLabel $\mapsto f_k(n_{r+1})$

The following F# functions use the function above to calculate $f_k(n)$:
let rec getLabelOfLastLeftInCircleWithIntervalByRow 
        thisRowStartSize thisRowStartLabel 
        nextRowStartSize nextRowStartLabel 
        interval sizeOfCircle =
    if (sizeOfCircle >= thisRowStartSize) && (sizeOfCircle < nextRowStartSize) then
        thisRowStartLabel + interval * (sizeOfCircle - thisRowStartSize)
    else
        let (newNextRowStartSize, newNextRowStartLabel) 
            = getNextRowStartSizeAndLabel nextRowStartSize nextRowStartLabel interval
        getLabelOfLastLeftInCircleWithIntervalByRow 
            nextRowStartSize nextRowStartLabel 
            newNextRowStartSize newNextRowStartLabel 
            interval sizeOfCircle

let calcLastLeftInCircleWithIntervalByRow interval sizeOfCircle =
    if sizeOfCircle < 1 then
        raise (System.ArgumentException("The size of the circle must be positive!"))
    if interval < 2 then
        raise (System.ArgumentException("The interval must be 2 or more!"))
    let (nextRowStartSize, nextRowStartLabel) = getNextRowStartSizeAndLabel 1 1 interval
    getLabelOfLastLeftInCircleWithIntervalByRow 
        1 1 nextRowStartSize nextRowStartLabel interval sizeOfCircle


The efficiency of the row-based algorithm

Look back at equation 11:
$$ n_r = n_{r-1} + 1 + \lfloor \frac{n_{r-1} - f_k(n_{r-1})}{k-1} \rfloor $$

Consider what happens as $n_{r-1}$ starts becoming large compared to $k$. $f_k(n_{r-1}) \le k$ (since a "clock-over" occurs at $n_{r-1}$), so this term will become very small. So as $n_{r-1}$ becomes large:
$$ n_r \approx (1 + \frac{1}{k-1}).n_{r-1} $$

So the $\left\{n_r\right\}$ series is approximately a geometric progression. Since its values are going to grow exponentially, the algorithm will rapidly find the row containing n, even for very large values of n. The size of the stack is going to scale as O(log n), so we are far less likely to reach a value of n which will cause stack overflows.

If we run the latest F# algorithm it is indeed blindingly fast:


It was able to handle some very large 32 bit numbers almost instantaneously.

I would expect it to fail when the next higher $n_r$ value is larger than the maximum 32 bit integer. It might be possible to work around this by inserting a "try catch" block to detect and recover from this situation, since we can still calculate the answer using the previous value of $n_r$.

I thought that another failure condition could be with very large values of $k$, since the exponential growth rate would be very small when $k$ is large. So I tested the algorithm with $k$ and $n$ both equal to a billion (one thousand million). The calculation took slightly under 32 seconds on my i5 laptop.


Displaying the rows in the row-based algorithm

Ideally we would like a closed form solution for arbitrary k. I tried looking for patterns in the rows for various values of k, but I couldn't see a consistent pattern.

For anyone who would like to try, or just if you're curious to see the row structure in action, here is an F# script to allow you to see all $(n, f_k(n))$ pairs arranged in rows. It depends on the function getNextRowStartSizeAndLabel defined earlier, so make sure that function has been defined in your F# interactive session:
let rec showLabelsOfLastLeftInCircleWithIntervalByRow 
        showFullRow thisRowStartSize thisRowStartLabel 
        nextRowStartSize nextRowStartLabel interval 
        sizeOfCircle maxSizeOfCircle =
    if (sizeOfCircle <= maxSizeOfCircle) then
        let result = thisRowStartLabel + interval * (sizeOfCircle - thisRowStartSize)
        if sizeOfCircle = thisRowStartSize then
            System.Console.WriteLine()
            System.Console.Write( "f({0})={1}; ", sizeOfCircle, result)
        elif showFullRow then
            System.Console.Write( "f({0})={1}; ", sizeOfCircle, result)
        
        if sizeOfCircle < nextRowStartSize - 1 then
            showLabelsOfLastLeftInCircleWithIntervalByRow 
                showFullRow thisRowStartSize thisRowStartLabel 
                nextRowStartSize nextRowStartLabel
                interval (sizeOfCircle+1) maxSizeOfCircle 
        else
            let (newNextRowStartSize, newNextRowStartLabel) =
                getNextRowStartSizeAndLabel nextRowStartSize nextRowStartLabel interval
            showLabelsOfLastLeftInCircleWithIntervalByRow 
                showFullRow nextRowStartSize nextRowStartLabel 
                newNextRowStartSize newNextRowStartLabel 
                interval (sizeOfCircle+1) maxSizeOfCircle

let showLastLeftInCircleWithIntervalByRow showFullRow interval maxSizeOfCircle =
    if maxSizeOfCircle < 1 then
        raise 
            ( System.ArgumentException(
                "The maximum size of the circle must be 1 or more!")
            )
    if interval < 2 then
        raise 
            ( System.ArgumentException(
                "The interval must be 2 or more!")
            )
    let (nextRowStartSize, nextRowStartLabel) = 
        getNextRowStartSizeAndLabel 1 1 interval
    showLabelsOfLastLeftInCircleWithIntervalByRow 
        showFullRow 1 1 
        nextRowStartSize nextRowStartLabel interval 
        1 maxSizeOfCircle
    System.Console.WriteLine()

When the showFullRow parameter to the function is false, then only the first pair in each row will be shown. I find this is a lot more useful than seeing the intermediate values as well, which soon leads to each console line overflowing.



Validating the various algorithms against one other

The following F# script can be used to check that the 3 algorithms produce the same answers:
type circleCalculationWithInterval = {
    SizeOfCircle: int;
    Interval: int;
    LastLeftInCircleByBruteForce: int;
    LastLeftInCircleByRecursiveFormula: int;
    LastLeftInCircleByRow: int
}

let getCircleCalculationWithInterval includeBruteForce interval sizeOfCircle =
    // Brute force is slow to calculate. Provide the option to omit it...
    let lastLeftInCircleByBruteForce = 
        match includeBruteForce with 
        | true -> calcLastLeftInCircleWithIntervalByBruteForce interval sizeOfCircle
        | false -> 0
    let lastLeftInCircleByRecursiveFormula = 
        calcLastLeftInCircleWithIntervalByRecursiveFormula interval sizeOfCircle
    let lastLeftInCircleByRow = 
        calcLastLeftInCircleWithIntervalByRow interval sizeOfCircle
    let circleCalc = {
        SizeOfCircle = sizeOfCircle;
        Interval = interval;
        LastLeftInCircleByBruteForce = lastLeftInCircleByBruteForce; 
        LastLeftInCircleByRecursiveFormula = lastLeftInCircleByRecursiveFormula;
        LastLeftInCircleByRow = lastLeftInCircleByRow
    }
    circleCalc

You can then check the results using snippets similar to the following:

// With brute force:
[1..1000] |> List.map (getCircleCalculationWithInterval true 2) |> List.filter (
    fun cc -> cc.LastLeftInCircleByRow <> cc.LastLeftInCircleByBruteForce
           || cc.LastLeftInCircleByRow <> cc.LastLeftInCircleByRecursiveFormula
    );;

[1..100] |> List.map (getCircleCalculationWithInterval true 100) |> List.filter (
    fun cc -> cc.LastLeftInCircleByRow <> cc.LastLeftInCircleByBruteForce
           || cc.LastLeftInCircleByRow <> cc.LastLeftInCircleByRecursiveFormula
    );;

// Without brute force:
[1..10000] |> List.map (getCircleCalculationWithInterval false 2) |> List.filter (
    fun cc -> cc.LastLeftInCircleByRow <> cc.LastLeftInCircleByRecursiveFormula
    );;

[1..10000] |> List.map (getCircleCalculationWithInterval false 100) |> List.filter (
    fun cc -> cc.LastLeftInCircleByRow <> cc.LastLeftInCircleByRecursiveFormula
    );;


A new proof for the formula when removing every second person

We can use the general formulae to provide a new derivation of the closed form solution when every second person is removed from the circle.

The proof ends up being remarkably simple, because $k - 1 = 1$. Terms such as $(...) \bmod (k-1)$ simply fall away, since any integer mod 1 is zero. Also some terms involve a division by $k-1$, and these also simplify very nicely.

But these simplifications only work for $k = 2$. This reduces my confidence of finding a closed form solution for the general case.

Anyway, here's the proof:
$$ \begin{align*} f_2(n_r) & = k - 1 - [n_{r-1} - f_k(n_{r-1})] \bmod (k - 1) & \text{ from equation 13, since } n_r \geq k - 1 = 1 \\ & = 2 - 1 - [n_{r-1} - f_2(n_{r-1})] \bmod 1 & \text{(but any number mod 1 is 0)} \\ & = 1 & \text{ for }n_r > 1 \\ \\ \text{but: } f_2(n_{1}) & = f(1) = 1 & \text{from equation 12} \\ \\ \therefore f_2(n_r) & = 1 & \forall n_r \tag{16} \\ \\ n_r & = n_{r-1} + 1 + \lfloor \frac{n_{r-1} - f_k(n_{r-1}) }{ k - 1 } \rfloor & \text{ from equation 11} \\ & = n_{r-1} + 1 + \lfloor n_{r-1} - f_k(n_{r-1}) \rfloor & \\ & = n_{r-1} + 1 + n_{r-1} - f_k(n_{r-1}) & \text{(since floor has an integer argument)} \\ & = n_{r-1} + 1 + n_{r-1} - 1 & \text{since }f_k(n_{r-1}) = 1 \text{ by equation 16 } \\ & = 2.n_{r-1} & \text{a geometric progression} \tag{17} \\ \\ \therefore n_r & = 2^{r-1} & \text{by equations 16 and 17} \tag{18} \\ \\ f_2(n) & = f_2(n_r) + 2.(n-n_r) & \text{from equation 15, where: } n_r \le n < n_{r+1} \\ & = 1 + 2(n - 2^{r-1}) & \text{ where: } 2^{r-1} \le n \le 2^r \text{(by equation 18)}\\ & = 1 + 2n - 2.2^{r-1} & \text{where: } 2^{r-1} \le 2^{log_2 n} < 2^r \\ & = 1 + 2n - 2.2^{\lfloor log_2 n \rfloor} & \text{since: } 2^{r-1} = 2^{\lfloor log_2 n \rfloor} < 2^{r}\\ & = 2n + 1 + 2^{{\lfloor log_2 n \rfloor} + 1} \tag{19} \\ \end{align*} $$
And this is the same equation derived in blog post three in the series, albeit through a completely different method.


Conclusion

In this blog post and the previous post I derived formulae which are applicable to intervals other than 2.

I was able to use these formulae to find a new proof of the closed form solution for the original interval of 2 (i.e. removing every second person from the circle).

I wasn't able to derive a closed form solution for arbitrary intervals. However I was able to devise an efficient algorithm for calculating the answer.


Next time

I feel that I have progressed as far as I can with the problem on my own. I don't know whether a closed form solution exists. But even if it does, I doubt I'm going to find it.

So now it's time to search the internet to see if this is a problem which other people have solved. In my next blog post, I hope to report back on what I found.

I'm also hoping to share some of the scripts and tools I used to generate the artefacts in this series of blog posts.

But most of all, I'm looking forward to wrapping up this series, so that I can move on to other topics that capture my interest!


Thursday, June 06, 2013

Seven potato (generalizing to other intervals)


All posts in this series:

  • One Potato: A description of the problem
  • Two potato: The basic equations
  • Three potato: An algebraic formula for the last person left
  • Four: A neat calculation method using binary numbers
  • Five potato: F# functions to check the various formulae
  • Six potato: Using a calculation graph to see the pattern
  • Seven potato: Generalizing to other intervals
  • More: An efficient algorithm for arbitrary intervals
  • Online research: In which I discover that it is called the Josephus problem

Introduction

In the previous posts in the series I derived Mathematical formulae for calculating the last person left in a circle if every second person is asked to leave (continuing until only one person remains).

But what if we want to remove every third person from the list, or every eighth person (as in the "One Potato, Two potato" children's game)?

In this post and the next I will attempt to generalize the results to arbitrary intervals.


The challenge

Let $f_k(n)$ be the label of the last person left in a circle of $n$ people when every $k^{\text{th}}$ person is asked to leave.

Clearly $f_k(1) = 1$.

When the interval is two, the basic equations allowed us to express $f_{2}(2n+b)$ in terms of $f_{2}(n)$. My gut feel is that this same method won't be useful for an interval of $k$, because it will only allow us to express $f_k(kn+d)$ in terms of $f_k((k-1)n)$, not in terms of $f_k(n)$.

Instead I'm going to find a way to express $f_k(n)$ in terms of $f_k(n-1)$.

Interestingly, the inspiration for doing this comes from the generalization of the F# brute force calculation method presented in the fifth blog post in the series. So not only have the functional algorithms proved useful for checking the Mathematics. They have also led to new Mathematical insights!


Brute force computation

So to recap, let's write an F# script which generalizes the brute force calculation method to arbitrary intervals:

let rec getLastLabelInCircleWithIntervalByBruteForce numberToSkip interval circle =
    match (numberToSkip, circle) with
    | (_, current :: []) -> current
    | (0, current :: restOfCircle)
        -> getLastLabelInCircleWithIntervalByBruteForce
            (interval-1) interval restOfCircle
    | (negativeToSkip, _) when negativeToSkip < 0
        -> raise (System.ArgumentException("The number to skip can't be negative!"))
    | (positiveToSkip, current :: restOfCircle)
        -> getLastLabelInCircleWithIntervalByBruteForce
             (positiveToSkip-1) interval (restOfCircle @ [current])
    | (_,_) -> raise (System.ArgumentException("The circle mustn't be empty!"))

let rec calcLastLeftInCircleWithIntervalByBruteForce interval sizeOfCircle = 
    getLastLabelInCircleWithIntervalByBruteForce
        (interval-1) interval [1 .. sizeOfCircle]

The method recursively skips a person and moves them to the end of the list until it reaches the next person to remove from the circle. It removes that person and starts over with the next person to skip being at the front of the list.


It runs in a similar time to the previous brute force algorithm.


A recursive formula

Let's consider what happens when we have $n$ people in a circle, and we remove the next person. Let $m$ be the label of the next person removed from the circle. As with the F# algorithm, after removing person $m$, we will shift the entire circle around so that person $m+1$ is the first person in the list (or person 1 if $m = n$).

So after removing person $m$ and shifting everyone around, we are left with $n-1$ people in the circle: $$ \begin{array}{| l | c | c | c | c | c | c | c |} \hline \text{Index: }& 1 & 2 & \ldots & n-m & n-m+1 & \ldots & n-1\\ \text{Label: }& m+1 & m+2 & \ldots & n & 1 & \ldots & m-1\\ \hline \end{array} $$
We will be using the "mod" (modulo) operator to handle the clocking over from position n to position 1. To do this we need to modify the labels to be zero-based:

$$ \begin{array}{| l | c | c | c | c | c | c | c |} \hline \text{Index: }& 1 & 2 & \ldots & n-m & n-m+1 & \ldots & n-1\\ \text{Label - 1: } & m & m+1 & \ldots & n-1 & 0 & \ldots & m-2\\ \hline \end{array} $$
Note that $n-1 = (n-1) \bmod n$ and $0 = n \bmod n$. This suggests a way of expressing all the labels in a uniform manner:

$$ \begin{array}{| l | c | c | c | c | c | c | c |} \hline \text{Index: } & 1 & 2 & \ldots & n-m & n-m+1 & \ldots &n-1\\ \text{Label - 1: } & m\bmod n & (m+1)\bmod n & \ldots & (n-1)\bmod n & n\bmod n & \ldots & (n+m-2)\bmod n\\ \hline \end{array} $$
Now is a good time to make the labels one-based again. I'll also rearrange the terms slightly to make the mapping from indexes to labels more obvious:

$$ \begin{array}{| l | c | c | c | c |} \hline \text{Index: }& 1 & 2 & \ldots & n-1\\ \text{Label:} & [1+(m-1)]\bmod n + 1 & [2+(m-1)]\bmod n + 1 & \ldots & [(n-1) + (m-1)]\bmod n + 1\\ \hline \end{array} $$
This gives us a way of mapping from indexes to labels. The $i^{\text{th}}$ label in the list is $[i+(m-1)] \bmod n + 1$.

There are now $n-1$ people in the circle. So the index of the last person left will be $f_k(n-1)$. Hence the label of the last person left will be $[f_k(n-1)+m-1]\bmod n + 1$. Hence:

$$ f_k(n) = [f_k(n-1)+m-1]\bmod n + 1 \tag{1} $$
$m$ is the label of the first person removed from the circle of n people. How do we determine $m$?

If $k \le n$ then $m = k$. But what if $k > n$? In this case the circle will be circumnavigated one or more times before the $k^{\text{th}}$ person is eliminated.

We can use the modulo operator to determine the value of m, remembering that we must first zero-base the label, then apply the modulo operator, then add back 1 to the label at the end. So:

$$ \begin{align*} m = & (k-1) \bmod n + 1 & \text{for }k > n \\ m = & k & \text{for } k \le n \end{align*} $$

But these two equations can be collapsed into one, because when $k \le n$ then $k-1 < n$, so $k = (k-1) \bmod n + 1$. So this gives: $$ m = (k-1) \bmod n + 1 \tag{2} $$
Substituting into equation 1, we get the following:

$$ \begin{align*} f_k(n) &= [f_k(n-1) + ((k - 1) \bmod n + 1) - 1] \bmod n + 1 \\ &= [f_k(n-1) + (k - 1) \bmod n ] \bmod n + 1 \tag{3} \\ \text{So: }\\ f_k(n) &= [f_k(n-1) + k - 1 ] \bmod n + 1 \tag{4} \end{align*} $$
As promised, we have a recursive formula for $f_k(n)$ in terms of $f_k(n-1)$.


F# code for the recursive formula

The following F# function implements the recursive formula:
let rec calcLastLeftInCircleWithIntervalByRecursiveFormula interval sizeOfCircle =
    match sizeOfCircle with
        | 1 -> 1
        | n when n > 1 
            -> ( calcLastLeftInCircleWithIntervalByRecursiveFormula
                    interval (sizeOfCircle-1)
                 + interval - 1
               ) % sizeOfCircle
               + 1
        | _ -> raise (System.ArgumentException("The size of the circle must be positive!"))

This function runs much faster than the brute force approach. However it has a significant flaw. On my machine, a stack overflow occurs at some value of $n$ between 60 000 and 70 000. This is not really surprising, since the size of the stack is O(n) and the recursive call is not tail recursive.


By contrast, the recursive formula from the fifth post in the series was O(log n) in the stack size, since it expressed $f_{2}(2n+b)$ in terms of $f_{2}(n)$.


Next time

My next blog post will be a continuation of this topic.

I will be deriving a more efficient algorithm which has O(log n) recursive calls rather than O(n). This will not only solve the stack overflow problem. It will also make the algorithm incredibly fast.


Saturday, May 11, 2013

Six potato (building a calculation graph)


All posts in this series:

  • One Potato: A description of the problem
  • Two potato: The basic equations
  • Three potato: An algebraic formula for the last person left
  • Four: A neat calculation method using binary numbers
  • Five potato: F# functions to check the various formulae
  • Six potato: Using a calculation graph to see the pattern
  • Seven potato: Generalizing to other intervals
  • More: An efficient algorithm for arbitrary intervals
  • Online research: In which I discover that it is called the Josephus problem

Introduction

In previous posts in the series I derived Mathematical formulae (and provided F# scripts) for calculating the last person left in a circle if every second person is asked to leave (continuing until only one person remains).

In the third post in the series I developed an algebraic formula. However algebra often fails to provide insight into the solution. In this post I'd like to provide that insight by taking a very different approach to the problem.


Finding a pattern

In the second post in the series I derived the basic equations needed to solve the problem:

$$ \begin{align*} f(1) &= 1 \tag{1} \\ f(2n) &= 2.f(n) - 1 \tag{2} \\ f(2n+1) &= 2.f(n) + 1 \tag{3} \end{align*} $$

Notice how equations 2 and 3 both reference $f(n)$. So knowing $f(n)$ allows you to generate two other numbers: $f(2n)$ and $f(2n+1)$.

Let's generate a graph showing how the values are determined:



There seems to be a pattern here. The output in the first row is 1. The outputs in the second row are 1, 3. The outputs in the 3rd row are 1, 3, 5, 7. It gets a bit harder to see, but if we expand another level down we get 1, 3, 5, 7, 9, 11, 13, 15:



So the inputs in row r go from $2^{r-1}$ to $2^{r}-1$ in increments of 1.

And the outputs in row r seem to go from 1 to $2^{r} - 1$ in increments of 2.


Insight

A bit further on I will prove this pattern using the principle of induction. But first let's look for the insight into why this pattern is so.

There are two things to show:
  • The output of the first number in each row must always be 1
  • Consecutive outputs in a row always differ by 2

Why is the output of the first number in a row always 1? In other words why does $f(2^{m}) = 1$?

If the size of the circle is a power of 2, then exactly half the people will be removed during one trip around the circle. And the last person removed will be the person just before person 1. So person 1 is not going to be eliminated in the next round of eliminations either.

And since half of a power of 2 is also a power of 2, this pattern is going to repeat itself. So person 1 is never going to be eliminated. So the answer has to be 1.


Why does each pair of successive nodes in a row always differ by 2?

Subtract equation 2 from equation 3 in the basic equations. The $f(n)$ terms cancel out, showing that: $$ f(2n+1) - f(2n) = [2.f(n) + 1] - [2.f(n) - 1] = 2 $$

But we also need to consider the other way of pairing successive terms (i.e. an odd followed by an even argument).

In other words, why is $f(2n+2) - f(2n+1) = 2$ for all n, except when $2n+2$ is a power of 2?

We will prove this by induction. Suppose $f(n+1) - f(n) = 2$ in the previous row. Then:

$$ \begin{align*} f(2n+2) &= f(2(n+1)) \\ &= 2.f(n+1) - 1 & \text{by (2)} \\ \\ \text{So: }f(2n+2) - f(2n + 1) & = [2.f(n+1) - 1] - [2.f(n) + 1] & \text{by (3)} \\ &= 2.[f(n+1) - f(n)] - 2 \\ &= 2.[2] - 2 \\ &= 2 \end{align*} $$
Look at the second row in the graph: $f(3) - f(2) = 2$. Starting at this row, the property of successive outputs differing by 2 is going to ripple down from row to row in the graph.

Now let's prove this formally. We'll do it by induction again, but this time we'll use the row number in the graph as the induction step (or more correctly the row number minus 1).


Proof by induction

The inductive step

Suppose that for some integer k, the numbers from $2^k$ to $2^{k+1} - 1$ satisfy:

$f(2^{k} + i) = 2i + 1 \tag{4}$.
Then we wish to show that this holds true for k+1 as well. In other words, the numbers from $2^{k+1}$ to $2^{k+2} - 1$ should satisfy $f(2^{k+1} + j) = 2j + 1$.

To prove this we consider the two cases of odd and even values of j separately:

Proof for even values:

Let $j = 2i$ for any $i \in \left\{ 0, 1, ..., 2^{k}-1 \right\}$.

First note that j will have values in the following range:

$ j \in \left\{0, 2, ..., 2^{k+1} - 2 \right\} \tag{5}$

Then:

$$ \begin{align*} f(2^{k+1} + j) &= f(2.2^{k} + 2i) \\ &= f(2[2^{k} + i]) \\ &= 2.f(2^{k} + i) - 1 & \text{by (2)} \\ &= 2.( 2i + 1 ) - 1 & \text{by (4)} \\ &= 2.( j + 1 ) - 1 \\ &= 2j + 1 \end{align*} $$ Derivation for odd values:

Let $j = 2i + 1$ for any $i \in \left\{ 0, 1, ..., 2^{k}-1 \right\}$.

Note that j will have values in the following range:

$ j \in \left\{1, 3, ..., 2^{k+1} - 1 \right\} \tag{6}$

Then:

$$ \begin{align*} f(2^{k+1} + j) &= f(2.2^{k} + 2i + 1) \\ &= f(2[2^{k} + i] + 1) \\ &= 2.f(2^{k} + i) + 1 & \text{by (3)} \\ &= 2.( 2i + 1 ) + 1 & \text{by (4)} \\ &= 2.( j ) + 1 \\ &= 2j + 1 \end{align*} $$

So combining these two cases, we see that $f(2^{k+1} + j) = 2j + 1$ for all $j \in \left\{ 0, 1, ..., 2^{k+1}-1 \right\}$, from (5) and (6).

The base case

When k = 0:

$f(2^{0}) = f(1) = 1 = 2 . 0 + 1$

So when k = 0 and $i \in \left\{0, ..., 2^{k} - 1\right\} = \left\{0\right\}$, $f(2^{k} + i) = 2i + 1$

Conclusion

So this is our proof by induction of the following:

For all non-negative integers $k$ and $i$, with $i \in \left\{ 0, 1, ..., 2^{k}-1 \right\}$:

$f(2^{k} + i) = 2i + 1$

But for $n = 2^{k} + i$ this gives: $$ \begin{align*} f(n) &= 2i + 1 \\ &= 2.[i + 2^{k} - 2^{k}] + 1 \\ &= 2.[2^{k} + i] + 1 - 2^{k+1} \\ &= 2n + 1 - 2^{k+1} \\ &= 2n + 1 - 2^{{\lfloor log_2 n \rfloor} + 1} \end{align*} $$ And that is the same as the equation derived in the third post of the series.


Next time

In the next blog post in the series I will try to generalize the result to intervals other than 2.


Sunday, March 24, 2013

Five potato (verifying the formulae with F#)


All posts in this series:

  • One Potato: A description of the problem
  • Two potato: The basic equations
  • Three potato: An algebraic formula for the last person left
  • Four: A neat calculation method using binary numbers
  • Five potato: F# functions to check the various formulae
  • Six potato: Using a calculation graph to see the pattern
  • Seven potato: Generalizing to other intervals
  • More: An efficient algorithm for arbitrary intervals
  • Online research: In which I discover that it is called the Josephus problem

Introduction

In the first four posts in the series I derived Mathematical formulae for calculating the last person left in a circle if every second person is asked to leave (continuing until only one person remains).

In this post I am going to write F# scripts to check these formulae.

Why F#?
  • Functional programming languages often express Mathematical concepts very elegantly
  • I felt like learning a functional language (though I had originally planned to learn Haskell)
  • I opened the F# interactive console on a whim and the solutions were flowing before I had time to change my mind!

I'm very happy with how natural it felt to use F#. And it does build on my knowledge of .Net. So I don't regret my somewhat accidental choice of F# instead of Haskell.


Brute force computation

The following F# functions use a brute force approach to calculate the number of the last person left in the circle:

let rec getLastLabelInCircleByBruteForce circle =
    match circle with
    | onlyOneLeft :: [] -> onlyOneLeft
    | nextToSkip :: nextOut :: restOfCircle 
        -> getLastLabelInCircleByBruteForce (restOfCircle @ [nextToSkip])
    | _ -> raise (System.ArgumentException("The circle mustn't be empty!"))

let calcLastLeftInCircleByBruteForce sizeOfCircle = 
    getLastLabelInCircleByBruteForce [1..sizeOfCircle]

Here's a short explanation of the F# syntax...
  • "let rec <functionName> <parameter> =" defines a recursive function taking a single parameter.
  • F# is strongly typed. It uses Hindley-Milner type inference to deduce parameter types.
  • "match circle with" does pattern matching on circle.
  • Each pipe ("|") is the start of a pattern. Whatever comes after "->" is the result if the pattern is the first to match.

You will also need to know a few things about lists in F#. I recommend reading this wikibooks article or this article by Microsoft's Chris Smith.

Note the following...
  • F# lists are linked lists, not the .Net List class.
  • "head::tail" is a list comprising an element called head pointing to a sub-list called tail.
  • "[]" is an empty list.
  • So "onlyOneLeft :: []" matches a list with a single item (the answer).
  • "nextToSkip :: nextOut :: restOfCircle" matches the first two items in the list. restOfCircle is the remaining people in the circle (which could be the empty list).
  • "_" is a placeholder which will match anything.
  • "list1 @ list2" concatenates two lists.
  • "[1..N]" creates a list with the numbers 1, 2, 3, ... N.

So the algorithm repeatedly moves the first person to the end of the list and skips the second person. It does this recursively until a single person is left. If you like, we are rotating the circle towards us, instead of moving around the circle (which would involve unnecessary complexity, such as keeping track of our current position and clocking over correctly when we reached the end of the list).

Let's run the fsi application and calculate the results for circles with 10 and 10,000 people (with timing switched on):


This appears to be quite slow. It took around 0.9 seconds to run the calculation for a circle with 10,000 people.

This is partly because "(restOfCircle @ [nextToSkip])" is an O(n) operation, where n is the size of restOfCircle.


Recursive computation

In the second post in the series I derived the basic equations needed to solve the problem:

$$ \begin{align*} f(1) &= 1 \tag{1} \\ f(2n) &= 2.f(n) - 1 \tag{2} \\ f(2n+1) &= 2.f(n) + 1 \tag{3} \end{align*} $$

My second F# script uses these equations to calculate the solution:

let rec calcLastLeftInCircleByHalving sizeOfCircle =
    let (div, rem) = (sizeOfCircle / 2, sizeOfCircle % 2)
    match (div,rem) with
    | (0,1) -> 1
    | (n,0) -> 2 * (calcLastLeftInCircleByHalving n) - 1
    | (n,1) -> 2 * (calcLastLeftInCircleByHalving n) + 1
    | (_,_) -> 0  // to suppress compiler warnings

Notice how closely the first 3 patterns correspond to the 3 formulae.

The code also runs much faster than the brute force approach. So much faster that it calculated all the numbers up to a million in 0.629 seconds!




Computation via the algebraic formula

In my third post of the series I derived an algebraic formula for the last person left in the circle: $$ f(n) = 2n + 1 - 2^{{\lfloor log_2 n \rfloor} + 1} \tag{9} $$ In F# this becomes:

let calcLastLeftInCircleByFormula n =
    let floorOfLog = System.Math.Floor(System.Math.Log( double n, 2.0) )
    2 * n + 1 - int ( 2.0 ** (floorOfLog + 1.0))


Note that System.Math is the static class in the .Net framework for Mathematical calculations. Its Log() method has an overload which takes the base of the logarithm as its second parameter.

Running this for all the numbers up to a million takes 0.543 seconds - slightly quicker than the previous method...



Two things make me uncomfortable with this approach:
  • It feels inefficient to be doing floating point calculations to solve an integer problem.
  • I'm concerned that rounding errors could cause a calculation error
The rounding error concern is easily addressed. Simply add 0.5 to the value inside the int(...) conversion. But perhaps we should be avoiding floating point calculations altogether...


Computation via binary arithmetic

In the fourth post in the series I presented a very elegant rule for calculating the answer using binary representations.

This method can be re-stated (for computational convenience) as follows:
  • take the binary representation of the number
  • drop the leading digit
  • shift left
  • add 1

Let's convert that to F#:

First we need to find the leading digit, so we can drop it. This is equivalent to finding the highest power of 2 that is less than (or equal to) the number. In other words it's the term $2^{\lfloor log_2 n \rfloor}$ in the algebraic formula.

This is the same as finding the left-most bit in the binary representation of the number and setting all remaining bits to zero.

We can do this easily by repeatedly shifting the number to the right (the final bit drops out) and simultaneously shifting the number 1 to the left each time (this multiplies by 2, giving the next higher power of 2). Repeat until the original number is 1, and the other number will be the power of 2 we are looking for.

In F# the ">>>" operator does a bitwise right shift, and "<<<" does a left shift. So the following pair of F# functions calculate $2^{\lfloor log_2 n \rfloor}$:

let rec calcFirstBinaryDigit powerOf2 shiftedNumber =
    match shiftedNumber with
    | 0 -> raise (System.ArgumentException("No power of 2 is less than zero"))
    | 1 -> powerOf2
    | _ -> calcFirstBinaryDigit (powerOf2 <<< 1) (shiftedNumber >>> 1)

let calcHighestPowerOfTwoNotGreaterThan = calcFirstBinaryDigit 1

With this building block in place, we can easily express the binary calculation method as follows:

let calcLastLeftInCircleByBinaryFormula n =
    ((n - calcHighestPowerOfTwoNotGreaterThan n) <<< 1) + 1


The binary formula took 0.465 seconds, which is marginally faster than the floating point formula:




Comparing the results of the various methods

Now let's return to the original purpose of writing these various F# scripts, which was to check the Mathematics.

First I'm going to create a data structure to hold all the answers from the various methods for a single size of circle...

type circleCalculation = { 
    SizeOfCircle: int;
    LastLeftInCircleByBruteForce: int;
    LastLeftInCircleByHalving: int;
    LastLeftInCircleByFormula: int;
    LastLeftInCircleByBinaryFormula: int
}


Then I'm going to create a function to perform the calculation of a single record. I want this function to have a parameter allowing me to switch off the brute force method if I want to, since I know it is much slower than the other methods:

let getCircleCalculation includeBruteForce sizeOfCircle =
    // Brute force is so slow to calculate. Provide the option to omit it...
    let lastLeftInCircleByBruteForce = 
        match includeBruteForce with 
        | true -> calcLastLeftInCircleByBruteForce sizeOfCircle
        | false -> 0
    let lastLeftInCircleByHalving = calcLastLeftInCircleByHalving sizeOfCircle
    let lastLeftInCircleByFormula = calcLastLeftInCircleByFormula sizeOfCircle
    let lastLeftInCircleByBinaryFormula = calcLastLeftInCircleByBinaryFormula sizeOfCircle
    let circleCalc = {
        SizeOfCircle = sizeOfCircle;
        LastLeftInCircleByBruteForce = lastLeftInCircleByBruteForce; 
        LastLeftInCircleByHalving = lastLeftInCircleByHalving;
        LastLeftInCircleByFormula = lastLeftInCircleByFormula;
        LastLeftInCircleByBinaryFormula = lastLeftInCircleByBinaryFormula
    }
    circleCalc


Something that really amazed me is that I didn't need to specify the circleCalculation type in the function. F# is able to infer the data type from the usage pattern. Very neat!

With this in place, I can use the following script to check that the methods all match up:

// With brute force:
[1..1000] |> List.map (getCircleCalculation true) |> List.filter (
    fun cc -> cc.LastLeftInCircleByBinaryFormula <> cc.LastLeftInCircleByBruteForce
           || cc.LastLeftInCircleByBinaryFormula <> cc.LastLeftInCircleByHalving
           || cc.LastLeftInCircleByBinaryFormula <> cc.LastLeftInCircleByFormula
    )

// Without brute force:
[1..1000000] |> List.map (getCircleCalculation false) |> List.filter (
    fun cc -> cc.LastLeftInCircleByBinaryFormula <> cc.LastLeftInCircleByHalving
           || cc.LastLeftInCircleByBinaryFormula <> cc.LastLeftInCircleByFormula
    )


Running this shows that all the methods produce the same result. Mission accomplished...



Next time

In the next blog post in the series, I will show a different way of solving the problem. Instead of using the basic equations to recursively calculate a large number from the formula for smaller numbers, I'm going to work in the opposite direction.

What will emerge from this is a pattern of answers for successive values of the circle size. You might like to have a look at some of the earlier screen captures to see if you can work out what the pattern is.