Problem 367: Bozo Sort

View on Project Euler

Project Euler Problem 367 Solution

EulerSolve provides an optimized solution for Project Euler Problem 367, Bozo Sort, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The implementation studies a Bozo-sort variant on permutations of \(\{1,2,\dots,n\}\). In one step we choose three distinct positions uniformly and then apply a uniformly random permutation to the values in those positions. The sorted permutation is the identity, and the goal is the expected number of steps needed to reach it from a uniformly random starting permutation. For the Project Euler instance the code evaluates this expectation for \(n=11\) and rounds the result to the nearest integer. Mathematical Approach A naive Markov chain on all \(n!\) permutations is already huge for \(n=11\). The key reduction in the local C++, Python, and Java solutions is that the transition rule depends only on cycle structure, so the chain can be compressed from individual permutations to conjugacy classes of \(S_n\), equivalently to integer partitions of \(n\). Step 1: Compress States by Cycle Type If a permutation has cycle type \(\lambda\), every conjugate permutation has the same cycle type. Write \(\mathcal C_\lambda\) for the conjugacy class corresponding to the partition \(\lambda \vdash n\). The program generates all partitions of \(n\), and for each partition builds one representative permutation consisting of disjoint cycles on consecutive labels. This compression is valid because the move set is itself invariant under conjugation....

Detailed mathematical approach

Problem Summary

An input sequence is one of the \(n!\) orderings of \(1,2,\dots,n\). Bozo sort checks whether it is sorted and, if not, swaps two randomly chosen elements, repeating until it is. [2] Project Euler 367 studies a variant: while the sequence is out of order we pick three elements at random and rearrange those three, all \(3!=6\) orderings equally likely. An already sorted sequence takes \(0\) steps. Over all \(4!\) inputs the statement gives \(24.75\) expected swaps for plain bozo sort and \(27.5\) expected shuffles for the variant. We must give the variant's expectation for \(n=11\), averaged over all \(11!\) inputs and rounded to the nearest integer. [1]

The whole solution in one sentence: a shuffle either changes nothing, swaps two entries or rotates three, so the expected effort depends on the input only through the lengths of its permutation cycles, which collapses the \(11!\) sequences to the \(56\) partitions of \(11\) and leaves \(55\) linear equations to solve and average.

Mathematical Approach

1. The process as a walk on the symmetric group

Start with \(n=3\), where the only three-element subset is the whole sequence, so a shuffle produces one of the six orders uniformly. Each attempt then succeeds with probability \(\tfrac16\) and the expected wait is \(6\); five of the six inputs need that wait and one is already sorted, giving the average \(5\).

For larger \(n\) we need permutations. The symmetric group \(S_n\) holds the permutations of \(n\) labels, composed as successive rearrangements. To read off the cycles of one, follow where each label goes: \((2,3,1,4)\) sends \(1\mapsto2\mapsto3\mapsto1\) and fixes \(4\), so its cycles are \((1\;2\;3)\) and \((4)\) and its cycle type is \((3,1)\); the identity has type \(1^n\).

Identify a sequence with the permutation \(\sigma\in S_n\) whose value at \(i\) is the entry in position \(i\); sorted means \(\sigma=\mathrm{id}\). Rearranging three entries touches only three positions, so the new sequence is \(\sigma\pi\) with \((\sigma\pi)(i)=\sigma(\pi(i))\), where \(\pi\) moves nothing else. Picking three values is the same as picking their three positions, the two corresponding one to one.

The process is thus \(\sigma_{t+1}=\sigma_t\pi_t\) with independent steps from one distribution \(Q\) on \(S_n\): a uniform \(3\)-element subset, then one of its \(6\) rearrangements. We stop when \(\sigma_t=\mathrm{id}\) and average that time over the \(n!\) starting points. It is finite: any permutation is sorted by at most \(n-1\) transpositions and a shuffle produces a prescribed transposition with probability \(1/[n(n-1)]\), so from every state the process finishes within \(n-1\) shuffles with probability at least \([n(n-1)]^{-(n-1)}\) and the tail decays geometrically.

2. What one shuffle really does

Of the \(6\) rearrangements of three fixed positions, one is the identity, three are transpositions (exchanging two labels) and two are three-cycles (sending three labels round one cycle), so a shuffle changes nothing with probability \(\tfrac16\). A transposition \((x\;y)\) lies in \(n-2\) of the subsets, so all transpositions are equally likely; a three-cycle determines its subset, so all three-cycles are too. Hence

$$Q=\tfrac16\,\delta_{\mathrm{id}}+\tfrac12\,U_2+\tfrac13\,U_3,$$

where \(\delta_{\mathrm{id}}\) gives probability \(1\) to the identity, \(U_2\) is uniform on the \(\binom n2\) transpositions and \(U_3\) on the \(2\binom n3\) three-cycles: \(55\) and \(330\) of them at \(n=11\).

Wasted shuffles count, and this is where a careless reading changes the answer: counting only shuffles that move something gives exactly \(\tfrac56\) of our value, each being wasted with probability \(\tfrac16\) independently. At \(n=4\) that is \(275/12=22.91\ldots\) rather than \(27.5\), so the checkpoint pins the convention down.

3. First-step analysis

Let \(h(\sigma)\) be the expected number of shuffles still needed when the sequence is \(\sigma\). Then \(h(\mathrm{id})=0\), and for \(\sigma\ne\mathrm{id}\) we spend one shuffle and land on \(\sigma\pi\): [5]

$$h(\sigma)=1+\sum_{\pi}Q(\pi)\,h(\sigma\pi).$$

Splitting \(Q\) as above and moving the \(\tfrac16 h(\sigma)\) term to the left gives, for every \(\sigma\ne\mathrm{id}\),

$$\tfrac56 h(\sigma)-\tfrac12\,\overline{h}_2(\sigma)-\tfrac13\,\overline{h}_3(\sigma)=1,$$

