-- |
-- 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.Array.Base
import Data.Array.ST
import Data.Bits
import Data.Coerce
import Data.Proxy
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

wSHFT :: (Bits a, Num a) => a
wSHFT :: forall a. (Bits a, Num a) => a
wSHFT = if Word -> Int
forall b. FiniteBits b => b -> Int
finiteBitSize (Word
0 :: Word) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
64 then a
6 else a
5

-- | Compact store of primality flags.
data PrimeSieve = PS !Integer {-# UNPACK #-} !(UArray Int Bool)

-- | 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 -> UArray Int Bool -> PrimeSieve
PS Integer
0 ((forall s. ST s (STUArray s Int Bool)) -> UArray Int Bool
forall i e. (forall s. ST s (STUArray s i e)) -> UArray i e
runSTUArray ((forall s. ST s (STUArray s Int Bool)) -> UArray Int Bool)
-> (forall s. ST s (STUArray s Int Bool)) -> UArray Int Bool
forall a b. (a -> b) -> a -> b
$ Integer -> ST s (STUArray s Int Bool)
forall s. Integer -> ST s (STUArray s Int Bool)
sieveTo Integer
bound)

-- | 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 UArray Int Bool
_)
  | 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 UArray Int Bool
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 (UArray Int Bool -> Int -> Bool
forall i. Ix i => UArray i Bool -> Int -> Bool
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> Int -> e
unsafeAt UArray Int Bool
bs) [Int
lo..Int
hi]
  where
    (Int
lo, Int
hi) = UArray Int Bool -> (Int, Int)
forall i. Ix i => UArray i Bool -> (i, i)
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> (i, i)
bounds UArray Int Bool
bs

-- | 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
-> UArray Int Word64
-> [PrimeSieve]
makeSieves Integer
plim Integer
sqlim Integer
0 Integer
0 UArray Int 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 :: UArray Int Word64
cache = (forall s. ST s (STUArray s Int Word64)) -> UArray Int Word64
forall i e. (forall s. ST s (STUArray s i e)) -> UArray i e
runSTUArray ((forall s. ST s (STUArray s Int Word64)) -> UArray Int Word64)
-> (forall s. ST s (STUArray s Int Word64)) -> UArray Int Word64
forall a b. (a -> b) -> a -> b
$ do
        STUArray s Int Bool
sieve <- Integer -> ST s (STUArray s Int Bool)
forall s. Integer -> ST s (STUArray s Int Bool)
sieveTo (Integer
4801 :: Integer)
        STUArray s Int Word64
new <- (Int, Int) -> ST s (STUArray s Int Word64)
forall i. Ix i => (i, i) -> ST s (STUArray s i Word64)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> m (a i e)
unsafeNewArray_ (Int
0,Int
1287) :: ST s (STUArray s Int Word64)
        let fill :: Int -> Int -> m (STUArray s Int Word64)
fill Int
j Int
indx
              | Int
1279 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
indx = STUArray s Int Word64 -> m (STUArray s Int Word64)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return STUArray s Int Word64
new    -- index of 4801 = 159*30 + 31 ~> 159*8+7
              | Bool
otherwise = do
                Bool
p <- STUArray s Int Bool -> Int -> m Bool
forall i. Ix i => STUArray s i Bool -> Int -> m Bool
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> m e
unsafeRead STUArray s Int Bool
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
                    STUArray s Int Word64 -> Int -> Word64 -> m ()
forall i. Ix i => STUArray s i Word64 -> Int -> Word64 -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> e -> m ()
unsafeWrite STUArray s Int Word64
new Int
j Word64
skip
                    STUArray s Int Word64 -> Int -> Word64 -> m ()
