preparation for interval arithmetic

MetaBrush.cabal

@@ -194,6 +194,8 @@ library splines
           ^>= 3.2.2.0
       , prim-instances
           ^>= 0.2
+      , rounded-hw
+          ^>= 0.3
       , vector
            >= 0.12.1.2   && < 0.14
 

cabal.project

@@ -1,7 +1,9 @@
 packages: .
 
 constraints:
-  acts -finitary
+  acts -finitary,
+  rounded-hw -pure-hs -c99 +avx512 +ghc-prim -x87-long-double
+
 
 -- Fix a severe bug in Waargonaut (no corresponding Hackage release???)
 source-repository-package

src/app/MetaBrush/Action.hs

@@ -1078,7 +1078,7 @@ instance HandleAction Scroll where
                     | dy > 0
                     = max 0.0078125 ( oldZoomFactor / sqrt 2 )
                     | otherwise
-                    = min 256 ( oldZoomFactor * sqrt 2 )
+                    = min 4096 ( oldZoomFactor * sqrt 2 )
                   newCenter :: ℝ 2
                   newCenter
                     = ( 1 - oldZoomFactor / newZoomFactor ) *^ ( oldCenter --> mousePos :: T ( ℝ 2 ) )

src/app/MetaBrush/Application.hs

@@ -200,8 +200,8 @@ runApplication application = do
   maxHistorySizeTVar   <- STM.newTVarIO @Int                                     1000
   fitParametersTVar    <- STM.newTVarIO @FitParameters
                             ( FitParameters
-                              { maxSubdiv  = 2 --3 -- 6
+                              { maxSubdiv  = 3 --2 --3 -- 6
-                              , nbSegments = 3 --6 -- 12
+                              , nbSegments = 40 --3 --6 -- 12
                               , dist_tol   = 0.1 -- 5e-3
                               , t_tol      = 0.1 -- 1e-4
                               , maxIters   = 2   -- 100

src/app/MetaBrush/UI/InfoBar.hs

@@ -172,7 +172,7 @@ updateInfoBar viewportDrawingArea ( InfoBar {..} ) ( Variables { mousePosTVar }
         GTK.labelSetText zoomText $ Text.pack ( fixed 5 2 ( 100 * zoomFactor ) <> "%" )
         case mbMousePos of
           Just ( ℝ2 mx my ) ->
-            GTK.labelSetText cursorPosText $ Text.pack ( "x: " <> fixed 6 2 mx <> "\ny: " <> fixed 6 2 my )
+            GTK.labelSetText cursorPosText $ Text.pack ( "x: " <> fixed 6 4 mx <> "\ny: " <> fixed 6 4 my )
           Nothing ->
             GTK.labelSetText cursorPosText $ Text.pack ( "x: " <> na <> "\ny: " <> na )
         GTK.labelSetText  topLeftPosText $ Text.pack ( "x: " <> fixed 6 2  l <> "\ny: " <> fixed 6 2  t )

src/metabrushes/MetaBrush/Asset/Brushes.hs

@@ -93,7 +93,7 @@ circleBrush =
     e_x = pure $ ℝ2 1 0
     e_y = pure $ ℝ2 0 1
 
-    kon = konst @( Record CircleBrushFields )
+    kon = konst @Double @( Record CircleBrushFields )
 
 ellipseBrush :: Record EllipseBrushFields ~> Spline 'Closed () ( ℝ 2 )
 ellipseBrush =
@@ -112,4 +112,4 @@ ellipseBrush =
     e_x = pure $ ℝ2 1 0
     e_y = pure $ ℝ2 0 1
 
-    kon = konst @( Record EllipseBrushFields )
+    kon = konst @Double @( Record EllipseBrushFields )

src/metabrushes/MetaBrush/Brush.hs

@@ -62,8 +62,8 @@ data Brush brushFields where
   BrushData
     :: forall brushFields
     .  ( KnownSymbols brushFields
-       , Representable ( ℝ ( Length brushFields) )
+       , Representable Double ( ℝ ( Length brushFields) )
-       , Diffy ( ℝ ( Length brushFields) )
+       , Diffy Double ( ℝ ( Length brushFields) )
        , Typeable brushFields )
     => { brushName     :: !Text
        , brushFunction :: BrushFunction brushFields
@@ -98,7 +98,7 @@ class ( KnownSymbols pointFields, Typeable pointFields
       , Show           ( Record pointFields )
       , NFData         ( Record pointFields )
       , Interpolatable ( Record pointFields )
-      , Representable  ( ℝ ( Length pointFields ) )
+      , Representable  Double ( ℝ ( Length pointFields ) )
       )
    => PointFields pointFields where { }
 instance ( KnownSymbols pointFields, Typeable pointFields
@@ -106,7 +106,7 @@ instance ( KnownSymbols pointFields, Typeable pointFields
          , Show           ( Record pointFields )
          , NFData         ( Record pointFields )
          , Interpolatable ( Record pointFields )
-         , Representable  ( ℝ ( Length pointFields ) )
+         , Representable  Double ( ℝ ( Length pointFields ) )
          )
       => PointFields pointFields where { }
 

src/metabrushes/MetaBrush/Records.hs

@@ -80,7 +80,7 @@ deriving newtype
         => NFData ( Record ks )
 
 -- | Show a record, using the given type-level field names.
-instance ( KnownSymbols ks, Representable ( ℝ ( Length ks ) ) )
+instance ( KnownSymbols ks, Representable Double ( ℝ ( Length ks ) ) )
        => Show ( Record ks ) where
   showsPrec p ( MkR r )
     = showParen ( p >= 11 )
@@ -117,11 +117,11 @@ instance ( Torsor ( T ( ℝ ( Length ks ) ) ) ( ℝ ( Length ks ) )
   MkR g --> MkR a = T $ MkR $ unT $ g --> a
 
 deriving newtype
-  instance Representable ( ℝ ( Length ks ) )
+  instance Representable r ( ℝ ( Length ks ) )
-        => Representable ( Record ks )
+        => Representable r ( Record ks )
 
 type instance D ( Record ks ) = D ( ℝ ( Length ks ) )
-deriving newtype instance Diffy ( ℝ ( Length ks ) ) => Diffy ( Record ks )
+deriving newtype instance Diffy Double ( ℝ ( Length ks ) ) => Diffy Double ( Record ks )
 
 --------------------------------------------------------------------------------
 
@@ -148,8 +148,8 @@ intersect :: forall r1 r2 l1 l2
           . ( Typeable r1, Typeable r2
             , KnownSymbols r1, KnownSymbols r2
             , l1 ~ Length r1, l2 ~ Length r2
-            , Representable ( ℝ l1 ), Representable ( ℝ l2 )
+            , Representable Double ( ℝ l1 ), Representable Double ( ℝ l2 )
-            , Interpolatable ( Record r1 ), Diffy ( ℝ l2 )
+            , Interpolatable ( Record r1 ), Diffy Double ( ℝ l2 )
             )
           => Intersection r1 r2
 intersect
@@ -172,8 +172,8 @@ data Intersection r1 r2 where
     :: forall r1r2 r1 r2 l12
     .  ( l12 ~ Length r1r2
        , KnownSymbols r1r2
-       , Representable ( ℝ l12 )
+       , Representable Double ( ℝ l12 )
-       , Diffy ( ℝ l12 )
+       , Diffy Double ( ℝ l12 )
        , Interpolatable ( Record r1r2 ) )
     => { project :: Record r1 -> Record r1r2
        , inject  :: Record r2 -> Record r1r2 ~> Record r2
@@ -185,12 +185,13 @@ doIntersection
   :: forall r1 r2 l1 l2 kont
   .  ( KnownSymbols r1, KnownSymbols r2
      , l1 ~ Length r1, l2 ~ Length r2
-     , Representable ( ℝ l1 ), Representable ( ℝ l2 )
+     , Representable Double ( ℝ l1 ), Representable Double ( ℝ l2 )
      )
   => ( forall r1r2 l12.
         ( r1r2 ~ Intersect r1 r2, l12 ~ Length r1r2
-        , Representable ( ℝ l12 ), Diffy ( ℝ l12 ), Interpolatable ( ℝ l12 )
+        , Representable Double ( ℝ l12 ), Diffy Double ( ℝ l12 )
-        , KnownSymbols r1r2, Representable ( ℝ ( Length r1r2 ) )
+        , Interpolatable ( ℝ l12 )
+        , KnownSymbols r1r2, Representable Double ( ℝ ( Length r1r2 ) )
         )
       => Proxy# r1r2 -> Vec l12 ( Fin l1 ) -> Vec l12 ( Fin l2 ) -> kont )
   -> kont

src/metabrushes/MetaBrush/Serialisable.hs

@@ -110,7 +110,8 @@ instance Serialisable ( T ( ℝ 2 ) ) where
       JSON.Encoder.atKey' "x" encoder x
     . JSON.Encoder.atKey' "y" encoder y
   decoder = V2 <$> JSON.Decoder.atKey "x" decoder <*> JSON.Decoder.atKey "y" decoder
-instance ( KnownSymbols ks, Representable ( ℝ ( Length ks ) ) ) => Serialisable ( Record ks ) where
+instance ( KnownSymbols ks, Representable Double ( ℝ ( Length ks ) ) )
+      => Serialisable ( Record ks ) where
   encoder = contramap encodeFields ( JSON.Encoder.keyValueTupleFoldable ( encoder @Double ) )
     where
       encodeFields :: Record ks -> [ ( Text, Double ) ]

src/splines/Math/Bezier/Stroke.hs

@@ -2,6 +2,8 @@
 {-# LANGUAGE QuantifiedConstraints #-}
 {-# LANGUAGE ScopedTypeVariables   #-}
 
+{-# LANGUAGE DuplicateRecordFields  #-}
+
 module Math.Bezier.Stroke
   ( Offset(..)
   , CachedStroke(..), discardCache, invalidateCache
@@ -106,9 +108,8 @@ import qualified Math.Bezier.Quadratic as Quadratic
 import Math.Epsilon
   ( epsilon )
 import Math.Module
-  ( Module(..), Inner((^.^)), Interpolatable
+  ( Module(..), Inner((^.^)), Cross(cross), Interpolatable
-  , lerp, cross
+  , lerp, convexCombination, strictlyParallel
-  , convexCombination, strictlyParallel
   )
 import Math.Orientation
   ( Orientation(..), splineOrientation
@@ -181,7 +182,7 @@ computeStrokeOutline ::
   forall ( clo :: SplineType ) usedParams brushParams crvData ptData s
   .  ( KnownSplineType clo
      , Interpolatable usedParams
-     , Diffy usedParams, Diffy brushParams
+     , Diffy Double usedParams, Diffy Double brushParams
      , HasType ( ℝ 2 ) ptData
      , HasType ( CachedStroke s ) crvData
      , NFData ptData, NFData crvData
@@ -414,7 +415,7 @@ computeStrokeOutline fitParams ptParams toBrushParams brushFn spline@( Spline {
 outlineFunction
   :: forall usedParams brushParams crvData ptData
   .  ( Interpolatable usedParams
-     , Diffy usedParams, Diffy brushParams
+     , Diffy Double usedParams, Diffy Double brushParams
      , HasType ( ℝ 2 ) ptData
      -- Debugging.
      , Show ptData, Show brushParams
@@ -456,8 +457,9 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv =
 
         D1 path_t path'_t _ = runD path t
         D1 params_t _ _     = runD usedParams t
-        brush_t             = value @brushParams
+        brush_t             = value @Double @brushParams
-                            $ runD brushFromParams ( fun toBrushParams params_t )
+                            $ runD brushFromParams
+                            $ fun toBrushParams params_t
 
   in fwdBwd
 
@@ -750,9 +752,10 @@ withTangent tgt_wanted spline@( Spline { splineStart } )
 --
 --   - \( p(t) \) is the path that the brush follows, and
 --   - \( b(t,s) \) is the brush shape, as it varies along the path.
-brushStroke :: D ( ℝ 1 ) ( ℝ 2 ) -- ^ stroke path \( p(t) \)
+brushStroke :: Module r ( T v )
-            -> D ( ℝ 2 ) ( ℝ 2 ) -- ^ brush \( b(t,s) \)
+            => D ( ℝ 1 ) v -- ^ stroke path \( p(t) \)
-            -> D ( ℝ 2 ) ( ℝ 2 )
+            -> D ( ℝ 2 ) v -- ^ brush \( b(t,s) \)
+            -> D ( ℝ 2 ) v
 brushStroke ( D1 p dpdt d2pdt2 ) ( D2 b dbdt dbds d2bdt2 d2bdtds d2bds2 ) =
   D2 ( unT $ T p ^+^ T b )
      -- c = p + b
@@ -770,27 +773,33 @@ brushStroke ( D1 p dpdt d2pdt2 ) ( D2 b dbdt dbds d2bdt2 d2bdtds d2bds2 ) =
 --
 -- \[ E = \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} = 0, ]
 --
--- as well as the vector
+-- together with the total derivative
 --
--- \[ \frac{\partial E}{\partial s} \frac{\mathrm{d} c}{\mathrm{d} t} \]
+-- \[ \frac{\mathrm{d} c}{\mathrm{d} t}, \]
 --
--- whose roots correspond to cusps in the envelope, and
+-- and the partial derivatives
 --
--- \[ \frac{\partial E}{\partial s}. \]
+-- \[ \frac{\partial E}{\partial s}, \qquad \frac{\partial E}{\partial s}. \]
-envelopeEquation :: D ( ℝ 2 ) ( ℝ 2 ) -> ( Double, T ( ℝ 2 ), Double, Double )
+--
+-- NB: if \( \frac{\partial E}{\partial s} \) is zero, the total derivative is ill-defined.
+envelopeEquation :: ( D ( i ( ℝ 2 ) ) ~ D ( ℝ 2 )
+                    , Cross ( i Double ) ( T ( i ( ℝ 2 ) ) )
+                    , Fractional ( i Double )
+                    )
+                 => D ( i ( ℝ 2 ) ) ( i ( ℝ 2 ) )
+                 -> ( i Double, T ( i ( ℝ 2 ) ), i Double, i Double )
 envelopeEquation ( D2 _ dcdt dcds d2cdt2 d2cdtds d2cds2 ) =
   let ee = dcdt `cross` dcds
       dEdt = d2cdt2  `cross` dcds + dcdt `cross` d2cdtds
       dEds = d2cdtds `cross` dcds + dcdt `cross` d2cds2
-      tot = dcdt ^-^ ( dEdt / dEds ) *^ dcds --dEds *^ dcdt ^-^ dEdt *^ dcds
+      tot = dcdt ^-^ ( dEdt / dEds ) *^ dcds
   in ( ee, tot, dEdt, dEds )
   -- Computation of total derivative dc/dt:
   --
   -- dc/dt = ∂c/∂t + ∂c/∂s ∂s/∂t
-  -- ∂s/∂t = - ( ∂E / ∂t ) / ( ∂E / ∂s )
+  -- ∂s/∂t = - ∂E/∂t / ∂E/∂s
   --
-  -- ( ∂E / ∂s ) dc/dt = ( ∂E / ∂s ) ∂c/∂t - ( ∂E / ∂t ) ∂c/∂s.
+  -- ∂E/∂s dc/dt = ∂E/∂s ∂c/∂t - ∂E/∂t ∂c/∂s.
-
 
 -- | Linear interpolation, as a differentiable function.
 line :: forall b. ( Module Double ( T b ), Torsor ( T b ) b )
@@ -906,9 +915,27 @@ solveEnvelopeEquations path_t path'_t ( fwdOffset, bwdOffset ) strokeData
                         Nothing -> ( False, initialGuess )
                         Just s0 -> ( True , s0 )
       in case f ( ℝ1 s ) of -- TODO: a bit redundant to have to compute this again...
-          StrokeDatum { dstroke, dcdt } ->
+          StrokeDatum { ee = _ee, dstroke, 𝛿E𝛿t = _𝛿E𝛿t, 𝛿E𝛿s, dcdt } ->
-            ( good, ℝ1 s, value @( ℝ 2 ) dstroke, dcdt )
+            -- The total derivative dc/dt is computed by dividing by ∂E/∂s,
-
+            -- so check it isn't zero first. This corresponds to cusps in the envelope.
+            let totDeriv
+                  | abs 𝛿E𝛿s < epsilon
+                  , let s' = if s >= 0.5 then s - 1e-9 else s + 1e-9
+                  = case f ( ℝ1 s' ) of { StrokeDatum { dcdt = dcdt_s' } -> dcdt_s' }
+                  | otherwise
+                  = dcdt
+            in --trace
+               --  ( unlines
+               --    [ "solveEnvelopeEquations"
+               --    , "      t = " ++ show t
+               --    , "      s = " ++ show s
+               --    , "      c = " ++ show dstroke
+               --    , "      E = " ++ show _ee
+               --    , "  ∂E/∂t = " ++ show _𝛿E𝛿t
+               --    , "  ∂E/∂s = " ++ show 𝛿E𝛿s
+               --    , "  dc/dt = " ++ show totDeriv
+               --    ] )
+               ( good, ℝ1 s, value @Double @( ℝ 2 ) dstroke, totDeriv )
 
     eqn :: ( ℝ 1 -> StrokeDatum ) -> ( Double -> ( Double, Double ) )
     eqn f s =
@@ -930,7 +957,7 @@ instance Applicative ZipSeq where
   liftA2 f ( ZipSeq xs ) ( ZipSeq ys ) = ZipSeq ( Seq.zipWith f xs ys )
 
 brushStrokeData :: forall brushParams
-                .  ( Diffy brushParams, Show brushParams )
+                .  ( Diffy Double brushParams, Show brushParams )
                 => ( ℝ 1 ~> ℝ 2 )                      -- ^ path
                 -> ( ℝ 1 ~> brushParams )              -- ^ brush parameters
                 -> ( brushParams ~> SplinePts Closed ) -- ^ brush from parameters
@@ -947,7 +974,7 @@ brushStrokeData path params brush =
         splines :: Seq ( D brushParams ( ℝ 1 ~> ℝ 2 ) )
         !splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns ) dbrush_params
         dbrushes_t :: Seq ( ℝ 1 -> D ( ℝ 2 ) ( ℝ 2 ) )
-        !dbrushes_t = force $ fmap ( uncurryD' . ( dparams_t `chain` ) ) splines
+        !dbrushes_t = force $ fmap ( uncurryD . ( dparams_t `chain` ) ) splines
 
     in fmap ( mkStrokeDatum dpath_t ) dbrushes_t
   where
@@ -958,10 +985,11 @@ brushStrokeData path params brush =
     mkStrokeDatum dpath_t dbrush_t s =
       let dbrush_t_s = dbrush_t s
           dstroke@( D2 _c _𝛿c𝛿t _𝛿c𝛿s _ _ _ ) = brushStroke dpath_t dbrush_t_s
-          ( ee, dcdt, _𝛿E𝛿t, 𝛿E𝛿s ) = envelopeEquation dstroke
+          ( ee, dcdt, 𝛿E𝛿t, 𝛿E𝛿s ) = coerce $ envelopeEquation @Identity $ coerce dstroke
       in  -- trace
           --   ( unlines
           --     [ "envelopeEquation:"
+          --     , "      t = " ++ show t
           --     , "      s = " ++ show s
           --     , "      c = " ++ show _c
           --     , "  ∂c/∂t = " ++ show _𝛿c𝛿t
@@ -974,7 +1002,7 @@ brushStrokeData path params brush =
             { dpath  = dpath_t
             , dbrush = dbrush_t_s
             , dstroke
-            , ee, dcdt, 𝛿E𝛿s }
+            , ee, dcdt, 𝛿E𝛿t, 𝛿E𝛿s }
 
 -- | The value and derivative of a brush stroke at a given coordinate
 -- \( (t_0, s_0) \), together with the value of the envelope equation at that
@@ -994,11 +1022,16 @@ data StrokeDatum
     --
     -- \[ E(t_0,s_0) = \left ( \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} \right )_{(t_0,s_0)}. \]
   , ee       :: Double
+    -- | \( \left ( \frac{\partial E}{\partial s} \right )_{(t_0,s_0)}. \)
+  , 𝛿E𝛿s     :: Double
+    -- | \( \left ( \frac{\partial E}{\partial t} \right )_{(t_0,s_0)}. \)
+  , 𝛿E𝛿t     :: Double
     -- | Total derivative
     --
     -- \[ \left ( \frac{\mathrm{d} c}{\mathrm{d} t} \right )_{(t_0,s_0)}. \]
-  , dcdt     :: T ( ℝ 2 )
+    --
-    -- | \( \left ( \frac{\partial E}{\partial s} \right )_{(t_0,s_0)}. \)
+    -- This is ill-defined when \( \frac{\partial E}{\partial s} = 0 \).
-  , 𝛿E𝛿s     :: Double
+  , dcdt :: T ( ℝ 2 )
+
   }
   deriving stock Show

src/splines/Math/Linear.hs

@@ -12,6 +12,9 @@ module Math.Linear
   , Fin(..), eqFin, MFin(..)
   , Dim, Representable(..), injection, projection
   , Vec(..), (!), find
+
+  -- * Intervals
+  , 𝕀, 𝕀ℝ
   ) where
 
 -- base
@@ -44,6 +47,10 @@ import Data.Group
 import Data.Group.Generics
   ( )
 
+-- rounded-hw
+import Numeric.Rounded.Hardware.Interval.NonEmpty
+  ( Interval )
+
 --------------------------------------------------------------------------------
 
 data Mat22 = Mat22 !Double !Double !Double !Double
@@ -171,24 +178,24 @@ type family Dim v
 
 type instance Dim ( ℝ n ) = n
 
-type Representable :: Type -> Constraint
+type Representable :: Type -> Type -> Constraint
-class Representable v where
+class Representable r v | v -> r where
-  tabulate :: ( Fin ( Dim v ) -> Double ) -> v
+  tabulate :: ( Fin ( Dim v ) -> r ) -> v
-  index    :: v -> Fin ( Dim v ) -> Double
+  index    :: v -> Fin ( Dim v ) -> r
 
-instance Representable ( ℝ 0 ) where
+instance Representable Double ( ℝ 0 ) where
   {-# INLINE tabulate #-}
   tabulate _ = ℝ0
   {-# INLINE index #-}
   index _ _ = 0
 
-instance Representable ( ℝ 1 ) where
+instance Representable Double ( ℝ 1 ) where
   {-# INLINE tabulate #-}
   tabulate f = ℝ1 ( f ( Fin 1## ) )
   {-# INLINE index #-}
   index ( ℝ1 x ) _ = x
 
-instance Representable ( ℝ 2 ) where
+instance Representable Double ( ℝ 2 ) where
   {-# INLINE tabulate #-}
   tabulate f = ℝ2 ( f ( Fin 1## ) ) ( f ( Fin 2## ) )
   {-# INLINE index #-}
@@ -196,7 +203,7 @@ instance Representable ( ℝ 2 ) where
     Fin 1## -> x
     _       -> y
 
-instance Representable ( ℝ 3 ) where
+instance Representable Double ( ℝ 3 ) where
   {-# INLINE tabulate #-}
   tabulate f = ℝ3 ( f ( Fin 1## ) ) ( f ( Fin 2## ) ) ( f ( Fin 3## ) )
   {-# INLINE index #-}
@@ -206,14 +213,14 @@ instance Representable ( ℝ 3 ) where
     _       -> z
 
 {-# INLINE projection #-}
-projection :: ( Representable u, Representable v )
+projection :: ( Representable r u, Representable r v )
            => ( Fin ( Dim v ) -> Fin ( Dim u ) )
            -> u -> v
 projection f = \ u ->
   tabulate \ i -> index u ( f i )
 
 {-# INLINE injection #-}
-injection :: ( Representable u, Representable v )
+injection :: ( Representable r u, Representable r v )
           => ( Fin ( Dim v ) -> MFin ( Dim u ) )
           -> u -> v -> v
 injection f = \ u v ->
@@ -242,3 +249,9 @@ find eq v b = MFin ( go 1## v )
       | otherwise
       = go ( j `plusWord#` 1## ) as
     go _ VZ = 0##
+
+--------------------------------------------------------------------------------
+-- Intervals.
+
+type 𝕀 = Interval
+type 𝕀ℝ n = 𝕀 ( ℝ n )

src/splines/Math/Linear/Dual.hs

@@ -11,6 +11,8 @@ import Control.Applicative
   ( liftA2 )
 import Data.Coerce
   ( coerce )
+import Data.Functor.Identity
+  ( Identity(..) )
 import Data.Kind
   ( Type, Constraint )
 import GHC.Generics
@@ -28,11 +30,12 @@ type (~>) :: Type -> Type -> Type
 newtype u ~> v = D { runD :: u -> D u v }
 deriving stock instance Functor ( D u ) => Functor ( (~>) u )
 
-instance ( Applicative ( D u ), Module ( D u Double ) ( D u v ) ) => Module Double ( T ( u ~> v ) ) where
+instance ( Applicative ( D u ), Module r ( T v ) )
-  origin = T $ D \ _ -> origin
+       => Module r ( T ( u ~> v ) ) where
-  T ( D f ) ^+^ T ( D g ) = T $ D \ t -> f t ^+^ g t
+  origin = T $ D \ _ -> pure $ coerce $ origin @r @( T v )
-  T ( D f ) ^-^ T ( D g ) = T $ D \ t -> f t ^-^ g t
+  T ( D f ) ^+^ T ( D g ) = T $ D \ t -> liftA2 ( coerce $ (^+^) @r @( T v ) ) ( f t ) ( g t )
-  a *^ T ( D f ) = T $ D \ t -> pure a *^ f t
+  T ( D f ) ^-^ T ( D g ) = T $ D \ t -> liftA2 ( coerce $ (^-^) @r @( T v ) ) ( f t ) ( g t )
+  a *^ T ( D f ) = T $ D \ t -> fmap ( coerce $ (*^) @r @( T v ) a ) $ f t
 
 -- | @D ( ℝ n ) v@ is \( \mathbb{R}[x_1, \ldots, x_n]/(x_1, \ldots, x_n)^3 \otimes_\mathbb{R} v \)
 type D :: Type -> Type -> Type
@@ -42,6 +45,9 @@ type instance D ( ℝ 1 ) = Dℝ1
 type instance D ( ℝ 2 ) = Dℝ2
 type instance D ( ℝ 3 ) = Dℝ3
 
+type instance D ( Identity a ) = D a
+type instance D ( 𝕀 u ) = D u
+
 newtype Dℝ0 v = D0 { v :: v }
   deriving stock   ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
   deriving newtype ( Num, Fractional, Floating )
@@ -67,7 +73,7 @@ instance Num ( Dℝ1 Double ) where
   (+) = liftA2 (+)
   (-) = liftA2 (-)
   negate = fmap negate
-  fromInteger = konst @( ℝ 1 ) . fromInteger
+  fromInteger = konst @Double @( ℝ 1 ) . fromInteger
 
   abs    = error "no"
   signum = error "no"
@@ -82,7 +88,7 @@ instance Num ( Dℝ2 Double ) where
   (+) = liftA2 (+)
   (-) = liftA2 (-)
   negate = fmap negate
-  fromInteger = konst @( ℝ 2 ) . fromInteger
+  fromInteger = konst @Double @( ℝ 2 ) . fromInteger
 
   abs    = error "no"
   signum = error "no"
@@ -101,7 +107,7 @@ instance Num ( Dℝ3 Double ) where
   (+) = liftA2 (+)
   (-) = liftA2 (-)
   negate = fmap negate
-  fromInteger = konst @( ℝ 3 ) . fromInteger
+  fromInteger = konst @Double @( ℝ 3 ) . fromInteger
 
   abs    = error "no"
   signum = error "no"
@@ -146,9 +152,9 @@ instance Module Double ( T v ) => Module ( Dℝ3 Double ) ( Dℝ3 v ) where
 
 instance Fractional ( Dℝ1 Double ) where
   (/) = error "I haven't yet defined (/) for Dℝ1"
-  fromRational = konst @( ℝ 1 ) . fromRational
+  fromRational = konst @Double @( ℝ 1 ) . fromRational
 instance Floating ( Dℝ1 Double ) where
-  pi = konst @( ℝ 1 ) pi
+  pi = konst @Double @( ℝ 1 ) pi
   sin ( D1 v ( T dx ) ( T ddx ) )
     = let !s = sin v
           !c = cos v
@@ -161,9 +167,9 @@ instance Floating ( Dℝ1 Double ) where
 
 instance Fractional ( Dℝ2 Double ) where
   (/) = error "I haven't yet defined (/) for Dℝ2"
-  fromRational = konst @( ℝ 2 ) . fromRational
+  fromRational = konst @Double @( ℝ 2 ) . fromRational
 instance Floating ( Dℝ2 Double ) where
-  pi = konst @( ℝ 2 ) pi
+  pi = konst @Double @( ℝ 2 ) pi
   sin ( D2 v ( T dx ) ( T dy ) ( T ddx ) ( T dxdy ) ( T ddy ) )
     = let !s = sin v
           !c = cos v
@@ -184,9 +190,9 @@ instance Floating ( Dℝ2 Double ) where
 
 instance Fractional ( Dℝ3 Double ) where
   (/) = error "I haven't yet defined (/) for Dℝ3"
-  fromRational = konst @( ℝ 3 ) . fromRational
+  fromRational = konst @Double @( ℝ 3 ) . fromRational
 instance Floating ( Dℝ3 Double ) where
-  pi = konst @( ℝ 3 ) pi
+  pi = konst @Double @( ℝ 3 ) pi
   sin ( D3 v ( T dx ) ( T dy ) ( T dz ) ( T ddx ) ( T dxdy ) ( T ddy ) ( T dxdz ) ( T dydz ) ( T ddz ) )
     = let !s = sin v
           !c = cos v
@@ -213,30 +219,24 @@ instance Floating ( Dℝ3 Double ) where
 
 --------------------------------------------------------------------------------
 
-uncurryD :: ( ℝ 1 ~> ℝ 1 ~> b ) -> ( ℝ 2 ~> b )
+uncurryD :: D a ~ D ( ℝ 1 )
-uncurryD ( D b ) = D \ ( ℝ2 t0 s0 ) -> uncurryD' ( b ( ℝ1 t0 ) ) ( ℝ1 s0 )
+         => D ( ℝ 1 ) ( a ~> b ) -> a -> D ( ℝ 2 ) b
-
+uncurryD ( D1 ( D b_t0 ) ( T ( D dbdt_t0 ) ) ( T ( D d2bdt2_t0 ) ) ) s0 =
-uncurryD' :: D ( ℝ 1 ) ( ℝ 1 ~> b ) -> ℝ 1 -> D ( ℝ 2 ) b
+  let !( D1 b_t0s0 dbds_t0s0 d2bds2_t0s0 ) = b_t0 s0
-uncurryD' ( D1 ( D b_t0 ) ( T ( D dbdt_t0 ) ) ( T ( D d2bdt2_t0 ) ) ) ( ℝ1 s0 ) =
+      !( D1 dbdt_t0s0 d2bdtds_t0s0 _ ) = dbdt_t0 s0
-  let !( D1 b_t0s0 dbds_t0s0 d2bds2_t0s0 ) = b_t0 ( ℝ1 s0 )
+      !( D1 d2bdt2_t0s0 _ _ ) = d2bdt2_t0 s0
-      !( D1 dbdt_t0s0 d2bdtds_t0s0 _ ) = dbdt_t0 ( ℝ1 s0 )
-      !( D1 d2bdt2_t0s0 _ _ ) = d2bdt2_t0 ( ℝ1 s0 )
   in D2 b_t0s0 ( T dbdt_t0s0 ) dbds_t0s0 ( T d2bdt2_t0s0 ) d2bdtds_t0s0 d2bds2_t0s0
 
 --------------------------------------------------------------------------------
 
-type Diffy :: Type -> Constraint
+type Diffy :: Type -> Type -> Constraint
-class ( Applicative ( D v )
+class ( Traversable ( D v ), Module r ( T v ) ) => Diffy r v where
-      , Traversable ( D v )
+  chain  :: ( Module r ( T w ) ) => D ( ℝ 1 ) v -> D v w -> D ( ℝ 1 ) w
-      , Module Double ( T v ) )
+  konst  :: Module r ( T w ) => w -> D v w
-   => Diffy v where
+  value  :: D v w -> w
-  chain :: ( Module Double ( T w ) )
+  linear :: Module r ( T w ) => ( v -> w ) -> ( v ~> w )
-        => D ( ℝ 1 ) v -> D v w -> D ( ℝ 1 ) w
-  konst :: Module Double ( T w ) => w -> D v w
-  value :: D v w -> w
-  linear :: Module Double ( T w ) => ( v -> w ) -> ( v ~> w )
 
-chainRule :: ( Diffy v, Module Double ( T w ) )
+chainRule :: ( Diffy r v, Module r ( T w ) )
           => ( ( ℝ 1 ) ~> v ) -> ( v ~> w ) -> ( ( ℝ 1 ) ~> w )
 chainRule ( D df ) ( D dg ) =
   D \ x ->
@@ -245,15 +245,16 @@ chainRule ( D df ) ( D dg ) =
         chain df_x ( dg f_x )
 
 -- | Recover the underlying function, discarding all infinitesimal information.
-fun :: forall v w. Diffy v => ( v ~> w ) -> ( v -> w )
+fun :: forall v w r. Diffy r v => ( v ~> w ) -> ( v -> w )
-fun ( D df ) = value @v . df
+fun ( D df ) = value @r @v . df
 {-# INLINE fun #-}
 
-var :: ( Representable u, Diffy u ) => Fin ( Dim u ) -> ( u ~> Double )
+var :: forall u r. ( Module r ( T r ), Representable r u, Diffy r u )
+    => Fin ( Dim u ) -> ( u ~> r )
 var i = linear ( `index` i )
 {-# INLINE var #-}
 
-instance Diffy ( ℝ 0 ) where
+instance Diffy Double ( ℝ 0 ) where
   chain _ ( D0 w ) = D1 w origin origin
   {-# INLINE chain #-}
   konst k = D0 k
@@ -263,7 +264,7 @@ instance Diffy ( ℝ 0 ) where
   linear f = D \ _ -> D0 ( f ℝ0 )
   {-# INLINE linear #-}
 
-instance Diffy ( ℝ 1 ) where
+instance Diffy Double ( ℝ 1 ) where
   chain ( D1 _ ( T ( ℝ1 x' ) ) ( T ( ℝ1 x'' ) ) ) ( D1 v g_x g_xx )
     = D1 v
          ( x' *^ g_x )
@@ -276,7 +277,7 @@ instance Diffy ( ℝ 1 ) where
   linear f = D \ u -> D1 ( f u ) ( T $ f u ) origin
   {-# INLINE linear #-}
 
-instance Diffy ( ℝ 2 ) where
+instance Diffy Double ( ℝ 2 ) where
   chain ( D1 _ ( T ( ℝ2 x' y' ) ) ( T ( ℝ2 x'' y'' ) ) ) ( D2 v g_x g_y g_xx g_xy g_yy )
     = D1 v
          ( x' *^ g_x ^+^ y' *^ g_y )
@@ -294,7 +295,7 @@ instance Diffy ( ℝ 2 ) where
        origin origin origin
   {-# INLINE linear #-}
 
-instance Diffy ( ℝ 3 ) where
+instance Diffy Double ( ℝ 3 ) where
   chain ( D1 _ ( T ( ℝ3 x' y' z' ) ) ( T ( ℝ3 x'' y'' z'' ) ) )
         ( D3 v g_x g_y g_z g_xx g_xy g_yy g_xz g_yz g_zz )
     = D1 v

src/splines/Math/Module.hs

@@ -4,12 +4,12 @@
 
 module Math.Module
   ( Module(..), lerp
-  , Inner(..)
+  , Inner(..), Cross(..)
   , Interpolatable
   , norm, squaredNorm, quadrance, distance
   , proj, projC, closestPointOnSegment
-  , cross
   , strictlyParallel, convexCombination
+  , 𝕀
   )
   where
 
@@ -18,6 +18,10 @@ import Control.Applicative
   ( liftA2 )
 import Control.Monad
   ( guard )
+import Data.Coerce
+  ( coerce )
+import Data.Functor.Identity
+  ( Identity(..) )
 import Data.Kind
   ( Type, Constraint )
 import Data.Monoid
@@ -31,6 +35,14 @@ import Data.Act
     ( (-->) )
   )
 
+-- rounded-hw
+import Numeric.Rounded.Hardware
+  ( Rounded(..) )
+import Numeric.Rounded.Hardware.Interval.NonEmpty
+  ( Interval(..) )
+import Numeric.Rounded.Hardware.Class
+  ( intervalAdd, intervalSub, intervalMul )
+
 -- MetaBrush
 import Math.Epsilon
   ( epsilon )
@@ -69,6 +81,9 @@ infixl 8 ^.^
 class Module r m => Inner r m where
   (^.^) :: m -> m -> r
 
+class Module r m => Cross r m where
+  cross :: m -> m -> r
+
 -- | Norm of a vector, computed using the inner product.
 norm :: forall m r. ( Floating r, Inner r m ) => m -> r
 norm = sqrt . squaredNorm
@@ -114,8 +129,8 @@ closestPointOnSegment c ( Segment p0 p1 )
 
 -- | A convenient constraint synonym for types that support interpolation.
 type Interpolatable :: Type -> Constraint
-class    ( Torsor ( T r ) r, Module Double ( T r ) ) => Interpolatable r
+class    ( Torsor ( T u ) u, Module Double ( T u ) ) => Interpolatable u
-instance ( Torsor ( T r ) r, Module Double ( T r ) ) => Interpolatable r
+instance ( Torsor ( T u ) u, Module Double ( T u ) ) => Interpolatable u
 
 --------------------------------------------------------------------------------
 
@@ -158,9 +173,18 @@ instance Module Double ( T ( ℝ 3 ) ) where
 instance Inner Double ( T ( ℝ 2 ) ) where
   V2 x1 y1 ^.^ V2 x2 y2 = x1 * x2 + y1 * y2
 
--- | Cross-product of two 2D vectors.
+instance Cross Double ( T ( ℝ 2 ) ) where
-cross :: T ( ℝ 2 ) -> T ( ℝ 2 ) -> Double
+  cross ( V2 x1 y1 ) ( V2 x2 y2 ) = x1 * y2 - x2 * y1
-cross ( V2 x1 y1 ) ( V2 x2 y2 ) = x1 * y2 - x2 * y1
+
+instance Module r ( T m ) => Module ( Identity r ) ( T ( Identity m ) ) where
+  origin = coerce $ origin @r @( T m )
+  (^+^)  = coerce $ (^+^)  @r @( T m )
+  (^-^)  = coerce $ (^-^)  @r @( T m )
+  (*^)   = coerce $ (*^)   @r @( T m )
+instance Inner  r ( T m ) => Inner  ( Identity r ) ( T ( Identity m ) ) where
+  (^.^)  = coerce $ (^.^)  @r @( T m )
+instance Cross  r ( T m ) => Cross  ( Identity r ) ( T ( Identity m ) ) where
+  cross  = coerce $ cross  @r @( T m )
 
 -- | Compute whether two vectors point in the same direction,
 -- that is, whether each vector is a (strictly) positive multiple of the other.
@@ -199,3 +223,39 @@ convexCombination v0 v1 u
     c0, c10 :: Double
     c0  = v0 `cross` u
     c10 = ( v0 ^-^ v1 ) `cross` u
+
+--------------------------------------------------------------------------------
+-- Interval arithmetic using rounded-hw library.
+
+instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ 2 ) ) where
+  origin = T ( I ( Rounded ( ℝ2 0 0 ) ) ( Rounded ( ℝ2 0 0 ) ) )
+  T ( I ( Rounded ( ℝ2 x1_lo y1_lo ) ) ( Rounded ( ℝ2 x1_hi y1_hi ) ) ) ^+^
+    T ( I ( Rounded ( ℝ2 x2_lo y2_lo ) ) ( Rounded ( ℝ2 x2_hi y2_hi ) ) )
+      = let !( Rounded x_lo, Rounded x_hi ) = intervalAdd ( Rounded x1_lo ) ( Rounded x1_hi ) ( Rounded x2_lo ) ( Rounded x2_hi )
+            !( Rounded y_lo, Rounded y_hi ) = intervalAdd ( Rounded y1_lo ) ( Rounded y1_hi ) ( Rounded y2_lo ) ( Rounded y2_hi )
+        in T ( I ( Rounded ( ℝ2 x_lo y_lo ) ) ( Rounded ( ℝ2 x_hi y_hi ) ) )
+  T ( I ( Rounded ( ℝ2 x1_lo y1_lo ) ) ( Rounded ( ℝ2 x1_hi y1_hi ) ) ) ^-^
+    T ( I ( Rounded ( ℝ2 x2_lo y2_lo ) ) ( Rounded ( ℝ2 x2_hi y2_hi ) ) )
+      = let !( Rounded x_lo, Rounded x_hi ) = intervalSub ( Rounded x1_lo ) ( Rounded x1_hi ) ( Rounded x2_lo ) ( Rounded x2_hi )
+            !( Rounded y_lo, Rounded y_hi ) = intervalSub ( Rounded y1_lo ) ( Rounded y1_hi ) ( Rounded y2_lo ) ( Rounded y2_hi )
+        in T ( I ( Rounded ( ℝ2 x_lo y_lo ) ) ( Rounded ( ℝ2 x_hi y_hi ) ) )
+  I ( Rounded k_lo ) ( Rounded k_hi ) *^ T ( I ( Rounded ( ℝ2 x1_lo y1_lo ) ) ( Rounded ( ℝ2 x1_hi y1_hi ) ) )
+    = let !( Rounded x_lo, Rounded x_hi ) = intervalMul ( Rounded k_lo ) ( Rounded k_hi ) ( Rounded x1_lo ) ( Rounded x1_hi )
+          !( Rounded y_lo, Rounded y_hi ) = intervalMul ( Rounded k_lo ) ( Rounded k_hi ) ( Rounded y1_lo ) ( Rounded y1_hi )
+      in T ( I ( Rounded ( ℝ2 x_lo y_lo ) ) ( Rounded ( ℝ2 x_hi y_hi ) ) )
+
+instance Inner  ( 𝕀 Double ) ( T ( 𝕀ℝ 2 ) ) where
+  T ( I ( Rounded ( ℝ2 x1_lo y1_lo ) ) ( Rounded ( ℝ2 x1_hi y1_hi ) ) ) ^.^
+    T ( I ( Rounded ( ℝ2 x2_lo y2_lo ) ) ( Rounded ( ℝ2 x2_hi y2_hi ) ) )
+      = let !( x_lo, x_hi ) = intervalMul ( Rounded x1_lo ) ( Rounded x1_hi ) ( Rounded x2_lo ) ( Rounded x2_hi )
+            !( y_lo, y_hi ) = intervalMul ( Rounded y1_lo ) ( Rounded y1_hi ) ( Rounded y2_lo ) ( Rounded y2_hi )
+            !( z_lo, z_hi ) = intervalAdd x_lo x_hi y_lo y_hi
+        in I z_lo z_hi
+
+instance Cross  ( 𝕀 Double ) ( T ( 𝕀ℝ 2 ) ) where
+  T ( I ( Rounded ( ℝ2 x1_lo y1_lo ) ) ( Rounded ( ℝ2 x1_hi y1_hi ) ) ) `cross`
+    T ( I ( Rounded ( ℝ2 x2_lo y2_lo ) ) ( Rounded ( ℝ2 x2_hi y2_hi ) ) )
+      = let !( x_lo, x_hi ) = intervalMul ( Rounded x1_lo ) ( Rounded x1_hi ) ( Rounded y2_lo ) ( Rounded y2_hi )
+            !( y_lo, y_hi ) = intervalMul ( Rounded x2_lo ) ( Rounded x2_hi ) ( Rounded y1_lo ) ( Rounded y1_hi )
+            !( z_lo, z_hi ) = intervalSub x_lo x_hi y_lo y_hi
+        in I z_lo z_hi