diff --git a/Changes.md b/Changes.md
--- a/Changes.md
+++ b/Changes.md
@@ -1,5 +1,23 @@
 # Change log for the `comfort-array` package
 
+## 0.5
+
+ * `Array.Boxed`.`map`, `zipWith`, `toList`: make lazy
+
+ * add `unified` methods to `Shape` classes:
+   `unifiedSize`, `unifiedOffset`, `unifiedSizeOffset`,
+   `uncheckedIndexFromOffset`.
+   They simplify to share code between checked and unchecked variants.
+   Actually, many implementations of these methods
+   recursively call themselves on part shapes.
+   However, the default methods have changed.
+
+ * `Shape.:+:` -> `Shape.::+`.
+   This resolves the name clash with the `:+:` operator from `tfp`.
+   It also highlights the right associativity and non-commutativity.
+
+ * `Shape.Simplex`
+
 ## 0.4.1
 
  * use `doctest-extract` for tests
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.4.1
+Version:          0.5
 License:          BSD3
 License-File:     LICENSE
 Author:           Henning Thielemann <haskell@henning-thielemann.de>
@@ -45,10 +45,10 @@
   * @Enumeration@:
     Arrays with indices like 'LT', 'EQ', 'GT' and a shape of fixed size.
   .
-  * @(:+:)@:
+  * @(::+)@:
     The Append type constructor allows to respresent block arrays,
     e.g. block matrices.
-    It also allows to represent non-empty arrays via @():+:sh@.
+    It also allows to represent non-empty arrays via @()::+sh@.
   .
   * @Set@: Use an arbitrary ordered set as index set.
   .
@@ -57,16 +57,32 @@
   * @Triangular@:
     A 2D array with the shape of a lower or upper triangular matrix.
   .
+  * @Simplex@:
+    Simplices of any dimension, where the dimension is encoded in the type.
+    An index is a tuple of monotonic ordered sub-indices.
+  .
   * @Square@: A 2D array where both dimensions always have equal size.
   .
   * @Cube@: A 3D array where all three dimensions always have equal size.
   .
   * @Tagged@: Statically distinguish shapes and indices that are isomorphic.
   .
-  The @lapack@ package defines even more fancy shapes
-  like tall rectangular matrices, triangular matrices and banded matrices.
+  With our @Array@ type you can perform
+  .
+  * Fast Linear Algebra using the package @lapack@.
+    The @lapack@ package defines even more fancy shapes
+    like tall rectangular matrices, triangular matrices and banded matrices.
+  .
+  * Fast Fourier Transforms using the package @comfort-fftw@
+  .
+  * Efficient Array Processing via LLVM Just-In-Time code generation
+    using the package @knead@.
+  .
+  See also @comfort-graph@ for a Graph data structure,
+  with non-Int node identifiers and flexible edge types.
 
-Tested-With:      GHC==7.4.2, GHC==7.8.4, GHC==8.2.2
+Tested-With:      GHC==7.4.2, GHC==7.8.4
+Tested-With:      GHC==8.2.2, GHC==8.6.5, GHC==8.10.4
 Cabal-Version:    1.14
 Build-Type:       Simple
 Extra-Source-Files:
@@ -74,7 +90,7 @@
   test-module.list
 
 Source-Repository this
-  Tag:         0.4.1
+  Tag:         0.5
   Type:        darcs
   Location:    https://hub.darcs.net/thielema/comfort-array/
 
@@ -88,6 +104,7 @@
 
 Library
   Build-Depends:
+    storablevector >=0.2 && <0.3,
     primitive >=0.6.4 && <0.8,
     guarded-allocation >=0.0.1 && <0.1,
     storable-record >=0.0.1 && <0.1,
@@ -122,11 +139,12 @@
     Data.Array.Comfort.Storable.Mutable.Unchecked
     Data.Array.Comfort.Storable.Mutable.Private
     Data.Array.Comfort.Boxed
+    Data.Array.Comfort.Boxed.Unchecked
     Data.Array.Comfort.Container
   Other-Modules:
     Data.Array.Comfort.Shape.Set
     Data.Array.Comfort.Shape.Utility
-    Data.Array.Comfort.Boxed.Unchecked
+    Data.Array.Comfort.Boxed.Strict.Unchecked
     Data.Array.Comfort.Storable.Memory
     Data.Array.Comfort.Check
 
diff --git a/set/0.4.0/Data/Array/Comfort/Shape/Set.hs b/set/0.4.0/Data/Array/Comfort/Shape/Set.hs
--- a/set/0.4.0/Data/Array/Comfort/Shape/Set.hs
+++ b/set/0.4.0/Data/Array/Comfort/Shape/Set.hs
@@ -5,24 +5,20 @@
 import qualified Data.Set as Set
 import Data.Set (Set)
 import Data.Tuple.HT (fst3)
+import Data.Maybe.HT (toMaybe)
 
 
-offset :: Ord a => Set a -> a -> Int
+offset :: Ord a => Set a -> a -> Maybe Int
 offset set ix =
    case Set.splitMember ix set of
-      (less, hit, _) ->
-         if hit
-            then Set.size less
-            else error "Shape.Set: array index not member of the index set"
+      (less, hit, _) -> toMaybe hit (Set.size less)
 
 uncheckedOffset :: Ord a => Set a -> a -> Int
 uncheckedOffset set = Set.size . fst3 . flip Set.splitMember set
 
-indexFromOffset :: Set a -> Int -> a
-indexFromOffset set k =
-   if 0<=k
-      then uncheckedIndexFromOffset set k
-      else errorIndexFromOffset "Set" k
+
+indexFromOffset :: Set a -> Int -> Maybe a
+indexFromOffset set k = toMaybe (0<=k) (uncheckedIndexFromOffset set k)
 
 uncheckedIndexFromOffset :: Set a -> Int -> a
 uncheckedIndexFromOffset set k =
diff --git a/set/0.5.4/Data/Array/Comfort/Shape/Set.hs b/set/0.5.4/Data/Array/Comfort/Shape/Set.hs
--- a/set/0.5.4/Data/Array/Comfort/Shape/Set.hs
+++ b/set/0.5.4/Data/Array/Comfort/Shape/Set.hs
@@ -2,12 +2,19 @@
 
 import qualified Data.Set as Set
 import Data.Set (Set)
+import Data.Maybe.HT (toMaybe)
 
 
-offset, uncheckedOffset :: Ord a => Set a -> a -> Int
-offset = flip Set.findIndex
-uncheckedOffset = offset
+offset :: Ord a => Set a -> a -> Maybe Int
+offset = flip Set.lookupIndex
 
-indexFromOffset, uncheckedIndexFromOffset :: Set a -> Int -> a
-indexFromOffset = flip Set.elemAt
-uncheckedIndexFromOffset = indexFromOffset
+uncheckedOffset :: Ord a => Set a -> a -> Int
+uncheckedOffset = flip Set.findIndex
+
+
+indexFromOffset :: Set a -> Int -> Maybe a
+indexFromOffset set k =
+   toMaybe (0<=k && k<Set.size set) (Set.elemAt k set)
+
+uncheckedIndexFromOffset :: Set a -> Int -> a
+uncheckedIndexFromOffset = flip Set.elemAt
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
@@ -3,7 +3,7 @@
    shape,
    reshape,
    mapShape,
-   (!),
+   accessMaybe, (!),
    Array.toList,
    Array.fromList,
    Array.vectorFromList,
@@ -40,6 +40,8 @@
 import Data.Map (Map)
 import Data.Set (Set)
 import Data.Foldable (forM_)
+import Data.Either.HT (maybeRight)
+import Data.Maybe (fromMaybe)
 
 import Prelude hiding (zipWith, replicate)
 
@@ -74,10 +76,14 @@
 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"
+(!) arr ix =
+   fromMaybe (error "Array.Comfort.Boxed.!: index out of bounds") $
+   accessMaybe arr ix
+
+accessMaybe :: (Shape.Indexed sh) => Array sh a -> Shape.Index sh -> Maybe a
+accessMaybe (Array sh arr) ix =
+   fmap (Prim.indexArray arr) $ maybeRight $
+   Shape.getChecked $ Shape.unifiedOffset sh ix
 
 
 zipWith ::
diff --git a/src/Data/Array/Comfort/Boxed/Strict/Unchecked.hs b/src/Data/Array/Comfort/Boxed/Strict/Unchecked.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Comfort/Boxed/Strict/Unchecked.hs
@@ -0,0 +1,35 @@
+{-# LANGUAGE TypeFamilies #-}
+module Data.Array.Comfort.Boxed.Strict.Unchecked where
+
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Boxed.Unchecked (Array(Array))
+
+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 Prelude hiding (map, zipWith)
+
+
+toList :: (Shape.C sh) => Array sh a -> [a]
+toList (Array sh arr) =
+   STStrict.runST (mapM (Prim.indexArrayM arr) $ take (Shape.size sh) [0..])
+
+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/Boxed/Unchecked.hs b/src/Data/Array/Comfort/Boxed/Unchecked.hs
--- a/src/Data/Array/Comfort/Boxed/Unchecked.hs
+++ b/src/Data/Array/Comfort/Boxed/Unchecked.hs
@@ -1,12 +1,23 @@
 {-# LANGUAGE TypeFamilies #-}
-module Data.Array.Comfort.Boxed.Unchecked where
+module Data.Array.Comfort.Boxed.Unchecked (
+   Array(..),
+   reshape,
+   mapShape,
+   (!),
+   toList,
+   fromList,
+   vectorFromList,
+   replicate,
+   map,
+   zipWith,
+   ) where
 
 import qualified Data.Array.Comfort.Shape as Shape
 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
+-- FixMe: In GHC-7.4.2 there is no instance PrimMonad (Lazy.ST s)
+-- import qualified Control.Monad.ST.Lazy as ST
+import qualified Control.Monad.ST.Strict as ST
 import Control.Monad (liftM)
 import Control.Applicative (Applicative, pure, (<*>), (<$>))
 import Control.DeepSeq (NFData, rnf)
@@ -29,11 +40,11 @@
          showString "BoxedArray.fromList " .
          showsPrec 11 (shape arr) .
          showChar ' ' .
-         shows (toListLazy arr)
+         shows (toList arr)
 
 
 instance (Shape.C sh, NFData sh, NFData a) => NFData (Array sh a) where
-   rnf a@(Array sh _arr) = rnf (sh, toListLazy a)
+   rnf a@(Array sh _arr) = rnf (sh, toList a)
 
 instance (Shape.C sh) => Functor (Array sh) where
    fmap = map
@@ -61,7 +72,6 @@
    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
 
@@ -74,13 +84,9 @@
 (!) :: (Shape.Indexed sh) => Array sh a -> Shape.Index sh -> a
 (!) (Array sh arr) ix = Prim.indexArray arr $ Shape.uncheckedOffset 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..])
+   List.map (Prim.indexArray arr) $ take (Shape.size sh) [0..]
 
 fromList :: (Shape.C sh) => sh -> [a] -> Array sh a
 fromList sh xs = Array sh $ Prim.fromListN (Shape.size sh) xs
@@ -93,21 +99,17 @@
 replicate :: (Shape.C sh) => sh -> a -> Array sh a
 replicate sh a =
    Array sh $
-   STStrict.runST (Prim.unsafeFreezeArray  =<< Prim.newArray (Shape.size sh) a)
+   ST.runST (Prim.unsafeFreezeArray  =<< Prim.newArray (Shape.size sh) a)
 
 map :: (Shape.C sh) => (a -> b) -> Array sh a -> Array sh b
-map f (Array sh arr) = Array sh $ Prim.mapArray' f arr
+map f (Array sh arr) = Array sh $
+   let n = Shape.size sh
+   in Prim.fromListN n $ List.map (f . Prim.indexArray arr) $ take n [0..]
 
 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)
+zipWith f (Array sha arra) (Array _shb arrb) = Array sha $
+   let n = Shape.size sha
+   in Prim.fromListN n $
+      List.map (\k -> f (Prim.indexArray arra k) (Prim.indexArray arrb k)) $
+      take n [0..]
diff --git a/src/Data/Array/Comfort/Container.hs b/src/Data/Array/Comfort/Container.hs
--- a/src/Data/Array/Comfort/Container.hs
+++ b/src/Data/Array/Comfort/Container.hs
@@ -49,7 +49,6 @@
 
 instance (C f) => Shape.C (Shape f) where
    size = shapeSize
-   uncheckedSize = shapeSize
 
 
 instance C [] where
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
@@ -1,11 +1,22 @@
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE EmptyDataDecls #-}
+{-# LANGUAGE GADTs #-}
 module Data.Array.Comfort.Shape (
    C(..),
    Indexed(..),
-   InvIndexed(..),
+   InvIndexed(..), messageIndexFromOffset, assertIndexFromOffset,
    Static(..),
 
+   requireCheck,
+   CheckSingleton(..),
+   Checking(..),
+   Result(..),
+   runChecked,
+   runUnchecked,
+   assert,
+   throwOrError,
+
    Zero(Zero),
    ZeroBased(..), zeroBasedSplit,
    OneBased(..),
@@ -15,7 +26,7 @@
    Enumeration(..),
    Deferred(..), DeferredIndex(..), deferIndex, revealIndex,
 
-   (:+:)(..),
+   (::+)(..),
 
    Square(..),
    Cube(..),
@@ -25,11 +36,23 @@
    lowerTriangular, upperTriangular,
    triangleSize, triangleRoot,
 
+   Simplex(..),
+   SimplexAscending, simplexAscending,
+   SimplexDescending, simplexDescending,
+   Ascending,
+   Descending,
+   SimplexOrder(..),
+   SimplexOrderC,
+   AllDistinct(..),
+   SomeRepetitive(..),
+   Collision(..),
+   CollisionC,
+
    Cyclic(..),
    ) where
 
 import qualified Data.Array.Comfort.Shape.Set as ShapeSet
-import Data.Array.Comfort.Shape.Utility (errorIndexFromOffset)
+import Data.Array.Comfort.Shape.Utility (messageIndexFromOffset, isRight)
 
 import qualified Foreign.Storable.Newtype as Store
 import Foreign.Storable
@@ -40,25 +63,32 @@
 
 import qualified Control.Monad.Trans.State as MS
 import qualified Control.Monad.HT as Monad
+import qualified Control.Applicative.HT as AppHT
 import qualified Control.Applicative.Backwards as Back
 import Control.DeepSeq (NFData, rnf)
-import Control.Applicative (Applicative, pure, liftA2, liftA3, (<*>))
+import Control.Monad (liftM)
+import Control.Applicative (Applicative, pure, liftA2, liftA3, (<*>), (<$>))
 import Control.Applicative (Const(Const, getConst))
+import Control.Functor.HT (void)
 
+import qualified Data.Functor.Classes as FunctorC
 import qualified Data.Traversable as Trav
 import qualified Data.Foldable as Fold
 import qualified Data.Map as Map
 import qualified Data.Set as Set
 import qualified Data.NonEmpty as NonEmpty
+import qualified Data.List.Match as Match
+import qualified Data.List.HT as ListHT
 import qualified Data.List as List
+import Data.Traversable (Traversable)
+import Data.Foldable (Foldable)
 import Data.Functor.Identity (Identity(Identity), runIdentity)
 import Data.Function.HT (compose2)
 import Data.Tagged (Tagged(Tagged, unTagged))
 import Data.Map (Map)
 import Data.Set (Set)
 import Data.List.HT (tails)
-import Data.Maybe (fromMaybe)
-import Data.Tuple.HT (mapSnd, mapPair, swap, fst3, snd3, thd3)
+import Data.Tuple.HT (mapFst, mapSnd, swap, fst3, snd3, thd3)
 import Data.Eq.HT (equating)
 
 
@@ -66,45 +96,153 @@
 >>> import qualified Data.Array.Comfort.Shape as Shape
 >>> import qualified Data.Map as Map
 >>> import qualified Data.Set as Set
->>> import Data.Array.Comfort.Shape ((:+:)((:+:)))
+>>> import Data.Array.Comfort.Shape ((::+)((::+)))
 -}
 
 
