diff --git a/knead.cabal b/knead.cabal
--- a/knead.cabal
+++ b/knead.cabal
@@ -1,5 +1,5 @@
 Name:             knead
-Version:          0.3
+Version:          0.4
 License:          BSD3
 License-File:     LICENSE
 Author:           Henning Thielemann <haskell@henning-thielemann.de>
@@ -49,7 +49,7 @@
 Build-Type:       Simple
 
 Source-Repository this
-  Tag:         0.3
+  Tag:         0.4
   Type:        darcs
   Location:    http://hub.darcs.net/thielema/knead/
 
@@ -62,6 +62,7 @@
     llvm-extra >=0.8 && <0.9,
     llvm-tf >=3.1.1 && <3.2,
     tfp >=1.0 && <1.1,
+    comfort-array >=0.3 && <0.5,
     storable-tuple >=0.0 && <0.1,
     storable-record >=0.0.3 && <0.1,
     bool8 >=0.0 && <0.1,
@@ -72,7 +73,7 @@
   GHC-Options:      -Wall
   Hs-Source-Dirs:   src
   Exposed-Modules:
-    Data.Array.Knead.Shape.Nested
+    Data.Array.Knead.Shape
     Data.Array.Knead.Shape.Cubic
     Data.Array.Knead.Shape.Cubic.Int
     Data.Array.Knead.Expression
diff --git a/src/Data/Array/Knead/Code.hs b/src/Data/Array/Knead/Code.hs
--- a/src/Data/Array/Knead/Code.hs
+++ b/src/Data/Array/Knead/Code.hs
@@ -1,7 +1,7 @@
 {-# LANGUAGE TypeFamilies #-}
 module Data.Array.Knead.Code where
 
-import qualified Data.Array.Knead.Shape.Nested as Shape
+import qualified Data.Array.Knead.Shape as Shape
 
 import qualified LLVM.Extra.Multi.Value as MultiValue
 
@@ -23,7 +23,7 @@
    MultiValue.T ix ->
    LLVM.CodeGenFunction r (LLVM.Value (Ptr a))
 getElementPtr sh ptr ix = do
-   n <- Shape.flattenIndex sh ix
+   n <- Shape.offset sh ix
    LLVM.getElementPtr ptr (n, ())
 
 
diff --git a/src/Data/Array/Knead/Expression.hs b/src/Data/Array/Knead/Expression.hs
--- a/src/Data/Array/Knead/Expression.hs
+++ b/src/Data/Array/Knead/Expression.hs
@@ -244,7 +244,8 @@
    patternA ->
    patternB ->
    (Decomposed Exp patternA -> Decomposed Exp patternB -> a) ->
-   Exp (PatternTuple patternA) -> Exp (PatternTuple patternB) -> Exp (Composed a)
+   Exp (PatternTuple patternA) ->
+   Exp (PatternTuple patternB) -> Exp (Composed a)
 modify2 pa pb f a b = compose $ f (decompose pa a) (decompose pb b)
 
 
@@ -277,7 +278,8 @@
 
 instance (Compose a, Compose b, Compose c) => Compose (a,b,c) where
    type Composed (a,b,c) = (Composed a, Composed b, Composed c)
-   compose = TupleHT.uncurry3 zip3 . TupleHT.mapTriple (compose, compose, compose)
+   compose =
+      TupleHT.uncurry3 zip3 . TupleHT.mapTriple (compose, compose, compose)
 
 instance
    (Decompose pa, Decompose pb, Decompose pc) =>
diff --git a/src/Data/Array/Knead/Parameterized/Physical.hs b/src/Data/Array/Knead/Parameterized/Physical.hs
--- a/src/Data/Array/Knead/Parameterized/Physical.hs
+++ b/src/Data/Array/Knead/Parameterized/Physical.hs
@@ -2,7 +2,7 @@
 {-# LANGUAGE ForeignFunctionInterface #-}
 module Data.Array.Knead.Parameterized.Physical (
    Phys.Array,
-   Phys.shape,
+   Array.shape,
    Phys.fromList,
    feed,
    the,
@@ -20,12 +20,14 @@
 import qualified Data.Array.Knead.Simple.Physical as Phys
 import qualified Data.Array.Knead.Simple.Private as Core
 import qualified Data.Array.Knead.Parameter as Param
-import qualified Data.Array.Knead.Shape.Nested as Shape
+import qualified Data.Array.Knead.Shape as Shape
 import qualified Data.Array.Knead.Expression as Expr
 import qualified Data.Array.Knead.Code as Code
 import Data.Array.Knead.Expression (Exp, unExp, )
 import Data.Array.Knead.Code (getElementPtr, compile, )
 
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
+
 import qualified LLVM.Extra.Multi.Value.Memory as MultiValueMemory
 import qualified LLVM.Extra.Multi.Value as MultiValue
 import qualified LLVM.Extra.Memory as Memory
@@ -50,7 +52,7 @@
     MultiValueMemory.C a) =>
    Param.T p (Phys.Array sh a) -> Sym.Array p sh a
 feed arr =
-   Param.withMulti (fmap Phys.shape arr) $ \getShape valueShape ->
+   Param.withMulti (fmap Array.shape arr) $ \getShape valueShape ->
    Sym.Array
       (\p ->
          case mapFst valueShape $ MultiValue.unzip p of
@@ -58,7 +60,7 @@
                Core.Array (Expr.lift0 sh) $
                   Memory.load <=< getElementPtr sh ptr)
       (\p ->
-         case Phys.buffer $ Param.get arr p of
+         case Array.buffer $ Param.get arr p of
             fptr ->
                withForeignPtr fptr $ \ptr ->
                   return (fptr, (getShape p, MultiValueMemory.castStructPtr ptr)))
@@ -110,7 +112,7 @@
            Core.Array esh _code -> do
               sh <- unExp esh
               MultiValueMemory.store sh resultPtr
-              Shape.sizeCode sh >>= LLVM.ret
+              Shape.size sh >>= LLVM.ret
    return $ \p ->
       bracket (create p) (delete . fst) $ \(_ctx, param) ->
       alloca $ \shptr ->
diff --git a/src/Data/Array/Knead/Parameterized/PhysicalHull.hs b/src/Data/Array/Knead/Parameterized/PhysicalHull.hs
--- a/src/Data/Array/Knead/Parameterized/PhysicalHull.hs
+++ b/src/Data/Array/Knead/Parameterized/PhysicalHull.hs
@@ -24,14 +24,15 @@
 
 import qualified Data.Array.Knead.Parameterized.Private as Sym
 import qualified Data.Array.Knead.Simple.PhysicalPrivate as Priv
-import qualified Data.Array.Knead.Simple.Physical as Phys
 import qualified Data.Array.Knead.Simple.Private as Core
-import qualified Data.Array.Knead.Shape.Nested as Shape
+import qualified Data.Array.Knead.Shape as Shape
 import qualified Data.Array.Knead.Expression as Expr
 import qualified Data.Array.Knead.Code as Code
 import Data.Array.Knead.Expression (Exp, unExp, )
 import Data.Array.Knead.Code (compile, )
 
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
+
 import qualified LLVM.Extra.Multi.Value.Memory as MultiValueMemory
 import qualified LLVM.Extra.Multi.Value as MultiValue
 import qualified LLVM.Extra.Memory as Memory
@@ -84,7 +85,7 @@
     LLVM.Value (Ptr (MultiValueMemory.Struct sh)) ->
     LLVM.Value (Ptr (MultiValueMemory.Struct a)) ->
     LLVM.CodeGenFunction () ()) ->
-   Sym.Hull p core -> IO (p -> IO (Phys.Array sh a))
+   Sym.Hull p core -> IO (p -> IO (Array sh a))
 materialize name shape fill (Sym.Hull core create delete) = do
    (fsh, farr) <-
       compile name $
@@ -94,7 +95,7 @@
             param <- Memory.load paramPtr
             sh <- unExp $ shape $ core param
             MultiValueMemory.store sh resultPtr
-            Shape.sizeCode sh >>= LLVM.ret)
+            Shape.size sh >>= LLVM.ret)
          (Code.createFunction callFill "fill" $
           \paramPtr shapePtr bufferPtr -> do
             param <- Memory.load paramPtr
@@ -110,7 +111,7 @@
          fptr <- mallocArray =<< fsh paramMVPtr shapeMVPtr
          withForeignMemPtr fptr $ farr paramMVPtr shapeMVPtr
          sh <- peek shptr
-         return (Phys.Array sh fptr)
+         return (Array sh fptr)
 
 
 foreign import ccall safe "dynamic" callFillExpArray ::
@@ -128,7 +129,7 @@
     LLVM.Value (Ptr (MultiValueMemory.Struct sh)) ->
     LLVM.Value (Ptr (MultiValueMemory.Struct a)) ->
     LLVM.CodeGenFunction () ()) ->
-   Sym.Hull p core -> IO (p -> IO (b, Phys.Array sh a))
+   Sym.Hull p core -> IO (p -> IO (b, Array sh a))
 materializeExpArray name shape fill (Sym.Hull core create delete) = do
    (fsh, farr) <-
       compile name $
@@ -138,7 +139,7 @@
             param <- Memory.load paramPtr
             sh <- unExp $ shape $ core param
             MultiValueMemory.store sh resultPtr
-            Shape.sizeCode sh >>= LLVM.ret)
+            Shape.size sh >>= LLVM.ret)
          (Code.createFunction callFillExpArray "fill" $
           \paramPtr finalPtr shapePtr bufferPtr -> do
             param <- Memory.load paramPtr
@@ -157,7 +158,7 @@
          withForeignMemPtr fptr $ farr paramMVPtr finalMVPtr shapeMVPtr
          sh <- peek shptr
          final <- peek finalPtr
-         return (final, Phys.Array sh fptr)
+         return (final, Array sh fptr)
 
 
 foreign import ccall safe "dynamic" callShaper2 ::
@@ -180,7 +181,7 @@
     (LLVM.Value (Ptr (MultiValueMemory.Struct shb)),
      LLVM.Value (Ptr (MultiValueMemory.Struct b))) ->
     LLVM.CodeGenFunction () ()) ->
