diff --git a/Math/NumberTheory/Primes.hs b/Math/NumberTheory/Primes.hs
index a60f5bc3f..7168f2b7e 100644
--- a/Math/NumberTheory/Primes.hs
+++ b/Math/NumberTheory/Primes.hs
@@ -20,6 +20,8 @@ module Math.NumberTheory.Primes
     , factorBack
     , -- * Old interface
       primes
+    , -- * Temporary
+      module Math.NumberTheory.Primes.Sieve.Atkin
     ) where
 
 import Data.Bits
@@ -31,6 +33,7 @@ import Numeric.Natural
 import Math.NumberTheory.Primes.Counting (nthPrime, primeCount)
 import qualified Math.NumberTheory.Primes.Factorisation.Montgomery as F (factorise)
 import qualified Math.NumberTheory.Primes.Testing.Probabilistic as T (isPrime)
+import Math.NumberTheory.Primes.Sieve.Atkin
 import Math.NumberTheory.Primes.Sieve.Eratosthenes (primes, sieveRange, primeList, psieveFrom, primeSieve)
 import Math.NumberTheory.Primes.Small
 import Math.NumberTheory.Primes.Types
diff --git a/Math/NumberTheory/Primes/Factorisation/Montgomery.hs b/Math/NumberTheory/Primes/Factorisation/Montgomery.hs
index aca75ffb1..5968c7f6b 100644
--- a/Math/NumberTheory/Primes/Factorisation/Montgomery.hs
+++ b/Math/NumberTheory/Primes/Factorisation/Montgomery.hs
@@ -35,7 +35,6 @@ module Math.NumberTheory.Primes.Factorisation.Montgomery
 
 import Control.Arrow
 import Control.Monad.Trans.State.Lazy
-import Data.Array.Base (bounds, unsafeAt)
 import Data.Bits
 import Data.IntMap (IntMap)
 import qualified Data.IntMap as IM