where \(\overline{h}_2(\sigma)\) averages \(h(\sigma t)\) over the transpositions \(t\) and \(\overline{h}_3(\sigma)\) averages \(h(\sigma c)\) over the three-cycles \(c\). That is \(n!-1\) equations in as many unknowns, which is \(39916799\) of each for \(n=11\): far too many.

4. Only the cycle type matters

Two permutations are conjugate exactly when they have the same cycle type, so the conjugacy classes of \(S_n\) correspond to the partitions of \(n\). [3] We claim \(h\) is constant on each class. Fix \(\rho\in S_n\) and relabel every state \(\sigma\) as \(\rho\sigma\rho^{-1}\). A step \(\sigma\mapsto\sigma\pi\) becomes

$$\rho\sigma\rho^{-1}\mapsto\rho\sigma\pi\rho^{-1} =(\rho\sigma\rho^{-1})(\rho\pi\rho^{-1}),$$

so the relabelled walk takes the step \(\rho\pi\rho^{-1}\). Conjugation maps the identity, the transpositions and the three-cycles onto themselves and \(Q\) is uniform on each, so \(\rho\pi\rho^{-1}\) again has distribution \(Q\), and \(\mathrm{id}\) is fixed too. The relabelled process is the same process, so \(h(\rho\sigma\rho^{-1})=h(\sigma)\), and likewise the cycle type of \(\sigma t\), counted over all transpositions \(t\), depends only on the class of \(\sigma\). One representative per class suffices.

5. The compressed linear system

Write \(h_\lambda\) for the common value on the class of type \(\lambda\), pick any representative \(\sigma\), and let

