diff --git a/changelog b/changelog --- a/changelog +++ b/changelog @@ -13,4 +13,11 @@ * 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 '(^)'. -0.0.0.3 + * 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. + diff --git a/factory.cabal b/factory.cabal --- a/factory.cabal +++ b/factory.cabal @@ -1,6 +1,6 @@ --Package-properties Name: factory -Version: 0.0.0.2 +Version: 0.1.0.0 Cabal-Version: >= 1.6 Copyright: (C) 2011 Dr. Alistair Ward License: GPL @@ -11,7 +11,7 @@ Build-Type: Simple Description: A library of number-theory functions, for; factorials, square-roots, Pi, primality-testing, prime-factorisation ... Category: Math, Number Theory -Tested-With: GHC == 6.12, GHC == 7.0 +Tested-With: GHC == 6.10, GHC == 6.12, GHC == 7.0 Homepage: http://functionalley.eu Maintainer: factory functionalley eu Bug-reports: factory functionalley eu @@ -30,8 +30,8 @@ hs-source-dirs: src Exposed-modules: - Factory.Data.Bounds Factory.Data.Exponential + Factory.Data.Interval Factory.Data.MonicPolynomial Factory.Data.Monomial Factory.Data.Polynomial @@ -43,11 +43,6 @@ Factory.Math.Factorial Factory.Math.Fibonacci Factory.Math.Implementations.Factorial - Factory.Math.Implementations.Primality - Factory.Math.Implementations.PrimeFactorisation - Factory.Math.Implementations.SquareRoot - Factory.Math.MultiplicativeOrder - Factory.Math.Pi Factory.Math.Implementations.Pi.AGM.Algorithm Factory.Math.Implementations.Pi.AGM.BrentSalamin Factory.Math.Implementations.Pi.BBP.Algorithm @@ -69,10 +64,16 @@ Factory.Math.Implementations.Pi.Spigot.RabinowitzWagon Factory.Math.Implementations.Pi.Spigot.Series Factory.Math.Implementations.Pi.Spigot.Spigot + Factory.Math.Implementations.Primality + Factory.Math.Implementations.PrimeFactorisation + Factory.Math.Implementations.SquareRoot + Factory.Math.MultiplicativeOrder + Factory.Math.Pi Factory.Math.Power Factory.Math.Precision Factory.Math.Primality Factory.Math.PrimeFactorisation + Factory.Math.Probability Factory.Math.Radix Factory.Math.SquareRoot Factory.Math.Statistics @@ -84,7 +85,8 @@ deepseq >= 1.1, containers, primes >= 0.1, - toolshed == 0.11.* + random, + toolshed >= 0.12 if flag(threaded) Build-depends: parallel >= 3.0 @@ -111,16 +113,17 @@ Factory.Test.Performance.SquareRoot Factory.Test.Performance.Statistics Factory.Test.QuickCheck.ArithmeticGeometricMean - Factory.Test.QuickCheck.Bounds Factory.Test.QuickCheck.Factorial + Factory.Test.QuickCheck.Interval Factory.Test.QuickCheck.MonicPolynomial Factory.Test.QuickCheck.Pi Factory.Test.QuickCheck.Polynomial Factory.Test.QuickCheck.Power Factory.Test.QuickCheck.Primality Factory.Test.QuickCheck.PrimeFactorisation - Factory.Test.QuickCheck.Radix + Factory.Test.QuickCheck.Probability Factory.Test.QuickCheck.QuickChecks + Factory.Test.QuickCheck.Radix Factory.Test.QuickCheck.SquareRoot Factory.Test.QuickCheck.Statistics Factory.Test.QuickCheck.Summation @@ -134,7 +137,10 @@ GHC-prof-options: -prof -auto-all -caf-all if flag(threaded) - GHC-options: -threaded -feager-blackholing + GHC-options: -threaded + + if impl(ghc >= 6.12) + GHC-options: -feager-blackholing if impl(ghc >= 7.0) GHC-options: -rtsopts diff --git a/makefile b/makefile --- a/makefile +++ b/makefile @@ -13,7 +13,7 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -.PHONY: all build clean configure copy haddock help hlint install prof sdist +.PHONY: all build check clean configure copy haddock help hlint install prof sdist all: install @@ -45,8 +45,12 @@ sdist: configure runhaskell Setup.hs $@ +check: sdist + cabal upload --check --verbose=3 dist/*.tar.gz; + clean: runhaskell Setup.hs $@ + find src -type f \( -name '*.o' -o -name '*.hi' \) -delete help: @grep '^[a-zA-Z].*:' makefile | sed -e 's/:.*//' diff --git a/src/Factory/Data/Bounds.hs b/src/Factory/Data/Bounds.hs deleted file mode 100644 --- a/src/Factory/Data/Bounds.hs +++ /dev/null @@ -1,170 +0,0 @@ -{-# LANGUAGE CPP #-} -{- - Copyright (C) 2011 Dr. Alistair Ward - - This program is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with this program. If not, see . --} -{- | - [@AUTHOR@] Dr. Alistair Ward - - [@DESCRIPTION@] - - * Describes a /bounded/ range of, typically integral, quantities. - - * Operations have been defined, on the list of /consecutive/ quantities delimited by these two bounds. - - * The point is that if the list is composed from /consecutive/ quantities, the intermediate values can be inferred, rather than physically represented. - - [@CAVEATS@] - - * The API was driven top-down by its caller's requirements, rather than a bottom-up attempt to provide a complete interface. - consequently there may be omissions from the view point of future callers. --} - -module Factory.Data.Bounds ( --- * Types --- ** Type-synonyms - Bounds, --- * Functions --- divideAndConquer, - elem', - length', - normalise, - product', - splitAt', - toList, --- ** Accessors - minBound', - maxBound' --- ** Predicates --- isReversed -) where - -import Control.Arrow((***)) -import qualified Data.Monoid -import qualified Data.Ratio - -#if MIN_VERSION_parallel(3,0,0) -import qualified Control.Parallel.Strategies -#endif - -#if MIN_VERSION_base(4,3,0) -import Data.Tuple(swap) -#else --- | Swap the components of a pair. -swap :: (a, b) -> (b, a) -swap (a, b) = (b, a) -#endif - --- | Defines a range of consecutive values, bracketed by /inclusive/ bounds. -type Bounds limit = (limit, limit) - --- | Accessor. -{-# INLINE minBound' #-} -minBound' :: Bounds a -> a -minBound' = fst - --- | Accessor. -{-# INLINE maxBound' #-} -maxBound' :: Bounds a -> a -maxBound' = snd - --- | 'True' if the specified value is within the inclusive 'Bounds'. -elem' :: Ord limit => limit -> Bounds limit -> Bool -elem' x = uncurry (&&) . ((<= x) *** (x <=)) - --- | 'True' if /minBound'/ exceeds /maxBound'/ extent. -isReversed :: Ord limit => Bounds limit -> Bool -isReversed = uncurry (>) - --- | Swap the limits where they were originally reversed, but otherwise do nothing. -normalise :: Ord limit => Bounds limit -> Bounds limit -normalise b - | isReversed b = swap b - | otherwise = b - --- | Bisect the bounds at the specified limit; which should be between the two existing limits. -splitAt' :: (Num limit, Ord limit) => limit -> Bounds limit -> (Bounds limit, Bounds limit) -splitAt' i bounds@(l, r) - | any ($ i) [(< l), (>= r)] = error $ "Factory.Data.Bounds.splitAt':\tunsuitable index=" ++ show i ++ " for bounds=" ++ show bounds ++ "." - | otherwise = ((l, i), (i + 1, r)) - --- | The length of 'toList'. -{-# INLINE length' #-} -length' :: (Num limit, Ord limit) => Bounds limit -> limit -length' (l, r) = r + 1 - l - --- | Converts 'Bounds' to a list by enumerating the values. -{-# INLINE toList #-} -toList :: Enum limit => Bounds limit -> [limit] -toList = uncurry enumFromTo - -{- | - * Reduces 'Bounds' to a single integral value encapsulated in a 'Data.Monoid.Monoid', - using a /divide-and-conquer/ strategy, - bisecting the /bounds/ and recursively evaluating each part; . - - * By choosing a 'ratio' other than @(1 % 2)@, the bisection can be made asymmetrical. - The specified ratio represents the length of the left-hand portion over the original list-length; - eg. @(1 % 3)@ results in the first part, half the length of the second. - - * This process of recursive bisection, is terminated beneath the specified minimum length, - after which the 'Bounds' are expanded into the corresponding list, and the /monoid/'s binary operator is directly /folded/ over it. - - * One can view this as a , - in which 'Bounds' is exploded into a binary tree-structure - (each leaf of which contains a list of up to 'minLength' integers, and each node of which contains an associative binary operator), - and then collapsed to a scalar, by application of the operators. --} -divideAndConquer :: (Integral i, Data.Monoid.Monoid monoid) - => (i -> monoid) -- ^ The monoid's constructor. - -> Data.Ratio.Ratio i -- ^ The ratio of the original span, at which to bisect the 'Bounds'. - -> i -- ^ For efficiency, the bounds will not be bisected, when it's length has been reduced to this value. - -> Bounds i - -> monoid -- ^ The resulting scalar. -divideAndConquer monoidConstructor ratio minLength - | any ($ ratio) [ - (< 0), - (>= 1) - ] = error $ "Factory.Data.Bounds.divideAndConquer:\tunsuitable ratio='" ++ show ratio ++ "'." - | minLength < 1 = error $ "Factory.Data.Bounds.divideAndConquer:\tunsuitable minLength=" ++ show minLength ++ "." - | otherwise = slave - where - slave bounds@(l, r) - | length' bounds <= minLength = Data.Monoid.mconcat . map monoidConstructor $ toList bounds --Fold the monoid's binary operator over the delimited list. - | otherwise = uncurry Data.Monoid.mappend . -#if MIN_VERSION_parallel(3,0,0) - Control.Parallel.Strategies.withStrategy ( - Control.Parallel.Strategies.parTuple2 Control.Parallel.Strategies.rseq Control.Parallel.Strategies.rseq - ) . -#endif - (slave *** slave) $ splitAt' ( - l + (r - l) * Data.Ratio.numerator ratio `div` Data.Ratio.denominator ratio --Use the ratio to generate the split-index. - ) bounds --Apply the monoid's binary operator to the two operands resulting from bisection. - -{- | - * Multiplies the consecutive sequence of integers within 'Bounds'. - - * Since the result can be large, 'divideAndConquer' is used to form operands of a similar order of magnitude, - thus improving the efficiency of the big-number multiplication. --} -product' :: Integral i - => Data.Ratio.Ratio i -- ^ The ratio at which to bisect the 'Bounds'. - -> i -- ^ For efficiency, the bounds will not be bisected, when it's length has been reduced to this value. - -> Bounds i - -> i -- ^ The resulting product. -product' ratio minLength bounds - | elem' 0 bounds = 0 - | otherwise = Data.Monoid.getProduct $ divideAndConquer Data.Monoid.Product ratio minLength bounds - diff --git a/src/Factory/Data/Interval.hs b/src/Factory/Data/Interval.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Data/Interval.hs @@ -0,0 +1,195 @@ +{-# LANGUAGE CPP #-} +{- + Copyright (C) 2011 Dr. Alistair Ward + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * Describes a bounded set of, typically integral, quantities. + + * Operations have been defined, on the list of /consecutive/ quantities delimited by these endpoints. + + * The point is that if the list is composed from /consecutive/ quantities, the intermediate values can be inferred, rather than physically represented. + + [@CAVEATS@] + + * The API was driven top-down by its caller's requirements, rather than a bottom-up attempt to provide a complete interface. + consequently there may be omissions from the view point of future callers. + + * Thought similar to the mathematical concept of an /interval/, the latter technically relates to /real/ numbers; . + + * No account has been made for /semi-closed/ or /open/ intervals. +-} + +module Factory.Data.Interval( +-- * Types +-- ** Type-synonyms + Interval, +-- * Constants + closedUnitInterval, +-- * Functions +-- divideAndConquer, + elem', + getLength, + normalise, + product', + shift, + splitAt', + toList, +-- ** Accessors + getMinBound, + getMaxBound, +-- ** Constructors + precisely, +-- ** Predicates + isReversed +) where + +import Control.Arrow((***), (&&&)) +import qualified Data.Monoid +import qualified Data.Ratio +import qualified ToolShed.Pair as Pair + +#if MIN_VERSION_parallel(3,0,0) +import qualified Control.Parallel.Strategies +#endif + +#if MIN_VERSION_base(4,3,0) +import Data.Tuple(swap) +#else +-- | Swap the components of a pair. +swap :: (a, b) -> (b, a) +swap (a, b) = (b, a) +#endif + +-- | Defines a closed (inclusive) interval of consecutive values. +type Interval endPoint = (endPoint, endPoint) + +-- | Accessor. +{-# INLINE getMinBound #-} +getMinBound :: Interval endPoint -> endPoint +getMinBound = fst + +-- | Accessor. +{-# INLINE getMaxBound #-} +getMaxBound :: Interval endPoint -> endPoint +getMaxBound = snd + +-- | Construct the interval from a single value. +precisely :: endPoint -> Interval endPoint +precisely = id &&& id + +-- | Construct the /closed unit-interval/; . +closedUnitInterval :: Num n => Interval n +closedUnitInterval = (0, 1) + +-- | Shift of both /end-points/ of the /interval/ by the specified amount. +shift :: Num endPoint + => endPoint -- ^ The magnitude of the require shift. + -> Interval endPoint -- ^ The interval to be shifted. + -> Interval endPoint +shift i = Pair.mirror (+ i) + +-- | 'True' if the specified value is within the inclusive bounds of the /interval/. +elem' :: Ord endPoint => endPoint -> Interval endPoint -> Bool +elem' x = Pair.both . ((<= x) *** (x <=)) + +-- | 'True' if 'getMinBound' exceeds 'getMaxBound' extent. +isReversed :: Ord endPoint => Interval endPoint -> Bool +isReversed = uncurry (>) + +-- | Swap the /end-points/ where they were originally reversed, but otherwise do nothing. +normalise :: Ord endPoint => Interval endPoint -> Interval endPoint +normalise b + | isReversed b = swap b + | otherwise = b + +-- | Bisect the /interval/ at the specified /end-point/; which should be between the two existing /end-points/. +splitAt' :: (Num endPoint, Ord endPoint) => endPoint -> Interval endPoint -> (Interval endPoint, Interval endPoint) +splitAt' i interval@(l, r) + | any ($ i) [(< l), (>= r)] = error $ "Factory.Data.Interval.splitAt':\tunsuitable index=" ++ show i ++ " for interval=" ++ show interval ++ "." + | otherwise = ((l, i), (i + 1, r)) + +-- | The length of 'toList'. +{-# INLINE getLength #-} +getLength :: (Num endPoint, Ord endPoint) => Interval endPoint -> endPoint +getLength (l, r) = r + 1 - l + +-- | Converts 'Interval' to a list by enumerating the values. +{-# INLINE toList #-} +toList :: Enum endPoint => Interval endPoint -> [endPoint] +toList = uncurry enumFromTo + +{- | + * Reduces 'Interval' to a single integral value encapsulated in a 'Data.Monoid.Monoid', + using a /divide-and-conquer/ strategy, + bisecting the /interval/ and recursively evaluating each part; . + + * By choosing a 'ratio' other than @(1 % 2)@, the bisection can be made asymmetrical. + The specified ratio represents the length of the left-hand portion over the original list-length; + eg. @(1 % 3)@ results in the first part, half the length of the second. + + * This process of recursive bisection, is terminated beneath the specified minimum length, + after which the 'Interval' are expanded into the corresponding list, and the /monoid/'s binary operator is directly /folded/ over it. + + * One can view this as a , + in which 'Interval' is exploded into a binary tree-structure + (each leaf of which contains a list of up to 'minLength' integers, and each node of which contains an associative binary operator), + and then collapsed to a scalar, by application of the operators. +-} +divideAndConquer :: (Integral i, Data.Monoid.Monoid monoid) + => (i -> monoid) -- ^ The monoid's constructor. + -> Data.Ratio.Ratio i -- ^ The ratio of the original span, at which to bisect the 'Interval'. + -> i -- ^ For efficiency, the /interval/ will not be bisected, when it's length has been reduced to this value. + -> Interval i + -> monoid -- ^ The resulting scalar. +divideAndConquer monoidConstructor ratio minLength + | any ($ ratio) [ + (< 0), + (>= 1) + ] = error $ "Factory.Data.Interval.divideAndConquer:\tunsuitable ratio='" ++ show ratio ++ "'." + | minLength < 1 = error $ "Factory.Data.Interval.divideAndConquer:\tunsuitable minLength=" ++ show minLength ++ "." + | otherwise = slave + where + slave interval@(l, r) + | getLength interval <= minLength = Data.Monoid.mconcat . map monoidConstructor $ toList interval --Fold the monoid's binary operator over the delimited list. + | otherwise = uncurry Data.Monoid.mappend . +#if MIN_VERSION_parallel(3,0,0) + Control.Parallel.Strategies.withStrategy ( + Control.Parallel.Strategies.parTuple2 Control.Parallel.Strategies.rseq Control.Parallel.Strategies.rseq + ) . +#endif + Pair.mirror slave $ splitAt' ( + l + (r - l) * Data.Ratio.numerator ratio `div` Data.Ratio.denominator ratio --Use the ratio to generate the split-index. + ) interval --Apply the monoid's binary operator to the two operands resulting from bisection. + +{- | + * Multiplies the consecutive sequence of integers within 'Interval'. + + * Since the result can be large, 'divideAndConquer' is used to form operands of a similar order of magnitude, + thus improving the efficiency of the big-number multiplication. +-} +product' :: Integral i + => Data.Ratio.Ratio i -- ^ The ratio at which to bisect the 'Interval'. + -> i -- ^ For efficiency, the /interval/ will not be bisected, when it's length has been reduced to this value. + -> Interval i + -> i -- ^ The resulting product. +product' ratio minLength interval + | elem' 0 interval = 0 + | otherwise = Data.Monoid.getProduct $ divideAndConquer Data.Monoid.Product ratio minLength interval + diff --git a/src/Factory/Data/MonicPolynomial.hs b/src/Factory/Data/MonicPolynomial.hs --- a/src/Factory/Data/MonicPolynomial.hs +++ b/src/Factory/Data/MonicPolynomial.hs @@ -32,7 +32,6 @@ mkMonicPolynomial ) where -import Control.Arrow((***)) import qualified Control.Arrow import qualified Factory.Data.Monomial as Data.Monomial import Factory.Data.Polynomial((*=)) @@ -40,6 +39,7 @@ import qualified Factory.Data.QuotientRing as Data.QuotientRing import Factory.Data.Ring((=*=), (=+=), (=-=)) import qualified Factory.Data.Ring as Data.Ring +import qualified ToolShed.Pair as Pair -- | A type of 'Data.Polynomial.Polynomial', in which the /leading term/ is required to have a /coefficient/ of one. newtype MonicPolynomial c e = MkMonicPolynomial { @@ -72,7 +72,7 @@ -- Since the /leading term/ of the /denominator/ is one, the /coefficient/ isn't required to implement 'Fractional'. instance (Num c, Num e, Ord e) => Data.QuotientRing.QuotientRing (MonicPolynomial c e) where - MkMonicPolynomial polynomialN `quotRem'` MkMonicPolynomial polynomialD = (MkMonicPolynomial *** MkMonicPolynomial) $ longDivide polynomialN where + MkMonicPolynomial polynomialN `quotRem'` MkMonicPolynomial polynomialD = Pair.mirror MkMonicPolynomial $ longDivide polynomialN where -- longDivide :: (Num c, Num e, Ord e) => Polynomial c e -> (Polynomial c e, Polynomial c e) longDivide numerator | Data.Polynomial.isZero numerator || Data.Monomial.getExponent quotient < 0 = (Data.Polynomial.zero, numerator) 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 @@ -25,7 +25,7 @@ * Implementations of this strategy have been provided for /addition/ and /multiplication/, though other associative binary operators, like 'gcd' or 'lcm' could also be used. - * Where the contents of the list are consecutive, a more efficient implementation is available in /Factory.Data.Bounds/. + * Where the contents of the list are consecutive, a more efficient implementation is available in /Factory.Data.Interval/. -} module Factory.Math.DivideAndConquer( 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 @@ -45,7 +45,7 @@ ) where import qualified Data.Numbers.Primes -import qualified Factory.Data.Bounds as Data.Bounds +import qualified Factory.Data.Interval as Data.Interval import qualified Factory.Data.PrimeFactors as Data.PrimeFactors import qualified Factory.Math.Factorial as Math.Factorial import qualified ToolShed.Defaultable as Defaultable @@ -54,7 +54,7 @@ -- | The algorithms by which /factorial/ has been implemented. data Algorithm = - Bisection -- ^ The integers from which the /factorial/ is composed, are multiplied using @Data.Bounds.product'@. + Bisection -- ^ The integers from which the /factorial/ is composed, are multiplied using @Data.Interval.product'@. | PrimeFactorisation -- ^ The /prime factors/ of the /factorial/ are extracted, then raised to the appropriate power, before multiplication. deriving (Eq, Read, Show) @@ -104,7 +104,7 @@ -> i -- ^ The result. risingFactorial _ 0 = 1 risingFactorial 0 _ = 0 -risingFactorial x n = Data.Bounds.product' (recip 2) 64 $ Data.Bounds.normalise (x, (x + n) - 1) +risingFactorial x n = Data.Interval.product' (recip 2) 64 $ Data.Interval.normalise (x, (x + n) - 1) -- | Returns the /falling factorial/; fallingFactorial :: Integral i @@ -113,7 +113,7 @@ -> i -- ^ The result. fallingFactorial _ 0 = 1 fallingFactorial 0 _ = 0 -fallingFactorial x n = Data.Bounds.product' (recip 2) 64 $ Data.Bounds.normalise (x, (x - n) + 1) +fallingFactorial x n = Data.Interval.product' (recip 2) 64 $ Data.Interval.normalise (x, (x - n) + 1) {- | * Returns the ratio of two factorials. @@ -134,5 +134,5 @@ | denominator <= 1 = fromIntegral $ Math.Factorial.factorial (Defaultable.defaultValue :: Algorithm) numerator | numerator == denominator = 1 | numerator < denominator = recip $ denominator !/! numerator --Recurse. - | otherwise = fromIntegral $ Data.Bounds.product' (recip 2) 64 (succ denominator, numerator) + | otherwise = fromIntegral $ Data.Interval.product' (recip 2) 64 (succ denominator, numerator) 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 @@ -17,7 +17,7 @@ {- | [@AUTHOR@] Dr. Alistair Ward - [@DESCRIPTION@] Defines the /Ramanujan/ series for /Pi/; . + [@DESCRIPTION@] Defines the /Ramanujan/ series for /Pi/; . -} module Factory.Math.Implementations.Pi.Ramanujan.Classic( diff --git a/src/Factory/Math/Implementations/Pi/Spigot/Spigot.hs b/src/Factory/Math/Implementations/Pi/Spigot/Spigot.hs --- a/src/Factory/Math/Implementations/Pi/Spigot/Spigot.hs +++ b/src/Factory/Math/Implementations/Pi/Spigot/Spigot.hs @@ -45,7 +45,6 @@ -- getRemainder ) where -import Data.Ratio((%)) import qualified Control.Arrow import qualified Data.Char import qualified Data.Ratio @@ -133,7 +132,7 @@ * Call 'processColumns'. -} mkRow :: Math.Implementations.Pi.Spigot.Series.Series I -> PreDigits -> Coefficients -> Pi -mkRow series preDigits = processColumns series preDigits . zip (1 % decimal : Math.Implementations.Pi.Spigot.Series.bases series) . map (* decimal) +mkRow series preDigits = processColumns series preDigits . zip (recip (fromIntegral decimal) : Math.Implementations.Pi.Spigot.Series.bases series) . map (* decimal) {- | * Initialises a /spigot/-table with the row of 'Math.Implementations.Pi.Spigot.Series.coefficients'. 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 @@ -39,7 +39,7 @@ -- factoriseByTrialDivision ) where -import Control.Arrow((&&&), (***)) +import Control.Arrow((&&&)) import qualified Control.Arrow import qualified Control.DeepSeq import qualified Data.Maybe @@ -50,6 +50,7 @@ import qualified Factory.Math.Power as Math.Power import qualified Factory.Math.PrimeFactorisation as Math.PrimeFactorisation import qualified ToolShed.Defaultable as Defaultable +import qualified ToolShed.Pair as Pair #if MIN_VERSION_parallel(3,0,0) import qualified Control.Parallel.Strategies @@ -108,9 +109,7 @@ Control.Parallel.Strategies.parTuple2 Control.Parallel.Strategies.rdeepseq Control.Parallel.Strategies.rdeepseq --CAVEAT: unproductive on the size of integers tested so far. ) . #endif - ( - factoriseByFermatsMethod *** factoriseByFermatsMethod --Divide and conquer. - ) $ head factors + Pair.mirror factoriseByFermatsMethod $ head factors where -- maybeSquareNumber :: Integral i => Maybe i maybeSquareNumber = Math.Power.maybeSquareNumber i diff --git a/src/Factory/Math/Probability.hs b/src/Factory/Math/Probability.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Probability.hs @@ -0,0 +1,181 @@ +{- + 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@] Miscellaneous functions for probability-distributions. +-} + +module Factory.Math.Probability( +-- * Types +-- ** Data-types + ContinuousDistribution(..), + DiscreteDistribution(..), +-- * Functions + boxMullerTransform, +-- minPositiveFloat, +-- reProfile, + generateStandardizedNormalDistribution, + generateContinuousPopulation, + generatePoissonDistribution, + generateDiscretePopulation +) where + +import qualified Control.Arrow +import Control.Arrow((***), (&&&)) +import qualified Factory.Data.Interval as Data.Interval +import qualified System.Random +import qualified ToolShed.ListPlus as ListPlus +import qualified ToolShed.Pair as Pair +import qualified ToolShed.SelfValidate as 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/; . + deriving (Eq, Read, Show) + +instance (Num a, Ord a) => SelfValidate.SelfValidator (ContinuousDistribution a) where + isValid (UniformDistribution interval) = not $ Data.Interval.isReversed interval + isValid (NormalDistribution _ v) = v >= 0 + +-- | Describes a /discrete probability-distribution/; . +data DiscreteDistribution f = PoissonDistribution f deriving (Eq, Read, Show) + +instance (Num f, Ord f) => SelfValidate.SelfValidator (DiscreteDistribution f) where + isValid (PoissonDistribution lambda) = lambda >= 0 + +{- | + * 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) + => (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 . Pair.both $ 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 + where + inSemiClosedUnitInterval :: (Num n, Ord n) => n -> Bool + inSemiClosedUnitInterval = Pair.both . ((> 0) &&& (<= 1)) + + polarToCartesianTransform :: Floating f => (f, f) -> (f, f) + 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. + + * , . +-} +generateStandardizedNormalDistribution :: (System.Random.RandomGen g, RealFloat f, System.Random.Random f) => g -> [f] +generateStandardizedNormalDistribution = ListPlus.linearise . uncurry (zipWith $ curry boxMullerTransform) . Pair.mirror ( + 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. + + * 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. +-} +generateContinuousPopulation :: (RealFloat f, System.Random.Random f, System.Random.RandomGen g) + => Int -- ^ number of items. + -> ContinuousDistribution f + -> g -- ^ A generator of /uniformly distributed/ random numbers. + -> [f] +generateContinuousPopulation 0 _ _ = [] +generateContinuousPopulation populationSize probabilityDistribution randomGen + | populationSize < 0 = error $ "Factory.Math.Probability.generateDiscretePopulation:\tinvalid population-size=" ++ show populationSize + | not $ SelfValidate.isValid probabilityDistribution = error $ "Factory.Math.Probability.generateContinuousPopulation:\tinvalid; '" ++ show probabilityDistribution ++ "'" + | otherwise = take populationSize $ ( + case probabilityDistribution of + UniformDistribution interval -> System.Random.randomRs interval + NormalDistribution requiredMean requiredVariance -> reProfile requiredMean (sqrt requiredVariance) . generateStandardizedNormalDistribution + ) 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). + + * . + + * 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/. +-} +generatePoissonDistribution :: ( + RealFloat lambda, + System.Random.Random lambda, + System.Random.RandomGen g, + Integral events + ) + => lambda -- ^ Defines the required approximate value of both /mean/ and /variance/. + -> g + -> [events] +generatePoissonDistribution lambda + | lambda < 0 = error $ "Factory.Math.Probability.generatePoissonDistribution:\tinvalid lambda=" ++ 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 + | otherwise = generator + where + generator = uncurry (:) . ( + fst . head . dropWhile ( + (> exp (negate lambda)) . snd --CAVEAT: underflows if lambda > (103 :: Float, 745 :: Double). + ) . scanl ( + \accumulator random -> succ *** (* random) $ accumulator + ) (negate 1, 1) . System.Random.randomRs (0, 1) *** generator {-recurse-} + ) . System.Random.split + +-- | Generates a random sample-population, with the specified discrete probability-distribution. +generateDiscretePopulation :: ( + Ord f, + RealFloat f, + System.Random.Random f, + System.Random.RandomGen g, + Integral events + ) + => Int -- ^ number of items. + -> DiscreteDistribution f + -> g -- ^ 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 + | not $ SelfValidate.isValid probabilityDistribution = error $ "Factory.Math.Probability.generateDiscretePopulation:\tinvalid; '" ++ show probabilityDistribution ++ "'" + | otherwise = take populationSize $ ( + case probabilityDistribution of + PoissonDistribution lambda -> generatePoissonDistribution lambda + ) randomGen + 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 @@ -22,7 +22,12 @@ module Factory.Math.Statistics( -- * Functions - mean, + getMean, +-- getDispersionFromMean, + getVariance, + getStandardDeviation, + getAverageAbsoluteDeviation, + getCoefficientOfVariance, nCr, nPr ) where @@ -30,21 +35,64 @@ import Control.Arrow((***)) import Control.Parallel(par, pseq) import qualified Data.List ---import qualified Factory.Data.PrimeFactors as Data.PrimeFactors ---import Factory.Data.PrimeFactors((>/<), (>*<)) +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 --- | Determines the of the supplied numbers. -mean :: (Real r, Fractional f) => [r] -> f -mean [] = error "Factory.Math.Statistics.mean:\tundefined result for specified null-list" -mean l = uncurry (/) . (realToFrac *** fromIntegral) $ foldr (\s -> (+ s) *** succ) (0, 0 :: Int) l +{- | + * Determines the /mean/ of the specified list of numbers; . + * Should the caller define the result-type as 'Data.Ratio.Rational', then it will be free from rounding-errors. +-} +getMean :: (Real r, Fractional result) => [r] -> result +getMean [] = error "Factory.Math.Statistics.getMean:\tundefined result for null-list." +getMean [x] = realToFrac x --Not necessary, but a shortcut for this special case. +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; . + + * 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 + +{- | + * Determines the exact /variance/ of the specified list of numbers; . + + * Should the caller define the result-type as 'Data.Ratio.Rational', then it will be free from rounding-errors. +-} +getVariance :: (Real r, Fractional result) => [r] -> result +getVariance = getDispersionFromMean Math.Power.square + +-- | Determines the /standard-deviation/ of the specified list of numbers; . +getStandardDeviation :: (Real r, Floating result) => [r] -> result +getStandardDeviation = sqrt . getVariance + +{- | + * Determines the /average absolute deviation/ of the specified list of numbers; . + + * Should the caller define the result-type as 'Data.Ratio.Rational', then it will be free from rounding-errors. +-} +getAverageAbsoluteDeviation :: (Real r, Fractional result) => [r] -> result +getAverageAbsoluteDeviation = getDispersionFromMean abs + +-- | Determines the /coefficient-of-variance/ of the specified list of numbers; . +getCoefficientOfVariance :: (Real r, Floating result) => [r] -> result +getCoefficientOfVariance l + | mean == 0 = error "Factory.Math.Statistics.getCoefficientOfVariance:\tundefined if mean is zero." + | otherwise = getStandardDeviation l / abs mean + where + mean = getMean l + -- | The number of unordered combinations of /r/ objects taken from /n/; . nCr :: (Math.Factorial.Algorithm factorialAlgorithm, Integral i) => factorialAlgorithm -> i -- ^ The total number of items from which to select. - -> i -- ^ The number of iterms in a sample. + -> i -- ^ The number of items in a sample. -> i -- ^ The number of combinations. nCr _ 0 _ = 1 nCr _ _ 0 = 1 @@ -52,12 +100,6 @@ | n < 0 = error $ "Factory.Math.Statistics.nCr:\tinvalid n; " ++ show n | r < 0 = error $ "Factory.Math.Statistics.nCr:\tinvalid r; " ++ show r | n < r = 0 -{- - | otherwise = uncurry div $ product' *** product' $ Math.Implementations.Factorial.primeFactors n >/< ( - Math.Implementations.Factorial.primeFactors r >*< Math.Implementations.Factorial.primeFactors (n - r) - ) where - product' = Data.PrimeFactors.product' (recip 2) 10 --} | otherwise = numerator `par` (denominator `pseq` numerator `div` denominator) where [smaller, bigger] = Data.List.sort [r, n - r] diff --git a/src/Factory/Test/QuickCheck/Bounds.hs b/src/Factory/Test/QuickCheck/Bounds.hs deleted file mode 100644 --- a/src/Factory/Test/QuickCheck/Bounds.hs +++ /dev/null @@ -1,43 +0,0 @@ -{- - Copyright (C) 2011 Dr. Alistair Ward - - This program is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with this program. If not, see . --} -{- | - [@AUTHOR@] Dr. Alistair Ward - - [@DESCRIPTION@] Defines /QuickCheck/-properties for "Data.Bounds". --} - -module Factory.Test.QuickCheck.Bounds( --- * Functions - quickChecks -) where - -import qualified Data.Ratio -import qualified Factory.Data.Bounds as Data.Bounds -import qualified Test.QuickCheck - --- | Defines invariant properties. -quickChecks :: IO () -quickChecks = Test.QuickCheck.quickCheckWith Test.QuickCheck.stdArgs {Test.QuickCheck.maxSuccess = 1000} prop_product where - prop_product :: Data.Ratio.Ratio Integer -> Integer -> Data.Bounds.Bounds Integer -> Test.QuickCheck.Property - prop_product ratio minLength bounds = Test.QuickCheck.label "prop_product" $ Data.Bounds.product' ratio' minLength' bounds' == product (Data.Bounds.toList bounds') where - bounds' = Data.Bounds.normalise bounds - minLength' = 1 + minLength `mod` 1000 - ratio' = if r > 1 - then recip r - else r - where - r = abs ratio diff --git a/src/Factory/Test/QuickCheck/Interval.hs b/src/Factory/Test/QuickCheck/Interval.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/QuickCheck/Interval.hs @@ -0,0 +1,43 @@ +{- + 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 "Data.Interval". +-} + +module Factory.Test.QuickCheck.Interval( +-- * Functions + quickChecks +) where + +import qualified Data.Ratio +import qualified Factory.Data.Interval as Data.Interval +import qualified Test.QuickCheck + +-- | Defines invariant properties. +quickChecks :: IO () +quickChecks = Test.QuickCheck.quickCheckWith Test.QuickCheck.stdArgs {Test.QuickCheck.maxSuccess = 1000} prop_product where + prop_product :: Data.Ratio.Ratio Integer -> Integer -> Data.Interval.Interval Integer -> Test.QuickCheck.Property + prop_product ratio minLength interval = Test.QuickCheck.label "prop_product" $ Data.Interval.product' ratio' minLength' interval' == product (Data.Interval.toList interval') where + interval' = Data.Interval.normalise interval + minLength' = 1 + minLength `mod` 1000 + ratio' = if r > 1 + then recip r + else r + where + r = abs ratio diff --git a/src/Factory/Test/QuickCheck/Probability.hs b/src/Factory/Test/QuickCheck/Probability.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/QuickCheck/Probability.hs @@ -0,0 +1,69 @@ +{- + 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.Probability". +-} + +module Factory.Test.QuickCheck.Probability( +-- * Functions + quickChecks +) where + +import Control.Arrow((&&&)) +import qualified Factory.Math.Probability as Math.Probability +import qualified Factory.Math.Statistics as Math.Statistics +import Factory.Test.QuickCheck.Factorial() +import qualified ToolShed.Pair as Pair +import qualified System.Random +import qualified Test.QuickCheck +import Test.QuickCheck((==>)) + +-- | 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 g => g -> (Double, Double) -> Test.QuickCheck.Property + prop_normalDistribution randomGen (mean, variance) = variance' /= 0 ==> Test.QuickCheck.label "prop_normalDistribution" . Pair.both . Pair.mirror ( + (< (0.05 :: Double)) . abs --Tolerance. + ) . ( + 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 + + prop_poissonDistribution :: System.Random.RandomGen g => g -> Int -> Test.QuickCheck.Property + prop_poissonDistribution randomGen lambda = lambda' /= 0 ==> Test.QuickCheck.label "prop_poissonDistribution" . Pair.both . Pair.mirror ( + (< (0.1 :: Double)) . abs --Tolerance. + ) . ( + 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 $ lambda `mod` 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 @@ -1,4 +1,3 @@ -{-# LANGUAGE CPP #-} {- Copyright (C) 2011 Dr. Alistair Ward @@ -27,14 +26,15 @@ ) where import qualified Factory.Test.QuickCheck.ArithmeticGeometricMean -import qualified Factory.Test.QuickCheck.Bounds import qualified Factory.Test.QuickCheck.Factorial +import qualified Factory.Test.QuickCheck.Interval import qualified Factory.Test.QuickCheck.MonicPolynomial import qualified Factory.Test.QuickCheck.Pi import qualified Factory.Test.QuickCheck.Polynomial import qualified Factory.Test.QuickCheck.Power import qualified Factory.Test.QuickCheck.Primality import qualified Factory.Test.QuickCheck.PrimeFactorisation +import qualified Factory.Test.QuickCheck.Probability import qualified Factory.Test.QuickCheck.Radix import qualified Factory.Test.QuickCheck.SquareRoot import qualified Factory.Test.QuickCheck.Statistics @@ -43,14 +43,15 @@ -- | Run the /quickChecks/-functions for modules supporting this feature. run :: IO () run = putStrLn "ArithmeticGeometricMean" >> Factory.Test.QuickCheck.ArithmeticGeometricMean.quickChecks - >> putStrLn "Bounds" >> Factory.Test.QuickCheck.Bounds.quickChecks >> putStrLn "Factorial" >> Factory.Test.QuickCheck.Factorial.quickChecks + >> putStrLn "Interval" >> Factory.Test.QuickCheck.Interval.quickChecks >> putStrLn "MonicPolynomial" >> Factory.Test.QuickCheck.MonicPolynomial.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 "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 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 @@ -27,7 +27,9 @@ import qualified Data.List import qualified Data.Numbers.Primes +import qualified Data.Ratio import qualified Factory.Math.Implementations.Factorial as Math.Implementations.Factorial +import qualified Factory.Math.Power as Math.Power import qualified Factory.Math.Statistics as Math.Statistics import Factory.Test.QuickCheck.Factorial() import qualified Test.QuickCheck @@ -35,34 +37,43 @@ -- | Defines invariant properties. quickChecks :: IO () -quickChecks = - Test.QuickCheck.quickCheck `mapM_` [prop_nC0, prop_nC1, prop_sum] +quickChecks = Test.QuickCheck.quickCheck `mapM_` [prop_nC0, prop_nC1, prop_sum] >> Test.QuickCheck.quickCheck `mapM_` [prop_symmetry, prop_prime] >> Test.QuickCheck.quickCheck `mapM_` [prop_nP0, prop_nP1] - >> Test.QuickCheck.quickCheck prop_balance where - prop_nC0, prop_nC1, prop_sum :: Math.Implementations.Factorial.Algorithm -> Integer -> Test.QuickCheck.Property - prop_nC0 algorithm n = Test.QuickCheck.label "prop_nC0" $ Math.Statistics.nCr algorithm (abs n) 0 == 1 + >> Test.QuickCheck.quickCheck `mapM_` [prop_zeroVariance, prop_zeroAverageAbsoluteDeviation] + >> Test.QuickCheck.quickCheck `mapM_` [prop_balance, prop_varianceRelocated, prop_varianceScaled, prop_varianceOrder, prop_equivalence] + where + prop_nC0, prop_nC1, prop_sum :: Math.Implementations.Factorial.Algorithm -> Integer -> Test.QuickCheck.Property + prop_nC0 algorithm n = Test.QuickCheck.label "prop_nC0" $ Math.Statistics.nCr algorithm (abs n) 0 == 1 - prop_nC1 algorithm i = Test.QuickCheck.label "prop_nC1" $ Math.Statistics.nCr algorithm n 1 == n where - n = 1 + abs i + prop_nC1 algorithm i = Test.QuickCheck.label "prop_nC1" $ Math.Statistics.nCr algorithm n 1 == n where + n = 1 + abs i - prop_sum algorithm i = Test.QuickCheck.label "prop_sum" $ sum (Math.Statistics.nCr algorithm n `map` [0 .. n]) == 2 ^ n where - n = 1 + abs i + prop_sum algorithm i = Test.QuickCheck.label "prop_sum" $ sum (Math.Statistics.nCr algorithm n `map` [0 .. n]) == 2 ^ n where + n = 1 + abs i - prop_symmetry, prop_prime :: Math.Implementations.Factorial.Algorithm -> (Integer, Integer) -> Test.QuickCheck.Property - prop_symmetry algorithm (i, j) = Test.QuickCheck.label "prop_symmetry" $ Math.Statistics.nCr algorithm n r == Math.Statistics.nCr algorithm n (n - r) where - [r, n] = Data.List.sort $ map abs [i, j] + prop_symmetry, prop_prime :: Math.Implementations.Factorial.Algorithm -> (Integer, Integer) -> Test.QuickCheck.Property + prop_symmetry algorithm (i, j) = Test.QuickCheck.label "prop_symmetry" $ Math.Statistics.nCr algorithm n r == Math.Statistics.nCr algorithm n (n - r) where + [r, n] = Data.List.sort $ map abs [i, j] - prop_prime algorithm (i, j) = r `notElem` [0, n] ==> Test.QuickCheck.label "prop_prime" $ (Math.Statistics.nCr algorithm n r `mod` n) == 0 where - n = Data.Numbers.Primes.primes !! fromIntegral (i `mod` 500000) - r = j `mod` n --Ensure r is smaller than n. + prop_prime algorithm (i, j) = r `notElem` [0, n] ==> Test.QuickCheck.label "prop_prime" $ (Math.Statistics.nCr algorithm n r `mod` n) == 0 where + n = Data.Numbers.Primes.primes !! fromIntegral (i `mod` 500000) + r = j `mod` n --Ensure r is smaller than n. - prop_nP0, prop_nP1 :: Integer -> Test.QuickCheck.Property - prop_nP0 n = Test.QuickCheck.label "prop_nP0" $ Math.Statistics.nPr (abs n) 0 == 1 + prop_nP0, prop_nP1 :: Integer -> Test.QuickCheck.Property + prop_nP0 n = Test.QuickCheck.label "prop_nP0" $ Math.Statistics.nPr (abs n) 0 == 1 - prop_nP1 i = Test.QuickCheck.label "prop_nP1" $ Math.Statistics.nPr n 1 == n where - n = 1 + abs i + prop_nP1 i = Test.QuickCheck.label "prop_nP1" $ Math.Statistics.nPr n 1 == n where + n = 1 + abs i - prop_balance :: [Integer] -> Test.QuickCheck.Property - prop_balance l = not (null l) ==> Test.QuickCheck.label "prop_balance" . (< 1e-11 {-rounding errors-}) . abs . sum $ map (\i -> fromIntegral i - (Math.Statistics.mean l :: Double)) l + 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_balance, prop_varianceRelocated, prop_varianceScaled, prop_varianceOrder, prop_equivalence :: [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) diff --git a/src/Main.hs b/src/Main.hs --- a/src/Main.hs +++ b/src/Main.hs @@ -98,7 +98,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, 2] [] + Distribution.Package.pkgVersion = Distribution.Version.Version [0, 0, 0, 3] [] } printUsage = System.IO.hPutStrLn System.IO.stderr usage >> System.exitWith System.ExitSuccess