diff --git a/changelog b/changelog --- a/changelog +++ b/changelog @@ -3,21 +3,33 @@ 0.0.0.1 * First version of the package. 0.0.0.2 - * Created modules; "Factory.Test.QuickCheck.Bounds", "Factory.Math.Implementations.Pi.Borwein.*" and "Factory.Test.Performance.Statistics". + * Created the modules; "Factory.Test.QuickCheck.Bounds", "Factory.Math.Implementations.Pi.Borwein.*" and "Factory.Test.Performance.Statistics". * Created a new module "Factory.Data.PrimeFactors", and migrated definitions from both "Factory.Math.PrimeFactorisation" and "Factory.Math.Implementations.PrimeFactorisation". - * Created class 'Factory.Math.Factorial.Factorial' and new module "Factory.Math.Implementations.Factorial". - Moved existing implementation (Bisection) into new module, with new implementation (PrimeFactorisation). - * Added function 'Factory.Math.Summation.sumR'. - * Added a parameter to functions 'Factory.Math.DivideAndConquer.divideAndConquer' and 'Factory.Data.Bounds.divideAndConquer' to permit asymmetric bisection. + * Created the class 'Factory.Math.Factorial.Factorial' and a new module "Factory.Math.Implementations.Factorial". + Moved existing implementation (Bisection) into the new module, with a new implementation (PrimeFactorisation). + * Added the function 'Factory.Math.Summation.sumR'. + * Added a parameter to the functions 'Factory.Math.DivideAndConquer.divideAndConquer' and 'Factory.Data.Bounds.divideAndConquer', to permit asymmetric bisection. * Added methods to class "Factory.Math.Pi.Algorithm" to permit the retrieval of /Pi/ as a 'Rational' or a 'String'. - * Renamed 'Factory.Math.Precision.capPrecision' to 'Factory.Math.Precision.simplify'. - * Removed module "Factory.Test.Performance.Exponential". - * Removed function 'Factory.Math.Power.raise', which was no more efficient than ghc's implementation of '(^)'. + * Renamed the function 'Factory.Math.Precision.capPrecision' to 'Factory.Math.Precision.simplify'. + * Removed the module "Factory.Test.Performance.Exponential". + * Removed the function 'Factory.Math.Power.raise', which was no more efficient than ghc's implementation of '(^)'. * Uploaded to . 0.1.0.0 * Amended 'factory.cabal' to more correctly specify the dependency on package 'toolshed'. - * Added module "Factory.Math.Probability". - * Renamed module "Factory.Data.Bounds" to "Factory.Data.Interval", - and added functions; 'Factory.Data.Interval.precisely', 'Factory.Data.Interval.shift', 'Factory.Data.Interval.closedUnitInterval'. - * Guarded 'eager-blackholing' flag in cabal file. + * Added the module "Factory.Math.Probability". + * Renamed the module "Factory.Data.Bounds" to "Factory.Data.Interval", + and added the functions; 'Factory.Data.Interval.precisely', 'Factory.Data.Interval.shift', 'Factory.Data.Interval.closedUnitInterval'. + * Guarded 'eager-blackholing' flag in /cabal/ file. +0.1.0.1 + * Renamed classes "Factory.Math.[Primality, Pi, Factorial, SquareRoot, PrimeFactorisation].Algorithm" to "Factory.Math.[Primality, Pi, Factorial, SquareRoot, PrimeFactorisation].Algorithmic", to distinguish them from the data-types which implement them. + * Added the modules "Factory.Math.Hyperoperation", "Factory.Test.QuickCheck.Hyperoperation" and "Factory.Test.Performance.Hyperoperation". + * Added the modules "Factory.Math.Primes", "Factory.Math.Implementation.Primes", "Factory.Test.Performance.Primes", "Factory.Test.QuickCheck.Primes" and "Factory.Data.PrimeWheel". + * Added the function 'Factory.Math.PrimeFactorisation.squareFree'. + * Added rewrite-rules to specialise 'Factory.Math.Power.isPerfectPower' for type-parameter='Int'. + * Recoded "Factory.Math.Radix" to the interface "Data.Array.IArray.IArray", rather than the data-type "Data.Array.Array". +0.1.0.2 + * Added 'Factory.Math.Primes.primorial'. + * Altered 'Factory.Math.Implementations.Primes.trialDivision' to take an integer defining the size of a 'Factory.Data.PrimeWheel', from which candidates are extracted. + * Removed the command-line option 'primesPerformanceGraph', which appears to memoise data from previous tests. + 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.0 +Version: 0.1.0.2 Cabal-Version: >= 1.6 Copyright: (C) 2011 Dr. Alistair Ward License: GPL @@ -36,12 +36,14 @@ Factory.Data.Monomial Factory.Data.Polynomial Factory.Data.PrimeFactors + Factory.Data.PrimeWheel Factory.Data.QuotientRing Factory.Data.Ring Factory.Math.ArithmeticGeometricMean Factory.Math.DivideAndConquer Factory.Math.Factorial Factory.Math.Fibonacci + Factory.Math.Hyperoperation Factory.Math.Implementations.Factorial Factory.Math.Implementations.Pi.AGM.Algorithm Factory.Math.Implementations.Pi.AGM.BrentSalamin @@ -66,6 +68,7 @@ Factory.Math.Implementations.Pi.Spigot.Spigot Factory.Math.Implementations.Primality Factory.Math.Implementations.PrimeFactorisation + Factory.Math.Implementations.Primes Factory.Math.Implementations.SquareRoot Factory.Math.MultiplicativeOrder Factory.Math.Pi @@ -73,6 +76,7 @@ Factory.Math.Precision Factory.Math.Primality Factory.Math.PrimeFactorisation + Factory.Math.Primes Factory.Math.Probability Factory.Math.Radix Factory.Math.SquareRoot @@ -107,13 +111,16 @@ Other-modules: Factory.Test.CommandOptions Factory.Test.Performance.Factorial + Factory.Test.Performance.Hyperoperation Factory.Test.Performance.Pi Factory.Test.Performance.Primality Factory.Test.Performance.PrimeFactorisation + Factory.Test.Performance.Primes Factory.Test.Performance.SquareRoot Factory.Test.Performance.Statistics Factory.Test.QuickCheck.ArithmeticGeometricMean Factory.Test.QuickCheck.Factorial + Factory.Test.QuickCheck.Hyperoperation Factory.Test.QuickCheck.Interval Factory.Test.QuickCheck.MonicPolynomial Factory.Test.QuickCheck.Pi @@ -121,6 +128,7 @@ Factory.Test.QuickCheck.Power Factory.Test.QuickCheck.Primality Factory.Test.QuickCheck.PrimeFactorisation + Factory.Test.QuickCheck.Primes Factory.Test.QuickCheck.Probability Factory.Test.QuickCheck.QuickChecks Factory.Test.QuickCheck.Radix diff --git a/makefile b/makefile --- a/makefile +++ b/makefile @@ -50,7 +50,7 @@ clean: runhaskell Setup.hs $@ - find src -type f \( -name '*.o' -o -name '*.hi' \) -delete + find src -type f \( -name '*.hc' -o -name '*.hcr' -o -name '*.hi' -o -name '*.o' \) -delete help: @grep '^[a-zA-Z].*:' makefile | sed -e 's/:.*//' diff --git a/src/Factory/Data/Monomial.hs b/src/Factory/Data/Monomial.hs --- a/src/Factory/Data/Monomial.hs +++ b/src/Factory/Data/Monomial.hs @@ -49,7 +49,6 @@ ) where import qualified Control.Arrow -import qualified Factory.Math.Power as Math.Power infix 4 <=> --Same as (==). infix 4 =~ --Same as (==). @@ -105,7 +104,7 @@ -- | Square the specified 'Monomial'. square :: (Num c, Num e) => Monomial c e -> Monomial c e -square (c, e) = (Math.Power.square c, 2 * e) +square (c, e) = (c ^ (2 :: Int), 2 * e) -- | Double the specified 'Monomial'. {-# INLINE double #-} diff --git a/src/Factory/Data/PrimeWheel.hs b/src/Factory/Data/PrimeWheel.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Data/PrimeWheel.hs @@ -0,0 +1,178 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines a /prime-wheel/, for use in prime-number generation; . +-} + +module Factory.Data.PrimeWheel( +-- * Types +-- ** Type-synonyms + Distance, + PrimeMultiples, +-- Repository, +-- ** Data-types + PrimeWheel(getPrimeComponents), +-- * Functions +-- findCoprimes, + generatePrimeMultiples, + roll, + rotate, +-- ** Constructors + mkPrimeWheel, +-- ** Query + getCircumference, + getSpokeCount +) where + +import Control.Arrow((&&&), (***)) +import qualified Data.IntMap +import qualified Data.List + +{- | + * A conceptual /wheel/, with irregularly spaced spokes; . + + * On being rolled, the trace of the spokes, identifies candidates which are /coprime/ to those primes from which the /wheel/ was composed. + + * One can alternatively view this as a set of vertical nested rings, each with a /prime circumference/, and touching at its lowest point. + Each has a single mark on its /circumference/, which when rolled identifies multiples of that /circumference/. + When the complete set is rolled, from the state where all marks are coincident, all multiples of the set of primes, are traced. + + * CAVEAT: The distance required to return this state, the /circumference/ grows rapidly, with the number of primes: + +> zip [0 ..] . scanl (*) 1 $ [2,3,5,7,11,13,17,19,23,29,31] +> [(0,1),(1,2),(2,6),(3,30),(4,210),(5,2310),(6,30030),(7,510510),(8,9699690),(9,223092870),(10,6469693230),(11,200560490130)] + + * The number of spokes also grows rapidly with the number of primes: + +> zip [0 ..] . scanl (*) 1 . map pred $ [2,3,5,7,11,13,17,19,23,29,31] +> [(0,1),(1,1),(2,2),(3,8),(4,48),(5,480),(6,5760),(7,92160),(8,1658880),(9,36495360),(10,1021870080),(11,30656102400)] +-} +data PrimeWheel i = MkPrimeWheel { + getPrimeComponents :: [i], -- ^ Accessor: the ordered sequence of initial primes, from which the /wheel/ was composed. + getSpokeGaps :: [i] -- ^ Accessor: the sequence of spoke-gaps, the sum of which equals its /circumference/. +} deriving Show + +-- | The /circumference/ of the specified 'PrimeWheel'. +getCircumference :: Num n => PrimeWheel n -> n +getCircumference = product . getPrimeComponents + +-- | The number of spokes in the specified 'PrimeWheel'. +getSpokeCount:: Integral i => PrimeWheel i -> i +getSpokeCount = foldr ((*) . pred) 1 . getPrimeComponents + +-- | An infinite increasing sequence, of the multiples of a specific prime. +type PrimeMultiples i = [i] + +-- | Defines a container for the 'PrimeMultiples' of an initial short sequence of primes. +type Repository = Data.IntMap.IntMap [PrimeMultiples 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. + + * 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 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 candidate found repository = case Data.IntMap.lookup candidate repository of + Just primeMultiplesList -> sieve' found $ foldr insert ( + Data.IntMap.delete candidate repository --Remove the matched prime-multiple. + ) primeMultiplesList + Nothing -> let + found' = succ found + (key : values) = iterate (+ gap * candidate) $ candidate ^ (2 :: Int) --Generate a sequence of prime-multiples, starting from its square. + in candidate : sieve' found' (if found' >= required then repository else Data.IntMap.insert key [values] repository) + where + gap :: Int + gap = 2 --For efficiency, only sieve odd integers. + + sieve' :: Int -> Repository -> [Int] + sieve' = sieve $ candidate + gap --Tail-recurse. + + insert :: PrimeMultiples Int -> Repository -> Repository + insert [] = error "Factory.Data.PrimeWheel.insert:\tnull list" + insert (key : values) = Data.IntMap.insertWith (++) key [values] + +{- | + * 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 /circumference/ should be approximately the /square-root/ of the number of integers it will be required to sieve. + + * The sequence of gaps between spokes on the /wheel/ is /symmetrical under reflection/; + though two values lie /on/ the axis, that aren't part of this symmetry. Eg: + +> nPrimes Gaps +> ====== ==== +> 0 [1] +> 1 [2] --The terminal gap for all subsequent wheels is '2'; [(succ circumference `mod` circumference) - (pred circumference `mod` circumference)]. +> 2 [4,2] --Both points are on the axis, so the symmetry isn't yet clear. +> 3 [6,4,2,4,2,4,6,2] +> 4 [10,2,4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,8,6,4,6,2,4,6,2,6,6,4,2,4,6,2,6,4,2,4,2,10,2] + + Exploitation of this property has proved counter-productive, probably because it requires /strict evaluation/, + exposing the user to the full cost of inadvertently choosing a /wheel/, which in practice, is rotated less than once. +-} +mkPrimeWheel :: Integral i => Int -> PrimeWheel i +mkPrimeWheel 0 = MkPrimeWheel [] [1] +mkPrimeWheel nPrimes + | nPrimes < 0 = error $ "Factory.Data.PrimeWheel.mkPrimeWheel: unable to construct a wheel from " ++ show nPrimes ++ " primes" + | otherwise = primeWheel + where + (primeComponents, coprimeCandidates) = (map fromIntegral *** map fromIntegral . Data.List.genericTake (getSpokeCount primeWheel)) $ findCoprimes nPrimes + primeWheel = MkPrimeWheel primeComponents $ zipWith (-) coprimeCandidates $ 1 : coprimeCandidates --Measure the gaps between candidate primes. + +-- | Couples a candidate prime with a /rolling wheel/, to define the distance rolled. +type Distance i = (i, [i]) + +-- | Generates a new candidate prime, from a /rolling wheel/, and the current candidate. +rotate :: Integral i => Distance i -> Distance i +rotate (candidate, rollingWheel) = (candidate +) . head &&& tail $ rollingWheel + +{-# INLINE rotate #-} + +-- | Generate an infinite, increasing sequence of candidate primes, from the specified /wheel/. +roll :: Integral i => PrimeWheel i -> [Distance i] +roll primeWheel = tail $ iterate rotate (1, cycle $ getSpokeGaps primeWheel) + +{- | + * Generates multiples of the specified prime, starting from its /square/, + skipping those multiples of the low primes from which the specified 'PrimeWheel' was composed, + and which therefore, the /wheel/ won't generate as candidates. Eg: + +> Prime Rotating PrimeWheel 3 Output +> ===== ===================== ====== +> 7 [4,2,4,2,4,6,2,6] [49,77,91,119,133,161,203,217,259 ..] +> 11 [2,4,2,4,6,2,6,4] [121,143,187,209,253,319,341,407 ..] +> 13 [4,2,4,6,2,6,4,2] [169,221,247,299,377,403,481,533,559 ..] +-} +generatePrimeMultiples :: Integral i + => i -- ^ The /prime/. + -> [i] -- ^ A /rolling wheel/, the track of which, delimits the gaps between /coprime/ candidates. + -> [i] +generatePrimeMultiples prime = scanl (\accumulator -> (+ accumulator) . (* prime)) (prime ^ (2 :: Int)) + +{-# INLINE generatePrimeMultiples #-} + 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 @@ -66,7 +66,7 @@ getGeometricMean = snd -- | Returns an infinite list which converges on the /Arithmetic-geometric mean/. -convergeToAGM :: Math.SquareRoot.Algorithm squareRootAlgorithm => squareRootAlgorithm -> Math.Precision.DecimalDigits -> AGM -> [AGM] +convergeToAGM :: Math.SquareRoot.Algorithmic squareRootAlgorithm => squareRootAlgorithm -> Math.Precision.DecimalDigits -> AGM -> [AGM] convergeToAGM squareRootAlgorithm decimalDigits agm | decimalDigits <= 0 = error $ "Factory.Math.ArithmeticGeometricMean.convergeToAGM:\tinvalid number of decimal digits; " ++ show decimalDigits | not $ isValid agm = error $ "Factory.Math.ArithmeticGeometricMean.convergeToAGM:\tboth means must be positive for a real geometric mean; " ++ show agm diff --git a/src/Factory/Math/Factorial.hs b/src/Factory/Math/Factorial.hs --- a/src/Factory/Math/Factorial.hs +++ b/src/Factory/Math/Factorial.hs @@ -28,10 +28,10 @@ module Factory.Math.Factorial( -- * Type-classes - Algorithm(..) + Algorithmic(..) ) where -- | Defines the methods expected of a /factorial/-algorithm. -class Algorithm algorithm where +class Algorithmic algorithm where factorial :: Integral i => algorithm -> i -> i diff --git a/src/Factory/Math/Hyperoperation.hs b/src/Factory/Math/Hyperoperation.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Hyperoperation.hs @@ -0,0 +1,113 @@ +{- + 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@] Provides various /hyperoperations/; . +-} + +module Factory.Math.Hyperoperation( +-- * Types +-- ** Type-synonyms + Base, + HyperExponent, +-- * Constants + succession, + addition, + multiplication, + exponentiation, + tetration, + pentation, + hexation, +-- * Functions + hyperoperation, + ackermannPeter, + powerTower, +-- ** Predicates + areCoincidental +) where + +import qualified Data.List + +{- | + * Merely to enhance self-documentation. + + * CAVEAT: whilst it may appear that 'Base' could be non-'Integral', the recursive definition for /hyper-exponents/ above 'tetration', prevents this. +-} +type Base = Integer + +{- | + * Merely to enhance self-documentation. + + * CAVEAT: whilst 'Base' and 'HyperExponent' can be independent types for both 'exponentiation' and 'tetration', they interact for other /hyper-exponents/. +-} +type HyperExponent = Base + +succession, addition, multiplication, exponentiation, tetration, pentation, hexation :: Int --Arbitrarily. +(succession : addition : multiplication : exponentiation : tetration : pentation : hexation : _) = [0 ..] + +{- | + * Returns the /power-tower/ of the specified /base/; . + + * A synonym for /tetration/; + , + . +-} +powerTower :: (Integral base, Integral hyperExponent) => base -> hyperExponent -> base +powerTower 0 hyperExponent + | odd hyperExponent = 0 + | otherwise = 1 +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 + | otherwise = Data.List.genericIndex (iterate (base ^) 1) hyperExponent + +-- | The /hyperoperation/-sequence; . +hyperoperation :: Integral rank => rank -> Base -> HyperExponent -> Base +hyperoperation rank base hyperExponent + | rank < fromIntegral succession = error $ "Factory.Math.Hyperoperation.hyperoperation:\tundefined for rank; " ++ show rank + | hyperExponent < 0 = error $ "Factory.Math.Hyperoperation.hyperoperation:\tundefined for hyper-exponent; " ++ show hyperExponent + | otherwise = rank ^# hyperExponent + where + (^#) :: Integral rank => rank -> HyperExponent -> Base + r ^# 0 = case r of + 1 {-addition-} -> base + 2 {-multiplication-} -> 0 + _ -> 1 + r ^# e = case r of + 0 {-succession-} -> succ {-fromIntegral-} e + 1 {-addition-} -> base + {-fromIntegral-} e + 2 {-multiplication-} -> base * {-fromIntegral-} e + 3 {-exponentiation-} -> base ^ e + 4 {-tetration-} -> base `powerTower` e + _ + | e' == e -> tetration ^# e' --To which it would otherwise be reduced by laborious recursion. + | otherwise -> pred r ^# e' + where + e' = {-fromIntegral $-} r ^# pred e + +-- | The /Ackermann-Peter/-function; . +ackermannPeter :: Integral rank => rank -> HyperExponent -> Base +ackermannPeter rank = (+ negate 3) . hyperoperation rank 2 {-base-} . (+ 3) + +-- | 'True' if @hyperoperation base hyperExponent@ has the same value for each specified 'rank'. +areCoincidental :: Integral rank => Base -> HyperExponent -> [rank] -> Bool +areCoincidental _ _ [] = True +areCoincidental _ _ [_] = True +areCoincidental base hyperExponent ranks = all (== h) hs where + (h : hs) = map (\rank -> hyperoperation rank base hyperExponent) ranks + 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 @@ -19,7 +19,7 @@ [@DESCRIPTION@] - * Provides implementations of the class 'Math.Factorial.Algorithm'. + * Provides implementations of the class 'Math.Factorial.Algorithmic'. * Provides additional functions related to /factorials/, but which depends on a specific implementation, and which therefore can't be accessed throught the class-interface. @@ -61,7 +61,7 @@ instance Defaultable.Defaultable Algorithm where defaultValue = Bisection -instance Math.Factorial.Algorithm Algorithm where +instance Math.Factorial.Algorithmic Algorithm where factorial algorithm n | n < 2 = 1 | otherwise = case algorithm of diff --git a/src/Factory/Math/Implementations/Pi/AGM/Algorithm.hs b/src/Factory/Math/Implementations/Pi/AGM/Algorithm.hs --- a/src/Factory/Math/Implementations/Pi/AGM/Algorithm.hs +++ b/src/Factory/Math/Implementations/Pi/AGM/Algorithm.hs @@ -37,6 +37,6 @@ instance Defaultable.Defaultable squareRootAlgorithm => Defaultable.Defaultable (Algorithm squareRootAlgorithm) where defaultValue = BrentSalamin Defaultable.defaultValue -instance Math.SquareRoot.Algorithm squareRootAlgorithm => Math.Pi.Algorithm (Algorithm squareRootAlgorithm) where +instance Math.SquareRoot.Algorithmic squareRootAlgorithm => Math.Pi.Algorithmic (Algorithm squareRootAlgorithm) where openR (BrentSalamin squareRootAlgorithm) = Math.Implementations.Pi.AGM.BrentSalamin.openR squareRootAlgorithm diff --git a/src/Factory/Math/Implementations/Pi/AGM/BrentSalamin.hs b/src/Factory/Math/Implementations/Pi/AGM/BrentSalamin.hs --- a/src/Factory/Math/Implementations/Pi/AGM/BrentSalamin.hs +++ b/src/Factory/Math/Implementations/Pi/AGM/BrentSalamin.hs @@ -56,7 +56,7 @@ > => 4*a[N]^2 / (1 - sum [2^(n+1) * (a[n-1]^2 - g[n-1]^2)]) -} -openR :: Math.SquareRoot.Algorithm squareRootAlgorithm => squareRootAlgorithm -> Math.Precision.DecimalDigits -> Data.Ratio.Rational +openR :: Math.SquareRoot.Algorithmic squareRootAlgorithm => squareRootAlgorithm -> Math.Precision.DecimalDigits -> Data.Ratio.Rational openR squareRootAlgorithm decimalDigits = uncurry (/) . ( Math.Power.square . uncurry (+) . last &&& negate . pred . sum . zipWith (*) (iterate (* 2) 1) . map (Math.Power.square . Math.ArithmeticGeometricMean.spread) ) . take ( diff --git a/src/Factory/Math/Implementations/Pi/BBP/Algorithm.hs b/src/Factory/Math/Implementations/Pi/BBP/Algorithm.hs --- a/src/Factory/Math/Implementations/Pi/BBP/Algorithm.hs +++ b/src/Factory/Math/Implementations/Pi/BBP/Algorithm.hs @@ -41,7 +41,7 @@ instance Defaultable.Defaultable Algorithm where defaultValue = Base65536 -instance Math.Pi.Algorithm Algorithm where +instance Math.Pi.Algorithmic Algorithm where openR Base65536 = Math.Implementations.Pi.BBP.Implementation.openR Math.Implementations.Pi.BBP.Base65536.series openR Bellard = Math.Implementations.Pi.BBP.Implementation.openR Math.Implementations.Pi.BBP.Bellard.series diff --git a/src/Factory/Math/Implementations/Pi/BBP/Base65536.hs b/src/Factory/Math/Implementations/Pi/BBP/Base65536.hs --- a/src/Factory/Math/Implementations/Pi/BBP/Base65536.hs +++ b/src/Factory/Math/Implementations/Pi/BBP/Base65536.hs @@ -31,7 +31,7 @@ -- | Defines the parameters of this specific series. series :: Math.Implementations.Pi.BBP.Series.Series series = Math.Implementations.Pi.BBP.Series.MkSeries { - Math.Implementations.Pi.BBP.Series.numerators = zipWith ($) (concat $ repeat [id, id, id, negate]) $ map (2 ^) [15 :: Integer, 14, 14, 12, 11, 10, 10, 8, 7, 6, 6, 4, 3, 2, 2, 0], + Math.Implementations.Pi.BBP.Series.numerators = zipWith ($) (cycle [id, id, id, negate]) $ map (2 ^) [15 :: Integer, 14, 14, 12, 11, 10, 10, 8, 7, 6, 6, 4, 3, 2, 2, 0], Math.Implementations.Pi.BBP.Series.getDenominators = \i -> map (32 * fromIntegral i +) [2, 3, 4, 7, 10, 11, 12, 15, 18, 19, 20, 23, 26, 27, 28, 31], Math.Implementations.Pi.BBP.Series.seriesScalingFactor = recip $ 2 ^ (13 :: Int), Math.Implementations.Pi.BBP.Series.base = 2 ^ (16 :: Int) diff --git a/src/Factory/Math/Implementations/Pi/Borwein/Algorithm.hs b/src/Factory/Math/Implementations/Pi/Borwein/Algorithm.hs --- a/src/Factory/Math/Implementations/Pi/Borwein/Algorithm.hs +++ b/src/Factory/Math/Implementations/Pi/Borwein/Algorithm.hs @@ -49,8 +49,8 @@ defaultValue = Borwein1993 Defaultable.defaultValue Defaultable.defaultValue instance ( - Math.SquareRoot.Algorithm squareRootAlgorithm, - Math.Factorial.Algorithm factorialAlgorithm - ) => Math.Pi.Algorithm (Algorithm squareRootAlgorithm factorialAlgorithm) where + Math.SquareRoot.Algorithmic squareRootAlgorithm, + Math.Factorial.Algorithmic factorialAlgorithm + ) => Math.Pi.Algorithmic (Algorithm squareRootAlgorithm factorialAlgorithm) where openR (Borwein1993 squareRootAlgorithm factorialAlgorithm) = Math.Implementations.Pi.Borwein.Implementation.openR Math.Implementations.Pi.Borwein.Borwein1993.series squareRootAlgorithm factorialAlgorithm 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 @@ -38,7 +38,7 @@ import qualified Factory.Math.SquareRoot as Math.SquareRoot -- | Defines the parameters of the /Borwein/ series. -series :: (Math.SquareRoot.Algorithm squareRootAlgorithm, Math.Factorial.Algorithm factorialAlgorithm) => Math.Implementations.Pi.Borwein.Series.Series squareRootAlgorithm factorialAlgorithm +series :: (Math.SquareRoot.Algorithmic squareRootAlgorithm, Math.Factorial.Algorithmic factorialAlgorithm) => Math.Implementations.Pi.Borwein.Series.Series squareRootAlgorithm factorialAlgorithm 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 diff --git a/src/Factory/Math/Implementations/Pi/Ramanujan/Algorithm.hs b/src/Factory/Math/Implementations/Pi/Ramanujan/Algorithm.hs --- a/src/Factory/Math/Implementations/Pi/Ramanujan/Algorithm.hs +++ b/src/Factory/Math/Implementations/Pi/Ramanujan/Algorithm.hs @@ -47,9 +47,9 @@ defaultValue = Chudnovsky Defaultable.defaultValue Defaultable.defaultValue instance ( - Math.SquareRoot.Algorithm squareRootAlgorithm, - Math.Factorial.Algorithm factorialAlgorithm - ) => Math.Pi.Algorithm (Algorithm squareRootAlgorithm factorialAlgorithm) where + Math.SquareRoot.Algorithmic squareRootAlgorithm, + Math.Factorial.Algorithmic factorialAlgorithm + ) => Math.Pi.Algorithmic (Algorithm squareRootAlgorithm factorialAlgorithm) where openR (Classic squareRootAlgorithm factorialAlgorithm) = Math.Implementations.Pi.Ramanujan.Implementation.openR Math.Implementations.Pi.Ramanujan.Classic.series squareRootAlgorithm factorialAlgorithm openR (Chudnovsky squareRootAlgorithm factorialAlgorithm) = Math.Implementations.Pi.Ramanujan.Implementation.openR Math.Implementations.Pi.Ramanujan.Chudnovsky.series squareRootAlgorithm factorialAlgorithm 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 @@ -37,8 +37,8 @@ -- | Defines the parameters of the /Chudnovsky/ series. series :: ( - Math.SquareRoot.Algorithm squareRootAlgorithm, - Math.Factorial.Algorithm factorialAlgorithm + Math.SquareRoot.Algorithmic squareRootAlgorithm, + Math.Factorial.Algorithmic factorialAlgorithm ) => Math.Implementations.Pi.Ramanujan.Series.Series squareRootAlgorithm factorialAlgorithm series = Math.Implementations.Pi.Ramanujan.Series.MkSeries { Math.Implementations.Pi.Ramanujan.Series.terms = \factorialAlgorithm -> zipWith ( 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 @@ -36,7 +36,7 @@ import qualified Factory.Math.SquareRoot as Math.SquareRoot -- | Defines the parameters of the /Ramanujan/ series. -series :: (Math.SquareRoot.Algorithm squareRootAlgorithm, Math.Factorial.Algorithm factorialAlgorithm) => Math.Implementations.Pi.Ramanujan.Series.Series squareRootAlgorithm factorialAlgorithm +series :: (Math.SquareRoot.Algorithmic squareRootAlgorithm, Math.Factorial.Algorithmic factorialAlgorithm) => Math.Implementations.Pi.Ramanujan.Series.Series squareRootAlgorithm factorialAlgorithm series = Math.Implementations.Pi.Ramanujan.Series.MkSeries { Math.Implementations.Pi.Ramanujan.Series.terms = \factorialAlgorithm -> let toFourthPower = (^ (4 :: Int)) 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 @@ -42,7 +42,7 @@ instance Defaultable.Defaultable Algorithm where defaultValue = Gosper -instance Math.Pi.Algorithm Algorithm where +instance Math.Pi.Algorithmic Algorithm where 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 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 @@ -65,7 +65,7 @@ instance Defaultable.Defaultable (Algorithm factorisationAlgorithm) where defaultValue = MillerRabin -instance Math.PrimeFactorisation.Algorithm factorisationAlgorithm => Math.Primality.Algorithm (Algorithm factorisationAlgorithm) where +instance Math.PrimeFactorisation.Algorithmic factorisationAlgorithm => Math.Primality.Algorithmic (Algorithm factorisationAlgorithm) where isPrime _ 2 = True --The only even prime. isPrime algorithm candidate | candidate < 2 || ( @@ -106,7 +106,7 @@ [@Vibhor Bhatt and G. K. Patra@] , -} -isPrimeByAKS :: (Math.PrimeFactorisation.Algorithm factorisationAlgorithm, Control.DeepSeq.NFData i, Integral i) => factorisationAlgorithm -> i -> Bool +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. Math.Primality.areCoprime n `all` filter (/= n) [2 .. r], --Step 3. 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 @@ -66,7 +66,7 @@ instance Defaultable.Defaultable Algorithm where defaultValue = TrialDivision -instance Math.PrimeFactorisation.Algorithm Algorithm where +instance Math.PrimeFactorisation.Algorithmic Algorithm where primeFactors algorithm = case algorithm of DixonsMethod -> factoriseByDixonsMethod FermatsMethod -> Data.PrimeFactors.reduce . factoriseByFermatsMethod diff --git a/src/Factory/Math/Implementations/Primes.hs b/src/Factory/Math/Implementations/Primes.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Primes.hs @@ -0,0 +1,223 @@ +{- + 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 Data.Sequence((|>)) +import qualified Data.Sequence +import qualified Data.Numbers.Primes +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 $ 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 + +-- | 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 primeMultiplesList -> sieve' $ Control.Arrow.second (\m -> foldr insert (Data.Map.delete candidate m) primeMultiplesList) repository --Re-insert subsequent multiples. + Nothing --Not a square-free composite. + | candidate == smallestPrimeSquare -> sieve' $ (tail' *** insert subsequentPrimeMultiples) repository --Migrate subsequent 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. + + insert :: Ord i => Data.PrimeWheel.PrimeMultiples i -> PrimeMultiplesMap i -> PrimeMultiplesMap i + insert [] = error "Factory.Math.Implementations.Primes.sieveOfEratosthenes.sieve.insert:\tnull list" + insert (key : values) = Data.Map.insertWith (++) key [values] --i.e. key => ([values] ++ [Data.PrimeWheel.PrimeMultiples i]) + +{-# 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; + 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 primeMultiplesList -> sieve' $ Control.Arrow.second (\m -> foldr insert (Data.IntMap.delete candidate m) primeMultiplesList) repository + Nothing + | candidate == smallestPrimeSquare -> sieve' $ (tail' *** insert 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 + + insert :: Data.PrimeWheel.PrimeMultiples Int -> PrimeMultiplesMapInt -> PrimeMultiplesMapInt + insert [] = error "Factory.Math.Implementations.Primes.sieveOfEratosthenesInt.sieve.insert:\tnull list" + insert (key : values) = Data.IntMap.insertWith (++) key [values] + 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 @@ -17,7 +17,7 @@ {- | [@AUTHOR@] Dr. Alistair Ward - [@DESCRIPTION@] Implements 'Math.SquareRoot.Algorithm' by a variety of methods. + [@DESCRIPTION@] Implements 'Math.SquareRoot.Algorithmic' by a variety of methods. [@CAVEAT@] @@ -70,7 +70,7 @@ -> operand -- ^ The value for which to find the /square-root/. -> Math.SquareRoot.Result -instance Math.SquareRoot.Algorithm Algorithm where +instance Math.SquareRoot.Algorithmic Algorithm where squareRootFrom _ _ _ 0 = 0 squareRootFrom _ _ _ 1 = 1 squareRootFrom algorithm estimate@(x, decimalDigits) requiredDecimalDigits y diff --git a/src/Factory/Math/MultiplicativeOrder.hs b/src/Factory/Math/MultiplicativeOrder.hs --- a/src/Factory/Math/MultiplicativeOrder.hs +++ b/src/Factory/Math/MultiplicativeOrder.hs @@ -42,7 +42,7 @@ * . -} -multiplicativeOrder :: (Math.PrimeFactorisation.Algorithm primeFactorisationAlgorithm, Control.DeepSeq.NFData i, Integral i) +multiplicativeOrder :: (Math.PrimeFactorisation.Algorithmic primeFactorisationAlgorithm, Control.DeepSeq.NFData i, Integral i) => primeFactorisationAlgorithm -> i -- ^ Base. -> i -- ^ Modulus. 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 @@ -22,7 +22,7 @@ module Factory.Math.Pi( -- * Type-classes - Algorithm(..), + Algorithmic(..), -- * Types -- ** Data-types Category(..) @@ -40,7 +40,7 @@ * Since representing /Pi/ as either a 'Rational' or promoted to an 'Integer', is inconvenient, an alternative decimal 'String'-representation is provided. -} -class Algorithm algorithm where +class Algorithmic algorithm where openR :: algorithm -> Math.Precision.DecimalDigits -> Data.Ratio.Rational -- ^ Returns the value of /Pi/ as a 'Rational'. openI :: algorithm -> Math.Precision.DecimalDigits -> Integer -- ^ Returns the value of /Pi/, promoted by the required precision to form an integer. @@ -75,12 +75,12 @@ defaultValue = BBP Defaultable.defaultValue instance ( - Algorithm agm, - Algorithm bbp, - Algorithm borwein, - Algorithm ramanujan, - Algorithm spigot - ) => Algorithm (Category agm bbp borwein ramanujan spigot) where + Algorithmic agm, + Algorithmic bbp, + Algorithmic borwein, + Algorithmic ramanujan, + Algorithmic spigot + ) => 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) 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 @@ -30,8 +30,10 @@ raiseModulo, -- ** Predicates isPerfectPower +-- isPerfectPowerInt ) where +import qualified Data.IntSet import qualified Data.Set -- | Mainly for convenience. @@ -49,7 +51,9 @@ * The initial value doesn't need to be either positive or integral. -} -squaresFrom :: Num n => n -> [(n, n)] +squaresFrom :: Num n + => n -- ^ Lower bound. + -> [(n, n)] -- ^ @ [(n, n^2)] @. squaresFrom from = iterate (\(x, y) -> (x + 1, y + 2 * x + 1)) (from, square from) -- | Just for convenience. @@ -107,7 +111,7 @@ (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, whereas 187598574531033121 is square. +-- ] && fromIntegral iSqrt == sqrt' = Just iSqrt --CAVEAT: erroneously True for 187598574531033120 (187598574531033121 is square). ] && square iSqrt == i = Just iSqrt | otherwise = Nothing where @@ -124,6 +128,8 @@ * . * . + + * A generalisation of the concept of /perfect squares/, in which only the exponent '2' is significant. -} isPerfectPower :: Integral i => i -> Bool isPerfectPower i @@ -131,7 +137,20 @@ | otherwise = i `Data.Set.member` foldr ( \n set -> if n `Data.Set.member` set then set --- else Data.Set.union set . Data.Set.fromList . takeWhile (<= i) . iterate (* n) $ square n --TODO: test relative speed. - else foldr Data.Set.insert set . takeWhile (<= i) . iterate (* n) $ square n +-- 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/Primality.hs b/src/Factory/Math/Primality.hs --- a/src/Factory/Math/Primality.hs +++ b/src/Factory/Math/Primality.hs @@ -26,7 +26,7 @@ module Factory.Math.Primality( -- * Type-classes - Algorithm(..), + Algorithmic(..), -- * Functions carmichaelNumbers, -- ** Predicates @@ -38,8 +38,8 @@ import qualified Control.DeepSeq import qualified Factory.Math.Power as Math.Power --- | Defines the methods expected of a primality-algorithm. -class Algorithm algorithm where +-- | Defines the methods expected of a primality-testing algorithm. +class Algorithmic algorithm where isPrime :: (Control.DeepSeq.NFData i, Integral i) => algorithm -> i -> Bool {- | @@ -79,7 +79,7 @@ * . -} -isCarmichaelNumber :: (Algorithm algorithm, Control.DeepSeq.NFData i, Integral i) => algorithm -> i -> Bool +isCarmichaelNumber :: (Algorithmic algorithm, Control.DeepSeq.NFData i, Integral i) => algorithm -> i -> Bool isCarmichaelNumber algorithm i = not $ or [ i <= 2, even i, @@ -88,5 +88,5 @@ ] -- | An ordered list of the /Carmichael/ numbers; . -carmichaelNumbers :: (Algorithm algorithm, Control.DeepSeq.NFData i, Integral i) => algorithm -> [i] +carmichaelNumbers :: (Algorithmic algorithm, Control.DeepSeq.NFData i, Integral i) => algorithm -> [i] carmichaelNumbers algorithm = isCarmichaelNumber algorithm `filter` [3, 5 ..] diff --git a/src/Factory/Math/PrimeFactorisation.hs b/src/Factory/Math/PrimeFactorisation.hs --- a/src/Factory/Math/PrimeFactorisation.hs +++ b/src/Factory/Math/PrimeFactorisation.hs @@ -24,7 +24,7 @@ * Exports a common interface to permit decomposition of positive integers, into the unique combination of /prime/-factors known to exist according to the /Fundamental Theorem of Arithmetic/; . - * Leveraging this abstract capability, it derives the /smoothness/, /power-smoothness/, and /omega/-numbers. + * Leveraging this abstract capability, it derives the /smoothness/, /power-smoothness/, /omega/-numbers and /square-free/ integers. * Filters the list of /regular-numbers/ from the list of /smoothness/. @@ -33,7 +33,7 @@ module Factory.Math.PrimeFactorisation( -- * Type-classes - Algorithm(..), + Algorithmic(..), -- * Functions maxBoundPrimeFactor, smoothness, @@ -41,7 +41,8 @@ regularNumbers, primePowerTotient, eulersTotient, - omega + omega, + squareFree ) where import qualified Control.DeepSeq @@ -50,7 +51,7 @@ import qualified Factory.Data.PrimeFactors as Data.PrimeFactors -- | Defines the methods expected of a /factorisation/-algorithm. -class Algorithm algorithm where +class Algorithmic algorithm where primeFactors :: (Control.DeepSeq.NFData base, Integral base) => algorithm -> base -- ^ The operand @@ -63,8 +64,10 @@ but though a prime-factor /greater/ than the /square-root/ of the number can exist, its smaller /cofactor/ decomposes to a prime which must be less than the /square-root/. - * NB: rather using @(primeFactor <= sqrt numerator)@ to filter the candidate prime-factors of a given numerator, + * NB: rather then using @(primeFactor <= sqrt numerator)@ to filter the candidate prime-factors of a given numerator, one can alternatively use @(numerator >= primeFactor ^ 2)@ to filter what can potentially be factored by a given prime-factor. + + * CAVEAT: suffers from rounding-errors, though no consequence has been witnessed. -} maxBoundPrimeFactor :: Integral i => i -> i maxBoundPrimeFactor = floor . (sqrt :: Double -> Double) . fromIntegral @@ -76,7 +79,7 @@ * . -} -smoothness :: (Algorithm algorithm, Control.DeepSeq.NFData base, Integral base) => algorithm -> [base] +smoothness :: (Algorithmic algorithm, Control.DeepSeq.NFData base, Integral base) => algorithm -> [base] smoothness algorithm = 0 : map (Data.Exponential.getBase . last . primeFactors algorithm) [1 ..] {- | @@ -84,7 +87,7 @@ * . -} -powerSmoothness :: (Algorithm algorithm, Control.DeepSeq.NFData base, Integral base) => algorithm -> [base] +powerSmoothness :: (Algorithmic algorithm, Control.DeepSeq.NFData base, Integral base) => algorithm -> [base] powerSmoothness algorithm = 0 : map (maximum . map Data.Exponential.evaluate . primeFactors algorithm) [1 ..] {- | @@ -92,7 +95,7 @@ * . -} -regularNumbers :: (Algorithm algorithm, Control.DeepSeq.NFData base, Integral base) => algorithm -> [base] +regularNumbers :: (Algorithmic algorithm, Control.DeepSeq.NFData base, Integral base) => algorithm -> [base] regularNumbers algorithm = map fst . filter ((<= (5 :: Integer)) . snd) . zip [1 ..] . tail $ smoothness algorithm {- | @@ -115,14 +118,14 @@ * AKA /EulerPhi/. -} -eulersTotient :: (Algorithm algorithm, Control.DeepSeq.NFData i, Integral i) => algorithm -> i -> i +eulersTotient :: (Algorithmic algorithm, Control.DeepSeq.NFData i, Integral i) => algorithm -> i -> i eulersTotient _ 1 = 1 eulersTotient algorithm i | i <= 0 = error $ "Factory.Math.PrimeFactorisation.eulersTotient:\tundefined for; " ++ show i | otherwise = product . map primePowerTotient $ primeFactors algorithm i {- | - * A constant, zero-indexed, conceptually infinite, list of the /small omega/ numbers, the number of /distinct/ prime factors; cf. /big omega/. + * A constant, zero-indexed, conceptually infinite, list of the /small omega/ numbers (i.e. the number of /distinct/ prime factors); cf. /big omega/. * . @@ -130,6 +133,14 @@ * . -} -omega :: (Algorithm algorithm, Control.DeepSeq.NFData i, Integral i) => algorithm -> [i] +omega :: (Algorithmic algorithm, Control.DeepSeq.NFData i, Integral i) => algorithm -> [i] omega algorithm = map (Data.List.genericLength . primeFactors algorithm) [0 :: Integer ..] + +{- | + * A constant, conceptually infinite, list of the /square-free/ numbers, i.e. those which aren't divisible by any /perfect square/. + + * . +-} +squareFree :: (Algorithmic algorithm, Control.DeepSeq.NFData i, Integral i) => algorithm -> [i] +squareFree algorithm = filter (all (== 1) . map Data.Exponential.getExponent . primeFactors algorithm) [1 ..] diff --git a/src/Factory/Math/Primes.hs b/src/Factory/Math/Primes.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Primes.hs @@ -0,0 +1,42 @@ +{- + 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 a common interface for implementations of /prime-number/ generators. +-} + +module Factory.Math.Primes( +-- * Types-classes + Algorithmic(..), +-- * Functions + primorial +) where + +-- | 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. + +{- | + * Returns the constant, infinite list, defining the /Primorial/. + + * . + + * . +-} +primorial :: (Integral i, Algorithmic algorithm) => 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 @@ -32,8 +32,8 @@ toBase ) where -import Data.Array((!)) -import qualified Data.Array +import Data.Array.IArray((!)) +import qualified Data.Array.IArray import qualified Data.Char import qualified Data.List import qualified Data.Maybe @@ -43,8 +43,8 @@ digits = ['0' .. '9'] ++ ['a' .. 'z'] -- | Constant random-access lookup for 'digits'. -encodes :: (Data.Array.Ix index, Integral index) => Data.Array.Array index Char -encodes = Data.Array.listArray (0, fromIntegral $ length digits - 1) 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 -- | Constant reverse-lookup for 'digits'. decodes :: Integral i => [(Char, i)] @@ -58,14 +58,14 @@ * The conversion to 'Char' can only succeed where printable and intelligible characters exist to represent all digits in the chosen base; which in practice means @(-36 <= base <= 36)@. -} -toBase :: (Integral base, Integral decimal) => base -> decimal -> String +toBase :: (Integral base, Data.Array.IArray.Ix decimal, Integral decimal) => base -> decimal -> String toBase 10 decimal = show decimal --Base unchanged. toBase _ 0 = "0" --Zero has the same representation in any base. toBase base decimal | abs base < 2 = error $ "Factory.Math.Radix.toBase:\tan arbitrary integer can't be represented in base " ++ show base | abs base > Data.List.genericLength digits = error $ "Factory.Math.Radix.toBase:\tunable to clearly represent the complete set of digits in base " ++ show base - | base > 0 && decimal < 0 = '-' : map (toDigit . fromIntegral) (fromDecimal (negate decimal) []) - | otherwise = (toDigit . fromIntegral) `map` fromDecimal decimal [] + | base > 0 && decimal < 0 = '-' : map toDigit (fromDecimal (negate decimal) []) + | otherwise = toDigit `map` fromDecimal decimal [] where fromDecimal 0 = id fromDecimal n @@ -74,14 +74,14 @@ where (quotient, remainder) = n `quotRem` fromIntegral base - toDigit :: Int -> Char + toDigit :: (Data.Array.IArray.Ix i, Integral i) => i -> Char toDigit n | n >&< encodes = encodes ! n | otherwise = error $ "Factory.Math.Radix.toBase.toDigit:\tno suitable character-representation for integer " ++ show n where - (>&<) :: Int -> Data.Array.Array Int Char -> Bool + (>&<) :: (Data.Array.IArray.Ix i, Integral i) => i -> Data.Array.IArray.Array i Char -> Bool index >&< array = ($ index) `all` [(>= lower), (<= upper)] where - (lower, upper) = Data.Array.bounds array + (lower, upper) = Data.Array.IArray.bounds array {- | * Convert the 'String'-representation of a number in the specified base, to a decimal integer. @@ -95,8 +95,8 @@ | abs base < 2 = error $ "Factory.Math.Radix.fromBase:\tan arbitrary integer can't be represented in base " ++ show base | abs base > Data.List.genericLength digits = error $ "Factory.Math.Radix.fromBase:\tunable to clearly represent the complete set of digits in base " ++ show base | base > 0 && head s == '-' = negate . fromBase base $ tail s --Recurse. - | otherwise = Data.List.foldl' (\l -> ((l * fromIntegral base) +) . fromIntegral . fromDigit) 0 s where - fromDigit :: Char -> Int + | otherwise = Data.List.foldl' (\l -> ((l * fromIntegral base) +) . fromDigit) 0 s where + fromDigit :: Integral i => Char -> i fromDigit c = case c `lookup` decodes of Just i | i >= abs (fromIntegral base) -> error $ "Factory.Math.Radix.fromBase.fromDigit:\tillegal char " ++ show c ++ ", for base " ++ show base @@ -108,11 +108,11 @@ * . -} -digitSum :: (Integral base, Integral decimal) => base -> decimal -> decimal +digitSum :: (Integral base, Data.Array.IArray.Ix decimal, Integral decimal) => base -> decimal -> decimal digitSum 10 = fromIntegral . foldr ((+) . Data.Char.digitToInt) 0 . show digitSum base = sum . Data.Maybe.mapMaybe (`lookup` decodes) . toBase base -- | . -digitalRoot :: Integral decimal => decimal -> decimal +digitalRoot :: (Data.Array.IArray.Ix decimal, Integral decimal) => decimal -> decimal digitalRoot = head . dropWhile (> 9) . iterate (digitSum (10 :: Int)) diff --git a/src/Factory/Math/SquareRoot.hs b/src/Factory/Math/SquareRoot.hs --- a/src/Factory/Math/SquareRoot.hs +++ b/src/Factory/Math/SquareRoot.hs @@ -26,7 +26,7 @@ module Factory.Math.SquareRoot( -- * Type-classes - Algorithm(..), + Algorithmic(..), Iterator(..), -- * Types -- ** Type-synonyms @@ -52,7 +52,7 @@ type Estimate = (Result, Math.Precision.DecimalDigits) -- | Defines the methods expected of a /square-root/ algorithm. -class Algorithm algorithm where +class Algorithmic algorithm where squareRootFrom :: Real operand => algorithm -> Estimate -- ^ An initial estimate from which to start. @@ -85,7 +85,7 @@ getEstimate :: Real operand => operand -> Estimate getEstimate y | y < 0 = error $ "Factory.Math.SquareRoot.getEstimate:\tthere's no real square-root of " ++ show y - | otherwise = (Math.Precision.simplify decimalDigits {-doubles performance by roughly halving number's length-} . toRational $ rSqrt y, decimalDigits) + | otherwise = (Math.Precision.simplify decimalDigits {-doubles performance by roughly length of the Rational representation-} . toRational $ rSqrt y, decimalDigits) where decimalDigits :: Math.Precision.DecimalDigits decimalDigits = 16 -- . 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 @@ -17,7 +17,7 @@ {- | [@AUTHOR@] Dr. Alistair Ward - [@DESCRIPTION@] Miscellaneous statistical functions. + [@DESCRIPTION@] Miscellaneous statistics functions. -} module Factory.Math.Statistics( @@ -51,14 +51,14 @@ getMean l = uncurry (/) . (realToFrac *** fromIntegral) $ foldr (\s -> (+ s) *** succ) (0, 0 :: Int) l {- | - * Measures the dispersion of a population of results from the mean value; . + * Measures the /dispersion/ of a /population/ of results from the /mean/ value; . * Should the caller define the result-type as 'Data.Ratio.Rational', then it will be free from rounding-errors. -} getDispersionFromMean :: (Real r, Fractional result) => (Data.Ratio.Rational -> Data.Ratio.Rational) -> [r] -> result -getDispersionFromMean _ [] = error "Factory.Math.Statistics.getDispersionFromMean:\tundefined result for null-list." -getDispersionFromMean _ [_] = 0 --Not necessary, but a shortcut for this special case. -getDispersionFromMean measure l = getMean $ map (measure . (+ negate (getMean l :: Data.Ratio.Rational)) . realToFrac) l +getDispersionFromMean _ [] = error "Factory.Math.Statistics.getDispersionFromMean:\tundefined result for null-list." +getDispersionFromMean _ [_] = 0 --Not necessary, but a shortcut for this special case. +getDispersionFromMean weighting l = getMean $ map (weighting . (+ negate (getMean l :: Data.Ratio.Rational)) . realToFrac) l {- | * Determines the exact /variance/ of the specified list of numbers; . @@ -88,8 +88,8 @@ where mean = getMean l --- | The number of unordered combinations of /r/ objects taken from /n/; . -nCr :: (Math.Factorial.Algorithm factorialAlgorithm, Integral i) +-- | The number of unordered /combinations/ of /r/ objects taken from /n/; . +nCr :: (Math.Factorial.Algorithmic factorialAlgorithm, Integral i) => factorialAlgorithm -> i -- ^ The total number of items from which to select. -> i -- ^ The number of items in a sample. @@ -106,7 +106,7 @@ numerator = Math.Implementations.Factorial.risingFactorial (bigger + 1) (n - bigger) denominator = Math.Factorial.factorial factorialAlgorithm smaller --- | The number of permutations of /r/ objects taken from /n/; . +-- | The number of /permutations/ of /r/ objects taken from /n/; . nPr :: Integral i => i -- ^ The total number of items from which to select. -> i -- ^ The number of items in a sample. diff --git a/src/Factory/Math/Summation.hs b/src/Factory/Math/Summation.hs --- a/src/Factory/Math/Summation.hs +++ b/src/Factory/Math/Summation.hs @@ -69,7 +69,7 @@ sum' _ = sum #endif -{- +{- | * Sums a list of /rational/ type numbers. * CAVEAT: though faster than 'Data.List.sum', this algorithm has poor space-complexity, making it unsuitable for unrestricted use. diff --git a/src/Factory/Test/Performance/Factorial.hs b/src/Factory/Test/Performance/Factorial.hs --- a/src/Factory/Test/Performance/Factorial.hs +++ b/src/Factory/Test/Performance/Factorial.hs @@ -17,7 +17,7 @@ {- | [@AUTHOR@] Dr. Alistair Ward - [@DESCRIPTION@] Times functions exported from module "Math.Factorial". + [@DESCRIPTION@] Times the methods exported from module "Math.Factorial". -} module Factory.Test.Performance.Factorial( @@ -34,7 +34,7 @@ import qualified ToolShed.TimePure as TimePure -- | Measures the CPU-time required by 'Math.Factorial.factorial'. -factorialPerformance :: (Math.Factorial.Algorithm algorithm, Control.DeepSeq.NFData i, Integral i) => algorithm -> i -> IO (Double, i) +factorialPerformance :: (Math.Factorial.Algorithmic algorithm, Control.DeepSeq.NFData i, Integral i) => algorithm -> i -> IO (Double, i) factorialPerformance algorithm = TimePure.getCPUSeconds . Math.Factorial.factorial algorithm -- | Measures the CPU-time required by a naive implementation. @@ -47,7 +47,7 @@ * CAVEAT: nothing is returned, since the result is printed ... and it never terminates. -} -factorialPerformanceGraph :: Math.Factorial.Algorithm algorithm => Bool -> algorithm -> IO () +factorialPerformanceGraph :: Math.Factorial.Algorithmic algorithm => Bool -> algorithm -> IO () factorialPerformanceGraph verbose algorithm = mapM_ ( \operand -> factorialPerformance algorithm operand >>= putStrLn . shows operand . showChar '\t' . ( if verbose diff --git a/src/Factory/Test/Performance/Hyperoperation.hs b/src/Factory/Test/Performance/Hyperoperation.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/Performance/Hyperoperation.hs @@ -0,0 +1,71 @@ +{- + 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@] Times functions exported from module "Math.Hyperoperation". +-} + +module Factory.Test.Performance.Hyperoperation( +-- * Functions + hyperoperationPerformance, + hyperoperationPerformanceGraphRank, + hyperoperationPerformanceGraphExponent +) where + +import qualified Factory.Math.Hyperoperation as Math.Hyperoperation +import qualified ToolShed.TimePure as TimePure + +-- | Measures the CPU-time required by 'Math.Hyperoperation.hyperoperation'. +hyperoperationPerformance :: Integral rank => rank -> Math.Hyperoperation.Base -> Math.Hyperoperation.HyperExponent -> IO (Double, Integer) +hyperoperationPerformance rank base = TimePure.getCPUSeconds . Math.Hyperoperation.hyperoperation rank base + +{- | + * Measure the CPU-time required by 'Math.Hyperoperation.hyperoperation', against a linearly increasing /rank/. + + * CAVEAT: nothing is returned, since the result is printed ... and it never terminates. +-} +hyperoperationPerformanceGraphRank :: + Bool -- ^ Verbose. + -> Math.Hyperoperation.Base + -> Math.Hyperoperation.HyperExponent + -> IO () +hyperoperationPerformanceGraphRank verbose base hyperExponent = mapM_ ( + \rank -> hyperoperationPerformance rank base hyperExponent >>= putStrLn . shows rank . showChar '\t' . ( + if verbose + then (`shows` "") + else (`shows` "") . fst + ) + ) [0 :: Int ..] + +{- | + * Measure the CPU-time required by 'Math.Hyperoperation.hyperoperation', against a linearly increasing /hyper-exponent/. + + * CAVEAT: nothing is returned, since the result is printed ... and it never terminates. +-} +hyperoperationPerformanceGraphExponent :: Integral rank + => Bool -- ^ Verbose. + -> rank + -> Math.Hyperoperation.Base + -> IO () +hyperoperationPerformanceGraphExponent verbose rank base = mapM_ ( + \hyperExponent -> hyperoperationPerformance rank base hyperExponent >>= putStrLn . shows hyperExponent . showChar '\t' . ( + if verbose + then (`shows` "") + else (`shows` "") . fst + ) + ) [0 ..] diff --git a/src/Factory/Test/Performance/Pi.hs b/src/Factory/Test/Performance/Pi.hs --- a/src/Factory/Test/Performance/Pi.hs +++ b/src/Factory/Test/Performance/Pi.hs @@ -17,7 +17,7 @@ {- | [@AUTHOR@] Dr. Alistair Ward - [@DESCRIPTION@] Times functions exported from module "Math.Pi". + [@DESCRIPTION@] Times the methods exported from module "Math.Pi". -} module Factory.Test.Performance.Pi( @@ -51,8 +51,8 @@ -- | Measures the CPU-time required to find Pi to the required precision. piPerformance :: ( - Math.SquareRoot.Algorithm squareRootAlgorithm, - Math.Factorial.Algorithm factorialAlgorithm + Math.SquareRoot.Algorithmic squareRootAlgorithm, + Math.Factorial.Algorithmic factorialAlgorithm ) => Category squareRootAlgorithm factorialAlgorithm -> Math.Precision.DecimalDigits -> IO (Double, String) piPerformance category = TimePure.getCPUSeconds . Math.Pi.openS category @@ -62,9 +62,9 @@ * CAVEAT: nothing is returned, since the result is printed ... and it never terminates. -} piPerformanceGraph :: ( - Math.SquareRoot.Algorithm squareRootAlgorithm, + Math.SquareRoot.Algorithmic squareRootAlgorithm, Show squareRootAlgorithm, - Math.Factorial.Algorithm factorialAlgorithm, + Math.Factorial.Algorithmic factorialAlgorithm, Show factorialAlgorithm ) => RealFrac i => Category squareRootAlgorithm factorialAlgorithm -- ^ The algorithm. diff --git a/src/Factory/Test/Performance/Primality.hs b/src/Factory/Test/Performance/Primality.hs --- a/src/Factory/Test/Performance/Primality.hs +++ b/src/Factory/Test/Performance/Primality.hs @@ -27,18 +27,19 @@ isPrimePerformanceGraph ) where +import qualified Control.DeepSeq import qualified Factory.Math.Fibonacci as Math.Fibonacci import qualified Factory.Math.Primality as Math.Primality import qualified ToolShed.TimePure as TimePure --- | Measures the CPU-time required to find the specified number of /Carmichael/ numbers, which is returned together with the requested list. -carmichaelNumbersPerformance :: Math.Primality.Algorithm primalityAlgorithm => primalityAlgorithm -> Int -> IO (Double, [Integer]) +-- | Measures the CPU-time required to find the specified number of /Carmichael/-numbers, which is returned together with the requested list. +carmichaelNumbersPerformance :: Math.Primality.Algorithmic primalityAlgorithm => primalityAlgorithm -> Int -> IO (Double, [Integer]) carmichaelNumbersPerformance primalityAlgorithm i | i < 0 = error $ "Factory.Test.Performance.Primality.carmichaelNumbersPerformance:\tnegative number; " ++ show i | otherwise = TimePure.getCPUSeconds . take i $ Math.Primality.carmichaelNumbers primalityAlgorithm -- | Measures the CPU-time required to determine whether the specified integer is prime, which is returned together with the Boolean result. -isPrimePerformance :: Math.Primality.Algorithm primalityAlgorithm => primalityAlgorithm -> Integer -> IO (Double, Bool) +isPrimePerformance :: (Control.DeepSeq.NFData i, Integral i) => Math.Primality.Algorithmic primalityAlgorithm => primalityAlgorithm -> i -> IO (Double, Bool) isPrimePerformance primalityAlgorithm = TimePure.getCPUSeconds . Math.Primality.isPrime primalityAlgorithm {- | @@ -46,8 +47,8 @@ * CAVEAT: nothing is returned, since the result is printed ... and it never terminates. -} -isPrimePerformanceGraph :: Math.Primality.Algorithm primalityAlgorithm => primalityAlgorithm -> IO () +isPrimePerformanceGraph :: Math.Primality.Algorithmic primalityAlgorithm => primalityAlgorithm -> IO () isPrimePerformanceGraph primalityAlgorithm = mapM_ ( \operand -> isPrimePerformance primalityAlgorithm operand >>= putStrLn . shows operand . showChar '\t' . (`shows` "") - ) Math.Fibonacci.primeIndexedFibonacci + ) (Math.Fibonacci.primeIndexedFibonacci :: [Integer]) diff --git a/src/Factory/Test/Performance/PrimeFactorisation.hs b/src/Factory/Test/Performance/PrimeFactorisation.hs --- a/src/Factory/Test/Performance/PrimeFactorisation.hs +++ b/src/Factory/Test/Performance/PrimeFactorisation.hs @@ -17,7 +17,7 @@ {- | [@AUTHOR@] Dr. Alistair Ward - [@DESCRIPTION@] Times functions exported by module "Math.PrimeFactorisation". + [@DESCRIPTION@] Times the methods exported by module "Math.PrimeFactorisation". -} module Factory.Test.Performance.PrimeFactorisation( @@ -32,7 +32,7 @@ import qualified ToolShed.TimePure as TimePure -- | Measures the CPU-time required to prime-factorise the specified integer, which is returned together with the resulting list of factors. -primeFactorsPerformance :: Math.PrimeFactorisation.Algorithm algorithm => algorithm -> Integer -> IO (Double, Data.PrimeFactors.Factors Integer Int) +primeFactorsPerformance :: Math.PrimeFactorisation.Algorithmic algorithm => algorithm -> Integer -> IO (Double, Data.PrimeFactors.Factors Integer Int) primeFactorsPerformance algorithm = TimePure.getCPUSeconds . Math.PrimeFactorisation.primeFactors algorithm {- | @@ -41,7 +41,7 @@ * CAVEAT: nothing is returned, since the result is printed ... and it never terminates. -} -primeFactorsPerformanceGraph :: Math.PrimeFactorisation.Algorithm algorithm => algorithm -> Int -> IO () +primeFactorsPerformanceGraph :: Math.PrimeFactorisation.Algorithmic algorithm => algorithm -> Int -> IO () primeFactorsPerformanceGraph algorithm tests | tests < 0 = error $ "Factory.Test.Performance.PrimeFactorisation.primeFactorsPerformanceGraph:\tnegative number; " ++ show tests | otherwise = mapM_ ( diff --git a/src/Factory/Test/Performance/Primes.hs b/src/Factory/Test/Performance/Primes.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/Performance/Primes.hs @@ -0,0 +1,34 @@ +{- + 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@] Measures the CPU-time required by "Math.Primes.primes". +-} + +module Factory.Test.Performance.Primes( +-- * Functions + primesPerformance +) where + +import qualified Control.DeepSeq +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 algorithm = TimePure.getCPUSeconds . (Math.Primes.primes algorithm !!) diff --git a/src/Factory/Test/Performance/SquareRoot.hs b/src/Factory/Test/Performance/SquareRoot.hs --- a/src/Factory/Test/Performance/SquareRoot.hs +++ b/src/Factory/Test/Performance/SquareRoot.hs @@ -17,7 +17,7 @@ {- | [@AUTHOR@] Dr. Alistair Ward - [@DESCRIPTION@] Times functions exported from module "Math.SquareRoot". + [@DESCRIPTION@] Measures the CPU-time required by the methods exported from module "Math.SquareRoot". -} module Factory.Test.Performance.SquareRoot( @@ -32,7 +32,7 @@ import qualified ToolShed.TimePure as TimePure -- | Measures the CPU-time required by 'Math.SquareRoot.squareRootFrom', which is returned together with the approximate rational result. -squareRootPerformance :: (Math.SquareRoot.Algorithm algorithm, Real operand) => algorithm -> operand -> Math.Precision.DecimalDigits -> IO (Double, Math.SquareRoot.Result) +squareRootPerformance :: (Math.SquareRoot.Algorithmic algorithm, Real operand) => algorithm -> operand -> Math.Precision.DecimalDigits -> IO (Double, Math.SquareRoot.Result) squareRootPerformance algorithm operand requiredDecimalDigits = TimePure.getCPUSeconds $ Math.SquareRoot.squareRoot algorithm requiredDecimalDigits operand {- | @@ -42,7 +42,7 @@ * CAVEAT: nothing is returned, since the result is printed ... and it never terminates. -} squareRootPerformanceGraph :: ( - Math.SquareRoot.Algorithm algorithm, + Math.SquareRoot.Algorithmic algorithm, Math.SquareRoot.Iterator algorithm, Show algorithm, Real operand diff --git a/src/Factory/Test/Performance/Statistics.hs b/src/Factory/Test/Performance/Statistics.hs --- a/src/Factory/Test/Performance/Statistics.hs +++ b/src/Factory/Test/Performance/Statistics.hs @@ -17,7 +17,7 @@ {- | [@AUTHOR@] Dr. Alistair Ward - [@DESCRIPTION@] Times functions exported from module "Math.Statistics". + [@DESCRIPTION@] Times the functions exported from module "Math.Statistics". -} module Factory.Test.Performance.Statistics( @@ -31,7 +31,7 @@ import qualified ToolShed.TimePure as TimePure -- | Measures the CPU-time required by 'Math.Statistics.nCr'. -nCrPerformance :: (Math.Factorial.Algorithm factorialAlgorithm, Control.DeepSeq.NFData i, Integral i) +nCrPerformance :: (Math.Factorial.Algorithmic factorialAlgorithm, Control.DeepSeq.NFData i, Integral i) => factorialAlgorithm -> i -- ^ The total number from which to select. -> i -- ^ The number of items in a sample. diff --git a/src/Factory/Test/QuickCheck/Hyperoperation.hs b/src/Factory/Test/QuickCheck/Hyperoperation.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/QuickCheck/Hyperoperation.hs @@ -0,0 +1,75 @@ +{- + 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.Hyperoperation". +-} + +module Factory.Test.QuickCheck.Hyperoperation( +-- * Functions + quickChecks +) where + +import qualified Factory.Math.Hyperoperation as Math.Hyperoperation +import qualified Test.QuickCheck + +type Rank = Int + +-- | Defines invariant properties. +quickChecks :: IO () +quickChecks = + Test.QuickCheck.quickCheck prop_rankCoincides + >> Test.QuickCheck.quickCheck prop_baseCoincides + >> Test.QuickCheck.quickCheck prop_hyperExponentCoincides + >> Test.QuickCheck.quickCheck `mapM_` [prop_succ, prop_addition, prop_multiplication, prop_exponentiation] where + 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) + + prop_baseCoincides :: Rank -> Integer -> Test.QuickCheck.Property + prop_baseCoincides rank base = Test.QuickCheck.label "prop_baseCoincides" $ Math.Hyperoperation.hyperoperation rank' base 1 == base where + rank' :: Rank + rank' = 2 + (rank `mod` 1000) + + prop_hyperExponentCoincides :: Rank -> Integer -> Test.QuickCheck.Property + prop_hyperExponentCoincides rank hyperExponent = Test.QuickCheck.label "prop_hyperExponentCoincides" $ Math.Hyperoperation.hyperoperation rank' 1 hyperExponent' == 1 where + rank' :: Rank + rank' = 3 + (rank `mod` 1000) + + hyperExponent' :: Math.Hyperoperation.HyperExponent + 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 + hyperExponent' :: Math.Hyperoperation.HyperExponent + hyperExponent' = abs hyperExponent + + prop_addition base hyperExponent = Test.QuickCheck.label "prop_addition" $ Math.Hyperoperation.hyperoperation Math.Hyperoperation.addition base hyperExponent' == base + fromIntegral hyperExponent' where + hyperExponent' :: Math.Hyperoperation.HyperExponent + hyperExponent' = abs hyperExponent + + prop_multiplication base hyperExponent = Test.QuickCheck.label "prop_multiplication" $ Math.Hyperoperation.hyperoperation Math.Hyperoperation.multiplication base hyperExponent' == base * fromIntegral hyperExponent' where + hyperExponent' :: Math.Hyperoperation.HyperExponent + hyperExponent' = abs hyperExponent + + prop_exponentiation base hyperExponent = Test.QuickCheck.label "prop_exponentiation" $ Math.Hyperoperation.hyperoperation Math.Hyperoperation.exponentiation base hyperExponent' == base ^ hyperExponent' where + hyperExponent' :: Math.Hyperoperation.HyperExponent + hyperExponent' = abs hyperExponent + + 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 @@ -49,7 +49,7 @@ instance ( Test.QuickCheck.Arbitrary squareRootAlgorithm, - Math.SquareRoot.Algorithm squareRootAlgorithm + Math.SquareRoot.Algorithmic squareRootAlgorithm ) => Test.QuickCheck.Arbitrary (Math.Implementations.Pi.AGM.Algorithm.Algorithm squareRootAlgorithm) where arbitrary = Math.Implementations.Pi.AGM.Algorithm.BrentSalamin <$> Test.QuickCheck.arbitrary #if !(MIN_VERSION_QuickCheck(2,1,0)) @@ -64,9 +64,9 @@ instance ( Test.QuickCheck.Arbitrary squareRootAlgorithm, - Math.SquareRoot.Algorithm squareRootAlgorithm, + Math.SquareRoot.Algorithmic squareRootAlgorithm, Test.QuickCheck.Arbitrary factorialAlgorithm, - Math.Factorial.Algorithm factorialAlgorithm + Math.Factorial.Algorithmic factorialAlgorithm ) => Test.QuickCheck.Arbitrary (Math.Implementations.Pi.Borwein.Algorithm.Algorithm squareRootAlgorithm factorialAlgorithm) where arbitrary = Test.QuickCheck.oneof [ Math.Implementations.Pi.Borwein.Algorithm.Borwein1993 <$> Test.QuickCheck.arbitrary <*> Test.QuickCheck.arbitrary @@ -77,9 +77,9 @@ instance ( Test.QuickCheck.Arbitrary squareRootAlgorithm, - Math.SquareRoot.Algorithm squareRootAlgorithm, + Math.SquareRoot.Algorithmic squareRootAlgorithm, Test.QuickCheck.Arbitrary factorialAlgorithm, - Math.Factorial.Algorithm factorialAlgorithm + Math.Factorial.Algorithmic factorialAlgorithm ) => Test.QuickCheck.Arbitrary (Math.Implementations.Pi.Ramanujan.Algorithm.Algorithm squareRootAlgorithm factorialAlgorithm) where arbitrary = Test.QuickCheck.oneof [ Math.Implementations.Pi.Ramanujan.Algorithm.Classic <$> Test.QuickCheck.arbitrary <*> Test.QuickCheck.arbitrary, 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 @@ -33,21 +33,28 @@ -- | Defines invariant properties. quickChecks :: IO () quickChecks = - Test.QuickCheck.quickCheck prop_maybeSquareNumber + Test.QuickCheck.quickCheck `mapM_` [prop_maybeSquareNumber, prop_rewriteRule] >> Test.QuickCheck.quickCheckWith Test.QuickCheck.stdArgs {Test.QuickCheck.maxSuccess = 10000} prop_notSquare - >> Test.QuickCheck.quickCheck prop_squaresFrom + >> Test.QuickCheck.quickCheck `mapM` [prop_squaresFrom, prop_isPerfectPower] >> Test.QuickCheck.quickCheck prop_raiseModulo where - prop_maybeSquareNumber, prop_notSquare :: Integer -> Test.QuickCheck.Property + 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 :: Integer -> Integer -> Test.QuickCheck.Property + prop_squaresFrom, prop_isPerfectPower :: 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 e' :: Integer e' = abs e + diff --git a/src/Factory/Test/QuickCheck/Primes.hs b/src/Factory/Test/QuickCheck/Primes.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/QuickCheck/Primes.hs @@ -0,0 +1,74 @@ +{-# LANGUAGE CPP #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +{- + 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@] Implements 'Test.QuickCheck.Arbitrary' and defines /QuickCheck/-properties for "Math.Implementations.Primes". +-} + +module Factory.Test.QuickCheck.Primes( +-- * Functions + quickChecks +-- isPrime +) where + +import Control.Applicative((<$>)) +import qualified Control.DeepSeq +import qualified Data.Set +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.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 + 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 + ] +#if !(MIN_VERSION_QuickCheck(2,1,0)) + coarbitrary = undefined --CAVEAT: stops warnings from ghc. +#endif + +isPrime :: (Control.DeepSeq.NFData i, Integral i) => i -> Bool +isPrime = Math.Primality.isPrime primalityAlgorithm where + primalityAlgorithm :: Math.Implementations.Primality.Algorithm Math.Implementations.PrimeFactorisation.Algorithm + primalityAlgorithm = Defaultable.defaultValue + +-- | Defines invariant properties. +quickChecks :: IO () +quickChecks = + Test.QuickCheck.quickCheckWith Test.QuickCheck.stdArgs {Test.QuickCheck.maxSuccess = 32} `mapM_` [prop_isPrime, prop_isComposite] + >> Test.QuickCheck.quickCheckWith Test.QuickCheck.stdArgs {Test.QuickCheck.maxSuccess = 32} 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 . take (i `mod` 4096) $ (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] + ) . Data.Set.fromList . takeWhile (<= upperBound) $ Math.Primes.primes algorithm where + upperBound :: Int + upperBound = i `mod` 32768 + + prop_consistency :: Math.Implementations.Primes.Algorithm -> Math.Implementations.Primes.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]) + 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 @@ -64,6 +64,6 @@ Math.Probability.generateDiscretePopulation 1000 (Math.Probability.PoissonDistribution lambda') randomGen :: [Int] ) where lambda' :: Double - lambda' = fromIntegral $ lambda `mod` 1000 + lambda' = fromIntegral $ mod lambda 1000 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 @@ -27,6 +27,7 @@ import qualified Factory.Test.QuickCheck.ArithmeticGeometricMean import qualified Factory.Test.QuickCheck.Factorial +import qualified Factory.Test.QuickCheck.Hyperoperation import qualified Factory.Test.QuickCheck.Interval import qualified Factory.Test.QuickCheck.MonicPolynomial import qualified Factory.Test.QuickCheck.Pi @@ -34,6 +35,7 @@ import qualified Factory.Test.QuickCheck.Power import qualified Factory.Test.QuickCheck.Primality import qualified Factory.Test.QuickCheck.PrimeFactorisation +import qualified Factory.Test.QuickCheck.Primes import qualified Factory.Test.QuickCheck.Probability import qualified Factory.Test.QuickCheck.Radix import qualified Factory.Test.QuickCheck.SquareRoot @@ -44,6 +46,7 @@ run :: IO () run = putStrLn "ArithmeticGeometricMean" >> Factory.Test.QuickCheck.ArithmeticGeometricMean.quickChecks >> putStrLn "Factorial" >> Factory.Test.QuickCheck.Factorial.quickChecks + >> putStrLn "Hyperoperation" >> Factory.Test.QuickCheck.Hyperoperation.quickChecks >> putStrLn "Interval" >> Factory.Test.QuickCheck.Interval.quickChecks >> putStrLn "MonicPolynomial" >> Factory.Test.QuickCheck.MonicPolynomial.quickChecks >> putStrLn "Pi" >> Factory.Test.QuickCheck.Pi.quickChecks @@ -51,6 +54,7 @@ >> putStrLn "Power" >> Factory.Test.QuickCheck.Power.quickChecks >> putStrLn "Primality" >> Factory.Test.QuickCheck.Primality.quickChecks >> putStrLn "PrimeFactorisation" >> Factory.Test.QuickCheck.PrimeFactorisation.quickChecks + >> putStrLn "Primes" >> Factory.Test.QuickCheck.Primes.quickChecks >> putStrLn "Probability" >> Factory.Test.QuickCheck.Probability.quickChecks >> putStrLn "Radix" >> Factory.Test.QuickCheck.Radix.quickChecks >> putStrLn "SquareRoot" >> Factory.Test.QuickCheck.SquareRoot.quickChecks diff --git a/src/Main.hs b/src/Main.hs --- a/src/Main.hs +++ b/src/Main.hs @@ -36,15 +36,19 @@ import qualified Distribution.Package import qualified Distribution.Text import qualified Distribution.Version +import qualified Factory.Math.Hyperoperation as Math.Hyperoperation 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.SquareRoot as Math.Implementations.SquareRoot import qualified Factory.Test.CommandOptions as Test.CommandOptions import qualified Factory.Test.Performance.Factorial as Test.Performance.Factorial +import qualified Factory.Test.Performance.Hyperoperation as Test.Performance.Hyperoperation import qualified Factory.Test.Performance.Pi as Test.Performance.Pi import qualified Factory.Test.Performance.Primality as Test.Performance.Primality import qualified Factory.Test.Performance.PrimeFactorisation as Test.Performance.PrimeFactorisation +import qualified Factory.Test.Performance.Primes as Test.Performance.Primes 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 @@ -79,14 +83,18 @@ G.Option "" ["factorialPerformance"] (factorialPerformance `G.ReqArg` "(Math.Implementations.Factorial.Algorithm, Integer)") "Test the performance of 'Math.Factorial.factorial'.", G.Option "" ["factorialPerformanceGraph"] (factorialPerformanceGraph `G.ReqArg` "Math.Implementations.Factorial.Algorithm") "Test the performance of 'Math.Factorial.factorial', with an exponentially increasing operand.", G.Option "" ["factorialPerformanceGraphControl"] (G.NoArg factorialPerformanceGraphControl) "Test the performance of a naive factorial-implementation, with an exponentially increasing operand.", + G.Option "" ["hyperoperationPerformance"] (hyperoperationPerformance `G.ReqArg` "(Integer, Math.Hyperoperation.Base, Math.Hyperoperation.HyperExponent)") "Test the performance of 'Math.Hyperoperation.hyperoperation', against the specified rank, base and hyper-exponent.", + G.Option "" ["hyperoperationPerformanceGraphRank"] (hyperoperationPerformanceGraphRank `G.ReqArg` "(Math.Hyperoperation.Base, Math.Hyperoperation.HyperExponent)") "Test the performance of 'Math.Hyperoperation.hyperoperation', for the specified base and hyper-exponent, and a linearly increasing rank.", + G.Option "" ["hyperoperationPerformanceGraphExponent"] (hyperoperationPerformanceGraphExponent `G.ReqArg` "(Integer, Math.Hyperoperation.Base)") "Test the performance of 'Math.Hyperoperation.hyperoperation', for the specified rank and base, and a linearly increasing hyper-exponent.", G.Option "" ["isPrimePerformance"] (isPrimePerformance `G.ReqArg` "(Math.Implementations.Primality.Algorithm, Integer)") "Test the performance of 'Math.Primality.isPrime'.", G.Option "" ["isPrimePerformanceGraph"] (isPrimePerformanceGraph `G.ReqArg` "Math.Implementations.Primality.Algorithm") "Test the performance of 'Math.Primality.isPrime', against the prime-indexed Fibonacci-numbers.", G.Option "" ["nCrPerformance"] (nCrPerformance `G.ReqArg` "(Math.Implementations.Factorial.Algorithm, Integer, Integer)") "Test the performance of 'Math.Factorial.factorial'.", G.Option "" ["piPerformance"] (piPerformance `G.ReqArg` "(Math.Pi.Category, Math.Precision.DecimalDigits)") "Test the performance of 'Math.Pi.openI'.", 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 odd integers from the Fibonacci sequence.", - G.Option "" ["squareRootPerformance"] (squareRootPerformance `G.ReqArg` "(Math.Implementations.SquareRoot.Algorithm, Data.Ratio.Rational, DecimalDigits)") "Test 'Math.SquareRoot.squareRoot'.", + 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 "" ["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) ++ "'."), G.Option "" ["version"] (G.NoArg $ const printVersion) "Print version-information & then exit.", @@ -98,7 +106,7 @@ packageIdentifier :: Distribution.Package.PackageIdentifier packageIdentifier = Distribution.Package.PackageIdentifier { Distribution.Package.pkgName = Distribution.Package.PackageName "factory", - Distribution.Package.pkgVersion = Distribution.Version.Version [0, 0, 0, 3] [] + Distribution.Package.pkgVersion = Distribution.Version.Version [0, 1, 0, 2] [] } printUsage = System.IO.hPutStrLn System.IO.stderr usage >> System.exitWith System.ExitSuccess @@ -107,7 +115,7 @@ factorialPerformanceGraphControl :: Test.CommandOptions.CommandOptions -> IO Test.CommandOptions.CommandOptions factorialPerformanceGraphControl commandOptions = Test.Performance.Factorial.factorialPerformanceGraphControl (Test.CommandOptions.verbose commandOptions) >> System.exitWith (System.ExitFailure 1) - carmichaelNumbersPerformance, factorialPerformance, factorialPerformanceGraph, isPrimePerformance, isPrimePerformanceGraph, piPerformance, piPerformanceGraph, primeFactorsPerformance, primeFactorsPerformanceGraph, squareRootPerformance, squareRootPerformanceGraph :: String -> CommandLineAction + carmichaelNumbersPerformance, factorialPerformance, factorialPerformanceGraph, hyperoperationPerformance, hyperoperationPerformanceGraphRank, hyperoperationPerformanceGraphExponent, isPrimePerformance, isPrimePerformanceGraph, piPerformance, piPerformanceGraph, primeFactorsPerformance, primesPerformance, squareRootPerformance, squareRootPerformanceGraph :: String -> CommandLineAction carmichaelNumbersPerformance arg _ = Test.Performance.Primality.carmichaelNumbersPerformance algorithm i >>= print >> System.exitWith System.ExitSuccess where algorithm :: PrimalityAlgorithm @@ -120,8 +128,25 @@ factorialPerformanceGraph arg commandOptions = Test.Performance.Factorial.factorialPerformanceGraph (Test.CommandOptions.verbose commandOptions) (read arg :: Math.Implementations.Factorial.Algorithm) >> System.exitWith (System.ExitFailure 1) + hyperoperationPerformance arg _ = Test.Performance.Hyperoperation.hyperoperationPerformance rank base hyperExponent >>= print >> System.exitWith System.ExitSuccess where + rank :: Integer + base :: Math.Hyperoperation.Base + hyperExponent :: Math.Hyperoperation.HyperExponent + (rank, base, hyperExponent) = read arg + + hyperoperationPerformanceGraphRank arg commandOptions = Test.Performance.Hyperoperation.hyperoperationPerformanceGraphRank (Test.CommandOptions.verbose commandOptions) base hyperExponent >> System.exitWith (System.ExitFailure 1) where + base :: Math.Hyperoperation.Base + hyperExponent :: Math.Hyperoperation.HyperExponent + (base, hyperExponent) = read arg + + hyperoperationPerformanceGraphExponent arg commandOptions = Test.Performance.Hyperoperation.hyperoperationPerformanceGraphExponent (Test.CommandOptions.verbose commandOptions) rank base >> System.exitWith (System.ExitFailure 1) where + rank :: Integer + base :: Math.Hyperoperation.Base + (rank, base) = read arg + isPrimePerformance arg _ = Test.Performance.Primality.isPrimePerformance algorithm i >>= print >> System.exitWith System.ExitSuccess where - algorithm :: PrimalityAlgorithm + algorithm :: PrimalityAlgorithm + i :: Integer (algorithm, i) = read arg isPrimePerformanceGraph arg _ = Test.Performance.Primality.isPrimePerformanceGraph (read arg :: Math.Implementations.Primality.Algorithm Math.Implementations.PrimeFactorisation.Algorithm) >> System.exitWith (System.ExitFailure 1) @@ -144,9 +169,28 @@ algorithm :: Math.Implementations.PrimeFactorisation.Algorithm (algorithm, i) = read arg - primeFactorsPerformanceGraph arg _ = Test.Performance.PrimeFactorisation.primeFactorsPerformanceGraph algorithm i >> System.exitWith (System.ExitFailure 1) where + primeFactorsPerformanceGraph arg _ = Test.Performance.PrimeFactorisation.primeFactorsPerformanceGraph algorithm index >> System.exitWith (System.ExitFailure 1) where algorithm :: Math.Implementations.PrimeFactorisation.Algorithm - (algorithm, i) = read arg + (algorithm, index) = read arg + + primesPerformance arg _ = ( + ( +{- + Hard-code specific algorithms, so the simplifier triggers rewrite-rules in "Math.Implementations.Primes", + ready for run-time definitions of 'algorithm' to exploit as appropriate. + CAVEAT: fragile. +-} + case algorithm of + Math.Implementations.Primes.SieveOfEratosthenes n -> Test.Performance.Primes.primesPerformance $ Math.Implementations.Primes.SieveOfEratosthenes n + _ -> Test.Performance.Primes.primesPerformance algorithm + ) index :: IO ( + Double, +-- Integer + Int --Exploits rewrite-rule in "Math.Implementations.Primes". + ) + ) >>= print >> System.exitWith System.ExitSuccess where + algorithm :: Math.Implementations.Primes.Algorithm + (algorithm, index) = read arg squareRootPerformance arg _ = Test.Performance.SquareRoot.squareRootPerformance algorithm operand decimalDigits >>= print >> System.exitWith System.ExitSuccess where algorithm :: Math.Implementations.SquareRoot.Algorithm