{-
	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 <http://www.gnu.org/licenses/>.
-}
{- |
 [@AUTHOR@]	Dr. Alistair Ward

 [@DESCRIPTION@]

	* Generates the constant, conceptually infinite, list of /prime-numbers/, using the /Sieve of Eratosthenes/; <https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes>.

	* Based on <http://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf>.

	* The implementation;
		has been optimised using a /wheel/ of static, but parameterised, size;
		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.SieveOfEratosthenes(
-- * Types
-- ** Type-synonyms
--	PrimeMultiplesQueue,
--	PrimeMultiplesMap,
--	Repository,
--	PrimeMultiplesMapInt,
--	RepositoryInt,
-- * Functions
--	head',
--	tail',
	sieveOfEratosthenes,
--	sieveOfEratosthenesInt
) where

import			Control.Arrow((&&&), (***))
import qualified	Control.Arrow
import qualified	Data.IntMap
import qualified	Data.Map
import			Data.Sequence((|>))
import qualified	Data.Sequence
import qualified	Factory.Data.PrimeWheel		as Data.PrimeWheel

-- | The 'Data.Sequence.Seq' counterpart to 'Data.List.head'.
head' :: Data.Sequence.Seq [a] -> [a]
head' :: Seq [a] -> [a]
head'	= (Seq [a] -> Int -> [a]
forall a. Seq a -> Int -> a
`Data.Sequence.index` Int
0)

{- |
	* The 'Data.Sequence.Seq' counterpart to 'Data.List.tail'.

	* CAVEAT: because @ Data.List.tail [] @ returns an error, whereas @ tail' Data.Sequence.empty @ returns 'Data.Sequence.empty',
	this function is for internal use only.
-}
tail' :: Data.Sequence.Seq [a] -> Data.Sequence.Seq [a]
tail' :: Seq [a] -> Seq [a]
tail'	= Int -> Seq [a] -> Seq [a]
forall a. Int -> Seq a -> Seq a
Data.Sequence.drop Int
1

-- | An ordered queue of the multiples of primes.
type PrimeMultiplesQueue i	= Data.Sequence.Seq (Data.PrimeWheel.PrimeMultiples i)

-- | A map of the multiples of primes.
type PrimeMultiplesMap i	= Data.Map.Map i (Data.PrimeWheel.PrimeMultiples i)

-- | Combine a /queue/, with a /map/, to form a repository to hold prime-multiples.
type Repository i	= (PrimeMultiplesQueue i, PrimeMultiplesMap i)

{- |
	* A refinement of the /Sieve Of Eratosthenes/, which pre-sieves candidates, selecting only those /coprime/ to the specified short sequence of low prime-numbers.

	* The short sequence of initial primes are represented by a 'Data.PrimeWheel.PrimeWheel',
	of parameterised, but static, size; <https://en.wikipedia.org/wiki/Wheel_factorization>.

	* The algorithm requires one to record multiples of previously discovered primes, allowing /composite/ candidates to be eliminated by comparison.

	* Because each /list/ of multiples, starts with the /square/ of the prime from which it was generated,
	the vast majority will be larger than the maximum prime ultimately demanded, and the effort of constructing and storing this list, is consequently wasted.
	Many implementations solve this, by requiring specification of the maximum prime required,
	thus allowing the construction of redundant lists of multiples to be avoided.

	* This implementation doesn't impose that constraint, leaving a requirement for /rapid/ storage,
	which is supported by /appending/ the /list/ of prime-multiples, to a /queue/.
	If a large enough candidate is ever generated, to match the /head/ of the /list/ of prime-multiples,
	at the /head/ of this /queue/, then the whole /list/ of prime-multiples is dropped from the /queue/,
	but the /tail/ of this /list/ of prime-multiples, for which there is now a high likelyhood of a subsequent match, must now be re-recorded.
	A /queue/ doesn't support efficient random /insertion/, so a 'Data.Map.Map' is used for these subsequent multiples.
	This solution is faster than just using a "Data.PQueue.Min".

	* CAVEAT: has linear /O(n)/ space-complexity.
-}
sieveOfEratosthenes :: Integral i
	=> Data.PrimeWheel.NPrimes
	-> [i]
