Problem 356: Largest Roots of Cubic Polynomials

Project Euler Problem 356 Solution

EulerSolve provides an optimized solution for Project Euler Problem 356, Largest Roots of Cubic Polynomials, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each $n\in\{1,\dots,30\}$, let $a_n$ be the largest real root of $$x^3-2^n x^2+n=0.$$ The task is to compute $$\sum_{n=1}^{30}\left\lfloor a_n^{987654321}\right\rfloor \pmod{10^8}.$$ A direct floating-point evaluation of such enormous powers would be fragile and unnecessary. The local C++, Python, and Java solutions instead turn the problem into an exact linear recurrence modulo $10^8$. Mathematical Approach Step 1: Locate the Three Real Roots Write $c=2^n$ and $f(x)=x^3-cx^2+n$. The implementations rely on the fact that the cubic has three real roots, with one large positive root and two small roots near the origin. Indeed, $$f(-1)=n-1-c<0,\qquad f(0)=n>0,$$ so there is a root $\gamma\in(-1,0)$. Also, $$f(1)=n+1-c\le 0,\qquad f(0)=n>0,$$ so there is a root $\beta\in(0,1]$, with $\beta=1$ only when $n=1$. Finally, $$f(c)=n>0,\qquad f(c-1)=n-(c-1)^2\le 0,$$ so the largest root $\alpha=a_n$ lies in $(c-1,c)$....

Detailed mathematical approach

Problem Summary

For an integer $n\ge1$ let $g(x)=x^3-2^nx^2+n$, and let $a_n$ be the largest real root of $g$. The statement gives $a_2=3.86619826\cdots$ and asks for the last eight digits of

$$\sum_{i=1}^{30}\lfloor a_i^{987654321}\rfloor,$$

where $\lfloor x\rfloor$ is the largest integer not exceeding $x$. ^[1]

Ordinary floating-point arithmetic cannot settle these floors. The largest root $a_{30}$ sits just below $2^{30}$, so the single term $i=30$ is a whole number with about $8.9$ billion decimal digits; the sum itself is never written down, only its residue modulo $10^8$.

The whole solution in one sentence: adding the $m$-th powers of all three roots of $g$ always gives a whole number, that whole number is exactly one more than the floor we want whenever $m$ is odd, and a three-term recurrence yields its last eight digits by repeated squaring.

Mathematical Approach

1. Three real roots, and where they lie

Write $c=2^n$ throughout. We first need the elementary inequality $2^n\ge n+1$, with equality only at $n=1$: it holds at $n=1$, and if $2^n\ge n+1$ then $2^{n+1}\ge2n+2\ge n+2$. A cubic has at most three real roots, and five values of $g$, all exact integers, locate every one of them.

$g(-1)=n-1-c\le-2$, since $c\ge n+1$.
$g(0)=n\gt0$.
$g(1)=n+1-c\le0$, with equality only at $n=1$. For $n=1,2,3,4$ the values are $0,-1,-4,-11$.
$g(c-1)=n-(c-1)^2\le0$: the inequality $c-1\ge n$ gives $(c-1)^2\ge n^2\ge n$. For $n=1,2,3,4$ the values are $0,-7,-46,-221$, so the inequality is strict from $n=2$ on.
$g(c)=n\gt0$.

A polynomial is continuous, so a sign change on an interval forces a root inside it; this is Bolzano's form of the intermediate value theorem. ^[2] For $n\ge2$ the three intervals $(-1,0)$, $(0,1)$ and $(c-1,c)$ are pairwise disjoint, because then $c-1\ge3$. Each carries a root, three distinct roots exhaust a cubic, and therefore all three roots are real and

$$\begin{aligned} \gamma&\in(-1,0),\\ \beta&\in(0,1],\\ \alpha&\in(c-1,c), \end{aligned}$$

with $\alpha=a_n$ the largest. The case $n=1$ is the one place where $g(1)=0$, and there the cubic factors:

$$x^3-2x^2+1=(x-1)(x^2-x-1),$$

so the roots are exactly $1$ and $\tfrac{1\pm\sqrt5}{2}$, that is $\gamma=-1/\varphi\approx-0.6180339887$, $\beta=1$ and $\alpha=\varphi\approx1.6180339887$, where $\varphi$ denotes the golden ratio: the same picture, with $\beta=1$ instead of $\beta\lt1$. Two further numerical examples: for $n=2$ the roots are $\gamma\approx-0.6554423815$, $\beta\approx0.7892441190$ and $\alpha\approx3.8661982625$, matching the $a_2$ of the statement; for $n=30$ they are $\gamma\approx-0.0001671517$, $\beta\approx0.0001671517$ and $\alpha\approx2^{30}-2.60\cdot10^{-17}$.

2. What Vieta's formulas give

Reading $g(x)=x^3-cx^2+0\cdot x+n$ coefficient by coefficient, Vieta's formulas ^[3] state that the three roots satisfy

$$\alpha+\beta+\gamma=c,\quad\alpha\beta+\alpha\gamma+\beta\gamma=0,$$

$$\alpha\beta\gamma=-n.$$

Three consequences are used below. First, $\beta+\gamma=c-\alpha$, and since $\alpha^3-c\alpha^2+n=0$ can be rewritten as $\alpha^2(c-\alpha)=n$,

$$\beta+\gamma=c-\alpha=\frac{n}{\alpha^2}\gt0 .$$

