changeset
AXPY as LinSpace attempts, difficulties with Pairs and nalgebra
Fri, 20 Dec 2024 16:14:17 -0500
- author
- Tuomo Valkonen <tuomov@iki.fi>
- date
- Fri, 20 Dec 2024 16:14:17 -0500
- branch
- dev
- changeset 79
- d63e40672dd6
- parent 78
- cebedc4a8331
AXPY as LinSpace attempts, difficulties with Pairs and nalgebra
nalgebra should allow various storages, so InstanceMut as &self, but that won't work.
| src/direct_product.rs | file | annotate | diff | comparison | revisions | |
| src/linops.rs | file | annotate | diff | comparison | revisions | |
| src/loc.rs | file | annotate | diff | comparison | revisions | |
| src/nalgebra_support.rs | file | annotate | diff | comparison | revisions | |
| src/norms.rs | file | annotate | diff | comparison | revisions |
--- a/src/direct_product.rs Sat Dec 21 14:27:14 2024 -0500 +++ b/src/direct_product.rs Fri Dec 20 16:14:17 2024 -0500 @@ -279,28 +279,27 @@ } } -impl<F, A, B, U, V> AXPY<F, Pair<U, V>> for Pair<A, B> +impl<F, A, B> AXPY for Pair<A, B> where - U : Space, - V : Space, - A : AXPY<F, U>, - B : AXPY<F, V>, - F : Num + A : AXPY<Field=F>, + B : AXPY<Field=F>, + F : Num + AXPY { + type Field = F; - fn axpy<I : Instance<Pair<U,V>>>(&mut self, α : F, x : I, β : F) { + fn axpy<I : Instance<Pair<A,B>>>(&mut self, α : F, x : I, β : F) { let Pair(u, v) = x.decompose(); self.0.axpy(α, u, β); self.1.axpy(α, v, β); } - fn copy_from<I : Instance<Pair<U,V>>>(&mut self, x : I) { + fn copy_from<I : Instance<Pair<A,B>>>(&mut self, x : I) { let Pair(u, v) = x.decompose(); self.0.copy_from(u); self.1.copy_from(v); } - fn scale_from<I : Instance<Pair<U,V>>>(&mut self, α : F, x : I) { + fn scale_from<I : Instance<Pair<A,B>>>(&mut self, α : F, x : I) { let Pair(u, v) = x.decompose(); self.0.scale_from(α, u); self.1.scale_from(α, v);
--- a/src/linops.rs Sat Dec 21 14:27:14 2024 -0500 +++ b/src/linops.rs Fri Dec 20 16:14:17 2024 -0500 @@ -12,53 +12,66 @@ use crate::norms::{NormExponent, PairNorm, L1, L2, Linfinity, Norm, HasDual}; /// Trait for linear operators on `X`. -pub trait Linear<X : Space> : Mapping<X> -{ } +pub trait Linear<X : AXPY> : Mapping<X, Codomain = Self::LinCodomain> { + type LinCodomain : AXPY; // = Self::Codomain, but with AXPY enforced without trait bound bloat. +} /// Efficient in-place summation. -#[replace_float_literals(F::cast_from(literal))] -pub trait AXPY<F, X = Self> : Space -where - F : Num, - X : Space, -{ +#[replace_float_literals(Self::Field::cast_from(literal))] +pub trait AXPY : Space { + type Field : Num; + /// Computes `y = βy + αx`, where `y` is `Self`. - fn axpy<I : Instance<X>>(&mut self, α : F, x : I, β : F); + fn axpy<I : Instance<Self>>(&mut self, α : Self::Field, x : I, β : Self::Field); /// Copies `x` to `self`. - fn copy_from<I : Instance<X>>(&mut self, x : I) { + fn copy_from<I : Instance<Self>>(&mut self, x : I) { self.axpy(1.0, x, 0.0) } /// Computes `y = αx`, where `y` is `Self`. - fn scale_from<I : Instance<X>>(&mut self, α : F, x : I) { + fn scale_from<I : Instance<Self>>(&mut self, α : Self::Field, x : I) { self.axpy(α, x, 0.0) } } +macro_rules! impl_axpy { + ($($type:ty)*) => { $( + impl AXPY for $type { + type Field = $type; + fn axpy<I : Instance<Self>>(&mut self, α : $type, x : I, β : $type) { + *self = *self * α + x.own() * β; + } + } + )* }; +} + +impl_axpy!(i8 i16 i32 i64 i128 isize f32 f64); + + /// Efficient in-place application for [`Linear`] operators. -#[replace_float_literals(F::cast_from(literal))] -pub trait GEMV<F : Num, X : Space, Y = <Self as Mapping<X>>::Codomain> : Linear<X> { +#[replace_float_literals(X::Field::cast_from(literal))] +pub trait GEMV<X : AXPY> : Linear<X> { /// Computes `y = αAx + βy`, where `A` is `Self`. - fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F); + fn gemv<I : Instance<X>>(&self, y : &mut Self::LinCodomain, α : <Self::LinCodomain as AXPY>::Field, x : I, β : <Self::LinCodomain as AXPY>::Field); /// Computes `y = Ax`, where `A` is `Self` - fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){ + fn apply_mut<I : Instance<X>>(&self, y : &mut Self::LinCodomain, x : I){ self.gemv(y, 1.0, x, 0.0) } /// Computes `y += Ax`, where `A` is `Self` - fn apply_add<I : Instance<X>>(&self, y : &mut Y, x : I){ + fn apply_add<I : Instance<X>>(&self, y : &mut Self::LinCodomain, x : I){ self.gemv(y, 1.0, x, 1.0) } } /// Bounded linear operators -pub trait BoundedLinear<X, XExp, CodExp, F = f64> : Linear<X> +pub trait BoundedLinear<X, XExp, CodExp> : Linear<X> where F : Num, - X : Space + Norm<F, XExp>, + X : AXP + Norm<F, XExp>, XExp : NormExponent, CodExp : NormExponent { @@ -66,7 +79,7 @@ /// This is not expected to be the norm, just any bound on it that can be /// reasonably implemented. The [`NormExponent`] `xexp` indicates the norm /// in `X`, and `codexp` in the codomain. - fn opnorm_bound(&self, xexp : XExp, codexp : CodExp) -> F; + fn opnorm_bound(&self, xexp : XExp, codexp : CodExp) -> X::Field; // TODO: Should pick real subfield } // Linear operator application into mutable target. The [`AsRef`] bound @@ -105,27 +118,26 @@ /// /// We do not set further restrictions on the spacds, to allow preadjointing when `X` /// is on the dual space of `Ypre`, but a subset of it. -pub trait Preadjointable<X : Space, Ypre : Space> : Linear<X> +pub trait Preadjointable<X : AXPY, Ypre : AXPY> : Linear<X> //where // Ypre : HasDual<F, DualSpace=Self::Codomain>, { - type PreadjointCodomain : Space; + type PreadjointCodomain : AXPY; type Preadjoint<'a> : Linear<Ypre, Codomain=Self::PreadjointCodomain> where Self : 'a; /// Form the adjoint operator of `self`. fn preadjoint(&self) -> Self::Preadjoint<'_>; } - /// The identity operator #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)] -pub struct IdOp<X> (PhantomData<X>); +pub struct IdOp<X : AXPY> (PhantomData<X>); -impl<X> IdOp<X> { +impl<X : AXPY> IdOp<X> { pub fn new() -> IdOp<X> { IdOp(PhantomData) } } -impl<X : Clone + Space> Mapping<X> for IdOp<X> { +impl<X : Clone + AXPY> Mapping<X> for IdOp<X> { type Codomain = X; fn apply<I : Instance<X>>(&self, x : I) -> X { @@ -133,32 +145,29 @@ } } -impl<X : Clone + Space> Linear<X> for IdOp<X> -{ } +impl<X : Clone + AXPY> Linear<X> for IdOp<X> { + type LinCodomain = 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 + Space -{ +impl<X : Clone + AXPY> GEMV<X> for IdOp<X> { // Computes `y = αAx + βy`, where `A` is `Self`. - fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F) { + fn gemv<I : Instance<X>>(&self, y : &mut Self::LinCodomain, α : X::Field, x : I, β : X::Field) { y.axpy(α, x, β) } - fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){ + fn apply_mut<I : Instance<X>>(&self, y : &mut X, x : I){ y.copy_from(x); } } -impl<F, X, E> BoundedLinear<X, E, E, F> for IdOp<X> +impl<X, E> BoundedLinear<X, E, E> for IdOp<X> where - X : Space + Clone + Norm<F, E>, - F : Num, + X : AXPY + Clone + Norm<F, E>, + F : Num + AXPY, E : NormExponent { - fn opnorm_bound(&self, _xexp : E, _codexp : E) -> F { F::ONE } + fn opnorm_bound(&self, _xexp : E, _codexp : E) -> X::Field { F::ONE } } impl<X : Clone + HasDual<F, DualSpace=X>, F : Float> Adjointable<X, F> for IdOp<X> { @@ -168,7 +177,7 @@ fn adjoint(&self) -> Self::Adjoint<'_> { IdOp::new() } } -impl<X : Clone + Space> Preadjointable<X, X> for IdOp<X> { +impl<X : Clone + AXPY> Preadjointable<X, X> for IdOp<X> { type PreadjointCodomain = X; type Preadjoint<'a> = IdOp<X> where X : 'a; @@ -178,10 +187,12 @@ impl<S, T, E, X> Linear<X> for Composition<S, T, E> where - X : Space, + X : AXPY, T : Linear<X>, - S : Linear<T::Codomain> -{ } + S : Linear<T::LinCodomain> +{ + type LinCodomain = S::LinCodomain; +} impl<F, S, T, E, X, Y> GEMV<F, X, Y> for Composition<S, T, E> where @@ -227,17 +238,16 @@ use std::ops::Add; -impl<A, B, S, T> Mapping<Pair<A, B>> for RowOp<S, T> +impl<A, B, S, T, Y> Mapping<Pair<A, B>> for RowOp<S, T> where A : Space, B : Space, + Y : Space, S : Mapping<A>, T : Mapping<B>, - S::Codomain : Add<T::Codomain>, - <S::Codomain as Add<T::Codomain>>::Output : Space, - + S::Codomain : Add<T::Codomain, Output = Y>, { - type Codomain = <S::Codomain as Add<T::Codomain>>::Output; + type Codomain = Y; fn apply<I : Instance<Pair<A, B>>>(&self, x : I) -> Self::Codomain { let Pair(a, b) = x.decompose(); @@ -245,25 +255,30 @@ } } -impl<A, B, S, T> Linear<Pair<A, B>> for RowOp<S, T> +impl<F, A, B, S, T, Y> Linear<Pair<A, B>> for RowOp<S, T> where - A : Space, - B : Space, + F : Num + AXPY, + A : AXPY<Field=F>, + B : AXPY<Field=F>, + Y : AXPY, S : Linear<A>, T : Linear<B>, - S::Codomain : Add<T::Codomain>, - <S::Codomain as Add<T::Codomain>>::Output : Space, -{ } + S::LinCodomain : Add<T::LinCodomain, Output=Y>, +{ + type LinCodomain = Y; +} -impl<'b, F, S, T, Y, U, V> GEMV<F, Pair<U, V>, Y> for RowOp<S, T> +impl<'b, F, S, T, Y, G, U, V> GEMV<Pair<U, V>> for RowOp<S, T> where - U : Space, - V : Space, - S : GEMV<F, U, Y>, - T : GEMV<F, V, Y>, - F : Num, - Self : Linear<Pair<U, V>, Codomain=Y> + U : AXPY<Field=G>, + V : AXPY<Field=G>, + Y : AXPY<Field=F>, + S : GEMV<U>, + T : GEMV<V>, + F : Num + AXPY, + G : Num + AXPY, + S::LinCodomain : Add<T::LinCodomain, Output=Y>, { fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Y, α : F, x : I, β : F) { let Pair(u, v) = x.decompose(); @@ -288,33 +303,39 @@ /// “Column operator” $(S; T)$; $(S; T)x=(Sx, Tx)$. pub struct ColOp<S, T>(pub S, pub T); -impl<A, S, T> Mapping<A> for ColOp<S, T> +impl<X, S, T> Mapping<X> for ColOp<S, T> where - A : Space, - S : Mapping<A>, - T : Mapping<A>, + X : AXPY, + S : Mapping<X>, + T : Mapping<X>, { type Codomain = Pair<S::Codomain, T::Codomain>; - fn apply<I : Instance<A>>(&self, a : I) -> Self::Codomain { + fn apply<I : Instance<X>>(&self, a : I) -> Self::Codomain { Pair(self.0.apply(a.ref_instance()), self.1.apply(a)) } } -impl<A, S, T> Linear<A> for ColOp<S, T> +impl<X, S, T, A, B, F> Linear<X> for ColOp<S, T> where - A : Space, - S : Mapping<A>, - T : Mapping<A>, -{ } + F : Num + AXPY, + X : AXPY, + S : Linear<X, LinCodomain=A>, + T : Linear<X, LinCodomain=B>, + A : AXPY<Field=F>, + B : AXPY<Field=F>, +{ + type LinCodomain = Pair<S::LinCodomain, T::LinCodomain>; +} -impl<F, S, T, A, B, X> GEMV<F, X, Pair<A, B>> for ColOp<S, T> +impl<F, S, T, A, B, X> GEMV<X> for ColOp<S, T> where - X : Space, - S : GEMV<F, X, A>, - T : GEMV<F, X, B>, - F : Num, - Self : Linear<X, Codomain=Pair<A, B>> + F : Num + AXPY, + X : AXPY, + S : GEMV<X, LinCodomain=A>, + T : GEMV<X, LinCodomain=B>, + A : AXPY<Field=F>, + B : AXPY<Field=F>, { fn gemv<I : Instance<X>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) { self.0.gemv(&mut y.0, α, x.ref_instance(), β); @@ -341,7 +362,7 @@ B : HasDual<F>, S : Adjointable<A, F>, T : Adjointable<B, F>, - Yʹ : Space, + Yʹ : AXPY, S :: Codomain : HasDual<F, DualSpace=Yʹ>, T :: Codomain : HasDual<F, DualSpace=Yʹ>, S::Codomain : Add<T::Codomain>, @@ -361,15 +382,15 @@ } } -impl<A, B, Yʹ, S, T> Preadjointable<Pair<A,B>, Yʹ> for RowOp<S, T> -where - A : Space, - B : Space, - Yʹ : Space, +impl<F, A, B, Y, Yʹ, S, T> Preadjointable<Pair<A,B>, Yʹ> for RowOp<S, T> +wherea + A : AXPY, + B : AXPY, + Yʹ : AXPY, S : Preadjointable<A, Yʹ>, T : Preadjointable<B, Yʹ>, S::Codomain : Add<T::Codomain>, - <S::Codomain as Add<T::Codomain>>::Output : Space, + <S::Codomain as Add<T::Codomain>>::Output : AXPY, //Self : Linear<Pair<A, B>, Codomain=Y>, //for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<Yʹ, F>, { @@ -407,10 +428,10 @@ impl<A, Aʹ, Xʹ, Yʹ, S, T> Preadjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T> where - A : Space, - Aʹ : Space, - Xʹ : Space, - Yʹ : Space, + A : AXPY, + Aʹ : AXPY, + Xʹ : AXPY, + Yʹ : AXPY, S : Preadjointable<A, Xʹ, PreadjointCodomain=Aʹ>, T : Preadjointable<A, Yʹ, PreadjointCodomain=Aʹ>, Aʹ : ClosedAdd, @@ -429,8 +450,8 @@ impl<A, B, S, T> Mapping<Pair<A, B>> for DiagOp<S, T> where - A : Space, - B : Space, + A : AXPY, + B : AXPY, S : Mapping<A>, T : Mapping<B>, { @@ -442,26 +463,32 @@ } } -impl<A, B, S, T> Linear<Pair<A, B>> for DiagOp<S, T> -where - A : Space, - B : Space, - S : Linear<A>, - T : Linear<B>, -{ } - -impl<F, S, T, A, B, U, V> GEMV<F, Pair<U, V>, Pair<A, B>> for DiagOp<S, T> +impl<F, A, B, S, T, G, U, V> Linear<Pair<A, B>> for DiagOp<S, T> where - A : Space, - B : Space, - U : Space, - V : Space, - S : GEMV<F, U, A>, - T : GEMV<F, V, B>, - F : Num, - Self : Linear<Pair<U, V>, Codomain=Pair<A, B>>, + F : Num + AXPY, + G : Num + AXPY, + A : AXPY<Field=F>, + B : AXPY<Field=F>, + U : AXPY<Field=G>, + V : AXPY<Field=G>, + S : Linear<A, LinCodomain=U>, + T : Linear<B, LinCodomain=V>, { - fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) { + type LinCodomain = Pair<U, V>; +} + +impl<F, A, B, S, T, G, U, V> GEMV<Pair<U, V>> for DiagOp<S, T> +where + F : Num + AXPY, + G : Num + AXPY, + A : AXPY<Field=F>, + B : AXPY<Field=F>, + U : AXPY<Field=G>, + V : AXPY<Field=G>, + S : GEMV<A, LinCodomain=U>, + T : GEMV<B, LinCodomain=V>, +{ + fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, α : G, x : I, β : G) { let Pair(u, v) = x.decompose(); self.0.gemv(&mut y.0, α, u, β); self.1.gemv(&mut y.1, α, v, β); @@ -501,10 +528,10 @@ impl<A, B, Xʹ, Yʹ, S, T> Preadjointable<Pair<A,B>, Pair<Xʹ,Yʹ>> for DiagOp<S, T> where - A : Space, - B : Space, - Xʹ : Space, - Yʹ : Space, + A : AXPY, + B : AXPY, + Xʹ : AXPY, + Yʹ : AXPY, S : Preadjointable<A, Xʹ>, T : Preadjointable<B, Yʹ>, { @@ -523,16 +550,16 @@ macro_rules! pairnorm { ($expj:ty) => { impl<F, A, B, S, T, ExpA, ExpB, ExpR> - BoundedLinear<Pair<A, B>, PairNorm<ExpA, ExpB, $expj>, ExpR, F> + BoundedLinear<Pair<A, B>, PairNorm<ExpA, ExpB, $expj>, ExpR> for RowOp<S, T> where F : Float, - A : Space + Norm<F, ExpA>, - B : Space + Norm<F, ExpB>, + A : AXPY + Norm<F, ExpA>, + B : AXPY + Norm<F, ExpB>, S : BoundedLinear<A, ExpA, ExpR, F>, T : BoundedLinear<B, ExpB, ExpR, F>, S::Codomain : Add<T::Codomain>, - <S::Codomain as Add<T::Codomain>>::Output : Space, + <S::Codomain as Add<T::Codomain>>::Output : AXPY, ExpA : NormExponent, ExpB : NormExponent, ExpR : NormExponent, @@ -541,7 +568,7 @@ &self, PairNorm(expa, expb, _) : PairNorm<ExpA, ExpB, $expj>, expr : ExpR - ) -> F { + ) -> <Pair<A, B> as AXPY>::Field { // An application of the triangle inequality bounds the norm by the maximum // of the individual norms. A simple observation shows this to be exact. let na = self.0.opnorm_bound(expa, expr); @@ -550,12 +577,12 @@ } } - impl<F, A, S, T, ExpA, ExpS, ExpT> - BoundedLinear<A, ExpA, PairNorm<ExpS, ExpT, $expj>, F> + impl<A, S, T, ExpA, ExpS, ExpT> + BoundedLinear<A, ExpA, PairNorm<ExpS, ExpT, $expj>> for ColOp<S, T> where F : Float, - A : Space + Norm<F, ExpA>, + A : AXPY + Norm<F, ExpA>, S : BoundedLinear<A, ExpA, ExpS, F>, T : BoundedLinear<A, ExpA, ExpT, F>, ExpA : NormExponent, @@ -566,7 +593,7 @@ &self, expa : ExpA, PairNorm(exps, expt, _) : PairNorm<ExpS, ExpT, $expj> - ) -> F { + ) -> A::Field { // This is based on the rule for RowOp and ‖A^*‖ = ‖A‖, hence, // for A=[S; T], ‖A‖=‖[S^*, T^*]‖ ≤ max{‖S^*‖, ‖T^*‖} = max{‖S‖, ‖T‖} let ns = self.0.opnorm_bound(expa, exps);
--- a/src/loc.rs Sat Dec 21 14:27:14 2024 -0500 +++ b/src/loc.rs Fri Dec 20 16:14:17 2024 -0500 @@ -709,9 +709,12 @@ } } -impl<F : Num, const N : usize> Linear<Loc<F, N>> for Loc<F, N> { } +impl<F : Num + AXPY, const N : usize> Linear<Loc<F, N>> for Loc<F, N> { + type LinCodomain = F; +} -impl<F : Num, const N : usize> AXPY<F, Loc<F, N>> for Loc<F, N> { +impl<F : Num + AXPY, const N : usize> AXPY for Loc<F, N> { + type Field = F; #[inline] fn axpy<I : Instance<Loc<F, N>>>(&mut self, α : F, x : I, β : F) { x.eval(|x̃| {
--- a/src/nalgebra_support.rs Sat Dec 21 14:27:14 2024 -0500 +++ b/src/nalgebra_support.rs Fri Dec 20 16:14:17 2024 -0500 @@ -61,11 +61,13 @@ DefaultAllocator : Allocator<N,K>, DefaultAllocator : Allocator<M,K>, DefaultAllocator : Allocator<N,M>, - DefaultAllocator : Allocator<M,N> { + DefaultAllocator : Allocator<M,N> +{ + type LinCodomain = OMatrix<E,N,K>; } -impl<SM,SV1,SV2,N,M,K,E> GEMV<E, 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> + Clone, SV2: StorageMut<E,N,K>, +impl<SM,SV1,N,M,K,E> GEMV<Matrix<E,M,K,SV1>> for Matrix<E,N,M,SM> +where SM: Storage<E,N,M>, SV1: Storage<E,M,K> + Clone, N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + Float, DefaultAllocator : Allocator<N,K>, DefaultAllocator : Allocator<M,K>, @@ -85,18 +87,19 @@ } } -impl<SM,SV1,M,E> AXPY<E, Vector<E,M,SV1>> for Vector<E,M,SM> -where SM: StorageMut<E,M> + Clone, SV1: Storage<E,M> + Clone, +impl<SM,M,E> AXPY for Vector<E,M,SM> +where SM: StorageMut<E,M> + Clone, M : Dim, E : Scalar + Zero + One + Float, DefaultAllocator : Allocator<M> { + type Field = E; #[inline] - fn axpy<I : Instance<Vector<E,M,SV1>>>(&mut self, α : E, x : I, β : E) { + fn axpy<I : Instance<Vector<E,M,SM>>>(&mut self, α : E, x : I, β : E) { x.eval(|x̃| Matrix::axpy(self, α, x̃, β)) } #[inline] - fn copy_from<I : Instance<Vector<E,M,SV1>>>(&mut self, y : I) { + fn copy_from<I : Instance<Vector<E,M,SM>>>(&mut self, y : I) { y.eval(|ỹ| Matrix::copy_from(self, ỹ)) } }
--- a/src/norms.rs Sat Dec 21 14:27:14 2024 -0500 +++ b/src/norms.rs Fri Dec 20 16:14:17 2024 -0500 @@ -7,6 +7,7 @@ use crate::types::*; use crate::euclidean::*; use crate::mapping::{Mapping, Space, Instance}; +use crate::linops::AXPY; // // Abstract norms @@ -202,7 +203,7 @@ } } -pub trait Normed<F : Num = f64> : Space + Norm<F, Self::NormExp> { +pub trait Normed<F : Num = f64> : AXPY + Norm<F, Self::NormExp> { type NormExp : NormExponent; fn norm_exponent(&self) -> Self::NormExp;
