Problem 367: Bozo Sort

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

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. If $\sigma'=\tau \sigma \tau^{-1}$ and $t$ is a transposition, then

$$\sigma' t = \tau \sigma (\tau^{-1} t \tau)\tau^{-1}.$$

As $\tau^{-1} t \tau$ runs through all transpositions again, the multiset of destination cycle types from $\sigma$ and from $\sigma'$ is identical. The same argument holds for oriented $3$-cycles. Therefore every permutation in the same class has the same transition probabilities between classes.

For $n=11$, the number of partitions is $p(11)=56$, so the compressed chain has only $56$ states instead of $11! = 39916800$.

Step 2: Decompose One Bozo Step

Fix three positions $i,j,k$. A uniform random permutation of those three entries has $3!=6$ possibilities. Relative to the current arrangement, these six permutations split into three types:

$$1 \text{ identity},\qquad 3 \text{ transpositions},\qquad 2 \text{ oriented }3\text{-cycles}.$$

Hence one Bozo step is equivalent to the mixture

$$\Pr(\text{no change})=\frac16,\qquad \Pr(\text{transposition})=\frac12,\qquad \Pr(\text{oriented }3\text{-cycle})=\frac13.$$

The code therefore enumerates all transpositions and all oriented $3$-cycles separately. Define

$$\mathcal T_n=\{(i\ j): 1 \le i \lt j \le n\},\qquad |\mathcal T_n|=\binom{n}{2},$$

$$\mathcal K_n=\{(i\ j\ k),(i\ k\ j): 1 \le i \lt j \lt k \le n\},\qquad |\mathcal K_n|=2\binom{n}{3}.$$

For $n=11$ this gives $55$ transpositions and $330$ oriented $3$-cycles per class representative.

Step 3: Transition Probabilities Between Classes

Choose a representative $\sigma\in \mathcal C_\lambda$. The implementations multiply on the right by every transposition and every oriented $3$-cycle, compute the cycle type of the result, and count how often each destination class $\mu$ appears. This yields the conditional class-to-class probabilities

$$P_2(\lambda,\mu)=\frac{\#\{t\in \mathcal T_n:\operatorname{type}(\sigma t)=\mu\}}{\binom{n}{2}},$$

$$P_3(\lambda,\mu)=\frac{\#\{c\in \mathcal K_n:\operatorname{type}(\sigma c)=\mu\}}{2\binom{n}{3}}.$$

Because of the conjugation argument above, these definitions do not depend on which representative of $\lambda$ is chosen.

Step 4: Linear Equations for Expected Hitting Time

Let $h_\lambda$ be the expected remaining number of Bozo steps needed to reach the identity class from any permutation of class $\lambda$. For the identity partition $1^n$ we have

$$h_{1^n}=0.$$

For every non-identity class, conditioning on the first step gives

$$h_\lambda=1+\frac16 h_\lambda+\frac12 \sum_{\mu} P_2(\lambda,\mu)h_\mu+\frac13 \sum_{\mu} P_3(\lambda,\mu)h_\mu.$$

Rearranging, exactly as in the code,

$$\boxed{\frac56 h_\lambda-\frac12 \sum_{\mu} P_2(\lambda,\mu)h_\mu-\frac13 \sum_{\mu} P_3(\lambda,\mu)h_\mu=1.}$$

There are $p(n)-1$ unknowns, since the identity state is fixed at zero. For $n=11$, the matrix dimension is $55$.

Step 5: Average Over All Permutations

Cycle types occur with different multiplicities, so the final answer is not the simple mean of the $h_\lambda$. If

$$\lambda = 1^{m_1}2^{m_2}3^{m_3}\cdots,$$

then the size of the conjugacy class is

$$|\mathcal C_\lambda|=\frac{n!}{\prod_{\ell \ge 1}\ell^{m_\ell} m_\ell!}.$$

Therefore the expectation for a uniformly random starting permutation is

$$\mathbb E[T_n]=\frac{1}{n!}\sum_{\lambda \vdash n} |\mathcal C_\lambda|\,h_\lambda.$$

This is the weighted sum computed at the end of all three language implementations.

Worked Example: $n=4$

The partitions of $4$ are $(4)$, $(3,1)$, $(2,2)$, $(2,1,1)$, and $(1,1,1,1)$. The C++ source uses $n=4$ as a checkpoint and verifies that the average expectation is $27.5$.

For the class $(2,1,1)$, the code finds

$$P_2((2,1,1),(3,1))=\frac46,\quad P_2((2,1,1),(2,2))=\frac16,\quad P_2((2,1,1),(1,1,1,1))=\frac16,$$

$$P_3((2,1,1),(4))=\frac48=\frac12,\quad P_3((2,1,1),(2,1,1))=\frac48=\frac12.$$

Substituting into the expectation equation gives

$$\frac23 h_{(2,1,1)}-\frac13 h_{(3,1)}-\frac1{12} h_{(2,2)}-\frac16 h_{(4)}=1.$$

Solving the full $4\times 4$ system yields

$$h_{(4)}=30,\qquad h_{(3,1)}=28.5,\qquad h_{(2,2)}=30,\qquad h_{(2,1,1)}=27.$$

Using class sizes $6,8,3,6,1$, the overall expectation is

$$\mathbb E[T_4]=\frac{6\cdot 30+8\cdot 28.5+3\cdot 30+6\cdot 27}{24}=27.5,$$

which matches the checkpoint in the implementation.

How the Code Works

The helper generate_partitions(n) lists all partitions of $n$. Then build_representative_permutation converts each partition into a concrete permutation with the requested cycle lengths, while cycle_type_partition recovers the sorted cycle lengths after a move. A map from partition to row index lets the program accumulate counts directly in the compressed state space.

For every non-identity class, the program enumerates all $55$ transpositions and all $330$ oriented $3$-cycles when $n=11$, fills one row of the augmented linear system, and solves the dense system by Gaussian elimination with partial pivoting. Finally class_size provides the weighting factor $ |\mathcal C_\lambda| $, the weighted mean is formed, and the returned Project Euler answer is the rounded value

$$\mathbb E[T_{11}] \approx 48271206.59545692,\qquad \text{answer}=48271207.$$

Complexity Analysis

Let $p(n)$ be the partition number. The compressed state space has $p(n)$ classes, so only $p(n)-1$ unknown expectations remain. Building one transition row requires enumerating $\binom{n}{2}$ transpositions and $2\binom{n}{3}$ oriented $3$-cycles, and each trial recomputes a cycle decomposition of a length-$n$ permutation. A faithful description of the implementation is therefore roughly

$$O\!\left(p(n)\left(\binom{n}{2}+2\binom{n}{3}\right)n\right)$$

for transition construction, followed by

$$O\!\left((p(n)-1)^3\right)$$

for dense Gaussian elimination. Memory usage is $O(p(n)^2)$ for the augmented matrix plus $O(p(n)\,n)$ for representatives and metadata. For the actual input $n=11$, this is completely practical because $p(11)=56$.

Footnotes and References

Problem page: https://projecteuler.net/problem=367
Symmetric group and conjugacy classes: Wikipedia — Symmetric group
Integer partitions: Wikipedia — Partition (number theory)
Expected hitting times in Markov chains: Wikipedia — Markov chain
Gaussian elimination: Wikipedia — Gaussian elimination

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());
    }
}