$$P_2(\lambda,\mu)=\frac{\#\{t:\operatorname{type}(\sigma t)=\mu\}}{\binom n2},$$

$$P_3(\lambda,\mu)=\frac{\#\{c:\operatorname{type}(\sigma c)=\mu\}}{2\binom n3}.$$

Section 3 becomes, for every \(\lambda\ne 1^n\),

$$\begin{aligned} &\tfrac56 h_\lambda-\tfrac12\sum_{\mu}P_2(\lambda,\mu)h_\mu\\ &\qquad-\tfrac13\sum_{\mu}P_3(\lambda,\mu)h_\mu=1, \end{aligned}$$

and the term \(\mu=1^n\) drops out because \(h_{1^n}=0\). There are \(p(n)-1\) unknowns, with \(p(n)\) the number of partitions of \(n\). [4] As \(p(11)=56\), the system for \(n=11\) has \(55\) equations: a reduction by a factor above \(700\,000\). It has one solution, since the geometric tail of Section 1 drives the powers of the restricted transition matrix \(R\) to \(0\), so \(I+R+R^2+\cdots\) converges and \(I-R\) is invertible. [5]

Only the transpositions and the three-cycles can reach the identity in one shuffle, those being the only non-identity elements \(Q\) produces: at \(n=11\) with probabilities \(\tfrac12\cdot\tfrac1{55}=\tfrac1{110}\) and \(\tfrac13\cdot\tfrac1{330}=\tfrac1{990}\) respectively, and \(0\) from every other class.

6. The worked example \(n=4\)

The five partitions of \(4\) are \((4)\), \((3,1)\), \((2,2)\), \((2,1,1)\) and \(1^4\), with class sizes \(6,8,3,6,1\); there are \(6\) transpositions and \(8\) three-cycles. Take \(\sigma=(1\;2\;3\;4)\) for the class \((4)\): \(\sigma\,(1\;2)=(1\;3\;4)(2)\) has type \((3,1)\) and \(\sigma\,(1\;3)=(1\;4)(2\;3)\) has type \((2,2)\), so the four swaps of labels adjacent on the cycle give \((3,1)\) and the two swaps of opposite labels give \((2,2)\). Counting the same way on all four representatives:

  • \((4)\): transpositions \(\to(3,1)\) four times, \((2,2)\) twice; three-cycles \(\to(4)\) and \((2,1,1)\) four times each.
  • \((3,1)\): transpositions \(\to(2,1,1)\) and \((4)\) three times each; three-cycles \(\to(3,1)\) four times, \((2,2)\) three times, \(1^4\) once.
  • \((2,2)\): transpositions \(\to(4)\) four times, \((2,1,1)\) twice; three-cycles \(\to(3,1)\) eight times.
  • \((2,1,1)\): transpositions \(\to(3,1)\) four times, \((2,2)\) once, \(1^4\) once; three-cycles \(\to(2,1,1)\) and \((4)\) four times each.

Abbreviating the unknowns as \(h_4,h_{31},h_{22},h_{211}\) and collecting coefficients:

$$\begin{aligned} \tfrac23 h_4-\tfrac13 h_{31}-\tfrac16 h_{22}-\tfrac16 h_{211}&=1,\\ \tfrac23 h_{31}-\tfrac14 h_4-\tfrac18 h_{22}-\tfrac14 h_{211}&=1,\\ \tfrac56 h_{22}-\tfrac13 h_4-\tfrac13 h_{31}-\tfrac16 h_{211}&=1,\\ \tfrac23 h_{211}-\tfrac16 h_4-\tfrac13 h_{31}-\tfrac1{12}h_{22}&=1. \end{aligned}$$

The leading coefficient stays \(\tfrac56\) only for a class that cannot return to itself, such as \((2,2)\); a three-cycle returns \((4)\) to \((4)\) in four of eight ways, moving \(\tfrac13\cdot\tfrac12=\tfrac16\) to the left. Subtracting the first row from the third gives \(h_{22}=h_4\) at once, and two more eliminations give the solution, which substitution also checks:

$$\begin{aligned} h_4=30,&\quad h_{31}=\tfrac{57}2,\\ h_{22}=30,&\quad h_{211}=27. \end{aligned}$$

7. Averaging over all inputs

A class is not a single state: if \(\lambda\) has \(m_\ell\) parts equal to \(\ell\), its size is

$$|\mathcal C_\lambda|=\frac{n!}{\prod_{\ell\ge1}\ell^{m_\ell}\,m_\ell!},$$

because we may list the \(n\) labels in any of \(n!\) orders, cut each list into blocks of the prescribed lengths and read every block as a cycle. Each permutation arises \(\ell\) times per length-\(\ell\) cycle, from the choice of starting label, and \(m_\ell!\) times from reordering the cycles of that length. [3] The quantity asked for is

$$E_n=\frac1{n!}\sum_{\lambda}|\mathcal C_\lambda|\,h_\lambda,$$

so at \(n=4\) the numerator is

$$\begin{aligned} &6\cdot30+8\cdot\tfrac{57}2\\ &\quad{}+3\cdot30+6\cdot27=660, \end{aligned}$$

giving the checkpoint \(E_4=660/24=\tfrac{55}2\). Both weightings matter: the unweighted mean of the four class values is \(28.875\), weighting by class size gives \(660/23=28.69\ldots\) over the unsorted inputs, and the one sorted input with value \(0\) lowers that to \(27.5\).

8. An independent closed form

The expectation also has a closed form, which checks everything above without a linear system. Let functions on \(S_n\) carry \(\langle f,g\rangle=\sum_\sigma f(\sigma)g(\sigma)\) and let \(P\) be the averaging operator \((Pf)(\sigma)=\sum_\pi Q(\pi)f(\sigma\pi)\), symmetric because each part of \(Q\) is closed under inverses, with rows and columns summing to \(1\). Section 3 says \((I-P)h\) equals \(1\) except at \(\mathrm{id}\), where it is \(1-c\) with \(c=1+\sum_\pi Q(\pi)h(\pi)\); summing over all \(\sigma\) gives both \(0\) and \(n!-c\), so \(c=n!\) and

$$(I-P)h=\mathbf 1-n!\,\delta_{\mathrm{id}}.$$

The irreducible representations of \(S_n\) — assignments of invertible matrices to permutations that respect composition and preserve no proper non-zero subspace — are indexed by the partitions \(\lambda\), with dimensions \(d_\lambda\), characters \(\chi_\lambda\) given by the traces, and \((n)\) trivial; their copies inside the functions on \(S_n\) span subspaces \(V_\lambda\) of dimension \(d_\lambda^2\). As \(Q\) is a class function it is central in the group algebra, so by Schur's lemma \(P\) acts on \(V_\lambda\) as one number: [6]

$$r_\lambda=\tfrac16+\tfrac12 a_\lambda+\tfrac13 b_\lambda,$$

where \(a_\lambda=\chi_\lambda(t)/d_\lambda\) and \(b_\lambda=\chi_\lambda(c)/d_\lambda\) for any transposition \(t\) and three-cycle \(c\). Column orthogonality with one argument \(\mathrm{id}\) gives \(\sum_\lambda d_\lambda\chi_\lambda(\pi)=n!\) at \(\pi=\mathrm{id}\) and \(0\) elsewhere, [6] so the \(V_\lambda\)-component of \(\delta_{\mathrm{id}}\) is \((d_\lambda/n!)\chi_\lambda\), of squared length \(d_\lambda^2/n!\). Projecting the identity above kills \(\mathbf 1\) for non-trivial \(\lambda\); pairing the result with \(\delta_{\mathrm{id}}\) and using \(h(\mathrm{id})=0\) leaves

$$E_n=\sum_{\lambda\ne(n)}\frac{d_\lambda^2}{1-r_\lambda}.$$

The denominators are positive: taking unitary matrices, \(\chi_\lambda(\pi)\le d_\lambda\) with equality only when the matrix is the identity, so \(r_\lambda\le1\), and \(r_\lambda=1\) would force \(a_\lambda=1\), making transpositions act trivially and hence the representation trivial.

Both ingredients come off the Young diagram of \(\lambda\), the left-aligned rows of its parts. The hook length \(H(i,j)\) counts cell \((i,j)\) together with the cells right of it in its row and below it in its column, and \(d_\lambda=n!/\prod H(i,j)\): for \((3,1)\) the hooks are \(4,2,1,1\), so \(d=3\). [7] Giving cell \((i,j)\) the content \(j-i\) and putting \(S_1=\sum(j-i)\), \(S_2=\sum(j-i)^2\),

$$a_\lambda=\frac{S_1}{\binom n2},\qquad b_\lambda=\frac{S_2-\binom n2}{2\binom n3}.$$

Both follow from the Jucys–Murphy elements \(X_k=\sum_{i\lt k}(i\;k)\), which multiply the basis vector of a standard Young tableau by the content of the cell holding \(k\). [8] Their sum is all the transpositions, acting by \(S_1\) with trace \(\binom n2\chi_\lambda(t)\), and

$$\sum_{k=2}^n X_k^2=\binom n2\,\mathrm{id}+\sum_c c,$$

since \((i\;k)^2=\mathrm{id}\) while the products \((i\;k)(j\;k)\) with \(i\ne j\) run once through both orientations of every three-cycle with largest label \(k\). The left side acts by \(S_2\), the right by \(\binom n2+2\binom n3 b_\lambda\).

At \(n=4\), where \(\binom42=6\) and \(2\binom43=8\), the four non-trivial \(\lambda\) have \((d,S_1,S_2)\) equal to \((3,2,6)\), \((2,0,2)\), \((3,-2,6)\) and \((1,-6,14)\), so the denominators are \(\tfrac23,1,1,1\) and the terms \(\tfrac{27}2,4,9,1\), adding to \(\tfrac{55}2=27.5\). Plain bozo sort is the same walk with \(Q=U_2\), so \(1-r_\lambda=1-a_\lambda\) and the terms become \(\tfrac{27}2,4,\tfrac{27}4,\tfrac12\), total \(\tfrac{99}4=24.75\): both numbers of the statement out of one formula. The variant is the dearer partly because a shuffle wastes a step one time in six, partly because a swap always flips the sequence's parity while a shuffle flips it only half the time, which in the sign representation is \(1-r=1\) against \(1-r=2\).

The same recipe gives \(E_5=\tfrac{2969}{20}=148.45\), \(E_6=\tfrac{225235}{252}=893.78\ldots\), \(E_7=\tfrac{3562386449}{574560}=6200.19\ldots\) and \(E_8=\tfrac{5392639268429}{109523700}=49237.19\ldots\). The large-dimensional representations carry nearly all of \(\sum_\lambda d_\lambda^2=n!\) and have small \(a_\lambda\) and \(b_\lambda\), so almost every denominator is near \(\tfrac56\) and \(E_n\) grows like a constant times \(n!\).

How the Code Works

All three build and solve the system of Section 5 inside solve_average_expectation(n), which Java and the helpers below spell in camel case. generate_partitions(n) is a depth-first search emitting every partition of \(n\) as a non-increasing list, largest remaining part first: \(56\) lists for \(n=11\). build_representative_permutation(n, partition) lays the cycles on consecutive labels, a part of length \(1\) becoming a fixed point and one of length \(\ell\ge2\) the cycle sending each label to the next and the last back to the first. cycle_type_partition(perm) inverts that, walking each cycle with a visited flag and sorting the lengths downwards to match the keys of generate_partitions; Python uses tuples for both, and class_index maps a key to a row.

One row is filled per non-identity class, the set in unknown_classes. For class \(\lambda\) the code applies every transposition to the representative by swapping two array entries, counting destinations in cnt2; swapping positions \(x\) and \(y\) in the array of \(\sigma\) yields the array of \(\sigma\,(x\;y)\), the right multiplication of Section 1. Both orientations of every three-element subset go into cnt3, the totals being trans_total\({}=\binom n2\) and cycle3_total\({}=n(n-1)(n-2)/3=2\binom n3\). Each entry becomes \(-\tfrac12P_2-\tfrac13P_3\), with \(\tfrac56\) on the diagonal and \(1\) on the right; the identity class gets no column, which is how \(h_{1^n}=0\) enters. Gaussian elimination with partial pivoting, taken to reduced row echelon form, leaves the solution in the last column. [9] Then class_size evaluates Section 7 in integer arithmetic and the weighted sum is divided by \(n!\) and rounded. The C++ run_checkpoints recomputes \(n=4\) and aborts unless it matches \(27.5\) to within \(10^{-10}\); --skip-checkpoints disables it, and the Python and Java mains have no equivalent test.

The three also differ in arithmetic and containers. C++ accumulates in long double and rounds with llround, Java uses double and Math.round, Python binary64 floats and round, which breaks exact ties to even rather than upward — harmless here, the value being far from a half-integer. The precision of long double is implementation-defined: GCC and Clang on x86-64 give it an 80-bit format, while on Apple Silicon and under Microsoft's compiler it is binary64, so that type alone guarantees nothing extra. Integers are safe at \(n=11\): C++ uses std::uint64_t and Java long for \(11!=39916800\) and the class sizes, whose denominators divide it, and Python has arbitrary precision, though neither fixed-width version checks overflow for larger \(n\). Java keys a HashMap with List<Integer> sorted by Collections.reverseOrder() and C++ an ordered std::map with std::vector<int>, both matching the generated partitions.

Complexity and Verification

Building the matrix costs \((p(n)-1)\bigl(\binom n2+2\binom n3\bigr)\) multiplications, each with an \(O(n)\) cycle walk, a sort of at most \(n\) lengths and a lookup: at \(n=11\), \(55\cdot385=21175\) cycle-type evaluations. The elimination is \(O(p(n)^3)\), of order \(55^3\approx1.7\cdot10^5\) operations, and the augmented matrix holds \(55\cdot56=3080\) numbers, about \(25\) kB in binary64. Measured here, C++ finished in under \(10\) ms and Python and Java in a few tens of milliseconds, all three printing the same integer.

The numerical margin is thinner than that size suggests. Write \(A=I-R\), with \(R\) the class-transition matrix without the identity row and column. Each row sum of \(A\) is the one-step probability of hitting the identity, so by Section 5 it is \(\tfrac1{110}\) or \(\tfrac1{990}\) for the two classes that can and \(0\) for the other \(53\) rows: \(A\) is close to \(I\) minus a stochastic matrix and nearly singular, with smallest eigenvalue about \(2.07\cdot10^{-8}\), of order \(1/n!\), and \(2\)-norm condition number about \(8.34\cdot10^7\). A condition number \(10^k\) costs up to \(k\) digits. [10] A binary64 run of this pipeline, against an exact rational solve of the same system, is off by \(0.084\), while the exact value sits \(0.275\) from the nearest rounding boundary: the answer rounds correctly, but only by a factor of about \(3.3\). The programs compute no residual or error bound, skip elimination factors below \(10^{-30}\) and check only \(n=4\), so none measures that margin.

These checks were run for this article, none of them part of the published solutions.

  • Enumerating every permutation and every allowed shuffle and solving the full chain exactly in rationals gives \(5\), \(\tfrac{55}2\) and \(\tfrac{2969}{20}\) for \(n=3,4,5\), each state value equalling its cycle-class value.
  • The same brute force on plain bozo sort gives \(\tfrac92\), \(\tfrac{99}4=24.75\) and \(\tfrac{525}4\), reproducing the other checkpoint.
  • An exact rational rerun of the published algorithm agrees with Section 8 for every \(n\) from \(3\) to \(11\), whether the character ratios come from Murnaghan–Nakayama or from contents.
  • The identity for \(\sum_k X_k^2\) holds term by term in the group algebra for \(n\le6\), the hook length dimensions satisfy \(\sum_\lambda d_\lambda^2=n!\) for \(n\le11\), and both content identities hold for all \(191\) partitions of the values \(3\le n\le11\).
  • Simulating the rule gives \(4.997\) over \(400000\) runs at \(n=3\) and \(27.571\) over \(200000\) at \(n=4\).
  • At \(n=4\) the time from \(\sigma\) to \(\mathrm{id}\) equals that from \(\mathrm{id}\) to \(\sigma^{-1}\) for all \(24\) permutations, and the average hitting time of a uniform random target is \(27.5\) from every starting state, so none is privileged. [5]
  • At \(n=11\) all \(55\) class values lie within \(1.2\%\) of one another.

Footnotes and References

The references supply the background facts; the decomposition of one shuffle, the class reduction, the compressed system, the closed form and the operation and error counts are derived here.

  1. Project Euler 367 — Bozo Sort. The official statement describes bozo sort, then the variant that picks three elements at random and shuffles them with all six orderings equally likely, says an already sorted sequence takes 0 steps, gives the averages 24.75 and 27.5 over all 4! inputs, and asks for the expected number of shuffles over all 11! inputs, rounded to the nearest integer.
  2. Wikipedia — Bogosort. Under "Related algorithms" it describes bozosort as picking two items at random, swapping them and then checking whether the list is sorted: the plain version used as a warm-up here.
  3. Wikipedia — Symmetric group. The conjugacy classes of Sₙ correspond to the cycle types of permutations and hence to the partitions of n, the basis of Section 4.
  4. OEIS A000041 — partition numbers. The number of partitions of n begins 1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42, 56, so p(11) = 56.
  5. Wikipedia — Absorbing Markov chain. The expected number of steps before absorption is N1 with fundamental matrix N = (I − Q)⁻¹, the matrix form of Section 3 and the reason the system has one solution. Wikipedia — Kemeny's constant adds that the expected number of steps to a stationary-distributed target is independent of the starting state; its diagonal term is the mean positive recurrence time, whereas this article uses 0 there, as the problem does, which lowers the constant by one and keeps the independence.
  6. Wikipedia — Schur's lemma. Its section on central characters notes that every element of the centre of the group algebra acts on an irreducible module as a scalar, which gives the number r. Wikipedia — Character theory gives the column orthogonality relation used for the point mass at the identity.
  7. Wikipedia — Hook length formula. The number of standard Young tableaux of shape λ, which equals the dimension of the matching irreducible representation of Sₙ, is n! divided by the product of the cells' hook lengths. Wikipedia — Murnaghan–Nakayama rule computes any character value of a symmetric group by removing border strips, each weighted by minus one to its height; it supplied the independent character values used above.
  8. Wikipedia — Jucys–Murphy element. It defines Xₖ as the sum of the transpositions (1 k), …, (k−1 k) and states that Xₖ acts on the basis vector of a standard Young tableau by the content b − a of the cell containing k.
  9. Wikipedia — Gaussian elimination. It describes row reduction to reduced row echelon form and partial pivoting, choosing the largest available pivot: the solver in all three programs.
  10. Wikipedia — Condition number. It bounds how far the relative error of a linear solve can exceed that of its data, with the rule of thumb that a value of 10ᵏ may cost up to k digits beyond the method's own loss.

Problem 367 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <map>
#include <numeric>
#include <string>
#include <vector>

namespace {

using u64 = std::uint64_t;

std::vector<std::vector<int>> generate_partitions(const int n) {
    std::vector<std::vector<int>> parts;
    std::vector<int> cur;

    const auto dfs = [&](const auto& self, int rem, int max_part) -> void {
        if (rem == 0) {
            parts.push_back(cur);
            return;
        }
        for (int x = std::min(rem, max_part); x >= 1; --x) {
            cur.push_back(x);
            self(self, rem - x, x);
            cur.pop_back();
        }
    };

    dfs(dfs, n, n);
    return parts;
}

std::vector<int> build_representative_permutation(const int n, const std::vector<int>& partition) {
    std::vector<int> perm(static_cast<std::size_t>(n), 0);
    int cur = 0;
    for (const int len : partition) {
        if (len == 1) {
            perm[static_cast<std::size_t>(cur)] = cur;
            ++cur;
            continue;
        }
        for (int i = 0; i < len - 1; ++i) {
            perm[static_cast<std::size_t>(cur + i)] = cur + i + 1;
        }
        perm[static_cast<std::size_t>(cur + len - 1)] = cur;
        cur += len;
    }
    return perm;
}

std::vector<int> cycle_type_partition(const std::vector<int>& perm) {
    const int n = static_cast<int>(perm.size());
    std::vector<char> vis(static_cast<std::size_t>(n), 0);
    std::vector<int> lengths;
    lengths.reserve(static_cast<std::size_t>(n));

    for (int i = 0; i < n; ++i) {
        if (vis[static_cast<std::size_t>(i)] != 0) {
            continue;
        }
        int j = i;
        int len = 0;
        while (vis[static_cast<std::size_t>(j)] == 0) {
            vis[static_cast<std::size_t>(j)] = 1;
            j = perm[static_cast<std::size_t>(j)];
            ++len;
        }
        lengths.push_back(len);
    }

    std::sort(lengths.begin(), lengths.end(), std::greater<int>());
    return lengths;
}

u64 factorial(const int n) {
    u64 v = 1;
    for (int i = 2; i <= n; ++i) {
        v *= static_cast<u64>(i);
    }
    return v;
}

u64 class_size(const int n, const std::vector<int>& partition) {
    const u64 num = factorial(n);
    std::vector<int> freq(static_cast<std::size_t>(n + 1), 0);
    for (const int len : partition) {
        ++freq[static_cast<std::size_t>(len)];
    }

    u64 den = 1;
    for (int len = 1; len <= n; ++len) {
        const int cnt = freq[static_cast<std::size_t>(len)];
        if (cnt == 0) {
            continue;
        }
        u64 p = 1;
        for (int i = 0; i < cnt; ++i) {
            p *= static_cast<u64>(len);
        }
        den *= p;
        den *= factorial(cnt);
    }

    return num / den;
}

long double solve_average_expectation(const int n) {
    const std::vector<std::vector<int>> partitions = generate_partitions(n);
    const int ccount = static_cast<int>(partitions.size());

    std::map<std::vector<int>, int> class_index;
    for (int i = 0; i < ccount; ++i) {
        class_index[partitions[static_cast<std::size_t>(i)]] = i;
    }

    std::vector<std::vector<int>> reps;
    reps.reserve(static_cast<std::size_t>(ccount));
    for (const auto& part : partitions) {
        reps.push_back(build_representative_permutation(n, part));
    }

    const std::vector<int> identity_partition(static_cast<std::size_t>(n), 1);
    const int id_idx = class_index[identity_partition];

    std::vector<int> unknown_classes;
    unknown_classes.reserve(static_cast<std::size_t>(ccount - 1));
    for (int i = 0; i < ccount; ++i) {
        if (i != id_idx) {
            unknown_classes.push_back(i);
        }
    }
    const int m = static_cast<int>(unknown_classes.size());

    const int trans_total = n * (n - 1) / 2;
    const int cycle3_total = n * (n - 1) * (n - 2) / 3;  // oriented 3-cycles: 2*C(n,3)

    std::vector<std::vector<long double>> a(static_cast<std::size_t>(m),
                                            std::vector<long double>(static_cast<std::size_t>(m + 1), 0.0L));

    for (int ri = 0; ri < m; ++ri) {
        const int cls = unknown_classes[static_cast<std::size_t>(ri)];
        const std::vector<int>& perm = reps[static_cast<std::size_t>(cls)];

        std::vector<int> cnt2(static_cast<std::size_t>(ccount), 0);
        std::vector<int> cnt3(static_cast<std::size_t>(ccount), 0);

        // Multiply on the right by transpositions.
        for (int x = 0; x < n; ++x) {
            for (int y = x + 1; y < n; ++y) {
                std::vector<int> q = perm;
                std::swap(q[static_cast<std::size_t>(x)], q[static_cast<std::size_t>(y)]);
                const int to = class_index[cycle_type_partition(q)];
                ++cnt2[static_cast<std::size_t>(to)];
            }
        }

        // Multiply on the right by oriented 3-cycles.
        for (int x = 0; x < n; ++x) {
            for (int y = x + 1; y < n; ++y) {
                for (int z = y + 1; z < n; ++z) {
                    {
                        std::vector<int> q = perm;
                        const int px = q[static_cast<std::size_t>(x)];
                        const int py = q[static_cast<std::size_t>(y)];
                        const int pz = q[static_cast<std::size_t>(z)];
                        q[static_cast<std::size_t>(x)] = py;
                        q[static_cast<std::size_t>(y)] = pz;
                        q[static_cast<std::size_t>(z)] = px;
                        const int to = class_index[cycle_type_partition(q)];
                        ++cnt3[static_cast<std::size_t>(to)];
                    }
                    {
                        std::vector<int> q = perm;
                        const int px = q[static_cast<std::size_t>(x)];
                        const int py = q[static_cast<std::size_t>(y)];
                        const int pz = q[static_cast<std::size_t>(z)];
                        q[static_cast<std::size_t>(x)] = pz;
                        q[static_cast<std::size_t>(y)] = px;
                        q[static_cast<std::size_t>(z)] = py;
                        const int to = class_index[cycle_type_partition(q)];
                        ++cnt3[static_cast<std::size_t>(to)];
                    }
                }
            }
        }

        for (int cj = 0; cj < m; ++cj) {
            const int cls_to = unknown_classes[static_cast<std::size_t>(cj)];
            const long double p2 =
                static_cast<long double>(cnt2[static_cast<std::size_t>(cls_to)]) / static_cast<long double>(trans_total);
            const long double p3 =
                static_cast<long double>(cnt3[static_cast<std::size_t>(cls_to)]) / static_cast<long double>(cycle3_total);

            long double coef = -0.5L * p2 - (1.0L / 3.0L) * p3;
            if (cls == cls_to) {
                coef += 5.0L / 6.0L;
            }
            a[static_cast<std::size_t>(ri)][static_cast<std::size_t>(cj)] = coef;
        }

        a[static_cast<std::size_t>(ri)][static_cast<std::size_t>(m)] = 1.0L;
    }

    // Gaussian elimination with partial pivoting.
    for (int col = 0; col < m; ++col) {
        int pivot = col;
        for (int r = col + 1; r < m; ++r) {
            if (std::fabsl(a[static_cast<std::size_t>(r)][static_cast<std::size_t>(col)]) >
                std::fabsl(a[static_cast<std::size_t>(pivot)][static_cast<std::size_t>(col)])) {
                pivot = r;
            }
        }
        std::swap(a[static_cast<std::size_t>(col)], a[static_cast<std::size_t>(pivot)]);

        const long double pv = a[static_cast<std::size_t>(col)][static_cast<std::size_t>(col)];
        for (int j = col; j <= m; ++j) {
            a[static_cast<std::size_t>(col)][static_cast<std::size_t>(j)] /= pv;
        }

        for (int r = 0; r < m; ++r) {
            if (r == col) {
                continue;
            }
            const long double factor = a[static_cast<std::size_t>(r)][static_cast<std::size_t>(col)];
            if (std::fabsl(factor) < 1e-30L) {
                continue;
            }
            for (int j = col; j <= m; ++j) {
                a[static_cast<std::size_t>(r)][static_cast<std::size_t>(j)] -=
                    factor * a[static_cast<std::size_t>(col)][static_cast<std::size_t>(j)];
            }
        }
    }

    std::vector<long double> h(static_cast<std::size_t>(ccount), 0.0L);
    for (int ri = 0; ri < m; ++ri) {
        const int cls = unknown_classes[static_cast<std::size_t>(ri)];
        h[static_cast<std::size_t>(cls)] = a[static_cast<std::size_t>(ri)][static_cast<std::size_t>(m)];
    }
    h[static_cast<std::size_t>(id_idx)] = 0.0L;

    const u64 total_perm = factorial(n);
    long double weighted_sum = 0.0L;
    for (int i = 0; i < ccount; ++i) {
        const u64 sz = class_size(n, partitions[static_cast<std::size_t>(i)]);
        weighted_sum += static_cast<long double>(sz) * h[static_cast<std::size_t>(i)];
    }

    return weighted_sum / static_cast<long double>(total_perm);
}

bool run_checkpoints() {
    const long double avg4 = solve_average_expectation(4);
    if (std::fabsl(avg4 - 27.5L) > 1e-10L) {
        std::cerr << "Checkpoint failed: n=4 average expectation\n";
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else {
            std::cerr << "Unknown argument: " << arg << '\n';
            return 1;
        }
    }

    if (!skip_checkpoints && !run_checkpoints()) {
        return 2;
    }

    const long double avg11 = solve_average_expectation(11);
    const long long answer = std::llround(avg11);
    std::cout << answer << '\n';
    return 0;
}

Python

import math
from collections import defaultdict

def generate_partitions(n):
    parts = []
    def dfs(rem, max_part, cur):
        if rem == 0:
            parts.append(list(cur))
            return
        for x in range(min(rem, max_part), 0, -1):
            cur.append(x)
            dfs(rem - x, x, cur)
            cur.pop()
    dfs(n, n, [])
    return parts

def build_representative_permutation(n, partition):
    perm = [0] * n
    cur = 0
    for length in partition:
        if length == 1:
            perm[cur] = cur
            cur += 1
            continue
        for i in range(length - 1):
            perm[cur + i] = cur + i + 1
        perm[cur + length - 1] = cur
        cur += length
    return perm

def cycle_type_partition(perm):
    n = len(perm)
    vis = [False] * n
    lengths = []
    for i in range(n):
        if vis[i]: continue
        j = i
        length = 0
        while not vis[j]:
            vis[j] = True
            j = perm[j]
            length += 1
        lengths.append(length)
    lengths.sort(reverse=True)
    return tuple(lengths)

def factorial(n):
    return math.factorial(n)

def class_size(n, partition):
    num = factorial(n)
    freq = defaultdict(int)
    for length in partition:
        freq[length] += 1
    
    den = 1
    for length, cnt in freq.items():
        if cnt == 0: continue
        den *= (length ** cnt)
        den *= factorial(cnt)
        
    return num // den

def solve_average_expectation(n):
    partitions = generate_partitions(n)
    ccount = len(partitions)
    
    class_index = {tuple(p): i for i, p in enumerate(partitions)}
    reps = [build_representative_permutation(n, p) for p in partitions]
    
    identity_partition = tuple([1] * n)
    id_idx = class_index[identity_partition]
    
    unknown_classes = [i for i in range(ccount) if i != id_idx]
    m = len(unknown_classes)
    
    trans_total = n * (n - 1) // 2
    cycle3_total = n * (n - 1) * (n - 2) // 3
    
    a = [[0.0] * (m + 1) for _ in range(m)]
    
    for ri in range(m):
        cls = unknown_classes[ri]
        perm = reps[cls]
        
        cnt2 = [0] * ccount
        cnt3 = [0] * ccount
        
        for x in range(n):
            for y in range(x + 1, n):
                q = list(perm)
                q[x], q[y] = q[y], q[x]
                to = class_index[cycle_type_partition(q)]
                cnt2[to] += 1
                
        for x in range(n):
            for y in range(x + 1, n):
                for z in range(y + 1, n):
                    q = list(perm)
                    px, py, pz = q[x], q[y], q[z]
                    q[x], q[y], q[z] = py, pz, px
                    to = class_index[cycle_type_partition(q)]
                    cnt3[to] += 1
                    
                    q = list(perm)
                    q[x], q[y], q[z] = pz, px, py
                    to = class_index[cycle_type_partition(q)]
                    cnt3[to] += 1
                    
        for cj in range(m):
            cls_to = unknown_classes[cj]
            p2 = cnt2[cls_to] / trans_total
            p3 = cnt3[cls_to] / cycle3_total
            
            coef = -0.5 * p2 - (1.0 / 3.0) * p3
            if cls == cls_to:
                coef += 5.0 / 6.0
            a[ri][cj] = coef
            
        a[ri][m] = 1.0
        
    for col in range(m):
        pivot = col
        for r in range(col + 1, m):
            if abs(a[r][col]) > abs(a[pivot][col]):
                pivot = r
        a[col], a[pivot] = a[pivot], a[col]
        
        pv = a[col][col]
        for j in range(col, m + 1):
            a[col][j] /= pv
            
        for r in range(m):
            if r == col: continue
            factor = a[r][col]
            if abs(factor) < 1e-30: continue
            for j in range(col, m + 1):
                a[r][j] -= factor * a[col][j]
                
    h = [0.0] * ccount
    for ri in range(m):
        cls = unknown_classes[ri]
        h[cls] = a[ri][m]
    h[id_idx] = 0.0
    
    total_perm = factorial(n)
    weighted_sum = 0.0
    for i in range(ccount):
        sz = class_size(n, partitions[i])
        weighted_sum += sz * h[i]
        
    return weighted_sum / total_perm

def solve():
    avg11 = solve_average_expectation(11)
    ans = round(avg11)
    return str(ans)

if __name__ == '__main__':
    print(solve())

Java

import java.util.*;

public class Euler367 {

    static void dfs(int rem, int maxPart, List<Integer> cur, List<List<Integer>> parts) {
        if (rem == 0) {
            parts.add(new ArrayList<>(cur));
            return;
        }
        for (int x = Math.min(rem, maxPart); x >= 1; x--) {
            cur.add(x);
            dfs(rem - x, x, cur, parts);
            cur.remove(cur.size() - 1);
        }
    }

    static List<List<Integer>> generatePartitions(int n) {
        List<List<Integer>> parts = new ArrayList<>();
        dfs(n, n, new ArrayList<>(), parts);
        return parts;
    }

    static int[] buildRepresentativePermutation(int n, List<Integer> partition) {
        int[] perm = new int[n];
        int cur = 0;
        for (int len : partition) {
            if (len == 1) {
                perm[cur] = cur;
                cur++;
                continue;
            }
            for (int i = 0; i < len - 1; i++) {
                perm[cur + i] = cur + i + 1;
            }
            perm[cur + len - 1] = cur;
            cur += len;
        }
        return perm;
    }

    static List<Integer> cycleTypePartition(int[] perm) {
        int n = perm.length;
        boolean[] vis = new boolean[n];
        List<Integer> lengths = new ArrayList<>();
        for (int i = 0; i < n; i++) {
            if (vis[i])
                continue;
            int j = i;
            int len = 0;
            while (!vis[j]) {
                vis[j] = true;
                j = perm[j];
                len++;
            }
            lengths.add(len);
        }
        lengths.sort(Collections.reverseOrder());
        return lengths;
    }

    static long factorial(int n) {
        long v = 1;
        for (int i = 2; i <= n; i++)
            v *= i;
        return v;
    }

    static long classSize(int n, List<Integer> partition) {
        long num = factorial(n);
        int[] freq = new int[n + 1];
        for (int len : partition)
            freq[len]++;

        long den = 1;
        for (int len = 1; len <= n; len++) {
            int cnt = freq[len];
            if (cnt == 0)
                continue;
            long p = 1;
            for (int i = 0; i < cnt; i++)
                p *= len;
            den *= p;
            den *= factorial(cnt);
        }
        return num / den;
    }

    static double solveAverageExpectation(int n) {
        List<List<Integer>> partitions = generatePartitions(n);
        int ccount = partitions.size();

        Map<List<Integer>, Integer> classIndex = new HashMap<>();
        for (int i = 0; i < ccount; i++) {
            classIndex.put(partitions.get(i), i);
        }

        List<int[]> reps = new ArrayList<>();
        for (List<Integer> part : partitions) {
            reps.add(buildRepresentativePermutation(n, part));
        }

        List<Integer> identityPartition = new ArrayList<>();
        for (int i = 0; i < n; i++)
            identityPartition.add(1);
        int idIdx = classIndex.get(identityPartition);

        List<Integer> unknownClasses = new ArrayList<>();
        for (int i = 0; i < ccount; i++) {
            if (i != idIdx)
                unknownClasses.add(i);
        }
        int m = unknownClasses.size();

        int transTotal = n * (n - 1) / 2;
        int cycle3Total = n * (n - 1) * (n - 2) / 3;

        double[][] a = new double[m][m + 1];

        for (int ri = 0; ri < m; ri++) {
            int cls = unknownClasses.get(ri);
            int[] perm = reps.get(cls);

            int[] cnt2 = new int[ccount];
            int[] cnt3 = new int[ccount];

            for (int x = 0; x < n; x++) {
                for (int y = x + 1; y < n; y++) {
                    int[] q = perm.clone();
                    int tmp = q[x];
                    q[x] = q[y];
                    q[y] = tmp;
                    int to = classIndex.get(cycleTypePartition(q));
                    cnt2[to]++;
                }
            }

            for (int x = 0; x < n; x++) {
                for (int y = x + 1; y < n; y++) {
                    for (int z = y + 1; z < n; z++) {
                        int[] q = perm.clone();
                        int px = q[x], py = q[y], pz = q[z];
                        q[x] = py;
                        q[y] = pz;
                        q[z] = px;
                        int to = classIndex.get(cycleTypePartition(q));
                        cnt3[to]++;

                        q = perm.clone();
                        q[x] = pz;
                        q[y] = px;
                        q[z] = py;
                        to = classIndex.get(cycleTypePartition(q));
                        cnt3[to]++;
                    }
                }
            }

            for (int cj = 0; cj < m; cj++) {
                int clsTo = unknownClasses.get(cj);
                double p2 = (double) cnt2[clsTo] / transTotal;
                double p3 = (double) cnt3[clsTo] / cycle3Total;

                double coef = -0.5 * p2 - (1.0 / 3.0) * p3;
                if (cls == clsTo)
                    coef += 5.0 / 6.0;
                a[ri][cj] = coef;
            }
            a[ri][m] = 1.0;
        }

        for (int col = 0; col < m; col++) {
            int pivot = col;
            for (int r = col + 1; r < m; r++) {
                if (Math.abs(a[r][col]) > Math.abs(a[pivot][col]))
                    pivot = r;
            }
            double[] temp = a[col];
            a[col] = a[pivot];
            a[pivot] = temp;

            double pv = a[col][col];
            for (int j = col; j <= m; j++)
                a[col][j] /= pv;

            for (int r = 0; r < m; r++) {
                if (r == col)
                    continue;
                double factor = a[r][col];
                if (Math.abs(factor) < 1e-30)
                    continue;
                for (int j = col; j <= m; j++) {
                    a[r][j] -= factor * a[col][j];
                }
            }
        }

        double[] h = new double[ccount];
        for (int ri = 0; ri < m; ri++) {
            int cls = unknownClasses.get(ri);
            h[cls] = a[ri][m];
        }
        h[idIdx] = 0.0;

        long totalPerm = factorial(n);
        double weightedSum = 0.0;
        for (int i = 0; i < ccount; i++) {
            long sz = classSize(n, partitions.get(i));
            weightedSum += sz * h[i];
        }

        return weightedSum / totalPerm;
    }

    static String solve() {
        double avg11 = solveAverageExpectation(11);
        long ans = Math.round(avg11);
        return Long.toString(ans);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}