Detailed mathematical approach

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\). Therefore the quotient \(P_{a,b,c}(n)/m\) is integer-valued for all \(n\) exactly when every coefficient in the binomial-basis expansion is divisible by \(m\). The maximal such divisor is the gcd of those coefficients:

$$\boxed{M(a,b,c)=\gcd\bigl(24,\ 36+6a,\ 14+6a+2b,\ 1+a+b+c\bigr).}$$

This explains the statement example immediately:

$$M(4,2,5)=\gcd(24,60,42,12)=6.$$

A tiny consistency check is

$$S(1)=M(1,1,1)=\gcd(24,42,22,4)=2.$$

3. Why \(S(N)\) becomes a cubic quasi-polynomial

Only residues modulo \(24\) matter, because \(\gcd(24,x)\) depends only on \(x \pmod{24}\). Write

$$N=24q+r,\qquad 0\le r<24.$$

For each residue \(t\in\{1,\dots,24\}\), the number of integers in \(\{1,\dots,N\}\) with that residue is

$$q+\varepsilon_t(r),\qquad \varepsilon_t(r)=\begin{cases}1,& t\le r,\\0,& t>r.\end{cases}$$

Hence

$$S(N)=\sum_{r_a=1}^{24}\sum_{r_b=1}^{24}\sum_{r_c=1}^{24} M(r_a,r_b,r_c)\bigl(q+\varepsilon_{r_a}\bigr)\bigl(q+\varepsilon_{r_b}\bigr)\bigl(q+\varepsilon_{r_c}\bigr).$$

Expanding the product yields

$$S(N)=C_3 q^3 + C_2(r) q^2 + C_1(r) q + C_0(r),$$

where \(C_3\) is constant and \(C_2,C_1,C_0\) depend only on the remainder \(r\). So \(S(N)\) is a degree-3 quasi-polynomial with period \(24\). The implementations precompute these four coefficient tables once by enumerating all \(24^3\) residue triples.

The built-in checkpoints are

$$S(10)=1972,\qquad S(10000)=2024258331114.$$

4. Evaluate the quasi-polynomial at Fibonacci arguments

Write each Fibonacci number as

$$F_k=24q_k+r_k,\qquad 0\le r_k<24.$$

Because Fibonacci numbers modulo \(24\) have Pisano period \(24\), the residue sequence \(r_k\) is periodic. From \(F_{k+2}=F_{k+1}+F_k\) we obtain

$$r_{k+2}\equiv r_{k+1}+r_k \pmod{24},$$

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

The carry \(\delta_k\) depends only on the residue pair \((r_k,r_{k+1})\), so it is periodic as well. This converts the huge-index problem into a fixed periodic affine recurrence for the quotients \(q_k\).

5. Linearize the recurrence with monomials up to degree \(3\)

Since \(S(F_k)\) is cubic in \(q_k\), it is enough to track all monomials in two consecutive quotient variables \(x=q_k\) and \(y=q_{k+1}\) up to degree \(3\):

$$1,\ x,\ y,\ x^2,\ xy,\ y^2,\ x^3,\ x^2y,\ xy^2,\ y^3,$$

plus one extra coordinate for the running total \(\sum_{j=2}^{k} S(F_j)\). Under the update

$$x'=y,\qquad y'=x+y+\delta_k,$$

every basis monomial becomes a linear combination of the same basis. For example,

$$x'y'=y(x+y+\delta_k),\qquad y'^2=(x+y+\delta_k)^2,\qquad y'^3=(x+y+\delta_k)^3.$$

Therefore one Fibonacci step is represented by an \(11\times 11\) matrix. Multiplying the \(24\) phase matrices of one full residue cycle gives a single block matrix, and binary exponentiation of that block reaches \(K=1234567890123\) in logarithmic time.

How the Code Works

The C++, Python, and Java implementations all use the same structure. They first precompute the quasi-polynomial coefficients of \(S(N)\) from the \(24^3\) residue kernel. They then build the \(24\) phase-dependent transition matrices determined by the periodic Fibonacci residues and carries. Finally they exponentiate the one-cycle block matrix, apply the remaining partial cycle, and read the accumulated-sum coordinate modulo \(10^9\).

Complexity Analysis

The residue precomputation is constant-size work, bounded by \(24^4\) elementary updates. The main stage is binary exponentiation of fixed \(11\times 11\) matrices, so the running time is \(O(\log K)\) and the memory usage is \(O(1)\) apart from fixed-size tables and matrices.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=402
2. Integer-valued polynomial: Wikipedia — Integer-valued polynomial
3. Binomial coefficient basis: Wikipedia — Binomial coefficient
4. Pisano periods: Wikipedia — Pisano period
5. Matrix exponentiation: Wikipedia — Matrix exponentiation
6. Quasi-polynomials: Wikipedia — Quasi-polynomial