experiment: FMA backend for interval arithmetic

also includes the start of a way to observe which equations
are being solved, which should help with improving the performance of
the interval Newton method

MetaBrush.cabal

@@ -29,6 +29,11 @@ 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:
@@ -53,7 +58,7 @@ common common
       , rounded-hw
           ^>= 0.3
       , transformers
-          ^>= 0.5.6.2
+           >= 0.5.6.2 && < 0.7
 
     default-extensions:
         BangPatterns
@@ -109,6 +114,14 @@ 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
 
@@ -201,6 +214,10 @@ library splines
       , Math.Module.Internal
       , TH.Utils
 
+    if flag(use-fma)
+      other-modules:
+        Math.Interval.FMA
+
     build-depends:
         bifunctors
            >= 5.5.4      && < 5.6
@@ -248,6 +265,19 @@ library metabrushes
       , bytestring
            >= 0.10.10.0  && < 0.12
 
+executable test-cusp-isolation
+
+    import:
+        common
+
+    hs-source-dirs:
+        src/cusps
+
+    build-depends:
+        splines
+
+    main-is:
+        Main.hs
 
 executable convert-metafont
 
@@ -261,7 +291,6 @@ executable convert-metafont
         Haskell2010
 
     main-is:
-
         Main.hs
 
     other-modules:

cabal.project

@@ -11,6 +11,10 @@ source-repository-package
   location: https://github.com/haskell-waargonaut/waargonaut
   tag: 5f838582a8c5aae1a198ecd4958729e53a6b03cf
 
+allow-newer:
+  *:base, *:ghc, *:ghc-prim, *:template-haskell,
+  *:text
+
 -------------
 -- GHC 9.2 --
 -------------
@@ -18,11 +22,8 @@ source-repository-package
 allow-newer:
   digit:lens,
   hedgehog:resourcet,
-  hoist-error:base,
   natural:lens,
-  prim-instances:base,
+  waargonaut:lens,
-  records-sop:ghc-prim,
-  waargonaut:base, waargonaut:lens,
   waargonaut:records-sop,
   waargonaut:semigroups, waargonaut:text,
   waargonaut:vector, waargonaut:witherable,
@@ -54,9 +55,6 @@ source-repository-package
 -- GHC 9.6 --
 -------------
 
-allow-newer:
+-------------
-  attoparsec:ghc-prim,
+-- GHC 9.8 --
-  generics-sop:base, generics-sop:ghc-prim,
+-------------
-  hw-prim:ghc-prim,
-  integer-logarithms:ghc-prim,
-  vector-stream:base, vector-stream:ghc-prim,

src/app/MetaBrush/Application.hs

@@ -170,8 +170,8 @@ runApplication application = do
                 { splineStart  = mkPoint ( ℝ2 10 -20 ) 2 1 0
                 , splineCurves = OpenCurves $ Seq.fromList
                   [ LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2  10  10 ) 10 5 ( pi / 4 ) ), curveData = invalidateCache undefined }
-                  , LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 -10  10 )  8 5 ( pi / 4 ) ), curveData = invalidateCache undefined }
+                  --, LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 -10  10 )  8 5 ( pi / 4 ) ), curveData = invalidateCache undefined }
-                  , LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 -10 -20 ) 10 7 ( pi / 2 ) ), curveData = invalidateCache undefined }
+                  --, LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 -10 -20 ) 10 7 ( pi / 2 ) ), curveData = invalidateCache undefined }
                   ]
                 }
             }

src/app/MetaBrush/Render/Document.hs

@@ -75,8 +75,6 @@ import Math.Bezier.Stroke
   ( Cusp(..), CachedStroke(..), invalidateCache
   , computeStrokeOutline
   )
-import Math.Interval
-  ( Extent(Point) )
 import Math.Linear
   ( ℝ(..), T(..) )
 import Math.Module
@@ -310,7 +308,7 @@ strokeRenderData fitParams
                       StrokeWithOutlineRenderData
                         { strokeDataSpline    = spline
                         , strokeOutlineData   = ( outline, fitPts, cusps )
-                        , strokeBrushFunction = fun @Double ( brushFn @Point proxy# id )
+                        , strokeBrushFunction = fun @Double ( brushFn @2 @() proxy# id )
                                               . embedUsedParams
                                               . toUsedParams
                         }

src/cusps/Main.hs

@@ -0,0 +1,209 @@
+{-# LANGUAGE AllowAmbiguousTypes #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+module Main where
+
+-- base
+import Prelude hiding ( Num(..), Fractional(..), Floating(..), (^) )
+import qualified Prelude ( Num(..), Fractional(..) )
+import Data.Coerce
+  ( Coercible )
+import GHC.Generics
+  ( Generic )
+
+-- acts
+import Data.Act
+  ( Torsor )
+
+-- containers
+import Data.Map
+  ( Map )
+import qualified Data.Map as Map
+  ( (!) )
+import qualified Data.Sequence as Seq
+  ( index )
+
+-- generic-lens
+import Data.Generics.Product.Typed
+  ( HasType )
+
+-- MetaBrush
+import Math.Algebra.Dual
+import Math.Bezier.Spline
+import Math.Bezier.Stroke
+  ( brushStrokeData, pathAndUsedParams )
+import Math.Bezier.Stroke.EnvelopeEquation
+  ( StrokeDatum(..) )
+import Math.Differentiable
+  ( I )
+import Math.Interval
+import Math.Linear
+  ( ℝ(..), T(..) )
+import Math.Module
+  ( Module, Cross )
+import Math.Ring
+  ( AbelianGroup(..), Ring(..), Field(..), Transcendental(..)
+  , ViaPrelude(..)
+  )
+
+--------------------------------------------------------------------------------
+
+data PointData params
+  = PointData
+  { pointCoords :: !( ℝ 2 )
+  , brushParams :: !params
+  }
+  deriving stock ( Show, Generic )
+
+
+{-
+
+newtype C k u v = D { runD :: u -> D k u v }
+
+type instance D k ( ℝ n ) = Dk𝔸n
+type instance D k ( 𝕀 v ) = D k v
+
+e.g.
+
+data D2𝔸2 v =
+  D22 { _D22_v :: !v
+      , _D22_dx, _D22_dy :: !( T v )
+      , _D22_dxdx, _D22_dxdy, _D22_dydy :: !( T v )
+      }
+
+-}
+
+
+outlineFunction
+  :: forall i brushParams
+  .  ( Show brushParams
+     , D 1 ( I i ( ℝ 2 ) ) ~ D 1 ( ℝ 2 )
+     , D 2 ( I i ( ℝ 2 ) ) ~ D 2 ( ℝ 2 )
+     , D 3 ( I i ( ℝ 1 ) ) ~ D 3 ( ℝ 1 )
+     , D 3 ( I i ( ℝ 2 ) ) ~ D 3 ( ℝ 2 )
+     , Coercible ( I i ( ℝ 1 ) ) ( I i Double )
+     , HasType ( ℝ 2 ) ( PointData brushParams )
+     , Cross  ( I i Double ) ( T ( I i ( ℝ 2 ) ) )
+     , Module ( I i Double ) ( T ( I i brushParams ) )
+     , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) )
+     , Torsor ( T ( I i brushParams ) ) ( I i brushParams )
+     , HasChainRule ( I i Double ) 3 ( I i ( ℝ 1 ) )
+     , HasChainRule ( I i Double ) 3 ( I i brushParams )
+     , Traversable ( D 3 brushParams )
+     , Traversable ( D 3 ( I i brushParams ) )
+     )
+  => ( forall a. a -> I i a )
+  -> C 3 ( I i brushParams ) ( Spline Closed () ( I i ( ℝ 2 ) ) )
+      -- ^ brush shape
+  -> Int -- ^ brush segment index to consider
+  -> PointData brushParams -> Curve Open () ( PointData brushParams )
+      -- ^ stroke path
+  -> ( I i ( ℝ 1 ) -> I i ( ℝ 1 ) -> StrokeDatum 3 i )
+outlineFunction single brush brush_index sp0 crv = strokeData
+
+  where
+    path   :: C 3 ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) )
+    params :: C 3 ( I i ( ℝ 1 ) ) ( I i brushParams )
+    (path, params) =
+      pathAndUsedParams @3 @i single brushParams sp0 crv
+
+    strokeData :: I i ( ℝ 1 ) -> I i ( ℝ 1 ) -> StrokeDatum 3 i
+    strokeData =
+      fmap ( `Seq.index` brush_index ) $
+        brushStrokeData @3 @brushParams @i path params brush
+
+main :: IO ()
+main = return ()
+
+--------------------------------------------------------------------------------
+-- EDSL for inspection
+
+type instance I AA a = AA a
+
+data AAVal a where
+  Val :: a -> AAVal a
+  Var :: String -> AAVal Double
+
+instance Show a => Show ( AAVal a ) where
+  show (Val v) = show v
+  show (Var v) = show v
+
+data AA a where
+  Pt     :: AAVal a -> AA a
+  Ival   :: AAVal a -> AAVal a -> AA a
+  (:+:)  :: AA Double -> AA Double -> AA Double
+  (:-:)  :: AA Double -> AA Double -> AA Double
+  Negate :: AA Double -> AA Double
+  (:*:)  :: AA Double -> AA Double -> AA Double
+  (:^:)  :: AA Double -> Word -> AA Double
+  Recip  :: AA Double -> AA Double
+  Cos    :: AA Double -> AA Double
+  Sin    :: AA Double -> AA Double
+  Pi     :: AA Double
+
+instance Show a => Show (AA a) where
+  showsPrec prec = \case
+    Pt v        -> showString "[" . showsPrec 0 v . showString "]"
+    Ival lo hi  -> showString "[" . showsPrec 0 lo . showString "," . showsPrec 0 hi . showString "]"
+    iv1 :+: iv2 -> showParen (prec >  6) $ showsPrec 6 iv1 . showString " + " . showsPrec 7 iv2
+    iv1 :-: iv2 -> showParen (prec >  6) $ showsPrec 6 iv1 . showString " - " . showsPrec 7 iv2
+    iv1 :*: iv2 -> showParen (prec >  7) $ showsPrec 7 iv1 . showString " - " . showsPrec 8 iv2
+    iv :^: n    -> showParen (prec >  8) $ showsPrec 9 iv  . showString " ^ " . showsPrec 8 n
+    Negate iv   -> showParen (prec > 10) $ showString "negate " . showsPrec 11 iv
+    Recip iv    -> showParen (prec > 10) $ showString "recip " . showsPrec 11 iv
+    Cos iv      -> showParen (prec > 10) $ showString "cos "   . showsPrec 11 iv
+    Sin iv      -> showParen (prec > 10) $ showString "sin "   . showsPrec 11 iv
+    Pi          -> showString "pi"
+
+infixl 6 :+:
+infixl 6 :-:
+infixl 7 :*:
+infixr 8 :^:
+
+val :: Map String Double -> AAVal a -> a
+val _    ( Val a ) = a
+val vars ( Var v ) = vars Map.! v
+
+eval :: Map String Double -> AA a -> 𝕀 a
+eval vars = \ case
+  Pt v        -> let f = val vars v in 𝕀 f f
+  Ival lo hi  -> 𝕀 (val vars lo) (val vars hi)
+  iv1 :+: iv2 -> eval vars iv1 + eval vars iv2
+  iv1 :-: iv2 -> eval vars iv1 - eval vars iv2
+  Negate iv   -> negate $ eval vars iv
+  iv1 :*: iv2 -> eval vars iv1 * eval vars iv2
+  iv :^: n    -> eval vars iv ^ n
+  Recip iv    -> recip $ eval vars iv
+  Cos iv      -> cos $ eval vars iv
+  Sin iv      -> sin $ eval vars iv
+  Pi          -> pi
+
+instance Prelude.Num ( AA Double ) where
+  (+) = (:+:)
+  (-) = (:-:)
+  negate = Negate
+  (*) = (:*:)
+  abs = error "No abs for abstract intervals"
+  signum = error "No signum for abstract intervals"
+  fromInteger = Pt . Val . fromInteger
+
+instance Prelude.Fractional ( AA Double ) where
+  recip = Recip
+  (/) = error "No division for abstract intervals"
+  fromRational = Pt . Val . fromRational
+
+instance Ring ( AA Double ) where
+  (*) = (:*:)
+  (^) = (:^:)
+
+deriving via ViaPrelude ( AA Double )
+  instance AbelianGroup ( AA Double )
+deriving via ViaPrelude ( AA Double )
+  instance AbelianGroup ( T ( AA Double ) )
+deriving via ViaPrelude ( AA Double )
+  instance Field ( AA Double )
+
+instance Transcendental ( AA Double ) where
+  pi  = Pi
+  cos = Cos
+  sin = Sin

