The next Project Euler problem I want to tackle is number 22:

What is the total of all the name scores in the file of first names?

The file in question contains data that looks like this:

"MARY","PATRICIA","LINDA","BARBARA",...

It's a nice list but I will need to clean it up to get it into a form that is easy to process. To that end I define a couple of helper functions. These allow me to use .NET's Regex class in a more functional manner:

open System.Text.RegularExpressions

let split pattern input =

    (new Regex(pattern)).Split(input)

let replace pattern (replacement:string) input =

    (new Regex(pattern)).Replace(input, replacement)

The value for a name is computed by adding the values for each letter ( a = 1, b = 2, c = 3, etc.). This is pretty easy in F# because a string can be treated as a sequence of characters. Thus the value of a name can be computed like this:

let value (name:string) =

    name |> Seq.map (fun c -> (int)c - (int)'A' + 1) |> Seq.sum

All the names in the file are already in uppercase so we don't have to do any normalisation there.

Once all the names have a value I need to compute a score for each name. This score is the name's value multiplied by its position in the alphabetically sorted list.

So, in words, here's what we do:

1. Read in the names.
2. Sanitize them by getting rid of the extraneous quotes.
3. Split the resulting string into an array of name strings.
4. Sort the array.
5. Compute the score for each name.
6. Sum the resulting scores.
7. Display the result.

Using the previously defined functions this algorithm translates directly into F# code:

open System.IO

File.ReadAllText(@"e:\Project Euler\names.txt")

    |> replace "\"" ""

    |> split ","

    |> Array.sort

    |> Array.mapi (fun i name -> (i + 1)*(value name))

    |> Array.sum

    |> printfn "total = %A"

The twenty-first Project Euler problem is:

Evaluate the sum of all amicable pairs under 10000.

An amicable pair is two numbers such that the sum of the proper divisors of the first number adds up to the second and vice versa. Each member of the pair is considered to be an amicable number. As always Wikipedia has a good article on the subject.

Following the suggestion from the statement of the problem I define a function d(n) that computes the sum of n's divisors:

let d n =

    {1..n/2} |> Seq.filter (fun i -> n%i = 0) |> Seq.sum

Now I need to test whether a given number is amicable:

let amicable n =

    let sum = d n

    (sum <> n) && (d sum = n)

Note that we exclude so called perfect numbers where the sum of n's divisors is n.

Finally, I can compute the sum of all the amicable numbers under 10,000 as follows:

{1..10000} |> Seq.filter (fun i -> amicable i) |> Seq.sum |> printfn "total = %A"

Numerically, the next Project Euler problem is number 20:

Find the sum of digits in 100!

Hmmm... what to do, what to do? Well, as always, lets start with the helper functions. 100! is a very big number. Thus we will have to use F#'s BigInt type. Here's a factorial function:

let factorial n = Seq.fold ( * ) 1I [1I .. n]

Of course, F# already has a factorial function for BigInts but since this is a learning exercise I felt that I had to write my own.

Next I want to turn a big number into a series of digits. This is the function I used for that purpose in my solution to problem number 16:

let rec digits n =

    seq {

        if n > 0I then

            yield n%10I

            yield! digits (n/10I)

    }

Finally we can put these together to create a solution:

factorial 100I |> digits |> Seq.sum |> printfn "sum = %A"

On my computer I get the answer in less time than it takes to threaten the monitor with my fist.

The next Project Euler problem on my list is number 19:

How many Sundays fell on the first of the month during the twentieth century?

For this problem it seemed to me that the easiest attack would be to create a sequence containing all the first days of each month in the given range and then filter it by testing whether the day in question was a Sunday. To make this work I would need a way of calculating the day of the week from the date. Fortunately, once again Wikipedia came to the rescue.

Following the algorithm from this article we define a few helper functions. To calculate the century we say:

let century y = 2 * (3 - (y/100)%4)

To determine whether we have a leap year:

