Tuesday, March 17, 2009

I’m back again with another post in my Yet Another Project Euler Series. As a reminder, my approach to these problems is to try to find the most beautiful solution that I can using F#. While performance is an important, I’m looking for the most elegant solutions that I can find.

Project Euler problem ten is another dealing prime numbers.

The sum of the primes below 10 is 2 + 3 + 5 + 7 = 17.

Find the sum of all the primes below two million.

As you can guess, there’s not too much to this problem. In fact, we’ll declare just one simple supporting function.

let lessThan y x = x < y

Using the above function and our existing prime number generator, our solution is short and simple.

primes
  |> Seq.take_while (lessThan 2000000L)
  |> Seq.sum

I feel a little guilty that this post is so short, but that’s the benefit of spending two whole blog posts working out a prime number generator for problem seven.

posted on Tuesday, March 17, 2009 2:13:48 AM (Pacific Standard Time, UTC-08:00)  #    Comments [4]

kick it on DotNetKicks.com
 Monday, March 09, 2009

Welcome back F# junkies! In this installment of YAPES1, we’re tackling Project Euler problem nine.

A Pythagorean triplet is a set of three natural numbers, a < b < c, for which,

a2 + b2 = c2

For example, 32 + 42 = 9 + 16 = 25 = 52.

There exists exactly one Pythagorean triplet for which a + b + c = 1000. Find the product abc.

Like we’ve done before, we’ll break this problem down into steps:

  1. Generate a stream of Pythagorean triples.
  2. Find the Pythagorean triple whose sum is equal to 1000.
  3. Calculate the product of that triple.

Sound good?

The first order of business is to determine how to generate Pythagorean triples. It turns out that there are many ways to do this. However, to keep it simple, we will use Euclid’s classic formula:

Given a pair of positive integers m and n with m > n

a = 2mn, b = m2n2, c = m2 + n2

Truthfully, this formula doesn’t actually produce every possible Pythagorean triple. However, since the solution is within the set that it generates, that really isn’t a problem.

NeRd Note
To produce every Pythagorean triple, we would have to introduce an additional parameter, k, into the formula.
a = 2kmn, b = k(m2n2), c = k(m2 + n2)
This is left as an exercise for you, dear reader!

First, we’ll declare an F# function that calculates a Pythagorean triple using Euclid’s formula.

let getTriple m n =
  let a = 2*m*n
  let b = m*m - n*n
  let c = m*m + n*n
  a,b,c

Using the above function, it’s easy to produce a stream of Pythagorean triples. Our strategy will be to start with the pair, m = 2 and n = 1. To generate the next pair, we’ll first check if n + 1 is less than m. If it is, we’ll use m and n + 1 as the next pair, otherwise we’ll use m + 1 and 1. Using Seq.unfold, we can use this algorithm to produce an infinite stream of Pythagorean triples:

let triples =
  Seq.unfold (fun (m,n) ->
                let res = getTriple m n
                let next = if n+1 < m then m,n+1 else m+1,1
                Some(res, next))
             (2,1)

With our generator in place, we can easily solve the problem.

triples
  |> Seq.find (fun (a,b,c) –> a+b+c = 1000)
  |> (fun (a,b,c) –> a*b*c)

While the above solution works well, we can abstract further to produce a more elegant solution. First, we’ll declare three functions to handle the three operations that we’re performing explicitly, sum, product and equals.

let sum (a,b,c) = a+b+c
let product (a,b,c) = a*b*c
let equals x y = x = y

Now, we can restate our solution using these functions with a bit of function composition to produce a truly beautiful solution.

triples
  |> Seq.find (sum >> equals 1000)
  |> product

As we’ve seen once before, function composition is beautiful.

1Yet Another Project Euler Series

posted on Monday, March 09, 2009 7:05:47 AM (Pacific Standard Time, UTC-08:00)  #    Comments [1]

kick it on DotNetKicks.com
 Saturday, March 07, 2009

Welcome coding ninjas!1 Continuing with Yet Another Project Euler Series (YAPES™), we come upon problem eight.

Find the greatest product of five consecutive digits in the 1000-digit number.

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

Being one of the earlier Project Euler problems, this is pretty trivial to solve. We’ll approach the problem by first breaking it into a set of simple steps:

  1. Transform the 1000-digit number into a list of digits.
  2. Transform the list of digits into a list of quintuplets—that is, a list containing tuples of each five consecutive elements.
  3. Transform the list of quintuplets into a list containing the product of each.
  4. Find the maximum element in the list.

Some or all of the steps above could be combined to produce a more optimal solution. However, the beauty of F#—and functional programming in general—is the ability to easily break a problem like this into a series of data transformations that run quickly enough for most situations. When a problem is broken down into simple operations, it is easier to optimize in other ways (e.g. parallelization).

The first step is to transform the 1000-digit number into a list of digits. To facilitate working with such a large number in code, we’ll represent it as a multi-line string. Like before, we’ll break the problem of converting a string into a list of digits into a set of simple steps:

  1. Transform the string into a list of chars.
  2. Because this was a multi-line string, some of the chars will be line breaks and whitespace. So, we’ll filter the list to include only the chars that represent digits.
  3. Transform the list of digit chars into a list containing the numerical value of each.

Let’s define a few library functions to to help with each step above.

module String =
  /// Takes a string and produces a list of chars.
  let toChars (s : string) =
    s.ToCharArray() |> Array.to_list

