Problem 402: Integer-valued Polynomials

Project Euler Problem 402 Solution

EulerSolve provides an optimized solution for Project Euler Problem 402, Integer-valued Polynomials, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Define $$P_{a,b,c}(n)=n^4+a n^3+b n^2+c n.$$ For each positive triple $(a,b,c)$, let $M(a,b,c)$ be the largest integer that divides $P_{a,b,c}(n)$ for every integer $n$. Then $$S(N)=\sum_{1\le a,b,c\le N} M(a,b,c).$$ The final target is $$T(K)=\sum_{k=2}^{K} S(F_k),\qquad K=1234567890123,$$ where $F_k$ denotes the Fibonacci sequence. The implementations return $T(K)\bmod 10^9$, i.e. the last nine digits. Mathematical Approach 1. Turn the divisibility condition into an integer-valued polynomial problem The crucial fact is that a polynomial is integer-valued on all integers if and only if it has integer coefficients in the binomial basis $\binom{n}{1},\binom{n}{2},\binom{n}{3},\binom{n}{4},\dots$. So instead of working with powers $n^j$, we rewrite $P_{a,b,c}(n)$ in that basis. The standard identities are $$n=\binom{n}{1},\qquad n^2=2\binom{n}{2}+\binom{n}{1},$$ $$n^3=6\binom{n}{3}+6\binom{n}{2}+\binom{n}{1},$$ $$n^4=24\binom{n}{4}+36\binom{n}{3}+14\binom{n}{2}+\binom{n}{1}.$$ Substituting these into $P_{a,b,c}(n)$ gives $$P_{a,b,c}(n)=24\binom{n}{4}+(36+6a)\binom{n}{3}+(14+6a+2b)\binom{n}{2}+(1+a+b+c)\binom{n}{1}.$$ 2. Closed form for $M(a,b,c)$ Each binomial polynomial $\binom{n}{j}$ is integer-valued for every integer $n$....

Detailed mathematical approach

Problem Summary

For integers $a,b,c$ put

$$P(n)=n^{4}+an^{3}+bn^{2}+cn.$$

Then $M(a,b,c)$ is the largest integer $m$ dividing $P(n)$ for every integer $n$; the statement's example is $M(4,2,5)=6$, so $n^4+4n^3+2n^2+5n$ is always a multiple of $6$ and of nothing larger. Next, $S(N)$ collects these values over all triples with $0\lt a,b,c\le N$,

$$S(N)=\sum_{1\le a,b,c\le N}M(a,b,c),$$

and the task is the last nine digits of

$$T=\sum_{k=2}^{K}S(F_k),\qquad K=1234567890123,$$

where $F_0=0$, $F_1=1$ and $F_k=F_{k-1}+F_{k-2}$. The statement supplies two checkpoints, $S(10)=1972$ and $S(10000)=2024258331114$. ^[1]

The scale is worth a moment: $F_K$ has about $2.58\cdot10^{11}$ decimal digits and the sum has over $10^{12}$ terms, so everything must be done without writing an $F_k$ down.

The whole solution in one sentence: $M(a,b,c)$ is the greatest common divisor of the four coefficients of $P$ in the basis of binomial coefficients, and that divisor depends only on $a,b,c$ modulo $24$, which makes $S(N)$ a cubic polynomial in $\lfloor N/24\rfloor$ with coefficients depending only on $N\bmod 24$; since the Fibonacci residues modulo $24$ repeat with period $24$, the quotients $\lfloor F_k/24\rfloor$ obey a fixed linear recurrence and the running total follows from one $11\times11$ matrix raised to a large power.

Mathematical Approach

1. The largest divisor is a greatest common divisor

The largest $m$ dividing every $P(n)$ is the greatest common divisor of them all:

$$M(a,b,c)=\gcd\{\,P(n):n\in\mathbb Z\,\}.$$

That already gives an upper bound. For $a=4$, $b=2$, $c=5$ the first values are $P(1)=12$, $P(2)=66$, $P(3)=222$, $P(4)=564$; as $\gcd(12,66)=6$ and $6$ divides $222$ and $564$, the answer is at most $6$. But no program loops over infinitely many values, so we need a finite description of the same quantity.

2. Newton's basis of binomial coefficients

Here $\binom nj$ always means the polynomial

$$\binom nj=\frac{n(n-1)\cdots(n-j+1)}{j!},\qquad\binom n0=1.$$

For $n\ge0$ this is the ordinary binomial coefficient; for $n=-u$ with $u\ge1$ the product gives $\binom{-u}{j}=(-1)^j\binom{u+j-1}{j}$, an integer too — which matters, since $P$ must be controlled at negative arguments. They form a basis of the rational polynomials, the one of degree $j$ having leading coefficient $1/j!$. ^[2] In that basis

$$\begin{aligned} n&=\binom n1,\\ n^2&=2\binom n2+\binom n1,\\ n^3&=6\binom n3+6\binom n2+\binom n1,\\ n^4&=24\binom n4+36\binom n3+14\binom n2+\binom n1. \end{aligned}$$

These are not coincidences. Let $\sigma(i,j)$, a Stirling number of the second kind, count the ways to split $i$ labelled positions into $j$ nonempty groups; grouping equal labels and then giving the groups distinct labels turns the $n^i$ labellings into $n^i=\sum_j \sigma(i,j)\,j!\binom nj$. So the coefficient of $\binom nj$ in $n^i$ is $j!\,\sigma(i,j)$, and the row $\sigma(4,1),\dots,\sigma(4,4)=1,7,6,1$ gives $1,14,36,24$. ^[3]

Substituting into $P$ and collecting terms gives

$$P(n)=c_4\binom n4+c_3\binom n3+c_2\binom n2+c_1\binom n1$$

with four integer coefficients

$$\begin{aligned} c_4&=24, & c_3&=6a+36,\\ c_2&=6a+2b+14, & c_1&=a+b+c+1. \end{aligned}$$

3. A closed form for $M(a,b,c)$

Claim. $\gcd\{P(n)\}=\gcd(c_1,c_2,c_3,c_4)$. ^[4]

Proof. One direction is immediate: every $\binom nj$ is an integer, so each $P(n)$ is an integer combination of $c_1,\dots,c_4$, hence divisible by their greatest common divisor. For the other, let $\Delta f(n)=f(n+1)-f(n)$ be the forward difference and, for $(a,b,c)=(4,2,5)$, subtract adjacent entries of $P(0),\dots,P(4)$ repeatedly:

