diff --git a/LICENSE.txt b/LICENSE.txt
new file mode 100644
--- /dev/null
+++ b/LICENSE.txt
@@ -0,0 +1,28 @@
+Copyright (c) 2017, Arthur Lee
+
+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 Arthur Lee 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/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,10 @@
+# The `exact-cover` package
+
+`exact-cover` is a Haskell library that implements a fast solver
+for [exact cover problems](https://en.wikipedia.org/wiki/Exact_cover) using
+Algorithm X as described in the
+paper [_Dancing Links_](https://arxiv.org/abs/cs/0011047), by Donald Knuth, in
+_Millennial Perspectives in Computer Science_, P159, 2000.
+
+See [`exact-cover` on Hackage](https://hackage.haskell.org/package/exact-cover)
+for more information.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/cbits/dbg.h b/cbits/dbg.h
new file mode 100644
--- /dev/null
+++ b/cbits/dbg.h
@@ -0,0 +1,34 @@
+/*
+ * Reference: http://c.learncodethehardway.org/book/ex20.html
+ */
+
+#ifndef __dbg_h__
+#define __dbg_h__
+
+#include <stdio.h>
+#include <errno.h>
+#include <string.h>
+
+#ifdef NDEBUG
+#define debug(M, ...)
+#else
+#define debug(M, ...) fprintf(stderr, "DEBUG %s:%d: " M "\n", __FILE__, __LINE__, ##__VA_ARGS__)
+#endif
+
+#define clean_errno() (errno == 0 ? "None" : strerror(errno))
+
+#define log_err(M, ...) fprintf(stderr, "[ERROR] (%s:%d: errno: %s) " M "\n", __FILE__, __LINE__, clean_errno(), ##__VA_ARGS__)
+
+#define log_warn(M, ...) fprintf(stderr, "[WARN] (%s:%d: errno: %s) " M "\n", __FILE__, __LINE__, clean_errno(), ##__VA_ARGS__)
+
+#define log_info(M, ...) fprintf(stderr, "[INFO] (%s:%d) " M "\n", __FILE__, __LINE__, ##__VA_ARGS__)
+
+#define check(A, M, ...) if(!(A)) { log_err(M, ##__VA_ARGS__); errno=0; goto error; }
+
+#define sentinel(M, ...)  { log_err(M, ##__VA_ARGS__); errno=0; goto error; }
+
+#define check_mem(A) check((A), "Out of memory.")
+
+#define check_debug(A, M, ...) if(!(A)) { debug(M, ##__VA_ARGS__); errno=0; goto error; }
+
+#endif
diff --git a/cbits/dlx.c b/cbits/dlx.c
new file mode 100644
--- /dev/null
+++ b/cbits/dlx.c
@@ -0,0 +1,334 @@
+#include "dbg.h"
+#include <stdint.h>
+#include "dlx.h"
+
+/*
+ * Creates empty DLX matrix.
+ */
+struct cell *create_empty_matrix() {
+  struct cell *h = malloc(sizeof(struct cell));
+  if (h != NULL) {
+    h->L = h;
+    h->R = h;
+    h->U = NULL;
+    h->D = NULL;
+    h->C = NULL;
+  }
+  return h;
+}
+
+/*
+ * Frees all memory associated with given matrix.
+ */
+void free_matrix(struct cell *h) {
+  if (h == NULL) return;
+
+  if (h->C != NULL) {  /* Not a header cell, so make h point to the header
+                        * cell. */
+    h = h->C;
+    while (h->C != NULL) h = h->R;
+  }
+
+  /* Cells are freed from left-most column to right-most column, and within each
+   * column, top to bottom. */
+  struct cell *c = h->R;
+  while (c != h) {
+    struct cell *r = c->D;
+    while (r != c) {
+      r = r->D;
+      free(r->U);
+    }
+    c = c->R;
+    free(c->L);
+  }
+  free(h);
+}
+
+/*
+ * Inserts column. No-op if column exists. Search always move rightwards from
+ * starting cell (this is intentional).
+ *
+ * Input: xp: a pointer to a pointer to any valid cell.
+ *
+ * Return code: 1 if column already exists, 0 if insertion is successful, -1 for
+ * any other error. Input pointer is modified to point to inserted cell.
+ *
+ * Warning: matrix is in an inconsistent state if function fails.
+ */
+int add_column(struct cell **xp, column_name new_N) {
+  struct cell *c = *xp;
+  if (c->C != NULL) c = c->C;  /* c to point to column header if it doesn't
+                                * already do so */
+
+  /* Step 1: check that c is to the left of the target location */
+  column_name current_N = (c->C == NULL) ? COLUMN_NAME_MIN : c->N;
+  if (current_N > new_N) {  /* too far to the right */
+    while (c->C != NULL) c = c->R;  /* Move back to the header */
+    current_N = COLUMN_NAME_MIN;
+  }
+  check_debug(current_N <= new_N, "assertion failed: current_N <= new_N");
+
+  /* Step 2: move c such that it just crosses or touches the target location */
+  while (current_N < new_N) {
+    c = c->R;
+    current_N = (c->C == NULL) ? COLUMN_NAME_MAX : c->N;
+  }
+
+  /* Step 3: check case where column already exists */
+  if (current_N == new_N) {
+    *xp = c;
+    return 1;
+  }
+  check_debug(((c->L->C == NULL) ? COLUMN_NAME_MIN : c->L->N) < new_N && new_N < current_N,
+              "assertion failed: previous_N < new_N && new_N < current_N");
+
+  /* Step 4: Insert col to the left of c */
+  struct cell *new_c = malloc(sizeof(struct cell));
+  check_mem(new_c);
+  new_c->L = c->L;
+  new_c->R = c;
+  new_c->U = new_c;
+  new_c->D = new_c;
+  new_c->C = new_c;
+  new_c->S = 0;
+  new_c->N = new_N;
+  c->L->R = new_c;
+  c->L = new_c;
+  *xp = new_c;
+  return 0;
+
+error:
+  return -1;
+}
+
+/*
+ * Inserts a set.
+ *
+ * Input: xp: a pointer to a pointer to any valid cell. set: array should (but
+ * not necessarily) be in circularly increasing order for efficiency. e.g. "1 2
+ * 3" or "3 1 2" is okay, but "3 2 1" or "1 3 2" is not.
+ *
+ * Output: 0 for success, -1 for error. Input pointer to the cell pointer is
+ * modified to point to the last added cell.
+ *
+ * Warning: matrix is in an inconsistent state if function fails.
+ */
+int add_set(struct cell **xp, column_name *set, size_t n_set_elem) {
+  struct cell *c = *xp;
+
+  struct cell *y_prev, *y;
+  for (size_t i = 0; i < n_set_elem; ++i) {
+    check(add_column(&c, set[i]) >= 0, "Failed to add column header.");
+    /* note: c points to the column header after add_column */
+
+    y = malloc(sizeof(struct cell));
+    check_mem(y);
+    if (i == 0) {
+      y->L = y;
+      y->R = y;
+    } else {
+      y->L = y_prev;
+      y->R = y_prev->R;
+      y_prev->R->L = y;
+      y_prev->R = y;
+    }
+    y->U = c->U;
+    y->D = c;
+    c->U->D = y;
+    c->U = y;
+    y->C = c;
+
+    ++(c->S);
+    y_prev = y;
+  }
+
+  *xp = y;
+  return 0;
+
+error:
+  return -1;
+}
+
+/*
+ * Find column with the smallest number of elements.
+ *
+ * Input: header cell.
+ */
+struct cell *find_smallest_col(const struct cell *h) {
+  check(h->C == NULL, "Not header cell.");
+  struct cell *c = h->R;
+  struct cell *min_c;
+  size_t min_S = SIZE_MAX;
+  while (c != h) {
+    if (c->S < min_S) {
+      min_S = c->S;
+      min_c = c;
+    }
+    c = c->R;
+  }
+  return min_c;
+
+error:
+  return NULL;
+}
+
+/*
+ * Covers the column that a given cell is in. Does not free memory.
+ */
+void cover_col(struct cell *x) {
+  struct cell *c = x->C;
+  c->R->L = c->L;
+  c->L->R = c->R;
+
+  struct cell *i = c->D;
+  while (i != c) {
+    /* For each cell i in the column, cover the row. */
+    struct cell *j = i->R;
+    while (j != i) {
+      j->D->U = j->U;
+      j->U->D = j->D;
+      --j->C->S;
+      j = j->R;
+    }
+
+    i = i->D;
+  }
+}
+
+/*
+ * Uncovers the column that a given cell is in.
+ */
+void uncover_col(struct cell *x) {
+  struct cell *c = x->C;
+
+  /* Programming note: the uncovering is done in exactly the reverse order by
+   * which the column was previously covered! */
+  struct cell *i = c->U;
+  while (i != c) {
+    /* For each cell i in the column, uncover the row. */
+    struct cell *j = i->L;
+    while (j != i) {
+      ++j->C->S;
+      j->D->U = j;
+      j->U->D = j;
+      j = j->L;
+    }
+
+    i = i->U;
+  }
+
+  c->R->L = c;
+  c->L->R = c;
+}
+
+/*
+ * The Dancing Links X algorithm. Note: allocates memory for output array.
+ *
+ * Input: h: header cell. k: always start at 0.
+ *
+ * Output: pointers to cells whose rows are part of the solution.
+ *
+ * Return code: 0 = good. 1 = no solution. -1 = error.
+ */
+int dlx_solve(const struct cell *h, depth_t k, struct cell ***output_arr, size_t *output_arr_size) {
+  check(h->C == NULL, "Not header cell.");
+
+  if (h->R == h) {
+    *output_arr = malloc(sizeof(struct cell *) * k);
+    check_mem(*output_arr);
+    *output_arr_size = k;
+    return 0;
+  }
+
+  int ret = 1;
+  struct cell *c = find_smallest_col(h);
+  cover_col(c);
+
+  struct cell *candidate_row = c->D;
+  struct cell *j;  /* variable to traverse candidate_row */
+  while (candidate_row != c) {
+    j = candidate_row->R;
+    while (j != candidate_row) {
+      cover_col(j);
+      j = j->R;
+    }
+
+    ret = dlx_solve(h, k+1, output_arr, output_arr_size);
+
+    j = candidate_row->L;
+    while (j != candidate_row) {
+      uncover_col(j);
+      j = j->L;
+    }
+
+    if (ret == 0) {
+      (*output_arr)[k] = candidate_row;
+      break;
+    }
+
+    candidate_row = candidate_row->D;
+  }
+
+  uncover_col(c);
+  return ret;
+
+error:
+  return -1;
+}
+
+/*
+ * Helper function that wraps dlx_solve(). Note: allocates memory for output
+ * array.
+ *
+ * Input: header cell and the most probable number of elements each solution set
+ * will have. The latter will be used to allocate memory (but will be expanded
+ * if there are more elements).
+ */
+int solve(const struct cell *h, size_t suggested_set_size,
+          int ***set_covers, size_t **set_cover_sizes, size_t *n_sets) {
+  struct cell **ans = NULL;
+  size_t ans_len = 0;
+  int ret_code = dlx_solve(h, 0, &ans, &ans_len);
+
+  if (ret_code == 0) {
+    *n_sets = ans_len;
+    *set_covers = malloc(sizeof(int *) * ans_len);
+    check_mem(*set_covers);
+    *set_cover_sizes = malloc(sizeof(size_t) * ans_len);
+    check_mem(*set_cover_sizes);
+
+    struct cell *x;
+    size_t j, set_max_size;
+    for (size_t i = 0; i < ans_len; ++i) {
+      x = ans[i];
+      set_max_size = suggested_set_size;
+      (*set_covers)[i] = malloc(sizeof(int) * set_max_size);
+
+      j = 0;
+      do {
+        if (j >= set_max_size) {
+          set_max_size *= 2;
+          int *tmp = realloc((*set_covers)[i], sizeof(int) * set_max_size);
+          check_mem(tmp);
+          (*set_covers)[i] = tmp;
+        }
+
+        (*set_covers)[i][j] = x->C->N;
+        x = x->R;
+        ++j;
+      } while (x != ans[i]);
+
+      (*set_cover_sizes)[i] = j;
+    }
+  }
+
+  return ret_code;
+
+error:
+  for (size_t i = 0; i < ans_len; ++i) {
+    free((*set_covers)[i]);
+  }
+  free(*set_cover_sizes);
+  free(*set_covers);
+  return -1;
+}
diff --git a/cbits/dlx.h b/cbits/dlx.h
new file mode 100644
--- /dev/null
+++ b/cbits/dlx.h
@@ -0,0 +1,31 @@
+#ifndef _dlx_h
+#define _dlx_h
+
+#include <stdlib.h>
+#include <limits.h>
+
+typedef unsigned int depth_t;
+typedef int column_name;
+#define COLUMN_NAME_MAX INT_MAX
+#define COLUMN_NAME_MIN INT_MIN
+
+/*
+ * Column and data object
+ */
+struct cell {
+  struct cell *L;
+  struct cell *R;
+  struct cell *U;
+  struct cell *D;
+  struct cell *C;
+  size_t S;
+  column_name N;
+};
+
+struct cell *create_empty_matrix();
+void free_matrix(struct cell *h);
+int add_set(struct cell **xp, column_name *set, size_t n_set_elem);
+int dlx_solve(const struct cell *h, depth_t k, struct cell ***output_arr, size_t *output_arr_size);
+int solve(const struct cell *h, size_t suggested_set_size, int ***set_covers, size_t **set_cover_sizes, size_t *n_sets);
+
+#endif
diff --git a/exact-cover.cabal b/exact-cover.cabal
new file mode 100644
--- /dev/null
+++ b/exact-cover.cabal
@@ -0,0 +1,91 @@
+name: exact-cover
+version: 0.1.0.0
+cabal-version: >=1.10
+build-type: Simple
+license: BSD3
+license-file: LICENSE.txt
+copyright: (c) 2017, Arthur Lee
+maintainer: me@arthur.li
+homepage: https://github.com/arthurl/exact-cover
+bug-reports: https://github.com/arthurl/exact-cover/issues
+synopsis: Efficient exact cover solver.
+description:
+    Fast solver for exact set cover problems
+    (<http://en.wikipedia.org/wiki/Exact_cover>) using Algorithm X as described in
+    the paper /Dancing Links/, by Donald Knuth, in
+    /Millennial Perspectives in Computer Science/, P159, 2000
+    (<https://arxiv.org/abs/cs/0011047>).
+    .
+    To get started, see the documentation for the "Math.ExactCover" module below.
+    .
+    Build examples with @cabal install -fbuildExamples@ or @stack build --flag
+    exact-cover:buildExamples@. Examples include a Sudoku solver.
+category: Math, Algorithms
+author: Arthur Lee
+tested-with: GHC ==7.10.3 GHC ==8.0.1 GHC ==8.0.2
+extra-source-files:
+    README.md
+    cbits/*.h
+
+source-repository head
+    type: git
+    location: https://github.com/arthurl/exact-cover
+
+flag buildexamples
+    description:
+        Build example executables.
+    default: False
+    manual: True
+
+library
+    exposed-modules:
+        Math.ExactCover
+        Math.ExactCover.Internal.DLX
+    build-depends:
+        base >=4.6 && <4.10,
+        containers ==0.5.*
+    cc-options: -std=c99
+    c-sources:
+        cbits/dlx.c
+    default-language: Haskell2010
+    default-extensions: OverloadedStrings ViewPatterns
+    other-extensions: ForeignFunctionInterface CPP EmptyDataDecls
+    include-dirs: cbits
+    hs-source-dirs: src
+    ghc-options: -funbox-strict-fields -Wall
+
+executable sudoku
+    
+    if !flag(buildexamples)
+        buildable: False
+    main-is: Main.hs
+    build-depends:
+        exact-cover <0.2,
+        base <4.10,
+        containers <0.6,
+        vector <0.12,
+        safe <0.4,
+        boxes <0.2
+    default-language: Haskell2010
+    default-extensions: OverloadedStrings ViewPatterns
+    hs-source-dirs: examples/sudoku
+    other-modules:
+        Sudoku
+        Sudoku.Grid
+    ghc-options: -Wall
+
+test-suite tasty
+    type: exitcode-stdio-1.0
+    main-is: Main.hs
+    build-depends:
+        exact-cover <0.2,
+        base <4.10,
+        tasty <0.12,
+        tasty-hunit <0.10,
+        containers <0.6
+    default-language: Haskell2010
+    default-extensions: OverloadedStrings ViewPatterns
+    hs-source-dirs: tests
+    other-modules:
+        Math.ExactCover.Tests
+    ghc-options: -Wall
diff --git a/examples/sudoku/Main.hs b/examples/sudoku/Main.hs
new file mode 100644
--- /dev/null
+++ b/examples/sudoku/Main.hs
@@ -0,0 +1,18 @@
+module Main where
+
+import Math.ExactCover
+import Sudoku
+
+import Control.Monad (forever)
+
+
+main :: IO ()
+main = do
+  let f sudokuPuzzleStr =
+        case solve . toConstraints . readLineFormat $ sudokuPuzzleStr of
+          [] -> putStrLn "No solution."
+          soln@(_:_) -> print . convert2SolvedGrid $ soln
+  forever $ do
+    putStrLn "Enter puzzle in single line format:"
+    putStrLn "(e.g. 8..........36......7..9.2...5...7.......457.....1...3...1....68..85...1..9....4..)"
+    getLine >>= f
diff --git a/examples/sudoku/Sudoku.hs b/examples/sudoku/Sudoku.hs
new file mode 100644
--- /dev/null
+++ b/examples/sudoku/Sudoku.hs
@@ -0,0 +1,156 @@
+-- |
+
+module Sudoku
+  (
+    -- * Types
+    Digit(..), allDigits
+  , SRow(..), allRows
+  , SCol(..), allColumns
+  , SBox(..), box
+  , SConstraint(..)
+
+    -- * Map Sudoku to exact-cover problem
+  , toConstraints
+  , convert2SolvedGrid
+
+    -- * Reading Sudokus
+  , readLineFormat
+  , readPEulerGrid
+  )
+where
+
+import Sudoku.Grid
+
+import Control.Arrow (second)
+import Safe (readMay, headMay)
+import qualified Data.Vector.Generic as V (generate)
+import Data.Map (Map)
+import qualified Data.Map.Strict as Map (lookup, fromList)
+import Data.Set (Set)
+import qualified Data.Set as Set (fromList)
+
+
+data Digit = D1 | D2 | D3 | D4 | D5 | D6 | D7 | D8 | D9
+  deriving (Eq, Ord, Enum)
+
+instance Show Digit where
+  show D1 = "1"
+  show D2 = "2"
+  show D3 = "3"
+  show D4 = "4"
+  show D5 = "5"
+  show D6 = "6"
+  show D7 = "7"
+  show D8 = "8"
+  show D9 = "9"
+
+instance Read Digit where
+  readsPrec _ ('1':rest) = [(D1, rest)]
+  readsPrec _ ('2':rest) = [(D2, rest)]
+  readsPrec _ ('3':rest) = [(D3, rest)]
+  readsPrec _ ('4':rest) = [(D4, rest)]
+  readsPrec _ ('5':rest) = [(D5, rest)]
+  readsPrec _ ('6':rest) = [(D6, rest)]
+  readsPrec _ ('7':rest) = [(D7, rest)]
+  readsPrec _ ('8':rest) = [(D8, rest)]
+  readsPrec _ ('9':rest) = [(D9, rest)]
+  readsPrec _ _ = []
+
+allDigits :: [Digit]
+allDigits = [toEnum x | x <- [0..8]]
+
+data SRow = R0 | R1 | R2 | R3 | R4 | R5 | R6 | R7 | R8
+  deriving (Show, Eq, Ord, Enum)
+
+allRows :: [SRow]
+allRows = [toEnum x | x <- [0..8]]
+
+data SCol = C0 | C1 | C2 | C3 | C4 | C5 | C6 | C7 | C8
+  deriving (Show, Eq, Ord, Enum)
+
+allColumns :: [SCol]
+allColumns = [toEnum x | x <- [0..8]]
+
+data SBox = B0 | B1 | B2 | B3 | B4 | B5 | B6 | B7 | B8
+  deriving (Show, Eq, Ord, Enum)
+
+-- | Given the row and column, get the box.
+box :: SRow -> SCol -> SBox
+box (fromEnum -> r) (fromEnum -> c) =
+  toEnum $ (r `div` 3) * 3 + (c `div` 3)
+
+data SConstraint = RowColumn !SRow !SCol
+                  | RowNumber !SRow !Digit
+                  | ColumnNumber !SCol !Digit
+                  | BoxNumber !SBox !Digit
+  deriving (Show, Eq, Ord)
+
+instance Enum SConstraint where
+  fromEnum (RowColumn r c) = fromEnum r * 9 + fromEnum c
+  fromEnum (RowNumber r n) = 81 + fromEnum r * 9 + fromEnum n
+  fromEnum (ColumnNumber c n) = 162 + fromEnum c * 9 + fromEnum n
+  fromEnum (BoxNumber b n) = 243 + fromEnum b * 9 + fromEnum n
+
+  toEnum x = case second (`divMod` 9) $ x `divMod` 81 of
+    (a, (b, c)) | a == 0 -> RowColumn (toEnum b) (toEnum c)
+                | a == 1 -> RowNumber (toEnum b) (toEnum c)
+                | a == 2 -> ColumnNumber (toEnum b) (toEnum c)
+                | a == 3 -> BoxNumber (toEnum b) (toEnum c)
+                | otherwise -> error "toEnum instance error."
+
+-- | Maps the Sudoku problem to a set of constraints following
+-- https://en.wikipedia.org/wiki/Exact_cover#Sudoku.
+toConstraints :: Grid (Maybe Digit)
+              -> Map (Set SConstraint) (SRow, SCol, Digit)
+toConstraints (toListList -> g) =
+  let expanded :: [(SRow, SCol, Maybe Digit)]
+      expanded = concatMap (\(r, xs) -> zip3 (repeat r) allColumns xs)
+                 . zip allRows $ g
+  in Map.fromList $ concatMap f expanded
+  where
+    f :: (SRow, SCol, Maybe Digit)
+      -> [((Set SConstraint), (SRow, SCol, Digit))]
+    f (r, c, Just n) = [clueExist r c n]
+    f (r, c, Nothing) = clueNotExist r c
+
+    clueExist :: SRow -> SCol -> Digit
+              -> ((Set SConstraint), (SRow, SCol, Digit))
+    clueExist r c n = let constraint = Set.fromList [ RowColumn r c
+                                                    , RowNumber r n
+                                                    , ColumnNumber c n
+                                                    , BoxNumber (box r c) n
+                                                    ]
+                      in (constraint, (r, c, n))
+
+    clueNotExist :: SRow -> SCol
+                 -> [((Set SConstraint), (SRow, SCol, Digit))]
+    clueNotExist r c = map (clueExist r c) allDigits
+
+convert2SolvedGrid :: [(SRow, SCol, Digit)] -> Grid Digit
+convert2SolvedGrid xs =
+  let m = Map.fromList . map (\(r,c,n) -> ((r,c),n)) $ xs
+  in Grid $ V.generate 9 $ \r ->
+            V.generate 9 $ \c ->
+              case Map.lookup (toEnum r, toEnum c) m of
+                Just n -> n
+                Nothing -> error "convert2SolvedGrid error."
+
+-- | Parse single line Sudoku format.
+readLineFormat :: String -> Grid (Maybe Digit)
+readLineFormat (headMay . lines -> Just str) | length str == 81 =
+  fromListList . chunksOf 9 . map (\c -> readMay [c]) $ str
+  where
+  chunksOf :: Int -> [a] -> [[a]]
+  chunksOf _ [] = [[]]
+  chunksOf n xs@(_:_) = let (y, rest) = splitAt n xs
+                        in y : chunksOf n rest
+readLineFormat _ = error "Invalid Sudoku grid."
+
+-- | Parse project Euler Sudoku format.
+readPEulerGrid :: String -> Grid (Maybe Digit)
+readPEulerGrid (lines -> strS) | length strS >= 9 =
+  fromListList $ map readPEulerLine strS
+  where
+    readPEulerLine :: String -> [Maybe Digit]
+    readPEulerLine (take 9 -> strLn) = map (\c -> readMay [c]) strLn
+readPEulerGrid _ = error "Invalid Sudoku grid."
diff --git a/examples/sudoku/Sudoku/Grid.hs b/examples/sudoku/Sudoku/Grid.hs
new file mode 100644
--- /dev/null
+++ b/examples/sudoku/Sudoku/Grid.hs
@@ -0,0 +1,70 @@
+-- |
+
+module Sudoku.Grid
+  (
+    -- * Types
+    Grid(..)
+
+    -- * Conversion
+  , fromListList
+  , toListList
+
+    -- * Printing
+  , prettyPrintGrid
+
+    -- * Access
+  , gridElem
+
+  )
+where
+
+import Data.Functor.Identity (Identity(..))
+import Data.List (transpose)
+import qualified Data.Vector as VB (Vector)
+import qualified Data.Vector.Generic as V (fromList, toList)
+import Data.Vector.Generic ((!))
+import qualified Text.PrettyPrint.Boxes as Box
+
+
+newtype Grid a = Grid (VB.Vector (VB.Vector a))
+
+instance (Show a) => Show (Grid a) where
+  show = runIdentity . prettyPrintGrid (Identity . show)
+
+fromListList :: [[a]] -> Grid a
+fromListList = Grid . V.fromList . map V.fromList
+
+toListList :: Grid a -> [[a]]
+toListList (Grid vv) = map V.toList . V.toList $ vv
+
+-- | Prints matrix.
+prettyPrintGrid :: (Monad m) => (a -> m String) -> Grid a -> m String
+prettyPrintGrid printCell (Grid g) = do
+  g2 <- mapM (mapM printCell) . transpose . map V.toList $ V.toList g
+  let g3 = map (Box.vcat Box.center1 . map Box.text) g2
+      colCount = map Box.cols g3
+      gStr = Box.render $ Box.hsep 1 Box.center1 $ intersperse3 verticleSep g3
+      gStr2 = unlines . intersperse3 (horizontalSep colCount) . lines $ gStr
+  pure gStr2
+  where
+    verticleSep :: Box.Box
+    verticleSep = Box.vcat Box.center1 . replicate 9 $ Box.text "|"
+    horizontalSep :: [Int] -> String
+    horizontalSep (c1:c2:c3:c4:c5:c6:c7:c8:c9:_) =
+      replicate (c1+c2+c3+3) '-'
+        ++ '+' : replicate (c4+c5+c6+4) '-'
+        ++ '+' : replicate (c7+c8+c9+3) '-'
+    horizontalSep _ = undefined
+
+    intersperse3 :: a -> [a] -> [a]
+    intersperse3 a (x1:x2:x3:xs@(_:_)) = x1:x2:x3:a:(intersperse3 a xs)
+    intersperse3 _ xs = xs
+
+data Coordinate = Coordinate !Int !Int
+
+instance Show Coordinate where
+  show (Coordinate r c) = '(' : show r ++ ',' : show c ++ ")"
+
+-- | Get the corresponding element of a grid.
+gridElem :: Grid a -> Coordinate -> a
+gridElem (Grid g) (Coordinate r c) = g ! r ! c
diff --git a/src/Math/ExactCover.hs b/src/Math/ExactCover.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/ExactCover.hs
@@ -0,0 +1,138 @@
+-- |
+--
+-- This module provides an efficient solver for exact set cover problems
+-- (<http://en.wikipedia.org/wiki/Exact_cover>) using Algorithm X as described
+-- in the paper /Dancing Links/, by Donald Knuth, in
+-- /Millennial Perspectives in Computer Science/, P159, 2000
+-- (<https://arxiv.org/abs/cs/0011047>).
+--
+-- For a quick start, go straight to the 'solve' function.
+
+module Math.ExactCover
+  (
+    -- * Mathematical definition
+    -- $def
+
+    -- * Simple interface
+    solve
+
+    -- * Types
+  , ExactCoverProblem
+
+    -- * Construction
+  , transform
+
+    -- * Solvers
+  , solveEC
+  )
+where
+
+import Math.ExactCover.Internal.DLX
+
+import Foreign
+import System.IO.Unsafe (unsafePerformIO)
+import Control.Concurrent (MVar, newMVar, withMVar)
+import Control.Monad (forM, forM_, when)
+import Data.Map (Map)
+import qualified Data.Map.Strict as Map (lookup, keys)
+import Data.Set (Set)
+import qualified Data.Set as Set (size, toList, fromList)
+
+
+-- | Basic type that represents the exact cover problem.
+data ExactCoverProblem setlabel a = ExactCoverProblem
+  { ec_sets :: !(Map (Set a) setlabel)  -- ^ Associates each set with a label.
+  , ec_dlx :: MVar (ForeignPtr DLXMatrix)  -- ^ Dancing links representation.
+  }
+
+-- PROGRAMMING NOTE: The Enum typeclass is needed as the set objects are
+-- internally represented as integers in the C portion of the library.
+
+-- | Constructs an 'ExactCoverProblem' given a collection of subsets \\(
+-- \\mathcal{S} \\) as represented by a 'Map' between each subset and its label.
+-- The set \\( \\mathcal{X} \\) over which the exact cover is to be found is
+-- assumed to be the union of the given collection of subsets \\( \\mathcal{S}
+-- \\).
+transform :: (Enum a)
+          => Map (Set a) setlabel
+          -> ExactCoverProblem setlabel a
+transform m = unsafePerformIO $ do
+  dlxHead <- newForeignPtr c_free_matrix =<< c_create_empty_matrix
+  withForeignPtr dlxHead $ \hPtr ->
+    forM_ (Map.keys m) $ \k ->
+      withArray (fromIntegral . fromEnum <$> Set.toList k) $ \constrainPtr ->
+        with hPtr $ \hPPtr -> do
+          ret <- c_add_set hPPtr constrainPtr (fromIntegral $ Set.size k)
+          when (ret /= 0) $ error "Could not add constraint."
+  dlxM <- newMVar dlxHead
+  pure $ ExactCoverProblem { ec_sets = m
+                           , ec_dlx = dlxM
+                           }
+
+-- | Solves the given 'ExactCoverProblem', returning the labels of the subsets
+-- that form the exact cover.
+solveEC :: (Enum a, Ord a)
+        => ExactCoverProblem setlabel a
+        -> [setlabel]
+solveEC ExactCoverProblem{ ec_sets = setLabels, ec_dlx = dlxM } = unsafePerformIO $
+  withMVar dlxM $ \dlxHead ->
+  withForeignPtr dlxHead $ \hPtr ->
+  alloca $ \setCoversPtrPtrPtr ->
+  alloca $ \setCoverSizesPtrPtr ->
+  alloca $ \nSetsPtr -> do
+    ret <- c_solve hPtr 4 setCoversPtrPtrPtr setCoverSizesPtrPtr nSetsPtr
+    case () of
+      _ | ret < 0 -> error ""
+        | ret == 0 -> do
+            -- Retrieve number of result sets.
+            nSets <- fromIntegral <$> peek nSetsPtr
+
+            -- Retrieve list of set sizes.
+            setCoverSizesPtr <- peek setCoverSizesPtrPtr
+            setCoverSizes <- pure . map fromIntegral
+                             =<< peekArray nSets setCoverSizesPtr
+            free setCoverSizesPtr
+
+            -- Retrieve list of pointers to the result sets.
+            setCoversPtrPtr <- peek setCoversPtrPtrPtr
+            setCoversPtr <- peekArray nSets setCoversPtrPtr
+            free setCoversPtrPtr
+
+            -- Retrieve the result sets.
+            forM (zip setCoverSizes setCoversPtr) $ \(setSize, setPtr) -> do
+              setC <- Set.fromList . map (toEnum . fromIntegral)
+                      <$> peekArray setSize setPtr
+              free setPtr
+              pure $ case Map.lookup setC setLabels of
+                Nothing -> error "Constrain set not in map."
+                Just label -> label
+
+        | ret == 1 -> pure mempty
+        | otherwise -> error "Unknown error code."
+
+-- | Given a collection of subsets \\( \\mathcal{S} \\), represented by a 'Map'
+-- between each subset (of type @'Set' a@) and its label, returns a list of
+-- labels that represents the exact cover \\( \\mathcal{S}^{*} \\).
+--
+-- Example: To find the exact cover of the collection of subsets
+-- \\( \\left\\{\\left\\{2,4,5\\right\\}, \\left\\{0,3,6\\right\\},
+-- \\left\\{1,2,5\\right\\}, \\left\\{0,3\\right\\}, \\left\\{1,6\\right\\},
+-- \\left\\{3,4,6\\right\\}\\right\\} \\),
+--
+-- > solve (Map.fromList [ (Set.fromList [2,4,5], 'A')
+-- >                     , (Set.fromList [0,3,6], 'B')
+-- >                     , (Set.fromList [1,2,5], 'C')
+-- >                     , (Set.fromList [0,3], 'D')
+-- >                     , (Set.fromList [1,6], 'E')
+-- >                     , (Set.fromList [3,4,6], 'F')
+-- >                     ] :: Map (Set Int) Char)
+-- > == "DAE"
+solve :: (Enum a, Ord a) => Map (Set a) setlabel -> [setlabel]
+solve = solveEC . transform
+
+-- $def
+--
+-- Given a collection \\( \\mathcal{S} \\) of subsets of a set \\( \\mathcal{X}
+-- \\), an exact cover is a subcollection \\( \\mathcal{S}^{*} \\) of \\(
+-- \\mathcal{S} \\) such that each element in \\( \\mathcal{X} \\) is contained
+-- in exactly one subset in \\( \\mathcal{S}^{*} \\) (from wikipedia).
diff --git a/src/Math/ExactCover/Internal/DLX.hsc b/src/Math/ExactCover/Internal/DLX.hsc
new file mode 100644
--- /dev/null
+++ b/src/Math/ExactCover/Internal/DLX.hsc
@@ -0,0 +1,31 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+{-# LANGUAGE CPP                      #-}
+{-# LANGUAGE EmptyDataDecls           #-}
+
+-- |
+-- = WARNING
+--
+-- This module is considered __internal__.
+--
+
+module Math.ExactCover.Internal.DLX where
+
+import Foreign
+import Foreign.C.Types
+
+
+#include <dlx.h>
+
+data DLXMatrix
+
+foreign import ccall unsafe "dlx.h solve"
+  c_solve :: Ptr DLXMatrix -> CSize -> Ptr (Ptr (Ptr CInt)) -> Ptr (Ptr CSize) -> Ptr CSize -> IO CInt
+
+foreign import ccall unsafe "dlx.h &free_matrix"
+  c_free_matrix :: FinalizerPtr DLXMatrix
+
+foreign import ccall unsafe "dlx.h create_empty_matrix"
+  c_create_empty_matrix :: IO (Ptr DLXMatrix)
+
+foreign import ccall unsafe "dlx.h add_set"
+  c_add_set :: Ptr (Ptr DLXMatrix) -> Ptr CInt -> CSize -> IO CInt
diff --git a/tests/Main.hs b/tests/Main.hs
new file mode 100644
--- /dev/null
+++ b/tests/Main.hs
@@ -0,0 +1,13 @@
+-- |
+
+module Main where
+
+import Test.Tasty (defaultMain, testGroup)
+
+import qualified Math.ExactCover.Tests
+
+
+main :: IO ()
+main = defaultMain $ testGroup "Tests"
+    [ Math.ExactCover.Tests.tests
+    ]
diff --git a/tests/Math/ExactCover/Tests.hs b/tests/Math/ExactCover/Tests.hs
new file mode 100644
--- /dev/null
+++ b/tests/Math/ExactCover/Tests.hs
@@ -0,0 +1,35 @@
+-- |
+
+module Math.ExactCover.Tests where
+
+import Math.ExactCover
+
+import Test.Tasty (TestTree, testGroup)
+import qualified Test.Tasty.HUnit as U
+import Test.Tasty.HUnit ((@?=))
+
+import qualified Data.Map.Strict as Map (fromList)
+import qualified Data.Set as Set (fromList)
+
+
+tests :: TestTree
+tests = testGroup "Math.ExactCover"
+  [ tests_solve
+  ]
+
+
+tests_solve :: TestTree
+tests_solve = testGroup "solve"
+  [ U.testCase "regular test case" case_solve_reg1
+  ]
+
+case_solve_reg1 :: U.Assertion
+case_solve_reg1 =
+  ( Set.fromList . solve $ Map.fromList [ (Set.fromList [2,4,5::Int], 'A')
+                                        , (Set.fromList [0,3,6], 'B')
+                                        , (Set.fromList [1,2,5], 'C')
+                                        , (Set.fromList [0,3], 'D')
+                                        , (Set.fromList [1,6], 'E')
+                                        , (Set.fromList [3,4,6], 'F')
+                                        ]
+  ) @?= Set.fromList "ADE"