Incoherent and disjointed opinionated drivel from somewhere near London

Problem 7

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?

Ah, what a nice, straightforward, unambiguous spec! If only business software specifications were so precise.

Way back in problem 3, I took a bit of a wander off-topic and built a prime generator in .Net using the Sieve of Eratosthenes. Armed with this, problem 7 should be easy, right? The sieve implementation generates an IEnumerable, which is non-indexable (i.e. I can’t just say Primes()[10001]), but I can take the first 10,001 and then ask for the last element, which will be the answer to the problem.

There’s a problem with this, however. The sieve requires an upper bound during initialisation. This means it’s great for solving problems like “generate all the primes less than 10,001″, but not so great at answering questions like “what is the 10,001st prime number?”, since it requires foreknowledge of the upper bound.

To illustrate the problem, I’ll take a wild guess at the upper bound. I’m going to guess that the 10,001st prime number is less than 99,999. What happens?

var sieve = new SieveOfEratosthenes(99999);
var result = sieve.Primes().Take(10001).Last();

This generates an answer of 99,991. If I enter this into the Project Euler website, however, it tells me the answer is wrong. Gah! What went wrong? A simple test reveals the problem:

var sieve = new SieveOfEratosthenes(99999);
var primes = sieve.Primes().Take(10001);
var count = primes.Count();

There’s only 9,592 primes generated! As the docs for Take() state (emphasis mine):

Take<TSource> enumerates source and yields elements until count elements have been yielded or source contains no more elements.

Damn. So, looks like my 99,999 guess was too small – with that as an upper bound, the sieve only finds 9,592 primes, and I need the 10,001st. OK, I’ll bump it up by an order of magnitude:

var sieve = new SieveOfEratosthenes(999999);
var result = sieve.Primes().Take(10001).Last();

This gives me the correct answer. Not exactly a wonderful solution though; the idea of having to guess the upper bound is pretty horrendous, and if this was real code it wouldn’t be particularly maintainable – what if the requirements changed and we had to find the nth prime, which happened to be >99,999? We’d have to guess again. Ugh.

Worse, the sieve algorithm precomputes all the primes up to the specified upper bound, meaning that in the above approach I’ve asked the sieve to generate primes up to 999,999 (all 78,498 of them!) despite only needing 10,001. Not very efficient.

Fortunately, the upper bound can be calculated separately. Where n>8601, as in this case, we can use the following equation:

p(n) < n (loge n + loge loge n - 0.9427)

where p(n) is the nth prime number.

Alternatively, for flexibility in handling n<8601, we can use the less accurate

p(n) < n loge loge n

which works for n>5. We can easily precompute the answers for n<=5, or simply calculate on demand.

The formula can be implemented on the sieve class, with a factory method to help when we want to use it:

public static SieveOfEratosthenes
    CreateSieveWithAtLeastNPrimes(int n)
{
    return new SieveOfEratosthenes((long)
            Math.Ceiling(UpperBoundEstimate(n)));
}
 
private static double UpperBoundEstimate(int n)
{
    return n * Ln(n) + n * (Ln(Ln(n)));
}
 
private static double Ln(double n)
{
    return Math.Log(n, Math.E);
}

This leaves us with an overall solution like so:

var sieve = SieveOfEratosthenes.CreateSieveWithAtLeastNPrimes(10001);
var result = sieve.Primes().Take(10001).Last();

This generates a total 10,018 primes, cutting the wasted effort from almost 70,000 superfluous primes to just 17, and takes around 20ms to execute on my machine. Plenty fast enough, I think.

On to the next Project Euler problem (after a bit of a hiatus)…

Problem 5

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?

In common with many of the other Euler problems, there’s a brute-force way to solve this, and a clean algorithmic way. And in common with my other Euler posts so far, I’ll start with the brute-force way

This problem can be tackled head-on with the following approach: Start from n=1 and increment in a loop. Test each value of n by attempting to divide it by all numbers m from 1 to 20. The first number to pass the test (i.e. n mod m is 0 for all values of m) is the answer.

private static long BruteForceSolver()
{
    long result;
    for (result = 1; !Check(result); ++result)
        ;
    return result;
}

private static bool Check(long result)
{
    for (int i = 1; i <= 20; ++i)
    {
        if (result % i != 0)
            return false;
    }

    return true;
}

This works, but it takes >12 seconds to execute on my PC, so it’s not what you’d call efficient (though it is well within the Euler execution time guidelines).

Some speed gains can be achieved by exploiting the information provided in the question itself. We are told that 2520 is the lowest number evenly divisible by all numbers from 1 to 10. Since the problem space (1 to 20) includes all these numbers, the answer must also be evenly divisible by 2520. This allows much bigger increments each loop – rather than incrementing by 1, why not increment by 2520? And since the answer must be greater than or equal to 2520, why not start the loop there instead of 1? Finally, since we already know that 1 to 10 divide evenly into 2520, each inner loop only needs to check numbers 11 to 20.

