diff --git a/MetaBrush.cabal b/MetaBrush.cabal
index 3f9edf9..641b871 100644
--- a/MetaBrush.cabal
+++ b/MetaBrush.cabal
@@ -187,49 +187,6 @@ library metabrushes
       , bytestring
            >= 0.10.10.0  && < 0.12
 
-executable cusps
-
-    import:
-        common
-
-    hs-source-dirs:
-        src/cusps
-
-    default-language:
-        Haskell2010
-
-    main-is:
-        Main.hs
-
-    other-modules:
-        Math.Interval.Abstract
-
-
-executable convert-metafont
-
-    import:
-        common, extras
-
-    hs-source-dirs:
-        src/convert
-
-    default-language:
-        Haskell2010
-
-    main-is:
-        Main.hs
-
-    other-modules:
-        MetaBrush.MetaFont.Convert
-
-    build-depends:
-        metabrushes,
-        diagrams-contrib,
-        diagrams-lib,
-        linear,
-        parsec
-
-
 executable MetaBrush
 
     import:
diff --git a/brush-strokes/bench/Main.hs b/brush-strokes/bench/Main.hs
deleted file mode 100644
index 929fb61..0000000
--- a/brush-strokes/bench/Main.hs
+++ /dev/null
@@ -1,370 +0,0 @@
-{-# LANGUAGE PolyKinds #-}
-{-# LANGUAGE ScopedTypeVariables #-}
-{-# LANGUAGE UndecidableInstances #-}
-
-module Main
-  ( main
-
-    -- Testing
-  , TestCase(..)
-  , testCases
-  , testCaseStrokeFunctions
-  , eval
-  , mkVal, mkBox
-  , potentialCusp
-  , dEdsdcdt
-  )
-  where
-
--- base
-import Control.Concurrent.MVar
-  ( newMVar )
-import Data.Coerce
-  ( coerce )
-import Data.Foldable
-  ( for_ )
-import Data.List
-  ( intercalate )
-import GHC.Exts
-  ( Proxy#, proxy# )
-import GHC.Generics
-  ( Generic )
-import GHC.TypeNats
-  ( type (-) )
-import Numeric
-  ( showFFloat )
-
--- containers
-import Data.Sequence
-  ( Seq )
-import qualified Data.Sequence as Seq
-  ( index )
-import Data.Tree
-  ( foldTree )
-
--- brush-strokes
-import Calligraphy.Brushes
-import Debug.Utils
-  ( logToFile )
-import Math.Algebra.Dual
-import Math.Bezier.Spline
-import Math.Bezier.Stroke
-import Math.Bezier.Stroke.EnvelopeEquation
-import Math.Differentiable
-import Math.Interval
-import Math.Linear
-import Math.Module
-import Math.Ring
-  ( Transcendental )
-
---------------------------------------------------------------------------------
-
-main :: IO ()
-main = for_ testCases $ \ testCase@( TestCase { testName, testAlgorithmParams } ) -> do
-  let ( _, testStrokeFnI ) = testCaseStrokeFunctions testCase
-      ( newtTrees, ( dunno, sols ) ) = computeCusps testAlgorithmParams testStrokeFnI
-      showedTrees = map ( uncurry showIntervalNewtonTree ) newtTrees
-  putStrLn $ unlines $
-    [ ""
-    , "Test case '" ++ testName ++ "':" ] ++
-    map ( "    " ++ )
-      [ " #sols: " ++ show (length sols)
-      , "#dunno: " ++ show (length dunno)
-      , "#trees: " ++ show @Int (sum @_ @Int $ map (foldTree ( \ _ bs -> 1 + sum bs )) showedTrees)
-      , " dunno: " ++ show dunno
-      , "  sols: " ++ show sols
-      ]
-  --logFileMVar <- newMVar "logs/trickyCusp.log"
-  --logToFile logFileMVar (unlines logLines)
-  --  `seq` return ()
-
-testCases :: [ TestCase ]
-testCases = [ ellipse , trickyCusp2 ]
-
---------------------------------------------------------------------------------
-
-data TestCase =
-  forall nbParams. ParamsCt nbParams =>
-  TestCase
-  { testName   :: !String
-  , testBrush  :: !( Brush nbParams )
-  , testStroke :: !( Point nbParams, Curve Open () ( Point nbParams ))
-  , testAlgorithmParams :: !CuspAlgorithmParams
-  }
-
-testCaseStrokeFunctions
-  :: TestCase
-  -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () )
-     , 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-testCaseStrokeFunctions ( TestCase { testStroke = ( sp0, crv ), testBrush } ) =
-  getStrokeFunctions testBrush sp0 crv
-
--- Utilities to use in GHCi to help debugging.
-
-eval
-  :: ( I i ( ℝ 1 ) -> Seq ( I i ( ℝ 1 ) -> StrokeDatum k i ) )
-  -> ( I i ( ℝ 1 ), Int, I i ( ℝ 1 ) )
-  -> StrokeDatum k i
-eval f ( t, i, s ) = ( f t `Seq.index` i ) s
-
-mkVal :: Double -> Int -> Double -> ( ℝ 1, Int, ℝ 1 )
-mkVal t i s = ( ℝ1 t, i, ℝ1 s )
-
-mkBox :: ( Double, Double ) -> Int -> ( Double, Double ) -> Box
-mkBox ( t_min, t_max ) i ( s_min, s_max ) =
-  ( 𝕀 ( ℝ1 t_min ) ( ℝ1 t_max ) , i, 𝕀 ( ℝ1 s_min ) ( ℝ1 s_max ) )
-
-potentialCusp :: StrokeDatum 3 𝕀 -> Bool
-potentialCusp
-  ( StrokeDatum
-    { ee = D22 { _D22_v = 𝕀 ( ℝ1 ee_min ) ( ℝ1 ee_max ) }
-    , 𝛿E𝛿sdcdt = D12 { _D12_v = T ( 𝕀 ( ℝ2 vx_min vy_min ) ( ℝ2 vx_max vy_max ) )}
-    }
-  ) = ee_min <= 0 && ee_max >= 0
-    && vx_min <= 0 && vx_max >= 0
-    && vy_min <= 0 && vy_max >= 0
-
-dEdsdcdt :: StrokeDatum k i -> D ( k - 2 ) ( I i ( ℝ 2 ) ) ( T ( I i ( ℝ 2 ) ) )
-dEdsdcdt ( StrokeDatum { 𝛿E𝛿sdcdt = v } ) = v
-
-{-
-let (f, fI) = testCaseStrokeFunctions trickyCusp2
-
-take 10 $ Data.List.sortOn ( \ ( _, ℝ1 e, v) -> abs e + norm v ) [ let { v = mkVal x 3 y; d = eval f v } in ( v, _D12_v $ ee d, _D0_v $ dEdsdcdt d ) | x <- [0.57,0.5701 .. 0.58], y <- [0.29,0.291..0.3] ]
-> ((ℝ1 0.5798800000000057,3,ℝ1 0.267980000000008),ℝ1 -2.8596965543670194e-4,V2 7.79559474412963e-2 2.0389671921293484e-2)
-
-potentialCusp $ eval fI $ mkBox (0.5798, 0.5799) 3 (0.26798, 0.26799)
-> True
-
-let nbPotentialSols b = let ( _newtTrees, ( dunno, sols ) ) = intervalNewtonGSFrom NoPreconditioning 1e-7 fI b in length dunno + length sols
-
-nbPotentialSols $ mkBox (0.5798, 0.5799) 3 (0.26798, 0.26799)
-1
-
-nbPotentialSols $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799)
-0
-
-let showTrees b = map ( uncurry showIntervalNewtonTree ) $ fst $ intervalNewtonGSFrom NoPreconditioning 1e-7 fI b
-
-putStrLn $ unlines $ map Data.Tree.View.showTree $ showTrees $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799)
-
-([ℝ1 0.5798, ℝ1 0.675],3,[ℝ1 0.26798, ℝ1 0.26799]) (area 0.000001) N []
- └─ ([ℝ1 0.5973000285624527, ℝ1 0.6750000000000002],3,[ℝ1 0.26798, ℝ1 0.26799000000000006]) (area 0.000001) NoSolution "ee" ([ℝ1 0.5973000285624527, ℝ1 0.6750000000000002],3,[ℝ1 0.26798, ℝ1 0.26799000000000006])
-
-eval fI $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799)
-> D12 {_D12_v = T[ℝ2 -10088.6674944889 -3281.3820867312834, ℝ2 4124.668381545453 4524.807156085763], _D12_dx = TT[ℝ2 -173746.97965005718 -33281.18494907289, ℝ2 298.2609121556852 23639.772884799597], _D12_dy = TT[ℝ2 -18454.27716258352 -28337.509817580823, ℝ2 1163.6949532017436 -13936.383137525536]}}
-i.e.
-> f = [ℝ2 -10088.6674944889 -3281.3820867312834, ℝ2 4124.668381545453 4524.807156085763]
-> f_t = [ℝ2 -173746.97965005718 -33281.18494907289, ℝ2 298.2609121556852 23639.772884799597]
-> f_s = [ℝ2 -18454.27716258352 -28337.509817580823, ℝ2 1163.6949532017436 -13936.383137525536]
-
-(f, fI) = testCaseStrokeFunctions trickyCusp2
-t = 𝕀 (ℝ1 0.5798) (ℝ1 0.675)
-s = 𝕀 (ℝ1 0.26798) (ℝ1 0.26799)
-t_mid = 0.5 * ( 0.5798 + 0.675 )
-s_mid = 0.5 * ( 0.26798 + 0.26799 )
-D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) = dEdsdcdt $ eval fI (t, 3, s)
-t' = coerce ( (-) @( 𝕀 Double ) ) t ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1
-s' = coerce ( (-) @( 𝕀 Double ) ) s ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1
-a = ( f_t, f_s )
-b = negV2 $ singleton $ midV2 f
-[((t2', s2'), isContr)] = gaussSeidel a b (t', s')
-t2 = coerce ( (+) @( 𝕀 Double ) ) t2' ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1
-s2 = coerce ( (+) @( 𝕀 Double ) ) s2' ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1
-
-t2
-> [ℝ1 0.6102365832093095, ℝ1 0.6750000000000002]
-s2
-> [ℝ1 0.26798, ℝ1 0.26799000000000006]
-
-
-let ( 𝕀 ( ℝ2 a11_lo a21_lo ) ( ℝ2 a11_hi a21_hi ), 𝕀 ( ℝ2 a12_lo a22_lo ) ( ℝ2 a12_hi a22_hi ) ) = a
-let ( 𝕀 ( ℝ2 b1_lo b2_lo ) ( ℝ2 b1_hi b2_hi ) ) = b
-let ( 𝕀 ( ℝ1 x1_lo ) ( ℝ1 x1_hi ), 𝕀 ( ℝ1 x2_lo ) ( ℝ1 x2_hi ) ) = ( t', s' )
-
-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
-
-( b1 - a12 * x2 )
-> [2981.90728508591, 2982.0918278575364]
-
-extendedRecip a11
-
--}
-
-negV2 :: 𝕀ℝ 2 -> 𝕀ℝ 2
-negV2 ( 𝕀 ( ℝ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 )
-
-midV2 :: 𝕀ℝ 2 -> ℝ 2
-midV2 ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) =
-  ℝ2 ( 0.5 * ( x_lo + x_hi ) ) ( 0.5 * ( y_lo + y_hi ) )
-
-logLines :: [ String ]
-logLines =
-  [ "f = dE/ds * dc/dt: f, df/dt, df/ds"
-  , "{" ++
-    (intercalate ","
-     [ "{" ++ showD (midPoint t) ++ "," ++ showD (midPoint s) ++ ",{" ++ intercalate "," vals ++ "}}"
-     | t <- map ( around 0.5798  ) [-0.05, -0.049.. 0.05]
-     , let i = 3
-     , s <- map ( around 0.26798 ) [-0.05, -0.049.. 0.05]
-     , let StrokeDatum
-             { 𝛿E𝛿sdcdt = D12 (T f) (T (T f_t)) (T (T f_s))
-             } = (curvesI t `Seq.index` i) s
-           ℝ2 vx vy = midPoint2 f
-           ℝ2 vx_t vy_t = midPoint2 f_t
-           ℝ2 vx_s vy_s = midPoint2 f_s
-           vals = [ "{" ++ showD vx ++ "," ++ showD vy ++ "}"
-                  , "{" ++ showD vx_t ++ "," ++ showD vy_t ++ "}"
-                  , "{" ++ showD vx_s ++ "," ++ showD vy_s ++ "}"
-                  ]
-     ]
-    ) ++ "}"
-  ]
-  where
-    around :: Double -> Double -> 𝕀ℝ 1
-    around z0 z = 𝕀 ( ℝ1 ( z + z0 - 1e-6 ) ) ( ℝ1 ( z + z0 + 1e-6 ) )
-    ( _, curvesI ) = testCaseStrokeFunctions trickyCusp2
-    midPoint (𝕀 (ℝ1 lo) (ℝ1 hi)) = 0.5 * ( lo + hi )
-    midPoint2 (𝕀 (ℝ2 lo_x lo_y) (ℝ2 hi_x hi_y))
-      = ℝ2 ( 0.5 * ( lo_x + hi_x ) ) ( 0.5 * ( lo_y + hi_y ) )
-
-showD :: Double -> String
-showD float = showFFloat (Just 6) float ""
-
---------------------------------------------------------------------------------
-
-ellipse :: TestCase
-ellipse =
-  TestCase
-    { testName = "ellipse"
-    , testBrush = ellipseBrush
-    , testStroke = ( p0, LineTo ( NextPoint p1 ) () )
-    , testAlgorithmParams =
-        CuspAlgorithmParams
-         { preconditioning = NoPreconditioning
-         , maxWidth = 1e-7
-         }
-    }
-  where
-    mkPt x y w h phi =
-      Point
-        { pointCoords = ℝ2 x y
-        , pointParams = Params $ ℝ3 w h phi
-        }
-    p0 = mkPt 0 0 10 25 0
-    p1 = mkPt 100 0 15 40 pi
-
-trickyCusp2 :: TestCase
-trickyCusp2 =
-  TestCase
-    { testName = "trickyCusp2"
-    , testBrush = circleBrush
-    , testStroke = ( p0, Bezier3To p1 p2 ( NextPoint p3 ) () )
-    , testAlgorithmParams =
-        CuspAlgorithmParams
-         { preconditioning = NoPreconditioning
-         , maxWidth = 1e-7
-         }
-    }
-  where
-    mkPt x y =
-      Point
-        { pointCoords = ℝ2 x y
-        , pointParams = Params $ ℝ1 5.0
-        }
-    p0 = mkPt 5e+1 -5e+1
-    p1 = mkPt -7.72994362904069e+1 -3.124468786098509e+1
-    p2 = mkPt -5.1505430313958364e+1 -3.9826386521527986e+1
-    p3 = mkPt -5e+1 -5e+1
-
---------------------------------------------------------------------------------
-
-type ParamsCt nbParams
-  = ( Show ( ℝ nbParams )
-    , HasChainRule Double 2 ( ℝ nbParams )
-    , HasChainRule ( 𝕀 Double ) 3 ( 𝕀 ( ℝ nbParams ) )
-    , Applicative ( D 2 ( ℝ nbParams ) )
-    , Applicative ( D 3 ( ℝ nbParams ) )
-    , Traversable ( D 2 ( ℝ nbParams ) )
-    , Traversable ( D 3 ( ℝ nbParams ) )
-    , Representable Double ( ℝ nbParams )
-    , Module Double ( T ( ℝ nbParams ) )
-    , Module ( 𝕀 Double ) ( T ( 𝕀 ( ℝ nbParams ) ) )
-    , Module ( D 2 ( ℝ nbParams ) Double ) ( D 2 ( ℝ nbParams ) ( ℝ 2 ) )
-    , Module ( D 3 ( ℝ nbParams ) ( 𝕀 Double ) ) ( D 3 ( ℝ nbParams ) ( 𝕀 ( ℝ 2 ) ) )
-    , Transcendental ( D 2 ( ℝ nbParams ) Double )
-    , Transcendental ( D 3 ( ℝ nbParams ) ( 𝕀 Double ) )
-    )
-
-newtype Params nbParams = Params { getParams :: ( ℝ nbParams ) }
-deriving newtype instance Show ( ℝ nbParams ) => Show ( Params nbParams )
-
-data Point nbParams =
-  Point
-    { pointCoords :: !( ℝ 2 )
-    , pointParams :: !( Params nbParams ) }
-  deriving stock Generic
-deriving stock instance Show ( ℝ nbParams ) => Show ( Point nbParams )
-
-data CuspAlgorithmParams =
-  CuspAlgorithmParams
-    { preconditioning :: !Preconditioner
-    , maxWidth :: !Double
-    }
-  deriving stock Show
-
-type Brush nbParams
-  = forall {t} k (i :: t)
-  .  DiffInterp k i ( ℝ nbParams )
-  => Proxy# i
-  -> ( forall a. a -> I i a )
-  -> C k ( I i ( ℝ nbParams ) )
-         ( Spline Closed () ( I i ( ℝ 2 ) ) )
-
-getStrokeFunctions
-  :: forall nbParams
-  .  ParamsCt nbParams
-  => Brush nbParams
-      -- ^ brush shape
-  -> Point nbParams
-      -- ^ start point
-  -> Curve Open () ( Point nbParams )
-      -- ^ curve points
-  -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () )
-     , 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-getStrokeFunctions brush sp0 crv =
-  let
-    usedParams :: C 2 ( ℝ 1 ) ( ℝ nbParams )
-    path :: C 2 ( ℝ 1 ) ( ℝ 2 )
-    ( path, usedParams ) =
-      pathAndUsedParams @2 @() coerce id ( getParams . pointParams )
-        sp0 crv
-    usedParamsI :: C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ nbParams )
-    pathI :: C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 )
-    ( pathI, usedParamsI ) =
-      pathAndUsedParams @3 @𝕀 coerce singleton ( getParams . pointParams )
-        sp0 crv
-  in ( brushStrokeData @2 @( ℝ nbParams ) coerce coerce
-        path usedParams $
-        brush @2 @() proxy# id
-     , brushStrokeData @3 @( ℝ nbParams ) coerce coerce
-        pathI usedParamsI $
-        brush @3 @𝕀 proxy# singleton )
-{-# INLINEABLE getStrokeFunctions #-}
-
-computeCusps
-  :: CuspAlgorithmParams
-  -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-  -> ( [ ( Box, IntervalNewtonTree Box ) ], ( [ Box ], [ Box ] ) )
-computeCusps params =
-  intervalNewtonGS ( preconditioning params ) ( maxWidth params )
diff --git a/brush-strokes/brush-strokes.cabal b/brush-strokes/brush-strokes.cabal
index e825794..9ff3b55 100644
--- a/brush-strokes/brush-strokes.cabal
+++ b/brush-strokes/brush-strokes.cabal
@@ -83,13 +83,23 @@ common common
           base
             >= 4.19
 
+common extra
+
+    build-depends:
+        acts
+          ^>= 0.3.1.0
+      , generic-lens
+           >= 2.2     && < 2.3
+      , groups
+          ^>= 0.5.3
+
 library
 
     import:
-        common
+        common, extra
 
     hs-source-dirs:
-        src
+        src/lib
 
     default-language:
         Haskell2010
@@ -129,9 +139,6 @@ library
     build-depends:
         template-haskell
            >= 2.18       && < 2.22
-
-      , acts
-          ^>= 0.3.1.0
       , bifunctors
            >= 5.5.4      && < 5.7
       , code-page
@@ -144,10 +151,6 @@ library
           ^>= 3.3.7.0
       , filepath
            >= 1.4     && < 1.6
-      , generic-lens
-           >= 2.2     && < 2.3
-      , groups
-          ^>= 0.5.3
       , groups-generic
           ^>= 0.3.1.0
       , parallel
@@ -161,13 +164,58 @@ library
       , transformers
            >= 0.5.6.2 && < 0.7
 
+--executable convert-metafont
+--
+--    import:
+--        common
+--
+--    hs-source-dirs:
+--        src/metafont
+--
+--    default-language:
+--        Haskell2010
+--
+--    main-is:
+--        Main.hs
+--
+--    other-modules:
+--        Calligraphy.MetaFont.Convert
+--
+--    build-depends:
+--        diagrams-contrib,
+--        diagrams-lib,
+--        linear,
+--        parsec
+
+executable inspect
+
+    import:
+        common, extra
+
+    hs-source-dirs:
+        src/cusps/inspect
+
+    default-language:
+        Haskell2010
+
+    main-is:
+        Main.hs
+
+    other-modules:
+        Math.Interval.Abstract
+
+    build-depends:
+      brush-strokes,
+      data-reify
+          ^>= 0.6.3
+
 benchmark cusps
 
   import:
     common
 
   hs-source-dirs:
-    bench
+    src/cusps/bench
 
   main-is:
     Main.hs
diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs
new file mode 100644
index 0000000..36e793c
--- /dev/null
+++ b/brush-strokes/src/cusps/bench/Main.hs
@@ -0,0 +1,616 @@
+{-# LANGUAGE PolyKinds #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+module Main
+  ( main
+
+    -- Testing
+  , TestCase(..)
+  , testCases
+  , BrushStroke(..)
+  , brushStrokeFunctions
+  , eval
+  , mkVal, mkBox
+  , potentialCusp
+  , dEdsdcdt
+  )
+  where
+
+-- base
+import Control.Concurrent.MVar
+  ( newMVar )
+import Data.Coerce
+  ( coerce )
+import Data.Foldable
+  ( for_ )
+import Data.List
+  ( intercalate )
+import GHC.Exts
+  ( Proxy#, proxy# )
+import GHC.Generics
+  ( Generic )
+import GHC.TypeNats
+  ( type (-) )
+import Numeric
+  ( showFFloat )
+
+-- containers
+import Data.Sequence
+  ( Seq )
+import qualified Data.Sequence as Seq
+  ( index )
+import Data.Tree
+  ( foldTree )
+
+-- tree-view
+import Data.Tree.View
+  ( showTree )
+
+-- brush-strokes
+import Calligraphy.Brushes
+import Debug.Utils
+  ( logToFile )
+import Math.Algebra.Dual
+import Math.Bezier.Spline
+import Math.Bezier.Stroke
+import Math.Bezier.Stroke.EnvelopeEquation
+import Math.Differentiable
+import Math.Interval
+import Math.Linear
+import Math.Module
+import Math.Ring
+  ( Transcendental )
+
+--------------------------------------------------------------------------------
+
+main :: IO ()
+main = do
+  for_ testCases $ \ testCase@( TestCase { testName, testBrushStroke, testAlgorithmParams, testStartBoxes } ) -> do
+    let ( _, testStrokeFnI ) = brushStrokeFunctions testBrushStroke
+        ( newtTrees, ( dunno, sols ) ) = computeCusps testAlgorithmParams testStrokeFnI testStartBoxes
+        showedTrees = map ( uncurry showIntervalNewtonTree ) newtTrees
+        testHeader =
+          [ "", "Test case '" ++ testName ++ "':" ]
+    putStrLn $ unlines $
+      testHeader ++
+      map ( "    " ++ )
+        [ " #sols: " ++ show (length sols)
+        , "#dunno: " ++ show (length dunno)
+        , "#trees: " ++ show @Int (sum @_ @Int $ map (foldTree ( \ _ bs -> 1 + sum bs )) showedTrees)
+        , " dunno: " ++ show dunno
+        , "  sols: " ++ show sols
+        ]
+  --  logFileMVar <- newMVar "logs/fnData.log"
+  --  logToFile logFileMVar (unlines logLines)
+  --    `seq` return ()
+
+testCases :: [ TestCase ]
+testCases = benchCases
+  -- [ --trickyCusp2
+  --   ellipse "full" (0,1) pi $ defaultStartBoxes [ 2 ]
+  -- ]
+  -- ++
+{-
+    [ ellipse ( "(k1, k2) = " ++ show (k1, k2) ) (k1, k2) pi $ defaultStartBoxes [ 2 ]
+    | (k1, k2) <-
+      [(0.5,0.6), (0.55, 0.56)]
+    ] ++
+    [ ellipse ( "'(k1, k2) = " ++ show (k1, k2) ) (0,1) pi [ mkBox (k1 + zero, k2 + zero) 2 (zero, one) ]
+    | (k1, k2) <-
+      [(0.5,0.6), (0.55, 0.56)]
+    ]
+-}
+
+benchCases :: [ TestCase ]
+benchCases = [ ellipseTestCase "full" ( 0, 1 ) pi $ defaultStartBoxes [ 0 .. 3 ] ]
+
+--------------------------------------------------------------------------------
+
+data BrushStroke =
+  forall nbParams. ParamsCt nbParams =>
+  BrushStroke
+  { brush  :: !( Brush nbParams )
+  , stroke :: !( Point nbParams, Curve Open () ( Point nbParams ) )
+  }
+
+
+data TestCase =
+  TestCase
+  { testName            :: String
+  , testBrushStroke     :: BrushStroke
+  , testAlgorithmParams :: CuspAlgorithmParams
+  , testStartBoxes      :: [ Box ]
+  }
+
+brushStrokeFunctions
+  :: BrushStroke
+  -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () )
+     , 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
+brushStrokeFunctions ( BrushStroke { stroke = ( sp0, crv ), brush } ) =
+  getStrokeFunctions brush sp0 crv
+
+-- Utilities to use in GHCi to help debugging.
+
+eval
+  :: ( I i ( ℝ 1 ) -> Seq ( I i ( ℝ 1 ) -> StrokeDatum k i ) )
+  -> ( I i ( ℝ 1 ), Int, I i ( ℝ 1 ) )
+  -> StrokeDatum k i
+eval f ( t, i, s ) = ( f t `Seq.index` i ) s
+
+mkVal :: Double -> Int -> Double -> ( ℝ 1, Int, ℝ 1 )
+mkVal t i s = ( ℝ1 t, i, ℝ1 s )
+
+mkBox :: ( Double, Double ) -> Int -> ( Double, Double ) -> Box
+mkBox ( t_min, t_max ) i ( s_min, s_max ) =
+  ( 𝕀 ( ℝ1 t_min ) ( ℝ1 t_max ) , i, 𝕀 ( ℝ1 s_min ) ( ℝ1 s_max ) )
+
+zero, one :: Double
+zero = 1e-6
+one = 1 - zero
+{-# INLINE zero #-}
+{-# INLINE one #-}
+
+potentialCusp :: StrokeDatum 3 𝕀 -> Bool
+potentialCusp
+  ( StrokeDatum
+    { ee = D22 { _D22_v = 𝕀 ( ℝ1 ee_min ) ( ℝ1 ee_max ) }
+    , 𝛿E𝛿sdcdt = D12 { _D12_v = T ( 𝕀 ( ℝ2 vx_min vy_min ) ( ℝ2 vx_max vy_max ) )}
+    }
+  ) =  ee_min <= 0 && ee_max >= 0
+    && vx_min <= 0 && vx_max >= 0
+    && vy_min <= 0 && vy_max >= 0
+
+dEdsdcdt :: StrokeDatum k i -> D ( k - 2 ) ( I i ( ℝ 2 ) ) ( T ( I i ( ℝ 2 ) ) )
+dEdsdcdt ( StrokeDatum { 𝛿E𝛿sdcdt = v } ) = v
+
+{-
+let (f, fI) = testCaseStrokeFunctions trickyCusp2
+
+take 10 $ Data.List.sortOn ( \ ( _, ℝ1 e, v) -> abs e + norm v ) [ let { v = mkVal x 3 y; d = eval f v } in ( v, _D12_v $ ee d, _D0_v $ dEdsdcdt d ) | x <- [0.57,0.5701 .. 0.58], y <- [0.29,0.291..0.3] ]
+> ((ℝ1 0.5798800000000057,3,ℝ1 0.267980000000008),ℝ1 -2.8596965543670194e-4,V2 7.79559474412963e-2 2.0389671921293484e-2)
+
+potentialCusp $ eval fI $ mkBox (0.5798, 0.5799) 3 (0.26798, 0.26799)
+> True
+
+let nbPotentialSols b = let ( _newtTrees, ( dunno, sols ) ) = intervalNewtonGSFrom NoPreconditioning 1e-7 fI b in length dunno + length sols
+
+nbPotentialSols $ mkBox (0.5798, 0.5799) 3 (0.26798, 0.26799)
+1
+
+nbPotentialSols $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799)
+0
+
+let showTrees b = map ( uncurry showIntervalNewtonTree ) $ fst $ intervalNewtonGSFrom NoPreconditioning 1e-7 fI b
+
+putStrLn $ unlines $ map Data.Tree.View.showTree $ showTrees $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799)
+
+([ℝ1 0.5798, ℝ1 0.675],3,[ℝ1 0.26798, ℝ1 0.26799]) (area 0.000001) N []
+ └─ ([ℝ1 0.5973000285624527, ℝ1 0.6750000000000002],3,[ℝ1 0.26798, ℝ1 0.26799000000000006]) (area 0.000001) NoSolution "ee" ([ℝ1 0.5973000285624527, ℝ1 0.6750000000000002],3,[ℝ1 0.26798, ℝ1 0.26799000000000006])
+
+eval fI $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799)
+> D12 {_D12_v = T[ℝ2 -10088.6674944889 -3281.3820867312834, ℝ2 4124.668381545453 4524.807156085763], _D12_dx = TT[ℝ2 -173746.97965005718 -33281.18494907289, ℝ2 298.2609121556852 23639.772884799597], _D12_dy = TT[ℝ2 -18454.27716258352 -28337.509817580823, ℝ2 1163.6949532017436 -13936.383137525536]}}
+i.e.
+> f = [ℝ2 -10088.6674944889 -3281.3820867312834, ℝ2 4124.668381545453 4524.807156085763]
+> f_t = [ℝ2 -173746.97965005718 -33281.18494907289, ℝ2 298.2609121556852 23639.772884799597]
+> f_s = [ℝ2 -18454.27716258352 -28337.509817580823, ℝ2 1163.6949532017436 -13936.383137525536]
+
+(f, fI) = testCaseStrokeFunctions trickyCusp2
+t = 𝕀 (ℝ1 0.5798) (ℝ1 0.675)
+s = 𝕀 (ℝ1 0.26798) (ℝ1 0.26799)
+t_mid = 0.5 * ( 0.5798 + 0.675 )
+s_mid = 0.5 * ( 0.26798 + 0.26799 )
+D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) = dEdsdcdt $ eval fI (t, 3, s)
+t' = coerce ( (-) @( 𝕀 Double ) ) t ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1
+s' = coerce ( (-) @( 𝕀 Double ) ) s ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1
+a = ( f_t, f_s )
+b = negV2 $ singleton $ midV2 f
+[((t2', s2'), isContr)] = gaussSeidel a b (t', s')
+t2 = coerce ( (+) @( 𝕀 Double ) ) t2' ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1
+s2 = coerce ( (+) @( 𝕀 Double ) ) s2' ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1
+
+t2
+> [ℝ1 0.6102365832093095, ℝ1 0.6750000000000002]
+s2
+> [ℝ1 0.26798, ℝ1 0.26799000000000006]
+
+
+let ( 𝕀 ( ℝ2 a11_lo a21_lo ) ( ℝ2 a11_hi a21_hi ), 𝕀 ( ℝ2 a12_lo a22_lo ) ( ℝ2 a12_hi a22_hi ) ) = a
+let ( 𝕀 ( ℝ2 b1_lo b2_lo ) ( ℝ2 b1_hi b2_hi ) ) = b
+let ( 𝕀 ( ℝ1 x1_lo ) ( ℝ1 x1_hi ), 𝕀 ( ℝ1 x2_lo ) ( ℝ1 x2_hi ) ) = ( t', s' )
+
+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
+
+( b1 - a12 * x2 )
+> [2981.90728508591, 2982.0918278575364]
+
+extendedRecip a11
+
+-----------------------
+
+fI k rot = snd $ testCaseStrokeFunctions $ ellipse k rot
+d k rot = width $ _D22_v $ ee $ eval (fI k rot) $ mkBox (1e-6, 1-1e-6) 3 (1e-6, 1-1e-6)
+
+-----------------------
+
+(f, fI) = testCaseStrokeFunctions $ ellipse 1 pi
+t = 𝕀 (ℝ1 0.001) (ℝ1 0.099)
+s = 𝕀 (ℝ1 0.001) (ℝ1 0.999)
+t_mid = 0.5 * ( 0.001 + 0.099 )
+s_mid = 0.5 * ( 0.001 + 0.999 )
+D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) = dEdsdcdt $ eval fI (t, 3, s)
+t' = coerce ( (-) @( 𝕀 Double ) ) t ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1
+s' = coerce ( (-) @( 𝕀 Double ) ) s ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1
+a = ( f_t, f_s )
+b = negV2 $ singleton $ midV2 f
+[((t2', s2'), isContr)] = gaussSeidel a b (t', s')
+t2 = coerce ( (+) @( 𝕀 Double ) ) t2' ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1
+s2 = coerce ( (+) @( 𝕀 Double ) ) s2' ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1
+
+> t = [0.001, 0.099]
+> t' = [-0.049, 0.049]
+
+(f, fI) = testCaseStrokeFunctions $ ellipse 0.1 pi
+t = 𝕀 (ℝ1 0.001) (ℝ1 0.999)
+s = 𝕀 (ℝ1 0.001) (ℝ1 0.999)
+t_mid = 0.5 * ( 0.001 + 0.999 )
+s_mid = 0.5 * ( 0.001 + 0.999 )
+D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) = dEdsdcdt $ eval fI (t, 3, s)
+t' = coerce ( (-) @( 𝕀 Double ) ) t ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1
+s' = coerce ( (-) @( 𝕀 Double ) ) s ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1
+a = ( f_t, f_s )
+b = negV2 $ singleton $ midV2 f
+[((t2', s2'), isContr)] = gaussSeidel a b (t', s')
+t2 = coerce ( (+) @( 𝕀 Double ) ) t2' ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1
+s2 = coerce ( (+) @( 𝕀 Double ) ) s2' ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1
+
+> t = [0.001, 0.999]
+> t' = [-0.499, 0.499]
+
+-}
+
+width :: 𝕀ℝ 1 -> Double
+width (𝕀 (ℝ1 lo) (ℝ1 hi)) = hi - lo
+
+negV2 :: 𝕀ℝ 2 -> 𝕀ℝ 2
+negV2 ( 𝕀 ( ℝ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 )
+
+midV2 :: 𝕀ℝ 2 -> ℝ 2
+midV2 ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) =
+  ℝ2 ( 0.5 * ( x_lo + x_hi ) ) ( 0.5 * ( y_lo + y_hi ) )
+
+
+logLines :: [ String ]
+logLines =
+  [ "E, dE/ds * dc/dt"
+  , "{" ++
+    (intercalate ","
+     [ "{" ++ showD t ++ "," ++ showD s ++ ",{" ++ intercalate "," vals ++ "}}"
+     | t <- [ 0.5484, 0.5484 + 0.00001 .. 0.5488 ]
+     , s <- [ 0.5479, 0.5479 + 0.00001 .. 0.5483 ]
+     , let StrokeDatum
+             { ee = D22 ee _ _ _ _ _
+             , 𝛿E𝛿sdcdt = D12 (T f) (T (T f_t)) (T (T f_s))
+             } = (curvesI (singleton (ℝ1 t)) `Seq.index` i) (singleton (ℝ1 s))
+           ℝ2 vx vy = midPoint2 f
+           --ℝ2 vx_t vy_t = midPoint2 f_t
+           --ℝ2 vx_s vy_s = midPoint2 f_s
+           vals = [ showD ( midPoint ee )
+                  , "{" ++ showD vx ++ "," ++ showD vy ++ "}"
+                --  , "{" ++ showD vx_t ++ "," ++ showD vy_t ++ "}"
+                --  , "{" ++ showD vx_s ++ "," ++ showD vy_s ++ "}"
+                  ]
+     ]
+    ) ++ "}"
+  ]
+  where
+    i = 2
+    ( curves, curvesI ) = brushStrokeFunctions $ ellipseBrushStroke ( 0, 1 ) pi
+    midPoint (𝕀 (ℝ1 lo) (ℝ1 hi)) = 0.5 * ( lo + hi )
+    midPoint2 (𝕀 (ℝ2 lo_x lo_y) (ℝ2 hi_x hi_y))
+      = ℝ2 ( 0.5 * ( lo_x + hi_x ) ) ( 0.5 * ( lo_y + hi_y ) )
+
+
+-- t = 0.5486102, s = 0.5480951
+
+
+bloo =
+  [ ( e * e + vx * vx + vy * vy, ( t, s ) )
+  | t <- [ 0.548609, 0.548609 + 0.0000001 .. 0.54862 ]
+  , s <- [ 0.548094, 0.548094 + 0.0000001 .. 0.548096 ]
+  , let StrokeDatum
+          { ee = D22 ee _ _ _ _ _
+          , 𝛿E𝛿sdcdt = D12 (T f) _ _
+          } = (curvesI (singleton (ℝ1 t)) `Seq.index` i) (singleton (ℝ1 s))
+        e = midPoint ee
+        ℝ2 vx vy = midPoint2 f
+        vals = [ showD ( midPoint ee )
+               , "{" ++ showD vx ++ "," ++ showD vy ++ "}"
+               ]
+  ]
+  where
+    i = 2
+    ( curves, curvesI ) = brushStrokeFunctions $ ellipseBrushStroke ( 0, 1 ) pi
+    midPoint (𝕀 (ℝ1 lo) (ℝ1 hi)) = 0.5 * ( lo + hi )
+    midPoint2 (𝕀 (ℝ2 lo_x lo_y) (ℝ2 hi_x hi_y))
+      = ℝ2 ( 0.5 * ( lo_x + hi_x ) ) ( 0.5 * ( lo_y + hi_y ) )
+
+{-
+
+(f, fI) = brushStrokeFunctions $ ellipseBrushStroke ( 0, 1 ) pi
+
+
+
+
+_D22_v $ ee $ eval fI $ mkBox (0.5, 0.6) 2 (0,1)
+> [ℝ1 -9531.427315889887, ℝ1 10135.074304695485]
+
+minimum $ map inf $ [ _D22_v $ ee $ eval fI $ mkBox (t, t + 0.01) 2 (0,1) | t <- [ 0.5, 0.51 .. 0.59 ] ]
+> ℝ1 -5718.905635365308
+
+maximum $ map sup $ [ _D22_v $ ee $ eval fI $ mkBox (t, t + 0.01) 2 (0,1) | t <- [ 0.5, 0.51 .. 0.59 ] ]
+> ℝ1 5099.008191092755
+
+minimum $ map inf $ [ _D22_v $ ee $ eval fI $ mkBox (t, t) 2 (s, s) | t <- [ 0.5, 0.501 .. 0.6 ], s <- [ zero, zero + 0.001 .. one ] ]
+> ℝ1 -675.9595496147458
+
+maximum $ map sup $ [ _D22_v $ ee $ eval fI $ mkBox (t, t) 2 (s, s) | t <- [ 0.5, 0.501 .. 0.6 ], s <- [ zero, zero + 0.001 .. one ] ]
+> ℝ1 2401.9644509525997
+
+_D12_v $ dEdsdcdt $ eval fI $ mkBox (0.5, 0.6) 2 (0,1)
+> T[ℝ2 -1.7300637136531524e7 -1.262824151868635e7, ℝ2 1.632868898735965e7 1.1869759856947478e7]
+
+minimum [ _x $ inf $ unT $ _D12_v $ dEdsdcdt $ eval fI $ mkBox (t, t + 0.01) 2 (0,1) | t <- [ 0.5, 0.51 .. 0.59 ] ]
+-5606615.948203902
+
+maximum [ _x $ sup $ unT $ _D12_v $ dEdsdcdt $ eval fI $ mkBox (t, t + 0.01) 2 (0,1) | t <- [ 0.5, 0.51 .. 0.59 ] ]
+4340858.832347277
+
+minimum [ _x $ inf $ unT $ _D12_v $ dEdsdcdt $ eval fI $ mkBox (t, t) 2 (s, s) | t <- [ 0.5, 0.501 .. 0.6 ], s <- [ zero, zero + 0.001 .. one ] ]
+-1785730.2396688666
+
+maximum [ _x $ sup $ unT $ _D12_v $ dEdsdcdt $ eval fI $ mkBox (t, t) 2 (s, s) | t <- [ 0.5, 0.501 .. 0.6 ], s <- [ zero, zero + 0.001 .. one ] ]
+974842.6547409865
+
+maximum [ _y $ sup $ unT $ _D12_v $ dEdsdcdt $ eval fI $ mkBox (t, t) 2 (s, s) | t <- [ 0.5, 0.501 .. 0.6 ], s <- [ zero, zero + 0.001 .. one ] ]
+845211.4833711373
+
+
+
+let showTrees b = map ( uncurry showIntervalNewtonTree ) $ fst $ intervalNewtonGSFrom NoPreconditioning 1e-7 fI b
+
+putStrLn $ unlines $ map Data.Tree.View.showTree $ showTrees $ mkBox (0.5486101933245248, 0.5486102071493622) 2 (0.548095036738487, 0.5480952)
+
+
+
+([ℝ1 0.5486101960904595, ℝ1 0.5486102071493623],2,[ℝ1 0.5480950771755867, ℝ1 0.5480952000000001])
+
+-}
+
+_x ( ℝ2 x _ ) = x
+_y ( ℝ2 _ y ) = y
+
+showD :: Double -> String
+showD float = showFFloat (Just 6) float ""
+
+--------------------------------------------------------------------------------
+
+ellipseTestCase :: String -> ( Double, Double ) -> Double -> [ Box ] -> TestCase
+ellipseTestCase str k0k1 rot startBoxes =
+  TestCase
+    { testName = "ellipse (" ++ str ++ ")"
+    , testBrushStroke = ellipseBrushStroke k0k1 rot
+    , testAlgorithmParams =
+        CuspAlgorithmParams
+         { preconditioning = InverseMidJacobian
+         , maxWidth = 1e-7
+         }
+    , testStartBoxes = startBoxes
+    }
+
+ellipseBrushStroke :: ( Double, Double ) -> Double -> BrushStroke
+ellipseBrushStroke ( k0, k1 ) rot =
+  BrushStroke
+    { brush = ellipseBrush
+    , stroke = ( p0, LineTo ( NextPoint p1 ) () ) }
+  where
+    mkPt x y w h phi =
+      Point
+        { pointCoords = ℝ2 x y
+        , pointParams = Params $ ℝ3 w h phi
+        }
+    l k = lerp @( T Double ) k
+    p k = mkPt ( l k 0 100 ) 0 ( l k 10 15 ) ( l k 25 40 ) ( l k 0 rot )
+    p0 = p k0
+    p1 = p k1
+
+trickyCusp2TestCase :: TestCase
+trickyCusp2TestCase =
+  TestCase
+    { testName = "trickyCusp2"
+    , testBrushStroke = trickyCusp2BrushStroke
+    , testAlgorithmParams =
+        CuspAlgorithmParams
+         { preconditioning = InverseMidJacobian
+         , maxWidth = 1e-7
+         }
+    , testStartBoxes = defaultStartBoxes [ 0 .. 3 ]
+    }
+
+trickyCusp2BrushStroke :: BrushStroke
+trickyCusp2BrushStroke =
+  BrushStroke
+    { brush = circleBrush
+    , stroke = ( p0, Bezier3To p1 p2 ( NextPoint p3 ) () )
+    }
+  where
+    mkPt x y =
+      Point
+        { pointCoords = ℝ2 x y
+        , pointParams = Params $ ℝ1 5.0
+        }
+    p0 = mkPt 5e+1 -5e+1
+    p1 = mkPt -7.72994362904069e+1 -3.124468786098509e+1
+    p2 = mkPt -5.1505430313958364e+1 -3.9826386521527986e+1
+    p3 = mkPt -5e+1 -5e+1
+
+--------------------------------------------------------------------------------
+
+type ParamsCt nbParams
+  = ( Show ( ℝ nbParams )
+    , HasChainRule Double 2 ( ℝ nbParams )
+    , HasChainRule ( 𝕀 Double ) 3 ( 𝕀 ( ℝ nbParams ) )
+    , Applicative ( D 2 ( ℝ nbParams ) )
+    , Applicative ( D 3 ( ℝ nbParams ) )
+    , Traversable ( D 2 ( ℝ nbParams ) )
+    , Traversable ( D 3 ( ℝ nbParams ) )
+    , Representable Double ( ℝ nbParams )
+    , Module Double ( T ( ℝ nbParams ) )
+    , Module ( 𝕀 Double ) ( T ( 𝕀 ( ℝ nbParams ) ) )
+    , Module ( D 2 ( ℝ nbParams ) Double ) ( D 2 ( ℝ nbParams ) ( ℝ 2 ) )
+    , Module ( D 3 ( ℝ nbParams ) ( 𝕀 Double ) ) ( D 3 ( ℝ nbParams ) ( 𝕀 ( ℝ 2 ) ) )
+    , Transcendental ( D 2 ( ℝ nbParams ) Double )
+    , Transcendental ( D 3 ( ℝ nbParams ) ( 𝕀 Double ) )
+    )
+
+newtype Params nbParams = Params { getParams :: ( ℝ nbParams ) }
+deriving newtype instance Show ( ℝ nbParams ) => Show ( Params nbParams )
+
+data Point nbParams =
+  Point
+    { pointCoords :: !( ℝ 2 )
+    , pointParams :: !( Params nbParams ) }
+  deriving stock Generic
+deriving stock instance Show ( ℝ nbParams ) => Show ( Point nbParams )
+
+data CuspAlgorithmParams =
+  CuspAlgorithmParams
+    { preconditioning :: !Preconditioner
+    , maxWidth :: !Double
+    }
+  deriving stock Show
+
+type Brush nbParams
+  = forall {t} k (i :: t)
+  .  DiffInterp k i ( ℝ nbParams )
+  => Proxy# i
+  -> ( forall a. a -> I i a )
+  -> C k ( I i ( ℝ nbParams ) )
+         ( Spline Closed () ( I i ( ℝ 2 ) ) )
+
+getStrokeFunctions
+  :: forall nbParams
+  .  ParamsCt nbParams
+  => Brush nbParams
+      -- ^ brush shape
+  -> Point nbParams
+      -- ^ start point
+  -> Curve Open () ( Point nbParams )
+      -- ^ curve points
+  -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () )
+     , 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
+getStrokeFunctions brush sp0 crv =
+  let
+    usedParams :: C 2 ( ℝ 1 ) ( ℝ nbParams )
+    path :: C 2 ( ℝ 1 ) ( ℝ 2 )
+    ( path, usedParams ) =
+      pathAndUsedParams @2 @() coerce id ( getParams . pointParams )
+        sp0 crv
+    usedParamsI :: C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ nbParams )
+    pathI :: C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 )
+    ( pathI, usedParamsI ) =
+      pathAndUsedParams @3 @𝕀 coerce singleton ( getParams . pointParams )
+        sp0 crv
+  in ( brushStrokeData @2 @( ℝ nbParams ) coerce coerce
+        path usedParams $
+        brush @2 @() proxy# id
+     , brushStrokeData @3 @( ℝ nbParams ) coerce coerce
+        pathI usedParamsI $
+        brush @3 @𝕀 proxy# singleton )
+{-# INLINEABLE getStrokeFunctions #-}
+
+computeCusps
+  :: CuspAlgorithmParams
+  -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
+  -> [ Box ]
+  -> ( [ ( Box, IntervalNewtonTree Box ) ], ( [ Box ], [ Box ] ) )
+computeCusps params eqs startBoxes =
+  foldMap
+    ( intervalNewtonGSFrom ( preconditioning params ) ( maxWidth params ) eqs )
+    startBoxes
+
+defaultStartBoxes :: [ Int ] -> [ Box ]
+defaultStartBoxes is =
+  [ mkBox (zero, one) i (zero, one) | i <- is ]
+
+getR1 (ℝ1 u) = u
+
+{-
+
+(f, fI) = brushStrokeFunctions $ ellipseBrushStroke (0,1) pi
+nbPotentialSols box = let ( _newtTrees, ( dunno, sols ) ) = intervalNewtonGSFrom NoPreconditioning 1e-7 fI box in length dunno + length sols
+showTrees box = map ( uncurry showIntervalNewtonTree ) $ fst $ intervalNewtonGSFrom NoPreconditioning 1e-7 fI box
+
+(t, i, s) = mkBox (0.548610200176363, 0.5486102071493623) 2 (0.5480950215354709, 0.5480952)
+putStrLn $ unlines $ map Data.Tree.View.showTree $ showTrees (t,i,s)
+
+
+t_mid = 0.5 * ( getR1 ( inf t ) + getR1 ( sup t ) )
+s_mid = 0.5 * ( getR1 ( inf s ) + getR1 ( sup s ) )
+D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) = dEdsdcdt $ eval fI (t, i, s)
+t' = coerce ( (-) @( 𝕀 Double ) ) t ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1
+s' = coerce ( (-) @( 𝕀 Double ) ) s ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1
+a = ( f_t, f_s )
+b = negV2 $ singleton $ midV2 f
+[((t2', s2'), isContr)] = gaussSeidel a b (t', s')
+t2 = coerce ( (+) @( 𝕀 Double ) ) t2' ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1
+s2 = coerce ( (+) @( 𝕀 Double ) ) s2' ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1
+
+t2
+> [ℝ1 0.548610200176363, ℝ1 0.5486102071493624]
+s2
+> [ℝ1 0.5480950911334656, ℝ1 0.5480952000000001]
+
+mkBox (0.548610200176363, 0.5486102071493624) i (0.5480950911334656, 0.5480952000000001)
+
+t inf (no change)
+
+t sup (no change)
+
+s_inf:
+0.5480950215354709
+0.5480950911334656
+
+s_sup (no change)
+
+ghci> potentialCusp $ eval fI $ mkBox (0.548610200176363, 0.5486102071493623) 2 (0.54809502, 0.5480952)
+True
+ghci> potentialCusp $ eval fI $ mkBox (0.548610200176363, 0.5486102071493623) 2 (0.54809503, 0.5480952)
+False
+
+
+
+let ( 𝕀 ( ℝ2 a11_lo a21_lo ) ( ℝ2 a11_hi a21_hi ), 𝕀 ( ℝ2 a12_lo a22_lo ) ( ℝ2 a12_hi a22_hi ) ) = a
+let ( 𝕀 ( ℝ2 b1_lo b2_lo ) ( ℝ2 b1_hi b2_hi ) ) = b
+let ( 𝕀 ( ℝ1 x1_lo ) ( ℝ1 x1_hi ), 𝕀 ( ℝ1 x2_lo ) ( ℝ1 x2_hi ) ) = ( t', s' )
+
+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
+
+-}
\ No newline at end of file
diff --git a/src/cusps/Main.hs b/brush-strokes/src/cusps/inspect/Main.hs
similarity index 99%
rename from src/cusps/Main.hs
rename to brush-strokes/src/cusps/inspect/Main.hs
index ee71b18..a876c12 100644
--- a/src/cusps/Main.hs
+++ b/brush-strokes/src/cusps/inspect/Main.hs
@@ -32,7 +32,7 @@ import qualified Data.Sequence as Seq
 import Data.Generics.Product.Typed
   ( HasType )
 
--- MetaBrush
+-- brush-strokes
 import Math.Algebra.Dual
 import qualified Math.Bezier.Quadratic as Quadratic
 import qualified Math.Bezier.Cubic     as Cubic
@@ -58,6 +58,7 @@ import Math.Ring
   , ifThenElse
   )
 
