diff --git a/brush-strokes/brush-strokes.cabal b/brush-strokes/brush-strokes.cabal
index 9ff3b55..08bf574 100644
--- a/brush-strokes/brush-strokes.cabal
+++ b/brush-strokes/brush-strokes.cabal
@@ -205,10 +205,38 @@ executable inspect
         Math.Interval.Abstract
 
     build-depends:
-      brush-strokes,
-      data-reify
+        brush-strokes
+      , data-reify
           ^>= 0.6.3
 
+test-suite tests
+
+    import:
+      common
+
+    type:
+      exitcode-stdio-1.0
+
+    hs-source-dirs:
+      src/test
+
+    default-language:
+      Haskell2010
+
+    main-is:
+      Main.hs
+
+    build-depends:
+        brush-strokes
+      , falsify
+          ^>= 0.2
+      , hspray
+          ^>= 0.1.3
+      , tasty
+          ^>= 1.5
+      , unordered-containers
+          >= 0.2.15 && < 0.3
+
 benchmark cusps
 
   import:
diff --git a/brush-strokes/cabal.project b/brush-strokes/cabal.project
index 51f427c..22cd7a7 100644
--- a/brush-strokes/cabal.project
+++ b/brush-strokes/cabal.project
@@ -4,10 +4,13 @@ constraints:
   acts -finitary,
   rounded-hw -pure-hs -c99 -avx512 +ghc-prim -x87-long-double
 
+tests: True
+
 allow-newer:
   acts:base, acts:deepseq,
   groups-generic:base,
   eigen:primitive,
+  falsify:base, falsify:tasty,
 
 --package brush-strokes
 profiling: True
diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs
index fc804e0..c58128f 100644
--- a/brush-strokes/src/cusps/bench/Main.hs
+++ b/brush-strokes/src/cusps/bench/Main.hs
@@ -608,3 +608,34 @@ x2  = 𝕀  x2_lo  x2_hi
 
 ( b1 - a12 * x2 ) `extendedDivide` a11
 -}
+
+{-
+SMOKING GUN:
+
+> (f, fI) = brushStrokeFunctions $ ellipseBrushStroke (0,1) pi
+> (t, i, s) = mkBox (0.54, 0.55) 2 (0.5480, 0.5481)
+> _D32_dxdxdx $ dbrush $ eval fI (t,i,s)
+> [ℝ2 54187.61174031432 -11626.47655724034, ℝ2 55257.54384135226 -9759.09092706883]
+
+In Mathematica:
+
+  FindMinimum[{D[b[vt, s, myPb0, myPb1][[1]], {vt, 3}] /. (vt -> t), 
+    0.54 <= t <= 0.55 && 0.548 <= s <= 0.5481}, {{t, 0.545}, {s, 
+     0.54805}}]
+> { -503.445, {t -> 0.54, s -> 0.548} }
+
+  FindMaximum[{D[b[vt, s, myPb0, myPb1][[1]], {vt, 3}] /. (vt -> t), 
+    0.54 <= t <= 0.55 && 0.548 <= s <= 0.5481}, {{t, 0.545}, {s, 
+     0.54805}}]
+> { -480.086, {t -> 0.55, s -> 0.5481} }
+
+so the actual range is [-503.445,-480.086] but we have computed a range of
+[54187.61174031432, 55257.54384135226] which does not contain the actual range!
+
+-}
+
+{-
+
+Math.Ring.cos  ( linearD @_ @3 (\i -> 𝕀 (getR1 $ inf i) (getR1 $ sup i)) (𝕀 (ℝ1 $ 0.548 * pi) (ℝ1 $ 0.5481 * pi)) :: D3𝔸1 ( 𝕀 Double ) )
+
+-}
\ No newline at end of file
diff --git a/brush-strokes/src/cusps/inspect/Main.hs b/brush-strokes/src/cusps/inspect/Main.hs
index a876c12..a412124 100644
--- a/brush-strokes/src/cusps/inspect/Main.hs
+++ b/brush-strokes/src/cusps/inspect/Main.hs
@@ -140,8 +140,8 @@ boo x i y z
 ellipseTest :: StrokeDatum 3 AI
 ellipseTest =
   boo ( ellipseBrush ( Pt . Val ) scaleI ) 1
