Tag Archives: fibonacci numbers

Project Euler in R: Problem 25

Problem 25:  What is the first term in the Fibonacci sequence to contain 1000 digits?

Commentary:  Well, I certainly hope you paid attention to all the boring crap in the solution to Problem 2, because you’re going to need it to understand my solution for Problem 25.

First, I’ll describe my first attempted solution, which was a “computery” approach.  My original plan was to use the closed form (discussed in the link above) to the Fibonacci Numbers, passing off all the “heavy lifting” to GNU bc.  This “works”, but is UNBELIEVABLY slow (I think I let it run for about 15 minutes before I decided to try to math my way out of things).  This problem is a good reminder (though we’ll definitely see others) of the inarguable fact that mathematicians are still kings of earth–or at least that the Project Euler people like to occasionally throw a bone to math-types instead of letting the computer science nerds have all the fun.

Ok, so, remember that weird “continued fraction” expansion of \phi that I gave in the solution to Problem 2?  If you stare at it for a few seconds, then it should be obvious that

\phi = 1+\displaystyle\frac{1}{\phi}

and hence

1-\phi = \displaystyle\frac{1}{\phi}

While we’re here, it’s worth pointing out that if you define \phi in terms of that continued fraction expansion, then you can deduce from that definition the numeric value of \phi , since from the first equation above, if we multiply both sides by \phi, then rearrange things and get

\phi^2-\phi-1=0

Using the good old quadratic formula gives two solutions (or in other words, two possibilities for \phi), namely \displaystyle\frac{1\pm\sqrt{5}}{2}.  One of these two possibilities is positive and the other negative.  From there it isn’t hard to see which is which.

Ok, so back to Project Euler.  With the above observation regarding 1-\phi, we can rewrite the closed form expression of the Fibonacci Numbers (using the Project Euler convention that F_1=F_2=1) as

F_n = \displaystyle\frac{\phi^n-\left(\frac{1}{\phi}\right)^n}{\sqrt{5}}

Then for sufficiently large n (and thinking about the problem at hand for 3 seconds should convince you that the n we’re after should do the trick), we have

F_n\approx\displaystyle\frac{\phi^n}{\sqrt{5}}

since of course \displaystyle\left(\frac{1}{\phi}\right)^n will be quite small.  It follows that we need only find the smallest natural number n satisfying

F_n = \displaystyle\frac{\phi^n}{\sqrt{5}}>10^{999}

And since \log is an increasing function, this is equivalent to needing to find the smallest n satisfying

\log_{10}(F_n)>999

The problem is now trivial.  Note that from the above observation regarding \displaystyle\left(\frac{1}{\phi}\right)^n and from basic properties of logarithms:

\log_{10}(F_n) = n\log_{10}(\phi)-\frac{1}{2}\log_{10}(5)

Whence, we wish to find the smallest natural number n satisfying

n>\displaystyle\frac{999+\frac{1}{2}\log_{10}(5)}{\log_{10}\left(\phi\right)}\approx 4781.859

It follows that n=4782.

Who says a masters in math is useless?

R Code:

ElapsedTime <- system.time({
##########################
phi<-(1+sqrt(5))/2
answer<-ceiling((999+1/2*log(5, 10))/log(phi, 10))
##########################
})[3]
ElapsedMins <- floor(ElapsedTime/60)
ElapsedSecs <- (ElapsedTime-ElapsedMins*60)
cat(sprintf("\nThe answer is:  %d\nTotal elapsed time:  %d minutes and %f seconds\n", answer, ElapsedMins, ElapsedSecs))


Output:
The answer is:  4782
Total elapsed time:  0 minutes and 0.000000 seconds


Project Euler in R: Problem 2

Problem 2: By considering the terms in the Fibonacci sequence whose values do not exceed four million, find the sum of the even-valued terms.

Commentary:  This problem again is not very tricky.  It is worth mentioning here that I am a mathematician by training, and I’m about to math the hell out.  If you’re not into that, you are encouraged to skip this entire section.

The Fibonacci numbers are governed by the recurrence relation:

F_n = F_{n-1} + F_{n-2}

with F_0 = F_1 = 1.  What’s not stated in the Project Euler explanation (which will become debatably useful in a later problem), is that there is a “closed form” to the Fibonacci Numbers, given by:

F_n = \displaystyle\frac{\left(\frac{1+\sqrt{5}}{2}\right)^{n+1} - \left(\frac{1-\sqrt{5}}{2}\right)^{n+1}}{\sqrt{5}}

(where the exponents “n+1” would be replaced by n if you index the sequence via F_1=F_2=1).  Demonstrating this is not particularly difficult if you’re at all comfortable with simple calculus (which now qualifies as high school mathematics).  To do this, we’ll use a generating function (there are other approaches to get this “closed form”, but this, in my opinion, is by far the simplest).  Let’s first construct a funny thing called S:

S:= \displaystyle\sum_{n=0}^\infty F_n x^n

You might be wondering if this is a sensible definition, i.e., does it even make sense to say that this arcane construction S even exists.  If so, congratulations, you think like a mathematician–yes, you too could be a poor man some day!  If you never bothered to wonder, then you lucked out, because that can ultimately be considered an irrelevant question in the first place.  But even the staunchest pedant can be convinced that this is a reasonable construction–which is a discussion I don’t really feel like having here.  Deal with it, nerd.

