The tenth problem in Project Euler revisits the prime numbers:

Calculate the sum of all the primes below two million.

Using the Sieve of Eratosthenes sequence developed earlier we can write the solution like this:

sieve 2L Map.empty |> 
    Seq.takeWhile (fun n -> n < 2000000L) |>
    Seq.sum |>
    printfn "sum = %A"

Arguably there is some inefficiency in generating a subsequence just consisting of the primes less than two million and then summing that sequence. I think the clarity of the resulting code makes up for the slight inefficiency.

Today I want to present my solution to problem nine in the Project Euler series:

Find the only Pythagorean triplet, {a, b, c}, for which a + b+ c = 1000.

A little bit of research on Wikipedia led to the page on Pythgorean triples. Using Euclid's formula we can generate a sequence of tuples which can be manipulated to find the desired answer:

let triples = 
    seq {
        for n in 1..100 do
            for m in n+1..100 do
                yield (2*m*n, m*m-n*n, m*m+n*n)
    }

Note that my selection of 100 as the endpoint for each of the loops is pretty arbitrary. I just needed a number big enough that the triple whose members summed to 1000 would be included in the sequence. Then to find the right one we can write:

triples |>
    Seq.find (fun (a, b, c) -> a+b+c = 1000) |>
    (fun (a, b, c) -> a*b*c) |>
    printfn "product = %A"

This sequence is lazily evaluated until we find the one we're looking for.

Yet another Project Euler post follows. This time I am looking at problem number eight:

Discover the largest product of five consecutive digits in the 1000-digit number.

The first thing to do was get the string of digits into a usable form. I copied them from the Project Euler web site, pasted them into a document and then did some manual labour to make them into a long string:

let chars = 
 "73167176531330624919225119674426574742355349194934"+
 "96983520312774506326239578318016984801869478851843"+
 "85861560789112949495459501737958331952853208805511"+
 "12540698747158523863050715693290963295227443043557"+
 "66896648950445244523161731856403098711121722383113"+
 "62229893423380308135336276614282806444486645238749"+
 "30358907296290491560440772390713810515859307960866"+
 "70172427121883998797908792274921901699720888093776"+
 "65727333001053367881220235421809751254540594752243"+
 "52584907711670556013604839586446706324415722155397"+
 "53697817977846174064955149290862569321978468622482"+
 "83972241375657056057490261407972968652414535100474"+
 "82166370484403199890008895243450658541227588666881"+
 "16427171479924442928230863465674813919123162824586"+
 "17866458359124566529476545682848912883142607690042"+
 "24219022671055626321111109370544217506941658960408"+
 "07198403850962455444362981230987879927244284909188"+
 "84580156166097919133875499200524063689912560717606"+
 "05886116467109405077541002256983155200055935729725"+
 "71636269561882670428252483600823257530420752963450"

I know this is kind of inefficient (each + operation will create a new immutable string object containing the contents of the previous operation with the new part appended) but it gets the job done. Now, in F#, a string is automatically a sequence so we can just process it like this:

chars |> 
    Seq.map (fun c -> int c - int '0') |>
    Seq.windowed 5 |>
    Seq.map (fun a -> Array.fold (*) 1 a) |>
    Seq.max |>
    printfn "max product = %A"

The idea here is that we convert the sequence of characters into a sequence of digits and then create a new sequence from it which consists of each set of five consecutive digits. Then we compute the product of each of these five digit sets and find the maximum.

This was one solution that I reworked extensively on revisiting it. My first attempt was not very idiomatic. I think this one is much better.

I will write more about music notation later as I explore LilyPond. In this post I want to explore my answer to problem #6 of Project Euler:

What is the difference between the sum of the squares and the square of the sums?

I first defined a helper function:

let square i = i*i

Then the sum of the squares can be expressed as:

let sum_of_squares = {1..100} |> Seq.map square |> Seq.sum

And the square of the sum is:

let square_of_sum = {1..100} |> Seq.sum |> square

So that the desired difference is:

square_of_sum - sum_of_squares |> printfn "diff = %A"

And there's really nothing much more to say.

Yesterday the Craftzine blog printed these words from Sarah Lohman's blog about historic cuisine:

Why bother deciphering a recipe over 150 years old?

You can take a collection of words and measurements written long ago, and turn it into a physical object. You can create something that looks, smells, and tastes just like it did hundreds of years ago. And that's the next best thing to time travel: it's a window to the past that lets you understand a little bit about another way of life.

