compute brush strokes using the envelope equation

MetaBrush.cabal

@@ -36,6 +36,8 @@ common common
            >= 4.17    && < 4.18
       , acts
           ^>= 0.3.1.0
+      , code-page
+          ^>= 0.2.1
       , containers
            >= 0.6.0.1 && < 0.7
       , deepseq
@@ -112,6 +114,7 @@ common common
 common extras
 
     build-depends:
+
         directory
            >= 1.3.4.0    && < 1.4
       , filepath
@@ -180,6 +183,7 @@ library splines
       , Math.Module
       , Math.Orientation
       , Math.Roots
+      , Debug.Utils
 
     build-depends:
         bifunctors

app/Main.hs

@@ -1,6 +1,4 @@
-{-# LANGUAGE BlockArguments    #-}
 {-# LANGUAGE OverloadedStrings #-}
-{-# LANGUAGE TypeApplications  #-}
 
 module Main
   ( main )
@@ -14,6 +12,10 @@ import System.Exit
 import GHC.Conc
   ( getNumProcessors, setNumCapabilities )
 
+-- code-page
+import System.IO.CodePage
+  ( withCP65001 )
+
 -- gi-gio
 import qualified GI.Gio as GIO
 
@@ -30,7 +32,7 @@ import MetaBrush.Application
 --------------------------------------------------------------------------------
 
 main :: IO ()
-main = do
+main = withCP65001 do
 
   procs <- getNumProcessors
   let

src/app/MetaBrush/Application.hs

@@ -84,7 +84,6 @@ import Math.Linear
 import MetaBrush.Action
   ( ActionOrigin(..) )
 import qualified MetaBrush.Asset.Brushes as Asset.Brushes
-  ( EllipseBrushFields, ellipse )
 import MetaBrush.Asset.Colours
   ( getColours )
 import MetaBrush.Asset.Logo
@@ -161,6 +160,12 @@ runApplication application = do
             , strokeUnique  = strokeUnique
             , strokeBrush   = Just Asset.Brushes.ellipse
             , strokeSpline  =
+--              Spline
+--                { splineStart  = mkPoint ( ℝ2 -20 -20 ) 5
+--                , splineCurves = OpenCurves $ Seq.fromList
+--                  [ LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 20 20 ) 5 ), curveData = invalidateCache undefined }
+--                  ]
+--                }
               Spline
                 { splineStart  = mkPoint ( ℝ2 10 -20 ) 2 1 0
                 , splineCurves = OpenCurves $ Seq.fromList
@@ -176,6 +181,8 @@ 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 -> PointData ( Record Asset.Brushes.CircleBrushFields )
+--          mkPoint pt r = PointData pt Normal ( MkR $ ℝ1 r )
 
   recomputeStrokesTVar <- STM.newTVarIO @Bool                                    False
   documentRenderTVar   <- STM.newTVarIO @( ( Int32, Int32 ) -> Cairo.Render () ) ( const $ pure () )
@@ -193,11 +200,11 @@ runApplication application = do
   maxHistorySizeTVar   <- STM.newTVarIO @Int                                     1000
   fitParametersTVar    <- STM.newTVarIO @FitParameters
                             ( FitParameters
-                              { maxSubdiv  = 6
+                              { maxSubdiv  = 2 --3 -- 6
-                              , nbSegments = 12
+                              , nbSegments = 3 --6 -- 12
-                              , dist_tol   = 5e-3
+                              , dist_tol   = 0.1 -- 5e-3
-                              , t_tol      = 1e-4
+                              , t_tol      = 0.1 -- 1e-4
-                              , maxIters   = 100
+                              , maxIters   = 2   -- 100
                               }
                             )
 

src/app/MetaBrush/Render/Document.hs

@@ -292,17 +292,18 @@ strokeRenderData fitParams
                  { inject
                  , project = toUsedParams :: Record pointFields -> Record usedFields }
                  -> do
-                    let embedUsedParams r = inject r brush_defaults
+                    let embedUsedParams = inject brush_defaults
 
                   -- Compute the outline using the brush function.
                     ( outline, fitPts ) <-
-                      computeStrokeOutline @( T ( Record usedFields) ) @clo
+                      computeStrokeOutline @clo fitParams
-                        fitParams ( toUsedParams . brushParams ) ( fun brushFn . embedUsedParams ) spline
+                        ( toUsedParams . brushParams ) embedUsedParams brushFn
+                        spline
                     pure $
                       StrokeWithOutlineRenderData
                         { strokeDataSpline    = spline
                         , strokeOutlineData   = ( outline, fitPts )
-                        , strokeBrushFunction = fun brushFn . embedUsedParams . toUsedParams
+                        , strokeBrushFunction = fun brushFn . fun embedUsedParams . toUsedParams
                         }
         _ -> pure $
                 StrokeRenderData
@@ -604,7 +605,7 @@ drawFitPoint _ zoom ( FitPoint { fitPoint = ℝ2 x y } ) = do
   lift do
     Cairo.save
     Cairo.translate x y
-    Cairo.arc 0 0 ( 2 / zoom ) 0 ( 2 * pi )
+    Cairo.arc 0 0 ( 6 / zoom ) 0 ( 2 * pi ) -- 2 / zoom
     Cairo.setSourceRGBA r g b 0.8
     Cairo.fill
     Cairo.restore
@@ -621,7 +622,7 @@ drawFitPoint _ zoom ( FitTangent { fitPoint = ℝ2 x y, fitTangent = V2 tx ty }
     Cairo.translate x y
     Cairo.moveTo 0 0
     Cairo.lineTo ( 0.05 * tx ) ( 0.05 * ty )
-    Cairo.setLineWidth ( 1 / zoom )
+    Cairo.setLineWidth ( 4 / zoom ) -- 1 / zoom
     Cairo.setSourceRGBA r g b 0.8
     Cairo.stroke
     Cairo.arc 0 0 ( 2 / zoom ) 0 ( 2 * pi )

src/metabrushes/MetaBrush/Asset/Brushes.hs

@@ -3,10 +3,6 @@
 
 module MetaBrush.Asset.Brushes where
 
--- base
-import Data.Coerce
-  ( coerce )
-
 -- containers
 import qualified Data.Sequence as Seq
   ( fromList )
@@ -24,9 +20,9 @@ import qualified Data.HashMap.Strict as HashMap
 -- MetaBrush
 import Math.Bezier.Spline
 import Math.Linear
-  ( ℝ(..), T(..) )
+  ( ℝ(..), Fin(..) )
 import Math.Linear.Dual
-  ( D, type (~>)(..), Var(var), konst )
+  ( D, type (~>)(..), var, konst )
 import Math.Module
   ( Module((^+^), (*^)) )
 import MetaBrush.Brush
@@ -86,17 +82,16 @@ circleBrush :: Record CircleBrushFields ~> Spline 'Closed () ( ℝ 2 )
 circleBrush =
   D \ params ->
     let r :: D ( Record CircleBrushFields ) Double
-        r   = runD ( var @1 ) params
+        r = runD ( var ( Fin 1## ) ) params
         mkPt :: Double -> Double -> D ( Record CircleBrushFields ) ( ℝ 2 )
         mkPt ( kon -> x ) ( kon -> y )
-          = fmap coerce
+          =   ( x * r ) *^ e_x
-          $   ( x * r ) *^ e_x
           ^+^ ( y * r ) *^ e_y
     in circleSpline @( Record CircleBrushFields ) mkPt
   where
-    e_x, e_y :: D ( Record CircleBrushFields ) ( T ( ℝ 2 ) )
+    e_x, e_y :: D ( Record CircleBrushFields ) ( ℝ 2 )
-    e_x = pure $ T $ ℝ2 1 0
+    e_x = pure $ ℝ2 1 0
-    e_y = pure $ T $ ℝ2 0 1
+    e_y = pure $ ℝ2 0 1
 
     kon = konst @( Record CircleBrushFields )
 
@@ -104,25 +99,17 @@ ellipseBrush :: Record EllipseBrushFields ~> Spline 'Closed () ( ℝ 2 )
 ellipseBrush =
   D \ params ->
     let a, b, phi :: D ( Record EllipseBrushFields ) Double
-        a   = runD ( var @1 ) params
+        a   = runD ( var ( Fin 1## ) ) params
-        b   = runD ( var @2 ) params
+        b   = runD ( var ( Fin 2## ) ) params
-        phi = runD ( var @3 ) params
+        phi = runD ( var ( Fin 3## ) ) params
         mkPt :: Double -> Double -> D ( Record EllipseBrushFields ) ( ℝ 2 )
         mkPt ( kon -> x ) ( kon -> y )
-          = fmap coerce
+          =   ( x * a * cos phi - y * b * sin phi ) *^ e_x
-          $   ( x * a * cos phi - y * b * sin phi ) *^ e_x
           ^+^ ( y * b * cos phi + x * a * sin phi ) *^ e_y
     in circleSpline @( Record EllipseBrushFields ) mkPt
   where
-    e_x, e_y :: D ( Record EllipseBrushFields ) ( T ( ℝ 2 ) )
+    e_x, e_y :: D ( Record EllipseBrushFields ) ( ℝ 2 )
-    e_x = pure $ T $ ℝ2 1 0
+    e_x = pure $ ℝ2 1 0
-    e_y = pure $ T $ ℝ2 0 1
+    e_y = pure $ ℝ2 0 1
 
     kon = konst @( Record EllipseBrushFields )
-
---ellipseArc :: ℝ 2 ~> ℝ 2
---ellipseArc = brushStroke ( linear myPath ) ( uncurryD $ fmap bezier3 myBrush )
-
---testing :: Double -> Double -> (# Double, T ( ℝ 2 ) #)
---testing :: Double -> Double -> (# Double, T (ℝ 2) #)
---testing t s = envelopeEquation ellipseArc t s

src/metabrushes/MetaBrush/Brush.hs

@@ -1,6 +1,7 @@
-{-# LANGUAGE AllowAmbiguousTypes  #-}
+{-# LANGUAGE AllowAmbiguousTypes   #-}
-{-# LANGUAGE ScopedTypeVariables  #-}
+{-# LANGUAGE QuantifiedConstraints #-}
-{-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE ScopedTypeVariables   #-}
+{-# LANGUAGE UndecidableInstances  #-}
 
 module MetaBrush.Brush
   ( Brush(..), SomeBrush(..), BrushFunction
@@ -38,12 +39,17 @@ import qualified Data.Text as Text
 
 -- MetaBrush
 import Math.Linear
+  ( ℝ, Representable )
 import Math.Linear.Dual
   ( Diffy )
+import Math.Module
+  ( Interpolatable )
 import Math.Bezier.Spline
   ( SplineType(Closed), SplinePts)
 import MetaBrush.Records
+  ( KnownSymbols, Length, Record, WithParams )
 import MetaBrush.Serialisable
+  ( Serialisable )
 
 --------------------------------------------------------------------------------
 

src/metabrushes/MetaBrush/Document/SubdivideStroke.hs

@@ -49,7 +49,7 @@ import Math.Bezier.Spline
 import Math.Bezier.Stroke
   ( CachedStroke(..), invalidateCache )
 import Math.Module
-  ( lerp, quadrance, closestPointOnSegment )
+  ( Interpolatable, lerp, quadrance, closestPointOnSegment )
 import Math.Linear
   ( Segment(..), ℝ(..), T(..) )
 import MetaBrush.Document
@@ -57,8 +57,6 @@ import MetaBrush.Document
   , PointData(..), DiffPointData(..)
   , coords, _strokeSpline
   )
-import MetaBrush.Records
-  ( Interpolatable )
 
 --------------------------------------------------------------------------------
 

src/metabrushes/MetaBrush/Records.hs

@@ -54,13 +54,6 @@ import Math.Module
 
 --------------------------------------------------------------------------------
 
--- | A convenient constraint synonym for types that support interpolation.
-type Interpolatable :: Type -> Constraint
-class    ( Torsor ( T r ) r, Module Double ( T r ) ) => Interpolatable r
-instance ( Torsor ( T r ) r, Module Double ( T r ) ) => Interpolatable r
-
---------------------------------------------------------------------------------
-
 -- | A function from a given record type, with provided default values
 -- that can be overridden.
 type WithParams :: [ Symbol ] -> Type -> Type
@@ -123,9 +116,11 @@ instance ( Torsor ( T ( ℝ ( Length ks ) ) ) ( ℝ ( Length ks ) )
         => Torsor ( T ( Record ks ) ) ( Record ks ) where
   MkR g --> MkR a = T $ MkR $ unT $ g --> a
 
-type instance D ( Record ks ) = D ( ℝ ( Length ks ) )
+deriving newtype
+  instance Representable ( ℝ ( Length ks ) )
+        => Representable ( Record ks )
 
-deriving newtype instance Var n ( ℝ ( Length ks ) ) => Var n ( Record ks )
+type instance D ( Record ks ) = D ( ℝ ( Length ks ) )
 deriving newtype instance Diffy ( ℝ ( Length ks ) ) => Diffy ( Record ks )
 
 --------------------------------------------------------------------------------
@@ -154,31 +149,35 @@ intersect :: forall r1 r2 l1 l2
             , KnownSymbols r1, KnownSymbols r2
             , l1 ~ Length r1, l2 ~ Length r2
             , Representable ( ℝ l1 ), Representable ( ℝ l2 )
-            , Interpolatable ( Record r1 )
+            , Interpolatable ( Record r1 ), Diffy ( ℝ l2 )
             )
           => Intersection r1 r2
 intersect
   -- Shortcut when the two rows are equal.
   | Just Refl <- eqT @r1 @r2
   , Refl <- ( unsafeCoerce Refl :: r1 :~: Intersect r1 r2 )
-  = Intersection { project = id, inject = const }
+  = Intersection { project = id, inject = \ _ -> linear id }
   | otherwise
   = doIntersection @r1 @r2 \ ( _ :: Proxy# r1r2 ) r1_idxs r2_idxs ->
       let
         project :: Record r1 -> Record r1r2
         project = \ ( MkR r1 ) -> MkR $ projection ( (!) r1_idxs ) r1
 
-        inject :: Record r1r2 -> Record r2 -> Record r2
+        inject :: Record r2 -> Record r1r2 ~> Record r2
-        inject  = \ ( MkR r1r2 ) ( MkR r2 ) -> MkR $ injection ( find eqFin r2_idxs ) r1r2 r2
+        inject = \ ( MkR r2 ) -> linear \ ( MkR r1r2 ) -> MkR $ injection ( find eqFin r2_idxs ) r1r2 r2
       in Intersection { project, inject }
 
 data Intersection r1 r2 where
   Intersection
-    :: forall r1r2 r1 r2
+    :: forall r1r2 r1 r2 l12
-    .  ( KnownSymbols r1r2, Representable ( ℝ ( Length r1r2 ) )
+    .  ( l12 ~ Length r1r2
+       , KnownSymbols r1r2
+       , Representable ( ℝ l12 )
+       , Diffy ( ℝ l12 )
        , Interpolatable ( Record r1r2 ) )
-    => { project :: Record r1   -> Record r1r2
+    => { project :: Record r1 -> Record r1r2
-       , inject  :: Record r1r2 -> Record r2 -> Record r2
+       , inject  :: Record r2 -> Record r1r2 ~> Record r2
+          -- ^ overrides the components of the first record with the second
        } -> Intersection r1 r2
 
 {-# INLINE doIntersection #-}
@@ -190,7 +189,7 @@ doIntersection
      )
   => ( forall r1r2 l12.
         ( r1r2 ~ Intersect r1 r2, l12 ~ Length r1r2
-        , Representable ( ℝ l12 ), Interpolatable ( ℝ l12 )
+        , Representable ( ℝ l12 ), Diffy ( ℝ l12 ), Interpolatable ( ℝ l12 )
         , KnownSymbols r1r2, Representable ( ℝ ( Length r1r2 ) )
         )
       => Proxy# r1r2 -> Vec l12 ( Fin l1 ) -> Vec l12 ( Fin l2 ) -> kont )

src/splines/Debug/Utils.hs

@@ -0,0 +1,26 @@
+module Debug.Utils ( trace ) where
+
+
+-- base
+import System.IO
+  ( BufferMode(..), hSetBuffering, hFlush, hPutStrLn, stdout )
+import System.IO.Unsafe
+  ( unsafePerformIO )
+
+-- deepseq
+import Control.DeepSeq
+  ( force )
+
+-- code-page
+import System.IO.CodePage
+  ( withCP65001 )
+
+--------------------------------------------------------------------------------
+
+trace :: String -> a -> a
+trace ( force -> !str ) expr =
+  unsafePerformIO $ withCP65001 do
+    hSetBuffering stdout ( BlockBuffering $ Just 50 )
+    hPutStrLn stdout str
+    hFlush stdout
+    return expr

src/splines/Math/Bezier/Cubic/Fit.hs

@@ -120,7 +120,7 @@ data FitPoint
 -- including the meaning of \( \texttt{t_tol} \) and \( \texttt{maxIters} \).
 fitSpline
   :: FitParameters
-  -> ( Double -> ( ℝ 2, T ( ℝ 2 ) ) ) -- ^ curve \( t \mapsto ( C(t), C'(t) ) \) to fit
+  -> ( ℝ 1 -> ( ℝ 2, T ( ℝ 2 ) ) ) -- ^ curve \( t \mapsto ( C(t), C'(t) ) \) to fit
   -> ( SplinePts Open, Seq FitPoint )
 fitSpline ( FitParameters {..} ) = go 0
   where
@@ -128,16 +128,16 @@ fitSpline ( FitParameters {..} ) = go 0
     dt = recip ( fromIntegral nbSegments )
     go
       :: Int
-      -> ( Double -> ( ℝ 2, T ( ℝ 2 ) ) )
+      -> ( ℝ 1 -> ( ℝ 2, T ( ℝ 2 ) ) )
       -> ( SplinePts Open, Seq FitPoint )
     go subdiv curveFn =
       let
         p, r :: ℝ 2
         tp, tr :: T ( ℝ 2 )
         qs :: [ ℝ 2 ]
-        (p, tp) = curveFn 0
+        (p, tp) = curveFn $ ℝ1 0
-        (r, tr) = curveFn 1
+        (r, tr) = curveFn $ ℝ1 1
-        qs = [ fst $ curveFn ( dt * fromIntegral j ) | j <- [ 1 .. nbSegments - 1 ] ]
+        qs = [ fst $ curveFn ( ℝ1 $ dt * fromIntegral j ) | j <- [ 1 .. nbSegments - 1 ] ]
       in
         case fitPiece dist_tol t_tol maxIters p tp qs r tr of
           ( bez, Max ( Arg max_sq_error t_split ) )
@@ -150,8 +150,8 @@ fitSpline ( FitParameters {..} ) = go 0
                 c1, c2 :: SplinePts Open
                 ps1, ps2 :: Seq FitPoint
                 ( ( c1, ps1 ), ( c2, ps2 ) )
-                  = ( go ( subdiv + 1 ) ( \ t -> curveFn $ t * t_split_eff )
+                  = ( go ( subdiv + 1 ) ( \ ( ℝ1 t ) -> curveFn $ ℝ1 $ t * t_split_eff )
-                    , go ( subdiv + 1 ) ( \ t -> curveFn $ t_split_eff + t * ( 1 - t_split_eff ) )
+                    , go ( subdiv + 1 ) ( \ ( ℝ1 t ) -> curveFn $ ℝ1 $ t_split_eff + t * ( 1 - t_split_eff ) )
                     ) `Parallel.Strategy.using`
                        ( Parallel.Strategy.parTuple2 Parallel.Strategy.rdeepseq Parallel.Strategy.rdeepseq )
             -> ( c1 <> c2, ps1 <> ps2 )

src/splines/Math/Bezier/Stroke.hs

@@ -1,5 +1,6 @@
-{-# LANGUAGE AllowAmbiguousTypes #-}
+{-# LANGUAGE AllowAmbiguousTypes   #-}
-{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE QuantifiedConstraints #-}
+{-# LANGUAGE ScopedTypeVariables   #-}
 
 module Math.Bezier.Stroke
   ( Offset(..)
@@ -9,10 +10,8 @@ module Math.Bezier.Stroke
 
     -- * Brush stroking
 
-    -- $brushes
   , brushStroke, envelopeEquation
-  , linear, bezier2, bezier3
+  , line, bezier2, bezier3
---  , uncurryD
 
   )
   where
@@ -22,6 +21,8 @@ import Prelude
   hiding ( unzip )
 import Control.Arrow
   ( first, (***) )
+import Control.Applicative
+  ( Applicative(..) )
 import Control.Monad
   ( unless )
 import Control.Monad.ST
@@ -30,6 +31,8 @@ import Data.Bifunctor
   ( Bifunctor(bimap) )
 import Data.Foldable
   ( for_ )
+import Data.Functor.Identity
+  ( Identity(..) )
 import Data.List.NonEmpty
   ( unzip )
 import Data.Maybe
@@ -53,11 +56,15 @@ import Data.Act
 import Data.Sequence
   ( Seq(..) )
 import qualified Data.Sequence as Seq
+  --( empty, reverse, singleton, index )
+import Data.Set
+  ( Set )
+import qualified Data.Set as Set
   ( singleton )
 
 -- deepseq
 import Control.DeepSeq
-  ( NFData(..), NFData1, deepseq )
+  ( NFData(..), NFData1, deepseq, force )
 
 -- generic-lens
 import Data.Generics.Product.Typed
@@ -65,10 +72,6 @@ import Data.Generics.Product.Typed
 import Data.Generics.Internal.VL
   ( set, view )
 
--- groups
-import Data.Group
-  ( Group )
-
 -- parallel
 import qualified Control.Parallel.Strategies as Strats
   ( rdeepseq, parTuple2, using )
@@ -101,8 +104,8 @@ import qualified Math.Bezier.Quadratic as Quadratic
 import Math.Epsilon
   ( epsilon )
 import Math.Module
-  ( Module(..), Inner((^.^))
+  ( Module(..), Inner((^.^)), Interpolatable
-  , lerp, squaredNorm, cross
+  , lerp, cross
   , convexCombination, strictlyParallel
   )
 import Math.Orientation
@@ -110,10 +113,12 @@ import Math.Orientation
   , between
   )
 import Math.Roots
-  ( solveQuadratic )
+  ( solveQuadratic, newtonRaphson )
 import Math.Linear
 import Math.Linear.Dual
 
+--import Debug.Utils
+
 --------------------------------------------------------------------------------
 
 data Offset
@@ -168,42 +173,42 @@ coords = view typed
 
 --------------------------------------------------------------------------------
 
+type OutlineFn = ℝ 1 -> ( ( ℝ 2, T ( ℝ 2 ) ), ( ℝ 2, T ( ℝ 2 ) ) )
+
 computeStrokeOutline ::
-  forall diffParams ( clo :: SplineType ) brushParams crvData ptData s
+  forall ( clo :: SplineType ) usedParams brushParams crvData ptData s
   .  ( KnownSplineType clo
-     , Group  diffParams
+     , Interpolatable usedParams
-     , Module Double diffParams
+     , Diffy usedParams, Diffy brushParams
-     , Torsor diffParams brushParams
      , HasType ( ℝ 2 ) ptData
      , HasType ( CachedStroke s ) crvData
      , NFData ptData, NFData crvData
      -- Debugging.
-     , Show ptData
+     , Show ptData, Show brushParams
      )
   => FitParameters
-  -> ( ptData -> brushParams )
+  -> ( ptData -> usedParams )
-  -> ( brushParams -> SplinePts Closed )
+  -> ( usedParams ~> brushParams )
+  -> ( brushParams ~> SplinePts Closed )
   -> Spline clo crvData ptData
   -> ST s
         ( Either ( SplinePts Closed ) ( SplinePts Closed, SplinePts Closed )
         , Seq FitPoint
         )
-computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart = spt0 } ) = case ssplineType @clo of
+computeStrokeOutline fitParams ptParams toBrushParams brushFn spline@( Spline { splineStart = spt0 } ) = case ssplineType @clo of
   -- Open brush path with at least one segment.
   -- Need to add caps at both ends of the path.
   SOpen
     | firstCurve :<| _ <- openCurves $ splineCurves spline
     , prevCurves :|> lastCurve <- openCurves $ splineCurves spline
-    , ( firstOutlineFwd, firstOutlineBwd ) :<| _ <- outlineFns
+    , firstOutlineFn :<| _ <- outlineFns
-    , _ :|> ( lastOutlineFwd, lastOutlineBwd )   <- outlineFns
+    , _ :|> lastOutlineFn  <- outlineFns
     , let
         endPt :: ptData
         endPt = openCurveEnd lastCurve
         startTgtFwd, startTgtBwd, endTgtFwd, endTgtBwd :: T ( ℝ 2 )
-        startTgtFwd =       snd ( firstOutlineFwd 0 )
+        ( ( _, startTgtFwd), ( _, startTgtBwd ) ) = firstOutlineFn $ ℝ1 0
-        startTgtBwd = -1 *^ snd ( firstOutlineBwd 1 )
+        ( ( _, endTgtFwd  ), ( _, endTgtBwd   ) ) =  lastOutlineFn $ ℝ1 1
-        endTgtFwd   =       snd (  lastOutlineFwd 1 )
-        endTgtBwd   = -1 *^ snd (  lastOutlineBwd 0 )
         startBrush, endBrush :: SplinePts Closed
         startBrush = brushShape spt0
         endBrush   = brushShape endPt
@@ -252,8 +257,8 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart =
   -- Add forward and backward caps at the start.
   SClosed
     | ClosedCurves prevCurves lastCurve <- splineCurves spline
-    , ( firstOutlineFwd, firstOutlineBwd ) :<| _ <- outlineFns
+    , firstOutlineFn :<| _ <- outlineFns
-    , _ :|> ( lastOutlineFwd, lastOutlineBwd )   <- outlineFns
+    , _ :|> lastOutlineFn  <- outlineFns
     , let
         startTgt, endTgt, startTgtFwd, startTgtBwd, endTgtFwd, endTgtBwd :: T ( ℝ 2 )
         startTgt = case prevCurves of
@@ -262,10 +267,8 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart =
         endTgt = case prevCurves of
           Empty          -> endTangent spt0 spt0                      lastCurve
           _ :|> lastPrev -> endTangent spt0 ( openCurveEnd lastPrev ) lastCurve
-        startTgtFwd =       snd ( firstOutlineFwd 0 )
+        ( ( _, startTgtFwd), ( _, startTgtBwd ) ) = firstOutlineFn $ ℝ1 0
-        startTgtBwd = -1 *^ snd ( firstOutlineBwd 1 )
+        ( ( _, endTgtFwd  ), ( _, endTgtBwd   ) ) =  lastOutlineFn $ ℝ1 1
-        endTgtFwd   =       snd (  lastOutlineFwd 1 )
-        endTgtBwd   = -1 *^ snd (  lastOutlineBwd 0 )
         fwdStartCap, bwdStartCap :: SplinePts Open
         TwoSided fwdStartCap bwdStartCap
           = fmap fst . snd . runWriter
@@ -284,26 +287,19 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart =
       )
   where
 
-    outlineFns
+    outlineFns :: Seq OutlineFn
-      :: Seq
-          ( Double -> ( ℝ 2, T ( ℝ 2 ) )
-          , Double -> ( ℝ 2, T ( ℝ 2 ) )
-          )
     outlineFns = go spt0 ( openCurves $ splineCurves ( adjustSplineType @Open spline ) )
       where
         go
           :: ptData
           -> Seq ( Curve Open crvData ptData )
-          -> Seq
+          -> Seq OutlineFn
-              ( Double -> ( ℝ 2, T ( ℝ 2 ) )
-              , Double -> ( ℝ 2, T ( ℝ 2 ) )
-              )
         go _ Empty = Empty
         go p0 ( crv :<| crvs ) =
-          outlineFunctions @diffParams ptParams brushFn p0 crv :<| go ( openCurveEnd crv ) crvs
+          outlineFunction ptParams toBrushParams brushFn p0 crv :<| go ( openCurveEnd crv ) crvs
 
     brushShape :: ptData -> SplinePts Closed
-    brushShape pt = brushFn ( ptParams pt )
+    brushShape pt = fun brushFn $ fun toBrushParams $ ptParams pt
 
     updateSpline :: ( T ( ℝ 2 ), T ( ℝ 2 ), T ( ℝ 2 ) ) -> ST s OutlineData
     updateSpline ( lastTgt, lastTgtFwd, lastTgtBwd )
@@ -313,17 +309,15 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart =
           ( \ ptData curve -> do
             ( prevTgt, prev_tgtFwd, prev_tgtBwd ) <- get
             let
-              fwd, bwd :: Double -> ( ℝ 2, T ( ℝ 2 ) )
+              fwdBwd :: OutlineFn
-              ( fwd, bwd ) = outlineFunctions @diffParams ptParams brushFn ptData curve
+              fwdBwd = outlineFunction ptParams toBrushParams brushFn ptData curve
               tgt, next_tgt, tgtFwd, next_tgtFwd, tgtBwd, next_tgtBwd :: T ( ℝ 2 )
-              tgt         = startTangent spt0 ptData curve
+              tgt      = startTangent spt0 ptData curve
-              next_tgt    =   endTangent spt0 ptData curve
+              next_tgt =   endTangent spt0 ptData curve
-              tgtFwd      =       snd ( fwd 0 )
+              ( ( _, tgtFwd      ), ( _, tgtBwd      ) ) = fwdBwd $ ℝ1 0
-              next_tgtFwd =       snd ( fwd 1 )
+              ( ( _, next_tgtFwd ), ( _, next_tgtBwd ) ) = fwdBwd $ ℝ1 1
-              tgtBwd      = -1 *^ snd ( bwd 1 )
-              next_tgtBwd = -1 *^ snd ( bwd 0 )
             lift $ tellBrushJoin ( prevTgt, prev_tgtFwd, tgtBwd ) ptData ( tgt, tgtFwd, prev_tgtBwd )
-            lift $ updateCurveData ( curveData curve ) fwd bwd
+            lift $ updateCurveData ( curveData curve ) fwdBwd
             put ( next_tgt, next_tgtFwd, next_tgtBwd )
           )
           ( const ( pure () ) )
@@ -331,10 +325,9 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart =
 
     updateCurveData
       :: crvData
-      -> ( Double -> ( ℝ 2, T ( ℝ 2 ) ) )
+      -> OutlineFn
-      -> ( Double -> ( ℝ 2, T ( ℝ 2 ) ) )
       -> WriterT OutlineData ( ST s ) ()
-    updateCurveData ( view ( typed @( CachedStroke s ) ) -> CachedStroke { cachedStrokeRef } ) fwd bwd = do
+    updateCurveData ( view ( typed @( CachedStroke s ) ) -> CachedStroke { cachedStrokeRef } ) fwdBwd = do
       mbOutline <- lift ( readSTRef cachedStrokeRef )
       case mbOutline of
         -- Cached fit data is available: use it.
@@ -346,12 +339,12 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart =
           let
             fwdData, bwdData :: ( SplinePts Open, Seq FitPoint )
             ( fwdData, bwdData ) =
-              ( fitSpline fitParams fwd, fitSpline fitParams bwd )
+              ( fitSpline fitParams ( fst . fwdBwd ), fitSpline fitParams ( snd . fwdBwd ) )
                 `Strats.using`
                   ( Strats.parTuple2 Strats.rdeepseq Strats.rdeepseq )
             outlineData :: OutlineData
-            outlineData = TwoSided fwdData bwdData
+            outlineData = TwoSided fwdData ( bimap reverseSpline Seq.reverse bwdData )
-          outlineData `deepseq` tell ( outlineData )
+          outlineData `deepseq` tell outlineData
           lift $ writeSTRef cachedStrokeRef ( Just outlineData )
 
     -- Connecting paths at a point of discontinuity of the tangent vector direction (G1 discontinuity).
@@ -416,107 +409,55 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart =
           $ joinWithBrush brush0 tgtBwd prevTgtBwd
 
 -- | Computes the forward and backward stroke outline functions for a single curve.
-outlineFunctions
+outlineFunction
-  :: forall diffParams brushParams crvData ptData
+  :: forall usedParams brushParams crvData ptData
-  .  ( Group  diffParams, Module Double diffParams
+  .  ( Interpolatable usedParams
-     , Torsor diffParams brushParams
+     , Diffy usedParams, Diffy brushParams
      , HasType ( ℝ 2 ) ptData
      -- Debugging.
-     , Show ptData
+     , Show ptData, Show brushParams
      )
-  => ( ptData -> brushParams )
+  => ( ptData -> usedParams )
-  -> ( brushParams -> SplinePts Closed )
+  -> ( usedParams ~> brushParams )
+  -> ( brushParams ~> SplinePts Closed )
   -> ptData
   -> Curve Open crvData ptData
-  -> ( Double -> ( ℝ 2, T ( ℝ 2 ) )
+  -> OutlineFn
-     , Double -> ( ℝ 2, T ( ℝ 2 ) )
+outlineFunction ptParams toBrushParams brushFromParams sp0 crv =
-     )
-outlineFunctions ptParams brushFn sp0 crv =
   let
-    p0 :: ℝ 2
+    usedParams :: ℝ 1 ~> usedParams
-    p0 = coords sp0
+    path :: ℝ 1 ~> ℝ 2
-    brush :: Double -> SplinePts Closed
+    ( path, usedParams ) =
-    f  :: Double -> ℝ 2
+      case crv of
-    f' :: Double -> T ( ℝ 2 )
+        LineTo    { curveEnd = NextPoint sp1 }
-    ( brush, f, f' ) = case crv of
+          | let seg = Segment sp0 sp1
-      LineTo    { curveEnd = NextPoint sp1 }
+          -> ( line ( fmap coords seg ), line ( fmap ptParams seg ) )
-        | let
+        Bezier2To { controlPoint = sp1, curveEnd = NextPoint sp2 }
-            p1 :: ℝ 2
+          | let bez2 = Quadratic.Bezier sp0 sp1 sp2
-            p1 = coords sp1
+          -> ( bezier2 ( fmap coords bez2 ), bezier2 ( fmap ptParams bez2 ) )
-            tgt :: T ( ℝ 2 )
+        Bezier3To { controlPoint1 = sp1, controlPoint2 = sp2, curveEnd = NextPoint sp3 }
-            tgt = p0 --> p1
+          | let bez3 = Cubic.Bezier sp0 sp1 sp2 sp3
-            brush1 :: Double -> SplinePts Closed
+          -> ( bezier3 ( fmap coords bez3 ), bezier3 ( fmap ptParams bez3 ) )
-            brush1 t = brushFn ( lerp @diffParams t ( ptParams sp0 ) ( ptParams sp1 ) )
+
-        -> ( brush1, \ t -> lerp @( T ( ℝ 2 ) ) t p0 p1, const tgt )
+    fwdBwd :: OutlineFn
-      Bezier2To { controlPoint = sp1, curveEnd = NextPoint sp2 }
+    fwdBwd t
-        | let
+      = solveEnvelopeEquations path_t path'_t ( fwdOffset, bwdOffset ) curves
-            p1, p2 :: ℝ 2
+    --  = ( ( offset fwdOffset • path_t,       path'_t )
-            p1 = coords sp1
+    --    , ( offset bwdOffset • path_t, -1 *^ path'_t ) )
-            p2 = coords sp2
-            bez :: Quadratic.Bezier ( ℝ 2 )
-            bez = Quadratic.Bezier {..}
-            brush2 :: Double -> SplinePts Closed
-            brush2 t =
-              brushFn $
-                Quadratic.bezier @diffParams
-                  ( Quadratic.Bezier ( ptParams sp0 ) ( ptParams sp1 ) ( ptParams sp2 ) ) t
-        -> ( brush2, Quadratic.bezier @( T ( ℝ 2 ) ) bez, Quadratic.bezier' bez )
-      Bezier3To { controlPoint1 = sp1, controlPoint2 = sp2, curveEnd = NextPoint sp3 }
-        | let
-            p1, p2, p3 :: ℝ 2
-            p1 = coords sp1
-            p2 = coords sp2
-            p3 = coords sp3
-            bez :: Cubic.Bezier ( ℝ 2 )
-            bez = Cubic.Bezier {..}
-            brush3 :: Double -> SplinePts Closed
-            brush3 t =
-              brushFn $
-                Cubic.bezier @diffParams
-                  ( Cubic.Bezier ( ptParams sp0 ) ( ptParams sp1 ) ( ptParams sp2 ) ( ptParams sp3 ) ) t
-        -> ( brush3, Cubic.bezier @( T ( ℝ 2 ) ) bez, Cubic.bezier' bez )
-    fwd, bwd :: Double -> ( ℝ 2, T ( ℝ 2 ) )
-    fwd t
-      = ( off t --offset ( withTangent ( fwd' t ) ( brush t ) ) • f t
-        , fwd' t
-        )
       where
-        off :: Double -> ℝ 2
+
-        off u = offset ( withTangent ( f' u ) ( brush u ) ) • f u
+        curves :: Seq ( ℝ 1 -> StrokeDatum )
-        offTgt :: Double -> T ( ℝ 2 )
+        curves = brushStrokeData path ( usedParams `chainRule` toBrushParams ) brushFromParams t
-        offTgt u
+
-          | u < 0.5
+        fwdOffset = withTangent         path'_t   brush_t
-          = 1e9 *^ ( off u --> off (u + 1e-9) )
+        bwdOffset = withTangent ( -1 *^ path'_t ) brush_t
-          | otherwise
+
-          = 1e9 *^ ( off (u - 1e-9) --> off u )
+        D1 path_t path'_t _ = runD path t
-        fwd' :: Double -> T ( ℝ 2 )
+        D1 params_t _ _     = runD usedParams t
-        fwd' u
+        brush_t             = value @brushParams
-          | squaredNorm ( offTgt u ) < epsilon
+                            $ runD brushFromParams ( fun toBrushParams params_t )
-          = f' u
+
-          | otherwise
+  in fwdBwd
-          = offTgt u
-    bwd t
-      = ( off s --offset ( withTangent ( -1 *^ bwd' s ) ( brush s ) ) • f s
-        , bwd' s
-        )
-      where
-        s :: Double
-        s = 1 - t
-        off :: Double -> ℝ 2
-        off u = offset ( withTangent ( -1 *^ f' u ) ( brush u ) ) • f u
-        offTgt :: Double -> T ( ℝ 2 )
-        offTgt u
-          | u < 0.5
-          = 1e9 *^ ( off u --> off (u + 1e-9) )
-          | otherwise
-          = 1e9 *^ ( off (u - 1e-9) --> off u )
-        bwd' :: Double -> T ( ℝ 2 )
-        bwd' u
-          | squaredNorm ( offTgt u ) < epsilon
-          = -1 *^ f' u
-          | otherwise
-          = offTgt u
-  in ( fwd, bwd )
 
 -----------------------------------
 -- Various utility functions
@@ -799,38 +740,6 @@ withTangent tgt_wanted spline@( Spline { splineStart } )
 
 --------------------------------------------------------------------------------
 
--- $brushes
---
--- You can compute the envelope equation for a brush stroke by using
--- the functions 'linear', 'bezier2' and 'bezier3' in conjunction with
--- the 'brushStroke' function, e.g.
---
--- > brushStroke ( bezier2 path ) ( uncurryD $ fmap bezier3 brush )
-
--- | Linear interpolation, as a differentiable function.
-linear :: forall b. ( Module Double ( T b ), Torsor ( T b ) b )
-       => Segment b -> ℝ 1 ~> b
-linear ( Segment a b ) = D \ ( ℝ1 t ) ->
-  D1 ( lerp @( T b ) t a b )
-     ( a --> b )
-     origin
-
--- | A quadratic Bézier curve, as a differentiable function.
-bezier2 :: forall b. ( Module Double ( T b ), Torsor ( T b ) b )
-        => Quadratic.Bezier b -> ℝ 1 ~> b
-bezier2 bez = D \ ( ℝ1 t ) ->
-  D1 ( Quadratic.bezier @( T b ) bez t )
-     ( Quadratic.bezier' bez t )
-     ( Quadratic.bezier'' bez )
-
--- | A cubic Bézier curve, as a differentiable function.
-bezier3 :: forall b. ( Module Double ( T b ), Torsor ( T b ) b )
-        => Cubic.Bezier b -> ℝ 1 ~> b
-bezier3 bez = D \ ( ℝ1 t ) ->
-  D1 ( Cubic.bezier @( T b ) bez t )
-     ( Cubic.bezier' bez t )
-     ( Cubic.bezier'' bez t )
-
 -- | A brush stroke, as described by the equation
 --
 -- \[ c(t,s) = p(t) + b(t,s) \]
@@ -839,15 +748,10 @@ bezier3 bez = D \ ( ℝ1 t ) ->
 --
 --   - \( p(t) \) is the path that the brush follows, and
 --   - \( b(t,s) \) is the brush shape, as it varies along the path.
-brushStroke :: ℝ 1 ~> ℝ 2 -- ^ stroke path \( p(t) \)
+brushStroke :: D ( ℝ 1 ) ( ℝ 2 ) -- ^ stroke path \( p(t) \)
-            -> ℝ 2 ~> ℝ 2 -- ^ brush \( b(t,s) \)
+            -> D ( ℝ 2 ) ( ℝ 2 ) -- ^ brush \( b(t,s) \)
-            -> ℝ 2 ~> ℝ 2
+            -> D ( ℝ 2 ) ( ℝ 2 )
-brushStroke ( D f_p ) ( D f_b ) = D \ ( ℝ2 t0 s0 ) ->
+brushStroke ( D1 p dpdt d2pdt2 ) ( D2 b dbdt dbds d2bdt2 d2bdtds d2bds2 ) =
-  let !( D1 p dpdt d2pdt2 )
-          = f_p ( ℝ1 t0 )
-      !( D2 b dbdt dbds d2bdt2 d2bdtds d2bds2 )
-          = f_b ( ℝ2 t0 s0 )
-  in
   D2 ( unT $ T p ^+^ T b )
      -- c = p + b
 
@@ -868,14 +772,16 @@ brushStroke ( D f_p ) ( D f_b ) = D \ ( ℝ2 t0 s0 ) ->
 --
 -- \[ \frac{\partial E}{\partial s} \frac{\mathrm{d} c}{\mathrm{d} t} \]
 --
--- whose roots correspond to cusps in the envelope.
+-- whose roots correspond to cusps in the envelope, and
-envelopeEquation :: ℝ 2 ~> ℝ 2 -> Double -> Double -> (# Double, T ( ℝ 2 ) #)
+--
-envelopeEquation ( D c ) t s =
+-- \[ \frac{\partial E}{\partial s}. \]
-  case c ( ℝ2 t s ) of
+envelopeEquation :: D ( ℝ 2 ) ( ℝ 2 ) -> ( Double, T ( ℝ 2 ), Double, Double )
-    D2 _ dcdt dcds d2cdt2 d2cdtds d2cds2 ->
+envelopeEquation ( D2 _ dcdt dcds d2cdt2 d2cdtds d2cds2 ) =
-      let dEdt = d2cdt2  `cross2` dcds + dcdt `cross2` d2cdtds
+  let ee = dcdt `cross` dcds
-          dEds = d2cdtds `cross2` dcds + dcdt `cross2` d2cds2
+      dEdt = d2cdt2  `cross` dcds + dcdt `cross` d2cdtds
-      in (# dcdt `cross2` dcds, dEds *^ dcdt ^-^ dEdt *^ dcds #)
+      dEds = d2cdtds `cross` dcds + dcdt `cross` d2cds2
+      tot = dcdt ^-^ ( dEdt / dEds ) *^ dcds --dEds *^ dcdt ^-^ dEdt *^ dcds
+  in ( ee, tot, dEdt, dEds )
   -- Computation of total derivative dc/dt:
   --
   -- dc/dt = ∂c/∂t + ∂c/∂s ∂s/∂t
@@ -883,7 +789,218 @@ envelopeEquation ( D c ) t s =
   --
   -- ( ∂E / ∂s ) dc/dt = ( ∂E / ∂s ) ∂c/∂t - ( ∂E / ∂t ) ∂c/∂s.
 
--- | Cross-product of two 2D vectors.
+
-cross2 :: T ( ℝ 2 ) -> T ( ℝ 2 ) -> Double
+-- | Linear interpolation, as a differentiable function.
-cross2 ( T ( ℝ2 x1 y1 ) ) ( T ( ℝ2 x2 y2 ) )
+line :: forall b. ( Module Double ( T b ), Torsor ( T b ) b )
-  = x1 * y2 - x2 * y1
+       => Segment b -> ℝ 1 ~> b
+line ( Segment a b ) = D \ ( ℝ1 t ) ->
+  D1 ( lerp @( T b ) t a b )
+     ( a --> b )
+     origin
+
+-- | A quadratic Bézier curve, as a differentiable function.
+bezier2 :: forall b. ( Module Double ( T b ), Torsor ( T b ) b )
+        => Quadratic.Bezier b -> ℝ 1 ~> b
+bezier2 bez = D \ ( ℝ1 t ) ->
+  D1 ( Quadratic.bezier @( T b ) bez t )
+     ( Quadratic.bezier'         bez t )
+     ( Quadratic.bezier''        bez   )
+
+-- | A cubic Bézier curve, as a differentiable function.
+bezier3 :: forall b. ( Module Double ( T b ), Torsor ( T b ) b )
+        => Cubic.Bezier b -> ℝ 1 ~> b
+bezier3 bez = D \ ( ℝ1 t ) ->
+  D1 ( Cubic.bezier @( T b ) bez t )
+     ( Cubic.bezier'         bez t )
+     ( Cubic.bezier''        bez t )
+
+splineCurveFns :: SplinePts Closed -> Seq ( ℝ 1 ~> ℝ 2 )
+splineCurveFns spls
+  = runIdentity
+  . bifoldSpline
+      ( \ pt crv -> Identity . Seq.singleton $ curveFn pt crv )
+      ( const $ Identity Seq.empty )
+  . adjustSplineType @Open
+  $ spls
+  where
+    curveFn :: ℝ 2 -> Curve Open () ( ℝ 2 ) -> ( ℝ 1 ~> ℝ 2 )
+    curveFn p0 = \case
+      LineTo    { curveEnd = NextPoint p1 }
+        -> line    $ Segment p0 p1
+      Bezier2To { controlPoint = p1, curveEnd = NextPoint p2 }
+        -> bezier2 $ Quadratic.Bezier p0 p1 p2
+      Bezier3To { controlPoint1 = p1, controlPoint2 = p2, curveEnd = NextPoint p3 }
+        -> bezier3 $ Cubic.Bezier p0 p1 p2 p3
+
+-- | Solve the envelope equations at a given point \( t = t_0 \), to find
+-- \( s_0 \) such that \( c(t_0, s_0) \) is on the envelope of the brush stroke.
+solveEnvelopeEquations :: ℝ 2
+                       -> T ( ℝ 2 )
+                       -> ( Offset, Offset )
+                       -> Seq ( ℝ 1 -> StrokeDatum )
+                       -> ( ( ℝ 2, T ( ℝ 2 ) ), ( ℝ 2, T ( ℝ 2 ) ) )
+solveEnvelopeEquations path_t path'_t ( fwdOffset, bwdOffset ) strokeData
+  = ( fwdSol, ( bwdPt, -1 *^ bwdTgt ) )
+
+  where
+
+    -- !_ = trace
+    --      ( unlines
+    --        [ "solveEnvelopeEquation"
+    --        , "  pt: " ++ show path_t
+    --        , "  tgt: "  ++ show path'_t
+    --        , "  fwdOffset: " ++ show fwdOffset
+    --        , "  bwdOffset: " ++ show bwdOffset ]
+    --      ) True
+
+    fwdSol = findSolFrom False fwdOffset
+    ( bwdPt, bwdTgt ) = findSolFrom True  bwdOffset
+
+    n :: Int
+    n = length strokeData
+
+    findSolFrom :: Bool -> Offset -> ( ℝ 2, T ( ℝ 2 ) )
+    findSolFrom goingBwds ( Offset { offsetIndex = i00, offsetParameter = s00, offset = off } )
+      = go ( Set.singleton i00 ) i00 ( fromMaybe 0.5 s00 )
+      where
+
+        plausibleTangent :: T ( ℝ 2 ) -> Bool
+        plausibleTangent tgt
+          | goingBwds
+          = path'_t ^.^ tgt < 0
+          | otherwise
+          = path'_t ^.^ tgt > 0
+
+        go :: Set Int -> Int -> Double -> ( ℝ 2, T ( ℝ 2 ) )
+        go seen i0 s0 =
+          case sol s0 ( strokeData `Seq.index` i0 ) of
+            ( good, s, pt, tgt )
+              | True --good && plausibleTangent tgt
+              -> ( pt, tgt )
+          --    | let ( i', s0' ) = mbNextPoint i0 ( unℝ1 s )
+          --    , not ( i' `Set.member` seen )
+          --    -> go ( Set.insert i' seen ) i' s0'
+              | otherwise
+              -> ( off • path_t, path'_t )
+
+        mbNextPoint :: Int -> Double -> ( Int, Double )
+        mbNextPoint i0 s0
+          | s0 <= 0.5
+          = ( prev i0, 0.9 )
+          | otherwise
+          = ( next i0, 0.1 )
+
+        prev, next :: Int -> Int
+        prev i
+          | i == 0
+          = n - 1
+          | otherwise
+          = i - 1
+        next i
+          | i == n - 1
+          = 0
+          | otherwise
+          = i + 1
+
+    sol :: Double -> ( ℝ 1 -> StrokeDatum ) -> ( Bool, ℝ 1, ℝ 2, T ( ℝ 2 ) )
+    sol initialGuess f =
+      let (good, s) = case newtonRaphson maxIters precision domain ( eqn f ) initialGuess of
+                        Nothing -> ( False, initialGuess )
+                        Just s0 -> ( True , s0 )
+      in case f ( ℝ1 s ) of -- TODO: a bit redundant to have to compute this again...
+          StrokeDatum { dstroke, dcdt } ->
+            ( good, ℝ1 s, value @( ℝ 2 ) dstroke, dcdt )
+
+
+    eqn :: ( ℝ 1 -> StrokeDatum ) -> ( Double -> ( Double, Double ) )
+    eqn f s =
+      case f ( ℝ1 s ) of
+        StrokeDatum { ee, 𝛿E𝛿s } ->
+          ( ee, 𝛿E𝛿s )
+
+    maxIters :: Word
+    maxIters = 5 --30
+    precision :: Int
+    precision = 10
+    domain :: ( Double, Double )
+    domain = ( 0, 1 )
+
+newtype ZipSeq a = ZipSeq { getZipSeq :: Seq a }
+  deriving stock Functor
+instance Applicative ZipSeq where
+  pure _ = error "no pure for ZipSeq"
+  liftA2 f ( ZipSeq xs ) ( ZipSeq ys ) = ZipSeq ( Seq.zipWith f xs ys )
+
+brushStrokeData :: forall brushParams
+                .  ( Diffy brushParams, Show brushParams )
+                => ( ℝ 1 ~> ℝ 2 )                      -- ^ path
+                -> ( ℝ 1 ~> brushParams )              -- ^ brush parameters
+                -> ( brushParams ~> SplinePts Closed ) -- ^ brush from parameters
+                -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum ) )
+brushStrokeData path params brush =
+  \ t ->
+    let
+        dpath_t :: D ( ℝ 1 ) ( ℝ 2 )
+        !dpath_t = runD path t
+        dparams_t :: D ( ℝ 1 ) brushParams
+        !dparams_t@( D1 { v = params_t } ) = runD params t
+        dbrush_params :: D brushParams ( SplinePts Closed )
+        !dbrush_params = runD brush params_t
+        splines :: Seq ( D brushParams ( ℝ 1 ~> ℝ 2 ) )
+        !splines = getZipSeq $ traverse @_ @ZipSeq ( ZipSeq . splineCurveFns ) dbrush_params
+        dbrushes_t :: Seq ( ℝ 1 -> D ( ℝ 2 ) ( ℝ 2 ) )
+        !dbrushes_t = force $ fmap ( uncurryD' . ( dparams_t `chain` ) ) splines
+
+    in fmap ( mkStrokeDatum dpath_t ) dbrushes_t
+  where
+
+    mkStrokeDatum :: D ( ℝ 1 ) ( ℝ 2 )
+                  -> ( ℝ 1 -> D ( ℝ 2 ) ( ℝ 2 ) )
+                  -> ( ℝ 1 -> StrokeDatum )
+    mkStrokeDatum dpath_t dbrush_t s =
+      let dbrush_t_s = dbrush_t s
+          dstroke@( D2 _c _𝛿c𝛿t _𝛿c𝛿s _ _ _ ) = brushStroke dpath_t dbrush_t_s
+          ( ee, dcdt, _𝛿E𝛿t, 𝛿E𝛿s ) = envelopeEquation dstroke
+      in  -- trace
+          --   ( unlines
+          --     [ "envelopeEquation:"
+          --     , "      s = " ++ show s
+          --     , "      c = " ++ show _c
+          --     , "  ∂c/∂t = " ++ show _𝛿c𝛿t
+          --     , "  ∂c/∂s = " ++ show _𝛿c𝛿s
+          --     , "      E = " ++ show ee
+          --     , "  ∂E/∂t = " ++ show _𝛿E𝛿t
+          --     , "  ∂E/∂s = " ++ show 𝛿E𝛿s
+          --     , "  dc/dt = " ++ show dcdt ] ) $
+          StrokeDatum
+            { dpath  = dpath_t
+            , dbrush = dbrush_t_s
+            , dstroke
+            , ee, dcdt, 𝛿E𝛿s }
+
+-- | The value and derivative of a brush stroke at a given coordinate
+-- \( (t_0, s_0) \), together with the value of the envelope equation at that
+-- point.
+data StrokeDatum
+  = StrokeDatum
+  { -- | Path \( p(t_0) \) (with its 1st and 2nd derivatives).
+    dpath    :: D ( ℝ 1 ) ( ℝ 2 )
+    -- | Brush shape \( b(t_0, s_0) \) (with its 1st and 2nd derivatives).
+  , dbrush   :: D ( ℝ 2 ) ( ℝ 2 )
+
+  -- Everything below can be computed in terms of the first two fields.
+
+    -- | Stroke \( c(t_0,s_0) = p(t_0) + b(t_0,s_0) \) (with its 1st and 2nd derivatives).
+  , dstroke  :: D ( ℝ 2 ) ( ℝ 2 )
+    -- | Envelope
+    --
+    -- \[ E(t_0,s_0) = \left ( \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} \right )_{(t_0,s_0)}. \]
+  , ee       :: Double
+    -- | Total derivative
+    --
+    -- \[ \left ( \frac{\mathrm{d} c}{\mathrm{d} t} \right )_{(t_0,s_0)}. \]
+  , dcdt     :: T ( ℝ 2 )
+    -- | \( \left ( \frac{\partial E}{\partial s} \right )_{(t_0,s_0)}. \)
+  , 𝛿E𝛿s     :: Double
+  }
+  deriving stock Show

src/splines/Math/Linear.hs

@@ -71,7 +71,7 @@ data family ℝ n
 data instance ℝ 0 = ℝ0
   deriving stock ( Show, Eq, Ord, Generic )
   deriving anyclass NFData
-newtype instance ℝ 1 = ℝ1 Double
+newtype instance ℝ 1 = ℝ1 { unℝ1 :: Double }
   deriving stock   ( Generic )
   deriving newtype ( Show, Eq, Ord, NFData )
 data    instance ℝ 2 = ℝ2 {-# UNPACK #-} !Double {-# UNPACK #-} !Double
@@ -86,8 +86,12 @@ data    instance ℝ 3 = ℝ3 {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-#
 -- | Tangent space to Euclidean space.
 type T :: Type -> Type
 newtype T e = T { unT :: e }
-  deriving stock ( Eq, Functor )
+  deriving stock ( Eq, Functor, Foldable, Traversable )
-  deriving newtype ( Show, NFData ) -- newtype Show instance for debugging...
+  deriving newtype NFData -- newtype Show instance for debugging...
+
+instance {-# OVERLAPPING #-} Show ( ℝ n ) => Show ( T ( ℝ n ) ) where
+  show ( T p ) = "V" ++ drop 1 ( show p )
+deriving stock instance Show v => Show ( T v )
 
 instance Applicative T where
   pure = T

src/splines/Math/Linear/Dual.hs

@@ -1,5 +1,6 @@
 {-# LANGUAGE AllowAmbiguousTypes   #-}
 {-# LANGUAGE DuplicateRecordFields #-}
+{-# LANGUAGE QuantifiedConstraints #-}
 {-# LANGUAGE ScopedTypeVariables   #-}
 {-# LANGUAGE UndecidableInstances  #-}
 
@@ -8,12 +9,12 @@ module Math.Linear.Dual where
 -- base
 import Control.Applicative
   ( liftA2 )
+import Data.Coerce
+  ( coerce )
 import Data.Kind
   ( Type, Constraint )
 import GHC.Generics
-  ( Generic, Generic1, Generically1(..) )
+  ( Generic, Generic1(..), Generically1(..) )
-import GHC.TypeNats
-  ( Nat )
 
 -- MetaBrush
 import Math.Module
@@ -24,8 +25,14 @@ import Math.Linear
 -- | Differentiable mappings between spaces.
 infixr 0 ~>
 type (~>) :: Type -> Type -> Type
-newtype a ~> b = D { runD :: a -> D a b }
+newtype u ~> v = D { runD :: u -> D u v }
-deriving stock instance Functor ( D a ) => Functor ( (~>) a )
+deriving stock instance Functor ( D u ) => Functor ( (~>) u )
+
+instance ( Applicative ( D u ), Module ( D u Double ) ( D u v ) ) => Module Double ( T ( u ~> v ) ) where
+  origin = T $ D \ _ -> origin
+  T ( D f ) ^+^ T ( D g ) = T $ D \ t -> f t ^+^ g t
+  T ( D f ) ^-^ T ( D g ) = T $ D \ t -> f t ^-^ g t
+  a *^ T ( D f ) = T $ D \ t -> pure a *^ f t
 
 -- | @D ( ℝ n ) v@ is \( \mathbb{R}[x_1, \ldots, x_n]/(x_1, \ldots, x_n)^3 \otimes_\mathbb{R} v \)
 type D :: Type -> Type -> Type
@@ -36,22 +43,25 @@ type instance D ( ℝ 2 ) = Dℝ2
 type instance D ( ℝ 3 ) = Dℝ3
 
 newtype Dℝ0 v = D0 { v :: v }
-  deriving stock   ( Show, Eq, Functor, Generic, Generic1 )
+  deriving stock   ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
   deriving newtype ( Num, Fractional, Floating )
   deriving Applicative
     via Generically1 Dℝ0
-data Dℝ1 v = D1 { v :: !v, dx :: !( T v ), ddx :: !( T v ) }
+data Dℝ1 v = D1 { v :: !v, df :: !( T v ), d²f :: !( T v ) }
-  deriving stock ( Show, Eq, Functor, Generic, Generic1 )
+  deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 )
   deriving Applicative
     via Generically1 Dℝ1
+deriving stock instance ( Show v, Show ( T v ) ) => Show ( Dℝ1 v )
 data Dℝ2 v = D2 { v :: !v, dx, dy :: !( T v ), ddx, dxdy, ddy :: !( T v ) }
-  deriving stock ( Show, Eq, Functor, Generic, Generic1 )
+  deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 )
   deriving Applicative
     via Generically1 Dℝ2
+deriving stock instance ( Show v, Show ( T v ) ) => Show ( Dℝ2 v )
 data Dℝ3 v = D3 { v :: !v, dx, dy, dz :: !( T v ), ddx, dxdy, ddy, dxdz, dydz, ddz :: !( T v ) }
-  deriving stock ( Show, Eq, Functor, Generic, Generic1 )
+  deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 )
   deriving Applicative
     via Generically1 Dℝ3
+deriving stock instance ( Show v, Show ( T v ) ) => Show ( Dℝ3 v )
 
 instance Num ( Dℝ1 Double ) where
   (+) = liftA2 (+)
@@ -110,29 +120,29 @@ instance Num ( Dℝ3 Double ) where
          ( T $ dz1 * dz2 + v1 * ddz2 + ddz1 * v2)
 
 
-instance Module Double v => Module ( Dℝ0 Double ) ( Dℝ0 v ) where
+instance Module Double ( T v ) => Module ( Dℝ0 Double ) ( Dℝ0 v ) where
-  (^+^)  = liftA2 (^+^)
+  (^+^)  = liftA2 ( coerce $ (^+^)  @Double @( T v ) )
-  (^-^)  = liftA2 (^-^)
+  (^-^)  = liftA2 ( coerce $ (^-^)  @Double @( T v ) )
-  origin = pure origin
+  origin = pure   ( coerce $ origin @Double @( T v ) )
-  (*^)   = liftA2 (*^)
+  (*^)   = liftA2 ( coerce $ (*^)   @Double @( T v ) )
 
-instance Module Double v => Module ( Dℝ1 Double ) ( Dℝ1 v ) where
+instance Module Double ( T v ) => Module ( Dℝ1 Double ) ( Dℝ1 v ) where
-  (^+^)  = liftA2 (^+^)
+  (^+^)  = liftA2 ( coerce $ (^+^)  @Double @( T v ) )
-  (^-^)  = liftA2 (^-^)
+  (^-^)  = liftA2 ( coerce $ (^-^)  @Double @( T v ) )
-  origin = pure origin
+  origin = pure   ( coerce $ origin @Double @( T v ) )
-  (*^)   = liftA2 (*^)
+  (*^)   = liftA2 ( coerce $ (*^)   @Double @( T v ) )
 
-instance Module Double v => Module ( Dℝ2 Double ) ( Dℝ2 v ) where
+instance Module Double ( T v ) => Module ( Dℝ2 Double ) ( Dℝ2 v ) where
-  (^+^)  = liftA2 (^+^)
+  (^+^)  = liftA2 ( coerce $ (^+^)  @Double @( T v ) )
-  (^-^)  = liftA2 (^-^)
+  (^-^)  = liftA2 ( coerce $ (^-^)  @Double @( T v ) )
-  origin = pure origin
+  origin = pure   ( coerce $ origin @Double @( T v ) )
-  (*^)   = liftA2 (*^)
+  (*^)   = liftA2 ( coerce $ (*^)   @Double @( T v ) )
 
-instance Module Double v => Module ( Dℝ3 Double ) ( Dℝ3 v ) where
+instance Module Double ( T v ) => Module ( Dℝ3 Double ) ( Dℝ3 v ) where
-  (^+^)  = liftA2 (^+^)
+  (^+^)  = liftA2 ( coerce $ (^+^)  @Double @( T v ) )
-  (^-^)  = liftA2 (^-^)
+  (^-^)  = liftA2 ( coerce $ (^-^)  @Double @( T v ) )
-  origin = pure origin
+  origin = pure   ( coerce $ origin @Double @( T v ) )
-  (*^)   = liftA2 (*^)
+  (*^)   = liftA2 ( coerce $ (*^)   @Double @( T v ) )
 
 instance Fractional ( Dℝ1 Double ) where
   (/) = error "I haven't yet defined (/) for Dℝ1"
@@ -204,48 +214,67 @@ instance Floating ( Dℝ3 Double ) where
 --------------------------------------------------------------------------------
 
 uncurryD :: ( ℝ 1 ~> ℝ 1 ~> b ) -> ( ℝ 2 ~> b )
-uncurryD ( D b ) = D \ ( ℝ2 t0 s0 ) ->
+uncurryD ( D b ) = D \ ( ℝ2 t0 s0 ) -> uncurryD' ( b ( ℝ1 t0 ) ) ( ℝ1 s0 )
-  let !(D1 ( D b_t0 ) ( T ( D dbdt_t0 ) ) ( T ( D d2bdt2_t0 ) ) ) = b ( ℝ1 t0 )
+
-      !(D1 b_t0s0 dbds_t0s0 d2bds2_t0s0 ) = b_t0 ( ℝ1 s0 )
+uncurryD' :: D ( ℝ 1 ) ( ℝ 1 ~> b ) -> ℝ 1 -> D ( ℝ 2 ) b
-      !(D1 dbdt_t0s0 d2bdtds_t0s0 _ ) = dbdt_t0 ( ℝ1 s0 )
+uncurryD' ( D1 ( D b_t0 ) ( T ( D dbdt_t0 ) ) ( T ( D d2bdt2_t0 ) ) ) ( ℝ1 s0 ) =
-      !(D1 d2bdt2_t0s0 _ _ ) = d2bdt2_t0 ( ℝ1 s0 )
+  let !( D1 b_t0s0 dbds_t0s0 d2bds2_t0s0 ) = b_t0 ( ℝ1 s0 )
+      !( D1 dbdt_t0s0 d2bdtds_t0s0 _ ) = dbdt_t0 ( ℝ1 s0 )
+      !( D1 d2bdt2_t0s0 _ _ ) = d2bdt2_t0 ( ℝ1 s0 )
   in D2 b_t0s0 ( T dbdt_t0s0 ) dbds_t0s0 ( T d2bdt2_t0s0 ) d2bdtds_t0s0 d2bds2_t0s0
 
 --------------------------------------------------------------------------------
 
+type Diffy :: Type -> Constraint
+class ( Applicative ( D v )
+      , Traversable ( D v )
+      , Module Double ( T v ) )
+   => Diffy v where
+  chain :: ( Module Double ( T w ) )
+        => D ( ℝ 1 ) v -> D v w -> D ( ℝ 1 ) w
+  konst :: Module Double ( T w ) => w -> D v w
+  value :: D v w -> w
+  linear :: Module Double ( T w ) => ( v -> w ) -> ( v ~> w )
+
 chainRule :: ( Diffy v, Module Double ( T w ) )
-          => ( ( ℝ 1 ) ~> v )
+          => ( ( ℝ 1 ) ~> v ) -> ( v ~> w ) -> ( ( ℝ 1 ) ~> w )
-          -> ( v ~> w )
+chainRule ( D df ) ( D dg ) =
-          -> ( ( ℝ 1 ) ~> w )
-chainRule ( D f ) ( D g ) =
   D \ x ->
-    case f x of
+    case df x of
-      df@( D1 { v = f_x } ) ->
+      df_x@( D1 { v = f_x } ) ->
-        chain df ( g f_x )
+        chain df_x ( dg f_x )
 
 -- | Recover the underlying function, discarding all infinitesimal information.
 fun :: forall v w. Diffy v => ( v ~> w ) -> ( v -> w )
-fun ( D f ) = value @v . f
+fun ( D df ) = value @v . df
+{-# INLINE fun #-}
 
-type Diffy :: Type -> Constraint
+var :: ( Representable u, Diffy u ) => Fin ( Dim u ) -> ( u ~> Double )
-class Diffy v where
+var i = linear ( `index` i )
-  chain :: ( Module Double ( T w ) )
+{-# INLINE var #-}
-        => Dℝ1 v -> D v w -> D ( ℝ 1 ) w
-  konst :: Module Double ( T w ) => w -> D v w
-  value :: D v w -> w
 
 instance Diffy ( ℝ 0 ) where
   chain _ ( D0 w ) = D1 w origin origin
+  {-# INLINE chain #-}
   konst k = D0 k
+  {-# INLINE konst #-}
   value ( D0 w ) = w
+  {-# INLINE value #-}
+  linear f = D \ _ -> D0 ( f ℝ0 )
+  {-# INLINE linear #-}
 
 instance Diffy ( ℝ 1 ) where
   chain ( D1 _ ( T ( ℝ1 x' ) ) ( T ( ℝ1 x'' ) ) ) ( D1 v g_x g_xx )
     = D1 v
          ( x' *^ g_x )
          ( x'' *^ g_x ^+^ ( x' * x' ) *^ g_xx )
+  {-# INLINE chain #-}
   konst k = D1 k origin origin
+  {-# INLINE konst #-}
   value ( D1 { v } ) = v
+  {-# INLINE value #-}
+  linear f = D \ u -> D1 ( f u ) ( T $ f u ) origin
+  {-# INLINE linear #-}
 
 instance Diffy ( ℝ 2 ) where
   chain ( D1 _ ( T ( ℝ2 x' y' ) ) ( T ( ℝ2 x'' y'' ) ) ) ( D2 v g_x g_y g_xx g_xy g_yy )
@@ -254,8 +283,16 @@ instance Diffy ( ℝ 2 ) where
          ( x'' *^ g_x ^+^ y'' *^ g_y
             ^+^ ( x' * x' ) *^ g_xx ^+^ ( y' * y' ) *^ g_yy
             ^+^ 2 *^ ( ( x' * y' ) *^ g_xy ) )
+  {-# INLINE chain #-}
   konst k = D2 k origin origin origin origin origin
+  {-# INLINE konst #-}
   value ( D2 { v } ) = v
+  {-# INLINE value #-}
+  linear f = D \ u ->
+    D2 ( f u )
+       ( T $ f ( ℝ2 1 0 ) ) ( T $ f ( ℝ2 0 1 ) )
+       origin origin origin
+  {-# INLINE linear #-}
 
 instance Diffy ( ℝ 3 ) where
   chain ( D1 _ ( T ( ℝ3 x' y' z' ) ) ( T ( ℝ3 x'' y'' z'' ) ) )
@@ -265,24 +302,13 @@ instance Diffy ( ℝ 3 ) where
         (       x'' *^ g_x ^+^ y'' *^ g_y ^+^ z'' *^ g_z
             ^+^ ( x' * x' ) *^ g_xx ^+^ ( y' * y' ) *^ g_yy ^+^ ( z' * z' ) *^ g_zz
             ^+^ 2 *^ ( ( x' * y' ) *^ g_xy ) ^+^ ( x' * z' ) *^ g_xz ^+^ ( y' * z' ) *^ g_yz )
+  {-# INLINE chain #-}
   konst k = D3 k origin origin origin origin origin origin origin origin origin
+  {-# INLINE konst #-}
   value ( D3 { v } ) = v
-
+  {-# INLINE value #-}
---------------------------------------------------------------------------------
+  linear f = D \ u ->
-
+    D3 ( f u )
-type Var :: Nat -> Type -> Constraint
+       ( T $ f ( ℝ3 1 0 0 ) ) ( T $ f ( ℝ3 0 1 0 ) ) ( T $ f ( ℝ3 0 0 1 ) )
-class Var n v where
+       origin origin origin origin origin origin
-  var :: v ~> Double
+  {-# INLINE linear #-}
-
-instance Var 1 ( ℝ 1 ) where
-  var = D \ ( ℝ1 x ) -> D1 x ( T 1 ) ( T 0 )
-instance Var 1 ( ℝ 2 ) where
-  var = D \ ( ℝ2 x _ ) -> D2 x ( T 1 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 )
-instance Var 2 ( ℝ 2 ) where
-  var = D \ ( ℝ2 _ y ) -> D2 y ( T 0 ) ( T 1 ) ( T 0 ) ( T 0 ) ( T 0 )
-instance Var 1 ( ℝ 3 ) where
-  var = D \ ( ℝ3 x _ _ ) -> D3 x ( T 1 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 )
-instance Var 2 ( ℝ 3 ) where
-  var = D \ ( ℝ3 _ y _ ) -> D3 y ( T 0 ) ( T 1 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 )
-instance Var 3 ( ℝ 3 ) where
-  var = D \ ( ℝ3 _ _ z ) -> D3 z ( T 0 ) ( T 0 ) ( T 1 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 )

src/splines/Math/Module.hs

@@ -5,6 +5,7 @@
 module Math.Module
   ( Module(..), lerp
   , Inner(..)
+  , Interpolatable
   , norm, squaredNorm, quadrance, distance
   , proj, projC, closestPointOnSegment
   , cross
@@ -17,6 +18,8 @@ import Control.Applicative
   ( liftA2 )
 import Control.Monad
   ( guard )
+import Data.Kind
+  ( Type, Constraint )
 import Data.Monoid
   ( Ap(..), Sum(..) )
 
@@ -109,6 +112,13 @@ closestPointOnSegment c ( Segment p0 p1 )
 
 --------------------------------------------------------------------------------
 
+-- | A convenient constraint synonym for types that support interpolation.
+type Interpolatable :: Type -> Constraint
+class    ( Torsor ( T r ) r, Module Double ( T r ) ) => Interpolatable r
+instance ( Torsor ( T r ) r, Module Double ( T r ) ) => Interpolatable r
+
+--------------------------------------------------------------------------------
+
 instance Num a => Module a ( Sum a ) where
 
   origin = Sum 0

src/splines/Math/Roots.hs

@@ -3,7 +3,18 @@
 
 {-# OPTIONS_GHC -Wno-name-shadowing #-}
 
-module Math.Roots where
+module Math.Roots
+  ( -- * Quadratic solver
+    solveQuadratic
+
+    -- * Laguerre method
+  , roots, realRoots
+  , laguerre, eval, derivative
+
+    -- * Newton–Raphson
+  , newtonRaphson
+  )
+  where
 
 -- base
 import Control.Monad
@@ -72,6 +83,9 @@ solveQuadratic a0 a1 a2
     disc :: a
     disc = a1 * a1 - 4 * a0 * a2
 
+--------------------------------------------------------------------------------
+-- Root finding using Laguerre's method
+
 -- | Find real roots of a polynomial with real coefficients.
 --
 -- Coefficients are given in order of decreasing degree, e.g.:
@@ -279,118 +293,140 @@ derivative p = do
   pure p'
 
 --------------------------------------------------------------------------------
-
+-- Newton–Raphson
-data NewtonRaphson
-  = NewtonRaphson
-    { f :: Double -> (# Double, Double #)
-    , maxIters     :: Word
-    , min_x, max_x :: Double
-    , digits       :: Int }
 
 -- Newton–Raphson implementation taken from Boost C++ library:
 -- https://github.com/boostorg/math/blob/0dc6a70caa6bbec2b6ae25eede36c430f0ccae13/include/boost/math/tools/roots.hpp#L217
-
+{-# SPECIALISE
-newtonRaphson :: NewtonRaphson
+  newtonRaphson
-              -> Double
+    :: Word -> Int -> ( Double, Double ) -> ( Double -> ( Double, Double ) ) -> Double -> Maybe Double
-              -> Maybe Double
+  #-}
-newtonRaphson ( NewtonRaphson { f, maxIters, min_x, max_x, digits } ) x0 =
+{-# INLINEABLE newtonRaphson #-}
-    go $
+newtonRaphson :: ( RealFloat r, Show r )
-      NewtonRaphsonState
+              => Word              -- ^ maximum number of iterations
-        { f_x_prev = 0
+              -> Int               -- ^ desired digits of precision
-        , x = x0
+              -> ( r, r )          -- ^ @(min_x, max_x)@.
-        , δ = maxDouble, δ1 = maxDouble
+              -> ( r -> ( r, r ) ) -- ^ function and its derivative
-        , iters = 0
+              -> r                 -- ^ initial guess
-        , min_x, max_x
+              -> Maybe r
-        , min_f_x = 0, max_f_x = 0 }
+newtonRaphson maxIters digits ( min_x, max_x ) f x0 =
+  doNewtonRaphson f maxIters factor min_x max_x 0 0 0 x0 maxRealFloat maxRealFloat
   where
-    factor = encodeFloat 1 ( 1 - digits )
+    !factor = encodeFloat 1 ( 1 - digits )
-    go ( NewtonRaphsonState { f_x_prev, x, δ = δ1, δ1 = δ2, iters, min_x, max_x, min_f_x, max_f_x } )
+
-      = case f x of
+doNewtonRaphson :: ( Fractional r, Ord r, Show r )
-          (# f_x, f'_x #)
+                => ( r -> (r, r) )
-            | f_x == 0
+                -> Word
+                -> r
+                -> r -> r
+                -> r -> r
+                -> r
+                -> r
+                -> r -> r
+                -> Maybe r
+doNewtonRaphson f maxIters factor min_x max_x min_f_x max_f_x f_x_prev x δ1 δ2 =
+  go min_x max_x min_f_x max_f_x f_x_prev 0 x δ1 δ2
+  where
+  go min_x max_x min_f_x max_f_x f_x_prev !iters !x !δ1 !δ2 =
+    case f x of
+      ( f_x, f'_x )
+        | f_x == 0
+        -> Just x
+        | ( new_x, δ, δ1 ) <- newtonRaphsonStep f min_x max_x f_x_prev x f_x f'_x δ1 δ2
+        -> if
+            |  abs δ <= abs ( new_x * factor )
+            || new_x == min_x || new_x == max_x
             -> Just x
-            | δ <-
+            | iters >= maxIters
-               if | f'_x == 0
+            -> Nothing
-                  -> handleZeroDerivative f_x_prev x f_x δ1 min_x max_x
+            | ( min_x, max_x, min_f_x, max_f_x ) <-
-                  | otherwise
+                 if δ > 0
-                  -> f_x / f'_x
+                 then ( min_x, x, min_f_x, f_x )
-            , (# δ, δ1 #) <-
+                 else ( x, max_x, f_x, max_f_x )
-                if | abs ( δ * δ ) > abs δ2
+            -> if min_f_x * max_f_x > 0
-                   , let shift = if δ > 0 then 0.5 * ( x - min_x ) else 0.5 * ( x - max_x )
+               then Nothing
-                   , δ <- if x /= 0 && ( abs shift > abs x )
+               else
-                          then signum δ * abs x * 1.1
+                 go min_x max_x min_f_x max_f_x f_x ( iters + 1 ) new_x δ δ1
-                          else shift
-                   -> (# δ, 3 * δ #)
-                   | otherwise
-                   -> (# δ, δ1 #)
-            , let new_x = x - δ
-            , (# δ, new_x #) <-
-               if | new_x <= min_x
-                  , δ <- 0.5 * ( x - min_x )
-                  -> (# δ, x - δ #)
-                  | new_x >= max_x
-                  , δ <- 0.5 * ( x - max_x )
-                  -> (# δ, x - δ #)
-                  | otherwise
-                  -> (# δ, new_x #)
-            -> if
-                |  abs δ <= abs ( new_x * factor )
-                || new_x == min_x || new_x == max_x
-                -> Just x
-                | iters >= maxIters
-                -> Nothing
-                | (# min_x, max_x, min_f_x, max_f_x #) <-
-                     if δ > 0
-                     then (# min_x, x, min_f_x, f_x #)
-                     else (# x, max_x, f_x, max_f_x #)
-                -> if min_f_x * max_f_x > 0
-                   then Nothing
-                   else
-                    go $ NewtonRaphsonState
-                         { f_x_prev = f_x, x = new_x
-                         , δ, δ1
-                         , iters = iters + 1
-                         , min_x, max_x
-                         , min_f_x, max_f_x
-                         }
 
-    handleZeroDerivative :: Double
+newtonRaphsonStep :: ( Fractional r, Ord r, Show r )
-                         -> Double -> Double
+                  => ( r -> ( r, r ) )
-                         -> Double
+                  -> r
-                         -> Double -> Double
+                  -> r -> r -> r
-                         -> Double
+                  -> r -> r
-    handleZeroDerivative f_x_prev x f_x δ min_x max_x
+                  -> r -> r
-      -- Handle zero derivative on first iteration.
+                  -> ( r, r, r )
-      | f_x_prev == 0
+newtonRaphsonStep f min_x max_x f_x_prev x f_x f'_x δ1 δ2
-      , x_prev <- if x <= 0.5 * ( min_x + max_x ) then max_x else min_x
+  | δ <-
-      , (# f_x_prev, _ #) <- f x_prev
+     if | f'_x == 0
-      , δ <- x_prev - x
+        -> zeroDerivativeStep f min_x max_x f_x_prev x f_x δ1
-      = finish f_x_prev δ
+        | otherwise
-      | otherwise
+        -> f_x / f'_x
-      = finish f_x_prev δ
+  , ( δ, δ1 ) <-
+     if | abs ( δ * δ ) > abs δ2
+        -> bisectStep min_x max_x x δ
+        | otherwise
+        -> ( δ, δ1 )
+  , let !new_x = x - δ
+  , ( δ, new_x ) <-
+     if | new_x <= min_x
+        , δ <- 0.5 * ( x - min_x )
+        -> ( δ, x - δ )
+        | new_x >= max_x
+        , δ <- 0.5 * ( x - max_x )
+        -> ( δ, x - δ )
+        | otherwise
+        -> ( δ, new_x )
+  = ( new_x, δ, δ1 )
 
-      where
+-- Handle \( f'(x_0) = 0 \).
-        finish f_x_prev δ
+zeroDerivativeStep :: ( Fractional r, Ord r, Show r )
-          | signum f_x_prev * signum f_x < 0
+                   => ( r -> ( r, r ) )
-          = if δ < 0 then 0.5 * ( x - min_x ) else 0.5 * ( x - max_x )
+                   -> r
-          | otherwise
+                   -> r -> r
-          = if δ < 0 then 0.5 * ( x - max_x ) else 0.5 * ( x - min_x )
+                   -> r
-
+                   -> r -> r
--- | Loop state for the 'newtonRaphson' function.
+                   -> r
-data NewtonRaphsonState =
+zeroDerivativeStep f min_x max_x f_x_prev x f_x δ
-  NewtonRaphsonState
+  -- Handle zero derivative on first iteration.
-    { f_x_prev         :: !Double
+  | f_x_prev == 0
-    , x                :: !Double
+  , x_prev <- if x <= 0.5 * ( min_x + max_x ) then max_x else min_x
-    , δ, δ1            :: !Double
+  , ( f_x_prev, _ ) <- f x_prev
-    , iters            :: !Word
+  , δ <- x_prev - x
-    , min_x, max_x     :: !Double
+  = finish f_x_prev δ
-    , min_f_x, max_f_x :: !Double }
+  | otherwise
-
+  = finish f_x_prev δ
-maxDouble :: Double
-maxDouble = encodeFloat m n
   where
-    b = floatRadix ( 0 :: Double )
+    finish f_x_prev δ
-    e = floatDigits ( 0 :: Double )
+      | signum f_x_prev * signum f_x < 0
-    (_, e') = floatRange ( 0 :: Double )
+      = if δ < 0 then 0.5 * ( x - min_x ) else 0.5 * ( x - max_x )
+      | otherwise
+      = if δ < 0 then 0.5 * ( x - max_x ) else 0.5 * ( x - min_x )
+
+-- Handle failure of convergence in the past two steps.
+bisectStep :: ( Fractional r, Ord r, Show r ) => r -> r -> r -> r -> ( r, r )
+bisectStep min_x max_x x δ
+  | let !shift = if δ > 0 then 0.5 * ( x - min_x ) else 0.5 * ( x - max_x )
+  , δ <- if x /= 0 && ( abs shift > abs x )
+         then signum δ * abs x * 1.1
+         else shift
+  = ( δ, 3 * δ )
+
+{-# SPECIALISE maxRealFloat :: Float #-}
+{-# SPECIALISE maxRealFloat :: Double #-}
+maxRealFloat :: forall r. RealFloat r => r
+maxRealFloat = encodeFloat m n
+  where
+    b = floatRadix @r 0
+    e = floatDigits @r 0
+    (_, e') = floatRange @r 0
     m = b ^ e - 1
     n = e' - e
+
+{-
+
+for some reason GHC isn't able to optimise maxRealFloat @Double into
+
+maxDouble :: Double
+maxDouble = D# 1.7976931348623157e308##
+
+-}