Detailed mathematical approach

Problem Summary

Let \(m=100\) and \(n=1+2+\cdots+m=5050\). The permutation in the problem is obtained by starting with a block permutation whose disjoint cycles have lengths \(1,2,\dots,100\), and then relabeling the points by a bijection induced by multiplication by \(10^9+7\) modulo \(5050\). The task is to sum the 1-indexed lexicographic ranks of the first \(100!\) powers of that permutation and report the result modulo \(10^9+7\).

A direct simulation would require handling \(100!\) different powers of a permutation on 5050 symbols, which is completely infeasible. The key observation is that lexicographic rank can be written as a weighted inversion count, and powers of a permutation behave periodically on each pair of cycles.

Mathematical Approach

Write the final permutation as \(\rho\), in one-line notation \((\rho(1),\rho(2),\dots,\rho(n))\). The entire solution comes from combining cycle structure, Lehmer-code weights, and a counting argument on pairs of cycles.

Step 1: Identify the cycle structure

The base permutation is a disjoint union of consecutive cycles of lengths \(1,2,\dots,100\). Conjugating by a relabeling permutation changes the labels written inside the cycles, but it does not change any cycle length. Therefore \(\rho\) still has exactly one cycle of each length from 1 through 100.

This immediately implies that the order of \(\rho\) divides \(100!\), because every cycle length divides \(100!\). So when we sum over the first \(100!\) powers, every cycle has completed an integer number of full turns, and every pair of cycles has completed an integer number of joint periods.

Step 2: Rewrite lexicographic rank as weighted inversions

For any permutation \(\alpha\) on \(\{1,\dots,n\}\), its 1-indexed lexicographic rank is

$\operatorname{rank}(\alpha)=1+\sum_{i=1}^{n} c_i(\alpha)\,(n-i)!,$

where

$c_i(\alpha)=\#\{j \gt i:\alpha(j)\lt \alpha(i)\}.$

This is exactly the factorial-number-system or Lehmer-code expansion. Summing over the first \(100!\) powers of \(\rho\) gives

$S=\sum_{k=1}^{100!}\operatorname{rank}(\rho^k)=100!+\sum_{1\le i\lt j\le n}(n-i)!\,N_{i,j},$

with

$N_{i,j}=\#\{1\le k\le 100!:\rho^k(j)\lt \rho^k(i)\}.$

So the whole problem reduces to one question: for each pair of positions \(i\lt j\), how often do their images appear in reversed numerical order as the exponent runs?

Step 3: Reduce one index pair to a pair of cycles

Suppose \(i\) belongs to a cycle

$C=(c_0,c_1,\dots,c_{s-1}),$

and \(j\) belongs to a cycle

$D=(d_0,d_1,\dots,d_{t-1}).$

If \(i\) starts at position \(u\) in \(C\) and \(j\) starts at position \(v\) in \(D\), then after \(k\) applications of \(\rho\) we have

$\rho^k(i)=c_{u+k \bmod s},\qquad \rho^k(j)=d_{v+k \bmod t}.$

Now define

$g=\gcd(s,t),\qquad L=\operatorname{lcm}(s,t)=\frac{st}{g}.$

The ordered pair of cycle positions repeats every \(L\) steps. Since \(L\mid 100!\), it is enough to count favorable exponents in one \(L\)-step period and then multiply by \(100!/L\).

Step 4: Use residue classes modulo the gcd

The relative phase between the two cycles matters only modulo \(g\). Set

$\delta\equiv v-u \pmod g.$

Partition each cycle according to position modulo \(g\):

$C_r=\{c_a:a\equiv r \pmod g\},\qquad D_r=\{d_b:b\equiv r \pmod g\}.$

During one full \(L\)-step period, the pair \((\rho^k(i),\rho^k(j))\) visits exactly the pairs of positions whose residues differ by \(\delta\). Each admissible position pair occurs once. Therefore the number of favorable exponents in one period is

$Q_\delta(C,D)=\sum_{r=0}^{g-1}\#\{(x,y)\in C_r\times D_{r+\delta}: y\lt x\}.$

Hence