@@ -56,12 +55,11 @@ import System.Random
 import Math.NumberTheory.Curves.Montgomery
 import Math.NumberTheory.Euclidean.Coprimes (splitIntoCoprimes, unCoprimes)
 import Math.NumberTheory.Logarithms (integerLogBase')
-import Math.NumberTheory.Roots
-import Math.NumberTheory.Primes.Sieve.Eratosthenes (PrimeSieve(..), psieveFrom)
-import Math.NumberTheory.Primes.Sieve.Indexing (toPrim)
+import Math.NumberTheory.Primes.Sieve.Atkin
 import Math.NumberTheory.Primes.Small
 import Math.NumberTheory.Primes.Testing.Probabilistic
-import Math.NumberTheory.Utils hiding (splitOff)
+import Math.NumberTheory.Roots
+import Math.NumberTheory.Utils (shiftToOddCount#, shiftToOddCountBigNat)
 import Math.NumberTheory.Utils.FromIntegral
 
 -- | @'factorise' n@ produces the prime factorisation of @n@. @'factorise' 0@ is
@@ -143,12 +141,12 @@ curveFactorisation primeBound primeTest prng seed mbdigs n
         fact :: Integer -> Int -> State g [(Integer, Word)]
         fact 1 _ = return mempty
         fact m digs = do
-          let (b1, b2, ct) = findParms digs
+          let (b1, b1Sieve, b2, ct) = findParms digs
           -- All factors (both @pfs@ and @cfs@), are pairwise coprime. This is
           -- because 'repFact' returns either a single factor, or output of 'workFact'.
           -- In its turn, 'workFact' returns either a single factor,
           -- or concats 'repFact's over coprime integers. Induction completes the proof.
-          Factors pfs cfs <- repFact m b1 b2 ct
+          Factors pfs cfs <- repFact m b1 b1Sieve b2 ct
           case cfs of
             [] -> return pfs
             _  -> do
@@ -156,23 +154,23 @@ curveFactorisation primeBound primeTest prng seed mbdigs n
                   map (second (* j)) <$> fact k (if null pfs then digs + 5 else digs)
               return $ mconcat (pfs : nfs)
 
-        repFact :: Integer -> Word -> Word -> Word -> State g Factors
-        repFact 1 _ _ _ = return mempty
-        repFact m b1 b2 count =
+        repFact :: Integer -> Word -> PrimeSieve -> Word -> Word -> State g Factors
+        repFact 1 _ _ _ _ = return mempty
+        repFact m b1 b1Sieve b2 count =
           case perfPw m of
-            (_, 1) -> workFact m b1 b2 count
+            (_, 1) -> workFact m b1 b1Sieve b2 count
             (b, e)
               | ptest b   -> return $ singlePrimeFactor b e
-              | otherwise -> modifyPowers (* e) <$> workFact b b1 b2 count
+              | otherwise -> modifyPowers (* e) <$> workFact b b1 b1Sieve b2 count
 
-        workFact :: Integer -> Word -> Word -> Word -> State g Factors
-        workFact 1 _ _ _ = return mempty
-        workFact m _ _ 0 = return $ singleCompositeFactor m 1
-        workFact m b1 b2 count = do
+        workFact :: Integer -> Word -> PrimeSieve -> Word -> Word -> State g Factors
+        workFact 1 _ _ _ _ = return mempty
+        workFact m _ _ _ 0 = return $ singleCompositeFactor m 1
+        workFact m b1 b1Sieve b2 count = do
           s <- rndR m
           case someNatVal (fromInteger m) of
-            SomeNat (_ :: Proxy t) -> case montgomeryFactorisation b1 b2 (fromInteger s :: Mod t) of
-              Nothing -> workFact m b1 b2 (count - 1)
+            SomeNat (_ :: Proxy t) -> case montgomeryFactorisation b1 b1Sieve b2 (fromInteger s :: Mod t) of
+              Nothing -> workFact m b1 b1Sieve b2 (count - 1)
               Just d  -> do
                 let cs = unCoprimes $ splitIntoCoprimes [(d, 1), (m `quot` d, 1)]
                 -- Since all @cs@ are coprime, we can factor each of
@@ -181,7 +179,7 @@ curveFactorisation primeBound primeTest prng seed mbdigs n
                 fmap mconcat $ forM cs $
                   \(x, xm) -> if ptest x
                               then pure $ singlePrimeFactor x xm
-                              else repFact x b1 b2 (count - 1)
+                              else repFact x b1 b1Sieve b2 (count - 1)
 
 data Factors = Factors
   { _primeFactors     :: [(Integer, Word)]
@@ -269,8 +267,8 @@ badPower mx n
 --   It is assumed that @n@ has no small prime factors.
 --
 --   The result is maybe a nontrivial divisor of @n@.
-montgomeryFactorisation :: KnownNat n => Word -> Word -> Mod n -> Maybe Integer
-montgomeryFactorisation b1 b2 s = case newPoint (toInteger (unMod s)) n of
+montgomeryFactorisation :: KnownNat n => Word -> PrimeSieve -> Word -> Mod n -> Maybe Integer
+montgomeryFactorisation b1 b1Sieve b2 s = case newPoint (toInteger (unMod s)) n of
   Nothing             -> Nothing
   Just (SomePoint p0) -> do
     -- Small step: for each prime p <= b1
@@ -286,9 +284,7 @@ montgomeryFactorisation b1 b2 s = case newPoint (toInteger (unMod s)) n of
       g -> Just g
   where
     n = toInteger (natVal s)
-    smallPowers
-      = map findPower
-      $ takeWhile (<= b1) (2 : 3 : 5 : list primeStore)
+    smallPowers = map (findPower . fromIntegral) (atkinPrimeList b1Sieve)
     findPower p = go p
       where
         go acc
@@ -311,7 +307,7 @@ bigStep q b1 b2 = rs
         us * (pointZ p * pointX pq - pointX p * pointZ pq) `rem` n
         ) ts qks) 1 qs
 
-wheel :: Word
+wheel :: Num a => a
 wheel = 210
 
 wheelCoprimes :: [Word]
@@ -336,15 +332,6 @@ enumAndMultiplyFromThenTo p from thn to = zip [from, thn .. to] progression
 
     progression = pFrom : pThen : zipWith (`add` pStep) progression (tail progression)
 
--- primes, compactly stored as a bit sieve
-primeStore :: [PrimeSieve]
-primeStore = psieveFrom 7
-
--- generate list of primes from arrays
-list :: [PrimeSieve] -> [Word]
-list sieves = concat [[off + toPrim i | i <- [0 .. li], unsafeAt bs i]
-                                | PS vO bs <- sieves, let { (_,li) = bounds bs; off = fromInteger vO; }]
-
 -- | @'smallFactors' n@ finds all prime divisors of @n > 1@ up to 2^16 by trial division and returns the
 --   list of these together with their multiplicities, and a possible remaining factor which may be composite.
 smallFactors :: Natural -> ([(Natural, Word)], Maybe Natural)
@@ -403,8 +390,10 @@ smallFactors = \case
 -- ("tier") return parameters B1, B2 and the number of curves to try
 -- before next "tier".
 -- Roughly based on http://www.mersennewiki.org/index.php/Elliptic_Curve_Method#Choosing_the_best_parameters_for_ECM
-testParms :: IntMap (Word, Word, Word)
-testParms = IM.fromList
+testParms :: IntMap (Word, PrimeSieve, Word, Word)
+testParms
+  = IM.fromList
+  $ map (\(k, (b1, b2, ct)) -> (k, (b1, atkinSieve 0 (fromIntegral b1), b2, ct)))
   [ (12, (       400,        40000,     10))
   , (15, (      2000,       200000,     25))
   , (20, (     11000,      1100000,     90))
@@ -420,5 +409,7 @@ testParms = IM.fromList
   , (70, (2900000000, 290000000000, 340000))
   ]
 
-findParms :: Int -> (Word, Word, Word)
-findParms digs = maybe (wheel, 1000, 7) snd (IM.lookupLT digs testParms)
+findParms :: Int -> (Word, PrimeSieve, Word, Word)
+findParms digs
+  = maybe (wheel, atkinSieve 0 wheel, 1000, 7) snd
+  $ IM.lookupLT digs testParms
diff --git a/Math/NumberTheory/Primes/Factorisation/TrialDivision.hs b/Math/NumberTheory/Primes/Factorisation/TrialDivision.hs
index d6616cace..b4e178d05 100644
--- a/Math/NumberTheory/Primes/Factorisation/TrialDivision.hs
+++ b/Math/NumberTheory/Primes/Factorisation/TrialDivision.hs
@@ -16,8 +16,8 @@ module Math.NumberTheory.Primes.Factorisation.TrialDivision
     , trialDivisionPrimeTo
     ) where
 
-import Math.NumberTheory.Primes.Sieve.Eratosthenes (primeList, primeSieve, psieveList)
 import Math.NumberTheory.Roots
