diff -r 000000000000 -r eb3c7813b67a src/kernels/mollifier.rs
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/kernels/mollifier.rs	Thu Dec 01 23:07:35 2022 +0200
@@ -0,0 +1,136 @@
+
+//! Implementation of the standard mollifier
+
+use rgsl::hypergeometric::hyperg_U;
+use float_extras::f64::{tgamma as gamma};
+use numeric_literals::replace_float_literals;
+use serde::Serialize;
+use alg_tools::types::*;
+use alg_tools::euclidean::Euclidean;
+use alg_tools::norms::*;
+use alg_tools::loc::Loc;
+use alg_tools::sets::Cube;
+use alg_tools::bisection_tree::{
+    Support,
+    Constant,
+    Bounds,
+    LocalAnalysis,
+    GlobalAnalysis
+};
+use alg_tools::mapping::Apply;
+use alg_tools::maputil::array_init;
+
+/// Reresentation of the (unnormalised) standard mollifier.
+///
+/// For the `width` parameter $ε>0$, this is
+/// <div>$$
+///     f(x)=\begin{cases}
+///         e^{\frac{ε^2}{\|x\|_2^2-ε^2}}, & \|x\|_2 < ε, \\
+///         0, & \text{otherwise}.
+///     \end{cases}
+/// $$</div>
+#[derive(Copy,Clone,Serialize,Debug,Eq,PartialEq)]
+pub struct Mollifier<C : Constant, const N : usize> {
+    /// The parameter $ε$ of the mollifier.
+    pub width : C,
+}
+
+#[replace_float_literals(C::Type::cast_from(literal))]
+impl<'a, C : Constant, const N : usize> Apply<&'a Loc<C::Type, N>> for Mollifier<C, N> {
+    type Output = C::Type;
+    #[inline]
+    fn apply(&self, x : &'a Loc<C::Type, N>) -> Self::Output {
+        let ε = self.width.value();
+        let ε2 = ε*ε;
+        let n2 = x.norm2_squared();
+        if n2 < ε2 {
+            (n2 / (n2 - ε2)).exp()
+        } else {
+            0.0
+        }
+    }
+}
+
+impl<C : Constant, const N : usize> Apply<Loc<C::Type, N>> for Mollifier<C, N> {
+    type Output = C::Type;
+    #[inline]
+    fn apply(&self, x : Loc<C::Type, N>) -> Self::Output {
+        self.apply(&x)
+    }
+}
+
+impl<'a, C : Constant, const N : usize> Support<C::Type, N> for Mollifier<C, N> {
+    #[inline]
+    fn support_hint(&self) -> Cube<C::Type,N> {
+        let ε = self.width.value();
+        array_init(|| [-ε, ε]).into()
+    }
+
+    #[inline]
+    fn in_support(&self, x : &Loc<C::Type,N>) -> bool {
+        x.norm2() < self.width.value()
+    }
+    
+    /*fn fully_in_support(&self, _cube : &Cube<C::Type,N>) -> bool {
+        todo!("Not implemented, but not used at the moment")
+    }*/
+}
+
+#[replace_float_literals(C::Type::cast_from(literal))]
+impl<'a, C : Constant, const N : usize> GlobalAnalysis<C::Type, Bounds<C::Type>>
+for Mollifier<C, N> {
+    #[inline]
+    fn global_analysis(&self) -> Bounds<C::Type> {
+        // The function is maximised/minimised where the 2-norm is minimised/maximised.
+        Bounds(0.0, 1.0)
+    }
+}
+
+impl<'a, C : Constant, const N : usize> LocalAnalysis<C::Type, Bounds<C::Type>, N>
+for Mollifier<C, N> {
+    #[inline]
+    fn local_analysis(&self, cube : &Cube<C::Type, N>) -> Bounds<C::Type> {
+        // The function is maximised/minimised where the 2-norm is minimised/maximised.
+        let lower = self.apply(cube.maxnorm_point());
+        let upper = self.apply(cube.minnorm_point());
+        Bounds(lower, upper)
+    }
+}
+
+/// Calculate integral of the standard mollifier of width 1 in $ℝ^n$.
+///
+/// This is based on the formula from
+/// [https://math.stackexchange.com/questions/4359683/integral-of-the-usual-mollifier-function-finding-its-necessary-constant]().
+///
+/// If `rescaled` is `true`, return the integral of the scaled mollifier that has value one at the
+/// origin.
+#[inline]
+pub fn mollifier_norm1(n_ : usize, rescaled : bool) -> f64 {
+    assert!(n_ > 0);
+    let n = n_ as f64;
+    let q = 2.0;
+    let p = 2.0;
+    let base = (2.0*gamma(1.0 + 1.0/p)).powi(n_ as i32)
+               /*/ gamma(1.0 + n / p)
+               * gamma(1.0 + n / q)*/
+               * hyperg_U(1.0 + n / q, 2.0, 1.0);
+    if rescaled { base } else { base / f64::E }
+}
+
+impl<'a, C : Constant, const N : usize> Norm<C::Type, L1>
+for Mollifier<C, N> {
+    #[inline]
+    fn norm(&self, _ : L1) -> C::Type {
+        let ε = self.width.value();
+        C::Type::cast_from(mollifier_norm1(N, true)) * ε.powi(N as i32)
+    }
+}
+
+#[replace_float_literals(C::Type::cast_from(literal))]
+impl<'a, C : Constant, const N : usize> Norm<C::Type, Linfinity>
+for Mollifier<C, N> {
+    #[inline]
+    fn norm(&self, _ : Linfinity) -> C::Type {
+        1.0
+    }
+}