Add mechanisms to log envelope equation data

brush-strokes/.gitignore

@@ -0,0 +1,3 @@
+dist-newstyle/
+logs/
+cabal.project.local

brush-strokes/brush-strokes.cabal

@@ -12,7 +12,7 @@ description:
 flag use-fma
   description: Use fused-muliply add instructions to implement interval arithmetic.
   default: True
-  manual:  True
+  manual:  False
 
 common common
 
@@ -27,6 +27,10 @@ common common
            >= 0.6.0.1 && < 0.8
       , deepseq
            >= 1.4.4.0 && < 1.6
+      , directory
+           >= 1.3.7.1 && < 1.4
+      , filepath
+           >= 1.4     && < 1.6
       , generic-lens
            >= 2.2     && < 2.3
       , groups
@@ -37,8 +41,12 @@ common common
           ^>= 0.9.0.0
       , rounded-hw
           ^>= 0.4
+      , time
+          ^>= 1.12.2
       , transformers
            >= 0.5.6.2 && < 0.7
+      , tree-view
+          ^>= 0.5
 
     default-extensions:
         BangPatterns

brush-strokes/src/Debug/Utils.hs

@@ -1,21 +1,44 @@
-module Debug.Utils ( trace ) where
+module Debug.Utils
+  ( trace
+  , logToFile
+  ) where
 
 -- base
+import Control.Concurrent.MVar
+  ( MVar, withMVarMasked )
+import Control.Monad
+  ( void )
 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 )
 
+-- deepseq
+import Control.DeepSeq
+  ( force )
+
+-- directory
+import qualified System.Directory as Directory
+  ( createDirectoryIfMissing )
+
+-- filepath
+import qualified System.FilePath as FilePath
+  ( takeDirectory )
+
+-- time
+import qualified Data.Time.Clock as Time
+  ( getCurrentTime )
+import qualified  Data.Time.Format as Time
+  ( defaultTimeLocale, formatTime )
+
 --------------------------------------------------------------------------------
 
+-- | Like 'Debug.Trace.trace', but using 'withCP65001` in order to handle
+-- Unicode without crashing on Windows.
 trace :: String -> a -> a
 trace ( force -> !str ) expr =
   unsafePerformIO $ withCP65001 do
@@ -23,3 +46,23 @@ trace ( force -> !str ) expr =
     hPutStrLn stdout str
     hFlush stdout
     return expr
+
+-- | Log the second argument to the file stored in the 'MVar' in the first
+-- argument.
+--
+-- The 'MVar' is used to avoid jumbled contents when attempting to write to
+-- the file concurrently.
+logToFile :: MVar FilePath -> String -> ()
+logToFile logFileMVar logContents =
+  unsafePerformIO $ withCP65001 $ void $ withMVarMasked logFileMVar $ \ logFile -> do
+    now <- Time.getCurrentTime
+    let timeString = Time.formatTime Time.defaultTimeLocale "%0Y-%m-%d (%Hh %Mm %Ss)" now
+        logContentsWithHeader = unlines
+          [ replicate 80 '='
+          , timeString
+          , logContents
+          ]
+    Directory.createDirectoryIfMissing True $
+      FilePath.takeDirectory logFile
+    appendFile logFile logContentsWithHeader
+    return logFile

brush-strokes/src/Math/Algebra/Dual.hs

@@ -33,7 +33,7 @@ import Data.Kind
 import GHC.TypeNats
   ( Nat )
 
--- MetaBrush
+-- brush-strokes
 import Math.Algebra.Dual.Internal
 import Math.Interval.Internal
 import Math.Linear

brush-strokes/src/Math/Algebra/Dual/Internal.hs

@@ -31,7 +31,7 @@ import Language.Haskell.TH
 import Language.Haskell.TH.Syntax
   ( Lift(..) )
 
