The quantity behind Problem 956 can be collapsed to a single integer
$$N_n=\prod_{t=1}^{n} t^{c_t},\qquad c_t=\binom{n-t+2}{2}=\frac{(n-t+1)(n-t+2)}{2}.$$
The required value is the sum of those divisors \(d\mid N_n\) for which the total exponent count
$$\Omega(d)=\sum_p \alpha_p$$
is divisible by \(m\), where \(d=\prod_p p^{\alpha_p}\). The implementations evaluate the case \(n=m=1000\) and reduce the answer modulo \(999999001\).
A brute-force divisor enumeration is impossible. Even after prime factorization, the exponent bounds \(E_p\) of \(N_{1000}\) are enormous, so the number of divisors is astronomically large. The workable approach is to separate the contribution of each prime and track only the residue of the exponent sum modulo \(m\).
The entire solution is a prime-by-prime generating-function computation. The key is that the condition \(\Omega(d)\equiv 0 \pmod m\) depends only on exponent residues, while the divisor value \(d\) factors multiplicatively over primes.
The code never expands the nested product directly. Instead it records how many times each integer \(t\) appears after everything is flattened. That multiplicity is the triangular number \(c_t\), so the whole object is exactly \(N_n=\prod t^{c_t}\).
For any prime \(p\), the exponent of \(p\) in \(N_n\) is therefore
$$E_p=v_p(N_n)=\sum_{t=1}^{n} c_t\,v_p(t).$$
Once all \(E_p\) are known, the original problem has been reduced to a divisor-sum problem on the prime factorization
$$N_n=\prod_p p^{E_p}.$$
Every divisor of \(N_n\) has the form
$$d=\prod_p p^{\alpha_p},\qquad 0\le \alpha_p\le E_p.$$
The quantity being summed is the value of the divisor itself, not just the number of admissible divisors. So a choice of exponent vector \((\alpha_p)\) contributes
$$\prod_p p^{\alpha_p}.$$
The only global constraint is
$$\Omega(d)=\sum_p \alpha_p\equiv 0 \pmod m.$$
This is exactly the kind of condition that can be handled by a residue-class convolution.
For a fixed prime \(p\), introduce
$$F_p(x)=\sum_{a=0}^{E_p} p^a x^a.$$
Multiplying these factors over all primes gives
$$F(x)=\prod_p F_p(x).$$
The coefficient of \(x^k\) in \(F(x)\) is precisely the sum of all divisor values \(d\mid N_n\) with \(\Omega(d)=k\), because choosing the term \(p^{\alpha_p}x^{\alpha_p}\) from each prime contributes the divisor value and records the total exponent count in the power of \(x\).
We do not need the exact exponent total \(k\); only its residue modulo \(m\) matters. So we can work in the quotient ring where \(x^m=1\). In that setting, each prime contributes only \(m\) residue classes:
$$F_p(x)\equiv \sum_{r=0}^{m-1} C_{p,r}x^r \pmod{x^m-1},$$
with
$$C_{p,r}=\sum_{\substack{0\le a\le E_p\\ a\equiv r \pmod m}} p^a.$$
The desired answer is then the coefficient of \(x^0\) in the product of these reduced polynomials.
If \(E_p<r\), then \(C_{p,r}=0\). Otherwise the admissible exponents are
$$a=r,\ r+m,\ r+2m,\ \dots,\ r+(N_{p,r}-1)m,$$
where
$$N_{p,r}=\left\lfloor\frac{E_p-r}{m}\right\rfloor+1.$$
Let \(q=p^m\). Then
$$C_{p,r}=p^r\sum_{j=0}^{N_{p,r}-1} q^j.$$
Modulo \(M=999999001\), this becomes a geometric series. If \(q\not\equiv 1 \pmod M\), then
$$C_{p,r}\equiv p^r\,(q^{N_{p,r}}-1)\,(q-1)^{-1}\pmod M.$$
If \(q\equiv 1 \pmod M\), the denominator would vanish, and the sum collapses to
$$C_{p,r}\equiv p^r\,N_{p,r}\pmod M.$$
This is the only subtle modular case in the entire computation.
Here
$$c_1=6,\qquad c_2=3,\qquad c_3=1,$$
so
$$N_3=1^6\cdot 2^3\cdot 3=24=2^3\cdot 3.$$
We want divisors whose total exponent count is a multiple of \(3\). For the prime \(2\) with exponent \(3\), the residue-class sums are
$$C_{2,0}=1+2^3=9,\qquad C_{2,1}=2,\qquad C_{2,2}=4.$$
For the prime \(3\) with exponent \(1\), they are
$$C_{3,0}=1,\qquad C_{3,1}=3,\qquad C_{3,2}=0.$$
To obtain total residue \(0\pmod 3\), we combine residue pairs \((0,0)\), \((2,1)\), and \((1,2)\). Therefore the required sum is
$$9\cdot 1+4\cdot 3+2\cdot 0=21.$$
Indeed the valid divisors are \(1\), \(8\), and \(12\), and their sum is \(21\). This small case is exactly the same convolution that the full solution performs for \(n=m=1000\).
After reducing exponents modulo \(m\), the coefficient vector for each prime has length \(m\). Multiplying in another prime means a cyclic convolution on those \(m\) residues. No larger state space is needed, because the problem never distinguishes between exponent totals that differ by a multiple of \(m\).
The C++, Python, and Java implementations first build a smallest-prime-factor table up to \(n\). Then they scan \(t=1,2,\dots,n\), factor each \(t\), and add \(c_t\,v_p(t)\) to the stored exponent of every prime \(p\) appearing in \(t\). At the end of this pass they know every nonzero \(E_p\) in the factorization of \(N_n\).
For each prime \(p\), the implementation computes the \(m\) numbers \(C_{p,0},C_{p,1},\dots,C_{p,m-1}\) modulo \(M\). Instead of summing powers one by one, it maintains the current factor \(p^r\) and evaluates the remaining geometric series in closed form. The branch \(p^m\equiv 1\pmod M\) is handled separately so that no invalid modular inverse is taken.
After processing any subset of the primes, the state vector satisfies the invariant
The state \(\mathrm{dp}[r]\) is the sum of divisor values formed from the processed primes whose total exponent residue is \(r\).
When a new prime is incorporated, the update is the cyclic convolution
$$\mathrm{next}[(r+t)\bmod m]\;{+}{=}\;\mathrm{dp}[r]\cdot C_{p,t}\pmod M.$$
After the final prime has been processed, \(\mathrm{dp}[0]\) is exactly the required answer. The C++ and Java versions also compare the modular routine against a tiny exact computation on a small sample, which is a good sanity check for the formulae.
The sieve used to build smallest prime factors is linear in \(n\). Factoring every \(t\le n\) by repeated division through that table is close to linear for these input sizes and is much cheaper than the final convolution phase.
The dominant cost is the residue DP. If \(\pi(n)\) denotes the number of primes up to \(n\), then the transition over all primes uses \(O(\pi(n)m^2)\) modular operations, because each prime contributes an \(m\times m\) cyclic convolution. Memory usage is \(O(n+m)\): the sieve and exponent arrays take \(O(n)\), and the algorithm needs only a constant number of DP buffers of length \(m\).
Die zentrale Größe von Problem 956 lässt sich zu einer einzigen Zahl zusammenfassen:
$$N_n=\prod_{t=1}^{n} t^{c_t},\qquad c_t=\binom{n-t+2}{2}=\frac{(n-t+1)(n-t+2)}{2}.$$
Gesucht ist die Summe aller Teiler \(d\mid N_n\), für die die gesamte Exponentenzahl
$$\Omega(d)=\sum_p \alpha_p$$
durch \(m\) teilbar ist, wenn \(d=\prod_p p^{\alpha_p}\) geschrieben wird. Die Implementierungen berechnen den Fall \(n=m=1000\) modulo \(999999001\).
Eine direkte Aufzählung der Teiler ist hoffnungslos. Schon nach der Primfaktorzerlegung sind die Exponenten \(E_p\) von \(N_{1000}\) riesig, also ist auch die Anzahl der Teiler astronomisch. Der gangbare Weg besteht darin, die Beiträge primeweise zu trennen und nur die Restklasse der gesamten Exponentensumme modulo \(m\) zu verfolgen.
Die ganze Lösung ist eine primeweise Generating-Function-Rechnung. Entscheidend ist, dass die Bedingung \(\Omega(d)\equiv 0 \pmod m\) nur von Restklassen der Exponenten abhängt, während der Wert des Teilers multiplikativ über die Primzahlen zerfällt.
Der Code expandiert das verschachtelte Produkt nie direkt. Stattdessen wird gezählt, wie oft jede Zahl \(t\) nach dem vollständigen Ausmultiplizieren vorkommt. Diese Häufigkeit ist die Dreieckszahl \(c_t\), also ist das gesamte Objekt genau \(N_n=\prod t^{c_t}\).
Für eine Primzahl \(p\) ist ihr Exponent in \(N_n\) damit
$$E_p=v_p(N_n)=\sum_{t=1}^{n} c_t\,v_p(t).$$
Sobald alle \(E_p\) bekannt sind, ist das ursprüngliche Problem auf eine Teilersumme über der Primfaktorzerlegung
$$N_n=\prod_p p^{E_p}$$
reduziert.
Jeder Teiler von \(N_n\) hat die Form
$$d=\prod_p p^{\alpha_p},\qquad 0\le \alpha_p\le E_p.$$
Summiert wird der Teilerwert selbst, nicht nur die Anzahl der zulässigen Teiler. Ein Exponentenvektor \((\alpha_p)\) trägt also
$$\prod_p p^{\alpha_p}$$
zur Summe bei. Die einzige globale Nebenbedingung lautet
$$\Omega(d)=\sum_p \alpha_p\equiv 0 \pmod m.$$
Genau eine solche Bedingung lässt sich mit einer Faltung über Restklassen effizient behandeln.
Für jede Primzahl \(p\) definieren wir
$$F_p(x)=\sum_{a=0}^{E_p} p^a x^a.$$
Das Produkt über alle Primzahlen ist
$$F(x)=\prod_p F_p(x).$$
Der Koeffizient von \(x^k\) in \(F(x)\) ist genau die Summe aller Teilerwerte \(d\mid N_n\) mit \(\Omega(d)=k\), denn aus jedem Primfaktor wählt man einen Term \(p^{\alpha_p}x^{\alpha_p}\), der zugleich den Teilerwert und die gesamte Exponentenzahl aufzeichnet.
Wir brauchen aber nicht den exakten Wert \(k\), sondern nur \(k\bmod m\). Deshalb arbeiten wir modulo \(x^m-1\). Dann liefert jede Primzahl nur noch \(m\) Restklassen:
$$F_p(x)\equiv \sum_{r=0}^{m-1} C_{p,r}x^r \pmod{x^m-1},$$
wobei
$$C_{p,r}=\sum_{\substack{0\le a\le E_p\\ a\equiv r \pmod m}} p^a.$$
Die gesuchte Antwort ist also der Koeffizient von \(x^0\) im Produkt dieser reduzierten Polynome.
Falls \(E_p<r\), gilt \(C_{p,r}=0\). Andernfalls sind die zulässigen Exponenten
$$a=r,\ r+m,\ r+2m,\ \dots,\ r+(N_{p,r}-1)m,$$
mit
$$N_{p,r}=\left\lfloor\frac{E_p-r}{m}\right\rfloor+1.$$
Setzt man \(q=p^m\), dann folgt
$$C_{p,r}=p^r\sum_{j=0}^{N_{p,r}-1} q^j.$$
Modulo \(M=999999001\) ist das eine geometrische Reihe. Für \(q\not\equiv 1 \pmod M\) erhält man
$$C_{p,r}\equiv p^r\,(q^{N_{p,r}}-1)\,(q-1)^{-1}\pmod M.$$
Im Sonderfall \(q\equiv 1 \pmod M\) verschwindet der Nenner, und die Summe vereinfacht sich zu
$$C_{p,r}\equiv p^r\,N_{p,r}\pmod M.$$
Das ist die einzige wirklich heikle Modulo-Stelle der gesamten Rechnung.
Hier gilt
$$c_1=6,\qquad c_2=3,\qquad c_3=1,$$
also
$$N_3=1^6\cdot 2^3\cdot 3=24=2^3\cdot 3.$$
Gesucht sind die Teiler, deren gesamte Exponentenzahl ein Vielfaches von \(3\) ist. Für die Primzahl \(2\) mit Exponent \(3\) ergeben sich
$$C_{2,0}=1+2^3=9,\qquad C_{2,1}=2,\qquad C_{2,2}=4.$$
Für die Primzahl \(3\) mit Exponent \(1\) gilt
$$C_{3,0}=1,\qquad C_{3,1}=3,\qquad C_{3,2}=0.$$
Um die Gesamtrestklasse \(0\pmod 3\) zu erhalten, kommen die Paare \((0,0)\), \((2,1)\) und \((1,2)\) infrage. Daher ist die gesuchte Summe
$$9\cdot 1+4\cdot 3+2\cdot 0=21.$$
Tatsächlich sind die gültigen Teiler \(1\), \(8\) und \(12\), und ihre Summe ist \(21\). Der volle Algorithmus für \(n=m=1000\) macht nichts anderes, nur mit viel mehr Primzahlen.
Nach der Reduktion modulo \(m\) hat jede Primzahl nur einen Koeffizientenvektor der Länge \(m\). Das Multiplizieren mit einer weiteren Primzahl ist daher nur eine zyklische Faltung über diese \(m\) Restklassen. Ein größerer Zustandsraum ist unnötig, weil das Problem Exponentensummen, die sich um ein Vielfaches von \(m\) unterscheiden, bewusst nicht auseinanderhält.
Die C++-, Python- und Java-Implementierungen bauen zunächst eine Tabelle des kleinsten Primteilers bis \(n\) auf. Dann laufen sie über \(t=1,2,\dots,n\), faktorisieren jedes \(t\) und addieren \(c_t\,v_p(t)\) zu jedem vorkommenden Primexponenten. Nach diesem Durchlauf sind alle nichtverschwindenden \(E_p\) der Zerlegung von \(N_n\) bekannt.
Für jedes \(p\) berechnet die Implementierung die \(m\) Werte \(C_{p,0},C_{p,1},\dots,C_{p,m-1}\) modulo \(M\). Statt Potenzen einzeln zu addieren, wird der aktuelle Faktor \(p^r\) fortgeschrieben und die restliche geometrische Reihe in geschlossener Form ausgewertet. Der Sonderfall \(p^m\equiv 1\pmod M\) wird separat behandelt, damit kein unzulässiges modulares Inverses entsteht.
Nach der Verarbeitung eines beliebigen Teilmengen von Primzahlen gilt die Invariante
Der Zustand \(\mathrm{dp}[r]\) ist die Summe der Teilerwerte aus den bereits verarbeiteten Primzahlen, deren gesamte Exponentensumme den Rest \(r\) hat.
Beim Hinzufügen einer neuen Primzahl lautet der Übergang
$$\mathrm{next}[(r+t)\bmod m]\;{+}{=}\;\mathrm{dp}[r]\cdot C_{p,t}\pmod M.$$
Nach der letzten Primzahl ist \(\mathrm{dp}[0]\) genau die gesuchte Antwort. Die C++- und Java-Versionen prüfen diese modulare Rechnung zusätzlich an einem kleinen exakten Beispiel, was die Herleitung gut absichert.
Die Konstruktion der Tabelle der kleinsten Primteiler ist mit dem verwendeten linearen Sieb \(O(n)\). Das Faktorisieren aller \(t\le n\) durch wiederholte Division über diese Tabelle ist für die hier vorkommenden Größen nahezu linear und deutlich günstiger als die abschließende Faltungsphase.
Dominant ist die DP über Restklassen. Bezeichnet \(\pi(n)\) die Anzahl der Primzahlen bis \(n\), dann benötigt der Übergang über alle Primzahlen \(O(\pi(n)m^2)\) modulare Operationen, weil jede Primzahl eine zyklische \(m\times m\)-Faltung beiträgt. Der Speicherbedarf ist \(O(n+m)\): Sieb und Exponentenfelder benötigen \(O(n)\), und dazu kommen nur konstant viele DP-Puffer der Länge \(m\).
Problem 956'nın merkezindeki nicelik tek bir sayıya indirgenebilir:
$$N_n=\prod_{t=1}^{n} t^{c_t},\qquad c_t=\binom{n-t+2}{2}=\frac{(n-t+1)(n-t+2)}{2}.$$
İstenen değer, \(d=\prod_p p^{\alpha_p}\) olmak üzere
$$\Omega(d)=\sum_p \alpha_p$$
toplamı \(m\) ile bölünebilen bütün \(d\mid N_n\) bölenlerinin değer toplamıdır. Uygulamalar \(n=m=1000\) durumunu hesaplayıp sonucu \(999999001\) modunda verir.
Bölenleri tek tek üretmek imkansızdır. \(N_{1000}\)'in asal üsleri \(E_p\) çok büyük olduğu için bölen sayısı astronomik boyuttadır. Çalışan fikir, her asalı bağımsız katkılar halinde ele almak ve yalnızca üs toplamının \(m\) modundaki kalanı izlemektir.
Bütün çözüm, asal asal ilerleyen bir üreteç fonksiyonu hesabıdır. Çünkü \(\Omega(d)\equiv 0 \pmod m\) koşulu yalnızca üslerin artık sınıflarına bakar; buna karşılık bölenin değeri olan \(d\) asal çarpanlara göre tamamen çarpımsal ayrışır.
Kod, iç içe geçmiş çarpımı doğrudan açmıyor. Bunun yerine her \(t\) sayısının tüm açılım sonunda kaç kez göründüğünü sayıyor. Bu çokluk üçgensel sayı \(c_t\) olduğu için elimizde tam olarak \(N_n=\prod t^{c_t}\) var.
Dolayısıyla herhangi bir asal \(p\) için \(N_n\) içindeki üs
$$E_p=v_p(N_n)=\sum_{t=1}^{n} c_t\,v_p(t)$$
olur. Tüm \(E_p\) değerleri elde edildiğinde başlangıçtaki problem,
$$N_n=\prod_p p^{E_p}$$
asal çarpan ayrışımı üzerinde tanımlı bir bölen-toplamı problemine dönüşmüş olur.
\(N_n\)'in her böleni
$$d=\prod_p p^{\alpha_p},\qquad 0\le \alpha_p\le E_p$$
şeklindedir. Toplanan şey yalnızca uygun bölen sayısı değildir; bölenin kendisidir. Bu yüzden \((\alpha_p)\) seçiminden gelen katkı
$$\prod_p p^{\alpha_p}$$
olur. Tek küresel koşul ise
$$\Omega(d)=\sum_p \alpha_p\equiv 0 \pmod m$$
şeklindedir. Bu tam olarak artık sınıf konvolüsyonu ile ele alınabilecek türden bir koşuldur.
Sabit bir asal \(p\) için
$$F_p(x)=\sum_{a=0}^{E_p} p^a x^a$$
tanımını yapalım. Tüm asallar üzerinde çarpınca
$$F(x)=\prod_p F_p(x)$$
elde edilir. \(F(x)\) içindeki \(x^k\) katsayısı, \(\Omega(d)=k\) olan bütün bölenlerin değerleri toplamıdır. Çünkü her asal için seçilen \(p^{\alpha_p}x^{\alpha_p}\) terimi hem bölenin sayısal değerini hem de toplam üs sayısını aynı anda taşır.
Fakat bize tam \(k\) değil, yalnızca \(k\bmod m\) gerekir. Bu yüzden \(x^m=1\) olacak şekilde çalışırız. Böylece her asal yalnızca \(m\) adet artık sınıf katkısı üretir:
$$F_p(x)\equiv \sum_{r=0}^{m-1} C_{p,r}x^r \pmod{x^m-1},$$
burada
$$C_{p,r}=\sum_{\substack{0\le a\le E_p\\ a\equiv r \pmod m}} p^a.$$
Aradığımız cevap, bu indirgenmiş polinomların çarpımındaki \(x^0\) katsayısıdır.
Eğer \(E_p<r\) ise \(C_{p,r}=0\) olur. Aksi halde uygun üsler
$$a=r,\ r+m,\ r+2m,\ \dots,\ r+(N_{p,r}-1)m$$
şeklindedir; burada
$$N_{p,r}=\left\lfloor\frac{E_p-r}{m}\right\rfloor+1.$$
\(q=p^m\) yazarsak
$$C_{p,r}=p^r\sum_{j=0}^{N_{p,r}-1} q^j$$
elde edilir. Modül \(M=999999001\) altında bu bir geometrik seri olur. \(q\not\equiv 1 \pmod M\) ise
$$C_{p,r}\equiv p^r\,(q^{N_{p,r}}-1)\,(q-1)^{-1}\pmod M,$$
\(q\equiv 1 \pmod M\) ise payda sıfıra düştüğünden seri uzunluğa indirgenir:
$$C_{p,r}\equiv p^r\,N_{p,r}\pmod M.$$
Hesabın modüler açıdan hassas olan tek noktası budur.
Bu durumda
$$c_1=6,\qquad c_2=3,\qquad c_3=1,$$
dolayısıyla
$$N_3=1^6\cdot 2^3\cdot 3=24=2^3\cdot 3.$$
Toplam üs sayısı \(3\)'ün katı olan bölenleri istiyoruz. \(2\) asalı için üs \(3\) olduğundan
$$C_{2,0}=1+2^3=9,\qquad C_{2,1}=2,\qquad C_{2,2}=4.$$
\(3\) asalı için ise
$$C_{3,0}=1,\qquad C_{3,1}=3,\qquad C_{3,2}=0.$$
Toplam kalanın \(0\pmod 3\) olması için \((0,0)\), \((2,1)\) ve \((1,2)\) çiftleri kullanılabilir. Bu yüzden istenen toplam
$$9\cdot 1+4\cdot 3+2\cdot 0=21$$
olur. Gerçekten de geçerli bölenler \(1\), \(8\) ve \(12\)'dir; toplamları \(21\) eder. Büyük girdide yapılan işlem de tam olarak bu konvolüsyonun daha büyük bir sürümüdür.
Üsleri \(m\) modunda indirdikten sonra her asalın katkısı uzunluğu \(m\) olan tek bir katsayı vektörüne dönüşür. Yeni bir asal eklemek de bu \(m\) artık sınıf üzerinde döngüsel bir konvolüsyon yapmaktan ibarettir. Daha büyük bir durum uzayı gereksizdir; çünkü problem, toplam üs sayıları arasındaki \(m\) katı farkları zaten ayırt etmez.
C++, Python ve Java uygulamaları önce \(n\)'e kadar en küçük asal bölen tablosu kurar. Sonra \(t=1,2,\dots,n\) üzerinden geçip her \(t\)'yi çarpanlarına ayırır ve görülen her asal için \(c_t\,v_p(t)\) kadar katkıyı ilgili üst toplama ekler. Bu tur bittiğinde \(N_n\)'in asal ayrışımındaki bütün sıfır olmayan \(E_p\) değerleri hazırdır.
Her asal \(p\) için uygulama \(C_{p,0},C_{p,1},\dots,C_{p,m-1}\) değerlerini mod \(M\) altında hesaplar. Kuvvetleri tek tek toplamaktansa, \(p^r\) çarpanı adım adım güncellenir ve geriye kalan geometrik seri kapalı formdan bulunur. \(p^m\equiv 1\pmod M\) özel durumu ayrıca ele alınır; böylece tanımsız bir modüler ters alınmaz.
Asalların herhangi bir alt kümesi işlendiğinde durum vektörü şu değişmezi sağlar:
\(\mathrm{dp}[r]\) durumu, işlenmiş asallardan kurulan ve toplam üs kalanı \(r\) olan bütün bölen değerlerinin toplamını tutar.
Yeni bir asal eklenirken geçiş
$$\mathrm{next}[(r+t)\bmod m]\;{+}{=}\;\mathrm{dp}[r]\cdot C_{p,t}\pmod M$$
şeklindedir. Son asal işlendiğinde \(\mathrm{dp}[0]\) doğrudan cevaptır. C++ ve Java sürümleri ayrıca küçük bir örnekte modüler hesabı tam sayı hesabıyla karşılaştırarak formülün doğru çalıştığını denetler.
En küçük asal bölen tablosu kullanılan doğrusal elek ile \(O(n)\) zamanda kurulur. Bu tablo üzerinden bütün \(t\le n\) değerlerini çarpanlarına ayırmak, bu girdi boyutlarında pratiğe yakın olarak doğrusal davranır ve son DP aşamasından daha ucuzdur.
Baskın maliyet artık-sınıf DP'sidir. \(\pi(n)\), \(n\)'e kadar olan asal sayısını göstersin. Tüm asallar üzerindeki geçişler \(O(\pi(n)m^2)\) modüler işlem gerektirir; çünkü her asal için \(m\times m\) boyutlu döngüsel bir konvolüsyon yapılır. Bellek kullanımı \(O(n+m)\) düzeyindedir: elek ve üs dizileri \(O(n)\) yer kaplar, buna ek olarak uzunluğu \(m\) olan sabit sayıda DP tamponu tutulur.
La cantidad central del Problema 956 puede reducirse a un solo entero:
$$N_n=\prod_{t=1}^{n} t^{c_t},\qquad c_t=\binom{n-t+2}{2}=\frac{(n-t+1)(n-t+2)}{2}.$$
Se pide la suma de todos los divisores \(d\mid N_n\) tales que, si escribimos \(d=\prod_p p^{\alpha_p}\), entonces
$$\Omega(d)=\sum_p \alpha_p$$
sea divisible por \(m\). Las implementaciones evalúan el caso \(n=m=1000\) y reducen el resultado módulo \(999999001\).
Enumerar divisores directamente es inviable. Incluso después de factorizar \(N_{1000}\), los exponentes \(E_p\) son enormes y el número total de divisores es descomunal. La idea útil es separar las contribuciones primo por primo y conservar solo el residuo de la suma total de exponentes módulo \(m\).
Toda la solución es un cálculo con funciones generadoras procesado primo a primo. La condición \(\Omega(d)\equiv 0 \pmod m\) depende solo de residuos de exponentes, mientras que el valor del divisor \(d\) se separa multiplicativamente por primos.
El código no expande el producto anidado de forma explícita. En cambio, cuenta cuántas veces aparece cada entero \(t\) al aplanar toda la expresión. Esa multiplicidad es el número triangular \(c_t\), de modo que el objeto completo es exactamente \(N_n=\prod t^{c_t}\).
Por tanto, para cualquier primo \(p\), su exponente en \(N_n\) es
$$E_p=v_p(N_n)=\sum_{t=1}^{n} c_t\,v_p(t).$$
Una vez conocidos todos los \(E_p\), el problema original se convierte en una suma de divisores sobre la factorización
$$N_n=\prod_p p^{E_p}.$$
Cada divisor de \(N_n\) tiene la forma
$$d=\prod_p p^{\alpha_p},\qquad 0\le \alpha_p\le E_p.$$
Lo que se suma no es el número de divisores válidos, sino el valor del divisor mismo. Así, una elección del vector \((\alpha_p)\) aporta
$$\prod_p p^{\alpha_p}.$$
La única condición global es
$$\Omega(d)=\sum_p \alpha_p\equiv 0 \pmod m.$$
Eso es exactamente el tipo de restricción que puede tratarse con una convolución por clases residuales.
Para un primo fijo \(p\), definimos
$$F_p(x)=\sum_{a=0}^{E_p} p^a x^a.$$
Multiplicando sobre todos los primos obtenemos
$$F(x)=\prod_p F_p(x).$$
El coeficiente de \(x^k\) en \(F(x)\) es precisamente la suma de los valores de todos los divisores \(d\mid N_n\) con \(\Omega(d)=k\). Cada elección del término \(p^{\alpha_p}x^{\alpha_p}\) registra simultáneamente el valor del divisor y la suma total de exponentes.
No necesitamos el valor exacto de \(k\), sino solo \(k\bmod m\). Por eso trabajamos módulo \(x^m-1\). En ese anillo, cada primo aporta solo \(m\) clases residuales:
$$F_p(x)\equiv \sum_{r=0}^{m-1} C_{p,r}x^r \pmod{x^m-1},$$
donde
$$C_{p,r}=\sum_{\substack{0\le a\le E_p\\ a\equiv r \pmod m}} p^a.$$
La respuesta buscada es entonces el coeficiente de \(x^0\) en el producto de estos polinomios reducidos.
Si \(E_p<r\), entonces \(C_{p,r}=0\). En caso contrario, los exponentes admisibles son
$$a=r,\ r+m,\ r+2m,\ \dots,\ r+(N_{p,r}-1)m,$$
con
$$N_{p,r}=\left\lfloor\frac{E_p-r}{m}\right\rfloor+1.$$
Si definimos \(q=p^m\), resulta
$$C_{p,r}=p^r\sum_{j=0}^{N_{p,r}-1} q^j.$$
Módulo \(M=999999001\), esto es una serie geométrica. Si \(q\not\equiv 1 \pmod M\), entonces
$$C_{p,r}\equiv p^r\,(q^{N_{p,r}}-1)\,(q-1)^{-1}\pmod M.$$
Si \(q\equiv 1 \pmod M\), el denominador desaparece y la suma se reduce a
$$C_{p,r}\equiv p^r\,N_{p,r}\pmod M.$$
Ese es el único caso modular delicado de todo el algoritmo.
Aquí tenemos
$$c_1=6,\qquad c_2=3,\qquad c_3=1,$$
por lo que
$$N_3=1^6\cdot 2^3\cdot 3=24=2^3\cdot 3.$$
Queremos los divisores cuya suma total de exponentes sea múltiplo de \(3\). Para el primo \(2\) con exponente \(3\), las sumas por residuo son
$$C_{2,0}=1+2^3=9,\qquad C_{2,1}=2,\qquad C_{2,2}=4.$$
Para el primo \(3\) con exponente \(1\), obtenemos
$$C_{3,0}=1,\qquad C_{3,1}=3,\qquad C_{3,2}=0.$$
Para que el residuo total sea \(0\pmod 3\), se combinan los pares \((0,0)\), \((2,1)\) y \((1,2)\). Por tanto la suma pedida es
$$9\cdot 1+4\cdot 3+2\cdot 0=21.$$
En efecto, los divisores válidos son \(1\), \(8\) y \(12\), cuya suma es \(21\). El caso grande repite exactamente esta misma convolución, pero con muchas más factorizaciones primas.
Después de reducir los exponentes módulo \(m\), cada primo queda representado por un vector de longitud \(m\). Incorporar otro primo significa hacer una convolución cíclica sobre esos \(m\) residuos. No hace falta un espacio de estados mayor porque el problema no distingue entre sumas de exponentes que difieren en un múltiplo de \(m\).
Las implementaciones en C++, Python y Java construyen primero una tabla de menor factor primo hasta \(n\). Luego recorren \(t=1,2,\dots,n\), factorizan cada \(t\) y añaden \(c_t\,v_p(t)\) al exponente acumulado de cada primo \(p\) que aparece. Al terminar este paso, ya conocen todos los \(E_p\) no nulos de la factorización de \(N_n\).
Para cada primo \(p\), la implementación calcula los \(m\) valores \(C_{p,0},C_{p,1},\dots,C_{p,m-1}\) módulo \(M\). En lugar de sumar potencias una por una, mantiene el factor actual \(p^r\) y evalúa la serie geométrica restante en forma cerrada. El caso especial \(p^m\equiv 1\pmod M\) se trata por separado para evitar invertir un denominador nulo.
Después de procesar cualquier subconjunto de los primos, el vector de estado cumple la invariante
El estado \(\mathrm{dp}[r]\) almacena la suma de los valores de divisores formados con los primos ya procesados cuyo residuo total de exponentes es \(r\).
Al incorporar un nuevo primo, la transición es la convolución cíclica
$$\mathrm{next}[(r+t)\bmod m]\;{+}{=}\;\mathrm{dp}[r]\cdot C_{p,t}\pmod M.$$
Al final, \(\mathrm{dp}[0]\) es exactamente la respuesta. Las versiones en C++ y Java también comparan la rutina modular con un cálculo exacto muy pequeño, lo que sirve como comprobación independiente de la derivación.
La criba usada para construir la tabla de menor factor primo es lineal en \(n\). Factorizar todos los \(t\le n\) mediante divisiones repetidas sobre esa tabla es casi lineal para estos tamaños de entrada y bastante más barato que la fase final de convoluciones.
El coste dominante es la programación dinámica por residuos. Si \(\pi(n)\) es el número de primos menores o iguales que \(n\), la transición total requiere \(O(\pi(n)m^2)\) operaciones modulares, porque cada primo aporta una convolución cíclica de tamaño \(m\times m\). El uso de memoria es \(O(n+m)\): la criba y los exponentes ocupan \(O(n)\), y solo se necesitan unos pocos vectores de DP de longitud \(m\).
第 956 题的核心对象可以压缩成一个单独的整数：
$$N_n=\prod_{t=1}^{n} t^{c_t},\qquad c_t=\binom{n-t+2}{2}=\frac{(n-t+1)(n-t+2)}{2}.$$
要求的是所有满足 \(d\mid N_n\) 且
$$\Omega(d)=\sum_p \alpha_p$$
能被 \(m\) 整除的除数 \(d=\prod_p p^{\alpha_p}\) 的数值之和。实现计算的是 \(n=m=1000\) 的情形，并对 \(999999001\) 取模。
如果从“把所有除数列出来再筛选”的角度出发，这个问题完全无法做。即使先把 \(N_{1000}\) 做素因子分解，各个素数的指数 \(E_p\) 也已经非常巨大，因此除数总数天文级别。真正可行的办法，是把每个素数的贡献分开处理，并且只跟踪总指数和在模 \(m\) 意义下的余数。
整个解法本质上是“按素数逐个处理的生成函数计算”。因为条件 \(\Omega(d)\equiv 0 \pmod m\) 只依赖指数的余数，而除数的值 \(d\) 本身又对不同素数完全乘法分离，这两个结构正好可以组合在一起。
代码并不会真的去展开那个层层嵌套的乘积，而是直接统计每个整数 \(t\) 在完全展开后会出现多少次。这个重数正是三角数 \(c_t\)，所以整体对象恰好可以写成 \(N_n=\prod t^{c_t}\)。
于是对任意素数 \(p\)，它在 \(N_n\) 中的指数就是
$$E_p=v_p(N_n)=\sum_{t=1}^{n} c_t\,v_p(t).$$
一旦所有 \(E_p\) 都算出来，原问题就已经被化成了基于
$$N_n=\prod_p p^{E_p}$$
这组素因子指数的除数求和问题。
\(N_n\) 的每个除数都可以写成
$$d=\prod_p p^{\alpha_p},\qquad 0\le \alpha_p\le E_p.$$
题目要求求和的不是“满足条件的除数个数”，而是这些除数的实际数值。所以一个指数向量 \((\alpha_p)\) 的贡献就是
$$\prod_p p^{\alpha_p}.$$
全局限制只有一个：
$$\Omega(d)=\sum_p \alpha_p\equiv 0 \pmod m.$$
这种“总和模 \(m\) 的限制”非常适合用按余数分类的卷积来处理。
对固定素数 \(p\)，定义
$$F_p(x)=\sum_{a=0}^{E_p} p^a x^a.$$
把所有素数的这类因子相乘：
$$F(x)=\prod_p F_p(x).$$
\(F(x)\) 中 \(x^k\) 的系数，恰好就是所有满足 \(\Omega(d)=k\) 的除数 \(d\mid N_n\) 的数值之和。原因很直接：对每个素数都选出一项 \(p^{\alpha_p}x^{\alpha_p}\)，乘起来以后，\(p^{\alpha_p}\) 部分给出除数的值，而 \(x\) 的指数把总的 \(\sum_p \alpha_p\) 记录下来。
但我们并不关心精确的 \(k\)，只关心 \(k\bmod m\)。因此可以在满足 \(x^m=1\) 的商环里工作。这样每个素数只需要保留 \(m\) 个余数类：
$$F_p(x)\equiv \sum_{r=0}^{m-1} C_{p,r}x^r \pmod{x^m-1},$$
其中
$$C_{p,r}=\sum_{\substack{0\le a\le E_p\\ a\equiv r \pmod m}} p^a.$$
最终答案就是这些约化后多项式乘积中的 \(x^0\) 系数。
如果 \(E_p<r\)，那么 \(C_{p,r}=0\)。否则允许的指数为
$$a=r,\ r+m,\ r+2m,\ \dots,\ r+(N_{p,r}-1)m,$$
其中
$$N_{p,r}=\left\lfloor\frac{E_p-r}{m}\right\rfloor+1.$$
令 \(q=p^m\)，则有
$$C_{p,r}=p^r\sum_{j=0}^{N_{p,r}-1} q^j.$$
在模 \(M=999999001\) 下，这就是一个等比数列。若 \(q\not\equiv 1 \pmod M\)，则
$$C_{p,r}\equiv p^r\,(q^{N_{p,r}}-1)\,(q-1)^{-1}\pmod M.$$
若 \(q\equiv 1 \pmod M\)，分母会变成零，此时应改写为
$$C_{p,r}\equiv p^r\,N_{p,r}\pmod M.$$
这就是整个算法里唯一需要特别小心的模运算分支。
此时
$$c_1=6,\qquad c_2=3,\qquad c_3=1,$$
所以
$$N_3=1^6\cdot 2^3\cdot 3=24=2^3\cdot 3.$$
我们需要的是总指数个数为 \(3\) 的倍数的除数。对素数 \(2\) 来说，指数上界是 \(3\)，于是余数类求和为
$$C_{2,0}=1+2^3=9,\qquad C_{2,1}=2,\qquad C_{2,2}=4.$$
对素数 \(3\) 来说，指数上界是 \(1\)，因此
$$C_{3,0}=1,\qquad C_{3,1}=3,\qquad C_{3,2}=0.$$
若总余数要等于 \(0\pmod 3\)，可用的组合只有 \((0,0)\)、\((2,1)\) 和 \((1,2)\)。因此所求和为
$$9\cdot 1+4\cdot 3+2\cdot 0=21.$$
这与直接枚举得到的结果一致：合法除数是 \(1\)、\(8\) 和 \(12\)，总和正好为 \(21\)。大规模输入下的程序做的只是同一种卷积，只不过素数更多、指数更大而已。
一旦把指数按模 \(m\) 归类，每个素数就只剩下一个长度为 \(m\) 的系数向量。再乘进一个新的素数，本质上只是对这 \(m\) 个余数做一次循环卷积。由于题目根本不会区分相差 \(m\) 的倍数的指数总和，所以没有必要保存更大的状态空间。
C++、Python 和 Java 实现都会先建立一个到 \(n\) 为止的最小素因子表。接着遍历 \(t=1,2,\dots,n\)，对每个 \(t\) 做分解，并把 \(c_t\,v_p(t)\) 加到对应素数 \(p\) 的累计指数上。完成这一轮以后，\(N_n\) 的所有非零素数指数 \(E_p\) 就都得到了。
对每个素数 \(p\)，实现会计算 \(C_{p,0},C_{p,1},\dots,C_{p,m-1}\) 这 \(m\) 个值。它不会逐项累加幂，而是维护当前的 \(p^r\)，再用闭式一次性算出后面的等比和。若出现 \(p^m\equiv 1\pmod M\) 的情况，则单独处理，避免去求一个零分母的模逆。
处理完任意一个素数子集之后，状态向量满足：
状态 \(\mathrm{dp}[r]\) 表示：只使用已经处理过的素数所形成、且总指数余数为 \(r\) 的所有除数值之和。
加入一个新的素数时，转移就是循环卷积
$$\mathrm{next}[(r+t)\bmod m]\;{+}{=}\;\mathrm{dp}[r]\cdot C_{p,t}\pmod M.$$
当最后一个素数处理完成，\(\mathrm{dp}[0]\) 就是答案。C++ 和 Java 版本还会把模运算结果与一个很小规模的精确计算进行比较，这为整个推导提供了额外的校验。
建立最小素因子表所用的筛法对 \(n\) 是线性的。利用这张表对所有 \(t\le n\) 反复整除分解，在本题的输入范围内也非常快，而且明显比最后的卷积阶段便宜。
真正的主导成本来自余数动态规划。若 \(\pi(n)\) 表示不超过 \(n\) 的素数个数，那么整体转移需要 \(O(\pi(n)m^2)\) 次模运算，因为每个素数都要做一次 \(m\times m\) 的循环卷积。空间复杂度为 \(O(n+m)\)：筛表和指数数组占 \(O(n)\)，再加上常数个长度为 \(m\) 的 DP 缓冲区。
Центральный объект в задаче 956 можно свернуть в одно число:
$$N_n=\prod_{t=1}^{n} t^{c_t},\qquad c_t=\binom{n-t+2}{2}=\frac{(n-t+1)(n-t+2)}{2}.$$
Нужно просуммировать все делители \(d\mid N_n\), для которых, если записать \(d=\prod_p p^{\alpha_p}\), то
$$\Omega(d)=\sum_p \alpha_p$$
делится на \(m\). Реализации вычисляют случай \(n=m=1000\) и берут ответ по модулю \(999999001\).
Прямой перебор делителей здесь невозможен. Даже после разложения \(N_{1000}\) на простые множители показатели \(E_p\) получаются огромными, а значит, число делителей становится астрономическим. Рабочая идея состоит в том, чтобы отделить вклад каждого простого числа и хранить только остаток суммы показателей по модулю \(m\).
Вся конструкция решения сводится к производящей функции, обрабатываемой по простым числам. Условие \(\Omega(d)\equiv 0 \pmod m\) зависит только от остатков показателей, тогда как значение делителя \(d\) раскладывается мультипликативно по простым множителям.
Код не раскрывает вложенное произведение буквально. Вместо этого он считает, сколько раз каждое число \(t\) встречается после полного распрямления выражения. Эта кратность равна треугольному числу \(c_t\), поэтому весь объект имеет вид \(N_n=\prod t^{c_t}\).
Следовательно, для любого простого \(p\) его показатель в \(N_n\) равен
$$E_p=v_p(N_n)=\sum_{t=1}^{n} c_t\,v_p(t).$$
Как только найдены все \(E_p\), исходная задача превращается в задачу о сумме делителей для разложения
$$N_n=\prod_p p^{E_p}.$$
Любой делитель числа \(N_n\) записывается в форме
$$d=\prod_p p^{\alpha_p},\qquad 0\le \alpha_p\le E_p.$$
Суммируются не количества подходящих делителей, а сами значения делителей. Поэтому вклад выбора \((\alpha_p)\) равен
$$\prod_p p^{\alpha_p}.$$
Единственное глобальное ограничение имеет вид
$$\Omega(d)=\sum_p \alpha_p\equiv 0 \pmod m.$$
Именно такое ограничение естественно обрабатывается сверткой по классам вычетов.
Для фиксированного простого \(p\) введем
$$F_p(x)=\sum_{a=0}^{E_p} p^a x^a.$$
Тогда произведение по всем простым равно
$$F(x)=\prod_p F_p(x).$$
Коэффициент при \(x^k\) в \(F(x)\) есть в точности сумма значений всех делителей \(d\mid N_n\), для которых \(\Omega(d)=k\). Действительно, выбор члена \(p^{\alpha_p}x^{\alpha_p}\) одновременно добавляет значение делителя и записывает суммарное число показателей в степень \(x\).
Но точное значение \(k\) нам не нужно, важен только остаток \(k\bmod m\). Поэтому достаточно работать по модулю \(x^m-1\). В этом кольце каждый простой множитель оставляет только \(m\) классов:
$$F_p(x)\equiv \sum_{r=0}^{m-1} C_{p,r}x^r \pmod{x^m-1},$$
где
$$C_{p,r}=\sum_{\substack{0\le a\le E_p\\ a\equiv r \pmod m}} p^a.$$
Искомый ответ является коэффициентом при \(x^0\) в произведении этих сокращенных полиномов.
Если \(E_p<r\), то \(C_{p,r}=0\). Иначе допустимые показатели имеют вид
$$a=r,\ r+m,\ r+2m,\ \dots,\ r+(N_{p,r}-1)m,$$
где
$$N_{p,r}=\left\lfloor\frac{E_p-r}{m}\right\rfloor+1.$$
Обозначив \(q=p^m\), получаем
$$C_{p,r}=p^r\sum_{j=0}^{N_{p,r}-1} q^j.$$
По модулю \(M=999999001\) это геометрическая прогрессия. Если \(q\not\equiv 1 \pmod M\), то
$$C_{p,r}\equiv p^r\,(q^{N_{p,r}}-1)\,(q-1)^{-1}\pmod M.$$
Если же \(q\equiv 1 \pmod M\), знаменатель обращается в нуль, и тогда сумма упрощается до
$$C_{p,r}\equiv p^r\,N_{p,r}\pmod M.$$
Это единственный нетривиальный модульный особый случай во всей задаче.
В этом случае
$$c_1=6,\qquad c_2=3,\qquad c_3=1,$$
поэтому
$$N_3=1^6\cdot 2^3\cdot 3=24=2^3\cdot 3.$$
Нужны делители, для которых суммарное число показателей кратно \(3\). Для простого \(2\) с показателем \(3\) имеем
$$C_{2,0}=1+2^3=9,\qquad C_{2,1}=2,\qquad C_{2,2}=4.$$
Для простого \(3\) с показателем \(1\) получаем
$$C_{3,0}=1,\qquad C_{3,1}=3,\qquad C_{3,2}=0.$$
Чтобы общий остаток был \(0\pmod 3\), подходят комбинации \((0,0)\), \((2,1)\) и \((1,2)\). Следовательно, искомая сумма равна
$$9\cdot 1+4\cdot 3+2\cdot 0=21.$$
И действительно, подходящие делители — это \(1\), \(8\) и \(12\), а их сумма равна \(21\). Полный алгоритм для \(n=m=1000\) делает ровно ту же свертку, только на значительно большем наборе простых чисел.
После сведения показателей по модулю \(m\) каждый простой множитель представлен вектором длины \(m\). Добавление следующего простого — это просто циклическая свертка по этим \(m\) остаткам. Более крупное пространство состояний не нужно, потому что задача не различает суммы показателей, отличающиеся на кратное \(m\).
Реализации на C++, Python и Java сначала строят таблицу наименьших простых делителей до \(n\). Затем они проходят по \(t=1,2,\dots,n\), раскладывают каждое \(t\) и добавляют \(c_t\,v_p(t)\) к накопленному показателю каждого встреченного простого \(p\). После этого прохода все ненулевые \(E_p\) уже известны.
Для каждого простого \(p\) программа вычисляет \(m\) чисел \(C_{p,0},C_{p,1},\dots,C_{p,m-1}\) по модулю \(M\). Вместо поэлементного суммирования степеней она поддерживает текущий множитель \(p^r\) и вычисляет оставшуюся геометрическую сумму по замкнутой формуле. Особый случай \(p^m\equiv 1\pmod M\) обрабатывается отдельно, чтобы не брать обратный элемент к нулю.
После обработки любого подмножества простых чисел вектор состояния удовлетворяет инварианту
Состояние \(\mathrm{dp}[r]\) хранит сумму значений делителей, составленных из уже обработанных простых чисел, у которых общий остаток суммы показателей равен \(r\).
При добавлении нового простого переход имеет вид
$$\mathrm{next}[(r+t)\bmod m]\;{+}{=}\;\mathrm{dp}[r]\cdot C_{p,t}\pmod M.$$
После обработки последнего простого \(\mathrm{dp}[0]\) и есть ответ. Версии на C++ и Java дополнительно сверяют модульную процедуру с маленьким точным вычислением, что служит хорошей внутренней проверкой формул.
Построение таблицы наименьших простых делителей выполняется за \(O(n)\) благодаря линейному решету. Разложение всех \(t\le n\) с помощью этой таблицы для данных размеров входа также работает очень быстро и заметно дешевле финальной фазы сверток.
Основная стоимость сосредоточена в DP по остаткам. Если \(\pi(n)\) — число простых, не превосходящих \(n\), то полный переход по всем простым требует \(O(\pi(n)m^2)\) модульных операций, потому что каждый простой вносит циклическую свертку размера \(m\times m\). Память равна \(O(n+m)\): решето и массив показателей занимают \(O(n)\), а дополнительно нужны лишь несколько DP-буферов длины \(m\).
يمكن اختزال الكمية المركزية في المسألة 956 إلى عدد واحد هو:
$$N_n=\prod_{t=1}^{n} t^{c_t},\qquad c_t=\binom{n-t+2}{2}=\frac{(n-t+1)(n-t+2)}{2}.$$
المطلوب هو جمع جميع القواسم \(d\mid N_n\) التي تحقق الشرط التالي عندما نكتب
$$d=\prod_p p^{\alpha_p}:$$
وهو أن
$$\Omega(d)=\sum_p \alpha_p$$
يكون من مضاعفات \(m\). التنفيذات تحسب الحالة \(n=m=1000\) ثم تأخذ النتيجة بترديد \(999999001\).
العد المباشر للقواسم غير ممكن عمليًا. حتى بعد تحليل \(N_{1000}\) إلى عوامله الأولية، تكون الأسس \(E_p\) كبيرة جدًا، وبالتالي يصبح عدد القواسم هائلًا. الفكرة الصحيحة هي فصل مساهمة كل عدد أولي على حدة، ثم تتبع باقي مجموع الأسس فقط modulo \(m\).
الحل كله مبني على دالة مولدة تُعالَج أوليًا بعد أولي. السبب هو أن الشرط \(\Omega(d)\equiv 0 \pmod m\) يعتمد فقط على بواقي الأسس، بينما قيمة القاسم \(d\) نفسها تتفكك ضربياً عبر الأعداد الأولية.
الكود لا يوسع الجداء المتداخل حرفيًا. بدلًا من ذلك، يحسب كم مرة يظهر كل عدد \(t\) بعد تسطيح التعبير بالكامل. هذه المضاعفية هي العدد المثلثي \(c_t\)، ولذلك يصبح الجسم كاملًا مساويًا لـ \(N_n=\prod t^{c_t}\).
ومن ثم، لأي عدد أولي \(p\)، يكون أسه في \(N_n\) مساويًا لـ
$$E_p=v_p(N_n)=\sum_{t=1}^{n} c_t\,v_p(t).$$
وبمجرد معرفة جميع القيم \(E_p\)، تتحول المسألة الأصلية إلى مسألة مجموع قواسم على التحليل
$$N_n=\prod_p p^{E_p}.$$
كل قاسم لـ \(N_n\) يمكن كتابته على الصورة
$$d=\prod_p p^{\alpha_p},\qquad 0\le \alpha_p\le E_p.$$
المجموع المطلوب ليس عدد القواسم المقبولة، بل مجموع قيمها الفعلية. لذلك فإن اختيار متجه الأسس \((\alpha_p)\) يضيف المقدار
$$\prod_p p^{\alpha_p}.$$
والقيد العالمي الوحيد هو
$$\Omega(d)=\sum_p \alpha_p\equiv 0 \pmod m.$$
وهذا بالضبط نوع القيود الذي يناسب الالتفاف على أصناف البواقي.
لكل أولي ثابت \(p\) نعرّف
$$F_p(x)=\sum_{a=0}^{E_p} p^a x^a.$$
وبضرب هذه العوامل على جميع الأعداد الأولية نحصل على
$$F(x)=\prod_p F_p(x).$$
معامل \(x^k\) في \(F(x)\) يساوي تمامًا مجموع قيم جميع القواسم \(d\mid N_n\) التي تحقق \(\Omega(d)=k\). فاختيار الحد \(p^{\alpha_p}x^{\alpha_p}\) من كل أولي يسجل في آن واحد قيمة القاسم ومجموع عدد الأسس في قوة \(x\).
لكننا لا نحتاج القيمة الدقيقة لـ \(k\)، بل نحتاج فقط \(k\bmod m\). لذلك نعمل في الحلقة التي تحقق \(x^m=1\). عندها يكفي لكل أولي أن نحتفظ بـ \(m\) أصناف بواقي فقط:
$$F_p(x)\equiv \sum_{r=0}^{m-1} C_{p,r}x^r \pmod{x^m-1},$$
حيث
$$C_{p,r}=\sum_{\substack{0\le a\le E_p\\ a\equiv r \pmod m}} p^a.$$
وعليه فإن الجواب المطلوب هو معامل \(x^0\) في حاصل ضرب هذه كثيرات الحدود المختزلة.
إذا كان \(E_p<r\)، فإن \(C_{p,r}=0\). وإلا فإن الأسس الممكنة هي
$$a=r,\ r+m,\ r+2m,\ \dots,\ r+(N_{p,r}-1)m,$$
حيث
$$N_{p,r}=\left\lfloor\frac{E_p-r}{m}\right\rfloor+1.$$
وبوضع \(q=p^m\) نحصل على
$$C_{p,r}=p^r\sum_{j=0}^{N_{p,r}-1} q^j.$$
وبترديد \(M=999999001\) يصبح ذلك متسلسلة هندسية. فإذا كان \(q\not\equiv 1 \pmod M\)، فإن
$$C_{p,r}\equiv p^r\,(q^{N_{p,r}}-1)\,(q-1)^{-1}\pmod M.$$
أما إذا كان \(q\equiv 1 \pmod M\)، فإن المقام يختفي، وتتبسط الصيغة إلى
$$C_{p,r}\equiv p^r\,N_{p,r}\pmod M.$$
وهذه هي الحالة المعيارية الحساسة الوحيدة في الخوارزمية كلها.
في هذه الحالة لدينا
$$c_1=6,\qquad c_2=3,\qquad c_3=1,$$
ومن ثم
$$N_3=1^6\cdot 2^3\cdot 3=24=2^3\cdot 3.$$
نريد القواسم التي يكون مجموع أسسها من مضاعفات \(3\). بالنسبة إلى الأولي \(2\) ذي الأس \(3\)، تكون مجاميع أصناف البواقي هي
$$C_{2,0}=1+2^3=9,\qquad C_{2,1}=2,\qquad C_{2,2}=4.$$
وبالنسبة إلى الأولي \(3\) ذي الأس \(1\)، نحصل على
$$C_{3,0}=1,\qquad C_{3,1}=3,\qquad C_{3,2}=0.$$
ولكي يكون الباقي الكلي مساويًا لـ \(0\pmod 3\)، فإن التركيبات الممكنة هي \((0,0)\) و\((2,1)\) و\((1,2)\). لذلك يكون المجموع المطلوب
$$9\cdot 1+4\cdot 3+2\cdot 0=21.$$
وفعلًا، القواسم المقبولة هي \(1\) و\(8\) و\(12\)، ومجموعها \(21\). وما يفعله البرنامج في الحالة الكبيرة هو نفس هذا الالتفاف تمامًا، لكن عبر عدد أكبر من العوامل الأولية.
بعد اختزال الأسس modulo \(m\)، يصبح لكل أولي متجه معاملات طوله \(m\). وعند إدخال أولي جديد لا نحتاج إلا إلى التفاف دوري على هذه البواقي \(m\). لا توجد حاجة إلى فضاء حالات أكبر، لأن المسألة أصلًا لا تميّز بين مجموعي أسس يختلفان بمضاعف من \(m\).
تنفيذات C++ وPython وJava تبدأ ببناء جدول لأصغر عامل أولي حتى \(n\). ثم تمر على \(t=1,2,\dots,n\)، وتفكك كل \(t\) إلى عوامله، وتضيف \(c_t\,v_p(t)\) إلى الأس المتراكم لكل أولي \(p\) يظهر في هذا العدد. في نهاية هذه المرحلة تكون جميع القيم غير الصفرية لـ \(E_p\) قد حُسِبت.
لكل أولي \(p\)، يحسب التنفيذ القيم \(C_{p,0},C_{p,1},\dots,C_{p,m-1}\) بترديد \(M\). وبدلًا من جمع القوى حدًا حدًا، يحافظ على العامل الحالي \(p^r\) ثم يقيم باقي المتسلسلة الهندسية بصيغة مغلقة. وتعالج حالة \(p^m\equiv 1\pmod M\) على نحو خاص حتى لا يتم أخذ معكوس معياري لمقام صفري.
بعد معالجة أي مجموعة جزئية من الأعداد الأولية، يحقق متجه الحالة الثابت التالي:
الحالة \(\mathrm{dp}[r]\) تخزن مجموع قيم القواسم المكوّنة من الأوليات التي عولجت بالفعل والتي يكون باقي مجموع أسسها الكلي مساويًا لـ \(r\).
وعند إضافة أولي جديد يصبح الانتقال
$$\mathrm{next}[(r+t)\bmod m]\;{+}{=}\;\mathrm{dp}[r]\cdot C_{p,t}\pmod M.$$
وبعد معالجة آخر أولي تكون \(\mathrm{dp}[0]\) هي الجواب تمامًا. كما أن نسختي C++ وJava تقارنان الروتين المعياري بحساب دقيق صغير الحجم، وهو فحص جيد لصحة الصيغ المستعملة.
المنخل المستخدم لبناء جدول أصغر عامل أولي خطي في \(n\). أما تفكيك جميع القيم \(t\le n\) باستخدام هذا الجدول فيبقى سريعًا جدًا لهذه الحدود، وهو أرخص بكثير من مرحلة الالتفاف النهائية.
الكلفة المسيطرة تأتي من البرمجة الديناميكية على البواقي. إذا رمزنا إلى عدد الأعداد الأولية حتى \(n\) بالرمز \(\pi(n)\)، فإن الانتقال الكلي عبر جميع الأوليات يحتاج \(O(\pi(n)m^2)\) عملية معيارية، لأن كل أولي يضيف التفافًا دوريًا بحجم \(m\times m\). استهلاك الذاكرة يساوي \(O(n+m)\): المنخل ومصفوفات الأسس يحتاجان \(O(n)\)، ويكفي فوق ذلك عدد ثابت من مخازن DP بطول \(m\).