diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs
index 7bd4867..82dfd88 100644
--- a/brush-strokes/src/cusps/bench/Main.hs
+++ b/brush-strokes/src/cusps/bench/Main.hs
@@ -35,7 +35,7 @@ import GHC.Exts
 import GHC.Generics
   ( Generic )
 import GHC.TypeNats
-  ( type (-) )
+  ( Nat, type (-) )
 import Numeric
   ( showFFloat )
 
@@ -126,8 +126,8 @@ benchTestCase ( TestCase { testName, testBrushStroke, testCuspOptions, testStart
           IntMap.toList $
             findCuspsIn testCuspOptions testStrokeFnI $
               IntMap.fromList
-                [ ( i, [ box ] )
-                | ( i, box ) <- testStartBoxes
+                [ ( i, boxes )
+                | ( i, boxes ) <- testStartBoxes
                 ]
   rnf dunno `seq` rnf sols `seq` return ()
   after <- getMonotonicTime
@@ -199,8 +199,8 @@ data TestCase =
   TestCase
   { testName        :: String
   , testBrushStroke :: BrushStroke
-  , testCuspOptions :: RootIsolationOptions 2 3
-  , testStartBoxes  :: [ ( Int, Box 2 ) ]
+  , testCuspOptions :: RootIsolationOptions N 3
+  , testStartBoxes  :: [ ( Int, [ Box 2 ] ) ]
   }
 
 brushStrokeFunctions
@@ -488,7 +488,7 @@ showD float = showFFloat (Just 6) float ""
 
 --------------------------------------------------------------------------------
 
-ellipseTestCase :: RootIsolationOptions 2 3 -> String -> ( Double, Double ) -> Double -> [ ( Int, Box 2 ) ] -> TestCase
+ellipseTestCase :: RootIsolationOptions N 3 -> String -> ( Double, Double ) -> Double -> [ ( Int, [ Box 2 ] ) ] -> TestCase
 ellipseTestCase opts str k0k1 rot startBoxes =
   TestCase
     { testName = "ellipse (" ++ str ++ ")"
@@ -607,11 +607,15 @@ getStrokeFunctions brush sp0 crv =
         brush @3 @𝕀 proxy# singleton )
 {-# INLINEABLE getStrokeFunctions #-}
 
-defaultStartBoxes :: [ Int ] -> [ ( Int, Box 2 ) ]
+defaultStartBoxes :: [ Int ] -> [ ( Int, [ Box 2 ] ) ]
 defaultStartBoxes is =
-  [ ( i, 𝕀 ( ℝ2 zero zero ) ( ℝ2 one one ) )
-  | i <- is
-  ]
+  [ ( i, [ 𝕀 ( ℝ2 zero zero ) ( ℝ2 one one ) ] ) | i <- is ]
+--  [ ( i, [ 𝕀 ( ℝ2 t s ) ( ℝ2 ( t + 0.099999 ) ( s + 0.099999 ) )
+--         | t <- [ 0.00001, 0.1 .. 0.9 ]
+--         , s <- [ 0.00001, 0.1 .. 0.9 ]
+--         ] )
+--  | i <- is
+--  ]
 
 getR1 (ℝ1 u) = u
 
diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke.hs b/brush-strokes/src/lib/Math/Bezier/Stroke.hs
index 60e6e64..cba086e 100644
--- a/brush-strokes/src/lib/Math/Bezier/Stroke.hs
+++ b/brush-strokes/src/lib/Math/Bezier/Stroke.hs
@@ -18,6 +18,9 @@ module Math.Bezier.Stroke
 
   -- ** Cusp finding
   , findCusps, findCuspsIn
+
+  -- TODO: hack for switching between 2 and 3 dim formulations of cusp-finding
+  , N
   )
   where
 
@@ -243,7 +246,7 @@ computeStrokeOutline ::
 
      )
   => RootSolvingAlgorithm
-  -> Maybe ( RootIsolationOptions 2 3 )
+  -> Maybe ( RootIsolationOptions N 3 )
   -> FitParameters
   -> ( ptData -> usedParams )
   -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing
@@ -522,7 +525,7 @@ outlineFunction
      , Show ptData, Show brushParams
      )
   => RootSolvingAlgorithm
-  -> Maybe ( RootIsolationOptions 2 3 )
+  -> Maybe ( RootIsolationOptions N 3 )
   -> ( ptData -> usedParams )
   -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing
   -> ( forall {t} k (i :: t)
@@ -1147,9 +1150,9 @@ 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 2 )
+           -> ( Int, Box N )
            -> Cusp
-cuspCoords eqs ( i, 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) )
+cuspCoords eqs ( i, box )
   | StrokeDatum
     { dpath
     , dstroke = D22 { _D22_v = stroke } }
@@ -1160,17 +1163,21 @@ cuspCoords eqs ( i, 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) )
   , cuspStrokeCoords = coerce stroke
   }
   where
