more WIP with TH now

MetaBrush.cabal

@@ -75,6 +75,7 @@ common common
         MagicHash
         MultiWayIf
         NamedFieldPuns
+        NoStarIsType
         PatternSynonyms
         RankNTypes
         RecordWildCards
@@ -92,7 +93,6 @@ common common
         -O1
         -fexpose-all-unfoldings
         -funfolding-use-threshold=16
-        -fexcess-precision
         -fspecialise-aggressively
         -optc-O3
         -Wall
@@ -173,20 +173,31 @@ library splines
         Haskell2010
 
     exposed-modules:
-        Math.Bezier.Cubic
+        Math.Algebra.Dual
+      , Math.Bezier.Cubic
       , Math.Bezier.Cubic.Fit
       , Math.Bezier.Quadratic
       , Math.Bezier.Spline
       , Math.Bezier.Stroke
+      , Math.Differentiable
       , Math.Epsilon
+      , Math.Interval
       , Math.Linear
-      , Math.Linear.Dual
       , Math.Linear.Solve
       , Math.Module
+      , Math.Monomial
       , Math.Orientation
+      , Math.Ring
       , Math.Roots
       , Debug.Utils
 
+    other-modules:
+        Math.Algebra.Dual.Internal
+      , Math.Interval.Internal
+      , Math.Linear.Internal
+      , Math.Module.Internal
+      , TH.Utils
+
     build-depends:
         bifunctors
            >= 5.5.4      && < 5.6
@@ -198,6 +209,8 @@ library splines
           ^>= 0.2
       , vector
            >= 0.12.1.2   && < 0.14
+      , template-haskell
+           >= 2.18       && < 2.20
 
 
 library metabrushes

cabal.project

@@ -28,10 +28,10 @@ allow-newer:
   waargonaut:vector, waargonaut:witherable,
 
 -- eigen
-source-repository-package
+--source-repository-package
-  type: git
+--  type: git
-  location: https://github.com/chessai/eigen
+--  location: https://github.com/chessai/eigen
-  tag: 8fff32a43df743c8c83428a86dd566a0936a4fba
+--  tag: 8fff32a43df743c8c83428a86dd566a0936a4fba
 
 -- records-sop
 source-repository-package

src/app/MetaBrush/Application.hs

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

src/app/MetaBrush/Render/Document.hs

@@ -56,6 +56,8 @@ import Control.Monad.Trans.State.Strict
   ( StateT, evalStateT, get, put )
 
 -- MetaBrush
+import Math.Algebra.Dual
+  ( fun )
 import qualified Math.Bezier.Cubic     as Cubic
   ( Bezier(..), fromQuadratic )
 import Math.Bezier.Cubic.Fit
@@ -73,10 +75,10 @@ import Math.Bezier.Stroke
   ( CachedStroke(..), invalidateCache
   , computeStrokeOutline
   )
+import Math.Interval
+  ( Extent(Point) )
 import Math.Linear
-  ( ℝ(..), T(..), Extent(Point) )
+  ( ℝ(..), T(..) )
-import Math.Linear.Dual
-  ( fun )
 import MetaBrush.Asset.Colours
   ( Colours, ColourRecord(..) )
 import MetaBrush.Brush
