diff --git a/src/test/test_vec.rs b/src/test/test_vec.rs
index 7b4dcf7..28abb7b 100644
--- a/src/test/test_vec.rs
+++ b/src/test/test_vec.rs
@@ -12,6 +12,8 @@ fn test_Vec2() {
     let f1 = 1.5f;
     let f2 = 0.5f;
     
+    let mut mut_a = a;
+    
     assert Vec2::new(1f, 2f) == a;
     // assert Vec2::from_value(1f32) == Vec2::new(1f32, 1f32);
     
@@ -20,10 +22,15 @@ fn test_Vec2() {
     // assert Vec2::unit_y()   == Vec2::new(0f, 1f);
     // assert Vec2::identity() == Vec2::new(1f, 1f);
     
-    let mut mut_a = a;
+    *mut_a.index_mut(0) = 42f;
+    *mut_a.index_mut(1) = 43f;
+    assert mut_a == Vec2::new(42f, 43f);
+    mut_a = a;
+    
     mut_a.swap(0, 1);
     assert mut_a[0] == a[1];
     assert mut_a[1] == a[0];
+    mut_a = a;
     
     assert a.x == 1f;
     assert a.y == 2f;
@@ -39,6 +46,26 @@ fn test_Vec2() {
     assert a.add_v(&b) == Vec2::new( 4f,  6f);
     assert a.sub_v(&b) == Vec2::new(-2f, -2f);
     
+    mut_a.neg_self();
+    assert mut_a == -a;
+    mut_a = a;
+    
+    mut_a.mul_self_t(&f1);
+    assert mut_a == a.mul_t(f1);
+    mut_a = a;
+    
+    mut_a.div_self_t(&f2);
+    assert mut_a == a.div_t(f2);
+    mut_a = a;
+    
+    mut_a.add_self_v(&b);
+    assert mut_a == a.add_v(&b);
+    mut_a = a;
+    
+    mut_a.sub_self_v(&b);
+    assert mut_a == a.sub_v(&b);
+    // mut_a = a;
+    
     // exact_eq
     // fuzzy_eq
     // eq
@@ -78,6 +105,8 @@ fn test_Vec3() {
     let f1 = 1.5f;
     let f2 = 0.5f;
     
+    let mut mut_a = a;
+    
     assert Vec3::new(1f, 2f, 3f) == a;
     // assert Vec3::from_value(1f32) == Vec3::new(1f32, 1f32, 1f32);
     
@@ -87,15 +116,21 @@ fn test_Vec3() {
     // assert Vec3::unit_z()   == Vec3::new(0f, 0f, 1f);
     // assert Vec3::identity() == Vec3::new(1f, 1f, 1f);
     
-    let mut mut_a = a;
+    *mut_a.index_mut(0) = 42f;
+    *mut_a.index_mut(1) = 43f;
+    *mut_a.index_mut(2) = 44f;
+    assert mut_a == Vec3::new(42f, 43f, 44f);
+    mut_a = a;
+    
     mut_a.swap(0, 2);
     assert mut_a[0] == a[2];
     assert mut_a[2] == a[0];
-    
     mut_a = a;
+    
     mut_a.swap(1, 2);
     assert mut_a[1] == a[2];
     assert mut_a[2] == a[1];
+    let mut mut_a = a;
     
     assert a.x == 1f;
     assert a.y == 2f;
@@ -115,6 +150,26 @@ fn test_Vec3() {
     assert a.add_v(&b) == Vec3::new( 5f,  7f,  9f);
     assert a.sub_v(&b) == Vec3::new(-3f, -3f, -3f);
     
+    mut_a.neg_self();
+    assert mut_a == -a;
+    mut_a = a;
+    
+    mut_a.mul_self_t(&f1);
+    assert mut_a == a.mul_t(f1);
+    mut_a = a;
+    
+    mut_a.div_self_t(&f2);
+    assert mut_a == a.div_t(f2);
+    mut_a = a;
+    
+    mut_a.add_self_v(&b);
+    assert mut_a == a.add_v(&b);
+    mut_a = a;
+    
+    mut_a.sub_self_v(&b);
+    assert mut_a == a.sub_v(&b);
+    // mut_a = a;
+    
     // exact_eq
     // fuzzy_eq
     // eq
@@ -154,18 +209,27 @@ fn test_Vec4() {
     let f1 = 1.5f;
     let f2 = 0.5f;
     
+    let mut mut_a = a;
+    
     assert Vec4::new(1f, 2f, 3f, 4f) == a;
     // assert Vec4::from_value(1f32) == Vec4::new(1f32, 1f32, 1f32, 1f32);
     
-    let mut mut_a = a;
+    *mut_a.index_mut(0) = 42f;
+    *mut_a.index_mut(1) = 43f;
+    *mut_a.index_mut(2) = 44f;
+    *mut_a.index_mut(3) = 45f;
+    assert mut_a == Vec4::new(42f, 43f, 44f, 45f);
+    mut_a = a;
+    
     mut_a.swap(0, 3);
     assert mut_a[0] == a[3];
     assert mut_a[3] == a[0];
-    
     mut_a = a;
+    
     mut_a.swap(1, 2);
     assert mut_a[1] == a[2];
     assert mut_a[2] == a[1];
+    mut_a = a;
     
     // assert Vec4::zero()     == Vec4::new(0f, 0f, 0f, 0f);
     // assert Vec4::unit_x()   == Vec4::new(1f, 0f, 0f, 0f);
@@ -194,6 +258,26 @@ fn test_Vec4() {
     
     assert a.dot(&b) == 70f;
     
+    mut_a.neg_self();
+    assert mut_a == -a;
+    mut_a = a;
+    
+    mut_a.mul_self_t(&f1);
+    assert mut_a == a.mul_t(f1);
+    mut_a = a;
+    
+    mut_a.div_self_t(&f2);
+    assert mut_a == a.div_t(f2);
+    mut_a = a;
+    
+    mut_a.add_self_v(&b);
+    assert mut_a == a.add_v(&b);
+    mut_a = a;
+    
+    mut_a.sub_self_v(&b);
+    assert mut_a == a.sub_v(&b);
+    // mut_a = a;
+    
     // exact_eq
     // fuzzy_eq
     // eq
diff --git a/src/vec.rs b/src/vec.rs
index 925d128..7220180 100644
--- a/src/vec.rs
+++ b/src/vec.rs
@@ -72,7 +72,7 @@ pub trait NumericVector<T>: Vector<T>, Neg<self> {
     pure fn mul_t(&self, value: T) -> self;
     
     /**
-     * Returns the scalar quotient of the vector and `value`
+     * Returns the scalar division of the vector and `value`
      */
     pure fn div_t(&self, value: T) -> self;
     
@@ -81,11 +81,6 @@ pub trait NumericVector<T>: Vector<T>, Neg<self> {
      */
     pure fn add_v(&self, other: &self) -> self;
     
-    /**
-     * Add the vector `other`
-     */
-    // fn add_v_self(&mut self, other: &self);
-    
     /**
      * Returns the difference between the vector and `other`
      */
@@ -97,6 +92,36 @@ pub trait NumericVector<T>: Vector<T>, Neg<self> {
     pure fn dot(&self, other: &self) -> T;
 }
 
+/**
+ * A mutable vector with numeric components
+ */
+pub trait MutableNumericVector<T>: MutableVector<&self/T>, NumericVector<T> {
+    /**
+     * Negate the vector
+     */
+    fn neg_self(&mut self);
+    
+    /**
+     * Multiply the vector by a scalar
+     */
+    fn mul_self_t(&mut self, value: T);
+    
+    /**
+     * Divide the vector by a scalar
+     */
+    fn div_self_t(&mut self, value: T);
+    
+    /**
+     * Add a vector
+     */
+    fn add_self_v(&mut self, other: &self);
+    
+    /**
+     * Subtract a vector
+     */
+    fn sub_self_v(&mut self, other: &self);
+}
+
 /**
  * A 2-dimensional vector with numeric components
  */
@@ -279,7 +304,39 @@ pub impl<T:Copy Number> Vec2<T>: Neg<Vec2<T>> {
         Vec2::new(-self[0], -self[1])
     }
 }
+
+pub impl<T:Copy Number> Vec2<T>: MutableNumericVector<&self/T> {
+    #[inline(always)]
+    fn neg_self(&mut self) {
+        *self.index_mut(0) = -*self.index_mut(0);
+        *self.index_mut(1) = -*self.index_mut(1);
+    }
     
+    #[inline(always)]
+    fn mul_self_t(&mut self, value: &T) {
+        *self.index_mut(0) *= (*value);
+        *self.index_mut(1) *= (*value);
+    }
+    
+    #[inline(always)]
+    fn div_self_t(&mut self, value: &T) {
+        *self.index_mut(0) /= (*value);
+        *self.index_mut(1) /= (*value);
+    }
+    
+    #[inline(always)]
+    fn add_self_v(&mut self, other: &Vec2<T>) {
+        *self.index_mut(0) += other[0];
+        *self.index_mut(1) += other[1];
+    }
+    
+    #[inline(always)]
+    fn sub_self_v(&mut self, other: &Vec2<T>) {
+        *self.index_mut(0) -= other[0];
+        *self.index_mut(1) -= other[1];
+    }
+}
+
 pub impl<T:Copy Number Exp> Vec2<T>: GeometricVector<T> {
     #[inline(always)]
     pure fn length2(&self) -> T {
@@ -477,6 +534,43 @@ pub impl<T:Copy Number> Vec3<T>: NumericVector3<T> {
     }
 }
 
+pub impl<T:Copy Number> Vec3<T>: MutableNumericVector<&self/T> {
+    #[inline(always)]
+    fn neg_self(&mut self) {
+        *self.index_mut(0) = -*self.index_mut(0);
+        *self.index_mut(1) = -*self.index_mut(1);
+        *self.index_mut(2) = -*self.index_mut(2);
+    }
+    
+    #[inline(always)]
+    fn mul_self_t(&mut self, value: &T) {
+        *self.index_mut(0) *= (*value);
+        *self.index_mut(1) *= (*value);
+        *self.index_mut(2) *= (*value);
+    }
+    
+    #[inline(always)]
+    fn div_self_t(&mut self, value: &T) {
+        *self.index_mut(0) /= (*value);
+        *self.index_mut(1) /= (*value);
+        *self.index_mut(2) /= (*value);
+    }
+    
+    #[inline(always)]
+    fn add_self_v(&mut self, other: &Vec3<T>) {
+        *self.index_mut(0) += other[0];
+        *self.index_mut(1) += other[1];
+        *self.index_mut(2) += other[2];
+    }
+    
+    #[inline(always)]
+    fn sub_self_v(&mut self, other: &Vec3<T>) {
+        *self.index_mut(0) -= other[0];
+        *self.index_mut(1) -= other[1];
+        *self.index_mut(2) -= other[2];
+    }
+}
+
 pub impl<T:Copy Number Exp> Vec3<T>: GeometricVector<T> {
     #[inline(always)]
     pure fn length2(&self) -> T {
@@ -675,6 +769,48 @@ pub impl<T:Copy Number> Vec4<T>: Neg<Vec4<T>> {
     }
 }
 
+pub impl<T:Copy Number> Vec4<T>: MutableNumericVector<&self/T> {
+    #[inline(always)]
+    fn neg_self(&mut self) {
+        *self.index_mut(0) = -*self.index_mut(0);
+        *self.index_mut(1) = -*self.index_mut(1);
+        *self.index_mut(2) = -*self.index_mut(2);
+        *self.index_mut(3) = -*self.index_mut(3);
+    }
+    
+    #[inline(always)]
+    fn mul_self_t(&mut self, value: &T) {
+        *self.index_mut(0) *= (*value);
+        *self.index_mut(1) *= (*value);
+        *self.index_mut(2) *= (*value);
+        *self.index_mut(3) *= (*value);
+    }
+    
+    #[inline(always)]
+    fn div_self_t(&mut self, value: &T) {
+        *self.index_mut(0) /= (*value);
+        *self.index_mut(1) /= (*value);
+        *self.index_mut(2) /= (*value);
+        *self.index_mut(3) /= (*value);
+    }
+    
+    #[inline(always)]
+    fn add_self_v(&mut self, other: &Vec4<T>) {
+        *self.index_mut(0) += other[0];
+        *self.index_mut(1) += other[1];
+        *self.index_mut(2) += other[2];
+        *self.index_mut(3) += other[3];
+    }
+    
+    #[inline(always)]
+    fn sub_self_v(&mut self, other: &Vec4<T>) {
+        *self.index_mut(0) -= other[0];
+        *self.index_mut(1) -= other[1];
+        *self.index_mut(2) -= other[2];
+        *self.index_mut(3) -= other[3];
+    }
+}
+
 pub impl<T:Copy Number Exp> Vec4<T>: GeometricVector<T> {
     #[inline(always)]
     pure fn length2(&self) -> T {