Project Euler

Integer-valued Polynomials

The polynomial n^4 + 4n^3 + 2n^2 + 5n is a multiple of 6 for every integer n, and 6 is the largest such integer. Define M(a, b, c) as the maximum m such that n^4 + an^3 + bn^2 + cn is a multiple of...

Source sync Apr 19, 2026

Problem #0402

Level Level 23

Solved By 504

Languages C++, Python

Answer 356019862

Length 454 words

modular_arithmeticlinear_algebranumber_theory

It can be shown that the polynomial $n^4 + 4n^3 + 2n^2 + 5n$ is a multiple of $6$ for every integer $n$. It can also be shown that $6$ is the largest integer satisfying this property.

Define $M(a, b, c)$ as the maximum $m$ such that $n^4 + an^3 + bn^2 + cn$ is a multiple of $m$ for all integers $n$. For example, $M(4, 2, 5) = 6$.

Also, define $S(N)$ as the sum of $M(a, b, c)$ for all $0 < a, b, c \leq N$.

We can verify that $S(10) = 1972$ and $S(10000) = 2024258331114$.

Let $F_k$ be the Fibonacci sequence: \begin {equation} \begin {cases} F_0 = 0 \\ F_1 = 1 \\ F_k = F_{k - 1} + F_{k - 2}, \quad \text {for } k \geq 2 \end {cases} \end {equation}

Find the last $9$ digits of $\sum S(F_k)$ for $2 \leq k \leq 1234567890123$.

Problem 402: Integer-valued Polynomials

Mathematical Foundation

Theorem 1 (Binomial basis characterization of $M(a,b,c)$ ). For integers $a, b, c$ ,

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

Proof. A polynomial $p(n) = \sum_{k=0}^{d} c_k \binom{n}{k}$ satisfies $m \mid p(n)$ for all $n \in \mathbb{Z}$ if and only if $m \mid c_k$ for every $k$ . This is because the binomial coefficients $\binom{n}{0}, \binom{n}{1}, \ldots$ form a $\mathbb{Z}$ -basis for the ring of integer-valued polynomials, and $\binom{n}{k} \in \mathbb{Z}$ for all $n \in \mathbb{Z}$ .

Lemma 1 (Newton forward-difference expansion). The monomials expand as:

n = \binom{n}{1}, \quad n^2 = 2\binom{n}{2} + \binom{n}{1}, \quad 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}.

Proof. By the Stirling number identity $n^k = \sum_{j=0}^{k} S(k,j) \cdot j! \cdot \binom{n}{j}$ , where $S(k,j)$ are Stirling numbers of the second kind. Direct computation gives the stated coefficients. $\square$

Applying the lemma, $p(n) = n^4 + an^3 + bn^2 + cn$ has binomial coefficients:

d_4 = 24, \quad d_3 = 36 + 6a, \quad d_2 = 14 + 6a + 2b, \quad d_1 = 1 + a + b + c.

The maximum universal divisor is $M(a,b,c) = \gcd(d_1, d_2, d_3, d_4)$ . $\square$

Theorem 2 (Euler totient decomposition). Since $M(a,b,c) \mid 24$ , we have

S(N) = \sum_{d \mid 24} \varphi(d) \cdot \#\{(a,b,c) \in [1,N]^3 : d \mid M(a,b,c)\}

where $\varphi$ is Euler’s totient function.

Proof. By the identity $n = \sum_{d \mid n} \varphi(d)$ , we have $M(a,b,c) = \sum_{d \mid M(a,b,c)} \varphi(d)$ . Summing over all triples and interchanging:

S(N) = \sum_{(a,b,c)} \sum_{d \mid M(a,b,c)} \varphi(d) = \sum_{d \mid 24} \varphi(d) \sum_{\substack{(a,b,c) \\ d \mid M(a,b,c)}} 1. \quad \square

Lemma 2 (Piecewise-cubic structure). For each residue $s = N \bmod 24$ , the function $S(N)$ is a cubic polynomial in $N$ :

288 \cdot S(N) = A_s N^3 + B_s N^2 + C_s N + D_s, \quad N \equiv s \pmod{24},

where $A_s = 583$ for all $s$ .

Proof. The condition $d \mid M(a,b,c)$ reduces to three congruences modulo $d$ on $(a, b, c)$ . For each valid residue class, the count of triples in $[1,N]^3$ is a product of three terms of the form $\lfloor (N - r)/d \rfloor$ , each piecewise-linear in $N$ with period $d$ . The product is piecewise-cubic with period $\operatorname{lcm}$ of all $d \mid 24$ , which is $24$ . The factor $288 = \operatorname{lcm}$ of the coefficient denominators ensures integrality. $\square$

