module Factory.Math.Implementations.Primes.SieveOfAtkin(
sieveOfAtkin,
) where
import qualified Control.DeepSeq
import qualified Control.Parallel.Strategies
import qualified Data.Array.IArray
import Data.Array.IArray((!))
import qualified Data.IntSet
import qualified Data.List
import qualified Data.Set
import qualified Factory.Data.PrimeWheel as Data.PrimeWheel
import qualified Factory.Math.Power as Math.Power
import qualified ToolShed.Data.List
data PolynomialType
= ModFour
| ModSix
| ModTwelve
| None
deriving PolynomialType -> PolynomialType -> Bool
(PolynomialType -> PolynomialType -> Bool)
-> (PolynomialType -> PolynomialType -> Bool) -> Eq PolynomialType
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: PolynomialType -> PolynomialType -> Bool
$c/= :: PolynomialType -> PolynomialType -> Bool
== :: PolynomialType -> PolynomialType -> Bool
$c== :: PolynomialType -> PolynomialType -> Bool
Eq
atkinsModulus :: Integral i => i
atkinsModulus :: i
atkinsModulus = (i -> i -> i) -> [i] -> i
forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldr1 i -> i -> i
forall a. Integral a => a -> a -> a
lcm [i
4, i
6, i
12]
inherentPrimes :: Integral i => [i]
inherentPrimes :: [i]
inherentPrimes = [i
2, i
3]
nInherentPrimes :: Int
nInherentPrimes :: Int
nInherentPrimes = [Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Int]
forall i. Integral i => [i]
inherentPrimes :: [Int])
getPrefactoredPrimes :: Integral i => Data.PrimeWheel.PrimeWheel i -> [i]
getPrefactoredPrimes :: PrimeWheel i -> [i]
getPrefactoredPrimes = [i] -> [i] -> [i]
forall a. Ord a => a -> a -> a
max [i]
forall i. Integral i => [i]
inherentPrimes ([i] -> [i]) -> (PrimeWheel i -> [i]) -> PrimeWheel i -> [i]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. PrimeWheel i -> [i]
forall i. PrimeWheel i -> [i]
Data.PrimeWheel.getPrimeComponents
polynomialTypeLookupPeriod :: Integral i => Data.PrimeWheel.PrimeWheel i -> i
polynomialTypeLookupPeriod :: PrimeWheel i -> i
polynomialTypeLookupPeriod = i -> i -> i
forall a. Integral a => a -> a -> a
lcm i
forall i. Integral i => i
atkinsModulus (i -> i) -> (PrimeWheel i -> i) -> PrimeWheel i -> i
forall b c a. (b -> c) -> (a -> b) -> a -> c
. PrimeWheel i -> i
forall i. Integral i => PrimeWheel i -> i
Data.PrimeWheel.getCircumference
polynomialTypeLookup :: (Data.Array.IArray.Ix i, Integral i)
=> Data.PrimeWheel.PrimeWheel i
-> i
-> Data.Array.IArray.Array i PolynomialType
polynomialTypeLookup :: PrimeWheel i -> i -> Array i PolynomialType
polynomialTypeLookup PrimeWheel i
primeWheel i
maxPrime = (i, i) -> [PolynomialType] -> Array i PolynomialType
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
(i, i) -> [e] -> a i e
Data.Array.IArray.listArray (i
0, i -> i
forall a. Enum a => a -> a
pred (PrimeWheel i -> i
forall i. Integral i => PrimeWheel i -> i
polynomialTypeLookupPeriod PrimeWheel i
primeWheel) i -> i -> i
forall a. Ord a => a -> a -> a
`min` i
maxPrime) ([PolynomialType] -> Array i PolynomialType)
-> [PolynomialType] -> Array i PolynomialType
forall a b. (a -> b) -> a -> b
$ (i -> PolynomialType) -> [i] -> [PolynomialType]
forall a b. (a -> b) -> [a] -> [b]
map i -> PolynomialType
select [i
0 ..] where
select :: i -> PolynomialType
select i
n
| (i -> Bool) -> [i] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (
(i -> i -> Bool
forall a. Eq a => a -> a -> Bool
== i
0) (i -> Bool) -> (i -> i) -> i -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (i
n i -> i -> i
forall a. Integral a => a -> a -> a
`rem`)
) [i]
primeComponents = PolynomialType
None
| i
r i -> [i] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [i
1, i
5] = PolynomialType
ModFour
| i
r i -> i -> Bool
forall a. Eq a => a -> a -> Bool
== i
7 = PolynomialType
ModSix
| i
r i -> i -> Bool
forall a. Eq a => a -> a -> Bool
== i
11 = PolynomialType
ModTwelve
| Bool
otherwise = PolynomialType
None
where
r :: i
r = i
n i -> i -> i
forall a. Integral a => a -> a -> a
`rem` i
forall i. Integral i => i
atkinsModulus
primeComponents :: [i]
primeComponents = Int -> [i] -> [i]
forall a. Int -> [a] -> [a]
drop Int
nInherentPrimes ([i] -> [i]) -> [i] -> [i]
forall a b. (a -> b) -> a -> b
$ PrimeWheel i -> [i]
forall i. PrimeWheel i -> [i]
Data.PrimeWheel.getPrimeComponents PrimeWheel i
primeWheel
squares :: Integral i => [i]
squares :: [i]
squares = ((i, i) -> i) -> [(i, i)] -> [i]
forall a b. (a -> b) -> [a] -> [b]
map (i, i) -> i
forall a b. (a, b) -> b
snd ([(i, i)] -> [i]) -> [(i, i)] -> [i]
forall a b. (a -> b) -> a -> b
$ i -> [(i, i)]
forall n. (Enum n, Num n) => n -> [(n, n)]
Math.Power.squaresFrom i
1
filterOddRepetitions :: Ord a => [a] -> [a]
filterOddRepetitions :: [a] -> [a]
filterOddRepetitions = Bool -> [a] -> [a]
forall a. Eq a => Bool -> [a] -> [a]
slave Bool
True ([a] -> [a]) -> ([a] -> [a]) -> [a] -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> [a]
forall a. Ord a => [a] -> [a]
Data.List.sort where
slave :: Bool -> [a] -> [a]
slave Bool
isOdd (a
one : remainder :: [a]
remainder@(a
two : [a]
_))
| a
one a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
two = Bool -> [a] -> [a]
slave (Bool -> Bool
not Bool
isOdd) [a]
remainder
| Bool
isOdd = a
one a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a]
beginSpan
| Bool
otherwise = [a]
beginSpan
where
beginSpan :: [a]
beginSpan = Bool -> [a] -> [a]
slave Bool
True [a]
remainder
slave Bool
True [a
singleton] = [a
singleton]
slave Bool
_ [a]
_ = []
findPolynomialSolutions :: (Control.DeepSeq.NFData i, Data.Array.IArray.Ix i, Integral i)
=> Data.PrimeWheel.PrimeWheel i
-> i
-> [i]
findPolynomialSolutions :: PrimeWheel i -> i -> [i]
findPolynomialSolutions PrimeWheel i
primeWheel i
maxPrime = ([i] -> [i] -> [i]) -> [[i]] -> [i]
forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldr1 [i] -> [i] -> [i]
forall a. Ord a => [a] -> [a] -> [a]
ToolShed.Data.List.merge ([[i]] -> [i]) -> [[i]] -> [i]
forall a b. (a -> b) -> a -> b
$ Strategy [[i]] -> [[i]] -> [[i]]
forall a. Strategy a -> a -> a
Control.Parallel.Strategies.withStrategy (
Strategy [i] -> Strategy [[i]]
forall a. Strategy a -> Strategy [a]
Control.Parallel.Strategies.parList Strategy [i]
forall a. NFData a => Strategy a
Control.Parallel.Strategies.rdeepseq
) [
{-# SCC "4x^2+y^2" #-} [i] -> [i]
forall a. Ord a => [a] -> [a]
filterOddRepetitions [
i
z |
i
x' <- (i -> Bool) -> [i] -> [i]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (i -> i -> Bool
forall a. Ord a => a -> a -> Bool
<= i -> i
forall a. Enum a => a -> a
pred i
maxPrime) ([i] -> [i]) -> [i] -> [i]
forall a b. (a -> b) -> a -> b
$ (i -> i) -> [i] -> [i]
forall a b. (a -> b) -> [a] -> [b]
map (i -> i -> i
forall a. Num a => a -> a -> a
* i
4) [i]
forall i. Integral i => [i]
squares,
i
z <- (i -> Bool) -> [i] -> [i]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (i -> i -> Bool
forall a. Ord a => a -> a -> Bool
<= i
maxPrime) ([i] -> [i]) -> [i] -> [i]
forall a b. (a -> b) -> a -> b
$ (i -> i) -> [i] -> [i]
forall a b. (a -> b) -> [a] -> [b]
map (i -> i -> i
forall a. Num a => a -> a -> a
+ i
x') [i]
oddSquares,
i -> PolynomialType
lookupPolynomialType i
z PolynomialType -> PolynomialType -> Bool
forall a. Eq a => a -> a -> Bool
== PolynomialType
ModFour
],
{-# SCC "3x^2+y^2" #-} [i] -> [i]
forall a. Ord a => [a] -> [a]
filterOddRepetitions [
i
z |
i
x' <- (i -> Bool) -> [i] -> [i]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (i -> i -> Bool
forall a. Ord a => a -> a -> Bool
<= i -> i
forall a. Enum a => a -> a
pred i
maxPrime) ([i] -> [i]) -> [i] -> [i]
forall a b. (a -> b) -> a -> b
$ (i -> i) -> [i] -> [i]
forall a b. (a -> b) -> [a] -> [b]
map (i -> i -> i
forall a. Num a => a -> a -> a
* i
3) [i]
forall i. Integral i => [i]
squares,
i
z <- (i -> Bool) -> [i] -> [i]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (i -> i -> Bool
forall a. Ord a => a -> a -> Bool
<= i
maxPrime) ([i] -> [i]) -> ([i] -> [i]) -> [i] -> [i]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (i -> i) -> [i] -> [i]
forall a b. (a -> b) -> [a] -> [b]
map (i -> i -> i
forall a. Num a => a -> a -> a
+ i
x') ([i] -> [i]) -> [i] -> [i]
forall a b. (a -> b) -> a -> b
$ if i -> Bool
forall a. Integral a => a -> Bool
even i
x' then [i]
oddSelection else [i]
evenSelection,
i -> PolynomialType
lookupPolynomialType i
z PolynomialType -> PolynomialType -> Bool
forall a. Eq a => a -> a -> Bool
== PolynomialType
ModSix
],
{-# SCC "3x^2-y^2" #-} [i] -> [i]
forall a. Ord a => [a] -> [a]
filterOddRepetitions [
i
z |
i
x2 <- (i -> Bool) -> [i] -> [i]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (i -> i -> Bool
forall a. Ord a => a -> a -> Bool
<= i
maxPrime i -> i -> i
forall a. Integral a => a -> a -> a
`div` i
2) [i]
forall i. Integral i => [i]
squares,
i
z <- (i -> Bool) -> [i] -> [i]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (i -> i -> Bool
forall a. Ord a => a -> a -> Bool
> i
maxPrime) ([i] -> [i]) -> ([i] -> [i]) -> [i] -> [i]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (i -> i) -> [i] -> [i]
forall a b. (a -> b) -> [a] -> [b]
map (i
3 i -> i -> i
forall a. Num a => a -> a -> a
* i
x2 i -> i -> i
forall a. Num a => a -> a -> a
-) ([i] -> [i]) -> ([i] -> [i]) -> [i] -> [i]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (i -> Bool) -> [i] -> [i]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (i -> i -> Bool
forall a. Ord a => a -> a -> Bool
< i
x2) ([i] -> [i]) -> [i] -> [i]
forall a b. (a -> b) -> a -> b
$ if i -> Bool
forall a. Integral a => a -> Bool
even i
x2 then [i]
oddSelection else [i]
evenSelection,
i -> PolynomialType
lookupPolynomialType i
z PolynomialType -> PolynomialType -> Bool
forall a. Eq a => a -> a -> Bool
== PolynomialType
ModTwelve
]
] where
([i]
evenSquares, [i]
oddSquares) = (i -> Bool) -> [i] -> ([i], [i])
forall a. (a -> Bool) -> [a] -> ([a], [a])
Data.List.partition i -> Bool
forall a. Integral a => a -> Bool
even [i]
forall i. Integral i => [i]
squares
evenSelection :: [i]
evenSelection = [i] -> [i]
forall a. [a] -> [a]
selection110 [i]
evenSquares where
selection110 :: [a] -> [a]
selection110 (a
x0 : a
x1 : a
_ : [a]
xs) = a
x0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a
x1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a]
selection110 [a]
xs
selection110 [a]
xs = [a]
xs
oddSelection :: [i]
oddSelection = [i] -> [i]
forall a. [a] -> [a]
selection101 [i]
oddSquares where
selection101 :: [a] -> [a]
selection101 (a
x0 : a
_ : a
x2 : [a]
xs) = a
x0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a
x2 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a]
selection101 [a]
xs
selection101 [a]
xs = [a]
xs
lookupPolynomialType :: i -> PolynomialType
lookupPolynomialType = (PrimeWheel i -> i -> Array i PolynomialType
forall i.
(Ix i, Integral i) =>
PrimeWheel i -> i -> Array i PolynomialType
polynomialTypeLookup PrimeWheel i
primeWheel i
maxPrime Array i PolynomialType -> i -> PolynomialType
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
!) (i -> PolynomialType) -> (i -> i) -> i -> PolynomialType
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (i -> i -> i
forall a. Integral a => a -> a -> a
`rem` PrimeWheel i -> i
forall i. Integral i => PrimeWheel i -> i
polynomialTypeLookupPeriod PrimeWheel i
primeWheel)
generateMultiplesOfSquareTo :: Integral i
=> Data.PrimeWheel.PrimeWheel i
-> i
-> i
-> [i]
generateMultiplesOfSquareTo :: PrimeWheel i -> i -> i -> [i]
generateMultiplesOfSquareTo PrimeWheel i
primeWheel i
prime i
max' = (i -> Bool) -> [i] -> [i]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (i -> i -> Bool
forall a. Ord a => a -> a -> Bool
<= i
max') ([i] -> [i]) -> ([i] -> [i]) -> [i] -> [i]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (i -> i -> i) -> i -> [i] -> [i]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (\i
accumulator -> (i -> i -> i
forall a. Num a => a -> a -> a
+ i
accumulator) (i -> i) -> (i -> i) -> i -> i
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (i -> i -> i
forall a. Num a => a -> a -> a
* i
prime2)) i
prime2 ([i] -> [i]) -> ([i] -> [i]) -> [i] -> [i]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [i] -> [i]
forall a. [a] -> [a]
cycle ([i] -> [i]) -> [i] -> [i]
forall a b. (a -> b) -> a -> b
$ PrimeWheel i -> [i]
forall i. PrimeWheel i -> [i]
Data.PrimeWheel.getSpokeGaps PrimeWheel i
primeWheel where
prime2 :: i
prime2 = i -> i
forall n. Num n => n -> n
Math.Power.square i
prime
sieveOfAtkin :: (Control.DeepSeq.NFData i, Data.Array.IArray.Ix i, Integral i)
=> Data.PrimeWheel.NPrimes
-> i
-> [i]
sieveOfAtkin :: Int -> i -> [i]
sieveOfAtkin Int
wheelSize i
maxPrime = ([i]
prefactoredPrimes [i] -> [i] -> [i]
forall a. [a] -> [a] -> [a]
++) ([i] -> [i]) -> ([i] -> [i]) -> [i] -> [i]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Set i -> [i] -> [i]
filterSquareFree Set i
forall a. Set a
Data.Set.empty ([i] -> [i]) -> ([i] -> [i]) -> [i] -> [i]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (i -> Bool) -> [i] -> [i]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (i -> i -> Bool
forall a. Ord a => a -> a -> Bool
<= [i] -> i
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [i]
prefactoredPrimes) ([i] -> [i]) -> [i] -> [i]
forall a b. (a -> b) -> a -> b
$ PrimeWheel i -> i -> [i]
forall i. (NFData i, Ix i, Integral i) => PrimeWheel i -> i -> [i]
findPolynomialSolutions PrimeWheel i
primeWheel i
maxPrime where
primeWheel :: PrimeWheel i
primeWheel = Int -> PrimeWheel i
forall i. Integral i => Int -> PrimeWheel i
Data.PrimeWheel.mkPrimeWheel Int
wheelSize
prefactoredPrimes :: [i]
prefactoredPrimes = PrimeWheel i -> [i]
forall i. Integral i => PrimeWheel i -> [i]
getPrefactoredPrimes PrimeWheel i
primeWheel
filterSquareFree :: Set i -> [i] -> [i]
filterSquareFree Set i
_ [] = []
filterSquareFree Set i
primeMultiples (i
candidate : [i]
candidates)
| i -> Set i -> Bool
forall a. Ord a => a -> Set a -> Bool
Data.Set.member i
candidate Set i
primeMultiples = {-# SCC "delete" #-} Set i -> [i] -> [i]
filterSquareFree (i -> Set i -> Set i
forall a. Ord a => a -> Set a -> Set a
Data.Set.delete i
candidate Set i
primeMultiples) [i]
candidates
| Bool
otherwise = {-# SCC "insert" #-} i
candidate i -> [i] -> [i]
forall a. a -> [a] -> [a]
: Set i -> [i] -> [i]
filterSquareFree (Set i -> Set i -> Set i
forall a. Ord a => Set a -> Set a -> Set a
Data.Set.union Set i
primeMultiples (Set i -> Set i) -> ([i] -> Set i) -> [i] -> Set i
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [i] -> Set i
forall a. [a] -> Set a
Data.Set.fromDistinctAscList ([i] -> Set i) -> [i] -> Set i
forall a b. (a -> b) -> a -> b
$ PrimeWheel i -> i -> i -> [i]
forall i. Integral i => PrimeWheel i -> i -> i -> [i]
generateMultiplesOfSquareTo PrimeWheel i
primeWheel i
candidate i
maxPrime) [i]
candidates
{-# NOINLINE sieveOfAtkin #-}
{-# RULES "sieveOfAtkin/Int" sieveOfAtkin = sieveOfAtkinInt #-}
sieveOfAtkinInt :: Data.PrimeWheel.NPrimes -> Int -> [Int]
sieveOfAtkinInt :: Int -> Int -> [Int]
sieveOfAtkinInt Int
wheelSize Int
maxPrime = ([Int]
prefactoredPrimes [Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++) ([Int] -> [Int]) -> ([Int] -> [Int]) -> [Int] -> [Int]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IntSet -> [Int] -> [Int]
filterSquareFree IntSet
Data.IntSet.empty ([Int] -> [Int]) -> ([Int] -> [Int]) -> [Int] -> [Int]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> Bool) -> [Int] -> [Int]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Int]
prefactoredPrimes) ([Int] -> [Int]) -> [Int] -> [Int]
forall a b. (a -> b) -> a -> b
$ PrimeWheel Int -> Int -> [Int]
forall i. (NFData i, Ix i, Integral i) => PrimeWheel i -> i -> [i]
findPolynomialSolutions PrimeWheel Int
primeWheel Int
maxPrime where
primeWheel :: PrimeWheel Int
primeWheel = Int -> PrimeWheel Int
forall i. Integral i => Int -> PrimeWheel i
Data.PrimeWheel.mkPrimeWheel Int
wheelSize
prefactoredPrimes :: [Int]
prefactoredPrimes = PrimeWheel Int -> [Int]
forall i. Integral i => PrimeWheel i -> [i]
getPrefactoredPrimes PrimeWheel Int
primeWheel
filterSquareFree :: Data.IntSet.IntSet -> [Int] -> [Int]
filterSquareFree :: IntSet -> [Int] -> [Int]
filterSquareFree IntSet
_ [] = []
filterSquareFree IntSet
primeMultiples (Int
candidate : [Int]
candidates)
| Int -> IntSet -> Bool
Data.IntSet.member Int
candidate IntSet
primeMultiples = IntSet -> [Int] -> [Int]
filterSquareFree (Int -> IntSet -> IntSet
Data.IntSet.delete Int
candidate IntSet
primeMultiples) [Int]
candidates
| Bool
otherwise = Int
candidate Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: IntSet -> [Int] -> [Int]
filterSquareFree (IntSet -> IntSet -> IntSet
Data.IntSet.union IntSet
primeMultiples (IntSet -> IntSet) -> ([Int] -> IntSet) -> [Int] -> IntSet
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Int] -> IntSet
Data.IntSet.fromDistinctAscList ([Int] -> IntSet) -> [Int] -> IntSet
forall a b. (a -> b) -> a -> b
$ PrimeWheel Int -> Int -> Int -> [Int]
forall i. Integral i => PrimeWheel i -> i -> i -> [i]
generateMultiplesOfSquareTo PrimeWheel Int
primeWheel Int
candidate Int
maxPrime) [Int]
candidates