Problem 903: Total Permutation Powers

Project Euler Problem 903 Solution

EulerSolve provides an optimized solution for Project Euler Problem 903, Total Permutation Powers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For $n \ge 1$, define $$Q(n)=\sum_{\pi\in S_n}\sum_{k=1}^{n!}\operatorname{rank}(\pi^k)\pmod{1\,000\,000\,007},$$ where $\operatorname{rank}(\pi^k)$ is the 1-based lexicographic rank of the permutation $\pi^k$ written in one-line notation. Problem 903 asks for the value at $n=10^6$. A direct approach would have to range over $n!$ permutations and, for each one, over $n!$ powers, which is hopeless even before ranking them. The actual solution never enumerates permutation powers. Instead, it computes the expected rank of a random power of a random permutation and then multiplies by $n!^2$. Mathematical Approach The key objects are a uniformly random permutation $\pi\in S_n$, a uniformly random exponent $i\in\{1,\dots,n!\}$, and the random power $\sigma=\pi^i$. Once we understand the pairwise comparison probabilities inside $\sigma$, the whole sum follows from the Lehmer-code formula for lexicographic rank. Random-power reformulation Choose $\pi$ uniformly from $S_n$ and choose $i$ uniformly from $\{1,\dots,n!\}$. Set $$\sigma=\pi^i.$$ Every ordered pair $(\pi,i)$ appears exactly once in the original double sum, so $$Q(n)=n!^2\,\mathbb{E}[\operatorname{rank}(\sigma)].$$ This is the first decisive simplification: the problem stops being a gigantic enumeration and becomes a single expectation problem....

Detailed mathematical approach

Problem Summary

For $n \ge 1$, define

$$Q(n)=\sum_{\pi\in S_n}\sum_{k=1}^{n!}\operatorname{rank}(\pi^k)\pmod{1\,000\,000\,007},$$

where $\operatorname{rank}(\pi^k)$ is the 1-based lexicographic rank of the permutation $\pi^k$ written in one-line notation. Problem 903 asks for the value at $n=10^6$.

A direct approach would have to range over $n!$ permutations and, for each one, over $n!$ powers, which is hopeless even before ranking them. The actual solution never enumerates permutation powers. Instead, it computes the expected rank of a random power of a random permutation and then multiplies by $n!^2$.

Mathematical Approach

The key objects are a uniformly random permutation $\pi\in S_n$, a uniformly random exponent $i\in\{1,\dots,n!\}$, and the random power $\sigma=\pi^i$. Once we understand the pairwise comparison probabilities inside $\sigma$, the whole sum follows from the Lehmer-code formula for lexicographic rank.

Random-power reformulation

Choose $\pi$ uniformly from $S_n$ and choose $i$ uniformly from $\{1,\dots,n!\}$. Set

$$\sigma=\pi^i.$$

Every ordered pair $(\pi,i)$ appears exactly once in the original double sum, so

$$Q(n)=n!^2\,\mathbb{E}[\operatorname{rank}(\sigma)].$$

This is the first decisive simplification: the problem stops being a gigantic enumeration and becomes a single expectation problem.

Lehmer digits reduce rank to pairwise comparisons

For a permutation $\tau$ on $\{1,\dots,n\}$, define its Lehmer digits by

$$a_j(\tau)=\#\{k>j:\tau(k)<\tau(j)\}.$$

The 1-based lexicographic rank is then

$$\operatorname{rank}(\tau)=1+\sum_{j=1}^{n} a_j(\tau)\,(n-j)!.$$

Applying this to $\sigma$ gives

$$\mathbb{E}[\operatorname{rank}(\sigma)]=1+\sum_{j=1}^{n}\mathbb{E}[a_j]\,(n-j)!.$$

For fixed $x<y$, write $d=y-x$ and define

$$f_d=\Pr\bigl(\sigma(x)<\sigma(y)\bigr).$$

Then

$$\mathbb{E}[a_j]=\sum_{d=1}^{n-j}\bigl(1-f_d\bigr),$$

because $a_j$ counts later positions whose values are smaller than $\sigma(j)$. So the entire problem now depends on understanding $f_d$.

Cycle probabilities for one or two marked points

The code introduces four basic probabilities for two distinct points $x\neq y$:

$$\alpha=\Pr(\sigma(x)=x),\qquad \beta=\Pr(\sigma(x)=y),$$

$$\gamma=\Pr(\sigma(x)=y,\ \sigma(y)=x),\qquad \delta=\Pr(\sigma(x)=x,\ \sigma(y)=y).$$

These numbers come directly from the cycle structure of a uniform random permutation.

A single marked point

