The Josephus Problem

Paul Wellin
Wolfram Research, Inc.

Flavius Josephus was a Jewish historian during the Roman-Jewish war of the first century.  Through his writings comes the following macabre story:

The Romans had chased a group of 10 Jews into a cave and were about to attack.  Rather than dying at the hands of their enemy, the group chose to commit suicide one by one.  Legend has it though, that they decided to go around their circle of 10 individuals and eliminate every other person until no one was left.  Who was the last to die?

A functional approach

Here is a list of 10 individuals:

Range[10]

First you must skip one,

RotateLeft[Range[10], 1]

then one must die,

Rest[RotateLeft[Range[10]]]

and repeat until only one person remains:

Nest[Rest[RotateLeft[#1, 1]] &, Range[10], 9]

Here is a listing of the survivors after each round of execution (beginning with the original orientation):

TableForm[NestList[Rest[RotateLeft[#1, 1]] &, Range[10], 9]]

The following functions generalize this problem by considering groups of n people, executing every m-th.

Josephus::"usage" = 
        "Josephus[n,m] determines who will be the last to live of a group \
of n people, numbered 1 to n if every remaining mth person were eliminated \
until only one survives.  The function outputs the lone survivor after the \
requisite executions.";

JosephusList::"usage" = 
        "JosephusList[n,m] determines who will be the last to live of a \
group of n people, numbered 1 to n if every remaining mth person were \
eliminated until only one survives.  The output consists of successive lists  \
containing the survivors after each execution.";

Josephus[n_, m_] := 
      Nest[Drop[RotateLeft[#1, 1], m - 1] &, Range[n], n - 1]

Josephus[10, 2]

JosephusList[n_, m_] := 
      TableForm[NestList[Drop[RotateLeft[#1, 1], m - 1] &, Range[n], n - 1]]

JosephusList[10, 2]

A recurrence relation

If you start out with an even number of people, after the first round of deaths, only odd-numbered people will remain.  Then the next to go will be person number 3, then 7, and so on.  But now this is just like the original group of n, except that each persons' number has been doubled and then decreased by 1.

In the case of an odd number of people, person number 1 is to die immediately after the nth.  A bit of thought should convince you that the following recurrence defines the Josephus function for n people, assuming that every second person is eliminated.

J[1] = 1;
J[n_?EvenQ] := 2  J[(n/2)] - 1 /; n ≥ 1
J[n_?OddQ] := 2  J[((n - 1)/2)] + 1 /; n ≥ 1

Here now is a table containing n and the Josephus number of n, 1  ≤ n ≤ 10:

TableForm[Table[{i, J[i]}, {i, 1, 10}]]

This is a very efficient recurrence because it reduces n by a factor of 2 each time through.  Note how quick it is in computing J[(10^9)] (how are you going to get a billion people to line up anyway especially for such an unpleasant experience?):

Timing[J[(10^9)]]

Looks like the 926,258,177th person will survive.

A closed form

A cursory look at a table of the Josephus numbers will convince you that powers of 2 are somehow involved.  See if you can use the following closed form of the Josephus function to figure out how powers of 2 are lurking behind the scenes.

Hint:  Consider any number n in the form n = 2 + q where 2 is the largest power of 2 less than or equal to n, and q is what's left over.

f[n_] := n - (2^Floor)[N[Log[2, n]]]

f[15]

J2[n_] := 2 f[n] + 1

J2[(10^9)]

Transpose[Table[{n, J2[n]}, {n, 1, 10}]]

Cyclic shifts

The Josephus number of 100 is 73.

Josephus[100, 2]

First put 100 in base 2

BaseForm[100, 2]

and then get a list of its digits:

IntegerDigits[100, 2]

Do a binary shift of one unit

l = RotateLeft[IntegerDigits[100, 2], 1]

and convert back from base 2 to base 10:

r = Reverse[l]

Dot[r, (2^Range)[0, Length[r] - 1]]

\(Josephus3[n_] := 
      Module[{p}, p = Reverse[RotateLeft[IntegerDigits[n, 2], 1]]; 
		
        Dot[p, 2\^Range[0, Length[p] - 1]]]\)

Timing[Josephus3[(10^9)]]

Functional programming the binary shift

Here is the ultimate functional approach to a binary shift.  If you look carefully, you will see Horner's method for multiplying polynomials lurking behind the scenes.  Fold is a nice little function!

Josephus4[n_Integer] := 
      Fold[2 #1 + #2 &, 0, RotateLeft[IntegerDigits[n, 2], 1]]

Timing[Josephus4[(10^9)]]

(Thanks to Paul Abbott, Richard Gaylord, Jerry Keiper, Shawn Sheridan, and Stephen Wolfram for their help with this last functional approach.)

Horner's method and the Fold function

If you are wondering what exactly the Fold function does, consider the following:

Horner[list_List, base_:  2] := Fold[base #1 + #2&, 0, list]

This implements Horner's method for fast polynomial multiplication.  Here is a fourth-order polynomial:

Horner[{a, b, c, d, e}, x]

Note there are 4 multiplies and 4 adds.  Compare this with evaluating a polynomial in the usual fashion:

FullForm[a (x^4) + b (x^3) + c (x^2) + d x + e]

Assuming for the moment that Mathematica doesn't implement Horner's rule for the Power function, you should count 10 multiplies and 4 adds.  In fact the number of multiplies for an nth degree polynomial evaluated in this manner is given by

Binomial[n + 1, 2]

This grows quadratically compared to the strictly linear growth using Horner's rule.