sieveOfEratosthenes :: Int -> [i]
sieveOfEratosthenes	= ([i] -> [i] -> [i]) -> ([i], [i]) -> [i]
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry [i] -> [i] -> [i]
forall a. [a] -> [a] -> [a]
(++) (([i], [i]) -> [i]) -> (Int -> ([i], [i])) -> Int -> [i]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (PrimeWheel i -> [i]
forall i. PrimeWheel i -> [i]
Data.PrimeWheel.getPrimeComponents (PrimeWheel i -> [i])
-> (PrimeWheel i -> [i]) -> PrimeWheel i -> ([i], [i])
forall (a :: * -> * -> *) b c c'.
Arrow a =>
a b c -> a b c' -> a b (c, c')
&&& [Distance i] -> [i]
forall i. Integral i => [Distance i] -> [i]
start ([Distance i] -> [i])
-> (PrimeWheel i -> [Distance i]) -> PrimeWheel i -> [i]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. PrimeWheel i -> [Distance i]
forall i. Integral i => PrimeWheel i -> [Distance i]
Data.PrimeWheel.roll) (PrimeWheel i -> ([i], [i]))
-> (Int -> PrimeWheel i) -> Int -> ([i], [i])
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> PrimeWheel i
forall i. Integral i => Int -> PrimeWheel i
Data.PrimeWheel.mkPrimeWheel	where
	start :: Integral i => [Data.PrimeWheel.Distance i] -> [i]
	start :: [Distance i] -> [i]
start ~((i
candidate, [i]
rollingWheel) : [Distance i]
distances)	= i
candidate i -> [i] -> [i]
forall a. a -> [a] -> [a]
: Distance i -> Repository i -> [i]
forall i. Integral i => Distance i -> Repository i -> [i]
sieve ([Distance i] -> Distance i
forall a. [a] -> a
head [Distance i]
distances) ([i] -> Seq [i]
forall a. a -> Seq a
Data.Sequence.singleton ([i] -> Seq [i]) -> [i] -> Seq [i]
forall a b. (a -> b) -> a -> b
$ i -> [i] -> [i]
forall i. Integral i => i -> [i] -> [i]
Data.PrimeWheel.generateMultiples i
candidate [i]
rollingWheel, Map i [i]
forall k a. Map k a
Data.Map.empty)

	sieve :: Integral i => Data.PrimeWheel.Distance i -> Repository i -> [i]
	sieve :: Distance i -> Repository i -> [i]
sieve distance :: Distance i
distance@(i
candidate, [i]
rollingWheel) repository :: Repository i
repository@(PrimeMultiplesQueue i
primeSquares, PrimeMultiplesMap i
squareFreePrimeMultiples)	= case i -> PrimeMultiplesMap i -> Maybe [i]
forall k a. Ord k => k -> Map k a -> Maybe a
Data.Map.lookup i
candidate PrimeMultiplesMap i
squareFreePrimeMultiples of
		Just [i]
primeMultiples	-> Repository i -> [i]
sieve' (Repository i -> [i]) -> Repository i -> [i]
forall a b. (a -> b) -> a -> b
$ (PrimeMultiplesMap i -> PrimeMultiplesMap i)
-> Repository i -> Repository i
forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (d, b) (d, c)
Control.Arrow.second ([i] -> PrimeMultiplesMap i -> PrimeMultiplesMap i
forall i.
Ord i =>
PrimeMultiples i -> PrimeMultiplesMap i -> PrimeMultiplesMap i
insertUniq [i]
primeMultiples (PrimeMultiplesMap i -> PrimeMultiplesMap i)
-> (PrimeMultiplesMap i -> PrimeMultiplesMap i)
-> PrimeMultiplesMap i
-> PrimeMultiplesMap i
forall b c a. (b -> c) -> (a -> b) -> a -> c
. i -> PrimeMultiplesMap i -> PrimeMultiplesMap i
forall k a. Ord k => k -> Map k a -> Map k a
Data.Map.delete i
candidate) Repository i
repository	-- Re-insert subsequent multiples.
		Maybe [i]
