Project Euler Problem 435 – Polynomials of Fibonacci Numbers

Spoiler Alert! This blog entry gives away the solution to problem 435 of Project Euler. Please don’t read any further if you have yet to attempt to solve the problem on your own. The information is intended for those who failed to solve the problem and are looking for hints or test data to help them track down bugs.</br>
First I provide references of the major theoretical background that you probably need to know to solve the problem and test data to validate your own algorithm. The last section presents the solution from start to finish.</br>
The post reflects my approach to the problem. Even though the final outcome was accepted by Project Euler doesn’t mean the information is correct or elegant. Algorithms won’t necessarily be the most efficient ones, but they are guaranteed to run within the time limit of one minute.

Test Data

Edit: This post originally listed wrong values. Thanks, Michael_Foo_Bar, for pointing out the errors.

 $f_{10^15} \text{ mod } 15!$ 36651874875 $\sum_{x=0}^{3} F_3(x)$ 92 $\sum_{x=0}^{10} F_{10}(x) \text{ mod } 27$ 14 $F_{10^15}(2) \text{ mod } 15!$ 960038235750

Solution

Simplify the polynomial

The generating function

can be simplified the same way the Sum of geometric progression is derived. First, multiply (1) by $x$ , and by $x^2$ respectively to obtain the following set of equations:

Now subtract the last two equations from the first to get:

Most of the terms like $(f_3 - f_2 - f_1)$ drop out because $f_n = f_{n-1} + f_{n-2}$ , and $F_n(x)$ is:

Calculate the Fibonacci terms

Equation (2) still requires Fibonacci terms for very large n. The closed-form won’t work. Instead, the matrix form turns out to be useful:

Exponentiation by squaring helps to calculate the nth power of the matrix – and the powers of x – fast.

Apply modular arithmetic

Since Fibonacci numbers grow large very fast we will only be able to calculate the modulus. Unfortunately we can’t just determine all values in (2) modulo 15! and end up with the desired result, because the following does not hold:

One way to do division is to use the modular multiplicative inverse, but those aren’t always defined in our instance. The Chinese Remainder Theorem would help fix that, but there is a simpler solution:

If we apply that to Equation (2):

The numerator only uses multiplication and addition which allows to calculate the Fibonacci terms as well as the powers of x modulo $m\cdot(x^2+x-1)$ .

Putting it all together

def mat_mult(A, B, m):
C = [[0,0],[0,0]]
for i in range(2):
for j in range(2):
for k in range(2):
C[i][j] += A[i][k]*B[k][j]
C[i][j] %= m
return C

def mat_pow(A, p, m):
if p == 1:
return A
if p % 2:
return mat_mult(A, mat_pow(A, p-1, m), m)
X = mat_pow(A, p//2, m)
return mat_mult(X, X, m)

def fib(n, m):
T = [[1,1], [1,0]]
if n == 1:
return 1
T = mat_pow(T, n-1, m)
return T[0][0]

def mpow(x, p, m):
if p == 1:
return x
if p % 2:
return x*mpow(x, p-1, m) % m
r = mpow(x, p//2, m)
return r*r % m

def F(n, x, m):
b = (x*x + x - 1)
f_nn = fib(n+1, m*b)
f_n = fib(n, m*b)
a = (f_nn*mpow(x, n+1, m*b) + f_n*mpow(x, n+2, m*b) - x) % (m*b)
return a/b

def calc(n,m,x_range):
R = 0
for x in x_range:
R += F(n, x, m)
R %= m
return R

if __name__ == "__main__":
n = pow(10,15)
m = 1307674368000
x_range = range(0, 101)
print(calc(n, m, x_range))