diff --git a/README.md b/README.md
index a986f1e4..528dc563 100644
--- a/README.md
+++ b/README.md
@@ -5,9 +5,9 @@ nalgebra
 
 **nalgebra** is a low-dimensional linear algebra library written for Rust targeting:
 
-* general-purpose linear algebra (still lacks a lot of features…).
-* real time computer graphics.
-* real time computer physics.
+* General-purpose linear algebra (still lacks a lot of features…)
+* Real time computer graphics.
+* Real time computer physics.
 
 An on-line version of this documentation is available [here](http://nalgebra.org/doc/nalgebra).
 
diff --git a/src/structs/dmat.rs b/src/structs/dmat.rs
index b99b25f1..310c7532 100644
--- a/src/structs/dmat.rs
+++ b/src/structs/dmat.rs
@@ -3,7 +3,7 @@
 use std::cmp;
 use std::mem;
 use std::iter::repeat;
-use std::ops::{Add, Sub, Mul, Div, Index, IndexMut};
+use std::ops::{Add, Sub, Mul, Div, AddAssign, SubAssign, MulAssign, Index, IndexMut};
 use std::fmt::{Debug, Formatter, Result};
 use rand::{self, Rand};
 use num::{Zero, One};
diff --git a/src/structs/dmat_macros.rs b/src/structs/dmat_macros.rs
index 8cd70f34..aaf02a2c 100644
--- a/src/structs/dmat_macros.rs
+++ b/src/structs/dmat_macros.rs
@@ -165,7 +165,13 @@ macro_rules! dmat_impl(
             }
         }
 