Theorem 3 (Matrix exponentiation for Fibonacci power sums). For each residue class $r \bmod 24$ , the subsequence $(F_{24i+r})$ satisfies a two-term linear recurrence. The sums $\sum F_{24i+r}^j$ for $j = 0, 1, 2, 3$ are computed via a $14$ -dimensional linear recurrence solved by matrix exponentiation.

Proof. The Pisano period $\pi(24) = 24$ , so $F_k \bmod 24$ depends only on $k \bmod 24$ . The matrix $A^{24}$ where $A = \bigl(\begin{smallmatrix}1&1\\1&0\end{smallmatrix}\bigr)$ advances the Fibonacci sequence by 24 steps: $(F_{k+24}, F_{k+25}) = M_{24}(F_k, F_{k+1})$ .

Define the state vector of dimension 14:

\mathbf{x} = (u^3, u^2v, uv^2, v^3, u^2, uv, v^2, u, v, 1, \Sigma u^3, \Sigma u^2, \Sigma u, \Sigma 1)

where $u = F_{24i+r}$ , $v = F_{24i+r+1}$ . The transition $(u,v) \to (F_{23}u + F_{24}v,\; F_{24}u + F_{25}v)$ is linear, and all monomials up to degree 3 in $(u,v)$ transform linearly. The accumulators track running sums. This yields a $14 \times 14$ transition matrix $T$ , and $T^{N}$ is computed via repeated squaring in $O(14^3 \log N)$ operations. $\square$

Editorial

The polynomial n^4 + 4n^3 + 2n^2 + 5n is a multiple of 6 for every integer n. Define M(a,b,c) as the maximum m such that n^4 + an^3 + bn^2 + cn is a multiple of m for all integers n. So M(4,2,5) = 6. Define S(N) = sum of M(a,b,c) for 0 < a,b,c <= N. Given: S(10) = 1972, S(10000) = 2024258331114. With Fibonacci F_0=0, F_1=1, F_k = F_{k-1}+F_{k-2}: Find last 9 digits of sum_{k=2}^{1234567890123} S(F_k). We precompute the 24 cubic polynomials (A_s, B_s, C_s, D_s). We then by fitting S(N) at 4 sample points per residue class s mod 24. Finally, iterate over each residue class r = 0..23 of k mod 24.

Pseudocode

Precompute the 24 cubic polynomials (A_s, B_s, C_s, D_s)
by fitting S(N) at 4 sample points per residue class s mod 24
For each residue class r = 0..23 of k mod 24:
Determine which polynomial to use: s = F_r mod 24 (Pisano period)
Build 14x14 transition matrix T for subsequence F_{24i+r}
Determine number of terms in subsequence with index k in [2, K], k ≡ r (mod 24)
Matrix exponentiation: T^{n_terms - 1} applied to initial state
Extract accumulated sums: Σu^3, Σu^2, Σu, Σ1
Divide by 288 and reduce modulo MOD

Complexity Analysis

Time: $O(14^3 \cdot 24 \cdot \log K) \approx O(10^6)$ for $K = 1.23 \times 10^{12}$ . The dominant cost is 24 matrix exponentiations of $14 \times 14$ matrices with exponents $\sim K/24$ .
Space: $O(14^2) = O(1)$ . Only a constant number of matrices and vectors.

Answer

\boxed{356019862}

C++ project_euler/problem_402/solution.cpp

#include <bits/stdc++.h>
using namespace std;

/*
 * Project Euler Problem 402: Integer-valued Polynomials
 *
 * M(a,b,c) = gcd(24, 36+6a, 14+6a+2b, 1+a+b+c)
 * S(N) = sum_{0<a,b,c<=N} M(a,b,c)
 * Find last 9 digits of sum_{k=2}^{1234567890123} S(F_k)
 *
 * Answer: 356019862
 */

typedef long long ll;
typedef __int128 lll;
typedef vector<vector<ll>> Matrix;

const ll MOD_FINAL = 1000000000LL;       // 10^9
const ll WORK_MOD  = 288LL * MOD_FINAL;  // 288 * 10^9

// =====================================================================
// Safe modular reduction for __int128
// =====================================================================

inline ll mmod(lll x) {
    ll r = (ll)(x % WORK_MOD);
    return r < 0 ? r + WORK_MOD : r;
}

// =====================================================================
// Matrix operations mod WORK_MOD
// =====================================================================

