{-
	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/ by a variety of different algorithms.

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

	* <http://www.haskell.org/haskellwiki/Prime_numbers>.

	* <http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.31.3936&rep=rep1&type=pdf>.

	* <http://larc.unt.edu/ian/pubs/sieve.pdf>.
-}

module Factory.Math.Implementations.Primes(
-- * Types
-- ** Type-synonyms
--	PrimeMultiplesQueue,
--	PrimeMultiplesMap,
--	Repository,
--	PrimeMultiplesMapInt,
--	RepositoryInt,
-- ** Data-types
	Algorithm(..),
-- * Functions
--	head',
--	tail',
--	turnersSieve,
--	trialDivision,
--	sieveOfEratosthenes,
--	sieveOfEratosthenesInt,
-- ** Predicates
--	isIndivisible
) where

import			Control.Arrow((&&&), (***))
import qualified	Control.Arrow
import qualified	Data.IntMap
import qualified	Data.List
import qualified	Data.Map
import qualified	Data.Numbers.Primes
import			Data.Sequence((|>))
import qualified	Data.Sequence
import qualified	Factory.Data.PrimeWheel		as Data.PrimeWheel
import qualified	Factory.Math.Power		as Math.Power
import qualified	Factory.Math.PrimeFactorisation	as Math.PrimeFactorisation
import qualified	Factory.Math.Primes		as Math.Primes
import qualified	ToolShed.Defaultable		as Defaultable