I play recorders with a group of like-minded people. We mostly play music of the Renaissance (on Baroque recorders...). Do I feel a connection to the musicians who played this music five hundred years ago? To the composers who wrote the music? To the audience who listened to it? Or danced to it? Well, not really, but I like the idea.

Adieu mes amours by Josquin des Prez in the Odhecaton.

Project Euler Problem 5:

What is the smallest number divisible by each of the numbers 1 to 20?

There are several ways to solve this problem. I chose to compute the Least Common Multiple of two numbers using their Greatest Common Divisor:

Using Euclid's algorithm, the F# code for gcd is:

let rec gcd a b =
    match a with
    | a when a = 0L -> b
    | a when a < b -> gcd (b-a) a
    | a -> gcd (a-b) a

And lcm is just as easy:

let lcm a b = (a*b)/gcd a b

Then we can write the code that solves the problem:

seq {1L..20L} |> Seq.fold lcm 1L |> printfn "lcm = %A"

The fold operation calculates a running lcm. At each step the result of the previous step is passed to the lcm function with the current element in the sequence of integers.

The downside of posting my solutions to the Project Euler problems is that I end up second guessing myself; "is that really the best way to solve this?" Then I have to spend some time trying alternative ideas just to convince myself that my first idea really was good enough. Often it is. Sometimes I tweak my solution to make some part of it more reusable. Reduce, recycle, reuse:

On to Project Euler problem number four:

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

I created a number of helper functions. The first takes a number and returns a sequence containing the number's digits:

let rec digits n =
    seq {
        if n > 0 then
            yield n%10
            yield! digits (n/10)
    }

The next one takes a number a reverses it (that is, 1234 would become 4321):

let reverse n =
    digits n |> Seq.fold (fun acc k -> 10*acc + k) 0

With these helpers we can easily check whether a number is a palindrome:

let palindrome n = 
    n = reverse n

So now we can create a sequence containing all the palindromes that result from multiplying two three-digit numbers together:

let palindromes = 
    seq { 
        for i in 1..999 do 
            for j in 1..999 do 
                if palindrome (i*j) then yield i*j 
    }

Finally we can find the largest palindrome:

palindromes |> Seq.max |> printfn "max = %A"

I think that's a pretty good solution.

Growing up I had a set of Pentominoes and a rectangular tray in which they could be placed. I may even still have it somewhere. It's got to be more than 35 years old. Hmmm... I guess it would be pretty miraculous to find it tucked away in a box after four or five major moves (including one from Canada to the U.S.).

Recently a friend gave me a puzzle set with pieces that include some pentominoes and some other things. It seems to have a lot of solutions. Here's one of them:

Problem 3 of the Euler Project is:

Find the largest prime factor of a composite number.

As previously mentioned, having a lazy prime generator is kind of nice when it's time to solve this problem. The idea I settled on is to create a prime factorization of the target number. From that list I can then select the largest prime factor.

Here's my prime factorization code:

let rec factors n = 
    let root (n:int64) =
        int64 (sqrt (float n))
    let rec exp e p n =
        match n%p with
        | 0L -> exp (e+1L) p (n/p)
        | _  -> e
    let rec pow e p =
        match e with
        | 0L -> 1L
        | e  -> p*pow (e-1L) p  
    let primes = sieve 2L Map.empty
    seq {
        if n > 1L then
            let divisor = primes |> Seq.takeWhile (fun k -> k <= root n) |> Seq.tryFind (fun k -> n%k = 0L)
            match divisor with
            | None      -> yield (n, 1L)
            | Some(k)   ->
                let e = exp 0L k n
                yield (k, e)
                yield! factors (n/pow e k)
    }

I think this looks more complicated than it actually is because of the embedded helper functions: the root function computes the square root of an integer, exp computes the number of times a prime number appears in the prime factorization of a composite number and pow raises an integer to a power. It's not clear to me whether it would be better to have these as standalone functions or as embedded ones.

The result of passing a number into factors is a sequence of tuples where the first entry is the prime factor and the second is the number of times that factor goes into the composite number. This function will come in useful in my solutions to later Project Euler problems.

With this factorization function the solution to the problem is:

factors 600851475143L |> Seq.map (fun k -> fst k) |> Seq.max |> printfn "%A"

The mapping step basically throws away the exponent and allows us just to focus on the prime factor.

