diff --git a/src/base/construction.rs b/src/base/construction.rs
index ae8c10d4..bb0a34c8 100644
--- a/src/base/construction.rs
+++ b/src/base/construction.rs
@@ -739,7 +739,7 @@ impl_constructors_from_data!(data; Dynamic, Dynamic;
  */
 impl<N, R: DimName, C: DimName> Zero for MatrixMN<N, R, C>
 where
-    N: Scalar + Zero + ClosedAdd,
+    N: Scalar + Zero,
     DefaultAllocator: Allocator<N, R, C>,
 {
     #[inline]
diff --git a/src/base/matrix.rs b/src/base/matrix.rs
index 783f8c19..5136eabe 100644
--- a/src/base/matrix.rs
+++ b/src/base/matrix.rs
@@ -16,7 +16,7 @@ use serde::{Deserialize, Deserializer, Serialize, Serializer};
 #[cfg(feature = "abomonation-serialize")]
 use abomonation::Abomonation;
 
-use simba::scalar::{ClosedAdd, ClosedMul, ClosedSub, Field, RealField};
+use simba::scalar::{ClosedAdd, ClosedMul, Field, RealField};
 use simba::simd::SimdPartialOrd;
 
 use crate::base::allocator::{Allocator, SameShapeAllocator, SameShapeC, SameShapeR};
@@ -28,7 +28,7 @@ use crate::base::iter::{
 use crate::base::storage::{
     ContiguousStorage, ContiguousStorageMut, Owned, SameShapeStorage, Storage, StorageMut,
 };
-use crate::base::{DefaultAllocator, MatrixMN, MatrixN, Scalar, Unit, VectorN};
+use crate::base::{DefaultAllocator, MatrixMN, MatrixN, Scalar, Unit, VectorN, ops::SimpleSub};
 use crate::SimdComplexField;
 
 /// A square matrix.
@@ -1572,7 +1572,7 @@ fn lower_exp() {
     )
 }
 
-impl<N: Scalar + ClosedAdd + ClosedSub + ClosedMul, R: Dim, C: Dim, S: Storage<N, R, C>>
+impl<N: Scalar + SimpleSub + ClosedMul, R: Dim, C: Dim, S: Storage<N, R, C>>
     Matrix<N, R, C, S>
 {
     /// The perpendicular product between two 2D column vectors, i.e. `a.x * b.y - a.y * b.x`.
@@ -1715,7 +1715,7 @@ impl<N: SimdComplexField, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S
     }
 }
 
-impl<N: Scalar + Zero + One + ClosedAdd + ClosedSub + ClosedMul, D: Dim, S: Storage<N, D>>
+impl<N: Scalar + Zero + One + ClosedAdd + SimpleSub + ClosedMul, D: Dim, S: Storage<N, D>>
     Vector<N, D, S>
 {
     /// Returns `self * (1.0 - t) + rhs * t`, i.e., the linear blend of the vectors x and y using the scalar value a.
diff --git a/src/base/matrix_alga.rs b/src/base/matrix_alga.rs
index 18182382..08908582 100644
--- a/src/base/matrix_alga.rs
+++ b/src/base/matrix_alga.rs
@@ -5,10 +5,11 @@ use num::{One, Zero};
 
 use alga::general::{
     AbstractGroup, AbstractGroupAbelian, AbstractLoop, AbstractMagma, AbstractModule,
-    AbstractMonoid, AbstractQuasigroup, AbstractSemigroup, Additive, ClosedAdd, ClosedMul,
-    ClosedNeg, ComplexField, Field, Identity, JoinSemilattice, Lattice, MeetSemilattice, Module,
+    AbstractMonoid, AbstractQuasigroup, AbstractSemigroup, Additive, ClosedMul, ClosedNeg,
+    ComplexField, Field, Identity, JoinSemilattice, Lattice, MeetSemilattice, Module,
     Multiplicative, RingCommutative, TwoSidedInverse,
 };
+use crate::base::{SimpleAdd, SimpleSub};
 use alga::linear::{
     FiniteDimInnerSpace, FiniteDimVectorSpace, InnerSpace, NormedSpace, VectorSpace,
 };
@@ -36,7 +37,7 @@ where
 
 impl<N, R: DimName, C: DimName> AbstractMagma<Additive> for MatrixMN<N, R, C>
 where
-    N: Scalar + ClosedAdd,
+    N: Scalar + SimpleAdd,
     DefaultAllocator: Allocator<N, R, C>,
 {
     #[inline]
@@ -71,12 +72,12 @@ macro_rules! inherit_additive_structure(
 );
 
 inherit_additive_structure!(
-    AbstractSemigroup<Additive> + ClosedAdd,
-    AbstractMonoid<Additive> + Zero + ClosedAdd,
-    AbstractQuasigroup<Additive> + ClosedAdd + ClosedNeg,
-    AbstractLoop<Additive> + Zero + ClosedAdd + ClosedNeg,
-    AbstractGroup<Additive> + Zero + ClosedAdd + ClosedNeg,
-    AbstractGroupAbelian<Additive> + Zero + ClosedAdd + ClosedNeg
+    AbstractSemigroup<Additive> + SimpleAdd,
+    AbstractMonoid<Additive> + Zero + SimpleAdd,
+    AbstractQuasigroup<Additive> + SimpleAdd + ClosedNeg,
+    AbstractLoop<Additive> + Zero + SimpleAdd + ClosedNeg,
+    AbstractGroup<Additive> + Zero + SimpleAdd + ClosedNeg,
+    AbstractGroupAbelian<Additive> + Zero + SimpleAdd + ClosedNeg
 );
 
 impl<N, R: DimName, C: DimName> AbstractModule for MatrixMN<N, R, C>
@@ -371,7 +372,7 @@ where
 
 impl<N, D: DimName> AbstractMagma<Multiplicative> for MatrixN<N, D>
 where
-    N: Scalar + Zero + One + ClosedAdd + ClosedMul,
+    N: Scalar + Zero + One + ClosedMul,
     DefaultAllocator: Allocator<N, D, D>,
 {
     #[inline]
@@ -383,7 +384,7 @@ where
 macro_rules! impl_multiplicative_structure(
     ($($marker: ident<$operator: ident> $(+ $bounds: ident)*),* $(,)*) => {$(
         impl<N, D: DimName> $marker<$operator> for MatrixN<N, D>
-            where N: Scalar + Zero + One + ClosedAdd + ClosedMul + $marker<$operator> $(+ $bounds)*,
+            where N: Scalar + Zero + One + ClosedMul + $marker<$operator> $(+ $bounds)*,
                   DefaultAllocator: Allocator<N, D, D> { }
     )*}
 );
diff --git a/src/base/mod.rs b/src/base/mod.rs
index 0e8ea72d..575076f8 100644
--- a/src/base/mod.rs
+++ b/src/base/mod.rs
@@ -8,6 +8,7 @@ pub mod default_allocator;
 pub mod dimension;
 pub mod iter;
 mod ops;
+pub use ops::{SimpleAdd, SimpleSub};
 pub mod storage;
 
 mod alias;
diff --git a/src/base/ops.rs b/src/base/ops.rs
index 12b26d1a..66ab1b88 100644
--- a/src/base/ops.rs
+++ b/src/base/ops.rs
@@ -132,9 +132,9 @@ where
  */
 
 macro_rules! componentwise_binop_impl(
-    ($Trait: ident, $method: ident, $bound: ident;
+    ($Trait: ident, $method: ident, $bound: ident, $bound_assign: ident;
      $TraitAssign: ident, $method_assign: ident, $method_assign_statically_unchecked: ident,
-     $method_assign_statically_unchecked_rhs: ident;
+     $method_assign_statically_unchecked_rhs: ident, $method_assign_statically_unchecked_lhs: ident;
      $method_to: ident, $method_to_statically_unchecked: ident) => {
 
         impl<N, R1: Dim, C1: Dim, SA: Storage<N, R1, C1>> Matrix<N, R1, C1, SA>
@@ -182,34 +182,24 @@ macro_rules! componentwise_binop_impl(
             }
 
 
+            /*
+             *
+             * Methods without dimension checking at compile-time.
+             * This is useful for code reuse because the sum representative system does not plays
+             * easily with static checks.
+             *
+             */
+            /// Equivalent to `self + rhs` but stores the result into `out` to avoid allocations.
             #[inline]
-            fn $method_assign_statically_unchecked<R2, C2, SB>(&mut self, rhs: &Matrix<N, R2, C2, SB>)
-                where R2: Dim,
-                      C2: Dim,
-                      SA: StorageMut<N, R1, C1>,
-                      SB: Storage<N, R2, C2> {
-                assert!(self.shape() == rhs.shape(), "Matrix addition/subtraction dimensions mismatch.");
-
-                // This is the most common case and should be deduced at compile-time.
-                // FIXME: use specialization instead?
-                if self.data.is_contiguous() && rhs.data.is_contiguous() {
-                    let arr1 = self.data.as_mut_slice();
-                    let arr2 = rhs.data.as_slice();
-                    for i in 0 .. arr2.len() {
-                        unsafe {
-                            arr1.get_unchecked_mut(i).$method_assign(arr2.get_unchecked(i).inlined_clone());
-                        }
-                    }
-                }
-                else {
-                    for j in 0 .. rhs.ncols() {
-                        for i in 0 .. rhs.nrows() {
-                            unsafe {
-                                self.get_unchecked_mut((i, j)).$method_assign(rhs.get_unchecked((i, j)).inlined_clone())
-                            }
-                        }
-                    }
-                }
+            pub fn $method_to<R2: Dim, C2: Dim, SB,
+                              R3: Dim, C3: Dim, SC>(&self,
+                                                    rhs: &Matrix<N, R2, C2, SB>,
+                                                    out: &mut Matrix<N, R3, C3, SC>)
+                where SB: Storage<N, R2, C2>,
+                      SC: StorageMut<N, R3, C3>,
+                      ShapeConstraint: SameNumberOfRows<R1, R2> + SameNumberOfColumns<C1, C2> +
+                                       SameNumberOfRows<R1, R3> + SameNumberOfColumns<C1, C3> {
+                self.$method_to_statically_unchecked(rhs, out)
             }
 
 
@@ -245,24 +235,71 @@ macro_rules! componentwise_binop_impl(
             }
 
 
-            /*
-             *
-             * Methods without dimension checking at compile-time.
-             * This is useful for code reuse because the sum representative system does not plays
-             * easily with static checks.
-             *
-             */
-            /// Equivalent to `self + rhs` but stores the result into `out` to avoid allocations.
             #[inline]
-            pub fn $method_to<R2: Dim, C2: Dim, SB,
-                              R3: Dim, C3: Dim, SC>(&self,
-                                                    rhs: &Matrix<N, R2, C2, SB>,
-                                                    out: &mut Matrix<N, R3, C3, SC>)
-                where SB: Storage<N, R2, C2>,
-                      SC: StorageMut<N, R3, C3>,
-                      ShapeConstraint: SameNumberOfRows<R1, R2> + SameNumberOfColumns<C1, C2> +
-                                       SameNumberOfRows<R1, R3> + SameNumberOfColumns<C1, C3> {
-                self.$method_to_statically_unchecked(rhs, out)
+            fn $method_assign_statically_unchecked_lhs<R2, C2, SB>(&mut self, rhs: &Matrix<N, R2, C2, SB>)
+                where R2: Dim,
+                      C2: Dim,
+                      SA: StorageMut<N, R1, C1>,
+                      SB: Storage<N, R2, C2> {
+                assert!(self.shape() == rhs.shape(), "Matrix addition/subtraction dimensions mismatch.");
+
+                // This is the most common case and should be deduced at compile-time.
+                // FIXME: use specialization instead?
+                if self.data.is_contiguous() && rhs.data.is_contiguous() {
+                    let arr1 = self.data.as_mut_slice();
+                    let arr2 = rhs.data.as_slice();
+                    for i in 0 .. arr2.len() {
+                        unsafe {
+                            let res = arr1.get_unchecked(i).inlined_clone().$method(arr2.get_unchecked(i).inlined_clone());
+                            *arr1.get_unchecked_mut(i) = res;
+                        }
+                    }
+                }
+                else {
+                    for j in 0 .. rhs.ncols() {
+                        for i in 0 .. rhs.nrows() {
+                            unsafe {
+                                let r = self.get_unchecked_mut((i, j));
+                                *r = r.inlined_clone().$method(r.inlined_clone())
+                            }
+                        }
+                    }
+                }
+            }
+        }
+
+
+        impl<N, R1: Dim, C1: Dim, SA: Storage<N, R1, C1>> Matrix<N, R1, C1, SA>
+            where N: Scalar + $bound_assign {
+
+            #[inline]
+            fn $method_assign_statically_unchecked<R2, C2, SB>(&mut self, rhs: &Matrix<N, R2, C2, SB>)
+                where R2: Dim,
+                      C2: Dim,
+                      SA: StorageMut<N, R1, C1>,
+                      SB: Storage<N, R2, C2> {
+                assert!(self.shape() == rhs.shape(), "Matrix addition/subtraction dimensions mismatch.");
+
+                // This is the most common case and should be deduced at compile-time.
+                // FIXME: use specialization instead?
+                if self.data.is_contiguous() && rhs.data.is_contiguous() {
+                    let arr1 = self.data.as_mut_slice();
+                    let arr2 = rhs.data.as_slice();
+                    for i in 0 .. arr2.len() {
+                        unsafe {
+                            arr1.get_unchecked_mut(i).$method_assign(arr2.get_unchecked(i).inlined_clone());
+                        }
+                    }
+                }
+                else {
+                    for j in 0 .. rhs.ncols() {
+                        for i in 0 .. rhs.nrows() {
+                            unsafe {
+                                self.get_unchecked_mut((i, j)).$method_assign(rhs.get_unchecked((i, j)).inlined_clone())
+                            }
+                        }
+                    }
+                }
             }
         }
 
@@ -279,7 +316,7 @@ macro_rules! componentwise_binop_impl(
             fn $method(self, rhs: &'b Matrix<N, R2, C2, SB>) -> Self::Output {
                 assert!(self.shape() == rhs.shape(), "Matrix addition/subtraction dimensions mismatch.");
                 let mut res = self.into_owned_sum::<R2, C2>();
-                res.$method_assign_statically_unchecked(rhs);
+                res.$method_assign_statically_unchecked_lhs(rhs);
                 res
             }
         }
@@ -342,7 +379,7 @@ macro_rules! componentwise_binop_impl(
 
         impl<'b, N, R1, C1, R2, C2, SA, SB> $TraitAssign<&'b Matrix<N, R2, C2, SB>> for Matrix<N, R1, C1, SA>
             where R1: Dim, C1: Dim, R2: Dim, C2: Dim,
-                  N: Scalar + $bound,
+                  N: Scalar + $bound_assign,
                   SA: StorageMut<N, R1, C1>,
                   SB: Storage<N, R2, C2>,
                   ShapeConstraint: SameNumberOfRows<R1, R2> + SameNumberOfColumns<C1, C2> {
@@ -355,7 +392,7 @@ macro_rules! componentwise_binop_impl(
 
         impl<N, R1, C1, R2, C2, SA, SB> $TraitAssign<Matrix<N, R2, C2, SB>> for Matrix<N, R1, C1, SA>
             where R1: Dim, C1: Dim, R2: Dim, C2: Dim,
-                  N: Scalar + $bound,
+                  N: Scalar + $bound_assign,
                   SA: StorageMut<N, R1, C1>,
                   SB: Storage<N, R2, C2>,
                   ShapeConstraint: SameNumberOfRows<R1, R2> + SameNumberOfColumns<C1, C2> {
@@ -368,16 +405,24 @@ macro_rules! componentwise_binop_impl(
     }
 );
 
-componentwise_binop_impl!(Add, add, ClosedAdd;
-                          AddAssign, add_assign, add_assign_statically_unchecked, add_assign_statically_unchecked_mut;
+/// Trait __alias__ for `Add` with result of type `Self`.
+pub trait SimpleAdd<Right = Self>: Sized + Add<Right, Output = Self> {}
+impl<T, Right> SimpleAdd<Right> for T where T: Add<Right, Output = T> {}
+
+/// Trait __alias__ for `Sub` with result of type `Self`.
+pub trait SimpleSub<Right = Self>: Sized + Sub<Right, Output = Self> {}
+impl<T, Right> SimpleSub<Right> for T where T: Sub<Right, Output = T> {}
+
+componentwise_binop_impl!(Add, add, SimpleAdd, ClosedAdd;
+                          AddAssign, add_assign, add_assign_statically_unchecked, add_assign_statically_unchecked_mut, add_assign_statically_unchecked_lhs;
                           add_to, add_to_statically_unchecked);
-componentwise_binop_impl!(Sub, sub, ClosedSub;
-                          SubAssign, sub_assign, sub_assign_statically_unchecked, sub_assign_statically_unchecked_mut;
+componentwise_binop_impl!(Sub, sub, SimpleSub, ClosedSub;
+                          SubAssign, sub_assign, sub_assign_statically_unchecked, sub_assign_statically_unchecked_mut, sub_assign_statically_unchecked_lhs;
                           sub_to, sub_to_statically_unchecked);
 
 impl<N, R: DimName, C: DimName> iter::Sum for MatrixMN<N, R, C>
 where
-    N: Scalar + ClosedAdd + Zero,
+    N: Scalar + SimpleAdd + Zero,
     DefaultAllocator: Allocator<N, R, C>,
 {
     fn sum<I: Iterator<Item = MatrixMN<N, R, C>>>(iter: I) -> MatrixMN<N, R, C> {
@@ -387,7 +432,7 @@ where
 
 impl<N, C: Dim> iter::Sum for MatrixMN<N, Dynamic, C>
 where
-    N: Scalar + ClosedAdd + Zero,
+    N: Scalar + SimpleAdd + Zero,
     DefaultAllocator: Allocator<N, Dynamic, C>,
 {
     /// # Example
@@ -417,7 +462,7 @@ where
 
 impl<'a, N, R: DimName, C: DimName> iter::Sum<&'a MatrixMN<N, R, C>> for MatrixMN<N, R, C>
 where
-    N: Scalar + ClosedAdd + Zero,
+    N: Scalar + SimpleAdd + Zero,
     DefaultAllocator: Allocator<N, R, C>,
 {
     fn sum<I: Iterator<Item = &'a MatrixMN<N, R, C>>>(iter: I) -> MatrixMN<N, R, C> {
@@ -427,7 +472,7 @@ where
 
 impl<'a, N, C: Dim> iter::Sum<&'a MatrixMN<N, Dynamic, C>> for MatrixMN<N, Dynamic, C>
 where
-    N: Scalar + ClosedAdd + Zero,
+    N: Scalar + SimpleAdd + Zero,
     DefaultAllocator: Allocator<N, Dynamic, C>,
 {
     /// # Example
diff --git a/src/geometry/point_ops.rs b/src/geometry/point_ops.rs
index e6547fc4..0676f15d 100644
--- a/src/geometry/point_ops.rs
+++ b/src/geometry/point_ops.rs
@@ -77,65 +77,66 @@ where
  *
  */
 
+use crate::base::{SimpleAdd, SimpleSub};
 // Point - Point
-add_sub_impl!(Sub, sub, ClosedSub;
+add_sub_impl!(Sub, sub, SimpleSub;
     (D, U1), (D, U1) for D: DimName;
     self: &'a Point<N, D>, right: &'b Point<N, D>, Output = VectorSum<N, D, D>;
     &self.coords - &right.coords; 'a, 'b);
 
-add_sub_impl!(Sub, sub, ClosedSub;
+add_sub_impl!(Sub, sub, SimpleSub;
     (D, U1), (D, U1) for D: DimName;
     self: &'a Point<N, D>, right: Point<N, D>, Output = VectorSum<N, D, D>;
     &self.coords - right.coords; 'a);
 
-add_sub_impl!(Sub, sub, ClosedSub;
+add_sub_impl!(Sub, sub, SimpleSub;
     (D, U1), (D, U1) for D: DimName;
     self: Point<N, D>, right: &'b Point<N, D>, Output = VectorSum<N, D, D>;
     self.coords - &right.coords; 'b);
 
-add_sub_impl!(Sub, sub, ClosedSub;
+add_sub_impl!(Sub, sub, SimpleSub;
     (D, U1), (D, U1) for D: DimName;
     self: Point<N, D>, right: Point<N, D>, Output = VectorSum<N, D, D>;
     self.coords - right.coords; );
 
 // Point - Vector
-add_sub_impl!(Sub, sub, ClosedSub;
+add_sub_impl!(Sub, sub, SimpleSub;
     (D1, U1), (D2, U1) -> (D1) for D1: DimName, D2: Dim, SB: Storage<N, D2>;
     self: &'a Point<N, D1>, right: &'b Vector<N, D2, SB>, Output = Point<N, D1>;
     Self::Output::from(&self.coords - right); 'a, 'b);
 
-add_sub_impl!(Sub, sub, ClosedSub;
+add_sub_impl!(Sub, sub, SimpleSub;
     (D1, U1), (D2, U1) -> (D1) for D1: DimName, D2: Dim, SB: Storage<N, D2>;
     self: &'a Point<N, D1>, right: Vector<N, D2, SB>, Output = Point<N, D1>;
     Self::Output::from(&self.coords - &right); 'a); // FIXME: should not be a ref to `right`.
 
-add_sub_impl!(Sub, sub, ClosedSub;
+add_sub_impl!(Sub, sub, SimpleSub;
     (D1, U1), (D2, U1) -> (D1) for D1: DimName, D2: Dim, SB: Storage<N, D2>;
     self: Point<N, D1>, right: &'b Vector<N, D2, SB>, Output = Point<N, D1>;
     Self::Output::from(self.coords - right); 'b);
 
-add_sub_impl!(Sub, sub, ClosedSub;
+add_sub_impl!(Sub, sub, SimpleSub;
     (D1, U1), (D2, U1) -> (D1) for D1: DimName, D2: Dim, SB: Storage<N, D2>;
     self: Point<N, D1>, right: Vector<N, D2, SB>, Output = Point<N, D1>;
     Self::Output::from(self.coords - right); );
 
 // Point + Vector
-add_sub_impl!(Add, add, ClosedAdd;
+add_sub_impl!(Add, add, SimpleAdd;
     (D1, U1), (D2, U1) -> (D1) for D1: DimName, D2: Dim, SB: Storage<N, D2>;
     self: &'a Point<N, D1>, right: &'b Vector<N, D2, SB>, Output = Point<N, D1>;
     Self::Output::from(&self.coords + right); 'a, 'b);
 
-add_sub_impl!(Add, add, ClosedAdd;
+add_sub_impl!(Add, add, SimpleAdd;
     (D1, U1), (D2, U1) -> (D1) for D1: DimName, D2: Dim, SB: Storage<N, D2>;
     self: &'a Point<N, D1>, right: Vector<N, D2, SB>, Output = Point<N, D1>;
     Self::Output::from(&self.coords + &right); 'a); // FIXME: should not be a ref to `right`.
 
-add_sub_impl!(Add, add, ClosedAdd;
+add_sub_impl!(Add, add, SimpleAdd;
     (D1, U1), (D2, U1) -> (D1) for D1: DimName, D2: Dim, SB: Storage<N, D2>;
     self: Point<N, D1>, right: &'b Vector<N, D2, SB>, Output = Point<N, D1>;
     Self::Output::from(self.coords + right); 'b);
 
-add_sub_impl!(Add, add, ClosedAdd;
+add_sub_impl!(Add, add, SimpleAdd;
     (D1, U1), (D2, U1) -> (D1) for D1: DimName, D2: Dim, SB: Storage<N, D2>;
     self: Point<N, D1>, right: Vector<N, D2, SB>, Output = Point<N, D1>;
     Self::Output::from(self.coords + right); );
diff --git a/src/geometry/translation_ops.rs b/src/geometry/translation_ops.rs
index c00b25bc..f0856724 100644
--- a/src/geometry/translation_ops.rs
+++ b/src/geometry/translation_ops.rs
@@ -9,45 +9,46 @@ use crate::base::{DefaultAllocator, Scalar};
 
 use crate::geometry::{Point, Translation};
 
+use crate::base::{SimpleAdd, SimpleSub};
 // Translation × Translation
-add_sub_impl!(Mul, mul, ClosedAdd;
+add_sub_impl!(Mul, mul, SimpleAdd;
     (D, U1), (D, U1) -> (D) for D: DimName;
     self: &'a Translation<N, D>, right: &'b Translation<N, D>, Output = Translation<N, D>;
     Translation::from(&self.vector + &right.vector); 'a, 'b);
 
-add_sub_impl!(Mul, mul, ClosedAdd;
+add_sub_impl!(Mul, mul, SimpleAdd;
     (D, U1), (D, U1) -> (D) for D: DimName;
     self: &'a Translation<N, D>, right: Translation<N, D>, Output = Translation<N, D>;
     Translation::from(&self.vector + right.vector); 'a);
 
-add_sub_impl!(Mul, mul, ClosedAdd;
+add_sub_impl!(Mul, mul, SimpleAdd;
     (D, U1), (D, U1) -> (D) for D: DimName;
     self: Translation<N, D>, right: &'b Translation<N, D>, Output = Translation<N, D>;
     Translation::from(self.vector + &right.vector); 'b);
 
-add_sub_impl!(Mul, mul, ClosedAdd;
+add_sub_impl!(Mul, mul, SimpleAdd;
     (D, U1), (D, U1) -> (D) for D: DimName;
     self: Translation<N, D>, right: Translation<N, D>, Output = Translation<N, D>;
     Translation::from(self.vector + right.vector); );
 
 // Translation ÷ Translation
 // FIXME: instead of calling inverse explicitly, could we just add a `mul_tr` or `mul_inv` method?
-add_sub_impl!(Div, div, ClosedSub;
+add_sub_impl!(Div, div, SimpleSub;
     (D, U1), (D, U1) -> (D) for D: DimName;
     self: &'a Translation<N, D>, right: &'b Translation<N, D>, Output = Translation<N, D>;
     Translation::from(&self.vector - &right.vector); 'a, 'b);
 
-add_sub_impl!(Div, div, ClosedSub;
+add_sub_impl!(Div, div, SimpleSub;
     (D, U1), (D, U1) -> (D) for D: DimName;
     self: &'a Translation<N, D>, right: Translation<N, D>, Output = Translation<N, D>;
     Translation::from(&self.vector - right.vector); 'a);
 
-add_sub_impl!(Div, div, ClosedSub;
+add_sub_impl!(Div, div, SimpleSub;
     (D, U1), (D, U1) -> (D) for D: DimName;
     self: Translation<N, D>, right: &'b Translation<N, D>, Output = Translation<N, D>;
     Translation::from(self.vector - &right.vector); 'b);
 
-add_sub_impl!(Div, div, ClosedSub;
+add_sub_impl!(Div, div, SimpleSub;
     (D, U1), (D, U1) -> (D) for D: DimName;
     self: Translation<N, D>, right: Translation<N, D>, Output = Translation<N, D>;
     Translation::from(self.vector - right.vector); );
@@ -55,22 +56,22 @@ add_sub_impl!(Div, div, ClosedSub;
 // Translation × Point
 // FIXME: we don't handle properly non-zero origins here. Do we want this to be the intended
 // behavior?
-add_sub_impl!(Mul, mul, ClosedAdd;
+add_sub_impl!(Mul, mul, SimpleAdd;
     (D, U1), (D, U1) -> (D) for D: DimName;
     self: &'a Translation<N, D>, right: &'b Point<N, D>, Output = Point<N, D>;
     right + &self.vector; 'a, 'b);
 
-add_sub_impl!(Mul, mul, ClosedAdd;
+add_sub_impl!(Mul, mul, SimpleAdd;
     (D, U1), (D, U1) -> (D) for D: DimName;
     self: &'a Translation<N, D>, right: Point<N, D>, Output = Point<N, D>;
     right + &self.vector; 'a);
 
-add_sub_impl!(Mul, mul, ClosedAdd;
+add_sub_impl!(Mul, mul, SimpleAdd;
     (D, U1), (D, U1) -> (D) for D: DimName;
     self: Translation<N, D>, right: &'b Point<N, D>, Output = Point<N, D>;
     right + self.vector; 'b);
 
-add_sub_impl!(Mul, mul, ClosedAdd;
+add_sub_impl!(Mul, mul, SimpleAdd;
     (D, U1), (D, U1) -> (D) for D: DimName;
     self: Translation<N, D>, right: Point<N, D>, Output = Point<N, D>;
     right + self.vector; );
diff --git a/tests/core/matrix.rs b/tests/core/matrix.rs
index 9e25db58..07988e95 100644
--- a/tests/core/matrix.rs
+++ b/tests/core/matrix.rs
@@ -1107,3 +1107,98 @@ fn partial_eq_different_types() {
     // assert_ne!(static_mat, typenum_static_mat);
     //assert_ne!(typenum_static_mat, static_mat);
 }
+
+#[test]
+fn add_without_add_assign() {
+    // Ensure adding matrices works without implementing AddAssign
+    #[derive(Clone, Copy, Debug, PartialEq)]
+    struct Value(f32);
+    impl std::ops::Add<&Value> for Value {
+        type Output = Self;
+        fn add(self, rhs: &Self) -> Self {
+            Value(self.0 + rhs.0)
+        }
+    }
+    impl std::ops::Add<Value> for Value {
+        type Output = Self;
+        fn add(self, rhs: Self) -> Self {
+            self.add(&rhs)
+        }
+    }
+    impl std::ops::Sub<&Value> for Value {
+        type Output = Self;
+        fn sub(self, rhs: &Self) -> Self {
+            Value(self.0 - rhs.0)
+        }
+    }
+    impl std::ops::Sub<Value> for Value {
+        type Output = Self;
+        fn sub(self, rhs: Self) -> Self {
+            self.sub(&rhs)
+        }
+    }
+
+    let a = Matrix2x3::new(
+        Value(1.0),
+        Value(2.0),
+        Value(3.0),
+        Value(4.0),
+        Value(5.0),
+        Value(6.0),
+    );
+
+    let b = Matrix2x3::new(
+        Value(10.0),
+        Value(20.0),
+        Value(30.0),
+        Value(40.0),
+        Value(50.0),
+        Value(60.0),
+    );
+    let c = DMatrix::from_row_slice(2, 3, &[
+         Value(10.0),
+         Value(20.0),
+         Value(30.0),
+         Value(40.0),
+         Value(50.0),
+         Value(60.0),
+    ]);
+
+    let expected_add = Matrix2x3::new(
+        Value(11.0),
+        Value(22.0),
+        Value(33.0),
+        Value(44.0),
+        Value(55.0),
+        Value(66.0)
+    );
+
+    let expected_sub = Matrix2x3::new(
+        Value(-9.0),
+        Value(-18.0),
+        Value(-27.0),
+        Value(-36.0),
+        Value(-45.0),
+        Value(-54.0)
+    );
+
+    assert_eq!(expected_add, &a + &b);
+    assert_eq!(expected_add, &a + b);
+    assert_eq!(expected_add, a + &b);
+    assert_eq!(expected_add, a + b);
+
+    // Sum of a static matrix with a dynamic one.
+    assert_eq!(expected_add, &a + &c);
+    assert_eq!(expected_add, a + &c);
+    assert_eq!(expected_add, &c + &a);
+    assert_eq!(expected_add, &c + a);
+
+    assert_eq!(expected_sub, &a - &b);
+    assert_eq!(expected_sub, &a - b);
+    assert_eq!(expected_sub, a - &b);
+    assert_eq!(expected_sub, a - b);
+
+    // Difference of a static matrix with a dynamic one.
+    assert_eq!(expected_sub, &a - &c);
+    assert_eq!(expected_sub, a - &c);
+}
diff --git a/tests/geometry/point.rs b/tests/geometry/point.rs
index 896a09d6..bba5b36c 100644
--- a/tests/geometry/point.rs
+++ b/tests/geometry/point.rs
@@ -93,6 +93,94 @@ fn to_homogeneous() {
     assert_eq!(a.to_homogeneous(), expected);
 }
 
+#[test]
+fn point_ops_without_assign() {
+    // Ensure adding matrices works without implementing AddAssign
+    #[derive(Clone, Copy, Debug, PartialEq)]
+    struct Value(f32);
+    impl std::ops::Add<&Value> for Value {
+        type Output = Self;
+        fn add(self, rhs: &Self) -> Self {
+            Value(self.0 + rhs.0)
+        }
+    }
+    impl std::ops::Add<Value> for Value {
+        type Output = Self;
+        fn add(self, rhs: Self) -> Self {
+            self.add(&rhs)
+        }
+    }
+    impl std::ops::Sub<&Value> for Value {
+        type Output = Self;
+        fn sub(self, rhs: &Self) -> Self {
+            Value(self.0 - rhs.0)
+        }
+    }
+    impl std::ops::Sub<Value> for Value {
+        type Output = Self;
+        fn sub(self, rhs: Self) -> Self {
+            self.sub(&rhs)
+        }
+    }
+    impl std::ops::Mul<Value> for f32 {
+        type Output = Value;
+        fn mul(self, rhs: Value) -> Self::Output {
+            Value(self * rhs.0)
+        }
+    }
+    impl std::ops::Mul<f32> for Value {
+        type Output = Value;
+        fn mul(self, rhs: f32) -> Self::Output {
+            Value(self.0 * rhs)
+        }
+    }
+    impl Zero for Value {
+        fn zero() -> Self {
+            Value(Zero::zero())
+        }
+
+        fn is_zero(&self) -> bool {
+            self.0.is_zero()
+        }
+    }
+
+    let a = Point3::new(
+        Value(1.0),
+        Value(2.0),
+        Value(3.0),
+    );
+    let b = Point3::new(
+        Value(1.0),
+        Value(2.0),
+        Value(3.0),
+    );
+    let c = Vector3::new(
+        Value(1.0),
+        Value(2.0),
+        Value(3.0),
+    );
+
+    assert_eq!(a - b, Vector3::zero());
+    assert_eq!(&a - &b, Vector3::zero());
+    assert_eq!(a - &b, Vector3::zero());
+    assert_eq!(&a - b, Vector3::zero());
+
+    assert_eq!(b - c, Point3::origin());
+    assert_eq!(&b - &c, Point3::origin());
+    assert_eq!(b - &c, Point3::origin());
+    assert_eq!(&b - c, Point3::origin());
+
+    let a2 = Point3::new(
+        Value(2.0),
+        Value(4.0),
+        Value(6.0),
+    );
+    assert_eq!(b + c, a2);
+    assert_eq!(&b + &c, a2);
+    assert_eq!(b + &c, a2);
+    assert_eq!(&b + c, a2);
+}
+
 #[cfg(feature = "arbitrary")]
 quickcheck!(
     fn point_sub(pt1: Point3<f64>, pt2: Point3<f64>) -> bool {