```{-
Copyright (C) 2011 Dr. Alistair Ward

This program is free software: you can redistribute it and/or modify
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@]

* Provides implementations of the class 'Math.Factorial.Algorithm'.

* Provides additional functions related to /factorials/, but which depends on a specific implementation,
and which therefore can't be accessed throught the class-interface.

* <http://en.wikipedia.org/wiki/Factorial>.

* <http://mathworld.wolfram.com/Factorial.html>.

* <http://www.luschny.de/math/factorial/FastFactorialFunctions.htm>.
-}

module Factory.Math.Implementations.Factorial(
-- * Types
-- ** Data-types
Algorithm(..),
-- * Functions
primeFactors,
--	primeMultiplicity,
risingFactorial,
fallingFactorial,
-- ** Operators
(!/!)
) where

import qualified	Data.Numbers.Primes
import qualified	Factory.Data.Bounds		as Data.Bounds
import qualified	Factory.Data.PrimeFactors	as Data.PrimeFactors
import qualified	Factory.Math.Factorial		as Math.Factorial
import qualified	ToolShed.Defaultable		as Defaultable

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.Bounds.product'@.
| PrimeFactorisation	-- ^ The /prime factors/ of the /factorial/ are extracted, then raised to the appropriate power, before multiplication.

instance Defaultable.Defaultable Algorithm	where
defaultValue	= Bisection

instance Math.Factorial.Algorithm Algorithm	where
factorial algorithm n
| n < 2		= 1
| otherwise	= case algorithm of
Bisection		-> risingFactorial 2 \$ n - 1
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.

* <http://en.wikipedia.org/wiki/Factorial#Number_theory>.

* 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/.

* <http://www.proofwiki.org/wiki/Multiplicity_of_Prime_Factor_in_Factorial>.
-}
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/; <http://mathworld.wolfram.com/RisingFactorial.html>
risingFactorial :: Integral 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.Bounds.product' (recip 2) 64 \$ Data.Bounds.normalise (x, (x + n) - 1)

-- | Returns the /falling factorial/; <http://mathworld.wolfram.com/FallingFactorial.html>
fallingFactorial :: Integral 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.Bounds.product' (recip 2) 64 \$ Data.Bounds.normalise (x, (x - n) + 1)

{- |
* 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)
=> i	-- ^ The /numerator/.
-> i	-- ^ The /denominator/.
-> f	-- ^ The resulting fraction.
numerator !/! denominator
| numerator <= 1		= recip . fromIntegral \$ Math.Factorial.factorial (Defaultable.defaultValue :: Algorithm) denominator
| denominator <= 1		= fromIntegral \$ Math.Factorial.factorial (Defaultable.defaultValue :: Algorithm) numerator
| numerator == denominator	= 1
| numerator < denominator	= recip \$ denominator !/! numerator	--Recurse.
| otherwise			= fromIntegral \$ Data.Bounds.product' (recip 2) 64 (succ denominator, numerator)

```