Re: To understand recursion, one must first understand recursion (was Re: apple iie overheat question)

[Scott Hemphill <hemphill@hemphills.net>]

"anoneds@netscape.net" <anoneds@netscape.net> writes:

> Scott Hemphill wrote:
> <snip>
> 
> >
> > I was hoping that somebody would bite on my direct calculation for
> the
> > Fibonacci sequence, wondering where my "magic constants" came from,
> and
> > why it works.  I guess most people either already know, or don't
> care.
> >
> 
> Hmm. Well, I care! :)
> Please, elaborate? I haven't had a chance to even try it, yet!

The first step is to find a closed form for the Fibonacci sequence.  Let's
call any sequence starts with two numbers a0 and a1 and then is further
defined by a2=a1+a0, a3=a2+a1, etc. an "additive sequence".  If you
multiply an additive sequence {a0,a1,a2,...} by a constant, then the
result is an additive sequence.  If you add two different additive
sequences {a0,a1,a2,...} and {b0,b2,b2,...} together, the result is also
an additive sequence.  Suppose we know a closed form for calculating two
different additive sequences.  Then we can "mix" them by multiplying the
first sequence by one constant, the second sequence by another constant,
and adding that together.  If we are careful about choosing the constants,
we can arrange that the first number of the resulting additive sequence is
zero and the second number is one, and the sequence will therefore be the
Fibonacci sequence.

Here's how this works:

Look at the sequence {1,g,g^2,g^3,...}.  If g^2=g+1, then the sequence
is additive because (multiplying both sides of the equation by g) g^3=g^2+g.
If you keep multiplying both sides of the equation by you wind up with
g^4=g^3+g^2, etc.

So what does "g" have to be to make this sequence additive?  The quadratic
equation g^2-g-1=0 has two roots, (1+sqrt(5))/2 and (1-sqrt(5))/2.
The first solution is the "golden ratio", 1.618... and the second is
-0.618....  From this point on, I'll use "g" to refer to the golden ratio,
and "1-g" to refer to the other solution of g^2-g-1=0.

So we have just come up with closed forms for two additive sequences:

  {1,g,g^2,...} and {1,1-g,(1-g)^2,...}.

We want to mix these by multiplying each by a constant and adding that
together in such a way as to form the Fibonacci sequence.

  c {1,g,g^2,...} + d {1,1-g,(1-g)^2,...} = {0,1,1,2,...}

If the first number of the sequence is going to be zero, then c + d = 0,
or d = -c.  Then we have:

  c {1,g,g^2,...} - c {1,1-g,(1-g)^2,...} = {0,1,1,2,...}

If the second number of the sequence is going to be one, then c*g-c*(1-g)=1,
or c=1/(2*g-1)=1/sqrt(5).

So we now have a closed form for the nth term of the Fibonacci sequence:

  (g^n-(1-g)^n)/sqrt(5)

Now look at what happens as n gets bigger and bigger.  "g" is bigger than
one, so g^n keeps growing.  However (1-g) has absolute value less than one,
so (1-g)^n keeps getting closer and closer to zero.  If n is large enough,
we can forget about the (1-g)^n term completely.  So lets look at the
sequence defined by g^n/sqrt(5):

  {0.447214, 0.723607, 1.17082, 1.89443, 3.06525, 4.95967, 8.02492, 
   12.9846, 21.0095, 33.9941, 55.0036, ...}

We can see that it gets closer and closer to the Fibonacci sequence.  The
contribution of the (1-g)^n is small enough that we can ignore it completely
for all n>=0.  All we have to do is to round the sequence to the nearest
integer, and it will be correct for all n.

  {0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, ...}

So the equation implementing the closed form solution looks something like:

  round(pow(g,n)/sqrt(5)), where g = 1.6180339887498948482045868343656
                       and sqrt(5) = 2.2360679774997896964091736687313

round is not a standard C function, but can be implemented as floor(0.5+ ...)

Since pow(g,n) might overflow in cases where pow(g,n)/sqrt(5) is still
within the representable range of numbers, and who knows, exp might be
faster than pow, I transformed pow(g,n)/sqrt(5):

  pow(g,n)/sqrt(5) = exp(log(pow(g,n)/sqrt(5))
                   = exp(n*log(g)-log(sqrt(5))

This is nice because the two logarithms are both constants.

  log(g) = 0.48121182505960344749775891342437
  log(sqrt(5)) = 0.80471895621705018730037966661309

So that should pretty much how I came up with this implementation

>  #include <math.h>
>  double F(int n)
>  {
>    const double c1 = 0.48121182505960344749775891342437; // magic constants 
>    const double c2 = 0.80471895621705018730037966661309; // :-)
>  // let n < 0 fall through.  function will return zero
>    return floor(0.5 + exp(n*c1-c2));
>  }

Scott
-- 
Scott Hemphill	hemphill@alumni.caltech.edu
"This isn't flying.  This is falling, with style."  -- Buzz Lightyear