Add Me!Close Menu Navigation

Calculating Modified Fibonacci Numbers, Quickly and Exactly

This is inspired by an article of similar name called Calculating Fibonacci Numbers, Quickly and Exactly. Consider the a modified fibonacci sequence given by the recurrence c_0 = 2, c_1 = 1, c_k = c_{k-1} + c_{k-2}. We shall look at various algebraic and analytic properties of this sequence, find two algorithms that computes c_k in O(\log(k)) time, and finally look at a motivating example where we need to compute this sequence quickly and exactly.

1. Properties of c_k

The sequence c_k is essentially the same as the Fibonacci sequence, with the exception that the initial conditions are altered slightly (so that c_0 \ne f_0, where f_k is the fibonacci sequence). This similarity lets us use many of the same algebraic machinery that can be used to tackle the fibonacci numbers.

Closed Form

I will defer the derivation until the appendix, but it can be quickly shown using a little bit of linear algebra that letting \phi,\psi = \frac{1\pm \sqrt{5}}{2} be the golden ratio (and it’s “conjugate” wrt to the root), then

    \[c_k = \phi^k + \psi^k\]

which looks to have the same form as the fibonacci numbers.

Analytic Property

The numerical value of \phi \approx 1.618\dots and \psi = -\phi^{-1}, therefore as k grows large, the contribution of \psi^k becomes insignificant, serving only to ensure that c_k is a whole number, therefore c_k = O(\phi^k), in fact, from simple observation, it seems that for k > 1, c_k = \left\lfloor\phi^k + \frac{1}{2}\right\rfloor (or \phi^k rounded to the nearest integer).

Computation by Reduction to Fibonacci Number

Unfortunately, we don’t have quite a good way to compute with irrational objects: if we use floating point arithmetic, we quickly lose information about the trailing terms, if we use some algebraic package, then simplifying expressions with irrational numbers may become computationally infeasible.

However, we can exploit the connection between c_k and f_k to help us out. Suppose we have some software that can efficiently compute f_k for any arbitrary k, can we possibly use that to our advantage in computing c_k?

It can be shown (similar to the steps used in the appendix) that \sqrt{5}f_k = \phi^k - \psi^k. Notice that

    \[\underbrace{\pa{\phi^k - \psi^k}}_{\sqrt{5}f_k}\underbrace{\pa{\phi^k + \psi^k}}_{c_k} = \phi^{2k} - \psi^{2k} = \sqrt{5}f_{2k}\]

therefore

    \[c_k = \frac{f_{2k}}{f_k} = 2f_{k+1} - f_k\]

Direct Computation

Alternatively, we can also compute this sequence directly in O(\log(k)) time. Recall that \phi\psi = -1, then

    \[c_k^2 = \phi^{2k} + 2\phi^k\psi^k + \psi^{2k} = c_{2k} + 2(-1)^k\]

so c_{2k} = c_k^2 - 2(-1)^k. Unfortunately, there’s no elegant solution for all odd terms, but recall that within the consecutive triple (k+1,k,k-1), at least one of these must be a multiple of three. Now, consider

    \begin{align*} c_k^3 &= c_k^2(\phi^k + \psi^k) \\  &= \pa{\phi^{2k} + \psi^{2k}}\pa{\phi^k + \psi^k} + 2(-1)^k c_k \\ &= \phi^{3k} + \psi^{3k} + \phi^{2k}\psi^k + \psi^{2k}\phi^k + 2(-1)^k c_k \\ &= c_{3k} + (-1)^k c_k + 2(-1)^k c_k \\ &= c_{3k} + 3(-1)^k c_k \\ c_{3k} &= c_k(c_k^2 - 3(-1)^k) \end{align*}

Therefore, we can compute c_k using

    \[c_0 = 2, c_1 = 1, c_k = \begin{cases} c_{k/2}^2 - 2(-1)^{k/2} & \mbox{if $k$ is even} \\ c_{k/3}(c_{k/3}^2 - 3(-1)^{k/3}) & \mbox{if $k$ is a multiple of three} \\ \underbrace{c_{k-1}}_{\g{even}} + \underbrace{c_{k-2}}_{\g{multiple of 3}} & \mbox{otherwise} \end{cases}\]

Memoizing gives a O(\log(k)) solution

C = {0:2,1:1}
d = lambda k: -1 if k % 2 else 1
def c(k):
    if k in C:
        return C[k]
    if k % 2 is 0:
        r = c(k/2)
        C[k] = r*r - 2*d(k/2)
    elif k % 3 is 0:
        r = c(k/3)
        C[k] = r*(r*r - 3*d(k/3))
    else:
        C[k] = c(k-1) + c(k-2)
    return C[k]