+import Math.NumberTheory.Primes.Sieve.Atkin
 import Math.NumberTheory.Primes.Types
 import Math.NumberTheory.Utils
 
@@ -32,7 +32,7 @@ trialDivisionWith prs n
     | otherwise = go n (integerSquareRoot n) prs
       where
         go !m !r (p:ps)
-            | r < p     = [(m,1)]
+            | r < p = [(m,1)]
             | otherwise =
                 case splitOff p m of
                   (0,_) -> go m r ps
@@ -45,11 +45,10 @@ trialDivisionWith prs n
 --   primes @<= bound@. If @n@ has prime divisors @> bound@, the last entry
 --   in the list is the product of all these. If @n <= bound^2@, this is a
 --   full factorisation, but very slow if @n@ has large prime divisors.
-trialDivisionTo :: Integer -> Integer -> [(Integer, Word)]
+trialDivisionTo :: Int -> Integer -> [(Integer, Word)]
 trialDivisionTo bd
-    | bd < 100      = trialDivisionTo 100
-    | bd < 10000000 = trialDivisionWith (map unPrime $ primeList $ primeSieve bd)
-    | otherwise     = trialDivisionWith (takeWhile (<= bd) $ map unPrime $ psieveList >>= primeList)
+    | bd < 100  = trialDivisionTo 100
+    | otherwise = trialDivisionWith (map (toInteger . unPrime) (atkinFromTo 0 bd))
 
 -- | Check whether a number is coprime to all of the numbers in the list
 --   (assuming that list contains only numbers > 1 and is ascending).
@@ -64,8 +63,7 @@ trialDivisionPrimeWith prs n
 
 -- | @'trialDivisionPrimeTo' bound n@ tests whether @n@ is coprime to all primes @<= bound@.
 --   If @n <= bound^2@, this is a full prime test, but very slow if @n@ has no small prime divisors.
-trialDivisionPrimeTo :: Integer -> Integer -> Bool
+trialDivisionPrimeTo :: Int -> Integer -> Bool
 trialDivisionPrimeTo bd
     | bd < 100      = trialDivisionPrimeTo 100
