Detailed mathematical approach

Problem Summary

Let \(r\) be the irrational base used in the problem. For each positive integer \(n\), we consider its greedy representation as a sum of distinct powers of \(r\), and we write \(\ell(n)\) for the number of selected powers. The quantity to compute is

$S(m)=\sum_{j=1}^{m}\ell(j^2),\qquad m=5{,}000{,}000.$

The challenge is that the base is irrational, so the implementation must make millions of precise comparisons and subtractions without letting floating-point error corrupt the greedy choices.

Mathematical Approach

The implementations rely on one algebraic identity for the base and one numerical device for comparing powers quickly.

Step 1: Identify the Irrational Base

The base is

$r=\frac{1+\sqrt[3]{\frac{29-3\sqrt{93}}{2}}+\sqrt[3]{\frac{29+3\sqrt{93}}{2}}}{3}\approx 1.465571231876768.$

This is the unique real root of

$x^3-x^2-1=0,$

so it satisfies

$r^3=r^2+1.$

Multiplying by \(r^{e-2}\) gives the recurrence

$r^{e+1}=r^e+r^{e-2}.$

That recurrence explains the combinatorial shape of the greedy expansion.

Step 2: Why the Exponents Are Spaced by at Least Three

Suppose the greedy process has chosen the largest admissible power \(r^e\le n\). If the residual is \(R=n-r^e\), then \(n<r^{e+1}\), hence

$0\le R<r^{e+1}-r^e=r^{e-2}.$

Therefore the next selected power cannot be \(r^{e-1}\) or \(r^{e-2}\); its exponent must be at most \(e-3\). Repeating the same argument at every step yields the greedy form used by the implementations:

$n=\sum_{t=1}^{\ell(n)} r^{e_t},\qquad e_1>e_2>\cdots,\qquad e_{t+1}\le e_t-3.$

So after one term is chosen, the next search can start three exponents lower and continue downward only if the residual is still smaller than that first candidate.

Step 3: Worked Example

Take \(n=4\). Since

$r^3\approx 3.147899<4<r^4\approx 4.613470,$

the first greedy term is \(r^3\). The remaining amount is about \(0.852101\), so the next admissible term is \(r^{-1}\approx 0.682328\). Continuing in the same way gives

$4=r^3+r^{-1}+r^{-5}+r^{-10}.$

Numerically,

$3.147899+0.682328+0.147899+0.021874\approx 4.$

The decimal values only illustrate the greedy choices; the equality itself is exact because all terms are powers of the algebraic number \(r\). Hence \(\ell(4)=4\).

Step 4: Encode Powers in a KaTeX-Friendly Fixed-Point Form

Directly comparing irrational numbers inside the inner loop would be both slow and fragile. Instead, the implementations precompute each power in the form

$r^e=a_e+\frac{b_e}{M}+\frac{c_e}{M^2}+O(M^{-3}),\qquad M=10^{17},$

where \(a_e\), \(b_e\), and \(c_e\) are integers obtained from a high-precision evaluation of \(r^e\). Residuals are stored in the same three-level format. To decide whether a power fits into the current residual, the implementation only needs a lexicographic comparison of the triples

$\left(a_e,b_e,c_e\right),\qquad \left(R_1,R_2,R_3\right).$

After each subtraction, borrow propagation normalizes the three components exactly as in arithmetic with base \(M\).

Step 5: Turn Individual Lengths into the Final Sum

For every square \(j^2\), the algorithm finds the largest precomputed power below it, repeatedly subtracts the largest admissible lower power, and counts how many terms were used. That count is \(\ell(j^2)\), so the final answer is simply

$S(m)=\sum_{j=1}^{m}\ell(j^2).$

The implementations verify two checkpoints before the full run:

$S(10)=61,\qquad S(1000)=19403.$

These values confirm that the greedy representation and the fixed-point comparisons are aligned.

How the Code Works

The C++ and Java implementations first evaluate the algebraic base with high decimal precision and build a table of powers over a finite exponent window wide enough for the target computation. Each stored power is split into three integer layers so that the inner greedy loop can avoid arbitrary-precision arithmetic.

For each square, the implementation locates the leading power, initializes the residual in the same three-layer format, then repeatedly searches downward for the next admissible term, subtracts it with borrow handling, and increments the representation length. Because the squares increase with \(j\), the search for the leading exponent only moves forward through the table. The Python implementation delegates to the same compiled algorithm, so all three language versions use the same greedy decomposition and the same checkpoints.

Complexity Analysis

Let \(W=E_{\max}-E_{\min}+1\) denote the size of the precomputed exponent window. Building the table costs \(O(W)\) high-precision steps and \(O(W)\) memory. Over the whole sweep \(j=1,\dots,m\), the forward search for the leading exponent is monotone, so it contributes \(O(W+m)\) total work. The dominant part is the greedy subtraction itself, which is proportional to the total number of selected terms:

$O\left(W+m+\sum_{j=1}^{m}\ell(j^2)\right).$

Since \(W\) is fixed in the implementations, memory usage is effectively constant with respect to \(m\), and the running time is essentially linear in the total output length of the greedy expansions.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=558
2. Beta expansion: Wikipedia — Beta expansion
3. Algebraic number: Wikipedia — Algebraic number
4. Fixed-point arithmetic: Wikipedia — Fixed-point arithmetic
5. Greedy algorithm: Wikipedia — Greedy algorithm