diff --git a/changelog b/changelog --- a/changelog +++ b/changelog @@ -35,4 +35,11 @@ 0.1.0.3 * Qualified 'Factory.Math.Implementations.Primes.trialDivision' with /NOINLINE/ pragma, to block optimization which conflicts with rewrite-rule for 'Factory.Math.Implementations.Primes.sieveOfEratosthenes' ! * Re-coded 'Factory.Data.PrimeWheel.coprimes' and 'Factory.Math.Implementations.Primes.sieveOfEratosthenes', to use a map of lists, rather than a map of lists of lists. - +0.2.0.0 + * Separately coded the special-case of a 'Factory.Data.PrimeWheel' of size zero, in 'Factory.Math.Implementations.Primes.trialDivision', to achieve better space-complexity. + * Added 'Factory.Data.PrimeWheel.estimateOptimalSize'. + * Split "Factory.Math.Implementations.Primes" into; "Factory.Math.Implementations.Primes.SieveOfEratosthenes", "Factory.Math.Implementations.Primes.TurnersSieve", "Factory.Math.Implementations.Primes.TrialDivision", and added a new module "Factory.Math.Implementations.Primes.SieveOfAtkin". This makes the rewrite-rules less fragile. + * Coded 'Factory.Math.Radix.digitalRoot' more concisely. + * Split "Factory.Math.Power" into an additional module "Factory.Math.PerfectPower". + * Replaced '(+ 1)' and '(- 1)' with the faster calls 'succ' and 'pred'. + * Used 'Paths_factory.version' in 'Main', rather than hard-coding it. diff --git a/factory.cabal b/factory.cabal --- a/factory.cabal +++ b/factory.cabal @@ -1,6 +1,6 @@ --Package-properties Name: factory -Version: 0.1.0.3 +Version: 0.2.0.0 Cabal-Version: >= 1.6 Copyright: (C) 2011 Dr. Alistair Ward License: GPL @@ -9,7 +9,7 @@ Stability: Unstable interface, incomplete features. Synopsis: Rational arithmetic in an irrational world. Build-Type: Simple -Description: A library of number-theory functions, for; factorials, square-roots, Pi, primality-testing, prime-factorisation ... +Description: A library of number-theory functions, for; factorials, square-roots, Pi and primes. Category: Math, Number Theory Tested-With: GHC == 6.10, GHC == 6.12, GHC == 7.0 Homepage: http://functionalley.eu @@ -68,9 +68,14 @@ Factory.Math.Implementations.Pi.Spigot.Spigot Factory.Math.Implementations.Primality Factory.Math.Implementations.PrimeFactorisation - Factory.Math.Implementations.Primes + Factory.Math.Implementations.Primes.Algorithm + Factory.Math.Implementations.Primes.SieveOfAtkin + Factory.Math.Implementations.Primes.SieveOfEratosthenes + Factory.Math.Implementations.Primes.TrialDivision + Factory.Math.Implementations.Primes.TurnersSieve Factory.Math.Implementations.SquareRoot Factory.Math.MultiplicativeOrder + Factory.Math.PerfectPower Factory.Math.Pi Factory.Math.Power Factory.Math.Precision @@ -123,6 +128,7 @@ Factory.Test.QuickCheck.Hyperoperation Factory.Test.QuickCheck.Interval Factory.Test.QuickCheck.MonicPolynomial + Factory.Test.QuickCheck.PerfectPower Factory.Test.QuickCheck.Pi Factory.Test.QuickCheck.Polynomial Factory.Test.QuickCheck.Power diff --git a/src/Factory/Data/Interval.hs b/src/Factory/Data/Interval.hs --- a/src/Factory/Data/Interval.hs +++ b/src/Factory/Data/Interval.hs @@ -120,15 +120,15 @@ | otherwise = b -- | Bisect the /interval/ at the specified /end-point/; which should be between the two existing /end-points/. -splitAt' :: (Num endPoint, Ord endPoint) => endPoint -> Interval endPoint -> (Interval endPoint, Interval endPoint) +splitAt' :: (Enum endPoint, Num endPoint, Ord endPoint) => endPoint -> Interval endPoint -> (Interval endPoint, Interval endPoint) splitAt' i interval@(l, r) | any ($ i) [(< l), (>= r)] = error $ "Factory.Data.Interval.splitAt':\tunsuitable index=" ++ show i ++ " for interval=" ++ show interval ++ "." - | otherwise = ((l, i), (i + 1, r)) + | otherwise = ((l, i), (succ i, r)) -- | The length of 'toList'. {-# INLINE getLength #-} -getLength :: (Num endPoint, Ord endPoint) => Interval endPoint -> endPoint -getLength (l, r) = r + 1 - l +getLength :: (Enum endPoint, Num endPoint) => Interval endPoint -> endPoint +getLength (l, r) = succ r - l -- | Converts 'Interval' to a list by enumerating the values. {-# INLINE toList #-} diff --git a/src/Factory/Data/PrimeFactors.hs b/src/Factory/Data/PrimeFactors.hs --- a/src/Factory/Data/PrimeFactors.hs +++ b/src/Factory/Data/PrimeFactors.hs @@ -32,7 +32,6 @@ -- * Functions insert', -- invert, --- merge, product', reduce, -- reduceSorted, @@ -90,7 +89,7 @@ * Preserves the sort-order. - * CAVEAT: this is tolerably efficient for the odd insertion; to insert a list, use '>*<'. + * CAVEAT: this is tolerably efficient for sporadic insertion; to insert a list, use '>*<'. -} insert' :: (Ord base, Num exponent) => Data.Exponential.Exponential base exponent -> Factors base exponent -> Factors base exponent insert' e [] = [e] diff --git a/src/Factory/Data/PrimeWheel.hs b/src/Factory/Data/PrimeWheel.hs --- a/src/Factory/Data/PrimeWheel.hs +++ b/src/Factory/Data/PrimeWheel.hs @@ -24,13 +24,15 @@ -- * Types -- ** Type-synonyms Distance, + NPrimes, PrimeMultiples, -- Repository, -- ** Data-types - PrimeWheel(getPrimeComponents), + PrimeWheel(getPrimeComponents, getSpokeGaps), -- * Functions + estimateOptimalSize, -- findCoprimes, - generatePrimeMultiples, + generateMultiples, roll, rotate, -- ** Constructors @@ -53,7 +55,7 @@ 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 this state, the /circumference/ grows rapidly, with the number of primes: + * 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)] @@ -69,7 +71,7 @@ } deriving Show -- | The /circumference/ of the specified 'PrimeWheel'. -getCircumference :: Num n => PrimeWheel n -> n +getCircumference :: Integral i => PrimeWheel i -> i getCircumference = product . getPrimeComponents -- | The number of spokes in the specified 'PrimeWheel'. @@ -82,43 +84,61 @@ -- | 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. + * 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 :: Int -> ([Int], [Int]) +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 -> Int -> Repository -> [Int] + 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 -> let + 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) + 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' :: Int -> Repository -> [Int] + 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 + 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) + {- | * Constructs a /wheel/ from the specified number of low primes. - * The optimum number of low primes from which to build the /wheel/, grows with the number of primes required; + * 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/; @@ -135,7 +155,7 @@ 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 => Int -> PrimeWheel i +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" @@ -168,11 +188,11 @@ > 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 ..] -} -generatePrimeMultiples :: Integral i - => i -- ^ The /prime/. +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] -generatePrimeMultiples prime = scanl (\accumulator -> (+ accumulator) . (* prime)) (prime ^ (2 :: Int)) +generateMultiples i = scanl (\accumulator -> (+ accumulator) . (* i)) (i ^ (2 :: Int)) -{-# INLINE generatePrimeMultiples #-} +{-# INLINE generateMultiples #-} diff --git a/src/Factory/Math/ArithmeticGeometricMean.hs b/src/Factory/Math/ArithmeticGeometricMean.hs --- a/src/Factory/Math/ArithmeticGeometricMean.hs +++ b/src/Factory/Math/ArithmeticGeometricMean.hs @@ -73,7 +73,7 @@ | spread agm == 0 = repeat agm | otherwise = let simplify :: Data.Ratio.Rational -> Data.Ratio.Rational - simplify = Math.Precision.simplify (decimalDigits - 1 {-ignore single integral digit-}) --This makes a gigantic difference to performance. + simplify = Math.Precision.simplify (pred decimalDigits {-ignore single integral digit-}) --This makes a gigantic difference to performance. findArithmeticMean :: AGM -> ArithmeticMean findArithmeticMean = (/ 2) . uncurry (+) diff --git a/src/Factory/Math/Hyperoperation.hs b/src/Factory/Math/Hyperoperation.hs --- a/src/Factory/Math/Hyperoperation.hs +++ b/src/Factory/Math/Hyperoperation.hs @@ -69,8 +69,8 @@ -} powerTower :: (Integral base, Integral hyperExponent) => base -> hyperExponent -> base powerTower 0 hyperExponent - | odd hyperExponent = 0 - | otherwise = 1 + | even hyperExponent = 1 + | otherwise = 0 powerTower _ (-1) = 0 --The only negative hyper-exponent for which there's a consistent result. powerTower base hyperExponent | base < 0 && hyperExponent > 1 = error $ "Factory.Math.Hyperoperation.powerTower:\tundefined for negative base; " ++ show base diff --git a/src/Factory/Math/Implementations/Factorial.hs b/src/Factory/Math/Implementations/Factorial.hs --- a/src/Factory/Math/Implementations/Factorial.hs +++ b/src/Factory/Math/Implementations/Factorial.hs @@ -65,7 +65,7 @@ factorial algorithm n | n < 2 = 1 | otherwise = case algorithm of - Bisection -> risingFactorial 2 $ n - 1 + Bisection -> risingFactorial 2 $ pred n PrimeFactorisation -> Data.PrimeFactors.product' (recip 5) {-empirical-} 10 {-empirical-} $ primeFactors n {- | @@ -104,7 +104,7 @@ -> i -- ^ The result. risingFactorial _ 0 = 1 risingFactorial 0 _ = 0 -risingFactorial x n = Data.Interval.product' (recip 2) 64 $ Data.Interval.normalise (x, (x + n) - 1) +risingFactorial x n = Data.Interval.product' (recip 2) 64 $ Data.Interval.normalise (x, pred $ x + n) -- | Returns the /falling factorial/; fallingFactorial :: Integral i @@ -113,7 +113,7 @@ -> i -- ^ The result. fallingFactorial _ 0 = 1 fallingFactorial 0 _ = 0 -fallingFactorial x n = Data.Interval.product' (recip 2) 64 $ Data.Interval.normalise (x, (x - n) + 1) +fallingFactorial x n = Data.Interval.product' (recip 2) 64 $ Data.Interval.normalise (x, succ $ x - n) {- | * Returns the ratio of two factorials. diff --git a/src/Factory/Math/Implementations/Pi/Borwein/Borwein1993.hs b/src/Factory/Math/Implementations/Pi/Borwein/Borwein1993.hs --- a/src/Factory/Math/Implementations/Pi/Borwein/Borwein1993.hs +++ b/src/Factory/Math/Implementations/Pi/Borwein/Borwein1993.hs @@ -42,7 +42,7 @@ series = Math.Implementations.Pi.Borwein.Series.MkSeries { Math.Implementations.Pi.Borwein.Series.terms = \squareRootAlgorithm factorialAlgorithm decimalDigits -> let simplify, squareRoot :: Data.Ratio.Rational -> Data.Ratio.Rational - simplify = Math.Precision.simplify (decimalDigits - 1 {-ignore single integral digit-}) --This makes a gigantic difference to performance. + simplify = Math.Precision.simplify $ pred decimalDigits {-ignore single integral digit-} --This makes a gigantic difference to performance. squareRoot = simplify . Math.SquareRoot.squareRoot squareRootAlgorithm decimalDigits sqrt5, a, b, c3 :: Data.Ratio.Rational @@ -64,7 +64,7 @@ ) -} \n power -> ( - Math.Implementations.Factorial.risingFactorial (3 * n + 1) (3 * n) % Math.Power.cube (Math.Factorial.factorial factorialAlgorithm n) + Math.Implementations.Factorial.risingFactorial (succ $ 3 * n) (3 * n) % Math.Power.cube (Math.Factorial.factorial factorialAlgorithm n) ) * ( (a + b * fromIntegral n) / power ) diff --git a/src/Factory/Math/Implementations/Pi/Ramanujan/Chudnovsky.hs b/src/Factory/Math/Implementations/Pi/Ramanujan/Chudnovsky.hs --- a/src/Factory/Math/Implementations/Pi/Ramanujan/Chudnovsky.hs +++ b/src/Factory/Math/Implementations/Pi/Ramanujan/Chudnovsky.hs @@ -52,7 +52,7 @@ ) -} \n power -> ( - Math.Implementations.Factorial.risingFactorial (3 * n + 1) (3 * n) % Math.Power.cube (Math.Factorial.factorial factorialAlgorithm n) + Math.Implementations.Factorial.risingFactorial (succ $ 3 * n) (3 * n) % Math.Power.cube (Math.Factorial.factorial factorialAlgorithm n) ) * ( (13591409 + 545140134 * n) % power ) -- CAVEAT: the order in which these terms are evaluated radically affects performance. diff --git a/src/Factory/Math/Implementations/Pi/Ramanujan/Classic.hs b/src/Factory/Math/Implementations/Pi/Ramanujan/Classic.hs --- a/src/Factory/Math/Implementations/Pi/Ramanujan/Classic.hs +++ b/src/Factory/Math/Implementations/Pi/Ramanujan/Classic.hs @@ -49,7 +49,7 @@ ) $ Math.Implementations.Factorial.primeFactors (4 * n) >/< Math.Implementations.Factorial.primeFactors n >^ 4 -} \n power -> ( - Math.Implementations.Factorial.risingFactorial (n + 1) (3 * n) % Math.Power.cube (Math.Factorial.factorial factorialAlgorithm n) + Math.Implementations.Factorial.risingFactorial (succ n) (3 * n) % Math.Power.cube (Math.Factorial.factorial factorialAlgorithm n) ) * ( (1103 + 26390 * n) % power ) -- CAVEAT: the order in which these terms are evaluated radically affects performance. diff --git a/src/Factory/Math/Implementations/Pi/Spigot/Algorithm.hs b/src/Factory/Math/Implementations/Pi/Spigot/Algorithm.hs --- a/src/Factory/Math/Implementations/Pi/Spigot/Algorithm.hs +++ b/src/Factory/Math/Implementations/Pi/Spigot/Algorithm.hs @@ -46,5 +46,5 @@ openI Gosper = Math.Implementations.Pi.Spigot.Spigot.openI Math.Implementations.Pi.Spigot.Gosper.series openI RabinowitzWagon = Math.Implementations.Pi.Spigot.Spigot.openI Math.Implementations.Pi.Spigot.RabinowitzWagon.series - openR algorithm decimalDigits = Math.Pi.openI algorithm decimalDigits % (10 ^ (decimalDigits - 1)) + openR algorithm decimalDigits = Math.Pi.openI algorithm decimalDigits % (10 ^ pred decimalDigits) diff --git a/src/Factory/Math/Implementations/Pi/Spigot/Gosper.hs b/src/Factory/Math/Implementations/Pi/Spigot/Gosper.hs --- a/src/Factory/Math/Implementations/Pi/Spigot/Gosper.hs +++ b/src/Factory/Math/Implementations/Pi/Spigot/Gosper.hs @@ -31,8 +31,8 @@ -- | Defines a series which converges to /Pi/. series :: Integral i => Math.Implementations.Pi.Spigot.Series.Series i series = Math.Implementations.Pi.Spigot.Series.MkSeries { - Math.Implementations.Pi.Spigot.Series.baseNumerators = map (\i -> i * (2 * i - 1)) [1 ..], - Math.Implementations.Pi.Spigot.Series.baseDenominators = map ((* 3) . (\i -> (i + 1) * (i + 2))) [3, 6 ..], + Math.Implementations.Pi.Spigot.Series.baseNumerators = map (\i -> i * pred (2 * i)) [1 ..], + Math.Implementations.Pi.Spigot.Series.baseDenominators = map ((* 3) . (\i -> succ i * (i + 2))) [3, 6 ..], Math.Implementations.Pi.Spigot.Series.coefficients = [3, 8 ..], --5n - 2 Math.Implementations.Pi.Spigot.Series.nTerms = Math.Precision.getTermsRequired $ 1 / 13 {-empirical convergence-rate-} } diff --git a/src/Factory/Math/Implementations/Pi/Spigot/Spigot.hs b/src/Factory/Math/Implementations/Pi/Spigot/Spigot.hs --- a/src/Factory/Math/Implementations/Pi/Spigot/Spigot.hs +++ b/src/Factory/Math/Implementations/Pi/Spigot/Spigot.hs @@ -108,9 +108,9 @@ -> [(Base, I)] -- ^ Data-row. -> Pi processColumns series preDigits l - | overflowMargin > 1 = preDigits ++ nextRow [digit] --There's neither overflow, nor risk of impact from subsequent overflow. - | overflowMargin == 1 = nextRow $ preDigits ++ [digit] --There's no overflow, but risk of impact from subsequent overflow. - | otherwise = map ((`mod` decimal) . (+ 1)) preDigits ++ nextRow [0] --Overflow => propagate the excess to previously withheld preDigits. + | overflowMargin > 1 = preDigits ++ nextRow [digit] --There's neither overflow, nor risk of impact from subsequent overflow. + | overflowMargin == 1 = nextRow $ preDigits ++ [digit] --There's no overflow, but risk of impact from subsequent overflow. + | otherwise = map ((`mod` decimal) . succ) preDigits ++ nextRow [0] --Overflow => propagate the excess to previously withheld preDigits. where results :: [QuotRem] results = init $ scanr carryAndDivide (0, undefined) l diff --git a/src/Factory/Math/Implementations/Primality.hs b/src/Factory/Math/Implementations/Primality.hs --- a/src/Factory/Math/Implementations/Primality.hs +++ b/src/Factory/Math/Implementations/Primality.hs @@ -47,6 +47,7 @@ import qualified Factory.Data.Polynomial as Data.Polynomial import qualified Factory.Data.QuotientRing as Data.QuotientRing import qualified Factory.Math.MultiplicativeOrder as Math.MultiplicativeOrder +import qualified Factory.Math.PerfectPower as Math.PerfectPower import qualified Factory.Math.Power as Math.Power import qualified Factory.Math.Primality as Math.Primality import qualified Factory.Math.PrimeFactorisation as Math.PrimeFactorisation @@ -108,7 +109,7 @@ -} isPrimeByAKS :: (Math.PrimeFactorisation.Algorithmic factorisationAlgorithm, Control.DeepSeq.NFData i, Integral i) => factorisationAlgorithm -> i -> Bool isPrimeByAKS factorisationAlgorithm n = and [ - not $ Math.Power.isPerfectPower n, --Step 1. + not $ Math.PerfectPower.isPerfectPower n, --Step 1. Math.Primality.areCoprime n `all` filter (/= n) [2 .. r], --Step 3. #if MIN_VERSION_parallel(3,0,0) and $ Control.Parallel.Strategies.parMap Control.Parallel.Strategies.rdeepseq --Benefits from '+RTS -H100M', which reduces garbage-collections. @@ -191,7 +192,7 @@ ) ( length binaryFactors --The number of times that 'two' can be factored-out from 'predecessor'. ) `any` testBases where - predecessor = primeCandidate - 1 + predecessor = pred primeCandidate binaryFactors = takeWhile ((== 0) . snd) . tail {-drop the original-} $ iterate ((`quotRem` 2) . fst) (predecessor, 0) --Factor-out powers of two. testBases | null fewestPrimeBases = let diff --git a/src/Factory/Math/Implementations/PrimeFactorisation.hs b/src/Factory/Math/Implementations/PrimeFactorisation.hs --- a/src/Factory/Math/Implementations/PrimeFactorisation.hs +++ b/src/Factory/Math/Implementations/PrimeFactorisation.hs @@ -47,6 +47,7 @@ import qualified Factory.Data.Exponential as Data.Exponential import Factory.Data.Exponential((<^)) import qualified Factory.Data.PrimeFactors as Data.PrimeFactors +import qualified Factory.Math.PerfectPower as Math.PerfectPower import qualified Factory.Math.Power as Math.Power import qualified Factory.Math.PrimeFactorisation as Math.PrimeFactorisation import qualified ToolShed.Defaultable as Defaultable @@ -112,19 +113,19 @@ Pair.mirror factoriseByFermatsMethod $ head factors where -- maybeSquareNumber :: Integral i => Maybe i - maybeSquareNumber = Math.Power.maybeSquareNumber i + maybeSquareNumber = Math.PerfectPower.maybeSquareNumber i -- factors :: Integral i => [i] factors = map ( ( - uncurry (+) &&& uncurry (-) --Construct the co-factors as the sum and difference of /larger/ and /smaller/. + uncurry (+) &&& uncurry (-) --Construct the co-factors as the sum and difference of /larger/ and /smaller/. ) . Control.Arrow.second Data.Maybe.fromJust ) . filter ( - Data.Maybe.isJust . snd --Search for a perfect square. + Data.Maybe.isJust . snd --Search for a perfect square. ) . map ( - Control.Arrow.second $ Math.Power.maybeSquareNumber {-hotspot-} . (+ negate i) --Associate the corresponding value of /smaller/. + Control.Arrow.second $ Math.PerfectPower.maybeSquareNumber {-hotspot-} . (+ negate i) --Associate the corresponding value of /smaller/. ) . takeWhile ( - (<= (i + 9) `div` 6) . fst --Terminate the search at the maximum value of /larger/. + (<= (i + 9) `div` 6) . fst --Terminate the search at the maximum value of /larger/. ) . Math.Power.squaresFrom {-hotspot-} . ceiling $ sqrt (fromIntegral i :: Double) --Start the search at the minimum value of /larger/. {- | diff --git a/src/Factory/Math/Implementations/Primes.hs b/src/Factory/Math/Implementations/Primes.hs deleted file mode 100644 --- a/src/Factory/Math/Implementations/Primes.hs +++ /dev/null @@ -1,229 +0,0 @@ -{- - 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@] - - * Generates the constant, conceptually infinite, list of /prime-numbers/ by a variety of different algorithms. - - * Based heavily on . - - * . - - * . - - * . --} - -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; . - | TrialDivision Int -- ^ For each candidate, confirm indivisibility, by all /primes/ smaller than its /square-root/, optimised using a 'Data.PrimeWheel.PrimeWheel' (). - | SieveOfEratosthenes Int -- ^ The /Sieve of Eratosthenes/ (), optimised using a 'Data.PrimeWheel.PrimeWheel' (). - | 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; . --} -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; . --} -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; . - - * 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 diff --git a/src/Factory/Math/Implementations/Primes/Algorithm.hs b/src/Factory/Math/Implementations/Primes/Algorithm.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Primes/Algorithm.hs @@ -0,0 +1,63 @@ +{- + 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@] + + * Generates the constant list of /prime-numbers/, by a variety of different algorithms. + + * . + + * . + + * . +-} + +module Factory.Math.Implementations.Primes.Algorithm( +-- * Types +-- ** Data-types + Algorithm(..) +) where + +import qualified Data.Numbers.Primes +import qualified Factory.Data.PrimeWheel as Data.PrimeWheel +import qualified Factory.Math.Implementations.Primes.SieveOfAtkin as Math.Implementations.Primes.SieveOfAtkin +import qualified Factory.Math.Implementations.Primes.SieveOfEratosthenes as Math.Implementations.Primes.SieveOfEratosthenes +import qualified Factory.Math.Implementations.Primes.TrialDivision as Math.Implementations.Primes.TrialDivision +import qualified Factory.Math.Implementations.Primes.TurnersSieve as Math.Implementations.Primes.TurnersSieve +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 + = SieveOfAtkin Integer -- ^ The /Sieve of Atkin/, optimised using a 'Data.PrimeWheel.PrimeWheel' of optimal size, for primes up to the specified maximum bound; . + | SieveOfEratosthenes Data.PrimeWheel.NPrimes -- ^ The /Sieve of Eratosthenes/ (), optimised using a 'Data.PrimeWheel.PrimeWheel'. + | TrialDivision Data.PrimeWheel.NPrimes -- ^ For each candidate, confirm indivisibility, by all /primes/ smaller than its /square-root/, optimised using a 'Data.PrimeWheel.PrimeWheel'. + | TurnersSieve -- ^ For each /prime/, the infinite list of candidates greater than its /square/, is filtered for indivisibility; . + | WheelSieve Int -- ^ 'Data.Numbers.Primes.wheelSieve'. + deriving (Eq, Read, Show) + +instance Defaultable.Defaultable Algorithm where + defaultValue = SieveOfEratosthenes 7 --Resulting in a wheel of circumference 510510. + +instance Math.Primes.Algorithmic Algorithm where + primes (SieveOfAtkin maxPrime) = Math.Implementations.Primes.SieveOfAtkin.sieveOfAtkin (Data.PrimeWheel.estimateOptimalSize maxPrime) $ fromIntegral maxPrime + primes (SieveOfEratosthenes wheelSize) = Math.Implementations.Primes.SieveOfEratosthenes.sieveOfEratosthenes wheelSize + primes (TrialDivision wheelSize) = Math.Implementations.Primes.TrialDivision.trialDivision wheelSize + primes TurnersSieve = Math.Implementations.Primes.TurnersSieve.turnersSieve + primes (WheelSieve wheelSize) = Data.Numbers.Primes.wheelSieve wheelSize --Has better space-complexity than 'SieveOfEratosthenes'. diff --git a/src/Factory/Math/Implementations/Primes/SieveOfAtkin.hs b/src/Factory/Math/Implementations/Primes/SieveOfAtkin.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Primes/SieveOfAtkin.hs @@ -0,0 +1,251 @@ +{-# LANGUAGE CPP #-} +{- + 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@] + + * Generates the constant /bounded/ list of /prime-numbers/, using the /Sieve of Atkin/; . + + * . + + * The implementation; + has been optimised using a /wheel/ of static, but parameterised, size; + has been parallelized; + is polymorphic, but with a specialisation for type 'Int'. + + [@CAVEAT@] The 'Int'-specialisation is implemented by a /rewrite-rule/, which is /very/ fragile. +-} + +module Factory.Math.Implementations.Primes.SieveOfAtkin( +-- * Types +-- ** Data-types +-- PolynomialType, +-- * Constants +-- atkinsModulus, +-- inherentPrimes, +-- nInherentPrimes, +-- squares, +-- * Functions +-- polynomialTypeLookupPeriod, +-- polynomialTypeLookup, +-- findPolynomialSolutions, +-- filterOddRepetitions, +-- generateMultiplesOfSquareTo, +-- getPrefactoredPrimes, + sieveOfAtkin, +-- sieveOfAtkinInt +) where + +import qualified Control.DeepSeq +import qualified Data.Array +import qualified Data.Array.IArray +import Data.Array.IArray((!)) +--import qualified Data.Array.Unboxed +import qualified Data.IntSet +import qualified Data.List +import qualified Data.Set +import qualified Factory.Data.PrimeWheel as Data.PrimeWheel +import qualified Factory.Math.Power as Math.Power +import qualified ToolShed.ListPlus as ListPlus + +#if MIN_VERSION_parallel(3,0,0) +import qualified Control.Parallel.Strategies +#endif + +-- | Defines the types of /quadratic/, available to test the potential primality of a candidate integer. +data PolynomialType + = ModFour -- ^ Suitable for primality-testing numbers meeting @(n `mod` 4 == 1)@. + | ModSix -- ^ Suitable for primality-testing numbers meeting @(n `mod` 6 == 1)@. + | ModTwelve -- ^ Suitable for primality-testing numbers meeting @(n `mod` 12 == 11)@. + | None -- ^ There's no polynomial which can assess primality, because the candidate is composite. + deriving Eq + +-- | The constant modulus used to select the appropriate quadratic for a prime candidate. +atkinsModulus :: Integral i => i +atkinsModulus = foldr1 lcm [4, 6, 12] --Sure, this is always '12', but this is the reason why. + +-- | The constant list of primes factored-out by the unoptimised algorithm. +inherentPrimes :: Integral i => [i] +inherentPrimes = [2, 3] + +-- | The constant number of primes factored-out by the unoptimised algorithm. +nInherentPrimes :: Int +nInherentPrimes = length (inherentPrimes :: [Int]) + +-- | Typically the set of primes which have been built into the specified /wheel/, but never fewer than 'inherentPrimes'. +getPrefactoredPrimes :: Integral i => Data.PrimeWheel.PrimeWheel i -> [i] +getPrefactoredPrimes = max inherentPrimes . Data.PrimeWheel.getPrimeComponents + +-- | The period over which the data returned by 'polynomialTypeLookup' repeats. +polynomialTypeLookupPeriod :: Integral i => Data.PrimeWheel.PrimeWheel i -> i +polynomialTypeLookupPeriod = lcm atkinsModulus . Data.PrimeWheel.getCircumference + +{- | + * Defines which, if any, of the three /quadratics/ is appropriate for the primality-test for each candidate. + + * Since this algorithm uses /modular arithmetic/, the /range/ of results repeat after a short /domain/ related to the /modulus/. + Thus one need calculate at most one period of this cycle, but fewer if the maximum prime required falls within the first cycle of results. + + * Because the results are /bounded/, they're returned in a zero-indexed /array/, to provide efficient random access; + the first few elements should never be required, but it makes query clearer. + + * . +-} +polynomialTypeLookup :: (Data.Array.IArray.Ix i, Integral i) + => Data.PrimeWheel.PrimeWheel i + -> i -- ^ The maximum prime required. +-- -> Data.Array.Unboxed.Array i PolynomialType --Changes neither execution-time nor space ?! + -> Data.Array.Array i PolynomialType +polynomialTypeLookup primeWheel maxPrime = Data.Array.IArray.listArray (0, pred (polynomialTypeLookupPeriod primeWheel) `min` maxPrime) $ map select [0 ..] where +-- select :: Integral i => i -> PolynomialType + select n + | any ( + (== 0) . (n `mod`) --Though this is merely /Trial Division/, it's only performed over a short bounded interval of numerators. + ) primeComponents = None + | r `elem` [1, 5] = ModFour --We actually require @(n `mod` 4 == 1)@, but this is the equivalent modulo 12, with @(r == 9)@ removed because they're all divisible by /3/. + | r == 7 = ModSix --We actually require @(n `mod` 6 == 1)@, but this is the equivalent modulo 12, where @(r == 1)@ has been accounted for above. + | r == 11 = ModTwelve --We require @(n `mod` 12 == 11)@. + | otherwise = None + where + r = n `mod` atkinsModulus + primeComponents = drop nInherentPrimes $ Data.PrimeWheel.getPrimeComponents primeWheel + +-- | The constant, infinite list of the /squares/, of integers increasing from /1/. +squares :: Integral i => [i] +squares = map snd $ Math.Power.squaresFrom 1 + +{- | + * Returns the /ordered/ list of those values with an /odd/ number of occurrences in the specified /unordered/ list. + + * CAVEAT: this is expensive in both execution-time and space. + The typical imperative-style implementation accumulates polynomial-solutions in a /mutable array/ indexed by the candidate integer. + This doesn't translate seamlessly to the /pure functional/ domain where /arrays/ naturally immutable, + so we /sort/ a /list/ of polynomial-solutions, then measure the length of the solution-spans, corresponding to viable candidates. + Regrettably, 'Data.List.sort' (implemented in /GHC/ by /mergesort/) has a time-complexity /O(n*log n)/ + which is greater than the theoretical /O(n)/ of the whole /Sieve of Atkin/; + /GHC/'s old /qsort/-implementation is even slower :( +-} +filterOddRepetitions :: Ord a => [a] -> [a] +--filterOddRepetitions = map head . filter (foldr (const not) False) . Data.List.group . Data.List.sort --Too slow. +filterOddRepetitions = slave True . Data.List.sort where + slave isOdd (one : remainder@(two : _)) + | one == two = slave (not isOdd) remainder + | isOdd = one : beginSpan + | otherwise = beginSpan + where + beginSpan = slave True remainder + slave True [singleton] = [singleton] + slave _ _ = [] + +{- | + * Returns the ordered list of solutions aggregated from each of three /bivariate quadratics/; @z = f(x, y)@. + + * For a candidate integer to be prime, it is necessary but insufficient, that there are an /odd/ number of solutions of value /candidate/. + + * At most one of these three polynomials is suitable for the validation of any specific candidate /z/, depending on 'lookupPolynomialType'. + so the three sets of solutions are mutually exclusive. + One coordinate @(x, y)@, can have solutions in more than one of the three polynomials. + + * This algorithm exhaustively traverses the domain @(x, y)@, for resulting /z/ of the required modulus. + Whilst it tightly constrains the bounds of the search-space, it searches the domain methodically rather than intelligently. +-} +findPolynomialSolutions :: (Control.DeepSeq.NFData i, Data.Array.IArray.Ix i, Integral i) + => Data.PrimeWheel.PrimeWheel i + -> i -- ^ The maximum prime-number required. + -> [i] +findPolynomialSolutions primeWheel maxPrime = foldr1 ListPlus.merge --The lists were previously sorted, as a side-effect, by 'filterOddRepetitions'. +#if MIN_VERSION_parallel(3,0,0) + $ Control.Parallel.Strategies.withStrategy (Control.Parallel.Strategies.parList Control.Parallel.Strategies.rdeepseq) +#endif + [ + {-# SCC "4x^2+y^2" #-} filterOddRepetitions [ + z | + x' <- takeWhile (<= pred maxPrime) $ map (* 4) squares, + z <- takeWhile (<= maxPrime) $ map (+ x') oddSquares, + lookupPolynomialType z == ModFour + ], --Twice the length of the other two lists. + {-# SCC "3x^2+y^2" #-} filterOddRepetitions [ + z | + x' <- takeWhile (<= pred maxPrime) $ map (* 3) squares, + z <- takeWhile (<= maxPrime) . map (+ x') $ if even x' then oddSelection else evenSelection, + lookupPolynomialType z == ModSix + ], + {-# SCC "3x^2-y^2" #-} filterOddRepetitions [ + z | + x2 <- takeWhile (<= maxPrime `div` 2) squares, + z <- dropWhile (> maxPrime) . map (3 * x2 -) . takeWhile (< x2) $ if even x2 then oddSelection else evenSelection, + lookupPolynomialType z == ModTwelve + ] + ] where + (evenSquares, oddSquares) = Data.List.partition even squares + +-- evenSelection, oddSelection :: Integral i => [i] + evenSelection = selection110 evenSquares where + selection110 (x0 : x1 : _ : xs) = x0 : x1 : selection110 xs --Effectively, those for meeting ((== 4) . (`mod` 6)). + selection110 xs = xs + oddSelection = selection101 oddSquares where + selection101 (x0 : _ : x2 : xs) = x0 : x2 : selection101 xs --Effectively, those for meeting ((== 1) . (`mod` 6)). + selection101 xs = xs + +-- lookupPolynomialType :: (Data.Array.IArray.Ix i, Integral i) => i -> PolynomialType + lookupPolynomialType = (polynomialTypeLookup primeWheel maxPrime !) . (`mod` polynomialTypeLookupPeriod primeWheel) + +-- | Generates the /bounded/ list of multiples, of the /square/ of the specified prime, skipping those which aren't required. +generateMultiplesOfSquareTo :: Integral i + => Data.PrimeWheel.PrimeWheel i -- ^ Used to generate the gaps between prime multiples of the square. + -> i -- ^ The /prime/. + -> i -- ^ The maximum bound. + -> [i] +generateMultiplesOfSquareTo primeWheel prime max' = takeWhile (<= max') . scanl (\accumulator -> (+ accumulator) . (* prime2)) prime2 . cycle $ Data.PrimeWheel.getSpokeGaps primeWheel where + prime2 = Math.Power.square prime + +{- | + * Generates the constant /bounded/ list of /prime-numbers/. + + * +-} +sieveOfAtkin :: (Control.DeepSeq.NFData i, Data.Array.IArray.Ix i, Integral i) + => Data.PrimeWheel.NPrimes -- ^ Other implementations effectively use a hard-coded value either /2/ or /3/, but /6/ seems better. + -> i -- ^ The maximum prime required. + -> [i] -- ^ The /bounded/ list of primes. +sieveOfAtkin wheelSize maxPrime = (prefactoredPrimes ++) . filterSquareFree Data.Set.empty . dropWhile (<= maximum prefactoredPrimes) $ findPolynomialSolutions primeWheel maxPrime where + primeWheel = Data.PrimeWheel.mkPrimeWheel wheelSize + prefactoredPrimes = getPrefactoredPrimes primeWheel + +-- filterSquareFree :: Integral i => Data.Set.Set i -> [i] -> [i] + filterSquareFree _ [] = [] + filterSquareFree primeMultiples (candidate : candidates) + | Data.Set.member candidate primeMultiples = {-# SCC "delete" #-} filterSquareFree (Data.Set.delete candidate primeMultiples) candidates --Tail-recurse. + | otherwise = {-# SCC "insert" #-} candidate : filterSquareFree (Data.Set.union primeMultiples . Data.Set.fromDistinctAscList $ generateMultiplesOfSquareTo primeWheel candidate maxPrime) candidates + +{-# NOINLINE sieveOfAtkin #-} +{-# RULES "sieveOfAtkin/Int" sieveOfAtkin = sieveOfAtkinInt #-} --CAVEAT: doesn't fire when built with profiling enabled. + +-- | A specialisation of 'sieveOfAtkin', which reduces both the execution-time and the space required. +sieveOfAtkinInt :: Data.PrimeWheel.NPrimes -> Int -> [Int] +sieveOfAtkinInt wheelSize maxPrime = (prefactoredPrimes ++) . filterSquareFree Data.IntSet.empty . dropWhile (<= maximum prefactoredPrimes) $ findPolynomialSolutions primeWheel maxPrime where + primeWheel = Data.PrimeWheel.mkPrimeWheel wheelSize + prefactoredPrimes = getPrefactoredPrimes primeWheel + + filterSquareFree :: Data.IntSet.IntSet -> [Int] -> [Int] + filterSquareFree _ [] = [] + filterSquareFree primeMultiples (candidate : candidates) + | Data.IntSet.member candidate primeMultiples = filterSquareFree (Data.IntSet.delete candidate primeMultiples) candidates + | otherwise = candidate : filterSquareFree (Data.IntSet.union primeMultiples . Data.IntSet.fromDistinctAscList $ generateMultiplesOfSquareTo primeWheel candidate maxPrime) candidates + diff --git a/src/Factory/Math/Implementations/Primes/SieveOfEratosthenes.hs b/src/Factory/Math/Implementations/Primes/SieveOfEratosthenes.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Primes/SieveOfEratosthenes.hs @@ -0,0 +1,157 @@ +{- + 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@] + + * Generates the constant, conceptually infinite, list of /prime-numbers/, using the /Sieve of Eratosthenes/; . + + * Based on . + + * The implementation; + has been optimised using a /wheel/ of static, but parameterised, size; + is polymorphic, but with a specialisation for type 'Int'. + + [@CAVEAT@] The 'Int'-specialisation is implemented by a /rewrite-rule/, which is /very/ fragile. +-} + +module Factory.Math.Implementations.Primes.SieveOfEratosthenes( +-- * Types +-- ** Type-synonyms +-- PrimeMultiplesQueue, +-- PrimeMultiplesMap, +-- Repository, +-- PrimeMultiplesMapInt, +-- RepositoryInt, +-- * Functions +-- head', +-- tail', + sieveOfEratosthenes, +-- sieveOfEratosthenesInt +) where + +import Control.Arrow((&&&), (***)) +import qualified Control.Arrow +import qualified Data.IntMap +import qualified Data.Map +import Data.Sequence((|>)) +import qualified Data.Sequence +import qualified Factory.Data.PrimeWheel as Data.PrimeWheel + +-- | 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 + +-- | 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; . + + * 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 ultimately demanded, and the effort of constructing and storing this list, is consequently 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 this /queue/, then the whole /list/ of prime-multiples is dropped from the /queue/, + 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 just using a "Data.PQueue.Min". + + * CAVEAT: has linear /O(n)/ space-complexity. +-} +sieveOfEratosthenes :: Integral i + => Data.PrimeWheel.NPrimes + -> [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.generateMultiples 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.generateMultiples 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.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 and reduces the space required. + + * CAVEAT: because the algorithm involves /squares/ of primes, + this implementation will overflow when finding primes greater than @2^16@ on a /32-bit/ machine. +-} +sieveOfEratosthenesInt :: Data.PrimeWheel.NPrimes -> [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.generateMultiples 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.generateMultiples 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.SieveOfEratosthenes.sieveOfEratosthenesInt.sieve.insertUniq.insert:\tnull list" + insert (key : values) = Data.IntMap.insert key values m diff --git a/src/Factory/Math/Implementations/Primes/TrialDivision.hs b/src/Factory/Math/Implementations/Primes/TrialDivision.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Primes/TrialDivision.hs @@ -0,0 +1,61 @@ +{- + 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@] Generates the constant, conceptually infinite, list of /prime-numbers/, using /Trial Division/. +-} + +module Factory.Math.Implementations.Primes.TrialDivision( +-- * Functions + trialDivision +-- ** Predicates +-- isIndivisibleBy +) where + +import qualified Control.Arrow +import qualified Data.List +import qualified Factory.Math.Power as Math.Power +import qualified Factory.Math.PrimeFactorisation as Math.PrimeFactorisation +import qualified Factory.Data.PrimeWheel as Data.PrimeWheel + +-- | Uses /Trial Division/, to determine whether the specified candidate is indivisible by all potential denominators from the specified list. +isIndivisibleBy :: Integral i + => i -- ^ The numerator. + -> [i] -- ^ The denominators of which it must not be a multiple. + -> Bool +isIndivisibleBy numerator = all ((/= 0) . (numerator `mod`)) . takeWhile (<= Math.PrimeFactorisation.maxBoundPrimeFactor numerator) + +{-# INLINE isIndivisibleBy #-} + +{- | + * 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; . +-} +trialDivision :: Integral prime => Data.PrimeWheel.NPrimes -> [prime] +trialDivision 0 = [2, 3] ++ filter (`isIndivisibleBy` trialDivision 0 {-recurse-}) [5 ..] --No faster than using 'Data.PrimeWheel.mkPrimeWheel 0', but apparently better space-complexity ?! +trialDivision wheelSize = Data.PrimeWheel.getPrimeComponents primeWheel ++ indivisible where + primeWheel = Data.PrimeWheel.mkPrimeWheel wheelSize + candidates = map fst $ Data.PrimeWheel.roll primeWheel + indivisible = uncurry (++) . Control.Arrow.second ( + filter (`isIndivisibleBy` 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 + diff --git a/src/Factory/Math/Implementations/Primes/TurnersSieve.hs b/src/Factory/Math/Implementations/Primes/TurnersSieve.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Primes/TurnersSieve.hs @@ -0,0 +1,48 @@ +{- + 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@] Generates the constant, conceptally infinite, list of /prime-numbers/, using /Turner's Sieve/; . +-} + +module Factory.Math.Implementations.Primes.TurnersSieve( +-- * Functions + turnersSieve +) where + +import qualified Factory.Math.Power as Math.Power + +{- | + * For each /prime/, the infinite list of candidates greater than its /square/, + is filtered for indivisibility; . + + * CAVEAT: though one can easily add a 'Data.PrimeWheel.PrimeWheel', it proved counterproductive. +-} +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 discovered. + ] + ) candidates + ) + diff --git a/src/Factory/Math/Implementations/SquareRoot.hs b/src/Factory/Math/Implementations/SquareRoot.hs --- a/src/Factory/Math/Implementations/SquareRoot.hs +++ b/src/Factory/Math/Implementations/SquareRoot.hs @@ -184,7 +184,7 @@ | otherwise = Math.Summation.sumR' . take terms . zipWith (*) taylorSeriesCoefficients $ iterate (* relativeError) x where relativeError :: Math.SquareRoot.Result - relativeError = (realToFrac y / Math.Power.square x) - 1 --Pedantically, this is the error in y, which is twice the magnitude of the error in x. + relativeError = pred $ realToFrac y / Math.Power.square x --Pedantically, this is the error in y, which is twice the magnitude of the error in x. -- | Iterates from the estimated value, towards the /square-root/, a sufficient number of times to achieve the required accuracy. squareRootByIteration :: Real operand => Algorithm -> ProblemSpecification operand diff --git a/src/Factory/Math/PerfectPower.hs b/src/Factory/Math/PerfectPower.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/PerfectPower.hs @@ -0,0 +1,100 @@ +{- + 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@] Exports functions related to /perfect powers/. +-} + +module Factory.Math.PerfectPower( +-- * Functions + maybeSquareNumber, +-- ** Predicates + isPerfectPower +-- isPerfectPowerInt +) where + +import qualified Data.IntSet +import qualified Data.Set +import qualified Factory.Math.Power as Math.Power + +{- | + * Returns @(Just . sqrt)@ if the specified integer is a /square number/ (AKA /perfect square/). + + * . + + * . + + * @(Math.Power.square . sqrt)@ is expensive, so the modulus of the operand is tested first, in an attempt to prove it isn't a /perfect square/. + The set of tests, and the valid moduli within each test, are ordered to maximize the rate of failure-detection. +-} +maybeSquareNumber :: Integral i => i -> Maybe i +maybeSquareNumber i +-- | i < 0 = Nothing --This function is performance-sensitive, but this test is neither strictly nor frequently required. + | all (\(modulus, valid) -> mod i modulus `elem` valid) [ +-- --Distribution of moduli amongst perfect squares Cumulative failure-detection. + (16, [0,1,4,9]), --All moduli are equally likely. 75% + (9, [0,1,4,7]), --Zero occurs 33%, the others only 22%. 88% + (17, [1,2,4,8,9,13,15,16,0]), --Zero only occurs 5.8%, the others 11.8%. 94% +-- These additional tests, aren't always cost-effective. + (13, [1,3,4,9,10,12,0]), --Zero only occurs 7.7%, the others 15.4%. 97% + (7, [1,2,4,0]), --Zero only occurs 14.3%, the others 28.6%. 98% + (5, [1,4,0]) --Zero only occurs 20%, the others 40%. 99% + +-- ] && fromIntegral iSqrt == sqrt' = Just iSqrt --CAVEAT: erroneously True for 187598574531033120 (187598574531033121 is square). + ] && Math.Power.square iSqrt == i = Just iSqrt + | otherwise = Nothing + where + sqrt' :: Double + sqrt' = sqrt $ fromIntegral i + + iSqrt = round sqrt' + +{- | + * An integer @(> 1)@ which can be expressed as an integral power @(> 1)@ of a smaller /natural/ number. + + * CAVEAT: /zero/ and /one/ are normally excluded from this set. + + * . + + * . + + * A generalisation of the concept of /perfect squares/, in which only the exponent '2' is significant. +-} +isPerfectPower :: Integral i => i -> Bool +isPerfectPower i + | i < Math.Power.square 2 = False + | otherwise = i `Data.Set.member` foldr ( + \n set -> if n `Data.Set.member` set + then set +-- else Data.Set.union set . Data.Set.fromDistinctAscList . takeWhile (<= i) . iterate (* n) $ Math.Power.square n + else foldr Data.Set.insert set . takeWhile (<= i) . iterate (* n) $ Math.Power.square n --Faster. + ) Data.Set.empty [2 .. round $ sqrt (fromIntegral i :: Double)] + +{-# NOINLINE isPerfectPower #-} +{-# RULES "isPerfectPower/Int" isPerfectPower = isPerfectPowerInt #-} + +-- | A specialisation of 'isPerfectPower'. +isPerfectPowerInt :: Int -> Bool +isPerfectPowerInt i + | i < Math.Power.square 2 = False + | otherwise = i `Data.IntSet.member` foldr ( + \n set -> if n `Data.IntSet.member` set + then set + else foldr Data.IntSet.insert set . takeWhile (<= i) . iterate (* n) $ Math.Power.square n + ) Data.IntSet.empty [2 .. round $ sqrt (fromIntegral i :: Double)] + diff --git a/src/Factory/Math/Pi.hs b/src/Factory/Math/Pi.hs --- a/src/Factory/Math/Pi.hs +++ b/src/Factory/Math/Pi.hs @@ -47,13 +47,13 @@ openI _ 1 = 3 openI algorithm decimalDigits | decimalDigits <= 0 = error $ "Factory.Math.Pi.openI:\tinsufficient decimalDigits=" ++ show decimalDigits - | otherwise = round . Math.Precision.promote (openR algorithm decimalDigits) $ decimalDigits - 1 + | otherwise = round . Math.Precision.promote (openR algorithm decimalDigits) $ pred decimalDigits openS :: algorithm -> Math.Precision.DecimalDigits -> String -- ^ Returns the value of /Pi/ as a decimal 'String'. openS _ 1 = "3" openS algorithm decimalDigits | decimalDigits <= 0 = "" - | decimalDigits <= 16 = take (decimalDigits + 1) $ show (pi :: Double) + | decimalDigits <= 16 = take (succ decimalDigits) $ show (pi :: Double) | otherwise = "3." ++ tail (show $ openI algorithm decimalDigits) --Insert a decimal point. -- | Categorises the various algorithms. @@ -83,7 +83,7 @@ ) => Algorithmic (Category agm bbp borwein ramanujan spigot) where openR algorithm decimalDigits | decimalDigits <= 0 = error $ "Factory.Math.Pi.openR:\tinsufficient decimalDigits=" ++ show decimalDigits - | decimalDigits <= 16 = Math.Precision.simplify (decimalDigits - 1) (pi :: Double) + | decimalDigits <= 16 = Math.Precision.simplify (pred decimalDigits) (pi :: Double) | otherwise = ( case algorithm of AGM agm -> openR agm @@ -97,5 +97,5 @@ openI (Spigot spigot) decimalDigits = openI spigot decimalDigits openI algorithm decimalDigits | decimalDigits <= 0 = error $ "Factory.Math.Pi.openI:\tinsufficient decimalDigits=" ++ show decimalDigits - | otherwise = round . Math.Precision.promote (openR algorithm decimalDigits) $ decimalDigits - 1 + | otherwise = round . Math.Precision.promote (openR algorithm decimalDigits) $ pred decimalDigits diff --git a/src/Factory/Math/Power.hs b/src/Factory/Math/Power.hs --- a/src/Factory/Math/Power.hs +++ b/src/Factory/Math/Power.hs @@ -24,23 +24,17 @@ -- * Functions square, squaresFrom, - maybeSquareNumber, cube, cubeRoot, - raiseModulo, --- ** Predicates - isPerfectPower --- isPerfectPowerInt + raiseModulo ) where -import qualified Data.IntSet -import qualified Data.Set - -- | Mainly for convenience. -{-# INLINE square #-} square :: Num n => n -> n square = (^ (2 :: Int)) +{-# INLINE square #-} + -- | Just for convenience. cube :: Num n => n -> n cube = (^ (3 :: Int)) @@ -51,10 +45,10 @@ * The initial value doesn't need to be either positive or integral. -} -squaresFrom :: Num n +squaresFrom :: (Enum n, Num n) => n -- ^ Lower bound. -> [(n, n)] -- ^ @ [(n, n^2)] @. -squaresFrom from = iterate (\(x, y) -> (x + 1, y + 2 * x + 1)) (from, square from) +squaresFrom from = iterate (\(x, y) -> (succ x, succ $ y + 2 * x)) (from, square from) -- | Just for convenience. cubeRoot :: Double -> Double @@ -77,7 +71,7 @@ raiseModulo _ _ 1 = 0 raiseModulo _ 0 modulus = 1 `mod` modulus raiseModulo base power modulus - | base < 0 = (`mod` modulus) . (if odd power then negate else id) $ raiseModulo (negate base) power modulus --Recurse. + | base < 0 = (`mod` modulus) . (if even power then id else negate) $ raiseModulo (negate base) power modulus --Recurse. | power < 0 = error $ "Factory.Math.Power.raiseModulo:\tnegative power; " ++ show power | first `elem` [0, 1] = first | otherwise = slave power @@ -87,70 +81,4 @@ slave 1 = first slave e = (`mod` modulus) . (if r == 0 {-even-} then id else (* base)) . square $ slave q {-recurse-} where (q, r) = e `quotRem` 2 - -{- | - * Returns @(Just . sqrt)@ if the specified integer is a /square number/ (AKA /perfect square/). - - * . - - * . - - * @(square . sqrt)@ is expensive, so the modulus of the operand is tested first, in an attempt to prove it isn't a /perfect square/. - The set of tests, and the valid moduli within each test, are ordered to maximize the rate of failure-detection. --} -maybeSquareNumber :: Integral i => i -> Maybe i -maybeSquareNumber i --- | i < 0 = Nothing --This function is performance-sensitive, but this test is neither strictly nor frequently required. - | all (\(modulus, valid) -> mod i modulus `elem` valid) [ --- --Distribution of moduli amongst perfect squares Cumulative failure-detection. - (16, [0,1,4,9]), --All moduli are equally likely. 75% - (9, [0,1,4,7]), --Zero occurs 33%, the others only 22%. 88% - (17, [1,2,4,8,9,13,15,16,0]), --Zero only occurs 5.8%, the others 11.8%. 94% --- These additional tests, aren't always cost-effective. - (13, [1,3,4,9,10,12,0]), --Zero only occurs 7.7%, the others 15.4%. 97% - (7, [1,2,4,0]), --Zero only occurs 14.3%, the others 28.6%. 98% - (5, [1,4,0]) --Zero only occurs 20%, the others 40%. 99% - --- ] && fromIntegral iSqrt == sqrt' = Just iSqrt --CAVEAT: erroneously True for 187598574531033120 (187598574531033121 is square). - ] && square iSqrt == i = Just iSqrt - | otherwise = Nothing - where - sqrt' :: Double - sqrt' = sqrt $ fromIntegral i - - iSqrt = round sqrt' - -{- | - * An integer @(> 1)@ which can be expressed as an integral power @(> 1)@ of a smaller /natural/ number. - - * CAVEAT: /zero/ and /one/ are normally excluded from this set. - - * . - - * . - - * A generalisation of the concept of /perfect squares/, in which only the exponent '2' is significant. --} -isPerfectPower :: Integral i => i -> Bool -isPerfectPower i - | i < square 2 = False - | otherwise = i `Data.Set.member` foldr ( - \n set -> if n `Data.Set.member` set - then set --- else Data.Set.union set . Data.Set.fromDistinctAscList . takeWhile (<= i) . iterate (* n) $ square n - else foldr Data.Set.insert set . takeWhile (<= i) . iterate (* n) $ square n --Faster. - ) Data.Set.empty [2 .. round $ sqrt (fromIntegral i :: Double)] - -{-# NOINLINE isPerfectPower #-} -{-# RULES "isPerfectPower/Int" isPerfectPower = isPerfectPowerInt #-} - --- | A specialisation of 'isPerfectPower'. -isPerfectPowerInt :: Int -> Bool -isPerfectPowerInt i - | i < square 2 = False - | otherwise = i `Data.IntSet.member` foldr ( - \n set -> if n `Data.IntSet.member` set - then set - else foldr Data.IntSet.insert set . takeWhile (<= i) . iterate (* n) $ square n - ) Data.IntSet.empty [2 .. round $ sqrt (fromIntegral i :: Double)] diff --git a/src/Factory/Math/Precision.hs b/src/Factory/Math/Precision.hs --- a/src/Factory/Math/Precision.hs +++ b/src/Factory/Math/Precision.hs @@ -114,5 +114,5 @@ => DecimalDigits -- ^ The number of places after the decimal point, which are required. -> operand -> Data.Ratio.Rational -simplify decimalDigits operand = Data.Ratio.approxRational operand . recip $ 4 * 10 ^ (decimalDigits + 1) --Tolerate any error less than half the least significant digit required. +simplify decimalDigits operand = Data.Ratio.approxRational operand . recip $ 4 * 10 ^ succ decimalDigits --Tolerate any error less than half the least significant digit required. diff --git a/src/Factory/Math/Primality.hs b/src/Factory/Math/Primality.hs --- a/src/Factory/Math/Primality.hs +++ b/src/Factory/Math/Primality.hs @@ -69,11 +69,11 @@ * TODO: confirm that all values must be tested. -} isFermatWitness :: Integral i => i -> Bool -isFermatWitness i = not . all isFermatPseudoPrime $ filter (areCoprime i) [2 .. i - 1] where - isFermatPseudoPrime base = Math.Power.raiseModulo base (i - 1) i == 1 --CAVEAT: a /Fermat Pseudo-prime/ must also be a /composite/ number. +isFermatWitness i = not . all isFermatPseudoPrime $ filter (areCoprime i) [2 .. pred i] where + isFermatPseudoPrime base = Math.Power.raiseModulo base (pred i) i == 1 --CAVEAT: a /Fermat Pseudo-prime/ must also be a /composite/ number. {- | - * A /Carmichael number/ is an odd /composite/ number which satisfies /Fermat's little theorem/. + * A /Carmichael number/ is an /odd/ /composite/ number which satisfies /Fermat's little theorem/. * . diff --git a/src/Factory/Math/Primes.hs b/src/Factory/Math/Primes.hs --- a/src/Factory/Math/Primes.hs +++ b/src/Factory/Math/Primes.hs @@ -27,16 +27,19 @@ primorial ) where +import qualified Control.DeepSeq +import qualified Data.Array.IArray + -- | Defines the methods expected of a /prime-number/ generator. class Algorithmic algorithm where - primes :: Integral i => algorithm -> [i] -- ^ Returns the constant, conceptually infinite, list of primes. + primes :: (Control.DeepSeq.NFData i, Data.Array.IArray.Ix i, Integral i) => algorithm -> [i] -- ^ Returns the constant, infinite, list of primes. {- | - * Returns the constant, infinite list, defining the /Primorial/. + * Returns the constant list, defining the /Primorial/. * . * . -} -primorial :: (Integral i, Algorithmic algorithm) => algorithm -> [i] +primorial :: (Algorithmic algorithm, Control.DeepSeq.NFData i, Data.Array.IArray.Ix i, Integral i) => algorithm -> [i] primorial = scanl (*) 1 . primes diff --git a/src/Factory/Math/Radix.hs b/src/Factory/Math/Radix.hs --- a/src/Factory/Math/Radix.hs +++ b/src/Factory/Math/Radix.hs @@ -44,7 +44,7 @@ -- | Constant random-access lookup for 'digits'. encodes :: (Data.Array.IArray.Ix index, Integral index) => Data.Array.IArray.Array index Char -encodes = Data.Array.IArray.listArray (0, fromIntegral $ length digits - 1) digits where +encodes = Data.Array.IArray.listArray (0, fromIntegral . pred $ length digits) digits where -- | Constant reverse-lookup for 'digits'. decodes :: Integral i => [(Char, i)] @@ -69,7 +69,7 @@ where fromDecimal 0 = id fromDecimal n - | remainder < 0 = fromDecimal (quotient + 1) . ((remainder - fromIntegral base) :) --This can only occur when base is negative; cf. 'divMod'. + | remainder < 0 = fromDecimal (succ quotient) . ((remainder - fromIntegral base) :) --This can only occur when base is negative; cf. 'divMod'. | otherwise = fromDecimal quotient . (remainder :) where (quotient, remainder) = n `quotRem` fromIntegral base @@ -114,5 +114,5 @@ -- | . digitalRoot :: (Data.Array.IArray.Ix decimal, Integral decimal) => decimal -> decimal -digitalRoot = head . dropWhile (> 9) . iterate (digitSum (10 :: Int)) +digitalRoot = until (<= 9) (digitSum (10 :: Int)) diff --git a/src/Factory/Math/Statistics.hs b/src/Factory/Math/Statistics.hs --- a/src/Factory/Math/Statistics.hs +++ b/src/Factory/Math/Statistics.hs @@ -103,7 +103,7 @@ | otherwise = numerator `par` (denominator `pseq` numerator `div` denominator) where [smaller, bigger] = Data.List.sort [r, n - r] - numerator = Math.Implementations.Factorial.risingFactorial (bigger + 1) (n - bigger) + numerator = Math.Implementations.Factorial.risingFactorial (succ bigger) (n - bigger) denominator = Math.Factorial.factorial factorialAlgorithm smaller -- | The number of /permutations/ of /r/ objects taken from /n/; . diff --git a/src/Factory/Test/Performance/Primes.hs b/src/Factory/Test/Performance/Primes.hs --- a/src/Factory/Test/Performance/Primes.hs +++ b/src/Factory/Test/Performance/Primes.hs @@ -26,9 +26,10 @@ ) where import qualified Control.DeepSeq +import qualified Data.Array.IArray import qualified Factory.Math.Primes as Math.Primes import qualified ToolShed.TimePure as TimePure -- | Measures the CPU-time required by 'Math.Primes.primes', to find the specified prime. -primesPerformance :: (Math.Primes.Algorithmic algorithm, Control.DeepSeq.NFData i, Integral i) => algorithm -> Int -> IO (Double, i) +primesPerformance :: (Math.Primes.Algorithmic algorithm, Control.DeepSeq.NFData i, Data.Array.IArray.Ix i, Integral i) => algorithm -> Int -> IO (Double, i) primesPerformance algorithm = TimePure.getCPUSeconds . (Math.Primes.primes algorithm !!) diff --git a/src/Factory/Test/QuickCheck/ArithmeticGeometricMean.hs b/src/Factory/Test/QuickCheck/ArithmeticGeometricMean.hs --- a/src/Factory/Test/QuickCheck/ArithmeticGeometricMean.hs +++ b/src/Factory/Test/QuickCheck/ArithmeticGeometricMean.hs @@ -55,11 +55,11 @@ ) ( Math.ArithmeticGeometricMean.convergeToAGM squareRootAlgorithm decimalDigits' $ swap agm ) where - decimalDigits' = 1 + (decimalDigits `mod` 64) - index' = 1 + (index `mod` 8) + decimalDigits' = succ $ decimalDigits `mod` 64 + index' = succ $ index `mod` 8 prop_bounds squareRootAlgorithm decimalDigits agm index = all ($ agm) [Math.ArithmeticGeometricMean.isValid, uncurry (/=)] ==> Test.QuickCheck.label "prop_bounds" . all (uncurry (>=)) . tail . take index' $ Math.ArithmeticGeometricMean.convergeToAGM squareRootAlgorithm decimalDigits' agm where decimalDigits' = 33 {-test is sensitive to rounding-errors-} + (decimalDigits `mod` 96) - index' = 1 + (index `mod` 5) + index' = succ $ index `mod` 5 diff --git a/src/Factory/Test/QuickCheck/Factorial.hs b/src/Factory/Test/QuickCheck/Factorial.hs --- a/src/Factory/Test/QuickCheck/Factorial.hs +++ b/src/Factory/Test/QuickCheck/Factorial.hs @@ -52,10 +52,10 @@ prop_equivalence x n = Test.QuickCheck.label "prop_equivalence" $ Math.Implementations.Factorial.risingFactorial x n == sign * Math.Implementations.Factorial.fallingFactorial (negate x) n && Math.Implementations.Factorial.fallingFactorial x n == sign * Math.Implementations.Factorial.risingFactorial (negate x) n where sign :: Integer sign - | odd n = negate 1 - | otherwise = 1 + | even n = 1 + | otherwise = negate 1 - prop_symmetry x n = Test.QuickCheck.label "prop_symmetry" $ Math.Implementations.Factorial.risingFactorial x n == Math.Implementations.Factorial.fallingFactorial (x + n - 1) n + prop_symmetry x n = Test.QuickCheck.label "prop_symmetry" $ Math.Implementations.Factorial.risingFactorial x n == Math.Implementations.Factorial.fallingFactorial (pred $ x + n) n prop_x0 x _ = Test.QuickCheck.label "prop_x0" $ all (== 1) $ map ($ 0) [Math.Implementations.Factorial.risingFactorial x, Math.Implementations.Factorial.fallingFactorial x] @@ -63,10 +63,10 @@ prop_ratio :: Math.Implementations.Factorial.Algorithm -> Integer -> Integer -> Test.QuickCheck.Property prop_ratio algorithm i j = Test.QuickCheck.label "prop_ratio" $ n !/! d == Math.Factorial.factorial algorithm n % Math.Factorial.factorial algorithm d where - n = (i `mod` 100000) - 1 - d = (j `mod` 100000) - 1 + n = pred $ i `mod` 100000 + d = pred $ j `mod` 100000 prop_consistency :: Math.Implementations.Factorial.Algorithm -> Math.Implementations.Factorial.Algorithm -> Integer -> Test.QuickCheck.Property prop_consistency l r i = l /= r ==> Test.QuickCheck.label "prop_consistency" $ Math.Factorial.factorial l n == Math.Factorial.factorial r n where - n = (i `mod` 100000) - 1 + n = pred $ i `mod` 100000 diff --git a/src/Factory/Test/QuickCheck/Hyperoperation.hs b/src/Factory/Test/QuickCheck/Hyperoperation.hs --- a/src/Factory/Test/QuickCheck/Hyperoperation.hs +++ b/src/Factory/Test/QuickCheck/Hyperoperation.hs @@ -40,7 +40,7 @@ prop_rankCoincides :: Rank -> Test.QuickCheck.Property prop_rankCoincides rank = Test.QuickCheck.label "prop_rankCoincides" $ Math.Hyperoperation.hyperoperation rank' 2 2 == 4 where rank' :: Rank - rank' = 1 + (rank `mod` 1000) + rank' = succ $ rank `mod` 1000 prop_baseCoincides :: Rank -> Integer -> Test.QuickCheck.Property prop_baseCoincides rank base = Test.QuickCheck.label "prop_baseCoincides" $ Math.Hyperoperation.hyperoperation rank' base 1 == base where @@ -56,7 +56,7 @@ hyperExponent' = abs hyperExponent prop_succ, prop_addition, prop_multiplication, prop_exponentiation :: Integer -> Integer -> Test.QuickCheck.Property - prop_succ base hyperExponent = Test.QuickCheck.label "prop_succ" $ Math.Hyperoperation.hyperoperation Math.Hyperoperation.succession base hyperExponent' == 1 + fromIntegral hyperExponent' where + prop_succ base hyperExponent = Test.QuickCheck.label "prop_succ" $ Math.Hyperoperation.hyperoperation Math.Hyperoperation.succession base hyperExponent' == succ (fromIntegral hyperExponent') where hyperExponent' :: Math.Hyperoperation.HyperExponent hyperExponent' = abs hyperExponent diff --git a/src/Factory/Test/QuickCheck/Interval.hs b/src/Factory/Test/QuickCheck/Interval.hs --- a/src/Factory/Test/QuickCheck/Interval.hs +++ b/src/Factory/Test/QuickCheck/Interval.hs @@ -35,9 +35,9 @@ prop_product :: Data.Ratio.Ratio Integer -> Integer -> Data.Interval.Interval Integer -> Test.QuickCheck.Property prop_product ratio minLength interval = Test.QuickCheck.label "prop_product" $ Data.Interval.product' ratio' minLength' interval' == product (Data.Interval.toList interval') where interval' = Data.Interval.normalise interval - minLength' = 1 + minLength `mod` 1000 - ratio' = if r > 1 - then recip r - else r + minLength' = succ $ minLength `mod` 1000 + ratio' + | r > 1 = recip r + | otherwise = r where r = abs ratio diff --git a/src/Factory/Test/QuickCheck/MonicPolynomial.hs b/src/Factory/Test/QuickCheck/MonicPolynomial.hs --- a/src/Factory/Test/QuickCheck/MonicPolynomial.hs +++ b/src/Factory/Test/QuickCheck/MonicPolynomial.hs @@ -47,7 +47,7 @@ arbitrary = do polynomial <- Test.QuickCheck.arbitrary - return . Data.MonicPolynomial.mkMonicPolynomial $ ((1, Data.Polynomial.getDegree polynomial + 1) :) `Data.Polynomial.lift` polynomial + return . Data.MonicPolynomial.mkMonicPolynomial $ ((1, succ $ Data.Polynomial.getDegree polynomial) :) `Data.Polynomial.lift` polynomial #if !(MIN_VERSION_QuickCheck(2,1,0)) coarbitrary = undefined --CAVEAT: stops warnings from ghc. #endif @@ -66,7 +66,7 @@ prop_perfectPower :: P -> Int -> Test.QuickCheck.Property prop_perfectPower polynomial power = Test.QuickCheck.label "prop_perfectPower" $ iterate (`Data.QuotientRing.quot'` polynomial) (polynomial =^ power') !! pred power' == polynomial where power' :: Int - power' = 1 + power `mod` 100 + power' = succ $ power `mod` 100 prop_isDivisibleBy :: [P] -> Test.QuickCheck.Property prop_isDivisibleBy monicPolynomials = Test.QuickCheck.label "prop_isDivisibleBy" $ all (Data.QuotientRing.isDivisibleBy (Data.Ring.product' (recip 2) {-TODO-} 10 monicPolynomials)) monicPolynomials diff --git a/src/Factory/Test/QuickCheck/PerfectPower.hs b/src/Factory/Test/QuickCheck/PerfectPower.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/QuickCheck/PerfectPower.hs @@ -0,0 +1,52 @@ +{- + 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 /QuickCheck/-properties for "Math.PerfectPower". +-} + +module Factory.Test.QuickCheck.PerfectPower( +-- * Functions + quickChecks +) where + +import qualified Factory.Math.PerfectPower as Math.PerfectPower +import qualified Factory.Math.Power as Math.Power +import qualified Test.QuickCheck +import Test.QuickCheck((==>)) + +-- | Defines invariant properties. +quickChecks :: IO () +quickChecks = + Test.QuickCheck.quickCheck `mapM_` [prop_maybeSquareNumber, prop_rewriteRule] + >> Test.QuickCheck.quickCheckWith Test.QuickCheck.stdArgs {Test.QuickCheck.maxSuccess = 10000} prop_notSquare + >> Test.QuickCheck.quickCheck prop_isPerfectPower + where + prop_maybeSquareNumber, prop_notSquare, prop_rewriteRule :: Integer -> Test.QuickCheck.Property + prop_maybeSquareNumber i = Test.QuickCheck.label "prop_maybeSquareNumber" $ Math.PerfectPower.maybeSquareNumber (Math.Power.square i) == Just (abs i) + + prop_notSquare i = abs i > 0 ==> Test.QuickCheck.label "prop_notSquare" $ Math.PerfectPower.maybeSquareNumber (succ $ i ^ (10 {-promote rounding-error using big number-} :: Int)) == Nothing + prop_rewriteRule i = Test.QuickCheck.label "prop_rewriteRule" $ Math.PerfectPower.isPerfectPower i' == Math.PerfectPower.isPerfectPower (fromIntegral i' :: Int) where + i' = abs i + + prop_isPerfectPower :: Integer -> Integer -> Test.QuickCheck.Property + prop_isPerfectPower b e = Test.QuickCheck.label "prop_isPerfectPower" . Math.PerfectPower.isPerfectPower $ b' ^ e' where + b' = 2 + (b `mod` 10) + e' = 2 + (e `mod` 8) + + diff --git a/src/Factory/Test/QuickCheck/Pi.hs b/src/Factory/Test/QuickCheck/Pi.hs --- a/src/Factory/Test/QuickCheck/Pi.hs +++ b/src/Factory/Test/QuickCheck/Pi.hs @@ -128,5 +128,5 @@ quickChecks = Test.QuickCheck.quickCheck prop_consistency where prop_consistency :: Testable prop_consistency l r decimalDigits = l /= r ==> Test.QuickCheck.label "prop_consistency" $ Math.Pi.openI l decimalDigits' - Math.Pi.openI r decimalDigits' <= 1 {-rounding error-} where - decimalDigits' = 1 + (decimalDigits `mod` 250) + decimalDigits' = succ $ decimalDigits `mod` 250 diff --git a/src/Factory/Test/QuickCheck/Polynomial.hs b/src/Factory/Test/QuickCheck/Polynomial.hs --- a/src/Factory/Test/QuickCheck/Polynomial.hs +++ b/src/Factory/Test/QuickCheck/Polynomial.hs @@ -88,30 +88,30 @@ prop_power, prop_perfectPower, prop_normalised :: Data.Polynomial.Polynomial Integer Integer -> Int -> Test.QuickCheck.Property prop_power polynomial power = Test.QuickCheck.label "prop_power" $ polynomial =^ power' == iterate (=*= polynomial) polynomial !! pred power' where power' :: Int - power' = 1 + power `mod` 100 + power' = succ $ power `mod` 100 prop_perfectPower polynomial power = polynomial' /= Data.Polynomial.zero ==> Test.QuickCheck.label "prop_perfectPower" $ iterate (`Data.QuotientRing.quot'` polynomial') (polynomial' =^ power') !! pred power' == polynomial' where polynomial' :: Data.Polynomial.Polynomial Data.Ratio.Rational Integer polynomial' = Data.Polynomial.realCoefficientsToFrac polynomial power' :: Int - power' = 1 + power `mod` 100 + power' = succ $ power `mod` 100 prop_normalised polynomial i = Test.QuickCheck.label "prop_normalised" $ all Data.Polynomial.isNormalised [ polynomial =^ power', polynomial `Data.Polynomial.mod'` modulus' ] where power' :: Int - power' = 1 + i `mod` 100 + power' = succ $ i `mod` 100 modulus' :: Integer - modulus' = 1 + fromIntegral i `mod` 100 + modulus' = succ $ fromIntegral i `mod` 100 prop_raiseModuloNormalised :: Data.Polynomial.Polynomial Integer Integer -> Integer -> Integer -> Test.QuickCheck.Property prop_raiseModuloNormalised polynomial power modulus = Test.QuickCheck.label "prop_raiseModuloNormalised" . Data.Polynomial.isNormalised $ Data.Polynomial.raiseModulo polynomial power' modulus' where power', modulus' :: Integer - power' = 1 + power `mod` 100 - modulus' = 1 + modulus `mod` 100 + power' = succ $ power `mod` 100 + modulus' = succ $ modulus `mod` 100 prop_integralDomain, prop_isDivisibleBy :: [Data.Polynomial.Polynomial Integer Integer] -> Test.QuickCheck.Property prop_integralDomain polynomials = Data.Polynomial.zero `notElem` polynomials ==> Test.QuickCheck.label "prop_integralDomain" $ Data.Ring.product' (recip 2) {-TODO-} 10 polynomials /= Data.Polynomial.zero diff --git a/src/Factory/Test/QuickCheck/Power.hs b/src/Factory/Test/QuickCheck/Power.hs --- a/src/Factory/Test/QuickCheck/Power.hs +++ b/src/Factory/Test/QuickCheck/Power.hs @@ -17,7 +17,7 @@ {- | [@AUTHOR@] Dr. Alistair Ward - [@DESCRIPTION@] Defines /QuickCheck/-properties for "Math.Power". + [@DESCRIPTION@] Defines /QuickCheck/-properties "Math.Power". -} module Factory.Test.QuickCheck.Power( @@ -32,25 +32,10 @@ -- | Defines invariant properties. quickChecks :: IO () -quickChecks = - Test.QuickCheck.quickCheck `mapM_` [prop_maybeSquareNumber, prop_rewriteRule] - >> Test.QuickCheck.quickCheckWith Test.QuickCheck.stdArgs {Test.QuickCheck.maxSuccess = 10000} prop_notSquare - >> Test.QuickCheck.quickCheck `mapM` [prop_squaresFrom, prop_isPerfectPower] - >> Test.QuickCheck.quickCheck prop_raiseModulo +quickChecks = Test.QuickCheck.quickCheck prop_squaresFrom >> Test.QuickCheck.quickCheck prop_raiseModulo where - prop_maybeSquareNumber, prop_notSquare, prop_rewriteRule :: Integer -> Test.QuickCheck.Property - prop_maybeSquareNumber i = Test.QuickCheck.label "prop_maybeSquareNumber" $ Math.Power.maybeSquareNumber (Math.Power.square i) == Just (abs i) - - prop_notSquare i = abs i > 0 ==> Test.QuickCheck.label "prop_notSquare" $ Math.Power.maybeSquareNumber (i ^ (10 {-promote rounding-error using big number-} :: Int) + 1) == Nothing - prop_rewriteRule i = Test.QuickCheck.label "prop_rewriteRule" $ Math.Power.isPerfectPower i' == Math.Power.isPerfectPower (fromIntegral i' :: Int) where - i' = abs i - - prop_squaresFrom, prop_isPerfectPower :: Integer -> Integer -> Test.QuickCheck.Property + prop_squaresFrom :: Integer -> Integer -> Test.QuickCheck.Property prop_squaresFrom from l = Test.QuickCheck.label "prop_squaresFrom" . (\(x, y) -> y == Math.Power.square x) . Data.List.genericIndex (Math.Power.squaresFrom from) $ abs l - - prop_isPerfectPower b e = Test.QuickCheck.label "prop_isPerfectPower" . Math.Power.isPerfectPower $ b' ^ e' where - b' = 2 + (b `mod` 10) - e' = 2 + (e `mod` 8) prop_raiseModulo :: Integer -> Integer -> Integer -> Test.QuickCheck.Property prop_raiseModulo b e m = m /= 0 ==> Test.QuickCheck.label "prop_raiseModulo" $ Math.Power.raiseModulo b e' m == (b ^ e') `mod` m where diff --git a/src/Factory/Test/QuickCheck/PrimeFactorisation.hs b/src/Factory/Test/QuickCheck/PrimeFactorisation.hs --- a/src/Factory/Test/QuickCheck/PrimeFactorisation.hs +++ b/src/Factory/Test/QuickCheck/PrimeFactorisation.hs @@ -58,30 +58,30 @@ prop_consistency :: Integer -> Test.QuickCheck.Property prop_consistency i = Test.QuickCheck.label "prop_consistency" $ (Math.PrimeFactorisation.primeFactors Math.Implementations.PrimeFactorisation.TrialDivision i' :: Data.PrimeFactors.Factors Integer Int) == Math.PrimeFactorisation.primeFactors Math.Implementations.PrimeFactorisation.FermatsMethod i' where i' :: Integer - i' = 1 + (i `mod` 1000000) + i' = succ $ i `mod` 1000000 prop_primeFactors, prop_smoothness, prop_eulersTotientP, prop_eulersTotientInequality :: Math.Implementations.PrimeFactorisation.Algorithm -> Integer -> Test.QuickCheck.Property prop_primeFactors algorithm i = Test.QuickCheck.label "prop_primeFactors" $ Data.PrimeFactors.product' (recip 2) {-TODO-} 10 (Math.PrimeFactorisation.primeFactors algorithm i') == i' where i' :: Integer - i' = 1 + (i `mod` 1000000) + i' = succ $ i `mod` 1000000 prop_smoothness algorithm i = Test.QuickCheck.label "prop_smoothness" $ (Math.PrimeFactorisation.smoothness algorithm !! (2 ^ i')) <= (2 :: Integer) where i' :: Integer i' = i `mod` 20 - prop_eulersTotientP algorithm i = Test.QuickCheck.label "prop_eulersTotient" $ Math.PrimeFactorisation.eulersTotient algorithm prime == prime - 1 where + prop_eulersTotientP algorithm i = Test.QuickCheck.label "prop_eulersTotient" $ Math.PrimeFactorisation.eulersTotient algorithm prime == pred prime where prime :: Integer prime = Data.List.genericIndex Data.Numbers.Primes.primes (i `mod` 10000) prop_eulersTotientInequality algorithm i = i `notElem` [2, 6] ==> Test.QuickCheck.label "prop_eulersTotientInequality" $ Math.PrimeFactorisation.eulersTotient algorithm i' >= floor (sqrt $ fromIntegral i' :: Double) where - i' = 1 + (i `mod` 100000) + i' = succ $ i `mod` 100000 prop_eulersTotient, prop_lagrange, prop_multiplicativeOrder, prop_perfectPower :: Math.Implementations.PrimeFactorisation.Algorithm -> Integer -> Integer -> Test.QuickCheck.Property - prop_eulersTotient algorithm i power = Test.QuickCheck.label "prop_eulersTotient" $ Math.PrimeFactorisation.eulersTotient algorithm (base ^ power') == (base ^ (power' - 1)) * (base - 1) where + prop_eulersTotient algorithm i power = Test.QuickCheck.label "prop_eulersTotient" $ Math.PrimeFactorisation.eulersTotient algorithm (base ^ power') == (base ^ pred power') * pred base where base :: Integer base = Data.List.genericIndex Data.Numbers.Primes.primes (i `mod` 8) - power' = 1 + (power `mod` 5) + power' = succ $ power `mod` 5 prop_lagrange algorithm base modulus = gcd base modulus' == 1 ==> Test.QuickCheck.label "prop_lagrange" $ (Math.PrimeFactorisation.eulersTotient algorithm modulus' `rem` Math.MultiplicativeOrder.multiplicativeOrder algorithm base modulus') == 0 where modulus' :: Integer diff --git a/src/Factory/Test/QuickCheck/Primes.hs b/src/Factory/Test/QuickCheck/Primes.hs --- a/src/Factory/Test/QuickCheck/Primes.hs +++ b/src/Factory/Test/QuickCheck/Primes.hs @@ -23,6 +23,8 @@ -} module Factory.Test.QuickCheck.Primes( +-- * Constants +-- defaultAlgorithm, -- * Functions quickChecks, -- isPrime, @@ -32,20 +34,21 @@ import Control.Applicative((<$>)) import qualified Control.DeepSeq import qualified Data.Set +import qualified Factory.Data.PrimeWheel as Data.PrimeWheel import qualified Factory.Math.Implementations.Primality as Math.Implementations.Primality import qualified Factory.Math.Implementations.PrimeFactorisation as Math.Implementations.PrimeFactorisation -import qualified Factory.Math.Implementations.Primes as Math.Implementations.Primes +import qualified Factory.Math.Implementations.Primes.Algorithm as Math.Implementations.Primes.Algorithm import qualified Factory.Math.Primality as Math.Primality import qualified Factory.Math.Primes as Math.Primes import qualified Test.QuickCheck import Test.QuickCheck((==>)) import qualified ToolShed.Defaultable as Defaultable -instance Test.QuickCheck.Arbitrary Math.Implementations.Primes.Algorithm where +instance Test.QuickCheck.Arbitrary Math.Implementations.Primes.Algorithm.Algorithm where arbitrary = Test.QuickCheck.oneof [ - return Math.Implementations.Primes.TurnersSieve, - Math.Implementations.Primes.TrialDivision . (`mod` 10) <$> Test.QuickCheck.arbitrary, - Math.Implementations.Primes.SieveOfEratosthenes . (`mod` 10) <$> Test.QuickCheck.arbitrary + return Math.Implementations.Primes.Algorithm.TurnersSieve, + Math.Implementations.Primes.Algorithm.TrialDivision . (`mod` 10) <$> Test.QuickCheck.arbitrary, + Math.Implementations.Primes.Algorithm.SieveOfEratosthenes . (`mod` 10) <$> Test.QuickCheck.arbitrary ] #if !(MIN_VERSION_QuickCheck(2,1,0)) coarbitrary = undefined --CAVEAT: stops warnings from ghc. @@ -56,23 +59,46 @@ primalityAlgorithm :: Math.Implementations.Primality.Algorithm Math.Implementations.PrimeFactorisation.Algorithm primalityAlgorithm = Defaultable.defaultValue -upperBound :: Math.Implementations.Primes.Algorithm -> Int -> Int -upperBound algorithm i = mod i $ if algorithm == Math.Implementations.Primes.TurnersSieve +upperBound :: Math.Implementations.Primes.Algorithm.Algorithm -> Int -> Int +upperBound algorithm i = mod i $ if algorithm == Math.Implementations.Primes.Algorithm.TurnersSieve then 8192 else 65536 +defaultAlgorithm :: Math.Implementations.Primes.Algorithm.Algorithm +defaultAlgorithm = Defaultable.defaultValue + -- | Defines invariant properties. quickChecks :: IO () quickChecks = Test.QuickCheck.quickCheck `mapM_` [prop_isPrime, prop_isComposite] - >> Test.QuickCheck.quickCheck prop_consistency where - prop_isPrime, prop_isComposite :: Math.Implementations.Primes.Algorithm -> Int -> Test.QuickCheck.Property - prop_isPrime algorithm i = Test.QuickCheck.label "prop_isPrime" . all isPrime . takeWhile (<= (upperBound algorithm i)) $ (Math.Primes.primes algorithm :: [Int]) - + >> Test.QuickCheck.quickCheck prop_consistency + >> Test.QuickCheck.quickCheck prop_rewriteRule + >> Test.QuickCheck.quickCheck `mapM_` [prop_sieveOfAtkin, prop_sieveOfAtkinRewrite] + where + prop_isPrime, prop_isComposite :: Math.Implementations.Primes.Algorithm.Algorithm -> Int -> Test.QuickCheck.Property + prop_isPrime algorithm i = Test.QuickCheck.label "prop_isPrime" . all isPrime . takeWhile (<= upperBound algorithm i) $ (Math.Primes.primes algorithm :: [Int]) prop_isComposite algorithm i = Test.QuickCheck.label "prop_isComposite" . not . any isPrime . Data.Set.toList . Data.Set.difference ( Data.Set.fromList [2 .. upperBound algorithm i] - ) . Data.Set.fromList . takeWhile (<= (upperBound algorithm i)) $ Math.Primes.primes algorithm + ) . Data.Set.fromList . takeWhile (<= upperBound algorithm i) $ Math.Primes.primes algorithm - prop_consistency :: Math.Implementations.Primes.Algorithm -> Math.Implementations.Primes.Algorithm -> Int -> Test.QuickCheck.Property + prop_consistency :: Math.Implementations.Primes.Algorithm.Algorithm -> Math.Implementations.Primes.Algorithm.Algorithm -> Int -> Test.QuickCheck.Property prop_consistency l r i = l /= r ==> Test.QuickCheck.label "prop_consistency" . and . take (i `mod` 4096) $ zipWith (==) (Math.Primes.primes l) (Math.Primes.primes r :: [Int]) + + prop_rewriteRule :: Data.PrimeWheel.NPrimes -> Int -> Test.QuickCheck.Property + prop_rewriteRule wheelSize i = Test.QuickCheck.label "prop_rewriteRule" $ toInteger (Math.Primes.primes (Math.Implementations.Primes.Algorithm.SieveOfEratosthenes wheelSize') !! index :: Int) == (Math.Primes.primes (Math.Implementations.Primes.Algorithm.SieveOfEratosthenes wheelSize') !! index :: Integer) where + wheelSize' = wheelSize `mod` 8 + index = i `mod` 131072 + + prop_sieveOfAtkin, prop_sieveOfAtkinRewrite :: Int -> Test.QuickCheck.Property + prop_sieveOfAtkin i = Test.QuickCheck.label "prop_sieveOfAtkin" $ Math.Primes.primes (Math.Implementations.Primes.Algorithm.SieveOfAtkin prime) !! index == prime where + index = i `mod` 131072 + + prime :: Integer + prime = Math.Primes.primes defaultAlgorithm !! index + + prop_sieveOfAtkinRewrite i = Test.QuickCheck.label "prop_sieveOfAtkin'" $ Math.Primes.primes (Math.Implementations.Primes.Algorithm.SieveOfAtkin $ fromIntegral prime) !! index == prime where + index = i `mod` 131072 + + prime :: Int + prime = Math.Primes.primes defaultAlgorithm !! index diff --git a/src/Factory/Test/QuickCheck/Probability.hs b/src/Factory/Test/QuickCheck/Probability.hs --- a/src/Factory/Test/QuickCheck/Probability.hs +++ b/src/Factory/Test/QuickCheck/Probability.hs @@ -45,7 +45,7 @@ ) where prop_normalDistribution :: System.Random.RandomGen g => g -> (Double, Double) -> Test.QuickCheck.Property prop_normalDistribution randomGen (mean, variance) = variance' /= 0 ==> Test.QuickCheck.label "prop_normalDistribution" . Pair.both . Pair.mirror ( - (< (0.05 :: Double)) . abs --Tolerance. + (< (0.1 :: Double)) . abs --Generous tolerance. ) . ( Math.Statistics.getMean &&& pred . Math.Statistics.getStandardDeviation ) . map ( diff --git a/src/Factory/Test/QuickCheck/QuickChecks.hs b/src/Factory/Test/QuickCheck/QuickChecks.hs --- a/src/Factory/Test/QuickCheck/QuickChecks.hs +++ b/src/Factory/Test/QuickCheck/QuickChecks.hs @@ -30,6 +30,7 @@ import qualified Factory.Test.QuickCheck.Hyperoperation import qualified Factory.Test.QuickCheck.Interval import qualified Factory.Test.QuickCheck.MonicPolynomial +import qualified Factory.Test.QuickCheck.PerfectPower import qualified Factory.Test.QuickCheck.Pi import qualified Factory.Test.QuickCheck.Polynomial import qualified Factory.Test.QuickCheck.Power @@ -49,6 +50,7 @@ >> putStrLn "Hyperoperation" >> Factory.Test.QuickCheck.Hyperoperation.quickChecks >> putStrLn "Interval" >> Factory.Test.QuickCheck.Interval.quickChecks >> putStrLn "MonicPolynomial" >> Factory.Test.QuickCheck.MonicPolynomial.quickChecks + >> putStrLn "PerfectPower" >> Factory.Test.QuickCheck.PerfectPower.quickChecks >> putStrLn "Pi" >> Factory.Test.QuickCheck.Pi.quickChecks >> putStrLn "Polynomial" >> Factory.Test.QuickCheck.Polynomial.quickChecks >> putStrLn "Power" >> Factory.Test.QuickCheck.Power.quickChecks diff --git a/src/Factory/Test/QuickCheck/Radix.hs b/src/Factory/Test/QuickCheck/Radix.hs --- a/src/Factory/Test/QuickCheck/Radix.hs +++ b/src/Factory/Test/QuickCheck/Radix.hs @@ -42,5 +42,5 @@ base = (b `mod` 73) - 36 prop_digitalRoot (_, n) = Test.QuickCheck.label "prop_digitalRoot" $ Math.Radix.digitalRoot n' == 9 where - n' = 9 * (1 + abs n) + n' = 9 * succ (abs n) diff --git a/src/Factory/Test/QuickCheck/SquareRoot.hs b/src/Factory/Test/QuickCheck/SquareRoot.hs --- a/src/Factory/Test/QuickCheck/SquareRoot.hs +++ b/src/Factory/Test/QuickCheck/SquareRoot.hs @@ -61,7 +61,7 @@ prop_accuracy, prop_factorable, prop_perfectSquare :: Testable prop_accuracy (algorithm, decimalDigits, operand) = Test.QuickCheck.label "prop_accuracy" . (>= requiredDecimalDigits) . Math.SquareRoot.getAccuracy operand' $ Math.SquareRoot.squareRoot algorithm requiredDecimalDigits operand' where requiredDecimalDigits :: Math.Precision.DecimalDigits - requiredDecimalDigits = 1 + (decimalDigits `mod` 1024) + requiredDecimalDigits = succ $ decimalDigits `mod` 1024 operand' :: Data.Ratio.Rational operand' = abs operand @@ -76,14 +76,14 @@ ) ) / Math.SquareRoot.squareRoot algorithm requiredDecimalDigits operand' where requiredDecimalDigits :: Math.Precision.DecimalDigits - requiredDecimalDigits = 1 + (decimalDigits `mod` 1024) + requiredDecimalDigits = succ $ decimalDigits `mod` 1024 operand' :: Data.Ratio.Rational - operand' = 1 + abs operand + operand' = succ $ abs operand prop_perfectSquare (algorithm, decimalDigits, operand) = Test.QuickCheck.label "prop_perfectSquare" . Math.SquareRoot.isPrecise perfectSquare $ Math.SquareRoot.squareRoot algorithm requiredDecimalDigits perfectSquare where requiredDecimalDigits :: Math.Precision.DecimalDigits - requiredDecimalDigits = 1 + (decimalDigits `mod` 32768) + requiredDecimalDigits = succ $ decimalDigits `mod` 32768 operand', perfectSquare :: Data.Ratio.Rational operand' = (abs (Data.Ratio.numerator operand) `min` (2 ^ (32 :: Int))) % (abs (Data.Ratio.denominator operand) `min` (2 ^ (32 :: Int))) --Avoid floating-point rounding-errors in 'Math.SquareRoot.rSqrt'. diff --git a/src/Factory/Test/QuickCheck/Statistics.hs b/src/Factory/Test/QuickCheck/Statistics.hs --- a/src/Factory/Test/QuickCheck/Statistics.hs +++ b/src/Factory/Test/QuickCheck/Statistics.hs @@ -47,10 +47,10 @@ prop_nC0 algorithm n = Test.QuickCheck.label "prop_nC0" $ Math.Statistics.nCr algorithm (abs n) 0 == 1 prop_nC1 algorithm i = Test.QuickCheck.label "prop_nC1" $ Math.Statistics.nCr algorithm n 1 == n where - n = 1 + abs i + n = succ $ abs i prop_sum algorithm i = Test.QuickCheck.label "prop_sum" $ sum (Math.Statistics.nCr algorithm n `map` [0 .. n]) == 2 ^ n where - n = 1 + abs i + n = succ $ abs i prop_symmetry, prop_prime :: Math.Implementations.Factorial.Algorithm -> (Integer, Integer) -> Test.QuickCheck.Property prop_symmetry algorithm (i, j) = Test.QuickCheck.label "prop_symmetry" $ Math.Statistics.nCr algorithm n r == Math.Statistics.nCr algorithm n (n - r) where @@ -64,7 +64,7 @@ prop_nP0 n = Test.QuickCheck.label "prop_nP0" $ Math.Statistics.nPr (abs n) 0 == 1 prop_nP1 i = Test.QuickCheck.label "prop_nP1" $ Math.Statistics.nPr n 1 == n where - n = 1 + abs i + n = succ $ abs i prop_zeroVariance, prop_zeroAverageAbsoluteDeviation :: Data.Ratio.Rational -> Test.QuickCheck.Property prop_zeroVariance x = Test.QuickCheck.label "prop_zeroVariance" $ Math.Statistics.getVariance (replicate 32 x) == (0 :: Data.Ratio.Rational) diff --git a/src/Main.hs b/src/Main.hs --- a/src/Main.hs +++ b/src/Main.hs @@ -25,7 +25,8 @@ -} module Main( --- * Type-classes +-- * Types +-- ** Type-synonyms -- CommandLineAction, -- * Functions main @@ -33,6 +34,7 @@ import qualified Data.List import qualified Data.Ratio +import qualified Data.Version import qualified Distribution.Package import qualified Distribution.Text import qualified Distribution.Version @@ -40,7 +42,7 @@ import qualified Factory.Math.Implementations.Factorial as Math.Implementations.Factorial import qualified Factory.Math.Implementations.Primality as Math.Implementations.Primality import qualified Factory.Math.Implementations.PrimeFactorisation as Math.Implementations.PrimeFactorisation -import qualified Factory.Math.Implementations.Primes as Math.Implementations.Primes +import qualified Factory.Math.Implementations.Primes.Algorithm as Math.Implementations.Primes.Algorithm import qualified Factory.Math.Implementations.SquareRoot as Math.Implementations.SquareRoot import qualified Factory.Test.CommandOptions as Test.CommandOptions import qualified Factory.Test.Performance.Factorial as Test.Performance.Factorial @@ -52,6 +54,7 @@ import qualified Factory.Test.Performance.SquareRoot as Test.Performance.SquareRoot import qualified Factory.Test.Performance.Statistics as Test.Performance.Statistics import qualified Factory.Test.QuickCheck.QuickChecks as Test.QuickCheck.QuickChecks +import qualified Paths_factory as Paths --Either local stub, or package-instance autogenerated by 'Setup.hs build'. import qualified System import qualified System.Console.GetOpt as G import qualified System.IO @@ -93,7 +96,7 @@ G.Option "" ["piPerformanceGraph"] (piPerformanceGraph `G.ReqArg` "(Math.Pi.Category, Double, Math.Precision.DecimalDigits)") "Test the performance of 'Math.Pi.openI', with an exponential precision-requirement (of the specified exponent), up to the specified limit.", G.Option "" ["primeFactorsPerformance"] (primeFactorsPerformance `G.ReqArg` "(Math.Implementations.PrimeFactorisation.Algorithm, Integer)") "Test the performance of 'Math.PrimeFactorisation.primeFactors'.", G.Option "" ["primeFactorsPerformanceGraph"] (primeFactorsPerformanceGraph `G.ReqArg` "(Math.Implementations.PrimeFactorisation.Algorithm, Int)") "Test the performance of 'Math.PrimeFactorisation.primeFactors', on the specified number of odd integers from the Fibonacci-sequence.", - G.Option "" ["primesPerformance"] (primesPerformance `G.ReqArg` "(Math.Implementations.Primes.Algorithm, Int)") "Test the performance of 'Math.Primes.primes'.", + G.Option "" ["primesPerformance"] (primesPerformance `G.ReqArg` "(Math.Implementations.Primes.Algorithm.Algorithm, Int)") "Test the performance of 'Math.Primes.primes'.", G.Option "" ["squareRootPerformance"] (squareRootPerformance `G.ReqArg` "(Math.Implementations.SquareRoot.Algorithm, Data.Ratio.Rational, DecimalDigits)") "Test the performance of 'Math.SquareRoot.squareRoot'.", G.Option "" ["squareRootPerformanceGraph"] (squareRootPerformanceGraph `G.ReqArg` "(Math.Implementations.SquareRoot.Algorithm, Data.Ratio.Rational)") "Test the performance of 'Math.SquareRoot.squareRoot', with an exponentially increasing precision-requirement.", G.Option "" ["verbose"] (G.NoArg $ return {-to IO-monad-} . Test.CommandOptions.setVerbose) ("Provide additional information where available; default '" ++ show (Test.CommandOptions.verbose Defaultable.defaultValue) ++ "'."), @@ -106,7 +109,7 @@ packageIdentifier :: Distribution.Package.PackageIdentifier packageIdentifier = Distribution.Package.PackageIdentifier { Distribution.Package.pkgName = Distribution.Package.PackageName "factory", - Distribution.Package.pkgVersion = Distribution.Version.Version [0, 1, 0, 3] [] + Distribution.Package.pkgVersion = Distribution.Version.Version (Data.Version.versionBranch Paths.version) [] } printUsage = System.IO.hPutStrLn System.IO.stderr usage >> System.exitWith System.ExitSuccess @@ -181,15 +184,16 @@ CAVEAT: fragile. -} case algorithm of - Math.Implementations.Primes.SieveOfEratosthenes n -> Test.Performance.Primes.primesPerformance $ Math.Implementations.Primes.SieveOfEratosthenes n - _ -> Test.Performance.Primes.primesPerformance algorithm + Math.Implementations.Primes.Algorithm.SieveOfEratosthenes wheelSize -> Test.Performance.Primes.primesPerformance $ Math.Implementations.Primes.Algorithm.SieveOfEratosthenes wheelSize + Math.Implementations.Primes.Algorithm.SieveOfAtkin maxPrime -> Test.Performance.Primes.primesPerformance $ Math.Implementations.Primes.Algorithm.SieveOfAtkin maxPrime + _ -> Test.Performance.Primes.primesPerformance algorithm ) index :: IO ( Double, -- Integer - Int --Exploits rewrite-rule in "Math.Implementations.Primes". + Int --Exploits rewrite-rules in "Math.Implementations.Primes.*". ) ) >>= print >> System.exitWith System.ExitSuccess where - algorithm :: Math.Implementations.Primes.Algorithm + algorithm :: Math.Implementations.Primes.Algorithm.Algorithm (algorithm, index) = read arg squareRootPerformance arg _ = Test.Performance.SquareRoot.squareRootPerformance algorithm operand decimalDigits >>= print >> System.exitWith System.ExitSuccess where