Problem 435: Polynomials of Fibonacci Numbers
View on Project EulerProject Euler Problem 435 Solution
EulerSolve provides an optimized solution for Project Euler Problem 435, Polynomials of Fibonacci Numbers, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Let the Fibonacci numbers be \(f_0=0\), \(f_1=1\), and \(f_n=f_{n-1}+f_{n-2}\). For each integer \(x\), define the truncated polynomial $$F_n(x)=\sum_{k=0}^{n} f_k x^k.$$ The task is to evaluate $$\sum_{x=0}^{100} F_n(x) \pmod{15!}$$ for the enormous index \(n=10^{15}\). A direct loop over all Fibonacci terms is impossible, so the key is to turn the polynomial sum into a fixed-size linear recurrence that can be exponentiated in logarithmic time. Mathematical Approach Step 1: Isolate the weighted Fibonacci term For a fixed value of \(x\), define $$a_k=f_k x^k.$$ Using the Fibonacci recurrence, $$a_{k+1}=f_{k+1}x^{k+1}=(f_k+f_{k-1})x^{k+1}=x a_k + x^2 a_{k-1}.$$ So once \(x\) is fixed, the sequence \(a_k\) satisfies a second-order linear recurrence with constant coefficients \(x\) and \(x^2\). Step 2: Add the partial sum to the state vector The polynomial value itself is the prefix sum of the \(a_k\): $$F_k(x)=\sum_{i=0}^{k} a_i.$$ Therefore $$F_{k+1}(x)=F_k(x)+a_{k+1}.$$ Now combine the recurrence for \(a_k\) with the update for the running sum....
Detailed mathematical approach
Problem Summary
Let the Fibonacci numbers be \(f_0=0\), \(f_1=1\), and \(f_n=f_{n-1}+f_{n-2}\). For each integer \(x\), define the truncated polynomial
$$F_n(x)=\sum_{k=0}^{n} f_k x^k.$$
The task is to evaluate
$$\sum_{x=0}^{100} F_n(x) \pmod{15!}$$
for the enormous index \(n=10^{15}\). A direct loop over all Fibonacci terms is impossible, so the key is to turn the polynomial sum into a fixed-size linear recurrence that can be exponentiated in logarithmic time.
Mathematical Approach
Step 1: Isolate the weighted Fibonacci term
For a fixed value of \(x\), define
$$a_k=f_k x^k.$$
Using the Fibonacci recurrence,
$$a_{k+1}=f_{k+1}x^{k+1}=(f_k+f_{k-1})x^{k+1}=x a_k + x^2 a_{k-1}.$$
So once \(x\) is fixed, the sequence \(a_k\) satisfies a second-order linear recurrence with constant coefficients \(x\) and \(x^2\).
Step 2: Add the partial sum to the state vector
The polynomial value itself is the prefix sum of the \(a_k\):
$$F_k(x)=\sum_{i=0}^{k} a_i.$$
Therefore
$$F_{k+1}(x)=F_k(x)+a_{k+1}.$$
Now combine the recurrence for \(a_k\) with the update for the running sum. With the state vector
$$v_k=\begin{bmatrix}a_k\\ a_{k-1}\\ F_k(x)\end{bmatrix},$$
we obtain the linear transition
$$v_{k+1}=T_x v_k,\qquad T_x=\begin{bmatrix} x & x^2 & 0\\ 1 & 0 & 0\\ x & x^2 & 1 \end{bmatrix}.$$
The first row reproduces \(a_{k+1}=x a_k+x^2 a_{k-1}\), the second row shifts \(a_k\) into the next position, and the third row accumulates the new term into the polynomial sum.
Step 3: Base state and matrix power
For \(k=1\),
$$a_1=f_1 x=x,\qquad a_0=f_0=0,\qquad F_1(x)=f_0+f_1 x=x.$$
Hence the base vector is
$$v_1=\begin{bmatrix}x\\ 0\\ x\end{bmatrix}.$$
For every \(n\ge 1\),
$$v_n=T_x^{\,n-1}v_1,$$
and the desired value \(F_n(x)\) is the third component of \(v_n\). The special case \(n=0\) is immediate because \(F_0(x)=f_0=0\).
Step 4: Why repeated squaring is the right tool
The matrix \(T_x\) has fixed size \(3\times 3\), so its \((n-1)\)-st power can be computed by binary exponentiation in \(O(\log n)\) matrix multiplications. Every operation is performed modulo
$$M=15!=1307674368000.$$
This keeps the numbers bounded while matching the modulus required by the problem.
Alternative closed form and why the implementation avoids division
The truncated generating function also satisfies
$$\left(1-x-x^2\right)F_n(x)=x-f_{n+1}x^{n+1}-f_n x^{n+2}.$$
Over the integers this identity is correct and useful, but modulo \(15!\) the factor \(1-x-x^2\) is not always invertible. For example, at \(x=2\) it equals \(-5\), and \(5\) shares a factor with \(15!\). The matrix formulation is therefore preferable: it is division-free and works uniformly for every \(x\in\{0,\dots,100\}\).
Worked checks
Several small cases confirm the recurrence.
First, \(F_n(0)=0\) for all \(n\), because every term contains either \(f_0=0\) or a positive power of \(0\).
Second, for \(n=7\) and \(x=11\),
$$F_7(11)=11+11^2+2\cdot 11^3+3\cdot 11^4+5\cdot 11^5+8\cdot 11^6+13\cdot 11^7=268357683.$$
This matches the value produced by the matrix method.
A larger checkpoint is
$$\sum_{x=0}^{10} F_{20}(x)\equiv 1044074802100 \pmod{15!},$$
which is exactly the validation target used by the implementations.
How the Code Works
The C++, Python, and Java implementations all follow the same structure. For each \(x\in\{0,\dots,100\}\), they build the transition matrix \(T_x\), raise it to the power \(n-1\) by repeated squaring, multiply by the base vector, and read the third component as \(F_n(x)\). These 101 values are then accumulated modulo \(15!\).
Because \(15!\approx 1.3\times 10^{12}\), an unreduced product of two residues can exceed ordinary 64-bit multiplication in some languages. The implementations therefore perform modular multiplication carefully before each reduction, but the mathematical algorithm itself remains the same in all three languages.
Complexity Analysis
A \(3\times 3\) matrix multiplication has constant cost, and binary exponentiation uses \(O(\log n)\) such multiplications. Therefore each single value \(F_n(x)\) is computed in \(O(\log n)\) time and \(O(1)\) extra space. Since the outer sum runs over exactly 101 values of \(x\), the full computation is still \(O(101\log n)=O(\log n)\) with a very small constant factor, and the memory usage stays \(O(1)\).
References
- Problem page: https://projecteuler.net/problem=435
- Fibonacci numbers: Wikipedia — Fibonacci number
- Matrix exponentiation / exponentiation by squaring: Wikipedia — Exponentiation by squaring
- Generating functions for Fibonacci numbers: Wikipedia — Fibonacci number, generating function section
Problem 435 source code
C++
#include <array>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>
namespace {
using u64 = std::uint64_t;
using u128 = __uint128_t;
constexpr u64 MOD = 1307674368000ULL; // 15!
struct Options {
u64 n = 1000000000000000ULL;
int x_max = 100;
bool run_checkpoints = true;
};
bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& 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::stoi(tail);
} catch (...) {
return false;
}
return true;
}
bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& value) {
if (arg.rfind(prefix, 0U) != 0U) {
return false;
}
const std::string tail = arg.substr(prefix.size());
if (tail.empty()) {
return false;
}
try {
value = static_cast<u64>(std::stoull(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_u64_after_prefix(arg, "--n=", options.n) ||
parse_int_after_prefix(arg, "--x-max=", options.x_max)) {
continue;
}
std::cerr << "Unknown argument: " << arg << '\n';
return false;
}
return options.x_max >= 0;
}
using Mat3 = std::array<std::array<u64, 3>, 3>;
using Vec3 = std::array<u64, 3>;
Mat3 mat_mul(const Mat3& a, const Mat3& b) {
Mat3 c{};
for (int i = 0; i < 3; ++i) {
for (int k = 0; k < 3; ++k) {
if (a[i][k] == 0ULL) {
continue;
}
const u64 aik = a[i][k];
for (int j = 0; j < 3; ++j) {
if (b[k][j] == 0ULL) {
continue;
}
c[i][j] = static_cast<u64>((c[i][j] + static_cast<u128>(aik) * b[k][j]) % MOD);
}
}
}
return c;
}
Vec3 mat_vec_mul(const Mat3& a, const Vec3& v) {
Vec3 out{};
for (int i = 0; i < 3; ++i) {
u128 sum = 0;
for (int j = 0; j < 3; ++j) {
if (a[i][j] == 0ULL || v[j] == 0ULL) {
continue;
}
sum += static_cast<u128>(a[i][j]) * v[j];
}
out[i] = static_cast<u64>(sum % MOD);
}
return out;
}
Mat3 mat_pow(Mat3 base, u64 exp) {
Mat3 result{};
for (int i = 0; i < 3; ++i) {
result[i][i] = 1ULL;
}
u64 e = exp;
while (e > 0ULL) {
if (e & 1ULL) {
result = mat_mul(base, result);
}
e >>= 1ULL;
if (e > 0ULL) {
base = mat_mul(base, base);
}
}
return result;
}
u64 F_value(u64 n, u64 x_raw) {
if (n == 0ULL) {
return 0ULL;
}
const u64 x = x_raw % MOD;
const u64 x2 = static_cast<u64>((static_cast<u128>(x) * x) % MOD);
// a_n = f_n * x^n obeys: a_{n+1} = x*a_n + x^2*a_{n-1}.
// State v_n = [a_n, a_{n-1}, F_n]^T with transition:
// [a_{n+1}] [x x^2 0] [a_n ]
// [a_n ] = [1 0 0] [a_{n-1}]
// [F_{n+1}] [x x^2 1] [F_n ].
const Mat3 t{{
{{x, x2, 0ULL}},
{{1ULL, 0ULL, 0ULL}},
{{x, x2, 1ULL}},
}};
const Vec3 v1{{x, 0ULL, x}}; // n=1
const Mat3 p = mat_pow(t, n - 1ULL);
const Vec3 vn = mat_vec_mul(p, v1);
return vn[2];
}
u64 brute_small(int n, int x) {
std::vector<u64> fib(static_cast<std::size_t>(n + 1), 0ULL);
if (n >= 1) {
fib[1] = 1ULL;
}
for (int i = 2; i <= n; ++i) {
fib[static_cast<std::size_t>(i)] = fib[static_cast<std::size_t>(i - 1)] + fib[static_cast<std::size_t>(i - 2)];
}
u64 p = 1ULL;
u64 sum = 0ULL;
for (int i = 0; i <= n; ++i) {
sum = (sum + static_cast<u128>(fib[static_cast<std::size_t>(i)] % MOD) * p) % MOD;
p = static_cast<u64>((static_cast<u128>(p) * static_cast<u64>(x)) % MOD);
}
return sum;
}
u64 solve_sum(u64 n, int x_max) {
u64 ans = 0ULL;
for (int x = 0; x <= x_max; ++x) {
ans += F_value(n, static_cast<u64>(x));
ans %= MOD;
}
return ans;
}
bool run_checkpoints() {
if (F_value(7ULL, 11ULL) != 268357683ULL) {
std::cerr << "Checkpoint failed: F_7(11)\n";
return false;
}
if (F_value(20ULL, 7ULL) != brute_small(20, 7)) {
std::cerr << "Checkpoint failed: matrix vs brute for n=20,x=7\n";
return false;
}
if (solve_sum(20ULL, 10) != 1044074802100ULL) {
std::cerr << "Checkpoint failed: sum_{x=0..10} F_20(x)\n";
return false;
}
return true;
}
} // namespace
int main(int argc, char** argv) {
Options options;
if (!parse_arguments(argc, argv, options)) {
return 1;
}
if (options.run_checkpoints && !run_checkpoints()) {
return 2;
}
std::cout << solve_sum(options.n, options.x_max) << '\n';
return 0;
}
Python
def solve():
MOD = 1307674368000
n = 1000000000000000
x_max = 100
def mat_mul(a, b):
c = [[0] * 3 for _ in range(3)]
for i in range(3):
for k in range(3):
if a[i][k] == 0:
continue
aik = a[i][k]
for j in range(3):
if b[k][j] == 0:
continue
c[i][j] = (c[i][j] + aik * b[k][j]) % MOD
return c
def mat_vec_mul(a, v):
out = [0] * 3
for i in range(3):
val = 0
for j in range(3):
if a[i][j] != 0 and v[j] != 0:
val = (val + a[i][j] * v[j]) % MOD
out[i] = val
return out
def mat_pow(base, exp):
result = [[0] * 3 for _ in range(3)]
for i in range(3):
result[i][i] = 1
e = exp
while e > 0:
if e & 1:
result = mat_mul(base, result)
e >>= 1
if e > 0:
base = mat_mul(base, base)
return result
def f_value(x_raw):
if n == 0:
return 0
x = x_raw % MOD
x2 = (x * x) % MOD
t = [
[x, x2, 0],
[1, 0, 0],
[x, x2, 1]
]
v1 = [x, 0, x]
p = mat_pow(t, n - 1)
vn = mat_vec_mul(p, v1)
return vn[2]
ans = 0
for x in range(x_max + 1):
ans = (ans + f_value(x)) % MOD
return str(ans)
if __name__ == '__main__':
print(solve())
Java
public class Euler435 {
private static final long MOD = 1307674368000L;
private static final long N_VAL = 1000000000000000L;
private static final int X_MAX = 100;
private static long[][] matMul(long[][] a, long[][] b) {
long[][] c = new long[3][3];
for (int i = 0; i < 3; i++) {
for (int k = 0; k < 3; k++) {
if (a[i][k] == 0)
continue;
// Since multiplication can exceed long capacity before modulo,
// we'll carefully handle it using big decimal or a safe multiplier,
// but MOD is ~10^12, so a[i][k] * b[k][j] can be ~10^24 which overflows long
// (max ~9*10^18).
// Use a safe multiply mod function or split math.
for (int j = 0; j < 3; j++) {
if (b[k][j] == 0)
continue;
long add = mulMod(a[i][k], b[k][j]);
c[i][j] = (c[i][j] + add) % MOD;
}
}
}
return c;
}
private static long mulMod(long a, long b) {
long q = (long) ((double) a * b / MOD);
long r = a * b - q * MOD;
while (r >= MOD)
r -= MOD;
while (r < 0)
r += MOD;
return r;
}
private static long[] matVecMul(long[][] a, long[] v) {
long[] out = new long[3];
for (int i = 0; i < 3; i++) {
long sum = 0;
for (int j = 0; j < 3; j++) {
if (a[i][j] == 0 || v[j] == 0)
continue;
sum = (sum + mulMod(a[i][j], v[j])) % MOD;
}
out[i] = sum;
}
return out;
}
private static long[][] matPow(long[][] base, long exp) {
long[][] result = new long[3][3];
for (int i = 0; i < 3; i++) {
result[i][i] = 1;
}
long e = exp;
while (e > 0) {
if ((e & 1) != 0) {
result = matMul(base, result);
}
e >>= 1;
if (e > 0) {
base = matMul(base, base);
}
}
return result;
}
private static long fValue(long xRaw) {
if (N_VAL == 0)
return 0;
long x = xRaw % MOD;
long x2 = mulMod(x, x);
long[][] t = {
{ x, x2, 0 },
{ 1, 0, 0 },
{ x, x2, 1 }
};
long[] v1 = { x, 0, x };
long[][] p = matPow(t, N_VAL - 1);
long[] vn = matVecMul(p, v1);
return vn[2];
}
public static String solve() {
long ans = 0;
for (int x = 0; x <= X_MAX; x++) {
ans = (ans + fValue(x)) % MOD;
}
return String.valueOf(ans);
}
public static void main(String[] args) {
System.out.println(solve());
}
}