[mynegation]

Very good demonstration of subsequent improvements of a naive algorithm. To me that was somewhat depreciated by the fact that you can actually calculate n-the Fibonacci number using Binet's closed form formula (http://en.wikipedia.org/wiki/Fibonacci_number#Closed-form_ex...). You will need arbitrary precision arithmetic starting with certain 'n' though, as IEEE 754 will not give you correct result.

[stralep]

Actually, you need just Z ring enriched (if this is right word) with sqrt(5). So instead of one float, you use pair of integers of arbitrary precision. So

(a,b)+(c,d) = (a+c,b+d)

(a,b)/sqrt(5) = (b,a/5)

(a,b)(c,d) = (ac+5bd,ad+bc)

2phi = (1,1)

F(n) = (2phi^n - (2-2phi)^n)/(2^n*sqrt(5))

[edit] As ot said, the right word is "extended"

[ot]

> Z ring enriched (if this is right word) with sqrt(5)

I think "extended" is the right word (see http://en.wikipedia.org/wiki/Ring_extension)

It is commonly denoted as Z[sqrt(5)]

[vog]

This is a great method! It avoids arbitrary floating point arithmetics[1] and thus makes it more comparable with the method presented in the article.

As far as I can see, this algorithm is very similar to the matrix operations of the article, at least in terms of complexity.

[1] Only the very last step might involve fp arithmetics, because you'll have to convert the result pair (x,y) back to x+y*sqrt(5). However, we already know that the result has to be an integer, so the its second part will always be 0, no fp arithmetics required.

[btilly]

And when you calculate this out, you internally wind up doing the same calculation as the matrix method.

[vog]

How do you know that? Have you actually done the math?

Although I agree that this whole stuff looks very similar, I'm not sure whether it really leads to the exact same calculations.

[onedognight]

If you compute the eigenvalue decomposition of

  A={{0,1},{1,1}}

i.e.

  A = V^-1 D V

where D is a diagonal matrix with the eigenvalues phi and 1/phi as the diagonals, then

  A^n = V^-1 D^n V

where D^n is diagonal with phi^n and 1/phi^n as the diagonals. This gives exactly the above well known formula.

[btilly]

Yes, I have actually done the math. On the surface it looks very different, but a lot of the same numbers show up in intermediate calculations.

[shasta]

Well, both have computed the Fibonacci terms at a given level, so how different could it be? Here's my implementation:

    def fib_fast2(n):
        assert n >= 0
        a, b = 2, 0 # invariant: a,b are components of 2(phi^n)
        for bit in bits(n):
            a, b = (a*a + 5*b*b)>>1, a*b
            if bit: a, b = (a + 5*b)>>1, (a+b)>>1
        return b

It's almost identical runtime as the one in the article - a hair slower (15.32s vs. 16.17s to compute fib 10M). They're probably related by some well known relation between Fibonacci numbers.

[robinhouston]

That’s very nice! And a little rearrangement will bring it from three to two “big multiplications”, making it faster than my routine:

  def fib_ring2(n):
    assert n >= 0
    a, b = 2, 0 # invariant: phi^n = (a + b*sqrt(5)) / 2
    for bit in bits(n):
        ab = a*b
        a, b = (a+b)*((a+5*b)//2) - 3*ab, ab
        if bit: a, b = (a + 5*b)//2, (a+b)//2
    return b

[Xk]

That method is, for me at least, really really fast.

Eleven seconds to compute the 1,000,000th Fibonacci number. Compared with the naive method which takes 143 seconds.

Code below:

    def add(a, b): return (a[0]+b[0], a[1]+b[1])
    def sub(a, b): return (a[0]-b[0], a[1]-b[1])
    def divsq5(a): return (a[1], a[0]/5)
    def mul(a, b): return (a[0]*b[0]+5*a[1]*b[1],a[0]*b[1]+a[1]*b[0])
    def subi(x, a): return (x-a[0],-a[1])
    def pow(b, e):
        r = (1, 0)
        while e > 0:
            if e & 1 == 1:
                result = mul(r, b)
            e = e >> 1
            b = mul(b, b)
        return r

    twophi = (1,1)

    def fib0(n):
        return divsq5(sub(pow(twophi, n), pow(subi(2,twophi), n)))[0]>>n

    def fib1(n):
        a = 0
        b = 1
        while n > 0:
            (a,b) = (b, a+b)
            n -= 1
        return a

[shasta]

Two optimizations you can apply:

1. Given how your big integer code is probably implemented, you probably want to distribute the divisions by 2 into the intermediate calculations, rather just shifting by n at the very end.

2. Consider that (a, -b)(c, -d) = (ac + bd, -(bc + ad)). That means that the powers of (1,1) and (1,-1) are always related by negating the second component. Thus you only need to actually compute the power of (1,1) and double it.

[Xk]

Yeah, there are a whole lot of things I could do, but I just wanted a really quick program to see what the relative speeds were. And I have to say I was impressed.

[waterhouse]

I did this once for a Project Euler question. I called the function "awful-thing". It worked.

  ;this calculates (a+b√sqr)^n, using exponentiation by squaring
  (def awful-thing (n a b sqr)
    (afnwith (n n i a r b it 1 rt 0)
      (if (is n 0)
          (list it rt)
          (even n)
          (self (/ n 2)
                (+ (square i) (* sqr (square r)))
                (* 2 i r)
                it
                rt)
          (self (- n 1)
                i
                r
                (+ (* i it) (* sqr r rt))
                (+ (* i rt) (* r it))))))

  (def fib (n)
    (with (pow2 (expt 2 n)
           hh   (map - (awful-thing n 1 1 5) (awful-thing n 1 -1 5)))
      ;now hh = '(a b) where (fib n) = (a + b*root5)/(root5*pow2)
      ;a should = 0
      (/ (cadr hh)
         pow2)))

[shasta]

Thank you, that's a nice observation

[shasta]

Given the precision needed in the arithmetic, I've always wondered if the matrix exponentiating method wouldn't be faster. I've never been motivated enough to actually research it, test it, or think hard about it though...

[robinhouston]

I’d love to see an exact algorithm using the closed form formula. It’s obviously possible to do — though it seems fairly complicated — but would it really be faster? My instinct is that it would be slower: if anyone wants to prove me wrong, I’d be thrilled!

(The asymptotic complexity is surely the same, in any case.)

[ephemient]

http://stackoverflow.com/questions/4327846/calculating-fibon...

It's straightforward (if messy) to calculate two rationals exactly, then approximate √5 to within certain error bounds.

I wrote some code which does this, just as an experiment... and it's unbearably slow. Just try it.

  import Control.Arrow
  import Data.Ratio

  newtype Matrix2 a = Matrix2 (a, a, a, a) deriving (Show, Eq)
  instance (Num a) => Num (Matrix2 a) where
      Matrix2 (a, b, c, d) * Matrix2 (e, f, g, h) =
          Matrix2 (a*e+b*g, a*f+b*h, c*e+d*g, c*f+d*h)
      fromInteger x = let y = fromInteger x in Matrix2 (y, 0, 0, y)
  fib n = let Matrix2 (_, x, _, _) = Matrix2 (1, 1, 1, 0) ^ n in x

  binom n =
      scanl (\a (b, c)-> a*b `div` c) 1 $
      takeWhile ((/=) 0 . fst) $ iterate (pred *** succ) (n, 1)
  evens (x:_:xs) = x : evens xs
  evens xs = xs
  odds (_:x:xs) = x : odds xs
  odds _ = []
  iterate' f x = x : (iterate f $! f x)
  powers b = iterate (b *) 1
  esqrt e n = x where
      (_, x):_ = dropWhile ((<=) e . abs . uncurry (-)) $ zip trials (tail trials)
      trials = iterate (\x -> (x + n / x) / 2) n
  fib' n = (a, b) where
      d = 2 ^ n
      a = sum (zipWith (*) (odds $ binom n) (powers 5)) % d
      b = sum (zipWith (*) (evens $ binom n) (powers 5)) % d
  fib2 n = numerator r `div` denominator r where
      (a, b) = fib' n
      l = lcm (denominator a) (denominator a)
      r = a + esqrt (1 % l) (b * b / 5) + 1 % 2

[vog]

I think the biggest issue in this code in the naive implementation of "powers". It runs in O(b) while it should take at most O(log(b)).

Also, the multiplication might be a bottleneck, there has been a lot of research in this area. (the original article mentions some)

It might be a good idea to use an arbitrary precision library like GMP, which provides super-fast multiplication and thus power() algorithms.

(It also has an implementation of Fib() which outperforms both approaches here, but that's not my point.)

[ephemient]

I have not profiled the code, but I expect that the Newtonian approximation of sqrt is the slowest part.

It would surprise me if "powers" was a big problem relative to the rest. And no shortcuts: each of [1, b b², b³, b⁴, b⁵, …] (up to n) is used in the binomial expansion.

[ephemient]

Hmm, I screwed up and wrote iterate instead of iterate' where I intended. Not that it really makes a difference.

In any case, the average time to print the first 20 Fibonacci numbers is 0.4157ms with matrix exponentiation and 1.707063s with exact integer math (yes, the units are different) on a T400.

[kenjackson]

(The asymptotic complexity is surely the same, in any case.)

The closed form complexity is not the same (not sure how it could be). Here's the code written in standard most language psuedocode:

   fib(n) {
      double c = 1.6180339887;
      return round((power(c, n) - power(1-c, n))/sqrt(5));
   }

Note: I assume your point wasn't that computing the golden ratio exactly is... well ummm... time consuming.

[vog]

> The closed form complexity is not the same (not sure how it could be)

I disagree. Both variants are more similar than it may seem at a first glance. In essence, the question here is which of the following calculations is more efficient for big n:

a) the power function for arbitrary precision floating point

b) the power function for 2x2 matrices of arbitrary precision integers

Assuming that we can ignore the initial effort to calculate the golden ratio (can we?), Binet's formula is based on a), while the article's best algorith is an application of b).

I'm absolutely undecided which one works better, mostly because both kinds of power functions are working essentially the same way, with a complexity of O(log(n)).

Also, the precisions required for both cases seem to be quite similar.

[kenjackson]

Indeed. I thought the OP was referring to the recursive algorithm, not the final algorithm.

[vog]

Note that you can optimize this by removing the "(1-c)^n" part, because its absolute value is smaller than 0.5 for big n, so it can't influence the (rounded) result.

This trick is especially funny when performed using a hand calculator. For instance, just calculate the powers of the golden ratio and observe numbers like "xxxxx,0000yyy" "xxxxx,9999yyy", which become closer and closer to integers.

[robinhouston]

Actually I think that roughly was my point. :-)

You can’t get an exact answer in general by using fixed-size floating-point numbers, as your pseudocode does. The number of bits of precision you need to get an exact answer is going to depend on the input value, presumably linearly.

[kenjackson]

When you round you'll always get the exact integer answer, even using the precision of the golden ratio I gave in the code.

[robinhouston]

That’s really not true. You don’t have to take my word for it: try it with fib(1000), say.

The answer should be:

  43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875

[kenjackson]

Oh, the old n=1000 trick. Well played.

[akronim]

can't wait to break this out next time it comes up in an interview!