-        impl<N: Copy + Mul<N, Output = N> + Add<N, Output = N> + Zero> Mul<$dmat<N>> for $dmat<N> {
+        /*
+         *
+         * Multiplications matrix/matrix.
+         *
+         */
+        impl<N> Mul<$dmat<N>> for $dmat<N>
+            where N: Copy + Mul<N, Output = N> + Add<N, Output = N> + Zero {
             type Output = $dmat<N>;
 
             #[inline]
@@ -174,7 +180,8 @@ macro_rules! dmat_impl(
             }
         }
 
-        impl<'a, N: Copy + Mul<N, Output = N> + Add<N, Output = N> + Zero> Mul<&'a $dmat<N>> for $dmat<N> {
+        impl<'a, N> Mul<&'a $dmat<N>> for $dmat<N>
+            where N: Copy + Mul<N, Output = N> + Add<N, Output = N> + Zero {
             type Output = $dmat<N>;
 
             #[inline]
@@ -183,7 +190,8 @@ macro_rules! dmat_impl(
             }
         }
 
-        impl<'a, N: Copy + Mul<N, Output = N> + Add<N, Output = N> + Zero> Mul<$dmat<N>> for &'a $dmat<N> {
+        impl<'a, N> Mul<$dmat<N>> for &'a $dmat<N>
+            where N: Copy + Mul<N, Output = N> + Add<N, Output = N> + Zero {
             type Output = $dmat<N>;
 
             #[inline]
@@ -192,7 +200,8 @@ macro_rules! dmat_impl(
             }
         }
 
-        impl<'a, 'b, N: Copy + Mul<N, Output = N> + Add<N, Output = N> + Zero> Mul<&'b $dmat<N>> for &'a $dmat<N> {
+        impl<'a, 'b, N> Mul<&'b $dmat<N>> for &'a $dmat<N>
+            where N: Copy + Mul<N, Output = N> + Add<N, Output = N> + Zero {
             type Output = $dmat<N>;
 
             #[inline]
@@ -201,14 +210,13 @@ macro_rules! dmat_impl(
 
                 let mut res = unsafe { $dmat::new_uninitialized(self.nrows, right.ncols) };
 
-                for i in 0..self.nrows {
-                    for j in 0..right.ncols {
+                for i in 0 .. self.nrows {
+                    for j in 0 .. right.ncols {
                         let mut acc: N = ::zero();
 
                         unsafe {
-                            for k in 0..self.ncols {
-                                acc = acc
-                                    + self.unsafe_at((i, k)) * right.unsafe_at((k, j));
+                            for k in 0 .. self.ncols {
+                                acc = acc + self.unsafe_at((i, k)) * right.unsafe_at((k, j));
                             }
 
                             res.unsafe_set((i, j), acc);
@@ -220,10 +228,64 @@ macro_rules! dmat_impl(
             }
         }
 
-        impl<N: Copy + Add<N, Output = N> + Mul<N, Output = N> + Zero> Mul<$dvec<N>> for $dmat<N> {
+        impl<N> MulAssign<$dmat<N>> for $dmat<N>
+            where N: Copy + Mul<N, Output = N> + Add<N, Output = N> + Zero {
+            #[inline]
+            fn mul_assign(&mut self, right: $dmat<N>) {
+                self.mul_assign(&right)
+            }
+        }
+
+        impl<'a, N> MulAssign<&'a $dmat<N>> for $dmat<N>
+            where N: Copy + Mul<N, Output = N> + Add<N, Output = N> + Zero {
+            #[inline]
+            fn mul_assign(&mut self, right: &'a $dmat<N>) {
+                assert!(self.ncols == right.nrows);
+
+                // FIXME: optimize when both matrices have the same layout.
+                let res = &*self * right;
+                *self = res;
+            }
+        }
+
+
+        /*
+         *
+         * Multiplication matrix/vector.
+         *
+         */
+        impl<N> Mul<$dvec<N>> for $dmat<N>
+            where N: Copy + Add<N, Output = N> + Mul<N, Output = N> + Zero {
             type Output = $dvec<N>;
 
             fn mul(self, right: $dvec<N>) -> $dvec<N> {
+                (&self) * (&right)
+            }
+        }
+
+        impl<'a, N> Mul<$dvec<N>> for &'a $dmat<N>
+            where N: Copy + Add<N, Output = N> + Mul<N, Output = N> + Zero {
+            type Output = $dvec<N>;
+
+            fn mul(self, right: $dvec<N>) -> $dvec<N> {
+                self * (&right)
+            }
+        }
+
+        impl<'a, N> Mul<&'a $dvec<N>> for $dmat<N>
+            where N: Copy + Add<N, Output = N> + Mul<N, Output = N> + Zero {
+            type Output = $dvec<N>;
+
+            fn mul(self, right: &'a $dvec<N>) -> $dvec<N> {
+                (&self) * right
+            }
+        }
+
+        impl<'a, 'b, N> Mul<&'b $dvec<N>> for &'a $dmat<N>
+            where N: Copy + Add<N, Output = N> + Mul<N, Output = N> + Zero {
+            type Output = $dvec<N>;
+
+            fn mul(self, right: &'b $dvec<N>) -> $dvec<N> {
                 assert!(self.ncols == right.len());
 
                 let mut res : $dvec<N> = unsafe { $dvec::new_uninitialized(self.nrows) };
@@ -246,11 +308,39 @@ macro_rules! dmat_impl(
             }
         }
 
-
-        impl<N: Copy + Add<N, Output = N> + Mul<N, Output = N> + Zero> Mul<$dmat<N>> for $dvec<N> {
+        impl<N> Mul<$dmat<N>> for $dvec<N>
+            where N: Copy + Add<N, Output = N> + Mul<N, Output = N> + Zero {
             type Output = $dvec<N>;
 
             fn mul(self, right: $dmat<N>) -> $dvec<N> {
+                (&self) * (&right)
+            }
+        }
+
+        impl<'a, N> Mul<$dmat<N>> for &'a $dvec<N>
+            where N: Copy + Add<N, Output = N> + Mul<N, Output = N> + Zero {
+            type Output = $dvec<N>;
+
+            fn mul(self, right: $dmat<N>) -> $dvec<N> {
+                self * (&right)
+            }
+        }
+
+        impl<'a, N> Mul<&'a $dmat<N>> for $dvec<N>
+            where N: Copy + Add<N, Output = N> + Mul<N, Output = N> + Zero {
+            type Output = $dvec<N>;
+
+            fn mul(self, right: &'a $dmat<N>) -> $dvec<N> {
+                (&self) * right
+            }
+        }
+
+
+        impl<'a, 'b, N> Mul<&'b $dmat<N>> for &'a $dvec<N>
+            where N: Copy + Add<N, Output = N> + Mul<N, Output = N> + Zero {
+            type Output = $dvec<N>;
+
+            fn mul(self, right: &'b $dmat<N>) -> $dvec<N> {
                 assert!(right.nrows == self.len());
 
                 let mut res : $dvec<N> = unsafe { $dvec::new_uninitialized(right.ncols) };
@@ -273,6 +363,206 @@ macro_rules! dmat_impl(
             }
         }
 
+        impl<N> MulAssign<$dmat<N>> for $dvec<N>
+            where N: Copy + Mul<N, Output = N> + Add<N, Output = N> + Zero {
+            #[inline]
+            fn mul_assign(&mut self, right: $dmat<N>) {
+                self.mul_assign(&right)
+            }
+        }
+
+        impl<'a, N> MulAssign<&'a $dmat<N>> for $dvec<N>
+            where N: Copy + Mul<N, Output = N> + Add<N, Output = N> + Zero {
+            #[inline]
+            fn mul_assign(&mut self, right: &'a $dmat<N>) {
+                assert!(right.nrows == self.len());
+
+                let res = &*self * right;
+                *self = res;
+            }
+        }
+
+        /*
+         *
+         * Addition matrix/matrix.
+         *
+         */
+        impl<N: Copy + Add<N, Output = N>> Add<$dmat<N>> for $dmat<N> {
+            type Output = $dmat<N>;
+
+            #[inline]
+            fn add(self, right: $dmat<N>) -> $dmat<N> {
+                self + (&right)
+            }
+        }
+
+        impl<'a, N: Copy + Add<N, Output = N>> Add<$dmat<N>> for &'a $dmat<N> {
+            type Output = $dmat<N>;
+
+            #[inline]
+            fn add(self, right: $dmat<N>) -> $dmat<N> {
+                right + self
+            }
+        }
+
+        impl<'a, N: Copy + Add<N, Output = N>> Add<&'a $dmat<N>> for $dmat<N> {
+            type Output = $dmat<N>;
+
+            #[inline]
+            fn add(self, right: &'a $dmat<N>) -> $dmat<N> {
+                let mut res = self;
+
+                for (mij, right_ij) in res.mij.iter_mut().zip(right.mij.iter()) {
+                    *mij = *mij + *right_ij;
+                }
+
+                res
+            }
+        }
+
+        impl<N: Copy + AddAssign<N>> AddAssign<$dmat<N>> for $dmat<N> {
+            #[inline]
+            fn add_assign(&mut self, right: $dmat<N>) {
+                self.add_assign(&right)
+            }
+        }
+
+        impl<'a, N: Copy + AddAssign<N>> AddAssign<&'a $dmat<N>> for $dmat<N> {
+            #[inline]
+            fn add_assign(&mut self, right: &'a $dmat<N>) {
+                assert!(self.nrows == right.nrows && self.ncols == right.ncols,
+                        "Unable to add matrices with different dimensions.");
+
+                for (mij, right_ij) in self.mij.iter_mut().zip(right.mij.iter()) {
+                    *mij += *right_ij;
+                }
+            }
+        }
+
+        /*
+         *
+         * Subtraction matrix/scalar.
+         *
+         */
+        impl<N: Copy + Sub<N, Output = N>> Sub<N> for $dmat<N> {
+            type Output = $dmat<N>;
+
+            #[inline]
+            fn sub(self, right: N) -> $dmat<N> {
+                let mut res = self;
+
+                for mij in res.mij.iter_mut() {
+                    *mij = *mij - right;
+                }
+
+                res
+            }
+        }
+
+        impl<'a, N: Copy + SubAssign<N>> SubAssign<N> for $dmat<N> {
+            #[inline]
+            fn sub_assign(&mut self, right: N) {
+                for mij in self.mij.iter_mut() {
+                    *mij -= right
+                }
+            }
+        }
+
+        impl Sub<$dmat<f32>> for f32 {
+            type Output = $dmat<f32>;
+
+            #[inline]
+            fn sub(self, right: $dmat<f32>) -> $dmat<f32> {
+                let mut res = right;
+
+                for mij in res.mij.iter_mut() {
+                    *mij = self - *mij;
+                }
+
+                res
+            }
+        }
+
+        impl Sub<$dmat<f64>> for f64 {
+            type Output = $dmat<f64>;
+
+            #[inline]
+            fn sub(self, right: $dmat<f64>) -> $dmat<f64> {
+                let mut res = right;
+
+                for mij in res.mij.iter_mut() {
+                    *mij = self - *mij;
+                }
+
+                res
+            }
+        }
+
+        /*
+         *
+         * Subtraction matrix/matrix.
+         *
+         */
+        impl<N: Copy + Sub<N, Output = N>> Sub<$dmat<N>> for $dmat<N> {
+            type Output = $dmat<N>;
+
+            #[inline]
+            fn sub(self, right: $dmat<N>) -> $dmat<N> {
+                self - (&right)
+            }
+        }
+
+        impl<'a, N: Copy + Sub<N, Output = N>> Sub<$dmat<N>> for &'a $dmat<N> {
+            type Output = $dmat<N>;
+
+            #[inline]
+            fn sub(self, right: $dmat<N>) -> $dmat<N> {
+                right - self
+            }
+        }
+
+        impl<'a, N: Copy + Sub<N, Output = N>> Sub<&'a $dmat<N>> for $dmat<N> {
+            type Output = $dmat<N>;
+
+            #[inline]
+            fn sub(self, right: &'a $dmat<N>) -> $dmat<N> {
+                assert!(self.nrows == right.nrows && self.ncols == right.ncols,
+                        "Unable to subtract matrices with different dimensions.");
+
+                let mut res = self;
+
+                for (mij, right_ij) in res.mij.iter_mut().zip(right.mij.iter()) {
+                    *mij = *mij - *right_ij;
+                }
+
+                res
+            }
+        }
+
+        impl<N: Copy + SubAssign<N>> SubAssign<$dmat<N>> for $dmat<N> {
+            #[inline]
+            fn sub_assign(&mut self, right: $dmat<N>) {
+                self.sub_assign(&right)
+            }
+        }
+
+        impl<'a, N: Copy + SubAssign<N>> SubAssign<&'a $dmat<N>> for $dmat<N> {
+            #[inline]
+            fn sub_assign(&mut self, right: &'a $dmat<N>) {
+                assert!(self.nrows == right.nrows && self.ncols == right.ncols,
+                        "Unable to subtract matrices with different dimensions.");
+
+                for (mij, right_ij) in self.mij.iter_mut().zip(right.mij.iter()) {
+                    *mij -= *right_ij;
+                }
+            }
+        }
+
+        /*
+         *
+         * Inversion.
+         *
+         */
         impl<N: BaseNum + Clone> Inv for $dmat<N> {
             #[inline]
             fn inv(&self) -> Option<$dmat<N>> {
@@ -619,6 +909,11 @@ macro_rules! dmat_impl(
             }
         }
 
+        /*
+         *
+         * Multpilication matrix/scalar.
+         *
+         */
         impl<N: Copy + Mul<N, Output = N>> Mul<N> for $dmat<N> {
             type Output = $dmat<N>;
 
@@ -664,6 +959,11 @@ macro_rules! dmat_impl(
             }
         }
 
+        /*
+         *
+         * Division matrix/scalar.
+         *
+         */
         impl<N: Copy + Div<N, Output = N>> Div<N> for $dmat<N> {
             type Output = $dmat<N>;
 
@@ -679,6 +979,12 @@ macro_rules! dmat_impl(
             }
         }
 
+
+        /*
+         *
+         * Addition matrix/scalar.
+         *
+         */
         impl<N: Copy + Add<N, Output = N>> Add<N> for $dmat<N> {
             type Output = $dmat<N>;
 
@@ -724,123 +1030,6 @@ macro_rules! dmat_impl(
             }
         }
 
-        impl<N: Copy + Add<N, Output = N>> Add<$dmat<N>> for $dmat<N> {
-            type Output = $dmat<N>;
-
-            #[inline]
-            fn add(self, right: $dmat<N>) -> $dmat<N> {
-                self + (&right)
-            }
-        }
-
-        impl<'a, N: Copy + Add<N, Output = N>> Add<$dmat<N>> for &'a $dmat<N> {
-            type Output = $dmat<N>;
-
-            #[inline]
-            fn add(self, right: $dmat<N>) -> $dmat<N> {
-                right + self
-            }
-        }
-
-        impl<'a, N: Copy + Add<N, Output = N>> Add<&'a $dmat<N>> for $dmat<N> {
-            type Output = $dmat<N>;
-
-            #[inline]
-            fn add(self, right: &'a $dmat<N>) -> $dmat<N> {
-                assert!(self.nrows == right.nrows && self.ncols == right.ncols,
-                        "Unable to add matrices with different dimensions.");
-
-                let mut res = self;
-
-                for (mij, right_ij) in res.mij.iter_mut().zip(right.mij.iter()) {
-                    *mij = *mij + *right_ij;
-                }
-
-                res
-            }
-        }
-
-        impl<N: Copy + Sub<N, Output = N>> Sub<N> for $dmat<N> {
-            type Output = $dmat<N>;
-
-            #[inline]
-            fn sub(self, right: N) -> $dmat<N> {
-                let mut res = self;
-
-                for mij in res.mij.iter_mut() {
-                    *mij = *mij - right;
-                }
-
-                res
-            }
-        }
-
-        impl Sub<$dmat<f32>> for f32 {
-            type Output = $dmat<f32>;
-
-            #[inline]
-            fn sub(self, right: $dmat<f32>) -> $dmat<f32> {
-                let mut res = right;
-
-                for mij in res.mij.iter_mut() {
-                    *mij = self - *mij;
-                }
-
-                res
-            }
-        }
-
-        impl Sub<$dmat<f64>> for f64 {
-            type Output = $dmat<f64>;
-
-            #[inline]
-            fn sub(self, right: $dmat<f64>) -> $dmat<f64> {
-                let mut res = right;
-
-                for mij in res.mij.iter_mut() {
-                    *mij = self - *mij;
-                }
-
-                res
-            }
-        }
-
-        impl<N: Copy + Sub<N, Output = N>> Sub<$dmat<N>> for $dmat<N> {
-            type Output = $dmat<N>;
-
-            #[inline]
-            fn sub(self, right: $dmat<N>) -> $dmat<N> {
-                self - (&right)
-            }
-        }
-
-        impl<'a, N: Copy + Sub<N, Output = N>> Sub<$dmat<N>> for &'a $dmat<N> {
-            type Output = $dmat<N>;
-
-            #[inline]
-            fn sub(self, right: $dmat<N>) -> $dmat<N> {
-                right - self
-            }
-        }
-
-        impl<'a, N: Copy + Sub<N, Output = N>> Sub<&'a $dmat<N>> for $dmat<N> {
-            type Output = $dmat<N>;
-
-            #[inline]
-            fn sub(self, right: &'a $dmat<N>) -> $dmat<N> {
-                assert!(self.nrows == right.nrows && self.ncols == right.ncols,
-                        "Unable to subtract matrices with different dimensions.");
-
-                let mut res = self;
-
-                for (mij, right_ij) in res.mij.iter_mut().zip(right.mij.iter()) {
-                    *mij = *mij - *right_ij;
-                }
-
-                res
-            }
-        }
-
         #[cfg(feature="arbitrary")]
         impl<N: Copy + Zero + Arbitrary> Arbitrary for $dmat<N> {
             fn arbitrary<G: Gen>(g: &mut G) -> $dmat<N> {
diff --git a/src/structs/dvec.rs b/src/structs/dvec.rs
index bc067e47..4f31a273 100644
--- a/src/structs/dvec.rs
+++ b/src/structs/dvec.rs
@@ -3,7 +3,7 @@
 use std::slice::{Iter, IterMut};
 use std::iter::{FromIterator, IntoIterator};
 use std::iter::repeat;
-use std::ops::{Add, Sub, Mul, Div, Neg, Index, IndexMut};
+use std::ops::{Add, Sub, Mul, Div, Neg, AddAssign, SubAssign, MulAssign, DivAssign, Index, IndexMut};
 use std::mem;
 use rand::{self, Rand};
 use num::{Zero, One};
diff --git a/src/structs/iso.rs b/src/structs/iso.rs
index a96584ef..f69c8eef 100644
--- a/src/structs/iso.rs
+++ b/src/structs/iso.rs
@@ -1,5 +1,5 @@
 use std::fmt;
-use std::ops::{Add, Sub, Mul, Neg};
+use std::ops::{Add, Sub, Mul, Neg, MulAssign};
 
 use rand::{Rand, Rng};
 use num::One;
diff --git a/src/structs/iso_macros.rs b/src/structs/iso_macros.rs
index 6c3fcc37..1ba52291 100644
--- a/src/structs/iso_macros.rs
+++ b/src/structs/iso_macros.rs
@@ -73,21 +73,29 @@ macro_rules! iso_mul_iso_impl(
                     self.rotation * right.rotation)
             }
         }
+
+        impl<N: BaseFloat> MulAssign<$t<N>> for $t<N> {
+            #[inline]
+            fn mul_assign(&mut self, right: $t<N>) {
+                self.translation += self.rotation * right.translation;
+                self.rotation    *= right.rotation;
+            }
+        }
     )
 );
 
 macro_rules! iso_mul_rot_impl(
-    ($t: ident, $tr: ident) => (
-        impl<N: BaseFloat> Mul<$tr<N>> for $t<N> {
+    ($t: ident, $rot: ident) => (
+        impl<N: BaseFloat> Mul<$rot<N>> for $t<N> {
             type Output = $t<N>;
 
             #[inline]
-            fn mul(self, right: $tr<N>) -> $t<N> {
+            fn mul(self, right: $rot<N>) -> $t<N> {
                 $t::new_with_rotmat(self.translation, self.rotation * right)
             }
         }
 
-        impl<N: BaseFloat> Mul<$t<N>> for $tr<N> {
+        impl<N: BaseFloat> Mul<$t<N>> for $rot<N> {
             type Output = $t<N>;
 
             #[inline]
@@ -97,6 +105,13 @@ macro_rules! iso_mul_rot_impl(
                     self * right.rotation)
             }
         }
+
+        impl<N: BaseFloat> MulAssign<$rot<N>> for $t<N> {
+            #[inline]
+            fn mul_assign(&mut self, right: $rot<N>) {
+                self.rotation *= right
+            }
+        }
     )
 );
 
@@ -110,9 +125,6 @@ macro_rules! iso_mul_pnt_impl(
                 self.rotation * right + self.translation
             }
         }
-
-        // NOTE: there is no viable pre-multiplication definition because of the translation
-        // component.
     )
 );
 
diff --git a/src/structs/mat.rs b/src/structs/mat.rs
index 9b580206..3285e539 100644
--- a/src/structs/mat.rs
+++ b/src/structs/mat.rs
@@ -3,7 +3,7 @@
 #![allow(missing_docs)] // we allow missing to avoid having to document the mij components.
 
 use std::fmt;
-use std::ops::{Add, Sub, Mul, Div, Index, IndexMut};
+use std::ops::{Add, Sub, Mul, Div, AddAssign, SubAssign, MulAssign, DivAssign, Index, IndexMut};
 use std::mem;
 use std::slice::{Iter, IterMut};
 use rand::{Rand, Rng};
diff --git a/src/structs/mat_macros.rs b/src/structs/mat_macros.rs
index 16a7d44b..a1fff8cd 100644
--- a/src/structs/mat_macros.rs
+++ b/src/structs/mat_macros.rs
@@ -92,6 +92,13 @@ macro_rules! add_impl(
                 $t::new($(self.$compN + right.$compN),+)
             }
         }
+
+        impl<N: AddAssign<N>> AddAssign<$t<N>> for $t<N> {
+            #[inline]
+            fn add_assign(&mut self, right: $t<N>) {
+                $( self.$compN += right.$compN; )+
+            }
+        }
     )
 );
 
@@ -105,6 +112,14 @@ macro_rules! sub_impl(
                 $t::new($(self.$compN - right.$compN),+)
             }
         }
+
+
+        impl<N: SubAssign<N>> SubAssign<$t<N>> for $t<N> {
+            #[inline]
+            fn sub_assign(&mut self, right: $t<N>) {
+                $( self.$compN -= right.$compN; )+
+            }
+        }
     )
 );
 
@@ -119,6 +134,13 @@ macro_rules! mat_mul_scalar_impl(
             }
         }
 
+        impl<N: MulAssign<N>> MulAssign<N> for $t<N> {
+            #[inline]
+            fn mul_assign(&mut self, right: N) {
+                $( self.$compN *= *right; )+
+            }
+        }
+
         impl Mul<$t<f32>> for f32 {
             type Output = $t<f32>;
 
@@ -149,6 +171,13 @@ macro_rules! mat_div_scalar_impl(
                 $t::new($(self.$compN / *right),+)
             }
         }
+
+        impl<N: DivAssign<N>> DivAssign<N> for $t<N> {
+            #[inline]
+            fn div_assign(&mut self, right: N) {
+                $( self.$compN /= *right; )+
+            }
+        }
     )
 );
 
@@ -163,6 +192,13 @@ macro_rules! mat_add_scalar_impl(
             }
         }
 
+        impl<N: AddAssign<N>> AddAssign<N> for $t<N> {
+            #[inline]
+            fn add_assign(&mut self, right: N) {
+                $( self.$compN += *right; )+
+            }
+        }
+
         impl Add<$t<f32>> for f32 {
             type Output = $t<f32>;
 
@@ -183,6 +219,44 @@ macro_rules! mat_add_scalar_impl(
     )
 );
 
+macro_rules! mat_sub_scalar_impl(
+    ($t: ident, $($compN: ident),+) => (
+        impl<N: Sub<N, Output = N>> Sub<N> for $t<N> {
+            type Output = $t<N>;
+
+            #[inline]
+            fn sub(self, right: &N) -> $t<N> {
+                $t::new($(self.$compN - *right),+)
+            }
+        }
+
+        impl<N: SubAssign<N>> SubAssign<N> for $t<N> {
+            #[inline]
+            fn sub_assign(&mut self, right: N) {
+                $( self.$compN -= *right; )+
+            }
+        }
+
+        impl Sub<f32> for $t<f32> {
+            type Output = $t<f32>;
+
+            #[inline]
+            fn sub(self, right: $t<f32>) -> $t<f32> {
+                $t::new($(self - right.$compN),+)
+            }
+        }
+
+        impl Sub<f64> for $t<f64> {
+            type Output = $t<f64>;
+
+            #[inline]
+            fn sub(self, right: $t<f64>) -> $t<f64> {
+                $t::new($(self - right.$compN),+)
+            }
+        }
+    )
+);
+
 
 macro_rules! eye_impl(
     ($t: ident, $dim: expr, $($comp_diagN: ident),+) => (
@@ -209,37 +283,6 @@ macro_rules! repeat_impl(
     )
 );
 
-macro_rules! mat_sub_scalar_impl(
-    ($t: ident, $($compN: ident),+) => (
-        impl<N: Sub<N, Output = N>> Sub<N> for $t<N> {
-            type Output = $t<N>;
-
-            #[inline]
-            fn sub(self, right: &N) -> $t<N> {
-                $t::new($(self.$compN - *right),+)
-            }
-        }
-
-        impl Sub<f32> for $t<f32> {
-            type Output = $t<f32>;
-
-            #[inline]
-            fn sub(self, right: $t<f32>) -> $t<f32> {
-                $t::new($(self - right.$compN),+)
-            }
-        }
-
-        impl Sub<f64> for $t<f64> {
-            type Output = $t<f64>;
-
-            #[inline]
-            fn sub(self, right: $t<f64>) -> $t<f64> {
-                $t::new($(self - right.$compN),+)
-            }
-        }
-    )
-);
-
 macro_rules! absolute_impl(
   ($t: ident, $($compN: ident),+) => (
     impl<N: Absolute<N>> Absolute<$t<N>> for $t<N> {
@@ -468,7 +511,7 @@ macro_rules! diag_impl(
             fn diag(&self) -> $tv<N> {
                 let mut diag: $tv<N> = ::zero();
 
-                for i in 0..$dim {
+                for i in 0 .. $dim {
                     unsafe { diag.unsafe_set(i, self.unsafe_at((i, i))) }
                 }
 
@@ -479,7 +522,7 @@ macro_rules! diag_impl(
         impl<N: Copy + Zero> DiagMut<$tv<N>> for $t<N> {
             #[inline]
             fn set_diag(&mut self, diag: &$tv<N>) {
-                for i in 0..$dim {
+                for i in 0 .. $dim {
                     unsafe { self.unsafe_set((i, i), diag.unsafe_at(i)) }
                 }
             }
@@ -493,15 +536,14 @@ macro_rules! mat_mul_mat_impl(
         type Output = $t<N>;
         #[inline]
         fn mul(self, right: $t<N>) -> $t<N> {
-            // careful! we need to comute other * self here (self is the rhs).
             let mut res: $t<N> = ::zero();
 
-            for i in 0..$dim {
-                for j in 0..$dim {
+            for i in 0 .. $dim {
+                for j in 0 .. $dim {
                     let mut acc: N = ::zero();
 
                     unsafe {
-                        for k in 0..$dim {
+                        for k in 0 .. $dim {
                             acc = acc + self.at_fast((i, k)) * right.at_fast((k, j));
                         }
 
@@ -513,6 +555,15 @@ macro_rules! mat_mul_mat_impl(
             res
         }
     }
+
+    impl<N: Copy + BaseNum> MulAssign<$t<N>> for $t<N> {
+        #[inline]
+        fn mul_assign(&mut self, right: $t<N>) {
+            // NOTE: there is probably not any useful optimization to perform here compaired to the
+            // version without assignment..
+            *self = *self * right
+        }
+    }
   )
 );
 
@@ -537,6 +588,15 @@ macro_rules! vec_mul_mat_impl(
             res
         }
     }
+
+    impl<N: Copy + BaseNum> MulAssign<$t<N>> for $v<N> {
+        #[inline]
+        fn mul_assign(&mut self, right: $t<N>) {
+            // NOTE: there is probably not any useful optimization to perform here compaired to the
+            // version without assignment..
+            *self = *self * right
+        }
+    }
   )
 );
 
@@ -549,8 +609,8 @@ macro_rules! mat_mul_vec_impl(
         fn mul(self, right: $v<N>) -> $v<N> {
             let mut res : $v<N> = $zero();
 
-            for i in 0..$dim {
-                for j in 0..$dim {
+            for i in 0 .. $dim {
+                for j in 0 .. $dim {
                     unsafe {
                         let val = res.at_fast(i) + self.at_fast((i, j)) * right.at_fast(j);
                         res.set_fast(i, val)
diff --git a/src/structs/pnt.rs b/src/structs/pnt.rs
index 2731cf33..2b5a9b74 100644
--- a/src/structs/pnt.rs
+++ b/src/structs/pnt.rs
@@ -4,7 +4,7 @@ use std::mem;
 use std::fmt;
 use std::slice::{Iter, IterMut};
 use std::iter::{Iterator, FromIterator, IntoIterator};
-use std::ops::{Add, Sub, Mul, Div, Neg, Index, IndexMut};
+use std::ops::{Add, Sub, Mul, Div, Neg, AddAssign, SubAssign, MulAssign, DivAssign, Index, IndexMut};
 use rand::{Rand, Rng};
 use num::{Zero, One};
 use traits::operations::{ApproxEq, POrd, POrdering, Axpy};
diff --git a/src/structs/pnt_macros.rs b/src/structs/pnt_macros.rs
index f5df77a9..5d294fc4 100644
--- a/src/structs/pnt_macros.rs
+++ b/src/structs/pnt_macros.rs
@@ -41,6 +41,13 @@ macro_rules! pnt_add_vec_impl(
                 $t::new($(self.$compN + right.$compN),+)
             }
         }
+
+        impl<N: Copy + AddAssign<N>> AddAssign<$tv<N>> for $t<N> {
+            #[inline]
+            fn add_assign(&mut self, right: $tv<N>) {
+                $( self.$compN += right.$compN; )+
+            }
+        }
     )
 );
 
@@ -54,6 +61,13 @@ macro_rules! pnt_sub_vec_impl(
                 $t::new($(self.$compN - right.$compN),+)
             }
         }
+
+        impl<N: Copy + SubAssign<N>> SubAssign<$tv<N>> for $t<N> {
+            #[inline]
+            fn sub_assign(&mut self, right: $tv<N>) {
+                $( self.$compN -= right.$compN; )+
+            }
+        }
     )
 );
 
diff --git a/src/structs/quat.rs b/src/structs/quat.rs
index 74db865d..a4d835d2 100644
--- a/src/structs/quat.rs
+++ b/src/structs/quat.rs
@@ -3,7 +3,7 @@
 use std::fmt;
 use std::mem;
 use std::slice::{Iter, IterMut};
-use std::ops::{Add, Sub, Mul, Div, Neg, Index, IndexMut};
+use std::ops::{Add, Sub, Mul, Div, Neg, AddAssign, SubAssign, MulAssign, DivAssign, Index, IndexMut};
 use std::iter::{FromIterator, IntoIterator};
 use rand::{Rand, Rng};
 use num::{Zero, One};
@@ -146,6 +146,14 @@ impl<N> Mul<Quat<N>> for Quat<N>
     }
 }
 
+impl<N> MulAssign<Quat<N>> for Quat<N>
+    where N: Copy + Mul<N, Output = N> + Sub<N, Output = N> + Add<N, Output = N> {
+    #[inline]
+    fn mul_assign(&mut self, right: Quat<N>) {
+        *self = *self * right;
+    }
+}
+
 impl<N: ApproxEq<N> + BaseFloat> Div<Quat<N>> for Quat<N> {
     type Output = Quat<N>;
 
@@ -155,6 +163,13 @@ impl<N: ApproxEq<N> + BaseFloat> Div<Quat<N>> for Quat<N> {
     }
 }
 
+impl<N: ApproxEq<N> + BaseFloat> DivAssign<Quat<N>> for Quat<N> {
+    #[inline]
+    fn div_assign(&mut self, right: Quat<N>) {
+        *self *= right.inv().expect("Unable to invert the denominator.")
+    }
+}
+
 impl<N: fmt::Display> fmt::Display for Quat<N> {
     fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
         write!(f, "Quaternion {} − ({}, {}, {})", self.w, self.i, self.j, self.k)
@@ -336,6 +351,13 @@ impl<N: BaseFloat + ApproxEq<N>> Div<UnitQuat<N>> for UnitQuat<N> {
     }
 }
 
+impl<N: BaseFloat + ApproxEq<N>> DivAssign<UnitQuat<N>> for UnitQuat<N> {
+    #[inline]
+    fn div_assign(&mut self, other: UnitQuat<N>) {
+        self.q /= other.q
+    }
+}
+
 impl<N: BaseNum> Mul<UnitQuat<N>> for UnitQuat<N> {
     type Output = UnitQuat<N>;
 
@@ -345,6 +367,13 @@ impl<N: BaseNum> Mul<UnitQuat<N>> for UnitQuat<N> {
     }
 }
 
+impl<N: BaseNum> MulAssign<UnitQuat<N>> for UnitQuat<N> {
+    #[inline]
+    fn mul_assign(&mut self, right: UnitQuat<N>) {
+        self.q *= right.q
+    }
+}
+
 impl<N: BaseNum> Mul<Vec3<N>> for UnitQuat<N> {
     type Output = Vec3<N>;
 
@@ -391,6 +420,20 @@ impl<N: BaseNum + Neg<Output = N>> Mul<UnitQuat<N>> for Pnt3<N> {
     }
 }
 
+impl<N: BaseNum + Neg<Output = N>> MulAssign<UnitQuat<N>> for Vec3<N> {
+    #[inline]
+    fn mul_assign(&mut self, right: UnitQuat<N>) {
+        *self = *self * right
+    }
+}
+
+impl<N: BaseNum + Neg<Output = N>> MulAssign<UnitQuat<N>> for Pnt3<N> {
+    #[inline]
+    fn mul_assign(&mut self, right: UnitQuat<N>) {
+        *self = *self * right
+    }
+}
+
 impl<N: BaseFloat> Rotation<Vec3<N>> for UnitQuat<N> {
     #[inline]
     fn rotation(&self) -> Vec3<N> {
diff --git a/src/structs/rot.rs b/src/structs/rot.rs
index 40ef54de..7b87a9d9 100644
--- a/src/structs/rot.rs
+++ b/src/structs/rot.rs
@@ -1,7 +1,7 @@
 //! Rotations matrices.
 
 use std::fmt;
-use std::ops::{Mul, Neg, Index};
+use std::ops::{Mul, Neg, MulAssign, Index};
 use rand::{Rand, Rng};
 use num::{Zero, One};
 use traits::geometry::{Rotate, Rotation, AbsoluteRotate, RotationMatrix, RotationTo, Transform,
diff --git a/src/structs/rot_macros.rs b/src/structs/rot_macros.rs
index 7c6573e5..043baf39 100644
--- a/src/structs/rot_macros.rs
+++ b/src/structs/rot_macros.rs
@@ -145,6 +145,13 @@ macro_rules! rot_mul_rot_impl(
                 $t { submat: self.submat * right.submat }
             }
         }
+
+        impl<N: Copy + BaseNum> MulAssign<$t<N>> for $t<N> {
+            #[inline]
+            fn mul_assign(&mut self, right: $t<N>) {
+                self.submat *= right.submat
+            }
+        }
     )
 );
 
@@ -177,6 +184,13 @@ macro_rules! vec_mul_rot_impl(
                 self * right.submat
             }
         }
+
+        impl<N: Copy + BaseNum> MulAssign<$t<N>> for $tv<N> {
+            #[inline]
+            fn mul_assign(&mut self, right: $t<N>) {
+                *self *= right.submat
+            }
+        }
     )
 );
 
diff --git a/src/structs/sim.rs b/src/structs/sim.rs
index e2f0a183..28345bcf 100644
--- a/src/structs/sim.rs
+++ b/src/structs/sim.rs
@@ -1,5 +1,5 @@
 use std::fmt;
-use std::ops::{Mul, Neg};
+use std::ops::{Mul, Neg, MulAssign};
 
 use rand::{Rand, Rng};
 use num::One;
diff --git a/src/structs/sim_macros.rs b/src/structs/sim_macros.rs
index 65f3e371..b82aaf0f 100644
--- a/src/structs/sim_macros.rs
+++ b/src/structs/sim_macros.rs
@@ -123,6 +123,15 @@ macro_rules! sim_mul_sim_impl(
                     self.scale * right.scale)
             }
         }
+
+        impl<N: BaseFloat> MulAssign<$t<N>> for $t<N> {
+            #[inline]
+            fn mul_assign(&mut self, right: $t<N>) {
+                self.isometry.translation += self.isometry.rotation * (right.isometry.translation * self.scale);
+                self.isometry.rotation    *= right.isometry.rotation;
+                self.scale                *= right.scale;
+            }
+        }
     )
 );
 
@@ -140,6 +149,14 @@ macro_rules! sim_mul_iso_impl(
             }
         }
 
+        impl<N: BaseFloat> MulAssign<$ti<N>> for $t<N> {
+            #[inline]
+            fn mul_assign(&mut self, right: $ti<N>) {
+                self.isometry.translation += self.isometry.rotation * (right.translation * self.scale);
+                self.isometry.rotation    *= right.rotation;
+            }
+        }
+
         impl<N: BaseFloat> Mul<$t<N>> for $ti<N> {
             type Output = $t<N>;
 
@@ -168,6 +185,13 @@ macro_rules! sim_mul_rot_impl(
             }
         }
 
+        impl<N: BaseFloat> MulAssign<$tr<N>> for $t<N> {
+            #[inline]
+            fn mul_assign(&mut self, right: $tr<N>) {
+                self.isometry.rotation *= right;
+            }
+        }
+
         impl<N: BaseFloat> Mul<$t<N>> for $tr<N> {
             type Output = $t<N>;
 
diff --git a/src/structs/spec/mat.rs b/src/structs/spec/mat.rs
index d4d80602..c5a8e4c3 100644
--- a/src/structs/spec/mat.rs
+++ b/src/structs/spec/mat.rs
@@ -1,4 +1,4 @@
-use std::ops::{Add, Mul, Neg};
+use std::ops::{Add, Mul, Neg, MulAssign};
 use structs::vec::{Vec2, Vec3};
 use structs::pnt::{Pnt2, Pnt3};
 use structs::mat::{Mat1, Mat2, Mat3};
@@ -345,3 +345,25 @@ impl<N: Copy + Mul<N, Output = N> + Add<N, Output = N>> Mul<Pnt2<N>> for Mat2<N>
         )
     }
 }
+
+
+macro_rules! impl_mul_assign_from_mul(
+    ($tleft: ident, $tright: ident) => (
+        impl<N: Copy + Mul<N, Output = N> + Add<N, Output = N>> MulAssign<$tright<N>> for $tleft<N> {
+            #[inline(always)]
+            fn mul_assign(&mut self, right: $tright<N>) {
+                // NOTE: there is probably no interesting optimization compared to the not-inplace
+                // operation.
+                *self = *self * right
+            }
+        }
+    )
+);
+
+impl_mul_assign_from_mul!(Mat3, Mat3);
+impl_mul_assign_from_mul!(Mat2, Mat2);
+
+impl_mul_assign_from_mul!(Vec3, Mat3);
+impl_mul_assign_from_mul!(Vec2, Mat2);
+impl_mul_assign_from_mul!(Pnt3, Mat3);
+impl_mul_assign_from_mul!(Pnt2, Mat2);
diff --git a/src/structs/vec.rs b/src/structs/vec.rs
index 37849331..93ec73b2 100644
--- a/src/structs/vec.rs
+++ b/src/structs/vec.rs
@@ -1,6 +1,6 @@
 //! Vectors with dimension known at compile-time.
 
-use std::ops::{Add, Sub, Mul, Div, Neg, Index, IndexMut};
+use std::ops::{Add, Sub, Mul, Div, Neg, AddAssign, SubAssign, MulAssign, DivAssign, Index, IndexMut};
 use std::mem;
 use std::slice::{Iter, IterMut};
 use std::iter::{Iterator, FromIterator, IntoIterator};
diff --git a/src/structs/vec_macros.rs b/src/structs/vec_macros.rs
index 0aa3cb46..58c2161c 100644
--- a/src/structs/vec_macros.rs
+++ b/src/structs/vec_macros.rs
@@ -372,6 +372,13 @@ macro_rules! add_impl(
                 $t::new($(self.$compN + right.$compN),+)
             }
         }
+
+        impl<N: AddAssign<N>> AddAssign<$t<N>> for $t<N> {
+            #[inline]
+            fn add_assign(&mut self, right: $t<N>) {
+                $( self.$compN += right.$compN; )+
+            }
+        }
     )
 );
 
@@ -387,6 +394,13 @@ macro_rules! scalar_add_impl(
             }
         }
 
