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.

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.

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.

Reference

• Sum of geometric progression
• Matrix form of Fibonacci sequence
• Exponentiation by squaring
• Modular arithmetic

Test Data

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

<td>
  36651874875
</td>

<td>
  92
</td>

<td>
  14
</td>

<td>
  960038235750
</td>

Solution

Simplify the polynomial

The generating function

$$ \begin{equation} F_n(x) = f_1 x + f_2 x^2 + f_3 x^3 + f_4 x^4 + \cdots + f_n x^n \end{equation} $$

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:

$$ \begin{align} F_n(x) &= f_1 x + f_2 x^2 + f_3 x^3 + f_4 x^4 + \cdots + f_n x^n \label{eq:a} \ x\cdot F_n(x) &= \phantom{f_1 x +}; f_1 x^2 + f_2 x^3 + f_3 x^4 + f_4 x^5 + \cdots + f_n x^{n+1} \label{eq:b} \ x^2\cdot F_n(x) &=\phantom{f_1 x + f_2 x^2 + }; f_1 x^3 + f_2 x^4 + f_3 x^5 + f_4 x^6 + \cdots + f_n x^{n+2} \label{eq:c} \end{align} $$

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

$$ \begin{equation} \begin{split} F_n(x)(1-x-x^2) =& f_1 x + (f_2-f_1) x^2 + (f_3 - f_2 - f_1) x^3 \ &+ (f_4 - f_3 - f_2) x^4 + \cdots + (f_n - f_{n-1} - f_{n-2}) x^n \ &- f_n x^{n+1} - f_{n-1} x^{n+1} - f_n x^{n+2} \end{split} \end{equation} $$

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:

$$ \begin{equation} F_n(x) = \frac{ f_n x^{n+2} + f_{n+1} x^{n+1} - x}{x^2+x-1} \end{equation} $$

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:

$$ \begin{equation} f_n = (1, 0) \begin{pmatrix}1 & 1 \1 & 0\end{pmatrix}^n \begin{pmatrix}1\0\end{pmatrix} \end{equation} $$

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:

$$ \begin{equation} \frac{a}{b} \text{ mod } m \neq \frac{a \text{ mod } m}{b \text{ mod } m} \end{equation} $$

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:

$$ \begin{equation} \frac{a}{b} \text{ mod } m = \frac{a \text{ mod } m\cdot b}{b}, \text{ if } a \text{ mod } b \equiv 0 \end{equation} $$

If we apply that to Equation (2):

$$ \begin{equation} \begin{split} F_n(x) \text{ mod } m = \frac{(f_n x^{n+2} + f_{n+1} x^{n+1} - x) \text{ mod } \big(m(x^2+x-1)\big)}{ (x^2+x-1)} \end{split} \end{equation} $$

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))

Archived Comments