factory 0.2.1.0 → 0.2.1.1

raw patch · 59 files changed

+371/−431 lines, 59 filesdep ~basedep ~parallelsetup-changedPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependency ranges changed: base, parallel

API changes (from Hackage documentation)

Files

• Setup.hs +2/−2
• changelog +9/−0
• factory.cabal +8/−14
• makefile +3/−3
• src/Factory/Data/Exponential.hs +3/−3
• src/Factory/Data/Interval.hs +10/−25
• src/Factory/Data/MonicPolynomial.hs +4/−4
• src/Factory/Data/Monomial.hs +8/−7
• src/Factory/Data/Polynomial.hs +10/−9
• src/Factory/Data/PrimeFactors.hs +5/−5
• src/Factory/Data/PrimeWheel.hs +7/−7
• src/Factory/Data/QuotientRing.hs +1/−1
• src/Factory/Data/Ring.hs +8/−8
• src/Factory/Math/ArithmeticGeometricMean.hs +3/−10
• src/Factory/Math/DivideAndConquer.hs +6/−14
• src/Factory/Math/Hyperoperation.hs +3/−3
• src/Factory/Math/Implementations/Factorial.hs +2/−2
• src/Factory/Math/Implementations/Pi/BBP/Implementation.hs +2/−2
• src/Factory/Math/Implementations/Pi/Borwein/Borwein1993.hs +6/−6
• src/Factory/Math/Implementations/Pi/Borwein/Implementation.hs +4/−10
• src/Factory/Math/Implementations/Pi/Borwein/Series.hs +2/−2
• src/Factory/Math/Implementations/Pi/Ramanujan/Chudnovsky.hs +4/−4
• src/Factory/Math/Implementations/Pi/Ramanujan/Classic.hs +4/−4
• src/Factory/Math/Implementations/Pi/Ramanujan/Implementation.hs +5/−11
• src/Factory/Math/Implementations/Pi/Spigot/Gosper.hs +1/−1
• src/Factory/Math/Implementations/Pi/Spigot/Spigot.hs +5/−5
• src/Factory/Math/Implementations/Primality.hs +20/−29
• src/Factory/Math/Implementations/PrimeFactorisation.hs +10/−18
• src/Factory/Math/Implementations/Primes/Algorithm.hs +2/−2
• src/Factory/Math/Implementations/Primes/SieveOfAtkin.hs +17/−23
• src/Factory/Math/Implementations/Primes/SieveOfEratosthenes.hs +5/−5
• src/Factory/Math/Implementations/Primes/TrialDivision.hs +2/−4
• src/Factory/Math/Implementations/Primes/TurnersSieve.hs +2/−2
• src/Factory/Math/Implementations/SquareRoot.hs +12/−12
• src/Factory/Math/MultiplicativeOrder.hs +1/−1
• src/Factory/Math/PerfectPower.hs +10/−10
• src/Factory/Math/Pi.hs +1/−1
• src/Factory/Math/Power.hs +2/−2
• src/Factory/Math/Precision.hs +1/−1
• src/Factory/Math/Primality.hs +1/−1
• src/Factory/Math/Primes.hs +1/−1
• src/Factory/Math/Probability.hs +8/−8
• src/Factory/Math/Radix.hs +6/−6
• src/Factory/Math/SquareRoot.hs +4/−4
• src/Factory/Math/Statistics.hs +57/−25
• src/Factory/Math/Summation.hs +5/−19
• src/Factory/Test/Performance/Factorial.hs +1/−1
• src/Factory/Test/QuickCheck/ArithmeticGeometricMean.hs +2/−10
• src/Factory/Test/QuickCheck/Factorial.hs +0/−4
• src/Factory/Test/QuickCheck/MonicPolynomial.hs +0/−4
• src/Factory/Test/QuickCheck/Pi.hs +2/−20
• src/Factory/Test/QuickCheck/Polynomial.hs +2/−7
• src/Factory/Test/QuickCheck/Primality.hs +3/−8
• src/Factory/Test/QuickCheck/PrimeFactorisation.hs +1/−5
• src/Factory/Test/QuickCheck/Primes.hs +3/−8
• src/Factory/Test/QuickCheck/Probability.hs +24/−9
• src/Factory/Test/QuickCheck/SquareRoot.hs +4/−9
• src/Factory/Test/QuickCheck/Statistics.hs +30/−3
• src/Main.hs +7/−7

@@ -1,5 +1,5 @@ #!/usr/bin/env runhaskell -import qualified	Distribution.Simple+import qualified Distribution.Simple -main	= Distribution.Simple.defaultMain+main = Distribution.Simple.defaultMain

@@ -76,4 +76,13 @@ 	* 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".+0.2.1.1+	* Added 'Factory.Test.QuickCheck.Probability.prop_logNormalDistributionEqual'.+	* Removed /INLINE/ pragma from 'Factory.Math.Implementations.Primes.TrialDivision.isIndivisibleBy', since to be effective it must be called with fully applied parameters (which it isn't).+	* Un eta-reduced 'Factory.Math.Power.square', since we want it to be inlined when called with one argument.+	* Tested with 'haskell-platform-2013.2.0.0'.+	* Replaced preprocessor-directives with 'build-depends' constraints in 'factory.cabal'.+	* Added function 'Factory.Math.Statistics.getWeightedMean' & corresponding tests in module "Factory.Test.QuickCheck.Statistics".+	* Since '(<$>)' is exported from the Prelude from 'base-4.8', imported "Prelude" hiding '(<*>)' into module "Factory.Data.Monomial", since this symbol is defined locally for other purposes.+	* Either replaced instances of '(<$>)' with 'fmap' to avoid ambiguity between "Control.Applicative" & "Prelude" which (from 'base-4.8') also exports this symbol, or hid the symbol when importing the "Prelude".. 

@@ -1,6 +1,6 @@---Package-properties+-- Package-properties Name:			factory-Version:		0.2.1.0+Version:		0.2.1.1 Cabal-Version:		>= 1.6 Copyright:		(C) 2011-2013 Dr. Alistair Ward License:		GPL@@ -11,12 +11,13 @@ Build-Type:		Simple Description:		A library of number-theory functions, for; factorials, square-roots, Pi and primes. Category:		Math, Number Theory-Tested-With:		GHC == 7.4+Tested-With:		GHC == 7.4, GHC == 7.6, GHC == 7.10 Homepage:		http://functionalley.eu Maintainer:		factory <at> functionalley <dot> eu Bug-reports:		factory <at> functionalley <dot> eu Extra-Source-Files:	changelog, copyright, makefile +-- Turn on using: 'runhaskell ./Setup.hs configure -f llvm'. flag llvm     Description:	Whether the 'llvm' compiler-backend has been installed and is required for code-generation.     manual:		True@@ -90,19 +91,15 @@      Build-depends:         array,-        base == 4.*,+        base >= 4.3 && < 5,         deepseq >= 1.1,         containers,+        parallel >= 3.0,         primes >= 0.1,         random,         toolshed >= 0.13 -    if flag(threaded)-        Build-depends:	parallel >= 3.0-    else-        Build-depends:	parallel--    GHC-options:	-Wall -O2+    GHC-options:	-Wall -O2 -fno-warn-tabs      if impl(ghc >= 7.4.1)         GHC-prof-options:	-prof -fprof-auto -fprof-cafs@@ -150,14 +147,11 @@         Cabal >= 1.6 && < 2,         QuickCheck >= 2.2 -    GHC-options:	-Wall -O2+    GHC-options:	-Wall -O2 -fno-warn-tabs     GHC-prof-options:	-prof -auto-all -caf-all      if flag(threaded)         GHC-options:	-threaded--        if impl(ghc >= 6.12)-            GHC-options:	-feager-blackholing      if impl(ghc >= 7.0)         GHC-options:	-rtsopts

@@ -1,4 +1,4 @@-# Copyright (C) 2011 Dr. Alistair Ward+# Copyright (C) 2011-2104 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@@ -40,10 +40,10 @@ 	PATH=~/.cabal/bin:$$PATH runhaskell Setup $@ --hyperlink-source	#Amend path to find 'HsColour', as required for 'hyperlink-source'.  hlint:-	@$@ -i 'Use &&' -i 'Reduce duplication' -i 'Redundant bracket' src/+	@$@ -i 'Use &&' -i 'Reduce duplication' -i 'Redundant bracket' src/ +RTS -N  sdist:-	runhaskell Setup $@+	TAR_OPTIONS='--format=ustar' runhaskell Setup $@  check: sdist 	cabal upload --check --verbose=3 dist/*.tar.gz;

