diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,26 @@
+Copyright (c) Henning Thielemann 2013
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+1. Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+2. 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.
+3. Neither the name of the University nor the names of its contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE REGENTS 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 REGENTS 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.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,3 @@
+#! /usr/bin/env runhaskell
+> import Distribution.Simple
+> main = defaultMain
diff --git a/example/Parallelism.hs b/example/Parallelism.hs
new file mode 100644
--- /dev/null
+++ b/example/Parallelism.hs
@@ -0,0 +1,24 @@
+module Parallelism where
+
+import qualified Control.Concurrent.MVar as MVar
+import Control.Concurrent (forkIO, getNumCapabilities)
+import Control.Exception (finally)
+
+import Control.Functor.HT (void)
+import Data.Foldable (forM_)
+
+
+schedule :: [IO ()] -> IO ()
+schedule acts = do
+   n <- getNumCapabilities
+   let (start, queue) = splitAt n acts
+   mvar <- MVar.newEmptyMVar
+   let newJob act = void $ forkIO $ finally act $ MVar.putMVar mvar ()
+   mapM_ newJob start
+   let loop [] = return ()
+       loop (act:remain) = do
+          MVar.takeMVar mvar
+          newJob act
+          loop remain
+   loop queue
+   forM_ start $ const $ MVar.takeMVar mvar
diff --git a/example/Queen8.hs b/example/Queen8.hs
new file mode 100644
--- /dev/null
+++ b/example/Queen8.hs
@@ -0,0 +1,52 @@
+module Main where
+
+import qualified Math.SetCover.Exact as ESC
+
+import Control.Monad (liftM2)
+
+import qualified Data.Array as Array
+import qualified Data.Set as Set
+import Data.Array (accumArray)
+import Data.Set (Set)
+import Data.List.HT (sliceVertical)
+import Data.List (intersperse)
+import Data.Maybe (catMaybes)
+
+
+n :: Int
+n = 8
+
+range :: [Int]
+range = [0 .. n-1]
+
+
+data X = Row Int | Column Int | Diag Int | Gaid Int
+         deriving (Eq, Ord, Show)
+
+type Assign = ESC.Assign (Maybe (Int, Int)) (Set X)
+
+assign :: Int -> Int -> Assign
+assign i j =
+   ESC.assign (Just (i,j)) $
+   Set.fromList [Row i, Column j, Diag (i+j), Gaid (i-j)]
+
+fill :: X -> Assign
+fill = ESC.assign Nothing . Set.singleton
+
+assigns :: [Assign]
+assigns =
+   liftM2 assign range range
+   ++
+   map (fill . Diag) [0 .. 2*(n-1)]
+   ++
+   map (fill . Gaid) [1-n .. n-1]
+
+format :: [Maybe (Int,Int)] -> String
+format =
+   unlines . map (intersperse ' ') . sliceVertical n . Array.elems .
+   accumArray (flip const) '.' ((0,0),(n-1,n-1)) .
+   map (flip (,) 'Q') . catMaybes
+
+
+main :: IO ()
+main = mapM_ (putStrLn . format) $ ESC.partitions assigns
diff --git a/example/Soma.hs b/example/Soma.hs
new file mode 100644
--- /dev/null
+++ b/example/Soma.hs
@@ -0,0 +1,128 @@
+module Main where
+
+import qualified Math.SetCover.Exact as ESC
+import qualified Math.SetCover.BitSet as BitSet
+import qualified Math.SetCover.Bit as Bit
+import Math.SetCover.Cuboid
+          (PackedCoords(PackedCoords), Coords, Size, forNestedCoords,
+           allPositions, allOrientations, packCoords, unpackCoords)
+
+import Data.Word (Word8, Word32)
+
+import qualified Data.Map as Map
+import qualified Data.Set as Set
+import qualified Data.List.Match as Match
+
+import Control.Applicative (pure)
+import Data.Foldable (foldMap)
+import Data.Function.HT (nest)
+import Data.List (minimumBy, intercalate)
+
+
+shapes :: [[PackedCoords]]
+shapes =
+   map (map PackedCoords) $
+   [[0,1,2,3], [0,1,3], [0,1,2,4], [0,1,4,5],
+    [0,1,3,9], [0,1,3,10], [0,1,3,12]]
+
+size :: Size
+size = pure 3
+
+
+newtype Brick = Brick Int deriving (Eq, Ord, Show)
+
+type Mask = Set.Set (Either Brick PackedCoords)
+
+type Assign = ESC.Assign (Map.Map PackedCoords Brick) Mask
+
+transformedBrickAssign :: Brick -> [PackedCoords] -> [Assign]
+transformedBrickAssign k =
+    map (brickAssign k) . concatMap (allPositions size) .
+    (if k == Brick 0 then \x->[x] else allOrientations) .
+    map (unpackCoords size)
+
+brickAssign :: Brick -> [Coords Int] -> Assign
+brickAssign k ts =
+   let xs = map (packCoords size) ts
+   in  ESC.assign (Map.fromList $ map (flip (,) k) xs) $
+       Set.fromList $ Left k : map Right xs
+
+allAssigns :: [Assign]
+allAssigns = concat $ zipWith transformedBrickAssign (map Brick [0 ..]) shapes
+
+allMasks :: [Mask]
+allMasks = map ESC.labeledSet allAssigns
+
+writeMasks :: IO ()
+writeMasks =
+   writeFile "somaA.txt" $ show allMasks
+
+format :: [Map.Map PackedCoords Brick] -> String
+format v =
+   let wuerfelx = Map.unions v
+   in  forNestedCoords
+          unlines (intercalate " | ") (intercalate " ")
+          (\c ->
+             maybe "." (\(Brick n) -> show n) $
+             Map.lookup (packCoords size c) wuerfelx)
+          size
+
+printMask :: [Map.Map PackedCoords Brick] -> IO ()
+printMask = putStrLn . format
+
+
+-- Setcovering
+omega :: Mask
+omega = Set.unions allMasks
+
+-- Erweiterungsoptionen für eine Lösung
+e :: [Assign] -> [Mask] -> Either Brick PackedCoords -> [Assign]
+e aktiv x o =
+   filter (ESC.disjoint (Set.unions x) . ESC.labeledSet) $
+   filter (Set.member o . ESC.labeledSet) $
+   aktiv
+
+-- Erweiterung der gesamten (partiellen) Lösungsmenge
+ew :: [Assign] -> [[Assign]]
+ew x =
+   map (:x) $
+   minimumBy Match.compareLength $
+   map (e allAssigns $ map ESC.labeledSet x) $
+   Set.toList $ Set.difference omega $ foldMap ESC.labeledSet x
+
+
+type BitMask = BitSet.Set (Bit.Sum Word8 Word32)
+
+packMask :: Mask -> BitMask
+packMask =
+   foldMap
+      (BitSet.Set .
+       either
+          (\(Brick x) -> Bit.bitLeft x)
+          (\(PackedCoords x) -> Bit.bitRight x))
+
+
+testme :: Brick -> IO ()
+testme b@(Brick n) =
+   mapM_ (printMask . (:[]) . ESC.label) $ transformedBrickAssign b (shapes!!n)
+
+main, mainBase, mainState, mainBits, testme0, testme1 :: IO ()
+testme0 = testme $ Brick 0
+testme1 = testme $ Brick 1
+
+mainBase = do
+   let lsg = map (map ESC.label) $ nest (length shapes) (concatMap ew) [[]]
+   mapM_ printMask lsg
+   print $ length lsg
+
+mainState = do
+   let lsg = ESC.partitions allAssigns
+   mapM_ printMask lsg
+   print $ length lsg
+
+mainBits = do
+   let lsg = ESC.partitions $ map (fmap packMask) allAssigns
+   mapM_ printMask lsg
+   print $ length lsg
+
+main = mainBits
diff --git a/example/Sudoku.hs b/example/Sudoku.hs
new file mode 100644
--- /dev/null
+++ b/example/Sudoku.hs
@@ -0,0 +1,104 @@
+module Main where
+
+import qualified Math.SetCover.BitSet as BitSet
+import qualified Math.SetCover.Bit as Bit
+import qualified Math.SetCover.Exact as ESC
+
+import Data.Word (Word32, Word64)
+
+import Control.Monad (liftM3, guard)
+
+import qualified Data.Array as Array
+import qualified Data.Set as Set
+import Data.Array (array)
+import Data.Set (Set)
+import Data.List.HT (sliceVertical)
+import Data.List (intersperse)
+
+
+data X = Pos Int Int | Row Int Int | Column Int Int | Square Int Int Int
+         deriving (Eq, Ord, Show)
+
+type Assign = ESC.Assign ((Int, Int), Int)
+
+assign :: Int -> Int -> Int -> Assign (Set X)
+assign k i j =
+   ESC.assign ((i,j), k) $
+   Set.fromList [Pos i j, Row k i, Column k j, Square k (div i 3) (div j 3)]
+
+assigns :: [Assign (Set X)]
+assigns = liftM3 assign [1..9] [0..8] [0..8]
+
+
+type Word81 = Bit.Sum Word64 Word32
+type Mask =
+   BitSet.Set
+      (Bit.Sum
+         (Bit.Sum Word81 Word81)
+         (Bit.Sum Word81 Word81))
+
+bit9x9 :: Int -> Int -> Word81
+bit9x9 i j =
+   let k = i*9+j
+   in  if k<64
+         then Bit.bitLeft k
+         else Bit.bitRight (k-64)
+
+bitAssign :: Int -> Int -> Int -> Assign Mask
+bitAssign k i j =
+   ESC.assign ((i,j), k) $
+   BitSet.Set $
+   Bit.Sum
+      (Bit.Sum (bit9x9 k i) (bit9x9 k j))
+      (Bit.Sum (bit9x9 i j) (bit9x9 k (div i 3 + 3 * div j 3)))
+
+bitAssigns :: [Assign Mask]
+bitAssigns = liftM3 bitAssign [1..9] [0..8] [0..8]
+
+
+format :: [((Int, Int), Int)] -> String
+format =
+   unlines . map (intersperse ' ') . sliceVertical 9 . Array.elems .
+   fmap (\n -> toEnum $ n + fromEnum '0') .
+   array ((0,0),(8,8))
+
+
+exampleHawiki1 :: [String]
+exampleHawiki1 =
+   "    6  8 " :
+   " 2       " :
+   "  1      " :
+   " 7    1 2" :
+   "5   3    " :
+   "      4  " :
+   "  42 1   " :
+   "3  7  6  " :
+   "       5 " :
+   []
+
+stateFromString ::
+   (ESC.Set set) =>
+   [Assign set] ->
+   (Int -> Int -> Int -> Assign set) ->
+   [String] -> ESC.State ((Int, Int), Int) set
+stateFromString asgns asgn css =
+   foldl (flip ESC.updateState) (ESC.initState asgns) $
+   do (i,cs) <- zip [0..] css
+      (j,c)  <- zip [0..] cs
+      guard $ c/=' '
+      return $ asgn (fromEnum c - fromEnum '0') i j
+
+
+main, mainAll, mainSolve, mainBit :: IO ()
+mainAll =
+   mapM_ (putStrLn . format) $ ESC.partitions bitAssigns
+
+mainSolve =
+   mapM_ (putStrLn . format) $ ESC.search $
+   stateFromString assigns assign exampleHawiki1
+
+mainBit =
+   mapM_ (putStrLn . format) $ ESC.search $
+   stateFromString bitAssigns bitAssign exampleHawiki1
+
+main = mainBit
diff --git a/example/TetrisCube.hs b/example/TetrisCube.hs
new file mode 100644
--- /dev/null
+++ b/example/TetrisCube.hs
@@ -0,0 +1,273 @@
+{-
+One solution:
+0 0 0 1 | 0 7 1 1 | 0 4 1 5 | 4 4 1 2
+7 B 0 8 | 7 7 5 5 | 7 3 6 5 | 4 3 2 2
+B B B 8 | 7 A 5 2 | 3 3 6 2 | 4 6 6 2
+A A B 8 | A A B 8 | 3 9 6 8 | 9 9 9 9
+
+Another one with colors:
+0 0 0 1 | 0 0 0 0 | 0 0 2 2 | 1 0 1 1
+0 2 0 1 | 0 0 3 2 | 3 1 1 2 | 1 1 1 2
+0 2 3 1 | 3 3 3 1 | 3 3 3 2 | 1 3 3 2
+0 2 2 2 | 3 3 2 1 | 3 3 3 2 | 1 3 2 2
+
+
+Algorithm by Helmut Podhaisky:
+It is a depth-first search where in each stage we choose a position
+where as few as possible bricks match. (see 'ew')
+The function 'ESC.step' is a slightly more efficient version
+that permanently manages the set of available bricks.
+
+dist/build/tetris-cube/tetris-cube +RTS -N4 -M500m
+-}
+module Main where
+
+import qualified Math.SetCover.Exact as ESC
+import qualified Math.SetCover.BitSet as BitSet
+import qualified Math.SetCover.Bit as Bit
+import Math.SetCover.Cuboid
+          (PackedCoords(PackedCoords), Coords(Coords), Size, forNestedCoords,
+           allPositions, allOrientations, packCoords)
+
+import Parallelism (schedule)
+
+import qualified Data.Map as Map
+import qualified Data.Set as Set
+import qualified Data.List.Match as Match
+import qualified Data.Foldable as Fold
+
+import Control.Applicative (pure)
+import Data.Function (on)
+import Data.Foldable (foldMap)
+import Data.List (intercalate, sortBy)
+import Data.Word (Word16, Word64)
+
+import qualified System.IO as IO
+import Text.Printf (printf)
+
+
+shapes, blueShapes, yellowShapes, redShapes :: [[String]]
+shapes = blueShapes ++ yellowShapes ++ redShapes
+
+blueShapes =
+   (
+   "..." :
+   ".  " :
+   ":  " :
+   [])
+   :
+   (
+   "... " :
+   "  .." :
+   [])
+   :
+   (
+   "..." :
+   ":  " :
+   ".  " :
+   [])
+   :
+   (
+   "':." :
+   "  ." :
+   [])
+   :
+   []
+
+
+yellowShapes =
+   (
+   "..." :
+   ":  " :
+   [])
+   :
+   (
+   ".. " :
+   " :'" :
+   [])
+   :
+   (
+   "..." :
+   " : " :
+   [])
+   :
+   (
+   " . " :
+   ".:." :
+   " ' " :
+   [])
+   :
+   []
+
+
+redShapes =
+   (
+   "..." :
+   ".  " :
+   ".  " :
+   [])
+   :
+   (
+   "...." :
+   " .  " :
+   [])
+   :
+   (
+   ":." :
+   ".." :
+   [])
+   :
+   (
+   ":  " :
+   "..." :
+   " . " :
+   [])
+   :
+   []
+
+
+numberOfAtoms :: Int
+numberOfAtoms =
+   Fold.sum $
+   Map.intersectionWith (*)
+      (Map.fromList [('.', 1), ('\'', 1), (':', 2)]) $
+   Map.fromListWith (+) $
+   map (flip (,) 1) $
+   concat $ concat shapes
+
+propNumberOfAtoms :: Bool
+propNumberOfAtoms = numberOfAtoms == 64
+
+
+coordsFromString :: [String] -> [Coords Int]
+coordsFromString ss = do
+   (rowN,row) <- zip [0..] ss
+   (colN,c) <- zip [0..] row
+   fmap (Coords rowN colN) $
+      case c of
+         ' '  -> []
+         '.'  -> [0]
+         '\'' -> [1]
+         ':'  -> [0, 1]
+         _ -> error "forbidden character"
+
+
+size :: Size
+size = pure 4
+
+
+data Color = Blue | Yellow | Red deriving (Eq, Ord, Enum, Show)
+
+type BrickId = (Color, Int)
+
+type Mask = Set.Set (Either BrickId PackedCoords)
+
+type Assign = ESC.Assign (Map.Map PackedCoords BrickId) Mask
+
+transformedBrickAssign :: BrickId -> [String] -> [Assign]
+transformedBrickAssign k =
+   map (brickAssign k) . concatMap (allPositions size) .
+   (if k==(Blue,0) then (:[]) else allOrientations) . coordsFromString
+
+brickAssign :: BrickId -> [Coords Int] -> Assign
+brickAssign k ts =
+   let xs = map (packCoords size) ts
+   in  ESC.assign (Map.fromList $ map (flip (,) k) xs) $
+       Set.fromList $ Left k : map Right xs
+
+allAssigns :: [Assign]
+allAssigns =
+   let gen color =
+          concat . zipWith transformedBrickAssign (map ((,) color) [0..])
+   in  gen Blue blueShapes ++
+       gen Yellow yellowShapes ++
+       gen Red redShapes
+
+allMasks :: [Mask]
+allMasks = map ESC.labeledSet allAssigns
+
+writeMasks :: IO ()
+writeMasks =
+   writeFile "tetriscube.txt" $ show allMasks
+
+formatBrickId :: BrickId -> String
+formatBrickId (color, num) =
+   case color of
+      Red -> "\ESC[31m"
+      Yellow -> "\ESC[33m"
+      Blue -> "\ESC[34m"
+   ++
+   show num
+   ++
+   "\ESC[m"
+
+
+hPutStrLnImmediate :: IO.Handle -> String -> IO ()
+hPutStrLnImmediate h str = do
+   IO.hPutStrLn h str
+   IO.hFlush h
+
+format :: [Map.Map PackedCoords BrickId] -> String
+format v =
+   let wuerfelx = Map.unions v
+   in  forNestedCoords
+          unlines (intercalate " | ") (intercalate " ")
+          (\c ->
+             maybe "." formatBrickId $
+             Map.lookup (packCoords size c) wuerfelx)
+          size
+
+printMask :: [Map.Map PackedCoords BrickId] -> IO ()
+printMask =
+   hPutStrLnImmediate IO.stdout . format
+
+
+type BitMask = BitSet.Set (Bit.Sum Word16 Word64)
+
+
+packMask :: Mask -> BitMask
+packMask =
+   foldMap
+      (BitSet.Set .
+       either
+          (\(color, n) -> Bit.bitLeft $ fromEnum color * 4 + n)
+          (\(PackedCoords x) -> Bit.bitRight x))
+
+
+testme :: BrickId -> IO ()
+testme b@(color, num) =
+   mapM_ (printMask . (:[]) . ESC.label) $
+   transformedBrickAssign b $ (!!num) $
+   case color of
+      Red -> redShapes
+      Blue -> blueShapes
+      Yellow -> yellowShapes
+
+
+main, mainState, mainBits, mainParallel, testme0, testme1 :: IO ()
+testme0 = testme (Blue, 0)
+testme1 = testme (Blue, 1)
+
+mainState = do
+   let lsg = ESC.partitions allAssigns
+   mapM_ printMask lsg
+   print $ length lsg
+
+mainBits = do
+   let lsg = ESC.partitions $ map (fmap packMask) allAssigns
+   mapM_ printMask lsg
+   print $ length lsg
+
+mainParallel =
+   schedule $ map snd $
+   sortBy (flip Match.compareLength `on` fst) $
+   let attempts =
+          ESC.step $ ESC.initState $ map (fmap packMask) allAssigns
+   in  (\f -> zipWith f [0..] attempts) $ \n attempt ->
+          let refinedAttempts = concatMap ESC.step $ ESC.step attempt
+          in  (refinedAttempts,
+               IO.withFile (printf "tetriscube%02d.txt" (n::Int)) IO.WriteMode $ \h ->
+                  mapM_ (hPutStrLnImmediate h . format) $
+                  concatMap ESC.search refinedAttempts)
+
+main = mainParallel
diff --git a/set-cover.cabal b/set-cover.cabal
new file mode 100644
--- /dev/null
+++ b/set-cover.cabal
@@ -0,0 +1,105 @@
+Name:             set-cover
+Version:          0.0
+License:          BSD3
+License-File:     LICENSE
+Author:           Henning Thielemann, Helmut Podhaisky
+Maintainer:       Henning Thielemann <haskell@henning-thielemann.de>
+Homepage:         http://code.haskell.org/~thielema/set-cover/
+Category:         Math, Algorithms
+Synopsis:         Solve exact set cover problems like Sudoku, 8 Queens, Soma Cube, Tetris Cube
+Description:
+  Solver for exact set cover problems.
+  Included examples: Sudoku, 8 Queens, Soma Cube, Tetris Cube.
+  Generic algorithm allows to choose between
+  slow but flexible @Set@ from @containers@ package
+  and fast but cumbersome bitvectors.
+  .
+  Build examples with @cabal install -fbuildExamples@.
+  .
+  The package needs only Haskell 98.
+Tested-With:      GHC==7.4.2
+Cabal-Version:    >=1.8
+Build-Type:       Simple
+
+Flag buildExamples
+  description: Build example executables
+  default:     False
+
+Source-Repository this
+  Tag:         0.0
+  Type:        darcs
+  Location:    http://code.haskell.org/~thielema/set-cover/
+
+Source-Repository head
+  Type:        darcs
+  Location:    http://code.haskell.org/~thielema/set-cover/
+
+Library
+  Build-Depends:
+    containers >=0.4 && <0.6,
+    utility-ht >=0.0.1 && <0.1,
+    base >=4 && <5
+
+  GHC-Options:      -Wall
+  Hs-Source-Dirs:   src
+  Exposed-Modules:
+    Math.SetCover.Bit
+    Math.SetCover.BitMap
+    Math.SetCover.BitSet
+    Math.SetCover.Exact
+    Math.SetCover.Cuboid
+
+Executable tetris-cube
+  If flag(buildExamples)
+    Build-Depends:
+      set-cover,
+      containers,
+      utility-ht,
+      base
+  Else
+    Buildable: False
+  GHC-Options:    -Wall -rtsopts -threaded
+  Hs-Source-Dirs: example
+  Main-Is: TetrisCube.hs
+  Other-Modules: Parallelism
+
+Executable soma-cube
+  If flag(buildExamples)
+    Build-Depends:
+      set-cover,
+      containers,
+      utility-ht,
+      base
+  Else
+    Buildable: False
+  GHC-Options:    -Wall -rtsopts -threaded
+  Hs-Source-Dirs: example
+  Main-Is: Soma.hs
+
+Executable queen8
+  If flag(buildExamples)
+    Build-Depends:
+      set-cover,
+      containers,
+      array >=0.1 && <0.5,
+      utility-ht,
+      base
+  Else
+    Buildable: False
+  GHC-Options:    -Wall -rtsopts -threaded
+  Hs-Source-Dirs: example
+  Main-Is: Queen8.hs
+
+Executable sudoku-setcover
+  If flag(buildExamples)
+    Build-Depends:
+      set-cover,
+      containers,
+      array >=0.1 && <0.5,
+      utility-ht,
+      base
+  Else
+    Buildable: False
+  GHC-Options:    -Wall -rtsopts -threaded
+  Hs-Source-Dirs: example
+  Main-Is: Sudoku.hs
diff --git a/src/Math/SetCover/Bit.hs b/src/Math/SetCover/Bit.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/SetCover/Bit.hs
@@ -0,0 +1,71 @@
+module Math.SetCover.Bit where
+
+import qualified Data.Bits as Bits
+import Data.Bits (Bits)
+import Data.Word (Word8, Word16, Word32, Word64)
+import Prelude hiding (null)
+
+
+infixl 7 .&.
+infixl 5 .|.
+
+class Eq bits => C bits where
+   empty :: bits
+   complement, keepMinimum :: bits -> bits
+   xor, (.&.), (.|.) :: bits -> bits -> bits
+
+instance C Word8 where
+   empty = 0
+   complement = Bits.complement
+   keepMinimum xs = xs .&. (-xs)
+   xor = Bits.xor
+   (.&.) = (Bits..&.)
+   (.|.) = (Bits..|.)
+
+instance C Word16 where
+   empty = 0
+   complement = Bits.complement
+   keepMinimum xs = xs .&. (-xs)
+   xor = Bits.xor
+   (.&.) = (Bits..&.)
+   (.|.) = (Bits..|.)
+
+instance C Word32 where
+   empty = 0
+   complement = Bits.complement
+   keepMinimum xs = xs .&. (-xs)
+   xor = Bits.xor
+   (.&.) = (Bits..&.)
+   (.|.) = (Bits..|.)
+
+instance C Word64 where
+   empty = 0
+   complement = Bits.complement
+   keepMinimum xs = xs .&. (-xs)
+   xor = Bits.xor
+   (.&.) = (Bits..&.)
+   (.|.) = (Bits..|.)
+
+
+{-
+cf. package largeword
+-}
+data Sum a b = Sum !a !b
+   deriving (Eq, Show)
+
+instance (C a, C b) => C (Sum a b) where
+   empty = Sum empty empty
+   complement (Sum l h) = Sum (complement l) (complement h)
+   xor (Sum xl xh) (Sum yl yh) = Sum (xor xl yl) (xor xh yh)
+   Sum xl xh .&. Sum yl yh = Sum (xl.&.yl) (xh.&.yh)
+   Sum xl xh .|. Sum yl yh = Sum (xl.|.yl) (xh.|.yh)
+   keepMinimum (Sum l h) =
+      if l == empty
+        then Sum empty (keepMinimum h)
+        else Sum (keepMinimum l) empty
+
+bitLeft :: (Bits a, C b) => Int -> Sum a b
+bitLeft n = Sum (Bits.bit n) empty
+
+bitRight :: (C a, Bits b) => Int -> Sum a b
+bitRight n = Sum empty (Bits.bit n)
diff --git a/src/Math/SetCover/BitMap.hs b/src/Math/SetCover/BitMap.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/SetCover/BitMap.hs
@@ -0,0 +1,53 @@
+module Math.SetCover.BitMap where
+
+import qualified Math.SetCover.BitSet as BitSet
+import qualified Math.SetCover.Bit as Bit
+import Math.SetCover.BitSet (Set(Set))
+import Math.SetCover.Bit (xor, (.|.), (.&.))
+
+import Data.Monoid (Monoid, mempty, mappend)
+
+
+{-
+Sliced representation of Map [0..bitSize-1] Integer.
+-}
+newtype Map bits = Map {unMap :: [bits]} deriving (Show)
+
+instance (Bit.C bits) => Monoid (Map bits) where
+   mempty = Map []
+   mappend = add
+
+
+fromSet :: Bit.C bits => Set bits -> Map bits
+fromSet (Set x) = Map [x]
+
+add :: Bit.C bits => Map bits -> Map bits -> Map bits
+add (Map xs0) (Map ys0) =
+   let go c xs [] = unMap $ inc (Set c) (Map xs)
+       go c [] ys = unMap $ inc (Set c) (Map ys)
+       go c (x:xs) (y:ys) =
+          xor c (xor x y) :
+          go (c.&.x .|. x.&.y .|. y.&.c) xs ys
+   in  Map $ go Bit.empty xs0 ys0
+
+inc :: Bit.C bits => Set bits -> Map bits -> Map bits
+inc (Set xs0) (Map ys0) =
+   let go c [] = if c==Bit.empty then [] else [c]
+       go c (x:xs) = xor c x : go (c .&. x) xs
+   in  Map $ go xs0 ys0
+
+
+{-
+Only elements from the base set are considered.
+This way we can distinguish between non-members and members with count zero.
+-}
+minimumSet :: Bit.C bits => Set bits -> Map bits -> Set bits
+minimumSet baseSet (Map xs) =
+   foldr
+      (\x mins ->
+         case BitSet.difference mins $ Set x of
+            newMins ->
+               if BitSet.null newMins
+                 then mins
+                 else newMins)
+      baseSet xs
diff --git a/src/Math/SetCover/BitSet.hs b/src/Math/SetCover/BitSet.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/SetCover/BitSet.hs
@@ -0,0 +1,28 @@
+module Math.SetCover.BitSet where
+
+import qualified Math.SetCover.Bit as Bit
+import Math.SetCover.Bit ((.|.), (.&.))
+
+import Data.Monoid (Monoid, mempty, mappend)
+
+
+newtype Set bits = Set bits deriving (Show)
+
+instance (Bit.C bits) => Monoid (Set bits) where
+   mempty = empty
+   mappend (Set x) (Set y) = Set $ x.|.y
+
+empty :: Bit.C bits => Set bits
+empty = Set Bit.empty
+
+null :: Bit.C bits => Set bits -> Bool
+null (Set xs)  =  xs == Bit.empty
+
+keepMinimum :: Bit.C bits => Set bits -> Set bits
+keepMinimum (Set xs) = Set $ Bit.keepMinimum xs
+
+disjoint :: Bit.C bits => Set bits -> Set bits -> Bool
+disjoint (Set xs) (Set ys)  =  xs.&.ys == Bit.empty
+
+difference :: Bit.C bits => Set bits -> Set bits -> Set bits
+difference (Set xs) (Set ys) = Set $ xs .&. Bit.complement ys
diff --git a/src/Math/SetCover/Cuboid.hs b/src/Math/SetCover/Cuboid.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/SetCover/Cuboid.hs
@@ -0,0 +1,86 @@
+module Math.SetCover.Cuboid where
+
+import qualified Data.Set as Set
+
+import qualified Data.Traversable as Trav
+import qualified Data.Foldable as Fold
+import Control.Applicative (Applicative, liftA2, liftA3, pure, (<*>))
+import Data.List (sort)
+
+
+data Coords a = Coords a a a
+   deriving (Eq, Ord, Show)
+
+instance Functor Coords where
+   fmap f (Coords x y z) = Coords (f x) (f y) (f z)
+
+instance Applicative Coords where
+   pure x = Coords x x x
+   Coords fx fy fz <*> Coords x y z = Coords (fx x) (fy y) (fz z)
+
+instance Fold.Foldable Coords where
+   foldMap = Trav.foldMapDefault
+
+instance Trav.Traversable Coords where
+   traverse f (Coords x y z) = liftA3 Coords (f x) (f y) (f z)
+
+
+forNestedCoords ::
+   (Enum a, Num a) =>
+   ([z] -> b) ->
+   ([y] -> z) ->
+   ([x] -> y) ->
+   (Coords a -> x) ->
+   Coords a -> b
+forNestedCoords fz fy fx f size =
+   case fmap (\k -> [0 .. k-1]) size of
+      Coords rx ry rz ->
+         fz $ flip map rz $ \z ->
+         fy $ flip map ry $ \y ->
+         fx $ flip map rx $ \x ->
+         f (Coords x y z)
+
+
+newtype PackedCoords = PackedCoords Int
+   deriving (Eq, Ord, Show)
+
+
+dx, dy, dz :: Num a => Coords a -> Coords a
+dx (Coords x y z) = Coords x (-z) y -- [1 0  0; 0 0 -1; 0 1 0]
+dy (Coords x y z) = Coords (-z) y x -- [0 0 -1; 0 1 0; 1 0 0]
+dz (Coords x y z) = Coords (-y) x z -- [0 -1 0; 1 0 0; 0 0 1]
+
+rotations :: Num a => [Coords a -> Coords a]
+rotations =
+   liftA2 (.)
+      [id, dx, dx.dx, dx.dx.dx]
+      [id, dz, dz.dz, dz.dz.dz, dy, dy.dy.dy]
+
+
+type Size = Coords Int
+
+unpackCoords :: Size -> PackedCoords -> Coords Int
+unpackCoords size (PackedCoords n) =
+   snd $ Trav.mapAccumL divMod n size
+
+packCoords :: Size -> Coords Int -> PackedCoords
+packCoords size =
+   PackedCoords . Fold.foldr (\(k,x) s -> k*s+x) 0 . liftA2 (,) size
+
+normalForm :: (Ord a, Num a) => [Coords a] -> [Coords a]
+normalForm ts = sort $ map (liftA2 subtract xyzm) ts
+   where xyzm = foldl1 (liftA2 min) ts
+
+allPositions :: Size -> [Coords Int] -> [[Coords Int]]
+allPositions size ts =
+   map (\displacement -> map (liftA2 (+) displacement) ts) $
+   Trav.sequence $
+   liftA2
+      (\k r -> [0 .. k-1-r])
+      size
+      (foldl1 (liftA2 max) ts)
+
+allOrientations :: (Num a, Ord a) => [Coords a] -> [[Coords a]]
+allOrientations ts =
+   Set.toList $ Set.fromList $
+   map (normalForm . flip map ts) rotations
diff --git a/src/Math/SetCover/Exact.hs b/src/Math/SetCover/Exact.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/SetCover/Exact.hs
@@ -0,0 +1,107 @@
+module Math.SetCover.Exact where
+
+import qualified Math.SetCover.BitMap as BitMap
+import qualified Math.SetCover.BitSet as BitSet
+import qualified Math.SetCover.Bit as Bit
+
+import qualified Data.Map as Map
+import qualified Data.Set as Set
+import qualified Data.List as List
+import qualified Data.List.Match as Match
+import qualified Data.Foldable as Fold
+
+import Prelude hiding (null)
+
+
+class Set set where
+   null :: set -> Bool
+   disjoint :: set -> set -> Bool
+   unions :: [set] -> set
+   difference :: set -> set -> set
+   minimize :: set -> [Assign label set] -> [Assign label set]
+
+instance (Ord a) => Set (Set.Set a) where
+   null = Set.null
+   disjoint x y = Set.null $ Set.intersection x y
+   unions = Set.unions
+   difference = Set.difference
+   minimize free =
+      Fold.minimumBy Match.compareLength . Map.unionsWith (++) .
+      (Fold.foldMap (flip Map.singleton []) free :) .
+      map (\a -> Fold.foldMap (flip Map.singleton [a]) $ labeledSet a)
+
+instance (Bit.C a) => Set (BitSet.Set a) where
+   null = BitSet.null
+   disjoint = BitSet.disjoint
+   unions = Fold.fold
+   difference = BitSet.difference
+   minimize free available =
+      let singleMin =
+             BitSet.keepMinimum $ BitMap.minimumSet free $
+             Fold.foldMap (BitMap.fromSet . labeledSet) available
+      in  filter (not . BitSet.disjoint singleMin . labeledSet) available
+
+
+data Assign label set =
+   Assign {
+      label :: label,
+      labeledSet :: set
+   }
+
+assign :: label -> set -> Assign label set
+assign = Assign
+
+
+data State label set =
+   State {
+      availableSubsets :: [Assign label set],
+      freeElements :: set,
+      usedSubsets :: [Assign label set]
+   }
+
+instance Functor (Assign label) where
+   fmap f (Assign lab set) = Assign lab (f set)
+
+instance Functor (State label) where
+   fmap f (State ab fp pb) =
+      State (map (fmap f) ab) (f fp) (map (fmap f) pb)
+
+initState :: Set set => [Assign label set] -> State label set
+initState subsets =
+   State {
+      availableSubsets = subsets,
+      freeElements = unions $ map labeledSet subsets,
+      usedSubsets = []
+   }
+
+{-# INLINE updateState #-}
+updateState :: Set set => Assign label set -> State label set -> State label set
+updateState attempt@(Assign _ attemptedSet) s =
+   State {
+      availableSubsets =
+         filter (disjoint attemptedSet . labeledSet) $
+         availableSubsets s,
+      freeElements = difference (freeElements s) attemptedSet,
+      usedSubsets = attempt : usedSubsets s
+   }
+
+
+{-# INLINE step #-}
+step :: Set set => State label set -> [State label set]
+step s =
+   if List.null (availableSubsets s) || null (freeElements s)
+     then []
+     else
+        map (flip updateState s) $
+        minimize (freeElements s) (availableSubsets s)
+
+{-# INLINE search #-}
+search :: Set set => State label set -> [[label]]
+search s =
+   if null (freeElements s)
+     then [map label $ usedSubsets s]
+     else step s >>= search
+
+{-# INLINE partitions #-}
+partitions :: Set set => [Assign label set] -> [[label]]
+partitions = search . initState