-   Sym.Hull p core -> IO (p -> IO (Phys.Array sha a, Phys.Array shb b))
+   Sym.Hull p core -> IO (p -> IO (Array sha a, Array shb b))
 materialize2 name shape fill (Sym.Hull core create delete) = do
    (fsh, farr) <-
       compile name $
@@ -192,9 +193,9 @@
             MultiValueMemory.store sha shapeAPtr
             MultiValueMemory.store shb shapeBPtr
             sizeAPtr <- LLVM.bitcast sizesPtr
-            flip LLVM.store sizeAPtr =<< Shape.sizeCode sha
+            flip LLVM.store sizeAPtr =<< Shape.size sha
             sizeBPtr <- A.advanceArrayElementPtr sizeAPtr
-            flip LLVM.store sizeBPtr =<< Shape.sizeCode shb
+            flip LLVM.store sizeBPtr =<< Shape.size shb
             LLVM.ret ())
          (Code.createFunction callFill2 "fill" $
           \paramPtr shapeAPtr bufferAPtr shapeBPtr bufferBPtr -> do
@@ -219,14 +220,14 @@
             farr paramMVPtr shapeAMVPtr aptr shapeBMVPtr bptr
          sha <- peek shaPtr
          shb <- peek shbPtr
-         return (Phys.Array sha afptr, Phys.Array shb bfptr)
+         return (Array sha afptr, Array shb bfptr)
 
 
 render ::
    (Shape.C sh, Shape.Index sh ~ ix,
     Storable sh, MultiValueMemory.C sh,
     Storable a, MultiValueMemory.C a) =>
-   Sym.Hull p (Core.Array sh a) -> IO (p -> IO (Phys.Array sh a))
+   Sym.Hull p (Core.Array sh a) -> IO (p -> IO (Array sh a))
 render =
    materialize "render" Core.shape
       (\(Core.Array esh code) shapePtr bufferPtr -> do
@@ -249,7 +250,7 @@
     Shape.C sh1, Shape.Index sh1 ~ ix1,
     Storable sh1, MultiValueMemory.C sh1,
     Storable a, MultiValueMemory.C a) =>
-   Sym.Hull p (Scatter sh0 sh1 a) -> IO (p -> IO (Phys.Array sh1 a))
+   Sym.Hull p (Scatter sh0 sh1 a) -> IO (p -> IO (Array sh1 a))
 scatter =
    materialize "scatter"
       (Core.shape . scatterInit)
@@ -270,7 +271,7 @@
     Shape.C sh1, Shape.Index sh1 ~ ix1,
     Storable sh1, MultiValueMemory.C sh1,
     Storable a, MultiValueMemory.C a) =>
-   Sym.Hull p (ScatterMaybe sh0 sh1 a) -> IO (p -> IO (Phys.Array sh1 a))
+   Sym.Hull p (ScatterMaybe sh0 sh1 a) -> IO (p -> IO (Array sh1 a))
 scatterMaybe =
    materialize "scatterMaybe"
       (Core.shape . scatterMaybeInit)
@@ -291,7 +292,7 @@
     MultiValue.C acc,
     Storable a, MultiValueMemory.C a,
     Storable b, MultiValueMemory.C b) =>
-   Sym.Hull p (MapAccumLSimple sh n acc a b) -> IO (p -> IO (Phys.Array (sh,n) b))
+   Sym.Hull p (MapAccumLSimple sh n acc a b) -> IO (p -> IO (Array (sh,n) b))
 mapAccumLSimple =
    materialize "mapAccumLSimple"
       (Core.shape . mapAccumLSimpleArray)
@@ -315,7 +316,7 @@
     Storable a, MultiValueMemory.C a,
     Storable b, MultiValueMemory.C b) =>
    Sym.Hull p (MapAccumLSequence n acc final a b) ->
-   IO (p -> IO (final, Phys.Array n b))
+   IO (p -> IO (final, Array n b))
 mapAccumLSequence =
    materializeExpArray "mapAccumLSequence"
       (Core.shape . mapAccumLSequenceArray)
@@ -340,7 +341,7 @@
     Storable a, MultiValueMemory.C a,
     Storable b, MultiValueMemory.C b) =>
    Sym.Hull p (MapAccumL sh n acc final a b) ->
