Fibonacci program makes matlab go busy forever

This is the function I ran for fibonacci numbers. It works well for small numbers but when you reach large ones matlab just turns "busy" and freezes up. I must use this program for the 10th, 50th, and 100th fibonacci numbers.

function f = fibonacci(n)

if(n==1)

      f=1;

elseif(n==2)

      f=1;
else
      f = fibonacci(n-1) + fibonacci(n-2);
end

Answer 1

[John D'Errico]

And, why are you surprised that it fails? :)

I go into a fair amount of detail on this in my Fibonacci submission on the File Exchange. But the gist of it is that recursion is a flat out terrible way to compute Fibonacci numbers. My guess is this is why you were given that homework problem.

Lets think about what happens here. the first call, fibonacci(100), then calls itself twice, trying to compute fibonacci(99) and fibonacci(98).

No obvious problem there. But now that call to fibonacci(99) calls itself twice, trying to compute f_98 and f_97. Each successive call spawns additional calls, and in fact, many of those calls are repeated. So there are two calls to get f_98, 3 times when it will evaluate f_97. Think of those recursive calls as a tree, with many branches that intertwine, crossing. But the problem is, fibonacci does not know that it has, someplace else in the tree, needed the value of a given fibonacci number.

How many calls are there? Here are the number of calls for each specific fibonacci number in that tree. See if you recognize the pattern:

f_100 - 1
f_99 - 1
f_98 - 2
f_97 - 3
f_96 - 5
f_95 - 8
f_94 - 13
...

Interesting. So the number of calls actually grows exactly as does the fibonacci sequence, but backwards. How many calls will there be for fibonacci(1)?

I can compute this, because I'm using my own code for Fibonacci numbers, that is quite efficient. As well, my own code uses my vpi tools, so the value given is exact.

fibonacci(100)
ans =
    354224848179261915075

DON'T TRY THAT WITH YOUR OWN CODE! My version uses a O(log(n)) scheme to compute the nth Fibonacci number.

And that is just the number of times that fibonacci(1) will be evaluated by your recursive code, so roughly 3e20 recursive calls!

vpi2english(fibonacci(100))
ans =
three hundred fifty four quintillion, two hundred twenty four quadrillion, eight hundred forty eight trillion, one hundred seventy nine billion, two hundred sixty one million, nine hundred fifteen thousand, seventy five

That is one big number. But the total number of recursive calls is even worse.

sum(fibonacci(1:100))
ans =
    927372692193078999175

Yes, that is 927 quintillion calls (with a few to spare.)

So there are good ways to compute large Fibonacci numbers, and there are obscenely bad ways. Straight recursion (as you have) is an example of an obscenely bad way. Again, my guess is this assignment was given to make a point.

Again, there are good ways to do that computation.

fibonacci(10000)
ans =
    33644764876431783266621612005107543310302148460680063906564769974680
081442166662368155595513633734025582065332680836159373734790483865268263
040892463056431887354544369559827491606602099884183933864652731300088830
269235673613135117579297437854413752130520504347701602264758318906527890
855154366159582987279682987510631200575428783453215515103870818298969791
613127856265033195487140214287532698187962046936097879900350962302291026
368131493195275630227837628441540360584402572114334961180023091208287046
088923962328835461505776583271252546093591128203925285393434620904245248
929403901706233888991085841065183173360437470737908552631764325733993712
871937587746897479926305837065742830161637408969178426378624212835258112
820516370298089332099905707920064367426202389783111470054074998459250360
633560933883831923386783056136435351892133279732908133732642652633989763
922723407882928177953580570993691049175470808931841056146322338217465637
321248226383092103297701648054726243842374862411453093812206564914032751
086643394517512161526545361333111314042436854805106765843493523836959653
428071768775328348234345557366719731392746273629108210679280784718035329
131176778924659089938635459327894523777674406192240337638674004021330343
297496902028328145933418826817683893072003634795623117103101291953169794
607632737589253530772552375943788434504067715555779056450443016640119462
580972216729758615026968443146952034614932291105970676243268515992834709
891284706740862008587135016260312071903172086094081298321581077282076353
186624611278245537208532365305775956430072517744315051539600905168603220
349163222640885248852433158051534849622434848299380905070483482449327453
732624567755879089187190803662058009594743150052402532709746995318770724
376825907419939632265984147498193609285223945039707165443156421328157688
908058783183404917434556270520223564846495196112460268313970975069382648
706613264507665074611512677522748621598642530711298441182622661057163515
069260029861704945425047491378115154139941550671256271197133252763631939
606902895650288268608362241082050562430701794976171121233066073310059947
366875