-    | bd < 10000000 = trialDivisionPrimeWith (map unPrime $ primeList $ primeSieve bd)
-    | otherwise     = trialDivisionPrimeWith (takeWhile (<= bd) $ map unPrime $ psieveList >>= primeList)
+    | otherwise     = trialDivisionPrimeWith (map (toInteger . unPrime) (atkinFromTo 0 bd))
diff --git a/Math/NumberTheory/Primes/Sieve/Atkin.hs b/Math/NumberTheory/Primes/Sieve/Atkin.hs
new file mode 100644
index 000000000..88d260d7b
--- /dev/null
+++ b/Math/NumberTheory/Primes/Sieve/Atkin.hs
@@ -0,0 +1,288 @@
+-- |
+-- Module:      Math.NumberTheory.Primes.Sieve.Atkin
+-- Copyright:   (c) 2019 Andrew Lelechenko
+-- Licence:     MIT
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+--
+-- Atkin sieve.
+--
+
+{-# LANGUAGE BangPatterns #-}
+
+module Math.NumberTheory.Primes.Sieve.Atkin
+  ( atkinPrimeList
+  , atkinSieve
+  , atkinFromTo
+  , PrimeSieve
+  ) where
+
+import Control.Monad
+import Control.Monad.ST
+import Data.Bit
+import Data.Bits
+import Data.Coerce
+import Data.Maybe
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Mutable as MU
+import Data.Word
+
+import Math.NumberTheory.Primes.Small
+import Math.NumberTheory.Primes.Types
+import Math.NumberTheory.Roots
+import Math.NumberTheory.Utils
+
+atkinFromTo :: Int -> Int -> [Prime Int]
+atkinFromTo low high = coerce $ atkinPrimeList $ atkinSieve low (high - low + 1)
+
+atkinPrimeList :: PrimeSieve -> [Int]
+atkinPrimeList (PrimeSieve low len segments)
+  | len60 == 0 = []
+  | otherwise = takeWhile (< high) $ dropWhile (< low) $ 2 : 3 : 5 : map ((+ low60 * 60) . fromWheel30) (listBits segments)
+  where
+    low60 = low `quot` 60
+    len60 = (low + len + 59) `quot` 60 - low60
+    high = low + len
+
+data PrimeSieve = PrimeSieve
+  { _psLowBound :: !Int
+  , _psLength   :: !Int
+  , _psSegments :: !(U.Vector Bit)
+  } deriving (Show)
+
+atkinSieve
+  :: Int
+  -> Int
+  -> PrimeSieve
+atkinSieve low len = PrimeSieve low len segments
+  where
+    low60 = low `quot` 60
+    len60 = (low + len + 59) `quot` 60 - low60
+    segments = sieveSegment low60 len60
+
+data SieveParams = SieveParams
+  { spDelta16  :: !Int -- [0..15]
+  , spDelta60  :: !Int -- [0..59]
+  -- ^ spDelta30 = fromWheel30 spDelta16
+  , spLowBound :: !Int
+  , spLength   :: !Int
+  } deriving (Show)
+
+-- | Take an index of the lowest period and number of periods.
+-- Return a bitmap for wheeled primes.
+sieveSegment
+  :: Int
+  -> Int
+  -> U.Vector Bit
+sieveSegment low60 len60 = runST $ do
+  vec <- MU.new (len60 `shiftL` 4)
+  forM_ fgs $ uncurry $ \delta16 -> traverse $
+    traverseLatticePoints (SieveParams delta16 (fromWheel30 delta16) low60 len60) vec
+  algo3steps456 low60 len60 vec
+  U.unsafeFreeze vec
+
+-- | Solutions of k * f^2 + l * g^2 = delta (mod 60)
+-- where (k, l) = (4, 1) for delta = 1 (mod 4)
+--              = (3, 1) for delta = 1 (mod 6)
+--              = (3,-1) for delta =11 (mod 12)
+-- Outer vector is indexed by delta=[0..15],
+-- inner vector is a list of pairs.
+fgs :: [(Int, [(Int, Int)])]
+fgs = map (\i -> (i, dispatch (fromWheel30 i))) [0..15]
+  where
+    dispatch delta
+      | delta `mod` 4 == 1
+      = [ (f, g) | f <- [1..15], g <- [1..30], (4*f*f + g*g - delta) `rem` 60 == 0]
+      | delta `mod` 6 == 1
+      = [ (f, g) | f <- [1..10], g <- [1..30], (3*f*f + g*g - delta) `rem` 60 == 0]
+      | delta `mod` 12 == 11
+      = [ (f, g) | f <- [1..10], g <- [1..30], (3*f*f - g*g - delta) `rem` 60 == 0]
+      | otherwise
+      = error "fgs: unexpected delta"
+
+traverseLatticePoints
+  :: SieveParams
+  -> MU.MVector s Bit
+  -> (Int, Int)
+  -> ST s ()
+traverseLatticePoints sp vec (x0, y0)
+  | spDelta60 sp `mod` 4 == 1
+  = traverseLatticePoints1 sp vec (x0, y0)
+  | spDelta60 sp `mod` 6 == 1
+  = traverseLatticePoints2 sp vec (x0, y0)
+  | spDelta60 sp `mod` 12 == 11
+  = traverseLatticePoints3 sp vec (x0, y0)
+  | otherwise
+  = error "traverseLatticePoints: unexpected delta"
+
+traverseLatticePoints1
+  :: SieveParams
+  -> MU.MVector s Bit
+  -> (Int, Int)
+  -> ST s ()
+traverseLatticePoints1 !sp vec (!x0, !y0) =
+  go kMax xMax y0
+  where
+    forwardY  (k, y) = (k +     y + 15, y + 30)
+    forwardX  (k, x) = (k + 2 * x + 15, x + 15)
+    backwardX (k, x) = (k - 2 * x + 15, x - 15)
+
+    -- Step 1
+    k0 = (4 * x0 * x0 + y0 * y0 - spDelta60 sp) `quot` 60 - spLowBound sp
+
+    -- Step 2
+    (kMax, xMax)
+      = backwardX
+      $ head
+      $ dropWhile (\(k, _) -> k < spLength sp)
+      $ iterate forwardX
+        (k0, x0)
+
+    -- Step 4
+    adjustY (!k, !y)
+      | k >= 0
+      = (k, y)
+      | otherwise
+      = adjustY $ forwardY (k, y)
+
+    -- Step 6
+    doActions (!k, !y)
+      | k < spLength sp
+      = unsafeFlipBit vec (k `shiftL` 4 + spDelta16 sp)
+        >> doActions (forwardY (k, y))
+      | otherwise
+      = pure ()
+
+    go !k !x !y
+      | x <= 0 = pure ()
+      | otherwise = do
+        let (k', y') = adjustY (k, y)
+        doActions (k', y')
+        let (k'', x') = backwardX (k', x)
+        go k'' x' y'
+
+traverseLatticePoints2
+  :: SieveParams
+  -> MU.MVector s Bit
+  -> (Int, Int)
+  -> ST s ()
+traverseLatticePoints2 sp vec (x0, y0) =
+  go kMax xMax y0
+  where
+    forwardY  (k, y) = (k + y + 15, y + 30)
+    forwardX  (k, x) = (k + x +  5, x + 10)
+    backwardX (k, x) = (k - x +  5, x - 10)
+
+    -- Step 1
+    k0 = (3 * x0 * x0 + y0 * y0 - spDelta60 sp) `quot` 60 - spLowBound sp
+
+    -- Step 2
+    (kMax, xMax)
+      = backwardX
+      $ head
+      $ dropWhile (\(k, _) -> k < spLength sp)
+      $ iterate forwardX
+        (k0, x0)
+
+    -- Step 4
+    adjustY (!k, !y)
+      | k >= 0
+      = (k, y)
+      | otherwise
+      = adjustY $ forwardY (k, y)
+
+    -- Step 6
+    doActions (!k, !y)
+      | k < spLength sp
+      = unsafeFlipBit vec (k `shiftL` 4 + spDelta16 sp)
+        >> doActions (forwardY (k, y))
+      | otherwise
+      = pure ()
+
+    go !k !x !y
+      | x <= 0 = pure ()
+      | otherwise = do
+        let (k', y') = adjustY (k, y)
+        doActions (k', y')
+        let (k'', x') = backwardX (k', x)
+        go k'' x' y'
+
+traverseLatticePoints3
+  :: SieveParams
+  -> MU.MVector s Bit
+  -> (Int, Int)
+  -> ST s ()
+traverseLatticePoints3 sp vec (x0, y0) =
+  go k0 x0 y0
+  where
+    forwardY  (k, y) = (k - y - 15, y + 30)
+    forwardX  (k, x) = (k + x +  5, x + 10)
+
+    -- Step 1
+    k0 = (3 * x0 * x0 - y0 * y0 - spDelta60 sp) `quot` 60 - spLowBound sp
+
+    -- Step 6
+    doActions (!k, !x, !y)
+      | k >= 0 && y < x
+      = unsafeFlipBit vec (k `shiftL` 4 + spDelta16 sp)
+        >> (let (k', y') = forwardY (k, y) in doActions (k', x, y'))
+      | otherwise
+      = pure ()
+
+    go !k !x !y
+      | k >= spLength sp
+      , x <= y
+      = pure ()
+      | k >= spLength sp
+      = let (k', y') = forwardY (k, y) in
+        go k' x y'
+      | otherwise
+      = do
+        doActions (k, x, y)
+        let (k', x') = forwardX (k, x)
+        go k' x' y
+
+-- | Perform steps 4-6 of Algorithm 3.X.
+algo3steps456
+  :: Int
+  -> Int
+  -> MU.MVector s Bit
+  -> ST s ()
+algo3steps456 _ 0 _ = pure ()
+algo3steps456 low60 len60 vec =
+  forM_ ps $ \p ->
+    crossMultiples low60 len60 vec (p * p)
+  where
+    low  = 7
+    high = integerSquareRoot (60 * (low60 + len60) - 1)
+    ps = case toIntegralSized high of
+      Just high' -> map fromIntegral (smallPrimesFromTo low high')
+      Nothing    -> map fromIntegral (smallPrimesFromTo low maxBound) ++ atkinPrimeList (atkinSieve low' (high - low' + 1))
+        where
+          low' = fromIntegral (maxBound :: Word16) + 1
+
+-- | Cross out multiples of the first argument
+-- in a given sieve.
+crossMultiples
+  :: Int
+  -> Int
+  -> MU.MVector s Bit
+  -> Int -- coprime with 60
+  -> ST s ()
+crossMultiples low60 len60 vec m =
+  forM_ [0..15] $ \i -> do
+        -- k0 is the smallest non-negative k such that 60k+delta = 0 (mod m)
+    let k0 = solveCongruence (fromWheel30 i) m
+        -- k1 = k0 (mod m), k1 >= lowBound
+        k1 = if r <= k0 then q * m + k0 else (q + 1) * m + k0
+    forM_ [k1, k1 + m .. (low60 + len60) - 1] $
+      \k -> MU.unsafeWrite vec ((k - low60) `shiftL` 4 + i) (Bit False)
+  where
+    (q, r) = low60 `quotRem` m
+
+-- Find the smallest k such that 60k+delta = 0 (mod m)
+-- Should be equal to
+-- (`quot` 60) $ fromJust $ chineseCoprime (delta, 60) (0, m)
+solveCongruence :: Int -> Int -> Int
+solveCongruence delta m = (- recip60 * delta) `mod` m
+  where
+    recip60 = fromInteger $ fromJust $ recipMod 60 (toInteger m)
diff --git a/Math/NumberTheory/Primes/Small.hs b/Math/NumberTheory/Primes/Small.hs
index 66c2f6ba1..3ba94d709 100644
--- a/Math/NumberTheory/Primes/Small.hs
+++ b/Math/NumberTheory/Primes/Small.hs
@@ -42,6 +42,7 @@ smallPrimesFromTo (W16# from#) (W16# to#) = go k0#
       = W16# p# : go (k# +# 1#)
       where
         p# = indexWord16OffAddr# smallPrimesAddr# k#
+{-# INLINE smallPrimesFromTo #-}
 
 -- length smallPrimes
 smallPrimesLength :: Int
diff --git a/Math/NumberTheory/Primes/Testing/Certified.hs b/Math/NumberTheory/Primes/Testing/Certified.hs
index e5a40b854..043af5a16 100644
--- a/Math/NumberTheory/Primes/Testing/Certified.hs
+++ b/Math/NumberTheory/Primes/Testing/Certified.hs
@@ -25,7 +25,7 @@ import Math.NumberTheory.Primes (unPrime)
 import Math.NumberTheory.Primes.Factorisation.TrialDivision (trialDivisionPrimeTo, trialDivisionTo, trialDivisionWith)
 import Math.NumberTheory.Primes.Factorisation.Montgomery (montgomeryFactorisation, smallFactors, findParms)
 import Math.NumberTheory.Primes.Testing.Probabilistic (bailliePSW, isPrime, isStrongFermatPP, lucasTest)
-import Math.NumberTheory.Primes.Sieve.Eratosthenes (primeList, primeSieve)
+import Math.NumberTheory.Primes.Sieve.Atkin
 import Math.NumberTheory.Utils (splitOff)
 
 -- | @'isCertifiedPrime' n@ tests primality of @n@, first trial division
@@ -51,7 +51,7 @@ data PrimalityProof
                   , _knownFactors :: ![(Integer, Word, Integer, PrimalityProof)]
                   }
     | TrialDivision { cprime :: !Integer        -- ^ The number whose primality is proved.
-                    , _tdLimit :: !Integer }
+                    , tdLimit :: !Int }
     | Trivial { cprime :: !Integer              -- ^ The number whose primality is proved.
               }
       deriving Show