+        impl<N: Copy + AddAssign<N>> AddAssign<N> for $t<N> {
+            #[inline]
+            fn add_assign(&mut self, right: N) {
+                $( self.$compN += right; )+
+            }
+        }
+
         impl Add<$t<f32>> for f32 {
             type Output = $t<f32>;
 
@@ -417,6 +431,13 @@ macro_rules! sub_impl(
                 $t::new($(self.$compN - right.$compN),+)
             }
         }
+
+        impl<N: SubAssign<N>> SubAssign<$t<N>> for $t<N> {
+            #[inline]
+            fn sub_assign(&mut self, right: $t<N>) {
+                $( self.$compN -= right.$compN; )+
+            }
+        }
     )
 );
 
@@ -431,6 +452,13 @@ macro_rules! scalar_sub_impl(
             }
         }
 
+        impl<N: Copy + SubAssign<N>> SubAssign<N> for $t<N> {
+            #[inline]
+            fn sub_assign(&mut self, right: N) {
+                $( self.$compN -= right; )+
+            }
+        }
+
         impl Sub<$t<f32>> for f32 {
             type Output = $t<f32>;
 
@@ -460,6 +488,13 @@ macro_rules! mul_impl(
                 $t::new($(self.$compN * right.$compN),+)
             }
         }
