Large integers in R: Fibonacci number with 1000 digits, Euler Problem 25

The Fibonacci Sequence in nature: The nautilus shell. Euler problem 25

The Fibonacci Sequence occurs in nature: The nautilus shell.

Euler Problem 25 takes us back to the Fibonacci sequence and the problems related to working with very large integers.

The Fibonacci sequence follows a simple mathematical rule but it can create things of great beauty. This pattern occurs quite often in nature, like to nautilus shell shown in the image. The video by Arthur Benjamin at the end of this post illustrates some of the magic of this sequence.

Large Integers in R

By default, numbers with more than 7 digits are shown in scientific notation in R, which reduces the accuracy of the calculation. You can change the precision of large integers with the options function but R struggles with integers with more than 22 digits. This example illustrates this issue.

factorial(24)
[1] 6.204484e+23
> options(digits=22)
> factorial(24)
[1] 620448401733239409999872

However, finding a number of 1000 digits is a problem with using special functions.

Euler Problem 25 Definition

The Fibonacci sequence is defined by the recurrence relation:

F_n = F_{n-1} + F_{n-2} , where F_1 = 1 and F_2 = 1 . The 12th term, F_{12} , is the first term to contain three digits.

What is the index of the first term in the Fibonacci sequence to contain 1000 digits?

Proposed Solutions

The first solution uses the GMP library to manage very large integers. This library also contains a function to generate Fibonacci numbers. This solution cycles through the Fibonacci sequence until it finds a number with 1000 digits.

library(gmp) # GNU Multiple Precision Arithmetic Library
n <- 1
fib <- 1
while (nchar(as.character(fib)) < 1000) {
   fib <- fibnum(n) # Determine next Fibonacci number
   n <- n + 1
}
answer <- n
print(answer)

This is a very fast solution but my aim is to solve the first 100 Project Euler problems using only base-R code. The big.add function I developed to solve Euler Problem 13.

t <- proc.time()
fib <- 1 # First Fibonaci number
cur <- 1 # Current number in sequence
pre <- 1 # Previous number in sequence
index <- 2
while (nchar(fib) < 1000) {
    fib <- big.add(cur, pre) # Determine next Fibonacci number
    pre <- cur
    cur <- fib
    index <- index + 1
}
answer <- index
print(answer)

This code is much slower than the GMP library but is was fun to develop.

The Magic of the Fibonacci Numbers

Euler Problem 20: Large Integer Factorials

Euler Problem 20Euler Problem 20 is the third problem that requires special consideration for working with very large integers. In this problem, we look at factorials. These numbers are useful in combinatorics if, for example, you like to know in how many ways you can arrange a deck of cards.

The problem with computing factorials is that they are mostly very large numbers, beyond the generic capabilities of computers to process. This problem can be solved using a specialised R package and using only base-R code.

Euler Problem 20 Definition

n! = n \times (n - 1) \times (n-2) \times \ldots \times 3 \times 2 \times 1 .

For example:  10! = 10 \times 9 \times \ldots \times 3 \times 2 \times 1 = 3628800 .

The sum of the digits in the number 10! is 3 + 6 + 2 + 8 + 8 + 0 + 0 = 27 .

Find the sum of the digits in the number 100!

Euler Problem 20 Solution

The factorial of the number 100 contains 158 digits, which is a lot more digits than a 64-bit operating system can accurately produce. Using the standard function: factorial(100) = 9.332622e+157. Without using a specialised algorithm, we cannot determine the sum of all digits. We need to deploy arbitrary-precision arithmetic to solve this problem.

Many computer languages, including R, have special libraries to deal with such large numbers. The GMP Multiple Precision Arithmetic package renders this problem almost trivial.

library(gmp)
digits <- factorialZ(100)
digits <- as.character(digits)
answer <- sum(as.numeric(unlist(strsplit(digits, ""))))

Base-R Code

The problem becomes more interesting when only using basic R code. I developed the big.add function to solve Euler Problem 13 through the addition of very large integers. We can extend this function to also calculate factorials. A factorial can be replaced by a series of additions, for example:

3! = 1 \times 2 \times 3 = (((1+1) + (1+1)) + (1+1))

This can be mimicked in R using the Reduce function. This function reduces a vector to a single value by recursively calling a function. Reduce(“+”, rep(4, 5)) is the same as:

4 \times 5 = ((((4 + 4) + 4) + 4) + 4) = 20

Using a loop, we can use the Reduce function to calculate a factorial, using only additions:

fact <- 1
x <- 100
for (i in 2:x) {
    fact <- Reduce("+", rep(fact, i))
}
print(fact)

The big.factorial function below implements this idea by combining the big.add and Reduce functions to calculate large integer factorials. The function returns a value of 1 when the factorial of 0 or 1 is requested. This function does not calculate the Gamma-function for fractions. For all other values, it goes through a loop from 2 to the requested factorial. The temporary values are stored in the bf variable. The code loops through the factorials by using the result of the previous Reduce call into the current one.

big.factorial <- function(x) {
    x <- floor(x)
    bf <- 1 if (x > 1) {
        for (i in 2:x) {
            bf <- Reduce(big.add, rep(bf,i))
        }
    }
    return (bf)
}

