27 December, 2009

Diving into Erlang with Project Euler, Part 1

Project who?

Project Euler is probably the largest (270 problems as of December 2009) and best collection of brainteasers specifically designed for programmers. Unlike other puzzles they are not only fun and challenging, but also help you discover and develop new patterns and programming techniques. I'm a big fan of similar computational mind-benders and use them both in job interviews to see how a candidate approaches a real problem and as a way to balance the boring work that's frequently necessary for even the most interesting projects.
I've been solving 1-2 problems every day and after completing number 25 today I got a nice message:

"Bravo, dragozov! Now that you have solved 25 problems you have achieved what 80.04% of members have failed to do and have advanced to level 1. Good luck as you continue."

So I thought I'd describe my experience, especially considering that I'm using Erlang, a language that I'm currently trying to master and which could use a little more sample code and demos online.

Why Erlang?

I'm working on a small side project and one of the decisions I made was to focus on learning new things and exploring technologies and problems that I find interesting, but are not useful or relevant to my day job. So I'm developing everything as a combination of Erlang and Python.
Of course I have more pragmatic reasons. Erlang's unique features make lots of sense for some parts of my app and the two languages complement each other quite well. I also wanted to experiment a little more with functional programming and I find the mental shift required to code in Erlang quite enjoyable. I guess the recent hype around the language played its role too.
Anyway I've been studying Erlang and OTP (the open source platform behind it) and use any chance to apply my new skills, so Project Euler was too good of an opportunity to miss.

The Idea

In this and the follow-up articles I'm showing how to use Erlang to solve computational problems. As I advance from the simpler questions to the harder ones, I'll try to also move from the language's basics to more advanced topics.
However this is not an introduction and I won't be explaining the syntax. There are a few good online tutorials and an excellent book, so you might want to check them first.
On the other hand I'll point out Erlang specific features and techniques and discuss their usage. I also keep the solutions as self-contained as possible, so in some cases there might be code implementing functionality that is already available in the standard libraries or functions copied between the modules.
I've tried to find the best solutions I could think of and optimize the code as much as possible, but I haven't searched the Internet or checked the forums, so there might be better/faster algorithms. There is rarely a single correct way of doing things in programming and this is not an exception.
I'll keep all solution sources in a Google Code project called Eulerl. The license is BSD, if you care, so you are free to use them for anything you like.

Level 1 (first 25 problems)

This one is really easy, it just sums two arithmetic progressions, for which the formula is like this:
SumM = 1*M + 2*M + ...+ K*M = M*(1+2+...K)=M*(K*(K+1)/2)
We add the sums of the multiples of 3 and 5 and subtract the sum of the multiples of 15 as they were added twice.
Not much to be seen here, just 2 simple function declarations and the use of the div operator for integer division.

Here we implement a function that iterates through the Fibonacci numbers and sums the even ones.
Although it's still a very simple problem, things get a little more interesting.
First we have an example of a tail-recursion, the most common (and efficient) way of implementing iterations in Erlang. I use the function's arguments, including an accumulator (a very common pattern for keeping track of the result) to keep the recursion's state between subsequent calls. At each step of the iteration the function receives the complete information and returns the whole result.
The code also demonstrates some basic usage of pattern matching and guard expressions. Although they don't make a big difference here, they can make the code for an algorithm shorter and more readable (my personal opinion, of course).

To find the prime factors of a number N iterate from 2 to N/2 checking if N is divisible by the current number. At each iteration collect the factor if testing positive and replace the number tested with the result of the division or, if negative, increment the factor to be tested.
Similarly to problem 3 I use tail-recursion and pattern matching to code the solution. Note the rem operator, which is Erlang's operator for modulo (the remainder of the division of one number by another).

Problem 4 is the first where the solution requires some more coding (at least the solution that I came up with). A palindrome is a number that reads the same both ways. We can generate all palindromes of length 2*N and 2*N - 1 from the numbers with N digits. For example 123321 and 12321 can both be generated from 123.
So to solve the problem I iterate backwards starting with 999 and ending with 100, generate all 6 and 5 digit palindromes and return the first that is a product of two 3 digit numbers.
On the Erlang side the most interesting is the use of functions as first-class elements of the language. You can see them passed as arguments to other functions and assigned to local variables.
Also note the syntax for declaring anonymous (fun) functions and how they have access to the enclosing context (the arguments of the factor/3 function for example). One of the quirks of an anonymous function in Erlang is that it can't call itself directly, so to implement recursion I pass the function as one of its own arguments (Funself is my naming convention, not a rule).
The palindrom_iterator function is interesting for a few reasons.
First it shows a common approach in functional programming. It defines some generic functionality and lets you change the behavior later by accepting a function as one of its arguments. In this case palindrom_iterator iterates through the palindromes and applies the function to each one. Then based on the result it either returns it or continues with the next palindrome.
Second it demonstrates how pattern matching can simplify complex logic. In this case the function iterates through all palindromes bellow a given size (although its not strictly a requirement for the solution). For example if we start with length 6, then we'll iterate through lengths 5, 4, 3 and 2. As it uses numbers half the length to generate the palindromes, there are two cases when it reaches the lowest number of a given length. If it was generating even length numbers, we need to run the same iteration but for odd length palindromes, else we need to reduce the length of the "template" numbers and switch to even palindromes. This is all encoded in only 4 lines at the 2nd and 3rd patterns of the function's declaration.