@@ -60,9 +60,12 @@ data PrimalityProof
 --   impossible to create invalid proofs by the public interface, this
 --   should never return 'False'.
 checkPrimalityProof :: PrimalityProof -> Bool
-checkPrimalityProof (Trivial n) = isTrivialPrime n
-checkPrimalityProof (TrialDivision p b) = p <= b*b && trialDivisionPrimeTo b p
-checkPrimalityProof (Pocklington p a b fcts) = b > 0 && a > b && a*b == pm1 && a == ppProd fcts && all verify fcts
+checkPrimalityProof (Trivial n) =
+  isTrivialPrime n
+checkPrimalityProof (TrialDivision p b) =
+  p <= toInteger b * toInteger b && trialDivisionPrimeTo b p
+checkPrimalityProof (Pocklington p a b fcts) =
+  b > 0 && a > b && a*b == pm1 && a == ppProd fcts && all verify fcts
   where
     pm1 = p-1
     ppProd pps = product [pf^e | (pf,e,_,_) <- pps]
@@ -87,7 +90,7 @@ trivialPrimes = [2,3,5,7,11,13,17,19,23,29]
 smallCert :: Integer -> PrimalityProof
 smallCert n
     | n < 30    = Trivial n
-    | otherwise = TrialDivision n (integerSquareRoot n + 1)
+    | otherwise = TrialDivision n (fromInteger (integerSquareRoot n + 1))
 
 -- | @'certify' n@ constructs, for @n > 1@, a proof of either
 --   primality or compositeness of @n@. This may take a very long