$$\begin{aligned} &0,\;12,\;66,\;222,\;564,\\ &12,\;54,\;156,\;342,\\ &42,\;102,\;186,\\ &60,\;84,\\ &24. \end{aligned}$$

The first entry of each new row is $c_1=12$, $c_2=42$, $c_3=60$, $c_4=24$ — no accident, since Pascal's identity gives $\Delta\binom nj=\binom n{j-1}$, so differencing the expansion of Section 2 and setting $n=0$ leaves $c_j=\Delta^jP(0)$. Iterating $\Delta$ gives

$$c_j=\sum_{i=0}^{j}(-1)^{j-i}\binom ji P(i).$$

Each $c_j$ is thus an integer combination of $P(0),\dots,P(j)$, so the two families generate the same integer combinations and their greatest common divisors agree. $\square$

Since $P(0)=0$, four values already decide everything:

$$M(a,b,c)=\gcd\bigl(P(1),P(2),P(3),P(4)\bigr).$$

Three are not enough: for $(a,b,c)=(1,4,1)$ the values are $7,42,147,388$, the first three with greatest common divisor $7$, yet $\gcd(7,388)=1=M(1,4,1)$. With Section 2's coefficients the claim becomes the formula the solution rests on:

$$M(a,b,c)=\gcd(24,\,6a+36,\,6a+2b+14,\,a+b+c+1).$$

It reproduces the statement at once: $\gcd(24,60,42,12)=6$, so $M(4,2,5)=6$; likewise $M(1,1,1)=\gcd(24,42,22,4)=2$, which is $S(1)$. Read differently, it is the standard fact about integer-valued polynomials: $m$ divides every $P(n)$ exactly when $P/m$ is integer-valued, that is when every $c_j/m$ is an integer. ^[2]

4. Only the residues modulo $24$ matter

Two consequences follow. First, $24$ is one of the four arguments, so $M(a,b,c)$ always divides $24$: the warm-up value $6$ is simply one of the eight divisors of $24$. Second, $\gcd(24,y)$ depends only on $y\bmod 24$, so if $a\equiv a'$, $b\equiv b'$ and $c\equiv c'$ modulo $24$, then $M(a,b,c)=M(a',b',c')$.

The first argument has extra structure: since

$$\gcd(24,6a+36)=6\gcd(4,a+2),$$

it is $24$ when $a\equiv2\pmod 4$, $12$ when $a\equiv0\pmod 4$ and $6$ when $a$ is odd — an odd $a$ alone caps $M$ at $6$.

5. Which values $M$ takes

By Section 4 it is enough to tabulate $M$ on the $24^3=13824$ residue triples with $1\le a,b,c\le24$. The value $1$ occurs $6144$ times, $2$ occurs $5376$, $3$ occurs $768$, $4$ and $6$ occur $672$ each, $8$ occurs $96$, $12$ occurs $84$, $24$ occurs $12$: total $27984$, average $27984/13824=583/288\approx2.0243$.

Two hand checks confirm the table. The arguments $24$, $6a+36$, $6a+2b+14$ are even, so $M$ is odd exactly when $c_1=a+b+c+1$ is odd, for half the triples, $6912$; the odd divisors of $24$ are $1$ and $3$, and $6144+768=6912$. Likewise $3$ divides $6a+36$ always, $6a+2b+14$ when $b\equiv2\pmod 3$ and $a+b+c+1$ when $a+b+c\equiv2\pmod 3$, leaving $24\cdot8\cdot8=1536$ triples, matching $768+672+84+12$. The maximum $24$ is reached by $12$ residue triples; among representatives with $1\le a,b,c\le24$ the lexicographically smallest is $(2,11,10)$.

6. $S(N)$ is a cubic quasi-polynomial

Write $N=24q+r$ with $0\le r\lt24$. For $1\le t\le24$, the integers $t,t+24,t+48,\dots$ not exceeding $N$ number

$$q+\varepsilon_t,\qquad \varepsilon_t=\begin{cases}1,&t\le r,\\0,&t\gt r,\end{cases}$$

because $N-t=24q+(r-t)$ when $t\le r$ and $N-t=24(q-1)+(24+r-t)$ otherwise. Group the triple sum by the residues $u,v,w$ of $a,b,c$, abbreviating $x=\varepsilon_u$, $y=\varepsilon_v$, $z=\varepsilon_w$:

$$\begin{aligned} S(N)=\sum_{u,v,w=1}^{24}M(u,v,w)\\ \times(q+x)(q+y)(q+z). \end{aligned}$$

The inner product expands with the elementary symmetric functions $e_1=x+y+z$, $e_2=xy+xz+yz$ and $e_3=xyz$:

$$(q+x)(q+y)(q+z)=q^3+e_1q^2+e_2q+e_3.$$

As $q$ is independent of the residues, the sum splits in four:

$$S(N)=A_3q^3+A_2(r)q^2+A_1(r)q+A_0(r),$$

$$A_{3-i}(r)=\sum_{u,v,w=1}^{24}M(u,v,w)\,e_i,$$

with $e_0=1$, so $A_3=27984$ is the constant of Section 5. Substituting $q=(N-r)/24$ turns this into a cubic polynomial in $N$, one per remainder, with coefficients repeating every $24$: such a polynomial with periodic coefficients is a quasi-polynomial, here of degree $3$. ^[5] Nothing below needs that period to be minimal. Taking $q=0$ gives $A_0(r)=S(r)$ for $0\le r\lt24$, and $S(N)$ grows like $(583/288)N^3$.

7. Worked example: the two checkpoints and $S(25)$

Each coefficient is a sum over the residue table. For $r=1$, $\varepsilon_t=1$ only at $t=1$, so $A_0(1)=M(1,1,1)=2$; $A_1(1)$ adds the one-variable sums $\sum_cM(1,1,c)$, $\sum_bM(1,b,1)$, $\sum_aM(a,1,1)$, each $36$, giving $108$; $A_2(1)$ adds the two-variable sums $1056$, $1008$, $1056$, giving $3120$. For $N=25$, $q=1$ and $r=1$:

$$S(25)=27984+3120+108+2=31214,$$

which a direct loop over $25^3$ triples confirms. The first checkpoint needs even less: $10=24\cdot0+10$, so $S(10)=A_0(10)=1972$. The second exercises all four terms: $10000=24\cdot416+16$, with $A_2(16)=55696$, $A_1(16)=36944$, $A_0(16)=S(16)=8170$, so the summands are $2014604427264$, $9638526976$, $15368704$ and $8170$:

$$\begin{aligned} S(10000)&=27984\cdot416^{3}+55696\cdot416^{2}\\ &\quad+36944\cdot416+8170\\ &=2024258331114, \end{aligned}$$

