diff --git a/bed-and-breakfast.cabal b/bed-and-breakfast.cabal
--- a/bed-and-breakfast.cabal
+++ b/bed-and-breakfast.cabal
@@ -1,5 +1,5 @@
 Name:           bed-and-breakfast
-Version:        0.3.2
+Version:        0.4
 Synopsis:       Efficient Matrix operations in 100% Haskell.
 Description:    Efficient Matrix operations in 100% Haskell.
                 .
@@ -48,6 +48,10 @@
                 .
                 [@v0.3.2@] Numeric.Matrix.Sugar was not mentioned in the
                     cabal file. Improved test suite. Improved documentation.
+                .
+                [@v0.4@] Fixed a bug regarding @empty@ and @fromList@.
+                    Use unsafe operations where it is safe for speed.
+                    Added RULES. Added an instance for binary.
 
 
 License:        MIT
@@ -70,6 +74,7 @@
     Build-Depends:      base >= 4.5 && < 5,
                         deepseq >= 1.3,
                         array >= 0.4,
+                        binary >= 0.5,
                         template-haskell >= 2.7
     Hs-Source-Dirs:     src
 
@@ -80,4 +85,5 @@
     build-depends:  base >= 4.5 && < 5,
                     bed-and-breakfast == 0.3.2,
                     QuickCheck >= 2.4.2
+
 
diff --git a/src/Numeric/Matrix.hs b/src/Numeric/Matrix.hs
--- a/src/Numeric/Matrix.hs
+++ b/src/Numeric/Matrix.hs
@@ -59,26 +59,22 @@
 import Control.Monad
 import Control.Monad.ST
 
---import Control.Concurrent
---import Control.Concurrent.MVar
---import qualified Control.Monad.Parallel as MPar
-
 import Data.Function (on)
 import Data.Ratio
 import Data.Complex
 import Data.Maybe
---import Data.Foldable (Foldable)
---import qualified Data.Foldable as F
+import Data.Int
+import Data.Word
 
 import qualified Data.List as L
 import Data.Array.IArray
 import Data.Array.MArray
 import Data.Array.Unboxed
 import Data.Array.ST
---import Data.Array.IO
+import Data.Array.Base (unsafeRead, unsafeWrite)
 import Data.STRef
+import Data.Binary
 
---import System.IO.Unsafe
 import qualified Data.Array.Unsafe as U
 
 import Data.Typeable
@@ -86,6 +82,7 @@
 import Prelude hiding (any, all, read, map)
 import qualified Prelude as P
 
+
 -- | Matrices are represented by a type which fits best the component type.
 -- For example a @Matrix Double@ is represented by unboxed arrays,
 -- @Matrix Integer@ by boxed arrays.
@@ -95,9 +92,9 @@
 -- example you can not compute the inverse matrix of a @Matrix Int@.
 --
 -- Every matrix (regardless of the component type) has instances for
--- 'Show', 'Read', 'Num', 'Fractional', 'Eq', 'Typeable', and 'NFData'.
--- This means that you can use arithmetic operations like '+', '*', and
--- '/', as well as functions like 'show', 'read', or 'typeOf'.
+-- 'Show', 'Read', 'Num', 'Fractional', 'Eq', 'Typeable', 'Binary',
+-- and 'NFData'. This means that you can use arithmetic operations like
+-- '+', '*', and '/', as well as functions like 'show', 'read', or 'typeOf'.
 --
 -- [@Show (Matrix e)@]
 -- Note that a Show instance for the component type @e@ must exist.
@@ -132,7 +129,11 @@
 -- monad (see the @monad-par@ package for details).
 --
 -- [@Typeable (Matrix e)@]
--- Allows you to use matrices as 'Data.Dynamic' values.
+-- Allows you to use matrices as 'Dynamic' values.
+--
+-- [@Binary (Matrix e)@]
+-- Serialize and unserialize matrices using the @binary@ package.
+-- See @encode@ and @decode@.
 data family Matrix e
 
 data instance Matrix Int
@@ -191,7 +192,25 @@
 instance (MatrixElement e) => NFData (Matrix e) where
     rnf matrix = matrix `deepseq` ()
 
+instance (MatrixElement e, Binary e) => Binary (Matrix e) where
+    
+    put m = do
+        let (rows, cols) = dimensions m
+        
+        put rows >> put cols
 
+        forM_ [1..rows] $ \i -> do
+            forM_ [1..cols] $ \j -> do
+                put (m `at` (i,j))
+
+    get = do
+        rows <- get :: Get Int
+        cols <- get :: Get Int
+
+        forM [1..rows] (const (forM [1..cols] (const get)))
+            >>= return . fromList
+
+
 (<|>) :: MatrixElement e => Matrix e -> Matrix e -> Matrix e
 -- ^ Joins two matrices horizontally.
 --
