{-
	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 {
	PrimeWheel i -> [i]
getPrimeComponents	:: [i],	-- ^ Accessor: the ordered sequence of initial primes, from which the /wheel/ was composed.
	PrimeWheel i -> [i]
getSpokeGaps		:: [i]	-- ^ Accessor: the sequence of spoke-gaps, the sum of which equals its /circumference/.
} deriving Int -> PrimeWheel i -> ShowS
[PrimeWheel i] -> ShowS
PrimeWheel i -> String
(Int -> PrimeWheel i -> ShowS)
-> (PrimeWheel i -> String)
-> ([PrimeWheel i] -> ShowS)
-> Show (PrimeWheel i)
forall i. Show i => Int -> PrimeWheel i -> ShowS
forall i. Show i => [PrimeWheel i] -> ShowS
forall i. Show i => PrimeWheel i -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [PrimeWheel i] -> ShowS
$cshowList :: forall i. Show i => [PrimeWheel i] -> ShowS
show :: PrimeWheel i -> String
$cshow :: forall i. Show i => PrimeWheel i -> String
showsPrec :: Int -> PrimeWheel i -> ShowS
$cshowsPrec :: forall i. Show i => Int -> PrimeWheel i -> ShowS
Show

-- | The /circumference/ of the specified 'PrimeWheel'.
getCircumference :: Integral i => PrimeWheel i -> i
getCircumference :: PrimeWheel i -> i
getCircumference	= [i] -> i
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ([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]
getPrimeComponents

-- | The number of spokes in the specified 'PrimeWheel'.
getSpokeCount :: Integral i => PrimeWheel i -> i
getSpokeCount :: PrimeWheel i -> i
getSpokeCount	= (i -> i -> i) -> i -> [i] -> i
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (i -> i -> i
forall a. Num a => a -> a -> a
(*) (i -> i -> i) -> (i -> i) -> i -> i -> i
forall b c a. (b -> c) -> (a -> b) -> a -> c
. i -> i
forall a. Enum a => a -> a
pred) i
1 ([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]
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 :: Int -> ([Int], [Int])
findCoprimes Int
0	= ([], [])
findCoprimes Int
required
	| Int
required Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0	= String -> ([Int], [Int])
forall a. HasCallStack => String -> a
error (String -> ([Int], [Int])) -> String -> ([Int], [Int])
forall a b. (a -> b) -> a -> b
$ String
"Factory.Data.PrimeWheel.findCoprimes: invalid number of coprimes; " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Int -> String
forall a. Show a => a -> String
show Int
required
	| Bool
otherwise	= Int -> [Int] -> ([Int], [Int])
forall a. Int -> [a] -> ([a], [a])
splitAt Int
required ([Int] -> ([Int], [Int])) -> [Int] -> ([Int], [Int])
forall a b. (a -> b) -> a -> b
$ Int
2 Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: Int -> Int -> Repository -> [Int]
sieve Int
3 Int
0 Repository
forall a. IntMap a
Data.IntMap.empty
	where
		sieve :: Int -> NPrimes -> Repository -> [Int]
		sieve :: Int -> Int -> Repository -> [Int]
sieve Int
candidate Int
found Repository
repository	= case Int -> Repository -> Maybe [Int]
forall a. Int -> IntMap a -> Maybe a
Data.IntMap.lookup Int
candidate Repository
repository of
			Just [Int]
primeMultiples	-> Int -> Repository -> [Int]
sieve' Int
found (Repository -> [Int])
-> (Repository -> Repository) -> Repository -> [Int]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Int] -> Repository -> Repository
insertUniq [Int]
primeMultiples (Repository -> [Int]) -> Repository -> [Int]
forall a b. (a -> b) -> a -> b
$ Int -> Repository -> Repository
forall a. Int -> IntMap a -> IntMap a
Data.IntMap.delete Int
candidate Repository
repository	-- Re-insert subsequent multiples.
			Maybe [Int]
Nothing {-prime-}	-> let
				found' :: Int
found'		= Int -> Int
forall a. Enum a => a -> a
succ Int
found
				(Int
key : [Int]
values)	= (Int -> Int) -> Int -> [Int]
forall a. (a -> a) -> a -> [a]
iterate (Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
gap Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
candidate) (Int -> [Int]) -> Int -> [Int]
forall a b. (a -> b) -> a -> b
$ Int
candidate Int -> Int -> Int
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)	-- Generate a sequence of prime-multiples, starting from its square.
			 in Int
candidate Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: Int -> Repository -> [Int]
sieve' Int
found' (
				if Int
found' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
required
					then Repository
repository
					else Int -> [Int] -> Repository -> Repository
forall a. Int -> a -> IntMap a -> IntMap a
Data.IntMap.insert Int
key [Int]
values Repository
repository
			 )
			where
				gap :: Int
				gap :: Int
