diff -r 000000000000 -r eb3c7813b67a src/kernels/ball_indicator.rs
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/kernels/ball_indicator.rs	Thu Dec 01 23:07:35 2022 +0200
@@ -0,0 +1,260 @@
+
+//! Implementation of the indicator function of a ball with respect to various norms.
+use float_extras::f64::{tgamma as gamma};
+use numeric_literals::replace_float_literals;
+use serde::Serialize;
+use alg_tools::types::*;
+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;
+use alg_tools::coefficients::factorial;
+
+use super::base::*;
+
+/// Representation of the indicator of the ball $𝔹_q = \\{ x ∈ ℝ^N \mid \\|x\\|\_q ≤ r \\}$,
+/// where $q$ is the `Exponent`, and $r$ is the radius [`Constant`] `C`.
+#[derive(Copy,Clone,Serialize,Debug,Eq,PartialEq)]
+pub struct BallIndicator<C : Constant, Exponent : NormExponent, const N : usize> {
+    /// The radius of the ball.
+    pub r : C,
+    /// The exponent $q$ of the norm creating the ball
+    pub exponent : Exponent,
+}
+
+/// Alias for the representation of the indicator of the $∞$-norm-ball
+/// $𝔹_∞ = \\{ x ∈ ℝ^N \mid \\|x\\|\_∞ ≤ c \\}$.
+pub type CubeIndicator<C, const N : usize> = BallIndicator<C, Linfinity, N>;
+
+#[replace_float_literals(C::Type::cast_from(literal))]
+impl<'a, F : Float, C : Constant<Type=F>, Exponent : NormExponent, const N : usize>
+Apply<&'a Loc<C::Type, N>>
+for BallIndicator<C, Exponent, N>
+where Loc<F, N> : Norm<F, Exponent> {
+    type Output = C::Type;
+    #[inline]
+    fn apply(&self, x : &'a Loc<C::Type, N>) -> Self::Output {
+        let r = self.r.value();
+        let n = x.norm(self.exponent);
+        if n <= r {
+            1.0
+        } else {
+            0.0
+        }
+    }
+}
+
+impl<F : Float, C : Constant<Type=F>, Exponent : NormExponent, const N : usize>
+Apply<Loc<C::Type, N>>
+for BallIndicator<C, Exponent, N>
+where Loc<F, N> : Norm<F, Exponent> {
+    type Output = C::Type;
+    #[inline]
+    fn apply(&self, x : Loc<C::Type, N>) -> Self::Output {
+        self.apply(&x)
+    }
+}
+
+
+impl<'a, F : Float, C : Constant<Type=F>, Exponent : NormExponent, const N : usize>
+Support<C::Type, N>
+for BallIndicator<C, Exponent, N>
+where Loc<F, N> : Norm<F, Exponent>,
+      Linfinity : Dominated<F, Exponent, Loc<F, N>> {
+
+    #[inline]
+    fn support_hint(&self) -> Cube<F,N> {
+        let r = Linfinity.from_norm(self.r.value(), self.exponent);
+        array_init(|| [-r, r]).into()
+    }
+
+    #[inline]
+    fn in_support(&self, x : &Loc<F,N>) -> bool {
+        let r = Linfinity.from_norm(self.r.value(), self.exponent);
+        x.norm(self.exponent) <= r
+    }
+
+    /// This can only really work in a reasonable fashion for N=1.
+    #[inline]
+    fn bisection_hint(&self, cube : &Cube<F, N>) -> [Option<F>; N] {
+        let r = Linfinity.from_norm(self.r.value(), self.exponent);
+        cube.map(|a, b| symmetric_interval_hint(r, a, b))
+    }
+}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F : Float, C : Constant<Type=F>, Exponent : NormExponent, const N : usize>
+GlobalAnalysis<F, Bounds<F>>
+for BallIndicator<C, Exponent, N>
+where Loc<F, N> : Norm<F, Exponent> {
+    #[inline]
+    fn global_analysis(&self) -> Bounds<F> {
+        Bounds(0.0, 1.0)
+    }
+}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F : Float, C : Constant<Type=F>, Exponent : NormExponent, const N : usize>
+Norm<F, Linfinity>
+for BallIndicator<C, Exponent, N>
+where Loc<F, N> : Norm<F, Exponent> {
+    #[inline]
+    fn norm(&self, _ : Linfinity) -> F {
+        1.0
+    }
+}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F : Float, C : Constant<Type=F>, const N : usize>
+Norm<F, L1>
+for BallIndicator<C, L1, N> {
+    #[inline]
+    fn norm(&self, _ : L1) -> F {
+        // Using https://en.wikipedia.org/wiki/Volume_of_an_n-ball#Balls_in_Lp_norms,
+        // we have V_N^1(r) = (2r)^N / N!
+        let r = self.r.value();
+        if N==1 {
+            2.0 * r
+        } else if N==2 {
+            r*r
+        } else {
+            (2.0 * r).powi(N as i32) * F::cast_from(factorial(N))
+        }
+    }
+}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F : Float, C : Constant<Type=F>, const N : usize>
+Norm<F, L1>
+for BallIndicator<C, L2, N> {
+    #[inline]
+    fn norm(&self, _ : L1) -> F {
+        // See https://en.wikipedia.org/wiki/Volume_of_an_n-ball#The_volume.
+        let r = self.r.value();
+        let π = F::PI;
+        if N==1 {
+            2.0 * r
+        } else if N==2 {
+            π * (r * r)
+        } else {
+            let ndiv2 = F::cast_from(N) / 2.0;
+            let γ = F::cast_from(gamma((ndiv2 + 1.0).as_()));
+            π.powf(ndiv2) / γ * r.powi(N as i32)
+        }
+    }
+}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F : Float, C : Constant<Type=F>, const N : usize>
+Norm<F, L1>
+for BallIndicator<C, Linfinity, N> {
+    #[inline]
+    fn norm(&self, _ : L1) -> F {
+        let two_r = 2.0 * self.r.value();
+        two_r.powi(N as i32)
+    }
+}
+
+
+macro_rules! indicator_local_analysis {
+    ($exponent:ident) => {
+        impl<'a, F : Float, C : Constant<Type=F>, const N : usize>
+        LocalAnalysis<F, Bounds<F>, N>
+        for BallIndicator<C, $exponent, N>
+        where Loc<F, N> : Norm<F, $exponent>,
+            Linfinity : Dominated<F, $exponent, Loc<F, N>> {
+            #[inline]
+            fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> {
+                // 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)
+            }
+        }
+    }
+}
+
+indicator_local_analysis!(L1);
+indicator_local_analysis!(L2);
+indicator_local_analysis!(Linfinity);
+
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F : Float, R, const N : usize> Apply<&'a Loc<F, N>>
+for AutoConvolution<CubeIndicator<R, N>>
+where R : Constant<Type=F> {
+    type Output = F;
+
+    #[inline]
+    fn apply(&self, y : &'a Loc<F, N>) -> F {
+        let two_r = 2.0 * self.0.r.value();
+        // This is just a product of one-dimensional versions
+        y.iter().map(|&x| {
+            0.0.max(two_r - x.abs())
+        }).product()
+    }
+}
+
+impl<F : Float, R, const N : usize> Apply<Loc<F, N>>
+for AutoConvolution<CubeIndicator<R, N>>
+where R : Constant<Type=F> {
+    type Output = F;
+
+    #[inline]
+    fn apply(&self, y : Loc<F, N>) -> F {
+        self.apply(&y)
+    }
+}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<F : Float, R, const N : usize> Support<F, N>
+for AutoConvolution<CubeIndicator<R, N>>
+where R : Constant<Type=F> {
+    #[inline]
+    fn support_hint(&self) -> Cube<F, N> {
+        let two_r = 2.0 * self.0.r.value();
+        array_init(|| [-two_r, two_r]).into()
+    }
+
+    #[inline]
+    fn in_support(&self, y : &Loc<F, N>) -> bool {
+        let two_r = 2.0 * self.0.r.value();
+        y.iter().all(|x| x.abs() <= two_r)
+    }
+
+    #[inline]
+    fn bisection_hint(&self, cube : &Cube<F, N>) -> [Option<F>; N] {
+        let two_r = 2.0 * self.0.r.value();
+        cube.map(|c, d| symmetric_interval_hint(two_r, c, d))
+    }
+}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<F : Float, R, const N : usize> GlobalAnalysis<F, Bounds<F>>
+for AutoConvolution<CubeIndicator<R, N>>
+where R : Constant<Type=F> {
+    #[inline]
+    fn global_analysis(&self) -> Bounds<F> {
+        Bounds(0.0, self.apply(Loc::ORIGIN))
+    }
+}
+
+impl<F : Float, R, const N : usize> LocalAnalysis<F, Bounds<F>, N>
+for AutoConvolution<CubeIndicator<R, N>>
+where R : Constant<Type=F> {
+    #[inline]
+    fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> {
+        // The function is maximised/minimised where the absolute value is minimised/maximised.
+        let lower = self.apply(cube.maxnorm_point());
+        let upper = self.apply(cube.minnorm_point());
+        Bounds(lower, upper)
+    }
+}