exactly the value in the statement. Evaluating $S$ now costs a handful of multiplications however large the argument is, provided we can reach $\lfloor N/24\rfloor$ and $N\bmod 24$.

8. Fibonacci numbers modulo $24$

The residues $r_k=F_k\bmod 24$ are easy. Modulo any $m$ the pair of consecutive residues has finitely many states, and its update $(x,y)\mapsto(y,x+y)$ is reversible by $(x,y)\mapsto(y-x,x)$, so the starting pair lies on a cycle and the sequence is periodic from the beginning. The cycle length is the Pisano period $\pi(m)$; for coprime $m,n$ a pair returns modulo $mn$ exactly when it returns modulo both factors, so $\pi(mn)=\operatorname{lcm}(\pi(m),\pi(n))$, a least common multiple and not a product. Thus $\pi(8)=12$ and $\pi(3)=8$ give $\pi(24)=24$. ^[6] ^[7] The cycle, for $k=0,1,\dots,23$, is

$$\begin{aligned} &0,1,1,2,3,5,8,13,21,10,7,17,\\ &0,17,17,10,3,13,16,5,21,2,23,1, \end{aligned}$$

and then it repeats. Note $r_{12}=0$: $F_{12}=144$ is a multiple of $24$, so $S(144)=27984\cdot6^3=6044544$, the three remainder-dependent terms vanishing.

9. The quotients obey a fixed recurrence

Write $F_k=24q_k+r_k$. Adding the decompositions of $F_k$ and $F_{k+1}$ and comparing with that of $F_{k+2}$ gives

$$24q_{k+2}+r_{k+2}=24(q_k+q_{k+1})+r_k+r_{k+1}.$$

As $r_k,r_{k+1}\lt24$, the sum $r_k+r_{k+1}$ is $r_{k+2}$ or $r_{k+2}+24$, so

$$q_{k+2}=q_{k+1}+q_k+\delta_k,\qquad \delta_k=\Bigl\lfloor\frac{r_k+r_{k+1}}{24}\Bigr\rfloor.$$

The carry $\delta_k$ is $0$ or $1$ and depends only on $k\bmod 24$, as the residues do. Nine of the $24$ phases carry:

$$k\equiv7,8,10,13,14,17,19,21,22\pmod{24}.$$

The starting values are $q_2=q_3=0$, since $F_2=1$ and $F_3=2$. The enormous Fibonacci numbers are thus replaced by a three-term recurrence with a periodic inhomogeneity — one a matrix can carry.

10. Linearising a cubic in two variables

Take the step $k=8$ first. From $F_8=21$ and $F_9=34$ the pair $(q_8,q_9)$ is $(0,1)$ and the remainders are $21$ and $10$, so the carry is $1$ and the next pair is $(q_9,q_{10})=(1,2)$, while the same step adds $S(21)=18560$, lifting the total from $5749$ to $24309$. Quotient powers and total advance together.

By Section 6, $S(F_k)$ is a cubic polynomial in $q_k$ with coefficients fixed by $k\bmod 24$, so the state is the pair $(q,q')=(q_k,q_{k+1})$ with enough of its powers, and the step is the affine map $(q,q')\mapsto(q',\,q+q'+\delta)$. Substituting an affine map into a monomial of degree at most $3$ gives a combination of such monomials, for instance

$$\begin{aligned} (q+q'+\delta)^3&=(q+q')^3+3\delta(q+q')^2\\ &\quad+3\delta^2(q+q')+\delta^3, \end{aligned}$$

and each piece expands into the monomials $q^3,q^2q',qq'^2,q'^3$, $q^2,qq',q'^2$, $q,q'$ and $1$. There are exactly $\binom{5}{2}=10$ monomials of degree at most $3$ in two variables, so ten coordinates suffice. One more holds the running total $\sum_{j=2}^{k-1}S(F_j)$, which is $0$ at $k=2$ and whose update $T\mapsto T+S(F_k)$ is linear in those ten and in $T$: dimension $11$. One step is thus a single $11\times11$ matrix $V_p$, one per phase $p=k\bmod 24$.

11. One cycle, one matrix power

The phase advances by one each step and repeats after $24$, so the ordered product of the $24$ matrices from phase $2$,

$$B=V_1V_0V_{23}\cdots V_3V_2,$$

advances the state by a whole cycle and returns it to phase $2$, so the same $B$ can be applied repeatedly. The steps $k=2,3,\dots,K$ number $K-1=1234567890122$, and

$$1234567890122=24\cdot51440328755+2,$$

so the computation is $B^{51440328755}$ followed by $V_2$ and $V_3$. The power comes from repeated squaring, one squaring per bit of the exponent and one extra multiplication per $1$ bit: here $36$ bits and $17$ one bits. ^[8] Everything is reduced modulo $10^9$ at every step, legitimate because the update uses only integer additions and multiplications and only the last nine digits are wanted.

12. Worked example: the first nine terms

For $k=2,\dots,10$ the Fibonacci numbers are $1,2,3,5,8,13,21,34,55$, the quotients $q_k$ are $0,0,0,0,0,0,0,1,2$ and the remainders $r_k$ are $1,2,3,5,8,13,21,10,7$, so the terms are

$$2,20,57,266,1056,4348,18560,79052,335311.$$

The eighth is a full four-term evaluation, as $F_9=34$ gives $q=1$ and $r=10$:

$$S(34)=27984+34736+14360+1972=79052.$$

Their total is $438672$, and the matrix method with $K=10$ prints $000438672$: the leading zeros belong to the answer format and must not be dropped.

13. Two routes that look shorter but are not

Enumerating a cycle looks tempting, since $S(N)\bmod 10^9$ is determined by $N\bmod 24\cdot10^9=2^{12}\cdot3\cdot5^9$, and $\pi(2^{12})=6144$, $\pi(3)=8$, $\pi(5^9)=4\cdot5^9$ make those residues cycle after $1.2\cdot10^{10}$ steps. That is a period of the summands but need not be the shortest: $A_3$ and every $A_2(r)$, $A_1(r)$ is divisible by $4$, so raising $q$ by $10^9/4$ leaves $S$ unchanged modulo $10^9$ and $\pi(6\cdot10^9)=3\cdot10^9$ is a period too — still billions of terms against a few dozen small matrix products. The identity $5F_n^3=F_{3n}-3(-1)^nF_n$ offers another route, turning $\sum q_k^3$ into Fibonacci values at tripled indices that fast doubling supplies in logarithmic time; but $A_2,A_1,A_0$ depend on the remainder, forcing a split into $24$ classes of $k$ — bookkeeping the single matrix power does by itself.