For $n=30$ this is about $2.60\cdot10^{-17}$, the gap between $\alpha$ and $2^{30}$ quoted above. Second, $\beta\gamma=-n/\alpha\lt0$, so the two small roots have opposite signs; for $n=2$ it is $\beta\gamma=-2/\alpha\approx-0.5173040450$. Third, because $\beta+\gamma$ is tiny while $\beta\gamma\approx-n/c$, both small roots are close to $\pm\sqrt{n/2^n}$ in size: for $n=30$ that estimate is $\sqrt{30}/2^{15}\approx0.0001671517$, which agrees with $\beta$ to the ten digits shown.

3. The power sums are integers

Define the power sums

$$S_m=\alpha^m+\beta^m+\gamma^m\qquad(m\ge0).$$

Every root $r$ of $g$ satisfies $r^3=cr^2-n$. Multiplying by $r^{m-3}$ gives $r^m=cr^{m-1}-nr^{m-3}$ for $m\ge3$, and adding the three resulting equations yields

$$S_m=cS_{m-1}-nS_{m-3}\qquad(m\ge3).$$

The three starting values come from Section 2: $S_0=3$ counts the roots, $S_1=c$, and

$$\begin{aligned} S_2&=(\alpha+\beta+\gamma)^2\\ &\quad-2(\alpha\beta+\alpha\gamma+\beta\gamma)\\ &=c^2 . \end{aligned}$$

Since $c=2^n$ and $n$ are integers and the recurrence has integer coefficients, induction on $m$ shows that every $S_m$ is an integer — even though for $n\ge2$ no root is rational at all, since a rational root of a monic polynomial with integer coefficients must be an integer and the three intervals above contain none. These are the power sums of Newton's identities, which for three variables and $m\gt3$ read $S_m=e_1S_{m-1}-e_2S_{m-2}+e_3S_{m-3}$ with $e_1=c$, $e_2=0$ and $e_3=-n$: the recurrence above. ^[4] The elementary derivation from $r^3=cr^2-n$ also covers $m=3$.

4. The key identity for odd exponents

Claim. For every $n\ge1$ and every odd $m\ge1$,

$$\lfloor a_n^m\rfloor=S_m-1 .$$

$$\rho_m=\beta^m-|\gamma|^m .$$

From $\beta\gt|\gamma|\gt0$ and the strict increase of $t\mapsto t^m$ on positive $t$ it follows that $\rho_m\gt0$; from $|\gamma|^m\gt0$ and $\beta^m\le1$ it follows that $\rho_m\lt1$. So $0\lt\rho_m\lt1$, and because $S_m$ is an integer,

$$\lfloor\alpha^m\rfloor=\lfloor S_m-\rho_m\rfloor=S_m-1 .$$

The exponent $987654321$ is odd, so the identity applies to all thirty terms of the required sum. Note also that $\alpha\gt1$ gives $S_m\ge2$ for $m\ge1$, so each $S_m-1$ is a positive integer. Nothing about the size of $n$ was used, so the identity holds for every $n\ge1$.

5. Worked example: $n=1$ and the Lucas numbers

For $n=1$ the recurrence is $S_m=2S_{m-1}-S_{m-3}$ with $S_0=3$, $S_1=2$, $S_2=4$. Running it by hand gives, for $m=0,\dots,9$,

$$3,\;2,\;4,\;5,\;8,\;12,\;19,\;30,\;48,\;77 .$$

These are the Lucas numbers $L_m=2,1,3,4,7,11,18,29,47,76$, each increased by one. ^[5] The reason is Section 1: the roots are $1$, $\varphi$ and $-1/\varphi$, so

$$\begin{aligned} S_m&=1^m+\varphi^m+(-1/\varphi)^m\\ &=1+L_m. \end{aligned}$$

The claim of Section 4 therefore says $\lfloor\varphi^m\rfloor=L_m$ for odd $m$, which a reader can check by hand: $\varphi\approx1.6180339887$ gives $\lfloor\varphi\rfloor=1=L_1$, $\varphi^3\approx4.2360679775$ gives $4=L_3$, and $\varphi^5\approx11.0901699437$ gives $11=L_5$. The same source records the general formula $L_m=\lfloor\varphi^m\rfloor+\tfrac{(-1)^m+1}2$, whose correction term is $0$ for odd $m$ and $1$ for even $m$ — a first warning that the parity in Section 4 matters.

6. Worked example: $n=2$, and why the exponent must be odd

For $n=2$ we have $c=4$, so $S_0=3$, $S_1=4$, $S_2=16$, and then

$$\begin{aligned} S_3&=4\cdot16-2\cdot3=58,\\ S_4&=4\cdot58-2\cdot4=224,\\ S_5&=4\cdot224-2\cdot16=864,\\ S_6&=4\cdot864-2\cdot58=3340 . \end{aligned}$$

Carrying $\alpha$ to high precision gives $\alpha^3\approx57.7899560201$ and $\alpha^5\approx863.8147322117$, so $\lfloor\alpha^3\rfloor=57=S_3-1$ and $\lfloor\alpha^5\rfloor=863=S_5-1$, as claimed.

The even exponents behave differently. Here $\alpha^2\approx14.9474890050$, so $\lfloor\alpha^2\rfloor=14$, while $S_2-1=15$. The reason is visible in the proof: for even $m$ both $\beta^m$ and $\gamma^m$ are positive, so $\rho_2=\beta^2+\gamma^2\approx1.0525109950$ exceeds $1$ and the floor drops by $2$ instead of $1$. For $n=1$ the same happens at every even $m$, since $\beta=1$ makes $\rho_m=1+|\gamma|^m\gt1$. The identity of Section 4 is therefore about odd exponents only.

7. How close the powers really are