Matrix mat_mul(const Matrix& A, const Matrix& B) {
    int n = A.size(), m = B[0].size(), k = B.size();
    Matrix C(n, vector<ll>(m, 0));
    for (int i = 0; i < n; i++)
        for (int j = 0; j < m; j++) {
            lll s = 0;
            for (int l = 0; l < k; l++)
                s += (lll)A[i][l] * B[l][j];
            C[i][j] = mmod(s);
        }
    return C;
}

Matrix mat_pow(Matrix A, ll p) {
    int n = A.size();
    Matrix result(n, vector<ll>(n, 0));
    for (int i = 0; i < n; i++) result[i][i] = 1;
    while (p > 0) {
        if (p & 1) result = mat_mul(result, A);
        A = mat_mul(A, A);
        p >>= 1;
    }
    return result;
}

vector<ll> mat_vec(const Matrix& A, const vector<ll>& v) {
    int n = A.size();
    vector<ll> result(n, 0);
    for (int i = 0; i < n; i++) {
        lll s = 0;
        for (int j = 0; j < (int)v.size(); j++)
            s += (lll)A[i][j] * v[j];
        result[i] = mmod(s);
    }
    return result;
}

// =====================================================================
// S(N) analytical computation (for verification)
// =====================================================================

int euler_totient(int n) {
    int result = n, temp = n;
    for (int p = 2; p * p <= temp; p++) {
        if (temp % p == 0) {
            while (temp % p == 0) temp /= p;
            result -= result / p;
        }
    }
    if (temp > 1) result -= result / temp;
    return result;
}

ll S_analytical(ll N) {
    int divs[] = {1, 2, 3, 4, 6, 8, 12, 24};
    ll total = 0;
    for (int d : divs) {
        int phi_d = euler_totient(d);
        ll cnt = 0;
        for (int ar = 0; ar < d; ar++) {
            if ((36 + 6 * ar) % d != 0) continue;
            int first_a = (ar >= 1) ? ar : d;
            if (first_a > N) continue;
            ll cnt_a = (N - first_a) / d + 1;
            for (int br = 0; br < d; br++) {
                if ((14 + 6 * ar + 2 * br) % d != 0) continue;
                int first_b = (br >= 1) ? br : d;
                if (first_b > N) continue;
                ll cnt_b = (N - first_b) / d + 1;
                int cr = ((-(1 + ar + br)) % d + d) % d;
                int first_c = (cr >= 1) ? cr : d;
                if (first_c > N) continue;
                ll cnt_c = (N - first_c) / d + 1;
                cnt += cnt_a * cnt_b * cnt_c;
            }
        }
        total += (ll)phi_d * cnt;
    }
    return total;
}

// =====================================================================
// Piecewise-cubic polynomial coefficients for 288 * S(N)
// 288*S(N) = A*N^3 + B*N^2 + C*N + D when N ≡ s (mod 24)
// =====================================================================

struct Poly { ll D, C, B, A; };

Poly POLYS[24] = {
    {0, 0, 0, 583},           // s=0
    {257, -75, -189, 583},    // s=1
    {520, 84, 102, 583},      // s=2
    {1215, -27, -51, 583},    // s=3
    {256, -64, -136, 583},    // s=4
    {2313, 229, 11, 583},     // s=5
    {360, -108, 14, 583},     // s=6
    {1063, -139, -139, 583},  // s=7
    {512, 128, 64, 583},      // s=8
    {225, 21, -125, 583},     // s=9
    {-1464, -140, -122, 583}, // s=10
    {4487, 245, 85, 583},     // s=11
    {4032, 0, 0, 583},        // s=12
    {4289, -75, -189, 583},   // s=13
    {1672, 84, 102, 583},     // s=14
    {1791, -27, -51, 583},    // s=15
    {832, -64, -136, 583},    // s=16
    {2889, 229, 11, 583},     // s=17
    {360, -108, 14, 583},     // s=18
    {1639, -139, -139, 583},  // s=19
    {1088, 128, 64, 583},     // s=20
    {801, 21, -125, 583},     // s=21
    {-1464, -140, -122, 583}, // s=22
    {455, 245, 85, 583}       // s=23
};

// Fibonacci mod 24 table (Pisano period = 24)
int fib_mod24[25];
void init_fib_mod24() {
    fib_mod24[0] = 0; fib_mod24[1] = 1;
    for (int i = 2; i < 25; i++)
        fib_mod24[i] = (fib_mod24[i-1] + fib_mod24[i-2]) % 24;
}

