# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1735654363 18000
# Node ID d2acaaddd9af19198a530d09ebd4116588f91122
# Parent  edb95d2b83ccb541d736e1742643068c95ef32f7
Try to have Field as member type in Mappings etc.

diff -r edb95d2b83cc -r d2acaaddd9af src/bisection_tree/bt.rs
--- a/src/bisection_tree/bt.rs	Sun Nov 10 09:02:57 2024 -0500
+++ b/src/bisection_tree/bt.rs	Tue Dec 31 09:12:43 2024 -0500
@@ -243,10 +243,11 @@
 
     /// Creates a new node branching structure, subdividing `domain` based on the
     /// [hint][Support::support_hint] of `support`.
-    pub(super) fn new_with<S : LocalAnalysis <F, A, N>>(
+    pub(super) fn new_with<S>(
         domain : &Cube<F,N>,
         support : &S
-    ) -> Self {
+    ) -> Self
+    where S : Support<N, RealField=F> + LocalAnalysis<A, Cube<F, N>> {
         let hint = support.bisection_hint(domain);
         let branch_at = map2(&hint, domain, |h, r| {
             h.unwrap_or_else(|| (r[0]+r[1])/F::TWO).max(r[0]).min(r[1])
@@ -321,14 +322,15 @@
     ///  * `support` is the [`Support`] that is used determine with which subcubes of `domain`
     ///     (at subdivision depth `new_leaf_depth`) the data `d` is to be associated with.
     ///
-    pub(super) fn insert<'refs, 'scope, M : Depth, S : LocalAnalysis<F, A, N>>(
+    pub(super) fn insert<'refs, 'scope, M : Depth, S>(
         &mut self,
         domain : &Cube<F,N>,
         d : D,
         new_leaf_depth : M,
         support : &S,
         task_budget : TaskBudget<'scope, 'refs>,
-    ) {
+    )
+    where S : Support<N, RealField=F> + LocalAnalysis<A, Cube<F, N>> {
         let support_hint = support.support_hint();
         self.recurse(domain, task_budget,
                      |_, subcube| support_hint.intersects(&subcube),
@@ -348,8 +350,8 @@
         domain : &Cube<F, N>
     ) -> Branches<F,D,ANew,N,P>
     where ANew : Aggregator,
-          G : SupportGenerator<F, N, Id=D>,
-          G::SupportType : LocalAnalysis<F, ANew, N> {
+          G : SupportGenerator<N, Id=D, RealField = F>,
+          G::SupportType : LocalAnalysis<ANew, Cube<F, N>> {
         let branch_at = self.branch_at;
         let subcube_iter = self.iter_subcubes(domain);
         let new_nodes = self.nodes.into_iter().zip(subcube_iter).map(|(node, subcube)| {
@@ -370,8 +372,8 @@
         generator : &G,
         domain : &Cube<F, N>,
         task_budget : TaskBudget<'scope, 'refs>,
-    ) where G : SupportGenerator<F, N, Id=D>,
-            G::SupportType : LocalAnalysis<F, A, N> {
+    ) where G : SupportGenerator<N, Id=D, RealField = F>,
+            G::SupportType : LocalAnalysis<A, Cube<F, N>> {
         self.recurse(domain, task_budget,
                      |_, _| true,
                      move |node, subcube, new_budget| node.refresh_aggregator(generator, subcube,
@@ -434,14 +436,15 @@
     /// If `self` is a [`NodeOption::Branches`], the data is passed to branches whose subcubes
     /// `support` intersects. If an [`NodeOption::Uninitialised`] node is encountered, a new leaf is
     /// created at a minimum depth of `new_leaf_depth`.
-    pub(super) fn insert<'refs, 'scope, M : Depth, S : LocalAnalysis <F, A, N>>(
+    pub(super) fn insert<'refs, 'scope, M : Depth, S>(
         &mut self,
         domain : &Cube<F,N>,
         d : D,
         new_leaf_depth : M,
         support : &S,
         task_budget : TaskBudget<'scope, 'refs>,
-    ) {
+    )
+    where S : Support<N, RealField = F> + LocalAnalysis<A, Cube<F, N>> {
         match &mut self.data {
             NodeOption::Uninitialised => {
                 // Replace uninitialised node with a leaf or a branch
@@ -493,8 +496,8 @@
         domain : &Cube<F, N>
     ) -> Node<F,D,ANew,N,P>
     where ANew : Aggregator,
-          G : SupportGenerator<F, N, Id=D>,
-          G::SupportType : LocalAnalysis<F, ANew, N> {
+          G : SupportGenerator<N, Id=D, RealField = F>,
+          G::SupportType : LocalAnalysis<ANew, Cube<F, N>> {
 
         // The mem::replace is needed due to the [`Drop`] implementation to extract self.data.
         match std::mem::replace(&mut self.data, NodeOption::Uninitialised) {
@@ -537,8 +540,8 @@
         generator : &G,
         domain : &Cube<F, N>,
         task_budget : TaskBudget<'scope, 'refs>,
-    ) where G : SupportGenerator<F, N, Id=D>,
-            G::SupportType : LocalAnalysis<F, A, N> {
+    ) where G : SupportGenerator<N, Id=D, RealField = F>,
+            G::SupportType : LocalAnalysis<A, Cube<F, N>> {
         match &mut self.data {
             NodeOption::Uninitialised => { },
             NodeOption::Leaf(v) => {
@@ -580,7 +583,7 @@
 /// Basic interface to a [`BT`] bisection tree.
 ///
 /// Further routines are provided by the [`BTSearch`][super::refine::BTSearch] trait.
-pub trait BTImpl<F : Float, const N : usize> : std::fmt::Debug + Clone + GlobalAnalysis<F, Self::Agg> {
+pub trait BTImpl<F : Float, const N : usize> : std::fmt::Debug + Clone + GlobalAnalysis<Self::Agg> {
     /// The data type stored in the tree
     type Data : 'static + Copy + Send + Sync;
     /// The depth type of the tree
@@ -595,11 +598,12 @@
     ///
     /// Every leaf node of the tree that intersects the `support` will contain a copy of
     /// `d`.
-    fn insert<S : LocalAnalysis<F, Self::Agg, N>>(
+    fn insert<S>(
         &mut self,
         d : Self::Data,
         support : &S
-    );
+    )
+    where S: Support<N, RealField=F> + LocalAnalysis<Self::Agg, Cube<F, N>>;
 
     /// Construct a new instance of the tree for a different aggregator
     ///
@@ -608,8 +612,8 @@
     fn convert_aggregator<ANew, G>(self, generator : &G)
     -> Self::Converted<ANew>
     where ANew : Aggregator,
-          G : SupportGenerator<F, N, Id=Self::Data>,
-          G::SupportType : LocalAnalysis<F, ANew, N>;
+          G : SupportGenerator< N, Id=Self::Data, RealField = F>,
+          G::SupportType : LocalAnalysis<ANew, Cube<F, N>>;
 
 
     /// Refreshes the aggregator of the three after possible changes to the support generator.
@@ -617,8 +621,8 @@
     /// The `generator` is used to convert the data of type [`Self::Data`] contained in the tree
     /// into corresponding [`Support`]s.
     fn refresh_aggregator<G>(&mut self, generator : &G)
-    where G : SupportGenerator<F, N, Id=Self::Data>,
-          G::SupportType : LocalAnalysis<F, Self::Agg, N>;
+    where G : SupportGenerator<N, Id=Self::Data, RealField = F>,
+          G::SupportType : LocalAnalysis<Self::Agg, Cube<F, N>>;
 
     /// Returns an iterator over all [`Self::Data`] items at the point `x` of the domain.
     fn iter_at(&self, x : &Loc<F,N>) -> std::slice::Iter<'_, Self::Data>;
@@ -672,11 +676,12 @@
             type Agg = A;
             type Converted<ANew> = BT<M,F,D,ANew,$n> where ANew : Aggregator;
 
-            fn insert<S: LocalAnalysis<F, A, $n>>(
+            fn insert<S>(
                 &mut self,
                 d : D,
                 support : &S
-            ) {
+            )
+            where S : Support<$n, RealField = F> + LocalAnalysis<A, Cube<F, $n>> {
                 with_task_budget(|task_budget|
                     self.topnode.insert(
                         &self.domain,
@@ -690,8 +695,8 @@
 
             fn convert_aggregator<ANew, G>(self, generator : &G) -> Self::Converted<ANew>
             where ANew : Aggregator,
-                  G : SupportGenerator<F, $n, Id=D>,
-                  G::SupportType : LocalAnalysis<F, ANew, $n> {
+                  G : SupportGenerator<$n, Id=D, RealField = F>,
+                  G::SupportType : LocalAnalysis<ANew, Cube<F, $n>> {
                 let topnode = self.topnode.convert_aggregator(generator, &self.domain);
 
                 BT {
@@ -702,8 +707,8 @@
             }
 
             fn refresh_aggregator<G>(&mut self, generator : &G)
-            where G : SupportGenerator<F, $n, Id=Self::Data>,
-                G::SupportType : LocalAnalysis<F, Self::Agg, $n> {
+            where G : SupportGenerator<$n, Id=Self::Data, RealField = F>,
+                  G::SupportType : LocalAnalysis<Self::Agg, Cube<F, $n>> {
                 with_task_budget(|task_budget|
                     self.topnode.refresh_aggregator(generator, &self.domain, task_budget)
                 )
@@ -726,7 +731,7 @@
             }
         }
 
-        impl<M,F,D,A> GlobalAnalysis<F,A> for BT<M,F,D,A,$n>
+        impl<M,F,D,A> GlobalAnalysis<A> for BT<M,F,D,A,$n>
         where M : Depth,
               F : Float,
               D : 'static + Copy + Send + Sync + std::fmt::Debug,
diff -r edb95d2b83cc -r d2acaaddd9af src/bisection_tree/btfn.rs
--- a/src/bisection_tree/btfn.rs	Sun Nov 10 09:02:57 2024 -0500
+++ b/src/bisection_tree/btfn.rs	Tue Dec 31 09:12:43 2024 -0500
@@ -1,9 +1,8 @@
 
 use numeric_literals::replace_float_literals;
 use std::iter::Sum;
-use std::marker::PhantomData;
 use std::sync::Arc;
-use crate::types::Float;
+use crate::types::{Float, HasScalarField, HasRealField};
 use crate::mapping::{Apply, Mapping, Differentiable};
 //use crate::linops::{Apply, Linear};
 use crate::sets::Set;
@@ -29,21 +28,23 @@
 /// for a [`PreBTFN`] that is only useful for vector space operations with a full [`BTFN`].
 #[derive(Clone,Debug)]
 pub struct BTFN<
-    F : Float,
-    G : SupportGenerator<F, N>,
+    G : SupportGenerator<N>,
     BT /*: BTImpl<F, N>*/,
     const N : usize
-> /*where G::SupportType : LocalAnalysis<F, A, N>*/ {
+> {
     bt : BT,
     generator : Arc<G>,
-    _phantoms : PhantomData<F>,
+}
+
+impl<G : SupportGenerator<N>, BT, const N : usize> HasScalarField for BTFN<G, BT, N> {
+    type Field = G::Field;
 }
 
 impl<F : Float, G, BT, const N : usize>
-BTFN<F, G, BT, N>
-where G : SupportGenerator<F, N, Id=BT::Data>,
-      G::SupportType : LocalAnalysis<F, BT::Agg, N>,
-      BT : BTImpl<F, N> {
+BTFN<G, BT, N>
+where G : SupportGenerator<N, RealField = F>,
+      G::SupportType : LocalAnalysis<BT::Agg, Cube<F, N>>,
+      BT : BTImpl<F, N, Data=G::Id> {
 
     /// Create a new BTFN from a support generator and a pre-initialised bisection tree.
     ///
@@ -60,7 +61,6 @@
         BTFN {
             bt : bt,
             generator : generator,
-            _phantoms : std::marker::PhantomData,
         }
     }
 
@@ -85,11 +85,11 @@
     /// The top node of the created [`BT`] will have the given `domain`.
     ///
     /// See the documentation for [`BTFN`] on the role of the `generator`.
-    pub fn construct(domain : Cube<F, N>, depth : BT::Depth, generator : G) -> Self {
+    pub fn construct(domain : Cube<G::Field, N>, depth : BT::Depth, generator : G) -> Self {
         Self::construct_arc(domain, depth, Arc::new(generator))
     }
 
-    fn construct_arc(domain : Cube<F, N>, depth : BT::Depth, generator : Arc<G>) -> Self {
+    fn construct_arc(domain : Cube<G::Field, N>, depth : BT::Depth, generator : Arc<G>) -> Self {
         let mut bt = BT::new(domain, depth);
         for (d, support) in generator.all_data() {
             bt.insert(d, &support);
@@ -101,10 +101,10 @@
     ///
     /// This will construct a [`BTFN`] with the same components and generator as the (consumed)
     /// `self`, but a new `BT` with [`Aggregator`]s of type `ANew`.
-    pub fn convert_aggregator<ANew>(self) -> BTFN<F, G, BT::Converted<ANew>, N>
+    pub fn convert_aggregator<ANew>(self) -> BTFN<G, BT::Converted<ANew>, N>
     where ANew : Aggregator,
-          G : SupportGenerator<F, N, Id=BT::Data>,
-          G::SupportType : LocalAnalysis<F, ANew, N> {
+          G : SupportGenerator<N>,
+          G::SupportType : LocalAnalysis<ANew, Cube<G::RealField, N>> {
         BTFN::new_arc(self.bt.convert_aggregator(&*self.generator), self.generator)
     }
 
@@ -121,16 +121,16 @@
 }
 
 impl<F : Float, G, BT, const N : usize>
-BTFN<F, G, BT, N>
-where G : SupportGenerator<F, N> {
+BTFN<G, BT, N>
+where G : SupportGenerator<N, RealField = F> {
     /// Change the [bisection tree][BTImpl] of the [`BTFN`] to a different one.
     ///
     /// This can be used to convert a [`PreBTFN`] to a full [`BTFN`], or the change
     /// the aggreagator; see also [`self.convert_aggregator`].
     pub fn instantiate<
         BTNew : BTImpl<F, N, Data=G::Id>,
-    > (self, domain : Cube<F, N>, depth : BTNew::Depth) -> BTFN<F, G, BTNew, N>
-    where G::SupportType : LocalAnalysis<F, BTNew::Agg, N>  {
+    > (self, domain : Cube<F, N>, depth : BTNew::Depth) -> BTFN<G, BTNew, N>
+    where G::SupportType : LocalAnalysis<BTNew::Agg, Cube<F, N>>  {
         BTFN::construct_arc(domain, depth, self.generator)
     }
 }
@@ -140,30 +140,29 @@
 /// Most BTFN methods are not available, but if a BTFN is going to be summed with another
 /// before other use, it will be more efficient to not construct an unnecessary bisection tree
 /// that would be shortly dropped.
-pub type PreBTFN<F, G, const N : usize> = BTFN<F, G, (), N>;
+pub type PreBTFN<G, const N : usize> = BTFN<G, (), N>;
 
-impl<F : Float, G, const N : usize> PreBTFN<F, G, N> where G : SupportGenerator<F, N> {
+impl<G, const N : usize> PreBTFN<G, N> where G : SupportGenerator<N> {
 
     /// Create a new [`PreBTFN`] with no bisection tree.
     pub fn new_pre(generator : G) -> Self {
         BTFN {
             bt : (),
             generator : Arc::new(generator),
-            _phantoms : std::marker::PhantomData,
         }
     }
 }
 
 impl<F : Float, G, BT, const N : usize>
-BTFN<F, G, BT, N>
-where G : SupportGenerator<F, N, Id=usize>,
-      G::SupportType : LocalAnalysis<F, BT::Agg, N>,
+BTFN<G, BT, N>
+where G : SupportGenerator<N, Id=usize, RealField = F>,
+      G::SupportType : LocalAnalysis<BT::Agg, Cube<F, N>>,
       BT : BTImpl<F, N, Data=usize> {
 
     /// Helper function for implementing [`std::ops::Add`].
-    fn add_another<G2>(&self, g2 : Arc<G2>) -> BTFN<F, BothGenerators<G, G2>, BT, N>
-    where G2 : SupportGenerator<F, N, Id=usize>,
-          G2::SupportType : LocalAnalysis<F, BT::Agg, N> {
+    fn add_another<G2>(&self, g2 : Arc<G2>) -> BTFN<BothGenerators<G, G2>, BT, N>
+    where G2 : SupportGenerator<N, Id=usize, RealField = F>,
+          G2::SupportType : LocalAnalysis<BT::Agg, Cube<F, N>> {
 
         let mut bt = self.bt.clone();
         let both = BothGenerators(Arc::clone(&self.generator), g2);
@@ -175,7 +174,6 @@
         BTFN {
             bt : bt,
             generator : Arc::new(both),
-            _phantoms : std::marker::PhantomData,
         }
     }
 }
@@ -183,54 +181,54 @@
 macro_rules! make_btfn_add {
     ($lhs:ty, $preprocess:path, $($extra_trait:ident)?) => {
         impl<'a, F : Float, G1, G2, BT1, BT2, const N : usize>
-        std::ops::Add<BTFN<F, G2, BT2, N>> for
+        std::ops::Add<BTFN<G2, BT2, N>> for
         $lhs
         where BT1 : BTImpl<F, N, Data=usize>,
-              G1 : SupportGenerator<F, N, Id=usize> + $($extra_trait)?,
-              G2 : SupportGenerator<F, N, Id=usize>,
-              G1::SupportType : LocalAnalysis<F, BT1::Agg, N>,
-              G2::SupportType : LocalAnalysis<F, BT1::Agg, N> {
-            type Output = BTFN<F, BothGenerators<G1, G2>, BT1, N>;
+              G1 : SupportGenerator<N, Id=usize, RealField = F> + $($extra_trait)?,
+              G2 : SupportGenerator<N, Id=usize, RealField = F>,
+              G1::SupportType : LocalAnalysis<BT1::Agg, Cube<F, N>>,
+              G2::SupportType : LocalAnalysis<BT1::Agg, Cube<F, N>> {
+            type Output = BTFN<BothGenerators<G1, G2>, BT1, N>;
             #[inline]
-            fn add(self, other : BTFN<F, G2, BT2, N>) -> Self::Output {
+            fn add(self, other : BTFN<G2, BT2, N>) -> Self::Output {
                 $preprocess(self).add_another(other.generator)
             }
         }
 
         impl<'a, 'b, F : Float, G1, G2,  BT1, BT2, const N : usize>
-        std::ops::Add<&'b BTFN<F, G2, BT2, N>> for
+        std::ops::Add<&'b BTFN<G2, BT2, N>> for
         $lhs
         where BT1 : BTImpl<F, N, Data=usize>,
-              G1 : SupportGenerator<F, N, Id=usize> + $($extra_trait)?,
-              G2 : SupportGenerator<F, N, Id=usize>,
-              G1::SupportType : LocalAnalysis<F, BT1::Agg, N>,
-              G2::SupportType : LocalAnalysis<F, BT1::Agg, N> {
+              G1 : SupportGenerator<N, Id=usize, RealField = F> + $($extra_trait)?,
+              G2 : SupportGenerator<N, Id=usize, RealField = F>,
+              G1::SupportType : LocalAnalysis<BT1::Agg, Cube<F, N>>,
+              G2::SupportType : LocalAnalysis<BT1::Agg, Cube<F, N>> {
 
-            type Output = BTFN<F, BothGenerators<G1, G2>, BT1, N>;
+            type Output = BTFN<BothGenerators<G1, G2>, BT1, N>;
             #[inline]
-            fn add(self, other : &'b BTFN<F, G2, BT2, N>) -> Self::Output {
+            fn add(self, other : &'b BTFN<G2, BT2, N>) -> Self::Output {
                 $preprocess(self).add_another(other.generator.clone())
             }
         }
     }
 }
 
-make_btfn_add!(BTFN<F, G1, BT1, N>, std::convert::identity, );
-make_btfn_add!(&'a BTFN<F, G1, BT1, N>, Clone::clone, );
+make_btfn_add!(BTFN<G1, BT1, N>, std::convert::identity, );
+make_btfn_add!(&'a BTFN<G1, BT1, N>, Clone::clone, );
 
 macro_rules! make_btfn_sub {
     ($lhs:ty, $preprocess:path, $($extra_trait:ident)?) => {
         impl<'a, F : Float, G1, G2, BT1, BT2, const N : usize>
-        std::ops::Sub<BTFN<F, G2, BT2, N>> for
+        std::ops::Sub<BTFN<G2, BT2, N>> for
         $lhs
         where BT1 : BTImpl<F, N, Data=usize>,
-              G1 : SupportGenerator<F, N, Id=usize> + $($extra_trait)?,
-              G2 : SupportGenerator<F, N, Id=usize>,
-              G1::SupportType : LocalAnalysis<F, BT1::Agg, N>,
-              G2::SupportType : LocalAnalysis<F, BT1::Agg, N> {
-            type Output = BTFN<F, BothGenerators<G1, G2>, BT1, N>;
+              G1 : SupportGenerator<N, Id=usize, RealField = F> + $($extra_trait)?,
+              G2 : SupportGenerator<N, Id=usize, RealField = F>,
+              G1::SupportType : LocalAnalysis<BT1::Agg, Cube<F, N>>,
+              G2::SupportType : LocalAnalysis<BT1::Agg, Cube<F, N>> {
+            type Output = BTFN<BothGenerators<G1, G2>, BT1, N>;
             #[inline]
-            fn sub(self, other : BTFN<F, G2, BT2, N>) -> Self::Output {
+            fn sub(self, other : BTFN<G2, BT2, N>) -> Self::Output {
                 $preprocess(self).add_another(Arc::new(
                     Arc::try_unwrap(other.generator)
                         .unwrap_or_else(|arc| (*arc).clone())
@@ -240,35 +238,35 @@
         }
 
         impl<'a, 'b, F : Float, G1, G2,  BT1, BT2, const N : usize>
-        std::ops::Sub<&'b BTFN<F, G2, BT2, N>> for
+        std::ops::Sub<&'b BTFN<G2, BT2, N>> for
         $lhs
         where BT1 : BTImpl<F, N, Data=usize>,
-              G1 : SupportGenerator<F, N, Id=usize> + $($extra_trait)?,
-              G2 : SupportGenerator<F, N, Id=usize> + Clone,
-              G1::SupportType : LocalAnalysis<F, BT1::Agg, N>,
-              G2::SupportType : LocalAnalysis<F, BT1::Agg, N>,
+              G1 : SupportGenerator<N, Id=usize, RealField = F> + $($extra_trait)?,
+              G2 : SupportGenerator<N, Id=usize, RealField = F> + Clone,
+              G1::SupportType : LocalAnalysis<BT1::Agg, Cube<F, N>>,
+              G2::SupportType : LocalAnalysis<BT1::Agg, Cube<F, N>>,
               &'b G2 : std::ops::Neg<Output=G2> {
 
-            type Output = BTFN<F, BothGenerators<G1, G2>, BT1, N>;
+            type Output = BTFN<BothGenerators<G1, G2>, BT1, N>;
             #[inline]
-            fn sub(self, other : &'b BTFN<F, G2, BT2, N>) -> Self::Output {
+            fn sub(self, other : &'b BTFN<G2, BT2, N>) -> Self::Output {
                 $preprocess(self).add_another(Arc::new((*other.generator).clone().neg()))
             }
         }
     }
 }
 
-make_btfn_sub!(BTFN<F, G1, BT1, N>, std::convert::identity, );
-make_btfn_sub!(&'a BTFN<F, G1, BT1, N>, std::convert::identity, );
+make_btfn_sub!(BTFN<G1, BT1, N>, std::convert::identity, );
+make_btfn_sub!(&'a BTFN<G1, BT1, N>, std::convert::identity, );
 
 macro_rules! make_btfn_scalarop_rhs {
     ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
         impl<F : Float, G, BT, const N : usize>
         std::ops::$trait_assign<F>
-        for BTFN<F, G, BT, N>
+        for BTFN<G, BT, N>
         where BT : BTImpl<F, N>,
-              G : SupportGenerator<F, N, Id=BT::Data>,
-              G::SupportType : LocalAnalysis<F, BT::Agg, N> {
+              G : SupportGenerator<N, Id=BT::Data, RealField = F>,
+              G::SupportType : LocalAnalysis<BT::Agg, Cube<F, N>> {
             #[inline]
             fn $fn_assign(&mut self, t : F) {
                 Arc::make_mut(&mut self.generator).$fn_assign(t);
@@ -278,10 +276,10 @@
 
         impl<F : Float, G, BT, const N : usize>
         std::ops::$trait<F>
-        for BTFN<F, G, BT, N>
+        for BTFN<G, BT, N>
         where BT : BTImpl<F, N>,
-              G : SupportGenerator<F, N, Id=BT::Data>,
-              G::SupportType : LocalAnalysis<F, BT::Agg, N> {
+              G : SupportGenerator<N, Id=BT::Data, RealField = F>,
+              G::SupportType : LocalAnalysis<BT::Agg, Cube<F, N>> {
             type Output = Self;
             #[inline]
             fn $fn(mut self, t : F) -> Self::Output {
@@ -293,12 +291,12 @@
 
         impl<'a, F : Float, G, BT, const N : usize>
         std::ops::$trait<F>
-        for &'a BTFN<F, G, BT, N>
+        for &'a BTFN<G, BT, N>
         where BT : BTImpl<F, N>,
-              G : SupportGenerator<F, N, Id=BT::Data>,
-              G::SupportType : LocalAnalysis<F, BT::Agg, N>,
+              G : SupportGenerator<N, Id=BT::Data, RealField = F>,
+              G::SupportType : LocalAnalysis<BT::Agg, Cube<F, N>>,
               &'a G : std::ops::$trait<F,Output=G> {
-            type Output = BTFN<F, G, BT, N>;
+            type Output = BTFN<G, BT, N>;
             #[inline]
             fn $fn(self, t : F) -> Self::Output {
                 self.new_generator(self.generator.$fn(t))
@@ -313,14 +311,14 @@
 macro_rules! make_btfn_scalarop_lhs {
     ($trait:ident, $fn:ident, $fn_assign:ident, $($f:ident)+) => { $(
         impl<G, BT, const N : usize>
-        std::ops::$trait<BTFN<$f, G, BT, N>>
+        std::ops::$trait<BTFN<G, BT, N>>
         for $f
         where BT : BTImpl<$f, N>,
-              G : SupportGenerator<$f, N, Id=BT::Data>,
-              G::SupportType : LocalAnalysis<$f, BT::Agg, N> {
-            type Output = BTFN<$f, G, BT, N>;
+              G : SupportGenerator<N, Id=BT::Data, RealField = $f>,
+              G::SupportType : LocalAnalysis<BT::Agg, Cube<$f, N>> {
+            type Output = BTFN<G, BT, N>;
             #[inline]
-            fn $fn(self, mut a : BTFN<$f, G, BT, N>) -> Self::Output {
+            fn $fn(self, mut a : BTFN<G, BT, N>) -> Self::Output {
                 Arc::make_mut(&mut a.generator).$fn_assign(self);
                 a.refresh_aggregator();
                 a
@@ -328,16 +326,16 @@
         }
 
         impl<'a, G, BT, const N : usize>
-        std::ops::$trait<&'a BTFN<$f, G, BT, N>>
+        std::ops::$trait<&'a BTFN<G, BT, N>>
         for $f
         where BT : BTImpl<$f, N>,
-              G : SupportGenerator<$f, N, Id=BT::Data> + Clone,
-              G::SupportType : LocalAnalysis<$f, BT::Agg, N>,
+              G : SupportGenerator<N, Id=BT::Data, RealField = $f> + Clone,
+              G::SupportType : LocalAnalysis<BT::Agg, Cube<$f, N>>,
               // FIXME: This causes compiler overflow
               /*&'a G : std::ops::$trait<$f,Output=G>*/ {
-            type Output = BTFN<$f, G, BT, N>;
+            type Output = BTFN<G, BT, N>;
             #[inline]
-            fn $fn(self, a : &'a BTFN<$f, G, BT, N>) -> Self::Output {
+            fn $fn(self, a : &'a BTFN<G, BT, N>) -> Self::Output {
                 let mut tmp = (*a.generator).clone();
                 tmp.$fn_assign(self);
                 a.new_generator(tmp)
@@ -355,10 +353,10 @@
     ($trait:ident, $fn:ident) => {
         impl<F : Float, G, BT, const N : usize>
         std::ops::$trait
-        for BTFN<F, G, BT, N>
+        for BTFN<G, BT, N>
         where BT : BTImpl<F, N>,
-              G : SupportGenerator<F, N, Id=BT::Data>,
-              G::SupportType : LocalAnalysis<F, BT::Agg, N> {
+              G : SupportGenerator<N, Id=BT::Data, RealField = F>,
+              G::SupportType : LocalAnalysis<BT::Agg, Cube<F, N>> {
             type Output = Self;
             #[inline]
             fn $fn(mut self) -> Self::Output {
@@ -391,10 +389,11 @@
 //
 
 impl<'a, F : Float, G, BT, V, const N : usize> Apply<&'a Loc<F, N>>
-for BTFN<F, G, BT, N>
+for BTFN<G, BT, N>
 where BT : BTImpl<F, N>,
-      G : SupportGenerator<F, N, Id=BT::Data>,
-      G::SupportType : LocalAnalysis<F, BT::Agg, N> + Apply<&'a Loc<F, N>, Output = V>,
+      G : SupportGenerator<N, Id=BT::Data, RealField = F>,
+      G::SupportType : LocalAnalysis<BT::Agg, Cube<F, N>>
+                       + Apply<&'a Loc<F, N>, Output = V>,
       V : Sum {
 
     type Output = V;
@@ -406,10 +405,11 @@
 }
 
 impl<F : Float, G, BT, V, const N : usize> Apply<Loc<F, N>>
-for BTFN<F, G, BT, N>
+for BTFN<G, BT, N>
 where BT : BTImpl<F, N>,
-      G : SupportGenerator<F, N, Id=BT::Data>,
-      G::SupportType : LocalAnalysis<F, BT::Agg, N> + Apply<Loc<F, N>, Output = V>,
+      G : SupportGenerator<N, Id=BT::Data, RealField = F>,
+      G::SupportType : LocalAnalysis<BT::Agg, Cube<F, N>>
+                       + Apply<Loc<F, N>, Output = V>,
       V : Sum {
 
     type Output = V;
@@ -420,11 +420,13 @@
     }
 }
 
-impl<'a, F : Float, G, BT, V, const N : usize> Differentiable<&'a Loc<F, N>>
-for BTFN<F, G, BT, N>
+impl<'a, F : Float, G, BT, V : HasRealField<RealField=F>, const N : usize>
+Differentiable<&'a Loc<F, N>>
+for BTFN<G, BT, N>
 where BT : BTImpl<F, N>,
-      G : SupportGenerator<F, N, Id=BT::Data>,
-      G::SupportType : LocalAnalysis<F, BT::Agg, N> + Differentiable<&'a Loc<F, N>, Derivative = V>,
+      G : SupportGenerator<N, Id=BT::Data, RealField = F>,
+      G::SupportType : LocalAnalysis<BT::Agg, Cube<F, N>>
+                       + Differentiable<&'a Loc<F, N>, Derivative = V>,
       V : Sum {
 
     type Derivative = V;
@@ -437,11 +439,12 @@
 }
 
 impl<F : Float, G, BT, V, const N : usize> Differentiable<Loc<F, N>>
-for BTFN<F, G, BT, N>
+for BTFN<G, BT, N>
 where BT : BTImpl<F, N>,
-      G : SupportGenerator<F, N, Id=BT::Data>,
-      G::SupportType : LocalAnalysis<F, BT::Agg, N> + Differentiable<Loc<F, N>, Derivative = V>,
-      V : Sum {
+      G : SupportGenerator<N, Id=BT::Data, RealField = F>,
+      G::SupportType : LocalAnalysis<BT::Agg, Cube<F, N>>
+                       + Differentiable<Loc<F, N>, Derivative = V>,
+      V : Sum + HasRealField<RealField=F> {
 
     type Derivative = V;
 
@@ -456,11 +459,11 @@
 // GlobalAnalysis
 //
 
-impl<F : Float, G, BT, const N : usize> GlobalAnalysis<F, BT::Agg>
-for BTFN<F, G, BT, N>
+impl<F : Float, G, BT, const N : usize> GlobalAnalysis<BT::Agg>
+for BTFN<G, BT, N>
 where BT : BTImpl<F, N>,
-      G : SupportGenerator<F, N, Id=BT::Data>,
-      G::SupportType : LocalAnalysis<F, BT::Agg, N> {
+      G : SupportGenerator<N, Id=BT::Data, RealField = F>,
+      G::SupportType : LocalAnalysis<BT::Agg, Cube<F, N>> {
 
     #[inline]
     fn global_analysis(&self) -> BT::Agg {
@@ -617,12 +620,12 @@
     how : T,
 }
 
-impl<F : Float, G, const N : usize> Refiner<F, Bounds<F>, G, N>
+impl<F : Float, G, const N : usize> Refiner<Bounds<F>, G, N>
 for P2Refiner<F, RefineMax>
 where Cube<F, N> : P2Minimise<Loc<F, N>, F>,
-      G : SupportGenerator<F, N>,
+      G : SupportGenerator<N, RealField = F>,
       G::SupportType : Mapping<Loc<F, N>, Codomain=F>
-                       + LocalAnalysis<F, Bounds<F>, N> {
+                       + LocalAnalysis<Bounds<F>, Cube<F, N>> {
     type Result = Option<(Loc<F, N>, F)>;
     type Sorting = UpperBoundSorting<F>;
 
@@ -674,12 +677,12 @@
 }
 
 
-impl<F : Float, G, const N : usize> Refiner<F, Bounds<F>, G, N>
+impl<F : Float, G, const N : usize> Refiner<Bounds<F>, G, N>
 for P2Refiner<F, RefineMin>
 where Cube<F, N> : P2Minimise<Loc<F, N>, F>,
-      G : SupportGenerator<F, N>,
+      G : SupportGenerator<N, RealField = F>,
       G::SupportType : Mapping<Loc<F, N>, Codomain=F>
-                       + LocalAnalysis<F, Bounds<F>, N> {
+                       + LocalAnalysis<Bounds<F>, Cube<F, N>> {
     type Result = Option<(Loc<F, N>, F)>;
     type Sorting = LowerBoundSorting<F>;
 
@@ -754,9 +757,9 @@
     how : T,
 }
 
-impl<F : Float, G, const N : usize> Refiner<F, Bounds<F>, G, N>
+impl<F : Float, G, const N : usize> Refiner<Bounds<F>, G, N>
 for BoundRefiner<F, RefineMax>
-where G : SupportGenerator<F, N> {
+where G : SupportGenerator<N, RealField = F> {
     type Result = bool;
     type Sorting = UpperBoundSorting<F>;
 
@@ -788,9 +791,9 @@
     }
 }
 
-impl<F : Float, G, const N : usize> Refiner<F, Bounds<F>, G, N>
+impl<F : Float, G, const N : usize> Refiner<Bounds<F>, G, N>
 for BoundRefiner<F, RefineMin>
-where G : SupportGenerator<F, N> {
+where G : SupportGenerator<N, RealField = F> {
     type Result = bool;
     type Sorting = UpperBoundSorting<F>;
 
@@ -837,11 +840,11 @@
 // there should be a result, or new nodes above the `glb` inserted into the queue. Then the waiting
 // threads can also continue processing. If, however, numerical inaccuracy destroyes the `glb`,
 // the queue may run out, and we get “Refiner failure”.
-impl<F : Float, G, BT, const N : usize> BTFN<F, G, BT, N>
+impl<F : Float, G, BT, const N : usize> BTFN<G, BT, N>
 where BT : BTSearch<F, N, Agg=Bounds<F>>,
-      G : SupportGenerator<F, N, Id=BT::Data>,
+      G : SupportGenerator<N, Id=BT::Data, RealField = F>,
       G::SupportType : Mapping<Loc<F, N>,Codomain=F>
-                       + LocalAnalysis<F, Bounds<F>, N>,
+                       + LocalAnalysis<Bounds<F>, Cube<F, N>>,
       Cube<F, N> : P2Minimise<Loc<F, N>, F> {
 
     /// Maximise the `BTFN` within stated value `tolerance`.
diff -r edb95d2b83cc -r d2acaaddd9af src/bisection_tree/either.rs
--- a/src/bisection_tree/either.rs	Sun Nov 10 09:02:57 2024 -0500
+++ b/src/bisection_tree/either.rs	Tue Dec 31 09:12:43 2024 -0500
@@ -29,12 +29,18 @@
     Right(B),
 }
 
+impl<A, B, F : Num> HasScalarField for EitherSupport<A, B>
+where A : HasScalarField<Field = F>,
+      B : HasScalarField<Field = F> {
+    type Field = F;
+}
+
 // We need type alias bounds to access associate types.
 #[allow(type_alias_bounds)]
 type BothAllDataIter<
-    'a, F,
-    G1 : SupportGenerator<F, N>,
-    G2 : SupportGenerator<F, N>,
+    'a,
+    G1 : SupportGenerator<N>,
+    G2 : SupportGenerator<N, RealField = G1::RealField>,
     const N : usize
 > = Chain<
     MapF<G1::AllDataIter<'a>, (usize, EitherSupport<G1::SupportType, G2::SupportType>)>,
@@ -46,8 +52,8 @@
     #[inline]
     fn map_left<F : Float, const N : usize>((d, support) : (G1::Id, G1::SupportType))
     -> (usize, EitherSupport<G1::SupportType, G2::SupportType>)
-    where G1 : SupportGenerator<F, N, Id=usize>,
-          G2 : SupportGenerator<F, N, Id=usize> {
+    where G1 : SupportGenerator<N, Id=usize, RealField =F>,
+          G2 : SupportGenerator<N, Id=usize, RealField = F> {
 
         let id : usize = d.into();
         (id.into(), EitherSupport::Left(support))
@@ -57,8 +63,8 @@
     #[inline]
     fn map_right<F : Float, const N : usize>(n0 : &usize, (d, support) : (G2::Id, G2::SupportType))
     -> (usize, EitherSupport<G1::SupportType, G2::SupportType>)
-    where G1 : SupportGenerator<F, N, Id=usize>,
-          G2 : SupportGenerator<F, N, Id=usize> {
+    where G1 : SupportGenerator<N, Id=usize, RealField = F>,
+          G2 : SupportGenerator<N, Id=usize, RealField = F> {
 
         let id : usize = d.into();
         ((n0+id).into(), EitherSupport::Right(support))
@@ -70,8 +76,8 @@
     #[inline]
     pub(super) fn all_left_data<F : Float, const N : usize>(&self)
     -> MapF<G1::AllDataIter<'_>, (usize, EitherSupport<G1::SupportType, G2::SupportType>)>
-    where G1 : SupportGenerator<F, N, Id=usize>,
-          G2 : SupportGenerator<F, N, Id=usize> {
+    where G1 : SupportGenerator<N, Id=usize, RealField = F>,
+          G2 : SupportGenerator<N, Id=usize, RealField = F> {
         self.0.all_data().mapF(Self::map_left)
     }
 
@@ -81,22 +87,29 @@
     #[inline]
     pub(super) fn all_right_data<F : Float, const N : usize>(&self)
     -> MapZ<G2::AllDataIter<'_>, usize, (usize, EitherSupport<G1::SupportType, G2::SupportType>)>
-    where G1 : SupportGenerator<F, N, Id=usize>,
-          G2 : SupportGenerator<F, N, Id=usize> {
+    where G1 : SupportGenerator<N, Id=usize, RealField=F>,
+          G2 : SupportGenerator<N, Id=usize, RealField=F> {
         let n0 = self.0.support_count();
         self.1.all_data().mapZ(n0, Self::map_right)
     }
 }
 
+impl<F : Float, G1, G2> HasScalarField for BothGenerators<G1, G2>
+where G1 : HasScalarField<Field=F>,
+      G2 : HasScalarField<Field=F> {
+    type Field = F;
+}
+
 impl<F : Float, G1, G2, const N : usize>
-SupportGenerator<F, N>
+SupportGenerator<N>
 for BothGenerators<G1, G2>
-where G1 : SupportGenerator<F, N, Id=usize>,
-      G2 : SupportGenerator<F, N, Id=usize> {
+where G1 : SupportGenerator<N, Id=usize> + HasRealField<RealField=F>,
+      G2 : SupportGenerator<N, Id=usize> + HasRealField<RealField=F>,
+      Self : HasRealField<RealField = F> {
 
     type Id = usize;
     type SupportType = EitherSupport<G1::SupportType, G2::SupportType>;
-    type AllDataIter<'a> = BothAllDataIter<'a, F, G1, G2, N> where G1 : 'a, G2 : 'a;
+    type AllDataIter<'a> = BothAllDataIter<'a, G1, G2, N> where G1 : 'a, G2 : 'a;
 
     #[inline]
     fn support_for(&self,  id : Self::Id)
@@ -120,9 +133,9 @@
     }
 }
 
-impl<F: Float, S1, S2, const N : usize> Support<F, N> for EitherSupport<S1, S2>
-where S1 : Support<F, N>,
-      S2 : Support<F, N> {
+impl<F: Float, S1, S2, const N : usize> Support<N> for EitherSupport<S1, S2>
+where S1 : Support<N, RealField = F>,
+      S2 : Support<N, RealField = F> {
 
     #[inline]
     fn support_hint(&self) -> Cube<F,N> {
@@ -149,10 +162,11 @@
     }
 }
 
-impl<F : Float, A, S1, S2, const N : usize> LocalAnalysis<F, A, N> for EitherSupport<S1, S2>
+impl<F : Float, A, S1, S2, const N : usize> LocalAnalysis<A, Cube<F, N>>
+for EitherSupport<S1, S2>
 where A : Aggregator,
-      S1 : LocalAnalysis<F, A, N>,
-      S2 : LocalAnalysis<F, A, N>, {
+      S1 : LocalAnalysis<A, Cube<F, N>>,
+      S2 : LocalAnalysis<A, Cube<F, N>> {
 
     #[inline]
     fn local_analysis(&self, cube : &Cube<F, N>) -> A {
@@ -163,10 +177,10 @@
     }
 }
 
-impl<F : Float, A, S1, S2> GlobalAnalysis<F, A> for EitherSupport<S1, S2>
+impl<A, S1, S2> GlobalAnalysis<A> for EitherSupport<S1, S2>
 where A : Aggregator,
-      S1 : GlobalAnalysis<F, A>,
-      S2 : GlobalAnalysis<F, A>, {
+      S1 : GlobalAnalysis<A>,
+      S2 : GlobalAnalysis<A>, {
 
     #[inline]
     fn global_analysis(&self) -> A {
@@ -190,9 +204,9 @@
     }
 }
 
-impl<F, S1, S2, X> Differentiable<X> for EitherSupport<S1, S2>
-where S1 : Differentiable<X, Derivative=F>,
-      S2 : Differentiable<X, Derivative=F> {
+impl<F : Float, S1, S2, X> Differentiable<X> for EitherSupport<S1, S2>
+where S1 : Differentiable<X, Derivative=F, RealField=F>,
+      S2 : Differentiable<X, Derivative=F, RealField=F> {
     type Derivative = F;
     #[inline]
     fn differential(&self, x : X) -> F {
diff -r edb95d2b83cc -r d2acaaddd9af src/bisection_tree/refine.rs
--- a/src/bisection_tree/refine.rs	Sun Nov 10 09:02:57 2024 -0500
+++ b/src/bisection_tree/refine.rs	Tue Dec 31 09:12:43 2024 -0500
@@ -92,10 +92,9 @@
 /// The `Refiner` is used to determine whether an [`Aggregator`] `A` stored in the [`BT`] is
 /// sufficiently refined within a [`Cube`], and in such a case, produce a desired result (e.g.
 /// a maximum value of a function).
-pub trait Refiner<F : Float, A, G, const N : usize> : Sync + Send + 'static
-where F : Num,
-      A : Aggregator,
-      G : SupportGenerator<F, N> {
+pub trait Refiner<A, G, const N : usize> : Sync + Send + 'static
+where A : Aggregator,
+      G : SupportGenerator<N> {
 
     /// The result type of the refiner
     type Result : std::fmt::Debug + Sync + Send + 'static;
@@ -125,7 +124,7 @@
     fn refine(
         &self,
         aggregator : &A,
-        domain : &Cube<F, N>,
+        domain : &Cube<G::Field, N>,
         data : &[G::Id],
         generator : &G,
         step : usize,
@@ -323,9 +322,9 @@
         container_arc : &'c Arc<Mutex<HeapContainer<'a, F, D, A, R::Sorting, R::Result, N, P>>>,
         step : usize
     ) -> Result<R::Result, MutexGuard<'c, HeapContainer<'a, F, D, A, R::Sorting, R::Result, N, P>>>
-    where R : Refiner<F, A, G, N>,
-          G : SupportGenerator<F, N, Id=D>,
-          G::SupportType : LocalAnalysis<F, A, N> {
+    where R : Refiner<A, G, N>,
+          G : SupportGenerator<N, Id=D, RealField = F>,
+          G::SupportType : LocalAnalysis<A, Cube<F, N>> {
 
         //drop(container);
 
@@ -435,9 +434,9 @@
         refiner : R,
         generator : &Arc<G>,
     ) -> Option<R::Result>
-    where R : Refiner<F, Self::Agg, G, N> + Sync + Send + 'static,
-          G : SupportGenerator<F, N, Id=Self::Data> + Sync + Send + 'static,
-          G::SupportType : LocalAnalysis<F, Self::Agg, N>;
+    where R : Refiner<Self::Agg, G, N> + Sync + Send + 'static,
+          G : SupportGenerator<N, Id=Self::Data, RealField = F> + Sync + Send + 'static,
+          G::SupportType : LocalAnalysis<Self::Agg, Cube<F, N>>;
 }
 
 fn refinement_loop<F : Float, D, A, R, G, const N : usize, const P : usize> (
@@ -446,9 +445,9 @@
     generator_arc : &Arc<G>,
     container_arc : &Arc<Mutex<HeapContainer<F, D, A, R::Sorting, R::Result, N, P>>>,
 ) where A : Aggregator,
-        R : Refiner<F, A, G, N>,
-        G : SupportGenerator<F, N, Id=D>,
-        G::SupportType : LocalAnalysis<F, A, N>,
+        R : Refiner<A, G, N>,
+        G : SupportGenerator<N, Id=D, RealField = F>,
+        G::SupportType : LocalAnalysis<A, Cube<F, N>>,
         Const<P> : BranchCount<N>,
         D : 'static + Copy + Sync + Send + std::fmt::Debug {
 
@@ -565,9 +564,9 @@
                 refiner : R,
                 generator : &Arc<G>,
             ) -> Option<R::Result>
-            where R : Refiner<F, A, G, $n>,
-                  G : SupportGenerator<F, $n, Id=D>,
-                  G::SupportType : LocalAnalysis<F, A, $n> {
+            where R : Refiner<A, G, $n>,
+                  G : SupportGenerator<$n, Id=D, RealField = F>,
+                  G::SupportType : LocalAnalysis<A, Cube<F, $n>> {
                 let mut init_container = HeapContainer {
                     heap : BinaryHeap::new(),
                     glb : R::Sorting::bottom(),
diff -r edb95d2b83cc -r d2acaaddd9af src/bisection_tree/support.rs
--- a/src/bisection_tree/support.rs	Sun Nov 10 09:02:57 2024 -0500
+++ b/src/bisection_tree/support.rs	Tue Dec 31 09:12:43 2024 -0500
@@ -10,18 +10,17 @@
 use crate::sets::Cube;
 use crate::loc::Loc;
 use super::aggregator::Bounds;
-use crate::norms::{Norm, L1, L2, Linfinity};
+use crate::norms::{Norm, L1, L2, Linfinity, HasScalarField, HasRealField};
 
 /// A trait for encoding constant [`Float`] values
-pub trait Constant : Copy + Sync + Send + 'static + std::fmt::Debug + Into<Self::Type> {
-    /// The type of the value
-    type Type : Float;
+pub trait Constant : Copy + Sync + Send + 'static
+                     + std::fmt::Debug + HasRealField
+                     + Into<Self::Field> {
     /// Returns the value of the constant
-    fn value(&self) -> Self::Type;
+    fn value(&self) -> Self::Field;
 }
 
 impl<F : Float> Constant for F {
-    type Type = F;
     #[inline]
     fn value(&self) -> F { *self }
 }
@@ -30,17 +29,17 @@
 /// A trait for working with the supports of [`Apply`]s.
 ///
 /// Apply is not a super-trait to allow more general use.
-pub trait Support<F : Num, const N : usize> : Sized + Sync + Send + 'static {
+pub trait Support<const N : usize> : Sized + Sync + Send + 'static + HasRealField {
     /// Return a cube containing the support of the function represented by `self`.
     ///
     /// The hint may be larger than the actual support, but must contain it.
-    fn support_hint(&self) -> Cube<F,N>;
+    fn support_hint(&self) -> Cube<Self::Field,N>;
 
     /// Indicate whether `x` is in the support of the function represented by `self`.
-    fn in_support(&self, x : &Loc<F,N>) -> bool;
+    fn in_support(&self, x : &Loc<Self::Field,N>) -> bool;
 
     // Indicate whether `cube` is fully in the support of the function represented by `self`.
-    //fn fully_in_support(&self, cube : &Cube<F,N>) -> bool;
+    //fn fully_in_support(&self, cube : &Cube<Self::Field,N>) -> bool;
 
     /// Return an optional hint for bisecting the support.
     ///
@@ -53,19 +52,19 @@
     /// The default implementation returns `[None; N]`.
     #[inline]
     #[allow(unused_variables)]
-    fn bisection_hint(&self, cube : &Cube<F, N>) -> [Option<F>; N] {
+    fn bisection_hint(&self, cube : &Cube<Self::Field, N>) -> [Option<Self::Field>; N] {
         [None; N]
     }
 
     /// Translate `self` by `x`.
     #[inline]
-    fn shift(self, x : Loc<F, N>) -> Shift<Self, F, N> {
+    fn shift(self, x : Loc<Self::Field, N>) -> Shift<Self, Self::Field, N> {
         Shift { shift : x, base_fn : self }
     }
 
     /// Multiply `self` by the scalar `a`.
     #[inline]
-    fn weigh<C : Constant<Type=F>>(self, a : C) -> Weighted<Self, C> {
+    fn weigh<C : Constant<RealField=Self::RealField>>(self, a : C) -> Weighted<Self, C> {
         Weighted { weight : a, base_fn : self }
     }
 }
@@ -74,7 +73,7 @@
 ///
 /// Typically `A` is an [`Aggregator`][super::aggregator::Aggregator] such as
 /// [`Bounds`][super::aggregator::Bounds].
-pub trait GlobalAnalysis<F : Num, A> {
+pub trait GlobalAnalysis<A> {
     /// Perform global analysis of the property `A` of `Self`.
     ///
     /// As an example, in the case of `A` being [`Bounds`][super::aggregator::Bounds],
@@ -83,26 +82,18 @@
     fn global_analysis(&self) -> A;
 }
 
-// default impl<F, A, N, L> GlobalAnalysis<F, A, N> for L
-// where L : LocalAnalysis<F, A, N> {
-//     #[inline]
-//     fn global_analysis(&self) -> Bounds<F> {
-//         self.local_analysis(&self.support_hint())
-//     }
-// }
-
 /// Trait for locally analysing a property `A` of a [`Apply`] (implementing [`Support`])
 /// within a [`Cube`].
 ///
 /// Typically `A` is an [`Aggregator`][super::aggregator::Aggregator] such as
 /// [`Bounds`][super::aggregator::Bounds].
-pub trait LocalAnalysis<F : Num, A, const N : usize> : GlobalAnalysis<F, A> + Support<F, N> {
+pub trait LocalAnalysis<A, LocalSet> : GlobalAnalysis<A> {
     /// Perform local analysis of the property `A` of `Self`.
     ///
     /// As an example, in the case of `A` being [`Bounds`][super::aggregator::Bounds],
     /// this function will return upper and lower bounds within `cube` for the mapping
     /// represented by `self`.
-    fn local_analysis(&self, cube : &Cube<F, N>) -> A;
+    fn local_analysis(&self, cube : &LocalSet) -> A;
 }
 
 /// Trait for determining the upper and lower bounds of an float-valued [`Apply`].
@@ -110,15 +101,15 @@
 /// This is a blanket-implemented alias for [`GlobalAnalysis`]`<F, Bounds<F>>`
 /// [`Apply`] is not a supertrait to allow flexibility in the implementation of either
 /// reference or non-reference arguments.
-pub trait Bounded<F : Float> : GlobalAnalysis<F, Bounds<F>> {
+pub trait Bounded : GlobalAnalysis<Bounds<Self::Field>> + HasScalarField {
     /// Return lower and upper bounds for the values of of `self`.
     #[inline]
-    fn bounds(&self) -> Bounds<F> {
+    fn bounds(&self) -> Bounds<Self::Field> {
         self.global_analysis()
     }
 }
 
-impl<F : Float, T : GlobalAnalysis<F, Bounds<F>>> Bounded<F> for T { }
+impl<T : GlobalAnalysis<Bounds<Self::Field>> + HasScalarField> Bounded for T { }
 
 /// Shift of [`Support`] and [`Apply`]; output of [`Support::shift`].
 #[derive(Copy,Clone,Debug,Serialize)] // Serialize! but not implemented by Loc.
@@ -127,6 +118,10 @@
     base_fn : T,
 }
 
+impl<T, F : Num, const N : usize>  HasScalarField for Shift<T, F, N> {
+    type Field = F;
+}
+
 impl<'a, T, V, F : Float, const N : usize> Apply<&'a Loc<F, N>> for Shift<T,F,N>
 where T : Apply<Loc<F, N>, Output=V> {
     type Output = V;
@@ -145,7 +140,8 @@
     }
 }
 
-impl<'a, T, V, F : Float, const N : usize> Differentiable<&'a Loc<F, N>> for Shift<T,F,N>
+impl<'a, T, V : HasRealField<RealField=F>, F : Float, const N : usize>
+Differentiable<&'a Loc<F, N>> for Shift<T,F,N>
 where T : Differentiable<Loc<F, N>, Derivative=V> {
     type Derivative = V;
     #[inline]
@@ -154,7 +150,8 @@
     }
 }
 
-impl<'a, T, V, F : Float, const N : usize> Differentiable<Loc<F, N>> for Shift<T,F,N>
+impl<'a, T, V : HasRealField<RealField=F>, F : Float, const N : usize>
+Differentiable<Loc<F, N>> for Shift<T,F,N>
 where T : Differentiable<Loc<F, N>, Derivative=V> {
     type Derivative = V;
     #[inline]
@@ -163,8 +160,8 @@
     }
 }
 
-impl<'a, T, F : Float, const N : usize> Support<F,N> for Shift<T,F,N>
-where T : Support<F, N> {
+impl<'a, T, F : Float, const N : usize> Support<N> for Shift<T,F,N>
+where T : Support<N, RealField = F> {
     #[inline]
     fn support_hint(&self) -> Cube<F,N> {
         self.base_fn.support_hint().shift(&self.shift)
@@ -188,16 +185,17 @@
 
 }
 
-impl<'a, T, F : Float, const N : usize> GlobalAnalysis<F, Bounds<F>> for Shift<T,F,N>
-where T : LocalAnalysis<F, Bounds<F>, N> {
+impl<'a, T, F : Float, const N : usize> GlobalAnalysis<Bounds<F>> for Shift<T,F,N>
+where T : GlobalAnalysis<Bounds<F>> {
     #[inline]
     fn global_analysis(&self) -> Bounds<F> {
         self.base_fn.global_analysis()
     }
 }
 
-impl<'a, T, F : Float, const N : usize> LocalAnalysis<F, Bounds<F>, N> for Shift<T,F,N>
-where T : LocalAnalysis<F, Bounds<F>, N> {
+impl<'a, T, F : Float, const N : usize> LocalAnalysis<Bounds<F>, Cube<F, N>>
+for Shift<T,F,N>
+where T : LocalAnalysis<Bounds<F>, Cube<F, N>> {
     #[inline]
     fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> {
         self.base_fn.local_analysis(&cube.shift(&(-self.shift)))
@@ -206,10 +204,10 @@
 
 macro_rules! impl_shift_norm {
     ($($norm:ident)*) => { $(
-        impl<'a, T, F : Float, const N : usize> Norm<F, $norm> for Shift<T,F,N>
-        where T : Norm<F, $norm> {
+        impl<'a, T, const N : usize> Norm<$norm> for Shift<T,T::Field,N>
+        where T : Norm<$norm> {
             #[inline]
-            fn norm(&self, n : $norm) -> F {
+            fn norm(&self, n : $norm) -> T::Field {
                 self.base_fn.norm(n)
             }
         }
@@ -228,10 +226,16 @@
     pub base_fn : T,
 }
 
+impl<T, C, F : Float> HasScalarField for Weighted<T, C>
+where T : HasRealField<RealField = F>,
+      C : Constant<RealField = F> {
+    type Field = F;
+}
+
 impl<'a, T, V, F : Float, C, const N : usize> Apply<&'a Loc<F, N>> for Weighted<T, C>
 where T : for<'b> Apply<&'b Loc<F, N>, Output=V>,
       V : std::ops::Mul<F,Output=V>,
-      C : Constant<Type=F> {
+      C : Constant<RealField=F> {
     type Output = V;
     #[inline]
     fn apply(&self, x : &'a Loc<F, N>) -> Self::Output {
@@ -242,7 +246,7 @@
 impl<'a, T, V, F : Float, C, const N : usize> Apply<Loc<F, N>> for Weighted<T, C>
 where T : Apply<Loc<F, N>, Output=V>,
       V : std::ops::Mul<F,Output=V>,
-      C : Constant<Type=F> {
+      C : Constant<RealField=F> {
     type Output = V;
     #[inline]
     fn apply(&self, x : Loc<F, N>) -> Self::Output {
@@ -250,10 +254,11 @@
     }
 }
 
-impl<'a, T, V, F : Float, C, const N : usize> Differentiable<&'a Loc<F, N>> for Weighted<T, C>
-where T : for<'b> Differentiable<&'b Loc<F, N>, Derivative=V>,
+impl<'a, T, V : HasRealField<RealField=F>, F : Float, C, const N : usize>
+Differentiable<&'a Loc<F, N>> for Weighted<T, C>
+where T : for<'b> Differentiable<&'b Loc<F, N>, Derivative=V, RealField=F>,
       V : std::ops::Mul<F, Output=V>,
-      C : Constant<Type=F> {
+      C : Constant<RealField=F> {
     type Derivative = V;
     #[inline]
     fn differential(&self, x : &'a Loc<F, N>) -> Self::Derivative {
@@ -261,10 +266,11 @@
     }
 }
 
-impl<'a, T, V, F : Float, C, const N : usize> Differentiable<Loc<F, N>> for Weighted<T, C>
-where T : Differentiable<Loc<F, N>, Derivative=V>,
+impl<'a, T, V : HasRealField<RealField = F>, F : Float, C, const N : usize>
+Differentiable<Loc<F, N>> for Weighted<T, C>
+where T : Differentiable<Loc<F, N>, Derivative=V, RealField = F>,
       V : std::ops::Mul<F, Output=V>,
-      C : Constant<Type=F> {
+      C : Constant<RealField=F> {
     type Derivative = V;
     #[inline]
     fn differential(&self, x : Loc<F, N>) -> Self::Derivative {
@@ -272,9 +278,9 @@
     }
 }
 
-impl<'a, T, F : Float, C, const N : usize> Support<F,N> for Weighted<T, C>
-where T : Support<F, N>,
-      C : Constant<Type=F> {
+impl<'a, T, F : Float, C, const N : usize> Support<N> for Weighted<T, C>
+where T : Support<N, RealField = F>,
+      C : Constant<RealField=F> {
 
     #[inline]
     fn support_hint(&self) -> Cube<F,N> {
@@ -296,9 +302,9 @@
     }
 }
 
-impl<'a, T, F : Float, C> GlobalAnalysis<F, Bounds<F>> for Weighted<T, C>
-where T : GlobalAnalysis<F, Bounds<F>>,
-      C : Constant<Type=F> {
+impl<'a, T, F : Float, C> GlobalAnalysis<Bounds<F>> for Weighted<T, C>
+where T : GlobalAnalysis<Bounds<F>>,
+      C : Constant<RealField=F> {
     #[inline]
     fn global_analysis(&self) -> Bounds<F> {
         let Bounds(lower, upper) = self.base_fn.global_analysis();
@@ -310,9 +316,10 @@
     }
 }
 
-impl<'a, T, F : Float, C, const N : usize> LocalAnalysis<F, Bounds<F>, N> for Weighted<T, C>
-where T : LocalAnalysis<F, Bounds<F>, N>,
-      C : Constant<Type=F> {
+impl<'a, T, F : Float, C, const N : usize> LocalAnalysis<Bounds<F>, Cube<F, N>>
+for Weighted<T, C>
+where T : LocalAnalysis<Bounds<F>, Cube<F, N>>,
+      C : Constant<RealField=F> {
     #[inline]
     fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> {
         let Bounds(lower, upper) = self.base_fn.local_analysis(cube);
@@ -358,8 +365,9 @@
 
 macro_rules! impl_weighted_norm {
     ($($norm:ident)*) => { $(
-        impl<'a, T, F : Float> Norm<F, $norm> for Weighted<T,F>
-        where T : Norm<F, $norm> {
+        impl<'a, T, F : Float> Norm<$norm> for Weighted<T,F>
+        where T : Norm<$norm, Field = F> {
+
             #[inline]
             fn norm(&self, n : $norm) -> F {
                 self.base_fn.norm(n) * self.weight.abs()
@@ -381,7 +389,7 @@
 );
 
 impl<'a, T, F : Float, const N : usize> Apply<&'a Loc<F, N>> for Normalised<T>
-where T : Norm<F, L1> + for<'b> Apply<&'b Loc<F, N>, Output=F> {
+where T : Norm<L1, Field=F> + for<'b> Apply<&'b Loc<F, N>, Output=F> {
     type Output = F;
     #[inline]
     fn apply(&self, x : &'a Loc<F, N>) -> Self::Output {
@@ -391,7 +399,7 @@
 }
 
 impl<'a, T, F : Float, const N : usize> Apply<Loc<F, N>> for Normalised<T>
-where T : Norm<F, L1> + Apply<Loc<F,N>, Output=F> {
+where T : Norm<L1, Field = F> + Apply<Loc<F,N>, Output=F> {
     type Output = F;
     #[inline]
     fn apply(&self, x : Loc<F, N>) -> Self::Output {
@@ -400,8 +408,8 @@
     }
 }
 
-impl<'a, T, F : Float, const N : usize> Support<F,N> for Normalised<T>
-where T : Norm<F, L1> + Support<F, N> {
+impl<'a, T, F : Float, const N : usize> Support<N> for Normalised<T>
+where T : Norm<L1> + Support<N, RealField = F> {
 
     #[inline]
     fn support_hint(&self) -> Cube<F,N> {
@@ -423,8 +431,8 @@
     }
 }
 
-impl<'a, T, F : Float> GlobalAnalysis<F, Bounds<F>> for Normalised<T>
-where T : Norm<F, L1> + GlobalAnalysis<F, Bounds<F>> {
+impl<'a, T, F : Float> GlobalAnalysis<Bounds<F>> for Normalised<T>
+where T : Norm<L1, Field = F> + GlobalAnalysis<Bounds<F>> {
     #[inline]
     fn global_analysis(&self) -> Bounds<F> {
         let Bounds(lower, upper) = self.0.global_analysis();
@@ -435,8 +443,9 @@
     }
 }
 
-impl<'a, T, F : Float, const N : usize> LocalAnalysis<F, Bounds<F>, N> for Normalised<T>
-where T : Norm<F, L1> + LocalAnalysis<F, Bounds<F>, N> {
+impl<'a, T, F : Float, const N : usize> LocalAnalysis<Bounds<F>, Cube<F, N>>
+for Normalised<T>
+where T : Norm<L1, Field = F> + LocalAnalysis<Bounds<F>, Cube<F, N>> {
     #[inline]
     fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> {
         let Bounds(lower, upper) = self.0.local_analysis(cube);
@@ -447,23 +456,26 @@
     }
 }
 
-impl<'a, T, F : Float> Norm<F, L1> for Normalised<T>
-where T : Norm<F, L1> {
+impl<'a, T> HasScalarField for Normalised<T> where T : HasScalarField {
+    type Field = T::Field;
+}
+
+impl<'a, T> Norm<L1> for Normalised<T> where T : Norm<L1> {
     #[inline]
-    fn norm(&self, _ : L1) -> F {
+    fn norm(&self, _ : L1) -> Self::Field {
         let w = self.0.norm(L1);
-        if w == F::ZERO { F::ZERO } else { F::ONE }
+        if w == Self::Field::ZERO { Self::Field::ZERO } else { Self::Field::ONE }
     }
 }
 
 macro_rules! impl_normalised_norm {
     ($($norm:ident)*) => { $(
-        impl<'a, T, F : Float> Norm<F, $norm> for Normalised<T>
-        where T : Norm<F, $norm> + Norm<F, L1> {
+        impl<'a, T> Norm<$norm> for Normalised<T> where T : Norm<$norm> + Norm<L1> {
+
             #[inline]
-            fn norm(&self, n : $norm) -> F {
+            fn norm(&self, n : $norm) -> Self::Field {
                 let w = self.0.norm(L1);
-                if w == F::ZERO { F::ZERO } else { self.0.norm(n) / w }
+                if w == Self::Field::ZERO { Self::Field::ZERO } else { self.0.norm(n) / w }
             }
         }
     )* }
@@ -471,26 +483,16 @@
 
 impl_normalised_norm!(L2 Linfinity);
 
-/*
-impl<F : Num, S : Support<F, N>, const N : usize> LocalAnalysis<F, NullAggregator, N> for S {
-    fn local_analysis(&self, _cube : &Cube<F, N>) -> NullAggregator { NullAggregator }
-}
-
-impl<F : Float, S : Bounded<F>, const N : usize> LocalAnalysis<F, Bounds<F>, N> for S {
-    #[inline]
-    fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> {
-        self.bounds(cube)
-    }
-}*/
 
 /// Generator of [`Support`]-implementing component functions based on low storage requirement
 /// [ids][`Self::Id`].
-pub trait SupportGenerator<F : Float, const N : usize>
-: MulAssign<F> + DivAssign<F> + Neg<Output=Self> + Clone + Sync + Send + 'static {
+pub trait SupportGenerator<const N : usize>
+: MulAssign<Self::Field> + DivAssign<Self::Field> + Neg<Output=Self> + HasRealField
+  + Clone + Sync + Send + 'static {
     /// The identification type
     type Id : 'static + Copy;
     /// The type of the [`Support`] (often also a [`Apply`]).
-    type SupportType : 'static + Support<F, N>;
+    type SupportType : 'static + Support<N, RealField=Self::RealField>;
     /// An iterator over all the [`Support`]s of the generator.
     type AllDataIter<'a> : Iterator<Item=(Self::Id, Self::SupportType)> where Self : 'a;
 
diff -r edb95d2b83cc -r d2acaaddd9af src/euclidean.rs
--- a/src/euclidean.rs	Sun Nov 10 09:02:57 2024 -0500
+++ b/src/euclidean.rs	Tue Dec 31 09:12:43 2024 -0500
@@ -4,13 +4,14 @@
 
 use crate::types::*;
 use std::ops::{Mul, MulAssign, Div, DivAssign, Add, Sub, AddAssign, SubAssign, Neg};
+pub use crate::types::{HasScalarField, HasRealField};
 
 /// Space (type) with a defined dot product.
 ///
 /// `U` is the space of the multiplier, and `F` the space of scalars.
 /// Since `U` ≠ `Self`, this trait can also implement dual products.
-pub trait Dot<U, F> {
-    fn dot(&self, other : &U) -> F;
+pub trait Dot<U> : HasScalarField {
+    fn dot(&self, other : U) -> Self::Field;
 }
 
 /// Space (type) with Euclidean and vector space structure
@@ -18,59 +19,64 @@
 /// The type should implement vector space operations (addition, subtraction, scalar
 /// multiplication and scalar division) along with their assignment versions, as well
 /// as the [`Dot`] product with respect to `Self`.
-pub trait Euclidean<F : Float> : Sized + Dot<Self,F>
-        + Mul<F, Output=<Self as Euclidean<F>>::Output> + MulAssign<F>
-        + Div<F, Output=<Self as Euclidean<F>>::Output> + DivAssign<F>
-        + Add<Self, Output=<Self as Euclidean<F>>::Output>
-        + Sub<Self, Output=<Self as Euclidean<F>>::Output>
-        + for<'b> Add<&'b Self, Output=<Self as Euclidean<F>>::Output>
-        + for<'b> Sub<&'b Self, Output=<Self as Euclidean<F>>::Output>
+pub trait Euclidean : Sized
+        + HasRealField
+        + Dot<Self> + for<'a> Dot<&'a Self>
+        + Mul<Self::Field, Output=<Self as Euclidean>::Output>
+        + MulAssign<Self::Field>
+        + Div<Self::Field, Output=<Self as Euclidean>::Output>
+        + DivAssign<Self::Field>
+        + Add<Self, Output=<Self as Euclidean>::Output>
+        + Sub<Self, Output=<Self as Euclidean>::Output>
+        + for<'b> Add<&'b Self, Output=<Self as Euclidean>::Output>
+        + for<'b> Sub<&'b Self, Output=<Self as Euclidean>::Output>
         + AddAssign<Self> + for<'b> AddAssign<&'b Self>
         + SubAssign<Self> + for<'b> SubAssign<&'b Self>
-        + Neg<Output=<Self as Euclidean<F>>::Output> {
-    type Output : Euclidean<F>;
+        + Neg<Output=<Self as Euclidean>::Output> {
+
+    type Output : Euclidean<Field = Self::Field>;
 
     /// Returns origin of same dimensions as `self`.
-    fn similar_origin(&self) -> <Self as Euclidean<F>>::Output;
+    fn similar_origin(&self) -> <Self as Euclidean>::Output;
 
     /// Calculate the square of the 2-norm, $\frac{1}{2}\\|x\\|_2^2$, where `self` is $x$.
     #[inline]
-    fn norm2_squared(&self) -> F {
+    fn norm2_squared(&self) -> Self::Field {
         self.dot(self)
     }
 
     /// Calculate the square of the 2-norm divided by 2, $\frac{1}{2}\\|x\\|_2^2$,
     /// where `self` is $x$.
     #[inline]
-    fn norm2_squared_div2(&self) -> F {
-        self.norm2_squared()/F::TWO
+    fn norm2_squared_div2(&self) -> Self::Field {
+        self.norm2_squared()/Self::Field::TWO
     }
 
     /// Calculate the 2-norm $‖x‖_2$, where `self` is $x$.
     #[inline]
-    fn norm2(&self) -> F {
+    fn norm2(&self) -> Self::Field {
         self.norm2_squared().sqrt()
     }
 
     /// Calculate the 2-distance squared $\\|x-y\\|_2^2$, where `self` is $x$.
-    fn dist2_squared(&self, y : &Self) -> F;
+    fn dist2_squared(&self, y : &Self) -> Self::Field;
 
     /// Calculate the 2-distance $\\|x-y\\|_2$, where `self` is $x$.
     #[inline]
-    fn dist2(&self, y : &Self) -> F {
+    fn dist2(&self, y : &Self) -> Self::Field {
         self.dist2_squared(y).sqrt()
     }
 
     /// Projection to the 2-ball.
     #[inline]
-    fn proj_ball2(mut self, ρ : F) -> Self {
+    fn proj_ball2(mut self, ρ : Self::Field) -> Self {
         self.proj_ball2_mut(ρ);
         self
     }
 
     /// In-place projection to the 2-ball.
     #[inline]
-    fn proj_ball2_mut(&mut self, ρ : F) {
+    fn proj_ball2_mut(&mut self, ρ : Self::Field) {
         let r = self.norm2();
         if r>ρ {
             *self *= ρ/r
@@ -79,7 +85,7 @@
 }
 
 /// Trait for [`Euclidean`] spaces with dimensions known at compile time.
-pub trait StaticEuclidean<F : Float> : Euclidean<F> {
+pub trait StaticEuclidean : Euclidean {
     /// Returns the origin
-    fn origin() -> <Self as Euclidean<F>>::Output;
+    fn origin() -> <Self as Euclidean>::Output;
 }
diff -r edb95d2b83cc -r d2acaaddd9af src/linops.rs
--- a/src/linops.rs	Sun Nov 10 09:02:57 2024 -0500
+++ b/src/linops.rs	Tue Dec 31 09:12:43 2024 -0500
@@ -10,7 +10,8 @@
 
 /// Trait for linear operators on `X`.
 pub trait Linear<X> : Apply<X, Output=Self::Codomain>
-                      + for<'a> Apply<&'a X, Output=Self::Codomain> {
+                      + for<'a> Apply<&'a X, Output=Self::Codomain>
+                      + HasScalarField {
     type Codomain;
 }
 
@@ -18,32 +19,32 @@
 #[replace_float_literals(F::cast_from(literal))]
 pub trait AXPY<F : Num, X = Self> {
     /// Computes  `y = βy + αx`, where `y` is `Self`.
-    fn axpy(&mut self, α : F, x : &X, β : F);
+    fn axpy(&mut self, α : F, x : X, β : F);
 
     /// Copies `x` to `self`.
-    fn copy_from(&mut self, x : &X) {
+    fn copy_from(&mut self, x : X) {
         self.axpy(1.0, x, 0.0)
     }
 
     /// Computes  `y = αx`, where `y` is `Self`.
-    fn scale_from(&mut self, α : F, x : &X) {
+    fn scale_from(&mut self, α : F, x : X) {
         self.axpy(α, x, 0.0)
     }
 }
 
 /// Efficient in-place application for [`Linear`] operators.
 #[replace_float_literals(F::cast_from(literal))]
-pub trait GEMV<F : Num, X, Y = <Self as Linear<X>>::Codomain> : Linear<X> {
+pub trait GEMV<F : Num, X, Y = <Self as Apply<X>>::Output> : Apply<X> {
     /// Computes  `y = αAx + βy`, where `A` is `Self`.
-    fn gemv(&self, y : &mut Y, α : F, x : &X, β : F);
+    fn gemv(&self, y : &mut Y, α : F, x : X, β : F);
 
     /// Computes `y = Ax`, where `A` is `Self`
-    fn apply_mut(&self, y : &mut Y, x : &X){
+    fn apply_mut(&self, y : &mut Y, x : X){
         self.gemv(y, 1.0, x, 0.0)
     }
 
     /// Computes `y += Ax`, where `A` is `Self`
-    fn apply_add(&self, y : &mut Y, x : &X){
+    fn apply_add(&self, y : &mut Y, x : X){
         self.gemv(y, 1.0, x, 1.0)
     }
 }
@@ -51,11 +52,10 @@
 
 /// Bounded linear operators
 pub trait BoundedLinear<X> : Linear<X> {
-    type FloatType : Float;
     /// A bound on the operator norm $\|A\|$ for the linear operator $A$=`self`.
     /// This is not expected to be the norm, just any bound on it that can be
     /// reasonably implemented.
-    fn opnorm_bound(&self) -> Self::FloatType;
+    fn opnorm_bound(&self) -> Self::Field;
 }
 
 // Linear operator application into mutable target. The [`AsRef`] bound
@@ -70,7 +70,7 @@
 
 /// Trait for forming the adjoint operator of `Self`.
 pub trait Adjointable<X,Yʹ> : Linear<X> {
-    type AdjointCodomain;
+    type AdjointCodomain : HasScalarField<Field=Self::Field>;
     type Adjoint<'a> : Linear<Yʹ, Codomain=Self::AdjointCodomain> where Self : 'a;
 
     /// Form the adjoint operator of `self`.
@@ -89,7 +89,7 @@
 /// domain of the adjointed operator. `Self::Preadjoint` should be
 /// [`Adjointable`]`<'a,Ypre,X>`.
 pub trait Preadjointable<X,Ypre> : Linear<X> {
-    type PreadjointCodomain;
+    type PreadjointCodomain : HasScalarField<Field=Self::Field>;
     type Preadjoint<'a> : Adjointable<Ypre, X, Codomain=Self::PreadjointCodomain> where Self : 'a;
 
     /// Form the preadjoint operator of `self`.
@@ -100,55 +100,49 @@
 pub trait SimplyAdjointable<X> : Adjointable<X,<Self as Linear<X>>::Codomain> {}
 impl<'a,X,T> SimplyAdjointable<X> for T where T : Adjointable<X,<Self as Linear<X>>::Codomain> {}
 
-/// The identity operator
+/// The identity operator on `X` with scalar field `F`.
 #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)]
-pub struct IdOp<X> (PhantomData<X>);
+pub struct IdOp<X, F : Num> (PhantomData<(X, F)>);
 
-impl<X> IdOp<X> {
-    fn new() -> IdOp<X> { IdOp(PhantomData) }
+impl<X, F : Num> HasScalarField for IdOp<X, F> {
+    type Field = F;
 }
 
-impl<X> Apply<X> for IdOp<X> {
+impl<X, F : Num> IdOp<X, F> {
+    fn new() -> IdOp<X, F> { IdOp(PhantomData) }
+}
+
+impl<X, F : Num, T : CloneIfNeeded<X>> Apply<T> for IdOp<X, F> {
     type Output = X;
 
-    fn apply(&self, x : X) -> X {
-        x
+    fn apply(&self, x : T) -> X {
+        x.clone_if_needed()
     }
 }
 
-impl<'a, X> Apply<&'a X> for IdOp<X> where X : Clone {
-    type Output = X;
-    
-    fn apply(&self, x : &'a X) -> X {
-        x.clone()
-    }
-}
-
-impl<X> Linear<X> for IdOp<X> where X : Clone {
+impl<X : Clone, F : Num> Linear<X> for IdOp<X, F> {
     type Codomain = X;
 }
 
-
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Num, X, Y> GEMV<F, X, Y> for IdOp<X> where Y : AXPY<F, X>, X : Clone {
+impl<F : Num, X : Clone, Y> GEMV<F, X, Y> for IdOp<X, F> where Y : AXPY<F, X> {
     // Computes  `y = αAx + βy`, where `A` is `Self`.
-    fn gemv(&self, y : &mut Y, α : F, x : &X, β : F) {
+    fn gemv(&self, y : &mut Y, α : F, x : X, β : F) {
         y.axpy(α, x, β)
     }
 
-    fn apply_mut(&self, y : &mut Y, x : &X){
+    fn apply_mut(&self, y : &mut Y, x : X){
         y.copy_from(x);
     }
 }
 
-impl<X> BoundedLinear<X> for IdOp<X> where X : Clone {
-    type FloatType = float;
-    fn opnorm_bound(&self) -> float { 1.0 }
+impl<X, F : Num> BoundedLinear<X> for IdOp<X, F> where X : Clone {
+    fn opnorm_bound(&self) -> F { F::ONE }
 }
 
-impl<X> Adjointable<X,X> for IdOp<X> where X : Clone {
+impl<X, F : Num> Adjointable<X,X> for IdOp<X, F> where X : Clone + HasScalarField<Field=F> {
     type AdjointCodomain=X;
-    type Adjoint<'a> = IdOp<X> where X : 'a;
+    type Adjoint<'a> = IdOp<X, F> where X : 'a;
     fn adjoint(&self) -> Self::Adjoint<'_> { IdOp::new() }
 }
 
diff -r edb95d2b83cc -r d2acaaddd9af src/loc.rs
--- a/src/loc.rs	Sun Nov 10 09:02:57 2024 -0500
+++ b/src/loc.rs	Tue Dec 31 09:12:43 2024 -0500
@@ -6,6 +6,7 @@
 use std::ops::{Add,Sub,AddAssign,SubAssign,Mul,Div,MulAssign,DivAssign,Neg,Index,IndexMut};
 use std::slice::{Iter,IterMut};
 use std::fmt::{Display, Formatter};
+use std::borrow::Borrow;
 use crate::types::{Float,Num,SignedNum};
 use crate::maputil::{FixedLength,FixedLengthMut,map1,map2,map1_mut,map2_mut};
 use crate::euclidean::*;
@@ -23,6 +24,10 @@
     pub [F; N]
 );
 
+impl<F : Num, const N : usize> HasScalarField for Loc<F, N> {
+    type Field = F;
+}
+
 impl<F : Display, const N : usize> Display for Loc<F, N>{
     // Required method
     fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
@@ -393,7 +398,7 @@
 
 macro_rules! domination {
     ($norm:ident, $dominates:ident) => {
-        impl<F : Float, const N : usize> Dominated<F, $dominates, Loc<F, N>> for $norm {
+        impl<F : Float, const N : usize> Dominated<$dominates, Loc<F, N>> for $norm {
             #[inline]
             fn norm_factor(&self, _p : $dominates) -> F {
                 F::ONE
@@ -405,7 +410,7 @@
         }
     };
     ($norm:ident, $dominates:ident, $fn:path) => {
-        impl<F : Float, const N : usize> Dominated<F, $dominates, Loc<F, N>> for $norm {
+        impl<F : Float, const N : usize> Dominated<$dominates, Loc<F, N>> for $norm {
             #[inline]
             fn norm_factor(&self, _p : $dominates) -> F {
                 $fn(F::cast_from(N))
@@ -426,18 +431,18 @@
 domination!(Linfinity, L2);
 domination!(L2, L1);
 
-impl<F : Num,const N : usize> Dot<Loc<F, N>,F> for Loc<F, N> {
+impl<F : Num, T : Borrow<Loc<F, N>>, const N : usize> Dot<T> for Loc<F, N> {
     /// This implementation is not stabilised as it's meant to be used for very small vectors.
     /// Use [`nalgebra`] for larger vectors.
     #[inline]
-    fn dot(&self, other : &Loc<F, N>) -> F {
+    fn dot(&self, other : T) -> F {
         self.0.iter()
-              .zip(other.0.iter())
+              .zip(other.borrow().0.iter())
               .fold(F::ZERO, |m, (&v, &w)| m + v * w)
     }
 }
 
-impl<F : Float,const N : usize> Euclidean<F> for Loc<F, N> {
+impl<F : Float,const N : usize> Euclidean for Loc<F, N> {
     type Output = Self;
 
     #[inline]
@@ -483,24 +488,24 @@
     pub const ORIGIN : Self = Loc([F::ZERO; N]);
 }
 
-impl<F : Float,const N : usize> StaticEuclidean<F> for Loc<F, N> {
+impl<F : Float,const N : usize> StaticEuclidean for Loc<F, N> {
     #[inline]
     fn origin() -> Self {
         Self::ORIGIN
     }
 }
 
-impl<F : Float, const N : usize> Norm<F, L2> for Loc<F, N> {
+impl<F : Float, const N : usize> Norm<L2> for Loc<F, N> {
     #[inline]
     fn norm(&self, _ : L2) -> F { self.norm2() }
 }
 
-impl<F : Float, const N : usize> Dist<F, L2> for Loc<F, N> {
+impl<F : Float, const N : usize> Dist<L2> for Loc<F, N> {
     #[inline]
     fn dist(&self, other : &Self, _ : L2) -> F { self.dist2(other) }
 }
 
-impl<F : Float, const N : usize> Norm<F, L1> for Loc<F, N> {
+impl<F : Float, const N : usize> Norm<L1> for Loc<F, N> {
     /// This implementation is not stabilised as it's meant to be used for very small vectors.
     /// Use [`nalgebra`] for larger vectors.
     #[inline]
@@ -509,7 +514,7 @@
     }
 }
 
-impl<F : Float, const N : usize> Dist<F, L1> for Loc<F, N> {
+impl<F : Float, const N : usize> Dist<L1> for Loc<F, N> {
     #[inline]
     fn dist(&self, other : &Self, _ : L1) -> F {
         self.iter()
@@ -518,14 +523,14 @@
     }
 }
 
-impl<F : Float, const N : usize> Projection<F, Linfinity> for Loc<F, N> {
+impl<F : Float, const N : usize> Projection<Linfinity> for Loc<F, N> {
     #[inline]
     fn proj_ball_mut(&mut self, ρ : F, _ : Linfinity) {
         self.iter_mut().for_each(|v| *v = num_traits::clamp(*v, -ρ, ρ))
     }
 }
 
-impl<F : Float, const N : usize> Norm<F, Linfinity> for Loc<F, N> {
+impl<F : Float, const N : usize> Norm<Linfinity> for Loc<F, N> {
     /// This implementation is not stabilised as it's meant to be used for very small vectors.
     /// Use [`nalgebra`] for larger vectors.
     #[inline]
@@ -534,7 +539,7 @@
     }
 }
 
-impl<F : Float, const N : usize> Dist<F, Linfinity> for Loc<F, N> {
+impl<F : Float, const N : usize> Dist<Linfinity> for Loc<F, N> {
     #[inline]
     fn dist(&self, other : &Self, _ : Linfinity) -> F {
         self.iter()
@@ -572,19 +577,18 @@
     }
 }
 
-impl<F : Num, const N : usize> AXPY<F, Loc<F, N>> for Loc<F, N> {
-
+impl<'a, T : Borrow<Loc<F, N>>, F : Num, const N : usize> AXPY<F, T> for Loc<F, N> {
     #[inline]
-    fn axpy(&mut self, α : F, x : &Loc<F, N>, β : F) {
+    fn axpy(&mut self, α : F, x : T, β : F) {
         if β == F::ZERO {
-            map2_mut(self, x, |yi, xi| { *yi = α * (*xi) })
+            map2_mut(self, x.borrow(), |yi, xi| { *yi = α * (*xi) })
         } else {
-            map2_mut(self, x, |yi, xi| { *yi = β * (*yi) + α * (*xi) })
+            map2_mut(self, x.borrow(), |yi, xi| { *yi = β * (*yi) + α * (*xi) })
         }
     }
 
     #[inline]
-    fn copy_from(&mut self, x : &Loc<F, N>) {
-        map2_mut(self, x, |yi, xi| *yi = *xi )
+    fn copy_from(&mut self, x : T) {
+        map2_mut(self, x.borrow(), |yi, xi| *yi = *xi )
     }
 }
diff -r edb95d2b83cc -r d2acaaddd9af src/mapping.rs
--- a/src/mapping.rs	Sun Nov 10 09:02:57 2024 -0500
+++ b/src/mapping.rs	Tue Dec 31 09:12:43 2024 -0500
@@ -3,9 +3,10 @@
 */
 
 use std::marker::PhantomData;
-use crate::types::Float;
+use crate::types::{Float, HasRealField, HasScalarField};
 use serde::Serialize;
 use crate::loc::Loc;
+use std::borrow::Cow;
 
 /// Trait for application of `Self` as a mathematical function or operator on `X`.
 pub trait Apply<X> {
@@ -44,43 +45,47 @@
 
 /// Automatically implemented shorthand for referring to [`Mapping`]s from [`Loc<F, N>`] to `F`.
 pub trait RealMapping<F : Float, const N : usize>
-: Mapping<Loc<F, N>, Codomain = F> {}
+: Mapping<Loc<F, N>, Codomain = F> + HasRealField<RealField = F> {}
 
 impl<F : Float, T, const N : usize> RealMapping<F, N> for T
-where T : Mapping<Loc<F, N>, Codomain = F> {}
+where T : Mapping<Loc<F, N>, Codomain = F> + HasRealField<RealField = F> {}
 
 /// Automatically implemented shorthand for referring to differentiable [`Mapping`]s from
 /// [`Loc<F, N>`] to `F`.
 pub trait DifferentiableRealMapping<F : Float, const N : usize>
-: DifferentiableMapping<Loc<F, N>, Codomain = F, DerivativeDomain=Loc<F, N>> {}
+: DifferentiableMapping<Loc<F, N>, Codomain = F, DerivativeDomain=Loc<F, N>>
+  + RealMapping<F, N>
+  + HasRealField<RealField=F> {}
 
 impl<F : Float, T, const N : usize> DifferentiableRealMapping<F, N> for T
-where T : DifferentiableMapping<Loc<F, N>, Codomain = F, DerivativeDomain=Loc<F, N>> {}
+where T : DifferentiableMapping<Loc<F, N>, Codomain = F, DerivativeDomain=Loc<F, N>>
+          + RealMapping<F, N>
+          + HasRealField<RealField=F> {}
 
 
 /// A helper trait alias for referring to [`Mapping`]s from [`Loc<F, N>`] to [`Loc<F, M>`].
 pub trait RealVectorField<F : Float, const N : usize, const M : usize>
-: Mapping<Loc<F, N>, Codomain = Loc<F, M>> {}
+: Mapping<Loc<F, N>, Codomain = Loc<F, M>> + HasRealField<RealField=F> {}
 
 impl<F : Float, T, const N : usize, const M : usize> RealVectorField<F, N, M> for T
-where T : Mapping<Loc<F, N>, Codomain = Loc<F, M>> {}
+where T : Mapping<Loc<F, N>, Codomain = Loc<F, M>> + HasRealField<RealField=F> {}
 
 
 /// Trait for calculation the differential of `Self` as a mathematical function on `X`.
-pub trait Differentiable<X> : Sized {
-    type Derivative;
+pub trait Differentiable<X> : Sized + HasRealField {
+    type Derivative : HasRealField<RealField = Self::RealField>;
 
     /// Compute the differential of `self` at `x`.
     fn differential(&self, x : X) -> Self::Derivative;
 }
 
-impl<'g, X, G : Differentiable<X>> Differentiable<X> for &'g G {
-    type Derivative = G::Derivative;
-    #[inline]
-    fn differential(&self, x : X) -> Self::Derivative {
-        (*self).differential(x)
-    }
-}
+// impl<'g, X, G : Differentiable<X>> Differentiable<X> for &'g G {
+//     type Derivative = G::Derivative;
+//     #[inline]
+//     fn differential(&self, x : X) -> Self::Derivative {
+//         (*self).differential(x)
+//     }
+// }
 
 /// A differentiable mapping from `Domain` to [`Mapping::Codomain`], with differentials
 /// `Differential`.
@@ -90,33 +95,34 @@
 : Mapping<Domain>
   + Differentiable<Domain, Derivative=Self::DerivativeDomain>
   + for<'a> Differentiable<&'a Domain, Derivative=Self::DerivativeDomain> {
-    type DerivativeDomain;
-    type Differential : Mapping<Domain, Codomain=Self::DerivativeDomain>;
-    type DifferentialRef<'b> : Mapping<Domain, Codomain=Self::DerivativeDomain> where Self : 'b;
+    type DerivativeDomain : HasRealField<RealField=Self::RealField>;
+    type Differential<'b> : Mapping<Domain, Codomain=Self::DerivativeDomain>
+                            + HasRealField<RealField=Self::RealField>
+                            + Clone where Self : 'b;
 
     /// Form the differential mapping of `self`.
-    fn diff(self) -> Self::Differential;
+    fn diff(self) -> Self::Differential<'static>;
     /// Form the differential mapping of `self`.
-    fn diff_ref(&self) -> Self::DifferentialRef<'_>;
+    fn diff_ref(&self) -> Self::Differential<'_>;
 }
 
 
-impl<Domain, Derivative, T> DifferentiableMapping<Domain> for T
-where T : Mapping<Domain>
+impl<Domain : Clone, Derivative, T> DifferentiableMapping<Domain> for T
+where T : Mapping<Domain> + Clone
           + Differentiable<Domain, Derivative=Derivative>
-          + for<'a> Differentiable<&'a Domain,Derivative=Derivative> {
+          + for<'a> Differentiable<&'a Domain, Derivative=Derivative>,
+      Derivative : HasRealField<RealField = T::RealField> {
     type DerivativeDomain = Derivative;
-    type Differential = Differential<Domain, Self>;
-    type DifferentialRef<'b> = Differential<Domain, &'b Self> where Self : 'b;
+    type Differential<'b> = Differential<'b, Domain, Self> where Self : 'b + Clone;
     
     /// Form the differential mapping of `self`.
-    fn diff(self) -> Self::Differential {
-        Differential{ g : self, _space : PhantomData }
+    fn diff(self) -> Self::Differential<'static> {
+        Differential{ g : Cow::Owned(self), _space : PhantomData }
     }
 
     /// Form the differential mapping of `self`.
-    fn diff_ref(&self) -> Self::DifferentialRef<'_> {
-        Differential{ g : self, _space : PhantomData }
+    fn diff_ref(&self) -> Self::Differential<'_> {
+        Differential{ g : Cow::Borrowed(self), _space : PhantomData }
     }
 }
 
@@ -127,6 +133,10 @@
     _domain : PhantomData<Domain>,
 }
 
+impl<Domain, M : Mapping<Domain> + HasScalarField> HasScalarField for Sum<Domain, M> {
+    type Field = M::Field;
+}
+
 impl<Domain, M : Mapping<Domain>> Sum<Domain, M> {
     /// Construct from an iterator.
     pub fn new<I : Iterator<Item = M>>(iter : I) -> Self {
@@ -174,18 +184,25 @@
 }
 
 /// Container for the differential [`Mapping`] of a [`Differentiable`] mapping.
-pub struct Differential<X, G : DifferentiableMapping<X>> {
-    g : G,
+#[derive(Clone, Debug)]
+pub struct Differential<'a, X, G : 'a + DifferentiableMapping<X> + Clone> {
+    g : Cow<'a, G>,
     _space : PhantomData<X>
 }
 
-impl<X, G : DifferentiableMapping<X>> Differential<X, G> {
-    pub fn base_fn(&self) -> &G {
-        &self.g
+impl<'a, X, G : DifferentiableMapping<X> + Clone> Differential<'a, X, G> {
+    pub fn base_fn(&'a self) -> &'a G {
+        &*self.g
     }
 }
 
-impl<X, G : DifferentiableMapping<X>> Apply<X> for Differential<X, G> {
+
+impl<'a, X, G> HasScalarField for Differential<'a, X, G>
+where G : 'a + DifferentiableMapping<X> + HasScalarField + Clone {
+    type Field = G::Field;
+}
+
+impl<'b, X, G : DifferentiableMapping<X> + Clone> Apply<X> for Differential<'b, X, G> {
     type Output = G::DerivativeDomain;
 
     #[inline]
@@ -194,7 +211,8 @@
     }
 }
 
-impl<'a, X, G : DifferentiableMapping<X>> Apply<&'a X> for Differential<X, G> {
+impl<'a, 'b, X, G : 'b + DifferentiableMapping<X> + Clone> Apply<&'a X>
+for Differential<'b, X, G> {
     type Output = G::DerivativeDomain;
 
     #[inline]
@@ -232,14 +250,23 @@
 
 
 /// Container for dimensional slicing [`Loc`]`<F, N>` codomain of a [`Mapping`] to `F`.
-pub struct SlicedCodomain<X, F, G : Mapping<X, Codomain=Loc<F, N>>, const N : usize> {
-    g : G,
+pub struct SlicedCodomain<
+    'a, X, F,
+    G : Mapping<X, Codomain=Loc<F, N>> + Clone,
+    const N : usize
+> {
+    g : Cow<'a, G>,
     slice : usize,
     _phantoms : PhantomData<(X, F)>
 }
 
-impl<X, F : Copy, G : Mapping<X, Codomain=Loc<F, N>>, const N : usize> Apply<X>
-for SlicedCodomain<X, F, G, N> {
+impl<'a, X, F, G, const N : usize> HasScalarField for SlicedCodomain<'a, X, F, G, N>
+where G : HasScalarField + Mapping<X, Codomain=Loc<F, N>> + Clone {
+    type Field = G::Field;
+}
+
+impl<'a, X, F : Copy, G : Mapping<X, Codomain=Loc<F, N>> + Clone, const N : usize> Apply<X>
+for SlicedCodomain<'a, X, F, G, N> {
     type Output = F;
 
     #[inline]
@@ -250,12 +277,13 @@
     }
 }
 
-impl<'a, X, F : Copy, G : Mapping<X, Codomain=Loc<F, N>>, const N : usize> Apply<&'a X>
-for SlicedCodomain<X, F, G, N> {
+impl<'a, 'b, X, F : Copy, G : Mapping<X, Codomain=Loc<F, N>> + Clone, const N : usize>
+Apply<&'b X>
+for SlicedCodomain<'a, X, F, G, N> {
     type Output = F;
 
     #[inline]
-    fn apply(&self, x : &'a X) -> Self::Output {
+    fn apply(&self, x : &'b X) -> Self::Output {
         let tmp : [F; N] = self.g.apply(x).into();
         // Safety: `slice_codomain` below checks the range.
         unsafe { *tmp.get_unchecked(self.slice) }
@@ -264,21 +292,21 @@
 
 /// An auto-trait for constructing a [`FlattenCodomain`] structure for
 /// flattening the codomain of a [`Mapping`] from [`Loc`]`<F, 1>` to `F`.
-pub trait SliceCodomain<X, F : Copy, const N : usize> : Mapping<X, Codomain=Loc<F, N>> + Sized {
+pub trait SliceCodomain<X, F : Copy, const N : usize> : Mapping<X, Codomain=Loc<F, N>> + Clone + Sized {
     /// Flatten the codomain from [`Loc`]`<F, 1>` to `F`.
-    fn slice_codomain(self, slice : usize) -> SlicedCodomain<X, F, Self, N> {
+    fn slice_codomain(self, slice : usize) -> SlicedCodomain<'static, X, F, Self, N> {
         assert!(slice < N);
-        SlicedCodomain{ g : self, slice, _phantoms : PhantomData }
+        SlicedCodomain{ g : Cow::Owned(self), slice, _phantoms : PhantomData }
     }
 
     /// Flatten the codomain from [`Loc`]`<F, 1>` to `F`.
-    fn slice_codomain_ref(&self, slice : usize) -> SlicedCodomain<X, F, &'_ Self, N> {
+    fn slice_codomain_ref(&self, slice : usize) -> SlicedCodomain<'_, X, F, Self, N> {
         assert!(slice < N);
-        SlicedCodomain{ g : self, slice, _phantoms : PhantomData }
+        SlicedCodomain{ g : Cow::Borrowed(self), slice, _phantoms : PhantomData }
     }
 }
 
-impl<X, F : Copy, G : Sized + Mapping<X, Codomain=Loc<F, N>>, const N : usize>
+impl<X, F : Copy, G : Sized + Mapping<X, Codomain=Loc<F, N>> + Clone, const N : usize>
 SliceCodomain<X, F, N>
 for G {}
 
diff -r edb95d2b83cc -r d2acaaddd9af src/nalgebra_support.rs
--- a/src/nalgebra_support.rs	Sun Nov 10 09:02:57 2024 -0500
+++ b/src/nalgebra_support.rs	Tue Dec 31 09:12:43 2024 -0500
@@ -56,9 +56,9 @@
     }
 }
 
-impl<'a, SM,SV,N,M,K,E> Linear<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM>
+impl<SM,SV,N,M,K,E> Linear<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM>
 where SM: Storage<E,N,M>, SV: Storage<E,M,K>,
-        N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One,
+        N : Dim, M : Dim, K : Dim, E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One,
         DefaultAllocator : Allocator<E,N,K>,
         DefaultAllocator : Allocator<E,M,K>,
         DefaultAllocator : Allocator<E,N,M>,
@@ -75,12 +75,31 @@
       DefaultAllocator : Allocator<E,M,N> {
 
     #[inline]
-    fn gemv(&self, y : &mut Matrix<E,N,K,SV2>, α : E, x : &Matrix<E,M,K,SV1>, β : E) {
+    fn gemv(&self, y : &mut Matrix<E,N,K,SV2>, α : E, x : Matrix<E,M,K,SV1>, β : E) {
+        Matrix::gemm(y, α, self, &x, β)
+    }
+
+    #[inline]
+    fn apply_mut(&self, y : &mut Matrix<E,N,K,SV2>, x : Matrix<E,M,K,SV1>) {
+        self.mul_to(&x, y)
+    }
+}
+
+impl<'a, SM,SV1,SV2,N,M,K,E> GEMV<E, &'a Matrix<E,M,K,SV1>, Matrix<E,N,K,SV2>> for Matrix<E,N,M,SM>
+where SM: Storage<E,N,M>, SV1: Storage<E,M,K>, SV2: StorageMut<E,N,K>,
+      N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + Float,
+      DefaultAllocator : Allocator<E,N,K>,
+      DefaultAllocator : Allocator<E,M,K>,
+      DefaultAllocator : Allocator<E,N,M>,
+      DefaultAllocator : Allocator<E,M,N> {
+
+    #[inline]
+    fn gemv(&self, y : &mut Matrix<E,N,K,SV2>, α : E, x : &'a Matrix<E,M,K,SV1>, β : E) {
         Matrix::gemm(y, α, self, x, β)
     }
 
     #[inline]
-    fn apply_mut<'a>(&self, y : &mut Matrix<E,N,K,SV2>, x : &Matrix<E,M,K,SV1>) {
+    fn apply_mut(&self, y : &mut Matrix<E,N,K,SV2>, x : &'a Matrix<E,M,K,SV1>) {
         self.mul_to(x, y)
     }
 }
@@ -91,17 +110,33 @@
       DefaultAllocator : Allocator<E,M> {
 
     #[inline]
-    fn axpy(&mut self, α : E, x : &Vector<E,M,SV1>, β : E) {
+    fn axpy(&mut self, α : E, x : Vector<E,M,SV1>, β : E) {
+        Matrix::axpy(self, α, &x, β)
+    }
+
+    #[inline]
+    fn copy_from(&mut self, y : Vector<E,M,SV1>) {
+        Matrix::copy_from(self, &y)
+    }
+}
+
+impl<'a, SM,SV1,M,E> AXPY<E, &'a Vector<E,M,SV1>> for Vector<E,M,SM>
+where SM: StorageMut<E,M>, SV1: Storage<E,M>,
+      M : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + Float,
+      DefaultAllocator : Allocator<E,M> {
+
+    #[inline]
+    fn axpy(&mut self, α : E, x : &'a Vector<E,M,SV1>, β : E) {
         Matrix::axpy(self, α, x, β)
     }
 
     #[inline]
-    fn copy_from(&mut self, y : &Vector<E,M,SV1>) {
+    fn copy_from(&mut self, y : &'a Vector<E,M,SV1>) {
         Matrix::copy_from(self, y)
     }
 }
 
-impl<SM,M,E> Projection<E, Linfinity> for Vector<E,M,SM>
+impl<SM,M,E> Projection<Linfinity> for Vector<E,M,SM>
 where SM: StorageMut<E,M>,
       M : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + Float + RealField,
       DefaultAllocator : Allocator<E,M> {
@@ -114,7 +149,8 @@
 impl<'own,SV1,SV2,SM,N,M,K,E> Adjointable<Matrix<E,M,K,SV1>,Matrix<E,N,K,SV2>>
 for Matrix<E,N,M,SM>
 where SM: Storage<E,N,M>, SV1: Storage<E,M,K>, SV2: Storage<E,N,K>,
-      N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + SimdComplexField,
+      N : Dim, M : Dim, K : Dim,
+      E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + SimdComplexField,
       DefaultAllocator : Allocator<E,N,K>,
       DefaultAllocator : Allocator<E,M,K>,
       DefaultAllocator : Allocator<E,N,M>,
@@ -128,7 +164,7 @@
     }
 }
 
-impl<E,M,S,Si> Dot<Vector<E,M,Si>,E>
+impl<E,M,S,Si> Dot<Vector<E,M,Si>>
 for Vector<E,M,S>
 where M : Dim,
       E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One,
@@ -137,7 +173,21 @@
       DefaultAllocator : Allocator<E,M> {
 
     #[inline]
-    fn dot(&self, other : &Vector<E,M,Si>) -> E {
+    fn dot(&self, other : Vector<E,M,Si>) -> E {
+        Vector::<E,M,S>::dot(self, &other)
+    }
+}
+
+impl<'a, E,M,S,Si> Dot<&'a Vector<E,M,Si>>
+for Vector<E,M,S>
+where M : Dim,
+      E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One,
+      S : Storage<E,M>,
+      Si : Storage<E,M>,
+      DefaultAllocator : Allocator<E,M> {
+
+    #[inline]
+    fn dot(&self, other : &'a Vector<E,M,Si>) -> E {
         Vector::<E,M,S>::dot(self, other)
     }
 }
@@ -167,7 +217,16 @@
 
 // TODO: should allow different input storages in `Euclidean`.
 
-impl<E,M,S> Euclidean<E>
+impl<E,M,K,S> HasScalarField
+for Matrix<E,M,K,S>
+where M : Dim, K : Dim,
+      S : Storage<E,M,K>,
+      E : Float + Scalar,
+      DefaultAllocator : Allocator<E,M,K> {
+    type Field = E;
+}
+
+impl<E,M,S> Euclidean
 for Vector<E,M,S>
 where M : Dim,
       S : StorageMut<E,M>,
@@ -192,7 +251,7 @@
     }
 }
 
-impl<E,M,S> StaticEuclidean<E>
+impl<E,M,S> StaticEuclidean
 for Vector<E,M,S>
 where M : DimName,
       S : StorageMut<E,M>,
@@ -205,7 +264,7 @@
     }
 }
 
-impl<E,M,S> Norm<E, L1>
+impl<E,M,S> Norm<L1>
 for Vector<E,M,S>
 where M : Dim,
       S : StorageMut<E,M>,
@@ -218,7 +277,7 @@
     }
 }
 
-impl<E,M,S> Dist<E, L1>
+impl<E,M,S> Dist<L1>
 for Vector<E,M,S>
 where M : Dim,
       S : StorageMut<E,M>,
@@ -230,7 +289,7 @@
     }
 }
 
-impl<E,M,S> Norm<E, L2>
+impl<E,M,S> Norm<L2>
 for Vector<E,M,S>
 where M : Dim,
       S : StorageMut<E,M>,
@@ -243,7 +302,7 @@
     }
 }
 
-impl<E,M,S> Dist<E, L2>
+impl<E,M,S> Dist<L2>
 for Vector<E,M,S>
 where M : Dim,
       S : StorageMut<E,M>,
@@ -255,7 +314,7 @@
     }
 }
 
-impl<E,M,S> Norm<E, Linfinity>
+impl<E,M,S> Norm<Linfinity>
 for Vector<E,M,S>
 where M : Dim,
       S : StorageMut<E,M>,
@@ -268,7 +327,7 @@
     }
 }
 
-impl<E,M,S> Dist<E, Linfinity>
+impl<E,M,S> Dist<Linfinity>
 for Vector<E,M,S>
 where M : Dim,
       S : StorageMut<E,M>,
diff -r edb95d2b83cc -r d2acaaddd9af src/norms.rs
--- a/src/norms.rs	Sun Nov 10 09:02:57 2024 -0500
+++ b/src/norms.rs	Tue Dec 31 09:12:43 2024 -0500
@@ -5,6 +5,7 @@
 use serde::Serialize;
 use crate::types::*;
 use crate::euclidean::*;
+pub use crate::types::{HasScalarField, HasRealField};
 
 //
 // Abstract norms
@@ -53,7 +54,6 @@
 pub struct HuberL21<F : Float>(pub F);
 impl<F : Float> NormExponent for HuberL21<F> {}
 
-
 /// A normed space (type) with exponent or other type `Exponent` for the norm.
 ///
 /// Use as
@@ -64,27 +64,28 @@
 ///
 /// println!("{}, {} {}", x.norm(L1), x.norm(L2), x.norm(Linfinity))
 /// ```
-pub trait Norm<F : Num, Exponent : NormExponent> {
+pub trait Norm<Exponent : NormExponent> : HasScalarField {
     /// Calculate the norm.
-    fn norm(&self, _p : Exponent) -> F;
+    fn norm(&self, _p : Exponent) -> Self::Field;
 }
 
 /// Indicates that the `Self`-[`Norm`] is dominated by the `Exponent`-`Norm` on the space
 /// `Elem` with the corresponding field `F`.
-pub trait Dominated<F : Num, Exponent : NormExponent, Elem> {
+pub trait Dominated<Exponent : NormExponent, Elem>
+where Elem : HasScalarField {
     /// Indicates the factor $c$ for the inequality $‖x‖ ≤ C ‖x‖_p$.
-    fn norm_factor(&self, p : Exponent) -> F;
+    fn norm_factor(&self, p : Exponent) -> Elem::Field;
     /// Given a norm-value $‖x‖_p$, calculates $C‖x‖_p$ such that $‖x‖ ≤ C‖x‖_p$
     #[inline]
-    fn from_norm(&self, p_norm : F, p : Exponent) -> F {
+    fn from_norm(&self, p_norm : Elem::Field, p : Exponent) -> Elem::Field {
         p_norm * self.norm_factor(p)
     }
 }
 
 /// Trait for distances with respect to a norm.
-pub trait Dist<F : Num, Exponent : NormExponent> : Norm<F, Exponent> {
+pub trait Dist<Exponent : NormExponent> : Norm<Exponent> {
     /// Calculate the distance
-    fn dist(&self, other : &Self, _p : Exponent) -> F;
+    fn dist(&self, other : &Self, _p : Exponent) -> Self::Field;
 }
 
 /// Trait for Euclidean projections to the `Exponent`-[`Norm`]-ball.
@@ -97,16 +98,15 @@
 ///
 /// println!("{:?}, {:?}", x.proj_ball(1.0, L2), x.proj_ball(0.5, Linfinity));
 /// ```
-pub trait Projection<F : Num, Exponent : NormExponent> : Norm<F, Exponent> + Euclidean<F>
-where F : Float {
+pub trait Projection<Exponent : NormExponent> : Norm<Exponent> + Euclidean {
     /// Projection of `self` to the `q`-norm-ball of radius ρ.
-    fn proj_ball(mut self, ρ : F, q : Exponent) -> Self {
+    fn proj_ball(mut self, ρ : Self::Field, q : Exponent) -> Self {
         self.proj_ball_mut(ρ, q);
         self
     }
 
     /// In-place projection of `self` to the `q`-norm-ball of radius ρ.
-    fn proj_ball_mut(&mut self, ρ : F, _q : Exponent);
+    fn proj_ball_mut(&mut self, ρ : Self::Field, _q : Exponent);
 }
 
 /*impl<F : Float, E : Euclidean<F>> Norm<F, L2> for E {
@@ -116,12 +116,12 @@
     fn dist(&self, other : &Self, _p : L2) -> F { self.dist2(other) }
 }*/
 
-impl<F : Float, E : Euclidean<F> + Norm<F, L2>> Projection<F, L2> for E {
+impl<E : Norm<L2> + Euclidean> Projection<L2> for E {
     #[inline]
-    fn proj_ball(self, ρ : F, _p : L2) -> Self { self.proj_ball2(ρ) }
+    fn proj_ball(self, ρ : Self::Field, _p : L2) -> Self { self.proj_ball2(ρ) }
 
     #[inline]
-    fn proj_ball_mut(&mut self, ρ : F, _p : L2) { self.proj_ball2_mut(ρ) }
+    fn proj_ball_mut(&mut self, ρ : Self::Field, _p : L2) { self.proj_ball2_mut(ρ) }
 }
 
 impl<F : Float> HuberL1<F> {
@@ -142,13 +142,13 @@
     }
 }
 
-impl<F : Float, E : Euclidean<F>> Norm<F, HuberL1<F>> for E {
+impl<F : Float, E : Euclidean<RealField=F>> Norm<HuberL1<F>> for E {
     fn norm(&self, huber : HuberL1<F>) -> F {
         huber.apply(self.norm2_squared())
     }
 }
 
-impl<F : Float, E : Euclidean<F>> Dist<F, HuberL1<F>> for E {
+impl<F : Float, E : Euclidean<RealField=F>> Dist<HuberL1<F>> for E {
     fn dist(&self, other : &Self, huber : HuberL1<F>) -> F {
         huber.apply(self.dist2_squared(other))
     }
diff -r edb95d2b83cc -r d2acaaddd9af src/sets.rs
--- a/src/sets.rs	Sun Nov 10 09:02:57 2024 -0500
+++ b/src/sets.rs	Tue Dec 31 09:12:43 2024 -0500
@@ -59,15 +59,19 @@
 /// vector and $t$ the offset.
 ///
 /// `U` is the element type, `F` the floating point number type, and `A` the type of the
-/// orthogonal (dual) vectors. They need implement [`Dot<U, F>`].
+/// orthogonal (dual) vectors. They need implement [`Dot<U, Output=F>`].
 #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)]
-pub struct Halfspace<A, F, U> where A : Dot<U, F>, F : Float {
+pub struct Halfspace<A, F, U>
+where A : Dot<U, Field=F> + for<'a> Dot<&'a U>,
+      F : Float {
     pub orthogonal : A,
     pub offset : F,
     _phantom : PhantomData<U>,
 }
 
-impl<A,F,U> Halfspace<A,F,U> where A : Dot<U,F>, F : Float {
+impl<A,F,U> Halfspace<A,F,U>
+where A : Dot<U, Field=F> + for<'a> Dot<&'a U>,
+      F : Float {
     #[inline]
     pub fn new(orthogonal : A, offset : F) -> Self {
         Halfspace{ orthogonal : orthogonal, offset : offset, _phantom : PhantomData }
@@ -77,7 +81,7 @@
 /// Trait for generating a halfspace spanned by another set `Self` of elements of type `U`.
 pub trait SpannedHalfspace<F, U> where  F : Float {
     /// Type of the orthogonal vector describing the halfspace.
-    type A : Dot<U, F>;
+    type A : Dot<U, Field=F> + for<'a> Dot<&'a U>;
     /// Returns the halfspace spanned by this set.
     fn spanned_halfspace(&self) -> Halfspace<Self::A, F, U>;
 }
@@ -103,7 +107,9 @@
     }
 }
 
-impl<A,U,F> Set<U> for Halfspace<A,F,U> where A : Dot<U,F>, F : Float {
+impl<A,U,F> Set<U> for Halfspace<A,F,U>
+where A : Dot<U, Field=F> + for<'a> Dot<&'a U>,
+      F : Float {
     #[inline]
     fn contains(&self, item : &U) -> bool {
         self.orthogonal.dot(item) >= self.offset
@@ -111,10 +117,14 @@
 }
 
 /// Polygons defined by `N` `Halfspace`s.
-#[derive(Clone,Copy,Debug,Eq,PartialEq)]
-pub struct NPolygon<A, F, U, const N : usize>(pub [Halfspace<A,F,U>; N]) where A : Dot<U,F>, F : Float;
+#[derive(Clone,Debug,Eq,PartialEq)]
+pub struct NPolygon<A, F, U, const N : usize>(pub [Halfspace<A,F,U>; N])
+where A : Dot<U, Field=F> + for<'a> Dot<&'a U>,
+      F : Float;
 
-impl<A,U,F,const N : usize> Set<U> for NPolygon<A,F,U,N> where A : Dot<U,F>, F : Float {
+impl<A,U,F,const N : usize> Set<U> for NPolygon<A,F,U,N>
+where A : Dot<U, Field=F> + for<'a> Dot<&'a U>,
+      F : Float {
     fn contains(&self, item : &U) -> bool {
         self.0.iter().all(|halfspace| halfspace.contains(item))
     }
diff -r edb95d2b83cc -r d2acaaddd9af src/types.rs
--- a/src/types.rs	Sun Nov 10 09:02:57 2024 -0500
+++ b/src/types.rs	Tue Dec 31 09:12:43 2024 -0500
@@ -57,7 +57,8 @@
                  + CastFrom<u128> + CastFrom<usize>
                  + CastFrom<i8>   + CastFrom<i16> + CastFrom<i32> + CastFrom<i64>
                  + CastFrom<i128> + CastFrom<isize>
-                 + CastFrom<f32>  + CastFrom<f64> {
+                 + CastFrom<f32>  + CastFrom<f64>
+                 + HasScalarField<Field = Self> {
 
     const ZERO : Self;
     const ONE : Self;
@@ -106,6 +107,9 @@
             const RANGE_MAX : Self = <$type>::MAX;
             const RANGE_MIN : Self = <$type>::MIN;
         }
+        impl HasScalarField for $type {
+            type Field = $type;
+        }
     )* }
 }
 
@@ -163,3 +167,54 @@
 }
 */
 
+
+/// Trait for on-demand cloning as of `Self` as `Target`.
+///
+/// Blanket implementations clone references, but do nothing if `Self` equals `Target`.
+///
+/// The idea is to allow other traits to be implemented generically, similar to now one
+/// can work with references and owned values through [`std::borrow::Borrow`], but, in
+/// contrast to the latter, this trait will provide an owned value instead of reference.
+pub trait CloneIfNeeded<Target> {
+    /// Clone `self` if needed as `Target`.
+    fn clone_if_needed(self) -> Target;
+}
+
+impl<'a, T : Clone> CloneIfNeeded<T> for &'a T {
+    #[inline]
+    /// Clone
+    fn clone_if_needed(self) -> T {
+        self.clone()
+    }
+}
+
+impl<'a, T : Clone> CloneIfNeeded<T> for &'a mut T {
+    #[inline]
+    /// Clone
+    fn clone_if_needed(self) -> T {
+        self.clone()
+    }
+}
+
+impl<'a, T> CloneIfNeeded<T> for T {
+    #[inline]
+    /// No clone needed
+    fn clone_if_needed(self) -> T {
+        self
+    }
+}
+
+/// Trait for objects that have an associated scalar field.
+pub trait HasScalarField {
+    type Field : Num;
+}
+
+/// Trait for objects that have an associated real field.
+pub trait HasRealField : HasScalarField<Field = Self::RealField> {
+    type RealField : Float;
+}
+
+impl<T> HasRealField for T where T : HasScalarField, T::Field : Float {
+    type RealField = T::Field;
+}
+