{-
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 .
-}
{- |
[@AUTHOR@] Dr. Alistair Ward
[@DESCRIPTION@] Defines a /prime-wheel/, for use in prime-number generation; .
-}
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; .
* 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/ (), 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 #-}