-   IO (p -> IO (Phys.Array sh final, Phys.Array (sh,n) b))
+   IO (p -> IO (Array sh final, Array (sh,n) b))
 mapAccumL =
    materialize2 "mapAccumL"
       (\core ->
@@ -363,7 +364,7 @@
    (Shape.C n, Storable n, MultiValueMemory.C n,
     Shape.C sh, Storable sh, MultiValueMemory.C sh,
     Storable a, MultiValueMemory.C a) =>
-   Sym.Hull p (FoldOuterL n sh a b) -> IO (p -> IO (Phys.Array sh a))
+   Sym.Hull p (FoldOuterL n sh a b) -> IO (p -> IO (Array sh a))
 foldOuterL =
    materialize "foldOuterL"
       (Core.shape . foldOuterLInit)
@@ -381,7 +382,7 @@
    (Shape.Sequence n,
     Storable n, MultiValueMemory.C n,
     Storable b, MultiValueMemory.C b) =>
-   Sym.Hull p (MapFilter n a b) -> IO (p -> IO (Phys.Array n b))
+   Sym.Hull p (MapFilter n a b) -> IO (p -> IO (Array n b))
 mapFilter =
    materialize "mapFilter"
       (Core.shape . mapFilterArray)
@@ -401,7 +402,7 @@
    (Shape.Sequence n, Storable n, MultiValueMemory.C n,
     Shape.C sh, Storable sh, MultiValueMemory.C sh,
     Storable a, MultiValueMemory.C a) =>
-   Sym.Hull p (FilterOuter n sh a) -> IO (p -> IO (Phys.Array (n,sh) a))
+   Sym.Hull p (FilterOuter n sh a) -> IO (p -> IO (Array (n,sh) a))
 filterOuter =
    materialize "filterOuter"
       (Core.shape . filterOuterArray)
@@ -421,7 +422,7 @@
    (Shape.C sh, Storable sh, MultiValueMemory.C sh,
     Shape.C n, Storable n, MultiValueMemory.C n,
     Storable b, MultiValueMemory.C b) =>
-   Sym.Hull p (AddDimension sh n a b) -> IO (p -> IO (Phys.Array (sh,n) b))
+   Sym.Hull p (AddDimension sh n a b) -> IO (p -> IO (Array (sh,n) b))
 addDimension =
    materialize "addDimension"
       (\r -> Expr.zip (Core.shape (addDimensionArray r)) (addDimensionSize r))
diff --git a/src/Data/Array/Knead/Parameterized/Private.hs b/src/Data/Array/Knead/Parameterized/Private.hs
--- a/src/Data/Array/Knead/Parameterized/Private.hs
+++ b/src/Data/Array/Knead/Parameterized/Private.hs
@@ -5,7 +5,7 @@
 import qualified Data.Array.Knead.Simple.Symbolic as Core
 
 import qualified Data.Array.Knead.Parameter as Param
-import qualified Data.Array.Knead.Shape.Nested as Shape
+import qualified Data.Array.Knead.Shape as Shape
 import qualified Data.Array.Knead.Expression as Expr
 import Data.Array.Knead.Expression (Exp, )
 
diff --git a/src/Data/Array/Knead/Parameterized/Render.hs b/src/Data/Array/Knead/Parameterized/Render.hs
--- a/src/Data/Array/Knead/Parameterized/Render.hs
+++ b/src/Data/Array/Knead/Parameterized/Render.hs
@@ -21,7 +21,7 @@
 import qualified Data.Array.Knead.Simple.Physical as Phys
 import qualified Data.Array.Knead.Simple.Private as Core
 import qualified Data.Array.Knead.Parameter as Param
-import qualified Data.Array.Knead.Shape.Nested as Shape
+import qualified Data.Array.Knead.Shape as Shape
 import Data.Array.Knead.Parameterized.PhysicalHull
          (MapFilter, FilterOuter,
           MapAccumLSimple, MapAccumLSequence, MapAccumL, FoldOuterL,
diff --git a/src/Data/Array/Knead/Parameterized/Slice.hs b/src/Data/Array/Knead/Parameterized/Slice.hs
--- a/src/Data/Array/Knead/Parameterized/Slice.hs
+++ b/src/Data/Array/Knead/Parameterized/Slice.hs
@@ -18,7 +18,7 @@
 import qualified Data.Array.Knead.Simple.Private as Core
 
 import qualified Data.Array.Knead.Shape.Cubic as Linear
-import qualified Data.Array.Knead.Shape.Nested as Shape
+import qualified Data.Array.Knead.Shape as Shape
 import qualified Data.Array.Knead.Parameter as Param
 import qualified Data.Array.Knead.Expression as Expr
 import Data.Array.Knead.Expression (Exp, )
diff --git a/src/Data/Array/Knead/Parameterized/Symbolic.hs b/src/Data/Array/Knead/Parameterized/Symbolic.hs
--- a/src/Data/Array/Knead/Parameterized/Symbolic.hs
+++ b/src/Data/Array/Knead/Parameterized/Symbolic.hs
@@ -22,7 +22,7 @@
 import qualified Data.Array.Knead.Simple.Symbolic as Core
 import Data.Array.Knead.Parameterized.Private (Array, gather, )
 
-import qualified Data.Array.Knead.Shape.Nested as Shape
+import qualified Data.Array.Knead.Shape as Shape
 import qualified Data.Array.Knead.Parameter as Param
 import qualified Data.Array.Knead.Expression as Expr
 import Data.Array.Knead.Expression (Exp, )
diff --git a/src/Data/Array/Knead/Shape.hs b/src/Data/Array/Knead/Shape.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Knead/Shape.hs
@@ -0,0 +1,498 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE Rank2Types #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+module Data.Array.Knead.Shape (
+   C(..), Index,
+   Size,
+   value,
+   paramWith,
+   load,
+   intersect,
+   offset,
+
+   ZeroBased(ZeroBased), zeroBased, zeroBasedSize,
+
+   Range(Range), range, rangeFrom, rangeTo,
+   Shifted(Shifted), shifted, shiftedOffset, shiftedSize,
+
+   Scalar(..),
+   Sequence(..),
+   ) where
+
+import qualified Data.Array.Knead.Expression as Expr
+import qualified Data.Array.Knead.Parameter as Param
+import Data.Array.Knead.Expression (Exp, )
+
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Shape
+         (Index, ZeroBased(ZeroBased), Range(Range), Shifted(Shifted))
+import Data.Ix (Ix)
+
+import qualified LLVM.Extra.Multi.Value.Memory as MultiMem
+import qualified LLVM.Extra.Multi.Value as MultiValue
+import qualified LLVM.Extra.Multi.Iterator as IterMV
+import qualified LLVM.Extra.Iterator as Iter
+import qualified LLVM.Extra.Arithmetic as A
+import LLVM.Extra.Multi.Value (atom)
+
+import qualified LLVM.Util.Loop as Loop
+import qualified LLVM.Core as LLVM
+
+import qualified Type.Data.Num.Decimal as TypeNum
+
+import Foreign.Storable (Storable)
+import Foreign.Ptr (Ptr)
+
+import Data.Word (Word8, Word16, Word32, Word64)
+import Data.Int (Int8, Int16, Int32, Int64)
+
+import qualified Control.Monad.HT as Monad
+import Control.Applicative ((<$>))
+
+
+type Size = Word64
+
+value :: (C sh, Expr.Value val) => sh -> val sh
+value = Expr.lift0 . MultiValue.cons
+
+paramWith ::
+   (Storable b, MultiMem.C b, Expr.Value val) =>
+   Param.T p b ->
+   (forall parameters.
+    (Storable parameters, MultiMem.C parameters) =>
+    (p -> parameters) ->
+    (MultiValue.T parameters -> val b) ->
+    a) ->
+   a
+paramWith p f =
+   Param.withMulti p (\get val -> f get (Expr.lift0 . val))
+
+load ::
+   (MultiMem.C sh) =>
+   f sh -> LLVM.Value (Ptr (MultiMem.Struct sh)) ->
+   LLVM.CodeGenFunction r (MultiValue.T sh)
+load _ = MultiMem.load
+
+intersect :: (C sh) => Exp sh -> Exp sh -> Exp sh
+intersect = Expr.liftM2 intersectCode
+
+offset ::
+   (C sh) =>
+   MultiValue.T sh -> MultiValue.T (Index sh) ->
+   LLVM.CodeGenFunction r (LLVM.Value Size)
+offset sh ix = ($ix) . snd =<< sizeOffset sh
+
+class (MultiValue.C sh, MultiValue.C (Index sh), Shape.Indexed sh) => C sh where
+   {-
+   It would be better to restrict zipWith to matching shapes
+   and turn shape intersection into a bound check.
+   -}
+   intersectCode ::
+      MultiValue.T sh -> MultiValue.T sh ->
+      LLVM.CodeGenFunction r (MultiValue.T sh)
+   size :: MultiValue.T sh -> LLVM.CodeGenFunction r (LLVM.Value Size)
+   {- |
+   Result is @(size, offset)@.
+   @size@ must equal the result of 'size'.
+   We use this for sharing intermediate results.
+   -}
+   sizeOffset ::
+      MultiValue.T sh ->
+      LLVM.CodeGenFunction r
+         (LLVM.Value Size,
+          MultiValue.T (Index sh) -> LLVM.CodeGenFunction r (LLVM.Value Size))
+   iterator :: (Index sh ~ ix) => MultiValue.T sh -> Iter.T r (MultiValue.T ix)
+   loop ::
+      (Index sh ~ ix, MultiValue.C ix, Loop.Phi state) =>
+      (MultiValue.T ix -> state -> LLVM.CodeGenFunction r state) ->
+      MultiValue.T sh -> state -> LLVM.CodeGenFunction r state
+   loop f sh = Iter.mapState_ f (iterator sh)
+
+
+instance C () where
+   intersectCode _ _ = return $ MultiValue.cons ()
+   size _ = return A.one
+   sizeOffset _ = return (A.one, \_ -> return A.zero)
+   iterator = Iter.singleton
+   loop = id
+
+
+class C sh => Scalar sh where
+   scalar :: (Expr.Value val) => val sh
+   zeroIndex :: (Expr.Value val) => f sh -> val (Index sh)
+
+instance Scalar () where
+   scalar = Expr.lift0 $ MultiValue.Cons ()
+   zeroIndex _ = Expr.lift0 $ MultiValue.Cons ()
+
+
+class
+   (C sh,
+    MultiValue.IntegerConstant (Index sh),
+    MultiValue.Additive (Index sh)) =>
+      Sequence sh where
+   sequenceShapeFromIndex ::
+      MultiValue.T (Index sh) -> LLVM.CodeGenFunction r (MultiValue.T sh)
+
+
+class
+   (MultiValue.Additive n, MultiValue.Real n, MultiValue.IntegerConstant n) =>
+      ToSize n where
+   toSize :: MultiValue.T n -> LLVM.CodeGenFunction r (LLVM.Value Size)
+
+instance ToSize Word8  where toSize (MultiValue.Cons n) = LLVM.ext n
+instance ToSize Word16 where toSize (MultiValue.Cons n) = LLVM.ext n
+instance ToSize Word32 where toSize (MultiValue.Cons n) = LLVM.ext n
+instance ToSize Word64 where toSize (MultiValue.Cons n) = return n
+instance ToSize Int8  where toSize (MultiValue.Cons n) = LLVM.zext n
+instance ToSize Int16 where toSize (MultiValue.Cons n) = LLVM.zext n
+instance ToSize Int32 where toSize (MultiValue.Cons n) = LLVM.zext n
+instance ToSize Int64 where toSize (MultiValue.Cons n) = LLVM.bitcast n
+
+
+unzipZeroBased :: MultiValue.T (ZeroBased n) -> ZeroBased (MultiValue.T n)
+unzipZeroBased (MultiValue.Cons (ZeroBased n)) = ZeroBased (MultiValue.Cons n)
+
+zeroBasedSize :: (Expr.Value val) => val (ZeroBased n) -> val n
+zeroBasedSize = Expr.lift1 $ Shape.zeroBasedSize . unzipZeroBased
+
+zeroBased :: (Expr.Value val) => val n -> val (ZeroBased n)
+zeroBased =
+   Expr.lift1 $ \(MultiValue.Cons n) -> MultiValue.Cons (ZeroBased n)
+
+instance (MultiValue.C n) => MultiValue.C (ZeroBased n) where
+   type Repr f (ZeroBased n) = ZeroBased (MultiValue.Repr f n)
+   cons (ZeroBased n) = zeroBased (MultiValue.cons n)
+   undef = zeroBased MultiValue.undef
+   zero = zeroBased MultiValue.zero
+   phis bb = Monad.lift zeroBased . MultiValue.phis bb . zeroBasedSize
+   addPhis bb a b = MultiValue.addPhis bb (zeroBasedSize a) (zeroBasedSize b)
+
+type instance
+   MultiValue.Decomposed f (ZeroBased pn) =
+      ZeroBased (MultiValue.Decomposed f pn)
+type instance
+   MultiValue.PatternTuple (ZeroBased pn) =
+      ZeroBased (MultiValue.PatternTuple pn)
+
+instance (MultiValue.Compose n) => MultiValue.Compose (ZeroBased n) where
+   type Composed (ZeroBased n) = ZeroBased (MultiValue.Composed n)
+   compose (ZeroBased n) = zeroBased (MultiValue.compose n)
+
+instance (MultiValue.Decompose pn) => MultiValue.Decompose (ZeroBased pn) where
+   decompose (ZeroBased p) sh =
+      MultiValue.decompose p <$> unzipZeroBased sh
+
+instance (Expr.Compose n) => Expr.Compose (ZeroBased n) where
+   type Composed (ZeroBased n) = ZeroBased (Expr.Composed n)
+   compose (ZeroBased n) = Expr.lift1 zeroBased (Expr.compose n)
+
+instance (Expr.Decompose pn) => Expr.Decompose (ZeroBased pn) where
+   decompose (ZeroBased p) = ZeroBased . Expr.decompose p . zeroBasedSize
+
+instance (MultiMem.C n) => MultiMem.C (ZeroBased n) where
+   type Struct (ZeroBased n) = MultiMem.Struct n
+   decompose = fmap zeroBased . MultiMem.decompose
+   compose = MultiMem.compose . zeroBasedSize
+
+{- |
+Array dimensions and indexes cannot be negative,
+but computations in indices may temporarily yield negative values
+or we want to add negative values to indices.
+
+So maybe, we would better have type Index (ZeroBased Word64) = Int64.
+This is not possible.
+Maybe we need an additional ZeroBased type for unsigned array sizes.
+-}
+instance
+      (Integral n, ToSize n, MultiValue.Comparison n) => C (ZeroBased n) where
+   intersectCode sha shb =
+      zeroBased <$> MultiValue.min (zeroBasedSize sha) (zeroBasedSize shb)
+   size = toSize . zeroBasedSize
+   sizeOffset sh = Monad.lift2 (,) (toSize $ zeroBasedSize sh) (return toSize)
+   iterator sh =
+      IterMV.take (zeroBasedSize sh) $
+      Iter.iterate MultiValue.inc MultiValue.zero
+
+instance
+   (Integral n, ToSize n, MultiValue.Comparison n) =>
+      Sequence (ZeroBased n) where
+   sequenceShapeFromIndex = return . zeroBased
+
+
+singletonRange :: n -> Range n
+singletonRange n = Range n n
+
+rangeSize ::
+   (ToSize n) =>
+   Range (MultiValue.T n) -> LLVM.CodeGenFunction r (LLVM.Value Size)
+rangeSize (Range from to) =
+   toSize =<< MultiValue.inc =<< MultiValue.sub to from
+
+unzipRange :: MultiValue.T (Range n) -> Range (MultiValue.T n)
+unzipRange (MultiValue.Cons (Range from to)) =
+   Range (MultiValue.Cons from) (MultiValue.Cons to)
+
+zipRange :: MultiValue.T n -> MultiValue.T n -> MultiValue.T (Range n)
+zipRange (MultiValue.Cons from) (MultiValue.Cons to) =
+   MultiValue.Cons (Range from to)
+
+
+rangeFrom :: (Expr.Value val) => val (Range n) -> val n
+rangeFrom = Expr.lift1 $ Shape.rangeFrom . unzipRange
+
+rangeTo :: (Expr.Value val) => val (Range n) -> val n
+rangeTo = Expr.lift1 $ Shape.rangeTo . unzipRange
+
+range :: (Expr.Value val) => val n -> val n -> val (Range n)
+range =
+   Expr.lift2 $
+      \(MultiValue.Cons from) (MultiValue.Cons to) ->
+         MultiValue.Cons (Range from to)
+
+
+instance (MultiValue.C n) => MultiValue.C (Range n) where
+   type Repr f (Range n) = Range (MultiValue.Repr f n)
+   cons (Range from to) = zipRange (MultiValue.cons from) (MultiValue.cons to)
+   undef = MultiValue.compose $ singletonRange MultiValue.undef
+   zero = MultiValue.compose $ singletonRange MultiValue.zero
+   phis bb a =
+      case unzipRange a of
+         Range a0 a1 ->
+            Monad.lift2 zipRange (MultiValue.phis bb a0) (MultiValue.phis bb a1)
+   addPhis bb a b =
+      case (unzipRange a, unzipRange b) of
+         (Range a0 a1, Range b0 b1) ->
+            MultiValue.addPhis bb a0 b0 >>
+            MultiValue.addPhis bb a1 b1
+
+type instance
+   MultiValue.Decomposed f (Range pn) =
+      Range (MultiValue.Decomposed f pn)
+type instance
+   MultiValue.PatternTuple (Range pn) =
+      Range (MultiValue.PatternTuple pn)
+
+instance (MultiValue.Compose n) => MultiValue.Compose (Range n) where
+   type Composed (Range n) = Range (MultiValue.Composed n)
+   compose (Range from to) =
+      zipRange (MultiValue.compose from) (MultiValue.compose to)
+
+instance (MultiValue.Decompose pn) => MultiValue.Decompose (Range pn) where
+   decompose (Range pfrom pto) rng =
+      case unzipRange rng of
+         Range from to ->
+            Range
+               (MultiValue.decompose pfrom from)
+               (MultiValue.decompose pto to)
+
+instance (MultiMem.C n) => MultiMem.C (Range n) where
+   type Struct (Range n) = PairStruct n
+   decompose = fmap (uncurry zipRange) . decomposeGen
+   compose x = case unzipRange x of Range n m -> composeGen n m
+
+instance (Ix n, ToSize n, MultiValue.Comparison n) => C (Range n) where
+   intersectCode =
+      MultiValue.modifyF2 (singletonRange atom) (singletonRange atom) $
+            \(Range fromN toN) (Range fromM toM) ->
+         Monad.lift2 Range (MultiValue.max fromN fromM) (MultiValue.min toN toM)
+   size = rangeSize . unzipRange
+   sizeOffset rngValue =
+      case unzipRange rngValue of
+         rng@(Range from _to) ->
+            Monad.lift2 (,) (rangeSize rng)
+               (return $ \i -> toSize =<< MultiValue.sub i from)
+   iterator rngValue =
+      case MultiValue.decompose (singletonRange atom) rngValue of
+         Range from to ->
+            IterMV.takeWhile (MultiValue.cmp LLVM.CmpGE to) $
+            Iter.iterate MultiValue.inc from
+
+
+singletonShifted :: n -> Shifted n
+singletonShifted n = Shifted n n
+
+
+unzipShifted :: MultiValue.T (Shifted n) -> Shifted (MultiValue.T n)
+unzipShifted (MultiValue.Cons (Shifted from to)) =
+   Shifted (MultiValue.Cons from) (MultiValue.Cons to)
+
+zipShifted :: MultiValue.T n -> MultiValue.T n -> MultiValue.T (Shifted n)
+zipShifted (MultiValue.Cons from) (MultiValue.Cons to) =
+   MultiValue.Cons (Shifted from to)
+
+
+shiftedOffset :: (Expr.Value val) => val (Shifted n) -> val n
+shiftedOffset = Expr.lift1 $ Shape.shiftedOffset . unzipShifted
+
+shiftedSize :: (Expr.Value val) => val (Shifted n) -> val n
+shiftedSize = Expr.lift1 $ Shape.shiftedSize . unzipShifted
+
+shifted :: (Expr.Value val) => val n -> val n -> val (Shifted n)
+shifted =
+   Expr.lift2 $
+      \(MultiValue.Cons from) (MultiValue.Cons to) ->
+         MultiValue.Cons (Shifted from to)
+
+
+instance (MultiValue.C n) => MultiValue.C (Shifted n) where
+   type Repr f (Shifted n) = Shifted (MultiValue.Repr f n)
+   cons (Shifted start len) =
+      zipShifted (MultiValue.cons start) (MultiValue.cons len)
+   undef = MultiValue.compose $ singletonShifted MultiValue.undef
+   zero = MultiValue.compose $ singletonShifted MultiValue.zero
+   phis bb a =
+      case unzipShifted a of
+         Shifted a0 a1 ->
+            Monad.lift2 zipShifted
+               (MultiValue.phis bb a0) (MultiValue.phis bb a1)
+   addPhis bb a b =
+      case (unzipShifted a, unzipShifted b) of
+         (Shifted a0 a1, Shifted b0 b1) ->
+            MultiValue.addPhis bb a0 b0 >>
+            MultiValue.addPhis bb a1 b1
+
+type instance
+   MultiValue.Decomposed f (Shifted pn) =
+      Shifted (MultiValue.Decomposed f pn)
+type instance
+   MultiValue.PatternTuple (Shifted pn) =
+      Shifted (MultiValue.PatternTuple pn)
+
+instance (MultiValue.Compose n) => MultiValue.Compose (Shifted n) where
+   type Composed (Shifted n) = Shifted (MultiValue.Composed n)
+   compose (Shifted start len) =
+      zipShifted (MultiValue.compose start) (MultiValue.compose len)
+
+instance (MultiValue.Decompose pn) => MultiValue.Decompose (Shifted pn) where
+   decompose (Shifted pstart plen) rng =
+      case unzipShifted rng of
+         Shifted start len ->
+            Shifted
+               (MultiValue.decompose pstart start)
+               (MultiValue.decompose plen len)
+
+instance (MultiMem.C n) => MultiMem.C (Shifted n) where
+   type Struct (Shifted n) = PairStruct n
+   decompose = fmap (uncurry zipShifted) . decomposeGen
+   compose x = case unzipShifted x of Shifted n m -> composeGen n m
+
+type PairStruct n = LLVM.Struct (MultiMem.Struct n, (MultiMem.Struct n, ()))
+
+decomposeGen ::
+   (MultiMem.C n) =>
+   LLVM.Value (PairStruct n) ->
+   LLVM.CodeGenFunction r (MultiValue.T n, MultiValue.T n)
+decomposeGen nm =
+   Monad.lift2 (,)
+      (MultiMem.decompose =<< LLVM.extractvalue nm TypeNum.d0)
+      (MultiMem.decompose =<< LLVM.extractvalue nm TypeNum.d1)
+
+composeGen ::
+   (MultiMem.C n) =>
+   MultiValue.T n -> MultiValue.T n ->
+   LLVM.CodeGenFunction r (LLVM.Value (PairStruct n))
+composeGen n m = do
+   sn <- MultiMem.compose n
+   sm <- MultiMem.compose m
+   rn <- LLVM.insertvalue (LLVM.value LLVM.undef) sn TypeNum.d0
+   LLVM.insertvalue rn sm TypeNum.d1
+
+
+instance (Integral n, ToSize n, MultiValue.Comparison n) => C (Shifted n) where
+   intersectCode =
+      MultiValue.modifyF2 (singletonShifted atom) (singletonShifted atom) $
+            \(Shifted startN lenN) (Shifted startM lenM) -> do
+         start <- MultiValue.max startN startM
+         endN <- MultiValue.add startN lenN
+         endM <- MultiValue.add startM lenM
+         end <- MultiValue.min endN endM
+         Shifted start <$> MultiValue.sub end start
+   size = toSize . shiftedSize
+   sizeOffset shapeValue =
+      case unzipShifted shapeValue of
+         Shifted start len ->
+            Monad.lift2 (,) (toSize len)
+               (return $ \i -> toSize =<< MultiValue.sub i start)
+   iterator rngValue =
+      case MultiValue.decompose (singletonShifted atom) rngValue of
+         Shifted from len ->
+            IterMV.take len $ Iter.iterate MultiValue.inc from
+
+
+
+instance (C n, C m) => C (n,m) where
+   intersectCode a b =
+      case (MultiValue.unzip a, MultiValue.unzip b) of
+         ((an,am), (bn,bm)) ->
+            Monad.lift2 MultiValue.zip
+               (intersectCode an bn)
+               (intersectCode am bm)
+   size nm =
+      case MultiValue.unzip nm of
+         (n,m) -> Monad.liftJoin2 A.mul (size n) (size m)
+   sizeOffset nm =
+      case MultiValue.unzip nm of
+         (n,m) -> do
+            (ns, iOffset) <- sizeOffset n
+            (ms, jOffset) <- sizeOffset m
+            sz <- A.mul ns ms
+            return
+               (sz,
+                \ij ->
+                  case MultiValue.unzip ij of
+                     (i,j) -> do
+                        il <- iOffset i
+                        jl <- jOffset j
+                        A.add jl =<< A.mul ms il)
+   iterator nm =
+      case MultiValue.unzip nm of
+         (n,m) ->
+            uncurry MultiValue.zip <$>
+            Iter.cartesian (iterator n) (iterator m)
+   loop code nm =
+      case MultiValue.unzip nm of
+         (n,m) -> loop (\i -> loop (\j -> code (MultiValue.zip i j)) m) n
+
+instance (C n, C m, C l) => C (n,m,l) where
+   intersectCode a b =
+      case (MultiValue.unzip3 a, MultiValue.unzip3 b) of
+         ((ai,aj,ak), (bi,bj,bk)) ->
+            Monad.lift3 MultiValue.zip3
+               (intersectCode ai bi)
+               (intersectCode aj bj)
+               (intersectCode ak bk)
+   size nml =
+      case MultiValue.unzip3 nml of
+         (n,m,l) ->
+            Monad.liftJoin2 A.mul (size n) $
+            Monad.liftJoin2 A.mul (size m) (size l)
+   sizeOffset nml =
+      case MultiValue.unzip3 nml of
+         (n,m,l) -> do
+            (ns, iOffset) <- sizeOffset n
+            (ms, jOffset) <- sizeOffset m
+            (ls, kOffset) <- sizeOffset l
+            sz <- A.mul ns =<< A.mul ms ls
+            return
+               (sz,
+                \ijk ->
+                  case MultiValue.unzip3 ijk of
+                     (i,j,k) -> do
+                        il <- iOffset i
+                        jl <- jOffset j
+                        kl <- kOffset k
+                        A.add kl =<< A.mul ls =<< A.add jl =<< A.mul ms il)
+   iterator nml =
+      case MultiValue.unzip3 nml of
+         (n,m,l) ->
+            fmap (\(a,(b,c)) -> MultiValue.zip3 a b c) $
+            Iter.cartesian (iterator n) $
+            Iter.cartesian (iterator m) (iterator l)
+   loop code nml =
+      case MultiValue.unzip3 nml of
+         (n,m,l) ->
+            loop (\i -> loop (\j -> loop (\k ->
+               code (MultiValue.zip3 i j k))
+            l) m) n
diff --git a/src/Data/Array/Knead/Shape/Cubic.hs b/src/Data/Array/Knead/Shape/Cubic.hs
--- a/src/Data/Array/Knead/Shape/Cubic.hs
+++ b/src/Data/Array/Knead/Shape/Cubic.hs
@@ -15,7 +15,7 @@
    constant,
    paramWith,
    tunnel,
-   flattenIndex,
+   offsetCode,
    peek,
    poke,
    computeSize,
@@ -34,13 +34,16 @@
    storeMultiValue,
    ) where
 
