Split off root-isolation algorithms into Math.Root.Isolation

brush-strokes/brush-strokes.cabal

@@ -124,6 +124,7 @@ library
       , Math.Orientation
       , Math.Ring
       , Math.Roots
+      , Math.Root.Isolation
       , Debug.Utils
 
     other-modules:

brush-strokes/src/cusps/bench/Main.hs

@@ -40,6 +40,8 @@ import Numeric
   ( showFFloat )
 
 -- containers
+import qualified Data.IntMap.Strict as IntMap
+  ( fromList, toList )
 import Data.Sequence
   ( Seq )
 import qualified Data.Sequence as Seq
@@ -72,6 +74,7 @@ import Math.Linear
 import Math.Module
 import Math.Ring
   ( Transcendental )
+import Math.Root.Isolation
 
 --------------------------------------------------------------------------------
 
@@ -120,7 +123,12 @@ benchTestCase ( TestCase { testName, testBrushStroke, testCuspOptions, testStart
   let ( _, testStrokeFnI ) = brushStrokeFunctions testBrushStroke
       ( dunno, sols ) =
         foldMap ( \ ( i, ( _trees, ( mbCusps, defCusps ) ) ) -> ( map ( i , ) mbCusps, map ( i, ) defCusps ) ) $
-          computeCusps testCuspOptions testStrokeFnI testStartBoxes
+          IntMap.toList $
+            findCuspsIn testCuspOptions testStrokeFnI $
+              IntMap.fromList
+                [ ( i, [ box ] )
+                | ( i, box ) <- testStartBoxes
+                ]
   rnf dunno `seq` rnf sols `seq` return ()
   after <- getMonotonicTime
   let dt = after - before
@@ -175,9 +183,7 @@ benchCases :: [ TestCase ]
 benchCases =
   [ ellipseTestCase opts "full" ( 0, 1 ) pi $ defaultStartBoxes [ 0 .. 3 ] ]
   where
-    opts = defaultCuspFindingOptions
+    opts = defaultRootIsolationOptions
-
-
 
 --------------------------------------------------------------------------------
 
@@ -193,7 +199,7 @@ data TestCase =
   TestCase
   { testName        :: String
   , testBrushStroke :: BrushStroke
-  , testCuspOptions :: CuspFindingOptions
+  , testCuspOptions :: RootIsolationOptions
   , testStartBoxes  :: [ ( Int, Box ) ]
   }
 
@@ -250,7 +256,7 @@ take 10 $ Data.List.sortOn ( \ ( _, ℝ1 e, v) -> abs e + norm v ) [ let { v = m
 potentialCusp $ eval fI $ mkBox (0.5798, 0.5799) 3 (0.26798, 0.26799)
 > True
 
-let nbPotentialSols b = let ( _newtTrees, ( dunno, sols ) ) = findCuspsFrom NoPreconditioning 1e-7 fI b in length dunno + length sols
+let nbPotentialSols b = let ( _newtTrees, ( dunno, sols ) ) = isolateRootsIn NoPreconditioning 1e-7 fI b in length dunno + length sols
 
 nbPotentialSols $ mkBox (0.5798, 0.5799) 3 (0.26798, 0.26799)
 1
@@ -258,7 +264,7 @@ nbPotentialSols $ mkBox (0.5798, 0.5799) 3 (0.26798, 0.26799)
 nbPotentialSols $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799)
 0
 
-let showTrees b = map ( uncurry showIntervalNewtonTree ) $ fst $ findCuspsFrom NoPreconditioning 1e-7 fI b
+let showTrees b = map ( uncurry showIntervalNewtonTree ) $ fst $ isolateRootsIn NoPreconditioning 1e-7 fI b
 
 putStrLn $ unlines $ map Data.Tree.View.showTree $ showTrees $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799)
 
@@ -464,7 +470,7 @@ maximum [ _y $ sup $ unT $ _D12_v $ dEdsdcdt $ eval fI $ mkBox (t, t) 2 (s, s) |
 
 
 
-let showTrees b = map ( uncurry showIntervalNewtonTree ) $ fst $ findCuspsFrom NoPreconditioning 1e-7 fI b
+let showTrees b = map ( uncurry showIntervalNewtonTree ) $ fst $ isolateRootsIn NoPreconditioning 1e-7 fI b
 
 putStrLn $ unlines $ map Data.Tree.View.showTree $ showTrees $ mkBox (0.5486101933245248, 0.5486102071493622) 2 (0.548095036738487, 0.5480952)
 
@@ -482,7 +488,7 @@ showD float = showFFloat (Just 6) float ""
 
 --------------------------------------------------------------------------------
 
-ellipseTestCase :: CuspFindingOptions -> String -> ( Double, Double ) -> Double -> [ ( Int, Box ) ] -> TestCase
+ellipseTestCase :: RootIsolationOptions -> String -> ( Double, Double ) -> Double -> [ ( Int, Box ) ] -> TestCase
 ellipseTestCase opts str k0k1 rot startBoxes =
   TestCase
     { testName = "ellipse (" ++ str ++ ")"
@@ -512,7 +518,7 @@ trickyCusp2TestCase =
   TestCase
     { testName = "trickyCusp2"
     , testBrushStroke = trickyCusp2BrushStroke
-    , testCuspOptions = defaultCuspFindingOptions
+    , testCuspOptions = defaultRootIsolationOptions
     , testStartBoxes = defaultStartBoxes [ 0 .. 3 ]
     }
 
@@ -601,19 +607,6 @@ getStrokeFunctions brush sp0 crv =
         brush @3 @𝕀 proxy# singleton )
 {-# INLINEABLE getStrokeFunctions #-}
 
-computeCusps
-  :: CuspFindingOptions
-  -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-  -> [ ( Int, Box ) ]
-  -> [ ( Int, ( [ ( Box, RootFindingTree Box ) ], ( [ Box ], [ Box ] ) ) ) ]
-computeCusps cuspOpts eqs startBoxes =
-  map ( \ ( i, box ) -> ( i , ) $ findCuspsFrom cuspOpts ( mkEqn i ) box ) startBoxes
-    where
-      mkEqn 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 )
-        in ( eqs t `Seq.index` i ) s
-
 defaultStartBoxes :: [ Int ] -> [ ( Int, Box ) ]
 defaultStartBoxes is =
   [ ( i, mkBox (zero, one) (zero, one) ) | i <- is ]
@@ -623,8 +616,8 @@ getR1 (ℝ1 u) = u
 {-
 
 (f, fI) = brushStrokeFunctions $ ellipseBrushStroke (0,1) pi
-nbPotentialSols box = let ( _newtTrees, ( dunno, sols ) ) = findCuspsFrom NoPreconditioning 1e-7 fI box in length dunno + length sols
+nbPotentialSols box = let ( _newtTrees, ( dunno, sols ) ) = isolateRootsIn NoPreconditioning 1e-7 fI box in length dunno + length sols
-showTrees box = putStrLn $ unlines $ map Data.Tree.View.showTree $ map ( uncurry showIntervalNewtonTree ) $ fst $ findCuspsFrom NoPreconditioning 1e-7 fI box
+showTrees box = putStrLn $ unlines $ map Data.Tree.View.showTree $ map ( uncurry showIntervalNewtonTree ) $ fst $ isolateRootsIn NoPreconditioning 1e-7 fI box
 
 sol_t = 0.5486100729150693796677845183880669025324233347060776339185 :: Double
 sol_s = 0.5480950141859386853197594577293968665598143630958601978245 :: Double

brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs

@@ -317,9 +317,9 @@ evaluateCubic bez t =
   let inf_bez = fmap inf bez
       sup_bez = fmap sup bez
       mins = fmap (Cubic.bezier @( T Double ) inf_bez)
-           $ inf t :| ( sup t : filter ( inside t ) ( Cubic.extrema inf_bez ) )
+           $ inf t :| ( sup t : filter ( `inside` t ) ( Cubic.extrema inf_bez ) )
       maxs = fmap (Cubic.bezier @( T Double ) sup_bez)
-           $ inf t :| ( sup t : filter ( inside t ) ( Cubic.extrema sup_bez ) )
+           $ inf t :| ( sup t : filter ( `inside` t ) ( Cubic.extrema sup_bez ) )
   in 𝕀 ( minimum mins ) ( maximum maxs )
 
 -- | Evaluate a quadratic Bézier curve, when both the coefficients and the
@@ -330,9 +330,9 @@ evaluateQuadratic bez t =
   let inf_bez = fmap inf bez
       sup_bez = fmap sup bez
       mins = fmap (Quadratic.bezier @( T Double ) inf_bez)
-           $ inf t :| ( sup t : filter ( inside t ) ( Quadratic.extrema inf_bez ) )
+           $ inf t :| ( sup t : filter ( `inside` t ) ( Quadratic.extrema inf_bez ) )
       maxs = fmap (Quadratic.bezier @( T Double ) sup_bez)
-           $ inf t :| ( sup t : filter ( inside t ) ( Quadratic.extrema sup_bez ) )
+           $ inf t :| ( sup t : filter ( `inside` t ) ( Quadratic.extrema sup_bez ) )
   in 𝕀 ( minimum mins ) ( maximum maxs )
 
 {-

brush-strokes/src/lib/Math/Interval.hs

@@ -59,8 +59,8 @@ width ( 𝕀 lo hi ) = hi - lo
 nonDecreasing :: ( a -> b ) -> 𝕀 a -> 𝕀 b
 nonDecreasing f ( 𝕀 lo hi ) = 𝕀 ( f lo ) ( f hi )
 
-inside :: Ord a => 𝕀 a -> a -> Bool
+inside :: Ord a => a -> 𝕀 a -> Bool
-inside ( 𝕀 lo hi ) x = x >= lo && x <= hi
+inside x ( 𝕀 lo hi ) = x >= lo && x <= hi
 
 deriving via ViaAbelianGroup ( T ( 𝕀 Double ) )
   instance Semigroup ( T ( 𝕀 Double ) )

brush-strokes/src/lib/Math/Linear.hs

@@ -12,7 +12,7 @@ module Math.Linear
   , ℝ(..), T(.., V2, V3, V4)
   , Fin(..), MFin(..)
   , RepDim, RepresentableQ(..)
-  , Representable(..), injection, projection
+  , Representable(..), set, injection, projection
   , Vec(..), (!), find, zipIndices, unzip
 
   ) where

brush-strokes/src/lib/Math/Linear/Internal.hs

@@ -6,7 +6,9 @@ module Math.Linear.Internal
   , Fin(..), MFin(..)
   , RepDim
   , RepresentableQ(..)
-  , Representable(..), projection, injection
+  , Representable(..)
+  , set
+  , projection, injection
   )
   where
 
@@ -111,6 +113,12 @@ class 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
 
 projection :: ( Representable r u, Representable r v )
            => ( Fin ( RepDim v ) -> Fin ( RepDim u ) )
@@ -138,20 +146,24 @@ instance RepresentableQ Double ( ℝ 0 ) where
 
 instance RepresentableQ Double ( ℝ 1 ) where
   tabulateQ f = [|| ℝ1 $$( f ( Fin 1 ) ) ||]
-  indexQ p _ = [|| unℝ1 $$p ||]
+  indexQ p = \ case
+    Fin 1 -> [|| unℝ1 $$p ||]
+    Fin i -> error $ "invalid index for ℝ 1: " ++ show i
 
 instance RepresentableQ Double ( ℝ 2 ) where
   tabulateQ f = [|| ℝ2 $$( f ( Fin 1 ) ) $$( f ( Fin 2 ) ) ||]
   indexQ p = \ case
     Fin 1 -> [|| _ℝ2_x $$p ||]
-    _     -> [|| _ℝ2_y $$p ||]
+    Fin 2 -> [|| _ℝ2_y $$p ||]
+    Fin i -> error $ "invalid index for ℝ 2: " ++ show i
 
 instance RepresentableQ Double ( ℝ 3 ) where
   tabulateQ f = [|| ℝ3 $$( f ( Fin 1 ) ) $$( f ( Fin 2 ) ) $$( f ( Fin 3 ) ) ||]
   indexQ p = \ case
     Fin 1 -> [|| _ℝ3_x $$p ||]
     Fin 2 -> [|| _ℝ3_y $$p ||]
-    _     -> [|| _ℝ3_z $$p ||]
+    Fin 3 -> [|| _ℝ3_z $$p ||]
+    Fin i -> error $ "invalid index for ℝ 3: " ++ show i
 
 instance RepresentableQ Double ( ℝ 4 ) where
   tabulateQ f = [|| ℝ4 $$( f ( Fin 1 ) ) $$( f ( Fin 2 ) ) $$( f ( Fin 3 ) ) $$( f ( Fin 4 ) ) ||]
@@ -159,7 +171,8 @@ instance RepresentableQ Double ( ℝ 4 ) where
     Fin 1 -> [|| _ℝ4_x $$p ||]
     Fin 2 -> [|| _ℝ4_y $$p ||]
     Fin 3 -> [|| _ℝ4_z $$p ||]
-    _     -> [|| _ℝ4_w $$p ||]
+    Fin 4 -> [|| _ℝ4_w $$p ||]
+    Fin i -> error $ "invalid index for ℝ 4: " ++ show i
 
 instance Representable Double ( ℝ 0 ) where
   tabulate _ = ℝ0
@@ -170,7 +183,9 @@ instance Representable Double ( ℝ 0 ) where
 instance Representable Double ( ℝ 1 ) where
   tabulate f = ℝ1 ( f ( Fin 1 ) )
   {-# INLINE tabulate #-}
-  index p _ = unℝ1 p
+  index p = \ case
+    Fin 1 -> unℝ1 p
+    Fin i -> error $ "invalid index for ℝ 1: " ++ show i
   {-# INLINE index #-}
 
 instance Representable Double ( ℝ 2 ) where
@@ -178,7 +193,8 @@ instance Representable Double ( ℝ 2 ) where
   {-# INLINE tabulate #-}
   index p = \ case
     Fin 1 -> _ℝ2_x p
-    _     -> _ℝ2_y p
+    Fin 2 -> _ℝ2_y p
+    Fin i -> error $ "invalid index for ℝ 2: " ++ show i
   {-# INLINE index #-}
 
 instance Representable Double ( ℝ 3 ) where
@@ -187,7 +203,8 @@ instance Representable Double ( ℝ 3 ) where
   index p = \ case
     Fin 1 -> _ℝ3_x p
     Fin 2 -> _ℝ3_y p
-    _     -> _ℝ3_z p
+    Fin 3 -> _ℝ3_z p
+    Fin i -> error $ "invalid index for ℝ 3: " ++ show i
   {-# INLINE index #-}
 
 instance Representable Double ( ℝ 4 ) where
@@ -197,5 +214,6 @@ instance Representable Double ( ℝ 4 ) where
     Fin 1 -> _ℝ4_x p
     Fin 2 -> _ℝ4_y p
     Fin 3 -> _ℝ4_z p
-    _     -> _ℝ4_w p
+    Fin 4 -> _ℝ4_w p
+    Fin i -> error $ "invalid index for ℝ 4: " ++ show i
   {-# INLINE index #-}

brush-strokes/src/lib/Math/Root/Isolation.hs

@@ -0,0 +1,746 @@
+module Math.Root.Isolation
+  ( -- * Root isolation using interval arithmetic
+    Box, BoxHistory
+  , isolateRootsIn
+
+    -- ** Configuration of the root isolation methods
+  , RootIsolationOptions(..), defaultRootIsolationOptions
+  , RootIsolationAlgorithm(..), defaultRootIsolationAlgorithms
+
+    -- *** Options for the bisection method
+  , BisectionOptions(..)
+
+    -- *** Options for the interval Newton method with Gauss–Seidel step
+  , GaussSeidelOptions(..), Preconditioner(..)
+
+    -- *** Options for the @box(1)@-consistency algorithm
+  , Box1Options(..)
+
+    -- *** Options for the @box(2)@-consistency algorithm
+  , Box2Options(..)
+
+    -- ** Trees recording search space of root isolation algorithms
+  , RootIsolationTree(..), showRootIsolationTree
+  , RootIsolationStep(..)
+  , RootIsolationLeaf(..)
+  )
+  where
+
+-- base
+import Control.Arrow
+  ( first )
+import Control.Monad
+  ( when )
+import Data.Bifunctor
+  ( Bifunctor(bimap) )
+import Data.List
+  ( nub, partition)
+import Data.List.NonEmpty
+  ( NonEmpty )
+import qualified Data.List.NonEmpty as NE
+  ( NonEmpty(..), last )
+import Numeric
+  ( showFFloat )
+
+-- containers
+import Data.Tree
+  ( Tree(..) )
+
+-- transformers
+import Control.Monad.Trans.State.Strict as State
+  ( State, evalState, get, put )
+import Control.Monad.Trans.Writer.CPS
+  ( Writer, runWriter, tell )
+
+-- MetaBrush
+import Math.Algebra.Dual
+  ( D, D1𝔸2(..) )
+import Math.Epsilon
+  ( nearZero )
+import Math.Float.Utils
+  ( succFP, prevFP )
+import Math.Interval
+import Math.Linear
+import Math.Module
+  ( Module(..) )
+import qualified Math.Ring as Ring
+
+--import Debug.Utils
+
+--------------------------------------------------------------------------------
+
+-- | A tree recording the steps taken when doing cusp finding.
+data RootIsolationTree d
+  = RootIsolationLeaf (RootIsolationLeaf d)
+  | RootIsolationStep RootIsolationStep [(d, RootIsolationTree d)]
+
+-- | A description of a step taken in cusp-finding.
+data RootIsolationStep
+  = GaussSeidelStep
+  | BisectionStep (String, Double)
+  | Box1Step
+  | Box2Step
+
+instance Show RootIsolationStep where
+  showsPrec _ GaussSeidelStep = showString "GS"
+  showsPrec _ ( BisectionStep (coord, w) )
+    = showString "bis "
+    . showParen True
+      ( showString coord
+      . showString " = "
+      . showsPrec 0 w
+      )
+  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 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 )
+
+showArea :: Double -> String
+showArea area = "(area " ++ showFFloat (Just 6) area "" ++ ")"
+
+type Box = 𝕀ℝ 2
+type BoxHistory = [ NE.NonEmpty ( RootIsolationStep, Box ) ]
+
+-- | Options for the root isolation methods in 'findCusps' and 'isolateRootsIn'.
+data RootIsolationOptions
+  = 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 )
+  }
+
+data RootIsolationAlgorithm
+   -- | Bisection step.
+  = Bisection !BisectionOptions
+  -- | Gauss–Seidel step with the given preconditioning method.
+  | GaussSeidel !GaussSeidelOptions
+  -- | @box(1)@-consistency.
+  | Box1 !Box1Options
+  -- | @box(2)@-consistency.
+  | Box2 !Box2Options
+
+-- | Options for the bisection method.
+data BisectionOptions =
+  BisectionOptions
+  { canHaveSols :: !( ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) ) -> Box -> Bool )
+  , fallbackBisectionDim :: !( [ ( RootIsolationStep, Box ) ] -> BoxHistory -> ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) ) -> ( Int, String ) )
+  }
+
+-- | Options for the interval Gauss–Seidel method.
+data GaussSeidelOptions =
+  GaussSeidelOptions
+    { gsPreconditioner :: !Preconditioner }
+
+-- | Options for the @box(1)@-consistency method.
+data Box1Options =
+  Box1Options
+    { box1EpsEq :: !Double }
+
+-- | Options for the @box(2)@-consistency method.
+data Box2Options =
+  Box2Options
+    { box2EpsEq :: !Double
+    , box2LambdaMin :: !Double
+    }
+
+defaultRootIsolationOptions :: RootIsolationOptions
+defaultRootIsolationOptions =
+  RootIsolationOptions
+    { minWidth
+    , cuspFindingAlgorithms = defaultRootIsolationAlgorithms minWidth narrowAbs
+    }
+  where
+    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 } ) ]
+  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 )
+
+defaultGaussSeidelOptions :: GaussSeidelOptions
+defaultGaussSeidelOptions = GaussSeidelOptions { gsPreconditioner = InverseMidJacobian }
+
+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
+
+            -- 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, "")
+      }
+
+-- | Use the following algorithms using interval arithmetic in order
+-- to isolate roots:
+--
+--  - interval Newton method with Gauss–Seidel step for inversion
+--    of the interval Jacobian,
+--  - coordinate-wise Newton method (@box(1)@-consistency algorithm),
+--  - @box(2)@-consistency algorithm,
+--  - bisection.
+--
+-- Returns @(tree, (dunno, sols))@ where @tree@ is the full tree (useful for debugging),
+-- @sols@ are boxes that contain a unique solution (and to which Newton's method
+-- 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 ] ) )
+isolateRootsIn
+  ( RootIsolationOptions
+    { minWidth
+    , cuspFindingAlgorithms
+    }
+  )
+  eqs initBox = runWriter $ map ( initBox , ) <$> go [] initBox
+
+  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.
+      = 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
+      where
+        D12 iRange _ _ = eqs $ 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi )
+
+    -- Run a round of cusp finding strategies, then recur.
+    doStrategies
+      :: BoxHistory
+      -> NonEmpty RootIsolationAlgorithm
+      -> Box
+      -> Writer ( [ Box ], [ Box ] )
+                [ RootIsolationTree Box ]
+    doStrategies prevRoundsHist = do_strats []
+      where
+        do_strats :: [ ( RootIsolationStep, Box ) ]
+                  -> NE.NonEmpty RootIsolationAlgorithm
+                  -> Box
+                  -> Writer ( [ Box ], [ Box ] )
+                            [ RootIsolationTree Box ]
+        do_strats thisRoundHist ( algo NE.:| algos ) box = do
+          -- Run one strategy in the round.
+          ( step, boxes ) <- doStrategy thisRoundHist prevRoundsHist eqs minWidth algo box
+          case algos of
+            -- If there are other algorithms to run in this round, run them next.
+            nextAlgo : otherAlgos ->
+              let thisRoundHist' = ( step, box ) : thisRoundHist
+              in recur step ( do_strats thisRoundHist' ( nextAlgo NE.:| otherAlgos ) ) boxes
+            -- Otherwise, we have done one full round of strategies.
+            -- Recur back to the top (calling 'go').
+            [] ->
+              let thisRoundHist' = ( step, box ) NE.:| thisRoundHist
+                  history = thisRoundHist' : prevRoundsHist
+              in recur step ( go history ) boxes
+
+    recur :: RootIsolationStep
+          -> ( Box -> Writer ( [ Box ], [ Box ] ) [ RootIsolationTree Box ] )
+          -> [ Box ]
+          -> Writer ( [Box], [Box] ) [ RootIsolationTree Box ]
+    recur step r cands = do
+      rest <- traverse ( \ c -> do { trees <- r c; return [ (c, tree) | tree <- trees ] } ) cands
+      return [ RootIsolationStep step $ concat rest ]
+
+
+-- | 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 roundHistory previousRoundsHistory eqs minWidth algo box =
+  case algo of
+    GaussSeidel ( GaussSeidelOptions { gsPreconditioner } ) -> do
+      boxes <- intervalGaussSeidel eqs gsPreconditioner 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 )
+
+-- | 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 ) )
+  -- 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
+
+        ( 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 )
+
+-- | 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 ]
+intervalGaussSeidel
+  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
+
+  = 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 )
+
+        -- 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 )
+    in
+         -- If the Gauss–Seidel step was a contraction, then the box
+         -- contains a unique solution (by the Banach fixed point theorem).
+         --
+         -- These boxes can thus be directly added to the solution set:
+         -- 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)
+               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.
+
+-- 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 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 )
+-- 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"
+
+-- | Compute the intersection of two intervals.
+--
+-- Returns whether the first interval is a strict subset of the second interval,
+-- or the intersection is a single point.
+intersect :: 𝕀 Double -> 𝕀 Double -> [ ( 𝕀 Double, Bool ) ]
+intersect ( 𝕀 lo1 hi1 ) ( 𝕀 lo2 hi2 )
+  | lo > hi
+  = [ ]
+  | otherwise
+  = [ ( 𝕀 lo hi, ( lo1 > lo2 && hi1 < hi2 ) || ( lo == hi ) ) ]
+  where
+    lo = max lo1 lo2
+    hi = min hi1 hi2
+
+-- | Preconditioner to use with the interval Gauss–Seidel method.
+data Preconditioner
+  = NoPreconditioning
+  | InverseMidJacobian
+  deriving stock Show
+
+-- | An implementation of "bc_enforce" from the paper
+-- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
+--
+-- See also
+-- "Presentation of a highly tuned multithreaded interval solver for underdetermined and well-determined nonlinear systems"
+makeBox1Consistent
+  :: ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) )
+  -> Double -> Double
+  -> Box -> [ Box ]
+makeBox1Consistent eqs minWidth epsEq x =
+  ( `State.evalState` False ) $
+    pipeFunctionsWhileTrue ( allNarrowingOperators epsEq minWidth eqs ) 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 ) )
+  -> Double -> Double -> Double
+  -> Box -> Box
+makeBox2Consistent eqs minWidth epsEq lambdaMin x0 = ( `State.evalState` False ) $ doLoop 0.25 x0
+  where
+    doLoop :: Double -> Box -> State Bool Box
+    doLoop lambda x = do
+      x'  <- boundConsistency get_t set_t lambda x
+      x'' <- boundConsistency get_s set_s lambda x'
+      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 getter setter lambda box = do
+      let x@( 𝕀 x_inf x_sup ) = getter box
+          c1 = ( 1 - lambda ) * x_inf + lambda * x_sup
+          c2 = lambda * x_inf + ( 1 - lambda ) * x_sup
+          x'_inf =
+            case makeBox1Consistent eqs minWidth epsEq ( setter ( 𝕀 x_inf c1 ) box ) of
+              []  -> c1
+              x's -> minimum $ map ( inf . getter ) x's
+          x'_sup =
+            case makeBox1Consistent eqs minWidth epsEq ( setter ( 𝕀 c2 x_sup ) box ) of
+              []  -> c2
+              x's -> maximum $ map ( sup . getter ) x's
+          x' = 𝕀 x'_inf x'_sup
+      when ( width x - width x' >= epsEq ) $
+        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 =
+  case meth of
+    NoPreconditioning
+      -> ( a, 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 )
+      | 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 )
+  where
+    𝕀 b1_lo b1_hi = scaleInterval s ( 𝕀 a1_lo a1_hi )
+    𝕀 b2_lo b2_hi = scaleInterval s ( 𝕀 a2_lo a2_hi )
+
+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
+
+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
+
+-- Use the univariate interval Newton method to narrow from the left
+-- a candidate interval.
+--
+-- See Algorithm 5 (Procedure left_narrow) in
+-- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
+-- by Bartłomiej Jacek Kubica, 2017
+leftNarrow :: Double
+           -> Double
+           -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) )
+           -> 𝕀 Double
+           -> [ 𝕀 Double ]
+leftNarrow eps_equal eps_bisection ff' = left_narrow
+  where
+    left_narrow ( 𝕀 x_inf x_sup ) =
+      go x_sup ( 𝕀 x_inf ( if x_inf == x_sup then x_inf else succFP x_inf ) )
+    go x_sup x_left =
+      let ( f_x_left, _f'_x_left ) = ff' x_left
+      in
+        if inf f_x_left <= 0 && sup f_x_left >= 0
+        then [ 𝕀 ( inf x_left ) x_sup ]
+        else
+          let x = 𝕀 ( sup x_left ) x_sup
+              ( _f_x, f'_x ) = ff' x
+              x's = do delta <- f_x_left `extendedDivide` f'_x --f'_x_left
+                       let x_new = x_left - delta
+                       map fst $ x_new `intersect` x
+          in
+            if | null x's
+               -> []
+               | ( width x - sum ( map width x's ) ) < eps_equal
+               -- TODO: do a check with the relative reduction in size?
+               -> x's
+               | otherwise
+               -> do
+                x' <- x's
+                if sup x' - inf x' < eps_bisection
+                then return x'
+                else left_narrow =<< bisectI x'
+
+-- TODO: de-duplicate with 'leftNarrow'?
+rightNarrow :: Double
+            -> Double
+            -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) )
+            -> 𝕀 Double
+            -> [ 𝕀 Double ]
+rightNarrow eps_equal eps_bisection ff' = right_narrow
+  where
+    right_narrow ( 𝕀 x_inf x_sup ) =
+      go x_inf ( 𝕀 ( if x_inf == x_sup then x_sup else prevFP x_sup ) x_sup )
+    go x_inf x_right =
+      let ( f_x_right, _f'_x_right ) = ff' x_right
+      in
+        if inf f_x_right <= 0 && sup f_x_right >= 0
+        then [ 𝕀 x_inf ( sup x_right ) ]
+        else
+          let x = 𝕀 x_inf ( inf x_right )
+              ( _f_x, f'_x ) = ff' x
+              x's = do delta <- f_x_right `extendedDivide` f'_x --f'_x_right
+                       let x_new = x_right - delta
+                       map fst $ x_new `intersect` x
+          in
+            if | null x's
+               -> []
+               | ( width x - sum ( map width x's ) ) < eps_equal
+               -- TODO: do a check with the relative reduction in size?
+               -> x's
+               | otherwise
+               -> do
+                x' <- x's
+                if sup x' - inf x' < eps_bisection
+                then return x'
+                else right_narrow =<< bisectI x'
+
+bisectI :: 𝕀 Double -> [ 𝕀 Double ]
+bisectI x@( 𝕀 x_inf x_sup )
+  | x_inf == x_sup
+  = [ x ]
+  | otherwise
+  = [ 𝕀 x_inf x_mid, 𝕀 x_mid x_sup ]
+  where x_mid = 0.5 * ( x_inf + x_sup )
+
+-- | Apply each function in turn, piping the output of one function into
+-- the next.
+--
+-- Once all the functions have been applied, check whether the Bool is True.
+-- If it is, go around again with all the functions; otherwise, stop.
+pipeFunctionsWhileTrue :: [ a -> State Bool [ a ] ] -> a -> State Bool [ a ]
+pipeFunctionsWhileTrue fns = go fns
+  where
+    go [] x = do
+      doAnotherRound <- State.get
+      if doAnotherRound
+      then do { State.put False ; go fns x }
+      else return [ x ]
+    go ( f : fs ) x = do
+      xs <- f x
+      concat <$> traverse ( go fs ) xs
+
+allNarrowingOperators
+  :: Double
+  -> Double
+  -> ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) )
+  -> [ Box -> State Bool [ Box ] ]
+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 )
+    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 ] ]
+      ]
+  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 )