Faster Recursive Fibonacci Numbers

[RJB]

I noticed some discussion of recursion..... the trick is to find a
formula where the arguments are divided, not decremented.
I've had a "divide-and-conquer" recursion for the Fibonacci numbers
for a couple of years in C++ but just for fun rewrote it
in Python. It was easy. Enjoy. And tell me how I can improve it!

def fibo(n):
"""A Faster recursive Fibonaci function
Use a formula from Knuth Vol 1 page 80, section 1.2.8:
If F[n] is the n'th Fibonaci number then
F[n+m] = F[m]*F[n+1] + F[m-1]*F[n].
First set m = n+1
F[ 2*n+1 ] = F[n+1]**2 + F[n]*2.

Then put m = n in Knuth's formula,
F[ 2*n ] = F[n]*F[n+1] + F[n-1]* F[n],
and replace F[n+1] by F[n]+F[n-1],
F[ 2*n ] = F[n]*(F[n] + 2*F[n-1]).
"""
if n<=0:
return 0
elif n<=2:
return 1
elif n%2==0:
half=n//2
f1=fibo(half)
f2=fibo(half-1)
return f1*(f1+2*f2)
else:
nearhalf=(n-1)//2
f1=fibo(nearhalf+1)
f2=fibo(nearhalf)
return f1*f1 + f2*f2


RJB the Lurker
http://www.csci.csusb.edu/dick/cs320/lab/10.html

 

[Ian Kelly]

I noticed some discussion of recursion..... the trick is to find a
formula where the arguments are divided, not decremented.
I've had a "divide-and-conquer" recursion for the Fibonacci numbers
for a couple of years in C++ but just for fun rewrote it
in Python.  It was easy.  Enjoy.  And tell me how I can improve it!

def fibo(n):
       """A Faster recursive Fibonaci function
Use a formula from Knuth Vol 1 page 80, section 1.2.8:
          If F[n] is the n'th Fibonaci number then
                  F[n+m] = F[m]*F[n+1] + F[m-1]*F[n].
 First set m = n+1
  F[ 2*n+1 ] = F[n+1]**2 + F[n]*2.

 Then put m = n in Knuth's formula,
          F[ 2*n ] = F[n]*F[n+1] + F[n-1]* F[n],
  and replace F[n+1] by F[n]+F[n-1],
          F[ 2*n ] = F[n]*(F[n] + 2*F[n-1]).
"""
       if n<=0:
               return 0
       elif n<=2:
               return 1
       elif n%2==0:
               half=n//2
               f1=fibo(half)
               f2=fibo(half-1)
               return f1*(f1+2*f2)
       else:
               nearhalf=(n-1)//2
               f1=fibo(nearhalf+1)
               f2=fibo(nearhalf)
               return f1*f1 + f2*f2

Thanks for posting! Actually, it looks like this is the same O(n)
algorithm that rusi posted. There was also a O(log n) algorithm
discussed that is based on vector math. You might want to take a
look.

Cheers,
Ian

 

[rusi]

I noticed some discussion of recursion..... the trick is to find a
formula where the arguments are divided, not decremented.
I've had a "divide-and-conquer" recursion for the Fibonacci numbers
for a couple of years in C++ but just for fun rewrote it
in Python.  It was easy.  Enjoy.  And tell me how I can improve it!

def fibo(n):
        """A Faster recursive Fibonaci function
Use a formula from Knuth Vol 1 page 80, section 1.2.8:
           If F[n] is the n'th Fibonaci number then
                   F[n+m] = F[m]*F[n+1] + F[m-1]*F[n].
  First set m = n+1
   F[ 2*n+1 ] = F[n+1]**2 + F[n]*2.

  Then put m = n in Knuth's formula,
           F[ 2*n ] = F[n]*F[n+1] + F[n-1]* F[n],
   and replace F[n+1] by F[n]+F[n-1],
           F[ 2*n ] = F[n]*(F[n] + 2*F[n-1]).
"""
        if n<=0:
                return 0
        elif n<=2:
                return 1
        elif n%2==0:
                half=n//2
                f1=fibo(half)
                f2=fibo(half-1)
                return f1*(f1+2*f2)
        else:
                nearhalf=(n-1)//2
                f1=fibo(nearhalf+1)
                f2=fibo(nearhalf)
                return f1*f1 + f2*f2

