4. The Rite of Spring by Igor Stravinsky
5. Playing the Fool by Gentle Giant

In the summer of 1978, probably August, I spent two or three weeks staying with an aunt and uncle who lived in Québec City. While there I purchased a couple of LP's. By extreme coincidence they were the two listed above. I'm not sure which orchestra did The Rite of Spring. I checked this morning and, sadly, neither disk remains in my collection. I do, however, have multiple versions of the former on my iPod as well as the reissued version of the latter.

In 1978 I was nineteen and had just completed my second year of university at UBC. My first year was at the Okanagan College in Vernon, BC. During that first year I took a music appreciation course which must have been where I was exposed to the The Rite of Spring. The teacher of that course suggested to my friends and me that we might like Gentle Giant. How had she even heard of them? But was she ever right (for me at least, my friends stuck with big band jazz). I loved Gentle Giant's brand of pretentious progressive rock. They referenced RD Laing. They played the harpsichord. They were gargantuan. They played recorders. What was not to like?

I never got to see Gentle Giant in concert. By the time I was interested in them they weren't touring North America and were pretty close to breaking up. One of the guys at my residence at UBC had seen them play live. I think they opened for Yes. He said something about a dragon.

A decade or so ago (in the late 1990's) I went up to San Francisco and saw the San Francisco Symphony Orchestra perform The Rite of Spring under Michael Tilson Thomas. I thought I knew that piece of music very well. I'd listened to it regularly for twenty-plus years. I knew that piece of music. Ha! Watching it performed live was a revelation. I had chosen a seat close to and in front of the cello section (I was learning to play the cello at the time). I remember my astonishment at seeing waves of melody passing through the cellos. Themes that I had never consciously heard before.

I suppose I should make some pithy observation about those who listen but do not hear. But who would benefit from that? Me? I think not.

3. London Calling by The Clash

I generally have my iPod cycle through albums at random. On 5 January 2010 I read a posting on Boing Boing about the 30th anniversary of the release of London Calling in the United States. That very day it started playing during my commute home from work. It came up again today. I am listening to it as I write this post. Personally I think it has aged very well. It sounds as fresh to me now as it did when I first heard it on CITR when I was a student at UBC in the early 1980's.

Update:

I forgot to mention that I actually got to see The Clash in concert. They opened for The Who during their first farewell tour in 1982. This was at the Kingdome in Seattle. A friend drove a bunch of us there for the concert. I was not a fan of The Who but my friends were. I had heard that The Clash were opening and that was good enough for me. We were miles away from the stage. Up in the cheap seats. Even the huge screen showing the performers was too small to be seen clearly. The first act was, interestingly, T-Bone Burnett. I had never heard of him before. Later on he became a favourite. In part because he married Sam Phillips whose music I loved. I don't really remember anything about the music. I believe I was underwhelmed by The Clash. I didn't know much about The Who. Their music washed over me and left me unchanged.

Since that was The Who's farewell tour and I was there I have refused to listen to anything by them ever since! Why no, I am not vindictive. Heh.

Final Update:

While I remember almost nothing of the music I heard at this concert I do have two very clear memories: one from before and one from after. My friends and I were living at a residence on UBC and the cook there made us all bag lunches to take along since we would be missing dinner that night. We waited outside the Kingdome for a while until they started letting people in. When it was our turn to enter, the security guards made us throw away our bag lunches. I guess they were worried we'd start throwing apples at people or something.

On the way home, being quite hungry, we stopped at a restaurant that was open late. Some kind of fast food place no doubt. In the parking lot was one of my friends from Vancouver who had also driven down for the concert.

Those two memories are very clear. Everything is kind of blurry.

I still like the Clash...

2. Six Suites for Unaccompanied Cello by J.S. Bach

Judging by the number of recordings I have of this music it is a huge favourite. Here's a partial list of artists I have playing one or more of these suites:

• Jaap Ter Linden
• Pieter Wispelwey
• Anner Bylsma
• Pablo Casals
• Janos Starker
• Yo-Yo Ma

When I studied cello (for about five years in the mid-1990's) I reached the point where I could play most of the first suite. I still have the score Karen, my teacher, and I marked up with fingerings. I don't think I will ever have the time or energy to pick up the cello again but I am happy that I was able play some of them.

My cello teacher played the first suite during my wedding reception at Hakone Gardens in Los Gatos in August of 1999. My wife and I had wanted her to play it while we were walking down the aisle but, well, I don't really remember why she wouldn't.

1. Historical Anthology of Music (volume I).

I have played the recorder off and on since the mid-1990's. I first encountered this book when my teacher gave me a photocopy of piece number 49, a madrigal by Jacopo da Bolgna. I believe I first played it on a Baroque alto in F.

When my son was born (premature and with issues) in 2000 I dropped almost all my extracurricular activities. I stopped all music lessons (voice, cello and recorder) and focussed on my wife and child. I started playing the recorder again about three years ago and restarted lessons (with the same teacher; thanks Carol!) about six months ago.

I recently revisited this piece. This time I am using it to learn to play the Renaissance alto in G. I also now have my own copy of this book and it is full of wonderful music. I suspect it will be a constant companion during 2011.

John Terauds of The Star proposes this challenge:

What are the 100 pieces of music -- including full operas -- that we simply can't live without in 2011?

Thinking about the music that I love has brought back memories of where I was and what I was doing when I first heard the piece. I'm going to try and capture some of those thoughts in a few posts about the music that is meaningful to me.

I'm going to talk about the music that I play (or have played) as an amateur musician as well as the music that I listen to a lot on my iPod. I will not claim that the music I discuss is significant to anyone other than me. Also, I will be listing pieces in the order in which they occur to me.

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. 

Many years ago, in the late 1980's when I lived in Vancouver, I took night school courses to learn conversational Cantonese. I stopped when I moved to California in 1990 and all the vocabulary I acquired then has been slowly eroding. However, my previous post on spelling out numbers reminded me of how one says the numbers in Cantonese. Here are the numbers with their Chinese character and Cantonese pronunciation:

one	1	yāt	一
two	2	yih	二
three	3	sàam	三
four	4	sei	四
five	5	nǵh	五
six	6	luhk	六
seven	7	chat	七
eight	8	baat	八
nine	9	gáu	九
ten	10	sahp	十

Subsequent numbers are built logically from these. For example, eleven (11) is 十一, twenty-one (21) is 二十一 and thirty-one (31) is 三十一. In Cantonese you would pronounce 三十一 as "sàam sahp yāt".


A program equivalent to the one that counted the number of letters in the English spelling of numbers would probably count the strokes of the Chinese numbers. It would have a lot fewer special cases than the program I presented earlier.


I use the Yale Romanization for Cantonese since that's what I was taught long back. There is a lot of information about Chinese Numerals online.

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"

Speaking of goofy pictures, I should mention that I used AutoDesk // LABS' Project // DRAW to create the images in my previous post.

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.

In the previous post I talked about computing the divisor count of an integer given its prime factorisation. I also said I didn't have to prove the formula. However, I find I just can't leave it alone. Why does that formula work? I will sketch out a brief rationale (that falls far short of a formal proof).

I think it is pretty clear that any number that divides into an integer has to be made up from some mixture of the integer's prime factors. In other words, if a prime p is not a factor of an integer n then p will not divide evenly into n. Also, if p is a factor of n and it divides into n α times then p^α+1 will not divide into n.

So we can build all divisors of n by multiplying all of n's prime factors and varying their exponent between 0 and allowed exponent for that factor. That is, for example, if n = p^α then the divisors of n are {p⁰, p¹, ... p^α}. This set contains α+1 elements.

If n has multiple prime factors then you can multiply them all together with all acceptable permutations of exponents to create the divisor set. In the field of combinatorics this is selection with replacement (pdf). If you have ω(n) positions to fill and each position has α_i+1 possible values then the total number of combinations is the product of the possible values. That is:

I feel much better now.

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.

The next book I want to classify is "Editing Early Music" by John Caldwell (ISBN 0-19-816544-7). The publisher, Oxford University Press, claims that it should be filed at 781.2 Elements of Music but this doesn't seem right to me. I'm not sure which edition of the Dewey Classification system they are using, perhaps this number made more sense before the great revision. Looking at the revised 780 schedule I found 780.149 Editing Music. This seems to be what I'm looking for. There is also a notation under 780.148 Musical notation, abbreviations, symbols that says "Class transcription from one form of notation to another in 780.149". This again indicates that I'm on the right track (this book is all about transforming early music notation into a form that can be understood by a modern musician).

At this point I could easily say I'm done but the Dewey system does allow further aspects or facets of the work to be captured. In this case I think you could add 780.902 500 - 1449 to the base number to indicate that the book is relevant to early music. The way you do this is still a little unclear to me. I believe that you add a facet indicator (which can be either 0 or 1) to the base number and then append the number after the decimal place in the facet. This would yield either 780.1490902 or 780.1491902. I think this would be fine in a large library with a ton of books on the same subject. For my library (with maybe 50 books and a couple hundred scores) I think it makes more sense to use the shorter classification number and then sort all the books in that class by Cutter number.

Update: I now believe that 0 is used to indicate a standard subdivision and 1 indicates a facet for a built number. So, in the example above, 780.1491902 would be the correct class number. Please correct me if I'm wrong.