How the Code Works

All three programs follow four stages: tabulate the coefficients, build the $24$ phase matrices, raise one cycle to a power, read one coordinate. The walkthrough uses the C++ names; Python shares them but hands the four coefficient tables to step_matrix separately, Java spells them buildCoefficients, stepMatrix, matMul, matPow and matVecMul, and only C++ carries the checkpoint helper s_mod. gcd4 runs Euclid's algorithm on $24$, 36 + 6*ra, 14 + 6*ra + 2*rb and 1 + ra + rb + rc — Section 3's formula transcribed literally.

build_coefficients is Section 6 in one pass: it loops ra, rb, rc over $1,\dots,24$, computes m once per triple, adds it to a3, then loops rem over $0,\dots,23$ to accumulate m times $e_1$, $e_2$, $e_3$ into a2[rem], a1[rem], a0[rem]. The symmetric functions appear verbatim as ea + eb + ec, ea*eb + ea*ec + eb*ec and ea*eb*ec, each ea an indicator $\varepsilon$, and all four tables stay modulo $10^9$. s_mod evaluates the quasi-polynomial at one argument, splitting it into q = n / 24 and rem = n % 24; only the checkpoints call it.

step_matrix(carry, rem, coef) builds one $V_p$. A comment records the coordinate order: 0:1, 1:q, 2:q1, 3:q^2, 4:q q1, 5:q1^2, 6:q^3, 7:q^2 q1, 8:q q1^2, 9:q1^3, 10:sum. Row $2$ is Section 9's recurrence $q'=q+q_1+c$, row $9$ holds the ten terms of $(q+q_1+c)^3$ from Section 10, and row $10$ is the accumulator with coef.a0[rem], coef.a1[rem], coef.a2[rem], coef.a3 and a $1$ on the diagonal. By phase a matrix has only $22$, $25$ or $40$ nonzero entries out of $121$, and mat_mul skips zero factors in both operands.

solve assembles the rest. It fills fib_mod24 with $F_i\bmod 24$ for $i\lt24$ and derives carry[i] as (fib_mod24[i] + fib_mod24[(i+1)%24]) / 24; the wrap-around at i = 23 reads $F_0\bmod 24$ for $F_{24}\bmod 24$, correct exactly because the period is $24$. It forms block as the ordered product of the phase matrices from start_phase = 2 and sets state to the vector whose only nonzero entry is the constant coordinate — the encoding of $q_2=q_3=0$ with an empty total. With steps = k_max - 1 it applies mat_pow(block, steps / 24), then steps % 24 single phase matrices, and returns coordinate $10$.

run_checkpoints guards the algebra: it compares s_mod(10, coef) with $1972$ and s_mod(10000, coef) with $2024258331114\bmod 10^9$, then recomputes $S(10000)$ exactly in __int128 arithmetic from the three residue counts. parse_arguments accepts --k-max= and --skip-checkpoints, making Section 12's small cases runnable, and main prints with std::setw(9) and std::setfill('0').

The languages differ only in arithmetic and packaging. In every multiply-and-add of mat_mul both factors and the accumulator are already reduced modulo $10^9$ — the reduction after each addition is what makes that true — so the value stays below $10^{18}$, within reach of C++'s std::uint64_t and of Java's signed long, which holds $2^{63}-1\approx9.22\cdot10^{18}$. The C++ code still promotes products to __uint128_t, a safety margin rather than a necessity, and Python's integers are unbounded. ^[9] ^[10] All three pad to nine digits — "{:09d}" in Python, String.format("%09d") in Java — since the last nine digits may begin with zeros, as $K=10$ shows. Python fixes k_max and Java kMax at $1234567890123$, neither running checkpoints; only C++ has run_checkpoints and the two flags.

Complexity and Verification

The coefficient stage is constant work: $24^3=13824$ evaluations of gcd4 and $24^4=331776$ inner iterations, each updating three coefficient tables, independent of $K$. The matrix stage is $O(\log K)$ products of a fixed $11\times11$ matrix, each at most $11^3=1331$ modular multiply-and-adds and fewer wherever zero factors are skipped. Concretely $24$ products build the cycle matrix, $36$ squarings and $17$ multiplications raise it to the power $51440328755$, and three matrix–vector products finish the state — $77$ matrix products in all, near $10^5$ modular operations. Memory is $O(1)$ in $K$ but not negligible — the $24$ phase matrices alone hold $2904$ entries, $23232$ bytes in C++ — and on an ordinary desktop all three programs finish in a fraction of a second.

Beyond the packaged checkpoints, separate programs verified the following for this article; none is part of the published sources.

Section 3's closed form equals the greatest common divisor of $P(n)$ over $0\le n\le59$ for all $64000$ triples with $1\le a,b,c\le40$, and equals $\gcd(P(1),P(2),P(3),P(4))$ for all $27000$ triples with $1\le a,b,c\le30$, while three values fail for $39$ of the $1728$ triples with $1\le a,b,c\le12$.
A direct triple loop over $1\le a,b,c\le N$ agrees with Section 6 for every $N\le199$.
A second evaluation of $S(N)$ by exact residue counts, without the $q$ and $\varepsilon$ split, agrees with Section 6 for every $N\le399$ and for $F_k$ with $2\le k\le69$, and gives $2024258331114$ at $N=10000$.
The recurrence $q_{k+2}=q_{k+1}+q_k+\delta_{k\bmod 24}$ agrees with exact Fibonacci quotients for all $k\le249$; the periods $\pi(5^9)=4\cdot5^9$ and $\pi(6\cdot10^9)=3\cdot10^9$ of Section 13, and the divisibility of $A_3,A_2,A_1$ by $4$, were confirmed directly.
The matrix method reproduces the exact big-integer running total modulo $10^9$ for every $K$ with $2\le K\le600$ and for $K=1000$; the C++ program with --k-max= matches that reference at twelve limits up to $1000$, giving $000438672$ at $K=10$ and $854524944$ at $K=60$. At the full $K$ all three programs print the same nine digits.

Footnotes and References

The references supply the background facts used above; the closed form for $M$, the residue tables, the quasi-polynomial, the carry recurrence and the eleven-coordinate linearisation are derived explicitly here.

