Let \(k\) be the number of decimal digits of \(d\). We want to count
$$F(n,d)=\#\left\{x\in \mathbb{Z}_{>0}:x\mid n!,\ x\equiv d \pmod{10^k}\right\},$$
with the final result reported modulo \(10^{16}+61\). For Project Euler 474 the concrete input is \(n=10^6\) and \(d=65432\). The difficulty is that \(n!\) has an enormous number of divisors, so the computation must work with prime exponents and modular residues rather than explicit divisor generation.
Every divisor of \(n!\) is determined by choosing one exponent for each prime \(p\le n\). The implementation uses that representation, then separates the forced powers of \(2\) and \(5\) coming from the target suffix, and finally runs a dynamic program on the unit group modulo a reduced power of \(10\).
If
$$n!=\prod_{p\le n} p^{e_p},\qquad e_p=v_p(n!)=\sum_{m\ge 1}\left\lfloor\frac{n}{p^m}\right\rfloor,$$
then every divisor \(x\mid n!\) has the form
$$x=\prod_{p\le n} p^{a_p},\qquad 0\le a_p\le e_p.$$
So the problem is equivalent to counting exponent vectors \((a_p)\) that place the resulting divisor in one prescribed residue class modulo \(10^k\).
Write
$$d=2^{\alpha}5^{\beta}\tau,\qquad \gcd(\tau,10)=1,$$
where \(\alpha=v_2(d)\) and \(\beta=v_5(d)\). In the main regime used for the actual Euler input, both \(\alpha\lt k\) and \(\beta\lt k\). Define
$$f=2^{\alpha}5^{\beta},\qquad M=\frac{10^k}{f}=2^{k-\alpha}5^{k-\beta},\qquad t=\frac{d}{f}.$$
If a divisor \(x\) satisfies \(x\equiv d \pmod{10^k}\), then \(x\) must contain exactly \(\alpha\) factors of \(2\) and exactly \(\beta\) factors of \(5\). Dividing the congruence by \(f\) gives
$$\frac{x}{f}\equiv t \pmod{M},$$
and \(t\) is coprime to \(10\), so \(x/f\) is automatically a unit modulo \(M\). Therefore the exponents of \(2\) and \(5\) are fixed, and only the primes \(p\ne 2,5\) remain free.
If \(\alpha>v_2(n!)\) or \(\beta>v_5(n!)\), the answer is immediately \(0\). The implementations also keep a separate tiny brute-force branch for the checkpoint \((n,d)=(12,12)\), because there \(\alpha=k\) and the unit reduction above does not apply.
After factoring out the fixed part \(f\), the remaining unit part of a divisor is
$$u=\prod_{\substack{p\le n\\ p\ne 2,5}} p^{a_p}\pmod{M}.$$
Only residues coprime to \(M\) can occur, so the state space has size
$$U=\varphi(M).$$
Define \(A(r)\) as the number of choices of exponents for the primes processed so far that produce residue \(r\) modulo \(M\). Initially only the empty product is present, so
$$A(1)=1,\qquad A(r)=0\text{ for }r\ne 1.$$
When processing one prime \(p\ne 2,5\) with exponent range \(0\le j\le e_p\), the transition is
$$A_{\mathrm{new}}(x)=\sum_{j=0}^{e_p} A_{\mathrm{old}}(x p^{-j}),$$
where the inverse is taken modulo \(M\). This is correct because choosing exponent \(j\) for \(p\) multiplies the current residue by \(p^j\).
Since \(\gcd(p,M)=1\), multiplication by \(p\) permutes the unit residues modulo \(M\). Hence the unit set splits into disjoint cycles
$$u_0,\ u_1,\ \dots,\ u_{L-1},\qquad u_{i+1}\equiv p\,u_i \pmod{M}.$$
On one cycle, the transition becomes a cyclic convolution:
$$A_{\mathrm{new}}(u_i)=\sum_{j=0}^{e_p} A_{\mathrm{old}}(u_{i-j}),$$
with indices interpreted modulo \(L\). Write
$$e_p+1=qL+r,\qquad 0\le r\lt L.$$
Then each new value contains \(q\) complete wraps around the cycle plus an extra window of length \(r\):
$$A_{\mathrm{new}}(u_i)=q\sum_{m=0}^{L-1}A_{\mathrm{old}}(u_m)+\sum_{j=0}^{r-1}A_{\mathrm{old}}(u_{i-j}).$$
The first term is constant on the whole cycle, and the second can be updated with a sliding window. That reduces the work for one prime from \(O(e_pL)\) on a cycle to \(O(L)\).
After all primes \(p\ne 2,5\) have been processed, the desired count is simply the state at the target unit residue \(t\):
$$F(n,d)\equiv A(t)\pmod{10^{16}+61}.$$
All arithmetic in the dynamic program is performed modulo \(10^{16}+61\), while the residue classes themselves live modulo \(M\).
Here \(k=5\) and
$$65432=2^3\cdot 8179.$$
So \(\alpha=3\), \(\beta=0\), \(f=8\), and
$$M=\frac{10^5}{8}=12500,\qquad t=\frac{65432}{8}=8179.$$
Because \(\gcd(8179,12500)=1\), the main unit-case algorithm applies directly. The counted divisors are exactly those divisors of \(n!\) whose exponent of \(2\) is \(3\), whose exponent of \(5\) is \(0\), and whose remaining unit part is congruent to \(8179\) modulo \(12500\). The dynamic program therefore runs on
$$U=\varphi(12500)=12500\left(1-\frac12\right)\left(1-\frac15\right)=5000$$
unit residues instead of on the enormous set of all divisors of \(10^6!\).
The C++, Python, and Java implementations follow the same numerical strategy. They first sieve all primes up to \(n\), compute \(v_p(n!)\) by Legendre's formula, and determine whether the target suffix falls into the generic unit-case or into the one small checkpoint handled separately by direct divisor enumeration.
In the unit-case, the implementation builds the list of residues modulo \(M\) that are coprime to \(M\), maps each such residue to an array position, and initializes a one-dimensional dynamic-programming table with count \(1\) at residue \(1\).
For each prime \(p\ne 2,5\), it reuses the cycle decomposition of the permutation \(r\mapsto pr \pmod{M}\). Inside each cycle it computes the full-cycle contribution once, then advances the remaining partial sum with a sliding window, so every state in that cycle is updated in constant amortized time.
Because the answer modulus is slightly above \(10^{16}\), the implementations also use overflow-safe modular arithmetic for the running totals. After all odd primes have been processed, the entry corresponding to \(t=d/f\) is the required count modulo \(10^{16}+61\).
Let \(U=\varphi(M)\), and let \(C\) be the number of distinct residue classes \(p\bmod M\) that occur among the primes \(p\le n\) with \(p\ne 2,5\). The prime sieve costs \(O(n\log\log n)\) time and \(O(n)\) memory. Building the cached cycle layouts costs \(O(CU)\) time and \(O(CU)\) memory in the current implementation. Once those layouts are available, each prime transition touches each unit residue only once, so the dynamic-programming work is \(O(\pi(n)\,U)\).
For the actual Euler input, \(M=12500\) and \(U=5000\), which is why this approach is practical even though \(10^6!\) itself has an astronomically large divisor set.
Sei \(k\) die Anzahl der Dezimalstellen von \(d\). Gesucht ist
$$F(n,d)=\#\left\{x\in \mathbb{Z}_{>0}:x\mid n!,\ x\equiv d \pmod{10^k}\right\},$$
wobei das Endergebnis modulo \(10^{16}+61\) ausgegeben wird. Für Problem 474 sind die konkreten Eingaben \(n=10^6\) und \(d=65432\). Die Zahl der Teiler von \(n!\) ist riesig, daher arbeitet die Lösung nicht mit expliziten Teilern, sondern nur mit Primexponenten und Restklassen.
Jeder Teiler von \(n!\) wird durch die Wahl eines Exponenten für jede Primzahl \(p\le n\) festgelegt. Die Implementierung nutzt genau diese Darstellung, trennt danach die vom Endblock erzwungenen Zweier- und Fünferpotenzen ab und führt schließlich eine dynamische Programmierung auf der Einheitengruppe modulo einer reduzierten Zehnerpotenz aus.
Ist
$$n!=\prod_{p\le n} p^{e_p},\qquad e_p=v_p(n!)=\sum_{m\ge 1}\left\lfloor\frac{n}{p^m}\right\rfloor,$$
dann hat jeder Teiler \(x\mid n!\) die Form
$$x=\prod_{p\le n} p^{a_p},\qquad 0\le a_p\le e_p.$$
Damit wird das Problem zu einer Zählung der Exponentenvektoren \((a_p)\), deren Produkt in genau eine vorgegebene Restklasse modulo \(10^k\) fällt.
Schreibe
$$d=2^{\alpha}5^{\beta}\tau,\qquad \gcd(\tau,10)=1,$$
wobei \(\alpha=v_2(d)\) und \(\beta=v_5(d)\) gilt. Im Hauptfall, der auch den eigentlichen Euler-Eingaben entspricht, sind sowohl \(\alpha\lt k\) als auch \(\beta\lt k\). Setze
$$f=2^{\alpha}5^{\beta},\qquad M=\frac{10^k}{f}=2^{k-\alpha}5^{k-\beta},\qquad t=\frac{d}{f}.$$
Erfüllt ein Teiler \(x\) die Kongruenz \(x\equiv d \pmod{10^k}\), dann enthält \(x\) genau \(\alpha\) Faktoren \(2\) und genau \(\beta\) Faktoren \(5\). Durch Division der Kongruenz durch \(f\) erhält man
$$\frac{x}{f}\equiv t \pmod{M},$$
und weil \(t\) teilerfremd zu \(10\) ist, ist \(x/f\) automatisch eine Einheit modulo \(M\). Damit sind die Exponenten von \(2\) und \(5\) fest vorgegeben; frei bleiben nur die Primzahlen \(p\ne 2,5\).
Falls \(\alpha>v_2(n!)\) oder \(\beta>v_5(n!)\) gilt, ist die Antwort sofort \(0\). Zusätzlich enthält die Implementierung einen separaten Brute-Force-Zweig für den kleinen Kontrollfall \((n,d)=(12,12)\), weil dort \(\alpha=k\) ist und die obige Einheitsreduktion nicht greift.
Nach dem Herausziehen des festen Faktors \(f\) bleibt vom Teiler nur noch der Einheitenteil
$$u=\prod_{\substack{p\le n\\ p\ne 2,5}} p^{a_p}\pmod{M}.$$
Es können nur zu \(M\) teilerfremde Reste auftreten, also hat der Zustandsraum Größe
$$U=\varphi(M).$$
Sei \(A(r)\) die Anzahl der bisherigen Exponentenwahlen, die den Rest \(r\) modulo \(M\) erzeugen. Am Anfang existiert nur das leere Produkt, also
$$A(1)=1,\qquad A(r)=0\text{ für }r\ne 1.$$
Beim Verarbeiten einer Primzahl \(p\ne 2,5\) mit Exponentenbereich \(0\le j\le e_p\) lautet der Übergang
$$A_{\mathrm{neu}}(x)=\sum_{j=0}^{e_p} A_{\mathrm{alt}}(x p^{-j}),$$
wobei das Inverse modulo \(M\) genommen wird. Das ist genau richtig, weil die Wahl des Exponenten \(j\) den bisherigen Rest mit \(p^j\) multipliziert.
Da \(\gcd(p,M)=1\) ist, permutiert die Multiplikation mit \(p\) die Einheiten modulo \(M\). Daher zerfällt die Einheitenmenge in disjunkte Zyklen
$$u_0,\ u_1,\ \dots,\ u_{L-1},\qquad u_{i+1}\equiv p\,u_i \pmod{M}.$$
Auf einem einzelnen Zyklus wird der Übergang zu einer zyklischen Faltung:
$$A_{\mathrm{neu}}(u_i)=\sum_{j=0}^{e_p} A_{\mathrm{alt}}(u_{i-j}),$$
wobei die Indizes modulo \(L\) gelesen werden. Schreibe
$$e_p+1=qL+r,\qquad 0\le r\lt L.$$
Dann enthält jeder neue Wert \(q\) vollständige Umläufe über den Zyklus plus ein zusätzliches Fenster der Länge \(r\):
$$A_{\mathrm{neu}}(u_i)=q\sum_{m=0}^{L-1}A_{\mathrm{alt}}(u_m)+\sum_{j=0}^{r-1}A_{\mathrm{alt}}(u_{i-j}).$$
Der erste Term ist auf dem ganzen Zyklus konstant, und der zweite lässt sich mit einem Sliding Window aktualisieren. Dadurch sinkt die Arbeit pro Primzahl auf einem Zyklus von \(O(e_pL)\) auf \(O(L)\).
Nachdem alle Primzahlen \(p\ne 2,5\) verarbeitet wurden, ist die gesuchte Anzahl einfach der Zustand am Zielrest \(t\):
$$F(n,d)\equiv A(t)\pmod{10^{16}+61}.$$
Die Restklassen leben also modulo \(M\), während die Zählwerte der dynamischen Programmierung stets modulo \(10^{16}+61\) geführt werden.
Hier ist \(k=5\) und
$$65432=2^3\cdot 8179.$$
Damit sind \(\alpha=3\), \(\beta=0\), \(f=8\) und
$$M=\frac{10^5}{8}=12500,\qquad t=\frac{65432}{8}=8179.$$
Weil \(\gcd(8179,12500)=1\) gilt, kann der Hauptalgorithmus direkt angewandt werden. Gezählt werden also genau die Teiler von \(n!\), deren Zweierexponent gleich \(3\), deren Fünferexponent gleich \(0\) und deren verbleibender Einheitenteil kongruent zu \(8179\) modulo \(12500\) ist. Die dynamische Programmierung läuft damit nur über
$$U=\varphi(12500)=12500\left(1-\frac12\right)\left(1-\frac15\right)=5000$$
Einheitsreste statt über die astronomisch große Gesamtmenge aller Teiler von \(10^6!\).
Die C++-, Python- und Java-Implementierungen folgen derselben numerischen Strategie. Zuerst werden alle Primzahlen bis \(n\) gesiebt, dann die Exponenten \(v_p(n!)\) mit der Legendre-Formel berechnet und anschließend wird entschieden, ob der Ziel-Endblock im allgemeinen Einheitsfall liegt oder in den einen kleinen Kontrollfall fällt, der separat per direkter Teilerenumeration behandelt wird.
Im Einheitsfall erstellt die Implementierung die Liste aller zu \(M\) teilerfremden Reste, ordnet jedem Rest eine Tabellenposition zu und initialisiert eine eindimensionale DP-Tabelle mit dem Wert \(1\) beim Rest \(1\).
Für jede Primzahl \(p\ne 2,5\) wird die Zyklendarstellung der Permutation \(r\mapsto pr \pmod{M}\) wiederverwendet. Innerhalb jedes Zyklus wird der konstante Vollzyklusbeitrag einmal berechnet; der verbleibende Teil wird dann mit einem Sliding Window fortgeschrieben, sodass jeder Zustand amortisiert in konstanter Zeit aktualisiert wird.
Weil der Ergebnismodulus etwas größer als \(10^{16}\) ist, verwenden die Implementierungen außerdem überlaufsichere modulare Arithmetik. Nach der Verarbeitung aller ungeraden Primzahlen ist der Eintrag zum Zielrest \(t=d/f\) genau die gesuchte Anzahl modulo \(10^{16}+61\).
Sei \(U=\varphi(M)\), und sei \(C\) die Anzahl der verschiedenen Restklassen \(p\bmod M\), die bei Primzahlen \(p\le n\) mit \(p\ne 2,5\) auftreten. Das Primzahlsieb benötigt \(O(n\log\log n)\) Zeit und \(O(n)\) Speicher. Das Cachen der Zykluslayouts kostet in der aktuellen Implementierung \(O(CU)\) Zeit und \(O(CU)\) Speicher. Sind diese Layouts vorhanden, berührt jeder Primzahl-Übergang jeden Einheitsrest genau einmal, also beträgt die eigentliche DP-Arbeit \(O(\pi(n)\,U)\).
Für die echte Euler-Eingabe gilt \(M=12500\) und \(U=5000\). Darum bleibt das Verfahren praktikabel, obwohl \(10^6!\) selbst eine unfassbar große Menge von Teilern besitzt.
\(k\), \(d\)'nin ondalık basamak sayısı olsun. Aradığımız değer
$$F(n,d)=\#\left\{x\in \mathbb{Z}_{>0}:x\mid n!,\ x\equiv d \pmod{10^k}\right\}$$
sayısıdır; sonuç \(10^{16}+61\) modunda verilir. Problem 474 için somut girdi \(n=10^6\) ve \(d=65432\)'dir. \(n!\)'ın bölen sayısı aşırı büyük olduğundan çözüm, bölenleri tek tek üretmek yerine asal üsleri ve modüler kalıntıları kullanır.
\(n!\)'ın her böleni, \(p\le n\) olan her asal için bir üs seçilerek belirlenir. Uygulama önce bu gösterimi kullanır, sonra hedef son blok tarafından zorlanan \(2\) ve \(5\) kuvvetlerini ayırır, en sonunda da indirgenmiş bir \(10\) kuvveti modunda birim kalıntılar üzerinde dinamik programlama yapar.
Eğer
$$n!=\prod_{p\le n} p^{e_p},\qquad e_p=v_p(n!)=\sum_{m\ge 1}\left\lfloor\frac{n}{p^m}\right\rfloor,$$
ise her \(x\mid n!\) böleni
$$x=\prod_{p\le n} p^{a_p},\qquad 0\le a_p\le e_p$$
şeklindedir. Dolayısıyla problem, çarpımı \(10^k\) modunda tam istenen kalıntıyı veren üs vektörlerini saymaya dönüşür.
Şöyle yazalım:
$$d=2^{\alpha}5^{\beta}\tau,\qquad \gcd(\tau,10)=1,$$
burada \(\alpha=v_2(d)\) ve \(\beta=v_5(d)\) olsun. Asıl Euler girdisinde olduğu gibi ana durumda hem \(\alpha\lt k\) hem de \(\beta\lt k\) geçerlidir. Şimdi
$$f=2^{\alpha}5^{\beta},\qquad M=\frac{10^k}{f}=2^{k-\alpha}5^{k-\beta},\qquad t=\frac{d}{f}$$
tanımlarını yapalım. Eğer bir bölen \(x\), \(x\equiv d \pmod{10^k}\) koşulunu sağlıyorsa, \(x\)'in içindeki \(2\) sayısı tam olarak \(\alpha\), \(5\) sayısı da tam olarak \(\beta\) olmak zorundadır. Çünkü kongruansı \(f\)'ye böldüğümüzde
$$\frac{x}{f}\equiv t \pmod{M}$$
elde edilir ve \(t\), \(10\) ile aralarında asaldır; yani \(x/f\), \(M\) modunda birimdir. Böylece \(2\) ve \(5\) üsleri sabitlenir, serbest kalanlar sadece \(p\ne 2,5\) asallarıdır.
Eğer \(\alpha>v_2(n!)\) ya da \(\beta>v_5(n!)\) ise cevap hemen \(0\) olur. Uygulamalar ayrıca \((n,d)=(12,12)\) küçük kontrolü için ayrı bir kaba kuvvet dalı tutar; çünkü orada \(\alpha=k\) olur ve yukarıdaki birim indirgemesi doğrudan kullanılamaz.
Sabit parça \(f\) dışarı alındıktan sonra bölenin kalan birim kısmı
$$u=\prod_{\substack{p\le n\\ p\ne 2,5}} p^{a_p}\pmod{M}$$
olur. Yalnızca \(M\) ile aralarında asal olan kalıntılar ortaya çıkabileceğinden durum uzayının büyüklüğü
$$U=\varphi(M)$$
kadardır. \(A(r)\), şimdiye kadar işlenen asallar için kalan kısmı \(r\) veren seçim sayısı olsun. Başlangıçta yalnızca boş çarpım bulunduğu için
$$A(1)=1,\qquad A(r)=0\text{, }r\ne 1$$
olur. \(p\ne 2,5\) olan bir asal ve \(0\le j\le e_p\) aralığındaki üsler için geçiş
$$A_{\mathrm{yeni}}(x)=\sum_{j=0}^{e_p} A_{\mathrm{eski}}(x p^{-j})$$
şeklindedir; burada ters alma \(M\) modunda yapılır. Bunun nedeni, \(p\)'nin üssü \(j\) seçildiğinde mevcut kalıntının \(p^j\) ile çarpılmasıdır.
\(\gcd(p,M)=1\) olduğu için \(p\) ile çarpmak, \(M\) modundaki birim kalıntıları bir permütasyon gibi dolaştırır. Bu yüzden birim kümesi ayrık döngülere bölünür:
$$u_0,\ u_1,\ \dots,\ u_{L-1},\qquad u_{i+1}\equiv p\,u_i \pmod{M}.$$
Tek bir döngü üzerinde geçiş, çevrimsel bir evrişim olur:
$$A_{\mathrm{yeni}}(u_i)=\sum_{j=0}^{e_p} A_{\mathrm{eski}}(u_{i-j}),$$
burada indisler \(L\)'ye göre mod alınır. Şimdi
$$e_p+1=qL+r,\qquad 0\le r\lt L$$
yazalım. O zaman her yeni değer, döngünün tamamının \(q\) kez toplanmasına ve uzunluğu \(r\) olan ek bir pencereye ayrılır:
$$A_{\mathrm{yeni}}(u_i)=q\sum_{m=0}^{L-1}A_{\mathrm{eski}}(u_m)+\sum_{j=0}^{r-1}A_{\mathrm{eski}}(u_{i-j}).$$
İlk terim döngü boyunca sabittir; ikinci terim ise kayan pencere ile güncellenebilir. Böylece tek bir asal için döngü başına maliyet \(O(e_pL)\)'den \(O(L)\)'ye iner.
Tüm \(p\ne 2,5\) asalları işlendiğinde aranan sayı, hedef birim kalıntı \(t\)'deki durumdur:
$$F(n,d)\equiv A(t)\pmod{10^{16}+61}.$$
Yani kalıntı uzayı \(M\) modundadır; fakat DP sayaçları boyunca bütün aritmetik \(10^{16}+61\) modunda tutulur.
Burada \(k=5\) ve
$$65432=2^3\cdot 8179$$
olur. Dolayısıyla \(\alpha=3\), \(\beta=0\), \(f=8\) ve
$$M=\frac{10^5}{8}=12500,\qquad t=\frac{65432}{8}=8179$$
elde edilir. \(\gcd(8179,12500)=1\) olduğundan ana birim-durum algoritması doğrudan uygulanır. Sayılan bölenler tam olarak \(2\) üssü \(3\), \(5\) üssü \(0\) olan ve geriye kalan birim kısmı \(12500\) modunda \(8179\)'a denk gelen bölenlerdir. Böylece dinamik programlama, \(10^6!\)'ın bütün bölenleri yerine yalnızca
$$U=\varphi(12500)=12500\left(1-\frac12\right)\left(1-\frac15\right)=5000$$
adet birim kalıntı üzerinde çalışır.
C++, Python ve Java uygulamaları aynı sayısal stratejiyi izler. Önce \(n\)'ye kadar olan asallar süzülür, ardından Legendre formülü ile \(v_p(n!)\) hesaplanır ve hedef son bloğun genel birim-durum rejimine mi yoksa küçük kaba kuvvet kontrolüne mi girdiği belirlenir.
Birim durumda uygulama, \(M\) modunda \(M\) ile aralarında asal olan bütün kalıntıları listeler, her birini bir dizi konumuna eşler ve tek boyutlu DP tablosunu \(1\) kalıntısında değer \(1\) ile başlatır.
Her \(p\ne 2,5\) asalı için \(r\mapsto pr \pmod{M}\) permütasyonunun döngü ayrışımı yeniden kullanılır. Her döngüde tam döngü katkısı bir kez hesaplanır; kalan kısım kayan pencere ile ilerletilir. Böylece döngü içindeki her durum sabit amortismanlı sürede güncellenir.
Cevap modu \(10^{16}\)'nın biraz üzerinde olduğu için uygulamalar taşmaya dayanıklı modüler aritmetik de kullanır. Tüm tek asallar işlendiğinde \(t=d/f\) kalıntısındaki tablo değeri aranan sonucu verir.
\(U=\varphi(M)\) olsun. Ayrıca \(C\), \(p\le n\) ve \(p\ne 2,5\) için ortaya çıkan farklı \(p\bmod M\) kalıntılarının sayısı olsun. Asal süzgeci \(O(n\log\log n)\) zaman ve \(O(n)\) bellek kullanır. Mevcut uygulamada döngü yerleşimlerini önbelleğe almak \(O(CU)\) zaman ve \(O(CU)\) bellek ister. Bu yerleşimler hazır olduktan sonra her asal geçişi her birim kalıntıya tam bir kez dokunduğu için DP işi \(O(\pi(n)\,U)\) olur.
Gerçek Euler girdisinde \(M=12500\) ve \(U=5000\) olduğundan bu yaklaşım pratiktir; oysa \(10^6!\)'ın bölenlerini doğrudan taramak imkansızdır.
Sea \(k\) el número de cifras decimales de \(d\). Queremos contar
$$F(n,d)=\#\left\{x\in \mathbb{Z}_{>0}:x\mid n!,\ x\equiv d \pmod{10^k}\right\},$$
y devolver el resultado módulo \(10^{16}+61\). En el problema 474 la entrada concreta es \(n=10^6\) y \(d=65432\). Como \(n!\) tiene una cantidad inmensa de divisores, el cálculo debe hacerse con exponentes primos y clases residuales, no enumerando divisores uno por uno.
Cada divisor de \(n!\) queda determinado al elegir un exponente para cada primo \(p\le n\). La implementación usa exactamente esa descripción, separa después las potencias de \(2\) y \(5\) forzadas por el sufijo objetivo, y finalmente ejecuta una programación dinámica sobre el grupo de unidades módulo una potencia reducida de \(10\).
Si
$$n!=\prod_{p\le n} p^{e_p},\qquad e_p=v_p(n!)=\sum_{m\ge 1}\left\lfloor\frac{n}{p^m}\right\rfloor,$$
entonces todo divisor \(x\mid n!\) puede escribirse como
$$x=\prod_{p\le n} p^{a_p},\qquad 0\le a_p\le e_p.$$
Así, el problema consiste en contar los vectores de exponentes \((a_p)\) cuyo producto cae en una clase residual concreta módulo \(10^k\).
Escribimos
$$d=2^{\alpha}5^{\beta}\tau,\qquad \gcd(\tau,10)=1,$$
donde \(\alpha=v_2(d)\) y \(\beta=v_5(d)\). En el régimen principal, que es precisamente el del dato real de Euler, se cumple \(\alpha\lt k\) y \(\beta\lt k\). Definimos
$$f=2^{\alpha}5^{\beta},\qquad M=\frac{10^k}{f}=2^{k-\alpha}5^{k-\beta},\qquad t=\frac{d}{f}.$$
Si un divisor \(x\) satisface \(x\equiv d \pmod{10^k}\), entonces \(x\) debe contener exactamente \(\alpha\) factores \(2\) y exactamente \(\beta\) factores \(5\). Al dividir la congruencia por \(f\) obtenemos
$$\frac{x}{f}\equiv t \pmod{M},$$
y como \(t\) es coprimo con \(10\), \(x/f\) es automáticamente una unidad módulo \(M\). Por tanto, los exponentes de \(2\) y \(5\) quedan fijados, y solo permanecen libres los primos \(p\ne 2,5\).
Si \(\alpha>v_2(n!)\) o \(\beta>v_5(n!)\), la respuesta es \(0\) de inmediato. Además, las implementaciones conservan una rama separada de fuerza bruta para el pequeño control \((n,d)=(12,12)\), porque allí \(\alpha=k\) y la reducción a unidades ya no vale.
Después de extraer la parte fija \(f\), la parte unitaria restante de un divisor es
$$u=\prod_{\substack{p\le n\\ p\ne 2,5}} p^{a_p}\pmod{M}.$$
Solo pueden aparecer residuos coprimos con \(M\), así que el espacio de estados tiene tamaño
$$U=\varphi(M).$$
Sea \(A(r)\) el número de elecciones de exponentes para los primos ya procesados que producen el residuo \(r\) módulo \(M\). Al principio solo existe el producto vacío, de modo que
$$A(1)=1,\qquad A(r)=0\text{ para }r\ne 1.$$
Al procesar un primo \(p\ne 2,5\) con rango de exponentes \(0\le j\le e_p\), la transición es
$$A_{\mathrm{nuevo}}(x)=\sum_{j=0}^{e_p} A_{\mathrm{antiguo}}(x p^{-j}),$$
donde el inverso se toma módulo \(M\). Esto es correcto porque escoger el exponente \(j\) para \(p\) multiplica el residuo actual por \(p^j\).
Como \(\gcd(p,M)=1\), multiplicar por \(p\) permuta las unidades módulo \(M\). Por eso el conjunto de unidades se divide en ciclos disjuntos
$$u_0,\ u_1,\ \dots,\ u_{L-1},\qquad u_{i+1}\equiv p\,u_i \pmod{M}.$$
Sobre un ciclo, la transición se convierte en una convolución cíclica:
$$A_{\mathrm{nuevo}}(u_i)=\sum_{j=0}^{e_p} A_{\mathrm{antiguo}}(u_{i-j}),$$
con índices tomados módulo \(L\). Escribimos
$$e_p+1=qL+r,\qquad 0\le r\lt L.$$
Entonces cada valor nuevo contiene \(q\) vueltas completas del ciclo más una ventana extra de longitud \(r\):
$$A_{\mathrm{nuevo}}(u_i)=q\sum_{m=0}^{L-1}A_{\mathrm{antiguo}}(u_m)+\sum_{j=0}^{r-1}A_{\mathrm{antiguo}}(u_{i-j}).$$
El primer término es constante en todo el ciclo, y el segundo puede actualizarse con una ventana deslizante. Así, el coste por ciclo baja de \(O(e_pL)\) a \(O(L)\).
Después de procesar todos los primos \(p\ne 2,5\), la cantidad buscada es simplemente el estado en el residuo objetivo \(t\):
$$F(n,d)\equiv A(t)\pmod{10^{16}+61}.$$
Las clases residuales viven módulo \(M\), mientras que la aritmética de conteo de la programación dinámica se lleva siempre módulo \(10^{16}+61\).
Aquí \(k=5\) y
$$65432=2^3\cdot 8179.$$
Por tanto \(\alpha=3\), \(\beta=0\), \(f=8\) y
$$M=\frac{10^5}{8}=12500,\qquad t=\frac{65432}{8}=8179.$$
Como \(\gcd(8179,12500)=1\), el algoritmo principal de unidades se aplica directamente. Los divisores contados son exactamente aquellos divisores de \(n!\) cuyo exponente de \(2\) es \(3\), cuyo exponente de \(5\) es \(0\) y cuya parte unitaria restante es congruente con \(8179\) módulo \(12500\). Así, la programación dinámica trabaja sobre
$$U=\varphi(12500)=12500\left(1-\frac12\right)\left(1-\frac15\right)=5000$$
residuos unitarios en lugar de sobre el conjunto gigantesco de todos los divisores de \(10^6!\).
Las implementaciones en C++, Python y Java siguen la misma estrategia numérica. Primero generan todos los primos hasta \(n\), calculan \(v_p(n!)\) mediante la fórmula de Legendre y deciden si el sufijo buscado cae en el caso general de unidades o en el único control pequeño que se trata aparte por enumeración directa de divisores.
En el caso unitario, la implementación construye la lista de residuos módulo \(M\) coprimos con \(M\), asigna a cada uno una posición en un arreglo e inicializa una tabla de programación dinámica unidimensional con valor \(1\) en el residuo \(1\).
Para cada primo \(p\ne 2,5\), reutiliza la descomposición en ciclos de la permutación \(r\mapsto pr \pmod{M}\). Dentro de cada ciclo calcula una vez la contribución de las vueltas completas y luego actualiza la parte restante con una ventana deslizante, de modo que cada estado del ciclo se procesa en tiempo amortizado constante.
Como el módulo de la respuesta está ligeramente por encima de \(10^{16}\), las implementaciones también usan aritmética modular segura frente a desbordamientos. Tras procesar todos los primos impares, la entrada correspondiente a \(t=d/f\) es exactamente la cuenta buscada módulo \(10^{16}+61\).
Sea \(U=\varphi(M)\), y sea \(C\) el número de clases residuales distintas \(p\bmod M\) que aparecen entre los primos \(p\le n\) con \(p\ne 2,5\). La criba de primos cuesta \(O(n\log\log n)\) tiempo y \(O(n)\) memoria. En la implementación actual, cachear los diseños de ciclos cuesta \(O(CU)\) tiempo y \(O(CU)\) memoria. Una vez disponibles, cada transición prima toca cada residuo unitario exactamente una vez, así que el trabajo de programación dinámica es \(O(\pi(n)\,U)\).
Para la entrada real de Euler se tiene \(M=12500\) y \(U=5000\), lo que hace viable el método aunque \(10^6!\) tenga un número astronómico de divisores.
设 \(k\) 是 \(d\) 的十进制位数。我们要求
$$F(n,d)=\#\left\{x\in \mathbb{Z}_{>0}:x\mid n!,\ x\equiv d \pmod{10^k}\right\},$$
并将结果对 \(10^{16}+61\) 取模。对 Project Euler 474 而言，具体输入是 \(n=10^6\) 与 \(d=65432\)。由于 \(n!\) 的因子数量极其巨大，算法不可能显式枚举所有因子，而必须完全依靠素因子指数和模意义下的剩余类来计数。
\(n!\) 的每个因子都对应于对每个素数 \(p\le n\) 选择一个指数。实现正是从这个表示出发，先把目标后缀强制规定的 \(2\) 和 \(5\) 的幂次剥离出来，再在一个较小模数下的单位群上做动态规划。
如果
$$n!=\prod_{p\le n} p^{e_p},\qquad e_p=v_p(n!)=\sum_{m\ge 1}\left\lfloor\frac{n}{p^m}\right\rfloor,$$
那么任意因子 \(x\mid n!\) 都可写成
$$x=\prod_{p\le n} p^{a_p},\qquad 0\le a_p\le e_p.$$
因此，原问题等价于统计所有指数向量 \((a_p)\) 中，有多少个会使对应的因子落入模 \(10^k\) 的某个指定剩余类。
把目标后缀写成
$$d=2^{\alpha}5^{\beta}\tau,\qquad \gcd(\tau,10)=1,$$
其中 \(\alpha=v_2(d)\)，\(\beta=v_5(d)\)。对本题真正需要求解的输入，满足 \(\alpha\lt k\) 且 \(\beta\lt k\)。定义
$$f=2^{\alpha}5^{\beta},\qquad M=\frac{10^k}{f}=2^{k-\alpha}5^{k-\beta},\qquad t=\frac{d}{f}.$$
若某个因子 \(x\) 满足 \(x\equiv d \pmod{10^k}\)，那么 \(x\) 中 \(2\) 的指数必须恰好等于 \(\alpha\)，\(5\) 的指数也必须恰好等于 \(\beta\)。原因是把同余式除以 \(f\) 后得到
$$\frac{x}{f}\equiv t \pmod{M},$$
而 \(t\) 与 \(10\) 互素，所以 \(x/f\) 自动是模 \(M\) 的单位。也就是说，\(2\) 和 \(5\) 的指数在这一阶段已经完全固定，真正还需要选择的只有所有 \(p\ne 2,5\) 的素数指数。
如果 \(\alpha>v_2(n!)\) 或 \(\beta>v_5(n!)\)，那么答案立刻就是 \(0\)。实现中还保留了一个单独的小规模暴力分支来处理检查点 \((n,d)=(12,12)\)，因为在那个例子里 \(\alpha=k\)，上面的单位化简不再适用。
把固定部分 \(f\) 提出去以后，一个因子的剩余单位部分为
$$u=\prod_{\substack{p\le n\\ p\ne 2,5}} p^{a_p}\pmod{M}.$$
可能出现的只会是与 \(M\) 互素的剩余类，因此状态总数为
$$U=\varphi(M).$$
令 \(A(r)\) 表示对于已经处理过的素数，能够得到剩余类 \(r\) 的指数选择个数。初始时只有空乘积，因此
$$A(1)=1,\qquad A(r)=0\text{，当 }r\ne 1.$$
现在处理一个素数 \(p\ne 2,5\)，其可选指数范围为 \(0\le j\le e_p\)。转移公式是
$$A_{\mathrm{new}}(x)=\sum_{j=0}^{e_p} A_{\mathrm{old}}(x p^{-j}),$$
其中逆元在模 \(M\) 下理解。这个公式很自然，因为若为 \(p\) 选定指数 \(j\)，现有剩余类就会被乘上 \(p^j\)。
由于 \(\gcd(p,M)=1\)，映射“乘以 \(p\)”会在模 \(M\) 的单位集合上给出一个置换。因此单位集合可以拆成若干个互不相交的循环：
$$u_0,\ u_1,\ \dots,\ u_{L-1},\qquad u_{i+1}\equiv p\,u_i \pmod{M}.$$
在单个循环上，转移变成循环卷积：
$$A_{\mathrm{new}}(u_i)=\sum_{j=0}^{e_p} A_{\mathrm{old}}(u_{i-j}),$$
其中下标对 \(L\) 取模。把
$$e_p+1=qL+r,\qquad 0\le r\lt L$$
写成整除形式后，可以看出每个新状态都由两部分组成：一部分是对整条循环完整绕行 \(q\) 次，另一部分是长度为 \(r\) 的额外窗口，即
$$A_{\mathrm{new}}(u_i)=q\sum_{m=0}^{L-1}A_{\mathrm{old}}(u_m)+\sum_{j=0}^{r-1}A_{\mathrm{old}}(u_{i-j}).$$
第一项在整条循环上是常数，第二项可以用滑动窗口维护。这样一来，一个素数在该循环上的更新时间就从朴素的 \(O(e_pL)\) 降为 \(O(L)\)。
当所有 \(p\ne 2,5\) 的素数都处理完成后，所求数量就是目标单位剩余类 \(t\) 处的状态值：
$$F(n,d)\equiv A(t)\pmod{10^{16}+61}.$$
也就是说，剩余类运算发生在模 \(M\) 下，而动态规划中的计数始终在模 \(10^{16}+61\) 下进行。
这里 \(k=5\)，而且
$$65432=2^3\cdot 8179.$$
因此 \(\alpha=3\)、\(\beta=0\)、\(f=8\)，从而得到
$$M=\frac{10^5}{8}=12500,\qquad t=\frac{65432}{8}=8179.$$
由于 \(\gcd(8179,12500)=1\)，主算法可以直接使用。也就是说，我们真正统计的是那些因子 \(x\mid n!\)，它们满足：\(2\) 的指数恰好为 \(3\)，\(5\) 的指数恰好为 \(0\)，并且去掉这部分以后剩余的单位部分在模 \(12500\) 下同余于 \(8179\)。于是动态规划只需要在
$$U=\varphi(12500)=12500\left(1-\frac12\right)\left(1-\frac15\right)=5000$$
个单位剩余类上运行，而不需要面对 \(10^6!\) 的天文数量级因子。
C++、Python 和 Java 实现遵循同一套数值流程。它们先筛出不超过 \(n\) 的全部素数，再用 Legendre 公式求出每个 \(v_p(n!)\)，随后判断目标后缀属于一般的单位情形，还是属于那个单独保留的小型检查点，需要直接枚举因子。
在单位情形下，实现会先构造所有与 \(M\) 互素的模 \(M\) 剩余类列表，并把每个剩余类映射到数组下标；然后用一个一维动态规划表记录计数，初值仅在剩余类 \(1\) 处为 \(1\)。
对每个 \(p\ne 2,5\) 的素数，实现都会复用置换 \(r\mapsto pr \pmod{M}\) 的循环分解。对于每个循环，先一次性求出完整绕行若干圈的常量贡献，再用滑动窗口维护剩余部分，因此循环中的每个状态都能以均摊常数时间完成更新。
由于答案模数略大于 \(10^{16}\)，实现还采用了防溢出的模加法和模乘法。全部奇素数处理完毕后，剩余类 \(t=d/f\) 对应的表项就是所求答案模 \(10^{16}+61\)。
记 \(U=\varphi(M)\)，再记 \(C\) 为所有满足 \(p\le n\)、\(p\ne 2,5\) 的素数中，不同剩余类 \(p\bmod M\) 的个数。素数筛部分需要 \(O(n\log\log n)\) 时间和 \(O(n)\) 空间。当前实现会缓存循环分解，因此这部分需要 \(O(CU)\) 时间和 \(O(CU)\) 额外空间。循环布局准备好以后，每个素数的转移只会扫描全部单位状态一次，所以动态规划主体的工作量是 \(O(\pi(n)\,U)\)。
对实际 Euler 输入来说，\(M=12500\)，\(U=5000\)。正因为状态空间被压缩到这个量级，算法才可以在 \(10^6!\) 这样庞大的对象上工作。
Пусть \(k\) обозначает число десятичных цифр в \(d\). Требуется найти
$$F(n,d)=\#\left\{x\in \mathbb{Z}_{>0}:x\mid n!,\ x\equiv d \pmod{10^k}\right\},$$
а итоговый результат вывести по модулю \(10^{16}+61\). В задаче 474 конкретные параметры равны \(n=10^6\) и \(d=65432\). Число делителей у \(n!\) колоссально, поэтому перебор самих делителей невозможен; считать нужно через показатели простых и классы вычетов.
Каждый делитель \(n!\) задается выбором одного показателя для каждого простого \(p\le n\). Реализация исходит именно из этого описания, затем отделяет степени \(2\) и \(5\), которые уже принудительно определены нужным десятичным суффиксом, и после этого запускает динамическое программирование по группе единиц modulo уменьшенного модуля.
Если
$$n!=\prod_{p\le n} p^{e_p},\qquad e_p=v_p(n!)=\sum_{m\ge 1}\left\lfloor\frac{n}{p^m}\right\rfloor,$$
то любой делитель \(x\mid n!\) имеет вид
$$x=\prod_{p\le n} p^{a_p},\qquad 0\le a_p\le e_p.$$
Следовательно, задача сводится к подсчету тех векторов показателей \((a_p)\), для которых произведение попадает в один заданный класс вычетов modulo \(10^k\).
Запишем
$$d=2^{\alpha}5^{\beta}\tau,\qquad \gcd(\tau,10)=1,$$
где \(\alpha=v_2(d)\), \(\beta=v_5(d)\). В основном режиме, который и используется для реального входа Euler, выполняется \(\alpha\lt k\) и \(\beta\lt k\). Положим
$$f=2^{\alpha}5^{\beta},\qquad M=\frac{10^k}{f}=2^{k-\alpha}5^{k-\beta},\qquad t=\frac{d}{f}.$$
Если делитель \(x\) удовлетворяет \(x\equiv d \pmod{10^k}\), то в \(x\) обязательно содержится ровно \(\alpha\) двоек и ровно \(\beta\) пятерок. После деления сравнения на \(f\) получаем
$$\frac{x}{f}\equiv t \pmod{M},$$
а число \(t\) взаимно просто с \(10\), поэтому \(x/f\) автоматически является единицей modulo \(M\). Значит, показатели при \(2\) и \(5\) уже полностью фиксированы, и свободный выбор остается только для простых \(p\ne 2,5\).
Если \(\alpha>v_2(n!)\) или \(\beta>v_5(n!)\), ответ сразу равен \(0\). Кроме того, в реализациях есть отдельная крошечная ветка полного перебора для контрольного случая \((n,d)=(12,12)\), потому что там \(\alpha=k\) и приведенная редукция к единицам уже не работает.
После вынесения фиксированного множителя \(f\) единичная часть делителя имеет вид
$$u=\prod_{\substack{p\le n\\ p\ne 2,5}} p^{a_p}\pmod{M}.$$
Появляться могут только остатки, взаимно простые с \(M\), поэтому размер пространства состояний равен
$$U=\varphi(M).$$
Пусть \(A(r)\) обозначает число способов выбрать показатели уже обработанных простых так, чтобы получить остаток \(r\) modulo \(M\). В начале имеется только пустое произведение, поэтому
$$A(1)=1,\qquad A(r)=0\text{ при }r\ne 1.$$
При обработке одного простого \(p\ne 2,5\) с диапазоном показателей \(0\le j\le e_p\) переход записывается так:
$$A_{\mathrm{new}}(x)=\sum_{j=0}^{e_p} A_{\mathrm{old}}(x p^{-j}),$$
где обратный элемент понимается modulo \(M\). Это верно, потому что выбор показателя \(j\) для \(p\) умножает текущий остаток на \(p^j\).
Так как \(\gcd(p,M)=1\), умножение на \(p\) переставляет все единицы modulo \(M\). Поэтому множество единиц распадается на непересекающиеся циклы
$$u_0,\ u_1,\ \dots,\ u_{L-1},\qquad u_{i+1}\equiv p\,u_i \pmod{M}.$$
На одном цикле переход превращается в циклическую свертку:
$$A_{\mathrm{new}}(u_i)=\sum_{j=0}^{e_p} A_{\mathrm{old}}(u_{i-j}),$$
где индексы читаются по модулю \(L\). Запишем
$$e_p+1=qL+r,\qquad 0\le r\lt L.$$
Тогда каждое новое значение состоит из \(q\) полных обходов цикла и еще одного окна длины \(r\):
$$A_{\mathrm{new}}(u_i)=q\sum_{m=0}^{L-1}A_{\mathrm{old}}(u_m)+\sum_{j=0}^{r-1}A_{\mathrm{old}}(u_{i-j}).$$
Первая часть постоянна на всем цикле, а вторую удобно поддерживать скользящим окном. За счет этого стоимость обновления на цикле падает с \(O(e_pL)\) до \(O(L)\).
После обработки всех простых \(p\ne 2,5\) искомое количество равно значению состояния в целевом остатке \(t\):
$$F(n,d)\equiv A(t)\pmod{10^{16}+61}.$$
То есть сами остатки живут modulo \(M\), а счетчики динамического программирования все время ведутся modulo \(10^{16}+61\).
Здесь \(k=5\), причем
$$65432=2^3\cdot 8179.$$
Значит \(\alpha=3\), \(\beta=0\), \(f=8\), а потому
$$M=\frac{10^5}{8}=12500,\qquad t=\frac{65432}{8}=8179.$$
Поскольку \(\gcd(8179,12500)=1\), основной алгоритм применяется напрямую. Мы считаем именно те делители \(n!\), у которых показатель при \(2\) равен \(3\), показатель при \(5\) равен \(0\), а оставшаяся единичная часть сравнима с \(8179\) modulo \(12500\). Поэтому динамическое программирование идет только по
$$U=\varphi(12500)=12500\left(1-\frac12\right)\left(1-\frac15\right)=5000$$
обратимым остаткам вместо астрономического множества всех делителей числа \(10^6!\).
Реализации на C++, Python и Java используют одну и ту же вычислительную схему. Сначала они строят список простых до \(n\), затем по формуле Лежандра находят все значения \(v_p(n!)\), после чего определяют, попадает ли нужный суффикс в общий единичный режим или в единственный маленький контрольный случай, где допустим прямой перебор делителей.
В единичном режиме реализация строит список всех остатков modulo \(M\), взаимно простых с \(M\), сопоставляет каждому остатку индекс массива и инициализирует одномерную таблицу динамического программирования значением \(1\) в остатке \(1\).
Для каждого простого \(p\ne 2,5\) повторно используется разложение перестановки \(r\mapsto pr \pmod{M}\) на циклы. На каждом цикле сначала вычисляется вклад полных обходов, а затем остаточная часть обновляется скользящим окном, поэтому каждое состояние обрабатывается за амортизированное постоянное время.
Так как модуль ответа немного больше \(10^{16}\), реализации также используют безопасную по переполнению модульную арифметику. После обработки всех нечетных простых значение в позиции \(t=d/f\) и есть искомый ответ modulo \(10^{16}+61\).
Пусть \(U=\varphi(M)\), а \(C\) обозначает число различных классов \(p\bmod M\), встречающихся среди простых \(p\le n\) с \(p\ne 2,5\). Решето простых требует \(O(n\log\log n)\) времени и \(O(n)\) памяти. Кэширование циклических разложений в текущей реализации требует \(O(CU)\) времени и \(O(CU)\) памяти. После этого переход для каждого простого проходит по каждому обратимому остатку ровно один раз, так что основная DP-часть работает за \(O(\pi(n)\,U)\).
Для настоящего входа Euler имеем \(M=12500\) и \(U=5000\). Именно это делает метод практичным, несмотря на чудовищное количество делителей у \(10^6!\).
ليكن \(k\) هو عدد الخانات العشرية في \(d\). نريد حساب
$$F(n,d)=\#\left\{x\in \mathbb{Z}_{>0}:x\mid n!,\ x\equiv d \pmod{10^k}\right\},$$
مع إرجاع النتيجة بترديد \(10^{16}+61\). في مسألة Euler 474 تكون القيم الفعلية \(n=10^6\) و\(d=65432\). وبما أن \(n!\) يملك عددًا هائلًا جدًا من القواسم، فلا يمكن تعدادها مباشرة، ولذلك يعتمد الحل على أسس العوامل الأولية وعلى البواقي المعيارية فقط.
كل قاسم لـ \(n!\) يتحدد باختيار أس واحد لكل عدد أولي \(p\le n\). تنطلق الشيفرة من هذا الوصف، ثم تفصل قوى \(2\) و\(5\) التي يفرضها اللاحق العشري المطلوب، وبعد ذلك تجري برمجة ديناميكية على مجموعة الوحدات modulo قوة مختزلة من \(10\).
إذا كان
$$n!=\prod_{p\le n} p^{e_p},\qquad e_p=v_p(n!)=\sum_{m\ge 1}\left\lfloor\frac{n}{p^m}\right\rfloor,$$
فإن كل قاسم \(x\mid n!\) يمكن كتابته على الصورة
$$x=\prod_{p\le n} p^{a_p},\qquad 0\le a_p\le e_p.$$
وبالتالي تصبح المسألة هي عدّ جميع متجهات الأسس \((a_p)\) التي تجعل القاسم الناتج يقع في فئة باقٍ محددة modulo \(10^k\).
نكتب
$$d=2^{\alpha}5^{\beta}\tau,\qquad \gcd(\tau,10)=1,$$
حيث \(\alpha=v_2(d)\) و\(\beta=v_5(d)\). في الحالة الرئيسية المستخدمة في الإدخال الحقيقي للمسألة يكون كل من \(\alpha\lt k\) و\(\beta\lt k\). نعرّف
$$f=2^{\alpha}5^{\beta},\qquad M=\frac{10^k}{f}=2^{k-\alpha}5^{k-\beta},\qquad t=\frac{d}{f}.$$
إذا حقق قاسم \(x\) العلاقة \(x\equiv d \pmod{10^k}\)، فلا بد أن يحتوي \(x\) على \(\alpha\) عوامل من \(2\) بالضبط وعلى \(\beta\) عوامل من \(5\) بالضبط. وبعد قسمة هذه العلاقة على \(f\) نحصل على
$$\frac{x}{f}\equiv t \pmod{M},$$
وبما أن \(t\) أولي نسبيًا مع \(10\)، فإن \(x/f\) يكون تلقائيًا عنصرًا قابلاً للعكس modulo \(M\). وهكذا تصبح أسس \(2\) و\(5\) ثابتة تمامًا، ولا يبقى حرًا إلا أس كل أولي \(p\ne 2,5\).
إذا كان \(\alpha>v_2(n!)\) أو \(\beta>v_5(n!)\)، فالجواب يساوي \(0\) مباشرة. كما تحتفظ التطبيقات بفرع صغير منفصل للعدّ المباشر في حالة الفحص \((n,d)=(12,12)\)، لأن \(\alpha=k\) هناك، وعندها لا يعود اختزال الوحدات صالحًا.
بعد إخراج الجزء الثابت \(f\)، يصبح الجزء الوحدوي المتبقي من القاسم هو
$$u=\prod_{\substack{p\le n\\ p\ne 2,5}} p^{a_p}\pmod{M}.$$
ولا يمكن أن تظهر إلا البواقي الأولية نسبيًا مع \(M\)، ولذلك يكون حجم فضاء الحالات
$$U=\varphi(M).$$
لتكن \(A(r)\) عدد اختيارات الأسس للأوليات التي تمت معالجتها حتى الآن والتي تعطي الباقي \(r\) modulo \(M\). في البداية لا يوجد إلا حاصل الضرب الفارغ، لذا
$$A(1)=1,\qquad A(r)=0\quad (r\ne 1).$$
وعند معالجة أولي \(p\ne 2,5\) مع مجال أسس \(0\le j\le e_p\)، تكون معادلة الانتقال
$$A_{\mathrm{new}}(x)=\sum_{j=0}^{e_p} A_{\mathrm{old}}(x p^{-j}),$$
حيث يؤخذ المعكوس modulo \(M\). وهذا صحيح لأن اختيار الأس \(j\) لهذا الأولي يعني ضرب الباقي الحالي في \(p^j\).
بما أن \(\gcd(p,M)=1\)، فإن الضرب في \(p\) يعطي تبديلًا على مجموعة الوحدات modulo \(M\). لذلك تنقسم هذه المجموعة إلى دورات منفصلة
$$u_0,\ u_1,\ \dots,\ u_{L-1},\qquad u_{i+1}\equiv p\,u_i \pmod{M}.$$
وعلى دورة واحدة يتحول الانتقال إلى التفاف دوري:
$$A_{\mathrm{new}}(u_i)=\sum_{j=0}^{e_p} A_{\mathrm{old}}(u_{i-j}),$$
مع قراءة المؤشرات modulo \(L\). وإذا كتبنا
$$e_p+1=qL+r,\qquad 0\le r\lt L,$$
فإن كل قيمة جديدة تتكون من \(q\) لفات كاملة حول الدورة، بالإضافة إلى نافذة إضافية طولها \(r\):
$$A_{\mathrm{new}}(u_i)=q\sum_{m=0}^{L-1}A_{\mathrm{old}}(u_m)+\sum_{j=0}^{r-1}A_{\mathrm{old}}(u_{i-j}).$$
الحد الأول ثابت على كامل الدورة، والحد الثاني يمكن تحديثه بنافذة منزلقة. وبهذا ينخفض العمل على الدورة الواحدة من \(O(e_pL)\) إلى \(O(L)\).
بعد معالجة جميع الأوليات \(p\ne 2,5\)، يكون العدد المطلوب ببساطة هو قيمة الحالة عند الباقي الهدف \(t\):
$$F(n,d)\equiv A(t)\pmod{10^{16}+61}.$$
أي أن فضاء البواقي نفسه يعمل modulo \(M\)، بينما تُحفظ أعداد الطرق كلها modulo \(10^{16}+61\).
لدينا هنا \(k=5\)، كما أن
$$65432=2^3\cdot 8179.$$
إذًا \(\alpha=3\) و\(\beta=0\) و\(f=8\)، ومن ثم
$$M=\frac{10^5}{8}=12500,\qquad t=\frac{65432}{8}=8179.$$
ولأن \(\gcd(8179,12500)=1\)، فإن الخوارزمية الرئيسية تنطبق مباشرة. وهذا يعني أن القواسم المعدودة هي بالضبط القواسم التي يساوي فيها أس \(2\) العدد \(3\)، وأُس \(5\) العدد \(0\)، ويكون الجزء الوحدوي الباقي فيها مطابقًا لـ \(8179\) modulo \(12500\). لذلك تعمل البرمجة الديناميكية على
$$U=\varphi(12500)=12500\left(1-\frac12\right)\left(1-\frac15\right)=5000$$
باقٍ وحدوي فقط، بدلًا من التعامل مع جميع قواسم \(10^6!\) ذات العدد الفلكي.
تتبع تطبيقات C++ وPython وJava الاستراتيجية العددية نفسها. فهي تبدأ بتوليد جميع الأوليات حتى \(n\)، ثم تحسب \(v_p(n!)\) بصيغة ليجاندر، وبعد ذلك تحدد هل يقع اللاحق المطلوب في الحالة العامة للوحدات أم في حالة الفحص الصغيرة الوحيدة التي تعالج بتعداد مباشر للقواسم.
في حالة الوحدات، تبني الشيفرة قائمة بكل البواقي modulo \(M\) التي تكون أولية نسبيًا مع \(M\)، وتربط كل باقٍ بموضع في مصفوفة، ثم تهيئ جدول برمجة ديناميكية أحادي البعد بحيث تكون القيمة \(1\) عند الباقي \(1\).
ولكل أولي \(p\ne 2,5\)، تعيد الشيفرة استخدام تفكيك التبديل \(r\mapsto pr \pmod{M}\) إلى دورات. داخل كل دورة يُحسب أولًا إسهام اللفات الكاملة، ثم تُحدّث البقية بنافذة منزلقة، ولذلك تُحدّث كل حالة في الدورة بزمن ثابت وسطي.
ولأن موديل الجواب أكبر قليلًا من \(10^{16}\)، تستخدم التطبيقات أيضًا حسابًا معياريًا آمنًا من ناحية الفيض. وبعد الانتهاء من كل الأوليات الفردية، تكون القيمة الموجودة عند \(t=d/f\) هي الجواب المطلوب modulo \(10^{16}+61\).
ليكن \(U=\varphi(M)\)، وليكن \(C\) عدد البواقي المختلفة \(p\bmod M\) التي تظهر بين الأوليات \(p\le n\) مع استثناء \(2\) و\(5\). يكلف منخل الأوليات \(O(n\log\log n)\) زمنًا و\(O(n)\) ذاكرة. وفي التطبيق الحالي يكلف تخزين تفكيكات الدورات مؤقتًا \(O(CU)\) زمنًا و\(O(CU)\) ذاكرة. وبعد تجهيز هذه التفكيكات، فإن كل انتقال خاص بأولي واحد يزور كل باقٍ وحدوي مرة واحدة فقط، لذا يكون عمل البرمجة الديناميكية \(O(\pi(n)\,U)\).
في الإدخال الحقيقي للمسألة نجد أن \(M=12500\) و\(U=5000\)، ولهذا يبقى الأسلوب عمليًا رغم أن عدد قواسم \(10^6!\) هائل جدًا.