{-
	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 <https://www.gnu.org/licenses/>.
-}
{- |
 [@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.

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

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

	* <https://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.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 (Algorithm -> Algorithm -> Bool
(Algorithm -> Algorithm -> Bool)
-> (Algorithm -> Algorithm -> Bool) -> Eq Algorithm
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Algorithm -> Algorithm -> Bool
$c/= :: Algorithm -> Algorithm -> Bool
== :: Algorithm -> Algorithm -> Bool
$c== :: Algorithm -> Algorithm -> Bool
Eq, ReadPrec [Algorithm]
ReadPrec Algorithm
Int -> ReadS Algorithm
ReadS [Algorithm]
(Int -> ReadS Algorithm)
-> ReadS [Algorithm]
-> ReadPrec Algorithm
-> ReadPrec [Algorithm]
-> Read Algorithm
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [Algorithm]
$creadListPrec :: ReadPrec [Algorithm]
readPrec :: ReadPrec Algorithm
$creadPrec :: ReadPrec Algorithm
readList :: ReadS [Algorithm]
$creadList :: ReadS [Algorithm]
readsPrec :: Int -> ReadS Algorithm
$creadsPrec :: Int -> ReadS Algorithm
Read, Int -> Algorithm -> ShowS
[Algorithm] -> ShowS
Algorithm -> String
(Int -> Algorithm -> ShowS)
-> (Algorithm -> String)
-> ([Algorithm] -> ShowS)
-> Show Algorithm
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Algorithm] -> ShowS
$cshowList :: [Algorithm] -> ShowS
show :: Algorithm -> String
$cshow :: Algorithm -> String
showsPrec :: Int -> Algorithm -> ShowS
$cshowsPrec :: Int -> Algorithm -> ShowS
Show)

instance Data.Default.Default Algorithm	where
	def :: Algorithm
def	= Algorithm
Bisection

instance Math.Factorial.Algorithmic Algorithm	where
	factorial :: Algorithm -> i -> i
factorial Algorithm
algorithm i
n
		| i
n i -> i -> Bool
forall a. Ord a => a -> a -> Bool
< i
2		= i
1
		| Bool
otherwise	= case Algorithm
algorithm of
			Algorithm
Bisection		-> i -> i -> i
forall i. (Integral i, Show i) => i -> i -> i
risingFactorial i
2 (i -> i) -> i -> i
forall a b. (a -> b) -> a -> b
$ i -> i
forall a. Enum a => a -> a
pred i
n
			Algorithm
PrimeFactorisation	-> BisectionRatio -> Int -> Factors i i -> i
forall base exponent.
(Num base, Integral exponent) =>
BisectionRatio -> Int -> Factors base exponent -> base
Data.PrimeFactors.product' (BisectionRatio -> BisectionRatio
forall a. Fractional a => a -> a
recip BisectionRatio
5) {-empirical-} Int
10 {-empirical-} (Factors i i -> i) -> Factors i i -> i
forall a b. (a -> b) -> a -> b
$ i -> Factors i i
forall base. Integral base => base -> Factors base base
primeFactors i
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.

	* <https://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 :: base -> Factors base base
primeFactors base
n	= ((base, base) -> Bool) -> Factors base base -> Factors base base
forall a. (a -> Bool) -> [a] -> [a]
takeWhile ((base -> base -> Bool
forall a. Ord a => a -> a -> Bool
> base
0) (base -> Bool) -> ((base, base) -> base) -> (base, base) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (base, base) -> base
forall a b. (a, b) -> b
snd) (Factors base base -> Factors base base)
-> Factors base base -> Factors base base
forall a b. (a -> b) -> a -> b
$ (base -> (base, base)) -> [base] -> Factors base base
forall a b. (a -> b) -> [a] -> [b]
map (\base
prime -> (base
prime, base -> base -> base
forall i. Integral i => i -> i -> i
primeMultiplicity base
prime base
n)) [base]
forall int. Integral int => [int]
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/.

	* <https://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 :: i -> i -> i
primeMultiplicity i
prime	= [i] -> i
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([i] -> i) -> (i -> [i]) -> i -> i
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (i -> Bool) -> [i] -> [i]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (i -> i -> Bool
forall a. Ord a => a -> a -> Bool
> i
0) ([i] -> [i]) -> (i -> [i]) -> i -> [i]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [i] -> [i]
forall a. [a] -> [a]
tail ([i] -> [i]) -> (i -> [i]) -> i -> [i]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (i -> i) -> i -> [i]
forall a. (a -> a) -> a -> [a]
iterate (i -> i -> i
forall i. Integral i => i -> i -> i
`div` i
prime)

