Let
$$M=11^8=214358881.$$
Problem 478 asks for a counting function \(E(n)\) modulo \(M\). The C++, Python, and Java implementations do not enumerate mixtures directly. Instead, they rewrite the problem in terms of primitive lattice data inside the box \(\{0,\dots,n\}^3\), together with a structured correction built from reduced two-dimensional boundary profiles. The checkpoints embedded in the implementations are
$$E(1)=103,\qquad E(2)=520447,\qquad E(10)=82608406,\qquad E(500)=13801403.$$
The solution is organized around four quantities: a primitive-triple count \(P(n)\), a totient-weighted boundary term \(B(n)\), an inner coprime count \(H_n(L)\), and a final correction sum \(T(n)\). Once these are known, \(E(n)\) follows from one modular identity.
Define
$$\mathcal{P}_n=\left\{(a,b,c)\in \{0,\dots,n\}^3\setminus\{(0,0,0)\}:\gcd(a,b,c)=1\right\}.$$
Let \(P(n)=|\mathcal{P}_n|\). For a fixed divisor \(d\), the number of triples whose three coordinates are all divisible by \(d\) is
$$\left(\left\lfloor\frac{n}{d}\right\rfloor+1\right)^3-1,$$
because each coordinate has \(\left\lfloor n/d\right\rfloor+1\) multiples of \(d\) between \(0\) and \(n\), and the all-zero triple must be excluded. Möbius inversion therefore gives
$$P(n)=\sum_{d=1}^{n}\mu(d)\left(\left(\left\lfloor\frac{n}{d}\right\rfloor+1\right)^3-1\right).$$
This is the large exponent that ultimately drives the term \(2^{P(n)}\).
For each \(L\ge 1\), the number of coprime positive pairs \((u,v)\) satisfying
$$u+v=L$$
is exactly \(\varphi(L)\). Indeed, choosing \(u\) with \(1\le u\le L\) and \(\gcd(u,L)=1\) uniquely determines \(v=L-u\), and then \(\gcd(u,v)=1\). The final correction formula treats each reduced pair with a sixfold symmetry factor, which explains the coefficient \(6\varphi(L)\). Hence the global boundary weight is
$$B(n)=6\sum_{L=1}^{n}\varphi(L).$$
Fix \(L\). The inner quantity used by the implementations is
$$H_n(L)=1+\#\left\{(x,y)\in \mathbb{Z}_{>0}^2:\gcd(x,y)=1,\ x+Ly\le n\right\}.$$
The initial \(1\) is explicit. To count the coprime pairs, apply Möbius inversion again. If \(d\mid x\) and \(d\mid y\), write \(x=da\), \(y=db\). Then
$$a+Lb\le \left\lfloor\frac{n}{d}\right\rfloor.$$
Set
$$q_{d,L}=\left\lfloor\frac{n}{Ld}\right\rfloor.$$
For each \(b=1,\dots,q_{d,L}\), the variable \(a\) has \(\left\lfloor n/d\right\rfloor-Lb\) positive choices. Summing those choices yields
$$\sum_{b=1}^{q_{d,L}}\left(\left\lfloor\frac{n}{d}\right\rfloor-Lb\right)=q_{d,L}\left\lfloor\frac{n}{d}\right\rfloor-L\frac{q_{d,L}(q_{d,L}+1)}{2}.$$
Therefore
$$H_n(L)=1+\sum_{d=1}^{\lfloor n/L\rfloor}\mu(d)\left(q_{d,L}\left\lfloor\frac{n}{d}\right\rfloor-L\frac{q_{d,L}(q_{d,L}+1)}{2}\right),\qquad q_{d,L}=\left\lfloor\frac{n}{Ld}\right\rfloor.$$
This matches the inner loop term for term.
Let
$$\Phi=\varphi(M)=\varphi(11^8)=10\cdot 11^7=194871710.$$
Since \(\gcd(2,M)=1\), Euler's theorem tells us that every power of \(2\) modulo \(M\) depends only on the exponent modulo \(\Phi\). However, the correction terms also need
$$\frac{P(n)-1}{2},$$
so parity must be preserved before dividing by \(2\). For that reason the implementations first compute \(P(n)\) modulo \(2\Phi\), subtract \(1\), verify that the result is even, and only then reduce the half-exponent modulo \(\Phi\).
With the reduced half-exponent understood modulo \(\Phi\), define
$$T(n)\equiv \sum_{L=1}^{n}6\varphi(L)\,2^{\left(\frac{P(n)-1}{2}-H_n(L)\right)\bmod \Phi}\pmod{M}.$$
Then the shared correction term is
$$F(n)\equiv 1+2^{\left(\frac{P(n)-1}{2}\right)\bmod \Phi}B(n)-T(n)\pmod{M},$$
and the required value is
$$\boxed{E(n)\equiv 2^{P(n)\bmod \Phi}-F(n)\pmod{M}.}$$
This is exactly the algebra used by the C++, Python, and Java implementations.
Here \(\mu(1)=1\), \(\mu(2)=-1\), and \(\varphi(1)=\varphi(2)=1\). First,
$$P(2)=\left((2+1)^3-1\right)-\left((1+1)^3-1\right)=26-7=19.$$
Hence
$$\frac{P(2)-1}{2}=9,\qquad B(2)=6(\varphi(1)+\varphi(2))=12.$$
For \(L=1\),
$$H_2(1)=1+\left(2\cdot 2-\frac{2\cdot 3}{2}\right)-\left(1\cdot 1-\frac{1\cdot 2}{2}\right)=2.$$
For \(L=2\),
$$H_2(2)=1+\left(1\cdot 2-2\cdot\frac{1\cdot 2}{2}\right)=1.$$
Therefore
$$T(2)=6\cdot 2^{9-2}+6\cdot 2^{9-1}=768+1536=2304,$$
$$F(2)=1+2^9\cdot 12-2304=3841,$$
and
$$E(2)=2^{19}-3841=520447.$$
This agrees with the checkpoint used by the implementations.
The implementation begins with a linear sieve that computes \(\mu(d)\) and \(\varphi(d)\) for every \(d\le n\). It also stores the quotient table \(\left\lfloor n/d\right\rfloor\), because both \(P(n)\) and \(H_n(L)\) reuse those values constantly.
Next, it evaluates \(P(n)\) modulo \(2\Phi\), accumulates \(B(n)=6\sum_{L\le n}\varphi(L)\) modulo \(M\), and precomputes modular powers of \(2\) by fast exponentiation. Then it loops over \(L=1,\dots,n\), computes \(H_n(L)\) from the Möbius-weighted inner sum, forms the exponent residue \(\left(\frac{P(n)-1}{2}-H_n(L)\right)\bmod \Phi\), and adds the contribution \(6\varphi(L)\cdot 2^{((P(n)-1)/2-H_n(L))\bmod\Phi}\) to \(T(n)\).
Finally, it combines \(2^{P(n)\bmod\Phi}\), \(B(n)\), and \(T(n)\) through the boxed formula. The C++ and Java implementations parallelize the outer loop over \(L\), while the Python implementation evaluates the same mathematics sequentially.
The sieve and the quotient table each take \(O(n)\) time and \(O(n)\) memory. The dominant cost is the double summation for \(H_n(L)\): for each \(L\), the inner loop runs to \(\left\lfloor n/L\right\rfloor\), so the total number of iterations is
$$\sum_{L=1}^{n}\left\lfloor\frac{n}{L}\right\rfloor=n\log n+O(n).$$
Thus the overall running time is \(O(n\log n)\) arithmetic operations with \(O(n)\) memory. Parallel execution improves wall-clock time but does not change the asymptotic complexity.
Sei
$$M=11^8=214358881.$$
Gesucht ist der Zählwert \(E(n)\) modulo \(M\). Die C++-, Python- und Java-Implementierungen enumerieren die Mischungen nicht direkt. Stattdessen wird das Problem auf primitive Gitterdaten im Würfel \(\{0,\dots,n\}^3\) und auf einen strukturierten Korrekturterm zurückgeführt, der aus reduzierten zweidimensionalen Randprofilen besteht. Die in den Implementierungen eingebauten Kontrollwerte sind
$$E(1)=103,\qquad E(2)=520447,\qquad E(10)=82608406,\qquad E(500)=13801403.$$
Die Lösung ist um vier Größen herum aufgebaut: die Zahl primitiver Tripel \(P(n)\), ein mit der Eulerschen Phi-Funktion gewichteter Randterm \(B(n)\), eine innere teilerfremde Zählfunktion \(H_n(L)\) und eine abschließende Korrektursumme \(T(n)\). Sind diese vier Größen bekannt, folgt \(E(n)\) aus einer einzigen modularen Identität.
Definiere
$$\mathcal{P}_n=\left\{(a,b,c)\in \{0,\dots,n\}^3\setminus\{(0,0,0)\}:\gcd(a,b,c)=1\right\}.$$
Sei \(P(n)=|\mathcal{P}_n|\). Für einen festen Teiler \(d\) ist die Anzahl der Tripel, deren drei Koordinaten alle durch \(d\) teilbar sind, gleich
$$\left(\left\lfloor\frac{n}{d}\right\rfloor+1\right)^3-1,$$
denn jede Koordinate besitzt \(\left\lfloor n/d\right\rfloor+1\) Vielfache von \(d\) zwischen \(0\) und \(n\), und das Nulltripel muss entfernt werden. Mit Möbius-Inversion erhält man daher
$$P(n)=\sum_{d=1}^{n}\mu(d)\left(\left(\left\lfloor\frac{n}{d}\right\rfloor+1\right)^3-1\right).$$
Genau dieser Wert liefert später den großen Exponenten in \(2^{P(n)}\).
Für jedes \(L\ge 1\) ist die Anzahl teilerfremder positiver Paare \((u,v)\) mit
$$u+v=L$$
genau \(\varphi(L)\). Wählt man nämlich \(u\) mit \(1\le u\le L\) und \(\gcd(u,L)=1\), dann ist \(v=L-u\) eindeutig bestimmt und automatisch \(\gcd(u,v)=1\). In der Endformel tritt zu jedem reduzierten Paar ein sechsfacher Symmetriefaktor auf; daher erklärt sich der Koeffizient \(6\varphi(L)\). Damit lautet das globale Randgewicht
$$B(n)=6\sum_{L=1}^{n}\varphi(L).$$
Fixiere \(L\). Die von den Implementierungen verwendete innere Größe ist
$$H_n(L)=1+\#\left\{(x,y)\in \mathbb{Z}_{>0}^2:\gcd(x,y)=1,\ x+Ly\le n\right\}.$$
Die führende \(1\) steht explizit in der Formel. Um die teilerfremden Paare zu zählen, verwendet man erneut Möbius-Inversion. Gilt \(d\mid x\) und \(d\mid y\), dann schreibt man \(x=da\), \(y=db\). Damit folgt
$$a+Lb\le \left\lfloor\frac{n}{d}\right\rfloor.$$
Setze
$$q_{d,L}=\left\lfloor\frac{n}{Ld}\right\rfloor.$$
Für jedes \(b=1,\dots,q_{d,L}\) besitzt \(a\) genau \(\left\lfloor n/d\right\rfloor-Lb\) positive Möglichkeiten. Die Summe dieser Möglichkeiten ist
$$\sum_{b=1}^{q_{d,L}}\left(\left\lfloor\frac{n}{d}\right\rfloor-Lb\right)=q_{d,L}\left\lfloor\frac{n}{d}\right\rfloor-L\frac{q_{d,L}(q_{d,L}+1)}{2}.$$
Daher gilt
$$H_n(L)=1+\sum_{d=1}^{\lfloor n/L\rfloor}\mu(d)\left(q_{d,L}\left\lfloor\frac{n}{d}\right\rfloor-L\frac{q_{d,L}(q_{d,L}+1)}{2}\right),\qquad q_{d,L}=\left\lfloor\frac{n}{Ld}\right\rfloor.$$
Das ist exakt dieselbe innere Summe wie im Programm.
Sei
$$\Phi=\varphi(M)=\varphi(11^8)=10\cdot 11^7=194871710.$$
Wegen \(\gcd(2,M)=1\) hängt jede Potenz von \(2\) modulo \(M\) nur vom Exponenten modulo \(\Phi\) ab. Die Korrekturterme benötigen jedoch auch
$$\frac{P(n)-1}{2},$$
also muss die Parität vor der Division durch \(2\) erhalten bleiben. Deshalb berechnen die Implementierungen zunächst \(P(n)\) modulo \(2\Phi\), ziehen \(1\) ab, prüfen die Geradheit und reduzieren erst danach den halbierten Exponenten modulo \(\Phi\).
Mit dem halbierten Exponenten modulo \(\Phi\) definiert man
$$T(n)\equiv \sum_{L=1}^{n}6\varphi(L)\,2^{\left(\frac{P(n)-1}{2}-H_n(L)\right)\bmod \Phi}\pmod{M}.$$
Dann ist der gemeinsame Korrekturterm
$$F(n)\equiv 1+2^{\left(\frac{P(n)-1}{2}\right)\bmod \Phi}B(n)-T(n)\pmod{M},$$
und der gesuchte Wert lautet
$$\boxed{E(n)\equiv 2^{P(n)\bmod \Phi}-F(n)\pmod{M}.}$$
Genau diese Algebra setzen die C++-, Python- und Java-Implementierungen um.
Hier gelten \(\mu(1)=1\), \(\mu(2)=-1\) und \(\varphi(1)=\varphi(2)=1\). Zuerst erhält man
$$P(2)=\left((2+1)^3-1\right)-\left((1+1)^3-1\right)=26-7=19.$$
Also
$$\frac{P(2)-1}{2}=9,\qquad B(2)=6(\varphi(1)+\varphi(2))=12.$$
Für \(L=1\) ist
$$H_2(1)=1+\left(2\cdot 2-\frac{2\cdot 3}{2}\right)-\left(1\cdot 1-\frac{1\cdot 2}{2}\right)=2.$$
Für \(L=2\) folgt
$$H_2(2)=1+\left(1\cdot 2-2\cdot\frac{1\cdot 2}{2}\right)=1.$$
Damit ist
$$T(2)=6\cdot 2^{9-2}+6\cdot 2^{9-1}=768+1536=2304,$$
$$F(2)=1+2^9\cdot 12-2304=3841,$$
und somit
$$E(2)=2^{19}-3841=520447.$$
Das stimmt mit dem eingebauten Kontrollwert überein.
Die Implementierung beginnt mit einem linearen Sieb für \(\mu(d)\) und \(\varphi(d)\) für alle \(d\le n\). Zusätzlich wird die Quotiententabelle \(\left\lfloor n/d\right\rfloor\) gespeichert, weil sowohl \(P(n)\) als auch \(H_n(L)\) diese Werte ständig benötigen.
Anschließend wird \(P(n)\) modulo \(2\Phi\) berechnet, parallel dazu \(B(n)=6\sum_{L\le n}\varphi(L)\) modulo \(M\) akkumuliert, und Potenzen von \(2\) werden per schneller modularer Exponentiation ausgewertet. Danach läuft eine Schleife über \(L=1,\dots,n\): Für jedes \(L\) wird \(H_n(L)\) aus der Möbius-gewichteten inneren Summe bestimmt, der Exponentenrest \(\left(\frac{P(n)-1}{2}-H_n(L)\right)\bmod \Phi\) gebildet und der Beitrag \(6\varphi(L)\cdot 2^{((P(n)-1)/2-H_n(L))\bmod\Phi}\) zu \(T(n)\) addiert.
Zum Schluss werden \(2^{P(n)\bmod\Phi}\), \(B(n)\) und \(T(n)\) durch die eingerahmte Formel kombiniert. Die C++- und Java-Implementierungen parallelisieren die äußere Schleife über \(L\), während die Python-Implementierung dieselbe Mathematik sequentiell ausführt.
Das Sieb und die Quotiententabelle benötigen jeweils \(O(n)\) Zeit und \(O(n)\) Speicher. Der dominante Aufwand entsteht in der Doppelsumme für \(H_n(L)\): Für jedes \(L\) läuft die innere Schleife bis \(\left\lfloor n/L\right\rfloor\), also insgesamt
$$\sum_{L=1}^{n}\left\lfloor\frac{n}{L}\right\rfloor=n\log n+O(n).$$
Damit beträgt die Gesamtlaufzeit \(O(n\log n)\) arithmetische Operationen bei \(O(n)\) Speicher. Parallelisierung verbessert nur die reale Laufzeit, nicht die asymptotische Schranke.
Şunu tanımlayalım:
$$M=11^8=214358881.$$
İstenen değer \(E(n)\)'dir ve sonuç \(M\) modunda alınır. C++, Python ve Java uygulamaları karışımları doğrudan tek tek saymaz. Bunun yerine problem, \(\{0,\dots,n\}^3\) kutusu içindeki primitif kafes verilerine ve indirgenmiş iki boyutlu sınır profillerinden gelen yapılandırılmış bir düzeltme terimine dönüştürülür. Uygulamaların içine yerleştirilmiş kontrol değerleri şunlardır:
$$E(1)=103,\qquad E(2)=520447,\qquad E(10)=82608406,\qquad E(500)=13801403.$$
Çözüm dört nicelik etrafında kuruludur: primitif üçlü sayısı \(P(n)\), \(\varphi\) ile ağırlıklandırılmış sınır terimi \(B(n)\), içteki aralarında asal sayım \(H_n(L)\) ve son düzeltme toplamı \(T(n)\). Bu dört parça elde edilince \(E(n)\) tek bir modüler bağıntıdan çıkar.
Şu kümeyi tanımlayalım:
$$\mathcal{P}_n=\left\{(a,b,c)\in \{0,\dots,n\}^3\setminus\{(0,0,0)\}:\gcd(a,b,c)=1\right\}.$$
\(P(n)=|\mathcal{P}_n|\) olsun. Sabit bir \(d\) böleni için, üç koordinatın da \(d\) ile bölündüğü üçlü sayısı
$$\left(\left\lfloor\frac{n}{d}\right\rfloor+1\right)^3-1$$
olur; çünkü her koordinatın \(0\) ile \(n\) arasında \(\left\lfloor n/d\right\rfloor+1\) tane uygun değeri vardır ve \((0,0,0)\) çıkarılmalıdır. Bu nedenle Möbius terslemesi
$$P(n)=\sum_{d=1}^{n}\mu(d)\left(\left(\left\lfloor\frac{n}{d}\right\rfloor+1\right)^3-1\right)$$
formülünü verir. Sonradan görünen \(2^{P(n)}\) teriminin büyük üssü tam olarak budur.
Her \(L\ge 1\) için
$$u+v=L$$
koşulunu sağlayan ve aralarında asal olan pozitif \((u,v)\) çiftlerinin sayısı tam olarak \(\varphi(L)\)'dir. Gerçekten, \(1\le u\le L\) aralığında \(\gcd(u,L)=1\) olan bir \(u\) seçildiğinde \(v=L-u\) tekil olarak belirlenir ve \(\gcd(u,v)=1\) olur. Son düzeltme toplamında her indirgenmiş çifte altılı bir simetri katsayısı uygulanır; bu da \(6\varphi(L)\) katsayısını açıklar. Böylece küresel sınır ağırlığı
$$B(n)=6\sum_{L=1}^{n}\varphi(L)$$
şeklindedir.
Sabit bir \(L\) alalım. Uygulamaların kullandığı iç nicelik
$$H_n(L)=1+\#\left\{(x,y)\in \mathbb{Z}_{>0}^2:\gcd(x,y)=1,\ x+Ly\le n\right\}$$
şeklindedir. Baştaki \(1\) doğrudan formülün bir parçasıdır. Aralarında asal çiftleri saymak için yine Möbius terslemesi uygulanır. Eğer \(d\mid x\) ve \(d\mid y\) ise \(x=da\), \(y=db\) yazılır. Bu durumda
$$a+Lb\le \left\lfloor\frac{n}{d}\right\rfloor$$
olur. Şimdi
$$q_{d,L}=\left\lfloor\frac{n}{Ld}\right\rfloor$$
tanımını yapalım. Her \(b=1,\dots,q_{d,L}\) için \(a\)'nın \(\left\lfloor n/d\right\rfloor-Lb\) tane pozitif seçeneği vardır. Bunların toplamı
$$\sum_{b=1}^{q_{d,L}}\left(\left\lfloor\frac{n}{d}\right\rfloor-Lb\right)=q_{d,L}\left\lfloor\frac{n}{d}\right\rfloor-L\frac{q_{d,L}(q_{d,L}+1)}{2}$$
eder. Dolayısıyla
$$H_n(L)=1+\sum_{d=1}^{\lfloor n/L\rfloor}\mu(d)\left(q_{d,L}\left\lfloor\frac{n}{d}\right\rfloor-L\frac{q_{d,L}(q_{d,L}+1)}{2}\right),\qquad q_{d,L}=\left\lfloor\frac{n}{Ld}\right\rfloor.$$
Bu, uygulamadaki iç döngünün terim terim aynı halidir.
Şunu yazalım:
$$\Phi=\varphi(M)=\varphi(11^8)=10\cdot 11^7=194871710.$$
\(\gcd(2,M)=1\) olduğu için Euler teoremi, \(M\) modunda her \(2\) kuvvetinin yalnızca üssün \(\Phi\) modundaki değerine bağlı olduğunu söyler. Fakat düzeltme terimleri ayrıca
$$\frac{P(n)-1}{2}$$
ifadesini gerektirir; yani \(2\)'ye bölmeden önce tek-çift bilgisi korunmalıdır. Bu yüzden uygulamalar önce \(P(n)\)'yi \(2\Phi\) modunda hesaplar, sonra \(1\) çıkarır, sonucun çift olduğunu doğrular ve ancak bundan sonra yarıya bölüp \(\Phi\) moduna indirger.
Yarılanmış üs \(\Phi\) modunda anlaşıldığında
$$T(n)\equiv \sum_{L=1}^{n}6\varphi(L)\,2^{\left(\frac{P(n)-1}{2}-H_n(L)\right)\bmod \Phi}\pmod{M}$$
tanımlanır. Ortak düzeltme terimi
$$F(n)\equiv 1+2^{\left(\frac{P(n)-1}{2}\right)\bmod \Phi}B(n)-T(n)\pmod{M}$$
olur ve istenen sonuç
$$\boxed{E(n)\equiv 2^{P(n)\bmod \Phi}-F(n)\pmod{M}.}$$
C++, Python ve Java uygulamalarının kullandığı cebir tam olarak budur.
Burada \(\mu(1)=1\), \(\mu(2)=-1\) ve \(\varphi(1)=\varphi(2)=1\). Önce
$$P(2)=\left((2+1)^3-1\right)-\left((1+1)^3-1\right)=26-7=19$$
bulunur. Dolayısıyla
$$\frac{P(2)-1}{2}=9,\qquad B(2)=6(\varphi(1)+\varphi(2))=12.$$
\(L=1\) için
$$H_2(1)=1+\left(2\cdot 2-\frac{2\cdot 3}{2}\right)-\left(1\cdot 1-\frac{1\cdot 2}{2}\right)=2,$$
\(L=2\) için ise
$$H_2(2)=1+\left(1\cdot 2-2\cdot\frac{1\cdot 2}{2}\right)=1.$$
Böylece
$$T(2)=6\cdot 2^{9-2}+6\cdot 2^{9-1}=768+1536=2304,$$
$$F(2)=1+2^9\cdot 12-2304=3841,$$
ve son olarak
$$E(2)=2^{19}-3841=520447.$$
Bu da uygulamalardaki kontrol değeriyle aynıdır.
Uygulama önce tüm \(d\le n\) değerleri için \(\mu(d)\) ve \(\varphi(d)\) üreten lineer bir sieve kurar. Ayrıca \(\left\lfloor n/d\right\rfloor\) bölüm tablosunu saklar; çünkü hem \(P(n)\) hem de \(H_n(L)\) bu değerlere sürekli ihtiyaç duyar.
Sonra \(P(n)\) değeri \(2\Phi\) modunda hesaplanır, \(B(n)=6\sum_{L\le n}\varphi(L)\) ifadesi \(M\) modunda biriktirilir ve \(2\)'nin modüler kuvvetleri hızlı üs alma ile hesaplanır. Ardından \(L=1,\dots,n\) üzerinde dolaşılır; her \(L\) için Möbius ağırlıklı iç toplamdan \(H_n(L)\) bulunur, \(\left(\frac{P(n)-1}{2}-H_n(L)\right)\bmod\Phi\) üs kalanı oluşturulur ve \(6\varphi(L)\cdot 2^{((P(n)-1)/2-H_n(L))\bmod\Phi}\) katkısı \(T(n)\)'ye eklenir.
Son aşamada \(2^{P(n)\bmod\Phi}\), \(B(n)\) ve \(T(n)\) kutu içindeki formülle birleştirilir. C++ ve Java sürümleri \(L\) üzerindeki dış döngüyü paralelleştirirken Python sürümü aynı matematiği sıralı biçimde yürütür.
Sieve ve bölüm tablosu ayrı ayrı \(O(n)\) zaman ve \(O(n)\) bellek gerektirir. Baskın maliyet \(H_n(L)\) için gereken çift toplamdan gelir: her \(L\) için iç döngü \(\left\lfloor n/L\right\rfloor\)'ye kadar gider, dolayısıyla toplam yineleme sayısı
$$\sum_{L=1}^{n}\left\lfloor\frac{n}{L}\right\rfloor=n\log n+O(n)$$
olur. Böylece toplam çalışma süresi \(O(n\log n)\) aritmetik işlem ve \(O(n)\) bellek düzeyindedir. Paralelleştirme yalnızca gerçek çalışma süresini azaltır; asimptotik sınırı değiştirmez.
Sea
$$M=11^8=214358881.$$
Hay que calcular \(E(n)\) módulo \(M\). Las implementaciones en C++, Python y Java no enumeran las mezclas de forma directa. En su lugar, reescriben el problema en términos de datos de retícula primitivos dentro del cubo \(\{0,\dots,n\}^3\), junto con una corrección estructurada construida a partir de perfiles de borde bidimensionales reducidos. Los valores de control incluidos en las implementaciones son
$$E(1)=103,\qquad E(2)=520447,\qquad E(10)=82608406,\qquad E(500)=13801403.$$
La solución gira alrededor de cuatro cantidades: el conteo de ternas primitivas \(P(n)\), un término de borde ponderado por \(\varphi\) llamado \(B(n)\), un conteo interior coprimo \(H_n(L)\) y una suma de corrección final \(T(n)\). Una vez conocidas esas piezas, \(E(n)\) sale de una sola identidad modular.
Definimos
$$\mathcal{P}_n=\left\{(a,b,c)\in \{0,\dots,n\}^3\setminus\{(0,0,0)\}:\gcd(a,b,c)=1\right\}.$$
Sea \(P(n)=|\mathcal{P}_n|\). Para un divisor fijo \(d\), el número de ternas cuyas tres coordenadas son divisibles por \(d\) es
$$\left(\left\lfloor\frac{n}{d}\right\rfloor+1\right)^3-1,$$
porque cada coordenada tiene \(\left\lfloor n/d\right\rfloor+1\) múltiplos de \(d\) entre \(0\) y \(n\), y la terna nula debe excluirse. Por inversión de Möbius se obtiene
$$P(n)=\sum_{d=1}^{n}\mu(d)\left(\left(\left\lfloor\frac{n}{d}\right\rfloor+1\right)^3-1\right).$$
Éste es el gran exponente que después aparece en \(2^{P(n)}\).
Para cada \(L\ge 1\), el número de pares positivos coprimos \((u,v)\) que cumplen
$$u+v=L$$
es exactamente \(\varphi(L)\). En efecto, elegir \(u\) con \(1\le u\le L\) y \(\gcd(u,L)=1\) determina de manera única \(v=L-u\), y entonces \(\gcd(u,v)=1\). En la suma de corrección final, cada par reducido recibe un factor de simetría igual a \(6\); por eso aparece el coeficiente \(6\varphi(L)\). El peso global de borde es así
$$B(n)=6\sum_{L=1}^{n}\varphi(L).$$
Fijemos \(L\). La cantidad interior usada por las implementaciones es
$$H_n(L)=1+\#\left\{(x,y)\in \mathbb{Z}_{>0}^2:\gcd(x,y)=1,\ x+Ly\le n\right\}.$$
El \(1\) inicial forma parte explícita de la fórmula. Para contar los pares coprimos, se aplica otra vez la inversión de Möbius. Si \(d\mid x\) y \(d\mid y\), escribimos \(x=da\), \(y=db\). Entonces
$$a+Lb\le \left\lfloor\frac{n}{d}\right\rfloor.$$
Definimos
$$q_{d,L}=\left\lfloor\frac{n}{Ld}\right\rfloor.$$
Para cada \(b=1,\dots,q_{d,L}\), la variable \(a\) tiene \(\left\lfloor n/d\right\rfloor-Lb\) opciones positivas. Sumando esas opciones se obtiene
$$\sum_{b=1}^{q_{d,L}}\left(\left\lfloor\frac{n}{d}\right\rfloor-Lb\right)=q_{d,L}\left\lfloor\frac{n}{d}\right\rfloor-L\frac{q_{d,L}(q_{d,L}+1)}{2}.$$
Por tanto,
$$H_n(L)=1+\sum_{d=1}^{\lfloor n/L\rfloor}\mu(d)\left(q_{d,L}\left\lfloor\frac{n}{d}\right\rfloor-L\frac{q_{d,L}(q_{d,L}+1)}{2}\right),\qquad q_{d,L}=\left\lfloor\frac{n}{Ld}\right\rfloor.$$
Esta expresión coincide término a término con el bucle interior de la implementación.
Sea
$$\Phi=\varphi(M)=\varphi(11^8)=10\cdot 11^7=194871710.$$
Como \(\gcd(2,M)=1\), el teorema de Euler implica que cualquier potencia de \(2\) módulo \(M\) sólo depende del exponente módulo \(\Phi\). Sin embargo, los términos de corrección también necesitan
$$\frac{P(n)-1}{2},$$
de modo que la paridad debe conservarse antes de dividir entre \(2\). Por eso las implementaciones calculan primero \(P(n)\) módulo \(2\Phi\), restan \(1\), comprueban que el resultado sea par y sólo entonces reducen el semiexponente módulo \(\Phi\).
Con el semiexponente entendido módulo \(\Phi\), definimos
$$T(n)\equiv \sum_{L=1}^{n}6\varphi(L)\,2^{\left(\frac{P(n)-1}{2}-H_n(L)\right)\bmod \Phi}\pmod{M}.$$
Entonces el término común de corrección es
$$F(n)\equiv 1+2^{\left(\frac{P(n)-1}{2}\right)\bmod \Phi}B(n)-T(n)\pmod{M},$$
y el valor buscado es
$$\boxed{E(n)\equiv 2^{P(n)\bmod \Phi}-F(n)\pmod{M}.}$$
Ésta es exactamente el álgebra usada por las implementaciones en C++, Python y Java.
Aquí \(\mu(1)=1\), \(\mu(2)=-1\) y \(\varphi(1)=\varphi(2)=1\). Primero,
$$P(2)=\left((2+1)^3-1\right)-\left((1+1)^3-1\right)=26-7=19.$$
Por tanto,
$$\frac{P(2)-1}{2}=9,\qquad B(2)=6(\varphi(1)+\varphi(2))=12.$$
Para \(L=1\),
$$H_2(1)=1+\left(2\cdot 2-\frac{2\cdot 3}{2}\right)-\left(1\cdot 1-\frac{1\cdot 2}{2}\right)=2.$$
Para \(L=2\),
$$H_2(2)=1+\left(1\cdot 2-2\cdot\frac{1\cdot 2}{2}\right)=1.$$
Entonces
$$T(2)=6\cdot 2^{9-2}+6\cdot 2^{9-1}=768+1536=2304,$$
$$F(2)=1+2^9\cdot 12-2304=3841,$$
y finalmente
$$E(2)=2^{19}-3841=520447.$$
Esto coincide con el valor de control de las implementaciones.
La implementación empieza con una criba lineal que produce \(\mu(d)\) y \(\varphi(d)\) para todo \(d\le n\). También almacena la tabla de cocientes \(\left\lfloor n/d\right\rfloor\), porque tanto \(P(n)\) como \(H_n(L)\) reutilizan constantemente esos valores.
Después evalúa \(P(n)\) módulo \(2\Phi\), acumula \(B(n)=6\sum_{L\le n}\varphi(L)\) módulo \(M\) y calcula potencias modulares de \(2\) mediante exponenciación rápida. Luego recorre \(L=1,\dots,n\), obtiene \(H_n(L)\) a partir de la suma interior ponderada por Möbius, forma el residuo exponencial \(\left(\frac{P(n)-1}{2}-H_n(L)\right)\bmod \Phi\) y añade la contribución \(6\varphi(L)\cdot 2^{((P(n)-1)/2-H_n(L))\bmod\Phi}\) a \(T(n)\).
Al final combina \(2^{P(n)\bmod\Phi}\), \(B(n)\) y \(T(n)\) mediante la fórmula enmarcada. Las versiones de C++ y Java paralelizan el bucle exterior sobre \(L\), mientras que la versión de Python evalúa la misma matemática de forma secuencial.
La criba y la tabla de cocientes cuestan \(O(n)\) tiempo y \(O(n)\) memoria. El coste dominante proviene de la doble suma usada para \(H_n(L)\): para cada \(L\), el bucle interior llega hasta \(\left\lfloor n/L\right\rfloor\), así que el número total de iteraciones es
$$\sum_{L=1}^{n}\left\lfloor\frac{n}{L}\right\rfloor=n\log n+O(n).$$
Por tanto, el tiempo total es \(O(n\log n)\) operaciones aritméticas con \(O(n)\) memoria. La paralelización mejora el tiempo real de ejecución, pero no la cota asintótica.
记
$$M=11^8=214358881.$$
题目要求计算 \(E(n)\) 对 \(M\) 取模后的值。C++、Python 和 Java 实现都没有直接暴力枚举所有“mixtures”。它们先把问题改写成 \(\{0,\dots,n\}^3\) 这个立方体中的本原格点统计，再加上一项由二维约化边界轮廓组成的结构化修正。实现中内置的校验点为
$$E(1)=103,\qquad E(2)=520447,\qquad E(10)=82608406,\qquad E(500)=13801403.$$
整个解法围绕四个量展开：本原三元组计数 \(P(n)\)、以欧拉函数加权的边界项 \(B(n)\)、内部互素计数 \(H_n(L)\)，以及最后的修正和 \(T(n)\)。只要把这四部分求出来，\(E(n)\) 就能通过一个模公式直接得到。
定义
$$\mathcal{P}_n=\left\{(a,b,c)\in \{0,\dots,n\}^3\setminus\{(0,0,0)\}:\gcd(a,b,c)=1\right\}.$$
记 \(P(n)=|\mathcal{P}_n|\)。对固定的除数 \(d\)，若三坐标都被 \(d\) 整除，那么这样的三元组个数是
$$\left(\left\lfloor\frac{n}{d}\right\rfloor+1\right)^3-1,$$
因为每个坐标在 \(0\) 到 \(n\) 之间都有 \(\left\lfloor n/d\right\rfloor+1\) 个 \(d\) 的倍数可选，但必须去掉全零三元组。于是利用 Möbius 反演可得
$$P(n)=\sum_{d=1}^{n}\mu(d)\left(\left(\left\lfloor\frac{n}{d}\right\rfloor+1\right)^3-1\right).$$
后面出现的 \(2^{P(n)}\) 中，那个巨大的指数正是这个 \(P(n)\)。
对每个 \(L\ge 1\)，满足
$$u+v=L$$
且 \((u,v)\) 互素的正整数对数量恰好等于 \(\varphi(L)\)。理由是：只要选定一个 \(1\le u\le L\) 且 \(\gcd(u,L)=1\) 的 \(u\)，就唯一确定了 \(v=L-u\)，同时自动有 \(\gcd(u,v)=1\)。在最终修正和里，每个约化后的整数对还要乘上一个六重对称因子，所以系数变成 \(6\varphi(L)\)。因此全局边界权重为
$$B(n)=6\sum_{L=1}^{n}\varphi(L).$$
固定 \(L\)。实现中使用的内部量是
$$H_n(L)=1+\#\left\{(x,y)\in \mathbb{Z}_{>0}^2:\gcd(x,y)=1,\ x+Ly\le n\right\}.$$
开头那个 \(1\) 是公式显式带上的常数项。为了统计互素整数对，再次使用 Möbius 反演。若 \(d\mid x\) 且 \(d\mid y\)，写成 \(x=da\)、\(y=db\)，则有
$$a+Lb\le \left\lfloor\frac{n}{d}\right\rfloor.$$
记
$$q_{d,L}=\left\lfloor\frac{n}{Ld}\right\rfloor.$$
对每个 \(b=1,\dots,q_{d,L}\)，变量 \(a\) 有 \(\left\lfloor n/d\right\rfloor-Lb\) 个正整数选择，因此这些选择总数为
$$\sum_{b=1}^{q_{d,L}}\left(\left\lfloor\frac{n}{d}\right\rfloor-Lb\right)=q_{d,L}\left\lfloor\frac{n}{d}\right\rfloor-L\frac{q_{d,L}(q_{d,L}+1)}{2}.$$
于是得到
$$H_n(L)=1+\sum_{d=1}^{\lfloor n/L\rfloor}\mu(d)\left(q_{d,L}\left\lfloor\frac{n}{d}\right\rfloor-L\frac{q_{d,L}(q_{d,L}+1)}{2}\right),\qquad q_{d,L}=\left\lfloor\frac{n}{Ld}\right\rfloor.$$
这与实现里的内部循环逐项完全一致。
令
$$\Phi=\varphi(M)=\varphi(11^8)=10\cdot 11^7=194871710.$$
由于 \(\gcd(2,M)=1\)，根据 Euler 定理，模 \(M\) 下任意 \(2\) 的幂只取决于指数对 \(\Phi\) 的余数。但修正项还需要
$$\frac{P(n)-1}{2},$$
也就是说，在除以 \(2\) 之前必须先保留奇偶性信息。因此实现先计算 \(P(n)\bmod 2\Phi\)，然后减去 \(1\)，确认结果为偶数，再把半指数降到 \(\Phi\) 模下。
在把半指数理解为模 \(\Phi\) 的量之后，定义
$$T(n)\equiv \sum_{L=1}^{n}6\varphi(L)\,2^{\left(\frac{P(n)-1}{2}-H_n(L)\right)\bmod \Phi}\pmod{M}.$$
然后公共修正项为
$$F(n)\equiv 1+2^{\left(\frac{P(n)-1}{2}\right)\bmod \Phi}B(n)-T(n)\pmod{M},$$
所求结果就是
$$\boxed{E(n)\equiv 2^{P(n)\bmod \Phi}-F(n)\pmod{M}.}$$
这正是 C++、Python 和 Java 实现所执行的代数过程。
此时 \(\mu(1)=1\)、\(\mu(2)=-1\)，并且 \(\varphi(1)=\varphi(2)=1\)。先算
$$P(2)=\left((2+1)^3-1\right)-\left((1+1)^3-1\right)=26-7=19.$$
因此
$$\frac{P(2)-1}{2}=9,\qquad B(2)=6(\varphi(1)+\varphi(2))=12.$$
当 \(L=1\) 时，
$$H_2(1)=1+\left(2\cdot 2-\frac{2\cdot 3}{2}\right)-\left(1\cdot 1-\frac{1\cdot 2}{2}\right)=2.$$
当 \(L=2\) 时，
$$H_2(2)=1+\left(1\cdot 2-2\cdot\frac{1\cdot 2}{2}\right)=1.$$
于是
$$T(2)=6\cdot 2^{9-2}+6\cdot 2^{9-1}=768+1536=2304,$$
$$F(2)=1+2^9\cdot 12-2304=3841,$$
最后得到
$$E(2)=2^{19}-3841=520447.$$
这与实现中的校验值完全一致。
实现首先建立线性筛，同时求出所有 \(d\le n\) 的 \(\mu(d)\) 和 \(\varphi(d)\)。另外还保存 \(\left\lfloor n/d\right\rfloor\) 的商表，因为 \(P(n)\) 与 \(H_n(L)\) 都会反复使用这些值。
接着，程序在模 \(2\Phi\) 下计算 \(P(n)\)，在模 \(M\) 下累加 \(B(n)=6\sum_{L\le n}\varphi(L)\)，并用快速幂计算模 \(M\) 的 \(2\) 次幂。然后对 \(L=1,\dots,n\) 逐个处理：先从 Möbius 加权的内部求和得到 \(H_n(L)\)，再形成指数余数 \(\left(\frac{P(n)-1}{2}-H_n(L)\right)\bmod \Phi\)，并把贡献 \(6\varphi(L)\cdot 2^{((P(n)-1)/2-H_n(L))\bmod\Phi}\) 加进 \(T(n)\)。
最后再把 \(2^{P(n)\bmod\Phi}\)、\(B(n)\) 与 \(T(n)\) 按照上面的盒装公式合并。C++ 与 Java 版本把外层 \(L\) 循环并行化，而 Python 版本则顺序执行同一套数学过程。
线性筛和商表都需要 \(O(n)\) 时间与 \(O(n)\) 空间。主导开销来自 \(H_n(L)\) 的双重求和：对每个 \(L\)，内部循环运行到 \(\left\lfloor n/L\right\rfloor\)，所以总迭代次数为
$$\sum_{L=1}^{n}\left\lfloor\frac{n}{L}\right\rfloor=n\log n+O(n).$$
因此总时间复杂度是 \(O(n\log n)\) 次算术操作，空间复杂度是 \(O(n)\)。并行化只能改善实际运行时间，不会改变渐近复杂度。
Положим
$$M=11^8=214358881.$$
Требуется вычислить \(E(n)\) по модулю \(M\). Реализации на C++, Python и Java не перебирают смеси напрямую. Вместо этого задача переписывается через подсчёт примитивных решёточных данных в кубе \(\{0,\dots,n\}^3\) и через структурированную поправку, построенную из приведённых двумерных граничных профилей. Контрольные значения, зашитые в реализации, таковы:
$$E(1)=103,\qquad E(2)=520447,\qquad E(10)=82608406,\qquad E(500)=13801403.$$
Решение организовано вокруг четырёх величин: числа примитивных троек \(P(n)\), граничного члена \(B(n)\) с весами \(\varphi\), внутреннего взаимно простого счёта \(H_n(L)\) и итоговой суммы поправок \(T(n)\). Когда эти четыре величины найдены, \(E(n)\) выражается одной модульной формулой.
Определим
$$\mathcal{P}_n=\left\{(a,b,c)\in \{0,\dots,n\}^3\setminus\{(0,0,0)\}:\gcd(a,b,c)=1\right\}.$$
Пусть \(P(n)=|\mathcal{P}_n|\). Для фиксированного делителя \(d\) число троек, у которых все три координаты делятся на \(d\), равно
$$\left(\left\lfloor\frac{n}{d}\right\rfloor+1\right)^3-1,$$
потому что каждая координата имеет \(\left\lfloor n/d\right\rfloor+1\) допустимых кратных между \(0\) и \(n\), а нулевую тройку нужно убрать. Тогда обращение Мёбиуса даёт
$$P(n)=\sum_{d=1}^{n}\mu(d)\left(\left(\left\lfloor\frac{n}{d}\right\rfloor+1\right)^3-1\right).$$
Именно это число затем становится большим показателем в \(2^{P(n)}\).
Для каждого \(L\ge 1\) число взаимно простых положительных пар \((u,v)\), удовлетворяющих
$$u+v=L,$$
равно точно \(\varphi(L)\). Действительно, выбор \(u\) при \(1\le u\le L\) и \(\gcd(u,L)=1\) однозначно определяет \(v=L-u\), после чего автоматически \(\gcd(u,v)=1\). В итоговой поправке каждая приведённая пара входит с шестикратным симметричным весом, откуда и появляется коэффициент \(6\varphi(L)\). Поэтому глобальный граничный вес равен
$$B(n)=6\sum_{L=1}^{n}\varphi(L).$$
Зафиксируем \(L\). Внутренняя величина, используемая реализациями, имеет вид
$$H_n(L)=1+\#\left\{(x,y)\in \mathbb{Z}_{>0}^2:\gcd(x,y)=1,\ x+Ly\le n\right\}.$$
Начальная единица входит в формулу явно. Чтобы посчитать взаимно простые пары, снова применяется обращение Мёбиуса. Если \(d\mid x\) и \(d\mid y\), то пишем \(x=da\), \(y=db\). Тогда
$$a+Lb\le \left\lfloor\frac{n}{d}\right\rfloor.$$
Обозначим
$$q_{d,L}=\left\lfloor\frac{n}{Ld}\right\rfloor.$$
Для каждого \(b=1,\dots,q_{d,L}\) переменная \(a\) имеет \(\left\lfloor n/d\right\rfloor-Lb\) положительных значений. Сумма этих вариантов равна
$$\sum_{b=1}^{q_{d,L}}\left(\left\lfloor\frac{n}{d}\right\rfloor-Lb\right)=q_{d,L}\left\lfloor\frac{n}{d}\right\rfloor-L\frac{q_{d,L}(q_{d,L}+1)}{2}.$$
Следовательно,
$$H_n(L)=1+\sum_{d=1}^{\lfloor n/L\rfloor}\mu(d)\left(q_{d,L}\left\lfloor\frac{n}{d}\right\rfloor-L\frac{q_{d,L}(q_{d,L}+1)}{2}\right),\qquad q_{d,L}=\left\lfloor\frac{n}{Ld}\right\rfloor.$$
Это в точности та же внутренняя сумма, что и в коде.
Пусть
$$\Phi=\varphi(M)=\varphi(11^8)=10\cdot 11^7=194871710.$$
Так как \(\gcd(2,M)=1\), по теореме Эйлера любая степень \(2\) по модулю \(M\) зависит только от показателя по модулю \(\Phi\). Но в поправке требуется ещё и величина
$$\frac{P(n)-1}{2},$$
поэтому перед делением на \(2\) нужно сохранить чётность. Именно поэтому реализации сначала вычисляют \(P(n)\) по модулю \(2\Phi\), затем вычитают \(1\), проверяют чётность результата и только после этого берут половину по модулю \(\Phi\).
Понимая половинный показатель по модулю \(\Phi\), определяем
$$T(n)\equiv \sum_{L=1}^{n}6\varphi(L)\,2^{\left(\frac{P(n)-1}{2}-H_n(L)\right)\bmod \Phi}\pmod{M}.$$
Тогда общий корректирующий член равен
$$F(n)\equiv 1+2^{\left(\frac{P(n)-1}{2}\right)\bmod \Phi}B(n)-T(n)\pmod{M},$$
а искомое значение имеет вид
$$\boxed{E(n)\equiv 2^{P(n)\bmod \Phi}-F(n)\pmod{M}.}$$
Именно эту алгебру используют реализации на C++, Python и Java.
Здесь \(\mu(1)=1\), \(\mu(2)=-1\), а \(\varphi(1)=\varphi(2)=1\). Сначала
$$P(2)=\left((2+1)^3-1\right)-\left((1+1)^3-1\right)=26-7=19.$$
Значит,
$$\frac{P(2)-1}{2}=9,\qquad B(2)=6(\varphi(1)+\varphi(2))=12.$$
Для \(L=1\)
$$H_2(1)=1+\left(2\cdot 2-\frac{2\cdot 3}{2}\right)-\left(1\cdot 1-\frac{1\cdot 2}{2}\right)=2.$$
Для \(L=2\)
$$H_2(2)=1+\left(1\cdot 2-2\cdot\frac{1\cdot 2}{2}\right)=1.$$
Следовательно,
$$T(2)=6\cdot 2^{9-2}+6\cdot 2^{9-1}=768+1536=2304,$$
$$F(2)=1+2^9\cdot 12-2304=3841,$$
и потому
$$E(2)=2^{19}-3841=520447.$$
Это совпадает с контрольным значением из реализаций.
Сначала строится линейное решето, которое вычисляет \(\mu(d)\) и \(\varphi(d)\) для всех \(d\le n\). Дополнительно сохраняется таблица частных \(\left\lfloor n/d\right\rfloor\), поскольку и \(P(n)\), и \(H_n(L)\) многократно используют эти значения.
Затем программа вычисляет \(P(n)\) по модулю \(2\Phi\), накапливает \(B(n)=6\sum_{L\le n}\varphi(L)\) по модулю \(M\) и считает степени двойки быстрым возведением в степень. После этого идёт цикл по \(L=1,\dots,n\): для каждого \(L\) по внутренней сумме с весами Мёбиуса вычисляется \(H_n(L)\), формируется показатель \(\left(\frac{P(n)-1}{2}-H_n(L)\right)\bmod \Phi\), и вклад \(6\varphi(L)\cdot 2^{((P(n)-1)/2-H_n(L))\bmod\Phi}\) добавляется к \(T(n)\).
В конце \(2^{P(n)\bmod\Phi}\), \(B(n)\) и \(T(n)\) объединяются по формуле в рамке. Реализации на C++ и Java распараллеливают внешний цикл по \(L\), а реализация на Python выполняет ту же математику последовательно.
Решето и таблица частных требуют по \(O(n)\) времени и \(O(n)\) памяти. Основная стоимость сосредоточена в двойной сумме для \(H_n(L)\): для каждого \(L\) внутренний цикл идёт до \(\left\lfloor n/L\right\rfloor\), так что суммарное число итераций равно
$$\sum_{L=1}^{n}\left\lfloor\frac{n}{L}\right\rfloor=n\log n+O(n).$$
Следовательно, общая трудоёмкость равна \(O(n\log n)\) арифметических операций при \(O(n)\) памяти. Параллельность улучшает только реальное время выполнения, но не асимптотику.
لنكتب
$$M=11^8=214358881.$$
المطلوب هو حساب \(E(n)\) بترديد \(M\). لا تقوم تطبيقات C++ وPython وJava بعدّ جميع الخلطات مباشرة. بدلًا من ذلك، تعيد صياغة المسألة بوصفها عدًا لبيانات شبكية بدائية داخل المكعب \(\{0,\dots,n\}^3\)، ثم تضيف حدًا تصحيحيًا منظمًا مبنيًا من ملفات حدية ثنائية الأبعاد في صورة مختزلة. نقاط التحقق المضمّنة في التنفيذات هي
$$E(1)=103,\qquad E(2)=520447,\qquad E(10)=82608406,\qquad E(500)=13801403.$$
الحل مبني حول أربع كميات: عدد الثلاثيات البدائية \(P(n)\)، وحد حدّي موزون بدالة أويلر \(B(n)\)، وعدّ داخلي للأزواج المتباينة أوليًا \(H_n(L)\)، ثم مجموع تصحيحي نهائي \(T(n)\). متى حُسبت هذه العناصر الأربع أمكن استخراج \(E(n)\) من هوية معيارية واحدة.
نعرّف
$$\mathcal{P}_n=\left\{(a,b,c)\in \{0,\dots,n\}^3\setminus\{(0,0,0)\}:\gcd(a,b,c)=1\right\}.$$
وليكن \(P(n)=|\mathcal{P}_n|\). إذا ثبتنا قاسمًا \(d\)، فإن عدد الثلاثيات التي تكون إحداثياتها الثلاثة جميعًا قابلة للقسمة على \(d\) يساوي
$$\left(\left\lfloor\frac{n}{d}\right\rfloor+1\right)^3-1,$$
لأن كل إحداثي يملك \(\left\lfloor n/d\right\rfloor+1\) مضاعفًا مناسبًا بين \(0\) و\(n\)، مع استبعاد الثلاثية الصفرية. ومن ثم تعطي صيغة Möbius العكسية
$$P(n)=\sum_{d=1}^{n}\mu(d)\left(\left(\left\lfloor\frac{n}{d}\right\rfloor+1\right)^3-1\right).$$
وهذا هو الأس الكبير الذي يظهر لاحقًا في الحد \(2^{P(n)}\).
لكل \(L\ge 1\)، فإن عدد الأزواج الموجبة المتباينة أوليًا \((u,v)\) التي تحقق
$$u+v=L$$
هو بالضبط \(\varphi(L)\). والسبب أن اختيار \(u\) بحيث \(1\le u\le L\) و\(\gcd(u,L)=1\) يحدد \(v=L-u\) بصورة وحيدة، وعندها يكون \(\gcd(u,v)=1\). في مجموع التصحيح النهائي تُعامل كل زوجة مختزلة مع عامل تماثل مقداره ستة، ومن هنا يأتي المعامل \(6\varphi(L)\). لذا فإن الوزن الحدّي الكلي هو
$$B(n)=6\sum_{L=1}^{n}\varphi(L).$$
ثبّت قيمة \(L\). الكمية الداخلية المستخدمة في التنفيذات هي
$$H_n(L)=1+\#\left\{(x,y)\in \mathbb{Z}_{>0}^2:\gcd(x,y)=1,\ x+Ly\le n\right\}.$$
إن الواحد في البداية جزء صريح من الصيغة. ولعدّ الأزواج المتباينة أوليًا نستخدم Möbius مرة أخرى. إذا كان \(d\mid x\) و\(d\mid y\)، نكتب \(x=da\) و\(y=db\). عندئذ
$$a+Lb\le \left\lfloor\frac{n}{d}\right\rfloor.$$
ولنعرّف
$$q_{d,L}=\left\lfloor\frac{n}{Ld}\right\rfloor.$$
لكل \(b=1,\dots,q_{d,L}\)، يملك المتغير \(a\) عددًا من الخيارات الموجبة يساوي \(\left\lfloor n/d\right\rfloor-Lb\). وبجمع هذه الخيارات نحصل على
$$\sum_{b=1}^{q_{d,L}}\left(\left\lfloor\frac{n}{d}\right\rfloor-Lb\right)=q_{d,L}\left\lfloor\frac{n}{d}\right\rfloor-L\frac{q_{d,L}(q_{d,L}+1)}{2}.$$
ومن ثم
$$H_n(L)=1+\sum_{d=1}^{\lfloor n/L\rfloor}\mu(d)\left(q_{d,L}\left\lfloor\frac{n}{d}\right\rfloor-L\frac{q_{d,L}(q_{d,L}+1)}{2}\right),\qquad q_{d,L}=\left\lfloor\frac{n}{Ld}\right\rfloor.$$
وهذا يطابق الحلقة الداخلية في التنفيذ بندًا بندًا.
لنضع
$$\Phi=\varphi(M)=\varphi(11^8)=10\cdot 11^7=194871710.$$
بما أن \(\gcd(2,M)=1\)، فإن مبرهنة أويلر تقول إن أي قوة للعدد \(2\) بترديد \(M\) تعتمد فقط على الأس بترديد \(\Phi\). لكن حدود التصحيح تحتاج أيضًا إلى
$$\frac{P(n)-1}{2},$$
ولذلك يجب الحفاظ على معلومة الزوجية قبل القسمة على \(2\). لهذا السبب تحسب التنفيذات أولًا \(P(n)\) بترديد \(2\Phi\)، ثم تطرح \(1\)، ثم تتحقق من أن الناتج زوجي، وبعد ذلك فقط تختزل نصف الأس بترديد \(\Phi\).
بعد فهم نصف الأس بترديد \(\Phi\)، نعرّف
$$T(n)\equiv \sum_{L=1}^{n}6\varphi(L)\,2^{\left(\frac{P(n)-1}{2}-H_n(L)\right)\bmod \Phi}\pmod{M}.$$
وعندئذ يكون حد التصحيح المشترك
$$F(n)\equiv 1+2^{\left(\frac{P(n)-1}{2}\right)\bmod \Phi}B(n)-T(n)\pmod{M},$$
والقيمة المطلوبة هي
$$\boxed{E(n)\equiv 2^{P(n)\bmod \Phi}-F(n)\pmod{M}.}$$
وهذا هو البناء الجبري نفسه الذي تستخدمه تطبيقات C++ وPython وJava.
هنا \(\mu(1)=1\) و\(\mu(2)=-1\)، كما أن \(\varphi(1)=\varphi(2)=1\). أولًا
$$P(2)=\left((2+1)^3-1\right)-\left((1+1)^3-1\right)=26-7=19.$$
إذًا
$$\frac{P(2)-1}{2}=9,\qquad B(2)=6(\varphi(1)+\varphi(2))=12.$$
عندما \(L=1\) نحصل على
$$H_2(1)=1+\left(2\cdot 2-\frac{2\cdot 3}{2}\right)-\left(1\cdot 1-\frac{1\cdot 2}{2}\right)=2.$$
وعندما \(L=2\) يكون
$$H_2(2)=1+\left(1\cdot 2-2\cdot\frac{1\cdot 2}{2}\right)=1.$$
ومن ثم
$$T(2)=6\cdot 2^{9-2}+6\cdot 2^{9-1}=768+1536=2304,$$
$$F(2)=1+2^9\cdot 12-2304=3841,$$
وأخيرًا
$$E(2)=2^{19}-3841=520447.$$
وهذا يطابق قيمة التحقق الموجودة في التنفيذات.
تبدأ الشيفرة ببناء غربال خطي ينتج \(\mu(d)\) و\(\varphi(d)\) لكل \(d\le n\). كما تحتفظ بجدول القيم \(\left\lfloor n/d\right\rfloor\)، لأن \(P(n)\) و\(H_n(L)\) يعيدان استخدام هذه الكسور الصحيحة مرات كثيرة.
بعد ذلك تُحسب \(P(n)\) بترديد \(2\Phi\)، ويُجمَّع \(B(n)=6\sum_{L\le n}\varphi(L)\) بترديد \(M\)، وتُحسب قوى العدد \(2\) بالرفع السريع. ثم تدور حلقة على \(L=1,\dots,n\): يُحسب \(H_n(L)\) من المجموع الداخلي الموزون بـ Möbius، ثم يُشكَّل الأس \(\left(\frac{P(n)-1}{2}-H_n(L)\right)\bmod \Phi\)، ثم تُضاف المساهمة \(6\varphi(L)\cdot 2^{((P(n)-1)/2-H_n(L))\bmod\Phi}\) إلى \(T(n)\).
في النهاية تُدمج \(2^{P(n)\bmod\Phi}\) و\(B(n)\) و\(T(n)\) وفق الصيغة المؤطرة أعلاه. وتقوم نسختا C++ وJava بموازاة الحلقة الخارجية على \(L\)، بينما تنفذ نسخة Python الرياضيات نفسها بصورة تسلسلية.
يحتاج الغربال وجدول القيم إلى \(O(n)\) زمنًا و\(O(n)\) ذاكرة. أما الكلفة الغالبة فتأتي من المجموع المزدوج الخاص بـ \(H_n(L)\): لكل \(L\) تمتد الحلقة الداخلية حتى \(\left\lfloor n/L\right\rfloor\)، ومن ثم يكون مجموع عدد التكرارات
$$\sum_{L=1}^{n}\left\lfloor\frac{n}{L}\right\rfloor=n\log n+O(n).$$
إذن الزمن الكلي هو \(O(n\log n)\) من العمليات الحسابية، مع ذاكرة \(O(n)\). الموازاة تحسن زمن التنفيذ الفعلي فقط، ولا تغيّر الرتبة التقاربية.