diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs
index 9e4366d..c408119 100644
--- a/brush-strokes/src/cusps/bench/Main.hs
+++ b/brush-strokes/src/cusps/bench/Main.hs
@@ -185,8 +185,7 @@ benchCases =
   , narrowAbs <- [ 5e-2 ]
   , let opts =
           RootIsolationOptions
-            { minWidth
-            , cuspFindingAlgorithms = defaultRootIsolationAlgorithms minWidth narrowAbs
+            { rootIsolationAlgorithms = defaultRootIsolationAlgorithms minWidth narrowAbs
             }
   ]
 
diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke.hs b/brush-strokes/src/lib/Math/Bezier/Stroke.hs
index 7c64a2e..737a371 100644
--- a/brush-strokes/src/lib/Math/Bezier/Stroke.hs
+++ b/brush-strokes/src/lib/Math/Bezier/Stroke.hs
@@ -229,6 +229,8 @@ data OutlineInfo =
     { outlineFn :: OutlineFn
     , outlineDefiniteCusps, outlinePotentialCusps :: [ Cusp ] }
 
+type N = 2
+
 computeStrokeOutline ::
   forall ( clo :: SplineType ) usedParams brushParams crvData ptData s
   .  ( KnownSplineType clo
diff --git a/brush-strokes/src/lib/Math/Linear.hs b/brush-strokes/src/lib/Math/Linear.hs
index 3d500f3..b246be4 100644
--- a/brush-strokes/src/lib/Math/Linear.hs
+++ b/brush-strokes/src/lib/Math/Linear.hs
@@ -10,7 +10,8 @@ module Math.Linear
 
   -- * Points and vectors (second version)
   , ℝ(..), T(.., V2, V3, V4)
-  , Fin(..), MFin(..), universe, coordinates
+  , Fin(..), MFin(..)
+  , universe, coordinates, choose
   , RepDim, RepresentableQ(..)
   , Representable(..), set, injection, projection
   , Vec(..), (!), find, zipIndices
@@ -34,7 +35,9 @@ import GHC.Generics
 import GHC.Stack
   ( HasCallStack )
 import GHC.TypeNats
-  ( Nat, KnownNat, natVal' )
+  ( Nat, KnownNat, type (<=)
+  , natVal'
+  )
 
 -- acts
 import Data.Act
@@ -136,6 +139,27 @@ coordinates :: forall r u. ( Representable r u ) => u -> Vec ( RepDim u ) r
 coordinates u = fmap ( index u ) $ universe @( RepDim u )
 {-# INLINEABLE coordinates #-}
 
+-- | Binomial coefficient: choose all subsets of size @k@ of the given set
+-- of size @n@.
+choose
+  :: forall n k
+   . ( KnownNat n, KnownNat k
+     , 1 <= n, 1 <= k, k <= n
+     )
+  => [ Vec k ( Fin n ) ]
+choose = coerce $ go ( fromIntegral $ natVal' @n proxy# )
+                     ( fromIntegral $ natVal' @k proxy# )
+  where
+    go :: Word -> Word -> [ [ Word ] ]
+    go n k
+      | k == 1
+      = [ [ i ] | i <- [ 1 .. n ] ]
+      | n == k
+      = [ [ 1 .. n ] ]
+    go n k = go ( n - 1 ) k
+          ++ ( map ( ++ [ n ] ) $ go ( n - 1 ) ( k - 1 ) )
+{-# INLINEABLE choose #-}
+
 infixl 9 !
 (!) :: forall l a. HasCallStack => Vec l a -> Fin l -> a
 ( Vec l ) ! Fin i = l !! ( fromIntegral i - 1 )
diff --git a/brush-strokes/src/lib/Math/Root/Isolation.hs b/brush-strokes/src/lib/Math/Root/Isolation.hs
index 5a8df0f..13e9114 100644
--- a/brush-strokes/src/lib/Math/Root/Isolation.hs
+++ b/brush-strokes/src/lib/Math/Root/Isolation.hs
@@ -24,10 +24,6 @@ module Math.Root.Isolation
     -- ** Trees recording search space of root isolation algorithms
   , RootIsolationTree(..), showRootIsolationTree
   , RootIsolationStep(..)
-
-    -- * Hack for changing between 2 and 3 d formulations
-    -- for my personal testing
-  , N
   )
   where
 
@@ -52,13 +48,19 @@ import Data.List.NonEmpty
   ( NonEmpty )
 import qualified Data.List.NonEmpty as NE
   ( NonEmpty(..), cons, filter, fromList, last, singleton, sort )
+import Data.Proxy
+  ( Proxy(..) )
 import Data.Semigroup
   ( Arg(..), Dual(..) )
+import Data.Type.Ord
+  ( OrderingI(..) )
 import Numeric
   ( showFFloat )
 import GHC.TypeNats
   ( Nat, KnownNat
-  , type (<=) )
+  , type (<=)
+  , cmpNat
+  )
 
 -- containers
 import Data.Tree
@@ -88,7 +90,7 @@ import Math.Linear
 import Math.Module
   ( Module(..) )
 import Math.Monomial
-  ( MonomialBasis(..)
+  ( MonomialBasis(..), Deg, Vars
   , linearMonomial, zeroMonomial
   )
 import qualified Math.Ring as Ring
@@ -137,10 +139,18 @@ showArea area = "(area " ++ showFFloat (Just 6) area "" ++ ")"
 
 type Box n = 𝕀ℝ n
 type BoxHistory n = [ NE.NonEmpty ( RootIsolationStep, Box n ) ]
-type N = 2
-type BoxCt n d = ( n ~ N, d ~ 3 )
-{-
-  ( Show ( 𝕀ℝ n ), Show ( ℝ n )
+
+-- | Dimension constraints for root isolation in a system of equations:
+--
+--   - @n@: number of variables
+--   - @d@: number of equations
+--
+-- NB: we require n <= d (no support for under-constrained systems).
+type BoxCt n d =
+  ( KnownNat n, KnownNat d
+  , 1 <= n, 1 <= d, n <= d
+
+  , Show ( 𝕀ℝ n ), Show ( ℝ n )
   , Eq ( ℝ n )
   , Representable Double ( ℝ n )
   , MonomialBasis ( D 1 ( ℝ n ) )
@@ -153,7 +163,6 @@ type BoxCt n d = ( n ~ N, d ~ 3 )
   , Module Double ( T ( ℝ d ) )
   , Representable Double ( ℝ d )
   )
--}
 
 -- | Options for the root isolation methods in 'isolateRootsIn'.
 type RootIsolationOptions :: Nat -> Nat -> Type
@@ -235,6 +244,7 @@ defaultRootIsolationOptions =
   where
     minWidth = 1e-5
     ε_eq = 5e-3
+{-# INLINEABLE defaultRootIsolationOptions #-}
 
 defaultRootIsolationAlgorithms
   :: forall n d
@@ -284,11 +294,17 @@ defaultRootIsolationAlgorithms minWidth ε_eq box history
       GaussSeidelOptions
         { gsPreconditioner = InverseMidJacobian
         , gsPickEqs =
-            \ ( 𝕀 ( ℝ3 a_lo b_lo c_lo ) ( ℝ3 a_hi b_hi c_hi ) ) ->
-              case length history `mod` 3 of
-                0 -> 𝕀 ( ℝ2 a_lo b_lo ) ( ℝ2 a_hi b_hi )
-                1 -> 𝕀 ( ℝ2 b_lo c_lo ) ( ℝ2 b_hi c_hi )
-                _ -> 𝕀 ( ℝ2 a_lo c_lo ) ( ℝ2 a_hi c_hi )
+            case cmpNat @n @d Proxy Proxy of
+              EQI -> id
+              LTI ->
+                -- If there are more equations (d) than variables (n),
+                -- pick a size n subset of the variables,
+                -- (go through all combinations cyclically).
+                let choices :: [ Vec n ( Fin d ) ]
+                    choices = choose @d @n
+                    choice :: Vec n ( Fin d )
+                    choice = choices !! ( length history `mod` length choices )
+                in \ u -> tabulate \ i -> index u ( choice ! i )
         }
 
     -- Did we reduce the box width by at least ε_eq