Solution

Answer: 356019862

Let

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

We seek the largest integer $m$ such that $P(n)\equiv 0\pmod m$ for all integers $n$.

1. Mathematical analysis

Step 1: Characterizing $M(a,b,c)$

A standard fact about integer-valued polynomials is:

The gcd of all values of a polynomial over the integers equals the gcd of its forward differences at $0$.

Compute the forward differences of $P(n)$:

$$\Delta P(0)=1+a+b+c,$$

$$\Delta^2 P(0)=14+6a+2b,$$

$$\Delta^3 P(0)=36+6a,$$

$$\Delta^4 P(0)=24.$$

Therefore

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

So $M(a,b,c)$ always divides $24$.

Hence $M$ depends only on $a,b,c \pmod{24}$.

Step 2: Periodicity of $M$

Because only residues modulo $24$ matter,

$$M(a+24,b,c)=M(a,b,c)$$

(and similarly for $b,c$).

Thus $S(N)$ is a cubic quasipolynomial with period $24$.

Write

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

For each residue class $t\in{1,\dots,24}$, the count of integers

$\le N$ congruent to $t\pmod{24}$ is either $q$ or $q+1$.

Summing over residue classes gives

$$S(N)=Aq^3+B_r q^2+C_r q+D_r,$$

where the constants are determined from the $24^3$ residue cube.

The full-period total is

$$A=\sum_{1\le a,b,c\le 24} M(a,b,c)=27984.$$

This reproduces the checks:

$$S(10)=1972,$$

$$S(10000)=2024258331114.$$

Step 3: Fibonacci input

We need

$$\sum_{k=2}^{1234567890123} S(F_k)\pmod{10^9}.$$

Since Fibonacci numbers modulo $24$ have Pisano period $24$,

$$F_{k+24}\equiv F_k\pmod{24},$$

the coefficients $B_r,C_r,D_r$ repeat with period $24$.

Thus $S(F_k)$ becomes a linear combination of

$$1,\quad F_k,\quad F_k^2,\quad F_k^3$$

with coefficients periodic modulo $24$.

These sequences satisfy linear recurrences, so the whole sequence can be evaluated with matrix exponentiation in $O(\log n)$.

2. Python implementation

from math import gcd
import numpy as np

MOD = 10**9
BIGMOD = 288 * MOD

# ------------------------------------------------------------
# Compute M(a,b,c)
# ------------------------------------------------------------

def M(a, b, c):
    return gcd(
        gcd(
            gcd(1 + a + b + c,
                14 + 6*a + 2*b),
            36 + 6*a),
        24
    )

# ------------------------------------------------------------
# Precompute the 24x24x24 residue cube
# ------------------------------------------------------------

cube = np.zeros((25, 25, 25), dtype=object)

for a in range(1, 25):
    for b in range(1, 25):
        for c in range(1, 25):
            cube[a, b, c] = M(a, b, c)

# 3D prefix sums
pref = cube.cumsum(0).cumsum(1).cumsum(2)

# ------------------------------------------------------------
# S(N)
# ------------------------------------------------------------

def S(N):
    q, r = divmod(N, 24)

    return (
        q**3 * pref[24,24,24]
        + q*q * (
            pref[r,24,24]
            + pref[24,r,24]
            + pref[24,24,r]
        )
        + q * (
            pref[r,r,24]
            + pref[r,24,r]
            + pref[24,r,r]
        )
        + pref[r,r,r]
    )

# ------------------------------------------------------------
# Verify examples
# ------------------------------------------------------------

assert S(10) == 1972
assert S(10000) == 2024258331114

# ------------------------------------------------------------
# Fibonacci residues mod 24
# ------------------------------------------------------------

fib24 = [0, 1]
for _ in range(22):
    fib24.append((fib24[-1] + fib24[-2]) % 24)

# ------------------------------------------------------------
# Build polynomial coefficients
# We work with 288*S(F_n) to avoid fractions.
# ------------------------------------------------------------

import sympy as sp

q = sp.Symbol('q')
x = sp.Symbol('x')

# Recover cubic polynomials for each residue
poly_by_r = {}