// =====================================================================
// Build 14x14 transition matrix
// State: [u^3, u^2v, uv^2, v^3, u^2, uv, v^2, u, v, 1,
//         sum_u^3, sum_u^2, sum_u, sum_1]
// Transition: u' = a*u + b*v, v' = c*u + d*v
// =====================================================================

Matrix build_transition(ll a, ll b, ll c, ll d) {
    Matrix T(14, vector<ll>(14, 0));

    // Degree-3 monomials
    T[0][0] = mmod((lll)a*a*a);
    T[0][1] = mmod((lll)3*a*a*b);
    T[0][2] = mmod((lll)3*a*b*b);
    T[0][3] = mmod((lll)b*b*b);

    T[1][0] = mmod((lll)a*a*c);
    T[1][1] = mmod((lll)a*a*d + (lll)2*a*b*c);
    T[1][2] = mmod((lll)2*a*b*d + (lll)b*b*c);
    T[1][3] = mmod((lll)b*b*d);

    T[2][0] = mmod((lll)a*c*c);
    T[2][1] = mmod((lll)2*a*c*d + (lll)b*c*c);
    T[2][2] = mmod((lll)a*d*d + (lll)2*b*c*d);
    T[2][3] = mmod((lll)b*d*d);

    T[3][0] = mmod((lll)c*c*c);
    T[3][1] = mmod((lll)3*c*c*d);
    T[3][2] = mmod((lll)3*c*d*d);
    T[3][3] = mmod((lll)d*d*d);

    // Degree-2 monomials
    T[4][4] = mmod((lll)a*a); T[4][5] = mmod((lll)2*a*b); T[4][6] = mmod((lll)b*b);
    T[5][4] = mmod((lll)a*c); T[5][5] = mmod((lll)a*d+(lll)b*c); T[5][6] = mmod((lll)b*d);
    T[6][4] = mmod((lll)c*c); T[6][5] = mmod((lll)2*c*d); T[6][6] = mmod((lll)d*d);

    // Degree-1
    T[7][7] = a; T[7][8] = b;
    T[8][7] = c; T[8][8] = d;

    // Constant
    T[9][9] = 1;

    // Accumulators: Su3' = Su3 + u'^3, etc.
    for (int j = 0; j < 4; j++) T[10][j] = T[0][j];
    T[10][10] = 1;

    for (int j = 4; j < 7; j++) T[11][j] = T[4][j];
    T[11][11] = 1;

    T[12][7] = T[7][7]; T[12][8] = T[7][8];
    T[12][12] = 1;

    T[13][9] = 1; T[13][13] = 1;

    return T;
}

// =====================================================================
// Main solver
// =====================================================================

ll solve(ll K) {
    init_fib_mod24();

    // Fibonacci constants: F_23=28657, F_24=46368, F_25=75025
    // Transition: u' = F23*u + F24*v, v' = F24*u + F25*v
    const ll F23 = 28657, F24 = 46368, F25 = 75025;

    Matrix T = build_transition(F23, F24, F24, F25);

    // Small Fibonacci values for initial conditions
    ll fib_small[30];
    fib_small[0] = 0; fib_small[1] = 1;
    for (int i = 2; i < 30; i++)
        fib_small[i] = fib_small[i-1] + fib_small[i-2];

    ll total_288 = 0;

    for (int r = 0; r < 24; r++) {
        int s = fib_mod24[r];
        Poly& poly = POLYS[s];

        // Polynomial coefficients mod WORK_MOD (handle negative values)
        ll D_c = mmod(poly.D);
        ll C_c = mmod(poly.C);
        ll B_c = mmod(poly.B);
        ll A_c = mmod(poly.A);

        // Index range for k = 24*i + r, with k in [2, K]
        ll i_min = (r >= 2) ? 0 : 1;

        // IMPORTANT: guard against negative (K - r) before division
        if (K < r) continue;
        ll i_max = (K - r) / 24;

        if (i_max < i_min) continue;
        ll n_terms = i_max - i_min + 1;

        // Initial Fibonacci values at k_start = 24*i_min + r
        ll k_start = 24 * i_min + r;
        ll u = fib_small[k_start] % WORK_MOD;
        ll v = fib_small[k_start + 1] % WORK_MOD;

        // Build initial state (with first term already in accumulators)
        vector<ll> state(14);
        state[0]  = mmod((lll)u*u*u);
        state[1]  = mmod((lll)u*u*v);
        state[2]  = mmod((lll)u*v*v);
        state[3]  = mmod((lll)v*v*v);
        state[4]  = mmod((lll)u*u);
        state[5]  = mmod((lll)u*v);
        state[6]  = mmod((lll)v*v);
        state[7]  = u;
        state[8]  = v;
        state[9]  = 1;
        state[10] = state[0]; // Su3 initialized with first term
        state[11] = state[4]; // Su2
        state[12] = u;        // Su
        state[13] = 1;        // count

        if (n_terms > 1) {
            Matrix T_pow = mat_pow(T, n_terms - 1);
            state = mat_vec(T_pow, state);
        }

        ll Su3 = state[10], Su2 = state[11], Su1 = state[12], S1 = state[13];

        lll contrib = (lll)A_c * Su3 + (lll)B_c * Su2
                    + (lll)C_c * Su1 + (lll)D_c * S1;
        total_288 = (total_288 + mmod(contrib)) % WORK_MOD;
    }

    return total_288 / 288;
}

