diff --git a/Data/Array/Repa/ExtShape.hs b/Data/Array/Repa/ExtShape.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Repa/ExtShape.hs
@@ -0,0 +1,41 @@
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE TypeOperators #-}
+
+-- | Additional functions on shapes
+
+module Data.Array.Repa.ExtShape where
+
+import Data.Array.Repa.Index
+import Data.Array.Repa.Shape
+
+
+
+-- | A number of additional operations that are useful together with
+-- 'PrimitiveArray's.
+
+class ExtShape sh where
+
+  -- | subtract the right coordinates from the left. Does not check if the
+  -- resulting shape make sense.
+
+  subDim :: sh -> sh -> sh
+
+  -- | Given an index and an extend, return a list of all indices. For
+  -- @rangeList (Z:.3) (Z:.2)@ this returns @[(Z:.3), (Z:.4), (Z:.5)]@.
+
+  rangeList :: sh -> sh -> [sh]
+
+
+
+instance ExtShape Z where
+  subDim _ _ = Z
+  {-# INLINE subDim #-}
+  rangeList _ _ = [Z]
+  {-# INLINE rangeList #-}
+
+instance ExtShape sh => ExtShape (sh:.Int) where
+  subDim (sh1:.n1) (sh2:.n2) = subDim sh1 sh2 :. (n1-n2)
+  {-# INLINE subDim #-}
+  rangeList (sh1:.n1) (sh2:.n2) = [sh:.n | sh <- rangeList sh1 sh2, n <- [n1 .. (n1+n2) ] ]
+  {-# INLINE rangeList #-}
+
diff --git a/Data/Array/Repa/Index.hs b/Data/Array/Repa/Index.hs
deleted file mode 100644
--- a/Data/Array/Repa/Index.hs
+++ /dev/null
@@ -1,171 +0,0 @@
-{-# LANGUAGE TypeOperators, FlexibleInstances, ScopedTypeVariables #-}
-
--- | Index types.
-module Data.Array.Repa.Index
-	(
-	-- * Index types
-	  Z	(..)
-	, (:.)	(..)
-
-	-- * Common dimensions.
-	, DIM0, DIM1, DIM2, DIM3, DIM4, DIM5
-        ,       ix1,  ix2,  ix3,  ix4,  ix5)
-where
-import Data.Array.Repa.Shape
-import GHC.Base 		(quotInt, remInt)
-
-stage	= "Data.Array.Repa.Index"
-
--- | An index of dimension zero
-data Z	= Z
-	deriving (Show, Read, Eq, Ord)
-
--- | Our index type, used for both shapes and indices.
-infixl 3 :.
-data tail :. head
-	= !tail :. !head
-	deriving (Show, Read, Eq, Ord)
-
--- Common dimensions
-type DIM0	= Z
-type DIM1	= DIM0 :. Int
-type DIM2	= DIM1 :. Int
-type DIM3	= DIM2 :. Int
-type DIM4	= DIM3 :. Int
-type DIM5	= DIM4 :. Int
-
-
--- | Helper for index construction.
---
---   Use this instead of explicit constructors like @(Z :. (x :: Int))@.
---   The this is sometimes needed to ensure that 'x' is constrained to 
---   be in @Int@.
-ix1 :: Int -> DIM1
-ix1 x = Z :. x
-{-# INLINE ix1 #-}
-
-ix2 :: Int -> Int -> DIM2
-ix2 y x = Z :. y :. x
-{-# INLINE ix2 #-}
-
-ix3 :: Int -> Int -> Int -> DIM3
-ix3 z y x = Z :. z :. y :. x
-{-# INLINE ix3 #-}
-
-ix4 :: Int -> Int -> Int -> Int -> DIM4
-ix4 a z y x = Z :. a :. z :. y :. x
-{-# INLINE ix4 #-}
-
-ix5 :: Int -> Int -> Int -> Int -> Int -> DIM5
-ix5 b a z y x = Z :. b :. a :. z :. y :. x
-{-# INLINE ix5 #-}
-
-
--- Shape ----------------------------------------------------------------------
-instance Shape Z where
-	{-# INLINE [1] rank #-}
-	rank _			= 0
-
-	{-# INLINE [1] zeroDim #-}
-	zeroDim		 	= Z
-
-	{-# INLINE [1] unitDim #-}
-	unitDim			= Z
-
-	{-# INLINE [1] intersectDim #-}
-	intersectDim _ _	= Z
-
-	{-# INLINE [1] addDim #-}
-	addDim _ _		= Z
-
-	{-# INLINE [1] size #-}
-	size _			= 1
-
-	{-# INLINE [1] sizeIsValid #-}
-	sizeIsValid _		= True
-
-
-	{-# INLINE [1] toIndex #-}
-	toIndex _ _		= 0
-
-	{-# INLINE [1] fromIndex #-}
-	fromIndex _ _		= Z
-
-
-	{-# INLINE [1] inShapeRange #-}
-	inShapeRange Z Z Z	= True
-
-        {-# NOINLINE listOfShape #-}
-	listOfShape _		= []
-
-        {-# NOINLINE shapeOfList #-}
-	shapeOfList []		= Z
-	shapeOfList _		= error $ stage ++ ".fromList: non-empty list when converting to Z."
-
-	{-# INLINE deepSeq #-}
-	deepSeq Z x		= x
-
-
-instance Shape sh => Shape (sh :. Int) where
-	{-# INLINE [1] rank #-}
-	rank   (sh  :. _)
-		= rank sh + 1
-
-	{-# INLINE [1] zeroDim #-}
-	zeroDim = zeroDim :. 0
-
-	{-# INLINE [1] unitDim #-}
-	unitDim = unitDim :. 1
-
-	{-# INLINE [1] intersectDim #-}
-	intersectDim (sh1 :. n1) (sh2 :. n2)
-		= (intersectDim sh1 sh2 :. (min n1 n2))
-
-	{-# INLINE [1] addDim #-}
-	addDim (sh1 :. n1) (sh2 :. n2)
-		= addDim sh1 sh2 :. (n1 + n2)
-
-	{-# INLINE [1] size #-}
-	size  (sh1 :. n)
-		= size sh1 * n
-
-	{-# INLINE [1] sizeIsValid #-}
-	sizeIsValid (sh1 :. n)
-		| size sh1 > 0
-		= n <= maxBound `div` size sh1
-
-		| otherwise
-		= False
-
-	{-# INLINE [1] toIndex #-}
-	toIndex (sh1 :. sh2) (sh1' :. sh2')
-		= toIndex sh1 sh1' * sh2 + sh2'
-
-	{-# INLINE [1] fromIndex #-}
-        fromIndex (ds :. d) n
-                = fromIndex ds (n `quotInt` d) :. r
-                where
-                -- If we assume that the index is in range, there is no point
-                -- in computing the remainder for the highest dimension since
-                -- n < d must hold. This saves one remInt per element access which
-                -- is quite a big deal.
-                r       | rank ds == 0  = n
-                        | otherwise     = n `remInt` d
-
-	{-# INLINE [1] inShapeRange #-}
-	inShapeRange (zs :. z) (sh1 :. n1) (sh2 :. n2)
-		= (n2 >= z) && (n2 < n1) && (inShapeRange zs sh1 sh2)
-
-        {-# NOINLINE listOfShape #-}
-       	listOfShape (sh :. n)
-	 = n : listOfShape sh
-
-        {-# NOINLINE shapeOfList #-}
-	shapeOfList xx
-	 = case xx of
-		[]	-> error $ stage ++ ".toList: empty list when converting to  (_ :. Int)"
-		x:xs	-> shapeOfList xs :. x
-
-	{-# INLINE deepSeq #-}
-	deepSeq (sh :. n) x = deepSeq sh (n `seq` x)
-
diff --git a/Data/Array/Repa/Index/Point.hs b/Data/Array/Repa/Index/Point.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Repa/Index/Point.hs
@@ -0,0 +1,115 @@
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE TypeOperators #-}
+
+-- | A Repa-compatible index of points in multi-dimensional space. A point
+-- represents the index of a left- or right-linear grammar.
+
+module Data.Array.Repa.Index.Point where
+
+import Control.Applicative
+import Control.DeepSeq
+import Data.Array.Repa.Index
+import Data.Array.Repa.Shape
+import GHC.Base (quotInt, remInt)
+import Test.QuickCheck
+import Test.QuickCheck.All
+
+import Data.Array.Repa.ExtShape
+
+
+
+stage = "Data.Array.Repa.Index.Point"
+
+-- |
+
+newtype Point = Point Int
+  deriving (Eq,Ord,Show)
+
+-- |
+
+instance Shape sh => Shape (sh :. Point) where
+  {-# INLINE [1] rank #-}
+  rank   (sh  :. _)
+    = rank sh + 1
+
+  {-# INLINE [1] zeroDim #-}
+  zeroDim = zeroDim :. Point 0
+
+  {-# INLINE [1] unitDim #-}
+  unitDim = unitDim :. Point 1
+
+  {-# INLINE [1] intersectDim #-}
+  intersectDim (sh1 :. Point n1) (sh2 :. Point n2)
+    = (intersectDim sh1 sh2 :. Point (min n1 n2))
+
+  {-# INLINE [1] addDim #-}
+  addDim (sh1 :. Point n1) (sh2 :. Point n2)
+    = addDim sh1 sh2 :. Point (n1 + n2)
+
+  {-# INLINE [1] size #-}
+  size  (sh1 :. Point n)
+    = size sh1 * n
+
+  {-# INLINE [1] sizeIsValid #-}
+  sizeIsValid (sh1 :. Point n)
+    | size sh1 > 0
+    = n <= maxBound `div` size sh1
+
+    | otherwise
+    = False
+
+  {-# INLINE [1] toIndex #-}
+  toIndex (sh1 :. Point sh2) (sh1' :. Point sh2')
+    = toIndex sh1 sh1' * sh2 + sh2'
+
+  {-# INLINE [1] fromIndex #-}
+  fromIndex (ds :. Point d) n
+          = fromIndex ds (n `quotInt` d) :. Point r
+          where
+          -- If we assume that the index is in range, there is no point
+          -- in computing the remainder for the highest dimension since
+          -- n < d must hold. This saves one remInt per element access which
+          -- is quite a big deal.
+          r       | rank ds == 0  = n
+                  | otherwise     = n `remInt` d
+
+  {-# INLINE [1] inShapeRange #-}
+  inShapeRange (zs :. Point z) (sh1 :. Point n1) (sh2 :. Point n2)
+    = (n2 >= z) && (n2 < n1) && (inShapeRange zs sh1 sh2)
+
+  {-# NOINLINE listOfShape #-}
+  listOfShape (sh :. Point n)
+   = n : listOfShape sh
+
+  {-# NOINLINE shapeOfList #-}
+  shapeOfList xx
+   = case xx of
+    []  -> error $ stage ++ ".toList: empty list when converting to  (_ :. Point)"
+    x:xs  -> shapeOfList xs :. Point x
+
+  {-# INLINE deepSeq #-}
+  deepSeq (sh :. Point n) x = deepSeq sh (n `seq` x)
+
+-- |
+
+instance ExtShape sh => ExtShape (sh:.Point) where
+  subDim (sh1:.Point n1) (sh2:.Point n2) = subDim sh1 sh2 :. Point (n1-n2)
+  {-# INLINE subDim #-}
+  rangeList (sh1:.Point n1) (sh2:.Point n2) = [sh:.Point n | sh <- rangeList sh1 sh2, n <- [n1 .. (n1+n2) ] ]
+  {-# INLINE rangeList #-}
+
+instance NFData Point where
+  rnf (Point i) = rnf i
+
+-- |
+
+instance Arbitrary Point where
+  arbitrary = Point <$> choose (0,100)
+  shrink (Point i)
+    | i>0       = [Point (i-1)]
+    | otherwise = []
+
+instance Arbitrary z => Arbitrary (z:.Point) where
+  arbitrary = (:.) <$> arbitrary <*> arbitrary
+  shrink (z:.s) = (:.) <$> shrink z <*> shrink s
+
diff --git a/Data/Array/Repa/Index/Subword.hs b/Data/Array/Repa/Index/Subword.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Repa/Index/Subword.hs
@@ -0,0 +1,186 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE TemplateHaskell #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE TypeOperators #-}
+
+-- | Subwords span upper triangular tables. A subword (i,j) is legal iff i<=j.
+--
+-- NOTE Using more complicated shapes has a number of benefits. We don't need
+-- to specify triangular or rectangular tables anymore. A rectangular
+-- one-dimensional table with a subword as shape actually /does/ create space
+-- as required for subwords.
+--
+-- TODO subword indexing is currently hardcoded to be zero-based. See
+-- 'subwordIndex'.
+--
+-- TODO consider replacing (`quot` 2) with (`shiftR` 1)
+--
+-- TODO all the QuickCheck stuff is missing
+
+module Data.Array.Repa.Index.Subword where
+
+import Control.Applicative
+import Control.DeepSeq
+import Data.Array.Repa.Index
+import Data.Array.Repa.Shape
+import GHC.Base (quotInt, remInt)
+import Test.QuickCheck
+import Test.QuickCheck.All
+import qualified Data.Vector.Unboxed as VU
+import Data.Vector.Unboxed.Deriving
+
+import Data.Array.Repa.ExtShape
+
+
+
+stage = "Data.Array.Repa.Index.Subword"
+
+-- | A subword wraps a simple pair.
+--
+-- Subwords always yield the upper-triangular part of a rect-angular array.
+-- This gives the quite curious effect that (0,N) points to the ``largest''
+-- index, while (0,0) and (N,N) both point to the smallest. We do, however, use
+-- (0,0) as the smallest as (0,k) gives successively smaller upper triangular
+-- parts.
+
+newtype Subword = Subword (Int:.Int)
+  deriving (Eq,Ord,Show)
+
+derivingUnbox "Subword"
+  [t| Subword -> (Int,Int) |]
+  [| \ (Subword (i:.j)) -> (i,j) |]
+  [| \ (i,j) -> Subword (i:.j) |]
+
+subword :: Int -> Int -> Subword
+subword i j = Subword (i:.j)
+{-# INLINE subword #-}
+
+-- | triangular numbers
+--
+-- A000217
+
+triangularNumber :: Int -> Int
+triangularNumber x = (x * (x+1)) `quot` 2
+{-# INLINE triangularNumber #-}
+
+-- | Size of an upper triangle starting at 'i' and ending at 'j'. "(0,N)" what
+-- be the normal thing to use.
+
+upperTri :: Subword -> Int
+upperTri (Subword (i:.j)) = triangularNumber $ j-i
+{-# INLINE upperTri #-}
+
+-- | Subword indexing. Given the longest subword and the current subword,
+-- calculate a linear index "[0,..]". "(l,n)" in this case means "l"ower bound,
+-- length "n". And "(i,j)" is the normal index.
+--
+-- TODO probably doesn't work right with non-zero base ?!
+
+subwordIndex :: Subword -> Subword -> Int
+subwordIndex (Subword (l:.n)) (Subword (i:.j)) = adr n (i,j) -- - adr n (l,n)
+  where
+    adr n (i,j) = n*i - triangularNumber i + j
+{-# INLINE subwordIndex #-}
+
+subwordFromIndex :: Subword -> Int -> Subword
+subwordFromIndex = error "not implemented"
+{-# INLINE subwordFromIndex #-}
+
+-- | Some weird things are going on here. Adding subwords (i,j) and (k,l)
+-- yields (i+k,j+l). Normally i==k==0 when calculating space requirements. If
+-- you have a subword (3,10) and want the next outer one add (-1,1) and you get
+-- what you want. We make NO(!) check that the final subword contains only
+-- non-negative indices.
+
+instance Shape sh => Shape (sh :. Subword) where
+  {-# INLINE [1] rank #-}
+  rank   (sh  :. _)
+    = rank sh + 1
+
+  {-# INLINE [1] zeroDim #-}
+  zeroDim = zeroDim :. Subword (0:.0)
+
+  {-# INLINE [1] unitDim #-}
+  unitDim = unitDim :. Subword (0:.1)
+
+  {-# INLINE [1] intersectDim #-}
+  intersectDim (sh1 :. Subword (i:.j)) (sh2 :. Subword (k:.l))
+    = (intersectDim sh1 sh2 :. Subword (max i k :. min j l))
+
+  {-# INLINE [1] addDim #-}
+  addDim (sh1 :. Subword (i:.j)) (sh2 :. Subword (k:.l))
+    = addDim sh1 sh2 :. Subword (i+k:.j+l)
+
+  {-# INLINE [1] size #-}
+  size  (sh1 :. sw) = size sh1 * upperTri sw
+
+  {-# INLINE [1] sizeIsValid #-}
+  sizeIsValid (sh1 :. Subword (i:.j))
+    | size sh1 > 0
+    = i>=0 && i<=j && j <= maxBound `div` size sh1
+    | otherwise
+    = False
+
+  {-# INLINE [1] toIndex #-}
+  toIndex (sh1 :. sh2) (sh1' :. sh2')
+    = toIndex sh1 sh1' * upperTri sh2 + subwordIndex sh2 sh2'
+
+  {-# INLINE [1] fromIndex #-}
+  fromIndex (ds :. d) n  = undefined -- fromIndex ds (n `quotInt` d) :. r
+    where
+      r = subwordFromIndex d n
+    -- If we assume that the index is in range, there is no point
+    -- in computing the remainder for the highest dimension since
+    -- n < d must hold. This saves one remInt per element access which
+    -- is quite a big deal.
+    {-
+    r       | rank ds == 0  = n
+            | otherwise     = n `remInt` d -}
+
+  -- | TODO fix for lower bounds check!
+  {-# INLINE [1] inShapeRange #-}
+  inShapeRange (zs :. Subword (_:._)) (sh1 :. Subword (l:.n)) (sh2 :. Subword (i:.j))
+    = i<=j && l<=i && j<n && (inShapeRange zs sh1 sh2)
+
+  {-# NOINLINE listOfShape #-}
+  listOfShape (sh :. Subword (i:.j)) = i : j : listOfShape sh
+
+  {-# NOINLINE shapeOfList #-}
+  shapeOfList xx
+   = case xx of
+    []     -> error $ stage ++ ".toList: empty list when converting to  (_ :. Int)"
+    [x]    -> error $ stage ++ ".toList: only single element remaining!"
+    i:j:xs -> shapeOfList xs :. Subword (i:.j)
+
+  {-# INLINE deepSeq #-}
+  deepSeq (sh :. n) x = deepSeq sh (n `seq` x)
+
+-- |
+
+instance ExtShape sh => ExtShape (sh:.Subword) where
+  subDim (sh1:.Subword (i:.j)) (sh2:.Subword (k:.l)) = subDim sh1 sh2 :. Subword (i-k:.j-l)
+  {-# INLINE subDim #-}
+  rangeList (sh1:.Subword (i:.j)) (sh2:.Subword (k:.l)) = error "not implemented" -- [sh:.Subword (m,n) | sh <- rangeList sh1 sh2, m <- [i .. [i+k], n <- [ n <- [n1 .. (n1+n2) ] ]
+  {-# INLINE rangeList #-}
+
+-- |
+
+instance NFData Subword where
+  rnf (Subword (i:.j)) = i `seq` rnf j
+
+-- |
+
+instance Arbitrary Subword where
+  arbitrary = do
+    a <- choose (0,100)
+    b <- choose (0,100)
+    return $ Subword (min a b :. max a b)
+  shrink (Subword (i:.j))
+    | i<j       = [Subword (i:.j-1)]
+    | otherwise = []
+
+instance Arbitrary z => Arbitrary (z:.Subword) where
+  arbitrary = (:.) <$> arbitrary <*> arbitrary
+  shrink (z:.s) = (:.) <$> shrink z <*> shrink s
+
diff --git a/Data/Array/Repa/Shape.hs b/Data/Array/Repa/Shape.hs
deleted file mode 100644
--- a/Data/Array/Repa/Shape.hs
+++ /dev/null
@@ -1,82 +0,0 @@
-{-# LANGUAGE RankNTypes #-}
-
--- | Class of types that can be used as array shapes and indices.
-module Data.Array.Repa.Shape
-	( Shape(..)
-        , inShape
-        , showShape )
-where
-
--- Shape ----------------------------------------------------------------------
--- | Class of types that can be used as array shapes and indices.
-class Eq sh => Shape sh where
-
-	-- | Get the number of dimensions in a shape.
-	rank	:: sh -> Int
-
-	-- | The shape of an array of size zero, with a particular dimensionality.
-	zeroDim	:: sh
-
-	-- | The shape of an array with size one, with a particular dimensionality.
-	unitDim :: sh
-
-	-- | Compute the intersection of two shapes.
-	intersectDim :: sh -> sh -> sh
-
-	-- | Add the coordinates of two shapes componentwise
-	addDim  :: sh -> sh -> sh
-
-	-- | Get the total number of elements in an array with this shape.
-	size	:: sh -> Int
-
-	-- | Check whether this shape is small enough so that its flat
-	--	indices an be represented as `Int`. If this returns `False` then your
-	--	array is too big. Mostly used for writing QuickCheck tests.
-	sizeIsValid :: sh -> Bool
-
-
-	-- | Convert an index into its equivalent flat, linear, row-major version.
-	toIndex :: sh	-- ^ Shape of the array.
-		-> sh 	-- ^ Index into the array.
-		-> Int
-
-	-- | Inverse of `toIndex`.
-	fromIndex
-		:: sh 	-- ^ Shape of the array.
-		-> Int 	-- ^ Index into linear representation.
-		-> sh
-
-	-- | Check whether an index is within a given shape.
-	inShapeRange
-		:: sh 	-- ^ Start index for range.
-		-> sh 	-- ^ Final index for range.
-		-> sh 	-- ^ Index to check for.
-		-> Bool
-
-	-- | Convert a shape into its list of dimensions.
-	listOfShape	:: sh -> [Int]
-
-	-- | Convert a list of dimensions to a shape
-	shapeOfList	:: [Int] -> sh
-
-	-- | Ensure that a shape is completely evaluated.
-	infixr 0 `deepSeq`
-	deepSeq :: sh -> a -> a
-
-
--- | Check whether an index is a part of a given shape.
-inShape :: forall sh
-	.  Shape sh
-	=> sh 		-- ^ Shape of the array.
-	-> sh		-- ^ Index.
-	-> Bool
-
-{-# INLINE inShape #-}
-inShape sh ix
-	= inShapeRange zeroDim sh ix
-
-
--- | Nicely format a shape as a string
-showShape :: Shape sh => sh -> String
-showShape = foldr (\sh str -> str ++ " :. " ++ show sh) "Z" . listOfShape
-
diff --git a/Data/ExtShape.hs b/Data/ExtShape.hs
deleted file mode 100644
--- a/Data/ExtShape.hs
+++ /dev/null
@@ -1,41 +0,0 @@
-{-# LANGUAGE FlexibleInstances #-}
-{-# LANGUAGE TypeOperators #-}
-
--- | Additional functions on shapes
-
-module Data.ExtShape where
-
-import Data.Array.Repa.Index
-import Data.Array.Repa.Shape
-
-
-
--- | A number of additional operations that are useful together with
--- 'PrimitiveArray's.
-
-class (Eq sh, Shape sh) => ExtShape sh where
-
-  -- | subtract the right coordinates from the left. Does not check if the
-  -- resulting shape make sense.
-
-  subDim :: sh -> sh -> sh
-
-  -- | Given an index and an extend, return a list of all indices. For
-  -- @rangeList (Z:.3) (Z:.2)@ this returns @[(Z:.3), (Z:.4), (Z:.5)]@.
-
-  rangeList :: sh -> sh -> [sh]
-
-
-
-instance ExtShape Z where
-  subDim _ _ = Z
-  {-# INLINE subDim #-}
-  rangeList _ _ = [Z]
-  {-# INLINE rangeList #-}
-
-instance ExtShape sh => ExtShape (sh:.Int) where
-  subDim (sh1:.n1) (sh2:.n2) = subDim sh1 sh2 :. (n1-n2)
-  {-# INLINE subDim #-}
-  rangeList (sh1:.n1) (sh2:.n2) = [sh:.n | sh <- rangeList sh1 sh2, n <- [n1 .. (n1+n2) ] ]
-  {-# INLINE rangeList #-}
-
diff --git a/Data/PrimitiveArray.hs b/Data/PrimitiveArray.hs
--- a/Data/PrimitiveArray.hs
+++ b/Data/PrimitiveArray.hs
@@ -6,70 +6,68 @@
 
 -- | Vastly extended primitive arrays. Some basic ideas are now modeled after
 -- the vector package, especially the monadic mutable / pure immutable array
--- system. There are eight flavors of arrays among three axes: mutable/pure +
--- boxed/unboxed + zero-based/lower-bound.
+-- system.
 --
 -- NOTE all operations in MPrimArrayOps and PrimArrayOps are highly unsafe. No
 -- bounds-checking is performed at all.
 
 module Data.PrimitiveArray where
 
+import Control.Exception (assert)
+import Control.Monad
+import Control.Monad.Primitive
+import Control.Monad.ST
 import Data.Array.Repa.Index
 import Data.Array.Repa.Shape
-import Data.Primitive.Types
 import Data.Primitive
-import Control.Monad.ST
-import Control.Monad
-import Control.Monad.Primitive
+import Data.Primitive.Types
+import Prelude as P
 import System.IO.Unsafe
-import Control.Exception (assert)
 
-import Data.ExtShape
+import Data.Array.Repa.ExtShape
 
 
 
--- | The core set of operations for monadic arrays.
+-- | Mutable version of an array.
 
-class (Shape sh, ExtShape sh) => MPrimArrayOps marr sh elm where
+data family MutArr (m :: * -> *) (arr :: *) :: *
 
+
+-- | The core set of operations for monadic arrays.
+class (Shape sh, ExtShape sh) => MPrimArrayOps arr sh elm where
+
   -- | Return the bounds of the array. All bounds are inclusive, as in
   -- @[lb..ub]@
 
-  boundsM :: marr s sh elm -> (sh,sh)
+  boundsM :: MutArr m (arr sh elm) -> (sh,sh)
 
   -- | Given lower and upper bounds and a list of /all/ elements, produce a
   -- mutable array.
 
-  fromListM :: PrimMonad m => sh -> sh -> [elm] -> m (marr (PrimState m) sh elm)
+  fromListM :: PrimMonad m => sh -> sh -> [elm] -> m (MutArr m (arr sh elm))
 
   -- | Creates a new array with the given bounds with each element within the
   -- array being in an undefined state.
 
-  newM :: PrimMonad m => sh -> sh -> m (marr (PrimState m) sh elm)
+  newM :: PrimMonad m => sh -> sh -> m (MutArr m (arr sh elm))
 
   -- | Creates a new array with all elements being equal to 'elm'.
 
-  newWithM :: PrimMonad m => sh -> sh -> elm -> m (marr (PrimState m) sh elm)
+  newWithM :: PrimMonad m => sh -> sh -> elm -> m (MutArr m (arr sh elm))
 
   -- | Reads a single element in the array.
 
-  readM :: PrimMonad m => marr (PrimState m) sh elm -> sh -> m elm
+  readM :: PrimMonad m => MutArr m (arr sh elm) -> sh -> m elm
 
   -- | Writes a single element in the array.
 
-  writeM :: PrimMonad m => marr (PrimState m) sh elm -> sh -> elm -> m ()
-
-
-
--- | Used to connect each immutable array with one mutable array.
-
-type family MutArray (v :: * -> * -> * ) :: * -> * -> * -> *
+  writeM :: PrimMonad m => MutArr m (arr sh elm) -> sh -> elm -> m ()
 
 
 
 -- | The core set of functions on immutable arrays.
 
-class (Shape sh, ExtShape sh, MPrimArrayOps (MutArray arr) sh elm) => PrimArrayOps arr sh elm where
+class (Shape sh, ExtShape sh) => PrimArrayOps arr sh elm where
 
   -- | Returns the bounds of an immutable array, again inclusive bounds: @ [lb..ub] @.
 
@@ -79,15 +77,21 @@
   -- is /O(1)/ and both arrays share the same memory. Do not use the mutable
   -- array afterwards.
 
-  freeze :: PrimMonad m => MutArray arr (PrimState m) sh elm -> m (arr sh elm)
+  freeze :: PrimMonad m => MutArr m (arr sh elm) -> m (arr sh elm)
 
   -- | Extract a single element from the array. Generally unsafe as not
   -- bounds-checking is performed.
 
   index :: arr sh elm -> sh -> elm
 
+class (Shape sh, ExtShape sh) => PrimArrayMap arr sh e e' where
 
+  -- | Map a function over each element, keeping the shape intact.
 
+  map :: (e -> e') -> arr sh e -> arr sh e'
+
+
+
 -- | Infix index operator. Performs minimal bounds-checking using assert in
 -- non-optimized code.
 
@@ -97,7 +101,7 @@
 
 -- | Returns true if the index is valid for the array.
 
-inBoundsM :: MPrimArrayOps marr sh elm => marr s sh elm -> sh -> Bool
+inBoundsM :: (Monad m, MPrimArrayOps arr sh elm) => MutArr m (arr sh elm) -> sh -> Bool
 inBoundsM marr idx = let (lb,ub) = boundsM marr in inShapeRange lb ub idx
 {-# INLINE inBoundsM #-}
 
@@ -111,16 +115,16 @@
 sliceEq :: (Eq elm, PrimArrayOps arr sh elm) => arr sh elm -> sh -> arr sh elm -> sh -> sh -> Bool
 sliceEq arr1 k1 arr2 k2 xtnd = assert ((inBounds arr1 k1) && (inBounds arr2 k2) && (inBounds arr1 $ k1 `addDim` xtnd) && (inBounds arr2 $ k2 `addDim` xtnd)) $ and res where
   res = zipWith (==) xs ys
-  xs = map (index arr1) $ rangeList k1 xtnd
-  ys = map (index arr2) $ rangeList k2 xtnd
+  xs = P.map (index arr1) $ rangeList k1 xtnd
+  ys = P.map (index arr2) $ rangeList k2 xtnd
 {-# INLINE sliceEq #-}
 
 -- | Construct a mutable primitive array from a lower and an upper bound, a
 -- default element, and a list of associations.
 
 fromAssocsM
-  :: (PrimMonad m, MPrimArrayOps marr sh elm)
-  => sh -> sh -> elm -> [(sh,elm)] -> m (marr (PrimState m) sh elm)
+  :: (PrimMonad m, MPrimArrayOps arr sh elm)
+  => sh -> sh -> elm -> [(sh,elm)] -> m (MutArr m (arr sh elm))
 fromAssocsM lb ub def xs = do
   ma <- newWithM lb ub def
   forM_ xs $ \(k,v) -> writeM ma k v
@@ -130,21 +134,21 @@
 -- | Return all associations from an array.
 
 assocs :: PrimArrayOps arr sh elm => arr sh elm -> [(sh,elm)]
-assocs arr = map (\k -> (k,index arr k)) $ rangeList lb (ub `subDim` lb) where
+assocs arr = P.map (\k -> (k,index arr k)) $ rangeList lb (ub `subDim` lb) where
   (lb,ub) = bounds arr
 {-# INLINE assocs #-}
 
 -- | Creates an immutable array from lower and upper bounds and a complete list
 -- of elements.
 
-fromList :: PrimArrayOps arr sh elm => sh -> sh -> [elm] -> arr sh elm
+fromList :: (PrimArrayOps arr sh elm, MPrimArrayOps arr sh elm) => sh -> sh -> [elm] -> arr sh elm
 fromList lb ub xs = runST $ fromListM lb ub xs >>= freeze
 {-# INLINE fromList #-}
 
 -- | Creates an immutable array from lower and upper bounds, a default element,
 -- and a list of associations.
 
-fromAssocs :: PrimArrayOps arr sh elm => sh -> sh -> elm -> [(sh,elm)] -> arr sh elm
+fromAssocs :: (PrimArrayOps arr sh elm, MPrimArrayOps arr sh elm) => sh -> sh -> elm -> [(sh,elm)] -> arr sh elm
 fromAssocs lb ub def xs = runST $ fromAssocsM lb ub def xs >>= freeze
 {-# INLINE fromAssocs #-}
 
@@ -157,6 +161,6 @@
 -- | Returns all elements of an immutable array as a list.
 
 toList :: PrimArrayOps arr sh elm =>  arr sh elm -> [elm]
-toList arr = let (lb,ub) = bounds arr in map ((!) arr) $ rangeList lb $ ub `subDim` lb
+toList arr = let (lb,ub) = bounds arr in P.map ((!) arr) $ rangeList lb $ ub `subDim` lb
 {-# INLINE toList #-}
 
diff --git a/Data/PrimitiveArray/FillTables.hs b/Data/PrimitiveArray/FillTables.hs
new file mode 100644
--- /dev/null
+++ b/Data/PrimitiveArray/FillTables.hs
@@ -0,0 +1,82 @@
+{-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE TypeOperators #-}
+
+-- | Operations to fill primitive arrays. Arrays are combined just like indices
+-- using 'Z' and '(:.)'. This allows filling an unlimited number of tables.
+--
+-- TODO make explicit in which order the tables are filled.
+
+module Data.PrimitiveArray.FillTables where
+
+import Control.Monad.Primitive
+import Data.Array.Repa.Index
+import Data.Array.Repa.Shape
+import Data.Vector.Fusion.Stream.Monadic as S
+
+import Data.PrimitiveArray
+import Data.Array.Repa.Index.Subword
+
+
+
+-- Upper triangular table filling. Right now, only a serial option 'upperTriS'
+-- is available.
+--
+-- TODO Using Repa, 'upperTriP' will soon become available.
+
+class UpperTriS m stack where
+  upperTriS :: stack -> m ()
+
+-- |
+--
+-- TODO Insert check that all extends are the same!
+
+instance
+  ( Monad m
+  , MPrimArrayOps arr Subword e
+  , Stack m Subword (xs:.MutArr m (arr Subword e))
+  ) => UpperTriS m (xs:.MutArr m (arr Subword e)) where
+  upperTriS xs@(_:.x) = do
+    -- TODO missing extends check
+    let (Subword (l:._),Subword (u:._)) = boundsM x
+    S.mapM_ (go xs) $ unfolder l u
+    where
+      -- Write all table values at a certain subword. Note that tables are
+      -- filled left to right in the stack order '(Z:.first:.next:.last)'
+      go xs (i,j) = writeStack xs (Subword (i:.j))
+      {-# INLINE go #-}
+      -- the unfolder steps through the diagonals from the main toward the
+      -- upper-right. It starts a the main-diagonal column and goes to the
+      -- right in each big step '(i==u)' and increases the row by one in each
+      -- small step.
+      unfolder l u = S.unfoldr step (l,l) where
+        step (j,i)
+          | j> u      = Nothing
+          | i==u      = Just ((i,j),(j+1,l))
+          | otherwise = Just ((i,j),(j,i+1))
+        {-# INLINE step #-}
+      {-# INLINE unfolder #-}
+  {-# INLINE upperTriS #-}
+
+
+
+-- | Defines how a single index in a stack of arrays + evaluation functions is
+-- handled.
+
+class Stack m sh xs where
+  writeStack :: xs -> sh -> m ()
+
+instance (Monad m) => Stack m sh Z where
+  writeStack _ _ = return ()
+  {-# INLINE writeStack #-}
+
+instance
+  ( PrimMonad m
+  , Stack m Subword xs
+  , MPrimArrayOps arr Subword e
+  ) => Stack m Subword (xs:.(MutArr m (arr Subword e),(Subword -> m e))) where
+  writeStack (xs:.(x,f)) i = writeStack xs i >> f i >>= writeM x i
+  {-# INLINE writeStack #-}
+
diff --git a/Data/PrimitiveArray/Zero.hs b/Data/PrimitiveArray/Zero.hs
--- a/Data/PrimitiveArray/Zero.hs
+++ b/Data/PrimitiveArray/Zero.hs
@@ -1,9 +1,14 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TemplateHaskell #-}
+{-# LANGUAGE StandaloneDeriving #-}
+{-# LANGUAGE DeriveDataTypeable #-}
+{-# LANGUAGE CPP #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE TypeFamilies #-}
 
--- | Boxed, primitive arrays. A good use-case is to store boxed or unboxed
--- vectors.
+-- | Primitive arrays where the lower index is zero (or the equivalent of zero
+-- for newtypes and enumerations).
 
 module Data.PrimitiveArray.Zero where
 
@@ -11,48 +16,44 @@
 import Data.Array.Repa.Index
 import Data.Array.Repa.Shape
 import Control.Exception (assert)
-import Data.Vector as VU hiding (forM_, length, zipWithM_)
-import Data.Vector.Mutable as VUM hiding (length)
+import qualified Data.Vector.Unboxed as VU hiding (forM_, length, zipWithM_)
+import qualified Data.Vector as V hiding (forM_, length, zipWithM_)
+import qualified Data.Vector.Unboxed.Mutable as VUM hiding (length)
+import Data.Vector.Generic as G hiding (forM_, length, zipWithM_, new)
+import Data.Vector.Generic.Mutable as GM hiding (length)
+import Control.Monad.Primitive (PrimState)
 
-import Data.ExtShape
+import Data.Array.Repa.ExtShape
 import Data.PrimitiveArray
 
 
 
--- | Monadic arrays of boxed type.
-
-data MArr0 s sh elm = MArr0 !sh !(MVector s elm)
-
--- | Immutable arrays of boxed type.
-
-data Arr0 sh elm = Arr0 !sh !(Vector elm)
-
-
+-- * Unboxed, multidimensional arrays.
 
-type instance MutArray Arr0 = MArr0
+data Unboxed sh e = Unboxed !sh !(VU.Vector e)
+  deriving (Read,Show,Eq)
 
--- NOTE inLb, inUb is including bound, while exUb is excluding upper bound.
--- Differentiates between largest included index, first excluded index.
+data instance MutArr m (Unboxed sh e) = MUnboxed !sh (VU.MVector (PrimState m) e)
 
-instance (Shape sh, ExtShape sh) => MPrimArrayOps MArr0 sh elm where
-  boundsM (MArr0 exUb _) = (zeroDim,exUb `subDim` unitDim)
+instance (Shape sh, ExtShape sh, VUM.Unbox elm) => MPrimArrayOps Unboxed sh elm where
+  boundsM (MUnboxed exUb _) = (zeroDim,exUb `subDim` unitDim)
   fromListM inLb inUb xs = do
     ma <- newM inLb inUb
     let exUb = inUb `addDim` unitDim
-    let (MArr0 _ mba) = ma
+    let (MUnboxed _ mba) = ma
     zipWithM_ (\k x -> assert (length xs == size exUb) $ unsafeWrite mba k x) [0.. toIndex exUb inUb] xs
     return ma
   newM inLb inUb = let exUb = inUb `addDim` unitDim in
     unless (inLb == zeroDim) (error "MArr0 lb/=zeroDim") >>
-    MArr0 exUb `liftM` new (size exUb)
+    MUnboxed exUb `liftM` new (size exUb)
   newWithM inLb inUb def = do
     let exUb = inUb `addDim` unitDim
     ma <- newM inLb inUb
-    let (MArr0 _ mba) = ma
+    let (MUnboxed _ mba) = ma
     forM_ [0 .. toIndex exUb inUb] $ \k -> unsafeWrite mba k def
     return ma
-  readM (MArr0 exUb mba) idx = assert (inShape exUb idx) $ unsafeRead mba (toIndex exUb idx)
-  writeM (MArr0 exUb mba) idx elm = assert (inShape exUb idx) $ unsafeWrite mba (toIndex exUb idx) elm
+  readM (MUnboxed exUb mba) idx = assert (inShape exUb idx) $ unsafeRead mba (toIndex exUb idx)
+  writeM (MUnboxed exUb mba) idx elm = assert (inShape exUb idx) $ unsafeWrite mba (toIndex exUb idx) elm
   {-# INLINE boundsM #-}
   {-# INLINE fromListM #-}
   {-# INLINE newM #-}
@@ -60,11 +61,46 @@
   {-# INLINE readM #-}
   {-# INLINE writeM #-}
 
-instance (Shape sh, ExtShape sh) => PrimArrayOps Arr0 sh elm where
-  bounds (Arr0 exUb _) = (zeroDim,exUb `subDim` unitDim)
-  freeze (MArr0 exUb mba) = Arr0 exUb `liftM` unsafeFreeze mba
-  index (Arr0 exUb ba) idx = assert (inShape exUb idx) $ unsafeIndex ba (toIndex exUb idx)
+instance (Shape sh, ExtShape sh, VUM.Unbox elm) => PrimArrayOps Unboxed sh elm where
+  bounds (Unboxed exUb _) = (zeroDim,exUb `subDim` unitDim)
+  freeze (MUnboxed exUb mba) = Unboxed exUb `liftM` unsafeFreeze mba
+  index (Unboxed exUb ba) idx = assert (inShape exUb idx) $ unsafeIndex ba (toIndex exUb idx)
   {-# INLINE bounds #-}
   {-# INLINE freeze #-}
   {-# INLINE index #-}
+
+instance (Shape sh, ExtShape sh, VUM.Unbox e, VUM.Unbox e') => PrimArrayMap Unboxed sh e e' where
+  map f (Unboxed sh xs) = Unboxed sh (VU.map f xs)
+  {-# INLINE map #-}
+
+
+
+-- * Boxed, multidimensional arrays.
+
+-- data Boxed sh e = Boxed !sh !(V.Vector e)
+
+
+{-
+ -
+ - This stuff tells us how to write efficient generics on large data
+ - constructors like the Turner and Vienna ctors.
+ -
+
+import qualified Data.Generics.TH as T
+
+data Unboxed sh e = Unboxed !sh !(VU.Vector e)
+  deriving (Show,Eq,Ord)
+
+data X e = X (Unboxed DIM1 e) (Unboxed DIM1 e)
+  deriving (Show,Eq,Ord)
+
+x :: X Int
+x = X z z where z = (Unboxed (Z:.10) (VU.fromList [ 0 .. 10] ))
+
+pot :: X Int -> X Double
+pot = $( T.thmapT (T.mkTs ['f]) [t| X Int |] ) where
+  f :: Unboxed DIM1 Int -> Unboxed DIM1 Double
+  f (Unboxed sh xs) = Unboxed sh (VU.map fromIntegral xs)
+
+-}
 
diff --git a/Data/PrimitiveArray/Zero/Unboxed.hs b/Data/PrimitiveArray/Zero/Unboxed.hs
deleted file mode 100644
--- a/Data/PrimitiveArray/Zero/Unboxed.hs
+++ /dev/null
@@ -1,71 +0,0 @@
-{-# LANGUAGE ScopedTypeVariables #-}
-{-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE FlexibleInstances #-}
-
--- | Strict, unboxed arrays of primitive type. Uses unboxed vectors internally
--- to provide tuple instances.
-
-module Data.PrimitiveArray.Zero.Unboxed where
-
-import Control.Monad
-import Data.Array.Repa.Index
-import Data.Array.Repa.Shape
-import Control.Exception (assert)
-import Data.Vector.Unboxed as VU hiding (forM_, length, zipWithM_)
-import Data.Vector.Unboxed.Mutable as VUM hiding (length)
-
-import Data.ExtShape
-import Data.PrimitiveArray
-
-
-
--- | Monadic arrays of primitive type.
-
-data MArr0 s sh elm = MArr0 !sh !(MVector s elm)
-
--- | Immutable arrays of primitive type.
-
-data Arr0 sh elm = Arr0 !sh !(Vector elm)
-
-
-
-type instance MutArray Arr0 = MArr0
-
--- NOTE inLb, inUb is including bound, while exUb is excluding upper bound.
--- Differentiates between largest included index, first excluded index.
-
-instance (Shape sh, ExtShape sh, VU.Unbox elm) => MPrimArrayOps MArr0 sh elm where
-  boundsM (MArr0 exUb _) = (zeroDim,exUb `subDim` unitDim)
-  fromListM inLb inUb xs = do
-    ma <- newM inLb inUb
-    let exUb = inUb `addDim` unitDim
-    let (MArr0 _ mba) = ma
-    zipWithM_ (\k x -> assert (length xs == size exUb) $ unsafeWrite mba k x) [0.. toIndex exUb inUb] xs
-    return ma
-  newM inLb inUb = let exUb = inUb `addDim` unitDim in
-    unless (inLb == zeroDim) (error "MArr0 lb/=zeroDim") >>
-    MArr0 exUb `liftM` new (size exUb)
-  newWithM inLb inUb def = do
-    let exUb = inUb `addDim` unitDim
-    ma <- newM inLb inUb
-    let (MArr0 _ mba) = ma
-    forM_ [0 .. toIndex exUb inUb] $ \k -> unsafeWrite mba k def
-    return ma
-  readM (MArr0 exUb mba) idx = assert (inShape exUb idx) $ unsafeRead mba (toIndex exUb idx)
-  writeM (MArr0 exUb mba) idx elm = assert (inShape exUb idx) $ unsafeWrite mba (toIndex exUb idx) elm
-  {-# INLINE boundsM #-}
-  {-# INLINE fromListM #-}
-  {-# INLINE newM #-}
-  {-# INLINE newWithM #-}
-  {-# INLINE readM #-}
-  {-# INLINE writeM #-}
-
-instance (Shape sh, ExtShape sh, VUM.Unbox elm) => PrimArrayOps Arr0 sh elm where
-  bounds (Arr0 exUb _) = (zeroDim,exUb `subDim` unitDim)
-  freeze (MArr0 exUb mba) = Arr0 exUb `liftM` unsafeFreeze mba
-  index (Arr0 exUb ba) idx = assert (inShape exUb idx) $ unsafeIndex ba (toIndex exUb idx)
-  {-# INLINE bounds #-}
-  {-# INLINE freeze #-}
-  {-# INLINE index #-}
-
diff --git a/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -1,4 +1,4 @@
-Copyright Christian Hoener zu Siederdissen 2010
+Copyright Christian Hoener zu Siederdissen 2010-2013
 
 All rights reserved.
 
diff --git a/PrimitiveArray.cabal b/PrimitiveArray.cabal
--- a/PrimitiveArray.cabal
+++ b/PrimitiveArray.cabal
@@ -1,10 +1,10 @@
 Name:           PrimitiveArray
-Version:        0.4.0.1
+Version:        0.5.0.0
 License:        BSD3
 License-file:   LICENSE
 Author:         Christian Hoener zu Siederdissen
 Maintainer:     choener@tbi.univie.ac.at
-Copyright:      Christian Hoener zu Siederdissen, 2010-2012
+Copyright:      Christian Hoener zu Siederdissen, 2010-2013
 Homepage:       http://www.tbi.univie.ac.at/~choener/
 Stability:      Experimental
 Category:       Data
@@ -23,24 +23,23 @@
                 zero-based/lower-bound arrays. Zero-based arrays save one
                 addition on each access if the lower bound or the array is
                 always zero.
-                .
-                We have forked two repa modules: Shape and Index.
 
 Library
   Exposed-modules:
-    Data.Array.Repa.Index
-    Data.Array.Repa.Shape
-    Data.ExtShape
+    Data.Array.Repa.ExtShape
+    Data.Array.Repa.Index.Point
+    Data.Array.Repa.Index.Subword
     Data.PrimitiveArray
+    Data.PrimitiveArray.FillTables
     Data.PrimitiveArray.Zero
-    Data.PrimitiveArray.Zero.Unboxed
---    Data.PrimitiveArray.UpperTriangular
---    Data.PrimitiveArray.UpperTriangular.Unboxed
---    Data.PrimitiveArray.FillTable
   Build-depends:
-    base      >= 4 && <5  ,
-    primitive >= 0.5.0.1  ,
-    vector    >= 0.10.0.1
+    base            >= 4 && <5  ,
+    deepseq         >= 1.3      ,
+    primitive       >= 0.5.0.1  ,
+    vector          >= 0.10.0.1 ,
+    vector-th-unbox >= 0.2      ,
+    repa            >= 3.2.3    ,
+    QuickCheck      >= 2.5
   ghc-options:
     -O2
     -funbox-strict-fields