From Werner KRAUTH

This is a Python3 program to sample random permutations P of N elements subject to the constraint that the lengths of the cycles in P can only be one or two. The algorithm is based on a recursive-sampling strategy discussed in my book and presented in my 2025 Oxford lectures.

Contents 1 Description 2 Program 3 Output 4 Version

Description

Program

import random
from sympy.combinatorics import Permutation
Y = {-1: 0, 0:1, 1:1}
N = 4
Stats = {}
for k in range(2, N + 1):
    Y[k] = Y[k - 1] + (k - 1) * Y[k - 2]
for iter in range(100000):
    Q = list(range(N))
    random.shuffle(Q)
    M = N
    P = []
    while M > 0:
        if random.uniform(0.0, Y[M]) < Y[M - 1]:
            P.append([Q[M - 1]])
            M -= 1
        else:
            P.append([Q[M - 1] , Q[M - 2]])
            M -= 2
    P = tuple(Permutation(P).array_form)
    Stats[P] = Stats.get(P, 0) + 1
print(Stats)

Output

The output is a dictionary of permutations (in array form) together with the number of times they were sampled. Here we see, that all the permutations are (roughly) equally frequent, that is, were drawn with equal probabilities. We may check that all these permutations have cycles of at most two.

{(0, 3, 2, 1): 9801, (0, 2, 1, 3): 10281, (0, 1, 2, 3): 10145, (2, 3, 0, 1): 9940, (3, 1, 2, 0): 9869, 
 (0, 1, 3, 2): 9950, (2, 1, 0, 3): 9830, (3, 2, 1, 0): 10154, (1, 0, 3, 2): 10016, (1, 0, 2, 3): 10014}

Version

See history for version information. For a more elaborate version of this program, see the HistoricDisks repository connected to our paper Li et al.: Hard-disk computer simulations -- a historic perspective (2022)