diff -r e2953ffd4e0b -r e9a460a0e638 src/kernels/hat_convolution.rs
--- a/src/kernels/hat_convolution.rs	Fri May 15 14:40:02 2026 -0500
+++ b/src/kernels/hat_convolution.rs	Sun Jul 19 07:34:39 2026 +0200
@@ -97,42 +97,30 @@
     #[inline]
     fn apply<I: Instance<Loc<N, S::Type>>>(&self, y: I) -> Self::Codomain {
         let σ = self.radius();
-        y.decompose().product_map(|x| self.value_1d_σ1(x / σ) / σ)
+        y.decompose().product_map(|x| self.value_1d_σ1(x / σ)) / σ.powi(N as i32)
     }
 }
 
 #[replace_float_literals(S::Type::cast_from(literal))]
-impl<S, const N: usize> Lipschitz<L1> for HatConv<S, N>
-where
-    S: Constant,
-{
-    type FloatType = S::Type;
-    #[inline]
-    fn lipschitz_factor(&self, L1: L1) -> DynResult<Self::FloatType> {
-        // For any ψ_i, we have
-        // ∏_{i=1}^N ψ_i(x_i) - ∏_{i=1}^N ψ_i(y_i)
-        // = [ψ_1(x_1)-ψ_1(y_1)] ∏_{i=2}^N ψ_i(x_i)
-        //   + ψ_1(y_1)[ ∏_{i=2}^N ψ_i(x_i) - ∏_{i=2}^N ψ_i(y_i)]
-        // = ∑_{j=1}^N [ψ_j(x_j)-ψ_j(y_j)]∏_{i > j} ψ_i(x_i) ∏_{i < j} ψ_i(y_i)
-        // Thus
-        // |∏_{i=1}^N ψ_i(x_i) - ∏_{i=1}^N ψ_i(y_i)|
-        // ≤ ∑_{j=1}^N |ψ_j(x_j)-ψ_j(y_j)| ∏_{j ≠ i} \max_j |ψ_j|
-        let σ = self.radius();
-        let l1d = self.lipschitz_1d_σ1() / (σ * σ);
-        let m1d = self.value_1d_σ1(0.0) / σ;
-        Ok(l1d * m1d.powi(N as i32 - 1))
-    }
-}
-
 impl<S, const N: usize> Lipschitz<L2> for HatConv<S, N>
 where
     S: Constant,
 {
     type FloatType = S::Type;
     #[inline]
-    fn lipschitz_factor(&self, L2: L2) -> DynResult<Self::FloatType> {
-        self.lipschitz_factor(L1)
-            .map(|l1| l1 * <S::Type>::cast_from(N).sqrt())
+    fn lipschitz_factor(&self, _L2: L2) -> DynResult<Self::FloatType> {
+        // For arbitrary ψ(x) = ∏_{i=1}^n ψ_i(x_i), we have
+        // ψ(x) - ψ(y) = ∑_i [ψ_i(x_i)-ψ_i(y_i)] ∏_{j < i} ψ_j(y_j) ∏_{j > i} ψ_j(x_j)
+        // by a simple recursive argument. In particular, if ψ_i=g for all i, j, we have
+        // |ψ(x) - ψ(y)| ≤  ∑_i L_g|x_i-y_i| M_g^{n-1} ≤  √n L_g M_g^{n-1} |x-y|₂
+        // where L_g is the Lipschitz factor of g, and M_g a bound on it.
+        let σ = self.radius();
+        let l1d = self.lipschitz_1d_σ1();
+        let m1d = self.value_1d_σ1(0.0);
+        Ok(
+            S::Type::cast_from(N).sqrt() * l1d * m1d.powi((N - 1) as i32)
+                / σ.powi((N - 1 + 2) as i32),
+        )
     }
 }
 
@@ -146,9 +134,9 @@
     fn differential_impl<I: Instance<Loc<N, S::Type>>>(&self, y0: I) -> Self::Derivative {
         let y = y0.decompose();
         let σ = self.radius();
-        let σ2 = σ * σ;
-        let vs = y.map(|x| self.value_1d_σ1(x / σ) / σ);
-        product_differential(&*y, &vs, |x| self.diff_1d_σ1(x / σ) / σ2)
+        let tmp = y.map(|x| x / σ);
+        let vs = tmp.map(|x| self.value_1d_σ1(x));
+        product_differential(&tmp, &vs, |x| self.diff_1d_σ1(x)) / σ.powi((N + 1) as i32)
     }
 }
 
