A tale of Project Euler

Hi guys.
Somebody just introduced me to a thing called "Project Euler". I gather
it's well known around here...
Anyway, I was a little bit concerned about problem #7. The problem is
basically "figure out what the 10,001st prime number is". Consider the
following two programs for solving this:
First, somebody else wrote this in C:
int n = 2 , m , primesFound = 0;
for( n=0;n < MAX_NUMBERS;n++ )
if( prime[n] )
{
primesFound++;
if( primesFound == 10001 )
cout << n << " is the 10001st prime." << endl;
for(m=2*n;m

Andrew Coppin wrote:
Also, I'm stuck with problem #10. (Find the sum of all primes less than 1 million.) I've let the program run for well over 15 minutes, and still no answer is forthcomming. It's implemented using the same primes function as above, with a simple filter and sum. (The type has to be changed to [Word64] to avoid overflowing the result though.) The other guy claims to have a C solution that takes 12 ms. There's a hell of a difference between 12 ms and over half an hour...(!) Clearly something is horribly wrong here. Uh... help?
I just let it run to completion. Took 25 minutes and 15 seconds to find the (correct) answer. I would have preferred it to take 12 ms...

On Nov 27, 2007 8:44 PM, Andrew Coppin
Andrew Coppin wrote:
Also, I'm stuck with problem #10. (Find the sum of all primes less than 1 million.) I've let the program run for well over 15 minutes, and still no answer is forthcomming. It's implemented using the same primes function as above, with a simple filter and sum. (The type has to be changed to [Word64] to avoid overflowing the result though.) The other guy claims to have a C solution that takes 12 ms. There's a hell of a difference between 12 ms and over half an hour...(!) Clearly something is horribly wrong here. Uh... help?
I just let it run to completion. Took 25 minutes and 15 seconds to find the (correct) answer. I would have preferred it to take 12 ms...
FWIW the following code took 0.77s on my machine: primes :: [Integer] primes = 2 : filter (null . primeFactors) [3,5..] primeFactors :: Integer-> [Integer] primeFactors n = factor n primes where factor m (p:ps) | p*p > m = [] | m `mod` p == 0 = p : factor (m `div` p) (p:ps) | otherwise = factor m ps main = print ( sum ( takeWhile (< 1000000) primes ) ) -- Sebastian Sylvan +44(0)7857-300802 UIN: 44640862

On Nov 27, 2007 2:34 PM, Andrew Coppin
Hi guys.
Somebody just introduced me to a thing called "Project Euler". I gather it's well known around here...
Anyway, I was a little bit concerned about problem #7. The problem is basically "figure out what the 10,001st prime number is". Consider the following two programs for solving this:
First, somebody else wrote this in C:
int n = 2 , m , primesFound = 0;
for( n=0;n < MAX_NUMBERS;n++ ) if( prime[n] ) { primesFound++;
if( primesFound == 10001 ) cout << n << " is the 10001st prime." << endl;
for(m=2*n;m
Second, my implementation in Haskell:
primes :: [Integer] primes = seive [2..] where seive (p:xs) = p : seive (filter (\x -> x `mod` p > 0) xs)
main = print (primes !! 10000)
Well, we know who's winning the beauty contest. But here's the worrying thing:
C program: 0.016 seconds Haskell program: 41 seconds
Eeeps! o_O That's Not Good(tm).
Replacing "primes :: [Integer]" with "primes :: [Word32]" speeds the Haskell version up to 18 seconds, which is much faster but still a joke compared to C.
Running the Haskell code on a 2.2 GHz AMD Athlon64 X2 instead of a 1.5 GHz AMD Athlon drops the execution time down to 3 seconds. (So clearly the box I was using in my lunch break at work is *seriously* slow... presumably the PC running the C code isn't.)
So, now I have a Haskell version that's "only" several hundred times slower. Neither program is especially optimised, yet the C version is drastically faster. This makes me sad. :-(
By the way... I tried using Data.List.Stream instead. This made the program about 40% slower. I also tried both -fasm and -fvia-c. The difference was statistically insignificant. (Hey guys, nice work on the native codegen!) The difference in compilation speed was fairly drastic though, needless to say. (The machine is kinda low on RAM, so trying to call GCC makes it thrash like hell. So does linking, actually...)
Also, I'm stuck with problem #10. (Find the sum of all primes less than 1 million.) I've let the program run for well over 15 minutes, and still no answer is forthcomming. It's implemented using the same primes function as above, with a simple filter and sum. (The type has to be changed to [Word64] to avoid overflowing the result though.) The other guy claims to have a C solution that takes 12 ms. There's a hell of a difference between 12 ms and over half an hour...(!) Clearly something is horribly wrong here. Uh... help?
The algorithm you use to compute primes is actually quite inefficient, and in particular, it is NOT the same algorithm that the C program is using, despite first appearances! The call to 'filter' in the sieve function works by checking *every* number for divisibility by p, and only keeping those that aren't divisible by p; whereas the C program simply marks multiples of p as non-prime, skipping over those which are not multiples. So the Haskell version basically searches for primes by doing trial division on every single integer, whereas the C version is actually using a sieve. For more information you might want to read this paper, which includes a much better primes implementation: www.cs.hmc.edu/~*oneill*/papers/*Sieve*-JFP.pdf In fact, this very same topic was discussed on the list a while back, you can read it here: http://thread.gmane.org/gmane.comp.lang.haskell.cafe/19699/focus=19798 -Brent

Brent Yorgey wrote:
The algorithm you use to compute primes is actually quite inefficient, and in particular, it is NOT the same algorithm that the C program is using, despite first appearances! The call to 'filter' in the sieve function works by checking *every* number for divisibility by p, and only keeping those that aren't divisible by p; whereas the C program simply marks multiples of p as non-prime, skipping over those which are not multiples. So the Haskell version basically searches for primes by doing trial division on every single integer, whereas the C version is actually using a sieve.
I had to reread that several times to figure out what you're actually saying. So the Haskell version tests all N elements to see if P divides them, whereas the C version loops over an array (roughly) N / P times flagging the composite numbers, which is why it's so much damn faster. (Since as P grows, N / P plummets.) OK, cool. So... now how do I implement this without using mutable arrays? :-} PS. The *original* prime number thing I had was much slower. Used an "is_prime" function to do trial division by *every* number below N. (!!) The version I showed as much faster than that, but still quite slow, as you say...

Brent Yorgey wrote:
For more information you might want to read this paper, which includes a much better primes implementation: www.cs.hmc.edu/~oneill /papers/Sieve-JFP.pdf In fact, this very same topic was discussed on the list a while back, you can read it here: http://thread.gmane.org/gmane.comp.lang.haskell.cafe/19699/focus=19798
While not nearly as good as O'Neil's sieve, I think this is a good compromise between beauty and speed: primes = 2 : 3 : sieve (tail primes) [5,7..] where sieve (p:ps) x = let (h, _:t) = span (< p*p) x in h ++ sieve ps (filter ((/=0).(`mod` p)) t) Often in Project Euler problems you need both a factorization function and a list of primes. In that case, I like this: primeFactors = pf primes where pf ps@(p:ps') n | p * p > n = [n] | r == 0 = p : pf ps q | otherwise = pf ps' n where (q, r) = n `quotRem` p isPrime = null . tail . primeFactors primes = 2 : filter isPrime [3,5..] -Yitz

On 11/27/07, Andrew Coppin
Hi guys.
Somebody just introduced me to a thing called "Project Euler". I gather it's well known around here...
Anyway, I was a little bit concerned about problem #7. The problem is basically "figure out what the 10,001st prime number is". Consider the following two programs for solving this:
First, somebody else wrote this in C:
int n = 2 , m , primesFound = 0;
for( n=0;n < MAX_NUMBERS;n++ ) if( prime[n] ) { primesFound++;
if( primesFound == 10001 ) cout << n << " is the 10001st prime." << endl;
for(m=2*n;m
Second, my implementation in Haskell:
primes :: [Integer] primes = seive [2..] where seive (p:xs) = p : seive (filter (\x -> x `mod` p > 0) xs)
main = print (primes !! 10000)
Well, we know who's winning the beauty contest. But here's the worrying thing:
C program: 0.016 seconds Haskell program: 41 seconds
Eeeps! o_O That's Not Good(tm).
Replacing "primes :: [Integer]" with "primes :: [Word32]" speeds the Haskell version up to 18 seconds, which is much faster but still a joke compared to C.
Running the Haskell code on a 2.2 GHz AMD Athlon64 X2 instead of a 1.5 GHz AMD Athlon drops the execution time down to 3 seconds. (So clearly the box I was using in my lunch break at work is *seriously* slow... presumably the PC running the C code isn't.)
So, now I have a Haskell version that's "only" several hundred times slower. Neither program is especially optimised, yet the C version is drastically faster. This makes me sad. :-(
By the way... I tried using Data.List.Stream instead. This made the program about 40% slower. I also tried both -fasm and -fvia-c. The difference was statistically insignificant. (Hey guys, nice work on the native codegen!) The difference in compilation speed was fairly drastic though, needless to say. (The machine is kinda low on RAM, so trying to call GCC makes it thrash like hell. So does linking, actually...)
Also, I'm stuck with problem #10. (Find the sum of all primes less than 1 million.) I've let the program run for well over 15 minutes, and still no answer is forthcomming. It's implemented using the same primes function as above, with a simple filter and sum. (The type has to be changed to [Word64] to avoid overflowing the result though.) The other guy claims to have a C solution that takes 12 ms. There's a hell of a difference between 12 ms and over half an hour...(!) Clearly something is horribly wrong here. Uh... help?
_______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Hi Andrew, I don't remember who I stole this prime computation from, but it is very fast (10001's prime in 0.06 sec with Int and 0.2 with Integer on my machine) and not overly complex: primes :: [Integer] primes = 2 : filter (l1 . primeFactors) [3,5..] where l1 (_:[]) = True l1 _ = False primeFactors :: Integer -> [Integer] primeFactors n = factor n primes where factor _ [] = [] factor m (p:ps) | p*p > m = [m] | m `mod` p == 0 = p : factor (m `div` p) (p:ps) | otherwise = factor m ps I used it a lot while playing with Euler Project. Many of the problems require prime calculation. Cheers, Olivier.

On Nov 27, 2007 8:54 PM, Olivier Boudry
Hi Andrew,
I don't remember who I stole this prime computation from, but it is very fast (10001's prime in 0.06 sec with Int and 0.2 with Integer on my machine) and not overly complex:
primes :: [Integer] primes = 2 : filter (l1 . primeFactors) [3,5..] where l1 (_:[]) = True l1 _ = False
primeFactors :: Integer -> [Integer] primeFactors n = factor n primes where factor _ [] = [] factor m (p:ps) | p*p > m = [m] | m `mod` p == 0 = p : factor (m `div` p) (p:ps) | otherwise = factor m ps
I used it a lot while playing with Euler Project. Many of the problems require prime calculation.
That is indeed a nice and clear version that's pretty fast. It's basically the same as the C version except "backwards" (i.e. examine a number and work backwards through its divisors, rather than filling in a "map" of all multiples of a known prime). Let me suggest the following slight modification (primeFactors in your version doesn't actually return prime factors - it returns prime factors *or* a list of just the number itself), primes :: [Integer] primes = 2 : filter (null . primeFactors) [3,5..] primeFactors :: Integer-> [Integer] primeFactors n = factor n primes where factor m (p:ps) | p*p > m = [] | m `mod` p == 0 = p : factor (m `div` p) (p:ps) | otherwise = factor m ps This is roughly 35% faster on my machine with GHC 6.7.20070730 too, but the point wasn't to make it faster, but clearer. -- Sebastian Sylvan +44(0)7857-300802 UIN: 44640862

On 11/27/07, Sebastian Sylvan
That is indeed a nice and clear version that's pretty fast. It's basically the same as the C version except "backwards" (i.e. examine a number and work backwards through its divisors, rather than filling in a "map" of all multiples of a known prime). Let me suggest the following slight modification (primeFactors in your version doesn't actually return prime factors - it returns prime factors *or* a list of just the number itself),
primes :: [Integer] primes = 2 : filter (null . primeFactors) [3,5..]
primeFactors :: Integer-> [Integer] primeFactors n = factor n primes where factor m (p:ps) | p*p > m = [] | m `mod` p == 0 = p : factor (m `div` p) (p:ps) | otherwise = factor m ps
This is roughly 35% faster on my machine with GHC 6.7.20070730 too, but the point wasn't to make it faster, but clearer. -- Sebastian Sylvan +44(0)7857-300802 UIN: 44640862
Great remark, it's even simpler like this. By the way I just found the article I stole this algorithm from: http://www.haskell.org/haskellwiki/99_questions/31_to_41 last one in problem 35. Cheers, Olivier.

Sebastian Sylvan:
primes :: [Integer] primes = 2 : filter (null . primeFactors) [3,5..]
primeFactors :: Integer-> [Integer] primeFactors n = factor n primes where factor m (p:ps) | p*p > m = [] | m `mod` p == 0 = p : factor (m `div` p) (p:ps) | otherwise = factor m ps
Your definition gives a strange meaning to primeFactors. I'd want that for all n, product (primeFactors n) == n. I think this law holds for the code posted by Olivier. Of course I'd beautify his definition slightly by writing primes = 2 : filter isPrime [3,5..] isPrime = null . drop 1 . primeFactors primeFactors n | n >= 2 = factor primes n factor pr@(p:ps) n | p*p > n = [n] | r == 0 = p : factor pr q | otherwise = factor ps n where (q,r) = quotRem n p Kalman ---------------------------------------------------------------------- Get a free email account with anti spam protection. http://www.bluebottle.com/tag/2

On the other hand, I must relay to you how much fun I had with certain other problems. For example, problem #12. I started with this: triangles = scanl1 (+) [1..] divisors n = length $ filter (\x -> n `mod` x == 0) [1..n] answer = head $ dropWhile (\n -> divisors n < 500) triangles Sadly, this is *absurdly* slow. I gave up after about 5 minutes of waiting. It had only scanned up to T[1,200]. Then I tried the following: triangles = scanl1 (+) [1..] primes :: [Word32] primes = seive [2..] where seive (p:xs) = p : seive (filter (\x -> x `mod` p > 0) xs) factors = work primes where work (p:ps) n = if p > n then [] else if n `mod` p == 0 then p : work (p:ps) (n `div` p) else work ps n count_factors n = let fs = factors n fss = group fs in product $ map ((1+) . length) fss answer = head $ dropWhile (\n -> count_factors n < 500) triangles By looking only at *prime* divisors and then figuring out how many divisors there are in total (you don't even have to *compute* what they are, just *count* them!), I was able to get the correct solution in UNDER ONE SECOND! o_O :-D :^] Now how about that for fast, eh? (Actually, having solved it I discovered there's an even better way - apparently there's a relationship between the number of divisors of consecutive triangular numbers. I'm too noobish to understand the number theory here...) Similarly, problem 24 neatly illustrates everything that is sweet and pure about Haskell: choose [x] = [(x,[])] choose (x:xs) = (x,xs) : map (\(y,ys) -> (y,x:ys)) (choose xs) permutations [x] = [[x]] permutations xs = do (x1,xs1) <- choose xs xs2 <- permutations xs1 return (x1 : xs2) answer = (permutations "0123456789") !! 999999 This finds the answer in less than 3 seconds. And it is beautiful to behold. Yes, there is apparently some more sophisticated algorithm than actually enumerating the permutations. But I love the way I threw code together in the most intuitive way that fell into my head, and got a fast, performant answer. Perfection! :-)

Andrew Coppin wrote:
On the other hand, I must relay to you how much fun I had with certain other problems.
You may want to look at: http://haskell.org/haskellwiki/Euler_problems and make some contributions. But be very careful what you peek at, so don't spoil your own fun. Regards, Yitz

andrewcoppin:
Hi guys.
Somebody just introduced me to a thing called "Project Euler". I gather it's well known around here...
Anyway, I was a little bit concerned about problem #7. The problem is basically "figure out what the 10,001st prime number is". Consider the following two programs for solving this:
First, somebody else wrote this in C:
int n = 2 , m , primesFound = 0;
for( n=0;n < MAX_NUMBERS;n++ ) if( prime[n] ) { primesFound++;
if( primesFound == 10001 ) cout << n << " is the 10001st prime." << endl;
for(m=2*n;m
Second, my implementation in Haskell:
primes :: [Integer] primes = seive [2..] where seive (p:xs) = p : seive (filter (\x -> x `mod` p > 0) xs)
main = print (primes !! 10000)
This is an FAQ. Unless you use the same algorithm and data types in benchmarks, you're not really benchmarking anything. And expecting one of the worst possible algorithms to be good is hoping for a little too much :) When you do use similar data types and algorithms, you get similar performance: http://shootout.alioth.debian.org/gp4/benchmark.php?test=nsievebits&lang=all I'll reiterate: use the same data structures and algorithms, if you want the same performance. -- Don

Don Stewart wrote:
This is an FAQ. Unless you use the same algorithm and data types in benchmarks, you're not really benchmarking anything. And expecting one of the worst possible algorithms to be good is hoping for a little too much :)
Well, if I was comparing my Haskell against some expert-optimised C library then yeah. What puzzled me is that both programs are unoptimised and appear to be using *the same* algorithm, and yet the performance is incomparable. Now that Brent has spoken, I understand the error of my ways. ;-)
When you do use similar data types and algorithms, you get similar performance.
One would hope so, yes... :-)

On Tue, 27 Nov 2007, Andrew Coppin wrote:
So, now I have a Haskell version that's "only" several hundred times slower. Neither program is especially optimised, yet the C version is drastically faster. This makes me sad. :-(
I think the C version is so much faster because it does not need allocations in the inner loop. If you rewrite your program using a mutable array your program should become faster - and more ugly. However it might be fast enough if you create a new array for each run over the array with a new prime. Then you do not need mutable arrays.
participants (8)
-
Andrew Coppin
-
Brent Yorgey
-
Don Stewart
-
Henning Thielemann
-
Kalman Noel
-
Olivier Boudry
-
Sebastian Sylvan
-
Yitzchak Gale