// =====================================================================
// Verification
// =====================================================================

void verify() {
    printf("Verification:\n");

    // M(4,2,5) = 6
    auto M = [](int a, int b, int c) -> int {
        return __gcd(__gcd(24, 36+6*a), __gcd(14+6*a+2*b, 1+a+b+c));
    };
    assert(M(4,2,5) == 6);
    printf("  [PASS] M(4,2,5) = 6\n");

    // S(10) = 1972
    assert(S_analytical(10) == 1972);
    printf("  [PASS] S(10) = 1972\n");

    // S(10000) = 2024258331114
    assert(S_analytical(10000) == 2024258331114LL);
    printf("  [PASS] S(10000) = 2024258331114\n");

    // Cross-validate K=20 with brute force
    ll fib[21]; fib[0] = 0; fib[1] = 1;
    for (int i = 2; i <= 20; i++) fib[i] = fib[i-1] + fib[i-2];
    ll bf_sum = 0;
    for (int k = 2; k <= 20; k++)
        bf_sum = (bf_sum + S_analytical(fib[k]) % MOD_FINAL) % MOD_FINAL;
    ll mat_sum = solve(20);
    printf("  K=20 brute-force: %lld, matrix: %lld\n", bf_sum, mat_sum);
    assert(mat_sum == bf_sum);
    printf("  [PASS] K=20 cross-validation\n");
}

// =====================================================================
// Main
// =====================================================================

int main() {
    printf("Project Euler Problem 402: Integer-valued Polynomials\n");
    printf("=====================================================\n\n");

    verify();
    printf("\n");

    auto t0 = chrono::steady_clock::now();
    ll answer = solve(1234567890123LL);
    auto t1 = chrono::steady_clock::now();
    double elapsed = chrono::duration<double>(t1 - t0).count();

    printf("Answer: %lld\n", answer);
    printf("Time:   %.3f s\n", elapsed);

    return 0;
}

Python project_euler/problem_402/solution.py

"""
Project Euler Problem 402: Integer-valued Polynomials

The polynomial n^4 + 4n^3 + 2n^2 + 5n is a multiple of 6 for every integer n.
Define M(a,b,c) as the maximum m such that n^4 + an^3 + bn^2 + cn is a multiple of m
for all integers n. So M(4,2,5) = 6.

Define S(N) = sum of M(a,b,c) for 0 < a,b,c <= N.
Given: S(10) = 1972, S(10000) = 2024258331114.

With Fibonacci F_0=0, F_1=1, F_k = F_{k-1}+F_{k-2}:
Find last 9 digits of sum_{k=2}^{1234567890123} S(F_k).

Answer: 356019862

Key insight:
  Using Newton's forward difference (binomial coefficient) basis:
    n^4 = 24*C(n,4) + 36*C(n,3) + 14*C(n,2) + C(n,1)
    n^3 = 6*C(n,3) + 6*C(n,2) + C(n,1)
    n^2 = 2*C(n,2) + C(n,1)
    n   = C(n,1)

  So p(n) = n^4 + a*n^3 + b*n^2 + c*n
          = 24*C(n,4) + (36+6a)*C(n,3) + (14+6a+2b)*C(n,2) + (1+a+b+c)*C(n,1)

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

  Using Euler's totient decomposition and modular counting, S(N) becomes
  piecewise-cubic in N with period 24 (i.e., for each residue class N mod 24,
  S(N) is a degree-3 polynomial in N).

  The Fibonacci sum is computed by grouping k by k mod 24 (the Pisano period of 24),
  then using 14x14 matrix exponentiation to track power sums of F_k along each
  arithmetic subsequence.
"""

import time
from math import gcd
from fractions import Fraction

