Mercurial > repos > pointsource_algs / file diff

--- a/src/kernels/base.rs	Fri May 15 14:40:02 2026 -0500
+++ b/src/kernels/base.rs	Sun Jul 19 07:34:39 2026 +0200
@@ -437,24 +437,28 @@
     dlip: F,
 ) -> DynResult<F> {
     // For arbitrary ψ(x) = ∏_{i=1}^n ψ_i(x_i), we have
-    // ψ(x) - ψ(y) = ∑_i [ψ_i(x_i)-ψ_i(y_i)] ∏_{j ≠ i} ψ_j(x_j)
+    // ψ(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 M_g^{n-1}|x-y|, where L_g is the Lipschitz factor of g, and
-    // M_g a bound on it.
+    // |ψ(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.
     //
-    // We also have in the general case ∇ψ(x) = ∑_i ∇ψ_i(x_i) ∏_{j ≠ i} ψ_j(x_j), whence
-    // using the previous formula for each i with f_i=∇ψ_i and f_j=ψ_j for j ≠ i, we get
-    //  ∇ψ(x) - ∇ψ(y) = ∑_i[ ∇ψ_i(x_i)∏_{j ≠ i} ψ_j(x_j) - ∇ψ_i(y_i)∏_{j ≠ i} ψ_j(y_j)]
-    //                = ∑_i[ [∇ψ_i(x_i) - ∇ψ_j(x_j)] ∏_{j ≠ i}ψ_j(x_j)
-    //                       + [∑_{k ≠ i} [ψ_k(x_k) - ∇ψ_k(x_k)] ∏_{j ≠ i, k}ψ_j(x_j)]∇ψ_i(x_i)].
-    // With $ψ_i=g for all i, j, it follows that
-    // |∇ψ(x) - ∇ψ(y)| ≤ ∑_i L_{∇g} M_g^{n-1} + ∑_{k ≠ i} L_g M_g^{n-2} M_{∇g}
-    //                 = n [L_{∇g} M_g^{n-1} + (n-1) L_g M_g^{n-2} M_{∇g}].
-    //                 = n M_g^{n-2}[L_{∇g} M_g + (n-1) L_g M_{∇g}].
-    if N >= 2 {
+    // We also have in the general case [∇ψ(x)]_i = ψ_i'(x_i) ∏_{j ≠ i} ψ_j(x_j), whence
+    // from above, the Lipschitz factor of [∇ψ(x)]_i is
+    // L' = √(L_{g'}^2M_g^{2(n-1)} + (n-1)L_g^2M_{g'}^2M_g^{2(n-2)}) if n > 2,
+    // in particular
+    // L' = √(L_{g'}^2M_g^{2(n-1)} + L_g^2M_{g'}^2) if n > 2, and
+    // Now the Lipschitz factor of ∇ψ is √n L' from
+    // ‖∇ψ(x)-∇ψ(y)‖ = √(∑_{i=1}^n [∇ψ(x)-∇ψ(y)]_i^2) ≤  √(∑_{i=1}^n (L')^2) ≤ √n L'
+    if N > 2 {
         Ok(F::cast_from(N)
-            * bound.powi((N - 2) as i32)
-            * (dlip * bound + F::cast_from(N - 1) * lip * dbound))
+            * (dlip.powi(2) * bound.powi(2 * (N - 1) as i32)
+                + F::cast_from(N - 1)
+                    * lip.powi(2)
+                    * dbound.powi(2)
+                    * bound.powi(2 * (N - 2) as i32))
+            .sqrt())
+    } else if N == 2 {
+        Ok((F::TWO * ((dlip * bound).powi(2) + (lip * dbound).powi(2))).sqrt())
     } else if N == 1 {
         Ok(dlip)
     } else {