-import qualified Data.Array.Knead.Shape.Nested as Shape
+import qualified Data.Array.Knead.Shape as Shape
 import qualified Data.Array.Knead.Shape.Cubic.Int as Index
 
 import qualified Data.Array.Knead.Parameter as Param
 import qualified Data.Array.Knead.Expression as Expr
 import Data.Array.Knead.Expression (Exp, )
 
+import qualified Data.Array.Comfort.Shape as ComfortShape
+import Data.Array.Comfort.Shape (ZeroBased(ZeroBased))
+
 import qualified LLVM.Extra.Multi.Value.Memory as MultiValueMemory
 import qualified LLVM.Extra.Multi.Value as MultiValue
 import qualified LLVM.Extra.Multi.Iterator as IterMV
@@ -87,7 +90,8 @@
       (runSwitchInt $ Index.switchSingle (SwitchInt cons0))
 
 
-newtype Op2 tag sh = Op2 {runOp2 :: Exp (T tag sh) -> Exp (T tag sh) -> Exp (T tag sh)}
+newtype Op2 tag sh =
+         Op2 {runOp2 :: Exp (T tag sh) -> Exp (T tag sh) -> Exp (T tag sh)}
 
 intersect :: C sh => Exp (Shape sh) -> Exp (Shape sh) -> Exp (Shape sh)
 intersect =
