🚀 go-pugleaf

Quoting quale@spock.cs.wisc.edu (Doug Quale) in article <1993Jan10.080153.10321@daffy.cs.wisc.edu>:
>In article <C0KG2y.HvE@dcs.ed.ac.uk> pdc@dcs.ed.ac.uk (Paul Crowley) writes:
>This C code does not solve the problem correctly.  It gives n**n
>possible outcomes, and in general the number of permutations of n
>objects, n!, does not divide n**n.  For example, it doesn't work
>for n = 3 (try it!).

Damn, you're right.  This is more than a little embarassing since I've
been standing by that code all last term.  However six does not divide
nine, so I guess I've got to modify it.  I think it'll work if I make a
slight change to the source array.  This time I'll prove it correct.

/* 
 * source and result are arrays of integers of length n.
 * For random permutations, source should be made up of a lot of
 * random numbers such that 0 <= source[i] <= i.
 */

void randpermute(source, result, n)
int *source, *result, n;
{
	int i;
	
	for(i = 0; i < n; i++) {
/*
 * Loop invariant: just before the test against n,
 * result[k] for 0 <= k < i is a permutation of the nat. numbers less than i
 * and the permutation is random if source is random as specified.
 * Proof below.
 */
		result[i] = result[source[i]];
		result[source[i]] = i;
	}
}

A permutation of length k is random if the probability of it being equal
to a given permutation gperm is 1/(k!).  Base case of the invariant: the
null permutation arises with probability 1 = 1/(0!).

Inductive step: supposing the invariant holds for i = k where k < n. 
After the execution of the next two source lines, is the subarray size
k+1 a permutation at all?  This divides into two subcases:

* source[k] = k

In this case, result[k] ends up being equal to k and the rest of the
array is unchanged; since the rest of the array is a permutation, this
is a permutation.

* source[k] < i

In this case, result[k] takes some value between 0 and k-1, and the
entry that used to hold that value now holds k.  Since the rest of the
array was a permutation, all the values between 0 and k are still
present and the new subarray is also a permutation.

Now, consider the element h in gperm such that gperm[h] = k.  In order
for the resulting subarray to equal gperm, two conditions must hold:
result[k] must equal h, and the old subarray before the loop stepped
around must be equal to gperm except that result[h] must equal gperm[k].
The probability of these two is 1/(k+1) and 1/(k!) respectively, so the
probability that the new subarray is gperm is 1/((k+1)!) as required.

Sorry if this isn't very clear but it's 3am where I am.  For people who
want to know what all this C is doing in comp.lang.functional, I asked
if the thing that this C function does in linear time can be done in a
pure functional style in linear time.  Except that I got the C wrong. 
Apologies, and thanks to Doug Quale for pointing it out.
  __                                   _____
\/ o\ Paul Crowley   pdc@dcs.ed.ac.uk  \\ //
/\__/ Trust me. I know what I'm doing.  \X/