src/metabrushes/MetaBrush/Asset/Brushes.hs

@@ -1,5 +1,6 @@
 {-# LANGUAGE AllowAmbiguousTypes #-}
 {-# LANGUAGE OverloadedStrings   #-}
+{-# LANGUAGE PolyKinds           #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 
 module MetaBrush.Asset.Brushes
@@ -32,8 +33,8 @@ import qualified Data.HashMap.Strict as HashMap
 -- MetaBrush
 import Math.Algebra.Dual
 import Math.Bezier.Spline
-import Math.Interval
+import Math.Differentiable
-
+  ( I )
 import Math.Linear
 import Math.Module
   ( Module((^+^), (*^)) )
@@ -81,7 +82,7 @@ ellipse = BrushData "ellipse" ( WithParams deflts ellipseBrush )
 --------------------------------------------------------------------------------
 -- Differentiable brushes.
 
-circleSpline :: forall i k u v
+circleSpline :: forall {t} (i :: t) k u v
              .  Applicative ( D k ( I i u ) )
              => ( Double -> Double -> D k ( I i u ) ( I i v ) )
              -> D k ( I i u ) ( Spline 'Closed () ( I i v ) )
@@ -98,7 +99,7 @@ circleSpline p = sequenceA $
         Bezier3To ( p  κ -1 ) ( p  1 -κ ) BackToStart ()
 {-# INLINE circleSpline #-}
 
-circleBrush :: forall i k irec
+circleBrush :: forall {t} (i :: t) k irec
             . ( irec ~ I i ( Record CircleBrushFields )
               , Module
                  ( D k irec ( I i Double ) )
@@ -128,7 +129,7 @@ circleBrush _ mkI =
     kon = konst @( I i Double ) @k @irec . mkI
 {-# INLINEABLE circleBrush #-}
 
-ellipseBrush :: forall i k irec
+ellipseBrush :: forall {t} (i :: t) k irec
              . ( irec ~ I i ( Record EllipseBrushFields )
                , Module
                   ( D k irec ( I i Double ) )

src/metabrushes/MetaBrush/Brush.hs

@@ -1,4 +1,5 @@
 {-# LANGUAGE AllowAmbiguousTypes   #-}
+{-# LANGUAGE PolyKinds             #-}
 {-# LANGUAGE QuantifiedConstraints #-}
 {-# LANGUAGE ScopedTypeVariables   #-}
 {-# LANGUAGE UndecidableInstances  #-}
@@ -44,9 +45,9 @@ import Math.Algebra.Dual
 import Math.Bezier.Spline
   ( SplineType(Closed), Spline )
 import Math.Differentiable
-  ( DiffInterp, ExtentOrder )
+  ( DiffInterp, I )
 import Math.Interval
-  ( type I, Extent(Point, Interval) )
+  ( 𝕀 )
 import Math.Linear
 import MetaBrush.Records
   ( KnownSymbols, Length, Record )
@@ -62,11 +63,11 @@ data WithParams params f =
   WithParams
     { defaultParams :: params
     , withParams
-        :: forall i
+        :: forall {t} k (i :: t)
-        .  ( DiffInterp i params )
+        .  ( DiffInterp k i params )
         => Proxy# i
         -> ( forall a. a -> I i a )
-        -> C ( ExtentOrder i ) ( I i params ) ( f ( I i ( ℝ 2 ) ) )
+        -> C k ( I i params ) ( f ( I i ( ℝ 2 ) ) )
     }
 
 --------------------------------------------------------------------------------
@@ -82,8 +83,8 @@ data Brush brushFields where
     :: forall brushFields
     .  ( KnownSymbols brushFields, Typeable brushFields
        , Representable Double ( ℝ ( Length brushFields ) )
-       , DiffInterp Point     ( ℝ ( Length brushFields ) )
+       , DiffInterp 2 () ( ℝ ( Length brushFields ) )
-       , DiffInterp Interval  ( ℝ ( Length brushFields ) )
+       , DiffInterp 3 𝕀  ( ℝ ( Length brushFields ) )
        )
     => { brushName     :: !Text
        , brushFunction :: BrushFunction brushFields
@@ -118,8 +119,8 @@ class ( KnownSymbols pointFields, Typeable pointFields
       , Show           ( Record pointFields )
       , NFData         ( Record pointFields )
       , Representable  Double   ( ℝ ( Length pointFields ) )
-      , DiffInterp Point    ( ℝ ( Length pointFields ) )
+      , DiffInterp 2 () ( ℝ ( Length pointFields ) )
-      , DiffInterp Interval ( ℝ ( Length pointFields ) )
+      , DiffInterp 3 𝕀  ( ℝ ( Length pointFields ) )
       )
    => PointFields pointFields where { }
 instance ( KnownSymbols pointFields, Typeable pointFields
@@ -127,8 +128,8 @@ instance ( KnownSymbols pointFields, Typeable pointFields
          , Show           ( Record pointFields )
          , NFData         ( Record pointFields )
          , Representable  Double   ( ℝ ( Length pointFields ) )
-         , DiffInterp Point    ( ℝ ( Length pointFields ) )
+         , DiffInterp 2 () ( ℝ ( Length pointFields ) )
-         , DiffInterp Interval ( ℝ ( Length pointFields ) )
+         , DiffInterp 3 𝕀  ( ℝ ( Length pointFields ) )
          )
       => PointFields pointFields where { }
 

src/metabrushes/MetaBrush/Records.hs

@@ -181,8 +181,8 @@ intersect :: forall r1 r2 l1 l2
             , l1 ~ Length r1, l2 ~ Length r2
             , Representable Double ( ℝ l1 )
             , Representable Double ( ℝ l2 )
-            , Differentiable 'Point    ( ℝ l2 )
+            , Differentiable 2 () ( ℝ l2 )
-            , Differentiable 'Interval ( ℝ l2 )
+            , Differentiable 3 𝕀  ( ℝ l2 )
             )
           => Intersection r1 r2
 intersect
@@ -206,8 +206,8 @@ data Intersection r1 r2 where
     .  ( l12 ~ Length r1r2
        , KnownSymbols r1r2
        , Representable Double ( ℝ l12 )
-       , Differentiable 'Point    ( ℝ l12 )
+       , Differentiable 2 () ( ℝ l12 )
-       , Differentiable 'Interval ( ℝ l12 )
+       , Differentiable 3 𝕀  ( ℝ l12 )
        )
     => { project :: Record r1 -> Record r1r2
           -- ^ project out fields present in both rows
@@ -226,8 +226,8 @@ doIntersection
   => ( forall r1r2 l12.
         ( r1r2 ~ Intersect r1 r2
         , KnownSymbols r1r2, l12 ~ Length r1r2
-        , Differentiable 'Point    ( ℝ l12 )
+        , Differentiable 2 () ( ℝ l12 )
-        , Differentiable 'Interval ( ℝ l12 )
+        , Differentiable 3 𝕀  ( ℝ l12 )
         , Representable Double ( ℝ l12 )
         )
       => Proxy# r1r2 -> Vec l12 ( Fin l1 ) -> Vec l12 ( Fin l2 ) -> kont )

src/splines/Math/Algebra/Dual.hs

@@ -6,7 +6,7 @@
 
 {-# OPTIONS_GHC -Wno-orphans -O2 #-}
 
-{-# OPTIONS_GHC -ddump-splices -ddump-simpl -ddump-to-file -dno-typeable-binds
+{-# OPTIONS_GHC -ddump-splices -ddump-simpl -ddump-cmm -ddump-to-file -dno-typeable-binds
   -dsuppress-unfoldings -dsuppress-coercions #-}
 
 module Math.Algebra.Dual

src/splines/Math/Algebra/Dual/Internal.hs

@@ -178,7 +178,7 @@ data D2𝔸4 v =
   deriving Applicative
     via Generically1 D2𝔸4
 
--- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^3 \).
+-- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^4 \).
 data D3𝔸4 v =
   D34 { _D34_v :: !v
       , _D34_dx, _D34_dy, _D34_dz, _D34_dw :: !( T v )

src/splines/Math/Bezier/Stroke.hs

@@ -1,4 +1,5 @@
 {-# LANGUAGE AllowAmbiguousTypes   #-}
+{-# LANGUAGE PolyKinds             #-}
 {-# LANGUAGE QuantifiedConstraints #-}
 {-# LANGUAGE ScopedTypeVariables   #-}
 {-# LANGUAGE UndecidableInstances  #-}
@@ -13,6 +14,7 @@ module Math.Bezier.Stroke
 
   , brushStroke, envelopeEquation
   , line, bezier2, bezier3
+  , brushStrokeData, pathAndUsedParams
 
   )
   where
@@ -45,9 +47,8 @@ import Data.Maybe
 import Data.Semigroup
   ( sconcat )
 import GHC.Exts
-  ( newMutVar#, runRW#
+  ( newMutVar#, runRW#, inline
   , Proxy#, proxy#
-  , inline
   )
 import GHC.STRef
   ( STRef(..), readSTRef, writeSTRef )
@@ -88,13 +89,6 @@ import Data.Generics.Internal.VL
 import qualified Control.Parallel.Strategies as Strats
   ( rdeepseq, parTuple2, using )
 
--- rounded-hw
-import Numeric.Rounded.Hardware
-  ( Rounded(..) )
-import Numeric.Rounded.Hardware.Interval.NonEmpty
-  ( Interval(..) )
-import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval
-
 -- transformers
 import Control.Monad.Trans.Class
   ( lift )
@@ -123,9 +117,7 @@ import Math.Bezier.Spline
 import qualified Math.Bezier.Quadratic as Quadratic
 import Math.Bezier.Stroke.EnvelopeEquation
 import Math.Differentiable
-  ( Differentiable, DiffInterp
+  ( Differentiable, DiffInterp, I )
-  , type ExtentOrder
-  )
 import Math.Epsilon
   ( epsilon, nearZero )
 import Math.Interval
@@ -236,13 +228,13 @@ computeStrokeOutline ::
      -- Differentiability.
      , Interpolatable Double usedParams
      , Interpolatable ( 𝕀 Double ) ( 𝕀 usedParams )
-     , DiffInterp 'Point    brushParams
+     , DiffInterp 2 () brushParams
-     , DiffInterp 'Interval brushParams
+     , DiffInterp 3 𝕀  brushParams
-     , HasChainRule Double ( ExtentOrder 'Point ) usedParams
+     , HasChainRule Double 2 usedParams
-     , HasChainRule ( 𝕀 Double ) ( ExtentOrder 'Interval ) ( 𝕀 usedParams )
+     , HasChainRule ( 𝕀 Double ) 3 ( 𝕀 usedParams )
-     , HasChainRule Double ( ExtentOrder 'Point ) brushParams
+     , HasChainRule Double 2 brushParams
-     , HasChainRule ( 𝕀 Double ) ( ExtentOrder 'Interval ) ( 𝕀 brushParams )
+     , HasChainRule ( 𝕀 Double ) 3 ( 𝕀 brushParams )
-     , Traversable ( D ( ExtentOrder 'Point ) brushParams )
+     , Traversable ( D 2 brushParams )
 
      -- Debugging.
      , Show ptData, Show brushParams
@@ -251,12 +243,11 @@ computeStrokeOutline ::
   => FitParameters
   -> ( ptData -> usedParams )
   -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing
-  -> ( forall i
+  -> ( forall {t} k (i :: t)
-       .  DiffInterp i brushParams
+       .  DiffInterp k i brushParams
        => Proxy# i
        -> ( forall a. a -> I i a )
-       -> C ( ExtentOrder i )
+       -> C k ( I i brushParams )
-            ( I i brushParams )
               ( Spline Closed () ( I i ( ℝ 2 ) ) )
       )
   -> Spline clo crvData ptData
@@ -381,7 +372,7 @@ computeStrokeOutline fitParams ptParams toBrushParams brushFn spline@( Spline {
           outlineInfo p0 crv :<| go ( openCurveEnd crv ) crvs
 
     brushShape :: ptData -> SplinePts Closed
-    brushShape pt = fun @Double ( brushFn @Point proxy# id ) $ toBrushParams $ ptParams pt
+    brushShape pt = fun @Double ( brushFn @2 @() proxy# id ) $ toBrushParams $ ptParams pt
 
     updateSpline :: ( T ( ℝ 2 ), T ( ℝ 2 ), T ( ℝ 2 ) ) -> ST s OutlineData
     updateSpline ( lastTgt, lastTgtFwd, lastTgtBwd )
@@ -514,25 +505,24 @@ outlineFunction
      -- Differentiability.
      , Interpolatable Double usedParams
      , Interpolatable ( 𝕀 Double ) ( 𝕀 usedParams )
-     , DiffInterp 'Point    brushParams
+     , DiffInterp 2 () brushParams
-     , DiffInterp 'Interval brushParams
+     , DiffInterp 3 𝕀  brushParams
-     , HasChainRule Double ( ExtentOrder 'Point ) usedParams
+     , HasChainRule Double 2 usedParams
-     , HasChainRule ( 𝕀 Double ) ( ExtentOrder 'Interval ) ( 𝕀 usedParams )
+     , HasChainRule ( 𝕀 Double ) 3 ( 𝕀 usedParams )
-     , HasChainRule Double ( ExtentOrder 'Point ) brushParams
+     , HasChainRule Double 2 brushParams
-     , HasChainRule ( 𝕀 Double ) ( ExtentOrder 'Interval ) ( 𝕀 brushParams )
+     , HasChainRule ( 𝕀 Double ) 3 ( 𝕀 brushParams )
-     , Traversable ( D ( ExtentOrder 'Point ) brushParams )
+     , Traversable ( D 2 brushParams )
 
      -- Debugging.
      , Show ptData, Show brushParams
      )
   => ( ptData -> usedParams )
   -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing
-  -> ( forall i
+  -> ( forall {t} k (i :: t)
-       .  DiffInterp i brushParams
+       .  DiffInterp k i brushParams
        => Proxy# i
        -> ( forall a. a -> I i a )
-       -> C ( ExtentOrder i )
+       -> C k ( I i brushParams )
-            ( I i brushParams )
               ( Spline Closed () ( I i ( ℝ 2 ) ) )
       )
   -> ptData
@@ -540,61 +530,36 @@ outlineFunction
   -> OutlineInfo
 outlineFunction ptParams toBrushParams brushFromParams = \ sp0 crv ->
   let
-    pathAndUsedParams :: forall i k arr
-                      .  ( k ~ ExtentOrder i, CurveOrder k
-                         , arr ~ C k
-                         , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 )
-                         , Coercible ( I i ( ℝ 1 ) ) ( I i Double )
-                         , Module ( I i Double ) ( T ( I i ( ℝ 2 ) ) )
-                         , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) )
-                         , Module ( I i Double ) ( T ( I i usedParams ) )
-                         , Torsor ( T ( I i usedParams ) ) ( I i usedParams )
-                         )
-                      => ( forall a. a -> I i a )
-                      -> ( I i ( ℝ 1 ) `arr` I i ( ℝ 2 ), I i ( ℝ 1 ) `arr` I i usedParams )
-    pathAndUsedParams toI =
-      case crv of
-        LineTo    { curveEnd = NextPoint sp1 }
-          | let seg = Segment sp0 sp1
-          -> ( line    @k @i ( fmap ( toI . coords   ) seg )
-             , line    @k @i ( fmap ( toI . ptParams ) seg ) )
-        Bezier2To { controlPoint = sp1, curveEnd = NextPoint sp2 }
-          | let bez2 = Quadratic.Bezier sp0 sp1 sp2
-          -> ( bezier2 @k @i ( fmap ( toI . coords   ) bez2 )
-             , bezier2 @k @i ( fmap ( toI . ptParams ) bez2 ) )
-        Bezier3To { controlPoint1 = sp1, controlPoint2 = sp2, curveEnd = NextPoint sp3 }
-          | let bez3 = Cubic.Bezier sp0 sp1 sp2 sp3
-          -> ( bezier3 @k @i ( fmap ( toI . coords   ) bez3 )
-             , bezier3 @k @i ( fmap ( toI . ptParams ) bez3 ) )
 
-    usedParams :: C ( ExtentOrder 'Point ) ( ℝ 1 ) usedParams
+    usedParams :: C 2 ( ℝ 1 ) usedParams
-    path :: C ( ExtentOrder 'Point ) ( ℝ 1 ) ( ℝ 2 )
+    path :: C 2 ( ℝ 1 ) ( ℝ 2 )
-    ( path, usedParams ) = pathAndUsedParams @Point id
+    ( path, usedParams )
+      = pathAndUsedParams @2 @() id ptParams sp0 crv
 
     curves :: ℝ 1 -- t
-           -> Seq ( ℝ 1 {- s -} -> StrokeDatum Point )
+           -> Seq ( ℝ 1 {- s -} -> StrokeDatum 2 () )
     curves =
-      brushStrokeData @Point @( ExtentOrder 'Point ) @brushParams
+      brushStrokeData @2 @brushParams
         path
-        ( chainRule @Double @( ExtentOrder 'Point )
+        ( chainRule @Double @2
             usedParams
             ( linear toBrushParams )
         )
-        ( brushFromParams @Point proxy# id )
+        ( brushFromParams @2 @() proxy# id )
 
     curvesI :: 𝕀ℝ 1 -- t
-            -> Seq ( 𝕀ℝ 1 {- s -} -> StrokeDatum 'Interval )
+            -> Seq ( 𝕀ℝ 1 {- s -} -> StrokeDatum 3 𝕀 )
-    curvesI = brushStrokeData @'Interval @( ExtentOrder 'Interval ) @brushParams
+    curvesI = brushStrokeData @3 @brushParams
                pathI
-               ( chainRule @( 𝕀 Double ) @( ExtentOrder 'Interval )
+               ( chainRule @( 𝕀 Double ) @3
                    usedParamsI
                    ( linear ( nonDecreasing toBrushParams ) )
                )
-               ( brushFromParams @'Interval proxy# singleton )
+               ( brushFromParams @3 @𝕀 proxy# singleton )
 
-    usedParamsI :: C ( ExtentOrder 'Interval ) ( 𝕀ℝ 1 ) ( 𝕀 usedParams )
+    usedParamsI :: C 3 ( 𝕀ℝ 1 ) ( 𝕀 usedParams )
-    pathI :: C ( ExtentOrder 'Interval ) ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 )
+    pathI :: C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 )
-    ( pathI, usedParamsI ) = pathAndUsedParams @'Interval singleton
+    ( pathI, usedParamsI ) = pathAndUsedParams @3 @𝕀 singleton ptParams sp0 crv
 
     fwdBwd :: OutlineFn
     fwdBwd t
@@ -608,10 +573,10 @@ outlineFunction ptParams toBrushParams brushFromParams = \ sp0 crv ->
         D21 path_t path'_t _ = runD path t
         D21 params_t _ _     = runD usedParams t
         brush_t              = value @Double @2 @brushParams
-                             $ runD ( brushFromParams @Point proxy# id )
+                             $ runD ( brushFromParams @2 @() proxy# id )
                              $ toBrushParams params_t
 
-    ( newtDunno, newtSols ) = intervalNewtonGS InverseMidJacobian 0.0001 curvesI
+    ( newtDunno, newtSols ) = intervalNewtonGS InverseMidJacobian 1e-7 curvesI
 
   in --trace
      --   ( unlines $
@@ -631,6 +596,37 @@ outlineFunction ptParams toBrushParams brushFromParams = \ sp0 crv ->
        }
 {-# INLINEABLE outlineFunction #-}
 
+pathAndUsedParams :: forall k i arr crvData ptData usedParams
+                  .  ( HasType ( ℝ 2 ) ptData
+                     , CurveOrder k
+                     , arr ~ C k
+                     , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 )
+                     , Coercible ( I i ( ℝ 1 ) ) ( I i Double )
+                     , Module ( I i Double ) ( T ( I i ( ℝ 2 ) ) )
+                     , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) )
+                     , Module ( I i Double ) ( T ( I i usedParams ) )
+                     , Torsor ( T ( I i usedParams ) ) ( I i usedParams )
+                     )
+                  => ( forall a. a -> I i a )
+                  -> ( ptData -> usedParams )
+                  -> ptData
+                  -> Curve Open crvData ptData
+                  -> ( I i ( ℝ 1 ) `arr` I i ( ℝ 2 ), I i ( ℝ 1 ) `arr` I i usedParams )
+pathAndUsedParams toI ptParams sp0 crv =
+  case crv of
+    LineTo    { curveEnd = NextPoint sp1 }
+      | let seg = Segment sp0 sp1
+      -> ( line    @k @i ( fmap ( toI . coords   ) seg )
+         , line    @k @i ( fmap ( toI . ptParams ) seg ) )
+    Bezier2To { controlPoint = sp1, curveEnd = NextPoint sp2 }
+      | let bez2 = Quadratic.Bezier sp0 sp1 sp2
+      -> ( bezier2 @k @i ( fmap ( toI . coords   ) bez2 )
+         , bezier2 @k @i ( fmap ( toI . ptParams ) bez2 ) )
+    Bezier3To { controlPoint1 = sp1, controlPoint2 = sp2, curveEnd = NextPoint sp3 }
+      | let bez3 = Cubic.Bezier sp0 sp1 sp2 sp3
+      -> ( bezier3 @k @i ( fmap ( toI . coords   ) bez3 )
+         , bezier3 @k @i ( fmap ( toI . ptParams ) bez3 ) )
+
 -----------------------------------
 -- Various utility functions
 -- used in the "stroke" function.
@@ -924,10 +920,8 @@ withTangent tgt_wanted spline@( Spline { splineStart } )
                 , offset          = T $ Cubic.bezier @( T ( ℝ 2 ) ) bez t
                 }
 
-
+splineCurveFns :: forall k i
-
+               .  ( CurveOrder k
-splineCurveFns :: forall i k
-               .  ( k ~ ExtentOrder i, CurveOrder k
                   , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 )
                   , Module ( I i Double ) ( T ( I i ( ℝ 2 ) ) )
                   , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) )
@@ -958,7 +952,7 @@ solveEnvelopeEquations :: ℝ 1 -- ^ @t@ (for debugging only)
                        -> ℝ 2
                        -> T ( ℝ 2 )
                        -> ( Offset, Offset )
-                       -> Seq ( ℝ 1 -> StrokeDatum Point )
+                       -> Seq ( ℝ 1 -> StrokeDatum 2 () )
                        -> ( ( ℝ 2, T ( ℝ 2 ) ), ( ℝ 2, T ( ℝ 2 ) ) )
 solveEnvelopeEquations _t path_t path'_t ( fwdOffset, bwdOffset ) strokeData
   = ( fwdSol, ( bwdPt, -1 *^ bwdTgt ) )
@@ -1020,7 +1014,7 @@ solveEnvelopeEquations _t path_t path'_t ( fwdOffset, bwdOffset ) strokeData
           | otherwise
           = i + 1
 
-    sol :: Double -> ( ℝ 1 -> StrokeDatum Point ) -> ( Bool, ℝ 1, ℝ 2, T ( ℝ 2 ) )
+    sol :: Double -> ( ℝ 1 -> StrokeDatum 2 () ) -> ( Bool, ℝ 1, ℝ 2, T ( ℝ 2 ) )
     sol initialGuess f =
       let (good, s) = case newtonRaphson maxIters precision domain ( eqn f ) initialGuess of
                         Nothing -> ( False, initialGuess )
@@ -1054,11 +1048,11 @@ solveEnvelopeEquations _t path_t path'_t ( fwdOffset, bwdOffset ) strokeData
                --    ] )
                 ( good, ℝ1 s, value @Double @2 @( ℝ 2 ) dstroke, dcdt )
 
-    eqn :: ( ℝ 1 -> StrokeDatum Point ) -> ( Double -> ( Double, Double ) )
+    eqn :: ( ℝ 1 -> StrokeDatum 2 () ) -> ( Double -> ( Double, Double ) )
     eqn f s =
       case f ( ℝ1 s ) of
-        StrokeDatum { ee = D12 ( ℝ1 ee ) _ ( T ( ℝ1 ee_s ) ) } ->
+        StrokeDatum { ee = D12 ee _ ee_s } ->
-          ( ee, ee_s )
+          coerce ( ee, ee_s )
 
     maxIters :: Word
     maxIters = 5 --30
@@ -1074,10 +1068,9 @@ instance Applicative ZipSeq where
   liftA2 f ( ZipSeq xs ) ( ZipSeq ys ) = ZipSeq ( Seq.zipWith f xs ys )
   {-# INLINE liftA2 #-}
 
-brushStrokeData :: forall i k brushParams arr
+brushStrokeData :: forall k brushParams i arr
-                .  ( k ~ ExtentOrder i, CurveOrder k, arr ~ C k
+                .  ( CurveOrder k, arr ~ C k
-                   , Differentiable i brushParams
+                   , Differentiable k i brushParams
-                   , Fractional ( I i Double )
                    , HasChainRule ( I i Double ) k ( I i ( ℝ 1 ) )
                    , Applicative ( D k ( ℝ 1 ) )
 
@@ -1097,7 +1090,7 @@ brushStrokeData :: forall i k brushParams arr
                      -- ^ brush parameters
                 -> ( I i brushParams `arr` Spline Closed () ( I i ( ℝ 2 ) ) )
                      -- ^ brush from parameters
-                -> ( I i ( ℝ 1 ) -> Seq ( I i ( ℝ 1 ) -> StrokeDatum i ) )
+                -> ( I i ( ℝ 1 ) -> Seq ( I i ( ℝ 1 ) -> StrokeDatum k i ) )
 brushStrokeData path params brush =
   \ t ->
     let
@@ -1108,7 +1101,7 @@ brushStrokeData path params brush =
         dbrush_params :: D k ( I i brushParams ) ( Spline Closed () ( I i ( ℝ 2 ) ) )
         !dbrush_params = runD brush $ value @( I i Double ) @k @( I i ( ℝ 1 ) ) dparams_t
         splines :: Seq ( D k ( I i brushParams ) ( I i ( ℝ 1 ) `arr` I i ( ℝ 2 ) ) )
-        !splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns @i @k ) dbrush_params
+        !splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns @k @i ) dbrush_params
         dbrushes_t :: Seq ( I i ( ℝ 1 ) -> D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) )
         !dbrushes_t = force $ fmap ( uncurryD @k . ( chain @(I i Double) @k dparams_t ) ) splines
           -- This is the crucial use of the chain rule.
@@ -1118,7 +1111,7 @@ brushStrokeData path params brush =
 
     mkStrokeDatum :: D k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) )
                   -> ( I i ( ℝ 1 ) -> D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) )
-                  -> ( I i ( ℝ 1 ) -> StrokeDatum i )
+                  -> ( I i ( ℝ 1 ) -> StrokeDatum k i )
     mkStrokeDatum dpath_t dbrush_t s =
       let dbrush_t_s = dbrush_t s
           dstroke = brushStroke @k dpath_t dbrush_t_s
@@ -1200,7 +1193,7 @@ extendedRecip x@( 𝕀 lo hi )
 
 -- | Computes the brush stroke coordinates of a cusp from
 -- the @(t,s)@ parameter values.
-cuspCoords :: ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 'Point ) )
+cuspCoords :: ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) )
            -> ( 𝕀ℝ 1, Int, 𝕀ℝ 1 )
            -> Cusp
 cuspCoords eqs ( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ), i, 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
@@ -1210,8 +1203,8 @@ cuspCoords eqs ( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ), i, 𝕀 ( ℝ1 s_lo ) ( ℝ1
       <- ( eqs ( ℝ1 t_mid ) `Seq.index` i ) ( ℝ1 s_mid )
   = Cusp
   { cuspParameters   = ℝ2 t_mid s_mid
-  , cuspPathCoords   = dpath
+  , cuspPathCoords   = coerce dpath
-  , cuspStrokeCoords = stroke
+  , cuspStrokeCoords = coerce stroke
   }
   where
     t_mid = 0.5 * ( t_lo + t_hi )
@@ -1228,14 +1221,15 @@ data Preconditioner
 --
 -- Returns @(dunno, sols)@ where @sols@ are boxes that contain a unique solution
 -- (and to which Newton's method will converge starting from anywhere inside
--- the box), and @dunno@ which are small boxes which might or might not
+-- the box), and @dunno@ are small boxes which might or might not
 -- contain solutions.
 intervalNewtonGS :: Preconditioner
                  -> Double
-                 -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 'Interval ) )
+                 -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
                  -> ( [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ], [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] )
 intervalNewtonGS precondMethod minWidth eqs =
-  go [ ( 𝕀 ( ℝ1 0 ) ( ℝ1 1 ), i, 𝕀 ( ℝ1 0 ) ( ℝ1 1 ) )
+  go (0,0)
+    [ ( 𝕀 ( ℝ1 0 ) ( ℝ1 1 ), i, 𝕀 ( ℝ1 0 ) ( ℝ1 1 ) )
     | i <- [ 0 .. length ( eqs ( 𝕀 ( ℝ1 0 ) ( ℝ1 1 ) ) ) - 1 ]
     ]
     []
@@ -1243,18 +1237,25 @@ intervalNewtonGS precondMethod minWidth eqs =
 
   where
 
-    go :: [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] -- boxes to work on
+    go :: ( Int, Int ) -- step counts (for debugging)
+       -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] -- boxes to work on
        -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] -- too small: don't shrink further
        -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] -- found solutions
        -> ( [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ], [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] )
-    go [] giveUp sols = ( giveUp, sols )
+    go ( bis, newt ) [] giveUp sols =
-    go ( cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) )
+      trace ( unlines [ "intervalNewtonGS done"
+                      , "   #bisections: " ++ show bis
+                      , "  #newtonSteps: " ++ show newt
+                      , "         #sols: " ++ show ( length sols )
+                      , "        #dunno: " ++ show ( length giveUp ) ] )
+        ( giveUp, sols )
+    go ( bis, newt ) ( cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) )
                    , i
                    , 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 cands ( cand : giveUp ) sols
+      = go ( bis, newt ) cands ( cand : giveUp ) sols
 
       | StrokeDatum { ee = D22 ee _ _ _ _ _
                     , 𝛿E𝛿sdcdt = D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) }
@@ -1262,14 +1263,14 @@ 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 | Interval.inf ( ival ee ) < Rounded ( ℝ1 0 )
+      = if | inf ee < ℝ1 0
-           , Interval.sup ( ival ee ) > Rounded ( ℝ1 0 )
+           , sup ee > ℝ1 0
-           , cmpℝ2 (<) ( getRounded ( Interval.inf $ ival f ) ) ( ℝ2 0 0 )
+           , cmpℝ2 (<) ( inf f ) ( ℝ2 0 0 )
-           , cmpℝ2 (>) ( getRounded ( Interval.sup $ ival 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
-                                ( f_t_mid, f_s_mid )
+                                ( f_t, f_s )
                                 ( f_t, f_s ) ( neg f_mid )
 
                   !gsGuesses = gaussSeidel a b
@@ -1284,7 +1285,7 @@ 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 ( todo ++ cands ) giveUp ( done ++ sols )
+                   in go ( bis, newt + 1 ) ( todo ++ cands ) giveUp ( done ++ sols )
                  else
                     -- Gauss–Seidel failed to shrink the boxes.
                     -- Bisect along the widest dimension instead.
@@ -1295,11 +1296,11 @@ intervalNewtonGS precondMethod minWidth eqs =
                          | otherwise
                          = [ ( t, i, 𝕀 ( ℝ1 s_lo  ) ( ℝ1 s_mid ) )
                            , ( t, i, 𝕀 ( ℝ1 s_mid ) ( ℝ1 s_hi  ) ) ]
-                   in go ( bisGuesses ++ cands ) giveUp sols
+                   in go ( bis + 1, newt ) ( bisGuesses ++ cands ) giveUp sols
 
            -- Box doesn't contain a solution: discard it.
            | otherwise
-           -> go cands giveUp sols
+           -> go ( bis, newt ) cands giveUp sols
       where
         t_mid = 0.5 * ( t_lo + t_hi )
         s_mid = 0.5 * ( s_lo + s_hi )
@@ -1336,7 +1337,7 @@ precondition meth jac_mid a@( a1, a2 ) b =
             !a22 = 0.5 * ( a22_lo + a22_hi )
             !d = a11 * a22 - a12 * a21
       , not ( nearZero d )
-      , let !precond = ( ℝ2 a22 -a21, ℝ2 -a12 a11 )
+      , let !precond = ( ℝ2 a22 -a12, ℝ2 -a21 a11 )
             !inv = recip d
             f x = scale inv $ matMulVec precond x
       -> ( ( f a1, f a2 ), f b )
@@ -1345,21 +1346,19 @@ precondition meth jac_mid a@( a1, a2 ) b =
 
 scale :: Double -> 𝕀ℝ 2 -> 𝕀ℝ 2
 scale s ( 𝕀 ( ℝ2 a1_lo a2_lo ) ( ℝ2 a1_hi a2_hi ) )
-  | I ( Rounded b1_lo ) ( Rounded b1_hi )
+  | 𝕀 b1_lo b1_hi <- scaleInterval s ( 𝕀 a1_lo a1_hi )
-     <- I ( Rounded s ) ( Rounded s ) * I ( Rounded a1_lo ) ( Rounded a1_hi )
+  , 𝕀 b2_lo b2_hi <- scaleInterval s ( 𝕀 a2_lo a2_hi )
-  , I ( Rounded b2_lo ) ( Rounded b2_hi )
-     <- I ( Rounded s ) ( Rounded s ) * I ( Rounded a2_lo ) ( Rounded a2_hi )
   = 𝕀 ( ℝ2 b1_lo b2_lo ) ( ℝ2 b1_hi b2_hi )
 
 matMulVec :: ( ℝ 2, ℝ 2 ) -> 𝕀ℝ 2 -> 𝕀ℝ 2
-matMulVec ( ℝ2 a11 a21, ℝ2 a12 a22 ) ( 𝕀 ( ℝ2 u_lo v_lo ) ( ℝ2 u_hi v_hi ) ) =
+matMulVec ( ℝ2 a11 a21, ℝ2 a12 a22 ) ( 𝕀 ( ℝ2 u_lo v_lo ) ( ℝ2 u_hi v_hi ) )
-  let !( I ( Rounded u'_lo ) ( Rounded u'_hi ) ) =
+  | 𝕀 u'_lo u'_hi <-
-         I ( Rounded a11 ) ( Rounded a11 ) * I ( Rounded u_lo ) ( Rounded u_hi )
+         𝕀 a11 a11 * 𝕀 u_lo u_hi
-       + I ( Rounded a12 ) ( Rounded a12 ) * I ( Rounded v_lo ) ( Rounded v_hi )
+       + 𝕀 a12 a12 * 𝕀 v_lo v_hi
-      !( I ( Rounded v'_lo ) ( Rounded v'_hi ) ) =
+  , 𝕀 v'_lo v'_hi <-
-         I ( Rounded a21 ) ( Rounded a21 ) * I ( Rounded u_lo ) ( Rounded u_hi )
+         𝕀 a21 a21 * 𝕀 u_lo u_hi
-       + I ( Rounded a22 ) ( Rounded a22 ) * I ( Rounded v_lo ) ( Rounded v_hi )
+       + 𝕀 a22 a22 * 𝕀 v_lo v_hi
-  in 𝕀 ( ℝ2 u'_lo v'_lo ) ( ℝ2 u'_hi v'_hi )
+  = 𝕀 ( ℝ2 u'_lo v'_lo ) ( ℝ2 u'_hi v'_hi )
 
 
 cmpℝ2 :: ( Double -> Double -> Bool ) -> ℝ 2 -> ℝ 2 -> Bool

src/splines/Math/Bezier/Stroke/EnvelopeEquation.hs

@@ -1,4 +1,5 @@
 {-# LANGUAGE AllowAmbiguousTypes #-}
+{-# LANGUAGE PolyKinds           #-}
 {-# LANGUAGE RebindableSyntax    #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 
@@ -11,8 +12,10 @@ module Math.Bezier.Stroke.EnvelopeEquation
 import Prelude hiding ( Num(..), (^) )
 import Data.Coerce
   ( Coercible, coerce )
+import Data.Functor.Identity
+  ( Identity(..) )
 import Data.Kind
-  ( Constraint )
+  ( Type, Constraint )
 import GHC.TypeNats
   ( Nat, type (-) )
 
@@ -25,7 +28,7 @@ import Math.Algebra.Dual
 import qualified Math.Bezier.Cubic     as Cubic
 import qualified Math.Bezier.Quadratic as Quadratic
 import Math.Differentiable
-  ( type ExtentOrder )
+  ( type I )
 import Math.Interval
 import Math.Linear
 import Math.Module
@@ -37,34 +40,35 @@ import Math.Ring
 -- | 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 i
+type StrokeDatum :: Nat -> k -> Type
+data StrokeDatum k i
   = StrokeDatum
   { -- | Path \( p(t_0) \).
-    dpath    :: D ( ExtentOrder i ) ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) )
+    dpath    :: D k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) )
     -- | Brush shape \( b(t_0, s_0) \).
-  , dbrush   :: D ( ExtentOrder i ) ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) )
+  , dbrush   :: D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 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) \).
-  , dstroke  :: D ( ExtentOrder i ) ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) )
+  , dstroke  :: D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) )
 
     -- | Envelope function
     --
     -- \[ E(t_0,s_0) = \left ( \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} \right )_{(t_0,s_0)}. \]
-  , ee       :: D ( ExtentOrder i - 1 ) ( I i ( ℝ 2 ) ) ( I i ( ℝ 1 ) )
+  , ee       :: D ( k - 1 ) ( I i ( ℝ 2 ) ) ( I i ( ℝ 1 ) )
 
    -- \[ \frac{\partial E}{\partial s} \frac{\mathrm{d} c}{\mathrm{d} t}, \]
    --
    -- where \( \frac{\mathrm{d} c}{\mathrm{d} t} \)
    --
    -- denotes a total derivative.
-  , 𝛿E𝛿sdcdt :: D ( ExtentOrder i - 2 ) ( I i ( ℝ 2 ) ) ( T ( I i ( ℝ 2 ) ) )
+  , 𝛿E𝛿sdcdt :: D ( k - 2 ) ( I i ( ℝ 2 ) ) ( T ( I i ( ℝ 2 ) ) )
 
   }
 
-deriving stock instance Show ( StrokeDatum 'Point    )
+deriving stock instance Show ( StrokeDatum 2 () )
-deriving stock instance Show ( StrokeDatum 'Interval )
+deriving stock instance Show ( StrokeDatum 3 𝕀 )
 
 --------------------------------------------------------------------------------
 
@@ -123,8 +127,7 @@ class CurveOrder k where
   --
   -- denotes a total derivative.
   envelopeEquation :: forall i
-                   .  ( k ~ ExtentOrder i
+                   .  ( D ( k - 2 ) ( I i ( ℝ 2 ) ) ~ D ( k - 2 ) ( ℝ 2 )
-                      , D ( k - 2 ) ( I i ( ℝ 2 ) ) ~ D ( k - 2 ) ( ℝ 2 )
                       , D ( k - 1 ) ( I i ( ℝ 2 ) ) ~ D ( k - 1 ) ( ℝ 2 )
                       , D k ( I i ( ℝ 2 ) ) ~ D k ( ℝ 2 )
                       , Cross ( I i Double ) ( T ( I i ( ℝ 2 ) ) )
@@ -134,7 +137,6 @@ class CurveOrder k where
                    -> ( D ( k - 1 ) ( I i ( ℝ 2 ) ) ( I i ( ℝ 1 ) )
                       , D ( k - 2 ) ( I i ( ℝ 2 ) ) ( T ( I i ( ℝ 2 ) ) ) )
 
-
 instance CurveOrder 2 where
   line ( Segment a b :: Segment b ) =
     D \ ( coerce -> t ) ->

src/splines/Math/Differentiable.hs

@@ -1,7 +1,8 @@
+{-# LANGUAGE PolyKinds            #-}
 {-# LANGUAGE UndecidableInstances #-}
 
 module Math.Differentiable
-  ( ExtentOrder, Differentiable, DiffInterp )
+  ( I, Differentiable, DiffInterp )
   where
 
 -- base
@@ -13,52 +14,53 @@ import GHC.TypeNats
 -- MetaBrush
 import Math.Algebra.Dual
   ( D, HasChainRule )
+import Math.Interval
+  ( 𝕀 )
 import Math.Linear
 import Math.Module
-import Math.Interval
-  ( Extent(..), type I )
 import Math.Ring
   ( Transcendental )
 
 --------------------------------------------------------------------------------
 
-type ExtentOrder :: Extent -> Nat
+-- | Type family to parametrise over both pointwise and interval computations.
-type family ExtentOrder e where
+--
-  ExtentOrder 'Point    = 2
+-- Use '()' parameter for points, and '𝕀' parameter for intervals.
-  ExtentOrder 'Interval = 3
+type I :: k -> Type -> Type
--- Currently we're doing order 2 derivatives for the brush stroke fitting,
+type family I i a
--- but order 3 derivatives for the interval Newton method to find cusps.
+type instance I () a = a
+type instance I 𝕀 a = 𝕀 a
 
-type Differentiable :: Extent -> Type -> Constraint
+type Differentiable :: Nat -> k -> Type -> Constraint
 class
   ( Module ( I i Double ) ( T ( I i u ) )
-  , HasChainRule ( I i Double ) ( ExtentOrder i ) ( I i u )
+  , HasChainRule ( I i Double ) k ( I i u )
-  , Traversable ( D ( ExtentOrder i ) ( I i u ) )
+  , Traversable ( D k ( I i u ) )
-  ) => Differentiable i u
+  ) => Differentiable k i u
 instance
   ( Module ( I i Double ) ( T ( I i u ) )
-  , HasChainRule ( I i Double ) ( ExtentOrder i ) ( I i u )
+  , HasChainRule ( I i Double ) k ( I i u )
-  , Traversable ( D ( ExtentOrder i ) ( I i u ) )
+  , Traversable ( D k ( I i u ) )
-  ) => Differentiable i u
+  ) => Differentiable k i u
 
-type DiffInterp :: Extent -> Type -> Constraint
+type DiffInterp :: Nat -> k -> Type -> Constraint
 class
-  ( Differentiable i u
+  ( Differentiable k i u
   , Interpolatable ( I i Double ) ( I i u )
   , Module ( I i Double ) ( T ( I i Double ) )
-  , Module ( D ( ExtentOrder i ) ( I i u ) ( I i Double ) )
+  , Module ( D k ( I i u ) ( I i Double ) )
-           ( D ( ExtentOrder i ) ( I i u ) ( I i ( ℝ 2 ) ) )
+           ( D k ( I i u ) ( I i ( ℝ 2 ) ) )
-  , Transcendental ( D ( ExtentOrder i ) ( I i u ) ( I i Double ) )
+  , Transcendental ( D k ( I i u ) ( I i Double ) )
-  , Applicative ( D ( ExtentOrder i ) ( I i u ) )
+  , Applicative ( D k ( I i u ) )
   , Representable ( I i Double ) ( I i u )
-  ) => DiffInterp i u
+  ) => DiffInterp k i u
 instance
-  ( Differentiable i u
+  ( Differentiable k i u
   , Interpolatable ( I i Double ) ( I i u )
   , Module ( I i Double ) ( T ( I i Double ) )
-  , Module ( D ( ExtentOrder i ) ( I i u ) ( I i Double ) )
+  , Module ( D k ( I i u ) ( I i Double ) )
-           ( D ( ExtentOrder i ) ( I i u ) ( I i ( ℝ 2 ) ) )
+           ( D k ( I i u ) ( I i ( ℝ 2 ) ) )
-  , Transcendental ( D ( ExtentOrder i ) ( I i u ) ( I i Double ) )
+  , Transcendental ( D k ( I i u ) ( I i Double ) )
-  , Applicative ( D ( ExtentOrder i ) ( I i u ) )
+  , Applicative ( D k ( I i u ) )
   , Representable ( I i Double ) ( I i u )
-  ) => DiffInterp i u
+  ) => DiffInterp k i u

src/splines/Math/Interval.hs

@@ -4,21 +4,19 @@
 
 {-# 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
-  , Extent(..), type I
+  , scaleInterval
+  , 𝕀ℝ
   , singleton, nonDecreasing
   )
   where
 
 -- base
 import Prelude hiding ( Num(..) )
-import Data.Coerce
-  ( coerce )
-import Data.Kind
-  ( Type )
-import Data.Monoid
-  ( Sum(..) )
 
 -- acts
 import Data.Act
@@ -37,7 +35,7 @@ import Math.Algebra.Dual
 import Math.Algebra.Dual.Internal
   ( chainRule1NQ )
 import Math.Interval.Internal
-  ( 𝕀(..) )
+  ( 𝕀(𝕀), inf, sup, scaleInterval )
 import Math.Linear
   ( ℝ(..), T(..)
   , RepresentableQ(..)
@@ -52,14 +50,6 @@ import Math.Ring
 type 𝕀ℝ n = 𝕀 ( ℝ n )
 type instance D k ( 𝕀 v ) = D k v
 
--- Handling points and intervals uniformly.
-data Extent = Point | Interval
-
-type I :: Extent -> Type -> Type
-type family I i a where
-  I 'Point    a = a
-  I 'Interval a = 𝕀 a
-
 singleton :: a -> 𝕀 a
 singleton a = 𝕀 a a
 
@@ -76,9 +66,9 @@ deriving via ViaAbelianGroup ( T ( 𝕀 Double ) )
   instance Group ( T ( 𝕀 Double ) )
 
 instance Act   ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where
-  T g • a = coerce ( Sum g • a )
+  T g • a = g + a
 instance Torsor ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where
-  a --> b = T $ getSum ( a --> b )
+  a --> b = T $ b - a
 
 -------------------------------------------------------------------------------
 
@@ -171,9 +161,9 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 1 ) where
                  | Just i <- isLinear mon
                  -> [|| f $$( tabulateQ \ j ->
                                 if | j == i
-                                   -> [|| 1 ||]
+                                   -> [|| 𝕀 1 1 ||]
                                    | otherwise
-                                   -> [|| 0 ||]
+                                   -> [|| 𝕀 0 0 ||]
                             ) ||]
                  | otherwise
                  -> [|| unT o ||]
@@ -206,9 +196,9 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 1 ) where
                  | Just i <- isLinear mon
                  -> [|| f $$( tabulateQ \ j ->
                                 if | j == i
-                                   -> [|| 1 ||]
+                                   -> [|| 𝕀 1 1 ||]
                                    | otherwise
-                                   -> [|| 0 ||]
+                                   -> [|| 𝕀 0 0 ||]
                             ) ||]
                  | otherwise
                  -> [|| unT o ||]
@@ -241,9 +231,9 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 2 ) where
                  | Just i <- isLinear mon
                  -> [|| f $$( tabulateQ \ j ->
                                 if | j == i
-                                   -> [|| 1 ||]
+                                   -> [|| 𝕀 1 1 ||]
                                    | otherwise
-                                   -> [|| 0 ||]
+                                   -> [|| 𝕀 0 0 ||]
                             ) ||]
                  | otherwise
                  -> [|| unT o ||]
@@ -276,9 +266,9 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 2 ) where
                  | Just i <- isLinear mon
                  -> [|| f $$( tabulateQ \ j ->
                                 if | j == i
-                                   -> [|| 1 ||]
+                                   -> [|| 𝕀 1 1 ||]
                                    | otherwise
-                                   -> [|| 0 ||]
+                                   -> [|| 𝕀 0 0 ||]
                             ) ||]
                  | otherwise
                  -> [|| unT o ||]
@@ -311,9 +301,9 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 3 ) where
                  | Just i <- isLinear mon
                  -> [|| f $$( tabulateQ \ j ->
                                 if | j == i
-                                   -> [|| 1 ||]
+                                   -> [|| 𝕀 1 1 ||]
                                    | otherwise
-                                   -> [|| 0 ||]
+                                   -> [|| 𝕀 0 0 ||]
                             ) ||]
                  | otherwise
                  -> [|| unT o ||]
@@ -346,9 +336,9 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 3 ) where
                  | Just i <- isLinear mon
                  -> [|| f $$( tabulateQ \ j ->
                                 if | j == i
-                                   -> [|| 1 ||]
+                                   -> [|| 𝕀 1 1 ||]
                                    | otherwise
-                                   -> [|| 0 ||]
+                                   -> [|| 𝕀 0 0 ||]
                             ) ||]
                  | otherwise
                  -> [|| unT o ||]
@@ -381,9 +371,9 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 4 ) where
                  | Just i <- isLinear mon
                  -> [|| f $$( tabulateQ \ j ->
                                 if | j == i
-                                   -> [|| 1 ||]
+                                   -> [|| 𝕀 1 1 ||]
                                    | otherwise
-                                   -> [|| 0 ||]
+                                   -> [|| 𝕀 0 0 ||]
                             ) ||]
                  | otherwise
                  -> [|| unT o ||]
@@ -416,9 +406,9 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 4 ) where
                  | Just i <- isLinear mon
                  -> [|| f $$( tabulateQ \ j ->
                                 if | j == i
-                                   -> [|| 1 ||]
+                                   -> [|| 𝕀 1 1 ||]
                                    | otherwise
-                                   -> [|| 0 ||]
+                                   -> [|| 𝕀 0 0 ||]
                             ) ||]
                  | otherwise
                  -> [|| unT o ||]

src/splines/Math/Interval/FMA.hs

@@ -0,0 +1,152 @@
+{-# LANGUAGE AllowAmbiguousTypes #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+module Math.Interval.FMA ( addI, subI, prodI, divI, posPowI ) where
+
+-- base
+import Data.Bits
+  ( Bits((.&.), shiftL) )
+import Data.Word
+  ( Word32, Word64 )
+import GHC.Exts
+  ( Double(D#), fmsubDouble#, fnmaddDouble# )
+import GHC.Float
+  ( castFloatToWord32 , castWord32ToFloat
+  , castDoubleToWord64, castWord64ToDouble
+  )
+
+--------------------------------------------------------------------------------
+
+class ( RealFloat f, Num b, Bits b ) => FPBits f b | f -> b, b -> f where
+  toBits   :: f -> b
+  fromBits :: b -> f
+  -- | Size in bytes.
+  sizeOf   :: Int
+
+instance FPBits Float Word32 where
+  toBits   = castFloatToWord32
+  fromBits = castWord32ToFloat
+  sizeOf   = 4
+
+instance FPBits Double Word64 where
+  toBits   = castDoubleToWord64
+  fromBits = castWord64ToDouble
+  sizeOf   = 8
+
+
+{-# SPECIALISE nextAfter :: Float -> Float -> Float #-}
+{-# SPECIALISE nextAfter :: Double -> Double -> Double #-}
+{-# INLINEABLE nextAfter #-}
+-- | @nextAfter a b@ computes the next floating-point value after @a@
+-- in the direction of @b@.
+nextAfter :: forall f b. FPBits f b => f -> f -> f
+nextAfter a b
+  | isNaN a
+  = a
+  | isNaN b
+  = b
+  | a == b
+  = b
+  | otherwise
+  = let !res_bits
+          | a == 0
+          , let !sgn_mask = 1 `shiftL` ( sizeOf @f * 8 - 1 )
+          = ( toBits b .&. sgn_mask ) + 1
+          | ( a < b ) == ( a > 0 )
+          = toBits a + 1
+          | otherwise
+          = toBits a - 1
+    in fromBits res_bits
+
+{-# SPECIALISE succFP :: Float -> Float #-}
+{-# SPECIALISE succFP :: Double -> Double #-}
+{-# INLINEABLE succFP #-}
+-- | The next floating-point number.
+succFP :: forall f b. FPBits f b => f -> f
+succFP x = nextAfter x (1/0)
+
+
+{-# SPECIALISE prevFP :: Float -> Float #-}
+{-# SPECIALISE prevFP :: Double -> Double #-}
+{-# INLINEABLE prevFP #-}
+-- | The previous floating-point number.
+prevFP :: forall f b. FPBits f b => f -> f
+prevFP x = nextAfter x (-1/0)
+
+--------------------------------------------------------------------------------
+
+{-# INLINE withError #-}
+withError :: Double -> Double -> ( Double, Double )
+withError f e =
+  case compare e 0 of
+    LT -> ( prevFP f, f )
+    EQ -> ( f, f )
+    GT -> ( f, succFP f )
+
+addI :: Double -> Double -> ( Double, Double )
+addI a b =
+  let !s  = a + b
+      !a' = s - b
+      !b' = s - a'
+      !da = a - a'
+      !db = b - b'
+      !e = da + db
+  in s `withError` e
+
+subI :: Double -> Double -> ( Double, Double )
+subI a b =
+  let !s  = a - b
+      !a' = s - b
+      !b' = s - a'
+      !da = a - a'
+      !db = b - b'
+      !e = da + db
+  in s `withError` e
+
+{-# INLINE fmsubDouble #-}
+fmsubDouble :: Double -> Double -> Double -> Double
+fmsubDouble ( D# x ) ( D# y ) ( D# z ) = D# ( fmsubDouble# x y z )
+{-# INLINE fnmaddDouble #-}
+fnmaddDouble :: Double -> Double -> Double -> Double
+fnmaddDouble ( D# x ) ( D# y ) ( D# z ) = D# ( fnmaddDouble# x y z )
+
+prodI :: Double -> Double -> ( Double, Double )
+prodI a b =
+  let !p = a * b
+      !e = fmsubDouble a b p
+  in p `withError` e
+
+divI :: Double -> Double -> ( Double, Double )
+divI a b =
+  let !r = a / b
+      !e = fnmaddDouble r b a
+  in r `withError` e
+
+-- | Power of a **non-negative** number to a natural power.
+posPowI :: Double -- ^ Assumed to be non-negative!
+        -> Word -> ( Double, Double )
+posPowI _ 0 = ( 1, 1 )
+posPowI f 1 = ( f, f )
+posPowI f 2 = prodI f f
+posPowI f n
+  | even n
+  , let m = n `quot` 2
+        ( f²_lo, f²_hi ) = prodI f f
+  = ( fst $ posPowI f²_lo m, snd $ posPowI f²_hi m )
+  | otherwise
+  , let m = n `quot` 2
+        ( f²_lo, f²_hi ) = prodI f f
+  = ( fst $ posPowAcc f²_lo m f, snd $ posPowAcc f²_hi m f )
+
+posPowAcc :: Double -> Word -> Double -> ( Double, Double )
+posPowAcc f 1 x = prodI f x
+posPowAcc f n x
+  | even n
+  , let m = n `quot` 2
+        ( f²_lo, f²_hi ) = prodI f f
+  = ( fst $ posPowAcc f²_lo m x, snd $ posPowAcc f²_hi m x )
+  | otherwise
+  , let m = n `quot` 2
+        ( f²_lo, f²_hi ) = prodI f f
+        ( y_lo, y_hi ) = prodI f x
+  = ( fst $ posPowAcc f²_lo m y_lo, snd $ posPowAcc f²_hi m y_hi )

src/splines/Math/Interval/Internal.hs

@@ -1,3 +1,4 @@
+{-# LANGUAGE CPP #-}
 {-# LANGUAGE ScopedTypeVariables  #-}
 {-# LANGUAGE TemplateHaskell      #-}
 {-# LANGUAGE UndecidableInstances #-}
@@ -5,7 +6,7 @@
 {-# OPTIONS_GHC -Wno-orphans #-}
 
 module Math.Interval.Internal
-  ( 𝕀(.., 𝕀) )
+  ( 𝕀(𝕀), inf, sup, scaleInterval )
   where
 
 -- base
@@ -18,17 +19,24 @@ import GHC.Show
 
 -- deepseq
 import Control.DeepSeq
-  ( NFData )
+  ( NFData(..) )
 
 -- rounded-hw
 import Numeric.Rounded.Hardware
   ( Rounded(..) )
-import Numeric.Rounded.Hardware.Interval.NonEmpty
-  ( Interval )
 import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval
   ( Interval(..) )
+
+#ifdef USE_FMA
+-- MetaBrush
+import Math.Interval.FMA
+  ( addI, subI, prodI, divI, posPowI )
+
+#else
+-- rounded-hw
 import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval
-  ( sup, inf, powInt )
+  ( powInt )
+#endif
 
 -- MetaBrush
 import Math.Linear
@@ -42,29 +50,104 @@ import Math.Ring
 --------------------------------------------------------------------------------
 -- Intervals.
 
-newtype 𝕀 a = MkI { ival :: Interval a }
+#ifdef USE_FMA
+data 𝕀 a = 𝕀 { inf, sup :: !a }
+
+instance NFData a => NFData ( 𝕀 a ) where
+  rnf ( 𝕀 lo hi ) = rnf lo `seq` rnf hi
+
+instance Prelude.Num ( 𝕀 Double ) where
+  𝕀 x_lo x_hi + 𝕀 y_lo y_hi
+    | let !z_lo = fst $ addI x_lo y_lo
+          !z_hi = snd $ addI x_hi y_hi
+    = 𝕀 z_lo z_hi
+  𝕀 x_lo x_hi - 𝕀 y_lo y_hi
+    | let !z_lo = fst $ subI x_lo y_hi
+          !z_hi = snd $ subI x_hi y_lo
+    = 𝕀 z_lo z_hi
+  negate (𝕀 lo hi) = 𝕀 (Prelude.negate hi) (Prelude.negate lo)
+  (*) = (*)
+  fromInteger i =
+    let !j = Prelude.fromInteger i
+    in 𝕀 j j
+  abs (𝕀 lo hi)
+    | 0 <= lo
+    = 𝕀 lo hi
+    | hi <= 0
+    = 𝕀 (Prelude.negate hi) (Prelude.negate lo)
+    | otherwise
+    = 𝕀 0 (max (Prelude.negate lo) hi)
+  signum _ = error "No implementation of signum for intervals"
+
+instance Ring ( 𝕀 Double ) where
+  𝕀 lo1 hi1 * 𝕀 lo2 hi2
+    | let !( x_min, x_max ) = prodI lo1 lo2
+          !( y_min, y_max ) = prodI lo1 hi2
+          !( z_min, z_max ) = prodI hi1 lo2
+          !( w_min, w_max ) = prodI hi1 hi2
+    = 𝕀 ( min ( min x_min y_min ) ( min z_min w_min ) )
+        ( max ( max x_max y_max ) ( max z_max w_max ) )
+  _ ^ 0 = 𝕀 1 1
+  iv ^ 1 = iv
+  𝕀 lo hi ^ n
+    | odd n || 0 <= lo
+    , let !lo' = fst $ posPowI lo n
+          !hi' = snd $ posPowI hi n
+    = 𝕀 lo' hi'
+    | hi <= 0
+    , let !lo' = fst $ posPowI (negate hi) n
+          !hi' = snd $ posPowI (negate lo) n
+    = 𝕀 lo' hi'
+    | otherwise
+    , let !hi1 = snd $ posPowI (negate lo) n
+    , let !hi2 = snd $ posPowI hi n
+    = 𝕀 0 ( max hi1 hi2 )
+
+instance Prelude.Fractional ( 𝕀 Double ) where
+  fromRational r =
+    let q = Prelude.fromRational r
+    in 𝕀 q q
+  recip (𝕀 lo hi)
+-- #ifdef ASSERTS
+    | lo >= 0 || hi <= 0
+-- #endif
+    = 𝕀 ( fst $ divI 1 hi ) ( snd $ divI 1 lo )
+-- #ifdef ASSERTS
+    | otherwise
+    = error "BAD interval recip; should use extendedRecip instead"
+-- #endif
+  _ / _ = error "TODO: interval division is not implemented"
+
+instance Transcendental ( 𝕀 Double ) where
+  pi = 𝕀 3.141592653589793 3.1415926535897936
+  cos = withHW Prelude.cos
+  sin = withHW Prelude.sin
+
+{-# INLINE withHW #-}
+withHW :: (Interval.Interval a -> Interval.Interval b) -> 𝕀 a -> 𝕀 b
+withHW f = \ ( 𝕀 lo hi ) ->
+  case f ( Interval.I ( Rounded lo ) ( Rounded hi ) ) of
+    Interval.I ( Rounded x ) ( Rounded y ) -> 𝕀 x y
+
+scaleInterval :: Double -> 𝕀 Double -> 𝕀 Double
+scaleInterval s ( 𝕀 lo hi ) =
+  case compare s 0 of
+    LT -> 𝕀 ( fst $ prodI s hi ) ( snd $ prodI s lo )
+    EQ -> 𝕀 0 0
+    GT -> 𝕀 ( fst $ prodI s lo ) ( snd $ prodI s hi )
+#else
+
+newtype 𝕀 a = MkI { ival :: Interval.Interval a }
   deriving newtype ( Prelude.Num, Prelude.Fractional, Prelude.Floating )
   deriving newtype NFData
 
-instance Eq a => Eq ( 𝕀 a ) where
-  𝕀 a b == 𝕀 c d =
-    a == c && b == d
-
 {-# COMPLETE 𝕀 #-}
 pattern 𝕀 :: a -> a -> 𝕀 a
 pattern 𝕀 x y = MkI ( Interval.I ( Rounded x ) ( Rounded y ) )
-instance Show a => Show ( 𝕀 a ) where
-  showsPrec _ ( 𝕀 x y )
-    = showString "["
-    . showsPrec 0 x
-    . showCommaSpace
-    . showsPrec 0 y
-    . showString "]"
 
-deriving via ViaPrelude ( 𝕀 Double )
+inf, sup :: 𝕀 a -> a
-  instance AbelianGroup  ( 𝕀 Double )
+inf (𝕀 a _) = a
-deriving via ViaPrelude ( 𝕀 Double )
+sup (𝕀 _ b) = b
-  instance AbelianGroup ( T ( 𝕀 Double ) )
 
 instance Ring ( 𝕀 Double ) where
   MkI i1 * MkI i2 = MkI $ i1 Prelude.* i2
@@ -74,16 +157,40 @@ instance Ring ( 𝕀 Double ) where
     -- accidentally use (^) from Prelude.
 
 deriving via ViaPrelude ( 𝕀 Double )
-  instance Field ( 𝕀 Double )
+  instance Transcendental ( 𝕀 Double )
+
+scaleInterval :: Double -> 𝕀 Double -> 𝕀 Double
+scaleInterval s iv = 𝕀 s s * iv -- TODO: could be better
+#endif
 
 deriving via ViaPrelude ( 𝕀 Double )
-  instance Transcendental ( 𝕀 Double )
+  instance AbelianGroup ( 𝕀 Double )
+deriving via ViaPrelude ( 𝕀 Double )
+  instance AbelianGroup ( T ( 𝕀 Double ) )
+deriving via ViaPrelude ( 𝕀 Double )
+  instance Field ( 𝕀 Double )
+
+--------------------------------------------------------------------------------
+
+instance Eq a => Eq ( 𝕀 a ) where
+  𝕀 a b == 𝕀 c d =
+    a == c && b == d
+
+instance Show a => Show ( 𝕀 a ) where
+  showsPrec _ ( 𝕀 x y )
+    = showString "["
+    . showsPrec 0 x
+    . showCommaSpace
+    . showsPrec 0 y
+    . showString "]"
+
+--------------------------------------------------------------------------------
 
 type instance RepDim ( 𝕀 u ) = RepDim u
 instance RepresentableQ r u => RepresentableQ ( 𝕀 r ) ( 𝕀 u ) where
   tabulateQ f =
-    let !lo = tabulateQ @r @u ( \ i -> [|| getRounded $ Interval.inf $ ival $$( f i ) ||] )
+    let !lo = tabulateQ @r @u ( \ i -> [|| inf $ $$( f i ) ||] )
-        !hi = tabulateQ @r @u ( \ i -> [|| getRounded $ Interval.sup $ ival $$( f i ) ||] )
+        !hi = tabulateQ @r @u ( \ i -> [|| sup $ $$( f i ) ||] )
     in [|| 𝕀 $$lo $$hi ||]
 
   indexQ d i =
@@ -95,8 +202,8 @@ instance RepresentableQ r u => RepresentableQ ( 𝕀 r ) ( 𝕀 u ) where
     ||]
 instance Representable r u => Representable ( 𝕀 r ) ( 𝕀 u ) where
   tabulate f =
-    let !lo = tabulate @r @u ( \ i -> getRounded $ Interval.inf $ ival ( f i ) )
+    let !lo = tabulate @r @u ( \ i -> inf $ f i )
-        !hi = tabulate @r @u ( \ i -> getRounded $ Interval.sup $ ival ( f i ) )
+        !hi = tabulate @r @u ( \ i -> sup $ f i )
     in 𝕀 lo hi
 
   index d i =