Detailed mathematical approach

Problem Summary

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.

Mathematical Approach

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\).

Step 1: Replace the Integer by Its Exponent Pattern

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.

Step 2: Count Chains with a Fixed Number of Steps

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.

Step 3: Use Inclusion-Exclusion to Remove Empty Steps

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.

Step 4: Sum over All Admissible Lengths

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.

Step 5: Convert \(g(n)=n\) into a Pattern Fixed-Point Test

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.

Step 6: Prune the Search with the Smallest Representative

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.

Worked Example

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.

How the Code Works

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.

Complexity Analysis

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.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=548
2. Stars and bars / weak compositions: Wikipedia — Stars and bars
3. Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
4. Miller-Rabin primality test: Wikipedia — Miller-Rabin primality test
5. Pollard-Rho factorization: Wikipedia — Pollard's rho algorithm