@@ -104,7 +108,8 @@
 _value = Expr.lift0 . MultiValue.cons
 
 
-newtype MakeValue val tag sh = MakeValue {runMakeValue :: T tag sh -> val (T tag sh)}
+newtype MakeValue val tag sh =
+         MakeValue {runMakeValue :: T tag sh -> val (T tag sh)}
 
 value :: (C sh, Expr.Value val) => T tag sh -> val (T tag sh)
 value =
@@ -285,16 +290,57 @@
    addPhis = addPhis
    phis = phis
 
-instance (tag ~ ShapeTag, C sh) => Shape.C (T tag sh) where
-   type Index (T tag sh) = Index sh
+instance (tag ~ ShapeTag, C sh) => ComfortShape.C (T tag sh) where
    size = fromIntegral . size . decons
-   sizeCode = computeSize
+
+instance (tag ~ ShapeTag, C sh) => ComfortShape.Indexed (T tag sh) where
+   type Index (T tag sh) = Index sh
+   indices (Cons ix) = map index $ indices ix
+   inBounds (Cons sh) (Cons ix) = inBounds sh ix
+   offset (Cons sh) (Cons ix) = offset sh ix
+
+
+newtype Indices sh = Indices {runIndices :: sh -> [sh]}
+
+indices :: (C sh) => sh -> [sh]
+indices =
+   runIndices $
+   switchInt
+      (Indices $ \Z -> [Z])
+      (Indices $ \(t :. Index.Int h) ->
+         liftM2 (:.) (indices t)
+            (map Index.Int $ ComfortShape.indices $ ZeroBased h))
+
+newtype InBounds sh = InBounds {runInBounds :: sh -> sh -> Bool}
+
+inBounds :: (C sh) => sh -> sh -> Bool
+inBounds =
+   runInBounds $
+   switchInt
+      (InBounds $ \Z Z -> True)
+      (InBounds $ \(sh :. Index.Int s) (ix :. Index.Int i) ->
+         inBounds sh ix && ComfortShape.inBounds (ZeroBased s) i)
+
+newtype Offset sh = Offset {runOffset :: sh -> sh -> Int}
+
+offset :: (C sh) => sh -> sh -> Int
+offset =
+   runOffset $
+   switchInt
+      (Offset $ \Z Z -> 0)
+      (Offset $ \(sh :. Index.Int s) (ix :. Index.Int i) ->
+         offset sh ix * fromIntegral s + fromIntegral i)
+
+
+
+instance (tag ~ ShapeTag, C sh) => Shape.C (T tag sh) where
+   size = computeSize
    intersectCode = Expr.unliftM2 intersect
-   flattenIndexRec sh ix =
-      -- a joint implementation would not be more efficient
+   sizeOffset sh =
+      -- would a joint implementation be more efficient?
       liftM2 (,)
          (computeSize sh)
-         (flattenIndex sh ix)
+         (return $ offsetCode sh)
    iterator = iterator
    loop = loop
 
@@ -327,25 +373,25 @@
       StructFieldsProp -> store x =<< castPtrValue ptr
 
 
-newtype FlattenIndex r sh =
-   FlattenIndex {
-      runFlattenIndex ::
+newtype OffsetCode r sh =
+   OffsetCode {
+      runOffsetCode ::
          MultiValue.T (Shape sh) -> MultiValue.T (Index sh) ->
          LLVM.CodeGenFunction r (LLVM.Value Shape.Size)
    }
 
-flattenIndex ::
+offsetCode ::
    (C sh) =>
    MultiValue.T (Shape sh) -> MultiValue.T (Index sh) ->
    LLVM.CodeGenFunction r (LLVM.Value Shape.Size)
-flattenIndex =
-   runFlattenIndex $
+offsetCode =
+   runOffsetCode $
    switchInt
-      (FlattenIndex $ \_zerosh _zeroix -> return A.zero)
-      (FlattenIndex $
+      (OffsetCode $ \_zerosh _zeroix -> return A.zero)
+      (OffsetCode $
          switchR $ \sh (MultiValue.Cons s) ->
          switchR $ \ix (MultiValue.Cons i) ->
-            A.add i =<< A.mul s =<< flattenIndex sh ix)
+            A.add i =<< A.mul s =<< offsetCode sh ix)
 
 
 newtype Rank sh = Rank {runRank :: sh -> Int}
@@ -432,10 +478,7 @@
          store sh =<< A.advanceArrayElementPtr ptr)
 
 
-newtype Size sh =
-   Size {
-      runSize :: sh -> Shape.Size
-   }
+newtype Size sh = Size {runSize :: sh -> Shape.Size}
 
 size :: (C sh) => sh -> Shape.Size
 size =