That should speed things up a bit:

private long BruteForceSolver()
{
    long result;
    for (result = 2520; !Check(result); result += 2520)
        ;
    return result;
}

private bool Check(long result)
{
    for (int i = 11; i <= 20; ++i)
    {
        if (result % i != 0)
            return false;
    }

    return true;
}

And indeed, on my machine this is now down to 150ms or so. It’s still not a very nice way to tackle the problem, though.

Thinking about it from a different angle yields an altogether smarter approach. Imagine we are looking for the lowest number evenly divisible by the numbers 1 to 2.

[1, 2]

Well that’s easy; since there are only two numbers we just find the lowest common multiple (LCM), which in this case is 2 (since 2 % 2 == 0, and 2 % 1 == 0). If we call this sequence s₁, we can say that LCM(s₁) = 2.

OK, now imagine we are solving the same problem for s₂, which contains the numbers 1 to 3.

[1, 2, 3]

You’ll notice that s₂ contains s₁ in its entirety. LCM(s₂) must therefore be a multiple of LCM(s₁), so we can rewrite s₂ as [LCM(s₁), 3], or [2, 3] (since we know LCM(s₁) = 2). Now we are down to two numbers again, so we can calculate the LCM of 2 and 3, which is 6, so LCM(s₂) = 6.

OK, now we solve the problem for the first 4 numbers (s₃).

[1, 2, 3, 4]

This sequence contains s₂, therefore LCM(s₃) is a multiple of LCM(s₂). We can rewrite s₃ as [LCM(s₂), 4], or [6, 4]. Thus, LCM(s₃) = 12.

This can be repeated as many times as necessary. Generally, we have s_n = [LCM(s_n-1), n+1] where n > 0.

This looks recursive, but a better way to think of it is as an excellent example of a fold. A fold is one of the fundamental tools of functional programming. In fact, it is perhaps the most fundamental, since map, filter etc can be implemented as right folds^[1].

I won’t inflict my pitiful Photoshop skills on anyone by trying to graphically represent a fold – try looking at this Wikipedia article if you want to try and visualise it.

Broadly, the behaviour of a fold is to apply a combining function to elements in a list (or other data structure) and accumulate the results. That’s exactly what we want here – our combining function is LCM, and our accumulating value is the LCM of the whole list. Effectively, for list s₃ above, we have

LCM(s₃) = LCM(LCM(LCM(1,2),3),4) = 12

Note how the result of the innermost LCM (applied to values 1 and 2) becomes a parameter to the next LCM, which in turn becomes a parameter to the outermost LCM which returns the result we want.

By using a fold, we can generalise. In Haskell, the whole problem is a one-liner:

foldl lcm 1 [1..20]

The 1 passed in as a parameter represents the terminating value to use when the end of the list is reached. It is common for this value to be the first element of the list, so Haskell provides a convenience function that removes the need to specify it as a parameter:

foldl1 lcm [1..20]

Not all languages and platforms provide an LCM function right out of the box, so to take this neat Haskell solution and port it to .Net, the LCM function needs to be implemented. This is easily done in terms of the greatest common divisor (GCD) like so:

.Net doesn’t provide a GCD function either, so I’ll implement it using Euclid’s Algorithm as an extension method on long ints:

public static long GCD(this long a, long b)
{
    while (b != 0)
    {
        long tmp = b;
        b = a % b;
        a = tmp;
    }

    return a;
}

With GCD defined, LCM can be implemented as above:

public static long LCM(this long a, long b)
{
    return (a * b) / a.GCD(b);
}

With this in place, it’s a simple matter to use .Net’s equivalent of fold – a method on IEnumerable called Aggregate – to get the answer^[2]:

return LongEnumerable.Range(1, 20)
    .Aggregate(1L, (curr, next) => curr.LCM(next));

And indeed, the same basic pattern can be used to solve the problem in a number of languages. In F#, given implementations of LCM and GCD as above, we have:

List.fold_left lcm 1 [1..20]

And in ruby:

require 'rational'
(1..20).inject { |c, n| c.lcm n }

Given that the right algorithm makes this problem a fairly trivial expression in all these languages, it’s pretty hard to identify which is the nicest. I think overall I’ll give the nod to Haskell, however, for not making me implement LCM and because I find ruby’s ‘inject’ a less intuitive function name than foldr (but that’s probably because I learned the technique in Haskell in the first place and am set in my ways…)

let filter p lst =
  List.fold_right (fun x xs -> if p x then x::xs else xs) lst []

let map f lst =
  List.fold_right (fun x xs -> f x :: xs) lst []

^[2]Note that in this code LongEnumerable is just a very simple partial reimplementation of Enumerable, using longs instead of ints