5 ms·

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

by stralep 15y ago

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"

vog 15y ago
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.
- deleted 15y ago
  [deleted]
ot 15y ago
> 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 http://en.wikipedia.org/wiki/Ring_extension) It is commonly denoted as Z[sqrt(5)]
shasta 15y ago
Thank you, that's a nice observation
btilly 15y ago
And when you calculate this out, you internally wind up doing the same calculation as the matrix method.
- deleted 15y ago
  [deleted]
- vog 15y ago
  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.
  - btilly 15y ago
    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 15y ago
      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 15y ago
        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
        shasta 15y ago
        Oops, I should have found that. I looked and thought there was some reason it wouldn't work. Interesting post, thanks for the fun.
Xk 15y ago
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 15y ago
  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 15y ago
    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 15y ago
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)))

by stralep 15y ago

vog 15y ago

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.

deleted 15y ago
[deleted]

ot 15y ago

> 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 http://en.wikipedia.org/wiki/Ring_extension) It is commonly denoted as Z[sqrt(5)]

shasta 15y ago

Thank you, that's a nice observation

btilly 15y ago

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

deleted 15y ago
[deleted]
vog 15y ago
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.
- btilly 15y ago
  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 15y ago
    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 15y ago
      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
      - shasta 15y ago
        Oops, I should have found that. I looked and thought there was some reason it wouldn't work. Interesting post, thanks for the fun.

Xk 15y ago

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 15y ago
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 15y ago
  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 15y ago

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