From 8aa7506b1dcea9121c2325e10d3c65961c658165 Mon Sep 17 00:00:00 2001
From: Bodigrim <andrew.lelechenko@gmail.com>
Date: Sun, 16 Jun 2019 19:01:15 +0100
Subject: [PATCH 01/18] Initial implementation of Atkin sieve

---
 Math/NumberTheory/Primes/Sieve/Atkin.hs | 273 ++++++++++++++++++++++++
 arithmoi.cabal                          |   3 +-
 2 files changed, 275 insertions(+), 1 deletion(-)
 create mode 100644 Math/NumberTheory/Primes/Sieve/Atkin.hs

diff --git a/Math/NumberTheory/Primes/Sieve/Atkin.hs b/Math/NumberTheory/Primes/Sieve/Atkin.hs
new file mode 100644
index 000000000..0ba3e70da
--- /dev/null
+++ b/Math/NumberTheory/Primes/Sieve/Atkin.hs
@@ -0,0 +1,273 @@
+-- |
+-- Module:      Math.NumberTheory.Primes.Sieve.Atkin
+-- Copyright:   (c) 2019 Andrew Lelechenko
+-- Licence:     MIT
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+--
+-- Atkin sieve.
+--
+
+module Math.NumberTheory.Primes.Sieve.Atkin
+  ( atkinPrimeList
+  , atkinSieve
+  ) where
+
+import Control.Monad
+import Control.Monad.ST
+import Data.Bit
+import Data.Foldable
+import Data.Maybe
+import qualified Data.Vector as V
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Mutable as MU
+
+import Math.NumberTheory.Moduli.Chinese
+import Math.NumberTheory.Roots
+import qualified Math.NumberTheory.Primes.Sieve.Eratosthenes as E
+import Math.NumberTheory.Primes.Types
+import Math.NumberTheory.Utils
+
+atkinPrimeList :: PrimeSieve -> [Int]
+atkinPrimeList (PrimeSieve low len segments)
+  | len60 == 0 = []
+  | len60 == 1 = takeWhile (< high) $ dropWhile (< low) $ 2 : 3 : 5 : doIx 0
+  | otherwise
+  =  dropWhile (< low) (2 : 3 : 5 : doIx 0)
+  ++ concatMap doIx [1 .. len60 - 2]
+  ++ takeWhile (< high) (doIx $ len60 - 1)
+  where
+    low60 = low `quot` 60
+    len60 = (low + len + 59) `quot` 60 - low60
+    high = low + len
+
+    js = map fromWheel30 [0..15]
+
+    doIx k
+      = map ((+ 60 * (low60 + k)) . fst)
+      $ filter (unBit . snd)
+      $ zip js
+      $ toList
+      $ V.map (U.! k) segments
+
+data PrimeSieve = PrimeSieve
+  { _psLowBound :: !Int
+  , _psLength   :: !Int
+  , _psSegments :: V.Vector (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
+    params = V.generate 16 (\i -> SieveParams (fromWheel30 i) low60 len60)
+    segments = V.map sieveSegment params
+
+data SieveParams = SieveParams
+  { spDelta    :: !Int
+  , spLowBound :: !Int
+  , spLength   :: !Int
+  } deriving (Show)
+
+spHighBound :: SieveParams -> Int
+spHighBound sp = spLowBound sp + spLength sp
+
+sieveSegment
+  :: SieveParams
+  -> U.Vector Bit
+sieveSegment sp = runST $ do
+  vec <- MU.new (spLength sp)
+  U.forM_ (fgs V.! toWheel30 (spDelta sp)) $
+    traverseLatticePoints sp (\k -> unsafeFlipBit vec (k - spLowBound sp))
+  algo3steps456 sp 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)
+fgs :: V.Vector (U.Vector (Int, Int))
+fgs = V.generate 16 (dispatch . fromWheel30)
+  where
+    dispatch delta
+      | delta `mod` 4 == 1
+      = U.fromList [ (f, g) | f <- [1..15], g <- [1..30], (4*f*f + g*g - delta) `rem` 60 == 0]
+      | delta `mod` 6 == 1
+      = U.fromList [ (f, g) | f <- [1..10], g <- [1..30], (3*f*f + g*g - delta) `rem` 60 == 0]
+      | delta `mod` 12 == 11
+      = U.fromList [ (f, g) | f <- [1..10], g <- [1..30], (3*f*f - g*g - delta) `rem` 60 == 0]
+      | otherwise
+      = error "fgs: unexpected delta"
+
+traverseLatticePoints
+  :: SieveParams
+  -> (Int -> ST s ())
+  -> (Int, Int)
+  -> ST s ()
+traverseLatticePoints sp action (x0, y0)
+  | spDelta sp `mod` 4 == 1
+  = traverseLatticePoints1 sp action (x0, y0)
+  | spDelta sp `mod` 6 == 1
+  = traverseLatticePoints2 sp action (x0, y0)
+  | spDelta sp `mod` 12 == 11
+  = traverseLatticePoints3 sp action (x0, y0)
+  | otherwise
+  = error "traverseLatticePoints: unexpected delta"
+
+traverseLatticePoints1
+  :: SieveParams
+  -> (Int -> ST s ())
+  -> (Int, Int)
+  -> ST s ()
+traverseLatticePoints1 sp action (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 - spDelta sp) `quot` 60
+
+    -- Step 2
+    (kMax, xMax)
+      = backwardX
+      $ head
+      $ dropWhile (\(k, _) -> k < spHighBound sp)
+      $ iterate forwardX
+      $ (k0, x0)
+
+    -- Step 4
+    adjustY
+      = head
+      . dropWhile (\(k, _) -> k < spLowBound sp)
+      . iterate forwardY
+
+    -- Step 6
+    doActions (k, y)
+      = traverse_ action
+      $ takeWhile (< spHighBound sp)
+      $ map fst
+      $ iterate forwardY
+      $ (k, y)
+
+    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
+  -> (Int -> ST s ())
+  -> (Int, Int)
+  -> ST s ()
+traverseLatticePoints2 sp action (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 - spDelta sp) `quot` 60
+
+    -- Step 2
+    (kMax, xMax)
+      = backwardX
+      $ head
+      $ dropWhile (\(k, _) -> k < spHighBound sp)
+      $ iterate forwardX
+      $ (k0, x0)
+
+    -- Step 4
+    adjustY
+      = head
+      . dropWhile (\(k, _) -> k < spLowBound sp)
+      . iterate forwardY
+
+    -- Step 6
+    doActions (k, y)
+      = traverse_ action
+      $ takeWhile (< spHighBound sp)
+      $ map fst
+      $ iterate forwardY
+      $ (k, y)
+
+    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
+  -> (Int -> ST s ())
+  -> (Int, Int)
+  -> ST s ()
+traverseLatticePoints3 sp action (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 - spDelta sp) `quot` 60
+
+    -- Step 6
+    doActions (k, x, y)
+      = traverse_ action
+      $ map fst
+      $ takeWhile (\(k', y') -> k' >= spLowBound sp && y' < x)
+      $ iterate forwardY
+      $ (k, y)
+
+    go k x y
+      | k >= spHighBound sp
+      , x <= y
+      = pure ()
+      | k >= spHighBound 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
+  :: SieveParams
+  -> MU.MVector s Bit
+  -> ST s ()
+algo3steps456 sp vec =
+  forM_ ps $ \p ->
+    crossMultiples sp vec (p * p)
+  where
+    low  = 7
+    high = integerSquareRoot (60 * spHighBound sp - 1)
+    ps   = takeWhile (<= high) $ dropWhile (< low) $ map unPrime E.primes
+
+-- | Cross out multiples of the first argument
+-- in a given sieve.
+crossMultiples
+  :: SieveParams
+  -> MU.MVector s Bit
+  -> Int -- coprime with 60
+  -> ST s ()
+crossMultiples sp vec m =
+  forM_ [k1, k1 + m .. spHighBound sp - 1] $
+    \k -> MU.unsafeWrite vec (k - spLowBound sp) (Bit False)
+  where
+    -- k0 is the smallest non-negative k such that 60k+delta = 0 (mod m)
+    k0 = (`quot` 60) $ (\k -> if k < 0 then k + 60 * m else k) $ fst $ fromJust $ chinese (spDelta sp, 60) (0, m)
+    -- k1 = k0 (mod m), k1 >= lowBound
+    (q, r) = spLowBound sp `quotRem` m
+    k1 = if r < k0 then q * m + k0 else (q + 1) * m + k0
diff --git a/arithmoi.cabal b/arithmoi.cabal
index 6b7e2d9b5..3abbc8dec 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

From a84fa5092dcd66d96a53c1e77f846a850d5e2ed9 Mon Sep 17 00:00:00 2001
From: Bodigrim <andrew.lelechenko@gmail.com>
Date: Fri, 26 Jul 2019 20:26:52 +0100
Subject: [PATCH 02/18] Benchmark and test Atkin sieve

---
 Math/NumberTheory/Primes.hs                   |  4 ++++
 benchmark/Math/NumberTheory/SequenceBench.hs  | 19 +++++++++++++++---
 .../Math/NumberTheory/Primes/SieveTests.hs    | 20 +++++++++++++------
 3 files changed, 34 insertions(+), 9 deletions(-)

diff --git a/Math/NumberTheory/Primes.hs b/Math/NumberTheory/Primes.hs
index a60f5bc3f..9c2530843 100644
--- a/Math/NumberTheory/Primes.hs
+++ b/Math/NumberTheory/Primes.hs
@@ -20,6 +20,9 @@ module Math.NumberTheory.Primes
     , factorBack
     , -- * Old interface
       primes
+    , -- * Temporary
+      atkinPrimeList
+    , atkinSieve
     ) where
 
 import Data.Bits