let leap y =

    if y%400 = 0

        then true

        else

            if y%100 = 0

                then false

                else

                    if y%4 = 0

                        then true

                        else false

For the month table I used a pattern match expression which takes the month number and whether it's a leap year or not:

let month l m =

    match m with

    |  0 -> if l then 6 else 0

    |  1 -> if l then 2 else 3

    |  2 -> 3

    |  3 -> 6

    |  4 -> 1

    |  5 -> 4

    |  6 -> 6

    |  7 -> 2

    |  8 -> 5

    |  9 -> 0

    | 10 -> 3

    | 11 -> 5

Note that this match generates a warning to the effect that there are cases which may not be covered:

      match m with
  ----------^
stdin(19,11): warning FS0025: Incomplete pattern matches on this expression. For example, the value '12' may indicate a case not covered by the pattern(s).

In this case I am controlling all the inputs and I know that I would never make a mistake so I decided not to worry about it.

Next I needed to compute the year value:

let year y = y%100 + (y%100)/4

With these helpers, given a date tuple in the form (y, m, d) the day of the week it falls on is given by:

let day (y, m, d) = (century y + year y + month (leap y) m + d)%7

In this scheme Sunday is day 0.

The list of dates I want to examine are given by this sequence:

let dates =

    seq {

        for y in 1901..2000 do

            for m in 0..11 do

                yield (y, m, 1) }

Putting the pieces together to solve the puzzle:

dates |> Seq.filter (fun x -> day x = 0) |> Seq.length |> printfn "count = %A"

Returning to my irregularly scheduled series on Project Euler solutions I'd like to look at problem 18:

Find the maximum sum travelling from the top of the triangle to the base.

Here is a little example:

1

3 2

4 1 6

I have marked the maximum sum with red numbers. For a small triangle of numbers like this you could easily brute force the solution and just check all paths from the top to the bottom. For a large triangle, that would be too labour intensive. A better solution, I think, is to coalesce the rows of the triangle working from the bottom up. So, in this simple example, you would first replace the second row with the maximum value possible by adding a number in that row to the two numbers below it in the triangle. In this case, 3 + 4 > 3 + 1 so you would replace 3 by 7 and 2 + 6 > 2 + 1 so you'd replace 2 by 8. That would reduce the triangle to this:

1

7 8

Performing the same step on the reduced triangle would give you the final answer (since 1 + 8 > 1 + 7).

9

To solve this problem I put the data into an array of arrays:

let data =

    [|

    [|04; 62; 98; 27; 23; 09; 70; 98; 73; 93; 38; 53; 60; 04; 23|];

    [|63; 66; 04; 68; 89; 53; 67; 30; 73; 16; 69; 87; 40; 31|];

    [|91; 71; 52; 38; 17; 14; 91; 43; 58; 50; 27; 29; 48|];

    [|70; 11; 33; 28; 77; 73; 17; 78; 39; 68; 17; 57|];

    [|53; 71; 44; 65; 25; 43; 91; 52; 97; 51; 14|];

    [|41; 48; 72; 33; 47; 32; 37; 16; 94; 29|];

    [|41; 41; 26; 56; 83; 40; 80; 70; 33|];

    [|99; 65; 04; 28; 06; 16; 70; 92|];

    [|88; 02; 77; 73; 07; 63; 67|];

    [|19; 01; 23; 75; 03; 34|];

    [|20; 04; 82; 47; 65|];

    [|18; 35; 87; 10|];

    [|17; 47; 82|];

    [|95; 64|];

    [|75|]

    |]

You will notice that I inverted the triangle so that the widest part is at the bottom. I wanted my algorithm to proceed from the start of the array to end and not vice versa.

Next I defined a helper function:

let max x y = if x < y then y else x

Now I want a function that will coalesce two arrays into one:

let coalesce (a:int[]) (b:int[]) =

    [| for i in 0..b.Length-1 -> max (a.[i]+b.[i]) (a.[i+1]+b.[i]) |]