Project Euler 402 — Integer-valued Polynomials. The official statement defines M(a, b, c) as the maximum m such that n⁴ + an³ + bn² + cn is a multiple of m for all integers n, notes M(4, 2, 5) = 6, sums it as S(N) over 0 < a, b, c ≤ N, gives S(10) = 1972 and S(10000) = 2024258331114, and asks for the last 9 digits of the sum of S(F_k) over 2 ≤ k ≤ 1234567890123.
Wikipedia — Integer-valued polynomial. The integer-valued polynomials form a free abelian group inside the rational polynomials with the binomial coefficients t(t−1)⋯(t−k+1)/k! as a basis, every one of them being an integer combination of that basis in exactly one way — the equivalence used in Sections 2 and 3.
Wikipedia — Stirling numbers of the second kind. It gives the identity summing these numbers, written S(n, k) there and σ(n, k) above, against the falling factorial (x)_k = x(x−1)⋯(x−k+1) to produce xⁿ, and lists the row 1, 7, 6, 1 for n = 4 behind the coefficients 1, 14, 36, 24.
Wikipedia — Finite difference. It defines the forward difference, gives the n-th difference as an alternating sum of binomial coefficients times function values, and states the Newton series, which writes a polynomial as the sum over k of C(x − a, k) times the k-th difference at a — both used in Section 3.
Wikipedia — Quasi-polynomial. A quasi-polynomial is a polynomial whose coefficients are periodic functions of the variable with integral period, its degree being the largest power whose coefficient is not identically zero. The equivalent description by constituent polynomials, one per residue class, is the form Section 6 produces.
Wikipedia — Pisano period. The Pisano period π(n) is the cycle length of the Fibonacci sequence modulo n; the article gives π(mn) = lcm(π(m), π(n)) for coprime m and n, π(2^k) = 3·2^(k−1), π(3) = 8, π(5^k) = 4·5^k and π(24) = 24, the formulas used in Sections 8 and 13.
OEIS A001175 — Pisano periods. The sequence lists the period of the Fibonacci numbers modulo n; its 24th term is 24, confirming the value used in Section 8.
Wikipedia — Exponentiation by squaring. The method needs about ⌊log₂ n⌋ squarings and, in the variant described there, one fewer multiplication than the number of 1 bits of the exponent; it works in non-commutative semigroups and is often used for powers of matrices. Section 11's counts follow the published loop instead, which starts from the identity and squares at every bit, the last included.
C++ — Euler402.cpp. The linked immutable revision holds build_coefficients, step_matrix, solve and run_checkpoints, and is the implementation record for this article; the extra comparisons above are not packaged in it.
GCC manual — 128-bit integers documents unsigned __int128 on targets with an integer mode wide enough for 128 bits. Java — Long gives MAX_VALUE as 2⁶³ − 1, and Python — numeric types states that Python integers have unlimited precision.

Problem 402 source code

C++

#include <array>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>

namespace {

using u64 = std::uint64_t;
using i64 = long long;

constexpr u64 kMod = 1000000000ULL;  // last 9 digits
constexpr int kDim = 11;

struct Options {
    i64 k_max = 1234567890123LL;
    bool run_checkpoints = true;
};

bool parse_i64_after_prefix(const std::string& arg, const std::string& prefix, i64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        value = std::stoll(tail);
    } catch (...) {
        return false;
    }
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_i64_after_prefix(arg, "--k-max=", options.k_max)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.k_max >= 2;
}

i64 gcd4(i64 a, i64 b, i64 c, i64 d) {
    auto g2 = [](i64 x, i64 y) {
        while (y != 0) {
            const i64 t = x % y;
            x = y;
            y = t;
        }
        return x < 0 ? -x : x;
    };
    return g2(g2(a, b), g2(c, d));
}

struct Coefs {
    u64 a3 = 0;  // coefficient of q^3
    std::array<u64, 24> a2{};  // coefficient of q^2 by remainder
    std::array<u64, 24> a1{};  // coefficient of q   by remainder
    std::array<u64, 24> a0{};  // constant by remainder
};

Coefs build_coefficients() {
    Coefs out;

    for (int ra = 1; ra <= 24; ++ra) {
        for (int rb = 1; rb <= 24; ++rb) {
            for (int rc = 1; rc <= 24; ++rc) {
                const i64 m = gcd4(
                    24,
                    36 + 6 * ra,
                    14 + 6 * ra + 2 * rb,
                    1 + ra + rb + rc);

                out.a3 = (out.a3 + static_cast<u64>(m)) % kMod;
                for (int rem = 0; rem < 24; ++rem) {
                    const int ea = (ra <= rem ? 1 : 0);
                    const int eb = (rb <= rem ? 1 : 0);
                    const int ec = (rc <= rem ? 1 : 0);
                    const u64 u = static_cast<u64>(m) % kMod;
                    out.a2[static_cast<std::size_t>(rem)] =
                        (out.a2[static_cast<std::size_t>(rem)] + u * static_cast<u64>(ea + eb + ec)) % kMod;
                    out.a1[static_cast<std::size_t>(rem)] =
                        (out.a1[static_cast<std::size_t>(rem)] + u * static_cast<u64>(ea * eb + ea * ec + eb * ec)) % kMod;
                    out.a0[static_cast<std::size_t>(rem)] =
                        (out.a0[static_cast<std::size_t>(rem)] + u * static_cast<u64>(ea * eb * ec)) % kMod;
                }
            }
        }
    }

    return out;
}

u64 s_mod(i64 n, const Coefs& coef) {
    const i64 q = n / 24;
    const int rem = static_cast<int>(n % 24);
    const u64 qq = static_cast<u64>(q % static_cast<i64>(kMod));
    const u64 q2 = static_cast<u64>((__uint128_t)qq * qq % kMod);
    const u64 q3 = static_cast<u64>((__uint128_t)q2 * qq % kMod);

    u64 ans = 0;
    ans = (ans + static_cast<u64>((__uint128_t)coef.a3 * q3 % kMod)) % kMod;
    ans = (ans + static_cast<u64>((__uint128_t)coef.a2[static_cast<std::size_t>(rem)] * q2 % kMod)) % kMod;
    ans = (ans + static_cast<u64>((__uint128_t)coef.a1[static_cast<std::size_t>(rem)] * qq % kMod)) % kMod;
    ans = (ans + coef.a0[static_cast<std::size_t>(rem)]) % kMod;
    return ans;
}

struct Mat {
    std::array<std::array<u64, kDim>, kDim> a{};
};

Mat mat_identity() {
    Mat m;
    for (int i = 0; i < kDim; ++i) {
        m.a[static_cast<std::size_t>(i)][static_cast<std::size_t>(i)] = 1ULL;
    }
    return m;
}

