diff --git a/changelog b/changelog --- a/changelog +++ b/changelog @@ -62,3 +62,18 @@ * Added 'Factory.Math.Primes.mersenneNumbers'. * Replaced use of 'mod' on positive integers, with the faster 'rem', in 'Factory.Math.Implementations.Pi.Spigot.Spigot.processColumns', 'Factory.Math.Implementations.Primality.witnessesCompositeness', 'Factory.Math.Implementations.Primes.TrialDivision.isIndivisibleBy', 'Factory.Math.Implementations.Primes.SieveOfAtkin.polynomialTypeLookup', 'Factory.Math.Implementations.Primes.SieveOfAtkin.findPolynomialSolutions', 'Factory.Math.Implementations.Primes.TurnersSieve.turnersSieve', 'Factory.Math.PerfectPower.maybeSquareNumber'. * Replaced calls to 'realToFrac' with 'toRational' in; "Factory.Math.Implementations.SquareRoot", 'Factory.Math.Statistics.getDispersionFromMean', 'Factory.Math.SquareRoot.getDiscrepancy', 'Factory.Math.SquareRoot.getAccuracy', to more clearly represent the required operation. +0.2.1.0 + * Refactored 'Factory.Test.QuickCheck.QuickChecks'. + * Remove redundant import of 'Data.Ratio' from many modules. + * Refactored 'Factory.Math.Radix.encodes' to make use of 'Data.List.genericLength', & removed empty 'where'. + * Explicitly closed standard-input in the executable. + * Replaced calls to 'error' from inside the IO-monad, with 'Control.Monad.fail'. + * Added function 'Factory.Math.Precision.roundTo'. + * Trapped command-line arguments to which garbage has been appended. + * Corrected the output of 'Main.main.optDescrList.printVersion'. + * Removed the integral population-size parameter from 'Factory.Math.Probability.generateContinuousPopulation' & 'Factory.Math.Probability.generateDiscretePopulation', making the result conceptually infinite. + * Created class 'Factory.Math.Probability.Distribution', to which data-types 'Factory.Math.Probability.ContinuousDistribution' & 'Factory.Math.Probability.DiscreteDistribution' conform. + * Added data-constructors 'Factory.Math.Probability.ExponentialDistribution', 'Factory.Math.Probability.ShiftedGeometricDistribution' & 'Factory.Math.Probability.LogNormal'. + * Added command-line option '--plotDiscreteDistribution' to "Main". + * Removed Preprocessor-check on the version of package 'toolshed', in "Factory/Math/Summation" & "Factory/Data/PrimeFactors". + diff --git a/copyright b/copyright --- a/copyright +++ b/copyright @@ -2,7 +2,7 @@ Dr. Alistair Ward . Copyright: - Copyright (C) 2011 Dr. Alistair Ward. All Rights Reserved. + Copyright (C) 2011-2013 Dr. Alistair Ward. All Rights Reserved. Home-page: http://functionalley.eu diff --git a/factory.cabal b/factory.cabal --- a/factory.cabal +++ b/factory.cabal @@ -1,8 +1,8 @@ --Package-properties Name: factory -Version: 0.2.0.5 +Version: 0.2.1.0 Cabal-Version: >= 1.6 -Copyright: (C) 2011 Dr. Alistair Ward +Copyright: (C) 2011-2013 Dr. Alistair Ward License: GPL License-file: LICENSE Author: Dr. Alistair Ward @@ -11,7 +11,7 @@ Build-Type: Simple Description: A library of number-theory functions, for; factorials, square-roots, Pi and primes. Category: Math, Number Theory -Tested-With: GHC == 6.10, GHC == 6.12, GHC == 7.0, GHC == 7.4 +Tested-With: GHC == 7.4 Homepage: http://functionalley.eu Maintainer: factory functionalley eu Bug-reports: factory functionalley eu diff --git a/src/Factory/Data/PrimeFactors.hs b/src/Factory/Data/PrimeFactors.hs --- a/src/Factory/Data/PrimeFactors.hs +++ b/src/Factory/Data/PrimeFactors.hs @@ -1,4 +1,3 @@ -{-# LANGUAGE CPP #-} {- Copyright (C) 2011 Dr. Alistair Ward @@ -49,10 +48,7 @@ import qualified Factory.Math.DivideAndConquer as Math.DivideAndConquer import qualified Factory.Data.Exponential as Data.Exponential import Factory.Data.Exponential((<^), (=~)) - -#if MIN_VERSION_toolshed(11,1,1) import qualified ToolShed.Data.List -#endif infixl 7 >/<, >*< --Same as (/). infixr 8 >^ --Same as (^). @@ -104,12 +100,7 @@ * Preserves the sort-order. -} (>*<) :: (Ord base, Num exponent, Ord exponent) => Factors base exponent -> Factors base exponent -> Factors base exponent -l >*< r = -#if MIN_VERSION_toolshed(11,1,1) - reduceSorted $ ToolShed.Data.List.merge l r -#else - reduce $ l ++ r --CAVEAT: concatenation disorders the list, necessitating a re-sort. -#endif +l >*< r = reduceSorted $ ToolShed.Data.List.merge l r -- | Invert the product of a list /prime factors/, by negating each of the /exponents/. invert :: Num exponent => Factors base exponent -> Factors base exponent 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 @@ -38,7 +38,6 @@ ) where import Control.Arrow((&&&)) -import qualified Data.Ratio import qualified Factory.Math.Precision as Math.Precision import qualified Factory.Math.SquareRoot as Math.SquareRoot @@ -47,10 +46,10 @@ #endif -- | The type of the /arithmetic mean/; . -type ArithmeticMean = Data.Ratio.Rational +type ArithmeticMean = Rational -- | The type of the /geometric mean/; . -type GeometricMean = Data.Ratio.Rational +type GeometricMean = Rational -- | Encapsulates both /arithmetic/ and /geometric/ means. type AGM = (ArithmeticMean, GeometricMean) @@ -72,7 +71,7 @@ | not $ isValid agm = error $ "Factory.Math.ArithmeticGeometricMean.convergeToAGM:\tboth means must be positive for a real geometric mean; " ++ show agm | spread agm == 0 = repeat agm | otherwise = let - simplify :: Data.Ratio.Rational -> Data.Ratio.Rational + simplify :: Rational -> Rational simplify = Math.Precision.simplify (pred decimalDigits {-ignore single integral digit-}) --This makes a gigantic difference to performance. findArithmeticMean :: AGM -> ArithmeticMean @@ -90,7 +89,7 @@ ) agm -- | Returns the bounds within which the 'AGM' has been constrained. -spread :: AGM -> Data.Ratio.Rational +spread :: AGM -> Rational spread = uncurry (-) -- | Checks that both /means/ are positive, as required for the /geometric mean/ to be consistently /real/. diff --git a/src/Factory/Math/DivideAndConquer.hs b/src/Factory/Math/DivideAndConquer.hs --- a/src/Factory/Math/DivideAndConquer.hs +++ b/src/Factory/Math/DivideAndConquer.hs @@ -1,6 +1,6 @@ {-# LANGUAGE CPP #-} {- - Copyright (C) 2010 Dr. Alistair Ward + 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 @@ -118,7 +118,7 @@ * Since the result can be large, 'divideAndConquer' is used in an attempt to form operands of a similar order of magnitude, which creates scope for the use of more efficient multiplication-algorithms. - /Multiplication/ is required for the /addition/ of 'Data.Ratio.Rational' numbers by cross-multiplication; + /Multiplication/ is required for the /addition/ of 'Rational' numbers by cross-multiplication; this function is unlikely to be useful for other numbers. -} sum' :: Num n 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 @@ -35,7 +35,6 @@ ) where import Control.Arrow((&&&)) -import qualified Data.Ratio import qualified Factory.Math.ArithmeticGeometricMean as Math.ArithmeticGeometricMean import qualified Factory.Math.Power as Math.Power import qualified Factory.Math.Precision as Math.Precision @@ -56,10 +55,10 @@ > => 4*a[N]^2 / (1 - sum [2^(n+1) * (a[n-1]^2 - g[n-1]^2)]) -} -openR :: Math.SquareRoot.Algorithmic squareRootAlgorithm => squareRootAlgorithm -> Math.Precision.DecimalDigits -> Data.Ratio.Rational +openR :: Math.SquareRoot.Algorithmic squareRootAlgorithm => squareRootAlgorithm -> Math.Precision.DecimalDigits -> 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 ( Math.Precision.getIterationsRequired Math.Precision.quadraticConvergence 1 decimalDigits - ) $ Math.ArithmeticGeometricMean.convergeToAGM squareRootAlgorithm decimalDigits (1, Math.SquareRoot.squareRoot squareRootAlgorithm decimalDigits (recip 2 :: Data.Ratio.Rational)) + ) $ Math.ArithmeticGeometricMean.convergeToAGM squareRootAlgorithm decimalDigits (1, Math.SquareRoot.squareRoot squareRootAlgorithm decimalDigits (recip 2 :: Rational)) diff --git a/src/Factory/Math/Implementations/Pi/BBP/Implementation.hs b/src/Factory/Math/Implementations/Pi/BBP/Implementation.hs --- a/src/Factory/Math/Implementations/Pi/BBP/Implementation.hs +++ b/src/Factory/Math/Implementations/Pi/BBP/Implementation.hs @@ -21,7 +21,7 @@ * Implements a /Bailey-Borwein-Plouffe/ formula; - * Surprisingly, because of the huge size of the 'Data.Ratio.Rational' quantities, + * Surprisingly, because of the huge size of the 'Rational' quantities, it is a /single/ call to @Factory.Math.Summation.sum'@, rather than the calculation of the many terms in the series, which is the performance-bottleneck. -} @@ -31,7 +31,6 @@ ) where import Data.Ratio((%)) -import qualified Data.Ratio import qualified Factory.Math.Implementations.Pi.BBP.Series as Math.Implementations.Pi.BBP.Series import qualified Factory.Math.Precision as Math.Precision import qualified Factory.Math.Summation as Math.Summation @@ -40,7 +39,7 @@ openR :: Math.Implementations.Pi.BBP.Series.Series -- ^ This /Pi/-algorithm is parameterised by the type of other algorithms to use. -> Math.Precision.DecimalDigits -- ^ The number of decimal digits required. - -> Data.Ratio.Rational + -> Rational openR Math.Implementations.Pi.BBP.Series.MkSeries { Math.Implementations.Pi.BBP.Series.numerators = numerators, Math.Implementations.Pi.BBP.Series.getDenominators = getDenominators, diff --git a/src/Factory/Math/Implementations/Pi/BBP/Series.hs b/src/Factory/Math/Implementations/Pi/BBP/Series.hs --- a/src/Factory/Math/Implementations/Pi/BBP/Series.hs +++ b/src/Factory/Math/Implementations/Pi/BBP/Series.hs @@ -26,13 +26,11 @@ Series(..) ) where -import qualified Data.Ratio - -- | Defines a series corresponding to a specific /BBP/-formula. data Series = MkSeries { numerators :: [Integer], -- ^ The constant numerators from which each term in the series is composed. getDenominators :: Int -> [Integer], -- ^ Generates the term-dependent denominators from which each term in the series is composed. - seriesScalingFactor :: Data.Ratio.Rational, -- ^ The ratio by which the sum to infinity of the series, must be scaled to result in /Pi/. + seriesScalingFactor :: Rational, -- ^ The ratio by which the sum to infinity of the series, must be scaled to result in /Pi/. base :: Integer -- ^ The geometric ratio, by which successive terms are scaled. } 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 @@ -26,7 +26,6 @@ ) where --import Control.Arrow((***)) -import qualified Data.Ratio import Data.Ratio((%)) --import Factory.Data.PrimeFactors((>*<), (>/<), (>^)) --import qualified Factory.Data.PrimeFactors as Data.PrimeFactors @@ -41,11 +40,11 @@ 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 + simplify, squareRoot :: Rational -> Rational simplify = Math.Precision.simplify $ pred decimalDigits {-ignore single integral digit-} --This makes a gigantic difference to performance. squareRoot = simplify . Math.SquareRoot.squareRoot squareRootAlgorithm decimalDigits - sqrt5, a, b, c3 :: Data.Ratio.Rational + sqrt5, a, b, c3 :: Rational sqrt5 = squareRoot 5 a = 63365028312971999585426220 + sqrt5 * (28337702140800842046825600 + 384 * squareRoot (10891728551171178200467436212395209160385656017 + 4870929086578810225077338534541688721351255040 * sqrt5)) diff --git a/src/Factory/Math/Implementations/Pi/Borwein/Implementation.hs b/src/Factory/Math/Implementations/Pi/Borwein/Implementation.hs --- a/src/Factory/Math/Implementations/Pi/Borwein/Implementation.hs +++ b/src/Factory/Math/Implementations/Pi/Borwein/Implementation.hs @@ -27,7 +27,6 @@ ) where import qualified Control.Arrow -import qualified Data.Ratio import qualified Factory.Math.Implementations.Pi.Borwein.Series as Math.Implementations.Pi.Borwein.Series import qualified Factory.Math.Precision as Math.Precision @@ -41,7 +40,7 @@ -> squareRootAlgorithm -- ^ The specific /square-root/ algorithm to apply to the above series. -> factorialAlgorithm -- ^ The specific /factorial/-algorithm to apply to the above series. -> Math.Precision.DecimalDigits -- ^ The number of decimal digits required. - -> Data.Ratio.Rational + -> Rational openR Math.Implementations.Pi.Borwein.Series.MkSeries { Math.Implementations.Pi.Borwein.Series.terms = terms, Math.Implementations.Pi.Borwein.Series.convergenceRate = convergenceRate diff --git a/src/Factory/Math/Implementations/Pi/Borwein/Series.hs b/src/Factory/Math/Implementations/Pi/Borwein/Series.hs --- a/src/Factory/Math/Implementations/Pi/Borwein/Series.hs +++ b/src/Factory/Math/Implementations/Pi/Borwein/Series.hs @@ -26,7 +26,6 @@ Series(..) ) where -import qualified Data.Ratio import qualified Factory.Math.Precision as Math.Precision -- | Defines a series corresponding to a specific /Borwein/-formula. @@ -36,8 +35,8 @@ -> factorialAlgorithm -> Math.Precision.DecimalDigits -> ( - Data.Ratio.Rational, --The factor into which the sum to infinity of the sequence, must be divided to result in /Pi/ - [Data.Ratio.Rational] --The sequence of terms, the sum to infinity of which defines the series. + Rational, --The factor into which the sum to infinity of the sequence, must be divided to result in /Pi/ + [Rational] --The sequence of terms, the sum to infinity of which defines the series. ), convergenceRate :: Math.Precision.ConvergenceRate -- ^ The expected number of digits of /Pi/, per term in the series. } diff --git a/src/Factory/Math/Implementations/Pi/Ramanujan/Implementation.hs b/src/Factory/Math/Implementations/Pi/Ramanujan/Implementation.hs --- a/src/Factory/Math/Implementations/Pi/Ramanujan/Implementation.hs +++ b/src/Factory/Math/Implementations/Pi/Ramanujan/Implementation.hs @@ -26,7 +26,6 @@ openR ) where -import qualified Data.Ratio import qualified Factory.Math.Implementations.Pi.Ramanujan.Series as Math.Implementations.Pi.Ramanujan.Series import qualified Factory.Math.Precision as Math.Precision import qualified Factory.Math.Summation as Math.Summation @@ -41,7 +40,7 @@ -> squareRootAlgorithm -- ^ The specific /square-root/ algorithm to apply to the above series. -> factorialAlgorithm -- ^ The specific /factorial/-algorithm to apply to the above series. -> Math.Precision.DecimalDigits -- ^ The number of decimal digits required. - -> Data.Ratio.Rational + -> Rational openR Math.Implementations.Pi.Ramanujan.Series.MkSeries { Math.Implementations.Pi.Ramanujan.Series.terms = terms, Math.Implementations.Pi.Ramanujan.Series.getSeriesScalingFactor = getSeriesScalingFactor, diff --git a/src/Factory/Math/Implementations/Pi/Ramanujan/Series.hs b/src/Factory/Math/Implementations/Pi/Ramanujan/Series.hs --- a/src/Factory/Math/Implementations/Pi/Ramanujan/Series.hs +++ b/src/Factory/Math/Implementations/Pi/Ramanujan/Series.hs @@ -26,13 +26,12 @@ Series(..) ) where -import qualified Data.Ratio import qualified Factory.Math.Precision as Math.Precision -- | Defines a series corresponding to a specific /Ramanujan/-formula. data Series squareRootAlgorithm factorialAlgorithm = MkSeries { - terms :: factorialAlgorithm -> [Data.Ratio.Rational], -- ^ The sequence of terms, the sum to infinity of which defines the series. - getSeriesScalingFactor :: squareRootAlgorithm -> Math.Precision.DecimalDigits -> Data.Ratio.Rational, -- ^ The ratio by which the sum to infinity of the sequence, must be scaled to result in /Pi/. - convergenceRate :: Math.Precision.ConvergenceRate -- ^ The expected number of digits of /Pi/, per term in the series. + terms :: factorialAlgorithm -> [Rational], -- ^ The sequence of terms, the sum to infinity of which defines the series. + getSeriesScalingFactor :: squareRootAlgorithm -> Math.Precision.DecimalDigits -> Rational, -- ^ The ratio by which the sum to infinity of the sequence, must be scaled to result in /Pi/. + convergenceRate :: Math.Precision.ConvergenceRate -- ^ The expected number of digits of /Pi/, per term in the series. } diff --git a/src/Factory/Math/Implementations/Primes/SieveOfEratosthenes.hs b/src/Factory/Math/Implementations/Primes/SieveOfEratosthenes.hs --- a/src/Factory/Math/Implementations/Primes/SieveOfEratosthenes.hs +++ b/src/Factory/Math/Implementations/Primes/SieveOfEratosthenes.hs @@ -57,7 +57,12 @@ head' :: Data.Sequence.Seq [a] -> [a] head' = (`Data.Sequence.index` 0) --- | The 'Data.Sequence.Seq' counterpart to 'Data.List.tail'. +{- | + * The 'Data.Sequence.Seq' counterpart to 'Data.List.tail'. + + * CAVEAT: because @ Data.List.tail [] @ returns an error, whereas @ tail' Data.Sequence.empty @ returns 'Data.Sequence.empty', + this function is for internal use only. +-} tail' :: Data.Sequence.Seq [a] -> Data.Sequence.Seq [a] tail' = Data.Sequence.drop 1 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 @@ -28,7 +28,6 @@ Category(..) ) where -import qualified Data.Ratio import qualified Factory.Math.Precision as Math.Precision import qualified ToolShed.Defaultable @@ -41,9 +40,9 @@ * Since representing /Pi/ as either a 'Rational' or promoted to an 'Integer', is inconvenient, an alternative decimal 'String'-representation is provided. -} class Algorithmic algorithm where - openR :: algorithm -> Math.Precision.DecimalDigits -> Data.Ratio.Rational -- ^ Returns the value of /Pi/ as a 'Rational'. + openR :: algorithm -> Math.Precision.DecimalDigits -> 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. + openI :: algorithm -> Math.Precision.DecimalDigits -> Integer -- ^ Returns the value of /Pi/, promoted by the required precision to form an integer. openI _ 1 = 3 openI algorithm decimalDigits | decimalDigits <= 0 = error $ "Factory.Math.Pi.openI:\tinsufficient decimalDigits=" ++ show decimalDigits diff --git a/src/Factory/Math/Precision.hs b/src/Factory/Math/Precision.hs --- a/src/Factory/Math/Precision.hs +++ b/src/Factory/Math/Precision.hs @@ -33,6 +33,7 @@ -- * Functions getIterationsRequired, getTermsRequired, + roundTo, promote, simplify ) where @@ -45,7 +46,7 @@ -- | The /rate of convergence/; . type ConvergenceRate = Double --- | A number of decimal digits. +-- | A number of decimal digits; presumably positive. type DecimalDigits = Int -- | /Linear/ convergence-rate; which may be qualified by the /rate of convergence/. @@ -99,20 +100,26 @@ | requiredDecimalDigits < 0 = error $ "Factory.Math.Precision.getTermsRequired:\t'requiredDecimalDigits' must be positive; " ++ show requiredDecimalDigits | otherwise = ceiling $ fromIntegral requiredDecimalDigits / negate (logBase 10 convergenceRate) --- | Promotes the specified number, by a number of 'DecimalDigits'. +-- | Rounds the specified number, to a positive number of 'DecimalDigits'. +roundTo :: (RealFrac a, Fractional f) => DecimalDigits -> a -> f +roundTo decimals = (/ fromInteger promotionFactor) . fromInteger . round . (* fromInteger promotionFactor) where + promotionFactor :: Integer + promotionFactor = 10 ^ decimals + +-- | Promotes the specified number, by a positive number of 'DecimalDigits'. promote :: Num n => n -> DecimalDigits -> n promote x = (* x) . (10 ^) {- | - * Reduces a 'Data.Ratio.Rational' to the minimal form required for the specified number of /fractional/ decimal places; + * Reduces a 'Rational' to the minimal form required for the specified number of /fractional/ decimal places; irrespective of the number of integral decimal places. - * A 'Data.Ratio.Rational' approximation to an irrational number, may be very long, and provide an unknown excess precision. + * A 'Rational' approximation to an irrational number, may be very long, and provide an unknown excess precision. Whilst this doesn't sound harmful, it costs in performance and memory-requirement, and being unpredictable isn't actually useful. -} simplify :: RealFrac operand => DecimalDigits -- ^ The number of places after the decimal point, which are required. -> operand - -> Data.Ratio.Rational + -> Rational simplify decimalDigits operand = Data.Ratio.approxRational operand . recip $ 4 * 10 ^ succ decimalDigits --Tolerate any error less than half the least significant digit required. diff --git a/src/Factory/Math/Probability.hs b/src/Factory/Math/Probability.hs --- a/src/Factory/Math/Probability.hs +++ b/src/Factory/Math/Probability.hs @@ -1,5 +1,5 @@ {- - Copyright (C) 2011 Dr. Alistair Ward + Copyright (C) 2011-2013 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 @@ -17,61 +17,135 @@ {- | [@AUTHOR@] Dr. Alistair Ward - [@DESCRIPTION@] Miscellaneous functions for probability-distributions. + [@DESCRIPTION@] Functions for probability-distributions. + + [@CAVEAT@] Because data-constructors are exposed, 'ToolShed.SelfValidate.isValid' need not be called. -} module Factory.Math.Probability( +-- * Type-classes + Distribution(..), -- * Types -- ** Data-types ContinuousDistribution(..), DiscreteDistribution(..), -- * Functions - boxMullerTransform, + maxPreciseInteger, -- minPositiveFloat, + boxMullerTransform, -- reProfile, generateStandardizedNormalDistribution, generateContinuousPopulation, - generatePoissonDistribution, +-- generatePoissonDistribution, generateDiscretePopulation ) where import qualified Control.Arrow import Control.Arrow((***), (&&&)) import qualified Factory.Data.Interval as Data.Interval +import qualified Factory.Math.Power as Math.Power import qualified System.Random import qualified ToolShed.Data.List import qualified ToolShed.Data.Pair import qualified ToolShed.SelfValidate --- | Describes a /continuous probability-distribution/; . -data ContinuousDistribution f - = UniformDistribution (Data.Interval.Interval f) -- ^ Defines a /Uniform/-distribution within a /closed interval/; . - | NormalDistribution f f -- ^ Defines a /Normal/-distribution with a particular /mean/ and /variance/; . +-- | The maximum integer which can be accurately represented as a Double. +maxPreciseInteger :: RealFloat a => a -> Integer +maxPreciseInteger = (2 ^) . floatDigits + +{- | + * Determines the minimum positive floating-point number, which can be represented by using the parameter's type. + + * Only the type of the parameter is relevant, not its value. +-} +minPositiveFloat :: RealFloat a => a -> a +minPositiveFloat = encodeFloat 1 . uncurry (-) . (fst . floatRange &&& floatDigits) + +-- | Describes /continuous probability-distributions/; . +data ContinuousDistribution parameter + = ExponentialDistribution parameter {-lambda-} -- ^ Defines an /Exponential/-distribution with a particular /lambda/; . + | LogNormalDistribution parameter {-location-} parameter {-scale2-} -- ^ Defines a distribution whose logarithm is normally distributed with a particular /mean/ & /variance/; . + | NormalDistribution parameter {-mean-} parameter {-variance-} -- ^ Defines a /Normal/-distribution with a particular /mean/ & /variance/; . + | UniformDistribution (Data.Interval.Interval parameter) -- ^ Defines a /Uniform/-distribution within a /closed interval/; . deriving (Eq, Read, Show) -instance (Num a, Ord a, Show a) => ToolShed.SelfValidate.SelfValidator (ContinuousDistribution a) where - getErrors distribution = ToolShed.SelfValidate.extractErrors $ case distribution of - UniformDistribution interval -> [(Data.Interval.isReversed interval, "Reversed interval='" ++ show interval ++ "'.")] - NormalDistribution _ v -> [(v < 0, "Negative variance=" ++ show v ++ ".")] +instance (Floating parameter, Ord parameter, Show parameter) => ToolShed.SelfValidate.SelfValidator (ContinuousDistribution parameter) where + getErrors probabilityDistribution = ToolShed.SelfValidate.extractErrors $ case probabilityDistribution of + ExponentialDistribution lambda -> [(lambda <= 0, "'lambda' must exceed zero; " ++ show probabilityDistribution ++ ".")] + LogNormalDistribution location scale2 -> let + maxParameter = log . fromInteger $ maxPreciseInteger (undefined :: Double) + in [ + (scale2 <= 0, "'scale' must exceed zero; " ++ show probabilityDistribution ++ "."), + (location > maxParameter || scale2 > maxParameter, "loss of precision will result from either 'location' or 'scale^2' exceeding '" ++ show maxParameter ++ "'; " ++ show probabilityDistribution ++ ".") + ] + NormalDistribution _ variance -> [(variance <= 0, "variance must exceed zero; " ++ show probabilityDistribution ++ ".")] + UniformDistribution interval -> [(Data.Interval.isReversed interval, "reversed interval='" ++ show probabilityDistribution ++ "'.")] --- | Describes a /discrete probability-distribution/; . -data DiscreteDistribution f = PoissonDistribution f deriving (Eq, Read, Show) +-- | Describes /discrete probability-distributions/; . +data DiscreteDistribution parameter + = PoissonDistribution parameter {-lambda-} -- ^ Defines an /Poisson/-distribution with a particular /lambda/; . + | ShiftedGeometricDistribution parameter {-probability-} -- ^ Defines an /Geometric/-distribution with a particular probability of success; . + deriving (Eq, Read, Show) -instance (Num f, Ord f, Show f) => ToolShed.SelfValidate.SelfValidator (DiscreteDistribution f) where - getErrors (PoissonDistribution lambda) = ToolShed.SelfValidate.extractErrors [(lambda < 0, "Negative lambda=" ++ show lambda ++ ".")] +instance (Num parameter, Ord parameter, Show parameter) => ToolShed.SelfValidate.SelfValidator (DiscreteDistribution parameter) where + getErrors probabilityDistribution = ToolShed.SelfValidate.extractErrors $ case probabilityDistribution of + PoissonDistribution lambda -> [(lambda <= 0, "'lambda' must exceed zero; " ++ show probabilityDistribution ++ ".")] + ShiftedGeometricDistribution probability -> [(any ($ probability) [(<= 0), (> 1)], "probability must be in the semi-closed unit-interval (0, 1]; " ++ show probabilityDistribution ++ ".")] +-- | Defines a common interface for probability-distributions. +class Distribution probabilityDistribution where + generatePopulation + :: (Fractional sample, System.Random.RandomGen randomGen) + => probabilityDistribution + -> randomGen -- ^ A generator of /uniformly distributed/ random numbers. + -> [sample] -- ^ CAVEAT: the integers generated for discrete distributions are represented by a fractional type; use 'generateDiscretePopulation' if this is a problem. + + getMean :: Fractional mean => probabilityDistribution -> mean -- ^ The theoretical mean. + + getStandardDeviation :: Floating standardDeviation => probabilityDistribution -> standardDeviation-- ^ The theoretical standard-deviation. + getStandardDeviation = sqrt . getVariance --Default implementation. + + getVariance :: Floating variance => probabilityDistribution -> variance -- ^ The theoretical variance. + getVariance = Math.Power.square . getStandardDeviation --Default implementation. + +instance (RealFloat parameter, Show parameter, System.Random.Random parameter) => Distribution (ContinuousDistribution parameter) where + generatePopulation probabilityDistribution = map realToFrac {-parameter -> sample-} . generateContinuousPopulation probabilityDistribution + + getMean (ExponentialDistribution lambda) = realToFrac $ recip lambda + getMean (LogNormalDistribution location scale2) = realToFrac . exp . (+ location) $ scale2 / 2 + getMean (NormalDistribution mean _) = realToFrac mean + getMean (UniformDistribution (minParameter, maxParameter)) = realToFrac $ (minParameter + maxParameter) / 2 + + getVariance (ExponentialDistribution lambda) = realToFrac . recip $ Math.Power.square lambda + getVariance (LogNormalDistribution location scale2) = realToFrac $ (exp scale2 - 1) * exp (2 * location + scale2) + getVariance (NormalDistribution _ variance) = realToFrac variance + getVariance (UniformDistribution (minParameter, maxParameter)) = realToFrac $ Math.Power.square (maxParameter - minParameter) / 12 + +instance (RealFloat parameter, Show parameter, System.Random.Random parameter) => Distribution (DiscreteDistribution parameter) where + generatePopulation probabilityDistribution = map fromInteger . generateDiscretePopulation probabilityDistribution + + getMean (PoissonDistribution lambda) = realToFrac lambda + getMean (ShiftedGeometricDistribution probability) = realToFrac $ recip probability + + getVariance (PoissonDistribution lambda) = realToFrac lambda + getVariance (ShiftedGeometricDistribution probability) = realToFrac $ (1 - probability) / Math.Power.square probability + {- | * Converts a pair of independent /uniformly distributed/ random numbers, within the /semi-closed unit interval/ /(0,1]/, to a pair of independent /normally distributed/ random numbers, of standardized /mean/=0, and /variance/=1. * . -} -boxMullerTransform :: (Floating f, Ord f, Show f) +boxMullerTransform :: ( + Floating f, + Ord f, + Show f + ) => (f, f) -- ^ Independent, /uniformly distributed/ random numbers, which must be within the /semi-closed unit interval/, /(0,1]/. -> (f, f) -- ^ Independent, /normally distributed/ random numbers, with standardized /mean/=0 and /variance/=1. boxMullerTransform cartesian | not . uncurry (&&) $ ToolShed.Data.Pair.mirror inSemiClosedUnitInterval cartesian = error $ "Factory.Math.Probability.boxMullerTransform:\tspecified Cartesian coordinates, must be within semi-closed unit-interval (0, 1]; " ++ show cartesian - | otherwise = polarToCartesianTransform $ (sqrt . negate . (* 2) . log *** (* 2) . (* pi)) cartesian + | otherwise = polarToCartesianTransform $ (sqrt . negate . (* 2) . log *** (* 2) . (* pi)) cartesian where inSemiClosedUnitInterval :: (Num n, Ord n) => n -> Bool inSemiClosedUnitInterval = uncurry (&&) . ((> 0) &&& (<= 1)) @@ -80,14 +154,6 @@ polarToCartesianTransform = uncurry (*) . Control.Arrow.second cos &&& uncurry (*) . Control.Arrow.second sin {- | - * Determines the minimum positive floating-point number, which can be represented by using the parameter's type. - - * Only the type of the parameter is relevant, not its value. --} -minPositiveFloat :: RealFloat a => a -> a -minPositiveFloat = encodeFloat 1 . uncurry (-) . (fst . floatRange &&& floatDigits) - -{- | * Uses the supplied random-number generator, to generate a conceptually infinite list, of /normally distributed/ random numbers, with standardized /mean/=0, and /variance/=1. @@ -103,50 +169,45 @@ System.Random.randomRs (minPositiveFloat undefined, 1) ) . System.Random.split --- | Stretches and shifts a /standardized normal distribution/ to achieve the required /mean/ and /standard-deviation/. -reProfile :: Num n => n -> n -> [n] -> [n] -reProfile mean standardDeviation = map ((+ mean) . (* standardDeviation)) - -{- | - * Generates a random sample-population, with the specified continuous probability-distribution. +-- | Stretches and shifts a /distribution/ to achieve the required /mean/ and /standard-deviation/. +reProfile :: (Distribution distribution, Floating n) => distribution -> [n] -> [n] +reProfile distribution = map ((+ getMean distribution) . (* getStandardDeviation distribution)) - * When a /Normal distribution/ is requested, - the generated population will only tend towards the requested /mean/ and /variance/ of, as the sample-size tends towards infinity. - Whilst one could arrange for these criteria to be precisely met for any sample-size, the sample would lose a degree of randomness as a result. --} +-- | Uses the supplied random-number generator, to generate a conceptually infinite population, with the specified continuous probability-distribution. generateContinuousPopulation :: ( RealFloat f, Show f, System.Random.Random f, System.Random.RandomGen randomGen ) - => Int -- ^ number of items. - -> ContinuousDistribution f + => ContinuousDistribution f -> randomGen -- ^ A generator of /uniformly distributed/ random numbers. -> [f] -generateContinuousPopulation 0 _ _ = [] -generateContinuousPopulation populationSize probabilityDistribution randomGen - | populationSize < 0 = error $ "Factory.Math.Probability.generateContinuousPopulation:\tinvalid population-size=" ++ show populationSize +generateContinuousPopulation probabilityDistribution randomGen | not $ ToolShed.SelfValidate.isValid probabilityDistribution = error $ "Factory.Math.Probability.generateContinuousPopulation:\t" ++ ToolShed.SelfValidate.getFirstError probabilityDistribution - | otherwise = take populationSize $ ( + | otherwise = ( case probabilityDistribution of - UniformDistribution interval -> System.Random.randomRs interval - NormalDistribution requiredMean requiredVariance -> reProfile requiredMean (sqrt requiredVariance) . generateStandardizedNormalDistribution + ExponentialDistribution lambda -> let + quantile = (/ lambda) . negate . log . (1 -) -- . + in map quantile . System.Random.randomRs (0, 1) + LogNormalDistribution location scale2 -> map ( + exp . (+ location) . (* sqrt scale2) --Stretch the standard-deviation & re-locate the mean to that specified for the log-space, then return to the original coordinates. + ) . generateStandardizedNormalDistribution + NormalDistribution _ _ -> reProfile probabilityDistribution . generateStandardizedNormalDistribution + UniformDistribution interval -> System.Random.randomRs interval ) randomGen {- | * Uses the supplied random-number generator, - to generate a conceptually infinite list, of random integers conforming to the /Poisson distribution/ (/mean/=lambda, /variance/=lambda). - - * . + to generate a conceptually infinite population, of random integers conforming to the /Poisson distribution/; . * CAVEAT: uses an algorithm by Knuth, which having a /linear time-complexity/ in /lambda/, can be intolerably slow; also, the term @exp $ negate lambda@, underflows for large /lambda/; - so for large /lambda/, this implementation returns the appropriate 'NormalDistribution', which is similar for large /lambda/. + so for large /lambda/, this implementation returns the appropriate 'NormalDistribution'. -} generatePoissonDistribution :: ( - Integral events, + Integral sample, RealFloat lambda, Show lambda, System.Random.Random lambda, @@ -154,12 +215,12 @@ ) => lambda -- ^ Defines the required approximate value of both /mean/ and /variance/. -> randomGen - -> [events] + -> [sample] generatePoissonDistribution lambda - | lambda < 0 = error $ "Factory.Math.Probability.generatePoissonDistribution:\tinvalid lambda=" ++ show lambda + | lambda <= 0 = error $ "Factory.Math.Probability.generatePoissonDistribution:\tlambda must exceed zero " ++ show lambda | lambda > ( negate . log $ minPositiveFloat lambda --Guard against underflow, in the user-defined type for lambda. - ) = filter (>= 0) . map round . reProfile lambda (sqrt lambda) . generateStandardizedNormalDistribution + ) = filter (>= 0) . map round . (reProfile (PoissonDistribution lambda) :: [Double] -> [Double]) . generateStandardizedNormalDistribution | otherwise = generator where generator = uncurry (:) . ( @@ -170,25 +231,25 @@ ) (negate 1, 1) . System.Random.randomRs (0, 1) *** generator {-recurse-} ) . System.Random.split --- | Generates a random sample-population, with the specified discrete probability-distribution. +-- | Uses the supplied random-number generator, to generate a conceptually infinite population, with the specified discrete probability-distribution. generateDiscretePopulation :: ( - Ord f, - RealFloat f, - Show f, - System.Random.Random f, - System.Random.RandomGen randomGen, - Integral events + Integral sample, + Ord parameter, + RealFloat parameter, + Show parameter, + System.Random.Random parameter, + System.Random.RandomGen randomGen ) - => Int -- ^ number of items. - -> DiscreteDistribution f + => DiscreteDistribution parameter -> randomGen -- ^ A generator of /uniformly distributed/ random numbers. - -> [events] -generateDiscretePopulation 0 _ _ = [] -generateDiscretePopulation populationSize probabilityDistribution randomGen - | populationSize < 0 = error $ "Factory.Math.Probability.generateDiscretePopulation:\tinvalid populationSize=" ++ show populationSize + -> [sample] +generateDiscretePopulation probabilityDistribution randomGen | not $ ToolShed.SelfValidate.isValid probabilityDistribution = error $ "Factory.Math.Probability.generateDiscretePopulation:\t" ++ ToolShed.SelfValidate.getFirstError probabilityDistribution - | otherwise = take populationSize $ ( + | otherwise = ( case probabilityDistribution of PoissonDistribution lambda -> generatePoissonDistribution lambda + ShiftedGeometricDistribution probability + | probability == 1 -> const $ repeat 1 --The first Bernoulli Trial is guaranteed to succeed. + | otherwise -> map ceiling {-minimum 1-} . (\x -> x :: [Rational]) . generatePopulation (ExponentialDistribution . negate $ log (1 - probability)) --The geometric distribution is a discrete version of the exponential distribution. ) randomGen diff --git a/src/Factory/Math/Radix.hs b/src/Factory/Math/Radix.hs --- a/src/Factory/Math/Radix.hs +++ b/src/Factory/Math/Radix.hs @@ -44,7 +44,7 @@ -- | Constant random-access lookup for 'digits'. encodes :: (Data.Array.IArray.Ix index, Integral index) => Data.Array.IArray.Array index Char -encodes = Data.Array.IArray.listArray (0, fromIntegral . pred $ length digits) digits where +encodes = Data.Array.IArray.listArray (0, pred $ Data.List.genericLength digits) digits -- | Constant reverse-lookup for 'digits'. decodes :: Integral i => [(Char, i)] 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 @@ -41,12 +41,11 @@ isPrecise ) where -import qualified Data.Ratio import qualified Factory.Math.Power as Math.Power import qualified Factory.Math.Precision as Math.Precision -- | The result-type; actually, only the concrete return-type of 'Math.Precision.simplify', stops it being a polymorphic instance of 'Fractional'. -type Result = Data.Ratio.Rational +type Result = Rational -- | Contains an estimate for the /square-root/ of a value, and its accuracy. type Estimate = (Result, Math.Precision.DecimalDigits) 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 @@ -36,7 +36,6 @@ import Control.Parallel(par, pseq) import qualified Data.Foldable import qualified Data.List -import qualified Data.Ratio import qualified Factory.Math.Factorial as Math.Factorial import qualified Factory.Math.Implementations.Factorial as Math.Implementations.Factorial import qualified Factory.Math.Power as Math.Power @@ -44,7 +43,7 @@ {- | * Determines the /mean/ of the specified numbers; . - * Should the caller define the result-type as 'Data.Ratio.Rational', then it will be free from rounding-errors. + * Should the caller define the result-type as 'Rational', then it will be free from rounding-errors. -} getMean :: (Data.Foldable.Foldable f, Real r, Fractional result) => f r -> result getMean x @@ -56,22 +55,22 @@ {- | * 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. + * Should the caller define the result-type as 'Rational', then it will be free from rounding-errors. -} getDispersionFromMean :: ( Data.Foldable.Foldable f, Fractional result, Functor f, Real r - ) => (Data.Ratio.Rational -> Data.Ratio.Rational) -> f r -> result + ) => (Rational -> Rational) -> f r -> result getDispersionFromMean weight x = getMean $ fmap (weight . (+ negate mean) . toRational) x where - mean :: Data.Ratio.Rational + mean :: Rational mean = getMean x {- | * Determines the exact /variance/ of the specified numbers; . - * Should the caller define the result-type as 'Data.Ratio.Rational', then it will be free from rounding-errors. + * Should the caller define the result-type as 'Rational', then it will be free from rounding-errors. -} getVariance :: ( Data.Foldable.Foldable f, @@ -93,7 +92,7 @@ {- | * Determines the /average absolute deviation/ of the specified numbers; . - * Should the caller define the result-type as 'Data.Ratio.Rational', then it will be free from rounding-errors. + * Should the caller define the result-type as 'Rational', then it will be free from rounding-errors. -} getAverageAbsoluteDeviation :: ( Data.Foldable.Foldable f, 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 @@ -1,6 +1,6 @@ {-# LANGUAGE CPP #-} {- - Copyright (C) 2010 Dr. Alistair Ward + 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 @@ -28,14 +28,14 @@ sumR ) where +import qualified Control.DeepSeq import qualified Data.List import qualified Data.Ratio import Data.Ratio((%)) -import qualified Control.DeepSeq +import qualified ToolShed.Data.List #if MIN_VERSION_parallel(3,0,0) import qualified Control.Parallel.Strategies -import qualified ToolShed.Data.List #endif {- | @@ -44,16 +44,12 @@ * Sparks the summation of @(list-length / chunk-size)@ chunks from the list, each of the specified size (thought the last chunk may be smaller), then recursively sums the list of results from each spark. - * CAVEAT: unless the numbers are large, 'Data.Ratio.Rational' (requiring /cross-multiplication/), or the list long, + * CAVEAT: unless the numbers are large, 'Rational' (requiring /cross-multiplication/), or the list long, 'sum' is too light-weight for sparking to be productive, therefore it is more likely to be the parallelised deep /evaluation/ of list-elements which saves time. -} sum' :: (Num n, Control.DeepSeq.NFData n) -#if MIN_VERSION_toolshed(11,1,1) => ToolShed.Data.List.ChunkLength -#else - => Int -- ^ The Chunk-length. -#endif -> [n] -> n #if MIN_VERSION_parallel(3,0,0) @@ -90,11 +86,7 @@ -} {-# INLINE sumR #-} --This makes a staggering difference to calls from other modules. sumR :: (Integral i, Control.DeepSeq.NFData i) -#if MIN_VERSION_toolshed(11,1,1) => ToolShed.Data.List.ChunkLength -#else - => Int -- ^ The Chunk-length. -#endif -> [Data.Ratio.Ratio i] -> Data.Ratio.Ratio i sumR chunkLength 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 @@ -35,7 +35,7 @@ -- | 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 + | i < 0 = fail $ "Factory.Test.Performance.Primality.carmichaelNumbersPerformance:\tnegative number; " ++ show i | otherwise = ToolShed.System.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. 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 @@ -43,7 +43,7 @@ -} primeFactorsPerformanceGraph :: Math.PrimeFactorisation.Algorithmic algorithm => algorithm -> Int -> IO () primeFactorsPerformanceGraph algorithm tests - | tests < 0 = error $ "Factory.Test.Performance.PrimeFactorisation.primeFactorsPerformanceGraph:\tnegative number; " ++ show tests + | tests < 0 = fail $ "Factory.Test.Performance.PrimeFactorisation.primeFactorsPerformanceGraph:\tnegative number; " ++ show tests | otherwise = mapM_ ( \operand -> primeFactorsPerformance algorithm operand >>= putStrLn . shows operand . showChar '\t' . (`shows` "") ) . take tests . dropWhile (< 2) $ Math.Fibonacci.fibonacci diff --git a/src/Factory/Test/Performance/Primes.hs b/src/Factory/Test/Performance/Primes.hs --- a/src/Factory/Test/Performance/Primes.hs +++ b/src/Factory/Test/Performance/Primes.hs @@ -43,5 +43,5 @@ -- | Measures the CPU-time required to find the specified number of /Mersenne/-numbers, which is returned together with the requested list. mersenneNumbersPerformance :: Math.Primes.Algorithmic algorithm => algorithm -> Int -> IO (Double, [Integer]) mersenneNumbersPerformance primalityAlgorithm i - | i < 0 = error $ "Factory.Test.Performance.Primes.mersenneNumbersPerformance:\tnegative number; " ++ show i + | i < 0 = fail $ "Factory.Test.Performance.Primes.mersenneNumbersPerformance:\tnegative number; " ++ show i | otherwise = ToolShed.System.TimePure.getCPUSeconds . take i $ Math.Primes.mersenneNumbers primalityAlgorithm diff --git a/src/Factory/Test/QuickCheck/Polynomial.hs b/src/Factory/Test/QuickCheck/Polynomial.hs --- a/src/Factory/Test/QuickCheck/Polynomial.hs +++ b/src/Factory/Test/QuickCheck/Polynomial.hs @@ -31,7 +31,6 @@ import Control.Arrow((***)) import Factory.Data.Ring((=*=), (=+=), (=-=), (=^)) import qualified Data.Numbers.Primes -import qualified Data.Ratio import qualified Factory.Data.Polynomial as Data.Polynomial import qualified Factory.Data.QuotientRing as Data.QuotientRing import qualified Factory.Data.Ring as Data.Ring @@ -65,14 +64,14 @@ prop_quotRem, prop_degree, prop_ringNormalised, prop_quotientRingNormalised :: Data.Polynomial.Polynomial Integer Integer -> Data.Polynomial.Polynomial Integer Integer -> Test.QuickCheck.Property prop_quotRem numerator denominator = denominator' /= Data.Polynomial.zero ==> Test.QuickCheck.label "prop_quotRem" $ numerator' == denominator' =*= quotient =+= remainder where - numerator', denominator' :: Data.Polynomial.Polynomial Data.Ratio.Rational Integer + numerator', denominator' :: Data.Polynomial.Polynomial Rational Integer numerator' = Data.Polynomial.realCoefficientsToFrac numerator denominator' = Data.Polynomial.realCoefficientsToFrac denominator (quotient, remainder) = numerator' `Data.QuotientRing.quotRem'` denominator' prop_degree numerator denominator = denominator' /= Data.Polynomial.zero ==> Test.QuickCheck.label "prop_degree" $ remainder == Data.Polynomial.zero || Data.Polynomial.getDegree remainder < Data.Polynomial.getDegree denominator' where - numerator', denominator' :: Data.Polynomial.Polynomial Data.Ratio.Rational Integer + numerator', denominator' :: Data.Polynomial.Polynomial Rational Integer numerator' = Data.Polynomial.realCoefficientsToFrac numerator denominator' = Data.Polynomial.realCoefficientsToFrac denominator @@ -81,7 +80,7 @@ prop_ringNormalised l r = Test.QuickCheck.label "prop_ringNormalised" $ all Data.Polynomial.isNormalised [l =*= r, l =+= r, l =-= r] prop_quotientRingNormalised numerator denominator = denominator' /= Data.Polynomial.zero ==> Test.QuickCheck.label "prop_quotientRingNormalised" $ all Data.Polynomial.isNormalised [numerator' `Data.QuotientRing.quot'` denominator', numerator' `Data.QuotientRing.rem'` denominator'] where - numerator', denominator' :: Data.Polynomial.Polynomial Data.Ratio.Rational Integer + numerator', denominator' :: Data.Polynomial.Polynomial Rational Integer numerator' = Data.Polynomial.realCoefficientsToFrac numerator denominator' = Data.Polynomial.realCoefficientsToFrac denominator @@ -91,7 +90,7 @@ power' = succ $ power `mod` 100 prop_perfectPower polynomial power = polynomial' /= Data.Polynomial.zero ==> Test.QuickCheck.label "prop_perfectPower" $ iterate (`Data.QuotientRing.quot'` polynomial') (polynomial' =^ power') !! pred power' == polynomial' where - polynomial' :: Data.Polynomial.Polynomial Data.Ratio.Rational Integer + polynomial' :: Data.Polynomial.Polynomial Rational Integer polynomial' = Data.Polynomial.realCoefficientsToFrac polynomial power' :: Int @@ -117,6 +116,6 @@ prop_integralDomain polynomials = Data.Polynomial.zero `notElem` polynomials ==> Test.QuickCheck.label "prop_integralDomain" $ Data.Ring.product' (recip 2) {-TODO-} 10 polynomials /= Data.Polynomial.zero prop_isDivisibleBy polynomials = Test.QuickCheck.label "prop_isDivisibleBy" . all (Data.QuotientRing.isDivisibleBy (Data.Ring.product' (recip 2) {-TODO-} 10 polynomials')) $ filter (/= Data.Polynomial.zero) polynomials' where - polynomials' :: [Data.Polynomial.Polynomial Data.Ratio.Rational Integer] + polynomials' :: [Data.Polynomial.Polynomial Rational Integer] polynomials' = map Data.Polynomial.realCoefficientsToFrac polynomials 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 @@ -1,5 +1,5 @@ {- - Copyright (C) 2011 Dr. Alistair Ward + Copyright (C) 2011-2013 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 @@ -22,44 +22,124 @@ module Factory.Test.QuickCheck.Probability( -- * Functions +-- normalise, quickChecks ) where import Control.Arrow((&&&)) -import qualified Factory.Math.Probability as Math.Probability -import qualified Factory.Math.Statistics as Math.Statistics +import qualified Data.List +import qualified Factory.Math.Probability as Math.Probability +import qualified Factory.Math.Statistics as Math.Statistics import Factory.Test.QuickCheck.Factorial() import qualified System.Random import qualified Test.QuickCheck import Test.QuickCheck((==>)) import qualified ToolShed.Data.Pair +-- | Re-profile a distribution to achieve a standard mean & variance. +normalise :: ( + Eq f, + Floating f, + Math.Probability.Distribution distribution + ) => distribution -> [f] -> [f] +normalise distribution + | variance == 0 = error "Factory.Test.Quick.Probability.normalise:\tzero variance => can't stretch to one." + | otherwise = map $ (/ sqrt variance) . (+ negate mean) + where + (mean, variance) = Math.Probability.getMean &&& Math.Probability.getVariance $ distribution + -- | Defines invariant properties. quickChecks :: IO () quickChecks = do randomGen <- System.Random.getStdGen - Test.QuickCheck.quickCheck (prop_normalDistribution randomGen) >> Test.QuickCheck.quickCheck (prop_poissonDistribution randomGen) where - prop_normalDistribution :: System.Random.RandomGen randomGen => randomGen -> (Double, Double) -> Test.QuickCheck.Property - prop_normalDistribution randomGen (mean, variance) = variance' /= 0 ==> Test.QuickCheck.label "prop_normalDistribution" . uncurry (&&) . ToolShed.Data.Pair.mirror ( - (< (0.1 :: Double)) . abs --Generous tolerance. + (Test.QuickCheck.quickCheck . ($ randomGen)) `mapM_` [ + prop_logNormalDistribution, + prop_logNormalDistribution', + prop_normalDistribution, + prop_uniformDistribution + ] >> (Test.QuickCheck.quickCheck . ($ randomGen)) `mapM_` [ + prop_exponentialDistribution, + prop_exponentialDistribution', + prop_poissonDistribution, + prop_poissonDistribution', + prop_shiftedGeometricDistribution, + prop_shiftedGeometricDistribution' + ] + where + isWithinTolerance :: Double -> Double -> Bool + isWithinTolerance i = (< recip i) . abs + + prop_logNormalDistribution, prop_logNormalDistribution', prop_normalDistribution, prop_uniformDistribution :: System.Random.RandomGen randomGen => randomGen -> Double -> Double -> Test.QuickCheck.Property + prop_logNormalDistribution randomGen location scale2 = scale2 /= 0 ==> Test.QuickCheck.label "prop_logNormalDistribution" . uncurry (&&) . ToolShed.Data.Pair.mirror (isWithinTolerance 1) . ( + Math.Statistics.getMean &&& pred . Math.Statistics.getStandardDeviation --Both of which, having been normalised, should be zero. ) . ( - Math.Statistics.getMean &&& pred . Math.Statistics.getStandardDeviation - ) . map ( - (/ sqrt variance') . (+ negate mean) --Standardize. - ) $ Math.Probability.generateContinuousPopulation 1000 (Math.Probability.NormalDistribution mean variance') randomGen where - variance' = abs variance + normalise distribution :: [Double] -> [Double] + ) . take 10000 $ Math.Probability.generatePopulation distribution randomGen where + maxParameter = log . fromInteger $ Math.Probability.maxPreciseInteger (undefined :: Double) + location' + | location >= 0 = maxParameter `min` location + | otherwise = negate maxParameter `max` location - prop_poissonDistribution :: System.Random.RandomGen randomGen => randomGen -> Int -> Test.QuickCheck.Property - prop_poissonDistribution randomGen lambda = lambda' /= 0 ==> Test.QuickCheck.label "prop_poissonDistribution" . uncurry (&&) . ToolShed.Data.Pair.mirror ( - (< (0.1 :: Double)) . abs --Tolerance. + distribution = Math.Probability.LogNormalDistribution location' . min maxParameter $ abs scale2 + + prop_logNormalDistribution' randomGen location scale2 = scale2 /= 0 ==> Test.QuickCheck.label "prop_logNormalDistribution'" . all ( + >= (0 :: Double) + ) . take 10 $ Math.Probability.generatePopulation (Math.Probability.LogNormalDistribution location' . min maxParameter $ abs scale2) randomGen where + maxParameter = log . fromInteger $ Math.Probability.maxPreciseInteger (undefined :: Double) + + location' + | location >= 0 = maxParameter `min` location + | otherwise = negate maxParameter `max` location + + prop_normalDistribution randomGen mean variance = variance /= 0 ==> Test.QuickCheck.label "prop_normalDistribution" . uncurry (&&) . ToolShed.Data.Pair.mirror (isWithinTolerance 10) . ( + Math.Statistics.getMean &&& pred . Math.Statistics.getStandardDeviation --Both of which, having been normalised, should be zero. ) . ( - Math.Statistics.getMean &&& pred . Math.Statistics.getStandardDeviation - ) $ map ( - (/ sqrt lambda') . (+ negate lambda') . fromIntegral --Standardize. - ) ( - Math.Probability.generateDiscretePopulation 1000 (Math.Probability.PoissonDistribution lambda') randomGen :: [Int] - ) where - lambda' :: Double - lambda' = fromIntegral $ mod lambda 1000 + normalise distribution :: [Double] -> [Double] + ) . take 1000 $ Math.Probability.generatePopulation distribution randomGen where + distribution = Math.Probability.NormalDistribution mean $ abs variance + + prop_uniformDistribution randomGen min' max' = min' /= max' ==> Test.QuickCheck.label "prop_uniformDistribution" . uncurry (&&) . ToolShed.Data.Pair.mirror (isWithinTolerance 10) . ( + Math.Statistics.getMean &&& pred . Math.Statistics.getStandardDeviation --Both of which, having been normalised, should be zero. + ) . ( + normalise distribution :: [Double] -> [Double] + ) . take 10000 $ Math.Probability.generatePopulation distribution randomGen where + [min'', max''] = Data.List.sort [min', max'] + distribution = Math.Probability.UniformDistribution (min'', max'') + + prop_exponentialDistribution, prop_exponentialDistribution', prop_poissonDistribution, prop_poissonDistribution', prop_shiftedGeometricDistribution, prop_shiftedGeometricDistribution' :: System.Random.RandomGen randomGen => randomGen -> Double -> Test.QuickCheck.Property + prop_exponentialDistribution randomGen lambda = Test.QuickCheck.label "prop_exponentialDistribution" . uncurry (&&) . ToolShed.Data.Pair.mirror (isWithinTolerance 10) . ( + Math.Statistics.getMean &&& pred . Math.Statistics.getStandardDeviation --Both of which, having been normalised, should be zero. + ) . ( + normalise distribution :: [Double] -> [Double] + ) . take 10000 $ Math.Probability.generatePopulation distribution randomGen where + distribution = Math.Probability.ExponentialDistribution . succ {-exclude zero-} $ abs lambda `max` 10 {-cap-} + + prop_exponentialDistribution' randomGen lambda = lambda /= 0 ==> Test.QuickCheck.label "prop_exponentialDistribution'" . all ( + >= (0 :: Double) + ) . take 10 $ Math.Probability.generatePopulation (Math.Probability.ExponentialDistribution $ abs lambda) randomGen + + prop_poissonDistribution randomGen lambda = Test.QuickCheck.label "prop_poissonDistribution" . uncurry (&&) . ToolShed.Data.Pair.mirror (isWithinTolerance 10) . ( + Math.Statistics.getMean &&& pred . Math.Statistics.getStandardDeviation --Both of which, having been normalised, should be zero. + ) . ( + normalise distribution :: [Double] -> [Double] + ) . take 1000 $ Math.Probability.generatePopulation distribution randomGen where + distribution = Math.Probability.PoissonDistribution . succ {-exclude zero-} $ abs lambda `max` 10 {-cap-} + + prop_poissonDistribution' randomGen lambda = lambda /= 0 ==> Test.QuickCheck.label "prop_poissonDistribution'" . all ( + >= (0 :: Double) + ) . take 10 $ Math.Probability.generatePopulation (Math.Probability.PoissonDistribution $ abs lambda) randomGen + + prop_shiftedGeometricDistribution randomGen probability = probability' /= 1 ==> Test.QuickCheck.label "prop_shiftedGeometricDistribution" . uncurry (&&) . ToolShed.Data.Pair.mirror (isWithinTolerance 10) . ( + Math.Statistics.getMean &&& pred . Math.Statistics.getStandardDeviation --Both of which, having been normalised, should be zero. + ) . ( + normalise distribution :: [Double] -> [Double] + ) . take 10000 $ Math.Probability.generatePopulation distribution randomGen where + probability' = recip . succ $ abs probability --Semi-closed unit-interval (0, 1]. + distribution = Math.Probability.ShiftedGeometricDistribution probability' + + prop_shiftedGeometricDistribution' randomGen probability = Test.QuickCheck.label "prop_shiftedGeometricDistribution'" . all ( + >= (1 :: Double) + ) . take 10 $ Math.Probability.generatePopulation (Math.Probability.ShiftedGeometricDistribution probability') randomGen where + probability' = recip . succ $ abs probability --Semi-closed unit-interval (0, 1]. 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 @@ -25,6 +25,7 @@ run ) where +import qualified Control.Arrow import qualified Factory.Test.QuickCheck.ArithmeticGeometricMean import qualified Factory.Test.QuickCheck.Factorial import qualified Factory.Test.QuickCheck.Hyperoperation @@ -45,22 +46,25 @@ -- | Run the /quickChecks/-functions for modules supporting this feature. 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 "PerfectPower" >> Factory.Test.QuickCheck.PerfectPower.quickChecks - >> putStrLn "Pi" >> Factory.Test.QuickCheck.Pi.quickChecks - >> putStrLn "Polynomial" >> Factory.Test.QuickCheck.Polynomial.quickChecks - >> putStrLn "Power" >> Factory.Test.QuickCheck.Power.quickChecks - >> 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 - >> putStrLn "Statistics" >> Factory.Test.QuickCheck.Statistics.quickChecks - >> putStrLn "Summation" >> Factory.Test.QuickCheck.Summation.quickChecks +run = mapM_ ( + uncurry (>>) . Control.Arrow.first putStrLn + ) [ + ("ArithmeticGeometricMean", Factory.Test.QuickCheck.ArithmeticGeometricMean.quickChecks), + ("Factorial", Factory.Test.QuickCheck.Factorial.quickChecks), + ("Hyperoperation", Factory.Test.QuickCheck.Hyperoperation.quickChecks), + ("Interval", Factory.Test.QuickCheck.Interval.quickChecks), + ("MonicPolynomial", Factory.Test.QuickCheck.MonicPolynomial.quickChecks), + ("PerfectPower", Factory.Test.QuickCheck.PerfectPower.quickChecks), + ("Pi", Factory.Test.QuickCheck.Pi.quickChecks), + ("Polynomial", Factory.Test.QuickCheck.Polynomial.quickChecks), + ("Power", Factory.Test.QuickCheck.Power.quickChecks), + ("Primality", Factory.Test.QuickCheck.Primality.quickChecks), + ("PrimeFactorisation", Factory.Test.QuickCheck.PrimeFactorisation.quickChecks), + ("Primes", Factory.Test.QuickCheck.Primes.quickChecks), + ("Probability", Factory.Test.QuickCheck.Probability.quickChecks), + ("Radix", Factory.Test.QuickCheck.Radix.quickChecks), + ("SquareRoot", Factory.Test.QuickCheck.SquareRoot.quickChecks), + ("Statistics", Factory.Test.QuickCheck.Statistics.quickChecks), + ("Summation", Factory.Test.QuickCheck.Summation.quickChecks) + ] diff --git a/src/Factory/Test/QuickCheck/SquareRoot.hs b/src/Factory/Test/QuickCheck/SquareRoot.hs --- a/src/Factory/Test/QuickCheck/SquareRoot.hs +++ b/src/Factory/Test/QuickCheck/SquareRoot.hs @@ -53,7 +53,7 @@ coarbitrary = undefined --CAVEAT: stops warnings from ghc. #endif -type Testable = (Math.Implementations.SquareRoot.Algorithm, Math.Precision.DecimalDigits, Data.Ratio.Rational) -> Test.QuickCheck.Property +type Testable = (Math.Implementations.SquareRoot.Algorithm, Math.Precision.DecimalDigits, Rational) -> Test.QuickCheck.Property -- | Defines invariant properties. quickChecks :: IO () @@ -63,7 +63,7 @@ requiredDecimalDigits :: Math.Precision.DecimalDigits requiredDecimalDigits = succ $ decimalDigits `mod` 1024 - operand' :: Data.Ratio.Rational + operand' :: Rational operand' = abs operand prop_factorable (algorithm, decimalDigits, operand) = Test.QuickCheck.label "prop_factorable" . (<= 5) . ( @@ -78,14 +78,14 @@ requiredDecimalDigits :: Math.Precision.DecimalDigits requiredDecimalDigits = succ $ decimalDigits `mod` 1024 - operand' :: Data.Ratio.Rational + operand' :: Rational operand' = succ $ abs operand prop_perfectSquare (algorithm, decimalDigits, operand) = Test.QuickCheck.label "prop_perfectSquare" . Math.SquareRoot.isPrecise perfectSquare $ Math.SquareRoot.squareRoot algorithm requiredDecimalDigits perfectSquare where requiredDecimalDigits :: Math.Precision.DecimalDigits requiredDecimalDigits = succ $ decimalDigits `mod` 32768 - operand', perfectSquare :: Data.Ratio.Rational + operand', perfectSquare :: Rational operand' = (abs (Data.Ratio.numerator operand) `min` (2 ^ (32 :: Int))) % (abs (Data.Ratio.denominator operand) `min` (2 ^ (32 :: Int))) --Avoid floating-point rounding-errors in 'Math.SquareRoot.rSqrt'. perfectSquare = Math.Power.square operand' diff --git a/src/Factory/Test/QuickCheck/Statistics.hs b/src/Factory/Test/QuickCheck/Statistics.hs --- a/src/Factory/Test/QuickCheck/Statistics.hs +++ b/src/Factory/Test/QuickCheck/Statistics.hs @@ -29,7 +29,6 @@ import qualified Data.List import qualified Data.Map import qualified Data.Numbers.Primes -import qualified Data.Ratio import qualified Data.Set import qualified Factory.Math.Implementations.Factorial as Math.Implementations.Factorial import qualified Factory.Math.Power as Math.Power @@ -69,18 +68,18 @@ prop_nP1 i = Test.QuickCheck.label "prop_nP1" $ Math.Statistics.nPr n 1 == n where n = succ $ abs i - prop_zeroVariance, prop_zeroAverageAbsoluteDeviation :: Data.Ratio.Rational -> Test.QuickCheck.Property - prop_zeroVariance x = Test.QuickCheck.label "prop_zeroVariance" $ Math.Statistics.getVariance (replicate 32 x) == (0 :: Data.Ratio.Rational) - prop_zeroAverageAbsoluteDeviation x = Test.QuickCheck.label "zeroAverageAbsoluteDeviation" $ Math.Statistics.getAverageAbsoluteDeviation (replicate 32 x) == (0 :: Data.Ratio.Rational) + prop_zeroVariance, prop_zeroAverageAbsoluteDeviation :: Rational -> Test.QuickCheck.Property + prop_zeroVariance x = Test.QuickCheck.label "prop_zeroVariance" $ Math.Statistics.getVariance (replicate 32 x) == (0 :: Rational) + prop_zeroAverageAbsoluteDeviation x = Test.QuickCheck.label "zeroAverageAbsoluteDeviation" $ Math.Statistics.getAverageAbsoluteDeviation (replicate 32 x) == (0 :: Rational) prop_balance, prop_varianceRelocated, prop_varianceScaled, prop_varianceOrder, prop_equivalence, prop_varianceOfMap, prop_meanOfSet, prop_varianceOfArray :: [Integer] -> Test.QuickCheck.Property - prop_balance l = not (null l) ==> Test.QuickCheck.label "prop_balance" . (== 0) . abs . sum $ map (\i -> fromIntegral i - (Math.Statistics.getMean l :: Data.Ratio.Rational)) l - prop_varianceRelocated l = not (null l) ==> Test.QuickCheck.label "prop_varianceRelocated" $ (Math.Statistics.getVariance l :: Data.Ratio.Rational) == Math.Statistics.getVariance (map succ l) - prop_varianceScaled l = not (null l) ==> Test.QuickCheck.label "prop_varianceScaled" $ (4 * Math.Statistics.getVariance l :: Data.Ratio.Rational) == Math.Statistics.getVariance (map (* 2) l) - prop_varianceOrder l = not (null l) ==> Test.QuickCheck.label "prop_varianceOrder" $ Math.Statistics.getVariance l == (Math.Statistics.getVariance (reverse l) :: Data.Ratio.Rational) - prop_equivalence l = not (null l) ==> Test.QuickCheck.label "prop_equivalence" $ Math.Statistics.getVariance l == Math.Statistics.getMean (map Math.Power.square l) - Math.Power.square (Math.Statistics.getMean l :: Data.Ratio.Rational) - prop_varianceOfArray l = not (null l) ==> Test.QuickCheck.label "prop_varianceOfArray" $ Math.Statistics.getVariance (Data.Array.array (1, length l) $ zip [1 ..] l) == (Math.Statistics.getVariance l :: Data.Ratio.Rational) - prop_varianceOfMap l = not (null l) ==> Test.QuickCheck.label "prop_varianceOfMap" $ Math.Statistics.getVariance (Data.Map.fromList $ zip [0 :: Int ..] l) == (Math.Statistics.getVariance l :: Data.Ratio.Rational) - prop_meanOfSet l = not (null l') ==> Test.QuickCheck.label "prop_meanOfSet" $ Math.Statistics.getMean (Data.Set.fromList l') == (Math.Statistics.getMean l' :: Data.Ratio.Rational) where + prop_balance l = not (null l) ==> Test.QuickCheck.label "prop_balance" . (== 0) . abs . sum $ map (\i -> fromIntegral i - (Math.Statistics.getMean l :: Rational)) l + prop_varianceRelocated l = not (null l) ==> Test.QuickCheck.label "prop_varianceRelocated" $ (Math.Statistics.getVariance l :: Rational) == Math.Statistics.getVariance (map succ l) + prop_varianceScaled l = not (null l) ==> Test.QuickCheck.label "prop_varianceScaled" $ (4 * Math.Statistics.getVariance l :: Rational) == Math.Statistics.getVariance (map (* 2) l) + prop_varianceOrder l = not (null l) ==> Test.QuickCheck.label "prop_varianceOrder" $ Math.Statistics.getVariance l == (Math.Statistics.getVariance (reverse l) :: Rational) + prop_equivalence l = not (null l) ==> Test.QuickCheck.label "prop_equivalence" $ Math.Statistics.getVariance l == Math.Statistics.getMean (map Math.Power.square l) - Math.Power.square (Math.Statistics.getMean l :: Rational) + prop_varianceOfArray l = not (null l) ==> Test.QuickCheck.label "prop_varianceOfArray" $ Math.Statistics.getVariance (Data.Array.array (1, length l) $ zip [1 ..] l) == (Math.Statistics.getVariance l :: Rational) + prop_varianceOfMap l = not (null l) ==> Test.QuickCheck.label "prop_varianceOfMap" $ Math.Statistics.getVariance (Data.Map.fromList $ zip [0 :: Int ..] l) == (Math.Statistics.getVariance l :: Rational) + prop_meanOfSet l = not (null l') ==> Test.QuickCheck.label "prop_meanOfSet" $ Math.Statistics.getMean (Data.Set.fromList l') == (Math.Statistics.getMean l' :: Rational) where l' = Data.List.nub l diff --git a/src/Factory/Test/QuickCheck/Summation.hs b/src/Factory/Test/QuickCheck/Summation.hs --- a/src/Factory/Test/QuickCheck/Summation.hs +++ b/src/Factory/Test/QuickCheck/Summation.hs @@ -25,7 +25,6 @@ quickChecks ) where -import qualified Data.Ratio import qualified Factory.Math.Summation as Math.Summation import qualified Test.QuickCheck import Test.QuickCheck((==>)) @@ -33,7 +32,7 @@ -- | Defines invariant properties. quickChecks :: IO () quickChecks = Test.QuickCheck.quickCheck `mapM_` [prop_sum, prop_sumR] where - prop_sum, prop_sumR :: Int -> [Data.Ratio.Rational] -> Test.QuickCheck.Property + prop_sum, prop_sumR :: Int -> [Rational] -> Test.QuickCheck.Property prop_sum chunkSize l = not (null l) ==> Test.QuickCheck.label "prop_sum" $ Math.Summation.sum' chunkSize' l == sum l where chunkSize' = 2 + (chunkSize `mod` length l) diff --git a/src/Main.hs b/src/Main.hs --- a/src/Main.hs +++ b/src/Main.hs @@ -1,5 +1,5 @@ {- - Copyright (C) 2011 Dr. Alistair Ward + Copyright (C) 2011-2013 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 @@ -24,18 +24,10 @@ * Facilitates testing. -} -module Main( --- * Types --- ** Type-synonyms --- CommandLineAction, --- * Functions --- read', --- readCommandArg, - main -) where +module Main(main) where +import qualified Data.Map import qualified Data.List -import qualified Data.Ratio import qualified Data.Version import qualified Distribution.Package import qualified Distribution.Text @@ -46,6 +38,7 @@ import qualified Factory.Math.Implementations.PrimeFactorisation as Math.Implementations.PrimeFactorisation import qualified Factory.Math.Implementations.Primes.Algorithm as Math.Implementations.Primes.Algorithm import qualified Factory.Math.Implementations.SquareRoot as Math.Implementations.SquareRoot +import qualified Factory.Math.Probability as Math.Probability 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 @@ -62,6 +55,7 @@ import qualified System.Exit import qualified System.IO import qualified System.IO.Error +import qualified System.Random import qualified ToolShed.Defaultable -- Local convenience definitions. @@ -74,7 +68,7 @@ -- | On failure to parse the specified string, returns an explanatory error. read' :: Read a => String -> String -> a read' errorMessage s = case reads s of - [(x, _)] -> x + [(x, "")] -> x _ -> error $ errorMessage ++ show s -- | On failure to parse a command-line argument, returns an explanatory error. @@ -84,6 +78,8 @@ -- | Parses the command-line arguments, to determine 'Test.CommandOptions.CommandOptions'. main :: IO () main = do + System.IO.hClose System.IO.stdin --Nothing is read from standard input. + progName <- System.Environment.getProgName let @@ -111,14 +107,15 @@ 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 "" ["plotDiscreteDistribution"] (plotDiscreteDistribution `G.ReqArg` "(Int, Math.Probability.DiscreteDistribution)") "Plot the Probability Mass function for the specified discrete distribution.", G.Option "" ["primeFactorsPerformance"] (primeFactorsPerformance `G.ReqArg` "(Math.Implementations.PrimeFactorisation.Algorithm, Integer)") "Test the performance of 'Math.PrimeFactorisation.primeFactors'.", G.Option "" ["primeFactorsPerformanceGraph"] (primeFactorsPerformanceGraph `G.ReqArg` "(Math.Implementations.PrimeFactorisation.Algorithm, Int)") "Test the performance of 'Math.PrimeFactorisation.primeFactors', on the specified number of odd integers from the Fibonacci-sequence.", G.Option "" ["primesPerformance"] (primesPerformance `G.ReqArg` "(Math.Implementations.Primes.Algorithm.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 "" ["squareRootPerformance"] (squareRootPerformance `G.ReqArg` "(Math.Implementations.SquareRoot.Algorithm, Rational, DecimalDigits)") "Test the performance of 'Math.SquareRoot.squareRoot'.", + G.Option "" ["squareRootPerformanceGraph"] (squareRootPerformanceGraph `G.ReqArg` "(Math.Implementations.SquareRoot.Algorithm, Rational)") "Test the performance of 'Math.SquareRoot.squareRoot', with an exponentially increasing precision-requirement." ] where printVersion, printUsage, runQuickChecks :: IO Test.CommandOptions.CommandOptions - printVersion = System.IO.hPutStrLn System.IO.stderr (Distribution.Text.display packageIdentifier ++ "\n\nCopyright (C) 2011 " ++ author ++ ".\nThis program comes with ABSOLUTELY NO WARRANTY.\nThis is free software, and you are welcome to redistribute it under certain conditions.\n\nWritten by " ++ author ++ ".") >> System.Exit.exitWith System.Exit.ExitSuccess where + printVersion = System.IO.hPutStrLn System.IO.stderr (Distribution.Text.display packageIdentifier ++ "\n\nCopyright (C) 2011-2013 " ++ author ++ ".\nThis program comes with ABSOLUTELY NO WARRANTY.\nThis is free software, and you are welcome to redistribute it under certain conditions.\n\nWritten by " ++ author ++ ".") >> System.Exit.exitWith System.Exit.ExitSuccess where packageIdentifier :: Distribution.Package.PackageIdentifier packageIdentifier = Distribution.Package.PackageIdentifier { Distribution.Package.pkgName = Distribution.Package.PackageName progName, --CAVEAT: coincidentally. @@ -135,7 +132,7 @@ factorialPerformanceGraphControl :: Test.CommandOptions.CommandOptions -> IO Test.CommandOptions.CommandOptions factorialPerformanceGraphControl commandOptions = Test.Performance.Factorial.factorialPerformanceGraphControl (Test.CommandOptions.verbose commandOptions) >> System.Exit.exitWith (System.Exit.ExitFailure 1) - carmichaelNumbersPerformance, factorialPerformance, factorialPerformanceGraph, hyperoperationPerformance, hyperoperationPerformanceGraphRank, hyperoperationPerformanceGraphExponent, isPrimePerformance, isPrimePerformanceGraph, mersenneNumbersPerformance, piPerformance, piPerformanceGraph, primeFactorsPerformance, primesPerformance, squareRootPerformance, squareRootPerformanceGraph :: String -> CommandLineAction + carmichaelNumbersPerformance, factorialPerformance, factorialPerformanceGraph, hyperoperationPerformance, hyperoperationPerformanceGraphRank, hyperoperationPerformanceGraphExponent, isPrimePerformance, isPrimePerformanceGraph, mersenneNumbersPerformance, piPerformance, piPerformanceGraph, plotDiscreteDistribution, primeFactorsPerformance, primesPerformance, squareRootPerformance, squareRootPerformanceGraph :: String -> CommandLineAction carmichaelNumbersPerformance arg _ = Test.Performance.Primality.carmichaelNumbersPerformance algorithm i >>= print >> System.Exit.exitWith System.Exit.ExitSuccess where algorithm :: PrimalityAlgorithm @@ -189,6 +186,14 @@ factor :: Double (category, factor, maxDecimalDigits) = readCommandArg arg + plotDiscreteDistribution arg _ = let + distribution :: Math.Probability.DiscreteDistribution Double + (n, distribution) = readCommandArg arg + in do + System.Random.getStdGen >>= print . Data.Map.toList . Data.Map.map ((/ (fromIntegral n :: Double)) . fromInteger) . Data.Map.fromListWith (+) . (`zip` repeat 1) . (take n :: [Integer] -> [Integer]) . Math.Probability.generateDiscretePopulation distribution + + System.Exit.exitWith System.Exit.ExitSuccess + primeFactorsPerformance arg _ = Test.Performance.PrimeFactorisation.primeFactorsPerformance algorithm i >>= print >> System.Exit.exitWith System.Exit.ExitSuccess where algorithm :: Math.Implementations.PrimeFactorisation.Algorithm (algorithm, i) = readCommandArg arg @@ -219,12 +224,12 @@ squareRootPerformance arg _ = Test.Performance.SquareRoot.squareRootPerformance algorithm operand decimalDigits >>= print >> System.Exit.exitWith System.Exit.ExitSuccess where algorithm :: Math.Implementations.SquareRoot.Algorithm - operand :: Data.Ratio.Rational + operand :: Rational (algorithm, operand, decimalDigits) = readCommandArg arg squareRootPerformanceGraph arg _ = Test.Performance.SquareRoot.squareRootPerformanceGraph algorithm operand >> System.Exit.exitWith (System.Exit.ExitFailure 1) where algorithm :: Math.Implementations.SquareRoot.Algorithm - operand :: Data.Ratio.Rational + operand :: Rational (algorithm, operand) = readCommandArg arg args <- System.Environment.getArgs