Least common multiples - the scenic route

I've been back to Edinburgh for a family reunion and got to spend a morning with my dad who's a maths teacher by trade — I always thought he might be interested in haskell but never got the time to really look at some stuff in a bit of detail. I'd been pointed at Project Euler and thought it would be an interesting thing to look at together.

We worked through 5. What is the smallest number divisible by each of the numbers 1 to 20?

We started off with factorial of 20, and then the product of the primes in the range. Which meant we needed a sieve. I'd written a sieve in Perl (using HOP::Stream

    sub pf {
        my $s = shift;
        my $head = drop $s;
        my $s2 = filter { $_[0] % $head } $s;
        return node($head, promise { pf($s2) } );
    }

I still don't find this technique (wrapping filtering coderefs around each other in a lazy stream) natural, and it took me a while to translate this haskellish technique back into haskell...

 > makeSieve  n = filter (\i -> i `mod` n /= )
 > primes [] = []
 > primes (x:xs) = x : (primes $ (makeSieve x) xs)

Which is all very nice, but of course this is not the right number:

 > *Main> primes [2..20]
 > [2,3,5,7,11,13,17,19]
 > *Main> product $ primes [2..20]
 > 9699690

The right answer is the product of all the numbers 1..20 except for the ones that are already included in another number. So for example 20*19*18*... but not including 10,5,4, as they're already factors of 20.

We had the idea that we needed to get the factors of each number, so that we could then throw out those whose factors have already been used.

So here's a nice prime factorization technique:

 > factorize n = factorize' n (primes [2..]) []
 > 
 > factorize' n (p:px) results
 >    |  n == 1
 >    =  results
 >    |  n `mod` p == 
 >    =  factorize' (n `quot` p) (p:px) (p:results)
 >    |  otherwise
 >    =  factorize' n px results

Used like this

 > *Main> factorize 18
 > [3,3,2]

But now we want to know the powers of each factor: e.g. in this case, 3^2 and 2^1.

 > powf n = map (\l -> (head l, length l) ) $ groupBy (==) n

It made sense for the result to come back as tuples of (value, power):

 > *Main> powf $ factorize 18
 > [(3,2),(2,1)]

And we can get the whole lot like so:

 > *Main>  map (powf.factorize) [1..20]
 > [[],[(2,1)],[(3,1)],[(2,2)],[(5,1)], ....

But now we just want the maximum power for each one. I tried to get my brain in gear to work out how to do this recursively, then gave up and decided to just sort them...

We sort by value and the descending power. It made sense to use the library function sortBy, which takes a sorting function:

 > cmp (a,a') (b,b') | a == b = compare (b') (a')
 >                   | a <  b    = LT
 >                   | otherwise = GT

(compare is the inbuilt function that compares two Ord values)

 > *Main> sortBy cmp $ concatMap (powf.factorize) [1..20]

Now we want the top power for each value. The way that came to mind was to group by the value:

 > groupBy (\(x,_) (y,_) -> x == y) $ 
 >   sortBy cmp $ concatMap (powf.factorize) [1..20]

Now that we've grouped we just want the maximum ord for each:

 > map (\l -> head l) $
 >   groupBy (\(x,_) (y,_) -> x == y) $ 
 >   sortBy cmp $ concatMap (powf.factorize) [1..20]

Giving us (value,power) pairs like so:

 > [(2,4),(3,2),(5,1),(7,1),(11,1),(13,1),(17,1),(19,1)]

And then of course we want to exponentiate these value:

 > let factors = map (\(n,p) -> n ^ p ) $ 
 >   map (\l -> head l) $
 >   groupBy (\(x,_) (y,_) -> x == y) $ 
 >   sortBy cmp $ concatMap (powf.factorize) [1..20]

Now it's just a matter of

 > *Main> product factors

And Bob's your uncle. What could be simpler? ;-)

I clicked through to the Project Euler discussion forums to see a nice Ruby version:

num = (1..20).inject(1) { |result, n| result.lcm n }

OK, Ruby still seems weird to me, but that's because I haven't learnt it. The point is that this seems rather more compact. As soon as I pointed this to my Dad, he said, “Of course it's a least common multiple problem”. Which of course is more or less the approach that we'd taken, but rather more verbosely than the Ruby version.

So I googled “haskell lcm” to check if there was a function with a similar name somewhere. There was indeed a function lcm... in the Prelude itself.

“Oh,” I said, “maybe this will work...” and typed the following into ghci:

 > foldl1 lcm [1..20]

This of course is equivalent to the above... but 20 times shorter. This is both humbling and rather interesting. I think it was worth going through the entire process to get to this point, as I learnt a lot. Also, now having an (inefficient) prime sieve and a prime factoriser, I can easily solve a number of other questions:

3: Find the largest prime factor of 317584931803.

 > *Main> factorize 317584931803

7: Find the 10001st prime number

 > *Main> show $ last (take 10001 $ primes [2..])

(This one takes quite a long time - about 12 seconds - to complete. Apparently there are much better sieves, but we're still comfortably within the project's “minute rule” so I'm not going to prematurely optimize yet!

Update: bluestorm points out that what I believed to be a Sieve of Eratosthenes actually isn't, citing The Genuine Sieve of Eratosthenes, a very readable paper on the topic.