forall i. Ix i => STUArray s i Word64 -> Int -> Word64 -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> e -> m ()
unsafeWrite STUArray s Int Word64
new (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Word64
ixes
                    Int -> Int -> m (STUArray s Int 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 (STUArray s Int Word64)
fill Int
j (Int
indxInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
        Int -> Int -> ST s (STUArray s Int Word64)
forall {m :: * -> *}.
(MArray (STUArray s) Bool m, MArray (STUArray s) Word64 m) =>
Int -> Int -> m (STUArray s Int Word64)
fill Int
0 Int
0

makeSieves :: Integer -> Integer -> Integer -> Integer -> UArray Int Word64 -> [PrimeSieve]
makeSieves :: Integer
-> Integer
-> Integer
-> Integer
-> UArray Int Word64
-> [PrimeSieve]
makeSieves Integer
plim Integer
sqlim Integer
bitOff Integer
valOff UArray Int Word64
cache
  | Integer
valOff' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
sqlim =
      let (UArray Int Word64
nc, UArray Int Bool
bs) = (forall s. ST s (UArray Int Word64, UArray Int Bool))
-> (UArray Int Word64, UArray Int Bool)
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s (UArray Int Word64, UArray Int Bool))
 -> (UArray Int Word64, UArray Int Bool))
-> (forall s. ST s (UArray Int Word64, UArray Int Bool))
-> (UArray Int Word64, UArray Int Bool)
forall a b. (a -> b) -> a -> b
$ do
            STUArray s Int Word64
cch <- UArray Int Word64 -> ST s (STUArray s Int Word64)
forall i (a :: * -> * -> *) e (b :: * -> * -> *) (m :: * -> *).
(Ix i, IArray a e, MArray b e m) =>
a i e -> m (b i e)
unsafeThaw UArray Int Word64
cache :: ST s (STUArray s Int Word64)
            STUArray s Int Bool
bs0 <- STUArray s Int Word64 -> ST s (STUArray s Int Bool)
forall s. STUArray s Int Word64 -> ST s (STUArray s Int Bool)
slice STUArray s Int Word64
cch
            UArray Int Word64
fcch <- STUArray s Int Word64 -> ST s (UArray Int Word64)
forall i (a :: * -> * -> *) e (m :: * -> *) (b :: * -> * -> *).
(Ix i, MArray a e m, IArray b e) =>
a i e -> m (b i e)
unsafeFreeze STUArray s Int Word64
cch
            UArray Int Bool
fbs0 <- STUArray s Int Bool -> ST s (UArray Int Bool)
forall i (a :: * -> * -> *) e (m :: * -> *) (b :: * -> * -> *).
(Ix i, MArray a e m, IArray b e) =>
a i e -> m (b i e)
unsafeFreeze STUArray s Int Bool
bs0
            (UArray Int Word64, UArray Int Bool)
-> ST s (UArray Int Word64, UArray Int Bool)
forall a. a -> ST s a
forall (m :: * -> *) a. Monad m => a -> m a
return (UArray Int Word64
fcch, UArray Int Bool
fbs0)
      in Integer -> UArray Int Bool -> PrimeSieve
PS Integer
valOff UArray Int Bool
bs PrimeSieve -> [PrimeSieve] -> [PrimeSieve]
forall a. a -> [a] -> [a]
: Integer
-> Integer
-> Integer
-> Integer
-> UArray Int Word64
-> [PrimeSieve]
makeSieves Integer
plim Integer
sqlim Integer
bitOff' Integer
valOff' UArray Int 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'
          (UArray Int Word64
nc,UArray Int Bool
bs) = (forall s. ST s (UArray Int Word64, UArray Int Bool))
-> (UArray Int Word64, UArray Int Bool)
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s (UArray Int Word64, UArray Int Bool))
 -> (UArray Int Word64, UArray Int Bool))
-> (forall s. ST s (UArray Int Word64, UArray Int Bool))
-> (UArray Int Word64, UArray Int Bool)
forall a b. (a -> b) -> a -> b
$ do
            STUArray s Int Word64
cch <- Integer
-> Integer -> UArray Int Word64 -> ST s (STUArray s Int Word64)
forall s.
Integer
-> Integer -> UArray Int Word64 -> ST s (STUArray s Int Word64)
growCache Integer
bitOff Integer
plim UArray Int Word64
cache
            STUArray s Int Bool
bs0 <- STUArray s Int Word64 -> ST s (STUArray s Int Bool)
forall s. STUArray s Int Word64 -> ST s (STUArray s Int Bool)
slice STUArray s Int Word64
cch
            UArray Int Word64
fcch <- STUArray s Int Word64 -> ST s (UArray Int Word64)
forall i (a :: * -> * -> *) e (m :: * -> *) (b :: * -> * -> *).
(Ix i, MArray a e m, IArray b e) =>
a i e -> m (b i e)
unsafeFreeze STUArray s Int Word64
cch
            UArray Int Bool
