Define
$$N=10^{18},\qquad K=10^9,\qquad B=\binom{N}{K}.$$
For every prime \(p\) with \(1000 \lt p \lt 5000\), we need the residue \(B \bmod p\). Then, for every triple of distinct primes \(p \lt q \lt r\) in that interval, we reconstruct the unique number \(x_{pqr}\) with
$$0 \le x_{pqr} \lt pqr,\qquad x_{pqr}\equiv B \pmod{p},\qquad x_{pqr}\equiv B \pmod{q},\qquad x_{pqr}\equiv B \pmod{r}.$$
The required answer is the sum of all these reconstructed values. The sieve used in the code finds exactly \(501\) primes in the interval, so the outer summation runs over
$$\binom{501}{3}=20833250$$
prime triples.
For a fixed prime \(p\), write \(N\) and \(K\) in base \(p\):
$$N=\sum_{t=0}^{s} N_t p^t,\qquad K=\sum_{t=0}^{s} K_t p^t,\qquad 0 \le N_t,K_t \lt p.$$
Lucas's theorem gives the congruence
$$\binom{N}{K}\equiv \prod_{t=0}^{s}\binom{N_t}{K_t}\pmod{p}.$$
If some digit satisfies \(K_t \gt N_t\), then \(\binom{N_t}{K_t}=0\), so the whole product is \(0\pmod p\). This is exactly why the implementation stops immediately and returns zero in that case.
The interval \(1000 \lt p \lt 5000\) makes the digit loop very short. Because \(p \gt 1000\), the number \(10^{18}\) has at most six base-\(p\) digits, and \(10^9\) has at most three. So each residue \(B \bmod p\) is computed from only a handful of digit-level binomials.
For \(0 \le b \le a \lt p\), we evaluate
$$\binom{a}{b}=\frac{a!}{b!(a-b)!}\pmod{p}.$$
Because \(p\) is prime and \(a,b \lt p\), the denominator is invertible modulo \(p\). The code therefore precomputes
$$\text{fact}[i]=i!\pmod p,\qquad \text{invFact}[i]=(i!)^{-1}\pmod p,$$
and then uses
$$\binom{a}{b}\equiv \text{fact}[a]\cdot \text{invFact}[b]\cdot \text{invFact}[a-b]\pmod p.$$
The modular inverses come from Fermat's little theorem:
$$x^{p-1}\equiv 1\pmod p\quad\Longrightarrow\quad x^{-1}\equiv x^{p-2}\pmod p,$$
which is why all three implementations contain a fast modular exponentiation routine.
The C++ code verifies Lucas's theorem on the smaller example \(\binom{30}{12}\). Take \(p=11\). In base \(11\),
$$30=2\cdot 11+8,\qquad 12=1\cdot 11+1.$$
Lucas gives
$$\binom{30}{12}\equiv \binom{8}{1}\binom{2}{1}=8\cdot 2=16\equiv 5\pmod{11}.$$
The exact value is \(\binom{30}{12}=86493225\), and indeed \(86493225 \equiv 5 \pmod{11}\). The checkpoint in the local source repeats this comparison for several primes.
Suppose we already know
$$a\equiv B\pmod p,\qquad b\equiv B\pmod q,\qquad c\equiv B\pmod r,$$
with \(p,q,r\) distinct primes. Since these moduli are pairwise coprime, the Chinese Remainder Theorem guarantees a unique residue modulo \(pqr\).
The code uses a Garner-style two-stage reconstruction. First write
$$x_{pq}=a+p\,t.$$
To make this also congruent to \(b\pmod q\), we need
$$a+p\,t\equiv b\pmod q,$$
so
$$t\equiv (b-a)\,p^{-1}\pmod q.$$
Hence
$$x_{pq}=a+p\left((b-a)\,p^{-1}\bmod q\right),\qquad 0 \le x_{pq} \lt pq.$$
Now lift once more by writing
$$x=x_{pq}+pq\,u.$$
Imposing \(x\equiv c\pmod r\) yields
$$u\equiv (c-x_{pq})(pq)^{-1}\pmod r,$$
and therefore
$$x=x_{pq}+pq\left((c-x_{pq})(pq)^{-1}\bmod r\right),\qquad 0 \le x \lt pqr.$$
This final \(x\) is the value denoted \(x_{pqr}\) in the statement.
The implementations store every pairwise inverse
$$\text{inv}[i][j]\equiv p_i^{-1}\pmod{p_j}.$$
Then the inverse of a product is obtained for free:
$$ (pq)^{-1}\equiv p^{-1}q^{-1}\pmod r.$$
So once the table of pairwise inverses has been built, each triple reconstruction uses only a few modular additions and multiplications. There is no extended Euclidean algorithm inside the cubic loop.
The local C++ checkpoints also test the CRT stage on \(\binom{20}{8}=125970\) and the tiny prime set \(\{7,11,13,17\}\). For the triple \((7,11,13)\), the exact residues are
$$125970\equiv 5\pmod 7,\qquad 125970\equiv 9\pmod{11},\qquad 125970\equiv 0\pmod{13}.$$
First combine \(7\) and \(11\):
$$t\equiv (9-5)\cdot 7^{-1}\equiv 4\cdot 8\equiv 10\pmod{11},$$
so
$$x_{7,11}=5+7\cdot 10=75.$$
Next combine with \(13\): since \(75\equiv 10\pmod{13}\) and \(77\equiv 12\pmod{13}\),
$$u\equiv (0-10)\cdot 12^{-1}\equiv 3\cdot 12\equiv 10\pmod{13},$$
which gives
$$x_{7,11,13}=75+77\cdot 10=845.$$
This matches the direct reduction \(125970\bmod(7\cdot 11\cdot 13)=845\). Summing the four tiny triples \((7,11,13)\), \((7,11,17)\), \((7,13,17)\), and \((11,13,17)\) gives the checkpoint total \(3803\) used by the source.
If \(\mathcal P\) denotes the set of primes between \(1000\) and \(5000\), the target quantity is
$$\boxed{S=\sum_{p \lt q \lt r,\; p,q,r\in\mathcal P} x_{pqr}.}$$
Running the supplied implementations yields
$$S=162619462356610313.$$
The three language versions follow the same structure. They first generate the prime list with a sieve, then compute every Lucas residue \(B \bmod p\), then precompute all pairwise inverses \(p_i^{-1}\bmod p_j\), and finally iterate over all prime triples. The C++ version accumulates the sum in unsigned __int128, the Java version uses BigInteger, and the Python version relies on Python's built-in arbitrary-precision integers.
Let \(m=501\) be the number of primes in the interval. Sieve construction up to \(5000\) is negligible, roughly \(O(5000\log\log 5000)\). The Lucas preprocessing costs
$$\sum_{p\in\mathcal P} O(p+\log_p N),$$
because each prime builds factorial and inverse-factorial tables of size \(p\), while the digit loop is tiny. Precomputing pairwise inverses costs \(O(m^2)\) time and memory. The dominant stage is the triple loop, which performs
$$\binom{m}{3}=\binom{501}{3}=20833250$$
constant-time CRT merges. Therefore the overall running time is \(O(m^3)\), and the memory usage is \(O(m^2)\) because of the inverse matrix.
Wir definieren
$$N=10^{18},\qquad K=10^9,\qquad B=\binom{N}{K}.$$
Für jedes Prim \(p\) mit \(1000 \lt p \lt 5000\) wird zunächst \(B \bmod p\) berechnet. Danach wird für jedes Primtripel \(p \lt q \lt r\) der eindeutige Rest \(x_{pqr}\) mit
$$0 \le x_{pqr} \lt pqr,\qquad x_{pqr}\equiv B \pmod p,\qquad x_{pqr}\equiv B \pmod q,\qquad x_{pqr}\equiv B \pmod r$$
rekonstruiert und schließlich aufsummiert. Das im Code verwendete Sieb findet genau \(501\) Primzahlen im Intervall, also insgesamt
$$\binom{501}{3}=20833250$$
Tripel.
Für ein festes Prim \(p\) schreiben wir \(N\) und \(K\) in Basis \(p\):
$$N=\sum_{t=0}^{s} N_t p^t,\qquad K=\sum_{t=0}^{s} K_t p^t,\qquad 0 \le N_t,K_t \lt p.$$
Dann gilt nach Lucas
$$\binom{N}{K}\equiv \prod_{t=0}^{s}\binom{N_t}{K_t}\pmod p.$$
Sobald für irgendeine Stelle \(K_t \gt N_t\) ist, wird der entsprechende kleine Binomialkoeffizient null, also auch das gesamte Produkt modulo \(p\). Genau dieses Verhalten sieht man in der Implementierung.
Da alle betrachteten Primzahlen größer als \(1000\) sind, hat \(10^{18}\) in Basis \(p\) höchstens sechs Stellen und \(10^9\) höchstens drei. Deshalb besteht jede Lucas-Berechnung nur aus wenigen Ziffernschritten.
Für \(0 \le b \le a \lt p\) verwenden wir
$$\binom{a}{b}=\frac{a!}{b!(a-b)!}\pmod p.$$
Weil \(p\) prim ist und alle Faktoren kleiner als \(p\) sind, ist der Nenner modulo \(p\) invertierbar. Daher berechnet der Code
$$\text{fact}[i]=i!\pmod p,\qquad \text{invFact}[i]=(i!)^{-1}\pmod p,$$
und setzt dann
$$\binom{a}{b}\equiv \text{fact}[a]\cdot \text{invFact}[b]\cdot \text{invFact}[a-b]\pmod p.$$
Die Inversen stammen aus dem kleinen Satz von Fermat, denn
$$x^{p-1}\equiv 1\pmod p\quad\Longrightarrow\quad x^{-1}\equiv x^{p-2}\pmod p.$$
Darum enthalten alle drei Sprachversionen eine schnelle modulare Potenzfunktion.
Die C++-Datei prüft Lucas an \(\binom{30}{12}\). Für \(p=11\) gilt in Basis \(11\)
$$30=2\cdot 11+8,\qquad 12=1\cdot 11+1.$$
Daraus folgt
$$\binom{30}{12}\equiv \binom{8}{1}\binom{2}{1}=16\equiv 5\pmod{11}.$$
Der exakte Wert \(\binom{30}{12}=86493225\) liefert ebenfalls Rest \(5\) modulo \(11\). Im lokalen Test wird derselbe Vergleich noch für weitere Primzahlen wiederholt.
Seien
$$a\equiv B\pmod p,\qquad b\equiv B\pmod q,\qquad c\equiv B\pmod r$$
die schon berechneten Reste. Da \(p,q,r\) verschiedene Primzahlen sind, sind die Moduli paarweise teilerfremd; also existiert genau eine Lösung modulo \(pqr\).
Der Code rekonstruiert diese Lösung in zwei Stufen. Zuerst setzen wir
$$x_{pq}=a+p\,t.$$
Aus der Bedingung \(x_{pq}\equiv b\pmod q\) folgt
$$t\equiv (b-a)\,p^{-1}\pmod q,$$
also
$$x_{pq}=a+p\left((b-a)\,p^{-1}\bmod q\right),\qquad 0 \le x_{pq} \lt pq.$$
Dann heben wir auf das dritte Modul an und schreiben
$$x=x_{pq}+pq\,u.$$
Die Forderung \(x\equiv c\pmod r\) ergibt
$$u\equiv (c-x_{pq})(pq)^{-1}\pmod r,$$
somit
$$x=x_{pq}+pq\left((c-x_{pq})(pq)^{-1}\bmod r\right),\qquad 0 \le x \lt pqr.$$
Genau dieses \(x\) ist der gesuchte Wert \(x_{pqr}\).
Die Programme speichern alle paarweisen Inversen
$$\text{inv}[i][j]\equiv p_i^{-1}\pmod{p_j}.$$
Dann gilt sofort
$$ (pq)^{-1}\equiv p^{-1}q^{-1}\pmod r.$$
Dadurch braucht die innere Dreifachschleife nur noch wenige modulare Additionen und Multiplikationen. Ein neuer Euklidischer Algorithmus pro Tripel wäre unnötig teuer.
Die lokale C++-Implementierung testet die CRT-Stufe außerdem an \(\binom{20}{8}=125970\) und den kleinen Primzahlen \(\{7,11,13,17\}\). Für das Tripel \((7,11,13)\) sind die Reste
$$125970\equiv 5\pmod 7,\qquad 125970\equiv 9\pmod{11},\qquad 125970\equiv 0\pmod{13}.$$
Zunächst kombiniert man \(7\) und \(11\):
$$t\equiv (9-5)\cdot 7^{-1}\equiv 4\cdot 8\equiv 10\pmod{11},$$
woraus
$$x_{7,11}=5+7\cdot 10=75$$
folgt. Dann mit \(13\): weil \(75\equiv 10\pmod{13}\) und \(77\equiv 12\pmod{13}\), erhält man
$$u\equiv (0-10)\cdot 12^{-1}\equiv 3\cdot 12\equiv 10\pmod{13},$$
also
$$x_{7,11,13}=75+77\cdot 10=845.$$
Das stimmt genau mit \(125970 \bmod (7\cdot 11\cdot 13)=845\) überein. Die Summe über die vier kleinen Tripel ergibt den Testwert \(3803\), den der Quellcode verifiziert.
Mit \(\mathcal P\) als Menge aller Primzahlen zwischen \(1000\) und \(5000\) lautet die gesuchte Größe
$$\boxed{S=\sum_{p \lt q \lt r,\; p,q,r\in\mathcal P} x_{pqr}.}$$
Die mitgelieferten Programme berechnen daraus
$$S=162619462356610313.$$
Alle drei Lösungen haben dieselbe Struktur: Sie erzeugen die Primliste per Sieb, berechnen für jedes Prim \(p\) den Lucas-Rest \(B \bmod p\), speichern alle paarweisen Inversen \(p_i^{-1}\bmod p_j\) und laufen danach über alle Primtripel. C++ akkumuliert in unsigned __int128, Java in BigInteger, und Python verwendet automatisch beliebig große Ganzzahlen.
Mit \(m=501\) Primzahlen ist das Sieb bis \(5000\) vernachlässigbar klein. Die Lucas-Vorberechnung kostet
$$\sum_{p\in\mathcal P} O(p+\log_p N),$$
weil für jedes Prim eine Faktorientabelle der Länge \(p\) aufgebaut wird, während die Ziffernzerlegung nur wenige Schritte braucht. Die Inversentabelle benötigt \(O(m^2)\) Zeit und Speicher. Dominierend ist die Dreifachschleife mit
$$\binom{m}{3}=\binom{501}{3}=20833250$$
konstant teuren CRT-Zusammenführungen. Insgesamt ist die Laufzeit also \(O(m^3)\), der Speicherbedarf \(O(m^2)\).
Şöyle tanımlayalım:
$$N=10^{18},\qquad K=10^9,\qquad B=\binom{N}{K}.$$
\(1000 \lt p \lt 5000\) aralığındaki her asal \(p\) için önce \(B \bmod p\) hesaplanır. Ardından bu aralıktaki her \(p \lt q \lt r\) asal üçlüsü için
$$0 \le x_{pqr} \lt pqr,\qquad x_{pqr}\equiv B \pmod p,\qquad x_{pqr}\equiv B \pmod q,\qquad x_{pqr}\equiv B \pmod r$$
koşullarını sağlayan tek değer yeniden kurulur ve tüm bu değerler toplanır. Kodun kullandığı eleme işlemi bu aralıkta tam \(501\) asal bulur; dolayısıyla toplam
$$\binom{501}{3}=20833250$$
üçlü incelenir.
Sabit bir asal \(p\) için \(N\) ve \(K\)'yi taban \(p\)'de yazalım:
$$N=\sum_{t=0}^{s} N_t p^t,\qquad K=\sum_{t=0}^{s} K_t p^t,\qquad 0 \le N_t,K_t \lt p.$$
Lucas teoremi şunu söyler:
$$\binom{N}{K}\equiv \prod_{t=0}^{s}\binom{N_t}{K_t}\pmod p.$$
Herhangi bir basamakta \(K_t \gt N_t\) olursa ilgili küçük binom katsayısı sıfırdır; bu yüzden tüm çarpım da \(0 \pmod p\) olur. Yerel çözümlerdeki erken çıkış mantığı tam olarak budur.
Buradaki önemli pratik nokta şudur: \(p \gt 1000\) olduğundan \(10^{18}\) sayısının taban-\(p\) gösteriminde en fazla altı basamak, \(10^9\) sayısında ise en fazla üç basamak vardır. Yani her Lucas hesabı yalnızca birkaç basamak üzerinde çalışır.
\(0 \le b \le a \lt p\) için
$$\binom{a}{b}=\frac{a!}{b!(a-b)!}\pmod p$$
ifadesini kullanırız. \(p\) asal olduğu ve tüm çarpanlar \(p\)'den küçük olduğu için paydadaki terimlerin mod \(p\) altında tersi vardır. Bu yüzden kod
$$\text{fact}[i]=i!\pmod p,\qquad \text{invFact}[i]=(i!)^{-1}\pmod p$$
tablolarını oluşturur ve
$$\binom{a}{b}\equiv \text{fact}[a]\cdot \text{invFact}[b]\cdot \text{invFact}[a-b]\pmod p$$
formülünü kullanır. Tersler de Fermat'nın küçük teoreminden gelir:
$$x^{p-1}\equiv 1\pmod p\quad\Longrightarrow\quad x^{-1}\equiv x^{p-2}\pmod p.$$
Bu nedenle C++, Python ve Java sürümlerinin hepsinde hızlı modüler üs alma yordamı vardır.
C++ kodu \(\binom{30}{12}\) örneğiyle Lucas hesabını doğrular. \(p=11\) için
$$30=2\cdot 11+8,\qquad 12=1\cdot 11+1$$
olduğundan
$$\binom{30}{12}\equiv \binom{8}{1}\binom{2}{1}=8\cdot 2=16\equiv 5\pmod{11}.$$
Gerçek değer \(\binom{30}{12}=86493225\) olup bu sayı da \(11\) ile bölümünden \(5\) kalanını verir. Kaynaktaki kontrol, aynı karşılaştırmayı başka asallar için de tekrar eder.
Şimdi
$$a\equiv B\pmod p,\qquad b\equiv B\pmod q,\qquad c\equiv B\pmod r$$
olsun. \(p,q,r\) farklı asallar olduğundan aralarında ikişer ikişer aralarında asaldırlar; dolayısıyla mod \(pqr\) altında tek bir çözüm vardır.
Kod Garner tarzı iki aşamalı bir birleştirme yapar. Önce
$$x_{pq}=a+p\,t$$
yazılır. Bunun aynı zamanda \(b \pmod q\) olması için
$$a+p\,t\equiv b\pmod q$$
ve dolayısıyla
$$t\equiv (b-a)\,p^{-1}\pmod q$$
gerekir. Böylece
$$x_{pq}=a+p\left((b-a)\,p^{-1}\bmod q\right),\qquad 0 \le x_{pq} \lt pq.$$
Sonra üçüncü modülü eklemek için
$$x=x_{pq}+pq\,u$$
yazarız. \(x\equiv c\pmod r\) koşulu
$$u\equiv (c-x_{pq})(pq)^{-1}\pmod r$$
sonucunu verir. Yani
$$x=x_{pq}+pq\left((c-x_{pq})(pq)^{-1}\bmod r\right),\qquad 0 \le x \lt pqr.$$
Bu son değer tam olarak \(x_{pqr}\)'dir.
Çözümler tüm çiftler için
$$\text{inv}[i][j]\equiv p_i^{-1}\pmod{p_j}$$
değerini saklar. Böylece çarpımın tersi doğrudan bulunur:
$$ (pq)^{-1}\equiv p^{-1}q^{-1}\pmod r.$$
Sonuç olarak kübik döngünün içinde yalnızca birkaç modüler toplama ve çarpma yapılır; her üçlüde yeniden genişletilmiş Öklid çalıştırmak gerekmez.
Yerel C++ kodu ayrıca \(\binom{20}{8}=125970\) ve \(\{7,11,13,17\}\) asal kümesiyle CRT kısmını sınar. \((7,11,13)\) üçlüsü için kalıntılar
$$125970\equiv 5\pmod 7,\qquad 125970\equiv 9\pmod{11},\qquad 125970\equiv 0\pmod{13}$$
olur. Önce \(7\) ve \(11\) birleştirilir:
$$t\equiv (9-5)\cdot 7^{-1}\equiv 4\cdot 8\equiv 10\pmod{11},$$
buradan
$$x_{7,11}=5+7\cdot 10=75$$
elde edilir. Sonra \(13\) eklenir. \(75\equiv 10\pmod{13}\) ve \(77\equiv 12\pmod{13}\) olduğundan
$$u\equiv (0-10)\cdot 12^{-1}\equiv 3\cdot 12\equiv 10\pmod{13},$$
ve sonuç
$$x_{7,11,13}=75+77\cdot 10=845$$
olur. Bu, \(125970 \bmod (7\cdot 11\cdot 13)=845\) ile tam uyumludur. Dört küçük üçlünün toplamı olan \(3803\) değeri de kodda doğrulanan kontrol sonucudur.
\(\mathcal P\), \(1000\) ile \(5000\) arasındaki asallar kümesi olmak üzere aranan nicelik
$$\boxed{S=\sum_{p \lt q \lt r,\; p,q,r\in\mathcal P} x_{pqr}.}$$
Yerel çözümler çalıştırıldığında
$$S=162619462356610313$$
elde edilir.
Üç dildeki çözüm aynı akışı izler: önce asallar eleme ile bulunur, sonra her asal için Lucas kalıntısı hesaplanır, ardından tüm ikili tersler \(p_i^{-1}\bmod p_j\) önceden tabloya yazılır ve son olarak tüm asal üçlüleri üzerinde dönülür. C++ toplamı unsigned __int128 ile, Java BigInteger ile, Python ise doğal olarak keyfi uzunlukta tamsayılarla taşır.
\(m=501\) olsun. \(5000\)'e kadar eleme işlemi çok küçüktür. Lucas önişlemi
$$\sum_{p\in\mathcal P} O(p+\log_p N)$$
maliyetindedir; çünkü her asal için uzunluğu \(p\) olan faktöriyel tabloları kurulur ve basamak döngüsü yalnızca birkaç kez çalışır. İkili tersler tablosu \(O(m^2)\) zaman ve bellek ister. Asıl baskın kısım
$$\binom{m}{3}=\binom{501}{3}=20833250$$
adet sabit-zamanlı CRT birleştirmesi yapan üçlü döngüdür. Bu nedenle toplam süre \(O(m^3)\), bellek kullanımı ise \(O(m^2)\)'dir.
Definimos
$$N=10^{18},\qquad K=10^9,\qquad B=\binom{N}{K}.$$
Para cada primo \(p\) con \(1000 \lt p \lt 5000\) hay que calcular \(B \bmod p\). Después, para cada triple de primos \(p \lt q \lt r\) del mismo intervalo, se reconstruye el único entero \(x_{pqr}\) tal que
$$0 \le x_{pqr} \lt pqr,\qquad x_{pqr}\equiv B \pmod p,\qquad x_{pqr}\equiv B \pmod q,\qquad x_{pqr}\equiv B \pmod r.$$
La respuesta final es la suma de todos esos valores. El cribado usado por el programa encuentra exactamente \(501\) primos entre \(1000\) y \(5000\), así que el sumatorio principal recorre
$$\binom{501}{3}=20833250$$
triples.
Fijado un primo \(p\), escribimos \(N\) y \(K\) en base \(p\):
$$N=\sum_{t=0}^{s} N_t p^t,\qquad K=\sum_{t=0}^{s} K_t p^t,\qquad 0 \le N_t,K_t \lt p.$$
Entonces Lucas afirma que
$$\binom{N}{K}\equiv \prod_{t=0}^{s}\binom{N_t}{K_t}\pmod p.$$
Si en alguna cifra ocurre \(K_t \gt N_t\), entonces el factor local es cero y todo el producto se anula módulo \(p\). Por eso las implementaciones devuelven inmediatamente \(0\) en ese caso.
Como todos los primos del problema satisfacen \(p \gt 1000\), el número \(10^{18}\) tiene a lo sumo seis cifras en base \(p\), y \(10^9\) tiene a lo sumo tres. Así, cada residuo \(B \bmod p\) se obtiene a partir de muy pocos coeficientes binomiales pequeños.
Para \(0 \le b \le a \lt p\), usamos
$$\binom{a}{b}=\frac{a!}{b!(a-b)!}\pmod p.$$
Como \(p\) es primo y todos esos factoriales son menores que \(p\), el denominador es invertible módulo \(p\). El código precalcula
$$\text{fact}[i]=i!\pmod p,\qquad \text{invFact}[i]=(i!)^{-1}\pmod p,$$
y luego evalúa
$$\binom{a}{b}\equiv \text{fact}[a]\cdot \text{invFact}[b]\cdot \text{invFact}[a-b]\pmod p.$$
Las inversas se obtienen con el pequeño teorema de Fermat:
$$x^{p-1}\equiv 1\pmod p\quad\Longrightarrow\quad x^{-1}\equiv x^{p-2}\pmod p.$$
De ahí sale la rutina de exponenciación modular rápida que aparece en C++, Python y Java.
El archivo C++ comprueba Lucas con \(\binom{30}{12}\). Para \(p=11\), las expansiones en base \(11\) son
$$30=2\cdot 11+8,\qquad 12=1\cdot 11+1.$$
Por Lucas,
$$\binom{30}{12}\equiv \binom{8}{1}\binom{2}{1}=16\equiv 5\pmod{11}.$$
El valor exacto \(\binom{30}{12}=86493225\) también deja resto \(5\) módulo \(11\). El programa repite esta verificación para varios primos pequeños.
Supongamos que ya conocemos
$$a\equiv B\pmod p,\qquad b\equiv B\pmod q,\qquad c\equiv B\pmod r.$$
Como \(p,q,r\) son primos distintos, son coprimos dos a dos, y el teorema chino del resto garantiza una única solución módulo \(pqr\).
La implementación usa una reconstrucción en dos pasos, al estilo de Garner. Primero se escribe
$$x_{pq}=a+p\,t.$$
Imponer \(x_{pq}\equiv b\pmod q\) da
$$t\equiv (b-a)\,p^{-1}\pmod q,$$
de modo que
$$x_{pq}=a+p\left((b-a)\,p^{-1}\bmod q\right),\qquad 0 \le x_{pq} \lt pq.$$
Luego se eleva al tercer módulo escribiendo
$$x=x_{pq}+pq\,u.$$
La condición \(x\equiv c\pmod r\) produce
$$u\equiv (c-x_{pq})(pq)^{-1}\pmod r,$$
y por tanto
$$x=x_{pq}+pq\left((c-x_{pq})(pq)^{-1}\bmod r\right),\qquad 0 \le x \lt pqr.$$
Ese valor es precisamente \(x_{pqr}\).
El programa almacena
$$\text{inv}[i][j]\equiv p_i^{-1}\pmod{p_j}.$$
Con eso, la inversa de un producto sale de inmediato:
$$ (pq)^{-1}\equiv p^{-1}q^{-1}\pmod r.$$
Así, cada triple se resuelve con unas pocas sumas y multiplicaciones modulares. No hace falta ejecutar Euclides extendido dentro del bucle cúbico.
El código C++ también verifica la fase CRT con \(\binom{20}{8}=125970\) y el conjunto \(\{7,11,13,17\}\). Para el triple \((7,11,13)\), los residuos exactos son
$$125970\equiv 5\pmod 7,\qquad 125970\equiv 9\pmod{11},\qquad 125970\equiv 0\pmod{13}.$$
Primero se combinan \(7\) y \(11\):
$$t\equiv (9-5)\cdot 7^{-1}\equiv 4\cdot 8\equiv 10\pmod{11},$$
de donde
$$x_{7,11}=5+7\cdot 10=75.$$
Después se incorpora \(13\). Como \(75\equiv 10\pmod{13}\) y \(77\equiv 12\pmod{13}\), obtenemos
$$u\equiv (0-10)\cdot 12^{-1}\equiv 3\cdot 12\equiv 10\pmod{13},$$
y por tanto
$$x_{7,11,13}=75+77\cdot 10=845.$$
Eso coincide con la reducción directa \(125970 \bmod (7\cdot 11\cdot 13)=845\). Sumando los cuatro triples pequeños se obtiene el valor de control \(3803\) usado en el repositorio.
Si \(\mathcal P\) es el conjunto de primos entre \(1000\) y \(5000\), el objetivo es
$$\boxed{S=\sum_{p \lt q \lt r,\; p,q,r\in\mathcal P} x_{pqr}.}$$
Ejecutando las implementaciones locales se obtiene
$$S=162619462356610313.$$
Las tres versiones siguen la misma estructura: generan la lista de primos con un cribado, calculan todos los residuos Lucas \(B \bmod p\), precalculan la tabla completa de inversas \(p_i^{-1}\bmod p_j\) y recorren todos los triples de primos. C++ acumula en unsigned __int128, Java usa BigInteger y Python aprovecha sus enteros de precisión arbitraria.
Sea \(m=501\) el número de primos. El cribado hasta \(5000\) es despreciable. El preprocesamiento de Lucas cuesta
$$\sum_{p\in\mathcal P} O(p+\log_p N),$$
porque para cada primo se construyen tablas de factoriales de tamaño \(p\), mientras que el recorrido por cifras es muy corto. La tabla de inversas requiere \(O(m^2)\) tiempo y memoria. La fase dominante es el bucle sobre triples, con
$$\binom{m}{3}=\binom{501}{3}=20833250$$
fusiones CRT de costo constante. En consecuencia, el tiempo total es \(O(m^3)\) y la memoria es \(O(m^2)\).
记
$$N=10^{18},\qquad K=10^9,\qquad B=\binom{N}{K}.$$
对每个满足 \(1000 \lt p \lt 5000\) 的素数 \(p\)，先求 \(B \bmod p\)。然后对区间内每个素数三元组 \(p \lt q \lt r\)，重建唯一的整数 \(x_{pqr}\)，使得
$$0 \le x_{pqr} \lt pqr,\qquad x_{pqr}\equiv B \pmod p,\qquad x_{pqr}\equiv B \pmod q,\qquad x_{pqr}\equiv B \pmod r.$$
题目要求把所有这样的 \(x_{pqr}\) 相加。代码中的筛法表明该区间内一共有 \(501\) 个素数，因此主循环一共处理
$$\binom{501}{3}=20833250$$
个三元组。
固定一个素数 \(p\)，把 \(N\) 和 \(K\) 写成 \(p\) 进制：
$$N=\sum_{t=0}^{s} N_t p^t,\qquad K=\sum_{t=0}^{s} K_t p^t,\qquad 0 \le N_t,K_t \lt p.$$
Lucas 定理给出
$$\binom{N}{K}\equiv \prod_{t=0}^{s}\binom{N_t}{K_t}\pmod p.$$
如果某一位出现 \(K_t \gt N_t\)，那么对应的小组合数就是 \(0\)，整个乘积也就变成 \(0 \pmod p\)。这正是实现里“若某位不合法就立刻返回零”的原因。
由于题目中的素数都满足 \(p \gt 1000\)，所以 \(10^{18}\) 的 \(p\) 进制长度最多只有六位，而 \(10^9\) 最多三位。也就是说，每个 \(B \bmod p\) 都只需要处理很少几个数字位。
当 \(0 \le b \le a \lt p\) 时，我们使用
$$\binom{a}{b}=\frac{a!}{b!(a-b)!}\pmod p.$$
因为 \(p\) 是素数，且这些数都小于 \(p\)，分母在模 \(p\) 下可逆。于是程序预先计算
$$\text{fact}[i]=i!\pmod p,\qquad \text{invFact}[i]=(i!)^{-1}\pmod p,$$
再用
$$\binom{a}{b}\equiv \text{fact}[a]\cdot \text{invFact}[b]\cdot \text{invFact}[a-b]\pmod p$$
得到小组合数。逆元来自费马小定理：
$$x^{p-1}\equiv 1\pmod p\quad\Longrightarrow\quad x^{-1}\equiv x^{p-2}\pmod p.$$
因此三种语言的实现里都出现了快速幂函数。
C++ 源码用 \(\binom{30}{12}\) 做 Lucas 检查。取 \(p=11\)，则
$$30=2\cdot 11+8,\qquad 12=1\cdot 11+1.$$
由 Lucas 定理可得
$$\binom{30}{12}\equiv \binom{8}{1}\binom{2}{1}=16\equiv 5\pmod{11}.$$
而精确值 \(\binom{30}{12}=86493225\) 对 \(11\) 取模同样等于 \(5\)。源码还对若干其他素数重复了这个比对。
设我们已经知道
$$a\equiv B\pmod p,\qquad b\equiv B\pmod q,\qquad c\equiv B\pmod r.$$
由于 \(p,q,r\) 是互不相同的素数，它们两两互素，所以模 \(pqr\) 的解唯一存在。
程序采用类似 Garner 的两步重建。先写成
$$x_{pq}=a+p\,t.$$
要求它同时满足 \(x_{pq}\equiv b\pmod q\)，就有
$$t\equiv (b-a)\,p^{-1}\pmod q,$$
因此
$$x_{pq}=a+p\left((b-a)\,p^{-1}\bmod q\right),\qquad 0 \le x_{pq} \lt pq.$$
再提升到第三个模数，写成
$$x=x_{pq}+pq\,u.$$
由 \(x\equiv c\pmod r\) 得
$$u\equiv (c-x_{pq})(pq)^{-1}\pmod r,$$
于是
$$x=x_{pq}+pq\left((c-x_{pq})(pq)^{-1}\bmod r\right),\qquad 0 \le x \lt pqr.$$
这个 \(x\) 就是题目中定义的 \(x_{pqr}\)。
实现中保存了所有
$$\text{inv}[i][j]\equiv p_i^{-1}\pmod{p_j}.$$
于是乘积的逆元可以直接写成
$$ (pq)^{-1}\equiv p^{-1}q^{-1}\pmod r.$$
这样在三重循环内部，只需要少量模加法和模乘法，不必对每个三元组重新做一次扩展欧几里得算法。
本地 C++ 代码还用 \(\binom{20}{8}=125970\) 和素数集合 \(\{7,11,13,17\}\) 检验 CRT 阶段。对于三元组 \((7,11,13)\)，有
$$125970\equiv 5\pmod 7,\qquad 125970\equiv 9\pmod{11},\qquad 125970\equiv 0\pmod{13}.$$
先合并 \(7\) 与 \(11\)：
$$t\equiv (9-5)\cdot 7^{-1}\equiv 4\cdot 8\equiv 10\pmod{11},$$
所以
$$x_{7,11}=5+7\cdot 10=75.$$
再合并 \(13\)。因为 \(75\equiv 10\pmod{13}\)，而 \(77\equiv 12\pmod{13}\)，于是
$$u\equiv (0-10)\cdot 12^{-1}\equiv 3\cdot 12\equiv 10\pmod{13},$$
从而
$$x_{7,11,13}=75+77\cdot 10=845.$$
这与直接计算 \(125970 \bmod (7\cdot 11\cdot 13)=845\) 完全一致。四个小三元组的总和为 \(3803\)，这也是源码里的检查值。
设 \(\mathcal P\) 为 \(1000\) 到 \(5000\) 之间的素数集合，则目标为
$$\boxed{S=\sum_{p \lt q \lt r,\; p,q,r\in\mathcal P} x_{pqr}.}$$
运行仓库中的实现可得到
$$S=162619462356610313.$$
三份代码结构完全一致：先筛出素数列表，再对每个素数计算 Lucas 残差 \(B \bmod p\)，然后预处理完整的两两逆元表 \(p_i^{-1}\bmod p_j\)，最后枚举所有素数三元组。C++ 用 unsigned __int128 累加总和，Java 用 BigInteger，Python 直接使用任意精度整数。
记素数个数为 \(m=501\)。筛到 \(5000\) 的成本很小。Lucas 预处理的代价为
$$\sum_{p\in\mathcal P} O(p+\log_p N),$$
因为每个素数都要建立长度为 \(p\) 的阶乘表，而 \(p\) 进制拆位只需很少几步。逆元矩阵需要 \(O(m^2)\) 的时间与空间。主导项是三重循环，它执行
$$\binom{m}{3}=\binom{501}{3}=20833250$$
次常数时间的 CRT 合并。因此总时间复杂度是 \(O(m^3)\)，空间复杂度是 \(O(m^2)\)。
Положим
$$N=10^{18},\qquad K=10^9,\qquad B=\binom{N}{K}.$$
Для каждого простого \(p\) из диапазона \(1000 \lt p \lt 5000\) нужно вычислить \(B \bmod p\). Затем для каждого тройного набора простых \(p \lt q \lt r\) восстанавливается единственное число \(x_{pqr}\), удовлетворяющее условиям
$$0 \le x_{pqr} \lt pqr,\qquad x_{pqr}\equiv B \pmod p,\qquad x_{pqr}\equiv B \pmod q,\qquad x_{pqr}\equiv B \pmod r.$$
Итоговый ответ равен сумме всех таких значений. Решето из программы находит в этом интервале ровно \(501\) простое число, поэтому внешний перебор охватывает
$$\binom{501}{3}=20833250$$
троек.
Для фиксированного простого \(p\) запишем \(N\) и \(K\) в системе счисления по основанию \(p\):
$$N=\sum_{t=0}^{s} N_t p^t,\qquad K=\sum_{t=0}^{s} K_t p^t,\qquad 0 \le N_t,K_t \lt p.$$
Тогда по теореме Лукаса
$$\binom{N}{K}\equiv \prod_{t=0}^{s}\binom{N_t}{K_t}\pmod p.$$
Если на каком-то разряде выполняется \(K_t \gt N_t\), то соответствующий локальный биномиальный коэффициент равен нулю, а значит и всё произведение равно \(0 \pmod p\). Именно поэтому локальная функция в исходниках сразу возвращает ноль при таком сравнении цифр.
Поскольку все рассматриваемые простые больше \(1000\), число \(10^{18}\) имеет в базе \(p\) не более шести цифр, а число \(10^9\) не более трёх. Значит, для каждого простого остаток вычисляется по очень короткому циклу.
Для \(0 \le b \le a \lt p\) используется формула
$$\binom{a}{b}=\frac{a!}{b!(a-b)!}\pmod p.$$
Так как \(p\) простое и все множители меньше \(p\), знаменатель обратим по модулю \(p\). Поэтому программа заранее строит таблицы
$$\text{fact}[i]=i!\pmod p,\qquad \text{invFact}[i]=(i!)^{-1}\pmod p,$$
после чего использует
$$\binom{a}{b}\equiv \text{fact}[a]\cdot \text{invFact}[b]\cdot \text{invFact}[a-b]\pmod p.$$
Обратные элементы получаются из малой теоремы Ферма:
$$x^{p-1}\equiv 1\pmod p\quad\Longrightarrow\quad x^{-1}\equiv x^{p-2}\pmod p.$$
Именно этим объясняется наличие быстрой модульной степени во всех трёх реализациях.
В C++-файле теорема Лукаса проверяется на \(\binom{30}{12}\). При \(p=11\)
$$30=2\cdot 11+8,\qquad 12=1\cdot 11+1.$$
Следовательно,
$$\binom{30}{12}\equiv \binom{8}{1}\binom{2}{1}=16\equiv 5\pmod{11}.$$
Точное значение \(\binom{30}{12}=86493225\) действительно даёт остаток \(5\) по модулю \(11\). В локальных проверках этот же тест повторяется для нескольких простых чисел.
Пусть уже известны остатки
$$a\equiv B\pmod p,\qquad b\equiv B\pmod q,\qquad c\equiv B\pmod r.$$
Поскольку \(p,q,r\) различны и попарно взаимно просты, существует единственное решение по модулю \(pqr\).
Код использует двухшаговую реконструкцию в стиле Гарнера. Сначала пишем
$$x_{pq}=a+p\,t.$$
Из условия \(x_{pq}\equiv b\pmod q\) следует
$$t\equiv (b-a)\,p^{-1}\pmod q,$$
и потому
$$x_{pq}=a+p\left((b-a)\,p^{-1}\bmod q\right),\qquad 0 \le x_{pq} \lt pq.$$
Далее поднимаем решение на третий модуль:
$$x=x_{pq}+pq\,u.$$
Условие \(x\equiv c\pmod r\) даёт
$$u\equiv (c-x_{pq})(pq)^{-1}\pmod r,$$
так что
$$x=x_{pq}+pq\left((c-x_{pq})(pq)^{-1}\bmod r\right),\qquad 0 \le x \lt pqr.$$
Это и есть искомое значение \(x_{pqr}\).
Реализации заранее вычисляют
$$\text{inv}[i][j]\equiv p_i^{-1}\pmod{p_j}.$$
После этого обратный элемент к произведению получается сразу:
$$ (pq)^{-1}\equiv p^{-1}q^{-1}\pmod r.$$
Значит, внутри кубического перебора остаются только несколько модульных сложений и умножений. Никакой новый расширенный алгоритм Евклида для каждой тройки не нужен.
Локальный C++-код также проверяет этап CRT на \(\binom{20}{8}=125970\) и множестве простых \(\{7,11,13,17\}\). Для тройки \((7,11,13)\) точные остатки таковы:
$$125970\equiv 5\pmod 7,\qquad 125970\equiv 9\pmod{11},\qquad 125970\equiv 0\pmod{13}.$$
Сначала объединяем \(7\) и \(11\):
$$t\equiv (9-5)\cdot 7^{-1}\equiv 4\cdot 8\equiv 10\pmod{11},$$
поэтому
$$x_{7,11}=5+7\cdot 10=75.$$
Затем добавляем модуль \(13\). Так как \(75\equiv 10\pmod{13}\), а \(77\equiv 12\pmod{13}\), имеем
$$u\equiv (0-10)\cdot 12^{-1}\equiv 3\cdot 12\equiv 10\pmod{13},$$
и значит
$$x_{7,11,13}=75+77\cdot 10=845.$$
Это в точности совпадает с прямым вычислением \(125970 \bmod (7\cdot 11\cdot 13)=845\). Сумма по четырём маленьким тройкам равна \(3803\), и именно этот контрольный результат проверяется в исходниках.
Если \(\mathcal P\) обозначает множество простых чисел между \(1000\) и \(5000\), то требуется вычислить
$$\boxed{S=\sum_{p \lt q \lt r,\; p,q,r\in\mathcal P} x_{pqr}.}$$
Запуск локальных решений даёт
$$S=162619462356610313.$$
Во всех трёх языках структура одинакова: сначала строится список простых через решето, затем для каждого простого считается остаток по Лукасу, потом заполняется полная таблица попарных обратных \(p_i^{-1}\bmod p_j\), и лишь после этого идёт перебор всех троек. Версия на C++ суммирует в unsigned __int128, версия на Java использует BigInteger, а Python полагается на встроенные целые произвольной длины.
Пусть \(m=501\). Решето до \(5000\) мало по стоимости. Предобработка Лукаса требует
$$\sum_{p\in\mathcal P} O(p+\log_p N),$$
потому что для каждого простого строятся таблицы факториалов длины \(p\), а разложение по цифрам очень короткое. Таблица обратных требует \(O(m^2)\) времени и памяти. Основной вклад даёт тройной цикл с
$$\binom{m}{3}=\binom{501}{3}=20833250$$
слияниями CRT константной стоимости. Следовательно, итоговая асимптотика равна \(O(m^3)\) по времени и \(O(m^2)\) по памяти.
لنعرّف
$$N=10^{18},\qquad K=10^9,\qquad B=\binom{N}{K}.$$
لكل عدد أولي \(p\) يحقق \(1000 \lt p \lt 5000\)، نحتاج أولًا إلى حساب \(B \bmod p\). بعد ذلك، ولكل ثلاثي أولي \(p \lt q \lt r\) في هذا المجال، نعيد بناء العدد الوحيد \(x_{pqr}\) الذي يحقق
$$0 \le x_{pqr} \lt pqr,\qquad x_{pqr}\equiv B \pmod p,\qquad x_{pqr}\equiv B \pmod q,\qquad x_{pqr}\equiv B \pmod r.$$
الجواب النهائي هو مجموع جميع هذه القيم. الغربال المستخدم في الشيفرة يجد بالضبط \(501\) عددًا أوليًا في المجال، ولذلك فإن الحلقة الرئيسية تمر على
$$\binom{501}{3}=20833250$$
ثلاثية أولية.
لعدد أولي ثابت \(p\)، نكتب \(N\) و\(K\) بالأساس \(p\):
$$N=\sum_{t=0}^{s} N_t p^t,\qquad K=\sum_{t=0}^{s} K_t p^t,\qquad 0 \le N_t,K_t \lt p.$$
عندئذ تعطي مبرهنة لوكاس
$$\binom{N}{K}\equiv \prod_{t=0}^{s}\binom{N_t}{K_t}\pmod p.$$
إذا تحقق في أي خانة أن \(K_t \gt N_t\)، فإن معامل تلك الخانة يساوي صفرًا، وبالتالي يصبح حاصل الضرب كله \(0 \pmod p\). ولهذا تعيد الدالة في التنفيذ القيمة صفرًا مباشرة عند ظهور هذه الحالة.
وبما أن كل الأعداد الأولية في المسألة أكبر من \(1000\)، فإن \(10^{18}\) يملك على الأكثر ست خانات في الأساس \(p\)، بينما يملك \(10^9\) على الأكثر ثلاث خانات. لذلك فإن حساب كل باقي \(B \bmod p\) يتم عبر عدد صغير جدًا من الخطوات.
عندما يكون \(0 \le b \le a \lt p\)، نستخدم الصيغة
$$\binom{a}{b}=\frac{a!}{b!(a-b)!}\pmod p.$$
وبما أن \(p\) أولي وجميع هذه القيم أصغر من \(p\)، فإن المقام قابل للعكس بترديد \(p\). لذلك تُنشئ الشيفرة الجدولين
$$\text{fact}[i]=i!\pmod p,\qquad \text{invFact}[i]=(i!)^{-1}\pmod p,$$
ثم تحسب
$$\binom{a}{b}\equiv \text{fact}[a]\cdot \text{invFact}[b]\cdot \text{invFact}[a-b]\pmod p.$$
أما المعكوسات فتعتمد على مبرهنة فيرما الصغرى:
$$x^{p-1}\equiv 1\pmod p\quad\Longrightarrow\quad x^{-1}\equiv x^{p-2}\pmod p.$$
ولهذا تحتوي جميع النسخ على دالة أسّية معيارية سريعة.
ملف C++ يتحقق من مبرهنة لوكاس باستخدام \(\binom{30}{12}\). عند \(p=11\) نحصل على
$$30=2\cdot 11+8,\qquad 12=1\cdot 11+1.$$
ومن ثم
$$\binom{30}{12}\equiv \binom{8}{1}\binom{2}{1}=16\equiv 5\pmod{11}.$$
والقيمة الدقيقة \(\binom{30}{12}=86493225\) تعطي أيضًا باقيًا يساوي \(5\) modulo \(11\). الشيفرة المحلية تكرر هذا الفحص لعدة أوليات أخرى.
لنفترض أننا عرفنا بالفعل
$$a\equiv B\pmod p,\qquad b\equiv B\pmod q,\qquad c\equiv B\pmod r.$$
بما أن \(p,q,r\) أوليات مختلفة فهي متباينة زوجيًا، وبالتالي تضمن مبرهنة البواقي الصينية وجود حل وحيد modulo \(pqr\).
التنفيذ يستخدم إعادة بناء من خطوتين على نمط Garner. أولًا نكتب
$$x_{pq}=a+p\,t.$$
ولكي يكون هذا العدد موافقًا لـ \(b \pmod q\) يجب أن يكون
$$t\equiv (b-a)\,p^{-1}\pmod q,$$
ومن هنا
$$x_{pq}=a+p\left((b-a)\,p^{-1}\bmod q\right),\qquad 0 \le x_{pq} \lt pq.$$
بعد ذلك نرفع الحل إلى الموديل الثالث بكتابة
$$x=x_{pq}+pq\,u.$$
ومن الشرط \(x\equiv c\pmod r\) نستنتج
$$u\equiv (c-x_{pq})(pq)^{-1}\pmod r,$$
أي أن
$$x=x_{pq}+pq\left((c-x_{pq})(pq)^{-1}\bmod r\right),\qquad 0 \le x \lt pqr.$$
وهذا هو بالضبط العدد \(x_{pqr}\) المطلوب.
التنفيذات تخزّن جميع القيم
$$\text{inv}[i][j]\equiv p_i^{-1}\pmod{p_j}.$$
وعندها يصبح معكوس حاصل الضرب متاحًا مباشرة:
$$ (pq)^{-1}\equiv p^{-1}q^{-1}\pmod r.$$
وهكذا تحتاج كل ثلاثية داخل الحلقة التكعيبية إلى عدد قليل من عمليات الجمع والضرب المعياريين فقط، من دون إعادة تشغيل خوارزمية إقليدس الموسعة في كل مرة.
شيفرة C++ المحلية تختبر أيضًا مرحلة CRT على \(\binom{20}{8}=125970\) ومع المجموعة \(\{7,11,13,17\}\). بالنسبة للثلاثية \((7,11,13)\) تكون البواقي الدقيقة
$$125970\equiv 5\pmod 7,\qquad 125970\equiv 9\pmod{11},\qquad 125970\equiv 0\pmod{13}.$$
أولًا ندمج \(7\) و\(11\):
$$t\equiv (9-5)\cdot 7^{-1}\equiv 4\cdot 8\equiv 10\pmod{11},$$
فنحصل على
$$x_{7,11}=5+7\cdot 10=75.$$
ثم نضيف \(13\). بما أن \(75\equiv 10\pmod{13}\) و\(77\equiv 12\pmod{13}\)، فإن
$$u\equiv (0-10)\cdot 12^{-1}\equiv 3\cdot 12\equiv 10\pmod{13},$$
ومن ثم
$$x_{7,11,13}=75+77\cdot 10=845.$$
وهذا يطابق مباشرة \(125970 \bmod (7\cdot 11\cdot 13)=845\). كما أن مجموع الأربع ثلاثيات الصغيرة يساوي \(3803\)، وهو رقم التحقق الموجود في المصدر.
إذا رمزنا إلى مجموعة الأوليات بين \(1000\) و\(5000\) بالرمز \(\mathcal P\)، فإن الكمية المطلوبة هي
$$\boxed{S=\sum_{p \lt q \lt r,\; p,q,r\in\mathcal P} x_{pqr}.}$$
وعند تشغيل الحلول المحلية نحصل على
$$S=162619462356610313.$$
النسخ الثلاث تتبع البنية نفسها: توليد قائمة الأوليات بغربال، ثم حساب بواقي Lucas لكل أولي، ثم بناء جدول كامل للمعكوسات \(p_i^{-1}\bmod p_j\)، وأخيرًا المرور على جميع الثلاثيات الأولية. نسخة C++ تجمع الناتج في unsigned __int128، ونسخة Java تستخدم BigInteger، أما Python فيعتمد على أعداد صحيحة ذات دقة غير محدودة.
لنرمز بعدد الأوليات في المجال بـ \(m=501\). الغربال حتى \(5000\) صغير الكلفة. أمّا تهيئة Lucas فتكلف
$$\sum_{p\in\mathcal P} O(p+\log_p N),$$
لأن كل أولي يبني جدول فواعل بطول \(p\)، بينما تفكيك الأرقام في الأساس \(p\) قصير جدًا. جدول المعكوسات يحتاج \(O(m^2)\) زمنًا وذاكرة. الجزء المهيمن هو الحلقة على جميع الثلاثيات، والتي تنفذ
$$\binom{m}{3}=\binom{501}{3}=20833250$$
عملية دمج CRT بزمن ثابت. لذلك يكون الزمن الكلي \(O(m^3)\)، واستهلاك الذاكرة \(O(m^2)\).