diff --git a/src/Data/Array/Knead/Shape/Nested.hs b/src/Data/Array/Knead/Shape/Nested.hs
deleted file mode 100644
--- a/src/Data/Array/Knead/Shape/Nested.hs
+++ /dev/null
@@ -1,539 +0,0 @@
-{-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE Rank2Types #-}
-{-# LANGUAGE FlexibleContexts #-}
-module Data.Array.Knead.Shape.Nested (
-   C(..),
-   Size,
-   value,
-   paramWith,
-   load,
-   intersect,
-   flattenIndex,
-
-   Range(..),
-   Shifted(..),
-
-   Scalar(..),
-   Sequence(..),
-   ) where
-
-import qualified Data.Array.Knead.Expression as Expr
-import qualified Data.Array.Knead.Parameter as Param
-import Data.Array.Knead.Expression (Exp, )
-
-import qualified LLVM.Extra.Multi.Value.Memory as MultiMem
-import qualified LLVM.Extra.Multi.Value as MultiValue
-import qualified LLVM.Extra.Multi.Iterator as IterMV
-import qualified LLVM.Extra.Iterator as Iter
-import qualified LLVM.Extra.Arithmetic as A
-import LLVM.Extra.Multi.Value (atom)
-
-import qualified LLVM.Util.Loop as Loop
-import qualified LLVM.Core as LLVM
-
-import qualified Type.Data.Num.Decimal as TypeNum
-
-import Foreign.Storable
-         (Storable, sizeOf, alignment, poke, peek, pokeElemOff, peekElemOff)
-import Foreign.Ptr (Ptr, castPtr)
-
-import Data.Word (Word8, Word16, Word32, Word64)
-import Data.Int (Int8, Int16, Int32, Int64)
-
-import qualified Control.Monad.HT as Monad
-import Control.Applicative ((<$>))
-
-
-type Size = Word64
-
-value :: (C sh, Expr.Value val) => sh -> val sh
-value = Expr.lift0 . MultiValue.cons
-
-paramWith ::
-   (Storable b, MultiMem.C b, Expr.Value val) =>
-   Param.T p b ->
-   (forall parameters.
-    (Storable parameters, MultiMem.C parameters) =>
-    (p -> parameters) ->
-    (MultiValue.T parameters -> val b) ->
-    a) ->
-   a
-paramWith p f =
-   Param.withMulti p (\get val -> f get (Expr.lift0 . val))
-
-load ::
-   (MultiMem.C sh) =>
-   f sh -> LLVM.Value (Ptr (MultiMem.Struct sh)) ->
-   LLVM.CodeGenFunction r (MultiValue.T sh)
-load _ = MultiMem.load
-
-intersect :: (C sh) => Exp sh -> Exp sh -> Exp sh
-intersect = Expr.liftM2 intersectCode
-
-flattenIndex ::
-   (C sh) =>
-   MultiValue.T sh -> MultiValue.T (Index sh) ->
-   LLVM.CodeGenFunction r (LLVM.Value Size)
-flattenIndex sh ix =
-   fmap snd $ flattenIndexRec sh ix
-
-class (MultiValue.C sh, MultiValue.C (Index sh)) => C sh where
-   type Index sh :: *
-   {-
-   It would be better to restrict zipWith to matching shapes
-   and turn shape intersection into a bound check.
-   -}
-   intersectCode ::
-      MultiValue.T sh -> MultiValue.T sh ->
-      LLVM.CodeGenFunction r (MultiValue.T sh)
-   sizeCode ::
-      MultiValue.T sh ->
-      LLVM.CodeGenFunction r (LLVM.Value Size)
-   size :: sh -> Int
-   {- |
-   Result is @(size, flattenedIndex)@.
-   @size@ must equal the result of 'sizeCode'.
-   We use this for sharing intermediate results.
-   -}
-   flattenIndexRec ::
-      MultiValue.T sh -> MultiValue.T (Index sh) ->
-      LLVM.CodeGenFunction r (LLVM.Value Size, LLVM.Value Size)
-   iterator :: (Index sh ~ ix) => MultiValue.T sh -> Iter.T r (MultiValue.T ix)
-   loop ::
-      (Index sh ~ ix, MultiValue.C ix, Loop.Phi state) =>
-      (MultiValue.T ix -> state -> LLVM.CodeGenFunction r state) ->
-      MultiValue.T sh -> state -> LLVM.CodeGenFunction r state
-   loop f sh = Iter.mapState_ f (iterator sh)
-
-
-instance C () where
-   type Index () = ()
-   intersectCode _ _ = return $ MultiValue.cons ()
-   sizeCode _ = return A.one
-   size _ = 1
-   flattenIndexRec _ _ = return (A.one, A.zero)
-   iterator = Iter.singleton
-   loop = id
-
-
-class C sh => Scalar sh where
-   scalar :: (Expr.Value val) => val sh
-   zeroIndex :: (Expr.Value val) => f sh -> val (Index sh)
-
-instance Scalar () where
-   scalar = Expr.lift0 $ MultiValue.Cons ()
-   zeroIndex _ = Expr.lift0 $ MultiValue.Cons ()
-
-
-class
-   (C sh,
-    MultiValue.IntegerConstant (Index sh),
-    MultiValue.Additive (Index sh)) =>
-      Sequence sh where
-   sequenceShapeFromIndex ::
-      MultiValue.T (Index sh) -> LLVM.CodeGenFunction r (MultiValue.T sh)
-
-
-iteratorPrimitive ::
-   (MultiValue.Repr LLVM.Value j ~ LLVM.Value j,
-    Num j, LLVM.IsConst j, LLVM.IsInteger j,
-    LLVM.CmpRet j, LLVM.CmpResult j ~ Bool,
-    MultiValue.Additive j, MultiValue.IntegerConstant j) =>
-   MultiValue.T j -> Iter.T r (MultiValue.T j)
-iteratorPrimitive (MultiValue.Cons n) = iteratorStart n MultiValue.zero
-
-iteratorStart ::
-   (Num j, LLVM.IsConst j, LLVM.IsInteger j,
-    LLVM.CmpRet j, LLVM.CmpResult j ~ Bool,
-    MultiValue.Additive i, MultiValue.IntegerConstant i) =>
-   LLVM.Value j -> MultiValue.T i -> Iter.T r (MultiValue.T i)
-iteratorStart n start = Iter.take n $ Iter.iterate MultiValue.inc start
-
-
-instance C Word32 where
-   type Index Word32 = Word32
-   intersectCode = MultiValue.min
-   sizeCode (MultiValue.Cons n) = LLVM.ext n
-   size = fromIntegral
-   flattenIndexRec (MultiValue.Cons n) (MultiValue.Cons i) =
-      Monad.lift2 (,) (LLVM.ext n) (LLVM.ext i)
-   iterator = iteratorPrimitive
-
-instance Sequence Word32 where
-   sequenceShapeFromIndex = return
-
-instance C Word64 where
-   type Index Word64 = Word64
-   intersectCode = MultiValue.min
-   sizeCode (MultiValue.Cons n) = return n
-   size = fromIntegral
-   flattenIndexRec (MultiValue.Cons n) (MultiValue.Cons i) = return (n, i)
-   iterator = iteratorPrimitive
-
-instance Sequence Word64 where
-   sequenceShapeFromIndex = return
-
-
-{- |
-Array dimensions and indexes cannot be negative,
-but computations in indices may temporarily yield negative values
-or we want to add negative values to indices.
-
-Maybe we should better have type Index Word64 = Int64?
--}
-unsigned8 :: LLVM.Value Int8 -> LLVM.CodeGenFunction r (LLVM.Value Word8)
-unsigned8 = LLVM.bitcast
-
-instance C Int8 where
-   type Index Int8 = Int8
-   intersectCode = MultiValue.min
-   sizeCode (MultiValue.Cons n) = LLVM.ext =<< unsigned8 n
-   size = fromIntegral
-   flattenIndexRec (MultiValue.Cons n) (MultiValue.Cons i) =
-      Monad.lift2 (,) (LLVM.ext =<< unsigned8 n) (LLVM.ext =<< unsigned8 i)
-   iterator = iteratorPrimitive
-
-instance Sequence Int8 where
-   sequenceShapeFromIndex = return
-
-unsigned16 :: LLVM.Value Int16 -> LLVM.CodeGenFunction r (LLVM.Value Word16)
-unsigned16 = LLVM.bitcast
-
-instance C Int16 where
-   type Index Int16 = Int16
-   intersectCode = MultiValue.min
-   sizeCode (MultiValue.Cons n) = LLVM.ext =<< unsigned16 n
-   size = fromIntegral
-   flattenIndexRec (MultiValue.Cons n) (MultiValue.Cons i) =
-      Monad.lift2 (,) (LLVM.ext =<< unsigned16 n) (LLVM.ext =<< unsigned16 i)
-   iterator = iteratorPrimitive
-
-instance Sequence Int16 where
-   sequenceShapeFromIndex = return
-
-instance C Int32 where
-   type Index Int32 = Int32
-   intersectCode = MultiValue.min
-   sizeCode (MultiValue.Cons n) = LLVM.zext n
-   size = fromIntegral
-   flattenIndexRec (MultiValue.Cons n) (MultiValue.Cons i) =
-      Monad.lift2 (,) (LLVM.zext n) (LLVM.zext i)
-   iterator = iteratorPrimitive
-
-instance Sequence Int32 where
-   sequenceShapeFromIndex = return
-
-instance C Int64 where
-   type Index Int64 = Int64
-   intersectCode = MultiValue.min
-   sizeCode (MultiValue.Cons n) = LLVM.bitcast n
-   size = fromIntegral
-   flattenIndexRec (MultiValue.Cons n) (MultiValue.Cons i) =
-      Monad.lift2 (,) (LLVM.bitcast n) (LLVM.bitcast i)
-   iterator = iteratorPrimitive
-
-instance Sequence Int64 where
-   sequenceShapeFromIndex = return
-
-
-{- |
-'Range' denotes an inclusive range like
-those of the Haskell 98 standard @Array@ type from the @array@ package.
-E.g. the shape type @(Range Int32, Range Int64)@
-is equivalent to the ix type @(Int32, Int64)@ for @Array@s.
--}
-data Range n = Range n n
-
-singletonRange :: n -> Range n
-singletonRange n = Range n n
-
-
-{-# INLINE castToElemPtr #-}
-castToElemPtr :: Ptr (f a) -> Ptr a
-castToElemPtr = castPtr
-
--- cf. sample-frame:Stereo
-instance Storable n => Storable (Range n) where
-   {-# INLINE sizeOf #-}
-   {-# INLINE alignment #-}
-   {-# INLINE peek #-}
-   {-# INLINE poke #-}
-   sizeOf ~(Range l r) = sizeOf l + mod (- sizeOf l) (alignment r) + sizeOf r
-   alignment ~(Range l _) = alignment l
-   poke p (Range l r) =
-      let q = castToElemPtr p
-      in  poke q l >> pokeElemOff q 1 r
-   peek p =
-      let q = castToElemPtr p
-      in  Monad.lift2 Range (peek q) (peekElemOff q 1)
-
-
-class
-   (MultiValue.Additive n, MultiValue.Real n, MultiValue.IntegerConstant n) =>
-      ToSize n where
-   toSize :: MultiValue.T n -> LLVM.CodeGenFunction r (LLVM.Value Size)
-
-instance ToSize Word32 where toSize (MultiValue.Cons n) = LLVM.ext n
-instance ToSize Word64 where toSize (MultiValue.Cons n) = return n
-instance ToSize Int32 where toSize (MultiValue.Cons n) = LLVM.zext n
-instance ToSize Int64 where toSize (MultiValue.Cons n) = LLVM.bitcast n
-
-rangeSize ::
-   (ToSize n) =>
-   Range (MultiValue.T n) -> LLVM.CodeGenFunction r (LLVM.Value Size)
-rangeSize (Range from to) =
-   toSize =<< MultiValue.inc =<< MultiValue.sub to from
-
-unzipRange :: MultiValue.T (Range n) -> Range (MultiValue.T n)
-unzipRange (MultiValue.Cons (Range from to)) =
-   Range (MultiValue.Cons from) (MultiValue.Cons to)
-
-zipRange :: MultiValue.T n -> MultiValue.T n -> MultiValue.T (Range n)
-zipRange (MultiValue.Cons from) (MultiValue.Cons to) =
-   MultiValue.Cons (Range from to)
-
-instance (MultiValue.C n) => MultiValue.C (Range n) where
-   type Repr f (Range n) = Range (MultiValue.Repr f n)
-   cons (Range from to) = zipRange (MultiValue.cons from) (MultiValue.cons to)
-   undef = MultiValue.compose $ singletonRange MultiValue.undef
-   zero = MultiValue.compose $ singletonRange MultiValue.zero
-   phis bb a =
-      case unzipRange a of
-         Range a0 a1 ->
-            Monad.lift2 zipRange (MultiValue.phis bb a0) (MultiValue.phis bb a1)
-   addPhis bb a b =
-      case (unzipRange a, unzipRange b) of
-         (Range a0 a1, Range b0 b1) ->
-            MultiValue.addPhis bb a0 b0 >>
-            MultiValue.addPhis bb a1 b1
-
-type instance
-   MultiValue.Decomposed f (Range pn) =
-      Range (MultiValue.Decomposed f pn)
-type instance
-   MultiValue.PatternTuple (Range pn) =
-      Range (MultiValue.PatternTuple pn)
-
-instance (MultiValue.Compose n) => MultiValue.Compose (Range n) where
-   type Composed (Range n) = Range (MultiValue.Composed n)
-   compose (Range from to) =
-      zipRange (MultiValue.compose from) (MultiValue.compose to)
-
-instance (MultiValue.Decompose pn) => MultiValue.Decompose (Range pn) where
-   decompose (Range pfrom pto) rng =
-      case unzipRange rng of
-         Range from to ->
-            Range
-               (MultiValue.decompose pfrom from)
-               (MultiValue.decompose pto to)
-
-instance (MultiMem.C n) => MultiMem.C (Range n) where
-   type Struct (Range n) = PairStruct n
-   decompose = fmap (uncurry zipRange) . decomposeGen
-   compose x = case unzipRange x of Range n m -> composeGen n m
-
-instance (Integral n, ToSize n, MultiValue.Comparison n) => C (Range n) where
-   type Index (Range n) = n
-   intersectCode =
-      MultiValue.modifyF2 (singletonRange atom) (singletonRange atom) $
-            \(Range fromN toN) (Range fromM toM) ->
-         Monad.lift2 Range (MultiValue.max fromN fromM) (MultiValue.min toN toM)
-   sizeCode = rangeSize . unzipRange
-   size (Range from to) = fromIntegral $ to-from+1
-   flattenIndexRec rngValue i =
-      case unzipRange rngValue of
-         rng@(Range from _to) ->
-            Monad.lift2 (,) (rangeSize rng) (toSize =<< MultiValue.sub i from)
-   iterator rngValue =
-      case MultiValue.decompose (singletonRange atom) rngValue of
-         Range from to ->
-            IterMV.takeWhile (MultiValue.cmp LLVM.CmpGE to) $
-            Iter.iterate MultiValue.inc from
-
-
-{- |
-'Shifted' denotes a range defined by the start index and the length.
--}
-data Shifted n = Shifted {shiftedOffset, shiftedSize :: n}
-
-singletonShifted :: n -> Shifted n
-singletonShifted n = Shifted n n
-
-
--- cf. sample-frame:Stereo
-instance Storable n => Storable (Shifted n) where
-   {-# INLINE sizeOf #-}
-   {-# INLINE alignment #-}
-   {-# INLINE peek #-}
-   {-# INLINE poke #-}
-   sizeOf ~(Shifted l n) = sizeOf l + mod (- sizeOf l) (alignment n) + sizeOf n
-   alignment ~(Shifted l _) = alignment l
-   poke p (Shifted l n) =
-      let q = castToElemPtr p
-      in  poke q l >> pokeElemOff q 1 n
-   peek p =
-      let q = castToElemPtr p
-      in  Monad.lift2 Shifted (peek q) (peekElemOff q 1)
-
-
-unzipShifted :: MultiValue.T (Shifted n) -> Shifted (MultiValue.T n)
-unzipShifted (MultiValue.Cons (Shifted from to)) =
-   Shifted (MultiValue.Cons from) (MultiValue.Cons to)
-
-zipShifted :: MultiValue.T n -> MultiValue.T n -> MultiValue.T (Shifted n)
-zipShifted (MultiValue.Cons from) (MultiValue.Cons to) =
-   MultiValue.Cons (Shifted from to)
-
-instance (MultiValue.C n) => MultiValue.C (Shifted n) where
-   type Repr f (Shifted n) = Shifted (MultiValue.Repr f n)
-   cons (Shifted offset len) =
-      zipShifted (MultiValue.cons offset) (MultiValue.cons len)
-   undef = MultiValue.compose $ singletonShifted MultiValue.undef
-   zero = MultiValue.compose $ singletonShifted MultiValue.zero
-   phis bb a =
-      case unzipShifted a of
-         Shifted a0 a1 ->
-            Monad.lift2 zipShifted
-               (MultiValue.phis bb a0) (MultiValue.phis bb a1)
-   addPhis bb a b =
-      case (unzipShifted a, unzipShifted b) of
-         (Shifted a0 a1, Shifted b0 b1) ->
-            MultiValue.addPhis bb a0 b0 >>
-            MultiValue.addPhis bb a1 b1
-
-type instance
-   MultiValue.Decomposed f (Shifted pn) =
-      Shifted (MultiValue.Decomposed f pn)
-type instance
-   MultiValue.PatternTuple (Shifted pn) =
-      Shifted (MultiValue.PatternTuple pn)
-
-instance (MultiValue.Compose n) => MultiValue.Compose (Shifted n) where
-   type Composed (Shifted n) = Shifted (MultiValue.Composed n)
-   compose (Shifted offset len) =
-      zipShifted (MultiValue.compose offset) (MultiValue.compose len)
-
-instance (MultiValue.Decompose pn) => MultiValue.Decompose (Shifted pn) where
-   decompose (Shifted poffset plen) rng =
-      case unzipShifted rng of
-         Shifted offset len ->
-            Shifted
-               (MultiValue.decompose poffset offset)
-               (MultiValue.decompose plen len)
-
-instance (MultiMem.C n) => MultiMem.C (Shifted n) where
-   type Struct (Shifted n) = PairStruct n
-   decompose = fmap (uncurry zipShifted) . decomposeGen
-   compose x = case unzipShifted x of Shifted n m -> composeGen n m
-
-type PairStruct n = LLVM.Struct (MultiMem.Struct n, (MultiMem.Struct n, ()))
-
-decomposeGen ::
-   (MultiMem.C n) =>
-   LLVM.Value (PairStruct n) ->
-   LLVM.CodeGenFunction r (MultiValue.T n, MultiValue.T n)
-decomposeGen nm =
-   Monad.lift2 (,)
-      (MultiMem.decompose =<< LLVM.extractvalue nm TypeNum.d0)
-      (MultiMem.decompose =<< LLVM.extractvalue nm TypeNum.d1)
-
-composeGen ::
-   (MultiMem.C n) =>
-   MultiValue.T n -> MultiValue.T n ->
-   LLVM.CodeGenFunction r (LLVM.Value (PairStruct n))
-composeGen n m = do
-   sn <- MultiMem.compose n
-   sm <- MultiMem.compose m
-   rn <- LLVM.insertvalue (LLVM.value LLVM.undef) sn TypeNum.d0
-   LLVM.insertvalue rn sm TypeNum.d1
-
-
-instance (Integral n, ToSize n, MultiValue.Comparison n) => C (Shifted n) where
-   type Index (Shifted n) = n
-   intersectCode =
-      MultiValue.modifyF2 (singletonShifted atom) (singletonShifted atom) $
-            \(Shifted offsetN lenN) (Shifted offsetM lenM) -> do
-         offset <- MultiValue.max offsetN offsetM
-         endN <- MultiValue.add offsetN lenN
-         endM <- MultiValue.add offsetM lenM
-         end <- MultiValue.min endN endM
-         Shifted offset <$> MultiValue.sub end offset
-   sizeCode = toSize . shiftedSize . unzipShifted
-   size (Shifted _offset len) = fromIntegral len
-   flattenIndexRec shapeValue i =
-      case unzipShifted shapeValue of
-         Shifted offset len ->
-            Monad.lift2 (,) (toSize len) (toSize =<< MultiValue.sub i offset)
-   iterator rngValue =
-      case MultiValue.decompose (singletonShifted atom) rngValue of
-         Shifted from len ->
-            IterMV.take len $ Iter.iterate MultiValue.inc from
-
-
-
-instance (C n, C m) => C (n,m) where
-   type Index (n,m) = (Index n, Index m)
-   intersectCode a b =
-      case (MultiValue.unzip a, MultiValue.unzip b) of
-         ((an,am), (bn,bm)) ->
-            Monad.lift2 MultiValue.zip
-               (intersectCode an bn)
-               (intersectCode am bm)
-   sizeCode nm =
-      case MultiValue.unzip nm of
-         (n,m) -> Monad.liftJoin2 A.mul (sizeCode n) (sizeCode m)
-   size (n,m) = size n * size m
-   flattenIndexRec nm ij =
-      case (MultiValue.unzip nm, MultiValue.unzip ij) of
-         ((n,m), (i,j)) -> do
-            (ns, il) <- flattenIndexRec n i
-            (ms, jl) <- flattenIndexRec m j
-            Monad.lift2 (,)
-               (A.mul ns ms)
-               (A.add jl =<< A.mul ms il)
-   iterator nm =
-      case MultiValue.unzip nm of
-         (n,m) ->
-            uncurry MultiValue.zip <$>
-            Iter.cartesian (iterator n) (iterator m)
-   loop code nm =
-      case MultiValue.unzip nm of
-         (n,m) -> loop (\i -> loop (\j -> code (MultiValue.zip i j)) m) n
-
-instance (C n, C m, C l) => C (n,m,l) where
-   type Index (n,m,l) = (Index n, Index m, Index l)
-   intersectCode a b =
-      case (MultiValue.unzip3 a, MultiValue.unzip3 b) of
-         ((ai,aj,ak), (bi,bj,bk)) ->
-            Monad.lift3 MultiValue.zip3
-               (intersectCode ai bi)
-               (intersectCode aj bj)
-               (intersectCode ak bk)
-   sizeCode nml =
-      case MultiValue.unzip3 nml of
-         (n,m,l) ->
-            Monad.liftJoin2 A.mul (sizeCode n) $
-            Monad.liftJoin2 A.mul (sizeCode m) (sizeCode l)
-   size (n,m,l) = size n * size m * size l
-   flattenIndexRec nml ijk =
-      case (MultiValue.unzip3 nml, MultiValue.unzip3 ijk) of
-         ((n,m,l), (i,j,k)) -> do
-            (ns, il) <- flattenIndexRec n i
-            (ms, jl) <- flattenIndexRec m j
-            x0 <- A.add jl =<< A.mul ms il
-            (ls, kl) <- flattenIndexRec l k
-            x1 <- A.add kl =<< A.mul ls x0
-            sz <- A.mul ns =<< A.mul ms ls
-            return (sz, x1)
-   iterator nml =
-      case MultiValue.unzip3 nml of
-         (n,m,l) ->
-            fmap (\(a,(b,c)) -> MultiValue.zip3 a b c) $
-            Iter.cartesian (iterator n) $
-            Iter.cartesian (iterator m) (iterator l)
-   loop code nml =
-      case MultiValue.unzip3 nml of
-         (n,m,l) ->
-            loop (\i -> loop (\j -> loop (\k ->
-               code (MultiValue.zip3 i j k))
-            l) m) n
diff --git a/src/Data/Array/Knead/Simple/Fold.hs b/src/Data/Array/Knead/Simple/Fold.hs
--- a/src/Data/Array/Knead/Simple/Fold.hs
+++ b/src/Data/Array/Knead/Simple/Fold.hs
@@ -22,7 +22,7 @@
 
 import qualified Data.Array.Knead.Shape.Cubic as Linear
 import qualified Data.Array.Knead.Shape.Cubic.Int as IndexInt
-import qualified Data.Array.Knead.Shape.Nested as Shape
+import qualified Data.Array.Knead.Shape as Shape
 import qualified Data.Array.Knead.Expression as Expr
 import Data.Array.Knead.Expression (Exp, unExp, )
 import Data.Array.Knead.Shape.Cubic ((#:.), (:.)((:.)), )
@@ -74,7 +74,8 @@
    (Val (Linear.Index (sh :. IndexInt.Int)) -> Code r a) ->
    (Val (Linear.Index sh) -> Code r a)
 fold1CodeLinear f nc code ix =
-   Core.fold1Code f (IndexInt.decons nc)
+   Core.fold1Code f
+      (Expr.lift1 (MultiValue.compose . Shape.ZeroBased) $ IndexInt.decons nc)
       (\j -> code (ix #:. IndexInt.cons j))
 
 fold ::
diff --git a/src/Data/Array/Knead/Simple/Physical.hs b/src/Data/Array/Knead/Simple/Physical.hs
--- a/src/Data/Array/Knead/Simple/Physical.hs
+++ b/src/Data/Array/Knead/Simple/Physical.hs
@@ -2,7 +2,8 @@
 {-# LANGUAGE Rank2Types #-}
 {-# LANGUAGE ForeignFunctionInterface #-}
 module Data.Array.Knead.Simple.Physical (
-   Array(Array, shape, buffer), -- data constructor intended for PhysicalParameterized
+   Array,
+   shape,
    toList,
    fromList,
    vectorFromList,
@@ -17,12 +18,17 @@
 
 import qualified Data.Array.Knead.Simple.PhysicalPrivate as Priv
 import qualified Data.Array.Knead.Simple.Private as Sym
-import qualified Data.Array.Knead.Shape.Nested as Shape
+import qualified Data.Array.Knead.Shape as Shape
 import qualified Data.Array.Knead.Expression as Expr
 import qualified Data.Array.Knead.Code as Code
 import Data.Array.Knead.Expression (Exp, unExp, )
 import Data.Array.Knead.Code (getElementPtr, compile, )
 
+import qualified Data.Array.Comfort.Storable.Mutable.Unchecked as MutArray
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
+import qualified Data.Array.Comfort.Shape as ComfortShape
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
+
 import qualified LLVM.Extra.Multi.Value.Memory as MultiValueMemory
 import qualified LLVM.Extra.Multi.Value as MultiValue
 import qualified LLVM.Extra.Arithmetic as A
@@ -31,52 +37,36 @@
 
 import qualified LLVM.Core as LLVM
 
-import Foreign.Marshal.Array (pokeArray, peekArray, )
 import Foreign.Marshal.Alloc (alloca, )
 import Foreign.Storable (Storable, peek, )
-import Foreign.ForeignPtr (ForeignPtr, withForeignPtr, mallocForeignPtrArray, )
+import Foreign.ForeignPtr (withForeignPtr, mallocForeignPtrArray, )
 import Foreign.Ptr (FunPtr, Ptr, )
 
-import Control.Monad.HT (void, )
+import Control.Monad.HT (void, (<=<), )
 import Control.Applicative (liftA2, (<$>), )
 
 import Prelude hiding (scanl1, )
 
 
-data Array sh a =
-   Array {
-      shape :: sh,
-      buffer :: ForeignPtr a
-   }
-
+shape :: Array sh a -> sh
+shape = Array.shape
 
 toList ::
    (Shape.C sh, Storable a) =>
    Array sh a -> IO [a]
-toList (Array sh fptr) =
-   withForeignPtr fptr $ peekArray (Shape.size sh)
+toList = MutArray.toList <=< MutArray.unsafeThaw
 
 fromList ::
    (Shape.C sh, Storable a) =>
    sh -> [a] -> IO (Array sh a)
-fromList sh xs = do
-   let size = Shape.size sh
-   fptr <- mallocForeignPtrArray size
-   withForeignPtr fptr $
-      \ptr ->
-         pokeArray ptr $
-         take size $
-         xs ++ repeat (error "Array.Knead.Physical.fromList: list too short for shape")
-   return (Array sh fptr)
+fromList sh = MutArray.unsafeFreeze <=< MutArray.fromList sh
 
 vectorFromList ::
-   (Shape.C sh, Num sh, Storable a) =>
-   [a] -> IO (Array sh a)
-vectorFromList xs = do
-   let size = length xs
-   fptr <- mallocForeignPtrArray size
-   withForeignPtr fptr $ flip pokeArray xs
-   return (Array (fromIntegral size) fptr)
+   (Num n, Storable a) =>
+   [a] -> IO (Array (ComfortShape.ZeroBased n) a)
+vectorFromList xs =
+   Array.mapShape (\(Shape.ZeroBased n) -> Shape.ZeroBased $ fromIntegral n) <$>
+   (MutArray.unsafeFreeze =<< MutArray.vectorFromList xs)
 
 
 {- |
@@ -122,7 +112,7 @@
             (Code.createFunction callShaper "shape" $ \ptr -> do
                sh <- unExp esh
                MultiValueMemory.store sh ptr
-               Shape.sizeCode sh >>= LLVM.ret)
+               Shape.size sh >>= LLVM.ret)
             (Code.createFunction callRenderer "fill"
                (\paramPtr arrayPtr -> code paramPtr arrayPtr >> LLVM.ret ()))
       let lshptr = MultiValueMemory.castStructPtr shptr
diff --git a/src/Data/Array/Knead/Simple/PhysicalPrivate.hs b/src/Data/Array/Knead/Simple/PhysicalPrivate.hs
--- a/src/Data/Array/Knead/Simple/PhysicalPrivate.hs
+++ b/src/Data/Array/Knead/Simple/PhysicalPrivate.hs
@@ -2,7 +2,7 @@
 module Data.Array.Knead.Simple.PhysicalPrivate where
 
 import qualified Data.Array.Knead.Simple.Private as Sym
-import qualified Data.Array.Knead.Shape.Nested as Shape
+import qualified Data.Array.Knead.Shape as Shape
 import qualified Data.Array.Knead.Expression as Expr
 import Data.Array.Knead.Expression (Exp, unExp)
 import Data.Array.Knead.Code (getElementPtr)
diff --git a/src/Data/Array/Knead/Simple/Private.hs b/src/Data/Array/Knead/Simple/Private.hs
--- a/src/Data/Array/Knead/Simple/Private.hs
+++ b/src/Data/Array/Knead/Simple/Private.hs
@@ -2,7 +2,7 @@
 {-# LANGUAGE TypeFamilies #-}
 module Data.Array.Knead.Simple.Private where
 
-import qualified Data.Array.Knead.Shape.Nested as Shape
+import qualified Data.Array.Knead.Shape as Shape
 import qualified Data.Array.Knead.Expression as Expr
 import Data.Array.Knead.Expression (Exp(Exp), )
 
diff --git a/src/Data/Array/Knead/Simple/ShapeDependent.hs b/src/Data/Array/Knead/Simple/ShapeDependent.hs
--- a/src/Data/Array/Knead/Simple/ShapeDependent.hs
+++ b/src/Data/Array/Knead/Simple/ShapeDependent.hs
@@ -4,7 +4,7 @@
 import qualified Data.Array.Knead.Simple.Private as Core
 import Data.Array.Knead.Simple.Private (Array(Array), )
 
-import qualified Data.Array.Knead.Shape.Nested as Shape
+import qualified Data.Array.Knead.Shape as Shape
 import qualified Data.Array.Knead.Expression as Expr
 import Data.Array.Knead.Expression (Exp, )
 
diff --git a/src/Data/Array/Knead/Simple/Slice.hs b/src/Data/Array/Knead/Simple/Slice.hs
--- a/src/Data/Array/Knead/Simple/Slice.hs
+++ b/src/Data/Array/Knead/Simple/Slice.hs
@@ -70,7 +70,7 @@
 import qualified Data.Array.Knead.Simple.Private as Core
 
 import qualified Data.Array.Knead.Shape.Cubic as Linear
-import qualified Data.Array.Knead.Shape.Nested as Shape
+import qualified Data.Array.Knead.Shape as Shape
 import qualified Data.Array.Knead.Expression as Expr
 import Data.Array.Knead.Shape.Cubic ((#:.), (:.)((:.)), )
 import Data.Array.Knead.Expression (Exp, )
diff --git a/src/Data/Array/Knead/Simple/Symbolic.hs b/src/Data/Array/Knead/Simple/Symbolic.hs
--- a/src/Data/Array/Knead/Simple/Symbolic.hs
+++ b/src/Data/Array/Knead/Simple/Symbolic.hs
@@ -31,7 +31,7 @@
 import qualified Data.Array.Knead.Simple.Private as Core
 import Data.Array.Knead.Simple.Private (Array, shape, gather, )
 
-import qualified Data.Array.Knead.Shape.Nested as Shape
+import qualified Data.Array.Knead.Shape as Shape
 import qualified Data.Array.Knead.Expression as Expr
 import Data.Array.Knead.Expression (Exp, )