diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. + diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,27 @@ +#!/usr/bin/env runhaskell +-- Cf. +-- + +{-# 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. diff --git a/combinatorics.cabal b/combinatorics.cabal new file mode 100644 --- /dev/null +++ b/combinatorics.cabal @@ -0,0 +1,51 @@ +---------------------------------------------------------------- +-- wren ng thornton ~ 2012.01.12 +---------------------------------------------------------------- + +Name: combinatorics +Version: 0.1.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 computation of common combinatoric functions. +Description: Efficient computation of common 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.Primes + , Math.Combinatorics.Factorial + , Math.Combinatorics.Binomial + -- Data.IntList + + if flag(base4) + Build-Depends: base >= 4 && < 5 + else + Build-Depends: base < 4 + +---------------------------------------------------------------- +----------------------------------------------------------- fin. diff --git a/src/Math/Combinatorics/Binomial.hs b/src/Math/Combinatorics/Binomial.hs new file mode 100644 --- /dev/null +++ b/src/Math/Combinatorics/Binomial.hs @@ -0,0 +1,139 @@ +{-# OPTIONS_GHC -Wall -fwarn-tabs #-} +---------------------------------------------------------------- +-- 2012.01.28 +-- | +-- Module : Math.Combinatorics.Binomial +-- Copyright : Copyright (c) 2011 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.Binomial (choose) where + +import Data.List (foldl') +import Math.Combinatorics.Primes (primes) + +{- + + +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= 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. +-- , +-- +-- +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. diff --git a/src/Math/Combinatorics/Factorial.hs b/src/Math/Combinatorics/Factorial.hs new file mode 100644 --- /dev/null +++ b/src/Math/Combinatorics/Factorial.hs @@ -0,0 +1,277 @@ +{-# OPTIONS_GHC -Wall -fwarn-tabs #-} +{-# LANGUAGE CPP #-} +---------------------------------------------------------------- +-- 2012.01.28 +-- | +-- Module : Math.Combinatorics.Factorial +-- Copyright : Copyright (c) 2011--2012 wren ng thornton +-- License : BSD +-- Maintainer : wren@community.haskell.org +-- Stability : provisional +-- Portability : Haskell98 + CPP +-- +-- The factorial numbers (). 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.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 + +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 ] + + +-- cf +-- cf +-- | 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. diff --git a/src/Math/Combinatorics/Primes.hs b/src/Math/Combinatorics/Primes.hs new file mode 100644 --- /dev/null +++ b/src/Math/Combinatorics/Primes.hs @@ -0,0 +1,61 @@ +{-# OPTIONS_GHC + -Wall + -fwarn-tabs + -fno-warn-incomplete-patterns + -fno-warn-name-shadowing + #-} +---------------------------------------------------------------- +-- 2011.12.07 +-- | +-- Module : Math.Combinatorics.Primes +-- Copyright : Copyright (c) 2011 wren ng thornton +-- License : BSD +-- Maintainer : wren@community.haskell.org +-- Stability : provisional +-- Portability : Haskell98 +-- +-- The prime numbers (). +---------------------------------------------------------------- +module Math.Combinatorics.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 +-- +-- +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.