exact-combinatorics-0.2.0: src/Math/Combinatorics/Exact/Factorial.hs
{-# OPTIONS_GHC -Wall -fwarn-tabs #-}
{-# LANGUAGE CPP #-}
----------------------------------------------------------------
-- 2012.02.02
-- |
-- Module : Math.Combinatorics.Exact.Factorial
-- Copyright : Copyright (c) 2011--2012 wren ng thornton
-- License : BSD
-- Maintainer : wren@community.haskell.org
-- Stability : provisional
-- 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.
#ifdef __HADDOCK__
import Data.Int (Int32, Int64)
#endif
import Data.Bits
{-
-- from <http://www.polyomino.f2s.com/david/haskell/hs/CombinatoricsCounting.hs.txt>
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
-- @math-functions:Numeric.SpecFunctions.factorial@ instead. The
-- 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 (`shiftR`1) or (`quot`2) ?
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.