euler_435

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

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

\[ \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

Note: I removed the Disqus integration in an effort to cut down on bloat. The following comments were retrieved with the export functionality of Disqus. If you have comments, please reach out to me by Twitter or email.

Michael_Foo_Bar Sep 14, 2013 19:43:21 UTC

Well, checking this with pencil and paper, Sum_x=0,...,3 F_3(x) seems not to be 26, but 92. And the next sum does not match my calculation, too. Can you please check your test data?

Johannes Sep 14, 2013 20:09:29 UTC

You are absolutely right. I got tricked by the upper limit of Python's range function. The value should now be correct.

essay writing Mar 08, 2016 05:43:48 UTC

It's actually good that you have shared about this kind of information in order for those people to know some effective strategies on how are they going to promote some good solution in case if they encounter that certain kind of problem.