+data Checked
+data Unchecked
+
+class Checking check where
+   data Result check a
+   switchCheck :: f Checked -> f Unchecked -> f check
+
+data CheckSingleton check where
+   Checked :: CheckSingleton Checked
+   Unchecked :: CheckSingleton Unchecked
+
+autoCheck :: (Checking check) => CheckSingleton check
+autoCheck = switchCheck Checked Unchecked
+
+checkFromResult :: (Checking check) => Result check a -> CheckSingleton check
+checkFromResult _ = autoCheck
+
+withCheck ::
+   (Checking check) =>
+   (CheckSingleton check -> Result check a) -> Result check a
+withCheck f = f autoCheck
+
+requireCheck :: CheckSingleton check -> Result check a -> Result check a
+requireCheck _ = id
+
+
+instance Checking Checked where
+   newtype Result Checked a = CheckedResult {getChecked :: Either String a}
+   switchCheck f _ = f
+
+runChecked :: String -> Result Checked a -> a
+runChecked name (CheckedResult m) =
+   either (error . (("Shape." ++ name ++ ": ") ++)) id m
+
+instance Checking Unchecked where
+   newtype Result Unchecked a = UncheckedResult {getUnchecked :: a}
+   switchCheck _ f = f
+
+runUnchecked :: Result Unchecked a -> a
+runUnchecked = getUnchecked
+
+
+throw :: String -> Result Checked a
+throw = CheckedResult . Left
+
+throwOrError :: (Checking check) => String -> Result check a
+throwOrError msg = withCheck $ \check ->
+   case check of
+      Checked -> throw msg
+      Unchecked -> error msg
+
+assert :: (Checking check) => String -> Bool -> Result check ()
+assert msg cond = withCheck $ \check ->
+   case check of
+      Unchecked -> UncheckedResult ()
+      Checked -> if cond then pure () else throw msg
+
+
+instance (Checking check, Eq a) => Eq (Result check a) where
+   a0 == b0 =
+      case (checkFromResult a0, a0, b0) of
+         (Checked, CheckedResult a, CheckedResult b)  ->  a==b
+         (Unchecked, UncheckedResult a, UncheckedResult b)  ->  a==b
+
+instance (Checking check) => Functor (Result check) where
+   fmap f m =
+      case (checkFromResult m, m) of
+         (Checked, CheckedResult e) -> CheckedResult $ fmap f e
+         (Unchecked, UncheckedResult a) -> UncheckedResult $ f a
+
+instance (Checking check) => Applicative (Result check) where
+   pure a = withCheck $ \check ->
+      case check of
+         Checked -> CheckedResult $ Right a
+         Unchecked -> UncheckedResult a
+   f<*>a =
+      case (checkFromResult a, f, a) of
+         (Checked, CheckedResult ff, CheckedResult fa) ->
+            CheckedResult $ ff<*>fa
+         (Unchecked, UncheckedResult xf, UncheckedResult xa) ->
+            UncheckedResult $ xf xa
+
+instance (Checking check) => Monad (Result check) where
+   return = pure
+   a >>= b =
+      case (checkFromResult a, a) of
+         (Checked, CheckedResult e) -> CheckedResult $ getChecked . b =<< e
+         (Unchecked, UncheckedResult x) -> b x
+
+
 class C sh where
-   -- Ix.rangeSize
+   {-
+   This is the counterpart to 'Ix.rangeSize'.
+   We do not support a counterpart to 'Ix.unsafeRangeSize' anymore.
+   First, there is hardly any speed advantage
+   of using 'Ix.unsafeRangeSize' instead of 'Ix.rangeSize'.
+   Second, I do not know of an 'Ix' instance
+   where 'Ix.rangeSize' and 'Ix.unsafeRangeSize' differ.
+   -}
    size :: sh -> Int
-   -- Ix.unsafeRangeSize
-   uncheckedSize :: sh -> Int
-   uncheckedSize = size
 
 class C sh => Indexed sh where