Nothing -- Not a square-free composite.
			| i
candidate i -> i -> Bool
forall a. Eq a => a -> a -> Bool
== i
smallestPrimeSquare	-> Repository i -> [i]
sieve' (Repository i -> [i]) -> Repository i -> [i]
forall a b. (a -> b) -> a -> b
$ (PrimeMultiplesQueue i -> PrimeMultiplesQueue i
forall a. Seq [a] -> Seq [a]
tail' (PrimeMultiplesQueue i -> PrimeMultiplesQueue i)
-> (PrimeMultiplesMap i -> PrimeMultiplesMap i)
-> Repository i
-> Repository i
forall (a :: * -> * -> *) b c b' c'.
Arrow a =>
a b c -> a b' c' -> a (b, b') (c, c')
*** [i] -> PrimeMultiplesMap i -> PrimeMultiplesMap i
forall i.
Ord i =>
PrimeMultiples i -> PrimeMultiplesMap i -> PrimeMultiplesMap i
insertUniq [i]
subsequentPrimeMultiples) Repository i
repository	-- Migrate subsequent prime-multiples, from 'primeSquares' to 'squareFreePrimeMultiples'.
			| Bool
otherwise {-prime-}			-> i
candidate i -> [i] -> [i]
forall a. a -> [a] -> [a]
: Repository i -> [i]
sieve' ((PrimeMultiplesQueue i -> PrimeMultiplesQueue i)
-> Repository i -> Repository i
forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (b, d) (c, d)
Control.Arrow.first (PrimeMultiplesQueue i -> [i] -> PrimeMultiplesQueue i
forall a. Seq a -> a -> Seq a
|> i -> [i] -> [i]
forall i. Integral i => i -> [i] -> [i]
Data.PrimeWheel.generateMultiples i
candidate [i]
rollingWheel) Repository i
repository)
			where
				(i
smallestPrimeSquare : [i]
subsequentPrimeMultiples)	= PrimeMultiplesQueue i -> [i]
forall a. Seq [a] -> [a]
head' PrimeMultiplesQueue i
primeSquares
		where
--			sieve' :: Repository i -> [i]
			sieve' :: Repository i -> [i]
sieve'	= Distance i -> Repository i -> [i]
forall i. Integral i => Distance i -> Repository i -> [i]
sieve (Distance i -> Repository i -> [i])
-> Distance i -> Repository i -> [i]
forall a b. (a -> b) -> a -> b
$ Distance i -> Distance i
forall i. Integral i => Distance i -> Distance i
Data.PrimeWheel.rotate Distance i
distance	-- Tail-recurse.

			insertUniq :: Ord i => Data.PrimeWheel.PrimeMultiples i -> PrimeMultiplesMap i -> PrimeMultiplesMap i
			insertUniq :: PrimeMultiples i -> PrimeMultiplesMap i -> PrimeMultiplesMap i
insertUniq PrimeMultiples i
l PrimeMultiplesMap i
m	= PrimeMultiples i -> PrimeMultiplesMap i
insert (PrimeMultiples i -> PrimeMultiplesMap i)
-> PrimeMultiples i -> PrimeMultiplesMap i
forall a b. (a -> b) -> a -> b
$ (i -> Bool) -> PrimeMultiples i -> PrimeMultiples i
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (i -> PrimeMultiplesMap i -> Bool
forall k a. Ord k => k -> Map k a -> Bool
`Data.Map.member` PrimeMultiplesMap i
m) PrimeMultiples i
l	where
--				insert :: Ord i => Data.PrimeWheel.PrimeMultiples i -> PrimeMultiplesMap i
				insert :: PrimeMultiples i -> PrimeMultiplesMap i
insert []		= [Char] -> PrimeMultiplesMap i
forall a. HasCallStack => [Char] -> a
error [Char]
"Factory.Math.Implementations.Primes.SieveOfEratosthenes.sieveOfEratosthenes.sieve.insertUniq.insert:\tnull list"
				insert (i
key : PrimeMultiples i
values)	= i -> PrimeMultiples i -> PrimeMultiplesMap i -> PrimeMultiplesMap i
forall k a. Ord k => k -> a -> Map k a -> Map k a
Data.Map.insert i
key PrimeMultiples i
values PrimeMultiplesMap i
m

