diff --git a/Changes.md b/Changes.md
new file mode 100644
--- /dev/null
+++ b/Changes.md
@@ -0,0 +1,20 @@
+# Change log for the `comfort-array` package
+
+## 0.2
+
+ * Add a monad parameter to the mutable `Storable` array type
+   and generalize functions to `PrimMonad`s.
+   This way the mutating functions can also be used in the `ST` monad.
+
+## 0.1.2
+
+ * Add immutable `Boxed` array type and mutable `Storable` array type.
+
+## 0.1
+
+ * Split `Shape.C` into `Shape.C` and `Shape.Indexed`.
+
+## 0.0
+
+ * Initial version featuring the `Shape.C` class with type function `Index`
+   and the immutable `Storable` array type.
diff --git a/comfort-array.cabal b/comfort-array.cabal
--- a/comfort-array.cabal
+++ b/comfort-array.cabal
@@ -1,5 +1,5 @@
 Name:             comfort-array
-Version:          0.1.2
+Version:          0.2
 License:          BSD3
 License-File:     LICENSE
 Author:           Henning Thielemann <haskell@henning-thielemann.de>
@@ -51,9 +51,11 @@
 Tested-With:      GHC==7.4.2, GHC==7.8.4, GHC==8.2.2
 Cabal-Version:    1.14
 Build-Type:       Simple
+Extra-Source-Files:
+  Changes.md
 
 Source-Repository this
-  Tag:         0.1.2
+  Tag:         0.2
   Type:        darcs
   Location:    http://hub.darcs.net/thielema/comfort-array/
 
@@ -66,8 +68,10 @@
     primitive >=0.6.4 && <0.7,
     guarded-allocation >=0.0 && <0.1,
     storable-record >=0.0.1 && <0.1,
+    deepseq >=1.3 && <1.5,
     QuickCheck >=2 && <3,
     transformers >=0.3 && <0.6,
+    non-empty >=0.3 && <0.4,
     utility-ht >=0.0.10 && <0.1,
     base >=4.5 && <5
 
@@ -79,10 +83,13 @@
     Data.Array.Comfort.Shape.Test
     Data.Array.Comfort.Storable
     Data.Array.Comfort.Storable.Internal
+    Data.Array.Comfort.Storable.Internal.Monadic
+    Data.Array.Comfort.Storable.Private
     Data.Array.Comfort.Storable.Mutable
+    Data.Array.Comfort.Storable.Mutable.Internal
+    Data.Array.Comfort.Storable.Mutable.Private
     Data.Array.Comfort.Boxed
   Other-Modules:
-    Data.Array.Comfort.Storable.Mutable.Internal
     Data.Array.Comfort.Boxed.Internal
 
 Test-Suite comfort-array-test
