Some combinatorics

.gitignore

@@ -0,0 +1 @@
+dist-newstyle

README.md

@@ -0,0 +1,12 @@
+Counting
+========
+
+Some combinatorics in Haskell. Currently has Bell numbers, binomial
+coefficients, and triangular numbers. And some number of orbits of
+sets with atoms (nominal sets).
+
+Use `cabal` to build:
+
+```
+cabal run counting
+```

Setup.hs

@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain

app/Main.hs

@@ -0,0 +1,27 @@
+module Main where
+
+import Bell
+import Matchings
+import Stirling
+import Control.Monad
+
+heightOfDistinctAtoms :: Int -> Integer
+heightOfDistinctAtoms n = toInteger $ n + 1
+
+heightOfAtoms :: Int -> Integer
+heightOfAtoms n = sum [ norbits k * heightOfDistinctAtoms k | k <- [0 .. n] ]
+                  + sum [ matchings o1 o2 | (o1, o2) <- cross ]
+  where
+    norbits k = secondStirling n k
+    orbits k = replicate (fromIntegral (norbits k)) k
+    allOrbits = concat [ orbits k | k <- [0..n] ]
+    cross = [ (allOrbits !! i, allOrbits !! j) | let m = length allOrbits, i <- [0 .. m-1], j <- [i+1 .. m-1] ]
+
+orbitsOf :: Int -> [Int]
+orbitsOf n = concat [ replicate (fromIntegral (secondStirling n k)) k | k <- [0..n] ]
+
+main :: IO ()
+main = do
+  let out = [ (i, heightOfAtoms i, bell (2*i)) | i <- [0..10] ]
+  forM_ out $ \(i, n, b) -> putStrLn ("h(A^" ++ show i ++ ") = " ++ show n ++ "   <   B_" ++ show (2*i) ++ " = " ++ show b)
+

counting.cabal

@@ -0,0 +1,28 @@
+cabal-version:       2.2
+name:                counting
+version:             0.1.0.0
+author:              Joshua Moerman
+maintainer:          joshua@cs.rwth-aachen.de
+build-type:          Simple
+
+common stuff
+  default-language: Haskell2010
+  build-depends:
+    base >=4.14 && <4.15,
+    containers,
+    MemoTrie
+
+library
+  import: stuff
+  hs-source-dirs: src
+  exposed-modules:
+    Bell,
+    Binomial,
+    Matchings,
+    Stirling
+
+executable counting
+  import: stuff
+  hs-source-dirs: app
+  main-is: Main.hs
+  build-depends: counting

src/Bell.hs

@@ -0,0 +1,23 @@
+module Bell where
+
+import Data.MemoTrie
+import Stirling
+
+--  A000110
+-- number of partitions of a set
+-- also the number of n-types in the structure (N, =)
+bell :: Int -> Integer
+bell = memo bell0 where
+  bell0 n = if n <= 12
+            then table !! n
+            else sum [ secondStirling n k | k <- [0 .. n] ]
+  table = [1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975, 678570, 4213597]
+
+-- number of n-types in the structure (N, <=)
+orderedBell :: Int -> Integer
+orderedBell = memo ob where
+  ob n = if n <= 10
+         then table !! n
+         else sum [ fac k * secondStirling n k | k <- [0 .. n] ]
+  table = [1, 1, 3, 13, 75, 541, 4683, 47293, 545835, 7087261, 102247563]
+  fac k = product [1 .. fromIntegral k]

src/Binomial.hs

@@ -0,0 +1,18 @@
+module Binomial where
+
+import Data.MemoTrie
+
+-- binomial coefficients
+binomial :: Int -> Int -> Integer
+binomial = memo2 b where
+    b n k = if k > n
+            then 0
+            else if n <= 5
+                 then table !! n !! k
+                 else binomial (n-1) k + binomial (n-1) (k-1)
+    table = [ [1]
+            , [1, 1]
+            , [1, 2, 1]
+            , [1, 3, 3, 1]
+            , [1, 4, 6, 4, 1]
+            , [1, 5, 10, 10, 5, 1] ]

src/Matchings.hs

@@ -0,0 +1,12 @@
+module Matchings where
+
+import Data.MemoTrie
+import Binomial
+
+-- number of ways to match up k elements with n elements
+-- matchings may be partial
+matchings :: Int -> Int -> Integer
+matchings = memo2 m where
+  m n k = if n > k
+          then matchings k n
+          else sum [ binomial n i * product [fromIntegral (k-i+1) .. fromIntegral k] | i <- [0 .. n]]

src/Stirling.hs

@@ -0,0 +1,20 @@
+module Stirling where
+
+import Data.MemoTrie
+
+-- Number of ways to partition n into k subsets
+secondStirling :: Int -> Int -> Integer
+secondStirling = memo2 ss where
+    ss n k = if k > n
+             then 0
+             else if n <= 6
+                  then table !! n !! k
+                  else (fromIntegral k) * secondStirling (n - 1) k + secondStirling n (k - 1)
+    table =
+        [ [1]
+        , [0, 1]
+        , [0, 1, 1]
+        , [0, 1, 3, 1]
+        , [0, 1, 7, 6, 1]
+        , [0, 1, 15, 25, 10, 1]
+        , [0, 1, 31, 90, 65, 15, 1] ]