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\).
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.
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.
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\).
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.
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}.$$
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)}.$$
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}.$$
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.
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.
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.
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\).
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.
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]\).
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.
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.
Für \(n \ge 1\) sei
$$Q(n)=\sum_{\pi\in S_n}\sum_{k=1}^{n!}\operatorname{rank}(\pi^k)\pmod{1\,000\,000\,007},$$
wobei \(\operatorname{rank}(\pi^k)\) den 1-basierten lexikographischen Rang der Potenz \(\pi^k\) in Einzeilen-Schreibweise bezeichnet. In Problem 903 soll der Wert für \(n=10^6\) bestimmt werden.
Ein direkter Ansatz wäre aussichtslos: Man müsste über \(n!\) Permutationen laufen und für jede davon noch \(n!\) Potenzen betrachten. Die wirkliche Lösung erzeugt keine Potenzen explizit. Stattdessen berechnet sie den erwarteten Rang einer zufälligen Potenz einer zufälligen Permutation und multipliziert am Ende mit \(n!^2\).
Die zentralen Objekte sind eine gleichverteilte Zufallspermutation \(\pi\in S_n\), ein gleichverteilter Exponent \(i\in\{1,\dots,n!\}\) und die Zufallspermutation \(\sigma=\pi^i\). Sobald die paarweisen Vergleichswahrscheinlichkeiten in \(\sigma\) bekannt sind, folgt die gesamte Summe aus der Lehmer-Code-Darstellung des lexikographischen Rangs.
Wähle \(\pi\) gleichverteilt aus \(S_n\) und \(i\) gleichverteilt aus \(\{1,\dots,n!\}\). Setze
$$\sigma=\pi^i.$$
Jedes geordnete Paar \((\pi,i)\) erscheint genau einmal in der ursprünglichen Doppelsumme. Daher gilt
$$Q(n)=n!^2\,\mathbb{E}[\operatorname{rank}(\sigma)].$$
Damit wird aus einer riesigen Enumeration ein einziges Erwartungswertproblem.
Für eine Permutation \(\tau\) auf \(\{1,\dots,n\}\) definiere
$$a_j(\tau)=\#\{k>j:\tau(k)<\tau(j)\}.$$
Dann ist der 1-basierte lexikographische Rang
$$\operatorname{rank}(\tau)=1+\sum_{j=1}^{n} a_j(\tau)\,(n-j)!.$$
Für \(\sigma\) folgt also
$$\mathbb{E}[\operatorname{rank}(\sigma)]=1+\sum_{j=1}^{n}\mathbb{E}[a_j]\,(n-j)!.$$
Für feste Indizes \(x<y\) schreiben wir \(d=y-x\) und definieren
$$f_d=\Pr\bigl(\sigma(x)<\sigma(y)\bigr).$$
Dann gilt
$$\mathbb{E}[a_j]=\sum_{d=1}^{n-j}\bigl(1-f_d\bigr),$$
denn \(a_j\) zählt genau die späteren Positionen mit kleinerem Bildwert. Das ganze Problem hängt also an einer expliziten Formel für \(f_d\).
Die Implementierungen verwenden für zwei verschiedene Punkte \(x\neq y\) vier Grundwahrscheinlichkeiten:
$$\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).$$
Alle vier Größen folgen unmittelbar aus der Zyklusstruktur einer gleichverteilten Zufallspermutation.
Die Zykluslänge eines markierten Punkts in einer gleichverteilten Permutation aus \(S_n\) ist gleichverteilt auf \(\{1,\dots,n\}\). Liegt der Punkt in einem Zyklus der Länge \(\ell\), dann wird er von \(\pi^i\) genau dann fixiert, wenn \(\ell\mid i\). Da jedes \(\ell\le n\) ein Teiler von \(n!\) ist, hat dieses Teilbarkeitsereignis für ein gleichverteiltes \(i\in\{1,\dots,n!\}\) die Wahrscheinlichkeit \(1/\ell\).
Daher
$$\alpha=\frac1n\sum_{\ell=1}^{n}\frac1\ell=\frac{H_n}{n},\qquad H_n=\sum_{r=1}^{n}\frac1r.$$
Ist \(x\) nicht fixiert, dann sind die übrigen \(n-1\) möglichen Bilder symmetrisch. Also
$$\beta=\frac{1-\alpha}{n-1}.$$
Die Vertauschungswahrscheinlichkeit \(\gamma\) kann nur auftreten, wenn \(x\) und \(y\) einander in einem geraden Zyklus der Länge \(2t\) genau gegenüberliegen. Die Zykluslänge von \(x\) muss also \(2t\) sein, \(y\) muss der eindeutig gegenüberliegende Punkt sein, und der Exponent muss \(i\equiv t\pmod{2t}\) erfüllen. Damit erhält man
$$\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)}.$$
Für \(\delta\) zerlegt man nach dem Fall, ob \(x\) und \(y\) im selben Zyklus liegen. Liegen sie in einem gemeinsamen Zyklus der Länge \(\ell\), dann liegt \(y\) nach Fixierung der Zykluslänge von \(x\) mit Wahrscheinlichkeit \((\ell-1)/(n-1)\) in demselben Zyklus, und beide Punkte werden mit Wahrscheinlichkeit \(1/\ell\) von \(\pi^i\) fixiert. Das liefert
$$\frac1n\sum_{\ell=1}^{n}\frac{\ell-1}{n-1}\cdot\frac1\ell=\frac{n-H_n}{n(n-1)}.$$
Liegen die beiden Punkte in verschiedenen Zyklen, gibt es ein nützliches Gleichverteilungsgesetz: Für jedes geordnete Paar \((a,b)\) mit \(a,b\ge 1\) und \(a+b\le n\) hat das Ereignis „der Zyklus von \(x\) hat Länge \(a\), der Zyklus von \(y\) hat Länge \(b\), und beide Zyklen sind verschieden“ die Wahrscheinlichkeit \(1/(n(n-1))\). Denn \(x\) hat Zykluslänge \(a\) mit Wahrscheinlichkeit \(1/n\); bedingt darauf liegt \(y\) mit Wahrscheinlichkeit \((n-a)/(n-1)\) außerhalb dieses Zyklus; und auf den verbleibenden \(n-a\) Punkten ist die Zykluslänge von \(y\) gleichverteilt auf \(\{1,\dots,n-a\}\). In diesem Fall werden beide Punkte genau dann fixiert, wenn \(i\) ein Vielfaches von \(\operatorname{lcm}(a,b)\) ist. Der Beitrag der verschiedenen Zyklen lautet also
$$\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)},$$
wobei
$$S(n)=\sum_{a=1}^{n-1}\sum_{b=1}^{n-a}\frac{1}{\operatorname{lcm}(a,b)},$$
und damit
$$\delta=\frac{n-H_n+S(n)}{n(n-1)}.$$
Fixiere nun \(x<y\) und schreibe \(d=y-x\). Zähle die geordnete Bildpaarung \((\sigma(x),\sigma(y))\) nach Kategorien. Das Paar \((x,y)\) selbst liefert den Beitrag \(\delta\), während die Vertauschung \((y,x)\) den Beitrag \(\gamma\) liefert und die Ungleichung nie erfüllt.
Wenn genau ein Bild fixiert ist, hat jede konkrete Wahl die Wahrscheinlichkeit \((\alpha-\delta)/(n-2)\). Im Fall \(\sigma(x)=x\) darf \(\sigma(y)\) jeder Wert größer als \(x\) außer \(y\) sein, also gibt es \(n-x-1\) günstige Möglichkeiten. Im Fall \(\sigma(y)=y\) darf \(\sigma(x)\) jeder Wert kleiner als \(y\) außer \(x\) sein, also \(y-2\) Möglichkeiten. Zusammen ergibt diese Kategorie
$$n-x-1+y-2=n+d-3.$$
Wenn genau ein Bild auf den anderen markierten Wert fällt, hat jede konkrete Wahl die Wahrscheinlichkeit \((\beta-\gamma)/(n-2)\). Im Fall \(\sigma(x)=y\) muss \(\sigma(y)\) größer als \(y\) sein, also gibt es \(n-y\) Möglichkeiten. Im Fall \(\sigma(y)=x\) muss \(\sigma(x)\) kleiner als \(x\) sein, also \(x-1\) Möglichkeiten. Diese Kategorie liefert daher
$$n-y+x-1=n-d-1.$$
Alle Bildpaare, die \(\{x,y\}\) ganz vermeiden, tragen wegen der Symmetrie genau die Hälfte ihrer Gesamtmasse bei. Kombiniert man diese vier Kategorien, so erhält man
$$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}.$$
Die konkrete Position spielt also nur noch über den Abstand \(d\) eine Rolle, und der Ausdruck ist affin:
$$f_d=A+Bd,$$
mit
$$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}.$$
Da \(a_j\) die Zahl späterer Werte zählt, die kleiner sind als \(\sigma(j)\), summieren wir \(1-f_d\) über alle späteren Abstände:
$$\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}.$$
Die Implementierungen verwenden die dazu äquivalente quadratische Form mit \(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}.$$
Genau dieser Ausdruck wird in der Endsumme mit \((n-j)!\) gewichtet.
Der einzige wirklich teure Term ist \(S(n)\). Eine direkte Doppelsumme wäre quadratisch, daher schreibt die Lösung
$$\frac{1}{\operatorname{lcm}(a,b)}=\frac{\gcd(a,b)}{ab}=\sum_{d\mid a,\ d\mid b}\frac{\varphi(d)}{ab},$$
wobei \(\gcd(a,b)=\sum_{d\mid \gcd(a,b)}\varphi(d)\) benutzt wird. Nach Vertauschen der Summationsreihenfolge folgt
$$S(n)=\sum_{d=1}^{n}\frac{\varphi(d)}{d^2}\,T\!\left(\left\lfloor\frac{n}{d}\right\rfloor\right),$$
mit
$$T(m)=\sum_{u=1}^{m-1}\sum_{v=1}^{m-u}\frac1{uv}.$$
Die innere Summe wird durch
$$\frac1{uv}=\frac{1}{u+v}\left(\frac1u+\frac1v\right)$$
auf die harmonische Identität
$$T(m)=H_m^2-H_m^{(2)},\qquad H_m^{(2)}=\sum_{r=1}^{m}\frac1{r^2}$$
zurückgeführt. Schließlich nimmt \(\left\lfloor n/d\right\rfloor\) nur \(O(\sqrt n)\) verschiedene Werte an, sodass die Summe über \(d\) blockweise nach gleichen Quotienten ausgewertet werden kann.
Für \(n=3\) gelten die harmonischen Werte
$$H_3=\frac{11}{6},\qquad H_{\lfloor 3/2\rfloor}=1,$$
und die kgV-Summe ist
$$S(3)=\frac1{\operatorname{lcm}(1,1)}+\frac1{\operatorname{lcm}(1,2)}+\frac1{\operatorname{lcm}(2,1)}=1+\frac12+\frac12=2.$$
Daraus folgen
$$\alpha=\frac{11}{18},\qquad \beta=\frac{7}{36},\qquad \gamma=\frac{1}{12},\qquad \delta=\frac{19}{36}.$$
In das affine Gesetz eingesetzt ergibt das
$$A=\frac34,\qquad B=-\frac1{36},\qquad f_d=\frac34-\frac{d}{36}.$$
Damit sind die beiden relevanten Abstände
$$f_1=\frac{13}{18},\qquad f_2=\frac{25}{36}.$$
Die erwarteten Lehmer-Ziffern lauten also
$$\mathbb{E}[a_1]=(1-f_1)+(1-f_2)=\frac{7}{12},\qquad \mathbb{E}[a_2]=1-f_1=\frac{5}{18}.$$
Folglich
$$\mathbb{E}[\operatorname{rank}(\sigma)]=1+2!\cdot\frac{7}{12}+1!\cdot\frac{5}{18}=\frac{22}{9},$$
und damit
$$Q(3)=3!^2\cdot\frac{22}{9}=88.$$
Genau dieser erste nichttriviale Prüfwert bestätigt, dass die geschlossene Herleitung mit der exakten Enumeration übereinstimmt.
Die C++-, Python- und Java-Implementierungen folgen alle derselben Herleitung. Jede rationale Größe wird modulo \(1\,000\,000\,007\) durch einen modularen Inversen dargestellt, was hier zulässig ist, weil \(n=10^6<1\,000\,000\,007\) gilt.
Zuerst erzeugt die Implementierung die modularen Inversen von \(1,2,\dots,n\). Daraus entstehen die Präfixarrays für \(H_m\) und \(H_m^{(2)}\). Zusätzlich werden die Werte der Eulerschen Phi-Funktion bis \(n\) mit einem linearen Sieb berechnet, Präfixsummen von \(\varphi(d)/d^2\) aufgebaut und die Fakultäten \(0!,1!,\dots,n!\) für die Lehmer-Gewichte vorbereitet.
Da im affinen Vergleichsgesetz der Faktor \(1/(n-2)\) auftaucht, werden die Kleinfälle \(n=1\) und \(n=2\) vor der allgemeinen Verzweigung separat behandelt. Kleine Prüfwerte wie \(Q(2)=5\) und \(Q(3)=88\) dienen zusätzlich als Konsistenztest der Herleitung.
Mit diesen Präfixdaten berechnet die Implementierung \(S(n)\) blockweise nach gleichen Werten von \(\left\lfloor n/d\right\rfloor\). Anschließend werden \(\alpha\), \(\beta\), \(\gamma\) und \(\delta\) bestimmt, die Einpunkt-Normierung \(\alpha+(n-1)\beta=1\) kontrolliert, daraus das affine Gesetz \(f_d=A+Bd\) aufgebaut und schließlich für jede Position die quadratische Formel für \(\mathbb{E}[a_j]\) ausgewertet.
Im letzten Schritt wird
$$1+\sum_{j=1}^{n}\mathbb{E}[a_j]\,(n-j)!$$
aufsummiert, um den erwarteten lexikographischen Rang zu erhalten. Danach wird mit \(n!^2\) multipliziert, um \(Q(n)\) zurückzugewinnen. Die C++- und Java-Versionen parallelisieren nur diese letzte gewichtete Summe über Positionsbereiche; die Python-Version führt dieselbe Rechnung seriell aus.
Das Berechnen der Inversen, harmonischen Präfixe, Totienten, Präfixsummen, Fakultäten und der finalen Erwartungswertsumme benötigt jeweils \(O(n)\) Zeit. Die Quotientenblock-Auswertung von \(S(n)\) kostet nach dem Aufbau der Präfixfelder nur \(O(\sqrt n)\), sodass die Gesamtlaufzeit \(O(n)\) bleibt.
Der Speicherverbrauch ist \(O(n)\), weil mehrere Arrays der Länge \(n+1\) gehalten werden. Für \(n=10^6\) ist das praktisch beherrschbar; ein direktes Durchlaufen aller Permutationen oder Potenzen wäre es nicht.
\(n \ge 1\) için
$$Q(n)=\sum_{\pi\in S_n}\sum_{k=1}^{n!}\operatorname{rank}(\pi^k)\pmod{1\,000\,000\,007}$$
tanımlansın. Burada \(\operatorname{rank}(\pi^k)\), \(\pi^k\) permütasyonunun tek satırlı gösterimde tüm \(n!\) permütasyon arasındaki 1 tabanlı leksikografik sırasıdır. Problem 903, \(n=10^6\) için bu değeri ister.
Doğrudan yaklaşım tamamen çıkmazdır: hem \(n!\) permütasyonu hem de her biri için \(n!\) kuvveti dolaşmak gerekir. Gerçek çözüm permütasyon kuvvetlerini hiç üretmez. Bunun yerine rastgele bir permütasyonun rastgele bir kuvvetinin beklenen leksikografik sırasını hesaplar ve sonucu \(n!^2\) ile çarpar.
Ana nesneler, düzgün rastgele seçilmiş bir \(\pi\in S_n\) permütasyonu, düzgün rastgele seçilmiş bir \(i\in\{1,\dots,n!\}\) üssü ve \(\sigma=\pi^i\) rastgele kuvvetidir. \(\sigma\) içindeki ikili karşılaştırma olasılıkları bulununca, toplamın tamamı Lehmer kodu formülünden çıkar.
\(\pi\) permütasyonunu \(S_n\) içinde düzgün rastgele, \(i\) değerini de \(\{1,\dots,n!\}\) aralığında düzgün rastgele seçelim. Tanım:
$$\sigma=\pi^i.$$
Orijinal çift toplamda her \((\pi,i)\) ikilisi tam bir kez yer aldığı için
$$Q(n)=n!^2\,\mathbb{E}[\operatorname{rank}(\sigma)]$$
elde edilir. Böylece devasa bir sayma problemi tek bir beklenti problemine indirgenir.
\(\{1,\dots,n\}\) üzerinde bir \(\tau\) permütasyonu için
$$a_j(\tau)=\#\{k>j:\tau(k)<\tau(j)\}$$
olsun. O zaman 1 tabanlı leksikografik sıra
$$\operatorname{rank}(\tau)=1+\sum_{j=1}^{n} a_j(\tau)\,(n-j)!$$
şeklindedir. Bunu \(\sigma\) için uygularsak
$$\mathbb{E}[\operatorname{rank}(\sigma)]=1+\sum_{j=1}^{n}\mathbb{E}[a_j]\,(n-j)!$$
bulunur. Şimdi sabit \(x<y\) için \(d=y-x\) yazalım ve
$$f_d=\Pr\bigl(\sigma(x)<\sigma(y)\bigr)$$
tanımını yapalım. O halde
$$\mathbb{E}[a_j]=\sum_{d=1}^{n-j}\bigl(1-f_d\bigr),$$
çünkü \(a_j\), daha sonraki konumlar arasında \(\sigma(j)\)'den küçük olanların sayısını verir. Böylece tüm problem \(f_d\) için açık bir formül bulmaya indirgenir.
Uygulamalar, iki farklı nokta \(x\neq y\) için dört temel olasılığı kullanır:
$$\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).$$
Bu dört nicelik düzgün rastgele bir permütasyonun çevrim yapısından doğrudan gelir.
Düzgün rastgele bir \(S_n\) permütasyonunda işaretli bir noktanın çevrim uzunluğu \(\{1,\dots,n\}\) üzerinde düzgündür. Eğer bu uzunluk \(\ell\) ise, \(\pi^i\) bu noktayı ancak \(\ell\mid i\) olduğunda sabit bırakır. Her \(\ell\le n\), \(n!\)'i böldüğü için, düzgün rastgele \(i\in\{1,\dots,n!\}\) altında bunun olasılığı \(1/\ell\)'dir.
Dolayısıyla
$$\alpha=\frac1n\sum_{\ell=1}^{n}\frac1\ell=\frac{H_n}{n},\qquad H_n=\sum_{r=1}^{n}\frac1r.$$
\(x\) sabit değilse kalan \(n-1\) olası görüntü tamamen simetriktir. Buradan
$$\beta=\frac{1-\alpha}{n-1}$$
çıkar.
\(\gamma\) olasılığı yalnızca \(x\) ile \(y\) bir \(2t\) uzunluklu çift çevrim üzerinde tam karşılıklı duruyorsa mümkündür. O halde \(x\)'in çevrim uzunluğu \(2t\) olmalı, \(y\) o çevrimin tek karşı noktası olmalı ve üs \(i\equiv t\pmod{2t}\) sağlamalıdır. Bu da
$$\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)}$$
sonucunu verir.
\(\delta\) için, \(x\) ile \(y\)'nin aynı çevrimde olup olmamasına göre ayrım yaparız. \(x\) ve \(y\), uzunluğu \(\ell\) olan aynı çevrimdeyse, \(x\)'in çevrim uzunluğunu sabitledikten sonra \(y\)'nin de o çevrimde olma olasılığı \((\ell-1)/(n-1)\) ve her ikisinin \(\pi^i\) altında sabit kalma olasılığı \(1/\ell\)'dir. Bu katkı
$$\frac1n\sum_{\ell=1}^{n}\frac{\ell-1}{n-1}\cdot\frac1\ell=\frac{n-H_n}{n(n-1)}$$
olur.
Farklı çevrimler için yararlı bir düzgünlük vardır: \(a,b\ge 1\) ve \(a+b\le n\) koşulunu sağlayan her sıralı \((a,b)\) çifti için, “\(x\)'in çevrimi \(a\), \(y\)'nin çevrimi \(b\) uzunluğunda ve bu iki çevrim farklıdır” olayının olasılığı \(1/(n(n-1))\)'dir. Çünkü \(x\)'in çevrim uzunluğu \(a\) olma olasılığı \(1/n\)'dir; buna koşullu olarak \(y\)'nin bu çevrimin dışında kalma olasılığı \((n-a)/(n-1)\)'dir; kalan \(n-a\) nokta üzerinde de \(y\)'nin çevrim uzunluğu \(\{1,\dots,n-a\}\) üzerinde düzgündür. Bu durumda iki noktanın birlikte sabit kalması tam olarak \(i\), \(\operatorname{lcm}(a,b)\)'nin katı olduğunda gerçekleşir. Dolayısıyla farklı-çevrim katkısı
$$\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)}$$
olur; burada
$$S(n)=\sum_{a=1}^{n-1}\sum_{b=1}^{n-a}\frac{1}{\operatorname{lcm}(a,b)}$$
ve dolayısıyla
$$\delta=\frac{n-H_n+S(n)}{n(n-1)}$$
elde edilir.
Şimdi \(x<y\) sabit olsun ve \(d=y-x\) yazalım. \(\sigma(x)<\sigma(y)\) olayını, sıralı görüntü çifti \((\sigma(x),\sigma(y))\)'ni kategorilere ayırarak sayarız. \((x,y)\) çifti doğrudan \(\delta\) katkısını verir; ters sıra \((y,x)\) ise \(\gamma\) katkısını verir ve eşitsizliği hiç sağlamaz.
Yalnızca bir görüntü sabitse, her somut seçim \((\alpha-\delta)/(n-2)\) olasılığına sahiptir. \(\sigma(x)=x\) ise \(\sigma(y)\), \(y\) hariç \(x\)'ten büyük herhangi bir değer olabilir; bu da \(n-x-1\) uygun seçenek verir. \(\sigma(y)=y\) ise \(\sigma(x)\), \(x\) hariç \(y\)'den küçük herhangi bir değer olabilir; bu da \(y-2\) seçenek verir. Toplam katkı
$$n-x-1+y-2=n+d-3$$
olur.
Yalnızca bir görüntü diğer işaretli değere gidiyorsa, her somut seçim \((\beta-\gamma)/(n-2)\) olasılığına sahiptir. \(\sigma(x)=y\) durumunda \(\sigma(y)\), \(y\)'den büyük olmalıdır; bu \(n-y\) seçenek verir. \(\sigma(y)=x\) durumunda \(\sigma(x)\), \(x\)'ten küçük olmalıdır; bu da \(x-1\) seçenek verir. Dolayısıyla bu kategori
$$n-y+x-1=n-d-1$$
katkısını üretir.
\(\{x,y\}\) kümesinden tamamen kaçınan bütün görüntü çiftleri ise simetri nedeniyle toplam kütlelerinin tam yarısını katkı olarak verir. Bu dört kategoriyi birleştirince
$$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}$$
elde edilir. Böylece somut konumların etkisi yalnızca fark \(d\) üzerinden kalır ve ifade doğrusal olur:
$$f_d=A+Bd,$$
burada
$$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}.$$
\(a_j\), \(\sigma(j)\)'den küçük olan sonraki değerleri saydığı için \(1-f_d\) terimleri toplanır:
$$\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}.$$
Uygulamalar bunun \(m=j-1\) yazılarak elde edilen eşdeğer ikinci dereceden biçimini kullanır:
$$\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}.$$
Son toplamda \((n-j)!\) ile çarpılan nicelik tam olarak budur.
Gerçekten pahalı tek terim \(S(n)\)'dir. Doğrudan çift toplam karesel maliyetlidir, bu yüzden çözüm
$$\frac{1}{\operatorname{lcm}(a,b)}=\frac{\gcd(a,b)}{ab}=\sum_{d\mid a,\ d\mid b}\frac{\varphi(d)}{ab}$$
yazımını kullanır; burada \(\gcd(a,b)=\sum_{d\mid\gcd(a,b)}\varphi(d)\) özdeşliği kullanılır. Toplamların sırası değiştirilince
$$S(n)=\sum_{d=1}^{n}\frac{\varphi(d)}{d^2}\,T\!\left(\left\lfloor\frac{n}{d}\right\rfloor\right)$$
çıkar; burada
$$T(m)=\sum_{u=1}^{m-1}\sum_{v=1}^{m-u}\frac1{uv}.$$
İç toplam
$$\frac1{uv}=\frac{1}{u+v}\left(\frac1u+\frac1v\right)$$
yardımıyla
$$T(m)=H_m^2-H_m^{(2)},\qquad H_m^{(2)}=\sum_{r=1}^{m}\frac1{r^2}$$
harmonik özdeşliğine indirgenir. Son olarak \(\left\lfloor n/d\right\rfloor\) yalnızca \(O(\sqrt n)\) farklı değer aldığı için, \(d\) toplamı eşit bölüm sonuçlarına göre bloklar halinde hesaplanabilir.
\(n=3\) için harmonik değerler
$$H_3=\frac{11}{6},\qquad H_{\lfloor 3/2\rfloor}=1$$
olur; EKOK toplamı da
$$S(3)=\frac1{\operatorname{lcm}(1,1)}+\frac1{\operatorname{lcm}(1,2)}+\frac1{\operatorname{lcm}(2,1)}=1+\frac12+\frac12=2$$
şeklindedir.
Böylece
$$\alpha=\frac{11}{18},\qquad \beta=\frac{7}{36},\qquad \gamma=\frac{1}{12},\qquad \delta=\frac{19}{36}$$
elde edilir.
Bunları doğrusal yasaya koyunca
$$A=\frac34,\qquad B=-\frac1{36},\qquad f_d=\frac34-\frac{d}{36}$$
çıkar. Dolayısıyla ilgili iki mesafe için
$$f_1=\frac{13}{18},\qquad f_2=\frac{25}{36}$$
olur.
Beklenen Lehmer basamakları
$$\mathbb{E}[a_1]=(1-f_1)+(1-f_2)=\frac{7}{12},\qquad \mathbb{E}[a_2]=1-f_1=\frac{5}{18}$$
şeklindedir. Buradan
$$\mathbb{E}[\operatorname{rank}(\sigma)]=1+2!\cdot\frac{7}{12}+1!\cdot\frac{5}{18}=\frac{22}{9}$$
ve son olarak
$$Q(3)=3!^2\cdot\frac{22}{9}=88$$
bulunur. Bu değer, kapalı formülün tam sayımla uyuştuğunu gösteren ilk ciddi kontrol noktasıdır.
C++, Python ve Java uygulamaları aynı türevi uygular. Tüm rasyonel nicelikler modüler tersler kullanılarak \(1\,000\,000\,007\) modunda temsil edilir; bu burada geçerlidir çünkü \(n=10^6<1\,000\,000\,007\).
Önce \(1,2,\dots,n\) için modüler tersler hesaplanır. Bunlardan \(H_m\) ve \(H_m^{(2)}\) önek dizileri oluşturulur. Ayrıca Euler totient değerleri lineer elek ile çıkarılır, \(\varphi(d)/d^2\) önek toplamları hazırlanır ve Lehmer ağırlıkları için \(0!,1!,\dots,n!\) faktöriyelleri hesaplanır.
Doğrusal karşılaştırma yasasında \(1/(n-2)\) çarpanı bulunduğu için, \(n=1\) ve \(n=2\) özel durumları genel dala girilmeden önce ayrı ele alınır. \(Q(2)=5\) ve \(Q(3)=88\) gibi küçük kontrol değerleri de türevin tutarlılığını sınar.
Önek veriler hazır olunca uygulama \(S(n)\)'yi \(\left\lfloor n/d\right\rfloor\) değerlerine göre bloklayarak hesaplar. Sonra \(\alpha\), \(\beta\), \(\gamma\), \(\delta\) bulunur, tek nokta normlaması \(\alpha+(n-1)\beta=1\) kontrol edilir, bunlardan \(f_d=A+Bd\) doğrusal yasası çıkarılır ve her \(\mathbb{E}[a_j]\) için ikinci dereceden ifade değerlendirilir.
Son aşama
$$1+\sum_{j=1}^{n}\mathbb{E}[a_j]\,(n-j)!$$
toplamını alarak beklenen leksikografik sırayı üretir; ardından \(Q(n)\)'i geri kazanmak için sonuç \(n!^2\) ile çarpılır. C++ ve Java sürümleri yalnızca bu son ağırlıklı toplamı konum aralıklarına bölerek paralelleştirir; Python sürümü aynı hesabı seri biçimde yapar.
Tersler, harmonik önekler, totientler, önek toplamlar, faktöriyeller ve son beklenti toplamı \(O(n)\) zamanda hesaplanır. \(S(n)\)'nin bölüm-blok değerlendirmesi, önek diziler hazırlandıktan sonra yalnızca \(O(\sqrt n)\) sürer; dolayısıyla toplam süre yine \(O(n)\) kalır.
Bellek kullanımı \(O(n)\)'dir; çünkü yöntem uzunluğu \(n+1\) olan birkaç dizi tutar. Bu, \(n=10^6\) için uygulanabilirdir; buna karşılık permütasyonları ya da kuvvetleri doğrudan saymak tamamen imkansızdır.
Para \(n \ge 1\), definimos
$$Q(n)=\sum_{\pi\in S_n}\sum_{k=1}^{n!}\operatorname{rank}(\pi^k)\pmod{1\,000\,000\,007},$$
donde \(\operatorname{rank}(\pi^k)\) es el rango lexicográfico basado en 1 de la potencia \(\pi^k\) escrita en notación de una sola fila. El Problema 903 pide el valor para \(n=10^6\).
Un enfoque directo sería imposible: habría que recorrer \(n!\) permutaciones y, para cada una, \(n!\) potencias. La solución real no enumera potencias de permutaciones. En cambio, calcula el rango esperado de una potencia aleatoria de una permutación aleatoria y luego multiplica por \(n!^2\).
Los objetos centrales son una permutación aleatoria uniforme \(\pi\in S_n\), un exponente uniforme \(i\in\{1,\dots,n!\}\) y la potencia aleatoria \(\sigma=\pi^i\). Una vez conocidas las probabilidades de comparación por pares dentro de \(\sigma\), toda la suma sale de la fórmula del código de Lehmer para el rango lexicográfico.
Elegimos \(\pi\) uniformemente en \(S_n\) y \(i\) uniformemente en \(\{1,\dots,n!\}\). Definimos
$$\sigma=\pi^i.$$
Cada par ordenado \((\pi,i)\) aparece exactamente una vez en la suma original, de modo que
$$Q(n)=n!^2\,\mathbb{E}[\operatorname{rank}(\sigma)].$$
Con esto, el problema deja de ser una enumeración gigantesca y pasa a ser un único problema de esperanza.
Para una permutación \(\tau\) sobre \(\{1,\dots,n\}\), definimos
$$a_j(\tau)=\#\{k>j:\tau(k)<\tau(j)\}.$$
Entonces el rango lexicográfico basado en 1 es
$$\operatorname{rank}(\tau)=1+\sum_{j=1}^{n} a_j(\tau)\,(n-j)!.$$
Aplicado a \(\sigma\), esto da
$$\mathbb{E}[\operatorname{rank}(\sigma)]=1+\sum_{j=1}^{n}\mathbb{E}[a_j]\,(n-j)!.$$
Para índices fijos \(x<y\), escribimos \(d=y-x\) y definimos
$$f_d=\Pr\bigl(\sigma(x)<\sigma(y)\bigr).$$
Entonces
$$\mathbb{E}[a_j]=\sum_{d=1}^{n-j}\bigl(1-f_d\bigr),$$
porque \(a_j\) cuenta las posiciones posteriores cuyo valor es menor que \(\sigma(j)\). Por tanto, todo el problema se reduce a describir \(f_d\).
La implementación utiliza cuatro probabilidades básicas para dos puntos distintos \(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).$$
Las cuatro salen directamente de la estructura de ciclos de una permutación uniforme aleatoria.
La longitud del ciclo de un punto marcado en una permutación uniforme de \(S_n\) es uniforme en \(\{1,\dots,n\}\). Si ese punto está en un ciclo de longitud \(\ell\), entonces \(\pi^i\) lo fija exactamente cuando \(\ell\mid i\). Como todo \(\ell\le n\) divide a \(n!\), para un exponente uniforme \(i\in\{1,\dots,n!\}\) esa condición tiene probabilidad \(1/\ell\).
Por lo tanto
$$\alpha=\frac1n\sum_{\ell=1}^{n}\frac1\ell=\frac{H_n}{n},\qquad H_n=\sum_{r=1}^{n}\frac1r.$$
Condicionado a que \(x\) no quede fijo, los otros \(n-1\) posibles valores de la imagen son simétricos, de modo que
$$\beta=\frac{1-\alpha}{n-1}.$$
La probabilidad de intercambio \(\gamma\) solo puede ocurrir cuando \(x\) e \(y\) están enfrentados en un ciclo par de longitud \(2t\). La longitud del ciclo de \(x\) debe ser \(2t\), el punto \(y\) debe ser el punto opuesto único de ese ciclo y el exponente debe satisfacer \(i\equiv t\pmod{2t}\). Así se obtiene
$$\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)}.$$
Para \(\delta\), conviene separar según \(x\) e \(y\) estén o no en el mismo ciclo. Si comparten un ciclo de longitud \(\ell\), entonces, una vez fijada la longitud del ciclo de \(x\), el punto \(y\) está en ese mismo ciclo con probabilidad \((\ell-1)/(n-1)\), y ambos quedan fijos por \(\pi^i\) con probabilidad \(1/\ell\). Eso aporta
$$\frac1n\sum_{\ell=1}^{n}\frac{\ell-1}{n-1}\cdot\frac1\ell=\frac{n-H_n}{n(n-1)}.$$
Si están en ciclos distintos, aparece además una ley uniforme muy útil: para cada par ordenado \((a,b)\) con \(a,b\ge 1\) y \(a+b\le n\), el evento “el ciclo de \(x\) tiene longitud \(a\), el ciclo de \(y\) tiene longitud \(b\) y ambos ciclos son distintos” tiene probabilidad \(1/(n(n-1))\). En efecto, \(x\) tiene ciclo de longitud \(a\) con probabilidad \(1/n\); condicionado a eso, \(y\) queda fuera de ese ciclo con probabilidad \((n-a)/(n-1)\); y sobre los \(n-a\) puntos restantes, la longitud del ciclo de \(y\) es uniforme en \(\{1,\dots,n-a\}\). En ese caso ambos puntos quedan fijos exactamente cuando \(i\) es múltiplo de \(\operatorname{lcm}(a,b)\), de modo que la contribución de ciclos distintos es
$$\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)},$$
donde
$$S(n)=\sum_{a=1}^{n-1}\sum_{b=1}^{n-a}\frac{1}{\operatorname{lcm}(a,b)},$$
y por tanto
$$\delta=\frac{n-H_n+S(n)}{n(n-1)}.$$
Fijemos ahora \(x<y\) y escribamos \(d=y-x\). El evento \(\sigma(x)<\sigma(y)\) se cuenta separando el par ordenado \((\sigma(x),\sigma(y))\) en categorías. El propio par \((x,y)\) aporta \(\delta\), mientras que el intercambio \((y,x)\) aporta \(\gamma\) y nunca satisface la desigualdad.
Cuando exactamente una imagen queda fija, cada elección concreta tiene probabilidad \((\alpha-\delta)/(n-2)\). Si \(\sigma(x)=x\), entonces \(\sigma(y)\) puede ser cualquier valor mayor que \(x\) salvo \(y\), lo que da \(n-x-1\) casos favorables. Si \(\sigma(y)=y\), entonces \(\sigma(x)\) puede ser cualquier valor menor que \(y\) salvo \(x\), lo que da \(y-2\) casos favorables. En total, esta categoría aporta
$$n-x-1+y-2=n+d-3.$$
Cuando exactamente una imagen cae en el otro valor marcado, cada elección concreta tiene probabilidad \((\beta-\gamma)/(n-2)\). Si \(\sigma(x)=y\), entonces \(\sigma(y)\) debe ser mayor que \(y\), y eso da \(n-y\) posibilidades favorables. Si \(\sigma(y)=x\), entonces \(\sigma(x)\) debe ser menor que \(x\), y eso da \(x-1\) posibilidades. Por tanto, esta categoría aporta
$$n-y+x-1=n-d-1.$$
Todos los pares de imágenes que evitan por completo \(\{x,y\}\) contribuyen exactamente la mitad de su masa total por simetría. Combinando las cuatro categorías se llega a
$$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}.$$
Toda la dependencia respecto a las posiciones concretas se reduce al salto \(d\), y la expresión es afín:
$$f_d=A+Bd,$$
con
$$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}.$$
Como \(a_j\) cuenta los valores posteriores menores que \(\sigma(j)\), sumamos \(1-f_d\) sobre todas las distancias posteriores:
$$\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}.$$
La implementación usa la forma cuadrática equivalente obtenida al escribir \(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}.$$
Ésa es exactamente la cantidad que luego se multiplica por \((n-j)!\) en la suma final.
El único término realmente costoso es \(S(n)\). Una doble suma directa sería cuadrática, así que la solución reescribe
$$\frac{1}{\operatorname{lcm}(a,b)}=\frac{\gcd(a,b)}{ab}=\sum_{d\mid a,\ d\mid b}\frac{\varphi(d)}{ab},$$
usando la identidad \(\gcd(a,b)=\sum_{d\mid\gcd(a,b)}\varphi(d)\). Al intercambiar el orden de la suma, se obtiene
$$S(n)=\sum_{d=1}^{n}\frac{\varphi(d)}{d^2}\,T\!\left(\left\lfloor\frac{n}{d}\right\rfloor\right),$$
donde
$$T(m)=\sum_{u=1}^{m-1}\sum_{v=1}^{m-u}\frac1{uv}.$$
La suma interior se simplifica con
$$\frac1{uv}=\frac{1}{u+v}\left(\frac1u+\frac1v\right),$$
lo que conduce a la identidad armónica
$$T(m)=H_m^2-H_m^{(2)},\qquad H_m^{(2)}=\sum_{r=1}^{m}\frac1{r^2}.$$
Además, \(\left\lfloor n/d\right\rfloor\) solo toma \(O(\sqrt n)\) valores distintos, así que la suma sobre \(d\) se puede evaluar por bloques de cocientes iguales después de precomputar los prefijos de \(\varphi(d)/d^2\).
Para \(n=3\), los valores armónicos son
$$H_3=\frac{11}{6},\qquad H_{\lfloor 3/2\rfloor}=1,$$
y la suma con mínimos comunes múltiplos vale
$$S(3)=\frac1{\operatorname{lcm}(1,1)}+\frac1{\operatorname{lcm}(1,2)}+\frac1{\operatorname{lcm}(2,1)}=1+\frac12+\frac12=2.$$
De ahí salen
$$\alpha=\frac{11}{18},\qquad \beta=\frac{7}{36},\qquad \gamma=\frac{1}{12},\qquad \delta=\frac{19}{36}.$$
Al sustituir en la ley afín obtenemos
$$A=\frac34,\qquad B=-\frac1{36},\qquad f_d=\frac34-\frac{d}{36}.$$
Por tanto, los dos saltos relevantes son
$$f_1=\frac{13}{18},\qquad f_2=\frac{25}{36}.$$
Los dígitos esperados del código de Lehmer quedan
$$\mathbb{E}[a_1]=(1-f_1)+(1-f_2)=\frac{7}{12},\qquad \mathbb{E}[a_2]=1-f_1=\frac{5}{18}.$$
En consecuencia,
$$\mathbb{E}[\operatorname{rank}(\sigma)]=1+2!\cdot\frac{7}{12}+1!\cdot\frac{5}{18}=\frac{22}{9},$$
y finalmente
$$Q(3)=3!^2\cdot\frac{22}{9}=88.$$
Éste es el primer valor de control no trivial que confirma que las fórmulas cerradas coinciden con la enumeración exacta.
Las implementaciones en C++, Python y Java siguen exactamente esta derivación. Toda cantidad racional se representa módulo \(1\,000\,000\,007\) con inversos modulares, lo cual es válido porque aquí \(n=10^6<1\,000\,000\,007\).
La implementación construye primero los inversos modulares de \(1,2,\dots,n\). Con ellos forma los arreglos prefijo de \(H_m\) y \(H_m^{(2)}\). También calcula los valores de la función totiente hasta \(n\) con un tamiz lineal, acumula prefijos de \(\varphi(d)/d^2\) y prepara los factoriales \(0!,1!,\dots,n!\) necesarios en la fórmula del rango de Lehmer.
Como la ley afín de comparación contiene el factor \(1/(n-2)\), los casos pequeños \(n=1\) y \(n=2\) se resuelven aparte antes de entrar en la rama general. Valores de control como \(Q(2)=5\) y \(Q(3)=88\) sirven además como prueba de consistencia de la derivación.
Con los prefijos listos, la implementación calcula \(S(n)\) agrupando intervalos con el mismo valor de \(\left\lfloor n/d\right\rfloor\). Después obtiene \(\alpha\), \(\beta\), \(\gamma\) y \(\delta\), comprueba la normalización de un punto \(\alpha+(n-1)\beta=1\), las convierte en la ley afín \(f_d=A+Bd\) y evalúa la fórmula cuadrática de cada dígito esperado \(\mathbb{E}[a_j]\).
La fase final acumula
$$1+\sum_{j=1}^{n}\mathbb{E}[a_j]\,(n-j)!$$
para obtener el rango lexicográfico esperado, y después multiplica por \(n!^2\) para recuperar \(Q(n)\). Las versiones en C++ y Java paralelizan solo esta última suma ponderada sobre rangos de posiciones; la versión en Python realiza la misma aritmética en serie.
Calcular inversos, prefijos armónicos, totientes, prefijos, factoriales y la suma final del rango esperado cuesta \(O(n)\) tiempo. La evaluación de \(S(n)\) por bloques de cociente solo requiere \(O(\sqrt n)\) una vez construidos los prefijos, así que el tiempo total sigue siendo \(O(n)\).
El uso de memoria es \(O(n)\), porque el método mantiene varios arreglos de longitud \(n+1\). Eso es práctico para \(n=10^6\); enumerar permutaciones o potencias directamente no lo es.
对任意 \(n \ge 1\)，定义
$$Q(n)=\sum_{\pi\in S_n}\sum_{k=1}^{n!}\operatorname{rank}(\pi^k)\pmod{1\,000\,000\,007},$$
其中 \(\operatorname{rank}(\pi^k)\) 表示排列 \(\pi^k\) 以单行表示写出后，在全部 \(n!\) 个排列中的 1 开始字典序排名。Problem 903 要求的是 \(n=10^6\) 时的这个值。
如果直接做，必须遍历 \(n!\) 个排列，并且对每个排列再处理 \(n!\) 个幂，这在计算上根本不可能。真正的解法完全不枚举这些幂，而是先求“随机排列的随机幂”的期望字典序排名，再把结果乘上 \(n!^2\)。
核心对象是一组随机变量：从 \(S_n\) 中等概率选出的排列 \(\pi\)，从 \(\{1,\dots,n!\}\) 中等概率选出的指数 \(i\)，以及对应的随机幂 \(\sigma=\pi^i\)。一旦知道 \(\sigma\) 中各个位置对的比较概率，整个总和就能通过 Lehmer 码公式完全恢复出来。
令 \(\pi\) 在 \(S_n\) 上均匀随机，令 \(i\) 在 \(\{1,\dots,n!\}\) 上均匀随机，并定义
$$\sigma=\pi^i.$$
原题中的每一个有序对 \((\pi,i)\) 都恰好对应双重求和中的一个项，因此
$$Q(n)=n!^2\,\mathbb{E}[\operatorname{rank}(\sigma)].$$
这一步是第一层关键化简：原题不再是一个不可承受的枚举问题，而变成了一个单独的期望问题。
对于 \(\{1,\dots,n\}\) 上的任意排列 \(\tau\)，定义
$$a_j(\tau)=\#\{k>j:\tau(k)<\tau(j)\}.$$
那么从 1 开始计数的字典序排名满足
$$\operatorname{rank}(\tau)=1+\sum_{j=1}^{n} a_j(\tau)\,(n-j)!.$$
应用到 \(\sigma\) 上，就得到
$$\mathbb{E}[\operatorname{rank}(\sigma)]=1+\sum_{j=1}^{n}\mathbb{E}[a_j]\,(n-j)!.$$
现在固定 \(x<y\)，记 \(d=y-x\)，并定义
$$f_d=\Pr\bigl(\sigma(x)<\sigma(y)\bigr).$$
由于 \(a_j\) 统计的是后面比 \(\sigma(j)\) 更小的值的个数，所以有
$$\mathbb{E}[a_j]=\sum_{d=1}^{n-j}\bigl(1-f_d\bigr).$$
因此，整道题的实质变成了：如何精确写出 \(f_d\)。
实现里引入了两个不同点 \(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).$$
这四个量都直接来自均匀随机排列的循环结构。
在均匀随机排列中，一个被标记点所在循环的长度在 \(\{1,\dots,n\}\) 上是均匀分布的。如果这个长度是 \(\ell\)，那么 \(\pi^i\) 固定该点，当且仅当 \(\ell\mid i\)。因为所有 \(\ell\le n\) 都整除 \(n!\)，所以对于均匀随机的 \(i\in\{1,\dots,n!\}\)，这个事件的概率正好是 \(1/\ell\)。
于是
$$\alpha=\frac1n\sum_{\ell=1}^{n}\frac1\ell=\frac{H_n}{n},\qquad H_n=\sum_{r=1}^{n}\frac1r.$$
而一旦 \(x\) 没有被固定，它可以落到其余 \(n-1\) 个值中的任何一个。由对称性可知这些目标等概率，所以
$$\beta=\frac{1-\alpha}{n-1}.$$
交换概率 \(\gamma\) 只可能出现在一种情形：\(x\) 和 \(y\) 位于一个长度为 \(2t\) 的偶数循环中，并且恰好彼此相对。也就是说，\(x\) 的循环长度必须是 \(2t\)，\(y\) 必须是那个循环上唯一与 \(x\) 相隔 \(t\) 步的点，同时还要求 \(i\equiv t\pmod{2t}\)。因此
$$\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)}.$$
对于 \(\delta\)，需要先按 \(x\) 和 \(y\) 是否位于同一个循环来分类。如果 \(x\) 和 \(y\) 位于同一个长度为 \(\ell\) 的循环中，那么在固定 \(x\) 的循环长度之后，\(y\) 落在该循环中的概率是 \((\ell-1)/(n-1)\)，而二者同时被 \(\pi^i\) 固定的概率是 \(1/\ell\)。这一部分贡献为
$$\frac1n\sum_{\ell=1}^{n}\frac{\ell-1}{n-1}\cdot\frac1\ell=\frac{n-H_n}{n(n-1)}.$$
如果它们位于不同循环中，还有一个很有用的均匀规律：对任意满足 \(a,b\ge 1\) 且 \(a+b\le n\) 的有序对 \((a,b)\)，事件“\(x\) 所在循环长度为 \(a\)，\(y\) 所在循环长度为 \(b\)，并且这两个循环不同”的概率恒为 \(1/(n(n-1))\)。原因是：\(x\) 的循环长度为 \(a\) 的概率是 \(1/n\)；在此条件下，\(y\) 落在该循环之外的概率是 \((n-a)/(n-1)\)；而在剩下的 \(n-a\) 个点上，\(y\) 的循环长度又在 \(\{1,\dots,n-a\}\) 上均匀分布。在这种情况下，二者同时被固定，当且仅当 \(i\) 是 \(\operatorname{lcm}(a,b)\) 的倍数，所以不同循环的贡献就是
$$\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)},$$
其中
$$S(n)=\sum_{a=1}^{n-1}\sum_{b=1}^{n-a}\frac{1}{\operatorname{lcm}(a,b)},$$
于是
$$\delta=\frac{n-H_n+S(n)}{n(n-1)}.$$
现在固定 \(x<y\)，记 \(d=y-x\)。事件 \(\sigma(x)<\sigma(y)\) 可以按有序像对 \((\sigma(x),\sigma(y))\) 的类别来计数。配对 \((x,y)\) 本身贡献 \(\delta\)，而交换配对 \((y,x)\) 贡献 \(\gamma\)，并且永远不会满足不等式。
当恰好有一个像被固定时，每一种具体选择的概率都是 \((\alpha-\delta)/(n-2)\)。若 \(\sigma(x)=x\)，那么 \(\sigma(y)\) 可以是所有大于 \(x\) 且不等于 \(y\) 的值，共有 \(n-x-1\) 个有利选择。若 \(\sigma(y)=y\)，那么 \(\sigma(x)\) 可以是所有小于 \(y\) 且不等于 \(x\) 的值，共有 \(y-2\) 个有利选择。因此这一类总共贡献
$$n-x-1+y-2=n+d-3.$$
当恰好有一个像落到另一个特殊值上时，每一种具体选择的概率都是 \((\beta-\gamma)/(n-2)\)。若 \(\sigma(x)=y\)，则 \(\sigma(y)\) 必须大于 \(y\)，共有 \(n-y\) 个有利选择。若 \(\sigma(y)=x\)，则 \(\sigma(x)\) 必须小于 \(x\)，共有 \(x-1\) 个有利选择。因此这一类贡献
$$n-y+x-1=n-d-1.$$
至于两个像都避开 \(\{x,y\}\) 的情形，由对称性可知它们恰好贡献其总概率质量的一半。把这四类合并起来，就得到
$$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}.$$
可以看到，具体的 \(x\) 和 \(y\) 已经不重要了，剩下的只是一维参数 \(d\)。因此它必然写成
$$f_d=A+Bd,$$
其中
$$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}.$$
因为 \(a_j\) 统计的是后面有多少项比 \(\sigma(j)\) 小，所以把所有后续距离上的 \(1-f_d\) 累加起来：
$$\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}.$$
实现里实际使用的是等价的二次形式。令 \(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}.$$
这正是最后要乘上 \((n-j)!\) 再求和的那个量。
整个推导里唯一真正昂贵的量是 \(S(n)\)。如果直接二重求和，复杂度会是二次的。因此解法先写出
$$\frac{1}{\operatorname{lcm}(a,b)}=\frac{\gcd(a,b)}{ab}=\sum_{d\mid a,\ d\mid b}\frac{\varphi(d)}{ab},$$
这里用到了 \(\gcd(a,b)=\sum_{d\mid\gcd(a,b)}\varphi(d)\)。交换求和顺序后得到
$$S(n)=\sum_{d=1}^{n}\frac{\varphi(d)}{d^2}\,T\!\left(\left\lfloor\frac{n}{d}\right\rfloor\right),$$
其中
$$T(m)=\sum_{u=1}^{m-1}\sum_{v=1}^{m-u}\frac1{uv}.$$
内层和式可以借助
$$\frac1{uv}=\frac{1}{u+v}\left(\frac1u+\frac1v\right)$$
化简成调和数恒等式
$$T(m)=H_m^2-H_m^{(2)},\qquad H_m^{(2)}=\sum_{r=1}^{m}\frac1{r^2}.$$
接下来再利用 \(\left\lfloor n/d\right\rfloor\) 只有 \(O(\sqrt n)\) 个不同取值这一事实，就可以按整除分块来计算 \(S(n)\)，而不必枚举所有 \((a,b)\) 对。
当 \(n=3\) 时，调和数取值为
$$H_3=\frac{11}{6},\qquad H_{\lfloor 3/2\rfloor}=1,$$
而最小公倍数求和为
$$S(3)=\frac1{\operatorname{lcm}(1,1)}+\frac1{\operatorname{lcm}(1,2)}+\frac1{\operatorname{lcm}(2,1)}=1+\frac12+\frac12=2.$$
因此
$$\alpha=\frac{11}{18},\qquad \beta=\frac{7}{36},\qquad \gamma=\frac{1}{12},\qquad \delta=\frac{19}{36}.$$
把这些值代入仿射规律后得到
$$A=\frac34,\qquad B=-\frac1{36},\qquad f_d=\frac34-\frac{d}{36}.$$
于是仅有的两个相关距离是
$$f_1=\frac{13}{18},\qquad f_2=\frac{25}{36}.$$
期望 Lehmer 码位变成
$$\mathbb{E}[a_1]=(1-f_1)+(1-f_2)=\frac{7}{12},\qquad \mathbb{E}[a_2]=1-f_1=\frac{5}{18}.$$
从而
$$\mathbb{E}[\operatorname{rank}(\sigma)]=1+2!\cdot\frac{7}{12}+1!\cdot\frac{5}{18}=\frac{22}{9},$$
最终得到
$$Q(3)=3!^2\cdot\frac{22}{9}=88.$$
这正是实现用来验证闭式推导的第一个非平凡检查值。
C++、Python 和 Java 三个实现遵循完全相同的数学路线。所有有理量都通过模逆元在 \(1\,000\,000\,007\) 这个素数模数下表示；由于这里 \(n=10^6<1\,000\,000\,007\)，所以这些分母都可逆。
实现首先求出 \(1,2,\dots,n\) 的模逆元，并据此构造 \(H_m\) 与 \(H_m^{(2)}\) 的前缀数组。随后用线性筛计算 \(1\) 到 \(n\) 的欧拉函数值，累加 \(\varphi(d)/d^2\) 的前缀和，同时准备 \(0!,1!,\dots,n!\) 这些阶乘，以便在 Lehmer 排名公式中使用。
因为仿射比较公式里含有 \(1/(n-2)\)，所以 \(n=1\) 和 \(n=2\) 这两个极小情形会在进入通用分支之前单独处理。像 \(Q(2)=5\)、\(Q(3)=88\) 这样的检验值也会用于确认整个推导没有偏差。
在前缀数据就绪后，实现按 \(\left\lfloor n/d\right\rfloor\) 相同的区间分块计算 \(S(n)\)。然后依次求出 \(\alpha\)、\(\beta\)、\(\gamma\)、\(\delta\)，检查单点分布的归一化关系 \(\alpha+(n-1)\beta=1\)，再把它们组合成仿射规律 \(f_d=A+Bd\)，并用上面的二次公式计算每一个期望 Lehmer 码位 \(\mathbb{E}[a_j]\)。
最后一步求和
$$1+\sum_{j=1}^{n}\mathbb{E}[a_j]\,(n-j)!$$
得到期望字典序排名，再乘上 \(n!^2\) 恢复出 \(Q(n)\)。C++ 和 Java 版本只对最后这段按位置加权求和做并行化；Python 版本执行相同的算术，只是串行完成。
逆元、调和前缀、欧拉函数、前缀和、阶乘以及最后的期望排名累加都需要 \(O(n)\) 时间。利用整除分块计算 \(S(n)\) 在前缀数组准备好之后只需要 \(O(\sqrt n)\)，因此整体时间复杂度仍然是 \(O(n)\)。
空间复杂度为 \(O(n)\)，因为算法维护了若干个长度为 \(n+1\) 的数组。对于实际参数 \(n=10^6\)，这完全可行；而直接枚举排列或幂则完全不可行。
Для \(n \ge 1\) определим
$$Q(n)=\sum_{\pi\in S_n}\sum_{k=1}^{n!}\operatorname{rank}(\pi^k)\pmod{1\,000\,000\,007},$$
где \(\operatorname{rank}(\pi^k)\) — это лексикографический ранг с нумерацией от 1 перестановки \(\pi^k\), записанной в однострочном виде. В задаче 903 требуется вычислить это значение при \(n=10^6\).
Прямой подход здесь невозможен: пришлось бы перебирать \(n!\) перестановок и для каждой ещё \(n!\) степеней. Реальное решение вообще не перечисляет степени перестановок. Вместо этого оно вычисляет ожидаемый ранг случайной степени случайной перестановки, а затем умножает результат на \(n!^2\).
Главные объекты — равномерно случайная перестановка \(\pi\in S_n\), равномерно случайный показатель \(i\in\{1,\dots,n!\}\) и случайная степень \(\sigma=\pi^i\). Как только известны попарные вероятности сравнения внутри \(\sigma\), вся сумма восстанавливается через формулу кода Лемера для лексикографического ранга.
Выберем \(\pi\) равновероятно из \(S_n\) и \(i\) равновероятно из \(\{1,\dots,n!\}\). Положим
$$\sigma=\pi^i.$$
Каждая упорядоченная пара \((\pi,i)\) появляется в исходной двойной сумме ровно один раз, поэтому
$$Q(n)=n!^2\,\mathbb{E}[\operatorname{rank}(\sigma)].$$
Тем самым задача перестаёт быть чудовищной по размеру суммой и превращается в одну задачу на ожидание.
Для перестановки \(\tau\) на \(\{1,\dots,n\}\) определим
$$a_j(\tau)=\#\{k>j:\tau(k)<\tau(j)\}.$$
Тогда лексикографический ранг с нумерацией от 1 равен
$$\operatorname{rank}(\tau)=1+\sum_{j=1}^{n} a_j(\tau)\,(n-j)!.$$
Для \(\sigma\) отсюда следует
$$\mathbb{E}[\operatorname{rank}(\sigma)]=1+\sum_{j=1}^{n}\mathbb{E}[a_j]\,(n-j)!.$$
Теперь зафиксируем \(x<y\), обозначим \(d=y-x\) и введём
$$f_d=\Pr\bigl(\sigma(x)<\sigma(y)\bigr).$$
Тогда
$$\mathbb{E}[a_j]=\sum_{d=1}^{n-j}\bigl(1-f_d\bigr),$$
поскольку \(a_j\) считает более поздние позиции, где значение меньше \(\sigma(j)\). Значит, ключевой вопрос — найти явную формулу для \(f_d\).
Для двух различных точек \(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).$$
Все они выводятся напрямую из структуры циклов равномерно случайной перестановки.
Длина цикла отмеченной точки в равномерно случайной перестановке из \(S_n\) равномерна на \(\{1,\dots,n\}\). Если длина цикла равна \(\ell\), то перестановка \(\pi^i\) фиксирует эту точку тогда и только тогда, когда \(\ell\mid i\). Поскольку каждое \(\ell\le n\) делит \(n!\), для равномерного \(i\in\{1,\dots,n!\}\) это происходит с вероятностью \(1/\ell\).
Следовательно,
$$\alpha=\frac1n\sum_{\ell=1}^{n}\frac1\ell=\frac{H_n}{n},\qquad H_n=\sum_{r=1}^{n}\frac1r.$$
Если \(x\) не фиксирован, то все остальные \(n-1\) возможных образов симметричны, поэтому
$$\beta=\frac{1-\alpha}{n-1}.$$
Вероятность обмена \(\gamma\) возникает только тогда, когда \(x\) и \(y\) находятся напротив друг друга в чётном цикле длины \(2t\). Значит, длина цикла точки \(x\) должна быть \(2t\), точка \(y\) должна быть единственной противоположной точкой этого цикла, а показатель должен удовлетворять сравнению \(i\equiv t\pmod{2t}\). Поэтому
$$\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)}.$$
Для \(\delta\) удобно разделить случаи по тому, лежат ли \(x\) и \(y\) в одном цикле. Если \(x\) и \(y\) лежат в одном цикле длины \(\ell\), то после фиксации длины цикла точки \(x\) точка \(y\) попадает в этот же цикл с вероятностью \((\ell-1)/(n-1)\), а обе точки фиксируются перестановкой \(\pi^i\) с вероятностью \(1/\ell\). Это даёт вклад
$$\frac1n\sum_{\ell=1}^{n}\frac{\ell-1}{n-1}\cdot\frac1\ell=\frac{n-H_n}{n(n-1)}.$$
Если же точки находятся в разных циклах, возникает полезный равномерный закон: для каждой упорядоченной пары \((a,b)\) с \(a,b\ge 1\) и \(a+b\le n\) событие «цикл точки \(x\) имеет длину \(a\), цикл точки \(y\) имеет длину \(b\), и это разные циклы» имеет вероятность \(1/(n(n-1))\). Действительно, длина цикла точки \(x\) равна \(a\) с вероятностью \(1/n\); при этом \(y\) лежит вне этого цикла с вероятностью \((n-a)/(n-1)\); а на оставшихся \(n-a\) точках длина цикла точки \(y\) равномерна на \(\{1,\dots,n-a\}\). В таком случае обе точки фиксируются тогда и только тогда, когда \(i\) кратно \(\operatorname{lcm}(a,b)\), так что вклад разных циклов равен
$$\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)},$$
где
$$S(n)=\sum_{a=1}^{n-1}\sum_{b=1}^{n-a}\frac{1}{\operatorname{lcm}(a,b)},$$
и в итоге
$$\delta=\frac{n-H_n+S(n)}{n(n-1)}.$$
Теперь зафиксируем \(x<y\) и положим \(d=y-x\). Событие \(\sigma(x)<\sigma(y)\) удобно считать, разбивая упорядоченную пару образов \((\sigma(x),\sigma(y))\) на категории. Пара \((x,y)\) сама по себе даёт вклад \(\delta\), а перестановка местами \((y,x)\) даёт вклад \(\gamma\) и никогда не удовлетворяет неравенству.
Когда ровно один образ фиксирован, каждая конкретная возможность имеет вероятность \((\alpha-\delta)/(n-2)\). Если \(\sigma(x)=x\), то \(\sigma(y)\) может быть любым значением больше \(x\), кроме \(y\), то есть имеется \(n-x-1\) благоприятных вариантов. Если \(\sigma(y)=y\), то \(\sigma(x)\) может быть любым значением меньше \(y\), кроме \(x\), то есть \(y-2\) вариантов. Итого эта категория даёт
$$n-x-1+y-2=n+d-3.$$
Когда ровно один образ попадает в другое отмеченное значение, каждая конкретная возможность имеет вероятность \((\beta-\gamma)/(n-2)\). Если \(\sigma(x)=y\), то \(\sigma(y)\) должно быть больше \(y\), что даёт \(n-y\) благоприятных вариантов. Если \(\sigma(y)=x\), то \(\sigma(x)\) должно быть меньше \(x\), что даёт \(x-1\) вариантов. Следовательно, эта категория даёт
$$n-y+x-1=n-d-1.$$
Все пары образов, полностью избегающие множества \(\{x,y\}\), по симметрии дают ровно половину своей суммарной массы. Объединяя эти четыре категории, получаем формулу
$$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}.$$
После этого зависимость от конкретных позиций исчезает: остаётся только расстояние \(d\), и выражение становится аффинным:
$$f_d=A+Bd,$$
где
$$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}.$$
Поскольку \(a_j\) считает, сколько более поздних значений меньше \(\sigma(j)\), нужно просуммировать \(1-f_d\) по всем допустимым расстояниям:
$$\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}.$$
В реализации используется эквивалентная квадратичная форма после замены \(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}.$$
Именно эта величина затем умножается на \((n-j)!\) в финальной сумме.
Единственный действительно дорогой объект — это \(S(n)\). Прямая двойная сумма имела бы квадратичную сложность, поэтому решение переписывает
$$\frac{1}{\operatorname{lcm}(a,b)}=\frac{\gcd(a,b)}{ab}=\sum_{d\mid a,\ d\mid b}\frac{\varphi(d)}{ab},$$
используя тождество \(\gcd(a,b)=\sum_{d\mid\gcd(a,b)}\varphi(d)\). После перестановки суммирования получается
$$S(n)=\sum_{d=1}^{n}\frac{\varphi(d)}{d^2}\,T\!\left(\left\lfloor\frac{n}{d}\right\rfloor\right),$$
где
$$T(m)=\sum_{u=1}^{m-1}\sum_{v=1}^{m-u}\frac1{uv}.$$
Внутренняя сумма упрощается через равенство
$$\frac1{uv}=\frac{1}{u+v}\left(\frac1u+\frac1v\right),$$
что приводит к гармонической формуле
$$T(m)=H_m^2-H_m^{(2)},\qquad H_m^{(2)}=\sum_{r=1}^{m}\frac1{r^2}.$$
Далее используется факт, что \(\left\lfloor n/d\right\rfloor\) принимает только \(O(\sqrt n)\) различных значений. Поэтому сумму по \(d\) можно вычислять блоками одинаковых частных, если заранее подготовлены префиксные суммы для \(\varphi(d)/d^2\).
При \(n=3\) гармонические величины равны
$$H_3=\frac{11}{6},\qquad H_{\lfloor 3/2\rfloor}=1,$$
а сумма по наименьшему общему кратному равна
$$S(3)=\frac1{\operatorname{lcm}(1,1)}+\frac1{\operatorname{lcm}(1,2)}+\frac1{\operatorname{lcm}(2,1)}=1+\frac12+\frac12=2.$$
Отсюда получаем
$$\alpha=\frac{11}{18},\qquad \beta=\frac{7}{36},\qquad \gamma=\frac{1}{12},\qquad \delta=\frac{19}{36}.$$
Подстановка в аффинный закон даёт
$$A=\frac34,\qquad B=-\frac1{36},\qquad f_d=\frac34-\frac{d}{36}.$$
Значит, для двух возможных расстояний имеем
$$f_1=\frac{13}{18},\qquad f_2=\frac{25}{36}.$$
Ожидаемые цифры кода Лемера равны
$$\mathbb{E}[a_1]=(1-f_1)+(1-f_2)=\frac{7}{12},\qquad \mathbb{E}[a_2]=1-f_1=\frac{5}{18}.$$
Следовательно,
$$\mathbb{E}[\operatorname{rank}(\sigma)]=1+2!\cdot\frac{7}{12}+1!\cdot\frac{5}{18}=\frac{22}{9},$$
и потому
$$Q(3)=3!^2\cdot\frac{22}{9}=88.$$
Это первое нетривиальное контрольное значение, подтверждающее совпадение закрытых формул с точным перебором.
Реализации на C++, Python и Java следуют одной и той же математической схеме. Все рациональные величины представляются по модулю простого числа \(1\,000\,000\,007\) с помощью модульных обратных; это корректно, потому что здесь \(n=10^6<1\,000\,000\,007\).
Сначала вычисляются модульные обратные для \(1,2,\dots,n\). По ним строятся префиксные массивы для \(H_m\) и \(H_m^{(2)}\). Затем линейным решетом вычисляются значения функции Эйлера до \(n\), накапливаются префиксные суммы \(\varphi(d)/d^2\), а также подготавливаются факториалы \(0!,1!,\dots,n!\), используемые в формуле ранга Лемера.
Поскольку в аффинном законе сравнения присутствует множитель \(1/(n-2)\), маленькие случаи \(n=1\) и \(n=2\) разбираются отдельно ещё до основной ветви. Контрольные значения вроде \(Q(2)=5\) и \(Q(3)=88\) дополнительно проверяют согласованность вывода.
После подготовки префиксных данных реализация вычисляет \(S(n)\) по блокам одинаковых значений \(\left\lfloor n/d\right\rfloor\). Затем находятся \(\alpha\), \(\beta\), \(\gamma\) и \(\delta\), проверяется нормировка для одной точки \(\alpha+(n-1)\beta=1\), из них строится аффинный закон \(f_d=A+Bd\), и по квадратичной формуле вычисляется каждое ожидаемое значение \(\mathbb{E}[a_j]\).
На последнем этапе суммируется
$$1+\sum_{j=1}^{n}\mathbb{E}[a_j]\,(n-j)!$$
для получения ожидаемого лексикографического ранга, после чего результат умножается на \(n!^2\), что и даёт \(Q(n)\). В версиях на C++ и Java распараллеливается только эта последняя взвешенная сумма по диапазонам позиций; версия на Python выполняет ту же арифметику последовательно.
Построение обратных, гармонических префиксов, значений \(\varphi\), префиксных сумм, факториалов и финального накопления ожидаемого ранга требует \(O(n)\) времени. Вычисление \(S(n)\) по блокам частного после этого занимает лишь \(O(\sqrt n)\), так что общая асимптотика остаётся \(O(n)\).
По памяти метод использует \(O(n)\), поскольку хранит несколько массивов длины \(n+1\). Для \(n=10^6\) это практично; прямой перебор перестановок или степеней абсолютно нереален.
لكل \(n \ge 1\)، نعرّف
$$Q(n)=\sum_{\pi\in S_n}\sum_{k=1}^{n!}\operatorname{rank}(\pi^k)\pmod{1\,000\,000\,007},$$
حيث تمثل \(\operatorname{rank}(\pi^k)\) الرتبة المعجمية المفهرسة من 1 للتبديل \(\pi^k\) عندما يُكتب بصيغة السطر الواحد. في المسألة 903 المطلوب هو قيمة هذه الكمية عند \(n=10^6\).
الحل المباشر غير ممكن عمليًا، لأنه يتطلب المرور على \(n!\) تبديلًا وعلى \(n!\) قوة لكل واحد منها. الحل الفعلي لا يولد هذه القوى أصلًا، بل يحسب الرتبة المعجمية المتوقعة لقوة عشوائية من تبديل عشوائي، ثم يضرب الناتج في \(n!^2\).
العناصر الأساسية هنا هي تبديل عشوائي منتظم \(\pi\in S_n\)، وأُسّ عشوائي منتظم \(i\in\{1,\dots,n!\}\)، والقوة العشوائية \(\sigma=\pi^i\). وما إن نعرف احتمالات المقارنة الثنائية داخل \(\sigma\)، حتى يصبح مجموع المسألة كله نتيجة مباشرة من صيغة ترميز ليهمر للرتبة المعجمية.
نختار \(\pi\) بتوزيع منتظم من \(S_n\)، ونختار \(i\) بتوزيع منتظم من \(\{1,\dots,n!\}\)، ثم نضع
$$\sigma=\pi^i.$$
كل زوج مرتب \((\pi,i)\) يظهر مرة واحدة بالضبط في المجموع الأصلي، ولذلك
$$Q(n)=n!^2\,\mathbb{E}[\operatorname{rank}(\sigma)].$$
وهكذا تتحول المسألة من تعداد ضخم جدًا إلى مسألة قيمة متوقعة واحدة.
لكل تبديل \(\tau\) على \(\{1,\dots,n\}\)، نعرّف
$$a_j(\tau)=\#\{k>j:\tau(k)<\tau(j)\}.$$
عندئذ تكون الرتبة المعجمية ذات الفهرسة من 1 هي
$$\operatorname{rank}(\tau)=1+\sum_{j=1}^{n} a_j(\tau)\,(n-j)!.$$
وبتطبيق ذلك على \(\sigma\) نحصل على
$$\mathbb{E}[\operatorname{rank}(\sigma)]=1+\sum_{j=1}^{n}\mathbb{E}[a_j]\,(n-j)!.$$
والآن نثبت \(x<y\)، ونكتب \(d=y-x\)، ثم نعرّف
$$f_d=\Pr\bigl(\sigma(x)<\sigma(y)\bigr).$$
ومن ثم
$$\mathbb{E}[a_j]=\sum_{d=1}^{n-j}\bigl(1-f_d\bigr),$$
لأن \(a_j\) يعدّ المواضع اللاحقة التي تحمل قيمًا أصغر من \(\sigma(j)\). إذن لبّ المسألة هو إيجاد صيغة دقيقة لـ \(f_d\).
تستخدم التطبيقات أربع احتمالات أساسية لنقطتين مختلفتين \(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).$$
هذه المقادير الأربعة تنشأ مباشرة من بنية الدورات في تبديل عشوائي منتظم.
طول الدورة التي تحتوي نقطة معلّمة في تبديل عشوائي منتظم من \(S_n\) موزع بانتظام على \(\{1,\dots,n\}\). وإذا كان طول الدورة هو \(\ell\)، فإن \(\pi^i\) يثبت هذه النقطة إذا وفقط إذا كان \(\ell\mid i\). وبما أن كل \(\ell\le n\) يقسم \(n!\)، فإن هذا يحدث مع أُسّ منتظم \(i\in\{1,\dots,n!\}\) باحتمال \(1/\ell\).
لذلك
$$\alpha=\frac1n\sum_{\ell=1}^{n}\frac1\ell=\frac{H_n}{n},\qquad H_n=\sum_{r=1}^{n}\frac1r.$$
وإذا لم تُثبَّت \(x\)، فإن الصور الأخرى الممكنة وعددها \(n-1\) تكون متناظرة تمامًا، ومن ثم
$$\beta=\frac{1-\alpha}{n-1}.$$
احتمال التبادل \(\gamma\) لا يمكن أن يحدث إلا عندما تقع \(x\) و\(y\) في موضعين متقابلين داخل دورة زوجية طولها \(2t\). وهذا يفرض أن يكون طول دورة \(x\) مساويًا لـ \(2t\)، وأن تكون \(y\) هي النقطة المقابلة الوحيدة، وأن يحقق الأس \(i\equiv t\pmod{2t}\). لذلك
$$\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)}.$$
أما \(\delta\) فمن الأفضل تقسيمها بحسب ما إذا كانت \(x\) و\(y\) في الدورة نفسها أم لا. إذا كانت \(x\) و\(y\) في الدورة نفسها وطولها \(\ell\)، فإن احتمال وجود \(y\) في تلك الدورة بعد تثبيت طول دورة \(x\) هو \((\ell-1)/(n-1)\)، واحتمال أن يثبتهما \(\pi^i\) معًا هو \(1/\ell\). وهذا يعطي المساهمة
$$\frac1n\sum_{\ell=1}^{n}\frac{\ell-1}{n-1}\cdot\frac1\ell=\frac{n-H_n}{n(n-1)}.$$
أما إذا كانت النقطتان في دورتين مختلفتين، فهناك انتظام مفيد جدًا: لكل زوج مرتب \((a,b)\) يحقق \(a,b\ge 1\) و\(a+b\le n\)، يكون الحدث «طول دورة \(x\) هو \(a\)، وطول دورة \(y\) هو \(b\)، والدورتان مختلفتان» ذا احتمال \(1/(n(n-1))\). والسبب أن طول دورة \(x\) يساوي \(a\) باحتمال \(1/n\)، ثم بافتراض ذلك تكون \(y\) خارج تلك الدورة باحتمال \((n-a)/(n-1)\)، وعلى النقاط \(n-a\) المتبقية يكون طول دورة \(y\) موزعًا بانتظام على \(\{1,\dots,n-a\}\). وفي هذه الحالة تثبت النقطتان معًا بالضبط عندما يكون \(i\) من مضاعفات \(\operatorname{lcm}(a,b)\)، ولذلك فإن مساهمة الدورتين المختلفتين هي
$$\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)},$$
حيث
$$S(n)=\sum_{a=1}^{n-1}\sum_{b=1}^{n-a}\frac{1}{\operatorname{lcm}(a,b)},$$
ومن ثم
$$\delta=\frac{n-H_n+S(n)}{n(n-1)}.$$
الآن نثبت \(x<y\) ونكتب \(d=y-x\). ويمكن عدّ الحدث \(\sigma(x)<\sigma(y)\) بعد تقسيم الزوج المرتب \((\sigma(x),\sigma(y))\) إلى فئات. فالزوج \((x,y)\) نفسه يعطي المساهمة \(\delta\)، في حين أن الزوج المتبادل \((y,x)\) يعطي \(\gamma\) ولا يحقق المتباينة أبدًا.
عندما تكون صورة واحدة فقط مثبتة، يكون احتمال كل اختيار ملموس مساويًا لـ \((\alpha-\delta)/(n-2)\). فإذا كان \(\sigma(x)=x\)، جاز أن تكون \(\sigma(y)\) أي قيمة أكبر من \(x\) ما عدا \(y\)، وهذا يعطي \(n-x-1\) اختيارًا مناسبًا. وإذا كان \(\sigma(y)=y\)، جاز أن تكون \(\sigma(x)\) أي قيمة أصغر من \(y\) ما عدا \(x\)، وهذا يعطي \(y-2\) اختيارًا مناسبًا. وبذلك تسهم هذه الفئة بعدد
$$n-x-1+y-2=n+d-3.$$
وعندما تقع صورة واحدة فقط على القيمة الخاصة الأخرى، يكون احتمال كل اختيار ملموس مساويًا لـ \((\beta-\gamma)/(n-2)\). فإذا كان \(\sigma(x)=y\)، وجب أن تكون \(\sigma(y)\) أكبر من \(y\)، وهذا يعطي \(n-y\) اختيارات مناسبة. وإذا كان \(\sigma(y)=x\)، وجب أن تكون \(\sigma(x)\) أصغر من \(x\)، وهذا يعطي \(x-1\) اختيارًا. ومن ثم تسهم هذه الفئة بعدد
$$n-y+x-1=n-d-1.$$
أما جميع الأزواج التي تتجنب المجموعة \(\{x,y\}\) بالكامل، فتعطي بالتماثل نصف كتلتها الاحتمالية بالضبط. وبجمع الفئات الأربع نحصل على
$$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}.$$
وهكذا يختفي اعتماد المسألة على الموضعين أنفسهما، ولا يبقى سوى الفاصل \(d\)، ولذلك يكتب القانون على صورة
$$f_d=A+Bd,$$
حيث
$$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}.$$
بما أن \(a_j\) يعدّ القيم اللاحقة الأصغر من \(\sigma(j)\)، فإننا نجمع \(1-f_d\) على جميع المسافات اللاحقة:
$$\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}.$$
وتستخدم التطبيقات الصورة التربيعية المكافئة بعد كتابة \(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}.$$
وهذه هي بالضبط الكمية التي تُضرب في \((n-j)!\) في التجميع النهائي.
الحد الوحيد المكلف فعليًا هو \(S(n)\). فالمجموع المزدوج المباشر سيكون تربيعيًا، ولذلك يُعاد كتابة
$$\frac{1}{\operatorname{lcm}(a,b)}=\frac{\gcd(a,b)}{ab}=\sum_{d\mid a,\ d\mid b}\frac{\varphi(d)}{ab},$$
باستخدام الهوية \(\gcd(a,b)=\sum_{d\mid\gcd(a,b)}\varphi(d)\). وبعد تبديل ترتيب الجمع نحصل على
$$S(n)=\sum_{d=1}^{n}\frac{\varphi(d)}{d^2}\,T\!\left(\left\lfloor\frac{n}{d}\right\rfloor\right),$$
حيث
$$T(m)=\sum_{u=1}^{m-1}\sum_{v=1}^{m-u}\frac1{uv}.$$
ويجري تبسيط هذا المجموع الداخلي عبر العلاقة
$$\frac1{uv}=\frac{1}{u+v}\left(\frac1u+\frac1v\right),$$
فنصل إلى الهوية التوافقية
$$T(m)=H_m^2-H_m^{(2)},\qquad H_m^{(2)}=\sum_{r=1}^{m}\frac1{r^2}.$$
وبما أن \(\left\lfloor n/d\right\rfloor\) لا يأخذ إلا \(O(\sqrt n)\) قيمة مختلفة، يمكن حساب مجموع \(d\) على شكل كتل ذات خارج قسمة ثابت بمجرد تجهيز المجاميع البادئة لـ \(\varphi(d)/d^2\).
عندما \(n=3\)، تكون القيم التوافقية
$$H_3=\frac{11}{6},\qquad H_{\lfloor 3/2\rfloor}=1,$$
ويكون مجموع المضاعف المشترك الأصغر
$$S(3)=\frac1{\operatorname{lcm}(1,1)}+\frac1{\operatorname{lcm}(1,2)}+\frac1{\operatorname{lcm}(2,1)}=1+\frac12+\frac12=2.$$
ومن ثم نحصل على
$$\alpha=\frac{11}{18},\qquad \beta=\frac{7}{36},\qquad \gamma=\frac{1}{12},\qquad \delta=\frac{19}{36}.$$
وبالتعويض في القانون الخطي-الثابت نجد
$$A=\frac34,\qquad B=-\frac1{36},\qquad f_d=\frac34-\frac{d}{36}.$$
إذن المسافتان الممكنتان تعطيان
$$f_1=\frac{13}{18},\qquad f_2=\frac{25}{36}.$$
وتصبح أرقام ليهمر المتوقعة
$$\mathbb{E}[a_1]=(1-f_1)+(1-f_2)=\frac{7}{12},\qquad \mathbb{E}[a_2]=1-f_1=\frac{5}{18}.$$
ولذلك
$$\mathbb{E}[\operatorname{rank}(\sigma)]=1+2!\cdot\frac{7}{12}+1!\cdot\frac{5}{18}=\frac{22}{9},$$
وفي النهاية
$$Q(3)=3!^2\cdot\frac{22}{9}=88.$$
وهذه أول قيمة غير بديهية تؤكد أن الصيغ المغلقة متفقة مع العدّ الدقيق.
تتبع تطبيقات C++ وPython وJava الاشتقاق نفسه تمامًا. وكل كمية كسرية تُمثَّل باستعمال مقلوبات معيارية بترديد العدد الأولي \(1\,000\,000\,007\)، وهذا صحيح هنا لأن \(n=10^6<1\,000\,000\,007\).
يبدأ التنفيذ بحساب المقلوبات المعيارية للأعداد \(1,2,\dots,n\). ومن هذه القيم تُبنى المصفوفات البادئة لـ \(H_m\) و\(H_m^{(2)}\). وبعد ذلك تُحسب قيم دالة أويلر حتى \(n\) بغربال خطي، وتُراكم بادئات \(\varphi(d)/d^2\)، كما تُحضَّر القيم \(0!,1!,\dots,n!\) اللازمة في صيغة رتبة ليهمر.
ولأن قانون المقارنة الخطي-الثابت يحتوي على العامل \(1/(n-2)\)، فإن الحالتين الصغيرتين \(n=1\) و\(n=2\) تُعالجان على نحو منفصل قبل البدء بالفرع العام. كما تُستخدم قيم فحص صغيرة مثل \(Q(2)=5\) و\(Q(3)=88\) للتأكد من اتساق الاشتقاق.
بعد تجهيز البيانات البادئة، يحسب التنفيذ \(S(n)\) على كتل متساوية في قيمة \(\left\lfloor n/d\right\rfloor\). ثم يستخرج \(\alpha\) و\(\beta\) و\(\gamma\) و\(\delta\)، ويتحقق من تطبيع التوزيع أحادي النقطة \(\alpha+(n-1)\beta=1\)، ثم يحولها إلى القانون \(f_d=A+Bd\)، ثم يقيّم الصيغة التربيعية الخاصة بكل \(\mathbb{E}[a_j]\).
في المرحلة الأخيرة يُجمع
$$1+\sum_{j=1}^{n}\mathbb{E}[a_j]\,(n-j)!$$
للحصول على الرتبة المعجمية المتوقعة، ثم يُضرب الناتج في \(n!^2\) لاستعادة \(Q(n)\). وتوازي نسختا C++ وJava فقط هذا الجمع الأخير الموزون على مجالات من المواضع، بينما تنفذ نسخة Python الحساب نفسه على نحو تسلسلي.
حساب المقلوبات، والمجاميع التوافقية البادئة، وقيم \(\varphi\)، والبادئات، والعوامل، والتجميع النهائي للرتبة المتوقعة، كلها تحتاج إلى \(O(n)\) زمنًا. أما تقييم \(S(n)\) بكتل خارج القسمة فيحتاج فقط إلى \(O(\sqrt n)\) بعد تجهيز هذه البادئات، ولذلك يبقى الزمن الكلي \(O(n)\).
أما الذاكرة فهي \(O(n)\)، لأن الطريقة تحتفظ بعدة مصفوفات طول كل منها \(n+1\). وهذا عملي عند \(n=10^6\)، في حين أن تعداد التبديلات أو القوى مباشرة غير عملي إطلاقًا.