+
+        impl<N: MulAssign<N>> MulAssign<$t<N>> for $t<N> {
+            #[inline]
+            fn mul_assign(&mut self, right: $t<N>) {
+                $( self.$compN *= right.$compN; )+
+            }
+        }
     )
 );
 
@@ -474,6 +509,13 @@ macro_rules! scalar_mul_impl(
             }
         }
 
+        impl<N: Copy + MulAssign<N>> MulAssign<N> for $t<N> {
+            #[inline]
+            fn mul_assign(&mut self, right: N) {
+                $( self.$compN *= right; )+
+            }
+        }
+
         impl Mul<$t<f32>> for f32 {
             type Output = $t<f32>;
 
@@ -504,6 +546,13 @@ macro_rules! div_impl(
                 $t::new($(self.$compN / right.$compN),+)
             }
         }
+
+        impl<N: DivAssign<N>> DivAssign<$t<N>> for $t<N> {
+            #[inline]
+            fn div_assign(&mut self, right: $t<N>) {
+                $( self.$compN /= right.$compN; )+
+            }
+        }
     )
 );
 
@@ -517,6 +566,13 @@ macro_rules! scalar_div_impl(
                 $t::new($(self.$compN / right),+)
             }
         }
+
+        impl<N: Copy + DivAssign<N>> DivAssign<N> for $t<N> {
+            #[inline]
+            fn div_assign(&mut self, right: N) {
+                $( self.$compN /= right; )+
+            }
+        }
     )
 );
 
