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.

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

|> 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.

Welcome back F# junkies! In this installment of YAPES^{1}, we’re tackling Project Euler problem nine.

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

a^{2}+b^{2}=c^{2}

For example, 3^{2}+ 4^{2}= 9 + 16 = 25 = 5^{2}.

There exists exactly one Pythagorean triplet for whicha+b+c= 1000. Find the productabc.

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

- Generate a stream of Pythagorean triples.
- Find the Pythagorean triple whose sum is equal to 1000.
- 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 integersmandnwithm>n

a= 2mn,b=m^{2}–n^{2},c=m^{2}+n^{2}

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**

*every*Pythagorean triple, we would have to introduce an additional parameter,

*k*, into the formula.

a= 2kmn,b=k(m^{2}–n^{2}),c=k(m^{2}+n^{2})

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

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:

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.

|> 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 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.

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

|> product

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

^{1}Yet Another Project Euler Series

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:

- Transform the 1000-digit number into a list of digits.
- Transform the list of digits into a list of quintuplets—that is, a list containing tuples of each five consecutive elements.
- Transform the list of quintuplets into a list containing the product of each.
- 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:

- Transform the string into a list of chars.
- 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.
- 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.

/// 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.

"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.

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**

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)

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

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

|> toQuintuplets

|> List.map product

|> List.max

It’s that simple!

^{1}Or pirates.

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 YAPES^{1} 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:

- Start with the natural numbers from 2 up to some maximum limit.
- 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.
- Take the next unmarked number and mark all of its multiples.
- Repeat step 3 until there are no more numbers left to test.
- 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 Eratosthenes^{2}, 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.

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.

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.

Let's skip ahead to some more interesting cases.

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

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.

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

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.

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 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:

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 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:

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!

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

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 6^{th}prime is 13.

What is the 10001^{st}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.

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.

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.

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

Putting the pieces together gives us a real working solution.

{ 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 that^{1} 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 limit = int64 (sqrt (float n))

{ 2L..limit } |> Seq.for_all (fun x -> n%x <> 0L)

That's a huge improvement! Finding the 10001^{st} 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.

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**

*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.

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. ^{2}F#'s yield! is similar to the stream-flattening concept of Cω. Perhaps this would be an useful extension to the C# language? Mads? What do you think?

Project Euler problem six is another easy one.

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

1^{2}+ 2^{2}+ ... + 10^{2}= 385

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

(1 + 2 + ... + 10)^{2}= 55^{2}= 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.

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.

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.

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.

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

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:

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.

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

"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.

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 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 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.

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.

|> 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:

160 = 2

^{5}* 5

^{1}

90 = 2

^{1}* 3

^{2}* 5

^{1}

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

lcm(160,90) = 2

^{5}* 3

^{2}* 5

^{1}= 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 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:

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

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

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#.

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.

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?

F# can produce truly beautiful code indeed!

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:

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

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

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 y in 100..999

when isPalindrome (x*y) -> x*y ] |> toLargest

Short and sweet—my favorite!

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.

The following F# function implements our algorithm.

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.

There's nothing tricky about isFactor. It simply abstracts the modulo operation that determines whether n is evenly divisible by 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?

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:

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.

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.

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

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

That's just lovely.

**NeRd Note**

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.

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.

{

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**

*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.

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)))))

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}

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.

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#.

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.

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!

|> 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 isLessThanOrEqualTo x = (fun n -> n <= x)

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

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

|> 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.

^{1}For 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.

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.

The man must be stopped before he goes too far.

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.

|> 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 sum list = List.fold_left (+) 0 list

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

That's beautiful.