@@ -268,19 +287,22 @@
 class Division e where
     divide :: e -> e -> e
 
-instance Division Int    where divide = quot
--- instance Division Int8   where divide = quot
--- instance Division Int16  where divide = quot
--- instance Division Int32  where divide = quot
--- instance Division Int64  where divide = quot
+instance Division Int     where divide = quot
+instance Division Int8    where divide = quot
+instance Division Int16   where divide = quot
+instance Division Int32   where divide = quot
+instance Division Int64   where divide = quot
+instance Division Word8   where divide = quot
+instance Division Word16  where divide = quot
+instance Division Word32  where divide = quot
+instance Division Word64  where divide = quot
 instance Division Integer where divide = quot
-instance Division Float  where divide = (/)
-instance Division Double where divide = (/)
+instance Division Float   where divide = (/)
+instance Division Double  where divide = (/)
 instance Integral a => Division (Ratio a) where divide = (/)
 instance RealFloat a => Division (Complex a) where divide = (/)
 
 
-
 class (Eq e, Num e) => MatrixElement e where
 
     -- | Creates a matrix of the given size using a generator
@@ -366,12 +388,14 @@
     -- Returns 0 if the matrix is not a square matrix.
     det       :: Matrix e -> e
 
-    -- | Flips rows and columns.
+    -- | Flip rows and columns.
     --
     -- > 1 8 9                1 2 3
     -- > 2 1 8  --transpose-> 8 1 2
     -- > 3 2 1                9 8 1 
     transpose :: Matrix e -> Matrix e
+
+    -- | Compute the rank of a matrix.
     rank      :: Matrix e -> e
     trace     :: Matrix e -> [e]
 
@@ -380,15 +404,15 @@
     adjugate :: MatrixElement e => Matrix e -> Matrix e
     minorMatrix :: MatrixElement e => Matrix e -> (Int, Int) -> Matrix e
 
-    -- | Applies a function on every component in the matrix.
+    -- | Apply a function on every component in the matrix.
     map :: MatrixElement f => (e -> f) -> Matrix e -> Matrix f
 
-    -- | Applies a predicate on every component in the matrix
-    -- and returns True if all components satisfy it.
+    -- | Apply a predicate on every component in the matrix
+    -- and returns True iff all components satisfy it.
     all :: (e -> Bool) -> Matrix e -> Bool
 
-    -- | Applies a predicate on every component in the matrix
-    -- and returns True if one or more components satisfy it.
+    -- | Apply a predicate on every component in the matrix
+    -- and return True if one or more components satisfy it.
     any :: (e -> Bool) -> Matrix e -> Bool
 
     mapWithIndex :: MatrixElement f => ((Int, Int) -> e -> f) -> Matrix e -> Matrix f
@@ -551,7 +575,7 @@
           => (Int -> Int -> a Int (u Int e) -> matrix e) -> [[e]] -> matrix e
 _fromList c xs =
     let lengths = P.map length xs
-        numCols = foldl1 min lengths
+        numCols = if null lengths then 0 else foldl1 min lengths
         numRows = length lengths
         
     in  c numRows numCols
@@ -602,21 +626,13 @@
          => (Int, Int) -> e -> ST s ((STUArray s) Int e)
 arraySTU = newArray
 
---arrayIO :: MArray IOArray e IO
---        => (Int, Int) -> e -> IO (IOArray Int e)
---arrayIO = newArray
 
---arrayIOU :: MArray IOUArray e IO
---         => (Int, Int) -> e -> IO (IOUArray Int e)
---arrayIOU = newArray
-
-
 tee :: Monad m => (b -> m a) -> b -> m b
 tee f x = f x >> return x
 
 read :: (MArray a1 b m, MArray a (a1 Int b) m) =>
                        a Int (a1 Int b) -> Int -> Int -> m b
-read a i j = readArray a i >>= flip readArray j
+read a i j = unsafeRead a i >>= flip unsafeRead j
 
 
 _inv :: (IArray a e, MArray (u s) e (ST s), Fractional e, Ord e, Show e)
@@ -628,15 +644,15 @@
         n = 2*m
 
         swap a i j = do
-            tmp <- readArray a i
-            readArray a j >>= writeArray a i
-            writeArray a j tmp
+            tmp <- unsafeRead a i
+            unsafeRead a j >>= unsafeWrite a i
+            unsafeWrite a j tmp
 
     okay <- newSTRef True
 
     a <- augment mkArrayST mat
 
-    flip mapM_ [1..m] $ \k -> do
+    forM_ [1..m] $ \k -> do
         iPivot <- zip [k..m] <$> mapM (\i -> abs <$> read a i k) [k..m]
                     >>= return . fst . L.maximumBy (compare `on` snd)
 
