{-# OPTIONS_GHC -Wall -fwarn-tabs #-}
{-# LANGUAGE CPP #-}
----------------------------------------------------------------
--                                                    2012.02.02
-- |
-- Module      :  Math.Combinatorics.Exact.Factorial
-- Stability   :  experimental
-- Portability :  Haskell98 + CPP
--
-- The factorial numbers (<http://oeis.org/A000142>). For negative
-- inputs, all functions return 0 (rather than throwing an exception
-- or using 'Maybe').
--
-- Notable limits:
--
-- * 12! is the largest factorial that can fit into 'Int32'.
--
-- * 20! is the largest factorial that can fit into 'Int64'.
--
-- * 170! is the largest factorial that can fit into 64-bit 'Double'.
----------------------------------------------------------------
module Math.Combinatorics.Exact.Factorial (factorial) where

-- N.B., we need a Custom cabal build-type for this to work.
import Data.Int  (Int32, Int64)
#endif
import Data.Bits

{-

fallingFactorial n k = product [n - fromInteger i | i <- [0..toInteger k - 1] ]
-- == factorial n div factorial (n-k)

risingFactorial n k = product [n + fromInteger i | i <- [0..toInteger k - 1] ]
-- == factorial (n+k) div factorial n

-- | A common under-approximation of the factorial numbers.
factorial_stirling :: (Integral a) => a -> a
{-# SPECIALIZE factorial_stirling ::
Integer -> Integer,
Int     -> Int,
Int32   -> Int32,
Int64   -> Int64
#-}
factorial_stirling n
| n < 0     = 0
| otherwise = ceiling (sqrt (2 * pi * n') * (n' / exp 1) ** n')
where
n' :: Double
n' = fromIntegral n
-}

----------------------------------------------------------------
{-
n!  = 2^{n - popCount n}
* \prod_{k \geq 1} \left(
\prod_{n/2^k < j \leq 2*n/2^k}
if odd j then j else 1
\right)^k
-}

-- | Exact factorial numbers. For a fast approximation see
-- naive definition of the factorial numbers is:
--
-- > factorial n
-- >     | n < 0     = 0
-- >     | otherwise = product [1..n]
--
-- However, we use a fast algorithm based on the split-recursive form:
--
-- > factorial n =
-- >     2^(n - popCount n) * product [(q k)^k | forall k, k >= 1]
-- >     where
-- >     q k = product [j | forall j, n*2^(-k) < j <= n*2^(-k+1), odd j]
--
factorial :: (Integral a, Bits a) => Int -> a
factorial n
| n < 0     = 0
| n < 2     = 1
| otherwise = go (highestBitPosition_Int n - 1) 0 0 1 1 1 1
where
-- lo  == n/2^(k+1)
-- lo' == n/2^k
-- qk  == product of odd @j@s for @k@ in [1..K]
-- p   == q1 * q2 * ... * qK
-- r   == (q1 ^ K) * (q2 ^ (K-1)) * ... * (qK ^ 1)
-- s   == 2^{n - popCount n}
-- go :: Int -> Int -> Int -> Int -> a -> a -> a -> a
go k lo s hi j p r
| k seq lo seq s seq hi seq j seq p seq r seq False = undefined
| k >= 0 =                     -- TODO: why did old version use lo/=n ?
let lo' = n shiftR k     -- TODO: use shiftRL#
hi' = (lo' - 1) .|. 1  -- if odd lo' then lo' else lo' - 1
len = (hi' - hi) div 2 -- TODO: why not (shiftR1) or (quot2) ?
in if len > 0
then let
(q, j') = partialProduct len j
p' = p * q
r' = r * p'
in go (k - 1) lo' (s + lo) hi' j' p' r'
else   go (k - 1) lo' (s + lo) hi' j  p  r
--
-- fromIntegral s /= fromIntegral n - popCount (fromIntegral n) = error "factorial_splitRecursive: bug in the computation of n - popCount n"
| otherwise = r shiftL s

-- | The product of odd @j@s between n/2^k and 2*n/2^k. @len@
-- is the count of @j@ terms to multiply, where the @j@ state
-- argument is the largest previously used term.
partialProduct :: (Integral a) => Int -> a -> (a,a)
partialProduct len j
| half == 0 = (,) <!>  (j+2)        <!> (j+2)
| len  == 2 = (,) <!> ((j+2)*(j+4)) <!> (j+4)
| otherwise =
let (qL, j' ) = partialProduct (len - half) j
(qR, j'') = partialProduct half         j'
in (,) <!> (qL*qR) <!> j''
where
half = len quot 2

(<!>) = ($!) -- fix associativity {- floorLog2 :: (Integral a, Bits a) => a -> Int floorLog2 n | n <= 0 = error "floorLog2: argument must be positive" | otherwise = highestBitPosition n - 1 highestBitPosition :: (Integral a, Bits a) => a -> Int {-# INLINE highestBitPosition #-} {-# SPECIALIZE highestBitPosition :: Int -> Int #-} highestBitPosition n0 | n0 < 0 = error _highestBitPosition_negative | n0 == 0 = 1 | otherwise = go 0 n0 where go d n | d seq n seq False = undefined | n > 0 = go (d+1) (n shiftR 1) | otherwise = d _highestBitPosition_negative :: String {-# NOINLINE _highestBitPosition_negative #-} _highestBitPosition_negative = "highestBitPosition: argument must be non-negative" floorLog2_Int :: Int -> Int floorLog2_Int n | n <= 0 = error "floorLog2_Int: argument must be positive" | otherwise = highestBitPosition_Int n - 1 -} highestBitPosition_Int :: Int -> Int highestBitPosition_Int w = if w < 1 shiftL 15 then if w < 1 shiftL 7 then if w < 1 shiftL 3 then if w < 1 shiftL 1 then if w < 1 shiftL 0 then if w < 0 then 32 else 0 -- N.B., Int semantics else 1 else if w < 1 shiftL 2 then 2 else 3 else if w < 1 shiftL 5 then if w < 1 shiftL 4 then 4 else 5 else if w < 1 shiftL 6 then 6 else 7 else if w < 1 shiftL 11 then if w < 1 shiftL 9 then if w < 1 shiftL 8 then 8 else 9 else if w < 1 shiftL 10 then 10 else 11 else if w < 1 shiftL 13 then if w < 1 shiftL 12 then 12 else 13 else if w < 1 shiftL 14 then 14 else 15 else if w < 1 shiftL 23 then if w < 1 shiftL 19 then if w < 1 shiftL 17 then if w < 1 shiftL 16 then 16 else 17 else if w < 1 shiftL 18 then 18 else 19 else if w < 1 shiftL 21 then if w < 1 shiftL 20 then 20 else 21 else if w < 1 shiftL 22 then 22 else 23 else if w < 1 shiftL 27 then if w < 1 shiftL 25 then if w < 1 shiftL 24 then 24 else 25 else if w < 1 shiftL 26 then 26 else 27 else if w < 1 shiftL 29 then if w < 1 shiftL 28 then 28 else 29 else if w < 1 shiftL 30 then 30 else 31 ---------------------------------------------------------------- {- factorial_primeSwing :: Int -> Integer factorial_primeSwing n0 | n0 < 0 = 0 | n0 < 20 = smallFactorials unsafeAt n0 | otherwise = go n0 shiftL (n0 - popCount n0) where go n | n < 2 = 1 | otherwise = (go (n div 2) ^ 2) * swing n swing n | n < 33 = smallOddSwing unsafeAt n | otherwise = let count = 0 rootN = floorSqrt n xs = primes 3 rootN ys = primes (rootN + 1) (n div 3) in forM_ xs$ \x -> do
let q = n
let p = 1
q := q div x
whileM_ (q > 0) $do when (q .&. 1 == 1) (p := p*x) q := q div x when (p > 1)$ do
primeList !! count := p
count := count+1
forM_ ys $\y -> do when ((n div y) .&. 1 == 1)$ do
primeList !! count := y
count := count+1
return
\$ primorial (n div 2 + 1) n
* xmathProduct primeList 0 count

-- With hsc2hs we can use #def to define these as static C-style arrays, and then use base:Foreign.Marshall.Array to access them. Instead of using array:Data.Array.Unboxed; Or we could try the Addr# trick used in Warp
smallOddSwing :: UArray Int Int32
smallOddSwing = listArray (0,32)
[ 1, 1, 1, 3, 3, 15, 5, 35, 35, 315, 63, 693, 231, 3003
, 429, 6435, 6435, 109395, 12155, 230945, 46189, 969969
, 88179, 2028117, 676039, 16900975, 1300075, 35102025
, 5014575, 145422675, 9694845, 300540195, 300540195 ]

smallFactorials :: UArray Int Int64
smallFactorials = listArray (0,20)
[ 1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800
, 39916800, 479001600, 6227020800, 87178291200, 1307674368000
, 20922789888000, 355687428096000, 6402373705728000
, 121645100408832000, 2432902008176640000 ]

-- Added to Bits class in base-4.5.0.0==ghc-7.4.1
-- cf <http://wiki.cs.pdx.edu/forge/popcount.html>
-- cf <http://en.wikipedia.org/wiki/Hamming_weight>
-- | The number of set bits.
popCount :: Int -> Int
popCount x0 =
let x1 = x0 - w2i ((w1 .&. i2w x0) shiftR 1)
x2 = (x1 .&. m2) + ((x1 shiftR 2) .&. m2)
x3 = (x2 + (x2 shiftR 4)) .&. m4
x4 = x3 + (x3 shiftR 8)
x5 = x4 + (x4 shiftR 16)
x6 = x5 + (x5 shiftR 32) -- for 64-bit platforms
in x6 .&. 0x7f
where
i2w :: Int -> Word
i2w = fromIntegral

w2i :: Word -> Int
w2i = fromIntegral

w1 = 0xaaaaaaaaaaaaaaaa    -- binary: 0101...
-- m1 = 0x5555555555555555 -- binary: 1010...
m2 = 0x3333333333333333    -- binary: 11001100...
m4 = 0x0f0f0f0f0f0f0f0f    -- binary: 11110000...

factorial_parallelPrimeSwing
-}
----------------------------------------------------------------
----------------------------------------------------------- fin.