If you just care about practical sieve sizes then a mutable sieve on an array might be much faster, e.g. https://gist.github.com/noughtmare/48b5755768ddec86ad50a3908e721410.
Otherwise you might even be able to use a fast primality test instead of trial division such as Miller or Baille-PSW primality tests. I recently wrote this Miller test which is the fastest primality test for 64 bit words in Haskell that I know of:
{-# LANGUAGE MagicHash, UnboxedTuples, BangPatterns #-}
import GHC.Exts
import GHC.IO
import System.IO.Unsafe
import GHC.Word
{-# INLINE mulMod #-}
mulMod :: Word# -> Word# -> Word# -> Word#
-- mulMod _ _ _ | trace "mulMod" False = undefined
mulMod a b mod =
case timesWord2# a b of { (# h, l #) ->
case quotRemWord2# h l mod of { (# _, r #) -> r }}
{-# INLINE powMod #-}
powMod :: Word# -> Word# -> Word# -> Word#
-- powMod _ _ _ | trace "powMod" False = undefined
powMod a b m = go a b 1## where
go _ 0## r = r
go a b r =
case case b `and#` 1## of 0## -> r; _ -> mulMod r a m of { r ->
case b `uncheckedShiftRL#` 1# of { b ->
case case b of 0## -> a; _ -> mulMod a a m of { a ->
go a b r }}}
sprp :: Word# -> Word# -> Int#
-- sprp _ _ | trace "sprp" False = undefined
sprp n a = go (n `minusWord#` 1##) 0## where
go :: Word# -> Word# -> Int#
go d s =
case d `and#` 0xff## of { 0## -> go (d `uncheckedShiftRL#` 8#) (s `plusWord#` 8##); _ ->
case case d `and#` 0xf## of 0## -> (# d `uncheckedShiftRL#` 4#, s `plusWord#` 4## #); _ -> (# d, s #) of { (# d, s #) ->
case case d `and#` 0x3## of 0## -> (# d `uncheckedShiftRL#` 2#, s `plusWord#` 2## #); _ -> (# d, s #) of { (# d, s #) ->
case case d `and#` 0x1## of 0## -> (# d `uncheckedShiftRL#` 1#, s `plusWord#` 1## #); _ -> (# d, s #) of { (# d, s #) ->
case powMod a d n of { b ->
case (b `eqWord#` 1##) `orI#` (b `eqWord#` (n `minusWord#` 1##)) of { 1# -> 1#; _ ->
go' b 1## s }}}}}}
go' :: Word# -> Word# -> Word# -> Int#
go' b r s = case r `ltWord#` s of
1# -> case mulMod b b n of { b ->
case b `leWord#` 1## of { 1# -> 0#; _ ->
case b `eqWord#` (n `minusWord#` 1##) of { 1# -> 1#; _ ->
go' b (r `plusWord#` 1##) s }}}
_ -> 0#
{-# NOINLINE isPrime #-}
{-# SCC isPrime #-}
isPrime :: Word# -> Word#
-- isPrime _ | trace "isPrime" False = undefined
isPrime n =
case n `ltWord#` 2## of { 1# -> 0##; _ ->
case n `ltWord#` 4## of { 1# -> n; _ ->
case n `and#` 1## of { 0## -> 0##; _ ->
go 0# }}}
where
ByteArray arr = smallPrimes
go i = case i <# 256# of
1# ->
case indexWord16Array# arr i of { p ->
case n `eqWord#` p of { 1# -> n; _ ->
case n `remWord#` p of 0## -> p; _ -> go (i +# 1#) }}
_ ->
case n `ltWord#` 2647129## of { 1# -> n; _ ->
-- the smallest prime we haven't checked yet is 1627
-- so the smallest composite that passes is its square
case sprp n 2## of { 0# -> 0##; _ ->
case sprp n 3## of { 0# -> 0##; _ ->
case sprp n 5## of { 0# -> 0##; _ ->
case n `ltWord#` 25326001## of { 1# -> n; _ ->
-- This bound was found by:
-- C. Pomerance, J. L. Selfridge and Wagstaff, Jr., S. S.,
-- "The pseudoprimes to 25 · 109,"
-- Math. Comp., 35:151 (1980) 1003-1026. MR 82g:10030
case sprp n 7## of { 0# -> 0##; _ ->
-- 3215031751 already caught by trial division
case n `ltWord#` 118670087467## of { 1# -> n; _ ->
case sprp n 11## of { 0# -> 0##; _ ->
case n `ltWord#` 2152302898747## of { 1# -> n; _ ->
case sprp n 13## of { 0# -> 0##; _ ->
-- I chose to ignore this bound because it is too close to the previous:
-- case n `ltWord#` 3474749660383## of { 1# -> n; _ ->
case sprp n 17## of { 0# -> 0##; _ ->
case n `ltWord#` 341550071728321## of { 1# -> n; _ ->
-- for bases 7 up to 17 these bounds were found by:
-- G. Jaeschke, "On strong pseudoprimes to several bases,"
-- Math. Comp., 61 (1993) 915-926. MR 94d:11004
case sprp n 19## of { 0# -> 0##; _ ->
case sprp n 23## of { 0# -> 0##; _ ->
case n `ltWord#` 3825123056546413051## of { 1# -> n; _ ->
-- This bound was conjectured by:
-- Z. Zhang and M. Tang,
-- "Finding strong pseudoprimes to several bases. II"
-- Math. Comp., 72 (2003) 2085-2097.
case sprp n 29## of { 0# -> 0##; _ ->
case sprp n 31## of { 0# -> 0##; _ ->
case sprp n 37## of { 0# -> 0##; _ ->
n }}}}}}}}}}}}}}}}}}
data ByteArray = ByteArray !ByteArray#
-- all primes under 2^16 except 2
{-# NOINLINE smallPrimes #-}
smallPrimes :: ByteArray
smallPrimes = unsafePerformIO $ IO $ \s ->
case newAlignedPinnedByteArray# (len *# 2#) 64# s of { (# s, arr #) ->
go arr 0# 3# (go0 arr 0# s) }
where
len = 65536#
go0 :: MutableByteArray# a -> Int# -> State# a -> State# a
go0 !arr !i !s
| isTrue# (i >=# len `quotInt#` 4#) = s
| otherwise = go0 arr (i +# 1#) (writeInt64Array# arr i 0# s)
go :: MutableByteArray# a -> Int# -> Int# -> State# a -> (# State# a, ByteArray #)
go !arr !i !j !s
| isTrue# (j ==# len) =
case shrinkMutableByteArray# arr (i *# 2#) s of { s ->
case unsafeFreezeByteArray# arr s of (# s, arr #) -> (# s, {- trace (show (I# i)) $ -} ByteArray arr #) }
| otherwise =
case readInt16Array# arr j s of { (# s, isComposite #) ->
case isComposite of
0# {- | trace ("write" ++ show (I# i, I# j)) True -} ->
case writeWord16Array# arr i (int2Word# j) s of { s ->
case go' arr (j +# j) j s of { s ->
go arr (i +# 1#) (j +# 1#) s }}
_ -> go arr i (j +# 1#) s }
go' :: MutableByteArray# a -> Int# -> Int# -> State# a -> State# a
go' !arr !i !j !s
| isTrue# (i >=# len) = s
| otherwise =
case writeInt16Array# arr i 1# s of s -> go' arr (i +# j) j s
(this is a very low level version; I am planning to make a higher level version)