--- MetaBrush
+-- brush-strokes
 import Math.Linear
   ( Vec(..), T(..)
   , RepresentableQ(..), RepDim

brush-strokes/src/Math/Bezier/Cubic.hs

@@ -50,7 +50,7 @@ import Data.Group.Generics
 import Data.Primitive.Types
   ( Prim )
 
--- MetaBrush
+-- brush-strokes
 import qualified Math.Bezier.Quadratic as Quadratic
   ( Bezier(..), bezier )
 import Math.Epsilon
@@ -99,6 +99,7 @@ fromQuadratic ( Quadratic.Bezier { p0 = q0, p1 = q1, p2 = q2 } ) = Bezier {..}
     p1 = lerp @v (2/3) q0 q1
     p2 = lerp @v (1/3) q1 q2
     p3 = q2
+{-# INLINEABLE fromQuadratic #-}
 
 -- | Cubic Bézier curve.
 bezier :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> p
@@ -106,17 +107,20 @@ bezier ( Bezier {..} ) t =
   lerp @v t
     ( Quadratic.bezier @v ( Quadratic.Bezier p0 p1 p2 ) t )
     ( Quadratic.bezier @v ( Quadratic.Bezier p1 p2 p3 ) t )
+{-# INLINEABLE bezier #-}
 
 -- | The derivative of a Cubic Bézier curve, as a quadratic Bézier curve.
 derivative :: ( Group v, Module r v ) => Bezier v -> Quadratic.Bezier v
 derivative ( Bezier {..} ) = ( Ring.fromInteger 3 *^ )
                           <$> Quadratic.Bezier ( p0 --> p1 ) ( p1 --> p2 ) ( p2 --> p3 )
+{-# INLINEABLE derivative #-}
 
 -- | Derivative of a cubic Bézier curve.
 bezier' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> v
 bezier' ( Bezier {..} )
   = ( Ring.fromInteger 3 *^ )
   . Quadratic.bezier @v ( Quadratic.Bezier ( p0 --> p1 ) ( p1 --> p2 ) ( p2 --> p3 ) )
+{-# INLINEABLE bezier' #-}
 
 -- | Second derivative of a cubic Bézier curve.
 bezier'' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> v
@@ -125,16 +129,19 @@ bezier'' ( Bezier {..} ) t
   $ lerp @v t
       ( p1 --> p0 ^+^ p1 --> p2 )
       ( p2 --> p1 ^+^ p2 --> p3 )
+{-# INLINEABLE bezier'' #-}
 
 -- | Third derivative of a cubic Bézier curve.
 bezier''' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> v
 bezier''' ( Bezier {..} )
   = ( Ring.fromInteger 6 *^ )
   $ ( ( p0 --> p3 ) ^+^ Ring.fromInteger 3 *^ ( p2 --> p1 ) )
+{-# INLINEABLE bezier''' #-}
 
 -- | Curvature of a cubic Bézier curve.
 curvature :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> r -> r
 curvature bez t = sqrt $ squaredCurvature @v bez t
+{-# INLINEABLE curvature #-}
 
 -- | Square of curvature of a cubic Bézier curve.
 squaredCurvature :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> r -> r
@@ -150,6 +157,7 @@ squaredCurvature bez t
     g'' = bezier'' @v bez t
     sq_nm_g' :: r
     sq_nm_g' = squaredNorm @v g'
+{-# INLINEABLE squaredCurvature #-}
 
 -- | Signed curvature of a planar cubic Bézier curve.
 signedCurvature :: Bezier ( ℝ 2 ) -> Double -> Double
@@ -170,6 +178,7 @@ subdivide ( Bezier {..} ) t = ( Bezier p0 q1 q2 pt, Bezier pt r1 r2 p3 )
     q2 = lerp @v t q1 s
     r1 = lerp @v t s  r2
     pt = lerp @v t q2 r1
+{-# INLINEABLE subdivide #-}
 
 -- | 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 ]
@@ -188,6 +197,7 @@ ddist ( Bezier {..} ) c = [ a5, a4, a3, a2, a1, a0 ]
     !a3 = 6 * squaredNorm v'' + 4 * v' ^.^ v'''
     !a4 = 5 * v'' ^.^ v'''
     !a5 = squaredNorm v'''
+{-# INLINEABLE ddist #-}
 
 -- | Finds the closest point to a given point on a cubic Bézier curve.
 closestPoint
@@ -214,6 +224,7 @@ closestPoint pts c = pickClosest ( 0 :| 1 : roots )
             p' = bezier @v pts t'
             nm' :: r
             nm' = squaredNorm ( c --> p' :: v )
+{-# INLINEABLE closestPoint #-}
 
 -- | Drag a cubic Bézier curve to pass through a given point.
 --
@@ -236,6 +247,7 @@ drag ( Bezier {..} ) t q = Bezier { p0, p1 = p1', p2 = p2', p3 }
     p1', p2' :: p
     p1' = ( ( 1 - t ) *^ delta ) • p1
     p2' = (     t     *^ delta ) • p2
+{-# INLINEABLE drag #-}
 
 -- | Compute parameter values for the self-intersection of a planar cubic Bézier curve, if such exist.
 --
@@ -270,3 +282,4 @@ extrema ( Bezier {..} ) = solveQuadratic c b a
     a = p3 - 3 * p2 + 3 * p1 - p0
     b = 2 * ( p0 - 2 * p1 + p2 )
     c = p1 - p0
+{-# INLINEABLE extrema #-}

brush-strokes/src/Math/Bezier/Cubic/Fit.hs

@@ -59,7 +59,7 @@ import Control.Monad.Trans.State.Strict
 import Control.Monad.Trans.Class
   ( lift )
 
--- MetaBrush
+-- brush-strokes
 import qualified Math.Bezier.Cubic as Cubic
   ( Bezier(..), bezier, ddist )
 import Math.Bezier.Spline

brush-strokes/src/Math/Bezier/Quadratic.hs

@@ -48,7 +48,7 @@ import Data.Group.Generics
 import Data.Primitive.Types
   ( Prim )
 
--- MetaBrush
+-- brush-strokes
 import Math.Epsilon
   ( epsilon )
 import Math.Module
@@ -90,18 +90,22 @@ instance Show p => Show (Bezier p) where
 -- | Quadratic Bézier curve.
 bezier :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> p
 bezier ( Bezier {..} ) t = lerp @v t ( lerp @v t p0 p1 ) ( lerp @v t p1 p2 )
+{-# INLINEABLE bezier #-}
 
 -- | Derivative of a quadratic Bézier curve.
 bezier' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> v
 bezier' ( Bezier {..} ) t = Ring.fromInteger 2 *^ lerp @v t ( p0 --> p1 ) ( p1 --> p2 )
+{-# INLINEABLE bezier' #-}
 
 -- | Second derivative of a quadratic Bézier curve.
 bezier'' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> v
 bezier'' ( Bezier {..} ) = Ring.fromInteger 2 *^ ( p1 --> p0 ^+^ p1 --> p2 )
+{-# INLINEABLE bezier'' #-}
 
 -- | Curvature of a quadratic Bézier curve.
 curvature :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> r -> r
 curvature bez t = sqrt $ squaredCurvature @v bez t
+{-# INLINEABLE curvature #-}
 
 -- | Square of curvature of a quadratic Bézier curve.
 squaredCurvature :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> r -> r
@@ -117,6 +121,7 @@ squaredCurvature bez t
     g'' = bezier'' @v bez
     sq_nm_g' :: r
     sq_nm_g' = squaredNorm @v g'
+{-# INLINEABLE squaredCurvature #-}
 
 -- | Signed curvature of a planar quadratic Bézier curve.
 signedCurvature :: Bezier ( ℝ 2 ) -> Double -> Double
@@ -134,6 +139,7 @@ subdivide ( Bezier {..} ) t = ( Bezier p0 q1 pt, Bezier pt r1 p2 )
     q1 = lerp @v t p0 p1
     r1 = lerp @v t p1 p2
     pt = lerp @v t q1 r1
+{-# INLINEABLE subdivide #-}
 
 -- | 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 ]
@@ -149,6 +155,7 @@ ddist ( Bezier {..} ) c = [ a3, a2, a1, a0 ]
     !a1 = v ^.^ v'' + 2 * squaredNorm v'
     !a2 = 3 * v' ^.^ v''
     !a3 = squaredNorm v''
+{-# INLINEABLE ddist #-}
 
 -- | Finds the closest point to a given point on a quadratic Bézier curve.
 closestPoint
@@ -175,6 +182,7 @@ closestPoint pts c = pickClosest ( 0 :| 1 : roots )
             p' = bezier @v pts t'
             nm' :: r
             nm' = squaredNorm ( c --> p' :: v )
+{-# INLINEABLE closestPoint #-}
 
 -- | Interpolation of a quadratic Bézier control point, given path points and Bézier parameter.
 --
@@ -189,9 +197,11 @@ interpolate p0 p2 t q = Bezier {..}
     p1 = (   ( 0.5 * ( t - 1 ) / t ) *^ ( q --> p0 :: v )
          ^+^ ( 0.5 * t / ( t - 1 ) ) *^ ( q --> p2 :: v )
          ) • q
+{-# INLINEABLE interpolate #-}
 
 -- | Extremal values of the Bézier parameter for a quadratic Bézier curve.
 extrema :: Fractional r => Bezier r -> [ r ]
 extrema ( Bezier {..} ) = [ t ]
   where
     t = ( p0 - p1 ) / ( p0 - 2 * p1 + p2 )
+{-# INLINEABLE extrema #-}

brush-strokes/src/Math/Bezier/Spline.hs

@@ -52,7 +52,7 @@ import Control.Monad.Trans.Class
 import Control.Monad.Trans.State.Strict
   ( StateT(runStateT), modify' )
 
--- MetaBrush
+-- brush-strokes
 import qualified Math.Bezier.Cubic as Cubic
   ( Bezier(..) )
 import Math.Linear

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

@@ -24,6 +24,8 @@ module Math.Bezier.Stroke
 -- base
 import Control.Arrow
   ( first, (***) )
+import Control.Concurrent.MVar
+  ( MVar, newMVar )
 import Control.Monad
   ( unless )
 import Control.Monad.ST
@@ -39,7 +41,7 @@ import Data.Foldable
 import Data.Functor.Identity
   ( Identity(..) )
 import Data.List
-  ( nub, partition, sort )
+  ( intercalate, nub, partition, sort )
 import Data.List.NonEmpty
   ( NonEmpty )
 import qualified Data.List.NonEmpty as NE
@@ -58,6 +60,10 @@ import GHC.Generics
   ( Generic, Generic1, Generically(..) )
 import GHC.TypeNats
   ( type (-) )
+import Numeric
+  ( showFFloat )
+import System.IO.Unsafe
+  ( unsafePerformIO )
 
 -- acts
 import Data.Act
@@ -72,6 +78,8 @@ import Data.Sequence
   ( Seq(..) )
 import qualified Data.Sequence as Seq
   ( empty, index, length, reverse, singleton, zipWith )
+import Data.Tree
+  ( Tree(..) )
 
 -- deepseq
 import Control.DeepSeq
@@ -95,7 +103,11 @@ import Control.Monad.Trans.Except
 import Control.Monad.Trans.State.Strict
   ( StateT, runStateT, evalStateT, get, put )
 import Control.Monad.Trans.Writer.CPS
-  ( WriterT, execWriterT, runWriter, tell )
+  ( Writer, WriterT, execWriterT, runWriter, tell )
+
+-- tree-view
+import Data.Tree.View
+  ( showTree )
 
 -- MetaBrush
 import Math.Algebra.Dual
@@ -426,9 +438,6 @@ computeStrokeOutline rootAlgo fitParams ptParams toBrushParams brushFn spline@(
               OutlineData
                 ( TwoSided fwdData ( bimap reverseSpline Seq.reverse bwdData ) )
                 cusps
-          trace ( "bwd at t = 0.58: " ++ show ( ( snd . outlineFn fwdBwd ) $ ℝ1 0.58 ) ) ( return () )
-          trace ( "bwd at t = 0.5966724346435021: " ++ show ( ( snd . outlineFn fwdBwd ) $ ℝ1 0.5966724346435021 ) ) ( return () )
-          trace ( "bwd at t = 0.60: " ++ show ( ( snd . outlineFn fwdBwd ) $ ℝ1 0.60 ) ) ( return () )
           outlineData `deepseq` tell outlineData
           lift $ writeSTRef cachedStrokeRef ( Just outlineData )
 
@@ -581,19 +590,35 @@ outlineFunction rootAlgo ptParams toBrushParams brushFromParams = \ sp0 crv ->
                              $ runD ( brushFromParams @2 @() proxy# id )
                              $ toBrushParams params_t
 
-    ( newtDunno, newtSols ) = intervalNewtonGS InverseMidJacobian 1e-7 curvesI
+    ( newtTrees, ( newtDunno, newtSols ) ) =
+      intervalNewtonGS
+        NoPreconditioning --InverseMidJacobian
+        1e-7
+        curvesI
 
-  in --trace
+    showD :: Double -> String
-     --   ( unlines $
+    showD float = showFFloat (Just 6) float ""
-     --     [ "newtonMethod: #(definite zeroes) = " ++ show ( length newtSols  )
+    functionDataLines =
-     --     , "newtonMethod: #(unknown)         = " ++ show ( length newtDunno )
+      [ "E, dE/ds, dE/ds * dc/dt:"
-     --     , ""
+      , "{" ++
-     --     , "definite solutions:"
+        (intercalate ","
-     --     , if null newtSols then "[]" else unlines $ map show newtSols
+         [ "{" ++ showD t ++ "," ++ showD (s + fromIntegral i) ++ ",{" ++ intercalate "," vals ++ "}}"
-     --     , ""
+         | t <- [0.001, 0.01.. 0.999]
-     --     , "unknown:"
+         , i <- [0 .. (fromIntegral $ length $ curves (ℝ1 0.5)) - 1]
-     --     , if null newtDunno then "[]" else unlines $ map show newtDunno ]
+         , s <- [0.001, 0.01.. 0.999]
-     --   ) $
+         , let StrokeDatum
+                 { ee = D12 (ℝ1 e) _dEdt (T (ℝ1 dEds))
+                 , 𝛿E𝛿sdcdt = D0 (T (ℝ2 vx vy))
+                 } = (curves (ℝ1 t) `Seq.index` i) (ℝ1 s)
+               vals = [showD e, showD dEds, "{" ++ showD vx ++ "," ++ showD vy ++ "}"]
+         ]
+        ) ++ "}"
+      ]
+    newtTreeLines = map (showTree . showIntervalNewtonTree 1) newtTrees
+
+    logContents = unlines $ functionDataLines ++ newtTreeLines
+
+  in logToFile cuspFindingMVar logContents `seq`
      OutlineInfo
        { outlineFn = fwdBwd
        , outlineDefiniteCusps  = map ( cuspCoords curves ) newtSols
@@ -601,6 +626,11 @@ outlineFunction rootAlgo ptParams toBrushParams brushFromParams = \ sp0 crv ->
        }
 {-# INLINEABLE outlineFunction #-}
 
+cuspFindingMVar :: MVar FilePath
+cuspFindingMVar = unsafePerformIO
+                $ newMVar "logs/cusps.txt"
+{-# NOINLINE cuspFindingMVar #-}
+
 pathAndUsedParams :: forall k i arr crvData ptData usedParams
                   .  ( HasType ( ℝ 2 ) ptData
                      , HasBézier k i
@@ -962,6 +992,11 @@ data RootSolvingAlgorithm
   -- | Use the modified Halley M2 method.
   | HalleyM2 { maxIters :: Word, precision :: Int }
 
+rootSolvingMVar :: MVar FilePath
+rootSolvingMVar = unsafePerformIO
+                $ newMVar "logs/envelope.txt"
+{-# NOINLINE rootSolvingMVar #-}
+
 -- | 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 :: RootSolvingAlgorithm
@@ -986,17 +1021,17 @@ solveEnvelopeEquations rootAlgo _t path_t path'_t ( fwdOffset, bwdOffset ) strok
     --        , "  bwdOffset: " ++ show bwdOffset ]
     --      ) True
 
-    fwdSol            = findSolFrom fwdOffset
+    fwdSol            = findSolFrom "fwd" fwdOffset
-    ( bwdPt, bwdTgt ) = findSolFrom bwdOffset
+    ( bwdPt, bwdTgt ) = findSolFrom "bwd" bwdOffset
 
-    findSolFrom :: Offset -> ( ℝ 2, T ( ℝ 2 ) )
+    findSolFrom :: String -> Offset -> ( ℝ 2, T ( ℝ 2 ) )
-    findSolFrom ( Offset { offsetIndex = i00, offsetParameter = s00, offset = off } )
+    findSolFrom desc ( Offset { offsetIndex = i00, offsetParameter = s00, offset = off } )
       = go ( fromIntegral i00 + fromMaybe 0.5 s00 )
       where
 
         go :: Double -> ( ℝ 2, T ( ℝ 2 ) )
         go is0 =
-          case sol strokeData is0 of
+          case sol desc strokeData is0 of
             ( goodSoln, pt, tgt )
               | goodSoln && plausibleTangent tgt
               -> ( pt, tgt )
@@ -1006,59 +1041,100 @@ solveEnvelopeEquations rootAlgo _t path_t path'_t ( fwdOffset, bwdOffset ) strok
         plausibleTangent :: T ( ℝ 2 ) -> Bool
         plausibleTangent tgt = path'_t ^.^ tgt > 0
 
-    sol :: Seq ( ℝ 1 -> StrokeDatum 2 () ) -> Double -> ( Bool, ℝ 2, T ( ℝ 2 ) )
+    sol :: String -> Seq ( ℝ 1 -> StrokeDatum 2 () ) -> Double -> ( Bool, ℝ 2, T ( ℝ 2 ) )
-    sol f is0 =
+    sol desc f is0 =
-      let ( good, is ) =
+      let (solRes, solSteps) = runSolveMethod ( eqn f ) is0
-            case runSolveMethod ( eqn f ) is0 of
+          ( good, is ) =
+            case solRes of
               Nothing  -> ( False, is0 )
               Just is1 -> ( True , is1 )
           ( ds, dcdt ) = finish f is
-      in ( good, ds, dcdt )
+          showD :: Double -> String
+          showD float = showFFloat (Just 6) float ""
+          (x_min, x_max) = domain
+          logContents = unlines
+              [ "method = " <> methodName
+              , desc
+              , "t = " ++ show _t
+              , "start (n,s) = " <> show (fromDomain is0)
+              , "steps = {"
+              , unlines $ map (\ (x, y, y') -> show (x, (y, y'))) solSteps
+              , "},"
+              , "E:"
+              , "{" ++
+                 (intercalate ","
+                   [ "{" ++ showD x ++ "," ++ showD y ++ "}"
+                   | x <- [x_min, x_min + 0.01 .. x_max]
+                   , let (# y, _y' #) = eqn f x
+                   ]
+                 ) ++ "}"
+              , "𝛿E𝛿s:"
+              , "{" ++
+                (intercalate ","
+                 [ "{" ++ showD x ++ "," ++ showD y' ++ "}"
+                 | x <- [x_min, x_min + 0.01 .. x_max]
+                 , let (# _y, y' #) = eqn f x
+                 ]
+                ) ++ "}"
+              , "𝛿E𝛿sdcdt (x):"
+              , "{" ++
+                (intercalate ","
+                 [ "{" ++ showD x ++ "," ++ showD 𝛿E𝛿sdcdt_x ++ "}"
+                 | x <- [x_min, x_min + 0.01 .. x_max]
+                 , let StrokeDatum { 𝛿E𝛿sdcdt = D0 (T (ℝ2 𝛿E𝛿sdcdt_x _)) } = evalStrokeDatum f x
+                 ]
+                )
+                ++ "}"
+              , "𝛿E𝛿sdcdt (y):"
+              , "{" ++
+                (intercalate ","
+                 [ "{" ++ showD x ++ "," ++ showD 𝛿E𝛿sdcdt_y ++ "}"
+                 | x <- [x_min, x_min + 0.01 .. x_max]
+                 , let StrokeDatum { 𝛿E𝛿sdcdt = D0 (T (ℝ2 _ 𝛿E𝛿sdcdt_y)) } = evalStrokeDatum f x
+                 ]
+                )
+                ++ "}"
+              ]
 
-    runSolveMethod = case rootAlgo of
+      in logToFile rootSolvingMVar logContents `seq`
+         ( good, ds, dcdt )
+
+    (runSolveMethod, methodName) = case rootAlgo of
       HalleyM2 { maxIters, precision } ->
-        halleyM2 maxIters precision
+        (halleyM2 maxIters precision, "HalleyM2")
       NewtonRaphson { maxIters, precision } ->
-        newtonRaphson maxIters precision domain
+        (newtonRaphson maxIters precision domain, "NewtonRaphson")
 
     finish :: Seq ( ℝ 1 -> StrokeDatum 2 () ) -> Double -> ( ℝ 2, T ( ℝ 2 ) )
-    finish f is =
+    finish fs is =
       let (i, s) = fromDomain is in
-      case ( f `Seq.index` i ) ( ℝ1 s ) of -- TODO: a bit redundant to have to compute this again...
+      case evalStrokeDatum fs is of -- TODO: a bit redundant to have to compute this again...
         StrokeDatum
           { dstroke
-          , ee = D12 ( ℝ1 _ee ) ( T ( ℝ1 _𝛿E𝛿t ) ) ( T ( ℝ1 ee_s ) )
+          , ee = D12 ( ℝ1 _ee ) _ ( T ( ℝ1 𝛿E𝛿s ) )
           , 𝛿E𝛿sdcdt = D0 𝛿E𝛿sdcdt
           } ->
           -- The total derivative dc/dt is computed by dividing by ∂E/∂s,
           -- so check it isn't zero first. This corresponds to cusps in the envelope.
           let dcdt
-                | abs ee_s < epsilon
+                | abs 𝛿E𝛿s < epsilon
-                , let s' = if s >= 0.5 then s - 1e-9 else s + 1e-9
+                , let s' = if s >= 0.5 then s - 1e-6 else s + 1e-6
-                = case ( f `Seq.index` i ) ( ℝ1 s' ) of
+                = case ( fs `Seq.index` i ) ( ℝ1 s' ) of
-                   StrokeDatum { ee = D12 _ _ ( T ( ℝ1 ee_s' ) ), 𝛿E𝛿sdcdt = D0 𝛿E𝛿sdcdt' }
+                   StrokeDatum { ee = D12 _ _ ( T ( ℝ1 𝛿E𝛿s' ) ), 𝛿E𝛿sdcdt = D0 𝛿E𝛿sdcdt' }
-                     -> recip ee_s' *^ 𝛿E𝛿sdcdt'
+                     -> recip 𝛿E𝛿s' *^ 𝛿E𝛿sdcdt'
                 | otherwise
-                = recip ee_s *^ 𝛿E𝛿sdcdt
+                = recip 𝛿E𝛿s *^ 𝛿E𝛿sdcdt
-          in --trace
+          in ( value @Double @2 @( ℝ 2 ) dstroke, dcdt )
-             --  ( unlines
+
-             --    [ "solveEnvelopeEquations"
+    evalStrokeDatum :: Seq ( ℝ 1 -> StrokeDatum 2 () ) -> ( Double -> StrokeDatum 2 () )
-             --    , "      t = " ++ show _t
+    evalStrokeDatum fs is =
-             --    , "      s = " ++ show s
+      let (i, s) = fromDomain is
-             --    , "      c = " ++ show dstroke
+      in ( fs `Seq.index` i ) ( ℝ1 $ max 1e-6 $ min (1 - 1e-6) $ s )
-             --    , "      E = " ++ show _ee
-             --    , "  ∂E/∂t = " ++ show _𝛿E𝛿t
-             --    , "  ∂E/∂s = " ++ show ee_s
-             --    , "  dc/dt = " ++ show dcdt
-             --    ] )
-              ( value @Double @2 @( ℝ 2 ) dstroke, dcdt )
 
     eqn :: Seq ( ℝ 1 -> StrokeDatum 2 () ) -> ( Double -> (# Double, Double #) )
     eqn fs is =
-      let (i, s) = fromDomain is
+      case evalStrokeDatum fs is of
-      in case ( fs `Seq.index` i ) ( ℝ1 s ) of
+        StrokeDatum { ee = D12 ee _ ee_s } -> coerce (# ee, ee_s #)
-           StrokeDatum { ee = D12 ee _ ee_s } ->
-             coerce (# ee, ee_s #)
 
     n :: Int
     n = Seq.length strokeData
@@ -1161,24 +1237,26 @@ gaussSeidel
   ( 𝕀 ( ℝ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 ) )
-  ( 𝕀 ( ℝ1 t0_lo ) ( ℝ1 t0_hi ), 𝕀 ( ℝ1 s0_lo ) ( ℝ1 s0_hi ) )
+  ( 𝕀 ( ℝ1 x1_lo ) ( ℝ1 x1_hi ), 𝕀 ( ℝ1 x2_lo ) ( ℝ1 x2_hi ) )
     = let !a11 = 𝕀 a11_lo a11_hi
           !a12 = 𝕀 a12_lo a12_hi
           !a21 = 𝕀 a21_lo a21_hi
           !a22 = 𝕀 a22_lo a22_hi
           !b1  = 𝕀  b1_lo  b1_hi
           !b2  = 𝕀  b2_lo  b2_hi
-          !t0  = 𝕀  t0_lo  t0_hi
+          !x1  = 𝕀  x1_lo  x1_hi
-          !s0  = 𝕀  s0_lo  s0_hi
+          !x2  = 𝕀  x2_lo  x2_hi
       in nub $ do
 
-        t' <- ( b1 - a12 * s0 ) `extendedDivide` a11
+        -- x1' = (b1 - a12 * x2) / a11
-        ( t@( 𝕀 t_lo t_hi ), sub_t ) <- t' `intersect` t0
+        x1'_0 <- ( b1 - a12 * x2 ) `extendedDivide` a11
-        s' <- ( b2 - a21 * t  ) `extendedDivide` a22
+        ( x1'@( 𝕀 x1'_lo x1'_hi ), sub_x1 ) <- x1'_0 `intersect` x1
-        ( 𝕀 s_lo s_hi, sub_s ) <- s' `intersect` s0
+        -- x2' = (b2 - a21 * x1') / a22
+        x2'_0 <- ( b2 - a21 * x1' ) `extendedDivide` a22
+        ( 𝕀 x2'_lo x2'_hi, sub_x2 ) <- x2'_0 `intersect` x2
 
-        return ( ( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ), 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
+        return ( ( 𝕀 ( ℝ1 x1'_lo ) ( ℝ1 x1'_hi ), 𝕀 ( ℝ1 x2'_lo ) ( ℝ1 x2'_hi ) )
-               , sub_t && sub_s )
+               , sub_x1 && sub_x2 )
 
 intersect :: 𝕀 Double -> 𝕀 Double -> [ ( 𝕀 Double, Bool ) ]
 intersect ( 𝕀 lo1 hi1 ) ( 𝕀 lo2 hi2 )
@@ -1227,46 +1305,87 @@ data Preconditioner
   | InverseMidJacobian
   deriving stock Show
 
+-- | A tree recording the steps taken with the interval Newton method.
+data IntervalNewtonTree d
+  = IntervalNewtonLeaf (IntervalNewtonLeaf d)
+  | IntervalNewtonStep (IntervalNewtonStep d) [(Double, IntervalNewtonTree d)]
+
+data IntervalNewtonStep d
+  = IntervalNewtonContraction [d] [d]
+  | IntervalNewtonBisection (String, Double)
+
+instance Show d => Show (IntervalNewtonStep d) where
+  showsPrec _ ( IntervalNewtonContraction v w )
+    = showString "N "
+    . showsPrec 0 v
+    . showString " "
+    . showsPrec 0 w
+  showsPrec _ ( IntervalNewtonBisection (coord, w) )
+    = showString "B "
+    . showParen True
+      ( showString coord
+      . showString " = "
+      . showsPrec 0 w
+      )
+
+data IntervalNewtonLeaf d
+  = TooSmall d
+  | NoSolution d
+  deriving stock Show
+
+showIntervalNewtonTree :: Show d => Double -> IntervalNewtonTree d -> Tree String
+showIntervalNewtonTree area (IntervalNewtonLeaf l) = Node (showArea area ++ " " ++ show l) []
+showIntervalNewtonTree area (IntervalNewtonStep s ts)
+  = Node (showArea area ++ " " ++ show s) $ map (\ (d,t) -> showIntervalNewtonTree d t) ts
+
+showArea :: Double -> [Char]
+showArea area = "(area " ++ showFFloat (Just 6) area "" ++ ")"
+
 -- | Interval Newton method with Gauss–Seidel step for inversion
 -- of the interval Jacobian.
 --
--- Returns @(dunno, sols)@ where @sols@ are boxes that contain a unique solution
+-- Returns @(tree, (dunno, sols))@ where @tree@ is the full tree (useful for debugging),
--- (and to which Newton's method will converge starting from anywhere inside
+-- @sols@ are boxes that contain a unique solution (and to which Newton's method
--- the box), and @dunno@ are small boxes which might or might not
+-- will converge starting from anywhere inside the box), and @dunno@ are small
--- contain solutions.
+-- boxes which might or might not contain solutions.
 intervalNewtonGS :: Preconditioner
                  -> Double
                  -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-                 -> ( [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ], [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] )
+                 -> ( [ IntervalNewtonTree ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ], ( [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ], [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] ) )
 intervalNewtonGS precondMethod minWidth eqs =
-  go (0,0)
+  first concat $ runWriter $
-    [ ( 𝕀 ( ℝ1 0 ) ( ℝ1 1 ), i, 𝕀 ( ℝ1 0 ) ( ℝ1 1 ) )
+    traverse go
-    | i <- [ 0 .. length ( eqs ( 𝕀 ( ℝ1 0 ) ( ℝ1 1 ) ) ) - 1 ]
+      [ ( 𝕀 ( ℝ1 zero ) ( ℝ1 one ), i, 𝕀 ( ℝ1 zero ) ( ℝ1 one ) )
-    ]
+      | i <- [ 0 .. length ( eqs ( 𝕀 ( ℝ1 zero ) ( ℝ1 one ) ) ) - 1 ]
-    []
+      ]
-    []
 
   where
+    zero, one :: Double
+    zero = 1e-6
+    one = 1 - zero
+    {-# INLINE zero #-}
+    {-# INLINE one #-}
 
-    go :: ( Int, Int ) -- step counts (for debugging)
+    recur f cands = do
-       -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] -- boxes to work on
+      rest <- traverse ( \ c -> do { trees <- go c; return [ (boxArea c, tree) | tree <- trees ] } ) cands
-       -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] -- too small: don't shrink further
+      return [ f $ concat rest ]
-       -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] -- found solutions
+
-       -> ( [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ], [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] )
+    boxArea :: ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) -> Double
-    go ( bis, newt ) [] giveUp sols =
+    boxArea ( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ), _, 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
-      trace ( unlines [ "intervalNewtonGS done"
+      = abs ( t_hi - t_lo ) * abs ( s_hi - s_lo )
-                      , "   #bisections: " ++ show bis
+
-                      , "  #newtonSteps: " ++ show newt
+    go :: ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) -- box to work on
-                      , "         #sols: " ++ show ( length sols )
+       -> Writer ( [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ], [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] )
-                      , "        #dunno: " ++ show ( length giveUp ) ] )
+                 [ IntervalNewtonTree ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ]
-        ( giveUp, sols )
+    go cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) )
-    go ( bis, newt ) ( cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) )
+            , i
-                   , i
+            , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
-                   , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
+            )
-                   ) : cands ) giveUp sols
       -- Box is small: stop processing it.
       | t_hi - t_lo < minWidth && s_hi - s_lo < minWidth
-      = go ( bis, newt ) cands ( cand : giveUp ) sols
+      = do let res = TooSmall cand
+           tell ( [ cand ], [] )
+           return [ IntervalNewtonLeaf res ]
 
       | StrokeDatum { ee = D22 ee _ _ _ _ _
                     , 𝛿E𝛿sdcdt = D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) }
@@ -1274,10 +1393,12 @@ intervalNewtonGS precondMethod minWidth eqs =
 
       , StrokeDatum { 𝛿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
-      = if | inf ee < ℝ1 0
+      = if | -- Check the envelope equation interval contains zero.
-           , sup ee > ℝ1 0
+             inf ee <= ℝ1 0
-           , cmpℝ2 (<) ( inf f ) ( ℝ2 0 0 )
+           , sup ee >= ℝ1 0
-           , cmpℝ2 (>) ( sup f ) ( ℝ2 0 0 )
+             -- Check the 𝛿E𝛿sdcdt interval contains zero.
+           , cmpℝ2 (<=) ( inf f ) ( ℝ2 0 0 )
+           , cmpℝ2 (>=) ( sup f ) ( ℝ2 0 0 )
            -> let -- Interval Newton method: take one Gauss–Seidel step
                   -- for the equation f'(X) v = - f(x_mid).
                   !( a, b ) = precondition precondMethod
@@ -1296,22 +1417,25 @@ intervalNewtonGS precondMethod minWidth eqs =
                    -- 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 go ( bis, newt + 1 ) ( todo ++ cands ) giveUp ( done ++ sols )
+                   in do tell ([], done)
+                         recur (IntervalNewtonStep (IntervalNewtonContraction done todo)) todo
                  else
                     -- Gauss–Seidel failed to shrink the boxes.
                     -- Bisect along the widest dimension instead.
-                   let bisGuesses
+                   let (bisGuesses, whatBis)
                          | t_hi - t_lo > s_hi - s_lo
-                         = [ ( 𝕀 ( ℝ1 t_lo  ) ( ℝ1 t_mid ), i, s )
+                         = ( [ ( 𝕀 ( ℝ1 t_lo  ) ( ℝ1 t_mid ), i, s )
-                           , ( 𝕀 ( ℝ1 t_mid ) ( ℝ1 t_hi  ), 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_lo  ) ( ℝ1 s_mid ) )
-                           , ( t, i, 𝕀 ( ℝ1 s_mid ) ( ℝ1 s_hi  ) ) ]
+                             , ( t, i, 𝕀 ( ℝ1 s_mid ) ( ℝ1 s_hi  ) ) ]
-                   in go ( bis + 1, newt ) ( bisGuesses ++ cands ) giveUp sols
+                           , ("s", s_mid) )
+                   in recur (IntervalNewtonStep (IntervalNewtonBisection whatBis)) bisGuesses
 
            -- Box doesn't contain a solution: discard it.
            | otherwise
-           -> go ( bis, newt ) cands giveUp sols
+           -> return [ IntervalNewtonLeaf $ NoSolution cand ]
       where
         t_mid = 0.5 * ( t_lo + t_hi )
         s_mid = 0.5 * ( s_lo + s_hi )

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

@@ -22,7 +22,7 @@ import GHC.TypeNats
 import Data.Act
   ( Torsor ((-->)) )
 
--- MetaBrush
+-- brush-strokes
 import Math.Algebra.Dual
 import qualified Math.Bezier.Cubic     as Cubic
 import qualified Math.Bezier.Quadratic as Quadratic
@@ -63,7 +63,6 @@ data StrokeDatum k i
    --
    -- denotes a total derivative.
   , 𝛿E𝛿sdcdt :: D ( k - 2 ) ( I i ( ℝ 2 ) ) ( T ( I i ( ℝ 2 ) ) )
-
   }
 
 deriving stock instance Show ( StrokeDatum 2 () )

brush-strokes/src/Math/Differentiable.hs

@@ -11,7 +11,7 @@ import Data.Kind
 import GHC.TypeNats
   ( Nat )
 
--- MetaBrush
+-- brush-strokes
 import Math.Algebra.Dual
   ( D, HasChainRule )
 import Math.Interval

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

@@ -28,7 +28,7 @@ import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval
   ( Interval(..) )
 
 #ifdef USE_FMA
--- MetaBrush
+-- brush-strokes
 import Math.Interval.FMA
   ( addI, subI, prodI, divI, posPowI )
 
@@ -38,7 +38,7 @@ import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval
   ( powInt )
 #endif
 
--- MetaBrush
+-- brush-strokes
 import Math.Linear
   ( T(..)
   , RepDim, RepresentableQ(..), Representable(..)

brush-strokes/src/Math/Linear.hs

@@ -46,7 +46,7 @@ import Data.Group
 import Data.Group.Generics
   ( )
 
--- MetaBrush
+-- brush-strokes
 import Math.Linear.Internal
 
 --------------------------------------------------------------------------------

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

@@ -14,7 +14,7 @@ import qualified Eigen.Matrix as Eigen
 import qualified Eigen.Solver.LA as Eigen
   ( Decomposition(..), solve )
 
--- MetaBrush
+-- brush-strokes
 import Math.Linear
   ( Mat22(..), ℝ(..), T(..) )
 

brush-strokes/src/Math/Module.hs

@@ -40,7 +40,7 @@ import Data.Act
 import Data.Group
   ( Group(..) )
 
--- MetaBrush
+-- brush-strokes
 import Math.Epsilon
   ( epsilon )
 import Math.Linear

brush-strokes/src/Math/Module/Internal.hs

@@ -10,7 +10,7 @@ import Data.Coerce
 import Data.Monoid
   ( Sum(..) )
 
--- MetaBrush
+-- brush-strokes
 import Math.Ring
   ( Ring )
 import qualified Math.Ring as Ring

brush-strokes/src/Math/Monomial.hs

@@ -36,7 +36,7 @@ import Unsafe.Coerce
 import Language.Haskell.TH
   ( CodeQ )
 
--- MetaBrush
+-- brush-strokes
 import Math.Linear
   ( Vec(..), Fin(..) )
 import TH.Utils

brush-strokes/src/Math/Orientation.hs

@@ -28,7 +28,7 @@ import Data.Generics.Product.Typed
 import Data.Generics.Internal.VL
   ( view )
 
--- MetaBrush
+-- brush-strokes
 import Math.Epsilon
   ( nearZero )
 import Math.Module

brush-strokes/src/Math/Roots.hs

@@ -46,7 +46,7 @@ import Data.Primitive.PrimArray
 import Data.Primitive.Types
   ( Prim )
 
--- MetaBrush
+-- brush-strokes
 import Math.Epsilon
   ( epsilon, nearZero )
 
@@ -298,7 +298,7 @@ derivative p = do
 -- https://github.com/boostorg/math/blob/0dc6a70caa6bbec2b6ae25eede36c430f0ccae13/include/boost/math/tools/roots.hpp#L217
 {-# SPECIALISE
   newtonRaphson
-    :: Word -> Int -> ( Double, Double ) -> ( Double -> (# Double, Double #) ) -> Double -> Maybe Double
+    :: Word -> Int -> ( Double, Double ) -> ( Double -> (# Double, Double #) ) -> Double -> (Maybe Double, [(Double,Double,Double)])
   #-}
 {-# INLINEABLE newtonRaphson #-}
 newtonRaphson :: ( RealFloat r, Show r )
@@ -307,7 +307,7 @@ newtonRaphson :: ( RealFloat r, Show r )
               -> ( r, r )           -- ^ @(min_x, max_x)@.
               -> ( r -> (# r, r #) ) -- ^ function and its derivative
               -> r                  -- ^ initial guess
-              -> Maybe r
+              -> (Maybe r, [(r,r,r)])
 newtonRaphson maxIters digits ( min_x, max_x ) f x0 =
   doNewtonRaphson f maxIters factor min_x max_x 0 0 0 x0 maxRealFloat maxRealFloat
   where
@@ -322,30 +322,31 @@ doNewtonRaphson :: ( Fractional r, Ord r, Show r )
                 -> r
                 -> r
                 -> r -> r
-                -> Maybe r
+                -> (Maybe r, [(r,r,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
+  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 =
+  go prev_acc min_x max_x min_f_x max_f_x f_x_prev !iters !x !δ1 !δ2 =
-    case f x of
+    let (# f_x, f'_x #) = f x
-      (# f_x, f'_x #)
+        acc = ((x, f_x, f'_x) : prev_acc)
+    in if
         | f_x == 0
-        -> Just x
+        -> (Just x, acc)
         | ( 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
+            -> (Just x, acc)
             | iters >= maxIters
-            -> Nothing
+            -> (Nothing, acc)
             | ( 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
+               then (Nothing, acc)
                else
-                 go min_x max_x min_f_x max_f_x f_x ( iters + 1 ) new_x δ δ1
+                 go acc min_x max_x min_f_x max_f_x f_x ( iters + 1 ) new_x δ δ1
 
 newtonRaphsonStep :: ( Fractional r, Ord r, Show r )
                   => ( r -> (# r, r #) )
@@ -444,9 +445,14 @@ halleyStep (# x, f, f', f'' #) =
 -- by A. Cordero, H. Ramos & J.R. Torregrosa, J Math Chem 58, 751–774 (2020).
 --
 -- @https://doi.org/10.1007/s10910-020-01108-3@
-halleyM2Step :: Fractional a => (# a, (# a, a #) #) -> (# a, (# a, a #) #) -> a
+halleyM2Step :: RealFloat a => (# a, (# a, a #) #) -> (# a, (# a, a #) #) -> a
-halleyM2Step (# x_nm1, (# f_nm1, f'_nm1 #) #) (# x_n, (# f_n, f'_n #) #) =
+halleyM2Step (# x_nm1, (# f_nm1, f'_nm1 #) #) (# x_n, (# f_n, f'_n #) #)
-  num / denom
+  | nearZero num && nearZero denom
+  = 0.1 * signum num * signum denom
+  | nearZero num
+  = num
+  | otherwise
+  = num / denom
   where
     u = f_n * f_nm1 * (f_n - f_nm1)
     dx = x_n - x_nm1
@@ -457,7 +463,7 @@ halleyM2Step (# x_nm1, (# f_nm1, f'_nm1 #) #) (# x_n, (# f_n, f'_n #) #) =
 
 {-# SPECIALISE
   halleyM2
-    :: Word -> Int -> ( Double -> (# Double, Double #) ) -> Double -> Maybe Double
+    :: Word -> Int -> ( Double -> (# Double, Double #) ) -> Double -> (Maybe Double, [(Double, Double, Double)])
   #-}
 {-# INLINEABLE halleyM2 #-}
 halleyM2 :: ( RealFloat r, Show r )
@@ -465,16 +471,17 @@ halleyM2 :: ( RealFloat r, Show r )
          -> Int               -- ^ desired digits of precision
          -> ( r -> (# r, r #) ) -- ^ function and its derivative
          -> r                 -- ^ initial guess
-         -> Maybe r
+         -> (Maybe r, [(r,r,r)])
 halleyM2 maxIters digits f x0 =
   let y0 = (# x0, f x0 #)
-  in go 0 y0 y0
+  in go [] 0 y0 y0
   where
     !factor = encodeFloat 1 ( 1 - digits )
-    go i y_nm1 y_n@(# x_n, _ #) =
+    go prev_acc i y_nm1 y_n@(# x_n, (# f_x_n, f'_x_n #) #) =
-      let x_np1 = halleyM2Step y_nm1 y_n
+      let acc = (x_n, f_x_n, f'_x_n) : prev_acc
+          x_np1 = halleyM2Step y_nm1 y_n
       in if | i >= maxIters
             || abs ( x_np1 - x_n ) < abs ( x_n * factor )
-            -> Just x_np1
+            -> (Just x_np1, acc)
             | otherwise
-            -> go (i+1) y_n (# x_n, f x_np1 #)
+            -> go acc (i+1) y_n (# x_np1, f x_np1 #)

brush-strokes/src/TH/Utils.hs

@@ -6,7 +6,7 @@ module TH.Utils where
 import Language.Haskell.TH
   ( CodeQ )
 
--- MetaBrush
+-- brush-strokes
 import Math.Ring ( Ring )
 import qualified Math.Ring as Ring
 

cabal.project

@@ -1,9 +1,8 @@
-packages: brush-strokes
+packages: ., brush-strokes
-        , .
 
 constraints:
   acts -finitary,
-  brush-strokes +use-fma,
+--  brush-strokes +use-fma,
   rounded-hw -pure-hs -c99 -avx512 +ghc-prim -x87-long-double,
   text -simdutf
 -- text +simdutf causes the "digit" package to fail to build with undefined symbol linker errors

hie.yaml

@@ -1,7 +1,7 @@
 cradle:
   cabal:
-    - path: "./src/splines"
+    - path: "./brush-strokes/src"
-      component: "lib:splines"
+      component: "lib:brush-strokes"
     - path: "./src/metabrushes"
       component: "lib:metabrushes"
     - path: "./src/convert"

src/app/MetaBrush/Application.hs

@@ -209,7 +209,9 @@ runApplication application = do
                               , maxIters   = 5   -- 100
                               }
   rootsAlgoTVar        <- STM.newTVarIO @RootSolvingAlgorithm $
-                            HalleyM2 { maxIters = 20, precision = 8 }
+                            --HalleyM2
+                            NewtonRaphson
+                              { maxIters = 20, precision = 8 }
 
   -- Put all these stateful variables in a record for conciseness.
   let