From c7cd6c2a1cb950a75d013d484dfe94b457bac02b Mon Sep 17 00:00:00 2001 From: sheaf Date: Mon, 16 Jan 2023 02:31:31 +0100 Subject: [PATCH] WIP on monomial bases --- MetaBrush.cabal | 1 - src/splines/Math/Bezier/Stroke.hs | 111 ++- src/splines/Math/Linear.hs | 85 ++- src/splines/Math/Linear/Dual.hs | 1153 ++++++++++++++++++++--------- 4 files changed, 944 insertions(+), 406 deletions(-) diff --git a/MetaBrush.cabal b/MetaBrush.cabal index 6304b85..84ba0ee 100644 --- a/MetaBrush.cabal +++ b/MetaBrush.cabal @@ -95,7 +95,6 @@ common common -fexcess-precision -fspecialise-aggressively -optc-O3 - -optc-ffast-math -Wall -Wcompat -fwarn-missing-local-signatures diff --git a/src/splines/Math/Bezier/Stroke.hs b/src/splines/Math/Bezier/Stroke.hs index f8c0277..6b6faa1 100644 --- a/src/splines/Math/Bezier/Stroke.hs +++ b/src/splines/Math/Bezier/Stroke.hs @@ -26,7 +26,7 @@ import Control.Arrow import Control.Applicative ( Applicative(..) ) import Control.Monad - ( unless ) + ( guard, unless ) import Control.Monad.ST ( RealWorld, ST ) import Data.Bifunctor @@ -34,7 +34,7 @@ import Data.Bifunctor import Data.Coerce ( Coercible, coerce ) import Data.Foldable - ( for_ ) + ( for_, toList ) import Data.Functor.Identity ( Identity(..) ) import Data.List.NonEmpty @@ -82,6 +82,13 @@ import Data.Generics.Internal.VL import qualified Control.Parallel.Strategies as Strats ( rdeepseq, parTuple2, using ) +-- rounded-hw +import Numeric.Rounded.Hardware + ( Rounded(..) ) +import Numeric.Rounded.Hardware.Interval.NonEmpty + ( Interval(I) ) +import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval + -- transformers import Control.Monad.Trans.Class ( lift ) @@ -122,7 +129,7 @@ import Math.Roots import Math.Linear import Math.Linear.Dual ---import Debug.Utils +import Debug.Utils -------------------------------------------------------------------------------- @@ -482,7 +489,6 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv = path :: ℝ 1 ~> ℝ 2 ( path, usedParams ) = pathAndUsedParams @Point id -{- curvesI :: 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 'Interval ) curvesI = brushStrokeData @'Interval @brushParams pathI @@ -492,7 +498,6 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv = usedParamsI :: 𝕀ℝ 1 ~> 𝕀 usedParams pathI :: 𝕀ℝ 1 ~> 𝕀ℝ 2 ( pathI, usedParamsI ) = pathAndUsedParams @'Interval singleton --} fwdBwd :: OutlineFn fwdBwd t @@ -517,7 +522,14 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv = $ runD ( brushFromParams @Point proxy# id ) $ toBrushParams params_t - in fwdBwd + bisSols = bisection 0.0001 curvesI + + in trace + ( unlines $ + ( "bisectionMethod: #(possible zeroes) = " ++ show ( length bisSols ) ) : + "" : + map show bisSols ) + fwdBwd ----------------------------------- -- Various utility functions @@ -844,13 +856,14 @@ envelopeEquation :: forall i , Fractional ( I i Double ) ) => D ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) - -> ( I i Double, T ( I i ( ℝ 2 ) ), I i Double, I i Double ) + -> ( I i Double, T ( I i ( ℝ 2 ) ), T ( I i ( ℝ 2 ) ), I i Double, I i Double ) envelopeEquation ( D2 _ dcdt dcds d2cdt2 d2cdtds d2cds2 ) = let ee = dcdt `cross` dcds dEdt = d2cdt2 `cross` dcds + dcdt `cross` d2cdtds dEds = d2cdtds `cross` dcds + dcdt `cross` d2cds2 - tot = dcdt ^-^ ( dEdt / dEds ) *^ dcds - in ( ee, tot, dEdt, dEds ) + tot = dcdt -- ^-^ ( dEdt / dEds ) *^ dcds + dEdsTot = dEds *^ dcdt ^-^ dEdt *^ dcds + in ( ee, tot, dEdsTot, dEdt, dEds ) -- Computation of total derivative dc/dt: -- -- dc/dt = ∂c/∂t + ∂c/∂s ∂s/∂t @@ -1072,7 +1085,7 @@ brushStrokeData path params brush = mkStrokeDatum dpath_t dbrush_t s = let dbrush_t_s = dbrush_t s dstroke@( D2 _c _𝛿c𝛿t _𝛿c𝛿s _ _ _ ) = brushStroke dpath_t dbrush_t_s - ( ee, dcdt, 𝛿E𝛿t, 𝛿E𝛿s ) = envelopeEquation @i dstroke + ( ee, dcdt, 𝛿E𝛿sdcdt, 𝛿E𝛿t, 𝛿E𝛿s ) = envelopeEquation @i dstroke in -- trace -- ( unlines -- [ "envelopeEquation:" @@ -1089,7 +1102,7 @@ brushStrokeData path params brush = { dpath = dpath_t , dbrush = dbrush_t_s , dstroke - , ee, dcdt, 𝛿E𝛿t, 𝛿E𝛿s } + , ee, dcdt, 𝛿E𝛿sdcdt, 𝛿E𝛿t, 𝛿E𝛿s } -- | The value and derivative of a brush stroke at a given coordinate @@ -1119,9 +1132,83 @@ data StrokeDatum i -- \[ \left ( \frac{\mathrm{d} c}{\mathrm{d} t} \right )_{(t_0,s_0)}. \] -- -- This is ill-defined when \( \frac{\partial E}{\partial s} = 0 \). - , dcdt :: T ( I i ( ℝ 2 ) ) + , dcdt, 𝛿E𝛿sdcdt :: T ( I i ( ℝ 2 ) ) + } deriving stock instance Show ( StrokeDatum 'Point ) deriving stock instance Show ( StrokeDatum 'Interval ) + +-------------------------------------------------------------------------------- + +bisection :: Double + -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 'Interval ) ) + -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1, 𝕀 Double, 𝕀ℝ 2 ) ] +bisection minWidth eqs = bisect initialCands [] [] + where + + bisect :: [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1, 𝕀 Double, 𝕀ℝ 2 ) ] -- have solutions, need bisection to refine + -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] -- have been bisected, don't know if they contain solutions yet + -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1, 𝕀 Double, 𝕀ℝ 2 ) ] -- have solutions, don't bisect further + -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1, 𝕀 Double, 𝕀ℝ 2 ) ] + bisect [] [] sols = sols + bisect cands ( ( t, i, s ) : toTry ) sols + | Just ( ee, 𝛿E𝛿sdcdt ) <- isCand t i s + = bisect ( ( t, i, s, ee, 𝛿E𝛿sdcdt ) : cands ) toTry sols + | otherwise + = bisect cands toTry sols + + bisect ( cand@( t@( I ( Rounded ( ℝ1 t_lo ) ) ( Rounded ( ℝ1 t_hi ) ) ) + , i + , s@( I ( Rounded ( ℝ1 s_lo ) ) ( Rounded ( ℝ1 s_hi ) ) ) + , _, _ + ) : cands ) + toTry + sols + -- If the box is small, don't bisect it further, and store it as a candidate solution. + | t_hi - t_lo < minWidth && s_hi - s_lo < minWidth + = trace ( "bisection sol: " ++ show cand ++ "\nnbCands = " ++ show ( length cands ) ++ "\nnbToTry = " ++ show ( length toTry ) ) + $ bisect cands toTry ( cand : sols ) + -- Otherwise, bisect in its longest direction and add the two resulting + -- boxes to the list of boxes to try. + | otherwise + = let newToTry + | t_hi - t_lo > s_hi - s_lo + , let t_mid = 0.5 * ( t_lo + t_hi ) + = ( I ( Rounded ( ℝ1 t_lo ) ) ( Rounded ( ℝ1 t_mid ) ), i, s ) + : ( I ( Rounded ( ℝ1 t_mid ) ) ( Rounded ( ℝ1 t_hi ) ), i, s ) + : toTry + | let s_mid = 0.5 * ( s_lo + s_hi ) + = ( t, i, I ( Rounded ( ℝ1 s_lo ) ) ( Rounded ( ℝ1 s_mid ) ) ) + : ( t, i, I ( Rounded ( ℝ1 s_mid ) ) ( Rounded ( ℝ1 s_hi ) ) ) + : toTry + in bisect cands newToTry sols + + initialCands = + getCands + ( I ( Rounded $ ℝ1 0 ) ( Rounded $ ℝ1 1 ) ) + ( I ( Rounded $ ℝ1 0 ) ( Rounded $ ℝ1 1 ) ) + + getCands t s = + [ (t, i, s, ee, 𝛿E𝛿sdcdt ) + | let !eqs_t = eqs t + , ( eq_t, i ) <- zip ( toList eqs_t ) ( [0,1..] :: [Int] ) + , let !( StrokeDatum { ee, 𝛿E𝛿sdcdt = T 𝛿E𝛿sdcdt } ) = eq_t s + , Interval.inf ee < 0 && Interval.sup ee > 0 + , cmpℝ2 (<) ( getRounded ( Interval.inf 𝛿E𝛿sdcdt ) ) ( ℝ2 0 0 ) + && cmpℝ2 (>) ( getRounded ( Interval.sup 𝛿E𝛿sdcdt ) ) ( ℝ2 0 0 ) + ] + + isCand :: 𝕀ℝ 1 -> Int -> 𝕀ℝ 1 -> Maybe ( 𝕀 Double, 𝕀ℝ 2 ) + isCand t i s = case ( ( eqs t ) `Seq.index` i ) s of + StrokeDatum { ee, 𝛿E𝛿sdcdt = T 𝛿E𝛿sdcdt } -> + do guard $ + Interval.inf ee < 0 && Interval.sup ee > 0 + && cmpℝ2 (<) ( getRounded ( Interval.inf 𝛿E𝛿sdcdt ) ) ( ℝ2 0 0 ) + && cmpℝ2 (>) ( getRounded ( Interval.sup 𝛿E𝛿sdcdt ) ) ( ℝ2 0 0 ) + return ( ee, 𝛿E𝛿sdcdt ) + +cmpℝ2 :: (Double -> Double -> Bool) -> ℝ 2 -> ℝ 2 -> Bool +cmpℝ2 cmp ( ℝ2 x1 y1 ) ( ℝ2 x2 y2 ) + = cmp x1 x2 && cmp y1 y2 diff --git a/src/splines/Math/Linear.hs b/src/splines/Math/Linear.hs index c0de5f8..a78bd73 100644 --- a/src/splines/Math/Linear.hs +++ b/src/splines/Math/Linear.hs @@ -10,10 +10,10 @@ module Math.Linear -- * Points and vectors (second version) , ℝ(..), T(.., V2, V3) - , Fin(..), eqFin, MFin(..) + , Fin(..), MFin(..) , Dim, Representable(..), ApRep(..) , injection, projection - , Vec(..), (!), find + , Vec(..), (!), find, zipIndices -- * Intervals , 𝕀, 𝕀ℝ, singleton, nonDecreasing @@ -27,10 +27,6 @@ import Data.Kind ( Type, Constraint ) import Data.Monoid ( Sum(..) ) -import GHC.Exts - ( TYPE, RuntimeRep(..) - , Word#, plusWord#, minusWord#, isTrue#, eqWord# - ) import GHC.Generics ( Generic, Generic1, Generically(..), Generically1(..) ) import GHC.TypeNats @@ -99,6 +95,11 @@ data instance ℝ 3 = ℝ3 {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-# 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 ) @@ -151,17 +152,13 @@ pattern V3 x y z = T ( ℝ3 x y z ) -------------------------------------------------------------------------------- -- | 1, ..., n -type Fin :: Nat -> TYPE WordRep -newtype Fin n = Fin Word# - -{-# INLINE eqFin #-} -eqFin :: Fin n -> Fin n -> Bool -eqFin ( Fin i ) ( Fin j ) = isTrue# ( i `eqWord#` j ) +type Fin :: Nat -> Type +newtype Fin n = Fin Word + deriving stock Eq -- | 0, ..., n -type MFin :: Nat -> TYPE WordRep -newtype MFin n = MFin Word# - +type MFin :: Nat -> Type +newtype MFin n = MFin Word type Dim :: k -> Nat type family Dim v @@ -181,26 +178,36 @@ instance Representable Double ( ℝ 0 ) where instance Representable Double ( ℝ 1 ) where {-# INLINE tabulate #-} - tabulate f = ℝ1 ( f ( Fin 1## ) ) + 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## ) ) + tabulate f = ℝ2 ( f ( Fin 1 ) ) ( f ( Fin 2 ) ) {-# INLINE index #-} index ( ℝ2 x y ) = \ case - Fin 1## -> x - _ -> y + Fin 1 -> x + _ -> y instance Representable Double ( ℝ 3 ) where {-# INLINE tabulate #-} - tabulate f = ℝ3 ( f ( Fin 1## ) ) ( f ( Fin 2## ) ) ( f ( Fin 3## ) ) + 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 + 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 ) @@ -215,30 +222,42 @@ injection :: ( Representable r u, Representable r v ) -> u -> v -> v injection f = \ u v -> tabulate \ i -> case f i of - MFin 0## -> index v i - MFin j -> index u ( Fin j ) + MFin 0 -> index v i + MFin j -> index u ( Fin j ) -type Vec :: Nat -> TYPE WordRep -> Type +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 +deriving stock instance Functor ( Vec n ) +deriving stock instance Foldable ( Vec n ) +deriving stock instance Traversable ( Vec n ) + infixl 9 ! (!) :: forall l a. Vec l a -> Fin l -> a -VS a _ ! Fin 1## = a -VS _ a ! Fin i = a ! Fin ( i `minusWord#` 1## ) -_ ! _ = error "impossible: Fin 0 is uninhabited" +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 eq v b = MFin ( go 1## v ) +find eq v b = MFin ( go 1 v ) where - go :: Word# -> Vec n a -> Word# + go :: Word -> Vec n a -> Word go j ( VS a as ) | a `eq` b = j | otherwise - = go ( j `plusWord#` 1## ) as - go _ VZ = 0## + = go ( j + 1 ) as + go _ VZ = 0 + +zipIndices :: forall n a. Vec n a -> [ ( Fin n, a ) ] +zipIndices = go 0 + where + go :: forall i. Word -> Vec i a -> [ ( Fin n, a ) ] + go _ VZ = [] + go i (a `VS` as) = ( Fin i, a ) : go ( i + 1 ) as -------------------------------------------------------------------------------- -- Instances in terms of representable. diff --git a/src/splines/Math/Linear/Dual.hs b/src/splines/Math/Linear/Dual.hs index 18d5c9b..6afba83 100644 --- a/src/splines/Math/Linear/Dual.hs +++ b/src/splines/Math/Linear/Dual.hs @@ -10,11 +10,18 @@ module Math.Linear.Dual where import Control.Applicative ( liftA2 ) import Data.Coerce - ( 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 @@ -33,385 +40,217 @@ import Math.Linear -------------------------------------------------------------------------------- --- | Differentiable mappings between spaces. -infixr 0 ~> -type (~>) :: Type -> Type -> Type -newtype u ~> v = D { runD :: u -> D u v } -deriving stock instance Functor ( D u ) => Functor ( (~>) u ) +-- | @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 ) -instance ( Applicative ( D u ), Module r ( T v ) ) - => Module r ( T ( 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 +-- | \( C^2 \)-differentiable mappings. +type (~>) = C 2 +-- | \( C^3 \)-differentiable mappings. +type (~~>) = C 3 --- | @D ( ℝ n ) v@ is \( \mathbb{R}[x_1, \ldots, x_n]/(x_1, \ldots, x_n)^3 \otimes_\mathbb{R} v \) -type D :: Type -> Type -> Type -type family D u -type instance D ( ℝ 0 ) = Dℝ0 -type instance D ( ℝ 1 ) = Dℝ1 -type instance D ( ℝ 2 ) = Dℝ2 -type instance D ( ℝ 3 ) = Dℝ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 ( 𝕀 u ) = D u +type instance D n ( ℝ 0 ) = D𝔸0 -newtype Dℝ0 v = D0 { v :: v } +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 newtype ( Num, Fractional, Floating ) deriving Applicative - via Generically1 Dℝ0 -data Dℝ1 v = D1 { v :: !v, df :: !( T v ), d²f :: !( T v ) } + 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 Dℝ1 -deriving stock instance ( Show v, Show ( T v ) ) => Show ( Dℝ1 v ) -data Dℝ2 v = D2 { v :: !v, dx, dy :: !( T v ), ddx, dxdy, ddy :: !( T v ) } + 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 Dℝ2 -deriving stock instance ( Show v, Show ( T v ) ) => Show ( Dℝ2 v ) -data Dℝ3 v = D3 { v :: !v, dx, dy, dz :: !( T v ), ddx, dxdy, ddy, dxdz, dydz, ddz :: !( T v ) } + 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 Dℝ3 -deriving stock instance ( Show v, Show ( T v ) ) => Show ( Dℝ3 v ) + 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 ) -instance ( Module s ( T r ), Num r ) => Num ( Dℝ1 r ) where - (+) = liftA2 (+) - (-) = liftA2 (-) - negate = fmap negate - fromInteger = konst @Double @( ℝ 1 ) . fromInteger +-- | \( \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 ) - abs = error "no" - signum = error "no" +-- | \( \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 ) - D1 v1 ( T dx1 ) ( T ddx1 ) - * D1 v2 ( T dx2 ) ( T ddx2 ) - = D1 ( v1 * v2 ) - ( T $ dx1 * v2 + v1 * dx2 ) - ( T $ dx1 * dx2 + v1 * ddx2 + ddx1 * v2 ) +-- | \( \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 ) -instance ( Num r, Module s ( T r ) ) => Num ( Dℝ2 r ) where - (+) = liftA2 (+) - (-) = liftA2 (-) - negate = fmap negate - fromInteger = konst @Double @( ℝ 2 ) . fromInteger +-- | \( \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 ) - abs = error "no" - signum = error "no" - - D2 v1 ( T dx1 ) ( T dy1 ) ( T ddx1 ) ( T dxdy1 ) ( T ddy1 ) - * D2 v2 ( T dx2 ) ( T dy2 ) ( T ddx2 ) ( T dxdy2 ) ( T ddy2 ) - = D2 ( v1 * v2 ) - ( T $ dx1 * v2 + v1 * dx2 ) - ( T $ dy1 * v2 + v1 * dy2 ) - ( T $ dx1 * dx2 + v1 * ddx2 + ddx1 * v2 ) - ( T $ dy1 * dx2 + dx1 * dy2 + v1 * dxdy2 + dxdy1 * v2 ) - ( T $ dy1 * dy2 + v1 * ddy2 + ddy1 * v2 ) - - -instance ( Module s ( T r ) , Num r ) => Num ( Dℝ3 r ) where - (+) = liftA2 (+) - (-) = liftA2 (-) - negate = fmap negate - fromInteger = konst @Double @( ℝ 3 ) . fromInteger - - abs = error "no" - signum = error "no" - - D3 v1 ( T dx1 ) ( T dy1 ) ( T dz1 ) ( T ddx1 ) ( T dxdy1 ) ( T ddy1 ) ( T dxdz1 ) ( T dydz1 ) ( T ddz1 ) - * D3 v2 ( T dx2 ) ( T dy2 ) ( T dz2 ) ( T ddx2 ) ( T dxdy2 ) ( T ddy2 ) ( T dxdz2 ) ( T dydz2 ) ( T ddz2 ) - = D3 ( v1 * v2 ) - ( T $ dx1 * v2 + v1 * dx2 ) - ( T $ dy1 * v2 + v1 * dy2 ) - ( T $ dz1 * v2 + v1 * dz2 ) - ( T $ dx1 * dx2 + ddx2 * v1 + ddx1 * v2 ) - ( T $ dy1 * dx2 + dx1 * dy2 + v1 * dxdy2 + dxdy1 * v2 ) - ( T $ dy1 * dy2 + v1 * ddy2 + ddy1 * v2 ) - ( T $ dz1 * dx2 + dx1 * dz2 + v1 * dxdz2 + dxdz1 * v2 ) - ( T $ dz1 * dy2 + dy1 * dz2 + v1 * dydz2 + dydz1 * v2 ) - ( T $ dz1 * dz2 + v1 * ddz2 + ddz1 * v2) - - -instance ( Num ( Dℝ0 r ), Module r ( T v ) ) => Module ( Dℝ0 r ) ( Dℝ0 v ) where - (^+^) = liftA2 ( coerce $ (^+^) @r @( T v ) ) - (^-^) = liftA2 ( coerce $ (^-^) @r @( T v ) ) - origin = pure ( coerce $ origin @r @( T v ) ) - (*^) = liftA2 ( coerce $ (*^) @r @( T v ) ) - -instance ( Num ( Dℝ1 r ), Module r ( T v ) ) => Module ( Dℝ1 r ) ( Dℝ1 v ) where - (^+^) = liftA2 ( coerce $ (^+^) @r @( T v ) ) - (^-^) = liftA2 ( coerce $ (^-^) @r @( T v ) ) - origin = pure ( coerce $ origin @r @( T v ) ) - (*^) = liftA2 ( coerce $ (*^) @r @( T v ) ) - -instance ( Num ( Dℝ2 r ), Module r ( T v ) ) => Module ( Dℝ2 r ) ( Dℝ2 v ) where - (^+^) = liftA2 ( coerce $ (^+^) @r @( T v ) ) - (^-^) = liftA2 ( coerce $ (^-^) @r @( T v ) ) - origin = pure ( coerce $ origin @r @( T v ) ) - (*^) = liftA2 ( coerce $ (*^) @r @( T v ) ) - -instance ( Num ( Dℝ3 r ), Module r ( T v ) ) => Module ( Dℝ3 r ) ( Dℝ3 v ) where - (^+^) = liftA2 ( coerce $ (^+^) @r @( T v ) ) - (^-^) = liftA2 ( coerce $ (^-^) @r @( T v ) ) - origin = pure ( coerce $ origin @r @( T v ) ) - (*^) = liftA2 ( coerce $ (*^) @r @( T v ) ) - -instance ( Module s ( T r ), Fractional r ) => Fractional ( Dℝ1 r ) where - (/) = error "I haven't yet defined (/) for Dℝ1" - fromRational = konst @Double @( ℝ 1 ) . fromRational -instance ( Module s ( T r ), Floating r ) => Floating ( Dℝ1 r ) where - pi = konst @Double @( ℝ 1 ) pi - sin ( D1 v ( T dx ) ( T ddx ) ) - = let !s = sin v - !c = cos v - in D1 s ( T $ c * dx ) ( T $ 2 * c * ddx - s * dx * dx ) - - cos ( D1 v ( T dx ) ( T ddx ) ) - = let !s = sin v - !c = cos v - in D1 c ( T $ -s * dx ) ( T $ -2 * s * ddx - c * dx * dx ) - -instance ( Module s ( T r ), Fractional r ) => Fractional ( Dℝ2 r ) where - (/) = error "I haven't yet defined (/) for Dℝ2" - fromRational = konst @Double @( ℝ 2 ) . fromRational -instance ( Module s ( T r ), Floating r ) => Floating ( Dℝ2 r ) where - pi = konst @Double @( ℝ 2 ) pi - sin ( D2 v ( T dx ) ( T dy ) ( T ddx ) ( T dxdy ) ( T ddy ) ) - = let !s = sin v - !c = cos v - in D2 s - ( T $ c * dx ) ( T $ c * dy ) - ( T $ 2 * c * ddx - s * dx * dx ) - ( T $ 2 * c * dxdy - 2 * s * dx * dy ) - ( T $ 2 * c * ddy - s * dy * dy ) - - cos ( D2 v ( T dx ) ( T dy ) ( T ddx ) ( T dxdy ) ( T ddy ) ) - = let !s = sin v - !c = cos v - in D2 c - ( T $ -s * dx ) ( T $ -s * dy ) - ( T $ -2 * s * ddx - c * dx * dx ) - ( T $ -2 * s * dxdy - 2 * c * dx * dy ) - ( T $ -2 * s * ddy - c * dy * dy ) - -instance ( Module s ( T r ), Fractional r ) => Fractional ( Dℝ3 r ) where - (/) = error "I haven't yet defined (/) for Dℝ3" - fromRational = konst @Double @( ℝ 3 ) . fromRational -instance ( Module s ( T r ), Floating r ) => Floating ( Dℝ3 r ) where - pi = konst @Double @( ℝ 3 ) pi - sin ( D3 v ( T dx ) ( T dy ) ( T dz ) ( T ddx ) ( T dxdy ) ( T ddy ) ( T dxdz ) ( T dydz ) ( T ddz ) ) - = let !s = sin v - !c = cos v - in D3 s - ( T $ c * dx ) ( T $ c * dy ) ( T $ c * dz ) - ( T $ 2 * c * ddx - s * dx * dx ) - ( T $ 2 * c * dxdy - 2 * s * dx * dy ) - ( T $ 2 * c * ddy - s * dy * dy ) - ( T $ 2 * c * dxdz - 2 * s * dx * dz ) - ( T $ 2 * c * dydz - 2 * s * dy * dz ) - ( T $ 2 * c * ddz - s * dz * dz ) - - cos ( D3 v ( T dx ) ( T dy ) ( T dz ) ( T ddx ) ( T dxdy ) ( T ddy ) ( T dxdz ) ( T dydz ) ( T ddz ) ) - = let !s = sin v - !c = cos v - in D3 c - ( T $ -s * dx ) ( T $ -s * dy ) ( T $ -s * dz ) - ( T $ -2 * s * ddx - c * dx * dx ) - ( T $ -2 * s * dxdy - 2 * c * dx * dy ) - ( T $ -2 * s * ddy - c * dy * dy ) - ( T $ -2 * s * dxdz - 2 * c * dx * dz ) - ( T $ -2 * s * dydz - 2 * c * dy * dz ) - ( T $ -2 * s * ddz - c * dz * dz ) - --------------------------------------------------------------------------------- - -uncurryD :: D a ~ D ( ℝ 1 ) - => D ( ℝ 1 ) ( a ~> b ) -> a -> D ( ℝ 2 ) b -uncurryD ( D1 ( D b_t0 ) ( T ( D dbdt_t0 ) ) ( T ( D d2bdt2_t0 ) ) ) s0 = - let !( D1 b_t0s0 dbds_t0s0 d2bds2_t0s0 ) = b_t0 s0 - !( D1 dbdt_t0s0 d2bdtds_t0s0 _ ) = dbdt_t0 s0 - !( D1 d2bdt2_t0s0 _ _ ) = d2bdt2_t0 s0 - in D2 b_t0s0 ( T dbdt_t0s0 ) dbds_t0s0 ( T d2bdt2_t0s0 ) d2bdtds_t0s0 d2bds2_t0s0 - --------------------------------------------------------------------------------- - -type Diffy :: Type -> Type -> Constraint -class ( Traversable ( D v ), Module r ( T v ) ) => Diffy r v where - chain :: ( Module r ( T w ) ) => D ( ℝ 1 ) v -> D v w -> D ( ℝ 1 ) w - konst :: Module s ( T w ) => w -> D v w - value :: D v w -> w - linear :: Module s ( T w ) => ( v -> w ) -> ( v ~> w ) - -chainRule :: ( Diffy r v, Module r ( T w ), D u ~ D ( ℝ 1 ) ) - => ( u ~> v ) -> ( v ~> w ) -> ( u ~> w ) -chainRule ( D df ) ( D dg ) = - D \ x -> - case df x of - df_x@( D1 { v = f_x } ) -> - chain df_x ( dg f_x ) - --- | Recover the underlying function, discarding all infinitesimal information. -fun :: forall v w r. Diffy r v => ( v ~> w ) -> ( v -> w ) -fun ( D df ) = value @r @v . df -{-# INLINE fun #-} - -var :: forall u r. ( Module r ( T r ), Representable r u, Diffy r u ) - => Fin ( Dim u ) -> ( u ~> r ) -var i = linear ( `index` i ) -{-# INLINE var #-} - -instance Diffy Double ( ℝ 0 ) where - chain _ ( D0 w ) = D1 w origin origin - {-# INLINE chain #-} - konst k = D0 k - {-# INLINE konst #-} - value ( D0 w ) = w - {-# INLINE value #-} - linear f = D \ _ -> D0 ( f ℝ0 ) - {-# INLINE linear #-} - -instance Diffy Double ( ℝ 1 ) where - chain ( D1 _ ( T ( ℝ1 x' ) ) ( T ( ℝ1 x'' ) ) ) ( D1 v g_x g_xx ) - = D1 v - ( x' *^ g_x ) - ( x'' *^ g_x ^+^ ( x' * x' ) *^ g_xx ) - {-# INLINE chain #-} - konst k = D1 k origin origin - {-# INLINE konst #-} - value ( D1 { v } ) = v - {-# INLINE value #-} - linear f = D \ u -> D1 ( f u ) ( T $ f u ) origin - {-# INLINE linear #-} - -instance Diffy Double ( ℝ 2 ) where - chain ( D1 _ ( T ( ℝ2 x' y' ) ) ( T ( ℝ2 x'' y'' ) ) ) ( D2 v g_x g_y g_xx g_xy g_yy ) - = D1 v - ( x' *^ g_x ^+^ y' *^ g_y ) - ( x'' *^ g_x ^+^ y'' *^ g_y - ^+^ ( x' * x' ) *^ g_xx ^+^ ( y' * y' ) *^ g_yy - ^+^ 2 *^ ( ( x' * y' ) *^ g_xy ) ) - {-# INLINE chain #-} - konst k = D2 k origin origin origin origin origin - {-# INLINE konst #-} - value ( D2 { v } ) = v - {-# INLINE value #-} - linear f = D \ u -> - D2 ( f u ) - ( T $ f ( ℝ2 1 0 ) ) ( T $ f ( ℝ2 0 1 ) ) - origin origin origin - {-# INLINE linear #-} - -instance Diffy Double ( ℝ 3 ) where - chain ( D1 _ ( T ( ℝ3 x' y' z' ) ) ( T ( ℝ3 x'' y'' z'' ) ) ) - ( D3 v g_x g_y g_z g_xx g_xy g_yy g_xz g_yz g_zz ) - = D1 v - ( x' *^ g_x ^+^ y' *^ g_y ^+^ z' *^ g_z ) - ( x'' *^ g_x ^+^ y'' *^ g_y ^+^ z'' *^ g_z - ^+^ ( x' * x' ) *^ g_xx ^+^ ( y' * y' ) *^ g_yy ^+^ ( z' * z' ) *^ g_zz - ^+^ 2 *^ ( ( x' * y' ) *^ g_xy ) ^+^ ( x' * z' ) *^ g_xz ^+^ ( y' * z' ) *^ g_yz ) - {-# INLINE chain #-} - konst k = D3 k origin origin origin origin origin origin origin origin origin - {-# INLINE konst #-} - value ( D3 { v } ) = v - {-# INLINE value #-} - linear f = D \ u -> - D3 ( f u ) - ( T $ f ( ℝ3 1 0 0 ) ) ( T $ f ( ℝ3 0 1 0 ) ) ( T $ f ( ℝ3 0 0 1 ) ) - origin origin origin origin origin origin - {-# INLINE linear #-} - --------------------------------------------------------------------------------- - --- TODO: avoid copying over code... - -instance Diffy ( 𝕀 Double ) ( 𝕀ℝ 0 ) where - chain _ ( D0 w ) = D1 w origin origin - {-# INLINE chain #-} - konst k = D0 k - {-# INLINE konst #-} - value ( D0 w ) = w - {-# INLINE value #-} - linear f = D \ _ -> D0 ( f $ singleton ℝ0 ) - {-# INLINE linear #-} - -instance Diffy ( 𝕀 Double ) ( 𝕀ℝ 1 ) where - chain ( D1 _ ( T ( I ( Rounded ( ℝ1 x'_lo ) ) ( Rounded ( ℝ1 x'_hi ) ) ) ) - ( T ( I ( Rounded ( ℝ1 x''_lo ) ) ( Rounded ( ℝ1 x''_hi ) ) ) ) ) - ( D1 v g_x g_xx ) - = let - !x' = I ( Rounded x'_lo ) ( Rounded x'_hi ) - !x'' = I ( Rounded x''_lo ) ( Rounded x''_hi ) - in D1 v - ( x' *^ g_x ) - ( x'' *^ g_x ^+^ ( x' * x' ) *^ g_xx ) - {-# INLINE chain #-} - konst k = D1 k origin origin - {-# INLINE konst #-} - value ( D1 { v } ) = v - {-# INLINE value #-} - linear f = D \ u -> D1 ( f u ) ( T $ f u ) origin - {-# INLINE linear #-} - -instance Diffy ( 𝕀 Double ) ( 𝕀ℝ 2 ) where - chain ( D1 _ ( T ( I ( Rounded ( ℝ2 x'_lo y'_lo ) ) ( Rounded ( ℝ2 x'_hi y'_hi ) ) ) ) - ( T ( I ( Rounded ( ℝ2 x''_lo y''_lo ) ) ( Rounded ( ℝ2 x''_hi y''_hi ) ) ) ) ) - ( D2 v g_x g_y g_xx g_xy g_yy ) - = let - !x' = I ( Rounded x'_lo ) ( Rounded x'_hi ) - !y' = I ( Rounded y'_lo ) ( Rounded y'_hi ) - !x'' = I ( Rounded x''_lo ) ( Rounded x''_hi ) - !y'' = I ( Rounded y''_lo ) ( Rounded y''_hi ) - in D1 v - ( x' *^ g_x ^+^ y' *^ g_y ) - ( x'' *^ g_x ^+^ y'' *^ g_y - ^+^ ( x' * x' ) *^ g_xx ^+^ ( y' * y' ) *^ g_yy - ^+^ 2 *^ ( ( x' * y' ) *^ g_xy ) ) - {-# INLINE chain #-} - konst k = D2 k origin origin origin origin origin - {-# INLINE konst #-} - value ( D2 { v } ) = v - {-# INLINE value #-} - linear f = D \ u -> - D2 ( f u ) - ( T $ f ( singleton $ ℝ2 1 0 ) ) ( T $ f ( singleton $ ℝ2 0 1 ) ) - origin origin origin - {-# INLINE linear #-} - -instance Diffy ( 𝕀 Double ) ( 𝕀ℝ 3 ) where - chain ( D1 _ ( T ( I ( Rounded ( ℝ3 x'_lo y'_lo z'_lo ) ) ( Rounded ( ℝ3 x'_hi y'_hi z'_hi ) ) ) ) - ( T ( I ( Rounded ( ℝ3 x''_lo y''_lo z''_lo ) ) ( Rounded ( ℝ3 x''_hi y''_hi z''_hi ) ) ) ) ) - ( D3 v g_x g_y g_z g_xx g_xy g_yy g_xz g_yz g_zz ) - = let - !x' = I ( Rounded x'_lo ) ( Rounded x'_hi ) - !y' = I ( Rounded y'_lo ) ( Rounded y'_hi ) - !z' = I ( Rounded z'_lo ) ( Rounded z'_hi ) - !x'' = I ( Rounded x''_lo ) ( Rounded x''_hi ) - !y'' = I ( Rounded y''_lo ) ( Rounded y''_hi ) - !z'' = I ( Rounded z''_lo ) ( Rounded z''_hi ) - in D1 v - ( x' *^ g_x ^+^ y' *^ g_y ^+^ z' *^ g_z ) - ( x'' *^ g_x ^+^ y'' *^ g_y ^+^ z'' *^ g_z - ^+^ ( x' * x' ) *^ g_xx ^+^ ( y' * y' ) *^ g_yy ^+^ ( z' * z' ) *^ g_zz - ^+^ 2 *^ ( ( x' * y' ) *^ g_xy ) ^+^ ( x' * z' ) *^ g_xz ^+^ ( y' * z' ) *^ g_yz ) - {-# INLINE chain #-} - konst k = D3 k origin origin origin origin origin origin origin origin origin - {-# INLINE konst #-} - value ( D3 { v } ) = v - {-# INLINE value #-} - linear f = D \ u -> - D3 ( f u ) - ( T $ f ( singleton $ ℝ3 1 0 0 ) ) ( T $ f ( singleton $ ℝ3 0 1 0 ) ) ( T $ f ( singleton $ ℝ3 0 0 1 ) ) - origin origin origin origin origin origin - {-# INLINE linear #-} +-- | \( \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 - ( Diffy ( I i Double ) ( I i u ) - , Module ( I i Double ) ( T ( I i Double ) ) + ( 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 ) @@ -420,11 +259,605 @@ class ) => Differentiable i u instance - ( Diffy ( I i Double ) ( I i u ) - , Module ( I i Double ) ( T ( I i Double ) ) + ( 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 #-}