exact-combinatorics (empty) → 0.2.0
raw patch · 6 files changed
+592/−0 lines, 6 filesdep +basebuild-type:Customsetup-changed
Dependencies added: base
Files
- LICENSE +33/−0
- Setup.hs +27/−0
- exact-combinatorics.cabal +51/−0
- src/Math/Combinatorics/Exact/Binomial.hs +139/−0
- src/Math/Combinatorics/Exact/Factorial.hs +278/−0
- src/Math/Combinatorics/Exact/Primes.hs +64/−0
+ LICENSE view
@@ -0,0 +1,33 @@+Copyright (c) 2011, 2012, wren ng thornton.+ALL RIGHTS RESERVED.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of the copyright holders nor the names of+ other contributors may be used to endorse or promote products+ derived from this software without specific prior written+ permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS+FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE+COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,+INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,+BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN+ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE+POSSIBILITY OF SUCH DAMAGE.+
+ Setup.hs view
@@ -0,0 +1,27 @@+#!/usr/bin/env runhaskell+-- Cf. <http://www.mail-archive.com/haskell-cafe@haskell.org/msg59984.html>+-- <http://www.haskell.org/pipermail/haskell-cafe/2008-December/051785.html>++{-# OPTIONS_GHC -Wall -fwarn-tabs -fno-warn-missing-signatures #-}+module Main (main) where+import Distribution.Simple+import Distribution.Simple.LocalBuildInfo (withPrograms)+import Distribution.Simple.Program (userSpecifyArgs)+----------------------------------------------------------------++-- | Define __HADDOCK__ when building documentation.+main :: IO ()+main = defaultMainWithHooks+ $ simpleUserHooks `modify_haddockHook` \oldHH pkg lbi hooks flags -> do+ + -- Call the old haddockHook with a modified LocalBuildInfo+ (\lbi' -> oldHH pkg lbi' hooks flags)+ $ lbi `modify_withPrograms` \oldWP ->+ userSpecifyArgs "haddock" ["--optghc=-D__HADDOCK__"] oldWP+++modify_haddockHook hooks f = hooks { haddockHook = f (haddockHook hooks) }+modify_withPrograms lbi f = lbi { withPrograms = f (withPrograms lbi) }++----------------------------------------------------------------+----------------------------------------------------------- fin.
+ exact-combinatorics.cabal view
@@ -0,0 +1,51 @@+----------------------------------------------------------------+-- wren ng thornton <wren@community.haskell.org> ~ 2012.02.02+----------------------------------------------------------------++Name: exact-combinatorics+Version: 0.2.0+Stability: provisional+Homepage: http://code.haskell.org/~wren/+Author: wren ng thornton+Maintainer: wren@community.haskell.org+Copyright: Copyright (c) 2011--2012 wren ng thornton+License: BSD3+License-File: LICENSE++Category: Statistics, Math+Synopsis: Efficient exact computation of combinatoric functions.+Description: Efficient exact computation of combinatoric functions.++-- By and large Cabal >=1.2 is fine; but >= 1.6 gives tested-with:+-- and source-repository:.+Cabal-Version: >= 1.6+-- We need a custom build in order to define __HADDOCK__+Build-Type: Custom+Tested-With: GHC == 6.12.1++Source-Repository head+ Type: darcs+ Location: http://community.haskell.org/~wren/combinatorics++----------------------------------------------------------------+Flag base4+ Default: True+ Description: base-4.0 emits "Prelude deprecated" messages in+ order to get people to be explicit about which+ version of base they use.++----------------------------------------------------------------+Library+ Hs-Source-Dirs: src+ Exposed-Modules: Math.Combinatorics.Exact.Primes+ , Math.Combinatorics.Exact.Factorial+ , Math.Combinatorics.Exact.Binomial+ -- Data.IntList+ + if flag(base4)+ Build-Depends: base >= 4 && < 5+ else+ Build-Depends: base < 4++----------------------------------------------------------------+----------------------------------------------------------- fin.
+ src/Math/Combinatorics/Exact/Binomial.hs view
@@ -0,0 +1,139 @@+{-# OPTIONS_GHC -Wall -fwarn-tabs #-}+----------------------------------------------------------------+-- 2012.02.02+-- |+-- Module : Math.Combinatorics.Exact.Binomial+-- Copyright : Copyright (c) 2011--2012 wren ng thornton+-- License : BSD+-- Maintainer : wren@community.haskell.org+-- Stability : provisional+-- Portability : Haskell98+--+-- Binomial coefficients, aka the count of possible combinations.+-- For negative inputs, all functions return 0 (rather than throwing+-- an exception or using 'Maybe').+----------------------------------------------------------------+module Math.Combinatorics.Exact.Binomial (choose) where++import Data.List (foldl')+import Math.Combinatorics.Exact.Primes (primes)++{-+<http://mathworld.wolfram.com/BinomialCoefficient.html>++Some identities, but not really material for RULES:+ n `choose` 0 = 1+ n `choose` 1 = n+ n `choose` 2 = 2*n*(n-1)+ n `choose` k = n `choose` (n-k) when 0<=k<=n+ n `choose` k = (-1)^k * ((k-n-1) `choose` k)+ n `choose` (k+1) = (n `choose` k) * ((n-k) / (k+1))+ (n+1) `choose` k = (n `choose` k) * (n `choose` (k-1))+ n `choose` j = ((n-1)`choose` j) + ((n-1)`choose`(j-1)) when 0<j<n++Regarding the prime factorization/carries thing, also cf:+ Kummer (1852);+ Graham et al. (1989), Exercise 5.36, p. 245;+ Ribenboim (1989);+ Vardi (1991), p. 68++To extend to negative arguments and to complex numbers, see (Kronenburg 2011):+ n `choose` k+ | k >= 0 = (-1)^k * ((-n+k-1) `choose` k)+ | k <= n = (-1)^(n-k) * ((-k-1) `choose` (n-k))+ | otherwise = 0++According to Grinstead&Snell, p.95, when using the naive implementation+if you alternatete the multiplications and divisions then all+intermediate values are integers and none of the intermediate values+exceeds the final value. This property is retained in the fast+implementation.+-}+++-- TODO: give a version that returns the prime-power factorization as [(Int,Int)]+++-- | Exact binomial coefficients. For a fast approximation see+-- @math-functions:Numeric.SpecFunctions.choose@ instead. The naive+-- definition of the binomial coefficients is:+--+-- > n `choose` k+-- > | k < 0 = 0+-- > | k > n = 0+-- > | otherwise = factorial n `div` (factorial k * factorial (n-k))+--+-- However, we use a fast implementation based on the prime-power+-- factorization of the result (Goetgheluck, 1987). Each time @n@+-- is larger than the previous calls, there will be some slowdown+-- as the prime numbers must be computed (though it is still much+-- faster than the naive implementation); however, subsequent calls+-- will be extremely fast, since we memoize the list of 'primes'.+-- Do note, however, that this will result in a space leak if you+-- call @choose@ for an extremely large @n@ and then don't need+-- that many primes in the future. Hopefully future versions will+-- correct this issue.+--+-- * P. Goetgheluck (1987)+-- /Computing Binomial Coefficients/,+-- American Mathematical Monthly, 94(4). pp.360--365.+-- <http://www.jstor.org/stable/2323099>,+-- <http://dl.acm.org/citation.cfm?id=26272>+--+choose :: (Integral a) => a -> a -> a+ -- The result type could be any (Num b) if desired.+{-# SPECIALIZE choose ::+ Integer -> Integer -> Integer,+ Int -> Int -> Int+ #-}+n `choose` k_+ | n `seq` k_`seq` False = undefined+ | 0 < k_ && k_ < n = + k `seq` nk `seq` sqrtN `seq`+ foldl'+ (\acc prime -> step acc (fromIntegral prime))+ 1+ (takeWhile (fromIntegral n >=) primes)+ -- BUG: 'takeWhile' isn't a good producer, so we shouldn't+ -- just @map fromIntegral@. In newer GHC my patch will make+ -- it in for it to be a good producer (and a good consumer).+ | 0 <= k_ && k_ <= n = 1 -- N.B., @binomial_naive 0 0 == 1@+ | otherwise = 0+ where+ -- TODO: since we know the second operands to quot/rem are+ -- positive, we should use quotInt/remInt directly to avoid the+ -- extra tests (the overflow errors are not optimized away).+ + k = fromIntegral $! if k_ > n `quot` 2 then n - k_ else k_+ nk = n - k+ sqrtN = floor (sqrt (fromIntegral n) :: Double) `asTypeOf` n++ step acc prime+ | acc `seq` prime `seq` False = undefined+ | prime > nk = acc * prime+ | prime > n `quot` 2 = acc+ | prime > sqrtN =+ if n `rem` prime < k `rem` prime+ then acc * prime+ else acc+ | otherwise = acc * go n k 0 1+ where+ go n' k' r p+ | n' `seq` k' `seq` r `seq` p `seq` False = undefined+ | n' <= 0 = p+ | n' `rem` prime < (k' `rem` prime) + r+ = go (n' `quot` prime) (k' `quot` prime) 1 $! p * prime+ | otherwise = go (n' `quot` prime) (k' `quot` prime) 0 p+ + {- -- BENCH: apparently this is an unreliable optimization.+ | otherwise = acc * (prime ^ go n k 0 0)+ where+ go n' k' r p+ | n' <= 0 = p `asTypeOf` acc+ | n' `rem` prime < (k' `rem` prime) + r+ = go (n' `quot` prime) (k' `quot` prime) 1 $! p+1+ | otherwise = go (n' `quot` prime) (k' `quot` prime) 0 p+ -}++----------------------------------------------------------------+----------------------------------------------------------- fin.
+ src/Math/Combinatorics/Exact/Factorial.hs view
@@ -0,0 +1,278 @@+{-# 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.
+ src/Math/Combinatorics/Exact/Primes.hs view
@@ -0,0 +1,64 @@+{-# OPTIONS_GHC+ -Wall+ -fwarn-tabs+ -fno-warn-incomplete-patterns+ -fno-warn-name-shadowing+ #-}+----------------------------------------------------------------+-- 2012.02.02+-- |+-- Module : Math.Combinatorics.Exact.Primes+-- Copyright : Copyright (c) 2011 wren ng thornton+-- License : BSD+-- Maintainer : wren@community.haskell.org+-- Stability : provisional+-- Portability : Haskell98+--+-- The prime numbers (<http://oeis.org/A000040>).+----------------------------------------------------------------+module Math.Combinatorics.Exact.Primes (primes) where+++data Wheel = Wheel {-# UNPACK #-}!Int ![Int]+++-- BUG: the CAF is nice for sharing, but what about when we want+-- fusion and to avoid sharing? Using Data.IntList seems to only+-- increase the overhead. I guess things aren't being memoized/freed+-- like they should...++-- | The prime numbers. Implemented with the algorithm in:+--+-- * Colin Runciman (1997)+-- /Lazy Wheel Sieves and Spirals of Primes/, Functional Pearl,+-- Journal of Functional Programming, 7(2). pp.219--225.+-- ISSN 0956-7968+-- <http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.55.7096>+--+primes :: [Int]+primes = seive wheels primes primeSquares+ where+ primeSquares = [p*p | p <- primes]+ + wheels = Wheel 1 [1] : zipWith nextSize wheels primes+ where+ nextSize (Wheel s ns) p =+ Wheel (s*p) [n' | o <- [0,s..(p-1)*s]+ , n <- ns+ , n' <- [n+o]+ , n' `mod` p > 0 ]+ + -- N.B., ps and qs must be lazy. Or else the circular program is _|_.+ seive (Wheel s ns : ws) ps qs =+ [ n' | o <- s : [2*s,3*s..(head ps-1)*s]+ , n <- ns+ , n' <- [n+o]+ , s <= 2 || noFactorIn ps qs n' ]+ ++ seive ws (tail ps) (tail qs)+ where+ -- noFactorIn :: [Int] -> [Int] -> Int -> Bool+ noFactorIn (p:ps) (q:qs) x =+ q > x || x `mod` p > 0 && noFactorIn ps qs x++----------------------------------------------------------------+----------------------------------------------------------- fin.