@@ -164,11 +152,14 @@
     fn diff_lipschitz_factor(&self, _l2: L2) -> DynResult<F> {
         let σ = self.radius();
         product_differential_lipschitz_factor::<F, N>(
-            self.value_1d_σ1(0.0) / σ,
-            self.lipschitz_1d_σ1() / (σ * σ),
-            self.maxabsdiff_1d_σ1() / (σ * σ),
-            self.lipschitz_diff_1d_σ1() / (σ * σ),
+            self.value_1d_σ1(0.0),
+            self.lipschitz_1d_σ1(),
+            self.maxabsdiff_1d_σ1(),
+            self.lipschitz_diff_1d_σ1(),
         )
+        .map(|l| l / σ.powi((N + 2) as i32))
+        // We have f(x)=f₀(x/σ)/σ^N, where f₀(z)=∏_{i=1}^N f₁(z_i).
+        // Thus |f'(x)-f'(y)| = |f₀'(x/σ)-f₀'(y/σ)|/σ^{N+1} ≤ L_{f₀'}/σ^{N+2}.
     }
 }
 
@@ -195,15 +186,15 @@
     /// Computes the differential of the kernel for $n=1$ with $σ=1$.
     #[inline]
     fn diff_1d_σ1(&self, x: F) -> F {
-        let y = x.abs();
-        if y >= 1.0 {
+        if x >= 1.0 || x <= -1.0 {
             0.0
-        } else if y > 0.5 {
-            -8.0 * (y - 1.0).powi(2)
-        } else
-        /* 0 ≤ y ≤ 0.5 */
-        {
-            (24.0 * y - 16.0) * y
+        } else if x > 0.5 {
+            -8.0 * (x - 1.0).powi(2)
+        } else if x < -0.5 {
+            8.0 * (x + 1.0).powi(2)
+        } else {
+            /* 0 ≤ y ≤ 0.5 */
+            (24.0 * x.abs() - 16.0) * x
         }
     }
 
@@ -221,22 +212,6 @@
         2.0
     }
 
-    /// Computes the second differential of the kernel for $n=1$ with $σ=1$.
-    #[inline]
-    #[allow(dead_code)]
-    fn diff2_1d_σ1(&self, x: F) -> F {
-        let y = x.abs();
-        if y >= 1.0 {
-            0.0
-        } else if y > 0.5 {
-            -16.0 * (y - 1.0)
-        } else
-        /* 0 ≤ y ≤ 0.5 */
-        {
-            48.0 * y - 16.0
-        }
-    }
-
     /// Computes the differential of the kernel for $n=1$ with $σ=1$.
     #[inline]
     fn lipschitz_diff_1d_σ1(&self) -> F {
@@ -332,7 +307,7 @@
             //      = ∫_{x-β}^{x+β} u_σ(z) d z
             //      = (1/σ)∫_{x-β}^{x+β} u_1(z/σ) d z
             //      = ∫_{(x-β)/σ}^{(x+β)/σ} u_1(z) d z
-            //      = [χ_{-β/σ, β/σ} * u_1](x/σ)
+            //      = (χ_{[-β/σ, β/σ]} * u_1)(x/σ)
             // $$
             self.value_1d_σ1(x / σ, β / σ)
         })
@@ -354,14 +329,14 @@
         let Convolution(ref ind, ref hatconv) = self;
         let β = ind.r.value();
         let σ = hatconv.radius();
-        let σ2 = σ * σ;
 
         let vs = y.map(|x| self.value_1d_σ1(x / σ, β / σ));
-        product_differential(&*y, &vs, |x| self.diff_1d_σ1(x / σ, β / σ) / σ2)
+        product_differential(&*y, &vs, |x| self.diff_1d_σ1(x / σ, β / σ)) / σ
     }
 }
 
-/// Integrate $f$, whose support is $[c, d]$, on $[a, b]$.
+/// Integrate $f'$, whose support is $[c, d]$, on $[a, b]$.
+/// The value $f$ is given, not the derivative.
 /// If $b > d$, add $g()$ to the result.
 #[inline]
 #[replace_float_literals(F::cast_from(literal))]
@@ -403,12 +378,6 @@
         let a = x - β;
         let b = x + β;
 
