{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# OPTIONS_GHC -fspec-constr-count=8 #-}
module Math.NumberTheory.Primes.Sieve.Eratosthenes
( primes
, psieveFrom
, PrimeSieve(..)
, psieveList
, primeList
, primeSieve
, sieveBits
, sieveRange
, sieveTo
) where
import Control.Monad (when)
import Control.Monad.ST
import Data.Bit
import Data.Bits
import Data.Coerce
import Data.Proxy
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Unboxed.Mutable as MU
import Data.Word
import Math.NumberTheory.Primes.Sieve.Indexing
import Math.NumberTheory.Primes.Types
import Math.NumberTheory.Roots
import Math.NumberTheory.Utils.FromIntegral
iXMASK :: Num a => a
iXMASK :: forall a. Num a => a
iXMASK = a
0xFFFFF
iXBITS :: Int
iXBITS :: Int
iXBITS = Int
20
iXJMASK :: Num a => a
iXJMASK :: forall a. Num a => a
iXJMASK = a
0x7FFFFF
iXJBITS :: Int
iXJBITS :: Int
iXJBITS = Int
23
jMASK :: Int
jMASK :: Int
jMASK = Int
7
jBITS :: Int
jBITS :: Int
jBITS = Int
3
sieveBytes :: Int
sieveBytes :: Int
sieveBytes = Int
128 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
1024
sieveBits :: Int
sieveBits :: Int
sieveBits = Int
8 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
sieveBytes
lastIndex :: Int
lastIndex :: Int
lastIndex = Int
sieveBits Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
sieveRange :: Int
sieveRange :: Int
sieveRange = Int
30 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
sieveBytes
data PrimeSieve = PS !Integer {-# UNPACK #-} !(U.Vector Bit)
primeSieve :: Integer -> PrimeSieve
primeSieve :: Integer -> PrimeSieve
primeSieve Integer
bound = Integer -> Vector Bit -> PrimeSieve
PS Integer
0 ((forall s. ST s (Vector Bit)) -> Vector Bit
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s (Vector Bit)) -> Vector Bit)
-> (forall s. ST s (Vector Bit)) -> Vector Bit
forall a b. (a -> b) -> a -> b
$ Integer -> ST s (MVector s Bit)
forall s. Integer -> ST s (MVector s Bit)
sieveTo Integer
bound ST s (MVector s Bit)
-> (MVector s Bit -> ST s (Vector Bit)) -> ST s (Vector Bit)
forall a b. ST s a -> (a -> ST s b) -> ST s b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= MVector s Bit -> ST s (Vector Bit)
MVector (PrimState (ST s)) Bit -> ST s (Vector Bit)
forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
MVector (PrimState m) a -> m (Vector a)
U.unsafeFreeze)
primeList :: forall a. Integral a => PrimeSieve -> [Prime a]
primeList :: forall a. Integral a => PrimeSieve -> [Prime a]
primeList ps :: PrimeSieve
ps@(PS Integer
v Vector Bit
_)
| Proxy a -> Integer -> Bool
forall a. Integral a => Proxy a -> Integer -> Bool
doesNotFit (Proxy a
forall {k} (t :: k). Proxy t
Proxy :: Proxy a) Integer
v
= []
| Integer
v Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 = ([a] -> [Prime a]
forall a b. Coercible a b => a -> b
coerce :: [a] -> [Prime a])
([a] -> [Prime a]) -> [a] -> [Prime a]
forall a b. (a -> b) -> a -> b
$ [a] -> [a]
forall a. Ord a => [a] -> [a]
takeWhileIncreasing ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ a
2 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a
3 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a
5 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: PrimeSieve -> [a]
forall a. Num a => PrimeSieve -> [a]
primeListInternal PrimeSieve
ps
| Bool
otherwise = ([a] -> [Prime a]
forall a b. Coercible a b => a -> b
coerce :: [a] -> [Prime a])
([a] -> [Prime a]) -> [a] -> [Prime a]
forall a b. (a -> b) -> a -> b
$ [a] -> [a]
forall a. Ord a => [a] -> [a]
takeWhileIncreasing ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ PrimeSieve -> [a]
forall a. Num a => PrimeSieve -> [a]
primeListInternal PrimeSieve
ps
primeListInternal :: Num a => PrimeSieve -> [a]
primeListInternal :: forall a. Num a => PrimeSieve -> [a]
primeListInternal (PS Integer
v0 Vector Bit
bs)
= (Int -> a) -> [Int] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map ((a -> a -> a
forall a. Num a => a -> a -> a
+ Integer -> a
forall a. Num a => Integer -> a
fromInteger Integer
v0) (a -> a) -> (Int -> a) -> Int -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> a
forall a. Num a => Int -> a
toPrim)
([Int] -> [a]) -> [Int] -> [a]
forall a b. (a -> b) -> a -> b
$ (Int -> Bool) -> [Int] -> [Int]
forall a. (a -> Bool) -> [a] -> [a]
filter (Bit -> Bool
unBit (Bit -> Bool) -> (Int -> Bit) -> Int -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Bit -> Int -> Bit
forall a. Unbox a => Vector a -> Int -> a
U.unsafeIndex Vector Bit
bs) [Int
lo..Int
hi]
where
(Int
lo, Int
hi) = (Int
0, Vector Bit -> Int
forall a. Unbox a => Vector a -> Int
U.length Vector Bit
bs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
doesNotFit :: forall a. Integral a => Proxy a -> Integer -> Bool
doesNotFit :: forall a. Integral a => Proxy a -> Integer -> Bool
doesNotFit Proxy a
_ Integer
v = a -> Integer
forall a. Integral a => a -> Integer
toInteger (Integer -> a
forall a. Num a => Integer -> a
fromInteger Integer
v :: a) Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer
v
takeWhileIncreasing :: Ord a => [a] -> [a]
takeWhileIncreasing :: forall a. Ord a => [a] -> [a]
takeWhileIncreasing = \case
[] -> []
a
x : [a]
xs -> a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (a -> (a -> [a]) -> a -> [a]) -> (a -> [a]) -> [a] -> a -> [a]
forall a b. (a -> b -> b) -> b -> [a] -> b
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr a -> (a -> [a]) -> a -> [a]
forall a. Ord a => a -> (a -> [a]) -> a -> [a]
go ([a] -> a -> [a]
forall a b. a -> b -> a
const []) [a]
xs a
x
where
go :: Ord a => a -> (a -> [a]) -> a -> [a]
go :: forall a. Ord a => a -> (a -> [a]) -> a -> [a]
go a
y a -> [a]
f a
z = if a
z a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
y then a
y a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a -> [a]
f a
y else []
primes :: Integral a => [Prime a]
primes :: forall a. Integral a => [Prime a]
primes
= ([a] -> [Prime a]
forall {a}. [a] -> [Prime a]
forall a b. Coercible a b => a -> b
coerce :: [a] -> [Prime a])
([a] -> [Prime a]) -> [a] -> [Prime a]
forall a b. (a -> b) -> a -> b
$ [a] -> [a]
forall a. Ord a => [a] -> [a]
takeWhileIncreasing ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ a
2 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a
3 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a
5 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (PrimeSieve -> [a]) -> [PrimeSieve] -> [a]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap PrimeSieve -> [a]
forall a. Num a => PrimeSieve -> [a]
primeListInternal [PrimeSieve]
psieveList
psieveList :: [PrimeSieve]
psieveList :: [PrimeSieve]
psieveList = Integer
-> Integer -> Integer -> Integer -> Vector Word64 -> [PrimeSieve]
makeSieves Integer
plim Integer
sqlim Integer
0 Integer
0 Vector Word64
cache
where
plim :: Integer
plim = Integer
4801
sqlim :: Integer
sqlim = Integer
plimInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
plim
cache :: Vector Word64
cache = (forall s. ST s (Vector Word64)) -> Vector Word64
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s (Vector Word64)) -> Vector Word64)
-> (forall s. ST s (Vector Word64)) -> Vector Word64
forall a b. (a -> b) -> a -> b
$ do
sieve <- Integer -> ST s (MVector s Bit)
forall s. Integer -> ST s (MVector s Bit)
sieveTo (Integer
4801 :: Integer)
new <- MU.unsafeNew 1288 :: ST s (MU.MVector s Word64)
let fill Int
j Int
indx
| Int
1279 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
indx = MVector s Word64 -> m (MVector s Word64)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return MVector s Word64
new
| Bool
otherwise = do
Bit p <- MVector (PrimState m) Bit -> Int -> m Bit
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
MU.unsafeRead MVector s Bit
MVector (PrimState m) Bit
sieve Int
indx
if p
then do
let !i = Int
indx Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
jMASK
k = Int
indx Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftR` Int
jBITS
strt1 = (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
30Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int -> Int
rho Int
i) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int -> Int
byte Int
i) Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftL` Int
jBITS Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int -> Int
idx Int
i
!strt = Int -> Word64
intToWord64 (Int
strt1 Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
forall a. Num a => a
iXMASK)
!skip = Int -> Word64
intToWord64 (Int
strt1 Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftR` Int
iXBITS)
!ixes = Int -> Word64
intToWord64 Int
indx Word64 -> Int -> Word64
forall a. Bits a => a -> Int -> a
`shiftL` Int
iXJBITS Word64 -> Word64 -> Word64
forall a. Num a => a -> a -> a
+ Word64
strt Word64 -> Int -> Word64
forall a. Bits a => a -> Int -> a
`shiftL` Int
jBITS Word64 -> Word64 -> Word64
forall a. Num a => a -> a -> a
+ Int -> Word64
intToWord64 Int
i
MU.unsafeWrite new j skip
MU.unsafeWrite new (j+1) ixes
fill (j+2) (indx+1)
else fill j (indx+1)
vec <- fill 0 0
U.unsafeFreeze vec
makeSieves :: Integer -> Integer -> Integer -> Integer -> U.Vector Word64 -> [PrimeSieve]
makeSieves :: Integer
-> Integer -> Integer -> Integer -> Vector Word64 -> [PrimeSieve]
makeSieves Integer
plim Integer
sqlim Integer
bitOff Integer
valOff Vector Word64
cache
| Integer
valOff' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
sqlim =
let (Vector Word64
nc, Vector Bit
bs) = (forall s. ST s (Vector Word64, Vector Bit))
-> (Vector Word64, Vector Bit)
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s (Vector Word64, Vector Bit))
-> (Vector Word64, Vector Bit))
-> (forall s. ST s (Vector Word64, Vector Bit))
-> (Vector Word64, Vector Bit)
forall a b. (a -> b) -> a -> b
$ do
cch <- Vector Word64 -> ST s (MVector (PrimState (ST s)) Word64)
forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
Vector a -> m (MVector (PrimState m) a)
U.unsafeThaw Vector Word64
cache :: ST s (MU.MVector s Word64)
bs0 <- slice cch
fcch <- U.unsafeFreeze cch
fbs0 <- U.unsafeFreeze bs0
return (fcch, fbs0)
in Integer -> Vector Bit -> PrimeSieve
PS Integer
valOff Vector Bit
bs PrimeSieve -> [PrimeSieve] -> [PrimeSieve]
forall a. a -> [a] -> [a]
: Integer
-> Integer -> Integer -> Integer -> Vector Word64 -> [PrimeSieve]
makeSieves Integer
plim Integer
sqlim Integer
bitOff' Integer
valOff' Vector Word64
nc
| Bool
otherwise =
let plim' :: Integer
plim' = Integer
plim Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
4800
sqlim' :: Integer
sqlim' = Integer
plim' Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
plim'
(Vector Word64
nc,Vector Bit
bs) = (forall s. ST s (Vector Word64, Vector Bit))
-> (Vector Word64, Vector Bit)
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s (Vector Word64, Vector Bit))
-> (Vector Word64, Vector Bit))
-> (forall s. ST s (Vector Word64, Vector Bit))
-> (Vector Word64, Vector Bit)
forall a b. (a -> b) -> a -> b
$ do
cch <- Integer -> Integer -> Vector Word64 -> ST s (MVector s Word64)
forall s.
Integer -> Integer -> Vector Word64 -> ST s (MVector s Word64)
growCache Integer
bitOff Integer
plim Vector Word64
cache
bs0 <- slice cch
fcch <- U.unsafeFreeze cch
fbs0 <- U.unsafeFreeze bs0
return (fcch, fbs0)
in Integer -> Vector Bit -> PrimeSieve
PS Integer
valOff Vector Bit
bs PrimeSieve -> [PrimeSieve] -> [PrimeSieve]
forall a. a -> [a] -> [a]
: Integer
-> Integer -> Integer -> Integer -> Vector Word64 -> [PrimeSieve]
makeSieves Integer
plim' Integer
sqlim' Integer
bitOff' Integer
valOff' Vector Word64
nc
where
valOff' :: Integer
valOff' = Integer
valOff Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
intToInteger Int
sieveRange
bitOff' :: Integer
bitOff' = Integer
bitOff Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
intToInteger Int
sieveBits
slice :: MU.MVector s Word64 -> ST s (MU.MVector s Bit)
slice :: forall s. MVector s Word64 -> ST s (MVector s Bit)
slice MVector s Word64
cache = do
let hi :: Int
hi = MVector s Word64 -> Int
forall a s. Unbox a => MVector s a -> Int
MU.length MVector s Word64
cache Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
sieve <- Int -> Bit -> ST s (MVector (PrimState (ST s)) Bit)
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> a -> m (MVector (PrimState m) a)
MU.replicate (Int
lastIndex Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) (Bool -> Bit
Bit Bool
True)
let treat Int
pr
| Int
hi Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
pr = MVector s Bit -> m (MVector s Bit)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return MVector s Bit
sieve
| Bool
otherwise = do
w <- MVector (PrimState m) Word64 -> Int -> m Word64
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
MU.unsafeRead MVector s Word64
MVector (PrimState m) Word64
cache Int
pr
if w /= 0
then MU.unsafeWrite cache pr (w-1)
else do
ixes <- MU.unsafeRead cache (pr+1)
let !stj = Word64 -> Int
word64ToInt Word64
ixes Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
forall a. Num a => a
iXJMASK
!ixw = Word64 -> Int
word64ToInt (Word64
ixes Word64 -> Int -> Word64
forall a. Bits a => a -> Int -> a
`shiftR` Int
iXJBITS)
!i = Int
ixw Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
jMASK
!k = Int
ixw Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
i
!o = Int
i Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftL` Int
jBITS
!j = Int
stj Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
jMASK
!s = Int
stj Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftR` Int
jBITS
(n, u) <- tick k o j s
let !skip = Int -> Word64
intToWord64 (Int
n Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftR` Int
iXBITS)
!strt = Int -> Word64
intToWord64 (Int
n Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
forall a. Num a => a
iXMASK)
MU.unsafeWrite cache pr skip
MU.unsafeWrite cache (pr+1) ((ixes .&. complement iXJMASK) .|. strt `shiftL` jBITS .|. intToWord64 u)
treat (pr+2)
tick Int
stp Int
off Int
j Int
ix
| Int
lastIndex Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
ix = (Int, Int) -> m (Int, Int)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return (Int
ix Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
sieveBits, Int
j)
| Bool
otherwise = do
Bit p <- MVector (PrimState m) Bit -> Int -> m Bit
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
MU.unsafeRead MVector s Bit
MVector (PrimState m) Bit
sieve Int
ix
when p (MU.unsafeWrite sieve ix (Bit False))
tick stp off ((j+1) .&. jMASK) (ix + stp*delta j + tau (off+j))
treat 0
sieveTo :: Integer -> ST s (MU.MVector s Bit)
sieveTo :: forall s. Integer -> ST s (MVector s Bit)
sieveTo Integer
bound = ST s (MVector s Bit)
ST s (MVector (PrimState (ST s)) Bit)
arr
where
(Int
bytes,Int
lidx) = Integer -> (Int, Int)
forall a. Integral a => a -> (Int, Int)
idxPr Integer
bound
!mxidx :: Int
mxidx = Int
8Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
bytesInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
lidx
mxval :: Integer
mxval :: Integer
mxval = Integer
30Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Int -> Integer
intToInteger Int
bytes Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
intToInteger (Int -> Int
rho Int
lidx)
!mxsve :: Integer
mxsve = Integer -> Integer
forall a. Integral a => a -> a
integerSquareRoot Integer
mxval
(Int
kr,Int
r) = Integer -> (Int, Int)
forall a. Integral a => a -> (Int, Int)
idxPr Integer
mxsve
!svbd :: Int
svbd = Int
8Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
krInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
r
arr :: ST s (MVector (PrimState (ST s)) Bit)
arr = do
ar <- Int -> Bit -> ST s (MVector (PrimState (ST s)) Bit)
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> a -> m (MVector (PrimState m) a)
MU.replicate (Int
mxidx Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) (Bool -> Bit
Bit Bool
True)
let start Int
k Int
i = Int
8Int -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
30Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int -> Int
rho Int
i) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int -> Int
byte Int
i) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int -> Int
idx Int
i
tick Int
stp Int
off Int
j Int
ix
| Int
mxidx Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
ix = () -> m ()
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return ()
| Bool
otherwise = do
Bit p <- MVector (PrimState m) Bit -> Int -> m Bit
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
MU.unsafeRead MVector (PrimState m) Bit
MVector (PrimState (ST s)) Bit
ar Int
ix
when p (MU.unsafeWrite ar ix (Bit False))
tick stp off ((j+1) .&. jMASK) (ix + stp*delta j + tau (off+j))
sift Int
ix
| Int
svbd Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
ix = MVector (PrimState (ST s)) Bit
-> m (MVector (PrimState (ST s)) Bit)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return MVector (PrimState (ST s)) Bit
ar
| Bool
otherwise = do
Bit p <- MVector (PrimState m) Bit -> Int -> m Bit
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
MU.unsafeRead MVector (PrimState m) Bit
MVector (PrimState (ST s)) Bit
ar Int
ix
when p (do let i = Int
ix Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
jMASK
k = Int
ix Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftR` Int
jBITS
!off = Int
i Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftL` Int
jBITS
!stp = Int
ix Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
i
tick stp off i (start k i))
sift (ix+1)
sift 0
growCache :: Integer -> Integer -> U.Vector Word64 -> ST s (MU.MVector s Word64)
growCache :: forall s.
Integer -> Integer -> Vector Word64 -> ST s (MVector s Word64)
growCache Integer
offset Integer
plim Vector Word64
old = do
let num :: Int
num = Vector Word64 -> Int
forall a. Unbox a => Vector a -> Int
U.length Vector Word64
old Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
(Int
bt,Int
ix) = Integer -> (Int, Int)
forall a. Integral a => a -> (Int, Int)
idxPr Integer
plim
!start :: Int
start = Int
8Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
btInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
ixInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
!nlim :: Integer
nlim = Integer
plimInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
4800
sieve <- Integer -> ST s (MVector s Bit)
forall s. Integer -> ST s (MVector s Bit)
sieveTo Integer
nlim
let hi = MVector s Bit -> Int
forall a s. Unbox a => MVector s a -> Int
MU.length MVector s Bit
sieve Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
more <- countFromToWd start hi sieve
new <- MU.unsafeNew (1 + num + 2 * more) :: ST s (MU.MVector s Word64)
let copy Int
i
| Int
num Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
i = () -> m ()
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return ()
| Bool
otherwise = do
MVector (PrimState m) Word64 -> Int -> Word64 -> m ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
MU.unsafeWrite MVector s Word64
MVector (PrimState m) Word64
new Int
i (Vector Word64
old Vector Word64 -> Int -> Word64
forall a. Unbox a => Vector a -> Int -> a
`U.unsafeIndex` Int
i)
Int -> m ()
copy (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
copy 0
let fill Int
j Int
indx
| Int
hi Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
indx = MVector s Word64 -> m (MVector s Word64)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return MVector s Word64
new
| Bool
otherwise = do
Bit p <- MVector (PrimState m) Bit -> Int -> m Bit
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
MU.unsafeRead MVector s Bit
MVector (PrimState m) Bit
sieve Int
indx
if p
then do
let !i = Int
indx Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
jMASK
k :: Integer
k = Int -> Integer
intToInteger (Int
indx Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftR` Int
jBITS)
strt0 = ((Integer
kInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*(Integer
30Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
k Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
intToInteger (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int -> Int
rho Int
i))
Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
intToInteger (Int -> Int
byte Int
i)) Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
`shiftL` Int
jBITS)
Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
intToInteger (Int -> Int
idx Int
i)
strt1 = Integer
strt0 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
offset
!strt = Integer -> Word64
integerToWord64 Integer
strt1 Word64 -> Word64 -> Word64
forall a. Bits a => a -> a -> a
.&. Word64
forall a. Num a => a
iXMASK
!skip = Integer -> Word64
integerToWord64 (Integer
strt1 Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
`shiftR` Int
iXBITS)
!ixes = Int -> Word64
intToWord64 Int
indx Word64 -> Int -> Word64
forall a. Bits a => a -> Int -> a
`shiftL` Int
iXJBITS Word64 -> Word64 -> Word64
forall a. Bits a => a -> a -> a
.|. Word64
strt Word64 -> Int -> Word64
forall a. Bits a => a -> Int -> a
`shiftL` Int
jBITS Word64 -> Word64 -> Word64
forall a. Bits a => a -> a -> a
.|. Int -> Word64
intToWord64 Int
i
MU.unsafeWrite new j skip
MU.unsafeWrite new (j+1) ixes
fill (j+2) (indx+1)
else fill j (indx+1)
fill (num+1) start
{-# INLINE countFromToWd #-}
countFromToWd :: Int -> Int -> MU.MVector s Bit -> ST s Int
countFromToWd :: forall s. Int -> Int -> MVector s Bit -> ST s Int
countFromToWd Int
start Int
end =
(Vector Bit -> Int) -> ST s (Vector Bit) -> ST s Int
forall a b. (a -> b) -> ST s a -> ST s b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap Vector Bit -> Int
countBits (ST s (Vector Bit) -> ST s Int)
-> (MVector s Bit -> ST s (Vector Bit))
-> MVector s Bit
-> ST s Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. MVector s Bit -> ST s (Vector Bit)
MVector (PrimState (ST s)) Bit -> ST s (Vector Bit)
forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
MVector (PrimState m) a -> m (Vector a)
U.unsafeFreeze (MVector s Bit -> ST s (Vector Bit))
-> (MVector s Bit -> MVector s Bit)
-> MVector s Bit
-> ST s (Vector Bit)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Int -> MVector s Bit -> MVector s Bit
forall a s. Unbox a => Int -> Int -> MVector s a -> MVector s a
MU.slice Int
start (Int
end Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
start Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
psieveFrom :: Integer -> [PrimeSieve]
psieveFrom :: Integer -> [PrimeSieve]
psieveFrom Integer
n = Integer
-> Integer -> Integer -> Integer -> Vector Word64 -> [PrimeSieve]
makeSieves Integer
plim Integer
sqlim Integer
bitOff Integer
valOff Vector Word64
cache
where
k0 :: Integer
k0 = ((Integer
n Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
`max` Integer
7) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
7) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`quot` Integer
30
valOff :: Integer
valOff = Integer
30Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
k0
bitOff :: Integer
bitOff = Integer
8Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
k0
start :: Integer
start = Integer
valOffInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
7
ssr :: Integer
ssr = Integer -> Integer
forall a. Integral a => a -> a
integerSquareRoot (Integer
startInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1
end1 :: Integer
end1 = Integer
start Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
6 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
intToInteger Int
sieveRange
plim0 :: Integer
plim0 = Integer -> Integer
forall a. Integral a => a -> a
integerSquareRoot Integer
end1
plim :: Integer
plim = Integer
plim0 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
4801 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- (Integer
plim0 Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
4800)
sqlim :: Integer
sqlim = Integer
plimInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
plim
cache :: Vector Word64
cache = (forall s. ST s (Vector Word64)) -> Vector Word64
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s (Vector Word64)) -> Vector Word64)
-> (forall s. ST s (Vector Word64)) -> Vector Word64
forall a b. (a -> b) -> a -> b
$ do
sieve <- Integer -> ST s (MVector s Bit)
forall s. Integer -> ST s (MVector s Bit)
sieveTo Integer
plim
let (lo,hi) = (0, MU.length sieve - 1)
pct <- countFromToWd lo hi sieve
new <- MU.unsafeNew (2 * pct) :: ST s (MU.MVector s Word64)
let fill Int
j Int
indx
| Int
hi Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
indx = MVector s Word64 -> m (MVector s Word64)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return MVector s Word64
new
| Bool
otherwise = do
Bit isPr <- MVector (PrimState m) Bit -> Int -> m Bit
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
MU.unsafeRead MVector s Bit
MVector (PrimState m) Bit
sieve Int
indx
if isPr
then do
let !i = Int
indx Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
jMASK
!moff = Int
i Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftL` Int
jBITS
k :: Integer
k = Int -> Integer
intToInteger (Int
indx Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftR` Int
jBITS)
p = Integer
30Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
kInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Int -> Integer
intToInteger (Int -> Int
rho Int
i)
q0 = (Integer
startInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`quot` Integer
p
(skp0,q1) = q0 `quotRem` intToInteger sieveRange
(b0,r0)
| q1 == 0 = (-1,6)
| q1 < 7 = (-1,7)
| otherwise = idxPr (integerToInt q1 :: Int)
(b1,r1) | r0 == 7 = (b0+1,0)
| otherwise = (b0,r0+1)
b2 = Integer
skp0Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Int -> Integer
intToInteger Int
sieveBytes Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
intToInteger Int
b1
strt0 = ((Integer
kInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*(Integer
30Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
b2 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
intToInteger (Int -> Int
rho Int
r1))
Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
b2 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
intToInteger (Int -> Int
rho Int
i)
Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
intToInteger (Int -> Int
mu (Int
moff Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
r1))) Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
`shiftL` Int
jBITS)
Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
intToInteger (Int -> Int
nu (Int
moff Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
r1))
strt1 = ((Integer
kInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*(Integer
30Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
k Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
intToInteger (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int -> Int
rho Int
i))
Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
intToInteger (Int -> Int
byte Int
i)) Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
`shiftL` Int
jBITS)
Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
intToInteger (Int -> Int
idx Int
i)
(strt2,r2)
| p < ssr = (strt0 - bitOff,r1)
| otherwise = (strt1 - bitOff, i)
!strt = Integer -> Word64
integerToWord64 Integer
strt2 Word64 -> Word64 -> Word64
forall a. Bits a => a -> a -> a
.&. Word64
forall a. Num a => a
iXMASK
!skip = Integer -> Word64
integerToWord64 (Integer
strt2 Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
`shiftR` Int
iXBITS)
!ixes = Int -> Word64
intToWord64 Int
indx Word64 -> Int -> Word64
forall a. Bits a => a -> Int -> a
`shiftL` Int
iXJBITS Word64 -> Word64 -> Word64
forall a. Bits a => a -> a -> a
.|. Word64
strt Word64 -> Int -> Word64
forall a. Bits a => a -> Int -> a
`shiftL` Int
jBITS Word64 -> Word64 -> Word64
forall a. Bits a => a -> a -> a
.|. Int -> Word64
intToWord64 Int
r2
MU.unsafeWrite new j skip
MU.unsafeWrite new (j+1) ixes
fill (j+2) (indx+1)
else fill j (indx+1)
vec <- fill 0 0
U.unsafeFreeze vec
{-# INLINE delta #-}
delta :: Int -> Int
delta :: Int -> Int
delta = Vector Int -> Int -> Int
forall a. Unbox a => Vector a -> Int -> a
U.unsafeIndex Vector Int
deltas
deltas :: U.Vector Int
deltas :: Vector Int
deltas = [Int] -> Vector Int
forall a. Unbox a => [a] -> Vector a
U.fromList [Int
4,Int
2,Int
4,Int
2,Int
4,Int
6,Int
2,Int
6]
{-# INLINE tau #-}
tau :: Int -> Int
tau :: Int -> Int
tau = Vector Int -> Int -> Int
forall a. Unbox a => Vector a -> Int -> a
U.unsafeIndex Vector Int
taus
taus :: U.Vector Int
taus :: Vector Int
taus = [Int] -> Vector Int
forall a. Unbox a => [a] -> Vector a
U.fromList
[ Int
7, Int
4, Int
7, Int
4, Int
7, Int
12, Int
3, Int
12
, Int
12, Int
6, Int
11, Int
6, Int
12, Int
18, Int
5, Int
18
, Int
14, Int
7, Int
13, Int
7, Int
14, Int
21, Int
7, Int
21
, Int
18, Int
9, Int
19, Int
9, Int
18, Int
27, Int
9, Int
27
, Int
20, Int
10, Int
21, Int
10, Int
20, Int
30, Int
11, Int
30
, Int
25, Int
12, Int
25, Int
12, Int
25, Int
36, Int
13, Int
36
, Int
31, Int
15, Int
31, Int
15, Int
31, Int
47, Int
15, Int
47
, Int
33, Int
17, Int
33, Int
17, Int
33, Int
49, Int
17, Int
49
]
{-# INLINE byte #-}
byte :: Int -> Int
byte :: Int -> Int
byte = Vector Int -> Int -> Int
forall a. Unbox a => Vector a -> Int -> a
U.unsafeIndex Vector Int
startByte
startByte :: U.Vector Int
startByte :: Vector Int
startByte = [Int] -> Vector Int
forall a. Unbox a => [a] -> Vector a
U.fromList [Int
1,Int
3,Int
5,Int
9,Int
11,Int
17,Int
27,Int
31]
{-# INLINE idx #-}
idx :: Int -> Int
idx :: Int -> Int
idx = Vector Int -> Int -> Int
forall a. Unbox a => Vector a -> Int -> a
U.unsafeIndex Vector Int
startIdx
startIdx :: U.Vector Int
startIdx :: Vector Int
startIdx = [Int] -> Vector Int
forall a. Unbox a => [a] -> Vector a
U.fromList [Int
4,Int
7,Int
4,Int
4,Int
7,Int
4,Int
7,Int
7]
{-# INLINE mu #-}
mu :: Int -> Int
mu :: Int -> Int
mu = Vector Int -> Int -> Int
forall a. Unbox a => Vector a -> Int -> a
U.unsafeIndex Vector Int
mArr
{-# INLINE nu #-}
nu :: Int -> Int
nu :: Int -> Int
nu = Vector Int -> Int -> Int
forall a. Unbox a => Vector a -> Int -> a
U.unsafeIndex Vector Int
nArr
mArr :: U.Vector Int
mArr :: Vector Int
mArr = [Int] -> Vector Int
forall a. Unbox a => [a] -> Vector a
U.fromList
[ Int
1, Int
2, Int
2, Int
3, Int
4, Int
5, Int
6, Int
7
, Int
2, Int
3, Int
4, Int
6, Int
6, Int
8, Int
10, Int
11
, Int
2, Int
4, Int
5, Int
7, Int
8, Int
9, Int
12, Int
13
, Int
3, Int
6, Int
7, Int
9, Int
10, Int
12, Int
16, Int
17
, Int
4, Int
6, Int
8, Int
10, Int
11, Int
14, Int
18, Int
19
, Int
5, Int
8, Int
9, Int
12, Int
14, Int
17, Int
22, Int
23
, Int
6, Int
10, Int
12, Int
16, Int
18, Int
22, Int
27, Int
29
, Int
7, Int
11, Int
13, Int
17, Int
19, Int
23, Int
29, Int
31
]
nArr :: U.Vector Int
nArr :: Vector Int
nArr = [Int] -> Vector Int
forall a. Unbox a => [a] -> Vector a
U.fromList
[ Int
4, Int
3, Int
7, Int
6, Int
2, Int
1, Int
5, Int
0
, Int
3, Int
7, Int
5, Int
0, Int
6, Int
2, Int
4, Int
1
, Int
7, Int
5, Int
4, Int
1, Int
0, Int
6, Int
3, Int
2
, Int
6, Int
0, Int
1, Int
4, Int
5, Int
7, Int
2, Int
3
, Int
2, Int
6, Int
0, Int
5, Int
7, Int
3, Int
1, Int
4
, Int
1, Int
2, Int
6, Int
7, Int
3, Int
4, Int
0, Int
5
, Int
5, Int
4, Int
3, Int
2, Int
1, Int
0, Int
7, Int
6
, Int
0, Int
1, Int
2, Int
3, Int
4, Int
5, Int
6, Int
7
]