@@ -645,38 +661,38 @@
 
             swap a iPivot k
 
-            flip mapM_ [k+1..m] $ \i -> do
-                a_i <- readArray a i
-                a_k <- readArray a k
-                flip mapM_ [k+1..n] $ \j -> do
-                    a_ij <- readArray a_i j
-                    a_kj <- readArray a_k j
-                    a_ik <- readArray a_i k
-                    writeArray a_i j (a_ij - a_kj * (a_ik / p))
-                writeArray a_i k 0
+            forM_ [k+1..m] $ \i -> do
+                a_i <- unsafeRead a i
+                a_k <- unsafeRead a k
+                forM_ [k+1..n] $ \j -> do
+                    a_ij <- unsafeRead a_i j
+                    a_kj <- unsafeRead a_k j
+                    a_ik <- unsafeRead a_i k
+                    unsafeWrite a_i j (a_ij - a_kj * (a_ik / p))
+                unsafeWrite a_i k 0
 
     invertible <- readSTRef okay
 
     if invertible then
       do
-        flip mapM_ [ m - v | v <- [0..m-1] ] $ \i -> do
-            a_i <- readArray a i
-            p   <- readArray a_i i
-            writeArray a_i i 1
-            flip mapM_ [i+1..n] $ \j -> do
-                readArray a_i j >>= writeArray a_i j . (/ p)
+        forM_ [ m - v | v <- [0..m-1] ] $ \i -> do
+            a_i <- unsafeRead a i
+            p   <- unsafeRead a_i i
+            unsafeWrite a_i i 1
+            forM_ [i+1..n] $ \j -> do
+                unsafeRead a_i j >>= unsafeWrite a_i j . (/ p)
 
             unless (i == m) $ do
-                flip mapM_ [i+1..m] $ \k -> do
-                    a_k <- readArray a k
-                    p   <- readArray a_i k
+                forM_ [i+1..m] $ \k -> do
+                    a_k <- unsafeRead a k
+                    p   <- unsafeRead a_i k
 
-                    flip mapM_ [k..n] $ \j -> do
-                        a_ij <- readArray a_i j
-                        a_kj <- readArray a_k j
-                        writeArray a_i j (a_ij - p * a_kj)
+                    forM_ [k..n] $ \j -> do
+                        a_ij <- unsafeRead a_i j
+                        a_kj <- unsafeRead a_k j
+                        unsafeWrite a_i j (a_ij - p * a_kj)
 
