Project Euler 367

发表于 分类于 阅读次数： 本文字数： 3.3k 阅读时长 ≈ 3 分钟

Project Euler 367

题目

Bozo sort

Bozo sort, not to be confused with the slightly less efficient bogo sort, consists out of checking if the input sequence is sorted and if not swapping randomly two elements. This is repeated until eventually the sequence is sorted.

If we consider all permutations of the first \(4\) natural numbers as input the expectation value of the number of swaps, averaged over all \(4!\) input sequences is \(24.75\).

The already sorted sequence takes \(0\) steps.

In this problem we consider the following variant on bozo sort.

If the sequence is not in order we pick three elements at random and shuffle these three elements randomly.

All \(3!=6\) permutations of those three elements are equally likely.

The already sorted sequence will take \(0\) steps.

If we consider all permutations of the first \(4\) natural numbers as input the expectation value of the number of shuffles, averaged over all \(4!\) input sequences is \(27.5\).

Consider as input sequences the permutations of the first \(11\) natural numbers.

Averaged over all \(11!\) input sequences, what is the expected number of shuffles this sorting algorithm will perform?

Give your answer rounded to the nearest integer.

解决方案

令\(n=11\)。输入序列是 \(1,2,\dots,n\) 的一个排列。用排列 \(\pi\in S_n\) 表示“位置 \(\to\) 值”的映射：$ (i)\(表示位置\)i$上的数。已经排好序当且仅当 \(\pi=\mathrm{Id}\)。其中\(\mathrm{Id}\)是指恒等置换。

一次操作：若未排序，随机选三个位点 \(i<j<k\)，并对这三个位点上的元素做等概率的 6 种重排。