There is a some inefficiency in this code. Every time factors is called recursively it starts checking prime factors from the very first prime. Yet we know that we've already checked some of them previously. I briefly messed around with fixing this by always pulling from the same sequence of primes rather than creating a new one every time. Then I decided it was more trouble than it was worth.

Overexposed.

I am going to skip ahead to Project Euler problem number 7.
Find the 10001^st prime.

The solution to this problem involves creating a prime number generator which will come in very handy when solving some of the earlier problems.

I have been writing procedural code for something like 25 years now in assembler, FORTRAN, C, C++, Java and C#. I think that my first attempt at a solution shows this background.

I decided to use the Sieve of Eratosthenes as my way of finding primes.

The first problem I faced was "how big is the 10001^st prime?" Poking around in Wikipedia led me to a page on the prime counting function. There's a lot of interesting ideas there but for me the most useful one is that if p_n is the n^th prime then:

for large n.

With that tidbit of information solving this problem is pretty straightforward. The idea is that we create a bit array big enough that we are guaranteed that it will contain the 10001st prime, apply the Sieve algorithm to it, then scan for the one we are interested in:

let prime m =
    let count n =
        let f = float n
        int (f*(log f + log (log f)))
    let flags = new System.Collections.BitArray(count m, true)
    let filter n = 
        match flags.[n] with
        | true  -> for i in n*n..n..flags.Count do flags.[i] <- false
        | false -> ()
    let root n =
        let f = float n
        int (sqrt f)   
    seq {2..root flags.Count} |> Seq.iter filter
    seq {
        for i in 2..flags.Count-1 do if flags.[i] then yield i
    } |> Seq.nth (m-1)

prime 10001 |> printfn "%A"

Note that the filter function applies the "square optimization": we really only need to mark multiples of the current prime as composite starting with its square; smaller multiples of it will be removed when marking multiples of the smaller primes. This lets us avoid marking some numbers as composite multiple times.

This code solves the problem just fine and does so quite quickly (it took about a third of a second on my computer) but I just don't like it very much. It's functional but not very functional.

During my research on the Sieve I came across a paper by Melissa O'Neill called The Genuine Sieve of Eratosthenes. In it she talks about implementing the sieve algorithm in a functional language (Haskell). The idea, as presented in section 3, is to mark numbers as composite "just in time". Rather than use a fixed size bit array to mark which numbers are prime and which are composite it uses a map to keep track of the next multiple of a particular prime.

let rec sieve n table =
    let rec update n p table =
        let n = n+p
        match Map.tryFind n table with
            | None          -> Map.add n p table
            | Some(prime)   -> update n p table
    seq {
        match Map.tryFind n table with
        | None              -> 
            yield n
            yield! sieve (n+1L) (table |> Map.add (n*n) n)
        | Some(prime)       ->
            yield! sieve (n+1L) (Map.remove n table |> update n prime)
    } 

We seed this routine with the first prime and an empty map. Each time we request the next element in the sequence the function checks to see if the current number is in the table (and hence is composite). If not, the number is returned and then the map is updated with the first composite number that is not a multiple of a previously discovered prime (i.e. the square of the current number). If the current number is in the table then it is composite so we update the table and then check the next number. Once a number has been discovered to be composite we don't need it in the table any more (the sequence delivers up primes one at a time with no backtracking allowed) so we remove it and add the next multiple of the underlying prime.

This algorithm has a couple of nice features:

1. The map only contains one entry for each prime we discover.
2. With a minor change to use BigInts instead of long ints we could generate primes until your computer runs out of memory.

At this point a small helper function will get us the 10001^stprime:

let primeprime n =
    sieve 2L Map.empty |> Seq.nth (n-1) 

primeprime 10001 |> printfn "%A"

Update: This routine runs a bit slower than the first one. It took about half a second to compute the answer. I believe this is because of the immutable map structure used. Every time an item is added or removed a revised copy of the map is created. After 10,000 primes it's doing a lot of work just in bookkeeping. I imagine you could improve performance by using a mutable data structure like System.Collections.Generic.Dictionary. I haven't cared enough about performance to make this change. 

The rain just would not stop yesterday. Our rainwater tanks are full. They hold about 6000 gallons. One storm equals full tanks!

Update: Apparently we got 11 inches of rain in one day.

I've been meaning to mention that I'm using MathCast to format the equations in this blog. I render them as images which I then embed in the post.