# ---------------------------------------------------------------------------
# Part 1: M(a,b,c) and S(N) via direct enumeration (for verification)
# ---------------------------------------------------------------------------

DIVISORS_24 = [1, 2, 3, 4, 6, 8, 12, 24]

def euler_totient(n):
    """Compute Euler's totient function phi(n)."""
    result = n
    p = 2
    temp = n
    while p * p <= temp:
        if temp % p == 0:
            while temp % p == 0:
                temp //= p
            result -= result // p
        p += 1
    if temp > 1:
        result -= result // temp
    return result

def M_value(a, b, c):
    """
    M(a,b,c) = max m such that m | (n^4 + a*n^3 + b*n^2 + c*n) for all n.

    Using the binomial basis representation:
    M(a,b,c) = gcd(24, 36+6a, 14+6a+2b, 1+a+b+c)
    """
    return gcd(gcd(24, 36 + 6 * a), gcd(14 + 6 * a + 2 * b, 1 + a + b + c))

def S_brute(N):
    """Compute S(N) by brute-force triple loop. O(N^3)."""
    return sum(M_value(a, b, c)
               for a in range(1, N + 1)
               for b in range(1, N + 1)
               for c in range(1, N + 1))

def S_analytical(N):
    """
    Compute S(N) in O(1) using Euler-totient decomposition.

    S(N) = sum_{d|24} phi(d) * #{(a,b,c) in [1,N]^3 : d | gcd(24, 36+6a, 14+6a+2b, 1+a+b+c)}

    For each d | 24, we enumerate O(d^2) residue classes and count via floor division.
    Since max(d) = 24, this is O(sum d^2 for d|24) = O(1).
    """
    total = 0
    for d in DIVISORS_24:
        phi_d = euler_totient(d)
        cnt = 0
        for ar in range(d):
            if (36 + 6 * ar) % d != 0:
                continue
            first_a = ar if ar >= 1 else d
            if first_a > N:
                continue
            cnt_a = (N - first_a) // d + 1
            for br in range(d):
                if (14 + 6 * ar + 2 * br) % d != 0:
                    continue
                first_b = br if br >= 1 else d
                if first_b > N:
                    continue
                cnt_b = (N - first_b) // d + 1
                cr = (-(1 + ar + br)) % d
                first_c = cr if cr >= 1 else d
                if first_c > N:
                    continue
                cnt_c = (N - first_c) // d + 1
                cnt += cnt_a * cnt_b * cnt_c
        total += phi_d * cnt
    return total

# ---------------------------------------------------------------------------
# Part 2: Piecewise-cubic representation of S(N)
# ---------------------------------------------------------------------------

def fit_cubic(points):
    """Fit cubic a0 + a1*x + a2*x^2 + a3*x^3 through 4 points."""
    xs = [Fraction(p[0]) for p in points]
    ys = [Fraction(p[1]) for p in points]
    n = 4
    A = [[x ** j for j in range(n)] for x in xs]
    for col in range(n):
        pivot = col
        while A[pivot][col] == 0:
            pivot += 1
        A[col], A[pivot] = A[pivot], A[col]
        ys[col], ys[pivot] = ys[pivot], ys[col]
        for row in range(n):
            if row == col:
                continue
            f = A[row][col] / A[col][col]
            for j in range(n):
                A[row][j] -= f * A[col][j]
            ys[row] -= f * ys[col]
    return [ys[i] / A[i][i] for i in range(n)]

def compute_piecewise_polys():
    """
    For each s in {0,...,23}, compute integer coefficients [D, C, B, A]
    such that 288 * S(N) = A*N^3 + B*N^2 + C*N + D whenever N = s (mod 24).

    Returns dict: s -> [D, C, B, A]
    """
    int_polys = {}
    for s in range(24):
        points = [(s + 24 * (m + 1), S_analytical(s + 24 * (m + 1))) for m in range(4)]
        coeffs = fit_cubic(points)
        ic = []
        for c in coeffs:
            v = c * 288
            assert v.denominator == 1, f"Non-integer coefficient for s={s}"
            ic.append(int(v))
        int_polys[s] = ic
    return int_polys

# ---------------------------------------------------------------------------
# Part 3: Matrix exponentiation for Fibonacci power sums
# ---------------------------------------------------------------------------

def mat_mul_mod(A, B, mod):
    """Multiply two square matrices mod `mod`."""
    n = len(A)
    m = len(B[0])
    k = len(B)
    C = [[0] * m for _ in range(n)]
    for i in range(n):
        for j in range(m):
            s = 0
            for l in range(k):
                s += A[i][l] * B[l][j]
            C[i][j] = s % mod
    return C