So, the idea here is that two one-dimensional arrays are passed in. The first one is assumed to be one longer than the second (in production code you'd probably want to enforce the precondition). This function creates a new one-dimensional array with the same length as the second input array. Each entry in the new array is the max of the entry in the second array with either of the numbers above it in the first array.

Finally, I want to apply the coalesce operation to the entire data array:

data |> Array.fold coalesce (Array.zeroCreate (data.[0].Length+1)) |> printfn "%A"

I have created an array with length one greater than the maximum line in the data array and I'm using that as the state accumulator passed to the fold function. At each step in the fold operation it will be replaced by an array one entry smaller until, when it's done, there will just be an array with a single element which will be the final answer.

I should note that one potential problem with this solution arises if you want to know the actual path that led to the answer. 

Project Euler's seventeenth problem wants to know:

How many letters would be needed to write all the numbers in words from 1 to 1000?

For example, in the British style, one would write out 256 as "two hundred and fifty-six". We are told not to count spaces and hyphens so 21 letters are needed for this number.

My approach is to decompose a number into a tuple containing the number of thousands, hundreds, tens and ones in the number. However, to make things simpler for me I chose to special case all the numbers under 20 since their spelling is a irregular. That is, instead of saying something like "ten plus one" we say "eleven". It is easier to break that out at this stage in the problem.

Here is my F# code to decompose a number:

let decompose x =

    let thousands x = (x/1000)*1000

    let hundreds  x = ((x - thousands x)/100)*100

    let tens      x =

        let temp = ((x - thousands x - hundreds x)/10)*10

        match temp with

        | _ when 0 < temp && temp < 20 -> 0

        | _                            -> temp

    let ones      x = (x - thousands x - hundreds x - tens x)

    (thousands x, hundreds x, tens x, ones x)

Just to show how this works, here is the decomposition of some example numbers:

> decompose 256;;

val it : int * int * int * int = (0, 200, 50, 6)

> decompose 512;;

val it : int * int * int * int = (0, 500, 0, 12)

The next piece of the puzzle is to take a section of the decomposed number and return the number of letters used to write it out:

let letters x =

    match x with

    |    1 ->  3    // one

    |    2 ->  3    // two

    |    3 ->  5    // three

    |    4 ->  4    // four

    |    5 ->  4    // five

    |    6 ->  3    // six

    |    7 ->  5    // seven

    |    8 ->  5    // eight

    |    9 ->  4    // nine

    |   10 ->  3    // ten

    |   11 ->  6    // eleven

    |   12 ->  6    // twelve

    |   13 ->  8    // thirteen

    |   14 ->  8    // fourteen

    |   15 ->  7    // fifteen

    |   16 ->  7    // sixteen

    |   17 ->  9    // seventeen

    |   18 ->  8    // eighteen

    |   19 ->  8    // nineteen

    |   20 ->  6    // twenty

    |   30 ->  6    // thirty

    |   40 ->  5    // forty

    |   50 ->  5    // fifty

    |   60 ->  5    // sixty

    |   70 ->  7    // seventy

    |   80 ->  6    // eighty

    |   90 ->  6    // ninety

    |  100 -> 10    // one hundred

    |  200 -> 10    // two hundred

    |  300 -> 12    // three hundred

    |  400 -> 11    // four hundred

    |  500 -> 11    // five hundred

    |  600 -> 10    // six hundred

    |  700 -> 12    // seven hundred

    |  800 -> 12    // eight hundred

    |  900 -> 11    // nine hundred

    | 1000 -> 11    // one thousand

    |    _ ->  0

Now, given a decomposed number, compute the total number of letters used to write it out:

let length (a, b, c, d) =

    match (a, b, c, d) with

    | (_, _, 0, 0) -> letters a + letters b

    | (0, 0, _, _) -> letters c + letters d

    | (_, _, _, _) -> letters a + letters b + letters c + letters d + 3

Note that the 3 in the last case is to account for the word "and".

Finally, to find out how many letters are needed to write out the numbers from one to 1000:

{1..1000} |> Seq.map decompose |> Seq.map length |> Seq.sum |> printfn "count = %A"

Problem 16 of Project Euler:

What is the sum of the digits of the number 2¹⁰⁰⁰?

This is a pretty straightforward problem to solve. I don't even get to draw any goofy pictures to make the logic clear.

The idea is to compute the number, split it into its individual digits and add them all together. Given a number n here's the F# code to extract each digit:

let rec digits n =

    seq {

        if n > 0I then

            yield n%10I

            yield! digits (n/10I)

    }

The desired answer is:

digits (2I ** 1000) |> Seq.sum |> printfn "%A"

Problem 15 of Project Euler asks:

Starting in the top left corner in a 20 by 20 grid, how many routes are there to the bottom right corner?

Without getting into complicated proofs I want to present some ideas that lead me to a solution for this problem.

First, suppose that we want to compute the number of paths from point A to point C and we know that all paths pass through point B. Say that a quick count shows there are m paths from A to B and n paths from B to C. Here is my goofy drawing showing the paths:

I think it's pretty obvious that there are mxn paths from A to C through B. Keep this in mind, it will be useful later on.

In general, when tackling a problem like this, I like to start by examining simpler versions of the same problem to see if I can come up with a general solution which I can then apply to the specific problem. I started by thinking about the problem on a 1 by 1 grid. After much thought it occurred to me that I could use a symmetry argument. If I consider all the midpoints along the line from the top right corner to the lower left corner then the number of paths from the start to the midpoint must be the same as the number of paths from the midpoint to the end. Here's a picture showing the paths to one of the midpoints in the grid:

In this diagram, the dotted line shows the axis of symmetry of the square. There are two midpoints (although I've only labelled the first). There is only one path from A to B and so there is one path from B to C. To find the total number of paths from A to C we add the number of paths from A to C through each midpoint. In this case it is 1² + 1² = 2. We are using the first idea mentioned above in this special, symmetrical, case.

Fair enough. Lets take a look at the 2 by 2 grid:

In this case we have three midpoints and I am just showing that there are two paths to the center one. By inspection we can see that there is only one path from A to the other midpoints and so the total number of path is 1² + 2² + 1² = 6. This, by the way, is a nice sanity check since it matches the answer given for a 2 by 2 grid in the statement of the problem.

I'm already starting to see a pattern here but just to really nail it down let's think about the 3 by 3 grid:

Now we have four midpoints to consider. As before, there is one path to each of the outside midpoints and we can see three paths to the each of the inside midpoints. We compute the total number of paths from A to C as 1² + 3² + 3² + 1² = 20.

Suppose we just look at the number of paths to the midpoints for each of the three cases I have presented here. Stacking them up in a provocative manner we see this:

   1 1 

  1 2 1

 1 3 3 1

To me this looks like Pascal's Triangle. And since I don't have to prove anything I'm just going to assume that it is. So, if I was to count the paths to each midpoint in a 4 by 4 grid I would see:

1 4 6 4 1

And so on. If I get the right answer then I'll know that my hunch is correct.

According to the Wikipedia article referenced above, the i^th entry in the n^th row of Pascal's triangle is:

So, with this formula, we have all the pieces necessary to solve this problem. To spell it out, the number of paths from the top left corner of the grid to the bottom right corner (with no backtracking) is sum of the squares of the coefficients of the n^th row of Pascal's triangle. To compute this number I use a couple of helper functions. Here is a factorial:

let factorial n = Seq.fold (*) 1I {1I..n}

Note that I could have used the built in factorial function but since I am using these problems as a learning experience it makes sense to write my own (or at least it seemed to make sense at the time!).

Here's how we can compute n choose i:

let choose n i = (factorial n)/((factorial i)*(factorial (n-i)))

With this we can create a sequence consisting of the entries in the n^th row of Pascal's triangle:

let pascal n =

    seq {

        for i in 0I..n -> choose n i

    }

The final step is to perform the actual computation:

pascal 20I |> Seq.map (fun i -> i*i) |> Seq.sum |> printfn "sum = %A"

Done!

Given the right insight the actual code was pretty easy.

Time for another Project Euler post. Problem fourteen asks you to:

Find the longest sequence using a starting number under one million.

Given a start number x_n the next entry in the sequence is given by:

You can find a nice discussion of this sequence in Wikipedia as the Collatz Problem. In F# this code will generate the sequence from a given starting value:

let rec sequence n =

    seq {

        yield n;

        match n with

        | 1ul -> ()

        | _ when n%2ul = 0ul -> yield! sequence (n/2ul)

        | _ when n%2ul = 1ul -> yield! sequence (3ul*n + 1ul)

    }

Then I can create a tuple containing the sequence starting value and the sequence length as follows:

let length n =

    (sequence n |> Seq.length, n)

So when I want to find the longest sequence I can write:

seq { 1ul..999999ul} |> Seq.map length |> Seq.max |> printfn "longest = %A"

Interestingly, the Seq.max function will just use the first entry in the tuple for its comparison. This means I don't have to clutter my code with special handling to extract the length information from the sequence of tuples. Obviously I got lucky when I put together the tuple and had the length come before the seed value.

The thirteenth Project Euler problem is:

Find the first ten digits of the sum of one-hundred 50-digit numbers.

I used F#'s built in BigInt data type. The array of 50 digit numbers looks like this:

let numbers = [|

    37107287533902102798797998220837590246510135740250I;

    46376937677490009712648124896970078050417018260538I;

    74324986199524741059474233309513058123726617309629I;

    91942213363574161572522430563301811072406154908250I;

    23067588207539346171171980310421047513778063246676I;

    89261670696623633820136378418383684178734361726757I;

    28112879812849979408065481931592621691275889832738I;

    44274228917432520321923589422876796487670272189318I;

    47451445736001306439091167216856844588711603153276I;

    70386486105843025439939619828917593665686757934951I;

    62176457141856560629502157223196586755079324193331I;

    64906352462741904929101432445813822663347944758178I;

    92575867718337217661963751590579239728245598838407I;

    58203565325359399008402633568948830189458628227828I;

    80181199384826282014278194139940567587151170094390I;

    35398664372827112653829987240784473053190104293586I;

    86515506006295864861532075273371959191420517255829I;

    71693888707715466499115593487603532921714970056938I;

    54370070576826684624621495650076471787294438377604I;

    53282654108756828443191190634694037855217779295145I;

    36123272525000296071075082563815656710885258350721I;

    45876576172410976447339110607218265236877223636045I;

    17423706905851860660448207621209813287860733969412I;

    81142660418086830619328460811191061556940512689692I;

    51934325451728388641918047049293215058642563049483I;

    62467221648435076201727918039944693004732956340691I;

    15732444386908125794514089057706229429197107928209I;

    55037687525678773091862540744969844508330393682126I;

    18336384825330154686196124348767681297534375946515I;

    80386287592878490201521685554828717201219257766954I;

    78182833757993103614740356856449095527097864797581I;

    16726320100436897842553539920931837441497806860984I;

    48403098129077791799088218795327364475675590848030I;

    87086987551392711854517078544161852424320693150332I;

    59959406895756536782107074926966537676326235447210I;

    69793950679652694742597709739166693763042633987085I;

    41052684708299085211399427365734116182760315001271I;

    65378607361501080857009149939512557028198746004375I;

    35829035317434717326932123578154982629742552737307I;

    94953759765105305946966067683156574377167401875275I;

    88902802571733229619176668713819931811048770190271I;

    25267680276078003013678680992525463401061632866526I;

    36270218540497705585629946580636237993140746255962I;

    24074486908231174977792365466257246923322810917141I;

    91430288197103288597806669760892938638285025333403I;

    34413065578016127815921815005561868836468420090470I;

    23053081172816430487623791969842487255036638784583I;

    11487696932154902810424020138335124462181441773470I;

    63783299490636259666498587618221225225512486764533I;

    67720186971698544312419572409913959008952310058822I;

    95548255300263520781532296796249481641953868218774I;

    76085327132285723110424803456124867697064507995236I;

    37774242535411291684276865538926205024910326572967I;

    23701913275725675285653248258265463092207058596522I;

    29798860272258331913126375147341994889534765745501I;

    18495701454879288984856827726077713721403798879715I;

    38298203783031473527721580348144513491373226651381I;

    34829543829199918180278916522431027392251122869539I;

    40957953066405232632538044100059654939159879593635I;

    29746152185502371307642255121183693803580388584903I;

    41698116222072977186158236678424689157993532961922I;

    62467957194401269043877107275048102390895523597457I;

    23189706772547915061505504953922979530901129967519I;

    86188088225875314529584099251203829009407770775672I;

    11306739708304724483816533873502340845647058077308I;

    82959174767140363198008187129011875491310547126581I;

    97623331044818386269515456334926366572897563400500I;

    42846280183517070527831839425882145521227251250327I;

    55121603546981200581762165212827652751691296897789I;

    32238195734329339946437501907836945765883352399886I;

    75506164965184775180738168837861091527357929701337I;

    62177842752192623401942399639168044983993173312731I;

    32924185707147349566916674687634660915035914677504I;

    99518671430235219628894890102423325116913619626622I;

    73267460800591547471830798392868535206946944540724I;

    76841822524674417161514036427982273348055556214818I;

    97142617910342598647204516893989422179826088076852I;

    87783646182799346313767754307809363333018982642090I;

    10848802521674670883215120185883543223812876952786I;

    71329612474782464538636993009049310363619763878039I;

    62184073572399794223406235393808339651327408011116I;

    66627891981488087797941876876144230030984490851411I;

    60661826293682836764744779239180335110989069790714I;

    85786944089552990653640447425576083659976645795096I;

    66024396409905389607120198219976047599490197230297I;

    64913982680032973156037120041377903785566085089252I;

    16730939319872750275468906903707539413042652315011I;

    94809377245048795150954100921645863754710598436791I;

    78639167021187492431995700641917969777599028300699I;

    15368713711936614952811305876380278410754449733078I;

    40789923115535562561142322423255033685442488917353I;

    44889911501440648020369068063960672322193204149535I;

    41503128880339536053299340368006977710650566631954I;

    81234880673210146739058568557934581403627822703280I;

    82616570773948327592232845941706525094512325230608I;

    22918802058777319719839450180888072429661980811197I;

    77158542502016545090413245809786882778948721859617I;

    72107838435069186155435662884062257473692284509516I;

    20849603980134001723930671666823555245252804609722I;

    53503534226472524250874054075591789781264330331690I;

    |]

Once we have the numbers in an array it is easy to add them up:

numbers |> Array.sum |> printfn "total = %A"

So easy it feels like cheating. Formatting the data was the most time consuming part of solving this problem.

I think I'll post my answer to Problem 12 of Project Euler next:

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

What exactly is a triangle number? It is a number given by this formula:

For this problem we actually want to generate a sequence of triangular numbers. It is most efficient to generate the n^th number by adding n to the (n-1)^th number. We can capture that in F# as follows:

let triangles =

    let rec inner m n =

        seq {

            let p = m + n

            yield p

            yield! inner p (n+1L)

        }

    seq {

        yield! inner 0L 1L

    }

The next question is how do we find the number of divisors of the number n? The brute force way would be to try all values up to n/2 and find the ones which divide evenly into the target number. I figured there had to be a better way. A little bit of research on the web showed me that if you know the prime factorization for a number then you can easily compute the number of divisors. Suppose that the a number n can be factored as follows:

where ω(n) is the number of distinct prime factors of n. Then the number of divisors of n is given by:

Fortunately, I'm no longer a math major and so I don't have to prove that this is true. For my purposes it's enough to just try it out and see if it seems reasonable. So, we can say that 12 = 2² x 3. According to the formula then the number of divisors should be (2 + 1) x (1 + 1) = 6. Listing the divisors we find the following set: {1, 2, 3, 4, 6, 12}. Good enough for me.

To solve this Project Euler problem I need to be able to factor an arbitrary number. The simplest thing is to try all the primes up to the square root of the target number and see if they divide into it evenly. As before I need a sequence of primes. Here's the code:

let primes =

    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)

        }

    sieve 2L Map.empty

