-- |
-- Module:      Math.NumberTheory.Primes.Sieve.Eratosthenes
-- Copyright:   (c) 2011 Daniel Fischer
-- Licence:     MIT
-- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>
--
-- Sieve
--
{-# 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

-- Sieve in 128K chunks.
-- Large enough to get something done per chunk
-- and hopefully small enough to fit in the cache.
sieveBytes :: Int
sieveBytes :: Int
sieveBytes = Int
128 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
1024

-- Number of bits per chunk.
sieveBits :: Int
sieveBits :: Int
sieveBits = Int
8 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
sieveBytes

-- Last index of chunk.
lastIndex :: Int
lastIndex :: Int
lastIndex = Int
sieveBits Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1

-- Range of a chunk.
sieveRange :: Int
sieveRange :: Int
sieveRange = Int
30 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
sieveBytes

-- | Compact store of primality flags.
data PrimeSieve = PS !Integer {-# UNPACK #-} !(U.Vector Bit)

-- | Sieve primes up to (and including) a bound (or 7, if bound is smaller).
--   For small enough bounds, this is more efficient than
--   using the segmented sieve.
--
--   Since arrays are 'Int'-indexed, overflow occurs when the sieve
--   size comes near @'maxBound' :: 'Int'@, that corresponds to an
--   upper bound near @15/8*'maxBound'@. On @32@-bit systems, that
--   is often within memory limits, so don't give bounds larger than
--   @8*10^9@ there.
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)

-- | Generate a list of primes for consumption from a
--   'PrimeSieve'.
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
              = [] -- has an overflow already happened?
  | 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)

-- | Returns true if integer is beyond representation range of type a.
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

-- | Extracts the longest strictly increasing prefix of the list
-- (possibly infinite).
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 []

-- | Ascending list of primes.
--
-- >>> take 10 primes
-- [Prime 2,Prime 3,Prime 5,Prime 7,Prime 11,Prime 13,Prime 17,Prime 19,Prime 23,Prime 29]
--
-- 'primes' is a polymorphic list, so the results of computations are not retained in memory.
-- Make it monomorphic to take advantages of memoization. Compare
--
-- >>> primes !! 1000000 :: Prime Int -- (5.32 secs, 6,945,267,496 bytes)
-- Prime 15485867
-- >>> primes !! 1000000 :: Prime Int -- (5.19 secs, 6,945,267,496 bytes)
-- Prime 15485867
--
-- against
--
-- >>> let primes' = primes :: [Prime Int]
-- >>> primes' !! 1000000 :: Prime Int -- (5.29 secs, 6,945,269,856 bytes)
-- Prime 15485867
-- >>> primes' !! 1000000 :: Prime Int -- (0.02 secs, 336,232 bytes)
-- Prime 15485867
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

-- | List of primes in the form of a list of 'PrimeSieve's, more compact than
--   'primes', thus it may be better to use @'psieveList' >>= 'primeList'@
--   than keeping the list of primes alive during the entire run.
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     -- prime #647, 644 of them to use
    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
        MVector s Bit
sieve <- Integer -> ST s (MVector s Bit)
forall s. Integer -> ST s (MVector s Bit)
sieveTo (Integer
4801 :: Integer)
        MVector s Word64
new <- Int -> ST s (MVector (PrimState (ST s)) Word64)
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> m (MVector (PrimState m) a)
MU.unsafeNew Int
1288 :: ST s (MU.MVector s Word64)
        let fill :: Int -> Int -> m (MVector s Word64)
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    -- index of 4801 = 159*30 + 31 ~> 159*8+7
              | Bool
otherwise = do
                Bit Bool
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 Bool
p
                  then do
                    let !i :: Int
i = Int
indx Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
jMASK
                        k :: Int
k = Int
indx Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftR` Int
jBITS
                        strt1 :: Int
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 :: Word64
strt = Int -> Word64
intToWord64 (Int
strt1 Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
forall a. Num a => a
iXMASK)
                        !skip :: Word64
skip = Int -> Word64
intToWord64 (Int
strt1 Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftR` Int
iXBITS)
                        !ixes :: Word64
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
                    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
j Word64
skip
                    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
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Word64
ixes
                    Int -> Int -> m (MVector s Word64)
