Earlier quoted context omitted.
And when you calculate this out, you internally wind up doing the same calculation as the matrix method.
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.
"The worst algorithm in the world?"
51–60 of 73 posts
Re: "The worst algorithm in the world?"
#52Very 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.
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) = (a c+5 b d,a d+b c) 2phi = (1,1) F(n) = (2phi^n - (2-2phi)^n)/(2^n*sqrt(5)) [edit] As ot said, the right word is "extended"
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 aRe: "The worst algorithm in the world?"
#53Earlier quoted context omitted.
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) = (a c+5 b d,a d+b c) 2phi = (1,1) F(n) = (2phi^n - (2-2phi)^n)/(2^n*sqrt(5)) [edit] As ot said, the right word is "extended"
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(…
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.
Re: "The worst algorithm in the world?"
#54It seems the author hasn't read SICP: http://mitpress.mit.edu/sicp/full-text/book/book-Z-H-11.html...
I don't think you've read the article, actually. The key point is that he is trying to calculate arbitrarily large fibonacci numbers, and so the addition is itself an O(n) operation. This means that even after memoization, the complexity is still O(n^2), and he uses a number of tricks to reduce that.
But it's a nice discussion in any case.
Re: "The worst algorithm in the world?"
#55Earlier quoted context omitted.
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(…
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 t…
Re: "The worst algorithm in the world?"
#56Earlier quoted context omitted.
When you round you'll always get the exact integer answer, even using the precision of the golden ratio I gave in the code.
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
Re: "The worst algorithm in the world?"
#57Earlier quoted context omitted.
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.
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.
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.Re: "The worst algorithm in the world?"
#58Earlier quoted context omitted.
And when you calculate this out, you internally wind up doing the same calculation as the matrix method.
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.
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.Re: "The worst algorithm in the world?"
#59Earlier quoted context omitted.
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.)
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) w…
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.
Re: "The worst algorithm in the world?"
#60Earlier quoted context omitted.
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) w…
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 h…
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.