@@ -43,8 +43,8 @@  import qualified	Control.Arrow -infix 4 =~	--Same as (==).-infixr 8 <^	--Same as (^).+infix 4 =~	-- Same as (==).+infixr 8 <^	-- Same as (^).  -- | Describes an /exponential/, in terms of its /base/ and /exponent/. type Exponential base exponent	= (base, exponent)@@ -70,7 +70,7 @@ -- | Evaluate the specified 'Exponential', returning the resulting number. {-# INLINE evaluate #-} evaluate :: (Num base, Integral exponent) => Exponential base exponent -> base-evaluate	= uncurry (^)+evaluate	= uncurry (^)	-- CAVEAT: in this eta-reduced form, it'll only be inlined when called without arguments.  -- | True if the /bases/ are equal. (=~) :: Eq base => Exponential base exponent -> Exponential base exponent -> Bool

@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {- 	Copyright (C) 2011 Dr. Alistair Ward @@ -62,22 +61,12 @@ ) where  import			Control.Arrow((***), (&&&))+import qualified	Control.Parallel.Strategies import qualified	Data.Monoid import qualified	Data.Ratio+import qualified	Data.Tuple import qualified	ToolShed.Data.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) @@ -121,7 +110,7 @@ -- | 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+	| isReversed b	= Data.Tuple.swap b 	| otherwise	= b  -- | Bisect the /interval/ at the specified /end-point/; which should be between the two existing /end-points/.@@ -150,11 +139,11 @@ {- | 	* Converts 'Interval' to a list by enumerating the values. -	* CAVEAT: produces rather odd results for 'Fractional' types, but no stranger than considering such types 'Enum'erable.+	* CAVEAT: produces rather odd results for 'Fractional' types, but no stranger than considering such types Enumerable in the first place. -} {-# INLINE toList #-} toList :: Enum endPoint => Interval endPoint -> [endPoint]-toList	= uncurry enumFromTo+toList	= uncurry enumFromTo	-- CAVEAT: in this eta-reduced form, it'll only be inlined when called without arguments.  {- | 	* Reduces 'Interval' to a single integral value encapsulated in a 'Data.Monoid.Monoid',@@ -188,16 +177,12 @@ 	| 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 (+			| 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 . Control.Parallel.Strategies.withStrategy ( 				Control.Parallel.Strategies.parTuple2 Control.Parallel.Strategies.rseq Control.Parallel.Strategies.rseq-			) .-#endif-			ToolShed.Data.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.+			) . ToolShed.Data.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'.

@@ -26,7 +26,7 @@ module Factory.Data.MonicPolynomial( -- * Types -- ** Data-types,-	MonicPolynomial(getPolynomial),	--Hide the data-constructor.+	MonicPolynomial(getPolynomial),	-- Hide the data-constructor. -- * Functions -- ** Constructors 	mkMonicPolynomial@@ -72,11 +72,11 @@ 	Show	e  ) => Data.Ring.Ring (MonicPolynomial c e)	where 	MkMonicPolynomial l =*= MkMonicPolynomial r	= MkMonicPolynomial $ l =*= r-	MkMonicPolynomial l =+= MkMonicPolynomial r	= mkMonicPolynomial $ l =+= r	--CAVEAT: potentially non-monic.---	additiveInverse (MkMonicPolynomial p)		= MkMonicPolynomial $ Data.Ring.additiveInverse p	--CAVEAT: not monic !+	MkMonicPolynomial l =+= MkMonicPolynomial r	= mkMonicPolynomial $ l =+= r	-- CAVEAT: potentially non-monic.+--	additiveInverse (MkMonicPolynomial p)		= MkMonicPolynomial $ Data.Ring.additiveInverse p	-- CAVEAT: not monic ! 	additiveInverse _				= error "Factory.Data.MonicPolynomial.additiveInverse:\tresult isn't monic" 	multiplicativeIdentity				= MkMonicPolynomial Data.Ring.multiplicativeIdentity-	additiveIdentity				= MkMonicPolynomial Data.Ring.additiveIdentity	--CAVEAT: not monic !+	additiveIdentity				= MkMonicPolynomial Data.Ring.additiveIdentity	-- CAVEAT: not monic !  -- Since the /leading term/ of the /denominator/ is one, the /coefficient/ isn't required to implement 'Fractional'. instance (

@@ -1,5 +1,5 @@ {--	Copyright (C) 2011 Dr. Alistair Ward+	Copyright (C) 2011-2015 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@@ -48,12 +48,13 @@ 	isMonomial ) where +import Prelude hiding ((<*>))	-- The "Prelude" from 'base-4.8' exports this symbol. import qualified	Control.Arrow -infix 4 <=>	--Same as (==).-infix 4 =~	--Same as (==).-infixl 7 </>	--Same as (/).-infixl 7 <*>	--Same as (*).+infix 4 <=>	-- Same as (==).+infix 4 =~	-- Same as (==).+infixl 7 </>	-- Same as (/).+infixl 7 <*>	-- Same as (*).  {- | 	* The type of an arbitrary monomial.@@ -117,7 +118,7 @@ 	=> Monomial c e 	-> c	-- ^ The magnitude of the shift. 	-> Monomial c e---m `shiftCoefficient` i	= Control.Arrow.first (+ i) m	--CAVEAT: Too slow.+-- m `shiftCoefficient` i	= Control.Arrow.first (+ i) m	-- CAVEAT: Too slow. (c, e) `shiftCoefficient` i	= (c + i, e)  -- | Shift the /exponent/, by the specified amount.@@ -126,7 +127,7 @@ 	=> Monomial c e 	-> e	-- ^ The magnitude of the shift. 	-> Monomial c e---m `shiftExponent` i	= Control.Arrow.second (+ i) m	--CAVEAT: Too slow.+-- m `shiftExponent` i	= Control.Arrow.second (+ i) m	-- CAVEAT: Too slow. (c, e) `shiftExponent` i	= (c, e + i)  -- | Negate the coefficient.

@@ -1,5 +1,5 @@ {--	Copyright (C) 2011 Dr. Alistair Ward+	Copyright (C) 2011-2015 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@@ -65,6 +65,7 @@ 	isZero ) where +import Prelude hiding ((<*>))	-- The "Prelude" from 'base-4.8' exports this symbol. import			Control.Arrow((&&&)) import qualified	Control.Arrow import qualified	Data.List@@ -74,7 +75,7 @@ import			Factory.Data.Ring((=*=), (=+=), (=-=)) import qualified	Factory.Data.Ring		as Data.Ring -infixl 7 *=	--Same as (*).+infixl 7 *=	-- Same as (*).  -- | The guts of a 'Polynomial'. type MonomialList coefficient exponent	= [Data.Monomial.Monomial coefficient exponent]@@ -108,8 +109,8 @@ 	MkPolynomial [] =*= _	= zero 	_ =*= MkPolynomial []	= zero 	polynomialL =*= polynomialR---		| polynomialL == one			= polynomialR	--Counterproductive.---		| polynomialR == one			= polynomialL	--Counterproductive.+--		| polynomialL == one			= polynomialR	-- Counterproductive.+--		| polynomialR == one			= polynomialL	-- Counterproductive. 		| terms polynomialL > terms polynomialR	= polynomialL `times` polynomialR 		| otherwise				= polynomialR `times` polynomialL 		where@@ -167,9 +168,9 @@ 	_ `quotRem'` MkPolynomial []		= error "Factory.Data.Polynomial.quotRem':\tzero denominator." 	polynomialN `quotRem'` polynomialD	= longDivide polynomialN	where --		longDivide :: (Fractional c, Num e, Ord e) => Polynomial c e -> (Polynomial c e, Polynomial c e)-		longDivide (MkPolynomial [])	= (zero, zero)	--Exactly divides.+		longDivide (MkPolynomial [])	= (zero, zero)	-- Exactly divides. 		longDivide numerator-			| Data.Monomial.getExponent quotient < 0	= (zero, numerator)	--Indivisible remainder.+			| Data.Monomial.getExponent quotient < 0	= (zero, numerator)	-- Indivisible remainder. 			| otherwise					= Control.Arrow.first (lift (quotient :)) $ longDivide (numerator =-= polynomialD *= quotient ) 			where --				quotient :: (Fractional c, Num e) => Data.Monomial.Monomial c e@@ -260,7 +261,7 @@ 	* <http://en.wikipedia.org/wiki/Monic_polynomial#Classifications>. -} isMonic :: (Eq c, Num c) => Polynomial c e -> Bool-isMonic (MkPolynomial [])	= False	--All coefficients are zero, and have therefore been removed.+isMonic (MkPolynomial [])	= False	-- All coefficients are zero, and have therefore been removed. isMonic p			= (== 1) . Data.Monomial.getCoefficient $ getLeadingTerm p  -- | True if there are zero terms.@@ -301,7 +302,7 @@ 	* <http://mathworld.wolfram.com/PolynomialDegree.html>. -} getDegree :: Num e => Polynomial c e -> e-getDegree (MkPolynomial [])	= -1	--CAVEAT: debatable, but makes some operations more robust and consistent.+getDegree (MkPolynomial [])	= -1	-- CAVEAT: debatable, but makes some operations more robust and consistent. getDegree p			= Data.Monomial.getExponent $ getLeadingTerm p  {- |@@ -329,7 +330,7 @@ raiseModulo _ 0 modulus			= mkConstant $ 1 `mod` modulus raiseModulo polynomial power modulus 	| power < 0			= error $ "Factory.Data.Polynomial.raiseModulo:\tthe result isn't guaranteed to be a polynomial, for power=" ++ show power-	| first `elem` [zero, one]	= first	--Eg 'raiseModulo (mkPolynomial [(3,1)]) 100 3' or 'raiseModulo (mkPolynomial [(3,1),(1,0)]) 100 3'.+	| first `elem` [zero, one]	= first	-- Eg 'raiseModulo (mkPolynomial [(3,1)]) 100 3' or 'raiseModulo (mkPolynomial [(3,1),(1,0)]) 100 3'. 	| otherwise			= slave power 	where --		first :: Integral c => Polynomial c e

@@ -50,8 +50,8 @@ import			Factory.Data.Exponential((<^), (=~)) import qualified	ToolShed.Data.List -infixl 7 >/<, >*<	--Same as (/).-infixr 8 >^		--Same as (^).+infixl 7 >/<, >*<	-- Same as (/).+infixr 8 >^		-- Same as (^).  {- | 	* Each element of this list represents one /prime-factor/, expressed as an /exponential/ with a /prime/ base, of the original integer.@@ -70,7 +70,7 @@  -- | Multiplies 'Data.Exponential.Exponential's of similar /base/. reduceSorted :: (Eq base, Num exponent) => Factors base exponent -> Factors base exponent---reduceSorted	= map (Data.Exponential.getBase . head &&& sumExponents) . Data.List.groupBy (=~)	--Slow+-- reduceSorted	= map (Data.Exponential.getBase . head &&& sumExponents) . Data.List.groupBy (=~)	-- Slow reduceSorted []	= [] reduceSorted (x : xs) 	| null matched	= x : reduceSorted remainder@@ -91,8 +91,8 @@ insert' e []		= [e] insert' e l@(x : xs)	= case Data.Ord.comparing Data.Exponential.getBase e x of 	LT	-> e : l-	GT	-> x : insert' e xs	--Recurse.-	_	-> Control.Arrow.second (+ Data.Exponential.getExponent e) x : xs	--Multiply by adding exponents.+	GT	-> x : insert' e xs	-- Recurse.+	_	-> Control.Arrow.second (+ Data.Exponential.getExponent e) x : xs	-- Multiply by adding exponents.  {- | 	* Multiplies two lists each representing a product of /prime factors/, and sorted by increasing /base/.

@@ -103,10 +103,10 @@ 	where 		sieve :: Int -> NPrimes -> Repository -> [Int] 		sieve candidate found repository	= case Data.IntMap.lookup candidate repository of-			Just primeMultiples	-> sieve' found . insertUniq primeMultiples $ Data.IntMap.delete candidate repository	--Re-insert subsequent multiples.+			Just primeMultiples	-> sieve' found . insertUniq primeMultiples $ Data.IntMap.delete candidate repository	-- Re-insert subsequent multiples. 			Nothing {-prime-}	-> let 				found'		= succ found-				(key : values)	= iterate (+ gap * candidate) $ candidate ^ (2 :: Int)	--Generate a sequence of prime-multiples, starting from its square.+				(key : values)	= iterate (+ gap * candidate) $ candidate ^ (2 :: Int)	-- Generate a sequence of prime-multiples, starting from its square. 			 in candidate : sieve' found' ( 				if found' >= required 					then repository@@ -114,10 +114,10 @@ 			 ) 			where 				gap :: Int-				gap	= 2	--For efficiency, only sieve odd integers.+				gap	= 2	-- For efficiency, only sieve odd integers.  				sieve' :: NPrimes -> Repository -> [Int]-				sieve'	= sieve $ candidate + gap	--Tail-recurse.+				sieve'	= sieve $ candidate + gap	-- Tail-recurse.  				insertUniq :: PrimeMultiples Int -> Repository -> Repository 				insertUniq l m	= insert $ dropWhile (`Data.IntMap.member` m) l	where@@ -147,8 +147,8 @@ >	nPrimes	Gaps >	======	==== >	0	[1]->	1	[2]	--The terminal gap for all subsequent wheels is '2'; [(succ circumference `mod` circumference) - (pred circumference `mod` circumference)].->	2	[4,2]	--Both points are on the axis, so the symmetry isn't yet clear.+>	1	[2]	-- The terminal gap for all subsequent wheels is '2'; [(succ circumference `mod` circumference) - (pred circumference `mod` circumference)].+>	2	[4,2]	-- Both points are on the axis, so the symmetry isn't yet clear. >	3	[6,4,2,4,2,4,6,2] >	4	[10,2,4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,8,6,4,6,2,4,6,2,6,6,4,2,4,6,2,6,4,2,4,2,10,2] @@ -162,7 +162,7 @@ 	| otherwise	= primeWheel 	where 		(primeComponents, coprimeCandidates)	= (map fromIntegral *** map fromIntegral . Data.List.genericTake (getSpokeCount primeWheel)) $ findCoprimes nPrimes-		primeWheel				= MkPrimeWheel primeComponents $ zipWith (-) coprimeCandidates $ 1 : coprimeCandidates	--Measure the gaps between candidate primes.+		primeWheel				= MkPrimeWheel primeComponents $ zipWith (-) coprimeCandidates $ 1 : coprimeCandidates	-- Measure the gaps between candidate primes.  -- | Couples a candidate prime with a /rolling wheel/, to define the distance rolled. type Distance i	= (i, [i])

@@ -67,7 +67,7 @@ 	-> q	-- ^ Modulus. 	-> Bool areCongruentModulo l r modulus-	| l == r	= True	--Only required for efficiency.+	| l == r	= True	-- Only required for efficiency. 	| otherwise	= (l =-= r) `isDivisibleBy` modulus  -- | True if the second operand /divides/ the first.

@@ -43,10 +43,10 @@ import qualified	Data.Monoid import qualified	Factory.Math.DivideAndConquer	as Math.DivideAndConquer -infixl 6 =+=	--Same as (+).-infixl 6 =-=	--Same as (-).-infixl 7 =*=	--Same as (*).-infixr 8 =^	--Same as (^).+infixl 6 =+=	-- Same as (+).+infixl 6 =-=	-- Same as (-).+infixl 7 =*=	-- Same as (*).+infixr 8 =^	-- Same as (^).  {- | 	* Define both the operations applicable to all members of the /ring/, and its mandatory members.@@ -61,10 +61,10 @@ 	additiveIdentity	:: r		-- ^ The /identity/-member under addition (AKA /zero/); <http://en.wikipedia.org/wiki/Additive_identity>.  	(=-=) :: r -> r -> r			-- ^ Subtract the two specified /ring/-members.-	l =-= r	= l =+= additiveInverse r	--Default implementation.+	l =-= r	= l =+= additiveInverse r	-- Default implementation.  	square :: r -> r			-- ^ Square the ring.-	square r	= r =*= r		--Default implementation; there may be a more efficient one.+	square r	= r =*= r		-- Default implementation; there may be a more efficient one.  {- | 	* Raise a /ring/-member to the specified positive integral power.@@ -99,7 +99,7 @@  -- | Returns the /product/ of the list of /ring/-members. product' :: Ring r => Math.DivideAndConquer.BisectionRatio -> Math.DivideAndConquer.MinLength -> [r] -> r---product' _ _			= getProduct . Data.Monoid.mconcat . map MkProduct+-- product' _ _			= getProduct . Data.Monoid.mconcat . map MkProduct product' ratio minLength	= getProduct . Math.DivideAndConquer.divideAndConquer ratio minLength . map MkProduct  -- | Does for 'Ring', what 'Data.Monoid.Sum' does for type 'Num', in that it makes it an instance of 'Data.Monoid.Monoid' under addition.@@ -113,6 +113,6 @@  -- | Returns the /sum/ of the list of /ring/-members. sum' :: Ring r => Math.DivideAndConquer.BisectionRatio -> Math.DivideAndConquer.MinLength -> [r] -> r---sum' _ _		= getSum . Data.Monoid.mconcat . map MkSum+-- sum' _ _		= getSum . Data.Monoid.mconcat . map MkSum sum' ratio minLength	= getSum . Math.DivideAndConquer.divideAndConquer ratio minLength . map MkSum 

@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {- 	Copyright (C) 2011 Dr. Alistair Ward @@ -38,13 +37,10 @@ ) where  import			Control.Arrow((&&&))+import qualified	Control.Parallel.Strategies import qualified	Factory.Math.Precision	as Math.Precision import qualified	Factory.Math.SquareRoot	as Math.SquareRoot -#if MIN_VERSION_parallel(3,0,0)-import qualified	Control.Parallel.Strategies-#endif- -- | The type of the /arithmetic mean/; <http://en.wikipedia.org/wiki/Arithmetic_mean>. type ArithmeticMean	= Rational @@ -72,7 +68,7 @@ 	| spread agm == 0	= repeat agm 	| otherwise		= let 		simplify :: Rational -> Rational-		simplify	= Math.Precision.simplify (pred decimalDigits {-ignore single integral digit-})	--This makes a gigantic difference to performance.+		simplify	= Math.Precision.simplify (pred decimalDigits {-ignore single integral digit-})	-- This makes a gigantic difference to performance.  		findArithmeticMean :: AGM -> ArithmeticMean 		findArithmeticMean	= (/ 2) . uncurry (+)@@ -80,12 +76,9 @@ 		findGeometricMean :: AGM -> GeometricMean 		findGeometricMean	= Math.SquareRoot.squareRoot squareRootAlgorithm decimalDigits . uncurry (*) 	in iterate (-#if MIN_VERSION_parallel(3,0,0) 		Control.Parallel.Strategies.withStrategy ( 			Control.Parallel.Strategies.parTuple2 Control.Parallel.Strategies.rdeepseq Control.Parallel.Strategies.rdeepseq-		) .-#endif-		(simplify . findArithmeticMean &&& simplify . findGeometricMean)+		) . (simplify . findArithmeticMean &&& simplify . findGeometricMean) 	) agm  -- | Returns the bounds within which the 'AGM' has been constrained.

@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {- 	Copyright (C) 2011 Dr. Alistair Ward @@ -40,13 +39,10 @@ ) where  import			Control.Arrow((***))+import qualified	Control.Parallel.Strategies import qualified	Data.Monoid import qualified	Data.Ratio -#if MIN_VERSION_parallel(3,0,0)-import qualified	Control.Parallel.Strategies-#endif- {- | 	* The ratio of the original list-length at which to bisect. @@ -81,8 +77,8 @@ 	-> monoid		-- ^ The resulting scalar. divideAndConquer bisectionRatio minLength l 	| any ($ apportion minLength) [-		(< 1),			--The left-hand list may be null.-		(> pred minLength)	--The right-hand list may be null.+		(< 1),			-- The left-hand list may be null.+		(> pred minLength)	-- The right-hand list may be null. 	]		= error $ "Factory.Math.DivideAndConquer.divideAndConquer:\tbisectionRatio='" ++ show bisectionRatio ++ "' is incompatible with minLength=" ++ show minLength ++ "." 	| otherwise	= slave (length l) l 	where@@ -90,14 +86,10 @@ 		apportion list	= (list * Data.Ratio.numerator bisectionRatio) `div` Data.Ratio.denominator bisectionRatio  		slave len list-			| len <= minLength	= Data.Monoid.mconcat list	--Fold the monoid's binary operator over the list.-			| otherwise		= uncurry Data.Monoid.mappend .-#if MIN_VERSION_parallel(3,0,0)-			Control.Parallel.Strategies.withStrategy (+			| len <= minLength	= Data.Monoid.mconcat list	-- Fold the monoid's binary operator over the list.+			| otherwise		= uncurry Data.Monoid.mappend . Control.Parallel.Strategies.withStrategy ( 				Control.Parallel.Strategies.parTuple2 Control.Parallel.Strategies.rseq Control.Parallel.Strategies.rseq-			) .-#endif-			(slave cut *** slave (len - cut)) $ splitAt cut list	where	--Apply the monoid's binary operator to the two operands resulting from bisection.+			) . (slave cut *** slave (len - cut)) $ splitAt cut list	where	-- Apply the monoid's binary operator to the two operands resulting from bisection. 				cut	= apportion len  {- |

@@ -57,7 +57,7 @@ -} type HyperExponent	= Base -succession, addition, multiplication, exponentiation, tetration, pentation, hexation :: Int	--Arbitrarily.+succession, addition, multiplication, exponentiation, tetration, pentation, hexation :: Int	-- Arbitrarily. (succession : addition : multiplication : exponentiation : tetration : pentation : hexation : _)	= [0 ..]  {- |@@ -71,7 +71,7 @@ powerTower 0 hyperExponent 	| even hyperExponent	= 1 	| otherwise		= 0-powerTower _ (-1)	= 0	--The only negative hyper-exponent for which there's a consistent result.+powerTower _ (-1)	= 0	-- The only negative hyper-exponent for which there's a consistent result. powerTower base hyperExponent 	| base < 0 && hyperExponent > 1	= error $ "Factory.Math.Hyperoperation.powerTower:\tundefined for negative base; " ++ show base 	| otherwise			= Data.List.genericIndex (iterate (base ^) 1) hyperExponent@@ -95,7 +95,7 @@ 			3 {-exponentiation-}	-> base ^ e 			4 {-tetration-}		-> base `powerTower` e 			_-				| e' == e	-> tetration ^# e'	--To which it would otherwise be reduced by laborious recursion.+				| e' == e	-> tetration ^# e'	-- To which it would otherwise be reduced by laborious recursion. 				| otherwise	-> pred r ^# e' 				where 					e'	= {-fromIntegral $-} r ^# pred e

@@ -50,7 +50,7 @@ import qualified	Factory.Math.Factorial		as Math.Factorial import qualified	ToolShed.Defaultable -infixl 7 !/!	--Same as (/).+infixl 7 !/!	-- Same as (/).  -- | The algorithms by which /factorial/ has been implemented. data Algorithm	=@@ -133,6 +133,6 @@ 	| numerator <= 1		= recip . fromIntegral $ Math.Factorial.factorial (ToolShed.Defaultable.defaultValue :: Algorithm) denominator 	| denominator <= 1		= fromIntegral $ Math.Factorial.factorial (ToolShed.Defaultable.defaultValue :: Algorithm) numerator 	| numerator == denominator	= 1-	| numerator < denominator	= recip $ denominator !/! numerator	--Recurse.+	| numerator < denominator	= recip $ denominator !/! numerator	-- Recurse. 	| otherwise			= fromIntegral $ Data.Interval.product' (recip 2) 64 (succ denominator, numerator) 

@@ -47,10 +47,10 @@ 	Math.Implementations.Pi.BBP.Series.base			= base } decimalDigits		= (seriesScalingFactor *) . Math.Summation.sum' 8 . take ( 	Math.Precision.getTermsRequired (-		recip . fromIntegral $ abs {-potentially negative-} base	--The convergence-rate.+		recip . fromIntegral $ abs {-potentially negative-} base	-- The convergence-rate. 	) decimalDigits  ) . zipWith (*) (-	iterate (/ fromIntegral base) 1	--Generate the scaling-ratio, between successive terms.+	iterate (/ fromIntegral base) 1	-- Generate the scaling-ratio, between successive terms.  ) $ map ( 	sum . zipWith (%) numerators . getDenominators  ) [0 ..]

@@ -25,10 +25,10 @@ 	series ) where ---import		Control.Arrow((***))+-- import		Control.Arrow((***)) import			Data.Ratio((%))---import		Factory.Data.PrimeFactors((>*<), (>/<), (>^))---import qualified	Factory.Data.PrimeFactors			as Data.PrimeFactors+-- import		Factory.Data.PrimeFactors((>*<), (>/<), (>^))+-- import qualified	Factory.Data.PrimeFactors			as Data.PrimeFactors import qualified	Factory.Math.Factorial				as Math.Factorial import qualified	Factory.Math.Implementations.Factorial		as Math.Implementations.Factorial import qualified	Factory.Math.Implementations.Pi.Borwein.Series	as Math.Implementations.Pi.Borwein.Series@@ -41,7 +41,7 @@ series = Math.Implementations.Pi.Borwein.Series.MkSeries { 	Math.Implementations.Pi.Borwein.Series.terms			= \squareRootAlgorithm factorialAlgorithm decimalDigits -> let 		simplify, squareRoot :: Rational -> Rational-		simplify	= Math.Precision.simplify $ pred decimalDigits {-ignore single integral digit-}	--This makes a gigantic difference to performance.+		simplify	= Math.Precision.simplify $ pred decimalDigits {-ignore single integral digit-}	-- This makes a gigantic difference to performance. 		squareRoot	= simplify . Math.SquareRoot.squareRoot squareRootAlgorithm decimalDigits  		sqrt5, a, b, c3 :: Rational@@ -51,7 +51,7 @@ 		b	= 7849910453496627210289749000 + 3510586678260932028965606400 * sqrt5 + 2515968 * squareRoot (3110 * (6260208323789001636993322654444020882161 + 2799650273060444296577206890718825190235 * sqrt5)) 		c3	= simplify . Math.Power.cube $ negate 214772995063512240 - sqrt5 * (96049403338648032 + 1296 * squareRoot (10985234579463550323713318473 + 4912746253692362754607395912 * sqrt5)) 	in (-		squareRoot $ negate c3,	--The factor into which the series must be divided, to yield Pi.+		squareRoot $ negate c3,	-- The factor into which the series must be divided, to yield Pi. 		zipWith ( {- 			\n power -> let@@ -69,5 +69,5 @@ 			) 		) [0 :: Integer ..] $ iterate (* c3) 1 	),-	Math.Implementations.Pi.Borwein.Series.convergenceRate		= 10 ** negate 50	--Empirical.+	Math.Implementations.Pi.Borwein.Series.convergenceRate		= 10 ** negate 50	-- Empirical. }

@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {- 	Copyright (C) 2011 Dr. Alistair Ward @@ -27,13 +26,10 @@ ) where  import qualified	Control.Arrow+import qualified	Control.Parallel.Strategies import qualified	Factory.Math.Implementations.Pi.Borwein.Series	as Math.Implementations.Pi.Borwein.Series import qualified	Factory.Math.Precision				as Math.Precision -#if MIN_VERSION_parallel(3,0,0)-import qualified	Control.Parallel.Strategies-#endif- -- | Returns /Pi/, accurate to the specified number of decimal digits. openR 	:: Math.Implementations.Pi.Borwein.Series.Series squareRootAlgorithm factorialAlgorithm	-- ^ This /Pi/-algorithm is parameterised by the type of other algorithms to use.@@ -44,11 +40,9 @@ openR Math.Implementations.Pi.Borwein.Series.MkSeries { 	Math.Implementations.Pi.Borwein.Series.terms		= terms, 	Math.Implementations.Pi.Borwein.Series.convergenceRate	= convergenceRate-} squareRootAlgorithm factorialAlgorithm decimalDigits	= uncurry (/)-#if MIN_VERSION_parallel(3,0,0)-	. Control.Parallel.Strategies.withStrategy (Control.Parallel.Strategies.parTuple2 Control.Parallel.Strategies.rdeepseq Control.Parallel.Strategies.rdeepseq)-#endif-	. Control.Arrow.second (+} squareRootAlgorithm factorialAlgorithm decimalDigits	= uncurry (/) . Control.Parallel.Strategies.withStrategy (+		Control.Parallel.Strategies.parTuple2 Control.Parallel.Strategies.rdeepseq Control.Parallel.Strategies.rdeepseq+	) . Control.Arrow.second ( 		sum . take ( 			Math.Precision.getTermsRequired convergenceRate decimalDigits 		)

@@ -35,8 +35,8 @@ 		-> factorialAlgorithm 		-> Math.Precision.DecimalDigits 		-> (-			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.+			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. }

@@ -25,10 +25,10 @@ 	series ) where ---import		Control.Arrow((***))+-- import		Control.Arrow((***)) import			Data.Ratio((%))---import		Factory.Data.PrimeFactors((>/<), (>*<), (>^))---import qualified	Factory.Data.PrimeFactors				as Data.PrimeFactors+-- import		Factory.Data.PrimeFactors((>/<), (>*<), (>^))+-- import qualified	Factory.Data.PrimeFactors				as Data.PrimeFactors import qualified	Factory.Math.Factorial					as Math.Factorial import qualified	Factory.Math.Implementations.Factorial			as Math.Implementations.Factorial import qualified	Factory.Math.Implementations.Pi.Ramanujan.Series	as Math.Implementations.Pi.Ramanujan.Series@@ -58,6 +58,6 @@ 		) -- CAVEAT: the order in which these terms are evaluated radically affects performance. 	) [0 ..] $ iterate (* (Math.Power.cube $ negate 640320 :: Integer)) 1, 	Math.Implementations.Pi.Ramanujan.Series.getSeriesScalingFactor	= \squareRootAlgorithm decimalDigits -> 426880 * Math.SquareRoot.squareRoot squareRootAlgorithm decimalDigits (10005 :: Integer),-	Math.Implementations.Pi.Ramanujan.Series.convergenceRate	= 10 ** negate 14.0	--Empirical.+	Math.Implementations.Pi.Ramanujan.Series.convergenceRate	= 10 ** negate 14.0	-- Empirical. } 