diff --git a/src/structs/vecn.rs b/src/structs/vecn.rs
index 05689a79..bb896d22 100644
--- a/src/structs/vecn.rs
+++ b/src/structs/vecn.rs
@@ -1,6 +1,6 @@
 use std::slice::{Iter, IterMut};
 use std::iter::{FromIterator, IntoIterator};
-use std::ops::{Add, Sub, Mul, Div, Neg, Index, IndexMut};
+use std::ops::{Add, Sub, Mul, Div, Neg, AddAssign, SubAssign, MulAssign, DivAssign, Index, IndexMut};
 use std::mem;
 use rand::{Rand, Rng};
 use num::{Zero, One};
diff --git a/src/structs/vecn_macros.rs b/src/structs/vecn_macros.rs
index dcdf7a7b..bdc80f3e 100644
--- a/src/structs/vecn_macros.rs
+++ b/src/structs/vecn_macros.rs
@@ -2,6 +2,11 @@
 
 macro_rules! vecn_dvec_common_impl(
     ($vecn: ident $(, $param: ident)*) => (
+        /*
+         *
+         * Zero.
+         *
+         */
         impl<N: Zero + Copy + Clone $(, $param: ArrayLength<N>)*> $vecn<N $(, $param)*> {
             /// Tests if all components of the vector are zeroes.
             #[inline]
@@ -10,6 +15,11 @@ macro_rules! vecn_dvec_common_impl(
             }
         }
 
+        /*
+         *
+         * AsRef/AsMut
+         *
+         */
         impl<N $(, $param: ArrayLength<N>)*> AsRef<[N]> for $vecn<N $(, $param)*> {
             #[inline]
             fn as_ref(&self) -> &[N] {
@@ -25,6 +35,11 @@ macro_rules! vecn_dvec_common_impl(
             }
         }
 
+        /*
+         *
+         * Shape.
+         *
+         */
         impl<N $(, $param: ArrayLength<N>)*> Shape<usize> for $vecn<N $(, $param)*> {
             #[inline]
             fn shape(&self) -> usize {
@@ -32,6 +47,11 @@ macro_rules! vecn_dvec_common_impl(
             }
         }
 
+        /*
+         *
+         * Index et. al.
+         *
+         */
         impl<N: Copy $(, $param : ArrayLength<N>)*> Indexable<usize, N> for $vecn<N $(, $param)*> {
             #[inline]
             fn swap(&mut self, i: usize, j: usize) {
@@ -66,6 +86,11 @@ macro_rules! vecn_dvec_common_impl(
             }
         }
 
+        /*
+         *
+         * Iterable et al.
+         *
+         */
         impl<N $(, $param : ArrayLength<N>)*> Iterable<N> for $vecn<N $(, $param)*> {
             #[inline]
             fn iter<'l>(&'l self) -> Iter<'l, N> {
@@ -80,6 +105,11 @@ macro_rules! vecn_dvec_common_impl(
             }
         }
 
+        /*
+         *
+         * Axpy
+         *
+         */
         impl<N: Copy + Add<N, Output = N> + Mul<N, Output = N> $(, $param : ArrayLength<N>)*>
             Axpy<N> for $vecn<N $(, $param)*> {
             fn axpy(&mut self, a: &N, x: &$vecn<N $(, $param)*>) {
@@ -94,6 +124,11 @@ macro_rules! vecn_dvec_common_impl(
             }
         }
 
+        /*
+         *
+         * Mul
+         *
+         */
         impl<N: Copy + Mul<N, Output = N> + Zero $(, $param : ArrayLength<N>)*>
             Mul<$vecn<N $(, $param)*>> for $vecn<N $(, $param)*> {
             type Output = $vecn<N $(, $param)*>;
@@ -112,142 +147,6 @@ macro_rules! vecn_dvec_common_impl(
             }
         }
 
-        impl<N: Copy + Div<N, Output = N> + Zero $(, $param : ArrayLength<N>)*>
-            Div<$vecn<N $(, $param)*>> for $vecn<N $(, $param)*> {
-            type Output = $vecn<N $(, $param)*>;
-
-            #[inline]
-            fn div(self, right: $vecn<N $(, $param)*>) -> $vecn<N $(, $param)*> {
-                assert!(self.len() == right.len());
-
-                let mut res = self;
-
-                for (left, right) in res.as_mut().iter_mut().zip(right.as_ref().iter()) {
-                    *left = *left / *right
-                }
-
-                res
-            }
-        }
-
-        impl<N: Copy + Add<N, Output = N> + Zero $(, $param : ArrayLength<N>)*>
-            Add<$vecn<N $(, $param)*>> for $vecn<N $(, $param)*> {
-            type Output = $vecn<N $(, $param)*>;
-
-            #[inline]
-            fn add(self, right: $vecn<N $(, $param)*>) -> $vecn<N $(, $param)*> {
-                assert!(self.len() == right.len());
-
-                let mut res = self;
-
-                for (left, right) in res.as_mut().iter_mut().zip(right.as_ref().iter()) {
-                    *left = *left + *right
-                }
-
-                res
-            }
-        }
-
-        impl<N: Copy + Sub<N, Output = N> + Zero $(, $param : ArrayLength<N>)*>
-            Sub<$vecn<N $(, $param)*>> for $vecn<N $(, $param)*> {
-            type Output = $vecn<N $(, $param)*>;
-
-            #[inline]
-            fn sub(self, right: $vecn<N $(, $param)*>) -> $vecn<N $(, $param)*> {
-                assert!(self.len() == right.len());
-
-                let mut res = self;
-
-                for (left, right) in res.as_mut().iter_mut().zip(right.as_ref().iter()) {
-                    *left = *left - *right
-                }
-
-                res
-            }
-        }
-
-        impl<N: Neg<Output = N> + Zero + Copy $(, $param : ArrayLength<N>)*> Neg for $vecn<N $(, $param)*> {
-            type Output = $vecn<N $(, $param)*>;
-
-            #[inline]
-            fn neg(mut self) -> $vecn<N $(, $param)*> {
-                for e in self.as_mut().iter_mut() {
-                    *e = -*e;
-                }
-
-                self
-            }
-        }
-
-        impl<N: BaseNum $(, $param : ArrayLength<N>)*> Dot<N> for $vecn<N $(, $param)*> {
-            #[inline]
-            fn dot(&self, other: &$vecn<N $(, $param)*>) -> N {
-                assert!(self.len() == other.len());
-                let mut res: N = ::zero();
-                for i in 0 .. self.len() {
-                    res = res + unsafe { self.unsafe_at(i) * other.unsafe_at(i) };
-                }
-                res
-            }
-        }
-
-        impl<N: BaseFloat $(, $param : ArrayLength<N>)*> Norm<N> for $vecn<N $(, $param)*> {
-            #[inline]
-            fn sqnorm(&self) -> N {
-                Dot::dot(self, self)
-            }
-
-            #[inline]
-            fn normalize(&self) -> $vecn<N $(, $param)*> {
-                let mut res : $vecn<N $(, $param)*> = self.clone();
-                let _ = res.normalize_mut();
-                res
-            }
-
-            #[inline]
-            fn normalize_mut(&mut self) -> N {
-                let l = Norm::norm(self);
-
-                for n in self.as_mut().iter_mut() {
-                    *n = *n / l;
-                }
-
-                l
-            }
-        }
-
-        impl<N: BaseFloat + Cast<f64> $(, $param : ArrayLength<N>)*> Mean<N> for $vecn<N $(, $param)*> {
-            #[inline]
-            fn mean(&self) -> N {
-                let normalizer = ::cast(1.0f64 / self.len() as f64);
-                self.iter().fold(::zero(), |acc, x| acc + *x * normalizer)
-            }
-        }
-
-        impl<N: ApproxEq<N> $(, $param : ArrayLength<N>)*> ApproxEq<N> for $vecn<N $(, $param)*> {
-            #[inline]
-            fn approx_epsilon(_: Option<$vecn<N $(, $param)*>>) -> N {
-                ApproxEq::approx_epsilon(None::<N>)
-            }
-
-            #[inline]
-            fn approx_ulps(_: Option<$vecn<N $(, $param)*>>) -> u32 {
-                ApproxEq::approx_ulps(None::<N>)
-            }
-
-            #[inline]
-            fn approx_eq_eps(&self, other: &$vecn<N $(, $param)*>, epsilon: &N) -> bool {
-                let mut zip = self.as_ref().iter().zip(other.as_ref().iter());
-                zip.all(|(a, b)| ApproxEq::approx_eq_eps(a, b, epsilon))
-            }
-
-            #[inline]
-            fn approx_eq_ulps(&self, other: &$vecn<N $(, $param)*>, ulps: u32) -> bool {
-                let mut zip = self.as_ref().iter().zip(other.as_ref().iter());
-                zip.all(|(a, b)| ApproxEq::approx_eq_ulps(a, b, ulps))
-            }
-        }
-
         impl<N: Copy + Mul<N, Output = N> + Zero $(, $param : ArrayLength<N>)*>
             Mul<N> for $vecn<N $(, $param)*> {
             type Output = $vecn<N $(, $param)*>;
@@ -264,6 +163,28 @@ macro_rules! vecn_dvec_common_impl(
             }
         }
 
+        impl<N> MulAssign<$vecn<N $(, $param)*>> for $vecn<N $(, $param)*>
+            where N: Copy + MulAssign<N> + Zero $(, $param : ArrayLength<N>)* {
+            #[inline]
+            fn mul_assign(&mut self, right: $vecn<N $(, $param)*>) {
+                assert!(self.len() == right.len());
+
+                for (left, right) in self.as_mut().iter_mut().zip(right.as_ref().iter()) {
+                    *left *= *right
+                }
+            }
+        }
+
+        impl<N> MulAssign<N> for $vecn<N $(, $param)*>
+            where N: Copy + MulAssign<N> + Zero $(, $param : ArrayLength<N>)* {
+            #[inline]
+            fn mul_assign(&mut self, right: N) {
+                for e in self.as_mut().iter_mut() {
+                    *e *= right
+                }
+            }
+        }
+
         impl<$($param : ArrayLength<N>),*> Mul<$vecn<f32 $(, $param)*>> for f32 {
             type Output = $vecn<f32 $(, $param)*>;
 
@@ -294,6 +215,29 @@ macro_rules! vecn_dvec_common_impl(
             }
         }
 
+        /*
+         *
+         * Div.
+         *
+         */
+        impl<N: Copy + Div<N, Output = N> + Zero $(, $param : ArrayLength<N>)*>
+            Div<$vecn<N $(, $param)*>> for $vecn<N $(, $param)*> {
+            type Output = $vecn<N $(, $param)*>;
+
+            #[inline]
+            fn div(self, right: $vecn<N $(, $param)*>) -> $vecn<N $(, $param)*> {
+                assert!(self.len() == right.len());
+
+                let mut res = self;
+
+                for (left, right) in res.as_mut().iter_mut().zip(right.as_ref().iter()) {
+                    *left = *left / *right
+                }
+
+                res
+            }
+        }
+
         impl<N: Copy + Div<N, Output = N> + Zero $(, $param : ArrayLength<N>)*> Div<N> for $vecn<N $(, $param)*> {
             type Output = $vecn<N $(, $param)*>;
 
@@ -309,6 +253,51 @@ macro_rules! vecn_dvec_common_impl(
             }
         }
 
+        impl<N> DivAssign<$vecn<N $(, $param)*>> for $vecn<N $(, $param)*>
+            where N: Copy + DivAssign<N> + Zero $(, $param : ArrayLength<N>)* {
+            #[inline]
+            fn div_assign(&mut self, right: $vecn<N $(, $param)*>) {
+                assert!(self.len() == right.len());
+
+                for (left, right) in self.as_mut().iter_mut().zip(right.as_ref().iter()) {
+                    *left /= *right
+                }
+            }
+        }
+
+        impl<N> DivAssign<N> for $vecn<N $(, $param)*>
+            where N: Copy + DivAssign<N> + Zero $(, $param : ArrayLength<N>)* {
+            #[inline]
+            fn div_assign(&mut self, right: N) {
+                for e in self.as_mut().iter_mut() {
+                    *e /= right
+                }
+            }
+        }
+
+        /*
+         *
+         * Add.
+         *
+         */
+        impl<N: Copy + Add<N, Output = N> + Zero $(, $param : ArrayLength<N>)*>
+            Add<$vecn<N $(, $param)*>> for $vecn<N $(, $param)*> {
+            type Output = $vecn<N $(, $param)*>;
+
+            #[inline]
+            fn add(self, right: $vecn<N $(, $param)*>) -> $vecn<N $(, $param)*> {
+                assert!(self.len() == right.len());
+
+                let mut res = self;
+
+                for (left, right) in res.as_mut().iter_mut().zip(right.as_ref().iter()) {
+                    *left = *left + *right
+                }
+
+                res
+            }
+        }
+
         impl<N: Copy + Add<N, Output = N> + Zero $(, $param : ArrayLength<N>)*> Add<N> for $vecn<N $(, $param)*> {
             type Output = $vecn<N $(, $param)*>;
 
@@ -324,6 +313,28 @@ macro_rules! vecn_dvec_common_impl(
             }
         }
 
+        impl<N> AddAssign<$vecn<N $(, $param)*>> for $vecn<N $(, $param)*>
+            where N: Copy + AddAssign<N> + Zero $(, $param : ArrayLength<N>)* {
+            #[inline]
+            fn add_assign(&mut self, right: $vecn<N $(, $param)*>) {
+                assert!(self.len() == right.len());
+
+                for (left, right) in self.as_mut().iter_mut().zip(right.as_ref().iter()) {
+                    *left += *right
+                }
+            }
+        }
+
+        impl<N> AddAssign<N> for $vecn<N $(, $param)*>
+            where N: Copy + AddAssign<N> + Zero $(, $param : ArrayLength<N>)* {
+            #[inline]
+            fn add_assign(&mut self, right: N) {
+                for e in self.as_mut().iter_mut() {
+                    *e += right
+                }
+            }
+        }
+
         impl<$($param : ArrayLength<f32>),*> Add<$vecn<f32 $(, $param)*>> for f32 {
             type Output = $vecn<f32 $(, $param)*>;
 
@@ -354,6 +365,29 @@ macro_rules! vecn_dvec_common_impl(
             }
         }
 
+        /*
+         *
+         * Sub.
+         *
+         */
+        impl<N: Copy + Sub<N, Output = N> + Zero $(, $param : ArrayLength<N>)*>
+            Sub<$vecn<N $(, $param)*>> for $vecn<N $(, $param)*> {
+            type Output = $vecn<N $(, $param)*>;
+
+            #[inline]
+            fn sub(self, right: $vecn<N $(, $param)*>) -> $vecn<N $(, $param)*> {
+                assert!(self.len() == right.len());
+
+                let mut res = self;
+
+                for (left, right) in res.as_mut().iter_mut().zip(right.as_ref().iter()) {
+                    *left = *left - *right
+                }
+
+                res
+            }
+        }
+
         impl<N: Copy + Sub<N, Output = N> + Zero $(, $param : ArrayLength<N>)*> Sub<N> for $vecn<N $(, $param)*> {
             type Output = $vecn<N $(, $param)*>;
 
@@ -369,6 +403,28 @@ macro_rules! vecn_dvec_common_impl(
             }
         }
 
+        impl<N> SubAssign<$vecn<N $(, $param)*>> for $vecn<N $(, $param)*>
+            where N: Copy + SubAssign<N> + Zero $(, $param : ArrayLength<N>)* {
+            #[inline]
+            fn sub_assign(&mut self, right: $vecn<N $(, $param)*>) {
+                assert!(self.len() == right.len());
+
+                for (left, right) in self.as_mut().iter_mut().zip(right.as_ref().iter()) {
+                    *left -= *right
+                }
+            }
+        }
+
+        impl<N> SubAssign<N> for $vecn<N $(, $param)*>
+            where N: Copy + SubAssign<N> + Zero $(, $param : ArrayLength<N>)* {
+            #[inline]
+            fn sub_assign(&mut self, right: N) {
+                for e in self.as_mut().iter_mut() {
+                    *e -= right
+                }
+            }
+        }
+
         impl<$($param : ArrayLength<f32>),*> Sub<$vecn<f32 $(, $param)*>> for f32 {
             type Output = $vecn<f32 $(, $param)*>;
 
@@ -398,5 +454,112 @@ macro_rules! vecn_dvec_common_impl(
                 res
             }
         }
+
+        /*
+         *
+         * Neg.
+         *
+         */
+        impl<N: Neg<Output = N> + Zero + Copy $(, $param : ArrayLength<N>)*> Neg for $vecn<N $(, $param)*> {
+            type Output = $vecn<N $(, $param)*>;
+
+            #[inline]
+            fn neg(mut self) -> $vecn<N $(, $param)*> {
+                for e in self.as_mut().iter_mut() {
+                    *e = -*e;
+                }
+
+                self
+            }
+        }
+
+        /*
+         *
+         * Dot.
+         *
+         */
+        impl<N: BaseNum $(, $param : ArrayLength<N>)*> Dot<N> for $vecn<N $(, $param)*> {
+            #[inline]
+            fn dot(&self, other: &$vecn<N $(, $param)*>) -> N {
+                assert!(self.len() == other.len());
+                let mut res: N = ::zero();
+                for i in 0 .. self.len() {
+                    res = res + unsafe { self.unsafe_at(i) * other.unsafe_at(i) };
+                }
+                res
+            }
+        }
+
+        /*
+         *
+         * Norm.
+         *
+         */
+        impl<N: BaseFloat $(, $param : ArrayLength<N>)*> Norm<N> for $vecn<N $(, $param)*> {
+            #[inline]
+            fn sqnorm(&self) -> N {
+                Dot::dot(self, self)
+            }
+
+            #[inline]
+            fn normalize(&self) -> $vecn<N $(, $param)*> {
+                let mut res : $vecn<N $(, $param)*> = self.clone();
+                let _ = res.normalize_mut();
+                res
+            }
+
+            #[inline]
+            fn normalize_mut(&mut self) -> N {
+                let l = Norm::norm(self);
+
+                for n in self.as_mut().iter_mut() {
+                    *n = *n / l;
+                }
+
+                l
+            }
+        }
+
+        /*
+         *
+         * Mean.
+         *
+         */
+        impl<N: BaseFloat + Cast<f64> $(, $param : ArrayLength<N>)*> Mean<N> for $vecn<N $(, $param)*> {
+            #[inline]
+            fn mean(&self) -> N {
+                let normalizer = ::cast(1.0f64 / self.len() as f64);
+                self.iter().fold(::zero(), |acc, x| acc + *x * normalizer)
+            }
+        }
+
+        /*
+         *
+         * ApproxEq
+         *
+         */
+        impl<N: ApproxEq<N> $(, $param : ArrayLength<N>)*> ApproxEq<N> for $vecn<N $(, $param)*> {
+            #[inline]
+            fn approx_epsilon(_: Option<$vecn<N $(, $param)*>>) -> N {
+                ApproxEq::approx_epsilon(None::<N>)
+            }
+
+            #[inline]
+            fn approx_ulps(_: Option<$vecn<N $(, $param)*>>) -> u32 {
+                ApproxEq::approx_ulps(None::<N>)
+            }
+
+            #[inline]
+            fn approx_eq_eps(&self, other: &$vecn<N $(, $param)*>, epsilon: &N) -> bool {
+                let mut zip = self.as_ref().iter().zip(other.as_ref().iter());
+                zip.all(|(a, b)| ApproxEq::approx_eq_eps(a, b, epsilon))
+            }
+
+            #[inline]
+            fn approx_eq_ulps(&self, other: &$vecn<N $(, $param)*>, ulps: u32) -> bool {
+                let mut zip = self.as_ref().iter().zip(other.as_ref().iter());
+                zip.all(|(a, b)| ApproxEq::approx_eq_ulps(a, b, ulps))
+            }
+        }
     )
 );
diff --git a/src/traits/structure.rs b/src/traits/structure.rs
index 468c0258..479a1a92 100644
--- a/src/traits/structure.rs
+++ b/src/traits/structure.rs
@@ -2,7 +2,9 @@
 
 use std::{f32, f64, i8, i16, i32, i64, u8, u16, u32, u64, isize, usize};
 use std::slice::{Iter, IterMut};
-use std::ops::{Add, Sub, Mul, Div, Rem, Index, IndexMut, Neg};
+use std::ops::{Add, Sub, Mul, Div, Rem,
+               AddAssign, SubAssign, MulAssign, DivAssign, RemAssign,
+               Index, IndexMut, Neg};
 use num::{Float, Zero, One};
 use traits::operations::{Axpy, Transpose, Inv, Absolute};
 use traits::geometry::{Dot, Norm, Orig};
@@ -11,8 +13,11 @@ use traits::geometry::{Dot, Norm, Orig};
 pub trait BaseNum: Copy + Zero + One +
                    Add<Self, Output = Self> + Sub<Self, Output = Self> +
                    Mul<Self, Output = Self> + Div<Self, Output = Self> +
-                   Rem<Self, Output = Self> + PartialEq +
-                   Absolute<Self> + Axpy<Self> {
+                   Rem<Self, Output = Self> +
+                   AddAssign<Self> + SubAssign<Self> +
+                   MulAssign<Self> + DivAssign<Self> +
+                   RemAssign<Self> +
+                   PartialEq + Absolute<Self> + Axpy<Self> {
 }
 
 /// Basic floating-point number numeric trait.
@@ -221,10 +226,19 @@ pub trait IterableMut<N> {
  * Vec related traits.
  */
 /// Trait grouping most common operations on vectors.
-pub trait NumVec<N>: Dim +
-                     Sub<Self, Output = Self> + Add<Self, Output = Self> +
+pub trait NumVec<N>: Add<Self, Output = Self> + Sub<Self, Output = Self> +
+                     Mul<Self, Output = Self> + Div<Self, Output = Self> +
+
+                     Add<N, Output = Self> + Sub<N, Output = Self> + 
                      Mul<N, Output = Self> + Div<N, Output = Self> + 
-                     Index<usize, Output = N> +
+
+                     AddAssign<Self> + SubAssign<Self> +
+                     MulAssign<Self> + DivAssign<Self> + 
+
+                     AddAssign<N> + SubAssign<N> + 
+                     MulAssign<N> + DivAssign<N> + 
+
+                     Dim + Index<usize, Output = N> +
                      Zero + PartialEq + Dot<N> + Axpy<N> {
 }
 
@@ -263,9 +277,13 @@ pub trait NumPnt<N>:
           PartialEq +
           Axpy<N> +
           Sub<Self, Output = <Self as PntAsVec>::Vec> +
-          Mul<N, Output = Self> +
-          Div<N, Output = Self> +
+
+          Mul<N, Output = Self> + Div<N, Output = Self> +
           Add<<Self as PntAsVec>::Vec, Output = Self> +
+
+          MulAssign<N> + DivAssign<N> +
+          AddAssign<<Self as PntAsVec>::Vec> +
+
           Index<usize, Output = N> { // FIXME: + Sub<V, Self>
 }
 
diff --git a/tests/op_assign.rs b/tests/op_assign.rs
new file mode 100644
index 00000000..e4d62657
--- /dev/null
+++ b/tests/op_assign.rs
@@ -0,0 +1,76 @@
+extern crate nalgebra as na;
+extern crate rand;
+
+use std::ops::{Mul, Div, Add, Sub, MulAssign, DivAssign, AddAssign, SubAssign};
+use rand::random;
+use na::{Pnt3, Vec3, Mat3, Rot3, Iso3, Sim3, Quat, UnitQuat};
+
+// NOTE: we test only the 3D version because the others share the same code anyway.
+
+macro_rules! test_op_vs_op_assign(
+    ($name: ident, $t1: ty, $t2: ty, $op: ident, $op_assign: ident) => (
+        #[test]
+        fn $name() {
+            for _ in 0usize .. 10000 {
+                let rand1 = random::<$t1>();
+                let rand2 = random::<$t2>();
+                let mut res = rand1;
+
+                res.$op_assign(rand2);
+
+                assert_eq!(rand1.$op(rand2), res)
+            }
+        }
+    )
+);
+
+// Multiplication.
+test_op_vs_op_assign!(test_vec3_f32_mul_assign, Vec3<f32>, f32, mul, mul_assign);
+test_op_vs_op_assign!(test_mat3_f32_mul_assign, Mat3<f32>, f32, mul, mul_assign);
+test_op_vs_op_assign!(test_quat_f32_mul_assign, Quat<f32>, f32, mul, mul_assign);
+
+test_op_vs_op_assign!(test_vec3_vec3_mul_assign, Vec3<f32>, Vec3<f32>, mul, mul_assign);
+test_op_vs_op_assign!(test_quat_quat_mul_assign, Quat<f32>, Quat<f32>, mul, mul_assign);
+test_op_vs_op_assign!(test_unit_quat_unit_quat_mul_assign, UnitQuat<f32>, UnitQuat<f32>, mul, mul_assign);
+
+test_op_vs_op_assign!(test_vec3_unit_quat_mul_assign, Vec3<f32>, UnitQuat<f32>, mul, mul_assign);
+test_op_vs_op_assign!(test_pnt3_unit_quat_mul_assign, Pnt3<f32>, UnitQuat<f32>, mul, mul_assign);
+
+test_op_vs_op_assign!(test_mat3_mat3_mul_assign, Mat3<f32>, Mat3<f32>, mul, mul_assign);
+test_op_vs_op_assign!(test_vec3_mat3_mul_assign, Vec3<f32>, Mat3<f32>, mul, mul_assign);
+test_op_vs_op_assign!(test_pnt3_mat3_mul_assign, Pnt3<f32>, Mat3<f32>, mul, mul_assign);
+
+test_op_vs_op_assign!(test_rot3_rot3_mul_assign, Rot3<f32>, Rot3<f32>, mul, mul_assign);
+test_op_vs_op_assign!(test_vec3_rot3_mul_assign, Vec3<f32>, Rot3<f32>, mul, mul_assign);
+test_op_vs_op_assign!(test_pnt3_rot3_mul_assign, Pnt3<f32>, Rot3<f32>, mul, mul_assign);
+
+test_op_vs_op_assign!(test_iso3_iso3_mul_assign, Iso3<f32>, Iso3<f32>, mul, mul_assign);
+test_op_vs_op_assign!(test_iso3_rot3_mul_assign, Iso3<f32>, Rot3<f32>, mul, mul_assign);
+
+test_op_vs_op_assign!(test_sim3_sim3_mul_assign, Sim3<f32>, Sim3<f32>, mul, mul_assign);
+test_op_vs_op_assign!(test_sim3_iso3_mul_assign, Sim3<f32>, Iso3<f32>, mul, mul_assign);
+test_op_vs_op_assign!(test_sim3_rot3_mul_assign, Sim3<f32>, Rot3<f32>, mul, mul_assign);
+
+// Division.
+test_op_vs_op_assign!(test_vec3_vec3_div_assign, Vec3<f32>, Vec3<f32>, div, div_assign);
+test_op_vs_op_assign!(test_quat_quat_div_assign, Quat<f32>, Quat<f32>, div, div_assign);
+test_op_vs_op_assign!(test_unit_quat_unit_quat_div_assign, UnitQuat<f32>, UnitQuat<f32>, div, div_assign);
+
+test_op_vs_op_assign!(test_vec3_f32_div_assign, Vec3<f32>, f32, div, div_assign);
+test_op_vs_op_assign!(test_mat3_f32_div_assign, Mat3<f32>, f32, div, div_assign);
+
+// Addition.
+test_op_vs_op_assign!(test_vec3_vec3_add_assign, Vec3<f32>, Vec3<f32>, add, add_assign);
+test_op_vs_op_assign!(test_mat3_mat3_add_assign, Mat3<f32>, Mat3<f32>, add, add_assign);
+test_op_vs_op_assign!(test_quat_quat_add_assign, Quat<f32>, Quat<f32>, add, add_assign);
+
+test_op_vs_op_assign!(test_vec3_f32_add_assign, Vec3<f32>, f32, add, add_assign);
+test_op_vs_op_assign!(test_mat3_f32_add_assign, Mat3<f32>, f32, add, add_assign);
+
+// Subtraction.
+test_op_vs_op_assign!(test_vec3_vec3_sub_assign, Vec3<f32>, Vec3<f32>, sub, sub_assign);
+test_op_vs_op_assign!(test_mat3_mat3_sub_assign, Mat3<f32>, Mat3<f32>, sub, sub_assign);
+test_op_vs_op_assign!(test_quat_quat_sub_assign, Quat<f32>, Quat<f32>, sub, sub_assign);
+
+test_op_vs_op_assign!(test_vec3_f32_sub_assign, Vec3<f32>, f32, sub, sub_assign);
+test_op_vs_op_assign!(test_mat3_f32_sub_assign, Mat3<f32>, f32, sub, sub_assign);