diff --git a/src/Data/Array/Comfort/Boxed.hs b/src/Data/Array/Comfort/Boxed.hs
--- a/src/Data/Array/Comfort/Boxed.hs
+++ b/src/Data/Array/Comfort/Boxed.hs
@@ -1,16 +1,82 @@
-{-# LANGUAGE TypeFamilies #-}
 module Data.Array.Comfort.Boxed (
-   Array.Array,
+   Array,
    shape,
-   (Array.!),
+   (!),
    Array.toList,
+   toAssociations,
    Array.fromList,
    Array.vectorFromList,
 
    Array.map,
+   zipWith,
+   (//),
+   accumulate,
+   fromAssociations,
    ) where
 
 import qualified Data.Array.Comfort.Boxed.Internal as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Boxed.Internal (Array(Array))
 
+import qualified Data.Primitive.Array as Prim
+
+import qualified Control.Monad.Primitive as PrimM
+import Control.Monad.ST (runST)
+import Control.Applicative ((<$>))
+
+import Data.Foldable (forM_)
+
+import Prelude hiding (zipWith)
+
+
 shape :: Array.Array sh a -> sh
 shape = Array.shape
+
+
+infixl 9 !
+
+(!) :: (Shape.Indexed sh) => Array sh a -> Shape.Index sh -> a
+(!) (Array sh arr) ix =
+   if Shape.inBounds sh ix
+      then Prim.indexArray arr $ Shape.offset sh ix
+      else error "Array.Comfort.Boxed.!: index out of bounds"
+
+
+zipWith ::
+   (Shape.C sh, Eq sh) =>
+   (a -> b -> c) -> Array sh a -> Array sh b -> Array sh c
+zipWith f a b =
+   if shape a == shape b
+      then Array.zipWith f a b
+      else error "zipWith: shapes mismatch"
+
+
+(//) ::
+   (Shape.Indexed sh) => Array sh a -> [(Shape.Index sh, a)] -> Array sh a
+(//) (Array sh arr) xs = runST (do
+   marr <- Prim.thawArray arr 0 (Shape.size sh)
+   forM_ xs $ \(ix,a) -> Prim.writeArray marr (Shape.offset sh ix) a
+   Array sh <$> Prim.unsafeFreezeArray marr)
+
+accumulate ::
+   (Shape.Indexed sh) =>
+   (a -> b -> a) -> Array sh a -> [(Shape.Index sh, b)] -> Array sh a
+accumulate f (Array sh arr) xs = runST (do
+   marr <- Prim.thawArray arr 0 (Shape.size sh)
+   forM_ xs $ \(ix,b) -> updateArray marr (Shape.offset sh ix) $ flip f b
+   Array sh <$> Prim.unsafeFreezeArray marr)
+
+updateArray ::
+   PrimM.PrimMonad m =>
+   Prim.MutableArray (PrimM.PrimState m) a -> Int -> (a -> a) -> m ()
+updateArray marr k f = Prim.writeArray marr k . f =<< Prim.readArray marr k
+
+toAssociations :: (Shape.Indexed sh) => Array sh a -> [(Shape.Index sh, a)]
+toAssociations arr = zip (Shape.indices $ shape arr) (Array.toList arr)
+
+fromAssociations ::
+   (Shape.Indexed sh) => sh -> a -> [(Shape.Index sh, a)] -> Array sh a
+fromAssociations sh a xs = runST (do
+   marr <- Prim.newArray (Shape.size sh) a
+   forM_ xs $ \(ix,x) -> Prim.writeArray marr (Shape.offset sh ix) x
+   Array sh <$> Prim.unsafeFreezeArray marr)
diff --git a/src/Data/Array/Comfort/Boxed/Internal.hs b/src/Data/Array/Comfort/Boxed/Internal.hs
--- a/src/Data/Array/Comfort/Boxed/Internal.hs
+++ b/src/Data/Array/Comfort/Boxed/Internal.hs
@@ -5,7 +5,15 @@
 import qualified Data.Primitive.Array as Prim
 
 import qualified Control.Monad.ST.Strict as STStrict
+import qualified Control.Monad.Trans.Class as MT
+import qualified Control.Monad.Trans.State as MS
+import Control.Monad (liftM)
+import Control.Applicative ((<$>))
+import Control.DeepSeq (NFData, rnf)
 
+import qualified Data.Traversable as Trav
+import qualified Data.Foldable as Fold
+import qualified Data.List as List
 import Prelude hiding (map, )
 
 
@@ -15,6 +23,33 @@
       buffer :: Prim.Array a
    }
 
+instance (Shape.C sh, Show sh, Show a) => Show (Array sh a) where
+   show arr =
+      "BoxedArray.fromList " ++
+      showsPrec 11 (shape arr) (' ' : show (toListLazy arr))
+
+
+instance (Shape.C sh, NFData sh, NFData a) => NFData (Array sh a) where
+   rnf a@(Array sh _arr) = rnf (sh, toListLazy a)
+
+instance (Shape.C sh) => Functor (Array sh) where
+   fmap = map
+
+instance (Shape.C sh) => Fold.Foldable (Array sh) where
+   fold = Fold.fold . buffer
+   foldMap f = Fold.foldMap f . buffer
+   foldl f a = Fold.foldl f a . buffer
+   foldr f a = Fold.foldr f a . buffer
+   foldl1 f = Fold.foldl1 f . buffer
+   foldr1 f = Fold.foldr1 f . buffer
+
+instance (Shape.C sh) => Trav.Traversable (Array sh) where
+   traverse f (Array sh arr) = Array sh <$> Trav.traverse f arr
+   sequenceA (Array sh arr) = Array sh <$> Trav.sequenceA arr
+   mapM f (Array sh arr) = liftM (Array sh) $ Trav.mapM f arr
+   sequence (Array sh arr) = liftM (Array sh) $ Trav.sequence arr
+
+
 -- add assertion, at least in an exposed version
 reshape :: sh1 -> Array sh0 a -> Array sh1 a
 reshape sh (Array _ arr) = Array sh arr
@@ -28,6 +63,10 @@
 (!) :: (Shape.Indexed sh) => Array sh a -> Shape.Index sh -> a
 (!) (Array sh arr) ix = Prim.indexArray arr $ Shape.offset sh ix
 
+toListLazy :: (Shape.C sh) => Array sh a -> [a]
+toListLazy (Array sh arr) =
+   List.map (Prim.indexArray arr) $ take (Shape.size sh) [0..]
+
 toList :: (Shape.C sh) => Array sh a -> [a]
 toList (Array sh arr) =
    STStrict.runST (mapM (Prim.indexArrayM arr) $ take (Shape.size sh) [0..])
@@ -42,3 +81,17 @@
 
 map :: (Shape.C sh) => (a -> b) -> Array sh a -> Array sh b
 map f (Array sh arr) = Array sh $ Prim.mapArray' f arr
+
+zipWith ::
+   (Shape.C sh) => (a -> b -> c) -> Array sh a -> Array sh b -> Array sh c
+zipWith f (Array sha arra) (Array _shb arrb) =
+   Array sha $
+   STStrict.runST
+      (flip MS.evalStateT 0 $
+       Prim.traverseArrayP
+         (\a -> do
+            k <- MS.get
+            b <- MT.lift $ Prim.indexArrayM arrb k
+            MS.put (k+1)
+            return $ f a b)
+         arra)
diff --git a/src/Data/Array/Comfort/Shape.hs b/src/Data/Array/Comfort/Shape.hs
--- a/src/Data/Array/Comfort/Shape.hs
+++ b/src/Data/Array/Comfort/Shape.hs
@@ -10,7 +10,11 @@
 
    Range(..),
    Shifted(..),
+   Enumeration(..),
+   Deferred(..), DeferredIndex, deferIndex, revealIndex,
+
    (:+:)(..),
+   Triangular(..), Lower(Lower), Upper(Upper), triangleSize, triangleRoot,
    ) where
 
 import qualified Foreign.Storable.Newtype as Store
@@ -23,10 +27,14 @@
 import qualified Control.Monad.Trans.State as MS
 import qualified Control.Monad.HT as Monad
 import qualified Control.Applicative.Backwards as Back
-import Control.Applicative (liftA2, liftA3)
+import Control.DeepSeq (NFData, rnf)
+import Control.Applicative (Applicative, pure, liftA2, liftA3, (<*>))
+import Control.Applicative (Const(Const, getConst))
 
 import Text.Printf (printf)
 
+import qualified Data.NonEmpty as NonEmpty
+import Data.List.HT (tails)
 import Data.Tuple.HT (mapFst, mapPair, swap)
 
 
@@ -96,6 +104,13 @@
 instance Functor ZeroBased where
    fmap f (ZeroBased n) = ZeroBased $ f n
 
+instance Applicative ZeroBased where
+   pure = ZeroBased
+   ZeroBased f <*> ZeroBased n = ZeroBased $ f n
+
+instance (NFData n) => NFData (ZeroBased n) where
+   rnf (ZeroBased n) = rnf n
+
 instance (Storable n) => Storable (ZeroBased n) where
    sizeOf = Store.sizeOf zeroBasedSize
    alignment = Store.alignment zeroBasedSize
@@ -131,6 +146,13 @@
 instance Functor OneBased where
    fmap f (OneBased n) = OneBased $ f n
 
+instance Applicative OneBased where
+   pure = OneBased
+   OneBased f <*> OneBased n = OneBased $ f n
+
+instance (NFData n) => NFData (OneBased n) where
+   rnf (OneBased n) = rnf n
+
 instance (Storable n) => Storable (OneBased n) where
    sizeOf = Store.sizeOf oneBasedSize
    alignment = Store.alignment oneBasedSize
@@ -169,6 +191,9 @@
 instance Functor Range where
    fmap f (Range from to) = Range (f from) (f to)
 
+instance (NFData n) => NFData (Range n) where
+   rnf (Range from to) = rnf (from,to)
+
 instance (Ix.Ix n) => C (Range n) where
    size (Range from to) = Ix.rangeSize (from,to)
    uncheckedSize (Range from to) = Ix.unsafeRangeSize (from,to)
@@ -213,6 +238,9 @@
 instance Functor Shifted where
    fmap f (Shifted from to) = Shifted (f from) (f to)
 
+instance (NFData n) => NFData (Shifted n) where
+   rnf (Shifted from to) = rnf (from,to)
+
 instance (Integral n) => C (Shifted n) where
    size (Shifted _offs len) = fromIntegral len
    uncheckedSize (Shifted _offs len) = fromIntegral len
@@ -267,6 +295,107 @@
 
 
 {- |
+'Enumeration' denotes a shape of fixed size
+that is defined by 'Enum' and 'Bounded' methods.
+For correctness it is necessary that the 'Enum' and 'Bounded'
+are properly implemented.
+Automatically derived instances are fine.
+-}
+data Enumeration n = Enumeration
+   deriving (Eq, Show)
+
+instance NFData (Enumeration n) where
+   rnf Enumeration = ()
+
+instance (Enum n, Bounded n) => C (Enumeration n) where
+   size = uncheckedSize
+   uncheckedSize sh = intFromEnum sh maxBound - intFromEnum sh minBound + 1
+
+instance (Enum n, Bounded n) => Indexed (Enumeration n) where
+   type Index (Enumeration n) = n
+   indices sh = [asEnumType sh minBound .. asEnumType sh maxBound]
+   offset = uncheckedOffset
+   uncheckedOffset sh ix = fromEnum ix - intFromEnum sh minBound
+   inBounds _sh _ix = True
+
+instance (Enum n, Bounded n) => InvIndexed (Enumeration n) where
+   indexFromOffset sh k =
+      if 0<=k && k <= intFromEnum sh maxBound - intFromEnum sh minBound
+         then uncheckedIndexFromOffset sh k
+         else errorIndexFromOffset "Enumeration" k
+   uncheckedIndexFromOffset sh k = toEnum $ intFromEnum sh minBound + k
+
+asEnumType :: Enumeration n -> n -> n
+asEnumType Enumeration = id
+
+intFromEnum :: (Enum n) => Enumeration n -> n -> Int
+intFromEnum Enumeration = fromEnum
+
+instance Storable (Enumeration n) where
+   {-# INLINE sizeOf #-}
+   {-# INLINE alignment #-}
+   {-# INLINE peek #-}
+   {-# INLINE poke #-}
+   sizeOf ~Enumeration = 0
+   alignment ~Enumeration = 1
+   poke _p Enumeration = return ()
+   peek _p = return Enumeration
+
+
+{- |
+This data type wraps another array shape.
+Its index type is a wrapped 'Int'.
+The advantages are:
+No conversion forth and back 'Int' and @Index sh@.
+You can convert once using 'deferIndex' and 'revealIndex'
+whenever you need your application specific index type.
+No need for e.g. @Storable (Index sh)@, because 'Int' is already 'Storable'.
+-}
+newtype Deferred sh = Deferred sh
+   deriving (Eq, Show)
+
+newtype DeferredIndex ix = DeferredIndex Int
+   deriving (Eq, Show)
+
+instance (NFData sh) => NFData (Deferred sh) where
+   rnf (Deferred sh) = rnf sh
+
+instance (C sh) => C (Deferred sh) where
+   size (Deferred sh) = size sh
+   uncheckedSize (Deferred sh) = uncheckedSize sh
+
+instance (C sh) => Indexed (Deferred sh) where
+   type Index (Deferred sh) = DeferredIndex (Index sh)
+   indices (Deferred sh) = map DeferredIndex $ take (size sh) [0 ..]
+   offset (Deferred sh) (DeferredIndex k) = offset (ZeroBased $ size sh) k
+   uncheckedOffset _ (DeferredIndex k) = k
+   inBounds (Deferred sh) (DeferredIndex k) =
+      inBounds (ZeroBased $ size sh) k
+
+instance (C sh) => InvIndexed (Deferred sh) where
+   indexFromOffset sh k =
+      DeferredIndex $ indexFromOffset (ZeroBased $ size sh) k
+   uncheckedIndexFromOffset _sh = DeferredIndex
+
+deferIndex :: (Indexed sh, Index sh ~ ix) => sh -> ix -> DeferredIndex ix
+deferIndex sh ix = DeferredIndex $ offset sh ix
+
+revealIndex :: (InvIndexed sh, Index sh ~ ix) => sh -> DeferredIndex ix -> ix
+revealIndex sh (DeferredIndex ix) = indexFromOffset sh ix
+
+instance Storable (DeferredIndex ix) where
+   {-# INLINE sizeOf #-}
+   {-# INLINE alignment #-}
+   {-# INLINE peek #-}
+   {-# INLINE poke #-}
+   sizeOf (DeferredIndex k) = sizeOf k
+   alignment (DeferredIndex k) = alignment k
+   poke p (DeferredIndex k) = poke (castPtr p) k
+   peek p = fmap DeferredIndex $ peek (castPtr p)
+
+
+
+{- |
 Row-major composition of two dimensions.
 -}
 instance (C sh0, C sh1) => C (sh0,sh1) where
@@ -368,10 +497,114 @@
 
 
 
+data Lower = Lower deriving (Eq, Show)
+data Upper = Upper deriving (Eq, Show)
+
+class TriangularPart part where
+   switchTriangularPart :: f Lower -> f Upper -> f part
+instance TriangularPart Lower where switchTriangularPart f _ = f
+instance TriangularPart Upper where switchTriangularPart _ f = f
+
+getConstAs :: c -> Const a c -> a
+getConstAs _ = getConst
+
+caseTriangularPart :: (TriangularPart part) => part -> a -> a -> a
+caseTriangularPart part lo up =
+   getConstAs part $ switchTriangularPart (Const lo) (Const up)
+
+data Triangular part size =
+   Triangular {
+      triangularPart :: part,
+      triangularSize :: size
+   } deriving (Eq, Show)
+
+-- cf. Data.Bifunctor.Flip
+newtype Flip f b a = Flip {getFlip :: f a b}
+
+instance
+      (TriangularPart part, NFData size) => NFData (Triangular part size) where
+   rnf (Triangular part sz) =
+      rnf
+         (flip getFlip part $
+            switchTriangularPart (Flip $ \Lower -> ()) (Flip $ \Upper -> ()),
+          sz)
+
+instance (TriangularPart part, C size) => C (Triangular part size) where
+   size (Triangular _part sz) = triangleSize $ size sz
+   uncheckedSize (Triangular _part sz) = triangleSize $ uncheckedSize sz
+
+instance
+   (TriangularPart part, Indexed size) =>
+      Indexed (Triangular part size) where
+   type Index (Triangular part size) = (Index size, Index size)
+
+   indices (Triangular part sz) =
+      let ixs = indices sz
+      in concat $
+         caseTriangularPart part
+            (zipWith (\cs r -> map ((,) r) cs)
+               (NonEmpty.tail $ NonEmpty.inits ixs) ixs)
+            (zipWith (\r cs -> map ((,) r) cs) ixs $ tails ixs)
+
+   uncheckedOffset sh ix = snd $ uncheckedSizeOffset sh ix
+
+   sizeOffset sh ix =
+      if inBounds sh ix
+        then uncheckedSizeOffset sh ix
+        else error "Shape.Triangular.sizeOffset: wrong array part"
+
+   uncheckedSizeOffset (Triangular part sz) (rs,cs) =
+      let (n,r) = uncheckedSizeOffset sz rs
+          c = uncheckedOffset sz cs
+      in (triangleSize n,
+          caseTriangularPart part
+            (triangleSize r + c)
+            (triangleSize n - triangleSize (n-r) + c-r))
+
+   inBounds (Triangular part sz) ix@(r,c) =
+      inBounds (sz,sz) ix
+      &&
+      caseTriangularPart part (>=) (<=) (offset sz r) (offset sz c)
+
+instance
+   (TriangularPart part, InvIndexed size) =>
+      InvIndexed (Triangular part size) where
+
+   indexFromOffset (Triangular part sz) k =
+      mapPair (indexFromOffset sz, indexFromOffset sz) $
+      caseTriangularPart part
+         (let r = floor (triangleRootDouble k)
+          in (r, k - triangleSize r))
+         (let n = size sz
+              triSize = triangleSize n
+              rr = ceiling (triangleRootDouble (triSize-k))
+              r = n - rr
+          in (r, k+r - (triSize - triangleSize rr)))
+
+triangleSize :: Int -> Int
+triangleSize n = div (n*(n+1)) 2
+
+{-
+n*(n+1)/2 = m
+n^2 + n - 2m = 0
+n = -1/2 + sqrt(1/4+2m)
+  = (sqrt(8m+1) - 1) / 2
+-}
+triangleRoot :: Floating a => a -> a
+triangleRoot sz = (sqrt (8*sz+1)-1)/2
+
+triangleRootDouble :: Int -> Double
+triangleRootDouble = triangleRoot . fromIntegral
+
+
+
 infixr 5 :+:
 
 data sh0:+:sh1 = sh0:+:sh1
    deriving (Eq, Show)
+
+instance (NFData sh0, NFData sh1) => NFData (sh0:+:sh1) where
+   rnf (sh0:+:sh1) = rnf (sh0,sh1)
 
 instance (C sh0, C sh1) => C (sh0:+:sh1) where
    size (sh0:+:sh1) = size sh0 + size sh1
diff --git a/src/Data/Array/Comfort/Storable.hs b/src/Data/Array/Comfort/Storable.hs
--- a/src/Data/Array/Comfort/Storable.hs
+++ b/src/Data/Array/Comfort/Storable.hs
@@ -1,17 +1,85 @@
-{-# LANGUAGE TypeFamilies #-}
 module Data.Array.Comfort.Storable (
-   Array.Array,
+   Array,
    shape,
-   (Array.!),
+   reshape,
+   mapShape,
+
+   (!),
    Array.toList,
    Array.fromList,
    Array.vectorFromList,
 
    Array.map,
-   (Array.//),
+   Array.mapWithIndex,
+   (//),
+   accumulate,
+   fromAssociations,
    ) where
 
+import qualified Data.Array.Comfort.Storable.Mutable.Internal as MutArray
 import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Internal (Array)
 
-shape :: Array.Array sh a -> sh
+import Foreign.Storable (Storable)
+
+import Control.Monad.ST (runST)
+
+import Data.Foldable (forM_)
+
+import Text.Printf (printf)
+
+import Prelude hiding (map)
+
+
+shape :: Array sh a -> sh
 shape = Array.shape
+
+reshape :: (Shape.C sh0, Shape.C sh1) => sh1 -> Array sh0 a -> Array sh1 a
+reshape sh1 arr =
+   let n0 = Shape.size $ shape arr
+       n1 = Shape.size sh1
+   in if n0 == n1
+         then Array.reshape sh1 arr
+         else error $
+              printf
+                 ("Array.Comfort.Storable.reshape: " ++
+                  "different sizes of old (%d) and new (%d) shape")
+                 n0 n1
+
+mapShape ::
+   (Shape.C sh0, Shape.C sh1) => (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
+mapShape f arr = reshape (f $ shape arr) arr
+
+
+infixl 9 !
+
+(!) :: (Shape.Indexed sh, Storable a) => Array sh a -> Shape.Index sh -> a
+(!) arr ix = runST (do
+   marr <- MutArray.unsafeThaw arr
+   MutArray.read marr ix)
+
+
+(//) ::
+   (Shape.Indexed sh, Storable a) =>
+   Array sh a -> [(Shape.Index sh, a)] -> Array sh a
+(//) arr xs = runST (do
+   marr <- MutArray.thaw arr
+   forM_ xs $ uncurry $ MutArray.write marr
+   MutArray.unsafeFreeze marr)
+
+accumulate ::
+   (Shape.Indexed sh, Storable a) =>
+   (a -> b -> a) -> Array sh a -> [(Shape.Index sh, b)] -> Array sh a
+accumulate f arr xs = runST (do
+   marr <- MutArray.thaw arr
+   forM_ xs $ \(ix,b) -> MutArray.update marr ix $ flip f b
+   MutArray.unsafeFreeze marr)
+
+fromAssociations ::
+   (Shape.Indexed sh, Storable a) =>
+   sh -> a -> [(Shape.Index sh, a)] -> Array sh a
+fromAssociations sh a xs = runST (do
+   marr <- MutArray.new sh a
+   forM_ xs $ uncurry $ MutArray.write marr
+   MutArray.unsafeFreeze marr)
diff --git a/src/Data/Array/Comfort/Storable/Internal.hs b/src/Data/Array/Comfort/Storable/Internal.hs
--- a/src/Data/Array/Comfort/Storable/Internal.hs
+++ b/src/Data/Array/Comfort/Storable/Internal.hs
@@ -1,90 +1,57 @@
 {-# LANGUAGE TypeFamilies #-}
+{- |
+The functions in this module miss any bound checking.
+-}
 module Data.Array.Comfort.Storable.Internal (
-   Array(Array, shape, buffer),
-   reshape,
-   mapShape,
+   Priv.Array(Array, shape, buffer),
+   Priv.reshape,
+   Priv.mapShape,
 
-   (!),
+   (Priv.!),
    unsafeCreate,
    unsafeCreateWithSize,
    unsafeCreateWithSizeAndResult,
-   toList,
-   fromList,
-   vectorFromList,
+   Priv.toList,
+   Priv.fromList,
+   Priv.vectorFromList,
 
    map,
-   copyIO,
-   (//),
-
-   createIO,
-   createWithSizeIO,
-   createWithSizeAndResultIO,
-   showIO,
-   readIO,
-   toListIO,
-   fromListIO,
-   vectorFromListIO,
+   mapWithIndex,
+   (Priv.//),
+   Priv.accumulate,
+   Priv.fromAssociations,
    ) where
 
+import qualified Data.Array.Comfort.Storable.Internal.Monadic as Monadic
+import qualified Data.Array.Comfort.Storable.Private as Priv
 import qualified Data.Array.Comfort.Shape as Shape
-
-import qualified Foreign.Marshal.Array.Guarded as Alloc
-import Foreign.Marshal.Array (copyArray, pokeArray, peekArray, advancePtr, )
-import Foreign.Storable (Storable, poke, pokeElemOff, peek, peekElemOff, )
-import Foreign.ForeignPtr (ForeignPtr, withForeignPtr, )
-import Foreign.Ptr (Ptr, )
-
-import System.IO.Unsafe (unsafePerformIO, )
-
-import Control.Applicative ((<$>))
-
-import Data.Tuple.HT (mapFst)
-
-import Prelude hiding (map, readIO, )
-
-
-data Array sh a =
-   Array {
-      shape :: sh,
-      buffer :: ForeignPtr a
-   }
+import Data.Array.Comfort.Storable.Private (Array(Array))
 
-instance (Shape.C sh, Show sh, Storable a, Show a) => Show (Array sh a) where
-   show = unsafePerformIO . showIO
+import Foreign.Marshal.Array (advancePtr)
+import Foreign.Storable (Storable, poke, peek)
+import Foreign.ForeignPtr (withForeignPtr)
+import Foreign.Ptr (Ptr)
 
-reshape :: sh1 -> Array sh0 a -> Array sh1 a
-reshape sh (Array _ fptr) = Array sh fptr
+import Control.Monad.ST (runST)
 
-mapShape :: (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
-mapShape f (Array sh fptr) = Array (f sh) fptr
+import Prelude hiding (map)
 
 
-infixl 9 !
-
 unsafeCreate ::
-   (Shape.C sh, Storable a) => sh -> (Ptr a -> IO ()) -> Array sh a
-unsafeCreate sh = unsafePerformIO . createIO sh
+   (Shape.C sh, Storable a) =>
+   sh -> (Ptr a -> IO ()) -> Array sh a
+unsafeCreate sh arr = runST (Monadic.unsafeCreate sh arr)
 
 unsafeCreateWithSize ::
-   (Shape.C sh, Storable a) => sh -> (Int -> Ptr a -> IO ()) -> Array sh a
-unsafeCreateWithSize sh = unsafePerformIO . createWithSizeIO sh
+   (Shape.C sh, Storable a) =>
+   sh -> (Int -> Ptr a -> IO ()) -> Array sh a
+unsafeCreateWithSize sh arr = runST (Monadic.unsafeCreateWithSize sh arr)
 
 unsafeCreateWithSizeAndResult ::
-   (Shape.C sh, Storable a) => sh -> (Int -> Ptr a -> IO b) -> (Array sh a, b)
-unsafeCreateWithSizeAndResult sh =
-   unsafePerformIO . createWithSizeAndResultIO sh
-
-(!) :: (Shape.Indexed sh, Storable a) => Array sh a -> Shape.Index sh -> a
-(!) arr = unsafePerformIO . readIO arr
-
-toList :: (Shape.C sh, Storable a) => Array sh a -> [a]
-toList = unsafePerformIO . toListIO
-
-fromList :: (Shape.C sh, Storable a) => sh -> [a] -> Array sh a
-fromList sh = unsafePerformIO . fromListIO sh
-
-vectorFromList :: (Storable a) => [a] -> Array (Shape.ZeroBased Int) a
-vectorFromList = unsafePerformIO . vectorFromListIO
+   (Shape.C sh, Storable a) =>
+   sh -> (Int -> Ptr a -> IO b) -> (Array sh a, b)
+unsafeCreateWithSizeAndResult sh arr =
+   runST (Monadic.unsafeCreateWithSizeAndResult sh arr)
 
 
 map ::
@@ -99,61 +66,15 @@
          (iterate (flip advancePtr 1) srcPtr)
          (iterate (flip advancePtr 1) dstPtr)
 
-copyIO :: (Shape.C sh, Storable a) => Array sh a -> IO (Array sh a)
-copyIO (Array sh srcFPtr) =
-   withForeignPtr srcFPtr $ \srcPtr ->
-   createWithSizeIO sh $ \n dstPtr ->
-      copyArray dstPtr srcPtr n
-
-(//) ::
-   (Shape.Indexed sh, Storable a) =>
-   Array sh a -> [(Shape.Index sh, a)] -> Array sh a
-(//) (Array sh fptr) xs =
-   unsafeCreateWithSize sh $ \n dstPtr ->
-   withForeignPtr fptr $ \srcPtr -> do
-      copyArray dstPtr srcPtr n
-      mapM_ (\(ix,a) -> pokeElemOff dstPtr (Shape.offset sh ix) a) xs
-
-
-
-createIO ::
-   (Shape.C sh, Storable a) => sh -> (Ptr a -> IO ()) -> IO (Array sh a)
-createIO sh f = createWithSizeIO sh $ const f
-
-createWithSizeIO ::
-   (Shape.C sh, Storable a) => sh -> (Int -> Ptr a -> IO ()) -> IO (Array sh a)
-createWithSizeIO sh f =
-   fst <$> createWithSizeAndResultIO sh f
-
-createWithSizeAndResultIO ::
-   (Shape.C sh, Storable a) =>
-   sh -> (Int -> Ptr a -> IO b) -> IO (Array sh a, b)
-createWithSizeAndResultIO sh f =
-   let size = Shape.size sh
-   in fmap (mapFst (Array sh)) $ Alloc.create size $ f size
-
-showIO :: (Shape.C sh, Show sh, Storable a, Show a) => Array sh a -> IO String
-showIO arr = do
-   xs <- toListIO arr
-   return $ "fromList " ++ showsPrec 11 (shape arr) (' ' : show xs)
-
-readIO :: (Shape.Indexed sh, Storable a) => Array sh a -> Shape.Index sh -> IO a
-readIO (Array sh fptr) ix =
-   withForeignPtr fptr $ flip peekElemOff (Shape.offset sh ix)
-
-toListIO :: (Shape.C sh, Storable a) => Array sh a -> IO [a]
-toListIO (Array sh fptr) =
-   withForeignPtr fptr $ peekArray (Shape.size sh)
-
-fromListIO ::
-   (Shape.C sh, Storable a) =>
-   sh -> [a] -> IO (Array sh a)
-fromListIO sh xs =
-   createWithSizeIO sh $ \size ptr ->
-      pokeArray ptr $ take size $
-      xs ++
-      repeat (error "Array.Comfort.Storable.fromList: list too short for shape")
-
-vectorFromListIO :: (Storable a) => [a] -> IO (Array (Shape.ZeroBased Int) a)
-vectorFromListIO xs =
-   createIO (Shape.ZeroBased $ length xs) $ flip pokeArray xs
+mapWithIndex ::
+   (Shape.Indexed sh, Shape.Index sh ~ ix, Storable a, Storable b) =>
+   (ix -> a -> b) -> Array sh a -> Array sh b
+mapWithIndex f (Array sh a) =
+   unsafeCreate sh $ \dstPtr ->
+   withForeignPtr a $ \srcPtr ->
+   sequence_ $
+      zipWith3
+         (\ix src dst -> poke dst . f ix =<< peek src)
+         (Shape.indices sh)
+         (iterate (flip advancePtr 1) srcPtr)
+         (iterate (flip advancePtr 1) dstPtr)
diff --git a/src/Data/Array/Comfort/Storable/Internal/Monadic.hs b/src/Data/Array/Comfort/Storable/Internal/Monadic.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Comfort/Storable/Internal/Monadic.hs
@@ -0,0 +1,30 @@
+module Data.Array.Comfort.Storable.Internal.Monadic where
+
+import qualified Data.Array.Comfort.Storable.Mutable.Private as MutArray
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Private (Array, unsafeFreeze)
+
+import Foreign.Storable (Storable, )
+import Foreign.Ptr (Ptr, )
+
+import Control.Monad.Primitive (PrimMonad)
+import Control.Monad.HT ((<=<))
+
+
+unsafeCreate ::
+   (PrimMonad m, Shape.C sh, Storable a) =>
+   sh -> (Ptr a -> IO ()) -> m (Array sh a)
+unsafeCreate sh = unsafeFreeze <=< MutArray.unsafeCreate sh
+
+unsafeCreateWithSize ::
+   (PrimMonad m, Shape.C sh, Storable a) =>
+   sh -> (Int -> Ptr a -> IO ()) -> m (Array sh a)
+unsafeCreateWithSize sh = unsafeFreeze <=< MutArray.unsafeCreateWithSize sh
+
+unsafeCreateWithSizeAndResult ::
+   (PrimMonad m, Shape.C sh, Storable a) =>
+   sh -> (Int -> Ptr a -> IO b) -> m (Array sh a, b)
+unsafeCreateWithSizeAndResult sh f = do
+   (arr, b) <- MutArray.unsafeCreateWithSizeAndResult sh f
+   marr <- unsafeFreeze arr
+   return (marr, b)
diff --git a/src/Data/Array/Comfort/Storable/Mutable.hs b/src/Data/Array/Comfort/Storable/Mutable.hs
--- a/src/Data/Array/Comfort/Storable/Mutable.hs
+++ b/src/Data/Array/Comfort/Storable/Mutable.hs
@@ -1,13 +1,67 @@
 module Data.Array.Comfort.Storable.Mutable (
-   Array.Array,
+   Array,
+   MutArray.STArray,
+   MutArray.IOArray,
    shape,
-   Array.read,
-   Array.write,
-   Array.thaw,
-   Array.freeze,
+
+   MutArray.new,
+   read,
+   write,
+   update,
+   toList,
+   fromList,
+   vectorFromList,
+
+   MutArray.thaw,
+   MutArray.freeze,
    ) where
 
-import qualified Data.Array.Comfort.Storable.Mutable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Mutable.Internal as MutArray
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Mutable.Internal (Array)
 
-shape :: Array.Array sh a -> sh
-shape = Array.shape
+import Foreign.Storable (Storable)
+
+import Control.Monad.Primitive (PrimMonad)
+
+import Prelude hiding (read)
+
+
+shape :: Array m sh a -> sh
+shape = MutArray.shape
+
+
+read ::
+   (PrimMonad m, Shape.Indexed sh, Storable a) =>
+   Array m sh a -> Shape.Index sh -> m a
+read arr ix =
+   if Shape.inBounds (shape arr) ix
+      then MutArray.read arr ix
+      else error "Array.Comfort.Storable.Mutable.read: index out of bounds"
+
+write ::
+   (PrimMonad m, Shape.Indexed sh, Storable a) =>
+   Array m sh a -> Shape.Index sh -> a -> m ()
+write arr ix a =
+   if Shape.inBounds (shape arr) ix
+      then MutArray.write arr ix a
+      else error "Array.Comfort.Storable.Mutable.write: index out of bounds"
+
+update ::
+   (PrimMonad m, Shape.Indexed sh, Storable a) =>
+   Array m sh a -> Shape.Index sh -> (a -> a) -> m ()
+update arr ix f =
+   if Shape.inBounds (shape arr) ix
+      then MutArray.update arr ix f
+      else error "Array.Comfort.Storable.Mutable.update: index out of bounds"
+
+toList :: (PrimMonad m, Shape.C sh, Storable a) => Array m sh a -> m [a]
+toList = MutArray.toList
+
+fromList ::
+   (PrimMonad m, Shape.C sh, Storable a) => sh -> [a] -> m (Array m sh a)
+fromList = MutArray.fromList
+
+vectorFromList ::
+   (PrimMonad m, Storable a) => [a] -> m (Array m (Shape.ZeroBased Int) a)
+vectorFromList = MutArray.vectorFromList
diff --git a/src/Data/Array/Comfort/Storable/Mutable/Internal.hs b/src/Data/Array/Comfort/Storable/Mutable/Internal.hs
--- a/src/Data/Array/Comfort/Storable/Mutable/Internal.hs
+++ b/src/Data/Array/Comfort/Storable/Mutable/Internal.hs
@@ -1,97 +1,28 @@
-{-# LANGUAGE TypeFamilies #-}
-module Data.Array.Comfort.Storable.Mutable.Internal where
-
-import qualified Data.Array.Comfort.Storable.Internal as Imm
-import qualified Data.Array.Comfort.Shape as Shape
-
-import qualified Foreign.Marshal.Array.Guarded as Alloc
-import Foreign.Marshal.Array (pokeArray, peekArray, )
-import Foreign.Storable (Storable, pokeElemOff, peekElemOff, )
-import Foreign.ForeignPtr (ForeignPtr, withForeignPtr, )
-import Foreign.Ptr (Ptr, )
-
-import System.IO.Unsafe (unsafePerformIO, )
-
-import Control.Applicative ((<$>))
-import Control.Monad ((<=<))
-
-import Data.Tuple.HT (mapFst)
-
-import Prelude hiding (map, read, )
-
-
-data Array sh a =
-   Array {
-      shape :: sh,
-      buffer :: ForeignPtr a
-   }
-
-instance (Shape.C sh, Show sh, Storable a, Show a) => Show (Array sh a) where
-   show = unsafePerformIO . showIO
-
-reshape :: sh1 -> Array sh0 a -> Array sh1 a
-reshape sh (Array _ fptr) = Array sh fptr
-
-mapShape :: (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
-mapShape f (Array sh fptr) = Array (f sh) fptr
-
-
-create ::
-   (Shape.C sh, Storable a) => sh -> (Ptr a -> IO ()) -> IO (Array sh a)
-create sh f = createWithSize sh $ const f
-
-createWithSize ::
-   (Shape.C sh, Storable a) => sh -> (Int -> Ptr a -> IO ()) -> IO (Array sh a)
-createWithSize sh f =
-   fst <$> createWithSizeAndResult sh f
-
-createWithSizeAndResult ::
-   (Shape.C sh, Storable a) =>
-   sh -> (Int -> Ptr a -> IO b) -> IO (Array sh a, b)
-createWithSizeAndResult sh f =
-   let size = Shape.size sh
-   in fmap (mapFst (Array sh)) $ Alloc.create size $ f size
-
-showIO :: (Shape.C sh, Show sh, Storable a, Show a) => Array sh a -> IO String
-showIO arr = do
-   xs <- toList arr
-   return $ "fromList " ++ showsPrec 11 (shape arr) (' ' : show xs)
-
-read :: (Shape.Indexed sh, Storable a) => Array sh a -> Shape.Index sh -> IO a
-read (Array sh fptr) ix =
-   withForeignPtr fptr $ flip peekElemOff (Shape.offset sh ix)
-
-write ::
-   (Shape.Indexed sh, Storable a) => Array sh a -> Shape.Index sh -> a -> IO ()
-write (Array sh fptr) ix a =
-   withForeignPtr fptr $ \ptr -> pokeElemOff ptr (Shape.offset sh ix) a
-
-toList :: (Shape.C sh, Storable a) => Array sh a -> IO [a]
-toList (Array sh fptr) =
-   withForeignPtr fptr $ peekArray (Shape.size sh)
-
-fromList :: (Shape.C sh, Storable a) => sh -> [a] -> IO (Array sh a)
-fromList sh xs =
-   createWithSize sh $ \size ptr ->
-      pokeArray ptr $ take size $
-      xs ++
-      repeat (error "Array.Comfort.Storable.fromList: list too short for shape")
-
-vectorFromList :: (Storable a) => [a] -> IO (Array (Shape.ZeroBased Int) a)
-vectorFromList xs =
-   create (Shape.ZeroBased $ length xs) $ flip pokeArray xs
-
-
-freeze :: (Shape.Indexed sh, Storable a) => Array sh a -> IO (Imm.Array sh a)
-freeze = Imm.copyIO <=< unsafeFreeze
-
-thaw :: (Shape.Indexed sh, Storable a) => Imm.Array sh a -> IO (Array sh a)
-thaw = unsafeThaw <=< Imm.copyIO
+{- |
+The functions in this module miss any bound checking.
+-}
+module Data.Array.Comfort.Storable.Mutable.Internal (
+   MutArray.Array(Array, shape, buffer),
+   MutArray.STArray,
+   MutArray.IOArray,
+   MutArray.new,
+   MutArray.copy,
+   MutArray.create,
+   MutArray.createWithSize,
+   MutArray.createWithSizeAndResult,
+   MutArray.unsafeCreate,
+   MutArray.unsafeCreateWithSize,
+   MutArray.unsafeCreateWithSizeAndResult,
+   MutArray.read,
+   MutArray.write,
+   MutArray.update,
+   MutArray.toList,
+   MutArray.fromList,
+   MutArray.vectorFromList,
 
-unsafeFreeze ::
-   (Shape.Indexed sh, Storable a) => Array sh a -> IO (Imm.Array sh a)
-unsafeFreeze (Array sh fptr) = return (Imm.Array sh fptr)
+   ImmArray.freeze, ImmArray.unsafeFreeze,
+   ImmArray.thaw, ImmArray.unsafeThaw,
+   ) where
 
-unsafeThaw ::
-   (Shape.Indexed sh, Storable a) => Imm.Array sh a -> IO (Array sh a)
-unsafeThaw (Imm.Array sh fptr) = return (Array sh fptr)
+import qualified Data.Array.Comfort.Storable.Private as ImmArray
+import qualified Data.Array.Comfort.Storable.Mutable.Private as MutArray
diff --git a/src/Data/Array/Comfort/Storable/Mutable/Private.hs b/src/Data/Array/Comfort/Storable/Mutable/Private.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Comfort/Storable/Mutable/Private.hs
@@ -0,0 +1,131 @@
+{-# LANGUAGE TypeFamilies #-}
+module Data.Array.Comfort.Storable.Mutable.Private where
+
+import qualified Data.Array.Comfort.Shape as Shape
+
+import qualified Foreign.Marshal.Array.Guarded as Alloc
+import Foreign.Marshal.Array (copyArray, pokeArray, peekArray)
+import Foreign.Storable (Storable, pokeElemOff, peekElemOff)
+import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)
+import Foreign.Ptr (Ptr)
+
+import Control.Monad.Primitive (PrimMonad, unsafeIOToPrim)
+import Control.Monad.ST (ST)
+import Control.Monad (liftM)
+import Control.Applicative ((<$>))
+
+import Data.Tuple.HT (mapFst)
+
+import qualified Prelude as P
+import Prelude hiding (read, show)
+
+
+data Array (m :: * -> *) sh a =
+   Array {
+      shape :: sh,
+      buffer :: ForeignPtr a
+   }
+
+type STArray s = Array (ST s)
+type IOArray = Array IO
+
+
+copy ::
+   (PrimMonad m, Shape.C sh, Storable a) =>
+   Array m sh a -> m (Array m sh a)
+copy (Array sh srcFPtr) =
+   unsafeCreateWithSize sh $ \n dstPtr ->
+   withForeignPtr srcFPtr $ \srcPtr ->
+      copyArray dstPtr srcPtr n
+
+
+create ::
+   (Shape.C sh, Storable a) =>
+   sh -> (Ptr a -> IO ()) -> IO (IOArray sh a)
+create sh f = createWithSize sh $ const f
+
+createWithSize ::
+   (Shape.C sh, Storable a) =>
+   sh -> (Int -> Ptr a -> IO ()) -> IO (IOArray sh a)
+createWithSize sh f =
+   fst <$> createWithSizeAndResult sh f
+
+createWithSizeAndResult ::
+   (Shape.C sh, Storable a) =>
+   sh -> (Int -> Ptr a -> IO b) -> IO (IOArray sh a, b)
+createWithSizeAndResult sh f =
+   let size = Shape.size sh
+   in fmap (mapFst (Array sh)) $ Alloc.create size $ f size
+
+
+unsafeCreate ::
+   (PrimMonad m, Shape.C sh, Storable a) =>
+   sh -> (Ptr a -> IO ()) -> m (Array m sh a)
+unsafeCreate sh f = unsafeCreateWithSize sh $ const f
+
+unsafeCreateWithSize ::
+   (PrimMonad m, Shape.C sh, Storable a) =>
+   sh -> (Int -> Ptr a -> IO ()) -> m (Array m sh a)
+unsafeCreateWithSize sh f =
+   liftM fst $ unsafeCreateWithSizeAndResult sh f
+
+unsafeCreateWithSizeAndResult ::
+   (PrimMonad m, Shape.C sh, Storable a) =>
+   sh -> (Int -> Ptr a -> IO b) -> m (Array m sh a, b)
+unsafeCreateWithSizeAndResult sh f =
+   unsafeIOToPrim $
+   fmap (mapFst unsafeArrayIOToPrim) $ createWithSizeAndResult sh f
+
+unsafeArrayIOToPrim :: (PrimMonad m) => IOArray sh a -> Array m sh a
+unsafeArrayIOToPrim (Array sh fptr) = Array sh fptr
+
+
+show ::
+   (PrimMonad m, Shape.C sh, Show sh, Storable a, Show a) =>
+   Array m sh a -> m String
+show arr = do
+   xs <- toList arr
+   return $
+      "StorableArray.fromList " ++ showsPrec 11 (shape arr) (' ' : P.show xs)
+
+read ::
+   (PrimMonad m, Shape.Indexed sh, Storable a) =>
+   Array m sh a -> Shape.Index sh -> m a
+read (Array sh fptr) ix =
+   unsafeIOToPrim $ withForeignPtr fptr $ flip peekElemOff (Shape.offset sh ix)
+
+write ::
+   (PrimMonad m, Shape.Indexed sh, Storable a) =>
+   Array m sh a -> Shape.Index sh -> a -> m ()
+write (Array sh fptr) ix a =
+   unsafeIOToPrim $
+   withForeignPtr fptr $ \ptr -> pokeElemOff ptr (Shape.offset sh ix) a
+
+update ::
+   (PrimMonad m, Shape.Indexed sh, Storable a) =>
+   Array m sh a -> Shape.Index sh -> (a -> a) -> m ()
+update (Array sh fptr) ix f =
+   unsafeIOToPrim $
+   let k = Shape.offset sh ix
+   in withForeignPtr fptr $ \ptr -> pokeElemOff ptr k . f =<< peekElemOff ptr k
+
+new :: (PrimMonad m, Shape.C sh, Storable a) => sh -> a -> m (Array m sh a)
+new sh x =
+   unsafeCreateWithSize sh $ \size ptr -> pokeArray ptr $ replicate size x
+
+toList :: (PrimMonad m, Shape.C sh, Storable a) => Array m sh a -> m [a]
+toList (Array sh fptr) =
+   unsafeIOToPrim $ withForeignPtr fptr $ peekArray (Shape.size sh)
+
+fromList ::
+   (PrimMonad m, Shape.C sh, Storable a) => sh -> [a] -> m (Array m sh a)
+fromList sh xs =
+   unsafeCreateWithSize sh $ \size ptr ->
+      pokeArray ptr $ take size $
+      xs ++
+      repeat (error "Array.Comfort.Storable.fromList: list too short for shape")
+
+vectorFromList ::
+   (PrimMonad m, Storable a) => [a] -> m (Array m (Shape.ZeroBased Int) a)
+vectorFromList xs =
+   unsafeCreate (Shape.ZeroBased $ length xs) $ flip pokeArray xs
diff --git a/src/Data/Array/Comfort/Storable/Private.hs b/src/Data/Array/Comfort/Storable/Private.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Comfort/Storable/Private.hs
@@ -0,0 +1,95 @@
+{-# LANGUAGE TypeFamilies #-}
+module Data.Array.Comfort.Storable.Private where
+
+import qualified Data.Array.Comfort.Storable.Mutable.Private as MutArray
+import qualified Data.Array.Comfort.Shape as Shape
+
+import Foreign.Storable (Storable, )
+import Foreign.ForeignPtr (ForeignPtr, )
+
+import Control.DeepSeq (NFData, rnf)
+import Control.Monad.Primitive (PrimMonad)
+import Control.Monad.ST (runST)
+import Control.Monad.HT ((<=<))
+
+import Data.Foldable (forM_)
+
+
+data Array sh a =
+   Array {
+      shape :: sh,
+      buffer :: ForeignPtr a
+   }
+
+instance (Shape.C sh, Show sh, Storable a, Show a) => Show (Array sh a) where
+   show arr = runST (MutArray.show =<< unsafeThaw arr)
+
+instance (NFData sh) => NFData (Array sh a) where
+   rnf (Array sh fptr) = seq fptr (rnf sh)
+
+reshape :: sh1 -> Array sh0 a -> Array sh1 a
+reshape sh (Array _ fptr) = Array sh fptr
+
+mapShape :: (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
+mapShape f arr = reshape (f $ shape arr) arr
+
+
+infixl 9 !
+
+(!) :: (Shape.Indexed sh, Storable a) => Array sh a -> Shape.Index sh -> a
+(!) arr ix = runST (flip MutArray.read ix =<< unsafeThaw arr)
+
+toList :: (Shape.C sh, Storable a) => Array sh a -> [a]
+toList arr = runST (MutArray.toList =<< unsafeThaw arr)
+
+fromList :: (Shape.C sh, Storable a) => sh -> [a] -> Array sh a
+fromList sh arr = runST (unsafeFreeze =<< MutArray.fromList sh arr)
+
+vectorFromList :: (Storable a) => [a] -> Array (Shape.ZeroBased Int) a
+vectorFromList arr = runST (unsafeFreeze =<< MutArray.vectorFromList arr)
+
+
+(//) ::
+   (Shape.Indexed sh, Storable a) =>
+   Array sh a -> [(Shape.Index sh, a)] -> Array sh a
+(//) arr xs = runST (do
+   marr <- thaw arr
+   forM_ xs $ uncurry $ MutArray.write marr
+   unsafeFreeze marr)
+
+accumulate ::
+   (Shape.Indexed sh, Storable a) =>
+   (a -> b -> a) -> Array sh a -> [(Shape.Index sh, b)] -> Array sh a
+accumulate f arr xs = runST (do
+   marr <- thaw arr
+   forM_ xs $ \(ix,b) -> MutArray.update marr ix $ flip f b
+   unsafeFreeze marr)
+
+fromAssociations ::
+   (Shape.Indexed sh, Storable a) =>
+   sh -> a -> [(Shape.Index sh, a)] -> Array sh a
+fromAssociations sh a xs = runST (do
+   marr <- MutArray.new sh a
+   forM_ xs $ uncurry $ MutArray.write marr
+   unsafeFreeze marr)
+
+
+freeze ::
+   (PrimMonad m, Shape.C sh, Storable a) =>
+   MutArray.Array m sh a -> m (Array sh a)
+freeze = unsafeFreeze <=< MutArray.copy
+
+thaw ::
+   (PrimMonad m, Shape.C sh, Storable a) =>
+   Array sh a -> m (MutArray.Array m sh a)
+thaw = MutArray.copy <=< unsafeThaw
+
+unsafeFreeze ::
+   (PrimMonad m, Shape.C sh, Storable a) =>
+   MutArray.Array m sh a -> m (Array sh a)
+unsafeFreeze (MutArray.Array sh fptr) = return (Array sh fptr)
+
+unsafeThaw ::
+   (PrimMonad m, Shape.C sh, Storable a) =>
+   Array sh a -> m (MutArray.Array m sh a)
+unsafeThaw (Array sh fptr) = return (MutArray.Array sh fptr)
diff --git a/test/Test/Shape.hs b/test/Test/Shape.hs
--- a/test/Test/Shape.hs
+++ b/test/Test/Shape.hs
@@ -28,6 +28,16 @@
       (ShapeTest.tests $
        liftA2 Shape.Shifted
          (QC.choose (-10,10::Int)) (QC.choose (0,10::Int))) ++
+   prefix "Enumeration Bool"
+      (ShapeTest.tests $
+       return (Shape.Enumeration :: Shape.Enumeration Bool)) ++
+   prefix "Enumeration Ordering"
+      (ShapeTest.tests $
+       return (Shape.Enumeration :: Shape.Enumeration Ordering)) ++
+   prefix "Deferred Shifted"
+      (ShapeTest.tests $ fmap Shape.Deferred $
+       liftA2 Shape.Shifted
+         (QC.choose (-10,10::Int)) (QC.choose (0,10::Int))) ++
    prefix "Pair"
       (ShapeTest.tests $ liftA2 (,) (genZeroBased 10) (genZeroBased 10)) ++
    prefix "Triple"
@@ -35,4 +45,10 @@
        liftA3 (,,) (genZeroBased 10) (genZeroBased 10) (genZeroBased 10)) ++
    prefix "Append"
       (ShapeTest.tests $ liftA2 (:+:) (genZeroBased 10) (genZeroBased 10)) ++
+   prefix "Triangular Lower"
+      (ShapeTest.tests $
+       fmap (Shape.Triangular Shape.Lower) (genZeroBased 10)) ++
+   prefix "Triangular Upper"
+      (ShapeTest.tests $
+       fmap (Shape.Triangular Shape.Upper) (genZeroBased 10)) ++
    []