On Wed, 2009-12-30 at 11:09 -0500, haskell-cafe-requ...@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 _______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe