"The worst algorithm in the world?"(bosker.wordpress.com) |
"The worst algorithm in the world?"(bosker.wordpress.com) |
(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"
I think "extended" is the right word (see http://en.wikipedia.org/wiki/Ring_extension)
It is commonly denoted as Z[sqrt(5)]
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.
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 ;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)))(The asymptotic complexity is surely the same, in any case.)
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 % 2The 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));
}
The code below consists of declaring that S is a sorted version of L as long as S is a permutation of L and S is sorted.
permutation_sort(L,S) :- permutation(L,S), sorted(S).
sorted([]).
sorted([_]).
sorted([X,Y|ZS]) :- X =< Y, sorted([Y|ZS]).
permutation([],[]).
permutation([X|XS],YS) :- permutation(XS,ZS),select(X,YS,ZS).Nonono, Bogosort is way worse than naive recursive fibonacci - the former doesn't even guarantee termination, recursive fibonacci still does.
If you want to calculate fibonacci numbers not as a misguided exercise in algorithms but actually efficiently, use an algebraic form: http://en.wikipedia.org/wiki/Fibonacci_number#Computation_by...
[1] http://en.wikipedia.org/wiki/Fibonacci_number#Closed-form_ex...
http://mitpress.mit.edu/sicp/full-text/book/book-Z-H-11.html...
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.
That's what I call a bad algorithm!
What's nice about it is that it is deterministic yet ridiculously slow.
It can also be really easily implemented in Prolog[1] where you simply define what a permutation and being sorted means. After that you just search for a sorted permutation.
[1] http://rosettacode.org/wiki/Sorting_algorithms/Permutation_s...
http://en.wikipedia.org/wiki/Quarantine_%28Greg_Egan_novel%2...
Since the Fibonacci numbers are defined by a linear recurrence, we can express the recurrence as a matrix, and it’s easy to verify that...
It's been way too long since my math minor for me to understand that.
I think it’s only the jargon that’s confusing you. As long as you can remember (or look up) the definition of matrix multiplication, then it genuinely is easy to verify.
[ fib(n-1) fib(n) ] x [0 1]
[ fib(n) fib(n+1) ] [1 1]
= [ fib(n) fib(n-1)+fib(n) ]
[ fib(n+1) fib(n)+fib(n+1) ]
= [ fib(n) fib(n+1) ]
[ fib(n+1) fib(n+2) ][Edit:] But there are some interesting data structures based on them: http://en.wikipedia.org/wiki/Fibonacci_heap
[Edited to add: thanks for the replies. I think I expressed myself poorly here. It’s not that I don’t understand the difference between “with probability 1” and “always”. What I mean is that I don’t understand why people sometimes make a big deal out of it, and say “that algorithm is not even guaranteed to terminate” as though that somehow means the algorithm is untrustworthy or useless. The trouble with the game “roll a die until you get 100 sixes in a row” is that it has an astronomical expected running time – not that it might _never_ end, but that it will almost certainly take a very long time.]
Unfortunately, it isn't. This is the reason why in mathematics, we say that "possibility 1" means that an event is "almost sure", which is different from a "sure" event.
For instance, you can play a game where you roll a dice over an over again. You win if you get 100 times a 6 in sequence. If you play this game without any time limit, your possibility to win is exactly 1, which means that you can be "almost sure" you'll win in the end. However, you can't be sure, because there is still the possibility that you don't get 100 times a 6 in an eternity.
Note that this only happens in infinite probability spaces. In finite spaces, "almost sure" and "sure" are equivalent.
Also note that the same holds for "probability 0" which means "almost impossible", not to be confused with "impossible".
Yes.
Particularly when done on real machines that use pseudo-random numbers that will repeat in finite time. For instance many use http://en.wikipedia.org/wiki/Mersenne_twister with a version that will repeat after 2^19937-1 steps. With a random array of 2500 elements, the odds are very good that bogosort will fall into a loop before it has successfully sorted the array.
http://en.wikipedia.org/wiki/Almost_surely
Bogosort almost surely terminates.
And the worst case of bogosort is not to terminate.
The probability that the game ends is 1, but it is not guaranteed to end. Yet no reasonable person will worry about this in practice.
Your example with 100 times 6 in a row is conceptually exactly the same, just the expected time until the game terminates is much, much larger.
Where is ZS defined? Bear in mind that I only took half a module's worth (6 modules in a year at my university) of Prolog and we really just used it for learning predicate logic.
Nevertheless, using Prolog (especially when my code worked!) blew my mind. In the assignment I had to write a program that mimicked part of an aircraft controller in that each aircraft in it's flightplan had to be separated from all other aircraft. Do you think Prolog will ever become popular (perhaps with the incoming Semantic web) or is it destined for academic work and system proofs only?
I don't believe Prolog will become any more popular than it is. The Japanese had some sort of a government project a few years ago, but I think that fell through. Now we have Erlang, Haskell and Java libraries that replicate most functionality, and then some.
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.
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.)
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.
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.
The actual notation for proportional asymptotic dominance of expectation is too convoluted for use, so people just don't use one. But to even leave out the word "expected" is just plain wrong.
AKA: (Ignoring the fact it takes more than one number to shuffle a deck). After the first cycle of 2^19937-1 you have ~1 in (2500!) that you repeat the initial position. After the second cycle it's ~2 in (2500!) etc. It may be possible that your odds of success are around (2^19937-2)/(2^19937-1).
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}}
A = V^-1 D V
A^n = V^-1 D^n V
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
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.
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.
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.
The answer should be:
43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875 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