diff --git a/src/aabb.rs b/src/aabb.rs
index 307c02a..c84c40a 100644
--- a/src/aabb.rs
+++ b/src/aabb.rs
@@ -82,8 +82,8 @@ pub trait Aabb<S: BaseNum, V: Vector<S>, P: Point<S, V>>: Sized {
 
     /// Multiply every point in the AABB by a vector, returning a new AABB.
     fn mul_v(&self, v: &V) -> Self {
-        let min : P = Point::from_vec(&self.min().to_vec().mul_v(v));
-        let max : P = Point::from_vec(&self.max().to_vec().mul_v(v));
+        let min = P::from_vec(&self.min().to_vec().mul_v(v));
+        let max = P::from_vec(&self.max().to_vec().mul_v(v));
         Aabb::new(min, max)
     }
 }
diff --git a/src/matrix.rs b/src/matrix.rs
index 5d7f290..433d384 100644
--- a/src/matrix.rs
+++ b/src/matrix.rs
@@ -227,7 +227,7 @@ impl<S: BaseFloat> Matrix4<S> {
     /// Create a transformation matrix that will cause a vector to point at
     /// `dir`, using `up` for orientation.
     pub fn look_at(eye: &Point3<S>, center: &Point3<S>, up: &Vector3<S>) -> Matrix4<S> {
-        let f = center.sub_p(eye).normalize();
+        let f = (center - eye).normalize();
         let s = f.cross(up).normalize();
         let u = s.cross(&f);
 
@@ -249,9 +249,18 @@ impl<S: Copy + Neg<Output = S>> Matrix4<S> {
     }
 }
 
-pub trait Matrix<S: BaseFloat, V: Clone + Vector<S> + 'static>: Array2<V, V, S>
-                                                              + ApproxEq<S>
-                                                              + Sized {
+pub trait Matrix<S: BaseFloat, V: Vector<S> + 'static>: Array2<V, V, S> + ApproxEq<S> + Sized // where
+    // FIXME: blocked by rust-lang/rust#20671
+    //
+    // for<'a, 'b> &'a Self: Add<&'b Self, Output = Self>,
+    // for<'a, 'b> &'a Self: Sub<&'b Self, Output = Self>,
+    // for<'a, 'b> &'a Self: Mul<&'b Self, Output = Self>,
+    // for<'a, 'b> &'a Self: Mul<&'b V, Output = V>,
+    //
+    // for<'a> &'a Self: Mul<S, Output = Self>,
+    // for<'a> &'a Self: Div<S, Output = Self>,
+    // for<'a> &'a Self: Rem<S, Output = Self>,
+{
     /// Create a new diagonal matrix using the supplied value.
     fn from_value(value: S) -> Self;
     /// Create a matrix from a non-uniform scale
@@ -350,117 +359,6 @@ pub trait Matrix<S: BaseFloat, V: Clone + Vector<S> + 'static>: Array2<V, V, S>
     fn is_symmetric(&self) -> bool;
 }
 
-impl<S: BaseFloat> Add for Matrix2<S> {
-    type Output = Matrix2<S>;
-
-    #[inline]
-    fn add(self, other: Matrix2<S>) -> Matrix2<S> { self.add_m(&other) }
-}
-
-impl<S: BaseFloat> Add for Matrix3<S> {
-    type Output = Matrix3<S>;
-
-    #[inline]
-    fn add(self, other: Matrix3<S>) -> Matrix3<S> { self.add_m(&other) }
-}
-
-impl<S: BaseFloat> Add for Matrix4<S> {
-    type Output = Matrix4<S>;
-
-    #[inline]
-    fn add(self, other: Matrix4<S>) -> Matrix4<S> { self.add_m(&other) }
-}
-
-impl<S: BaseFloat> Sub for Matrix2<S> {
-    type Output = Matrix2<S>;
-
-    #[inline]
-    fn sub(self, other: Matrix2<S>) -> Matrix2<S> { self.sub_m(&other) }
-}
-
-impl<S: BaseFloat> Sub for Matrix3<S> {
-    type Output = Matrix3<S>;
-
-    #[inline]
-    fn sub(self, other: Matrix3<S>) -> Matrix3<S> { self.sub_m(&other) }
-}
-
-impl<S: BaseFloat> Sub for Matrix4<S> {
-    type Output = Matrix4<S>;
-
-    #[inline]
-    fn sub(self, other: Matrix4<S>) -> Matrix4<S> { self.sub_m(&other) }
-}
-
-impl<S: Neg<Output = S>> Neg for Matrix2<S> {
-    type Output = Matrix2<S>;
-
-    #[inline]
-    fn neg(self) -> Matrix2<S> {
-        Matrix2::from_cols(self.x.neg(), self.y.neg())
-    }
-}
-
-impl<S: Neg<Output = S>> Neg for Matrix3<S> {
-    type Output = Matrix3<S>;
-
-    #[inline]
-    fn neg(self) -> Matrix3<S> {
-        Matrix3::from_cols(self.x.neg(), self.y.neg(), self.z.neg())
-    }
-}
-
-impl<S: Neg<Output = S>> Neg for Matrix4<S> {
-    type Output = Matrix4<S>;
-
-    #[inline]
-    fn neg(self) -> Matrix4<S> {
-        Matrix4::from_cols(self.x.neg(), self.y.neg(), self.z.neg(), self.w.neg())
-    }
-}
-
-impl<S: BaseFloat> Mul for Matrix2<S> {
-    type Output = Matrix2<S>;
-
-    #[inline]
-    fn mul(self, other: Matrix2<S>) -> Matrix2<S> { self.mul_m(&other) }
-}
-
-impl<S: BaseFloat> Mul<S> for Matrix2<S> {
-    type Output = Matrix2<S>;
-
-    #[inline]
-    fn mul(self, other: S) -> Matrix2<S> { self.mul_s(other) }
-}
-
-impl<S: BaseFloat> Mul for Matrix3<S> {
-    type Output = Matrix3<S>;
-
-    #[inline]
-    fn mul(self, other: Matrix3<S>) -> Matrix3<S> { self.mul_m(&other) }
-}
-
-impl<S: BaseFloat> Mul<S> for Matrix3<S> {
-    type Output = Matrix3<S>;
-
-    #[inline]
-    fn mul(self, other: S) -> Matrix3<S> { self.mul_s(other) }
-}
-
-impl<S: BaseFloat> Mul for Matrix4<S> {
-    type Output = Matrix4<S>;
-
-    #[inline]
-    fn mul(self, other: Matrix4<S>) -> Matrix4<S> { self.mul_m(&other) }
-}
-
-impl<S: BaseFloat> Mul<S> for Matrix4<S> {
-    type Output = Matrix4<S>;
-
-    #[inline]
-    fn mul(self, other: S) -> Matrix4<S> { self.mul_s(other) }
-}
-
 impl<S: Copy + 'static> Array2<Vector2<S>, Vector2<S>, S> for Matrix2<S> {
     #[inline]
     fn row(&self, r: usize) -> Vector2<S> {
@@ -522,46 +420,13 @@ impl<S: BaseFloat> Matrix<S, Vector2<S>> for Matrix2<S> {
                      S::zero(), value.y)
     }
 
-    #[inline]
-    fn mul_s(&self, s: S) -> Matrix2<S> {
-        Matrix2::from_cols(self[0].mul_s(s),
-                           self[1].mul_s(s))
-    }
-
-    #[inline]
-    fn div_s(&self, s: S) -> Matrix2<S> {
-        Matrix2::from_cols(self[0].div_s(s),
-                           self[1].div_s(s))
-    }
-
-    #[inline]
-    fn rem_s(&self, s: S) -> Matrix2<S> {
-        Matrix2::from_cols(self[0].rem_s(s),
-                           self[1].rem_s(s))
-    }
-
-    #[inline]
-    fn add_m(&self, m: &Matrix2<S>) -> Matrix2<S> {
-        Matrix2::from_cols(self[0].add_v(&m[0]),
-                           self[1].add_v(&m[1]))
-    }
-
-    #[inline]
-    fn sub_m(&self, m: &Matrix2<S>) -> Matrix2<S> {
-        Matrix2::from_cols(self[0].sub_v(&m[0]),
-                           self[1].sub_v(&m[1]))
-    }
-
-    #[inline]
-    fn mul_v(&self, v: &Vector2<S>) -> Vector2<S> {
-        Vector2::new(self.row(0).dot(v),
-                     self.row(1).dot(v))
-    }
-
-    fn mul_m(&self, other: &Matrix2<S>) -> Matrix2<S> {
-        Matrix2::new(self.row(0).dot(&other[0]), self.row(1).dot(&other[0]),
-                     self.row(0).dot(&other[1]), self.row(1).dot(&other[1]))
-    }
+    #[inline] fn mul_s(&self, s: S) -> Matrix2<S> { self * s }
+    #[inline] fn div_s(&self, s: S) -> Matrix2<S> { self / s }
+    #[inline] fn rem_s(&self, s: S) -> Matrix2<S> { self % s }
+    #[inline] fn add_m(&self, m: &Matrix2<S>) -> Matrix2<S> { self + m }
+    #[inline] fn sub_m(&self, m: &Matrix2<S>) -> Matrix2<S> { self - m }
+    fn mul_m(&self, other: &Matrix2<S>) -> Matrix2<S> { self * other }
+    #[inline] fn mul_v(&self, v: &Vector2<S>) -> Vector2<S> { self * v }
 
     #[inline]
     fn mul_self_s(&mut self, s: S) {
@@ -654,53 +519,13 @@ impl<S: BaseFloat> Matrix<S, Vector3<S>> for Matrix3<S> {
                      S::zero(), S::zero(), value.z)
     }
 
-    #[inline]
-    fn mul_s(&self, s: S) -> Matrix3<S> {
-        Matrix3::from_cols(self[0].mul_s(s),
-                           self[1].mul_s(s),
-                           self[2].mul_s(s))
-    }
-
-    #[inline]
-    fn div_s(&self, s: S) -> Matrix3<S> {
-        Matrix3::from_cols(self[0].div_s(s),
-                           self[1].div_s(s),
-                           self[2].div_s(s))
-    }
-
-    #[inline]
-    fn rem_s(&self, s: S) -> Matrix3<S> {
-        Matrix3::from_cols(self[0].rem_s(s),
-                           self[1].rem_s(s),
-                           self[2].rem_s(s))
-    }
-
-    #[inline]
-    fn add_m(&self, m: &Matrix3<S>) -> Matrix3<S> {
-        Matrix3::from_cols(self[0].add_v(&m[0]),
-                           self[1].add_v(&m[1]),
-                           self[2].add_v(&m[2]))
-    }
-
-    #[inline]
-    fn sub_m(&self, m: &Matrix3<S>) -> Matrix3<S> {
-        Matrix3::from_cols(self[0].sub_v(&m[0]),
-                           self[1].sub_v(&m[1]),
-                           self[2].sub_v(&m[2]))
-    }
-
-    #[inline]
-    fn mul_v(&self, v: &Vector3<S>) -> Vector3<S> {
-        Vector3::new(self.row(0).dot(v),
-                     self.row(1).dot(v),
-                     self.row(2).dot(v))
-    }
-
-    fn mul_m(&self, other: &Matrix3<S>) -> Matrix3<S> {
-        Matrix3::new(self.row(0).dot(&other[0]),self.row(1).dot(&other[0]),self.row(2).dot(&other[0]),
-                     self.row(0).dot(&other[1]),self.row(1).dot(&other[1]),self.row(2).dot(&other[1]),
-                     self.row(0).dot(&other[2]),self.row(1).dot(&other[2]),self.row(2).dot(&other[2]))
-    }
+    #[inline] fn mul_s(&self, s: S) -> Matrix3<S> { self * s }
+    #[inline] fn div_s(&self, s: S) -> Matrix3<S> { self / s }
+    #[inline] fn rem_s(&self, s: S) -> Matrix3<S> { self % s }
+    #[inline] fn add_m(&self, m: &Matrix3<S>) -> Matrix3<S> { self + m }
+    #[inline] fn sub_m(&self, m: &Matrix3<S>) -> Matrix3<S> { self - m }
+    fn mul_m(&self, other: &Matrix3<S>) -> Matrix3<S> { self * other }
+    #[inline] fn mul_v(&self, v: &Vector3<S>) -> Vector3<S> { self * v}
 
     #[inline]
     fn mul_self_s(&mut self, s: S) {
@@ -766,9 +591,9 @@ impl<S: BaseFloat> Matrix<S, Vector3<S>> for Matrix3<S> {
     fn invert(&self) -> Option<Matrix3<S>> {
         let det = self.determinant();
         if det.approx_eq(&S::zero()) { None } else {
-            Some(Matrix3::from_cols(self[1].cross(&self[2]).div_s(det),
-                                    self[2].cross(&self[0]).div_s(det),
-                                    self[0].cross(&self[1]).div_s(det)).transpose())
+            Some(Matrix3::from_cols(&self[1].cross(&self[2]) / det,
+                                    &self[2].cross(&self[0]) / det,
+                                    &self[0].cross(&self[1]) / det).transpose())
         }
     }
 
@@ -795,18 +620,6 @@ impl<S: BaseFloat> Matrix<S, Vector3<S>> for Matrix3<S> {
     }
 }
 
-// Using self.row(0).dot(other[0]) like the other matrix multiplies
-// causes the LLVM to miss identical loads and multiplies. This optimization
-// causes the code to be auto vectorized properly increasing the performance
-// around ~4 times.
-macro_rules! dot_matrix4(
-    ($A:expr, $B:expr, $I:expr, $J:expr) => (
-        ($A[0][$I]) * ($B[$J][0]) +
-        ($A[1][$I]) * ($B[$J][1]) +
-        ($A[2][$I]) * ($B[$J][2]) +
-        ($A[3][$I]) * ($B[$J][3])
-));
-
 impl<S: BaseFloat> Matrix<S, Vector4<S>> for Matrix4<S> {
     #[inline]
     fn from_value(value: S) -> Matrix4<S> {
@@ -824,60 +637,13 @@ impl<S: BaseFloat> Matrix<S, Vector4<S>> for Matrix4<S> {
                      S::zero(), S::zero(), S::zero(), value.w)
     }
 
-    #[inline]
-    fn mul_s(&self, s: S) -> Matrix4<S> {
-        Matrix4::from_cols(self[0].mul_s(s),
-                           self[1].mul_s(s),
-                           self[2].mul_s(s),
-                           self[3].mul_s(s))
-    }
-
-    #[inline]
-    fn div_s(&self, s: S) -> Matrix4<S> {
-        Matrix4::from_cols(self[0].div_s(s),
-                           self[1].div_s(s),
-                           self[2].div_s(s),
-                           self[3].div_s(s))
-    }
-
-    #[inline]
-    fn rem_s(&self, s: S) -> Matrix4<S> {
-        Matrix4::from_cols(self[0].rem_s(s),
-                           self[1].rem_s(s),
-                           self[2].rem_s(s),
-                           self[3].rem_s(s))
-    }
-
-    #[inline]
-    fn add_m(&self, m: &Matrix4<S>) -> Matrix4<S> {
-        Matrix4::from_cols(self[0].add_v(&m[0]),
-                           self[1].add_v(&m[1]),
-                           self[2].add_v(&m[2]),
-                           self[3].add_v(&m[3]))
-    }
-
-    #[inline]
-    fn sub_m(&self, m: &Matrix4<S>) -> Matrix4<S> {
-        Matrix4::from_cols(self[0].sub_v(&m[0]),
-                           self[1].sub_v(&m[1]),
-                           self[2].sub_v(&m[2]),
-                           self[3].sub_v(&m[3]))
-    }
-
-    #[inline]
-    fn mul_v(&self, v: &Vector4<S>) -> Vector4<S> {
-        Vector4::new(self.row(0).dot(v),
-                     self.row(1).dot(v),
-                     self.row(2).dot(v),
-                     self.row(3).dot(v))
-    }
-
-    fn mul_m(&self, other: &Matrix4<S>) -> Matrix4<S> {
-        Matrix4::new(dot_matrix4!(self, other, 0, 0), dot_matrix4!(self, other, 1, 0), dot_matrix4!(self, other, 2, 0), dot_matrix4!(self, other, 3, 0),
-                     dot_matrix4!(self, other, 0, 1), dot_matrix4!(self, other, 1, 1), dot_matrix4!(self, other, 2, 1), dot_matrix4!(self, other, 3, 1),
-                     dot_matrix4!(self, other, 0, 2), dot_matrix4!(self, other, 1, 2), dot_matrix4!(self, other, 2, 2), dot_matrix4!(self, other, 3, 2),
-                     dot_matrix4!(self, other, 0, 3), dot_matrix4!(self, other, 1, 3), dot_matrix4!(self, other, 2, 3), dot_matrix4!(self, other, 3, 3))
-    }
+    #[inline] fn mul_s(&self, s: S) -> Matrix4<S> { self * s }
+    #[inline] fn div_s(&self, s: S) -> Matrix4<S> { self / s }
+    #[inline] fn rem_s(&self, s: S) -> Matrix4<S> { self % s }
+    #[inline] fn add_m(&self, m: &Matrix4<S>) -> Matrix4<S> { self + m }
+    #[inline] fn sub_m(&self, m: &Matrix4<S>) -> Matrix4<S> { self - m }
+    fn mul_m(&self, other: &Matrix4<S>) -> Matrix4<S> { self * other }
+    #[inline] fn mul_v(&self, v: &Vector4<S>) -> Vector4<S> { self * v }
 
     #[inline]
     fn mul_self_s(&mut self, s: S) {
@@ -965,24 +731,16 @@ impl<S: BaseFloat> Matrix<S, Vector4<S>> for Matrix4<S> {
 
     fn invert(&self) -> Option<Matrix4<S>> {
         let det = self.determinant();
-        if !det.approx_eq(&S::zero()) {
+        if det.approx_eq(&S::zero()) { None } else {
             let inv_det = S::one() / det;
             let t = self.transpose();
             let cf = |i, j| {
                 let mat = match i {
-                    0 => Matrix3::from_cols(t.y.truncate_n(j),
-                                            t.z.truncate_n(j),
-                                            t.w.truncate_n(j)),
-                    1 => Matrix3::from_cols(t.x.truncate_n(j),
-                                            t.z.truncate_n(j),
-                                            t.w.truncate_n(j)),
-                    2 => Matrix3::from_cols(t.x.truncate_n(j),
-                                            t.y.truncate_n(j),
-                                            t.w.truncate_n(j)),
-                    3 => Matrix3::from_cols(t.x.truncate_n(j),
-                                            t.y.truncate_n(j),
-                                            t.z.truncate_n(j)),
-                    _ => panic!("out of range")
+                    0 => Matrix3::from_cols(t.y.truncate_n(j), t.z.truncate_n(j), t.w.truncate_n(j)),
+                    1 => Matrix3::from_cols(t.x.truncate_n(j), t.z.truncate_n(j), t.w.truncate_n(j)),
+                    2 => Matrix3::from_cols(t.x.truncate_n(j), t.y.truncate_n(j), t.w.truncate_n(j)),
+                    3 => Matrix3::from_cols(t.x.truncate_n(j), t.y.truncate_n(j), t.z.truncate_n(j)),
+                    _ => panic!("out of range"),
                 };
                 let sign = if (i + j) & 1 == 1 { -S::one() } else { S::one() };
                 mat.determinant() * sign * inv_det
@@ -992,9 +750,6 @@ impl<S: BaseFloat> Matrix<S, Vector4<S>> for Matrix4<S> {
                               cf(1, 0), cf(1, 1), cf(1, 2), cf(1, 3),
                               cf(2, 0), cf(2, 1), cf(2, 2), cf(2, 3),
                               cf(3, 0), cf(3, 1), cf(3, 2), cf(3, 3)))
-
-        } else {
-            None
         }
     }
 
@@ -1062,6 +817,150 @@ impl<S: BaseFloat> ApproxEq<S> for Matrix4<S> {
     }
 }
 
+impl<S: Neg<Output = S>> Neg for Matrix2<S> {
+    type Output = Matrix2<S>;
+
+    #[inline]
+    fn neg(self) -> Matrix2<S> {
+        Matrix2::from_cols(-self.x, -self.y)
+    }
+}
+
+impl<S: Neg<Output = S>> Neg for Matrix3<S> {
+    type Output = Matrix3<S>;
+
+    #[inline]
+    fn neg(self) -> Matrix3<S> {
+        Matrix3::from_cols(-self.x, -self.y, -self.z)
+    }
+}
+
+impl<S: Neg<Output = S>> Neg for Matrix4<S> {
+    type Output = Matrix4<S>;
+
+    #[inline]
+    fn neg(self) -> Matrix4<S> {
+        Matrix4::from_cols(-self.x, -self.y, -self.z, -self.w)
+    }
+}
+
+macro_rules! impl_scalar_binary_operator {
+    ($Binop:ident :: $binop:ident, $MatrixN:ident { $($field:ident),+ }) => {
+        impl<'a, S: BaseNum> $Binop<S> for &'a $MatrixN<S> {
+            type Output = $MatrixN<S>;
+
+            #[inline]
+            fn $binop(self, s: S) -> $MatrixN<S> {
+                $MatrixN { $($field: self.$field.$binop(s)),+ }
+            }
+        }
+    }
+}
+
+impl_scalar_binary_operator!(Mul::mul, Matrix2 { x, y });
+impl_scalar_binary_operator!(Mul::mul, Matrix3 { x, y, z });
+impl_scalar_binary_operator!(Mul::mul, Matrix4 { x, y, z, w });
+impl_scalar_binary_operator!(Div::div, Matrix2 { x, y });
+impl_scalar_binary_operator!(Div::div, Matrix3 { x, y, z });
+impl_scalar_binary_operator!(Div::div, Matrix4 { x, y, z, w });
+impl_scalar_binary_operator!(Rem::rem, Matrix2 { x, y });
+impl_scalar_binary_operator!(Rem::rem, Matrix3 { x, y, z });
+impl_scalar_binary_operator!(Rem::rem, Matrix4 { x, y, z, w });
+
+macro_rules! impl_binary_operator {
+    ($Binop:ident :: $binop:ident, $MatrixN:ident { $($field:ident),+ }) => {
+        impl<'a, 'b, S: BaseNum> $Binop<&'a $MatrixN<S>> for &'b $MatrixN<S> {
+            type Output = $MatrixN<S>;
+
+            #[inline]
+            fn $binop(self, other: &'a $MatrixN<S>) -> $MatrixN<S> {
+                $MatrixN { $($field: self.$field.$binop(&other.$field)),+ }
+            }
+        }
+    }
+}
+
+impl_binary_operator!(Add::add, Matrix2 { x, y });
+impl_binary_operator!(Add::add, Matrix3 { x, y, z });
+impl_binary_operator!(Add::add, Matrix4 { x, y, z, w });
+impl_binary_operator!(Sub::sub, Matrix2 { x, y });
+impl_binary_operator!(Sub::sub, Matrix3 { x, y, z });
+impl_binary_operator!(Sub::sub, Matrix4 { x, y, z, w });
+
+impl<'a, 'b, S: BaseNum> Mul<&'a Vector2<S>> for &'b Matrix2<S> {
+    type Output = Vector2<S>;
+
+    fn mul(self, v: &'a Vector2<S>) -> Vector2<S> {
+        Vector2::new(self.row(0).dot(v),
+                     self.row(1).dot(v))
+    }
+}
+
+impl<'a, 'b, S: BaseNum> Mul<&'a Vector3<S>> for &'b Matrix3<S> {
+    type Output = Vector3<S>;
+
+    fn mul(self, v: &'a Vector3<S>) -> Vector3<S> {
+        Vector3::new(self.row(0).dot(v),
+                     self.row(1).dot(v),
+                     self.row(2).dot(v))
+    }
+}
+
+impl<'a, 'b, S: BaseNum> Mul<&'a Vector4<S>> for &'b Matrix4<S> {
+    type Output = Vector4<S>;
+
+    fn mul(self, v: &'a Vector4<S>) -> Vector4<S> {
+        Vector4::new(self.row(0).dot(v),
+                     self.row(1).dot(v),
+                     self.row(2).dot(v),
+                     self.row(3).dot(v))
+    }
+}
+
+impl<'a, 'b, S: BaseNum> Mul<&'a Matrix2<S>> for &'b Matrix2<S> {
+    type Output = Matrix2<S>;
+
+    fn mul(self, other: &'a Matrix2<S>) -> Matrix2<S> {
+        Matrix2::new(self.row(0).dot(&other[0]), self.row(1).dot(&other[0]),
+                     self.row(0).dot(&other[1]), self.row(1).dot(&other[1]))
+    }
+}
+
+impl<'a, 'b, S: BaseNum> Mul<&'a Matrix3<S>> for &'b Matrix3<S> {
+    type Output = Matrix3<S>;
+
+    fn mul(self, other: &'a Matrix3<S>) -> Matrix3<S> {
+        Matrix3::new(self.row(0).dot(&other[0]),self.row(1).dot(&other[0]),self.row(2).dot(&other[0]),
+                     self.row(0).dot(&other[1]),self.row(1).dot(&other[1]),self.row(2).dot(&other[1]),
+                     self.row(0).dot(&other[2]),self.row(1).dot(&other[2]),self.row(2).dot(&other[2]))
+    }
+}
+
+impl<'a, 'b, S: BaseNum> Mul<&'a Matrix4<S>> for &'b Matrix4<S> {
+    type Output = Matrix4<S>;
+
+    fn mul(self, other: &'a Matrix4<S>) -> Matrix4<S> {
+        // Using self.row(0).dot(other[0]) like the other matrix multiplies
+        // causes the LLVM to miss identical loads and multiplies. This optimization
+        // causes the code to be auto vectorized properly increasing the performance
+        // around ~4 times.
+        macro_rules! dot_matrix4 {
+            ($A:expr, $B:expr, $I:expr, $J:expr) => {
+                ($A[0][$I]) * ($B[$J][0]) +
+                ($A[1][$I]) * ($B[$J][1]) +
+                ($A[2][$I]) * ($B[$J][2]) +
+                ($A[3][$I]) * ($B[$J][3])
+            };
+        };
+
+        Matrix4::new(dot_matrix4!(self, other, 0, 0), dot_matrix4!(self, other, 1, 0), dot_matrix4!(self, other, 2, 0), dot_matrix4!(self, other, 3, 0),
+                     dot_matrix4!(self, other, 0, 1), dot_matrix4!(self, other, 1, 1), dot_matrix4!(self, other, 2, 1), dot_matrix4!(self, other, 3, 1),
+                     dot_matrix4!(self, other, 0, 2), dot_matrix4!(self, other, 1, 2), dot_matrix4!(self, other, 2, 2), dot_matrix4!(self, other, 3, 2),
+                     dot_matrix4!(self, other, 0, 3), dot_matrix4!(self, other, 1, 3), dot_matrix4!(self, other, 2, 3), dot_matrix4!(self, other, 3, 3))
+
+    }
+}
+
 macro_rules! index_operators {
     ($MatrixN:ident<$S:ident>, $n:expr, $Output:ty, $I:ty) => {
         impl<$S> Index<$I> for $MatrixN<$S> {
diff --git a/src/point.rs b/src/point.rs
index e9ebf96..94c8ce7 100644
--- a/src/point.rs
+++ b/src/point.rs
@@ -68,7 +68,16 @@ impl<S: BaseNum> Point3<S> {
 }
 
 /// Specifies the numeric operations for point types.
-pub trait Point<S: BaseNum, V: Vector<S>>: Array1<S> + Clone {
+pub trait Point<S: BaseNum, V: Vector<S>>: Array1<S> + Clone // where
+    // FIXME: blocked by rust-lang/rust#20671
+    //
+    // for<'a, 'b> &'a Self: Add<&'b V, Output = Self>,
+    // for<'a, 'b> &'a Self: Sub<&'b Self, Output = V>,
+    //
+    // for<'a> &'a Self: Mul<S, Output = Self>,
+    // for<'a> &'a Self: Div<S, Output = Self>,
+    // for<'a> &'a Self: Rem<S, Output = Self>,
+{
     /// Create a point at the origin.
     fn origin() -> Self;
 
@@ -128,39 +137,14 @@ impl<S: BaseNum> Point<S, Vector2<S>> for Point2<S> {
 
     #[inline]
     fn to_vec(&self) -> Vector2<S> {
-        Vector2::new(self.x,
-                     self.y)
+        Vector2::new(self.x, self.y)
     }
 
-    #[inline]
-    fn mul_s(&self, s: S) -> Point2<S> {
-        Point2::new(self.x * s,
-                    self.y * s)
-    }
-
-    #[inline]
-    fn div_s(&self, s: S) -> Point2<S> {
-        Point2::new(self.x / s,
-                    self.y / s)
-    }
-
-    #[inline]
-    fn rem_s(&self, s: S) -> Point2<S> {
-        Point2::new(self.x % s,
-                    self.y % s)
-    }
-
-    #[inline]
-    fn add_v(&self, v: &Vector2<S>) -> Point2<S> {
-        Point2::new(self.x + v.x,
-                    self.y + v.y)
-    }
-
-    #[inline]
-    fn sub_p(&self, p: &Point2<S>) -> Vector2<S> {
-        Vector2::new(self.x - p.x,
-                     self.y - p.y)
-    }
+    #[inline] fn mul_s(&self, s: S) -> Point2<S> { self * s }
+    #[inline] fn div_s(&self, s: S) -> Point2<S> { self / s }
+    #[inline] fn rem_s(&self, s: S) -> Point2<S> { self % s }
+    #[inline] fn add_v(&self, v: &Vector2<S>) -> Point2<S> { self + v }
+    #[inline] fn sub_p(&self, p: &Point2<S>) -> Vector2<S> { self - p }
 
     #[inline]
     fn mul_self_s(&mut self, s: S) {
@@ -194,14 +178,12 @@ impl<S: BaseNum> Point<S, Vector2<S>> for Point2<S> {
 
     #[inline]
     fn min(&self, p: &Point2<S>) -> Point2<S> {
-        Point2::new(self.x.partial_min(p.x),
-                    self.y.partial_min(p.y))
+        Point2::new(self.x.partial_min(p.x), self.y.partial_min(p.y))
     }
 
     #[inline]
     fn max(&self, p: &Point2<S>) -> Point2<S> {
-        Point2::new(self.x.partial_max(p.x),
-                    self.y.partial_max(p.y))
+        Point2::new(self.x.partial_max(p.x), self.y.partial_max(p.y))
     }
 }
 
@@ -228,45 +210,14 @@ impl<S: BaseNum> Point<S, Vector3<S>> for Point3<S> {
 
     #[inline]
     fn to_vec(&self) -> Vector3<S> {
-        Vector3::new(self.x,
-                     self.y,
-                     self.z)
+        Vector3::new(self.x, self.y, self.z)
     }
 
-    #[inline]
-    fn mul_s(&self, s: S) -> Point3<S> {
-        Point3::new(self.x * s,
-                    self.y * s,
-                    self.z * s)
-    }
-
-    #[inline]
-    fn div_s(&self, s: S) -> Point3<S> {
-        Point3::new(self.x / s,
-                    self.y / s,
-                    self.z / s)
-    }
-
-    #[inline]
-    fn rem_s(&self, s: S) -> Point3<S> {
-        Point3::new(self.x % s,
-                    self.y % s,
-                    self.z % s)
-    }
-
-    #[inline]
-    fn add_v(&self, v: &Vector3<S>) -> Point3<S> {
-        Point3::new(self.x + v.x,
-                    self.y + v.y,
-                    self.z + v.z)
-    }
-
-    #[inline]
-    fn sub_p(&self, p: &Point3<S>) -> Vector3<S> {
-        Vector3::new(self.x - p.x,
-                     self.y - p.y,
-                     self.z - p.z)
-    }
+    #[inline] fn mul_s(&self, s: S) -> Point3<S> { self * s }
+    #[inline] fn div_s(&self, s: S) -> Point3<S> { self / s }
+    #[inline] fn rem_s(&self, s: S) -> Point3<S> { self % s }
+    #[inline] fn add_v(&self, v: &Vector3<S>) -> Point3<S> { self + v }
+    #[inline] fn sub_p(&self, p: &Point3<S>) -> Vector3<S> { self - p }
 
     #[inline]
     fn mul_self_s(&mut self, s: S) {
@@ -305,16 +256,12 @@ impl<S: BaseNum> Point<S, Vector3<S>> for Point3<S> {
 
     #[inline]
     fn min(&self, p: &Point3<S>) -> Point3<S> {
-        Point3::new(self.x.partial_min(p.x),
-                    self.y.partial_min(p.y),
-                    self.z.partial_min(p.z))
+        Point3::new(self.x.partial_min(p.x), self.y.partial_min(p.y), self.z.partial_min(p.z))
     }
 
     #[inline]
     fn max(&self, p: &Point3<S>) -> Point3<S> {
-        Point3::new(self.x.partial_max(p.x),
-                    self.y.partial_max(p.y),
-                    self.z.partial_max(p.z))
+        Point3::new(self.x.partial_max(p.x), self.y.partial_max(p.y), self.z.partial_max(p.z))
     }
 }
 
@@ -327,6 +274,59 @@ impl<S: BaseFloat> ApproxEq<S> for Point3<S> {
     }
 }
 
+
+macro_rules! impl_operators {
+    ($PointN:ident { $($field:ident),+ }, $VectorN:ident) => {
+        impl<'a, S: BaseNum> Mul<S> for &'a $PointN<S> {
+            type Output = $PointN<S>;
+
+            #[inline]
+            fn mul(self, s: S) -> $PointN<S> {
+                $PointN::new($(self.$field * s),+)
+            }
+        }
+
+        impl<'a, S: BaseNum> Div<S> for &'a $PointN<S> {
+            type Output = $PointN<S>;
+
+            #[inline]
+            fn div(self, s: S) -> $PointN<S> {
+                $PointN::new($(self.$field / s),+)
+            }
+        }
+
+        impl<'a, S: BaseNum> Rem<S> for &'a $PointN<S> {
+            type Output = $PointN<S>;
+
+            #[inline]
+            fn rem(self, s: S) -> $PointN<S> {
+                $PointN::new($(self.$field % s),+)
+            }
+        }
+
+        impl<'a, 'b, S: BaseNum> Add<&'a $VectorN<S>> for &'b $PointN<S> {
+            type Output = $PointN<S>;
+
+            #[inline]
+            fn add(self, v: &'a $VectorN<S>) -> $PointN<S> {
+                $PointN::new($(self.$field + v.$field),+)
+            }
+        }
+
+        impl<'a, 'b, S: BaseNum> Sub<&'a $PointN<S>> for &'b $PointN<S> {
+            type Output = $VectorN<S>;
+
+            #[inline]
+            fn sub(self, p: &'a $PointN<S>) -> $VectorN<S> {
+                $VectorN::new($(self.$field - p.$field),+)
+            }
+        }
+    }
+}
+
+impl_operators!(Point2 { x, y }, Vector2);
+impl_operators!(Point3 { x, y, z }, Vector3);
+
 macro_rules! fixed_array_conversions {
     ($PointN:ident <$S:ident> { $($field:ident : $index:expr),+ }, $n:expr) => {
         impl<$S> Into<[$S; $n]> for $PointN<$S> {
diff --git a/src/quaternion.rs b/src/quaternion.rs
index dceec36..ed8fd8a 100644
--- a/src/quaternion.rs
+++ b/src/quaternion.rs
@@ -68,79 +68,6 @@ impl<S: BaseFloat> Quaternion<S> {
         Quaternion::from_sv(S::one(), Vector3::zero())
     }
 
-    /// The result of multiplying the quaternion a scalar
-    #[inline]
-    pub fn mul_s(&self, value: S) -> Quaternion<S> {
-        Quaternion::from_sv(self.s * value, self.v.mul_s(value))
-    }
-
-    /// The result of dividing the quaternion a scalar
-    #[inline]
-    pub fn div_s(&self, value: S) -> Quaternion<S> {
-        Quaternion::from_sv(self.s / value, self.v.div_s(value))
-    }
-
-    /// The result of multiplying the quaternion by a vector
-    #[inline]
-    pub fn mul_v(&self, vec: &Vector3<S>) -> Vector3<S>  {
-        let tmp = self.v.cross(vec).add_v(&vec.mul_s(self.s.clone()));
-        self.v.cross(&tmp).mul_s(cast(2i8).unwrap()).add_v(vec)
-    }
-
-    /// The sum of this quaternion and `other`
-    #[inline]
-    pub fn add_q(&self, other: &Quaternion<S>) -> Quaternion<S> {
-        Quaternion::from_sv(self.s + other.s, self.v + other.v)
-    }
-
-    /// The difference between this quaternion and `other`
-    #[inline]
-    pub fn sub_q(&self, other: &Quaternion<S>) -> Quaternion<S> {
-        Quaternion::from_sv(self.s - other.s, self.v - other.v)
-    }
-
-    /// The result of multipliplying the quaternion by `other`
-    pub fn mul_q(&self, other: &Quaternion<S>) -> Quaternion<S> {
-        Quaternion::new(self.s * other.s - self.v.x * other.v.x - self.v.y * other.v.y - self.v.z * other.v.z,
-                        self.s * other.v.x + self.v.x * other.s + self.v.y * other.v.z - self.v.z * other.v.y,
-                        self.s * other.v.y + self.v.y * other.s + self.v.z * other.v.x - self.v.x * other.v.z,
-                        self.s * other.v.z + self.v.z * other.s + self.v.x * other.v.y - self.v.y * other.v.x)
-    }
-
-    /// Multiply this quaternion by a scalar, in-place.
-    #[inline]
-    pub fn mul_self_s(&mut self, s: S) {
-        self.s = self.s * s;
-        self.v.mul_self_s(s);
-    }
-
-    /// Divide this quaternion by a scalar, in-place.
-    #[inline]
-    pub fn div_self_s(&mut self, s: S) {
-        self.s = self.s / s;
-        self.v.div_self_s(s);
-    }
-
-    /// Add this quaternion by another, in-place.
-    #[inline]
-    pub fn add_self_q(&mut self, q: &Quaternion<S>) {
-        self.s = self.s + q.s;
-        self.v.add_self_v(&q.v);
-    }
-
-    /// Subtract another quaternion from this one, in-place.
-    #[inline]
-    pub fn sub_self_q(&mut self, q: &Quaternion<S>) {
-        self.s = self.s - q.s;
-        self.v.sub_self_v(&q.v);
-    }
-
-    /// Multiply this quaternion by another, in-place.
-    #[inline]
-    pub fn mul_self_q(&mut self, q: &Quaternion<S>) {
-        *self = self.mul_q(q);
-    }
-
     /// The dot product of the quaternion and `q`.
     #[inline]
     pub fn dot(&self, q: &Quaternion<S>) -> S {
@@ -176,12 +103,70 @@ impl<S: BaseFloat> Quaternion<S> {
     /// Normalize this quaternion, returning the new quaternion.
     #[inline]
     pub fn normalize(&self) -> Quaternion<S> {
-        self.mul_s(S::one() / self.magnitude())
+        self * (S::one() / self.magnitude())
     }
 
     /// Do a normalized linear interpolation with `other`, by `amount`.
     pub fn nlerp(&self, other: &Quaternion<S>, amount: S) -> Quaternion<S> {
-        self.mul_s(S::one() - amount).add_q(&other.mul_s(amount)).normalize()
+        (&(self * (S::one() - amount)) + &(other * amount)).normalize()
+    }
+}
+
+impl<'a, S: BaseFloat> Mul<S> for &'a Quaternion<S> {
+    type Output = Quaternion<S>;
+
+    #[inline]
+    fn mul(self, value: S) -> Quaternion<S> {
+        Quaternion::from_sv(self.s * value, &self.v * value)
+    }
+}
+
+impl<'a, S: BaseFloat> Div<S> for &'a Quaternion<S> {
+    type Output = Quaternion<S>;
+
+    #[inline]
+    fn div(self, value: S) -> Quaternion<S> {
+        Quaternion::from_sv(self.s / value, &self.v / value)
+    }
+}
+
+impl<'a, 'b, S: BaseFloat> Mul<&'b Vector3<S>> for &'a Quaternion<S> {
+    type Output = Vector3<S>;
+
+    #[inline]
+    fn mul(self, vec: &'b Vector3<S>) -> Vector3<S>  {
+        let tmp = self.v.cross(vec).add_v(&vec.mul_s(self.s.clone()));
+        self.v.cross(&tmp).mul_s(cast(2i8).unwrap()).add_v(vec)
+    }
+}
+
+impl<'a, 'b, S: BaseFloat> Add<&'b Quaternion<S>> for &'a Quaternion<S> {
+    type Output = Quaternion<S>;
+
+    #[inline]
+    fn add(self, other: &'b Quaternion<S>) -> Quaternion<S> {
+        Quaternion::from_sv(self.s + other.s, &self.v + &other.v)
+    }
+}
+
+impl<'a, 'b, S: BaseFloat> Sub<&'b Quaternion<S>> for &'a Quaternion<S> {
+    type Output = Quaternion<S>;
+
+    #[inline]
+    fn sub(self, other: &'b Quaternion<S>) -> Quaternion<S> {
+        Quaternion::from_sv(self.s - other.s, &self.v - &other.v)
+    }
+}
+
+impl<'a, 'b, S: BaseFloat> Mul<&'b Quaternion<S>> for &'a Quaternion<S> {
+    type Output = Quaternion<S>;
+
+    #[inline]
+    fn mul(self, other: &'b Quaternion<S>) -> Quaternion<S>  {
+        Quaternion::new(self.s * other.s - self.v.x * other.v.x - self.v.y * other.v.y - self.v.z * other.v.z,
+                        self.s * other.v.x + self.v.x * other.s + self.v.y * other.v.z - self.v.z * other.v.y,
+                        self.s * other.v.y + self.v.y * other.s + self.v.z * other.v.x - self.v.x * other.v.z,
+                        self.s * other.v.z + self.v.z * other.s + self.v.x * other.v.y - self.v.y * other.v.x)
     }
 }
 
@@ -233,9 +218,7 @@ impl<S: BaseFloat> Quaternion<S> {
             let scale1 = sin(theta.mul_s(S::one() - amount));
             let scale2 = sin(theta.mul_s(amount));
 
-            self.mul_s(scale1)
-                .add_q(&other.mul_s(scale2))
-                .mul_s(sin(theta).recip())
+            &(&(self * scale1) + &(other * scale2)) * sin(theta).recip()
         }
     }
 
@@ -251,10 +234,10 @@ impl<S: BaseFloat> Quaternion<S> {
         let one: S = cast(1f64).unwrap();
 
         let (qw, qx, qy, qz) = (self.s, self.v.x, self.v.y, self.v.z);
-        let (sqw, sqx, sqy, sqz) = (qw*qw, qx*qx, qy*qy, qz*qz);
+        let (sqw, sqx, sqy, sqz) = (qw * qw, qx * qx, qy * qy, qz * qz);
 
         let unit = sqx + sqy + sqz + sqw;
-        let test = qx*qy + qz*qw;
+        let test = qx * qy + qz * qw;
 
         if test > sig * unit {
             (
@@ -271,9 +254,9 @@ impl<S: BaseFloat> Quaternion<S> {
             )
         } else {
             (
-                rad((two * (qy*qw - qx*qz)).atan2(one - two*(sqy + sqz))),
-                rad((two * (qx*qy + qz*qw)).asin()),
-                rad((two * (qx*qw - qy*qz)).atan2(one - two*(sqx + sqz))),
+                rad((two * (qy * qw - qx * qz)).atan2(one - two * (sqy + sqz))),
+                rad((two * (qx * qy + qz * qw)).asin()),
+                rad((two * (qx * qw - qy * qz)).atan2(one - two * (sqx + sqz))),
             )
         }
     }
@@ -373,16 +356,16 @@ impl<S: BaseFloat + 'static> Rotation<S, Vector3<S>, Point3<S>> for Quaternion<S
     }
 
     #[inline]
-    fn rotate_vector(&self, vec: &Vector3<S>) -> Vector3<S> { self.mul_v(vec) }
+    fn rotate_vector(&self, vec: &Vector3<S>) -> Vector3<S> { self * vec }
 
     #[inline]
-    fn concat(&self, other: &Quaternion<S>) -> Quaternion<S> { self.mul_q(other) }
+    fn concat(&self, other: &Quaternion<S>) -> Quaternion<S> { self * other }
 
     #[inline]
-    fn concat_self(&mut self, other: &Quaternion<S>) { self.mul_self_q(other); }
+    fn concat_self(&mut self, other: &Quaternion<S>) { *self = &*self * other; }
 
     #[inline]
-    fn invert(&self) -> Quaternion<S> { self.conjugate().div_s(self.magnitude2()) }
+    fn invert(&self) -> Quaternion<S> { &self.conjugate() / self.magnitude2() }
 
     #[inline]
     fn invert_self(&mut self) { *self = self.invert() }
diff --git a/src/transform.rs b/src/transform.rs
index 8ef9ee8..cb8d930 100644
--- a/src/transform.rs
+++ b/src/transform.rs
@@ -262,7 +262,7 @@ impl<
     R: Rotation<S, V, P> + Clone,
 > ToComponents<S, V, P, R> for Decomposed<S, V, R> {
     fn decompose(&self) -> (V, R, V) {
-        (V::one().mul_s(self.scale), self.rot.clone(), self.disp.clone())
+        (V::identity().mul_s(self.scale), self.rot.clone(), self.disp.clone())
     }
 }
 
diff --git a/src/vector.rs b/src/vector.rs
index 8a10b8d..7ec191d 100644
--- a/src/vector.rs
+++ b/src/vector.rs
@@ -22,27 +22,26 @@
 //! vector are also provided:
 //!
 //! ```rust
-//! use cgmath::{Vector, Vector2, Vector3, Vector4, zero, vec2, vec3};
+//! use cgmath::{Vector, Vector2, Vector3, Vector4, vec2, vec3};
 //!
 //! assert_eq!(Vector2::new(1.0f64, 0.0f64), Vector2::unit_x());
-//! assert_eq!(vec3(0.0f64, 0.0f64, 0.0f64), zero());
+//! assert_eq!(vec3(0.0f64, 0.0f64, 0.0f64), Vector3::zero());
 //! assert_eq!(Vector2::from_value(1.0f64), vec2(1.0, 1.0));
 //! ```
 //!
 //! Vectors can be manipulated with typical mathematical operations (addition,
 //! subtraction, element-wise multiplication, element-wise division, negation)
-//! using the built-in operators. The additive and multiplicative neutral
-//! elements (zero and one) are also provided by this library
+//! using the built-in operators.
 //!
 //! ```rust
-//! use cgmath::{Vector2, Vector3, Vector4, one, zero};
+//! use cgmath::{Vector, Vector2, Vector3, Vector4};
 //!
 //! let a: Vector2<f64> = Vector2::new(3.0, 4.0);
 //! let b: Vector2<f64> = Vector2::new(-3.0, -4.0);
 //!
-//! assert_eq!(a + b, zero());
-//! assert_eq!(-(a * b), Vector2::new(9.0f64, 16.0f64));
-//! assert_eq!(a / one(), a);
+//! assert_eq!(&a + &b, Vector2::zero());
+//! assert_eq!(-(&a * &b), Vector2::new(9.0f64, 16.0f64));
+//! assert_eq!(&a / &Vector2::identity(), a);
 //!
 //! // As with Rust's `int` and `f32` types, Vectors of different types cannot
 //! // be added and so on with impunity. The following will fail to compile:
@@ -50,10 +49,10 @@
 //!
 //! // Instead, we need to convert the Vector2 to a Vector3 by "extending" it
 //! // with the value for the last coordinate:
-//! let c: Vector3<f64> = a.extend(0.0) + Vector3::new(1.0, 0.0, 2.0);
+//! let c: Vector3<f64> = &a.extend(0.0) + &Vector3::new(1.0, 0.0, 2.0);
 //!
 //! // Similarly, we can "truncate" a Vector4 down to a Vector3:
-//! let d: Vector3<f64> = c + Vector4::unit_x().truncate();
+//! let d: Vector3<f64> = &c + &Vector4::unit_x().truncate();
 //!
 //! assert_eq!(d, Vector3::new(5.0f64, 4.0f64, 2.0f64));
 //! ```
@@ -64,12 +63,12 @@
 //! and [cross products](http://en.wikipedia.org/wiki/Cross_product).
 //!
 //! ```rust
-//! use cgmath::{Vector, Vector2, Vector3, Vector4, dot, zero};
+//! use cgmath::{Vector, Vector2, Vector3, Vector4, dot};
 //!
 //! // All vectors implement the dot product as a method:
 //! let a: Vector2<f64> = Vector2::new(3.0, 6.0);
 //! let b: Vector2<f64> = Vector2::new(-2.0, 1.0);
-//! assert_eq!(a.dot(&b), zero());
+//! assert_eq!(a.dot(&b), 0.0);
 //!
 //! // But there is also a top-level function:
 //! assert_eq!(a.dot(&b), dot(a, b));
@@ -112,9 +111,32 @@ use num::{BaseNum, BaseFloat};
 /// A trait that specifies a range of numeric operations for vectors. Not all
 /// of these make sense from a linear algebra point of view, but are included
 /// for pragmatic reasons.
-pub trait Vector<S: BaseNum>: Array1<S> + Zero + One {
+pub trait Vector<S: BaseNum>: Array1<S> + Clone // where
+    // FIXME: blocked by rust-lang/rust#20671
+    //
+    // for<'a, 'b> &'a Self: Add<&'b Self, Output = Self>,
+    // for<'a, 'b> &'a Self: Sub<&'b Self, Output = Self>,
+    // for<'a, 'b> &'a Self: Mul<&'b Self, Output = Self>,
+    // for<'a, 'b> &'a Self: Div<&'b Self, Output = Self>,
+    // for<'a, 'b> &'a Self: Rem<&'b Self, Output = Self>,
+    // for<'a, 'b> &'a Self: Sub<&'b Self, Output = Self>,
+    //
+    // for<'a> &'a Self: Add<S, Output = Self>,
+    // for<'a> &'a Self: Sub<S, Output = Self>,
+    // for<'a> &'a Self: Mul<S, Output = Self>,
+    // for<'a> &'a Self: Div<S, Output = Self>,
+    // for<'a> &'a Self: Rem<S, Output = Self>,
+{
     /// Construct a vector from a single value, replicating it.
     fn from_value(s: S) -> Self;
+
+    /// The zero vector (with all components set to zero)
+    #[inline]
+    fn zero() -> Self { Self::from_value(S::zero()) }
+    /// The identity vector (with all components set to one)
+    #[inline]
+    fn identity() -> Self { Self::from_value(S::one()) }
+
     /// Add a scalar to this vector, returning a new vector.
     #[must_use]
     fn add_s(&self, s: S) -> Self;
@@ -215,19 +237,6 @@ macro_rules! vec {
             $Self_::new($($field),+)
         }
 
-        impl<$S: Zero + BaseNum> Zero for $Self_<$S> {
-            #[inline]
-            fn zero() -> $Self_<S> { $Self_ { $($field: $S::zero()),+ } }
-
-            #[inline]
-            fn is_zero(&self) -> bool { $((self.$field.is_zero()) )&&+ }
-        }
-
-        impl<$S: One + BaseNum> One for $Self_<$S> {
-            #[inline]
-            fn one() -> $Self_<$S> { $Self_ { $($field: S::one()),+ } }
-        }
-
         impl<$S: NumCast + Copy> $Self_<$S> {
             /// Component-wise casting to another type
             #[inline]
@@ -240,29 +249,30 @@ macro_rules! vec {
 
         impl<S: BaseNum> Vector<S> for $Self_<S> {
             #[inline] fn from_value(s: S) -> $Self_<S> { $Self_ { $($field: s),+ } }
-            #[inline] fn add_s(&self, s: S) -> $Self_<S> { $Self_::new($(self.$field + s),+) }
-            #[inline] fn sub_s(&self, s: S) -> $Self_<S> { $Self_::new($(self.$field - s),+) }
-            #[inline] fn mul_s(&self, s: S) -> $Self_<S> { $Self_::new($(self.$field * s),+) }
-            #[inline] fn div_s(&self, s: S) -> $Self_<S> { $Self_::new($(self.$field / s),+) }
-            #[inline] fn rem_s(&self, s: S) -> $Self_<S> { $Self_::new($(self.$field % s),+) }
 
-            #[inline] fn add_v(&self, v: &$Self_<S>) -> $Self_<S> { $Self_::new($(self.$field + v.$field),+) }
-            #[inline] fn sub_v(&self, v: &$Self_<S>) -> $Self_<S> { $Self_::new($(self.$field - v.$field),+) }
-            #[inline] fn mul_v(&self, v: &$Self_<S>) -> $Self_<S> { $Self_::new($(self.$field * v.$field),+) }
-            #[inline] fn div_v(&self, v: &$Self_<S>) -> $Self_<S> { $Self_::new($(self.$field / v.$field),+) }
-            #[inline] fn rem_v(&self, v: &$Self_<S>) -> $Self_<S> { $Self_::new($(self.$field % v.$field),+) }
+            #[inline] fn add_s(&self, s: S) -> $Self_<S> { self + s }
+            #[inline] fn sub_s(&self, s: S) -> $Self_<S> { self - s }
+            #[inline] fn mul_s(&self, s: S) -> $Self_<S> { self * s }
+            #[inline] fn div_s(&self, s: S) -> $Self_<S> { self / s }
+            #[inline] fn rem_s(&self, s: S) -> $Self_<S> { self % s }
 
-            #[inline] fn add_self_s(&mut self, s: S) { $(self.$field = self.$field + s;)+ }
-            #[inline] fn sub_self_s(&mut self, s: S) { $(self.$field = self.$field - s;)+ }
-            #[inline] fn mul_self_s(&mut self, s: S) { $(self.$field = self.$field * s;)+ }
-            #[inline] fn div_self_s(&mut self, s: S) { $(self.$field = self.$field / s;)+ }
-            #[inline] fn rem_self_s(&mut self, s: S) { $(self.$field = self.$field % s;)+ }
+            #[inline] fn add_v(&self, v: &$Self_<S>) -> $Self_<S> { self + v }
+            #[inline] fn sub_v(&self, v: &$Self_<S>) -> $Self_<S> { self - v }
+            #[inline] fn mul_v(&self, v: &$Self_<S>) -> $Self_<S> { self * v }
+            #[inline] fn div_v(&self, v: &$Self_<S>) -> $Self_<S> { self / v }
+            #[inline] fn rem_v(&self, v: &$Self_<S>) -> $Self_<S> { self % v }
 
-            #[inline] fn add_self_v(&mut self, v: &$Self_<S>) { $(self.$field = self.$field + v.$field;)+ }
-            #[inline] fn sub_self_v(&mut self, v: &$Self_<S>) { $(self.$field = self.$field - v.$field;)+ }
-            #[inline] fn mul_self_v(&mut self, v: &$Self_<S>) { $(self.$field = self.$field * v.$field;)+ }
-            #[inline] fn div_self_v(&mut self, v: &$Self_<S>) { $(self.$field = self.$field / v.$field;)+ }
-            #[inline] fn rem_self_v(&mut self, v: &$Self_<S>) { $(self.$field = self.$field % v.$field;)+ }
+            #[inline] fn add_self_s(&mut self, s: S) { *self = &*self + s; }
+            #[inline] fn sub_self_s(&mut self, s: S) { *self = &*self - s; }
+            #[inline] fn mul_self_s(&mut self, s: S) { *self = &*self * s; }
+            #[inline] fn div_self_s(&mut self, s: S) { *self = &*self / s; }
+            #[inline] fn rem_self_s(&mut self, s: S) { *self = &*self % s; }
+
+            #[inline] fn add_self_v(&mut self, v: &$Self_<S>) { *self = &*self + v; }
+            #[inline] fn sub_self_v(&mut self, v: &$Self_<S>) { *self = &*self - v; }
+            #[inline] fn mul_self_v(&mut self, v: &$Self_<S>) { *self = &*self * v; }
+            #[inline] fn div_self_v(&mut self, v: &$Self_<S>) { *self = &*self / v; }
+            #[inline] fn rem_self_v(&mut self, v: &$Self_<S>) { *self = &*self % v; }
 
             #[inline] fn comp_add(&self) -> S { fold!(add, { $(self.$field),+ }) }
             #[inline] fn comp_mul(&self) -> S { fold!(mul, { $(self.$field),+ }) }
@@ -270,20 +280,6 @@ macro_rules! vec {
             #[inline] fn comp_max(&self) -> S { fold!(partial_max, { $(self.$field),+ }) }
         }
 
-        impl<S: BaseNum> Add for $Self_<S> {
-            type Output = $Self_<S>;
-
-            #[inline]
-            fn add(self, v: $Self_<S>) -> $Self_<S> { self.add_v(&v) }
-        }
-
-        impl<S: BaseNum> Sub for $Self_<S> {
-            type Output = $Self_<S>;
-
-            #[inline]
-            fn sub(self, v: $Self_<S>) -> $Self_<S> { self.sub_v(&v) }
-        }
-
         impl<S: Neg<Output = S>> Neg for $Self_<S> {
             type Output = $Self_<S>;
 
@@ -291,27 +287,6 @@ macro_rules! vec {
             fn neg(self) -> $Self_<S> { $Self_::new($(-self.$field),+) }
         }
 
-        impl<S: BaseNum> Mul for $Self_<S> {
-            type Output = $Self_<S>;
-
-            #[inline]
-            fn mul(self, v: $Self_<S>) -> $Self_<S> { self.mul_v(&v) }
-        }
-
-        impl<S: BaseNum> Div for $Self_<S> {
-            type Output = $Self_<S>;
-
-            #[inline]
-            fn div(self, v: $Self_<S>) -> $Self_<S> { self.div_v(&v) }
-        }
-
-        impl<S: BaseNum> Rem for $Self_<S> {
-            type Output = $Self_<S>;
-
-            #[inline]
-            fn rem(self, v: $Self_<S>) -> $Self_<S> { self.rem_v(&v) }
-        }
-
         impl<S: BaseFloat> ApproxEq<S> for $Self_<S> {
             #[inline]
             fn approx_eq_eps(&self, other: &$Self_<S>, epsilon: &S) -> bool {
@@ -328,6 +303,44 @@ macro_rules! vec {
     }
 }
 
+macro_rules! impl_binary_operator {
+    ($Binop:ident :: $binop:ident, $VectorN:ident { $($field:ident),+ }) => {
+        impl<'a, S: BaseNum> $Binop<S> for &'a $VectorN<S> {
+            type Output = $VectorN<S>;
+
+            #[inline]
+            fn $binop(self, s: S) -> $VectorN<S> {
+                $VectorN::new($(self.$field.$binop(s)),+)
+            }
+        }
+
+        impl<'a, 'b, S: BaseNum> $Binop<&'a $VectorN<S>> for &'b $VectorN<S> {
+            type Output = $VectorN<S>;
+
+            #[inline]
+            fn $binop(self, other: &'a $VectorN<S>) -> $VectorN<S> {
+                $VectorN::new($(self.$field.$binop(other.$field)),+)
+            }
+        }
+    }
+}
+
+impl_binary_operator!(Add::add, Vector2 { x, y });
+impl_binary_operator!(Add::add, Vector3 { x, y, z });
+impl_binary_operator!(Add::add, Vector4 { x, y, z, w });
+impl_binary_operator!(Sub::sub, Vector2 { x, y });
+impl_binary_operator!(Sub::sub, Vector3 { x, y, z });
+impl_binary_operator!(Sub::sub, Vector4 { x, y, z, w });
+impl_binary_operator!(Mul::mul, Vector2 { x, y });
+impl_binary_operator!(Mul::mul, Vector3 { x, y, z });
+impl_binary_operator!(Mul::mul, Vector4 { x, y, z, w });
+impl_binary_operator!(Div::div, Vector2 { x, y });
+impl_binary_operator!(Div::div, Vector3 { x, y, z });
+impl_binary_operator!(Div::div, Vector4 { x, y, z, w });
+impl_binary_operator!(Rem::rem, Vector2 { x, y });
+impl_binary_operator!(Rem::rem, Vector3 { x, y, z });
+impl_binary_operator!(Rem::rem, Vector4 { x, y, z, w });
+
 macro_rules! fold {
     (&$method:ident, { $x:expr, $y:expr })                   => { $x.$method(&$y) };
     (&$method:ident, { $x:expr, $y:expr, $z:expr })          => { $x.$method(&$y).$method(&$z) };
diff --git a/tests/matrix.rs b/tests/matrix.rs
index d087a02..881d615 100644
--- a/tests/matrix.rs
+++ b/tests/matrix.rs
@@ -134,7 +134,7 @@ fn test_add_m() {
     let mut mut_a = matrix2::A;
     mut_a.add_self_m(&matrix2::B);
     assert_eq!(mut_a, matrix2::A.add_m(&matrix2::B));
-    assert_eq!(mut_a, matrix2::A + matrix2::B);
+    assert_eq!(mut_a, &matrix2::A + &matrix2::B);
 
     // Matrix3
     assert_eq!(matrix3::A.add_m(&matrix3::B),
@@ -144,7 +144,7 @@ fn test_add_m() {
     let mut mut_a = matrix3::A;
     mut_a.add_self_m(&matrix3::B);
     assert_eq!(mut_a, matrix3::A.add_m(&matrix3::B));
-    assert_eq!(mut_a, matrix3::A + matrix3::B);
+    assert_eq!(mut_a, &matrix3::A + &matrix3::B);
 
     // Matrix4
     assert_eq!(matrix4::A.add_m(&matrix4::B),
@@ -155,7 +155,7 @@ fn test_add_m() {
     let mut mut_a = matrix4::A;
     mut_a.add_self_m(&matrix4::B);
     assert_eq!(mut_a, matrix4::A.add_m(&matrix4::B));
-    assert_eq!(mut_a, matrix4::A + matrix4::B);
+    assert_eq!(mut_a, &matrix4::A + &matrix4::B);
 }
 
 #[test]
@@ -167,7 +167,7 @@ fn test_sub_m() {
     let mut mut_a = matrix2::A;
     mut_a.sub_self_m(&matrix2::B);
     assert_eq!(mut_a, matrix2::A.sub_m(&matrix2::B));
-    assert_eq!(matrix2::A.sub_m(&matrix2::B), matrix2::A - matrix2::B);
+    assert_eq!(matrix2::A.sub_m(&matrix2::B), &matrix2::A - &matrix2::B);
 
     // Matrix3
     assert_eq!(matrix3::A.sub_m(&matrix3::B),
@@ -177,7 +177,7 @@ fn test_sub_m() {
     let mut mut_a = matrix3::A;
     mut_a.sub_self_m(&matrix3::B);
     assert_eq!(mut_a, matrix3::A.sub_m(&matrix3::B));
-    assert_eq!(matrix3::A.sub_m(&matrix3::B), matrix3::A - matrix3::B);
+    assert_eq!(matrix3::A.sub_m(&matrix3::B), &matrix3::A - &matrix3::B);
 
     // Matrix4
     assert_eq!(matrix4::A.sub_m(&matrix4::B),
@@ -188,7 +188,7 @@ fn test_sub_m() {
     let mut mut_a = matrix4::A;
     mut_a.sub_self_m(&matrix4::B);
     assert_eq!(mut_a, matrix4::A.sub_m(&matrix4::B));
-    assert_eq!(matrix4::A.sub_m(&matrix4::B), matrix4::A - matrix4::B);
+    assert_eq!(matrix4::A.sub_m(&matrix4::B), &matrix4::A - &matrix4::B);
 }
 
 #[test]
@@ -213,9 +213,9 @@ fn test_mul_m() {
                             120.0f64, 280.0f64, 440.0f64, 600.0f64,
                             130.0f64, 306.0f64, 482.0f64, 658.0f64));
 
-    assert_eq!(matrix2::A.mul_m(&matrix2::B), matrix2::A * matrix2::B);
-    assert_eq!(matrix3::A.mul_m(&matrix3::B), matrix3::A * matrix3::B);
-    assert_eq!(matrix4::A.mul_m(&matrix4::B), matrix4::A * matrix4::B);
+    assert_eq!(matrix2::A.mul_m(&matrix2::B), &matrix2::A * &matrix2::B);
+    assert_eq!(matrix3::A.mul_m(&matrix3::B), &matrix3::A * &matrix3::B);
+    assert_eq!(matrix4::A.mul_m(&matrix4::B), &matrix4::A * &matrix4::B);
 }
 
 #[test]