def mat_pow_mod(A, p, mod):
    """Matrix exponentiation: A^p mod `mod`."""
    n = len(A)
    result = [[1 if i == j else 0 for j in range(n)] for i in range(n)]
    while p > 0:
        if p & 1:
            result = mat_mul_mod(result, A, mod)
        A = mat_mul_mod(A, A, mod)
        p >>= 1
    return result

def build_transition_matrix(AA, BB, CC, DD, mod):
    """
    Build 14x14 transition matrix for the monomial + accumulator state.

    State vector: [u^3, u^2*v, u*v^2, v^3, u^2, u*v, v^2, u, v, 1,
                   sum_u^3, sum_u^2, sum_u, sum_1]

    Recurrence: u' = AA*u + BB*v, v' = CC*u + DD*v
    Accumulators updated after transition: Su3' = Su3 + u'^3, etc.
    """
    T = [[0] * 14 for _ in range(14)]
    a, b, c, d = AA % mod, BB % mod, CC % mod, DD % mod

    # Row 0: u'^3 = (a*u + b*v)^3
    T[0][0] = a * a % mod * a % mod
    T[0][1] = 3 * a % mod * a % mod * b % mod
    T[0][2] = 3 * a * b % mod * b % mod
    T[0][3] = b * b % mod * b % mod

    # Row 1: u'^2 * v'
    T[1][0] = a * a % mod * c % mod
    T[1][1] = (a * a % mod * d % mod + 2 * a * b % mod * c % mod) % mod
    T[1][2] = (2 * a * b % mod * d % mod + b * b % mod * c % mod) % mod
    T[1][3] = b * b % mod * d % mod

    # Row 2: u' * v'^2
    T[2][0] = a * c % mod * c % mod
    T[2][1] = (2 * a * c % mod * d % mod + b * c % mod * c % mod) % mod
    T[2][2] = (a * d % mod * d % mod + 2 * b * c % mod * d % mod) % mod
    T[2][3] = b * d % mod * d % mod

    # Row 3: v'^3 = (c*u + d*v)^3
    T[3][0] = c * c % mod * c % mod
    T[3][1] = 3 * c % mod * c % mod * d % mod
    T[3][2] = 3 * c * d % mod * d % mod
    T[3][3] = d * d % mod * d % mod

    # Row 4: u'^2
    T[4][4] = a * a % mod
    T[4][5] = 2 * a * b % mod
    T[4][6] = b * b % mod

    # Row 5: u' * v'
    T[5][4] = a * c % mod
    T[5][5] = (a * d + b * c) % mod
    T[5][6] = b * d % mod

    # Row 6: v'^2
    T[6][4] = c * c % mod
    T[6][5] = 2 * c * d % mod
    T[6][6] = d * d % mod

    # Row 7: u'
    T[7][7] = a
    T[7][8] = b

    # Row 8: v'
    T[8][7] = c
    T[8][8] = d

    # Row 9: constant 1
    T[9][9] = 1

    # Row 10: Su3' = Su3 + u'^3
    for j in range(4):
        T[10][j] = T[0][j]
    T[10][10] = 1

    # Row 11: Su2' = Su2 + u'^2
    for j in range(4, 7):
        T[11][j] = T[4][j]
    T[11][11] = 1

    # Row 12: Su' = Su + u'
    T[12][7] = T[7][7]
    T[12][8] = T[7][8]
    T[12][12] = 1

    # Row 13: S1' = S1 + 1
    T[13][9] = 1
    T[13][13] = 1

    for i in range(14):
        for j in range(14):
            T[i][j] %= mod
    return T

# ---------------------------------------------------------------------------
# Part 4: Main solver
# ---------------------------------------------------------------------------