An algebraic integer is a root of a monic polynomial with integer coefficients, so $\alpha$ is one; its conjugates are the other roots of its minimal polynomial, the monic rational polynomial of least degree with $\alpha$ as a root. That minimal polynomial divides $g$, so every conjugate of $\alpha$ besides $\alpha$ itself is $\beta$ or $\gamma$. For $n\ge2$, Section 1 gives $0\lt\beta\lt1$ and $-1\lt\gamma\lt0$ while $\alpha\gt1$, so $\alpha$ is a Pisot–Vijayaraghavan number: a real algebraic integer greater than $1$ whose other conjugates all have absolute value below $1$. ^[6] For $n=1$ the factor $x-1$ splits off and $\alpha=\varphi$ has the single conjugate $-1/\varphi$, so the golden ratio is one too. The defining property of these numbers is exactly the mechanism of Section 4: the sum of the $m$-th powers of such a number and its conjugates is an integer, and the conjugate powers shrink geometrically, so $\alpha^m$ approaches an integer at an exponential rate.

Here that approach is extreme. For $n=1$ and odd $m$ we have $\rho_m=1-\varphi^{-m}$, so $\alpha^m$ lies above the integer $S_m-1$ by $\varphi^{-m}$, a distance whose base-ten logarithm at $m=987654321$ is about $-2.1\cdot10^{8}$. For $n=30$, where $\beta\approx1.67\cdot10^{-4}$, the number $\alpha^m$ lies below the integer $S_m$ by $\rho_m=\beta^m-|\gamma|^m$, whose base-ten logarithm is about $-3.7\cdot10^{9}$. No hardware floating-point format comes anywhere near either gap, and the strict inequalities of Section 4 settle the floor without approximating these powers at all.

8. The matrix form and binary exponentiation

The recurrence is linear of order $3$, so three consecutive power sums determine the next triple. Collect them newest first in the column vector

$$v_m=\begin{pmatrix}S_m\\S_{m-1}\\S_{m-2}\end{pmatrix},$$

so that the recurrence reads $v_m=Mv_{m-1}$ with

$$M=\begin{pmatrix}c&0&-n\\1&0&0\\0&1&0\end{pmatrix}.$$

Its first row is the recurrence itself; the other two shift the window. Iterating from the starting triple gives, for $m\ge2$,

$$v_m=M^{m-2}\begin{pmatrix}S_2\\S_1\\S_0\end{pmatrix},$$

so $S_m$ is the first entry, that is the first row of $M^{m-2}$ applied to $(c^2,c,3)$. This $M$ is the companion matrix of $g$ written for a window ordered newest first ^[7]: expanding $\det(xI-M)$ along its first row gives $(x-c)x^2+n=g(x)$, so by the Cayley–Hamilton theorem $M^3=cM^2-nI$, which is the recurrence once more.

A small instance can be checked by hand. For $n=2$, squaring $M$ gives the rows $(16,-2,-8)$, $(4,0,-2)$, $(1,0,0)$, and one further product gives $M^3$ with rows $(62,-8,-32)$, $(16,-2,-8)$, $(4,0,-2)$. Applying $M^3$ to $(16,4,3)$ yields $992-32-96=864$, then $224$ and $58$ — exactly $(S_5,S_4,S_3)$ of Section 6.

The power $M^{m-2}$ is built by repeated squaring. ^[8] Start an accumulator at the identity matrix and a running power at $M$, and read the bits of the exponent from the least significant upwards: at bit $j$ the running power is $M^{2^j}$, so multiply it into the accumulator when that bit is $1$, then square it for the next bit. The selected powers $2^j$ add up to $m-2$, so the accumulator ends as $M^{m-2}$. In the hand example $m-2=3=11_2$ selects $M^1$ and $M^2$, whose product is $M^3$. For the required exponent, $m-2=987654319$ has $30$ binary digits of which $19$ are ones, and the published loops square on every pass, the last one included, and start from the identity: $30$ squarings and $19$ accumulator products, that is $49$ matrix products for one value of $n$ and $1470$ for all thirty. Two of those are avoidable, the multiplication by the identity and the final unused square, so a tighter loop would use $47$.

9. Working modulo $10^8$

Only additions, subtractions and multiplications occur, and congruence modulo a fixed $10^8$ is compatible with all three operations, so every intermediate value may be reduced immediately. ^[9] No modular inverse is ever needed, so no division has to be justified: since $10^8=2^8\cdot5^8$, a residue is invertible modulo $10^8$ exactly when it is coprime to $10$, which most residues are not. Three details complete the method.

The entry $-n$ of $M$ is stored as its non-negative representative $10^8-n$, so all nine entries stay in $\{0,\dots,10^8-1\}$.
Each term is $S_m-1$, a positive integer by Section 4. If $r=S_m\bmod10^8$, its residue is $(r+10^8-1)\bmod10^8$. Subtracting before or after reducing gives the same residue; adding the modulus first only keeps the intermediate value non-negative, which matters when $r=0$ in unsigned arithmetic.
The thirty residues are added and reduced one at a time. The exact sum is a positive integer, so its last eight digits are this residue, padded with leading zeros if it is below $10^7$.

Finally, only $c=2^n\bmod10^8$ is ever needed, never $2^n$ itself.

How the Code Works

For each $n$, all three programs reduce $c=2^n$ modulo $10^8$, raise the matrix of Section 8 to the power $m-2$, combine its first row with $(S_2,S_1,S_0)$, subtract $1$ modularly and accumulate.

