{-
	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@]	Defines a /prime-wheel/, for use in prime-number generation; <https://en.wikipedia.org/wiki/Wheel_factorization>.
-}

module Factory.Data.PrimeWheel(
-- * Types
-- ** Type-synonyms
	Distance,
	NPrimes,
	PrimeMultiples,
--	Repository,
-- ** Data-types
	PrimeWheel(getPrimeComponents, getSpokeGaps),
-- * Functions
	estimateOptimalSize,
--	findCoprimes,
	generateMultiples,
	roll,
	rotate,
-- ** Constructor
	mkPrimeWheel,
-- ** Query
	getCircumference,
	getSpokeCount
) where

import			Control.Arrow((&&&), (***))
import qualified	Data.IntMap
import qualified	Data.List

{- |
	* A conceptual /wheel/, with irregularly spaced spokes; <http://www.haskell.org/haskellwiki/Prime_numbers_miscellaneous#Prime_Wheels>.

	* On being rolled, the trace of the spokes, identifies candidates which are /coprime/ to those primes from which the /wheel/ was composed.

	* One can alternatively view this as a set of vertical nested rings, each with a /prime circumference/, and touching at its lowest point.
	Each has a single mark on its /circumference/, which when rolled identifies multiples of that /circumference/.
	When the complete set is rolled, from the state where all marks are coincident, all multiples of the set of primes, are traced.

	* CAVEAT: The distance required to return to this state (the wheel's /circumference/), grows rapidly with the number of primes:

>	zip [0 ..] . scanl (*) 1 $ [2,3,5,7,11,13,17,19,23,29,31]
>	[(0,1),(1,2),(2,6),(3,30),(4,210),(5,2310),(6,30030),(7,510510),(8,9699690),(9,223092870),(10,6469693230),(11,200560490130)]

	* The number of spokes also grows rapidly with the number of primes:

>	zip [0 ..] . scanl (*) 1 . map pred $ [2,3,5,7,11,13,17,19,23,29,31]
>	[(0,1),(1,1),(2,2),(3,8),(4,48),(5,480),(6,5760),(7,92160),(8,1658880),(9,36495360),(10,1021870080),(11,30656102400)]
-}
data PrimeWheel i	= MkPrimeWheel {
	getPrimeComponents	:: [i],	-- ^ Accessor: the ordered sequence of initial primes, from which the /wheel/ was composed.
	getSpokeGaps		:: [i]	-- ^ Accessor: the sequence of spoke-gaps, the sum of which equals its /circumference/.
} deriving Show

-- | The /circumference/ of the specified 'PrimeWheel'.
getCircumference :: Integral i => PrimeWheel i -> i
getCircumference	= product . getPrimeComponents

-- | The number of spokes in the specified 'PrimeWheel'.
getSpokeCount :: Integral i => PrimeWheel i -> i
getSpokeCount	= foldr ((*) . pred) 1 . getPrimeComponents

-- | An infinite increasing sequence, of the multiples of a specific prime.
type PrimeMultiples i	= [i]

-- | Defines a container for the 'PrimeMultiples'.
type Repository	= Data.IntMap.IntMap (PrimeMultiples Int)

-- | The size of the /wheel/, measured by the number of primes from which it is composed.
type NPrimes	= Int

{- |
	* Uses a /Sieve of Eratosthenes/ (<https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes>), to generate an initial sequence of primes.

	* Also generates an infinite sequence of candidate primes, each of which is /coprime/ to the primes just found, e.g.:
	@filter ((== 1) . (gcd (2 * 3 * 5 * 7))) [11 ..] = [11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,121 ..]@; NB /121/ isn't prime.

	* CAVEAT: the use, for efficiency, of "Data.IntMap", limits the maximum bound of this sequence, though not to a significant extent.
-}
findCoprimes :: NPrimes -> ([Int], [Int])
findCoprimes 0	= ([], [])
findCoprimes required
	| required < 0	= error $ "Factory.Data.PrimeWheel.findCoprimes: invalid number of coprimes; " ++ show required
	| otherwise	= splitAt required $ 2 : sieve 3 0 Data.IntMap.empty
	where
		sieve :: Int -> NPrimes -> Repository -> [Int]
		sieve candidate found repository	= case Data.IntMap.lookup candidate repository of
			Just primeMultiples	-> sieve' found . insertUniq primeMultiples $ Data.IntMap.delete candidate repository	-- Re-insert subsequent multiples.
			Nothing {-prime-}	-> let
				found'		= succ found
				(key : values)	= iterate (+ gap * candidate) $ candidate ^ (2 :: Int)	-- Generate a sequence of prime-multiples, starting from its square.
			 in candidate : sieve' found' (
				if found' >= required
					then repository
					else Data.IntMap.insert key values repository
			 )
			where
				gap :: Int
				gap	= 2	-- For efficiency, only sieve odd integers.

				sieve' :: NPrimes -> Repository -> [Int]
				sieve'	= sieve $ candidate + gap	-- Tail-recurse.

				insertUniq :: PrimeMultiples Int -> Repository -> Repository
				insertUniq l m	= insert $ dropWhile (`Data.IntMap.member` m) l	where
					insert :: PrimeMultiples Int -> Repository
					insert []		= error "Factory.Data.PrimeWheel.findCoprimes.sieve.insertUniq.insert:\tnull list"
					insert (key : values)	= Data.IntMap.insert key values m
{- |
	* The optimal number of low primes from which to build the /wheel/, grows with the number of primes required;
	the /circumference/ should be approximately the /square-root/ of the number of integers it will be required to sieve.

	* CAVEAT: one greater than this is returned, which empirically seems better.
-}
estimateOptimalSize :: Integral i => i -> NPrimes
estimateOptimalSize maxPrime	= succ . length . takeWhile (<= optimalCircumference) . scanl1 (*) {-circumference-} . map fromIntegral {-prevent overflow-} . fst {-primes-} $ findCoprimes 10 {-arbitrary maximum bound-}	where
	optimalCircumference :: Integer
	optimalCircumference	= round (sqrt $ fromIntegral maxPrime :: Double)

{- |
	Smart constructor for a /wheel/ from the specified number of low primes.

	* The optimal number of low primes from which to build the /wheel/, grows with the number of primes required;
	the /circumference/ should be approximately the /square-root/ of the number of integers it will be required to sieve.

	* The sequence of gaps between spokes on the /wheel/ is /symmetrical under reflection/;
	though two values lie /on/ the axis, that aren't part of this symmetry. Eg:

>	nPrimes	Gaps
>	======	====
>	0	[1]
>	1	[2]	-- The terminal gap for all subsequent wheels is '2'; [(succ circumference `mod` circumference) - (pred circumference `mod` circumference)].
>	2	[4,2]	-- Both points are on the axis, so the symmetry isn't yet clear.
>	3	[6,4,2,4,2,4,6,2]
>	4	[10,2,4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,8,6,4,6,2,4,6,2,6,6,4,2,4,6,2,6,4,2,4,2,10,2]

	Exploitation of this property has proved counter-productive, probably because it requires /strict evaluation/,
	exposing the user to the full cost of inadvertently choosing a /wheel/, which in practice, is rotated less than once.
-}
mkPrimeWheel :: Integral i => NPrimes -> PrimeWheel i
mkPrimeWheel 0	= MkPrimeWheel [] [1]
mkPrimeWheel nPrimes
	| nPrimes < 0	= error $ "Factory.Data.PrimeWheel.mkPrimeWheel: unable to construct from " ++ show nPrimes ++ " primes"
	| otherwise	= primeWheel
	where
		(primeComponents, coprimeCandidates)	= (map fromIntegral *** map fromIntegral . Data.List.genericTake (getSpokeCount primeWheel)) $ findCoprimes nPrimes
		primeWheel				= MkPrimeWheel primeComponents $ zipWith (-) coprimeCandidates $ 1 : coprimeCandidates	-- Measure the gaps between candidate primes.

-- | Couples a candidate prime with a /rolling wheel/, to define the distance rolled.
type Distance i	= (i, [i])

-- | Generates a new candidate prime, from a /rolling wheel/, and the current candidate.
rotate :: Integral i => Distance i -> Distance i
rotate (candidate, rollingWheel)	= (candidate +) . head &&& tail $ rollingWheel

{-# INLINE rotate #-}

-- | Generate an infinite, increasing sequence of candidate primes, from the specified /wheel/.
roll :: Integral i => PrimeWheel i -> [Distance i]
roll primeWheel	= tail $ iterate rotate (1, cycle $ getSpokeGaps primeWheel)

{- |
	* Generates multiples of the specified prime, starting from its /square/,
	skipping those multiples of the low primes from which the specified 'PrimeWheel' was composed,
	and which therefore, the /wheel/ won't generate as candidates. Eg:

>	Prime	Rotating PrimeWheel 3	Output
>	=====	=====================	======
>	7	[4,2,4,2,4,6,2,6]	[49,77,91,119,133,161,203,217,259 ..]
>	11	[2,4,2,4,6,2,6,4]	[121,143,187,209,253,319,341,407 ..]
>	13	[4,2,4,6,2,6,4,2]	[169,221,247,299,377,403,481,533,559 ..]
-}
generateMultiples :: Integral i
	=> i	-- ^ The number to square and multiply
	-> [i]	-- ^ A /rolling wheel/, the track of which, delimits the gaps between /coprime/ candidates.
	-> [i]
generateMultiples i	= scanl (\accumulator -> (+ accumulator) . (* i)) (i ^ (2 :: Int))

{-# INLINE generateMultiples #-}