For each triple \((a,b,c)\) with \(0 \le a,b,c \le 87\), define \(X=a^3\), \(Y=b^3\), and \(Z=c^3\). The quantity \(N(X,Y,Z)\) counts words made from \(X\) copies of \(i\), \(Y\) copies of \(j\), and \(Z\) copies of \(k\) whose product in the quaternion group \(Q_8\) lands on the required central sign. The final task is
$$A=\sum_{a=0}^{87}\sum_{b=0}^{87}\sum_{c=0}^{87} N(a^3,b^3,c^3)\pmod{888888883}.$$
The outer loop over \(88^3\) cube triples is not the hard part. The real issue is to evaluate \(N(X,Y,Z)\) without enumerating the raw multinomial family of words, whose size is
$$\binom{X+Y+Z}{X,Y,Z}=\frac{(X+Y+Z)!}{X!\,Y!\,Z!}.$$
The implementations solve that counting problem exactly by turning quaternion multiplication into a statement about inversion parity in multiset permutations.
Write \(S=X+Y+Z\). Every word with those letter counts can be viewed as a permutation of the multiset containing \(X\) letters \(i\), \(Y\) letters \(j\), and \(Z\) letters \(k\). The key is to compare such a word with the canonical ordered word \(i^Xj^Yk^Z\).
Fix the alphabet order \(i \lt j \lt k\). For a word \(w\), let \(\operatorname{inv}(w)\) be the number of pairs that are out of that order, equivalently the total number of occurrences of \(j\) before \(i\), plus \(k\) before \(i\), plus \(k\) before \(j\).
Each adjacent swap of two distinct generators flips the sign, because
$$ji=-ij,\qquad kj=-jk,\qquad ik=-ki.$$
Therefore reordering any word into \(i^Xj^Yk^Z\) gives
$$P(w)=(-1)^{\operatorname{inv}(w)}\,i^Xj^Yk^Z.$$
This is the central invariant behind the whole solution: the word value depends only on the fixed ordered product \(i^Xj^Yk^Z\) and on the inversion parity of the word.
The quaternion group has central elements only \(\pm1\). So we first ask when \(i^Xj^Yk^Z\) is already central.
If the parities of \(X\), \(Y\), and \(Z\) are not all the same, then one of the basis directions survives, and the ordered product is one of \(\pm i\), \(\pm j\), or \(\pm k\). In that case no word can contribute to \(N(X,Y,Z)\), so
$$N(X,Y,Z)=0.$$
If \(X\), \(Y\), and \(Z\) are all even, say \(X=2x\), \(Y=2y\), \(Z=2z\), then
$$i^Xj^Yk^Z=(i^2)^x(j^2)^y(k^2)^z=(-1)^{x+y+z}=(-1)^{S/2}.$$
If they are all odd, say \(X=2x+1\), \(Y=2y+1\), \(Z=2z+1\), then
$$i^Xj^Yk^Z=(-1)^{x+y+z}ijk=(-1)^{x+y+z+1}=(-1)^{\lfloor S/2\rfloor}.$$
So whenever the three parities agree, the ordered product simplifies to the compact formula
$$i^Xj^Yk^Z=(-1)^{\lfloor S/2\rfloor}.$$
The condition encoded by the implementations is that the quaternion word must land on the central sign \((-1)^S\). Combining that with the previous identity gives
$$(-1)^{\operatorname{inv}(w)}(-1)^{\lfloor S/2\rfloor}=(-1)^S.$$
Equivalently,
$$\operatorname{inv}(w)\equiv \left\lceil\frac{S}{2}\right\rceil \pmod 2.$$
So \(N(X,Y,Z)\) is not an arbitrary group-theoretic count any more. It is exactly the number of multiset permutations with a prescribed inversion parity, provided the parity filter from the previous subsection is satisfied.
Let
$$T=\binom{S}{X,Y,Z},$$
and let \(E\) and \(O\) denote the numbers of words with even and odd inversion counts. Then
$$E+O=T.$$
To separate the two halves, introduce the inversion generating function
$$\sum_w q^{\operatorname{inv}(w)}=\binom{S}{X,Y,Z}_q,$$
the \(q\)-multinomial over the multiset words \(w\). Evaluating at \(q=-1\) gives
$$E-O=\binom{S}{X,Y,Z}_{q=-1}.$$
For this problem, the only cases that survive the central-sign filter are especially simple:
$$E-O=0 \quad \text{when } X,Y,Z \text{ are all odd},$$
and if \(X,Y,Z\) are all even, with
$$H=\binom{S/2}{X/2,Y/2,Z/2},$$
then
$$E-O=H.$$
That half-scale multinomial \(H\) is exactly the correction term used by the implementations.
Because the required inversion parity is \(\lceil S/2\rceil \bmod 2\), the final count is
$$\boxed{ N(X,Y,Z)= \begin{cases} 0, & \text{if } X,Y,Z \text{ do not have the same parity},\\[6pt] \dfrac{T}{2}, & \text{if } X,Y,Z \text{ are all odd},\\[10pt] \dfrac{T+H}{2}, & \text{if } X,Y,Z \text{ are all even and } S/2 \text{ is even},\\[10pt] \dfrac{T-H}{2}, & \text{if } X,Y,Z \text{ are all even and } S/2 \text{ is odd}. \end{cases}}$$
This is precisely the branch structure shared by the C++, Python, and Java implementations.
The small exact values used in the implementations are good sanity checks for the derivation.
For \((X,Y,Z)=(2,2,2)\), we have \(S=6\),
$$T=\binom{6}{2,2,2}=90,\qquad H=\binom{3}{1,1,1}=6.$$
Since \(S/2=3\) is odd, the admissible words are the odd-inversion ones, so
$$N(2,2,2)=\frac{90-6}{2}=42.$$
For \((X,Y,Z)=(8,8,8)\), we have \(S=24\),
$$T=\binom{24}{8,8,8}=9465511770,\qquad H=\binom{12}{4,4,4}=34650.$$
Now \(S/2=12\) is even, so the admissible words are the even-inversion ones, giving
$$N(8,8,8)=\frac{9465511770+34650}{2}=4732773210.$$
Those are the exact checkpoints verified before the main modular sweep.
Once the closed form for \(N(X,Y,Z)\) is known, the Project Euler quantity is just
$$A=\sum_{a=0}^{87}\sum_{b=0}^{87}\sum_{c=0}^{87} N(a^3,b^3,c^3)\pmod{888888883}.$$
Because cubing preserves parity, a triple contributes only when \(a\), \(b\), and \(c\) are all even or all odd. Every surviving term is then obtained by a constant-time multinomial calculation.
The implementations first build the list of cubes \(0^3,1^3,\dots,87^3\) and precompute factorials and inverse factorials up to
$$3\cdot 87^3=1975509.$$
Since the modulus \(M=888888883\) is prime, every inverse factorial is obtained with Fermat's little theorem:
$$a^{-1}\equiv a^{M-2}\pmod M.$$
This makes every multinomial evaluation \(O(1)\) after preprocessing:
$$\binom{n}{a,b,c}\equiv n!\,(a!)^{-1}(b!)^{-1}(c!)^{-1}\pmod M.$$
For one triple \((X,Y,Z)\), the implementation first applies the parity filter. If the parities do not match, the contribution is zero immediately. Otherwise it computes the full multinomial \(T\), and in the all-even case it also computes the half-scale multinomial \(H\). Dividing by 2 modulo \(M\) is done by multiplying with \(2^{-1}=(M+1)/2\).
Before the main accumulation, the implementations also confirm the exact integer values \(N(2,2,2)=42\) and \(N(8,8,8)=4732773210\) using wide exact arithmetic. After that validation, they run a straightforward serial triple loop over all \(88^3\) cube triples and add the modular contribution of each one.
The factorial and inverse-factorial tables dominate memory, and they are sized up to \(3\cdot 87^3\). That preprocessing phase costs \(O(87^3)\) time and \(O(87^3)\) memory.
The final sweep visits exactly \(88^3\) triples. Each visit performs only constant-time arithmetic once the tables exist, so that phase is \(O(88^3)\). Overall, the algorithm runs in
$$O(87^3+88^3)$$
time and uses \(O(87^3)\) additional memory. The exact checkpoint calculations are tiny compared with those two main phases.
Für jedes Tripel \((a,b,c)\) mit \(0 \le a,b,c \le 87\) setzen wir \(X=a^3\), \(Y=b^3\) und \(Z=c^3\). Die Größe \(N(X,Y,Z)\) zählt Wörter aus \(X\) Kopien von \(i\), \(Y\) Kopien von \(j\) und \(Z\) Kopien von \(k\), deren Produkt in der Quaternionengruppe \(Q_8\) auf das geforderte zentrale Vorzeichen fällt. Gesucht ist also
$$A=\sum_{a=0}^{87}\sum_{b=0}^{87}\sum_{c=0}^{87} N(a^3,b^3,c^3)\pmod{888888883}.$$
Der äußere Durchlauf über \(88^3\) Kubiktripel ist nicht das eigentliche Problem. Entscheidend ist eine geschlossene Formel für \(N(X,Y,Z)\), denn die rohe Anzahl aller Wörter ist der Multinomialkoeffizient
$$\binom{X+Y+Z}{X,Y,Z}=\frac{(X+Y+Z)!}{X!\,Y!\,Z!},$$
und explizites Auflisten dieser Wörter ist unmöglich. Die Implementierungen umgehen das vollständig, indem sie Quaternionenmultiplikation auf eine Aussage über die Inversionsparität von Multimengen-Permutationen zurückführen.
Schreibe \(S=X+Y+Z\). Jedes Wort mit diesen Häufigkeiten kann man als Permutation der Multimenge mit \(X\) Buchstaben \(i\), \(Y\) Buchstaben \(j\) und \(Z\) Buchstaben \(k\) auffassen. Der Schlüssel ist der Vergleich mit dem kanonisch geordneten Wort \(i^Xj^Yk^Z\).
Fixiere die Alphabetordnung \(i \lt j \lt k\). Für ein Wort \(w\) bezeichne \(\operatorname{inv}(w)\) die Anzahl der Paare, die zu dieser Ordnung invertiert sind; also die Gesamtzahl der Vorkommen von \(j\) vor \(i\), von \(k\) vor \(i\) und von \(k\) vor \(j\).
Jeder benachbarte Tausch zweier verschiedener Generatoren ändert das Vorzeichen, denn
$$ji=-ij,\qquad kj=-jk,\qquad ik=-ki.$$
Ordnet man ein beliebiges Wort zu \(i^Xj^Yk^Z\) um, so erhält man daher
$$P(w)=(-1)^{\operatorname{inv}(w)}\,i^Xj^Yk^Z.$$
Das ist die zentrale Invariante der gesamten Lösung: Der Gruppenwert des Wortes hängt nur vom festen geordneten Produkt \(i^Xj^Yk^Z\) und von der Parität der Inversionszahl ab.
Die Quaternionengruppe besitzt als zentrale Elemente nur \(\pm1\). Deshalb fragt man zuerst, wann \(i^Xj^Yk^Z\) überhaupt zentral sein kann.
Wenn die Paritäten von \(X\), \(Y\) und \(Z\) nicht alle gleich sind, bleibt eine der Richtungen \(i\), \(j\) oder \(k\) übrig, und das geordnete Produkt ist eines von \(\pm i\), \(\pm j\) oder \(\pm k\). Dann kann kein Wort zu \(N(X,Y,Z)\) beitragen, also
$$N(X,Y,Z)=0.$$
Sind \(X\), \(Y\) und \(Z\) alle gerade, also \(X=2x\), \(Y=2y\), \(Z=2z\), dann gilt
$$i^Xj^Yk^Z=(i^2)^x(j^2)^y(k^2)^z=(-1)^{x+y+z}=(-1)^{S/2}.$$
Sind sie alle ungerade, also \(X=2x+1\), \(Y=2y+1\), \(Z=2z+1\), dann folgt
$$i^Xj^Yk^Z=(-1)^{x+y+z}ijk=(-1)^{x+y+z+1}=(-1)^{\lfloor S/2\rfloor}.$$
Damit erhält man in allen zulässigen Fällen die kompakte Formel
$$i^Xj^Yk^Z=(-1)^{\lfloor S/2\rfloor}.$$
Die von den Implementierungen kodierte Bedingung lautet, dass das Wort auf das zentrale Vorzeichen \((-1)^S\) fallen muss. Zusammen mit der vorigen Identität ergibt das
$$(-1)^{\operatorname{inv}(w)}(-1)^{\lfloor S/2\rfloor}=(-1)^S.$$
Äquivalent dazu ist
$$\operatorname{inv}(w)\equiv \left\lceil\frac{S}{2}\right\rceil \pmod 2.$$
Damit ist \(N(X,Y,Z)\) kein unbestimmter gruppentheoretischer Ausdruck mehr, sondern genau die Anzahl der Multimengen-Permutationen mit vorgeschriebener Inversionsparität, vorausgesetzt der Paritätsfilter aus dem vorherigen Abschnitt ist erfüllt.
Setze
$$T=\binom{S}{X,Y,Z},$$
und bezeichne mit \(E\) bzw. \(O\) die Anzahl der Wörter mit gerader bzw. ungerader Inversionszahl. Dann gilt
$$E+O=T.$$
Zur Trennung der beiden Teile betrachtet man die Inversions-Erzeugende
$$\sum_w q^{\operatorname{inv}(w)}=\binom{S}{X,Y,Z}_q,$$
also den \(q\)-Multinomialkoeffizienten über alle Wörter \(w\). Setzt man \(q=-1\), so erhält man
$$E-O=\binom{S}{X,Y,Z}_{q=-1}.$$
Für dieses Problem bleiben nach dem Zentralitätsfilter nur zwei Fälle übrig:
$$E-O=0 \quad \text{falls } X,Y,Z \text{ alle ungerade sind},$$
und im vollständig geraden Fall mit
$$H=\binom{S/2}{X/2,Y/2,Z/2}$$
gilt
$$E-O=H.$$
Genau dieser halbierte Multinomialkoeffizient \(H\) ist der Korrekturterm der Implementierungen.
Da die verlangte Inversionsparität gleich \(\lceil S/2\rceil \bmod 2\) ist, lautet die Endformel
$$\boxed{ N(X,Y,Z)= \begin{cases} 0, & \text{wenn } X,Y,Z \text{ nicht dieselbe Parität haben},\\[6pt] \dfrac{T}{2}, & \text{wenn } X,Y,Z \text{ alle ungerade sind},\\[10pt] \dfrac{T+H}{2}, & \text{wenn } X,Y,Z \text{ alle gerade sind und } S/2 \text{ gerade ist},\\[10pt] \dfrac{T-H}{2}, & \text{wenn } X,Y,Z \text{ alle gerade sind und } S/2 \text{ ungerade ist}. \end{cases}}$$
Genau diese Fallunterscheidung wird in C++, Python und Java implementiert.
Die kleinen exakten Werte aus den Implementierungen sind ideale Prüfsteine für die Herleitung.
Für \((X,Y,Z)=(2,2,2)\) ist \(S=6\), also
$$T=\binom{6}{2,2,2}=90,\qquad H=\binom{3}{1,1,1}=6.$$
Weil \(S/2=3\) ungerade ist, sind die zulässigen Wörter genau die mit ungerader Inversionszahl. Also
$$N(2,2,2)=\frac{90-6}{2}=42.$$
Für \((X,Y,Z)=(8,8,8)\) gilt \(S=24\),
$$T=\binom{24}{8,8,8}=9465511770,\qquad H=\binom{12}{4,4,4}=34650.$$
Nun ist \(S/2=12\) gerade, also zählen die Wörter mit gerader Inversionszahl, und damit
$$N(8,8,8)=\frac{9465511770+34650}{2}=4732773210.$$
Diese beiden exakten Werte werden vor der eigentlichen Modulo-Summe überprüft.
Nachdem die geschlossene Formel für \(N(X,Y,Z)\) feststeht, ist die gesuchte Euler-Größe einfach
$$A=\sum_{a=0}^{87}\sum_{b=0}^{87}\sum_{c=0}^{87} N(a^3,b^3,c^3)\pmod{888888883}.$$
Weil Potenzieren mit Exponent 3 die Parität erhält, kann ein Tripel nur dann beitragen, wenn \(a\), \(b\) und \(c\) alle gerade oder alle ungerade sind. Jeder verbleibende Term ist dann eine \(O(1)\)-Multinomialauswertung.
Die Implementierungen bauen zuerst die Liste der Kuben \(0^3,1^3,\dots,87^3\) auf und berechnen Fakultäten sowie inverse Fakultäten bis
$$3\cdot 87^3=1975509.$$
Da der Modul \(M=888888883\) prim ist, lassen sich inverse Fakultäten mit dem kleinen Satz von Fermat berechnen:
$$a^{-1}\equiv a^{M-2}\pmod M.$$
Danach kostet jede Multinomialauswertung nur noch \(O(1)\):
$$\binom{n}{a,b,c}\equiv n!\,(a!)^{-1}(b!)^{-1}(c!)^{-1}\pmod M.$$
Für ein einzelnes Tripel \((X,Y,Z)\) prüft die Implementierung zuerst die Paritäten. Stimmen sie nicht überein, ist der Beitrag sofort null. Andernfalls berechnet sie den vollen Multinomialkoeffizienten \(T\) und im vollständig geraden Fall zusätzlich den halb skalierten Koeffizienten \(H\). Die Division durch 2 im Modulo-Ring wird als Multiplikation mit \(2^{-1}=(M+1)/2\) ausgeführt.
Vor der Hauptsumme bestätigen die Implementierungen außerdem die exakten Werte \(N(2,2,2)=42\) und \(N(8,8,8)=4732773210\) mit exakter Ganzzahlarithmetik. Danach läuft eine direkte serielle Dreifachschleife über alle \(88^3\) Kubiktripel und addiert die modulare Version jedes Beitrags.
Die Fakultäts- und Invers-Fakultäts-Tabellen dominieren den Speicherbedarf; ihre Größe reicht bis \(3\cdot 87^3\). Diese Vorverarbeitung kostet \(O(87^3)\) Zeit und \(O(87^3)\) Speicher.
Die Endsumme besucht genau \(88^3\) Tripel. Jeder Besuch braucht nach der Vorverarbeitung nur konstante Zeit, also kostet diese Phase \(O(88^3)\). Insgesamt erhält man
$$O(87^3+88^3)$$
Zeit und \(O(87^3)\) zusätzlichen Speicher. Die exakten Kontrollrechnungen sind gegenüber diesen beiden Hauptphasen vernachlässigbar.
Her \(0 \le a,b,c \le 87\) üçlüsü için \(X=a^3\), \(Y=b^3\), \(Z=c^3\) tanımlanıyor. \(N(X,Y,Z)\), kuaterniyon grubu \(Q_8\) içinde çarpımı istenen merkezi işarete giden, \(X\) tane \(i\), \(Y\) tane \(j\) ve \(Z\) tane \(k\) harfinden oluşan kelimelerin sayısıdır. Aranan toplam ise
$$A=\sum_{a=0}^{87}\sum_{b=0}^{87}\sum_{c=0}^{87} N(a^3,b^3,c^3)\pmod{888888883}.$$
Zor olan kısım \(88^3\) üçlü üzerinde dolaşmak değildir. Asıl mesele,
$$\binom{X+Y+Z}{X,Y,Z}=\frac{(X+Y+Z)!}{X!\,Y!\,Z!}$$
kadar büyük bir kelime ailesini tek tek üretmeden \(N(X,Y,Z)\) için kapalı form bulmaktır. Uygulamalar bunu, kuaterniyon çarpımını çoklu-küme permütasyonlarının terslik paritesine indirgeyerek yapar.
\(S=X+Y+Z\) olsun. Bu harf sayılarıyla kurulan her kelime, \(X\) adet \(i\), \(Y\) adet \(j\) ve \(Z\) adet \(k\) içeren çoklu kümenin bir permütasyonu olarak düşünülebilir. Ana fikir, bu kelimeyi kanonik sıralı kelime \(i^Xj^Yk^Z\) ile karşılaştırmaktır.
Alfabe sırasını \(i \lt j \lt k\) olarak sabitleyelim. Bir kelime \(w\) için \(\operatorname{inv}(w)\), bu sıraya ters düşen çiftlerin sayısı olsun; yani \(j\)'nin \(i\)'den önce gelmeleri, \(k\)'nin \(i\)'den önce gelmeleri ve \(k\)'nin \(j\)'den önce gelmelerinin toplamı.
Farklı iki üretecin yan yana yer değiştirmesi işareti ters çevirir; çünkü
$$ji=-ij,\qquad kj=-jk,\qquad ik=-ki.$$
Bu yüzden her kelime kanonik sıraya getirildiğinde
$$P(w)=(-1)^{\operatorname{inv}(w)}\,i^Xj^Yk^Z$$
elde edilir. Bütün çözümün temel değişmezi budur: kelimenin grup değeri yalnızca sabit sıralı çarpıma ve terslik sayısının paritesine bağlıdır.
Kuaterniyon grubunun merkezinde yalnızca \(\pm1\) vardır. Dolayısıyla önce \(i^Xj^Yk^Z\) ifadesinin ne zaman merkezi olabileceğini belirlemek gerekir.
\(X\), \(Y\), \(Z\) pariteleri aynı değilse geriye \(i\), \(j\) veya \(k\) yönlerinden biri kalır; sıralı çarpım bu durumda \(\pm i\), \(\pm j\) ya da \(\pm k\) olur. O zaman hiçbir kelime \(N(X,Y,Z)\)'ye katkı veremez ve
$$N(X,Y,Z)=0$$
olur.
\(X\), \(Y\), \(Z\) hep çiftse, yani \(X=2x\), \(Y=2y\), \(Z=2z\) ise
$$i^Xj^Yk^Z=(i^2)^x(j^2)^y(k^2)^z=(-1)^{x+y+z}=(-1)^{S/2}.$$
Hepsi tekse, yani \(X=2x+1\), \(Y=2y+1\), \(Z=2z+1\) ise
$$i^Xj^Yk^Z=(-1)^{x+y+z}ijk=(-1)^{x+y+z+1}=(-1)^{\lfloor S/2\rfloor}.$$
Demek ki üç parite aynı olduğunda sıralı çarpım tek bir formülle yazılabilir:
$$i^Xj^Yk^Z=(-1)^{\lfloor S/2\rfloor}.$$
Uygulamaların kodladığı koşul, kelimenin merkezi işaret olarak \((-1)^S\) değerine gitmesidir. Önceki kimlikle birleştirince
$$(-1)^{\operatorname{inv}(w)}(-1)^{\lfloor S/2\rfloor}=(-1)^S$$
elde edilir. Bu da
$$\operatorname{inv}(w)\equiv \left\lceil\frac{S}{2}\right\rceil \pmod 2$$
demektir. Böylece \(N(X,Y,Z)\), belirsiz bir grup sayımı olmaktan çıkar; uygun parite filtresi geçtiği sürece, doğrudan belirli bir terslik paritesine sahip çoklu-küme permütasyonlarının sayısı olur.
Şimdi
$$T=\binom{S}{X,Y,Z}$$
olsun; \(E\) çift terslikli, \(O\) tek terslikli kelime sayısı olsun. O zaman
$$E+O=T.$$
Bu iki parçayı ayırmak için terslik üreteç polinomunu kullanırız:
$$\sum_w q^{\operatorname{inv}(w)}=\binom{S}{X,Y,Z}_q.$$
\(q=-1\) koyulunca
$$E-O=\binom{S}{X,Y,Z}_{q=-1}$$
elde edilir. Bu problemde merkez filtresini geçen iki durum vardır:
$$E-O=0 \quad \text{eğer } X,Y,Z \text{ hep tekse},$$
ve \(X\), \(Y\), \(Z\) hep çift olduğunda
$$H=\binom{S/2}{X/2,Y/2,Z/2}$$
tanımıyla
$$E-O=H$$
olur. Kodun kullandığı yarım ölçekli düzeltme terimi tam olarak bu \(H\)'dir.
Gerekli terslik paritesi \(\lceil S/2\rceil \bmod 2\) olduğundan son formül şöyledir:
$$\boxed{ N(X,Y,Z)= \begin{cases} 0, & \text{eğer } X,Y,Z \text{ aynı pariteye sahip değilse},\\[6pt] \dfrac{T}{2}, & \text{eğer } X,Y,Z \text{ hep tekse},\\[10pt] \dfrac{T+H}{2}, & \text{eğer } X,Y,Z \text{ hep çift ve } S/2 \text{ çiftse},\\[10pt] \dfrac{T-H}{2}, & \text{eğer } X,Y,Z \text{ hep çift ve } S/2 \text{ tekse}. \end{cases}}$$
C++, Python ve Java uygulamalarındaki dal yapısı aynen budur.
Uygulamaların doğruladığı küçük tam sayılı örnekler, türetimin doğru olduğunu hızlıca gösterir.
\((X,Y,Z)=(2,2,2)\) için \(S=6\) ve
$$T=\binom{6}{2,2,2}=90,\qquad H=\binom{3}{1,1,1}=6.$$
\(S/2=3\) tek olduğu için geçerli kelimeler tek terslikli olanlardır. Dolayısıyla
$$N(2,2,2)=\frac{90-6}{2}=42.$$
\((X,Y,Z)=(8,8,8)\) için ise \(S=24\),
$$T=\binom{24}{8,8,8}=9465511770,\qquad H=\binom{12}{4,4,4}=34650.$$
Bu kez \(S/2=12\) çift olduğu için çift terslikli kelimeler sayılır ve
$$N(8,8,8)=\frac{9465511770+34650}{2}=4732773210$$
elde edilir. Ana modüler toplamdan önce tam olarak bu iki değer kontrol edilir.
\(N(X,Y,Z)\) için kapalı form bilindikten sonra Euler toplamı doğrudan
$$A=\sum_{a=0}^{87}\sum_{b=0}^{87}\sum_{c=0}^{87} N(a^3,b^3,c^3)\pmod{888888883}$$
haline gelir. Küp alma işlemi pariteyi koruduğu için yalnızca \(a\), \(b\), \(c\)'nin ya hep çift ya da hep tek olduğu üçlüler katkı verir. Böylece hayatta kalan her terim \(O(1)\) maliyetli bir multinom hesabına indirgenir.
Uygulamalar önce \(0^3,1^3,\dots,87^3\) küp listesini kurar ve faktöriyel ile ters faktöriyel tablolarını
$$3\cdot 87^3=1975509$$
değerine kadar hazırlar. Modül \(M=888888883\) asal olduğu için ters alma Fermat'ın küçük teoremiyle yapılır:
$$a^{-1}\equiv a^{M-2}\pmod M.$$
Bundan sonra her multinom hesabı \(O(1)\) olur:
$$\binom{n}{a,b,c}\equiv n!\,(a!)^{-1}(b!)^{-1}(c!)^{-1}\pmod M.$$
Bir \((X,Y,Z)\) üçlüsü için uygulama önce parite filtresini uygular. Pariteler uyuşmuyorsa katkı anında sıfırdır. Aksi halde tam multinom \(T\) hesaplanır; tümü çift durumunda ayrıca yarım ölçekli multinom \(H\) de hesaplanır. Modüler ortamda 2'ye bölme, \(2^{-1}=(M+1)/2\) ile çarpma olarak yapılır.
Ana toplama geçmeden önce, uygulamalar tam sayı aritmetiğiyle \(N(2,2,2)=42\) ve \(N(8,8,8)=4732773210\) değerlerini doğrular. Ardından bütün \(88^3\) küp üçlüsü üzerinde basit bir seri üçlü döngü çalıştırılır ve her katkının modüler değeri toplama eklenir.
Faktöriyel ve ters faktöriyel tabloları bellek kullanımının ana kaynağıdır; boyutları \(3\cdot 87^3\)'e kadar çıkar. Bu ön işlem aşaması \(O(87^3)\) zaman ve \(O(87^3)\) bellek gerektirir.
Son tarama tam olarak \(88^3\) üçlü ziyaret eder. Tablolar hazır olduktan sonra her ziyaret sabit zamanlıdır; dolayısıyla bu aşama \(O(88^3)\)'tür. Toplam karmaşıklık
$$O(87^3+88^3)$$
zaman ve \(O(87^3)\) ek bellek olur. Küçük tam sayı doğrulama adımları bu iki ana fazın yanında ihmal edilebilir.
Para cada triple \((a,b,c)\) con \(0 \le a,b,c \le 87\), se define \(X=a^3\), \(Y=b^3\) y \(Z=c^3\). La cantidad \(N(X,Y,Z)\) cuenta las palabras formadas por \(X\) copias de \(i\), \(Y\) copias de \(j\) y \(Z\) copias de \(k\) cuyo producto en el grupo cuaterniónico \(Q_8\) cae en el signo central exigido. El valor buscado es
$$A=\sum_{a=0}^{87}\sum_{b=0}^{87}\sum_{c=0}^{87} N(a^3,b^3,c^3)\pmod{888888883}.$$
El bucle exterior sobre \(88^3\) ternas cúbicas no es el verdadero obstáculo. El punto difícil es evaluar \(N(X,Y,Z)\) sin enumerar la enorme familia de palabras posible, cuyo tamaño bruto es
$$\binom{X+Y+Z}{X,Y,Z}=\frac{(X+Y+Z)!}{X!\,Y!\,Z!}.$$
Las implementaciones resuelven ese conteo convirtiendo el producto de cuaterniones en una condición sobre la paridad del número de inversiones de una permutación con repetición.
Escribamos \(S=X+Y+Z\). Toda palabra con esas cantidades de letras puede verse como una permutación del multiconjunto que contiene \(X\) letras \(i\), \(Y\) letras \(j\) y \(Z\) letras \(k\). La idea clave es compararla con la palabra ordenada canónica \(i^Xj^Yk^Z\).
Fijemos el orden alfabético \(i \lt j \lt k\). Para una palabra \(w\), definimos \(\operatorname{inv}(w)\) como el número de pares que violan ese orden; es decir, cuántas veces aparece \(j\) antes que \(i\), más cuántas veces aparece \(k\) antes que \(i\), más cuántas veces aparece \(k\) antes que \(j\).
Cada intercambio adyacente entre dos generadores distintos cambia el signo, porque
$$ji=-ij,\qquad kj=-jk,\qquad ik=-ki.$$
Por tanto, al reordenar cualquier palabra hasta \(i^Xj^Yk^Z\) se obtiene
$$P(w)=(-1)^{\operatorname{inv}(w)}\,i^Xj^Yk^Z.$$
Ese es el invariante fundamental de toda la solución: el valor del producto depende solo del producto ordenado fijo \(i^Xj^Yk^Z\) y de la paridad del número de inversiones.
En el grupo cuaterniónico, los únicos elementos centrales son \(\pm1\). Por eso lo primero es decidir cuándo \(i^Xj^Yk^Z\) puede ser central.
Si las paridades de \(X\), \(Y\) y \(Z\) no coinciden todas, queda una dirección \(i\), \(j\) o \(k\) sin cancelar, y el producto ordenado es uno de \(\pm i\), \(\pm j\) o \(\pm k\). En ese caso ninguna palabra contribuye a \(N(X,Y,Z)\), y por tanto
$$N(X,Y,Z)=0.$$
Si \(X\), \(Y\) y \(Z\) son todos pares, digamos \(X=2x\), \(Y=2y\), \(Z=2z\), entonces
$$i^Xj^Yk^Z=(i^2)^x(j^2)^y(k^2)^z=(-1)^{x+y+z}=(-1)^{S/2}.$$
Si son todos impares, es decir \(X=2x+1\), \(Y=2y+1\), \(Z=2z+1\), entonces
$$i^Xj^Yk^Z=(-1)^{x+y+z}ijk=(-1)^{x+y+z+1}=(-1)^{\lfloor S/2\rfloor}.$$
Así, siempre que las tres paridades coincidan, el producto ordenado se reduce a
$$i^Xj^Yk^Z=(-1)^{\lfloor S/2\rfloor}.$$
La condición implementada en el código es que la palabra termine en el signo central \((-1)^S\). Combinando eso con la identidad anterior se obtiene
$$(-1)^{\operatorname{inv}(w)}(-1)^{\lfloor S/2\rfloor}=(-1)^S.$$
De manera equivalente,
$$\operatorname{inv}(w)\equiv \left\lceil\frac{S}{2}\right\rceil \pmod 2.$$
Por lo tanto, \(N(X,Y,Z)\) deja de ser un conteo grupal abstracto y se convierte exactamente en el número de permutaciones del multiconjunto con una paridad de inversiones prescrita, siempre que antes se haya superado el filtro de paridad central.
Sea
$$T=\binom{S}{X,Y,Z},$$
y llamemos \(E\) y \(O\) a los números de palabras con número de inversiones par e impar, respectivamente. Entonces
$$E+O=T.$$
Para distinguir ambas mitades se usa la función generadora de inversiones
$$\sum_w q^{\operatorname{inv}(w)}=\binom{S}{X,Y,Z}_q,$$
el \(q\)-multinomial de todas las palabras \(w\). Al evaluar en \(q=-1\), aparece la diferencia
$$E-O=\binom{S}{X,Y,Z}_{q=-1}.$$
En este problema, los casos que sobreviven al filtro central son especialmente simples:
$$E-O=0 \quad \text{si } X,Y,Z \text{ son todos impares},$$
y si \(X,Y,Z\) son todos pares, con
$$H=\binom{S/2}{X/2,Y/2,Z/2},$$
se cumple
$$E-O=H.$$
Ese multinomial a media escala \(H\) es exactamente el término corrector usado por las implementaciones.
Como la paridad de inversiones exigida es \(\lceil S/2\rceil \bmod 2\), el conteo final queda
$$\boxed{ N(X,Y,Z)= \begin{cases} 0, & \text{si } X,Y,Z \text{ no tienen la misma paridad},\\[6pt] \dfrac{T}{2}, & \text{si } X,Y,Z \text{ son todos impares},\\[10pt] \dfrac{T+H}{2}, & \text{si } X,Y,Z \text{ son todos pares y } S/2 \text{ es par},\\[10pt] \dfrac{T-H}{2}, & \text{si } X,Y,Z \text{ son todos pares y } S/2 \text{ es impar}. \end{cases}}$$
Esa es exactamente la estructura de casos que comparten las versiones en C++, Python y Java.
Los valores pequeños comprobados por las implementaciones sirven como validación directa de la derivación.
Para \((X,Y,Z)=(2,2,2)\), se tiene \(S=6\),
$$T=\binom{6}{2,2,2}=90,\qquad H=\binom{3}{1,1,1}=6.$$
Como \(S/2=3\) es impar, cuentan las palabras con inversión impar, y por tanto
$$N(2,2,2)=\frac{90-6}{2}=42.$$
Para \((X,Y,Z)=(8,8,8)\), se obtiene \(S=24\),
$$T=\binom{24}{8,8,8}=9465511770,\qquad H=\binom{12}{4,4,4}=34650.$$
Ahora \(S/2=12\) es par, así que cuentan las palabras con inversión par, lo que da
$$N(8,8,8)=\frac{9465511770+34650}{2}=4732773210.$$
Esos son exactamente los puntos de control verificados antes del barrido modular principal.
Una vez conocida la fórmula cerrada de \(N(X,Y,Z)\), la cantidad de Project Euler es simplemente
$$A=\sum_{a=0}^{87}\sum_{b=0}^{87}\sum_{c=0}^{87} N(a^3,b^3,c^3)\pmod{888888883}.$$
Como elevar al cubo conserva la paridad, solo contribuyen los triples en los que \(a\), \(b\) y \(c\) son todos pares o todos impares. Cada término superviviente se obtiene entonces con una sola evaluación multinomial de costo constante.
Las implementaciones construyen primero la lista de cubos \(0^3,1^3,\dots,87^3\) y precalculan factoriales e inversos factoriales hasta
$$3\cdot 87^3=1975509.$$
Como el módulo \(M=888888883\) es primo, los inversos se obtienen mediante el pequeño teorema de Fermat:
$$a^{-1}\equiv a^{M-2}\pmod M.$$
Así, cada multinomial se evalúa en \(O(1)\) una vez hechas las tablas:
$$\binom{n}{a,b,c}\equiv n!\,(a!)^{-1}(b!)^{-1}(c!)^{-1}\pmod M.$$
Para un triple \((X,Y,Z)\), la implementación aplica primero el filtro de paridad. Si las paridades no coinciden, la contribución es cero de inmediato. En caso contrario calcula el multinomial completo \(T\), y en el caso totalmente par también el multinomial a media escala \(H\). Dividir por 2 en módulo \(M\) equivale a multiplicar por \(2^{-1}=(M+1)/2\).
Antes de acumular la suma principal, las implementaciones validan con aritmética exacta los valores \(N(2,2,2)=42\) y \(N(8,8,8)=4732773210\). Después ejecutan un triple bucle serial sobre las \(88^3\) ternas cúbicas y añaden la contribución modular de cada una.
Las tablas de factoriales e inversos factoriales dominan la memoria, y llegan hasta \(3\cdot 87^3\). Esa fase de precálculo cuesta \(O(87^3)\) tiempo y \(O(87^3)\) memoria.
El barrido final visita exactamente \(88^3\) triples. Cada visita requiere solo tiempo constante tras la precalculación, de modo que esa fase es \(O(88^3)\). En conjunto, el algoritmo usa
$$O(87^3+88^3)$$
tiempo y \(O(87^3)\) memoria adicional. Las comprobaciones exactas son diminutas comparadas con esas dos fases principales.
对每个满足 \(0 \le a,b,c \le 87\) 的三元组 \((a,b,c)\)，定义 \(X=a^3\)、\(Y=b^3\)、\(Z=c^3\)。记 \(N(X,Y,Z)\) 为这样一种字符串数量：它由 \(X\) 个 \(i\)、\(Y\) 个 \(j\)、\(Z\) 个 \(k\) 组成，并且在四元数群 \(Q_8\) 中的乘积恰好落到题目要求的中心符号上。最终要计算的是
$$A=\sum_{a=0}^{87}\sum_{b=0}^{87}\sum_{c=0}^{87} N(a^3,b^3,c^3)\pmod{888888883}.$$
真正的难点不是外层的 \(88^3\) 次循环，而是如何在不枚举全部字符串的情况下求出 \(N(X,Y,Z)\)。因为原始总数已经是
$$\binom{X+Y+Z}{X,Y,Z}=\frac{(X+Y+Z)!}{X!\,Y!\,Z!},$$
这个规模很快就会大到无法直接处理。三个实现共同采用的办法，是把四元数乘法问题转写成“多重集排列的逆序数奇偶性”问题。
记 \(S=X+Y+Z\)。所有满足这些字母计数的字符串，都可以看成是多重集 \(\{i^X,j^Y,k^Z\}\) 的排列。核心思路是：把任意字符串同规范顺序串 \(i^Xj^Yk^Z\) 进行比较。
固定字母顺序 \(i \lt j \lt k\)。对一个字符串 \(w\)，定义 \(\operatorname{inv}(w)\) 为相对这个顺序的逆序对个数。具体地说，就是“\(j\) 出现在 \(i\) 前面”的次数，加上“\(k\) 出现在 \(i\) 前面”的次数，再加上“\(k\) 出现在 \(j\) 前面”的次数。
每交换一次相邻且不同的生成元，符号都会翻转，因为
$$ji=-ij,\qquad kj=-jk,\qquad ik=-ki.$$
因此把任意字符串整理成 \(i^Xj^Yk^Z\) 之后，乘积必然满足
$$P(w)=(-1)^{\operatorname{inv}(w)}\,i^Xj^Yk^Z.$$
这就是整个解法的基本不变量：字符串的群值，只取决于两个因素，一是固定的有序乘积 \(i^Xj^Yk^Z\)，二是逆序数的奇偶性。
四元数群 \(Q_8\) 的中心只有 \(\pm1\)。所以第一步要先判断 \(i^Xj^Yk^Z\) 什么时候可能是中心元素。
如果 \(X\)、\(Y\)、\(Z\) 的奇偶性不全相同，那么最终一定会剩下某个方向 \(i\)、\(j\) 或 \(k\)，于是有序乘积只能是 \(\pm i\)、\(\pm j\)、\(\pm k\) 中的一个，不可能落在中心。于是这种情况下直接有
$$N(X,Y,Z)=0.$$
如果 \(X\)、\(Y\)、\(Z\) 全是偶数，记作 \(X=2x\)、\(Y=2y\)、\(Z=2z\)，那么
$$i^Xj^Yk^Z=(i^2)^x(j^2)^y(k^2)^z=(-1)^{x+y+z}=(-1)^{S/2}.$$
如果 \(X\)、\(Y\)、\(Z\) 全是奇数，记作 \(X=2x+1\)、\(Y=2y+1\)、\(Z=2z+1\)，那么
$$i^Xj^Yk^Z=(-1)^{x+y+z}ijk=(-1)^{x+y+z+1}=(-1)^{\lfloor S/2\rfloor}.$$
所以只要三者奇偶一致，就都可以统一写成
$$i^Xj^Yk^Z=(-1)^{\lfloor S/2\rfloor}.$$
从三个实现的分支结构可以看出，被计入 \(N(X,Y,Z)\) 的字符串，最终必须落在中心符号 \((-1)^S\) 上。把这个条件和上一节合并，就得到
$$(-1)^{\operatorname{inv}(w)}(-1)^{\lfloor S/2\rfloor}=(-1)^S.$$
也就是
$$\operatorname{inv}(w)\equiv \left\lceil\frac{S}{2}\right\rceil \pmod 2.$$
这样一来，\(N(X,Y,Z)\) 就不再是一个模糊的群论计数了。它精确地等于“满足给定逆序奇偶性”的多重集排列数量，前提是上一节的中心过滤条件已经通过。
记
$$T=\binom{S}{X,Y,Z},$$
再令 \(E\) 表示逆序数为偶的字符串个数，\(O\) 表示逆序数为奇的字符串个数，那么
$$E+O=T.$$
为了把这两部分拆开，可以看逆序生成函数
$$\sum_w q^{\operatorname{inv}(w)}=\binom{S}{X,Y,Z}_q,$$
也就是所有字符串对应的 \(q\)-多项式系数。在 \(q=-1\) 处取值，就得到
$$E-O=\binom{S}{X,Y,Z}_{q=-1}.$$
对本题来说，能通过中心过滤的情况非常整齐：
$$E-O=0 \quad \text{当 } X,Y,Z \text{ 全为奇数时},$$
而在全偶数情况下，如果记
$$H=\binom{S/2}{X/2,Y/2,Z/2},$$
那么就有
$$E-O=H.$$
这个半尺度多项式系数 \(H\)，正是实现里出现的修正项。
因为要求的逆序奇偶性是 \(\lceil S/2\rceil \bmod 2\)，所以最终计数可以直接写成
$$\boxed{ N(X,Y,Z)= \begin{cases} 0, & \text{如果 } X,Y,Z \text{ 的奇偶性不一致},\\[6pt] \dfrac{T}{2}, & \text{如果 } X,Y,Z \text{ 全为奇数},\\[10pt] \dfrac{T+H}{2}, & \text{如果 } X,Y,Z \text{ 全为偶数且 } S/2 \text{ 为偶数},\\[10pt] \dfrac{T-H}{2}, & \text{如果 } X,Y,Z \text{ 全为偶数且 } S/2 \text{ 为奇数}. \end{cases}}$$
这正好对应三个实现里完全一致的分支逻辑。
实现中使用了两个小规模但非常有代表性的精确检查值。
当 \((X,Y,Z)=(2,2,2)\) 时，\(S=6\)，有
$$T=\binom{6}{2,2,2}=90,\qquad H=\binom{3}{1,1,1}=6.$$
因为 \(S/2=3\) 是奇数，所以应该取逆序数为奇的那一半，于是
$$N(2,2,2)=\frac{90-6}{2}=42.$$
当 \((X,Y,Z)=(8,8,8)\) 时，\(S=24\)，于是
$$T=\binom{24}{8,8,8}=9465511770,\qquad H=\binom{12}{4,4,4}=34650.$$
此时 \(S/2=12\) 是偶数，所以应取逆序数为偶的那一半，从而
$$N(8,8,8)=\frac{9465511770+34650}{2}=4732773210.$$
这两个数都会在主模运算求和之前先被验证一次。
一旦 \(N(X,Y,Z)\) 的闭式已经得到，题目的总量就只是
$$A=\sum_{a=0}^{87}\sum_{b=0}^{87}\sum_{c=0}^{87} N(a^3,b^3,c^3)\pmod{888888883}.$$
由于立方运算保持奇偶性不变，所以只有 \(a,b,c\) 同奇或同偶时才可能有贡献。所有存活下来的项都可以在常数时间内用多项式系数公式得到。
三个实现都会先构造立方表 \(0^3,1^3,\dots,87^3\)，然后把阶乘和逆阶乘预处理到
$$3\cdot 87^3=1975509.$$
因为模数 \(M=888888883\) 是素数，所以逆元可以用费马小定理得到：
$$a^{-1}\equiv a^{M-2}\pmod M.$$
这样预处理完成后，每个多项式系数都能在 \(O(1)\) 时间内求出：
$$\binom{n}{a,b,c}\equiv n!\,(a!)^{-1}(b!)^{-1}(c!)^{-1}\pmod M.$$
对某个 \((X,Y,Z)\)，实现先做奇偶过滤。如果三者奇偶不一致，贡献立刻就是 0。否则先计算完整多项式系数 \(T\)；如果三者全偶，还会再计算半尺度的 \(H\)。在模 \(M\) 下除以 2，就是乘以 \(2^{-1}=(M+1)/2\)。
在主循环开始前，三个实现都会先用精确整数运算验证 \(N(2,2,2)=42\) 和 \(N(8,8,8)=4732773210\)。验证通过后，再用一个直接的串行三重循环遍历全部 \(88^3\) 个立方三元组并累加模贡献。
内存主要花在阶乘与逆阶乘表上，它们的长度达到 \(3\cdot 87^3\)。这一预处理阶段需要 \(O(87^3)\) 时间和 \(O(87^3)\) 空间。
最终求和一共访问 \(88^3\) 个三元组。由于表已经准备好，每次访问只做常数次运算，所以这一阶段是 \(O(88^3)\)。总体复杂度为
$$O(87^3+88^3)$$
时间和 \(O(87^3)\) 额外空间。前面的两个精确检查只占很小的一点开销。
Для каждого тройного индекса \((a,b,c)\) при \(0 \le a,b,c \le 87\) задаются значения \(X=a^3\), \(Y=b^3\) и \(Z=c^3\). Величина \(N(X,Y,Z)\) обозначает число слов, состоящих из \(X\) копий \(i\), \(Y\) копий \(j\) и \(Z\) копий \(k\), произведение которых в группе кватернионов \(Q_8\) попадает в требуемый центральный знак. Итоговая величина равна
$$A=\sum_{a=0}^{87}\sum_{b=0}^{87}\sum_{c=0}^{87} N(a^3,b^3,c^3)\pmod{888888883}.$$
Сложность задачи не в самом внешнем проходе по \(88^3\) тройкам. Настоящая трудность состоит в вычислении \(N(X,Y,Z)\) без перебора всех слов, число которых уже равно
$$\binom{X+Y+Z}{X,Y,Z}=\frac{(X+Y+Z)!}{X!\,Y!\,Z!}.$$
Решения обходят прямой перебор, переводя умножение в \(Q_8\) в условие на четность числа инверсий у перестановок мультимножества.
Обозначим \(S=X+Y+Z\). Любое слово с такими кратностями можно рассматривать как перестановку мультимножества из \(X\) букв \(i\), \(Y\) букв \(j\) и \(Z\) букв \(k\). Основная идея состоит в сравнении этого слова с канонически упорядоченным словом \(i^Xj^Yk^Z\).
Зафиксируем порядок \(i \lt j \lt k\). Для слова \(w\) пусть \(\operatorname{inv}(w)\) обозначает число пар, нарушающих этот порядок; то есть число случаев, когда \(j\) стоит раньше \(i\), плюс число случаев, когда \(k\) стоит раньше \(i\), плюс число случаев, когда \(k\) стоит раньше \(j\).
Каждая соседняя перестановка двух различных генераторов меняет знак, поскольку
$$ji=-ij,\qquad kj=-jk,\qquad ik=-ki.$$
Следовательно, после упорядочивания любого слова в вид \(i^Xj^Yk^Z\) получаем
$$P(w)=(-1)^{\operatorname{inv}(w)}\,i^Xj^Yk^Z.$$
Это и есть главный инвариант решения: значение слова определяется только фиксированным упорядоченным произведением \(i^Xj^Yk^Z\) и четностью числа инверсий.
В центре группы кватернионов лежат только элементы \(\pm1\). Поэтому сначала надо понять, когда \(i^Xj^Yk^Z\) вообще может быть центральным.
Если четности \(X\), \(Y\), \(Z\) не совпадают, то после сокращений остается одна из направленных компонент, и упорядоченное произведение оказывается равным одному из \(\pm i\), \(\pm j\), \(\pm k\). Тогда ни одно слово не дает вклад в \(N(X,Y,Z)\), а значит
$$N(X,Y,Z)=0.$$
Если \(X\), \(Y\), \(Z\) все четные, то есть \(X=2x\), \(Y=2y\), \(Z=2z\), имеем
$$i^Xj^Yk^Z=(i^2)^x(j^2)^y(k^2)^z=(-1)^{x+y+z}=(-1)^{S/2}.$$
Если \(X\), \(Y\), \(Z\) все нечетные, то есть \(X=2x+1\), \(Y=2y+1\), \(Z=2z+1\), то
$$i^Xj^Yk^Z=(-1)^{x+y+z}ijk=(-1)^{x+y+z+1}=(-1)^{\lfloor S/2\rfloor}.$$
Тем самым во всех допустимых случаях можно написать
$$i^Xj^Yk^Z=(-1)^{\lfloor S/2\rfloor}.$$
Условие, зашитое в реализациях, состоит в том, что слово должно давать центральный знак \((-1)^S\). Совместив это с предыдущей формулой, получаем
$$(-1)^{\operatorname{inv}(w)}(-1)^{\lfloor S/2\rfloor}=(-1)^S.$$
То есть
$$\operatorname{inv}(w)\equiv \left\lceil\frac{S}{2}\right\rceil \pmod 2.$$
Следовательно, \(N(X,Y,Z)\) есть в точности число перестановок мультимножества с заданной четностью числа инверсий, если только предварительный фильтр по центральности пройден.
Положим
$$T=\binom{S}{X,Y,Z},$$
а через \(E\) и \(O\) обозначим числа слов с четным и нечетным числом инверсий. Тогда
$$E+O=T.$$
Чтобы различить эти два количества, используют производящую функцию по инверсиям
$$\sum_w q^{\operatorname{inv}(w)}=\binom{S}{X,Y,Z}_q,$$
то есть \(q\)-мультиномиальный коэффициент для всех слов \(w\). Подстановка \(q=-1\) дает
$$E-O=\binom{S}{X,Y,Z}_{q=-1}.$$
В нашей задаче после фильтра по центральности остаются только два простых сценария:
$$E-O=0 \quad \text{если } X,Y,Z \text{ все нечетные},$$
а если \(X,Y,Z\) все четные и
$$H=\binom{S/2}{X/2,Y/2,Z/2},$$
то
$$E-O=H.$$
Именно этот полуразмерный мультиномиальный коэффициент \(H\) и выступает в коде как поправочный член.
Так как требуемая четность инверсий равна \(\lceil S/2\rceil \bmod 2\), окончательная формула такова:
$$\boxed{ N(X,Y,Z)= \begin{cases} 0, & \text{если } X,Y,Z \text{ не имеют одинаковой четности},\\[6pt] \dfrac{T}{2}, & \text{если } X,Y,Z \text{ все нечетные},\\[10pt] \dfrac{T+H}{2}, & \text{если } X,Y,Z \text{ все четные и } S/2 \text{ четно},\\[10pt] \dfrac{T-H}{2}, & \text{если } X,Y,Z \text{ все четные и } S/2 \text{ нечетно}. \end{cases}}$$
Эта схема ветвления в точности совпадает с тем, что делают реализации на C++, Python и Java.
Небольшие точные значения, используемые в реализациях, хорошо подтверждают вывод.
Для \((X,Y,Z)=(2,2,2)\) имеем \(S=6\),
$$T=\binom{6}{2,2,2}=90,\qquad H=\binom{3}{1,1,1}=6.$$
Так как \(S/2=3\) нечетно, нужно брать слова с нечетным числом инверсий, и потому
$$N(2,2,2)=\frac{90-6}{2}=42.$$
Для \((X,Y,Z)=(8,8,8)\) получаем \(S=24\),
$$T=\binom{24}{8,8,8}=9465511770,\qquad H=\binom{12}{4,4,4}=34650.$$
Теперь \(S/2=12\) четно, значит нужны слова с четным числом инверсий, и
$$N(8,8,8)=\frac{9465511770+34650}{2}=4732773210.$$
Именно эти два значения проверяются до запуска основной модульной суммы.
После получения явной формулы для \(N(X,Y,Z)\) искомая величина задачи Project Euler принимает вид
$$A=\sum_{a=0}^{87}\sum_{b=0}^{87}\sum_{c=0}^{87} N(a^3,b^3,c^3)\pmod{888888883}.$$
Поскольку возведение в куб сохраняет четность, вклад возможен только тогда, когда \(a\), \(b\), \(c\) одновременно четны или одновременно нечетны. Каждый оставшийся член вычисляется затем за \(O(1)\) по формуле с мультиномиальными коэффициентами.
Сначала реализации строят список кубов \(0^3,1^3,\dots,87^3\), а затем предвычисляют факториалы и обратные факториалы вплоть до
$$3\cdot 87^3=1975509.$$
Так как модуль \(M=888888883\) прост, обратные элементы вычисляются по малой теореме Ферма:
$$a^{-1}\equiv a^{M-2}\pmod M.$$
После этого любой мультиномиальный коэффициент считается за \(O(1)\):
$$\binom{n}{a,b,c}\equiv n!\,(a!)^{-1}(b!)^{-1}(c!)^{-1}\pmod M.$$
Для одной тройки \((X,Y,Z)\) код сначала применяет фильтр по четности. Если четности различны, вклад немедленно равен нулю. Иначе вычисляется полный мультиномиальный коэффициент \(T\), а в полностью четном случае также полуразмерный коэффициент \(H\). Деление на 2 по модулю выполняется как умножение на \(2^{-1}=(M+1)/2\).
Перед основным накоплением реализации отдельно проверяют точные целочисленные равенства \(N(2,2,2)=42\) и \(N(8,8,8)=4732773210\). Затем запускается обычный последовательный тройной цикл по всем \(88^3\) кубическим тройкам, и модульный вклад каждой из них добавляется к ответу.
Основной объем памяти занимают таблицы факториалов и обратных факториалов, размер которых доходит до \(3\cdot 87^3\). Этап предвычисления стоит \(O(87^3)\) по времени и \(O(87^3)\) по памяти.
Финальный проход посещает ровно \(88^3\) троек. После предвычисления каждая тройка обрабатывается за константное время, значит эта фаза имеет сложность \(O(88^3)\). В сумме алгоритм работает за
$$O(87^3+88^3)$$
времени и использует \(O(87^3)\) дополнительной памяти. Точные проверочные вычисления по сравнению с этими двумя фазами ничтожны.
لكل ثلاثي \((a,b,c)\) حيث \(0 \le a,b,c \le 87\)، نعرّف \(X=a^3\) و\(Y=b^3\) و\(Z=c^3\). وتمثل الكمية \(N(X,Y,Z)\) عدد الكلمات المؤلفة من \(X\) نسخة من \(i\) و\(Y\) نسخة من \(j\) و\(Z\) نسخة من \(k\) بحيث يكون حاصل ضربها في مجموعة الكواتيرنيون \(Q_8\) هو الإشارة المركزية المطلوبة. والمطلوب في النهاية هو
$$A=\sum_{a=0}^{87}\sum_{b=0}^{87}\sum_{c=0}^{87} N(a^3,b^3,c^3)\pmod{888888883}.$$
الصعوبة الحقيقية ليست في المرور على \(88^3\) ثلاثيات، بل في حساب \(N(X,Y,Z)\) من غير تعداد جميع الكلمات الممكنة، لأن العدد الخام لهذه الكلمات يساوي
$$\binom{X+Y+Z}{X,Y,Z}=\frac{(X+Y+Z)!}{X!\,Y!\,Z!}.$$
الحل يحول هذه المسألة من حساب مباشر في المجموعة إلى مسألة عن زوجية عدد الانعكاسات في تبديلات متعددة العناصر.
لنكتب \(S=X+Y+Z\). كل كلمة بهذه الأعداد من الحروف يمكن النظر إليها على أنها تبديل للمتعدد الذي يحتوي على \(X\) حرفًا من \(i\) و\(Y\) حرفًا من \(j\) و\(Z\) حرفًا من \(k\). الفكرة الأساسية هي مقارنة أي كلمة بالكلمة المرتبة ترتيبا قانونيا \(i^Xj^Yk^Z\).
نثبت الترتيب \(i \lt j \lt k\). ولأي كلمة \(w\)، نرمز بـ \(\operatorname{inv}(w)\) إلى عدد الأزواج التي تخالف هذا الترتيب؛ أي عدد مرات ظهور \(j\) قبل \(i\)، مضافًا إليه عدد مرات ظهور \(k\) قبل \(i\)، ثم عدد مرات ظهور \(k\) قبل \(j\).
كل تبادل مجاور بين مولدين مختلفين يقلب الإشارة لأن
$$ji=-ij,\qquad kj=-jk,\qquad ik=-ki.$$
لذلك عندما نعيد ترتيب أي كلمة إلى الشكل \(i^Xj^Yk^Z\) نحصل على
$$P(w)=(-1)^{\operatorname{inv}(w)}\,i^Xj^Yk^Z.$$
هذه هي الثابتة الأساسية في الحل كله: قيمة الكلمة في المجموعة تعتمد فقط على الجداء المرتب الثابت \(i^Xj^Yk^Z\) وعلى زوجية عدد الانعكاسات.
العناصر المركزية الوحيدة في \(Q_8\) هي \(\pm1\). لذلك يجب أولًا أن نعرف متى يمكن أن يكون \(i^Xj^Yk^Z\) عنصرًا مركزيًا أصلًا.
إذا لم تكن زوجيات \(X\) و\(Y\) و\(Z\) متطابقة، فإن أحد الاتجاهات \(i\) أو \(j\) أو \(k\) يبقى حاضرًا، ويصبح الجداء المرتب واحدًا من \(\pm i\) أو \(\pm j\) أو \(\pm k\). عندها لا يمكن لأي كلمة أن تسهم في \(N(X,Y,Z)\)، وبالتالي
$$N(X,Y,Z)=0.$$
أما إذا كانت \(X\) و\(Y\) و\(Z\) كلها زوجية، فنكتب \(X=2x\) و\(Y=2y\) و\(Z=2z\)، وعندئذ
$$i^Xj^Yk^Z=(i^2)^x(j^2)^y(k^2)^z=(-1)^{x+y+z}=(-1)^{S/2}.$$
وإذا كانت كلها فردية، أي \(X=2x+1\) و\(Y=2y+1\) و\(Z=2z+1\)، فإن
$$i^Xj^Yk^Z=(-1)^{x+y+z}ijk=(-1)^{x+y+z+1}=(-1)^{\lfloor S/2\rfloor}.$$
إذن في كل حالة تبقى فيها الزوجيات متساوية، نستطيع كتابة الجداء المرتب بالشكل المختصر
$$i^Xj^Yk^Z=(-1)^{\lfloor S/2\rfloor}.$$
الشرط الذي تطبقه الحلول هو أن تكون قيمة الكلمة هي الإشارة المركزية \((-1)^S\). وبدمج هذا مع العلاقة السابقة نحصل على
$$(-1)^{\operatorname{inv}(w)}(-1)^{\lfloor S/2\rfloor}=(-1)^S.$$
وهذا يكافئ
$$\operatorname{inv}(w)\equiv \left\lceil\frac{S}{2}\right\rceil \pmod 2.$$
وهكذا تصبح \(N(X,Y,Z)\) ببساطة عدد تبديلات المتعدد التي تملك زوجية انعكاسات محددة سلفًا، بشرط أن تكون حالة المركزية منطبقة أصلًا.
لنضع
$$T=\binom{S}{X,Y,Z},$$
ولنرمز بـ \(E\) إلى عدد الكلمات ذات عدد الانعكاسات الزوجي، وبـ \(O\) إلى عدد الكلمات ذات عدد الانعكاسات الفردي. إذن
$$E+O=T.$$
ولفصل هذين العددين نستخدم الدالة المولدة للانعكاسات
$$\sum_w q^{\operatorname{inv}(w)}=\binom{S}{X,Y,Z}_q,$$
أي المعامل متعدد الحدود من نوع \(q\). وعند التعويض بـ \(q=-1\) نحصل على
$$E-O=\binom{S}{X,Y,Z}_{q=-1}.$$
في هذه المسألة، الحالات التي تتجاوز مرشح المركزية سهلة جدًا:
إذا كانت \(X\) و\(Y\) و\(Z\) كلها فردية، فإن الفرق بين العدّين يساوي
$$E-O=0.$$
أما إذا كانت كلها زوجية، ومع تعريف
$$H=\binom{S/2}{X/2,Y/2,Z/2},$$
فإن
$$E-O=H.$$
وهذا المعامل نصفي القياس \(H\) هو بالضبط حد التصحيح الذي تستعمله التطبيقات.
بما أن زوجية الانعكاسات المطلوبة هي \(\lceil S/2\rceil \bmod 2\)، فإن العدد النهائي يساوي
إذا لم تكن \(X\) و\(Y\) و\(Z\) متطابقة الزوجية، فلدينا
$$N(X,Y,Z)=0.$$
إذا كانت كلها فردية، فلدينا
$$N(X,Y,Z)=\frac{T}{2}.$$
إذا كانت كلها زوجية وكان \(S/2\) زوجيًا، فلدينا
$$N(X,Y,Z)=\frac{T+H}{2}.$$
إذا كانت كلها زوجية وكان \(S/2\) فرديًا، فلدينا
$$N(X,Y,Z)=\frac{T-H}{2}.$$
وهذه هي تمامًا البنية الفرعية نفسها الموجودة في تطبيقات C++ وPython وJava.
القيم الصغيرة التي تتحقق منها التطبيقات تعطي اختبارًا واضحًا لصحة الاشتقاق.
بالنسبة إلى \((X,Y,Z)=(2,2,2)\)، لدينا \(S=6\)، ومن ثم
$$T=\binom{6}{2,2,2}=90,\qquad H=\binom{3}{1,1,1}=6.$$
وبما أن \(S/2=3\) فردي، فنحن نعد الكلمات ذات الانعكاسات الفردية، وبالتالي
$$N(2,2,2)=\frac{90-6}{2}=42.$$
وبالنسبة إلى \((X,Y,Z)=(8,8,8)\)، نحصل على \(S=24\)، ولذلك
$$T=\binom{24}{8,8,8}=9465511770,\qquad H=\binom{12}{4,4,4}=34650.$$
وهنا \(S/2=12\) زوجي، فنعد الكلمات ذات الانعكاسات الزوجية، فنصل إلى
$$N(8,8,8)=\frac{9465511770+34650}{2}=4732773210.$$
هاتان هما القيمتان الدقيقتان اللتان يتم التحقق منهما قبل بدء الجمع المعياري الرئيسي.
بعد الحصول على الصيغة المغلقة لـ \(N(X,Y,Z)\)، تصبح كمية Project Euler المطلوبة هي
$$A=\sum_{a=0}^{87}\sum_{b=0}^{87}\sum_{c=0}^{87} N(a^3,b^3,c^3)\pmod{888888883}.$$
ولأن التكعيب يحافظ على الزوجية، فلا تسهم إلا الثلاثيات التي يكون فيها \(a\) و\(b\) و\(c\) كلها زوجية أو كلها فردية. وكل حد ناجٍ بعد ذلك يُحسب في زمن ثابت باستعمال معاملات متعددة الحدود.
تبدأ التطبيقات ببناء قائمة المكعبات \(0^3,1^3,\dots,87^3\)، ثم تحسب مسبقًا جدولي المضروبات ومعكوسات المضروبات حتى
$$3\cdot 87^3=1975509.$$
ولأن المودولو \(M=888888883\) عدد أولي، فإن المعكوسات تُحسب باستعمال مبرهنة فيرما الصغرى:
$$a^{-1}\equiv a^{M-2}\pmod M.$$
وبعد هذه التهيئة يصبح حساب أي معامل متعدد الحدود مسألة \(O(1)\):
$$\binom{n}{a,b,c}\equiv n!\,(a!)^{-1}(b!)^{-1}(c!)^{-1}\pmod M.$$
بالنسبة إلى ثلاثي واحد \((X,Y,Z)\)، تطبق الشيفرة أولًا مرشح الزوجية. إذا لم تتطابق الزوجيات فالمساهمة تساوي صفرًا فورًا. وإلا يُحسب المعامل الكامل \(T\)، وفي الحالة الزوجية الكاملة يُحسب أيضًا المعامل نصفي القياس \(H\). أما القسمة على 2 في الحساب المعياري فتنفذ بالضرب في \(2^{-1}=(M+1)/2\).
وقبل بدء الجمع الرئيسي، تتحقق التطبيقات بالحساب الصحيح من القيمتين \(N(2,2,2)=42\) و\(N(8,8,8)=4732773210\). وبعد ذلك تُشغَّل حلقة ثلاثية متسلسلة على جميع الثلاثيات المكعبة \(88^3\)، وتُضاف المساهمة المعيارية لكل حالة إلى الجواب.
تهيمن جداول المضروبات ومعكوساتها على استهلاك الذاكرة، وهي تمتد حتى \(3\cdot 87^3\). لذلك تكلف مرحلة التهيئة المسبقة \(O(87^3)\) زمنًا و\(O(87^3)\) ذاكرة.
أما المرحلة النهائية فتمر على \(88^3\) ثلاثيات بالضبط. وبعد اكتمال التهيئة تصبح كلفة كل ثلاثي ثابتة، فتكون هذه المرحلة \(O(88^3)\). بالتالي فإن التعقيد الكلي هو
$$O(87^3+88^3)$$
زمنًا مع \(O(87^3)\) ذاكرة إضافية. أما حسابات التحقق الدقيقة فهي صغيرة جدًا مقارنة بهاتين المرحلتين الرئيسيتين.