{-
	Copyright (C) 2011 Dr. Alistair Ward

	This program is free software: you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation, either version 3 of the License, or
	(at your option) any later version.

	This program is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with this program.  If not, see <https://www.gnu.org/licenses/>.
-}
{- |
 [@AUTHOR@]	Dr. Alistair Ward

 [@DESCRIPTION@]

	* Generates the constant /bounded/ list of /prime-numbers/, using the /Sieve of Atkin/; <https://en.wikipedia.org/wiki/Sieve_of_Atkin>.

	* <cr.yp.to/papers/primesieves-19990826.pdf>.

	* The implementation;
		has been optimised using a /wheel/ of static, but parameterised, size;
		has been parallelized;
		is polymorphic, but with a specialisation for type 'Int'.

 [@CAVEAT@] The 'Int'-specialisation is implemented by a /rewrite-rule/, which is /very/ fragile.
-}

module Factory.Math.Implementations.Primes.SieveOfAtkin(
-- * Types
-- ** Data-types
--	PolynomialType,
-- * Constants
--	atkinsModulus,
--	inherentPrimes,
--	nInherentPrimes,
--	squares,
-- * Functions
--	polynomialTypeLookupPeriod,
--	polynomialTypeLookup,
--	findPolynomialSolutions,
--	filterOddRepetitions,
--	generateMultiplesOfSquareTo,
--	getPrefactoredPrimes,
	sieveOfAtkin,
--	sieveOfAtkinInt
) 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

-- | Defines the types of /quadratic/, available to test the potential primality of a candidate integer.
data PolynomialType
	= ModFour	-- ^ Suitable for primality-testing numbers meeting @(n `mod` 4 == 1)@.
	| ModSix	-- ^ Suitable for primality-testing numbers meeting @(n `mod` 6 == 1)@.
	| ModTwelve	-- ^ Suitable for primality-testing numbers meeting @(n `mod` 12 == 11)@.
	| None		-- ^ There's no polynomial which can assess primality, because the candidate is composite.
	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

-- | The constant modulus used to select the appropriate quadratic for a prime candidate.
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]	-- Sure, this is always '12', but this is the reason why.

-- | The constant list of primes factored-out by the unoptimised algorithm.
inherentPrimes :: Integral i => [i]
inherentPrimes :: [i]
inherentPrimes	= [i
2, i
3]

-- | The constant number of primes factored-out by the unoptimised algorithm.
nInherentPrimes :: Int
nInherentPrimes :: Int
nInherentPrimes	= [Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Int]
forall i. Integral i => [i]
inherentPrimes :: [Int])

-- | Typically the set of primes which have been built into the specified /wheel/, but never fewer than 'inherentPrimes'.
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

-- | The period over which the data returned by 'polynomialTypeLookup' repeats.
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