gap	= Int
2	-- For efficiency, only sieve odd integers.

				sieve' :: NPrimes -> Repository -> [Int]
				sieve' :: Int -> Repository -> [Int]
sieve'	= Int -> Int -> Repository -> [Int]
sieve (Int -> Int -> Repository -> [Int])
-> Int -> Int -> Repository -> [Int]
forall a b. (a -> b) -> a -> b
$ Int
candidate Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
gap	-- Tail-recurse.

				insertUniq :: PrimeMultiples Int -> Repository -> Repository
				insertUniq :: [Int] -> Repository -> Repository
insertUniq [Int]
l Repository
m	= [Int] -> Repository
insert ([Int] -> Repository) -> [Int] -> Repository
forall a b. (a -> b) -> a -> b
$ (Int -> Bool) -> [Int] -> [Int]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (Int -> Repository -> Bool
forall a. Int -> IntMap a -> Bool
`Data.IntMap.member` Repository
m) [Int]
l	where
					insert :: PrimeMultiples Int -> Repository
					insert :: [Int] -> Repository
insert []		= String -> Repository
forall a. HasCallStack => String -> a
error String
"Factory.Data.PrimeWheel.findCoprimes.sieve.insertUniq.insert:\tnull list"
					insert (Int
key : [Int]
values)	= Int -> [Int] -> Repository -> Repository
forall a. Int -> a -> IntMap a -> IntMap a
Data.IntMap.insert Int
key [Int]
values Repository
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 :: i -> Int
estimateOptimalSize i
maxPrime	= Int -> Int
forall a. Enum a => a -> a
succ (Int -> Int) -> (([Int], [Int]) -> Int) -> ([Int], [Int]) -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Integer] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Integer] -> Int)
-> (([Int], [Int]) -> [Integer]) -> ([Int], [Int]) -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
optimalCircumference) ([Integer] -> [Integer])
-> (([Int], [Int]) -> [Integer]) -> ([Int], [Int]) -> [Integer]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer -> Integer -> Integer) -> [Integer] -> [Integer]
forall a. (a -> a -> a) -> [a] -> [a]
scanl1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
(*) {-circumference-} ([Integer] -> [Integer])
-> (([Int], [Int]) -> [Integer]) -> ([Int], [Int]) -> [Integer]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> Integer) -> [Int] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral {-prevent overflow-} ([Int] -> [Integer])
-> (([Int], [Int]) -> [Int]) -> ([Int], [Int]) -> [Integer]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([Int], [Int]) -> [Int]
forall a b. (a, b) -> a
fst {-primes-} (([Int], [Int]) -> Int) -> ([Int], [Int]) -> Int
forall a b. (a -> b) -> a -> b
$ Int -> ([Int], [Int])
findCoprimes Int
10 {-arbitrary maximum bound-}	where
	optimalCircumference :: Integer
	optimalCircumference :: Integer