+    𝕀 t_lo t_hi = box `index` Fin 1
+    𝕀 s_lo s_hi = box `index` Fin 2
     t_mid = 0.5 * ( t_lo + t_hi )
     s_mid = 0.5 * ( s_lo + s_hi )
 
+type N = 2
+
 -- | Find cusps in the envelope for values of the parameters in
 -- \( 0 \leqslant t, s \leqslant 1 \), using interval arithmetic.
 --
 -- See 'isolateRootsIn' for details of the algorithms used.
 findCusps
-  :: RootIsolationOptions 2 3
+  :: RootIsolationOptions N 3
   -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-  -> IntMap ( [ ( Box 2, RootIsolationTree ( Box 2 ) ) ], ( [ Box 2 ], [ Box 2 ] ) )
+  -> IntMap ( [ ( Box N, RootIsolationTree ( Box N ) ) ], ( [ Box N ], [ Box N ] ) )
 findCusps opts boxStrokeData =
   findCuspsIn opts boxStrokeData $
     IntMap.fromList
@@ -1187,10 +1194,10 @@ findCusps opts boxStrokeData =
 -- | Like 'findCusps', but parametrised over the initial boxes for the
 -- root isolation method.
 findCuspsIn
-  :: RootIsolationOptions 2 3
+  :: RootIsolationOptions N 3
   -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
   -> IntMap [ Box 2 ]
-  -> IntMap ( [ ( Box 2, RootIsolationTree ( Box 2 ) ) ], ( [ Box 2 ], [ Box 2 ] ) )
+  -> IntMap ( [ ( Box N, RootIsolationTree ( Box N ) ) ], ( [ Box N ], [ Box N ] ) )
 findCuspsIn opts boxStrokeData initBoxes =
   IntMap.mapWithKey ( \ i -> foldMap ( isolateRootsIn opts ( eqnPiece i ) ) ) initBoxes
     where
@@ -1214,13 +1221,42 @@ findCuspsIn opts boxStrokeData initBoxes =
 
 {-
 findCuspsIn
-  :: RootIsolationOptions 3 3
+  :: RootIsolationOptions N 3
   -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-  -> IntMap [ Box 3 ]
+  -> IntMap [ Box 2 ]
   -> IntMap ( [ ( Box 3, RootIsolationTree ( Box 3 ) ) ], ( [ Box 3 ], [ Box 3 ] ) )
 findCuspsIn opts boxStrokeData initBoxes =
-  IntMap.mapWithKey ( \ i -> foldMap ( isolateRootsIn opts ( eqnPiece i ) ) ) initBoxes
+  IntMap.mapWithKey ( \ i boxes -> foldMap ( isolateRootsIn opts ( eqnPiece i ) ) $ concatMap ( mkInitBox i ) boxes ) initBoxes
     where
+      mkInitBox :: Int -> Box 2 -> [ Box 3 ]
+      mkInitBox i ( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) =
+        let t = 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi )
+            s = 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi )
+            StrokeDatum
+              { dstroke =
+                D32
+                  { _D32_dx = T ( 𝕀 ( ℝ2 cx_t_lo cy_t_lo ) ( ℝ2 cx_t_hi cy_t_hi ) )
+                  , _D32_dy = T ( 𝕀 ( ℝ2 cx_s_lo cy_s_lo ) ( ℝ2 cx_s_hi cy_s_hi ) )}
+              , ee =
+                D22
+                  { _D22_dx = T ( 𝕀 ( ℝ1 ee_t_lo ) ( ℝ1 ee_t_hi ) )
+                  , _D22_dy = T ( 𝕀 ( ℝ1 ee_s_lo ) ( ℝ1 ee_s_hi ) ) }
+              } = ( boxStrokeData t `Seq.index` i ) s
+        -- λ = ∂E/∂t / ∂E/∂s
+            λ1 = 𝕀 ee_t_lo ee_t_hi `extendedDivide` 𝕀 ee_s_lo ee_s_hi
+        -- λ = ∂c/∂t / ∂c/∂s
+            λ2 = 𝕀 cx_t_lo cx_t_hi `extendedDivide` 𝕀 cx_s_lo cx_s_hi
+            λ3 = 𝕀 cy_t_lo cy_t_hi `extendedDivide` 𝕀 cy_s_lo cy_s_hi
+            λ = [ 𝕀 ( recip -0 ) ( recip 0 ) ]
+                `intersectMany` λ1
+                `intersectMany` λ2
+                `intersectMany` λ3
+        in
+          let boxes = [ 𝕀 ( ℝ3 t_lo s_lo λ_lo ) ( ℝ3 t_hi s_hi λ_hi )
+                      | 𝕀 λ_lo' λ_hi' <- λ
+                      , let λ_lo = max λ_lo' ( min ( λ_hi' - 10 ) -100 )
+                            λ_hi = min λ_hi' ( max ( λ_lo' + 10 )  100 ) ]
+          in boxes
       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 )