module Char =
  /// Determines whether a char represents a digit.
  let isDigit (c : char) =
    System.Char.IsDigit(c)

  /// Converts a char representing a digit into its numerical value.
  let toNumber (c : char) =
    int c - int '0'

Armed with these, we can perform the 3 steps above in a declarative fashion.

let number =
  "73167176531330624919225119674426574742355349194934
   96983520312774506326239578318016984801869478851843
   85861560789112949495459501737958331952853208805511
   12540698747158523863050715693290963295227443043557
   66896648950445244523161731856403098711121722383113
   62229893423380308135336276614282806444486645238749
   30358907296290491560440772390713810515859307960866
   70172427121883998797908792274921901699720888093776
   65727333001053367881220235421809751254540594752243
   52584907711670556013604839586446706324415722155397
   53697817977846174064955149290862569321978468622482
   83972241375657056057490261407972968652414535100474
   82166370484403199890008895243450658541227588666881
   16427171479924442928230863465674813919123162824586
   17866458359124566529476545682848912883142607690042
   24219022671055626321111109370544217506941658960408
   07198403850962455444362981230987879927244284909188
   84580156166097919133875499200524063689912560717606
   05886116467109405077541002256983155200055935729725
   71636269561882670428252483600823257530420752963450"
   |> String.toChars
   |> List.filter Char.isDigit
   |> List.map Char.toNumber

The next step is to transform the list of digits into a list of quintuplets. To achieve this, we’ll write a simple recursive function.

let rec toQuintuplets l =
  match l with
  | x1::(x2::x3::x4::x5::_ as t) -> (x1,x2,x3,x4,x5)::(toQuintuplets t)
  | _ -> []

Notice that the first pattern match above is slightly more complicated as we bind the list starting with the next element in the list to t, to make the recursion cleaner.

NeRd Note
The toQuintuplets function defined above will work well for this particular problem but will stack overflow if passed too large of a list because it is head-recursive. There are a couple of ways to make it tail-recursive and avoid blowing the stack, but the most elegant is to employ a continuation.
let toQuintuplets l =
  let rec loop l cont =
    match l with
    | x1::(x2::x3::x4::x5::_ as t) -> loop t (fun l –> cont ((x1,x2,x3,x4,x5)::l))
    | _ -> cont []

  loop l (fun x -> x)
Of course, that unnecessarily complicates this solution by turning the logic inside out!

The last function we’ll define is a small helper to produce the product of all of the values in a quintuplet.

let product (x1,x2,x3,x4,x5) = x1*x2*x3*x4*x5

All of the pieces are now in place, and the steps can be combined to solve Project Euler problem eight.

number
  |> toQuintuplets
  |> List.map product
  |> List.max

It’s that simple!

1Or pirates.

posted on Saturday, March 07, 2009 1:09:11 PM (Pacific Standard Time, UTC-08:00)  #    Comments [2]

kick it on DotNetKicks.com
 Monday, January 19, 2009

Last time, we took a brute force approach to solving Project Euler problem seven. Unfortunately, the resulting prime number generator turned out to be fairly ugly and not really efficient enough to handle the needs of future YAPES1 problems. This time, we’ll approach the problem using a tried-and-true algorithm that should give us the desired performance and produce a beautiful solution. I can already hear some of you rolling your eyes, because we’ll be looking at...

The Sieve of Eratosthenes

Amazingly, the Sieve of Eratosthenes is perhaps the oldest known algorithm for generating prime numbers, and yet it remains reasonably efficient by today’s standards. It’s true that there are faster methods available, but this should be fast enough. The basic algorithm is nothing earth-shattering—it’s first semester computer science stuff. However, we’ll add a twist to keep things interesting.

The Sieve of Eratosthenes works by marking composite numbers:

  1. Start with the natural numbers from 2 up to some maximum limit.
  2. Take the first unmarked number (2) and mark all multiples of that number. Note that 2 remains unmarked. Only its multiples (e.g. 4,6,8,…) are marked.
  3. Take the next unmarked number and mark all of its multiples.
  4. Repeat step 3 until there are no more numbers left to test.
  5. When finished, the unmarked numbers are the primes.

As described, this algorithm doesn’t quite work for our purposes. It is certainly easy to implement, but because of the starting requirement (a big list of numbers), it doesn’t lend itself well to generating a, potentially, infinite stream of primes. However, with a small modification, we can make it work.

The adaptation we’ll implement comes from from the excellent paper, The Genuine Sieve of Eratosthenes2, by Melissa E. O’Neill. The idea is to turn the algorithm inside out. Instead of keeping track of which numbers are prime, we’ll maintain a table of composite numbers along with their known factors. As we iterate the natural numbers, the table will be updated with new information about which numbers are known to be composite and what factors have been found.

Confused? Let’s walk through the algorithm, step by step.

Tracking Composites

We start with 2 as the first natural number and check to see if it exists in the composite table. Since the table is empty, 2 isn’t present, so we add a new list containing 2 to the table with the square of 2 (4) as the key. Finally, since 2 was not found in the composite table, we know it is prime.

Sieve Diagram 1 (after 2)

We take the next natural number, 3. Like before, 3 is not found in the composite table, so we add it to the table using its square (9) as the key and return 3 as prime.

Sieve Diagram 2 (after 3)