把这次“重排三个位点”看成对位置施加一个排列 \(\tau\in S_n\)：\(\tau\) 在 \({i,j,k}\) 上是某个 \(S_3\) 元素，在其他位置不动。由于“新位置上的值来自旧位置”，更新是\(\pi'=\pi\circ\tau.\)

因此整个算法就是：在 \(\pi\neq\mathrm{Id}\) 时，做一次随机右乘 \(\tau\)；在 \(\pi=\mathrm{Id}\) 时吸收（0 步结束）。

排列 \(\pi\) 的不交循环分解给出若干个循环，其循环长度构成一个多重集（例如 \((5,4,2)\)）。该多重集在共轭下保持不变：若 \(\pi_2=g^{-1}\pi_1 g\)，则 \(\pi_1,\pi_2\) 的循环长度多重集相同。反过来，循环长度多重集也恰好刻画了 \(S_n\) 的共轭类。由于这些长度之和为 \(n\)，并且不计顺序，因此它等价于 \(n\) 的一个整数分拆，记作 \(\lambda\vdash n\)。我们称 \(\lambda\) 为 \(\pi\) 的循环类型。

关键在于：本题每一步随机作用的 \(\tau\)（对随机选取的三个位点做等概率重排）在分布上是共轭不变的。直观地说，这个随机操作对位置的重标号完全对称：

  • 任意固定的换位（2-循环）出现的概率相同；
  • 任意固定的 3-循环出现的概率相同；
  • 恒等置换出现的概率也相同。

因此对任意 \(g\in S_n\)，都有\(g^{-1}\tau g \stackrel{d}{=} \tau\)（其中 \(\stackrel{d}{=}\) 表示两者作为随机置换具有相同的概率分布）。

从而若 \(\pi_2=g^{-1}\pi_1 g\)（即 \(\pi_1,\pi_2\) 同循环类型），则

\[\pi_2\circ\tau= g^{-1}\pi_1 g\circ\tau= g^{-1}(\pi_1\circ(g\tau g^{-1}))g. \] 由于 \(g\tau g^{-1}\stackrel{d}{=}\tau\)，可知 \(\pi_1\circ\tau\) 与 \(\pi_2\circ\tau\) 的分布互为共轭，因此它们的循环类型分布完全相同。

因此有结论：下一步的循环类型只依赖于当前的循环类型，与同类型内选取的具体排列无关。于是我们可以把状态空间从 \(n!\) 个排列压缩为“循环类型”的个数，也就是 \(n\) 的整数分拆数，记为 \(p(n)\)。当 \(n=11\) 时，\(p(n)=56\)，因此只需在 56 个分拆类型上建立转移并求期望。

设分拆 \(\lambda\) 中长度为 \(k\) 的循环有 \(m_k\) 个（满足 \(\sum_k k m_k=n\)）。循环类型为 \(\lambda\) 的排列个数为：\(N(\lambda)=\dfrac{n!}{\prod_{k\ge 1} k^{m_k}m_k!}.\)

因为题目“对所有 \(n!\) 个输入排列等权平均”，初始在类型 \(\lambda\) 的概率为\(w(\lambda)=\dfrac{N(\lambda)}{n!}.\)已排序态对应 \(\lambda=(1,1,\dots,1)\)，其期望步数为 0。

对每个分拆 \(\lambda\vdash n\)，构造一个代表排列 \(\pi_\lambda\)：把 \(0,1,2,\dots,n-1\) 按分块组成若干个循环，每个循环用“向后轮换 1 位”实现（例如长度 \(L\) 的块形成 \((s,s+1,\dots,s+L-1)\)）。

然后对所有三元组 \({i,j,k}\)以及 \(S_3\) 的 6 种重排，生成\(\pi'=\pi_\lambda\circ\tau,\)计算其循环类型 \(\mu\)，计数累加。

总的等可能一步操作数为\(6\dbinom{n}{3}.\)因此在“循环类型”这一层，若从类型 \(\lambda\) 出发，经一步操作落到类型 \(\mu\) 的次数为 \(\#(\lambda\to\mu)\)，对应的转移概率为

\[P_{\mu,\lambda}=\frac{\#(\lambda\to\mu)}{6\binom{n}{3}}.\]

这里采用的矩阵约定是：列索引表示出发类型 \(\lambda\)，行索引表示到达类型 \(\mu\)。

令 \(t(\lambda)\) 为“当前循环类型为 \(\lambda\) 时到吸收态的期望步数”。对吸收态 \(\lambda_\star=(1^{n})\)，有 \(t(\lambda_\star)=0\)。对其他类型：

\[ t(\lambda)=1+\sum_{\mu\neq \lambda_\star} P_{\mu,\lambda}t(\mu). \]

把 \(p(n)-1\) 个非吸收态组成向量 \(t\)，令 \(Q\) 为对应的 \((p(n)-1)\times (p(n)-1)\) 子转移矩阵（排除吸收态），则

\[ (I-Q^T)t=\mathbf{1}, \]

解出 \(t\) 后，总平均期望为

\[ \mathbb{E}=\sum_{\lambda\neq \lambda_\star} w(\lambda)t(\lambda). \]

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
from fractions import Fraction
from itertools import combinations, permutations
from tools import *
from collections import Counter

N = 11
H = 3

def gen_partitions(n, max_part=None):
    """Generate all integer partitions of n as non-increasing lists."""
    if max_part is None or max_part > n:
        max_part = n
    if n == 0:
        yield []
        return
    for first in range(min(max_part, n), 0, -1):
        for rest in gen_partitions(n - first, first):
            yield [first] + rest

def rep_perm_from_partition(part, n):
    """
    Build a representative permutation π for a given cycle-type partition.
    0-indexed; permutation is represented as a list p where p[i] = π(i).
    For each block of length L>1, make a cycle rotating by 1.
    """
    p = list(range(n))
    s = 0
    for L in part:
        if L == 1:
            s += 1
            continue
        for i in range(L - 1):
            p[s + i] = s + i + 1
        p[s + L - 1] = s
        s += L
    return p

def cycle_type(p):
    """Return the cycle type (sorted tuple of cycle lengths) of permutation p."""
    n = len(p)
    vis = [False] * n
    lens = []
    for i in range(n):
        if not vis[i]:
            cur = i
            l = 0
            while not vis[cur]:
                vis[cur] = True
                cur = p[cur]
                l += 1
            lens.append(l)
    lens.sort()
    return tuple(lens)

def count_perms_of_partition(part, n):
    """Number of permutations in S_n with cycle type given by partition 'part'."""
    c = Counter(part)
    denom = 1
    for L, m in c.items():
        denom *= (L ** m) * fac(m)
    return fac(n) // denom

def gauss_solve(A, b):
    """Solve A x = b over Fractions with Gauss-Jordan elimination."""
    m = len(A)
    M = [A[i][:] + [b[i]] for i in range(m)]

    for col in range(m):
        pivot = None
        for r in range(col, m):
            if M[r][col] != 0:
                pivot = r
                break
        if pivot is None:
            raise ValueError("Singular matrix")

        if pivot != col:
            M[col], M[pivot] = M[pivot], M[col]

        pv = M[col][col]
        if pv != 1:
            inv = Fraction(1, 1) / pv
            for c in range(col, m + 1):
                M[col][c] *= inv

        for r in range(m):
            if r == col:
                continue
            factor = M[r][col]
            if factor == 0:
                continue
            for c in range(col, m + 1):
                M[r][c] -= factor * M[col][c]

    return [M[i][m] for i in range(m)]

def solve_expected_shuffles(n,h):
    # List all cycle-type partitions (as sorted tuples)
    parts = sorted({tuple(sorted(p)) for p in gen_partitions(n)}, key=lambda x: (len(x), x))
    idx = {p: i for i, p in enumerate(parts)}

    id_part = (1,) * n
    id_idx = idx[id_part]

    reps = [rep_perm_from_partition(list(p), n) for p in parts]

    total_outcomes = comb(n, h) * fac(h)  # = 6*C(n,3)

    # P[to][from] as Fractions
    P = [[Fraction(0, 1) for _ in parts] for __ in parts]

    # Build transition probabilities by enumeration on representatives
    for fr, part in enumerate(parts):
        if fr == id_idx:
            P[fr][fr] = Fraction(1, 1)  # absorbing
            continue

        base = reps[fr]
        counts = [0] * len(parts)

        for co in combinations(range(n), h):
            co = list(co)
            for perm_src in permutations(co):
                newp = base[:]  # copy
                for pos, src in zip(co, perm_src):
                    newp[pos] = base[src]
                ct = cycle_type(newp)
                counts[idx[ct]] += 1

        for to, c in enumerate(counts):
            if c:
                P[to][fr] = Fraction(c, total_outcomes)

    # Transient states exclude identity
    trans = [i for i in range(len(parts)) if i != id_idx]
    m = len(trans)

    # Solve t(fr) = 1 + sum_to P[to][fr] * t(to), for fr in transient
    # => sum_to (delta(fr,to) - P[to][fr]) * t(to) = 1
    A = [[Fraction(0, 1) for _ in range(m)] for __ in range(m)]
    b = [Fraction(1, 1) for _ in range(m)]

    for r, fr in enumerate(trans):
        for c, to in enumerate(trans):
            A[r][c] = (Fraction(1, 1) if fr == to else Fraction(0, 1)) - P[to][fr]

    t = gauss_solve(A, b)

    # Initial distribution over cycle types under uniform random permutation
    total_fact = fac(n)
    w = []
    for fr in trans:
        w.append(Fraction(count_perms_of_partition(parts[fr], n), total_fact))

    E = sum(w[i] * t[i] for i in range(m))  # exact Fraction

    # Round to nearest integer: floor(E + 1/2) since E>0
    rounded = (2 * E.numerator + E.denominator) // (2 * E.denominator)
    return E, int(rounded)

if __name__ == "__main__":
    E, rounded = solve_expected_shuffles(N,H)
    print(rounded)