fbs0 <- STUArray s Int Bool -> ST s (UArray Int Bool)
forall i (a :: * -> * -> *) e (m :: * -> *) (b :: * -> * -> *).
(Ix i, MArray a e m, IArray b e) =>
a i e -> m (b i e)
unsafeFreeze STUArray s Int Bool
bs0
            (UArray Int Word64, UArray Int Bool)
-> ST s (UArray Int Word64, UArray Int Bool)
forall a. a -> ST s a
forall (m :: * -> *) a. Monad m => a -> m a
return (UArray Int Word64
fcch, UArray Int Bool
fbs0)
      in Integer -> UArray Int Bool -> PrimeSieve
PS Integer
valOff UArray Int Bool
bs PrimeSieve -> [PrimeSieve] -> [PrimeSieve]
forall a. a -> [a] -> [a]
: Integer
-> Integer
-> Integer
-> Integer
-> UArray Int Word64
-> [PrimeSieve]
makeSieves Integer
plim' Integer
sqlim' Integer
bitOff' Integer
valOff' UArray Int 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 :: STUArray s Int Word64 -> ST s (STUArray s Int Bool)
slice :: forall s. STUArray s Int Word64 -> ST s (STUArray s Int Bool)
slice STUArray s Int Word64
cache = do
    Int