-- | Returns the /rising factorial/; <https://mathworld.wolfram.com/RisingFactorial.html>
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 :: i -> i -> i
risingFactorial i
_ i
0	= i
1
risingFactorial i
0 i
_	= i
0
risingFactorial i
x i
n	= Ratio i -> i -> Interval i -> i
forall i. (Integral i, Show i) => Ratio i -> i -> Interval i -> i
Data.Interval.product' (Ratio i -> Ratio i
forall a. Fractional a => a -> a
recip Ratio i
2) i
64 (Interval i -> i) -> Interval i -> i
forall a b. (a -> b) -> a -> b
$ Interval i -> Interval i
forall endPoint.
Ord endPoint =>
Interval endPoint -> Interval endPoint
Data.Interval.normalise (i
x, i -> i
forall a. Enum a => a -> a
pred (i -> i) -> i -> i
forall a b. (a -> b) -> a -> b
$ i
x i -> i -> i
forall a. Num a => a -> a -> a
+ i
n)

-- | Returns the /falling factorial/; <https://mathworld.wolfram.com/FallingFactorial.html>
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 :: i -> i -> i
fallingFactorial i
_ i
0	= i
1
fallingFactorial i
0 i
_	= i
0
fallingFactorial i
x i
n	= Ratio i -> i -> Interval i -> i
forall i. (Integral i, Show i) => Ratio i -> i -> Interval i -> i
Data.Interval.product' (Ratio i -> Ratio i
forall a. Fractional a => a -> a
recip Ratio i
2) i
64 (Interval i -> i) -> Interval i -> i
forall a b. (a -> b) -> a -> b
$ Interval i -> Interval i
forall endPoint.
Ord endPoint =>
Interval endPoint -> Interval endPoint
Data.Interval.normalise (i
x, i -> i
forall a. Enum a => a -> a
succ (i -> i) -> i -> i
forall a b. (a -> b) -> a -> b
$ i
x i -> i -> i
forall a. Num a => a -> a -> a
- i
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.
i
numerator !/! :: i -> i -> f
!/! i
denominator
	| i
numerator i -> i -> Bool
forall a. Ord a => a -> a -> Bool
<= i
1		= f -> f
forall a. Fractional a => a -> a
recip (f -> f) -> (i -> f) -> i -> f
forall b c a. (b -> c) -> (a -> b) -> a -> c
. i -> f
forall a b. (Integral a, Num b) => a -> b
fromIntegral (i -> f) -> i -> f
forall a b. (a -> b) -> a -> b
$ Algorithm -> i -> i
forall algorithm i.
(Algorithmic algorithm, Integral i, Show i) =>
algorithm -> i -> i
Math.Factorial.factorial (Algorithm
forall a. Default a => a
Data.Default.def :: Algorithm) i
denominator
	| i
denominator i -> i -> Bool
forall a. Ord a => a -> a -> Bool
<= i
1		= i -> f
forall a b. (Integral a, Num b) => a -> b
fromIntegral (i -> f) -> i -> f
forall a b. (a -> b) -> a -> b
$ Algorithm -> i -> i
forall algorithm i.
(Algorithmic algorithm, Integral i, Show i) =>
algorithm -> i -> i
Math.Factorial.factorial (Algorithm
forall a. Default a => a
Data.Default.def :: Algorithm) i
numerator
	| i
numerator i -> i -> Bool
forall a. Eq a => a -> a -> Bool
== i
denominator	= f
1
	| i
numerator i -> i -> Bool
forall a. Ord a => a -> a -> Bool
< i
denominator	= f -> f
forall a. Fractional a => a -> a
recip (f -> f) -> f -> f
forall a b. (a -> b) -> a -> b
$ i
denominator i -> i -> f
forall i f. (Integral i, Fractional f, Show i) => i -> i -> f
!/! i
numerator	-- Recurse.
	| Bool
otherwise			= i -> f
forall a b. (Integral a, Num b) => a -> b
fromIntegral (i -> f) -> i -> f
forall a b. (a -> b) -> a -> b
$ Ratio i -> i -> Interval i -> i
forall i. (Integral i, Show i) => Ratio i -> i -> Interval i -> i
Data.Interval.product' (Ratio i -> Ratio i
forall a. Fractional a => a -> a
recip Ratio i
2) i
64 (i -> i
forall a. Enum a => a -> a
succ i
denominator, i
numerator)