3 ms·

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 fa

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

kenjackson 15y ago
(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.
- robinhouston 15y ago
  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 15y ago
    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 15y ago
      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 15y ago
        Oh, the old n=1000 trick. Well played.
- vog 15y ago
  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.
- akronim 15y ago
  can't wait to break this out next time it comes up in an interview!
- vog 15y ago
  > 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 15y ago
    Indeed. I thought the OP was referring to the recursive algorithm, not the final algorithm.
ephemient 15y ago
http://stackoverflow.com/questions/4327846/calculating-fibonacci/4328194#4328194 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 15y ago
  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 15y ago
    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.

kenjackson 15y ago

(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.

robinhouston 15y ago
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 15y ago
  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 15y ago
    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 15y ago
      Oh, the old n=1000 trick. Well played.
vog 15y ago
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.
akronim 15y ago
can't wait to break this out next time it comes up in an interview!
vog 15y ago
> 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 15y ago
  Indeed. I thought the OP was referring to the recursive algorithm, not the final algorithm.

ephemient 15y ago

http://stackoverflow.com/questions/4327846/calculating-fibonacci/4328194#4328194 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 15y ago
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 15y ago
  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.