-        #[inline]
-        fn pow4<F: Float>(x: F) -> F {
-            let y = x * x;
-            y * y
-        }
-
         // Observe the factor 1/6 at the front from the antiderivatives below.
         // The factor 4 is from normalisation of the original function.
         (4.0 / 6.0)
@@ -419,37 +388,28 @@
                 -0.5,
                 // (2/3) (y+1)^3  on  -1 < y ≤ -1/2
                 // The antiderivative is  (2/12)(y+1)^4 = (1/6)(y+1)^4
-                |y| pow4(y + 1.0),
+                |y| (y + 1.0).powi(4),
                 || {
                     i(
                         a,
                         b,
                         -0.5,
-                        0.0,
+                        0.5,
                         // -2 y^3 - 2 y^2 + 1/3  on  -1/2 < y ≤ 0
                         // The antiderivative is -1/2 y^4 - 2/3 y^3 + 1/3 y
-                        |y| y * (-y * y * (y * 3.0 + 4.0) + 2.0),
+                        // 2 y^3 - 2 y^2 + 1/3 on 0 < y < 1/2
+                        // The antiderivative is 1/2 y^4 - 2/3 y^3 + 1/3 y
+                        |y| y * (y * y * (y.abs() * 3.0 - 4.0) + 2.0),
                         || {
                             i(
                                 a,
                                 b,
-                                0.0,
                                 0.5,
-                                // 2 y^3 - 2 y^2 + 1/3 on 0 < y < 1/2
-                                // The antiderivative is 1/2 y^4 - 2/3 y^3 + 1/3 y
-                                |y| y * (y * y * (y * 3.0 - 4.0) + 2.0),
-                                || {
-                                    i(
-                                        a,
-                                        b,
-                                        0.5,
-                                        1.0,
-                                        // -(2/3) (y-1)^3  on  1/2 < y ≤ 1
-                                        // The antiderivative is  -(2/12)(y-1)^4 = -(1/6)(y-1)^4
-                                        |y| -pow4(y - 1.0),
-                                        || 0.0,
-                                    )
-                                },
+                                1.0,
+                                // -(2/3) (y-1)^3  on  1/2 < y ≤ 1
+                                // The antiderivative is  -(2/12)(y-1)^4 = -(1/6)(y-1)^4
+                                |y| -(y - 1.0).powi(4),
+                                || 0.0,
                             )
                         },
                     )
@@ -466,41 +426,30 @@
         let a = x - β;
         let b = x + β;
 
-        // The factor 4 is from normalisation of the original function.
-        4.0 * i(
+        i(
             a,
             b,
             -1.0,
             -0.5,
             // (2/3) (y+1)^3  on  -1 < y ≤ -1/2
-            |y| (2.0 / 3.0) * (y + 1.0).powi(3),
+            |y| (8.0 / 3.0) * (y + 1.0).powi(3),
             || {
                 i(
                     a,
                     b,
                     -0.5,
-                    0.0,
-                    // -2 y^3 - 2 y^2 + 1/3  on  -1/2 < y ≤ 0
-                    |y| -2.0 * (y + 1.0) * y * y + (1.0 / 3.0),
+                    0.5,
+                    // 2 y^3 - 2 y^2 + 1/3 on 0 < y < 1/2
+                    |y| 8.0 * (y.abs() - 1.0) * y * y + (4.0 / 3.0),
                     || {
                         i(
                             a,
                             b,
-                            0.0,
                             0.5,
-                            // 2 y^3 - 2 y^2 + 1/3 on 0 < y < 1/2
-                            |y| 2.0 * (y - 1.0) * y * y + (1.0 / 3.0),
-                            || {
-                                i(
-                                    a,
-                                    b,
-                                    0.5,
-                                    1.0,
-                                    // -(2/3) (y-1)^3  on  1/2 < y ≤ 1
-                                    |y| -(2.0 / 3.0) * (y - 1.0).powi(3),
-                                    || 0.0,
-                                )
-                            },
+                            1.0,
+                            // -(2/3) (y-1)^3  on  1/2 < y ≤ 1
+                            |y| -(8.0 / 3.0) * (y - 1.0).powi(3),
+                            || 0.0,
                         )
                     },
                 )
@@ -525,6 +474,7 @@
 }
 */
 
+#[replace_float_literals(F::cast_from(literal))]
 impl<F: Float, R, C, const N: usize> Convolution<CubeIndicator<R, N>, HatConv<C, N>>
 where
     R: Constant<Type = F>,