In the C++ program ^[11] the defaults are the modulus $10^8$, the exponent $987654321$ and the range $1\le n\le30$. parse_arguments accepts --n-min=, --n-max=, --exp= and --skip-checks, and it rejects an even or zero exponent — the hypothesis of Section 4 is enforced at the command line rather than assumed. pow_mod_u64 computes $c$ by repeated squaring. Matrix3 holds nine 64-bit residues; multiply_matrix3 accumulates the three products of a row and a column in an unsigned __int128 before reducing, and power_matrix3 is the loop of Section 8, starting from make_identity_matrix3. compute_newton_sum_mod returns $S_0$, $S_1$ or $S_2$ directly for the exponents $0$, $1$, $2$, and otherwise builds the transition matrix, writing its top-right entry as (mod - (n % mod)) % mod, which reduces $n$ first and equals $10^8-n$ for every $n$ this problem needs, and then forms the first-row combination. compute_term_mod adds $10^8-1$ and reduces, giving the residue of $S_m-1$. solve_range walks over the requested values of $n$ and keeps the running total in range with a conditional subtraction instead of a remainder operation; its counter is a 32-bit unsigned integer, so the extreme upper bound $2^{32}-1$ that the parser still accepts would make that counter wrap round instead of ending the loop, a limit the required range $1\le n\le30$ never approaches. main prints the total through std::setfill('0') and std::setw(8), so the eight digits appear even with leading zeros.

run_validations, switched off by --skip-checks, contains the self-tests. largest_real_root brackets $\alpha$ in $[c-1,c]$ as in Section 1 and halves the interval $220$ times, and floor_power multiplies the result out for a small exponent. The first check requires $\lfloor a_2^2\rfloor=14$, the value implied by the $a_2$ of the statement. The second compares compute_newton_sum_mod with compute_newton_sum_mod_iterative, which applies the recurrence one step at a time, for every $n\le12$ and every exponent $k\le220$. The third uses compute_newton_sum_exact_small, an exact 128-bit iteration with no modulus, to confirm $\lfloor a_n^m\rfloor=S_m-1$ against the floating-point value for $n\le3$ and $m\in\{1,3,5,7,9,11\}$. A fourth check, when the execution environment allows it, evaluates the enlarged range $n=1,\dots,400$ twice with the selected exponent and requires both totals to agree; otherwise it reports itself as skipped.

The Python and Java versions carry the same arithmetic under names adapted to each language: pow_mod and powMod, multiply_matrix and multiplyMatrix, power_matrix and powerMatrix, compute_newton_sum_mod and computeNewtonSumMod. Python keeps a separate compute_term_mod, whereas the Java solve performs that subtraction itself; both use the fixed range and exponent of the statement, and neither version carries self-tests or argument parsing. Java holds the matrix entries, the powers and the running total in long variables and the index $n$ in an int: a product of two residues is below $10^{16}$, and the three products inside a matrix entry sum to less than $3\cdot10^{16}$, far below $2^{63}-1\approx9.22\cdot10^{18}$, so nothing overflows. ^[10] Python's integers are unbounded, so no bound has to be checked there. One difference remains: C++ and Java pad the result to eight digits, while Python returns str(total), which would show fewer characters if the residue were below $10^7$.

Complexity and Verification

For one value of $n$ the work is the reduction of $c=2^n$, at most five squarings and four products for $n\le30$, plus the $49$ matrix products of Section 8. Each $3\times3$ product performs $27$ integer multiplications and nine remainder operations, one per entry of the result, so the matrix powers alone cost $1323$ multiplications for one value of $n$ and $39690$ for all thirty; the remaining operations are a handful. With the modulus fixed, the cost per value of $n$ is $O(\log n+\log m)$ multiplications, that is $O(\log m)$ over the range of this problem. A matrix power keeps a constant number of $3\times3$ arrays, so the working storage is $O(1)$. All three programs finish in a few hundredths of a second, the C++ self-tests included. A one-step iteration of the recurrence would instead need about $3.0\cdot10^{10}$ steps in all, which is why compute_newton_sum_mod_iterative serves only the checks.

The following facts were recomputed for this article, with exact integers where possible and otherwise with $200$ or more decimal digits.

The five integer values of Section 1 have the stated signs for every $n$ from $1$ to $30$, with equality only at $n=1$.
The roots to $200$ digits satisfy $\beta\le1$, with equality exactly at $n=1$, as well as $|\gamma|\lt1$, $\beta\gt|\gamma|$ and $\beta+\gamma=n/\alpha^2$, for every $n\le30$.
$\lfloor a_n^m\rfloor=S_m-1$ and $0\lt\rho_m\lt1$ hold for every $n\le30$ and every odd $m\le61$, with the working precision raised with $m$ so that $\rho_m$ stays resolvable.
The matrix route reproduces the exact integer recurrence reduced modulo $10^8$ for every $n\le30$ and every $m\le200$.
The full-size total was recomputed by two further methods: multiplying polynomials modulo $x^3-cx^2+n$, which uses Cayley–Hamilton instead of matrices, and running the computation modulo $2^8$ and modulo $5^8$ separately and recombining the two residues with the Chinese remainder theorem. ^[9] Both agree with the matrix method, and the three published programs print the same eight digits.

One limitation deserves naming. The third self-test compares an exact integer with a repeated product of long double values, and the precision of that type depends on the platform: it need not be the IEEE binary64 format. ^[10] In the build used for this article its significand has $53$ bits, and with the bisection and repeated-product routines of the program the comparison first fails at $m=71$ for $n=1$, at $m=27$ for $n=2$ and at $m=15$ for $n=3$ — observations for that build, not portable thresholds. To settle a floor numerically one needs an error bound below the distance to the nearer neighbouring integer, namely $\min(\rho_m,1-\rho_m)$; for $n=1$ and odd $m$ the smaller distance is $1-\rho_m=\varphi^{-m}$, not $\rho_m$. The packaged checks stop at $m=11$, well inside the reliable range, and every required exponent is covered by the proof of Section 4 instead.