optimalCircumference	= Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round (Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ i -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral i
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 :: Int -> PrimeWheel i
mkPrimeWheel Int
0	= [i] -> [i] -> PrimeWheel i
forall i. [i] -> [i] -> PrimeWheel i
MkPrimeWheel [] [i
1]
mkPrimeWheel Int
nPrimes
	| Int
nPrimes Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0	= String -> PrimeWheel i
forall a. HasCallStack => String -> a
error (String -> PrimeWheel i) -> String -> PrimeWheel i
forall a b. (a -> b) -> a -> b
$ String
"Factory.Data.PrimeWheel.mkPrimeWheel: unable to construct from " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Int -> String
forall a. Show a => a -> String
show Int
nPrimes String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
" primes"
	| Bool
otherwise	= PrimeWheel i
primeWheel
	where
		([i]
primeComponents, [i]
coprimeCandidates)	= ((Int -> i) -> [Int] -> [i]
forall a b. (a -> b) -> [a] -> [b]
map Int -> i
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Int] -> [i]) -> ([Int] -> [i]) -> ([Int], [Int]) -> ([i], [i])
forall (a :: * -> * -> *) b c b' c'.
Arrow a =>
a b c -> a b' c' -> a (b, b') (c, c')
*** (Int -> i) -> [Int] -> [i]
forall a b. (a -> b) -> [a] -> [b]
map Int -> i
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Int] -> [i]) -> ([Int] -> [Int]) -> [Int] -> [i]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. i -> [Int] -> [Int]
forall i a. Integral i => i -> [a] -> [a]
Data.List.genericTake (PrimeWheel i -> i
forall i. Integral i => PrimeWheel i -> i
getSpokeCount PrimeWheel i
primeWheel)) (([Int], [Int]) -> ([i], [i])) -> ([Int], [Int]) -> ([i], [i])
forall a b. (a -> b) -> a -> b
$ Int -> ([Int], [Int])
findCoprimes Int
nPrimes
		primeWheel :: PrimeWheel i
primeWheel				= [i] -> [i] -> PrimeWheel i
forall i. [i] -> [i] -> PrimeWheel i
MkPrimeWheel [i]
primeComponents ([i] -> PrimeWheel i) -> [i] -> PrimeWheel i
forall a b. (a -> b) -> a -> b
$ (i -> i -> i) -> [i] -> [i] -> [i]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [i]
coprimeCandidates ([i] -> [i]) -> [i] -> [i]
forall a b. (a -> b) -> a -> b
$ i
1 i -> [i] -> [i]
forall a. a -> [a] -> [a]
: [i]
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 :: Distance i -> Distance i
rotate (i
candidate, [i]
rollingWheel)	= (i
candidate 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] -> i
forall a. [a] -> a
head ([i] -> i) -> ([i] -> [i]) -> [i] -> Distance i
forall (a :: * -> * -> *) b c c'.
Arrow a =>
a b c -> a b c' -> a b (c, c')
&&& [i] -> [i]
forall a. [a] -> [a]
tail ([i] -> Distance i) -> [i] -> Distance i
forall a b. (a -> b) -> a -> b
$ [i]
rollingWheel

{-# INLINE rotate #-}

-- | Generate an infinite, increasing sequence of candidate primes, from the specified /wheel/.
roll :: Integral i => PrimeWheel i -> [Distance i]
roll :: PrimeWheel i -> [Distance i]
roll PrimeWheel i
primeWheel	= [Distance i] -> [Distance i]
forall a. [a] -> [a]
tail ([Distance i] -> [Distance i]) -> [Distance i] -> [Distance i]
forall a b. (a -> b) -> a -> b
$ (Distance i -> Distance i) -> Distance i -> [Distance i]
forall a. (a -> a) -> a -> [a]
iterate Distance i -> Distance i
forall i. Integral i => Distance i -> Distance i
rotate (i
1, [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]
getSpokeGaps PrimeWheel i
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 -> [i] -> [i]
generateMultiples i
i	= (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
i)) (i
i i -> Int -> i
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int))

{-# INLINE generateMultiples #-}