@@ -1253,5 +1289,20 @@ findCuspsIn opts boxStrokeData initBoxes =
                ( 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 ) )
 
+intersectMany :: [𝕀 Double] -> [𝕀 Double] -> [𝕀 Double]
+intersectMany _ [] = []
+intersectMany is (j : js) = intersectOne is j ++ intersectMany is js
 
+intersectOne :: [ 𝕀 Double ] -> 𝕀 Double -> [ 𝕀 Double ]
+intersectOne is i = concatMap ( intersect i ) is
+
+intersect :: 𝕀 Double -> 𝕀 Double -> [ 𝕀 Double ]
+intersect ( 𝕀 lo1 hi1 ) ( 𝕀 lo2 hi2 )
+  | lo > hi
+  = [ ]
+  | otherwise
+  = [ 𝕀 lo hi ]
+  where
+    lo = max lo1 lo2
+    hi = min hi1 hi2
 -}
\ 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 0ddbe84..39d7c0d 100644
--- a/brush-strokes/src/lib/Math/Interval.hs
+++ b/brush-strokes/src/lib/Math/Interval.hs
@@ -18,6 +18,8 @@ module Math.Interval
 
 -- base
 import Prelude hiding ( Num(..), Fractional(..) )
+import Data.List
+  ( nub )
 
 -- acts
 import Data.Act
@@ -78,7 +80,7 @@ instance Torsor ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where
 -- Extended division
 
 extendedDivide :: 𝕀 Double -> 𝕀 Double -> [ 𝕀 Double ]
-extendedDivide x y = map ( x * ) ( extendedRecip y )
+extendedDivide x y = nub $ map ( x * ) ( extendedRecip y )
 
 extendedRecip :: 𝕀 Double -> [ 𝕀 Double ]
 extendedRecip x@( 𝕀 lo hi )
diff --git a/brush-strokes/src/lib/Math/Root/Isolation.hs b/brush-strokes/src/lib/Math/Root/Isolation.hs
index 5d18c44..7613a52 100644
--- a/brush-strokes/src/lib/Math/Root/Isolation.hs
+++ b/brush-strokes/src/lib/Math/Root/Isolation.hs
@@ -41,7 +41,7 @@ import Data.Kind
 import Data.Foldable
  ( toList )
 import Data.List
-  ( nub, partition, sort )
+  ( partition, sort )
 import Data.List.NonEmpty
   ( NonEmpty )
 import qualified Data.List.NonEmpty as NE
@@ -81,7 +81,7 @@ import Math.Linear
 import Math.Module
   ( Module(..) )
 import Math.Monomial
-  ( MonomialBasis(..), Deg, Vars
+  ( MonomialBasis(..)
   , linearMonomial, zeroMonomial
   )
 import qualified Math.Ring as Ring
@@ -202,7 +202,7 @@ data Box2Options =
     , box2LambdaMin :: !Double
     }
 
-defaultRootIsolationOptions :: RootIsolationOptions 2 3
+defaultRootIsolationOptions :: BoxCt n d => RootIsolationOptions n d
 defaultRootIsolationOptions =
   RootIsolationOptions
     { minWidth
@@ -251,9 +251,8 @@ 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 )
+    , gsDims = \ ( 𝕀 ( ℝ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
@@ -576,7 +575,7 @@ gaussSeidelStep
   -> 𝕀ℝ n           -- ^ \( B \)
   -> T ( 𝕀ℝ n )     -- ^ initial box \( X \)
   -> [ ( T ( 𝕀ℝ n ), Bool ) ]
-gaussSeidelStep as b ( T x0 ) = coerce $ nub $
+gaussSeidelStep as b ( T x0 ) = coerce $
   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 ] )