-    ( mkPoint ( ℝ2 0 0 ) 1.0111 1.0222 0 )
-    ( LineTo { curveEnd = NextPoint $ mkPoint ( ℝ2 100 0 ) 10.0333 10.0444 ( pi Prelude./ 2 )
+    ( mkPoint ( ℝ2 0 0 ) 10 25 0 )
+    ( LineTo { curveEnd = NextPoint $ mkPoint ( ℝ2 100 0 ) 15 40 pi
              , curveData = () } )
   where
     mkPoint :: ℝ 2 -> Double -> Double -> Double -> PointData ( ℝ 3 )
@@ -258,7 +258,7 @@ instance HasBézier 3 AI where
           ( Cubic.bezier''' bez )
 
 --------------------------------------------------------------------------------
--- Brushes (TODO copied from MetaBrush.Asset.Brushes)
+-- Brushes (TODO copied from Calligraphy.Brushes)
 
 κ :: Double
 κ = 0.5519150244935105707435627227925
@@ -304,13 +304,13 @@ ellipseBrush mkI scaleI =
         mkPt x y
           = let !x' = a `scaledBy` x
                 !y' = b `scaledBy` y
-{-
+-- {-
                 !c  = cos phi
                 !s  = sin phi
             in   ( x' * c - y' * s ) *^ e_x
              ^+^ ( x' * s + y' * c ) *^ e_y
--}
--- {-
+-- -}
+{-
                 !r  = sqrt $ x' ^ 2 + y' ^ 2
                 !arctgt = atan ( y' / x' )
                 -- a and b are always strictly positive, so we can compute
@@ -326,7 +326,7 @@ ellipseBrush mkI scaleI =
             in
                 ( r * cos phi' ) *^ e_x
             ^+^ ( r * sin phi' ) *^ e_y
--- -}
+-}
 
     in circleSpline @AI @k @( ℝ 3 ) @( ℝ 2 ) mkPt
   where
diff --git a/brush-strokes/src/lib/Math/Algebra/Dual/Internal.hs b/brush-strokes/src/lib/Math/Algebra/Dual/Internal.hs
index 181ec9f..2d9c446 100644
--- a/brush-strokes/src/lib/Math/Algebra/Dual/Internal.hs
+++ b/brush-strokes/src/lib/Math/Algebra/Dual/Internal.hs
@@ -244,7 +244,7 @@ chainRule1NQ zero_w sum_w scale_w df dg =
                   ]
               ) ||]
 
--- | The chain rule for a composition \( \mathbb{R}^n \to \mathbb{R} \to \mathbb{R} \)
+-- | The chain rule for a composition \( \mathbb{R}^n \to \mathbb{R}^1 \to \mathbb{R}^1 \)
 --
 -- (To be spliced in using Template Haskell.)
 chainRuleN1Q :: forall du dr1 r
@@ -291,7 +291,7 @@ type instance Vars D𝔸0 = 0
 instance MonomialBasis D𝔸0 where
   monTabulate f =
     [|| let
-          !_D0_v = $$( f $ Mon VZ )
+          _D0_v = $$( f $ Mon VZ )
         in D0 { .. }
       ||]
 
@@ -302,8 +302,8 @@ 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 ) )
+          _D11_v    =   $$( f $ Mon ( 0 `VS` VZ ) )
+          _D11_dx   = T $$( f $ Mon ( 1 `VS` VZ ) )
         in D11 { .. }
       ||]
 
@@ -316,9 +316,9 @@ 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 ) )
+          _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 { .. }
       ||]
 
@@ -332,10 +332,10 @@ 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 ) )
+          _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 { .. }
       ||]
 
@@ -350,9 +350,9 @@ 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 ) )
+          _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 { .. }
       ||]
 
@@ -366,12 +366,12 @@ 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 ) )
+          _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 { .. }
       ||]
 
@@ -388,16 +388,16 @@ 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 ) )
+          _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
@@ -435,16 +435,16 @@ 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 ) )
+          _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 { .. }
       ||]
 
@@ -466,26 +466,26 @@ 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 ) )
+          _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
@@ -505,6 +505,7 @@ instance MonomialBasis D3𝔸3 where
     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` 2 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dydydz $$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 ||]
@@ -514,11 +515,11 @@ 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 ) )
+          _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
@@ -533,21 +534,21 @@ 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 ) )
+          _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
@@ -573,41 +574,41 @@ 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 ) )
+          _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
diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs b/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs
index fcbb7d6..b05f051 100644
--- a/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs
+++ b/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs
@@ -180,10 +180,10 @@ instance HasEnvelopeEquation 2 where
        -- ∂²c/∂s²  = ∂²b/∂s²
 
   envelopeEquation co ( D22 _ c_t c_s c_tt c_ts c_ss ) =
-    let !ee   = c_t × c_s
-        !ee_t = c_tt × c_s + c_t × c_ts
-        !ee_s = c_ts × c_s + c_t × c_ss
-        !𝛿E𝛿sdcdt = ee_s *^ c_t ^-^ ee_t *^ c_s
+    let ee   = c_t × c_s
+        ee_t = c_tt × c_s + c_t × c_ts
+        ee_s = c_ts × c_s + c_t × c_ss
+        𝛿E𝛿sdcdt = ee_s *^ c_t ^-^ ee_t *^ c_s
           -- TODO: we get c_t * c_t and c_s * c_s terms...
           -- These could be squares (better with interval arithmetic)?
     in ( D12 ( co ee ) ( T $ co ee_t ) ( T $ co ee_s )
@@ -235,24 +235,24 @@ instance HasEnvelopeEquation 3 where
     ( D32 _ c_t c_s
           c_tt c_ts c_ss
           c_ttt c_tts c_tss c_sss )
-    = let !ee    = c_t × c_s
-          !ee_t  = c_tt × c_s + c_t × c_ts
-          !ee_s  = c_ts × c_s + c_t × c_ss
-          !ee_tt = c_ttt × c_s
-                 + c_tt  × c_ts * 2
-                 + c_t   × c_tts
-          !ee_ts = c_tts × c_s
-                 + c_tt  × c_ss
-              -- + c_ts  × c_ts -- cancels out
-                 + c_t   × c_tss
-          !ee_ss = c_tss × c_s
-                 + c_ts  × c_ss * 2
-                 + c_t   × c_sss
-          !𝛿E𝛿sdcdt   = ee_s *^ c_t ^-^ ee_t *^ c_s
-          !𝛿E𝛿sdcdt_t =    ee_ts *^ c_t ^+^ ee_s *^ c_tt
-                     ^-^ ( ee_tt *^ c_s ^+^ ee_t *^ c_ts )
-          !𝛿E𝛿sdcdt_s =   ee_ss *^ c_t ^+^ ee_s *^ c_ts
-                    ^-^ ( ee_ts *^ c_s ^+^ ee_t *^ c_ss )
+    = let ee    = c_t × c_s
+          ee_t  = c_tt × c_s + c_t × c_ts
+          ee_s  = c_ts × c_s + c_t × c_ss
+          ee_tt = c_ttt × c_s
+                + c_tt  × c_ts * 2
+                + c_t   × c_tts
+          ee_ts = c_tts × c_s
+                + c_tt  × c_ss
+             -- + c_ts  × c_ts -- cancels out
+                + c_t   × c_tss
+          ee_ss = c_tss × c_s
+                + c_ts  × c_ss * 2
+                + c_t   × c_sss
+          𝛿E𝛿sdcdt   = ee_s *^ c_t ^-^ ee_t *^ c_s
+          𝛿E𝛿sdcdt_t =    ee_ts *^ c_t ^+^ ee_s *^ c_tt
+                    ^-^ ( ee_tt *^ c_s ^+^ ee_t *^ c_ts )
+          𝛿E𝛿sdcdt_s =   ee_ss *^ c_t ^+^ ee_s *^ c_ts
+                   ^-^ ( ee_ts *^ c_s ^+^ ee_t *^ c_ss )
       in ( D22
             ( co ee )
             ( T $ co ee_t ) ( T $ co ee_s )
diff --git a/brush-strokes/src/lib/Math/Linear.hs b/brush-strokes/src/lib/Math/Linear.hs
index 77344df..e0260a2 100644
--- a/brush-strokes/src/lib/Math/Linear.hs
+++ b/brush-strokes/src/lib/Math/Linear.hs
@@ -122,7 +122,9 @@ data Vec n a where
     -- can't be strict, otherwise we can't conveniently
     -- unsafeCoerce from lists
 
-deriving stock instance Show a => Show ( Vec n a )
+--deriving stock instance Show a => Show ( Vec n a )
+instance Show a => Show ( Vec n a ) where
+  showsPrec p v = showsPrec p ( unsafeCoerce v :: [ a ] )
 
 deriving stock instance Functor     ( Vec n )
 deriving stock instance Foldable    ( Vec n )
@@ -137,9 +139,7 @@ instance Ord a => Ord ( Vec n a ) where
 
 infixl 9 !
 (!) :: forall l a. Vec l a -> Fin l -> a
-VS a _ ! Fin 1 = a
-VS _ a ! Fin i = a ! Fin ( i - 1 )
-_      ! _     = error "impossible: Fin 0 is uninhabited"
+v ! Fin i = ( unsafeCoerce v :: [ a ] ) !! fromIntegral i
 
 find :: forall l a. ( a -> Bool ) -> Vec l a -> MFin l
 find f v = MFin ( find_ 1 v )
diff --git a/brush-strokes/src/test/Main.hs b/brush-strokes/src/test/Main.hs
new file mode 100644
index 0000000..b85fc22
--- /dev/null
+++ b/brush-strokes/src/test/Main.hs
@@ -0,0 +1,319 @@
+{-# LANGUAGE AllowAmbiguousTypes #-}
+{-# LANGUAGE NumericUnderscores  #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TemplateHaskell #-}
+
+module Main (main) where
+
+-- base
+import Prelude hiding
+  ( Num(..), (^) )
+import Data.Foldable
+  ( toList )
+import Data.List.NonEmpty
+  ( NonEmpty(..) )
+import Data.Maybe
+  ( catMaybes )
+import Data.Traversable
+  ( for )
+import Unsafe.Coerce
+  ( unsafeCoerce )
+
+-- brush-strokes
+import Math.Algebra.Dual
+import Math.Linear
+import Math.Module
+import Math.Monomial
+  ( multiSubsetSum, multiSubsetsSum
+  , MonomialBasis ( monTabulate, monIndex )
+  )
+import Math.Ring
+
+-- hspray
+import Math.Algebra.Hspray
+  ( Spray )
+import qualified Math.Algebra.Hspray as Spray
+
+-- falsify
+import Test.Tasty.Falsify
+import qualified Test.Falsify.Generator as Falsify
+  ( Gen )
+import qualified Test.Falsify.Generator as Falsify.Gen
+import Test.Falsify.Predicate
+  ( (.$) )
+import qualified Test.Falsify.Predicate as Falsify.Prop
+import qualified Test.Falsify.Property as Falsify
+  ( Property
+  , assert
+  , discard
+  , gen, genWith
+  )
+import qualified Test.Falsify.Range as Falsify
+
+-- tasty
+import qualified Test.Tasty as Tasty
+
+-- unordered-containers
+import qualified Data.HashMap.Lazy as HashMap
+
+--------------------------------------------------------------------------------
+
+
+main :: IO ()
+main =
+  Tasty.defaultMain $
+    Tasty.testGroup "brush-strokes property tests"
+    [ Tasty.testGroup "Automatic differentiation"
+      [ Tasty.testGroup "Monomial basis"
+        [ testProperty "Round trip D33" testMonomialBasisD33
+        ]
+      , Tasty.testGroup "Monomials"
+        [ Tasty.testGroup "multiSubsetSum"
+          [ testProperty "multiSubsetSum valid" testMultiSubsetSumValid
+          , testProperty "multiSubsetSum exhaustive" testMultiSubsetSumExhaustive
+          ]
+        -- , Tasty.testGroup "multiSubsetsSum"
+        --   [ testProperty "multiSubsetsSum exhaustive" testMultiSubsetsSumExhaustive
+        --   ]
+        ]
+      , Tasty.testGroup "chainRule1NQ"
+        [ testProperty "chainRule1NQ_1" testChainRule1NQ_1
+        , testProperty "chainRule1NQ_2" testChainRule1NQ_2
+        , testProperty "chainRule1NQ_3" testChainRule1NQ_3
+        ]
+      ]
+    ]
+
+-- | Check that the 'multiSubsetSum' function returns valid answers, i.e.
+-- all returned multisubsets have the desired size and sum.
+testMultiSubsetSumValid :: Falsify.Property ()
+testMultiSubsetSumValid = do
+  rg  <- Falsify.genWith (\ rg  -> Just $ "range = " ++ show rg ) $ Falsify.Gen.inRange $ Falsify.between ( 1, 6 )
+  sz  <- Falsify.genWith (\ sz  -> Just $  "size = " ++ show sz ) $ Falsify.Gen.inRange $ Falsify.between ( 0, 20 )
+  tot <- Falsify.genWith (\ tot -> Just $   "tot = " ++ show tot) $ Falsify.Gen.inRange $ Falsify.between ( sz, sz * rg )
+  let range = [ 1 .. rg ]
+      mss = multiSubsetSum sz tot range
+  case mss of
+    [] -> Falsify.discard
+    r:rs -> do
+      ms <- Falsify.gen $ Falsify.Gen.elem ( r :| rs )
+      Falsify.assert
+         $ Falsify.Prop.eq
+        .$ ("(sz, tot)", (sz, tot) )
+        .$ ("computed (sz, tot)", (size ms, total ms))
+  where
+    size, total :: [ ( Word, Word ) ] -> Word
+    size [] = 0
+    size ((_,n):ins) = n + size ins
+    total [] = 0
+    total ((i,n):ins) = i * n + total ins
+
+-- | Check that the 'multiSubsetSum' function returns all multisubsets of
+-- the given set, by generating a random multisubset, computing its size, and
+-- checking it belongs to the output of the 'multiSubsetSum' function.
+testMultiSubsetSumExhaustive :: Falsify.Property ()
+testMultiSubsetSumExhaustive = do
+  rg <- Falsify.genWith (\ rg -> Just $ "range = " ++ show rg) $ Falsify.Gen.inRange $ Falsify.between ( 1, 6 )
+  sz <- Falsify.genWith (\ sz -> Just $  "size = " ++ show sz) $ Falsify.Gen.inRange $ Falsify.between ( 0, 10 )
+  let range = [ 1 .. rg ]
+  (multiSubset, tot) <- Falsify.genWith (\ ms -> Just $ "multisubset = " ++ show ms) $ genMultiSubset range sz
+  Falsify.assert
+     $ Falsify.Prop.elem
+    .$ ("all multisubsets", multiSubsetSum sz tot range )
+    .$ ("random multisubset", multiSubset)
+
+genMultiSubset :: [ Word ] -> Word -> Falsify.Gen ( [ ( Word, Word ) ] , Word )
+genMultiSubset [i] sz =
+  return $
+    if sz == 0
+    then ( [], 0 )
+    else ( [ ( i, sz ) ], i * sz )
+genMultiSubset (i:is) sz = do
+  nb <- Falsify.Gen.inRange $ Falsify.between ( 0, sz )
+  (rest, tot) <- genMultiSubset is ( sz - nb )
+  return $ ( if nb == 0 then rest else ( i, nb ) : rest, tot + nb * i )
+genMultiSubset [] _ = error "impossible"
+
+coerceVec1 :: [ a ] -> Vec n a
+coerceVec1 = unsafeCoerce
+
+coerceVec2 :: Vec n a -> [ a ]
+coerceVec2 = toList
+
+-- | Check that the 'multiSubsetSums' function returns all collections of
+-- multisubsets of the given set (see 'testMultiSubsetSumExhaustive').
+testMultiSubsetsSumExhaustive :: Falsify.Property ()
+testMultiSubsetsSumExhaustive = do
+  rg <- Falsify.genWith (\ rg -> Just $ "range = " ++ show rg) $ Falsify.Gen.inRange $ Falsify.between ( 1, 5 )
+  let range = [ 1 .. rg ]
+  n  <- Falsify.genWith (\ n  -> Just $     "n = " ++ show n ) $ Falsify.Gen.inRange $ Falsify.between ( 1, 10 )
+  multiSubsets <- for ( [ 0 .. n - 1 ] :: [ Word ] ) \ i -> do
+    sz <- Falsify.gen $ Falsify.Gen.inRange $ Falsify.between ( 0, 5 )
+    ( ms, tot ) <- Falsify.genWith ( \ ms -> Just $ "ms_" ++ show i ++ " = " ++ show ms ) $ genMultiSubset range sz
+    return ( ms, sz, tot )
+  let mss = map ( \ (ms, _,_) -> ms ) multiSubsets
+      szs = map ( \ (_,sz,_) -> sz) multiSubsets
+      tot = sum $ map ( \(_,_,t) -> t) multiSubsets
+  Falsify.assert
+     $ Falsify.Prop.elem
+    .$ ("all multisubsets", map coerceVec2 $ multiSubsetsSum range tot $ coerceVec1 szs )
+    .$ ("random multisubset", mss)
+
+
+testRoundTrip
+  :: ( Show a, Eq a )
+  => Falsify.Gen a
+  -> ( a -> a )
+  -> Falsify.Property ()
+testRoundTrip g roundTrip = do
+  d <- Falsify.gen g
+  Falsify.assert
+     $ Falsify.Prop.eq
+    .$ ("value", d )
+    .$ ("round tripped", roundTrip d )
+
+testMonomialBasisD33 :: Falsify.Property ()
+testMonomialBasisD33 =
+  testRoundTrip genD33 \ d -> $$( monTabulate \ mon -> monIndex [|| d ||] mon )
+  where
+    genD33 :: Falsify.Gen ( D3𝔸3 Double )
+    genD33 =
+      D33 <$> (unT <$> g)
+          <*> g <*> g <*> g
+          <*> g <*> g <*> g <*> g <*> g <*> g
+          <*> g <*> g <*> g <*> g <*> g <*> g <*> g <*> g <*> g <*> g
+
+    g :: Falsify.Gen ( T Double )
+    g = T . fromIntegral <$> Falsify.Gen.inRange ( Falsify.withOrigin ( -100, 100 ) ( 0 :: Int ) )
+
+-- | Test the Faà di Bruno formula on polynomials, with a composition
+-- \( g(f_1(x), f_2(x), .., f_n(x)) \).
+testChainRule1NQ_1 :: Falsify.Property ()
+testChainRule1NQ_1 = do
+  f <- genSpray "f" 1
+  g <- genSpray "g" 1
+  let gof_spray = Spray.composeSpray g [f]
+      gof_chain =
+        chain @_ @3 @( ℝ 1 ) ( ℝ1 <$> fromSpray @3 @( ℝ 1 ) f ) ( fromSpray @3 @( ℝ 1 ) g )
+  Falsify.assert
+     $ Falsify.Prop.eq
+    .$ ("direct", fromSpray @3 @( ℝ 1 ) gof_spray )
+    .$ ("chain rule", gof_chain )
+
+-- | Test the Faà di Bruno formula on polynomials, with a composition
+-- \( g(f_1(x), f_2(x), .., f_n(x)) \).
+testChainRule1NQ_2 :: Falsify.Property ()
+testChainRule1NQ_2 = do
+  f1 <- genSpray "f1" 1
+  f2 <- genSpray "f2" 1
+  g  <- genSpray "g"  2
+  let gof_spray = Spray.composeSpray g [f1, f2]
+      f = ℝ2 <$> fromSpray @3 @( ℝ 1 ) f1
+             <*> fromSpray @3 @( ℝ 1 ) f2
+      gof_chain =
+        chain @_ @3 @( ℝ 2 ) f ( fromSpray @3 @( ℝ 2 ) g )
+  Falsify.assert
+     $ Falsify.Prop.eq
+    .$ ("direct", fromSpray @3 @( ℝ 1 ) gof_spray )
+    .$ ("chain rule", gof_chain )
+
+-- | Test the Faà di Bruno formula on polynomials, with a composition
+-- \( g(f_1(x), f_2(x), .., f_n(x)) \).
+testChainRule1NQ_3 :: Falsify.Property ()
+testChainRule1NQ_3 = do
+  f1 <- genSpray "f1" 1
+  f2 <- genSpray "f2" 1
+  f3 <- genSpray "f3" 1
+  g  <- genSpray "g"  3
+  let gof_spray = Spray.composeSpray g [f1, f2, f3]
+      f = ℝ3 <$> fromSpray @3 @( ℝ 1 ) f1
+             <*> fromSpray @3 @( ℝ 1 ) f2
+             <*> fromSpray @3 @( ℝ 1 ) f3
+      gof_chain =
+        chain @_ @3 @( ℝ 3 ) f ( fromSpray @3 @( ℝ 3 ) g )
+  Falsify.assert
+     $ Falsify.Prop.eq
+    .$ ("direct", fromSpray @3 @( ℝ 1 ) gof_spray )
+    .$ ("chain rule", gof_chain )
+
+class FromSpray v where
+  varFn :: Int -> v
+  linFn :: v -> Int -> Double
+
+instance FromSpray ( ℝ 1 ) where
+  varFn = \case
+    0 -> ℝ1 1
+    i -> error $ "fromSpray in 1d but variable " ++ show i
+  linFn ( ℝ1 x ) = \case
+    0 -> x
+    i -> error $ "fromSpray in 1d but variable " ++ show i
+instance FromSpray ( ℝ 2 ) where
+  varFn = \case
+    0 -> ℝ2 1 0
+    1 -> ℝ2 0 1
+    i -> error $ "fromSpray in 2d but variable " ++ show i
+  linFn ( ℝ2 x y ) = \case
+    0 -> x
+    1 -> y
+    i -> error $ "fromSpray in 2d but variable " ++ show i
+instance FromSpray ( ℝ 3 ) where
+  varFn = \case
+    0 -> ℝ3 1 0 0
+    1 -> ℝ3 0 1 0
+    2 -> ℝ3 0 0 1
+    i -> error $ "fromSpray in 3d but variable " ++ show i
+  linFn ( ℝ3 x y z ) = \case
+    0 -> x
+    1 -> y
+    2 -> z
+    i -> error $ "fromSpray in 3d but variable " ++ show i
+
+genSpray :: String -> Word -> Falsify.Property ( Spray Double )
+genSpray lbl nbVars = Falsify.genWith (\ p -> Just $ lbl ++ " = " ++ Spray.prettySpray show "x" p) $ do
+  deg <- Falsify.Gen.inRange $ Falsify.between ( 0, 10 )
+  let mons = allMonomials deg nbVars
+  coeffs <-
+    for mons $ \ mon -> do
+      if all (== 0) mon
+      then return Nothing
+      else do
+        nonZero <- Falsify.Gen.bool False
+        if nonZero
+        then return Nothing
+        else do
+          -- Just use (small) integral values in tests for now,
+          -- to avoid errors arising from rounding.
+          c <- Falsify.Gen.inRange $ Falsify.withOrigin ( -100, 100 ) ( 0 :: Int )
+          return $ Just ( map fromIntegral mon, fromIntegral c )
+  return $ Spray.fromList $ catMaybes coeffs
+
+allMonomials :: Word -> Word -> [ [ Word ] ]
+allMonomials k _ | k < 0 = []
+allMonomials _ 0 = [ [] ]
+allMonomials 0 n = [ replicate ( fromIntegral n ) 0 ]
+allMonomials k n = [ i : is | i <- reverse [ 0 .. k ], is <- allMonomials ( k - i ) ( n - 1 ) ]
+
+-- | Convert a multivariate polynomial from the @hspray@ library to the dual algebra.
+fromSpray
+  :: forall k v
+  .  ( HasChainRule Double k v
+     , Module Double (T v)
+     , Applicative ( D k v )
+     , Ring ( D k v Double )
+     , FromSpray v
+     )
+  => Spray Double
+  -> D k v Double
+fromSpray coeffs = HashMap.foldlWithKey' addMonomial ( konst @Double @k @v $ HashMap.lookupDefault 0 (Spray.Powers mempty 0) coeffs ) coeffs
+  where
+    addMonomial :: D k v Double -> Spray.Powers -> Double -> D k v Double
+    addMonomial a xs c = a + monomial c ( toList $ Spray.exponents xs )
+    monomial :: Double -> [ Int ] -> D k v Double
+    monomial _ [] = konst @Double @k @v 0
+    monomial c is = fmap ( c * ) $ go 0 is
+    go :: Int -> [ Int ] -> D k v Double
+    go _ [] = konst @Double @k @v 1
+    go d (i : is) = pow d i * go ( d + 1 ) is
+    pow :: Int -> Int -> D k v Double
+    pow _ 0 = konst @Double @k @v 1
+    pow d i = linearD @Double @k @v ( \ x -> linFn @v x d ) ( unT origin :: v ) ^ ( fromIntegral i )