Footnotes and References

The references below supply the background facts used above; the root locations, the integrality of the power sums, the identity for odd exponents and the operation counts are derived explicitly in this article.

Project Euler 356 — Largest Roots of Cubic Polynomials. The official statement defines aₙ as the largest real root of g(x) = x³ − 2ⁿ·x² + n, gives a₂ = 3.86619826…, notes that ⌊a⌋ is the floor function, and asks for the last eight digits of the sum of ⌊aᵢ^987654321⌋ for i = 1 to 30.
Wikipedia — Intermediate value theorem. The article states Bolzano's corollary: a continuous function with values of opposite sign on an interval has a root there. Section 1 applies it to the five integer values of g.
Wikipedia — Vieta's formulas. For a monic cubic x³ + bx² + cx + d with roots r₁, r₂, r₃ the article gives r₁ + r₂ + r₃ = −b, r₁r₂ + r₁r₃ + r₂r₃ = c and r₁r₂r₃ = −d; it attributes the relations to François Viète (1540–1603). Section 2 uses the three identities with b = −2ⁿ, c = 0 and d = n.
Wikipedia — Newton's identities. The Girard–Newton formulae relate power sums to elementary symmetric polynomials; for k greater than the number of variables n they reduce to a recurrence in the n previous power sums with no extra term, which for three variables is the relation of Section 3.
OEIS A000032 — Lucas numbers. L(0) = 2, L(1) = 1, L(n) = L(n−1) + L(n−2), with terms 2, 1, 3, 4, 7, 11, 18, 29, 47, 76; the entry gives the closed form L(n) = φⁿ + (1−φ)ⁿ and the formula L(n) = ⌊φⁿ⌋ + ((−1)ⁿ + 1)/2. Section 5 uses both to check the case n = 1 by hand.
Wikipedia — Pisot–Vijayaraghavan number. A real algebraic integer greater than 1 whose other conjugates all have absolute value below 1; the sum of the nth powers of such a number and its conjugates is an integer while the conjugate powers tend to 0, so its powers approach integers at an exponential rate. Section 7 identifies our roots as examples.
Wikipedia — Companion matrix. The companion matrix of a monic polynomial carries ones on one diagonal and the negated coefficients in one row or column, its characteristic polynomial is that polynomial, and it advances the state vector of the corresponding linear recurrence. Section 8 uses the variant for a window ordered newest first and verifies the characteristic polynomial directly.
Wikipedia — Exponentiation by squaring. Repeated squaring needs about ⌊log₂ e⌋ squarings for the exponent e and one multiplication fewer than the number of ones in its binary expansion; the article states that the method is often used for powers of matrices and for large exponents in modular arithmetic. Those minimal counts assume that the multiplication by the identity and the unused final square are omitted; the published loops keep both, which gives the 30 squarings and 19 products counted in Section 8.
Wikipedia — Modular arithmetic. Congruence modulo m is compatible with addition, subtraction and multiplication, a residue is invertible modulo m exactly when it is coprime to m, and the Chinese remainder theorem determines a value modulo a product of coprime moduli from its residues. Section 9 relies on the first two statements and the verification above on the third.
Wikipedia — Double-precision floating-point format documents the 53-bit significand of IEEE binary64. cppreference \u2014 Fundamental types describes long double as an extended-precision type that does not necessarily map to an IEEE-754 format, so its precision depends on the platform. Java API — Long gives Long.MAX_VALUE as 2⁶³ − 1, the bound used above, and Python — numeric types states that its integers have unlimited precision.
C++ — Euler356.cpp. The linked immutable revision contains power_matrix3, compute_newton_sum_mod, compute_newton_sum_exact_small, largest_real_root, solve_range and run_validations. It is the implementation record for this article; the wider cross-checks described above were run separately and are not part of that file.

Problem 356 source code

C++

#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <vector>