We take the next natural number, 4. This time is different—4 is found in the composite table, so we know it’s not prime. In this case, we take the list found in the table under 4 (the known factors), and for each factor, we add 4 and insert the result into the table. Thus, we insert 6 into the table along with its known factor, 2. (Interestingly, this algorithm doesn’t tell us that 3 is a factor of 6). Finally, we remove 4 from the table.

Sieve Diagram 3 (after 4)

Let's skip ahead to some more interesting cases.

After 9 has been processed, the table looks like so:

Sieve Diagram 4 (after 9)

The next natural number, 10, is found in the table, so we know that it’s composite and not prime. Again, we take the list of known factors, and for each factor, we add 10 and insert the result into the table. This time, the value (12) already exists in the table, so we add it to the list of known factors for that value. Lastly, we remove 10 from the table.

Sieve Diagram 5 (after 10)

The next natural number, 11, is not found in the table, so we insert its square and return it as prime.

Sieve Diagram 6 (after 11)

We find the next natural number, 12, in the table. However, because there is more than one known factor in the table, we add two new composites, 12+2=14 and 12+3=15, before finally removing 12.

Sieve Diagram 7 (after 12)

As you can see, this take on the classic algorithm allows us to determine which numbers are prime and composite by using a minimal amount of memory. Obviously, the table can grow large as we find larger primes, but the memory usage is quite reasonable relative to the traditional sieve where a large list of natural numbers would have to be generated before even the first prime is returned.

Show Me Some Code!

The F# code nearly writes itself.

let reinsert x table prime =
   let comp = x+prime
   match Map.tryfind comp table with
   | None        -> table |> Map.add comp [prime]
   | Some(facts) -> table |> Map.add comp (prime::facts)

let rec sieve x table =
  seq {
    match Map.tryfind x table with
    | None ->
        yield x
        yield! sieve (x+1L) (table |> Map.add (x*x) [x])
    | Some(factors) ->
        yield! sieve (x+1L) (factors |> List.fold_left (reinsert x) (table |> Map.remove x)) }

let primes =
  sieve 2L Map.empty

The meat of the algorithm is in the sieve function. Here, we try to locate x in the table. If it is not found, we yield x (because it’s prime) and recursively call sieve passing x+1 and a table with the square of x added to it. (Note that Map is an immutable data structure, so “mutating” functions like Map.add actually return a brand new Map.) If x is found in the table, we recursively call sieve passing x+1 and a table from which x has been removed and each known factor has been reinserted. The reinsert function performs the trivial work of reinserting a known factor into the table.

A couple of other interesting details about this code are worth pointing out:

  • sieve returns a sequence expression, denoted by the seq keyword and the curly braces. Sequence expressions are lazily evaluated, meaning that the next prime won't be calculated until requested. In addition, the sequence expression is recursive due to the use of the yield! keyword, which can be used to return another sequence expression as part of a sequence.
  • The call to List.fold_left might look a bit exotic because we are folding over the list of known factors to produce a new table with each factor reinserted.

Mutation Isn't All Bad

Now that we have a working prime number generator, let’s take a moment to optimize. Until now, we’ve been sticking with the immutable F# Map type, but we could employ a mutable data structure to improve performance. To replace the data structure in a way that preserves the beauty of the algorithm, we’ll define a set of APIs for Dictionary<TKey,TValue> that mimic the Map APIs:

