We must evaluate
$$D(m,n)=\sum_{d \mid m}\sum_{k=1}^{n}\sigma_0(kd),$$
for \(m=200!\) and \(n=10^{12}\), modulo \(10^9+7\). Here \(\sigma_0=\tau\) is the divisor-counting function. A direct double loop is hopeless, so the solution rewrites the divisor sum into an inclusion-exclusion over the primes of \(200!\), plus fast queries to the summatory divisor function.
Write \(m=f!\) with \(f=200\), and rename the upper limit as \(N=n\). Let \(\mathbb{P}\) denote the set of prime numbers and define
$$\mathcal{P}=\{p\in\mathbb{P}: p\le f\}.$$
For each \(p\in\mathcal{P}\), the exponent of \(p\) in \(f!\) is
$$e_p=v_p(f!)=\sum_{j\ge 1}\left\lfloor\frac{f}{p^j}\right\rfloor.$$
The implementations exploit this prime-power structure very directly.
Every divisor of \(f!\) has the form
$$d=\prod_{p\in\mathcal{P}}p^{a_p},\qquad 0\le a_p\le e_p.$$
For a fixed integer \(k\), write \(b_p=v_p(k)\) for primes \(p\in\mathcal{P}\). Since \(\tau\) is multiplicative on prime powers, the sum over all divisors of \(f!\) separates into local contributions:
$$\sum_{d \mid f!}\tau(kd)=\left(\prod_{q\notin\mathcal{P}}(v_q(k)+1)\right)\prod_{p\in\mathcal{P}}\sum_{a=0}^{e_p}(a+b_p+1).$$
So the entire problem reduces to understanding the one-prime expression \(\sum_{a=0}^{e}(a+b+1)\).
Define the triangular-number function
$$T(t)=\frac{t(t+1)}{2}.$$
Then
$$\sum_{a=0}^{e}(a+b+1)=(e+1)(b+1)+\frac{e(e+1)}{2}=T(e+1)(b+1)-T(e)b.$$
This identity is the key algebraic step. The term \(T(e+1)(b+1)\) keeps the current exponent of \(p\) inside \(k\), while the correction term \(T(e)b\) corresponds to removing one copy of \(p\) when \(p\mid k\).
Now define
$$K=\prod_{p\in\mathcal{P}}T(e_p+1),\qquad \beta_p=\frac{T(e_p)}{T(e_p+1)}.$$
Choosing the correction term at a prime \(p\) introduces a factor \(-\beta_p\) and can only happen when \(p\) already divides \(k\). If \(Q\subseteq\mathcal{P}\) is the set of primes where we choose that correction, and
$$P_Q=\prod_{p\in Q}p,$$
then one copy of every prime in \(Q\) is removed from \(k\). Therefore, for each fixed \(k\),
$$\sum_{d \mid f!}\tau(kd)=K\sum_{Q\subseteq\mathcal{P},\ P_Q\mid k}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)\tau\left(\frac{k}{P_Q}\right).$$
This is exactly the inclusion-exclusion that the implementations realize through a recursive traversal of prime subsets.
Now sum over \(1\le k\le N\). For a fixed subset \(Q\), the condition \(P_Q\mid k\) lets us write \(k=P_Qr\), so
$$D(f!,N)=K\sum_{Q\subseteq\mathcal{P}}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)\sum_{r\le N/P_Q}\tau(r).$$
Introduce the summatory divisor function
$$S_\tau(x)=\sum_{r=1}^{x}\tau(r).$$
Then the whole problem becomes
$$D(f!,N)=K\sum_{Q\subseteq\mathcal{P}}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)S_\tau\left(\left\lfloor\frac{N}{P_Q}\right\rfloor\right).$$
If \(P_Q>N\), the corresponding term is zero, which explains the branch pruning in the recursive search.
The identity
$$S_\tau(x)=\sum_{d=1}^{x}\left\lfloor\frac{x}{d}\right\rfloor$$
counts factor pairs \((a,b)\) with \(ab\le x\). Splitting those lattice points across the hyperbola \(ab=x\) gives the standard Dirichlet-hyperbola formula
$$S_\tau(x)=2\sum_{d=1}^{\lfloor\sqrt{x}\rfloor}\left\lfloor\frac{x}{d}\right\rfloor-\lfloor\sqrt{x}\rfloor^2.$$
This is the large-\(x\) formula used by the implementation. Small values are handled by a precomputed prefix table.
Here \(\mathcal{P}=\{2,3\}\) and \(e_2=e_3=1\). Hence
$$T(2)=3,\qquad T(1)=1,\qquad K=3\cdot 3=9,\qquad \beta_2=\beta_3=\frac{1}{3}.$$
The subset formula becomes
$$D(6,5)=9S_\tau(5)-3S_\tau(2)-3S_\tau(1)+S_\tau(0).$$
Now
$$S_\tau(5)=1+2+2+3+2=10,\qquad S_\tau(2)=3,\qquad S_\tau(1)=1,\qquad S_\tau(0)=0,$$
so
$$D(6,5)=9\cdot 10-3\cdot 3-3\cdot 1=78.$$
A direct check agrees:
$$\sum_{k=1}^{5}\tau(k)=10,\qquad \sum_{k=1}^{5}\tau(2k)=17,\qquad \sum_{k=1}^{5}\tau(3k)=19,\qquad \sum_{k=1}^{5}\tau(6k)=32,$$
and therefore
$$10+17+19+32=78.$$
This small case shows exactly how the inclusion-exclusion reproduces the brute-force sum.
The C++, Python, and Java implementations first enumerate the primes up to \(200\) and compute each factorial exponent \(e_p\) with Legendre's formula. From those exponents they build the global factor \(K\) and the prime-specific ratios \(\beta_p\), always working modulo \(10^9+7\).
Next, they precompute \(S_\tau(x)\) for all \(x<10^6\) with a divisor sieve followed by prefix sums. For larger arguments they use the hyperbola formula above and store the results in a cache, so repeated requests for the same value are answered immediately.
The main recursive traversal walks through the primes in increasing order, maintains the current squarefree product \(P_Q\), the alternating sign, and the multiplicative weight \(K\prod_{p\in Q}\beta_p\). At each visited subset it adds the corresponding term
$$K\left(\prod_{p\in Q}\beta_p\right)S_\tau\left(\left\lfloor\frac{N}{P_Q}\right\rfloor\right)$$
with the appropriate sign. Because multiplying by another prime only increases \(P_Q\), the recursion stops as soon as the next product would exceed \(N\).
Let \(B=10^6\). Building the small prefix table costs \(O(B\log B)\) time and \(O(B)\) memory, because every divisor updates all of its multiples. The prime list and factorial exponents up to \(200\) are tiny by comparison.
Let \(R\) be the number of squarefree products of primes \(\le 200\) that do not exceed \(N\). The subset traversal uses \(O(R)\) arithmetic steps. If \(X\) is the set of distinct large arguments passed to \(S_\tau\), then the uncached large-\(x\) cost is
$$O\left(\sum_{x\in X}\sqrt{x}\right),$$
because each such value is computed once with the hyperbola identity and then memoized. In practice, the pruning \(P_Q\le N\) and the cache are what make the computation feasible.
Gesucht ist der Wert von
$$D(m,n)=\sum_{d \mid m}\sum_{k=1}^{n}\sigma_0(kd),$$
für \(m=200!\) und \(n=10^{12}\), modulo \(10^9+7\). Dabei ist \(\sigma_0=\tau\) die Teileranzahlfunktion. Eine direkte doppelte Iteration ist völlig unpraktikabel, daher formt die Lösung die Summe in eine Inklusions-Exklusions-Form über die Primzahlen von \(200!\) plus schnelle Abfragen der summatorischen Teilerfunktion um.
Schreibe \(m=f!\) mit \(f=200\) und nenne die obere Grenze \(N=n\). Sei \(\mathbb{P}\) die Menge aller Primzahlen und setze
$$\mathcal{P}=\{p\in\mathbb{P}: p\le f\}.$$
Für jedes \(p\in\mathcal{P}\) ist der Exponent von \(p\) in \(f!\)
$$e_p=v_p(f!)=\sum_{j\ge 1}\left\lfloor\frac{f}{p^j}\right\rfloor.$$
Genau diese Primzahlpotenzstruktur wird in den Implementierungen ausgenutzt.
Jeder Teiler von \(f!\) hat die Form
$$d=\prod_{p\in\mathcal{P}}p^{a_p},\qquad 0\le a_p\le e_p.$$
Für ein festes \(k\) schreiben wir \(b_p=v_p(k)\) für \(p\in\mathcal{P}\). Weil \(\tau\) auf Primzahlpotenzen multiplikativ ist, zerfällt die Summe über alle Teiler von \(f!\) in lokale Faktoren:
$$\sum_{d \mid f!}\tau(kd)=\left(\prod_{q\notin\mathcal{P}}(v_q(k)+1)\right)\prod_{p\in\mathcal{P}}\sum_{a=0}^{e_p}(a+b_p+1).$$
Damit reduziert sich das Problem auf die lokale Summe \(\sum_{a=0}^{e}(a+b+1)\).
Definiere
$$T(t)=\frac{t(t+1)}{2}.$$
Dann gilt
$$\sum_{a=0}^{e}(a+b+1)=(e+1)(b+1)+\frac{e(e+1)}{2}=T(e+1)(b+1)-T(e)b.$$
Das ist der entscheidende algebraische Schritt. Der Term \(T(e+1)(b+1)\) behält den vorhandenen \(p\)-Exponent in \(k\) bei, während der Korrekturterm \(T(e)b\) genau einem Entfernen einer Kopie von \(p\) entspricht, falls \(p\mid k\).
Setze
$$K=\prod_{p\in\mathcal{P}}T(e_p+1),\qquad \beta_p=\frac{T(e_p)}{T(e_p+1)}.$$
Wählt man an einer Primzahl \(p\) den Korrekturterm, so entsteht ein Faktor \(-\beta_p\), und dies ist nur möglich, wenn \(p\) bereits \(k\) teilt. Ist \(Q\subseteq\mathcal{P}\) die Menge aller so gewählten Primzahlen und
$$P_Q=\prod_{p\in Q}p,$$
dann wird aus \(k\) genau eine Kopie jeder Primzahl aus \(Q\) herausdividiert. Daher gilt für jedes feste \(k\)
$$\sum_{d \mid f!}\tau(kd)=K\sum_{Q\subseteq\mathcal{P},\ P_Q\mid k}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)\tau\left(\frac{k}{P_Q}\right).$$
Dies ist exakt die Inklusions-Exklusions-Formel, die die Implementierungen durch eine rekursive Traversierung der Primzahlmengen auswerten.
Nun summieren wir über \(1\le k\le N\). Für eine feste Menge \(Q\) erlaubt die Bedingung \(P_Q\mid k\) die Substitution \(k=P_Qr\). Damit folgt
$$D(f!,N)=K\sum_{Q\subseteq\mathcal{P}}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)\sum_{r\le N/P_Q}\tau(r).$$
Führe die summatorische Teilerfunktion ein:
$$S_\tau(x)=\sum_{r=1}^{x}\tau(r).$$
Dann wird das Problem zu
$$D(f!,N)=K\sum_{Q\subseteq\mathcal{P}}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)S_\tau\left(\left\lfloor\frac{N}{P_Q}\right\rfloor\right).$$
Falls \(P_Q>N\), verschwindet der Beitrag. Genau das erklärt das Abschneiden von Zweigen in der rekursiven Suche.
Aus dem Zählen von Faktorenpaaren folgt
$$S_\tau(x)=\sum_{d=1}^{x}\left\lfloor\frac{x}{d}\right\rfloor.$$
Durch Spiegelung der Gitterpunkte an der Hyperbel \(ab=x\) erhält man die Dirichlet-Hyperbel-Formel
$$S_\tau(x)=2\sum_{d=1}^{\lfloor\sqrt{x}\rfloor}\left\lfloor\frac{x}{d}\right\rfloor-\lfloor\sqrt{x}\rfloor^2.$$
Genau diese Formel verwenden die Implementierungen für große \(x\). Kleine Werte kommen aus einer vorberechneten Präfixtabelle.
Hier ist \(\mathcal{P}=\{2,3\}\) und \(e_2=e_3=1\). Also
$$T(2)=3,\qquad T(1)=1,\qquad K=3\cdot 3=9,\qquad \beta_2=\beta_3=\frac{1}{3}.$$
Die Teilmengenformel lautet dann
$$D(6,5)=9S_\tau(5)-3S_\tau(2)-3S_\tau(1)+S_\tau(0).$$
Mit
$$S_\tau(5)=10,\qquad S_\tau(2)=3,\qquad S_\tau(1)=1,\qquad S_\tau(0)=0$$
erhält man
$$D(6,5)=9\cdot 10-3\cdot 3-3\cdot 1=78.$$
Eine direkte Kontrolle bestätigt das:
$$\sum_{k=1}^{5}\tau(k)=10,\qquad \sum_{k=1}^{5}\tau(2k)=17,\qquad \sum_{k=1}^{5}\tau(3k)=19,\qquad \sum_{k=1}^{5}\tau(6k)=32,$$
also insgesamt
$$10+17+19+32=78.$$
Das Beispiel zeigt transparent, dass die Inklusions-Exklusions-Formel denselben Wert wie die direkte Auswertung liefert.
Die C++-, Python- und Java-Implementierungen erzeugen zuerst alle Primzahlen bis \(200\) und bestimmen mit der üblichen Fakultäts-Exponentenformel den Wert \(e_p\) für jede dieser Primzahlen. Daraus werden der globale Faktor \(K\) und die primspezifischen Verhältnisse \(\beta_p\) modulo \(10^9+7\) aufgebaut.
Anschließend wird \(S_\tau(x)\) für alle \(x<10^6\) über ein Teilersieb und anschließende Präfixsummen vorbereitet. Für größere Argumente wird die Hyperbel-Formel berechnet und das Ergebnis zwischengespeichert, damit identische Anfragen nicht erneut die Wurzel-Schleife durchlaufen.
Die Hauptrekursion geht die Primzahlen in aufsteigender Reihenfolge durch und hält das aktuelle quadratfreie Produkt \(P_Q\), das Vorzeichen und das Gewicht \(K\prod_{p\in Q}\beta_p\) fest. Für jede besuchte Teilmenge wird der passende Term addiert oder subtrahiert, und sobald das nächste Produkt größer als \(N\) würde, endet der betreffende Zweig sofort.
Sei \(B=10^6\). Das Aufbauen der kleinen Präfixtabelle benötigt \(O(B\log B)\) Zeit und \(O(B)\) Speicher, weil jeder Teiler alle seine Vielfachen aktualisiert. Die Primzahlliste und die Fakultäts-Exponenten bis \(200\) sind dagegen vernachlässigbar klein.
Sei \(R\) die Anzahl quadratfreier Produkte von Primzahlen \(\le 200\), die \(N\) nicht überschreiten. Die Teilmengentraversierung braucht \(O(R)\) arithmetische Schritte. Ist \(X\) die Menge der verschiedenen großen Argumente von \(S_\tau\), so kostet der nicht gecachte Anteil
$$O\left(\sum_{x\in X}\sqrt{x}\right),$$
weil jeder solche Wert genau einmal mit der Hyperbel-Identität berechnet und danach gespeichert wird. Praktisch entscheidend sind also das Abschneiden \(P_Q\le N\) und die Memoisierung.
Hesaplanması gereken ifade
$$D(m,n)=\sum_{d \mid m}\sum_{k=1}^{n}\sigma_0(kd),$$
ve burada \(m=200!\), \(n=10^{12}\), sonuç da \(10^9+7\) modunda istenir. \(\sigma_0=\tau\) pozitif bölen sayısı fonksiyonudur. Doğrudan tüm \(d\mid 200!\) ve tüm \(k\le n\) değerlerini gezmek mümkün olmadığından çözüm, toplamı \(200!\)'ın asal çarpanları üzerinde bir dahil et-çıkar toplamına ve hızlı bir bölen-toplamı ön-ek fonksiyonuna dönüştürür.
\(m=f!\) olacak şekilde \(f=200\) yazalım ve üst sınırı \(N=n\) ile gösterelim. \(\mathbb{P}\) asal sayılar kümesi olsun ve
$$\mathcal{P}=\{p\in\mathbb{P}: p\le f\}.$$
Her \(p\in\mathcal{P}\) için \(f!\) içindeki üs
$$e_p=v_p(f!)=\sum_{j\ge 1}\left\lfloor\frac{f}{p^j}\right\rfloor$$
olur. Uygulamanın kullandığı bütün yapı doğrudan bu asal kuvvet ayrışımından gelir.
\(f!\)'ın her böleni
$$d=\prod_{p\in\mathcal{P}}p^{a_p},\qquad 0\le a_p\le e_p$$
şeklindedir. Sabit bir \(k\) için \(p\in\mathcal{P}\) olduğunda \(b_p=v_p(k)\) diyelim. \(\tau\) asal kuvvetlerde çarpımsal davrandığı için, \(f!\)'ın tüm bölenleri üzerindeki toplam yerel çarpanlara ayrılır:
$$\sum_{d \mid f!}\tau(kd)=\left(\prod_{q\notin\mathcal{P}}(v_q(k)+1)\right)\prod_{p\in\mathcal{P}}\sum_{a=0}^{e_p}(a+b_p+1).$$
Dolayısıyla ana problem tek bir asal için çıkan \(\sum_{a=0}^{e}(a+b+1)\) ifadesini anlamaya indirgenir.
Üçgensel sayı fonksiyonunu
$$T(t)=\frac{t(t+1)}{2}$$
olarak tanımlayalım. O zaman
$$\sum_{a=0}^{e}(a+b+1)=(e+1)(b+1)+\frac{e(e+1)}{2}=T(e+1)(b+1)-T(e)b.$$
Asıl cebirsel kırılma noktası budur. \(T(e+1)(b+1)\) terimi \(k\) içindeki mevcut \(p\) üssünü korur; \(T(e)b\) düzeltmesi ise \(p\mid k\) olduğunda bir tane \(p\)'yi eksiltmeye karşılık gelir.
Şimdi
$$K=\prod_{p\in\mathcal{P}}T(e_p+1),\qquad \beta_p=\frac{T(e_p)}{T(e_p+1)}$$
tanımlarını yapalım. Bir asal \(p\) için düzeltme terimini seçmek \(-\beta_p\) katsayısı getirir ve bu ancak \(p\) zaten \(k\)'yi bölüyorsa mümkündür. Eğer \(Q\subseteq\mathcal{P}\) bu düzeltmenin seçildiği asallar kümesi ise ve
$$P_Q=\prod_{p\in Q}p,$$
o zaman \(Q\) içindeki her asal için \(k\)'den bir kopya çıkarılmış olur. Böylece sabit bir \(k\) için
$$\sum_{d \mid f!}\tau(kd)=K\sum_{Q\subseteq\mathcal{P},\ P_Q\mid k}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)\tau\left(\frac{k}{P_Q}\right)$$
elde edilir. Uygulamadaki özyinelemeli asal-altküme dolaşımı tam olarak bu dahil et-çıkar açılımını hesaplar.
Artık \(1\le k\le N\) üzerinde toplayabiliriz. Sabit bir \(Q\) için \(P_Q\mid k\) koşulu \(k=P_Qr\) yazmamızı sağlar. Böylece
$$D(f!,N)=K\sum_{Q\subseteq\mathcal{P}}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)\sum_{r\le N/P_Q}\tau(r)$$
olur. Şimdi bölen sayısı fonksiyonunun kümülatif toplamını
$$S_\tau(x)=\sum_{r=1}^{x}\tau(r)$$
diye tanımlarsak, sonuç
$$D(f!,N)=K\sum_{Q\subseteq\mathcal{P}}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)S_\tau\left(\left\lfloor\frac{N}{P_Q}\right\rfloor\right)$$
şeklini alır. \(P_Q>N\) olduğunda katkı sıfırdır; bu yüzden aramadaki dallar erkenden budanabilir.
Çarpan çiftlerini sayarak
$$S_\tau(x)=\sum_{d=1}^{x}\left\lfloor\frac{x}{d}\right\rfloor$$
eşitliğini elde ederiz. \(ab=x\) hiperbolünün iki tarafını eşleştirince standart Dirichlet-hiperbol özdeşliği gelir:
$$S_\tau(x)=2\sum_{d=1}^{\lfloor\sqrt{x}\rfloor}\left\lfloor\frac{x}{d}\right\rfloor-\lfloor\sqrt{x}\rfloor^2.$$
Büyük \(x\) için uygulamanın kullandığı formül budur. Küçük değerler ise önceden hazırlanmış bir ön-ek tablosundan okunur.
Bu durumda \(\mathcal{P}=\{2,3\}\) ve \(e_2=e_3=1\) olur. Dolayısıyla
$$T(2)=3,\qquad T(1)=1,\qquad K=3\cdot 3=9,\qquad \beta_2=\beta_3=\frac{1}{3}.$$
Alt küme formülü
$$D(6,5)=9S_\tau(5)-3S_\tau(2)-3S_\tau(1)+S_\tau(0)$$
haline gelir. Burada
$$S_\tau(5)=10,\qquad S_\tau(2)=3,\qquad S_\tau(1)=1,\qquad S_\tau(0)=0,$$
yani
$$D(6,5)=9\cdot 10-3\cdot 3-3\cdot 1=78.$$
Doğrudan kontrol de aynı sonucu verir:
$$\sum_{k=1}^{5}\tau(k)=10,\qquad \sum_{k=1}^{5}\tau(2k)=17,\qquad \sum_{k=1}^{5}\tau(3k)=19,\qquad \sum_{k=1}^{5}\tau(6k)=32,$$
ve toplam
$$10+17+19+32=78$$
olur. Böylece dahil et-çıkar formülünün küçük bir örnekte kaba kuvvet sonucu aynen verdiği görülür.
C++, Python ve Java uygulamaları önce \(200\)'e kadar tüm asalları çıkarır ve her asal için \(200!\) içindeki üs değerini Legendre formülüyle hesaplar. Ardından bu üslerden genel katsayı \(K\) ve asal başına oranlar \(\beta_p\) oluşturulur; bütün aritmetik \(10^9+7\) modunda yapılır.
Sonra \(x<10^6\) için \(S_\tau(x)\) değerleri bir bölen eleği ve ön-ek toplamlarıyla önceden hesaplanır. Daha büyük \(x\) değerleri için hiperbol formülü uygulanır ve sonuçlar önbelleğe alınır; böylece aynı sorgu tekrar geldiğinde karekök döngüsü yeniden çalışmaz.
Ana özyinelemeli dolaşım asalları artan sırayla gezer; mevcut karesiz çarpım \(P_Q\), işaret ve ağırlık \(K\prod_{p\in Q}\beta_p\) birlikte taşınır. Her ziyaret edilen altküme için uygun terim eklenir veya çıkarılır. Bir sonraki asal çarpımı \(N\)'yi aşacaksa dal hemen sonlandırılır.
\(B=10^6\) olsun. Küçük ön-ek tablosunu kurmak \(O(B\log B)\) zaman ve \(O(B)\) bellek ister; çünkü her bölen bütün katlarını günceller. \(200\)'e kadar asal listesi ve faktöriyel üsleri bunun yanında çok küçüktür.
\(R\), \(200\)'den küçük ya da eşit asalların \(N\)'yi aşmayan karesiz çarpımlarının sayısı olsun. Alt küme dolaşımı \(O(R)\) aritmetik adım kullanır. Eğer \(X\), \(S_\tau\) fonksiyonuna gönderilen farklı büyük argümanların kümesi ise, önbellekte olmayan büyük-\(x\) maliyeti
$$O\left(\sum_{x\in X}\sqrt{x}\right)$$
olur; çünkü her böyle değer hiperbol özdeşliğiyle bir kez hesaplanır ve sonra saklanır. Pratikte hesaplamayı mümkün kılan iki nokta, \(P_Q\le N\) budaması ve önbelleklemedir.
Hay que calcular
$$D(m,n)=\sum_{d \mid m}\sum_{k=1}^{n}\sigma_0(kd),$$
con \(m=200!\) y \(n=10^{12}\), módulo \(10^9+7\). Aquí \(\sigma_0=\tau\) es la función que cuenta divisores. Un recorrido directo sobre todos los divisores de \(200!\) y todos los \(k\le n\) es inviable, así que la solución transforma la suma doble en una inclusión-exclusión sobre los primos de \(200!\), apoyada por consultas rápidas a la función sumatoria de divisores.
Escribimos \(m=f!\) con \(f=200\) y denotamos por \(N=n\) el límite superior. Sea \(\mathbb{P}\) el conjunto de los números primos y definamos
$$\mathcal{P}=\{p\in\mathbb{P}: p\le f\}.$$
Para cada \(p\in\mathcal{P}\), el exponente de \(p\) en \(f!\) viene dado por
$$e_p=v_p(f!)=\sum_{j\ge 1}\left\lfloor\frac{f}{p^j}\right\rfloor.$$
Toda la derivación parte de esa descomposición en potencias primas.
Cada divisor de \(f!\) se escribe como
$$d=\prod_{p\in\mathcal{P}}p^{a_p},\qquad 0\le a_p\le e_p.$$
Fijado un entero \(k\), escribimos \(b_p=v_p(k)\) para \(p\in\mathcal{P}\). Como \(\tau\) es multiplicativa en potencias primas, la suma sobre todos los divisores se factoriza:
$$\sum_{d \mid f!}\tau(kd)=\left(\prod_{q\notin\mathcal{P}}(v_q(k)+1)\right)\prod_{p\in\mathcal{P}}\sum_{a=0}^{e_p}(a+b_p+1).$$
Por tanto, el problema entero se reduce a estudiar la suma local \(\sum_{a=0}^{e}(a+b+1)\).
Definimos
$$T(t)=\frac{t(t+1)}{2}.$$
Entonces
$$\sum_{a=0}^{e}(a+b+1)=(e+1)(b+1)+\frac{e(e+1)}{2}=T(e+1)(b+1)-T(e)b.$$
Éste es el paso algebraico central. El término \(T(e+1)(b+1)\) conserva el exponente actual de \(p\) dentro de \(k\), mientras que el término correctivo \(T(e)b\) equivale a quitar una copia de \(p\) cuando \(p\mid k\).
Introducimos
$$K=\prod_{p\in\mathcal{P}}T(e_p+1),\qquad \beta_p=\frac{T(e_p)}{T(e_p+1)}.$$
Escoger el término correctivo en un primo \(p\) aporta un factor \(-\beta_p\) y sólo es posible si \(p\) ya divide a \(k\). Si \(Q\subseteq\mathcal{P}\) es el conjunto de primos donde se hace esa elección y
$$P_Q=\prod_{p\in Q}p,$$
entonces se elimina una copia de cada primo de \(Q\) del valor \(k\). Así, para todo \(k\) fijo,
$$\sum_{d \mid f!}\tau(kd)=K\sum_{Q\subseteq\mathcal{P},\ P_Q\mid k}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)\tau\left(\frac{k}{P_Q}\right).$$
Ésa es exactamente la identidad que las implementaciones evalúan mediante un recorrido recursivo de subconjuntos de primos.
Ahora sumamos sobre \(1\le k\le N\). Para un subconjunto fijo \(Q\), la condición \(P_Q\mid k\) permite escribir \(k=P_Qr\). Entonces
$$D(f!,N)=K\sum_{Q\subseteq\mathcal{P}}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)\sum_{r\le N/P_Q}\tau(r).$$
Definimos la función sumatoria de divisores
$$S_\tau(x)=\sum_{r=1}^{x}\tau(r).$$
Con ello obtenemos
$$D(f!,N)=K\sum_{Q\subseteq\mathcal{P}}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)S_\tau\left(\left\lfloor\frac{N}{P_Q}\right\rfloor\right).$$
Si \(P_Q>N\), el término vale cero. Por eso la búsqueda puede podar inmediatamente esas ramas.
Contando pares de factores se obtiene
$$S_\tau(x)=\sum_{d=1}^{x}\left\lfloor\frac{x}{d}\right\rfloor.$$
Al emparejar los puntos de la retícula a ambos lados de la hipérbola \(ab=x\), se llega a la identidad clásica de Dirichlet:
$$S_\tau(x)=2\sum_{d=1}^{\lfloor\sqrt{x}\rfloor}\left\lfloor\frac{x}{d}\right\rfloor-\lfloor\sqrt{x}\rfloor^2.$$
Ésta es la fórmula usada para valores grandes de \(x\). Para valores pequeños, la implementación consulta una tabla prefijo precalculada.
Aquí \(\mathcal{P}=\{2,3\}\) y \(e_2=e_3=1\). Por tanto
$$T(2)=3,\qquad T(1)=1,\qquad K=3\cdot 3=9,\qquad \beta_2=\beta_3=\frac{1}{3}.$$
La fórmula por subconjuntos queda
$$D(6,5)=9S_\tau(5)-3S_\tau(2)-3S_\tau(1)+S_\tau(0).$$
Como
$$S_\tau(5)=10,\qquad S_\tau(2)=3,\qquad S_\tau(1)=1,\qquad S_\tau(0)=0,$$
obtenemos
$$D(6,5)=9\cdot 10-3\cdot 3-3\cdot 1=78.$$
La comprobación directa coincide:
$$\sum_{k=1}^{5}\tau(k)=10,\qquad \sum_{k=1}^{5}\tau(2k)=17,\qquad \sum_{k=1}^{5}\tau(3k)=19,\qquad \sum_{k=1}^{5}\tau(6k)=32,$$
y por tanto
$$10+17+19+32=78.$$
Este ejemplo pequeño muestra con claridad que la fórmula transformada reproduce exactamente la suma ingenua.
Las implementaciones en C++, Python y Java generan primero todos los primos hasta \(200\) y calculan, para cada uno, su exponente en \(200!\) mediante la fórmula de Legendre. A partir de esos exponentes construyen el factor global \(K\) y las razones \(\beta_p\), siempre módulo \(10^9+7\).
Después precalculan \(S_\tau(x)\) para todos los \(x<10^6\) usando un tamiz de divisores seguido de sumas prefijo. Para argumentos mayores aplican la identidad hiperbólica y guardan los resultados en una caché, de modo que una consulta repetida no vuelve a ejecutar el bucle hasta \(\sqrt{x}\).
El recorrido principal avanza por los primos en orden creciente y mantiene el producto libre de cuadrados \(P_Q\), el signo alternante y el peso \(K\prod_{p\in Q}\beta_p\). En cada subconjunto visitado agrega o resta el término correspondiente, y se detiene en cuanto el siguiente producto excedería \(N\).
Sea \(B=10^6\). Construir la tabla prefijo pequeña cuesta \(O(B\log B)\) tiempo y \(O(B)\) memoria, porque cada divisor actualiza todos sus múltiplos. La preparación de primos y exponentes hasta \(200\) es muy pequeña en comparación.
Sea \(R\) el número de productos libres de cuadrados de primos \(\le 200\) que no superan \(N\). El recorrido de subconjuntos usa \(O(R)\) pasos aritméticos. Si \(X\) es el conjunto de argumentos grandes distintos enviados a \(S_\tau\), el coste no cacheado para esos valores es
$$O\left(\sum_{x\in X}\sqrt{x}\right),$$
porque cada uno se calcula una sola vez con la identidad hiperbólica y luego queda memoizado. En la práctica, la poda \(P_Q\le N\) y la caché son las dos razones por las que el método funciona.
题目要求计算
$$D(m,n)=\sum_{d \mid m}\sum_{k=1}^{n}\sigma_0(kd),$$
其中 \(m=200!\)，\(n=10^{12}\)，结果对 \(10^9+7\) 取模。这里 \(\sigma_0=\tau\) 表示约数个数函数。若直接枚举 \(200!\) 的所有约数，再对每个 \(k\le n\) 逐项计算，规模完全不可行，所以解法必须先把双重求和改写成更有结构的形式。
实际做法是把 \(200!\) 的所有素因子贡献拆开，先得到一个按素数子集进行的容斥公式，再把内部求和压缩成“约数个数前缀和”函数的查询。这样，问题的核心就从暴力枚举转成了两个部分：一是如何构造正确的容斥权重，二是如何快速计算前缀和。
设 \(m=f!\)，这里 \(f=200\)，并把上界记为 \(N=n\)。记 \(\mathbb{P}\) 为全体素数的集合，并定义
$$\mathcal{P}=\{p\in\mathbb{P}: p\le f\}.$$
对于每个 \(p\in\mathcal{P}\)，它在 \(f!\) 中的指数是
$$e_p=v_p(f!)=\sum_{j\ge 1}\left\lfloor\frac{f}{p^j}\right\rfloor.$$
这正是 Legendre 公式。后面的所有推导都建立在这个素数幂分解上。
\(f!\) 的任意一个约数都可以写成
$$d=\prod_{p\in\mathcal{P}}p^{a_p},\qquad 0\le a_p\le e_p.$$
固定一个整数 \(k\)，对每个 \(p\in\mathcal{P}\) 记 \(b_p=v_p(k)\)。由于 \(\tau\) 在素数幂上是乘法型展开的，所以对所有 \(d\mid f!\) 求和时，可以把来自不同素数的影响分离出来：
$$\sum_{d \mid f!}\tau(kd)=\left(\prod_{q\notin\mathcal{P}}(v_q(k)+1)\right)\prod_{p\in\mathcal{P}}\sum_{a=0}^{e_p}(a+b_p+1).$$
这说明真正需要处理的只是单个素数上的局部和 \(\sum_{a=0}^{e}(a+b+1)\)。一旦这个局部表达式被整理好，整个公式就会自然地变成乘积结构。
定义三角数函数
$$T(t)=\frac{t(t+1)}{2}.$$
那么有
$$\sum_{a=0}^{e}(a+b+1)=(e+1)(b+1)+\frac{e(e+1)}{2}=T(e+1)(b+1)-T(e)b.$$
这是整个推导里最关键的代数整理。前一部分 \(T(e+1)(b+1)\) 可以理解为“保持 \(k\) 中当前的 \(p\) 次幂不变”，后一部分 \(T(e)b\) 则对应“如果 \(p\) 已经整除 \(k\)，就从 \(k\) 里去掉一个 \(p\)”。正是这个“保留”与“去掉一份”的二项结构，使容斥展开成为可能。
定义
$$K=\prod_{p\in\mathcal{P}}T(e_p+1),\qquad \beta_p=\frac{T(e_p)}{T(e_p+1)}.$$
当我们在某个素数 \(p\) 上选取修正项时，会额外得到一个因子 \(-\beta_p\)，而且这一项只有在 \(p\mid k\) 时才有意义，因为它对应于从 \(k\) 中拿掉一个 \(p\)。如果把所有选取修正项的素数组成集合 \(Q\subseteq\mathcal{P}\)，并记
$$P_Q=\prod_{p\in Q}p,$$
那么 \(Q\) 中每个素数都恰好从 \(k\) 里消去一次，于是得到
$$\sum_{d \mid f!}\tau(kd)=K\sum_{Q\subseteq\mathcal{P},\ P_Q\mid k}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)\tau\left(\frac{k}{P_Q}\right).$$
这一步就是代码所实现的核心恒等式。程序递归枚举素数子集，本质上就是在对这个容斥式逐项求值。
现在对 \(1\le k\le N\) 再求和。对固定的子集 \(Q\)，因为必须满足 \(P_Q\mid k\)，所以可以把 \(k\) 写成 \(k=P_Qr\)。于是
$$D(f!,N)=K\sum_{Q\subseteq\mathcal{P}}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)\sum_{r\le N/P_Q}\tau(r).$$
定义约数个数函数的前缀和
$$S_\tau(x)=\sum_{r=1}^{x}\tau(r).$$
那么总公式就变成
$$D(f!,N)=K\sum_{Q\subseteq\mathcal{P}}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)S_\tau\left(\left\lfloor\frac{N}{P_Q}\right\rfloor\right).$$
如果某个子集的乘积 \(P_Q\) 已经大于 \(N\)，那么这一项显然为零，因此递归时可以立刻剪枝。这正是实现中性能优化最直接的一点。
由“计数满足 \(ab\le x\) 的因子对”这一观点，可以得到
$$S_\tau(x)=\sum_{d=1}^{x}\left\lfloor\frac{x}{d}\right\rfloor.$$
接着把格点按双曲线 \(ab=x\) 对称配对，就得到经典的 Dirichlet 双曲线恒等式：
$$S_\tau(x)=2\sum_{d=1}^{\lfloor\sqrt{x}\rfloor}\left\lfloor\frac{x}{d}\right\rfloor-\lfloor\sqrt{x}\rfloor^2.$$
这正是实现处理大参数时所使用的公式。对于较小的 \(x\)，实现会先用筛法预处理出全部前缀值，随后直接查表；对于较大的 \(x\)，则用上式计算并把结果缓存起来，避免重复求值。
此时 \(\mathcal{P}=\{2,3\}\)，并且 \(e_2=e_3=1\)。所以
$$T(2)=3,\qquad T(1)=1,\qquad K=3\cdot 3=9,\qquad \beta_2=\beta_3=\frac{1}{3}.$$
总公式退化成
$$D(6,5)=9S_\tau(5)-3S_\tau(2)-3S_\tau(1)+S_\tau(0).$$
又因为
$$S_\tau(5)=10,\qquad S_\tau(2)=3,\qquad S_\tau(1)=1,\qquad S_\tau(0)=0,$$
所以
$$D(6,5)=9\cdot 10-3\cdot 3-3\cdot 1=78.$$
如果直接枚举 \(6\) 的约数 \(1,2,3,6\)，也会得到同样的值：
$$\sum_{k=1}^{5}\tau(k)=10,\qquad \sum_{k=1}^{5}\tau(2k)=17,\qquad \sum_{k=1}^{5}\tau(3k)=19,\qquad \sum_{k=1}^{5}\tau(6k)=32,$$
因此
$$10+17+19+32=78.$$
这个小例子非常适合说明公式为什么成立，也能直观看出容斥写法并不是近似，而是与暴力求和完全等价的重组。
C++、Python 和 Java 实现都先枚举不超过 \(200\) 的素数，再用 Legendre 公式计算每个素数在 \(200!\) 中的指数。基于这些指数，程序构造出整体乘子 \(K\) 以及每个素数对应的比例 \(\beta_p\)，整个过程中始终在 \(10^9+7\) 模下运算。
接下来，程序先对所有 \(x<10^6\) 预处理 \(S_\tau(x)\)。做法是先用一个约数筛统计每个整数的约数个数，再做前缀和。对于更大的参数，则使用上面的双曲线公式计算，并把结果放入缓存，因此相同的查询不会重复进入平方根级别的循环。
最后，主递归按素数递增顺序遍历所有可能的子集，维护当前的平方自由乘积 \(P_Q\)、当前符号以及权重 \(K\prod_{p\in Q}\beta_p\)。每访问到一个子集，就把对应项加到答案里或从答案里减掉；一旦继续乘上下一个素数会使乘积超过 \(N\)，这一支就直接停止。
设 \(B=10^6\)。构建小范围前缀表的时间复杂度是 \(O(B\log B)\)，空间复杂度是 \(O(B)\)，因为每个约数都会更新它的所有倍数。相比之下，求不超过 \(200\) 的素数以及对应指数的代价很小。
设 \(R\) 表示所有不超过 \(N\) 的、由素数 \(\le 200\) 组成的平方自由乘积数量。子集递归本身需要 \(O(R)\) 次基本运算。再设 \(X\) 为所有传给 \(S_\tau\) 的不同大参数集合，那么这些未命中缓存的查询总代价为
$$O\left(\sum_{x\in X}\sqrt{x}\right),$$
因为每个这样的值只会用双曲线恒等式计算一次，之后就直接复用缓存结果。实际运行中，真正让程序可行的是两点：一是 \(P_Q\le N\) 的剪枝，二是对大参数前缀值的记忆化。
Требуется вычислить
$$D(m,n)=\sum_{d \mid m}\sum_{k=1}^{n}\sigma_0(kd),$$
где \(m=200!\), \(n=10^{12}\), а ответ нужен по модулю \(10^9+7\). Здесь \(\sigma_0=\tau\) обозначает функцию числа делителей. Перебор всех \(d\mid 200!\) и всех \(k\le n\) невозможен, поэтому решение переводит исходную двойную сумму в гораздо более компактную формулу включений-исключений по простым числам, входящим в \(200!\).
Положим \(m=f!\), где \(f=200\), а верхнюю границу обозначим \(N=n\). Пусть \(\mathbb{P}\) обозначает множество всех простых чисел, и положим
$$\mathcal{P}=\{p\in\mathbb{P}: p\le f\}.$$
Для каждого \(p\in\mathcal{P}\) показатель этого простого в \(f!\) равен
$$e_p=v_p(f!)=\sum_{j\ge 1}\left\lfloor\frac{f}{p^j}\right\rfloor.$$
Это стандартная формула Лежандра, и именно она задает всю локальную структуру задачи.
Любой делитель числа \(f!\) имеет вид
$$d=\prod_{p\in\mathcal{P}}p^{a_p},\qquad 0\le a_p\le e_p.$$
Для фиксированного \(k\) обозначим через \(b_p=v_p(k)\) показатели при простых \(p\in\mathcal{P}\). Так как \(\tau\) раскладывается по простым степеням, сумма по всем делителям \(f!\) факторизуется:
$$\sum_{d \mid f!}\tau(kd)=\left(\prod_{q\notin\mathcal{P}}(v_q(k)+1)\right)\prod_{p\in\mathcal{P}}\sum_{a=0}^{e_p}(a+b_p+1).$$
Значит, достаточно понять локальную сумму \(\sum_{a=0}^{e}(a+b+1)\) для одного простого числа.
Введем функцию треугольных чисел
$$T(t)=\frac{t(t+1)}{2}.$$
Тогда
$$\sum_{a=0}^{e}(a+b+1)=(e+1)(b+1)+\frac{e(e+1)}{2}=T(e+1)(b+1)-T(e)b.$$
Это главный алгебраический шаг. Слагаемое \(T(e+1)(b+1)\) сохраняет текущую степень \(p\) внутри \(k\), а поправка \(T(e)b\) соответствует снятию одной копии \(p\), если \(p\) уже делит \(k\).
Определим
$$K=\prod_{p\in\mathcal{P}}T(e_p+1),\qquad \beta_p=\frac{T(e_p)}{T(e_p+1)}.$$
Если для простого \(p\) берется поправочный член, появляется множитель \(-\beta_p\), и это возможно только тогда, когда \(p\mid k\). Пусть \(Q\subseteq\mathcal{P}\) состоит из всех простых, на которых выбрана поправка, и
$$P_Q=\prod_{p\in Q}p.$$
Тогда из \(k\) убирается по одному множителю каждого простого из \(Q\), и для любого фиксированного \(k\) получается формула
$$\sum_{d \mid f!}\tau(kd)=K\sum_{Q\subseteq\mathcal{P},\ P_Q\mid k}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)\tau\left(\frac{k}{P_Q}\right).$$
Именно эту формулу реализации и вычисляют, перебирая подмножества простых рекурсивно.
Теперь просуммируем по всем \(1\le k\le N\). Для фиксированного подмножества \(Q\) условие \(P_Q\mid k\) позволяет написать \(k=P_Qr\). Поэтому
$$D(f!,N)=K\sum_{Q\subseteq\mathcal{P}}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)\sum_{r\le N/P_Q}\tau(r).$$
Введем сумматор функции числа делителей:
$$S_\tau(x)=\sum_{r=1}^{x}\tau(r).$$
Тогда окончательная форма имеет вид
$$D(f!,N)=K\sum_{Q\subseteq\mathcal{P}}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)S_\tau\left(\left\lfloor\frac{N}{P_Q}\right\rfloor\right).$$
Если \(P_Q>N\), соответствующий вклад равен нулю. Именно поэтому рекурсия может немедленно обрывать такие ветви.
Подсчет пар множителей дает тождество
$$S_\tau(x)=\sum_{d=1}^{x}\left\lfloor\frac{x}{d}\right\rfloor.$$
Если сгруппировать точки решетки по разные стороны гиперболы \(ab=x\), получаем стандартную формулу Дирихле:
$$S_\tau(x)=2\sum_{d=1}^{\lfloor\sqrt{x}\rfloor}\left\lfloor\frac{x}{d}\right\rfloor-\lfloor\sqrt{x}\rfloor^2.$$
Именно так реализации считают большие значения. Для малых аргументов используется заранее построенная префиксная таблица.
Здесь \(\mathcal{P}=\{2,3\}\) и \(e_2=e_3=1\). Следовательно,
$$T(2)=3,\qquad T(1)=1,\qquad K=3\cdot 3=9,\qquad \beta_2=\beta_3=\frac{1}{3}.$$
Формула по подмножествам принимает вид
$$D(6,5)=9S_\tau(5)-3S_\tau(2)-3S_\tau(1)+S_\tau(0).$$
Так как
$$S_\tau(5)=10,\qquad S_\tau(2)=3,\qquad S_\tau(1)=1,\qquad S_\tau(0)=0,$$
получаем
$$D(6,5)=9\cdot 10-3\cdot 3-3\cdot 1=78.$$
Прямая проверка совпадает:
$$\sum_{k=1}^{5}\tau(k)=10,\qquad \sum_{k=1}^{5}\tau(2k)=17,\qquad \sum_{k=1}^{5}\tau(3k)=19,\qquad \sum_{k=1}^{5}\tau(6k)=32,$$
поэтому
$$10+17+19+32=78.$$
Этот небольшой пример наглядно показывает, что формула включений-исключений не приближает ответ, а в точности переупорядочивает исходную сумму.
Реализации на C++, Python и Java сначала строят список всех простых чисел не больше \(200\) и для каждого такого простого вычисляют показатель в \(200!\) по формуле Лежандра. Затем из этих показателей составляются общий множитель \(K\) и отношения \(\beta_p\), причем все операции выполняются по модулю \(10^9+7\).
После этого заранее вычисляются значения \(S_\tau(x)\) для всех \(x<10^6\): сначала решетом считается число делителей каждого значения, затем строятся префиксные суммы. Для больших аргументов применяется гиперболическая формула, а результаты складываются в кэш, чтобы одинаковые запросы не пересчитывались заново.
Основной рекурсивный обход идет по простым в возрастающем порядке и хранит текущий квадратсвободный продукт \(P_Q\), знак и вес \(K\prod_{p\in Q}\beta_p\). Для каждого достигнутого подмножества к ответу добавляется или из него вычитается соответствующий вклад. Если умножение на следующий простой уже превышает \(N\), ветвь сразу прекращается.
Пусть \(B=10^6\). Построение малой префиксной таблицы занимает \(O(B\log B)\) времени и \(O(B)\) памяти, потому что каждый делитель обновляет все свои кратные. Подготовка простых и показателей до \(200\) на этом фоне мала.
Обозначим через \(R\) число квадратсвободных произведений простых \(\le 200\), не превосходящих \(N\). Тогда перебор подмножеств требует \(O(R)\) арифметических шагов. Если \(X\) - множество различных больших аргументов, переданных в \(S_\tau\), то стоимость некэшированных вычислений равна
$$O\left(\sum_{x\in X}\sqrt{x}\right),$$
поскольку каждое такое значение считается по гиперболической формуле только один раз, а затем переиспользуется из памяти. На практике решающими являются отсечение по условию \(P_Q\le N\) и мемоизация.
نريد حساب
$$D(m,n)=\sum_{d \mid m}\sum_{k=1}^{n}\sigma_0(kd),$$
عندما يكون \(m=200!\) و\(n=10^{12}\)، مع أخذ الناتج بترديد \(10^9+7\). هنا \(\sigma_0=\tau\) هي دالة عدد القواسم. الحساب المباشر عبر جميع القواسم \(d\mid 200!\) وجميع القيم \(k\le n\) غير عملي تمامًا، لذلك تُعاد صياغة المجموع المزدوج إلى صيغة شمول واستبعاد على أوليات \(200!\)، مع طريقة سريعة لحساب مجموع دالة عدد القواسم حتى حد معين.
نكتب \(m=f!\) حيث \(f=200\)، ونرمز للحد الأعلى بـ \(N=n\). ولنرمز بمجموعة الأعداد الأولية إلى \(\mathbb{P}\)، ثم نعرّف
$$\mathcal{P}=\{p\in\mathbb{P}: p\le f\}.$$
ولكل \(p\in\mathcal{P}\) فإن أس \(p\) في \(f!\) هو
$$e_p=v_p(f!)=\sum_{j\ge 1}\left\lfloor\frac{f}{p^j}\right\rfloor.$$
هذه هي صيغة Legendre، ومنها تنطلق البنية الحسابية كلها.
كل قاسم من قواسم \(f!\) يمكن كتابته على الصورة
$$d=\prod_{p\in\mathcal{P}}p^{a_p},\qquad 0\le a_p\le e_p.$$
ولعدد ثابت \(k\)، نكتب \(b_p=v_p(k)\) لكل \(p\in\mathcal{P}\). وبما أن \(\tau\) تتفكك على قوى الأوليات، فإن الجمع على جميع القواسم ينفصل إلى عوامل محلية:
$$\sum_{d \mid f!}\tau(kd)=\left(\prod_{q\notin\mathcal{P}}(v_q(k)+1)\right)\prod_{p\in\mathcal{P}}\sum_{a=0}^{e_p}(a+b_p+1).$$
إذًا جوهر المسألة هو فهم التعبير المحلي \(\sum_{a=0}^{e}(a+b+1)\) عند أولي واحد فقط.
نعرّف دالة الأعداد المثلثية:
$$T(t)=\frac{t(t+1)}{2}.$$
وعندئذ
$$\sum_{a=0}^{e}(a+b+1)=(e+1)(b+1)+\frac{e(e+1)}{2}=T(e+1)(b+1)-T(e)b.$$
هذه هي الخطوة الجبرية الأساسية. الحد \(T(e+1)(b+1)\) يبقي أس \(p\) الموجود داخل \(k\) كما هو، أما حد التصحيح \(T(e)b\) فيمثل نزع نسخة واحدة من \(p\) إذا كان \(p\) يقسم \(k\).
نعرّف
$$K=\prod_{p\in\mathcal{P}}T(e_p+1),\qquad \beta_p=\frac{T(e_p)}{T(e_p+1)}.$$
اختيار حد التصحيح عند أولي \(p\) يضيف العامل \(-\beta_p\)، ولا يحدث ذلك إلا إذا كان \(p\mid k\). إذا كانت \(Q\subseteq\mathcal{P}\) هي مجموعة الأوليات التي اخترنا عندها هذا التصحيح، وكتبنا
$$P_Q=\prod_{p\in Q}p,$$
فهذا يعني أننا نحذف عاملًا واحدًا من كل أولي في \(Q\) من داخل \(k\). لذلك نحصل، لكل \(k\) ثابت، على
$$\sum_{d \mid f!}\tau(kd)=K\sum_{Q\subseteq\mathcal{P},\ P_Q\mid k}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)\tau\left(\frac{k}{P_Q}\right).$$
وهذه هي بالضبط صيغة شمول-استبعاد التي تنفذها التطبيقات عبر استعراض تفرعي لمجموعات الأوليات.
بعد ذلك نجمع على جميع \(1\le k\le N\). عند تثبيت مجموعة \(Q\)، فإن الشرط \(P_Q\mid k\) يسمح بكتابة \(k=P_Qr\). ومن ثم
$$D(f!,N)=K\sum_{Q\subseteq\mathcal{P}}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)\sum_{r\le N/P_Q}\tau(r).$$
نعرّف الآن المجموع التراكمي لدالة عدد القواسم:
$$S_\tau(x)=\sum_{r=1}^{x}\tau(r).$$
فتصبح الصيغة النهائية
$$D(f!,N)=K\sum_{Q\subseteq\mathcal{P}}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)S_\tau\left(\left\lfloor\frac{N}{P_Q}\right\rfloor\right).$$
إذا كان \(P_Q>N\) فإن الحد المقابل يساوي صفرًا، ولهذا تقطع الخوارزمية هذا الفرع فورًا.
من عدّ أزواج العوامل نحصل على الهوية
$$S_\tau(x)=\sum_{d=1}^{x}\left\lfloor\frac{x}{d}\right\rfloor.$$
وبتقسيم نقاط الشبكة حول القطع الزائد \(ab=x\) نحصل على صيغة Dirichlet المعروفة:
$$S_\tau(x)=2\sum_{d=1}^{\lfloor\sqrt{x}\rfloor}\left\lfloor\frac{x}{d}\right\rfloor-\lfloor\sqrt{x}\rfloor^2.$$
هذه هي الصيغة التي تُستخدم للقيم الكبيرة. أما القيم الصغيرة فتُجلب من جدول تراكمي محسوب مسبقًا.
هنا \(\mathcal{P}=\{2,3\}\) و\(e_2=e_3=1\)، ولذلك
$$T(2)=3,\qquad T(1)=1,\qquad K=3\cdot 3=9,\qquad \beta_2=\beta_3=\frac{1}{3}.$$
فتصبح صيغة المجموعات
$$D(6,5)=9S_\tau(5)-3S_\tau(2)-3S_\tau(1)+S_\tau(0).$$
ومع
$$S_\tau(5)=10,\qquad S_\tau(2)=3,\qquad S_\tau(1)=1,\qquad S_\tau(0)=0,$$
نحصل على
$$D(6,5)=9\cdot 10-3\cdot 3-3\cdot 1=78.$$
والتحقق المباشر يعطي النتيجة نفسها:
$$\sum_{k=1}^{5}\tau(k)=10,\qquad \sum_{k=1}^{5}\tau(2k)=17,\qquad \sum_{k=1}^{5}\tau(3k)=19,\qquad \sum_{k=1}^{5}\tau(6k)=32,$$
ومن ثم
$$10+17+19+32=78.$$
هذا المثال الصغير يوضح أن الصيغة الجديدة ليست تقريبًا، بل إعادة تنظيم دقيقة للمجموع الأصلي.
تبدأ تطبيقات C++ وPython وJava بحصر جميع الأوليات حتى \(200\)، ثم تحسب أس كل أولي في \(200!\) باستخدام صيغة Legendre. بعد ذلك تُبنى الكمية العامة \(K\) والنسب \(\beta_p\) لكل أولي، وكل الحسابات تتم بترديد \(10^9+7\).
ثم تُجهَّز قيم \(S_\tau(x)\) لكل \(x<10^6\) بواسطة غربال للقواسم يتبعه جمع تراكمي. أما للقيم الأكبر، فتُستخدم صيغة القطع الزائد السابقة، وتُخزَّن النتائج في ذاكرة مؤقتة حتى لا يُعاد حساب الحلقة حتى \(\sqrt{x}\) أكثر من مرة لنفس القيمة.
أما الجزء الرئيس فيسير على الأوليات بترتيب تصاعدي، ويحمل حاصل الضرب الخالي من المربعات \(P_Q\)، والإشارة، والوزن \(K\prod_{p\in Q}\beta_p\). وعند كل مجموعة مزارة يضيف أو يطرح الحد المناسب. وإذا كان ضرب الأولي التالي سيجعل الناتج أكبر من \(N\)، فإن ذلك الفرع يتوقف فورًا.
لنفرض \(B=10^6\). بناء الجدول التراكمي الصغير يحتاج إلى زمن \(O(B\log B)\) وذاكرة \(O(B)\)، لأن كل قاسم يحدث جميع مضاعفاته. أما تجهيز قائمة الأوليات وأسسها حتى \(200\) فتكلفته صغيرة مقارنة بذلك.
إذا رمزنا بـ \(R\) إلى عدد الحواصل الخالية من المربعات المكوّنة من أوليات \(\le 200\) والتي لا تتجاوز \(N\)، فإن استعراض المجموعات يحتاج إلى \(O(R)\) خطوة حسابية. وإذا كانت \(X\) هي مجموعة القيم الكبيرة المختلفة التي تُمرَّر إلى \(S_\tau\)، فإن تكلفة الحسابات غير الموجودة في الذاكرة المؤقتة هي
$$O\left(\sum_{x\in X}\sqrt{x}\right),$$
لأن كل قيمة من هذه القيم تُحسب مرة واحدة فقط ثم تُخزن. عمليًا، العاملان الحاسمان هما القطع المبكر عند الشرط \(P_Q\le N\) والتخزين المؤقت لنتائج الدالة التراكمية.