hi <- (Int, Int) -> Int
forall a b. (a, b) -> b
snd ((Int, Int) -> Int) -> ST s (Int, Int) -> 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` STUArray s Int Word64 -> ST s (Int, Int)
forall i. Ix i => STUArray s i Word64 -> ST s (i, i)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds STUArray s Int Word64
cache
    STUArray s Int Bool
sieve <- (Int, Int) -> Bool -> ST s (STUArray s Int Bool)
forall i. Ix i => (i, i) -> Bool -> ST s (STUArray s i Bool)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> e -> m (a i e)
newArray (Int
0,Int
lastIndex) Bool
True
    let treat :: Int -> m (STUArray s Int Bool)
treat Int
pr
          | Int
hi Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
pr     = STUArray s Int Bool -> m (STUArray s Int Bool)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return STUArray s Int Bool
sieve
          | Bool
otherwise   = do
            Word64
w <- STUArray s Int Word64 -> Int -> m Word64
forall i. Ix i => STUArray s i Word64 -> Int -> m Word64
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> m e
unsafeRead STUArray s Int Word64
cache Int
pr
            if Word64
w Word64 -> Word64 -> Bool
forall a. Eq a => a -> a -> Bool
/= Word64
0
              then STUArray s Int Word64 -> Int -> Word64 -> m ()
forall i. Ix i => STUArray s i Word64 -> Int -> Word64 -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> e -> m ()
unsafeWrite STUArray s Int Word64
cache Int
pr (Word64
wWord64 -> Word64 -> Word64
forall a. Num a => a -> a -> a
-Word64
1)
              else do
                Word64
ixes <- STUArray s Int Word64 -> Int -> m Word64
forall i. Ix i => STUArray s i Word64 -> Int -> m Word64
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> m e
unsafeRead STUArray s Int 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 :: * -> *}.
MArray (STUArray s) Bool 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)
                STUArray s Int Word64 -> Int -> Word64 -> m ()
forall i. Ix i => STUArray s i Word64 -> Int -> Word64 -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> e -> m ()
unsafeWrite STUArray s Int Word64
cache Int
pr Word64
skip
                STUArray s Int Word64 -> Int -> Word64 -> m ()
forall i. Ix i => STUArray s i Word64 -> Int -> Word64 -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> e -> m ()
unsafeWrite STUArray s Int 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 (STUArray s Int Bool)
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
            Bool
p <- STUArray s Int Bool -> Int -> m Bool
forall i. Ix i => STUArray s i Bool -> Int -> m Bool
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> m e
unsafeRead STUArray s Int Bool
sieve Int
ix
            Bool -> m () -> m ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when Bool
p (STUArray s Int Bool -> Int -> Bool -> m ()
forall i. Ix i => STUArray s i Bool -> Int -> Bool -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> e -> m ()
unsafeWrite STUArray s Int Bool
sieve Int
ix 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 (STUArray s Int Bool)
forall {m :: * -> *}.
(MArray (STUArray s) Bool m, MArray (STUArray s) Word64 m) =>
Int -> m (STUArray s Int Bool)
treat Int
0

-- | Sieve up to bound in one go.
sieveTo :: Integer -> ST s (STUArray s Int Bool)
sieveTo :: forall s. Integer -> ST s (STUArray s Int Bool)
sieveTo Integer
bound = ST s (STUArray s Int Bool)
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 (STUArray s Int Bool)
arr = do
        STUArray s Int Bool
ar <- (Int, Int) -> Bool -> ST s (STUArray s Int Bool)
forall i. Ix i => (i, i) -> Bool -> ST s (STUArray s i Bool)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> e -> m (a i e)
newArray (Int
0,Int
mxidx) 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
                Bool
p <- STUArray s Int Bool -> Int -> m Bool
forall i. Ix i => STUArray s i Bool -> Int -> m Bool
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> m e
unsafeRead STUArray s Int Bool
ar Int
ix
                Bool -> m () -> m ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when Bool
p (STUArray s Int Bool -> Int -> Bool -> m ()
forall i. Ix i => STUArray s i Bool -> Int -> Bool -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> e -> m ()
unsafeWrite STUArray s Int Bool
ar Int
ix 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 (STUArray s Int Bool)
sift Int
ix
              | Int
svbd Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
ix = STUArray s Int Bool -> m (STUArray s Int Bool)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return STUArray s Int Bool
ar
              | Bool
otherwise = do
                Bool
p <- STUArray s Int Bool -> Int -> m Bool
forall i. Ix i => STUArray s i Bool -> Int -> m Bool
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> m e
unsafeRead STUArray s Int Bool
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 :: * -> *}.
MArray (STUArray s) Bool m =>
Int -> Int -> Int -> Int -> m ()
tick Int
stp Int
off Int
i (Int -> Int -> Int
start Int
k Int
i))
                Int -> m (STUArray s Int Bool)
sift (Int
ixInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
        Int -> ST s (STUArray s Int Bool)
forall {m :: * -> *}.
MArray (STUArray s) Bool m =>
Int -> m (STUArray s Int Bool)
sift Int
0

growCache :: Integer -> Integer -> UArray Int Word64 -> ST s (STUArray s Int Word64)
growCache :: forall s.
Integer
-> Integer -> UArray Int Word64 -> ST s (STUArray s Int Word64)
growCache Integer
offset Integer
plim UArray Int Word64
old = do
    let (Int
_,Int
num) = UArray Int Word64 -> (Int, Int)
forall i. Ix i => UArray i Word64 -> (i, i)
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> (i, i)
bounds UArray Int Word64
old
        (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
    STUArray s Int Bool
sieve <- Integer -> ST s (STUArray s Int Bool)
forall s. Integer -> ST s (STUArray s Int Bool)
sieveTo Integer
nlim       -- Implement SieveFromTo for this, it's pretty wasteful when nlim isn't
    (Int
_,Int
hi) <- STUArray s Int Bool -> ST s (Int, Int)
forall i. Ix i => STUArray s i Bool -> ST s (i, i)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds STUArray s Int Bool
sieve   -- very small anymore
    Int
more <- Int -> Int -> STUArray s Int Bool -> ST s Int
forall s. Int -> Int -> STUArray s Int Bool -> ST s Int
countFromToWd Int
start Int
hi STUArray s Int Bool
sieve
    STUArray s Int Word64
new <- (Int, Int) -> ST s (STUArray s Int Word64)
forall i. Ix i => (i, i) -> ST s (STUArray s i Word64)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> m (a i e)
unsafeNewArray_ (Int
0,Int
numInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
more) :: ST s (STUArray s Int 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
            STUArray s Int Word64 -> Int -> Word64 -> m ()
forall i. Ix i => STUArray s i Word64 -> Int -> Word64 -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> e -> m ()
unsafeWrite STUArray s Int Word64
new Int
i (UArray Int Word64
old UArray Int Word64 -> Int -> Word64
forall i. Ix i => UArray i Word64 -> Int -> Word64
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> Int -> e
`unsafeAt` Int
i)
            Int -> m ()
copy (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
    Int -> ST s ()
forall {m :: * -> *}. MArray (STUArray s) Word64 m => Int -> m ()
copy Int
0
    let fill :: Int -> Int -> m (STUArray s Int Word64)
fill Int
j Int
indx
          | Int
hi Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
indx = STUArray s Int Word64 -> m (STUArray s Int Word64)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return STUArray s Int Word64
new
          | Bool
otherwise = do
            Bool
p <- STUArray s Int Bool -> Int -> m Bool
forall i. Ix i => STUArray s i Bool -> Int -> m Bool
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> m e
unsafeRead STUArray s Int Bool
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
                STUArray s Int Word64 -> Int -> Word64 -> m ()
forall i. Ix i => STUArray s i Word64 -> Int -> Word64 -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> e -> m ()
unsafeWrite STUArray s Int Word64
new Int
j Word64
skip
                STUArray s Int Word64 -> Int -> Word64 -> m ()
forall i. Ix i => STUArray s i Word64 -> Int -> Word64 -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> e -> m ()
unsafeWrite STUArray s Int Word64
new (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Word64
ixes
                Int -> Int -> m (STUArray s Int 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 (STUArray s Int Word64)
fill Int
j (Int
indxInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
    Int -> Int -> ST s (STUArray s Int Word64)
forall {m :: * -> *}.
(MArray (STUArray s) Bool m, MArray (STUArray s) Word64 m) =>
Int -> Int -> m (STUArray s Int Word64)
fill (Int
numInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int
start

-- Danger: relies on start and end being the first resp. last
-- index in a Word
-- Do not use except in growCache and psieveFrom
{-# INLINE countFromToWd #-}
countFromToWd :: Int -> Int -> STUArray s Int Bool -> ST s Int
countFromToWd :: forall s. Int -> Int -> STUArray s Int Bool -> ST s Int
countFromToWd Int
start Int
end STUArray s Int Bool
ba = do
    STUArray s Int Word
wa <- (STUArray s Int Bool -> ST s (STUArray s Int Word)
forall {s}. STUArray s Int Bool -> ST s (STUArray s Int Word)
forall s ix a b. STUArray s ix a -> ST s (STUArray s ix b)
castSTUArray :: STUArray s Int Bool -> ST s (STUArray s Int Word)) STUArray s Int Bool
ba
    let !sb :: Int
sb = Int
start Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftR` Int
forall a. (Bits a, Num a) => a
wSHFT
        !eb :: Int