+-- brush-strokes:inspect-cusps
 import Math.Interval.Abstract
 
 --------------------------------------------------------------------------------
diff --git a/src/cusps/Math/Interval/Abstract.hs b/brush-strokes/src/cusps/inspect/Math/Interval/Abstract.hs
similarity index 99%
rename from src/cusps/Math/Interval/Abstract.hs
rename to brush-strokes/src/cusps/inspect/Math/Interval/Abstract.hs
index 0d54cb4..2d4baa6 100644
--- a/src/cusps/Math/Interval/Abstract.hs
+++ b/brush-strokes/src/cusps/inspect/Math/Interval/Abstract.hs
@@ -43,7 +43,7 @@ import qualified Data.Map.Strict as Map
 import Data.Group
   ( Group(..) )
 
--- MetaBrush
+-- brush-strokes
 import Math.Algebra.Dual
 import qualified Math.Bezier.Quadratic as Quadratic
 import qualified Math.Bezier.Cubic     as Cubic
diff --git a/brush-strokes/src/Calligraphy/Brushes.hs b/brush-strokes/src/lib/Calligraphy/Brushes.hs
similarity index 100%
rename from brush-strokes/src/Calligraphy/Brushes.hs
rename to brush-strokes/src/lib/Calligraphy/Brushes.hs
diff --git a/brush-strokes/src/Debug/Utils.hs b/brush-strokes/src/lib/Debug/Utils.hs
similarity index 98%
rename from brush-strokes/src/Debug/Utils.hs
rename to brush-strokes/src/lib/Debug/Utils.hs
index 1687a61..0b32b91 100644
--- a/brush-strokes/src/Debug/Utils.hs
+++ b/brush-strokes/src/lib/Debug/Utils.hs
@@ -66,3 +66,5 @@ logToFile logFileMVar logContents =
       FilePath.takeDirectory logFile
     appendFile logFile logContentsWithHeader
     return logFile