fill (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) (Int
indxInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
                  else Int -> Int -> m (MVector s Word64)
fill Int
j (Int
indxInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
        MVector s Word64
vec <- Int -> Int -> ST s (MVector s Word64)
forall {m :: * -> *}.
(PrimState m ~ s, PrimMonad m) =>
Int -> Int -> m (MVector s Word64)
fill Int
0 Int
0
        MVector (PrimState (ST s)) Word64 -> ST s (Vector Word64)
forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
MVector (PrimState m) a -> m (Vector a)
U.unsafeFreeze MVector s Word64
MVector (PrimState (ST s)) Word64
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
            MVector s Word64
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)
            MVector s Bit
bs0 <- MVector s Word64 -> ST s (MVector s Bit)
forall s. MVector s Word64 -> ST s (MVector s Bit)
slice MVector s Word64
cch
            Vector Word64
fcch <- MVector (PrimState (ST s)) Word64 -> ST s (Vector Word64)
forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
MVector (PrimState m) a -> m (Vector a)
U.unsafeFreeze MVector s Word64
MVector (PrimState (ST s)) Word64
cch
            Vector Bit
fbs0 <- 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
MVector (PrimState (ST s)) Bit
bs0
            (Vector Word64, Vector Bit) -> ST s (Vector Word64, Vector Bit)
forall a. a -> ST s a
forall (m :: * -> *) a. Monad m => a -> m a
return (Vector Word64
fcch, Vector Bit
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
            MVector s Word64
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
            MVector s Bit
bs0 <- MVector s Word64 -> ST s (MVector s Bit)
forall s. MVector s Word64 -> ST s (MVector s Bit)
slice MVector s Word64
cch
            Vector Word64
fcch <- MVector (PrimState (ST s)) Word64 -> ST s (Vector Word64)
forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
MVector (PrimState m) a -> m (Vector a)
U.unsafeFreeze MVector s Word64
MVector (PrimState (ST s)) Word64
cch
            Vector Bit
fbs0 <- 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
MVector (PrimState (ST s)) Bit
bs0
            (Vector Word64, Vector Bit) -> ST s (Vector Word64, Vector Bit)
forall a. a -> ST s a
forall (m :: * -> *) a. Monad m => a -> m a
return (Vector Word64
fcch, Vector Bit
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
    MVector s Bit
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 -> m (MVector s Bit)
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
            Word64
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 Word64
w Word64 -> Word64 -> Bool
forall a. Eq a => a -> a -> Bool
/= Word64
0
              then 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
cache Int
pr (Word64
wWord64 -> Word64 -> Word64
forall a. Num a => a -> a -> a
-Word64
1)
              else do
                Word64
ixes <- 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
prInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
                let !stj :: Int
stj = Word64 -> Int
word64ToInt Word64
ixes Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
forall a. Num a => a
iXJMASK   -- position of multiple and index of cofactor
                    !ixw :: Int
ixw = Word64 -> Int
word64ToInt (Word64
ixes Word64 -> Int -> Word64
forall a. Bits a => a -> Int -> a
`shiftR` Int
iXJBITS)  -- prime data, up to 41 bits
                    !i :: Int
i = Int
ixw Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
jMASK
                    !k :: Int
k = Int
ixw Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
i        -- On 32-bits, k > 44717396 means overflow is possible in tick
                    !o :: Int
o = Int
i Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftL` Int
jBITS
                    !j :: Int
j = Int
stj Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
jMASK          -- index of cofactor
                    !s :: Int
s = Int
stj Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftR` Int
jBITS     -- index of first multiple to tick off
                (Int
n, Int
u) <- Int -> Int -> Int -> Int -> m (Int, Int)
forall {m :: * -> *}.
(PrimState m ~ s, PrimMonad m) =>
Int -> Int -> Int -> Int -> m (Int, Int)
tick Int
k Int
o Int
j Int
s
                let !skip :: Word64
skip = Int -> Word64
intToWord64 (Int
n Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftR` Int
iXBITS)
                    !strt :: Word64
strt = Int -> Word64
intToWord64 (Int
n Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
forall a. Num a => a
iXMASK)
                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
cache Int
pr Word64
skip
                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
cache (Int
prInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) ((Word64
ixes Word64 -> Word64 -> Word64
forall a. Bits a => a -> a -> a
.&. Word64 -> Word64
forall a. Bits a => a -> a
complement Word64
forall a. Num a => a
iXJMASK) 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
u)
            Int -> m (MVector s Bit)
treat (Int
prInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2)
        tick :: Int -> Int -> Int -> Int -> m (Int, Int)
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 Bool
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
            Bool -> m () -> m ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when Bool
p (MVector (PrimState m) Bit -> Int -> Bit -> m ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
MU.unsafeWrite MVector s Bit
MVector (PrimState m) Bit
sieve Int
ix (Bool -> Bit
Bit Bool
False))
            Int -> Int -> Int -> Int -> m (Int, Int)
tick Int
stp Int
off ((Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
jMASK) (Int
ix Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
stpInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int -> Int
delta Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int -> Int
tau (Int
offInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
j))
    Int -> ST s (MVector s Bit)
forall {m :: * -> *}.
(PrimState m ~ s, PrimMonad m) =>
Int -> m (MVector s Bit)
treat Int
0

-- | Sieve up to bound in one go.
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
        MVector (PrimState (ST s)) Bit
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 -> Int -> Int
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 -> Int -> Int -> Int -> m ()
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 Bool
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
                Bool -> m () -> m ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when Bool
p (MVector (PrimState m) Bit -> Int -> Bit -> m ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
MU.unsafeWrite MVector (PrimState m) Bit
MVector (PrimState (ST s)) Bit
ar Int
ix (Bool -> Bit
Bit Bool
False))
                Int -> Int -> Int -> Int -> m ()
tick Int
stp Int
off ((Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
jMASK) (Int
ix Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
stpInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int -> Int
delta Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int -> Int
tau (Int
offInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
j))
            sift :: Int -> m (MVector (PrimState (ST s)) Bit)
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 Bool
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
                Bool -> m () -> m ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when Bool
p  (do let i :: Int
i = Int
ix Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
jMASK
                                k :: Int
k = Int
ix Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftR` Int
jBITS
                                !off :: Int
off = Int
i Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftL` Int
jBITS
                                !stp :: Int
stp = Int
ix Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
i
                            Int -> Int -> Int -> Int -> m ()
forall {m :: * -> *}.
(PrimState m ~ PrimState (ST s), PrimMonad m) =>
Int -> Int -> Int -> Int -> m ()
tick Int
stp Int
off Int
i (Int -> Int -> Int
start Int
k Int
i))
                Int -> m (MVector (PrimState (ST s)) Bit)
sift (Int
ixInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
        Int -> ST s (MVector (PrimState (ST s)) Bit)
forall {m :: * -> *}.
(PrimState m ~ PrimState (ST s), PrimMonad m) =>
Int -> m (MVector (PrimState (ST s)) Bit)
sift Int
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
    MVector s Bit
sieve <- Integer -> ST s (MVector s Bit)
forall s. Integer -> ST s (MVector s Bit)
sieveTo Integer
nlim       -- Implement SieveFromTo for this, it's pretty wasteful when nlim isn't
    let hi :: Int
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
    Int
more <- Int -> Int -> MVector s Bit -> ST s Int
forall s. Int -> Int -> MVector s Bit -> ST s Int
countFromToWd Int
start Int
hi MVector s Bit
sieve
    MVector s Word64
new <- Int -> ST s (MVector (PrimState (ST s)) Word64)
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> m (MVector (PrimState m) a)
MU.unsafeNew (Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
num Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
more) :: ST s (MU.MVector s Word64)
    let copy :: Int -> m ()
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)
    Int -> ST s ()
forall {m :: * -> *}. (PrimState m ~ s, PrimMonad m) => Int -> m ()
copy Int
0
    let fill :: Int -> Int -> m (MVector s Word64)
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 Bool
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 Bool
p
              then do
                let !i :: Int
i = Int
indx Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
jMASK
                    k :: Integer
                    k :: Integer
k = Int -> Integer
intToInteger (Int
indx Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftR` Int
jBITS)
                    strt0 :: Integer
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
strt1 = Integer
strt0 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
offset
                    !strt :: Word64
strt = Integer -> Word64
integerToWord64 Integer
strt1 Word64 -> Word64 -> Word64
forall a. Bits a => a -> a -> a
.&. Word64
forall a. Num a => a
iXMASK
                    !skip :: Word64
skip = Integer -> Word64
integerToWord64 (Integer
strt1 Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
`shiftR` Int
iXBITS)
                    !ixes :: Word64
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
                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
j Word64
skip
                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
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Word64
ixes
                Int -> Int -> m (MVector s Word64)
fill (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) (Int
indxInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
              else Int -> Int -> m (MVector s Word64)
fill Int
j (Int
indxInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
    Int -> Int -> ST s (MVector s Word64)
forall {m :: * -> *}.
(PrimState m ~ s, PrimMonad m) =>
Int -> Int -> m (MVector s Word64)
fill (Int
numInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int
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' n@ creates the list of 'PrimeSieve's starting roughly
--   at @n@. Due to the organisation of the sieve, the list may contain
--   a few primes less than @n@.
--   This form uses less memory than @['Integer']@, hence it may be preferable
--   to use this if it is to be reused.
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 -- beware arithmetic underflow
      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
          MVector s Bit
sieve <- Integer -> ST s (MVector s Bit)
forall s. Integer -> ST s (MVector s Bit)
sieveTo Integer
plim
          let (Int
lo,Int
hi) = (Int
0, 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)
          Int
pct <- Int -> Int -> MVector s Bit -> ST s Int
forall s. Int -> Int -> MVector s Bit -> ST s Int
countFromToWd Int
lo Int
hi MVector s Bit
sieve
          MVector s Word64
new <- Int -> ST s (MVector (PrimState (ST s)) Word64)
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> m (MVector (PrimState m) a)
MU.unsafeNew (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
pct) ::  ST s (MU.MVector s Word64)
          let fill :: Int -> Int -> m (MVector s Word64)
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 Bool
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 Bool
isPr
                    then do
                      let !i :: Int
i = Int
indx Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
jMASK
                          !moff :: Int
moff = Int
i Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftL` Int
jBITS
                          k :: Integer
                          k :: Integer
k = Int -> Integer
intToInteger (Int
indx Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftR` Int
jBITS)
                          p :: Integer
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
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
                          (Integer
skp0,Integer
q1) = Integer
q0 Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
`quotRem` Int -> Integer
intToInteger Int
sieveRange
                          (Int
b0,Int
r0)
                              | Integer
q1 Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0   = (-Int
1,Int
6)
                              | Integer
q1 Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
7    = (-Int
1,Int
7)
                              | Bool
otherwise = Int -> (Int, Int)
forall a. Integral a => a -> (Int, Int)
idxPr (Integer -> Int
integerToInt Integer
q1 :: Int)
                          (Int
b1,Int
r1) | Int
r0 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
7 = (Int
b0Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1,Int
0)
                                  | Bool
otherwise = (Int
b0,Int
r0Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
                          b2 :: Integer
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
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
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)
                          (Integer
strt2,Int
r2)
                              | Integer
p Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
ssr   = (Integer
strt0 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
bitOff,Int
r1)
                              | Bool
otherwise = (Integer
strt1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
bitOff, Int
i)
                          !strt :: Word64
strt = Integer -> Word64
integerToWord64 Integer
strt2 Word64 -> Word64 -> Word64
forall a. Bits a => a -> a -> a
.&. Word64
forall a. Num a => a
iXMASK
                          !skip :: Word64
skip = Integer -> Word64
integerToWord64 (Integer
strt2 Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
`shiftR` Int
iXBITS)
                          !ixes :: Word64
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
                      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
j Word64
skip
                      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
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Word64
ixes
                      Int -> Int -> m (MVector s Word64)
fill (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) (Int
indxInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
                    else Int -> Int -> m (MVector s Word64)
fill Int
j (Int
indxInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
          MVector s Word64
vec <- Int -> Int -> ST s (MVector s Word64)
forall {m :: * -> *}.
(PrimState m ~ s, PrimMonad m) =>
Int -> Int -> m (MVector s Word64)
fill Int
0 Int
0
          MVector (PrimState (ST s)) Word64 -> ST s (Vector Word64)
forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
MVector (PrimState m) a -> m (Vector a)
U.unsafeFreeze MVector s Word64
MVector (PrimState (ST s)) Word64
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
        ]