RJB the Lurkerhttp://www.csci.csusb.edu/dick/cs320/lab/10.html

 

[rusi]

I noticed some discussion of recursion..... the trick is to find a
formula where the arguments are divided, not decremented.
I've had a "divide-and-conquer" recursion for the Fibonacci numbers
for a couple of years in C++ but just for fun rewrote it
in Python.  It was easy.  Enjoy.  And tell me how I can improve it!

def fibo(n):
        """A Faster recursive Fibonaci function
Use a formula from Knuth Vol 1 page 80, section 1.2.8:
           If F[n] is the n'th Fibonaci number then
                   F[n+m] = F[m]*F[n+1] + F[m-1]*F[n].
  First set m = n+1
   F[ 2*n+1 ] = F[n+1]**2 + F[n]*2.

  Then put m = n in Knuth's formula,
           F[ 2*n ] = F[n]*F[n+1] + F[n-1]* F[n],
   and replace F[n+1] by F[n]+F[n-1],
           F[ 2*n ] = F[n]*(F[n] + 2*F[n-1]).
"""
        if n<=0:
                return 0
        elif n<=2:
                return 1
        elif n%2==0:
                half=n//2
                f1=fibo(half)
                f2=fibo(half-1)
                return f1*(f1+2*f2)
        else:
                nearhalf=(n-1)//2
                f1=fibo(nearhalf+1)
                f2=fibo(nearhalf)
                return f1*f1 + f2*f2

RJB the Lurkerhttp://www.csci.csusb.edu/dick/cs320/lab/10.html

-------------------------------------------------------------
Its an interesting problem and you are 75% there.
You see the halving gives you logarithmic behavior and the double
calls give exponential behavior.

So how to get rid of double calls? Its quite simple: Just define your
function in terms of return pairs of adjacent pairs ie (fib(n), fib(n
+1)) for some n rather then a single number fib(n)

Here's a straightforward linear function:

def fp(n): #fibpair
if n==1:
return (1,1)
else:
a,b = fp(n-1)
return (b, a+b)

def fib(n):
a,b = fp(n)
return a

---------------
Now use this (pairing) idea with your (halving) identities and you
should get a logarithmic algo.

[If you cant do it ask again but yes its fun to work out so do
try ]

 

[Jussi Piitulainen]

or O(1):

φ = (1 + sqrt(5)) / 2
def fib(n):
numerator = (φ**n) - (1 - φ)**n
denominator = sqrt(5)
return round(numerator/denominator)

Testing indicates that it's faster somewhere around 7 or so.

And increasingly inaccurate from 71 on.

 

[geremy condra]

And increasingly inaccurate from 71 on.

Yup. That's floating point for you. For larger values you could just
add a linear search at the bottom using the 5f**2 +/- 4 rule, which
would still be quite fast out to about 10 times that. The decimal
module gets you a tiny bit further, and after that it's time to just
use Dijkstra's, like rusi suggested. In any event, I still think this
formulation is the most fun .

Geremy Condra

 

[Wolfram Hinderer]

Yup. That's floating point for you. For larger values you could just
add a linear search at the bottom using the 5f**2 +/- 4 rule, which
would still be quite fast out to about 10 times that. The decimal
module gets you a tiny bit further, and after that it's time to just
use Dijkstra's, like rusi suggested. In any event, I still think this
formulation is the most fun .

I think you can write it even more funny

def fib(n):
return round(((.5 + .5 * 5 ** .5) ** n - (.5 - .5 * 5 ** .5) **
n) * 5 ** -.5)

;-)

 

[geremy condra]

I think you can write it even more funny