2. Example Problem: Fibonacci Graph

Suppose we have the Fibonacci sequence

    \[0, 1, 1, 2, 3, 5, 8, 13, 21, \dots\]

defined by

    \[f_0 = 0 , f_1 = 1 , f_{k+1} = f_k + f_{k-1}.\]

Let’s construct a sequence of points P_k = \pa{x_k,y_k} such that

    \[P_0 = (0,0), P_k = (f_{k-1},f_k)\]

so that it looks like the sequence

    \[(0,0), (0,1), (1,1), (1,2), (2,3), (3,5), (5,8), \dots\]

Let’s connect all of the points in (P_k)_k together (so that it forms a jagged line) and look at the area under the curve; we will call this quantity B_n.

Give an algorithm that computes B_n \mod 10^{10} in O(\log(n)) time and compute the area for the case where n = 1234567891011121314.
splittingplot

Solution Sketch

If you look at the figure above, you’ll notice that the total region is made up of these trapezoids. Let’s call A_k the area of the trapezoid defined by P_k and P_{k+1}, then as illustrated below,

Rendered by QuickLaTeX.com

we know that

    \[A_k = (x_{k+1} - x_k) y_k + \frac{(y_{k+1} - y_k)(x_{k+1} - x_k)}{2},\]

which is just the area of the rectangle on the bottom plus the triangle on the top. From the relation x_{k+1} = x_k + x_{k-1}, we know that x_{k+1} - x_k = x_{k-1}, and similarly for y, so that we have

    \[A_k = x_{k-1} \pa{y_k + \frac{y_{k-1}}{2}} = f_{k-2}\pa{f_k + \frac{f_{k-1}}{2}}\]

Let’s make a connection to the exponential form of the Fibonacci sequence:

    \[f_k = \frac{\phi^k - \psi^k}{\sqrt{5}}\]

where \phi,\psi = \frac{1 \pm \sqrt{5}}{2} are the golden ratio and its conjugate respectively. If we actually substitute this equation into the definition of A_k and expand everything out, we will find that

    \[A_k = \frac{1}{10} \phi^{2k-3} + \frac{1}{5}\phi^{2k-2} - \frac{1}{10} \phi^{k-1}\psi^{k-2} - \frac{1}{5} \phi^k\psi^{k-2} - \frac{1}{10} \phi^{k-2}\psi^{k-1} - \frac{1}{5}\phi^{k-2}\psi^k + \frac{1}{10} \psi^{2k-3} + \frac{1}{5}\psi^{2k-2}\]

Since \psi is \phi‘s conjugate, we have \phi\psi = -1, so each of the \phi^{k-2+a}\psi^{k-2+b} term can be rewritten as (-1)^{k-2} \phi^a\psi^b = (-1)^k \phi^a\psi^b. Grouping related terms together, we get

    \[A_k = \phi^{2k} \underbrace{\pa{\frac{1}{10} \phi^{-3} + \frac{1}{5} \phi^{-2}}}_{10^{-1}} + \psi^{2k} \underbrace{\pa{\frac{1}{10} \psi^{-3} + \frac{1}{5} \psi^{-2}}}_{10^{-1}} - (-1)^k \pa{\frac{\overbrace{\phi + \psi}^{1}}{10} + \frac{\overbrace{\phi^2 + \psi^2}^{3}}{5}},\]

simplifying, we find that

    \[A_k = \frac{\phi^{2k} + \psi^{2k} - 7 (-1)^k }{10}, A_0 = 0\]