The cycle length of a marked point in a uniformly random permutation of $S_n$ is uniform on $\{1,\dots,n\}$. If the marked point lies on a cycle of length $\ell$, then $\pi^i$ fixes that point exactly when $\ell\mid i$. Because every $\ell\le n$ divides $n!$, the uniform exponent $i\in\{1,\dots,n!\}$ satisfies this with probability $1/\ell$.

Therefore

$$\alpha=\frac1n\sum_{\ell=1}^{n}\frac1\ell=\frac{H_n}{n},\qquad H_n=\sum_{r=1}^{n}\frac1r.$$

Once $x$ is not fixed, symmetry among the remaining $n-1$ possible images gives

$$\beta=\frac{1-\alpha}{n-1}.$$

Two marked points on the same cycle or on different cycles

The swap probability $\gamma$ can happen only when $x$ and $y$ lie opposite each other on an even cycle of length $2t$. The cycle length of $x$ must be $2t$, the point $y$ must be the unique opposite point on that cycle, and the exponent must satisfy $i\equiv t\pmod{2t}$. That yields

$$\gamma=\sum_{t=1}^{\lfloor n/2\rfloor}\frac1n\cdot\frac1{n-1}\cdot\frac1{2t}=\frac{H_{\lfloor n/2\rfloor}}{2n(n-1)}.$$

For $\delta$, split according to whether $x$ and $y$ lie on the same cycle. If they share a cycle of length $\ell$, then after conditioning on the cycle length of $x$, the point $y$ lands on that same cycle with probability $(\ell-1)/(n-1)$, and both points are fixed by $\pi^i$ with probability $1/\ell$. This contributes

$$\frac1n\sum_{\ell=1}^{n}\frac{\ell-1}{n-1}\cdot\frac1\ell=\frac{n-H_n}{n(n-1)}.$$

If they lie on different cycles, there is a useful uniform law: for every ordered pair $(a,b)$ with $a,b\ge 1$ and $a+b\le n$, the event “the cycle of $x$ has length $a$, the cycle of $y$ has length $b$, and the cycles are distinct” has probability $1/(n(n-1))$. Indeed, $x$ has cycle length $a$ with probability $1/n$; conditional on that, $y$ lies outside that cycle with probability $(n-a)/(n-1)$; and among the remaining $n-a$ points, the cycle length of $y$ is uniform on $\{1,\dots,n-a\}$. In that case both points are fixed exactly when $i$ is a multiple of $\operatorname{lcm}(a,b)$, so the different-cycle contribution is

$$\frac1{n(n-1)}\sum_{a=1}^{n-1}\sum_{b=1}^{n-a}\frac{1}{\operatorname{lcm}(a,b)}=\frac{S(n)}{n(n-1)},$$

where

$$S(n)=\sum_{a=1}^{n-1}\sum_{b=1}^{n-a}\frac{1}{\operatorname{lcm}(a,b)},$$

and therefore

$$\delta=\frac{n-H_n+S(n)}{n(n-1)}.$$

Why the comparison law is affine in the distance

Now fix $x<y$ and write $d=y-x$. Count the ordered image pair $(\sigma(x),\sigma(y))$ by categories. The pair $(x,y)$ itself contributes $\delta$, while the swapped pair $(y,x)$ contributes $\gamma$ and never satisfies the inequality.

When exactly one image is fixed, each concrete choice has probability $(\alpha-\delta)/(n-2)$. If $\sigma(x)=x$, then $\sigma(y)$ may be any value larger than $x$ except $y$, giving $n-x-1$ favorable choices. If $\sigma(y)=y$, then $\sigma(x)$ may be any value smaller than $y$ except $x$, giving $y-2$ favorable choices. Altogether this category contributes

$$n-x-1+y-2=n+d-3.$$

When exactly one image hits the other marked value, each concrete choice has probability $(\beta-\gamma)/(n-2)$. If $\sigma(x)=y$, then $\sigma(y)$ must be larger than $y$, giving $n-y$ favorable choices. If $\sigma(y)=x$, then $\sigma(x)$ must be smaller than $x$, giving $x-1$ favorable choices. So this category contributes

$$n-y+x-1=n-d-1.$$

All image pairs that avoid $\{x,y\}$ altogether contribute exactly half of their total mass by symmetry. Combining the four categories gives

$$f_d=\delta+\frac{n+d-3}{n-2}(\alpha-\delta)+\frac{n-d-1}{n-2}(\beta-\gamma)+\frac{1-2\alpha-2\beta+\delta+\gamma}{2}.$$

All dependence on the actual positions has collapsed to the gap $d$, and the expression is affine:

$$f_d=A+Bd,$$