I think this will be the last time I display this code. From now on just assume that when needed I can generate a list of primes with a lazy sequence.

I also need a few helper functions. Here's an integer square root:

let root (n:int64) =

    int64 (sqrt (float n))

And an integer power function (tail recursive just for fun):

let power p e  =

    let rec inner p e acc =

        if e = 0L then

            acc

        else

            inner p (e-1L) p*acc

    inner p e 1L

Finally here's a function to compute the number of times a prime divides into a number:

let exponent p n =

    let rec inner p n acc =

        if n%p = 0L then

            inner p (n/p) (acc+1L)

        else

            acc

    inner p n 0L

The next step is to compute a sequence which consists of the prime factors of a number along with each factor's exponent. Here's what I came up with:

let rec factors n =

    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 = exponent k n

                yield (k, e)

                yield! factors (n/power k e)

    }

There's some inefficiency here in that each time we recurse we start the sequence of primes over again. I'm not losing any sleep over it...

Given all this we can now compute the number of divisors of a number:

let divisors n =

    factors n

        |> Seq.map (fun k -> snd k)

        |> Seq.fold (fun acc e -> acc*(e+1L)) 1L

The final step is to find the first entry in the triangular number sequence which has more than 500 divisors.

triangles

    |> Seq.tryFind (fun n -> divisors n > 500L)

    |> printfn "%A"