namespace {

using u32 = std::uint32_t;
using u64 = std::uint64_t;
using u128 = unsigned __int128;

constexpr u64 kMod = 100'000'000ULL;
constexpr u32 kDefaultNMin = 1U;
constexpr u32 kDefaultNMax = 30U;
constexpr u64 kDefaultExponent = 987'654'321ULL;
constexpr std::size_t kAutoThreadThreshold = 96ULL;

struct Options {
    u32 n_min = kDefaultNMin;
    u32 n_max = kDefaultNMax;
    u64 exponent = kDefaultExponent;
    bool run_checks = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
};

bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(c - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_u32_after_prefix(const std::string& arg, const char* prefix, u32& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(std::numeric_limits<u32>::max())) {
        return false;
    }
    value = static_cast<u32>(parsed);
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg, const char* prefix, unsigned& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
        return false;
    }
    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(const int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

        if (arg == "--skip-checks") {
            options.run_checks = false;
            continue;
        }
        if (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }

        u32 parsed_u32 = 0U;
        if (parse_u32_after_prefix(arg, "--n-min=", parsed_u32)) {
            options.n_min = parsed_u32;
            continue;
        }
        if (parse_u32_after_prefix(arg, "--n-max=", parsed_u32)) {
            options.n_max = parsed_u32;
            continue;
        }

        u64 parsed_u64 = 0ULL;
        if (parse_u64_after_prefix(arg, "--exp=", parsed_u64)) {
            options.exponent = parsed_u64;
            continue;
        }

        unsigned parsed_unsigned = 0U;
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
            options.requested_threads = parsed_unsigned;
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    if (options.n_min == 0U || options.n_max == 0U) {
        std::cerr << "n bounds must be >= 1.\n";
        return false;
    }
    if (options.n_min > options.n_max) {
        std::cerr << "Invalid range: --n-min must be <= --n-max.\n";
        return false;
    }
    if (options.exponent == 0ULL || (options.exponent & 1ULL) == 0ULL) {
        std::cerr << "Exponent must be a positive odd integer.\n";
        return false;
    }

    return true;
}

unsigned choose_thread_count(const bool allow_multithreading,
                             const unsigned requested_threads,
                             const std::size_t workload_units) {
    if (!allow_multithreading || workload_units < 2ULL) {
        return 1U;
    }

    const unsigned max_threads = static_cast<unsigned>(
        std::min<std::size_t>(workload_units, static_cast<std::size_t>(std::numeric_limits<unsigned>::max())));

    if (requested_threads > 0U) {
        return std::max(1U, std::min(requested_threads, max_threads));
    }

    if (workload_units < kAutoThreadThreshold) {
        return 1U;
    }

    unsigned hardware_threads = std::thread::hardware_concurrency();
    if (hardware_threads == 0U) {
        hardware_threads = 1U;
    }
    return std::max(1U, std::min(hardware_threads, max_threads));
}

u64 pow_mod_u64(u64 base, u64 exponent, const u64 mod) {
    if (mod == 1ULL) {
        return 0ULL;
    }

    u64 result = 1ULL % mod;
    u64 cur = base % mod;
    u64 exp = exponent;

    while (exp > 0ULL) {
        if ((exp & 1ULL) != 0ULL) {
            result = static_cast<u64>((static_cast<u128>(result) * cur) % mod);
        }
        cur = static_cast<u64>((static_cast<u128>(cur) * cur) % mod);
        exp >>= 1ULL;
    }

    return result;
}

struct Matrix3 {
    u64 v[3][3] = {{0ULL, 0ULL, 0ULL}, {0ULL, 0ULL, 0ULL}, {0ULL, 0ULL, 0ULL}};
};

Matrix3 make_identity_matrix3() {
    Matrix3 out;
    for (int i = 0; i < 3; ++i) {
        out.v[i][i] = 1ULL;
    }
    return out;
}

Matrix3 multiply_matrix3(const Matrix3& a, const Matrix3& b, const u64 mod) {
    Matrix3 out;
    for (int i = 0; i < 3; ++i) {
        for (int j = 0; j < 3; ++j) {
            u128 sum = 0;
            for (int k = 0; k < 3; ++k) {
                sum += static_cast<u128>(a.v[i][k]) * static_cast<u128>(b.v[k][j]);
            }
            out.v[i][j] = static_cast<u64>(sum % mod);
        }
    }
    return out;
}

Matrix3 power_matrix3(Matrix3 base, u64 exponent, const u64 mod) {
    Matrix3 result = make_identity_matrix3();
    u64 exp = exponent;

    while (exp > 0ULL) {
        if ((exp & 1ULL) != 0ULL) {
            result = multiply_matrix3(result, base, mod);
        }
        base = multiply_matrix3(base, base, mod);
        exp >>= 1ULL;
    }

    return result;
}

u64 compute_newton_sum_mod(const u32 n, const u64 exponent, const u64 mod) {
    const u64 c = pow_mod_u64(2ULL, static_cast<u64>(n), mod);
    const u64 s0 = 3ULL % mod;
    const u64 s1 = c;
    const u64 s2 = static_cast<u64>((static_cast<u128>(c) * c) % mod);

    if (exponent == 0ULL) {
        return s0;
    }
    if (exponent == 1ULL) {
        return s1;
    }
    if (exponent == 2ULL) {
        return s2;
    }

    Matrix3 transition;
    transition.v[0][0] = c;
    transition.v[0][1] = 0ULL;
    transition.v[0][2] = (mod - (static_cast<u64>(n) % mod)) % mod;
    transition.v[1][0] = 1ULL % mod;
    transition.v[1][1] = 0ULL;
    transition.v[1][2] = 0ULL;
    transition.v[2][0] = 0ULL;
    transition.v[2][1] = 1ULL % mod;
    transition.v[2][2] = 0ULL;

    const Matrix3 p = power_matrix3(transition, exponent - 2ULL, mod);
    const u128 top = static_cast<u128>(p.v[0][0]) * s2 +
                     static_cast<u128>(p.v[0][1]) * s1 +
                     static_cast<u128>(p.v[0][2]) * s0;
    return static_cast<u64>(top % mod);
}

u64 compute_newton_sum_mod_iterative(const u32 n, const u64 exponent, const u64 mod) {
    const u64 c = pow_mod_u64(2ULL, static_cast<u64>(n), mod);
    const u64 n_mod = static_cast<u64>(n) % mod;
    const u64 s0 = 3ULL % mod;
    const u64 s1 = c;
    const u64 s2 = static_cast<u64>((static_cast<u128>(c) * c) % mod);

    if (exponent == 0ULL) {
        return s0;
    }
    if (exponent == 1ULL) {
        return s1;
    }
    if (exponent == 2ULL) {
        return s2;
    }

    u64 a = s0;
    u64 b = s1;
    u64 cur = s2;
    for (u64 k = 3ULL; k <= exponent; ++k) {
        const u64 lhs = static_cast<u64>((static_cast<u128>(c) * cur) % mod);
        const u64 rhs = static_cast<u64>((static_cast<u128>(n_mod) * a) % mod);
        u64 next = lhs;
        if (next >= rhs) {
            next -= rhs;
        } else {
            next += mod - rhs;
        }
        a = b;
        b = cur;
        cur = next;
    }
    return cur;
}

u128 compute_newton_sum_exact_small(const u32 n, const u32 exponent) {
    u128 c = 1;
    for (u32 i = 0; i < n; ++i) {
        c *= 2U;
    }

    const u128 s0 = 3U;
    const u128 s1 = c;
    const u128 s2 = c * c;

    if (exponent == 0U) {
        return s0;
    }
    if (exponent == 1U) {
        return s1;
    }
    if (exponent == 2U) {
        return s2;
    }

    u128 a = s0;
    u128 b = s1;
    u128 cur = s2;
    for (u32 k = 3U; k <= exponent; ++k) {
        const u128 next = c * cur - static_cast<u128>(n) * a;
        a = b;
        b = cur;
        cur = next;
    }
    return cur;
}

long double evaluate_polynomial(const long double x, const long double c, const u32 n) {
    return x * x * x - c * x * x + static_cast<long double>(n);
}

long double largest_real_root(const u32 n) {
    const long double c = std::ldexp(1.0L, static_cast<int>(n));
    long double low = std::max(0.0L, c - 1.0L);
    long double high = c;

    while (evaluate_polynomial(low, c, n) > 0.0L && low > 0.0L) {
        low = std::max(0.0L, low - 1.0L);
    }

    for (int iter = 0; iter < 220; ++iter) {
        const long double mid = (low + high) * 0.5L;
        if (evaluate_polynomial(mid, c, n) <= 0.0L) {
            low = mid;
        } else {
            high = mid;
        }
    }
    return (low + high) * 0.5L;
}

u64 floor_power(const long double base, const u32 exponent) {
    long double value = 1.0L;
    for (u32 i = 0; i < exponent; ++i) {
        value *= base;
    }
    const long double floored = std::floor(value + 1e-18L);
    return (floored <= 0.0L) ? 0ULL : static_cast<u64>(floored);
}

u64 compute_term_mod(const u32 n, const u64 exponent) {
    const u64 sum_mod = compute_newton_sum_mod(n, exponent, kMod);
    return (sum_mod + kMod - 1ULL) % kMod;
}

u64 solve_range(const u32 n_min,
                const u32 n_max,
                const u64 exponent,
                const bool allow_multithreading,
                const unsigned requested_threads) {
    const std::size_t workload = static_cast<std::size_t>(n_max) - static_cast<std::size_t>(n_min) + 1ULL;
    const unsigned threads = choose_thread_count(allow_multithreading, requested_threads, workload);

    if (threads <= 1U) {
        u64 total = 0ULL;
        for (u32 n = n_min; n <= n_max; ++n) {
            total += compute_term_mod(n, exponent);
            if (total >= kMod) {
                total -= kMod;
            }
        }
        return total;
    }

    std::atomic<u32> next_n(n_min);
    std::vector<u64> partial(static_cast<std::size_t>(threads), 0ULL);
    std::vector<std::thread> pool;
    pool.reserve(static_cast<std::size_t>(threads));

    auto worker = [&](const unsigned tid) {
        u64 local = 0ULL;
        while (true) {
            const u32 n = next_n.fetch_add(1U, std::memory_order_relaxed);
            if (n > n_max) {
                break;
            }
            local += compute_term_mod(n, exponent);
            if (local >= kMod) {
                local -= kMod;
            }
        }
        partial[static_cast<std::size_t>(tid)] = local;
    };

    for (unsigned t = 0U; t < threads; ++t) {
        pool.emplace_back(worker, t);
    }
    for (std::thread& th : pool) {
        th.join();
    }

    u64 total = 0ULL;
    for (const u64 part : partial) {
        total += part;
        if (total >= kMod) {
            total -= kMod;
        }
    }
    return total;
}

u64 solve(const Options& options) {
    return solve_range(options.n_min,
                       options.n_max,
                       options.exponent,
                       options.allow_multithreading,
                       options.requested_threads);
}

bool run_validations(const Options& options) {
    // Checkpoint 1: statement example floor(a_2^2) = 14.
    {
        const long double a2 = largest_real_root(2U);
        const u64 floor_a2_sq = floor_power(a2, 2U);
        if (floor_a2_sq != 14ULL) {
            std::cerr << "Validation failed (1/4): floor(a_2^2) expected 14, got "
                      << floor_a2_sq << '\n';
            return false;
        }
        std::cout << "Validation 1/4 [PASS]\n";
    }

    // Checkpoint 2: matrix recurrence matches direct iterative recurrence.
    {
        for (u32 n = 1U; n <= 12U; ++n) {
            for (u64 k = 0ULL; k <= 220ULL; ++k) {
                const u64 by_matrix = compute_newton_sum_mod(n, k, kMod);
                const u64 by_iter = compute_newton_sum_mod_iterative(n, k, kMod);
                if (by_matrix != by_iter) {
                    std::cerr << "Validation failed (2/4): recurrence mismatch at n=" << n
                              << ", k=" << k << " (matrix=" << by_matrix
                              << ", iterative=" << by_iter << ")\n";
                    return false;
                }
            }
        }
        std::cout << "Validation 2/4 [PASS]\n";
    }

    // Checkpoint 3: floor(a_n^m) == S_m(n)-1 for small odd-exponent spot checks.
    {
        const u32 odd_exponents[] = {1U, 3U, 5U, 7U, 9U, 11U};
        for (u32 n = 1U; n <= 3U; ++n) {
            const long double alpha = largest_real_root(n);
            for (const u32 m : odd_exponents) {
                const u64 floor_numeric = floor_power(alpha, m);
                const u128 exact_sum = compute_newton_sum_exact_small(n, m);
                const u64 floor_via_identity = static_cast<u64>(exact_sum - 1U);
                if (floor_numeric != floor_via_identity) {
                    std::cerr << "Validation failed (3/4): floor identity mismatch at n=" << n
                              << ", m=" << m << " (numeric=" << floor_numeric
                              << ", identity=" << floor_via_identity << ")\n";
                    return false;
                }
            }
        }
        std::cout << "Validation 3/4 [PASS]\n";
    }

    // Checkpoint 4: thread consistency on a larger reduced workload.
    if (options.allow_multithreading) {
        unsigned hardware_threads = std::thread::hardware_concurrency();
        if (hardware_threads == 0U) {
            hardware_threads = 1U;
        }

        if (hardware_threads > 1U) {
            const u32 check_n_min = 1U;
            const u32 check_n_max = 400U;
            const u64 single = solve_range(check_n_min, check_n_max, options.exponent, false, 1U);
            const unsigned compare_threads =
                (options.requested_threads > 1U) ? options.requested_threads : 2U;
            const u64 multi =
                solve_range(check_n_min, check_n_max, options.exponent, true, compare_threads);
            if (single != multi) {
                std::cerr << "Validation failed (4/4): single-thread vs multi-thread mismatch"
                          << " (single=" << single << ", multi=" << multi << ")\n";
                return false;
            }
            std::cout << "Validation 4/4 [PASS]\n";
        } else {
            std::cout << "Validation 4/4 [SKIP: single-core environment]\n";
        }
    } else {
        std::cout << "Validation 4/4 [SKIP: --single-thread enabled]\n";
    }

    return true;
}

void print_usage() {
    std::cerr << "Usage: ./Euler356 [--skip-checks] [--single-thread] [--threads=<unsigned>] "
                 "[--n-min=<u32>] [--n-max=<u32>] [--exp=<odd-u64>]\n";
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        print_usage();
        return 1;
    }