with

$$B=\frac{\alpha-\delta-\beta+\gamma}{n-2},$$

$$A=\delta+\frac{n-3}{n-2}(\alpha-\delta)+\frac{n-1}{n-2}(\beta-\gamma)+\frac{1-2\alpha-2\beta+\delta+\gamma}{2}.$$

Expected Lehmer digits

Since $a_j$ counts how many later entries are smaller than $\sigma(j)$, we sum $1-f_d$ over all later distances:

$$\mathbb{E}[a_j]=\sum_{d=1}^{n-j}\bigl(1-f_d\bigr)=(n-j)(1-A)-B\frac{(n-j)(n-j+1)}{2}.$$

The implementations use an algebraically equivalent quadratic form obtained by writing $m=j-1$:

$$\mathbb{E}[a_j]=\frac{n-H_n}{2}+m\left(\frac{H_n-1}{n-1}-A\right)-B\frac{m(m+1)}{2}.$$

That is exactly the quantity multiplied by $(n-j)!$ in the final accumulation.

Collapsing the least-common-multiple sum

The only genuinely expensive term is $S(n)$. A direct double sum would be quadratic, so the solution rewrites

$$\frac{1}{\operatorname{lcm}(a,b)}=\frac{\gcd(a,b)}{ab}=\sum_{d\mid a,\ d\mid b}\frac{\varphi(d)}{ab},$$

using the identity $\gcd(a,b)=\sum_{d\mid \gcd(a,b)}\varphi(d)$. Swapping the order of summation gives

$$S(n)=\sum_{d=1}^{n}\frac{\varphi(d)}{d^2}\,T\!\left(\left\lfloor\frac{n}{d}\right\rfloor\right),$$

where

$$T(m)=\sum_{u=1}^{m-1}\sum_{v=1}^{m-u}\frac1{uv}.$$

The inner sum is simplified by

$$\frac1{uv}=\frac{1}{u+v}\left(\frac1u+\frac1v\right),$$

which leads to the harmonic identity

$$T(m)=H_m^2-H_m^{(2)},\qquad H_m^{(2)}=\sum_{r=1}^{m}\frac1{r^2}.$$

Finally, $\left\lfloor n/d\right\rfloor$ takes only $O(\sqrt n)$ distinct values, so the sum over $d$ can be evaluated by quotient blocks once prefix sums of $\varphi(d)/d^2$ are available.

Worked example: $n=3$

For $n=3$, the harmonic values are

$$H_3=\frac{11}{6},\qquad H_{\lfloor 3/2\rfloor}=1,$$

and the least-common-multiple sum is

$$S(3)=\frac1{\operatorname{lcm}(1,1)}+\frac1{\operatorname{lcm}(1,2)}+\frac1{\operatorname{lcm}(2,1)}=1+\frac12+\frac12=2.$$

Therefore

$$\alpha=\frac{11}{18},\qquad \beta=\frac{7}{36},\qquad \gamma=\frac{1}{12},\qquad \delta=\frac{19}{36}.$$

Substituting into the affine law gives

$$A=\frac34,\qquad B=-\frac1{36},\qquad f_d=\frac34-\frac{d}{36}.$$

So the only two relevant gaps are

$$f_1=\frac{13}{18},\qquad f_2=\frac{25}{36}.$$

The expected Lehmer digits become

$$\mathbb{E}[a_1]=(1-f_1)+(1-f_2)=\frac{7}{12},\qquad \mathbb{E}[a_2]=1-f_1=\frac{5}{18}.$$

Hence

$$\mathbb{E}[\operatorname{rank}(\sigma)]=1+2!\cdot\frac{7}{12}+1!\cdot\frac{5}{18}=\frac{22}{9},$$

and finally

$$Q(3)=3!^2\cdot\frac{22}{9}=88.$$

This is the first nontrivial checkpoint used to confirm that the closed formulas agree with exact enumeration.

How the Code Works

The C++, Python, and Java implementations all follow the same derivation. Every rational quantity is represented modulo $1\,000\,000\,007$ by using modular inverses, which is valid here because the problem parameter satisfies $n=10^6<1\,000\,000\,007$.

Precompute the modular tables

The implementation first builds modular inverses of $1,2,\dots,n$. From them it forms the prefix arrays for $H_m$ and $H_m^{(2)}$. It also computes Euler's totient values up to $n$ with a linear sieve, accumulates prefix sums of $\varphi(d)/d^2$, and prepares the factorials $0!,1!,\dots,n!$ needed in the Lehmer-rank formula.