It was nice to revisit this solution while putting together this post. My original solution was a bit hard to read. I think it is much clearer now.

Problem 11 of the Project Euler:

What is the greatest product of four numbers on the same straight line in the 20 by 20 grid?

The problem description lays out an array of two digit integers where any group of 4 numbers either horizontally or vertically or diagonally must be evaluated to find the maximum product. My algorithms kept banging into the edges of the array (which caused ugly code) until I realised that I could embed their array inside a larger one and then just look at a subset of it. The four numbers I would look at in each step could include the padding and, if the pad number was zero, the maximum product would not be affected. I'd have to examine a few more products than if I handled the edges correctly but that's just a minor inefficiency.

With that in mind I defined the array as follows. The actual data starts at (3,3):

let data =

    [|

      [|00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00|];

      [|00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00|];

      [|00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00|];

      [|00; 00; 00; 08; 02; 22; 97; 38; 15; 00; 40; 00; 75; 04; 05; 07; 78; 52; 12; 50; 77; 91; 08; 00; 00; 00|];

      [|00; 00; 00; 49; 49; 99; 40; 17; 81; 18; 57; 60; 87; 17; 40; 98; 43; 69; 48; 04; 56; 62; 00; 00; 00; 00|];

      [|00; 00; 00; 81; 49; 31; 73; 55; 79; 14; 29; 93; 71; 40; 67; 53; 88; 30; 03; 49; 13; 36; 65; 00; 00; 00|];

      [|00; 00; 00; 52; 70; 95; 23; 04; 60; 11; 42; 69; 24; 68; 56; 01; 32; 56; 71; 37; 02; 36; 91; 00; 00; 00|];

      [|00; 00; 00; 22; 31; 16; 71; 51; 67; 63; 89; 41; 92; 36; 54; 22; 40; 40; 28; 66; 33; 13; 80; 00; 00; 00|];

      [|00; 00; 00; 24; 47; 32; 60; 99; 03; 45; 02; 44; 75; 33; 53; 78; 36; 84; 20; 35; 17; 12; 50; 00; 00; 00|];

      [|00; 00; 00; 32; 98; 81; 28; 64; 23; 67; 10; 26; 38; 40; 67; 59; 54; 70; 66; 18; 38; 64; 70; 00; 00; 00|];

      [|00; 00; 00; 67; 26; 20; 68; 02; 62; 12; 20; 95; 63; 94; 39; 63; 08; 40; 91; 66; 49; 94; 21; 00; 00; 00|];

      [|00; 00; 00; 24; 55; 58; 05; 66; 73; 99; 26; 97; 17; 78; 78; 96; 83; 14; 88; 34; 89; 63; 72; 00; 00; 00|];

      [|00; 00; 00; 21; 36; 23; 09; 75; 00; 76; 44; 20; 45; 35; 14; 00; 61; 33; 97; 34; 31; 33; 95; 00; 00; 00|];

      [|00; 00; 00; 78; 17; 53; 28; 22; 75; 31; 67; 15; 94; 03; 80; 04; 62; 16; 14; 09; 53; 56; 92; 00; 00; 00|];

      [|00; 00; 00; 16; 39; 05; 42; 96; 35; 31; 47; 55; 58; 88; 24; 00; 17; 54; 24; 36; 29; 85; 57; 00; 00; 00|];

      [|00; 00; 00; 86; 56; 00; 48; 35; 71; 89; 07; 05; 44; 44; 37; 44; 60; 21; 58; 51; 54; 17; 58; 00; 00; 00|];

      [|00; 00; 00; 19; 80; 81; 68; 05; 94; 47; 69; 28; 73; 92; 13; 86; 52; 17; 77; 04; 89; 55; 40; 00; 00; 00|];

      [|00; 00; 00; 04; 52; 08; 83; 97; 35; 99; 16; 07; 97; 57; 32; 16; 26; 26; 79; 33; 27; 98; 66; 00; 00; 00|];

      [|00; 00; 00; 88; 36; 68; 87; 57; 62; 20; 72; 03; 46; 33; 67; 46; 55; 12; 32; 63; 93; 53; 69; 00; 00; 00|];

      [|00; 00; 00; 04; 42; 16; 73; 38; 25; 39; 11; 24; 94; 72; 18; 08; 46; 29; 32; 40; 62; 76; 36; 00; 00; 00|];

      [|00; 00; 00; 20; 69; 36; 41; 72; 30; 23; 88; 34; 62; 99; 69; 82; 67; 59; 85; 74; 04; 36; 16; 00; 00; 00|];

      [|00; 00; 00; 20; 73; 35; 29; 78; 31; 90; 01; 74; 31; 49; 71; 48; 86; 81; 16; 23; 57; 05; 54; 00; 00; 00|];

      [|00; 00; 00; 01; 70; 54; 71; 83; 51; 54; 69; 16; 92; 33; 48; 61; 43; 52; 01; 89; 19; 67; 48; 00; 00; 00|];

      [|00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00|];

      [|00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00|];

      [|00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00; 00|];

    |]

Now that I have the array data in place I can create a sequence consisting of the products of four numbers in various directions. I chose to look at the four numbers to the left of the current number, down diagonally to the left, straight down and down diagonally to the right. See the coloured numbers in the above array. That's all that's needed to get every possibly four digit product. I created a sequence from the array data that looked like this:

let products =

    seq {

        for i in 3..22 do

            for j in 3..22 do

                yield (data.[i].[j] * data.[i-1].[j  ] * data.[i-2].[j  ] * data.[i-3].[j  ])

                yield (data.[i].[j] * data.[i-1].[j+1] * data.[i-2].[j+2] * data.[i-3].[j+3])

                yield (data.[i].[j] * data.[i  ].[j+1] * data.[i  ].[j+2] * data.[i  ].[j+3])

                yield (data.[i].[j] * data.[i+1].[j+1] * data.[i+2].[j+2] * data.[i+3].[j+3])

    }

A simple chain of operations finds the maximum product:

products |> Seq.max |> printfn "maximum = %A"

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.

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.