A gozinta chain for \(n\) is a sequence
$$1=d_0 \lt d_1 \lt \cdots \lt d_k=n$$
in which every term divides the next one. Let \(g(n)\) be the number of such chains. The task is to sum all integers \(n\le 10^{16}\) satisfying
$$g(n)=n.$$
The decisive observation is that \(g(n)\) depends only on the sorted prime-exponent pattern of \(n\), not on the actual primes. That turns a hopeless search over all integers up to \(10^{16}\) into a manageable search over exponent patterns.
Write the prime factorization in non-increasing exponent order:
$$n=\prod_{i=1}^{r} p_i^{a_i},\qquad a_1\ge a_2\ge \cdots \ge a_r \gt 0.$$
Let \(\mathbf a=(a_1,\dots,a_r)\) and \(\omega=a_1+\cdots+a_r\).
Every divisor of \(n\) has the form
$$\prod_{i=1}^{r} p_i^{b_i},\qquad 0\le b_i\le a_i.$$
So a gozinta chain is equivalent to a coordinatewise non-decreasing sequence of exponent vectors starting at \((0,\dots,0)\) and ending at \((a_1,\dots,a_r)\).
Because that description uses only the exponents, the chain count is determined entirely by \(\mathbf a\). We may therefore write
$$g(n)=G(\mathbf a),$$
where \(G\) is a function on exponent patterns. This is the key reduction used by the implementation.
Suppose the chain has exactly \(k\) nontrivial divisibility steps:
$$1=d_0 \lt d_1 \lt \cdots \lt d_k=n.$$
Define the step ratios \(q_t=d_t/d_{t-1}\). Then each \(q_t \gt 1\) and
$$q_1q_2\cdots q_k=n.$$
For one prime \(p_i^{a_i}\), distribute the exponent \(a_i\) across those \(k\) ratios:
$$e_{i,1}+e_{i,2}+\cdots+e_{i,k}=a_i,\qquad e_{i,t}\ge 0.$$
The number of weak compositions of \(a_i\) into \(k\) parts is
$$\binom{a_i+k-1}{k-1}.$$
Independence across primes gives
$$T_k(\mathbf a)=\prod_{i=1}^{r}\binom{a_i+k-1}{k-1}.$$
This counts ordered \(k\)-step decompositions, but it still allows some ratios to equal \(1\), so it overcounts genuine chains.
A genuine gozinta chain has no empty step. If exactly \(j\) positions among \(k\) are active, the exponent assignments are counted by \(T_j(\mathbf a)\), and there are
$$\binom{k}{j}$$
ways to choose those active positions. Inclusion-exclusion over the empty positions yields the number of chains with exactly \(k\) nontrivial steps:
$$F_k(\mathbf a)=\sum_{j=1}^{k}(-1)^{k-j}\binom{k}{j}T_j(\mathbf a).$$
This alternating sum is the core combinatorial formula in the C++, Python, and Java implementations.
Each nontrivial step increases the total exponent sum by at least \(1\), so the number of steps cannot exceed
$$\omega=a_1+\cdots+a_r.$$
Therefore the total number of gozinta chains for pattern \(\mathbf a\) is
$$G(\mathbf a)=\sum_{k=1}^{\omega}F_k(\mathbf a).$$
At this point the original divisor-chain problem has been reduced to exact finite sums of binomial coefficients.
Let
$$m=G(\mathbf a).$$
If the exponent pattern of \(m\) is again \(\mathbf a\), then \(m\) itself has pattern \(\mathbf a\), so
$$g(m)=G(\mathbf a)=m.$$
Conversely, every fixed point \(n\) produces its own exponent pattern and passes this test. So the search does not need to enumerate all integers; it only needs to enumerate exponent patterns, compute \(G(\mathbf a)\), factor the resulting integer, and compare patterns.
Among all integers having exponent pattern \(\mathbf a\), the smallest one is obtained by matching the largest exponent with the smallest prime, the next largest exponent with the next prime, and so on:
$$n_{\min}(\mathbf a)=2^{a_1}3^{a_2}5^{a_3}\cdots.$$
If
$$n_{\min}(\mathbf a)\gt 10^{16},$$
then no integer with that pattern can contribute and the entire branch is discarded. There is a second useful rejection test:
$$G(\mathbf a)\lt n_{\min}(\mathbf a)\implies \text{pattern}(G(\mathbf a))\ne \mathbf a.$$
So factorization is only attempted when the candidate is large enough even to have the required pattern.
Take \(\mathbf a=(4,1)\), the exponent pattern of \(48=2^4\cdot 3\). Here \(\omega=5\). First compute
$$\begin{aligned} T_1&=\binom{4}{0}\binom{1}{0}=1,\\ T_2&=\binom{5}{1}\binom{2}{1}=10,\\ T_3&=\binom{6}{2}\binom{3}{2}=45,\\ T_4&=\binom{7}{3}\binom{4}{3}=140,\\ T_5&=\binom{8}{4}\binom{5}{4}=350. \end{aligned}$$
Then inclusion-exclusion gives
$$\begin{aligned} F_1&=1,\\ F_2&=-2T_1+T_2=8,\\ F_3&=3T_1-3T_2+T_3=18,\\ F_4&=-4T_1+6T_2-4T_3+T_4=16,\\ F_5&=5T_1-10T_2+10T_3-5T_4+T_5=5. \end{aligned}$$
So
$$G(4,1)=1+8+18+16+5=48.$$
Since \(48\) has exponent pattern \((4,1)\), it is a true fixed point and must be included in the final sum. By contrast, \((2,1)\) gives \(G(2,1)=8\), whose pattern is \((3)\), so \(12\) is not a fixed point.
The C++, Python, and Java implementations recursively enumerate exponent patterns in non-increasing order. The search starts with a safe upper bound on the first exponent and extends the pattern one prime at a time. Every recursive step updates the smallest possible integer with that pattern, so branches exceeding \(10^{16}\) are cut immediately.
For each surviving pattern, the implementation builds the needed binomial coefficients exactly, evaluates the product terms \(T_k\), and applies the inclusion-exclusion formula to obtain \(F_k\) and then \(G(\mathbf a)\). The running total is stopped as soon as it exceeds the global limit, because larger values can never satisfy the problem condition.
Only candidates \(m\) with \(m\ge n_{\min}(\mathbf a)\) are factored. Their prime exponents are extracted, sorted in non-increasing order, and compared with the original pattern. Previously analyzed candidates are cached, valid fixed points are deduplicated, and the final sum is accumulated with wide integer arithmetic.
Let \(P\) be the number of exponent patterns whose smallest representative does not exceed \(10^{16}\). For a pattern \(\mathbf a\) with weight \(\omega=\sum a_i\) and length \(r\), evaluating \(G(\mathbf a)\) costs \(O(\omega^2+r\omega)\) arithmetic operations: \(O(\omega^2)\) for the binomial table and the inclusion-exclusion layers, and \(O(r\omega)\) for the products \(T_k\).
The recursive enumeration is practical because the lower-bound pruning keeps \(P\) small. The factorization stage uses deterministic Miller-Rabin primality testing for 64-bit integers together with Pollard-Rho splitting. That stage has no simple closed worst-case formula, but below \(10^{16}\) it is fast in practice and runs only on the comparatively small set of surviving candidates. Memory usage is modest: \(O(\omega^2)\) while processing one pattern, plus the cache of analyzed candidates and the set of discovered fixed points.
Eine Gozinta-Kette zu \(n\) ist eine Folge
$$1=d_0 \lt d_1 \lt \cdots \lt d_k=n$$
mit der Eigenschaft, dass jedes Glied das folgende teilt. Sei \(g(n)\) die Anzahl dieser Ketten. Gesucht ist die Summe aller Zahlen \(n\le 10^{16}\) mit
$$g(n)=n.$$
Die entscheidende Beobachtung lautet: \(g(n)\) hängt nur vom sortierten Muster der Primfaktor-Exponenten ab, nicht von den konkreten Primzahlen. Damit wird die Suche über alle Zahlen bis \(10^{16}\) auf eine viel kleinere Suche im Raum der Exponentenmuster reduziert.
Schreibe die Primfaktorzerlegung mit nichtzunehmenden Exponenten:
$$n=\prod_{i=1}^{r} p_i^{a_i},\qquad a_1\ge a_2\ge \cdots \ge a_r \gt 0.$$
Sei \(\mathbf a=(a_1,\dots,a_r)\) und \(\omega=a_1+\cdots+a_r\).
Jeder Teiler von \(n\) hat die Form
$$\prod_{i=1}^{r} p_i^{b_i},\qquad 0\le b_i\le a_i.$$
Eine Gozinta-Kette entspricht also genau einer koordinatenweise nichtfallenden Folge von Exponentenvektoren, die bei \((0,\dots,0)\) beginnt und bei \((a_1,\dots,a_r)\) endet.
Da diese Beschreibung nur die Exponenten benutzt, wird die Kettenzahl vollständig durch \(\mathbf a\) bestimmt. Deshalb schreiben wir
$$g(n)=G(\mathbf a).$$
Das ist die zentrale Vereinfachung der gesamten Lösung.
Angenommen, die Kette besitzt genau \(k\) nichttriviale Teilbarkeitsschritte:
$$1=d_0 \lt d_1 \lt \cdots \lt d_k=n.$$
Definiere die Quotienten \(q_t=d_t/d_{t-1}\). Dann gilt \(q_t \gt 1\) und
$$q_1q_2\cdots q_k=n.$$
Für einen Primfaktor \(p_i^{a_i}\) verteilen wir den Exponenten \(a_i\) auf diese \(k\) Quotienten:
$$e_{i,1}+e_{i,2}+\cdots+e_{i,k}=a_i,\qquad e_{i,t}\ge 0.$$
Die Anzahl schwacher Kompositionen von \(a_i\) in \(k\) Teile ist
$$\binom{a_i+k-1}{k-1}.$$
Wegen der Unabhängigkeit der Primzahlen folgt
$$T_k(\mathbf a)=\prod_{i=1}^{r}\binom{a_i+k-1}{k-1}.$$
Diese Größe zählt geordnete Zerlegungen in \(k\) Schritte, erlaubt aber noch Quotienten gleich \(1\) und zählt daher zu viel.
Eine echte Gozinta-Kette darf keinen leeren Schritt enthalten. Sind unter \(k\) Positionen genau \(j\) aktiv, dann werden die Exponentenzuweisungen durch \(T_j(\mathbf a)\) gezählt, und die aktiven Positionen können auf
$$\binom{k}{j}$$
Arten gewählt werden. Inklusion-Exklusion über die leeren Positionen liefert daher die Anzahl der Ketten mit genau \(k\) nichttrivialen Schritten:
$$F_k(\mathbf a)=\sum_{j=1}^{k}(-1)^{k-j}\binom{k}{j}T_j(\mathbf a).$$
Genau diese alternierende Summe wird in den drei Implementierungen ausgewertet.
Jeder nichttriviale Schritt erhöht die gesamte Exponentensumme mindestens um \(1\), also kann die Zahl der Schritte höchstens
$$\omega=a_1+\cdots+a_r$$
betragen. Damit ist die Gesamtzahl der Gozinta-Ketten zum Muster \(\mathbf a\)
$$G(\mathbf a)=\sum_{k=1}^{\omega}F_k(\mathbf a).$$
Damit ist das ursprüngliche Divisorkettenproblem auf endliche Summen aus Binomialkoeffizienten reduziert.
Setze
$$m=G(\mathbf a).$$
Hat \(m\) wieder das Exponentenmuster \(\mathbf a\), dann besitzt \(m\) selbst dieses Muster und somit
$$g(m)=G(\mathbf a)=m.$$
Umgekehrt erzeugt jeder Fixpunkt \(n\) sein eigenes Exponentenmuster und besteht genau diesen Test. Deshalb muss die Suche nicht alle Zahlen durchlaufen; es genügt, Exponentenmuster zu erzeugen, \(G(\mathbf a)\) zu berechnen, \(m\) zu faktorisieren und die Muster zu vergleichen.
Unter allen Zahlen mit Muster \(\mathbf a\) ist die kleinste diejenige, bei der der größte Exponent an die kleinste Primzahl, der zweitgrößte an die zweitkleinste Primzahl usw. gekoppelt wird:
$$n_{\min}(\mathbf a)=2^{a_1}3^{a_2}5^{a_3}\cdots.$$
Wenn
$$n_{\min}(\mathbf a)\gt 10^{16},$$
dann kann keine Zahl mit diesem Muster beitragen und der ganze Ast wird sofort verworfen. Zusätzlich gilt:
$$G(\mathbf a)\lt n_{\min}(\mathbf a)\implies \text{Muster}(G(\mathbf a))\ne \mathbf a.$$
In diesem Fall ist also nicht einmal eine Faktorisierung nötig.
Betrachte \(\mathbf a=(4,1)\), also das Exponentenmuster von \(48=2^4\cdot 3\). Hier ist \(\omega=5\). Zunächst erhält man
$$\begin{aligned} T_1&=\binom{4}{0}\binom{1}{0}=1,\\ T_2&=\binom{5}{1}\binom{2}{1}=10,\\ T_3&=\binom{6}{2}\binom{3}{2}=45,\\ T_4&=\binom{7}{3}\binom{4}{3}=140,\\ T_5&=\binom{8}{4}\binom{5}{4}=350. \end{aligned}$$
Mit Inklusion-Exklusion folgt
$$\begin{aligned} F_1&=1,\\ F_2&=-2T_1+T_2=8,\\ F_3&=3T_1-3T_2+T_3=18,\\ F_4&=-4T_1+6T_2-4T_3+T_4=16,\\ F_5&=5T_1-10T_2+10T_3-5T_4+T_5=5. \end{aligned}$$
Daher
$$G(4,1)=1+8+18+16+5=48.$$
Da \(48\) selbst das Muster \((4,1)\) hat, ist es ein echter Fixpunkt und gehört in die Endsumme. Dagegen liefert \((2,1)\) den Wert \(G(2,1)=8\); dessen Muster ist \((3)\), also ist \(12\) kein Fixpunkt.
Die C++-, Python- und Java-Implementierungen erzeugen Exponentenmuster rekursiv in nichtzunehmender Reihenfolge. Die Suche startet mit einer sicheren Obergrenze für den ersten Exponenten und erweitert das Muster Primzahl für Primzahl. Bei jedem Schritt wird die kleinste mögliche Zahl mit diesem Muster fortgeschrieben; überschreitet sie \(10^{16}\), wird der Ast sofort abgeschnitten.
Für jedes überlebende Muster berechnet die Implementierung die benötigten Binomialkoeffizienten exakt, bildet daraus die Produkte \(T_k\) und wertet anschließend die Inklusion-Exklusion-Formel für \(F_k\) und schließlich \(G(\mathbf a)\) aus. Sobald die laufende Summe über die globale Schranke hinauswächst, wird abgebrochen, weil größere Werte ohnehin nicht mehr als Lösungen infrage kommen.
Nur Kandidaten \(m\) mit \(m\ge n_{\min}(\mathbf a)\) werden faktorisiert. Deren Primexponenten werden extrahiert, absteigend sortiert und mit dem ursprünglichen Muster verglichen. Bereits analysierte Kandidaten werden zwischengespeichert, gültige Fixpunkte werden dedupliziert und die Endsumme wird mit breiter Ganzzahlarithmetik akkumuliert.
Sei \(P\) die Anzahl der Exponentenmuster, deren kleinster Repräsentant \(10^{16}\) nicht überschreitet. Für ein Muster \(\mathbf a\) mit Gewicht \(\omega=\sum a_i\) und Länge \(r\) kostet die Auswertung von \(G(\mathbf a)\) insgesamt \(O(\omega^2+r\omega)\) arithmetische Operationen: \(O(\omega^2)\) für die Binomialtabelle und die Inklusion-Exklusion, sowie \(O(r\omega)\) für die Produktterme \(T_k\).
Die rekursive Enumeration bleibt praktikabel, weil das Abschneiden über die Untergrenze \(P\) stark reduziert. Die Faktorisierung verwendet einen deterministischen Miller-Rabin-Primzahltest für 64-Bit-Zahlen zusammen mit Pollard-Rho-Zerlegung. Dafür gibt es keine einfache geschlossene Worst-Case-Formel, doch unterhalb von \(10^{16}\) ist das Verfahren in der Praxis schnell und läuft nur auf einer vergleichsweise kleinen Kandidatenmenge. Der Speicherbedarf bleibt moderat: \(O(\omega^2)\) während der Auswertung eines Musters, plus Cache für analysierte Kandidaten und Menge der gefundenen Fixpunkte.
Bir \(n\) sayısı için gozinta zinciri,
$$1=d_0 \lt d_1 \lt \cdots \lt d_k=n$$
şeklinde, her terimin bir sonrakini böldüğü bir dizidir. \(g(n)\), bu zincirlerin sayısı olsun. Aranan şey,
$$g(n)=n$$
koşulunu sağlayan tüm \(n\le 10^{16}\) değerlerinin toplamıdır.
Çözümün temel gözlemi şudur: \(g(n)\), kullanılan asal sayıların kendisine değil, yalnızca \(n\)'nin sıralanmış asal üs desenine bağlıdır. Böylece \(10^{16}\)'ya kadar tüm sayıları taramak yerine çok daha küçük bir üs deseni uzayında arama yapılır.
Asal çarpanlara ayrılışı üsler azalmayan sırada yazalım:
$$n=\prod_{i=1}^{r} p_i^{a_i},\qquad a_1\ge a_2\ge \cdots \ge a_r \gt 0.$$
\(\mathbf a=(a_1,\dots,a_r)\) ve \(\omega=a_1+\cdots+a_r\) olsun.
\(n\)'nin her böleni
$$\prod_{i=1}^{r} p_i^{b_i},\qquad 0\le b_i\le a_i$$
şeklindedir. Bu yüzden bir gozinta zinciri, \((0,\dots,0)\)'dan başlayıp \((a_1,\dots,a_r)\)'de biten, koordinat bazında azalmayan üs vektörleri dizisine denktir.
Bu tanım yalnızca üslere baktığı için zincir sayısı tamamen \(\mathbf a\) tarafından belirlenir. Dolayısıyla
$$g(n)=G(\mathbf a)$$
yazabiliriz. Uygulamanın bütün gücü bu indirgemeden gelir.
Zincirin tam olarak \(k\) tane önemsiz olmayan bölünebilme adımı içerdiğini varsayalım:
$$1=d_0 \lt d_1 \lt \cdots \lt d_k=n.$$
Adım oranlarını \(q_t=d_t/d_{t-1}\) ile tanımlayalım. Her \(q_t \gt 1\) olur ve
$$q_1q_2\cdots q_k=n.$$
Tek bir asal kuvvet \(p_i^{a_i}\) için, \(a_i\) üssünü bu \(k\) orana dağıtalım:
$$e_{i,1}+e_{i,2}+\cdots+e_{i,k}=a_i,\qquad e_{i,t}\ge 0.$$
\(a_i\)'nin \(k\) parçaya zayıf bileşim sayısı
$$\binom{a_i+k-1}{k-1}$$
olur. Asallar birbirinden bağımsız olduğu için
$$T_k(\mathbf a)=\prod_{i=1}^{r}\binom{a_i+k-1}{k-1}$$
elde edilir. Fakat bu sayı bazı oranların \(1\) olmasına izin verdiği için gerçek zincirleri fazla sayar.
Gerçek bir gozinta zincirinde boş adım olamaz. \(k\) konumun tam olarak \(j\) tanesi aktifse, üs atamaları \(T_j(\mathbf a)\) ile sayılır ve aktif konumlar
$$\binom{k}{j}$$
farklı biçimde seçilebilir. Boş konumlar üzerinde dahil et-çıkar uygularsak, tam olarak \(k\) önemsiz olmayan adım içeren zincirlerin sayısı
$$F_k(\mathbf a)=\sum_{j=1}^{k}(-1)^{k-j}\binom{k}{j}T_j(\mathbf a)$$
olur. Üç dildeki uygulamanın merkezindeki dönüşümlü toplam budur.
Önemsiz olmayan her adım, toplam üs sayısını en az \(1\) artırır. Bu yüzden adım sayısı en fazla
$$\omega=a_1+\cdots+a_r$$
olabilir. Böylece \(\mathbf a\) deseni için toplam gozinta zinciri sayısı
$$G(\mathbf a)=\sum_{k=1}^{\omega}F_k(\mathbf a)$$
şeklindedir. Yani problem, açıkça bölenleri üretmeden sonlu binom toplamlarına indirgenmiş olur.
Şimdi
$$m=G(\mathbf a)$$
olsun. Eğer \(m\)'nin üs deseni yine \(\mathbf a\) ise, \(m\) zaten bu desene sahip bir sayıdır ve
$$g(m)=G(\mathbf a)=m$$
elde edilir. Tersi de geçerlidir: her sabit nokta kendi üs desenini üretir ve bu testi geçer. Dolayısıyla tüm sayıları denemek yerine üs desenleri üretilir, \(G(\mathbf a)\) hesaplanır, çıkan sayı çarpanlara ayrılır ve desenler karşılaştırılır.
\(\mathbf a\) desenine sahip sayılar arasında en küçük olan, en büyük üssü en küçük asala, ikinci büyük üssü ikinci küçük asala eşleyerek elde edilir:
$$n_{\min}(\mathbf a)=2^{a_1}3^{a_2}5^{a_3}\cdots.$$
Eğer
$$n_{\min}(\mathbf a)\gt 10^{16},$$
ise bu desenle hiçbir çözüm üretilemez ve dal tamamen kesilir. Bir ikinci hızlı eleme daha vardır:
$$G(\mathbf a)\lt n_{\min}(\mathbf a)\implies \text{desen}(G(\mathbf a))\ne \mathbf a.$$
Bu durumda çarpanlara ayırma yapmaya bile gerek kalmaz.
\(48=2^4\cdot 3\) sayısının deseni olan \(\mathbf a=(4,1)\) için \(\omega=5\) olur. Önce
$$\begin{aligned} T_1&=\binom{4}{0}\binom{1}{0}=1,\\ T_2&=\binom{5}{1}\binom{2}{1}=10,\\ T_3&=\binom{6}{2}\binom{3}{2}=45,\\ T_4&=\binom{7}{3}\binom{4}{3}=140,\\ T_5&=\binom{8}{4}\binom{5}{4}=350. \end{aligned}$$
hesaplanır. Dahil et-çıkar sonrası
$$\begin{aligned} F_1&=1,\\ F_2&=-2T_1+T_2=8,\\ F_3&=3T_1-3T_2+T_3=18,\\ F_4&=-4T_1+6T_2-4T_3+T_4=16,\\ F_5&=5T_1-10T_2+10T_3-5T_4+T_5=5. \end{aligned}$$
bulunur. O halde
$$G(4,1)=1+8+18+16+5=48.$$
\(48\)'in üs deseni gerçekten \((4,1)\) olduğu için bu sayı çözümdür ve toplamda yer alır. Buna karşılık \((2,1)\) için \(G(2,1)=8\) çıkar; \(8\)'in deseni \((3)\) olduğu için \(12\) sabit nokta değildir.
C++, Python ve Java uygulamaları üs desenlerini azalmayan biçimde özyineli olarak üretir. Arama, ilk üs için güvenli bir üst sınırla başlar ve desen her adımda bir sonraki asalın kuvveti eklenerek büyütülür. O anki desenin alabileceği en küçük sayı \(10^{16}\)'yı geçerse dal hemen kapatılır.
Hayatta kalan her desen için uygulama gerekli binom katsayılarını tam olarak üretir, \(T_k\) çarpımlarını hesaplar ve ardından dahil et-çıkar formülüyle önce \(F_k\), sonra da \(G(\mathbf a)\) değerini bulur. Kısmi toplam sınırı aştığı anda işlem durdurulur; daha büyük değerlerin çözüm olma ihtimali yoktur.
Yalnızca \(m\ge n_{\min}(\mathbf a)\) koşulunu sağlayan adaylar çarpanlara ayrılır. Elde edilen asal üsleri büyükten küçüğe sıralanır ve başlangıç deseniyle karşılaştırılır. Daha önce analiz edilmiş adaylar önbelleğe alınır, geçerli sabit noktalar tekilleştirilir ve son toplam geniş tamsayı aritmetiğiyle biriktirilir.
\(P\), en küçük temsilcisi \(10^{16}\)'yı aşmayan üs deseni sayısı olsun. Ağırlığı \(\omega=\sum a_i\), uzunluğu \(r\) olan bir \(\mathbf a\) deseni için \(G(\mathbf a)\)'nın hesaplanması \(O(\omega^2+r\omega)\) aritmetik işlem gerektirir: binom tablosu ve dahil et-çıkar katmanları için \(O(\omega^2)\), \(T_k\) çarpımları için ise \(O(r\omega)\).
Özyineli arama pratikte etkilidir çünkü alt sınır budaması \(P\)'yi küçük tutar. Çarpanlara ayırma aşaması, 64 bit tamsayılar için deterministik Miller-Rabin asal testi ile Pollard-Rho ayrıştırmasını birlikte kullanır. Bunun basit bir kapalı en kötü durum formülü yoktur; ancak \(10^{16}\) altında oldukça hızlıdır ve yalnızca ayakta kalan az sayıdaki aday üzerinde çalışır. Bellek kullanımı düşüktür: tek bir desen işlenirken \(O(\omega^2)\), ayrıca analiz edilen adayların önbelleği ve bulunan sabit noktaların kümesi kadar ek alan gerekir.
Una cadena gozinta de \(n\) es una sucesión
$$1=d_0 \lt d_1 \lt \cdots \lt d_k=n$$
en la que cada término divide al siguiente. Sea \(g(n)\) el número de esas cadenas. Hay que sumar todos los enteros \(n\le 10^{16}\) tales que
$$g(n)=n.$$
La observación decisiva es que \(g(n)\) depende solo del patrón ordenado de exponentes primos de \(n\), no de cuáles sean los primos concretos. Eso convierte una búsqueda imposible sobre todos los enteros hasta \(10^{16}\) en una búsqueda manejable sobre patrones de exponentes.
Escribamos la factorización prima con exponentes en orden no creciente:
$$n=\prod_{i=1}^{r} p_i^{a_i},\qquad a_1\ge a_2\ge \cdots \ge a_r \gt 0.$$
Sea \(\mathbf a=(a_1,\dots,a_r)\) y \(\omega=a_1+\cdots+a_r\).
Cada divisor de \(n\) tiene la forma
$$\prod_{i=1}^{r} p_i^{b_i},\qquad 0\le b_i\le a_i.$$
Por tanto, una cadena gozinta equivale a una sucesión de vectores de exponentes no decrecientes coordenada a coordenada, desde \((0,\dots,0)\) hasta \((a_1,\dots,a_r)\).
Como esta descripción solo usa los exponentes, el número de cadenas queda completamente determinado por \(\mathbf a\). Podemos escribir entonces
$$g(n)=G(\mathbf a).$$
Esa es la reducción fundamental que explota la implementación.
Supongamos que la cadena tiene exactamente \(k\) pasos de divisibilidad no triviales:
$$1=d_0 \lt d_1 \lt \cdots \lt d_k=n.$$
Definimos las razones \(q_t=d_t/d_{t-1}\). Entonces cada \(q_t \gt 1\) y
$$q_1q_2\cdots q_k=n.$$
Para un primo \(p_i^{a_i}\), distribuimos el exponente \(a_i\) entre esas \(k\) razones:
$$e_{i,1}+e_{i,2}+\cdots+e_{i,k}=a_i,\qquad e_{i,t}\ge 0.$$
El número de composiciones débiles de \(a_i\) en \(k\) partes es
$$\binom{a_i+k-1}{k-1}.$$
Como los distintos primos actúan de forma independiente, obtenemos
$$T_k(\mathbf a)=\prod_{i=1}^{r}\binom{a_i+k-1}{k-1}.$$
Sin embargo, \(T_k(\mathbf a)\) todavía permite que algunas razones sean \(1\), así que cuenta de más.
Una cadena gozinta auténtica no puede tener pasos vacíos. Si entre \(k\) posiciones exactamente \(j\) están activas, las asignaciones de exponentes se cuentan con \(T_j(\mathbf a)\), y hay
$$\binom{k}{j}$$
formas de elegir esas posiciones activas. Aplicando inclusión-exclusión sobre las posiciones vacías se obtiene el número de cadenas con exactamente \(k\) pasos no triviales:
$$F_k(\mathbf a)=\sum_{j=1}^{k}(-1)^{k-j}\binom{k}{j}T_j(\mathbf a).$$
Esta suma alternante es la pieza combinatoria central en las tres implementaciones.
Cada paso no trivial incrementa la suma total de exponentes al menos en \(1\), así que el número de pasos no puede superar
$$\omega=a_1+\cdots+a_r.$$
Por tanto, el número total de cadenas gozinta del patrón \(\mathbf a\) es
$$G(\mathbf a)=\sum_{k=1}^{\omega}F_k(\mathbf a).$$
Con ello, el problema original de cadenas de divisores queda reducido a sumas finitas exactas de coeficientes binomiales.
Sea
$$m=G(\mathbf a).$$
Si el patrón de exponentes de \(m\) vuelve a ser \(\mathbf a\), entonces \(m\) mismo tiene ese patrón y por tanto
$$g(m)=G(\mathbf a)=m.$$
Recíprocamente, todo punto fijo \(n\) genera su propio patrón y supera exactamente esta prueba. Así, la búsqueda se hace en el espacio de patrones: se genera \(\mathbf a\), se calcula \(G(\mathbf a)\), se factoriza el entero obtenido y se comparan los patrones.
Entre todos los enteros con patrón \(\mathbf a\), el menor se obtiene asignando el mayor exponente al primo más pequeño, el segundo mayor al segundo primo más pequeño, etc.:
$$n_{\min}(\mathbf a)=2^{a_1}3^{a_2}5^{a_3}\cdots.$$
Si
$$n_{\min}(\mathbf a)\gt 10^{16},$$
entonces ningún entero con ese patrón puede contribuir y toda la rama se descarta. Además,
$$G(\mathbf a)\lt n_{\min}(\mathbf a)\implies \text{patrón}(G(\mathbf a))\ne \mathbf a.$$
En ese caso ni siquiera hace falta factorizar el candidato.
Tomemos \(\mathbf a=(4,1)\), el patrón de exponentes de \(48=2^4\cdot 3\). Aquí \(\omega=5\). Primero se calcula
$$\begin{aligned} T_1&=\binom{4}{0}\binom{1}{0}=1,\\ T_2&=\binom{5}{1}\binom{2}{1}=10,\\ T_3&=\binom{6}{2}\binom{3}{2}=45,\\ T_4&=\binom{7}{3}\binom{4}{3}=140,\\ T_5&=\binom{8}{4}\binom{5}{4}=350. \end{aligned}$$
Luego, por inclusión-exclusión,
$$\begin{aligned} F_1&=1,\\ F_2&=-2T_1+T_2=8,\\ F_3&=3T_1-3T_2+T_3=18,\\ F_4&=-4T_1+6T_2-4T_3+T_4=16,\\ F_5&=5T_1-10T_2+10T_3-5T_4+T_5=5. \end{aligned}$$
Por tanto,
$$G(4,1)=1+8+18+16+5=48.$$
Como \(48\) tiene precisamente el patrón \((4,1)\), es un punto fijo real y entra en la suma final. En cambio, \((2,1)\) produce \(G(2,1)=8\), cuyo patrón es \((3)\), así que \(12\) no es un punto fijo.
Las implementaciones en C++, Python y Java enumeran recursivamente los patrones de exponentes en orden no creciente. La búsqueda comienza con una cota superior segura para el primer exponente y extiende el patrón usando potencias de primos sucesivos. Cada paso actualiza el menor entero posible compatible con el patrón actual, de modo que cualquier rama que ya exceda \(10^{16}\) se corta inmediatamente.
Para cada patrón que sobrevive, la implementación construye exactamente los coeficientes binomiales necesarios, calcula los productos \(T_k\) y aplica la fórmula de inclusión-exclusión para obtener \(F_k\) y después \(G(\mathbf a)\). Si la suma parcial ya supera el límite global, se detiene en ese punto porque los valores mayores ya no pueden ser soluciones.
Solo se factorizan candidatos \(m\) con \(m\ge n_{\min}(\mathbf a)\). Sus exponentes primos se extraen, se ordenan de mayor a menor y se comparan con el patrón original. Los candidatos ya analizados se guardan en caché, los puntos fijos válidos se deduplican y la suma final se acumula con aritmética entera amplia.
Sea \(P\) el número de patrones de exponentes cuyo representante mínimo no supera \(10^{16}\). Para un patrón \(\mathbf a\) de peso \(\omega=\sum a_i\) y longitud \(r\), evaluar \(G(\mathbf a)\) cuesta \(O(\omega^2+r\omega)\) operaciones aritméticas: \(O(\omega^2)\) para la tabla binomial y las capas de inclusión-exclusión, y \(O(r\omega)\) para los productos \(T_k\).
La enumeración recursiva es viable porque la poda mediante la cota inferior mantiene pequeño el valor de \(P\). La fase de factorización usa la prueba de primalidad determinista de Miller-Rabin para enteros de 64 bits junto con descomposición Pollard-Rho. Esa parte no tiene una cota cerrada simple en el peor caso, pero para números por debajo de \(10^{16}\) es rápida en la práctica y solo se aplica al conjunto relativamente pequeño de candidatos supervivientes. La memoria utilizada es moderada: \(O(\omega^2)\) al procesar un patrón, más la caché de candidatos analizados y el conjunto de puntos fijos descubiertos.
对整数 \(n\)，gozinta 链是一个序列
$$1=d_0 \lt d_1 \lt \cdots \lt d_k=n$$
其中每一项都整除下一项。记 \(g(n)\) 为这样的链的条数。题目要求求出所有满足
$$g(n)=n$$
且 \(n\le 10^{16}\) 的整数之和。
本题的关键并不是直接枚举整数，而是先看素因子指数模式。实现所利用的核心事实是：\(g(n)\) 只依赖于 \(n\) 的素因子指数按降序排列后的模式，而不依赖于具体用了哪些素数。这样就把原本几乎无法处理的整数搜索，变成了对指数模式的递归枚举。
把 \(n\) 的素因子分解写成指数非增的形式：
$$n=\prod_{i=1}^{r} p_i^{a_i},\qquad a_1\ge a_2\ge \cdots \ge a_r \gt 0.$$
记指数模式为 \(\mathbf a=(a_1,\dots,a_r)\)，再记
$$\omega=a_1+\cdots+a_r.$$
\(n\) 的任意一个因子都可以写成
$$\prod_{i=1}^{r} p_i^{b_i},\qquad 0\le b_i\le a_i.$$
因此，一条 gozinta 链本质上就是一串指数向量：从 \((0,\dots,0)\) 出发，按坐标逐项不减地走到 \((a_1,\dots,a_r)\)。
这个描述完全没有用到具体素数的数值，只用到了指数，所以链的数量只由 \(\mathbf a\) 决定。于是可以写成
$$g(n)=G(\mathbf a).$$
这一步是整个解法最重要的降维：从“按整数搜索”变成“按指数模式搜索”。
设一条链恰好有 \(k\) 个真正的整除步骤：
$$1=d_0 \lt d_1 \lt \cdots \lt d_k=n.$$
定义相邻两项的比值
$$q_t=d_t/d_{t-1}.$$
那么每个 \(q_t \gt 1\)，并且满足
$$q_1q_2\cdots q_k=n.$$
对某个固定的素因子 \(p_i^{a_i}\)，把指数 \(a_i\) 分配到这 \(k\) 个比值上：
$$e_{i,1}+e_{i,2}+\cdots+e_{i,k}=a_i,\qquad e_{i,t}\ge 0.$$
这就是把 \(a_i\) 分成 \(k\) 个允许为零的部分，其方案数为
$$\binom{a_i+k-1}{k-1}.$$
不同素因子之间相互独立，所以得到
$$T_k(\mathbf a)=\prod_{i=1}^{r}\binom{a_i+k-1}{k-1}.$$
但是 \(T_k(\mathbf a)\) 还允许某些比值 \(q_t=1\)，也就是某些步骤实际上什么都没发生，因此它统计的是“带空步骤”的过大计数。
真正的 gozinta 链要求每一步都必须非平凡。若在总共 \(k\) 个位置中只有 \(j\) 个位置是真正活跃的，那么相应的指数分配由 \(T_j(\mathbf a)\) 给出，而活跃位置本身可以从 \(k\) 个位置里选出
$$\binom{k}{j}$$
种。对“空步骤”的集合应用容斥原理，可得恰好含 \(k\) 个非平凡步骤的链数：
$$F_k(\mathbf a)=\sum_{j=1}^{k}(-1)^{k-j}\binom{k}{j}T_j(\mathbf a).$$
这就是三份实现共同使用的核心交错求和公式。
每一个非平凡步骤至少会让总指数和增加 \(1\)，因此步骤数不可能超过
$$\omega=a_1+\cdots+a_r.$$
于是模式 \(\mathbf a\) 的 gozinta 链总数为
$$G(\mathbf a)=\sum_{k=1}^{\omega}F_k(\mathbf a).$$
到这里，原问题已经被改写成有限个二项式系数和容斥求和，没有显式枚举因子，也没有显式枚举链。
令
$$m=G(\mathbf a).$$
如果 \(m\) 的指数模式恰好仍然是 \(\mathbf a\)，那么 \(m\) 自身就具有模式 \(\mathbf a\)，从而
$$g(m)=G(\mathbf a)=m.$$
反过来，任何满足 \(g(n)=n\) 的固定点 \(n\) 也一定会产生自己的指数模式并通过这个判定。所以搜索过程可以完全在“模式空间”里进行：枚举一个模式 \(\mathbf a\)，算出 \(m=G(\mathbf a)\)，再把 \(m\) 分解质因数并比较其指数模式是否回到 \(\mathbf a\)。
在所有具有模式 \(\mathbf a\) 的整数里，最小的那个一定是把最大的指数配给最小的素数、第二大的指数配给第二小的素数，依此类推：
$$n_{\min}(\mathbf a)=2^{a_1}3^{a_2}5^{a_3}\cdots.$$
如果
$$n_{\min}(\mathbf a)\gt 10^{16},$$
那么具有该模式的任何整数都不可能成为答案，整棵递归分支可以直接丢弃。还有一个同样重要的快速过滤：
$$G(\mathbf a)\lt n_{\min}(\mathbf a)\implies \text{pattern}(G(\mathbf a))\ne \mathbf a.$$
也就是说，若候选值甚至比该模式的最小代表还小，那么它绝不可能拥有这个模式，这时就没有必要再做质因数分解。
取模式 \(\mathbf a=(4,1)\)，它对应于 \(48=2^4\cdot 3\) 这样的整数。此时 \(\omega=5\)。先计算
$$\begin{aligned} T_1&=\binom{4}{0}\binom{1}{0}=1,\\ T_2&=\binom{5}{1}\binom{2}{1}=10,\\ T_3&=\binom{6}{2}\binom{3}{2}=45,\\ T_4&=\binom{7}{3}\binom{4}{3}=140,\\ T_5&=\binom{8}{4}\binom{5}{4}=350. \end{aligned}$$
再做容斥：
$$\begin{aligned} F_1&=1,\\ F_2&=-2T_1+T_2=8,\\ F_3&=3T_1-3T_2+T_3=18,\\ F_4&=-4T_1+6T_2-4T_3+T_4=16,\\ F_5&=5T_1-10T_2+10T_3-5T_4+T_5=5. \end{aligned}$$
因此
$$G(4,1)=1+8+18+16+5=48.$$
因为 \(48\) 的指数模式确实就是 \((4,1)\)，所以它是一个真正的固定点，必须计入最终答案。作为对照，\((2,1)\) 会得到 \(G(2,1)=8\)，而 \(8\) 的模式是 \((3)\)，因此 \(12\) 不是固定点。
C++、Python 和 Java 实现都会按指数非增的顺序递归生成模式。搜索从一个足够安全的首项指数上界开始，每向下一层递归就把下一种素数的若干次幂并入当前模式。与此同时，程序持续维护该模式能够对应的最小整数；一旦这个下界已经超过 \(10^{16}\)，整条分支立刻终止。
对于每个保留下来的模式，实现会精确构造所需的二项式系数，计算各个 \(T_k\) 的乘积，再套用容斥公式得到 \(F_k\)，最后累加成 \(G(\mathbf a)\)。如果中间的累计值已经超过全局上限，就立即停止，因为更大的值不可能再满足题意。
只有满足 \(m\ge n_{\min}(\mathbf a)\) 的候选 \(m\) 才会进入质因数分解阶段。分解后提取出的指数按非增顺序排序，并与原模式比较。实现还会缓存已经分析过的候选值，最后对所有有效固定点去重，并用宽整数类型完成求和。
设 \(P\) 为所有最小代表不超过 \(10^{16}\) 的指数模式数量。对于一个权重为 \(\omega=\sum a_i\)、长度为 \(r\) 的模式 \(\mathbf a\)，计算 \(G(\mathbf a)\) 需要 \(O(\omega^2+r\omega)\) 次算术操作：其中 \(O(\omega^2)\) 来自二项式表与容斥层，\(O(r\omega)\) 来自各个 \(T_k\) 乘积的构造。
递归搜索之所以可行，是因为最小代表下界的剪枝把 \(P\) 控制得很小。质因数分解阶段使用适用于 64 位整数的确定性 Miller-Rabin 素性测试与 Pollard-Rho 分裂方法。这个阶段没有一个简单的闭式最坏情况公式，但在 \(10^{16}\) 以内运行非常快，而且只会作用于少量通过前面筛选的候选。空间开销也较温和：处理单个模式时为 \(O(\omega^2)\)，再加上候选缓存和已发现固定点集合的存储。
Gozinta-цепочка для числа \(n\) - это последовательность
$$1=d_0 \lt d_1 \lt \cdots \lt d_k=n$$
такая, что каждый член делит следующий. Обозначим через \(g(n)\) количество таких цепочек. Нужно найти сумму всех чисел \(n\le 10^{16}\), удовлетворяющих условию
$$g(n)=n.$$
Ключевое наблюдение состоит в том, что \(g(n)\) зависит не от самих простых чисел, входящих в разложение \(n\), а только от отсортированного шаблона их показателей. Поэтому вместо перебора всех чисел до \(10^{16}\) можно перебирать лишь шаблоны показателей.
Запишем разложение на простые множители в порядке невозрастания показателей:
$$n=\prod_{i=1}^{r} p_i^{a_i},\qquad a_1\ge a_2\ge \cdots \ge a_r \gt 0.$$
Обозначим шаблон через \(\mathbf a=(a_1,\dots,a_r)\), а сумму показателей через
$$\omega=a_1+\cdots+a_r.$$
Любой делитель числа \(n\) имеет вид
$$\prod_{i=1}^{r} p_i^{b_i},\qquad 0\le b_i\le a_i.$$
Следовательно, gozinta-цепочка эквивалентна последовательности векторов показателей, которые покоординатно не убывают и идут от \((0,\dots,0)\) к \((a_1,\dots,a_r)\).
В этом описании важны только показатели, а не сами простые числа. Поэтому число цепочек определяется только шаблоном \(\mathbf a\), и можно писать
$$g(n)=G(\mathbf a).$$
Именно это наблюдение делает задачу вычислимой.
Пусть цепочка имеет ровно \(k\) нетривиальных шагов делимости:
$$1=d_0 \lt d_1 \lt \cdots \lt d_k=n.$$
Введем отношения соседних членов \(q_t=d_t/d_{t-1}\). Тогда каждое \(q_t \gt 1\) и
$$q_1q_2\cdots q_k=n.$$
Для одного простого множителя \(p_i^{a_i}\) распределим показатель \(a_i\) по этим \(k\) отношениям:
$$e_{i,1}+e_{i,2}+\cdots+e_{i,k}=a_i,\qquad e_{i,t}\ge 0.$$
Количество слабых композиций числа \(a_i\) на \(k\) частей равно
$$\binom{a_i+k-1}{k-1}.$$
Поскольку простые множители независимы, получаем
$$T_k(\mathbf a)=\prod_{i=1}^{r}\binom{a_i+k-1}{k-1}.$$
Но величина \(T_k(\mathbf a)\) все еще допускает шаги с отношением \(1\), то есть считает также цепочки с пустыми шагами.
В настоящей gozinta-цепочке каждый шаг обязан быть нетривиальным. Если среди \(k\) позиций активны ровно \(j\), то распределения показателей считаются величиной \(T_j(\mathbf a)\), а сами активные позиции можно выбрать
$$\binom{k}{j}$$
способами. Применяя принцип включения-исключения по пустым позициям, получаем число цепочек с ровно \(k\) нетривиальными шагами:
$$F_k(\mathbf a)=\sum_{j=1}^{k}(-1)^{k-j}\binom{k}{j}T_j(\mathbf a).$$
Именно эта знакопеременная сумма реализована во всех трех версиях решения.
Каждый нетривиальный шаг увеличивает суммарный показатель как минимум на \(1\), поэтому число шагов не может превосходить
$$\omega=a_1+\cdots+a_r.$$
Следовательно, полное число gozinta-цепочек для шаблона \(\mathbf a\) равно
$$G(\mathbf a)=\sum_{k=1}^{\omega}F_k(\mathbf a).$$
Таким образом, задача о цепочках делителей сводится к конечным суммам биномиальных коэффициентов.
Положим
$$m=G(\mathbf a).$$
Если шаблон показателей числа \(m\) снова равен \(\mathbf a\), то само число \(m\) имеет этот шаблон, а значит
$$g(m)=G(\mathbf a)=m.$$
Обратно, любой фиксированный пункт \(n\) порождает свой шаблон и проходит именно эту проверку. Поэтому поиск можно вести целиком в пространстве шаблонов: сгенерировать \(\mathbf a\), вычислить \(G(\mathbf a)\), разложить получившееся число на множители и сравнить шаблоны.
Среди всех чисел с шаблоном \(\mathbf a\) наименьшим будет то, где наибольший показатель присвоен наименьшему простому, второй по величине - второму простому и так далее:
$$n_{\min}(\mathbf a)=2^{a_1}3^{a_2}5^{a_3}\cdots.$$
Если
$$n_{\min}(\mathbf a)\gt 10^{16},$$
то ни одно число с таким шаблоном не может быть ответом, и вся ветвь рекурсии отбрасывается сразу. Есть и второе быстрое условие:
$$G(\mathbf a)\lt n_{\min}(\mathbf a)\implies \text{шаблон}(G(\mathbf a))\ne \mathbf a.$$
В таком случае даже факторизация кандидата не нужна.
Возьмем \(\mathbf a=(4,1)\), то есть шаблон числа \(48=2^4\cdot 3\). Здесь \(\omega=5\). Сначала получаем
$$\begin{aligned} T_1&=\binom{4}{0}\binom{1}{0}=1,\\ T_2&=\binom{5}{1}\binom{2}{1}=10,\\ T_3&=\binom{6}{2}\binom{3}{2}=45,\\ T_4&=\binom{7}{3}\binom{4}{3}=140,\\ T_5&=\binom{8}{4}\binom{5}{4}=350. \end{aligned}$$
После применения включения-исключения имеем
$$\begin{aligned} F_1&=1,\\ F_2&=-2T_1+T_2=8,\\ F_3&=3T_1-3T_2+T_3=18,\\ F_4&=-4T_1+6T_2-4T_3+T_4=16,\\ F_5&=5T_1-10T_2+10T_3-5T_4+T_5=5. \end{aligned}$$
Значит,
$$G(4,1)=1+8+18+16+5=48.$$
Поскольку \(48\) действительно имеет шаблон \((4,1)\), это настоящий фиксированный пункт, который входит в искомую сумму. Для сравнения, \((2,1)\) дает \(G(2,1)=8\), а шаблон числа \(8\) равен \((3)\), поэтому \(12\) фиксированным пунктом не является.
Реализации на C++, Python и Java рекурсивно перечисляют шаблоны показателей в невозрастающем порядке. Поиск начинается с безопасной верхней границы для первого показателя и затем расширяет шаблон по одному простому множителю за шаг. На каждом уровне поддерживается минимально возможное число с текущим шаблоном; если эта нижняя граница уже превысила \(10^{16}\), ветвь немедленно обрывается.
Для каждого выжившего шаблона реализация точно строит нужные биномиальные коэффициенты, вычисляет произведения \(T_k\), затем применяет формулу включения-исключения для получения \(F_k\), а после этого суммирует их в \(G(\mathbf a)\). Если промежуточная сумма уже превысила глобальный предел, дальнейшие вычисления прекращаются, потому что более крупные значения не могут быть решениями.
Факторизация выполняется только для кандидатов \(m\), удовлетворяющих \(m\ge n_{\min}(\mathbf a)\). Из числа \(m\) извлекаются показатели простых множителей, сортируются по невозрастанию и сравниваются с исходным шаблоном. Уже проанализированные кандидаты кэшируются, найденные фиксированные пункты дедуплицируются, а итоговая сумма накапливается в широком целочисленном типе.
Пусть \(P\) обозначает количество шаблонов показателей, у которых минимальный представитель не превосходит \(10^{16}\). Для шаблона \(\mathbf a\) веса \(\omega=\sum a_i\) и длины \(r\) вычисление \(G(\mathbf a)\) требует \(O(\omega^2+r\omega)\) арифметических операций: \(O(\omega^2)\) уходит на таблицу биномиальных коэффициентов и слои включения-исключения, а \(O(r\omega)\) - на произведения \(T_k\).
Рекурсивный перебор остается практичным, потому что отсечение по нижней границе делает \(P\) небольшим. Этап факторизации использует детерминированный тест простоты Миллера-Рабина для 64-битных чисел вместе с разложением Полларда rho. Для этой части нет простой закрытой формулы на худший случай, но для чисел ниже \(10^{16}\) она очень быстра на практике и применяется лишь к небольшому числу кандидатов. Память расходуется умеренно: \(O(\omega^2)\) при обработке одного шаблона, плюс кэш разобранных кандидатов и множество найденных фиксированных пунктов.
سلسلة gozinta للعدد \(n\) هي متتالية من الشكل
$$1=d_0 \lt d_1 \lt \cdots \lt d_k=n$$
بحيث يقسم كل حد الحد الذي يليه. ولتكن \(g(n)\) عدد هذه السلاسل. المطلوب هو جمع كل الأعداد \(n\le 10^{16}\) التي تحقق
$$g(n)=n.$$
الفكرة الحاسمة هي أن \(g(n)\) لا تعتمد على القيم الفعلية للأعداد الأولية في تحليل \(n\)، بل تعتمد فقط على نمط أسسها بعد ترتيبها تنازليًا. بهذه الملاحظة تتحول المسألة من بحث هائل على جميع الأعداد حتى \(10^{16}\) إلى بحث أصغر بكثير في فضاء أنماط الأسس.
نكتب التحليل الأولي بالشكل ذي الأسس غير المتزايدة:
$$n=\prod_{i=1}^{r} p_i^{a_i},\qquad a_1\ge a_2\ge \cdots \ge a_r \gt 0.$$
ولنرمز إلى النمط بالمتجه \(\mathbf a=(a_1,\dots,a_r)\)، وإلى مجموع الأسس بالرمز
$$\omega=a_1+\cdots+a_r.$$
كل قاسم للعدد \(n\) يكتب على الصورة
$$\prod_{i=1}^{r} p_i^{b_i},\qquad 0\le b_i\le a_i.$$
إذن سلسلة gozinta تكافئ تمامًا متتالية من متجهات الأسس، تبدأ من \((0,\dots,0)\) وتنتهي عند \((a_1,\dots,a_r)\)، مع عدم تناقص كل إحداثي على امتداد السلسلة.
وبما أن هذا الوصف لا يستخدم إلا الأسس نفسها، فإن عدد السلاسل يتحدد بالكامل بالنمط \(\mathbf a\). لذلك يمكن كتابة
$$g(n)=G(\mathbf a).$$
هذه هي الخطوة الجوهرية التي تبني عليها الشيفرات الثلاث الحل كله.
لنفرض أن السلسلة تحتوي بالضبط على \(k\) خطوات غير تافهة:
$$1=d_0 \lt d_1 \lt \cdots \lt d_k=n.$$
نعرّف النسب بين الحدود المتتالية بالصيغة \(q_t=d_t/d_{t-1}\). عندئذ يكون كل \(q_t \gt 1\) ولدينا
$$q_1q_2\cdots q_k=n.$$
بالنسبة إلى عامل أولي واحد \(p_i^{a_i}\)، نوزع الأس \(a_i\) على هذه النسب \(k\):
$$e_{i,1}+e_{i,2}+\cdots+e_{i,k}=a_i,\qquad e_{i,t}\ge 0.$$
عدد التركيبات الضعيفة للعدد \(a_i\) إلى \(k\) أجزاء هو
$$\binom{a_i+k-1}{k-1}.$$
وبسبب استقلال العوامل الأولية نحصل على
$$T_k(\mathbf a)=\prod_{i=1}^{r}\binom{a_i+k-1}{k-1}.$$
لكن \(T_k(\mathbf a)\) ما زالت تسمح ببعض النسب المساوية لـ \(1\)، أي بخطوات فارغة، ولذلك فهي تعد أكثر من السلاسل الحقيقية.
السلسلة الصحيحة يجب أن تكون كل خطوة فيها غير تافهة. إذا كانت \(j\) مواضع فقط من أصل \(k\) مواضع نشطة، فإن توزيعات الأسس تُحصى بواسطة \(T_j(\mathbf a)\)، كما أن اختيار المواضع النشطة يتم بعدد
$$\binom{k}{j}$$
من الطرق. بتطبيق مبدأ الشمول والاستبعاد على مواضع الفراغ نحصل على عدد السلاسل التي تحتوي بالضبط على \(k\) خطوات غير تافهة:
$$F_k(\mathbf a)=\sum_{j=1}^{k}(-1)^{k-j}\binom{k}{j}T_j(\mathbf a).$$
وهذا هو المجموع المتناوب المركزي الذي تنفذه تطبيقات C++ وPython وJava.
كل خطوة غير تافهة تزيد مجموع الأسس بمقدار لا يقل عن \(1\)، لذلك لا يمكن أن يتجاوز عدد الخطوات القيمة
$$\omega=a_1+\cdots+a_r.$$
ومن ثم يكون العدد الكلي لسلاسل gozinta للنمط \(\mathbf a\) هو
$$G(\mathbf a)=\sum_{k=1}^{\omega}F_k(\mathbf a).$$
بهذا نكون قد حوّلنا مسألة سلاسل القواسم إلى مجاميع منتهية دقيقة من معاملات ثنائية الحدود.
لتكن
$$m=G(\mathbf a).$$
إذا كان نمط أسس \(m\) هو نفسه \(\mathbf a\)، فهذا يعني أن \(m\) يحمل هذا النمط أصلًا، وبالتالي
$$g(m)=G(\mathbf a)=m.$$
وبالعكس، كل نقطة ثابتة \(n\) تولد نمطها الخاص وستجتاز هذا الاختبار. لذلك يمكن إجراء البحث كله في فضاء الأنماط: نولد النمط \(\mathbf a\)، نحسب \(G(\mathbf a)\)، نحلل العدد الناتج إلى عوامله الأولية، ثم نقارن نمط الأسس بالنمط الأصلي.
أصغر عدد يملك النمط \(\mathbf a\) ينتج من إسناد أكبر أس إلى أصغر عدد أولي، وثاني أكبر أس إلى ثاني أصغر عدد أولي، وهكذا:
$$n_{\min}(\mathbf a)=2^{a_1}3^{a_2}5^{a_3}\cdots.$$
إذا كان
$$n_{\min}(\mathbf a)\gt 10^{16},$$
فلا يوجد أي عدد بهذا النمط يمكن أن يساهم في الجواب، ولذلك تُقص هذه الشعبة مباشرة. وهناك فحص سريع ثانٍ:
$$G(\mathbf a)\lt n_{\min}(\mathbf a)\implies \text{pattern}(G(\mathbf a))\ne \mathbf a.$$
أي إذا كانت القيمة المحسوبة أصغر من أصغر عدد ممكن لهذا النمط، فلا داعي أصلًا لتحليلها إلى عوامل أولية.
خذ النمط \(\mathbf a=(4,1)\)، وهو نمط العدد \(48=2^4\cdot 3\). هنا \(\omega=5\). نحسب أولًا
$$\begin{aligned} T_1&=\binom{4}{0}\binom{1}{0}=1,\\ T_2&=\binom{5}{1}\binom{2}{1}=10,\\ T_3&=\binom{6}{2}\binom{3}{2}=45,\\ T_4&=\binom{7}{3}\binom{4}{3}=140,\\ T_5&=\binom{8}{4}\binom{5}{4}=350. \end{aligned}$$
ثم بعد تطبيق الشمول والاستبعاد نحصل على
$$\begin{aligned} F_1&=1,\\ F_2&=-2T_1+T_2=8,\\ F_3&=3T_1-3T_2+T_3=18,\\ F_4&=-4T_1+6T_2-4T_3+T_4=16,\\ F_5&=5T_1-10T_2+10T_3-5T_4+T_5=5. \end{aligned}$$
إذن
$$G(4,1)=1+8+18+16+5=48.$$
وبما أن \(48\) نفسه يملك النمط \((4,1)\)، فهو نقطة ثابتة حقيقية ويجب أن يدخل في المجموع النهائي. أما النمط \((2,1)\) فيعطي \(G(2,1)=8\)، ونمط \(8\) هو \((3)\)، لذا فإن \(12\) ليست نقطة ثابتة.
تقوم تطبيقات C++ وPython وJava بتوليد أنماط الأسس بشكل عودي بترتيب غير متزايد. يبدأ البحث بحد أعلى آمن للأس الأول، ثم يمد النمط باستخدام قوى الأعداد الأولية المتتالية. وفي كل خطوة تُحدَّث أصغر قيمة ممكنة تمثل هذا النمط؛ فإذا تجاوز هذا الحد \(10^{16}\)، أُغلقت الشعبة فورًا.
لكل نمط ينجو من هذا القطع، تبني الشيفرة معاملات ثنائية الحدود اللازمة بدقة، وتحسب جداءات \(T_k\)، ثم تطبق صيغة الشمول والاستبعاد للحصول على \(F_k\) ومن ثم \(G(\mathbf a)\). وإذا تجاوز المجموع الجزئي الحد الكلي تتوقف العملية عند ذلك، لأن القيم الأكبر لا يمكن أن تحقق شرط المسألة.
لا تُحلَّل إلى عوامل أولية إلا القيم \(m\) التي تحقق \(m\ge n_{\min}(\mathbf a)\). بعد ذلك تُستخرج أسس العوامل الأولية وتُرتب ترتيبًا تنازليًا، ثم تُقارن بالنمط الأصلي. كما تُخزَّن القيم التي سبق تحليلها، وتُزال التكرارات من نقاط التثبيت الصحيحة، ويُجمع الناتج النهائي باستخدام حساب صحيح واسع السعة.
لنرمز بـ \(P\) إلى عدد أنماط الأسس التي لا يتجاوز أصغر ممثل لها \(10^{16}\). بالنسبة إلى نمط \(\mathbf a\) ذي الوزن \(\omega=\sum a_i\) والطول \(r\)، فإن حساب \(G(\mathbf a)\) يتطلب \(O(\omega^2+r\omega)\) عملية حسابية: \(O(\omega^2)\) لبناء جدول معاملات ثنائية الحدود ولمراحل الشمول والاستبعاد، و\(O(r\omega)\) لبناء جداءات \(T_k\).
يبقى التعداد العودي عمليًا لأن القطع بواسطة الحد الأدنى يجعل \(P\) صغيرًا. أما مرحلة التحليل إلى عوامل أولية فتستخدم اختبار Miller-Rabin الحتمي للأعداد الصحيحة ذات 64 بت مع طريقة Pollard-Rho للتفكيك. لا توجد لها صيغة مغلقة بسيطة لأسوأ حالة، لكنها سريعة جدًا عمليًا للأعداد الأصغر من \(10^{16}\)، كما أنها لا تعمل إلا على مجموعة صغيرة نسبيًا من المرشحين الذين نجوا من المراحل السابقة. واستهلاك الذاكرة معتدل: \(O(\omega^2)\) أثناء معالجة نمط واحد، إضافة إلى مخزن القيم المحللة ومجموعة نقاط التثبيت المكتشفة.