Because the affine comparison law contains the factor $1/(n-2)$, the tiny cases $n=1$ and $n=2$ are handled separately before the general branch starts. Small checkpoint values such as $Q(2)=5$ and $Q(3)=88$ provide an additional sanity test for the derivation.

Evaluate the probabilistic parameters

With the prefix data available, the implementation computes $S(n)$ by grouping equal values of $\left\lfloor n/d\right\rfloor$. It then evaluates $\alpha$, $\beta$, $\gamma$, and $\delta$, checks the one-point normalization $\alpha+(n-1)\beta=1$, converts the probabilities into the affine law $f_d=A+Bd$, and derives the quadratic formula for each expected Lehmer digit $\mathbb{E}[a_j]$.

Assemble the final answer

The last phase sums

$$1+\sum_{j=1}^{n}\mathbb{E}[a_j]\,(n-j)!$$

to obtain the expected lexicographic rank, then multiplies by $n!^2$ to recover $Q(n)$. The C++ and Java implementations parallelize only this final weighted summation across ranges of positions; the Python implementation performs the same arithmetic serially.

Complexity Analysis

Computing inverses, harmonic prefixes, totients, prefix sums, factorials, and the final expected-rank accumulation all costs $O(n)$ time. The quotient-block evaluation of $S(n)$ is only $O(\sqrt n)$ after the prefix arrays have been built, so the overall running time remains $O(n)$.

The memory usage is $O(n)$, because the method stores several arrays of length $n+1$. This is practical for $n=10^6$, whereas any attempt to enumerate permutations or powers would be completely infeasible.

Footnotes and References

Problem page: https://projecteuler.net/problem=903
Permutation and cycle notation: Wikipedia - Permutation
Lexicographic order: Wikipedia - Lexicographic order
Lehmer code: Wikipedia - Lehmer code
Harmonic number: Wikipedia - Harmonic number
Euler's totient function: Wikipedia - Euler's totient function
Least common multiple: Wikipedia - Least common multiple

Problem 903 source code

C++

#include <algorithm>
#include <cassert>
#include <iostream>
#include <numeric>
#include <stdexcept>
#include <string>
#include <thread>
#include <unordered_map>
#include <vector>
using namespace std;

static constexpr int MOD = 1'000'000'007;

static inline int addmod(int a, int b) { a += b; if (a >= MOD) a -= MOD; return a; }
static inline int submod(int a, int b) { a -= b; if (a < 0) a += MOD; return a; }
static inline int mulmod(long long a, long long b) { return int((a * b) % MOD); }

static int mod_pow(long long a, long long e) {
    long long r = 1 % MOD;
    a %= MOD;
    while (e > 0) {
        if (e & 1) r = (r * a) % MOD;
        a = (a * a) % MOD;
        e >>= 1;
    }
    return (int)r;
}
static int mod_inv(int a) { return mod_pow(a, MOD - 2); }

// inv[i] for 1..n (n < MOD)
static vector<int> compute_inverses(int n) {
    vector<int> inv(n + 1, 0);
    if (n >= 1) inv[1] = 1;
    for (int i = 2; i <= n; i++) {
        inv[i] = int(MOD - (long long)(MOD / i) * inv[MOD % i] % MOD);
    }
    return inv;
}

struct Harmonics {
    vector<int> H;   // H[k] = sum_{i=1..k} 1/i
    vector<int> H2;  // H2[k] = sum_{i=1..k} 1/i^2
};

static Harmonics compute_harmonics(int n, const vector<int>& inv) {
    Harmonics out;
    out.H.assign(n + 1, 0);
    out.H2.assign(n + 1, 0);
    for (int i = 1; i <= n; i++) {
        out.H[i]  = addmod(out.H[i - 1], inv[i]);
        out.H2[i] = addmod(out.H2[i - 1], mulmod(inv[i], inv[i]));
    }
    return out;
}

// Euler phi via linear sieve
static vector<int> compute_phi(int n) {
    vector<int> phi(n + 1, 0);
    vector<int> primes;
    primes.reserve(n / 10);
    vector<bool> isComp(n + 1, false);

    phi[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!isComp[i]) {
            primes.push_back(i);
            phi[i] = i - 1;
        }
        for (int p : primes) {
            long long v = 1LL * p * i;
            if (v > n) break;
            isComp[(int)v] = true;
            if (i % p == 0) {
                phi[(int)v] = phi[i] * p;
                break;
            } else {
                phi[(int)v] = phi[i] * (p - 1);
            }
        }
    }
    return phi;
}