@@ -305,7 +307,7 @@ strokeRenderData fitParams
                       StrokeWithOutlineRenderData
                         { strokeDataSpline    = spline
                         , strokeOutlineData   = ( outline, fitPts )
-                        , strokeBrushFunction = fun ( brushFn @Point proxy# id )
+                        , strokeBrushFunction = fun @Double ( brushFn @Point proxy# id )
                                               . embedUsedParams
                                               . toUsedParams
                         }

src/metabrushes/MetaBrush/Asset/Brushes.hs

@@ -1,10 +1,14 @@
 {-# LANGUAGE AllowAmbiguousTypes #-}
 {-# LANGUAGE OverloadedStrings   #-}
 {-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TemplateHaskell     #-}
 
 module MetaBrush.Asset.Brushes where
 
 -- base
+import Prelude
+  hiding
+    ( Num(..), Floating(..) )
 import GHC.Exts
   ( Proxy# )
 
@@ -23,12 +27,16 @@ import qualified Data.HashMap.Strict as HashMap
   ( fromList, lookup )
 
 -- MetaBrush
+import Math.Algebra.Dual
 import Math.Bezier.Spline
+import Math.Differentiable
+  ( DiffInterp )
+import Math.Interval
+  ( type I )
 import Math.Linear
-import Math.Linear.Dual
-  ( D, type (~>)(..), Differentiable, Diffy(konst), var )
 import Math.Module
   ( Module((^+^), (*^)) )
+import Math.Ring
 import MetaBrush.Brush
   ( Brush(..), SomeBrush(..), WithParams(..) )
 import MetaBrush.Records
@@ -70,62 +78,62 @@ ellipse = BrushData "ellipse" ( WithParams deflts ellipseBrush )
 --------------------------------------------------------------------------------
 -- Differentiable brushes.
 
-circleSpline :: forall i u v
+circleSpline :: forall i k u v
-             .  Applicative ( D ( I i u ) )
+             .  Applicative ( D k ( I i u ) )
-             => ( Double -> Double -> D ( I i u ) ( I i v ) )
+             => ( Double -> Double -> D k ( I i u ) ( I i v ) )
-             -> D ( I i u ) ( Spline 'Closed () ( I i v ) )
+             -> D k ( I i u ) ( Spline 'Closed () ( I i v ) )
 circleSpline p = sequenceA $
   Spline { splineStart  = p 1 0
          , splineCurves = ClosedCurves crvs lastCrv }
   where
     crvs = Seq.fromList
-      [ Bezier3To ( p  1  κ ) ( p  κ  1 ) ( NextPoint (p  0  1) ) ()
+      [ Bezier3To ( p  1  κ ) ( p  κ  1 ) ( NextPoint ( p  0  1 ) ) ()
-      , Bezier3To ( p -κ  1 ) ( p -1  κ ) ( NextPoint (p -1  0) ) ()
+      , Bezier3To ( p -κ  1 ) ( p -1  κ ) ( NextPoint ( p -1  0 ) ) ()
-      , Bezier3To ( p -1 -κ ) ( p -κ -1 ) ( NextPoint (p  0 -1) ) ()
+      , Bezier3To ( p -1 -κ ) ( p -κ -1 ) ( NextPoint ( p  0 -1 ) ) ()
       ]
     lastCrv =
         Bezier3To ( p  κ -1 ) ( p  1 -κ ) BackToStart ()
 
-circleBrush :: forall i
+circleBrush :: forall i k
-            .  ( Differentiable i ( Record CircleBrushFields ) )
+            .  ( k ~ 2, DiffInterp i ( Record CircleBrushFields ) )
             => Proxy# i
             -> ( forall a. a -> I i a )
-            -> I i ( Record CircleBrushFields ) ~> Spline 'Closed () ( I i ( ℝ 2 ) )
+            -> C k ( I i ( Record CircleBrushFields ) ) ( Spline 'Closed () ( I i ( ℝ 2 ) ) )
 circleBrush _ mkI =
   D \ params ->
-    let r :: D ( I i ( Record CircleBrushFields ) ) ( I i Double )
+    let r :: D k ( I i ( Record CircleBrushFields ) ) ( I i Double )
-        r = runD ( var ( Fin 1## ) ) params
+        r = runD ( var @_ @k ( Fin 1 ) ) params
-        mkPt :: Double -> Double -> D ( I i ( Record CircleBrushFields ) ) ( I i ( ℝ 2 ) )
+        mkPt :: Double -> Double -> D k ( I i ( Record CircleBrushFields ) ) ( I i ( ℝ 2 ) )
         mkPt ( kon -> x ) ( kon -> y )
           =   ( x * r ) *^ e_x
           ^+^ ( y * r ) *^ e_y
-    in circleSpline @i @( Record CircleBrushFields ) @( ℝ 2 ) mkPt
+    in circleSpline @i @k @( Record CircleBrushFields ) @( ℝ 2 ) mkPt
   where
-    e_x, e_y :: D ( I i ( Record CircleBrushFields ) ) ( I i ( ℝ 2 ) )
+    e_x, e_y :: D k ( I i ( Record CircleBrushFields ) ) ( I i ( ℝ 2 ) )
     e_x = pure $ mkI $ ℝ2 1 0
     e_y = pure $ mkI $ ℝ2 0 1
 
-    kon = konst @( I i Double ) @( I i ( Record CircleBrushFields ) ) . mkI
+    kon = konst @( I i Double ) @k @( I i ( Record CircleBrushFields ) ) . mkI
 
-ellipseBrush :: forall i
+ellipseBrush :: forall i k
-             .  ( Differentiable i ( Record EllipseBrushFields ) )
+             .  ( k ~ 2, DiffInterp i ( Record EllipseBrushFields ) )
              => Proxy# i
              -> ( forall a. a -> I i a )
-             -> I i ( Record EllipseBrushFields ) ~> Spline 'Closed () ( I i ( ℝ 2 ) )
+             -> C k ( I i ( Record EllipseBrushFields ) ) ( Spline 'Closed () ( I i ( ℝ 2 ) ) )
 ellipseBrush _ mkI =
   D \ params ->
-    let a, b, phi :: D ( I i ( Record EllipseBrushFields ) ) ( I i Double )
+    let a, b, phi :: D k ( I i ( Record EllipseBrushFields ) ) ( I i Double )
-        a   = runD ( var ( Fin 1## ) ) params
+        a   = runD ( var @_ @k ( Fin 1 ) ) params
-        b   = runD ( var ( Fin 2## ) ) params
+        b   = runD ( var @_ @k ( Fin 2 ) ) params
-        phi = runD ( var ( Fin 3## ) ) params
+        phi = runD ( var @_ @k ( Fin 3 ) ) params
-        mkPt :: Double -> Double -> D ( I i ( Record EllipseBrushFields ) ) ( I i ( ℝ 2 ) )
+        mkPt :: Double -> Double -> D k ( I i ( Record EllipseBrushFields ) ) ( I i ( ℝ 2 ) )
         mkPt ( kon -> x ) ( kon -> y )
           =   ( x * a * cos phi - y * b * sin phi ) *^ e_x
           ^+^ ( y * b * cos phi + x * a * sin phi ) *^ e_y
-    in circleSpline @i @( Record EllipseBrushFields ) @( ℝ 2 ) mkPt
+    in circleSpline @i @k @( Record EllipseBrushFields ) @( ℝ 2 ) mkPt
   where
-    e_x, e_y :: D ( I i ( Record EllipseBrushFields ) ) ( I i ( ℝ 2 ) )
+    e_x, e_y :: D k ( I i ( Record EllipseBrushFields ) ) ( I i ( ℝ 2 ) )
     e_x = pure $ mkI $ ℝ2 1 0
     e_y = pure $ mkI $ ℝ2 0 1
 
-    kon = konst @( I i Double ) @( I i ( Record EllipseBrushFields ) ) . mkI
+    kon = konst @( I i Double ) @k @( I i ( Record EllipseBrushFields ) ) . mkI

src/metabrushes/MetaBrush/Brush.hs

@@ -39,12 +39,16 @@ import qualified Data.Text as Text
   ( unpack )
 
 -- MetaBrush
+import Math.Algebra.Dual
+  ( type (~>) )
 import Math.Linear
-  ( ℝ, type I, Extent(Point, Interval) )
+import Math.Interval
-import Math.Linear.Dual
+  ( type I, Extent(Point, Interval) )
-  ( type (~>), Differentiable )
+
 import Math.Bezier.Spline
   ( SplineType(Closed), Spline )
+import Math.Differentiable
+  ( DiffInterp )
 import MetaBrush.Records
   ( KnownSymbols, Length, Record )
 import MetaBrush.Serialisable
@@ -54,30 +58,33 @@ import MetaBrush.Serialisable
 
 -- | A differentiable function from a given record type,
 -- with provided default values that can be overridden.
-type WithParams :: [ Symbol ] -> ( Type -> Type ) -> Type
+type WithParams :: Type -> ( Type -> Type ) -> Type
 data WithParams params f =
   WithParams
-    { defaultParams :: Record params
+    { defaultParams :: params
-    , withParams    :: forall i
+    , withParams
-                    .  ( Differentiable i ( Record params ) )
+        :: forall i
-                    => Proxy# i
+        .  ( DiffInterp i params )
-                    -> ( forall a. a -> I i a )
+        => Proxy# i
-                    -> I i ( Record params ) ~> f ( I i ( ℝ 2 ) )
+        -> ( forall a. a -> I i a )
+        -> I i params ~> f ( I i ( ℝ 2 ) )
     }
 
 --------------------------------------------------------------------------------
 
 -- | A brush function: a function from a record of parameters to a closed spline.
 type BrushFunction :: [ Symbol ] -> Type
-type BrushFunction brushFields = WithParams brushFields ( Spline Closed () )
+type BrushFunction brushFields =
+  WithParams ( Record brushFields ) ( Spline Closed () )
 
 type Brush :: [ Symbol ] -> Type
 data Brush brushFields where
   BrushData
     :: forall brushFields
     .  ( KnownSymbols brushFields, Typeable brushFields
-       , Differentiable Point    ( ℝ ( Length brushFields ) )
+       , Representable Double ( ℝ ( Length brushFields ) )
-       , Differentiable Interval ( ℝ ( Length brushFields ) )
+       , DiffInterp Point     ( ℝ ( Length brushFields ) )
+       , DiffInterp Interval  ( ℝ ( Length brushFields ) )
        )
     => { brushName     :: !Text
        , brushFunction :: BrushFunction brushFields
@@ -111,21 +118,25 @@ class ( KnownSymbols pointFields, Typeable pointFields
       , Serialisable   ( Record pointFields )
       , Show           ( Record pointFields )
       , NFData         ( Record pointFields )
-      , Differentiable Point    ( ℝ ( Length pointFields ) )
+      , Representable  Double   ( ℝ ( Length pointFields ) )
-      , Differentiable Interval ( ℝ ( Length pointFields ) )
+      , DiffInterp Point    ( ℝ ( Length pointFields ) )
+      , DiffInterp Interval ( ℝ ( Length pointFields ) )
       )
    => PointFields pointFields where { }
 instance ( KnownSymbols pointFields, Typeable pointFields
          , Serialisable   ( Record pointFields )
          , Show           ( Record pointFields )
          , NFData         ( Record pointFields )
-         , Differentiable Point    ( ℝ ( Length pointFields ) )
+         , Representable  Double   ( ℝ ( Length pointFields ) )
-         , Differentiable Interval ( ℝ ( Length pointFields ) )
+         , DiffInterp Point    ( ℝ ( Length pointFields ) )
+         , DiffInterp Interval ( ℝ ( Length pointFields ) )
          )
       => PointFields pointFields where { }
 
 -- | Assumes the input has no duplicates (doesn't check.)
-provePointFields :: [ Text ] -> ( forall pointFields. PointFields pointFields => Proxy# pointFields -> r ) -> r
+provePointFields :: [ Text ]
+                 -> ( forall pointFields. PointFields pointFields => Proxy# pointFields -> r )
+                 -> r
 provePointFields fieldNames k =
   case fieldNames of
     []

src/metabrushes/MetaBrush/Records.hs

@@ -17,7 +17,7 @@ import Data.Typeable
 import Data.Type.Equality
   ( (:~:)(Refl) )
 import GHC.Exts
-  ( Word(W#), Proxy#, proxy# )
+  ( Proxy#, proxy# )
 import GHC.Show
   ( showCommaSpace )
 import GHC.TypeLits
@@ -54,9 +54,10 @@ import qualified Data.Text as Text
   ( pack, unpack )
 
 -- MetaBrush
+import Math.Algebra.Dual
+import Math.Differentiable
+import Math.Interval
 import Math.Linear
-import Math.Linear.Dual
-  ( type (~>)(..), D, Diffy, Differentiable )
 import Math.Module
   ( Module )
 
@@ -87,7 +88,7 @@ instance ( KnownSymbols ks, Representable Double ( ℝ ( Length ks ) ) )
       where
         fields :: [ ShowS ]
         fields =
-          zip [ 1.. ] ( knownSymbols @ks ) <&> \ ( W# i, fld ) ->
+          zip [ 1.. ] ( knownSymbols @ks ) <&> \ ( i, fld ) ->
             let v = index r ( Fin i )
             in showString ( Text.unpack fld ) . showString " = " . showsPrec 0 v
 
@@ -143,17 +144,18 @@ instance ( Torsor ( T ( 𝕀ℝ ( Length ks ) ) ) ( 𝕀ℝ ( Length ks ) )
           T ( I ( Rounded c_lo ) ( Rounded c_hi ) ) ->
             T ( I ( Rounded ( MkR c_lo ) ) ( Rounded ( MkR c_hi ) ) )
 
+type instance RepDim ( Record ks ) = Length ks
 deriving newtype
   instance Representable r ( ℝ ( Length ks ) )
         => Representable r ( Record ks )
 
-type instance D ( Record ks ) = D ( ℝ ( Length ks ) )
+type instance D k ( Record ks ) = D k ( ℝ ( Length ks ) )
-deriving newtype instance Diffy Double ( ℝ ( Length ks ) )
+deriving newtype instance HasChainRule Double 2 ( ℝ ( Length ks ) )
-                       => Diffy Double ( Record ks )
+                       => HasChainRule Double 2 ( Record ks )
 
 deriving via 𝕀ℝ ( Length ks )
-  instance Diffy ( 𝕀 Double ) ( 𝕀ℝ ( Length ks ) )
+  instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ ( Length ks ) )
-        => Diffy ( 𝕀 Double ) ( 𝕀 ( Record ks ) )
+        => HasChainRule ( 𝕀 Double ) 2 ( 𝕀 ( Record ks ) )
 
 --------------------------------------------------------------------------------
 
@@ -181,6 +183,7 @@ intersect :: forall r1 r2 l1 l2
             , KnownSymbols r1, KnownSymbols r2
             , l1 ~ Length r1, l2 ~ Length r2
             , Representable Double ( ℝ l1 )
+            , Representable Double ( ℝ l2 )
             , Differentiable 'Point    ( ℝ l2 )
             , Differentiable 'Interval ( ℝ l2 )
             )
@@ -197,7 +200,7 @@ intersect
         project = \ ( MkR r1 ) -> MkR $ projection ( (!) r1_idxs ) r1
 
         inject :: Record r2 -> Record r1r2 -> Record r2
-        inject = \ ( MkR r2 ) -> \ ( MkR r1r2 ) -> MkR $ injection ( find eqFin r2_idxs ) r1r2 r2
+        inject = \ ( MkR r2 ) -> \ ( MkR r1r2 ) -> MkR $ injection ( \ i -> find ( == i ) r2_idxs ) r1r2 r2
       in Intersection { project, inject }
 
 data Intersection r1 r2 where
@@ -205,6 +208,7 @@ data Intersection r1 r2 where
     :: forall r1r2 r1 r2 l12
     .  ( l12 ~ Length r1r2
        , KnownSymbols r1r2
+       , Representable Double ( ℝ l12 )
        , Differentiable 'Point    ( ℝ l12 )
        , Differentiable 'Interval ( ℝ l12 )
        )
@@ -221,13 +225,13 @@ doIntersection
   :: forall r1 r2 l1 l2 kont
   .  ( KnownSymbols r1, KnownSymbols r2
      , l1 ~ Length r1, l2 ~ Length r2
-     , Representable Double ( ℝ l1 )
-     , Representable Double ( ℝ l2 )
      )
   => ( forall r1r2 l12.
-        ( r1r2 ~ Intersect r1 r2, KnownSymbols r1r2, l12 ~ Length r1r2
+        ( r1r2 ~ Intersect r1 r2
+        , KnownSymbols r1r2, l12 ~ Length r1r2
         , Differentiable 'Point    ( ℝ l12 )
         , Differentiable 'Interval ( ℝ l12 )
+        , Representable Double ( ℝ l12 )
         )
       => Proxy# r1r2 -> Vec l12 ( Fin l1 ) -> Vec l12 ( Fin l2 ) -> kont )
   -> kont
@@ -237,11 +241,11 @@ doIntersection k =
     [ ]
       | ( _ :: Proxy# r1r2 ) <- proxy# @'[ ]
       , Refl <- ( unsafeCoerce Refl :: r1r2 :~: Intersect r1 r2 )
-      -> k @'[] proxy#
+      -> k @r1r2 proxy#
            VZ
            VZ
 
-    [   ( f1, W# r1_i1, W# r2_i1 ) ]
+    [   ( f1, r1_i1, r2_i1 ) ]
       | SomeSymbol @f1 _     <- someSymbolVal ( Text.unpack f1 )
       , ( _ :: Proxy# r1r2 ) <- proxy# @'[ f1 ]
       , Refl <- ( unsafeCoerce Refl :: r1r2 :~: Intersect r1 r2 )
@@ -249,8 +253,8 @@ doIntersection k =
            ( VS ( Fin r1_i1 ) VZ )
            ( VS ( Fin r2_i1 ) VZ )
 
-    [   ( f1, W# r1_i1, W# r2_i1 )
+    [   ( f1, r1_i1, r2_i1 )
-      , ( f2, W# r1_i2, W# r2_i2 ) ]
+      , ( f2, r1_i2, r2_i2 ) ]
       | SomeSymbol @f1 _     <- someSymbolVal ( Text.unpack f1 )
       , SomeSymbol @f2 _     <- someSymbolVal ( Text.unpack f2 )
       , ( _ :: Proxy# r1r2 ) <- proxy# @'[ f1, f2 ]
@@ -259,9 +263,9 @@ doIntersection k =
            ( VS ( Fin r1_i1 ) $ VS ( Fin r1_i2 ) VZ )
            ( VS ( Fin r2_i1 ) $ VS ( Fin r2_i2 ) VZ )
 
-    [   ( f1, W# r1_i1, W# r2_i1 )
+    [   ( f1, r1_i1, r2_i1 )
-      , ( f2, W# r1_i2, W# r2_i2 )
+      , ( f2, r1_i2, r2_i2 )
-      , ( f3, W# r1_i3, W# r2_i3 ) ]
+      , ( f3, r1_i3, r2_i3 ) ]
       | SomeSymbol @f1 _     <- someSymbolVal ( Text.unpack f1 )
       , SomeSymbol @f2 _     <- someSymbolVal ( Text.unpack f2 )
       , SomeSymbol @f3 _     <- someSymbolVal ( Text.unpack f3 )
@@ -271,6 +275,20 @@ doIntersection k =
            ( VS ( Fin r1_i1 ) $ VS ( Fin r1_i2 ) $ VS ( Fin r1_i3 ) VZ )
            ( VS ( Fin r2_i1 ) $ VS ( Fin r2_i2 ) $ VS ( Fin r2_i3 ) VZ )
 
+    [   ( f1, r1_i1, r2_i1 )
+      , ( f2, r1_i2, r2_i2 )
+      , ( f3, r1_i3, r2_i3 )
+      , ( f4, r1_i4, r2_i4 ) ]
+      | SomeSymbol @f1 _     <- someSymbolVal ( Text.unpack f1 )
+      , SomeSymbol @f2 _     <- someSymbolVal ( Text.unpack f2 )
+      , SomeSymbol @f3 _     <- someSymbolVal ( Text.unpack f3 )
+      , SomeSymbol @f4 _     <- someSymbolVal ( Text.unpack f4 )
+      , ( _ :: Proxy# r1r2 ) <- proxy# @'[ f1, f2, f3, f4 ]
+      , Refl <- ( unsafeCoerce Refl :: r1r2 :~: Intersect r1 r2 )
+      -> k @r1r2 proxy#
+           ( VS ( Fin r1_i1 ) $ VS ( Fin r1_i2 ) $ VS ( Fin r1_i3 ) $ VS ( Fin r1_i4 ) VZ )
+           ( VS ( Fin r2_i1 ) $ VS ( Fin r2_i2 ) $ VS ( Fin r2_i3 ) $ VS ( Fin r2_i4 ) VZ )
+
     other -> error $ "Intersection not defined in dimension " ++ show ( length other )
 
 ------

src/metabrushes/MetaBrush/Serialisable.hs

@@ -29,8 +29,6 @@ import Data.STRef
   ( newSTRef )
 import Data.Traversable
   ( for )
-import GHC.Exts
-  ( Word(W#) )
 
 -- containers
 import Data.Map.Strict
@@ -116,14 +114,14 @@ instance ( KnownSymbols ks, Representable Double ( ℝ ( Length ks ) ) )
     where
       encodeFields :: Record ks -> [ ( Text, Double ) ]
       encodeFields ( MkR r ) =
-        zip [1..] ( knownSymbols @ks ) <&> \ ( W# i, fld ) ->
+        zip [1..] ( knownSymbols @ks ) <&> \ ( i, fld ) ->
           ( fld, index r ( Fin i ) )
 
   decoder = fmap decodeFields $ for ( knownSymbols @ks ) \ k -> JSON.Decoder.atKey k ( decoder @Double )
     where
       decodeFields :: [ Double ] -> Record ks
-      decodeFields coords = MkR $ tabulate \ ( Fin i# ) ->
+      decodeFields coords = MkR $ tabulate \ ( Fin i ) ->
-        coords !! ( fromIntegral ( W# i# ) - 1 )
+        coords !! ( fromIntegral i - 1 )
 
 --------------------------------------------------------------------------------
 

src/splines/Math/Algebra/Dual.hs

@@ -0,0 +1,761 @@
+{-# LANGUAGE AllowAmbiguousTypes  #-}
+{-# LANGUAGE RebindableSyntax     #-}
+{-# LANGUAGE ScopedTypeVariables  #-}
+{-# LANGUAGE TemplateHaskell      #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+{-# OPTIONS_GHC -Wno-orphans -O2 #-}
+
+{-# OPTIONS_GHC -ddump-splices -ddump-simpl -ddump-to-file -dno-typeable-binds
+  -dsuppress-unfoldings -dsuppress-coercions #-}
+
+module Math.Algebra.Dual
+  ( C(..), D, type (~>), type (~~>)
+  , HasChainRule(..), chainRule
+  , uncurryD2, uncurryD3
+  , linear, fun, var
+
+  , D𝔸0(..)
+  , D1𝔸1(..), D2𝔸1(..), D3𝔸1(..)
+  , D1𝔸2(..), D2𝔸2(..), D3𝔸2(..)
+  , D1𝔸3(..), D2𝔸3(..), D3𝔸3(..)
+  , D1𝔸4(..), D2𝔸4(..), D3𝔸4(..)
+  ) where
+
+-- base
+import Prelude hiding ( Num(..), Floating(..), (^) )
+import Control.Applicative
+  ( liftA2 )
+import Data.Coerce
+  ( coerce )
+import Data.Kind
+  ( Type )
+import Data.Monoid
+  ( Ap(..) )
+import GHC.TypeNats
+  ( Nat )
+
+-- MetaBrush
+import Math.Algebra.Dual.Internal
+import Math.Linear
+import Math.Module
+import Math.Monomial
+import Math.Ring
+
+--------------------------------------------------------------------------------
+
+-- | @C n u v@ is the space of @C^k@-differentiable maps from @u@ to @v@.
+type C :: Nat -> Type -> Type -> Type
+newtype C k u v = D { runD :: u -> D k u v }
+deriving stock instance Functor ( D k u ) => Functor ( C k u )
+
+-- | \( C^2 \)-differentiable mappings.
+type (~>)  = C 2
+-- | \( C^3 \)-differentiable mappings.
+type (~~>) = C 3
+
+-- | @D k u v@ is the space of @k@-th order germs of functions from @u@ to @v@,
+-- represented by the algebra:
+--
+-- \[ \mathbb{Z}[x_1, \ldots, x_n]/(x_1, \ldots, x_n)^(k+1) \otimes_\mathbb{Z} v \]
+--
+-- when @u@ is of dimension @n@.
+type D :: Nat -> Type -> Type -> Type
+type family D k u
+
+type instance D k ( ℝ 0 ) = D𝔸0
+
+type instance D 1 ( ℝ 1 ) = D1𝔸1
+type instance D 1 ( ℝ 2 ) = D1𝔸2
+type instance D 1 ( ℝ 3 ) = D1𝔸3
+type instance D 1 ( ℝ 4 ) = D1𝔸4
+
+type instance D 2 ( ℝ 1 ) = D2𝔸1
+type instance D 2 ( ℝ 2 ) = D2𝔸2
+type instance D 2 ( ℝ 3 ) = D2𝔸3
+type instance D 2 ( ℝ 4 ) = D2𝔸4
+
+type instance D 3 ( ℝ 1 ) = D3𝔸1
+type instance D 3 ( ℝ 2 ) = D3𝔸2
+type instance D 3 ( ℝ 3 ) = D3𝔸3
+type instance D 3 ( ℝ 4 ) = D3𝔸4
+
+--------------------------------------------------------------------------------
+
+-- Weird instance needed in just one place;
+-- see use of chain in 'Math.Bezier.Stroke.brushStrokeData'.
+instance ( Applicative ( D k u ), Module r ( T v ) )
+       => Module r ( T ( C k u v ) ) where
+  origin = T $ D \ _ -> pure $ coerce $ origin @r @( T v )
+  T ( D f ) ^+^ T ( D g ) = T $ D \ t -> liftA2 ( coerce $ (^+^) @r @( T v ) ) ( 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 -> fmap ( coerce $ (*^) @r @( T v ) a ) $ f t
+
+--------------------------------------------------------------------------------
+
+-- | @HasChainRule r k v@ means we have a chain rule
+-- with @D k v w@ in the middle, for any @r@-module @w@.
+class HasChainRule r k v where
+  chain   :: Module r ( T w )
+          => D k ( ℝ 1 ) v -> D k v w -> D k ( ℝ 1 ) w
+  konst   :: AbelianGroup w => w -> D k v w
+  value   :: D k v w -> w
+  linearD :: Module r ( T w ) => ( v -> w ) -> v -> D k v w
+
+linear :: forall k r v w
+       .  ( HasChainRule r k v, Module r ( T w ) )
+       => ( v -> w ) -> C k v w
+linear f = D \ x -> linearD @r @k @v @w f x
+
+chainRule :: forall r k u v w
+          .  ( HasChainRule r k v, Module r ( T w )
+             , D k u ~ D k ( ℝ 1 ), HasChainRule r k u
+             )
+          => C k u v -> C k v w -> C k u w
+chainRule ( D df ) ( D dg ) =
+  D \ x ->
+    case df x of
+      df_x ->
+        chain @r @k @v df_x ( dg $ value @r @k @u df_x )
+
+
+uncurryD2 :: D 2 a ~ D 2 ( ℝ 1 )
+          => D 2 ( ℝ 1 ) ( C 2 a b ) -> a -> D 2 ( ℝ 2 ) b
+uncurryD2 ( D21 ( D b_t0 ) ( T ( D dbdt_t0 ) ) ( T ( D d2bdt2_t0 ) ) ) s0 =
+  let !( D21 b_t0s0 dbds_t0s0 d2bds2_t0s0 ) = b_t0 s0
+      !( D21 dbdt_t0s0 d2bdtds_t0s0 _ ) = dbdt_t0 s0
+      !( D21 d2bdt2_t0s0 _ _ ) = d2bdt2_t0 s0
+  in D22
+       b_t0s0
+       ( T dbdt_t0s0 ) dbds_t0s0
+       ( T d2bdt2_t0s0 ) d2bdtds_t0s0 d2bds2_t0s0
+
+uncurryD3 :: D 3 a ~ D 3 ( ℝ 1 )
+          => D 3 ( ℝ 1 ) ( C 3 a b ) -> a -> D 3 ( ℝ 2 ) b
+uncurryD3 ( D31 ( D b_t0 ) ( T ( D dbdt_t0 ) ) ( T ( D d2bdt2_t0 ) ) ( T ( D d3bdt3_t0 ) ) ) s0 =
+  let !( D31 b_t0s0 dbds_t0s0 d2bds2_t0s0 d3bds3_t0s0 ) = b_t0 s0
+      !( D31 dbdt_t0s0 d2bdtds_t0s0 d3bdtds2_t0s0 _ ) = dbdt_t0 s0
+      !( D31 d2bdt2_t0s0 d3bdt2ds_t0s0 _ _ ) = d2bdt2_t0 s0
+      !( D31 d3bdt3_t0s0 _ _ _ ) = d3bdt3_t0 s0
+  in D32
+       b_t0s0
+       ( T dbdt_t0s0 ) dbds_t0s0
+       ( T d2bdt2_t0s0 ) d2bdtds_t0s0 d2bds2_t0s0
+       ( T d3bdt3_t0s0 ) d3bdt2ds_t0s0 d3bdtds2_t0s0 d3bds3_t0s0
+
+-- | Recover the underlying function, discarding all infinitesimal information.
+fun :: forall r k v w. HasChainRule r k v => C k v w -> ( v -> w )
+fun ( D df ) = value @r @k @v . df
+{-# INLINE fun #-}
+
+-- | The differentiable germ of a coordinate variable.
+var :: forall r k v
+    .  ( Module r ( T r ), Representable r v, HasChainRule r k v )
+    => Fin ( RepDim v ) -> C k v r
+var i = D $ linearD @r @k @v ( `index` i )
+{-# INLINE var #-}
+
+--------------------------------------------------------------------------------
+
+-- | Newtype for the module instance @Module r v => Module ( dr r ) ( dr v )@.
+newtype ApAp r dr v = ApAp { unApAp :: dr v }
+
+instance ( Ring ( dr r ), Module r ( T v ), Applicative dr )
+       => Module ( dr r ) ( ApAp r dr v ) where
+  ApAp !u ^+^ ApAp !v = ApAp $ liftA2 ( coerce $ (^+^) @r @( T v ) ) u v
+  ApAp !u ^-^ ApAp !v = ApAp $ liftA2 ( coerce $ (^-^) @r @( T v ) ) u v
+  origin = ApAp $ pure $ coerce $ origin @r @( T v )
+  !k *^ ApAp !u = ApAp $ liftA2 ( coerce $ (*^) @r @( T v ) ) k u
+
+
+deriving via ApAp r D𝔸0 v
+  instance Module r ( T v ) => Module ( D𝔸0  r ) ( D𝔸0  v )
+--deriving via ApAp r D1𝔸1 v
+--  instance Module r ( T v ) => Module ( D1𝔸1 r ) ( D1𝔸1 v )
+deriving via ApAp r D2𝔸1 v
+  instance Module r ( T v ) => Module ( D2𝔸1 r ) ( D2𝔸1 v )
+deriving via ApAp r D3𝔸1 v
+  instance Module r ( T v ) => Module ( D3𝔸1 r ) ( D3𝔸1 v )
+--deriving via ApAp r D1𝔸2 v
+--  instance Module r ( T v ) => Module ( D1𝔸2 r ) ( D1𝔸2 v )
+deriving via ApAp r D2𝔸2 v
+  instance Module r ( T v ) => Module ( D2𝔸2 r ) ( D2𝔸2 v )
+deriving via ApAp r D3𝔸2 v
+  instance Module r ( T v ) => Module ( D3𝔸2 r ) ( D3𝔸2 v )
+--deriving via ApAp r D1𝔸3 v
+--  instance Module r ( T v ) => Module ( D1𝔸3 r ) ( D1𝔸3 v )
+deriving via ApAp r D2𝔸3 v
+  instance Module r ( T v ) => Module ( D2𝔸3 r ) ( D2𝔸3 v )
+deriving via ApAp r D3𝔸3 v
+  instance Module r ( T v ) => Module ( D3𝔸3 r ) ( D3𝔸3 v )
+--deriving via ApAp r D1𝔸4 v
+--  instance Module r ( T v ) => Module ( D1𝔸4 r ) ( D1𝔸4 v )
+deriving via ApAp r D2𝔸4 v
+  instance Module r ( T v ) => Module ( D2𝔸4 r ) ( D2𝔸4 v )
+deriving via ApAp r D3𝔸4 v
+  instance Module r ( T v ) => Module ( D3𝔸4 r ) ( D3𝔸4 v )
+
+--------------------------------------------------------------------------------
+-- AbelianGroup instances
+
+newtype ApAp2 k u r = ApAp2 { unApAp2 :: D k u r }
+
+instance ( Applicative ( D k u )
+         , AbelianGroup r
+         , HasChainRule Double k u
+         ) => AbelianGroup ( ApAp2 k u r ) where
+  ApAp2 !x + ApAp2 !y = ApAp2 $ liftA2 ( (+) @r ) x y
+  ApAp2 !x - ApAp2 !y = ApAp2 $ liftA2 ( (-) @r ) x y
+  negate ( ApAp2 !x ) = ApAp2 $ fmap ( negate @r ) x
+
+  -- DO NOT USE PURE!!
+  fromInteger !i = ApAp2 $ konst @Double @k @u ( fromInteger @r i )
+
+deriving newtype instance AbelianGroup r => AbelianGroup ( D𝔸0 r )
+
+deriving via ApAp2 2 ( ℝ 1 ) r
+  instance AbelianGroup r => AbelianGroup ( D2𝔸1 r )
+deriving via ApAp2 3 ( ℝ 1 ) r
+  instance AbelianGroup r => AbelianGroup ( D3𝔸1 r )
+deriving via ApAp2 2 ( ℝ 2 ) r
+  instance AbelianGroup r => AbelianGroup ( D2𝔸2 r )
+deriving via ApAp2 3 ( ℝ 2 ) r
+  instance AbelianGroup r => AbelianGroup ( D3𝔸2 r )
+deriving via ApAp2 2 ( ℝ 3 ) r
+  instance AbelianGroup r => AbelianGroup ( D2𝔸3 r )
+deriving via ApAp2 3 ( ℝ 3 ) r
+  instance AbelianGroup r => AbelianGroup ( D3𝔸3 r )
+deriving via ApAp2 2 ( ℝ 4 ) r
+  instance AbelianGroup r => AbelianGroup ( D2𝔸4 r )
+deriving via ApAp2 3 ( ℝ 4 ) r
+  instance AbelianGroup r => AbelianGroup ( D3𝔸4 r )
+
+--------------------------------------------------------------------------------
+-- Ring instances.
+
+-- TODO: should also implement a power operation specifically,
+-- as this is important for interval arithmethic.
+
+deriving newtype instance Ring r => Ring ( D𝔸0 r )
+
+--instance Ring r => Ring ( D1𝔸1 r ) where
+--  !dr1 * !dr2 =
+--    let
+--      o :: r
+--      o = fromInteger 0
+--      p :: r -> r -> r
+--      p = (+) @r
+--      m :: r -> r -> r
+--      m = (*) @r
+--   in
+--       $$( prodRuleQ
+--             [|| o ||] [|| p ||] [|| m ||]
+--             [|| dr1 ||] [|| dr2 ||] )
+
+instance Ring r => Ring ( D2𝔸1 r ) where
+  !dr1 * !dr2 =
+    let
+      o :: r
+      o = fromInteger 0
+      p :: r -> r -> r
+      p = (+) @r
+      m :: r -> r -> r
+      m = (*) @r
+   in
+       $$( prodRuleQ
+             [|| o ||] [|| p ||] [|| m ||]
+             [|| dr1 ||] [|| dr2 ||] )
+
+instance Ring r => Ring ( D3𝔸1 r ) where
+  !dr1 * !dr2 =
+    let
+      o :: r
+      o = fromInteger 0
+      p :: r -> r -> r
+      p = (+) @r
+      m :: r -> r -> r
+      m = (*) @r
+   in
+       $$( prodRuleQ
+             [|| o ||] [|| p ||] [|| m ||]
+             [|| dr1 ||] [|| dr2 ||] )
+
+--instance Ring r => Ring ( D1𝔸2 r ) where
+--  !dr1 * !dr2 =
+--    let
+--      o :: r
+--      o = fromInteger 0
+--      p :: r -> r -> r
+--      p = (+) @r
+--      m :: r -> r -> r
+--      m = (*) @r
+--   in
+--       $$( prodRuleQ
+--             [|| o ||] [|| p ||] [|| m ||]
+--             [|| dr1 ||] [|| dr2 ||] )
+
+instance Ring r => Ring ( D2𝔸2 r ) where
+  !dr1 * !dr2 =
+    let
+      o :: r
+      o = fromInteger 0
+      p :: r -> r -> r
+      p = (+) @r
+      m :: r -> r -> r
+      m = (*) @r
+   in
+       $$( prodRuleQ
+             [|| o ||] [|| p ||] [|| m ||]
+             [|| dr1 ||] [|| dr2 ||] )
+
+instance Ring r => Ring ( D3𝔸2 r ) where
+  !dr1 * !dr2 =
+    let
+      o :: r
+      o = fromInteger 0
+      p :: r -> r -> r
+      p = (+) @r
+      m :: r -> r -> r
+      m = (*) @r
+   in
+       $$( prodRuleQ
+             [|| o ||] [|| p ||] [|| m ||]
+             [|| dr1 ||] [|| dr2 ||] )
+
+--instance Ring r => Ring ( D1𝔸3 r ) where
+--  !dr1 * !dr2 =
+--    let
+--      o :: r
+--      o = fromInteger 0
+--      p :: r -> r -> r
+--      p = (+) @r
+--      m :: r -> r -> r
+--      m = (*) @r
+--   in
+--       $$( prodRuleQ
+--             [|| o ||] [|| p ||] [|| m ||]
+--             [|| dr1 ||] [|| dr2 ||] )
+
+instance Ring r => Ring ( D2𝔸3 r ) where
+  !dr1 * !dr2 =
+    let
+      o :: r
+      o = fromInteger 0
+      p :: r -> r -> r
+      p = (+) @r
+      m :: r -> r -> r
+      m = (*) @r
+   in
+       $$( prodRuleQ
+             [|| o ||] [|| p ||] [|| m ||]
+             [|| dr1 ||] [|| dr2 ||] )
+
+instance Ring r => Ring ( D3𝔸3 r ) where
+  !dr1 * !dr2 =
+    let
+      o :: r
+      o = fromInteger 0
+      p :: r -> r -> r
+      p = (+) @r
+      m :: r -> r -> r
+      m = (*) @r
+   in
+       $$( prodRuleQ
+             [|| o ||] [|| p ||] [|| m ||]
+             [|| dr1 ||] [|| dr2 ||] )
+
+--instance Ring r => Ring ( D1𝔸4 r ) where
+--  !dr1 * !dr2 =
+--    let
+--      o :: r
+--      o = fromInteger 0
+--      p :: r -> r -> r
+--      p = (+) @r
+--      m :: r -> r -> r
+--      m = (*) @r
+--   in
+--       $$( prodRuleQ
+--             [|| o ||] [|| p ||] [|| m ||]
+--             [|| dr1 ||] [|| dr2 ||] )
+
+instance Ring r => Ring ( D2𝔸4 r ) where
+  !dr1 * !dr2 =
+    let
+      o :: r
+      o = fromInteger 0
+      p :: r -> r -> r
+      p = (+) @r
+      m :: r -> r -> r
+      m = (*) @r
+   in
+       $$( prodRuleQ
+             [|| o ||] [|| p ||] [|| m ||]
+             [|| dr1 ||] [|| dr2 ||] )
+
+instance Ring r => Ring ( D3𝔸4 r ) where
+  !dr1 * !dr2 =
+    let
+      o :: r
+      o = fromInteger 0
+      p :: r -> r -> r
+      p = (+) @r
+      m :: r -> r -> r
+      m = (*) @r
+   in
+       $$( prodRuleQ
+             [|| o ||] [|| p ||] [|| m ||]
+             [|| dr1 ||] [|| dr2 ||] )
+
+--------------------------------------------------------------------------------
+-- Field & transcendental instances
+
+-- TODO!!
+
+deriving newtype instance Field r => Field ( D𝔸0 r )
+--instance Field r => Field ( D1𝔸1 r ) where
+--instance Field r => Field ( D1𝔸2 r ) where
+--instance Field r => Field ( D1𝔸3 r ) where
+--instance Field r => Field ( D1𝔸4 r ) where
+instance Field r => Field ( D2𝔸1 r ) where
+instance Field r => Field ( D2𝔸2 r ) where
+instance Field r => Field ( D2𝔸3 r ) where
+instance Field r => Field ( D2𝔸4 r ) where
+instance Field r => Field ( D3𝔸1 r ) where
+instance Field r => Field ( D3𝔸2 r ) where
+instance Field r => Field ( D3𝔸3 r ) where
+instance Field r => Field ( D3𝔸4 r ) where
+
+
+deriving newtype instance Transcendental r => Transcendental ( D𝔸0 r )
+--instance Transcendental r => Transcendental ( D1𝔸1 r ) where
+--instance Transcendental r => Transcendental ( D1𝔸2 r ) where
+--instance Transcendental r => Transcendental ( D1𝔸3 r ) where
+--instance Transcendental r => Transcendental ( D1𝔸4 r ) where
+instance Transcendental r => Transcendental ( D2𝔸1 r ) where
+instance Transcendental r => Transcendental ( D2𝔸2 r ) where
+instance Transcendental r => Transcendental ( D2𝔸3 r ) where
+  pi = konst @Double @2 @( ℝ 3 ) pi
+  sin ( D23 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
+      in D23 s
+            ( T $ c * dx ) ( T $ c * dy ) ( T $ c * dz )
+            ( T $ 2 * c * ddx - s * ( dx ^ 2 ) )
+            ( T $ 2 * c * dxdy - 2 * s * dx * dy )
+            ( T $ 2 * c * ddy - s * ( dy ^ 2 )  )
+            ( T $ 2 * c * dxdz - 2 * s * dx * dz )
+            ( T $ 2 * c * dydz - 2 * s * dy * dz )
+            ( T $ 2 * c * ddz - s * ( dz ^ 2 ) )
+
+  cos ( D23 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
+      in D23 c
+            ( T $ -s * dx ) ( T $ -s * dy ) ( T $ -s * dz )
+            ( T $ -2 * s * ddx - c * ( dx ^ 2 ) )
+            ( T $ -2 * s * dxdy - 2 * c * dx * dy )
+            ( T $ -2 * s * ddy - c * ( dy ^ 2 ) )
+            ( T $ -2 * s * dxdz - 2 * c * dx * dz )
+            ( T $ -2 * s * dydz - 2 * c * dy * dz )
+            ( T $ -2 * s * ddz - c * ( dz ^ 2 ) )
+
+instance Transcendental r => Transcendental ( D2𝔸4 r ) where
+instance Transcendental r => Transcendental ( D3𝔸1 r ) where
+instance Transcendental r => Transcendental ( D3𝔸2 r ) where
+instance Transcendental r => Transcendental ( D3𝔸3 r ) where
+instance Transcendental r => Transcendental ( D3𝔸4 r ) where
+
+--------------------------------------------------------------------------------
+-- HasChainRule instances.
+
+instance HasChainRule Double 2 ( ℝ 0 ) where
+  konst w = D0 w
+  linearD f v = D0 ( f v )
+  value ( D0 v ) = v
+  chain _ ( D0 gfx ) = D21 gfx origin origin
+
+instance HasChainRule Double 3 ( ℝ 0 ) where
+  konst w = D0 w
+  linearD f v = D0 ( f v )
+  value ( D0 v ) = v
+  chain _ ( D0 gfx ) = D31 gfx origin origin origin
+
+instance HasChainRule Double 2 ( ℝ 1 ) where
+
+  konst   :: forall w. AbelianGroup w => w -> D2𝔸1 w
+  konst w =
+    let !o = fromInteger @w 0
+    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+
+  value df = $$( monIndex [|| df ||] zeroMonomial )
+
+  linearD :: forall w. Module Double ( T w ) => ( ℝ 1 -> w ) -> ℝ 1 -> D2𝔸1 w
+  linearD f v =
+    let !o = origin @Double @( T w )
+    in $$( monTabulate \ mon ->
+              if | isZeroMonomial mon
+                 -> [|| f v ||]
+                 | Just i <- isLinear mon
+                 -> [|| f $$( tabulateQ \ j ->
+                                if | j == i
+                                   -> [|| 1 ||]
+                                   | otherwise
+                                   -> [|| 0 ||]
+                            ) ||]
+                 | otherwise
+                 -> [|| unT o ||]
+         )
+
+  chain :: forall w. Module Double ( T w ) => D2𝔸1 ( ℝ 1 ) -> D2𝔸1 w -> D2𝔸1 w
+  chain !df !dg =
+    let !o = origin @Double @( T w )
+        !p = (^+^)  @Double @( T w )
+        !s = (^*)   @Double @( T w )
+    in $$( chainRuleQ
+         [|| o ||] [|| p ||] [|| s ||]
+         [|| df ||] [|| dg ||] )
+
+instance HasChainRule Double 3 ( ℝ 1 ) where
+
+  konst   :: forall w. AbelianGroup w => w -> D3𝔸1 w
+  konst w =
+    let !o = fromInteger @w 0
+    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+
+  value df = $$( monIndex [|| df ||] zeroMonomial )
+
+  linearD :: forall w. Module Double ( T w ) => ( ℝ 1 -> w ) -> ℝ 1 -> D3𝔸1 w
+  linearD f v =
+    let !o = origin @Double @( T w )
+    in $$( monTabulate \ mon ->
+              if | isZeroMonomial mon
+                 -> [|| f v ||]
+                 | Just i <- isLinear mon
+                 -> [|| f $$( tabulateQ \ j ->
+                                if | j == i
+                                   -> [|| 1 ||]
+                                   | otherwise
+                                   -> [|| 0 ||]
+                            ) ||]
+                 | otherwise
+                 -> [|| unT o ||]
+         )
+
+  chain :: forall w. Module Double ( T w ) => D3𝔸1 ( ℝ 1 ) -> D3𝔸1 w -> D3𝔸1 w
+  chain !df !dg =
+    let !o = origin @Double @( T w )
+        !p = (^+^)  @Double @( T w )
+        !s = (^*)   @Double @( T w )
+    in $$( chainRuleQ
+         [|| o ||] [|| p ||] [|| s ||]
+         [|| df ||] [|| dg ||] )
+
+instance HasChainRule Double 2 ( ℝ 2 ) where
+
+  konst   :: forall w. AbelianGroup w => w -> D2𝔸2 w
+  konst w =
+    let !o = fromInteger @w 0
+    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+
+  value df = $$( monIndex [|| df ||] zeroMonomial )
+
+  linearD :: forall w. Module Double ( T w ) => ( ℝ 2 -> w ) -> ℝ 2 -> D2𝔸2 w
+  linearD f v =
+    let !o = origin @Double @( T w )
+    in $$( monTabulate \ mon ->
+              if | isZeroMonomial mon
+                 -> [|| f v ||]
+                 | Just i <- isLinear mon
+                 -> [|| f $$( tabulateQ \ j ->
+                                if | j == i
+                                   -> [|| 1 ||]
+                                   | otherwise
+                                   -> [|| 0 ||]
+                            ) ||]
+                 | otherwise
+                 -> [|| unT o ||]
+         )
+
+  chain :: forall w. Module Double ( T w ) => D2𝔸1 ( ℝ 2 ) -> D2𝔸2 w -> D2𝔸1 w
+  chain !df !dg =
+    let !o = origin @Double @( T w )
+        !p = (^+^)  @Double @( T w )
+        !s = (^*)   @Double @( T w )
+    in $$( chainRuleQ
+         [|| o ||] [|| p ||] [|| s ||]
+         [|| df ||] [|| dg ||] )
+
+instance HasChainRule Double 3 ( ℝ 2 ) where
+
+  konst   :: forall w. AbelianGroup w => w -> D3𝔸2 w
+  konst w =
+    let !o = fromInteger @w 0
+    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+
+  value df = $$( monIndex [|| df ||] zeroMonomial )
+
+  linearD :: forall w. Module Double ( T w ) => ( ℝ 2 -> w ) -> ℝ 2 -> D3𝔸2 w
+  linearD f v =
+    let !o = origin @Double @( T w )
+    in $$( monTabulate \ mon ->
+              if | isZeroMonomial mon
+                 -> [|| f v ||]
+                 | Just i <- isLinear mon
+                 -> [|| f $$( tabulateQ \ j ->
+                                if | j == i
+                                   -> [|| 1 ||]
+                                   | otherwise
+                                   -> [|| 0 ||]
+                            ) ||]
+                 | otherwise
+                 -> [|| unT o ||]
+         )
+
+  chain :: forall w. Module Double ( T w ) => D3𝔸1 ( ℝ 2 ) -> D3𝔸2 w -> D3𝔸1 w
+  chain !df !dg =
+    let !o = origin @Double @( T w )
+        !p = (^+^)  @Double @( T w )
+        !s = (^*)   @Double @( T w )
+    in $$( chainRuleQ
+         [|| o ||] [|| p ||] [|| s ||]
+         [|| df ||] [|| dg ||] )
+
+instance HasChainRule Double 2 ( ℝ 3 ) where
+
+  konst   :: forall w. AbelianGroup w => w -> D2𝔸3 w
+  konst w =
+    let !o = fromInteger @w 0
+    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+
+  value df = $$( monIndex [|| df ||] zeroMonomial )
+
+  linearD :: forall w. Module Double ( T w ) => ( ℝ 3 -> w ) -> ℝ 3 -> D2𝔸3 w
+  linearD f v =
+    let !o = origin @Double @( T w )
+    in $$( monTabulate \ mon ->
+              if | isZeroMonomial mon
+                 -> [|| f v ||]
+                 | Just i <- isLinear mon
+                 -> [|| f $$( tabulateQ \ j ->
+                                if | j == i
+                                   -> [|| 1 ||]
+                                   | otherwise
+                                   -> [|| 0 ||]
+                            ) ||]
+                 | otherwise
+                 -> [|| unT o ||]
+         )
+
+  chain :: forall w. Module Double ( T w ) => D2𝔸1 ( ℝ 3 ) -> D2𝔸3 w -> D2𝔸1 w
+  chain !df !dg =
+    let !o = origin @Double @( T w )
+        !p = (^+^)  @Double @( T w )
+        !s = (^*)   @Double @( T w )
+    in $$( chainRuleQ
+         [|| o ||] [|| p ||] [|| s ||]
+         [|| df ||] [|| dg ||] )
+
+instance HasChainRule Double 3 ( ℝ 3 ) where
+
+  konst   :: forall w. AbelianGroup w => w -> D3𝔸3 w
+  konst w =
+    let !o = fromInteger @w 0
+    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+
+  value df = $$( monIndex [|| df ||] zeroMonomial )
+
+  linearD :: forall w. Module Double ( T w ) => ( ℝ 3 -> w ) -> ℝ 3 -> D3𝔸3 w
+  linearD f v =
+    let !o = origin @Double @( T w )
+    in $$( monTabulate \ mon ->
+              if | isZeroMonomial mon
+                 -> [|| f v ||]
+                 | Just i <- isLinear mon
+                 -> [|| f $$( tabulateQ \ j ->
+                                if | j == i
+                                   -> [|| 1 ||]
+                                   | otherwise
+                                   -> [|| 0 ||]
+                            ) ||]
+                 | otherwise
+                 -> [|| unT o ||]
+         )
+
+  chain :: forall w. Module Double ( T w ) => D3𝔸1 ( ℝ 3 ) -> D3𝔸3 w -> D3𝔸1 w
+  chain !df !dg =
+    let !o = origin @Double @( T w )
+        !p = (^+^)  @Double @( T w )
+        !s = (^*)   @Double @( T w )
+    in $$( chainRuleQ
+         [|| o ||] [|| p ||] [|| s ||]
+         [|| df ||] [|| dg ||] )
+
+instance HasChainRule Double 2 ( ℝ 4 ) where
+
+  konst   :: forall w. AbelianGroup w => w -> D2𝔸4 w
+  konst w =
+    let !o = fromInteger @w 0
+    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+
+  value df = $$( monIndex [|| df ||] zeroMonomial )
+
+  linearD :: forall w. Module Double ( T w ) => ( ℝ 4 -> w ) -> ℝ 4 -> D2𝔸4 w
+  linearD f v =
+    let !o = origin @Double @( T w )
+    in $$( monTabulate \ mon ->
+              if | isZeroMonomial mon
+                 -> [|| f v ||]
+                 | Just i <- isLinear mon
+                 -> [|| f $$( tabulateQ \ j ->
+                                if | j == i
+                                   -> [|| 1 ||]
+                                   | otherwise
+                                   -> [|| 0 ||]
+                            ) ||]
+                 | otherwise
+                 -> [|| unT o ||]
+         )
+
+  chain :: forall w. Module Double ( T w ) => D2𝔸1 ( ℝ 4 ) -> D2𝔸4 w -> D2𝔸1 w
+  chain !df !dg =
+    let !o = origin @Double @( T w )
+        !p = (^+^)  @Double @( T w )
+        !s = (^*)   @Double @( T w )
+    in $$( chainRuleQ
+         [|| o ||] [|| p ||] [|| s ||]
+         [|| df ||] [|| dg ||] )
+
+instance HasChainRule Double 3 ( ℝ 4 ) where
+
+  konst   :: forall w. AbelianGroup w => w -> D3𝔸4 w
+  konst w =
+    let !o = fromInteger @w 0
+    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+
+  value df = $$( monIndex [|| df ||] zeroMonomial )
+
+  linearD :: forall w. Module Double ( T w ) => ( ℝ 4 -> w ) -> ℝ 4 -> D3𝔸4 w
+  linearD f v =
+    let !o = origin @Double @( T w )
+    in $$( monTabulate \ mon ->
+              if | isZeroMonomial mon
+                 -> [|| f v ||]
+                 | Just i <- isLinear mon
+                 -> [|| f $$( tabulateQ \ j ->
+                                if | j == i
+                                   -> [|| 1 ||]
+                                   | otherwise
+                                   -> [|| 0 ||]
+                            ) ||]
+                 | otherwise
+                 -> [|| unT o ||]
+         )
+
+  chain :: forall w. Module Double ( T w ) => D3𝔸1 ( ℝ 4 ) -> D3𝔸4 w -> D3𝔸1 w
+  chain !df !dg =
+    let !o = origin @Double @( T w )
+        !p = (^+^)  @Double @( T w )
+        !s = (^*)   @Double @( T w )
+    in $$( chainRuleQ
+         [|| o ||] [|| p ||] [|| s ||]
+         [|| df ||] [|| dg ||] )

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

@@ -0,0 +1,588 @@
+
+{-# LANGUAGE AllowAmbiguousTypes   #-}
+{-# LANGUAGE QuantifiedConstraints #-}
+{-# LANGUAGE ScopedTypeVariables   #-}
+{-# LANGUAGE TemplateHaskell       #-}
+{-# LANGUAGE UndecidableInstances  #-}
+
+module Math.Algebra.Dual.Internal
+  ( D𝔸0(..)
+  , D1𝔸1(..), D2𝔸1(..), D3𝔸1(..)
+  , D1𝔸2(..), D2𝔸2(..), D3𝔸2(..)
+  , D1𝔸3(..), D2𝔸3(..), D3𝔸3(..)
+  , D1𝔸4(..), D2𝔸4(..), D3𝔸4(..)
+
+  , chainRuleQ
+  ) where
+
+-- base
+import GHC.Generics
+  ( Generic, Generic1, Generically1(..) )
+import GHC.TypeNats
+  ( KnownNat )
+
+-- template-haskell
+import Language.Haskell.TH
+  ( CodeQ )
+import Language.Haskell.TH.Syntax
+  ( Lift(..) )
+
+-- MetaBrush
+import Math.Linear
+  ( Vec(..), T(..)
+  , RepresentableQ(..), RepDim
+  , zipIndices
+  )
+import Math.Monomial
+  ( Mon(..), MonomialBasis(..), Vars, Deg
+  , mons, faà, multiSubsetsSum, zeroMonomial
+  )
+import Math.Ring
+  ( Ring )
+import qualified Math.Ring as Ring
+import TH.Utils
+  ( foldQ, powQ )
+
+--------------------------------------------------------------------------------
+
+-- | \( \mathbb{Z} \).
+newtype D𝔸0 v = D0 { _D0_v :: v }
+  deriving stock   ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving Applicative
+    via Generically1 D𝔸0
+
+-- | \( \mathbb{Z}[x]/(x)^2 \).
+data D1𝔸1 v =
+  D11 { _D11_v :: !v
+      , _D11_dx :: !( T v )
+      }
+  deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving Applicative
+    via Generically1 D1𝔸1
+
+-- | \( \mathbb{Z}[x]/(x)^3 \).
+data D2𝔸1 v =
+  D21 { _D21_v :: !v
+      , _D21_dx :: !( T v )
+      , _D21_dxdx :: !( T v )
+      }
+  deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving Applicative
+    via Generically1 D2𝔸1
+
+-- | \( \mathbb{Z}[x]/(x)^4 \).
+data D3𝔸1 v =
+  D31 { _D31_v :: !v
+      , _D31_dx :: !( T v )
+      , _D31_dxdx :: !( T v )
+      , _D31_dxdxdx :: !( T v )
+      }
+  deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving Applicative
+    via Generically1 D3𝔸1
+
+-- | \( \mathbb{Z}[x, y]/(x, y)^2 \).
+data D1𝔸2 v =
+  D12 { _D12_v :: !v
+      , _D12_dx, _D12_dy :: !( T v )
+      }
+  deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving Applicative
+    via Generically1 D1𝔸2
+
+-- | \( \mathbb{Z}[x, y]/(x, y)^3 \).
+data D2𝔸2 v =
+  D22 { _D22_v :: !v
+      , _D22_dx, _D22_dy :: !( T v )
+      , _D22_dxdx, _D22_dxdy, _D22_dydy :: !( T v )
+      }
+  deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving Applicative
+    via Generically1 D2𝔸2
+
+-- | \( \mathbb{Z}[x, y]/(x, y)^4 \).
+data D3𝔸2 v =
+  D32 { _D32_v :: !v
+      , _D32_dx, _D32_dy :: !( T v )
+      , _D32_dxdx, _D32_dxdy, _D32_dydy :: !( T v )
+      , _D32_dxdxdx, _D32_dxdxdy, _D32_dxdydy, _D32_dydydy :: !( T v )
+      }
+  deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving Applicative
+    via Generically1 D3𝔸2
+
+-- | \( \mathbb{Z}[x, y, z]/(x, y, z)^2 \).
+data D1𝔸3 v =
+  D13 { _D13_v :: !v
+      , _D13_dx, _D13_dy, _D13_dz :: !( T v )
+      }
+  deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving Applicative
+    via Generically1 D1𝔸3
+
+-- | \( \mathbb{Z}[x, y, z]/(x, y, z)^3 \).
+data D2𝔸3 v =
+  D23 { _D23_v :: !v
+      , _D23_dx, _D23_dy, _D23_dz :: !( T v )
+      , _D23_dxdx, _D23_dxdy, _D23_dydy, _D23_dxdz, _D23_dydz, _D23_dzdz :: !( T v )
+      }
+  deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving Applicative
+    via Generically1 D2𝔸3
+
+-- | \( \mathbb{Z}[x, y, z]/(x, y, z)^4 \).
+data D3𝔸3 v =
+  D33 { _D33_v :: !v
+      , _D33_dx, _D33_dy, _D33_dz :: !( T v )
+      , _D33_dxdx, _D33_dxdy, _D33_dydy, _D33_dxdz, _D33_dydz, _D33_dzdz :: !( T v )
+      , _D33_dxdxdx, _D33_dxdxdy, _D33_dxdydy, _D33_dydydy
+         , _D33_dxdxdz, _D33_dxdydz, _D33_dxdzdz, _D33_dydydz, _D33_dydzdz, _D33_dzdzdz :: !( T v )
+      }
+  deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving Applicative
+    via Generically1 D3𝔸3
+
+-- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^2 \).
+data D1𝔸4 v =
+  D14 { _D14_v :: !v
+      , _D14_dx, _D14_dy, _D14_dz, _D14_dw :: !( T v )
+      }
+  deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving Applicative
+    via Generically1 D1𝔸4
+
+-- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^3 \).
+data D2𝔸4 v =
+  D24 { _D24_v :: !v
+      , _D24_dx, _D24_dy, _D24_dz, _D24_dw :: !( T v )
+      , _D24_dxdx, _D24_dxdy, _D24_dydy, _D24_dxdz
+         , _D24_dydz, _D24_dzdz, _D24_dxdw, _D24_dydw, _D24_dzdw, _D24_dwdw :: !( T v )
+      }
+  deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving Applicative
+    via Generically1 D2𝔸4
+
+-- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^3 \).
+data D3𝔸4 v =
+  D34 { _D34_v :: !v
+      , _D34_dx, _D34_dy, _D34_dz, _D34_dw :: !( T v )
+      , _D34_dxdx, _D34_dxdy, _D34_dydy, _D34_dxdz, _D34_dydz, _D34_dzdz
+        , _D34_dxdw, _D34_dydw, _D34_dzdw, _D34_dwdw :: !( T v )
+      , _D34_dxdxdx, _D34_dxdxdy, _D34_dxdydy, _D34_dydydy,
+         _D34_dxdxdz, _D34_dxdydz, _D34_dxdzdz, _D34_dydzdz, _D34_dydydz, _D34_dzdzdz,
+         _D34_dxdxdw, _D34_dxdydw, _D34_dydydw, _D34_dxdzdw, _D34_dydzdw, _D34_dzdzdw,
+         _D34_dxdwdw, _D34_dydwdw, _D34_dzdwdw, _D34_dwdwdw :: !( T v )
+      }
+  deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving Applicative
+    via Generically1 D3𝔸4
+
+--------------------------------------------------------------------------------
+
+-- | The chain rule, to be spliced in using Template Haskell.
+chainRuleQ :: forall dr1 dv v r w d
+           .  ( Ring r, RepresentableQ r v
+              , MonomialBasis dr1, Vars dr1 ~ 1
+              , MonomialBasis dv , Vars dv  ~ RepDim v
+              , Deg dr1 ~ Deg dv
+              , d ~ Vars dv, KnownNat d
+              )
+           => CodeQ ( T w )               -- Module r ( T w )
+           -> CodeQ ( T w -> T w -> T w ) --
+           -> CodeQ ( T w -> r -> T w )   -- (circumvent TH constraint issue)
+           -> CodeQ ( dr1 v )
+           -> CodeQ ( dv  w )
+           -> CodeQ ( dr1 w )
+chainRuleQ zero_w sum_w scale_w df dg =
+    monTabulate @dr1 \ ( Mon ( k `VS` _ ) ) ->
+      case k of
+        -- Set the value of the composition separately,
+        -- as that isn't handled by the Faà di Bruno formula.
+        0 -> monIndex @dv dg zeroMonomial
+        _ ->
+          [|| unT $ $$( foldQ sum_w zero_w
+                  [ [|| $$scale_w ( T $$( monIndex @dv dg m_g ) )
+                          $$( foldQ [|| (Ring.+) ||] [|| Ring.fromInteger ( 0 :: Integer ) ||]
+                                [ [|| Ring.fromInteger $$( liftTyped $ fromIntegral $ faà k is ) Ring.*
+                                     $$( foldQ [|| (Ring.*) ||] [|| Ring.fromInteger ( 1 :: Integer ) ||]
+                                        [ foldQ [|| (Ring.*) ||] [|| Ring.fromInteger ( 1 :: Integer ) ||]
+                                          [ ( indexQ @r @v ( monIndex @dr1 df ( Mon $ f_deg `VS` VZ ) ) v_index )
+                                                `powQ` pow
+                                          | ( f_deg, pow ) <- pows
+                                          ]
+                                        | ( v_index, pows ) <- zipIndices is
+                                        ]
+                                      ) ||]
+                                | is <- mss
+                                ]
+                           )
+                      ||]
+                  | m_g <- mons @_ @d k
+                  , let mss = multiSubsetsSum [ 1 .. k ] k ( monDegs m_g )
+                  , not ( null mss ) -- avoid terms of the form x * 0
+                  ]
+              ) ||]
+
+--------------------------------------------------------------------------------
+-- MonomialBasis instances follow (nothing else).
+
+type instance Deg  D𝔸0 = 0
+type instance Vars D𝔸0 = 0
+instance MonomialBasis D𝔸0 where
+  monTabulate f =
+    [|| let
+          !_D0_v = $$( f $ Mon VZ )
+        in D0 { .. }
+      ||]
+
+  monIndex d _ = [|| _D0_v $$d ||]
+
+type instance Deg  D1𝔸1 = 1
+type instance Vars D1𝔸1 = 1
+instance MonomialBasis D1𝔸1 where
+  monTabulate f =
+    [|| let
+          !_D11_v    =   $$( f $ Mon ( 0 `VS` VZ ) )
+          !_D11_dx   = T $$( f $ Mon ( 1 `VS` VZ ) )
+        in D11 { .. }
+      ||]
+
+  monIndex d = \ case
+    Mon ( 1 `VS` VZ ) -> [|| unT $ _D11_dx   $$d ||]
+    _                 -> [||       _D11_v    $$d ||]
+
+type instance Deg  D2𝔸1 = 2
+type instance Vars D2𝔸1 = 1
+instance MonomialBasis D2𝔸1 where
+  monTabulate f =
+    [|| let
+          !_D21_v    =   $$( f $ Mon ( 0 `VS` VZ ) )
+          !_D21_dx   = T $$( f $ Mon ( 1 `VS` VZ ) )
+          !_D21_dxdx = T $$( f $ Mon ( 2 `VS` VZ ) )
+        in D21 { .. }
+      ||]
+
+  monIndex d = \ case
+    Mon ( 1 `VS` VZ ) -> [|| unT $ _D21_dx   $$d ||]
+    Mon ( 2 `VS` VZ ) -> [|| unT $ _D21_dxdx $$d ||]
+    _                 -> [||       _D21_v    $$d ||]
+
+type instance Deg  D3𝔸1 = 3
+type instance Vars D3𝔸1 = 1
+instance MonomialBasis D3𝔸1 where
+  monTabulate f =
+    [|| let
+          !_D31_v      =   $$( f $ Mon ( 0 `VS` VZ ) )
+          !_D31_dx     = T $$( f $ Mon ( 1 `VS` VZ ) )
+          !_D31_dxdx   = T $$( f $ Mon ( 2 `VS` VZ ) )
+          !_D31_dxdxdx = T $$( f $ Mon ( 3 `VS` VZ ) )
+        in D31 { .. }
+      ||]
+
+  monIndex d = \ case
+    Mon ( 1 `VS` VZ ) -> [|| unT $ _D31_dx     $$d ||]
+    Mon ( 2 `VS` VZ ) -> [|| unT $ _D31_dxdx   $$d ||]
+    Mon ( 3 `VS` VZ ) -> [|| unT $ _D31_dxdxdx $$d ||]
+    _                 -> [||       _D31_v      $$d ||]
+
+type instance Deg  D1𝔸2 = 1
+type instance Vars D1𝔸2 = 2
+instance MonomialBasis D1𝔸2 where
+  monTabulate f =
+    [|| let
+          !_D12_v    =   $$( f $ Mon ( 0 `VS` 0 `VS` VZ ) )
+          !_D12_dx   = T $$( f $ Mon ( 1 `VS` 0 `VS` VZ ) )
+          !_D12_dy   = T $$( f $ Mon ( 0 `VS` 1 `VS` VZ ) )
+        in D12 { .. }
+      ||]
+
+  monIndex d = \ case
+    Mon ( 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D12_dx $$d ||]
+    Mon ( 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D12_dy $$d ||]
+    _                        -> [||       _D12_v  $$d ||]
+
+type instance Deg  D2𝔸2 = 2
+type instance Vars D2𝔸2 = 2
+instance MonomialBasis D2𝔸2 where
+  monTabulate f =
+    [|| let
+          !_D22_v    =   $$( f $ Mon ( 0 `VS` 0 `VS` VZ ) )
+          !_D22_dx   = T $$( f $ Mon ( 1 `VS` 0 `VS` VZ ) )
+          !_D22_dy   = T $$( f $ Mon ( 0 `VS` 1 `VS` VZ ) )
+          !_D22_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` VZ ) )
+          !_D22_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` VZ ) )
+          !_D22_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` VZ ) )
+        in D22 { .. }
+      ||]
+
+  monIndex d = \ case
+    Mon ( 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D22_dx   $$d ||]
+    Mon ( 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D22_dy   $$d ||]
+    Mon ( 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D22_dxdx $$d ||]
+    Mon ( 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D22_dxdy $$d ||]
+    Mon ( 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D22_dydy $$d ||]
+    _                        -> [||       _D22_v    $$d ||]
+
+type instance Deg  D3𝔸2 = 3
+type instance Vars D3𝔸2 = 2
+instance MonomialBasis D3𝔸2 where
+  monTabulate f =
+    [|| let
+          !_D32_v      =   $$( f $ Mon ( 0 `VS` 0 `VS` VZ ) )
+          !_D32_dx     = T $$( f $ Mon ( 1 `VS` 0 `VS` VZ ) )
+          !_D32_dy     = T $$( f $ Mon ( 0 `VS` 1 `VS` VZ ) )
+          !_D32_dxdx   = T $$( f $ Mon ( 2 `VS` 0 `VS` VZ ) )
+          !_D32_dxdy   = T $$( f $ Mon ( 1 `VS` 1 `VS` VZ ) )
+          !_D32_dydy   = T $$( f $ Mon ( 0 `VS` 2 `VS` VZ ) )
+          !_D32_dxdxdx = T $$( f $ Mon ( 3 `VS` 0 `VS` VZ ) )
+          !_D32_dxdxdy = T $$( f $ Mon ( 2 `VS` 1 `VS` VZ ) )
+          !_D32_dxdydy = T $$( f $ Mon ( 1 `VS` 2 `VS` VZ ) )
+          !_D32_dydydy = T $$( f $ Mon ( 0 `VS` 3 `VS` VZ ) )
+       in D32 { .. } ||]
+
+  monIndex d = \ case
+    Mon ( 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D32_dx     $$d ||]
+    Mon ( 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D32_dy     $$d ||]
+    Mon ( 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D32_dxdx   $$d ||]
+    Mon ( 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D32_dxdy   $$d ||]
+    Mon ( 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D32_dydy   $$d ||]
+    Mon ( 3 `VS` 0 `VS` VZ ) -> [|| unT $ _D32_dxdxdx $$d ||]
+    Mon ( 2 `VS` 1 `VS` VZ ) -> [|| unT $ _D32_dxdxdy $$d ||]
+    Mon ( 1 `VS` 2 `VS` VZ ) -> [|| unT $ _D32_dxdydy $$d ||]
+    Mon ( 0 `VS` 3 `VS` VZ ) -> [|| unT $ _D32_dydydy $$d ||]
+    _                        -> [||       _D32_v      $$d ||]
+
+type instance Deg  D1𝔸3 = 1
+type instance Vars D1𝔸3 = 3
+instance MonomialBasis D1𝔸3 where
+  monTabulate f =
+    [|| let
+          !_D13_v    =   $$( f ( Mon ( 0 `VS` 0 `VS` 0 `VS` VZ ) ) )
+          !_D13_dx   = T $$( f ( Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) ) )
+          !_D13_dy   = T $$( f ( Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) ) )
+          !_D13_dz   = T $$( f ( Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) ) )
+        in D13 { .. }
+      ||]
+
+  monIndex d = \ case
+    Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D13_dx $$d ||]
+    Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D13_dy $$d ||]
+    Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D13_dz $$d ||]
+    _                               -> [||       _D13_v  $$d ||]
+
+type instance Deg  D2𝔸3 = 2
+type instance Vars D2𝔸3 = 3
+instance MonomialBasis D2𝔸3 where
+  monTabulate f =
+    [|| let
+          !_D23_v    =   $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D23_dx   = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D23_dy   = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) )
+          !_D23_dz   = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) )
+          !_D23_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D23_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) )
+          !_D23_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) )
+          !_D23_dxdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) )
+          !_D23_dydz = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) )
+          !_D23_dzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) )
+        in D23 { .. }
+      ||]
+
+  monIndex d = \ case
+    Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D23_dx   $$d ||]
+    Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D23_dy   $$d ||]
+    Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D23_dz   $$d ||]
+    Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D23_dxdx $$d ||]
+    Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D23_dxdy $$d ||]
+    Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D23_dydy $$d ||]
+    Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D23_dxdz $$d ||]
+    Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D23_dydz $$d ||]
+    Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D23_dzdz $$d ||]
+    _                               -> [||       _D23_v    $$d ||]
+
+
+type instance Deg  D3𝔸3 = 3
+type instance Vars D3𝔸3 = 3
+instance MonomialBasis D3𝔸3 where
+  monTabulate f =
+    [|| let
+          !_D33_v      =   $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D33_dx     = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D33_dy     = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) )
+          !_D33_dz     = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) )
+          !_D33_dxdx   = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D33_dxdy   = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) )
+          !_D33_dydy   = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) )
+          !_D33_dxdz   = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) )
+          !_D33_dydz   = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) )
+          !_D33_dzdz   = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) )
+          !_D33_dxdxdx = T $$( f $ Mon ( 3 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D33_dxdxdy = T $$( f $ Mon ( 2 `VS` 1 `VS` 0 `VS` VZ ) )
+          !_D33_dxdydy = T $$( f $ Mon ( 1 `VS` 2 `VS` 0 `VS` VZ ) )
+          !_D33_dydydy = T $$( f $ Mon ( 0 `VS` 3 `VS` 0 `VS` VZ ) )
+          !_D33_dxdxdz = T $$( f $ Mon ( 2 `VS` 0 `VS` 1 `VS` VZ ) )
+          !_D33_dxdydz = T $$( f $ Mon ( 1 `VS` 1 `VS` 1 `VS` VZ ) )
+          !_D33_dxdzdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 2 `VS` VZ ) )
+          !_D33_dydydz = T $$( f $ Mon ( 0 `VS` 2 `VS` 1 `VS` VZ ) )
+          !_D33_dydzdz = T $$( f $ Mon ( 0 `VS` 1 `VS` 2 `VS` VZ ) )
+          !_D33_dzdzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 3 `VS` VZ ) )
+       in D33 { .. } ||]
+
+  monIndex d = \ case
+    Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dx     $$d ||]
+    Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dy     $$d ||]
+    Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dz     $$d ||]
+    Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dxdx   $$d ||]
+    Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dxdy   $$d ||]
+    Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dydy   $$d ||]
+    Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dxdz   $$d ||]
+    Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dydz   $$d ||]
+    Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D33_dzdz   $$d ||]
+    Mon ( 3 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dxdxdx $$d ||]
+    Mon ( 2 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dxdxdy $$d ||]
+    Mon ( 1 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dxdydy $$d ||]
+    Mon ( 0 `VS` 3 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dydydy $$d ||]
+    Mon ( 2 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dxdxdz $$d ||]
+    Mon ( 1 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dxdydz $$d ||]
+    Mon ( 1 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D33_dxdzdz $$d ||]
+    Mon ( 0 `VS` 1 `VS` 2 `VS` VZ ) -> [|| unT $ _D33_dydzdz $$d ||]
+    Mon ( 0 `VS` 0 `VS` 3 `VS` VZ ) -> [|| unT $ _D33_dzdzdz $$d ||]
+    _                               -> [||       _D33_v      $$d ||]
+
+type instance Deg  D1𝔸4 = 1
+type instance Vars D1𝔸4 = 4
+instance MonomialBasis D1𝔸4 where
+  monTabulate f =
+    [|| let
+          !_D14_v    =   $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D14_dx   = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D14_dy   = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D14_dz   = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) )
+          !_D14_dw   = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) )
+        in D14 { .. } ||]
+
+  monIndex d = \ case
+    Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D14_dx   $$d ||]
+    Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D14_dy   $$d ||]
+    Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D14_dz   $$d ||]
+    Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D14_dw   $$d ||]
+    _                                      -> [||       _D14_v    $$d ||]
+
+type instance Deg  D2𝔸4 = 2
+type instance Vars D2𝔸4 = 4
+instance MonomialBasis D2𝔸4 where
+  monTabulate f =
+    [|| let
+          !_D24_v    =   $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D24_dx   = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D24_dy   = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D24_dz   = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) )
+          !_D24_dw   = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) )
+          !_D24_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D24_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D24_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D24_dxdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) )
+          !_D24_dydz = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) )
+          !_D24_dzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) )
+          !_D24_dxdw = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) )
+          !_D24_dydw = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) )
+          !_D24_dzdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) )
+          !_D24_dwdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) )
+        in D24 { .. } ||]
+
+  monIndex d = \ case
+    Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dx   $$d ||]
+    Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dy   $$d ||]
+    Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dz   $$d ||]
+    Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D24_dw   $$d ||]
+    Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dxdx $$d ||]
+    Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dxdy $$d ||]
+    Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dydy $$d ||]
+    Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dxdz $$d ||]
+    Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dydz $$d ||]
+    Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dzdz $$d ||]
+    Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D24_dxdw $$d ||]
+    Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D24_dydw $$d ||]
+    Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D24_dzdw $$d ||]
+    Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D24_dwdw $$d ||]
+    _                                      -> [||       _D24_v    $$d ||]
+
+
+type instance Deg  D3𝔸4 = 3
+type instance Vars D3𝔸4 = 4
+instance MonomialBasis D3𝔸4 where
+  monTabulate f =
+    [|| let
+          !_D34_v      =   $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D34_dx     = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D34_dy     = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D34_dz     = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) )
+          !_D34_dw     = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) )
+          !_D34_dxdx   = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D34_dxdy   = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D34_dydy   = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D34_dxdz   = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) )
+          !_D34_dydz   = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) )
+          !_D34_dzdz   = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) )
+          !_D34_dxdw   = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) )
+          !_D34_dydw   = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) )
+          !_D34_dzdw   = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) )
+          !_D34_dwdw   = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) )
+          !_D34_dxdxdx = T $$( f $ Mon ( 3 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D34_dxdxdy = T $$( f $ Mon ( 2 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D34_dxdydy = T $$( f $ Mon ( 1 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D34_dydydy = T $$( f $ Mon ( 0 `VS` 3 `VS` 0 `VS` 0 `VS` VZ ) )
+          !_D34_dxdxdz = T $$( f $ Mon ( 2 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) )
+          !_D34_dxdydz = T $$( f $ Mon ( 1 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) )
+          !_D34_dxdzdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) )
+          !_D34_dydydz = T $$( f $ Mon ( 0 `VS` 2 `VS` 1 `VS` 0 `VS` VZ ) )
+          !_D34_dydzdz = T $$( f $ Mon ( 0 `VS` 1 `VS` 2 `VS` 0 `VS` VZ ) )
+          !_D34_dzdzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 3 `VS` 0 `VS` VZ ) )
+          !_D34_dxdxdw = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) )
+          !_D34_dxdydw = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) )
+          !_D34_dydydw = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 1 `VS` VZ ) )
+          !_D34_dxdzdw = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) )
+          !_D34_dydzdw = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 1 `VS` VZ ) )
+          !_D34_dzdzdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 1 `VS` VZ ) )
+          !_D34_dxdwdw = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) )
+          !_D34_dydwdw = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 2 `VS` VZ ) )
+          !_D34_dzdwdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 2 `VS` VZ ) )
+          !_D34_dwdwdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 3 `VS` VZ ) )
+        in D34 { .. } ||]
+
+  monIndex d = \ case
+    Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dx     $$d ||]
+    Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dy     $$d ||]
+    Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dz     $$d ||]
+    Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dw     $$d ||]
+    Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdx   $$d ||]
+    Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdy   $$d ||]
+    Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dydy   $$d ||]
+    Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdz   $$d ||]
+    Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dydz   $$d ||]
+    Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dzdz   $$d ||]
+    Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dxdw   $$d ||]
+    Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dydw   $$d ||]
+    Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dzdw   $$d ||]
+    Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D34_dwdw   $$d ||]
+    Mon ( 3 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdxdx $$d ||]
+    Mon ( 2 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdxdy $$d ||]
+    Mon ( 1 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdydy $$d ||]
+    Mon ( 0 `VS` 3 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dydydy $$d ||]
+    Mon ( 2 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdxdz $$d ||]
+    Mon ( 1 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdydz $$d ||]
+    Mon ( 1 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdzdz $$d ||]
+    Mon ( 0 `VS` 2 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dydydz $$d ||]
+    Mon ( 0 `VS` 1 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dydzdz $$d ||]
+    Mon ( 0 `VS` 0 `VS` 3 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dzdzdz $$d ||]
+    Mon ( 2 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dxdxdw $$d ||]
+    Mon ( 1 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dxdydw $$d ||]
+    Mon ( 0 `VS` 2 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dydydw $$d ||]
+    Mon ( 1 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dxdzdw $$d ||]
+    Mon ( 0 `VS` 1 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dydzdw $$d ||]
+    Mon ( 0 `VS` 0 `VS` 2 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dzdzdw $$d ||]
+    Mon ( 1 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D34_dxdwdw $$d ||]
+    Mon ( 0 `VS` 1 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D34_dydwdw $$d ||]
+    Mon ( 0 `VS` 0 `VS` 1 `VS` 2 `VS` VZ ) -> [|| unT $ _D34_dzdwdw $$d ||]
+    Mon ( 0 `VS` 0 `VS` 0 `VS` 3 `VS` VZ ) -> [|| unT $ _D34_dwdwdw $$d ||]
+    _                                      -> [||       _D34_v      $$d ||]

src/splines/Math/Bezier/Cubic.hs

@@ -63,6 +63,7 @@ import Math.Roots
   ( realRoots, solveQuadratic )
 import Math.Linear
   ( ℝ(..), T(..) )
+import qualified Math.Ring as Ring
 
 --------------------------------------------------------------------------------
 
@@ -107,13 +108,13 @@ bezier ( Bezier {..} ) t =
 -- | Derivative of a cubic Bézier curve.
 bezier' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> v
 bezier' ( Bezier {..} )
-  = ( 3 *^ )
+  = ( Ring.fromInteger 3 *^ )
   . Quadratic.bezier @v ( Quadratic.Bezier ( p0 --> p1 ) ( p1 --> p2 ) ( p2 --> p3 ) )
 
 -- | Second derivative of a cubic Bézier curve.
 bezier'' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> v
 bezier'' ( Bezier {..} ) t
-  = ( 6 *^ )
+  = ( Ring.fromInteger 6 *^ )
   $ lerp @v t
       ( p1 --> p0 ^+^ p1 --> p2 )
       ( p2 --> p1 ^+^ p2 --> p3 )

src/splines/Math/Bezier/Quadratic.hs

@@ -60,6 +60,7 @@ import Math.Roots
   ( realRoots )
 import Math.Linear
   ( ℝ(..), T(..) )
+import qualified Math.Ring as Ring
 
 --------------------------------------------------------------------------------
 
@@ -91,11 +92,11 @@ bezier ( Bezier {..} ) t = lerp @v t ( lerp @v t p0 p1 ) ( lerp @v t p1 p2 )
 
 -- | Derivative of a quadratic Bézier curve.
 bezier' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> v
-bezier' ( Bezier {..} ) t = 2 *^ lerp @v t ( p0 --> p1 ) ( p1 --> p2 )
+bezier' ( Bezier {..} ) t = Ring.fromInteger 2 *^ lerp @v t ( p0 --> p1 ) ( p1 --> p2 )
 
 -- | Second derivative of a quadratic Bézier curve.
 bezier'' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> v
-bezier'' ( Bezier {..} ) = 2 *^ ( p1 --> p0 ^+^ p1 --> p2 )
+bezier'' ( Bezier {..} ) = Ring.fromInteger 2 *^ ( p1 --> p0 ^+^ p1 --> p2 )
 
 -- | Curvature of a quadratic Bézier curve.
 curvature :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> r -> r

src/splines/Math/Bezier/Stroke.hs

@@ -100,6 +100,7 @@ import Control.Monad.Trans.Writer.CPS
   ( WriterT, execWriterT, runWriter, tell )
 
 -- MetaBrush
+import Math.Algebra.Dual
 import qualified Math.Bezier.Cubic     as Cubic
 import Math.Bezier.Cubic.Fit
   ( FitPoint, FitParameters, fitSpline )
@@ -114,8 +115,12 @@ import Math.Bezier.Spline
   , showSplinePoints
   )
 import qualified Math.Bezier.Quadratic as Quadratic
+import Math.Differentiable
+  ( Differentiable, DiffInterp )
 import Math.Epsilon
   ( epsilon )
+import Math.Interval
+import Math.Linear
 import Math.Module
   ( Module(..), Inner((^.^)), Cross(cross), Interpolatable
   , lerp, convexCombination, strictlyParallel
@@ -126,8 +131,7 @@ import Math.Orientation
   )
 import Math.Roots
   ( solveQuadratic, newtonRaphson )
-import Math.Linear
+
-import Math.Linear.Dual
 
 import Debug.Utils
 
@@ -195,12 +199,15 @@ computeStrokeOutline ::
      , NFData ptData, NFData crvData
 
      -- Differentiability.
-     , Differentiable 'Point    brushParams
+     , DiffInterp 'Point    brushParams
-     , Differentiable 'Interval brushParams
+     , DiffInterp 'Interval brushParams
      , Interpolatable Double usedParams
      , Interpolatable ( 𝕀 Double ) ( 𝕀 usedParams )
-     , Diffy Double usedParams
+     , HasChainRule Double 2 brushParams
-     , Diffy ( 𝕀 Double ) ( 𝕀 usedParams )
+     , HasChainRule ( 𝕀 Double ) 2 ( 𝕀 brushParams )
+     , HasChainRule Double 2 usedParams
+     , HasChainRule ( 𝕀 Double ) 2 ( 𝕀 usedParams )
+     , Traversable ( D 2 brushParams )
 
      -- Debugging.
      , Show ptData, Show brushParams
@@ -209,7 +216,8 @@ computeStrokeOutline ::
   => FitParameters
   -> ( ptData -> usedParams )
   -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing
-  -> ( forall i. Differentiable i brushParams
+  -> ( forall i
+       .  DiffInterp i brushParams
        => Proxy# i
        -> ( forall a. a -> I i a )
        -> I i brushParams ~> Spline Closed () ( I i ( ℝ 2 ) )
@@ -323,7 +331,7 @@ computeStrokeOutline fitParams ptParams toBrushParams brushFn spline@( Spline {
           outlineFunction ptParams toBrushParams brushFn p0 crv :<| go ( openCurveEnd crv ) crvs
 
     brushShape :: ptData -> SplinePts Closed
-    brushShape pt = fun ( brushFn @Point proxy# id ) $ toBrushParams $ ptParams pt
+    brushShape pt = fun @Double ( brushFn @Point proxy# id ) $ toBrushParams $ ptParams pt
 
     updateSpline :: ( T ( ℝ 2 ), T ( ℝ 2 ), T ( ℝ 2 ) ) -> ST s OutlineData
     updateSpline ( lastTgt, lastTgtFwd, lastTgtBwd )
@@ -438,19 +446,26 @@ outlineFunction
   .  ( HasType ( ℝ 2 ) ptData
 
      -- Differentiability.
-     , Differentiable 'Point    brushParams
-     , Differentiable 'Interval brushParams
      , Interpolatable Double usedParams
      , Interpolatable ( 𝕀 Double ) ( 𝕀 usedParams )
-     , Diffy Double usedParams
+     , DiffInterp 'Point    brushParams
-     , Diffy ( 𝕀 Double ) ( 𝕀 usedParams )
+     , DiffInterp 'Interval brushParams
+     , HasChainRule Double 2 usedParams
+     , HasChainRule ( 𝕀 Double ) 2 ( 𝕀 usedParams )
+     , HasChainRule Double 2 brushParams
+     , HasChainRule ( 𝕀 Double ) 2 ( 𝕀 brushParams )
+     , Traversable ( D 2 brushParams )
+
+  --   , Diffy Double usedParams
+  --   , Diffy ( 𝕀 Double ) ( 𝕀 usedParams )
 
      -- Debugging.
      , Show ptData, Show brushParams
      )
   => ( ptData -> usedParams )
   -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing
-  -> ( forall i. Differentiable i brushParams
+  -> ( forall i
+       .  DiffInterp i brushParams
        => Proxy# i
        -> ( forall a. a -> I i a )
        -> I i brushParams ~> Spline Closed () ( I i ( ℝ 2 ) )
@@ -461,7 +476,7 @@ outlineFunction
 outlineFunction ptParams toBrushParams brushFromParams sp0 crv =
   let
     pathAndUsedParams :: forall i
-                      .  ( D ( I i ( ℝ 1 ) ) ~ D ( ℝ 1 )
+                      .  ( D 2 ( I i ( ℝ 1 ) ) ~ D 2 ( ℝ 1 )
                          , Coercible ( I i ( ℝ 1 ) ) ( I i Double )
                          , Module ( I i Double ) ( T ( I i ( ℝ 2 ) ) )
                          , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) )
@@ -492,7 +507,7 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv =
     curvesI :: 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 'Interval )
     curvesI = brushStrokeData @'Interval @brushParams
                pathI
-               ( usedParamsI `chainRule` linear ( nonDecreasing toBrushParams ) )
+               ( chainRule @( 𝕀 Double ) @2 usedParamsI $ linear ( nonDecreasing toBrushParams ) )
                ( brushFromParams @'Interval proxy# singleton )
 
     usedParamsI :: 𝕀ℝ 1 ~> 𝕀 usedParams
@@ -501,7 +516,7 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv =
 
     fwdBwd :: OutlineFn
     fwdBwd t
-      = solveEnvelopeEquations path_t path'_t ( fwdOffset, bwdOffset ) curves
+      = solveEnvelopeEquations t path_t path'_t ( fwdOffset, bwdOffset ) curves
     --  = ( ( offset fwdOffset • path_t,       path'_t )
     --    , ( offset bwdOffset • path_t, -1 *^ path'_t ) )
       where
@@ -509,26 +524,26 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv =
         curves :: Seq ( ℝ 1 -> StrokeDatum Point )
         curves = brushStrokeData @Point @brushParams
                    path
-                   ( usedParams `chainRule` linear toBrushParams )
+                   ( chainRule @Double usedParams $ linear toBrushParams )
                    ( brushFromParams @Point proxy# id )
                    t
 
         fwdOffset = withTangent         path'_t   brush_t
         bwdOffset = withTangent ( -1 *^ path'_t ) brush_t
 
-        D1 path_t path'_t _ = runD path t
+        D21 path_t path'_t _ = runD path t
-        D1 params_t _ _     = runD usedParams t
+        D21 params_t _ _     = runD usedParams t
-        brush_t             = value @Double @brushParams
+        brush_t              = value @Double @2 @brushParams
-                            $ runD ( brushFromParams @Point proxy# id )
+                             $ runD ( brushFromParams @Point proxy# id )
-                            $ toBrushParams params_t
+                             $ toBrushParams params_t
 
     bisSols = bisection 0.0001 curvesI
 
-  in trace
+  in --trace
-         ( unlines $
+     --    ( unlines $
-           ( "bisectionMethod: #(possible zeroes) = " ++ show ( length bisSols ) ) :
+     --      ( "bisectionMethod: #(possible zeroes) = " ++ show ( length bisSols ) ) :
-           "" :
+     --      "" :
-           map show bisSols )
+     --      map show bisSols )
          fwdBwd
 
 -----------------------------------
@@ -821,11 +836,11 @@ 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 :: Module r ( T v )
-            => D ( ℝ 1 ) v -- ^ stroke path \( p(t) \)
+            => D 2 ( ℝ 1 ) v -- ^ stroke path \( p(t) \)
-            -> D ( ℝ 2 ) v -- ^ brush \( b(t,s) \)
+            -> D 2 ( ℝ 2 ) v -- ^ brush \( b(t,s) \)
-            -> D ( ℝ 2 ) v
+            -> D 2 ( ℝ 2 ) v
-brushStroke ( D1 p dpdt d2pdt2 ) ( D2 b dbdt dbds d2bdt2 d2bdtds d2bds2 ) =
+brushStroke ( D21 p dpdt d2pdt2 ) ( D22 b dbdt dbds d2bdt2 d2bdtds d2bds2 ) =
-  D2 ( unT $ T p ^+^ T b )
+  D22 ( unT $ T p ^+^ T b )
      -- c = p + b
 
    ( dpdt ^+^ dbdt ) dbds
@@ -851,13 +866,13 @@ brushStroke ( D1 p dpdt d2pdt2 ) ( D2 b dbdt dbds d2bdt2 d2bdtds d2bds2 ) =
 --
 -- NB: if \( \frac{\partial E}{\partial s} \) is zero, the total derivative is ill-defined.
 envelopeEquation :: forall i
-                 .  ( D ( I i ( ℝ 2 ) ) ~ D ( ℝ 2 )
+                 .  ( D 2 ( I i ( ℝ 2 ) ) ~ D 2 ( ℝ 2 )
                     , Cross ( I i Double ) ( T ( I i ( ℝ 2 ) ) )
                     , Fractional ( I i Double )
                     )
-                 => D ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) )
+                 => D 2 ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) )
                  -> ( I i Double, T ( I i ( ℝ 2 ) ), T ( I i ( ℝ 2 ) ), I i Double, I i Double )
-envelopeEquation ( D2 _ dcdt dcds d2cdt2 d2cdtds d2cds2 ) =
+envelopeEquation ( D22 _ dcdt dcds d2cdt2 d2cdtds d2cds2 ) =
   let ee = dcdt `cross` dcds
       dEdt = d2cdt2  `cross` dcds + dcdt `cross` d2cdtds
       dEds = d2cdtds `cross` dcds + dcdt `cross` d2cds2
@@ -874,41 +889,41 @@ envelopeEquation ( D2 _ dcdt dcds d2cdt2 d2cdtds d2cds2 ) =
 -- | Linear interpolation, as a differentiable function.
 line :: forall i b
      .  ( Module ( I i Double ) ( T b ), Torsor ( T b ) b
-        , D ( I i ( ℝ 1 ) ) ~ D ( ℝ 1 )
+        , D 2 ( I i ( ℝ 1 ) ) ~ D 2 ( ℝ 1 )
         , Coercible ( I i ( ℝ 1 ) ) ( I i Double )
         )
      => Segment b -> I i ( ℝ 1 ) ~> b
 line ( Segment a b ) = D \ ( coerce -> t ) ->
-  D1 ( lerp @( T b ) t a b )
+  D21 ( lerp @( T b ) t a b )
-     ( a --> b )
+      ( a --> b )
-     origin
+      origin
 
 -- | A quadratic Bézier curve, as a differentiable function.
 bezier2 :: forall i b
         .  ( Module ( I i Double ) ( T b ), Torsor ( T b ) b
-           , D ( I i ( ℝ 1 ) ) ~ D ( ℝ 1 )
+           , D 2 ( I i ( ℝ 1 ) ) ~ D 2 ( ℝ 1 )
            , Coercible ( I i ( ℝ 1 ) ) ( I i Double )
            )
         => Quadratic.Bezier b -> I i ( ℝ 1 ) ~> b
 bezier2 bez = D \ ( coerce -> t ) ->
-  D1 ( Quadratic.bezier @( T b ) bez t )
+  D21 ( Quadratic.bezier @( T b ) bez t )
-     ( Quadratic.bezier'         bez t )
+      ( Quadratic.bezier'         bez t )
-     ( Quadratic.bezier''        bez   )
+      ( Quadratic.bezier''        bez   )
 
 -- | A cubic Bézier curve, as a differentiable function.
 bezier3 :: forall i b
         .  ( Module ( I i Double ) ( T b ), Torsor ( T b ) b
-           , D ( I i ( ℝ 1 ) ) ~ D ( ℝ 1 )
+           , D 2 ( I i ( ℝ 1 ) ) ~ D 2 ( ℝ 1 )
            , Coercible ( I i ( ℝ 1 ) ) ( I i Double )
            )
         => Cubic.Bezier b -> I i ( ℝ 1 ) ~> b
 bezier3 bez = D \ ( coerce -> t ) ->
-  D1 ( Cubic.bezier @( T b ) bez t )
+  D21 ( Cubic.bezier @( T b ) bez t )
-     ( Cubic.bezier'         bez t )
+      ( Cubic.bezier'         bez t )
-     ( Cubic.bezier''        bez t )
+      ( Cubic.bezier''        bez t )
 
 splineCurveFns :: forall i
-               .  ( D ( I i ( ℝ 1 ) ) ~ D ( ℝ 1 )
+               .  ( D 2 ( I i ( ℝ 1 ) ) ~ D 2 ( ℝ 1 )
                   , Module ( I i Double ) ( T ( I i ( ℝ 2 ) ) )
                   , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) )
                   , Coercible ( I i ( ℝ 1 ) ) ( I i Double ) )
@@ -934,12 +949,13 @@ splineCurveFns spls
 
 -- | Solve the envelope equations at a given point \( t = t_0 \), to find
 -- \( s_0 \) such that \( c(t_0, s_0) \) is on the envelope of the brush stroke.
-solveEnvelopeEquations :: ℝ 2
+solveEnvelopeEquations :: ℝ 1 -- ^ @t@ (for debugging only)
+                       -> ℝ 2
                        -> T ( ℝ 2 )
                        -> ( Offset, Offset )
                        -> Seq ( ℝ 1 -> StrokeDatum Point )
                        -> ( ( ℝ 2, T ( ℝ 2 ) ), ( ℝ 2, T ( ℝ 2 ) ) )
-solveEnvelopeEquations path_t path'_t ( fwdOffset, bwdOffset ) strokeData
+solveEnvelopeEquations _t path_t path'_t ( fwdOffset, bwdOffset ) strokeData
   = ( fwdSol, ( bwdPt, -1 *^ bwdTgt ) )
 
   where
@@ -947,6 +963,7 @@ solveEnvelopeEquations path_t path'_t ( fwdOffset, bwdOffset ) strokeData
     -- !_ = trace
     --      ( unlines
     --        [ "solveEnvelopeEquation"
+    --        , "  t: " ++ show _t
     --        , "  pt: " ++ show path_t
     --        , "  tgt: "  ++ show path'_t
     --        , "  fwdOffset: " ++ show fwdOffset
@@ -1016,7 +1033,7 @@ solveEnvelopeEquations path_t path'_t ( fwdOffset, bwdOffset ) strokeData
             in --trace
                --  ( unlines
                --    [ "solveEnvelopeEquations"
-               --    , "      t = " ++ show t
+               --    , "      t = " ++ show _t
                --    , "      s = " ++ show s
                --    , "      c = " ++ show dstroke
                --    , "      E = " ++ show _ee
@@ -1024,7 +1041,8 @@ solveEnvelopeEquations path_t path'_t ( fwdOffset, bwdOffset ) strokeData
                --    , "  ∂E/∂s = " ++ show 𝛿E𝛿s
                --    , "  dc/dt = " ++ show totDeriv
                --    ] )
-               ( good, ℝ1 s, value @Double @( ℝ 2 ) dstroke, totDeriv )
+               ( good, ℝ1 s, value @Double @2 @( ℝ 2 ) dstroke
+                 , totDeriv )
 
     eqn :: ( ℝ 1 -> StrokeDatum Point ) -> ( Double -> ( Double, Double ) )
     eqn f s =
@@ -1046,10 +1064,10 @@ instance Applicative ZipSeq where
   liftA2 f ( ZipSeq xs ) ( ZipSeq ys ) = ZipSeq ( Seq.zipWith f xs ys )
 
 brushStrokeData :: forall i brushParams
-                .  ( Diffy ( I i Double ) ( I i brushParams )
+                .  ( Differentiable i brushParams
                    , Fractional ( I i Double )
-                   , D ( I i ( ℝ 1 ) ) ~ D ( ℝ 1 )
+                   , D 2 ( I i ( ℝ 1 ) ) ~ D 2 ( ℝ 1 )
-                   , D ( I i ( ℝ 2 ) ) ~ D ( ℝ 2 )
+                   , D 2 ( I i ( ℝ 2 ) ) ~ D 2 ( ℝ 2 )
                    , Cross ( I i Double ) ( T ( I i ( ℝ 2 ) ) )
                    , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) )
                    , Coercible ( I i ( ℝ 1 ) ) ( I i Double )
@@ -1065,26 +1083,27 @@ brushStrokeData :: forall i brushParams
 brushStrokeData path params brush =
   \ t ->
     let
-        dpath_t :: D ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) )
+        dpath_t :: D 2 ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) )
         !dpath_t = runD path t
-        dparams_t :: D ( I i ( ℝ 1 ) ) ( I i brushParams )
+        dparams_t :: D 2 ( I i ( ℝ 1 ) ) ( I i brushParams )
-        !dparams_t@( D1 { v = params_t } ) = runD params t
+        !dparams_t@( D21 { _D21_v = params_t } ) = runD params t
-        dbrush_params :: D ( I i brushParams ) ( Spline Closed () ( I i ( ℝ 2 ) ) )
+        dbrush_params :: D 2 ( I i brushParams ) ( Spline Closed () ( I i ( ℝ 2 ) ) )
         !dbrush_params = runD brush params_t
-        splines :: Seq ( D ( I i brushParams ) ( I i ( ℝ 1 ) ~> I i ( ℝ 2 ) ) )
+        splines :: Seq ( D 2 ( I i brushParams ) ( I i ( ℝ 1 ) ~> I i ( ℝ 2 ) ) )
         !splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns @i ) dbrush_params
-        dbrushes_t :: Seq ( I i ( ℝ 1 ) -> D ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) )
+        dbrushes_t :: Seq ( I i ( ℝ 1 ) -> D 2 ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) )
-        !dbrushes_t = force $ fmap ( uncurryD . ( dparams_t `chain` ) ) splines
+        !dbrushes_t = force $ fmap ( uncurryD2 . ( chain @(I i Double) @2 dparams_t ) ) splines
+          -- This is the crucial use of the chain rule.
 
     in fmap ( mkStrokeDatum dpath_t ) dbrushes_t
   where
 
-    mkStrokeDatum :: D ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) )
+    mkStrokeDatum :: D 2 ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) )
-                  -> ( I i ( ℝ 1 ) -> D ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) )
+                  -> ( I i ( ℝ 1 ) -> D 2 ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) )
                   -> ( I i ( ℝ 1 ) -> StrokeDatum i )
     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
+          dstroke@( D22 _c _𝛿c𝛿t _𝛿c𝛿s _ _ _ ) = brushStroke dpath_t dbrush_t_s
           ( ee, dcdt, 𝛿E𝛿sdcdt, 𝛿E𝛿t, 𝛿E𝛿s ) = envelopeEquation @i dstroke
       in  -- trace
           --   ( unlines
@@ -1111,14 +1130,14 @@ brushStrokeData path params brush =
 data StrokeDatum i
   = StrokeDatum
   { -- | Path \( p(t_0) \) (with its 1st and 2nd derivatives).
-    dpath    :: D ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) )
+    dpath    :: D 2 ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) )
     -- | Brush shape \( b(t_0, s_0) \) (with its 1st and 2nd derivatives).
-  , dbrush   :: D ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) )
+  , dbrush   :: D 2 ( 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) \) (with its 1st and 2nd derivatives).
-  , dstroke  :: D ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) )
+  , dstroke  :: D 2 ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) )
     -- | Envelope
     --
     -- \[ E(t_0,s_0) = \left ( \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} \right )_{(t_0,s_0)}. \]
@@ -1209,6 +1228,6 @@ bisection minWidth eqs = bisect initialCands [] []
               && cmpℝ2 (>) ( getRounded ( Interval.sup 𝛿E𝛿sdcdt ) ) ( ℝ2 0 0 )
            return ( ee, 𝛿E𝛿sdcdt )
 
-cmpℝ2 :: (Double -> Double -> Bool) -> ℝ 2 -> ℝ 2 -> Bool
+cmpℝ2 :: ( Double -> Double -> Bool ) -> ℝ 2 -> ℝ 2 -> Bool
 cmpℝ2 cmp ( ℝ2 x1 y1 ) ( ℝ2 x2 y2 )
   = cmp x1 x2 && cmp y1 y2

src/splines/Math/Differentiable.hs

@@ -0,0 +1,66 @@
+{-# LANGUAGE UndecidableInstances #-}
+
+module Math.Differentiable
+  ( Differentiable, DiffInterp )
+  where
+
+-- base
+import Data.Kind
+  ( Type, Constraint )
+import GHC.TypeNats
+  ( Nat )
+
+-- MetaBrush
+import Math.Algebra.Dual
+  ( D, HasChainRule )
+import Math.Linear
+import Math.Module
+import Math.Interval
+  ( Extent(..), type I )
+import Math.Ring
+  ( Transcendental )
+
+--------------------------------------------------------------------------------
+
+type ExtentOrder :: Extent -> Nat
+type family ExtentOrder e where
+  ExtentOrder i = 2
+  --ExtentOrder 'Point    = 2
+  --ExtentOrder 'Interval = 2
+-- Currently we're doing order 2 derivatives for the brush stroke fitting,
+-- but order 3 derivatives for the interval Newton method to find cusps.
+-- TODO: using 2 for both until migration finishes.
+
+type Differentiable :: Extent -> Type -> Constraint
+class
+  ( Module ( I i Double ) ( T ( I i u ) )
+  , HasChainRule ( I i Double ) ( ExtentOrder i ) ( I i u )
+  , Traversable ( D ( ExtentOrder i ) ( I i u ) )
+  ) => Differentiable i u
+instance
+  ( Module ( I i Double ) ( T ( I i u ) )
+  , HasChainRule ( I i Double ) ( ExtentOrder i ) ( I i u )
+  , Traversable ( D ( ExtentOrder i ) ( I i u ) )
+  ) => Differentiable i u
+
+type DiffInterp :: Extent -> Type -> Constraint
+class
+  ( Differentiable 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 ) )
+           ( D ( ExtentOrder i ) ( I i u ) ( I i ( ℝ 2 ) ) )
+  , Transcendental ( D ( ExtentOrder i ) ( I i u ) ( I i Double ) )
+  , Applicative ( D ( ExtentOrder i ) ( I i u ) )
+  , Representable ( I i Double ) ( I i u )
+  ) => DiffInterp i u
+instance
+  ( Differentiable 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 ) )
+           ( D ( ExtentOrder i ) ( I i u ) ( I i ( ℝ 2 ) ) )
+  , Transcendental ( D ( ExtentOrder i ) ( I i u ) ( I i Double ) )
+  , Applicative ( D ( ExtentOrder i ) ( I i u ) )
+  , Representable ( I i Double ) ( I i u )
+  ) => DiffInterp i u

src/splines/Math/Interval.hs

@@ -0,0 +1,441 @@
+{-# LANGUAGE ScopedTypeVariables  #-}
+{-# LANGUAGE TemplateHaskell      #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+{-# OPTIONS_GHC -Wno-orphans #-}
+
+module Math.Interval
+  ( 𝕀, 𝕀ℝ
+  , Extent(..), type I
+  , singleton, nonDecreasing
+  )
+  where
+
+-- base
+import Prelude hiding ( Num(..) )
+import Data.Coerce
+  ( coerce )
+import Data.Kind
+  ( Type )
+import Data.Monoid
+  ( Sum(..) )
+
+-- acts
+import Data.Act
+  ( Act((•)), Torsor((-->)) )
+
+-- groups
+import Data.Group
+  ( Group(..) )
+
+-- groups-generic
+import Data.Group.Generics
+  ( )
+
+-- rounded-hw
+import Numeric.Rounded.Hardware
+  ( Rounded(..) )
+import Numeric.Rounded.Hardware.Interval.NonEmpty
+  ( Interval(..) )
+
+-- splines
+import Math.Algebra.Dual
+import Math.Algebra.Dual.Internal
+  ( chainRuleQ )
+import Math.Interval.Internal
+  ( type 𝕀 )
+import Math.Linear
+  ( ℝ(..), T(..)
+  , RepresentableQ(..)
+  )
+import Math.Module
+import Math.Monomial
+import Math.Ring
+
+--------------------------------------------------------------------------------
+-- Interval arithmetic using rounded-hw library.
+
+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 = I ( Rounded a ) ( Rounded a )
+
+-- | Turn a non-decreasing function into a function on intervals.
+nonDecreasing :: ( a -> b ) -> 𝕀 a -> 𝕀 b
+nonDecreasing f ( I ( Rounded lo ) ( Rounded hi ) ) =
+  I ( Rounded $ f lo ) ( Rounded $ f hi )
+
+
+deriving via ViaAbelianGroup ( T ( 𝕀 Double ) )
+  instance Semigroup ( T ( 𝕀 Double ) )
+deriving via ViaAbelianGroup ( T ( 𝕀 Double ) )
+  instance Monoid ( T ( 𝕀 Double ) )
+deriving via ViaAbelianGroup ( T ( 𝕀 Double ) )
+  instance Group ( T ( 𝕀 Double ) )
+
+instance Act   ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where
+  T g • a = coerce ( Sum g • a )
+instance Torsor ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where
+  a --> b = T $ getSum ( a --> b )
+
+-------------------------------------------------------------------------------
+
+instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ 0 ) ) where
+  origin      = T $$( tabulateQ \ _ -> [|| unT $ origin ||] )
+  T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] )
+  T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] )
+  k *^ T a    = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] )
+
+instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ 1 ) ) where
+  origin      = T $$( tabulateQ \ _ -> [|| unT $ origin ||] )
+  T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] )
+  T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] )
+  k *^ T a    = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] )
+
+instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ 2 ) ) where
+  origin      = T $$( tabulateQ \ _ -> [|| unT $ origin ||] )
+  T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] )
+  T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] )
+  k *^ T a    = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] )
+
+instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ 3 ) ) where
+  origin      = T $$( tabulateQ \ _ -> [|| unT $ origin ||] )
+  T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] )
+  T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] )
+  k *^ T a    = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] )
+
+instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ 4 ) ) where
+  origin      = T $$( tabulateQ \ _ -> [|| unT $ origin ||] )
+  T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] )
+  T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] )
+  k *^ T a    = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] )
+
+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 !x1x2 = I ( Rounded x1_lo ) ( Rounded x1_hi ) * I ( Rounded x2_lo ) ( Rounded x2_hi )
+            !y1y2 = I ( Rounded y1_lo ) ( Rounded y1_hi ) * I ( Rounded y2_lo ) ( Rounded y2_hi )
+        in x1x2 + y1y2
+
+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 !x1y2 = I ( Rounded x1_lo ) ( Rounded x1_hi ) * I ( Rounded y2_lo ) ( Rounded y2_hi )
+            !y2x1 = I ( Rounded x2_lo ) ( Rounded x2_hi ) * I ( Rounded y1_lo ) ( Rounded y1_hi )
+        in x1y2 - y2x1
+
+deriving via ViaModule ( 𝕀 Double ) ( T ( 𝕀ℝ n ) )
+  instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) => Semigroup ( T ( 𝕀ℝ n ) )
+deriving via ViaModule ( 𝕀 Double ) ( T ( 𝕀ℝ n ) )
+  instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) => Monoid ( T ( 𝕀ℝ n ) )
+deriving via ViaModule ( 𝕀 Double ) ( T ( 𝕀ℝ n ) )
+  instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) => Group ( T ( 𝕀ℝ n ) )
+
+deriving via ViaModule ( 𝕀 Double ) ( 𝕀ℝ n )
+  instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) => Act ( T ( 𝕀ℝ n ) ) ( 𝕀ℝ n )
+deriving via ( ViaModule ( 𝕀 Double ) ( 𝕀ℝ n ) )
+  instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) => Torsor ( T ( 𝕀ℝ n ) ) ( 𝕀ℝ n )
+
+--------------------------------------------------------------------------------
+-- HasChainRule instances.
+
+instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 0 ) where
+  konst w = D0 w
+  linearD f v = D0 ( f v )
+  value ( D0 v ) = v
+  chain _ ( D0 gfx ) = D21 gfx origin origin
+
+instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 0 ) where
+  konst w = D0 w
+  linearD f v = D0 ( f v )
+  value ( D0 v ) = v
+  chain _ ( D0 gfx ) = D31 gfx origin origin origin
+
+instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 1 ) where
+
+  konst   :: forall w. AbelianGroup w => w -> D2𝔸1 w
+  konst w =
+    let !o = fromInteger @w 0
+    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+
+  value df = $$( monIndex [|| df ||] zeroMonomial )
+
+  linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 1 -> w ) -> 𝕀ℝ 1 -> D2𝔸1 w
+  linearD f v =
+    let !o = origin @( 𝕀 Double ) @( T w )
+    in $$( monTabulate \ mon ->
+              if | isZeroMonomial mon
+                 -> [|| f v ||]
+                 | Just i <- isLinear mon
+                 -> [|| f $$( tabulateQ \ j ->
+                                if | j == i
+                                   -> [|| 1 ||]
+                                   | otherwise
+                                   -> [|| 0 ||]
+                            ) ||]
+                 | otherwise
+                 -> [|| unT o ||]
+         )
+
+  chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D2𝔸1 ( 𝕀ℝ 1 ) -> D2𝔸1 w -> D2𝔸1 w
+  chain !df !dg =
+    let !o = origin @( 𝕀 Double ) @( T w )
+        !p = (^+^)  @( 𝕀 Double ) @( T w )
+        !s = (^*)   @( 𝕀 Double ) @( T w )
+    in $$( chainRuleQ
+         [|| o ||] [|| p ||] [|| s ||]
+         [|| df ||] [|| dg ||] )
+
+instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 1 ) where
+
+  konst   :: forall w. AbelianGroup w => w -> D3𝔸1 w
+  konst w =
+    let !o = fromInteger @w 0
+    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+
+  value df = $$( monIndex [|| df ||] zeroMonomial )
+
+  linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 1 -> w ) -> 𝕀ℝ 1 -> D3𝔸1 w
+  linearD f v =
+    let !o = origin @( 𝕀 Double ) @( T w )
+    in $$( monTabulate \ mon ->
+              if | isZeroMonomial mon
+                 -> [|| f v ||]
+                 | Just i <- isLinear mon
+                 -> [|| f $$( tabulateQ \ j ->
+                                if | j == i
+                                   -> [|| 1 ||]
+                                   | otherwise
+                                   -> [|| 0 ||]
+                            ) ||]
+                 | otherwise
+                 -> [|| unT o ||]
+         )
+
+  chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D3𝔸1 ( 𝕀ℝ 1 ) -> D3𝔸1 w -> D3𝔸1 w
+  chain !df !dg =
+    let !o = origin @( 𝕀 Double ) @( T w )
+        !p = (^+^)  @( 𝕀 Double ) @( T w )
+        !s = (^*)   @( 𝕀 Double ) @( T w )
+    in $$( chainRuleQ
+         [|| o ||] [|| p ||] [|| s ||]
+         [|| df ||] [|| dg ||] )
+
+instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 2 ) where
+
+  konst   :: forall w. AbelianGroup w => w -> D2𝔸2 w
+  konst w =
+    let !o = fromInteger @w 0
+    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+
+  value df = $$( monIndex [|| df ||] zeroMonomial )
+
+  linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 2 -> w ) -> 𝕀ℝ 2 -> D2𝔸2 w
+  linearD f v =
+    let !o = origin @( 𝕀 Double ) @( T w )
+    in $$( monTabulate \ mon ->
+              if | isZeroMonomial mon
+                 -> [|| f v ||]
+                 | Just i <- isLinear mon
+                 -> [|| f $$( tabulateQ \ j ->
+                                if | j == i
+                                   -> [|| 1 ||]
+                                   | otherwise
+                                   -> [|| 0 ||]
+                            ) ||]
+                 | otherwise
+                 -> [|| unT o ||]
+         )
+
+  chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D2𝔸1 ( 𝕀ℝ 2 ) -> D2𝔸2 w -> D2𝔸1 w
+  chain !df !dg =
+    let !o = origin @( 𝕀 Double ) @( T w )
+        !p = (^+^)  @( 𝕀 Double ) @( T w )
+        !s = (^*)   @( 𝕀 Double ) @( T w )
+    in $$( chainRuleQ
+         [|| o ||] [|| p ||] [|| s ||]
+         [|| df ||] [|| dg ||] )
+
+instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 2 ) where
+
+  konst   :: forall w. AbelianGroup w => w -> D3𝔸2 w
+  konst w =
+    let !o = fromInteger @w 0
+    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+
+  value df = $$( monIndex [|| df ||] zeroMonomial )
+
+  linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 2 -> w ) -> 𝕀ℝ 2 -> D3𝔸2 w
+  linearD f v =
+    let !o = origin @( 𝕀 Double ) @( T w )
+    in $$( monTabulate \ mon ->
+              if | isZeroMonomial mon
+                 -> [|| f v ||]
+                 | Just i <- isLinear mon
+                 -> [|| f $$( tabulateQ \ j ->
+                                if | j == i
+                                   -> [|| 1 ||]
+                                   | otherwise
+                                   -> [|| 0 ||]
+                            ) ||]
+                 | otherwise
+                 -> [|| unT o ||]
+         )
+
+  chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D3𝔸1 ( 𝕀ℝ 2 ) -> D3𝔸2 w -> D3𝔸1 w
+  chain !df !dg =
+    let !o = origin @( 𝕀 Double ) @( T w )
+        !p = (^+^)  @( 𝕀 Double ) @( T w )
+        !s = (^*)   @( 𝕀 Double ) @( T w )
+    in $$( chainRuleQ
+         [|| o ||] [|| p ||] [|| s ||]
+         [|| df ||] [|| dg ||] )
+
+instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 3 ) where
+
+  konst   :: forall w. AbelianGroup w => w -> D2𝔸3 w
+  konst w =
+    let !o = fromInteger @w 0
+    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+
+  value df = $$( monIndex [|| df ||] zeroMonomial )
+
+  linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 3 -> w ) -> 𝕀ℝ 3 -> D2𝔸3 w
+  linearD f v =
+    let !o = origin @( 𝕀 Double ) @( T w )
+    in $$( monTabulate \ mon ->
+              if | isZeroMonomial mon
+                 -> [|| f v ||]
+                 | Just i <- isLinear mon
+                 -> [|| f $$( tabulateQ \ j ->
+                                if | j == i
+                                   -> [|| 1 ||]
+                                   | otherwise
+                                   -> [|| 0 ||]
+                            ) ||]
+                 | otherwise
+                 -> [|| unT o ||]
+         )
+
+  chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D2𝔸1 ( 𝕀ℝ 3 ) -> D2𝔸3 w -> D2𝔸1 w
+  chain !df !dg =
+    let !o = origin @( 𝕀 Double ) @( T w )
+        !p = (^+^)  @( 𝕀 Double ) @( T w )
+        !s = (^*)   @( 𝕀 Double ) @( T w )
+    in $$( chainRuleQ
+         [|| o ||] [|| p ||] [|| s ||]
+         [|| df ||] [|| dg ||] )
+
+instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 3 ) where
+
+  konst   :: forall w. AbelianGroup w => w -> D3𝔸3 w
+  konst w =
+    let !o = fromInteger @w 0
+    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+
+  value df = $$( monIndex [|| df ||] zeroMonomial )
+
+  linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 3 -> w ) -> 𝕀ℝ 3 -> D3𝔸3 w
+  linearD f v =
+    let !o = origin @( 𝕀 Double ) @( T w )
+    in $$( monTabulate \ mon ->
+              if | isZeroMonomial mon
+                 -> [|| f v ||]
+                 | Just i <- isLinear mon
+                 -> [|| f $$( tabulateQ \ j ->
+                                if | j == i
+                                   -> [|| 1 ||]
+                                   | otherwise
+                                   -> [|| 0 ||]
+                            ) ||]
+                 | otherwise
+                 -> [|| unT o ||]
+         )
+
+  chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D3𝔸1 ( 𝕀ℝ 3 ) -> D3𝔸3 w -> D3𝔸1 w
+  chain !df !dg =
+    let !o = origin @( 𝕀 Double ) @( T w )
+        !p = (^+^)  @( 𝕀 Double ) @( T w )
+        !s = (^*)   @( 𝕀 Double ) @( T w )
+    in $$( chainRuleQ
+         [|| o ||] [|| p ||] [|| s ||]
+         [|| df ||] [|| dg ||] )
+
+instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 4 ) where
+
+  konst   :: forall w. AbelianGroup w => w -> D2𝔸4 w
+  konst w =
+    let !o = fromInteger @w 0
+    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+
+  value df = $$( monIndex [|| df ||] zeroMonomial )
+
+  linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 4 -> w ) -> 𝕀ℝ 4 -> D2𝔸4 w
+  linearD f v =
+    let !o = origin @( 𝕀 Double ) @( T w )
+    in $$( monTabulate \ mon ->
+              if | isZeroMonomial mon
+                 -> [|| f v ||]
+                 | Just i <- isLinear mon
+                 -> [|| f $$( tabulateQ \ j ->
+                                if | j == i
+                                   -> [|| 1 ||]
+                                   | otherwise
+                                   -> [|| 0 ||]
+                            ) ||]
+                 | otherwise
+                 -> [|| unT o ||]
+         )
+
+  chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D2𝔸1 ( 𝕀ℝ 4 ) -> D2𝔸4 w -> D2𝔸1 w
+  chain !df !dg =
+    let !o = origin @( 𝕀 Double ) @( T w )
+        !p = (^+^)  @( 𝕀 Double ) @( T w )
+        !s = (^*)   @( 𝕀 Double ) @( T w )
+    in $$( chainRuleQ
+         [|| o ||] [|| p ||] [|| s ||]
+         [|| df ||] [|| dg ||] )
+
+instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 4 ) where
+
+  konst   :: forall w. AbelianGroup w => w -> D3𝔸4 w
+  konst w =
+    let !o = fromInteger @w 0
+    in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
+
+  value df = $$( monIndex [|| df ||] zeroMonomial )
+
+  linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 4 -> w ) -> 𝕀ℝ 4 -> D3𝔸4 w
+  linearD f v =
+    let !o = origin @( 𝕀 Double ) @( T w )
+    in $$( monTabulate \ mon ->
+              if | isZeroMonomial mon
+                 -> [|| f v ||]
+                 | Just i <- isLinear mon
+                 -> [|| f $$( tabulateQ \ j ->
+                                if | j == i
+                                   -> [|| 1 ||]
+                                   | otherwise
+                                   -> [|| 0 ||]
+                            ) ||]
+                 | otherwise
+                 -> [|| unT o ||]
+         )
+
+  chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D3𝔸1 ( 𝕀ℝ 4 ) -> D3𝔸4 w -> D3𝔸1 w
+  chain !df !dg =
+    let !o = origin @( 𝕀 Double ) @( T w )
+        !p = (^+^)  @( 𝕀 Double ) @( T w )
+        !s = (^*)   @( 𝕀 Double ) @( T w )
+    in $$( chainRuleQ
+         [|| o ||] [|| p ||] [|| s ||]
+         [|| df ||] [|| dg ||] )

src/splines/Math/Interval/Internal.hs

@@ -0,0 +1,82 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TemplateHaskell      #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+{-# OPTIONS_GHC -Wno-orphans #-}
+
+module Math.Interval.Internal
+  ( type 𝕀 )
+  where
+
+-- base
+import Data.Monoid
+  ( Sum(..) )
+
+-- rounded-hw
+import Numeric.Rounded.Hardware
+  ( Rounded(..) )
+import Numeric.Rounded.Hardware.Interval.NonEmpty
+  ( Interval(..) )
+import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval
+  ( sup, inf, powInt )
+
+-- MetaBrush
+import Math.Linear
+  ( T(..)
+  , RepDim, RepresentableQ(..), Representable(..)
+  )
+import Math.Module
+  ( Module(..) )
+import Math.Ring
+
+--------------------------------------------------------------------------------
+-- Intervals.
+
+type 𝕀 = Interval
+
+deriving via ViaPrelude ( 𝕀 Double )
+  instance AbelianGroup  ( 𝕀 Double )
+deriving via ViaPrelude ( 𝕀 Double )
+  instance AbelianGroup ( T ( 𝕀 Double ) )
+
+instance Ring ( 𝕀 Double ) where
+  (*) = (Prelude.*)
+  x ^ n = Interval.powInt x ( Prelude.fromIntegral n )
+    -- This is very important, as x^2 is not the same as x * x
+    -- in interval arithmetic. This ensures we don't
+    -- accidentally use (^) from Prelude.
+
+deriving via ViaPrelude ( 𝕀 Double )
+  instance Field ( 𝕀 Double )
+
+deriving via ViaPrelude ( 𝕀 Double )
+  instance Transcendental ( 𝕀 Double )
+
+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 $$( f i ) ||] )
+        !hi = tabulateQ @r @u ( \ i -> [|| getRounded $ Interval.sup $$( f i ) ||] )
+    in [|| I ( Rounded $$lo ) ( Rounded $$hi ) ||]
+
+  indexQ d i =
+    [|| case $$d of
+          I ( Rounded lo ) ( Rounded hi ) ->
+            let !lo_i = $$( indexQ @r @u [|| lo ||] i )
+                !hi_i = $$( indexQ @r @u [|| hi ||] i )
+            in I ( Rounded lo_i ) ( Rounded hi_i )
+    ||]
+instance Representable r u => Representable ( 𝕀 r ) ( 𝕀 u ) where
+  tabulate f =
+    let !lo = tabulate @r @u ( \ i -> getRounded $ Interval.inf ( f i ) )
+        !hi = tabulate @r @u ( \ i -> getRounded $ Interval.sup ( f i ) )
+    in I ( Rounded lo ) ( Rounded hi )
+
+  index d i =
+    case d of
+      I ( Rounded lo ) ( Rounded hi ) ->
+        let !lo_i = index @r @u lo i
+            !hi_i = index @r @u hi i
+        in I ( Rounded lo_i ) ( Rounded hi_i )
+
+deriving via Sum ( 𝕀 Double ) instance Module ( 𝕀 Double ) ( T ( 𝕀 Double ) )

src/splines/Math/Linear.hs

@@ -1,8 +1,8 @@
 {-# LANGUAGE AllowAmbiguousTypes  #-}
 {-# LANGUAGE PolyKinds            #-}
 {-# LANGUAGE ScopedTypeVariables  #-}
+{-# LANGUAGE TemplateHaskell      #-}
 {-# LANGUAGE UndecidableInstances #-}
-{-# LANGUAGE UnliftedNewtypes     #-}
 
 module Math.Linear
   ( -- * Points and vectors
@@ -11,22 +11,17 @@ module Math.Linear
   -- * Points and vectors (second version)
   , ℝ(..), T(.., V2, V3)
   , Fin(..), MFin(..)
-  , Dim, Representable(..), ApRep(..)
+  , RepDim, RepresentableQ(..)
-  , injection, projection
+  , Representable(..), injection, projection
   , Vec(..), (!), find, zipIndices
 
-  -- * Intervals
-  , 𝕀, 𝕀ℝ, singleton, nonDecreasing
-  , Extent(..), I
   ) where
 
 -- base
 import Data.Coerce
-  ( Coercible, coerce )
+  ( coerce )
 import Data.Kind
-  ( Type, Constraint )
+  ( Type )
-import Data.Monoid
-  ( Sum(..) )
 import GHC.Generics
   ( Generic, Generic1, Generically(..), Generically1(..) )
 import GHC.TypeNats
@@ -48,18 +43,12 @@ import Data.Group
 import Data.Group.Generics
   ( )
 
--- rounded-hw
+-- MetaBrush
-import Numeric.Rounded.Hardware
+import Math.Linear.Internal
-  ( Rounded(..) )
-import Numeric.Rounded.Hardware.Interval.NonEmpty
-  ( Interval(..) )
-import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval
-  ( sup, inf )
 
 --------------------------------------------------------------------------------
 
 data Mat22 = Mat22 !Double !Double !Double !Double
-
 data Segment p =
   Segment
     { segmentStart :: !p
@@ -72,45 +61,11 @@ data Segment p =
     via Generically1 Segment
   deriving anyclass ( NFData, NFData1 )
 
-instance Show p => Show (Segment p) where
+instance Show p => Show ( Segment p ) where
-  show (Segment s e) = show s ++ " -> " ++ show e
+  show ( Segment s e ) = show s ++ " -> " ++ show e
 
 --------------------------------------------------------------------------------
 
--- | Euclidean space \( \mathbb{R}^n \).
-type ℝ :: Nat -> Type
-data family ℝ n
-data instance ℝ 0 = ℝ0
-  deriving stock ( Show, Eq, Ord, Generic )
-  deriving anyclass NFData
-newtype instance ℝ 1 = ℝ1 { unℝ1 :: Double }
-  deriving stock   ( Generic )
-  deriving newtype ( Show, Eq, Ord, NFData )
-data    instance ℝ 2 = ℝ2 {-# UNPACK #-} !Double {-# UNPACK #-} !Double
-  deriving stock    Generic
-  deriving anyclass NFData
-  deriving stock    ( Show, Eq, Ord )
-data    instance ℝ 3 = ℝ3 {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-# UNPACK #-} !Double
-  deriving stock    Generic
-  deriving anyclass NFData
-  deriving stock    ( Show, Eq, Ord )
-
-data    instance ℝ 4 = ℝ4 {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-# UNPACK #-} !Double
-  deriving stock    Generic
-  deriving anyclass NFData
-  deriving stock    ( Show, Eq, Ord )
-
-deriving via ApRep ( Sum Double ) ( ℝ n )
-  instance Representable Double ( ℝ n ) => Semigroup ( T ( ℝ n ) )
-deriving via ApRep ( Sum Double ) ( ℝ n )
-  instance Representable Double ( ℝ n ) => Monoid    ( T ( ℝ n ) )
-deriving via ApRep ( Sum Double ) ( ℝ n )
-  instance Representable Double ( ℝ n ) => Group     ( T ( ℝ n ) )
-deriving via ApRep Double ( ℝ n )
-  instance Representable Double ( ℝ n ) => Act    ( T ( ℝ n ) ) ( ℝ n )
-deriving via ApRep Double ( ℝ n )
-  instance Representable Double ( ℝ n ) => Torsor ( T ( ℝ n ) ) ( ℝ n )
-
 -- | Tangent space to Euclidean space.
 type T :: Type -> Type
 newtype T e = T { unT :: e }
@@ -130,9 +85,8 @@ instance Act ( T Double ) Double where
 instance Torsor ( T Double ) Double where
   a --> b = T ( b - a )
 
-
+--instance {-# OVERLAPPING #-} Show ( ℝ n ) => Show ( T ( ℝ n ) ) where
-instance {-# OVERLAPPING #-} Show ( ℝ n ) => Show ( T ( ℝ n ) ) where
+--  show ( T p ) = "V" ++ drop 1 ( show p )
-  show ( T p ) = "V" ++ drop 1 ( show p )
 deriving stock instance Show v => Show ( T v )
 
 instance Applicative T where
@@ -151,85 +105,15 @@ pattern V3 x y z = T ( ℝ3 x y z )
 
 --------------------------------------------------------------------------------
 
--- | 1, ..., n
-type Fin :: Nat -> Type
-newtype Fin n = Fin Word
-  deriving stock Eq
-
--- | 0, ..., n
-type MFin :: Nat -> Type
-newtype MFin n = MFin Word
-
-type Dim :: k -> Nat
-type family Dim v
-
-type instance Dim ( ℝ n ) = n
-
-type Representable :: Type -> Type -> Constraint
-class Representable r v | v -> r where
-  tabulate :: ( Fin ( Dim v ) -> r ) -> v
-  index    :: v -> Fin ( Dim v ) -> r
-
-instance Representable Double ( ℝ 0 ) where
-  {-# INLINE tabulate #-}
-  tabulate _ = ℝ0
-  {-# INLINE index #-}
-  index _ _ = 0
-
-instance Representable Double ( ℝ 1 ) where
-  {-# INLINE tabulate #-}
-  tabulate f = ℝ1 ( f ( Fin 1 ) )
-  {-# INLINE index #-}
-  index ( ℝ1 x ) _ = x
-
-instance Representable Double ( ℝ 2 ) where
-  {-# INLINE tabulate #-}
-  tabulate f = ℝ2 ( f ( Fin 1 ) ) ( f ( Fin 2 ) )
-  {-# INLINE index #-}
-  index ( ℝ2 x y ) = \ case
-    Fin 1 -> x
-    _     -> y
-
-instance Representable Double ( ℝ 3 ) where
-  {-# INLINE tabulate #-}
-  tabulate f = ℝ3 ( f ( Fin 1 ) ) ( f ( Fin 2 ) ) ( f ( Fin 3 ) )
-  {-# INLINE index #-}
-  index ( ℝ3 x y z ) = \ case
-    Fin 1 -> x
-    Fin 2 -> y
-    _     -> z
-
-instance Representable Double ( ℝ 4 ) where
-  {-# INLINE tabulate #-}
-  tabulate f = ℝ4 ( f ( Fin 1 ) ) ( f ( Fin 2 ) ) ( f ( Fin 3 ) ) ( f ( Fin 4 ) )
-  {-# INLINE index #-}
-  index ( ℝ4 x y z w ) = \ case
-    Fin 1 -> x
-    Fin 2 -> y
-    Fin 3 -> z
-    _     -> w
-
-{-# INLINE projection #-}
-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 r u, Representable r v )
-          => ( Fin ( Dim v ) -> MFin ( Dim u ) )
-          -> u -> v -> v
-injection f = \ u v ->
-  tabulate \ i -> case f i of
-    MFin 0 -> index v i
-    MFin j -> index u ( Fin j )
-
 infixr 5 `VS`
 type Vec :: Nat -> Type -> Type
 data Vec n a where
   VZ :: Vec 0 a
   VS :: a -> Vec n a -> Vec ( 1 + n ) a
+    -- can't be strict, otherwise we can't conveniently
+    -- unsafeCoerce from lists
+
+deriving stock instance Show a => Show ( Vec n a )
 
 deriving stock instance Functor     ( Vec n )
 deriving stock instance Foldable    ( Vec n )
@@ -241,106 +125,20 @@ VS a _ ! Fin 1 = a
 VS _ a ! Fin i = a ! Fin ( i - 1 )
 _      ! _     = error "impossible: Fin 0 is uninhabited"
 
-find :: forall l a. ( a -> a -> Bool ) -> Vec l a -> a -> MFin l
+find :: forall l a. ( a -> Bool ) -> Vec l a -> MFin l
-find eq v b = MFin ( go 1 v )
+find f v = MFin ( find_ 1 v )
   where
-    go :: Word -> Vec n a -> Word
+    find_ :: Word -> Vec n a -> Word
-    go j ( VS a as )
+    find_ j ( VS a as )
-      | a `eq` b
+      | f a
       = j
       | otherwise
-      = go ( j + 1 ) as
+      = find_ ( j + 1 ) as
-    go _ VZ = 0
+    find_ _ VZ = 0
 
 zipIndices :: forall n a. Vec n a -> [ ( Fin n, a ) ]
-zipIndices = go 0
+zipIndices = zipIndices_ 1
   where
-    go :: forall i. Word -> Vec i a -> [ ( Fin n, a ) ]
+    zipIndices_ :: forall i. Word -> Vec i a -> [ ( Fin n, a ) ]
-    go _ VZ = []
+    zipIndices_ _ VZ = []
-    go i (a `VS` as) = ( Fin i, a ) : go ( i + 1 ) as
+    zipIndices_ i (a `VS` as) = ( Fin i, a ) : zipIndices_ ( i + 1 ) as
-
---------------------------------------------------------------------------------
--- Instances in terms of representable.
-
--- | A newtype to hang off instances for representable functors.
-newtype ApRep r u = ApRep { unApRep :: u }
-
-instance ( Representable r u, Coercible r m, Semigroup m ) => Semigroup ( ApRep m u ) where
-  ApRep a <> ApRep b = ApRep $ tabulate @r @u \ i ->
-    coerce $ (<>) @m ( coerce ( index @r @u a i ) ) ( coerce ( index @r @u b i ) )
-  {-# INLINE (<>) #-}
-instance ( Representable r u, Coercible r m, Monoid m ) => Monoid ( ApRep m u ) where
-  mempty = ApRep $ tabulate @r @u \ _ -> coerce $ mempty @m
-  {-# INLINE mempty #-}
-instance ( Representable r u, Coercible r m, Group m ) => Group ( ApRep m u ) where
-  invert ( ApRep a ) = ApRep $ tabulate @r @u \ i ->
-    coerce $ invert @m $ coerce ( index @r @u a i )
-  {-# INLINE invert #-}
-instance ( Act ( T r ) r , Semigroup ( T u ), Representable r u ) => Act ( T u ) ( ApRep r u ) where
-  T g • ApRep a = ApRep $ tabulate @r @u \ i ->
-    coerce $ (•) @(T r) @r ( coerce $ index @r @u g i ) ( coerce ( index @r @u a i ) )
-  {-# INLINE (•) #-}
-instance ( Torsor ( T r ) r , Group ( T u ), Representable r u ) => Torsor ( T u ) ( ApRep r u ) where
-  ApRep a --> ApRep b = T $ tabulate @r @u \ i ->
-    coerce $ (-->) @(T r) @r ( coerce $ index @r @u a i ) ( coerce ( index @r @u b i ) )
-  {-# INLINE (-->) #-}
-
---------------------------------------------------------------------------------
--- Intervals.
-
-type 𝕀 = Interval
-type 𝕀ℝ n = 𝕀 ( ℝ n )
-
--- 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 = I ( Rounded a ) ( Rounded a )
-
--- | Turn a non-decreasing function into a function on intervals.
-nonDecreasing :: ( a -> b ) -> 𝕀 a -> 𝕀 b
-nonDecreasing f ( I ( Rounded lo ) ( Rounded hi ) ) =
-  I ( Rounded $ f lo ) ( Rounded $ f hi )
-
-type instance Dim ( 𝕀 u ) = Dim u
-
-instance Representable r u => Representable ( 𝕀 r ) ( 𝕀 u ) where
-  tabulate f =
-    let !lo = tabulate @r @u ( \ i -> getRounded $ Interval.inf ( f i ) )
-        !hi = tabulate @r @u ( \ i -> getRounded $ Interval.sup ( f i ) )
-    in I ( Rounded lo ) ( Rounded hi )
-  {-# INLINE tabulate #-}
-  index ( I ( Rounded lo ) ( Rounded hi ) ) i =
-    let !lo_i = index @r @u lo i
-        !hi_i = index @r @u hi i
-    in I ( Rounded lo_i ) ( Rounded hi_i )
-  {-# INLINE index #-}
-
-deriving via ApRep ( Sum ( 𝕀 Double ) ) ( 𝕀ℝ n )
-  instance Representable ( 𝕀 Double ) ( 𝕀ℝ n ) => Semigroup ( T ( 𝕀ℝ n ) )
-deriving via ApRep ( Sum ( 𝕀 Double ) ) ( 𝕀ℝ n )
-  instance Representable ( 𝕀 Double ) ( 𝕀ℝ n ) => Monoid    ( T ( 𝕀ℝ n ) )
-deriving via ApRep ( Sum ( 𝕀 Double ) ) ( 𝕀ℝ n )
-  instance Representable ( 𝕀 Double ) ( 𝕀ℝ n ) => Group     ( T ( 𝕀ℝ n ) )
-
-deriving via Sum ( 𝕀 Double )
-  instance Semigroup ( T ( 𝕀 Double ) )
-deriving via Sum ( 𝕀 Double )
-  instance Monoid ( T ( 𝕀 Double ) )
-deriving via Sum ( 𝕀 Double )
-  instance Group ( T ( 𝕀 Double ) )
-
-instance Act ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where
-  T g • a = coerce ( Sum g • a )
-instance Torsor ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where
-  a --> b = T $ getSum ( a --> b )
-
-deriving via ApRep ( 𝕀 Double ) ( 𝕀ℝ n )
-  instance Representable ( 𝕀 Double ) ( 𝕀ℝ n ) => Act    ( T ( 𝕀ℝ n ) ) ( 𝕀ℝ n )
-deriving via ApRep ( 𝕀 Double ) ( 𝕀ℝ n )
-  instance Representable ( 𝕀 Double ) ( 𝕀ℝ n ) => Torsor ( T ( 𝕀ℝ n ) ) ( 𝕀ℝ n )

src/splines/Math/Linear/Dual.hs

@@ -1,863 +0,0 @@
-{-# LANGUAGE AllowAmbiguousTypes   #-}
-{-# LANGUAGE DuplicateRecordFields #-}
-{-# LANGUAGE QuantifiedConstraints #-}
-{-# LANGUAGE ScopedTypeVariables   #-}
-{-# LANGUAGE UndecidableInstances  #-}
-
-module Math.Linear.Dual where
-
--- base
-import Control.Applicative
-  ( liftA2 )
-import Data.Coerce
-  ( Coercible, coerce )
-import Data.Foldable
-  ( toList )
-import Data.Kind
-  ( Type, Constraint )
-import GHC.Exts
-  ( Any, Word#, proxy#, plusWord# )
-import GHC.Generics
-  ( Generic, Generic1(..), Generically1(..) )
-import GHC.TypeNats
-import Unsafe.Coerce
-  ( unsafeCoerce )
-
--- acts
-import Data.Act
-  ( Torsor )
-
--- rounded-hw
-import Numeric.Rounded.Hardware
-  ( Rounded(..) )
-import Numeric.Rounded.Hardware.Interval.NonEmpty
-  ( Interval(..) )
-
--- MetaBrush
-import Math.Module
-  ( Module(..) )
-import Math.Linear
-
---------------------------------------------------------------------------------
-
--- | @C n u v@ is the space of @C^k@-differentiable maps from @u@ to @v@.
-type C :: Nat -> Type -> Type -> Type
-newtype C k u v = D { runD :: u -> D k u v }
-deriving stock instance Functor ( D k u ) => Functor ( C k u )
-
--- | \( C^2 \)-differentiable mappings.
-type (~>)  = C 2
--- | \( C^3 \)-differentiable mappings.
-type (~~>) = C 3
-
--- | @D k u v@ is the space of @k@-th order germs of functions from @u@ to @v@,
--- represented by the algebra:
---
--- \[ \mathbb{Z}[x_1, \ldots, x_n]/(x_1, \ldots, x_n)^(k+1) \otimes_\mathbb{Z} v \]
---
--- when @u@ is of dimension @n@.
-type D :: Nat -> Type -> Type -> Type
-type family D k u
-
-type instance D n ( ℝ 0 ) = D𝔸0
-
-type instance D 2 ( ℝ 1 ) = D2𝔸1
-type instance D 2 ( ℝ 2 ) = D2𝔸2
-type instance D 2 ( ℝ 3 ) = D2𝔸3
-type instance D 2 ( ℝ 4 ) = D2𝔸4
-
-type instance D 3 ( ℝ 1 ) = D3𝔸1
-type instance D 3 ( ℝ 2 ) = D3𝔸2
-type instance D 3 ( ℝ 3 ) = D3𝔸3
-type instance D 3 ( ℝ 4 ) = D3𝔸3
-
--- | \( \mathbb{Z} \otimes_\mathbb{Z} v \).
-newtype D𝔸0 v = D0 { v :: v }
-  deriving stock   ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
-  deriving Applicative
-    via Generically1 D𝔸0
-deriving via MonApRep Double D𝔸0
-  instance HasChainRule 2 Double ( ℝ 0 ) D𝔸0
-deriving via MonApRep ( 𝕀 Double ) D𝔸0
-  instance HasChainRule 2 ( 𝕀 Double ) ( 𝕀ℝ 0 ) D𝔸0
-deriving via MonApRep r D𝔸0 r
-  instance Num r => Num ( D𝔸0 r )
-deriving via ( ApAp r D𝔸0 v )
-  instance ( Num r, Module r ( T v ) ) => Module ( D𝔸0 r ) ( D𝔸0 v )
-
--- | \( \mathbb{Z}[x]/(x)^3 \otimes_\mathbb{Z} v \).
-data D2𝔸1 v =
-  D21 { v :: !v
-      , dx :: !( T v )
-      , dxdx :: !( T v )
-      }
-  deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 )
-  deriving Applicative
-    via Generically1 D2𝔸1
-deriving stock instance ( Show v, Show ( T v ) ) => Show ( D2𝔸1 v )
-deriving via MonApRep Double D2𝔸1
-  instance HasChainRule 2 Double ( ℝ 1 ) D2𝔸1
-deriving via MonApRep ( 𝕀 Double ) D2𝔸1
-  instance HasChainRule 2 ( 𝕀 Double ) ( 𝕀ℝ 1 ) D2𝔸1
-deriving via MonApRep r D2𝔸1 r
-  instance Num r => Num ( D2𝔸1 r )
-deriving via ( ApAp r D2𝔸1 v )
-  instance ( Num r, Module r ( T v ) ) => Module ( D2𝔸1 r ) ( D2𝔸1 v )
-
--- | \( \mathbb{Z}[x]/(x)^4 \otimes_\mathbb{Z} v \).
-data D3𝔸1 v =
-  D31 { v :: !v
-      , dx :: !( T v )
-      , dxdx :: !( T v )
-      , dxdxdx :: !( T v )
-      }
-  deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 )
-  deriving Applicative
-    via Generically1 D3𝔸1
-deriving stock instance ( Show v, Show ( T v ) ) => Show ( D3𝔸1 v )
-deriving via MonApRep Double D3𝔸1
-  instance HasChainRule 3 Double ( ℝ 1 ) D3𝔸1
-deriving via MonApRep ( 𝕀 Double ) D3𝔸1
-  instance HasChainRule 3 ( 𝕀 Double ) ( 𝕀ℝ 1 ) D3𝔸1
-deriving via MonApRep r D3𝔸1 r
-  instance Num r => Num ( D3𝔸1 r )
-deriving via ( ApAp r D3𝔸1 v )
-  instance ( Num r, Module r ( T v ) ) => Module ( D3𝔸1 r ) ( D3𝔸1 v )
-
--- | \( \mathbb{Z}[x, y]/(x, y)^3 \otimes_\mathbb{Z} v \).
-data D2𝔸2 v =
-  D22 { v :: !v
-      , dx, dy :: !( T v )
-      , dxdx, dxdy, dydy :: !( T v )
-      }
-  deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 )
-  deriving Applicative
-    via Generically1 D2𝔸2
-deriving stock instance ( Show v, Show ( T v ) ) => Show ( D2𝔸2 v )
-deriving via MonApRep Double D2𝔸2
-  instance HasChainRule 2 Double ( ℝ 2 ) D2𝔸2
-deriving via MonApRep ( 𝕀 Double ) D2𝔸2
-  instance HasChainRule 2 ( 𝕀 Double ) ( 𝕀ℝ 2 ) D2𝔸2
-deriving via MonApRep r D2𝔸2 r
-  instance Num r => Num ( D2𝔸2 r )
-deriving via ( ApAp r D2𝔸2 v )
-  instance ( Num r, Module r ( T v ) ) => Module ( D2𝔸2 r ) ( D2𝔸2 v )
-
--- | \( \mathbb{Z}[x, y]/(x, y)^4 \otimes_\mathbb{Z} v \).
-data D3𝔸2 v =
-  D32 { v :: !v
-      , dx, dy :: !( T v )
-      , dxdx, dxdy, dydy :: !( T v )
-      , dxdxdx, dxdxdy, dxdydy, dydydy :: !( T v )
-      }
-  deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 )
-  deriving Applicative
-    via Generically1 D3𝔸2
-deriving stock instance ( Show v, Show ( T v ) ) => Show ( D3𝔸2 v )
-deriving via MonApRep Double D3𝔸2
-  instance HasChainRule 3 Double ( ℝ 2 ) D3𝔸2
-deriving via MonApRep ( 𝕀 Double ) D3𝔸2
-  instance HasChainRule 3 ( 𝕀 Double ) ( 𝕀ℝ 2 ) D3𝔸2
-deriving via MonApRep r D3𝔸2 r
-  instance Num r => Num ( D3𝔸2 r )
-deriving via ( ApAp r D3𝔸2 v )
-  instance ( Num r, Module r ( T v ) ) => Module ( D3𝔸2 r ) ( D3𝔸2 v )
-
--- | \( \mathbb{Z}[x, y, z]/(x, y, z)^3 \otimes_\mathbb{Z} v \).
-data D2𝔸3 v =
-  D23 { v :: !v
-      , dx, dy, dz :: !( T v )
-      , dxdx, dxdy, dydy, dxdz, dydz, dzdz :: !( T v )
-      }
-  deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 )
-  deriving Applicative
-    via Generically1 D2𝔸3
-deriving stock instance ( Show v, Show ( T v ) ) => Show ( D2𝔸3 v )
-deriving via MonApRep Double D2𝔸3
-  instance HasChainRule 2 Double ( ℝ 3 ) D2𝔸3
-deriving via MonApRep ( 𝕀 Double ) D2𝔸3
-  instance HasChainRule 2 ( 𝕀 Double ) ( 𝕀ℝ 3 ) D2𝔸3
-deriving via MonApRep r D2𝔸3 r
-  instance Num r => Num ( D2𝔸3 r )
-deriving via ( ApAp r D2𝔸3 v )
-  instance ( Num r, Module r ( T v ) ) => Module ( D2𝔸3 r ) ( D2𝔸3 v )
-
--- | \( \mathbb{Z}[x, y, z]/(x, y, z)^4 \otimes_\mathbb{Z} v \).
-data D3𝔸3 v =
-  D33 { v :: !v
-      , dx, dy, dz :: !( T v )
-      , dxdx, dxdy, dydy, dxdz, dydz, dzdz :: !( T v )
-      , dxdxdx, dxdxdy, dxdydy, dydydy,
-         dxdxdz, dxdydz, dxdzdz, dydydz, dydzdz, dzdzdz :: !( T v )
-      }
-  deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 )
-  deriving Applicative
-    via Generically1 D3𝔸3
-deriving stock instance ( Show v, Show ( T v ) ) => Show ( D3𝔸3 v )
-deriving via MonApRep Double D3𝔸3
-  instance HasChainRule 3 Double ( ℝ 3 ) D3𝔸3
-deriving via MonApRep ( 𝕀 Double ) D3𝔸3
-  instance HasChainRule 3 ( 𝕀 Double ) ( 𝕀ℝ 3 ) D3𝔸3
-deriving via MonApRep r D3𝔸3 r
-  instance Num r => Num ( D3𝔸3 r )
-deriving via ( ApAp r D3𝔸3 v )
-  instance ( Num r, Module r ( T v ) ) => Module ( D3𝔸3 r ) ( D3𝔸3 v )
-
--- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^3 \otimes_\mathbb{Z} v \).
-data D2𝔸4 v =
-  D24 { v :: !v
-      , dx, dy, dz, dw :: !( T v )
-      , dxdx, dxdy, dydy, dxdz, dydz, dzdz, dxdw, dydw, dzdw, dwdw :: !( T v )
-      }
-  deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 )
-  deriving Applicative
-    via Generically1 D2𝔸4
-deriving stock instance ( Show v, Show ( T v ) ) => Show ( D2𝔸4 v )
-deriving via MonApRep Double D2𝔸4
-  instance HasChainRule 2 Double ( ℝ 4 ) D2𝔸4
-deriving via MonApRep ( 𝕀 Double ) D2𝔸4
-  instance HasChainRule 2 ( 𝕀 Double ) ( 𝕀ℝ 4 ) D2𝔸4
-deriving via MonApRep r D2𝔸4 r
-  instance Num r => Num ( D2𝔸4 r )
-deriving via ( ApAp r D2𝔸4 v )
-  instance ( Num r, Module r ( T v ) ) => Module ( D2𝔸4 r ) ( D2𝔸4 v )
-
--- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^3 \otimes_\mathbb{Z} v \).
-data D3𝔸4 v =
-  D34 { v :: !v
-      , dx, dy, dz, dw :: !( T v )
-      , dxdx, dxdy, dydy, dxdz, dydz, dzdz, dxdw, dydw, dzdw, dwdw :: !( T v )
-      , dxdxdx, dxdxdy, dxdydy, dydydy,
-         dxdxdz, dxdydz, dxdzdz, dydzdz, dydydz, dzdzdz,
-         dxdxdw, dxdydw, dydydw, dxdzdw, dydzdw, dzdzdw
-      , dxdwdw, dydwdw, dzdwdw, dwdwdw :: !( T v )
-      }
-  deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 )
-  deriving Applicative
-    via Generically1 D3𝔸4
-deriving stock instance ( Show v, Show ( T v ) ) => Show ( D3𝔸4 v )
-deriving via MonApRep Double D3𝔸4
-  instance HasChainRule 3 Double ( ℝ 4 ) D3𝔸4
-deriving via MonApRep ( 𝕀 Double ) D3𝔸4
-  instance HasChainRule 3 ( 𝕀 Double ) ( 𝕀ℝ 4 ) D3𝔸4
-deriving via MonApRep r D3𝔸4 r
-  instance Num r => Num ( D3𝔸4 r )
-deriving via ( ApAp r D3𝔸4 v )
-  instance ( Num r, Module r ( T v ) ) => Module ( D3𝔸4 r ) ( D3𝔸4 v )
-
---------------------------------------------------------------------------------
-
-{-
-type Differentiable :: Extent -> Type -> Constraint
-class
-  ( Module ( I i Double ) ( T ( I i Double ) )
-  , Torsor ( T ( I i u ) ) ( I i u )
-  , Module ( D ( I i u ) ( I i Double ) ) ( D ( I i u ) ( I i ( ℝ 2 ) ) )
-  , Representable ( I i Double ) ( I i u )
-  , Floating ( D ( I i u ) ( I i Double ) )
-  , Applicative ( D ( I i u ) )
-  ) => Differentiable i u
-
-instance
-  ( Module ( I i Double ) ( T ( I i Double ) )
-  , Torsor ( T ( I i u ) ) ( I i u )
-  , Module ( D ( I i u ) ( I i Double ) ) ( D ( I i u ) ( I i ( ℝ 2 ) ) )
-  , Representable ( I i Double ) ( I i u )
-  , Floating ( D ( I i u ) ( I i Double ) )
-  , Applicative ( D ( I i u ) )
-  ) => Differentiable i u
--}
-
---------------------------------------------------------------------------------
-
--- | @Mon k n@ is the set of monomials in @n@ variables of degree less than or equal to @k@.
-type Mon :: Nat -> Nat -> Type
-newtype Mon k n = Mon { monDegs :: Vec n Word } -- sum <= k
-
-zeroMonomial :: forall n k. KnownNat n => Mon k n
-zeroMonomial = Mon
-             $ unsafeCoerce @[ Word ] @( Vec n Word )
-             $ replicate ( fromIntegral $ natVal' @n proxy# ) 0
-
-isZeroMonomial :: Mon k n -> Bool
-isZeroMonomial ( Mon ds ) = isZeroVector ds
-
-
-isZeroVector :: forall i. Vec i Word -> Bool
-isZeroVector VZ = True
-isZeroVector ( 0 `VS` ds ) = isZeroVector ds
-isZeroVector _ = False
-
-isLinear :: Mon k n -> Maybe Word
-isLinear = go 1 . monDegs
-  where
-    go :: forall i. Word -> Vec i Word -> Maybe Word
-    go _ VZ = Nothing
-    go i ( 1 `VS` ds )
-      | isZeroVector ds
-      = Just i
-      | otherwise
-      = Nothing
-    go i ( 0 `VS` ds )
-      = go ( i + 1 ) ds
-    go _ _ = Nothing
-
-type Deg :: ( Type -> Type ) -> Nat
-type Vars :: ( Type -> Type ) -> Nat
-type family Deg  f
-type family Vars f
-
--- | @'MonomialBasis' f@ exhibits @f u@ as a free @r@-module with basis the
--- monomials in @Vars u@ variables, of degree up to (and including) @Deg u@.
---
--- This is used as an accessory class to derive all the algebra, but calls to
--- @monTabulate@ and @monIndex@ should not remain in the optimised Core.
-type MonomialBasis :: ( Type -> Type ) -> Constraint
-class MonomialBasis f where
-  monTabulate :: ( Mon ( Deg f ) ( Vars f ) -> u ) -> f u
-  monIndex :: f u -> Mon ( Deg f ) ( Vars f ) -> u
-
--- | A newtype to hang off instances for functors representable
--- by a monomial basis (see 'MonomialBasis').
-newtype MonApRep r dv w = MonApRep { unMonApRep :: dv w }
-
-
-instance ( Num v, Applicative dr, MonomialBasis dr ) => Num ( MonApRep r dr v ) where
-  (+) = coerce $ liftA2 @dr ( (+) @v )
-  (-) = coerce $ liftA2 @dr ( (-) @v )
-  fromInteger i = MonApRep $ pure @dr $ fromInteger @v i
-
-  MonApRep df1 * MonApRep df2 =
-    MonApRep $
-      monTabulate \ mon ->
-        sum [ monIndex df1 m1 * monIndex df2 m2 | (m1, m2) <- split mon ]
-
-  abs    = error "no abs please"
-  signum = error "no signum please"
-
--- | Newtype for the module instance @Module r v => Module ( dr r ) ( dr v )@.
-newtype ApAp r dr v = ApAp { unApAp :: dr v }
-
-instance ( Num ( dr r ), Module r ( T v ), Applicative dr ) => Module ( dr r ) ( ApAp r dr v ) where
-  ApAp u ^+^ ApAp v = ApAp $ liftA2 ( coerce $ (^+^) @r @( T v ) ) u v
-  ApAp u ^-^ ApAp v = ApAp $ liftA2 ( coerce $ (^-^) @r @( T v ) ) u v
-  origin = ApAp $ pure $ coerce $ origin @r @( T v )
-  k *^ ApAp u = ApAp $ liftA2 ( coerce $ (*^) @r @( T v ) ) k u
-
---------------------------------------------------------------------------------
--- Chain rule and Faà di Bruno formula.
-
-class HasChainRule n r v dn_v where
-  chain   :: Module r ( T w )
-          => D n ( ℝ 1 ) v -> dn_v w -> D n ( ℝ 1 ) w
-  konst   :: Module s ( T w ) => w -> dn_v w
-  value   :: dn_v w -> w
-  linearD :: Module r ( T w ) => ( v -> w ) -> v -> dn_v w
-
-chainRule2 :: forall r v w
-           .  ( HasChainRule 2 r v ( D 2 v ), Module r ( T w ) )
-           => C 2 ( ℝ 1 ) v -> C 2 v w -> C 2 ( ℝ 1 ) w
-chainRule2 ( D df ) ( D dg ) =
-  D \ x ->
-    case df x of
-      df_x@( D21 { v = f_x } ) ->
-        chain @2 @r @v @( D 2 v ) df_x ( dg f_x )
-
-
-uncurryD2 :: D 2 a ~ D 2 ( ℝ 1 )
-          => D 2 ( ℝ 1 ) ( C 2 a b ) -> a -> D 2 ( ℝ 2 ) b
-uncurryD2 ( D21 ( D b_t0 ) ( T ( D dbdt_t0 ) ) ( T ( D d2bdt2_t0 ) ) ) s0 =
-  let !( D21 b_t0s0 dbds_t0s0 d2bds2_t0s0 ) = b_t0 s0
-      !( D21 dbdt_t0s0 d2bdtds_t0s0 _ ) = dbdt_t0 s0
-      !( D21 d2bdt2_t0s0 _ _ ) = d2bdt2_t0 s0
-  in D22 b_t0s0 ( T dbdt_t0s0 ) dbds_t0s0 ( T d2bdt2_t0s0 ) d2bdtds_t0s0 d2bds2_t0s0
-
-
-
--- | Recover the underlying function, discarding all infinitesimal information.
-fun :: forall n r v w. HasChainRule n r v ( D n v ) => C n v w -> ( v -> w )
-fun ( D df ) = value @n @r @v @( D n v ) . df
-{-# INLINE fun #-}
-
--- | The differentiable germ of a coordinate variable.
-var :: forall n r v
-    .  ( Module r ( T r ), Representable r v, HasChainRule n r v ( D n v ) )
-    => Fin ( Dim v ) -> C n v r
-var i = D $ linearD @n @r @v @( D n v ) ( `index` i )
-{-# INLINE var #-}
-
-
-instance forall r dn_v n v
-      .  ( Representable r v
-         , MonomialBasis dn_v
-         , Deg dn_v ~ n
-         , MonomialBasis ( D n ( ℝ 1 ) )
-         , Deg  ( D n ( ℝ 1 ) ) ~ n
-         , Vars ( D n ( ℝ 1 ) ) ~ 1
-         , KnownNat n, KnownNat ( Dim v )
-         , Vars dn_v ~ Dim v
-         ) => HasChainRule n r v ( MonApRep r dn_v ) where
-  chain :: forall w
-          .  Module r ( T w )
-          => D n ( ℝ 1 ) v -> MonApRep r dn_v w -> D n ( ℝ 1 ) w
-  chain = coerce $ monChain @r @v @w @( D n ( ℝ 1 ) ) @dn_v
-  value ( MonApRep u ) = u `monIndex` zeroMonomial
-  {-# INLINE value #-}
-  konst :: forall s w . Module s ( T w ) => w -> MonApRep r dn_v w
-  konst k = MonApRep $ monTabulate \ mon ->
-    if | isZeroMonomial mon
-       -> k
-       | otherwise
-       -> unT $ origin @s @( T w )
-  {-# INLINE konst #-}
-  linearD :: forall w. Module r ( T w ) => ( v -> w ) -> v -> MonApRep r dn_v w
-  linearD f v = MonApRep $ monTabulate \ mon ->
-    if | isZeroMonomial mon
-       -> f v
-       | Just i <- isLinear mon
-       -> f $ tabulate @r @v \ ( Fin j ) ->
-                if | j == i
-                   -> 1 :: r
-                   | otherwise
-                   -> 0 :: r
-       | otherwise
-       -> unT $ origin @r @( T w )
-  {-# INLINE linearD #-}
-
-
-monChain :: forall r v w dr1 dv
-      .  ( Module r ( T w ), Representable r v
-         , MonomialBasis dr1, Vars dr1 ~ 1
-         , MonomialBasis dv , Vars dv  ~ Dim v
-         , Deg dr1 ~ Deg dv
-         , KnownNat ( Dim v )
-         , KnownNat ( Deg dr1 )
-         )
-      => dr1 v -> dv w -> dr1 w
-monChain df dg = monTabulate @dr1
-               $ doChain @r @w ( index @r @v ) ( monIndex @dr1 df ) ( monIndex @dv dg )
-{-# INLINE monChain #-}
-
--- | Compute the chain rule for the composition \( g(f_1(t), .., f_d(t)) \),
--- for derivatives up to order @k@.
-doChain :: forall r w k d v
-        .  ( KnownNat k, KnownNat d, Module r ( T w ) )
-        => ( v -> Fin d -> r )
-        -> ( Mon k 1 -> v )
-        -> ( Mon k d -> w )
-        -> ( Mon k 1 -> w )
-doChain index_v df dg = \ ( Mon ( deg `VS` _ ) ) ->
-  unT $ moduleSum
-    [ T ( dg m_g ) ^*
-    ( sum [ fromIntegral ( faà k is )
-          * product
-              [ product
-                [ ( df ( Mon $ f_deg `VS` VZ ) `index_v` v_index ) ^ pow
-                | ( f_deg, pow ) <- pows ]
-              | ( v_index, pows ) <- zipIndices is ]
-            | is <- multiSubsetsSum [ 1 .. k ] k ( monDegs m_g ) ]
-      )
-    | m <- monomials @_ @d deg
-    , m_g <- subs m
-    ]
-  where
-    k = word @k
-{-# INLINE doChain #-}
-
-moduleSum :: Module r u => [ u ] -> u
-moduleSum = ( foldr (^+^) ) origin
-
--- | Faà di Bruno coefficient.
-faà :: Word -> Vec n [ ( Word, Word ) ] -> Word
-faà k multisubsets =
-  factorial k `div`
-    product [ factorial p * ( factorial i ) ^ p
-            | multisubset <- toList multisubsets
-            , ( i, p ) <- multisubset ]
-
-factorial :: Word -> Word
-factorial i = product [ 1 .. i ]
-
--- | All monomials of degree **exactly** @k@ in @n@ variables,
--- in lexicographic order.
-monomials :: forall k n
-          .  ( KnownNat n )
-          => Word -- ^ degree
-          -> [ Mon k n ]
-monomials k = unsafeCoerce @[ [ Word ] ] @[ Mon k n ]
-            $ monomials' k ( word @n )
-
-monomials' :: Word -> Word -> [ [ Word ] ]
-monomials' k n | k < 0 || n <= 0 = []
-monomials' 0 n = [ replicate ( fromIntegral n ) 0 ]
-monomials' k n = addOne      ( monomials' ( k - 1 ) n )
-              ++ map ( 0 : ) ( monomials' k ( n - 1 ) )
-  where
-    addOne :: [ [ Word ] ] -> [ [ Word ] ]
-    addOne []                   = []
-    addOne ( []         : _   ) = []
-    addOne ( ( d : ds ) : dss ) = ( ( d + 1 ) : ds ) : addOne dss
-
--- | All monomials less than or equal to a given monomial.
-subs :: Mon k n -> [ Mon k n ]
-subs = unsafeCoerce subs'
-
-subs' :: [ Word ] -> [ [ Word ] ]
-subs' [] = [ [] ]
-subs' ( d : ds ) =
-  [ i : as
-  | i <- [ 0 .. d ]
-  , as <- subs' ds ]
-
-split :: Mon k n -> [ ( Mon k n, Mon k n ) ]
-split = unsafeCoerce split'
-
-split' :: [ Word ] -> [ ( [ Word ], [ Word ] ) ]
-split' [] = [ ( [], [] ) ]
-split' ( d : ds ) =
-  [ ( i : as, d - i : bs )
-  | i <- [ 0 .. d ]
-  , (as, bs) <- split' ds ]
-
--- | @multiSubsetsSum is s ns@ computes all collection of multisubsets of @is@,
--- with sizes specifeid by @ns@, such that the total sum is @s@.
-multiSubsetsSum :: forall n
-                .  [ Word ]   -- ^ set to pick from
-                -> Word       -- ^ desired total sum
-                -> Vec n Word -- ^ sizes of each multisets
-                -> [ Vec n [ ( Word, Word ) ] ]
-multiSubsetsSum is = go
-  where
-    go :: forall i. Word -> Vec i Word -> [ Vec i [ ( Word, Word ) ] ]
-    go 0 VZ = [ VZ ]
-    go _ VZ = [ ]
-    go s (n `VS` ns) =
-      [ multi `VS` rest
-      | s_i <- [ n * i_min .. s ]
-      , multi <- multiSubsetSum n s_i is
-      , rest <- go ( s - s_i ) ns ]
-    i_min = case is of
-      [] -> 0
-      _  -> max 0 $ minimum is
-
--- | Computes the multisubsets of the given set which have the specified sum
--- and number of elements.
-multiSubsetSum :: Word -- ^ size of multisubset
-               -> Word -- ^ desired sum
-               -> [ Word ] -- ^ set to pick from
-               -> [ [ ( Word, Word ) ] ]
-multiSubsetSum 0 0 _  = [ [] ]
-multiSubsetSum 0 _ _  = []
-multiSubsetSum _ _ [] = []
-multiSubsetSum n s ( i : is ) =
-  [ if j == 0 then js else ( i, j ) : js
-  | j <- [ 0 .. n ]
-  , js <- multiSubsetSum ( n - j ) ( s - i * j ) is
-  ]
-
-word :: forall n. KnownNat n => Word
-word = fromIntegral $ natVal' @n proxy#
-
-
---------------------------------------------------------------------------------
--- MonomialBasis instances.
-
-type instance Deg  D𝔸0 = 2
-type instance Vars D𝔸0 = 0
-instance MonomialBasis D𝔸0 where
-  monTabulate f = D0 { v }
-     where
-      v = f $ Mon VZ
-  {-# INLINE monTabulate #-}
-  monIndex ( D0 { v } ) = \ _ -> v
-  {-# INLINE monIndex #-}
-
-type instance Deg  D2𝔸1 = 2
-type instance Vars D2𝔸1 = 1
-instance MonomialBasis D2𝔸1 where
-  monTabulate f =
-    D21 { v, dx, dxdx }
-     where
-      v    =     f $ Mon ( 0 `VS` VZ )
-      dx   = T $ f $ Mon ( 1 `VS` VZ )
-      dxdx = T $ f $ Mon ( 2 `VS` VZ )
-  {-# INLINE monTabulate #-}
-
-  monIndex ( D21 { .. } ) = \ case
-    Mon ( 1 `VS` VZ ) -> unT dx
-    Mon ( 2 `VS` VZ ) -> unT dxdx
-    _                 -> v
-  {-# INLINE monIndex #-}
-
-type instance Deg  D3𝔸1 = 3
-type instance Vars D3𝔸1 = 1
-instance MonomialBasis D3𝔸1 where
-  monTabulate f =
-    D31 { v, dx, dxdx, dxdxdx }
-     where
-      v      =     f $ Mon ( 0 `VS` VZ )
-      dx     = T $ f $ Mon ( 1 `VS` VZ )
-      dxdx   = T $ f $ Mon ( 2 `VS` VZ )
-      dxdxdx = T $ f $ Mon ( 3 `VS` VZ )
-  {-# INLINE monTabulate #-}
-
-  monIndex ( D31 { .. } ) = \ case
-    Mon ( 1 `VS` VZ ) -> unT dx
-    Mon ( 2 `VS` VZ ) -> unT dxdx
-    Mon ( 3 `VS` VZ ) -> unT dxdxdx
-    _                 -> v
-  {-# INLINE monIndex #-}
-
-
-type instance Deg  D2𝔸2 = 2
-type instance Vars D2𝔸2 = 2
-instance MonomialBasis D2𝔸2 where
-  monTabulate f = D22 { .. }
-     where
-      v    =     f $ Mon ( 0 `VS` 0 `VS` VZ )
-      dx   = T $ f $ Mon ( 1 `VS` 0 `VS` VZ )
-      dy   = T $ f $ Mon ( 0 `VS` 1 `VS` VZ )
-      dxdx = T $ f $ Mon ( 2 `VS` 0 `VS` VZ )
-      dxdy = T $ f $ Mon ( 1 `VS` 1 `VS` VZ )
-      dydy = T $ f $ Mon ( 0 `VS` 2 `VS` VZ )
-  {-# INLINE monTabulate #-}
-
-  monIndex ( D22 { .. } ) = \ case
-    Mon ( 1 `VS` 0 `VS` VZ ) -> unT dx
-    Mon ( 0 `VS` 1 `VS` VZ ) -> unT dy
-    Mon ( 2 `VS` 0 `VS` VZ ) -> unT dxdx
-    Mon ( 1 `VS` 1 `VS` VZ ) -> unT dxdy
-    Mon ( 0 `VS` 2 `VS` VZ ) -> unT dydy
-    _                        -> v
-  {-# INLINE monIndex #-}
-
-type instance Deg  D3𝔸2 = 3
-type instance Vars D3𝔸2 = 2
-instance MonomialBasis D3𝔸2 where
-  monTabulate f = D32 { .. }
-     where
-      v      =     f $ Mon ( 0 `VS` 0 `VS` VZ )
-      dx     = T $ f $ Mon ( 1 `VS` 0 `VS` VZ )
-      dy     = T $ f $ Mon ( 0 `VS` 1 `VS` VZ )
-      dxdx   = T $ f $ Mon ( 2 `VS` 0 `VS` VZ )
-      dxdy   = T $ f $ Mon ( 1 `VS` 1 `VS` VZ )
-      dydy   = T $ f $ Mon ( 0 `VS` 2 `VS` VZ )
-      dxdxdx = T $ f $ Mon ( 3 `VS` 0 `VS` VZ )
-      dxdxdy = T $ f $ Mon ( 2 `VS` 1 `VS` VZ )
-      dxdydy = T $ f $ Mon ( 1 `VS` 2 `VS` VZ )
-      dydydy = T $ f $ Mon ( 0 `VS` 3 `VS` VZ )
-  {-# INLINE monTabulate #-}
-
-  monIndex ( D32 { .. } ) = \ case
-    Mon ( 1 `VS` 0 `VS` VZ ) -> unT dx
-    Mon ( 0 `VS` 1 `VS` VZ ) -> unT dy
-    Mon ( 2 `VS` 0 `VS` VZ ) -> unT dxdx
-    Mon ( 1 `VS` 1 `VS` VZ ) -> unT dxdy
-    Mon ( 0 `VS` 2 `VS` VZ ) -> unT dydy
-    Mon ( 3 `VS` 0 `VS` VZ ) -> unT dxdxdx
-    Mon ( 2 `VS` 1 `VS` VZ ) -> unT dxdxdy
-    Mon ( 1 `VS` 2 `VS` VZ ) -> unT dxdydy
-    Mon ( 0 `VS` 3 `VS` VZ ) -> unT dydydy
-    _                        -> v
-  {-# INLINE monIndex #-}
-
-type instance Deg  D2𝔸3 = 2
-type instance Vars D2𝔸3 = 3
-instance MonomialBasis D2𝔸3 where
-  monTabulate f = D23 { .. }
-     where
-      v    =     f $ Mon ( 0 `VS` 0 `VS` 0 `VS` VZ )
-      dx   = T $ f $ Mon ( 1 `VS` 0 `VS` 0 `VS` VZ )
-      dy   = T $ f $ Mon ( 0 `VS` 1 `VS` 0 `VS` VZ )
-      dz   = T $ f $ Mon ( 0 `VS` 0 `VS` 1 `VS` VZ )
-      dxdx = T $ f $ Mon ( 2 `VS` 0 `VS` 0 `VS` VZ )
-      dxdy = T $ f $ Mon ( 1 `VS` 1 `VS` 0 `VS` VZ )
-      dydy = T $ f $ Mon ( 0 `VS` 2 `VS` 0 `VS` VZ )
-      dxdz = T $ f $ Mon ( 1 `VS` 0 `VS` 1 `VS` VZ )
-      dydz = T $ f $ Mon ( 0 `VS` 1 `VS` 1 `VS` VZ )
-      dzdz = T $ f $ Mon ( 0 `VS` 0 `VS` 2 `VS` VZ )
-  {-# INLINE monTabulate #-}
-
-
-  monIndex ( D23 { .. } ) = \ case
-    Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) -> unT dx
-    Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) -> unT dy
-    Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) -> unT dz
-    Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) -> unT dxdx
-    Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) -> unT dxdy
-    Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) -> unT dydy
-    Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) -> unT dxdz
-    Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) -> unT dydz
-    Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) -> unT dzdz
-    _                               -> v
-  {-# INLINE monIndex #-}
-
-type instance Deg  D3𝔸3 = 3
-type instance Vars D3𝔸3 = 3
-instance MonomialBasis D3𝔸3 where
-  monTabulate f = D33 { .. }
-     where
-      v      =     f $ Mon ( 0 `VS` 0 `VS` 0 `VS` VZ )
-      dx     = T $ f $ Mon ( 1 `VS` 0 `VS` 0 `VS` VZ )
-      dy     = T $ f $ Mon ( 0 `VS` 1 `VS` 0 `VS` VZ )
-      dz     = T $ f $ Mon ( 0 `VS` 0 `VS` 1 `VS` VZ )
-      dxdx   = T $ f $ Mon ( 2 `VS` 0 `VS` 0 `VS` VZ )
-      dxdy   = T $ f $ Mon ( 1 `VS` 1 `VS` 0 `VS` VZ )
-      dydy   = T $ f $ Mon ( 0 `VS` 2 `VS` 0 `VS` VZ )
-      dxdz   = T $ f $ Mon ( 1 `VS` 0 `VS` 1 `VS` VZ )
-      dydz   = T $ f $ Mon ( 0 `VS` 1 `VS` 1 `VS` VZ )
-      dzdz   = T $ f $ Mon ( 0 `VS` 0 `VS` 2 `VS` VZ )
-      dxdxdx = T $ f $ Mon ( 3 `VS` 0 `VS` 0 `VS` VZ )
-      dxdxdy = T $ f $ Mon ( 2 `VS` 1 `VS` 0 `VS` VZ )
-      dxdydy = T $ f $ Mon ( 1 `VS` 2 `VS` 0 `VS` VZ )
-      dydydy = T $ f $ Mon ( 0 `VS` 3 `VS` 0 `VS` VZ )
-      dxdxdz = T $ f $ Mon ( 2 `VS` 0 `VS` 1 `VS` VZ )
-      dxdydz = T $ f $ Mon ( 1 `VS` 1 `VS` 1 `VS` VZ )
-      dxdzdz = T $ f $ Mon ( 1 `VS` 0 `VS` 2 `VS` VZ )
-      dydydz = T $ f $ Mon ( 0 `VS` 2 `VS` 1 `VS` VZ )
-      dydzdz = T $ f $ Mon ( 0 `VS` 1 `VS` 2 `VS` VZ )
-      dzdzdz = T $ f $ Mon ( 0 `VS` 0 `VS` 3 `VS` VZ )
-  {-# INLINE monTabulate #-}
-
-  monIndex ( D33 { .. } ) = \ case
-    Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) -> unT dx
-    Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) -> unT dy
-    Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) -> unT dz
-    Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) -> unT dxdx
-    Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) -> unT dxdy
-    Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) -> unT dydy
-    Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) -> unT dxdz
-    Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) -> unT dydz
-    Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) -> unT dzdz
-    Mon ( 3 `VS` 0 `VS` 0 `VS` VZ ) -> unT dxdxdx
-    Mon ( 2 `VS` 1 `VS` 0 `VS` VZ ) -> unT dxdxdy
-    Mon ( 1 `VS` 2 `VS` 0 `VS` VZ ) -> unT dxdydy
-    Mon ( 0 `VS` 3 `VS` 0 `VS` VZ ) -> unT dydydy
-    Mon ( 2 `VS` 0 `VS` 1 `VS` VZ ) -> unT dxdxdz
-    Mon ( 1 `VS` 1 `VS` 1 `VS` VZ ) -> unT dxdydz
-    Mon ( 1 `VS` 0 `VS` 2 `VS` VZ ) -> unT dxdzdz
-    Mon ( 0 `VS` 1 `VS` 2 `VS` VZ ) -> unT dydzdz
-    Mon ( 0 `VS` 0 `VS` 3 `VS` VZ ) -> unT dzdzdz
-    _                               -> v
-  {-# INLINE monIndex #-}
-
-type instance Deg  D2𝔸4 = 2
-type instance Vars D2𝔸4 = 4
-instance MonomialBasis D2𝔸4 where
-  monTabulate f = D24 { .. }
-     where
-      v    =     f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 0 `VS` VZ )
-      dx   = T $ f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ )
-      dy   = T $ f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ )
-      dz   = T $ f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ )
-      dw   = T $ f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ )
-      dxdx = T $ f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ )
-      dxdy = T $ f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ )
-      dydy = T $ f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ )
-      dxdz = T $ f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ )
-      dydz = T $ f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ )
-      dzdz = T $ f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ )
-      dxdw = T $ f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ )
-      dydw = T $ f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ )
-      dzdw = T $ f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ )
-      dwdw = T $ f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ )
-  {-# INLINE monTabulate #-}
-
-  monIndex ( D24 { .. } ) = \ case
-    Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> unT dx
-    Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> unT dy
-    Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> unT dz
-    Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> unT dw
-    Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> unT dxdx
-    Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> unT dxdy
-    Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) -> unT dydy
-    Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> unT dxdz
-    Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) -> unT dydz
-    Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) -> unT dzdz
-    Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> unT dxdw
-    Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) -> unT dydw
-    Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) -> unT dzdw
-    Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) -> unT dwdw
-    _                                      -> v
-  {-# INLINE monIndex #-}
-
-type instance Deg  D3𝔸4 = 3
-type instance Vars D3𝔸4 = 4
-instance MonomialBasis D3𝔸4 where
-  monTabulate f = D34 { .. }
-     where
-      v      =     f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 0 `VS` VZ )
-      dx     = T $ f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ )
-      dy     = T $ f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ )
-      dz     = T $ f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ )
-      dw     = T $ f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ )
-      dxdx   = T $ f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ )
-      dxdy   = T $ f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ )
-      dydy   = T $ f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ )
-      dxdz   = T $ f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ )
-      dydz   = T $ f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ )
-      dzdz   = T $ f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ )
-      dxdw   = T $ f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ )
-      dydw   = T $ f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ )
-      dzdw   = T $ f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ )
-      dwdw   = T $ f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ )
-      dxdxdx = T $ f $ Mon ( 3 `VS` 0 `VS` 0 `VS` 0 `VS` VZ )
-      dxdxdy = T $ f $ Mon ( 2 `VS` 1 `VS` 0 `VS` 0 `VS` VZ )
-      dxdydy = T $ f $ Mon ( 1 `VS` 2 `VS` 0 `VS` 0 `VS` VZ )
-      dydydy = T $ f $ Mon ( 0 `VS` 3 `VS` 0 `VS` 0 `VS` VZ )
-      dxdxdz = T $ f $ Mon ( 2 `VS` 0 `VS` 1 `VS` 0 `VS` VZ )
-      dxdydz = T $ f $ Mon ( 1 `VS` 1 `VS` 1 `VS` 0 `VS` VZ )
-      dxdzdz = T $ f $ Mon ( 1 `VS` 0 `VS` 2 `VS` 0 `VS` VZ )
-      dydydz = T $ f $ Mon ( 0 `VS` 2 `VS` 1 `VS` 0 `VS` VZ )
-      dydzdz = T $ f $ Mon ( 0 `VS` 1 `VS` 2 `VS` 0 `VS` VZ )
-      dzdzdz = T $ f $ Mon ( 0 `VS` 0 `VS` 3 `VS` 0 `VS` VZ )
-      dxdxdw = T $ f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 1 `VS` VZ )
-      dxdydw = T $ f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 1 `VS` VZ )
-      dydydw = T $ f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 1 `VS` VZ )
-      dxdzdw = T $ f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 1 `VS` VZ )
-      dydzdw = T $ f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 1 `VS` VZ )
-      dzdzdw = T $ f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 1 `VS` VZ )
-      dxdwdw = T $ f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 2 `VS` VZ )
-      dydwdw = T $ f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 2 `VS` VZ )
-      dzdwdw = T $ f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 2 `VS` VZ )
-      dwdwdw = T $ f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 3 `VS` VZ )
-
-  {-# INLINE monTabulate #-}
-
-  monIndex ( D34 { .. } ) = \ case
-    Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> unT dx
-    Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> unT dy
-    Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> unT dz
-    Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> unT dw
-    Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> unT dxdx
-    Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> unT dxdy
-    Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) -> unT dydy
-    Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> unT dxdz
-    Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) -> unT dydz
-    Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) -> unT dzdz
-    Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> unT dxdw
-    Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) -> unT dydw
-    Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) -> unT dzdw
-    Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) -> unT dwdw
-    Mon ( 3 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> unT dxdxdx
-    Mon ( 2 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> unT dxdxdy
-    Mon ( 1 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) -> unT dxdydy
-    Mon ( 0 `VS` 3 `VS` 0 `VS` 0 `VS` VZ ) -> unT dydydy
-    Mon ( 2 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> unT dxdxdz
-    Mon ( 1 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) -> unT dxdydz
-    Mon ( 1 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) -> unT dxdzdz
-    Mon ( 0 `VS` 2 `VS` 1 `VS` 0 `VS` VZ ) -> unT dydydz
-    Mon ( 0 `VS` 1 `VS` 2 `VS` 0 `VS` VZ ) -> unT dydzdz
-    Mon ( 0 `VS` 0 `VS` 3 `VS` 0 `VS` VZ ) -> unT dzdzdz
-    Mon ( 2 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> unT dxdxdw
-    Mon ( 1 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) -> unT dxdydw
-    Mon ( 0 `VS` 2 `VS` 0 `VS` 1 `VS` VZ ) -> unT dydydw
-    Mon ( 1 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) -> unT dxdzdw
-    Mon ( 0 `VS` 1 `VS` 1 `VS` 1 `VS` VZ ) -> unT dydzdw
-    Mon ( 0 `VS` 0 `VS` 2 `VS` 1 `VS` VZ ) -> unT dzdzdw
-    Mon ( 1 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) -> unT dxdwdw
-    Mon ( 0 `VS` 1 `VS` 0 `VS` 2 `VS` VZ ) -> unT dydwdw
-    Mon ( 0 `VS` 0 `VS` 1 `VS` 2 `VS` VZ ) -> unT dzdwdw
-    Mon ( 0 `VS` 0 `VS` 0 `VS` 3 `VS` VZ ) -> unT dwdwdw
-    _                                      -> v
-  {-# INLINE monIndex #-}

src/splines/Math/Linear/Internal.hs

@@ -0,0 +1,154 @@
+{-# LANGUAGE PolyKinds       #-}
+{-# LANGUAGE TemplateHaskell #-}
+
+
+module Math.Linear.Internal
+  ( ℝ(..)
+  , Fin(..), MFin(..)
+  , RepDim
+  , RepresentableQ(..)
+  , Representable(..), projection, injection
+  )
+  where
+
+-- base
+import Data.Kind
+  ( Type, Constraint )
+import GHC.Generics
+  ( Generic )
+import GHC.TypeNats
+  ( Nat )
+
+-- template-haskell
+import Language.Haskell.TH
+  ( CodeQ )
+
+-- deepseq
+import Control.DeepSeq
+  ( NFData)
+
+--------------------------------------------------------------------------------
+
+-- | Euclidean space \( \mathbb{R}^n \).
+type ℝ :: Nat -> Type
+data family ℝ n
+
+data instance ℝ 0 = ℝ0
+  deriving stock ( Show, Eq, Ord, Generic )
+  deriving anyclass NFData
+newtype instance ℝ 1 = ℝ1 { unℝ1 :: Double }
+  deriving stock   ( Generic )
+  deriving newtype ( Show, Eq, Ord, NFData )
+data    instance ℝ 2 = ℝ2 { _ℝ2_x, _ℝ2_y :: {-# UNPACK #-} !Double }
+  deriving stock    Generic
+  deriving anyclass NFData
+  deriving stock    ( Show, Eq, Ord )
+data    instance ℝ 3 = ℝ3 { _ℝ3_x, _ℝ3_y, _ℝ3_z :: {-# UNPACK #-} !Double }
+  deriving stock    Generic
+  deriving anyclass NFData
+  deriving stock    ( Show, Eq, Ord )
+data    instance ℝ 4 = ℝ4 { _ℝ4_x, _ℝ4_y, _ℝ4_z, _ℝ4_w :: {-# UNPACK #-} !Double }
+  deriving stock    Generic
+  deriving anyclass NFData
+  deriving stock    ( Show, Eq, Ord )
+
+--------------------------------------------------------------------------------
+
+-- | 1, ..., n
+type Fin :: Nat -> Type
+newtype Fin n = Fin Word
+  deriving stock Eq
+
+-- | 0, ..., n
+type MFin :: Nat -> Type
+newtype MFin n = MFin Word
+
+type RepDim :: k -> Nat
+type family RepDim v
+
+type RepresentableQ :: Type -> Type -> Constraint
+class RepresentableQ r v | v -> r where
+  tabulateQ :: ( Fin ( RepDim v ) -> CodeQ r ) -> CodeQ v
+  indexQ    :: CodeQ v -> Fin ( RepDim v ) -> CodeQ r
+
+type Representable :: Type -> Type -> Constraint
+class Representable r v | v -> r where
+  tabulate :: ( Fin ( RepDim v ) -> r ) -> v
+  index    :: v -> Fin ( RepDim v ) -> r
+
+
+projection :: ( Representable r u, Representable r v )
+           => ( Fin ( RepDim v ) -> Fin ( RepDim u ) )
+           -> u -> v
+projection f = \ u ->
+  tabulate \ i -> index u ( f i )
+
+injection :: ( Representable r u, Representable r v )
+          => ( Fin ( RepDim v ) -> MFin ( RepDim u ) )
+          -> u -> v -> v
+injection f = \ u v ->
+  tabulate \ i -> case f i of
+    MFin 0 -> index v i
+    MFin j -> index u ( Fin j )
+
+--------------------------------------------------------------------------------
+
+type instance RepDim ( ℝ n ) = n
+
+instance RepresentableQ Double ( ℝ 0 ) where
+  tabulateQ _ = [|| ℝ0 ||]
+  indexQ _ _ = [|| 0 ||]
+
+instance RepresentableQ Double ( ℝ 1 ) where
+  tabulateQ f = [|| ℝ1 $$( f ( Fin 1 ) ) ||]
+  indexQ p _ = [|| unℝ1 $$p ||]
+
+instance RepresentableQ Double ( ℝ 2 ) where
+  tabulateQ f = [|| ℝ2 $$( f ( Fin 1 ) ) $$( f ( Fin 2 ) ) ||]
+  indexQ p = \ case
+    Fin 1 -> [|| _ℝ2_x $$p ||]
+    _     -> [|| _ℝ2_y $$p ||]
+
+instance RepresentableQ Double ( ℝ 3 ) where
+  tabulateQ f = [|| ℝ3 $$( f ( Fin 1 ) ) $$( f ( Fin 2 ) ) $$( f ( Fin 3 ) ) ||]
+  indexQ p = \ case
+    Fin 1 -> [|| _ℝ3_x $$p ||]
+    Fin 2 -> [|| _ℝ3_y $$p ||]
+    _     -> [|| _ℝ3_z $$p ||]
+
+instance RepresentableQ Double ( ℝ 4 ) where
+  tabulateQ f = [|| ℝ4 $$( f ( Fin 1 ) ) $$( f ( Fin 2 ) ) $$( f ( Fin 3 ) ) $$( f ( Fin 4 ) ) ||]
+  indexQ p = \ case
+    Fin 1 -> [|| _ℝ4_x $$p ||]
+    Fin 2 -> [|| _ℝ4_y $$p ||]
+    Fin 3 -> [|| _ℝ4_z $$p ||]
+    _     -> [|| _ℝ4_w $$p ||]
+
+instance Representable Double ( ℝ 0 ) where
+  tabulate _ = ℝ0
+  index _ _ = 0
+
+instance Representable Double ( ℝ 1 ) where
+  tabulate f = ℝ1 ( f ( Fin 1 ) )
+  index p _ = unℝ1 p
+
+instance Representable Double ( ℝ 2 ) where
+  tabulate f = ℝ2 ( f ( Fin 1 ) ) ( f ( Fin 2 ) )
+  index p = \ case
+    Fin 1 -> _ℝ2_x p
+    _     -> _ℝ2_y p
+
+instance Representable Double ( ℝ 3 ) where
+  tabulate f = ℝ3 ( f ( Fin 1 ) ) ( f ( Fin 2 ) ) ( f ( Fin 3 ) )
+  index p = \ case
+    Fin 1 -> _ℝ3_x p
+    Fin 2 -> _ℝ3_y p
+    _     -> _ℝ3_z p
+
+instance Representable Double ( ℝ 4 ) where
+  tabulate f = ℝ4 ( f ( Fin 1 ) ) ( f ( Fin 2 ) ) ( f ( Fin 3 ) ) ( f ( Fin 4 ) )
+  index p = \ case
+    Fin 1 -> _ℝ4_x p
+    Fin 2 -> _ℝ4_y p
+    Fin 3 -> _ℝ4_z p
+    _     -> _ℝ4_w p

src/splines/Math/Module.hs

@@ -1,7 +1,10 @@
 {-# LANGUAGE AllowAmbiguousTypes  #-}
 {-# LANGUAGE ScopedTypeVariables  #-}
+{-# LANGUAGE TemplateHaskell      #-}
 {-# LANGUAGE UndecidableInstances #-}
 
+{-# OPTIONS_GHC -Wno-orphans #-}
+
 module Math.Module
   ( Module(..), lerp
   , Inner(..), Cross(..)
@@ -9,6 +12,8 @@ module Math.Module
   , norm, squaredNorm, quadrance, distance
   , proj, projC, closestPointOnSegment
   , strictlyParallel, convexCombination
+
+  , ViaModule(..)
   )
   where
 
@@ -18,7 +23,7 @@ import Control.Applicative
 import Control.Monad
   ( guard )
 import Data.Coerce
-  ( Coercible, coerce )
+  ( coerce )
 import Data.Kind
   ( Type, Constraint )
 import Data.Monoid
@@ -32,36 +37,22 @@ import Data.Act
     ( (-->) )
   )
 
--- rounded-hw
+-- groups
-import Numeric.Rounded.Hardware
+import Data.Group
-  ( Rounded(..) )
+  ( Group(..) )
-import Numeric.Rounded.Hardware.Interval.NonEmpty
-  ( Interval(..) )
 
 -- MetaBrush
 import Math.Epsilon
   ( epsilon )
 import Math.Linear
+import Math.Module.Internal
+  ( Module(..), Inner(..) )
+import Math.Ring
+  ( Ring )
+import qualified Math.Ring as Ring
 
 --------------------------------------------------------------------------------
 
-infixl 6 ^+^, ^-^
-infix  9 ^*, *^
-
-class Num r => Module r m | m -> r where
-
-  {-# MINIMAL origin, (^+^), ( (^*) | (*^) ) #-}
-
-  origin :: m
-  (^+^) :: m -> m -> m
-  (^-^) :: m -> m -> m
-  (*^) :: r -> m -> m
-  (^*) :: m -> r -> m
-
-  (*^) = flip (^*)
-  (^*) = flip (*^)
-  m ^-^ n = m ^+^ -1 *^ n
-
 instance ( Applicative f, Module r m ) => Module r ( Ap f m ) where
   origin = pure origin
   (^+^)  = liftA2 (^+^)
@@ -71,11 +62,6 @@ instance ( Applicative f, Module r m ) => Module r ( Ap f m ) where
 lerp :: forall v r p. ( Module r v, Torsor v p ) => r -> p -> p -> p
 lerp t p0 p1 = ( t *^ ( p0 --> p1 :: v ) ) • p0
 
-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
 
@@ -129,42 +115,14 @@ instance ( Torsor ( T u ) u, Module r ( T u ) ) => Interpolatable r u
 
 --------------------------------------------------------------------------------
 
-instance ( Representable r u, Coercible r m, Module r m ) => Module r ( ApRep m u ) where
+instance Ring a => Inner a ( Sum a ) where
-  origin = ApRep $ tabulate @r @u \ _ -> coerce $ origin @r @m
+  Sum a ^.^ Sum b = a Ring.* b
-  {-# INLINE origin #-}
-  ApRep a ^+^ ApRep b = ApRep $ tabulate @r @u \ i ->
-    coerce $ (^+^) @r @m ( coerce ( index @r @u a i ) ) ( coerce ( index @r @u b i ) )
-  {-# INLINE (^+^) #-}
-  ApRep a ^-^ ApRep b = ApRep $ tabulate @r @u \ i ->
-    coerce $ (^-^) @r @m ( coerce ( index @r @u a i ) ) ( coerce ( index @r @u b i ) )
-  {-# INLINE (^-^) #-}
-  k *^ ApRep a = ApRep $ tabulate @r @u \ i ->
-    coerce $ (*^) @r @m k ( coerce ( index @r @u a i ) )
-  {-# INLINE (*^) #-}
-
-deriving via ( ApRep ( Sum Double ) ( ℝ n ) )
-  instance Representable Double ( ℝ n ) => Module Double ( T ( ℝ n ) )
-
-instance Num a => Module a ( Sum a ) where
-
-  origin = Sum 0
-
-  (^+^) = (<>)
-  Sum x ^-^ Sum y = Sum ( x - y )
-
-  c *^ Sum x = Sum ( c * x )
-  Sum x ^* c = Sum ( x * c )
-
-instance Num a => Inner a ( Sum a ) where
-  Sum a ^.^ Sum b = a * b
-
-deriving via Sum Double instance Module Double ( T Double )
 
 instance Inner Double ( T ( ℝ 2 ) ) where
-  V2 x1 y1 ^.^ V2 x2 y2 = x1 * x2 + y1 * y2
+  V2 x1 y1 ^.^ V2 x2 y2 = x1 Ring.* x2 + y1 Ring.* y2
 
 instance Cross Double ( T ( ℝ 2 ) ) where
-  cross ( V2 x1 y1 ) ( V2 x2 y2 ) = x1 * y2 - x2 * y1
+  cross ( V2 x1 y1 ) ( V2 x2 y2 ) = x1 Ring.* y2 Ring.- x2 Ring.* y1
 
 -- | Compute whether two vectors point in the same direction,
 -- that is, whether each vector is a (strictly) positive multiple of the other.
@@ -205,23 +163,62 @@ convexCombination v0 v1 u
     c10 = ( v0 ^-^ v1 ) `cross` u
 
 --------------------------------------------------------------------------------
--- Interval arithmetic using rounded-hw library.
+-- Not sure how to set things up to automate the following...
 
-deriving via Sum ( 𝕀 Double ) instance Module ( 𝕀 Double ) ( T ( 𝕀 Double ) )
+instance Module Double ( T ( ℝ 0 ) ) where
+  origin      = T $$( tabulateQ \ _ -> [|| unT $ origin ||] )
+  T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] )
+  T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] )
+  k *^ T a    = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] )
 
-deriving via ApRep ( Sum ( 𝕀 Double ) ) ( 𝕀ℝ n )
+instance Module Double ( T ( ℝ 1 ) ) where
-  instance Representable ( 𝕀 Double ) ( 𝕀ℝ n ) => Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) )
+  origin      = T $$( tabulateQ \ _ -> [|| unT $ origin ||] )
+  T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] )
+  T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] )
+  k *^ T a    = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] )
 
-instance Inner  ( 𝕀 Double ) ( T ( 𝕀ℝ 2 ) ) where
+instance Module Double ( T ( ℝ 2 ) ) where
-  T ( I ( Rounded ( ℝ2 x1_lo y1_lo ) ) ( Rounded ( ℝ2 x1_hi y1_hi ) ) ) ^.^
+  origin      = T $$( tabulateQ \ _ -> [|| unT $ origin ||] )
-    T ( I ( Rounded ( ℝ2 x2_lo y2_lo ) ) ( Rounded ( ℝ2 x2_hi y2_hi ) ) )
+  T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] )
-      = let !x1x2 = I ( Rounded x1_lo ) ( Rounded x1_hi ) * I ( Rounded x2_lo ) ( Rounded x2_hi )
+  T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] )
-            !y1y2 = I ( Rounded y1_lo ) ( Rounded y1_hi ) * I ( Rounded y2_lo ) ( Rounded y2_hi )
+  k *^ T a    = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] )
-        in x1x2 + y1y2
 
-instance Cross  ( 𝕀 Double ) ( T ( 𝕀ℝ 2 ) ) where
+instance Module Double ( T ( ℝ 3 ) ) where
-  T ( I ( Rounded ( ℝ2 x1_lo y1_lo ) ) ( Rounded ( ℝ2 x1_hi y1_hi ) ) ) `cross`
+  origin      = T $$( tabulateQ \ _ -> [|| unT $ origin ||] )
-    T ( I ( Rounded ( ℝ2 x2_lo y2_lo ) ) ( Rounded ( ℝ2 x2_hi y2_hi ) ) )
+  T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] )
-      = let !x1y2 = I ( Rounded x1_lo ) ( Rounded x1_hi ) * I ( Rounded y2_lo ) ( Rounded y2_hi )
+  T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] )
-            !y2x1 = I ( Rounded x2_lo ) ( Rounded x2_hi ) * I ( Rounded y1_lo ) ( Rounded y1_hi )
+  k *^ T a    = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] )
-        in x1y2 - y2x1
+
+instance Module Double ( T ( ℝ 4 ) ) where
+  origin      = T $$( tabulateQ \ _ -> [|| unT $ origin ||] )
+  T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] )
+  T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] )
+  k *^ T a    = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] )
+
+deriving via ViaModule Double ( T ( ℝ n ) )
+  instance Module Double ( T ( ℝ n ) ) => Semigroup ( T ( ℝ n ) )
+deriving via ViaModule Double ( T ( ℝ n ) )
+  instance Module Double ( T ( ℝ n ) ) => Monoid ( T ( ℝ n ) )
+deriving via ViaModule Double ( T ( ℝ n ) )
+  instance Module Double ( T ( ℝ n ) ) => Group ( T ( ℝ n ) )
+
+deriving via ViaModule Double ( ℝ n )
+  instance Module Double ( T ( ℝ n ) ) => Act ( T ( ℝ n ) ) ( ℝ n )
+deriving via ( ViaModule Double ( ℝ n ) )
+  instance Module Double ( T ( ℝ n ) ) => Torsor ( T ( ℝ n ) ) ( ℝ n )
+
+--------------------------------------------------------------------------------
+
+newtype ViaModule r m = ViaModule { unViaModule :: m }
+
+instance Module r m => Semigroup ( ViaModule r m ) where
+  (<>) = coerce $ (^+^) @r @m
+instance Module r m => Monoid ( ViaModule r m ) where
+  mempty = coerce $ origin @r @m
+instance Module r m => Group ( ViaModule r m ) where
+  invert = coerce $ ( \ x -> (^-^) @r @m ( origin @r @m ) x )
+
+instance ( Semigroup ( T m ), Module r ( T m ) ) => Act    ( T m ) ( ViaModule r m ) where
+  g • ViaModule m = ViaModule ( unT $ g ^+^ T m )
+instance ( Group     ( T m ), Module r ( T m ) ) => Torsor ( T m ) ( ViaModule r m ) where
+  ViaModule a --> ViaModule b = T ( unT $ T b ^-^ T a )

src/splines/Math/Module/Internal.hs

@@ -0,0 +1,90 @@
+{-# LANGUAGE ScopedTypeVariables  #-}
+{-# LANGUAGE TemplateHaskell      #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+module Math.Module.Internal where
+
+-- base
+import Data.Coerce
+  ( coerce )
+import Data.Monoid
+  ( Sum(..) )
+
+-- MetaBrush
+import Math.Ring
+  ( Ring )
+import qualified Math.Ring as Ring
+import Math.Linear
+
+--------------------------------------------------------------------------------
+
+infixl 6 ^+^, ^-^
+infix  9 ^*, *^
+
+class Ring r => Module r m | m -> r where
+
+  {-# MINIMAL origin, (^+^), ( (^*) | (*^) ) #-}
+
+  origin :: m
+  (^+^) :: m -> m -> m
+  (^-^) :: m -> m -> m
+  (*^) :: r -> m -> m
+  (^*) :: m -> r -> m
+
+  (*^) = flip (^*)
+  (^*) = flip (*^)
+  m ^-^ n = m ^+^ Ring.fromInteger ( -1 :: Integer ) *^ n
+
+infixl 8 ^.^
+
+class Module r m => Inner r m where
+  (^.^) :: m -> m -> r
+
+--------------------------------------------------------------------------------
+
+instance Ring a => Module a ( Sum a ) where
+
+  origin = Sum $ Ring.fromInteger ( 0 :: Integer )
+
+  (^+^) = coerce $ (Ring.+) @a
+  (^-^) = coerce $ (Ring.-) @a
+
+  c *^ Sum x = Sum ( c Ring.* x )
+  Sum x ^* c = Sum ( x Ring.* c )
+
+deriving via Sum Double instance Module Double ( T Double )
+
+{-
+moduleViaRepresentable :: Q TH.Type -> Q TH.Type -> TH.DecsQ
+moduleViaRepresentable r m = do
+  i <- TH.newName "i"
+  let tab b = ( ( TH.varE 'tabulateQ `TH.appTypeE` r ) `TH.appTypeE` m ) `TH.appE` TH.lamE [ TH.varP i ] b
+      idx a = ( ( TH.varE 'indexQ `TH.appTypeE` r ) `TH.appTypeE` m ) `TH.appE` TH.varE i
+
+  boo1 <- tab ( TH.varE 'origin `TH.appTypeE` r `TH.appTypeE` m )
+
+
+  let
+
+      meths =
+        [ mkMeth 'origin $ pure boo1 -- $ tab ( TH.varE 'origin `TH.appTypeE` r `TH.appTypeE` m )
+        --, mkMeth '(^+^) $ TH.varE '(^+^) `TH.appE` ( TH ) `TH.appE` ()
+        ]
+
+  i1 <- TH.instanceD ( pure [] ) ( TH.conT ''Module `TH.appT` r `TH.appT` m ) meths
+  return [i1]
+
+  where
+    mkMeth nm body = TH.funD nm [ TH.clause [] ( TH.normalB body ) [] ]
+
+
+moduleViaRepresentable :: ( ( i -> a ) -> TH.CodeQ b ) -> ( b -> i -> TH.CodeQ a ) -> Q TH.Type -> Q TH.Type -> TH.DecsQ
+moduleViaRepresentable tab idx r m = do
+  [d|
+    instance Module $r $m where
+       origin  = $( TH.unTypeCode $ tab \ _ -> origin )
+       a ^+^ b = $( TH.unTypeCode $ tab \ i -> idx [|| a ||] i ^+^ idx [|| b ||] i )
+       a ^-^ b = $( TH.unTypeCode $ tab \ i -> idx [|| a ||] i ^-^ idx [|| b ||] i )
+       k *^ a  = $( TH.unTypeCode $ tab \ i -> [|| {- k *^ -} $$( idx [|| a ||] i ) ||] )
+    |]
+-}

src/splines/Math/Monomial.hs

@@ -0,0 +1,174 @@
+{-# LANGUAGE AllowAmbiguousTypes   #-}
+{-# LANGUAGE QuantifiedConstraints #-}
+{-# LANGUAGE ScopedTypeVariables   #-}
+{-# LANGUAGE TemplateHaskell       #-}
+{-# LANGUAGE UndecidableInstances  #-}
+
+{-# OPTIONS_GHC -ddump-splices -ddump-to-file #-}
+
+module Math.Monomial
+  ( Mon(..)
+  , MonomialBasis(..), Deg, Vars
+  , zeroMonomial, isZeroMonomial, isLinear
+
+  , split, mons
+
+  , faà, multiSubsetsSum, multiSubsetSum
+
+  , prodRuleQ
+
+  ) where
+
+-- base
+import Data.Foldable
+  ( toList )
+import Data.Kind
+  ( Type, Constraint )
+import GHC.Exts
+  ( proxy# )
+import GHC.TypeNats
+  ( KnownNat, Nat, natVal' )
+import Unsafe.Coerce
+  ( unsafeCoerce )
+
+-- template-haskell
+import Language.Haskell.TH
+  ( CodeQ )
+
+-- MetaBrush
+import Math.Linear
+  ( Vec(..), Fin(..) )
+import TH.Utils
+  ( foldQ )
+
+--------------------------------------------------------------------------------
+
+-- | @Mon k n@ is the set of monomials in @n@ variables of degree less than or equal to @k@.
+type Mon :: Nat -> Nat -> Type
+newtype Mon k n = Mon { monDegs :: Vec n Word } -- sum <= k
+  deriving stock Show
+
+type Deg  :: ( Type -> Type ) -> Nat
+type Vars :: ( Type -> Type ) -> Nat
+type family Deg  f
+type family Vars f
+
+zeroMonomial :: forall k n. KnownNat n => Mon k n
+zeroMonomial = Mon $ unsafeCoerce @[ Word ] @( Vec n Word )
+                   $ replicate ( fromIntegral $ word @n ) 0
+
+isZeroMonomial :: Mon k n -> Bool
+isZeroMonomial ( Mon VZ ) = True
+isZeroMonomial ( Mon ( i `VS` is ) )
+  | i == 0
+  = isZeroMonomial ( Mon is )
+  | otherwise
+  = False
+
+isLinear :: Mon k n -> Maybe ( Fin n )
+isLinear = fmap Fin . go 1
+  where
+    go :: forall k' n'. Word -> Mon k' n' -> Maybe Word
+    go _ ( Mon VZ )
+      = Nothing
+    go j ( Mon ( i `VS` is ) )
+      | i == 0
+      = go ( j + 1 ) ( Mon is )
+      | i == 1 && isZeroMonomial ( Mon is )
+      = Just j
+      | otherwise
+      = Nothing
+
+-- | Comultiplication of monomials.
+split :: Mon k n -> [ ( Mon k n, Mon k n ) ]
+split ( Mon VZ ) = [ ( Mon VZ, Mon VZ ) ]
+split ( Mon ( d `VS` ds ) ) =
+      [ ( Mon ( i `VS` as ), Mon ( ( d - i ) `VS` bs ) )
+      | i <- [ 0 .. d ]
+      , ( Mon as, Mon bs ) <- ( split ( Mon ds ) )
+      ]
+
+-- | All monomials of degree less than or equal to @k@ in @n@ variables,
+-- in lexicographic order.
+mons :: forall k n. ( KnownNat n ) => Word -> [ Mon k n ]
+mons k = unsafeCoerce ( mons' k ( word @n ) )
+
+mons' :: Word -> Word -> [ [ Word ] ]
+mons' k _ | k < 0 = []
+mons' _ 0 = [ [] ]
+mons' 0 n = [ replicate ( fromIntegral n ) 0 ]
+mons' k n = [ i : is | i <- reverse [ 0 .. k ], is <- mons' ( k - i ) ( n - 1 ) ]
+
+word :: forall n. KnownNat n => Word
+word = fromIntegral $ natVal' @n proxy#
+
+-- | Faà di Bruno coefficient (naive implementation).
+faà :: Word -> Vec n [ ( Word, Word ) ] -> Word
+faà k multisubsets =
+  factorial k `div`
+    product [ factorial p * ( factorial i ) ^ p
+            | multisubset <- toList multisubsets
+            , ( i, p ) <- multisubset ]
+
+-- | Computes the multisubsets of the given set which have the specified sum
+-- and number of elements.
+multiSubsetSum :: Word -- ^ size of multisubset
+               -> Word -- ^ desired sum
+               -> [ Word ] -- ^ set to pick from
+               -> [ [ ( Word, Word ) ] ]
+multiSubsetSum 0 0 _  = [ [] ]
+multiSubsetSum 0 _ _  = []
+multiSubsetSum _ _ [] = []
+multiSubsetSum n s ( i : is ) =
+  [ if j == 0 then js else ( i, j ) : js
+  | j <- [ 0 .. n ]
+  , js <- multiSubsetSum ( n - j ) ( s - i * j ) is
+  ]
+
+-- | @multiSubsetsSum is s ns@ computes all collection of multisubsets of @is@,
+-- with sizes specified by @ns@, such that the total sum is @s@.
+multiSubsetsSum :: forall n
+                .  [ Word ]   -- ^ set to pick from
+                -> Word       -- ^ desired total sum
+                -> Vec n Word -- ^ sizes of each multisets
+                -> [ Vec n [ ( Word, Word ) ] ]
+multiSubsetsSum is = goMSS
+  where
+    goMSS :: forall i. Word -> Vec i Word -> [ Vec i [ ( Word, Word ) ] ]
+    goMSS 0 VZ = [ VZ ]
+    goMSS _ VZ = [ ]
+    goMSS s (n `VS` ns) =
+      [ multi `VS` rest
+      | s_i <- [ n * i_min .. s ]
+      , multi <- multiSubsetSum n s_i is
+      , rest <- goMSS ( s - s_i ) ns ]
+    i_min = case is of
+      [] -> 0
+      _  -> max 0 $ minimum is
+
+-- | The factorial function \( n! = n \cdot (n-1) \cdot `ldots` \cdot 2 `cdot` 1 \).
+factorial :: Word -> Word
+factorial i = product [ 1 .. i ]
+
+--------------------------------------------------------------------------------
+
+-- | @'MonomialBasis' f@ exhibits @f u@ as a free @r@-module with basis the
+-- monomials in @Vars u@ variables, of degree up to (and including) @Deg u@.
+type MonomialBasis :: ( Type -> Type ) -> Constraint
+class MonomialBasis f where
+  monTabulate :: ( ( Mon ( Deg f ) ( Vars f ) ) -> CodeQ u ) -> CodeQ ( f u )
+  monIndex :: CodeQ ( f u ) -> Mon ( Deg f ) ( Vars f ) -> CodeQ u
+
+--------------------------------------------------------------------------------
+
+prodRuleQ :: forall f r. MonomialBasis f
+          => CodeQ r -> CodeQ ( r -> r -> r ) -> CodeQ ( r -> r -> r )
+            -- Ring r constraint (circumvent TH constraint problem)
+          -> CodeQ ( f r ) -> CodeQ ( f r ) -> CodeQ ( f r )
+prodRuleQ zero plus times df1 df2 =
+  monTabulate @f \ mon ->
+    [|| $$( foldQ plus zero
+           [ [|| $$times $$( monIndex @f df1 m1 ) $$( monIndex @f df2 m2 ) ||]
+           | (m1, m2) <- split mon
+           ] )
+    ||]

src/splines/Math/Ring.hs

@@ -0,0 +1,177 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+
+module Math.Ring
+  ( AbelianGroup(..), Signed(..), Ring(..), Field(..), Transcendental(..)
+
+  , ViaPrelude(..), ViaAbelianGroup(..)
+
+  , ifThenElse
+  )
+  where
+
+-- base
+import Prelude ( Num, Fractional )
+import Prelude hiding ( Num(..), Fractional(..) )
+import qualified Prelude
+import Control.Applicative
+  ( liftA2 )
+import Data.Coerce
+  ( coerce )
+import Data.Monoid
+  ( Ap(..) )
+
+-- groups
+import Data.Group
+  ( Group(invert) )
+
+--------------------------------------------------------------------------------
+
+infixl 6 +, -
+infixl 7 *, /
+infixr 8 ^
+
+class AbelianGroup a where
+  {-# MINIMAL (+), fromInteger, ( (-) | negate ) #-}
+
+  (+) :: a -> a -> a
+  (-) :: a -> a -> a
+  negate :: a -> a
+  fromInteger :: Integer -> a
+
+  a - b = a + negate b
+  negate b = fromInteger 0 - b
+
+class AbelianGroup a => Signed a where
+  abs    :: a -> a
+  signum :: a -> a
+
+class AbelianGroup a => Ring a where
+  (*) :: a -> a -> a
+  (^) :: a -> Word -> a
+  (^) = defaultPow
+
+class Ring a => Field a where
+  {-# MINIMAL fromRational, ( recip | (/) ) #-}
+
+  fromRational :: Rational -> a
+  (/) :: a -> a -> a
+  recip :: a -> a
+
+  recip x = fromInteger 1 / x
+  x / y = x * recip y
+
+class Field a => Transcendental a where
+  pi   :: a
+  cos  :: a -> a
+  sin  :: a -> a
+--  sqrt :: a -> a
+
+--------------------------------------------------------------------------------
+
+ifThenElse :: Bool -> a -> a -> a
+ifThenElse b x y = if b then x else y
+
+--------------------------------------------------------------------------------
+-- A default power function, taken from base (GHC.Real).
+--
+-- DO NOT use with interval arithmetic!
+
+{-# RULES
+"defaultPow 2" forall x. defaultPow x 2 = x*x
+"defaultPow 3" forall x. defaultPow x 3 = x*x*x
+"defaultPow 4" forall x. defaultPow x 4 = let u = x*x in u*u
+"defaultPow 5" forall x. defaultPow x 5 = let u = x*x in u*u*x
+  #-}
+
+{-# INLINE [1] defaultPow #-}
+defaultPow :: Ring a => a -> Word -> a
+defaultPow x0 y0
+  | y0 < 0    = errorWithoutStackTrace "Negative exponent"
+  | y0 == 0   = fromInteger 1
+  | otherwise = powImpl x0 y0
+
+powImpl :: Ring a => a -> Word -> a
+-- powImpl : x0 ^ y0 = (x ^ y)
+powImpl x y
+  | even y    = powImpl (x * x) (y `quot` 2)
+  | y == 1    = x
+  | otherwise = powImplAcc (x * x) (y `quot` 2) x
+
+{-# INLINABLE powImplAcc #-}
+powImplAcc :: Ring a => a -> Word -> a -> a
+-- powImplAcc : x0 ^ y0 = (x ^ y) * z
+powImplAcc x y z
+  | even y    = powImplAcc (x * x) (y `quot` 2) z
+  | y == 1    = x * z
+  | otherwise = powImplAcc (x * x) (y `quot` 2) (x * z)
+
+--------------------------------------------------------------------------------
+
+newtype ViaPrelude a = ViaPrelude { viaPrelude :: a }
+
+instance Num a => AbelianGroup ( ViaPrelude a ) where
+  (+) = coerce $ (Prelude.+) @a
+  (-) = coerce $ (Prelude.-) @a
+  fromInteger = coerce $ Prelude.fromInteger @a
+
+instance Num a => Signed ( ViaPrelude a ) where
+  abs    = coerce $ Prelude.abs    @a
+  signum = coerce $ Prelude.signum @a
+
+instance Num a => Ring ( ViaPrelude a ) where
+  (*) = coerce $ (Prelude.*) @a
+  (^) = coerce $ (Prelude.^) @a @Word
+
+instance Fractional a => Field ( ViaPrelude a ) where
+  fromRational = coerce $ Prelude.fromRational @a
+  (/) = coerce $ (Prelude./) @a
+
+instance Floating a => Transcendental ( ViaPrelude a ) where
+  pi = coerce $ Prelude.pi @a
+  sin = coerce $ Prelude.sin @a
+  cos = coerce $ Prelude.cos @a
+--  sqrt = coerce $ Prelude.sqrt @a
+
+--------------------------------------------------------------------------------
+
+newtype ViaAbelianGroup a = ViaAbelianGroup { viaAbelianGroup :: a }
+
+instance AbelianGroup a => Semigroup ( ViaAbelianGroup a ) where
+  (<>) = coerce $ (+) @a
+instance AbelianGroup a => Monoid    ( ViaAbelianGroup a ) where
+  mempty = ViaAbelianGroup $ fromInteger 0
+instance AbelianGroup a => Group     ( ViaAbelianGroup a ) where
+  invert = coerce $ negate @a
+
+--------------------------------------------------------------------------------
+
+deriving via ViaPrelude Integer instance AbelianGroup   Integer
+deriving via ViaPrelude Integer instance Ring           Integer
+deriving via ViaPrelude Integer instance Signed         Integer
+
+deriving via ViaPrelude Int     instance AbelianGroup   Int
+deriving via ViaPrelude Int     instance Ring           Int
+deriving via ViaPrelude Int     instance Signed         Int
+
+deriving via ViaPrelude Word    instance AbelianGroup   Word
+deriving via ViaPrelude Word    instance Ring           Word
+
+deriving via ViaPrelude Float   instance AbelianGroup   Float
+deriving via ViaPrelude Float   instance Ring           Float
+deriving via ViaPrelude Float   instance Signed         Float
+deriving via ViaPrelude Float   instance Field          Float
+deriving via ViaPrelude Float   instance Transcendental Float
+
+deriving via ViaPrelude Double  instance AbelianGroup   Double
+deriving via ViaPrelude Double  instance Ring           Double
+deriving via ViaPrelude Double  instance Signed         Double
+deriving via ViaPrelude Double  instance Field          Double
+deriving via ViaPrelude Double  instance Transcendental Double
+
+--------------------------------------------------------------------------------
+
+instance ( AbelianGroup r, Applicative f ) => AbelianGroup ( Ap f r ) where
+  (+) = liftA2 $ (+) @r
+  (-) = liftA2 $ (-) @r
+  negate = fmap $ negate @r
+  fromInteger = pure . fromInteger @r

src/splines/TH/Utils.hs

@@ -0,0 +1,23 @@
+{-# LANGUAGE TemplateHaskell #-}
+
+module TH.Utils where
+
+-- template-haskell
+import Language.Haskell.TH
+  ( CodeQ )
+
+-- MetaBrush
+import Math.Ring ( Ring )
+import qualified Math.Ring as Ring
+
+--------------------------------------------------------------------------------
+
+foldQ :: CodeQ ( a -> a -> a ) -> CodeQ a -> [ CodeQ a ] -> CodeQ a
+foldQ _ a0 []  = a0
+foldQ _ _  [a] = a
+foldQ f a0 (a:as) = [|| $$f $$a $$( foldQ f a0 as ) ||]
+
+powQ :: Ring a => CodeQ a -> Word -> CodeQ a
+powQ _ 0 = [|| Ring.fromInteger ( 1 :: Integer ) ||]
+powQ x 1 = x
+powQ x n = [|| $$x Ring.^ ( n :: Word ) ||]