module Dict =
  open System.Collections.Generic

  let empty() = new Dictionary<_,_>()

  let add k v (d : #IDictionary<'a,'b>) =
    d.[k] <- v; d

  let remove (k : 'a) (d : #IDictionary<'a,'b>) =
    d.Remove(k) |> ignore; d

  let tryfind k (d : #IDictionary<'a,'b>) =
    match d.TryGetValue(k) with
    | true, v  -> Some(v)
    | false, _ -> None

Now, it’s a simple matter of replacing references to Map with Dict, and our algorithm immediately benefits from the performance boost of using a mutable table.

let reinsert x table prime =
   let comp = x+prime
   match Dict.tryfind comp table with
   | None        -> table |> Dict.add comp [prime]
   | Some(facts) -> table |> Dict.add comp (prime::facts)

let rec sieve x table =
  seq {
    match Dict.tryfind x table with
    | None ->
        yield x
        yield! sieve (x+1L) (table |> Dict.add (x*x) [x])
    | Some(factors) ->
        yield! sieve (x+1L) (factors |> List.fold_left (reinsert x) (table |> Dict.remove x)) }

let primes =
  sieve 2L (Dict.empty())

With our prime number generator in hand, solving Project Euler problem seven is trivial:

primes |> Seq.nth 10000

At last, we’ve produced a prime number generator with a good balance of performance and beauty. If we’ve learned anything from this article or its predecessor, I hope it’s that finding an optimal solution often requires finding an appropriate algorithm.

Take that Chris!

1Yet Another Project Euler Series
2O'Neill's paper appears in The Journal of Functional Programming, Volume 19, Issue 01.

posted on Monday, January 19, 2009 11:44:09 PM (Pacific Standard Time, UTC-08:00)  #    Comments [6]

kick it on DotNetKicks.com
 Wednesday, September 17, 2008

Prime numbers.

If there’s one mathematical curiosity that appears more often than any other in the Project Euler problems, it’s prime numbers. To be fair, we've dealt with primes before, but problem seven is the first that requires a prime number generator as part of its solution.

By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.

What is the 10001st prime number?

I must admit that I’ve been riffing on this problem for quite a while now. There are so many variations to prime number generation that I’ve had difficulty choosing one with the right balance of elegance and efficiency. Because I’ll be reusing my generator for future problems, I must be certain that it’s fast enough. However, as always my primary goal is to produce the most beautiful solution that I can.

Spoiler alert! I'm revealing the problem solution early.

primes |> Seq.nth 10000

Trivial, right? Of course, the challenge of this problem is in declaring the magic behind primes. We have a couple of options available to us, but since there’s a great deal that can be learned from using brute force to solve a problem, we’ll first try the…

The Naïve Approach

The most straightforward way to generate prime numbers is to simply test every number, starting from 2—the first prime—with some primality test. Perhaps something like the code below.

{ 2L..System.Int64.MaxValue } |> Seq.filter isPrime

For the primality test, we can test the candidate prime against every smaller number. If it is evenly divisible by any smaller number other than 1, it isn't prime.

let isPrime n =
  { 2L..n-1L } |> Seq.for_all (fun x -> n%x <> 0L)

Putting the pieces together gives us a real working solution.

let isPrime n =
  { 2L..n-1L } |> Seq.for_all (fun x -> n%x <> 0L)

let primes =
  { 2L..System.Int64.MaxValue } |> Seq.filter isPrime

primes |> Seq.nth 10000

Are we finished? Hardly! While simple, this prime number generator takes a whopping 30 seconds on my hefty-yet-quickly-obsoleting machine! That might fall within Project Euler's "one-minute rule," but we can certainly do better.

Optimizing Naïvely

The obvious optimization is to reduce the set of numbers that isPrime tests. Observe that1 the largest factor of a number (other than itself) is its square root. Armed with this knowledge, we can improve the primality test by testing just the natural numbers from 2 through the square root of n.

let isPrime n =
  let limit = int64 (sqrt (float n))
  { 2L..limit } |> Seq.for_all (fun x -> n%x <> 0L)

That's a huge improvement! Finding the 10001st prime now only takes about .25 seconds.

We can do better yet. If 2 is the only even prime, why are we bothering to test any other even numbers? We can save more time by testing only the odds.

let odds limit =
  Seq.unfold (fun n ->
    if n > limit then None
    else Some(n,n+2L)) 3L

let isPrime n =
  let limit = int64 (sqrt (float n))
  odds limit |> Seq.for_all (fun x -> n%x <> 0L)

let primes =
  seq { yield 2L
        yield! odds System.Int64.MaxValue |> Seq.filter isPrime }

That brings the solution down to approximately .2 seconds.

NeRd Note
Curious about the use of yield and yield! in the example above? The difference between these keywords is simple yet powerful. yield simply returns a single element of a sequence expression, much like the yield return of C# iterators. However, yield! returns another sequence expression as part of the sequence.2 This is an extraordinarily powerful feature that offers a lot of flexibility. In part 2, we'll use yield! to produce an elegant recursive sequence expression.

Before moving on, let's make one last optimization to our naïve algorithm. We can take advantage of the fact that every prime, after 2 and 3, is of the form 6k ± 1. By reducing the set of numbers used by our primality test to those of this form, we can eke out a tiny bit more speed.

let inline next k i =
  if i = -1L then (k,1L)
  else ((k+1L),-1L)

let candidates limit =
  Seq.unfold (fun (k,i) ->
    let v = 6L*k + i
    if (v > limit) then None
    else Some(v, (next k i))) (1L,-1L)

let isPrime n =
  let limit = int64 (sqrt (float n))
  candidates limit |> Seq.for_all (fun x -> n%x <> 0L)

let primes =
  seq { yield! [2L;3L]
        yield! candidates System.Int64.MaxValue |> Seq.filter isPrime }

Using the prime number generator above, our solution takes around .15 seconds. Not too shabby! To be fair, this problem deals with reasonably small prime numbers. Future Project Euler problems (like problem ten) will benefit from a more efficient algorithm. Next time we’ll take a look at another well-known algorithm for generating prime numbers.

 

1"Observe that"? Clearly, I've been reading too many academic papers lately.
2F#'s yield! is similar to the stream-flattening concept of . Perhaps this would be an useful extension to the C# language? Mads? What do you think?

posted on Wednesday, September 17, 2008 1:33:05 PM (Pacific Standard Time, UTC-08:00)  #    Comments [1]

kick it on DotNetKicks.com
 Tuesday, May 06, 2008

Project Euler problem six is another easy one.

The sum of the squares of the first ten natural numbers is,

12 + 22 + ... + 102 = 385

The square of the sum of the first ten natural numbers is,

(1 + 2 + ... + 10)2 = 552 = 3025

Hence the difference between the sum of the squares of the first ten natural numbers and the square of the sum is 3025 − 385 = 2640.

Find the difference between the sum of the squares of the first one hundred natural numbers and the square of the sum.

The solution to this problem boils down to a few folding operations and a map. The one-liner is below.

List.fold_left (+) 0 [1..100] * List.fold_left (+) 0 [1..100] - List.fold_left (+) 0 (List.map (fun x -> x * x) [1..100])

Pretty nasty, eh? Quite a bit of code duplication can be removed. Since they're identical, let's generalize all of the folds first by extracting them to a sum function.

let sum lst = List.fold_left (+) 0 lst

sum [1..100] * sum [1..100] - sum (List.map (fun x -> x * x) [1..100])

That already looks a lot better.

Next, we can generalize the multiplication operations. Each time multiplication occurs in the solution above, it's simply squaring a value. So, we can extract those operations into a square function.

let square x = x * x

square (sum [1..100]) - sum (List.map (fun x -> square x) [1..100])

We can simplify that even further. Because the anonymous function passed to List.map just applies its argument to the square function, we can pass square directly.

square (sum [1..100]) - sum (List.map square [1..100])

Next, let's generalize the call to List.map that produces a list of squares by moving it to a new function, squares.

let squares lst = List.map square lst

square (sum [1..100]) - sum (squares [1..100])

At this point, we have a perfectly acceptable solution. It states the problem almost like natural English: "The square of the sum of 1 to 100 minus the sum of the squares of 1 to 100." So, why are there a few more inches left in this article? Well, I'd like to take this a step further.

Thinking more abstractly, what does our solution do? It computes the difference of two calculations that are based on the same list. We can extract this general process to a new function like so:

let difference f1 f2 lst = f1 lst - f2 lst

difference (fun l -> square (sum l)) (fun l -> sum (squares l)) [1..100]

It turns out that we can simplify these anonymous functions in the same way that we did with the square function earlier. However, because there are two functions involved in each calculation, we must compose the functions together. In F#, there are two operators used to perform function composition: the >> operator, which applies the functions from left to right, and the << operator, which applies the functions from right to left. Obviously, we need the latter.

difference (square << sum) (sum << squares) [1..100]

After using the forward pipe operator to move the list to the front, we're finished.

[1..100] |> difference (square << sum) (sum << squares)

"Take the numbers 1 to 100 and find the difference of the square of the sum and the sum of the squares."

Function composition is beautiful.

posted on Tuesday, May 06, 2008 3:21:26 AM (Pacific Standard Time, UTC-08:00)  #    Comments [3]

kick it on DotNetKicks.com

A few days ago, I presented a solution for Project Euler problem four that I didn't really like. The challenge of problem four is to write a function that determines whether a number is a palindrome, that is, whether it reads the same backward as forward. When presented with that challenge, I took an approach that I feel is a bit of a cop-out: converting the number to a string, reversing the string and comparing the result. This felt somehow wrong because I'm not really solving the problem in a mathematical way. So, I'm declaring a mulligan. Below is a new function which properly performs the math necessary to reverse a base-10 number.

let reverse n =
  let rec loop x res =
    if x = 0 then res
    else loop (x/10) (res*10 + (x%10))

  loop n 0

let isPalindrome n =
  n = reverse n

Our original list comprehension below still works properly with the new isPalindrome function.

[ for x in 100..999
    for y in 100..999
      when isPalindrome(x*y) -> x*y ] |> toLargest

This solution is twice as fast as the original string-based solution. In addition, I'd argue that the tail-recursive loop is at least four times as beautiful. :-)

posted on Tuesday, May 06, 2008 3:21:13 AM (Pacific Standard Time, UTC-08:00)  #    Comments [0]

kick it on DotNetKicks.com
 Monday, May 05, 2008

At first glance, Project Euler problem five looks like a walk in the park:

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

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

Sounds easy! The most straightforward solution is to take the sequence of all natural numbers, filter those that are evenly divisible by 1 through 20, and pop off the first element of the sequence (the head). Something like the code below would do the trick.

{ 1L .. Int64.max_int }
  |> Seq.filter (fun n ->
       [1L .. 20L] |> List.for_all (fun d -> n % d = 0L))
  |> Seq.hd

Unfortunately, that solution, while direct, falls far outside of Project Euler's "one-minute rule." It eventually calculates the correct answer but takes as much as 10 minutes on my machine!

OK, let's take a step back. What exactly is the problem asking us to find? Stating it differently, "What is the least common multiple of all of the numbers from 1 to 20?"

The least common multiple (LCM) of two numbers is the smallest number that is evenly divisible by each. Still not familiar? Think about how fractions are added. The first step in adding fractions is to find the least common denominator, which is simply the LCM of the denominators. For example, given the fractions 1/8 and 1/12, we would find the LCM of 8 and 12. Then, the fractions would be rewritten with the LCM (24) as their denominators. Once this is done, we can easily add the fractions 3/24 and 2/24 to get the answer, 5/24.

So, how should we go about calculating the LCM of two numbers? It turns out that there are many well-known possibilities. One of the most popular methods involves finding the prime factors of both numbers. It goes something like this:

Suppose we wanted to find the least common multiple of 160 and 90. First, we would write out the prime factors of each:

160 = 25 * 51
90 = 21 * 32 * 51

The least common multiple can be computed by multiplying the highest power of each unique factor.

lcm(160,90) = 25 * 32 * 51 = 1440.

Many people have chosen this method when working through problem five, and I was tempted to take this road as well because it would allow us to reuse our primeFactors function from problem three. However, the code would be fairly complex.

let countItems lst =
  let incrCount m i =
    match Map.tryfind i m with
    | Some(c) -> Map.add i (c+1) m
    | None -> Map.add i 1 m

  lst |> List.fold_left incrCount Map.empty |> Map.to_list

let lcm x y =
  let rec updateMap m t =
    let i,c = t
    match Map.tryfind i m with
    | Some(v) when v < c -> Map.add i c m
    | None -> Map.add i c m
    | _ -> m

  let factors =
    [x; y]
    |> List.map primeFactors
    |> List.map countItems
    |> List.fold_left (List.fold_left updateMap) Map.empty

  Map.fold (fun i c res -> res * int64 (float i ** float c)) factors 1L

Personally, I feel a sense of accomplishment at writing all of those folds—particularly the double-fold near the end, that's really cool. :-) However, it's pretty far below my standard for code beauty. If you recall, I'm trying to present the most beautiful solution that I can. So, I'm rejecting this solution, even though it's efficient enough to meet Project Euler's requirements. Admittedly, there's a certain beauty in the list transformations, but there's a much better method.

The least common multiple of two numbers can be calculated quite simply using their greatest common divisor (GCD), or the largest number that divides evenly into both numbers. The GCD can be computed easily with the Euclidean algorithm. Here's how it works:

  1. Start with 2 natural numbers, x and y
  2. If y is equal to zero, the answer is x.
  3. If not, set x to the value of y, and y to the remainder of dividing x by y.
  4. Go back to step 2.

For the more visual among you, a flowchart of the Euclidean algorithm is pictured below.

EuclideanAlgorithm

Once we have the GCD, calculating the LCM is easy. Simply divide x by the GCD of x and y, and multiply the result by y. These two algorithms can be implemented quite beautifully in F#.

let rec gcd x y =
  if y = 0 then x
  else gcd y (x%y)

let lcm x y =
  if x = 0 or y = 0 then 0
  else (x / (gcd x y)) * y

However, the F# libraries already supply a function to calculate the GCD of two numbers. The greatest common denominator also goes by another name, highest common factor (HCF), and there is an HCF function in the Microsoft.FSharp.Math.BigInt module. It's a simple matter to rewrite lcm using BigInt.hcf.

open Microsoft.FSharp.Math

let lcm x y =
  if x = 0I or y = 0I then 0I
  else (x / (BigInt.hcf x y)) * y

With lcm in place, would you believe that our solution looks like this?

[1I..20I] |> List.reduce_left lcm

F# can produce truly beautiful code indeed!

posted on Monday, May 05, 2008 5:12:44 AM (Pacific Standard Time, UTC-08:00)  #    Comments [6]

kick it on DotNetKicks.com
 Friday, May 02, 2008

Yet Another Project Euler Series (YAPES) continues with problem four:

A palindromic number reads the same both ways. The largest palindrome made from the product of two 2-digit numbers is 9009 = 91 × 99.

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

The most straightforward way to determine if a number is palindromic is to convert it to a string and compare that string with its reverse. Sound easy? It is!

One minor snag is the lack of a library function in F# for reversing strings, but that's easily defined like so:

module String =
  let rev (s : string) = new string(s.ToCharArray() |> Array.rev)

With String.rev in place, writing an isPalindrome function is trivial.

let isPalindrome n =
  let text = Int32.to_string n
  text = String.rev text

Using a list comprehension, we can generate all of the palindromes that are products of 3-digit numbers. Once we have this list, producing the result is as simple as passing it to the toLargest function that we defined for Problem Three.

[ for x in 100..999
    for y in 100..999
      when isPalindrome (x*y) -> x*y ] |> toLargest

Short and sweet—my favorite!

posted on Friday, May 02, 2008 4:49:17 AM (Pacific Standard Time, UTC-08:00)  #    Comments [1]

kick it on DotNetKicks.com
 Thursday, May 01, 2008

Project Euler problem three is first of many to deal with prime numbers.

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

What is the largest prime factor of the number 600851475143?

Eventually, we'll need a prime number generator to solve some of the more advanced problems, but this problem can be solved efficiently without one. The number in question is small enough (just 12 digits) that the divide-and-conquer method that many of us learned in elementary school will suffice.

Consider how we might use this process to find the prime factors of 140.

140 Is 140 evenly divisible by 2? Yes! Remember 2 and divide 140 by 2.
2 * 70 Is 70 evenly divisible by 2? Yes! Remember 2 and divide 70 by 2.
2 * 2 * 35 Is 35 evenly divisible by 2? No, how about 3? No. 4? Nope. 5? Yes! Remember 5 and divide 35 by 5.
2 * 2 * 5 * 7 And we're done!

This method isn't rocket science, but it gets the job done. In fact, it's pretty fast for reasonably small numbers. After all, we're not trying to find the factors of RSA-200. :-)

The basic algorithm is pictured as a flowchart below.

PrimeFactorization

The following F# function implements our algorithm.

let primeFactors n =
  let inline isFactor n d = n % d = 0L

  let rec nextFactor n d =
    let x = if d = 2L then 3L else d+2L
    if isFactor n x then x else nextFactor n x

  let rec findFactors n d acc =
    if isFactor n d then
      findFactors (n/d) d (d::acc)
    elif n > d then
      findFactors n (nextFactor n d) acc
    else
      acc

  findFactors n 2L [] |> List.rev

To the uninitiated, that function might look pretty complex. In reality, it's extremely simple, but three other functions are nested inside of it. Let's look at each nested function in turn.

let inline isFactor n d = n % d = 0L

There's nothing tricky about isFactor. It simply abstracts the modulo operation that determines whether n is evenly divisible by d.

let rec nextFactor n d =
  let x = if d = 2L then 3L else d+2L
  if isFactor n x then x else nextFactor n x

nextFactor recursively determines the next value of d to be used in the algorithm. There is a small optimization here: nextFactor only produces odd numbers. Since 2 is the only even prime, why bother checking any other evens?

let rec findFactors n d acc =
  if isFactor n d then
    findFactors (n/d) d (d::acc)
  elif n > d then
    findFactors n (nextFactor n d) acc
  else
    acc

The meat of the algorithm is handled by findFactors. Any factors found are cons'd up with the accumulator variable, acc. Note that both findFactors and nextFactor are written tail-recursively, so they can be optimized by the compiler to conserve stack space.

The real body of primeFactors kicks off the recursion:

findFactors n 2L [] |> List.rev.

The result of findFactors is passed to List.rev to return the prime factors in a more logical order (smallest to largest).

A simple test in the F# Interactive Environment shows that primeFactors works as expected.

> primeFactors 140L;;

val it : int64 list = [2L; 2L; 5L; 7L]

Almost done.

Project Euler Problem Three asks, "What is the largest prime factor of the number 600851475143?" That's just a matter of folding the list of prime factors with the max function (from the F# libraries) to get the answer.

primeFactors 600851475143L |> List.fold1_left max

We can generalize the folding logic above with a new function...

let toLargest l = List.fold1_left max l

...And now we can write the following solution.

primeFactors 600851475143L |> toLargest

That's just lovely.

NeRd Note
Eagle-eyed readers might have noticed that the problem could have been solved several inches ago. If primeFactors didn't reorder its results from smallest to largest, the solution to the problem would be in the head of the result list!
primeFactors 600851475143L |> List.hd
However, that solution has some very real consequences. First of all, primeFactors won't return its results in the most logical order, which limits its reusability. Secondly, the intent of the code isn't stated as clearly. And finally, it's a leaky abstraction because the solution relies upon intimate knowledge of how primeFactors returns its results. If primeFactors were changed later, the solution would be broken!
posted on Thursday, May 01, 2008 6:55:56 AM (Pacific Standard Time, UTC-08:00)  #    Comments [2]

kick it on DotNetKicks.com
 Friday, April 25, 2008

Today, I'm tackling Project Euler problem two in F#:

Each new term in the Fibonacci sequence is generated by adding the previous two terms. By starting with 1 and 2, the first 10 terms will be:

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

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

Like problem one, this is pretty trivial. It's a simple matter of filtering even numbers, taking only the numbers less than some value and folding to get the sum. However, a couple of supporting cast members need to be introduced.

First, I need a way to generate the Fibonacci sequence. There are several ways to calculate Fibonacci numbers (including a clever method that takes advantage of the relationship between Fibonacci numbers and the Golden Ratio). However, I'll heed the heavy hints of the Project Euler problem description above and go with the iterative approach.

My first stab at writing a Fibonacci sequence generator in F# is simply an imperative approach, similar to what I might write in C#, wrapped in an F# sequence expression.

let fibs =
  seq {
        let i = ref 1
        let j = ref 2
        while true do
          let x = !i
          yield x
          do i := !j
          do j := !j + x
      }

In fact, that looks remarkably similar to how I might write a Fibonacci generator as a C# iterator.

static IEnumerable<int> Fibs
{
  get
  {
    int i = 1;
    int j = 2;
    while (true)
    {
      int x = i;
      yield return x;
      i = j;
      j = j + x;
    }
  }
}

The F# and C# Fibonacci generators above are functionally equivalent. The obvious syntactic difference is that the F# version uses reference cells to support mutable variables. Because it uses reference cells, the F# version inherits a bit of operator baggage that might looks strange to C# developers. Most confusing is the unary ! operator, which is used to retrieve the value from a reference cell (i.e., i is a reference cell containing an int, and !i is the int contained within). This will likely look bizarre to many programmers used to C-style syntax where the unary ! operator is used to negate its operand.

NeRd Note
While the C# and F# Fibonacci sequence generators above look essentially the same, they're implemented very differently under the covers. The C# iterator compiles to a class implementing IEnumerable<int> that works like a little state machine. However, the F# sequence expression is compiled as a series of continuations.
let fibs =
  Seq.delay (fun () ->
    let i = ref 1
    let j = ref 2
    Seq.generated
      (fun () -> true)
      (Seq.delay (fun () ->
        let x = !i
        Seq.append
          (Seq.singleton x)
          (Seq.delay (fun () ->
            i := !j
            j := !j + x
            Seq.empty)))))
It's OK if your brain hurts.

I dislike the F# sequence expression approach above for one reason: it seems like a cop-out. It works fine, but it just feels wrong. There has to be a more declarative way to generate the Fibonacci sequence. Fortunately, there is! I can use the Seq.unfold function like so:1

let fibs =
  Seq.unfold (fun (i,j) -> Some(i,(j,i+j))) (1,2)

The generator function passed to Seq.unfold uses a tuple to represent both values needed to iterate the Fibonacci sequence. We can verify that this works properly using the F# Interative Environment.

> fibs |> Seq.take 10;;

val it : int list = [1; 2; 3; 5; 8; 13; 21; 34; 55; 89]

OK. Almost done. I just need a way to take values from the Fibonacci sequence that are less than or equal to four million, and then stop. Effectively, I need something like the LINQ TakeWhile query operator. If I had that, I could use it similarly to the following C#.

foreach (int n in Fibs.TakeWhile(n => n <= 4000000)
  Console.WriteLine(n);

I looked through the F# libraries for a function like TakeWhile but couldn't find one. Either I missed it, or it just isn't there (Don?). Fortunately, this function is trivial to define in F# with a sequence expression. In fact, it's the perfect opportunity to use a sequence expression because the code must interact with an IEnumerator<int>, which has an inherently imperative programming model.

module Seq =
  let
takeWhile f (ie : #seq<'a>) =
    seq { use e = ie.GetEnumerator()
          while e.MoveNext() && f e.Current do
            yield e.Current }

It took a little while to get here, but now I'm ready to solve Project Euler problem two. To restate, we're looking for the sum of the even-valued terms from the Fibonacci sequence that are less than or equal to four million. No problem!

fibs
  |> Seq.filter (fun n -> n % 2 = 0)
  |> Seq.takeWhile (fun n -> n <= 4000000)
  |> Seq.fold (+) 0

As stated last time, I want to present the most beautiful solution that I can. To me, that means the solution should be concise, and it should read like natural language. As before, we can achieve this by abstracting some of the logic into functions whose names better indicate the intent.

let isEven n = n % 2 = 0

let isLessThanOrEqualTo x = (fun n -> n <= x)

let sum s = Seq.fold (+) 0 s

With these functions defined, we can rewrite the solution like so:

fibs
  |> Seq.filter isEven
  |> Seq.takeWhile (isLessThanOrEqualTo 4000000)
  |> sum

The beauty of this code is that it simply states the original problem:

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

F# takes care of the rest.

1For those curious about how Seq.unfold works, check out my Apples and Oranges post. For fun, try generating the Fibonacci sequence in C# using the Unfold function presented in that article.

posted on Friday, April 25, 2008 9:44:24 AM (Pacific Standard Time, UTC-08:00)  #    Comments [2]

kick it on DotNetKicks.com
 Thursday, April 24, 2008

Jason Follas has done it again! He has trumped all of our feeble geeky efforts to firmly establish himself as the reigning "King Nerd." Over the years, Jason repeatedly has proven his worthiness. Here are a few highlights.

  • Built his own Donkey Kong machine. Not a MAME cabinet. A real, full-blown Donkey Kong machine. Sold at auction.
  • Played the 1984 Hitchhiker's Guide to the Galaxy text adventure game on SQL Server 2005 via SQLCLR.
  • Presented a public demo of the SQL Server 2005 XML Data Type by importing World of Warcraft guild data via web services.

Today, I began a modest blog post series of F# solutions to the Project Euler problems. Not to be outdone, Jason has started a similar series. However, Jason has far loftier, or rather, nerdier goals for his series. You see, Jason is solving the problems with LUA, the scripting language beneath World of Warcraft. In fact, he's using WoW as his test-bed.

Don't believe me? Check out the screenshot below.

WoWEuler

The man must be stopped before he goes too far.

posted on Thursday, April 24, 2008 7:17:59 PM (Pacific Standard Time, UTC-08:00)  #    Comments [1]

kick it on DotNetKicks.com

For the past several months, I've been using F# to solve at least two Project Euler problems each week. I find this is a great way to sharpen my math skills and my F# skills simultaneously. If you're looking for a way to flex your programming muscles, you really should check out Project Euler.

Last week at the MVP Summit, my friend, Bill Wagner, pressured me to suggested that I post some of my solutions. Now, there are already plenty of smart people posting F# solutions to the Project Euler problems. That's why I've resisted starting my own series: I'm not certain that I have anything new to say on the topic. However, Bill was very convincing (especially when he mentioned that I would be starting a series to a couple hundred C# MVPs).

So, here's the deal. I will try to present the most beautiful solution that I can. Naturally, beauty is in the eye of the beholder, so you might disagree with me. That's OK. Just make certain to let me know why you disagree so that I can grow as a developer. If anything, this about learning to be a better programmer.

Let's get started.

Problem One

If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of these multiples is 23.

Find the sum of all the multiples of 3 or 5 below 1000.

Obviously, Project Euler starts with some easy problems and moves up to harder fare. The first problem is pretty trivial. It's even reminiscent of the famous FizzBuzz question.

In F#, a list of all multiples of 3 or 5 less than 1000 can be created using a list comprehension. Once we have the list, it can be summed by folding with the addition operator.

[ for x in 1 .. 999 when x % 3 = 0 or x % 5 = 0 -> x ]
  |> List.fold_left (+) 0

The power of F# (and functional programming in general) is its ability to express problems in a more declarative way. This lends the language to mathematical problems very naturally. Looking at the solution above, there are some obvious changes that could make it more succinct. First, the duplicated modulo operation can be abstracted away by declaring a new operator. (Did I mention that F# allows you to declare new operators?) Second, we can extract the folding logic to a new function that better describes its intent.

let inline (/:) x y = x % y = 0 // evenly-divisible by...

let sum list = List.fold_left (+) 0 list

With these defined, we can express the solution more cleanly.

[ for x in 1 .. 999 when x /: 3 or x /: 5 -> x ] |> sum

That's beautiful.

posted on Thursday, April 24, 2008 12:37:46 PM (Pacific Standard Time, UTC-08:00)  #    Comments [3]

kick it on DotNetKicks.com