timeit(@() fibonacci(10000))
ans =
          0.027053912189

So 0.027 seconds using the O(log(n)) scheme, and it is exact, down to the last digit. Had I tried to use the basic recursive scheme on this number, it would have brought the world's largest supercomputer to its knees, for the life of the universe. Sometimes the algorithm you use matters. A lot.

[Andrew Newell]

Rather than compete with such an excellent discussion, I'd just like to add that replacing recursion by a simple loop will greatly improve the speed.

[John D'Errico]

Yes, very true. And I admit that my all time favorite Fibonacci scheme is this one:

F = @(n) det(eye(n) + i*diag(ones(n-1,1),1) + i*diag(ones(n-1,1),-1));

F(1)
ans =
   1
F(2)
ans =
   2
F(3)
ans =
   3
F(4)
ans =
   5
F(5)
ans =
   8
F(6)
ans =
  13
...

That one was worth a paper or two because of where it takes you.

Regardless, even with a boring loop you still have the problem that beyond a limited number, the computation will still fail unless you use a higher precision type.

As it turns out, a double will be limited to about F_76 or so, and F_100 will exceed even the capabilities of uint64.

log2(fibonacci(76))
ans =
        51.6014213884832

log2(fibonacci(93))
ans =
        63.4035339202037

log2(fibonacci(100))
ans =
        68.2632273156181

[Andrew Newell]

What do you mean by "was worth" a paper or two? Have they been published?

[John D'Errico]

This is one of them, in that great read - The Fibonacci Quarterly. I'm sure we all keep a copy at hand.

http://www.fq.math.ca/Scanned/41-1/cahill.pdf

It has been a while since we did this, but the general idea is you can write that determinant recursively, as is possible with determinants. When you do so, lo and behold, it generates a similar recurrence relation to that which generates the Fibonacci numbers. Then you might notice a very interesting feature about the eigenvalues of those matrices. They look like the roots of a set of Chebyshev polynomials, which, to no surprise, ALSO can be generated recursively, using a similar recurrence relation. A nice thing about those polynomials is you can write down their roots simply, so since the determinant will be simply derived from the product of the eigenvalues, we come full circle to be able to derive a formula for the nth Fibonacci number, in a quite different looking form.

Another one here.

[Andrew Newell]

I seem to have misplaced that issue of the Fibonacci Quarterly, so I appreciate the link. Seriously, that is a very intriguing result!

[John D'Errico]

If you think about it, the 2 term recurrence relation for orthogonal polynomials looks amazingly like the relation for Fibonacci numbers, IF you choose the right parameters. We tripped over the idea when we heard about the determinant trick for Fibonacci numbers. Then it just leads you down a natural path to the solution, as long as you recognize the connections.

This is one of the fruits of a sufficiently broad education in mathematics, that one starts to see connections between disparate areas. Here the trick is to connect orthogonal polynomials, determinants, and Fibonacci numbers.

Answer 2

[Jos (10584)]

The Nth matrix power of the square matrix [1 1 ; 1 0] will give you three Fibonacci numbers at once, namely the (N+1)-th, the N-th (twice) and the (N-1)-th.

[1 1 ; 1 0]^10
% -> 89    55
%    55    34

Answer 3

[Roger Stafford]

You can also compute the n-th term in the Fibonacci series with this formula:

   f(n) = round((((1+sqrt(5))/2)^n-((1-sqrt(5))/2)^n)/sqrt(5));

Using double precision it is valid up to about the 70-th term at which point the round-off error becomes too large for 'round' to correct it properly. With the symbolic toolbox it would be valid to however many decimal places the user sets.