{- |
	* Defines which, if any, of the three /quadratics/ is appropriate for the primality-test for each candidate.

	* Since this algorithm uses /modular arithmetic/, the /range/ of results repeat after a short /domain/ related to the /modulus/.
	Thus one need calculate at most one period of this cycle, but fewer if the maximum prime required falls within the first cycle of results.

	* Because the results are /bounded/, they're returned in a zero-indexed /array/, to provide efficient random access;
	the first few elements should never be required, but it makes query clearer.

	* <https://en.wikipedia.org/wiki/Sieve_of_Atkin>.
-}
polynomialTypeLookup :: (Data.Array.IArray.Ix i, Integral i)
	=> Data.PrimeWheel.PrimeWheel i
	-> i	-- ^ The maximum prime required.
	-> 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 :: Integral i => i -> PolynomialType
	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`)		-- Though this is merely /Trial Division/, it's only performed over a short bounded interval of numerators.
		) [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	-- We actually require @(n `mod` 4 == 1)@, but this is the equivalent modulo 12, with @(r == 9)@ removed because they're all divisible by /3/.
		| i
r i -> i -> Bool
forall a. Eq a => a -> a -> Bool
== i
7		= PolynomialType
ModSix	-- We actually require @(n `mod` 6 == 1)@, but this is the equivalent modulo 12, where @(r == 1)@ has been accounted for above.
		| i
r i -> i -> Bool
forall a. Eq a => a -> a -> Bool
== i
11		= PolynomialType
ModTwelve	-- We require @(n `mod` 12 == 11)@.
		| 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

-- | The constant, infinite list of the /squares/, of integers increasing from /1/.
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

{- |
	* Returns the /ordered/ list of those values with an /odd/ number of occurrences in the specified /unordered/ list.

	* CAVEAT: this is expensive in both execution-time and space.
	The typical imperative-style implementation accumulates polynomial-solutions in a /mutable array/ indexed by the candidate integer.
	This doesn't translate seamlessly to the /pure functional/ domain where /arrays/ naturally immutable,
	so we /sort/ a /list/ of polynomial-solutions, then measure the length of the solution-spans, corresponding to viable candidates.
	Regrettably, 'Data.List.sort' (implemented in /GHC/ by /mergesort/) has a time-complexity /O(n*log n)/
	which is greater than the theoretical /O(n)/ of the whole /Sieve of Atkin/;
	/GHC/'s old /qsort/-implementation is even slower :(
-}
filterOddRepetitions :: Ord a => [a] -> [a]
-- filterOddRepetitions	= map head . filter (foldr (const not) False) . Data.List.group . Data.List.sort	-- Too slow.
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]
_		= []

{- |
	* Returns the ordered list of solutions aggregated from each of three /bivariate quadratics/; @z = f(x, y)@.

	* For a candidate integer to be prime, it is necessary but insufficient, that there are an /odd/ number of solutions of value /candidate/.

	* At most one of these three polynomials is suitable for the validation of any specific candidate /z/, depending on 'lookupPolynomialType'.
	so the three sets of solutions are mutually exclusive.
	One coordinate @(x, y)@, can have solutions in more than one of the three polynomials.

	* This algorithm exhaustively traverses the domain @(x, y)@, for resulting /z/ of the required modulus.
	Whilst it tightly constrains the bounds of the search-space, it searches the domain methodically rather than intelligently.
-}
findPolynomialSolutions :: (Control.DeepSeq.NFData i, Data.Array.IArray.Ix i, Integral i)
	=> Data.PrimeWheel.PrimeWheel i
	-> i	-- ^ The maximum prime-number required.
	-> [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 {-The lists were previously sorted, as a side-effect, by 'filterOddRepetitions'-} ([[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
		], -- List-comprehension. Twice the length of the other two lists.
		{-# 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
		], -- List-comprehension.
		{-# 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
		] -- List-comprehension.
	] 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, oddSelection :: Integral i => [i]
		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	-- Effectively, those for meeting ((== 4) . (`mod` 6)).
			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	-- Effectively, those for meeting ((== 1) . (`mod` 6)).
			selection101 [a]
xs			= [a]
xs

--		lookupPolynomialType :: (Data.Array.IArray.Ix i, Integral i) => i -> PolynomialType
		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)

-- | Generates the /bounded/ list of multiples, of the /square/ of the specified prime, skipping those which aren't required.
generateMultiplesOfSquareTo :: Integral i
	=> Data.PrimeWheel.PrimeWheel i	-- ^ Used to generate the gaps between prime multiples of the square.
	-> i				-- ^ The /prime/.
	-> i				-- ^ The maximum bound.
	-> [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

{- |
	* Generates the constant /bounded/ list of /prime-numbers/.

	* <https://cr.yp.to/papers/primesieves-19990826.pdf>
-}
sieveOfAtkin :: (Control.DeepSeq.NFData i, Data.Array.IArray.Ix i, Integral i)
	=> Data.PrimeWheel.NPrimes	-- ^ Other implementations effectively use a hard-coded value either /2/ or /3/, but /6/ seems better.
	-> i				-- ^ The maximum prime required.
	-> [i]				-- ^ The /bounded/ list of primes.
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 :: Integral i => Data.Set.Set i -> [i] -> [i]
	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	-- Tail-recurse.
		| 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 #-}	-- CAVEAT: doesn't fire when built with profiling enabled.

-- | A specialisation of 'sieveOfAtkin', which reduces both the execution-time and the space required.
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