@@ -99,7 +102,7 @@ certify n
                     ((p,_):_) | p < n     -> Nothing
                               | otherwise -> Just (Trivial n)
                     _ -> error "Impossible"
-    | n < billi = let r2 = integerSquareRoot n + 2 in
+    | n < billi = let r2 = fromInteger (integerSquareRoot n + 2) in
                   case trialDivisionTo r2 n of
                     ((p,_):_) | p < n       -> Nothing
                               | otherwise   -> Just (TrialDivision n r2)
@@ -109,7 +112,7 @@ certify n
                                  | not (lucasTest n) -> Nothing
                                  | otherwise -> Just (certifyBPSW n)       -- if it isn't we error and ask for a report.
                     ((toInteger -> p,_):_, _)
-                      | p == n -> Just (TrialDivision n (min 100000 n))
+                      | p == n -> Just (TrialDivision n (fromInteger (min 100000 n)))
                       | otherwise -> Nothing
                     _ -> error ("***Error factorising " ++ show n ++ "! Please report this to maintainer of arithmoi.")
       where
@@ -154,8 +157,8 @@ findDecomposition :: Integer -> (Integer, [(Integer, Word, Bool)], Integer)
 findDecomposition n = go 1 n [] prms
   where
     sr = integerSquareRoot n
-    pbd = min 1000000 (sr+20)
-    prms = map unPrime $ primeList (primeSieve pbd)
+    pbd = fromInteger (min 1000000 (sr+20))
+    prms = map (toInteger . unPrime) (atkinFromTo 0 pbd)
     go a b afs (p:ps)
         | a > b     = (a,afs,b)
         | otherwise = case splitOff p b of
@@ -174,20 +177,20 @@ findDecomposition n = go 1 n [] prms
 --   starting with curve s. Actually, this may loop infinitely, but the
 --   loop should not be entered before the heat death of the universe.
 findFactor :: Integer -> Int -> Integer -> Integer
-findFactor n digits s = case findLoop n lo hi count s of
+findFactor n digits s = case findLoop n lo loSieve hi count s of
                           Left t  -> findFactor n (digits+5) t
                           Right f -> f
   where
-    (lo,hi,count) = findParms digits
+    (lo, loSieve, hi, count) = findParms digits
 
 -- | Find a factor or say with which curve to continue.
-findLoop :: Integer -> Word -> Word -> Word -> Integer -> Either Integer Integer
-findLoop _ _  _  0  s = Left s
-findLoop n lo hi ct s
+findLoop :: Integer -> Word -> PrimeSieve -> Word -> Word -> Integer -> Either Integer Integer
+findLoop _ _ _ _ 0 s = Left s
+findLoop n lo loSieve hi ct s
     | n <= s+2  = Left 6
     | otherwise = case someNatVal (fromInteger n) of
-                    SomeNat (_ :: Proxy t) -> case montgomeryFactorisation lo hi (fromInteger s :: Mod t) of
-                      Nothing -> findLoop n lo hi (ct-1) (s+1)
+                    SomeNat (_ :: Proxy t) -> case montgomeryFactorisation lo loSieve hi (fromInteger s :: Mod t) of
+                      Nothing -> findLoop n lo loSieve hi (ct-1) (s+1)
                       Just fct
                         | bailliePSW fct -> Right fct
                         | otherwise -> Right (findFactor fct 8 (s+1))
diff --git a/Math/NumberTheory/Utils.hs b/Math/NumberTheory/Utils.hs
index d951a4672..ee60bce7b 100644
--- a/Math/NumberTheory/Utils.hs
+++ b/Math/NumberTheory/Utils.hs
@@ -30,6 +30,7 @@ module Math.NumberTheory.Utils
     , recipMod
 
     , toWheel30
+    , toWheel30Int
     , fromWheel30
     , withSomeKnown
     , intVal
@@ -202,16 +203,45 @@ bigNatToNatural bn
 
 -------------------------------------------------------------------------------
 -- Helpers for mapping to rough numbers and back.
--- Copypasted from Data.BitStream.WheelMapping
+-- Copypasted from Data.Chimera.WheelMapping
 