Mat mat_mul(const Mat& x, const Mat& y) {
    Mat z;
    for (int i = 0; i < kDim; ++i) {
        for (int k = 0; k < kDim; ++k) {
            const u64 xik = x.a[static_cast<std::size_t>(i)][static_cast<std::size_t>(k)];
            if (xik == 0ULL) {
                continue;
            }
            for (int j = 0; j < kDim; ++j) {
                const u64 ykj = y.a[static_cast<std::size_t>(k)][static_cast<std::size_t>(j)];
                if (ykj == 0ULL) {
                    continue;
                }
                z.a[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)] =
                    (z.a[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)] +
                     static_cast<u64>((__uint128_t)xik * ykj % kMod)) %
                    kMod;
            }
        }
    }
    return z;
}

Mat mat_pow(Mat base, i64 exp) {
    Mat res = mat_identity();
    i64 e = exp;
    while (e > 0) {
        if (e & 1LL) {
            res = mat_mul(base, res);
        }
        base = mat_mul(base, base);
        e >>= 1LL;
    }
    return res;
}

std::array<u64, kDim> mat_vec_mul(const Mat& m, const std::array<u64, kDim>& v) {
    std::array<u64, kDim> out{};
    for (int i = 0; i < kDim; ++i) {
        u64 s = 0ULL;
        for (int j = 0; j < kDim; ++j) {
            const u64 mij = m.a[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)];
            if (mij == 0ULL) {
                continue;
            }
            s = (s + static_cast<u64>((__uint128_t)mij * v[static_cast<std::size_t>(j)] % kMod)) % kMod;
        }
        out[static_cast<std::size_t>(i)] = s;
    }
    return out;
}

Mat step_matrix(int carry, int rem, const Coefs& coef) {
    Mat m{};
    const u64 c = static_cast<u64>(carry) % kMod;
    const u64 c2 = static_cast<u64>((__uint128_t)c * c % kMod);
    const u64 c3 = static_cast<u64>((__uint128_t)c2 * c % kMod);

    // Basis:
    // 0:1, 1:q, 2:q1, 3:q^2, 4:q q1, 5:q1^2, 6:q^3, 7:q^2 q1, 8:q q1^2, 9:q1^3, 10:sum
    m.a[0][0] = 1;

    m.a[1][2] = 1;                // q' = q1
    m.a[2][0] = c;                // q1' = q + q1 + c
    m.a[2][1] = 1;
    m.a[2][2] = 1;

    m.a[3][5] = 1;                // q'^2 = q1^2

    m.a[4][2] = c;                // q' q1' = q q1 + q1^2 + c q1
    m.a[4][4] = 1;
    m.a[4][5] = 1;

    m.a[5][0] = c2;               // q1'^2
    m.a[5][1] = (2 * c) % kMod;
    m.a[5][2] = (2 * c) % kMod;
    m.a[5][3] = 1;
    m.a[5][4] = 2;
    m.a[5][5] = 1;

    m.a[6][9] = 1;                // q'^3 = q1^3

    m.a[7][5] = c;                // q'^2 q1'
    m.a[7][8] = 1;
    m.a[7][9] = 1;

    m.a[8][2] = c2;               // q' q1'^2
    m.a[8][4] = (2 * c) % kMod;
    m.a[8][5] = (2 * c) % kMod;
    m.a[8][7] = 1;
    m.a[8][8] = 2;
    m.a[8][9] = 1;

    m.a[9][0] = c3;               // q1'^3
    m.a[9][1] = (3 * c2) % kMod;
    m.a[9][2] = (3 * c2) % kMod;
    m.a[9][3] = (3 * c) % kMod;
    m.a[9][4] = (6 * c) % kMod;
    m.a[9][5] = (3 * c) % kMod;
    m.a[9][6] = 1;
    m.a[9][7] = 3;
    m.a[9][8] = 3;
    m.a[9][9] = 1;

    // sum' = sum + S(F_k), and S(F_k)=a3*q^3 + a2(rem)*q^2 + a1(rem)*q + a0(rem)
    m.a[10][0] = coef.a0[static_cast<std::size_t>(rem)];
    m.a[10][1] = coef.a1[static_cast<std::size_t>(rem)];
    m.a[10][3] = coef.a2[static_cast<std::size_t>(rem)];
    m.a[10][6] = coef.a3;
    m.a[10][10] = 1;

    return m;
}

u64 solve(const i64 k_max, const Coefs& coef) {
    if (k_max < 2) {
        return 0;
    }

    std::array<int, 24> fib_mod24{};
    fib_mod24[0] = 0;
    fib_mod24[1] = 1;
    for (int i = 2; i < 24; ++i) {
        fib_mod24[static_cast<std::size_t>(i)] =
            (fib_mod24[static_cast<std::size_t>(i - 1)] + fib_mod24[static_cast<std::size_t>(i - 2)]) % 24;
    }

    std::array<int, 24> carry{};
    for (int i = 0; i < 24; ++i) {
        carry[static_cast<std::size_t>(i)] =
            (fib_mod24[static_cast<std::size_t>(i)] + fib_mod24[static_cast<std::size_t>((i + 1) % 24)]) / 24;
    }

    std::array<Mat, 24> step{};
    for (int phase = 0; phase < 24; ++phase) {
        step[static_cast<std::size_t>(phase)] =
            step_matrix(carry[static_cast<std::size_t>(phase)],
                        fib_mod24[static_cast<std::size_t>(phase)],
                        coef);
    }

    // We apply steps for k = 2..k_max, so number of steps is k_max-1.
    const i64 steps = k_max - 1;
    const int start_phase = 2 % 24;

    Mat block = mat_identity();
    for (int t = 0; t < 24; ++t) {
        const int phase = (start_phase + t) % 24;
        block = mat_mul(step[static_cast<std::size_t>(phase)], block);
    }

    std::array<u64, kDim> state{};
    // F_2 = 1 = 24*0 + 1, F_3 = 2 = 24*0 + 2 => q_2=0, q_3=0.
    state[0] = 1;  // constant term
    state[1] = 0;
    state[2] = 0;
    state[3] = 0;
    state[4] = 0;
    state[5] = 0;
    state[6] = 0;
    state[7] = 0;
    state[8] = 0;
    state[9] = 0;
    state[10] = 0;  // accumulated sum

    const i64 full_blocks = steps / 24;
    const int rem_steps = static_cast<int>(steps % 24);

    if (full_blocks > 0) {
        const Mat block_pow = mat_pow(block, full_blocks);
        state = mat_vec_mul(block_pow, state);
    }

    for (int t = 0; t < rem_steps; ++t) {
        const int phase = (start_phase + t) % 24;
        state = mat_vec_mul(step[static_cast<std::size_t>(phase)], state);
    }

    return state[10] % kMod;
}

