Two cycles.py

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 Lecture 7 of my 2025 Oxford lectures.

1 Description
2 Program
3 Output
4 Version

[edit]

Description

To sample permutations of N elements which only contain cycles of length 1 or 2, we derived in Lecture 7 a recursion formula for the number Y_N of such permutations (eq. (7.6)), namely

Y_N = Y_{N-1} + (N - 1) Y_{N-2}

The first term on the right-hand site contains the number of permutations such that the last element is itself in a cycle of length 1, while the second term corresponds to permutations of length 2. It thus suffices to first compute the sequence of the Y_0, Y_1, Y_2, .... Y_N, then sample a random number between 0 and Y_N, and then check whether it is smaller than Y_{N-1}. If it is, the last element is in a cycle of length 1, else in one of length 2. We can continue our sampling either with N - 1 particles left in a permutation that we know nothing about, or else with N - 2 particles. All this may sound a bit mysterious, but it certainly works out OK!

The reason we discuss this algorithm is that a very similar algorithm is used for sampling cycles in the Bose-gas problem of ideal particles (for example, in a harmonic trap in three dimensions).

[edit]

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)

[edit]

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}

[edit]