diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs
index 3fa5fcd..7bd4867 100644
--- a/brush-strokes/src/cusps/bench/Main.hs
+++ b/brush-strokes/src/cusps/bench/Main.hs
@@ -199,8 +199,8 @@ data TestCase =
   TestCase
   { testName        :: String
   , testBrushStroke :: BrushStroke
-  , testCuspOptions :: RootIsolationOptions
-  , testStartBoxes  :: [ ( Int, Box ) ]
+  , testCuspOptions :: RootIsolationOptions 2 3
+  , testStartBoxes  :: [ ( Int, Box 2 ) ]
   }
 
 brushStrokeFunctions
@@ -224,7 +224,7 @@ mkVal t i s = ( ℝ1 t, i, ℝ1 s )
 mkI :: ( Double, Double ) -> 𝕀 Double
 mkI ( lo, hi ) = 𝕀 lo hi
 
-mkBox :: ( Double, Double ) -> ( Double, Double ) -> Box
+mkBox :: ( Double, Double ) -> ( Double, Double ) -> Box 2
 mkBox ( t_min, t_max ) ( s_min, s_max ) =
   ( 𝕀 ( ℝ2 t_min s_min ) ( ℝ2 t_max s_max ) )
 
@@ -488,7 +488,7 @@ showD float = showFFloat (Just 6) float ""
 
 --------------------------------------------------------------------------------
 