eb = Int
end Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftR` Int
forall a. (Bits a, Num a) => a
wSHFT
        count :: Int -> Int -> m Int
count !Int
acc Int
i
          | Int
eb Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
i    = Int -> m Int
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return Int
acc
          | Bool
otherwise = do
            Word
w <- STUArray s Int Word -> Int -> m Word
forall i. Ix i => STUArray s i Word -> Int -> m Word
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> m e
unsafeRead STUArray s Int Word
wa Int
i
            Int -> Int -> m Int
count (Int
acc Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Word -> Int
forall a. Bits a => a -> Int
popCount Word
w) (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
    Int -> Int -> ST s Int
forall {m :: * -> *}.
MArray (STUArray s) Word m =>
Int -> Int -> m Int
count Int
0 Int
sb

-- | @'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
-> UArray Int Word64
-> [PrimeSieve]
makeSieves Integer
plim Integer
sqlim Integer
bitOff Integer
valOff UArray Int 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 :: UArray Int Word64
cache = (forall s. ST s (STUArray s Int Word64)) -> UArray Int Word64
forall i e. (forall s. ST s (STUArray s i e)) -> UArray i e
runSTUArray ((forall s. ST s (STUArray s Int Word64)) -> UArray Int Word64)
-> (forall s. ST s (STUArray s Int Word64)) -> UArray Int Word64
forall a b. (a -> b) -> a -> b
$ do
          STUArray s Int Bool
sieve <- Integer -> ST s (STUArray s Int Bool)
forall s. Integer -> ST s (STUArray s Int Bool)
sieveTo Integer
plim
          (Int
lo,Int
hi) <- STUArray s Int Bool -> ST s (Int, Int)
forall i. Ix i => STUArray s i Bool -> ST s (i, i)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds STUArray s Int Bool
sieve
          Int
pct <- Int -> Int -> STUArray s Int Bool -> ST s Int
forall s. Int -> Int -> STUArray s Int Bool -> ST s Int
countFromToWd Int
lo Int
hi STUArray s Int Bool
sieve
          STUArray s Int Word64
new <- (Int, Int) -> ST s (STUArray s Int Word64)
forall i. Ix i => (i, i) -> ST s (STUArray s i Word64)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> m (a i e)
unsafeNewArray_ (Int
0,Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
pctInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) ::  ST s (STUArray s Int Word64)
          let fill :: Int -> Int -> m (STUArray s Int Word64)
fill Int
j Int
indx
                | Int
hi Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
indx = STUArray s Int Word64 -> m (STUArray s Int Word64)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return STUArray s Int Word64
new
                | Bool
otherwise = do
                  Bool
isPr <- STUArray s Int Bool -> Int -> m Bool
forall i. Ix i => STUArray s i Bool -> Int -> m Bool
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> m e
unsafeRead STUArray s Int Bool
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
                      STUArray s Int Word64 -> Int -> Word64 -> m ()
forall i. Ix i => STUArray s i Word64 -> Int -> Word64 -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> e -> m ()
unsafeWrite STUArray s Int Word64
new Int
j Word64
skip
                      STUArray s Int Word64 -> Int -> Word64 -> m ()
forall i. Ix i => STUArray s i Word64 -> Int -> Word64 -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> Int -> e -> m ()
unsafeWrite STUArray s Int Word64
new (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Word64
ixes
                      Int -> Int -> m (STUArray s Int 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 (STUArray s Int Word64)
fill Int
j (Int
indxInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
          Int -> Int -> ST s (STUArray s Int Word64)
forall {m :: * -> *}.
(MArray (STUArray s) Bool m, MArray (STUArray s) Word64 m) =>
Int -> Int -> m (STUArray s Int Word64)
fill Int
0 Int
0


{-# INLINE delta #-}
delta :: Int -> Int
delta :: Int -> Int
delta = UArray Int Int -> Int -> Int
forall i. Ix i => UArray i Int -> Int -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> Int -> e
unsafeAt UArray Int Int
deltas

deltas :: UArray Int Int
deltas :: UArray Int Int
deltas = (Int, Int) -> [Int] -> UArray Int Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
(i, i) -> [e] -> a i e
listArray (Int
0,Int
7) [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 = UArray Int Int -> Int -> Int
forall i. Ix i => UArray i Int -> Int -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> Int -> e
unsafeAt UArray Int Int
taus

taus :: UArray Int Int
taus :: UArray Int Int
taus = (Int, Int) -> [Int] -> UArray Int Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
(i, i) -> [e] -> a i e
listArray (Int
0,Int
63)
        [  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 = UArray Int Int -> Int -> Int
forall i. Ix i => UArray i Int -> Int -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> Int -> e
unsafeAt UArray Int Int
startByte

startByte :: UArray Int Int
startByte :: UArray Int Int
startByte = (Int, Int) -> [Int] -> UArray Int Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
(i, i) -> [e] -> a i e
listArray (Int
0,Int
7) [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 = UArray Int Int -> Int -> Int
forall i. Ix i => UArray i Int -> Int -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> Int -> e
unsafeAt UArray Int Int
startIdx

startIdx :: UArray Int Int
startIdx :: UArray Int Int
startIdx = (Int, Int) -> [Int] -> UArray Int Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
(i, i) -> [e] -> a i e
listArray (Int
0,Int
7) [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 = UArray Int Int -> Int -> Int
forall i. Ix i => UArray i Int -> Int -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> Int -> e
unsafeAt UArray Int Int
mArr

{-# INLINE nu #-}
nu :: Int -> Int
nu :: Int -> Int
nu = UArray Int Int -> Int -> Int
forall i. Ix i => UArray i Int -> Int -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> Int -> e
unsafeAt UArray Int Int
nArr

mArr :: UArray Int Int
mArr :: UArray Int Int
mArr = (Int, Int) -> [Int] -> UArray Int Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
(i, i) -> [e] -> a i e
listArray (Int
0,Int
63)
        [ 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 :: UArray Int Int
nArr :: UArray Int Int
nArr = (Int, Int) -> [Int] -> UArray Int Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
(i, i) -> [e] -> a i e
listArray (Int
0,Int
63)
        [ 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
        ]