def fib(n):
   return round(((.5 + .5 * 5 ** .5) ** n -  (.5 - .5 * 5 ** .5) **
n) * 5 ** -.5)

;-)

Ok, that's amusing. It does hide the interaction with the golden ratio
though, which is what I find so fascinating about the earlier one.

Geremy Condra

 

[rusi]

I think you can write it even more funny

def fib(n):
    return round(((.5 + .5 * 5 ** .5) ** n -  (.5 - .5 * 5 ** .5) **
n) * 5 ** -.5)

;-)

VOW!
Tour de Force - Thanks
[I am going to trouble some class of students with that ]

 

[Steven D'Aprano]

or O(1):

φ = (1 + sqrt(5)) / 2
def fib(n):
numerator = (φ**n) - (1 - φ)**n
denominator = sqrt(5)
return round(numerator/denominator)

I'd just like to point out that, strictly speaking, it's only O(1) if you
assume that exponentiation is O(1). Computer scientists often like to
make this simplifying assumption, and it might even be true for floats,
but for long ints and any numeric data types with unlimited precision, it
won't be.

 

[Chris Angelico]

I'd just like to point out that, strictly speaking, it's only O(1) if you
assume that exponentiation is O(1). Computer scientists often like to
make this simplifying assumption, and it might even be true for floats,
but for long ints and any numeric data types with unlimited precision, it
won't be.

Python doesn't have arbitrary precision non-integers, does it? So this
is going to be done with floats.

Chris Angelico

 

[geremy condra]

I'd just like to point out that, strictly speaking, it's only O(1) if you
assume that exponentiation is O(1). Computer scientists often like to
make this simplifying assumption, and it might even be true for floats,
but for long ints and any numeric data types with unlimited precision, it
won't be.

Great point, and something I should have paid attention to. Thanks.

Geremy Condra

 

[Steven D'Aprano]

Python doesn't have arbitrary precision non-integers, does it? So this
is going to be done with floats.

Sure, which is why the above fib() function will become increasing
inaccurate beyond some given n, by memory about n=71 or so. Er, at least
the fib() function that *was* above until you deleted most of it

If you want an *accurate* fib() function using exponentiation of φ, you
need arbitrary precision non-integers.

Nevertheless, at some point you will hit the limit of floats, which
thanks to the exponential growth of the Fibonacci sequence won't take
that long: it takes roughly 1475 iterations to exceed the range of Python
floats.

 

[Chris Angelico]

... until you deleted most of it

Minimalist quoting practice!

If you want an *accurate* fib() function using exponentiation of ö, you
need arbitrary precision non-integers.

I believe the 'bc' command-line calculator can do a-p non-i, and I
know REXX can, but it seems to be quite an unusual thing. Is it that
much harder than a-p integer that it's just not worthwhile? It seems
strange to smoothly slide from native integer to long integer and just
keep on going, and yet to be unable to do the same if there's a
fractional part on it.

Chris Angelico

 

[harrismh777]

I believe the 'bc' command-line calculator can do a-p non-i, and I
know REXX can

Yes, bc is wonderful in this regard. Actually, bc does this sort of
thing in 'circles' around Python. This is one of Python's weaknesses for
some problem solving... no arbitrary precision. And its not just that bc
does arbitrary precision--- its that it does it fast!

Actually, it should be relatively easy to incorporate parts of bc into
Python as C extensions. On the other hand, when needing specialized math
work from bc, its probably just better to use bc and leave Python alone.

On the other hand, most of the time (and I mean 99.999% of the time)
floats are going to work just fine... usually folks don't even need
doubles....

With fifteen or twenty digits of PI we can calculate the circumference
of the visible universe to within the width of a proton... er, I mean
hadron... and you know what... how many folks need to do that anyway??


Don't get me wrong... I absolutely love playing around with bignums, but
then, I'm a math geek... ;-)



kind regards,
m harris

 

[Chris Angelico]

Actually, it should be relatively easy to incorporate parts of bc into
Python as C extensions. On the other hand, when needing specialized math
work from bc, its probably just better to use bc and leave Python alone.

If someone has time to kill (as if!), it'd be awesome to get a new
numeric type that uses bc's code; any other numeric type (int, long,
float) could autopromote to it, removing the dilemma of which to
promote out of long and float. Hmm... Python 4.0, 'bc' is the new
default integer and everything else is a performance optimization?
Heh.

On the other hand, most of the time (and I mean 99.999% of the time) floats
are going to work just fine... usually folks don't even need doubles....

99.9% of the time int will work fine, too. Most people don't need
arbitrary precision OR floating point.

Don't get me wrong... I absolutely love playing around with bignums, but
then, I'm a math geek...    ;-)

Absolutely. Bring on the geekiness.

I've used bignums for things other than straight arithmetic, actually.
In REXX, where everything is a string, I've done some fascinating (and
completely useless) analyses that involve taking internal digits from
a number, performing arithmetic on them, getting huge numbers back,
and then searching for substrings (ie digit strings) in the result.
What's the lowest power of 2 that has 5 consecutive digits in it? All
10 digits? (That is, it has '0123456789' somewhere in its decimal
representation.)