+{-# NOINLINE logToFile #-}
+
diff --git a/brush-strokes/src/Math/Algebra/Dual.hs b/brush-strokes/src/lib/Math/Algebra/Dual.hs
similarity index 100%
rename from brush-strokes/src/Math/Algebra/Dual.hs
rename to brush-strokes/src/lib/Math/Algebra/Dual.hs
diff --git a/brush-strokes/src/Math/Algebra/Dual/Internal.hs b/brush-strokes/src/lib/Math/Algebra/Dual/Internal.hs
similarity index 95%
rename from brush-strokes/src/Math/Algebra/Dual/Internal.hs
rename to brush-strokes/src/lib/Math/Algebra/Dual/Internal.hs
index 5b2a1cd..181ec9f 100644
--- a/brush-strokes/src/Math/Algebra/Dual/Internal.hs
+++ b/brush-strokes/src/lib/Math/Algebra/Dual/Internal.hs
@@ -60,8 +60,8 @@ newtype D𝔸0 v = D0 { _D0_v :: v }
 
 -- | \( \mathbb{Z}[x]/(x)^2 \).
 data D1𝔸1 v =
-  D11 { _D11_v :: !v
-      , _D11_dx :: !( T v )
+  D11 { _D11_v :: v
+      , _D11_dx :: ( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
   deriving anyclass NFData
@@ -70,9 +70,9 @@ data D1𝔸1 v =
 
 -- | \( \mathbb{Z}[x]/(x)^3 \).
 data D2𝔸1 v =
-  D21 { _D21_v :: !v
-      , _D21_dx :: !( T v )
-      , _D21_dxdx :: !( T v )
+  D21 { _D21_v :: v
+      , _D21_dx :: ( T v )
+      , _D21_dxdx :: ( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
   deriving anyclass NFData
@@ -81,10 +81,10 @@ data D2𝔸1 v =
 
 -- | \( \mathbb{Z}[x]/(x)^4 \).
 data D3𝔸1 v =
-  D31 { _D31_v :: !v
-      , _D31_dx :: !( T v )
-      , _D31_dxdx :: !( T v )
-      , _D31_dxdxdx :: !( T v )
+  D31 { _D31_v :: v
+      , _D31_dx :: ( T v )
+      , _D31_dxdx :: ( T v )
+      , _D31_dxdxdx :: ( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
   deriving anyclass NFData
@@ -93,8 +93,8 @@ data D3𝔸1 v =
 
 -- | \( \mathbb{Z}[x, y]/(x, y)^2 \).
 data D1𝔸2 v =
-  D12 { _D12_v :: !v
-      , _D12_dx, _D12_dy :: !( T v )
+  D12 { _D12_v :: v
+      , _D12_dx, _D12_dy :: ( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
   deriving anyclass NFData
@@ -103,9 +103,9 @@ data D1𝔸2 v =
 
 -- | \( \mathbb{Z}[x, y]/(x, y)^3 \).
 data D2𝔸2 v =
-  D22 { _D22_v :: !v
-      , _D22_dx, _D22_dy :: !( T v )
-      , _D22_dxdx, _D22_dxdy, _D22_dydy :: !( T v )
+  D22 { _D22_v :: v
+      , _D22_dx, _D22_dy :: ( T v )
+      , _D22_dxdx, _D22_dxdy, _D22_dydy :: ( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
   deriving anyclass NFData
@@ -114,10 +114,10 @@ data D2𝔸2 v =
 
 -- | \( \mathbb{Z}[x, y]/(x, y)^4 \).
 data D3𝔸2 v =
-  D32 { _D32_v :: !v
-      , _D32_dx, _D32_dy :: !( T v )
-      , _D32_dxdx, _D32_dxdy, _D32_dydy :: !( T v )
-      , _D32_dxdxdx, _D32_dxdxdy, _D32_dxdydy, _D32_dydydy :: !( T v )
+  D32 { _D32_v :: v
+      , _D32_dx, _D32_dy :: ( T v )
+      , _D32_dxdx, _D32_dxdy, _D32_dydy :: ( T v )
+      , _D32_dxdxdx, _D32_dxdxdy, _D32_dxdydy, _D32_dydydy :: ( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
   deriving anyclass NFData
@@ -126,8 +126,8 @@ data D3𝔸2 v =
 
 -- | \( \mathbb{Z}[x, y, z]/(x, y, z)^2 \).
 data D1𝔸3 v =
-  D13 { _D13_v :: !v
-      , _D13_dx, _D13_dy, _D13_dz :: !( T v )
+  D13 { _D13_v :: v
+      , _D13_dx, _D13_dy, _D13_dz :: ( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
   deriving anyclass NFData
@@ -136,8 +136,8 @@ data D1𝔸3 v =
 
 -- | \( \mathbb{Z}[x, y, z]/(x, y, z)^3 \).
 data D2𝔸3 v =
-  D23 { _D23_v :: !v
-      , _D23_dx, _D23_dy, _D23_dz :: !( T v )
+  D23 { _D23_v :: v
+      , _D23_dx, _D23_dy, _D23_dz :: ( T v )
       , _D23_dxdx, _D23_dxdy, _D23_dydy, _D23_dxdz, _D23_dydz, _D23_dzdz :: !( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
@@ -147,8 +147,8 @@ data D2𝔸3 v =
 
 -- | \( \mathbb{Z}[x, y, z]/(x, y, z)^4 \).
 data D3𝔸3 v =
-  D33 { _D33_v :: !v
-      , _D33_dx, _D33_dy, _D33_dz :: !( T v )
+  D33 { _D33_v :: v
+      , _D33_dx, _D33_dy, _D33_dz :: ( T v )
       , _D33_dxdx, _D33_dxdy, _D33_dydy, _D33_dxdz, _D33_dydz, _D33_dzdz :: !( T v )
       , _D33_dxdxdx, _D33_dxdxdy, _D33_dxdydy, _D33_dydydy
          , _D33_dxdxdz, _D33_dxdydz, _D33_dxdzdz, _D33_dydydz, _D33_dydzdz, _D33_dzdzdz :: !( T v )
@@ -160,8 +160,8 @@ data D3𝔸3 v =
 
 -- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^2 \).
 data D1𝔸4 v =
-  D14 { _D14_v :: !v
-      , _D14_dx, _D14_dy, _D14_dz, _D14_dw :: !( T v )
+  D14 { _D14_v :: v
+      , _D14_dx, _D14_dy, _D14_dz, _D14_dw :: ( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
   deriving anyclass NFData
@@ -170,10 +170,10 @@ data D1𝔸4 v =
 
 -- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^3 \).
 data D2𝔸4 v =
-  D24 { _D24_v :: !v
-      , _D24_dx, _D24_dy, _D24_dz, _D24_dw :: !( T v )
+  D24 { _D24_v :: v
+      , _D24_dx, _D24_dy, _D24_dz, _D24_dw :: ( T v )
       , _D24_dxdx, _D24_dxdy, _D24_dydy, _D24_dxdz
-         , _D24_dydz, _D24_dzdz, _D24_dxdw, _D24_dydw, _D24_dzdw, _D24_dwdw :: !( T v )
+         , _D24_dydz, _D24_dzdz, _D24_dxdw, _D24_dydw, _D24_dzdw, _D24_dwdw :: ( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
   deriving anyclass NFData
@@ -182,14 +182,14 @@ data D2𝔸4 v =
 
 -- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^4 \).
 data D3𝔸4 v =
-  D34 { _D34_v :: !v
-      , _D34_dx, _D34_dy, _D34_dz, _D34_dw :: !( T v )
+  D34 { _D34_v :: v
+      , _D34_dx, _D34_dy, _D34_dz, _D34_dw :: ( T v )
       , _D34_dxdx, _D34_dxdy, _D34_dydy, _D34_dxdz, _D34_dydz, _D34_dzdz
-        , _D34_dxdw, _D34_dydw, _D34_dzdw, _D34_dwdw :: !( T v )
+        , _D34_dxdw, _D34_dydw, _D34_dzdw, _D34_dwdw :: ( T v )
       , _D34_dxdxdx, _D34_dxdxdy, _D34_dxdydy, _D34_dydydy,
          _D34_dxdxdz, _D34_dxdydz, _D34_dxdzdz, _D34_dydzdz, _D34_dydydz, _D34_dzdzdz,
          _D34_dxdxdw, _D34_dxdydw, _D34_dydydw, _D34_dxdzdw, _D34_dydzdw, _D34_dzdzdw,
-         _D34_dxdwdw, _D34_dydwdw, _D34_dzdwdw, _D34_dwdwdw :: !( T v )
+         _D34_dxdwdw, _D34_dydwdw, _D34_dzdwdw, _D34_dwdwdw :: ( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
   deriving anyclass NFData
diff --git a/brush-strokes/src/Math/Bezier/Cubic.hs b/brush-strokes/src/lib/Math/Bezier/Cubic.hs
similarity index 93%
rename from brush-strokes/src/Math/Bezier/Cubic.hs
rename to brush-strokes/src/lib/Math/Bezier/Cubic.hs
index 173a01e..ba4fc6d 100644
--- a/brush-strokes/src/Math/Bezier/Cubic.hs
+++ b/brush-strokes/src/lib/Math/Bezier/Cubic.hs
@@ -8,7 +8,7 @@ module Math.Bezier.Cubic
   , bezier, bezier', bezier'', bezier'''
   , derivative
   , curvature, squaredCurvature, signedCurvature
-  , subdivide
+  , subdivide, restrict
   , ddist, closestPoint
   , drag, selfIntersectionParameters
   , extrema
@@ -180,6 +180,17 @@ subdivide ( Bezier {..} ) t = ( Bezier p0 q1 q2 pt, Bezier pt r1 r2 p3 )
     pt = lerp @v t q2 r1
 {-# INLINEABLE subdivide #-}
 
+-- | Restrict a cubic Bézier curve to a sub-interval, re-parametrising
+-- to \( [0,1] \).
+restrict :: forall v r p. ( Torsor v p, Ring.Field r, Module r v ) => Bezier p -> ( r , r ) -> Bezier p
+restrict bez ( a, b ) = fst $ ( flip ( subdivide @v ) b' ) $ snd $ subdivide @v bez a
+  where
+    b' = ( b Ring.- a ) Ring./ ( Ring.fromInteger 1 Ring.- a )
+  -- TODO: this could be made more efficient.
+  -- See e.g. "https://math.stackexchange.com/questions/4172835/cubic-b%C3%A9zier-spline-multiple-split"
+  -- or the paper "On the numerical condition of Bernstein-Bézier subdivision process".
+{-# INLINEABLE restrict #-}
+
 -- | Polynomial coefficients of the derivative of the distance to a cubic Bézier curve.
 ddist :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> p -> [ r ]
 ddist ( Bezier {..} ) c = [ a5, a4, a3, a2, a1, a0 ]
diff --git a/brush-strokes/src/Math/Bezier/Cubic/Fit.hs b/brush-strokes/src/lib/Math/Bezier/Cubic/Fit.hs
similarity index 100%
rename from brush-strokes/src/Math/Bezier/Cubic/Fit.hs
rename to brush-strokes/src/lib/Math/Bezier/Cubic/Fit.hs
diff --git a/brush-strokes/src/Math/Bezier/Quadratic.hs b/brush-strokes/src/lib/Math/Bezier/Quadratic.hs
similarity index 90%
rename from brush-strokes/src/Math/Bezier/Quadratic.hs
rename to brush-strokes/src/lib/Math/Bezier/Quadratic.hs
index 8f5c8d5..4846036 100644
--- a/brush-strokes/src/Math/Bezier/Quadratic.hs
+++ b/brush-strokes/src/lib/Math/Bezier/Quadratic.hs
@@ -6,7 +6,7 @@ module Math.Bezier.Quadratic
   ( Bezier(..)
   , bezier, bezier', bezier''
   , curvature, squaredCurvature, signedCurvature
-  , subdivide
+  , subdivide, restrict
   , ddist, closestPoint
   , interpolate
   , extrema
@@ -141,6 +141,17 @@ subdivide ( Bezier {..} ) t = ( Bezier p0 q1 pt, Bezier pt r1 p2 )
     pt = lerp @v t q1 r1
 {-# INLINEABLE subdivide #-}
 
+-- | Restrict a quadratic Bézier curve to a sub-interval, re-parametrising
+-- to \( [0,1] \).
+restrict :: forall v r p. ( Torsor v p, Ring.Field r, Module r v ) => Bezier p -> ( r , r ) -> Bezier p
+restrict bez ( a, b ) = fst $ ( flip ( subdivide @v ) b' ) $ snd $ subdivide @v bez a
+  where
+    b' = ( b Ring.- a ) Ring./ ( Ring.fromInteger 1 Ring.- a )
+  -- TODO: this could be made more efficient.
+  -- See e.g. "https://math.stackexchange.com/questions/4172835/cubic-b%C3%A9zier-spline-multiple-split"
+  -- or the paper "On the numerical condition of Bernstein-Bézier subdivision process".
+{-# INLINEABLE restrict #-}
+
 -- | Polynomial coefficients of the derivative of the distance to a quadratic Bézier curve.
 ddist :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> p -> [ r ]
 ddist ( Bezier {..} ) c = [ a3, a2, a1, a0 ]
diff --git a/brush-strokes/src/Math/Bezier/Spline.hs b/brush-strokes/src/lib/Math/Bezier/Spline.hs
similarity index 100%
rename from brush-strokes/src/Math/Bezier/Spline.hs
rename to brush-strokes/src/lib/Math/Bezier/Spline.hs
diff --git a/brush-strokes/src/Math/Bezier/Stroke.hs b/brush-strokes/src/lib/Math/Bezier/Stroke.hs
similarity index 82%
rename from brush-strokes/src/Math/Bezier/Stroke.hs
rename to brush-strokes/src/lib/Math/Bezier/Stroke.hs
index 0db8ed5..6113b74 100644
--- a/brush-strokes/src/Math/Bezier/Stroke.hs
+++ b/brush-strokes/src/lib/Math/Bezier/Stroke.hs
@@ -148,6 +148,7 @@ import Math.Orientation
   ( Orientation(..), splineOrientation
   , between
   )
+import qualified Math.Ring as Ring
 import Math.Roots
 
 import Debug.Utils
@@ -600,7 +601,7 @@ outlineFunction rootAlgo ptParams toBrushParams brushFromParams = \ sp0 crv ->
 
     ( newtTrees, ( newtDunno, newtSols ) ) =
       intervalNewtonGS
-        NoPreconditioning --InverseMidJacobian
+        InverseMidJacobian
         1e-7
         curvesI
 
@@ -634,7 +635,7 @@ outlineFunction rootAlgo ptParams toBrushParams brushFromParams = \ sp0 crv ->
     logContents = unlines $ functionDataLines ++ treeLines
 
   in trace (unlines solLines) $
-     logToFile cuspFindingMVar logContents `seq`
+     --logToFile cuspFindingMVar logContents `seq`
      OutlineInfo
        { outlineFn = fwdBwd
        , outlineDefiniteCusps  = map ( cuspCoords curves ) newtSols
@@ -1112,7 +1113,7 @@ solveEnvelopeEquations rootAlgo _t path_t path'_t ( fwdOffset, bwdOffset ) strok
                 ++ "}"
               ]
 
-      in logToFile rootSolvingMVar logContents `seq`
+      in --logToFile rootSolvingMVar logContents `seq`
          ( good, ds, dcdt )
 
     (runSolveMethod, methodName) = case rootAlgo of
@@ -1206,7 +1207,7 @@ brushStrokeData co1 co2 path params brush =
         splines :: Seq ( D k ( I i brushParams ) ( I i ( ℝ 1 ) `arr` I i ( ℝ 2 ) ) )
         !splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns @k @i co2 ) dbrush_params
         dbrushes_t :: Seq ( I i ( ℝ 1 ) -> D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) )
-        !dbrushes_t = force $ fmap ( uncurryD @k . chain @(I i Double) @k dparams_t ) splines
+        !dbrushes_t = force $ fmap ( uncurryD @k . chain @( I i Double ) @k dparams_t ) splines
           -- This is the crucial use of the chain rule.
 
     in fmap ( mkStrokeDatum dpath_t ) dbrushes_t
@@ -1274,6 +1275,115 @@ gaussSeidel
         return ( ( 𝕀 ( ℝ1 x1'_lo ) ( ℝ1 x1'_hi ), 𝕀 ( ℝ1 x2'_lo ) ( ℝ1 x2'_hi ) )
                , sub_x1 && sub_x2 )
 
+{-
+gaussSeidel2 :: Int
+             -> Double
+             -> Double
+             -> ( 𝕀ℝ 2, 𝕀ℝ 2 ) -- ^ columns of \( A \)
+             -> 𝕀ℝ 2           -- ^ \( B \)
+             -> ( 𝕀ℝ 1, 𝕀ℝ 1 ) -- ^ initial box \( X \)
+             -> [ ( ( 𝕀ℝ 1, 𝕀ℝ 1 ), Bool ) ]
+gaussSeidel2 maxIters eps_abs eps_rel
+  ( 𝕀 ( ℝ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 ) )
+  x0
+    = 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
+
+          -- See "Algorithm 2" in
+          --  "Using interval unions to solve linear systems of equations with uncertainties"
+          iter ( 𝕀 ( ℝ1 x1_lo ) ( ℝ1 x1_hi ), 𝕀 ( ℝ1 x2_lo ) ( ℝ1 x2_hi ) ) = do
+            let
+              !x1  = 𝕀 x1_lo x1_hi
+              !x2  = 𝕀 x2_lo x2_hi
+              blah1 = do
+                let s = b1 - a12 * x2
+                let s1 = s `monus` ( a11 * x1 )
+                y1 <-
+                  if not $ containsZero ( s1 - a11 * x1 )
+                  then []
+                  else
+                    if containsZero s1 && containsZero a11
+                    then [ ( x1, False ) ]
+                    else do
+                      x1'' <- s1 `extendedDivide` a11
+                      x1'' `intersect` x1
+
+                let s2 = s `monus` ( a12 * x2 )
+                y2 <-
+                  if not $ containsZero ( s2 - a11 * x1 )
+                  then []
+                  else
+                    if containsZero s2 && containsZero a11
+                    then [ ( x1, False ) ]
+                    else do
+                      x1'' <- s2 `extendedDivide` a12
+                      x1'' `intersect` x1
+                return ( y1 `cart` y2 )
+
+              blah2 = do
+                let s = b2 - a21 * x1
+                let s1 = s `monus` ( a21 * x1 )
+                y1 <-
+                  if not $ containsZero ( s1 - a22 * x2 )
+                  then []
+                  else
+                    if containsZero s1 && containsZero a22
+                    then [ ( x1, False ) ]
+                    else do
+                      x1'' <- s1 `extendedDivide` a21
+                      x1'' `intersect` x1
+
+                let s2 = s `monus` ( a12 * x2 )
+                y2 <-
+                  if not $ containsZero ( s2 - a22 * x2 )
+                  then []
+                  else
+                    if containsZero s2 && containsZero a22
+                    then [ ( x1, False ) ]
+                    else do
+                      x1'' <- s2 `extendedDivide` a22
+                      x1'' `intersect` x1
+                return ( y1 `cart` y2 )
+
+            blah1 ++ blah2
+
+          go :: Int -> ( 𝕀ℝ 1, 𝕀ℝ 1 ) -> [ ( ( 𝕀ℝ 1, 𝕀ℝ 1 ), Bool ) ]
+          go !i x
+            = do { nxt@( x', sub ) <- iter x
+                 ; let dw_abs = maxWidth x - maxWidth x'
+                       dw_rel = 1 - ( maxWidth x' / maxWidth x )
+                 ; if sub
+                   || i >= maxIters
+                   || dw_abs < eps_abs
+                   || dw_rel < eps_rel
+                   then return nxt
+                   else go ( i + 1 ) x'
+                 }
+
+      in go 1 x0
+  where
+    maxWidth :: ( 𝕀ℝ 1, 𝕀ℝ 1 ) -> Double
+    maxWidth ( 𝕀 ( ℝ1 x1_lo ) ( ℝ1 x1_hi ), 𝕀 ( ℝ1 x2_lo ) ( ℝ1 x2_hi ) )
+      = max ( x1_hi - x1_lo ) ( x2_hi - x2_lo )
+    containsZero :: 𝕀 Double -> Bool
+    containsZero ( 𝕀 lo hi ) = lo <= 0 && hi >= 0
+    monus :: 𝕀 Double -> 𝕀 Double -> 𝕀 Double
+    monus ( 𝕀 lo1 hi1 ) ( 𝕀 lo2 hi2 )
+      | hi1 - lo1 >= hi2 - lo2
+      = 𝕀 ( lo1 - lo2 ) ( hi1 - hi2 )
+      | otherwise
+      = 𝕀 ( hi1 - hi2 ) ( lo1 - lo2 )
+    cart :: ( 𝕀 Double, Bool ) -> ( 𝕀 Double, Bool ) -> ( ( 𝕀ℝ 1, 𝕀ℝ 1 ), Bool )
+    cart ( 𝕀 lo1 hi1, sub1 ) ( 𝕀 lo2 hi2, sub2 ) =
+      ( ( 𝕀 ( ℝ1 lo1 ) ( ℝ1 hi1 ), 𝕀 ( ℝ1 lo2 ) ( ℝ1 hi2 ) ), sub1 && sub2 )
+-}
+
 -- | Compute the intersection of two intervals.
 --
 -- Returns whether the first interval is a strict subset of the second interval,
@@ -1342,7 +1452,7 @@ data IntervalNewtonLeaf d
 showIntervalNewtonTree :: Box -> IntervalNewtonTree Box -> Tree String
 showIntervalNewtonTree cand (IntervalNewtonLeaf l) = Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ show l) []
 showIntervalNewtonTree cand (IntervalNewtonStep s ts)
-  = Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ show s) $ map (\ (c,t) -> showIntervalNewtonTree c t) ts
+  = Node (show cand ++ " abc " ++ showArea (boxArea cand) ++ " " ++ show s) $ map (\ (c,t) -> showIntervalNewtonTree c t) ts
 
 boxArea :: Box -> Double
 boxArea ( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ), _, 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
@@ -1383,96 +1493,152 @@ intervalNewtonGSFrom
   -> ( [ ( Box, IntervalNewtonTree Box ) ], ( [ Box ], [ Box ] ) )
 intervalNewtonGSFrom precondMethod minWidth eqs initBox =
   runWriter $
-    map ( initBox , ) <$> go initBox
+    map ( initBox , ) <$> evalStrokeDataAndGo initBox
 
   where
 
-    recur f cands = do
-      rest <- traverse ( \ c -> do { trees <- go c; return [ (c, tree) | tree <- trees ] } ) cands
+    recur :: ( cand -> Writer ( [ Box ], [ Box ] ) [ IntervalNewtonTree Box ] )
+          -> ( [ ( cand, IntervalNewtonTree Box ) ] -> IntervalNewtonTree Box )
+          -> [ cand ]
+          -> Writer ( [Box], [Box] ) [ IntervalNewtonTree Box ]
+    recur r f cands = do
+      rest <- traverse ( \ c -> do { trees <- r c; return [ (c, tree) | tree <- trees ] } ) cands
       return [ f $ concat rest ]
 
-    go :: Box -- box to work on
+    evalStrokeDataAndGo :: Box -> Writer ( [Box], [Box] ) [ IntervalNewtonTree Box ]
+    evalStrokeDataAndGo box@( t, i, s ) = go ( box, ( eqs t `Seq.index` i ) s )
+
+    go :: ( Box, StrokeDatum 3 𝕀  ) -- box to work on
        -> Writer ( [ Box ], [ Box ] )
                  [ IntervalNewtonTree Box ]
-    go cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) )
-            , i
-            , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
-            )
+    go ( cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) )
+              , i
+              , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
+              )
+       , sd@( StrokeDatum { ee = D22 _ ( T ee_t ) (T ee_s ) _ _ _
+                          , 𝛿E𝛿sdcdt = D12 _ ( T ( T f_t ) ) ( T ( T f_s ) ) })
+       )
       -- Box is small: stop processing it.
       | t_hi - t_lo < minWidth && s_hi - s_lo < minWidth
       = do let res = TooSmall cand
            tell ( [ cand ], [] )
            return [ IntervalNewtonLeaf res ]
-
-      | StrokeDatum { ee = D22 ee ( T _ee_t ) ( T _ee_s ) _ _ _
-                    , 𝛿E𝛿sdcdt = D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) }
-        <- ( eqs t `Seq.index` i ) s
-
+      | -- Check the envelope equation interval contains zero.
+        ee_potential_zero sd
+        -- Check the 𝛿E𝛿sdcdt interval contains zero.
+      , 𝛿E𝛿sdcdt_potential_zero sd
       , StrokeDatum { ee = D22 _ee_mid _ _ _ _ _
                     , 𝛿E𝛿sdcdt = D12 ( T f_mid ) ( T ( T _f_t_mid ) ) ( T ( T _f_s_mid ) ) }
         <- ( eqs i_t_mid `Seq.index` i ) i_s_mid
-      , let ee_potential_zero = inf ee <= ℝ1 0 && sup ee >= ℝ1 0
-            𝛿E𝛿sdcdt_potential_zero = cmpℝ2 (<=) ( inf f ) ( ℝ2 0 0 ) && cmpℝ2 (>=) ( sup f ) ( ℝ2 0 0 )
-      = if | -- Check the envelope equation interval contains zero.
-             ee_potential_zero
-             -- Check the 𝛿E𝛿sdcdt interval contains zero.
-           , 𝛿E𝛿sdcdt_potential_zero
-           -> 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 )
-                  --(a, b)
-                  --  | 𝕀 (ℝ1 ee_lo) (ℝ1 ee_hi) <- ee_mid
-                  --  , 𝕀 (ℝ1 ee_t_lo) (ℝ1 ee_t_hi) <- ee_t
-                  --  , 𝕀 (ℝ1 ee_s_lo) (ℝ1 ee_s_hi) <- ee_s
-                  --  , 𝕀 (ℝ2 fx_lo fy_lo) (ℝ2 fx_hi fy_hi) <- f_mid
-                  --  , 𝕀 (ℝ2 f_tx_lo f_ty_lo) (ℝ2 f_tx_hi f_ty_hi) <- f_t
-                  --  , 𝕀 (ℝ2 f_sx_lo f_sy_lo) (ℝ2 f_sx_hi f_sy_hi) <- f_s
-                  --  = ( ( 𝕀 (ℝ2 f_tx_lo ee_t_lo) (ℝ2 f_tx_hi ee_t_hi)
-                  --      , 𝕀 (ℝ2 f_sx_lo ee_s_lo) (ℝ2 f_sx_hi ee_s_hi)
-                  --      )
-                  --    , neg $ 𝕀 (ℝ2 fx_lo ee_lo) (ℝ2 fx_hi ee_hi)
-                  --    )
+      = 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 )
+            --(a, b)
+            --  | 𝕀 (ℝ1 ee_lo) (ℝ1 ee_hi) <- ee_mid
+            --  , 𝕀 (ℝ1 ee_t_lo) (ℝ1 ee_t_hi) <- ee_t
+            --  , 𝕀 (ℝ1 ee_s_lo) (ℝ1 ee_s_hi) <- ee_s
+            --  , 𝕀 (ℝ2 fx_lo fy_lo) (ℝ2 fx_hi fy_hi) <- f_mid
+            --  , 𝕀 (ℝ2 f_tx_lo f_ty_lo) (ℝ2 f_tx_hi f_ty_hi) <- f_t
+            --  , 𝕀 (ℝ2 f_sx_lo f_sy_lo) (ℝ2 f_sx_hi f_sy_hi) <- f_s
+            --  = ( ( 𝕀 (ℝ2 f_tx_lo ee_t_lo) (ℝ2 f_tx_hi ee_t_hi)
+            --      , 𝕀 (ℝ2 f_sx_lo ee_s_lo) (ℝ2 f_sx_hi ee_s_hi)
+            --      )
+            --    , neg $ 𝕀 (ℝ2 fx_lo ee_lo) (ℝ2 fx_hi ee_hi)
+            --    )
 
-                  !gsGuesses = gaussSeidel a b
-                                 ( coerce ( (-) @( 𝕀 Double ) ) t i_t_mid
-                                 , coerce ( (-) @( 𝕀 Double ) ) s i_s_mid )
-              in if all ( smaller . fst ) gsGuesses
-                 then
-                   -- 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 ( mkGuess . fst ) ) ( map ( mkGuess . fst ) )
-                                     $ partition snd gsGuesses
-                   in do tell ([], done)
-                         case todo of
-                           [] -> return [ IntervalNewtonLeaf $ NoSolution "GaussSeidel" cand ]
-                           _ -> recur (IntervalNewtonStep (IntervalNewtonContraction done)) todo
-                 else
-                    -- Gauss–Seidel failed to shrink the boxes.
-                    -- Bisect along the widest dimension instead.
-                   let (bisGuesses, whatBis)
-                         | t_hi - t_lo > s_hi - s_lo
-                         = ( [ ( 𝕀 ( ℝ1 t_lo  ) ( ℝ1 t_mid ), i, s )
-                             , ( 𝕀 ( ℝ1 t_mid ) ( ℝ1 t_hi  ), i, s ) ]
-                           , ("t", t_mid) )
-                         | otherwise
-                         = ( [ ( t, i, 𝕀 ( ℝ1 s_lo  ) ( ℝ1 s_mid ) )
-                             , ( t, i, 𝕀 ( ℝ1 s_mid ) ( ℝ1 s_hi  ) ) ]
-                           , ("s", s_mid) )
-                   in recur (IntervalNewtonStep (IntervalNewtonBisection whatBis)) bisGuesses
+            !gsGuesses = gaussSeidel a b
+                           ( coerce ( (-) @( 𝕀 Double ) ) t i_t_mid
+                           , coerce ( (-) @( 𝕀 Double ) ) s i_s_mid )
+        in if any ( smaller . fst ) gsGuesses
+           then
+             -- 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 ( mkGuess . fst ) ) ( map ( mkGuess . fst ) )
+                               $ partition snd gsGuesses
+             in do tell ([], done)
+                   case todo of
+                     [] -> return [ IntervalNewtonLeaf $ NoSolution "GaussSeidel" cand ]
+                     _ -> recur evalStrokeDataAndGo ( IntervalNewtonStep ( IntervalNewtonContraction done ) )
+                            todo
+           else
+              -- Gauss–Seidel failed to shrink the boxes, so bisect instead.
+              -- We have to pick along which dimension to bisect:
+              --  - if bisecting along a particular dimension discards one of
+              --    the boxes, do that;
+              --  - otherwise, bisect along the dimension j that maximises
+              --    width(x_j) * || J_j ||.
+             let l_t = 𝕀 ( ℝ1 t_lo  ) ( ℝ1 t_mid )
+                 r_t = 𝕀 ( ℝ1 t_mid ) ( ℝ1 t_hi )
+                 d_s = 𝕀 ( ℝ1 s_lo  ) ( ℝ1 s_mid )
+                 u_s = 𝕀 ( ℝ1 s_mid ) ( ℝ1 s_hi  )
+                 l = ( l_t, i, s )
+                 r = ( r_t, i, s )
+                 d = ( t, i, d_s )
+                 u = ( t, i, u_s )
+                 l_dat = ( eqs l_t `Seq.index` i ) s
+                 r_dat = ( eqs r_t `Seq.index` i ) s
+                 d_dat = ( eqs t   `Seq.index` i ) d_s
+                 u_dat = ( eqs t   `Seq.index` i ) u_s
+                 l_skip =
+                     not ( ee_potential_zero       l_dat )
+                  || not ( 𝛿E𝛿sdcdt_potential_zero l_dat )
+                 r_skip =
+                      not ( ee_potential_zero       r_dat )
+                   || not ( 𝛿E𝛿sdcdt_potential_zero r_dat )
+                 d_skip =
+                      not ( ee_potential_zero       d_dat )
+                   || not ( 𝛿E𝛿sdcdt_potential_zero d_dat )
+                 u_skip =
+                      not ( ee_potential_zero       u_dat )
+                   || not ( 𝛿E𝛿sdcdt_potential_zero u_dat )
+                 tJWidth = ( t_hi - t_lo ) * normVI3 ee_t f_t
+                 sJWidth = ( s_hi - s_lo ) * normVI3 ee_s f_s
+                 ( bisGuesses, whatBis )
+                   | l_skip && r_skip
+                   = ( [], ( "lr", t_mid ) )
+                   | d_skip && u_skip
+                   = ( [], ( "du", s_mid ) )
+                   | l_skip
+                   = ( [ ( r, r_dat ) ], ( "r", t_mid ) )
+                   | r_skip
+                   = ( [ ( l, l_dat ) ], ( "l", t_mid ) )
+                   | d_skip
+                   = ( [ ( u, u_dat ) ], ( "u", s_mid ) )
+                   | u_skip
+                   = ( [ ( d, d_dat ) ], ( "d", s_mid ) )
+                   | tJWidth >= sJWidth
+                 --  t_hi - t_lo >= s_hi - s_lo
+                   = ( [ ( l, l_dat ), ( r, r_dat ) ], ( "t", t_mid ) )
+                   | otherwise
+                   = ( [ ( d, d_dat ), ( u, u_dat ) ], ( "s", s_mid ) )
+             in recur go ( IntervalNewtonStep ( IntervalNewtonBisection whatBis ) . map (first fst) )
+                  bisGuesses
 
            -- Box doesn't contain a solution: discard it.
            | otherwise
-           -> return [ IntervalNewtonLeaf $ NoSolution (if ee_potential_zero then "dc/dt" else "ee") cand ]
+           = return [ IntervalNewtonLeaf $ NoSolution ( if ee_potential_zero sd then "dc/dt" else "ee" ) cand ]
       where
         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 :: 𝕀ℝ 1 -> 𝕀ℝ 2 -> Double
+        normVI3 ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 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
         t_mid = 0.5 * ( t_lo + t_hi )
         s_mid = 0.5 * ( s_lo + s_hi )
         i_t_mid = singleton ( ℝ1 t_mid )
@@ -1490,6 +1656,16 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox =
                 !( 𝕀 y'_lo y'_hi ) = negate $ 𝕀 y_lo y_hi
             in 𝕀 ( ℝ2 x'_lo y'_lo ) ( ℝ2 x'_hi y'_hi )
 
+ee_potential_zero :: StrokeDatum 3 𝕀 -> Bool
+ee_potential_zero dat =
+     inf ( _D22_v $ ee dat ) <= ℝ1 0
+  && sup ( _D22_v $ ee dat ) >= ℝ1 0
+𝛿E𝛿sdcdt_potential_zero :: StrokeDatum 3 𝕀 -> Bool
+𝛿E𝛿sdcdt_potential_zero dat =
+     cmpℝ2 (<=) ( inf $ unT $ _D12_v $ 𝛿E𝛿sdcdt dat ) ( ℝ2 0 0 )
+  && cmpℝ2 (>=) ( sup $ unT $ _D12_v $ 𝛿E𝛿sdcdt dat ) ( ℝ2 0 0 )
+
+
 zero, one :: Double
 zero = 1e-6
 one = 1 - zero
diff --git a/brush-strokes/src/Math/Bezier/Stroke/EnvelopeEquation.hs b/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs
similarity index 89%
rename from brush-strokes/src/Math/Bezier/Stroke/EnvelopeEquation.hs
rename to brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs
index 9c53198..d5a8adb 100644
--- a/brush-strokes/src/Math/Bezier/Stroke/EnvelopeEquation.hs
+++ b/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs
@@ -332,3 +332,31 @@ evaluateQuadratic bez t =
       maxs = fmap (Quadratic.bezier @( T Double ) sup_bez)
            $ inf t :| ( sup t : filter ( inside t ) ( Quadratic.extrema sup_bez ) )
   in 𝕀 ( minimum mins ) ( maximum maxs )
+
+{-
+
+evaluateCubic :: Cubic.Bezier ( 𝕀 Double ) -> 𝕀 Double -> 𝕀 Double
+evaluateCubic bez t =
+  -- assert (inf t >= 0 && sup t <= 1) "evaluateCubic: t ⊊ [0,1]" $ -- Requires t ⊆ [0,1]
+  let inf_bez = Cubic.restrict @( T Double ) ( fmap inf bez ) ( inf t, sup t )
+      sup_bez = Cubic.restrict @( T Double ) ( fmap sup bez ) ( inf t, sup t )
+      mins = fmap (Cubic.bezier @( T Double ) inf_bez)
+           $ 0 :| ( 1 : Cubic.extrema inf_bez )
+      maxs = fmap (Cubic.bezier @( T Double ) sup_bez)
+           $ 0 :| ( 1 : Cubic.extrema sup_bez )
+  in 𝕀 ( minimum mins ) ( maximum maxs )
+
+-- | Evaluate a quadratic Bézier curve, when both the coefficients and the
+-- parameter are intervals.
+evaluateQuadratic :: Quadratic.Bezier ( 𝕀 Double ) -> 𝕀 Double -> 𝕀 Double
+evaluateQuadratic bez t =
+  -- assert (inf t >= 0 && sup t <= 1) "evaluateCubic: t ⊊ [0,1]" $ -- Requires t ⊆ [0,1]
+  let inf_bez = Quadratic.restrict @( T Double ) ( fmap inf bez ) ( inf t, sup t )
+      sup_bez = Quadratic.restrict @( T Double ) ( fmap sup bez ) ( inf t, sup t )
+      mins = fmap (Quadratic.bezier @( T Double ) inf_bez)
+           $ 0 :| ( 1 : Quadratic.extrema inf_bez )
+      maxs = fmap (Quadratic.bezier @( T Double ) sup_bez)
+           $ 0 :| ( 1 : Quadratic.extrema sup_bez )
+  in 𝕀 ( minimum mins ) ( maximum maxs )
+
+-}
\ No newline at end of file
diff --git a/brush-strokes/src/Math/Differentiable.hs b/brush-strokes/src/lib/Math/Differentiable.hs
similarity index 100%
rename from brush-strokes/src/Math/Differentiable.hs
rename to brush-strokes/src/lib/Math/Differentiable.hs
diff --git a/brush-strokes/src/Math/Epsilon.hs b/brush-strokes/src/lib/Math/Epsilon.hs
similarity index 100%
rename from brush-strokes/src/Math/Epsilon.hs
rename to brush-strokes/src/lib/Math/Epsilon.hs
diff --git a/brush-strokes/src/Math/Interval.hs b/brush-strokes/src/lib/Math/Interval.hs
similarity index 100%
rename from brush-strokes/src/Math/Interval.hs
rename to brush-strokes/src/lib/Math/Interval.hs
diff --git a/brush-strokes/src/Math/Interval/FMA.hs b/brush-strokes/src/lib/Math/Interval/FMA.hs
similarity index 100%
rename from brush-strokes/src/Math/Interval/FMA.hs
rename to brush-strokes/src/lib/Math/Interval/FMA.hs
diff --git a/brush-strokes/src/Math/Interval/Internal.hs b/brush-strokes/src/lib/Math/Interval/Internal.hs
similarity index 94%
rename from brush-strokes/src/Math/Interval/Internal.hs
rename to brush-strokes/src/lib/Math/Interval/Internal.hs
index b170745..d58e8e6 100644
--- a/brush-strokes/src/Math/Interval/Internal.hs
+++ b/brush-strokes/src/lib/Math/Interval/Internal.hs
@@ -120,19 +120,7 @@ instance Prelude.Fractional ( 𝕀 Double ) where
     | otherwise
     = error "BAD interval recip; should use extendedRecip instead"
 -- #endif
-  𝕀 x_lo x_hi / 𝕀 y_lo y_hi
--- #ifdef ASSERTS
-    | y_lo == 0
-    = 𝕀 ( fst $ divI x_lo y_hi ) ( 1 / 0 )
-    | y_hi == 0
-    = 𝕀 ( -1 / 0 ) ( snd $ divI x_hi y_lo )
-    | y_lo > 0 || y_hi < 0
--- #endif
-    = 𝕀 ( fst $ divI x_lo y_hi ) ( snd $ divI x_hi y_lo )
--- #ifdef ASSERTS
-    | otherwise
-    = error "BAD interval division; should use extendedRecip instead"
--- #endif
+  p / q = p * recip q
 
 instance Floating ( 𝕀 Double ) where
   sqrt = withHW Prelude.sqrt
diff --git a/brush-strokes/src/Math/Linear.hs b/brush-strokes/src/lib/Math/Linear.hs
similarity index 100%
rename from brush-strokes/src/Math/Linear.hs
rename to brush-strokes/src/lib/Math/Linear.hs
diff --git a/brush-strokes/src/Math/Linear/Internal.hs b/brush-strokes/src/lib/Math/Linear/Internal.hs
similarity index 100%
rename from brush-strokes/src/Math/Linear/Internal.hs
rename to brush-strokes/src/lib/Math/Linear/Internal.hs
diff --git a/brush-strokes/src/Math/Linear/Solve.hs b/brush-strokes/src/lib/Math/Linear/Solve.hs
similarity index 100%
rename from brush-strokes/src/Math/Linear/Solve.hs
rename to brush-strokes/src/lib/Math/Linear/Solve.hs
diff --git a/brush-strokes/src/Math/Module.hs b/brush-strokes/src/lib/Math/Module.hs
similarity index 100%
rename from brush-strokes/src/Math/Module.hs
rename to brush-strokes/src/lib/Math/Module.hs
diff --git a/brush-strokes/src/Math/Module/Internal.hs b/brush-strokes/src/lib/Math/Module/Internal.hs
similarity index 100%
rename from brush-strokes/src/Math/Module/Internal.hs
rename to brush-strokes/src/lib/Math/Module/Internal.hs
diff --git a/brush-strokes/src/Math/Monomial.hs b/brush-strokes/src/lib/Math/Monomial.hs
similarity index 100%
rename from brush-strokes/src/Math/Monomial.hs
rename to brush-strokes/src/lib/Math/Monomial.hs
diff --git a/brush-strokes/src/Math/Orientation.hs b/brush-strokes/src/lib/Math/Orientation.hs
similarity index 100%
rename from brush-strokes/src/Math/Orientation.hs
rename to brush-strokes/src/lib/Math/Orientation.hs
diff --git a/brush-strokes/src/Math/Ring.hs b/brush-strokes/src/lib/Math/Ring.hs
similarity index 100%
rename from brush-strokes/src/Math/Ring.hs
rename to brush-strokes/src/lib/Math/Ring.hs
diff --git a/brush-strokes/src/Math/Roots.hs b/brush-strokes/src/lib/Math/Roots.hs
similarity index 100%
rename from brush-strokes/src/Math/Roots.hs
rename to brush-strokes/src/lib/Math/Roots.hs
diff --git a/brush-strokes/src/TH/Utils.hs b/brush-strokes/src/lib/TH/Utils.hs
similarity index 100%
rename from brush-strokes/src/TH/Utils.hs
rename to brush-strokes/src/lib/TH/Utils.hs
diff --git a/src/convert/Main.hs b/brush-strokes/src/metafont/Main.hs
similarity index 100%
rename from src/convert/Main.hs
rename to brush-strokes/src/metafont/Main.hs
diff --git a/src/convert/MetaBrush/MetaFont/Convert.hs b/brush-strokes/src/metafont/MetaBrush/MetaFont/Convert.hs
similarity index 100%
rename from src/convert/MetaBrush/MetaFont/Convert.hs
rename to brush-strokes/src/metafont/MetaBrush/MetaFont/Convert.hs
diff --git a/src/app/MetaBrush/Application.hs b/src/app/MetaBrush/Application.hs
index 5f62bed..253626f 100644
--- a/src/app/MetaBrush/Application.hs
+++ b/src/app/MetaBrush/Application.hs
@@ -158,7 +158,7 @@ runApplication application = do
             { strokeName    = "Stroke 1"
             , strokeVisible = True
             , strokeUnique  = strokeUnique
-            , strokeBrush   = Just Asset.Brushes.tearDrop
+            , strokeBrush   = Just Asset.Brushes.ellipse --tearDrop
             , strokeSpline  =
 --              Spline
 --                { splineStart  = mkPoint ( ℝ2 -20 -20 ) 5
@@ -169,7 +169,7 @@ runApplication application = do
               Spline
                 { splineStart  = mkPoint ( ℝ2 0 0 ) 10 25 0
                 , splineCurves = OpenCurves $ Seq.fromList
-                  [ LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 100 0 ) 15 40 pi ), curveData = invalidateCache undefined }
+                  [ LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 100 0 ) 15 40 (0.1 * pi) ), curveData = invalidateCache undefined }
                   --, LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 -10  10 )  8 5 ( pi / 4 ) ), curveData = invalidateCache undefined }
                   --, LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 -10 -20 ) 10 7 ( pi / 2 ) ), curveData = invalidateCache undefined }
                   ]
@@ -179,12 +179,12 @@ runApplication application = do
         )
       ]
         where
-          --mkPoint :: ℝ 2 -> Double -> Double -> Double -> PointData ( Record Asset.Brushes.EllipseBrushFields )
-          --mkPoint pt a b phi = PointData pt Normal ( MkR $ ℝ3 a b phi )
+          mkPoint :: ℝ 2 -> Double -> Double -> Double -> PointData ( Record Asset.Brushes.EllipseBrushFields )
+          mkPoint pt a b phi = PointData pt Normal ( MkR $ ℝ3 a b phi )
           --mkPoint :: ℝ 2 -> Double -> PointData ( Record Asset.Brushes.CircleBrushFields )
           --mkPoint pt r = PointData pt Normal ( MkR $ ℝ1 r )
-          mkPoint ::  ℝ 2 -> Double -> Double -> Double -> PointData ( Record Asset.Brushes.TearDropBrushFields )
-          mkPoint pt w h phi = PointData pt Normal ( MkR $ ℝ3 w h phi )
+          --mkPoint :: ℝ 2 -> Double -> Double -> Double -> PointData ( Record Asset.Brushes.TearDropBrushFields )
+          --mkPoint pt w h phi = PointData pt Normal ( MkR $ ℝ3 w h phi )
 
   recomputeStrokesTVar <- STM.newTVarIO @Bool                                    False
   documentRenderTVar   <- STM.newTVarIO @( ( Int32, Int32 ) -> Cairo.Render () ) ( const $ pure () )
diff --git a/src/metabrushes/MetaBrush/Records.hs b/src/metabrushes/MetaBrush/Records.hs
index 0d86207..d154ac2 100644
--- a/src/metabrushes/MetaBrush/Records.hs
+++ b/src/metabrushes/MetaBrush/Records.hs
@@ -139,12 +139,8 @@ instance ( Torsor ( T ( 𝕀ℝ ( Length ks ) ) ) ( 𝕀ℝ ( Length ks ) )
             T ( 𝕀 ( MkR c_lo ) ( MkR c_hi ) )
 
 type instance RepDim ( Record ks ) = Length ks
-instance Representable r ( ℝ ( Length ks ) )
-      => Representable r ( Record ks ) where
-  tabulate f = Record $ tabulate f
-  {-# INLINE tabulate #-}
-  index f (Record r) = index f r
-  {-# INLINE index #-}
+deriving newtype instance Representable r ( ℝ ( Length ks ) )
+                       => Representable r ( Record ks )
 
 type instance D k ( Record ks ) = D k ( ℝ ( Length ks ) )
 deriving newtype instance HasChainRule Double 2 ( ℝ ( Length ks ) )