@@ -25,10 +25,10 @@ 	series ) where ---import		Control.Arrow((***))+-- import		Control.Arrow((***)) import			Data.Ratio((%))---import		Factory.Data.PrimeFactors((>/<), (>^))---import qualified	Factory.Data.PrimeFactors				as Data.PrimeFactors+-- import		Factory.Data.PrimeFactors((>/<), (>^))+-- import qualified	Factory.Data.PrimeFactors				as Data.PrimeFactors import qualified	Factory.Math.Factorial					as Math.Factorial import qualified	Factory.Math.Implementations.Factorial			as Math.Implementations.Factorial import qualified	Factory.Math.Implementations.Pi.Ramanujan.Series	as Math.Implementations.Pi.Ramanujan.Series@@ -55,6 +55,6 @@ 		) -- CAVEAT: the order in which these terms are evaluated radically affects performance. 	) [0 ..] $ iterate (* toFourthPower 396) 1, 	Math.Implementations.Pi.Ramanujan.Series.getSeriesScalingFactor	= \squareRootAlgorithm decimalDigits -> 9801 / Math.SquareRoot.squareRoot squareRootAlgorithm decimalDigits (8 :: Integer),-	Math.Implementations.Pi.Ramanujan.Series.convergenceRate	= 10 ** negate 7.9	--Empirical.+	Math.Implementations.Pi.Ramanujan.Series.convergenceRate	= 10 ** negate 7.9	-- Empirical. } 

@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {- 	Copyright (C) 2011 Dr. Alistair Ward @@ -26,14 +25,11 @@ 	openR ) where +import qualified	Control.Parallel.Strategies 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 -#if MIN_VERSION_parallel(3,0,0)-import qualified	Control.Parallel.Strategies-#endif- -- | Returns /Pi/, accurate to the specified number of decimal digits. openR 	:: Math.Implementations.Pi.Ramanujan.Series.Series squareRootAlgorithm factorialAlgorithm	-- ^ This /Pi/-algorithm is parameterised by the type of other algorithms to use.@@ -45,14 +41,12 @@ 	Math.Implementations.Pi.Ramanujan.Series.terms			= terms, 	Math.Implementations.Pi.Ramanujan.Series.getSeriesScalingFactor	= getSeriesScalingFactor, 	Math.Implementations.Pi.Ramanujan.Series.convergenceRate	= convergenceRate-} squareRootAlgorithm factorialAlgorithm decimalDigits	= uncurry (/)-#if MIN_VERSION_parallel(3,0,0)-	$ Control.Parallel.Strategies.withStrategy (Control.Parallel.Strategies.parTuple2 Control.Parallel.Strategies.rdeepseq Control.Parallel.Strategies.rdeepseq)-#endif-	(+} squareRootAlgorithm factorialAlgorithm decimalDigits	= uncurry (/) $ Control.Parallel.Strategies.withStrategy (+		Control.Parallel.Strategies.parTuple2 Control.Parallel.Strategies.rdeepseq Control.Parallel.Strategies.rdeepseq+	) ( 		getSeriesScalingFactor squareRootAlgorithm decimalDigits, 		Math.Summation.sumR 64 . take ( 			Math.Precision.getTermsRequired convergenceRate decimalDigits 		) $ terms factorialAlgorithm-	)+	) -- Pair. 

@@ -33,7 +33,7 @@ series	= Math.Implementations.Pi.Spigot.Series.MkSeries { 	Math.Implementations.Pi.Spigot.Series.baseNumerators	= map (\i -> i * pred (2 * i)) [1 ..], 	Math.Implementations.Pi.Spigot.Series.baseDenominators	= map ((* 3) . (\i -> succ i * (i + 2))) [3, 6 ..],-	Math.Implementations.Pi.Spigot.Series.coefficients	= [3, 8 ..],	--5n - 2+	Math.Implementations.Pi.Spigot.Series.coefficients	= [3, 8 ..],	-- 5n - 2 	Math.Implementations.Pi.Spigot.Series.nTerms		= Math.Precision.getTermsRequired $ 1 / 13 {-empirical convergence-rate-} } 

