A random permutation of n elements can be done by randomly selecting one of the elements for the first position, one of the remaining randomly for the second position and so on. In a program, selecting from the remaining is easily done, if we store a start permutation say

\(1,\ldots,n\)

and exchange the n with a randomly chosen position, including the n-th position. Then the first n elements in the array contain the „remaining“ elements, and we can continue with the rest of the array in the same way until only one element is left.

The code for this in Euler is the following.

>function overwrite shuffle (n) ...
  p=1:n;
  for i=n to 2 step -1;
    j=intrandom(i);
    p[[i,j]]=p[[j,i]];
  end;
  return p;
endfunction

The line, where p[i] and p[j] are exchanged is typical matrix language. Probably it would be better to do it the following way.

    h=p[j]; p[j]=p[i]; p[i]=h;

Euler is not very fast in this kind of loops, and therefore there is a shuffle(n) function predefined in Euler.

But I did not want to talk about this algorithm. Rather, one might get the idea to make things a bit simpler. What if we exchange all elements with any other element randomly? The innocent change in the code above would be intrandom(n) instead of intrandom(i).

In fact, this does not work.

For an investigation of this faulty algorithm, assume

\(p_1,\ldots,p_n\)

are the probabilities of the number k to be at the positions 1 to n in one step of the algorithm. Then the probabilities after exchanging the number at i with one other position randomly are

\(\tilde p_j = \dfrac{1}{n} \, p_i + \left( 1 – \dfrac{1}{n} \right) p_j, \quad j \ne i,\)

and

\(\tilde p_i = \dfrac{1}{n}.\)

For k will remain at i, if we select the position i from the n positions with a chance of 1/n, and it can switch to position j with a chance of 1/n. Additionally, it remains at position j with a chance of (1-1/n).

Take n=3 and k=2. We start with 1,2,3 and thus with the vector of probabilities

\(0,1,0.\)

In the first step, the third position is exchanged with any of the previous. The chances for the number 2 to be at each position are then

\(0,\dfrac{2}{3},\dfrac{1}{3}.\)

In the second step, the 2 can be at all positions. The probabilities are

\(\dfrac{2}{9},\dfrac{1}{3},\dfrac{4}{9}\)

After the last step, the probabilities are

\(\dfrac{1}{3},\dfrac{8}{27},\dfrac{10}{27}\)

So we see that the chances are not evenly distributed for our number 2.

It is interesting to think about the correct algorithm in the same way. Here, we have

\(\tilde p_j = \dfrac{1}{i} \, p_i + \left( 1 – \dfrac{1}{i} \right) p_j, \quad j < i,\)

and

\(\tilde p_i = \dfrac{1}{i} (p_1+\ldots+p_i),\)

and

\(\tilde p_j = p_j, \quad j>i\)

By induction from n down to 2, we see

\(p_1+\ldots+p_i = \dfrac{i}{n}, \quad p_{i+1}=\ldots=p_n=\dfrac{1}{n}\)

 So we get the same chance for k to be at each position once we arrive at i=2.