    if (options.run_checks && !run_validations(options)) {
        return 1;
    }

    const u64 answer = solve(options);
    std::cout << std::setfill('0') << std::setw(8) << answer << '\n';
    return 0;
}

Python

def pow_mod(base, exp, mod):
    if mod == 1: return 0
    res = 1
    cur = base % mod
    while exp > 0:
        if exp & 1:
            res = (res * cur) % mod
        cur = (cur * cur) % mod
        exp >>= 1
    return res

def multiply_matrix(a, b, mod):
    out = [[0]*3 for _ in range(3)]
    for i in range(3):
        for j in range(3):
            s = 0
            for k in range(3):
                s += a[i][k] * b[k][j]
            out[i][j] = s % mod
    return out

def power_matrix(base, exp, mod):
    res = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
    while exp > 0:
        if exp & 1:
            res = multiply_matrix(res, base, mod)
        base = multiply_matrix(base, base, mod)
        exp >>= 1
    return res

def compute_newton_sum_mod(n, exp, mod):
    c = pow_mod(2, n, mod)
    s0 = 3 % mod
    s1 = c
    s2 = (c * c) % mod
    
    if exp == 0: return s0
    if exp == 1: return s1
    if exp == 2: return s2
    
    transition = [
        [c, 0, (mod - (n % mod)) % mod],
        [1, 0, 0],
        [0, 1, 0]
    ]
    p = power_matrix(transition, exp - 2, mod)
    top = p[0][0] * s2 + p[0][1] * s1 + p[0][2] * s0
    return top % mod