$N_{i,j}=Q_\delta(C,D)\cdot\frac{100!}{L}.$

The comparison problem is now finite and static: sort each residue class once, and for every phase \(\delta\) count how many elements of the second class are smaller than each element of the first.

Step 5: Assemble the final summation

Substituting the cycle-pair count into the weighted inversion formula yields

$\boxed{S=100!+\sum_{1\le i\lt j\le n}(n-i)!\,Q_{\delta(i,j)}(C(i),C(j))\frac{100!}{\operatorname{lcm}(|C(i)|,|C(j)|)} \pmod{10^9+7}.}$

Here \(C(i)\) and \(C(j)\) are the cycles containing \(i\) and \(j\), and \(\delta(i,j)\) is the phase difference of their positions modulo the relevant gcd. This is exactly the quantity accumulated by the implementations.

Worked Example: One Cycle Pair

Take two sample cycles

$C=(9,1,8,2),\qquad D=(6,3,7,4,5,10).$

Then \(s=4\), \(t=6\), so

$g=\gcd(4,6)=2,\qquad L=\operatorname{lcm}(4,6)=12.$

Split them by position parity:

$C_0=\{9,8\},\ C_1=\{1,2\},\qquad D_0=\{6,7,5\},\ D_1=\{3,4,10\}.$

After sorting,

$C_0=\{8,9\},\ C_1=\{1,2\},\qquad D_0=\{5,6,7\},\ D_1=\{3,4,10\}.$

If the starting phase difference is \(\delta=1\), then

$Q_1(C,D)=\#\{(x,y)\in C_0\times D_1:y\lt x\}+\#\{(x,y)\in C_1\times D_0:y\lt x\}=4+0=4.$

So there are exactly 4 favorable exponents in each 12-step joint period. Across the first \(100!\) powers this becomes

$N=4\cdot\frac{100!}{12}=\frac{100!}{3}.$

This miniature example is the same mechanism used for every pair of indices in the full problem.

How the Code Works

The C++, Python, and Java implementations build the relabeled permutation, decompose it into cycles, and record for every value both its cycle index and its position inside that cycle. They also precompute the factorial weights \((n-i)!\) used by the lexicographic-rank formula, together with \(100! \bmod (10^9+7)\).

For each ordered pair of cycles, the implementation computes \(g=\gcd(s,t)\), \(L=\operatorname{lcm}(s,t)\), and the phase tables \(Q_\delta(C,D)\) for \(\delta=0,1,\dots,g-1\). Because the modulus is prime and \(L<10^9+7\), the factor \(100!/L\) is represented modulo \(10^9+7\) as \(100!\cdot L^{-1}\).

After that preprocessing, the final accumulation is straightforward: for each \(i<j\), look up the two cycles, recover the correct phase \(\delta\), read the precomputed table entry, multiply by \((n-i)!\) and by the modular version of \(100!/L\), and add the result to the running sum. The C++ and Java implementations parallelize this last double loop across several threads; the Python implementation performs the same arithmetic serially.

Complexity Analysis

The permutation size is fixed at \(n=5050\), so the dominant cost is the sweep over all \(\binom{5050}{2}\) index pairs in the final accumulation. In asymptotic terms, if one writes \(n=1+2+\cdots+m=\Theta(m^2)\), then the main phase is \(O(n^2)=O(m^4)\).

The cycle-pair preprocessing is smaller: there are only 100 cycles, and each pair of cycles needs a table of size \(\gcd(s,t)\) plus sorted residue classes. Memory usage is \(O(n)\) for the permutation metadata and factorial weights, plus the cycle-pair tables, which are much smaller than the final pair sweep. In practice the algorithm is fast because it replaces gigantic exponentiation over \(100!\) powers by a single precompute-and-count pass.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=902
2. Permutation: Wikipedia - Permutation
3. Cycle notation and decomposition: Wikipedia - Cycle notation
4. Lehmer code: Wikipedia - Lehmer code
5. Lexicographic order: Wikipedia - Lexicographic order
6. Greatest common divisor: Wikipedia - Greatest common divisor
7. Least common multiple: Wikipedia - Least common multiple