Re: FASTER primes

31 Dec 2009


      On Wed, 2009-12-30 at 11:09 -0500, haskell-cafe-request@haskell.org
wrote:
...
Am Mittwoch 30 Dezember 2009 01:23:32 schrieb Will Ness:
...
Daniel Fischer  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

Steve

tags

participants (1)