-ellipseTestCase :: RootIsolationOptions -> String -> ( Double, Double ) -> Double -> [ ( Int, Box ) ] -> TestCase
+ellipseTestCase :: RootIsolationOptions 2 3 -> String -> ( Double, Double ) -> Double -> [ ( Int, Box 2 ) ] -> TestCase
 ellipseTestCase opts str k0k1 rot startBoxes =
   TestCase
     { testName = "ellipse (" ++ str ++ ")"
@@ -607,9 +607,11 @@ getStrokeFunctions brush sp0 crv =
         brush @3 @𝕀 proxy# singleton )
 {-# INLINEABLE getStrokeFunctions #-}
 
-defaultStartBoxes :: [ Int ] -> [ ( Int, Box ) ]
+defaultStartBoxes :: [ Int ] -> [ ( Int, Box 2 ) ]
 defaultStartBoxes is =
-  [ ( i, mkBox (zero, one) (zero, one) ) | i <- is ]
+  [ ( i, 𝕀 ( ℝ2 zero zero ) ( ℝ2 one one ) )
+  | i <- is
+  ]
 
 getR1 (ℝ1 u) = u
 
diff --git a/brush-strokes/src/lib/Math/Algebra/Dual.hs b/brush-strokes/src/lib/Math/Algebra/Dual.hs
index 4f73d74..d32f2cd 100644
--- a/brush-strokes/src/lib/Math/Algebra/Dual.hs
+++ b/brush-strokes/src/lib/Math/Algebra/Dual.hs
@@ -1061,14 +1061,14 @@ instance HasChainRule Double 2 ( ℝ 1 ) where
   konst   :: forall w. AbelianGroup w => w -> D2𝔸1 w
   konst w =
     let !o = fromInteger @w 0
-    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+    in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
 
-  value df = $$( monIndex [|| df ||] zeroMonomial )
+  value df = $$( monIndexQ [|| df ||] zeroMonomial )
 
   linearD :: forall w. Module Double ( T w ) => ( ℝ 1 -> w ) -> ℝ 1 -> D2𝔸1 w
   linearD f v =
     let !o = origin @Double @( T w )
-    in $$( monTabulate \ mon ->
+    in $$( monTabulateQ \ mon ->
               if | isZeroMonomial mon
                  -> [|| f v ||]
                  | Just i <- isLinear mon
@@ -1100,14 +1100,14 @@ instance HasChainRule Double 3 ( ℝ 1 ) where
   konst   :: forall w. AbelianGroup w => w -> D3𝔸1 w
   konst w =
     let !o = fromInteger @w 0
-    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+    in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
 
-  value df = $$( monIndex [|| df ||] zeroMonomial )
+  value df = $$( monIndexQ [|| df ||] zeroMonomial )
 
   linearD :: forall w. Module Double ( T w ) => ( ℝ 1 -> w ) -> ℝ 1 -> D3𝔸1 w
   linearD f v =
     let !o = origin @Double @( T w )
-    in $$( monTabulate \ mon ->
+    in $$( monTabulateQ \ mon ->
               if | isZeroMonomial mon
                  -> [|| f v ||]
                  | Just i <- isLinear mon
@@ -1139,14 +1139,14 @@ instance HasChainRule Double 2 ( ℝ 2 ) where
   konst   :: forall w. AbelianGroup w => w -> D2𝔸2 w
   konst w =
     let !o = fromInteger @w 0
-    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+    in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
 
-  value df = $$( monIndex [|| df ||] zeroMonomial )
+  value df = $$( monIndexQ [|| df ||] zeroMonomial )
 
   linearD :: forall w. Module Double ( T w ) => ( ℝ 2 -> w ) -> ℝ 2 -> D2𝔸2 w
   linearD f v =
     let !o = origin @Double @( T w )
-    in $$( monTabulate \ mon ->
+    in $$( monTabulateQ \ mon ->
               if | isZeroMonomial mon
                  -> [|| f v ||]
                  | Just i <- isLinear mon
@@ -1178,14 +1178,14 @@ instance HasChainRule Double 3 ( ℝ 2 ) where
   konst   :: forall w. AbelianGroup w => w -> D3𝔸2 w
   konst w =
     let !o = fromInteger @w 0
-    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+    in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
 
-  value df = $$( monIndex [|| df ||] zeroMonomial )
+  value df = $$( monIndexQ [|| df ||] zeroMonomial )
 
   linearD :: forall w. Module Double ( T w ) => ( ℝ 2 -> w ) -> ℝ 2 -> D3𝔸2 w
   linearD f v =
     let !o = origin @Double @( T w )
-    in $$( monTabulate \ mon ->
+    in $$( monTabulateQ \ mon ->
               if | isZeroMonomial mon
                  -> [|| f v ||]
                  | Just i <- isLinear mon
@@ -1217,14 +1217,14 @@ instance HasChainRule Double 2 ( ℝ 3 ) where
   konst   :: forall w. AbelianGroup w => w -> D2𝔸3 w
   konst w =
     let !o = fromInteger @w 0
-    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+    in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
 
-  value df = $$( monIndex [|| df ||] zeroMonomial )
+  value df = $$( monIndexQ [|| df ||] zeroMonomial )
 
   linearD :: forall w. Module Double ( T w ) => ( ℝ 3 -> w ) -> ℝ 3 -> D2𝔸3 w
   linearD f v =
     let !o = origin @Double @( T w )
-    in $$( monTabulate \ mon ->
+    in $$( monTabulateQ \ mon ->
               if | isZeroMonomial mon
                  -> [|| f v ||]
                  | Just i <- isLinear mon
@@ -1256,14 +1256,14 @@ instance HasChainRule Double 3 ( ℝ 3 ) where
   konst   :: forall w. AbelianGroup w => w -> D3𝔸3 w
   konst w =
     let !o = fromInteger @w 0
-    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+    in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
 
-  value df = $$( monIndex [|| df ||] zeroMonomial )
+  value df = $$( monIndexQ [|| df ||] zeroMonomial )
 
   linearD :: forall w. Module Double ( T w ) => ( ℝ 3 -> w ) -> ℝ 3 -> D3𝔸3 w
   linearD f v =
     let !o = origin @Double @( T w )
-    in $$( monTabulate \ mon ->
+    in $$( monTabulateQ \ mon ->
               if | isZeroMonomial mon
                  -> [|| f v ||]
                  | Just i <- isLinear mon
@@ -1295,14 +1295,14 @@ instance HasChainRule Double 2 ( ℝ 4 ) where
   konst   :: forall w. AbelianGroup w => w -> D2𝔸4 w
   konst w =
     let !o = fromInteger @w 0
-    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+    in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
 
-  value df = $$( monIndex [|| df ||] zeroMonomial )
+  value df = $$( monIndexQ [|| df ||] zeroMonomial )
 
   linearD :: forall w. Module Double ( T w ) => ( ℝ 4 -> w ) -> ℝ 4 -> D2𝔸4 w
   linearD f v =
     let !o = origin @Double @( T w )
-    in $$( monTabulate \ mon ->
+    in $$( monTabulateQ \ mon ->
               if | isZeroMonomial mon
                  -> [|| f v ||]
                  | Just i <- isLinear mon
@@ -1334,14 +1334,14 @@ instance HasChainRule Double 3 ( ℝ 4 ) where
   konst   :: forall w. AbelianGroup w => w -> D3𝔸4 w
   konst w =
     let !o = fromInteger @w 0
-    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+    in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
 
-  value df = $$( monIndex [|| df ||] zeroMonomial )
+  value df = $$( monIndexQ [|| df ||] zeroMonomial )
 
   linearD :: forall w. Module Double ( T w ) => ( ℝ 4 -> w ) -> ℝ 4 -> D3𝔸4 w
   linearD f v =
     let !o = origin @Double @( T w )
-    in $$( monTabulate \ mon ->
+    in $$( monTabulateQ \ mon ->
               if | isZeroMonomial mon
                  -> [|| f v ||]
                  | Just i <- isLinear mon
diff --git a/brush-strokes/src/lib/Math/Algebra/Dual/Internal.hs b/brush-strokes/src/lib/Math/Algebra/Dual/Internal.hs
index 2d9c446..787c03f 100644
--- a/brush-strokes/src/lib/Math/Algebra/Dual/Internal.hs
+++ b/brush-strokes/src/lib/Math/Algebra/Dual/Internal.hs
@@ -38,7 +38,8 @@ import Math.Linear
   , zipIndices
   )
 import Math.Monomial
-  ( Mon(..), MonomialBasis(..), Vars, Deg
+  ( Mon(..), MonomialBasisQ(..), MonomialBasis(..)
+  , Vars, Deg
   , mons, zeroMonomial, isZeroMonomial, totalDegree
   , multiSubsetSumFaà, multiSubsetsSum
   , partitionFaà, partitions
@@ -203,8 +204,8 @@ data D3𝔸4 v =
 -- (To be spliced in using Template Haskell.)
 chainRule1NQ :: forall dr1 dv v r w d
              .  ( Ring r, RepresentableQ r v
-                , MonomialBasis dr1, Vars dr1 ~ 1
-                , MonomialBasis dv , Vars dv  ~ RepDim v
+                , MonomialBasisQ dr1, Vars dr1 ~ 1
+                , MonomialBasisQ dv , Vars dv  ~ RepDim v
                 , Deg dr1 ~ Deg dv
                 , d ~ Vars dv, KnownNat d
                 )
@@ -215,19 +216,19 @@ chainRule1NQ :: forall dr1 dv v r w d
              -> CodeQ ( dv  w )
              -> CodeQ ( dr1 w )
 chainRule1NQ zero_w sum_w scale_w df dg =
-    monTabulate @dr1 \ ( Mon ( k `VS` _ ) ) ->
-      case k of
+    monTabulateQ @dr1 \ mon ->
+      case totalDegree mon of
         -- Set the value of the composition separately,
         -- as that isn't handled by the Faà di Bruno formula.
-        0 -> monIndex @dv dg zeroMonomial
-        _ ->
+        0 -> monIndexQ @dv dg zeroMonomial
+        k ->
           [|| unT $ $$( foldQ sum_w zero_w
-                  [ [|| $$scale_w ( T $$( monIndex @dv dg m_g ) )
+                  [ [|| $$scale_w ( T $$( monIndexQ @dv dg m_g ) )
                           $$( foldQ [|| (Ring.+) ||] [|| Ring.fromInteger ( 0 :: Integer ) ||]
                                 [ [|| Ring.fromInteger $$( liftTyped $ fromIntegral $ multiSubsetSumFaà k is ) Ring.*
                                      $$( foldQ [|| (Ring.*) ||] [|| Ring.fromInteger ( 1 :: Integer ) ||]
                                         [ foldQ [|| (Ring.*) ||] [|| Ring.fromInteger ( 1 :: Integer ) ||]
-                                          [ ( indexQ @r @v ( monIndex @dr1 df ( Mon $ f_deg `VS` VZ ) ) v_index )
+                                          [ ( indexQ @r @v ( monIndexQ @dr1 df ( Mon $ Vec [ f_deg ] ) ) v_index )
                                                 `powQ` pow
                                           | ( f_deg, pow ) <- pows
                                           ]
@@ -248,8 +249,8 @@ chainRule1NQ zero_w sum_w scale_w df dg =
 --
 -- (To be spliced in using Template Haskell.)
 chainRuleN1Q :: forall du dr1 r
-             .  ( MonomialBasis du
-                , MonomialBasis dr1, Vars dr1 ~ 1
+             .  ( MonomialBasisQ du
+                , MonomialBasisQ dr1, Vars dr1 ~ 1
                 )
              => CodeQ ( Integer -> r )   -- ^ fromInteger (circumvent TH constraint issue)
              -> CodeQ ( r -> r -> r )    -- ^ (+) (circumvent TH constraint issue)
@@ -259,19 +260,19 @@ chainRuleN1Q :: forall du dr1 r
              -> CodeQ ( dr1 r )
              -> CodeQ ( du  r )
 chainRuleN1Q fromInt add times pow df dg =
-    monTabulate @du \ mon ->
+    monTabulateQ @du \ mon ->
       if
         | isZeroMonomial mon
         -- Set the value of the composition separately,
         -- as that isn't handled by the Faà di Bruno formula.
-        -> monIndex @dr1 dg zeroMonomial
+        -> monIndexQ @dr1 dg zeroMonomial
         | otherwise
         -> foldQ add [|| $$fromInt ( 0 :: Integer ) ||]
-            [ [|| $$times $$( monIndex @dr1 dg ( Mon ( k `VS` VZ ) ) )
+            [ [|| $$times $$( monIndexQ @dr1 dg ( Mon ( Vec [ k ] ) ) )
                   $$( foldQ add [|| $$fromInt ( 0 :: Integer ) ||]
                       [ [|| $$times ( $$fromInt $$( liftTyped $ fromIntegral $ partitionFaà mon is ) )
                             $$( foldQ times [|| $$fromInt ( 1 :: Integer ) ||]
-                               [ [|| $$pow $$( monIndex @du df f_mon ) p ||]
+                               [ [|| $$pow $$( monIndexQ @du df f_mon ) p ||]
                                | ( f_mon, p ) <- is
                                ]
                              ) ||]
@@ -284,366 +285,395 @@ chainRuleN1Q fromInt add times pow df dg =
             ]
 
 --------------------------------------------------------------------------------
--- MonomialBasis instances follow (nothing else).
+-- MonomialBasisQ instances follow (nothing else).
 
 type instance Deg  D𝔸0 = 0
 type instance Vars D𝔸0 = 0
-instance MonomialBasis D𝔸0 where
-  monTabulate f =
+instance MonomialBasisQ D𝔸0 where
+  monTabulateQ f =
     [|| let
-          _D0_v = $$( f $ Mon VZ )
+          _D0_v = $$( f $ Mon ( Vec [ ] ) )
         in D0 { .. }
       ||]
 
-  monIndex d _ = [|| _D0_v $$d ||]
+  monIndexQ d _ = [|| _D0_v $$d ||]
 
 type instance Deg  D1𝔸1 = 1
 type instance Vars D1𝔸1 = 1
-instance MonomialBasis D1𝔸1 where
-  monTabulate f =
+instance MonomialBasisQ D1𝔸1 where
+  monTabulateQ f =
     [|| let
-          _D11_v    =   $$( f $ Mon ( 0 `VS` VZ ) )
-          _D11_dx   = T $$( f $ Mon ( 1 `VS` VZ ) )
+          _D11_v    =   $$( f $ Mon ( Vec [ 0 ] ) )
+          _D11_dx   = T $$( f $ Mon ( Vec [ 1 ] ) )
         in D11 { .. }
       ||]
 
-  monIndex d = \ case
-    Mon ( 1 `VS` VZ ) -> [|| unT $ _D11_dx   $$d ||]
+  monIndexQ d = \ case
+    Mon ( Vec [ 1 ] ) -> [|| unT $ _D11_dx   $$d ||]
     _                 -> [||       _D11_v    $$d ||]
 
 type instance Deg  D2𝔸1 = 2
 type instance Vars D2𝔸1 = 1
-instance MonomialBasis D2𝔸1 where
-  monTabulate f =
+instance MonomialBasisQ D2𝔸1 where
+  monTabulateQ f =
     [|| let
-          _D21_v    =   $$( f $ Mon ( 0 `VS` VZ ) )
-          _D21_dx   = T $$( f $ Mon ( 1 `VS` VZ ) )
-          _D21_dxdx = T $$( f $ Mon ( 2 `VS` VZ ) )
+          _D21_v    =   $$( f $ Mon ( Vec [ 0 ] ) )
+          _D21_dx   = T $$( f $ Mon ( Vec [ 1 ] ) )
+          _D21_dxdx = T $$( f $ Mon ( Vec [ 2 ] ) )
         in D21 { .. }
       ||]
 
-  monIndex d = \ case
-    Mon ( 1 `VS` VZ ) -> [|| unT $ _D21_dx   $$d ||]
-    Mon ( 2 `VS` VZ ) -> [|| unT $ _D21_dxdx $$d ||]
+  monIndexQ d = \ case
+    Mon ( Vec [ 1 ] ) -> [|| unT $ _D21_dx   $$d ||]
+    Mon ( Vec [ 2 ] ) -> [|| unT $ _D21_dxdx $$d ||]
     _                 -> [||       _D21_v    $$d ||]
 
 type instance Deg  D3𝔸1 = 3
 type instance Vars D3𝔸1 = 1
-instance MonomialBasis D3𝔸1 where
-  monTabulate f =
+instance MonomialBasisQ D3𝔸1 where
+  monTabulateQ f =
     [|| let
-          _D31_v      =   $$( f $ Mon ( 0 `VS` VZ ) )
-          _D31_dx     = T $$( f $ Mon ( 1 `VS` VZ ) )
-          _D31_dxdx   = T $$( f $ Mon ( 2 `VS` VZ ) )
-          _D31_dxdxdx = T $$( f $ Mon ( 3 `VS` VZ ) )
+          _D31_v      =   $$( f $ Mon ( Vec [ 0 ] ) )
+          _D31_dx     = T $$( f $ Mon ( Vec [ 1 ] ) )
+          _D31_dxdx   = T $$( f $ Mon ( Vec [ 2 ] ) )
+          _D31_dxdxdx = T $$( f $ Mon ( Vec [ 3 ] ) )
         in D31 { .. }
       ||]
 
-  monIndex d = \ case
-    Mon ( 1 `VS` VZ ) -> [|| unT $ _D31_dx     $$d ||]
-    Mon ( 2 `VS` VZ ) -> [|| unT $ _D31_dxdx   $$d ||]
-    Mon ( 3 `VS` VZ ) -> [|| unT $ _D31_dxdxdx $$d ||]
+  monIndexQ d = \ case
+    Mon ( Vec [ 1 ] ) -> [|| unT $ _D31_dx     $$d ||]
+    Mon ( Vec [ 2 ] ) -> [|| unT $ _D31_dxdx   $$d ||]
+    Mon ( Vec [ 3 ] ) -> [|| unT $ _D31_dxdxdx $$d ||]
     _                 -> [||       _D31_v      $$d ||]
 
 type instance Deg  D1𝔸2 = 1
 type instance Vars D1𝔸2 = 2
-instance MonomialBasis D1𝔸2 where
-  monTabulate f =
+instance MonomialBasisQ D1𝔸2 where
+  monTabulateQ f =
     [|| let
-          _D12_v    =   $$( f $ Mon ( 0 `VS` 0 `VS` VZ ) )
-          _D12_dx   = T $$( f $ Mon ( 1 `VS` 0 `VS` VZ ) )
-          _D12_dy   = T $$( f $ Mon ( 0 `VS` 1 `VS` VZ ) )
+          _D12_v    =   $$( f $ Mon ( Vec [ 0, 0 ] ) )
+          _D12_dx   = T $$( f $ Mon ( Vec [ 1, 0 ] ) )
+          _D12_dy   = T $$( f $ Mon ( Vec [ 0, 1 ] ) )
         in D12 { .. }
       ||]
 
+  monIndexQ d = \ case
+    Mon ( Vec [ 1, 0 ] ) -> [|| unT $ _D12_dx $$d ||]
+    Mon ( Vec [ 0, 1 ] ) -> [|| unT $ _D12_dy $$d ||]
+    _                    -> [||       _D12_v  $$d ||]
+
+
+instance MonomialBasis D1𝔸2 where
+  monTabulate f =
+    let
+      _D12_v    =     f $ Mon ( Vec [ 0, 0 ] )
+      _D12_dx   = T $ f $ Mon ( Vec [ 1, 0 ] )
+      _D12_dy   = T $ f $ Mon ( Vec [ 0, 1 ] )
+    in D12 { .. }
+
   monIndex d = \ case
-    Mon ( 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D12_dx $$d ||]
-    Mon ( 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D12_dy $$d ||]
-    _                        -> [||       _D12_v  $$d ||]
+    Mon ( Vec [ 1, 0 ] ) -> unT $ _D12_dx d
+    Mon ( Vec [ 0, 1 ] ) -> unT $ _D12_dy d
+    _                    ->       _D12_v  d
 
 type instance Deg  D2𝔸2 = 2
 type instance Vars D2𝔸2 = 2
-instance MonomialBasis D2𝔸2 where
-  monTabulate f =
+instance MonomialBasisQ D2𝔸2 where
+  monTabulateQ f =
     [|| let
-          _D22_v    =   $$( f $ Mon ( 0 `VS` 0 `VS` VZ ) )
-          _D22_dx   = T $$( f $ Mon ( 1 `VS` 0 `VS` VZ ) )
-          _D22_dy   = T $$( f $ Mon ( 0 `VS` 1 `VS` VZ ) )
-          _D22_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` VZ ) )
-          _D22_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` VZ ) )
-          _D22_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` VZ ) )
+          _D22_v    =   $$( f $ Mon ( Vec [ 0, 0 ] ) )
+          _D22_dx   = T $$( f $ Mon ( Vec [ 1, 0 ] ) )
+          _D22_dy   = T $$( f $ Mon ( Vec [ 0, 1 ] ) )
+          _D22_dxdx = T $$( f $ Mon ( Vec [ 2, 0 ] ) )
+          _D22_dxdy = T $$( f $ Mon ( Vec [ 1, 1 ] ) )
+          _D22_dydy = T $$( f $ Mon ( Vec [ 0, 2 ] ) )
         in D22 { .. }
       ||]
 
-  monIndex d = \ case
-    Mon ( 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D22_dx   $$d ||]
-    Mon ( 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D22_dy   $$d ||]
-    Mon ( 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D22_dxdx $$d ||]
-    Mon ( 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D22_dxdy $$d ||]
-    Mon ( 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D22_dydy $$d ||]
-    _                        -> [||       _D22_v    $$d ||]
+  monIndexQ d = \ case
+    Mon ( Vec [ 1, 0 ] ) -> [|| unT $ _D22_dx   $$d ||]
+    Mon ( Vec [ 0, 1 ] ) -> [|| unT $ _D22_dy   $$d ||]
+    Mon ( Vec [ 2, 0 ] ) -> [|| unT $ _D22_dxdx $$d ||]
+    Mon ( Vec [ 1, 1 ] ) -> [|| unT $ _D22_dxdy $$d ||]
+    Mon ( Vec [ 0, 2 ] ) -> [|| unT $ _D22_dydy $$d ||]
+    _                    -> [||       _D22_v    $$d ||]
 
 type instance Deg  D3𝔸2 = 3
 type instance Vars D3𝔸2 = 2
-instance MonomialBasis D3𝔸2 where
-  monTabulate f =
+instance MonomialBasisQ D3𝔸2 where
+  monTabulateQ f =
     [|| let
-          _D32_v      =   $$( f $ Mon ( 0 `VS` 0 `VS` VZ ) )
-          _D32_dx     = T $$( f $ Mon ( 1 `VS` 0 `VS` VZ ) )
-          _D32_dy     = T $$( f $ Mon ( 0 `VS` 1 `VS` VZ ) )
-          _D32_dxdx   = T $$( f $ Mon ( 2 `VS` 0 `VS` VZ ) )
-          _D32_dxdy   = T $$( f $ Mon ( 1 `VS` 1 `VS` VZ ) )
-          _D32_dydy   = T $$( f $ Mon ( 0 `VS` 2 `VS` VZ ) )
-          _D32_dxdxdx = T $$( f $ Mon ( 3 `VS` 0 `VS` VZ ) )
-          _D32_dxdxdy = T $$( f $ Mon ( 2 `VS` 1 `VS` VZ ) )
-          _D32_dxdydy = T $$( f $ Mon ( 1 `VS` 2 `VS` VZ ) )
-          _D32_dydydy = T $$( f $ Mon ( 0 `VS` 3 `VS` VZ ) )
+          _D32_v      =   $$( f $ Mon ( Vec [ 0, 0 ] ) )
+          _D32_dx     = T $$( f $ Mon ( Vec [ 1, 0 ] ) )
+          _D32_dy     = T $$( f $ Mon ( Vec [ 0, 1 ] ) )
+          _D32_dxdx   = T $$( f $ Mon ( Vec [ 2, 0 ] ) )
+          _D32_dxdy   = T $$( f $ Mon ( Vec [ 1, 1 ] ) )
+          _D32_dydy   = T $$( f $ Mon ( Vec [ 0, 2 ] ) )
+          _D32_dxdxdx = T $$( f $ Mon ( Vec [ 3, 0 ] ) )
+          _D32_dxdxdy = T $$( f $ Mon ( Vec [ 2, 1 ] ) )
+          _D32_dxdydy = T $$( f $ Mon ( Vec [ 1, 2 ] ) )
+          _D32_dydydy = T $$( f $ Mon ( Vec [ 0, 3 ] ) )
        in D32 { .. } ||]
 
-  monIndex d = \ case
-    Mon ( 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D32_dx     $$d ||]
-    Mon ( 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D32_dy     $$d ||]
-    Mon ( 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D32_dxdx   $$d ||]
-    Mon ( 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D32_dxdy   $$d ||]
-    Mon ( 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D32_dydy   $$d ||]
-    Mon ( 3 `VS` 0 `VS` VZ ) -> [|| unT $ _D32_dxdxdx $$d ||]
-    Mon ( 2 `VS` 1 `VS` VZ ) -> [|| unT $ _D32_dxdxdy $$d ||]
-    Mon ( 1 `VS` 2 `VS` VZ ) -> [|| unT $ _D32_dxdydy $$d ||]
-    Mon ( 0 `VS` 3 `VS` VZ ) -> [|| unT $ _D32_dydydy $$d ||]
-    _                        -> [||       _D32_v      $$d ||]
+  monIndexQ d = \ case
+    Mon ( Vec [ 1, 0 ] ) -> [|| unT $ _D32_dx     $$d ||]
+    Mon ( Vec [ 0, 1 ] ) -> [|| unT $ _D32_dy     $$d ||]
+    Mon ( Vec [ 2, 0 ] ) -> [|| unT $ _D32_dxdx   $$d ||]
+    Mon ( Vec [ 1, 1 ] ) -> [|| unT $ _D32_dxdy   $$d ||]
+    Mon ( Vec [ 0, 2 ] ) -> [|| unT $ _D32_dydy   $$d ||]
+    Mon ( Vec [ 3, 0 ] ) -> [|| unT $ _D32_dxdxdx $$d ||]
+    Mon ( Vec [ 2, 1 ] ) -> [|| unT $ _D32_dxdxdy $$d ||]
+    Mon ( Vec [ 1, 2 ] ) -> [|| unT $ _D32_dxdydy $$d ||]
+    Mon ( Vec [ 0, 3 ] ) -> [|| unT $ _D32_dydydy $$d ||]
+    _                    -> [||       _D32_v      $$d ||]
 
 type instance Deg  D1𝔸3 = 1
 type instance Vars D1𝔸3 = 3
-instance MonomialBasis D1𝔸3 where
-  monTabulate f =
+instance MonomialBasisQ D1𝔸3 where
+  monTabulateQ f =
     [|| let
-          !_D13_v    =   $$( f ( Mon ( 0 `VS` 0 `VS` 0 `VS` VZ ) ) )
-          !_D13_dx   = T $$( f ( Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) ) )
-          !_D13_dy   = T $$( f ( Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) ) )
-          !_D13_dz   = T $$( f ( Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) ) )
+          !_D13_v    =   $$( f ( Mon ( Vec [ 0, 0, 0 ] ) ) )
+          !_D13_dx   = T $$( f ( Mon ( Vec [ 1, 0, 0 ] ) ) )
+          !_D13_dy   = T $$( f ( Mon ( Vec [ 0, 1, 0 ] ) ) )
+          !_D13_dz   = T $$( f ( Mon ( Vec [ 0, 0, 1 ] ) ) )
         in D13 { .. }
       ||]
 
+  monIndexQ d = \ case
+    Mon ( Vec [ 1, 0, 0 ] ) -> [|| unT $ _D13_dx $$d ||]
+    Mon ( Vec [ 0, 1, 0 ] ) -> [|| unT $ _D13_dy $$d ||]
+    Mon ( Vec [ 0, 0, 1 ] ) -> [|| unT $ _D13_dz $$d ||]
+    _                       -> [||       _D13_v  $$d ||]
+
+instance MonomialBasis D1𝔸3 where
+  monTabulate f =
+    let
+      !_D13_v    =     f ( Mon ( Vec [ 0, 0, 0 ] ) )
+      !_D13_dx   = T $ f ( Mon ( Vec [ 1, 0, 0 ] ) )
+      !_D13_dy   = T $ f ( Mon ( Vec [ 0, 1, 0 ] ) )
+      !_D13_dz   = T $ f ( Mon ( Vec [ 0, 0, 1 ] ) )
+    in D13 { .. }
+
   monIndex d = \ case
-    Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D13_dx $$d ||]
-    Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D13_dy $$d ||]
-    Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D13_dz $$d ||]
-    _                               -> [||       _D13_v  $$d ||]
+    Mon ( Vec [ 1, 0, 0 ] ) -> unT $ _D13_dx d
+    Mon ( Vec [ 0, 1, 0 ] ) -> unT $ _D13_dy d
+    Mon ( Vec [ 0, 0, 1 ] ) -> unT $ _D13_dz d
+    _                       ->       _D13_v  d
 
 type instance Deg  D2𝔸3 = 2
 type instance Vars D2𝔸3 = 3
-instance MonomialBasis D2𝔸3 where
-  monTabulate f =
+instance MonomialBasisQ D2𝔸3 where
+  monTabulateQ f =
     [|| let
-          _D23_v    =   $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D23_dx   = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D23_dy   = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) )
-          _D23_dz   = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) )
-          _D23_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D23_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) )
-          _D23_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) )
-          _D23_dxdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) )
-          _D23_dydz = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) )
-          _D23_dzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) )
+          _D23_v    =   $$( f $ Mon ( Vec [ 0, 0, 0 ] ) )
+          _D23_dx   = T $$( f $ Mon ( Vec [ 1, 0, 0 ] ) )
+          _D23_dy   = T $$( f $ Mon ( Vec [ 0, 1, 0 ] ) )
+          _D23_dz   = T $$( f $ Mon ( Vec [ 0, 0, 1 ] ) )
+          _D23_dxdx = T $$( f $ Mon ( Vec [ 2, 0, 0 ] ) )
+          _D23_dxdy = T $$( f $ Mon ( Vec [ 1, 1, 0 ] ) )
+          _D23_dydy = T $$( f $ Mon ( Vec [ 0, 2, 0 ] ) )
+          _D23_dxdz = T $$( f $ Mon ( Vec [ 1, 0, 1 ] ) )
+          _D23_dydz = T $$( f $ Mon ( Vec [ 0, 1, 1 ] ) )
+          _D23_dzdz = T $$( f $ Mon ( Vec [ 0, 0, 2 ] ) )
         in D23 { .. }
       ||]
 
-  monIndex d = \ case
-    Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D23_dx   $$d ||]
-    Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D23_dy   $$d ||]
-    Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D23_dz   $$d ||]
-    Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D23_dxdx $$d ||]
-    Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D23_dxdy $$d ||]
-    Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D23_dydy $$d ||]
-    Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D23_dxdz $$d ||]
-    Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D23_dydz $$d ||]
-    Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D23_dzdz $$d ||]
-    _                               -> [||       _D23_v    $$d ||]
+  monIndexQ d = \ case
+    Mon ( Vec [ 1, 0, 0 ] ) -> [|| unT $ _D23_dx   $$d ||]
+    Mon ( Vec [ 0, 1, 0 ] ) -> [|| unT $ _D23_dy   $$d ||]
+    Mon ( Vec [ 0, 0, 1 ] ) -> [|| unT $ _D23_dz   $$d ||]
+    Mon ( Vec [ 2, 0, 0 ] ) -> [|| unT $ _D23_dxdx $$d ||]
+    Mon ( Vec [ 1, 1, 0 ] ) -> [|| unT $ _D23_dxdy $$d ||]
+    Mon ( Vec [ 0, 2, 0 ] ) -> [|| unT $ _D23_dydy $$d ||]
+    Mon ( Vec [ 1, 0, 1 ] ) -> [|| unT $ _D23_dxdz $$d ||]
+    Mon ( Vec [ 0, 1, 1 ] ) -> [|| unT $ _D23_dydz $$d ||]
+    Mon ( Vec [ 0, 0, 2 ] ) -> [|| unT $ _D23_dzdz $$d ||]
+    _                       -> [||       _D23_v    $$d ||]
 
 
 type instance Deg  D3𝔸3 = 3
 type instance Vars D3𝔸3 = 3
-instance MonomialBasis D3𝔸3 where
-  monTabulate f =
+instance MonomialBasisQ D3𝔸3 where
+  monTabulateQ f =
     [|| let
-          _D33_v      =   $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D33_dx     = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D33_dy     = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) )
-          _D33_dz     = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) )
-          _D33_dxdx   = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D33_dxdy   = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) )
-          _D33_dydy   = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) )
-          _D33_dxdz   = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) )
-          _D33_dydz   = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) )
-          _D33_dzdz   = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) )
-          _D33_dxdxdx = T $$( f $ Mon ( 3 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D33_dxdxdy = T $$( f $ Mon ( 2 `VS` 1 `VS` 0 `VS` VZ ) )
-          _D33_dxdydy = T $$( f $ Mon ( 1 `VS` 2 `VS` 0 `VS` VZ ) )
-          _D33_dydydy = T $$( f $ Mon ( 0 `VS` 3 `VS` 0 `VS` VZ ) )
-          _D33_dxdxdz = T $$( f $ Mon ( 2 `VS` 0 `VS` 1 `VS` VZ ) )
-          _D33_dxdydz = T $$( f $ Mon ( 1 `VS` 1 `VS` 1 `VS` VZ ) )
-          _D33_dxdzdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 2 `VS` VZ ) )
-          _D33_dydydz = T $$( f $ Mon ( 0 `VS` 2 `VS` 1 `VS` VZ ) )
-          _D33_dydzdz = T $$( f $ Mon ( 0 `VS` 1 `VS` 2 `VS` VZ ) )
-          _D33_dzdzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 3 `VS` VZ ) )
+          _D33_v      =   $$( f $ Mon ( Vec [ 0, 0, 0 ] ) )
+          _D33_dx     = T $$( f $ Mon ( Vec [ 1, 0, 0 ] ) )
+          _D33_dy     = T $$( f $ Mon ( Vec [ 0, 1, 0 ] ) )
+          _D33_dz     = T $$( f $ Mon ( Vec [ 0, 0, 1 ] ) )
+          _D33_dxdx   = T $$( f $ Mon ( Vec [ 2, 0, 0 ] ) )
+          _D33_dxdy   = T $$( f $ Mon ( Vec [ 1, 1, 0 ] ) )
+          _D33_dydy   = T $$( f $ Mon ( Vec [ 0, 2, 0 ] ) )
+          _D33_dxdz   = T $$( f $ Mon ( Vec [ 1, 0, 1 ] ) )
+          _D33_dydz   = T $$( f $ Mon ( Vec [ 0, 1, 1 ] ) )
+          _D33_dzdz   = T $$( f $ Mon ( Vec [ 0, 0, 2 ] ) )
+          _D33_dxdxdx = T $$( f $ Mon ( Vec [ 3, 0, 0 ] ) )
+          _D33_dxdxdy = T $$( f $ Mon ( Vec [ 2, 1, 0 ] ) )
+          _D33_dxdydy = T $$( f $ Mon ( Vec [ 1, 2, 0 ] ) )
+          _D33_dydydy = T $$( f $ Mon ( Vec [ 0, 3, 0 ] ) )
+          _D33_dxdxdz = T $$( f $ Mon ( Vec [ 2, 0, 1 ] ) )
+          _D33_dxdydz = T $$( f $ Mon ( Vec [ 1, 1, 1 ] ) )
+          _D33_dxdzdz = T $$( f $ Mon ( Vec [ 1, 0, 2 ] ) )
+          _D33_dydydz = T $$( f $ Mon ( Vec [ 0, 2, 1 ] ) )
+          _D33_dydzdz = T $$( f $ Mon ( Vec [ 0, 1, 2 ] ) )
+          _D33_dzdzdz = T $$( f $ Mon ( Vec [ 0, 0, 3 ] ) )
        in D33 { .. } ||]
 
-  monIndex d = \ case
-    Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dx     $$d ||]
-    Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dy     $$d ||]
-    Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dz     $$d ||]
-    Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dxdx   $$d ||]
-    Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dxdy   $$d ||]
-    Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dydy   $$d ||]
-    Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dxdz   $$d ||]
-    Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dydz   $$d ||]
-    Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D33_dzdz   $$d ||]
-    Mon ( 3 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dxdxdx $$d ||]
-    Mon ( 2 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dxdxdy $$d ||]
-    Mon ( 1 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dxdydy $$d ||]
-    Mon ( 0 `VS` 3 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dydydy $$d ||]
-    Mon ( 2 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dxdxdz $$d ||]
-    Mon ( 1 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dxdydz $$d ||]
-    Mon ( 1 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D33_dxdzdz $$d ||]
-    Mon ( 0 `VS` 2 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dydydz $$d ||]
-    Mon ( 0 `VS` 1 `VS` 2 `VS` VZ ) -> [|| unT $ _D33_dydzdz $$d ||]
-    Mon ( 0 `VS` 0 `VS` 3 `VS` VZ ) -> [|| unT $ _D33_dzdzdz $$d ||]
-    _                               -> [||       _D33_v      $$d ||]
+  monIndexQ d = \ case
+    Mon ( Vec [ 1, 0, 0 ] ) -> [|| unT $ _D33_dx     $$d ||]
+    Mon ( Vec [ 0, 1, 0 ] ) -> [|| unT $ _D33_dy     $$d ||]
+    Mon ( Vec [ 0, 0, 1 ] ) -> [|| unT $ _D33_dz     $$d ||]
+    Mon ( Vec [ 2, 0, 0 ] ) -> [|| unT $ _D33_dxdx   $$d ||]
+    Mon ( Vec [ 1, 1, 0 ] ) -> [|| unT $ _D33_dxdy   $$d ||]
+    Mon ( Vec [ 0, 2, 0 ] ) -> [|| unT $ _D33_dydy   $$d ||]
+    Mon ( Vec [ 1, 0, 1 ] ) -> [|| unT $ _D33_dxdz   $$d ||]
+    Mon ( Vec [ 0, 1, 1 ] ) -> [|| unT $ _D33_dydz   $$d ||]
+    Mon ( Vec [ 0, 0, 2 ] ) -> [|| unT $ _D33_dzdz   $$d ||]
+    Mon ( Vec [ 3, 0, 0 ] ) -> [|| unT $ _D33_dxdxdx $$d ||]
+    Mon ( Vec [ 2, 1, 0 ] ) -> [|| unT $ _D33_dxdxdy $$d ||]
+    Mon ( Vec [ 1, 2, 0 ] ) -> [|| unT $ _D33_dxdydy $$d ||]
+    Mon ( Vec [ 0, 3, 0 ] ) -> [|| unT $ _D33_dydydy $$d ||]
+    Mon ( Vec [ 2, 0, 1 ] ) -> [|| unT $ _D33_dxdxdz $$d ||]
+    Mon ( Vec [ 1, 1, 1 ] ) -> [|| unT $ _D33_dxdydz $$d ||]
+    Mon ( Vec [ 1, 0, 2 ] ) -> [|| unT $ _D33_dxdzdz $$d ||]
+    Mon ( Vec [ 0, 2, 1 ] ) -> [|| unT $ _D33_dydydz $$d ||]
+    Mon ( Vec [ 0, 1, 2 ] ) -> [|| unT $ _D33_dydzdz $$d ||]
+    Mon ( Vec [ 0, 0, 3 ] ) -> [|| unT $ _D33_dzdzdz $$d ||]
+    _                       -> [||       _D33_v      $$d ||]
 
 type instance Deg  D1𝔸4 = 1
 type instance Vars D1𝔸4 = 4
-instance MonomialBasis D1𝔸4 where
-  monTabulate f =
+instance MonomialBasisQ D1𝔸4 where
+  monTabulateQ f =
     [|| let
-          _D14_v    =   $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D14_dx   = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D14_dy   = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D14_dz   = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) )
-          _D14_dw   = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) )
+          _D14_v    =   $$( f $ Mon ( Vec [ 0, 0, 0, 0 ] ) )
+          _D14_dx   = T $$( f $ Mon ( Vec [ 1, 0, 0, 0 ] ) )
+          _D14_dy   = T $$( f $ Mon ( Vec [ 0, 1, 0, 0 ] ) )
+          _D14_dz   = T $$( f $ Mon ( Vec [ 0, 0, 1, 0 ] ) )
+          _D14_dw   = T $$( f $ Mon ( Vec [ 0, 0, 0, 1 ] ) )
         in D14 { .. } ||]
 
-  monIndex d = \ case
-    Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D14_dx   $$d ||]
-    Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D14_dy   $$d ||]
-    Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D14_dz   $$d ||]
-    Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D14_dw   $$d ||]
-    _                                      -> [||       _D14_v    $$d ||]
+  monIndexQ d = \ case
+    Mon ( Vec [ 1, 0, 0, 0 ] ) -> [|| unT $ _D14_dx   $$d ||]
+    Mon ( Vec [ 0, 1, 0, 0 ] ) -> [|| unT $ _D14_dy   $$d ||]
+    Mon ( Vec [ 0, 0, 1, 0 ] ) -> [|| unT $ _D14_dz   $$d ||]
+    Mon ( Vec [ 0, 0, 0, 1 ] ) -> [|| unT $ _D14_dw   $$d ||]
+    _                          -> [||       _D14_v    $$d ||]
 
 type instance Deg  D2𝔸4 = 2
 type instance Vars D2𝔸4 = 4
-instance MonomialBasis D2𝔸4 where
-  monTabulate f =
+instance MonomialBasisQ D2𝔸4 where
+  monTabulateQ f =
     [|| let
-          _D24_v    =   $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D24_dx   = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D24_dy   = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D24_dz   = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) )
-          _D24_dw   = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) )
-          _D24_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D24_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D24_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D24_dxdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) )
-          _D24_dydz = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) )
-          _D24_dzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) )
-          _D24_dxdw = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) )
-          _D24_dydw = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) )
-          _D24_dzdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) )
-          _D24_dwdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) )
+          _D24_v    =   $$( f $ Mon ( Vec [ 0, 0, 0, 0 ] ) )
+          _D24_dx   = T $$( f $ Mon ( Vec [ 1, 0, 0, 0 ] ) )
+          _D24_dy   = T $$( f $ Mon ( Vec [ 0, 1, 0, 0 ] ) )
+          _D24_dz   = T $$( f $ Mon ( Vec [ 0, 0, 1, 0 ] ) )
+          _D24_dw   = T $$( f $ Mon ( Vec [ 0, 0, 0, 1 ] ) )
+          _D24_dxdx = T $$( f $ Mon ( Vec [ 2, 0, 0, 0 ] ) )
+          _D24_dxdy = T $$( f $ Mon ( Vec [ 1, 1, 0, 0 ] ) )
+          _D24_dydy = T $$( f $ Mon ( Vec [ 0, 2, 0, 0 ] ) )
+          _D24_dxdz = T $$( f $ Mon ( Vec [ 1, 0, 1, 0 ] ) )
+          _D24_dydz = T $$( f $ Mon ( Vec [ 0, 1, 1, 0 ] ) )
+          _D24_dzdz = T $$( f $ Mon ( Vec [ 0, 0, 2, 0 ] ) )
+          _D24_dxdw = T $$( f $ Mon ( Vec [ 1, 0, 0, 1 ] ) )
+          _D24_dydw = T $$( f $ Mon ( Vec [ 0, 1, 0, 1 ] ) )
+          _D24_dzdw = T $$( f $ Mon ( Vec [ 0, 0, 1, 1 ] ) )
+          _D24_dwdw = T $$( f $ Mon ( Vec [ 0, 0, 0, 2 ] ) )
         in D24 { .. } ||]
 
-  monIndex d = \ case
-    Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dx   $$d ||]
-    Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dy   $$d ||]
-    Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dz   $$d ||]
-    Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D24_dw   $$d ||]
-    Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dxdx $$d ||]
-    Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dxdy $$d ||]
-    Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dydy $$d ||]
-    Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dxdz $$d ||]
-    Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dydz $$d ||]
-    Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dzdz $$d ||]
-    Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D24_dxdw $$d ||]
-    Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D24_dydw $$d ||]
-    Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D24_dzdw $$d ||]
-    Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D24_dwdw $$d ||]
-    _                                      -> [||       _D24_v    $$d ||]
+  monIndexQ d = \ case
+    Mon ( Vec [ 1, 0, 0, 0 ] ) -> [|| unT $ _D24_dx   $$d ||]
+    Mon ( Vec [ 0, 1, 0, 0 ] ) -> [|| unT $ _D24_dy   $$d ||]
+    Mon ( Vec [ 0, 0, 1, 0 ] ) -> [|| unT $ _D24_dz   $$d ||]
+    Mon ( Vec [ 0, 0, 0, 1 ] ) -> [|| unT $ _D24_dw   $$d ||]
+    Mon ( Vec [ 2, 0, 0, 0 ] ) -> [|| unT $ _D24_dxdx $$d ||]
+    Mon ( Vec [ 1, 1, 0, 0 ] ) -> [|| unT $ _D24_dxdy $$d ||]
+    Mon ( Vec [ 0, 2, 0, 0 ] ) -> [|| unT $ _D24_dydy $$d ||]
+    Mon ( Vec [ 1, 0, 1, 0 ] ) -> [|| unT $ _D24_dxdz $$d ||]
+    Mon ( Vec [ 0, 1, 1, 0 ] ) -> [|| unT $ _D24_dydz $$d ||]
+    Mon ( Vec [ 0, 0, 2, 0 ] ) -> [|| unT $ _D24_dzdz $$d ||]
+    Mon ( Vec [ 1, 0, 0, 1 ] ) -> [|| unT $ _D24_dxdw $$d ||]
+    Mon ( Vec [ 0, 1, 0, 1 ] ) -> [|| unT $ _D24_dydw $$d ||]
+    Mon ( Vec [ 0, 0, 1, 1 ] ) -> [|| unT $ _D24_dzdw $$d ||]
+    Mon ( Vec [ 0, 0, 0, 2 ] ) -> [|| unT $ _D24_dwdw $$d ||]
+    _                          -> [||       _D24_v    $$d ||]
 
 
 type instance Deg  D3𝔸4 = 3
 type instance Vars D3𝔸4 = 4
-instance MonomialBasis D3𝔸4 where
-  monTabulate f =
+instance MonomialBasisQ D3𝔸4 where
+  monTabulateQ f =
     [|| let
-          _D34_v      =   $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D34_dx     = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D34_dy     = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D34_dz     = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) )
-          _D34_dw     = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) )
-          _D34_dxdx   = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D34_dxdy   = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D34_dydy   = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D34_dxdz   = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) )
-          _D34_dydz   = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) )
-          _D34_dzdz   = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) )
-          _D34_dxdw   = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) )
-          _D34_dydw   = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) )
-          _D34_dzdw   = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) )
-          _D34_dwdw   = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) )
-          _D34_dxdxdx = T $$( f $ Mon ( 3 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D34_dxdxdy = T $$( f $ Mon ( 2 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D34_dxdydy = T $$( f $ Mon ( 1 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D34_dydydy = T $$( f $ Mon ( 0 `VS` 3 `VS` 0 `VS` 0 `VS` VZ ) )
-          _D34_dxdxdz = T $$( f $ Mon ( 2 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) )
-          _D34_dxdydz = T $$( f $ Mon ( 1 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) )
-          _D34_dxdzdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) )
-          _D34_dydydz = T $$( f $ Mon ( 0 `VS` 2 `VS` 1 `VS` 0 `VS` VZ ) )
-          _D34_dydzdz = T $$( f $ Mon ( 0 `VS` 1 `VS` 2 `VS` 0 `VS` VZ ) )
-          _D34_dzdzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 3 `VS` 0 `VS` VZ ) )
-          _D34_dxdxdw = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) )
-          _D34_dxdydw = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) )
-          _D34_dydydw = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 1 `VS` VZ ) )
-          _D34_dxdzdw = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) )
-          _D34_dydzdw = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 1 `VS` VZ ) )
-          _D34_dzdzdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 1 `VS` VZ ) )
-          _D34_dxdwdw = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) )
-          _D34_dydwdw = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 2 `VS` VZ ) )
-          _D34_dzdwdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 2 `VS` VZ ) )
-          _D34_dwdwdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 3 `VS` VZ ) )
+          _D34_v      =   $$( f $ Mon ( Vec [ 0, 0, 0, 0 ] ) )
+          _D34_dx     = T $$( f $ Mon ( Vec [ 1, 0, 0, 0 ] ) )
+          _D34_dy     = T $$( f $ Mon ( Vec [ 0, 1, 0, 0 ] ) )
+          _D34_dz     = T $$( f $ Mon ( Vec [ 0, 0, 1, 0 ] ) )
+          _D34_dw     = T $$( f $ Mon ( Vec [ 0, 0, 0, 1 ] ) )
+          _D34_dxdx   = T $$( f $ Mon ( Vec [ 2, 0, 0, 0 ] ) )
+          _D34_dxdy   = T $$( f $ Mon ( Vec [ 1, 1, 0, 0 ] ) )
+          _D34_dydy   = T $$( f $ Mon ( Vec [ 0, 2, 0, 0 ] ) )
+          _D34_dxdz   = T $$( f $ Mon ( Vec [ 1, 0, 1, 0 ] ) )
+          _D34_dydz   = T $$( f $ Mon ( Vec [ 0, 1, 1, 0 ] ) )
+          _D34_dzdz   = T $$( f $ Mon ( Vec [ 0, 0, 2, 0 ] ) )
+          _D34_dxdw   = T $$( f $ Mon ( Vec [ 1, 0, 0, 1 ] ) )
+          _D34_dydw   = T $$( f $ Mon ( Vec [ 0, 1, 0, 1 ] ) )
+          _D34_dzdw   = T $$( f $ Mon ( Vec [ 0, 0, 1, 1 ] ) )
+          _D34_dwdw   = T $$( f $ Mon ( Vec [ 0, 0, 0, 2 ] ) )
+          _D34_dxdxdx = T $$( f $ Mon ( Vec [ 3, 0, 0, 0 ] ) )
+          _D34_dxdxdy = T $$( f $ Mon ( Vec [ 2, 1, 0, 0 ] ) )
+          _D34_dxdydy = T $$( f $ Mon ( Vec [ 1, 2, 0, 0 ] ) )
+          _D34_dydydy = T $$( f $ Mon ( Vec [ 0, 3, 0, 0 ] ) )
+          _D34_dxdxdz = T $$( f $ Mon ( Vec [ 2, 0, 1, 0 ] ) )
+          _D34_dxdydz = T $$( f $ Mon ( Vec [ 1, 1, 1, 0 ] ) )
+          _D34_dxdzdz = T $$( f $ Mon ( Vec [ 1, 0, 2, 0 ] ) )
+          _D34_dydydz = T $$( f $ Mon ( Vec [ 0, 2, 1, 0 ] ) )
+          _D34_dydzdz = T $$( f $ Mon ( Vec [ 0, 1, 2, 0 ] ) )
+          _D34_dzdzdz = T $$( f $ Mon ( Vec [ 0, 0, 3, 0 ] ) )
+          _D34_dxdxdw = T $$( f $ Mon ( Vec [ 2, 0, 0, 1 ] ) )
+          _D34_dxdydw = T $$( f $ Mon ( Vec [ 1, 1, 0, 1 ] ) )
+          _D34_dydydw = T $$( f $ Mon ( Vec [ 0, 2, 0, 1 ] ) )
+          _D34_dxdzdw = T $$( f $ Mon ( Vec [ 1, 0, 1, 1 ] ) )
+          _D34_dydzdw = T $$( f $ Mon ( Vec [ 0, 1, 1, 1 ] ) )
+          _D34_dzdzdw = T $$( f $ Mon ( Vec [ 0, 0, 2, 1 ] ) )
+          _D34_dxdwdw = T $$( f $ Mon ( Vec [ 1, 0, 0, 2 ] ) )
+          _D34_dydwdw = T $$( f $ Mon ( Vec [ 0, 1, 0, 2 ] ) )
+          _D34_dzdwdw = T $$( f $ Mon ( Vec [ 0, 0, 1, 2 ] ) )
+          _D34_dwdwdw = T $$( f $ Mon ( Vec [ 0, 0, 0, 3 ] ) )
         in D34 { .. } ||]
 
-  monIndex d = \ case
-    Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dx     $$d ||]
-    Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dy     $$d ||]
-    Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dz     $$d ||]
-    Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dw     $$d ||]
-    Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdx   $$d ||]
-    Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdy   $$d ||]
-    Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dydy   $$d ||]
-    Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdz   $$d ||]
-    Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dydz   $$d ||]
-    Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dzdz   $$d ||]
-    Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dxdw   $$d ||]
-    Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dydw   $$d ||]
-    Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dzdw   $$d ||]
-    Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D34_dwdw   $$d ||]
-    Mon ( 3 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdxdx $$d ||]
-    Mon ( 2 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdxdy $$d ||]
-    Mon ( 1 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdydy $$d ||]
-    Mon ( 0 `VS` 3 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dydydy $$d ||]
-    Mon ( 2 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdxdz $$d ||]
-    Mon ( 1 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdydz $$d ||]
-    Mon ( 1 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdzdz $$d ||]
-    Mon ( 0 `VS` 2 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dydydz $$d ||]
-    Mon ( 0 `VS` 1 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dydzdz $$d ||]
-    Mon ( 0 `VS` 0 `VS` 3 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dzdzdz $$d ||]
-    Mon ( 2 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dxdxdw $$d ||]
-    Mon ( 1 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dxdydw $$d ||]
-    Mon ( 0 `VS` 2 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dydydw $$d ||]
-    Mon ( 1 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dxdzdw $$d ||]
-    Mon ( 0 `VS` 1 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dydzdw $$d ||]
-    Mon ( 0 `VS` 0 `VS` 2 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dzdzdw $$d ||]
-    Mon ( 1 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D34_dxdwdw $$d ||]
-    Mon ( 0 `VS` 1 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D34_dydwdw $$d ||]
-    Mon ( 0 `VS` 0 `VS` 1 `VS` 2 `VS` VZ ) -> [|| unT $ _D34_dzdwdw $$d ||]
-    Mon ( 0 `VS` 0 `VS` 0 `VS` 3 `VS` VZ ) -> [|| unT $ _D34_dwdwdw $$d ||]
-    _                                      -> [||       _D34_v      $$d ||]
+  monIndexQ d = \ case
+    Mon ( Vec [ 1, 0, 0, 0 ] ) -> [|| unT $ _D34_dx     $$d ||]
+    Mon ( Vec [ 0, 1, 0, 0 ] ) -> [|| unT $ _D34_dy     $$d ||]
+    Mon ( Vec [ 0, 0, 1, 0 ] ) -> [|| unT $ _D34_dz     $$d ||]
+    Mon ( Vec [ 0, 0, 0, 1 ] ) -> [|| unT $ _D34_dw     $$d ||]
+    Mon ( Vec [ 2, 0, 0, 0 ] ) -> [|| unT $ _D34_dxdx   $$d ||]
+    Mon ( Vec [ 1, 1, 0, 0 ] ) -> [|| unT $ _D34_dxdy   $$d ||]
+    Mon ( Vec [ 0, 2, 0, 0 ] ) -> [|| unT $ _D34_dydy   $$d ||]
+    Mon ( Vec [ 1, 0, 1, 0 ] ) -> [|| unT $ _D34_dxdz   $$d ||]
+    Mon ( Vec [ 0, 1, 1, 0 ] ) -> [|| unT $ _D34_dydz   $$d ||]
+    Mon ( Vec [ 0, 0, 2, 0 ] ) -> [|| unT $ _D34_dzdz   $$d ||]
+    Mon ( Vec [ 1, 0, 0, 1 ] ) -> [|| unT $ _D34_dxdw   $$d ||]
+    Mon ( Vec [ 0, 1, 0, 1 ] ) -> [|| unT $ _D34_dydw   $$d ||]
+    Mon ( Vec [ 0, 0, 1, 1 ] ) -> [|| unT $ _D34_dzdw   $$d ||]
+    Mon ( Vec [ 0, 0, 0, 2 ] ) -> [|| unT $ _D34_dwdw   $$d ||]
+    Mon ( Vec [ 3, 0, 0, 0 ] ) -> [|| unT $ _D34_dxdxdx $$d ||]
+    Mon ( Vec [ 2, 1, 0, 0 ] ) -> [|| unT $ _D34_dxdxdy $$d ||]
+    Mon ( Vec [ 1, 2, 0, 0 ] ) -> [|| unT $ _D34_dxdydy $$d ||]
+    Mon ( Vec [ 0, 3, 0, 0 ] ) -> [|| unT $ _D34_dydydy $$d ||]
+    Mon ( Vec [ 2, 0, 1, 0 ] ) -> [|| unT $ _D34_dxdxdz $$d ||]
+    Mon ( Vec [ 1, 1, 1, 0 ] ) -> [|| unT $ _D34_dxdydz $$d ||]
+    Mon ( Vec [ 1, 0, 2, 0 ] ) -> [|| unT $ _D34_dxdzdz $$d ||]
+    Mon ( Vec [ 0, 2, 1, 0 ] ) -> [|| unT $ _D34_dydydz $$d ||]
+    Mon ( Vec [ 0, 1, 2, 0 ] ) -> [|| unT $ _D34_dydzdz $$d ||]
+    Mon ( Vec [ 0, 0, 3, 0 ] ) -> [|| unT $ _D34_dzdzdz $$d ||]
+    Mon ( Vec [ 2, 0, 0, 1 ] ) -> [|| unT $ _D34_dxdxdw $$d ||]
+    Mon ( Vec [ 1, 1, 0, 1 ] ) -> [|| unT $ _D34_dxdydw $$d ||]
+    Mon ( Vec [ 0, 2, 0, 1 ] ) -> [|| unT $ _D34_dydydw $$d ||]
+    Mon ( Vec [ 1, 0, 1, 1 ] ) -> [|| unT $ _D34_dxdzdw $$d ||]
+    Mon ( Vec [ 0, 1, 1, 1 ] ) -> [|| unT $ _D34_dydzdw $$d ||]
+    Mon ( Vec [ 0, 0, 2, 1 ] ) -> [|| unT $ _D34_dzdzdw $$d ||]
+    Mon ( Vec [ 1, 0, 0, 2 ] ) -> [|| unT $ _D34_dxdwdw $$d ||]
+    Mon ( Vec [ 0, 1, 0, 2 ] ) -> [|| unT $ _D34_dydwdw $$d ||]
+    Mon ( Vec [ 0, 0, 1, 2 ] ) -> [|| unT $ _D34_dzdwdw $$d ||]
+    Mon ( Vec [ 0, 0, 0, 3 ] ) -> [|| unT $ _D34_dwdwdw $$d ||]
+    _                          -> [||       _D34_v      $$d ||]
diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke.hs b/brush-strokes/src/lib/Math/Bezier/Stroke.hs
index 93b3140..60e6e64 100644
--- a/brush-strokes/src/lib/Math/Bezier/Stroke.hs
+++ b/brush-strokes/src/lib/Math/Bezier/Stroke.hs
@@ -243,7 +243,7 @@ computeStrokeOutline ::
 
      )
   => RootSolvingAlgorithm
-  -> Maybe RootIsolationOptions
+  -> Maybe ( RootIsolationOptions 2 3 )
   -> FitParameters
   -> ( ptData -> usedParams )
   -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing
@@ -522,7 +522,7 @@ outlineFunction
      , Show ptData, Show brushParams
      )
   => RootSolvingAlgorithm
-  -> Maybe RootIsolationOptions
+  -> Maybe ( RootIsolationOptions 2 3 )
   -> ( ptData -> usedParams )
   -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing
   -> ( forall {t} k (i :: t)
@@ -1147,7 +1147,7 @@ solveEnvelopeEquations rootAlgo _t path_t path'_t ( fwdOffset, bwdOffset ) strok
 -- TODO: use Newton's method starting at the midpoint of the box,
 -- instead of just taking the midpoint of the box.
 cuspCoords :: ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) )
-           -> ( Int, Box )
+           -> ( Int, Box 2 )
            -> Cusp
 cuspCoords eqs ( i, 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) )
   | StrokeDatum
@@ -1168,13 +1168,13 @@ cuspCoords eqs ( i, 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) )
 --
 -- See 'isolateRootsIn' for details of the algorithms used.
 findCusps
-  :: RootIsolationOptions
+  :: RootIsolationOptions 2 3
   -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-  -> IntMap ( [ ( Box, RootIsolationTree Box ) ], ( [ Box ], [ Box ] ) )
+  -> IntMap ( [ ( Box 2, RootIsolationTree ( Box 2 ) ) ], ( [ Box 2 ], [ Box 2 ] ) )
 findCusps opts boxStrokeData =
   findCuspsIn opts boxStrokeData $
     IntMap.fromList
-      [ ( i, [ 𝕀 ( ℝ2 zero zero ) ( ℝ2 one one ) ] )
+      [ ( i, [ 𝕀 ( ℝ2 zero zero ) ( ℝ2 one one ) ] ) --[ 𝕀 ( ℝ3 zero zero -1e6 ) ( ℝ3 one one 1e6 ) ] )
       | i <- [ 0 .. length ( boxStrokeData ( 𝕀 ( ℝ1 zero ) ( ℝ1 one ) ) ) - 1 ]
       ]
   where
@@ -1187,10 +1187,10 @@ findCusps opts boxStrokeData =
 -- | Like 'findCusps', but parametrised over the initial boxes for the
 -- root isolation method.
 findCuspsIn
-  :: RootIsolationOptions
+  :: RootIsolationOptions 2 3
   -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-  -> IntMap [ Box ]
-  -> IntMap ( [ ( Box, RootIsolationTree Box ) ], ( [ Box ], [ Box ] ) )
+  -> IntMap [ Box 2 ]
+  -> IntMap ( [ ( Box 2, RootIsolationTree ( Box 2 ) ) ], ( [ Box 2 ], [ Box 2 ] ) )
 findCuspsIn opts boxStrokeData initBoxes =
   IntMap.mapWithKey ( \ i -> foldMap ( isolateRootsIn opts ( eqnPiece i ) ) ) initBoxes
     where
@@ -1211,3 +1211,47 @@ findCuspsIn opts boxStrokeData initBoxes =
         in D12 ( 𝕀 ( ℝ3 ee_lo f_lo g_lo ) ( ℝ3 ee_hi f_hi g_hi ) )
                ( T $ 𝕀 ( ℝ3 ee_t_lo f_t_lo g_t_lo ) ( ℝ3 ee_t_hi f_t_hi g_t_hi ) )
                ( T $ 𝕀 ( ℝ3 ee_s_lo f_s_lo g_s_lo ) ( ℝ3 ee_s_hi f_s_hi g_s_hi ) )
+
+{-
+findCuspsIn
+  :: RootIsolationOptions 3 3
+  -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
+  -> IntMap [ Box 3 ]
+  -> IntMap ( [ ( Box 3, RootIsolationTree ( Box 3 ) ) ], ( [ Box 3 ], [ Box 3 ] ) )
+findCuspsIn opts boxStrokeData initBoxes =
+  IntMap.mapWithKey ( \ i -> foldMap ( isolateRootsIn opts ( eqnPiece i ) ) ) initBoxes
+    where
+      eqnPiece :: Int -> 𝕀 ( ℝ 3 ) -> D1𝔸3 ( 𝕀 ( ℝ 3 ) )
+      eqnPiece i ( 𝕀 ( ℝ3 t_lo s_lo λ_lo ) ( ℝ3 t_hi s_hi λ_hi ) ) =
+        let t = 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi )
+            s = 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi )
+            λ = 𝕀 λ_lo λ_hi
+
+            StrokeDatum
+              { dstroke =
+                 D32 { _D32_dx = T c_t, _D32_dy = T c_s
+                     , _D32_dxdx = T c_tt, _D32_dxdy = T c_ts, _D32_dydy = T c_ss }
+              , ee =
+                D22 { _D22_dx = T ee_t, _D22_dy = T ee_s
+                    , _D22_dxdx = T ee_tt, _D22_dxdy = T ee_ts, _D22_dydy = T ee_ss }
+              } = ( boxStrokeData t `Seq.index` i ) s
+
+            -- λ ∂E/∂s - ∂E/∂t = 0
+            𝕀 ( ℝ1 f1_lo   ) ( ℝ1 f1_hi   ) = unT $ λ *^ T ee_s  ^-^ T ee_t
+            𝕀 ( ℝ1 f1_t_lo ) ( ℝ1 f1_t_hi ) = unT $ λ *^ T ee_ts ^-^ T ee_tt
+            𝕀 ( ℝ1 f1_s_lo ) ( ℝ1 f1_s_hi ) = unT $ λ *^ T ee_ss ^-^ T ee_ts
+            𝕀 ( ℝ1 f1_λ_lo ) ( ℝ1 f1_λ_hi ) = ee_s
+
+            -- λ ∂c/∂s - ∂c/∂t = 0
+            𝕀 ( ℝ2 f2_lo   f3_lo   ) ( ℝ2 f2_hi   f3_hi   ) = unT $ λ *^ T c_s  ^-^ T c_t
+            𝕀 ( ℝ2 f2_t_lo f3_t_lo ) ( ℝ2 f2_t_hi f3_t_hi ) = unT $ λ *^ T c_ts ^-^ T c_tt
+            𝕀 ( ℝ2 f2_s_lo f3_s_lo ) ( ℝ2 f2_s_hi f3_s_hi ) = unT $ λ *^ T c_ss ^-^ T c_ts
+            𝕀 ( ℝ2 f2_λ_lo f3_λ_lo ) ( ℝ2 f2_λ_hi f3_λ_hi ) = c_s
+
+        in D13 ( 𝕀 ( ℝ3 f1_lo f2_lo f3_lo ) ( ℝ3 f1_hi f2_hi f3_hi ) )
+               ( T $ 𝕀 ( ℝ3 f1_t_lo f2_t_lo f3_t_lo ) ( ℝ3 f1_t_hi f2_t_hi f3_t_hi ) )
+               ( T $ 𝕀 ( ℝ3 f1_s_lo f2_s_lo f3_s_lo ) ( ℝ3 f1_s_hi f2_s_hi f3_s_hi ) )
+               ( T $ 𝕀 ( ℝ3 f1_λ_lo f2_λ_lo f3_λ_lo ) ( ℝ3 f1_λ_hi f2_λ_hi f3_λ_hi ) )
+
+
+-}
\ No newline at end of file
diff --git a/brush-strokes/src/lib/Math/Interval.hs b/brush-strokes/src/lib/Math/Interval.hs
index 1ca0582..0ddbe84 100644
--- a/brush-strokes/src/lib/Math/Interval.hs
+++ b/brush-strokes/src/lib/Math/Interval.hs
@@ -179,14 +179,14 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 1 ) where
   konst   :: forall w. AbelianGroup w => w -> D2𝔸1 w
   konst w =
     let !o = fromInteger @w 0
-    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+    in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
 
-  value df = $$( monIndex [|| df ||] zeroMonomial )
+  value df = $$( monIndexQ [|| df ||] zeroMonomial )
 
   linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 1 -> w ) -> 𝕀ℝ 1 -> D2𝔸1 w
   linearD f v =
     let !o = origin @( 𝕀 Double ) @( T w )
-    in $$( monTabulate \ mon ->
+    in $$( monTabulateQ \ mon ->
               if | isZeroMonomial mon
                  -> [|| f v ||]
                  | Just i <- isLinear mon
@@ -214,14 +214,14 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 1 ) where
   konst   :: forall w. AbelianGroup w => w -> D3𝔸1 w
   konst w =
     let !o = fromInteger @w 0
-    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+    in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
 
-  value df = $$( monIndex [|| df ||] zeroMonomial )
+  value df = $$( monIndexQ [|| df ||] zeroMonomial )
 
   linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 1 -> w ) -> 𝕀ℝ 1 -> D3𝔸1 w
   linearD f v =
     let !o = origin @( 𝕀 Double ) @( T w )
-    in $$( monTabulate \ mon ->
+    in $$( monTabulateQ \ mon ->
               if | isZeroMonomial mon
                  -> [|| f v ||]
                  | Just i <- isLinear mon
@@ -249,14 +249,14 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 2 ) where
   konst   :: forall w. AbelianGroup w => w -> D2𝔸2 w
   konst w =
     let !o = fromInteger @w 0
-    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+    in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
 
-  value df = $$( monIndex [|| df ||] zeroMonomial )
+  value df = $$( monIndexQ [|| df ||] zeroMonomial )
 
   linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 2 -> w ) -> 𝕀ℝ 2 -> D2𝔸2 w
   linearD f v =
     let !o = origin @( 𝕀 Double ) @( T w )
-    in $$( monTabulate \ mon ->
+    in $$( monTabulateQ \ mon ->
               if | isZeroMonomial mon
                  -> [|| f v ||]
                  | Just i <- isLinear mon
@@ -284,14 +284,14 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 2 ) where
   konst   :: forall w. AbelianGroup w => w -> D3𝔸2 w
   konst w =
     let !o = fromInteger @w 0
-    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+    in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
 
-  value df = $$( monIndex [|| df ||] zeroMonomial )
+  value df = $$( monIndexQ [|| df ||] zeroMonomial )
 
   linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 2 -> w ) -> 𝕀ℝ 2 -> D3𝔸2 w
   linearD f v =
     let !o = origin @( 𝕀 Double ) @( T w )
-    in $$( monTabulate \ mon ->
+    in $$( monTabulateQ \ mon ->
               if | isZeroMonomial mon
                  -> [|| f v ||]
                  | Just i <- isLinear mon
@@ -319,14 +319,14 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 3 ) where
   konst   :: forall w. AbelianGroup w => w -> D2𝔸3 w
   konst w =
     let !o = fromInteger @w 0
-    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+    in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
 
-  value df = $$( monIndex [|| df ||] zeroMonomial )
+  value df = $$( monIndexQ [|| df ||] zeroMonomial )
 
   linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 3 -> w ) -> 𝕀ℝ 3 -> D2𝔸3 w
   linearD f v =
     let !o = origin @( 𝕀 Double ) @( T w )
-    in $$( monTabulate \ mon ->
+    in $$( monTabulateQ \ mon ->
               if | isZeroMonomial mon
                  -> [|| f v ||]
                  | Just i <- isLinear mon
@@ -354,14 +354,14 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 3 ) where
   konst   :: forall w. AbelianGroup w => w -> D3𝔸3 w
   konst w =
     let !o = fromInteger @w 0
-    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+    in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
 
-  value df = $$( monIndex [|| df ||] zeroMonomial )
+  value df = $$( monIndexQ [|| df ||] zeroMonomial )
 
   linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 3 -> w ) -> 𝕀ℝ 3 -> D3𝔸3 w
   linearD f v =
     let !o = origin @( 𝕀 Double ) @( T w )
-    in $$( monTabulate \ mon ->
+    in $$( monTabulateQ \ mon ->
               if | isZeroMonomial mon
                  -> [|| f v ||]
                  | Just i <- isLinear mon
@@ -389,14 +389,14 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 4 ) where
   konst   :: forall w. AbelianGroup w => w -> D2𝔸4 w
   konst w =
     let !o = fromInteger @w 0
-    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+    in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
 
-  value df = $$( monIndex [|| df ||] zeroMonomial )
+  value df = $$( monIndexQ [|| df ||] zeroMonomial )
 
   linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 4 -> w ) -> 𝕀ℝ 4 -> D2𝔸4 w
   linearD f v =
     let !o = origin @( 𝕀 Double ) @( T w )
-    in $$( monTabulate \ mon ->
+    in $$( monTabulateQ \ mon ->
               if | isZeroMonomial mon
                  -> [|| f v ||]
                  | Just i <- isLinear mon
@@ -424,14 +424,14 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 4 ) where
   konst   :: forall w. AbelianGroup w => w -> D3𝔸4 w
   konst w =
     let !o = fromInteger @w 0
-    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+    in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
 
-  value df = $$( monIndex [|| df ||] zeroMonomial )
+  value df = $$( monIndexQ [|| df ||] zeroMonomial )
 
   linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 4 -> w ) -> 𝕀ℝ 4 -> D3𝔸4 w
   linearD f v =
     let !o = origin @( 𝕀 Double ) @( T w )
-    in $$( monTabulate \ mon ->
+    in $$( monTabulateQ \ mon ->
               if | isZeroMonomial mon
                  -> [|| f v ||]
                  | Just i <- isLinear mon
diff --git a/brush-strokes/src/lib/Math/Linear.hs b/brush-strokes/src/lib/Math/Linear.hs
index d49c687..c5f0554 100644
--- a/brush-strokes/src/lib/Math/Linear.hs
+++ b/brush-strokes/src/lib/Math/Linear.hs
@@ -10,25 +10,29 @@ module Math.Linear
 
   -- * Points and vectors (second version)
   , ℝ(..), T(.., V2, V3, V4)
-  , Fin(..), MFin(..)
+  , Fin(..), MFin(..), universe, coordinates
   , RepDim, RepresentableQ(..)
   , Representable(..), set, injection, projection
-  , Vec(..), (!), find, zipIndices, unzip
+  , Vec(..), (!), find, zipIndices
 
   ) where
 
 -- base
 import Prelude hiding ( unzip )
+import Control.Applicative
+  ( ZipList(..) )
 import Data.Coerce
   ( coerce )
 import Data.Kind
   ( Type )
+import GHC.Exts
+  ( proxy# )
 import GHC.Generics
   ( Generic, Generic1, Generically(..), Generically1(..) )
+import GHC.Stack
+  ( HasCallStack )
 import GHC.TypeNats
-  ( Nat, type (+) )
-import Unsafe.Coerce
-  ( unsafeCoerce )
+  ( Nat, KnownNat, natVal' )
 
 -- acts
 import Data.Act
@@ -114,53 +118,36 @@ pattern V4 x y z w = T ( ℝ4 x y z w )
 
 --------------------------------------------------------------------------------
 
-infixr 5 `VS`
 type Vec :: Nat -> Type -> Type
-data Vec n a where
-  VZ :: Vec 0 a
-  VS :: forall n a. a -> Vec n a -> Vec ( 1 + n ) a
-    -- can't be strict, otherwise we can't conveniently
-    -- unsafeCoerce from lists
+newtype Vec n a = Vec { vecList :: [ a ] }
+  deriving newtype ( Show, Eq, Ord, Functor, Foldable )
+  deriving Applicative via ZipList
 
---deriving stock instance Show a => Show ( Vec n a )
-instance Show a => Show ( Vec n a ) where
-  showsPrec p v = showsPrec p ( unsafeCoerce v :: [ a ] )
+universe :: forall n. KnownNat n => Vec n ( Fin n )
+universe = Vec [ Fin i | i <- [ 1 .. fromIntegral ( natVal' @n proxy# ) ] ]
 
-deriving stock instance Functor     ( Vec n )
-deriving stock instance Foldable    ( Vec n )
-deriving stock instance Traversable ( Vec n )
-
-instance Eq a => Eq ( Vec n a ) where
-  (==) = unsafeCoerce $ (==) @[a]
-
-instance Ord a => Ord ( Vec n a ) where
-  compare = unsafeCoerce $ compare @[a]
-  (<=) = unsafeCoerce $ (<=) @[a]
+coordinates :: forall r u. ( Representable r u ) => u -> Vec ( RepDim u ) r
+coordinates u = fmap ( index u ) $ universe @( RepDim u )
+{-# INLINEABLE coordinates #-}
 
 infixl 9 !
-(!) :: forall l a. Vec l a -> Fin l -> a
-v ! Fin i = ( unsafeCoerce v :: [ a ] ) !! fromIntegral i
+(!) :: forall l a. HasCallStack => Vec l a -> Fin l -> a
+( Vec l ) ! Fin i = l !! ( fromIntegral i - 1 )
 
 find :: forall l a. ( a -> Bool ) -> Vec l a -> MFin l
-find f v = MFin ( find_ 1 v )
+find f ( Vec v ) = MFin ( find_ 1 v )
   where
-    find_ :: Word -> Vec n a -> Word
-    find_ j ( VS a as )
+    find_ :: Word -> [ a ] -> Word
+    find_ j ( a : as )
       | f a
       = j
       | otherwise
       = find_ ( j + 1 ) as
-    find_ _ VZ = 0
+    find_ _ [] = 0
 
 zipIndices :: forall n a. Vec n a -> [ ( Fin n, a ) ]
-zipIndices = zipIndices_ 1
+zipIndices ( Vec v ) = zipIndices_ 1 v
   where
-    zipIndices_ :: forall i. Word -> Vec i a -> [ ( Fin n, a ) ]
-    zipIndices_ _ VZ = []
-    zipIndices_ i (a `VS` as) = ( Fin i, a ) : zipIndices_ ( i + 1 ) as
-
-unzip :: Vec n ( a, b ) -> ( Vec n a, Vec n b )
-unzip VZ = ( VZ, VZ )
-unzip ( ( a, b ) `VS` ab ) =
-  case unzip ab of
-    ( as, bs ) -> ( a `VS` as, b `VS` bs )
+    zipIndices_ ::  Word -> [ a ] -> [ ( Fin n, a ) ]
+    zipIndices_ _ [] = []
+    zipIndices_ i (a : as) = ( Fin i, a ) : zipIndices_ ( i + 1 ) as
diff --git a/brush-strokes/src/lib/Math/Linear/Internal.hs b/brush-strokes/src/lib/Math/Linear/Internal.hs
index f0b1bbf..d01d724 100644
--- a/brush-strokes/src/lib/Math/Linear/Internal.hs
+++ b/brush-strokes/src/lib/Math/Linear/Internal.hs
@@ -1,5 +1,6 @@
-{-# LANGUAGE PolyKinds       #-}
-{-# LANGUAGE TemplateHaskell #-}
+{-# LANGUAGE PolyKinds           #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TemplateHaskell     #-}
 
 module Math.Linear.Internal
   ( ℝ(..)
@@ -22,7 +23,7 @@ import GHC.Generics
 import GHC.Show
   ( showSpace )
 import GHC.TypeNats
-  ( Nat )
+  ( Nat, KnownNat )
 
 -- template-haskell
 import Language.Haskell.TH
@@ -94,7 +95,7 @@ instance Show ( ℝ 4 ) where
 -- | 1, ..., n
 type Fin :: Nat -> Type
 newtype Fin n = Fin Word
-  deriving stock Eq
+  deriving stock ( Eq, Ord, Show )
 
 -- | 0, ..., n
 type MFin :: Nat -> Type
@@ -109,16 +110,16 @@ class RepresentableQ r v | v -> r where
   indexQ    :: CodeQ v -> Fin ( RepDim v ) -> CodeQ r
 
 type Representable :: Type -> Type -> Constraint
-class Representable r v | v -> r where
+class KnownNat ( RepDim v ) => Representable r v | v -> r where
   tabulate :: ( Fin ( RepDim v ) -> r ) -> v
   index    :: v -> Fin ( RepDim v ) -> r
 
-{-# INLINEABLE set #-}
 set :: Representable r v => Fin ( RepDim v ) -> r -> v -> v
 set i r u = tabulate \ j ->
   if i == j
   then r
   else index u j
+{-# INLINEABLE set #-}
 
 projection :: ( Representable r u, Representable r v )
            => ( Fin ( RepDim v ) -> Fin ( RepDim u ) )
diff --git a/brush-strokes/src/lib/Math/Monomial.hs b/brush-strokes/src/lib/Math/Monomial.hs
index ee2be74..2356202 100644
--- a/brush-strokes/src/lib/Math/Monomial.hs
+++ b/brush-strokes/src/lib/Math/Monomial.hs
@@ -7,9 +7,9 @@
 
 module Math.Monomial
   ( Mon(..)
-  , MonomialBasis(..), Deg, Vars
+  , MonomialBasis(..), MonomialBasisQ(..), Deg, Vars
   , zeroMonomial, isZeroMonomial
-  , totalDegree, isLinear
+  , totalDegree, isLinear, linearMonomial
 
   , split, mons
 
@@ -57,40 +57,44 @@ type family Vars f
 zeroMonomial :: forall k n. KnownNat n => Mon k n
 zeroMonomial = Mon $ unsafeCoerce @[ Word ] @( Vec n Word )
                    $ replicate ( fromIntegral $ word @n ) 0
+{-# INLINEABLE zeroMonomial #-}
 
 isZeroMonomial :: Mon k n -> Bool
-isZeroMonomial ( Mon VZ ) = True
-isZeroMonomial ( Mon ( i `VS` is ) )
-  | i == 0
-  = isZeroMonomial ( Mon is )
-  | otherwise
-  = False
+isZeroMonomial ( Mon ( Vec pows ) ) = all ( == 0 ) pows
 
 totalDegree :: Mon k n -> Word
-totalDegree ( Mon VZ ) = 0
-totalDegree ( Mon ( i `VS` is ) ) = i + totalDegree ( Mon is )
+totalDegree ( Mon ( Vec pows ) ) = sum pows
+
+linearMonomial :: forall k n. ( KnownNat n, 1 <= k ) => Fin n -> Mon k n
+linearMonomial ( Fin i' ) =
+  Mon $ unsafeCoerce @[ Word ] @( Vec n Word )
+      $ replicate ( i - 1 ) 0 ++ [ 1 ] ++ replicate ( n - i ) 0
+  where
+    n, i :: Int
+    n = fromIntegral $ word @n
+    i = fromIntegral i'
+{-# INLINEABLE linearMonomial #-}
 
 isLinear :: Mon k n -> Maybe ( Fin n )
-isLinear = fmap Fin . go 1
+isLinear ( Mon ( Vec pows ) ) = fmap Fin $ go 1 pows
   where
-    go :: forall k' n'. Word -> Mon k' n' -> Maybe Word
-    go _ ( Mon VZ )
-      = Nothing
-    go j ( Mon ( i `VS` is ) )
+    go :: Word -> [ Word ] -> Maybe Word
+    go _ [] = Nothing
+    go j ( i : is )
       | i == 0
-      = go ( j + 1 ) ( Mon is )
-      | i == 1 && isZeroMonomial ( Mon is )
+      = go ( j + 1 ) is
+      | i == 1 && all ( == 0 ) is
       = Just j
       | otherwise
       = Nothing
 
--- | Comultiplication of monomials.
+-- | Co-multiplication of monomials.
 split :: Mon k n -> [ ( Mon k n, Mon k n ) ]
-split ( Mon VZ ) = [ ( Mon VZ, Mon VZ ) ]
-split ( Mon ( d `VS` ds ) ) =
-      [ ( Mon ( i `VS` as ), Mon ( ( d - i ) `VS` bs ) )
+split ( Mon ( Vec [] ) ) = [ ( Mon ( Vec [] ), Mon ( Vec [] ) ) ]
+split ( Mon ( Vec ( d : ds ) ) ) =
+      [ ( Mon ( Vec ( i : as ) ), Mon ( Vec ( ( ( d - i ) : bs ) ) ) )
       | i <- [ 0 .. d ]
-      , ( Mon as, Mon bs ) <- ( split ( Mon ds ) )
+      , ( Mon ( Vec as ), Mon ( Vec bs ) ) <- ( split ( Mon ( Vec ds ) ) )
       ]
 
 -- | All monomials of degree less than or equal to @k@ in @n@ variables,
@@ -105,12 +109,12 @@ mons' 0 n = [ replicate ( fromIntegral n ) 0 ]
 mons' k n = [ i : is | i <- reverse [ 0 .. k ], is <- mons' ( k - i ) ( n - 1 ) ]
 
 subs :: Mon k n -> [ ( Mon k n, Word ) ]
-subs ( Mon VZ ) = [ ( Mon VZ, maxBound ) ]
-subs ( Mon ( i `VS` is ) )
-  = [ ( Mon ( j `VS` js )
+subs ( Mon ( Vec []) ) = [ ( Mon ( Vec [] ), maxBound ) ]
+subs ( Mon ( Vec ( i : is ) ) )
+  = [ ( Mon ( Vec ( j : js ) )
       , if j == 0 then mult else min ( i `quot` j ) mult )
     | j <- [ 0 .. i ]
-    , ( Mon js, mult ) <- subs ( Mon is )
+    , ( Mon ( Vec js ), mult ) <- subs ( Mon ( Vec is ) )
     ]
 
 word :: forall n. KnownNat n => Word
@@ -121,8 +125,8 @@ factorial :: Word -> Word
 factorial i = product [ 1 .. i ]
 
 vecFactorial :: Vec n Word -> Word
-vecFactorial VZ = 1
-vecFactorial ( i `VS` is ) = factorial i * vecFactorial is
+vecFactorial ( Vec [] ) = 1
+vecFactorial ( Vec ( i : is ) ) = factorial i * vecFactorial ( Vec is )
 
 --------------------------------------------------------------------------------
 -- Computations for the chain rule R^1 -> R^n -> R^m
@@ -160,13 +164,13 @@ multiSubsetsSum :: forall n
 multiSubsetsSum is = goMSS
   where
     goMSS :: forall i. Word -> Vec i Word -> [ Vec i [ ( Word, Word ) ] ]
-    goMSS 0 VZ = [ VZ ]
-    goMSS _ VZ = [ ]
-    goMSS s (n `VS` ns) =
-      [ multi `VS` rest
+    goMSS 0 ( Vec [] ) = [ Vec [] ]
+    goMSS _ ( Vec [] ) = [ ]
+    goMSS s ( Vec ( n : ns ) ) =
+      [ Vec ( multi : rest )
       | s_i <- [ n * i_min .. s ]
       , multi <- multiSubsetSum n s_i is
-      , rest <- goMSS ( s - s_i ) ns ]
+      , Vec rest <- goMSS ( s - s_i ) ( Vec ns ) ]
     i_min = case is of
       [] -> 0
       _  -> max 0 $ minimum is
@@ -224,19 +228,27 @@ subPart' m = zipWith ( \ i j -> i - m * j )
 -- monomials in @Vars u@ variables, of degree up to (and including) @Deg u@.
 type MonomialBasis :: ( Type -> Type ) -> Constraint
 class MonomialBasis f where
-  monTabulate :: ( ( Mon ( Deg f ) ( Vars f ) ) -> CodeQ u ) -> CodeQ ( f u )
-  monIndex :: CodeQ ( f u ) -> Mon ( Deg f ) ( Vars f ) -> CodeQ u
+  monTabulate :: ( ( Mon ( Deg f ) ( Vars f ) ) -> u ) -> f u
+  monIndex :: f u -> Mon ( Deg f ) ( Vars f ) -> u
+
+-- | @'MonomialBasisQ' f@ exhibits @f u@ as a free @r@-module with basis the
+-- monomials in @Vars u@ variables, of degree up to (and including) @Deg u@,
+-- at compile-time (to be spliced in using TemplateHaskell).
+type MonomialBasisQ :: ( Type -> Type ) -> Constraint
+class MonomialBasisQ f where
+  monTabulateQ :: ( ( Mon ( Deg f ) ( Vars f ) ) -> CodeQ u ) -> CodeQ ( f u )
+  monIndexQ :: CodeQ ( f u ) -> Mon ( Deg f ) ( Vars f ) -> CodeQ u
 
 --------------------------------------------------------------------------------
 
-prodRuleQ :: forall f r. MonomialBasis f
+prodRuleQ :: forall f r. MonomialBasisQ f
           => CodeQ r -> CodeQ ( r -> r -> r ) -> CodeQ ( r -> r -> r )
             -- Ring r constraint (circumvent TH constraint problem)
           -> CodeQ ( f r ) -> CodeQ ( f r ) -> CodeQ ( f r )
 prodRuleQ zero plus times df1 df2 =
-  monTabulate @f \ mon ->
+  monTabulateQ @f \ mon ->
     [|| $$( foldQ plus zero
-           [ [|| $$times $$( monIndex @f df1 m1 ) $$( monIndex @f df2 m2 ) ||]
+           [ [|| $$times $$( monIndexQ @f df1 m1 ) $$( monIndexQ @f df2 m2 ) ||]
            | (m1, m2) <- split mon
            ] )
     ||]
diff --git a/brush-strokes/src/lib/Math/Root/Isolation.hs b/brush-strokes/src/lib/Math/Root/Isolation.hs
index c603b08..5d18c44 100644
--- a/brush-strokes/src/lib/Math/Root/Isolation.hs
+++ b/brush-strokes/src/lib/Math/Root/Isolation.hs
@@ -1,3 +1,5 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+
 module Math.Root.Isolation
   ( -- * Root isolation using interval arithmetic
     Box, BoxHistory
@@ -22,7 +24,6 @@ module Math.Root.Isolation
     -- ** Trees recording search space of root isolation algorithms
   , RootIsolationTree(..), showRootIsolationTree
   , RootIsolationStep(..)
-  , RootIsolationLeaf(..)
   )
   where
 
@@ -33,19 +34,35 @@ import Control.Monad
   ( when )
 import Data.Bifunctor
   ( Bifunctor(bimap) )
+import Data.Coerce
+  ( coerce )
+import Data.Kind
+  ( Type )
+import Data.Foldable
+ ( toList )
 import Data.List
-  ( nub, partition)
+  ( nub, partition, sort )
 import Data.List.NonEmpty
   ( NonEmpty )
 import qualified Data.List.NonEmpty as NE
-  ( NonEmpty(..), last )
+  ( NonEmpty(..), last, singleton )
+import Data.Semigroup
+  ( Arg(..), Dual(..) )
 import Numeric
   ( showFFloat )
+import GHC.TypeNats
+  ( Nat, KnownNat )
 
 -- containers
 import Data.Tree
   ( Tree(..) )
 
+-- eigen
+import qualified Eigen.Matrix as Eigen
+  ( Matrix
+  , determinant, generate, inverse, unsafeCoeff
+  )
+
 -- transformers
 import Control.Monad.Trans.State.Strict as State
   ( State, evalState, get, put )
@@ -54,7 +71,7 @@ import Control.Monad.Trans.Writer.CPS
 
 -- MetaBrush
 import Math.Algebra.Dual
-  ( D, D1𝔸2(..) )
+  ( D )
 import Math.Epsilon
   ( nearZero )
 import Math.Float.Utils
@@ -63,6 +80,10 @@ import Math.Interval
 import Math.Linear
 import Math.Module
   ( Module(..) )
+import Math.Monomial
+  ( MonomialBasis(..), Deg, Vars
+  , linearMonomial, zeroMonomial
+  )
 import qualified Math.Ring as Ring
 
 --import Debug.Utils
@@ -71,19 +92,19 @@ import qualified Math.Ring as Ring
 
 -- | A tree recording the steps taken when doing cusp finding.
 data RootIsolationTree d
-  = RootIsolationLeaf (RootIsolationLeaf d)
+  = RootIsolationLeaf String d
   | RootIsolationStep RootIsolationStep [(d, RootIsolationTree d)]
 
 -- | A description of a step taken in cusp-finding.
 data RootIsolationStep
   = GaussSeidelStep
-  | BisectionStep (String, Double)
+  | BisectionStep String Double
   | Box1Step
   | Box2Step
 
 instance Show RootIsolationStep where
   showsPrec _ GaussSeidelStep = showString "GS"
-  showsPrec _ ( BisectionStep (coord, w) )
+  showsPrec _ ( BisectionStep coord w )
     = showString "bis "
     . showParen True
       ( showString coord
@@ -93,56 +114,81 @@ instance Show RootIsolationStep where
   showsPrec _ Box1Step = showString "box(1)"
   showsPrec _ Box2Step = showString "box(2)"
 
-data RootIsolationLeaf d
-  = TooSmall d
-  deriving stock Show
-
-showRootIsolationTree :: Box -> RootIsolationTree Box -> Tree String
-showRootIsolationTree cand (RootIsolationLeaf l) = Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ show l) []
+showRootIsolationTree
+  :: ( Representable Double ( ℝ n ), Show ( Box n ) )
+  => Box n -> RootIsolationTree ( Box n ) -> Tree String
+showRootIsolationTree cand (RootIsolationLeaf why l) = Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ why ++ " " ++ show l) []
 showRootIsolationTree cand (RootIsolationStep s ts)
   = Node (show cand ++ " abc " ++ showArea (boxArea cand) ++ " " ++ show s) $ map (\ (c,t) -> showRootIsolationTree c t) ts
 
-boxArea :: Box -> Double
-boxArea ( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) )
-  = abs ( t_hi - t_lo ) * abs ( s_hi - s_lo )
+boxArea :: Representable Double ( ℝ n ) => Box n -> Double
+boxArea ( 𝕀 lo hi ) =
+  product ( ( \ l h -> abs ( h - l ) ) <$> coordinates lo <*> coordinates hi )
 
 showArea :: Double -> String
 showArea area = "(area " ++ showFFloat (Just 6) area "" ++ ")"
 
-type Box = 𝕀ℝ 2
-type BoxHistory = [ NE.NonEmpty ( RootIsolationStep, Box ) ]
+type Box n = 𝕀ℝ n
+type BoxHistory n = [ NE.NonEmpty ( RootIsolationStep, Box n ) ]
+type BoxCt n d = ( n ~ 2, d ~ 3 )
+{-
+  ( Show ( 𝕀ℝ n ), Show ( ℝ n )
+  , Eq ( ℝ n )
+  , Representable Double ( ℝ n )
+  , MonomialBasis ( D 1 ( ℝ n ) )
+  , Deg ( D 1 ( ℝ n ) ) ~ 1
+  , Vars ( D 1 ( ℝ n ) ) ~ n
+  , Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) )
+  , Applicative ( Vec n )
+
+  , Ord ( ℝ d )
+  , Module Double ( T ( ℝ d ) )
+  , Representable Double ( ℝ d )
+  )
+-}
 
 -- | Options for the root isolation methods in 'findCusps' and 'isolateRootsIn'.
-data RootIsolationOptions
+type RootIsolationOptions :: Nat -> Nat -> Type
+data RootIsolationOptions n d
   = RootIsolationOptions
   { -- | Minimum width of a box.
     minWidth :: !Double
     -- | Given a box and its history, return a round of cusp-finding strategies
     -- to run, in sequence, on this box.
-  , cuspFindingAlgorithms :: !( Box -> BoxHistory -> NonEmpty RootIsolationAlgorithm )
+    --
+    -- Returning an empty list means: give up on this box.
+  , cuspFindingAlgorithms :: !( Box n -> BoxHistory n -> Either String ( NE.NonEmpty ( RootIsolationAlgorithm n d ) ) )
   }
 
-data RootIsolationAlgorithm
+type RootIsolationAlgorithm :: Nat -> Nat -> Type
+data RootIsolationAlgorithm n d
    -- | Bisection step.
-  = Bisection !BisectionOptions
+  = Bisection !( BisectionOptions n d )
   -- | Gauss–Seidel step with the given preconditioning method.
-  | GaussSeidel !GaussSeidelOptions
+  | GaussSeidel !( GaussSeidelOptions n d )
   -- | @box(1)@-consistency.
   | Box1 !Box1Options
   -- | @box(2)@-consistency.
   | Box2 !Box2Options
 
 -- | Options for the bisection method.
-data BisectionOptions =
+type BisectionOptions :: Nat -> Nat -> Type
+data BisectionOptions n d =
   BisectionOptions
-  { canHaveSols :: !( ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) ) -> Box -> Bool )
-  , fallbackBisectionDim :: !( [ ( RootIsolationStep, Box ) ] -> BoxHistory -> ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) ) -> ( Int, String ) )
+  { canHaveSols :: !( ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) -> Box n -> Bool )
+  , fallbackBisectionDim :: !( [ ( RootIsolationStep, Box n ) ] -> BoxHistory n -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) -> ( Fin n, String ) )
   }
 
 -- | Options for the interval Gauss–Seidel method.
-data GaussSeidelOptions =
+type GaussSeidelOptions :: Nat -> Nat -> Type
+data GaussSeidelOptions n d =
   GaussSeidelOptions
-    { gsPreconditioner :: !Preconditioner }
+    { -- | Which preconditioner to user?
+      gsPreconditioner :: !Preconditioner
+      -- | Function that projects over the equations we will consider
+      -- (the identity for a well-determined problem, or a project for
+      -- an overdetermined system).
+    , gsDims :: ( 𝕀ℝ d -> 𝕀ℝ n ) }
 
 -- | Options for the @box(1)@-consistency method.
 data Box1Options =
@@ -156,7 +202,7 @@ data Box2Options =
     , box2LambdaMin :: !Double
     }
 
-defaultRootIsolationOptions :: RootIsolationOptions
+defaultRootIsolationOptions :: RootIsolationOptions 2 3
 defaultRootIsolationOptions =
   RootIsolationOptions
     { minWidth
@@ -166,68 +212,119 @@ defaultRootIsolationOptions =
     minWidth = 1e-5
     narrowAbs = 5e-3
 
-defaultRootIsolationAlgorithms :: Double -> Double -> Box -> BoxHistory -> NonEmpty RootIsolationAlgorithm
-defaultRootIsolationAlgorithms minWidth narrowAbs box history =
-  case history of
-    lastRoundBoxes : _
-      -- If, in the last round of strategies, we didn't try bisection...
-      | any ( \case { BisectionStep {} -> False; _ -> True } . fst ) lastRoundBoxes
-      , (_, lastRoundFirstBox ) <- NE.last lastRoundBoxes
-      -- ...and the last round didn't sufficiently reduce the size of the box...
-      , not $ box `sufficientlySmallerThan` lastRoundFirstBox
-      -- ...then try bisecting the box.
-      -> Bisection ( defaultBisectionOptions minWidth narrowAbs box ) NE.:| []
-
-    -- Otherwise, do a normal round.
-    -- Currently: we try an interval Gauss–Seidel step followed by box(1)-consistency.
-    _ -> GaussSeidel defaultGaussSeidelOptions
-          NE.:| [ Box1 ( Box1Options { box1EpsEq = narrowAbs } ) ]
+defaultRootIsolationAlgorithms
+  :: forall n d
+  .  BoxCt n d
+  => Double -> Double -> Box n -> BoxHistory n -> Either String ( NE.NonEmpty ( RootIsolationAlgorithm n d ) )
+defaultRootIsolationAlgorithms minWidth narrowAbs box history
+  -- All the box widths are lower than the minimum width: give up.
+  | and $ ( \ cd -> width cd <= minWidth ) <$> coordinates box
+  = Left $ "sizes <= " ++ show minWidth
+  -- Give up after depth 50
+  -- | length history >= 50
+  -- = Left "depth >= 50"
+  | otherwise
+  = Right $ case history of
+      lastRoundBoxes : _
+        -- If, in the last round of strategies, we didn't try bisection...
+        | any ( \case { BisectionStep {} -> False; _ -> True } . fst ) lastRoundBoxes
+        , ( _step, lastRoundFirstBox ) <- NE.last lastRoundBoxes
+        -- ...and the last round didn't sufficiently reduce the size of the box...
+        , not $ box `sufficientlySmallerThan` lastRoundFirstBox
+        -- ...then try bisecting the box.
+        -> NE.singleton $ Bisection ( defaultBisectionOptions minWidth narrowAbs box )
+      -- Otherwise, do a normal round.
+      -- Currently: we try an interval Gauss–Seidel step followed by box(1)-consistency.
+      _ -> GaussSeidel defaultGaussSeidelOptions
+           NE.:| [ Box1 ( Box1Options { box1EpsEq = narrowAbs } ) ]
   where
-    sufficientlySmallerThan :: Box -> Box -> Bool
-    sufficientlySmallerThan
-      ( 𝕀 ( ℝ2 t1_lo s1_lo ) ( ℝ2 t1_hi s1_hi ) )
-      ( 𝕀 ( ℝ2 t2_lo s2_lo ) ( ℝ2 t2_hi s2_hi ) ) =
-        ( ( t1_hi - t1_lo ) - ( t2_hi - t2_lo ) > narrowAbs )
-          ||
-        ( ( s1_hi - s1_lo ) - ( s2_hi - s2_lo ) > narrowAbs )
+    -- Did we reduce the box width by at least "narrowAbs" in at least one of the dimensions?
+    sufficientlySmallerThan :: Box n -> Box n -> Bool
+    b1 `sufficientlySmallerThan` b2 =
+      or $
+        ( \ cd1 cd2 -> width cd2 - width cd1 > narrowAbs )
+          <$> coordinates b1
+          <*> coordinates b2
+{-# INLINEABLE defaultRootIsolationAlgorithms #-}
 
-defaultGaussSeidelOptions :: GaussSeidelOptions
-defaultGaussSeidelOptions = GaussSeidelOptions { gsPreconditioner = InverseMidJacobian }
+defaultGaussSeidelOptions :: GaussSeidelOptions 2 3
+defaultGaussSeidelOptions =
+  GaussSeidelOptions
+    { gsPreconditioner = InverseMidJacobian
+    , gsDims = --id
+         \ ( 𝕀 ( ℝ3 _a_lo b_lo c_lo ) ( ℝ3 _a_hi b_hi c_hi ) )
+        -> 𝕀 ( ℝ2 b_lo c_lo ) ( ℝ2 b_hi c_hi )
+}
 
-defaultBisectionOptions :: Double -> Double -> Box -> BisectionOptions
 defaultBisectionOptions
-  _minWidth
-  _narrowAbs
-  box@( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) )
-  = BisectionOptions
-      { canHaveSols =
-          \ eqs box' ->
-            -- box(0)-consistency
-            let D12 iRange _ _ = eqs box'
-            in ℝ3 0 0 0 `inside` iRange
+  :: forall n d
+  .  BoxCt n d
+  => Double -> Double
+  -> Box n -> BisectionOptions n d
+defaultBisectionOptions _minWidth _narrowAbs box =
+  BisectionOptions
+    { canHaveSols =
+        \ eqs box' ->
+          -- box(0)-consistency
+          let iRange' :: Box d
+              iRange' = ( `monIndex` zeroMonomial ) $ eqs box'
+          in unT ( origin @Double ) `inside` iRange'
 
-            -- box(1)-consistency
-            --not $ null $ makeBox1Consistent eqs _minWidth _narrowAbs box'
+          -- box(1)-consistency
+          --not $ null $ makeBox1Consistent eqs _minWidth _narrowAbs box'
 
-            -- box(2)-consistency
-            --let box'' = makeBox2Consistent eqs _minWidth _narrowAbs 0.2 box'
-            --    D12 iRange'' = eqs box''
-            --in ℝ3 0 0 0 `inside` iRange''
-      , fallbackBisectionDim =
-          \ _roundHist _prevRoundsHist eqs ->
-            let D12 _ ( T f_t ) ( T f_s ) = eqs box
-                wd_t = t_hi - t_lo
-                wd_s = s_hi - s_lo
-                tJWidth = wd_t * normVI3 f_t
-                sJWidth = wd_s * normVI3 f_s
-            in if tJWidth >= sJWidth
-                    -- bisect along dimension that maximises  width(x_j) * || J_j || ...
-                    -- ... but don't allow thin boxes
-                    || ( wd_t >= 10 * wd_s )
-                    && not ( wd_s >= 10 * wd_t )
-               then (0, "")
-               else (1, "")
-      }
+          -- box(2)-consistency
+          --let box'' = makeBox2Consistent eqs _minWidth _narrowAbs 0.2 box'
+          --    iRange'' :: Box d
+          --    iRange'' = value @Double @1 @( Box n ) $ eqs box''
+          --in origin @Double `inside` iRange''
+    , fallbackBisectionDim =
+        \ _roundHist _prevRoundsHist eqs ->
+          let df = eqs box
+              datPerCoord =
+                [ CoordBisectionData
+                   { coordIndex = i
+                   , coordInterval = box `index` i
+                   , coordJacobianColumn = df `monIndex` ( linearMonomial i )
+                   }
+                | i <- toList $ universe @n ]
+
+              -- First, check if the largest dimension is over 10 times larger
+              -- than the smallest dimension; if so bisect along that coordinate.
+          in case sortOnArg ( width . coordInterval ) datPerCoord of
+                [] -> error "dimension 0"
+                [Arg _ d] -> (coordIndex d, "")
+                Arg w0 _ : ds ->
+                  let Arg w1 d1 = last ds
+                  in if w1 >= 10 * w0
+                     then ( coordIndex d1, "tooWide")
+                     -- Otherwise, pick the dimension with the largest spread = width * Jacobian column norm
+                     else case sortOnArg ( Dual . spread ) datPerCoord of
+                       [] -> error "dimension 0"
+                       Arg _ d : _ ->
+                        (coordIndex d, "spread")
+    }
+{-# INLINEABLE defaultBisectionOptions #-}
+
+sortOnArg :: Ord b => (a -> b) -> [a] -> [Arg b a]
+sortOnArg f = sort . map ( \ a -> Arg ( f a ) a )
+{-# INLINEABLE sortOnArg #-}
+
+type CoordBisectionData :: Nat -> Nat -> Type
+data CoordBisectionData n d =
+  CoordBisectionData
+    { coordIndex :: !( Fin n )
+    , coordInterval :: !( 𝕀 Double )
+    , coordJacobianColumn :: !( 𝕀ℝ d )
+    }
+deriving stock instance Show ( ℝ d )
+                     => Show ( CoordBisectionData n d )
+
+spread :: ( BoxCt n d, Representable Double ( ℝ d ) )
+       => CoordBisectionData n d -> Double
+spread ( CoordBisectionData { coordInterval = cd, coordJacobianColumn = j_cd } )
+  = width cd * normVI j_cd
+{-# INLINEABLE spread #-}
 
 -- | Use the following algorithms using interval arithmetic in order
 -- to isolate roots:
@@ -243,10 +340,12 @@ defaultBisectionOptions
 -- will converge starting from anywhere inside the box), and @dunno@ are small
 -- boxes which might or might not contain solutions.
 isolateRootsIn
-  :: RootIsolationOptions
-  -> ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) )
-  -> Box
-  -> ( [ ( Box, RootIsolationTree Box ) ], ( [ Box ], [ Box ] ) )
+  :: forall n d
+  .  BoxCt n d
+  => RootIsolationOptions n d
+  -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
+  -> Box n
+  -> ( [ ( Box n, RootIsolationTree ( Box n ) ) ], ( [ Box n ], [ Box n ] ) )
 isolateRootsIn
   ( RootIsolationOptions
     { minWidth
@@ -257,40 +356,41 @@ isolateRootsIn
 
   where
 
-    go :: BoxHistory
-       -> Box -- box to work on
-       -> Writer ( [ Box ], [ Box ] )
-                 [ RootIsolationTree Box ]
-    go history
-       cand@( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) )
-      | -- Check the envelope equation interval contains zero.
-        not $ ( ℝ3 0 0 0 `inside` iRange )
-      -- Box doesn't contain a solution: discard it.
+    go :: BoxHistory n
+       -> Box n -- box to work on
+       -> Writer ( [ Box n ], [ Box n ] )
+                 [ RootIsolationTree ( Box n ) ]
+    go history cand
+      | -- Check the range of the equations contains zero.
+        not $ ( unT ( origin @Double ) `inside` iRange )
+        -- Box doesn't contain a solution: discard it.
       = return []
-      -- Box is small: stop processing it.
-      | t_hi - t_lo < minWidth && s_hi - s_lo < minWidth
-      = do let res = TooSmall cand
-           tell ( [ cand ], [] )
-           return [ RootIsolationLeaf res ]
       | otherwise
-      = doStrategies history ( cuspFindingAlgorithms cand history ) cand
+      = case cuspFindingAlgorithms cand history of
+          Right strats -> doStrategies history strats cand
+          Left whyStop -> do
+            -- We are giving up on this box (e.g. because it is too small,
+            -- or we have reached an iteration depth).
+            tell ( [ cand ], [] )
+            return [ RootIsolationLeaf whyStop cand ]
       where
-        D12 iRange _ _ = eqs $ 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi )
+        iRange :: Box d
+        iRange = ( `monIndex` zeroMonomial ) $ eqs cand
 
     -- Run a round of cusp finding strategies, then recur.
     doStrategies
-      :: BoxHistory
-      -> NonEmpty RootIsolationAlgorithm
-      -> Box
-      -> Writer ( [ Box ], [ Box ] )
-                [ RootIsolationTree Box ]
+      :: BoxHistory n
+      -> NonEmpty ( RootIsolationAlgorithm n d )
+      -> Box n
+      -> Writer ( [ Box n ], [ Box n ] )
+                [ RootIsolationTree ( Box n ) ]
     doStrategies prevRoundsHist = do_strats []
       where
-        do_strats :: [ ( RootIsolationStep, Box ) ]
-                  -> NE.NonEmpty RootIsolationAlgorithm
-                  -> Box
-                  -> Writer ( [ Box ], [ Box ] )
-                            [ RootIsolationTree Box ]
+        do_strats :: [ ( RootIsolationStep, Box n ) ]
+                  -> NE.NonEmpty ( RootIsolationAlgorithm n d )
+                  -> Box n
+                  -> Writer ( [ Box n ], [ Box n ] )
+                            [ RootIsolationTree ( Box n )]
         do_strats thisRoundHist ( algo NE.:| algos ) box = do
           -- Run one strategy in the round.
           ( step, boxes ) <- doStrategy thisRoundHist prevRoundsHist eqs minWidth algo box
@@ -307,116 +407,137 @@ isolateRootsIn
               in recur step ( go history ) boxes
 
     recur :: RootIsolationStep
-          -> ( Box -> Writer ( [ Box ], [ Box ] ) [ RootIsolationTree Box ] )
-          -> [ Box ]
-          -> Writer ( [Box], [Box] ) [ RootIsolationTree Box ]
+          -> ( Box n -> Writer ( [ Box n ], [ Box n ] ) [ RootIsolationTree ( Box n ) ] )
+          -> [ Box n ]
+          -> Writer ( [ Box n ], [ Box n ] ) [ RootIsolationTree ( Box n ) ]
     recur step r cands = do
       rest <- traverse ( \ c -> do { trees <- r c; return [ (c, tree) | tree <- trees ] } ) cands
       return [ RootIsolationStep step $ concat rest ]
-
+{-# INLINEABLE isolateRootsIn #-}
 
 -- | Execute a cusp-finding strategy, replacing the input box with
 -- (hopefully smaller) output boxes.
-doStrategy :: [ ( RootIsolationStep, Box ) ]
-           -> BoxHistory
-           -> ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) )
-           -> Double
-           -> RootIsolationAlgorithm
-           -> Box
-           -> Writer ( [ Box ], [ Box ] )
-                     ( RootIsolationStep, [ Box ] )
+doStrategy
+  :: BoxCt n d
+  => [ ( RootIsolationStep, Box n ) ]
+  -> BoxHistory n
+  -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
+  -> Double
+  -> RootIsolationAlgorithm n d
+  -> Box n
+  -> Writer ( [ Box n ], [ Box n ] )
+            ( RootIsolationStep, [ Box n ] )
 doStrategy roundHistory previousRoundsHistory eqs minWidth algo box =
   case algo of
-    GaussSeidel ( GaussSeidelOptions { gsPreconditioner } ) -> do
-      boxes <- intervalGaussSeidel eqs gsPreconditioner box
+    GaussSeidel gsOptions -> do
+      boxes <- intervalGaussSeidel gsOptions eqs box
       return ( GaussSeidelStep, boxes )
     Box1 ( Box1Options { box1EpsEq } ) ->
       return ( Box1Step, makeBox1Consistent eqs minWidth box1EpsEq box )
     Box2 ( Box2Options { box2LambdaMin, box2EpsEq } ) ->
       return ( Box2Step, [ makeBox2Consistent eqs minWidth box2EpsEq box2LambdaMin box ] )
     Bisection ( BisectionOptions { canHaveSols, fallbackBisectionDim } ) -> do
-      let ( boxes, whatBis ) = bisect ( canHaveSols eqs ) ( fallbackBisectionDim roundHistory previousRoundsHistory eqs ) box
-      return ( BisectionStep whatBis, boxes )
+      let ( boxes, ( whatBis, mid ) ) = bisect ( canHaveSols eqs ) ( fallbackBisectionDim roundHistory previousRoundsHistory eqs ) box
+      return ( BisectionStep whatBis mid, boxes )
+{-# INLINEABLE doStrategy #-}
 
 -- | Bisect the given box.
 --
 -- (The difficult part lies in determining along which dimension to bisect.)
-bisect :: ( Box -> Bool )
-          -- ^ how to check whether a box contains solutions
-       -> ( Int, String )
-           -- ^ fall-back choice of dimension (and "why" we chose it)
-       -> Box
-       -> ( [ Box ], ( String, Double ) )
-bisect canHaveSols ( dim, why )
-  ( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) )
+bisect
+  :: forall n
+  .  ( KnownNat n, Representable Double ( ℝ n ) )
+  => ( Box n -> Bool )
+     -- ^ how to check whether a box contains solutions
+  -> ( Fin n, String )
+      -- ^ fallback choice of dimension (and "why" we chose it)
+  -> Box n
+  -> ( [ Box n ], ( String, Double ) )
+bisect canHaveSols ( fallbackDim, why ) box =
   -- We try to bisect along a dimension which eliminates zeros from one of the
   -- sub-regions.
   --
-  -- If this fails, we fall-back to the provided dimension choice.
-  = let t_mid = 0.5 * ( t_lo + t_hi )
-        s_mid = 0.5 * ( s_lo + s_hi )
-        l = 𝕀 ( ℝ2 t_lo  s_lo  ) ( ℝ2 t_mid s_hi  )
-        r = 𝕀 ( ℝ2 t_mid s_lo  ) ( ℝ2 t_hi  s_hi  )
-        d = 𝕀 ( ℝ2 t_lo  s_lo  ) ( ℝ2 t_hi  s_mid )
-        u = 𝕀 ( ℝ2 t_lo  s_mid ) ( ℝ2 t_hi  s_hi  )
-        l_skip = not $ canHaveSols l
-        r_skip = not $ canHaveSols r
-        d_skip = not $ canHaveSols d
-        u_skip = not $ canHaveSols u
+  -- If this fails, we fall back to the provided dimension choice.
+  case findFewestSols solsList of
+    Just ( Arg _nbSols ( i, mid, oks ) ) ->
+      ( oks, ("cd = " ++ show i, mid ) )
+    Nothing ->
+      case bisectInCoord box fallbackDim of
+        ( mid, (lo, hi) ) ->
+          ( [ lo, hi ], ( why, mid ) )
+  where
+    solsList =
+      [ Arg ( fromIntegral $ length oks ) ( i, mid, oks )
+      | i <- toList $ universe @n
+      , let (mid, (lo, hi)) = bisectInCoord box i
+            lo_ok = canHaveSols lo
+            hi_ok = canHaveSols hi
+            oks = [ lo | lo_ok ] ++ [ hi | hi_ok ]
+      ]
+{-# INLINEABLE bisect #-}
 
-        ( bisGuesses, whatBis )
-          | l_skip && r_skip
-          = ( [], ( "lr", t_mid ) )
-          | d_skip && u_skip
-          = ( [], ( "du", s_mid ) )
-          | l_skip
-          = ( [ r ], ( "r", t_mid ) )
-          | r_skip
-          = ( [ l ], ( "l", t_mid ) )
-          | d_skip
-          = ( [ u ], ( "u", s_mid ) )
-          | u_skip
-          = ( [ d ], ( "d", s_mid ) )
-          | otherwise
-          = let why' = case why of
-                  "" -> ""
-                  _ -> " (" ++ why ++ ")"
-            in case dim of
-                 0 -> ( [ l, r ], ( "t" ++ why', t_mid ) )
-                 1 -> ( [ d, u ], ( "s" ++ why', s_mid ) )
-                 _ -> error "bisect: fall-back bisection dimension should be either 0 or 1"
-         in ( bisGuesses, whatBis )
+-- | Find any element with the least argument.
+--
+-- NB: this function shortcuts as soon as it finds an element with argument 0.
+findFewestSols :: forall a. [ Arg Word a ] -> Maybe ( Arg Word a )
+findFewestSols [] = Nothing
+findFewestSols ( arg@( Arg nbSols _ ) : args )
+  | nbSols == 0
+  = Just arg
+  | otherwise
+  = Just $ go arg args
+  where
+    go :: Arg Word a -> [ Arg Word a ] -> Arg Word a
+    go bestSoFar [] = bestSoFar
+    go bestSoFar@( Arg bestNbSolsSoFar _ ) ( arg'@( Arg nbSols' _ ) : args' )
+      | nbSols' == 0
+      = arg'
+      | nbSols' < bestNbSolsSoFar
+      = go arg' args'
+      | otherwise
+      = go bestSoFar args'
+
+bisectInCoord
+  :: Representable Double ( ℝ n )
+  => Box n -> Fin n -> ( Double, ( Box n, Box n ) )
+bisectInCoord box i =
+  let 𝕀 z_lo z_hi = box `index` i
+      z_mid = 0.5 * ( z_lo + z_hi )
+  in ( z_mid
+     , ( set i ( 𝕀 z_lo  z_mid ) box
+       , set i ( 𝕀 z_mid z_hi  ) box ) )
+{-# INLINEABLE bisectInCoord #-}
 
 -- | Interval Newton method with Gauss–Seidel step.
 intervalGaussSeidel
-  :: ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) )
-     -- ^ equations
-  -> Preconditioner
-     -- ^ preconditioner to use for the Gauss–Seidel step
-  -> 𝕀ℝ 2
-  -> Writer ( [ 𝕀ℝ 2 ], [ 𝕀ℝ 2 ] ) [ 𝕀ℝ 2 ]
+  :: forall n d
+  .  BoxCt n d
+  => GaussSeidelOptions n d
+  -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
+      -- ^ equations
+  -> 𝕀ℝ n
+     -- ^ box
+  -> Writer ( [ 𝕀ℝ n ], [ 𝕀ℝ n ] ) [ 𝕀ℝ n ]
 intervalGaussSeidel
+  ( GaussSeidelOptions { gsPreconditioner = precondMeth, gsDims = projectDims } )
   eqs
-  precondMethod
-  ts@( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) )
-  | D12 _f ( T ee_f_t ) ( T ee_f_s )
-    <- eqs ts
-  , let f_t = take2 ee_f_t
-        f_s = take2 ee_f_s
-  , D12 ee_f_mid ( T _ee_f_t_mid ) ( T _ee_f_s_mid )
-    <- eqs ts_mid
-  , let f_mid = take2 ee_f_mid
+  box
+  | let boxMid = singleton $ midpoint box
+        df = eqs box
+        f' :: Vec n ( 𝕀ℝ n )
+        f' = fmap ( \ i -> projectDims $ df `monIndex` linearMonomial i ) ( universe @n )
+        f_mid = projectDims $ eqs boxMid `monIndex` zeroMonomial
 
   = let -- Interval Newton method: take one Gauss–Seidel step
-        -- for the equation f'(x) v = - f(x_mid),
-        -- where f = 𝛿E/𝛿s * dc/dt
-        ( a, b ) = precondition precondMethod
-                      ( midI f_t, midI f_s )
-                      ( f_t, f_s ) ( neg f_mid )
+        -- for the system of equations f'(x) ( x - x_mid ) = - f(x_mid).
+        ( a, b ) = precondition precondMeth
+                      ( fmap midpoint f' )
+                      f' ( negateBox $ singleton $ midpoint $ f_mid )
 
-        -- NB: we have to re-center around zero to take a Gauss–Seidel step.
-        gsGuesses = map ( first ( \ ts' -> unT $ T ts' ^+^ T ts_mid ) )
-                  $ gaussSeidelStep a b ( unT $ T ts ^-^ T ts_mid )
+        -- NB: we have to change coordinates, putting the midpoint of the box
+        -- at the origin, in order to take a Gauss–Seidel step.
+        gsGuesses = map ( first ( \ box' -> unT $ box' ^+^ T boxMid ) )
+                  $ gaussSeidelStep a b ( T box ^-^ T boxMid )
     in
          -- If the Gauss–Seidel step was a contraction, then the box
          -- contains a unique solution (by the Banach fixed point theorem).
@@ -425,57 +546,58 @@ intervalGaussSeidel
          -- Newton's method is guaranteed to converge to the unique solution.
          let !(done, todo) = bimap ( map fst ) ( map fst )
                            $ partition snd gsGuesses
-         in do tell ([], done)
+         in do tell ( [], done )
                return todo
   where
-    t_mid = 0.5 * ( t_lo + t_hi )
-    s_mid = 0.5 * ( s_lo + s_hi )
-    ts_mid = singleton ( ℝ2 t_mid s_mid )
-    neg ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) )
-      = let !( 𝕀 x'_lo x'_hi ) = negate $ 𝕀 x_lo x_hi
-            !( 𝕀 y'_lo y'_hi ) = negate $ 𝕀 y_lo y_hi
-        in 𝕀 ( ℝ2 x'_lo y'_lo ) ( ℝ2 x'_hi y'_hi )
-    take2 :: 𝕀ℝ 3 -> 𝕀ℝ 2
-    take2 ( 𝕀 ( ℝ3 _a_lo b_lo c_lo ) ( ℝ3 _a_hi b_hi c_hi ) )
-      = 𝕀 ( ℝ2 b_lo c_lo ) ( ℝ2 b_hi c_hi )
-      -- TODO: always using the last 2 coordinates, instead of any
-      -- pair of two coordinates.
+{-# INLINEABLE intervalGaussSeidel #-}
 
--- Take one interval Gauss–Seidel step for the equation \( A X = B \),
--- refining the initial guess box for \( X \) into up to four (disjoint) new boxes.
+-- | The midpoint of a box.
+midpoint :: Representable Double ( ℝ n ) => 𝕀ℝ n -> ℝ n
+midpoint box =
+  tabulate $ \ i ->
+    let 𝕀 z_lo z_hi = index box i
+        z_mid = 0.5 * ( z_lo + z_hi )
+    in z_mid
+{-# INLINEABLE midpoint #-}
+
+-- | Negate a box.
+negateBox :: Representable Double ( ℝ n ) => 𝕀ℝ n -> 𝕀ℝ n
+negateBox v = tabulate $ \ i -> negate ( index v i )
+{-# INLINEABLE negateBox #-}
+
+-- | Take one interval Gauss–Seidel step for the equation \( A X = B \),
+-- refining the initial guess box for \( X \) into up to \( 2^n \) (disjoint) new boxes.
 --
 -- The boolean indicates whether the Gauss–Seidel step was a contraction.
-gaussSeidelStep :: ( 𝕀ℝ 2, 𝕀ℝ 2 ) -- ^ columns of \( A \)
-                -> 𝕀ℝ 2           -- ^ \( B \)
-                -> 𝕀ℝ 2           -- ^ initial box \( X \)
-                -> [ ( 𝕀ℝ 2, Bool ) ]
 gaussSeidelStep
-  ( 𝕀 ( ℝ2 a11_lo a21_lo ) ( ℝ2 a11_hi a21_hi )
-  , 𝕀 ( ℝ2 a12_lo a22_lo ) ( ℝ2 a12_hi a22_hi ) )
-  ( 𝕀 ( ℝ2 b1_lo b2_lo ) ( ℝ2 b1_hi b2_hi ) )
-  ( 𝕀 ( ℝ2 x1_lo x2_lo ) ( ℝ2 x1_hi x2_hi ) )
-    = let !a11 = 𝕀 a11_lo a11_hi
-          !a12 = 𝕀 a12_lo a12_hi
-          !a21 = 𝕀 a21_lo a21_hi
-          !a22 = 𝕀 a22_lo a22_hi
-          !b1  = 𝕀  b1_lo  b1_hi
-          !b2  = 𝕀  b2_lo  b2_hi
-          !x1  = 𝕀  x1_lo  x1_hi
-          !x2  = 𝕀  x2_lo  x2_hi
-      in nub $ do
-
-        -- x1' = ( ( b1 - a12 * x2 ) / a11 ) ∩ x1
-        x1'_0 <- ( b1 - a12 * x2 ) `extendedDivide` a11
-        ( x1'@( 𝕀 x1'_lo x1'_hi ), sub_x1 ) <- x1'_0 `intersect` x1
-        -- x2' = ( ( b2 - a21 * x1' ) / a22 ) ∩ x2
-        x2'_0 <- ( b2 - a21 * x1' ) `extendedDivide` a22
-        ( 𝕀 x2'_lo x2'_hi, sub_x2 ) <- x2'_0 `intersect` x2
-
-        return ( ( 𝕀 ( ℝ2 x1'_lo x2'_lo) ( ℝ2 x1'_hi x2'_hi ) )
-               , sub_x1 && sub_x2 )
+  :: forall n
+  .  ( Representable Double ( ℝ n ), Eq ( ℝ n ) )
+  => Vec n ( 𝕀ℝ n ) -- ^ columns of \( A \)
+  -> 𝕀ℝ n           -- ^ \( B \)
+  -> T ( 𝕀ℝ n )     -- ^ initial box \( X \)
+  -> [ ( T ( 𝕀ℝ n ), Bool ) ]
+gaussSeidelStep as b ( T x0 ) = coerce $ nub $
+  forEachDim @n ( x0, True ) $ \ i ( x, contraction ) -> do
+    -- x_i' = ( b_i - sum { j /= i } a_ij * x_j ) / a_ii
+    x_i'0 <- ( b `index` i - sum [ ( as ! j ) `index` i * x `index` j | j <- toList ( universe @n ), j /= i ] )
+          `extendedDivide` ( ( as ! i ) `index` i )
+    ( x_i', sub_i ) <- x_i'0 `intersect` ( x `index` i )
+    return $ ( set i x_i' x, sub_i && contraction )
 -- TODO: try implementing the complete interval union Gauss–Seidel algorithm.
 -- See "Algorithm 2" in
 -- "Using interval unions to solve linear systems of equations with uncertainties"
+{-# INLINEABLE gaussSeidelStep #-}
+
+-- | Run an effectful computation several times in sequence, piping its output
+-- into the next input, once for each coordinate dimension.
+forEachDim :: forall n a m. ( KnownNat n, Monad m ) => a -> ( Fin n -> a -> m a ) -> m a
+forEachDim a0 f = go ( toList $ universe @n ) a0
+  where
+    go [] a = return a
+    go ( i : is ) a = do
+      a' <- f i a
+      go is a'
+{-# INLINEABLE forEachDim #-}
 
 -- | Compute the intersection of two intervals.
 --
@@ -503,9 +625,10 @@ data Preconditioner
 -- See also
 -- "Presentation of a highly tuned multithreaded interval solver for underdetermined and well-determined nonlinear systems"
 makeBox1Consistent
-  :: ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) )
+  :: BoxCt n d
+  => ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
   -> Double -> Double
-  -> Box -> [ Box ]
+  -> Box n -> [ Box n ]
 makeBox1Consistent eqs minWidth epsEq x =
   ( `State.evalState` False ) $
     pipeFunctionsWhileTrue ( allNarrowingOperators epsEq minWidth eqs ) x
@@ -513,24 +636,26 @@ makeBox1Consistent eqs minWidth epsEq x =
 -- | An implementation of "bound-consistency" from the paper
 -- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
 makeBox2Consistent
-  :: ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) )
+  :: forall n d
+  .  BoxCt n d
+  => ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
   -> Double -> Double -> Double
-  -> Box -> Box
+  -> Box n -> Box n
 makeBox2Consistent eqs minWidth epsEq lambdaMin x0 = ( `State.evalState` False ) $ doLoop 0.25 x0
   where
-    doLoop :: Double -> Box -> State Bool Box
+    doLoop :: Double -> Box n -> State Bool ( Box n )
     doLoop lambda x = do
-      x'  <- boundConsistency get_t set_t lambda x
-      x'' <- boundConsistency get_s set_s lambda x'
+      x'' <- forEachDim @n x $ \ i ->
+        boundConsistency ( `index` i ) ( set i ) lambda
       modified <- State.get
       let lambda' = if modified then lambda else 0.5 * lambda
       if lambda' < lambdaMin
       then return x''
       else do { State.put False ; doLoop lambda' x'' }
 
-    boundConsistency :: ( Box -> 𝕀 Double )
-                     -> ( 𝕀 Double -> Box -> Box )
-                     -> Double -> Box -> State Bool Box
+    boundConsistency :: ( Box n -> 𝕀 Double )
+                     -> ( 𝕀 Double -> Box n -> Box n )
+                     -> Double -> Box n -> State Bool ( Box n )
     boundConsistency getter setter lambda box = do
       let x@( 𝕀 x_inf x_sup ) = getter box
           c1 = ( 1 - lambda ) * x_inf + lambda * x_sup
@@ -548,62 +673,63 @@ makeBox2Consistent eqs minWidth epsEq lambdaMin x0 = ( `State.evalState` False )
         State.put True
       return $ setter x' box
 
-precondition :: Preconditioner
-             -> ( 𝕀ℝ 2, 𝕀ℝ 2 )
-             -> ( 𝕀ℝ 2, 𝕀ℝ 2 )
-             -> 𝕀ℝ 2
-             -> ( ( 𝕀ℝ 2, 𝕀ℝ 2 ), 𝕀ℝ 2 )
-precondition meth jac_mid a@( a1, a2 ) b =
+-- | Pre-condition the system \( AX = B \).
+precondition
+  :: forall n
+  .  ( KnownNat n, Representable Double ( ℝ n ) )
+  => Preconditioner -- ^ pre-conditioning method to use
+  -> Vec n ( ℝ n )  -- ^ entry-wise midpoint matrix of the interval Jacobian matrix
+  -> Vec n ( 𝕀ℝ n ) -- ^ columns of \( A \)
+  -> 𝕀ℝ n           -- ^ \( B \)
+  -> ( Vec n ( 𝕀ℝ n ), 𝕀ℝ n )
+precondition meth jac_mid as b =
   case meth of
     NoPreconditioning
-      -> ( a, b )
+      -> ( as, b )
     InverseMidJacobian
-      | ( 𝕀 ( ℝ2 a11_lo a21_lo ) ( ℝ2 a11_hi a21_hi )
-        , 𝕀 ( ℝ2 a12_lo a22_lo ) ( ℝ2 a12_hi a22_hi ) ) <- jac_mid
-      , let !a11 = 0.5 * ( a11_lo + a11_hi )
-            !a12 = 0.5 * ( a12_lo + a12_hi )
-            !a21 = 0.5 * ( a21_lo + a21_hi )
-            !a22 = 0.5 * ( a22_lo + a22_hi )
-            !d = a11 * a22 - a12 * a21
-      , not ( nearZero d )
-      , let !precond = ( ℝ2 a22 -a12, ℝ2 -a21 a11 )
-            !inv = recip d
-            f x = scale inv $ matMulVec precond x
-      -> ( ( f a1, f a2 ), f b )
+      | let mat = toEigen jac_mid
+            det = Eigen.determinant mat
+      , not $ nearZero det
+      , let precond = Eigen.inverse mat
+            doPrecond = matMulVec ( fromEigen precond )
+      -> ( fmap doPrecond as, doPrecond b )
       | otherwise
-      -> ( a, b )
-
-scale :: Double -> 𝕀ℝ 2 -> 𝕀ℝ 2
-scale s ( 𝕀 ( ℝ2 a1_lo a2_lo ) ( ℝ2 a1_hi a2_hi ) )
-  = 𝕀 ( ℝ2 b1_lo b2_lo ) ( ℝ2 b1_hi b2_hi )
+      -> ( as, b )
   where
-    𝕀 b1_lo b1_hi = scaleInterval s ( 𝕀 a1_lo a1_hi )
-    𝕀 b2_lo b2_hi = scaleInterval s ( 𝕀 a2_lo a2_hi )
+    toEigen :: Vec n ( ℝ n ) -> Eigen.Matrix n n Double
+    toEigen cols =
+      Eigen.generate $ \ r c ->
+        ( cols ! Fin ( fromIntegral c + 1 ) ) `index` ( Fin ( fromIntegral r + 1 ) )
 
-matMulVec :: ( ℝ 2, ℝ 2 ) -> 𝕀ℝ 2 -> 𝕀ℝ 2
-matMulVec ( ℝ2 a11 a21, ℝ2 a12 a22 ) ( 𝕀 ( ℝ2 u_lo v_lo ) ( ℝ2 u_hi v_hi ) )
-  = 𝕀 ( ℝ2 u'_lo v'_lo ) ( ℝ2 u'_hi v'_hi )
-  where
-    u = 𝕀 u_lo u_hi
-    v = 𝕀 v_lo v_hi
-    𝕀 u'_lo u'_hi = scaleInterval a11 u + scaleInterval a12 v
-    𝕀 v'_lo v'_hi = scaleInterval a21 u + scaleInterval a22 v
+    fromEigen :: Eigen.Matrix n n Double -> Vec n ( ℝ n )
+    fromEigen mat =
+      fmap
+        ( \ ( Fin c ) ->
+          tabulate $ \ ( Fin r ) ->
+            Eigen.unsafeCoeff ( fromIntegral r - 1 ) ( fromIntegral c - 1 ) mat
+        )
+        ( universe @n )
 
-midI :: 𝕀ℝ 2 -> 𝕀ℝ 2
-midI ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) =
-  singleton $ ℝ2 ( 0.5 * ( x_lo + x_hi ) ) ( 0.5 * ( y_lo + y_hi ) )
---width :: 𝕀ℝ 1 -> Double
---width ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) = hi - lo
---normI :: 𝕀ℝ 1 -> Double
---normI ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) = sqrt $ sup $ ( 𝕀 lo hi ) Ring.^ 2
---normVI :: 𝕀ℝ 2 -> Double
---normVI ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) =
---  sqrt $ sup $ ( 𝕀 x_lo x_hi ) Ring.^ 2 + ( 𝕀 y_lo y_hi ) Ring.^ 2
-normVI3 :: 𝕀ℝ 3 -> Double
-normVI3 ( 𝕀 ( ℝ3 lo x_lo y_lo ) ( ℝ3 hi x_hi y_hi ) ) =
-  sqrt $ max ( abs lo   ) ( abs hi   ) Ring.^ 2
-       + max ( abs x_lo ) ( abs x_hi ) Ring.^ 2
-       + max ( abs y_lo ) ( abs y_hi ) Ring.^ 2
+-- | Matrix multiplication \( A v \).
+matMulVec
+  :: forall n m
+  .  ( Representable Double ( ℝ n ), Representable Double ( ℝ m ) )
+  => Vec m ( ℝ n ) -- ^ columns of the matrix \( A )
+  -> 𝕀ℝ m          -- ^ vector \( v \)
+  -> 𝕀ℝ n
+matMulVec as v = tabulate $ \ r ->
+  sum [ scaleInterval ( ( as ! c ) `index` r ) ( index v c )
+      | c <- toList ( universe @m )
+      ]
+{-# INLINEABLE matMulVec #-}
+
+normVI :: ( Applicative ( Vec d ), Representable Double ( ℝ d ) ) => 𝕀ℝ d -> Double
+normVI ( 𝕀 los his ) =
+  sqrt $ sum ( nm1 <$> coordinates los <*> coordinates his )
+    where
+      nm1 :: Double -> Double -> Double
+      nm1 lo hi = max ( abs lo ) ( abs hi ) Ring.^ 2
+{-# INLINEABLE normVI #-}
 
 -- Use the univariate interval Newton method to narrow from the left
 -- a candidate interval.
@@ -704,43 +830,36 @@ pipeFunctionsWhileTrue fns = go fns
       concat <$> traverse ( go fs ) xs
 
 allNarrowingOperators
-  :: Double
+  :: forall n d
+  .  BoxCt n d
+  => Double
   -> Double
-  -> ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) )
-  -> [ Box -> State Bool [ Box ] ]
+  -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
+  -> [ Box n -> State Bool [ Box n ] ]
 allNarrowingOperators eps_eq eps_bis eqs =
   [ \ cand ->
     let getter = ( `index` coordIndex )
         setter = set coordIndex
-        newCands = map ( `setter` cand ) $ narrowFn eps_eq eps_bis ( \ x0 -> ff' $ setter x0 cand ) ( getter cand )
+        newCands = map ( `setter` cand ) $ narrowFn eps_eq eps_bis ( \ x0 -> fn $ setter x0 cand ) ( getter cand )
     in do
       when ( ( width ( getter cand ) - sum ( map ( width . getter ) newCands ) ) > eps_eq ) $
         -- TODO: do a check with the relative reduction in size?
         State.put True
       return newCands
   | narrowFn <- [ leftNarrow, rightNarrow ]
-  , ( coordIndex, ff' ) <-
-      [ ( Fin 1, ff'_t i ) | i <- [ Fin 1, Fin 2, Fin 3 ] ]
-        ++
-      [ ( Fin 2, ff'_s i ) | i <- [ Fin 1, Fin 2, Fin 3 ] ]
+  , ( coordIndex, fn ) <-
+      [ ( i, ff' i d )
+      | i <- toList $ universe @n
+      , d <- toList $ universe @d
       ]
+  ]
   where
 
-    ff'_t i ts =
-      let D12 { _D12_v = f
-              , _D12_dx = T f_t
-              } = eqs ts
-      in ( f `index` i, f_t `index` i )
-    ff'_s i ts =
-      let D12 { _D12_v = f
-              , _D12_dy = T f_t
-              } = eqs ts
-      in ( f `index` i, f_t `index` i )
-
-get_t, get_s :: Box -> 𝕀 Double
-get_t = ( `index` ( Fin 1 ) )
-get_s = ( `index` ( Fin 2 ) )
-
-set_t, set_s :: 𝕀 Double -> Box -> Box
-set_t = set ( Fin 1 )
-set_s = set ( Fin 2 )
+    ff' :: Fin n -> Fin d -> 𝕀ℝ n -> ( 𝕀 Double, 𝕀 Double )
+    ff' i d ts =
+      let df = eqs ts
+          f, f' :: Box d
+          f = df `monIndex` zeroMonomial
+          f' = df `monIndex` linearMonomial i
+      in ( f `index` d, f' `index` d )
+{-# INLINEABLE allNarrowingOperators #-}
diff --git a/brush-strokes/src/test/Main.hs b/brush-strokes/src/test/Main.hs
index b85fc22..734ef35 100644
--- a/brush-strokes/src/test/Main.hs
+++ b/brush-strokes/src/test/Main.hs
@@ -25,7 +25,7 @@ import Math.Linear
 import Math.Module
 import Math.Monomial
   ( multiSubsetSum, multiSubsetsSum
-  , MonomialBasis ( monTabulate, monIndex )
+  , MonomialBasisQ ( monTabulateQ, monIndexQ )
   )
 import Math.Ring
 
@@ -65,7 +65,7 @@ main =
     Tasty.testGroup "brush-strokes property tests"
     [ Tasty.testGroup "Automatic differentiation"
       [ Tasty.testGroup "Monomial basis"
-        [ testProperty "Round trip D33" testMonomialBasisD33
+        [ testProperty "Round trip D33" testMonomialBasisQD33
         ]
       , Tasty.testGroup "Monomials"
         [ Tasty.testGroup "multiSubsetSum"
@@ -172,9 +172,9 @@ testRoundTrip g roundTrip = do
     .$ ("value", d )
     .$ ("round tripped", roundTrip d )
 
-testMonomialBasisD33 :: Falsify.Property ()
-testMonomialBasisD33 =
-  testRoundTrip genD33 \ d -> $$( monTabulate \ mon -> monIndex [|| d ||] mon )
+testMonomialBasisQD33 :: Falsify.Property ()
+testMonomialBasisQD33 =
+  testRoundTrip genD33 \ d -> $$( monTabulateQ \ mon -> monIndexQ [|| d ||] mon )
   where
     genD33 :: Falsify.Gen ( D3𝔸3 Double )
     genD33 =