// S(n) = sum_{a=1..n-1} sum_{b=1..n-a} 1/lcm(a,b)   (mod MOD)
// Uses identity:
//   S(n) = sum_{d=1..n} phi(d)/d^2 * ( H_{floor(n/d)}^2 - H2_{floor(n/d)} )
static int compute_S(int n,
                     const vector<int>& inv,
                     const vector<int>& H,
                     const vector<int>& H2,
                     const vector<int>& phi)
{
    vector<int> pref(n + 1, 0);
    for (int i = 1; i <= n; i++) {
        int inv_i2 = mulmod(inv[i], inv[i]);
        int term = mulmod(phi[i], inv_i2); // phi(i)/i^2
        pref[i] = addmod(pref[i - 1], term);
    }

    auto Tval = [&](int m) -> int {
        long long h = H[m];
        long long t = (h * h) % MOD;
        t = (t - H2[m]) % MOD;
        if (t < 0) t += MOD;
        return (int)t;
    };

    long long S = 0;
    for (int l = 1; l <= n; ) {
        int q = n / l;
        int r = n / q;
        int sumF = submod(pref[r], pref[l - 1]);
        S = (S + 1LL * sumF * Tval(q)) % MOD;
        l = r + 1;
    }
    return (int)S;
}

// Fast solver for given n (requires n < MOD).
static int solve_fast(int n) {
    if (n == 1) return 1;
    if (n == 2) return 5;
    if (n >= MOD) throw runtime_error("n must be < MOD (so denominators are invertible).");

    const auto inv = compute_inverses(n);
    const auto harm = compute_harmonics(n, inv);
    const auto phi = compute_phi(n);
    const int S = compute_S(n, inv, harm.H, harm.H2, phi);

    const int Hn = harm.H[n];
    const int inv2 = inv[2];

    const int inv_n   = inv[n];
    const int inv_nm1 = inv[n - 1];
    const int inv_d   = inv[n - 2]; // n>=3

    // Probabilities (in mod field) for σ = π^i where π uniform in S_n, i uniform in [1..n!]
    const int pF  = mulmod(Hn, inv_n);                                   // P(σ fixes a given point)
    const int pO  = mulmod(submod(1, pF), inv_nm1);                      // P(σ sends x to a specific y!=x)
    const int pSW = mulmod(mulmod(mulmod(harm.H[n / 2], inv2), inv_n), inv_nm1); // P(σ swaps two given points)
    const int pFF = mulmod(mulmod(addmod(submod(n % MOD, Hn), S), inv_n), inv_nm1); // P(σ fixes two given points)

#ifndef NDEBUG
    // Check: for fixed x, P(σ(x)=x) + sum_{y!=x} P(σ(x)=y) == 1
    long long chk = (pF + 1LL * (n - 1) % MOD * pO) % MOD;
    assert(chk == 1);
#endif

    // f(d) = P( σ(a) < σ(a+d) ) is affine in d:
    // f(d) = K0 + K1*d.
    const int K1 = mulmod(addmod(submod(submod(pF, pFF), pO), pSW), inv_d);

    int K0 = pFF;
    K0 = addmod(K0, mulmod(submod(pF, pFF), mulmod((n - 3) % MOD, inv_d)));
    K0 = addmod(K0, mulmod(submod(pO, pSW), mulmod((n - 1) % MOD, inv_d)));
    K0 = addmod(K0, inv2);
    K0 = submod(K0, pF);
    K0 = submod(K0, pO);
    K0 = addmod(K0, mulmod(addmod(pFF, pSW), inv2)); // +(pFF+pSW)/2

    // Factorials
    vector<int> fact(n + 1, 1);
    for (int i = 1; i <= n; i++) fact[i] = mulmod(fact[i - 1], i);

    // Expected Lehmer digit for position j: let m=j-1
    // E[a_{j}] = (n - Hn)/2 + m*((Hn-1)/(n-1) - K0) - K1*m(m+1)/2
    const int term1_const = mulmod(submod(n % MOD, Hn), inv2);          // (n - Hn)/2
    const int A_const     = mulmod(submod(Hn, 1), inv_nm1);             // (Hn - 1)/(n-1)
    const int C1          = submod(A_const, K0);

    // Parallel sum: Sum_{j=1..n} E[a_j] * (n-j)!  where (n-j)! = fact[n-j]
    unsigned T = max(1u, thread::hardware_concurrency());
    if (n < 100'000) T = 1; // avoid overhead for small n

    vector<long long> partial(T, 0);
    vector<thread> pool;
    pool.reserve(T);

    const int block = (n + (int)T - 1) / (int)T;

    for (unsigned t = 0; t < T; t++) {
        int L = (int)t * block + 1;
        int R = min(n, (int)(t + 1) * block);
        if (L > R) continue;

        pool.emplace_back([&, t, L, R]() {
            long long sum = 0;
            for (int j = L; j <= R; j++) {
                long long m = (long long)j - 1;
                long long mm = m % MOD;
                long long tri = mm * ((mm + 1) % MOD) % MOD;
                tri = tri * inv2 % MOD; // m(m+1)/2

                long long Ea = term1_const;
                Ea = (Ea + mm * C1) % MOD;
                Ea = (Ea - 1LL * K1 * tri) % MOD;
                if (Ea < 0) Ea += MOD;

                sum += Ea * fact[n - j] % MOD;
                if (sum > (1LL << 62)) sum %= MOD; // overflow guard
            }
            partial[t] = sum % MOD;
        });
    }

    for (auto &th : pool) th.join();

    long long sum = 0;
    for (auto v : partial) sum = (sum + v) % MOD;

    const int E_rank = addmod(1, (int)sum);
    const long long ans = 1LL * fact[n] * fact[n] % MOD * E_rank % MOD;
    return (int)ans;
}

static int brute_Q_mod(int n) {
    // Only intended for small n (<= 8).
    vector<int> p(n);
    iota(p.begin(), p.end(), 1);
    vector<vector<int>> perms;
    do perms.push_back(p);
    while (next_permutation(p.begin(), p.end()));

    int fact = 1;
    for (int i = 2; i <= n; i++) fact *= i;

    unordered_map<long long, int> rank;
    rank.reserve(perms.size() * 2);

    auto key = [&](const vector<int>& perm) -> long long {
        // n<=8 => digits are single-digit; base-10 concatenation is unique.
        long long k = 0;
        for (int x : perm) k = k * 10 + x;
        return k;
    };
    for (int i = 0; i < (int)perms.size(); i++) rank[key(perms[i])] = i + 1;

    auto compose = [&](const vector<int>& a, const vector<int>& b) -> vector<int> {
        vector<int> c(n);
        for (int i = 0; i < n; i++) c[i] = a[b[i] - 1];
        return c;
    };
    auto lcm_int = [&](int a, int b) -> int { return a / std::gcd(a, b) * b; };
    auto order = [&](const vector<int>& perm) -> int {
        vector<char> vis(n, 0);
        int ord = 1;
        for (int i = 0; i < n; i++) if (!vis[i]) {
            int len = 0, j = i;
            while (!vis[j]) {
                vis[j] = 1;
                j = perm[j] - 1;
                len++;
            }
            ord = lcm_int(ord, len);
        }
        return ord;
    };

    long long total = 0;
    for (const auto &perm : perms) {
        int m = order(perm);
        vector<int> cur = perm;
        long long s = 0;
        for (int k = 1; k <= m; k++) {
            s += rank[key(cur)];
            cur = compose(perm, cur);
        }
        total += 1LL * (fact / m) * s;
    }
    return (int)(total % MOD);
}

static void self_test() {
    // Given checkpoints
    assert(solve_fast(2)  == 5);
    assert(solve_fast(3)  == 88);
    assert(solve_fast(6)  == 133103808);
    assert(solve_fast(10) == 468421536);

    // Brute vs fast for small n
    for (int n = 1; n <= 8; n++) {
        int fast = solve_fast(n);
        int brute = brute_Q_mod(n);
        assert(fast == brute);
    }
}

int main(int argc, char** argv) {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    if (argc > 1 && string(argv[1]) == "--test") {
        self_test();
        return 0;
    }

#ifndef NDEBUG
    // Lightweight correctness checkpoints in debug builds
    self_test();
#endif

    int n = 1'000'000;
    if (argc > 1) n = stoi(argv[1]);

    cout << solve_fast(n) << "\n";
    return 0;
}

Python

def solve():
    MOD = 1000000007; n = 1000000

    def am(a, b):
        s = a+b; return s-MOD if s>=MOD else s
    def sm(a, b): return a-b if a>=b else a+MOD-b
    def mm(a, b): return a*b%MOD
    def mp(a, e):
        r = 1; a %= MOD
        while e > 0:
            if e&1: r=r*a%MOD
            a=a*a%MOD; e >>= 1
        return r
    def mi(a): return mp(a, MOD-2)

    # Inverses
    inv = [0]*(n+1)
    if n >= 1: inv[1] = 1
    for i in range(2, n+1): inv[i] = (MOD - MOD//i * inv[MOD%i] % MOD) % MOD

    # Harmonics
    H = [0]*(n+1); H2 = [0]*(n+1)
    for i in range(1, n+1):
        H[i] = am(H[i-1], inv[i]); H2[i] = am(H2[i-1], mm(inv[i], inv[i]))

    # Phi via linear sieve
    phi = [0]*(n+1); isComp = [False]*(n+1); primes = []; phi[1] = 1
    for i in range(2, n+1):
        if not isComp[i]: primes.append(i); phi[i] = i-1
        for p in primes:
            if p*i > n: break
            isComp[p*i] = True
            if i%p == 0: phi[p*i] = phi[i]*p; break
            else: phi[p*i] = phi[i]*(p-1)

    # Prefix sums of phi(i)/i^2
    pref = [0]*(n+1)
    for i in range(1, n+1):
        inv_i2 = mm(inv[i], inv[i])
        pref[i] = am(pref[i-1], mm(phi[i], inv_i2))

    def Tval(m):
        h = H[m]; t = mm(h, h); return sm(t, H2[m])

    # S(n) = sum_{d} phi(d)/d^2 * T(n/d) via block decomposition
    S = 0; l = 1
    while l <= n:
        q = n//l; r = n//q
        sumF = sm(pref[r], pref[l-1])
        S = (S + mm(sumF, Tval(q))) % MOD
        l = r+1

    Hn = H[n]; inv2 = inv[2]; inv_n = inv[n]; inv_nm1 = inv[n-1]; inv_d = inv[n-2]

    pF = mm(Hn, inv_n)
    pO = mm(sm(1, pF), inv_nm1)
    pSW = mm(mm(mm(H[n//2], inv2), inv_n), inv_nm1)
    pFF = mm(mm(am(sm(n%MOD, Hn), S), inv_n), inv_nm1)

    K1 = mm(am(sm(sm(pF, pFF), pO), pSW), inv_d)
    K0 = pFF
    K0 = am(K0, mm(sm(pF, pFF), mm((n-3)%MOD, inv_d)))
    K0 = am(K0, mm(sm(pO, pSW), mm((n-1)%MOD, inv_d)))
    K0 = am(K0, inv2); K0 = sm(K0, pF); K0 = sm(K0, pO)
    K0 = am(K0, mm(am(pFF, pSW), inv2))

    fact = [1]*(n+1)
    for i in range(1, n+1): fact[i] = mm(fact[i-1], i)

    term1_const = mm(sm(n%MOD, Hn), inv2)
    A_const = mm(sm(Hn, 1), inv_nm1)
    C1 = sm(A_const, K0)

    total = 0
    for j in range(1, n+1):
        m = j-1; mM = m%MOD
        tri = mM*(mM+1)%MOD*inv2%MOD
        Ea = (term1_const + mM*C1 - K1*tri) % MOD
        if Ea < 0: Ea += MOD
        total += Ea*fact[n-j]%MOD
        if total >= MOD*1000: total %= MOD

    total %= MOD
    E_rank = am(1, total)
    ans = fact[n]*fact[n]%MOD*E_rank%MOD
    return str(ans)

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

Java

import java.util.ArrayList;
import java.util.List;

public class Euler903 {
    static final int MOD = 1000000007;

    static int addmod(int a, int b) {
        a += b;
        if (a >= MOD)
            a -= MOD;
        return a;
    }

    static int submod(int a, int b) {
        a -= b;
        if (a < 0)
            a += MOD;
        return a;
    }

    static int mulmod(long a, long b) {
        return (int) ((a * b) % MOD);
    }

    static int modPow(long a, long e) {
        long r = 1 % MOD;
        a %= MOD;
        while (e > 0) {
            if ((e & 1) != 0)
                r = (r * a) % MOD;
            a = (a * a) % MOD;
            e >>= 1;
        }
        return (int) r;
    }

    static int modInv(int a) {
        return modPow(a, MOD - 2);
    }

    static int[] computeInverses(int n) {
        int[] inv = new int[n + 1];
        if (n >= 1)
            inv[1] = 1;
        for (int i = 2; i <= n; i++) {
            inv[i] = (int) (MOD - (long) (MOD / i) * inv[MOD % i] % MOD);
        }
        return inv;
    }

    static class Harmonics {
        int[] H;
        int[] H2;
    }

    static Harmonics computeHarmonics(int n, int[] inv) {
        Harmonics out = new Harmonics();
        out.H = new int[n + 1];
        out.H2 = new int[n + 1];
        for (int i = 1; i <= n; i++) {
            out.H[i] = addmod(out.H[i - 1], inv[i]);
            out.H2[i] = addmod(out.H2[i - 1], mulmod(inv[i], inv[i]));
        }
        return out;
    }

    static int[] computePhi(int n) {
        int[] phi = new int[n + 1];
        List<Integer> primes = new ArrayList<>();
        boolean[] isComp = new boolean[n + 1];

        phi[1] = 1;
        for (int i = 2; i <= n; i++) {
            if (!isComp[i]) {
                primes.add(i);
                phi[i] = i - 1;
            }
            for (int p : primes) {
                long v = 1L * p * i;
                if (v > n)
                    break;
                isComp[(int) v] = true;
                if (i % p == 0) {
                    phi[(int) v] = phi[i] * p;
                    break;
                } else {
                    phi[(int) v] = phi[i] * (p - 1);
                }
            }
        }
        return phi;
    }

    static int computeS(int n, int[] inv, int[] H, int[] H2, int[] phi) {
        int[] pref = new int[n + 1];
        for (int i = 1; i <= n; i++) {
            int invI2 = mulmod(inv[i], inv[i]);
            int term = mulmod(phi[i], invI2);
            pref[i] = addmod(pref[i - 1], term);
        }

        long S = 0;
        for (int l = 1; l <= n;) {
            int q = n / l;
            int r = n / q;
            int sumF = submod(pref[r], pref[l - 1]);
            long h = H[q];
            long t = (h * h) % MOD;
            t = (t - H2[q]) % MOD;
            if (t < 0)
                t += MOD;
            S = (S + 1L * sumF * t) % MOD;
            l = r + 1;
        }
        return (int) S;
    }

    static int solveFast(int n) {
        if (n == 1)
            return 1;
        if (n == 2)
            return 5;

        int[] inv = computeInverses(n);
        Harmonics harm = computeHarmonics(n, inv);
        int[] phi = computePhi(n);
        int S = computeS(n, inv, harm.H, harm.H2, phi);

        int Hn = harm.H[n];
        int inv2 = inv[2];

        int invN = inv[n];
        int invNm1 = inv[n - 1];
        int invD = inv[n - 2];

        int pF = mulmod(Hn, invN);
        int pO = mulmod(submod(1, pF), invNm1);
        int pSW = mulmod(mulmod(mulmod(harm.H[n / 2], inv2), invN), invNm1);
        int pFF = mulmod(mulmod(addmod(submod(n % MOD, Hn), S), invN), invNm1);

        int K1 = mulmod(addmod(submod(submod(pF, pFF), pO), pSW), invD);

        int K0 = pFF;
        K0 = addmod(K0, mulmod(submod(pF, pFF), mulmod((n - 3) % MOD, invD)));
        K0 = addmod(K0, mulmod(submod(pO, pSW), mulmod((n - 1) % MOD, invD)));
        K0 = addmod(K0, inv2);
        K0 = submod(K0, pF);
        K0 = submod(K0, pO);
        K0 = addmod(K0, mulmod(addmod(pFF, pSW), inv2));

        int[] fact = new int[n + 1];
        fact[0] = 1;
        for (int i = 1; i <= n; i++)
            fact[i] = mulmod(fact[i - 1], i);

        int term1Const = mulmod(submod(n % MOD, Hn), inv2);
        int AConst = mulmod(submod(Hn, 1), invNm1);
        int C1 = submod(AConst, K0);

        int T = Math.max(1, Runtime.getRuntime().availableProcessors());
        if (n < 100000)
            T = 1;

        long[] partial = new long[T];
        Thread[] pool = new Thread[T];
        int block = (n + T - 1) / T;

        for (int t = 0; t < T; t++) {
            final int tIdx = t;
            final int L = t * block + 1;
            final int R = Math.min(n, (t + 1) * block);
            if (L > R)
                continue;

            pool[t] = new Thread(() -> {
                long sum = 0;
                for (int j = L; j <= R; j++) {
                    long m = j - 1;
                    long mm = m % MOD;
                    long tri = mm * ((mm + 1) % MOD) % MOD;
                    tri = (tri * inv2) % MOD;

                    long Ea = term1Const;
                    Ea = (Ea + mm * C1) % MOD;
                    Ea = (Ea - 1L * K1 * tri) % MOD;
                    if (Ea < 0)
                        Ea += MOD;

                    sum += Ea * fact[n - j] % MOD;
                    if (sum > (1L << 60))
                        sum %= MOD;
                }
                partial[tIdx] = sum % MOD;
            });
            pool[t].start();
        }

        try {
            for (int t = 0; t < T; t++) {
                if (pool[t] != null)
                    pool[t].join();
            }
        } catch (InterruptedException e) {
            e.printStackTrace();
        }

        long sum = 0;
        for (long v : partial)
            sum = (sum + v) % MOD;

        int ERank = addmod(1, (int) sum);
        long ans = 1L * fact[n] * fact[n] % MOD * ERank % MOD;
        return (int) ans;
    }

    public static String solve() {
        return Integer.toString(solveFast(1000000));
    }

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