To find the smallest number that is divided by every number in a list we need to collect the prime factors of all numbers in the list, then multiply them raised by their maximum powers. For example if we have 6, 8 and 18, all prime factors are 2 and 3 (6 = 2*3, 8 = 23 and 18 = 2*(32)) with maximum powers 3 and 2 respectively. So the solution would be (23)*(32)=72.
There are a few new concepts introduced in the Erlang code:
- The use of a look-up table (in the max_count functions) implemented using the ets module, an efficient library for storage and access of data in a table form.
- List comprehensions (in the solution/1 method), Erlang's syntax for "functional" creation of lists.
- The syntax for splitting a list into a "head" element and a "tail" sublist or constructing a list in a similar fashion ([H|T] = List and List = [H|T]).
- The lists module, one of the most used standard libraries.

The square of a sum can be written as:
(a+b+c +c...+y+z)2 = a2+b2+c2...+y2 + z2 + 2*a*b + 2*a*c...+2*a*y+2*a*z + 2*b*c + ...+ 2*b*z + ...2*y*z
(a+b+c+...+y+z)2 - (a2 + b2 + c2 +...+ y2 + z2) = 2*(a*(b+c+...+y+z) + b*(c+...+y+z) + ...+y*z)

The actual implementation is another example of a tail-recursive function, where the state of the computation is kept in the function's arguments.

To solve this I implement a simple check for primality where I keep a list of all prime numbers found until now and test new candidates against the list.

Iterate through the digits, keeping a first-in-last-out list of the last N digits and check if their product is the maximum until now.

The solution is an example of a nested iteration. We have 2 conditions:
(1) A2 + B2 = C2
(2) A + B + C = N
If we postulate that A is the smallest and C the biggest number, we can extract two additional constraints:
(3) A < N/3
(4) B < N/2
I run an iteration through all possible A, for each run another nested iteration through all possible B and check if the resulting triplet {A, B, C} meets all the conditions and constraints (1,2,3 and 4).
The code also demonstrates the use of tuples to return and unpack multiples values as a result of a function.

Similar to problem 7, solved with a simple primality check and a tail-recursive iteration until we've collected the sum of the first N primes.

My solution is to generate all possible subsquares of a given size, calculate all the different products of each subsqare (4 horizontal, 4 vertical and 2 diagonals in this case) and find the maximum one.
The code introduces the array module, providing functions for working with arrays in Erlang (0 based unlike most other data structures). It's also an example of extensive use of anonymous (lambda) functions to simplify the implementation of complex logic.