for r in range(24):
    vals = [S(24*k + r) for k in range(4)]
    poly = sp.interpolate(
        [(k, vals[k]) for k in range(4)],
        q
    )
    poly_by_r[r] = sp.expand(poly)

coeffs = []

for r in range(24):
    expr = sp.expand(
        288 * poly_by_r[r].subs(q, (x-r)/24)
    )
    P = sp.Poly(expr, x)

    # const, x, x^2, x^3
    coeffs.append([
        int(P.nth(0)),
        int(P.nth(1)),
        int(P.nth(2)),
        int(P.nth(3)),
    ])

# ------------------------------------------------------------
# State vector:
#
# [F_n^3,
#  F_n^2 F_{n+1},
#  F_n F_{n+1}^2,
#  F_{n+1}^3,
#  F_n^2,
#  F_n F_{n+1},
#  F_{n+1}^2,
#  F_n,
#  F_{n+1},
#  1,
#  accumulated_sum]
# ------------------------------------------------------------

mons = [
    (3,0),(2,1),(1,2),(0,3),
    (2,0),(1,1),(0,2),
    (1,0),(0,1),(0,0)
]

idx = {m:i for i,m in enumerate(mons)}

T = np.zeros((10,10), dtype=object)

# Fibonacci transition:
# (x,y) -> (y, x+y)

from math import comb

for j, (a,b) in enumerate(mons):
    for k in range(b+1):
        coef = comb(b, k)
        mon = (k, a+b-k)
        i = idx[mon]
        T[i,j] += coef

# ------------------------------------------------------------
# Phase-dependent matrices
# ------------------------------------------------------------

def make_matrix(t):
    A = np.zeros((11,11), dtype=object)

    A[:10,:10] = T
    A[10,10] = 1

    c = coeffs[(t+1) % 24]

    # depends only on F_{n+1}
    A[10,3] = c[3]
    A[10,6] = c[2]
    A[10,8] = c[1]
    A[10,9] = c[0]

    return A % BIGMOD

phase = [make_matrix(t) for t in range(24)]

# ------------------------------------------------------------
# Multiply matrices
# ------------------------------------------------------------

def matmul(A, B):
    return (A.dot(B)) % BIGMOD

# 24-step block
B = np.identity(11, dtype=object)

for t in range(24):
    B = matmul(phase[t], B)

# ------------------------------------------------------------
# Fast exponentiation
# ------------------------------------------------------------

def mpow(M, n):
    R = np.identity(11, dtype=object)

    while n:
        if n & 1:
            R = matmul(M, R)

        M = matmul(M, M)
        n >>= 1

    return R

# ------------------------------------------------------------
# Initial state at n=0
# ------------------------------------------------------------

state0 = np.array([
    0,0,0,1,
    0,0,1,
    0,1,
    1,
    0
], dtype=object)

K = 1234567890123

blocks, rem = divmod(K, 24)

state = mpow(B, blocks).dot(state0) % BIGMOD

for t in range(rem):
    state = phase[t].dot(state) % BIGMOD

# state[10] = 288 * sum_{k=1}^K S(F_k)

answer = ((state[10] // 288) - 2) % MOD

print(answer)

3. Code walkthrough

Computing $M(a,b,c)$

We use the finite-difference characterization:

gcd(1+a+b+c,
    14+6a+2b,
    36+6a,
    24)

This immediately implies periodicity modulo $24$.

Computing $S(N)$

Since $M$ depends only on residues mod $24$, we precompute the entire $24^3$ cube once.

Then for

$$N=24q+r,$$

we count how many times each residue appears.

This reduces the whole problem to evaluating a cubic quasipolynomial.

Fibonacci acceleration

The dependence on $F_k$ is cubic, so we track all monomials up to degree $3$:

$$F_k^3,\quad F_k^2F_{k+1},\quad \dots$$

The Fibonacci transition

$$(F_k,F_{k+1})\to(F_{k+1},F_k+F_{k+1})$$

induces a linear transformation on these monomials.

We then use fast matrix exponentiation to jump directly to

$$k=1234567890123.$$

4. Final verification

The program reproduces the given checks exactly: