Restructure project & update bounds

MetaBrush.cabal

@@ -29,34 +29,27 @@ flag asserts
   default: False
   manual:  True
 
-flag use-fma
-  description: Use fused-muliply add instructions to implement interval arithmetic.
-  default: False
-  manual:  True
-
 common common
 
     build-depends:
-        base
+        brush-strokes
-           >= 4.17    && < 4.19
+          ^>= 0.1.0.0
+      , base
+           >= 4.17    && < 5
       , acts
           ^>= 0.3.1.0
       , code-page
           ^>= 0.2.1
       , containers
-           >= 0.6.0.1 && < 0.7
+           >= 0.6.0.1 && < 0.8
       , deepseq
-           >= 1.4.4.0 && < 1.5
+           >= 1.4.4.0 && < 1.6
       , generic-lens
            >= 2.2     && < 2.3
       , groups
           ^>= 0.5.3
-      , groups-generic
-          ^>= 0.3.1.0
       , primitive
-          ^>= 0.7.1.0
+          ^>= 0.9.0.0
-      , rounded-hw
-          ^>= 0.3
       , transformers
            >= 0.5.6.2 && < 0.7
 
@@ -94,14 +87,11 @@ common common
         UnboxedTuples
         ViewPatterns
 
+    if impl(ghc >= 9.8)
+      default-extensions:
+        TypeAbstractions
+
     ghc-options:
-        -O2
-        -fexpose-all-unfoldings
-      --  -funfolding-use-threshold=1000
-        -fspecialise-aggressively
-        -flate-dmd-anal
-        -fmax-worker-args=200
-        -optc-O3
         -Wall
         -Wcompat
         -fwarn-missing-local-signatures
@@ -114,14 +104,6 @@ common common
       cpp-options:
           -DASSERTS
 
-    if flag(use-fma)
-      cpp-options:
-          -DUSE_FMA
-      ghc-options:
-          -mfma
-      if impl(ghc < 9.7)
-        buildable: False
-
     autogen-modules:
         Paths_MetaBrush
 
@@ -131,11 +113,10 @@ common common
 common extras
 
     build-depends:
-
         directory
            >= 1.3.4.0    && < 1.4
       , filepath
-          ^>= 1.4.2.1
+           >= 1.4.2.1    && < 1.6
       , hashable
            >= 1.3.0.0    && < 1.5
       , lens
@@ -147,9 +128,9 @@ common extras
       , stm
           ^>= 2.5.0.0
       , text
-           >= 1.2.3.1    && < 2.1
+          ^>= 2.1.1
       , unordered-containers
-           >= 0.2.11     && < 0.2.20
+           >= 0.2.11     && < 0.3
       , waargonaut
           ^>= 0.8.0.2
 
@@ -175,64 +156,6 @@ common gtk
       , haskell-gi-base
            >= 0.26       && < 0.27
 
-
-library splines
-
-    import:
-        common
-
-    hs-source-dirs:
-        src/splines
-
-    default-language:
-        Haskell2010
-
-    exposed-modules:
-        Math.Algebra.Dual
-      , Math.Bezier.Cubic
-      , Math.Bezier.Cubic.Fit
-      , Math.Bezier.Quadratic
-      , Math.Bezier.Spline
-      , Math.Bezier.Stroke
-      , Math.Bezier.Stroke.EnvelopeEquation
-      , Math.Differentiable
-      , Math.Epsilon
-      , Math.Interval
-      , Math.Linear
-      , Math.Linear.Solve
-      , Math.Module
-      , Math.Monomial
-      , Math.Orientation
-      , Math.Ring
-      , Math.Roots
-      , Debug.Utils
-
-    other-modules:
-        Math.Algebra.Dual.Internal
-      , Math.Interval.Internal
-      , Math.Linear.Internal
-      , Math.Module.Internal
-      , TH.Utils
-
-    if flag(use-fma)
-      other-modules:
-        Math.Interval.FMA
-
-    build-depends:
-        bifunctors
-           >= 5.5.4      && < 5.6
-      , eigen
-          ^>= 3.3.7.0
-      , parallel
-          ^>= 3.2.2.0
-      , prim-instances
-          ^>= 0.2
-      , vector
-           >= 0.12.1.2   && < 0.14
-      , template-haskell
-           >= 2.18       && < 2.20
-
-
 library metabrushes
 
     import:
@@ -259,8 +182,7 @@ library metabrushes
       , MetaBrush.Util
 
     build-depends:
-        splines
+        atomic-file-ops
-      , atomic-file-ops
           ^>= 0.3.0.0
       , bytestring
            >= 0.10.10.0  && < 0.12
@@ -282,9 +204,6 @@ executable cusps
     other-modules:
         Math.Interval.Abstract
 
-    build-depends:
-        splines
-
 
 executable convert-metafont
 
@@ -304,7 +223,6 @@ executable convert-metafont
         MetaBrush.MetaFont.Convert
 
     build-depends:
-        splines,
         metabrushes,
         diagrams-contrib,
         diagrams-lib,
@@ -358,7 +276,6 @@ executable MetaBrush
         -threaded -rtsopts
 
     build-depends:
-        splines
+        metabrushes
-      , metabrushes
       , tardis
            >= 0.4.2.0    && < 0.5

brush-strokes/brush-strokes.cabal

@@ -0,0 +1,152 @@
+cabal-version:  3.0
+name:           brush-strokes
+version:        0.1.0.0
+synopsis:       Stroking brush paths
+category:       Calligraphy, Geometry, Graphics
+license:        BSD-3-Clause
+homepage:       https://gitlab.com/sheaf/MetaBrush
+build-type:     Simple
+description:
+  Computing brush strokes using Bézier curves.
+
+flag use-fma
+  description: Use fused-muliply add instructions to implement interval arithmetic.
+  default: True
+  manual:  True
+
+common common
+
+    build-depends:
+        base
+           >= 4.17    && < 4.20
+      , acts
+          ^>= 0.3.1.0
+      , code-page
+          ^>= 0.2.1
+      , containers
+           >= 0.6.0.1 && < 0.8
+      , deepseq
+           >= 1.4.4.0 && < 1.6
+      , generic-lens
+           >= 2.2     && < 2.3
+      , groups
+          ^>= 0.5.3
+      , groups-generic
+          ^>= 0.3.1.0
+      , primitive
+          ^>= 0.9.0.0
+      , rounded-hw
+          ^>= 0.4
+      , transformers
+           >= 0.5.6.2 && < 0.7
+
+    default-extensions:
+        BangPatterns
+        BlockArguments
+        ConstraintKinds
+        DataKinds
+        DeriveAnyClass
+        DeriveTraversable
+        DeriveGeneric
+        DerivingVia
+        FlexibleContexts
+        FlexibleInstances
+        FunctionalDependencies
+        GADTs
+        GeneralisedNewtypeDeriving
+        InstanceSigs
+        LambdaCase
+        LexicalNegation
+        MagicHash
+        MultiWayIf
+        NamedFieldPuns
+        NoStarIsType
+        PatternSynonyms
+        RankNTypes
+        RecordWildCards
+        RoleAnnotations
+        StandaloneDeriving
+        StandaloneKindSignatures
+        TupleSections
+        TypeApplications
+        TypeFamilies
+        TypeOperators
+        UnboxedTuples
+        ViewPatterns
+
+    ghc-options:
+        -O2
+        -fexpose-all-unfoldings
+      --  -funfolding-use-threshold=1000
+        -fspecialise-aggressively
+        -flate-dmd-anal
+        -fmax-worker-args=200
+        -optc-O3
+        -Wall
+        -Wcompat
+        -fwarn-missing-local-signatures
+        -fwarn-incomplete-patterns
+        -fwarn-incomplete-uni-patterns
+        -fwarn-missing-deriving-strategies
+        -fno-warn-unticked-promoted-constructors
+
+    if flag(use-fma)
+      cpp-options:
+          -DUSE_FMA
+      ghc-options:
+          -mfma
+      build-depends:
+          base
+            >= 4.19
+
+library
+
+    import:
+        common
+
+    hs-source-dirs:
+        src
+
+    default-language:
+        Haskell2010
+
+    exposed-modules:
+        Math.Algebra.Dual
+      , Math.Bezier.Cubic
+      , Math.Bezier.Cubic.Fit
+      , Math.Bezier.Quadratic
+      , Math.Bezier.Spline
+      , Math.Bezier.Stroke
+      , Math.Bezier.Stroke.EnvelopeEquation
+      , Math.Differentiable
+      , Math.Epsilon
+      , Math.Interval
+      , Math.Linear
+      , Math.Linear.Solve
+      , Math.Module
+      , Math.Monomial
+      , Math.Orientation
+      , Math.Ring
+      , Math.Roots
+      , Debug.Utils
+
+    other-modules:
+        Math.Algebra.Dual.Internal
+      , Math.Interval.Internal
+      , Math.Linear.Internal
+      , Math.Module.Internal
+      , TH.Utils
+
+    if flag(use-fma)
+      other-modules:
+        Math.Interval.FMA
+
+    build-depends:
+        bifunctors
+           >= 5.5.4      && < 5.7
+      , eigen
+          ^>= 3.3.7.0
+      , parallel
+          ^>= 3.2.2.0
+      , template-haskell
+           >= 2.18       && < 2.22

brush-strokes/cabal.project

@@ -0,0 +1,28 @@
+packages: .
+
+constraints:
+  acts -finitary,
+  rounded-hw -pure-hs -c99 -avx512 +ghc-prim -x87-long-double
+
+allow-newer:
+  acts:base, acts:deepseq,
+  groups-generic:base,
+  eigen:primitive,
+
+-------------
+-- GHC 9.4 --
+-------------
+
+-- eigen
+source-repository-package
+  type: git
+  location: https://github.com/chessai/eigen
+  tag: 1790fdf9138970dde0dbabf8b270698145a4a88c
+
+-------------
+-- GHC 9.6 --
+-------------
+
+-------------
+-- GHC 9.8 --
+-------------

brush-strokes/src/Debug/Utils.hs

@@ -1,6 +1,5 @@
 module Debug.Utils ( trace ) where
 
-
 -- base
 import System.IO
   ( BufferMode(..), hSetBuffering, hFlush, hPutStrLn, stdout )

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

@@ -6,9 +6,6 @@
 
 {-# OPTIONS_GHC -Wno-orphans -O2 #-}
 
-{-# OPTIONS_GHC -ddump-splices -ddump-simpl -ddump-cmm -ddump-to-file -dno-typeable-binds
-  -dsuppress-unfoldings -dsuppress-coercions #-}
-
 module Math.Algebra.Dual
   ( C(..), D
   , HasChainRule(..), chainRule

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

@@ -46,7 +46,12 @@ import qualified Control.Parallel.Strategies as Parallel.Strategy
 
 -- primitive
 import Data.Primitive.PrimArray
-  ( primArrayFromListN, unsafeThawPrimArray )
+  ( MutablePrimArray
+  , generatePrimArray
+  , primArrayFromListN
+  , unsafeThawPrimArray
+  , readPrimArray, writePrimArray
+  )
 
 -- transformers
 import Control.Monad.Trans.State.Strict
@@ -54,14 +59,6 @@ import Control.Monad.Trans.State.Strict
 import Control.Monad.Trans.Class
   ( lift )
 
--- vector
-import qualified Data.Vector.Unboxed.Mutable as Unboxed
-  ( MVector )
-import qualified Data.Vector.Unboxed.Mutable as Unboxed.MVector
-  ( unsafeRead, unsafeWrite )
-import qualified Data.Vector.Unboxed         as Unboxed.Vector
-  ( unsafeThaw, generate )
-
 -- MetaBrush
 import qualified Math.Bezier.Cubic as Cubic
   ( Bezier(..), bezier, ddist )
@@ -202,7 +199,7 @@ fitPiece
 fitPiece dist_tol t_tol maxIters p tp qs r tr =
   runST do
     -- Initialise the parameter values to a uniform subdivision.
-    ts <- Unboxed.Vector.unsafeThaw ( Unboxed.Vector.generate n uniform )
+    ts <- unsafeThawPrimArray $ generatePrimArray n uniform
     loop ts 0
   where
     n :: Int
@@ -216,12 +213,12 @@ fitPiece dist_tol t_tol maxIters p tp qs r tr =
     f2 t = h0 t *^ ( T p )
     f3 t = h3 t *^ ( T r )
 
-    loop :: forall s. Unboxed.MVector s Double -> Int -> ST s ( Cubic.Bezier ( ℝ 2 ), ArgMax Double Double )
+    loop :: forall s. MutablePrimArray s Double -> Int -> ST s ( Cubic.Bezier ( ℝ 2 ), ArgMax Double Double )
     loop ts count = do
       let
         hermiteParameters :: Mat22 -> T ( ℝ 2 ) -> Int -> [ ℝ 2 ] -> ST s ( T ( ℝ 2 ) )
         hermiteParameters ( Mat22 a11 a12 _ a22 ) ( V2 b1 b2 ) i ( q : rest ) = do
-          ti <- Unboxed.MVector.unsafeRead ts i
+          ti <- readPrimArray ts i
           let
             f0i, f1i, f2i, f3i :: T ( ℝ 2 )
             f0i = f0 ti
@@ -254,7 +251,7 @@ fitPiece dist_tol t_tol maxIters p tp qs r tr =
       -- so that t_i' is a better approximation of the parameter
       -- at which the curve is closest to q_i.
       ( dts_changed, argmax_sq_dist ) <- ( `execStateT` ( False, Max ( Arg 0 0 ) ) ) $ for_ ( zip qs [ 0 .. ] ) \( q, i ) -> do
-        ti <- lift ( Unboxed.MVector.unsafeRead ts i )
+        ti <- lift ( readPrimArray ts i )
         let
           laguerreStepResult :: Complex Double
           laguerreStepResult = runST do
@@ -277,7 +274,7 @@ fitPiece dist_tol t_tol maxIters p tp qs r tr =
           sq_dist :: Double
           sq_dist = quadrance @( T ( ℝ 2 ) ) q ( Cubic.bezier @( T ( ℝ 2 ) ) bez ti' )
         modify' ( second ( <> Max ( Arg sq_dist ti' ) ) )
-        lift ( Unboxed.MVector.unsafeWrite ts i ti' )
+        lift ( writePrimArray ts i ti' )
 
       case argmax_sq_dist of
         Max ( Arg max_sq_dist _ )

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

@@ -10,6 +10,8 @@ module Math.Bezier.Stroke
   , computeStrokeOutline, joinWithBrush
   , withTangent
 
+  , RootSolvingAlgorithm(..) -- TODO: move this?
+
     -- * Brush stroking
 
   , brushStroke, envelopeEquation
@@ -127,8 +129,6 @@ import Math.Orientation
   , between
   )
 import Math.Roots
-  ( solveQuadratic, newtonRaphson )
-
 
 import Debug.Utils
 
@@ -237,7 +237,8 @@ computeStrokeOutline ::
      , Show ptData, Show brushParams
 
      )
-  => FitParameters
+  => RootSolvingAlgorithm
+  -> FitParameters
   -> ( ptData -> usedParams )
   -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing
   -> ( forall {t} k (i :: t)
@@ -253,7 +254,7 @@ computeStrokeOutline ::
         , Seq FitPoint
         , [ Cusp ]
         )
-computeStrokeOutline fitParams ptParams toBrushParams brushFn spline@( Spline { splineStart = spt0 } ) = case ssplineType @clo of
+computeStrokeOutline rootAlgo 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
@@ -355,7 +356,7 @@ computeStrokeOutline fitParams ptParams toBrushParams brushFn spline@( Spline {
   where
 
     outlineInfo :: ptData -> Curve Open crvData ptData -> OutlineInfo
-    outlineInfo = inline ( outlineFunction ptParams toBrushParams brushFn )
+    outlineInfo = inline ( outlineFunction rootAlgo ptParams toBrushParams brushFn )
 
     outlineFns :: Seq OutlineInfo
     outlineFns = go spt0 ( openCurves $ splineCurves ( adjustSplineType @Open spline ) )
@@ -517,7 +518,8 @@ outlineFunction
      -- Debugging.
      , Show ptData, Show brushParams
      )
-  => ( ptData -> usedParams )
+  => RootSolvingAlgorithm
+  -> ( ptData -> usedParams )
   -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing
   -> ( forall {t} k (i :: t)
        .  DiffInterp k i brushParams
@@ -529,7 +531,7 @@ outlineFunction
   -> ptData
   -> Curve Open crvData ptData
   -> OutlineInfo
-outlineFunction ptParams toBrushParams brushFromParams = \ sp0 crv ->
+outlineFunction rootAlgo ptParams toBrushParams brushFromParams = \ sp0 crv ->
   let
 
     usedParams :: C 2 ( ℝ 1 ) usedParams
@@ -566,7 +568,7 @@ outlineFunction ptParams toBrushParams brushFromParams = \ sp0 crv ->
 
     fwdBwd :: OutlineFn
     fwdBwd t
-      = solveEnvelopeEquations t path_t path'_t ( fwdOffset, bwdOffset )
+      = solveEnvelopeEquations rootAlgo t path_t path'_t ( fwdOffset, bwdOffset )
           ( curves t )
       where
 
@@ -953,15 +955,23 @@ splineCurveFns co spls
       Bezier3To { controlPoint1 = p1, controlPoint2 = p2, curveEnd = NextPoint p3 }
         -> bezier3 @k @i @( ℝ 2 ) co $ Cubic.Bezier p0 p1 p2 p3
 
+-- | Which method to use to solve the envelope equation?
+data RootSolvingAlgorithm
+  -- | Use the Newton–Raphson method.
+  = NewtonRaphson { maxIters :: Word, precision :: Int }
+  -- | Use the modified Halley M2 method.
+  | HalleyM2 { maxIters :: Word, precision :: Int }
+
 -- | 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 :: ℝ 1 -- ^ @t@ (for debugging only)
+solveEnvelopeEquations :: RootSolvingAlgorithm
+                       -> ℝ 1 -- ^ @t@ (for debugging only)
                        -> ℝ 2
                        -> T ( ℝ 2 )
                        -> ( Offset, Offset )
                        -> Seq ( ℝ 1 -> StrokeDatum 2 () )
                        -> ( ( ℝ 2, T ( ℝ 2 ) ), ( ℝ 2, T ( ℝ 2 ) ) )
-solveEnvelopeEquations _t path_t path'_t ( fwdOffset, bwdOffset ) strokeData
+solveEnvelopeEquations rootAlgo _t path_t path'_t ( fwdOffset, bwdOffset ) strokeData
   = ( fwdSol, ( bwdPt, -1 *^ bwdTgt ) )
 
   where
@@ -999,12 +1009,18 @@ solveEnvelopeEquations _t path_t path'_t ( fwdOffset, bwdOffset ) strokeData
     sol :: Seq ( ℝ 1 -> StrokeDatum 2 () ) -> Double -> ( Bool, ℝ 2, T ( ℝ 2 ) )
     sol f is0 =
       let ( good, is ) =
-            case newtonRaphson maxIters precision domain ( eqn f ) is0 of
+            case runSolveMethod ( eqn f ) is0 of
               Nothing  -> ( False, is0 )
               Just is1 -> ( True , is1 )
           ( ds, dcdt ) = finish f is
       in ( good, ds, dcdt )
 
+    runSolveMethod = case rootAlgo of
+      HalleyM2 { maxIters, precision } ->
+        halleyM2 maxIters precision
+      NewtonRaphson { maxIters, precision } ->
+        newtonRaphson maxIters precision domain
+
     finish :: Seq ( ℝ 1 -> StrokeDatum 2 () ) -> Double -> ( ℝ 2, T ( ℝ 2 ) )
     finish f is =
       let (i, s) = fromDomain is in
@@ -1037,17 +1053,13 @@ solveEnvelopeEquations _t path_t path'_t ( fwdOffset, bwdOffset ) strokeData
              --    ] )
               ( value @Double @2 @( ℝ 2 ) dstroke, dcdt )
 
-    eqn :: Seq ( ℝ 1 -> StrokeDatum 2 () ) -> ( Double -> ( Double, Double ) )
+    eqn :: Seq ( ℝ 1 -> StrokeDatum 2 () ) -> ( Double -> (# Double, Double #) )
     eqn fs is =
       let (i, s) = fromDomain is
       in case ( fs `Seq.index` i ) ( ℝ1 s ) of
            StrokeDatum { ee = D12 ee _ ee_s } ->
-             coerce ( ee, ee_s )
+             coerce (# ee, ee_s #)
 
-    maxIters :: Word
-    maxIters = 20
-    precision :: Int
-    precision = 8
     n :: Int
     n = Seq.length strokeData
     domain :: ( Double, Double )

brush-strokes/src/Math/Interval.hs

@@ -4,9 +4,6 @@
 
 {-# OPTIONS_GHC -Wno-orphans #-}
 
-{-# OPTIONS_GHC -ddump-splices -ddump-simpl -ddump-cmm -ddump-to-file -dno-typeable-binds
-  -dsuppress-unfoldings -dsuppress-coercions #-}
-
 module Math.Interval
   ( 𝕀(𝕀), inf, sup
   , scaleInterval
@@ -32,7 +29,7 @@ import Data.Group
 import Data.Group.Generics
   ( )
 
--- splines
+-- brush-strokes
 import Math.Algebra.Dual
 import Math.Interval.Internal
   ( 𝕀(𝕀), inf, sup, scaleInterval )

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

@@ -136,6 +136,7 @@ instance Transcendental ( 𝕀 Double ) where
   atan = withHW Prelude.atan
 
 {-# INLINE withHW #-}
+-- | Internal function: use @rounded-hw@ to define a function on intervals.
 withHW :: (Interval.Interval a -> Interval.Interval b) -> 𝕀 a -> 𝕀 b
 withHW f = \ ( 𝕀 lo hi ) ->
   case f ( Interval.I ( Rounded lo ) ( Rounded hi ) ) of

brush-strokes/src/Math/Roots.hs

@@ -13,6 +13,9 @@ module Math.Roots
 
     -- * Newton–Raphson
   , newtonRaphson
+
+    -- * Modified Halley's method M2
+  , halleyStep, halleyM2
   )
   where
 
@@ -43,10 +46,6 @@ import Data.Primitive.PrimArray
 import Data.Primitive.Types
   ( Prim )
 
--- prim-instances
-import Data.Primitive.Instances
-  () -- instance Prim a => Prim ( Complex a )
-
 -- MetaBrush
 import Math.Epsilon
   ( epsilon, nearZero )
@@ -299,14 +298,14 @@ 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
   #-}
 {-# INLINEABLE newtonRaphson #-}
 newtonRaphson :: ( RealFloat r, Show r )
               => Word               -- ^ maximum number of iterations
               -> Int                -- ^ desired digits of precision
               -> ( r, r )           -- ^ @(min_x, max_x)@.
-              -> ( r -> ( r, r ) ) -- ^ function and its derivative
+              -> ( r -> (# r, r #) ) -- ^ function and its derivative
               -> r                  -- ^ initial guess
               -> Maybe r
 newtonRaphson maxIters digits ( min_x, max_x ) f x0 =
@@ -315,7 +314,7 @@ newtonRaphson maxIters digits ( min_x, max_x ) f x0 =
     !factor = encodeFloat 1 ( 1 - digits )
 
 doNewtonRaphson :: ( Fractional r, Ord r, Show r )
-                => ( r -> (r, r) )
+                => ( r -> (# r, r #) )
                 -> Word
                 -> r
                 -> r -> r
@@ -329,7 +328,7 @@ doNewtonRaphson f maxIters factor min_x max_x min_f_x max_f_x f_x_prev 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, 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
@@ -349,7 +348,7 @@ 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 ( iters + 1 ) new_x δ δ1
 
 newtonRaphsonStep :: ( Fractional r, Ord r, Show r )
-                  => ( r -> ( r, r ) )
+                  => ( r -> (# r, r #) )
                   -> r
                   -> r -> r -> r
                   -> r -> r
@@ -380,7 +379,7 @@ newtonRaphsonStep f min_x max_x f_x_prev x f_x f'_x δ1 δ2
 
 -- Handle \( f'(x_0) = 0 \).
 zeroDerivativeStep :: ( Fractional r, Ord r, Show r )
-                   => ( r -> ( r, r ) )
+                   => ( r -> (# r, r #) )
                    -> r
                    -> r -> r
                    -> r
@@ -390,7 +389,7 @@ zeroDerivativeStep f min_x max_x f_x_prev x f_x δ
   -- Handle zero derivative on first iteration.
   | f_x_prev == 0
   , x_prev <- if x <= 0.5 * ( min_x + max_x ) then max_x else min_x
-  , ( f_x_prev, _ ) <- f x_prev
+  , (# f_x_prev, _ #) <- f x_prev
   , δ <- x_prev - x
   = finish f_x_prev δ
   | otherwise
@@ -430,3 +429,52 @@ maxDouble :: Double
 maxDouble = D# 1.7976931348623157e308##
 
 -}
+
+--------------------------------------------------------------------------------
+-- Halley's method (WIP).
+
+-- | Take a single step with Halley's method.
+halleyStep :: Fractional a => (# a, a, a, a #) -> a
+halleyStep (# x, f, f', f'' #) =
+  x - 2 * f * f' / ( 2 * f' ^ ( 2 :: Int ) - f * f'' )
+
+-- | Take a single step in the M2 modified Halley method.
+--
+-- Taken from @Some variants of Halley’s method with memory and their applications for solving several chemical problems@
+-- 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 (# x_nm1, (# f_nm1, f'_nm1 #) #) (# x_n, (# f_n, f'_n #) #) =
+  num / denom
+  where
+    u = f_n * f_nm1 * (f_n - f_nm1)
+    dx = x_n - x_nm1
+    g1 = f_nm1 ^ ( 2 :: Int ) * f'_n
+    g2 = f_n ^ ( 2 :: Int ) * f'_nm1
+    num = (x_n + x_nm1) * u - dx * ( g1 * x_n + g2 * x_nm1)
+    denom = 2 * u - dx * ( g1 + g2 )
+
+{-# SPECIALISE
+  halleyM2
+    :: Word -> Int -> ( Double -> (# Double, Double #) ) -> Double -> Maybe Double
+  #-}
+{-# INLINEABLE halleyM2 #-}
+halleyM2 :: ( RealFloat r, Show r )
+         => Word              -- ^ maximum number of iterations
+         -> Int               -- ^ desired digits of precision
+         -> ( r -> (# r, r #) ) -- ^ function and its derivative
+         -> r                 -- ^ initial guess
+         -> Maybe r
+halleyM2 maxIters digits f x0 =
+  let y0 = (# x0, f x0 #)
+  in go 0 y0 y0
+  where
+    !factor = encodeFloat 1 ( 1 - digits )
+    go i y_nm1 y_n@(# x_n, _ #) =
+      let x_np1 = halleyM2Step y_nm1 y_n
+      in if | i >= maxIters
+            || abs ( x_np1 - x_n ) < abs ( x_n * factor )
+            -> Just x_np1
+            | otherwise
+            -> go (i+1) y_n (# x_n, f x_np1 #)

cabal.project

@@ -1,9 +1,12 @@
-packages: .
+packages: brush-strokes
+        , .
 
 constraints:
   acts -finitary,
-  rounded-hw -pure-hs -c99 -avx512 +ghc-prim -x87-long-double
+  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
 
 -- Fix a severe bug in Waargonaut (no corresponding Hackage release???)
 source-repository-package
@@ -12,27 +15,23 @@ source-repository-package
   tag: 5f838582a8c5aae1a198ecd4958729e53a6b03cf
 
 allow-newer:
-  *:base, *:ghc, *:ghc-prim, *:template-haskell,
+  *:base, *:template-haskell, *:ghc-prim,
-  *:text
+  acts:deepseq,
-
--------------
--- GHC 9.2 --
--------------
-
-allow-newer:
   digit:lens,
+  eigen:primitive,
+  eigen:transformers,
+  gi-cairo-connector:mtl,
   hedgehog:resourcet,
+  JuicyPixels:zlib,
   natural:lens,
+  records-sop:deepseq,
+  waargonaut:bifunctors,
   waargonaut:lens,
   waargonaut:records-sop,
-  waargonaut:semigroups, waargonaut:text,
+  waargonaut:semigroups,
-  waargonaut:vector, waargonaut:witherable,
+  waargonaut:text,
-
+  waargonaut:vector,
--- records-sop
+  waargonaut:witherable,
-source-repository-package
-  type: git
-  location: https://github.com/kosmikus/records-sop
-  tag: abab99b4b870fce55e81dd03d4e41fb50502ca4e
 
 -------------
 -- GHC 9.4 --
@@ -44,13 +43,6 @@ source-repository-package
   location: https://github.com/chessai/eigen
   tag: 1790fdf9138970dde0dbabf8b270698145a4a88c
 
--- rounded-hw
-source-repository-package
-  type: git
-  location: https://github.com/minoki/haskell-floating-point
-  subdir: rounded-hw
-  tag: c141bc8435c556e94d51254479a8213698679261
-
 -------------
 -- GHC 9.6 --
 -------------

src/app/MetaBrush/Application.hs

@@ -78,7 +78,7 @@ import Math.Bezier.Cubic.Fit
 import Math.Bezier.Spline
   ( Spline(..), Curves(..), Curve(..), NextPoint(..) )
 import Math.Bezier.Stroke
-  ( invalidateCache )
+  ( RootSolvingAlgorithm(..), invalidateCache )
 import Math.Linear
   ( ℝ(..) )
 import MetaBrush.Action
@@ -200,15 +200,16 @@ runApplication application = do
   fileBarTabsTVar      <- STM.newTVarIO @( Map Unique FileBarTab )               Map.empty
   showGuidesTVar       <- STM.newTVarIO @Bool                                    True
   maxHistorySizeTVar   <- STM.newTVarIO @Int                                     1000
-  fitParametersTVar    <- STM.newTVarIO @FitParameters
+  fitParametersTVar    <- STM.newTVarIO @FitParameters $
-                            ( FitParameters
+                            FitParameters
                               { maxSubdiv  = 5 --2 --3 -- 6
                               , nbSegments = 3 --3 --6 -- 12
                               , dist_tol   = 0.1 -- 5e-3
                               , t_tol      = 0.1 -- 1e-4
                               , maxIters   = 5   -- 100
                               }
-                            )
+  rootsAlgoTVar        <- STM.newTVarIO @RootSolvingAlgorithm $
+                            HalleyM2 { maxIters = 20, precision = 8 }
 
   -- Put all these stateful variables in a record for conciseness.
   let
@@ -333,6 +334,7 @@ runApplication application = do
                 showGuides     <- STM.readTVar showGuidesTVar
                 debug          <- STM.readTVar debugTVar
                 fitParameters  <- STM.readTVar fitParametersTVar
+                rootsAlgo      <- STM.readTVar rootsAlgoTVar
                 STM.writeTVar recomputeStrokesTVar False
                 let
                   addRulers :: ( ( Int32, Int32 ) -> Cairo.Render () ) -> ( ( Int32, Int32 ) -> Cairo.Render () ) 
@@ -344,7 +346,7 @@ runApplication application = do
                       doc
                 pure
                   ( addRulers <$> getDocumentRender
-                        colours fitParameters mode debug
+                      colours rootsAlgo fitParameters mode debug
                       modifiers mbMousePos mbHoldAction mbPartialPath
                       doc
                   )

src/app/MetaBrush/Context.hs

@@ -41,6 +41,8 @@ import Data.HashMap.Strict
 -- MetaBrush
 import Math.Bezier.Cubic.Fit
   ( FitParameters )
+import Math.Bezier.Stroke
+  ( RootSolvingAlgorithm )
 import Math.Linear
   ( ℝ(..) )
 import {-# SOURCE #-} MetaBrush.Action
@@ -99,6 +101,7 @@ data Variables
   , showGuidesTVar       :: !( STM.TVar Bool                                    )
   , maxHistorySizeTVar   :: !( STM.TVar Int                                     )
   , fitParametersTVar    :: !( STM.TVar FitParameters                           )
+  , rootsAlgoTVar        :: !( STM.TVar RootSolvingAlgorithm                    )
   }
 
 --------------------------------------------------------------------------------

src/app/MetaBrush/Render/Document.hs

@@ -74,6 +74,7 @@ import Math.Bezier.Spline
 import Math.Bezier.Stroke
   ( Cusp(..), CachedStroke(..), invalidateCache
   , computeStrokeOutline
+  , RootSolvingAlgorithm
   )
 import Math.Linear
   ( ℝ(..), T(..) )
@@ -147,12 +148,12 @@ blankRender :: Colours -> Cairo.Render ()
 blankRender _ = pure ()
 
 getDocumentRender
-  :: Colours -> FitParameters -> Mode -> Bool
+  :: Colours -> RootSolvingAlgorithm -> FitParameters -> Mode -> Bool
   -> Set Modifier -> Maybe ( ℝ 2 ) -> Maybe HoldAction -> Maybe PartialPath
   -> Document
   -> ST RealWorld ( ( Int32, Int32 ) -> Cairo.Render () )
 getDocumentRender
-  cols fitParams mode debug
+  cols rootAlgo fitParams mode debug
   modifiers mbMousePos mbHoldEvent mbPartialPath
   doc@( Document { viewportCenter = ℝ2 cx cy, zoomFactor, documentContent = content } )
     = do
@@ -214,7 +215,7 @@ getDocumentRender
               -> previewStroke :<| foldMap visibleStrokes ( strokes content )
           _ -> foldMap visibleStrokes ( strokes content )
 
-      strokesRenderData <- traverseMaybe ( sequenceA . strokeRenderData fitParams ) modifiedStrokes
+      strokesRenderData <- traverseMaybe ( sequenceA . strokeRenderData rootAlgo fitParams ) modifiedStrokes
 
       let
         renderSelectionRect :: Cairo.Render ()
@@ -275,8 +276,8 @@ instance NFData StrokeRenderData where
 --    - the computed outline (using fitting algorithm),
 --    - the brush shape function.
 -- - Otherwise, this consists of the underlying spline path only.
-strokeRenderData :: FitParameters -> Stroke -> Maybe ( ST RealWorld StrokeRenderData )
+strokeRenderData :: RootSolvingAlgorithm -> FitParameters -> Stroke -> Maybe ( ST RealWorld StrokeRenderData )
-strokeRenderData fitParams
+strokeRenderData rootAlgo fitParams
   ( Stroke
     { strokeSpline = spline :: StrokeSpline clo ( Record pointFields )
     , strokeBrush  = ( strokeBrush :: Maybe ( Brush brushFields ) )
@@ -301,7 +302,7 @@ strokeRenderData fitParams
 
                   -- Compute the outline using the brush function.
                     ( outline, fitPts, cusps ) <-
-                      computeStrokeOutline @clo fitParams
+                      computeStrokeOutline @clo rootAlgo fitParams
                         ( toUsedParams . brushParams ) embedUsedParams brushFn
                         spline
                     pure $

src/app/MetaBrush/UI/InfoBar.hs

@@ -179,7 +179,8 @@ updateInfoBar viewportDrawingArea ( InfoBar {..} ) ( Variables { mousePosTVar }
         GTK.labelSetText botRightPosText $ Text.pack ( "x: " <> fixed 6 2  r <> "\ny: " <> fixed 6 2  b )
 
 fixed :: Int -> Int -> Double -> String
-fixed digitsBefore digitsAfter x = case second tail . break ( == '.' ) $ showFFloat ( Just digitsAfter ) x "" of
+fixed digitsBefore digitsAfter x =
+  case second (drop 1) . break ( == '.' ) $ showFFloat ( Just digitsAfter ) x "" of
     ( as, bs ) ->
       let
         l, r :: Int

src/convert/Main.hs

@@ -34,7 +34,7 @@ import qualified Linear.V2 as Linear
 import Control.Monad.Trans.Reader
   ( runReaderT )
 
--- splines
+-- brush-strokes
 import Math.Bezier.Spline
   ( Spline, SplineType(..) )
 import Math.Bezier.Stroke

src/convert/MetaBrush/MetaFont/Convert.hs

@@ -56,7 +56,7 @@ import qualified Text.Parsec.Error as Parsec
 import Data.Text
   ( Text )
 
--- splines
+-- brush-strokes
 import Math.Bezier.Spline
   ( Spline(..), SplineType(..), KnownSplineType
   , Curves(..), Curve(..), NextPoint(..)