Project Euler in R: Problem 25

Solutions in R to Project Euler problems are generally hard to find, so I've recently started posting R solutions on this site.  Here are the previous problems: problem 22, problem 23 and problem 24

Now let's look at Problem 25 from Project Euler.  Here is the problem statement:

The Fibonacci sequence is defined by the recurrence relation: Fn = Fn1 + Fn2, where F1 = 1 and F2 = 1.
Hence the first 12 terms will be:
F1 = 1
F2 = 1
...
F11 = 89
F12 = 144
The 12th term, F12, is the first term to contain three digits. What is the first term in the Fibonacci sequence to contain 1000 digits?

This is an easy problem to solve for a small number of digits. It gets interesting once we increase the required number of digits to 1000. The R integer type can hold 10 digits at most, and the type double can hold 309 digits.  There is no Big Integer class in R that can hold larger numbers, so to compute the 1000 digit Fibonacci number we will have to be creative.

Computing Fibonacci numbers iteratively

Here is a basic iterative Fibonacci function, that works for Fibonacci numbers with 309 or fewer digits.

fib <- function(n) {
  a = 0
  b = 1
  for (i in 1:n) {
    tmp = b
    b = a
    a = a + tmp
  }
  return (a)
}

The version that works with larger numbers is similar, with the following changes

• Integer vectors are used to store each Fibonacci number.
• Each element of the vector is at most nine digits long.  
• Once all the nine digits are used, a new element is added to the vector and the Fibonacci number computation is carried over to the new element.
• The variable numberDigits counts the number of digits in the latest Fibonacci number.  When numberDigits exceeds 999, the function exits, returning the Fibonacci number index.

fib.bigint <- function(lim) {
  a = c(0)
  b = c(1)
  n = 0           # Fibonacci number index
  integerSize <- 1000000000  ## How big is each integer
  while (1) {
    n = n + 1     # Compute next Fibonacci number
    tmp = c(rep(0,length(a)-length(b)), b)  # Pad beginning with zero
    b = a
 
    # Add a and tmp
    carry = 0
    for (i in length(a):1) {
      sum.i = tmp[i] + a[i] + carry
      a[i] = sum.i %% integerSize
      carry = floor (sum.i / integerSize)
    }
    if (carry > 0) {a <- c(carry, a)}
    numberDigits <- (length(a) - 1) * 9 + countDigits(a[1])
    if (numberDigits > (lim - 1)) { return (n) }
  }
}

countDigits <- function(n) {
  count <- 1
  for (i in 1:9){
    if ((n/10) > 1) {
      count <- count + 1
      n = n / 10
    } else {
      return (count)
    }
  }
}

On my Intel Core 2 Duo 1.6Ghz laptop running Linux, this function less than 4 seconds to find the Fibonacci number with 1000 digits.

The analytical solution
Fibonacci numbers have a closed form solution, which leads to a simpler approximation by rounding (from Wikipedia),

F_n = the closest integer to (golden ratio)^n / sqrt(5) 

We need the first Fibonacci number with 1000 digits.  This number will be greater than or equal to 10^999.

F_n = the closest integer to (golden ratio)^n / sqrt(5) >= 10^999

We can drop the closest integer function because if F_n rounds off to the floor of the ratio, the resulting Fibonacci number will have less than the required digits.  So we need smallest n satisfying

(golden ratio)^n / sqrt(5) >= 10^999 

Or, n >= (999 log(10) + log(sqrt(5)))/log(golden ratio) 

Or, n >= 4781.86

So we need the 4782nd Fibonacci number to get the required 1000 digits.

This problem was particularly interesting to R users because it runs into an R-specific limitation. Many other languages don't have this problem and can just use a simple Fibonacci generator to get the answer. For example, in Scheme (or other Lisp variants), a naive Fibonacci implementation produces all the digits.

Aside: Having started writing about R, I have found an excellent online resource, R-bloggers -  an aggregator of content from various blogs about R. You can search for other R Project Euler solutions here.  Be sure to check out the top articles box on R-bloggers landing page if you're looking for new and interesting applications of R.