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

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