diff -r 1f19c6bbf07b -r 3868555d135c src/direct_product.rs
--- a/src/direct_product.rs	Sun Apr 27 20:29:43 2025 -0500
+++ b/src/direct_product.rs	Fri May 15 14:46:30 2026 -0500
@@ -6,11 +6,11 @@
 */
 
 use crate::euclidean::Euclidean;
-use crate::instance::{Decomposition, DecompositionMut, Instance, InstanceMut, MyCow};
-use crate::linops::AXPY;
+use crate::instance::{Decomposition, DecompositionMut, Instance, InstanceMut, MyCow, Ownable};
+use crate::linops::{VectorSpace, AXPY};
 use crate::loc::Loc;
 use crate::mapping::Space;
-use crate::norms::{HasDual, Norm, NormExponent, Normed, PairNorm, L2};
+use crate::norms::{Dist, HasDual, Norm, NormExponent, Normed, PairNorm, L2};
 use crate::types::{Float, Num};
 use core::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign};
 use serde::{Deserialize, Serialize};
@@ -39,239 +39,290 @@
     }
 }
 
-macro_rules! impl_binop {
-    (($a : ty, $b : ty), $trait : ident, $fn : ident, $refl:ident, $refr:ident) => {
-        impl_binop!(@doit: $a, $b, $trait, $fn;
-                           maybe_lifetime!($refl, &'l Pair<$a,$b>),
-                           (maybe_lifetime!($refl, &'l $a),
-                            maybe_lifetime!($refl, &'l $b));
-                           maybe_lifetime!($refr, &'r Pair<Ai,Bi>),
-                           (maybe_lifetime!($refr, &'r Ai),
-                            maybe_lifetime!($refr, &'r Bi));
-                           $refl, $refr);
-    };
-
-    (@doit: $a:ty, $b:ty,
-            $trait:ident, $fn:ident;
-            $self:ty, ($aself:ty, $bself:ty);
-            $in:ty, ($ain:ty, $bin:ty);
-            $refl:ident, $refr:ident) => {
-        impl<'l, 'r, Ai, Bi> $trait<$in>
-        for $self
-        where $aself: $trait<$ain>,
-              $bself: $trait<$bin> {
-            type Output = Pair<<$aself as $trait<$ain>>::Output,
-                               <$bself as $trait<$bin>>::Output>;
-
-            #[inline]
-            fn $fn(self, y : $in) -> Self::Output {
-                Pair(maybe_ref!($refl, self.0).$fn(maybe_ref!($refr, y.0)),
-                     maybe_ref!($refl, self.1).$fn(maybe_ref!($refr, y.1)))
+macro_rules! impl_unary {
+    ($trait:ident, $fn:ident) => {
+        impl<A, B> $trait for Pair<A, B>
+        where
+            A: $trait,
+            B: $trait,
+        {
+            type Output = Pair<A::Output, B::Output>;
+            fn $fn(self) -> Self::Output {
+                let Pair(a, b) = self;
+                Pair(a.$fn(), b.$fn())
             }
         }
+
+        // Compiler overflow
+        // impl<'a, A, B> $trait for &'a Pair<A, B>
+        // where
+        //     &'a A: $trait,
+        //     &'a B: $trait,
+        // {
+        //     type Output = Pair<<&'a A as $trait>::Output, <&'a B as $trait>::Output>;
+        //     fn $fn(self) -> Self::Output {
+        //         let Pair(ref a, ref b) = self;
+        //         Pair(a.$fn(), b.$fn())
+        //     }
+        // }
     };
 }
 
-macro_rules! impl_assignop {
-    (($a : ty, $b : ty), $trait : ident, $fn : ident, $refr:ident) => {
-        impl_assignop!(@doit: $a, $b,
-                              $trait, $fn;
-                              maybe_lifetime!($refr, &'r Pair<Ai,Bi>),
-                              (maybe_lifetime!($refr, &'r Ai),
-                               maybe_lifetime!($refr, &'r Bi));
-                              $refr);
-    };
-    (@doit: $a : ty, $b : ty,
-            $trait:ident, $fn:ident;
-            $in:ty, ($ain:ty, $bin:ty);
-            $refr:ident) => {
-        impl<'r, Ai, Bi> $trait<$in>
-        for Pair<$a,$b>
-        where $a: $trait<$ain>,
-              $b: $trait<$bin> {
-            #[inline]
-            fn $fn(&mut self, y : $in) -> () {
-                self.0.$fn(maybe_ref!($refr, y.0));
-                self.1.$fn(maybe_ref!($refr, y.1));
+impl_unary!(Neg, neg);
+
+macro_rules! impl_binary {
+    ($trait:ident, $fn:ident) => {
+        impl<A, B, C, D> $trait<Pair<C, D>> for Pair<A, B>
+        where
+            A: $trait<C>,
+            B: $trait<D>,
+        {
+            type Output = Pair<A::Output, B::Output>;
+            fn $fn(self, Pair(c, d): Pair<C, D>) -> Self::Output {
+                let Pair(a, b) = self;
+                Pair(a.$fn(c), b.$fn(d))
             }
         }
-    }
-}
 
-macro_rules! impl_scalarop {
-    (($a : ty, $b : ty), $field : ty, $trait : ident, $fn : ident, $refl:ident) => {
-        impl_scalarop!(@doit: $field,
-                              $trait, $fn;
-                              maybe_lifetime!($refl, &'l Pair<$a,$b>),
-                              (maybe_lifetime!($refl, &'l $a),
-                               maybe_lifetime!($refl, &'l $b));
-                              $refl);
-    };
-    (@doit: $field : ty,
-            $trait:ident, $fn:ident;
-            $self:ty, ($aself:ty, $bself:ty);
-            $refl:ident) => {
-        // Scalar as Rhs
-        impl<'l> $trait<$field>
-        for $self
-        where $aself: $trait<$field>,
-              $bself: $trait<$field> {
-            type Output = Pair<<$aself as $trait<$field>>::Output,
-                               <$bself as $trait<$field>>::Output>;
-            #[inline]
-            fn $fn(self, a : $field) -> Self::Output {
-                Pair(maybe_ref!($refl, self.0).$fn(a),
-                     maybe_ref!($refl, self.1).$fn(a))
+        impl<'a, A, B, C, D> $trait<Pair<C, D>> for &'a Pair<A, B>
+        where
+            &'a A: $trait<C>,
+            &'a B: $trait<D>,
+        {
+            type Output = Pair<<&'a A as $trait<C>>::Output, <&'a B as $trait<D>>::Output>;
+            fn $fn(self, Pair(c, d): Pair<C, D>) -> Self::Output {
+                let Pair(ref a, ref b) = self;
+                Pair(a.$fn(c), b.$fn(d))
             }
         }
-    }
-}
 
-// Not used due to compiler overflow
-#[allow(unused_macros)]
-macro_rules! impl_scalarlhs_op {
-    (($a : ty, $b : ty), $field : ty, $trait:ident, $fn:ident, $refr:ident) => {
-        impl_scalarlhs_op!(@doit: $trait, $fn,
-                                  maybe_lifetime!($refr, &'r Pair<$a,$b>),
-                                  (maybe_lifetime!($refr, &'r $a),
-                                   maybe_lifetime!($refr, &'r $b));
-                                  $refr, $field);
-    };
-    (@doit: $trait:ident, $fn:ident,
-            $in:ty, ($ain:ty, $bin:ty);
-            $refr:ident, $field:ty) => {
-        impl<'r> $trait<$in>
-        for $field
-        where $field : $trait<$ain>
-                     + $trait<$bin> {
-            type Output = Pair<<$field as $trait<$ain>>::Output,
-                               <$field as $trait<$bin>>::Output>;
-            #[inline]
-            fn $fn(self, x : $in) -> Self::Output {
-                Pair(self.$fn(maybe_ref!($refr, x.0)),
-                     self.$fn(maybe_ref!($refr, x.1)))
+        impl<'a, 'b, A, B, C, D> $trait<&'b Pair<C, D>> for &'a Pair<A, B>
+        where
+            &'a A: $trait<&'b C>,
+            &'a B: $trait<&'b D>,
+        {
+            type Output = Pair<<&'a A as $trait<&'b C>>::Output, <&'a B as $trait<&'b D>>::Output>;
+            fn $fn(self, Pair(ref c, ref d): &'b Pair<C, D>) -> Self::Output {
+                let Pair(ref a, ref b) = self;
+                Pair(a.$fn(c), b.$fn(d))
+            }
+        }
+
+        impl<'b, A, B, C, D> $trait<&'b Pair<C, D>> for Pair<A, B>
+        where
+            A: $trait<&'b C>,
+            B: $trait<&'b D>,
+        {
+            type Output = Pair<<A as $trait<&'b C>>::Output, <B as $trait<&'b D>>::Output>;
+            fn $fn(self, Pair(ref c, ref d): &'b Pair<C, D>) -> Self::Output {
+                let Pair(a, b) = self;
+                Pair(a.$fn(c), b.$fn(d))
             }
         }
     };
 }
 
-macro_rules! impl_scalar_assignop {
-    (($a : ty, $b : ty), $field : ty, $trait : ident, $fn : ident) => {
-        impl<'r> $trait<$field> for Pair<$a, $b>
+impl_binary!(Add, add);
+impl_binary!(Sub, sub);
+
+macro_rules! impl_scalar {
+    ($trait:ident, $fn:ident) => {
+        impl<A, B, F: Num> $trait<F> for Pair<A, B>
+        where
+            A: $trait<F>,
+            B: $trait<F>,
+        {
+            type Output = Pair<A::Output, B::Output>;
+            fn $fn(self, t: F) -> Self::Output {
+                let Pair(a, b) = self;
+                Pair(a.$fn(t), b.$fn(t))
+            }
+        }
+
+        impl<'a, A, B, F: Num> $trait<F> for &'a Pair<A, B>
         where
-            $a: $trait<$field>,
-            $b: $trait<$field>,
+            &'a A: $trait<F>,
+            &'a B: $trait<F>,
+        {
+            type Output = Pair<<&'a A as $trait<F>>::Output, <&'a B as $trait<F>>::Output>;
+            fn $fn(self, t: F) -> Self::Output {
+                let Pair(ref a, ref b) = self;
+                Pair(a.$fn(t), b.$fn(t))
+            }
+        }
+
+        // impl<'a, 'b, A, B> $trait<&'b $F> for &'a Pair<A, B>
+        // where
+        //     &'a A: $trait<&'b $F>,
+        //     &'a B: $trait<&'b $F>,
+        // {
+        //     type Output =
+        //         Pair<<&'a A as $trait<&'b $F>>::Output, <&'a B as $trait<&'b $F>>::Output>;
+        //     fn $fn(self, t: &'b $F) -> Self::Output {
+        //         let Pair(ref a, ref b) = self;
+        //         Pair(a.$fn(t), b.$fn(t))
+        //     }
+        // }
+
+        // impl<'b, A, B> $trait<&'b $F> for Pair<A, B>
+        // where
+        //     A: $trait<&'b $F>,
+        //     B: $trait<&'b $F>,
+        // {
+        //     type Output = Pair<<A as $trait<&'b $F>>::Output, <B as $trait<&'b $F>>::Output>;
+        //     fn $fn(self, t: &'b $F) -> Self::Output {
+        //         let Pair(a, b) = self;
+        //         Pair(a.$fn(t), b.$fn(t))
+        //     }
+        // }
+    };
+}
+
+impl_scalar!(Mul, mul);
+impl_scalar!(Div, div);
+
+macro_rules! impl_scalar_lhs {
+    ($trait:ident, $fn:ident, $F:ty) => {
+        impl<A, B> $trait<Pair<A, B>> for $F
+        where
+            $F: $trait<A> + $trait<B>,
         {
-            #[inline]
-            fn $fn(&mut self, a: $field) -> () {
-                self.0.$fn(a);
-                self.1.$fn(a);
+            type Output = Pair<<$F as $trait<A>>::Output, <$F as $trait<B>>::Output>;
+            fn $fn(self, Pair(a, b): Pair<A, B>) -> Self::Output {
+                Pair(self.$fn(a), self.$fn(b))
+            }
+        }
+
+        // Compiler overflow:
+        //
+        // impl<'a, A, B> $trait<&'a Pair<A, B>> for $F
+        // where
+        //     $F: $trait<&'a A> + $trait<&'a B>,
+        // {
+        //     type Output = Pair<<$F as $trait<&'a A>>::Output, <$F as $trait<&'a B>>::Output>;
+        //     fn $fn(self, Pair(a, b): &'a Pair<A, B>) -> Self::Output {
+        //         Pair(self.$fn(a), self.$fn(b))
+        //     }
+        // }
+    };
+}
+
+impl_scalar_lhs!(Mul, mul, f32);
+impl_scalar_lhs!(Mul, mul, f64);
+impl_scalar_lhs!(Div, div, f32);
+impl_scalar_lhs!(Div, div, f64);
+
+macro_rules! impl_binary_mut {
+    ($trait:ident, $fn:ident) => {
+        impl<'a, A, B, C, D> $trait<Pair<C, D>> for Pair<A, B>
+        where
+            A: $trait<C>,
+            B: $trait<D>,
+        {
+            fn $fn(&mut self, Pair(c, d): Pair<C, D>) {
+                let Pair(ref mut a, ref mut b) = self;
+                a.$fn(c);
+                b.$fn(d);
+            }
+        }
+
+        impl<'a, 'b, A, B, C, D> $trait<&'b Pair<C, D>> for Pair<A, B>
+        where
+            A: $trait<&'b C>,
+            B: $trait<&'b D>,
+        {
+            fn $fn(&mut self, Pair(ref c, ref d): &'b Pair<C, D>) {
+                let Pair(ref mut a, ref mut b) = self;
+                a.$fn(c);
+                b.$fn(d);
             }
         }
     };
 }
 
-macro_rules! impl_unaryop {
-    (($a : ty, $b : ty), $trait:ident, $fn:ident, $refl:ident) => {
-        impl_unaryop!(@doit: $trait, $fn;
-                             maybe_lifetime!($refl, &'l Pair<$a,$b>),
-                             (maybe_lifetime!($refl, &'l $a),
-                              maybe_lifetime!($refl, &'l $b));
-                             $refl);
-    };
-    (@doit: $trait:ident, $fn:ident;
-            $self:ty, ($aself:ty, $bself:ty);
-            $refl : ident) => {
-        impl<'l> $trait
-        for $self
-        where $aself: $trait,
-              $bself: $trait {
-            type Output = Pair<<$aself as $trait>::Output,
-                               <$bself as $trait>::Output>;
-            #[inline]
-            fn $fn(self) -> Self::Output {
-                Pair(maybe_ref!($refl, self.0).$fn(),
-                     maybe_ref!($refl, self.1).$fn())
+impl_binary_mut!(AddAssign, add_assign);
+impl_binary_mut!(SubAssign, sub_assign);
+
+macro_rules! impl_scalar_mut {
+    ($trait:ident, $fn:ident) => {
+        impl<'a, A, B, F: Num> $trait<F> for Pair<A, B>
+        where
+            A: $trait<F>,
+            B: $trait<F>,
+        {
+            fn $fn(&mut self, t: F) {
+                let Pair(ref mut a, ref mut b) = self;
+                a.$fn(t);
+                b.$fn(t);
             }
         }
+    };
+}
+
+impl_scalar_mut!(MulAssign, mul_assign);
+impl_scalar_mut!(DivAssign, div_assign);
+
+/// Trait for ownable-by-consumption objects
+impl<A, B> Ownable for Pair<A, B>
+where
+    A: Ownable,
+    B: Ownable,
+{
+    type OwnedVariant = Pair<A::OwnedVariant, B::OwnedVariant>;
+
+    #[inline]
+    fn into_owned(self) -> Self::OwnedVariant {
+        Pair(self.0.into_owned(), self.1.into_owned())
+    }
+
+    #[inline]
+    fn clone_owned(&self) -> Self::OwnedVariant {
+        Pair(self.0.clone_owned(), self.1.clone_owned())
+    }
+
+    #[inline]
+    fn cow_owned<'b>(self) -> MyCow<'b, Self::OwnedVariant>
+    where
+        Self: 'b,
+    {
+        MyCow::Owned(self.into_owned())
+    }
+
+    #[inline]
+    fn ref_cow_owned<'b>(&'b self) -> MyCow<'b, Self::OwnedVariant>
+    where
+        Self: 'b,
+    {
+        MyCow::Owned(self.clone_owned())
     }
 }
 
-#[macro_export]
-macro_rules! impl_pair_vectorspace_ops {
-    (($a:ty, $b:ty), $field:ty) => {
-        impl_pair_vectorspace_ops!(@binary, ($a, $b), Add, add);
-        impl_pair_vectorspace_ops!(@binary, ($a, $b), Sub, sub);
-        impl_pair_vectorspace_ops!(@assign, ($a, $b), AddAssign, add_assign);
-        impl_pair_vectorspace_ops!(@assign, ($a, $b), SubAssign, sub_assign);
-        impl_pair_vectorspace_ops!(@scalar, ($a, $b), $field, Mul, mul);
-        impl_pair_vectorspace_ops!(@scalar, ($a, $b), $field, Div, div);
-        // Compiler overflow
-        // $(
-        //     impl_pair_vectorspace_ops!(@scalar_lhs, ($a, $b), $field, $impl_scalarlhs_op, Mul, mul);
-        // )*
-        impl_pair_vectorspace_ops!(@scalar_assign, ($a, $b), $field, MulAssign, mul_assign);
-        impl_pair_vectorspace_ops!(@scalar_assign, ($a, $b), $field, DivAssign, div_assign);
-        impl_pair_vectorspace_ops!(@unary, ($a, $b), Neg, neg);
-    };
-    (@binary, ($a : ty, $b : ty), $trait : ident, $fn : ident) => {
-        impl_binop!(($a, $b), $trait, $fn, ref, ref);
-        impl_binop!(($a, $b), $trait, $fn, ref, noref);
-        impl_binop!(($a, $b), $trait, $fn, noref, ref);
-        impl_binop!(($a, $b), $trait, $fn, noref, noref);
-    };
-    (@assign, ($a : ty, $b : ty), $trait : ident, $fn :ident) => {
-        impl_assignop!(($a, $b), $trait, $fn, ref);
-        impl_assignop!(($a, $b), $trait, $fn, noref);
-    };
-    (@scalar, ($a : ty, $b : ty), $field : ty, $trait : ident, $fn :ident) => {
-        impl_scalarop!(($a, $b), $field, $trait, $fn, ref);
-        impl_scalarop!(($a, $b), $field, $trait, $fn, noref);
-    };
-    (@scalar_lhs, ($a : ty, $b : ty), $field : ty, $trait : ident, $fn : ident) => {
-        impl_scalarlhs_op!(($a, $b), $field, $trait, $fn, ref);
-        impl_scalarlhs_op!(($a, $b), $field, $trait, $fn, noref);
-    };
-    (@scalar_assign, ($a : ty, $b : ty), $field : ty, $trait : ident, $fn : ident) => {
-        impl_scalar_assignop!(($a, $b), $field, $trait, $fn);
-    };
-    (@unary, ($a : ty, $b : ty), $trait : ident, $fn :  ident) => {
-        impl_unaryop!(($a, $b), $trait, $fn, ref);
-        impl_unaryop!(($a, $b), $trait, $fn, noref);
-    };
-}
-
-impl_pair_vectorspace_ops!((f32, f32), f32);
-impl_pair_vectorspace_ops!((f64, f64), f64);
-
-type PairOutput<F, A, B> = Pair<<A as Euclidean<F>>::Output, <B as Euclidean<F>>::Output>;
-
-impl<A, B, F> Euclidean<F> for Pair<A, B>
+/// We only support 'closed' `Euclidean` `Pair`s, as more general ones cause
+/// compiler overflows.
+impl<A, B, F: Float> Euclidean<F> for Pair<A, B>
 where
     A: Euclidean<F>,
     B: Euclidean<F>,
-    F: Float,
-    PairOutput<F, A, B>: Euclidean<F>,
-    Self: Sized
-        + Mul<F, Output = PairOutput<F, A, B>>
-        + MulAssign<F>
-        + Div<F, Output = PairOutput<F, A, B>>
-        + DivAssign<F>
-        + Add<Self, Output = PairOutput<F, A, B>>
-        + Sub<Self, Output = PairOutput<F, A, B>>
-        + for<'b> Add<&'b Self, Output = PairOutput<F, A, B>>
-        + for<'b> Sub<&'b Self, Output = PairOutput<F, A, B>>
-        + AddAssign<Self>
-        + for<'b> AddAssign<&'b Self>
-        + SubAssign<Self>
-        + for<'b> SubAssign<&'b Self>
-        + Neg<Output = PairOutput<F, A, B>>,
+    // //Pair<A, B>: Euclidean<F>,
+    // Self: Sized
+    //     + Mul<F, Output = Self::OwnedEuclidean>
+    //     + MulAssign<F>
+    //     + Div<F, Output = Self::OwnedEuclidean>
+    //     + DivAssign<F>
+    //     + Add<Self, Output = Self::OwnedEuclidean>
+    //     + Sub<Self, Output = Self::OwnedEuclidean>
+    //     + for<'b> Add<&'b Self, Output = Self::OwnedEuclidean>
+    //     + for<'b> Sub<&'b Self, Output = Self::OwnedEuclidean>
+    //     + AddAssign<Self>
+    //     + for<'b> AddAssign<&'b Self>
+    //     + SubAssign<Self>
+    //     + for<'b> SubAssign<&'b Self>
+    //     + Neg<Output = Self::OwnedEuclidean>,
 {
-    type Output = PairOutput<F, A, B>;
+    type PrincipalE = Pair<A::PrincipalE, B::PrincipalE>;
 
     fn dot<I: Instance<Self>>(&self, other: I) -> F {
-        let Pair(u, v) = other.decompose();
-        self.0.dot(u) + self.1.dot(v)
+        other.eval_decompose(|Pair(u, v)| self.0.dot(u) + self.1.dot(v))
     }
 
     fn norm2_squared(&self) -> F {
@@ -279,45 +330,59 @@
     }
 
     fn dist2_squared<I: Instance<Self>>(&self, other: I) -> F {
-        let Pair(u, v) = other.decompose();
-        self.0.dist2_squared(u) + self.1.dist2_squared(v)
+        other.eval_decompose(|Pair(u, v)| self.0.dist2_squared(u) + self.1.dist2_squared(v))
     }
 }
 
-impl<F, A, B, U, V> AXPY<F, Pair<U, V>> for Pair<A, B>
+impl<F, A, B> VectorSpace for Pair<A, B>
+where
+    A: VectorSpace<Field = F>,
+    B: VectorSpace<Field = F>,
+    F: Num,
+{
+    type Field = F;
+    type PrincipalV = Pair<A::PrincipalV, B::PrincipalV>;
+
+    /// Return a similar zero as `self`.
+    fn similar_origin(&self) -> Self::PrincipalV {
+        Pair(self.0.similar_origin(), self.1.similar_origin())
+    }
+
+    // #[inline]
+    // fn into_owned(self) -> Self::Owned {
+    //     Pair(self.0.into_owned(), self.1.into_owned())
+    // }
+}
+
+impl<F, A, B, U, V> AXPY<Pair<U, V>> for Pair<A, B>
 where
     U: Space,
     V: Space,
-    A: AXPY<F, U>,
-    B: AXPY<F, V>,
+    A: AXPY<U, Field = F>,
+    B: AXPY<V, Field = F>,
     F: Num,
-    Self: MulAssign<F>,
-    Pair<A, B>: MulAssign<F>,
-    Pair<A::Owned, B::Owned>: AXPY<F, Pair<U, V>>,
+    // Self: MulAssign<F> + DivAssign<F>,
+    // Pair<A, B>: MulAssign<F> + DivAssign<F>,
 {
-    type Owned = Pair<A::Owned, B::Owned>;
-
     fn axpy<I: Instance<Pair<U, V>>>(&mut self, α: F, x: I, β: F) {
-        let Pair(u, v) = x.decompose();
-        self.0.axpy(α, u, β);
-        self.1.axpy(α, v, β);
+        x.eval_decompose(|Pair(u, v)| {
+            self.0.axpy(α, u, β);
+            self.1.axpy(α, v, β);
+        })
     }
 
     fn copy_from<I: Instance<Pair<U, V>>>(&mut self, x: I) {
-        let Pair(u, v) = x.decompose();
-        self.0.copy_from(u);
-        self.1.copy_from(v);
+        x.eval_decompose(|Pair(u, v)| {
+            self.0.copy_from(u);
+            self.1.copy_from(v);
+        })
     }
 
     fn scale_from<I: Instance<Pair<U, V>>>(&mut self, α: F, x: I) {
-        let Pair(u, v) = x.decompose();
-        self.0.scale_from(α, u);
-        self.1.scale_from(α, v);
-    }
-
-    /// Return a similar zero as `self`.
-    fn similar_origin(&self) -> Self::Owned {
-        Pair(self.0.similar_origin(), self.1.similar_origin())
+        x.eval_decompose(|Pair(u, v)| {
+            self.0.scale_from(α, u);
+            self.1.scale_from(α, v);
+        })
     }
 
     /// Set self to zero.
@@ -332,6 +397,7 @@
 pub struct PairDecomposition<D, Q>(D, Q);
 
 impl<A: Space, B: Space> Space for Pair<A, B> {
+    type Principal = Pair<A::Principal, B::Principal>;
     type Decomp = PairDecomposition<A::Decomp, B::Decomp>;
 }
 
@@ -367,25 +433,16 @@
     V: Instance<B, Q>,
 {
     #[inline]
-    fn decompose<'b>(
-        self,
-    ) -> <PairDecomposition<D, Q> as Decomposition<Pair<A, B>>>::Decomposition<'b>
+    fn eval_ref<'b, R>(&'b self, f: impl FnOnce(Pair<D::Reference<'b>, Q::Reference<'b>>) -> R) -> R
     where
+        Pair<A, B>: 'b,
         Self: 'b,
-        Pair<A, B>: 'b,
     {
-        Pair(self.0.decompose(), self.1.decompose())
+        self.0.eval_ref(|a| self.1.eval_ref(|b| f(Pair(a, b))))
     }
 
     #[inline]
-    fn ref_instance(
-        &self,
-    ) -> <PairDecomposition<D, Q> as Decomposition<Pair<A, B>>>::Reference<'_> {
-        Pair(self.0.ref_instance(), self.1.ref_instance())
-    }
-
-    #[inline]
-    fn cow<'b>(self) -> MyCow<'b, Pair<A, B>>
+    fn cow<'b>(self) -> MyCow<'b, Pair<A::Principal, B::Principal>>
     where
         Self: 'b,
     {
@@ -393,9 +450,17 @@
     }
 
     #[inline]
-    fn own(self) -> Pair<A, B> {
+    fn own(self) -> Pair<A::Principal, B::Principal> {
         Pair(self.0.own(), self.1.own())
     }
+
+    #[inline]
+    fn decompose<'b>(self) -> Pair<D::Decomposition<'b>, Q::Decomposition<'b>>
+    where
+        Self: 'b,
+    {
+        Pair(self.0.decompose(), self.1.decompose())
+    }
 }
 
 impl<'a, A, B, U, V, D, Q> Instance<Pair<A, B>, PairDecomposition<D, Q>> for &'a Pair<U, V>
@@ -409,29 +474,16 @@
     &'a U: Instance<A, D>,
     &'a V: Instance<B, Q>,
 {
-    #[inline]
-    fn decompose<'b>(
-        self,
-    ) -> <PairDecomposition<D, Q> as Decomposition<Pair<A, B>>>::Decomposition<'b>
+    fn eval_ref<'b, R>(&'b self, f: impl FnOnce(Pair<D::Reference<'b>, Q::Reference<'b>>) -> R) -> R
     where
+        Pair<A, B>: 'b,
         Self: 'b,
-        Pair<A, B>: 'b,
     {
-        Pair(
-            D::lift(self.0.ref_instance()),
-            Q::lift(self.1.ref_instance()),
-        )
+        self.0.eval_ref(|a| self.1.eval_ref(|b| f(Pair(a, b))))
     }
 
     #[inline]
-    fn ref_instance(
-        &self,
-    ) -> <PairDecomposition<D, Q> as Decomposition<Pair<A, B>>>::Reference<'_> {
-        Pair(self.0.ref_instance(), self.1.ref_instance())
-    }
-
-    #[inline]
-    fn cow<'b>(self) -> MyCow<'b, Pair<A, B>>
+    fn cow<'b>(self) -> MyCow<'b, Pair<A::Principal, B::Principal>>
     where
         Self: 'b,
     {
@@ -439,10 +491,19 @@
     }
 
     #[inline]
-    fn own(self) -> Pair<A, B> {
+    fn own(self) -> Pair<A::Principal, B::Principal> {
         let Pair(ref u, ref v) = self;
         Pair(u.own(), v.own())
     }
+
+    #[inline]
+    fn decompose<'b>(self) -> Pair<D::Decomposition<'b>, Q::Decomposition<'b>>
+    where
+        Self: 'b,
+    {
+        let Pair(u, v) = self;
+        Pair(u.decompose(), v.decompose())
+    }
 }
 
 impl<A, B, D, Q> DecompositionMut<Pair<A, B>> for PairDecomposition<D, Q>
@@ -492,18 +553,63 @@
     }
 }
 
-impl<F, A, B, ExpA, ExpB, ExpJ> Norm<F, PairNorm<ExpA, ExpB, ExpJ>> for Pair<A, B>
+impl<F, A, B, ExpA, ExpB, ExpJ> Norm<PairNorm<ExpA, ExpB, ExpJ>, F> for Pair<A, B>
+where
+    F: Num,
+    ExpA: NormExponent,
+    ExpB: NormExponent,
+    ExpJ: NormExponent,
+    A: Norm<ExpA, F>,
+    B: Norm<ExpB, F>,
+    Loc<2, F>: Norm<ExpJ, F>,
+{
+    fn norm(&self, PairNorm(expa, expb, expj): PairNorm<ExpA, ExpB, ExpJ>) -> F {
+        Loc([self.0.norm(expa), self.1.norm(expb)]).norm(expj)
+    }
+}
+
+impl<F, A, B, ExpA, ExpB, ExpJ> Dist<PairNorm<ExpA, ExpB, ExpJ>, F> for Pair<A, B>
 where
     F: Num,
     ExpA: NormExponent,
     ExpB: NormExponent,
     ExpJ: NormExponent,
-    A: Norm<F, ExpA>,
-    B: Norm<F, ExpB>,
-    Loc<F, 2>: Norm<F, ExpJ>,
+    A: Dist<ExpA, F>,
+    B: Dist<ExpB, F>,
+    Loc<2, F>: Norm<ExpJ, F>,
 {
-    fn norm(&self, PairNorm(expa, expb, expj): PairNorm<ExpA, ExpB, ExpJ>) -> F {
-        Loc([self.0.norm(expa), self.1.norm(expb)]).norm(expj)
+    fn dist<I: Instance<Self>>(
+        &self,
+        x: I,
+        PairNorm(expa, expb, expj): PairNorm<ExpA, ExpB, ExpJ>,
+    ) -> F {
+        x.eval_decompose(|Pair(x1, x2)| {
+            Loc([self.0.dist(x1, expa), self.1.dist(x2, expb)]).norm(expj)
+        })
+    }
+}
+
+impl<F, A, B> Norm<L2, F> for Pair<A, B>
+where
+    F: Num,
+    A: Norm<L2, F>,
+    B: Norm<L2, F>,
+    Loc<2, F>: Norm<L2, F>,
+{
+    fn norm(&self, _: L2) -> F {
+        Loc([self.0.norm(L2), self.1.norm(L2)]).norm(L2)
+    }
+}
+
+impl<F, A, B> Dist<L2, F> for Pair<A, B>
+where
+    F: Num,
+    A: Dist<L2, F>,
+    B: Dist<L2, F>,
+    Loc<2, F>: Norm<L2, F>,
+{
+    fn dist<I: Instance<Self>>(&self, x: I, _: L2) -> F {
+        x.eval_decompose(|Pair(x1, x2)| Loc([self.0.dist(x1, L2), self.1.dist(x2, L2)]).norm(L2))
     }
 }
 
@@ -531,4 +637,72 @@
     B: HasDual<F>,
 {
     type DualSpace = Pair<A::DualSpace, B::DualSpace>;
+
+    fn dual_origin(&self) -> <Self::DualSpace as VectorSpace>::PrincipalV {
+        Pair(self.0.dual_origin(), self.1.dual_origin())
+    }
 }
+
+#[cfg(feature = "pyo3")]
+mod python {
+    use super::Pair;
+    use pyo3::conversion::FromPyObject;
+    use pyo3::types::{PyAny, PyTuple};
+    use pyo3::{Borrowed, Bound, IntoPyObject, PyErr, Python};
+
+    impl<'py, A, B> IntoPyObject<'py> for Pair<A, B>
+    where
+        A: IntoPyObject<'py>,
+        B: IntoPyObject<'py>,
+    {
+        type Target = PyTuple;
+        type Error = PyErr;
+        type Output = Bound<'py, Self::Target>;
+
+        fn into_pyobject(self, py: Python<'py>) -> Result<Self::Output, Self::Error> {
+            (self.0, self.1).into_pyobject(py)
+        }
+    }
+
+    /*
+    impl<'a, 'py, A, B> IntoPyObject<'py> for &'a mut Pair<A, B>
+    where
+        &'a mut A: IntoPyObject<'py>,
+        &'a mut B: IntoPyObject<'py>,
+    {
+        type Target = PyTuple;
+        type Error = PyErr;
+        type Output = Bound<'py, Self::Target>;
+
+        fn into_pyobject(self, py: Python<'py>) -> Result<Self::Output, Self::Error> {
+            (&mut self.0, &mut self.1).into_pyobject(py)
+        }
+    }
+
+    impl<'a, 'py, A, B> IntoPyObject<'py> for &'a Pair<A, B>
+    where
+        &'a A: IntoPyObject<'py>,
+        &'a B: IntoPyObject<'py>,
+    {
+        type Target = PyTuple;
+        type Error = PyErr;
+        type Output = Bound<'py, Self::Target>;
+
+        fn into_pyobject(self, py: Python<'py>) -> Result<Self::Output, Self::Error> {
+            (&self.0, &self.1).into_pyobject(py)
+        }
+    }
+    */
+
+    impl<'a, 'py, A, B> FromPyObject<'a, 'py> for Pair<A, B>
+    where
+        A: Clone + FromPyObject<'a, 'py>,
+        B: Clone + FromPyObject<'a, 'py>,
+    {
+        type Error = PyErr;
+
+        fn extract(ob: Borrowed<'a, 'py, PyAny>) -> Result<Self, PyErr> {
+            FromPyObject::extract(ob).map(|(a, b)| Pair(a, b))
+        }
+    }
+}