def compute_term_mod(n, exp, mod):
    sum_mod = compute_newton_sum_mod(n, exp, mod)
    return (sum_mod + mod - 1) % mod

def solve():
    mod = 100000000
    total = 0
    for n in range(1, 31):
        total = (total + compute_term_mod(n, 987654321, mod)) % mod
    return str(total)

if __name__ == '__main__':
    print(solve())

Java

public class Euler356 {
    static long powMod(long base, long exp, long mod) {
        if (mod == 1)
            return 0;
        long res = 1;
        long cur = base % mod;
        while (exp > 0) {
            if ((exp & 1) != 0)
                res = (res * cur) % mod;
            cur = (cur * cur) % mod;
            exp >>= 1;
        }
        return res;
    }

    static long[][] multiplyMatrix(long[][] a, long[][] b, long mod) {
        long[][] out = new long[3][3];
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++) {
                long s = 0;
                for (int k = 0; k < 3; k++) {
                    s += a[i][k] * b[k][j];
                }
                out[i][j] = s % mod;
            }
        }
        return out;
    }

    static long[][] powerMatrix(long[][] base, long exp, long mod) {
        long[][] res = { { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } };
        while (exp > 0) {
            if ((exp & 1) != 0)
                res = multiplyMatrix(res, base, mod);
            base = multiplyMatrix(base, base, mod);
            exp >>= 1;
        }
        return res;
    }

    static long computeNewtonSumMod(int n, long exp, long mod) {
        long c = powMod(2, n, mod);
        long s0 = 3 % mod;
        long s1 = c;
        long s2 = (c * c) % mod;

        if (exp == 0)
            return s0;
        if (exp == 1)
            return s1;
        if (exp == 2)
            return s2;

        long[][] transition = {
                { c, 0, (mod - (n % mod)) % mod },
                { 1, 0, 0 },
                { 0, 1, 0 }
        };
        long[][] p = powerMatrix(transition, exp - 2, mod);
        long top = p[0][0] * s2 + p[0][1] * s1 + p[0][2] * s0;
        return top % mod;
    }

    public static String solve() {
        long mod = 100000000;
        long total = 0;
        for (int n = 1; n <= 30; n++) {
            long term = (computeNewtonSumMod(n, 987654321, mod) + mod - 1) % mod;
            total = (total + term) % mod;
        }
        return String.format("%08d", total);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}