-   {-# MINIMAL indices, (sizeOffset|offset), inBounds #-}
-   type Index sh :: *
+   {-# MINIMAL indices, (unifiedOffset|unifiedSizeOffset) #-}
+   type Index sh
    -- Ix.range
    indices :: sh -> [Index sh]
    -- Ix.index
    offset :: sh -> Index sh -> Int
-   offset sh = snd $ sizeOffset sh
+   offset sh = runChecked "offset" . unifiedOffset sh
    -- Ix.unsafeIndex
    uncheckedOffset :: sh -> Index sh -> Int
-   uncheckedOffset = offset
+   uncheckedOffset sh = getUnchecked . unifiedOffset sh
+   unifiedOffset :: (Checking check) => sh -> Index sh -> Result check Int
+   unifiedOffset sh = snd $ unifiedSizeOffset sh
    -- Ix.inRange
    inBounds :: sh -> Index sh -> Bool
+   inBounds sh = isRight . getChecked . unifiedOffset sh
 
    sizeOffset :: sh -> (Int, Index sh -> Int)
    sizeOffset sh = (size sh, offset sh)
    uncheckedSizeOffset :: sh -> (Int, Index sh -> Int)
-   uncheckedSizeOffset sh = (uncheckedSize sh, uncheckedOffset sh)
+   uncheckedSizeOffset sh = (size sh, uncheckedOffset sh)
+   unifiedSizeOffset ::
+      (Checking check) => sh -> (Int, Index sh -> Result check Int)
+   unifiedSizeOffset sh = (size sh, unifiedOffset sh)
 
 class Indexed sh => InvIndexed sh where
+   {-# MINIMAL unifiedIndexFromOffset #-}
    {- |
    It should hold @indexFromOffset sh k == indices sh !! k@,
    but 'indexFromOffset' should generally be faster.
    -}
    indexFromOffset :: sh -> Int -> Index sh
+   indexFromOffset sh = runChecked "indexFromOffset" . unifiedIndexFromOffset sh
    uncheckedIndexFromOffset :: sh -> Int -> Index sh
-   uncheckedIndexFromOffset = indexFromOffset
+   uncheckedIndexFromOffset sh = getUnchecked . unifiedIndexFromOffset sh
+   unifiedIndexFromOffset ::
+      (Checking check) => sh -> Int -> Result check (Index sh)
 
+assertIndexFromOffset ::
+   (Checking check) => String -> Int -> Bool -> Result check ()
+assertIndexFromOffset name k cond = assert (messageIndexFromOffset name k) cond
+
 class (C sh, Eq sh) => Static sh where
    static :: sh
 
@@ -114,7 +252,6 @@
 
 instance C Zero where
    size Zero = 0
-   uncheckedSize Zero = 0
 
 instance Static Zero where
    static = Zero
@@ -122,7 +259,6 @@
 
 instance C () where
    size () = 1
-   uncheckedSize () = 1
 
 {- |
 >>> Shape.indices ()
@@ -131,14 +267,11 @@
 instance Indexed () where
    type Index () = ()
    indices () = [()]
-   offset () () = 0
-   uncheckedOffset () () = 0
+   unifiedOffset () () = pure 0
    inBounds () () = True
 
 instance InvIndexed () where
-   indexFromOffset () 0 = ()
-   indexFromOffset () k = errorIndexFromOffset "()" k
-   uncheckedIndexFromOffset () _ = ()
+   unifiedIndexFromOffset () k = assertIndexFromOffset "()" k (k==0)
 
 instance Static () where
    static = ()
@@ -171,28 +304,24 @@
 
 instance (Integral n) => C (ZeroBased n) where
    size (ZeroBased len) = fromIntegral len
-   uncheckedSize (ZeroBased len) = fromIntegral len
 
 instance (Integral n) => Indexed (ZeroBased n) where
    type Index (ZeroBased n) = n
    indices (ZeroBased len) = indices $ Shifted 0 len
-   offset (ZeroBased len) = offset $ Shifted 0 len
-   uncheckedOffset _ ix = fromIntegral ix
+   unifiedOffset (ZeroBased len) = unifiedOffset $ Shifted 0 len
    inBounds (ZeroBased len) ix = 0<=ix && ix<len
 
 instance (Integral n) => InvIndexed (ZeroBased n) where
-   indexFromOffset (ZeroBased len) k0 =
+   unifiedIndexFromOffset (ZeroBased len) k0 = do
       let k = fromIntegral k0
-      in  if 0<=k && k<len
-            then k
-            else errorIndexFromOffset "ZeroBased" k0
-   uncheckedIndexFromOffset _ k = fromIntegral k
+      assertIndexFromOffset "ZeroBased" k0 $ 0<=k && k<len
+      pure k
 
-zeroBasedSplit :: (Real n) => n -> ZeroBased n -> ZeroBased n :+: ZeroBased n
+zeroBasedSplit :: (Real n) => n -> ZeroBased n -> ZeroBased n ::+ ZeroBased n
 zeroBasedSplit n (ZeroBased m) =
    if n<0
       then error "Shape.zeroBasedSplit: negative number of elements"
-      else let k = min n m in ZeroBased k :+: ZeroBased (m-k)
+      else let k = min n m in ZeroBased k ::+ ZeroBased (m-k)
 
 
 {- |
@@ -222,22 +351,18 @@
 
 instance (Integral n) => C (OneBased n) where
    size (OneBased len) = fromIntegral len
-   uncheckedSize (OneBased len) = fromIntegral len
 
 instance (Integral n) => Indexed (OneBased n) where
    type Index (OneBased n) = n
    indices (OneBased len) = indices $ Shifted 1 len
-   offset (OneBased len) = offset $ Shifted 1 len
-   uncheckedOffset _ ix = fromIntegral ix - 1
+   unifiedOffset (OneBased len) = unifiedOffset $ Shifted 1 len
    inBounds (OneBased len) ix = 0<ix && ix<=len
 
 instance (Integral n) => InvIndexed (OneBased n) where
-   indexFromOffset (OneBased len) k0 =
+   unifiedIndexFromOffset (OneBased len) k0 = do
       let k = fromIntegral k0
-      in  if 0<=k && k<len
-            then 1+k
-            else errorIndexFromOffset "OneBased" k0
-   uncheckedIndexFromOffset _ k = 1 + fromIntegral k
+      assertIndexFromOffset "OneBased" k0 $ 0<=k && k<len
+      pure $ 1+k
 
 
 {- |
@@ -262,22 +387,22 @@
 
 instance (Ix.Ix n) => C (Range n) where
    size (Range from to) = Ix.rangeSize (from,to)
-   uncheckedSize (Range from to) = Ix.unsafeRangeSize (from,to)
 
 instance (Ix.Ix n) => Indexed (Range n) where
    type Index (Range n) = n
    indices (Range from to) = Ix.range (from,to)
    offset (Range from to) ix = Ix.index (from,to) ix
    uncheckedOffset (Range from to) ix = Ix.unsafeIndex (from,to) ix
+   unifiedOffset (Range from to) ix = do
+      assert "Shape.Range: index out of range" $ Ix.inRange (from,to) ix
+      return $ Ix.unsafeIndex (from,to) ix
    inBounds (Range from to) ix = Ix.inRange (from,to) ix
 
 -- pretty inefficient when we rely solely on Ix
 instance (Ix.Ix n) => InvIndexed (Range n) where
-   indexFromOffset (Range from to) k =
-      if 0<=k && k < Ix.rangeSize (from,to)
-         then Ix.range (from,to) !! k
-         else errorIndexFromOffset "Range" k
-   uncheckedIndexFromOffset (Range from to) k = Ix.range (from,to) !! k
+   unifiedIndexFromOffset (Range from to) k = do
+      assertIndexFromOffset "Range" k $ 0<=k && k < Ix.rangeSize (from,to)
+      return $ Ix.range (from,to) !! k
 
 -- cf. sample-frame:Stereo
 instance Storable n => Storable (Range n) where
@@ -312,7 +437,6 @@
 
 instance (Integral n) => C (Shifted n) where
    size (Shifted _offs len) = fromIntegral len
-   uncheckedSize (Shifted _offs len) = fromIntegral len
 
 instance (Integral n) => Indexed (Shifted n) where
    type Index (Shifted n) = n
@@ -322,24 +446,18 @@
       zip
          (iterate (subtract 1) len)
          (iterate (1+) offs)
-   offset (Shifted offs len) ix =
-      if ix<offs
-        then error "Shape.Shifted: array index too small"
-        else
-          let k = ix-offs
-          in  if k<len
-                then fromIntegral k
-                else error "Shape.Shifted: array index too big"
-   uncheckedOffset (Shifted offs _len) ix = fromIntegral $ ix-offs
+   unifiedOffset (Shifted offs len) ix = do
+      assert "Shape.Shifted: array index too small" $ ix>=offs
+      let k = ix-offs
+      assert "Shape.Shifted: array index too big" $ k<len
+      return $ fromIntegral k
    inBounds (Shifted offs len) ix = offs <= ix && ix < offs+len
 
 instance (Integral n) => InvIndexed (Shifted n) where
-   indexFromOffset (Shifted offs len) k0 =
+   unifiedIndexFromOffset (Shifted offs len) k0 = do
       let k = fromIntegral k0
-      in  if 0<=k && k<len
-            then offs+k
-            else errorIndexFromOffset "Shifted" k0
-   uncheckedIndexFromOffset (Shifted offs _len) k = offs + fromIntegral k
+      assertIndexFromOffset "Shifted" k0 $ 0<=k && k<len
+      return $ offs+k
 
 -- cf. sample-frame:Stereo
 instance Storable n => Storable (Shifted n) where
@@ -380,22 +498,20 @@
    rnf Enumeration = ()
 
 instance (Enum n, Bounded n) => C (Enumeration n) where
-   size = uncheckedSize
-   uncheckedSize sh = intFromEnum sh maxBound - intFromEnum sh minBound + 1
+   size 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
+   unifiedOffset sh ix = pure $ 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
+   unifiedIndexFromOffset sh k = do
+      let minBnd = intFromEnum sh minBound
+      assertIndexFromOffset "Enumeration" k $
+         0<=k && k <= intFromEnum sh maxBound - minBnd
+      return $ toEnum $ minBnd + k
 
 asEnumType :: Enumeration n -> n -> n
 asEnumType Enumeration = id
@@ -418,26 +534,41 @@
 
 
 instance (Ord n) => C (Set n) where
-   size = uncheckedSize
-   uncheckedSize = Set.size
+   size = Set.size
 
 {- |
 You can use an arbitrary 'Set' of indices as shape.
 The array elements are ordered according to the index order in the 'Set'.
 
+Disadvantage is that combinators of different Set indexed arrays
+have to compare whole sets.
+However, the Set implementation may have low-level optimizations
+for pointer equality.
+
 >>> Shape.indices (Set.fromList "comfort")
 "cfmort"
 -}
 instance (Ord n) => Indexed (Set n) where
    type Index (Set n) = n
    indices = Set.toAscList
-   offset = ShapeSet.offset
-   uncheckedOffset = ShapeSet.uncheckedOffset
+   unifiedOffset sh ix = withCheck $ \check ->
+      case check of
+         Unchecked -> pure $ ShapeSet.uncheckedOffset sh ix
+         Checked ->
+            case ShapeSet.offset sh ix of
+               Just k -> pure k
+               Nothing ->
+                  throw "Shape.Set: array index not member of the index set"
    inBounds = flip Set.member
 
 instance (Ord n) => InvIndexed (Set n) where
-   indexFromOffset = ShapeSet.indexFromOffset
-   uncheckedIndexFromOffset = ShapeSet.uncheckedIndexFromOffset
+   unifiedIndexFromOffset sh k = withCheck $ \check ->
+      case check of
+         Unchecked -> pure $ ShapeSet.uncheckedIndexFromOffset sh k
+         Checked ->
+            case ShapeSet.indexFromOffset sh k of
+               Just ix -> pure ix
+               Nothing -> throw $ messageIndexFromOffset "Set" k
 
 
 {- |
@@ -445,7 +576,6 @@
 -}
 instance (Ord k, C shape) => C (Map k shape) where
    size = Fold.sum . Map.map size
-   uncheckedSize = Fold.sum . Map.map uncheckedSize
 
 {- |
 The implementations of 'offset' et.al.
@@ -458,54 +588,45 @@
    type Index (Map k shape) = (k, Index shape)
    indices =
       Fold.fold . Map.mapWithKey (\k shape -> map ((,) k) $ indices shape)
-   offset m =
-      let mu = snd $ Trav.mapAccumL (\l sh -> (l + size sh, (l,sh))) 0 m
-      in \(k,ix) ->
-         case Map.lookup k mu of
-            Nothing -> error "Shape.Map.offset: unknown key"
-            Just (l,sh) -> l + offset sh ix
-   uncheckedOffset m =
-      let mu =
-            snd $ Trav.mapAccumL (\l sh -> (l + uncheckedSize sh, (l,sh))) 0 m
+   unifiedOffset m =
+      let ms = fmap unifiedSizeOffset m
+          mu = snd $
+            Trav.mapAccumL (\l (sz,getOffset) -> (l + sz, (l,getOffset))) 0 ms
       in \(k,ix) ->
          case Map.lookup k mu of
-            Nothing -> error "Shape.Map.uncheckedOffset: unknown key"
-            Just (l,sh) -> l + uncheckedOffset sh ix
+            Nothing -> throwOrError "Shape.Map.offset: unknown key"
+            Just (l,getOffset) -> (l+) <$> getOffset ix
    inBounds m (k,ix) = Fold.any (flip inBounds ix) $ Map.lookup k m
 
-   sizeOffset = mapSizeOffset . Map.map sizeOffset
-   uncheckedSizeOffset = mapSizeOffset . Map.map uncheckedSizeOffset
+   unifiedSizeOffset = mapSizeOffset . fmap unifiedSizeOffset
 
 {-# INLINE mapSizeOffset #-}
-mapSizeOffset :: (Ord k, Num i) => Map k (i, ix -> i) -> (i, (k, ix) -> i)
+mapSizeOffset ::
+   (Checking check, Ord k, Num i) =>
+   Map k (i, ix -> Result check i) -> (i, (k, ix) -> Result check i)
 mapSizeOffset ms =
    (Fold.sum $ Map.map fst ms,
-    let mu = snd $ Trav.mapAccumL (\l (sz,offs) -> (l + sz, (l+) . offs)) 0 ms
+    let mu = snd $
+         Trav.mapAccumL (\l (sz,offs) -> (l + sz, fmap (l+) . offs)) 0 ms
     in \(k,ix) ->
-         fromMaybe (error "Shape.Map.sizeOffset: unknown key")
-            (Map.lookup k mu) ix)
+         maybe
+            (throwOrError "Shape.Map.sizeOffset: unknown key")
+            ($ix)
+            (Map.lookup k mu))
 
 instance (Ord k, InvIndexed shape) => InvIndexed (Map k shape) where
-   indexFromOffset m i =
+   unifiedIndexFromOffset m i =
       (\xs ->
          case xs of
             (_u,ix):_ -> ix
-            [] -> errorIndexFromOffset "Map" i) $
+            [] -> throwOrError $ messageIndexFromOffset "Map" i) $
       dropWhile (\(u,_ix) -> u<=i) $ snd $
       List.mapAccumL
          (\l (k,sh) ->
             let u = l + size sh
-            in (u, (u, (k, indexFromOffset sh (i-l))))) 0 $
+            in (u, (u, (,) k <$> unifiedIndexFromOffset sh (i-l)))) 0 $
       Map.toAscList m
 
-   uncheckedIndexFromOffset m i =
-      (\((_u,ix):_) -> ix) $
-      dropWhile (\(u,_ix) -> u<=i) $ snd $
-      List.mapAccumL
-         (\l (k,sh) ->
-            let u = l + size sh
-            in (u, (u, (k, uncheckedIndexFromOffset sh (i-l))))) 0 $
-      Map.toAscList m
 
 
 {- |
@@ -554,19 +675,17 @@
 
 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 sh
    indices (Deferred sh) = map DeferredIndex $ take (size sh) [0 ..]
-   offset (Deferred sh) (DeferredIndex k) = offset (ZeroBased $ size sh) k
-   uncheckedOffset _ (DeferredIndex k) = k
-   sizeOffset (Deferred sh) =
+   unifiedOffset (Deferred sh) (DeferredIndex k) = withCheck $ \check ->
+      case check of
+         Checked -> unifiedOffset (ZeroBased $ size sh) k
+         Unchecked -> pure k
+   unifiedSizeOffset (Deferred sh) =
       mapSnd (\offs (DeferredIndex k) -> offs k) $
-      sizeOffset (ZeroBased $ size sh)
-   uncheckedSizeOffset (Deferred sh) =
-      mapSnd (\ _offs (DeferredIndex k) -> k) $
-      uncheckedSizeOffset (ZeroBased $ size sh)
+      unifiedSizeOffset (ZeroBased $ size sh)
    inBounds (Deferred sh) (DeferredIndex k) =
       inBounds (ZeroBased $ size sh) k
 
@@ -574,6 +693,11 @@
    indexFromOffset (Deferred sh) k =
       DeferredIndex $ indexFromOffset (ZeroBased $ size sh) k
    uncheckedIndexFromOffset _sh = DeferredIndex
+   unifiedIndexFromOffset (Deferred sh) k = withCheck $ \check ->
+      case check of
+         Unchecked -> pure $ DeferredIndex k
+         Checked ->
+            DeferredIndex <$> unifiedIndexFromOffset (ZeroBased $ size sh) k
 
 deferIndex :: (Indexed sh, Index sh ~ ix) => sh -> ix -> DeferredIndex sh
 deferIndex sh ix = DeferredIndex $ offset sh ix
@@ -598,22 +722,18 @@
 
 instance (C sh) => C (Tagged s sh) where
    size (Tagged sh) = size sh
-   uncheckedSize (Tagged sh) = uncheckedSize sh
 
 instance (Indexed sh) => Indexed (Tagged s sh) where
    type Index (Tagged s sh) = Tagged s (Index sh)
    indices (Tagged sh) = map Tagged $ indices sh
-   offset (Tagged sh) (Tagged k) = offset sh k
-   uncheckedOffset (Tagged sh) (Tagged k) = uncheckedOffset sh k
-   sizeOffset (Tagged sh) = mapSnd (. unTagged) $ sizeOffset sh
-   uncheckedSizeOffset (Tagged sh) =
-      mapSnd (. unTagged) $ uncheckedSizeOffset sh
+   unifiedOffset (Tagged sh) = unifiedOffset sh . unTagged
+   unifiedSizeOffset (Tagged sh) =
+      mapSnd (. unTagged) $ unifiedSizeOffset sh
    inBounds (Tagged sh) (Tagged k) = inBounds sh k
 
 instance (InvIndexed sh) => InvIndexed (Tagged s sh) where
-   indexFromOffset (Tagged sh) k = Tagged $ indexFromOffset sh k
-   uncheckedIndexFromOffset (Tagged sh) k =
-      Tagged $ uncheckedIndexFromOffset sh k
+   unifiedIndexFromOffset (Tagged sh) k =
+      Tagged <$> unifiedIndexFromOffset sh k
 
 instance (Static sh) => Static (Tagged s sh) where
    static = Tagged static
@@ -622,7 +742,6 @@
 
 instance (C sh0, C sh1) => C (sh0,sh1) where
    size (sh0,sh1) = size sh0 * size sh1
-   uncheckedSize (sh0,sh1) = uncheckedSize sh0 * uncheckedSize sh1
 
 {- |
 Row-major composition of two dimensions.
@@ -633,29 +752,22 @@
 instance (Indexed sh0, Indexed sh1) => Indexed (sh0,sh1) where
    type Index (sh0,sh1) = (Index sh0, Index sh1)
    indices (sh0,sh1) = Monad.lift2 (,) (indices sh0) (indices sh1)
-   offset (sh0,sh1) =
-      offset sh0 . fst
-      `combineOffset`
-      mapSnd (.snd) (sizeOffset sh1)
-   uncheckedOffset (sh0,sh1) =
-      uncheckedOffset sh0 . fst
+   unifiedOffset (sh0,sh1) =
+      (unifiedOffset sh0 . fst)
       `combineOffset`
-      mapSnd (.snd) (uncheckedSizeOffset sh1)
-   sizeOffset (sh0,sh1) =
-      mapSnd (.fst) (sizeOffset sh0)
-      `combineSizeOffset`
-      mapSnd (.snd) (sizeOffset sh1)
-   uncheckedSizeOffset (sh0,sh1) =
-      mapSnd (.fst) (uncheckedSizeOffset sh0)
+      (mapSnd (.snd) $ unifiedSizeOffset sh1)
+   unifiedSizeOffset (sh0,sh1) =
+      (mapSnd (.fst) $ unifiedSizeOffset sh0)
       `combineSizeOffset`
-      mapSnd (.snd) (uncheckedSizeOffset sh1)
+      (mapSnd (.snd) $ unifiedSizeOffset sh1)
    inBounds (sh0,sh1) (ix0,ix1) = inBounds sh0 ix0 && inBounds sh1 ix1
 
 instance (InvIndexed sh0, InvIndexed sh1) => InvIndexed (sh0,sh1) where
-   indexFromOffset (sh0,sh1) k =
-      runInvIndex k $ liftA2 (,) (pickLastIndex sh0) (pickIndex sh1)
-   uncheckedIndexFromOffset (sh0,sh1) k =
-      runInvIndex k $ liftA2 (,) (uncheckedPickLastIndex sh0) (pickIndex sh1)
+   unifiedIndexFromOffset (sh0,sh1) k = do
+      let (rix0,ix1) =
+            runInvIndex k $ liftA2 (,) (pickLastIndex sh0) (pickIndex sh1)
+      ix0 <- rix0
+      return (ix0,ix1)
 
 instance (Static sh0, Static sh1) => Static (sh0,sh1) where
    static = (static, static)
@@ -663,43 +775,35 @@
 
 instance (C sh0, C sh1, C sh2) => C (sh0,sh1,sh2) where
    size (sh0,sh1,sh2) = size sh0 * size sh1 * size sh2
-   uncheckedSize (sh0,sh1,sh2) =
-      uncheckedSize sh0 * uncheckedSize sh1 * uncheckedSize sh2
 
 instance (Indexed sh0, Indexed sh1, Indexed sh2) => Indexed (sh0,sh1,sh2) where
    type Index (sh0,sh1,sh2) = (Index sh0, Index sh1, Index sh2)
    indices (sh0,sh1,sh2) =
       Monad.lift3 (,,) (indices sh0) (indices sh1) (indices sh2)
-   uncheckedOffset (sh0,sh1,sh2) =
-      uncheckedOffset sh0 . fst3
+   unifiedOffset (sh0,sh1,sh2) =
+      (unifiedOffset sh0 . fst3)
       `combineOffset`
-      mapSnd (.snd3) (uncheckedSizeOffset sh1)
-      `combineSizeOffset`
-      mapSnd (.thd3) (uncheckedSizeOffset sh2)
-   sizeOffset (sh0,sh1,sh2) =
-      mapSnd (.fst3) (sizeOffset sh0)
-      `combineSizeOffset`
-      mapSnd (.snd3) (sizeOffset sh1)
+      (mapSnd (.snd3) $ unifiedSizeOffset sh1)
       `combineSizeOffset`
-      mapSnd (.thd3) (sizeOffset sh2)
-   uncheckedSizeOffset (sh0,sh1,sh2) =
-      mapSnd (.fst3) (uncheckedSizeOffset sh0)
+      (mapSnd (.thd3) $ unifiedSizeOffset sh2)
+   unifiedSizeOffset (sh0,sh1,sh2) =
+      (mapSnd (.fst3) $ unifiedSizeOffset sh0)
       `combineSizeOffset`
-      mapSnd (.snd3) (uncheckedSizeOffset sh1)
+      (mapSnd (.snd3) $ unifiedSizeOffset sh1)
       `combineSizeOffset`
-      mapSnd (.thd3) (uncheckedSizeOffset sh2)
+      (mapSnd (.thd3) $ unifiedSizeOffset sh2)
    inBounds (sh0,sh1,sh2) (ix0,ix1,ix2) =
       inBounds sh0 ix0 && inBounds sh1 ix1 && inBounds sh2 ix2
 
 instance
    (InvIndexed sh0, InvIndexed sh1, InvIndexed sh2) =>
       InvIndexed (sh0,sh1,sh2) where
-   indexFromOffset (sh0,sh1,sh2) k =
-      runInvIndex k $
-      liftA3 (,,) (pickLastIndex sh0) (pickIndex sh1) (pickIndex sh2)
-   uncheckedIndexFromOffset (sh0,sh1,sh2) k =
-      runInvIndex k $
-      liftA3 (,,) (uncheckedPickLastIndex sh0) (pickIndex sh1) (pickIndex sh2)
+   unifiedIndexFromOffset (sh0,sh1,sh2) k = do
+      let (rix0,ix1,ix2) =
+            runInvIndex k $
+            liftA3 (,,) (pickLastIndex sh0) (pickIndex sh1) (pickIndex sh2)
+      ix0 <- rix0
+      return (ix0,ix1,ix2)
 
 instance (Static sh0, Static sh1, Static sh2) => Static (sh0,sh1,sh2) where
    static = (static, static, static)
@@ -708,14 +812,10 @@
 runInvIndex k = flip MS.evalState k . Back.forwards
 
 pickLastIndex ::
-   (InvIndexed sh) => sh -> Back.Backwards (MS.State Int) (Index sh)
+   (Checking check, InvIndexed sh) =>
+   sh -> Back.Backwards (MS.State Int) (Result check (Index sh))
 pickLastIndex sh =
-   Back.Backwards $ MS.gets $ indexFromOffset sh
-
-uncheckedPickLastIndex ::
-   (InvIndexed sh) => sh -> Back.Backwards (MS.State Int) (Index sh)
-uncheckedPickLastIndex sh =
-   Back.Backwards $ MS.gets $ uncheckedIndexFromOffset sh
+   Back.Backwards $ MS.gets $ unifiedIndexFromOffset sh
 
 pickIndex :: (InvIndexed sh) => sh -> Back.Backwards (MS.State Int) (Index sh)
 pickIndex sh =
@@ -727,13 +827,18 @@
 infixr 7 `combineOffset`, `combineSizeOffset`
 
 {-# INLINE combineOffset #-}
-combineOffset :: Num a => (ix -> a) -> (a, ix -> a) -> ix -> a
-combineOffset offset0 (size1,offset1) ix = offset0 ix * size1 + offset1 ix
+combineOffset ::
+   (Applicative f, Num a) =>
+   (ix -> f a) -> (a, ix -> f a) -> (ix -> f a)
+combineOffset offset0 (size1,offset1) ix =
+   offset0 ix |* size1 |+| offset1 ix
 
 {-# INLINE combineSizeOffset #-}
-combineSizeOffset :: Num a => (a, ix -> a) -> (a, ix -> a) -> (a, ix -> a)
+combineSizeOffset ::
+   (Applicative f, Num a) =>
+   (a, ix -> f a) -> (a, ix -> f a) -> (a, ix -> f a)
 combineSizeOffset (size0,offset0) (size1,offset1) =
-   (size0*size1, \ix -> offset0 ix * size1 + offset1 ix)
+   (size0*size1, \ix -> offset0 ix |* size1 |+| offset1 ix)
 
 
 
@@ -765,24 +870,18 @@
 
 instance (C sh) => C (Square sh) where
    size (Square sh) = size sh ^ (2::Int)
-   uncheckedSize (Square sh) = uncheckedSize sh ^ (2::Int)
 
 instance (Indexed sh) => Indexed (Square sh) where
    type Index (Square sh) = (Index sh, Index sh)
    indices (Square sh) = indices (sh,sh)
-   offset (Square sh) = offset (sh,sh)
-   uncheckedOffset (Square sh) = uncheckedOffset (sh,sh)
-   sizeOffset (Square sh) =
-      let szo = sizeOffset sh
-      in mapSnd (.fst) szo `combineSizeOffset` mapSnd (.snd) szo
-   uncheckedSizeOffset (Square sh) =
-      let szo = uncheckedSizeOffset sh
+   unifiedSizeOffset (Square sh) =
+      let szo = unifiedSizeOffset sh
       in mapSnd (.fst) szo `combineSizeOffset` mapSnd (.snd) szo
    inBounds (Square sh) = inBounds (sh,sh)
 
 instance (InvIndexed sh) => InvIndexed (Square sh) where
-   indexFromOffset (Square sh) = indexFromOffset (sh,sh)
-   uncheckedIndexFromOffset (Square sh) = uncheckedIndexFromOffset (sh,sh)
+   unifiedIndexFromOffset (Square sh) =
+      unifiedIndexFromOffset (sh,sh)
 
 
 
@@ -814,22 +913,12 @@
 
 instance (C sh) => C (Cube sh) where
    size (Cube sh) = size sh ^ (3::Int)
-   uncheckedSize (Cube sh) = uncheckedSize sh ^ (3::Int)
 
 instance (Indexed sh) => Indexed (Cube sh) where
    type Index (Cube sh) = (Index sh, Index sh, Index sh)
    indices (Cube sh) = indices (sh,sh,sh)
-   offset (Cube sh) = offset (sh,sh,sh)
-   uncheckedOffset (Cube sh) = uncheckedOffset (sh,sh,sh)
-   sizeOffset (Cube sh) =
-      let szo = sizeOffset sh
-      in mapSnd (.fst3) szo
-         `combineSizeOffset`
-         mapSnd (.snd3) szo
-         `combineSizeOffset`
-         mapSnd (.thd3) szo
-   uncheckedSizeOffset (Cube sh) =
-      let szo = uncheckedSizeOffset sh
+   unifiedSizeOffset (Cube sh) =
+      let szo = unifiedSizeOffset sh
       in mapSnd (.fst3) szo
          `combineSizeOffset`
          mapSnd (.snd3) szo
@@ -838,8 +927,8 @@
    inBounds (Cube sh) = inBounds (sh,sh,sh)
 
 instance (InvIndexed sh) => InvIndexed (Cube sh) where
-   indexFromOffset (Cube sh) = indexFromOffset (sh,sh,sh)
-   uncheckedIndexFromOffset (Cube sh) = uncheckedIndexFromOffset (sh,sh,sh)
+   unifiedIndexFromOffset (Cube sh) =
+      unifiedIndexFromOffset (sh,sh,sh)
 
 
 
@@ -900,7 +989,6 @@
 
 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) =>
@@ -915,21 +1003,14 @@
                (NonEmpty.tail $ NonEmpty.inits ixs) ixs)
             (zipWith (\r cs -> map ((,) r) cs) ixs $ tails ixs)
 
-   uncheckedOffset sh = snd $ uncheckedSizeOffset sh
-
-   sizeOffset (Triangular part sz) =
-      let (n, getOffset) = sizeOffset sz
-      in (triangleSize n, \(rs,cs) ->
-            let r = getOffset rs
-                c = getOffset cs
-            in if compareIndices part r c
-                  then triangleOffset part n (r,c)
-                  else error "Shape.Triangular.sizeOffset: wrong array part")
-
-   uncheckedSizeOffset (Triangular part sz) =
-      let (n, getOffset) = uncheckedSizeOffset sz
-      in (triangleSize n, \(rs,cs) ->
-            triangleOffset part n (getOffset rs, getOffset cs))
+   unifiedSizeOffset (Triangular part sz) =
+      let (n, getOffset) = unifiedSizeOffset sz
+      in (triangleSize n, \(rs,cs) -> do
+         r <- getOffset rs
+         c <- getOffset cs
+         assert "Shape.Triangular.sizeOffset: wrong array part" $
+            compareIndices part r c
+         return $ triangleOffset part n (r,c))
 
    inBounds (Triangular part sz) ix@(r,c) =
       inBounds (sz,sz) ix
@@ -950,13 +1031,13 @@
    (TriangularPart part, InvIndexed size) =>
       InvIndexed (Triangular part size) where
 
-   indexFromOffset (Triangular part sz) k =
-      mapPair (indexFromOffset sz, indexFromOffset sz) $
-      caseTriangularPart part
+   unifiedIndexFromOffset (Triangular part sz) k =
+      let n = size sz in
+      AppHT.mapPair (unifiedIndexFromOffset sz, unifiedIndexFromOffset sz) $
+       caseTriangularPart part
          (let r = floor (triangleRootDouble k)
           in (r, k - triangleSize r))
-         (let n = size sz
-              triSize = triangleSize n
+         (let triSize = triangleSize n
               rr = ceiling (triangleRootDouble (triSize-k))
               r = n - rr
           in (r, k+r - (triSize - triangleSize rr)))
@@ -988,6 +1069,265 @@
 
 
 {- |
+Simplex is a generalization of 'Triangular' to more than two dimensions.
+Indices are tuples of fixed size
+with elements ordered in ascending, strictly ascending,
+descending or strictly descending order.
+\"Order\" refers to the index order in 'indices'.
+In order to avoid confusion we suggest that the order of 'indices'
+is consistent with '<='.
+
+Obviously, 'offset' implements ranking
+and 'indexFromOffset' implements unranking
+of combinations (in the combinatorial sense)
+with or without repetitions.
+
+>>> Shape.indices $ Shape.simplexAscending (replicate 3 Shape.AllDistinct) $ Shape.ZeroBased (4::Int)
+[[0,1,2],[0,1,3],[0,2,3],[1,2,3]]
+>>> Shape.indices $ Shape.simplexAscending (replicate 3 Shape.SomeRepetitive) $ Shape.ZeroBased (3::Int)
+[[0,0,0],[0,0,1],[0,0,2],[0,1,1],[0,1,2],[0,2,2],[1,1,1],[1,1,2],[1,2,2],[2,2,2]]
+>>> Shape.indices $ Shape.simplexAscending [Shape.Repetitive,Shape.Distinct,Shape.Repetitive] $ Shape.ZeroBased (4::Int)
+[[0,0,1],[0,0,2],[0,0,3],[0,1,2],[0,1,3],[0,2,3],[1,1,2],[1,1,3],[1,2,3],[2,2,3]]
+>>> Shape.indices $ Shape.simplexAscending [Shape.Repetitive,Shape.Distinct,Shape.Distinct] $ Shape.ZeroBased (4::Int)
+[[0,0,1],[0,0,2],[0,0,3],[0,1,2],[0,1,3],[0,2,3],[1,1,2],[1,1,3],[1,2,3],[2,2,3]]
+
+>>> Shape.indices $ Shape.simplexDescending (replicate 3 Shape.AllDistinct) $ Shape.ZeroBased (4::Int)
+[[2,1,0],[3,1,0],[3,2,0],[3,2,1]]
+>>> Shape.indices $ Shape.simplexDescending (replicate 3 Shape.SomeRepetitive) $ Shape.ZeroBased (3::Int)
+[[0,0,0],[1,0,0],[1,1,0],[1,1,1],[2,0,0],[2,1,0],[2,1,1],[2,2,0],[2,2,1],[2,2,2]]
+>>> Shape.indices $ Shape.simplexDescending [Shape.Repetitive,Shape.Distinct,Shape.Repetitive] $ Shape.ZeroBased (4::Int)
+[[1,1,0],[2,1,0],[2,2,0],[2,2,1],[3,1,0],[3,2,0],[3,2,1],[3,3,0],[3,3,1],[3,3,2]]
+>>> Shape.indices $ Shape.simplexDescending [Shape.Repetitive,Shape.Distinct,Shape.Distinct] $ Shape.ZeroBased (4::Int)
+[[1,1,0],[2,1,0],[2,2,0],[2,2,1],[3,1,0],[3,2,0],[3,2,1],[3,3,0],[3,3,1],[3,3,2]]
+-}
+data Simplex order coll f size =
+   Simplex {
+      simplexOrder :: SimplexOrder order,
+      simplexDimension :: f coll,
+      simplexSize :: size
+   }
+
+data Ascending
+data Descending
+data SimplexOrder order where
+   Ascending :: SimplexOrder Ascending
+   Descending :: SimplexOrder Descending
+
+instance Eq (SimplexOrder order) where
+   Ascending == Ascending = True
+   Descending == Descending = True
+
+instance Show (SimplexOrder order) where
+   show Ascending = "Ascending"
+   show Descending = "Descending"
+
+type SimplexAscending = Simplex Ascending
+type SimplexDescending = Simplex Descending
+
+simplexAscending :: f coll -> size -> SimplexAscending coll f size
+simplexAscending = Simplex Ascending
+
+simplexDescending :: f coll -> size -> SimplexDescending coll f size
+simplexDescending = Simplex Descending
+
+isAscending :: SimplexOrder order -> Bool
+isAscending Ascending = True
+isAscending Descending = False
+
+class SimplexOrderC order where
+instance SimplexOrderC Ascending where
+instance SimplexOrderC Descending where
+
+data AllDistinct = AllDistinct deriving (Show, Eq)
+data SomeRepetitive = SomeRepetitive deriving (Show, Eq)
+data Collision = Distinct | Repetitive deriving (Show, Eq, Ord, Enum)
+
+class CollisionC coll where repetitionAllowed :: coll -> Bool
+instance CollisionC AllDistinct where repetitionAllowed AllDistinct = False
+instance CollisionC SomeRepetitive where repetitionAllowed SomeRepetitive = True
+instance CollisionC Collision where
+   repetitionAllowed Distinct = False
+   repetitionAllowed Repetitive = True
+
+instance
+   (SimplexOrderC order, Show coll, FunctorC.Show1 f, Show size) =>
+      Show (Simplex order coll f size) where
+   showsPrec p (Simplex order d sz) =
+      showParen (p>10) $
+         showString "Simplex " .
+         shows order .
+         showString " " .
+         FunctorC.showsPrec1 11 d .
+         showString " " .
+         showsPrec 11 sz
+
+instance
+   (SimplexOrderC order, CollisionC coll, Traversable f, C size) =>
+      C (Simplex order coll f size) where
+   size (Simplex _order d sz) =
+      let ds = Fold.toList d
+          rep = length $ filter repetitionAllowed $ laxInit ds
+      in simplexLayoutSize (length ds) (size sz + rep)
+
+laxInit :: [a] -> [a]
+laxInit xs = Match.take (drop 1 xs) xs
+
+simplexLayoutSize :: Integral i => Int -> i -> i
+simplexLayoutSize d n =
+   case drop d $ binomials n of
+      [] -> 0
+      m:_ -> m
+
+-- cf. package combinatorial
+binomials :: Integral a => a -> [a]
+binomials n =
+   scanl (\acc (num,den) -> div (acc*num) den) 1
+         (zip [n, pred n ..] [1..n])
+
+foldLength :: (Foldable f) => f a -> Int
+foldLength = length . Fold.toList
+
+instance
+   (SimplexOrderC order, CollisionC coll,
+    Traversable f, FunctorC.Eq1 f, Indexed size) =>
+      Indexed (Simplex order coll f size) where
+   type Index (Simplex order coll f size) = f (Index size)
+   indices (Simplex order d sz) =
+      flip MS.evalStateT (indices sz) $
+      Trav.traverse
+         (if isAscending order
+             then chooseIndexAscending
+             else chooseIndexDescending)
+         d
+   inBounds (Simplex order d sz) =
+      let getOffset = offset sz in \ix ->
+      let ixs = Fold.toList ix in
+         all (inBounds sz) ixs &&
+         FunctorC.eq1 (void d) (void ix) &&
+         isMonotonic order (Fold.toList d) (map getOffset ixs)
+   unifiedSizeOffset (Simplex order d sz) =
+      let (n, getOffset) = unifiedSizeOffset sz in
+      let dInt = foldLength d
+          prep = prepareSimplexIndexingOrder order d n in
+      (simplexLayoutSize dInt (fst prep),
+          -- cf. Combinatorics.chooseRank
+          \ixf -> do
+            ks <- Trav.traverse getOffset $ Fold.toList ixf
+            assert
+               "Shape.Simplex.offset: simplex and index structure mismatch"
+               (FunctorC.eq1 (void d) (void ixf))
+            assert
+               "Shape.Simplex.offset: index elements not monotonic"
+               (isMonotonic order (Fold.toList d) ks)
+            return $
+               simplexOffset order dInt
+                  (mapSnd (map snd . Fold.toList) prep) ks)
+
+simplexOffset ::
+   (Integral i) => SimplexOrder order -> Int -> (i, [(Int, i)]) -> [i] -> i
+simplexOffset order d (nsum,cis) ks =
+   case order of
+      Ascending ->
+         simplexLayoutSize d nsum - 1
+         -
+         sum (zipWith (\k (x,y) -> simplexLayoutSize x (y-k)) ks cis)
+      Descending ->
+         sum (zipWith (\k (x,y) -> simplexLayoutSize x (y+k)) ks cis)
+
+isMonotonic ::
+   (CollisionC coll) => SimplexOrder order -> [coll] -> [Int] -> Bool
+isMonotonic order cs =
+   and
+   .
+   (if isAscending order
+      then
+         ListHT.mapAdjacent
+            (\(c,x) (_,y) -> if repetitionAllowed c then x<=y else x<y)
+      else
+         ListHT.mapAdjacent
+            (\(c,x) (_,y) -> if repetitionAllowed c then x>=y else x>y))
+   .
+   zip cs
+
+chooseIndexAscending, chooseIndexDescending ::
+   (CollisionC coll) => coll -> MS.StateT [a] [] a
+
+chooseIndexAscending coll =
+   MS.StateT $ \as -> zip as $
+      (if repetitionAllowed coll then NonEmpty.flatten else NonEmpty.tail) $
+      NonEmpty.tails as
+
+chooseIndexDescending coll =
+   MS.StateT $ \as -> zip as $
+      (if repetitionAllowed coll then NonEmpty.tail else NonEmpty.flatten) $
+      NonEmpty.inits as
+
+instance
+   (SimplexOrderC order, CollisionC coll,
+    Traversable f, FunctorC.Eq1 f, InvIndexed size) =>
+      InvIndexed (Simplex order coll f size) where
+   unifiedIndexFromOffset (Simplex order d sh) =
+      let n = size sh in
+      let (nSum,deco) = prepareSimplexIndexingOrder order d n in
+      let dInt = foldLength d in \k ->
+      maybe
+         (throwOrError $ messageIndexFromOffset "Simplex" k)
+         (Trav.traverse (unifiedIndexFromOffset sh) . snd) $
+      if isAscending order
+         then
+            mapAccumLM
+               (\(a,k0) (db,(x,y)) ->
+                  case dropWhile ((<0) . snd) $
+                        map (\bi -> (bi, k0 - simplexLayoutSize x (y-bi))) $
+                        takeWhile (<n) $ iterate (1+) a of
+                     [] -> Nothing
+                     (b,k1):_ -> Just ((b+db, k1), b))
+               (0, simplexLayoutSize dInt nSum - 1 - k)
+               deco
+         else
+            mapAccumLM
+               (\(a,k0) (db,(x,y)) ->
+                  case dropWhile ((<0) . snd) $
+                        map (\bi -> (bi, k0 - simplexLayoutSize x (y+bi))) $
+                        takeWhile (>=0) $ iterate (subtract 1) a of
+                     [] -> Nothing
+                     (b,k1):_ -> Just ((b-db, k1), b))
+               (n,k)
+               deco
+
+mapAccumLM ::
+   (Traversable t, Monad m) => (a -> b -> m (a, c)) -> a -> t b -> m (a, t c)
+mapAccumLM f a0 xs =
+   liftM swap $
+   MS.runStateT
+      (Trav.mapM (\b -> MS.StateT $ \a -> liftM swap $ f a b) xs) a0
+
+
+prepareSimplexIndexingOrder ::
+   (Traversable t, Num i, CollisionC coll) =>
+   SimplexOrder order -> t coll -> Int -> (Int, t (Int, (i, Int)))
+prepareSimplexIndexingOrder order d n =
+   if isAscending order
+      then mapFst (1+) $ prepareSimplexIndexing d (n-1)
+      else mapFst (n+) $ prepareSimplexIndexing d 0
+
+prepareSimplexIndexing ::
+   (Traversable t, Num i, CollisionC coll) =>
+   t coll -> Int -> (Int, t (Int, (i, Int)))
+prepareSimplexIndexing d n =
+   let ((_,(_,nSum)), deco) =
+         Trav.mapAccumR
+            (\(c0,(x,y)) ci ->
+               let c1 = fromEnum (ci&&c0)
+                   p = (x+1,y+c1)
+               in ((True,p),(1-c1,p)))
+            (False,(0,n))
+            (fmap repetitionAllowed d)
+   in (nSum, deco)
+
+
+
+{- |
 'Cyclic' is a shape, where the indices wrap around at the array boundaries.
 E.g.
 
@@ -1020,75 +1360,71 @@
 
 instance (Integral n) => C (Cyclic n) where
    size (Cyclic len) = fromIntegral len
-   uncheckedSize (Cyclic len) = fromIntegral len
 
 instance (Integral n) => Indexed (Cyclic n) where
    type Index (Cyclic n) = n
    indices (Cyclic len) = indices $ ZeroBased len
-   offset = uncheckedOffset
-   uncheckedOffset (Cyclic len) ix = fromIntegral $ mod ix len
+   unifiedOffset (Cyclic len) ix = pure $ fromIntegral $ mod ix len
    inBounds (Cyclic len) _ix = len>0
 
 instance (Integral n) => InvIndexed (Cyclic n) where
-   indexFromOffset (Cyclic len) k0 =
+   unifiedIndexFromOffset (Cyclic len) k0 = do
       let k = fromIntegral k0
-      in  if 0<=k && k<len
-            then k
-            else errorIndexFromOffset "Cyclic" k0
-   uncheckedIndexFromOffset _ k = fromIntegral k
+      assertIndexFromOffset "Cyclic" k0 $ 0<=k && k<len
+      return k
 
 
 
-infixr 5 :+:
+infixr 5 ::+
 
 {- |
 Row-major composition of two dimensions.
 
->>> Shape.indices (Shape.ZeroBased (3::Int) :+: Shape.Range 'a' 'c')
+>>> Shape.indices (Shape.ZeroBased (3::Int) ::+ Shape.Range 'a' 'c')
 [Left 0,Left 1,Left 2,Right 'a',Right 'b',Right 'c']
 -}
-data sh0:+:sh1 = sh0:+:sh1
+data sh0::+sh1 = sh0::+sh1
    deriving (Eq, Show)
 
-instance (NFData sh0, NFData sh1) => NFData (sh0:+:sh1) where
-   rnf (sh0:+:sh1) = rnf (sh0,sh1)
+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
-   uncheckedSize (sh0:+:sh1) = uncheckedSize sh0 + uncheckedSize sh1
+instance (C sh0, C sh1) => C (sh0::+sh1) where
+   size (sh0::+sh1) = size sh0 + size sh1
 
-instance (Indexed sh0, Indexed sh1) => Indexed (sh0:+:sh1) where
-   type Index (sh0:+:sh1) = Either (Index sh0) (Index sh1)
-   indices (sh0:+:sh1) = map Left (indices sh0) ++ map Right (indices sh1)
-   offset (sh0:+:sh1) ix =
-      case ix of
-         Left ix0 -> offset sh0 ix0
-         Right ix1 -> size sh0 + offset sh1 ix1
-   uncheckedOffset (sh0:+:sh1) ix =
-      case ix of
-         Left ix0 -> uncheckedOffset sh0 ix0
-         Right ix1 -> uncheckedSize sh0 + uncheckedOffset sh1 ix1
-   sizeOffset (sh0:+:sh1) =
-      let (n0, getOffset0) = sizeOffset sh0
-          (n1, getOffset1) = sizeOffset sh1
-      in (n0+n1, either getOffset0 ((n0+) . getOffset1))
-   uncheckedSizeOffset (sh0:+:sh1) =
-      let (n0, getOffset0) = uncheckedSizeOffset sh0
-          (n1, getOffset1) = uncheckedSizeOffset sh1
-      in (n0+n1, either getOffset0 ((n0+) . getOffset1))
-   inBounds (sh0:+:sh1) = either (inBounds sh0) (inBounds sh1)
+instance (Indexed sh0, Indexed sh1) => Indexed (sh0::+sh1) where
+   type Index (sh0::+sh1) = Either (Index sh0) (Index sh1)
+   indices (sh0::+sh1) = map Left (indices sh0) ++ map Right (indices sh1)
+   unifiedOffset (sh0::+sh1) =
+      let (n0,getOffset0) = unifiedSizeOffset sh0
+          getOffset1 = unifiedOffset sh1
+      in \ix ->
+         case ix of
+            Left ix0 -> getOffset0 ix0
+            Right ix1 -> (n0 +) <$> getOffset1 ix1
+   unifiedSizeOffset (sh0::+sh1) =
+      let (n0, getOffset0) = unifiedSizeOffset sh0
+          (n1, getOffset1) = unifiedSizeOffset sh1
+      in (n0+n1, either getOffset0 (fmap (n0+) . getOffset1))
+   inBounds (sh0::+sh1) = either (inBounds sh0) (inBounds sh1)
 
-instance (InvIndexed sh0, InvIndexed sh1) => InvIndexed (sh0:+:sh1) where
-   indexFromOffset (sh0:+:sh1) k =
-      let pivot = size sh0
-      in if k < pivot
-            then Left $ indexFromOffset sh0 k
-            else Right $ indexFromOffset sh1 $ k-pivot
-   uncheckedIndexFromOffset (sh0:+:sh1) k =
+instance (InvIndexed sh0, InvIndexed sh1) => InvIndexed (sh0::+sh1) where
+   unifiedIndexFromOffset (sh0::+sh1) =
       let pivot = size sh0
-      in if k < pivot
-            then Left $ uncheckedIndexFromOffset sh0 k
-            else Right $ uncheckedIndexFromOffset sh1 $ k-pivot
+      in \k ->
+         if k < pivot
+            then Left <$> unifiedIndexFromOffset sh0 k
+            else Right <$> unifiedIndexFromOffset sh1 (k-pivot)
 
-instance (Static sh0, Static sh1) => Static (sh0:+:sh1) where
-   static = static:+:static
+instance (Static sh0, Static sh1) => Static (sh0::+sh1) where
+   static = static::+static
+
+
+infixl 7 |*
+infixl 6 |+|
+
+(|*) :: (Functor f, Num a) => f a -> a -> f a
+f|*a = fmap (*a) f
+
+(|+|) :: (Applicative f, Num a) => f a -> f a -> f a
+(|+|) = liftA2 (+)
diff --git a/src/Data/Array/Comfort/Shape/Test.hs b/src/Data/Array/Comfort/Shape/Test.hs
--- a/src/Data/Array/Comfort/Shape/Test.hs
+++ b/src/Data/Array/Comfort/Shape/Test.hs
@@ -1,15 +1,16 @@
 {-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE GADTs #-}
 module Data.Array.Comfort.Shape.Test (tests) where
 
 import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Shape.Utility (isRight)
+
+import Control.Applicative (pure)
 import Data.Tuple.HT (mapSnd)
 
 import qualified Test.QuickCheck as QC
 
 
-uncheckedSize :: (Shape.C sh) => sh -> Bool
-uncheckedSize sh  =  Shape.size sh == Shape.uncheckedSize sh
-
 inBounds :: (Shape.Indexed sh) => sh -> Bool
 inBounds sh  =  all (Shape.inBounds sh) $ Shape.indices sh
 
@@ -21,6 +22,21 @@
    let ixs = Shape.indices sh
    in not (null ixs)  QC.==>  QC.forAll (QC.elements ixs) f
 
+-- ToDo: we need to check for indices outside of bounds, too
+inBoundsOffset ::
+   (Shape.Indexed sh, Shape.Index sh ~ ix, Show ix) => sh -> QC.Property
+inBoundsOffset sh =
+   forAllIndices sh $ \ix ->
+      Shape.inBounds sh ix ==
+      isRight (Shape.getChecked (Shape.unifiedOffset sh ix))
+
+inBoundsSizeOffset ::
+   (Shape.Indexed sh, Shape.Index sh ~ ix, Show ix) => sh -> QC.Property
+inBoundsSizeOffset sh =
+   forAllIndices sh $ \ix ->
+      Shape.inBounds sh ix ==
+      isRight (Shape.getChecked (snd (Shape.unifiedSizeOffset sh) ix))
+
 sizeOffset ::
    (Shape.Indexed sh, Shape.Index sh ~ ix, Show ix) => sh -> QC.Property
 sizeOffset sh =
@@ -34,14 +50,45 @@
 uncheckedSizeOffset sh =
    forAllIndices sh $ \ix ->
       mapSnd ($ix) (Shape.uncheckedSizeOffset sh) ==
-         (Shape.uncheckedSize sh, Shape.uncheckedOffset sh ix)
+         (Shape.size sh, Shape.uncheckedOffset sh ix)
 
+unifiedSizeOffsetA ::
+   (Shape.Checking check, Shape.Indexed sh, Shape.Index sh ~ ix, Show ix) =>
+   Shape.CheckSingleton check -> sh -> QC.Property
+unifiedSizeOffsetA check sh =
+   forAllIndices sh $ \ix ->
+      mapSnd ($ix) (Shape.unifiedSizeOffset sh) ==
+         (Shape.size sh, Shape.requireCheck check $ Shape.unifiedOffset sh ix)
+
+unifiedSizeOffsetB ::
+   (Shape.Checking check, Shape.Indexed sh, Shape.Index sh ~ ix, Show ix) =>
+   Shape.CheckSingleton check -> sh -> QC.Property
+unifiedSizeOffsetB check sh =
+   forAllIndices sh $ \ix ->
+      (mapSnd (Shape.requireCheck check . ($ix)) $ Shape.unifiedSizeOffset sh)
+      ==
+      case check of
+         Shape.Checked ->
+            mapSnd (pure . ($ix)) (Shape.sizeOffset sh)
+         Shape.Unchecked ->
+            mapSnd (pure . ($ix)) (Shape.uncheckedSizeOffset sh)
+
 uncheckedOffset ::
    (Shape.Indexed sh, Shape.Index sh ~ ix, Show ix) => sh -> QC.Property
 uncheckedOffset sh =
    forAllIndices sh $ \ix ->
       Shape.offset sh ix == Shape.uncheckedOffset sh ix
 
+unifiedOffset ::
+   (Shape.Checking check, Shape.Indexed sh, Shape.Index sh ~ ix, Show ix) =>
+   Shape.CheckSingleton check -> sh -> QC.Property
+unifiedOffset check sh =
+   forAllIndices sh $ \ix ->
+      Shape.requireCheck check (Shape.unifiedOffset sh ix) ==
+      case check of
+         Shape.Checked -> pure $ Shape.offset sh ix
+         Shape.Unchecked -> pure $ Shape.uncheckedOffset sh ix
+
 lengthIndices :: (Shape.Indexed sh) => sh -> Bool
 lengthIndices sh  =  length (Shape.indices sh) == Shape.size sh
 
@@ -60,13 +107,50 @@
    Shape.indices sh ==
    map (Shape.uncheckedIndexFromOffset sh) (take (Shape.size sh) [0..])
 
+unifiedInvIndicesA ::
+   (Shape.Checking check, Shape.InvIndexed sh, Shape.Index sh ~ ix, Eq ix) =>
+   Shape.CheckSingleton check -> sh -> Bool
+unifiedInvIndicesA check sh =
+   map pure (Shape.indices sh) ==
+   map (Shape.requireCheck check . Shape.unifiedIndexFromOffset sh)
+      (take (Shape.size sh) [0..])
 
+unifiedInvIndicesB ::
+   (Shape.Checking check, Shape.InvIndexed sh, Shape.Index sh ~ ix, Eq ix) =>
+   Shape.CheckSingleton check -> sh -> QC.Property
+unifiedInvIndicesB check sh =
+   let n = Shape.size sh in n>0 QC.==>
+   QC.forAll (QC.choose (0, n-1)) $ \k ->
+   Shape.requireCheck check (Shape.unifiedIndexFromOffset sh k) ==
+   case check of
+      Shape.Checked -> pure $ Shape.indexFromOffset sh k
+      Shape.Unchecked -> pure $ Shape.uncheckedIndexFromOffset sh k
+
+
+unifiedTests ::
+   (Shape.Checking check,
+    Shape.InvIndexed sh, Show sh, Shape.Index sh ~ ix, Eq ix, Show ix) =>
+   Shape.CheckSingleton check -> QC.Gen sh -> [(String, QC.Property)]
+unifiedTests check gen =
+   ("unifiedSizeOffsetA", QC.forAll gen (unifiedSizeOffsetA check)) :
+   ("unifiedSizeOffsetB", QC.forAll gen (unifiedSizeOffsetB check)) :
+   ("unifiedOffset", QC.forAll gen (unifiedOffset check)) :
+   ("unifiedInvIndicesA", QC.forAll gen (unifiedInvIndicesA check)) :
+   ("unifiedInvIndicesB", QC.forAll gen (unifiedInvIndicesB check)) :
+   []
+
+-- cf. Test.Utility
+prefix :: String -> [(String, test)] -> [(String, test)]
+prefix msg =
+   map (\(str,test) -> (msg ++ "." ++ str, test))
+
 tests ::
    (Shape.InvIndexed sh, Show sh, Shape.Index sh ~ ix, Eq ix, Show ix) =>
    QC.Gen sh -> [(String, QC.Property)]
 tests gen =
-   ("uncheckedSize", QC.forAll gen uncheckedSize) :
    ("inBounds", QC.forAll gen inBounds) :
+   ("inBoundsOffset", QC.forAll gen inBoundsOffset) :
+   ("inBoundsSizeOffset", QC.forAll gen inBoundsSizeOffset) :
    ("sizeOffset", QC.forAll gen sizeOffset) :
    ("uncheckedSizeOffset", QC.forAll gen uncheckedSizeOffset) :
    ("uncheckedOffset", QC.forAll gen uncheckedOffset) :
@@ -74,4 +158,6 @@
    ("indexOffsets", QC.forAll gen indexOffsets) :
    ("invIndices", QC.forAll gen invIndices) :
    ("uncheckedInvIndices", QC.forAll gen uncheckedInvIndices) :
+   prefix "Checked" (unifiedTests Shape.Checked gen) ++
+   prefix "Unchecked" (unifiedTests Shape.Unchecked gen) ++
    []
diff --git a/src/Data/Array/Comfort/Shape/Utility.hs b/src/Data/Array/Comfort/Shape/Utility.hs
--- a/src/Data/Array/Comfort/Shape/Utility.hs
+++ b/src/Data/Array/Comfort/Shape/Utility.hs
@@ -3,6 +3,13 @@
 import Text.Printf (printf)
 
 
+messageIndexFromOffset :: String -> Int -> String
+messageIndexFromOffset name k =
+   printf "indexFromOffset (%s): index %d out of range" name k
+
 errorIndexFromOffset :: String -> Int -> a
-errorIndexFromOffset name k =
-   error $ printf "indexFromOffset (%s): index %d out of range" name k
+errorIndexFromOffset name = error . messageIndexFromOffset name
+
+
+isRight :: Either a b -> Bool
+isRight = either (const False) (const True)
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
@@ -5,7 +5,7 @@
    reshape,
    mapShape,
 
-   (!),
+   accessMaybe, (!),
    Array.toList,
    Array.vectorFromList,
    toAssociations,
@@ -16,6 +16,8 @@
    sample,
    fromBoxed,
    toBoxed,
+   fromStorableVector,
+   toStorableVector,
 
    Array.map,
    Array.mapWithIndex,
@@ -24,6 +26,9 @@
    accumulate,
    fromAssociations,
 
+   pick,
+   toRowArray,
+   fromRowArray,
    Array.singleton,
    Array.append,
    Array.take, Array.drop,
@@ -48,21 +53,26 @@
 import qualified Data.Array.Comfort.Check as Check
 import qualified Data.Array.Comfort.Shape as Shape
 import Data.Array.Comfort.Storable.Unchecked (Array(Array))
+import Data.Array.Comfort.Shape ((::+)((::+)))
 
 import System.IO.Unsafe (unsafePerformIO)
+import Foreign.Marshal.Array (copyArray, advancePtr)
 import Foreign.Storable (Storable)
 import Foreign.ForeignPtr (withForeignPtr)
 
 import Control.Monad.ST (runST)
 
+import qualified Data.StorableVector.Base as SVB
 import qualified Data.Map as Map
 import qualified Data.Set as Set
+import qualified Data.Traversable as Trav
 import qualified Data.Foldable as Fold
 import qualified Data.List as List
 import qualified Data.Tuple.Strict as StrictTuple
 import Data.Map (Map)
 import Data.Set (Set)
 import Data.Foldable (forM_)
+import Data.Maybe (fromMaybe)
 import Data.Semigroup
          (Semigroup, (<>), Min(Min,getMin), Max(Max,getMax), Arg(Arg))
 
@@ -73,12 +83,12 @@
 {- $setup
 >>> import qualified Data.Array.Comfort.Storable as Array
 >>> import qualified Data.Array.Comfort.Shape as Shape
->>> import Data.Array.Comfort.Storable (Array)
+>>> import Data.Array.Comfort.Storable (Array, (!))
 >>>
 >>> import qualified Test.QuickCheck as QC
 >>> import Test.ChasingBottoms.IsBottom (isBottom)
 >>>
->>> import Control.Applicative ((<$>))
+>>> import Control.Applicative ((<$>), (<*>))
 >>>
 >>> import Data.Word (Word16)
 >>>
@@ -87,6 +97,23 @@
 >>> genArray :: QC.Gen (Array ShapeInt Word16)
 >>> genArray = Array.vectorFromList <$> QC.arbitrary
 >>>
+>>> genArray2 :: QC.Gen (Array (ShapeInt,ShapeInt) Word16)
+>>> genArray2 = do
+>>>    xs <- QC.arbitrary
+>>>    let n = length xs
+>>>    (k,m) <-
+>>>       if n == 0
+>>>          then QC.elements [(,) 0, flip (,) 0] <*> QC.choose (1,n)
+>>>          else fmap (\m -> (div n m, m)) $ QC.choose (1,n)
+>>>    return $ Array.fromList (Shape.ZeroBased k, Shape.ZeroBased m) xs
+>>>
+>>> genNonEmptyArray2 :: QC.Gen (Array (ShapeInt,ShapeInt) Word16)
+>>> genNonEmptyArray2 = do
+>>>    xs <- QC.getNonEmpty <$> QC.arbitrary
+>>>    let n = length xs
+>>>    m <- QC.choose (1,n)
+>>>    return $ Array.fromList (Shape.ZeroBased (div n m), Shape.ZeroBased m) xs
+>>>
 >>> infix 4 ==?
 >>> (==?) :: a -> a -> (a,a)
 >>> (==?) = (,)
@@ -141,17 +168,42 @@
 toBoxed :: (Shape.C sh, Storable a) => Array sh a -> BoxedArray.Array sh a
 toBoxed arr = BoxedArray.fromList (Array.shape arr) $ Array.toList arr
 
+
+fromStorableVector ::
+   (Storable a) => SVB.Vector a -> Array (Shape.ZeroBased Int) a
+fromStorableVector xs =
+   case SVB.toForeignPtr xs of
+      (fptr,0,n) -> Array (Shape.ZeroBased n) fptr
+      (fptr,s,n) ->
+         Array.takeRight $
+         Array (Shape.ZeroBased s ::+ Shape.ZeroBased n) fptr
+
+toStorableVector :: (Shape.C sh, Storable a) => Array sh a -> SVB.Vector a
+toStorableVector (Array sh fptr) =
+   SVB.fromForeignPtr fptr $ Shape.size sh
+
+
 toAssociations ::
    (Shape.Indexed sh, Storable a) => Array sh a -> [(Shape.Index sh, a)]
 toAssociations arr = zip (Shape.indices $ shape arr) (Array.toList arr)
 
 
+errorArray :: String -> String -> a
+errorArray name msg =
+   error ("Array.Comfort.Storable." ++ name ++ ": " ++ msg)
+
 infixl 9 !
 
 (!) :: (Shape.Indexed sh, Storable a) => Array sh a -> Shape.Index sh -> a
-(!) arr ix = runST (do
+(!) arr ix =
+   fromMaybe (errorArray "!" "index out of bounds") $
+   accessMaybe arr ix
+
+accessMaybe ::
+   (Shape.Indexed sh, Storable a) => Array sh a -> Shape.Index sh -> Maybe a
+accessMaybe arr ix = runST (do
    marr <- MutArrayNC.unsafeThaw arr
-   MutArray.read marr ix)
+   Trav.sequenceA $ MutArray.readMaybe marr ix)
 
 
 zipWith ::
@@ -160,7 +212,7 @@
 zipWith f a b =
    if shape a == shape b
       then Array.zipWith f a b
-      else error "zipWith: shapes mismatch"
+      else errorArray "zipWith" "shapes mismatch"
 
 (//) ::
    (Shape.Indexed sh, Storable a) =>
@@ -186,6 +238,38 @@
    forM_ xs $ uncurry $ MutArray.write marr
    MutArrayNC.unsafeFreeze marr)
 
+
+{- |
+prop> QC.forAll genNonEmptyArray2 $ \xs -> QC.forAll (QC.elements $ Shape.indices $ Array.shape xs) $ \(ix0,ix1) -> Array.pick xs ix0 ! ix1 == xs!(ix0,ix1)
+-}
+pick ::
+   (Shape.Indexed sh0, Shape.C sh1, Storable a) =>
+   Array (sh0,sh1) a -> Shape.Index sh0 -> Array sh1 a
+pick (Array (sh0,sh1) x) ix0 =
+   Array.unsafeCreateWithSize sh1 $ \k yPtr ->
+   withForeignPtr x $ \xPtr ->
+      copyArray yPtr (advancePtr xPtr (Shape.offset sh0 ix0 * Shape.size sh1)) k
+
+toRowArray ::
+   (Shape.Indexed sh0, Shape.C sh1, Storable a) =>
+   Array (sh0,sh1) a -> BoxedArray.Array sh0 (Array sh1 a)
+toRowArray x = fmap (pick x) $ BoxedArray.indices $ fst $ Array.shape x
+
+{- |
+It is a checked error if a row width differs from the result array width.
+
+prop> QC.forAll genArray2 $ \xs -> xs == Array.fromRowArray (snd $ Array.shape xs) (Array.toRowArray xs)
+-}
+fromRowArray ::
+   (Shape.C sh0, Shape.C sh1, Eq sh1, Storable a) =>
+   sh1 -> BoxedArray.Array sh0 (Array sh1 a) -> Array (sh0,sh1) a
+fromRowArray sh1 x =
+   Array.unsafeCreate (BoxedArray.shape x, sh1) $ \yPtr ->
+   let k = Shape.size sh1 in
+   forM_ (zip [0,k..] (BoxedArray.toList x)) $ \(j, Array sh1i row) ->
+   if sh1 == sh1i
+      then withForeignPtr row $ \xPtr -> copyArray (advancePtr yPtr j) xPtr k
+      else errorArray "fromRowArray" "mismatching row width"
 
 
 {- |
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
@@ -5,7 +5,7 @@
    shape,
 
    MutArray.new,
-   read,
+   read, MutArray.readMaybe,
    write,
    update,
    toList,
@@ -19,6 +19,7 @@
 import qualified Data.Array.Comfort.Storable.Mutable.Unchecked as MutArray
 import qualified Data.Array.Comfort.Shape as Shape
 import Data.Array.Comfort.Storable.Mutable.Unchecked (Array)
+import Data.Maybe (fromMaybe)
 
 import Foreign.Marshal.Array (pokeArray)
 import Foreign.Storable (Storable)
@@ -36,9 +37,9 @@
    (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"
+   fromMaybe
+      (error "Array.Comfort.Storable.Mutable.read: index out of bounds")
+      (MutArray.readMaybe arr ix)
 
 write ::
    (PrimMonad m, Shape.Indexed sh, Storable a) =>
diff --git a/src/Data/Array/Comfort/Storable/Mutable/Private.hs b/src/Data/Array/Comfort/Storable/Mutable/Private.hs
--- a/src/Data/Array/Comfort/Storable/Mutable/Private.hs
+++ b/src/Data/Array/Comfort/Storable/Mutable/Private.hs
@@ -13,6 +13,7 @@
 import Control.Monad (liftM)
 import Control.Applicative ((<$>))
 
+import Data.Either.HT (maybeRight)
 import Data.Tuple.HT (mapFst)
 
 import qualified Prelude as P
@@ -100,6 +101,13 @@
    Array m sh a -> Shape.Index sh -> m a
 read (Array sh fptr) ix =
    withArrayPtr fptr $ flip peekElemOff (Shape.uncheckedOffset sh ix)
+
+readMaybe ::
+   (PrimMonad m, Shape.Indexed sh, Storable a) =>
+   Array m sh a -> Shape.Index sh -> Maybe (m a)
+readMaybe (Array sh fptr) ix =
+   fmap (withArrayPtr fptr . flip peekElemOff) $ maybeRight $
+   Shape.getChecked $ Shape.unifiedOffset sh ix
 
 write ::
    (PrimMonad m, Shape.Indexed sh, Storable a) =>
diff --git a/src/Data/Array/Comfort/Storable/Mutable/Unchecked.hs b/src/Data/Array/Comfort/Storable/Mutable/Unchecked.hs
--- a/src/Data/Array/Comfort/Storable/Mutable/Unchecked.hs
+++ b/src/Data/Array/Comfort/Storable/Mutable/Unchecked.hs
@@ -15,6 +15,7 @@
    MutArray.unsafeCreateWithSizeAndResult,
    MutArray.withPtr,
    MutArray.read,
+   MutArray.readMaybe,
    MutArray.write,
    MutArray.update,
    MutArray.toList,
diff --git a/src/Data/Array/Comfort/Storable/Unchecked.hs b/src/Data/Array/Comfort/Storable/Unchecked.hs
--- a/src/Data/Array/Comfort/Storable/Unchecked.hs
+++ b/src/Data/Array/Comfort/Storable/Unchecked.hs
@@ -38,7 +38,7 @@
 import qualified Data.Array.Comfort.Storable.Memory as Memory
 import qualified Data.Array.Comfort.Shape as Shape
 import Data.Array.Comfort.Storable.Private (Array(Array), mapShape)
-import Data.Array.Comfort.Shape ((:+:)((:+:)))
+import Data.Array.Comfort.Shape ((::+)((::+)))
 
 import System.IO.Unsafe (unsafePerformIO)
 import Foreign.Marshal.Array (copyArray, advancePtr)
@@ -136,11 +136,14 @@
 singleton :: (Storable a) => a -> Array () a
 singleton a = unsafeCreate () $ flip poke a
 
+
+infixr 5 `append`
+
 append ::
    (Shape.C shx, Shape.C shy, Storable a) =>
-   Array shx a -> Array shy a -> Array (shx:+:shy) a
+   Array shx a -> Array shy a -> Array (shx::+shy) a
 append (Array shX x) (Array shY y) =
-   unsafeCreate (shX:+:shY) $ \zPtr ->
+   unsafeCreate (shX::+shY) $ \zPtr ->
    withForeignPtr x $ \xPtr ->
    withForeignPtr y $ \yPtr -> do
       let sizeX = Shape.size shX
@@ -160,7 +163,7 @@
 splitN ::
    (Integral n, Storable a) =>
    n -> Array (Shape.ZeroBased n) a ->
-   Array (Shape.ZeroBased n :+: Shape.ZeroBased n) a
+   Array (Shape.ZeroBased n ::+ Shape.ZeroBased n) a
 splitN n = mapShape (Shape.zeroBasedSplit n)
 
 {- |
@@ -168,19 +171,19 @@
 -}
 takeLeft ::
    (Shape.C sh0, Shape.C sh1, Storable a) =>
-   Array (sh0:+:sh1) a -> Array sh0 a
+   Array (sh0::+sh1) a -> Array sh0 a
 takeLeft =
-   takeCenter . mapShape (\(sh0 :+: sh1) -> (Shape.Zero :+: sh0 :+: sh1))
+   takeCenter . mapShape (\(sh0 ::+ sh1) -> (Shape.Zero ::+ sh0 ::+ sh1))
 
 takeRight ::
    (Shape.C sh0, Shape.C sh1, Storable a) =>
-   Array (sh0:+:sh1) a -> Array sh1 a
+   Array (sh0::+sh1) a -> Array sh1 a
 takeRight =
-   takeCenter . mapShape (\(sh0 :+: sh1) -> (sh0 :+: sh1 :+: Shape.Zero))
+   takeCenter . mapShape (\(sh0 ::+ sh1) -> (sh0 ::+ sh1 ::+ Shape.Zero))
 
 split ::
    (Shape.C sh0, Shape.C sh1, Storable a) =>
-   Array (sh0:+:sh1) a -> (Array sh0 a, Array sh1 a)
+   Array (sh0::+sh1) a -> (Array sh0 a, Array sh1 a)
 split x = (takeLeft x, takeRight x)
 
 {- |
@@ -188,8 +191,8 @@
 -}
 takeCenter ::
    (Shape.C sh0, Shape.C sh1, Shape.C sh2, Storable a) =>
-   Array (sh0:+:sh1:+:sh2) a -> Array sh1 a
-takeCenter (Array (sh0:+:sh1:+:_sh2) x) =
+   Array (sh0::+sh1::+sh2) a -> Array sh1 a
+takeCenter (Array (sh0::+sh1::+_sh2) x) =
    unsafeCreateWithSize sh1 $ \k yPtr ->
    withForeignPtr x $ \xPtr ->
       copyArray yPtr (advancePtr xPtr (Shape.size sh0)) k
diff --git a/test/DocTest/Data/Array/Comfort/Shape.hs b/test/DocTest/Data/Array/Comfort/Shape.hs
--- a/test/DocTest/Data/Array/Comfort/Shape.hs
+++ b/test/DocTest/Data/Array/Comfort/Shape.hs
@@ -1,77 +1,77 @@
 -- Do not edit! Automatically created with doctest-extract from src/Data/Array/Comfort/Shape.hs
-{-# LINE 65 "src/Data/Array/Comfort/Shape.hs" #-}
+{-# LINE 95 "src/Data/Array/Comfort/Shape.hs" #-}
 
 module DocTest.Data.Array.Comfort.Shape where
 
 import Test.DocTest.Base
 import qualified Test.DocTest.Driver as DocTest
 
-{-# LINE 66 "src/Data/Array/Comfort/Shape.hs" #-}
+{-# LINE 96 "src/Data/Array/Comfort/Shape.hs" #-}
 import     qualified Data.Array.Comfort.Shape as Shape
 import     qualified Data.Map as Map
 import     qualified Data.Set as Set
-import     Data.Array.Comfort.Shape ((:+:)((:+:)))
+import     Data.Array.Comfort.Shape ((::+)((::+)))
 
 test :: DocTest.T ()
 test = do
- DocTest.printPrefix "Data.Array.Comfort.Shape:128: "
-{-# LINE 128 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Shape:264: "
+{-# LINE 264 "src/Data/Array/Comfort/Shape.hs" #-}
  DocTest.example
-{-# LINE 128 "src/Data/Array/Comfort/Shape.hs" #-}
+{-# LINE 264 "src/Data/Array/Comfort/Shape.hs" #-}
    (Shape.indices ())
   [ExpectedLine [LineChunk "[()]"]]
- DocTest.printPrefix "Data.Array.Comfort.Shape:150: "
-{-# LINE 150 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Shape:283: "
+{-# LINE 283 "src/Data/Array/Comfort/Shape.hs" #-}
  DocTest.example
-{-# LINE 150 "src/Data/Array/Comfort/Shape.hs" #-}
+{-# LINE 283 "src/Data/Array/Comfort/Shape.hs" #-}
    (Shape.indices (Shape.ZeroBased (7::Int)))
   [ExpectedLine [LineChunk "[0,1,2,3,4,5,6]"]]
- DocTest.printPrefix "Data.Array.Comfort.Shape:201: "
-{-# LINE 201 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Shape:330: "
+{-# LINE 330 "src/Data/Array/Comfort/Shape.hs" #-}
  DocTest.example
-{-# LINE 201 "src/Data/Array/Comfort/Shape.hs" #-}
+{-# LINE 330 "src/Data/Array/Comfort/Shape.hs" #-}
    (Shape.indices (Shape.OneBased (7::Int)))
   [ExpectedLine [LineChunk "[1,2,3,4,5,6,7]"]]
- DocTest.printPrefix "Data.Array.Comfort.Shape:249: "
-{-# LINE 249 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Shape:374: "
+{-# LINE 374 "src/Data/Array/Comfort/Shape.hs" #-}
  DocTest.example
-{-# LINE 249 "src/Data/Array/Comfort/Shape.hs" #-}
+{-# LINE 374 "src/Data/Array/Comfort/Shape.hs" #-}
    (Shape.indices (Shape.Range (-5) (5::Int)))
   [ExpectedLine [LineChunk "[-5,-4,-3,-2,-1,0,1,2,3,4,5]"]]
- DocTest.printPrefix "Data.Array.Comfort.Shape:251: "
-{-# LINE 251 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Shape:376: "
+{-# LINE 376 "src/Data/Array/Comfort/Shape.hs" #-}
  DocTest.example
-{-# LINE 251 "src/Data/Array/Comfort/Shape.hs" #-}
+{-# LINE 376 "src/Data/Array/Comfort/Shape.hs" #-}
    (Shape.indices (Shape.Range (-1,-1) (1::Int,1::Int)))
   [ExpectedLine [LineChunk "[(-1,-1),(-1,0),(-1,1),(0,-1),(0,0),(0,1),(1,-1),(1,0),(1,1)]"]]
- DocTest.printPrefix "Data.Array.Comfort.Shape:301: "
-{-# LINE 301 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Shape:426: "
+{-# LINE 426 "src/Data/Array/Comfort/Shape.hs" #-}
  DocTest.example
-{-# LINE 301 "src/Data/Array/Comfort/Shape.hs" #-}
+{-# LINE 426 "src/Data/Array/Comfort/Shape.hs" #-}
    (Shape.indices (Shape.Shifted (-4) (8::Int)))
   [ExpectedLine [LineChunk "[-4,-3,-2,-1,0,1,2,3]"]]
- DocTest.printPrefix "Data.Array.Comfort.Shape:373: "
-{-# LINE 373 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Shape:491: "
+{-# LINE 491 "src/Data/Array/Comfort/Shape.hs" #-}
  DocTest.example
-{-# LINE 373 "src/Data/Array/Comfort/Shape.hs" #-}
+{-# LINE 491 "src/Data/Array/Comfort/Shape.hs" #-}
    (Shape.indices (Shape.Enumeration :: Shape.Enumeration Ordering))
   [ExpectedLine [LineChunk "[LT,EQ,GT]"]]
- DocTest.printPrefix "Data.Array.Comfort.Shape:428: "
-{-# LINE 428 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Shape:548: "
+{-# LINE 548 "src/Data/Array/Comfort/Shape.hs" #-}
  DocTest.example
-{-# LINE 428 "src/Data/Array/Comfort/Shape.hs" #-}
+{-# LINE 548 "src/Data/Array/Comfort/Shape.hs" #-}
    (Shape.indices (Set.fromList "comfort"))
   [ExpectedLine [LineChunk "\"cfmort\""]]
- DocTest.printPrefix "Data.Array.Comfort.Shape:454: "
-{-# LINE 454 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Shape:584: "
+{-# LINE 584 "src/Data/Array/Comfort/Shape.hs" #-}
  DocTest.example
-{-# LINE 454 "src/Data/Array/Comfort/Shape.hs" #-}
+{-# LINE 584 "src/Data/Array/Comfort/Shape.hs" #-}
    (Shape.indices $ fmap Shape.ZeroBased $ Map.fromList [('b', (0::Int)), ('a', 5), ('c', 2)])
   [ExpectedLine [LineChunk "[('a',0),('a',1),('a',2),('a',3),('a',4),('c',0),('c',1)]"]]
- DocTest.printPrefix "Data.Array.Comfort.Shape:528: "
-{-# LINE 528 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Shape:649: "
+{-# LINE 649 "src/Data/Array/Comfort/Shape.hs" #-}
  DocTest.example
-{-# LINE 528 "src/Data/Array/Comfort/Shape.hs" #-}
+{-# LINE 649 "src/Data/Array/Comfort/Shape.hs" #-}
    (
   let sh2 = (Shape.ZeroBased (2::Int), Shape.ZeroBased (2::Int)) in
   let sh3 = (Shape.ZeroBased (3::Int), Shape.ZeroBased (3::Int)) in
@@ -80,50 +80,98 @@
   Shape.indexFromOffset (Shape.Deferred sh2) 3)
   )
   [ExpectedLine [LineChunk "(4,3)"]]
- DocTest.printPrefix "Data.Array.Comfort.Shape:630: "
-{-# LINE 630 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Shape:749: "
+{-# LINE 749 "src/Data/Array/Comfort/Shape.hs" #-}
  DocTest.example
-{-# LINE 630 "src/Data/Array/Comfort/Shape.hs" #-}
+{-# LINE 749 "src/Data/Array/Comfort/Shape.hs" #-}
    (Shape.indices (Shape.ZeroBased (3::Int), Shape.ZeroBased (3::Int)))
   [ExpectedLine [LineChunk "[(0,0),(0,1),(0,2),(1,0),(1,1),(1,2),(2,0),(2,1),(2,2)]"]]
- DocTest.printPrefix "Data.Array.Comfort.Shape:744: "
-{-# LINE 744 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Shape:849: "
+{-# LINE 849 "src/Data/Array/Comfort/Shape.hs" #-}
  DocTest.example
-{-# LINE 744 "src/Data/Array/Comfort/Shape.hs" #-}
+{-# LINE 849 "src/Data/Array/Comfort/Shape.hs" #-}
    (Shape.indices $ Shape.Square $ Shape.ZeroBased (3::Int))
   [ExpectedLine [LineChunk "[(0,0),(0,1),(0,2),(1,0),(1,1),(1,2),(2,0),(2,1),(2,2)]"]]
- DocTest.printPrefix "Data.Array.Comfort.Shape:793: "
-{-# LINE 793 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Shape:892: "
+{-# LINE 892 "src/Data/Array/Comfort/Shape.hs" #-}
  DocTest.example
-{-# LINE 793 "src/Data/Array/Comfort/Shape.hs" #-}
+{-# LINE 892 "src/Data/Array/Comfort/Shape.hs" #-}
    (Shape.indices $ Shape.Cube $ Shape.ZeroBased (2::Int))
   [ExpectedLine [LineChunk "[(0,0,0),(0,0,1),(0,1,0),(0,1,1),(1,0,0),(1,0,1),(1,1,0),(1,1,1)]"]]
- DocTest.printPrefix "Data.Array.Comfort.Shape:862: "
-{-# LINE 862 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Shape:951: "
+{-# LINE 951 "src/Data/Array/Comfort/Shape.hs" #-}
  DocTest.example
-{-# LINE 862 "src/Data/Array/Comfort/Shape.hs" #-}
+{-# LINE 951 "src/Data/Array/Comfort/Shape.hs" #-}
    (Shape.indices $ Shape.Triangular Shape.Upper $ Shape.ZeroBased (3::Int))
   [ExpectedLine [LineChunk "[(0,0),(0,1),(0,2),(1,1),(1,2),(2,2)]"]]
- DocTest.printPrefix "Data.Array.Comfort.Shape:864: "
-{-# LINE 864 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Shape:953: "
+{-# LINE 953 "src/Data/Array/Comfort/Shape.hs" #-}
  DocTest.example
-{-# LINE 864 "src/Data/Array/Comfort/Shape.hs" #-}
+{-# LINE 953 "src/Data/Array/Comfort/Shape.hs" #-}
    (Shape.indices $ Shape.Triangular Shape.Lower $ Shape.ZeroBased (3::Int))
   [ExpectedLine [LineChunk "[(0,0),(1,0),(1,1),(2,0),(2,1),(2,2)]"]]
- DocTest.printPrefix "Data.Array.Comfort.Shape:994: "
-{-# LINE 994 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Shape:1085: "
+{-# LINE 1085 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.example
+{-# LINE 1085 "src/Data/Array/Comfort/Shape.hs" #-}
+   (Shape.indices $ Shape.simplexAscending (replicate 3 Shape.AllDistinct) $ Shape.ZeroBased (4::Int))
+  [ExpectedLine [LineChunk "[[0,1,2],[0,1,3],[0,2,3],[1,2,3]]"]]
+ DocTest.printPrefix "Data.Array.Comfort.Shape:1087: "
+{-# LINE 1087 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.example
+{-# LINE 1087 "src/Data/Array/Comfort/Shape.hs" #-}
+   (Shape.indices $ Shape.simplexAscending (replicate 3 Shape.SomeRepetitive) $ Shape.ZeroBased (3::Int))
+  [ExpectedLine [LineChunk "[[0,0,0],[0,0,1],[0,0,2],[0,1,1],[0,1,2],[0,2,2],[1,1,1],[1,1,2],[1,2,2],[2,2,2]]"]]
+ DocTest.printPrefix "Data.Array.Comfort.Shape:1089: "
+{-# LINE 1089 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.example
+{-# LINE 1089 "src/Data/Array/Comfort/Shape.hs" #-}
+   (Shape.indices $ Shape.simplexAscending [Shape.Repetitive,Shape.Distinct,Shape.Repetitive] $ Shape.ZeroBased (4::Int))
+  [ExpectedLine [LineChunk "[[0,0,1],[0,0,2],[0,0,3],[0,1,2],[0,1,3],[0,2,3],[1,1,2],[1,1,3],[1,2,3],[2,2,3]]"]]
+ DocTest.printPrefix "Data.Array.Comfort.Shape:1091: "
+{-# LINE 1091 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.example
+{-# LINE 1091 "src/Data/Array/Comfort/Shape.hs" #-}
+   (Shape.indices $ Shape.simplexAscending [Shape.Repetitive,Shape.Distinct,Shape.Distinct] $ Shape.ZeroBased (4::Int))
+  [ExpectedLine [LineChunk "[[0,0,1],[0,0,2],[0,0,3],[0,1,2],[0,1,3],[0,2,3],[1,1,2],[1,1,3],[1,2,3],[2,2,3]]"]]
+ DocTest.printPrefix "Data.Array.Comfort.Shape:1094: "
+{-# LINE 1094 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.example
+{-# LINE 1094 "src/Data/Array/Comfort/Shape.hs" #-}
+   (Shape.indices $ Shape.simplexDescending (replicate 3 Shape.AllDistinct) $ Shape.ZeroBased (4::Int))
+  [ExpectedLine [LineChunk "[[2,1,0],[3,1,0],[3,2,0],[3,2,1]]"]]
+ DocTest.printPrefix "Data.Array.Comfort.Shape:1096: "
+{-# LINE 1096 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.example
+{-# LINE 1096 "src/Data/Array/Comfort/Shape.hs" #-}
+   (Shape.indices $ Shape.simplexDescending (replicate 3 Shape.SomeRepetitive) $ Shape.ZeroBased (3::Int))
+  [ExpectedLine [LineChunk "[[0,0,0],[1,0,0],[1,1,0],[1,1,1],[2,0,0],[2,1,0],[2,1,1],[2,2,0],[2,2,1],[2,2,2]]"]]
+ DocTest.printPrefix "Data.Array.Comfort.Shape:1098: "
+{-# LINE 1098 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.example
+{-# LINE 1098 "src/Data/Array/Comfort/Shape.hs" #-}
+   (Shape.indices $ Shape.simplexDescending [Shape.Repetitive,Shape.Distinct,Shape.Repetitive] $ Shape.ZeroBased (4::Int))
+  [ExpectedLine [LineChunk "[[1,1,0],[2,1,0],[2,2,0],[2,2,1],[3,1,0],[3,2,0],[3,2,1],[3,3,0],[3,3,1],[3,3,2]]"]]
+ DocTest.printPrefix "Data.Array.Comfort.Shape:1100: "
+{-# LINE 1100 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.example
+{-# LINE 1100 "src/Data/Array/Comfort/Shape.hs" #-}
+   (Shape.indices $ Shape.simplexDescending [Shape.Repetitive,Shape.Distinct,Shape.Distinct] $ Shape.ZeroBased (4::Int))
+  [ExpectedLine [LineChunk "[[1,1,0],[2,1,0],[2,2,0],[2,2,1],[3,1,0],[3,2,0],[3,2,1],[3,3,0],[3,3,1],[3,3,2]]"]]
+ DocTest.printPrefix "Data.Array.Comfort.Shape:1334: "
+{-# LINE 1334 "src/Data/Array/Comfort/Shape.hs" #-}
  DocTest.property
-{-# LINE 994 "src/Data/Array/Comfort/Shape.hs" #-}
+{-# LINE 1334 "src/Data/Array/Comfort/Shape.hs" #-}
      (let shape = Shape.Cyclic (10::Int) in Shape.offset shape (-1) == Shape.offset shape 9)
- DocTest.printPrefix "Data.Array.Comfort.Shape:999: "
-{-# LINE 999 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Shape:1339: "
+{-# LINE 1339 "src/Data/Array/Comfort/Shape.hs" #-}
  DocTest.example
-{-# LINE 999 "src/Data/Array/Comfort/Shape.hs" #-}
+{-# LINE 1339 "src/Data/Array/Comfort/Shape.hs" #-}
    (Shape.indices (Shape.Cyclic (7::Int)))
   [ExpectedLine [LineChunk "[0,1,2,3,4,5,6]"]]
- DocTest.printPrefix "Data.Array.Comfort.Shape:1047: "
-{-# LINE 1047 "src/Data/Array/Comfort/Shape.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Shape:1383: "
+{-# LINE 1383 "src/Data/Array/Comfort/Shape.hs" #-}
  DocTest.example
-{-# LINE 1047 "src/Data/Array/Comfort/Shape.hs" #-}
-   (Shape.indices (Shape.ZeroBased (3::Int) :+: Shape.Range 'a' 'c'))
+{-# LINE 1383 "src/Data/Array/Comfort/Shape.hs" #-}
+   (Shape.indices (Shape.ZeroBased (3::Int) ::+ Shape.Range 'a' 'c'))
   [ExpectedLine [LineChunk "[Left 0,Left 1,Left 2,Right 'a',Right 'b',Right 'c']"]]
diff --git a/test/DocTest/Data/Array/Comfort/Storable.hs b/test/DocTest/Data/Array/Comfort/Storable.hs
--- a/test/DocTest/Data/Array/Comfort/Storable.hs
+++ b/test/DocTest/Data/Array/Comfort/Storable.hs
@@ -1,19 +1,19 @@
 -- Do not edit! Automatically created with doctest-extract from src/Data/Array/Comfort/Storable.hs
-{-# LINE 73 "src/Data/Array/Comfort/Storable.hs" #-}
+{-# LINE 83 "src/Data/Array/Comfort/Storable.hs" #-}
 
 module DocTest.Data.Array.Comfort.Storable where
 
 import qualified Test.DocTest.Driver as DocTest
 
-{-# LINE 74 "src/Data/Array/Comfort/Storable.hs" #-}
+{-# LINE 84 "src/Data/Array/Comfort/Storable.hs" #-}
 import     qualified Data.Array.Comfort.Storable as Array
 import     qualified Data.Array.Comfort.Shape as Shape
-import     Data.Array.Comfort.Storable (Array)
+import     Data.Array.Comfort.Storable (Array, (!))
 
 import     qualified Test.QuickCheck as QC
 import     Test.ChasingBottoms.IsBottom (isBottom)
 
-import     Control.Applicative ((<$>))
+import     Control.Applicative ((<$>), (<*>))
 
 import     Data.Word (Word16)
 
@@ -22,6 +22,23 @@
 genArray     :: QC.Gen (Array ShapeInt Word16)
 genArray     = Array.vectorFromList <$> QC.arbitrary
 
+genArray2     :: QC.Gen (Array (ShapeInt,ShapeInt) Word16)
+genArray2     = do
+       xs <- QC.arbitrary
+       let n = length xs
+       (k,m) <-
+          if n == 0
+             then QC.elements [(,) 0, flip (,) 0] <*> QC.choose (1,n)
+             else fmap (\m -> (div n m, m)) $ QC.choose (1,n)
+       return $ Array.fromList (Shape.ZeroBased k, Shape.ZeroBased m) xs
+
+genNonEmptyArray2     :: QC.Gen (Array (ShapeInt,ShapeInt) Word16)
+genNonEmptyArray2     = do
+       xs <- QC.getNonEmpty <$> QC.arbitrary
+       let n = length xs
+       m <- QC.choose (1,n)
+       return $ Array.fromList (Shape.ZeroBased (div n m), Shape.ZeroBased m) xs
+
 infix     4 ==?
 (==?)     :: a -> a -> (a,a)
 (==?)     = (,)
@@ -37,18 +54,28 @@
 
 test :: DocTest.T ()
 test = do
- DocTest.printPrefix "Data.Array.Comfort.Storable:194: "
-{-# LINE 194 "src/Data/Array/Comfort/Storable.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Storable:243: "
+{-# LINE 243 "src/Data/Array/Comfort/Storable.hs" #-}
  DocTest.property
-{-# LINE 194 "src/Data/Array/Comfort/Storable.hs" #-}
+{-# LINE 243 "src/Data/Array/Comfort/Storable.hs" #-}
+     (QC.forAll genNonEmptyArray2 $ \xs -> QC.forAll (QC.elements $ Shape.indices $ Array.shape xs) $ \(ix0,ix1) -> Array.pick xs ix0 ! ix1 == xs!(ix0,ix1))
+ DocTest.printPrefix "Data.Array.Comfort.Storable:261: "
+{-# LINE 261 "src/Data/Array/Comfort/Storable.hs" #-}
+ DocTest.property
+{-# LINE 261 "src/Data/Array/Comfort/Storable.hs" #-}
+     (QC.forAll genArray2 $ \xs -> xs == Array.fromRowArray (snd $ Array.shape xs) (Array.toRowArray xs))
+ DocTest.printPrefix "Data.Array.Comfort.Storable:278: "
+{-# LINE 278 "src/Data/Array/Comfort/Storable.hs" #-}
+ DocTest.property
+{-# LINE 278 "src/Data/Array/Comfort/Storable.hs" #-}
      (forAllNonEmpty $ \xs -> Array.minimum xs ==? minimum (Array.toList xs))
- DocTest.printPrefix "Data.Array.Comfort.Storable:202: "
-{-# LINE 202 "src/Data/Array/Comfort/Storable.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Storable:286: "
+{-# LINE 286 "src/Data/Array/Comfort/Storable.hs" #-}
  DocTest.property
-{-# LINE 202 "src/Data/Array/Comfort/Storable.hs" #-}
+{-# LINE 286 "src/Data/Array/Comfort/Storable.hs" #-}
      (forAllNonEmpty $ \xs -> Array.maximum xs ==? maximum (Array.toList xs))
- DocTest.printPrefix "Data.Array.Comfort.Storable:214: "
-{-# LINE 214 "src/Data/Array/Comfort/Storable.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Storable:298: "
+{-# LINE 298 "src/Data/Array/Comfort/Storable.hs" #-}
  DocTest.property
-{-# LINE 214 "src/Data/Array/Comfort/Storable.hs" #-}
+{-# LINE 298 "src/Data/Array/Comfort/Storable.hs" #-}
      (forAllNonEmpty $ \xs -> Array.limits xs ==? (Array.minimum xs, Array.maximum xs))
diff --git a/test/DocTest/Data/Array/Comfort/Storable/Unchecked.hs b/test/DocTest/Data/Array/Comfort/Storable/Unchecked.hs
--- a/test/DocTest/Data/Array/Comfort/Storable/Unchecked.hs
+++ b/test/DocTest/Data/Array/Comfort/Storable/Unchecked.hs
@@ -29,28 +29,28 @@
  DocTest.property
 {-# LINE 134 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
      (\x  ->  Array.singleton x ! () == (x::Word16))
- DocTest.printPrefix "Data.Array.Comfort.Storable.Unchecked:152: "
-{-# LINE 152 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Storable.Unchecked:155: "
+{-# LINE 155 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
  DocTest.property
-{-# LINE 152 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
+{-# LINE 155 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
      (\(QC.NonNegative n) (Array16 x)  ->  x == Array.mapShape (Shape.ZeroBased . Shape.size) (Array.append (Array.take n x) (Array.drop n x)))
- DocTest.printPrefix "Data.Array.Comfort.Storable.Unchecked:167: "
-{-# LINE 167 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Storable.Unchecked:170: "
+{-# LINE 170 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
  DocTest.property
-{-# LINE 167 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
+{-# LINE 170 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
      (\(Array16 x) (Array16 y) -> let xy = Array.append x y in x == Array.takeLeft xy  &&  y == Array.takeRight xy)
- DocTest.printPrefix "Data.Array.Comfort.Storable.Unchecked:187: "
-{-# LINE 187 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Storable.Unchecked:190: "
+{-# LINE 190 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
  DocTest.property
-{-# LINE 187 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
+{-# LINE 190 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
      (\(Array16 x) (Array16 y) (Array16 z) -> let xyz = Array.append x $ Array.append y z in y == Array.takeCenter xyz)
- DocTest.printPrefix "Data.Array.Comfort.Storable.Unchecked:200: "
-{-# LINE 200 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Storable.Unchecked:203: "
+{-# LINE 203 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
  DocTest.property
-{-# LINE 200 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
+{-# LINE 203 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
      (\(Array16 xs)  ->  Array.sum xs == sum (Array.toList xs))
- DocTest.printPrefix "Data.Array.Comfort.Storable.Unchecked:206: "
-{-# LINE 206 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
+ DocTest.printPrefix "Data.Array.Comfort.Storable.Unchecked:209: "
+{-# LINE 209 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
  DocTest.property
-{-# LINE 206 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
+{-# LINE 209 "src/Data/Array/Comfort/Storable/Unchecked.hs" #-}
      (\(Array16 xs)  ->  Array.product xs == product (Array.toList xs))
diff --git a/test/Main.hs b/test/Main.hs
--- a/test/Main.hs
+++ b/test/Main.hs
@@ -1,4 +1,3 @@
-{-# LANGUAGE TypeFamilies #-}
 module Main where
 
 import qualified DocTest.Main as DocTestMain
diff --git a/test/Test/Shape.hs b/test/Test/Shape.hs
--- a/test/Test/Shape.hs
+++ b/test/Test/Shape.hs
@@ -2,12 +2,12 @@
 
 import qualified Data.Array.Comfort.Shape.Test as ShapeTest
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Shape ((:+:)((:+:)))
+import Data.Array.Comfort.Shape ((::+)((::+)))
 
 import qualified Test.QuickCheck as QC
 import Test.Utility (prefix)
 
-import Control.Applicative (liftA2, liftA3)
+import Control.Applicative (liftA2, liftA3, pure, (<$>))
 
 import qualified Data.Map as Map
 import qualified Data.Set as Set
@@ -20,6 +20,27 @@
 tag :: sh -> Tagged Ordering sh
 tag = Tagged
 
+simplex ::
+   (Shape.SimplexOrderC order) =>
+   Shape.SimplexOrder order -> [(String, QC.Property)]
+simplex order =
+   prefix "Mixed"
+      (ShapeTest.tests $
+       liftA2 (Shape.Simplex order)
+         (take 4 <$> QC.listOf (QC.elements [Shape.Distinct,Shape.Repetitive]))
+         (genZeroBased 10)) ++
+   prefix "Strict"
+      (ShapeTest.tests $
+       liftA2 (Shape.Simplex order)
+         (take 4 <$> QC.listOf (pure Shape.AllDistinct))
+         (genZeroBased 10)) ++
+   prefix "Weak"
+      (ShapeTest.tests $
+       liftA2 (Shape.Simplex order)
+         (take 4 <$> QC.listOf (pure Shape.SomeRepetitive))
+         (genZeroBased 10)) ++
+   []
+
 tests :: [(String, QC.Property)]
 tests =
    prefix "ZeroBased"
@@ -62,7 +83,7 @@
       (ShapeTest.tests $
        liftA3 (,,) (genZeroBased 10) (genZeroBased 10) (genZeroBased 10)) ++
    prefix "Append"
-      (ShapeTest.tests $ liftA2 (:+:) (genZeroBased 10) (genZeroBased 10)) ++
+      (ShapeTest.tests $ liftA2 (::+) (genZeroBased 10) (genZeroBased 10)) ++
    prefix "Square"
       (ShapeTest.tests $ fmap Shape.Square $ genZeroBased 10) ++
    prefix "Cube"
@@ -73,6 +94,9 @@
    prefix "Triangular Upper"
       (ShapeTest.tests $
        fmap (Shape.Triangular Shape.Upper) (genZeroBased 10)) ++
+   prefix "Simplex"
+      (prefix "Upper" (simplex Shape.Ascending) ++
+       prefix "Lower" (simplex Shape.Descending)) ++
    prefix "Cyclic"
       (ShapeTest.tests $ fmap Shape.Cyclic $ QC.choose (0,10::Int)) ++
    []