[Haskell-cafe] Re: FASTER primes

[Steve]

On Wed, 2009-12-30 at 11:09 -0500, haskell-cafe-request at haskell.org
wrote:
> Am Mittwoch 30 Dezember 2009 01:23:32 schrieb Will Ness:
> > Daniel Fischer <daniel.is.fischer <at> web.de> writes:

....

> No, it's my own code. Nothing elaborate, just sieving numbers 6k1,
> twice as fast as the 
> haskellwiki code (here) and uses only 1/3 the memory. For the record:
> 
> ----------------------------------------------------------------------
> {-# LANGUAGE BangPatterns #-}
> module SievePrimes where
> 
> import Data.Array.Base (unsafeRead, unsafeWrite, unsafeAt)
> import Data.Array.ST
> import Data.Array.Unboxed
> import Control.Monad (when)
> import Data.Bits
> 
> primesUpTo :: Integer -> [Integer]
> primesUpTo bound
>     = 2:3:[fromIntegral (3*i + (if testBit i 0 then 4 else 5))
>             | i <- [0 .. bd], par `unsafeAt` i]
>       where
>         bd0 = 2*(bound `div` 6) - case bound `mod` 6 of { 0 -> 2; 5 ->
> 0; _ -> 1 }
>         bd = fromInteger bd0
>         sr = floor (sqrt (fromIntegral bound))
>         sbd = 2*(sr `div` 6) - case sr `mod` 6 of { 0 -> 2; 5 -> 0; _
> -> 1 }
>         par :: UArray Int Bool
>         par = runSTUArray $ do
>                 ar <- newArray (0,bd) True
>                 let tick i s1 s2
>                         | bd < i    = return ()
>                         | otherwise = do
>                             p <- unsafeRead ar i
>                             when p (unsafeWrite ar i False)
>                             tick (i+s1) s2 s1
>                     sift i
>                         | i > sbd   = return ar
>                         | otherwise = do
>                             p <- unsafeRead ar i
>                             when p (let (!start,!s1,!s2) = if even i
> then 
> (i*(3*i+10)+7,2*i+3,4*i+7) else (i*(3*i+8)+4,4*i+5,2*i+3) in tick
> start s1 s2)
>                             sift (i+1)
>                 sift 0
> ----------------------------------------------------------------------
> 
> The index magic is admittedly not obvious, but if you take that on
> faith, the rest is 
> pretty clear. To find the n-th prime,
> 
> ----------------------------------------------------------------------
> module Main (main) where
> 
> import SievePrimes (primesUpTo)
> import System.Environment (getArgs)
> 
> printNthPrime :: Int -> IO ()
> printNthPrime n = print (n, primesUpTo (estim n) !! (n-1))
> 
> main = do
>     args <- getArgs
>     printNthPrime $ read $ head args
> 
> estim :: Int -> Integer
> estim k
>     | n < 13    = 37
>     | n < 70    = 349
>     | otherwise = round x
>       where
>         n = fromIntegral k :: Double
>         x = n * (log n + log (log n) - 1 + 1.05 * log (log n) / log n)
> ----------------------------------------------------------------------
> 
> If I don't construct the prime list but count directly on the array,
> it's ~6% faster.

The speed difference varies with the value of 'bound' - the difference
between counting directly on the array and creating a list from the
array increases as 'bound' increases. I tested it for 10^8 and 10^9
(using Int not Integer) and using the array without creating a list is
50+% faster.
It looks like a very fast method for creating the 'set' of primes up to
bound.

Steve