Link to Mathematica script is outdated (says an advertising pest)

[Joerg Arndt]

On 6/9/06, Max <maxale at gmail.com> wrote:
> On 6/9/06, Graeme McRae <g_m at mcraefamily.com> wrote:
>
> > It took me quite some time to understand how Ranier Rosenthal is counting
> > rounds, before I obtained A112088 as the number of rounds necessary to kill
> > all but one of n players of the Josephus Game with every third man out.  I
> > first tried using Hugo Pfoertner's idea where the final survivor counts the
> > times he was passed by the executioner including the final rendezvous.  But
> > that method gives an entirely different sequence (A005428:
> > 1,2,3,4,6,9,14,21,31,47,70,105,158,...).
>
> I've also ended up with this exactly sequence in my computations.

I have found an old draft message where I was explaining the details
of my computation.
I reproduce them below in hope that someone will find them useful.

Suppose that we have x persons at the circle (enumerated from 0 to
x-1) and the executioner stands at the person number s. After one
round we will eliminate 1+[(x-s-1)/3] persons and will end up at the
person number (s-x) mod 3. This observation leads to the following
algorithm (stated in the form of PARI/GP program) counting the
required number of rounds:

{ a(n) = local(count = 0, s = 0, x = n, x_new);
while(x>1,
   x_new = x - 1 - (x-s-1)\3;
   s = (s-x) % 3;
   x = x_new;
   count++;
);
count
}

For n=1..50, PARI/GP gives:

? vector(50,n,a(n))
%1 = [0, 1, 2, 3, 3, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7,
7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9,
9, 9, 9, 9, 9, 9]

This sequence (counting the number of rounds) seems to be missing in OEIS.

Now we want to "reverse" this algorithm, i.e., starting with a single
person make the given number of rounds r and output n. This involves
solving of the following system of equations:

x_new = x_old - 1 - [(x_old - s_old - 1)/3];
s_new = (s_old - x_old) mod 3;

with respect to x_old and s_old.
It is useful to notice that
[(x_old - s_old - 1)/3] = (x_old - s_old + s_new)/3 - 1,
implying
x_new = x_old - (x_old - s_old + s_new)/3.

Therefore,

x_old = (3*x_new + s_new)/2 - s_old/2.

 From this equation we get value of s_old modulo 2 (it must make x_old
integer). However, that does not uniquely define s_old since it can
take 3 values: 0, 1, and 2. But we can follow greedy strategy here
(though I'm not going to prove that) making x_old the smallest
possible on each round. In other words, we can limit values of s_old
to 1 and 2, and state that

x_old = [(3*x_new + s_new - 1)/2].
s_old = (s_new + x_old) mod 3;

That leads to the following PARI/GP implementation:

{ b(n) = local( x=1, s=0 );
for(r=1,n,
  x = (3*x + s)\2;
  s = (s + x) % 3;
);
x
}

The first 50 values are:

? vector(50,n,b(n))
%1 = [1, 2, 3, 4, 6, 9, 14, 21, 31, 47, 70, 105, 158, 237, 355, 533,
799, 1199, 1798, 2697, 4046, 6069, 9103, 13655, 20482, 30723, 46085,
69127, 103691, 155536, 233304, 349956, 524934, 787401, 1181102,
1771653, 2657479, 3986219, 5979328, 8968992, 13453488, 20180232,
30270348, 45405522, 68108283, 102162425, 153243637, 229865456,
344798184, 517197276]

That is exactly the sequence A005428. Based on our method of
computing, we can give an alternative description of A005428:

Define a mapping f over vectors with two integer components as follows:
f: (x,s)   -->   ( [(3*x + s) / 2], (s + x) mod 3 )

Then A005428(n) is the first component of the n-th iteration of f on
the vector (1,0).

Regards,
Max