#!/usr/bin/perl use strict; sub d_l_rec { # Calculates the pair (d_n, d_{n-1}) # # Why calculate the pair? It's O(n) instead of O(2^n) # for the 'obvious' recursive algorithm. # # Why recursive? No need to do it iteratively. my (\$n) = @_; return 1 if \$n < 1; return (0, 1) if \$n == 1; my (\$d1, \$d2) = d_l_rec(\$n-1); my \$d = (\$n-1) * (\$d1 + \$d2); return (\$d, \$d1); } sub random_local_derangement { # Returns a randomly-chosen local derangement of # (0..(\$n-1)). A local derangement is a derangement # *except* that the last place may be a fixed point. # # It's 'local' in the sense that, given \$n people and a # hat of \$n tickets you can generate a local derangement # by each person in turn pulling tickets out of the hat # until they get one that isn't theirs, ensuring that # (local to their draw) the permutation looks like it's # going to be a derangement. This is how 'Secret Santa' # derangements often seem to be organised, and this # process sometimes leaves the last person with their # own ticket. # # Unfortunately the 'Secret Santa' approach definitely # does not give uniform probabilities to each outcome. # This function, on the other hand, does. [At least, # it's definitely uniform for \$n <= 12 and merely a very # good approximation for larger \$n.] my (\$n) = @_; if (\$n == 0) { return []; } my (\$i, \$threshold); # A local \$n-derangement is either a full 'total' # \$n-derangement or else it is a (\$n-1)-derangement with # a fixed point added at the end. We must choose between # these options with appropriate probability weighting. # Note that this means that l_k = d_k + d_{k-1} where # l_k is the number of local k-derangements and d_k is # the number of total k-derangements. # # Note that d_{12} < 2^31 < d_{13} so we have to be # careful of overflow for \$n > 13. Fortunately there's a # good approximation that can be used. if (\$n <= 12) { # Calculate d_{\$n} and d_{\$n-1} and a random value # \$i in [0, d_{\$n} + d_{\$n-1}) to decide which of # the two options to use. my (\$dn, \$dn1) = d_l_rec(\$n); \$i = int(rand(\$dn+\$dn1)); \$threshold = \$dn1; } else { # If \$n is large then d_\$n/d_{\$n-1} is very close to # \$n. Therefore it'll do to pick a random \$i in the # range [0, \$n] and see if it is 0 or not. \$i = int(rand(\$n+1)); \$threshold = 1; # I think this is ok, but it just might contain an # off-by-one error. The upshot of such an error # would be a degree of bias in the results that is # going to be hard to detect - you may have to run # it literally trillions of times to pick up a # statistically significant result. } if (\$i < \$threshold) { # Case 1 - pick a properly local derangement my \$d = random_derangement(\$n-1); push @\$d, \$n-1; return \$d; } else { # Case 2 - pick a total derangement my \$d = random_derangement(\$n); return \$d; } } sub random_derangement { # Returns a randomly-chosen (total) derangement of # (0..(\$n-1)), uniformly-chosen amongst all possible # derangements. my (\$n) = @_; if (\$n == 0) { return []; } # There are (n-1) l_{n-1} of them, so pick a (uniformly) # random local (\$n-1)-derangement and a random \$m in the # range [0, \$n-1). my \$ld = random_local_derangement(\$n-1); my \$m = int(rand(\$n-1)); # If L_k is the set of all local k-derangements and D_k # is the set of all total k-derangements then the code # below encodes the proof that (n-1) l_{n-1} = d_n in a # bijection between [0, \$n-1) x L_{n-1} and D_{n}. # # Since the pair (\$m, \$ld) are chosen uniformly, this # shows that the resulting derangement is also uniformly # chosen. if (\$n-2 == \$ld->[\$n-2]) { # \$ld is properly local. Therefore the desired # derangement swaps the \$m'th and last places and # uses \$ld to derange the other places. my \$j = \$n-1; while (\$j--) { my \$k = \$j < \$m ? \$j : \$j-1; \$ld->[\$j] = \$ld->[\$k] < \$m ? \$ld->[\$k] : \$ld->[\$k]+1; } \$ld->[\$n-1] = \$m; \$ld->[\$m] = \$n-1; return \$ld; } else { # \$ld is total. Therefore put the \$m'th entry at the # end and put \$n-1 in the \$m'th place. \$ld->[\$n-1] = \$ld->[\$m]; \$ld->[\$m] = \$n-1; return \$ld; } } sub check_derangement { # Check that we have really generated a derangement my (\$n, \$d) = @_; my \$s = join ', ', @\$d; die "Wrong length: \$s (\$n)" unless (\$n == @\$d); for (my \$i = 0; \$i < \$n; \$i++) { die "Not a derangement: \$s" if (\$i == \$d->[\$i]); die "Illegal value: \$d->[\$i] in \$s" if (\$d->[\$i] < 0 || \$d->[\$i] >= \$n); } eval { my @check_unique = sort { \$a <=> \$b || undef } @\$d; }; die "Uniqueness check failed: \$s" if \$@; } my \$n = \$ARGV[0]; my %f = (); my \$c = 0; for (my \$i = 0; \$i < 1e6; \$i++) { my \$d = random_derangement(\$n); check_derangement(\$n, \$d); my \$s = join ',', @\$d; \$f{\$s} += 1; \$c += 1; } for my \$key (sort {\$f{\$a} <=> \$f{\$b}} keys %f) { printf "%s: %0.2f%% (expected %+0.2f%%)\n", \$key, 100.0*(\$f{\$key}/\$c), 100.0*((\$f{\$key}/\$c)-(1.0/(scalar keys %f))); } printf "Total %d (%d runs)\n", (scalar keys %f), \$c; ##```## \$ ./derangements.pl 4 2,3,0,1: 11.08% (expected -0.03%) 1,0,3,2: 11.09% (expected -0.02%) 2,3,1,0: 11.10% (expected -0.01%) 1,3,0,2: 11.11% (expected -0.01%) 2,0,3,1: 11.11% (expected -0.00%) 1,2,3,0: 11.12% (expected +0.01%) 3,0,1,2: 11.12% (expected +0.01%) 3,2,1,0: 11.13% (expected +0.02%) 3,2,0,1: 11.14% (expected +0.03%) Total 9 (1000000 runs) ```