Ok, so weird thing in hand, let’s play with it.  Notice that

S = F_0 + F_1x + \displaystyle\sum_{n=2}^\infty F_n x^n\\= F_0 + F_1x + \displaystyle\sum_{n=2}^\infty \left(F_{n-1} + F_{n-2}\right) x^n\\=F_0 + F_1x + \displaystyle\sum_{n=2}^\infty F_{n-1} x^n + \displaystyle\sum_{n=2}^\infty F_{n-2} x^n\\ = F_0 + F_1x + x\displaystyle\sum_{n=2}^\infty F_{n-1} x^{n-1} + x^2\displaystyle\sum_{n=2}^\infty F_{n-2} x^{n-2}

and after re-kajiggering the indices (that’s the technical way to say it), and doing a little grouping, we have:

S = F_0 + F_1x + x\displaystyle\sum_{n=1}^\infty F_{n} x^{n} + x^2\displaystyle\sum_{n=0}^\infty F_{n} x^{n}\\ = F_0 + \left(F_1x + x\displaystyle\sum_{n=1}^\infty F_{n} x^{n}\right) + x^2\displaystyle\sum_{n=0}^\infty F_{n} x^{n}\\= F_0 + xS + x^2S\\=1+xS+x^2S

So that

S(1-x-x^2)=1

Hence

S = \displaystyle\frac{1}{1-x-x^2}

The last bit just devolves into a bunch of boring arithmetic, so I’ll keep the details fairly simple from here on.  Using the method of partial fractions decomposition, we can further write

S = \displaystyle\frac{1}{1-x-x^2}\\[.3cm]=\displaystyle\frac{\alpha_1}{1-\beta_1x}+\displaystyle\frac{\alpha_2}{1-\beta_2x}

and with only a little elbow grease, you can show that \beta_1=\frac{1+\sqrt{5}}{2}, \beta_2=\frac{1-\sqrt{5}}{2}, \alpha_1 = \frac{1}{\sqrt{5}}\beta_1, and \alpha_2 = -\frac{1}{\sqrt{5}}\beta_2.  Of course, these numbers are themselves quite well studied–typically \beta_1 is written as \phi, and given the interesting name “the golden ratio” (more on that in a moment).  So then

S=\displaystyle\sum_{n=0}^\infty \left( \alpha_1\beta_1^{n+1}+\alpha_2\beta_2^{n+1} \right)x^n

and by merely comparing coefficients among these two series, we have

F_n = \alpha_1\beta_1^n + \alpha_2\beta_2^n

which is exactly what we wanted to show!

Of course, the Golden Ratio has many cute properties and ties to the Fibonacci Numbers.  For example,

\displaystyle\lim_{n\rightarrow\infty}\frac{F_{n+1}}{F_n} = \phi

Actually, as it turns out, this limit shows up more commonly than you might expect.  But it does give some insight into the funny way some people define \phi–in terms of its so-called continued fraction expansion.  Starting with an example calculation, let’s look at F_5 and F_4.  Notice that

\displaystyle\frac{F_5}{F_4} = \displaystyle\frac{8}{5}\\=1+\displaystyle\frac{3}{5}\\=1+\displaystyle\frac{1}{\frac{5}{3}}\\=1+\displaystyle\frac{1}{1+\frac{2}{3}}\\=1+\displaystyle\frac{1}{1+\frac{1}{\frac{3}{2}}}\\=1+\displaystyle\frac{1}{1+\frac{1}{1+\frac{1}{2}}}\\=1+\displaystyle\frac{1}{1+\frac{1}{1+\frac{1}{1+\frac{1}{1}}}}

So you might be inclined to think that

\phi = 1+\displaystyle\frac{1}{1+\frac{1}{1+\frac{1}{1+\frac{1}{\ddots}}}}

And you would be right!

Now where were we?  Oh right; that R thing.  This closed form is easy enough to implement in R to give you the Fibonacci Number of your choice without having to compute all prior ones–but it’s not terribly efficient about it:

GetNthFib <- function(n){
  1/sqrt(5)*(((1+sqrt(5))/2)^{n+1} - ((1-sqrt(5))/2)^{n+1})
}

One final word before the R Project Euler solution: the Project Euler people use the indexing F_1=F_2=1, and then define the recurrence relation for all n>2. This is fine (and all R code beyond this point is a reflection of their choice), but it is worth noting that combinatorists–the people who have the most claim to that dumb list of numbers, generally agree that it should be indexed as F_0=F_1=1.

R Code:

ElapsedTime<-system.time({
##########################
### Function to get a vector of all Fibonacci #'s below n
FibsBelow <- function(n){
  fib <- c(1, 1)
  i <- 2
  while (sum(tail(fib, 2)) < n){
    i <- i+1
    fib[i] <- fib[i-1]+fib[i-2]
  }
  fib
}
###
n <- 4*10^6

fib.numbers<-FibsBelow(n)
answer<-sum(fib.numbers[-which(fib.numbers%%2 == 0)])
##########################
})[3]
ElapsedMins<-floor(ElapsedTime/60)
ElapsedSecs<-(ElapsedTime-ElapsedMins*60)
cat(sprintf("\nThe answer is:  %d\nTotal elapsed time:  %d minutes and %f seconds\n",answer, ElapsedMins, ElapsedSecs))

Output:
The answer is: 4613732
Total elapsed time: 0 minutes and 0.001000 seconds