[Haskell-cafe] A tale of Project Euler

[Daniel Fischer]

Am Mittwoch, 28. November 2007 22:31 schrieb Andrew Coppin:
> There are problems for which it's important to know how many times a
> given prime factor occurs. And there are other problems where it is
> merely necessary to know which primes are factors. I would say it's
> useful to have *both* functions available...

Yup
>
>
>
> By the way, I'm now using a new and much uglier prime seive:
>
> size = 1000000 :: Word64
>
> calc_primes :: IO [Word64]
> calc_primes = do
>     grid <- newArray (2,size) True
>     seive 2 grid
>   where
>     seive :: Word64 -> IOUArray Word64 Bool -> IO [Word64]
>     seive p g = do
>       mapM_ (\n -> writeArray g n False) [p, 2*p .. size]
>       mp' <- next (p+1) g
>       case mp' of
>         Nothing -> return [p]
>         Just p' -> do
>           ps <- seive p' g
>           return (p:ps)
>
>     next :: Word64 -> IOUArray Word64 Bool -> IO (Maybe Word64)
>     next p g = do
>       if p == size
>         then return Nothing
>         else do
>           t <- readArray g p
>           if t
>             then return (Just p)
>             else next (p+1) g
>
> Seems moderately fast. At least, I can find the 10,001st prime in under
> 1 second... 

One thing: since You check the array bounds, the system needn't check them 
again, use unsafeWrite and unsafeRead. And use Int for the index, that would 
be MUCH faster.
Another thing: you can stop sieving when p*p > size, another speedup
Third thing: In my experience mapM_ is relatively slow, hand coded loops are 
faster (though much uglier), but YMMV
Fourth thing: sieve only odd numbers
Fifth thing: better use an STUArray, don't drag IO in if it's not necessary.

> The obvious downside of course is that you must know how big
> to make the array at the start of your program, and you must compute the
> entire array before you can use any of it. Both limitations not precent
> in the lazy recursive list implementation...
>
The size of the array can easily be estimated by the prime number theorem,
pi(x) ~ x/log x, where pi(x) is the number of primes not exceeding x, so for 
10,001 primes, find x with x/log x about 10,000, add a bit for safety, a 
bound of 120,000 will do.
If you want the n-th prime, you don't want to use the smaller primes anyway, 
if you need all primes, chances are that prime generation will take 
negligible time in comparison to your main algo anyway.

Keep enjoying Project Euler,
Daniel