{-# NOINLINE sieveOfEratosthenes #-}
{-# RULES "sieveOfEratosthenes/Int" sieveOfEratosthenes = sieveOfEratosthenesInt #-}	-- CAVEAT: doesn't fire when built with profiling enabled.

-- | A specialisation of 'PrimeMultiplesMap'.
type PrimeMultiplesMapInt	= Data.IntMap.IntMap (Data.PrimeWheel.PrimeMultiples Int)

-- | A specialisation of 'Repository'.
type RepositoryInt	= (PrimeMultiplesQueue Int, PrimeMultiplesMapInt)

{- |
	* A specialisation of 'sieveOfEratosthenes', which approximately /doubles/ the speed and reduces the space required.

	* CAVEAT: because the algorithm involves /squares/ of primes,
	this implementation will overflow when finding primes greater than @2^16@ on a /32-bit/ machine.
-}
sieveOfEratosthenesInt :: Data.PrimeWheel.NPrimes -> [Int]
sieveOfEratosthenesInt :: Int -> [Int]
sieveOfEratosthenesInt	= ([Int] -> [Int] -> [Int]) -> ([Int], [Int]) -> [Int]
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry [Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
(++) (([Int], [Int]) -> [Int])
-> (Int -> ([Int], [Int])) -> Int -> [Int]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (PrimeWheel Int -> [Int]
forall i. PrimeWheel i -> [i]
Data.PrimeWheel.getPrimeComponents (PrimeWheel Int -> [Int])
-> (PrimeWheel Int -> [Int]) -> PrimeWheel Int -> ([Int], [Int])
forall (a :: * -> * -> *) b c c'.
Arrow a =>
a b c -> a b c' -> a b (c, c')
&&& [Distance Int] -> [Int]
start ([Distance Int] -> [Int])
-> (PrimeWheel Int -> [Distance Int]) -> PrimeWheel Int -> [Int]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. PrimeWheel Int -> [Distance Int]
forall i. Integral i => PrimeWheel i -> [Distance i]
Data.PrimeWheel.roll) (PrimeWheel Int -> ([Int], [Int]))
-> (Int -> PrimeWheel Int) -> Int -> ([Int], [Int])
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> PrimeWheel Int
forall i. Integral i => Int -> PrimeWheel i
Data.PrimeWheel.mkPrimeWheel	where
	start :: [Data.PrimeWheel.Distance Int] -> [Int]
	start :: [Distance Int] -> [Int]
start ~((Int
candidate, [Int]
rollingWheel) : [Distance Int]
distances)	= Int
candidate Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: Distance Int -> RepositoryInt -> [Int]
sieve ([Distance Int] -> Distance Int
forall a. [a] -> a
head [Distance Int]
distances) ([Int] -> Seq [Int]
forall a. a -> Seq a
Data.Sequence.singleton ([Int] -> Seq [Int]) -> [Int] -> Seq [Int]
forall a b. (a -> b) -> a -> b
$ Int -> [Int] -> [Int]
forall i. Integral i => i -> [i] -> [i]
Data.PrimeWheel.generateMultiples Int
candidate [Int]
rollingWheel, IntMap [Int]
forall a. IntMap a
Data.IntMap.empty)

	sieve :: Data.PrimeWheel.Distance Int -> RepositoryInt -> [Int]
	sieve :: Distance Int -> RepositoryInt -> [Int]
sieve distance :: Distance Int
distance@(Int
candidate, [Int]
rollingWheel) repository :: RepositoryInt
repository@(Seq [Int]
primeSquares, IntMap [Int]
squareFreePrimeMultiples)	= case Int -> IntMap [Int] -> Maybe [Int]
forall a. Int -> IntMap a -> Maybe a
Data.IntMap.lookup Int
candidate IntMap [Int]
squareFreePrimeMultiples of
		Just [Int]
primeMultiples	-> RepositoryInt -> [Int]
sieve' (RepositoryInt -> [Int]) -> RepositoryInt -> [Int]
forall a b. (a -> b) -> a -> b
$ (IntMap [Int] -> IntMap [Int]) -> RepositoryInt -> RepositoryInt
forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (d, b) (d, c)
Control.Arrow.second ([Int] -> IntMap [Int] -> IntMap [Int]
insertUniq [Int]
primeMultiples (IntMap [Int] -> IntMap [Int])
-> (IntMap [Int] -> IntMap [Int]) -> IntMap [Int] -> IntMap [Int]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> IntMap [Int] -> IntMap [Int]
forall a. Int -> IntMap a -> IntMap a
Data.IntMap.delete Int
candidate) RepositoryInt
repository
		Maybe [Int]
Nothing
			| Int
candidate Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
smallestPrimeSquare	-> RepositoryInt -> [Int]
sieve' (RepositoryInt -> [Int]) -> RepositoryInt -> [Int]
forall a b. (a -> b) -> a -> b
$ (Seq [Int] -> Seq [Int]
forall a. Seq [a] -> Seq [a]
tail' (Seq [Int] -> Seq [Int])
-> (IntMap [Int] -> IntMap [Int]) -> RepositoryInt -> RepositoryInt
forall (a :: * -> * -> *) b c b' c'.
Arrow a =>
a b c -> a b' c' -> a (b, b') (c, c')
*** [Int] -> IntMap [Int] -> IntMap [Int]
insertUniq [Int]
subsequentPrimeMultiples) RepositoryInt
repository
			| Bool
otherwise				-> Int
candidate Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: RepositoryInt -> [Int]
sieve' ((Seq [Int] -> Seq [Int]) -> RepositoryInt -> RepositoryInt
forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (b, d) (c, d)
Control.Arrow.first (Seq [Int] -> [Int] -> Seq [Int]
forall a. Seq a -> a -> Seq a
|> Int -> [Int] -> [Int]
forall i. Integral i => i -> [i] -> [i]
Data.PrimeWheel.generateMultiples Int
candidate [Int]
rollingWheel) RepositoryInt
repository)
			where
				(Int
smallestPrimeSquare : [Int]
subsequentPrimeMultiples)	= Seq [Int] -> [Int]
forall a. Seq [a] -> [a]
head' Seq [Int]
primeSquares
		where
			sieve' :: RepositoryInt -> [Int]
			sieve' :: RepositoryInt -> [Int]
sieve'	= Distance Int -> RepositoryInt -> [Int]
sieve (Distance Int -> RepositoryInt -> [Int])
-> Distance Int -> RepositoryInt -> [Int]
forall a b. (a -> b) -> a -> b
$ Distance Int -> Distance Int
forall i. Integral i => Distance i -> Distance i
Data.PrimeWheel.rotate Distance Int
distance

			insertUniq :: Data.PrimeWheel.PrimeMultiples Int -> PrimeMultiplesMapInt -> PrimeMultiplesMapInt
			insertUniq :: [Int] -> IntMap [Int] -> IntMap [Int]
insertUniq [Int]
l IntMap [Int]
m	= [Int] -> IntMap [Int]
insert ([Int] -> IntMap [Int]) -> [Int] -> IntMap [Int]
forall a b. (a -> b) -> a -> b
$ (Int -> Bool) -> [Int] -> [Int]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (Int -> IntMap [Int] -> Bool
forall a. Int -> IntMap a -> Bool
`Data.IntMap.member` IntMap [Int]
m) [Int]
l	where
				insert :: Data.PrimeWheel.PrimeMultiples Int -> PrimeMultiplesMapInt
				insert :: [Int] -> IntMap [Int]
insert []		= [Char] -> IntMap [Int]
forall a. HasCallStack => [Char] -> a
error [Char]
"Factory.Math.Implementations.Primes.SieveOfEratosthenes.sieveOfEratosthenesInt.sieve.insertUniq.insert:\tnull list"
				insert (Int
key : [Int]
values)	= Int -> [Int] -> IntMap [Int] -> IntMap [Int]
forall a. Int -> a -> IntMap a -> IntMap a
Data.IntMap.insert Int
key [Int]
values IntMap [Int]
m