+bits :: Int
+bits = finiteBitSize (0 :: Word)
+
+-- | Left inverse for 'fromWheel30'. Monotonically non-decreasing function.
+--
+-- prop> toWheel30 . fromWheel30 == id
 toWheel30 :: (Integral a, Bits a) => a -> a
 toWheel30 i = q `shiftL` 3 + (r + r `shiftR` 4) `shiftR` 2
   where
     (q, r) = i `P.quotRem` 30
+{-# INLINE toWheel30 #-}
+
+toWheel30Int :: Int -> Int
+toWheel30Int i@(I# i#) = q `shiftL` 3 + (r + r `shiftR` 4) `shiftR` 2
+  where
+    (q, r) = case bits of
+      64 -> (q64, r64)
+      _  -> i `quotRem` 30
 
+    m# = 9838263505978427529## -- (2^67+7) / 15
+    !(# z1#, _ #) = timesWord2# m# (int2Word# i#)
+    q64 = I# (word2Int# z1#) `shiftR` 4
+    r64 = i - q64 `shiftL` 5 + q64 `shiftL` 1
+
+{-# INLINE toWheel30Int #-}
+
+-- | 'fromWheel30' n is the (n+1)-th positive number, not divisible by 2, 3 or 5.
+-- Sequence <https://oeis.org/A007775 A007775>.
+--
+-- prop> map fromWheel30 [0..] == [ n | n <- [0..], n `gcd` 30 == 1 ]
+--
+-- > > map fromWheel30 [0..9]
+-- > [1,7,11,13,17,19,23,29,31,37]
 fromWheel30 :: (Num a, Bits a) => a -> a
 fromWheel30 i = ((i `shiftL` 2 - i `shiftR` 2) .|. 1)
               + ((i `shiftL` 1 - i `shiftR` 1) .&. 2)
+{-# INLINE fromWheel30 #-}
 
 -------------------------------------------------------------------------------
 -- Helpers for dealing with data types parametrised by natural numbers.
diff --git a/app/Atkin.hs b/app/Atkin.hs
new file mode 100644
index 000000000..f7b46b069
--- /dev/null
+++ b/app/Atkin.hs
@@ -0,0 +1,9 @@
+module Main where
+
+import Math.NumberTheory.Primes
+
+atkin :: (Int, Int) -> Int
+atkin (p, q) = sum $ atkinPrimeList $ atkinSieve p q
+
+main :: IO ()
+main = print $ atkin (10000000000,100000000)
diff --git a/arithmoi.cabal b/arithmoi.cabal
index 6b7e2d9b5..0abafcdad 100644
--- a/arithmoi.cabal
+++ b/arithmoi.cabal
@@ -30,7 +30,7 @@ library
   build-depends:
     base >=4.10 && <5,
     array >=0.5 && <0.6,
-    bitvec,
+    bitvec >=0.2,
     containers >=0.5.8 && <0.7,
     chimera >=0.3,
     constraints,
@@ -90,6 +90,7 @@ library
     Math.NumberTheory.Primes.Counting.Impl
     Math.NumberTheory.Primes.Factorisation.Montgomery
     Math.NumberTheory.Primes.Factorisation.TrialDivision
+    Math.NumberTheory.Primes.Sieve.Atkin
     Math.NumberTheory.Primes.Sieve.Eratosthenes
     Math.NumberTheory.Primes.Sieve.Indexing
     Math.NumberTheory.Primes.Small
@@ -237,3 +238,14 @@ executable arithmoi-quadratic-sieve
   default-language: Haskell2010
   hs-source-dirs: app
   ghc-options: -Wall -Widentities -Wcompat -rtsopts
+
+executable atkin
+  build-depends:
+    base,
+    arithmoi,
+    containers
+  buildable: True
+  main-is: Atkin.hs
+  hs-source-dirs: app
+  default-language: Haskell2010
+  ghc-options: -Wall -Widentities -Wcompat
diff --git a/benchmark/Math/NumberTheory/SequenceBench.hs b/benchmark/Math/NumberTheory/SequenceBench.hs
index aee78fc15..1c92765d4 100644
--- a/benchmark/Math/NumberTheory/SequenceBench.hs
+++ b/benchmark/Math/NumberTheory/SequenceBench.hs
@@ -1,3 +1,5 @@
+{-# LANGUAGE TupleSections #-}
+
 {-# OPTIONS_GHC -fno-warn-type-defaults #-}
 
 module Math.NumberTheory.SequenceBench
@@ -8,35 +10,38 @@ import Gauge.Main
 
 import Data.Array.Unboxed
 import Data.Bits
+import Data.Maybe
 
-import Math.NumberTheory.Primes (Prime(..), nextPrime, precPrime)
-import Math.NumberTheory.Primes.Testing
+import Math.NumberTheory.Primes
 
 filterIsPrime :: (Integer, Integer) -> Integer
-filterIsPrime (p, q) = sum $ takeWhile (<= q) $ dropWhile (< p) $ filter isPrime (map toPrim [toIdx p .. toIdx q])
+filterIsPrime (p, q) = sum $ takeWhile (<= p + q) $ dropWhile (< p) $ filter (isJust . isPrime) (map toPrim [toIdx p .. toIdx (p + q)])
 
 eratosthenes :: (Integer, Integer) -> Integer
-eratosthenes (p, q) = sum (map unPrime [nextPrime p .. precPrime q])
+eratosthenes (p, q) = sum (map unPrime [nextPrime p .. precPrime (p + q)])
+
+atkin :: (Integer, Integer) -> Integer
+atkin (p, q) = toInteger $ sum $ atkinPrimeList $ atkinSieve (fromInteger p) (fromInteger q)
 
 filterIsPrimeBench :: Benchmark
 filterIsPrimeBench = bgroup "filterIsPrime" $
-  [ bench (show (10^x, 10^y)) $ nf filterIsPrime (10^x, 10^x + 10^y)
+  [ bench (show (10^x, 10^y)) $ nf filterIsPrime (10^x, 10^x)
   | x <- [5..8]
   , y <- [3..x-1]
   ]
 
-eratosthenesBench :: Benchmark
-eratosthenesBench = bgroup "eratosthenes" $
-  [ bench (show (10^x, 10^y)) $ nf eratosthenes (10^x, 10^x + 10^y)
-  | x <- [10..17]
-  , y <- [6..x-1]
-  , x == 10 || y == 7
+sieveBench :: Benchmark
+sieveBench = bgroup "sieve" $ concat
+  [ [ bench ("eratosthenes/" ++ show (10^x, 10^y)) $ nf eratosthenes (10^x, 10^y)
+    , bench ("atkin/"        ++ show (10^x, 10^y)) $ nf atkin        (10^x, 10^y)
+    ]
+  | (x, y) <- map (10,) [6..9] ++ map (,7) [10..16]
   ]
 
 benchSuite :: Benchmark
 benchSuite = bgroup "Sequence"
-    [ filterIsPrimeBench
-    , eratosthenesBench
+    [ sieveBench
+    , filterIsPrimeBench
     ]
 
 -------------------------------------------------------------------------------
diff --git a/cabal.project b/cabal.project
new file mode 100644
index 000000000..5a9152e11
--- /dev/null
+++ b/cabal.project
@@ -0,0 +1,7 @@
+packages: .
+tests: True
+benchmarks: True
+constraints:
+  optparse-applicative -process,
+  tasty -unix,
+  quickcheck-classes -aeson -semigroupois
diff --git a/test-suite/Math/NumberTheory/Primes/SieveTests.hs b/test-suite/Math/NumberTheory/Primes/SieveTests.hs
index 4e9776452..b6b6931a3 100644
--- a/test-suite/Math/NumberTheory/Primes/SieveTests.hs
+++ b/test-suite/Math/NumberTheory/Primes/SieveTests.hs
@@ -21,15 +21,16 @@ import Prelude hiding (words)
 
 import Test.Tasty
 import Test.Tasty.HUnit
+import Test.Tasty.QuickCheck
 
 import Data.Bits
 import Data.Int
+import Data.Maybe
 import Data.Proxy (Proxy(..))
 import Data.Word
 import Numeric.Natural (Natural)
 
-import Math.NumberTheory.Primes (Prime, unPrime, primes, nextPrime, precPrime, UniqueFactorisation)
-import Math.NumberTheory.Primes.Testing
+import Math.NumberTheory.Primes
 import Math.NumberTheory.TestUtils
 
 lim1 :: Num a => a
@@ -42,12 +43,21 @@ lim2 = 100000
 primesProperty1 :: forall a. (Integral a, Show a) => Proxy a -> Assertion
 primesProperty1 _ = assertEqual "primes matches isPrime"
   (takeWhile (<= lim1) (map unPrime primes) :: [a])
-  (filter (isPrime . toInteger) [1..lim1])
+  (filter (isJust . isPrime . toInteger) [1..lim1])
 
 primesProperty2 :: forall a. (Integral a, Bounded a, Show a) => Proxy a -> Assertion
 primesProperty2 _ = assertEqual "primes matches isPrime"
   (map unPrime primes :: [a])
-  (filter (isPrime . toInteger) [1..maxBound])
+  (filter (isJust . isPrime . toInteger) [1..maxBound])
+
+atkinPrimesProperty1 :: Large Word -> Large Word -> Property
+atkinPrimesProperty1 (Large x) (Large y) = actual === expected
+  where
+    lim = 1000000
+    from = fromIntegral $ x `mod` lim
+    to   = fromIntegral $ x `mod` lim + y `mod` lim
+    expected = mapMaybe isPrime [from..to]
+    actual   = atkinFromTo from to
 
 -- | Check that 'primeList' from 'primeSieve' matches truncated 'primes'.
 primeSieveProperty1 :: AnySign Integer -> Bool
@@ -66,11 +76,12 @@ psieveListProperty1 _ = assertEqual "primes == primeList . psieveList"
 psieveListProperty2 :: forall a. (Integral a, Bounded a, Show a) => Proxy a -> Assertion
 psieveListProperty2 _ = assertEqual "primes == primeList . psieveList"
   (map unPrime primes :: [a])
-  (filter (isPrime . toInteger) [0..maxBound])
+  (filter (isJust . isPrime . toInteger) [0..maxBound])
 
 testSuite :: TestTree
 testSuite = testGroup "Sieve"
-  [ testGroup "primes"
+  [ testProperty "atkinPrimes" atkinPrimesProperty1
+  , testGroup "primes"
     [ testCase "Int"     (primesProperty1 (Proxy :: Proxy Int))
     , testCase "Word"    (primesProperty1 (Proxy :: Proxy Word))
     , testCase "Integer" (primesProperty1 (Proxy :: Proxy Integer))