Project Euler Problem #15

Problem 15 on Project Euler asks us to find the number of distinct routes between the top left and bottom right corners in a 20×20 grid, with no backtracking allowed.

I originally saw this type of problem tackled in the book Notes On Introductory Combinatorics, by George Polya amongst others. This book is hard to find now, but it is a really clear intro to combinatoric math.

The solution can be paraphrased as follows: if the grid is of size 20×20, and it takes 2 movements to navigate a single square in the grid, then we must make a total of 40 movements to get from the top right to the bottom left. Exactly half of these movements will be left-to-right, and the other half will be up-down. The total number of distinct routes is the number of ways that we can choose 20 of each type of move from the 40 total moves required. So we need the combinatoric construct n-choose-k, or how many ways k items can be selected from n total items. This is represented as \(n\choose k\).

In R, calculating \[{40\choose 20}\] is just:

choose(40, 20)

Learning Lilypond : Luciana’s Theme

There is a whole bunch of obscure and wonderful music out there that I have been meaning to transcribe and arrange for the piano, and so I have been looking around for decent music transcription software, and I recently came across lilypond. This is to music engraving what TeX is to mathematical typesetting, and the results are really very nice. The downside of lilypond (as with TeX) is the slightly steep initial learning curve. To put it to the test, I have transcribed a piece of music that I came across a while ago – it’s called “Tema Par Luciana”, and is from an old Italian movie called C’eravamo tanto amati. I haven’t seen the movie, but I came across the theme in a mix created by a French DJ called elektrosonik. The piece is simple enough that it was pretty straightforward to transcribe, but also has enough variety that it is a good starting point for learning lilypond techniques.

Here is the PDF output (Lilypond will also produce a MIDI file if asked):

Tema Par Luciana (PDF)

And here is the lilypond source:

Project Euler Problem #13

Problem 13 on Project Euler asks us to sum 100 50-digit numbers and give the first 10 digits of the result. This is pretty easy. Note we are using R’s integer division operator %/% to discard the remainder of the large summed integer and just gives us the first 10 digits of the result.

## Problem 13
problem13 <- function() {
    nums <- scan("problem13.dat")
    s <- sum(nums)
    s %/% 10^(floor(log10(s))9)

Project Euler Problem #14

Problem 14 on the Project Euler site asks us to find the longest chain under 1 million created using the Collatz mapping. This is fairly straightforward, although performance again is not great:

## Problem 14
# Collatz conjecture
problem14 <- function(N) {
        maxChain <- 0
        chains <- rep(0,N)
        x <- 1
        for (i in 1:N) {
                n <- i
                chain <- 0
                while(n > 1) {
                        n <- ifelse(n %% 2 == 0, n/2, 3*n+1)
                        chain <- chain + 1
                        if (n < N && chains[n] > 0) {
                            chain <- chain + chains[n]
                chains[i] <- chain
                if (chain > maxChain) {
                        maxChain <- chain
                        x <- i

Project Euler Problem #12

Problem 12 on the Project Euler site asks:

What is the value of the first triangle number to have over five hundred divisors?

A triangular number T(n) is defined as \(T(n) = \frac{n(n+1)}{2}\). The R code below consists of a solution, which involves the fact that the number of proper divisors of an integer n can be calculated by first computing a prime factorisation of the number n, e.g. if <\(n = p^aq^b\), where p,q are prime, then the number of proper divisors of n can be calculated as \(d(n) = (a+1)(b+1)\). This solution is extremely slow (mainly due to the naive prime sieving algorithm), and could be speeded up dramatically with a little effort. # Sieve of Eratosthenes
prime.sieve <- function(n) {
 a <-,n)
 p <- 1
 M <- as.integer(sqrt(n))
 while ((p <- p + 1) <= M) {
  if (a[p] != 0)
    a[*p, n, p)] <- 0
 a[a>1 & a>0]

# Trial Division
# Returns the exponents of the prime
# factors of n
# e.g. if n = p^a*q^b
# tdiv(n) will return (a,b)
tdiv <- function(n) {
 primes <- prime.sieve(n)
 factors <- c()
 i <- 1
 curr <- 0
 for (p in primes) {
  while (n %% p == 0) {
   curr <- curr + 1
   n <- n %/% p
  factors[i] <- curr
  i <- i + 1
  curr <- 0
 factors[factors > 0]

# Compute nth triangular number
T <- function(n) {
## Problem 12
# This is a slooow solution
problem12 <- function(N) {
 n <- 0 # current triangular number Tn
 i <- 5 # \sum_{i=1}^n{i}
 while (TRUE) { 
  n <- T(i)
  factors <- tdiv(n) 
  if (prod(factors+1) >= N) {
  i <- i + 1

Project Euler Problem #19

Problem 19 on the Project Euler website asks the user, given some initial information:

How many Sundays fell on the first of the month during the twentieth century (1 Jan 1901 to 31 Dec 2000)?

The obvious (but longer) way is to calculate the sum of the days between 1901 and 2000, given the number of days in each month, and a helper function to determine whether a year is a leap year or not:

is.leap <- function(year) {
        return (year %% 4 == 0 || (year %% 100 == 0 && year && 4 == 0))

# Problem 19
problem19 <- function() {
        daycount <- 1
        daylist <- list()
        i <- 1
        for (year in 1900:2000) {
                months <- c(0,31,28,31,30,31,30,31,31,30,31,30,31)
                if (is.leap(year)) {
                        months[3] <- 29
                        days <- daycount + cumsum(months)
                        daycount <- days[length(days)]
                        daylist[[i]] <- (days[(length(days))])
                        i <- i + 1
                sum(unlist(lapply(daylist[1], function(x) {sum(x %% 7==1)} )))

However, with the aid of R’s chron library, there is a much easier way:

# Problem 19, method 2
sum(weekdays(seq.dates(”01/01/1901″, â€œ12/31/2000″, by=”months”))==”Sun”)