Like I said, completely useless... but how many of you immediately
pondered which language to implement the search in?

Chris Angelico

 

[Ian Kelly]

Sure, which is why the above fib() function will become increasing
inaccurate beyond some given n, by memory about n=71 or so. Er, at least
the fib() function that *was* above until you deleted most of it

If you want an *accurate* fib() function using exponentiation of ö, you
need arbitrary precision non-integers.

This seems to work for arbitrary n, although I only tested it out to
about n=2000.

from decimal import Decimal, localcontext
from math import sqrt

def fib(n):
if n <= 70:
return fib_float(n)
else:
return fib_decimal(n)

phi_float = (1 + sqrt(5)) / 2
def fib_float(n):
numerator = (phi_float ** n) - (1 - phi_float) ** n
denominator = sqrt(5)
return int(round(numerator / denominator))

def fib_decimal(n):
with localcontext() as ctx:
ctx.prec = n // 4 + 1
sqrt_5 = Decimal('5').sqrt()
phi = (1 + sqrt_5) / 2
numerator = (phi ** n) - (1 - phi) ** n
return int((numerator / sqrt_5).to_integral_value())

Cheers,
Ian

 

[Ian Kelly]

def fib_decimal(n):
   with localcontext() as ctx:
       ctx.prec = n // 4 + 1
       sqrt_5 = Decimal('5').sqrt()
       phi = (1 + sqrt_5) / 2
       numerator = (phi ** n) - (1 - phi) ** n
       return int((numerator / sqrt_5).to_integral_value())

Note that I believe this is O(n) since the necessary precision grows
linearly with n.

 

[SigmundV]

There is a nice matrix representation of consecutive Fibonacci
numbers: [[1, 1], [1, 0]] ** n = [[F_n+1, F_n], [F_n, F_n-1]]. Using
the third party mpmath module, which uses arbitrary precision floating
point arithmetic, we can calculate the n'th Fibonacci number for an
arbitrary n as follows:

import mpmath
A = mpmath.matrix([[1, 1], [1, 0]])
F = A ** n

The n'th Fibonacci number is then found as the elements [0, 1] and [1,
0] in the matrix F.

This is more expensive than the formula involving the golden ratio,
but I like the compact representation.

 

[Steven D'Aprano]

If someone has time to kill (as if!), it'd be awesome to get a new
numeric type that uses bc's code; any other numeric type (int, long,
float) could autopromote to it, removing the dilemma of which to promote
out of long and float. Hmm... Python 4.0, 'bc' is the new default
integer and everything else is a performance optimization? Heh.

The problem is, it isn't *just* a performance optimization, there is a
semantic difference too. Consider:
True

But using something with more precision:
False


So you get different behaviour between floats and arbitrary precision
numbers.