The N-th triangular number is equal to the sum of a simple arithmetic progression:
(1) Tn = n*(n+1)/2
n and n-1 can be written as products of prime factors:
(2) n = 2p2*3p3...*NpN
(3) n+1 = 2q2*3q3...*MpM
so Tn:
(4) Tn = 2(p2 + q2 - 1)*3(p3 + q3)*...*K(pK+qK), K is the max of M and N and pi,qi >= 0
If the prime factors of x are:
(5) x = 2s2*3s3*...*LsL
then all possible divisors of x have the form:
(6) D = 2t2*3t3*...*LtL, where 0 <= ti <= si
So the total number of unique divisors of x is:
(7) Divisors(x) = (s2 + 1)*(s3+1)*...*(sL+1), adding 1 for the case where ti = 0 (see 6 above)
From 4 and 7:
(8) Divisors(Tn) = (p2 + q2)*(p3+q3+1)*...*(pk+qK+1)
First I implement a method that returns all the prime factors and their powers. Then I iterate through the triangular numbers and calculate expression 8 for each until I finds the first that meets the problem's condition.
The code is another example of a recursive function, with the optimization of reusing (through the function's arguments) the calculated factors and their powers for each number in two consecutive calls, as the (n+1) in one triangular number is the n in the next.

Yet another example of a recursive function and the use of pattern matching to implement a solution. Just sum the digits from right to left, keeping the addition carry in the function's arguments and storing the resulting digits in an "accumulator" list.

My solution is to iterate through all numbers bellow 1 Million and calculate their chain lengths, but using a look-up table to reuse the results in longer chains that contain previously computed shorter ones.

An N x N grid, contains (N+1)x(N+1) nodes. The number of paths from a node to the bottom-right corner is equal to the sum of the number of paths from its right and bottom neighbor nodes (if they exist). As some nodes have to be counted multiple times a look-up table is used (the ets module to the rescue again).

I implement a function that takes a number as a list of its digits and returns a list of the digits of the number multiplied by 2. Then getting the digits of any power of 2 means, starting with 1 (20), repeatedly calling the function with the result of its previous call. Nothing much else to be seen here.

This is an example of using Erlang's pattern matching syntax to implement a set of facts. The solution then uses these facts to iterate through all numbers and calculate the sum.

To solve this problem I'm building a tree-like data structure with Erlang (it's not a real tree as each node except the root has two parents).
Every node is represented by a tuple of four items. First the value of the node itself, then the tuples of the left and right children (or empty tuples if it's the bottom row) and then the "weight" of the node - the sum of the its value with the maximum weight of its children (as a matter of fact we only need to keep the weights). Building the tree from the bottom up gives the solution as the weight of its root.

Another example of using pattern matching to encode the facts of the problem. Then iterate through all the 1-st days of the months of the 20th century's years and check if they were Sundays.

Very similar to problem 16, but this time I implement a method that multiplies an arbitrary number to another and returns the result as a list of its digits.

From the expression for an arbitrary divisor of a number (problem 12, formula 6) we can calculate the sum of all possible divisors:
(1) DivisorsSum(x) = (p20 +p21 +...+p2k2)*(p30 +p31 +...+p3k3)*...(pN0 +pN1 +...+pNkN), where k2,k3, ...kN are the powers of each prime factor of x
Using this we can iterate through all numbers and check if they are amicable and meet the other conditions of the problem, then sum them into the result.

The solution is very straightforward on a technical level.
However the code introduces some important libraries and concepts:
- Working with files is usually done using the file module (and a few others usually prefixed with file_).
- The binary syntax and the pattern matching are a very powerful combination frequently used as a replacement for the standard (and somewhat limited) representation of strings as lists of integers. They make it really easy to parse and construct streams of bytes, while keeping the code short and readable. A personal favorite of mine.
I could have probably used standard libraries for a few of the functions, but chose to implement them myself:
- The qsort function demonstrates a very cool (but not optimal) one-line implementation of the QSort algorithm using list comprehensions, guard expressions and recursion.
- The a2i function demonstrates the syntax for retrieving the integer value of a character by prefixing it with a '$'.
- The split function demonstrates the binary syntax as discussed above.

The sum of all numbers that CAN'T be represented as a sum of 2 abundant ones in a given range is equal to the sum of all numbers in the range minus the sum of the numbers that CAN be represented like this:
(1) SumNA(Max) = Sum(Max) - SumA(Max) = Max*(Max + 1) /2 - SumA(Max)
If we know all abundant numbers in the range their sum can be calculated like this:
(2) SumA(Max) = (2*A1 + 2*A2 + ...+ 2*AN + (A1 + A2) + (A1 + A3) + ...+ (A1+AN) + (A2+A3) +...), where A1,...,AN are all abundant numbers bellow Max
So my solution is to iterate through all numbers in the range, check if they are abundant and update their sum according to formula 2 above. Then the solution is computed using formula 1.
The code is a combination of previous techniques, so nothing much to be explained.

Basic combinatorics. The number of permutations of N unique items is N!. Using this I compute every digit of the I-th permutation noting that the digit at position P changes at every (P-1)! permutations.

Another Fibonacci numbers problem solved with a tail recursive iteration until the problem's conditions are met.

Conclusion of Level 1

This was the first post about my experience as I move through the problems. It's a demonstration of the sequential style of programming in Erlang. In the future I intend to focus on the parallel and distributed concepts of the language, which is what makes it so different. I hope that the problems will get a little harder and less repetitive so that I'll be also able to introduce and experiment a little more with the OTP modules.


http://projecteuler.net/ - the home page of Project Euler.
http://erlang.org/ - the home page of the language.
http://code.google.com/p/eulerl/ - my project at Google Code, where I keep all the sources.

13 December, 2009

Random Subset of Combinations

I was solving an interesting problem today and thought I'd write about it and post some nice code.
Imagine you have a list of fields with a discrete set of possible values each. Example:

one of [John, Paul, Marie, ….]
one of [Turner, Conner, Lennon…]
a date between 01 January and 31 December
one of [USA, Canada, Germany, France...]

The number of possible combinations is equal to the multiplication of the number of possible values for each field and grows very fast.
So now imagine that you want a way to randomly generate a small number of combinations, where there are no repetitions (like 100 random, unique people from the example above) and you want it done in an optimal way.
My solution (I have the source in Python so don't freak out from my attempt at explaining):
1. Map every combination to an integer in the range between 0 and the total number of combinations (so the first possible person has id of 0 and the last one has id of TOTAL - 1).
2. Generating a random combination is equivalent to generating a random id between 0 and TOTAL - 1 and then building the corresponding combination.
3. The mapping is done by representing the combinations as points in a multidimensional space using each field as a dimension. The index of the value of the field in the set of all possible value is the coordinate in this dimension.
4. Lets say that we have N fields each with a list of possible values, L1, …, LN and sizes Size1, …, SizeN.
The number of all combinations in a space is the multiplication of all the sizes:
TOTAL = Size1*Size2…*SizeN
5. If we remove the Nth field we are left with a space with a lower dimension and a total of:
TotalN-1 = Size1*Size2*…SizeN-1
We can call this the Subtotal of N or SubN. With some calculations it's clear that
SubN = SizeN-1*SubN-1
and we set Sub1 = 1
6. Lets name the index of the value of the i-th field of a combination Xi. Then the id of a combination can be calculated like this:
ID = XN*SubN + XN-1*SubN-1 + …+ X1 * Sub1
7. If we have the ID we can calculate X1, …XN like this:
XN = ID/SubN (/ is integer division, i.e we ignore the remainder)
XN-1 = (ID % SubN)/SubN-1 (% is modulo i.e. the remainder of the division)

X1 = ((ID % SubN)%SubN-1)…%Sub2
8. So this gives us a way to quickly generate a random combination, but still doesn't solve the initial problem of getting a whole subset of non-repeating, random combinations.
This is equivalent to generating a list of random, non-repeating numbers in the range 0 - (TOTAL - 1).
9. I wasn't sure if this can be done optimally, but after some reading (Great article on k-permutations) it seemed that there is a very efficient algorithm.
For a subset of size K from N numbers it is possible to do it in O(K) both for space and speed. After I implemented the algorithm I found out that Python has a method that can do this (random.sample) so I'm using it in my code, but my implementation is also there (shuffle_range) if you need to do it in a different language or would like to appreciate how short and beautiful it is ;-)
10. The idea of the algorithm is to take all the numbers from 1 to N (bare with me, it gets optimized later) and shuffle them. Then get the first K numbers and throw the rest.
The shuffling is done by putting the numbers in a big array, iterating through them and for each position generating a random number between 0 and n-1 then swapping the value at the current position with the value at the random position.
The optimization that I mention allows this to be reduced from an O(N) to an O(K) algorithm both for space and speed.
Instead of creating an array of size N we use a hash table with the rule that if the table contains a value for a position this is the value at that position, else the value is equal to the index of the position. Then we only iterate from 0 to K-1 and generate the result using the above rule.

And after this long explanation, which will either put you to sleep or make you have nightmares in the next few nights, here is the source code, which I hope will clarify things a bit:

"""Generation of a random subset of all possible combinations for a list
of list of options.

This code is released in the public domain, meaning it is free for any
2009 Plamen Dragozov
import random

class Combinator(object):
"""Manages a list of lists and generates combinations of the elements of these lists.
def __init__(self, lists):
"""Initialize with the given lists and pre-calculate the subtotals.

Subtotals are the number of combinations of all lists without the last one,
so subtotal3 is equal to Size2*Size1.
Total is the number of all possible combinations.
subtotals = []
total = 1
for l in lists:
total = total * len(l)

self.total = total
self.lists = lists
self.subtotals = subtotals

def combination(self, n):
""" A combination is a list where the element at position 'i' is one of the elements in list 'i'.

All possible combinations (with a count of Total) can be mapped to the first Total integers.
This method returns the n-th combination, calculating all element indices in the lists.
result = []
for i in range(len(self.lists) - 1, -1, -1):#iterate backwards, from the last to the 1st
subtotal = self.subtotals[i]
pos = n / subtotal
n = n % subtotal
return result

def random(self, count):
"""Generate 'count' random combinations without repetition.
ids = random.sample(xrange(self.total), min(count, self.total))
return [self.combination(i) for i in ids]

def shuffle_range(n, k = None):
"""Creates a random sub-permutation of size k for the first n integers.

This is pretty much equivalent to random.sample.
if (not k) or (k > n):
k = n
hashT = {}
for i in xrange(k):
j = random.randint(0, n - 1)
#swap whats at i with what's at j and vice versa
hashT[i], hashT[j] = hashT.get(j, j), hashT.get(i, i)

return [hashT.get(i, i) for i in xrange(k)]