@@ -87,12 +87,12 @@ -} carryAndDivide :: (Base, I) -> QuotRem -> QuotRem carryAndDivide (base, lhs) rhs-	| n < d		= (0, n)	--In some degenerate cases, the result of the subsequent calculation can be more simply determined.+	| n < d		= (0, n)	-- In some degenerate cases, the result of the subsequent calculation can be more simply determined. 	| otherwise	= Control.Arrow.first (* Data.Ratio.numerator base) $ n `quotRem` d 	where 		d, n :: I 		d	= Data.Ratio.denominator base-		n	= lhs + getQuotient rhs	--Carry numerator from the column to the right and add it to the current digit.+		n	= lhs + getQuotient rhs	-- Carry numerator from the column to the right and add it to the current digit.  {- | 	* Fold 'carryAndDivide', from right to left, over the columns of a row in the spigot-table, continuously checking for overflow.@@ -109,9 +109,9 @@ 	-> [(Base, I)]	-- ^ Data-row. 	-> Pi processColumns series preDigits l-	| overflowMargin > 1	= preDigits ++ nextRow [digit]				--There's neither overflow, nor risk of impact from subsequent overflow.-	| overflowMargin == 1	= nextRow $ preDigits ++ [digit]			--There's no overflow, but risk of impact from subsequent overflow.-	| otherwise		= map ((`rem` decimal) . succ) preDigits ++ nextRow [0]	--Overflow => propagate the excess to previously withheld preDigits.+	| overflowMargin > 1	= preDigits ++ nextRow [digit]				-- There's neither overflow, nor risk of impact from subsequent overflow.+	| overflowMargin == 1	= nextRow $ preDigits ++ [digit]			-- There's no overflow, but risk of impact from subsequent overflow.+	| otherwise		= map ((`rem` decimal) . succ) preDigits ++ nextRow [0]	-- Overflow => propagate the excess to previously withheld preDigits. 	where 		results :: [QuotRem] 		results	= init $ scanr carryAndDivide (0, undefined) l

@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {- 	Copyright (C) 2011 Dr. Alistair Ward @@ -42,6 +41,7 @@  import			Control.Arrow((&&&)) import qualified	Control.DeepSeq+import qualified	Control.Parallel.Strategies import qualified	Data.Numbers.Primes import qualified	Factory.Data.MonicPolynomial		as Data.MonicPolynomial import qualified	Factory.Data.Polynomial			as Data.Polynomial@@ -53,10 +53,6 @@ import qualified	Factory.Math.PrimeFactorisation		as Math.PrimeFactorisation import qualified	ToolShed.Defaultable -#if MIN_VERSION_parallel(3,0,0)-import qualified	Control.Parallel.Strategies-#endif- -- | The algorithms by which /primality/-testing has been implemented. data Algorithm factorisationAlgorithm	= 	AKS factorisationAlgorithm	-- ^ <http://en.wikipedia.org/wiki/AKS_primality_test>.@@ -67,14 +63,14 @@ 	defaultValue	= MillerRabin  instance Math.PrimeFactorisation.Algorithmic factorisationAlgorithm => Math.Primality.Algorithmic (Algorithm factorisationAlgorithm)	where-	isPrime _ 2	= True	--The only even prime.+	isPrime _ 2	= True	-- The only even prime. 	isPrime algorithm candidate 		| candidate < 2 || ( 			any (-				(== 0) . (candidate `rem`)			--The candidate has a small prime-factor, and is therefore composite.+				(== 0) . (candidate `rem`)			-- The candidate has a small prime-factor, and is therefore composite. 			) . filter (-				(candidate >=) . (* 2)				--The candidate must be at least double the small prime, for it to be a potential factor.-			) . take 5 {-arbitrarily-} $ Data.Numbers.Primes.primes	--Excludes even numbers, provided at least the 1st prime is tested.+				(candidate >=) . (* 2)				-- The candidate must be at least double the small prime, for it to be a potential factor.+			) . take 5 {-arbitrarily-} $ Data.Numbers.Primes.primes	-- Excludes even numbers, provided at least the 1st prime is tested. 		)		= False 		| otherwise	= ( 			case algorithm of@@ -114,14 +110,9 @@ 	Show					i  ) => factorisationAlgorithm -> i -> Bool isPrimeByAKS factorisationAlgorithm n	= and [-	not $ Math.PerfectPower.isPerfectPower n,	--Step 1.-	Math.Primality.areCoprime n `all` filter (/= n) [2 .. r],	--Step 3.-#if MIN_VERSION_parallel(3,0,0)-	and $ Control.Parallel.Strategies.parMap Control.Parallel.Strategies.rdeepseq	--Benefits from '+RTS -H100M', which reduces garbage-collections.-#else-	all-#endif-	(+	not $ Math.PerfectPower.isPerfectPower n,	-- Step 1.+	Math.Primality.areCoprime n `all` filter (/= n) [2 .. r],	-- Step 3.+	and $ Control.Parallel.Strategies.parMap Control.Parallel.Strategies.rdeepseq	{-Benefits from '+RTS -H100M', which reduces garbage-collections-} ( 		\a	-> let --			lhs, rhs :: Data.Polynomial.Polynomial i i 			lhs	= Data.Polynomial.raiseModulo (Data.Polynomial.mkLinear 1 a) n {-power-} n {-modulus-}@@ -135,7 +126,7 @@ 		) -- Because all these polynomials are /monic/, one can establish /congruence/ using /integer/-division. 	) [ 		1 .. floor . (* lg) . sqrt $ fromIntegral r-	] --Step 4; (x + a)^n ~ x^n + a mod (x^r - 1, n).+	] -- Step 4; (x + a)^n ~ x^n + a mod (x^r - 1, n).  ] where 	lg :: Double 	lg	= logBase 2 $ fromIntegral n@@ -145,7 +136,7 @@ 		(<= floor (Math.Power.square lg)) . snd 	 ) . map ( 		id &&& Math.MultiplicativeOrder.multiplicativeOrder factorisationAlgorithm n-	 ) $ Math.Primality.areCoprime n `filter` [2 ..]	--Step 2.+	 ) $ Math.Primality.areCoprime n `filter` [2 ..]	-- Step 2.  --	modulus :: Data.Polynomial.Polynomial i i 	modulus	= Data.Polynomial.mkPolynomial [(1, r), (negate 1, 0)]@@ -167,10 +158,10 @@ 	-> i	-- ^ Base. 	-> Bool witnessesCompositeness candidate oddRemainder nPowersOfTwo base	= all (-	$ ((`rem` candidate) . Math.Power.square) `iterate` Math.Power.raiseModulo base oddRemainder candidate	--Repeatedly modulo-square.+	$ ((`rem` candidate) . Math.Power.square) `iterate` Math.Power.raiseModulo base oddRemainder candidate	-- Repeatedly modulo-square.  ) [-	(/= 1) . head,					--Check whether the zeroeth modulo-power is incongruent to one.-	notElem (pred candidate) . take nPowersOfTwo	--Check whether any modulo-power is incongruent to -1.+	(/= 1) . head,					-- Check whether the zeroeth modulo-power is incongruent to one.+	notElem (pred candidate) . take nPowersOfTwo	-- Check whether any modulo-power is incongruent to -1.  ]  {- |@@ -193,12 +184,12 @@ -} isPrimeByMillerRabin :: (Integral i, Show i) => i -> Bool isPrimeByMillerRabin primeCandidate	= not $ witnessesCompositeness primeCandidate (-	fst $ last binaryFactors	--Odd-remainder.+	fst $ last binaryFactors	-- Odd-remainder.  ) (-	length binaryFactors	--The number of times that 'two' can be factored-out from 'predecessor'.+	length binaryFactors	-- The number of times that 'two' can be factored-out from 'predecessor'.  ) `any` testBases	where 	predecessor	= pred primeCandidate-	binaryFactors	= takeWhile ((== 0) . snd) . tail {-drop the original-} $ iterate ((`quotRem` 2) . fst) (predecessor, 0)	--Factor-out powers of two.+	binaryFactors	= takeWhile ((== 0) . snd) . tail {-drop the original-} $ iterate ((`quotRem` 2) . fst) (predecessor, 0)	-- Factor-out powers of two. 	testBases 		| null fewestPrimeBases	= let 			millersTestSet	= floor . (* 2 {-Eric Bach-}) . Math.Power.square . toRational {-avoid premature rounding-} $ log (fromIntegral primeCandidate :: Double {-overflows at 10^851-})@@ -206,21 +197,21 @@ 		| otherwise		= head fewestPrimeBases `take` Data.Numbers.Primes.primes 		where 			fewestPrimeBases	= map fst $ dropWhile ((primeCandidate >=) . snd) [-				(0,	9),			--All odd integers less this, are prime, and require no further verification.+				(0,	9),			-- All odd integers less this, are prime, and require no further verification. 				(1,	2047), 				(2,	1373653), 				(3,	25326001), 				(4,	3215031751),-				(5,	2152302898747),		--Jaeschke ...+				(5,	2152302898747),		-- Jaeschke ... 				(6,	3474749660383), 				(8,	341550071728321),-				(11,	3825123056546413051),	--Zhang ...+				(11,	3825123056546413051),	-- Zhang ... 				(12,	318665857834031151167461), 				(13,	3317044064679887385961981), 				(14,	6003094289670105800312596501), 				(15,	59276361075595573263446330101), 				(17,	564132928021909221014087501701), 				(19,	1543267864443420616877677640751301),-				(20,	10 ^ (36 :: Int))	--At least.+				(20,	10 ^ (36 :: Int))	-- At least. 			 ] 

@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {- 	Copyright (C) 2011 Dr. Alistair Ward @@ -42,6 +41,7 @@ import			Control.Arrow((&&&)) import qualified	Control.Arrow import qualified	Control.DeepSeq+import qualified	Control.Parallel.Strategies import qualified	Data.Maybe import qualified	Data.Numbers.Primes import qualified	Factory.Data.Exponential	as Data.Exponential@@ -53,10 +53,6 @@ import qualified	ToolShed.Data.Pair import qualified	ToolShed.Defaultable -#if MIN_VERSION_parallel(3,0,0)-import qualified	Control.Parallel.Strategies-#endif- -- | The algorithms by which prime-factorisation has been implemented. data Algorithm 	= DixonsMethod	-- ^ <http://en.wikipedia.org/wiki/Dixon%27s_factorization_method>.@@ -106,14 +102,10 @@ 	| i <= 3				= [Data.Exponential.rightIdentity i] 	| even i				= Data.Exponential.rightIdentity 2 : factoriseByFermatsMethod (i `div` 2) {-recurse-} 	| Data.Maybe.isJust maybeSquareNumber	= (<^ 2) `map` factoriseByFermatsMethod (Data.Maybe.fromJust maybeSquareNumber) {-recurse-}-	| null factors				= [Data.Exponential.rightIdentity i]	--Prime.-	| otherwise				= uncurry (++) .-#if MIN_VERSION_parallel(3,0,0)-	Control.Parallel.Strategies.withStrategy (-		Control.Parallel.Strategies.parTuple2 Control.Parallel.Strategies.rdeepseq Control.Parallel.Strategies.rdeepseq	--CAVEAT: unproductive on the size of integers tested so far.-	) .-#endif-	ToolShed.Data.Pair.mirror factoriseByFermatsMethod $ head factors+	| null factors				= [Data.Exponential.rightIdentity i]	-- Prime.+	| otherwise				= uncurry (++) . Control.Parallel.Strategies.withStrategy (+		Control.Parallel.Strategies.parTuple2 Control.Parallel.Strategies.rdeepseq Control.Parallel.Strategies.rdeepseq	-- CAVEAT: unproductive on the size of integers tested so far.+	) . ToolShed.Data.Pair.mirror factoriseByFermatsMethod $ head factors 	where --		maybeSquareNumber :: Integral i => Maybe i 		maybeSquareNumber	= Math.PerfectPower.maybeSquareNumber i@@ -121,15 +113,15 @@ --		factors :: Integral i => [i] 		factors	= map ( 			(-				uncurry (+) &&& uncurry (-)	--Construct the co-factors as the sum and difference of /larger/ and /smaller/.+				uncurry (+) &&& uncurry (-)	-- Construct the co-factors as the sum and difference of /larger/ and /smaller/. 			) . Control.Arrow.second Data.Maybe.fromJust 		 ) . filter (-			Data.Maybe.isJust . snd	--Search for a perfect square.+			Data.Maybe.isJust . snd	-- Search for a perfect square. 		 ) . map (-			Control.Arrow.second $ Math.PerfectPower.maybeSquareNumber {-hotspot-} . (+ negate i)	--Associate the corresponding value of /smaller/.+			Control.Arrow.second $ Math.PerfectPower.maybeSquareNumber {-hotspot-} . (+ negate i)	-- Associate the corresponding value of /smaller/. 		 ) . takeWhile (-			(<= (i + 9) `div` 6) . fst	--Terminate the search at the maximum value of /larger/.-		 ) . Math.Power.squaresFrom {-hotspot-} . ceiling $ sqrt (fromIntegral i :: Double)	--Start the search at the minimum value of /larger/.+			(<= (i + 9) `div` 6) . fst	-- Terminate the search at the maximum value of /larger/.+		 ) . Math.Power.squaresFrom {-hotspot-} . ceiling $ sqrt (fromIntegral i :: Double)	-- Start the search at the minimum value of /larger/.  {- | 	* Decomposes the specified integer, into a product of /prime/-factors,

@@ -53,11 +53,11 @@ 	deriving (Eq, Read, Show)  instance ToolShed.Defaultable.Defaultable Algorithm	where-	defaultValue	= SieveOfEratosthenes 7	--Resulting in a wheel of circumference 510510.+	defaultValue	= SieveOfEratosthenes 7	-- Resulting in a wheel of circumference 510510.  instance Math.Primes.Algorithmic Algorithm	where 	primes (SieveOfAtkin maxPrime)		= Math.Implementations.Primes.SieveOfAtkin.sieveOfAtkin (Data.PrimeWheel.estimateOptimalSize maxPrime) $ fromIntegral maxPrime 	primes (SieveOfEratosthenes wheelSize)	= Math.Implementations.Primes.SieveOfEratosthenes.sieveOfEratosthenes wheelSize 	primes (TrialDivision wheelSize)	= Math.Implementations.Primes.TrialDivision.trialDivision wheelSize 	primes TurnersSieve			= Math.Implementations.Primes.TurnersSieve.turnersSieve-	primes (WheelSieve wheelSize)		= Data.Numbers.Primes.wheelSieve wheelSize	--Has better space-complexity than 'SieveOfEratosthenes'.+	primes (WheelSieve wheelSize)		= Data.Numbers.Primes.wheelSieve wheelSize	-- Has better space-complexity than 'SieveOfEratosthenes'.

@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {- 	Copyright (C) 2011 Dr. Alistair Ward @@ -53,6 +52,7 @@ ) where  import qualified	Control.DeepSeq+import qualified	Control.Parallel.Strategies import qualified	Data.Array.IArray import			Data.Array.IArray((!)) import qualified	Data.IntSet@@ -62,10 +62,6 @@ import qualified	Factory.Math.Power	as Math.Power import qualified	ToolShed.Data.List -#if MIN_VERSION_parallel(3,0,0)-import qualified	Control.Parallel.Strategies-#endif- -- | Defines the types of /quadratic/, available to test the potential primality of a candidate integer. data PolynomialType 	= ModFour	-- ^ Suitable for primality-testing numbers meeting @(n `mod` 4 == 1)@.@@ -76,7 +72,7 @@  -- | The constant modulus used to select the appropriate quadratic for a prime candidate. atkinsModulus :: Integral i => i-atkinsModulus	= foldr1 lcm [4, 6, 12]	--Sure, this is always '12', but this is the reason why.+atkinsModulus	= foldr1 lcm [4, 6, 12]	-- Sure, this is always '12', but this is the reason why.  -- | The constant list of primes factored-out by the unoptimised algorithm. inherentPrimes :: Integral i => [i]@@ -113,11 +109,11 @@ --	select :: Integral i => i -> PolynomialType 	select n 		| any (-			(== 0) . (n `rem`)		--Though this is merely /Trial Division/, it's only performed over a short bounded interval of numerators.+			(== 0) . (n `rem`)		-- Though this is merely /Trial Division/, it's only performed over a short bounded interval of numerators. 		) primeComponents	= None-		| r `elem` [1, 5]	= ModFour	--We actually require @(n `mod` 4 == 1)@, but this is the equivalent modulo 12, with @(r == 9)@ removed because they're all divisible by /3/.-		| r == 7		= ModSix	--We actually require @(n `mod` 6 == 1)@, but this is the equivalent modulo 12, where @(r == 1)@ has been accounted for above.-		| r == 11		= ModTwelve	--We require @(n `mod` 12 == 11)@.+		| r `elem` [1, 5]	= ModFour	-- We actually require @(n `mod` 4 == 1)@, but this is the equivalent modulo 12, with @(r == 9)@ removed because they're all divisible by /3/.+		| r == 7		= ModSix	-- We actually require @(n `mod` 6 == 1)@, but this is the equivalent modulo 12, where @(r == 1)@ has been accounted for above.+		| r == 11		= ModTwelve	-- We require @(n `mod` 12 == 11)@. 		| otherwise		= None 		where 			r		= n `rem` atkinsModulus@@ -139,7 +135,7 @@ 	/GHC/'s old /qsort/-implementation is even slower :( -} filterOddRepetitions :: Ord a => [a] -> [a]---filterOddRepetitions	= map head . filter (foldr (const not) False) . Data.List.group . Data.List.sort	--Too slow.+-- filterOddRepetitions	= map head . filter (foldr (const not) False) . Data.List.group . Data.List.sort	-- Too slow. filterOddRepetitions	= slave True . Data.List.sort where 	slave isOdd (one : remainder@(two : _)) 		| one == two	= slave (not isOdd) remainder@@ -166,38 +162,36 @@ 	=> Data.PrimeWheel.PrimeWheel i 	-> i	-- ^ The maximum prime-number required. 	-> [i]-findPolynomialSolutions primeWheel maxPrime	= foldr1 ToolShed.Data.List.merge --The lists were previously sorted, as a side-effect, by 'filterOddRepetitions'.-#if MIN_VERSION_parallel(3,0,0)-	$ Control.Parallel.Strategies.withStrategy (Control.Parallel.Strategies.parList Control.Parallel.Strategies.rdeepseq)-#endif-	[+findPolynomialSolutions primeWheel maxPrime	= foldr1 ToolShed.Data.List.merge {-The lists were previously sorted, as a side-effect, by 'filterOddRepetitions'-} $ Control.Parallel.Strategies.withStrategy (+		Control.Parallel.Strategies.parList Control.Parallel.Strategies.rdeepseq+	 ) [ 		{-# SCC "4x^2+y^2" #-} filterOddRepetitions [ 			z | 				x'	<- takeWhile (<= pred maxPrime) $ map (* 4) squares, 				z	<- takeWhile (<= maxPrime) $ map (+ x') oddSquares, 				lookupPolynomialType z == ModFour-		], --List-comprehension. Twice the length of the other two lists.+		], -- List-comprehension. Twice the length of the other two lists. 		{-# SCC "3x^2+y^2" #-} filterOddRepetitions [ 			z | 				x'	<- takeWhile (<= pred maxPrime) $ map (* 3) squares, 				z	<- takeWhile (<= maxPrime) . map (+ x') $ if even x' then oddSelection else evenSelection, 				lookupPolynomialType z == ModSix-		], --List-comprehension.+		], -- List-comprehension. 		{-# SCC "3x^2-y^2" #-} filterOddRepetitions [ 			z | 				x2	<- takeWhile (<= maxPrime `div` 2) squares, 				z	<- dropWhile (> maxPrime) . map (3 * x2 -) . takeWhile (< x2) $ if even x2 then oddSelection else evenSelection, 				lookupPolynomialType z == ModTwelve-		] --List-comprehension.+		] -- List-comprehension. 	] where 		(evenSquares, oddSquares)	= Data.List.partition even squares  --		evenSelection, oddSelection :: Integral i => [i] 		evenSelection	= selection110 evenSquares	where-			selection110 (x0 : x1 : _ : xs)	= x0 : x1 : selection110 xs	--Effectively, those for meeting ((== 4) . (`mod` 6)).+			selection110 (x0 : x1 : _ : xs)	= x0 : x1 : selection110 xs	-- Effectively, those for meeting ((== 4) . (`mod` 6)). 			selection110 xs			= xs 		oddSelection	= selection101 oddSquares	where-			selection101 (x0 : _ : x2 : xs)	= x0 : x2 : selection101 xs	--Effectively, those for meeting ((== 1) . (`mod` 6)).+			selection101 (x0 : _ : x2 : xs)	= x0 : x2 : selection101 xs	-- Effectively, those for meeting ((== 1) . (`mod` 6)). 			selection101 xs			= xs  --		lookupPolynomialType :: (Data.Array.IArray.Ix i, Integral i) => i -> PolynomialType@@ -228,11 +222,11 @@ --	filterSquareFree :: Integral i => Data.Set.Set i -> [i] -> [i] 	filterSquareFree _ []	= [] 	filterSquareFree primeMultiples (candidate : candidates)-		| Data.Set.member candidate primeMultiples	= {-# SCC "delete" #-} filterSquareFree (Data.Set.delete candidate primeMultiples) candidates	--Tail-recurse.+		| Data.Set.member candidate primeMultiples	= {-# SCC "delete" #-} filterSquareFree (Data.Set.delete candidate primeMultiples) candidates	-- Tail-recurse. 		| otherwise					= {-# SCC "insert" #-} candidate : filterSquareFree (Data.Set.union primeMultiples . Data.Set.fromDistinctAscList $ generateMultiplesOfSquareTo primeWheel candidate maxPrime) candidates  {-# NOINLINE sieveOfAtkin #-}-{-# RULES "sieveOfAtkin/Int" sieveOfAtkin = sieveOfAtkinInt #-}	--CAVEAT: doesn't fire when built with profiling enabled.+{-# RULES "sieveOfAtkin/Int" sieveOfAtkin = sieveOfAtkinInt #-}	-- CAVEAT: doesn't fire when built with profiling enabled.  -- | A specialisation of 'sieveOfAtkin', which reduces both the execution-time and the space required. sieveOfAtkinInt :: Data.PrimeWheel.NPrimes -> Int -> [Int]

@@ -107,15 +107,15 @@  	sieve :: Integral i => Data.PrimeWheel.Distance i -> Repository i -> [i] 	sieve distance@(candidate, rollingWheel) repository@(primeSquares, squareFreePrimeMultiples)	= case Data.Map.lookup candidate squareFreePrimeMultiples of-		Just primeMultiples	-> sieve' $ Control.Arrow.second (insertUniq primeMultiples . Data.Map.delete candidate) repository	--Re-insert subsequent multiples.-		Nothing --Not a square-free composite.-			| candidate == smallestPrimeSquare	-> sieve' $ (tail' *** insertUniq subsequentPrimeMultiples) repository	--Migrate subsequent prime-multiples, from 'primeSquares' to 'squareFreePrimeMultiples'.+		Just primeMultiples	-> sieve' $ Control.Arrow.second (insertUniq primeMultiples . Data.Map.delete candidate) repository	-- Re-insert subsequent multiples.+		Nothing -- Not a square-free composite.+			| candidate == smallestPrimeSquare	-> sieve' $ (tail' *** insertUniq subsequentPrimeMultiples) repository	-- Migrate subsequent prime-multiples, from 'primeSquares' to 'squareFreePrimeMultiples'. 			| otherwise {-prime-}			-> candidate : sieve' (Control.Arrow.first (|> Data.PrimeWheel.generateMultiples candidate rollingWheel) repository) 			where 				(smallestPrimeSquare : subsequentPrimeMultiples)	= head' primeSquares 		where --			sieve' :: Repository i -> [i]-			sieve'	= sieve $ Data.PrimeWheel.rotate distance	--Tail-recurse.+			sieve'	= sieve $ Data.PrimeWheel.rotate distance	-- Tail-recurse.  			insertUniq :: Ord i => Data.PrimeWheel.PrimeMultiples i -> PrimeMultiplesMap i -> PrimeMultiplesMap i 			insertUniq l m	= insert $ dropWhile (`Data.Map.member` m) l	where@@ -124,7 +124,7 @@ 				insert (key : values)	= Data.Map.insert key values m  {-# NOINLINE sieveOfEratosthenes #-}-{-# RULES "sieveOfEratosthenes/Int" sieveOfEratosthenes = sieveOfEratosthenesInt #-}	--CAVEAT: doesn't fire when built with profiling enabled.+{-# RULES "sieveOfEratosthenes/Int" sieveOfEratosthenes = sieveOfEratosthenesInt #-}	-- CAVEAT: doesn't fire when built with profiling enabled.  -- | A specialisation of 'PrimeMultiplesMap'. type PrimeMultiplesMapInt	= Data.IntMap.IntMap (Data.PrimeWheel.PrimeMultiples Int)

@@ -40,8 +40,6 @@ 	-> Bool isIndivisibleBy numerator	= all ((/= 0) . (numerator `rem`)) . takeWhile (<= Math.PrimeFactorisation.maxBoundPrimeFactor numerator) -{-# INLINE isIndivisibleBy #-}- {- | 	* For each candidate, confirm indivisibility, by all /primes/ smaller than its /square-root/. @@ -49,13 +47,13 @@ 	of parameterised, but static, size; <http://en.wikipedia.org/wiki/Wheel_factorization>. -} trialDivision :: Integral prime => Data.PrimeWheel.NPrimes -> [prime]-trialDivision 0	= [2, 3] ++ filter (`isIndivisibleBy` trialDivision 0 {-recurse-}) [5 ..]	--No faster than using 'Data.PrimeWheel.mkPrimeWheel 0', but apparently better space-complexity ?!+trialDivision 0	= [2, 3] ++ filter (`isIndivisibleBy` trialDivision 0 {-recurse-}) [5 ..]	-- No faster than using 'Data.PrimeWheel.mkPrimeWheel 0', but apparently better space-complexity ?! trialDivision wheelSize	= Data.PrimeWheel.getPrimeComponents primeWheel ++ indivisible	where 	primeWheel	= Data.PrimeWheel.mkPrimeWheel wheelSize 	candidates	= map fst $ Data.PrimeWheel.roll primeWheel 	indivisible	= uncurry (++) . Control.Arrow.second ( 		filter (`isIndivisibleBy` indivisible {-recurse-}) 	 ) $ Data.List.span (-		< Math.Power.square (head candidates)	--The first composite candidate, is the square of the next prime after the wheel's constituent ones.+		< Math.Power.square (head candidates)	-- The first composite candidate, is the square of the next prime after the wheel's constituent ones. 	 ) candidates 

@@ -40,8 +40,8 @@ 	sieve (prime : candidates)	= prime : sieve ( 		filter ( 			\candidate	-> any ($ candidate) [-				(< Math.Power.square prime),	--Unconditionally admit any candidate smaller than the square of the last prime.-				(/= 0) . (`rem` prime)		--Ensure indivisibility, of all subsequent candidates, by the last prime discovered.+				(< Math.Power.square prime),	-- Unconditionally admit any candidate smaller than the square of the last prime.+				(/= 0) . (`rem` prime)		-- Ensure indivisibility, of all subsequent candidates, by the last prime discovered. 			] 		) candidates 	 )

@@ -85,14 +85,14 @@  instance Math.SquareRoot.Iterator Algorithm where 	step BakhshaliApproximation y x-		| dy == 0	= x		--The estimate was precise.-		| otherwise	= x' - dx'	--Correct the estimate.+		| dy == 0	= x		-- The estimate was precise.+		| otherwise	= x' - dx'	-- Correct the estimate. 		where 			dy, dydx, dx, x', dydx', dx' :: Math.SquareRoot.Result 			dy	= Math.SquareRoot.getDiscrepancy y x 			dydx	= 2 * x 			dx	= dy / dydx-			x'	= x + dx	--Identical to Newton-Raphson iteration.+			x'	= x + dx	-- Identical to Newton-Raphson iteration. 			dydx'	= 2 * x' 			dx'	= Math.Power.square dx / dydx' @@ -104,17 +104,17 @@ >			=> Xn - 1 / [2Xn / (Xn^2 - Y) - 1 / 2Xn] -} 	step HalleysMethod y x-		| dy == 0	= x		--The estimate was precise.-		| otherwise	= x - dx	--Correct the estimate.+		| dy == 0	= x		-- The estimate was precise.+		| otherwise	= x - dx	-- Correct the estimate. 		where 			dy, dydx, dx :: Math.SquareRoot.Result-			dy	= negate $ Math.SquareRoot.getDiscrepancy y x	--Use the estimate to determine the error in 'y'.-			dydx	= 2 * x						--The gradient, at the estimated value 'x'.+			dy	= negate $ Math.SquareRoot.getDiscrepancy y x	-- Use the estimate to determine the error in 'y'.+			dydx	= 2 * x						-- The gradient, at the estimated value 'x'. 			dx	= recip $ dydx / dy - recip dydx ---	step NewtonRaphsonIteration y x	= (x + toRational y / x) / 2		--This is identical to the /Babylonian Method/.---	step NewtonRaphsonIteration y x	= x / 2 + toRational y / (2 * x)	--Faster.-	step NewtonRaphsonIteration y x	= x / 2 + (toRational y / 2) / x	--Faster still.+--	step NewtonRaphsonIteration y x	= (x + toRational y / x) / 2		-- This is identical to the /Babylonian Method/.+--	step NewtonRaphsonIteration y x	= x / 2 + toRational y / (2 * x)	-- Faster.+	step NewtonRaphsonIteration y x	= x / 2 + (toRational y / 2) / x	-- Faster still.  	step (TaylorSeries terms) y x	= squareRootByTaylorSeries terms y x @@ -124,7 +124,7 @@ 	convergenceOrder ContinuedFraction	= Math.Precision.linearConvergence 	convergenceOrder HalleysMethod		= Math.Precision.cubicConvergence 	convergenceOrder NewtonRaphsonIteration	= Math.Precision.quadraticConvergence-	convergenceOrder (TaylorSeries terms)	= terms	--The order of convergence, per iteration, equals the number of terms in the series on each iteration.+	convergenceOrder (TaylorSeries terms)	= terms	-- The order of convergence, per iteration, equals the number of terms in the series on each iteration.  {- | 	* Uses /continued-fractions/, to iterate towards the principal /square-root/ of the specified positive integer;@@ -184,7 +184,7 @@ 	| otherwise	= Math.Summation.sumR' . take terms . zipWith (*) taylorSeriesCoefficients $ iterate (* relativeError) x 	where 		relativeError :: Math.SquareRoot.Result-		relativeError	= pred $ toRational y / Math.Power.square x	--Pedantically, this is the error in y, which is twice the magnitude of the error in x.+		relativeError	= pred $ toRational y / Math.Power.square x	-- Pedantically, this is the error in y, which is twice the magnitude of the error in x.  -- | Iterates from the estimated value, towards the /square-root/, a sufficient number of times to achieve the required accuracy. squareRootByIteration :: Real operand => Algorithm -> ProblemSpecification operand

@@ -50,7 +50,7 @@ multiplicativeOrder primeFactorisationAlgorithm base modulus 	| modulus < 2					= error $ "Factory.Math.MultiplicativeOrder.multiplicativeOrder:\tinvalid modulus; " ++ show modulus 	| not $ Math.Primality.areCoprime base modulus	= error $ "Factory.Math.MultiplicativeOrder.multiplicativeOrder:\targuments aren't coprime; " ++ show (base, modulus)-	| otherwise					= foldr (lcm . multiplicativeOrder') 1 $ Math.PrimeFactorisation.primeFactors primeFactorisationAlgorithm modulus	--Combine the /multiplicative order/ of the constituent /prime-factors/.+	| otherwise					= foldr (lcm . multiplicativeOrder') 1 $ Math.PrimeFactorisation.primeFactors primeFactorisationAlgorithm modulus	-- Combine the /multiplicative order/ of the constituent /prime-factors/. 	where --		multiplicativeOrder' :: (Control.DeepSeq.NFData i, Integral i) => Data.Exponential.Exponential i -> i 		multiplicativeOrder' e	= product . map (

@@ -44,18 +44,18 @@ -} maybeSquareNumber :: Integral i => i -> Maybe i maybeSquareNumber i---	| i < 0					= Nothing	--This function is performance-sensitive, but this test is neither strictly nor frequently required.+--	| i < 0					= Nothing	-- This function is performance-sensitive, but this test is neither strictly nor frequently required. 	| all (\(modulus, valid) -> rem i modulus `elem` valid) [---							--Distribution of moduli amongst perfect squares	Cumulative failure-detection.-		(16,	[0,1,4,9]),			--All moduli are equally likely.			75%-		(9,	[0,1,4,7]),			--Zero occurs 33%, the others only 22%.			88%-		(17,	[1,2,4,8,9,13,15,16,0]),	--Zero only occurs 5.8%, the others 11.8%.		94%+--							-- Distribution of moduli amongst perfect squares	Cumulative failure-detection.+		(16,	[0,1,4,9]),			-- All moduli are equally likely.			75%+		(9,	[0,1,4,7]),			-- Zero occurs 33%, the others only 22%.			88%+		(17,	[1,2,4,8,9,13,15,16,0]),	-- Zero only occurs 5.8%, the others 11.8%.		94% -- These additional tests, aren't always cost-effective.-		(13,	[1,3,4,9,10,12,0]),		--Zero only occurs 7.7%, the others 15.4%.		97%-		(7,	[1,2,4,0]),			--Zero only occurs 14.3%, the others 28.6%.		98%-		(5,	[1,4,0])			--Zero only occurs 20%, the others 40%.			99%+		(13,	[1,3,4,9,10,12,0]),		-- Zero only occurs 7.7%, the others 15.4%.		97%+		(7,	[1,2,4,0]),			-- Zero only occurs 14.3%, the others 28.6%.		98%+		(5,	[1,4,0])			-- Zero only occurs 20%, the others 40%.			99% ---	] && fromIntegral iSqrt == sqrt'	= Just iSqrt	--CAVEAT: erroneously True for 187598574531033120 (187598574531033121 is square).+--	] && fromIntegral iSqrt == sqrt'	= Just iSqrt	-- CAVEAT: erroneously True for 187598574531033120 (187598574531033121 is square). 	] && Math.Power.square iSqrt == i	= Just iSqrt 	| otherwise				= Nothing 	where@@ -82,7 +82,7 @@ 		\n set	-> if n `Data.Set.member` set 			then set --			else Data.Set.union set . Data.Set.fromDistinctAscList . takeWhile (<= i) . iterate (* n) $ Math.Power.square n-			else foldr Data.Set.insert set . takeWhile (<= i) . iterate (* n) $ Math.Power.square n	--Faster.+			else foldr Data.Set.insert set . takeWhile (<= i) . iterate (* n) $ Math.Power.square n	-- Faster. 	) Data.Set.empty [2 .. round $ sqrt (fromIntegral i :: Double)]  {-# NOINLINE isPerfectPower #-}

@@ -53,7 +53,7 @@ 	openS algorithm decimalDigits 		| decimalDigits <= 0	= "" 		| decimalDigits <= 16	= take (succ decimalDigits) $ show (pi :: Double)-		| otherwise		= "3." ++ tail (show $ openI algorithm decimalDigits)	--Insert a decimal point.+		| otherwise		= "3." ++ tail (show $ openI algorithm decimalDigits)	-- Insert a decimal point.  -- | Categorises the various algorithms. data Category agm bbp borwein ramanujan spigot

@@ -31,7 +31,7 @@  -- | Mainly for convenience. square :: Num n => n -> n-square	= (^ (2 :: Int))+square x	= x ^ (2 :: Int)	-- CAVEAT: this could be eta-reduced, but it won't then inline when called with a single argument.  {-# INLINE square #-} @@ -71,7 +71,7 @@ raiseModulo _ _ 1	= 0 raiseModulo _ 0 modulus	= 1 `mod` modulus raiseModulo base power modulus-	| base < 0		= (`mod` modulus) . (if even power then id else negate) $ raiseModulo (negate base) power modulus	--Recurse.+	| base < 0		= (`mod` modulus) . (if even power then id else negate) $ raiseModulo (negate base) power modulus	-- Recurse. 	| power < 0		= error $ "Factory.Math.Power.raiseModulo:\tnegative power; " ++ show power 	| first `elem` [0, 1]	= first 	| otherwise		= slave power

@@ -121,5 +121,5 @@ 	=> DecimalDigits	-- ^ The number of places after the decimal point, which are required. 	-> operand 	-> Rational-simplify decimalDigits operand	= Data.Ratio.approxRational operand . recip $ 4 * 10 ^ succ decimalDigits	--Tolerate any error less than half the least significant digit required.+simplify decimalDigits operand	= Data.Ratio.approxRational operand . recip $ 4 * 10 ^ succ decimalDigits	-- Tolerate any error less than half the least significant digit required. 

@@ -70,7 +70,7 @@ -} isFermatWitness :: (Integral i, Show i) => i -> Bool isFermatWitness i	= not . all isFermatPseudoPrime $ filter (areCoprime i) [2 .. pred i]	where-	isFermatPseudoPrime base	= Math.Power.raiseModulo base (pred i) i == 1	--CAVEAT: a /Fermat Pseudo-prime/ must also be a /composite/ number.+	isFermatPseudoPrime base	= Math.Power.raiseModulo base (pred i) i == 1	-- CAVEAT: a /Fermat Pseudo-prime/ must also be a /composite/ number.  {- | 	* A /Carmichael number/ is an /odd/ /composite/ number which satisfies /Fermat's little theorem/.

@@ -60,5 +60,5 @@ 	* <http://mathworld.wolfram.com/MersenneNumber.html> -} mersenneNumbers :: (Algorithmic algorithm, Integral i) => algorithm -> [i]-mersenneNumbers algorithm	= map (pred . (2 ^)) (primes algorithm :: [Int])	--Whilst the exponentiation could be parallelised, not all values are known to be required.+mersenneNumbers algorithm	= map (pred . (2 ^)) (primes algorithm :: [Int])	-- Whilst the exponentiation could be parallelised, not all values are known to be required. 

@@ -103,10 +103,10 @@ 	getMean :: Fractional mean => probabilityDistribution -> mean	-- ^ The theoretical mean.  	getStandardDeviation :: Floating standardDeviation => probabilityDistribution -> standardDeviation-- ^ The theoretical standard-deviation.-	getStandardDeviation	= sqrt . getVariance	--Default implementation.+	getStandardDeviation	= sqrt . getVariance	-- Default implementation.  	getVariance :: Floating variance => probabilityDistribution -> variance	-- ^ The theoretical variance.-	getVariance	= Math.Power.square . getStandardDeviation	--Default implementation.+	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@@ -117,7 +117,7 @@ 	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 (LogNormalDistribution location scale2)		= realToFrac $ (exp scale2 - 1) * exp (2 * location + scale2)	-- NB: for standard-deviation == mean, use scale^2 == ln 2. 	getVariance (NormalDistribution _ variance)			= realToFrac variance 	getVariance (UniformDistribution (minParameter, maxParameter))	= realToFrac $ Math.Power.square (maxParameter - minParameter) / 12 @@ -191,7 +191,7 @@ 				quantile	= (/ lambda) . negate . log . (1 -)	-- <http://en.wikipedia.org/wiki/Quantile_function>. 			 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.+				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@@ -219,13 +219,13 @@ generatePoissonDistribution 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.+		negate . log $ minPositiveFloat lambda	-- Guard against underflow, in the user-defined type for lambda. 	)		= filter (>= 0) . map round . (reProfile (PoissonDistribution lambda) :: [Double] -> [Double]) . generateStandardizedNormalDistribution 	| otherwise	= generator 	where 		generator	= uncurry (:) . ( 			fst . head . dropWhile (-				(> exp (negate lambda)) . snd	--CAVEAT: underflows if lambda > (103 :: Float, 745 :: Double).+				(> 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-}@@ -249,7 +249,7 @@ 		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.+				| 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 

@@ -65,8 +65,8 @@ 	Show			base, 	Show			decimal  ) => base -> decimal -> String-toBase 10 decimal	= show decimal	--Base unchanged.-toBase _ 0		= "0"		--Zero has the same representation in any base.+toBase 10 decimal	= show decimal	-- Base unchanged.+toBase _ 0		= "0"		-- Zero has the same representation in any base. toBase base decimal 	| abs base < 2					= error $ "Factory.Math.Radix.toBase:\tan arbitrary integer can't be represented in base " ++ show base 	| abs base > Data.List.genericLength digits	= error $ "Factory.Math.Radix.toBase:\tunable to clearly represent the complete set of digits in base " ++ show base@@ -75,7 +75,7 @@ 	where 		fromDecimal 0		= id 		fromDecimal n-			| remainder < 0	= fromDecimal (succ quotient) . ((remainder - fromIntegral base) :)	--This can only occur when base is negative; cf. 'divMod'.+			| remainder < 0	= fromDecimal (succ quotient) . ((remainder - fromIntegral base) :)	-- This can only occur when base is negative; cf. 'divMod'. 			| otherwise	= fromDecimal quotient . (remainder :) 			where 				(quotient, remainder)	= n `quotRem` fromIntegral base@@ -100,12 +100,12 @@ 	Read		decimal, 	Show		base  ) => base -> String -> decimal-fromBase 10 s	= read s	--Base unchanged.-fromBase _ "0"	= 0		--Zero has the same representation in any base.+fromBase 10 s	= read s	-- Base unchanged.+fromBase _ "0"	= 0		-- Zero has the same representation in any base. fromBase base s 	| abs base < 2					= error $ "Factory.Math.Radix.fromBase:\tan arbitrary integer can't be represented in base " ++ show base 	| abs base > Data.List.genericLength digits	= error $ "Factory.Math.Radix.fromBase:\tunable to clearly represent the complete set of digits in base " ++ show base-	| base > 0 && head s == '-'			= negate . fromBase base $ tail s	--Recurse.+	| base > 0 && head s == '-'			= negate . fromBase base $ tail s	-- Recurse. 	| otherwise					= Data.List.foldl' (\l -> ((l * fromIntegral base) +) . fromDigit) 0 s	where 		fromDigit :: Integral i => Char -> i 		fromDigit c	= case c `lookup` decodes of

@@ -64,7 +64,7 @@ 		-> Math.Precision.DecimalDigits	-- ^ The required precision. 		-> operand			-- ^ The value for which to find the /square-root/. 		-> Result			-- ^ Returns an estimate of the /square-root/, found using the specified algorithm, accurate to at least the required number of decimal digits.-	squareRoot algorithm decimalDigits operand	= squareRootFrom algorithm (getEstimate operand) decimalDigits operand	--Default implementation+	squareRoot algorithm decimalDigits operand	= squareRootFrom algorithm (getEstimate operand) decimalDigits operand	-- Default implementation  -- | The interface required to iterate, from an estimate of the required value, to the next approximation. class Iterator algorithm where@@ -111,10 +111,10 @@ -} getAccuracy :: Real operand => operand -> Result -> Math.Precision.DecimalDigits getAccuracy y x-	| absoluteError == 0	= maxBound	--Bodge.---	| otherwise		= length . takeWhile (< 1) $ iterate (* 10) relativeError	--CAVEAT: too slow.+	| absoluteError == 0	= maxBound	-- Bodge.+--	| otherwise		= length . takeWhile (< 1) $ iterate (* 10) relativeError	-- CAVEAT: too slow. 	| otherwise		= length $ show (round $ toRational y / absoluteError :: Integer) 	where 		absoluteError :: Result-		absoluteError	= abs (getDiscrepancy y x) / 2	--NB: the magnitude of the error in 'y', is twice the error in its square-root, 'x'.+		absoluteError	= abs (getDiscrepancy y x) / 2	-- NB: the magnitude of the error in 'y', is twice the error in its square-root, 'x'. 

@@ -23,6 +23,7 @@ module Factory.Math.Statistics( -- * Functions 	getMean,+	getWeightedMean, --	getDispersionFromMean, 	getVariance, 	getStandardDeviation,@@ -45,27 +46,58 @@  	* 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+getMean :: (+	Data.Foldable.Foldable	foldable,+	Fractional		result,+	Real			value+ )+	=> foldable value+	-> result+getMean foldable 	| denominator == 0	= error "Factory.Math.Statistics.getMean:\tno data => undefined result." 	| otherwise		= realToFrac numerator / fromIntegral denominator 	where-		(numerator, denominator)	= Data.Foldable.foldr (\s -> (+ s) *** succ) (0, 0 :: Int) x+		(numerator, denominator)	= Data.Foldable.foldr (\s -> (+ s) *** succ) (0, 0 :: Int) foldable  {- |+	* Determines the /weighted mean/ of the specified numbers; <http://en.wikipedia.org/wiki/Weighted_arithmetic_mean>.++	* The specified value is only evaluated if the corresponding weight is non-zero.++	* Should the caller define the result-type as 'Rational', then it will be free from rounding-errors.+-}+getWeightedMean :: (+	Data.Foldable.Foldable	foldable,+	Fractional		result,+	Real			value,+	Real			weight+ )+	=> foldable (value, weight)	-- ^ Each pair consists of a value & the corresponding weight.+	-> result+getWeightedMean foldable+	| denominator == 0	= error "Factory.Math.Statistics.getWeightedMean:\tzero weight => undefined result."+	| otherwise		= numerator / realToFrac denominator+	where+		(numerator, denominator)	= Data.Foldable.foldr (+			\(value, weight)	-> if weight == 0+				then id	--Avoid unnecessarily evaluation.+				else (+ realToFrac value * realToFrac weight) *** (+ weight)+		 ) (0, 0) foldable++{- | 	* Measures the /dispersion/ of a /population/ of results from the /mean/ value; <http://en.wikipedia.org/wiki/Statistical_dispersion>.  	* Should the caller define the result-type as 'Rational', then it will be free from rounding-errors. -} getDispersionFromMean :: (-	Data.Foldable.Foldable	f,+	Data.Foldable.Foldable	foldable, 	Fractional		result,-	Functor			f,-	Real			r- ) => (Rational -> Rational) -> f r -> result-getDispersionFromMean weight x	= getMean $ fmap (weight . (+ negate mean) . toRational) x	where+	Functor			foldable,+	Real			value+ ) => (Rational -> Rational) -> foldable value -> result+getDispersionFromMean weight foldable	= getMean $ fmap (weight . (+ negate mean) . toRational) foldable	where 	mean :: Rational-	mean	= getMean x+	mean	= getMean foldable  {- | 	* Determines the exact /variance/ of the specified numbers; <http://en.wikipedia.org/wiki/Variance>.@@ -73,20 +105,20 @@ 	* Should the caller define the result-type as 'Rational', then it will be free from rounding-errors. -} getVariance :: (-	Data.Foldable.Foldable	f,+	Data.Foldable.Foldable	foldable, 	Fractional		variance,-	Functor			f,-	Real			r- ) => f r -> variance+	Functor			foldable,+	Real			value+ ) => foldable value -> variance getVariance	= getDispersionFromMean Math.Power.square  -- | Determines the /standard-deviation/ of the specified numbers; <http://en.wikipedia.org/wiki/Standard_deviation>. getStandardDeviation :: (-	Data.Foldable.Foldable	f,+	Data.Foldable.Foldable	foldable, 	Floating		result,-	Functor			f,-	Real			r- ) => f r -> result+	Functor			foldable,+	Real			value+ ) => foldable value -> result getStandardDeviation	= sqrt . getVariance  {- |@@ -95,21 +127,21 @@ 	* Should the caller define the result-type as 'Rational', then it will be free from rounding-errors. -} getAverageAbsoluteDeviation :: (-	Data.Foldable.Foldable	f,+	Data.Foldable.Foldable	foldable, 	Fractional		result,-	Functor			f,-	Real			r- ) => f r -> result+	Functor			foldable,+	Real			value+ ) => foldable value -> result getAverageAbsoluteDeviation	= getDispersionFromMean abs  -- | Determines the /coefficient-of-variance/ of the specified numbers; <http://en.wikipedia.org/wiki/Coefficient_of_variation>. getCoefficientOfVariance :: (-	Data.Foldable.Foldable	f,+	Data.Foldable.Foldable	foldable, 	Eq			result, 	Floating		result,-	Functor			f,-	Real			r- ) => f r -> result+	Functor			foldable,+	Real			value+ ) => foldable value -> result getCoefficientOfVariance l 	| mean == 0	= error "Factory.Math.Statistics.getCoefficientOfVariance:\tundefined if mean is zero." 	| otherwise	= getStandardDeviation l / abs mean

@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {- 	Copyright (C) 2011 Dr. Alistair Ward @@ -29,15 +28,12 @@ ) where  import qualified	Control.DeepSeq+import qualified	Control.Parallel.Strategies import qualified	Data.List import qualified	Data.Ratio import			Data.Ratio((%)) import qualified	ToolShed.Data.List -#if MIN_VERSION_parallel(3,0,0)-import qualified	Control.Parallel.Strategies-#endif- {- | 	* Sums a list of numbers of arbitrary type. @@ -52,7 +48,6 @@ 	=> ToolShed.Data.List.ChunkLength 	-> [n] 	-> n-#if MIN_VERSION_parallel(3,0,0) sum' chunkLength 	| chunkLength <= 1	= error $ "Factory.Math.Summation.sum':\tinvalid chunk-size; " ++ show chunkLength 	| otherwise		= slave@@ -61,20 +56,17 @@ 		slave []	= 0 		slave [x]	= x 		slave l		= slave {-recurse-} . Control.Parallel.Strategies.parMap Control.Parallel.Strategies.rdeepseq sum $ ToolShed.Data.List.chunk chunkLength l-#else-sum' _	= sum-#endif  {- | 	* Sums a list of /rational/ type numbers.  	* CAVEAT: though faster than 'Data.List.sum', this algorithm has poor space-complexity, making it unsuitable for unrestricted use. -}-{-# INLINE sumR' #-}	--This makes a staggering difference.+{-# INLINE sumR' #-}	-- This makes a staggering difference. sumR' :: Integral i => [Data.Ratio.Ratio i] -> Data.Ratio.Ratio i sumR' l	= foldr (\ratio -> ((Data.Ratio.numerator ratio * (commonDenominator `div` Data.Ratio.denominator ratio)) +)) 0 l % commonDenominator	where --	commonDenominator	= foldr (lcm . Data.Ratio.denominator) 1 l-	commonDenominator	= Data.List.foldl' (\multiple -> lcm multiple . Data.Ratio.denominator) 1 l	--Slightly faster.+	commonDenominator	= Data.List.foldl' (\multiple -> lcm multiple . Data.Ratio.denominator) 1 l	-- Slightly faster.  {- | 	* Sums a list of /rational/ numbers.@@ -84,7 +76,7 @@  	* CAVEAT: memory-use is proportional to chunk-size. -}-{-# INLINE sumR #-}	--This makes a staggering difference to calls from other modules.+{-# INLINE sumR #-}	-- This makes a staggering difference to calls from other modules. sumR :: (Integral i, Control.DeepSeq.NFData i) 	=> ToolShed.Data.List.ChunkLength 	-> [Data.Ratio.Ratio i]@@ -96,10 +88,4 @@ 		slave :: (Integral i, Control.DeepSeq.NFData i) => [Data.Ratio.Ratio i] -> Data.Ratio.Ratio i 		slave l 			| length l <= chunkLength	= sumR' l-			| otherwise		= slave {-recurse-} .-#if MIN_VERSION_parallel(3,0,0)-				Control.Parallel.Strategies.parMap Control.Parallel.Strategies.rdeepseq-#else-				map-#endif-				sumR' $ ToolShed.Data.List.chunk chunkLength l+			| otherwise			= slave {-recurse-} . Control.Parallel.Strategies.parMap Control.Parallel.Strategies.rdeepseq sumR' $ ToolShed.Data.List.chunk chunkLength l

@@ -44,7 +44,7 @@  -- | Measures the CPU-time required by a naive implementation. factorialPerformanceControl :: (Control.DeepSeq.NFData i, Integral i) => i -> IO (Double, i)---factorialPerformanceControl i	= ToolShed.System.TimePure.getCPUSeconds $ product [1 .. i]	--CAVEAT: too lazy.+-- factorialPerformanceControl i	= ToolShed.System.TimePure.getCPUSeconds $ product [1 .. i]	-- CAVEAT: too lazy. factorialPerformanceControl i	= ToolShed.System.TimePure.getCPUSeconds $ Data.List.foldl' (*) 1 [2 .. i]  {- |

@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {- 	Copyright (C) 2011 Dr. Alistair Ward @@ -29,6 +28,7 @@ 	quickChecks ) where +import qualified	Data.Tuple import qualified	Factory.Math.ArithmeticGeometricMean	as Math.ArithmeticGeometricMean import qualified	Factory.Math.Implementations.SquareRoot	as Math.Implementations.SquareRoot import qualified	Factory.Math.Precision			as Math.Precision@@ -36,14 +36,6 @@ import qualified	Test.QuickCheck import			Test.QuickCheck((==>)) -#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- type Testable	= Math.Implementations.SquareRoot.Algorithm -> Math.Precision.DecimalDigits -> Math.ArithmeticGeometricMean.AGM -> Int -> Test.QuickCheck.Property  -- | Defines invariant properties.@@ -53,7 +45,7 @@ 	prop_symmetrical squareRootAlgorithm decimalDigits agm index	= Math.ArithmeticGeometricMean.isValid agm ==> Test.QuickCheck.label "prop_symmetrical" . and . tail . take index' $ zipWith (==) ( 		Math.ArithmeticGeometricMean.convergeToAGM squareRootAlgorithm decimalDigits' agm 	 ) (-		Math.ArithmeticGeometricMean.convergeToAGM squareRootAlgorithm decimalDigits' $ swap agm+		Math.ArithmeticGeometricMean.convergeToAGM squareRootAlgorithm decimalDigits' $ Data.Tuple.swap agm 	 ) where 		decimalDigits'	= succ $ decimalDigits `mod` 64 		index'		= succ $ index `mod` 8

@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {-# OPTIONS_GHC -fno-warn-orphans #-} {- 	Copyright (C) 2011 Dr. Alistair Ward@@ -39,9 +38,6 @@  instance Test.QuickCheck.Arbitrary Math.Implementations.Factorial.Algorithm	where 	arbitrary	= Test.QuickCheck.elements [Math.Implementations.Factorial.Bisection, Math.Implementations.Factorial.PrimeFactorisation]-#if !(MIN_VERSION_QuickCheck(2,1,0))-	coarbitrary	= undefined	--CAVEAT: stops warnings from ghc.-#endif  type Testable	= Integer -> Integer -> Test.QuickCheck.Property 

@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {-# OPTIONS_GHC -fno-warn-orphans #-} {- 	Copyright (C) 2011 Dr. Alistair Ward@@ -50,9 +49,6 @@ 		polynomial	<- Test.QuickCheck.arbitrary  		return {-to Gen-monad-} . Data.MonicPolynomial.mkMonicPolynomial $ ((1, succ $ Data.Polynomial.getDegree polynomial) :) `Data.Polynomial.lift` polynomial-#if !(MIN_VERSION_QuickCheck(2,1,0))-	coarbitrary	= undefined	--CAVEAT: stops warnings from ghc.-#endif  type P	= Data.MonicPolynomial.MonicPolynomial Integer Integer 

@@ -1,7 +1,6 @@-{-# LANGUAGE CPP #-} {-# OPTIONS_GHC -fno-warn-orphans #-} {--	Copyright (C) 2011 Dr. Alistair Ward+	Copyright (C) 2011-2015 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@@ -30,6 +29,7 @@ 	quickChecks ) where +import Prelude hiding ((<*>))	-- The "Prelude" from 'base-4.8' exports this symbol. import			Control.Applicative((<$>), (<*>)) import			Factory.Test.QuickCheck.Factorial() import			Factory.Test.QuickCheck.SquareRoot()@@ -52,15 +52,9 @@ 	Math.SquareRoot.Algorithmic	squareRootAlgorithm  ) => Test.QuickCheck.Arbitrary (Math.Implementations.Pi.AGM.Algorithm.Algorithm squareRootAlgorithm)	where 	arbitrary	= Math.Implementations.Pi.AGM.Algorithm.BrentSalamin <$> Test.QuickCheck.arbitrary-#if !(MIN_VERSION_QuickCheck(2,1,0))-	coarbitrary	= undefined	--CAVEAT: stops warnings from ghc.-#endif  instance Test.QuickCheck.Arbitrary Math.Implementations.Pi.BBP.Algorithm.Algorithm	where 	arbitrary	= Test.QuickCheck.elements [Math.Implementations.Pi.BBP.Algorithm.Bellard, Math.Implementations.Pi.BBP.Algorithm.Base65536]-#if !(MIN_VERSION_QuickCheck(2,1,0))-	coarbitrary	= undefined	--CAVEAT: stops warnings from ghc.-#endif  instance ( 	Test.QuickCheck.Arbitrary	squareRootAlgorithm,@@ -71,9 +65,6 @@ 	arbitrary	= Test.QuickCheck.oneof [ 		Math.Implementations.Pi.Borwein.Algorithm.Borwein1993 <$> Test.QuickCheck.arbitrary <*> Test.QuickCheck.arbitrary 	 ]-#if !(MIN_VERSION_QuickCheck(2,1,0))-	coarbitrary	= undefined	--CAVEAT: stops warnings from ghc.-#endif  instance ( 	Test.QuickCheck.Arbitrary	squareRootAlgorithm,@@ -85,15 +76,9 @@ 		Math.Implementations.Pi.Ramanujan.Algorithm.Classic <$> Test.QuickCheck.arbitrary <*> Test.QuickCheck.arbitrary, 		Math.Implementations.Pi.Ramanujan.Algorithm.Chudnovsky <$> Test.QuickCheck.arbitrary <*> Test.QuickCheck.arbitrary 	 ]-#if !(MIN_VERSION_QuickCheck(2,1,0))-	coarbitrary	= undefined	--CAVEAT: stops warnings from ghc.-#endif  instance Test.QuickCheck.Arbitrary Math.Implementations.Pi.Spigot.Algorithm.Algorithm	where 	arbitrary	= Test.QuickCheck.elements [Math.Implementations.Pi.Spigot.Algorithm.RabinowitzWagon, Math.Implementations.Pi.Spigot.Algorithm.Gosper]-#if !(MIN_VERSION_QuickCheck(2,1,0))-	coarbitrary	= undefined	--CAVEAT: stops warnings from ghc.-#endif  instance ( 	Test.QuickCheck.Arbitrary agm,@@ -109,9 +94,6 @@ 		Math.Pi.Ramanujan <$> Test.QuickCheck.arbitrary, 		Math.Pi.Spigot <$> Test.QuickCheck.arbitrary 	 ]-#if !(MIN_VERSION_QuickCheck(2,1,0))-	coarbitrary	= undefined	--CAVEAT: stops warnings from ghc.-#endif  type Category	= Math.Pi.Category ( 	Math.Implementations.Pi.AGM.Algorithm.Algorithm Math.Implementations.SquareRoot.Algorithm

@@ -1,7 +1,6 @@-{-# LANGUAGE CPP #-} {-# OPTIONS_GHC -fno-warn-orphans #-} {--	Copyright (C) 2011 Dr. Alistair Ward+	Copyright (C) 2011-2015 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@@ -27,7 +26,6 @@ 	quickChecks ) where -import			Control.Applicative((<$>)) import			Control.Arrow((***)) import			Factory.Data.Ring((=*=), (=+=), (=-=), (=^)) import qualified	Data.Numbers.Primes@@ -43,10 +41,7 @@ 	Test.QuickCheck.Arbitrary	e, 	Integral			e  ) => Test.QuickCheck.Arbitrary (Data.Polynomial.Polynomial c e)	where-	arbitrary	= Data.Polynomial.mkPolynomial . map ((+ negate 4) . (`mod` 8) *** (`mod` 8)) <$> Test.QuickCheck.arbitrary-#if !(MIN_VERSION_QuickCheck(2,1,0))-	coarbitrary	= undefined	--CAVEAT: stops warnings from ghc.-#endif+	arbitrary	= (Data.Polynomial.mkPolynomial . map ((+ negate 4) . (`mod` 8) *** (`mod` 8))) `fmap` Test.QuickCheck.arbitrary  -- | Defines invariant properties. quickChecks :: IO ()

@@ -1,7 +1,6 @@-{-# LANGUAGE CPP #-} {-# OPTIONS_GHC -fno-warn-orphans #-} {--	Copyright (C) 2011 Dr. Alistair Ward+	Copyright (C) 2011-2015 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@@ -27,7 +26,6 @@ 	quickChecks ) where -import			Control.Applicative((<$>)) import			Factory.Test.QuickCheck.PrimeFactorisation() import qualified	Data.List import qualified	Data.Numbers.Primes@@ -39,12 +37,9 @@  instance Test.QuickCheck.Arbitrary factorisationAlgorithm => Test.QuickCheck.Arbitrary (Math.Implementations.Primality.Algorithm factorisationAlgorithm)	where 	arbitrary	= Test.QuickCheck.oneof [-		Math.Implementations.Primality.AKS <$> Test.QuickCheck.arbitrary,+		Math.Implementations.Primality.AKS `fmap` Test.QuickCheck.arbitrary, 		return {-to Gen-monad-} Math.Implementations.Primality.MillerRabin 	 ]-#if !(MIN_VERSION_QuickCheck(2,1,0))-	coarbitrary	= undefined	--CAVEAT: stops warnings from ghc.-#endif  -- | Defines invariant properties. quickChecks :: IO ()@@ -56,7 +51,7 @@ 		prop_prime :: Math.Implementations.Primality.Algorithm Math.Implementations.PrimeFactorisation.Algorithm -> Integer -> Test.QuickCheck.Property 		prop_prime primalityAlgorithm i	= Test.QuickCheck.label "prop_prime" $ Math.Primality.isPrime primalityAlgorithm prime	where 			normalise n-				| primalityAlgorithm == Math.Implementations.Primality.MillerRabin	= n `mod` 1000000	--Limited by the efficiency of 'Data.Numbers.Primes.primes'.+				| primalityAlgorithm == Math.Implementations.Primality.MillerRabin	= n `mod` 1000000	-- Limited by the efficiency of 'Data.Numbers.Primes.primes'. 				| otherwise								= n `mod` 59  			prime :: Integer

@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {-# OPTIONS_GHC -fno-warn-orphans #-} {- 	Copyright (C) 2011 Dr. Alistair Ward@@ -44,9 +43,6 @@ 			Math.Implementations.PrimeFactorisation.FermatsMethod 		] 	 ]-#if !(MIN_VERSION_QuickCheck(2,1,0))-	coarbitrary	= undefined	--CAVEAT: stops warnings from ghc.-#endif  -- | Defines invariant properties. quickChecks :: IO ()@@ -69,7 +65,7 @@ 			i' :: Integer 			i'	= i `mod` 20 -		prop_eulersTotientP algorithm i	= Test.QuickCheck.label "prop_eulersTotient" $ Math.PrimeFactorisation.eulersTotient algorithm prime == pred prime	where+		prop_eulersTotientP algorithm i	= Test.QuickCheck.label "prop_eulersTotientP" $ Math.PrimeFactorisation.eulersTotient algorithm prime == pred prime	where 			prime :: Integer 			prime	= Data.List.genericIndex Data.Numbers.Primes.primes (i `mod` 10000) 

@@ -1,7 +1,6 @@-{-# LANGUAGE CPP #-} {-# OPTIONS_GHC -fno-warn-orphans #-} {--	Copyright (C) 2011 Dr. Alistair Ward+	Copyright (C) 2011-2015 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@@ -31,7 +30,6 @@ 	upperBound ) where -import			Control.Applicative((<$>)) import qualified	Control.DeepSeq import qualified	Data.Set import qualified	Factory.Data.PrimeWheel				as Data.PrimeWheel@@ -47,12 +45,9 @@ instance Test.QuickCheck.Arbitrary Math.Implementations.Primes.Algorithm.Algorithm	where 	arbitrary	= Test.QuickCheck.oneof [ 		return {-to Gen-monad-} Math.Implementations.Primes.Algorithm.TurnersSieve,-		Math.Implementations.Primes.Algorithm.TrialDivision . (`mod` 10) <$> Test.QuickCheck.arbitrary,-		Math.Implementations.Primes.Algorithm.SieveOfEratosthenes . (`mod` 10) <$> Test.QuickCheck.arbitrary+		(Math.Implementations.Primes.Algorithm.TrialDivision . (`mod` 10)) `fmap` Test.QuickCheck.arbitrary,+		(Math.Implementations.Primes.Algorithm.SieveOfEratosthenes . (`mod` 10)) `fmap` Test.QuickCheck.arbitrary 	 ]-#if !(MIN_VERSION_QuickCheck(2,1,0))-	coarbitrary	= undefined	--CAVEAT: stops warnings from ghc.-#endif  isPrime :: (Control.DeepSeq.NFData i, Integral i, Show i) => i -> Bool isPrime	= Math.Primality.isPrime primalityAlgorithm	where

@@ -53,7 +53,7 @@ quickChecks	= do 	randomGen	<- System.Random.getStdGen -	(Test.QuickCheck.quickCheck . ($ randomGen)) `mapM_` [+	Test.QuickCheck.quickCheck (prop_logNormalDistributionEqual randomGen) >> (Test.QuickCheck.quickCheck . ($ randomGen)) `mapM_` [ 		prop_logNormalDistribution, 		prop_logNormalDistribution', 		prop_normalDistribution,@@ -72,7 +72,7 @@  		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	-- Both of which, having been normalised, should be zero. 		 ) . ( 			normalise distribution :: [Double] -> [Double] 		 ) . take 10000 $ Math.Probability.generatePopulation distribution randomGen	where@@ -92,15 +92,30 @@ 				| location >= 0	= maxParameter `min` location 				| otherwise	= negate maxParameter `max` location +-- The mean & standard-deviation are equal when scale^2 == ln 2, but this seems to break-down when the mean is close to zero.+		prop_logNormalDistributionEqual :: System.Random.RandomGen randomGen => randomGen -> Double -> Test.QuickCheck.Property+		prop_logNormalDistributionEqual randomGen location	= location' > 1 ==> Test.QuickCheck.label "prop_logNormalDistributionEqual" . (+			< (recip 1000000 :: Double)+		 ) . pred . abs . uncurry (/) . (+			Math.Statistics.getMean &&& Math.Statistics.getStandardDeviation+		 ) $ take 10000 (+			Math.Probability.generatePopulation (Math.Probability.LogNormalDistribution location' $ log 2) randomGen :: [Double]+		 ) 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	-- Both of which, having been normalised, should be zero. 		 ) . ( 			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.+			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@@ -109,7 +124,7 @@  		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.+			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@@ -120,7 +135,7 @@ 		 ) . 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.+			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@@ -131,15 +146,15 @@ 		 ) . 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.+			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].+			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].+			probability'	= recip . succ $ abs probability	-- Semi-closed unit-interval (0, 1]. 

@@ -1,7 +1,6 @@-{-# LANGUAGE CPP #-} {-# OPTIONS_GHC -fno-warn-orphans #-} {--	Copyright (C) 2011 Dr. Alistair Ward+	Copyright (C) 2011-2015 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@@ -30,7 +29,6 @@ 	quickChecks ) where -import			Control.Applicative((<$>)) import			Data.Ratio((%)) import qualified	Data.Ratio import qualified	Factory.Math.Implementations.SquareRoot	as Math.Implementations.SquareRoot@@ -47,11 +45,8 @@ 			Math.Implementations.SquareRoot.HalleysMethod, 			Math.Implementations.SquareRoot.NewtonRaphsonIteration 		],-		Math.Implementations.SquareRoot.TaylorSeries <$> Test.QuickCheck.elements [2 .. 32]+		Math.Implementations.SquareRoot.TaylorSeries `fmap` Test.QuickCheck.elements [2 .. 32] 	 ]-#if !(MIN_VERSION_QuickCheck(2,1,0))-	coarbitrary	= undefined	--CAVEAT: stops warnings from ghc.-#endif  type Testable	= (Math.Implementations.SquareRoot.Algorithm, Math.Precision.DecimalDigits, Rational) -> Test.QuickCheck.Property @@ -67,7 +62,7 @@ 		operand'	= abs operand  	prop_factorable (algorithm, decimalDigits, operand)	= Test.QuickCheck.label "prop_factorable" . (<= 5) . (-		* 10 ^ requiredDecimalDigits	--Promote the relative error.+		* 10 ^ requiredDecimalDigits	-- Promote the relative error. 	 ) . abs $ 1 - ( 		Math.SquareRoot.squareRoot algorithm requiredDecimalDigits ( 			toRational $ Data.Ratio.numerator operand'@@ -86,6 +81,6 @@ 		requiredDecimalDigits	= succ $ decimalDigits `mod` 32768  		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'.+		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' 

@@ -1,5 +1,5 @@ {--	Copyright (C) 2011 Dr. Alistair Ward+	Copyright (C) 2011-2014 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@@ -44,6 +44,9 @@ 	>> Test.QuickCheck.quickCheck `mapM_` [prop_nP0, prop_nP1] 	>> Test.QuickCheck.quickCheck `mapM_` [prop_zeroVariance, prop_zeroAverageAbsoluteDeviation] 	>> Test.QuickCheck.quickCheck `mapM_` [prop_balance, prop_varianceRelocated, prop_varianceScaled, prop_varianceOrder, prop_equivalence, prop_varianceOfArray, prop_varianceOfMap, prop_meanOfSet]+	>> Test.QuickCheck.quickCheck prop_weightedMeanRational+	>> Test.QuickCheck.quickCheck prop_weightedMeanInteger+	>> Test.QuickCheck.quickCheck prop_weightedMeanUniformDenominator  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@@ -60,7 +63,7 @@  	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.+		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@@ -73,7 +76,7 @@ 	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 :: Rational)) l+	prop_balance l			= not (null l)	==> Test.QuickCheck.label "prop_balance" . (== 0) . abs . sum $ map (\i -> toRational i - Math.Statistics.getMean l) 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)@@ -82,4 +85,28 @@ 	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++	prop_weightedMeanRational :: [(Rational, Rational)] -> Test.QuickCheck.Property+	prop_weightedMeanRational assoc	= (denominator /= 0) ==> Test.QuickCheck.label "prop_weightedMeanRational" $ Math.Statistics.getWeightedMean assoc == (+		sum (map (uncurry (*)) assoc) / denominator+	 ) where+		denominator	= sum $ map snd assoc+++	prop_weightedMeanInteger :: [(Integer, Integer)] -> Test.QuickCheck.Property+	prop_weightedMeanInteger assoc	= (denominator /= 0) ==> Test.QuickCheck.label "prop_weightedMeanInteger" $ Math.Statistics.getWeightedMean assoc == (+		toRational (+			sum $ map (+				uncurry (*)+			) assoc+		) / toRational denominator+	 ) where+		denominator	= sum $ map snd assoc++	prop_weightedMeanUniformDenominator :: [Rational] -> Integer -> Test.QuickCheck.Property+	prop_weightedMeanUniformDenominator numerators i	= (not (null numerators) && i /= 0) ==> Test.QuickCheck.label "prop_weightedMeanUniformDenominator" $ Math.Statistics.getWeightedMean (+		zip numerators $ repeat i+	 ) == (+		Math.Statistics.getMean numerators :: Rational+	 ) 

@@ -49,7 +49,7 @@ import qualified	Factory.Test.Performance.SquareRoot		as Test.Performance.SquareRoot import qualified	Factory.Test.Performance.Statistics		as Test.Performance.Statistics import qualified	Factory.Test.QuickCheck.QuickChecks		as Test.QuickCheck.QuickChecks-import qualified	Paths_factory					as Paths	--Either local stub, or package-instance autogenerated by 'Setup.hs build'.+import qualified	Paths_factory					as Paths	-- Either local stub, or package-instance autogenerated by 'Setup.hs build'. import qualified	System.Console.GetOpt				as G import qualified	System.Environment import qualified	System.Exit@@ -63,7 +63,7 @@ type PiCategory			= Test.Performance.Pi.Category Math.Implementations.SquareRoot.Algorithm Math.Implementations.Factorial.Algorithm  -- | Used to thread user-defined command-line options, though the list of functions which implement them.-type CommandLineAction	= Test.CommandOptions.CommandOptions -> IO Test.CommandOptions.CommandOptions	--Supplied as the type-argument to 'G.OptDescr'.+type CommandLineAction	= Test.CommandOptions.CommandOptions -> IO Test.CommandOptions.CommandOptions	-- Supplied as the type-argument to 'G.OptDescr'.  -- | On failure to parse the specified string, returns an explanatory error. read' :: Read a => String -> String -> a@@ -78,7 +78,7 @@ -- | 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.+	System.IO.hClose System.IO.stdin	-- Nothing is read from standard input.  	progName	<- System.Environment.getProgName @@ -115,10 +115,10 @@ 			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-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+			printVersion	= System.IO.hPutStrLn System.IO.stderr (Distribution.Text.display packageIdentifier ++ "\n\nCopyright (C) 2011-2015 " ++ 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.+					Distribution.Package.pkgName	= Distribution.Package.PackageName progName,	-- CAVEAT: coincidentally. 					Distribution.Package.pkgVersion	= Distribution.Version.Version (Data.Version.versionBranch Paths.version) [] 				} @@ -216,7 +216,7 @@ 				) index :: IO ( 					Double, --					Integer-					Int	--Exploits rewrite-rules in "Math.Implementations.Primes.*".+					Int	-- Exploits rewrite-rules in "Math.Implementations.Primes.*". 				) 			 ) >>= print >> System.Exit.exitWith System.Exit.ExitSuccess	where 				algorithm :: Math.Implementations.Primes.Algorithm.Algorithm@@ -237,5 +237,5 @@ --	G.getOpt :: G.ArgOrder CommandLineAction -> [G.OptDescr Action] -> [String] -> ([Action], [String], [String]) 	case G.getOpt G.RequireOrder optDescrList args of 		(commandLineActions, _, [])	-> Data.List.foldl' (>>=) (return {-to IO-monad-} ToolShed.Defaultable.defaultValue) commandLineActions	>> System.Exit.exitWith System.Exit.ExitSuccess-		(_, _, errors)			-> System.IO.Error.ioError . System.IO.Error.userError $ concat errors ++ usageMessage	--Throw.+		(_, _, errors)			-> System.IO.Error.ioError . System.IO.Error.userError $ concat errors ++ usageMessage	-- Throw.