bool run_checkpoints(const Coefs& coef) {
    if (s_mod(10, coef) != 1972ULL) {
        std::cerr << "Checkpoint failed: S(10)\n";
        return false;
    }
    if (s_mod(10000, coef) != (2024258331114ULL % kMod)) {
        std::cerr << "Checkpoint failed: S(10000) mod 1e9\n";
        return false;
    }
    // Exact S(10000) is given; recompute directly in 128-bit for full check.
    {
        const i64 n = 10000;
        const i64 q = n / 24;
        const int rem = static_cast<int>(n % 24);
        __int128 total = 0;
        for (int a = 1; a <= 24; ++a) {
            const i64 ca = q + (a <= rem ? 1 : 0);
            for (int b = 1; b <= 24; ++b) {
                const i64 cb = q + (b <= rem ? 1 : 0);
                for (int c = 1; c <= 24; ++c) {
                    const i64 cc = q + (c <= rem ? 1 : 0);
                    const i64 m = gcd4(
                        24,
                        36 + 6 * a,
                        14 + 6 * a + 2 * b,
                        1 + a + b + c);
                    total += static_cast<__int128>(ca) * cb * cc * m;
                }
            }
        }
        const i64 exact = static_cast<i64>(total);
        if (exact != 2024258331114LL) {
            std::cerr << "Checkpoint failed: exact S(10000)\n";
            return false;
        }
    }
    return true;
}

}  // namespace

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

    const Coefs coef = build_coefficients();
    if (options.run_checkpoints && !run_checkpoints(coef)) {
        return 2;
    }

    const u64 answer = solve(options.k_max, coef);
    std::cout << std::setfill('0') << std::setw(9) << answer << '\n';
    return 0;
}

Python

import math

MOD = 1000000000
DIM = 11

def gcd4(a, b, c, d):
    g2 = math.gcd
    return g2(g2(a, b), g2(c, d))

def build_coefficients():
    a3 = 0
    a2 = [0] * 24
    a1 = [0] * 24
    a0 = [0] * 24
    for ra in range(1, 25):
        for rb in range(1, 25):
            for rc in range(1, 25):
                m = gcd4(24, 36 + 6 * ra, 14 + 6 * ra + 2 * rb, 1 + ra + rb + rc)
                a3 = (a3 + m) % MOD
                for rem in range(24):
                    ea = 1 if ra <= rem else 0
                    eb = 1 if rb <= rem else 0
                    ec = 1 if rc <= rem else 0
                    u = m % MOD
                    a2[rem] = (a2[rem] + u * (ea + eb + ec)) % MOD
                    a1[rem] = (a1[rem] + u * (ea * eb + ea * ec + eb * ec)) % MOD
                    a0[rem] = (a0[rem] + u * (ea * eb * ec)) % MOD
    return a3, a2, a1, a0

def mat_identity():
    return [[1 if i == j else 0 for j in range(DIM)] for i in range(DIM)]

def mat_mul(x, y):
    z = [[0] * DIM for _ in range(DIM)]
    for i in range(DIM):
        for k in range(DIM):
            if x[i][k] == 0:
                continue
            for j in range(DIM):
                if y[k][j] == 0:
                    continue
                z[i][j] = (z[i][j] + x[i][k] * y[k][j]) % MOD
    return z

def mat_pow(base, exp):
    res = mat_identity()
    e = exp
    while e > 0:
        if e & 1:
            res = mat_mul(base, res)
        base = mat_mul(base, base)
        e >>= 1
    return res

def mat_vec_mul(m, v):
    out = [0] * DIM
    for i in range(DIM):
        s = 0
        for j in range(DIM):
            if m[i][j] != 0:
                s = (s + m[i][j] * v[j]) % MOD
        out[i] = s
    return out

def step_matrix(carry, rem, a3, a2, a1, a0):
    m = [[0] * DIM for _ in range(DIM)]
    c = carry % MOD
    c2 = (c * c) % MOD
    c3 = (c2 * c) % MOD
    
    m[0][0] = 1
    m[1][2] = 1
    m[2][0] = c
    m[2][1] = 1
    m[2][2] = 1
    
    m[3][5] = 1
    m[4][2] = c
    m[4][4] = 1
    m[4][5] = 1
    
    m[5][0] = c2
    m[5][1] = (2 * c) % MOD
    m[5][2] = (2 * c) % MOD
    m[5][3] = 1
    m[5][4] = 2
    m[5][5] = 1
    
    m[6][9] = 1
    m[7][5] = c
    m[7][8] = 1
    m[7][9] = 1
    
    m[8][2] = c2
    m[8][4] = (2 * c) % MOD
    m[8][5] = (2 * c) % MOD
    m[8][7] = 1
    m[8][8] = 2
    m[8][9] = 1
    
    m[9][0] = c3
    m[9][1] = (3 * c2) % MOD
    m[9][2] = (3 * c2) % MOD
    m[9][3] = (3 * c) % MOD
    m[9][4] = (6 * c) % MOD
    m[9][5] = (3 * c) % MOD
    m[9][6] = 1
    m[9][7] = 3
    m[9][8] = 3
    m[9][9] = 1
    
    m[10][0] = a0[rem]
    m[10][1] = a1[rem]
    m[10][3] = a2[rem]
    m[10][6] = a3
    m[10][10] = 1
    
    return m

def solve():
    k_max = 1234567890123
    a3, a2, a1, a0 = build_coefficients()
    
    fib_mod24 = [0] * 24
    fib_mod24[0] = 0
    fib_mod24[1] = 1
    for i in range(2, 24):
        fib_mod24[i] = (fib_mod24[i-1] + fib_mod24[i-2]) % 24
        
    carry = [0] * 24
    for i in range(24):
        carry[i] = (fib_mod24[i] + fib_mod24[(i+1)%24]) // 24
        
    step = []
    for phase in range(24):
        step.append(step_matrix(carry[phase], fib_mod24[phase], a3, a2, a1, a0))
        
    steps = k_max - 1
    start_phase = 2 % 24
    
    block = mat_identity()
    for t in range(24):
        phase = (start_phase + t) % 24
        block = mat_mul(step[phase], block)
        
    state = [0] * DIM
    state[0] = 1
    
    full_blocks = steps // 24
    rem_steps = steps % 24
    
    if full_blocks > 0:
        block_pow = mat_pow(block, full_blocks)
        state = mat_vec_mul(block_pow, state)
        
    for t in range(rem_steps):
        phase = (start_phase + t) % 24
        state = mat_vec_mul(step[phase], state)
        
    return "{:09d}".format(state[10] % MOD)

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