-        mapM (\i -> readArray a i >>= getElems
+        mapM (\i -> unsafeRead a i >>= getElems
                         >>= return . listArray (1, m) . drop m) [1..m]
             >>= return . Just . listArray (1, m)
 
@@ -691,16 +707,16 @@
         n = snd $ bounds (mat ! 1)
 
         swap a i j = do
-            tmp <- readArray a i
-            readArray a j >>= writeArray a i
-            writeArray a j tmp
+            tmp <- unsafeRead a i
+            unsafeRead a j >>= unsafeWrite a i
+            unsafeWrite a j tmp
 
     a <- thaws mat >>= arrays
 
     ixPivot <- newSTRef 1
     prevR   <- newSTRef 1
 
-    flip mapM_ [1..n] $ \k -> do
+    forM_ [1..n] $ \k -> do
         pivotRow <- readSTRef ixPivot
 
         switchRow <- mapM (\i -> read a i k) [pivotRow .. m]
@@ -710,17 +726,17 @@
             let ix = fromJust switchRow + pivotRow
             when (pivotRow /= ix) (swap a pivotRow ix)
 
-            a_p   <- readArray a k
-            pivot <- readArray a_p k
+            a_p   <- unsafeRead a k
+            pivot <- unsafeRead a_p k
             prev  <- readSTRef prevR
             
-            flip mapM_ [pivotRow+1..m] $ \i -> do
-                a_i <- readArray a i
-                flip mapM_ [k+1..n] $ \j -> do
-                    a_ij <- readArray a_i j
-                    a_ik <- readArray a_i k
-                    a_pj <- readArray a_p j
-                    writeArray a_i j ((pivot * a_ij - a_ik * a_pj)
+            forM_ [pivotRow+1..m] $ \i -> do
+                a_i <- unsafeRead a i
+                forM_ [k+1..n] $ \j -> do
+                    a_ij <- unsafeRead a_i j
+                    a_ik <- unsafeRead a_i k
+                    a_pj <- unsafeRead a_p j
+                    unsafeWrite a_i j ((pivot * a_ij - a_ik * a_pj)
                                         `divide` prev)
 
             writeSTRef ixPivot (pivotRow + 1)
@@ -742,7 +758,7 @@
     signR  <- newSTRef 1
     pivotR <- newSTRef 1
 
-    flip mapM_ [1..size] $ \k -> do
+    forM_ [1..size] $ \k -> do
         sign <- readSTRef signR
         unless (sign == 0) $ do
 
@@ -750,7 +766,7 @@
             pivot <- read a k k >>= tee (writeSTRef pivotR)
 
             when (pivot == 0) $ do
-                s <- flip mapM [(k+1)..size] $ \r -> do
+                s <- forM [(k+1)..size] $ \r -> do
                     a_rk <- read a r k
                     if a_rk == 0 then return 0 else return r
                 let sf = filter (>0) s
@@ -758,9 +774,9 @@
                 when (not $ null sf) $ do
                     let sw = head sf
 
-                    row <- readArray a sw
-                    readArray a k >>= writeArray a sw
-                    writeArray a k row
+                    row <- unsafeRead a sw
+                    unsafeRead a k >>= unsafeWrite a sw
+                    unsafeWrite a k row
 
                     read a k k >>= writeSTRef pivotR
                     readSTRef signR >>= writeSTRef signR . negate
@@ -770,13 +786,13 @@
             sign' <- readSTRef signR
             unless (sign' == 0) $ do
                 pivot' <- readSTRef pivotR
-                flip mapM_ [(k+1)..size] $ \i -> do
-                    a_i <- readArray a i
-                    flip mapM [(k+1)..size] $ \j -> do
-                        a_ij <- readArray a_i j
-                        a_ik <- readArray a_i k
+                forM_ [(k+1)..size] $ \i -> do
+                    a_i <- unsafeRead a i
+                    forM [(k+1)..size] $ \j -> do
+                        a_ij <- unsafeRead a_i j
+                        a_ik <- unsafeRead a_i k
                         a_kj <- read a k j
-                        writeArray a_i j ((pivot' * a_ij - a_ik * a_kj) `divide` prev)
+                        unsafeWrite a_i j ((pivot' * a_ij - a_ik * a_kj) `divide` prev)
 
     liftM2 (*) (readSTRef pivotR) (readSTRef signR)
 
@@ -787,21 +803,7 @@
                 colsB = numCols b
             in  matrix (rowsA, colsB) (\(i,j) -> L.foldl' (+) 0 [a `at` (i, k) * b `at` (k, j) | k <- [1..rowsB]])
 
-{- Old function: TODO benchmark against new one
 
-_matrix_ :: (IArray a (u Int e), IArray u e)
-        => (Int -> Int -> (a Int (u Int e)) -> matrix e)
-        -> (Int, Int)
-        -> ((Int, Int) -> e)
-        -> matrix e
-_matrix_ c (numRows, numCols) generator =
-    c numRows numCols
-      $ array (1, numRows)
-      $ [ (i, array (1, numCols) [(j, generator (i, j))
-                                 | j <- [1..numCols]])
-          | i <- [1..numRows] ]
--}
-
 _matrix :: (IArray a1 (u Int e), IArray u e,
             MArray a2 (u Int e) (ST s), MArray a3 e (ST s),
             Num e)
@@ -813,38 +815,13 @@
         -> ST s matrix
 _matrix c newArray newArrayU (m, n) g = do
     rows <- newArray (1, m) undefined
-    flip mapM_ [1..m] $ \i -> do
+    forM_ [1..m] $ \i -> do
         cols <- newArrayU (1, n) 0
-        flip mapM_ [1..n] $ \j -> do
-            writeArray cols j (g (i,j))
-        U.unsafeFreeze cols >>= writeArray rows i
+        forM_ [1..n] $ \j -> do
+            unsafeWrite cols j (g (i,j))
+        U.unsafeFreeze cols >>= unsafeWrite rows i
     U.unsafeFreeze rows >>= return . c m n
 
-{-
-_matrixIO c newArray newArrayU (m, n) g = do
-    rows <- newArray (1, m) undefined
-
-    let half = m `div` 2
-
-    wait <- newEmptyMVar
-    forkOS $ do
-        flip mapM_ [1..half] $ \i -> do
-            cols <- newArrayU (1, n) 0
-            flip mapM_ [1..n] $ \j -> do
-                writeArray cols j (g (i,j))
-            U.unsafeFreeze cols >>= writeArray rows i
-        putMVar wait ()
-
-    flip mapM_ [half+1..m] $ \i -> do
-        cols <- newArrayU (1, n) 0
-        flip mapM_ [1..n] $ \j -> do
-            writeArray cols j (g (i,j))
-        U.unsafeFreeze cols >>= writeArray rows i
-
-    takeMVar wait
-
-    U.unsafeFreeze rows >>= return . c m n
--}
 
 {-# RULES