@@ -31,6 +34,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/benchmark/Math/NumberTheory/SequenceBench.hs b/benchmark/Math/NumberTheory/SequenceBench.hs
index aee78fc15..c843ef46c 100644
--- a/benchmark/Math/NumberTheory/SequenceBench.hs
+++ b/benchmark/Math/NumberTheory/SequenceBench.hs
@@ -8,16 +8,19 @@ 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 (<= q) $ dropWhile (< p) $ filter (isJust . isPrime) (map toPrim [toIdx p .. toIdx q])
 
 eratosthenes :: (Integer, Integer) -> Integer
 eratosthenes (p, q) = sum (map unPrime [nextPrime p .. precPrime q])
 
+atkin :: (Integer, Integer) -> Integer
+atkin (p, q) = toInteger $ sum $ atkinPrimeList $ atkinSieve (fromInteger p) (fromInteger $ q - p)
+
 filterIsPrimeBench :: Benchmark
 filterIsPrimeBench = bgroup "filterIsPrime" $
   [ bench (show (10^x, 10^y)) $ nf filterIsPrime (10^x, 10^x + 10^y)
@@ -33,10 +36,20 @@ eratosthenesBench = bgroup "eratosthenes" $
   , x == 10 || y == 7
   ]
 
+atkinBench :: Benchmark
+atkinBench = bgroup "atkin" $
+  map (\(x, y) -> bench (show (x, y)) $ nf atkin (x, x + y))
+  [ (10 ^ x, 10 ^ y)
+  | x <- [10..17]
+  , y <- [6..x-1]
+  , x == 10 || y == 7
+  ]
+
 benchSuite :: Benchmark
 benchSuite = bgroup "Sequence"
     [ filterIsPrimeBench
     , eratosthenesBench
+    , atkinBench
     ]
 
 -------------------------------------------------------------------------------
diff --git a/test-suite/Math/NumberTheory/Primes/SieveTests.hs b/test-suite/Math/NumberTheory/Primes/SieveTests.hs
index 4e9776452..3effa5ef5 100644
--- a/test-suite/Math/NumberTheory/Primes/SieveTests.hs
+++ b/test-suite/Math/NumberTheory/Primes/SieveTests.hs
@@ -24,12 +24,12 @@ import Test.Tasty.HUnit
 
 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 +42,19 @@ 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 :: Assertion
+atkinPrimesProperty1 = assertEqual "atkinPrimes matches isPrime" expected actual
+  where
+    lim = 1000000
+    expected = filter (isJust . isPrime . toInteger) [1..lim]
+    actual   = atkinPrimeList $ atkinSieve 1 lim
 
 -- | Check that 'primeList' from 'primeSieve' matches truncated 'primes'.
 primeSieveProperty1 :: AnySign Integer -> Bool