digits <- big.factorial(100)
answer <- sum(as.numeric(unlist(strsplit(as.character(digits), ""))))
print(answer)

This function is most certainly not as fast as the GMP package but it was fun to write and to learn about the mechanics behind arbitrary precision arithmetic at work.

If you like to know how factorials can be used to determine the number of ways a deck can be shuffled the watch this video.

Euler Problem 16: Power Digit Sum

Euler Problem 16: Power Digit SumEuler Problem 16 is reminiscent of the famous fable of wheat and chess. Lahur Sessa invented the game of chess for King Iadava. The king was very pleased with the game and asked Lahur to name his reward.

Lahur asked the king to place one grain of rice on the first square of a chessboard, two on the next square, four on the third square and so on until the board is filled. The king was happy with his humble request until his mathematicians worked out that it would take millions of tonnes of grain. Assuming there are 25 grains of wheat in a gramme, the last field will contain more than 461,168,602,000 tonnes of grain.

Euler Problem 16 Definition

2^{15} = 32768 and the sum of its digits is 3 + 2 + 7 + 6 + 8 = 26 . What is the sum of the digits of the number 2^{1000} ?

Solution

The most straightforward solution uses the GMP package for Multiple Precision Arithmetic to calculate big integers. The as.bigz function results in a special class of arbitrarily large integer numbers

# Raise 2 to the power 1000
library(gmp)
digits <- as.bigz(2^1000) # Define number
# Sum all digits
answer <- sum(as.numeric(unlist(strsplit(as.character(digits), ""))))
print(answer)

We can also solve this problem in base-r with the bigg.add function which I developed for Euler Problem 13. This function uses basic string operations to add to arbitrarily large numbers. Raising a number to the power of two can also be written as a series of additions:

2^4 = 2 \times 2 \times 2 \times 2 = ((2+2)+(2+2)) + ((2+2)+(2+2))

The solution to this problem is to add 2 + 2 then add the outcome of that equation to itself, and so on. Repeat this one thousand times to raise the number two to the power of one thousand.

# Raise 2 to the power 1000
pow <- 2
for (i in 2:1000)
    pow <- big.add(pow, pow)
# Sum all digits
answer <- sum(as.numeric(unlist(strsplit(pow, ""))))
print(answer)

Euler Problem 13: Large Sum of 1000 Digits

Euler Problem 13 asks to add one hundred numbers with fifty digits. This seems like a simple problem where it not that most computers are not designed to deal with numbers with a lot of integers. For example:

 2^{64} = 18446744073709551616

When asking R to compute this value we get 1.844674e+19, losing most of the digits and limiting the accuracy of the results. Computers solve this problem using Arbitrary-precision Arithmetic. There are many software libraries that can process long integers without loosing accuracy. Euler Problem 13 requires this type of approach.

Euler Problem 13 Definition

Work out the first ten digits of the sum of the following one-hundred 50-digit numbers.

Solution

The easy way to solve this problem is to use the gmp package for working with very large integers. This package uses a special number types such as Big Rational and Big Integer. The number of digits in these number types is only limited by the size of the memory.

library(gmp)
numbers <- readLines("Euler/p013_numbers.txt")
digits <- sum(as.bigz(numbers))
answer <- substr(as.character(digits),1,10)

Using Base-R

To find the solution to this problem using only base R, I wrote a function to add numbers using strings instead of integers. The function adds leading zeros to the smallest number to make them both the same length. The function then proceeds to add numbers in the same way we were taught in primary school. This function can in principle be used for several other Euler Problems using large integers.

# Add numbers with many digits
big.add <- function(a, b) {
    # Add leading zeros to smallest numer
    if (nchar(a) < nchar(b)) 
        a <- paste0(paste(rep(0, nchar(b) - nchar(a)), collapse = ""), a) 
    if (nchar(a) > nchar(b)) 
        b <- paste0(paste(rep(0, nchar(a) - nchar(b)), collapse = ""), b)
    solution <- vector()
    remainder <- 0
    for (i in nchar(b):1) {
        p <- as.numeric(substr(a, i, i))
        q <- as.numeric(substr(b, i, i))
        r <- p + q + remainder 
        if (r >= 10 & i!=1) {
            solution <- c(solution, r %% 10)
            remainder <- (r - (r %% 10))/10
        } else {
            solution <- c(solution, r)
            remainder <- 0
        }
    }
    return(paste(rev(solution), collapse = ""))
}

With this function, the problem is easy to solve. The second part of the code runs this function over the one hundred numbers provided on the Euler Problem page and calculates the answer.

numbers <- readLines("Euler/p013_numbers.txt")
for (i in numbers) {
    answer <- big.add(answer, i)
}
answer <- substr(answer, 1, 10)

Multiplying Big Numbers

You can expand this function to multiply a very large number with a smaller number using the Reduce function. This function adds the number a to itself, using the big.add function. The outcome of the addition is used in the next iteration until it has been repeated b times. The number b in this function needs to be a ‘low’ number because it uses a vector of the length b.

big.mult <- function(a, b) {
    Reduce(big.add, rep(a, as.numeric(b)))
}