{-
	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 #-}