Monday, 28 April 2008

Project Euler Problem 8

Problem 8 in the the Project Euler series is pure number crunching. The site gives a 1000 digit number, and asks us to find the maximum product of 5 consecutive digits within this number. No number theory required here then.

My idea for solving this problem is to create a sliding window into the big number. The window will be 5 digits wide. We'll position it over the first digit, record the product of all the visible digits, then slide it along to the next position. When we can't slide it any further (when we're 5 digits from the end of the number), then we'll look over all the products we've recorded, and report the maximum.

The first problem in writing the code is to know how to deal with the 1000 digit number. C# cannot deal with such a beast by treating it as a number, and in fact that wouldn't be helpful anyway. We need to look at the digits individually. So what we'll do is supply the number as a string, then convert it into a list of integers, each integer representing one digit. This is what the code looks like:

const string InputDigits = @"73167176531330624919225119674426574742355349194934
                            96983520312774506326239578318016984801869478851843
                            85861560789112949495459501737958331952853208805511
                            12540698747158523863050715693290963295227443043557
                            66896648950445244523161731856403098711121722383113
                            62229893423380308135336276614282806444486645238749
                            30358907296290491560440772390713810515859307960866
                            70172427121883998797908792274921901699720888093776
                            65727333001053367881220235421809751254540594752243
                            52584907711670556013604839586446706324415722155397
                            53697817977846174064955149290862569321978468622482
                            83972241375657056057490261407972968652414535100474
                            82166370484403199890008895243450658541227588666881
                            16427171479924442928230863465674813919123162824586
                            17866458359124566529476545682848912883142607690042
                            24219022671055626321111109370544217506941658960408
                            07198403850962455444362981230987879927244284909188
                            84580156166097919133875499200524063689912560717606
                            05886116467109405077541002256983155200055935729725
                            7163626956188267042825248360082325753042075296345";

var digits = InputDigits
   .ToCharArray()
   .Where(c => char.IsDigit(c))
   .Select(c=> int.Parse(c.ToString()))
   .ToList();

Notice that I've used an "@" quoted string to allow it to span multiple lines, so as to preserve the presentation of the number. To create the List digits we first break up the string into individual characters with ToCharArray(). The Where clause filters out any new line or space characters from the sequence (they're there because of the way I've formatted the string). Then the Select clause converts each character into its corresponding number. Finally the ToList() method causes the sequence of integers to be evaluated and put into a List.

Now that we have our list of digits, the next step is to run our sliding window over it. This will create a list of quintuplets for us - mini-sequences, each containing 5 digits. I'm going to show two ways of doing this: one using my Lazy Trampoline from last time, and a second method that, in this case, evaluates much quicker.

The idea with using the Lazy Trampoline is to process the sequence of digits recursively. Remember that the Lazy Trampoline allows a tail recursive function to return a result after each iteration. Given a sequence, we'll return a sequence consisting of the first 5 elements; then we'll drop the first element of the sequence, and recurse, passing the remainder of the sequence as the parameter. The recursion stops when there are fewer than 5 elements left. Here's the code:

var quintupletsSelector = Trampoline.MakeLazyTrampoline((IEnumerable<int> digitsSequence) =>
{
   // if there are five elements remaining in the sequence
   if (digitsSequence.Skip(5).Any())
   {
       return Trampoline.YieldAndRecurse(digitsSequence.Take(5), digitsSequence.Skip(1));
   }
   else
   {
       return Trampoline.YieldBreak<IEnumerable<int>, IEnumerable<int>>();
   }
});

Executing this statement will cause quintupletsSelector to be initialised as a delegate that takes a sequence of integers, and returns a sequence of quintuplets. This is quite a nice demonstration of how you can recursively process a sequence. The problem is that this implementation is fairly slow (on my machine it takes about 20 seconds to produce all the quintuplets). The reason for that is hidden away behind the implementation of the Skip method. We're using it to create a new sequence that starts from the second member of the first sequence. The problem is that Skip works on sequences, not lists. This means that it can't just jump to the index you give it, it has to enumerate all the preceding items, so that it can ignore them. It may look like we're only asking it to Skip one item each time, but remember that it's working recursively. So we start off working with digitSequence. In the next iteration we're working with digitSequence.Skip(1), in the second digitSequence.Skip(1).Skip(1), and so on. Before any real work can be done, we have to run through the sequence, ignoring any digits that we've already processed. Not very efficient! In a real functional language, we would avoid this by using a data structure that actually allowed us to drop the first item in the list, and just work with the remainder in the next iteration. I might investigate such data structures soon, but I don't want to complicate things just yet.

So, in the shower this morning I thought up a faster way of producing the quintuplets. We'll treat the list as a list, rather than pretending its just a sequence. This will allow us to jump in at a particular index, rather than having to scan through from the beginning of the sequence each time.

I created an extension method RangeFrom that will extract a range from a List, starting from a given index:

public static IEnumerable<T> RangeFrom<T>(this IList<T> list, int startIndex, int count)
{
   for (var i = startIndex; i < startIndex + count; i++)
   {
       yield return list[i];
   }
}

With this extension method, we can now create a query expression to create the quintuplets. This completes almost instantaneously, as you would hope

var quintupletsSelector2 = from index in 0.To(digits.Count - 5)
                          select digits.RangeFrom(index, 5);

The final step is easy (since I've also defined an extension method to calculate the Product of a sequence):

quintupletsSelector2
   .Select(quintuplet => quintuplet.Product())
   .Max()
   .DisplayAndPause();

Or if you want to try out the Lazy Trampoline method:

quintupletsSelector(digits)
   .Select(quintuplet => quintuplet.Product())
   .Max()
   .DisplayAndPause();

Here's the complete sample (you'll also need the Trampoline code from last time)

public class Problem8
{
   public static void Main()
   {
       const string InputDigits = @"73167176531330624919225119674426574742355349194934
                                96983520312774506326239578318016984801869478851843
                                85861560789112949495459501737958331952853208805511
                                12540698747158523863050715693290963295227443043557
                                66896648950445244523161731856403098711121722383113
                                62229893423380308135336276614282806444486645238749
                                30358907296290491560440772390713810515859307960866
                                70172427121883998797908792274921901699720888093776
                                65727333001053367881220235421809751254540594752243
                                52584907711670556013604839586446706324415722155397
                                53697817977846174064955149290862569321978468622482
                                83972241375657056057490261407972968652414535100474
                                82166370484403199890008895243450658541227588666881
                                16427171479924442928230863465674813919123162824586
                                17866458359124566529476545682848912883142607690042
                                24219022671055626321111109370544217506941658960408
                                07198403850962455444362981230987879927244284909188
                                84580156166097919133875499200524063689912560717606
                                05886116467109405077541002256983155200055935729725
                                7163626956188267042825248360082325753042075296345";
       var digits = InputDigits
           .ToCharArray()
           .Where(c => char.IsDigit(c))
           .Select(c => int.Parse(c.ToString()))
           .ToList();

       var quintupletsSelector = Trampoline.MakeLazyTrampoline((IEnumerable<int> digitsSequence) =>
       {
           // if there are five elements remaining in the sequence
           if (digitsSequence.Skip(5).Any())
           {
               return Trampoline.YieldAndRecurse(digitsSequence.Take(5), digitsSequence.Skip(1));
           }
           else
           {
               return Trampoline.YieldBreak<IEnumerable<int>, IEnumerable<int>>();
           }
       });

       var quintupletsSelector2 = from index in 0.To(digits.Count - 5)
                                  select digits.RangeFrom(index, 5);

       // use the Trampoline method:
       quintupletsSelector(digits)
           .Select(quintuplet => quintuplet.Product())
           .Max()
           .DisplayAndPause();

       // use the fast method:
       quintupletsSelector2
           .Select(quintuplet => quintuplet.Product())
           .Max()
           .DisplayAndPause();
   }
}

public static class Extensions
{

   public static IEnumerable<int> To(this int start, int end)
   {
       for (int i = start; i <= end; i++)
       {
           yield return i;
       }
   }

   public static IEnumerable<T> RangeFrom<T>(this IList<T> list, int startIndex, int count)
   {
       for (var i = startIndex; i < startIndex + count; i++)
       {
           yield return list[i];
       }
   }

   public static int Product(this IEnumerable<int> factors)
   {
       int product = 1;

       foreach (var factor in factors)
       {
           product *= factor;
       }

       return product;
   }

   public static void DisplayAndPause(this object result)
   {
       Console.WriteLine(result);
       Console.ReadLine();
   }
}

Friday, 25 April 2008

Commenting

I've not had many comments on my posts: 5 in total so far. That could be because nobody's reading the blog; or maybe I'm not writing anything interesting enough to comment on.

But perhaps it's because commenting is too hard: up till now, Blogger has only been accepting comments from registered users. In the interests of getting some feedback, I'm now allowing annoymous comments. But don't get excited, spammers - I will be moderating everything before it's posted.

Thursday, 24 April 2008

Lazy Trampolining

A couple of epsiodes ago, I showed how you can avoid stack overflow in certain kinds of recursive function using a technique called trampolining. Just mentioning such a word made me feel out of breath, so I went on to develop a technique that I've called lazy trampolining, or lazy recursion.

The idea is to combine the trampoline that I created last time with an Iterator. This allows you to create functions that return a result at each iteration before recursing. Because it uses an Iterator, and Iterators are Lazy, you get to decide how long the recursion should carry on for. I'm still trying to work out how valuable this idea is, but since I've already come up with a couple of practical uses, I figured it would be worth sharing with the world.

The simplest example I've come up with so far is a Fibonacci sequence generator. You may remember that we defined one for Project Euler Problem 2 that used a while loop. Here's an equivalent generator using my Lazy Trampoline:

public static IEnumerable<int> FibonacciTerms(int a, int b)
{
   var recursiveFibonacciFunction = Trampoline.MakeLazyTrampoline((int x, int y) =>
                Trampoline.YieldAndRecurse(x + y, y, x + y)
                );

  return recursiveFibonacciFunction(a, b);
}

You can see that I'm calling the Trampoline.YieldAndRecurse method to indicate what I want to happen in the next iteration. I pass first the result from the current iteration that I want included in the output sequence, then I specify the parameters to be used for the next recursive call.

Now for a slightly longer example. I've taken the algorithm that I had for finding the largest prime factor of a number, and modified it so that it returns all the factors of the number as a sequence. In each iteration of the function, a new factor is found, and this is returned as a PrimeFactor object which contains both the prime number, and its multiplicity (how many times the factor divides into the target number).

        private static IEnumerable<PrimeFactor> PrimeFactors(long prime)
        {
            var function = Trampoline.MakeLazyTrampoline((long previousPrime, long remainder) =>
                {
                    // find next factor of number
                    long nextFactor = (previousPrime + 1).To(remainder)
                        .SkipWhile(x => remainder % x > 0)
                        .FirstOrDefault();

                    if (nextFactor == remainder)
                    {
                        return Trampoline.YieldBreak<long, long, PrimeFactor>(new PrimeFactor { Prime = remainder, Multiplicity = 1 });
                    }
                    else
                    {
                        // find its multiplicity
                        long multiplicity = Enumerable.Range(1, Int32.MaxValue)
                            .TakeWhile(x => remainder % (long)Math.Pow(nextFactor, x) == 0).Last();
                        long quotient = remainder / (long)Math.Pow(nextFactor, multiplicity);

                        PrimeFactor nextPrimeFactor = new PrimeFactor { Prime = nextFactor, Multiplicity = multiplicity };

                        if (quotient == 1)
                        {
                            return Trampoline.YieldBreak<long, long, PrimeFactor>(nextPrimeFactor);
                        }
                        else
                        {
                            return Trampoline.YieldAndRecurse<long, long, PrimeFactor>(nextPrimeFactor, nextFactor, quotient);
                        }
                    }
                }
                );

            return function(1, prime);
        }

Here I use the YieldBreak method: this is used to return one final result, and then terminate the recursion. In the code shown below, I've also included a YieldBreak method that just terminates without returning a final result.

This Lazy Trampoline technique appears to me to be a generalisation of the Unfold operation that we discussed last time. Unfold generates a sequence, but at each iteration you're only given the previous item in the sequence from which to generate the next item. With this technique, you get to decide what parameters to pass through to the next iteration.

This is what the code for the Lazy Trampoline looks like.

public static class Trampoline
{

    public static Func<T1, T2, IEnumerable<TResult>> MakeLazyTrampoline<T1, T2, TResult>(this Func<T1, T2, Bounce<T1, T2, TResult>> function)
    {
        return (T1 param1, T2 param2) => LazyTrampoline(function, param1, param2);
    }

    private static IEnumerable<TResult> LazyTrampoline<T1, T2, TResult>(Func<T1, T2, Bounce<T1, T2, TResult>> function, T1 param1, T2 param2)
    {
        T1 currentParam1 = param1;
        T2 currentParam2 = param2;

        while (true)
        {
            Bounce<T1, T2, TResult> result = function(currentParam1, currentParam2);

            if (result.HasResult)
            {
                yield return result.Result;
            }

            if (!result.Recurse)
            {
                yield break;
            }

            currentParam1 = result.Param1;
            currentParam2 = result.Param2;
        }
    }

    public static Bounce<T1, T2, TResult> YieldAndRecurse<T1, T2, TResult>(TResult result, T1 arg1, T2 arg2)
    {
        return new Bounce<T1, T2, TResult>(arg1, arg2, result);
    }

    public static Bounce<T1, T2, TResult> YieldBreak<T1, T2, TResult>()
    {
        return new Bounce<T1, T2, TResult>();
    }

    public static Bounce<T1, T2, TResult> YieldBreak<T1, T2, TResult>(TResult result)
    {
        return new Bounce<T1, T2, TResult>(result);
    }
}

As you can see, I had to be slightly creative in generating the recursive function in MakeLazyTrampoline. This is because the C# compiler won't allow you to create Iterators as lambda expressions (for very good reasons, no doubt). So I define the Iterator as its own method, then the lambda function that I actually return calls through to the Iterator to generate the sequence. As before, I've only shown the two-parameter version of the code. In the version I've included for download, there's a one-parameter version - you can easily extend it to any number of parameters.

You can download the Trampoline code, and the two examples shown above, here.

In the next episode, I'll show you how I used the Lazy Trampoline to solve Project Euler Problem 8.

Monday, 21 April 2008

Project Euler Problem 7 (and 10)

This post tackles Problem 7 in the Project Euler series.

A couple of episodes ago we introduced the Aggregate method. What this does is work its way through a sequence, summarising it in some way, until finally it produces a single result representing all the other items. I gave examples like finding the sum of a sequence, or computing the average value. "Aggregate" is what we call it if we speak C#: in other languages the same operation goes under the names "reduce", "accumulate", "compress" and perhaps most commonly, "fold".

These all convey the idea of working through a data structure (often a list or sequence, but not always), systematically combining the elements; because we're talking functional languages here, when you use one of these Fold operations you get to specify a function saying how the elements are to be combined.

Today, I'm going to introduce the opposite of the Fold operation, Unfold. (Aside: did anybody else have the thought that these ideas where invented by a gal keeping her mind off the tedium of doing the laundry?)

Where Fold takes a sequence and reduces it to a single element, Unfold starts with a single element, and amplifies it into a sequence. Unfold doesn't exist as such in the .Net class libraries, so I've had to code it up myself. It was a lot of work, as you can see:

    public static IEnumerable<T> Unfold<T>(this T seed, Func<T, T> generator)
    {
        // include seed in the sequence
        yield return seed;

        T current = seed;

        // now continue the sequence
        while (true)
        {
            current = generator(current);
            yield return current;
        }
    }

What I've done is to define an Iterator. For its parameters it takes a seed value and a function that will generate the next item in the sequence given the previous item. Then it's just a case of calling the generator function in a loop, yielding up the result each time, then feeding that result back to the generator function. The OddNumbersGreaterThan method shown below is an example of how you might use this Unfold method.

But what has any of this to do with finding the 10001st Prime, as Problem 7 requests?

Well, there are any ways of calculating the sequence of Prime numbers (the Sieve of Eratosthenes being probably the most famous). But the one that works most naturally (if not most efficiently) in a Functional language involves Unfolding.

What you do is start with a small list of the first few primes - this will be our seed list. Then you use an unfold operation that adds a new prime to the list at each stage. It does that by trial division: starting with the next odd number after the biggest prime in the list, check whether any of the known primes divide into it. If not, you've found your new prime; if it is divisible by any, move on to the next candidate.

Obviously, as the list of known primes gets bigger, finding the next prime in this way will get slower and slower: this isn't the most efficient algorithm, as I said. But we can make it fast enough to be acceptable if we remember again that when we're checking a particular candidate, we don't need to worry about checking divisibility by primes bigger than the square root of the candidate: if it is divisible by a prime bigger than the square root, it must also be divisible by a prime smaller than the square root.

This, then is what the solution looks like:

public class Problem7
{
    public static void Main()
    {
        var firstPrimes = new long[] { 2, 3, 5, 7, 11 };
        //
        var primes = firstPrimes.Unfold(priorPrimes =>
                        priorPrimes
                            .Last()
                            .OddNumbersGreaterThan()
                            .SkipWhile(
                                candidate => priorPrimes
                                  .TakeWhile(prime => prime * prime <= candidate)
                                  .Any(prime => candidate.IsDivisibleBy(prime)))
                            .First()
                            );
        //
        primes.Skip(10000).First().DisplayAndPause();
    }
}

The var keyword, by the way, is not defining a variant, or anything loosely-typed like that. It's just telling the compiler to declare a variable but to figure out what the type should be and substitute that type before it continues the compilation.

First up, we create a list of the first few primes, to get the whole thing started. Then, in line 7, we call the Unfold extension method (I'll show you that code in a minute) to generate the sequence of primes by extending the list firstPrimes. We pass the Unfold method a lambda expression ( lines 8 to 15) which tells it how to use the list of all previously found primes, priorPrimes, to generates the next prime.

And this is how it does it. The Last() method in line 9 finds the last item in this list - the biggest prime; this is taken up by the OddNumbersGreaterThan() method which generates a sequence of odd numbers bigger than that prime. The SkipWhile clause causes any of these odd numbers which are divisible by any of the known primes to be ignored; and finally the First() expression returns the first candidate which is not divisible by any of the previous primes - this is the next prime number: and thus we've extended the sequence by one item.

Just to explain that SkipWhile clause in more detail: we use it to ignore odd numbers from the sequence which are not prime. It does this by running through the sequence of known primes smaller than the Square root of candidate (the TakeWhile clause creates the subsequence for us). The Any method checks whether there are any items in a sequence meeting a condition: in this case, are there any primes which divide into our candidate?

It just remains, in line 18, to work our way through the sequence, ignoring the first 10000 items, finally displaying the 10001st.

As a bonus, now you've seen that, you have everything you need to go and solve Problem 10 as well: summing all primes less than 2 million. All it takes is a slight modification to line 18, but I'll leave that to you.

Thanks to Jacob Carpenter for inspiration for this post.

The complete code is shown below:

public class Problem7
{
    public static void Main()
    {
        var firstPrimes = new long[] { 2, 3, 5, 7, 11 };
        //
        var primes = firstPrimes.Unfold(priorPrimes =>
                        priorPrimes
                            .Last()
                            .OddNumbersGreaterThan()
                            .SkipWhile(
                                candidate => priorPrimes
                                  .TakeWhile(prime => prime * prime <= candidate)
                                  .Any(prime => candidate.IsDivisibleBy(prime)))
                            .First()
                            );
        //
        primes.Skip(10000).First().DisplayAndPause();
    }
}

public static class Extensions
{
    public static bool IsDivisibleBy(this long number, long factor)
    {
        return number % factor == 0;
    }

    public static IEnumerable<long> OddNumbersGreaterThan(this long prime)
    {
        return (prime + 2).Unfold(item => item + 2);
    }

    public static void DisplayAndPause(this object result)
    {
        Console.WriteLine(result);
        Console.ReadLine();
    }
}

public static class Functional
{
    public static IEnumerable<T> Unfold<T>(this IList<T> seedList, Func<IList<T>, T> generator)
    {
        List<T> previousItems = new List<T>(seedList);

        // enumerate all the items in the seed list
        foreach (T item in seedList)
        {
            yield return item;
        }

        // now extend the list
        while (true)
        {
            T newItem = generator(previousItems);
            previousItems.Add(newItem);
            yield return newItem;
        }
    }

    public static IEnumerable<T> Unfold<T>(this T seed, Func<T, T> generator)
    {
        // include seed in the sequence
        yield return seed;

        T current = seed;

        // now continue the sequence
        while (true)
        {
            current = generator(current);
            yield return current;
        }
    }
}