diff --git a/Data/Array/Repa/ExtShape.hs b/Data/Array/Repa/ExtShape.hs
deleted file mode 100644
--- a/Data/Array/Repa/ExtShape.hs
+++ /dev/null
@@ -1,41 +0,0 @@
-{-# 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/Outside.hs b/Data/Array/Repa/Index/Outside.hs
deleted file mode 100644
--- a/Data/Array/Repa/Index/Outside.hs
+++ /dev/null
@@ -1,173 +0,0 @@
-{-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE TemplateHaskell #-}
-{-# LANGUAGE FlexibleInstances #-}
-{-# LANGUAGE TypeOperators #-}
-
--- | 'Outside' covers subwords for outside calculations. These types of
--- calculations requires quite "weird" index movements if you want to stay with
--- usual grammars. This remains true if grammars are transformed to Chomsky
--- normal form, only that in said form it is easier to write down the
--- recursions.
---
--- TODO basically untested!
-
-module Data.Array.Repa.Index.Outside where
-
-import           Control.Applicative
-import           Control.DeepSeq
-import           Data.Array.Repa.Index
-import           Data.Array.Repa.Shape
-import           Data.Vector.Unboxed.Deriving
-import           GHC.Base (quotInt, remInt)
-import qualified Data.Vector.Generic.Base
-import qualified Data.Vector.Generic.Mutable
-import qualified Data.Vector.Unboxed as VU
-import           Test.QuickCheck
-import           Test.QuickCheck.All
-
-import           Data.Array.Repa.ExtShape
-import           Data.Array.Repa.Index.Subword hiding (upperTri, subwordIndex, subwordFromIndex,stage)
-import qualified Data.Array.Repa.Index.Subword as SW
-
-
-
-stage = "Data.Array.Repa.Index.Outside"
-
--- | 'Outside' inverts the usual subword (i,j) system.
---
--- TODO do I need to store N ?
-
-newtype Outside = Outside (Int:.Int)
-  deriving (Eq,Ord,Show)
-
-derivingUnbox "Outside"
-  [t| Outside -> (Int,Int) |]
-  [| \ (Outside (i:.j)) -> (i,j) |]
-  [| \ (i,j) -> Outside (i:.j) |]
-
-outside :: Int -> Int -> Outside
-outside i j = Outside (i:.j)
-{-# INLINE outside #-}
-
--- | Size of an upper triangle starting at 'i' and ending at 'j'. "(0,N)" what
--- be the normal thing to use. Internally, we stell upper triangular matrices.
-
-upperTri :: Outside -> Int
-upperTri (Outside (i:.j)) = triangularNumber $ j-i
-{-# INLINE upperTri #-}
-
--- | Outside 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 :: Outside -> Outside -> Int
-subwordIndex (Outside (l:.n)) (Outside (i:.j)) = adr n (i,j) -- - adr n (l,n)
-  where
-    adr n (i,j) = n*i - triangularNumber i + j
-{-# INLINE subwordIndex #-}
-
-subwordFromIndex :: Outside -> Int -> Outside
-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 :. Outside) where
-  {-# INLINE [1] rank #-}
-  rank   (sh  :. _)
-    = rank sh + 1
-  {-# INLINE [1] zeroDim #-}
-  zeroDim = zeroDim :. Outside (0:.0)
-
-  {-# INLINE [1] unitDim #-}
-  unitDim = unitDim :. Outside (0:.1)
-
-  {-# INLINE [1] intersectDim #-}
-  intersectDim (sh1 :. Outside (i:.j)) (sh2 :. Outside (k:.l))
-    = (intersectDim sh1 sh2 :. Outside (max i k :. min j l))
-
-  {-# INLINE [1] addDim #-}
-  addDim (sh1 :. Outside (i:.j)) (sh2 :. Outside (k:.l))
-    = addDim sh1 sh2 :. Outside (i+k:.j+l)
-
-  {-# INLINE [1] size #-}
-  size  (sh1 :. sw) = size sh1 * upperTri sw
-
-  {-# INLINE [1] sizeIsValid #-}
-  sizeIsValid (sh1 :. Outside (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 :. Outside (_:._)) (sh1 :. Outside (l:.n)) (sh2 :. Outside (i:.j))
-    = i<=j && l<=i && j<n && (inShapeRange zs sh1 sh2)
-
-  {-# NOINLINE listOfShape #-}
-  listOfShape (sh :. Outside (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 :. Outside (i:.j)
-
-  {-# INLINE deepSeq #-}
-  deepSeq (sh :. n) x = deepSeq sh (n `seq` x)
-
--- |
-
-instance ExtShape sh => ExtShape (sh:.Outside) where
-  subDim (sh1:.Outside (i:.j)) (sh2:.Outside (k:.l)) = subDim sh1 sh2 :. Outside (i-k:.j-l)
-  {-# INLINE subDim #-}
-  rangeList (sh1:.Outside (i:.j)) (sh2:.Outside (k:.l)) = error "not implemented" -- [sh:.Outside (m,n) | sh <- rangeList sh1 sh2, m <- [i .. [i+k], n <- [ n <- [n1 .. (n1+n2) ] ]
-  {-# INLINE rangeList #-}
-
--- |
-
-instance NFData Outside where
-  rnf (Outside (i:.j)) = i `seq` rnf j
-
--- |
-
-instance Arbitrary Outside where
-  arbitrary = do
-    a <- choose (0,100)
-    b <- choose (0,100)
-    return $ Outside (min a b :. max a b)
-  shrink (Outside (i:.j))
-    | i<j       = [Outside (i:.j-1)]
-    | otherwise = []
-
-instance Arbitrary z => Arbitrary (z:.Outside) where
-  arbitrary = (:.) <$> arbitrary <*> arbitrary
-  shrink (z:.s) = (:.) <$> shrink z <*> shrink s
-
diff --git a/Data/Array/Repa/Index/Point.hs b/Data/Array/Repa/Index/Point.hs
deleted file mode 100644
--- a/Data/Array/Repa/Index/Point.hs
+++ /dev/null
@@ -1,115 +0,0 @@
-{-# 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/Points.hs b/Data/Array/Repa/Index/Points.hs
deleted file mode 100644
--- a/Data/Array/Repa/Index/Points.hs
+++ /dev/null
@@ -1,239 +0,0 @@
-{-# LANGUAGE FlexibleInstances #-}
-{-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE TemplateHaskell #-}
-{-# LANGUAGE TypeOperators #-}
-
--- | Index structures for left- and right-linear grammars. Do not use this
--- index for general linear- or context-free grammars.
---
--- Internally, both 'PointL' and 'PointR' work a lot like 'Subword's, but in
--- non-terminals we only store the left- or right part.
-
-module Data.Array.Repa.Index.Points where
-
-import           Control.Applicative
-import           Control.DeepSeq
-import           Data.Array.Repa.Index
-import           Data.Array.Repa.Shape
-import           Data.Vector.Unboxed.Deriving
-import           GHC.Base (quotInt, remInt)
-import qualified Data.Vector.Generic.Base
-import qualified Data.Vector.Generic.Mutable
-import qualified Data.Vector.Unboxed as VU
-import           Test.QuickCheck
-import           Test.QuickCheck.All
-
-import           Data.Array.Repa.ExtShape
-import           Data.Array.Repa.Index.Subword
-
-
-
--- | A point in left-linear grammars. In @(i:.j)@, @j@ is the non-terminal
--- storage point, @i==0@ always for the non-terminal, while @i>=0@ for
--- terminals, which are on the right of the non-terminal. (This is why
--- left-linear grammars are called left-linear: they recurse on the left).
---
--- PS: all this left/right talk deals with the RHS of a production rule, the
--- LHS is always a non-terminal ;-)
-
-newtype PointL = PointL (Int:.Int)
-  deriving (Eq,Read,Show)
-
-pointL :: Int -> Int -> PointL
-pointL i j = PointL (i:.j)
-{-# INLINE pointL #-}
-
--- | A point in right-linear grammars.
-
-newtype PointR = PointR (Int:.Int)
-  deriving (Eq,Read,Show)
-
-pointR :: Int -> Int -> PointR
-pointR i j = PointR (i:.j)
-{-# INLINE pointR #-}
-
-
-
--- * Instances: PointL
-
-derivingUnbox "PointL"
-  [t| PointL -> (Int,Int) |]
-  [| \ (PointL (i:.j)) -> (i,j) |]
-  [| \ (i,j) -> PointL (i:.j) |]
-
-instance Shape sh => Shape (sh :. PointL) where
-  {-# INLINE [1] rank #-}
-  rank   (sh  :. _)
-    = rank sh + 1
-
-  {-# INLINE [1] zeroDim #-}
-  zeroDim = zeroDim :. PointL (0:.0)
-
-  {-# INLINE [1] unitDim #-}
-  unitDim = unitDim :. PointL (0:.1)
-
-  {-# INLINE [1] intersectDim #-}
-  intersectDim (sh1 :. PointL (i:.j)) (sh2 :. PointL (k:.l))
-    = (intersectDim sh1 sh2 :. PointL (max i k :. min j l))
-
-  {-# INLINE [1] addDim #-}
-  addDim (sh1 :. PointL (i:.j)) (sh2 :. PointL (k:.l))
-    = addDim sh1 sh2 :. PointL (i+k:.j+l)
-
-  -- NOTE size is calculated NOT as upper-triangular, but linear!
-  {-# INLINE [1] size #-}
-  size  (sh1 :. PointL (i:.j)) = size sh1 * (j-i)
-
-  {-# INLINE [1] sizeIsValid #-}
-  sizeIsValid (sh1 :. PointL (i:.j))
-    | size sh1 > 0
-    = i>=0 && i<=j && j <= maxBound `div` size sh1
-    | otherwise
-    = False
-
-  -- NOTE only the @j@ coordinate is used for indexing NTs, @i@ is just for
-  -- convenience. @l@ however restricts the NT to some value @>0@ if desired.
-  {-# INLINE [1] toIndex #-}
-  toIndex (sh1 :. PointL(l:.r)) (sh1' :. PointL(i:.j))
-    = toIndex sh1 sh1' * (r-l) + (j-l)
-
-  {-# INLINE [1] fromIndex #-}
-  fromIndex (ds :. d) n  = undefined -- fromIndex ds (n `quotInt` d) :. r
-    where
-      r = undefined
-
-  -- | TODO fix for lower bounds check!
-  {-# INLINE [1] inShapeRange #-}
-  inShapeRange (zs :. PointL (_:._)) (sh1 :. PointL (l:.n)) (sh2 :. PointL (i:.j))
-    = i<=j && l<=i && j<n && (inShapeRange zs sh1 sh2)
-
-  {-# NOINLINE listOfShape #-}
-  listOfShape (sh :. PointL (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 :. PointL (i:.j)
-
-  {-# INLINE deepSeq #-}
-  deepSeq (sh :. n) x = deepSeq sh (n `seq` x)
-
-instance ExtShape sh => ExtShape (sh:.PointL) where
-  {-# INLINE [1] subDim #-}
-  subDim (sh1:.PointL (i:.j)) (sh2:.PointL (k:.l)) = subDim sh1 sh2 :. PointL (i-k:.j-l)
-  {-# INLINE [1] rangeList #-}
-  rangeList _ _ = error "PointL:rangeList not implemented"
-
-instance NFData PointL where
-  rnf (PointL (i:.j)) = i `seq` rnf j
-  {-# INLINE rnf #-}
-
--- TODO maybe vary the left border, too? Since this invalidates that @i==0@ in
--- @PointL (i:.j)@, we would need to make sure that the memoizers for NTs get
--- notified ...
-
-instance Arbitrary PointL where
-  arbitrary = do
-    b <- choose (0,100)
-    return $ pointL 0 b
-  shrink (PointL (i:.j))
-    | i<j = [pointL i $ j-1]
-    | otherwise = []
-
-instance Arbitrary z => Arbitrary (z:.PointL) where
-  arbitrary = (:.) <$> arbitrary <*> arbitrary
-  shrink (z:.s) = (:.) <$> shrink z <*> shrink s
-
-
-
--- * Instances: PointR
-
-derivingUnbox "PointR"
-  [t| PointR -> (Int,Int) |]
-  [| \ (PointR (i:.j)) -> (i,j) |]
-  [| \ (i,j) -> PointR (i:.j) |]
-
-instance Shape sh => Shape (sh :. PointR) where
-  {-# INLINE [1] rank #-}
-  rank   (sh  :. _)
-    = rank sh + 1
-
-  {-# INLINE [1] zeroDim #-}
-  zeroDim = zeroDim :. PointR (0:.0)
-
-  {-# INLINE [1] unitDim #-}
-  unitDim = unitDim :. PointR (0:.1)
-
-  {-# INLINE [1] intersectDim #-}
-  intersectDim (sh1 :. PointR (i:.j)) (sh2 :. PointR (k:.l))
-    = (intersectDim sh1 sh2 :. PointR (max i k :. min j l))
-
-  {-# INLINE [1] addDim #-}
-  addDim (sh1 :. PointR (i:.j)) (sh2 :. PointR (k:.l))
-    = addDim sh1 sh2 :. PointR (i+k:.j+l)
-
-  -- NOTE size is calculated NOT as upper-triangular, but linear!
-  {-# INLINE [1] size #-}
-  size  (sh1 :. PointR (i:.j)) = size sh1 * (j-i)
-
-  {-# INLINE [1] sizeIsValid #-}
-  sizeIsValid (sh1 :. PointR (i:.j))
-    | size sh1 > 0
-    = i>=0 && i<=j && j <= maxBound `div` size sh1
-    | otherwise
-    = False
-
-  -- NOTE only the @i@ coordinate is used for indexing NTs, @j@ is just for
-  -- convenience. @l@ however restricts the NT to some value @>0@ if desired.
-  {-# INLINE [1] toIndex #-}
-  toIndex (sh1 :. PointR(l:.r)) (sh1' :. PointR(i:.j))
-    = toIndex sh1 sh1' * (r-l) + i
-
-  {-# INLINE [1] fromIndex #-}
-  fromIndex (ds :. d) n  = undefined -- fromIndex ds (n `quotInt` d) :. r
-    where
-      r = undefined
-
-  -- | TODO fix for lower bounds check!
-  {-# INLINE [1] inShapeRange #-}
-  inShapeRange (zs :. PointR (_:._)) (sh1 :. PointR (l:.n)) (sh2 :. PointR (i:.j))
-    = i<=j && l<=i && j<n && (inShapeRange zs sh1 sh2)
-
-  {-# NOINLINE listOfShape #-}
-  listOfShape (sh :. PointR (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 :. PointR (i:.j)
-
-  {-# INLINE deepSeq #-}
-  deepSeq (sh :. n) x = deepSeq sh (n `seq` x)
-
-instance ExtShape sh => ExtShape (sh:.PointR) where
-  {-# INLINE [1] subDim #-}
-  subDim (sh1:.PointR (i:.j)) (sh2:.PointR (k:.l)) = subDim sh1 sh2 :. PointR (i-k:.j-l)
-  {-# INLINE [1] rangeList #-}
-  rangeList _ _ = error "PointR:rangeList not implemented"
-
-instance NFData PointR where
-  rnf (PointR (i:.j)) = i `seq` rnf j
-  {-# INLINE rnf #-}
-
-instance Arbitrary PointR where
-  arbitrary = do
-    b <- choose (0,100)
-    return $ pointR b 100
-  shrink (PointR (i:.j))
-    | i<j = [pointR (i+1) $ j]
-    | otherwise = []
-
-instance Arbitrary z => Arbitrary (z:.PointR) 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
deleted file mode 100644
--- a/Data/Array/Repa/Index/Subword.hs
+++ /dev/null
@@ -1,188 +0,0 @@
-{-# 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           Data.Vector.Unboxed.Deriving
-import           GHC.Base (quotInt, remInt)
-import qualified Data.Vector.Generic.Base
-import qualified Data.Vector.Generic.Mutable
-import qualified Data.Vector.Unboxed as VU
-import           Test.QuickCheck
-import           Test.QuickCheck.All
-
-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/PrimitiveArray.hs b/Data/PrimitiveArray.hs
--- a/Data/PrimitiveArray.hs
+++ b/Data/PrimitiveArray.hs
@@ -1,166 +1,13 @@
-{-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE ScopedTypeVariables #-}
-{-# LANGUAGE FlexibleContexts #-}
-{-# LANGUAGE FlexibleInstances #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
 
--- | Vastly extended primitive arrays. Some basic ideas are now modeled after
--- the vector package, especially the monadic mutable / pure immutable array
--- 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
-import Data.Primitive.Types
-import Prelude as P
-import System.IO.Unsafe
-
-import Data.Array.Repa.ExtShape
-
-
-
--- | Mutable version of an array.
-
-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 :: 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 (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 (MutArr m (arr sh elm))
-
-  -- | Creates a new array with all elements being equal to 'elm'.
-
-  newWithM :: PrimMonad m => sh -> sh -> elm -> m (MutArr m (arr sh elm))
-
-  -- | Reads a single element in the array.
-
-  readM :: PrimMonad m => MutArr m (arr sh elm) -> sh -> m elm
-
-  -- | Writes a single element in the array.
-
-  writeM :: PrimMonad m => MutArr m (arr sh elm) -> sh -> elm -> m ()
-
-
-
--- | The core set of functions on immutable arrays.
-
-class (Shape sh, ExtShape sh) => PrimArrayOps arr sh elm where
-
-  -- | Returns the bounds of an immutable array, again inclusive bounds: @ [lb..ub] @.
-
-  bounds :: arr sh elm -> (sh,sh)
-
-  -- | Freezes a mutable array an returns its immutable version. This operation
-  -- is /O(1)/ and both arrays share the same memory. Do not use the mutable
-  -- array afterwards.
-
-  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.
-
-(!) :: PrimArrayOps arr sh elm => arr sh elm -> sh -> elm
-(!) arr idx = assert (inBounds arr idx) $ index arr idx
-{-# INLINE (!) #-}
-
--- | Returns true if the index is valid for the array.
-
-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 #-}
-
--- | Given two arrays with the same dimensionality, their respective starting
--- index, and how many steps to go in each dimension (in terms of a dimension
--- again), determine if the multidimensional slices have the same value at
--- all positions
---
--- TODO specialize for DIM1 (and maybe higher dim's) to use memcmp
-
-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 = 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 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
-  return ma
-{-# INLINE fromAssocsM #-}
-
--- | Return all associations from an array.
-
-assocs :: PrimArrayOps arr sh elm => arr sh elm -> [(sh,elm)]
-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, 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, 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 #-}
-
--- | Determines if an index is valid for a given immutable array.
-
-inBounds :: PrimArrayOps arr sh elm => arr sh elm -> sh -> Bool
-inBounds arr idx = let (lb,ub) = bounds arr in inShapeRange lb (ub `addDim` unitDim) idx
-{-# INLINE inBounds #-}
-
--- | Returns all elements of an immutable array as a list.
+module Data.PrimitiveArray 
+  ( module Data.PrimitiveArray.Class
+  , module Data.PrimitiveArray.Dense
+  , module Data.PrimitiveArray.FillTables
+  , module Data.PrimitiveArray.Index
+  ) where
 
-toList :: PrimArrayOps arr sh elm =>  arr sh elm -> [elm]
-toList arr = let (lb,ub) = bounds arr in P.map ((!) arr) $ rangeList lb $ ub `subDim` lb
-{-# INLINE toList #-}
+import Data.PrimitiveArray.Class
+import Data.PrimitiveArray.Dense
+import Data.PrimitiveArray.FillTables
+import Data.PrimitiveArray.Index
 
diff --git a/Data/PrimitiveArray/Class.hs b/Data/PrimitiveArray/Class.hs
new file mode 100644
--- /dev/null
+++ b/Data/PrimitiveArray/Class.hs
@@ -0,0 +1,188 @@
+
+-- | Vastly extended primitive arrays. Some basic ideas are now modeled after
+-- the vector package, especially the monadic mutable / pure immutable array
+-- system.
+--
+-- NOTE all operations in MPrimArrayOps and PrimArrayOps are highly unsafe. No
+-- bounds-checking is performed at all.
+
+module Data.PrimitiveArray.Class where
+
+import           Control.Applicative (Applicative, pure, (<$>), (<*>))
+import           Control.Exception (assert)
+import           Control.Monad (forM_)
+import           Control.Monad.Primitive (PrimMonad)
+import           Control.Monad.ST (runST)
+import           Prelude as P
+import qualified Data.Vector.Fusion.Stream as S
+
+import Data.PrimitiveArray.Index
+
+
+
+-- | Mutable version of an array.
+
+data family MutArr (m :: * -> *) (arr :: *) :: *
+
+
+-- | The core set of operations for monadic arrays.
+
+class (Index sh) => MPrimArrayOps arr sh elm where
+
+  -- | Return the bounds of the array. All bounds are inclusive, as in
+  -- @[lb..ub]@
+
+  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 (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 (MutArr m (arr sh elm))
+
+  -- | Creates a new array with all elements being equal to 'elm'.
+
+  newWithM :: PrimMonad m => sh -> sh -> elm -> m (MutArr m (arr sh elm))
+
+  -- | Reads a single element in the array.
+
+  readM :: PrimMonad m => MutArr m (arr sh elm) -> sh -> m elm
+
+  -- | Writes a single element in the array.
+
+  writeM :: PrimMonad m => MutArr m (arr sh elm) -> sh -> elm -> m ()
+
+
+
+-- | The core set of functions on immutable arrays.
+
+class (Index sh) => PrimArrayOps arr sh elm where
+
+  -- | Returns the bounds of an immutable array, again inclusive bounds: @ [lb..ub] @.
+
+  bounds :: arr sh elm -> (sh,sh)
+
+  -- | Freezes a mutable array an returns its immutable version. This operation
+  -- is /O(1)/ and both arrays share the same memory. Do not use the mutable
+  -- array afterwards.
+
+  unsafeFreeze :: PrimMonad m => MutArr m (arr sh elm) -> m (arr sh elm)
+
+  -- | Thaw an immutable array into a mutable one. Both versions share
+  -- memory.
+
+  unsafeThaw :: PrimMonad m => arr sh elm -> m (MutArr m (arr sh elm))
+
+  -- | Extract a single element from the array. Generally unsafe as not
+  -- bounds-checking is performed.
+
+  unsafeIndex :: arr sh elm -> sh -> elm
+
+  -- | Savely transform the shape space of a table.
+
+  transformShape :: (Index sh') => (sh -> sh') -> arr sh elm -> arr sh' elm
+
+class (Index 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.
+
+(!) :: PrimArrayOps arr sh elm => arr sh elm -> sh -> elm
+(!) arr idx = assert (uncurry inBounds (bounds arr) idx) $ unsafeIndex arr idx
+{-# INLINE (!) #-}
+
+-- | Returns true if the index is valid for the array.
+
+inBoundsM :: (Monad m, MPrimArrayOps arr sh elm) => MutArr m (arr sh elm) -> sh -> Bool
+inBoundsM marr idx = let (lb,ub) = boundsM marr in inBounds lb ub idx
+{-# INLINE inBoundsM #-}
+
+-- -- | Given two arrays with the same dimensionality, their respective starting
+-- -- index, and how many steps to go in each dimension (in terms of a dimension
+-- -- again), determine if the multidimensional slices have the same value at
+-- -- all positions
+-- --
+-- -- TODO specialize for DIM1 (and maybe higher dim's) to use memcmp
+-- 
+-- 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 = P.map (unsafeIndex arr1) $ rangeList k1 xtnd
+--   ys = P.map (unsafeIndex 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 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
+  return ma
+{-# INLINE fromAssocsM #-}
+
+-- | Return all associations from an array.
+
+assocs :: (IndexStream sh, PrimArrayOps arr sh elm) => arr sh elm -> [(sh,elm)]
+assocs arr = P.map (\k -> (k,unsafeIndex arr k)) . S.toList $ streamUp lb ub 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, MPrimArrayOps arr sh elm) => sh -> sh -> [elm] -> arr sh elm
+fromList lb ub xs = runST $ fromListM lb ub xs >>= unsafeFreeze
+{-# INLINE fromList #-}
+
+-- | Creates an immutable array from lower and upper bounds, a default element,
+-- and a list of associations.
+
+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 >>= unsafeFreeze
+{-# INLINE fromAssocs #-}
+
+-- -- | Determines if an index is valid for a given immutable array.
+-- 
+-- inBounds :: PrimArrayOps arr sh elm => arr sh elm -> sh -> Bool
+-- inBounds arr idx = let (lb,ub) = bounds arr in inShapeRange lb (ub `addDim` unitDim) idx
+-- {-# INLINE inBounds #-}
+
+-- | Returns all elements of an immutable array as a list.
+
+toList :: (IndexStream sh, PrimArrayOps arr sh elm) => arr sh elm -> [elm]
+toList arr = let (lb,ub) = bounds arr in P.map ((!) arr) . S.toList $ streamUp lb ub
+{-# INLINE toList #-}
+
+
+
+-- * Freeze an inductive stack of tables with a 'Z' at the bottom.
+
+-- | 'freezeTables' freezes a stack of tables.
+
+class FreezeTables m t where
+    type Frozen t :: *
+    freezeTables :: t -> m (Frozen t)
+
+instance Applicative m => FreezeTables m Z where
+    type Frozen Z = Z
+    freezeTables Z = pure Z
+    {-# INLINE freezeTables #-}
+
+instance (Functor m, Applicative m, Monad m, PrimMonad m, FreezeTables m ts, PrimArrayOps arr sh elm) => FreezeTables m (ts:.MutArr m (arr sh elm)) where
+    type Frozen (ts:.MutArr m (arr sh elm)) = Frozen ts :. arr sh elm
+    freezeTables (ts:.t) = (:.) <$> freezeTables ts <*> unsafeFreeze t
+    {-# INLINE freezeTables #-}
+
diff --git a/Data/PrimitiveArray/Dense.hs b/Data/PrimitiveArray/Dense.hs
new file mode 100644
--- /dev/null
+++ b/Data/PrimitiveArray/Dense.hs
@@ -0,0 +1,178 @@
+
+-- | Dense primitive arrays where the lower index is zero (or the
+-- equivalent of zero for newtypes and enumerations).
+--
+-- Actual @write@s to data structures use a more safe @write@ instead of
+-- the unsafe @unsafeWrite@. Writes also tend to occur much less in DP
+-- algorithms (say, N^2 writes for an N^3 time algorithm -- mostly reads
+-- are being executed).
+--
+-- TODO consider if we want to force the lower index to be zero, or allow
+-- non-zero lower indices. Will have to be considered together with the
+-- @Index.Class@ module!
+
+module Data.PrimitiveArray.Dense where
+
+import           Control.DeepSeq
+import           Control.Exception (assert)
+import           Control.Monad (liftM, forM_, zipWithM_)
+import           Control.Monad.Primitive (PrimState)
+import           Data.Aeson (ToJSON,FromJSON)
+import           Data.Binary (Binary)
+import           Data.Serialize (Serialize)
+import           Data.Vector.Binary
+import           Data.Vector.Cereal
+import           Data.Vector.Generic.Mutable as GM hiding (length)
+import           Data.Vector.Unboxed.Mutable (Unbox)
+import           GHC.Generics (Generic)
+import qualified Data.Vector as V hiding (forM_, length, zipWithM_)
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed as VU hiding (forM_, length, zipWithM_)
+
+
+import           Data.PrimitiveArray.Class
+import           Data.PrimitiveArray.Index
+
+
+
+-- * Unboxed, multidimensional arrays.
+
+data Unboxed sh e = Unboxed !sh !sh !(VU.Vector e)
+  deriving (Read,Show,Eq,Generic)
+
+instance (Binary    sh, Binary    e, Unbox e) => Binary    (Unboxed sh e)
+instance (Serialize sh, Serialize e, Unbox e) => Serialize (Unboxed sh e)
+instance (ToJSON    sh, ToJSON    e, Unbox e) => ToJSON    (Unboxed sh e)
+instance (FromJSON  sh, FromJSON  e, Unbox e) => FromJSON  (Unboxed sh e)
+
+instance NFData (Unboxed sh e) where
+  rnf !_ = ()
+
+data instance MutArr m (Unboxed sh e) = MUnboxed !sh !sh !(VU.MVector (PrimState m) e)
+
+instance NFData (MutArr m (Unboxed sh e)) where
+  rnf !_ = ()
+
+instance (Index sh, Unbox elm) => MPrimArrayOps Unboxed sh elm where
+  boundsM (MUnboxed l h _) = (l,h)
+  fromListM l h xs = do
+    ma <- newM l h
+    let (MUnboxed _ _ mba) = ma
+    zipWithM_ (\k x -> assert (length xs == size l h) $ unsafeWrite mba k x) [0.. size l h -1] xs
+    return ma
+  newM l h = MUnboxed l h `liftM` new (size l h)
+  newWithM l h def = do
+    ma <- newM l h
+    let (MUnboxed _ _ mba) = ma
+    forM_ [0 .. size l h -1] $ \k -> unsafeWrite mba k def
+    return ma
+  readM  (MUnboxed l h mba) idx     = assert (inBounds l h idx) $ unsafeRead  mba (linearIndex l h idx)
+  writeM (MUnboxed l h mba) idx elm = write mba (linearIndex l h idx) elm
+  {-# INLINE boundsM #-}
+  {-# INLINE fromListM #-}
+  {-# INLINE newM #-}
+  {-# INLINE newWithM #-}
+  {-# INLINE readM #-}
+  {-# INLINE writeM #-}
+
+instance (Index sh, Unbox elm) => PrimArrayOps Unboxed sh elm where
+  bounds (Unboxed l h _) = (l,h)
+  unsafeFreeze (MUnboxed l h mba) = Unboxed l h `liftM` G.unsafeFreeze mba
+  unsafeThaw   (Unboxed  l h ba) = MUnboxed l h `liftM` G.unsafeThaw ba
+  unsafeIndex  (Unboxed  l h ba) idx = {- assert (inShape exUb idx) $ -} G.unsafeIndex ba (linearIndex l h idx)
+  transformShape tr (Unboxed l h ba) = Unboxed (tr l) (tr h) ba
+  {-# INLINE bounds #-}
+  {-# INLINE unsafeFreeze #-}
+  {-# INLINE unsafeThaw #-}
+  {-# INLINE unsafeIndex #-}
+  {-# INLINE transformShape #-}
+
+instance (Index sh, Unbox e, Unbox e') => PrimArrayMap Unboxed sh e e' where
+  map f (Unboxed l h xs) = Unboxed l h (VU.map f xs)
+  {-# INLINE map #-}
+
+
+
+-- * Boxed, multidimensional arrays.
+
+data Boxed sh e = Boxed !sh !sh !(V.Vector e)
+  deriving (Read,Show,Eq,Generic)
+
+instance (Binary    sh, Binary    e)  => Binary    (Boxed sh e)
+instance (Serialize sh, Serialize e)  => Serialize (Boxed sh e)
+instance (ToJSON    sh, ToJSON    e)  => ToJSON    (Boxed sh e)
+instance (FromJSON  sh, FromJSON  e)  => FromJSON  (Boxed sh e)
+
+instance NFData (Boxed sh e) where
+  rnf !_ = ()
+
+data instance MutArr m (Boxed sh e) = MBoxed !sh !sh !(V.MVector (PrimState m) e)
+
+instance NFData (MutArr m (Boxed sh e)) where
+  rnf !_ = ()
+
+instance (Index sh) => MPrimArrayOps Boxed sh elm where
+  boundsM (MBoxed l h _) = (l,h)
+  fromListM l h xs = do
+    ma <- newM l h
+    let (MBoxed _ _ mba) = ma
+    zipWithM_ (\k x -> assert (length xs == size l h) $ unsafeWrite mba k x) [0 .. size l h - 1] xs
+    return ma
+  newM l h =
+    MBoxed l h `liftM` new (size l h)
+  newWithM l h def = do
+    ma <- newM l h
+    let (MBoxed _ _ mba) = ma
+    forM_ [0 .. size l h -1] $ \k -> unsafeWrite mba k def
+    return ma
+  readM  (MBoxed l h mba) idx     = assert (inBounds l h idx) $ GM.unsafeRead mba (linearIndex l h idx)
+  writeM (MBoxed l h mba) idx elm = assert (inBounds l h idx) $ GM.write mba (linearIndex l h idx) elm
+  {-# INLINE boundsM #-}
+  {-# INLINE fromListM #-}
+  {-# INLINE newM #-}
+  {-# INLINE newWithM #-}
+  {-# INLINE readM #-}
+  {-# INLINE writeM #-}
+
+instance (Index sh, Unbox elm) => PrimArrayOps Boxed sh elm where
+  bounds (Boxed l h _) = (l,h)
+  unsafeFreeze (MBoxed l h mba) = Boxed l h `liftM` G.unsafeFreeze mba
+  unsafeThaw   (Boxed l h ba) = MBoxed l h `liftM` G.unsafeThaw ba
+  unsafeIndex (Boxed l h ba) idx = {- assert (inShape exUb idx) $ -} G.unsafeIndex ba (linearIndex l h idx)
+  transformShape tr (Boxed l h ba) = Boxed (tr l) (tr h) ba
+  {-# INLINE bounds #-}
+  {-# INLINE unsafeFreeze #-}
+  {-# INLINE unsafeThaw #-}
+  {-# INLINE unsafeIndex #-}
+  {-# INLINE transformShape #-}
+
+instance (Index sh) => PrimArrayMap Boxed sh e e' where
+  map f (Boxed l h xs) = Boxed l h (V.map f xs)
+  {-# INLINE map #-}
+
+
+
+{-
+ -
+ - 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/FillTables.hs b/Data/PrimitiveArray/FillTables.hs
--- a/Data/PrimitiveArray/FillTables.hs
+++ b/Data/PrimitiveArray/FillTables.hs
@@ -1,103 +1,64 @@
-{-# 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.
+-- | Operations to fill primitive arrays. Arrays are combined just like
+-- indices using 'Z' and '(:.)'. This allows filling an unlimited number of
+-- tables. 'ExtShape' provides the 'rangeStream' function with generates
+-- a stream of indices in (generally) the right order.
 
 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 Control.Monad (when)
+import Data.Vector.Fusion.Stream as S
+import Data.Vector.Fusion.Stream.Monadic as M
+import Data.Vector.Fusion.Stream.Size
 
-import Data.PrimitiveArray
-import Data.Array.Repa.Index.Subword
+import Data.PrimitiveArray.Class
+import Data.PrimitiveArray.Index
 
 
 
--- * Driver classes for table filling system.
+-- * High-level table filling system.
 
--- Upper triangular table filling. Right now, only a serial option 'upperTriS'
--- is available.
+-- | Run the forward phase of algorithms. Is *really* unsafe for now if
+-- tables have different sizes, as in its broken.
 --
--- TODO Using Repa, 'upperTriP' will soon become available.
-
-class UpperTriS m stack where
-  upperTriS :: stack -> m ()
+-- TODO Need to run min/max on the bounds for all tables, not just the last
+-- table. Otherwise we don't really need the distinction between save and
+-- unsafe. This will have to be in @runFillTables@.
 
--- | Defines how a single index in a stack of arrays + evaluation functions is
--- handled. The instances *should* work for any index @ix@.
+unsafeRunFillTables
+  :: ( Index sh, IndexStream sh
+     , WriteCell m (tail :. (MutArr m (arr sh elm), t)) sh
+     , MPrimArrayOps arr sh elm
+     , Monad m
+     , PrimMonad m
+     )
+  => (tail :. (MutArr m (arr sh elm), t)) -> m ()
 
-class Stack m sh xs where
-  writeStack :: xs -> sh -> m ()
+unsafeRunFillTables (ts:.(t,f)) = M.mapM_ (unsafeWriteCell (ts:.(t,f))) $ streamUp from to where -- generateIndices from to where
+  (from,to) = boundsM t -- TODO min/max over all tables [for the safe version, the unsafe version *always* assumes equal-size tables; we still should check this during runtime]
+{-# INLINE unsafeRunFillTables #-}
 
 
 
--- * Instances
-
--- ** 1-tape grammars with 'Subword' indices.
-
--- TODO Insert check that all extends are the same!
-
-instance
-  ( Monad m
-  , MPrimArrayOps arr (Z:.Subword) e
-  , Stack m Subword (xs :. SubwordNonTerminal m arr e)
-  ) => UpperTriS m  (xs :. SubwordNonTerminal m arr e) where
-  upperTriS xs@(_:.(x,f)) = do
-    -- TODO missing extends check
-    let (Z:.Subword (l:._),Z:.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 #-}
-
-instance
-  ( PrimMonad m
-  , Stack m Subword xs
-  , MPrimArrayOps arr (Z:.Subword) e
-  ) => Stack m Subword (xs :. SubwordNonTerminal m arr e) where
-  writeStack (xs:.(x,f)) i = writeStack xs i >> f i >>= writeM x (Z:.i)
-  {-# INLINE writeStack #-}
-
--- ** Multi-tape indices.
-
-instance (Monad m) => Stack m sh Z where
-  writeStack _ _ = return ()
-  {-# INLINE writeStack #-}
-
-instance
-  ( PrimMonad m
-  , Stack m ix xs
-  , MPrimArrayOps arr ix e
-  ) => Stack m ix (xs :. GeneralNonTerminal m arr ix e) where
-  writeStack (xs:.(x,f)) i = writeStack xs i >> f i >>= writeM x i
-  {-# INLINE writeStack #-}
+-- * Write to individuel cells.
 
+-- | 'WriteCell' provides methods to fill all cells with a specific index
+-- @sh@ in a stack of non-terminal tables @c@.
 
+class (Monad m) => WriteCell m c sh where
+    unsafeWriteCell :: c -> sh -> m ()
+    writeCell       :: c -> sh -> m ()
 
--- Wrap non-terminal symbol type, corresponding rule type.
+instance (Monad m) => WriteCell m Z sh where
+    unsafeWriteCell _ _ = return ()
+    writeCell _ _ = return ()
+    {-# INLINE unsafeWriteCell #-}
+    {-# INLINE writeCell #-}
 
-type SubwordNonTerminal m arr e = (MutArr m (arr (Z:.Subword) e), Subword -> m e)
-type GeneralNonTerminal m arr ix e = (MutArr m (arr ix e), ix -> m e)
+instance (WriteCell m cs sh, Monad m, MPrimArrayOps arr sh a, PrimMonad m) => WriteCell m (cs:.(MutArr m (arr sh a), sh -> m a)) sh where
+    unsafeWriteCell (cs:.(t,f)) sh = unsafeWriteCell cs sh >> (f sh >>= writeM t sh)
+    writeCell (cs:.(t,f)) sh = writeCell cs sh >> (when (inBoundsM t sh) (f sh >>= writeM t sh))
+    {-# INLINE unsafeWriteCell #-}
+    {-# INLINE writeCell #-}
 
diff --git a/Data/PrimitiveArray/Index.hs b/Data/PrimitiveArray/Index.hs
new file mode 100644
--- /dev/null
+++ b/Data/PrimitiveArray/Index.hs
@@ -0,0 +1,19 @@
+
+module Data.PrimitiveArray.Index
+  ( module Data.PrimitiveArray.Index.Class
+  , module Data.PrimitiveArray.Index.Complement
+  , module Data.PrimitiveArray.Index.Int
+  , module Data.PrimitiveArray.Index.Outside
+  , module Data.PrimitiveArray.Index.Point
+  , module Data.PrimitiveArray.Index.Set
+  , module Data.PrimitiveArray.Index.Subword
+  ) where
+
+import Data.PrimitiveArray.Index.Class
+import Data.PrimitiveArray.Index.Complement
+import Data.PrimitiveArray.Index.Int
+import Data.PrimitiveArray.Index.Outside
+import Data.PrimitiveArray.Index.Point
+import Data.PrimitiveArray.Index.Set
+import Data.PrimitiveArray.Index.Subword
+
diff --git a/Data/PrimitiveArray/Index/Class.hs b/Data/PrimitiveArray/Index/Class.hs
new file mode 100644
--- /dev/null
+++ b/Data/PrimitiveArray/Index/Class.hs
@@ -0,0 +1,204 @@
+
+module Data.PrimitiveArray.Index.Class where
+
+import           Control.Applicative
+import           Control.DeepSeq (NFData(..))
+import           Control.Monad (liftM2)
+import           Data.Aeson
+import           Data.Binary
+import           Data.Serialize
+import           Data.Vector.Fusion.Stream.Monadic (Stream)
+import           Data.Vector.Unboxed.Deriving
+import           Data.Vector.Unboxed (Unbox(..))
+import           GHC.Generics
+import qualified Data.Vector.Fusion.Stream.Monadic as SM
+import           Test.QuickCheck
+
+
+
+infixl 3 :.
+
+-- | Strict pairs -- as in @repa@.
+
+data a :. b = !a :. !b
+  deriving (Eq,Ord,Show,Generic)
+
+derivingUnbox "StrictPair"
+  [t| forall a b . (Unbox a, Unbox b) => (a:.b) -> (a,b) |]
+  [| \(a:.b) -> (a, b) |]
+  [| \(a,b)  -> (a:.b) |]
+
+instance (Binary    a, Binary    b) => Binary    (a:.b)
+instance (Serialize a, Serialize b) => Serialize (a:.b)
+instance (ToJSON    a, ToJSON    b) => ToJSON    (a:.b)
+instance (FromJSON  a, FromJSON  b) => FromJSON  (a:.b)
+
+deriving instance (Read a, Read b) => Read (a:.b)
+
+instance (NFData a, NFData b) => NFData (a:.b) where
+  rnf (a:.b) = rnf a `seq` rnf b
+  {-# Inline rnf #-}
+
+instance (Arbitrary a, Arbitrary b) => Arbitrary (a :. b) where
+  arbitrary     = liftM2 (:.) arbitrary arbitrary
+  shrink (a:.b) = [ (a':.b) | a' <- shrink a ] ++ [ (a:.b') | b' <- shrink b ]
+
+
+
+infixl 3 :>
+
+-- | A different version of strict pairs. Makes for simpler type inference in
+-- multi-tape grammars. We use @:>@ when we have special needs, like
+-- non-recursive instances on inductives tuples, as used for set indices.
+
+data a :> b = !a :> !b
+  deriving (Eq,Ord,Show,Generic)
+
+derivingUnbox "StrictIxPair"
+  [t| forall a b . (Unbox a, Unbox b) => (a:>b) -> (a,b) |]
+  [| \(a:>b) -> (a, b) |]
+  [| \(a,b)  -> (a:>b) |]
+
+instance (Binary    a, Binary    b) => Binary    (a:>b)
+instance (Serialize a, Serialize b) => Serialize (a:>b)
+instance (ToJSON    a, ToJSON    b) => ToJSON    (a:>b)
+instance (FromJSON  a, FromJSON  b) => FromJSON  (a:>b)
+
+deriving instance (Read a, Read b) => Read (a:>b)
+
+instance (NFData a, NFData b) => NFData (a:>b) where
+  rnf (a:>b) = rnf a `seq` rnf b
+  {-# Inline rnf #-}
+
+--instance (Arbitrary a, Arbitrary b) => Arbitrary (a :> b) where
+--  arbitrary = (:>) <$> arbitrary <*> arbitrary
+--  shrink (a:>b) = (:>) <$> shrink a <*> shrink b
+
+
+
+-- | Base data constructor for multi-dimensional indices.
+
+data Z = Z
+  deriving (Eq,Ord,Read,Show,Generic)
+
+derivingUnbox "Z"
+  [t| Z -> () |]
+  [| const () |]
+  [| const Z  |]
+
+instance Binary    Z
+instance Serialize Z
+instance ToJSON    Z
+instance FromJSON  Z
+
+instance Arbitrary Z where
+  arbitrary = return Z
+
+instance NFData Z where
+  rnf Z = ()
+  {-# Inline rnf #-}
+
+
+
+-- | Index structures for complex, heterogeneous indexing. Mostly designed for
+-- indexing in DP grammars, where the indices work for linear and context-free
+-- grammars on one or more tapes, for strings, sets, later on tree structures.
+
+class Index i where
+
+  -- | Given a minimal size, a maximal size, and a current index, calculate
+  -- the linear index.
+
+  linearIndex :: i -> i -> i -> Int
+
+  -- | Given an index element from the smallest subset, calculate the
+  -- highest linear index that is *not* stored.
+
+  smallestLinearIndex :: i -> Int -- LH i
+
+  -- | Given an index element from the largest subset, calculate the
+  -- highest linear index that *is* stored.
+
+  largestLinearIndex :: i -> Int -- LH i
+
+  -- | Given smallest and largest index, return the number of cells
+  -- required for storage.
+
+  size :: i -> i -> Int
+
+  -- | Check if an index is within the bounds.
+
+  inBounds :: i -> i -> i -> Bool
+
+
+
+-- | Generate a stream of indices in correct order for dynamic programming.
+-- Since the stream generators require @concatMap@ / @flatten@ we have to
+-- write more specialized code for @(z:.IX)@ stuff.
+
+class IndexStream i where
+
+  -- | This generates an index stream suitable for @forward@ structure filling.
+  -- The first index is the smallest (or the first indices considered are all
+  -- equally small in partially ordered sets). Larger indices follow up until
+  -- the largest one.
+
+  streamUp   :: Monad m => i -> i -> Stream m i
+  default streamUp :: (Monad m, IndexStream (Z:.i)) => i -> i -> Stream m i
+  streamUp l h = SM.map (\(Z:.i) -> i) $ streamUp (Z:.l) (Z:.h)
+  {-# INLINE streamUp #-}
+
+  -- | If 'streamUp' generates indices from smallest to largest, then
+  -- 'streamDown' generates indices from largest to smallest. Outside grammars
+  -- make implicit use of this. Asking for an axiom in backtracking requests
+  -- the first element from this stream.
+
+  streamDown :: Monad m => i -> i -> Stream m i
+  default streamDown :: (Monad m, IndexStream (Z:.i)) => i -> i -> Stream m i
+  streamDown l h = SM.map (\(Z:.i) -> i) $ streamDown (Z:.l) (Z:.h)
+  {-# INLINE streamDown #-}
+
+
+
+instance Index Z where
+  linearIndex _ _ _ = 0
+  {-# INLINE linearIndex #-}
+  smallestLinearIndex _ = 0
+  {-# INLINE smallestLinearIndex #-}
+  largestLinearIndex _ = 0
+  {-# INLINE largestLinearIndex #-}
+  size _ _ = 1
+  {-# INLINE size #-}
+  inBounds _ _ _ = True
+  {-# INLINE inBounds #-}
+
+instance IndexStream Z where
+  streamUp   Z Z = SM.singleton Z
+  {-# INLINE streamUp #-}
+  streamDown Z Z = SM.singleton Z
+  {-# INLINE streamDown #-}
+
+instance (Index zs, Index z) => Index (zs:.z) where
+  linearIndex (ls:.l) (hs:.h) (zs:.z) = linearIndex ls hs zs * (largestLinearIndex h + 1) + linearIndex l h z
+  {-# INLINE linearIndex #-}
+  smallestLinearIndex (ls:.l) = smallestLinearIndex ls * smallestLinearIndex l
+  {-# INLINE smallestLinearIndex #-}
+  largestLinearIndex (hs:.h) = largestLinearIndex hs * largestLinearIndex h
+  {-# INLINE largestLinearIndex #-}
+  size (ls:.l) (hs:.h) = size ls hs * (size l h)
+  {-# INLINE size #-}
+  inBounds (ls:.l) (hs:.h) (zs:.z) = inBounds ls hs zs && inBounds l h z
+  {-# INLINE inBounds #-}
+
+instance (Index zs, Index z) => Index (zs:>z) where
+  linearIndex (ls:>l) (hs:>h) (zs:>z) = linearIndex ls hs zs * (largestLinearIndex h + 1) + linearIndex l h z
+  {-# INLINE linearIndex #-}
+  smallestLinearIndex (ls:>l) = smallestLinearIndex ls * smallestLinearIndex l
+  {-# INLINE smallestLinearIndex #-}
+  largestLinearIndex (hs:>h) = largestLinearIndex hs * largestLinearIndex h
+  {-# INLINE largestLinearIndex #-}
+  size (ls:>l) (hs:>h) = size ls hs * (size l h)
+  {-# INLINE size #-}
+  inBounds (ls:>l) (hs:>h) (zs:>z) = inBounds ls hs zs && inBounds l h z
+  {-# INLINE inBounds #-}
+
diff --git a/Data/PrimitiveArray/Index/Complement.hs b/Data/PrimitiveArray/Index/Complement.hs
new file mode 100644
--- /dev/null
+++ b/Data/PrimitiveArray/Index/Complement.hs
@@ -0,0 +1,60 @@
+
+module Data.PrimitiveArray.Index.Complement where
+
+import Control.Applicative
+import Control.DeepSeq (NFData(..))
+import Data.Aeson
+import Data.Binary
+import Data.Serialize
+import Data.Vector.Unboxed.Deriving
+import Data.Vector.Unboxed (Unbox(..))
+import GHC.Generics
+import Test.QuickCheck
+
+import Data.PrimitiveArray.Index.Class
+
+
+
+-- | A special index wrapper -- like @Outside@. @Complement@ allows combining
+-- inside and outside symbols which complement each other. This then yields
+-- ensemble results for each index (you need @ADPfusion@ for this).
+
+newtype Complement z = C { unC :: z }
+  deriving (Eq,Ord,Read,Show,Generic)
+
+derivingUnbox "Complement"
+  [t| forall z . Unbox z => Complement z -> z |]
+  [| unC |]
+  [| C   |]
+
+instance Binary    z => Binary    (Complement z)
+instance Serialize z => Serialize (Complement z)
+instance ToJSON    z => ToJSON    (Complement z)
+instance FromJSON  z => FromJSON  (Complement z)
+
+instance NFData z => NFData (Complement z) where
+  rnf (C z) = rnf z
+  {-# Inline rnf #-}
+
+instance Index i => Index (Complement i) where
+  linearIndex (C l) (C h) (C i) = linearIndex l h i
+  {-# INLINE linearIndex #-}
+  smallestLinearIndex (C i) = smallestLinearIndex i
+  {-# INLINE smallestLinearIndex #-}
+  largestLinearIndex (C i) = largestLinearIndex i
+  {-# INLINE largestLinearIndex #-}
+  size (C l) (C h) = size l h
+  {-# INLINE size #-}
+  inBounds (C l) (C h) (C z) = inBounds l h z
+  {-# INLINE inBounds #-}
+
+instance IndexStream i => IndexStream (Complement i) where
+  streamUp   (C l) (C h) = fmap C $ streamUp l h
+  {-# INLINE streamUp #-}
+  streamDown (C l) (C h) = fmap C $ streamDown l h
+  {-# INLINE streamDown #-}
+
+instance Arbitrary z => Arbitrary (Complement z) where
+  arbitrary    = C <$> arbitrary
+  shrink (C z) = C <$> shrink z
+
diff --git a/Data/PrimitiveArray/Index/Int.hs b/Data/PrimitiveArray/Index/Int.hs
new file mode 100644
--- /dev/null
+++ b/Data/PrimitiveArray/Index/Int.hs
@@ -0,0 +1,47 @@
+
+module Data.PrimitiveArray.Index.Int where
+
+import Data.Vector.Fusion.Stream.Monadic (flatten,map,Step(..))
+import Data.Vector.Fusion.Stream.Size
+import Prelude hiding (map)
+
+import Data.PrimitiveArray.Index.Class
+
+
+
+instance Index Int where
+  linearIndex _ _ k = k
+  {-# Inline linearIndex #-}
+  smallestLinearIndex _ = error "still needed?"
+  {-# Inline smallestLinearIndex #-}
+  largestLinearIndex h = h
+  {-# Inline largestLinearIndex #-}
+  size _ h = h+1
+  {-# Inline size #-}
+  inBounds l h k = l <= k && k <= h
+  {-# Inline inBounds #-}
+
+instance IndexStream z => IndexStream (z:.Int) where
+  streamUp (ls:.l) (hs:.h) = flatten mk step Unknown $ streamUp ls hs
+    where mk z = return (z,l)
+          step (z,k)
+            | k > h     = return $ Done
+            | otherwise = return $ Yield (z:.k) (z,k+1)
+          {-# Inline [0] mk   #-}
+          {-# Inline [0] step #-}
+  {-# Inline streamUp #-}
+  streamDown (ls:.l) (hs:.h) = flatten mk step Unknown $ streamDown ls hs
+    where mk z = return (z,h)
+          step (z,k)
+            | k < l     = return $ Done
+            | otherwise = return $ Yield (z:.k) (z,k-1)
+          {-# Inline [0] mk   #-}
+          {-# Inline [0] step #-}
+  {-# Inline streamDown #-}
+
+instance IndexStream Int where
+  streamUp l h = map (\(Z:.k) -> k) $ streamUp (Z:.l) (Z:.h)
+  {-# Inline streamUp #-}
+  streamDown l h = map (\(Z:.k) -> k) $ streamDown (Z:.l) (Z:.h)
+  {-# Inline streamDown #-}
+
diff --git a/Data/PrimitiveArray/Index/Outside.hs b/Data/PrimitiveArray/Index/Outside.hs
new file mode 100644
--- /dev/null
+++ b/Data/PrimitiveArray/Index/Outside.hs
@@ -0,0 +1,62 @@
+
+module Data.PrimitiveArray.Index.Outside where
+
+import Control.Applicative
+import Control.DeepSeq (NFData(..))
+import Data.Aeson
+import Data.Binary
+import Data.Serialize
+import Data.Vector.Unboxed.Deriving
+import Data.Vector.Unboxed (Unbox(..))
+import GHC.Generics
+import Test.QuickCheck
+
+import Data.PrimitiveArray.Index.Class
+
+
+
+-- | The 'Outside' wrapper takes an index structure, and provides
+-- 'IndexStream' functions 'streamUp' and 'streamDown' that work the other
+-- way around. In particular, for @Outside z@ @streamUp (Outside z) = fmap
+-- Outside $ streamDown z@ and vice versa. @Index@ functions are unwrapped
+-- but otherwise work as before.
+
+newtype Outside z = O { unO :: z }
+  deriving (Eq,Ord,Read,Show,Generic)
+
+derivingUnbox "Outside"
+  [t| forall z . Unbox z => Outside z -> z |]
+  [| unO |]
+  [| O   |]
+
+instance Binary    z => Binary    (Outside z)
+instance Serialize z => Serialize (Outside z)
+instance ToJSON    z => ToJSON    (Outside z)
+instance FromJSON  z => FromJSON  (Outside z)
+
+instance NFData z => NFData (Outside z) where
+  rnf (O z) = rnf z
+  {-# Inline rnf #-}
+
+instance Index i => Index (Outside i) where
+  linearIndex (O l) (O h) (O i) = linearIndex l h i
+  {-# INLINE linearIndex #-}
+  smallestLinearIndex (O i) = smallestLinearIndex i
+  {-# INLINE smallestLinearIndex #-}
+  largestLinearIndex (O i) = largestLinearIndex i
+  {-# INLINE largestLinearIndex #-}
+  size (O l) (O h) = size l h
+  {-# INLINE size #-}
+  inBounds (O l) (O h) (O z) = inBounds l h z
+  {-# INLINE inBounds #-}
+
+instance IndexStream i => IndexStream (Outside i) where
+  streamUp (O l) (O h) = fmap O $ streamDown l h
+  {-# INLINE streamUp #-}
+  streamDown (O l) (O h) = fmap O $ streamUp l h
+  {-# INLINE streamDown #-}
+
+instance Arbitrary z => Arbitrary (Outside z) where
+  arbitrary = O <$> arbitrary
+  shrink (O z) = O <$> shrink z
+
diff --git a/Data/PrimitiveArray/Index/Point.hs b/Data/PrimitiveArray/Index/Point.hs
new file mode 100644
--- /dev/null
+++ b/Data/PrimitiveArray/Index/Point.hs
@@ -0,0 +1,119 @@
+
+-- | @Point@ index structures are used for left- and right-linear grammars.
+-- Such grammars have at most one syntactic symbol on each r.h.s. of a rule.
+-- The syntactic symbol needs to be in an outermost position.
+
+module Data.PrimitiveArray.Index.Point where
+
+import           Control.Applicative
+import           Control.DeepSeq (NFData(..))
+import           Data.Aeson
+import           Data.Binary
+import           Data.Bits
+import           Data.Bits.Extras (Ranked)
+import           Data.Serialize
+import           Data.Vector.Fusion.Stream.Size
+import           Data.Vector.Unboxed.Deriving
+import           Data.Vector.Unboxed (Unbox(..))
+import           GHC.Generics
+import qualified Data.Vector.Fusion.Stream.Monadic as SM
+import qualified Data.Vector.Unboxed as VU
+import           Test.QuickCheck
+
+import           Data.PrimitiveArray.Index.Class
+
+
+
+-- | A point in a left-linear grammar. The syntactic symbol is in left-most
+-- position.
+
+newtype PointL = PointL {fromPointL :: Int}
+  deriving (Eq,Read,Show,Generic)
+
+-- | A point in a right-linear grammars.
+
+newtype PointR = PointR {fromPointR :: Int}
+  deriving (Eq,Read,Show,Generic)
+
+
+
+derivingUnbox "PointL"
+  [t| PointL -> Int    |]
+  [| \ (PointL i) -> i |]
+  [| \ i -> PointL i   |]
+
+instance Binary    PointL
+instance Serialize PointL
+instance FromJSON  PointL
+instance ToJSON    PointL
+
+instance NFData PointL where
+  rnf (PointL l) = rnf l
+  {-# Inline rnf #-}
+
+instance Index PointL where
+  linearIndex _ _ (PointL z) = z
+  {-# INLINE linearIndex #-}
+  smallestLinearIndex (PointL l) = error "still needed?"
+  {-# INLINE smallestLinearIndex #-}
+  largestLinearIndex (PointL h) = h
+  {-# INLINE largestLinearIndex #-}
+  size (_) (PointL h) = h + 1
+  {-# INLINE size #-}
+  inBounds (_) (PointL h) (PointL x) = 0<=x && x<=h
+  {-# INLINE inBounds #-}
+
+instance IndexStream z => IndexStream (z:.PointL) where
+  streamUp (ls:.PointL lf) (hs:.PointL ht) = SM.flatten mk step Unknown $ streamUp ls hs
+    where mk z = return (z,lf)
+          step (z,k)
+            | k > ht    = return $ SM.Done
+            | otherwise = return $ SM.Yield (z:.PointL k) (z,k+1)
+          {-# Inline [0] mk   #-}
+          {-# Inline [0] step #-}
+  {-# Inline streamUp #-}
+  streamDown (ls:.PointL lf) (hs:.PointL ht) = SM.flatten mk step Unknown $ streamDown ls hs
+    where mk z = return (z,ht)
+          step (z,k)
+            | k < lf    = return $ SM.Done
+            | otherwise = return $ SM.Yield (z:.PointL k) (z,k-1)
+          {-# Inline [0] mk   #-}
+          {-# Inline [0] step #-}
+  {-# Inline streamDown #-}
+
+instance IndexStream PointL
+
+instance Arbitrary PointL where
+  arbitrary = do
+    b <- choose (0,100)
+    return $ PointL b
+  shrink (PointL j)
+    | 0<j = [PointL $ j-1]
+    | otherwise = []
+
+
+
+derivingUnbox "PointR"
+  [t| PointR -> Int    |]
+  [| \ (PointR i) -> i |]
+  [| \ i -> PointR i   |]
+
+instance Binary    PointR
+instance Serialize PointR
+instance FromJSON  PointR
+instance ToJSON    PointR
+
+instance NFData PointR where
+  rnf (PointR l) = rnf l
+  {-# Inline rnf #-}
+
+instance Index PointR where
+  linearIndex l _ (PointR z) = undefined
+  {-# INLINE linearIndex #-}
+  smallestLinearIndex = undefined
+  {-# INLINE smallestLinearIndex #-}
+  largestLinearIndex = undefined
+  {-# INLINE largestLinearIndex #-}
+  size = undefined
+  {-# INLINE size #-}
+
diff --git a/Data/PrimitiveArray/Index/Set.hs b/Data/PrimitiveArray/Index/Set.hs
new file mode 100644
--- /dev/null
+++ b/Data/PrimitiveArray/Index/Set.hs
@@ -0,0 +1,508 @@
+
+-- | Set with and without interfaces. We provide instances for sets, and
+-- sets with one or two interfaces. The @First@ and @Last@ annotation is
+-- purely cosmetical (apart from introducing type safety).
+
+module Data.PrimitiveArray.Index.Set where
+
+import           Control.Applicative ((<$>),(<*>))
+import           Control.DeepSeq (NFData(..))
+import           Data.Aeson (FromJSON,ToJSON)
+import           Data.Binary (Binary)
+import           Data.Bits
+import           Data.Bits.Extras
+import           Data.Serialize (Serialize)
+import           Data.Vector.Fusion.Stream.Size
+import           Data.Vector.Unboxed.Deriving
+import           Data.Vector.Unboxed (Unbox(..))
+import           Debug.Trace
+import           GHC.Generics
+import qualified Data.Vector.Fusion.Stream.Monadic as SM
+import qualified Data.Vector.Unboxed as VU
+import           Test.QuickCheck (Arbitrary(..), choose, elements)
+
+import           Data.Bits.Ordered
+import           Data.PrimitiveArray.Index.Class
+
+
+
+-- * @newtype@s, @data@ types, @class@es.
+
+
+
+-- | Certain sets have an interface, a particular element with special
+-- meaning. In this module, certain ``meanings'' are already provided.
+-- These include a @First@ element and a @Last@ element. We phantom-type
+-- these to reduce programming overhead.
+
+newtype Interface t = Iter { getIter :: Int }
+  deriving (Eq,Ord,Read,Show,Generic,Num)
+
+-- | Declare the interface to be the start of a path.
+
+data First
+
+-- | Declare the interface to be the end of a path.
+
+data Last
+
+-- | Declare the interface to match anything.
+--
+-- TODO needed? want to use later in ADPfusion
+
+data Any
+
+-- | Newtype for a bitset. We'd use @Word@s but that requires more shape
+-- instances.
+--
+-- TODO can we use @Word@s now?
+
+newtype BitSet = BitSet { getBitSet :: Int }
+  deriving (Eq,Ord,Read,Generic,FiniteBits,Ranked,Num,Bits)
+
+-- | A bitset with one interface.
+
+type BS1I i = BitSet:>Interface i
+
+-- | A bitset with two interfaces.
+
+type BS2I i j = BitSet:>Interface i:>Interface j
+
+-- | Successor and Predecessor for sets. Designed as a class to accomodate
+-- sets with interfaces and without interfaces with one function.
+--
+-- The functions are not written recursively, as we currently only have
+-- three cases, and we do not want to "reset" while generating successors
+-- and predecessors.
+--
+-- Note that sets have a partial order. Within the group of element with
+-- the same @popCount@, we use @popPermutation@ which has the same stepping
+-- order for both, @setSucc@ and @setPred@.
+
+class SetPredSucc s where
+  -- | Set successor. The first argument is the lower set limit, the second
+  -- the upper set limit, the third the current set.
+  setSucc :: s -> s -> s -> Maybe s
+  -- | Set predecessor. The first argument is the lower set limit, the
+  -- second the upper set limit, the third the current set.
+  setPred :: s -> s -> s -> Maybe s
+
+-- | Masks are used quite often for different types of bitsets. We liberate
+-- them as a type family.
+
+type family Mask s :: *
+
+-- | @Fixed@ allows us to fix some or all bits of a bitset, thereby
+-- providing @succ/pred@ operations which are only partially free.
+--
+-- The mask is lazy, this allows us to have @undefined@ for @l@ and @h@.
+--
+-- @f = getFixedMask .&. getFixed@ are the fixed bits.
+-- @n = getFixed .&. complement getFixedMask@ are the free bits.
+-- @to = complement getFixed@ is the to move mask
+-- @n' = popShiftR to n@ yields the population after the move
+-- @p = popPermutation undefined n'@ yields the new population permutation
+-- @p' = popShiftL to p@ yields the population moved back
+-- @final = p' .|. f@
+
+data Fixed t = Fixed { getFixedMask :: (Mask t) , getFixed :: !t }
+
+-- | Assuming a bitset on bits @[0 .. highbit]@, we can apply a mask that
+-- stretches out those bits over @[0 .. higherBit]@ with @highbit <=
+-- higherBit@. Any active interfaces are correctly set as well.
+
+class ApplyMask s where
+  applyMask :: Mask s -> s -> s
+
+
+
+-- * Instances
+
+
+
+derivingUnbox "Interface"
+  [t| forall t . Interface t -> Int |]
+  [| \(Iter i) -> i            |]
+  [| Iter                      |]
+
+instance Binary    (Interface t)
+instance Serialize (Interface t)
+instance ToJSON    (Interface t)
+instance FromJSON  (Interface t)
+
+instance NFData (Interface t) where
+  rnf (Iter i) = rnf i
+  {-# Inline rnf #-}
+
+instance Index (Interface i) where
+  linearIndex l _ (Iter z) = z - smallestLinearIndex l
+  {-# INLINE linearIndex #-}
+  smallestLinearIndex (Iter l) = l
+  {-# INLINE smallestLinearIndex #-}
+  largestLinearIndex (Iter h) = h
+  {-# INLINE largestLinearIndex #-}
+  size (Iter l) (Iter h) = h - l + 1
+  {-# INLINE size #-}
+  inBounds l h z = l <= z && z <= h
+  {-# INLINE inBounds #-}
+
+
+
+derivingUnbox "BitSet"
+  [t| BitSet     -> Int |]
+  [| \(BitSet s) -> s   |]
+  [| BitSet             |]
+
+instance Show BitSet where
+  show (BitSet s) = "<" ++ (show $ activeBitsL s) ++ ">(" ++ show s ++ ")"
+
+instance Binary    BitSet
+instance Serialize BitSet
+instance ToJSON    BitSet
+instance FromJSON  BitSet
+
+instance NFData BitSet where
+  rnf (BitSet s) = rnf s
+  {-# Inline rnf #-}
+
+instance Index BitSet where
+  linearIndex l _ (BitSet z) = z - smallestLinearIndex l -- (2 ^ popCount l - 1)
+  {-# INLINE linearIndex #-}
+  smallestLinearIndex l = 2 ^ popCount l - 1
+  {-# INLINE smallestLinearIndex #-}
+  largestLinearIndex h = 2 ^ popCount h - 1
+  {-# INLINE largestLinearIndex #-}
+  size l h = 2 ^ popCount h - 2 ^ popCount l + 1
+  {-# INLINE size #-}
+  inBounds l h z = popCount l <= popCount z && popCount z <= popCount h
+  {-# INLINE inBounds #-}
+
+
+
+instance IndexStream z => IndexStream (z:.BitSet) where
+  streamUp (ls:.l) (hs:.h) = SM.flatten mk step Unknown $ streamUp ls hs
+    where mk z = return (z , (if l <= h then Just l else Nothing))
+          step (z , Nothing) = return $ SM.Done
+          step (z , Just t ) = return $ SM.Yield (z:.t) (z , setSucc l h t)
+          {-# Inline [0] mk   #-}
+          {-# Inline [0] step #-}
+  {-# Inline streamUp   #-}
+  streamDown (ls:.l) (hs:.h) = SM.flatten mk step Unknown $ streamDown ls hs
+    where mk z = return (z :. (if l <= h then Just h else Nothing))
+          step (z :. Nothing) = return $ SM.Done
+          step (z :. Just t ) = return $ SM.Yield (z:.t) (z :. setPred l h t)
+          {-# Inline [0] mk   #-}
+          {-# Inline [0] step #-}
+  {-# Inline streamDown #-}
+
+instance IndexStream z => IndexStream (z:.(BitSet:>Interface i)) where
+  streamUp (ls:.l@(sl:>_)) (hs:.h@(sh:>_)) = SM.flatten mk step Unknown $ streamUp ls hs
+    where mk z = return (z, (if sl<=sh then Just (sl:>(Iter . max 0 $ lsbZ sl)) else Nothing))
+          step (z , Nothing) = return $ SM.Done
+          step (z,  Just t ) = return $ SM.Yield (z:.t) (z , setSucc l h t)
+          {-# Inline [0] mk   #-}
+          {-# Inline [0] step #-}
+  {-# Inline streamUp #-}
+  streamDown (ls:.l@(sl:>_)) (hs:.h@(sh:>_)) = SM.flatten mk step Unknown $ streamDown ls hs
+    where mk z = return (z, (if sl<=sh then Just (sh:>(Iter . max 0 $ lsbZ sh)) else Nothing))
+          step (z , Nothing) = return $ SM.Done
+          step (z , Just t ) = return $ SM.Yield (z:.t) (z , setPred l h t)
+          {-# Inline [0] mk   #-}
+          {-# Inline [0] step #-}
+  {-# Inline streamDown #-}
+
+instance IndexStream z => IndexStream (z:.(BitSet:>Interface i:>Interface j)) where
+  streamUp (ls:.l@(sl:>_:>_)) (hs:.h@(sh:>_:>_)) = SM.flatten mk step Unknown $ streamUp ls hs
+    where mk z | sl > sh   = return (z , Nothing)
+               | cl == 0   = return (z , Just (0:>0:>0))
+               | cl == 1   = let i = lsbZ sl
+                             in  return (z , Just (sl :> Iter i :> Iter i))
+               | otherwise = let i = lsbZ sl; j = lsbZ (sl `clearBit` i)
+                             in  return (z , Just (sl :> Iter i :> Iter j))
+               where cl = popCount sl
+          step (z , Nothing) = return $ SM.Done
+          step (z , Just t ) = return $ SM.Yield (z:.t) (z , setSucc l h t)
+          {-# Inline [0] mk   #-}
+          {-# Inline [0] step #-}
+  {-# Inline streamUp #-}
+  streamDown (ls:.l@(sl:>_:>_)) (hs:.h@(sh:>_:>_)) = SM.flatten mk step Unknown $ streamDown ls hs
+    where mk z | sl > sh   = return (z , Nothing)
+               | ch == 0   = return (z , Just (0:>0:>0))
+               | ch == 1   = let i = lsbZ sh
+                             in  return (z , Just (sh :> Iter i :> Iter i))
+               | otherwise = let i = lsbZ sh; j = lsbZ sh
+                             in  return (z , Just (sh :> Iter i :> Iter j))
+               where ch = popCount sh
+          step (z , Nothing) = return $ SM.Done
+          step (z , Just t ) = return $ SM.Yield (z:.t) (z , setPred l h t)
+          {-# Inline [0] mk   #-}
+          {-# Inline [0] step #-}
+  {-# Inline streamDown #-}
+
+
+
+instance SetPredSucc BitSet where
+  setSucc l h s
+    | cs > ch                        = Nothing
+    | Just s' <- popPermutation ch s = Just s'
+    | cs >= ch                       = Nothing
+    | cs < ch                        = Just . BitSet $ 2^(cs+1) -1
+    where ch = popCount h
+          cs = popCount s
+  {-# Inline setSucc #-}
+  setPred l h s
+    | cs < cl                        = Nothing
+    | Just s' <- popPermutation ch s = Just s'
+    | cs <= cl                       = Nothing
+    | cs > cl                        = Just . BitSet $ 2^(cs-1) -1
+    where cl = popCount l
+          ch = popCount h
+          cs = popCount s
+  {-# Inline setPred #-}
+
+instance SetPredSucc (BitSet:>Interface i) where
+  setSucc (l:>il) (h:>ih) (s:>Iter is)
+    | cs > ch                         = Nothing
+    | Just is' <- maybeNextActive is s     = Just (s:>Iter is')
+    | Just s'  <- popPermutation ch s = Just (s':>Iter (lsbZ s'))
+    | cs >= ch                        = Nothing
+    | cs < ch                         = let s' = BitSet $ 2^(cs+1)-1 in Just (s' :> Iter (lsbZ s'))
+    where ch = popCount h
+          cs = popCount s
+  {-# Inline setSucc #-}
+  setPred (l:>il) (h:>ih) (s:>Iter is)
+    | cs < cl                         = Nothing
+    | Just is' <- maybeNextActive is s     = Just (s:>Iter is')
+    | Just s'  <- popPermutation ch s = Just (s':>Iter (lsbZ s'))
+    | cs <= cl                        = Nothing
+    | cs > cl                         = let s' = BitSet $ 2^(cs-1)-1 in Just (s' :> Iter (max 0 $ lsbZ s'))
+    where cl = popCount l
+          ch = popCount h
+          cs = popCount s
+  {-# Inline setPred #-}
+
+instance SetPredSucc (BitSet:>Interface i:>Interface j) where
+  setSucc (l:>il:>jl) (h:>ih:>jh) (s:>Iter is:>Iter js)
+    -- early termination
+    | cs > ch                         = Nothing
+    -- in case nothing was set, set initial set @1@ with both interfaces
+    -- pointing to the same element
+    | cs == 0                         = Just (1:>0:>0)
+    -- when only a single element is set, we just permute the population
+    -- and set the single interface
+    | cs == 1
+    , Just s'  <- popPermutation ch s
+    , let is' = lsbZ s'          = Just (s':>Iter is':>Iter is')
+    -- try advancing only one of the interfaces, doesn't collide with @is@
+    | Just js' <- maybeNextActive js (s `clearBit` is) = Just (s:>Iter is:>Iter js')
+    -- advance other interface, 
+    | Just is' <- maybeNextActive is s
+    , let js' = lsbZ (s `clearBit` is')      = Just (s:>Iter is':>Iter js')
+    -- find another permutation of the population
+    | Just s'  <- popPermutation ch s
+    , let is' = lsbZ s'
+    , Just js' <- maybeNextActive is' s'   = Just (s':>Iter is':>Iter js')
+    -- increasing the population forbidden by upper limit
+    | cs >= ch                        = Nothing
+    -- increase population
+    | cs < ch
+    , let s' = BitSet $ 2^(cs+1)-1
+    , let is' = lsbZ s'
+    , Just js' <- maybeNextActive is' s'   = Just (s':>Iter is':>Iter js')
+    where ch = popCount h
+          cs = popCount s
+  {-# Inline setSucc #-}
+  setPred (l:>il:>jl) (h:>ih:>jh) (s:>Iter is:>Iter js)
+    -- early termination
+    | cs < cl                         = Nothing
+    -- in case nothing was set, set initial set @1@ with both interfaces
+    -- pointing to the same element
+    | cs == 0                         = Nothing
+    -- when only a single element is set, we just permute the population
+    -- and set the single interface
+    | cs == 1
+    , Just s'  <- popPermutation ch s
+    , let is' = lsbZ s'          = Just (s':>Iter is':>Iter is')
+    -- return the single @0@ set
+    | cs == 1                         = Just (0:>0:>0)
+    -- try advancing only one of the interfaces, doesn't collide with @is@
+    | Just js' <- maybeNextActive js (s `clearBit` is) = Just (s:>Iter is:>Iter js')
+    -- advance other interface, 
+    | Just is' <- maybeNextActive is s
+    , let js' = lsbZ (s `clearBit` is')      = Just (s:>Iter is':>Iter js')
+    -- find another permutation of the population
+    | Just s'  <- popPermutation ch s
+    , let is' = lsbZ s'
+    , Just js' <- maybeNextActive is' s'   = Just (s':>Iter is':>Iter js')
+    -- decreasing the population forbidden by upper limit
+    | cs <= cl                        = Nothing
+    -- decrease population
+    | cs > cl && cs > 2
+    , let s' = BitSet $ 2^(cs-1)-1
+    , let is' = lsbZ s'
+    , Just js' <- maybeNextActive is' s'   = Just (s':>Iter is':>Iter js')
+    -- decrease population to single-element sets
+    | cs > cl && cs == 2              = Just (1:>0:>0)
+    where cl = popCount l
+          ch = popCount h
+          cs = popCount s
+  {-# Inline setPred #-}
+
+
+
+type instance Mask BitSet = BitSet
+
+type instance Mask (BitSet :> Interface i) = BitSet
+
+type instance Mask (BitSet :> Interface i :> Interface j) = BitSet
+
+
+
+derivingUnbox "Fixed"
+  [t| forall t . (Unbox t, Unbox (Mask t)) => Fixed t -> (Mask t, t) |]
+  [| \(Fixed m s) -> (m,s)              |]
+  [| uncurry Fixed                      |]
+
+deriving instance (Eq t     , Eq      (Mask t)) => Eq      (Fixed t)
+deriving instance (Ord t    , Ord     (Mask t)) => Ord     (Fixed t)
+deriving instance (Read t   , Read    (Mask t)) => Read    (Fixed t)
+deriving instance (Show t   , Show    (Mask t)) => Show    (Fixed t)
+deriving instance (Generic t, Generic (Mask t)) => Generic (Fixed t)
+
+instance (Generic t, Generic (Mask t), Binary t   , Binary    (Mask t)) => Binary    (Fixed t)
+instance (Generic t, Generic (Mask t), Serialize t, Serialize (Mask t)) => Serialize (Fixed t)
+{- -- TODO do json instances work automatically here?
+instance ToJSON    (Interface t)
+instance FromJSON  (Interface t)
+-}
+
+instance NFData (Fixed t) where
+  rnf (Fixed m s) = m `seq` s `seq` ()
+
+-- TODO we need to be careful here, that we actually fix all bits that are
+-- fixed AND that during permutations / increases in popCount we do not set
+-- an already fixed bit -- as otherwise we lose one in popCount.
+
+testBsS :: BitSet -> Maybe (Fixed BitSet)
+testBsS k = setSucc (Fixed 0 0) (Fixed 0 7) (Fixed 4 k)
+{-# NoInline testBsS #-}
+
+instance SetPredSucc (Fixed BitSet) where
+  setPred (Fixed _ l) (Fixed _ h) (Fixed !m s) = Fixed m <$> setPred l h (s .&. complement m)
+  {-# Inline setPred #-}
+  --setSucc (Fixed _ l) (Fixed _ h) (Fixed !m s) = Fixed m <$> setSucc l h (s .&. complement m)
+  --setSucc (Fixed _ l) (Fixed _ h) (Fixed !m' s) = (Fixed m . (.|. f)) <$> p -- return population, now again including the fixed part @f@
+  --  where m = m' .&. h            -- constrain the mask to just the bits until @h@
+  --        f = s .&. m             -- these bits are fixed to @1@
+  --        n = s .&. complement m  -- these bits are free to be @0@ or @1@ and may move around; this means that @n `subset` complement m@
+  --        to = complement m       -- once we have calculated our permutation, we move it to the correct places via @to@
+  --        n' = popShiftR to n     -- population without holes. all primes denote that we are in hole-free space.
+  --        p' = popPermutation (popCount $ h .&. to) n'  -- permutate the shifted population
+  --        p  = popShiftL to <$> p'  -- undo the shift
+  setSucc (Fixed _ l) (Fixed _ h) (Fixed !m' s) = traceShow (h,m,s,' ',fb0,fb1,' ',p',p'',p) $ (Fixed m . (.|. fb1)) <$> p
+    where m   = m' .&. h
+          fb0 = m  .&. complement s
+          fb1 = m  .&. s
+          p'  = popShiftR m s
+          p'' = setSucc (popShiftR m l) (popShiftR m h) p'
+          p   = popShiftL m <$> p''
+  {-# Inline setSucc #-}
+
+instance SetPredSucc (Fixed (BitSet:>Interface i)) where
+  setPred (Fixed _ (l:>li)) (Fixed _ (h:>hi)) (Fixed !m (s:>i))
+    | s `testBit` getIter i = (Fixed m . (:> i) . ( `setBit` getIter i)) <$> setPred l h (s .&. complement m)
+    | otherwise             = (Fixed m) <$> setPred (l:>li) (h:>hi) ((s .&. complement m):>i)
+  {-# Inline setPred #-}
+  setSucc (Fixed _ (l:>li)) (Fixed _ (h:>hi)) (Fixed !m (s:>i))
+    | s `testBit` getIter i = (Fixed m . (:> i) . ( `setBit` getIter i)) <$> setSucc l h (s .&. complement m)
+    | otherwise             = (Fixed m) <$> setSucc (l:>li) (h:>hi) ((s .&. complement m):>i)
+  {-# Inline setSucc #-}
+
+instance SetPredSucc (Fixed (BitSet:>Interface i:>Interface j)) where
+  setPred (Fixed _ (l:>li:>lj)) (Fixed _ (h:>hi:>hj)) (Fixed !m (s:>i:>j))
+    | s `testBit` getIter i && s `testBit` getIter j
+    = (Fixed m . (\z       -> (z `setBit` getIter i `setBit` getIter j:>i:>j ))) <$> setPred l h (s .&. complement m)
+    | s `testBit` getIter i
+    = (Fixed m . (\(z:>j') -> (z `setBit` getIter i                   :>i:>j'))) <$> setPred (l:>lj) (h:>hj) (s .&. complement m :>j)
+    | s `testBit` getIter j
+    = (Fixed m . (\(z:>i') -> (z `setBit` getIter j                   :>i':>j))) <$> setPred (l:>li) (h:>hi) (s .&. complement m :>i)
+  {-# Inline setPred #-}
+  setSucc (Fixed _ (l:>li:>lj)) (Fixed _ (h:>hi:>hj)) (Fixed !m (s:>i:>j))
+    | s `testBit` getIter i && s `testBit` getIter j
+    = (Fixed m . (\z       -> (z `setBit` getIter i `setBit` getIter j:>i:>j ))) <$> setSucc l h (s .&. complement m)
+    | s `testBit` getIter i
+    = (Fixed m . (\(z:>j') -> (z `setBit` getIter i                   :>i:>j'))) <$> setSucc (l:>lj) (h:>hj) (s .&. complement m :>j)
+    | s `testBit` getIter j
+    = (Fixed m . (\(z:>i') -> (z `setBit` getIter j                   :>i':>j))) <$> setSucc (l:>li) (h:>hi) (s .&. complement m :>i)
+  {-# Inline setSucc #-}
+
+
+
+instance ApplyMask BitSet where
+  applyMask = popShiftL
+  {-# Inline applyMask #-}
+
+instance ApplyMask (BitSet :> Interface i) where
+  applyMask m (s:>i)
+    | popCount s == 0 = 0:>0
+    | otherwise       = popShiftL m s :> (Iter . getBitSet . popShiftL m . BitSet $ 2 ^ getIter i)
+  {-# Inline applyMask #-}
+
+instance ApplyMask (BitSet :> Interface i :> Interface j) where
+  applyMask m (s:>i:>j)
+    | popCount s == 0 = 0:>0:>0
+    | popCount s == 1 = s' :> i' :> Iter (getIter i')
+    | otherwise       = s' :> i' :> j'
+    where s' = popShiftL m s
+          i' = Iter . getBitSet . popShiftL m . BitSet $ 2 ^ getIter i
+          j' = Iter . getBitSet . popShiftL m . BitSet $ 2 ^ getIter j
+  {-# Inline applyMask #-}
+
+
+
+arbitraryBitSetMax = 6
+
+instance (Arbitrary t, Arbitrary (Mask t)) => Arbitrary (Fixed t) where
+  arbitrary = Fixed <$> arbitrary <*> arbitrary
+  shrink (Fixed m s) = [ Fixed m' s' | m' <- shrink m, s' <- shrink s ]
+
+instance Arbitrary BitSet where
+  arbitrary = BitSet <$> choose (0,2^arbitraryBitSetMax-1)
+  shrink s = let s' = [ s `clearBit` a | a <- activeBitsL s ]
+             in  s' ++ concatMap shrink s'
+
+instance Arbitrary (BitSet:>Interface i) where
+  arbitrary = do
+    s <- arbitrary
+    if s==0
+      then return (s:>Iter 0)
+      else do i <- elements $ activeBitsL s
+              return (s:>Iter i)
+  shrink (s:>i) =
+    let s' = [ (s `clearBit` a:>i)
+             | a <- activeBitsL s
+             , Iter a /= i ]
+             ++ [ 0 :> Iter 0 | popCount s == 1 ]
+    in  s' ++ concatMap shrink s'
+
+instance Arbitrary (BitSet:>Interface i:>Interface j) where
+  arbitrary = do
+    s <- arbitrary
+    case (popCount s) of
+      0 -> return (s:>Iter 0:>Iter 0)
+      1 -> do i <- elements $ activeBitsL s
+              return (s:>Iter i:>Iter i)
+      _ -> do i <- elements $ activeBitsL s
+              j <- elements $ activeBitsL (s `clearBit` i)
+              return (s:>Iter i:>Iter j)
+  shrink (s:>i:>j) =
+    let s' = [ (s `clearBit` a:>i:>j)
+             | a <- activeBitsL s
+             , Iter a /= i, Iter a /= j ]
+             ++ [ 0 `setBit` a :> Iter a :> Iter a
+                | popCount s == 2
+                , a <- activeBitsL s ]
+             ++ [ 0 :> Iter 0 :> Iter 0
+                | popCount s == 1 ]
+    in  s' ++ concatMap shrink s'
+
diff --git a/Data/PrimitiveArray/Index/Subword.hs b/Data/PrimitiveArray/Index/Subword.hs
new file mode 100644
--- /dev/null
+++ b/Data/PrimitiveArray/Index/Subword.hs
@@ -0,0 +1,132 @@
+
+-- | Index structure for context-free grammars on strings. A @Subword@ captures
+-- a pair @(i,j)@ with @i<=j@.
+
+module Data.PrimitiveArray.Index.Subword where
+
+import Control.DeepSeq (NFData(..))
+import Data.Aeson (FromJSON,ToJSON)
+import Data.Binary (Binary)
+import Data.Serialize (Serialize)
+import Data.Vector.Fusion.Stream.Monadic (Step(..), flatten, map)
+import Data.Vector.Fusion.Stream.Size
+import Data.Vector.Unboxed.Deriving
+import GHC.Generics (Generic)
+import Test.QuickCheck (Arbitrary(..), choose)
+import Prelude hiding (map)
+
+import Data.PrimitiveArray.Index.Class
+
+
+
+-- | A subword wraps a pair of @Int@ indices @i,j@ with @i<=j@.
+--
+-- 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) ... (1,1) ... (k,k) ... (N,N)@ 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 {fromSubword :: (Int:.Int)}
+  deriving (Eq,Ord,Show,Generic,Read)
+
+derivingUnbox "Subword"
+  [t| Subword -> (Int,Int) |]
+  [| \ (Subword (i:.j)) -> (i,j) |]
+  [| \ (i,j) -> Subword (i:.j) |]
+
+instance Binary    Subword
+instance Serialize Subword
+instance FromJSON  Subword
+instance ToJSON    Subword
+
+instance NFData Subword where
+  rnf (Subword (i:.j)) = i `seq` rnf j
+  {-# Inline rnf #-}
+
+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+1
+{-# 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+1)*i - triangularNumber i + j
+{-# INLINE subwordIndex #-}
+
+subwordFromIndex :: Subword -> Int -> Subword
+subwordFromIndex = error "subwordFromIndex not implemented"
+{-# INLINE subwordFromIndex #-}
+
+
+
+instance Index Subword where
+  linearIndex _ h i = subwordIndex h i
+  {-# Inline linearIndex #-}
+  smallestLinearIndex _ = error "still needed?"
+  {-# Inline smallestLinearIndex #-}
+  largestLinearIndex = upperTri
+  {-# Inline largestLinearIndex #-}
+  size _ h = upperTri h
+  {-# Inline size #-}
+  inBounds _ (Subword (_:.h)) (Subword (i:.j)) = 0<=i && i<=j && j<=h
+  {-# Inline inBounds #-}
+
+instance IndexStream z => IndexStream (z:.Subword) where
+  streamUp (ls:.Subword (l:._)) (hs:.Subword (_:.h)) = flatten mk step Unknown $ streamUp ls hs
+    where mk z = return (z,h,h)
+          step (z,i,j)
+            | i < l     = return $ Done
+            | j > h     = return $ Skip (z,i-1,i-1)
+            | otherwise = return $ Yield (z:.subword i j) (z,i,j+1)
+          {-# Inline [0] mk   #-}
+          {-# Inline [0] step #-}
+  {-# Inline streamUp #-}
+  streamDown (ls:.Subword (l:._)) (hs:.Subword (_:.h)) = flatten mk step Unknown $ streamDown ls hs
+    where mk z = return (z,l,h)
+          step (z,i,j)
+            | i > h     = return $ Done
+            | j < i     = return $ Skip (z,i+1,h)
+            | otherwise = return $ Yield (z:.subword i j) (z,i,j-1)
+          {-# Inline [0] mk   #-}
+          {-# Inline [0] step #-}
+  {-# Inline streamDown #-}
+
+-- Default methods don't inline in a good way!
+
+instance IndexStream Subword where
+  streamUp l h = map (\(Z:.i) -> i) $ streamUp (Z:.l) (Z:.h)
+  {-# INLINE streamUp #-}
+  streamDown l h = map (\(Z:.i) -> i) $ streamDown (Z:.l) (Z:.h)
+  {-# INLINE streamDown #-}
+
+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), Subword (i+1:.j)]
+    | otherwise = []
+
diff --git a/Data/PrimitiveArray/QuickCheck.hs b/Data/PrimitiveArray/QuickCheck.hs
deleted file mode 100644
--- a/Data/PrimitiveArray/QuickCheck.hs
+++ /dev/null
@@ -1,5 +0,0 @@
-
-module Data.PrimitiveArray.QuickCheck where
-
-
-
diff --git a/Data/PrimitiveArray/QuickCheck/Index/Set.hs b/Data/PrimitiveArray/QuickCheck/Index/Set.hs
new file mode 100644
--- /dev/null
+++ b/Data/PrimitiveArray/QuickCheck/Index/Set.hs
@@ -0,0 +1,31 @@
+
+module Data.PrimitiveArray.QuickCheck.Index.Set where
+
+import Control.Applicative
+import Data.Bits
+import Data.Word (Word)
+import Debug.Trace
+import Test.QuickCheck hiding (Fixed(..), (.&.))
+
+import Data.PrimitiveArray.Index.Set
+
+
+
+-- TODO what exactly does the mask fix? Only bits already @1@, or every bit
+-- as it is? The mask should actually freeze-fix those bits, where we are
+-- set to @1@!
+
+prop_Fixed_BitSet_setSucc (u :: Word, Fixed m s :: Fixed BitSet) = traceShow (tgo, tsu) $ tgo == tsu
+  where tgo = go s
+        tsu = (getFixed <$> setSucc (Fixed 0 0) (Fixed 0 h) (Fixed m s))
+        fb1 = m .&. s -- fixed bits to 1
+        fb0 = m .&. complement s  -- fixed bits to 0
+        h   = bit (fromIntegral $ u `mod` 8) - 1
+        go x -- continue creating successors, until the mask criterion is met (again).
+          | Nothing <- ssx = Nothing
+          | Just x' <- ssx
+          , fb0 == m .&. complement x'
+          , fb1 == m .&. x' = traceShow ('j',fb0,fb1,m,x,x') $ Just x'
+          | Just x' <- ssx  = traceShow ('g',fb0,fb1,m,x,x') $ go x'
+          where ssx = setSucc 0 h x
+
diff --git a/Data/PrimitiveArray/Zero.hs b/Data/PrimitiveArray/Zero.hs
deleted file mode 100644
--- a/Data/PrimitiveArray/Zero.hs
+++ /dev/null
@@ -1,147 +0,0 @@
-{-# LANGUAGE ScopedTypeVariables #-}
-{-# LANGUAGE TemplateHaskell #-}
-{-# LANGUAGE StandaloneDeriving #-}
-{-# LANGUAGE DeriveDataTypeable #-}
-{-# LANGUAGE CPP #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE FlexibleInstances #-}
-{-# LANGUAGE TypeFamilies #-}
-
--- | Primitive arrays where the lower index is zero (or the equivalent of zero
--- for newtypes and enumerations).
-
-module Data.PrimitiveArray.Zero where
-
-import           Control.Exception (assert)
-import           Control.Monad
-import           Control.Monad.Primitive (PrimState)
-import           Data.Array.Repa.Index
-import           Data.Array.Repa.Shape
-import           Data.Vector.Generic as G hiding (forM_, length, zipWithM_, new)
-import           Data.Vector.Generic.Mutable as GM hiding (length)
-import qualified Data.Vector as V hiding (forM_, length, zipWithM_)
-import qualified Data.Vector.Unboxed as VU hiding (forM_, length, zipWithM_)
-import qualified Data.Vector.Unboxed.Mutable as VUM hiding (length)
-
-import Data.Array.Repa.ExtShape
-import Data.PrimitiveArray
-
-
-
--- * Unboxed, multidimensional arrays.
-
-data Unboxed sh e = Unboxed !sh !(VU.Vector e)
-  deriving (Read,Show,Eq)
-
-data instance MutArr m (Unboxed sh e) = MUnboxed !sh !(VU.MVector (PrimState m) e)
-
-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 (MUnboxed _ mba) = ma
-    zipWithM_ (\k x -> assert (length xs == size exUb) $ unsafeWrite mba k x) [0.. size exUb -1] xs
-    return ma
-  newM inLb inUb = let exUb = inUb `addDim` unitDim in
-    unless (inLb == zeroDim) (error "MArr0 lb/=zeroDim") >>
-    MUnboxed exUb `liftM` new (size exUb)
-  newWithM inLb inUb def = do
-    let exUb = inUb `addDim` unitDim
-    ma <- newM inLb inUb
-    let (MUnboxed _ mba) = ma
-    forM_ [0 .. size exUb -1] $ \k -> unsafeWrite mba k def
-    return ma
-  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 #-}
-  {-# INLINE newWithM #-}
-  {-# INLINE readM #-}
-  {-# INLINE writeM #-}
-
-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)
-  deriving (Read,Show,Eq)
-
-data instance MutArr m (Boxed sh e) = MBoxed !sh !(V.MVector (PrimState m) e)
-
-instance (Shape sh, ExtShape sh, VUM.Unbox elm) => MPrimArrayOps Boxed sh elm where
-  boundsM (MBoxed exUb _) = (zeroDim,exUb `subDim` unitDim)
-  fromListM inLb inUb xs = do
-    ma <- newM inLb inUb
-    let exUb = inUb `addDim` unitDim
-    let (MBoxed _ mba) = ma
-    zipWithM_ (\k x -> assert (length xs == size exUb) $ unsafeWrite mba k x) [0 .. size exUb - 1] xs -- [0.. toIndex exUb inUb] xs
-    return ma
-  newM inLb inUb = let exUb = inUb `addDim` unitDim in
-    unless (inLb == zeroDim) (error "MArr0 lb/=zeroDim") >>
-    MBoxed exUb `liftM` new (size exUb)
-  newWithM inLb inUb def = do
-    let exUb = inUb `addDim` unitDim
-    ma <- newM inLb inUb
-    let (MBoxed _ mba) = ma
-    forM_ [0 .. size exUb -1] $ \k -> unsafeWrite mba k def
-    return ma
-  readM (MBoxed exUb mba) idx = assert (inShape exUb idx) $ unsafeRead mba (toIndex exUb idx)
-  writeM (MBoxed 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 Boxed sh elm where
-  bounds (Boxed exUb _) = (zeroDim,exUb `subDim` unitDim)
-  freeze (MBoxed exUb mba) = Boxed exUb `liftM` unsafeFreeze mba
-  index (Boxed exUb ba) idx = assert (inShape exUb idx) $ unsafeIndex ba (toIndex exUb idx)
-  {-# INLINE bounds #-}
-  {-# INLINE freeze #-}
-  {-# INLINE index #-}
-
-instance (Shape sh, ExtShape sh) => PrimArrayMap Boxed sh e e' where
-  map f (Boxed sh xs) = Boxed sh (V.map f xs)
-  {-# INLINE map #-}
-
-
-{-
- -
- - 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/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -1,4 +1,4 @@
-Copyright Christian Hoener zu Siederdissen 2010-2013
+Copyright Christian Hoener zu Siederdissen 2010-2015
 
 All rights reserved.
 
diff --git a/PrimitiveArray.cabal b/PrimitiveArray.cabal
--- a/PrimitiveArray.cabal
+++ b/PrimitiveArray.cabal
@@ -1,49 +1,109 @@
 Name:           PrimitiveArray
-Version:        0.5.4.0
+Version:        0.6.0.0
 License:        BSD3
 License-file:   LICENSE
-Author:         Christian Hoener zu Siederdissen
-Maintainer:     choener@tbi.univie.ac.at
-Copyright:      Christian Hoener zu Siederdissen, 2010-2014
-Homepage:       http://www.tbi.univie.ac.at/~choener/
+Maintainer:     choener@bioinf.uni-leipzig.de
+author:         Christian Hoener zu Siederdissen, 2011-2015
+copyright:      Christian Hoener zu Siederdissen, 2011-2015
+homepage:       http://www.bioinf.uni-leipzig.de/Software/gADP/
 Stability:      Experimental
 Category:       Data
 Build-type:     Simple
-Cabal-version:  >=1.6
-Synopsis:
-                Efficient multidimensional arrays
+Cabal-version:  >=1.10.0
+tested-with:    GHC == 7.8.4, GHC == 7.10.1
+Synopsis:       Efficient multidimensional arrays
 Description:
-                This library provides efficient multidimensional arrays.
+                This library provides efficient multidimensional arrays. Import
+                @Data.PrimitiveArray@ for indices, lenses, and arrays.
                 .
+                For ADPfusion users, the library also provides the machinary to
+                fill tables in the correct order required by usual CYK-style
+                parsers, or regular grammars (used e.g. in alignment
+                algorithms). This means that unless your grammar require a
+                strange order in which parsing is to be performed, it will
+                mostly "just work".
+                .
                 In general all operations are (highly) unsafe, no
                 bounds-checking or other sanity-checking is performed.
                 Operations are aimed toward efficiency as much as possible.
 
+
+
 extra-source-files:
-  changelog
+  README.md
+  changelog.md
 
+
+
 Library
   Exposed-modules:
-    Data.Array.Repa.ExtShape
-    Data.Array.Repa.Index.Outside
-    Data.Array.Repa.Index.Point
-    Data.Array.Repa.Index.Points
-    Data.Array.Repa.Index.Subword
     Data.PrimitiveArray
+    Data.PrimitiveArray.Class
+    Data.PrimitiveArray.Dense
     Data.PrimitiveArray.FillTables
-    Data.PrimitiveArray.QuickCheck
-    Data.PrimitiveArray.Zero
-  Build-depends:
-    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
+    Data.PrimitiveArray.Index
+    Data.PrimitiveArray.Index.Class
+    Data.PrimitiveArray.Index.Complement
+    Data.PrimitiveArray.Index.Int
+    Data.PrimitiveArray.Index.Outside
+    Data.PrimitiveArray.Index.Point
+    Data.PrimitiveArray.Index.Set
+    Data.PrimitiveArray.Index.Subword
+    Data.PrimitiveArray.QuickCheck.Index.Set
+  build-depends: base                     >= 4.7      && < 4.9
+               , aeson                    == 0.8.*
+               , binary                   == 0.7.*
+               , bits                     == 0.4.*
+               , cereal                   == 0.4.*
+               , deepseq                  >= 1.3      && < 1.5
+               , OrderedBits              == 0.0.0.*
+               , primitive                >= 0.5.4    && < 0.7
+               , QuickCheck               >= 2.7      && < 2.9
+               , vector                   == 0.10.*
+               , vector-binary-instances  == 0.2.*
+               , vector-th-unbox          == 0.2.*
+  default-extensions: BangPatterns
+                    , DefaultSignatures
+                    , DeriveGeneric
+                    , FlexibleContexts
+                    , FlexibleInstances
+                    , GeneralizedNewtypeDeriving
+                    , MultiParamTypeClasses
+                    , RankNTypes
+                    , ScopedTypeVariables
+                    , StandaloneDeriving
+                    , TemplateHaskell
+                    , TypeFamilies
+                    , TypeOperators
+                    , UndecidableInstances
+  default-language:
+    Haskell2010
   ghc-options:
     -O2
     -funbox-strict-fields
+
+
+
+test-suite properties
+  type:
+    exitcode-stdio-1.0
+  main-is:
+    properties.hs
+  ghc-options:
+    -threaded -rtsopts -with-rtsopts=-N
+  hs-source-dirs:
+    tests
+  default-language:
+    Haskell2010
+  default-extensions: TemplateHaskell
+  build-depends: base
+               , PrimitiveArray
+               , QuickCheck
+               , test-framework               >= 0.8  && < 0.9
+               , test-framework-quickcheck2   >= 0.3  && < 0.4
+               , test-framework-th            >= 0.2  && < 0.3
+
+
 
 source-repository head
   type: git
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,18 @@
+# PrimitiveArray
+
+[![Build Status](https://travis-ci.org/choener/PrimitiveArray.svg?branch=master)](https://travis-ci.org/choener/PrimitiveArray)
+
+PrimitiveArray provides operations on multi-dimensional arrays. Internally, the
+representation is based on the vector library, while the multi-dimensional
+indexing follows repa.
+
+Primitive arrays are designed to be used together with ADPfusion.
+
+
+
+#### Contact
+
+Christian Hoener zu Siederdissen
+choener@bioinf.uni-leipzig.de
+http://www.bioinf.uni-leipzig.de/~choener/
+
diff --git a/changelog b/changelog
deleted file mode 100644
--- a/changelog
+++ /dev/null
@@ -1,9 +0,0 @@
-0.5.4.0
-
-- actually implemented PointR
-
-- added the rather important strictness annotation for mutable arrays in .Zero
-
-0.5.3.0
-
-- fixed vector-th-unbox problem
diff --git a/changelog.md b/changelog.md
new file mode 100644
--- /dev/null
+++ b/changelog.md
@@ -0,0 +1,37 @@
+0.6.0.0
+-------
+
+- moved primitive array classes to Data.PrimitiveArray.Class
+- added _from / _to lenses
+- Field1 .. Field6 lenses for indices (Z:.a:.b...) (with Z being Field0)
+  - lens stuff currently commented out; aiming to have an extra package [lens
+    is fairly heavy]
+- FillTables should work now (with PointL, Subword)
+- freezing of whole stacks of (Z:.mutarr:.mutarr:. ...) tables
+- explicit 'Shape Subword'; this allows for simpler code in a number of places
+  and is especially useful for CYK-style algorithms that have a
+  single-dimensional upper-triangular matrix.
+- rangeStream of Extshape is new and used by the FillTables module
+- Binary, Cereal, Aeson instances for indices and immutable tables
+- orphan instances of Binary, Cereal, Aeson for Z, and (:.)
+- topmostIndex returns the final index position for CYK-style (bottom to top)
+  parsing
+- removed Data.Array.Repa.Index.Point (we have PointL, PointR in Points.hs)
+- added   Data.Array.Repa.Index.Set (for sets with an interface, used by
+  Hamiltonian path problems)
+- Data.Array.Repa.Index.Outside is now just a newtype wrapped around other
+  Index types. We want to be able to say "a Subword, but for Outside
+  algorithms"
+- travis-ci integration
+
+0.5.4.0
+-------
+
+- actually implemented PointR
+
+- added the rather important strictness annotation for mutable arrays in .Zero
+
+0.5.3.0
+-------
+
+- fixed vector-th-unbox problem
diff --git a/tests/properties.hs b/tests/properties.hs
new file mode 100644
--- /dev/null
+++ b/tests/properties.hs
@@ -0,0 +1,15 @@
+
+module Main where
+
+import           Test.Framework.Providers.QuickCheck2
+import           Test.Framework.TH
+
+-- import qualified Data.Bits.Ordered.QuickCheck as QC
+
+
+
+-- prop_PopCountSet = QC.prop_PopCountSet
+
+main :: IO ()
+main = $(defaultMainGenerator)
+