{-
Copyright (C) 2011-2015 Dr. Alistair Ward
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see .
-}
{- |
[@AUTHOR@] Dr. Alistair Ward
[@DESCRIPTION@]
* Provides implementations of the class 'Math.Factorial.Algorithmic'.
* Provides additional functions related to /factorials/, but which depends on a specific implementation,
and which therefore can't be accessed throught the class-interface.
* .
* .
* .
-}
module Factory.Math.Implementations.Factorial(
-- * Types
-- ** Data-types
Algorithm(..),
-- * Functions
primeFactors,
-- primeMultiplicity,
risingFactorial,
fallingFactorial,
-- ** Operators
(!/!)
) where
import qualified Data.Default
import qualified Data.Numbers.Primes
import qualified Factory.Data.Interval as Data.Interval
import qualified Factory.Data.PrimeFactors as Data.PrimeFactors
import qualified Factory.Math.Factorial as Math.Factorial
infixl 7 !/! -- Same as (/).
-- | The algorithms by which /factorial/ has been implemented.
data Algorithm
= Bisection -- ^ The integers from which the /factorial/ is composed, are multiplied using @Data.Interval.product'@.
| PrimeFactorisation -- ^ The /prime factors/ of the /factorial/ are extracted, then raised to the appropriate power, before multiplication.
deriving (Eq, Read, Show)
instance Data.Default.Default Algorithm where
def = Bisection
instance Math.Factorial.Algorithmic Algorithm where
factorial algorithm n
| n < 2 = 1
| otherwise = case algorithm of
Bisection -> risingFactorial 2 $ pred n
PrimeFactorisation -> Data.PrimeFactors.product' (recip 5) {-empirical-} 10 {-empirical-} $ primeFactors n
{- |
* Returns the /prime factors/, of the /factorial/ of the specifed integer.
* Precisely all the primes less than or equal to the specified integer /n/, are included in /n!/;
only the multiplicity of each of these known prime components need be determined.
* .
* CAVEAT: currently a hotspot.
-}
primeFactors :: Integral base
=> base -- ^ The number, whose /factorial/ is to be factorised.
-> Data.PrimeFactors.Factors base base -- ^ The /base/ and /exponent/ of each /prime factor/ in the /factorial/, ordered by increasing /base/ (and decreasing /exponent/).
primeFactors n = takeWhile ((> 0) . snd) $ map (\prime -> (prime, primeMultiplicity prime n)) Data.Numbers.Primes.primes
{- |
* The number of times a specific /prime/, can be factored from the /factorial/ of the specified integer.
* General purpose /prime-factorisation/ has /exponential time-complexity/,
so use /Legendre's Theorem/, which relates only to the /prime factors/ of /factorials/.
* .
-}
primeMultiplicity :: Integral i
=> i -- ^ A prime number.
-> i -- ^ The integer, the factorial of which the prime is a factor.
-> i -- ^ The number of times the prime occurs in the factorial.
primeMultiplicity prime = sum . takeWhile (> 0) . tail . iterate (`div` prime)
-- | Returns the /rising factorial/;
risingFactorial :: (Integral i, Show i)
=> i -- ^ The lower bound of the integer-range, whose product is returned.
-> i -- ^ The number of integers in the range above.
-> i -- ^ The result.
risingFactorial _ 0 = 1
risingFactorial 0 _ = 0
risingFactorial x n = Data.Interval.product' (recip 2) 64 $ Data.Interval.normalise (x, pred $ x + n)
-- | Returns the /falling factorial/;
fallingFactorial :: (Integral i, Show i)
=> i -- ^ The upper bound of the integer-range, whose product is returned.
-> i -- ^ The number of integers in the range beneath.
-> i -- ^ The result.
fallingFactorial _ 0 = 1
fallingFactorial 0 _ = 0
fallingFactorial x n = Data.Interval.product' (recip 2) 64 $ Data.Interval.normalise (x, succ $ x - n)
{- |
* Returns the ratio of two factorials.
* It is more efficient than evaluating both factorials, and then dividing.
* For more complex combinations of factorials, such as in the /Binomial coefficient/,
extract the /prime factors/ using 'primeFactors'
then manipulate them using the module "Data.PrimeFactors",
and evaluate it using by /Data.PrimeFactors.product'/.
-}
(!/!) :: (Integral i, Fractional f, Show i)
=> i -- ^ The /numerator/.
-> i -- ^ The /denominator/.
-> f -- ^ The resulting fraction.
numerator !/! denominator
| numerator <= 1 = recip . fromIntegral $ Math.Factorial.factorial (Data.Default.def :: Algorithm) denominator
| denominator <= 1 = fromIntegral $ Math.Factorial.factorial (Data.Default.def :: Algorithm) numerator
| numerator == denominator = 1
| numerator < denominator = recip $ denominator !/! numerator -- Recurse.
| otherwise = fromIntegral $ Data.Interval.product' (recip 2) 64 (succ denominator, numerator)