Month: July 2008

Problem 5 on the Project Euler list is the following:

2520 is the smallest number that can be divided by each of the numbers from 1 to 10 without any remainder.

What is the smallest number that is evenly divisible by all of the numbers from 1 to 20?

This is asking us to compute the lowest common multiple (lcm) from 1 to 20. The lcm of two integers a and b is defined as

\[ \frac{ab}{\mbox{gcd}(a,b)} \]

In order to compute the gcd, we can use a simple algorithm, like a Euclidian reduction:

# Iterative Euclidean algorithm
gcd <- function(a,b) {
    t <- 0
    while (b != 0) {
        t <- b
        b <- a %% b
        a <- t
    }
    a
}

Defining the lcm then becomes:

lcm <- function(a,b) {
    abs(a * (b/(gcd(a,b))))
}

Now we have defined pretty much everything we need to solve the problem. We can write an iterative lcm routine that calculates the least common denominator from 1 to x:

lcm_ <- function(xs) {
    l <- 1
    for (i in 1:(length(xs))) {
        l <- lcm(l,xs[i])
    }
    l
}

And we can calculate the answer:

> lcm_(c(1:20))
[1] 232792560

However, I also experimented with writing a reduce() function (similar to Haskell’s foldr operation). We can write reduce() as follows:

reduce <- function(op,x,xs) {
    y <- x
    for (i in 1:(length(xs)-1)) {
        y <- do.call(op, list(y, xs[i]))
    } 
    y
}

And then invoke:

>reduce(lcm, 1, c(1:20))
[1] 232792560

Note that R actually has inbuilt Map/Reduce operations. Using the native Reduce, this is:

> Reduce(f=lcm, x=c(1:20), right=1)
[1] 232792560

Problem 4 in on Project Euler’s list is the following:

Find the largest palindrome made from the product of two 3-digit numbers.

We can attack this as follows: firstly, we know that for an n-digit number, the minimum n-digit number is \(10^{n-1}\) and the maximum is \(9\sum_{i=0}^{n-1}10^{i}\). So we need the cartesian product of all numbers in this range, which we can then test for the palindrome property, and finally filter for the maximum element.

Testing for the palindrome property requires that we check each digit in the number for equality on either side. We do this first by calculating the number of digits in the product (using \(\lfloor\log_{10}(x)\rfloor\) and then just comparing the digits by using the modulus and shifting. Note a couple of points in the R code below:

• %/% is integer division;
• Note the vectorization at work in statements like (p %/% 10^x): this will divide p by successive powers of 10, where the powers are defined in the x vector;
• Note the use of the logical all() function to test whether all palindromic contenders are equal.

pal < - function(N) {
 palindromes <- c()
 min_no <- 1*10^(N-1)
 max_no <- sum(9*10^(seq(0,N-1)))
 n <- 1
 for (i in seq(max_no,min_no,-1)) {
  for (j in seq(i,min_no,-1)) {
   p <- i * j
   digits <- floor(log10(p))
   x <- seq(0, digits %/% 2)
    if (all( ((p %/% 10^x) %% 10)==((p %/% 10^(digits-x)) %% 10) )) {
     palindromes[n] <- p
     n <- n + 1
    }
  }
 }
 max(palindromes)
}

The third problem in Project Euler is an interesting one:

The prime factors of 13195 are 5, 7, 13 and 29.

What is the largest prime factor of the number 600851475143 ?

In order to compute prime factors, we can use a sieving method. The simplest and oldest sieving algorithm is the sieve of Erastothenes, which we can implement in R as follows:

sieve <- function(n) {
    if (n <= 1) {
        return(numeric(0))
    }
    if (n == 2) {
        return(c(2))
    }

    a <- seq(2,n)
    b <- c()
    b[1] <- 2

    i <- 1
    while (length(a) > 0) {
        b[i] <- a[1]
        i <- i + 1
        a <- a[!(a %% a[1] == 0)]
    }
    b
}

So now we need to compute the largest prime factor of 600851475143. An upper bound for the factorisation is \(p \leq \lfloor\sqrt{x}\rfloor\), so we can halt our search at this point (alternatively, we may be able to begin our search here and work downwards). Let’s do it the obvious (and inefficient) way, starting at 1 and working up to \(\lfloor\sqrt{x}\rfloor\), building up a list of prime factors along the way:

pfactor <- function(x) {
    factors <- c()
    j <- 1
    primes <- sieve( floor(sqrt(x)) )
    for (prime in primes) {
        if (x %% prime == 0) {
            factors[j] <- prime
            j <- j + 1
    }
    }
    factors
}

Now to compute the maximum factor, we just take the maximum of this list:

max( pfactor(600851475143) )