-- | The implemented methods by which the primes may be generated.
data Algorithm
	= TurnersSieve			-- ^ For each /prime/, the infinite list of candidates greater than its /square/, is filtered for indivisibility; <http://www.haskell.org/haskellwiki/Prime_numbers#Turner.27s_sieve_-_Trial_division>.
	| TrialDivision Int		-- ^ For each candidate, confirm indivisibility, by all /primes/ smaller than its /square-root/, optimised using a 'Data.PrimeWheel.PrimeWheel' (<http://en.wikipedia.org/wiki/Wheel_factorization>).
	| SieveOfEratosthenes Int	-- ^ The /Sieve of Eratosthenes/ (<http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes>), optimised using a 'Data.PrimeWheel.PrimeWheel' (<http://en.wikipedia.org/wiki/Wheel_factorization>).
	| WheelSieve Int		-- ^ 'Data.Numbers.Primes.wheelSieve'.
	deriving (Eq, Read, Show)

instance Defaultable.Defaultable Algorithm	where
	defaultValue	= SieveOfEratosthenes 7	--Resulting in wheel-circumference=510510.

instance Math.Primes.Algorithmic Algorithm	where
	primes TurnersSieve		= turnersSieve
	primes (TrialDivision n)	= trialDivision n
	primes (SieveOfEratosthenes n)	= sieveOfEratosthenes n			--When (n == 0), this degenerates to the unoptimised classic form.
	primes (WheelSieve n)		= Data.Numbers.Primes.wheelSieve n	--Has better space-complexity than 'SieveOfEratosthenes'.

-- | Uses /Trial Division/, to determine whether the specified numerator is indivisible by all the specified denominators.
isIndivisible :: Integral i => i -> [i] -> Bool
isIndivisible numerator	= all ((/= 0) . (numerator `mod`))

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

-- | The 'Data.Sequence.Seq' counterpart to 'Data.List.tail'.
tail' :: Data.Sequence.Seq [a] -> Data.Sequence.Seq [a]
tail'	= Data.Sequence.drop 1

{- |
	* Generates the constant, conceptually infinite, list of /prime-numbers/.

	* For each /prime/, the infinite list of candidates greater than its /square/,
	is filtered for indivisibility; <http://www.haskell.org/haskellwiki/Prime_numbers#Turner.27s_sieve_-_Trial_division>.
-}
turnersSieve :: Integral prime => [prime]
turnersSieve	= 2 : sieve [3, 5 ..]	where
	sieve :: Integral i => [i] -> [i]
	sieve []			= []
	sieve (prime : candidates)	= prime : sieve (
		filter (
			\candidate	-> any ($ candidate) [
				(< Math.Power.square prime),	--Unconditionally admit any candidate smaller than the square of the last prime.
				(/= 0) . (`mod` prime)		--Ensure indivisibility, of all subsequent candidates, by the last prime.
			]
		) candidates
	 )

{- |
	* Generates the constant, conceptually infinite, list of /prime-numbers/.

	* For each candidate, confirm indivisibility, by all /primes/ smaller than its /square-root/.

	* The candidates to sieve, are generated by a 'Data.PrimeWheel.PrimeWheel',
	of parameterised, but static, size; <http://en.wikipedia.org/wiki/Wheel_factorization>.
-}
trialDivision :: Integral prime => Int -> [prime]
trialDivision n	= Data.PrimeWheel.getPrimeComponents primeWheel ++ indivisible where
	primeWheel	= Data.PrimeWheel.mkPrimeWheel n
	candidates	= map fst $ Data.PrimeWheel.roll primeWheel
	indivisible	= uncurry (++) . Control.Arrow.second (
--		filter (\candidate -> isIndivisible candidate . zipWith const indivisible . takeWhile (<= candidate) $ map Math.Power.square indivisible)
		filter (\candidate -> isIndivisible candidate $ takeWhile (<= Math.PrimeFactorisation.maxBoundPrimeFactor candidate) indivisible {-recurse-})
	 ) $ Data.List.span (
		< Math.Power.square (head candidates)	--The first composite candidate, is the square of the next prime after the wheel's constituent ones.
	 ) candidates

{-# NOINLINE trialDivision #-}	--Required to prevent optimization prior to firing of rewrite-rule ?!

-- | 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; <http://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 required, and the effort of constructing and storing this list, is 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 the /queue/, then the whole /list/ of prime-multiples is dropped,
	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 the same algorithm using "Data.PQueue.Min".

	* CAVEAT: has linear /O(primes)/ space-complexity.
-}
sieveOfEratosthenes :: Integral i => Int -> [i]
sieveOfEratosthenes	= uncurry (++) . (Data.PrimeWheel.getPrimeComponents &&& start . Data.PrimeWheel.roll) . Data.PrimeWheel.mkPrimeWheel	where
	start :: Integral i => [Data.PrimeWheel.Distance i] -> [i]
	start ~((candidate, rollingWheel) : distances)	= candidate : sieve (head distances) (Data.Sequence.singleton $ Data.PrimeWheel.generatePrimeMultiples candidate rollingWheel, Data.Map.empty)

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

			insertUniq :: Ord i => Data.PrimeWheel.PrimeMultiples i -> PrimeMultiplesMap i -> PrimeMultiplesMap i
			insertUniq l m	= insert (dropWhile (`Data.Map.member` m) l)	where
--				insert :: Ord i => Data.PrimeWheel.PrimeMultiples i -> PrimeMultiplesMap i
				insert []		= error "Factory.Math.Implementations.Primes.sieveOfEratosthenes.sieve.insertUniq.insert:\tnull list"
				insert (key : values)	= Data.Map.insert key values 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.

	* CAVEAT: because the algorithm involves /squares/ of primes,
	this implementation will overflow when finding primes greater than @ 2^16 @ on a /32-bit/ machine;
	but it will exhaust the memory before that anyway.
-}
sieveOfEratosthenesInt :: Int -> [Int]
sieveOfEratosthenesInt	= uncurry (++) . (Data.PrimeWheel.getPrimeComponents &&& start . Data.PrimeWheel.roll) . Data.PrimeWheel.mkPrimeWheel	where
	start :: [Data.PrimeWheel.Distance Int] -> [Int]
	start ~((candidate, rollingWheel) : distances)	= candidate : sieve (head distances) (Data.Sequence.singleton $ Data.PrimeWheel.generatePrimeMultiples candidate rollingWheel, Data.IntMap.empty)

	sieve :: Data.PrimeWheel.Distance Int -> RepositoryInt -> [Int]
	sieve distance@(candidate, rollingWheel) repository@(primeSquares, squareFreePrimeMultiples)	= case Data.IntMap.lookup candidate squareFreePrimeMultiples of
		Just primeMultiples	-> sieve' $ Control.Arrow.second (insertUniq primeMultiples . Data.IntMap.delete candidate) repository
		Nothing
			| candidate == smallestPrimeSquare	-> sieve' $ (tail' *** insertUniq subsequentPrimeMultiples) repository
			| otherwise				-> candidate : sieve' (Control.Arrow.first (|> Data.PrimeWheel.generatePrimeMultiples candidate rollingWheel) repository)
			where
				(smallestPrimeSquare : subsequentPrimeMultiples)	= head' primeSquares
		where
			sieve' :: RepositoryInt -> [Int]
			sieve'	= sieve $ Data.PrimeWheel.rotate distance

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