Continuing on with the Project Euler solutions I'm now going to look at the second problem:


Find the sum of all the even-valued terms in the Fibonacci sequence which do not exceed four million.


Back when I studied math (B.Sc. Math, UBC 1980) the definition of the Fibonacci sequence was usually given as something like:

The sequence defined by this recurrence relation is:

1, 1, 2, 3, 5, 8, 13, 21, ...

For the purposes of this problem the series is defined a bit differently:

So that the sequence is:

1, 2, 3, 5, 8, 13, 21, 34, ...

If I used the first definition my answer would be off by one. There are a lot of ways to implement this in F#. After briefly flirting with a matrix multiplication formula (in a foolish attempt to avoid recursion) I came up with this:

let fib m n =
    let rec inner a b =
        seq {
            let c = a + b
            yield c
            yield! inner b c
        }
    seq {
        yield m
        yield n
        yield! inner m n
    }

For the purposes of this sequence the two parameters of the function are the first two values in the series. I coded it this way so that type inference would allow it to be used with any integer type. This is definitely not the simplest F# implementation of the Fibonacci series. A naive implementation might look like this:

let rec fib_naive n =
    match n with
    | 1|2 -> n
    | k -> fib_naive (k-1) + fib_naive (k-2)

The problem with this code is that in the computation of the n-1^th element requires the computation of the n-2^nd element. But then you turn around and recompute the  n-2^nd element all over again. In other words you do the same computation twice which, for large values of n, is just wasteful. Also, to use this function you would need to know the approximate value of n for which fib_naive n is greater than 4 million.

The part I like about the first function is that at each step you already have the previous value of the series computed and you can just use it. Plus the first function is really a generator. You can get the next element from it as many times as you want (or as is allowed by the underlying representation of an integer).

So, given the Fibonacci number generator, the solution to this problem can be written as:

fib 1 2 |> 
    Seq.takeWhile (fun n -> n < 4000000) |> 
    Seq.filter (fun n -> n%2 = 0) |> 
    Seq.sum |> 
    printfn "sum = %A"

To put it in words, take the Fibonacci series and select the values from it that are less than 4 million. Then select all the even numbers in the new subset and add them together.

I find myself somewhat obsessed with Project Euler these days. It contains a series of math challenges that (mostly) require a computer program of some sort to find the right answer. If you start looking around the web there seems to be a tradition of using these problems as a way to learn a new programming language. And that's what I'm doing too. My language of choice is F#.

The first problem in the project is:

Add all the natural numbers below one thousand that are multiples of 3 or 5.

I have come up with four solutions. Here's number one:

seq {for i in 0..999 do if i%3=0 or i%5=0 then yield i} |> Seq.sum |> printfn "total = %A"

Solution number two is to sum the sequence of integer multiples of 3 and 5 and then subtract the overlap (the sum of the multiples of 15). That would look like this:

({3..3..999} |> Seq.sum) + ({5..5..999} |> Seq.sum) - ({15..15..999} |> Seq.sum) |> printfn "total = %A"

The third solution is another variation on the theme of summing a sequence:

seq {1..999} |> Seq.filter (fun k -> k%3=0 or k%5=0) |> Seq.sum |> printfn "total = %A"

Finally, solution number four comes from noting that the sum from k to n of the multiples of k can be written like this:

With some tweaking to get the end of the series correct I translated the math into F# as follows:

let sigma k n = 
    let m = k*((n-1)/k)
    (m*(m+k))/(2*k)

In this case m is the largest integer less than n that needs to be included in the sum (this is not the same as n since the problem states that we are looking at numbers below one thousand). With this helper function the solution to the problem becomes:

(sigma 3 1000) + (sigma 5 1000) - (sigma 15 1000) |> printfn "total = %A"

To time these solutions I used the #time functionality of F# Interactive.

Solution 1 took 00:00:00.010.
Solution 2 took 00:00:00.011.
Solution 3 took 00:00:00.024.
Solution 4 took 00:00:00.003.

It's not at all clear to me why the third solution is two times slower than the first two. As expected the fourth solution is by far the fastest.

And that's way too much analysis for such a simple problem.

In case you were wondering, the outfit Mallory is wearing in the previous picture was not included in the package when we took possession. She came to us with a dress and pair of squeaky shoes. The squeaky shoes were cute for about a minute. Then we disabled them.