@@ -66,11 +73,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"
+  [ testCase "atkinPrimes" atkinPrimesProperty1
+  , testGroup "primes"
     [ testCase "Int"     (primesProperty1 (Proxy :: Proxy Int))
     , testCase "Word"    (primesProperty1 (Proxy :: Proxy Word))
     , testCase "Integer" (primesProperty1 (Proxy :: Proxy Integer))

From 76bf177aa1ebb85d720b476c2a791c54810d200d Mon Sep 17 00:00:00 2001
From: Bodigrim <andrew.lelechenko@gmail.com>
Date: Fri, 26 Jul 2019 20:27:28 +0100
Subject: [PATCH 03/18] Use manual recursion in traverseLatticePoints{1,2,3}

---
 Math/NumberTheory/Primes/Sieve/Atkin.hs | 67 +++++++++++++------------
 1 file changed, 34 insertions(+), 33 deletions(-)

diff --git a/Math/NumberTheory/Primes/Sieve/Atkin.hs b/Math/NumberTheory/Primes/Sieve/Atkin.hs
index 0ba3e70da..7cad670f5 100644
--- a/Math/NumberTheory/Primes/Sieve/Atkin.hs
+++ b/Math/NumberTheory/Primes/Sieve/Atkin.hs
@@ -7,6 +7,8 @@
 -- Atkin sieve.
 --
 
+{-# LANGUAGE BangPatterns #-}
+
 module Math.NumberTheory.Primes.Sieve.Atkin
   ( atkinPrimeList
   , atkinSieve
@@ -122,7 +124,7 @@ traverseLatticePoints1
   -> (Int -> ST s ())
   -> (Int, Int)
   -> ST s ()
-traverseLatticePoints1 sp action (x0, y0) =
+traverseLatticePoints1 !sp action (!x0, !y0) =
   go kMax xMax y0
   where
     forwardY  (k, y) = (k +     y + 15, y + 30)
@@ -141,20 +143,20 @@ traverseLatticePoints1 sp action (x0, y0) =
       $ (k0, x0)
 
     -- Step 4
-    adjustY
-      = head
-      . dropWhile (\(k, _) -> k < spLowBound sp)
-      . iterate forwardY
+    adjustY (!k, !y)
+      | k >= spLowBound sp
+      = (k, y)
+      | otherwise
+      = adjustY $ forwardY (k, y)
 
     -- Step 6
-    doActions (k, y)
-      = traverse_ action
-      $ takeWhile (< spHighBound sp)
-      $ map fst
-      $ iterate forwardY
-      $ (k, y)
-
-    go k x y
+    doActions (!k, !y)
+      | k < spHighBound sp
+      = action k >> doActions (forwardY (k, y))
+      | otherwise
+      = pure ()
+
+    go !k !x !y
       | x <= 0 = pure ()
       | otherwise = do
         let (k', y') = adjustY (k, y)
@@ -186,20 +188,20 @@ traverseLatticePoints2 sp action (x0, y0) =
       $ (k0, x0)
 
     -- Step 4
-    adjustY
-      = head
-      . dropWhile (\(k, _) -> k < spLowBound sp)
-      . iterate forwardY
+    adjustY (!k, !y)
+      | k >= spLowBound sp
+      = (k, y)
+      | otherwise
+      = adjustY $ forwardY (k, y)
 
     -- Step 6
-    doActions (k, y)
-      = traverse_ action
-      $ takeWhile (< spHighBound sp)
-      $ map fst
-      $ iterate forwardY
-      $ (k, y)
-
-    go k x y
+    doActions (!k, !y)
+      | k < spHighBound sp
+      = action k >> doActions (forwardY (k, y))
+      | otherwise
+      = pure ()
+
+    go !k !x !y
       | x <= 0 = pure ()
       | otherwise = do
         let (k', y') = adjustY (k, y)
@@ -222,14 +224,13 @@ traverseLatticePoints3 sp action (x0, y0) =
     k0 = (3 * x0 * x0 - y0 * y0 - spDelta sp) `quot` 60
 
     -- Step 6
-    doActions (k, x, y)
-      = traverse_ action
-      $ map fst
-      $ takeWhile (\(k', y') -> k' >= spLowBound sp && y' < x)
-      $ iterate forwardY
-      $ (k, y)
-
-    go k x y
+    doActions (!k, !x, !y)
+      | k >= spLowBound sp && y < x
+      = action k >> (let (k', y') = forwardY (k, y) in doActions (k', x, y'))
+      | otherwise
+      = pure ()
+
+    go !k !x !y
       | k >= spHighBound sp
       , x <= y
       = pure ()

From 7639b5e49c8fc2e302052f8a0f3efdb1388db91d Mon Sep 17 00:00:00 2001
From: Bodigrim <andrew.lelechenko@gmail.com>
Date: Fri, 21 Jun 2019 12:21:11 +0200
Subject: [PATCH 04/18] Shift k to start from 0, not from spLowBound

---
 Math/NumberTheory/Primes/Sieve/Atkin.hs | 54 ++++++++++++-------------
 1 file changed, 27 insertions(+), 27 deletions(-)

diff --git a/Math/NumberTheory/Primes/Sieve/Atkin.hs b/Math/NumberTheory/Primes/Sieve/Atkin.hs
index 7cad670f5..10079ee1e 100644
--- a/Math/NumberTheory/Primes/Sieve/Atkin.hs
+++ b/Math/NumberTheory/Primes/Sieve/Atkin.hs
@@ -83,7 +83,7 @@ sieveSegment
 sieveSegment sp = runST $ do
   vec <- MU.new (spLength sp)
   U.forM_ (fgs V.! toWheel30 (spDelta sp)) $
-    traverseLatticePoints sp (\k -> unsafeFlipBit vec (k - spLowBound sp))
+    traverseLatticePoints sp vec
   algo3steps456 sp vec
   U.unsafeFreeze vec
 
@@ -106,25 +106,25 @@ fgs = V.generate 16 (dispatch . fromWheel30)
 
 traverseLatticePoints
   :: SieveParams
-  -> (Int -> ST s ())
+  -> MU.MVector s Bit
   -> (Int, Int)
   -> ST s ()
-traverseLatticePoints sp action (x0, y0)
+traverseLatticePoints sp vec (x0, y0)
   | spDelta sp `mod` 4 == 1
-  = traverseLatticePoints1 sp action (x0, y0)
+  = traverseLatticePoints1 sp vec (x0, y0)
   | spDelta sp `mod` 6 == 1
-  = traverseLatticePoints2 sp action (x0, y0)
+  = traverseLatticePoints2 sp vec (x0, y0)
   | spDelta sp `mod` 12 == 11
-  = traverseLatticePoints3 sp action (x0, y0)
+  = traverseLatticePoints3 sp vec (x0, y0)
   | otherwise
   = error "traverseLatticePoints: unexpected delta"
 
 traverseLatticePoints1
   :: SieveParams
-  -> (Int -> ST s ())
+  -> MU.MVector s Bit
   -> (Int, Int)
   -> ST s ()
-traverseLatticePoints1 !sp action (!x0, !y0) =
+traverseLatticePoints1 !sp vec (!x0, !y0) =
   go kMax xMax y0
   where
     forwardY  (k, y) = (k +     y + 15, y + 30)
@@ -132,27 +132,27 @@ traverseLatticePoints1 !sp action (!x0, !y0) =
     backwardX (k, x) = (k - 2 * x + 15, x - 15)
 
     -- Step 1
-    k0 = (4 * x0 * x0 + y0 * y0 - spDelta sp) `quot` 60
+    k0 = (4 * x0 * x0 + y0 * y0 - spDelta sp) `quot` 60 - spLowBound sp
 
     -- Step 2
     (kMax, xMax)
       = backwardX
       $ head
-      $ dropWhile (\(k, _) -> k < spHighBound sp)
+      $ dropWhile (\(k, _) -> k < spLength sp)
       $ iterate forwardX
       $ (k0, x0)
 
     -- Step 4
     adjustY (!k, !y)
-      | k >= spLowBound sp
+      | k >= 0
       = (k, y)
       | otherwise
       = adjustY $ forwardY (k, y)
 
     -- Step 6
     doActions (!k, !y)
-      | k < spHighBound sp
-      = action k >> doActions (forwardY (k, y))
+      | k < spLength sp
+      = unsafeFlipBit vec k >> doActions (forwardY (k, y))
       | otherwise
       = pure ()
 
@@ -166,10 +166,10 @@ traverseLatticePoints1 !sp action (!x0, !y0) =
 
 traverseLatticePoints2
   :: SieveParams
-  -> (Int -> ST s ())
+  -> MU.MVector s Bit
   -> (Int, Int)
   -> ST s ()
-traverseLatticePoints2 sp action (x0, y0) =
+traverseLatticePoints2 sp vec (x0, y0) =
   go kMax xMax y0
   where
     forwardY  (k, y) = (k + y + 15, y + 30)
@@ -177,27 +177,27 @@ traverseLatticePoints2 sp action (x0, y0) =
     backwardX (k, x) = (k - x +  5, x - 10)
 
     -- Step 1
-    k0 = (3 * x0 * x0 + y0 * y0 - spDelta sp) `quot` 60
+    k0 = (3 * x0 * x0 + y0 * y0 - spDelta sp) `quot` 60 - spLowBound sp
 
     -- Step 2
     (kMax, xMax)
       = backwardX
       $ head
-      $ dropWhile (\(k, _) -> k < spHighBound sp)
+      $ dropWhile (\(k, _) -> k < spLength sp)
       $ iterate forwardX
       $ (k0, x0)
 
     -- Step 4
     adjustY (!k, !y)
-      | k >= spLowBound sp
+      | k >= 0
       = (k, y)
       | otherwise
       = adjustY $ forwardY (k, y)
 
     -- Step 6
     doActions (!k, !y)
-      | k < spHighBound sp
-      = action k >> doActions (forwardY (k, y))
+      | k < spLength sp
+      = unsafeFlipBit vec k >> doActions (forwardY (k, y))
       | otherwise
       = pure ()
 
@@ -211,30 +211,30 @@ traverseLatticePoints2 sp action (x0, y0) =
 
 traverseLatticePoints3
   :: SieveParams
-  -> (Int -> ST s ())
+  -> MU.MVector s Bit
   -> (Int, Int)
   -> ST s ()
-traverseLatticePoints3 sp action (x0, y0) =
+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 - spDelta sp) `quot` 60
+    k0 = (3 * x0 * x0 - y0 * y0 - spDelta sp) `quot` 60 - spLowBound sp
 
     -- Step 6
     doActions (!k, !x, !y)
-      | k >= spLowBound sp && y < x
-      = action k >> (let (k', y') = forwardY (k, y) in doActions (k', x, y'))
+      | k >= 0 && y < x
+      = unsafeFlipBit vec k >> (let (k', y') = forwardY (k, y) in doActions (k', x, y'))
       | otherwise
       = pure ()
 
     go !k !x !y
-      | k >= spHighBound sp
+      | k >= spLength sp
       , x <= y
       = pure ()
-      | k >= spHighBound sp
+      | k >= spLength sp
       = let (k', y') = forwardY (k, y) in
         go k' x y'
       | otherwise

From 9b98e57c6f8c6f45af41eb8430af0edde653bdf0 Mon Sep 17 00:00:00 2001
From: Bodigrim <andrew.lelechenko@gmail.com>
Date: Fri, 21 Jun 2019 14:07:24 +0200
Subject: [PATCH 05/18] Avoid chineseCoprime in algo3steps456

---
 Math/NumberTheory/Primes/Sieve/Atkin.hs | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/Math/NumberTheory/Primes/Sieve/Atkin.hs b/Math/NumberTheory/Primes/Sieve/Atkin.hs
index 10079ee1e..3ec4c5e8a 100644
--- a/Math/NumberTheory/Primes/Sieve/Atkin.hs
+++ b/Math/NumberTheory/Primes/Sieve/Atkin.hs
@@ -23,7 +23,6 @@ import qualified Data.Vector as V
 import qualified Data.Vector.Unboxed as U
 import qualified Data.Vector.Unboxed.Mutable as MU
 
-import Math.NumberTheory.Moduli.Chinese
 import Math.NumberTheory.Roots
 import qualified Math.NumberTheory.Primes.Sieve.Eratosthenes as E
 import Math.NumberTheory.Primes.Types
@@ -268,7 +267,15 @@ crossMultiples sp vec m =
     \k -> MU.unsafeWrite vec (k - spLowBound sp) (Bit False)
   where
     -- k0 is the smallest non-negative k such that 60k+delta = 0 (mod m)
-    k0 = (`quot` 60) $ (\k -> if k < 0 then k + 60 * m else k) $ fst $ fromJust $ chinese (spDelta sp, 60) (0, m)
+    k0 = solveCongruence (spDelta sp) m
     -- k1 = k0 (mod m), k1 >= lowBound
     (q, r) = spLowBound sp `quotRem` m
     k1 = if r < k0 then q * m + k0 else (q + 1) * m + k0
+
+-- 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)

From fa1fe2ac68267feca6dce45f8927b44766d63780 Mon Sep 17 00:00:00 2001
From: Bodigrim <andrew.lelechenko@gmail.com>
Date: Fri, 21 Jun 2019 14:34:42 +0200
Subject: [PATCH 06/18] Use listBits to explode PrimeSieve into a list of
 primes

---
 Math/NumberTheory/Primes/Sieve/Atkin.hs | 55 ++++++++++++++++++-------
 1 file changed, 41 insertions(+), 14 deletions(-)

diff --git a/Math/NumberTheory/Primes/Sieve/Atkin.hs b/Math/NumberTheory/Primes/Sieve/Atkin.hs
index 3ec4c5e8a..0711a6b49 100644
--- a/Math/NumberTheory/Primes/Sieve/Atkin.hs
+++ b/Math/NumberTheory/Primes/Sieve/Atkin.hs
@@ -17,7 +17,6 @@ module Math.NumberTheory.Primes.Sieve.Atkin
 import Control.Monad
 import Control.Monad.ST
 import Data.Bit
-import Data.Foldable
 import Data.Maybe
 import qualified Data.Vector as V
 import qualified Data.Vector.Unboxed as U
@@ -31,24 +30,52 @@ import Math.NumberTheory.Utils
 atkinPrimeList :: PrimeSieve -> [Int]
 atkinPrimeList (PrimeSieve low len segments)
   | len60 == 0 = []
-  | len60 == 1 = takeWhile (< high) $ dropWhile (< low) $ 2 : 3 : 5 : doIx 0
-  | otherwise
-  =  dropWhile (< low) (2 : 3 : 5 : doIx 0)
-  ++ concatMap doIx [1 .. len60 - 2]
-  ++ takeWhile (< high) (doIx $ len60 - 1)
+  | otherwise = takeWhile (< high) $ dropWhile (< low) $ 2 : 3 : 5 : merge l0 l1
   where
+    list00 = map (\k -> 60 * (low60 + k) + fromWheel30  0) (listBits $ segments V.!  0)
+    list01 = map (\k -> 60 * (low60 + k) + fromWheel30  1) (listBits $ segments V.!  1)
+    list02 = map (\k -> 60 * (low60 + k) + fromWheel30  2) (listBits $ segments V.!  2)
+    list03 = map (\k -> 60 * (low60 + k) + fromWheel30  3) (listBits $ segments V.!  3)
+    list04 = map (\k -> 60 * (low60 + k) + fromWheel30  4) (listBits $ segments V.!  4)
+    list05 = map (\k -> 60 * (low60 + k) + fromWheel30  5) (listBits $ segments V.!  5)
+    list06 = map (\k -> 60 * (low60 + k) + fromWheel30  6) (listBits $ segments V.!  6)
+    list07 = map (\k -> 60 * (low60 + k) + fromWheel30  7) (listBits $ segments V.!  7)
+    list08 = map (\k -> 60 * (low60 + k) + fromWheel30  8) (listBits $ segments V.!  8)
+    list09 = map (\k -> 60 * (low60 + k) + fromWheel30  9) (listBits $ segments V.!  9)
+    list10 = map (\k -> 60 * (low60 + k) + fromWheel30 10) (listBits $ segments V.! 10)
+    list11 = map (\k -> 60 * (low60 + k) + fromWheel30 11) (listBits $ segments V.! 11)
+    list12 = map (\k -> 60 * (low60 + k) + fromWheel30 12) (listBits $ segments V.! 12)
+    list13 = map (\k -> 60 * (low60 + k) + fromWheel30 13) (listBits $ segments V.! 13)
+    list14 = map (\k -> 60 * (low60 + k) + fromWheel30 14) (listBits $ segments V.! 14)
+    list15 = map (\k -> 60 * (low60 + k) + fromWheel30 15) (listBits $ segments V.! 15)
+
+    lst0 = merge list00 list01
+    lst1 = merge list02 list03
+    lst2 = merge list04 list05
+    lst3 = merge list06 list07
+    lst4 = merge list08 list09
+    lst5 = merge list10 list11
+    lst6 = merge list12 list13
+    lst7 = merge list14 list15
+
+    ls0 = merge lst0 lst1
+    ls1 = merge lst2 lst3
+    ls2 = merge lst4 lst5
+    ls3 = merge lst6 lst7
+
+    l0 = merge ls0 ls1
+    l1 = merge ls2 ls3
+
     low60 = low `quot` 60
     len60 = (low + len + 59) `quot` 60 - low60
     high = low + len
 
-    js = map fromWheel30 [0..15]
-
-    doIx k
-      = map ((+ 60 * (low60 + k)) . fst)
-      $ filter (unBit . snd)
-      $ zip js
-      $ toList
-      $ V.map (U.! k) segments
+merge :: Ord a => [a] -> [a] -> [a]
+merge [] ys = ys
+merge xs [] = xs
+merge xs@(x:xs') ys@(y:ys')
+  | x < y     = x : merge xs' ys
+  | otherwise = y : merge xs ys'
 
 data PrimeSieve = PrimeSieve
   { _psLowBound :: !Int

From 0edc954f309043fe52cc51e77fb5baeed5deb31b Mon Sep 17 00:00:00 2001
From: Bodigrim <andrew.lelechenko@gmail.com>
Date: Thu, 20 Jun 2019 22:09:57 +0200
Subject: [PATCH 07/18] Add application for profiling

---
 app/Atkin.hs   |  9 +++++++++
 arithmoi.cabal | 11 +++++++++++
 2 files changed, 20 insertions(+)
 create mode 100644 app/Atkin.hs

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 3abbc8dec..0abafcdad 100644
--- a/arithmoi.cabal
+++ b/arithmoi.cabal
@@ -238,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

From 77f2814bf4403742ec1eeec70365280fe62013c1 Mon Sep 17 00:00:00 2001
From: Bodigrim <andrew.lelechenko@gmail.com>
Date: Sat, 7 Dec 2019 18:56:29 +0000
Subject: [PATCH 08/18] Unify subsieves

---
 Math/NumberTheory/Primes/Sieve/Atkin.hs | 104 ++++++++----------------
 1 file changed, 34 insertions(+), 70 deletions(-)

diff --git a/Math/NumberTheory/Primes/Sieve/Atkin.hs b/Math/NumberTheory/Primes/Sieve/Atkin.hs
index 0711a6b49..38b2edec0 100644
--- a/Math/NumberTheory/Primes/Sieve/Atkin.hs
+++ b/Math/NumberTheory/Primes/Sieve/Atkin.hs
@@ -17,6 +17,7 @@ module Math.NumberTheory.Primes.Sieve.Atkin
 import Control.Monad
 import Control.Monad.ST
 import Data.Bit
+import Data.Bits
 import Data.Maybe
 import qualified Data.Vector as V
 import qualified Data.Vector.Unboxed as U
@@ -30,57 +31,16 @@ import Math.NumberTheory.Utils
 atkinPrimeList :: PrimeSieve -> [Int]
 atkinPrimeList (PrimeSieve low len segments)
   | len60 == 0 = []
-  | otherwise = takeWhile (< high) $ dropWhile (< low) $ 2 : 3 : 5 : merge l0 l1
+  | otherwise = takeWhile (< high) $ dropWhile (< low) $ 2 : 3 : 5 : map fromWheel30 (listBits segments)
   where
-    list00 = map (\k -> 60 * (low60 + k) + fromWheel30  0) (listBits $ segments V.!  0)
-    list01 = map (\k -> 60 * (low60 + k) + fromWheel30  1) (listBits $ segments V.!  1)
-    list02 = map (\k -> 60 * (low60 + k) + fromWheel30  2) (listBits $ segments V.!  2)
-    list03 = map (\k -> 60 * (low60 + k) + fromWheel30  3) (listBits $ segments V.!  3)
-    list04 = map (\k -> 60 * (low60 + k) + fromWheel30  4) (listBits $ segments V.!  4)
-    list05 = map (\k -> 60 * (low60 + k) + fromWheel30  5) (listBits $ segments V.!  5)
-    list06 = map (\k -> 60 * (low60 + k) + fromWheel30  6) (listBits $ segments V.!  6)
-    list07 = map (\k -> 60 * (low60 + k) + fromWheel30  7) (listBits $ segments V.!  7)
-    list08 = map (\k -> 60 * (low60 + k) + fromWheel30  8) (listBits $ segments V.!  8)
-    list09 = map (\k -> 60 * (low60 + k) + fromWheel30  9) (listBits $ segments V.!  9)
-    list10 = map (\k -> 60 * (low60 + k) + fromWheel30 10) (listBits $ segments V.! 10)
-    list11 = map (\k -> 60 * (low60 + k) + fromWheel30 11) (listBits $ segments V.! 11)
-    list12 = map (\k -> 60 * (low60 + k) + fromWheel30 12) (listBits $ segments V.! 12)
-    list13 = map (\k -> 60 * (low60 + k) + fromWheel30 13) (listBits $ segments V.! 13)
-    list14 = map (\k -> 60 * (low60 + k) + fromWheel30 14) (listBits $ segments V.! 14)
-    list15 = map (\k -> 60 * (low60 + k) + fromWheel30 15) (listBits $ segments V.! 15)
-
-    lst0 = merge list00 list01
-    lst1 = merge list02 list03
-    lst2 = merge list04 list05
-    lst3 = merge list06 list07
-    lst4 = merge list08 list09
-    lst5 = merge list10 list11
-    lst6 = merge list12 list13
-    lst7 = merge list14 list15
-
-    ls0 = merge lst0 lst1
-    ls1 = merge lst2 lst3
-    ls2 = merge lst4 lst5
-    ls3 = merge lst6 lst7
-
-    l0 = merge ls0 ls1
-    l1 = merge ls2 ls3
-
     low60 = low `quot` 60
     len60 = (low + len + 59) `quot` 60 - low60
     high = low + len
 
-merge :: Ord a => [a] -> [a] -> [a]
-merge [] ys = ys
-merge xs [] = xs
-merge xs@(x:xs') ys@(y:ys')
-  | x < y     = x : merge xs' ys
-  | otherwise = y : merge xs ys'
-
 data PrimeSieve = PrimeSieve
   { _psLowBound :: !Int
   , _psLength   :: !Int
-  , _psSegments :: V.Vector (U.Vector Bit)
+  , _psSegments :: !(U.Vector Bit)
   } deriving (Show)
 
 atkinSieve
@@ -91,8 +51,7 @@ atkinSieve low len = PrimeSieve low len segments
   where
     low60 = low `quot` 60
     len60 = (low + len + 59) `quot` 60 - low60
-    params = V.generate 16 (\i -> SieveParams (fromWheel30 i) low60 len60)
-    segments = V.map sieveSegment params
+    segments = sieveSegment low60 len60
 
 data SieveParams = SieveParams
   { spDelta    :: !Int
@@ -100,17 +59,16 @@ data SieveParams = SieveParams
   , spLength   :: !Int
   } deriving (Show)
 
-spHighBound :: SieveParams -> Int
-spHighBound sp = spLowBound sp + spLength sp
-
 sieveSegment
-  :: SieveParams
+  :: Int
+  -> Int
   -> U.Vector Bit
-sieveSegment sp = runST $ do
-  vec <- MU.new (spLength sp)
-  U.forM_ (fgs V.! toWheel30 (spDelta sp)) $
-    traverseLatticePoints sp vec
-  algo3steps456 sp vec
+sieveSegment low60 len60 = runST $ do
+  vec <- MU.new (len60 `shiftL` 4)
+  forM_ [0..15] $ \i ->
+    U.forM_ (fgs V.! i) $
+      traverseLatticePoints (SieveParams (fromWheel30 i) low60 len60) vec
+  algo3steps456 low60 len60 vec
   U.unsafeFreeze vec
 
 -- | Solutions of k * f^2 + l * g^2 = delta (mod 60)
@@ -178,7 +136,8 @@ traverseLatticePoints1 !sp vec (!x0, !y0) =
     -- Step 6
     doActions (!k, !y)
       | k < spLength sp
-      = unsafeFlipBit vec k >> doActions (forwardY (k, y))
+      = unsafeFlipBit vec (k `shiftL` 4 + toWheel30 (spDelta sp))
+        >> doActions (forwardY (k, y))
       | otherwise
       = pure ()
 
@@ -223,7 +182,8 @@ traverseLatticePoints2 sp vec (x0, y0) =
     -- Step 6
     doActions (!k, !y)
       | k < spLength sp
-      = unsafeFlipBit vec k >> doActions (forwardY (k, y))
+      = unsafeFlipBit vec (k `shiftL` 4 + toWheel30 (spDelta sp))
+        >> doActions (forwardY (k, y))
       | otherwise
       = pure ()
 
@@ -252,7 +212,8 @@ traverseLatticePoints3 sp vec (x0, y0) =
     -- Step 6
     doActions (!k, !x, !y)
       | k >= 0 && y < x
-      = unsafeFlipBit vec k >> (let (k', y') = forwardY (k, y) in doActions (k', x, y'))
+      = unsafeFlipBit vec (k `shiftL` 4 + toWheel30 (spDelta sp))
+        >> (let (k', y') = forwardY (k, y) in doActions (k', x, y'))
       | otherwise
       = pure ()
 
@@ -271,33 +232,36 @@ traverseLatticePoints3 sp vec (x0, y0) =
 
 -- | Perform steps 4-6 of Algorithm 3.X.
 algo3steps456
-  :: SieveParams
+  :: Int
+  -> Int
   -> MU.MVector s Bit
   -> ST s ()
-algo3steps456 sp vec =
+algo3steps456 low60 len60 vec =
   forM_ ps $ \p ->
-    crossMultiples sp vec (p * p)
+    crossMultiples low60 len60 vec (p * p)
   where
     low  = 7
-    high = integerSquareRoot (60 * spHighBound sp - 1)
+    high = integerSquareRoot (60 * (low60 + len60) - 1)
     ps   = takeWhile (<= high) $ dropWhile (< low) $ map unPrime E.primes
 
 -- | Cross out multiples of the first argument
 -- in a given sieve.
 crossMultiples
-  :: SieveParams
+  :: Int
+  -> Int
   -> MU.MVector s Bit
   -> Int -- coprime with 60
   -> ST s ()
-crossMultiples sp vec m =
-  forM_ [k1, k1 + m .. spHighBound sp - 1] $
-    \k -> MU.unsafeWrite vec (k - spLowBound sp) (Bit False)
+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
-    -- k0 is the smallest non-negative k such that 60k+delta = 0 (mod m)
-    k0 = solveCongruence (spDelta sp) m
-    -- k1 = k0 (mod m), k1 >= lowBound
-    (q, r) = spLowBound sp `quotRem` m
-    k1 = if r < k0 then q * m + k0 else (q + 1) * m + k0
+    (q, r) = low60 `quotRem` m
 
 -- Find the smallest k such that 60k+delta = 0 (mod m)
 -- Should be equal to

From f18b8134a3121f37a633129f7686ef5c01214a8a Mon Sep 17 00:00:00 2001
From: Bodigrim <andrew.lelechenko@gmail.com>
Date: Sat, 7 Dec 2019 20:42:43 +0000
Subject: [PATCH 09/18] Upgrade toWheel30, backporting from Chimera

---
 Math/NumberTheory/Primes/Sieve/Atkin.hs |  6 ++---
 Math/NumberTheory/Utils.hs              | 32 ++++++++++++++++++++++++-
 2 files changed, 34 insertions(+), 4 deletions(-)

diff --git a/Math/NumberTheory/Primes/Sieve/Atkin.hs b/Math/NumberTheory/Primes/Sieve/Atkin.hs
index 38b2edec0..ff24ea72f 100644
--- a/Math/NumberTheory/Primes/Sieve/Atkin.hs
+++ b/Math/NumberTheory/Primes/Sieve/Atkin.hs
@@ -136,7 +136,7 @@ traverseLatticePoints1 !sp vec (!x0, !y0) =
     -- Step 6
     doActions (!k, !y)
       | k < spLength sp
-      = unsafeFlipBit vec (k `shiftL` 4 + toWheel30 (spDelta sp))
+      = unsafeFlipBit vec (k `shiftL` 4 + toWheel30Int (spDelta sp))
         >> doActions (forwardY (k, y))
       | otherwise
       = pure ()
@@ -182,7 +182,7 @@ traverseLatticePoints2 sp vec (x0, y0) =
     -- Step 6
     doActions (!k, !y)
       | k < spLength sp
-      = unsafeFlipBit vec (k `shiftL` 4 + toWheel30 (spDelta sp))
+      = unsafeFlipBit vec (k `shiftL` 4 + toWheel30Int (spDelta sp))
         >> doActions (forwardY (k, y))
       | otherwise
       = pure ()
@@ -212,7 +212,7 @@ traverseLatticePoints3 sp vec (x0, y0) =
     -- Step 6
     doActions (!k, !x, !y)
       | k >= 0 && y < x
-      = unsafeFlipBit vec (k `shiftL` 4 + toWheel30 (spDelta sp))
+      = unsafeFlipBit vec (k `shiftL` 4 + toWheel30Int (spDelta sp))
         >> (let (k', y') = forwardY (k, y) in doActions (k', x, y'))
       | otherwise
       = pure ()
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.

From 8825b8c4c2b4506867ee51bfe37af6ae93a48a44 Mon Sep 17 00:00:00 2001
From: Bodigrim <andrew.lelechenko@gmail.com>
Date: Tue, 7 Jan 2020 00:11:01 +0000
Subject: [PATCH 10/18] Make Atkin sieve bootstrappable

---
 Math/NumberTheory/Primes/Sieve/Atkin.hs | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/Math/NumberTheory/Primes/Sieve/Atkin.hs b/Math/NumberTheory/Primes/Sieve/Atkin.hs
index ff24ea72f..12b754667 100644
--- a/Math/NumberTheory/Primes/Sieve/Atkin.hs
+++ b/Math/NumberTheory/Primes/Sieve/Atkin.hs
@@ -24,8 +24,7 @@ import qualified Data.Vector.Unboxed as U
 import qualified Data.Vector.Unboxed.Mutable as MU
 
 import Math.NumberTheory.Roots
-import qualified Math.NumberTheory.Primes.Sieve.Eratosthenes as E
-import Math.NumberTheory.Primes.Types
+import Math.NumberTheory.Primes.Small
 import Math.NumberTheory.Utils
 
 atkinPrimeList :: PrimeSieve -> [Int]
@@ -242,7 +241,9 @@ algo3steps456 low60 len60 vec =
   where
     low  = 7
     high = integerSquareRoot (60 * (low60 + len60) - 1)
-    ps   = takeWhile (<= high) $ dropWhile (< low) $ map unPrime E.primes
+    ps = case toIntegralSized high of
+      Just high' -> map fromIntegral $ smallPrimesFromTo (fromIntegral low) high'
+      Nothing    -> atkinPrimeList $ atkinSieve low (high - low + 1)
 
 -- | Cross out multiples of the first argument
 -- in a given sieve.

From ff517ac6568cf3392563a8ea9b669b4788d5a062 Mon Sep 17 00:00:00 2001
From: Bodigrim <andrew.lelechenko@gmail.com>
Date: Tue, 7 Jan 2020 19:44:59 +0000
Subject: [PATCH 11/18] Migrate TrialDivision to Atkin sieve

---
 .../Primes/Factorisation/TrialDivision.hs     | 16 ++++++-------
 Math/NumberTheory/Primes/Sieve/Atkin.hs       |  9 +++++++-
 Math/NumberTheory/Primes/Testing/Certified.hs | 23 +++++++++++--------
 3 files changed, 28 insertions(+), 20 deletions(-)

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
index 12b754667..3d40b4664 100644
--- a/Math/NumberTheory/Primes/Sieve/Atkin.hs
+++ b/Math/NumberTheory/Primes/Sieve/Atkin.hs
@@ -8,25 +8,32 @@
 --
 
 {-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE MagicHash    #-}
 
 module Math.NumberTheory.Primes.Sieve.Atkin
   ( atkinPrimeList
   , atkinSieve
+  , atkinFromTo
   ) where
 
 import Control.Monad
 import Control.Monad.ST
 import Data.Bit
 import Data.Bits
+import Data.Coerce
 import Data.Maybe
 import qualified Data.Vector as V
 import qualified Data.Vector.Unboxed as U
 import qualified Data.Vector.Unboxed.Mutable as MU
 
-import Math.NumberTheory.Roots
 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 = []
diff --git a/Math/NumberTheory/Primes/Testing/Certified.hs b/Math/NumberTheory/Primes/Testing/Certified.hs
index e5a40b854..fb825d11b 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

From a4a226580dd1f8477df26b0bb0ef9c1e95bf433b Mon Sep 17 00:00:00 2001
From: Bodigrim <andrew.lelechenko@gmail.com>
Date: Tue, 7 Jan 2020 20:25:06 +0000
Subject: [PATCH 12/18] Migrate Montgomery to Atkin sieve

---
 .../Primes/Factorisation/Montgomery.hs        | 67 ++++++++-----------
 Math/NumberTheory/Primes/Sieve/Atkin.hs       |  1 +
 Math/NumberTheory/Primes/Testing/Certified.hs | 14 ++--
 3 files changed, 37 insertions(+), 45 deletions(-)

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/Sieve/Atkin.hs b/Math/NumberTheory/Primes/Sieve/Atkin.hs
index 3d40b4664..f915369bf 100644
--- a/Math/NumberTheory/Primes/Sieve/Atkin.hs
+++ b/Math/NumberTheory/Primes/Sieve/Atkin.hs
@@ -14,6 +14,7 @@ module Math.NumberTheory.Primes.Sieve.Atkin
   ( atkinPrimeList
   , atkinSieve
   , atkinFromTo
+  , PrimeSieve
   ) where
 
 import Control.Monad
diff --git a/Math/NumberTheory/Primes/Testing/Certified.hs b/Math/NumberTheory/Primes/Testing/Certified.hs
index fb825d11b..043af5a16 100644
--- a/Math/NumberTheory/Primes/Testing/Certified.hs
+++ b/Math/NumberTheory/Primes/Testing/Certified.hs
@@ -177,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))

From 1751b19f5998d31f869026e30d580f0a68853e70 Mon Sep 17 00:00:00 2001
From: Bodigrim <andrew.lelechenko@gmail.com>
Date: Sat, 4 Jul 2020 15:43:20 +0100
Subject: [PATCH 13/18] Fix hlint suggestions

---
 Math/NumberTheory/Primes/Sieve/Atkin.hs      | 5 ++---
 benchmark/Math/NumberTheory/SequenceBench.hs | 3 +--
 2 files changed, 3 insertions(+), 5 deletions(-)

diff --git a/Math/NumberTheory/Primes/Sieve/Atkin.hs b/Math/NumberTheory/Primes/Sieve/Atkin.hs
index f915369bf..7ea9f7336 100644
--- a/Math/NumberTheory/Primes/Sieve/Atkin.hs
+++ b/Math/NumberTheory/Primes/Sieve/Atkin.hs
@@ -8,7 +8,6 @@
 --
 
 {-# LANGUAGE BangPatterns #-}
-{-# LANGUAGE MagicHash    #-}
 
 module Math.NumberTheory.Primes.Sieve.Atkin
   ( atkinPrimeList
@@ -131,7 +130,7 @@ traverseLatticePoints1 !sp vec (!x0, !y0) =
       $ head
       $ dropWhile (\(k, _) -> k < spLength sp)
       $ iterate forwardX
-      $ (k0, x0)
+        (k0, x0)
 
     -- Step 4
     adjustY (!k, !y)
@@ -177,7 +176,7 @@ traverseLatticePoints2 sp vec (x0, y0) =
       $ head
       $ dropWhile (\(k, _) -> k < spLength sp)
       $ iterate forwardX
-      $ (k0, x0)
+        (k0, x0)
 
     -- Step 4
     adjustY (!k, !y)
diff --git a/benchmark/Math/NumberTheory/SequenceBench.hs b/benchmark/Math/NumberTheory/SequenceBench.hs
index c843ef46c..2a5edde6b 100644
--- a/benchmark/Math/NumberTheory/SequenceBench.hs
+++ b/benchmark/Math/NumberTheory/SequenceBench.hs
@@ -38,8 +38,7 @@ eratosthenesBench = bgroup "eratosthenes" $
 
 atkinBench :: Benchmark
 atkinBench = bgroup "atkin" $
-  map (\(x, y) -> bench (show (x, y)) $ nf atkin (x, x + y))
-  [ (10 ^ x, 10 ^ y)
+  [ bench (show (10^x, 10^y)) $ nf atkin (10^x, 10^x + 10^y)
   | x <- [10..17]
   , y <- [6..x-1]
   , x == 10 || y == 7

From 03c49a98b6d33da05a5833567f2b98653ce225ff Mon Sep 17 00:00:00 2001
From: Bodigrim <andrew.lelechenko@gmail.com>
Date: Sun, 19 Jul 2020 19:09:07 +0100
Subject: [PATCH 14/18] Improve benchmarks of sieves

---
 benchmark/Math/NumberTheory/SequenceBench.hs | 35 ++++++++------------
 1 file changed, 14 insertions(+), 21 deletions(-)

diff --git a/benchmark/Math/NumberTheory/SequenceBench.hs b/benchmark/Math/NumberTheory/SequenceBench.hs
index 2a5edde6b..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
@@ -13,42 +15,33 @@ import Data.Maybe
 import Math.NumberTheory.Primes
 
 filterIsPrime :: (Integer, Integer) -> Integer
-filterIsPrime (p, q) = sum $ takeWhile (<= q) $ dropWhile (< p) $ filter (isJust . 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 - p)
+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
-  ]
-
-atkinBench :: Benchmark
-atkinBench = bgroup "atkin" $
-  [ bench (show (10^x, 10^y)) $ nf atkin (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
-    , atkinBench
+    [ sieveBench
+    , filterIsPrimeBench
     ]
 
 -------------------------------------------------------------------------------

From 265d1d077f2b2850f51401787c4968c4a4e2ae77 Mon Sep 17 00:00:00 2001
From: Bodigrim <andrew.lelechenko@gmail.com>
Date: Tue, 29 Dec 2020 21:56:20 +0000
Subject: [PATCH 15/18] Add cabal.project

---
 cabal.project | 7 +++++++
 1 file changed, 7 insertions(+)
 create mode 100644 cabal.project

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

From 610003c72892f0a8018c07b1933d75d6a33d9f16 Mon Sep 17 00:00:00 2001
From: Bodigrim <andrew.lelechenko@gmail.com>
Date: Tue, 29 Dec 2020 21:56:39 +0000
Subject: [PATCH 16/18] Minor tweaks

---
 Math/NumberTheory/Primes/Sieve/Atkin.hs | 38 +++++++++++++++----------
 Math/NumberTheory/Primes/Small.hs       |  1 +
 2 files changed, 24 insertions(+), 15 deletions(-)

diff --git a/Math/NumberTheory/Primes/Sieve/Atkin.hs b/Math/NumberTheory/Primes/Sieve/Atkin.hs
index 7ea9f7336..4904e1ad4 100644
--- a/Math/NumberTheory/Primes/Sieve/Atkin.hs
+++ b/Math/NumberTheory/Primes/Sieve/Atkin.hs
@@ -25,6 +25,7 @@ import Data.Maybe
 import qualified Data.Vector as V
 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
@@ -60,20 +61,23 @@ atkinSieve low len = PrimeSieve low len segments
     segments = sieveSegment low60 len60
 
 data SieveParams = SieveParams
-  { spDelta    :: !Int
+  { spDelta16  :: !Int
+  , spDelta60  :: !Int
+  -- ^ 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_ [0..15] $ \i ->
-    U.forM_ (fgs V.! i) $
-      traverseLatticePoints (SieveParams (fromWheel30 i) low60 len60) vec
+  flip V.imapM_ fgs $ \delta16 -> U.mapM_ $
+    traverseLatticePoints (SieveParams delta16 (fromWheel30 delta16) low60 len60) vec
   algo3steps456 low60 len60 vec
   U.unsafeFreeze vec
 
@@ -81,6 +85,8 @@ sieveSegment low60 len60 = runST $ do
 -- 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 :: V.Vector (U.Vector (Int, Int))
 fgs = V.generate 16 (dispatch . fromWheel30)
   where
@@ -100,11 +106,11 @@ traverseLatticePoints
   -> (Int, Int)
   -> ST s ()
 traverseLatticePoints sp vec (x0, y0)
-  | spDelta sp `mod` 4 == 1
+  | spDelta60 sp `mod` 4 == 1
   = traverseLatticePoints1 sp vec (x0, y0)
-  | spDelta sp `mod` 6 == 1
+  | spDelta60 sp `mod` 6 == 1
   = traverseLatticePoints2 sp vec (x0, y0)
-  | spDelta sp `mod` 12 == 11
+  | spDelta60 sp `mod` 12 == 11
   = traverseLatticePoints3 sp vec (x0, y0)
   | otherwise
   = error "traverseLatticePoints: unexpected delta"
@@ -122,7 +128,7 @@ traverseLatticePoints1 !sp vec (!x0, !y0) =
     backwardX (k, x) = (k - 2 * x + 15, x - 15)
 
     -- Step 1
-    k0 = (4 * x0 * x0 + y0 * y0 - spDelta sp) `quot` 60 - spLowBound sp
+    k0 = (4 * x0 * x0 + y0 * y0 - spDelta60 sp) `quot` 60 - spLowBound sp
 
     -- Step 2
     (kMax, xMax)
@@ -142,7 +148,7 @@ traverseLatticePoints1 !sp vec (!x0, !y0) =
     -- Step 6
     doActions (!k, !y)
       | k < spLength sp
-      = unsafeFlipBit vec (k `shiftL` 4 + toWheel30Int (spDelta sp))
+      = unsafeFlipBit vec (k `shiftL` 4 + spDelta16 sp)
         >> doActions (forwardY (k, y))
       | otherwise
       = pure ()
@@ -168,7 +174,7 @@ traverseLatticePoints2 sp vec (x0, y0) =
     backwardX (k, x) = (k - x +  5, x - 10)
 
     -- Step 1
-    k0 = (3 * x0 * x0 + y0 * y0 - spDelta sp) `quot` 60 - spLowBound sp
+    k0 = (3 * x0 * x0 + y0 * y0 - spDelta60 sp) `quot` 60 - spLowBound sp
 
     -- Step 2
     (kMax, xMax)
@@ -188,7 +194,7 @@ traverseLatticePoints2 sp vec (x0, y0) =
     -- Step 6
     doActions (!k, !y)
       | k < spLength sp
-      = unsafeFlipBit vec (k `shiftL` 4 + toWheel30Int (spDelta sp))
+      = unsafeFlipBit vec (k `shiftL` 4 + spDelta16 sp)
         >> doActions (forwardY (k, y))
       | otherwise
       = pure ()
@@ -213,12 +219,12 @@ traverseLatticePoints3 sp vec (x0, y0) =
     forwardX  (k, x) = (k + x +  5, x + 10)
 
     -- Step 1
-    k0 = (3 * x0 * x0 - y0 * y0 - spDelta sp) `quot` 60 - spLowBound sp
+    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 + toWheel30Int (spDelta sp))
+      = unsafeFlipBit vec (k `shiftL` 4 + spDelta16 sp)
         >> (let (k', y') = forwardY (k, y) in doActions (k', x, y'))
       | otherwise
       = pure ()
@@ -249,8 +255,10 @@ algo3steps456 low60 len60 vec =
     low  = 7
     high = integerSquareRoot (60 * (low60 + len60) - 1)
     ps = case toIntegralSized high of
-      Just high' -> map fromIntegral $ smallPrimesFromTo (fromIntegral low) high'
-      Nothing    -> atkinPrimeList $ atkinSieve low (high - low + 1)
+      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.
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

From c057bfa3f9d5a4b6d4cc2dd28838da0bbc965ba7 Mon Sep 17 00:00:00 2001
From: Bodigrim <andrew.lelechenko@gmail.com>
Date: Tue, 29 Dec 2020 22:58:00 +0000
Subject: [PATCH 17/18] Write proper test and debug it

---
 Math/NumberTheory/Primes.hs                       |  3 +--
 Math/NumberTheory/Primes/Sieve/Atkin.hs           |  5 +++--
 test-suite/Math/NumberTheory/Primes/SieveTests.hs | 13 ++++++++-----
 3 files changed, 12 insertions(+), 9 deletions(-)

diff --git a/Math/NumberTheory/Primes.hs b/Math/NumberTheory/Primes.hs
index 9c2530843..7168f2b7e 100644
--- a/Math/NumberTheory/Primes.hs
+++ b/Math/NumberTheory/Primes.hs
@@ -21,8 +21,7 @@ module Math.NumberTheory.Primes
     , -- * Old interface
       primes
     , -- * Temporary
-      atkinPrimeList
-    , atkinSieve
+      module Math.NumberTheory.Primes.Sieve.Atkin
     ) where
 
 import Data.Bits
diff --git a/Math/NumberTheory/Primes/Sieve/Atkin.hs b/Math/NumberTheory/Primes/Sieve/Atkin.hs
index 4904e1ad4..2bf2f753c 100644
--- a/Math/NumberTheory/Primes/Sieve/Atkin.hs
+++ b/Math/NumberTheory/Primes/Sieve/Atkin.hs
@@ -38,7 +38,7 @@ 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 fromWheel30 (listBits segments)
+  | 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
@@ -248,6 +248,7 @@ algo3steps456
   -> Int
   -> MU.MVector s Bit
   -> ST s ()
+algo3steps456 _ 0 _ = pure ()
 algo3steps456 low60 len60 vec =
   forM_ ps $ \p ->
     crossMultiples low60 len60 vec (p * p)
@@ -273,7 +274,7 @@ crossMultiples low60 len60 vec m =
         -- 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
+        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
diff --git a/test-suite/Math/NumberTheory/Primes/SieveTests.hs b/test-suite/Math/NumberTheory/Primes/SieveTests.hs
index 3effa5ef5..b6b6931a3 100644
--- a/test-suite/Math/NumberTheory/Primes/SieveTests.hs
+++ b/test-suite/Math/NumberTheory/Primes/SieveTests.hs
@@ -21,6 +21,7 @@ import Prelude hiding (words)
 
 import Test.Tasty
 import Test.Tasty.HUnit
+import Test.Tasty.QuickCheck
 
 import Data.Bits
 import Data.Int
@@ -49,12 +50,14 @@ primesProperty2 _ = assertEqual "primes matches isPrime"
   (map unPrime primes :: [a])
   (filter (isJust . isPrime . toInteger) [1..maxBound])
 
-atkinPrimesProperty1 :: Assertion
-atkinPrimesProperty1 = assertEqual "atkinPrimes matches isPrime" expected actual
+atkinPrimesProperty1 :: Large Word -> Large Word -> Property
+atkinPrimesProperty1 (Large x) (Large y) = actual === expected
   where
     lim = 1000000
-    expected = filter (isJust . isPrime . toInteger) [1..lim]
-    actual   = atkinPrimeList $ atkinSieve 1 lim
+    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
@@ -77,7 +80,7 @@ psieveListProperty2 _ = assertEqual "primes == primeList . psieveList"
 
 testSuite :: TestTree
 testSuite = testGroup "Sieve"
-  [ testCase "atkinPrimes" atkinPrimesProperty1
+  [ testProperty "atkinPrimes" atkinPrimesProperty1
   , testGroup "primes"
     [ testCase "Int"     (primesProperty1 (Proxy :: Proxy Int))
     , testCase "Word"    (primesProperty1 (Proxy :: Proxy Word))

From 2361ae1bca33ea0ba3379bd1832f1f15d6096e84 Mon Sep 17 00:00:00 2001
From: Bodigrim <andrew.lelechenko@gmail.com>
Date: Sun, 24 Jan 2021 23:24:52 +0000
Subject: [PATCH 18/18] Minor changes

---
 Math/NumberTheory/Primes/Sieve/Atkin.hs | 17 ++++++++---------
 1 file changed, 8 insertions(+), 9 deletions(-)

diff --git a/Math/NumberTheory/Primes/Sieve/Atkin.hs b/Math/NumberTheory/Primes/Sieve/Atkin.hs
index 2bf2f753c..88d260d7b 100644
--- a/Math/NumberTheory/Primes/Sieve/Atkin.hs
+++ b/Math/NumberTheory/Primes/Sieve/Atkin.hs
@@ -22,7 +22,6 @@ import Data.Bit
 import Data.Bits
 import Data.Coerce
 import Data.Maybe
-import qualified Data.Vector as V
 import qualified Data.Vector.Unboxed as U
 import qualified Data.Vector.Unboxed.Mutable as MU
 import Data.Word
@@ -61,8 +60,8 @@ atkinSieve low len = PrimeSieve low len segments
     segments = sieveSegment low60 len60
 
 data SieveParams = SieveParams
-  { spDelta16  :: !Int
-  , spDelta60  :: !Int
+  { spDelta16  :: !Int -- [0..15]
+  , spDelta60  :: !Int -- [0..59]
   -- ^ spDelta30 = fromWheel30 spDelta16
   , spLowBound :: !Int
   , spLength   :: !Int
@@ -76,7 +75,7 @@ sieveSegment
   -> U.Vector Bit
 sieveSegment low60 len60 = runST $ do
   vec <- MU.new (len60 `shiftL` 4)
-  flip V.imapM_ fgs $ \delta16 -> U.mapM_ $
+  forM_ fgs $ uncurry $ \delta16 -> traverse $
     traverseLatticePoints (SieveParams delta16 (fromWheel30 delta16) low60 len60) vec
   algo3steps456 low60 len60 vec
   U.unsafeFreeze vec
@@ -87,16 +86,16 @@ sieveSegment low60 len60 = runST $ do
 --              = (3,-1) for delta =11 (mod 12)
 -- Outer vector is indexed by delta=[0..15],
 -- inner vector is a list of pairs.
-fgs :: V.Vector (U.Vector (Int, Int))
-fgs = V.generate 16 (dispatch . fromWheel30)
+fgs :: [(Int, [(Int, Int)])]
+fgs = map (\i -> (i, dispatch (fromWheel30 i))) [0..15]
   where
     dispatch delta
       | delta `mod` 4 == 1
-      = U.fromList [ (f, g) | f <- [1..15], g <- [1..30], (4*f*f + g*g - delta) `rem` 60 == 0]
+      = [ (f, g) | f <- [1..15], g <- [1..30], (4*f*f + g*g - delta) `rem` 60 == 0]
       | delta `mod` 6 == 1
-      = U.fromList [ (f, g) | f <- [1..10], g <- [1..30], (3*f*f + g*g - delta) `rem` 60 == 0]
+      = [ (f, g) | f <- [1..10], g <- [1..30], (3*f*f + g*g - delta) `rem` 60 == 0]
       | delta `mod` 12 == 11
-      = U.fromList [ (f, g) | f <- [1..10], g <- [1..30], (3*f*f - g*g - delta) `rem` 60 == 0]
+      = [ (f, g) | f <- [1..10], g <- [1..30], (3*f*f - g*g - delta) `rem` 60 == 0]
       | otherwise
       = error "fgs: unexpected delta"