Java

public class Euler402 {
    private static final long MOD = 1000000000L;
    private static final int DIM = 11;

    private static long gcd(long x, long y) {
        while (y != 0) {
            long t = x % y;
            x = y;
            y = t;
        }
        return x < 0 ? -x : x;
    }

    private static long gcd4(long a, long b, long c, long d) {
        return gcd(gcd(a, b), gcd(c, d));
    }

    private static class Coefs {
        long a3 = 0;
        long[] a2 = new long[24];
        long[] a1 = new long[24];
        long[] a0 = new long[24];
    }

    private static Coefs buildCoefficients() {
        Coefs out = new Coefs();
        for (int ra = 1; ra <= 24; ++ra) {
            for (int rb = 1; rb <= 24; ++rb) {
                for (int rc = 1; rc <= 24; ++rc) {
                    long m = gcd4(24, 36 + 6L * ra, 14 + 6L * ra + 2L * rb, 1L + ra + rb + rc);
                    out.a3 = (out.a3 + m) % MOD;
                    for (int rem = 0; rem < 24; ++rem) {
                        int ea = (ra <= rem ? 1 : 0);
                        int eb = (rb <= rem ? 1 : 0);
                        int ec = (rc <= rem ? 1 : 0);
                        long u = m % MOD;
                        out.a2[rem] = (out.a2[rem] + u * (ea + eb + ec)) % MOD;
                        out.a1[rem] = (out.a1[rem] + u * (ea * eb + ea * ec + eb * ec)) % MOD;
                        out.a0[rem] = (out.a0[rem] + u * (ea * eb * ec)) % MOD;
                    }
                }
            }
        }
        return out;
    }

    private static class Mat {
        long[][] a = new long[DIM][DIM];
    }

    private static Mat matIdentity() {
        Mat m = new Mat();
        for (int i = 0; i < DIM; ++i) {
            m.a[i][i] = 1;
        }
        return m;
    }

    private static Mat matMul(Mat x, Mat y) {
        Mat z = new Mat();
        for (int i = 0; i < DIM; ++i) {
            for (int k = 0; k < DIM; ++k) {
                if (x.a[i][k] == 0)
                    continue;
                for (int j = 0; j < DIM; ++j) {
                    if (y.a[k][j] == 0)
                        continue;
                    z.a[i][j] = (z.a[i][j] + x.a[i][k] * y.a[k][j]) % MOD;
                }
            }
        }
        return z;
    }

    private static Mat matPow(Mat base, long exp) {
        Mat res = matIdentity();
        long e = exp;
        while (e > 0) {
            if ((e & 1) != 0) {
                res = matMul(base, res);
            }
            base = matMul(base, base);
            e >>= 1;
        }
        return res;
    }

    private static long[] matVecMul(Mat m, long[] v) {
        long[] out = new long[DIM];
        for (int i = 0; i < DIM; ++i) {
            long s = 0;
            for (int j = 0; j < DIM; ++j) {
                if (m.a[i][j] == 0)
                    continue;
                s = (s + m.a[i][j] * v[j]) % MOD;
            }
            out[i] = s;
        }
        return out;
    }

    private static Mat stepMatrix(int carry, int rem, Coefs coef) {
        Mat m = new Mat();
        long c = carry % MOD;
        long c2 = (c * c) % MOD;
        long c3 = (c2 * c) % MOD;

        m.a[0][0] = 1;
        m.a[1][2] = 1;
        m.a[2][0] = c;
        m.a[2][1] = 1;
        m.a[2][2] = 1;

        m.a[3][5] = 1;
        m.a[4][2] = c;
        m.a[4][4] = 1;
        m.a[4][5] = 1;

        m.a[5][0] = c2;
        m.a[5][1] = (2 * c) % MOD;
        m.a[5][2] = (2 * c) % MOD;
        m.a[5][3] = 1;
        m.a[5][4] = 2;
        m.a[5][5] = 1;

        m.a[6][9] = 1;
        m.a[7][5] = c;
        m.a[7][8] = 1;
        m.a[7][9] = 1;

        m.a[8][2] = c2;
        m.a[8][4] = (2 * c) % MOD;
        m.a[8][5] = (2 * c) % MOD;
        m.a[8][7] = 1;
        m.a[8][8] = 2;
        m.a[8][9] = 1;

        m.a[9][0] = c3;
        m.a[9][1] = (3 * c2) % MOD;
        m.a[9][2] = (3 * c2) % MOD;
        m.a[9][3] = (3 * c) % MOD;
        m.a[9][4] = (6 * c) % MOD;
        m.a[9][5] = (3 * c) % MOD;
        m.a[9][6] = 1;
        m.a[9][7] = 3;
        m.a[9][8] = 3;
        m.a[9][9] = 1;

        m.a[10][0] = coef.a0[rem];
        m.a[10][1] = coef.a1[rem];
        m.a[10][3] = coef.a2[rem];
        m.a[10][6] = coef.a3;
        m.a[10][10] = 1;

        return m;
    }

    public static String solve() {
        long kMax = 1234567890123L;
        Coefs coef = buildCoefficients();

        int[] fibMod24 = new int[24];
        fibMod24[0] = 0;
        fibMod24[1] = 1;
        for (int i = 2; i < 24; ++i) {
            fibMod24[i] = (fibMod24[i - 1] + fibMod24[i - 2]) % 24;
        }

        int[] carry = new int[24];
        for (int i = 0; i < 24; ++i) {
            carry[i] = (fibMod24[i] + fibMod24[(i + 1) % 24]) / 24;
        }

        Mat[] step = new Mat[24];
        for (int phase = 0; phase < 24; ++phase) {
            step[phase] = stepMatrix(carry[phase], fibMod24[phase], coef);
        }

        long steps = kMax - 1;
        int startPhase = 2 % 24;

        Mat block = matIdentity();
        for (int t = 0; t < 24; ++t) {
            int phase = (startPhase + t) % 24;
            block = matMul(step[phase], block);
        }

        long[] state = new long[DIM];
        state[0] = 1;

        long fullBlocks = steps / 24;
        int remSteps = (int) (steps % 24);

        if (fullBlocks > 0) {
            Mat blockPow = matPow(block, fullBlocks);
            state = matVecMul(blockPow, state);
        }

        for (int t = 0; t < remSteps; ++t) {
            int phase = (startPhase + t) % 24;
            state = matVecMul(step[phase], state);
        }

        return String.format("%09d", state[10] % MOD);
    }

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