diff -r bd13c2ae3450 -r 56c8adc32b09 src/kernels/gaussian.rs
--- a/src/kernels/gaussian.rs	Fri Apr 28 13:15:19 2023 +0300
+++ b/src/kernels/gaussian.rs	Tue Dec 31 09:34:24 2024 -0500
@@ -17,9 +17,10 @@
     Weighted,
     Bounded,
 };
-use alg_tools::mapping::Apply;
+use alg_tools::mapping::{Apply, Differentiable};
 use alg_tools::maputil::array_init;
 
+use crate::types::Lipschitz;
 use crate::fourier::Fourier;
 use super::base::*;
 use super::ball_indicator::CubeIndicator;
@@ -75,14 +76,45 @@
 impl<S, const N : usize> Apply<Loc<S::Type, N>> for Gaussian<S, N>
 where S : Constant {
     type Output = S::Type;
-    // This is not normalised to neither to have value 1 at zero or integral 1
-    // (unless the cut-off ε=0).
     #[inline]
     fn apply(&self, x : Loc<S::Type, N>) -> Self::Output {
         self.apply(&x)
     }
 }
 
+#[replace_float_literals(S::Type::cast_from(literal))]
+impl<'a, S, const N : usize> Differentiable<&'a Loc<S::Type, N>> for Gaussian<S, N>
+where S : Constant {
+    type Output = Loc<S::Type, N>;
+    #[inline]
+    fn differential(&self, x : &'a Loc<S::Type, N>) -> Self::Output {
+        x * (self.apply(x) / self.variance.value())
+    }
+}
+
+impl<S, const N : usize> Differentiable<Loc<S::Type, N>> for Gaussian<S, N>
+where S : Constant {
+    type Output = Loc<S::Type, N>;
+    // This is not normalised to neither to have value 1 at zero or integral 1
+    // (unless the cut-off ε=0).
+    #[inline]
+    fn differential(&self, x : Loc<S::Type, N>) -> Self::Output {
+        x * (self.apply(&x) / self.variance.value())
+    }
+}
+
+#[replace_float_literals(S::Type::cast_from(literal))]
+impl<S, const N : usize> Lipschitz<L2> for Gaussian<S, N>
+where S : Constant {
+    type FloatType = S::Type;
+    fn lipschitz_factor(&self, L2 : L2) -> Option<Self::FloatType> {
+        // f(x)=f_1(‖x‖_2/σ) * √(2π) / √(2πσ)^N, where f_1 is one-dimensional Gaussian with
+        // variance 1. The Lipschitz factor of f_1 is e^{-1/2}/√(2π), see, e.g.,
+        // https://math.stackexchange.com/questions/3630967/is-the-gaussian-density-lipschitz-continuous
+        // Thus the Lipschitz factor we want is e^{-1/2} / (√(2πσ)^N * σ).
+        Some((-0.5).exp() / (self.scale() * self.variance.value().sqrt()))
+    }
+}
 
 #[replace_float_literals(S::Type::cast_from(literal))]
 impl<'a, S, const N : usize> Gaussian<S, N>
@@ -169,8 +201,8 @@
                                                                        Gaussian<S, N>>;
 
 
-/// This implements $χ\_{[-b, b]^n} \* (f χ\_{[-a, a]^n})$
-/// where $a,b>0$ and $f$ is a gaussian kernel on $ℝ^n$.
+/// This implements $g := χ\_{[-b, b]^n} \* (f χ\_{[-a, a]^n})$ where $a,b>0$ and $f$ is
+/// a gaussian kernel on $ℝ^n$. For an expression for $g$, see Lemma 3.9 in the manuscript.
 #[replace_float_literals(F::cast_from(literal))]
 impl<'a, F : Float, R, C, S, const N : usize> Apply<&'a Loc<F, N>>
 for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>
@@ -224,6 +256,89 @@
     }
 }
 
+/// This implements the differential of $g := χ\_{[-b, b]^n} \* (f χ\_{[-a, a]^n})$ where $a,b>0$
+/// and $f$ is a gaussian kernel on $ℝ^n$. For an expression for the value of $g$, from which the
+/// derivative readily arises (at points of differentiability), see Lemma 3.9 in the manuscript.
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F : Float, R, C, S, const N : usize> Differentiable<&'a Loc<F, N>>
+for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>
+where R : Constant<Type=F>,
+      C : Constant<Type=F>,
+      S : Constant<Type=F> {
+
+    type Output = Loc<F, N>;
+
+    #[inline]
+    fn differential(&self, y : &'a Loc<F, N>) -> Loc<F, N> {
+        let Convolution(ref ind,
+                        SupportProductFirst(ref cut,
+                                            ref gaussian)) = self;
+        let a = cut.r.value();
+        let b = ind.r.value();
+        let σ = gaussian.variance.value().sqrt();
+        let π = F::PI;
+        let t = F::SQRT_2 * σ;
+        let c = σ * (8.0/π).sqrt();
+        let cd = (8.0).sqrt(); // σ * (8.0/π).sqrt() / t * (√2/π)
+        
+        // Calculate the values for all component functions of the
+        // product. This is just the loop from apply above.
+        let unscaled_vs = y.map(|x| {
+            let c1 = -(a.min(b + x)); //(-a).max(-x-b);
+            let c2 = a.min(b - x);
+            if c1 >= c2 {
+                0.0
+            } else {
+                let e1 = F::cast_from(erf((c1 / t).as_()));
+                let e2 = F::cast_from(erf((c2 / t).as_()));
+                debug_assert!(e2 >= e1);
+                c * (e2 - e1)
+            }
+        });
+        // This computes the gradient for each coordinate
+        product_differential(y, &unscaled_vs, |x| {
+            let c1 = -(a.min(b + x)); //(-a).max(-x-b);
+            let c2 = a.min(b - x);
+            if c1 >= c2 {
+                0.0
+            } else {
+                // erf'(z) = (2/√π)*exp(-z^2), and we get extra factor -1/√(2*σ) = -1/t
+                // from the chain rule
+                let de1 = (-(c1/t).powi(2)).exp();
+                let de2 = (-(c2/t).powi(2)).exp();
+                cd * (de1 - de2)
+            }
+        }) / gaussian.scale()
+    }
+}
+
+impl<F : Float, R, C, S, const N : usize> Differentiable<Loc<F, N>>
+for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>
+where R : Constant<Type=F>,
+      C : Constant<Type=F>,
+      S : Constant<Type=F> {
+
+    type Output = Loc<F, N>;
+
+    #[inline]
+    fn differential(&self, y : Loc<F, N>) -> Loc<F, N> {
+        self.differential(&y)
+    }
+}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F : Float, R, C, S, const N : usize> Lipschitz<L2>
+for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>
+where R : Constant<Type=F>,
+      C : Constant<Type=F>,
+      S : Constant<Type=F> {
+    type FloatType = F;
+
+    fn lipschitz_factor(&self, L2 : L2) -> Option<F> {
+        todo!("This requirement some error function work.")
+    }
+}
+
 impl<F : Float, R, C, S, const N : usize>
 Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>
 where R : Constant<Type=F>,