Problem 402: Integer-valued Polynomials
View on Project EulerProject 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
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
- Problem page: https://projecteuler.net/problem=402
- Integer-valued polynomial: Wikipedia — Integer-valued polynomial
- Binomial coefficient basis: Wikipedia — Binomial coefficient
- Pisano periods: Wikipedia — Pisano period
- Matrix exponentiation: Wikipedia — Matrix exponentiation
- Quasi-polynomials: Wikipedia — Quasi-polynomial
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());
}
}