This is awesome! But we need to find the sum of B_k = A_0 + A_1 + \cdots + A_{k-1}, so let’s do just that

    \begin{align*} B_{k+1} &= \sum_{j=0}^{k} A_j \\ &= \sum_{j=1}^{k-1} \frac{\phi^{2j} + \psi^{2j} - 7 (-1)^j }{10} \\ \intertext{using the geometric series} &= \frac{1}{10} \pa{ \frac{\phi^{2k+2} - \phi^2}{\phi^2 - 1} + \frac{\psi^2 - \psi^{2k+2}}{1 - \psi^2} + 7 \frac{1 - (-1)^k}{2} } \\ \intertext{let's define $\gamma_k = \begin{cases}0 & k \g{ even} \\ 1 & \g{otherwise}\end{cases}$}  &= \frac{1}{10} \pa{ \frac{\phi^{2k+2} - \phi^2}{\phi^2 - 1} + \frac{\psi^2 - \psi^{2k+2}}{1 - \psi^2} + 7\gamma_k } \\ \intertext{furthermore, notice that $\phi = \phi^2 - 1$ and likewise with $\psi$, so} &= \frac{1}{10} \pa{ \frac{\phi^{2k+2} - \phi^2}{\phi} + \frac{\psi^2 - \psi^{2k+2}}{- \psi} + 7 \gamma_k} \\ &= \frac{ \phi^{2k+1} + \psi^{2k+1} - \overbrace{\pa{\phi + \psi}}^{-1} +  7 \gamma_k}{10} \\ &= \frac{\phi^{2k+1} + \psi^{2k+1} - 1 + 7 \gamma_k}{10} \\ &= \frac{c_{2k+1} - 1 + 7\gamma_k}{10} \end{align*}

Sweet!

Since I was a little vague on specifying the range of trapezoids to consider for B_n, I will accept both 7997480508.5 (the literalist interpretation) and 7203060060 (the more casual interpretation).




3. Appendix

Deriving Closed Form of c_k

Consider the sequence (c_k), we can represent this sequence as an infinite dimensional vector called c (so that indexing c is the same as finding the k^{th} element of the sequence). Now, consider the vector-space in which c inhabits, we can define a transformation R (that I will call the right-shift operator) that basically takes in x and outputs x shifted one index over, so that (Rc)_k = c_{k+1}. e.g.

    \[R\mat{x_0\\x_1\\x_2\\\vdots} = \mat{x_1\\x_2\\x_3\\\vdots}\]

It’s clear that this is a linear transform, because we can inductively define R as the limit of the matrix-valued sequence of matrices whose first upper diagonal is constant 1:

    \[R_k = \mat{0 & I_{k-1} \\ 0 & 0} = \mat{ 0 & 1 & 0 & 0 & 0 \\0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & \ddots & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0}\]

Consider the subspace spanned by all of the shifted c sequences:

    \[\mathcal{C} = \g{span}\pa{c, Rc, R^2c, \cdots} = \g{span}(c,Rc)\]

the second equality comes from the fact that because c_{k+2} = c_{k+1} + c => R^2c = Rc + c. Now, because we don’t actually know what c is (only that it exists), we can alternatively look at it as the kernel of the system

    \[(R - \phi)(R-\psi)c = 0\]

(because we’ve already established above that R is a linear operator). Furthermore, R and all scalar matrices commute. Appealing to a little bit of algebra, we know that, letting \rho(x) = x^2 - x - 1 be the characteristic and \phi,\psi be scalar roots of the quadratic x^2-x-1 such that \rho(\phi) = \rho(\psi) = 0, then

    \[\rho(R) = (R - \phi)(R-\psi) = (R-\psi)(R-\phi)\]

Therefore, let (R-\phi)a = 0 and (R-\psi)b = 0, then the solution c must be some linear combination of the solutions a and b:

    \[c = \alpha a + \beta b\]

for some scalar \alpha, \beta, because then

    \begin{align*} \rho(R)c &= (R^2 - R - 1) c \\ &= \alpha(R - \psi)\underbrace{(R-\phi)a}_0 + \beta (R-\phi)\underbrace{(R-\psi)b}_0 \\ &= 0 \end{align*}

as expected.

Now, we know that Ra = \phi a, but this is the same as saying that letting a_0 = 1, then a_{k+1} = \phi a_k, so a_k = \phi^k. Using a similar argument, we can also argue that b_k = \psi^k, so

    \[c_k = \alpha \phi^k + \beta \psi^k\]

We can then use our initial values that c_0 = 2, c_1 = 1 to solve for \alpha, \beta

    \[\mat{1 & 1\\ \phi & \psi} \mat{\alpha\\\beta} = \mat{2\\1} => \mat{\alpha\\\beta} = \frac{1}{\psi-\phi}\mat{\psi & -1 \\ -\phi & 1}\mat{2\\1} = \mat{1\\1}\]

so

    \[c_k = \phi^k + \psi^k\]

Posted By Lee

Woosh

  • AM

    Nice post; a few additional notes:

    1. These “modified” Fibonacci numbers are often called the Lucas numbers

    2. A lot of the above analysis also falls out of the matrix form of the Fibonacci recurrence, as seen here:
    http://en.wikipedia.org/wiki/Fibonacci_number#Matrix_form

    The relevant exponential bases fall out of the eigenvalues of the 2×2 matrix there, and the O(log k)-arithmetic-operation algorithm can be interpreted as exponentiation-by-repeated-squaring.

    • http://phailed.me/ Phailure

      Thanks AM, I meant to add that in as an addendum but forgot. Have a great day now! :)