def solve(K=1234567890123):
    """
    Compute sum_{k=2}^{K} S(F_k) mod 10^9.
    """
    MOD_FINAL = 10 ** 9
    WORK_MOD = 288 * MOD_FINAL  # to handle the 1/288 factor exactly

    # Precompute piecewise-cubic polynomials (288 * S(N))
    int_polys = compute_piecewise_polys()

    # Fibonacci mod 24 table (Pisano period of 24 is 24)
    fib_mod24 = [0, 1]
    for i in range(2, 25):
        fib_mod24.append((fib_mod24[-1] + fib_mod24[-2]) % 24)

    # Map each k mod 24 to the polynomial for its Fibonacci residue class
    class_poly = {r: int_polys[fib_mod24[r]] for r in range(24)}

    # Fibonacci constants for 24-step advancement
    # A^24 = [[F_25, F_24], [F_24, F_23]]
    F23, F24, F25 = 28657, 46368, 75025

    # Build 14x14 transition matrix
    T = build_transition_matrix(F23, F24, F24, F25, WORK_MOD)

    # Small Fibonacci values for initial conditions
    fib_small = [0, 1]
    for i in range(2, 30):
        fib_small.append(fib_small[-1] + fib_small[-2])

    total_288 = 0  # accumulate sum of 288*S(F_k), mod WORK_MOD

    for r in range(24):
        poly = class_poly[r]
        D_c = poly[0] % WORK_MOD
        C_c = poly[1] % WORK_MOD
        B_c = poly[2] % WORK_MOD
        A_c = poly[3] % WORK_MOD

        # Index range: k = 24*i + r, k in [2, K]
        i_min = 0 if r >= 2 else 1
        i_max = (K - r) // 24
        if i_max < i_min:
            continue
        n_terms = i_max - i_min + 1

        # Initial Fibonacci values
        k_start = 24 * i_min + r
        u = fib_small[k_start] % WORK_MOD
        v = fib_small[k_start + 1] % WORK_MOD

        # Build initial state (including first term in sums)
        state = [
            u * u % WORK_MOD * u % WORK_MOD,
            u * u % WORK_MOD * v % WORK_MOD,
            u * v % WORK_MOD * v % WORK_MOD,
            v * v % WORK_MOD * v % WORK_MOD,
            u * u % WORK_MOD,
            u * v % WORK_MOD,
            v * v % WORK_MOD,
            u, v, 1,
            u * u % WORK_MOD * u % WORK_MOD,  # Su3 starts with first term
            u * u % WORK_MOD,                   # Su2
            u,                                   # Su
            1,                                   # count
        ]
        state = [x % WORK_MOD for x in state]

        if n_terms > 1:
            T_pow = mat_pow_mod(T, n_terms - 1, WORK_MOD)
            new_state = [0] * 14
            for i in range(14):
                s = 0
                for j in range(14):
                    s += T_pow[i][j] * state[j]
                new_state[i] = s % WORK_MOD
            state = new_state

        Su3, Su2, Su1, S1 = state[10], state[11], state[12], state[13]
        contrib = (A_c * Su3 + B_c * Su2 + C_c * Su1 + D_c * S1) % WORK_MOD
        total_288 = (total_288 + contrib) % WORK_MOD

    return total_288 // 288

# ---------------------------------------------------------------------------
# Part 5: Verification and output
# ---------------------------------------------------------------------------

def verify():
    """Run verification checks."""
    print("=" * 60)
    print("VERIFICATION")
    print("=" * 60)

    # Check M(a,b,c)
    assert M_value(4, 2, 5) == 6, "M(4,2,5) should be 6"
    print("[PASS] M(4,2,5) = 6")

    # Check S(10) and S(10000)
    assert S_analytical(10) == 1972, "S(10) should be 1972"
    print("[PASS] S(10) = 1972")

    assert S_analytical(10000) == 2024258331114, "S(10000) should be 2024258331114"
    print("[PASS] S(10000) = 2024258331114")

    # Brute-force cross-check for S(10)
    assert S_brute(10) == 1972, "S_brute(10) should be 1972"
    print("[PASS] S_brute(10) = 1972 (brute-force confirmed)")

    # Verify polynomial representation for all N in [0, 100]
    int_polys = compute_piecewise_polys()
    for N in range(101):
        s = N % 24
        poly = int_polys[s]
        computed = poly[0] + poly[1] * N + poly[2] * N ** 2 + poly[3] * N ** 3
        actual = 288 * S_analytical(N)
        assert computed == actual, f"Polynomial mismatch at N={N}"
    print("[PASS] Polynomial representation verified for N=0..100")

    # Cross-validate with brute-force Fibonacci sum for K=20
    MOD = 10 ** 9
    fib = [0, 1]
    for i in range(2, 21):
        fib.append(fib[-1] + fib[-2])
    bf_sum = sum(S_analytical(fib[k]) for k in range(2, 21)) % MOD
    matrix_sum = solve(K=20)
    assert matrix_sum == bf_sum, f"K=20 mismatch: matrix={matrix_sum}, brute={bf_sum}"
    print(f"[PASS] K=20 cross-validation: {bf_sum}")

    print()

def generate_visualization():
    """Generate visualization and save as visualization.png."""