packages feed

factory 0.1.0.3 → 0.2.0.0

raw patch · 49 files changed

+938/−452 lines, 49 files

Files

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