diff --git a/src/ndim/dmat.rs b/src/dmat.rs
similarity index 97%
rename from src/ndim/dmat.rs
rename to src/dmat.rs
index fdeba343..87f7cb7e 100644
--- a/src/ndim/dmat.rs
+++ b/src/dmat.rs
@@ -7,7 +7,7 @@ use traits::inv::Inv;
 use traits::division_ring::DivisionRing;
 use traits::transpose::Transpose;
 use traits::rlmul::{RMul, LMul};
-use ndim::dvec::{DVec, zero_vec_with_dim};
+use dvec::{DVec, zero_vec_with_dim};
 
 #[deriving(Eq, ToStr, Clone)]
 pub struct DMat<N>
@@ -133,13 +133,13 @@ LMul<DVec<N>> for DMat<N>
   }
 }
 
-impl<N: Clone + Copy + Eq + DivisionRing>
+impl<N: Copy + Eq + DivisionRing>
 Inv for DMat<N>
 {
   #[inline]
   fn inverse(&self) -> DMat<N>
   {
-    let mut res : DMat<N> = self.clone();
+    let mut res : DMat<N> = copy *self;
 
     res.invert();
 
diff --git a/src/ndim/dvec.rs b/src/dvec.rs
similarity index 100%
rename from src/ndim/dvec.rs
rename to src/dvec.rs
diff --git a/src/mat.rs b/src/mat.rs
index d1e0995d..8286f9e3 100644
--- a/src/mat.rs
+++ b/src/mat.rs
@@ -13,330 +13,7 @@ use traits::transpose::Transpose;
 use traits::rlmul::{RMul, LMul};
 use traits::transformation::Transform;
 
-macro_rules! mat_impl(
-  ($t: ident, $dim: expr) => (
-    impl<N> $t<N>
-    {
-      #[inline]
-      pub fn new(mij: [N, ..$dim * $dim]) -> $t<N>
-      { $t { mij: mij } }
-    }
-  )
-)
-
-macro_rules! one_impl(
-  ($t: ident, [ $($value: ident)|+ ] ) => (
-    impl<N: Copy + One + Zero> One for $t<N>
-    {
-      #[inline]
-      fn one() -> $t<N>
-      {
-        let (_0, _1) = (Zero::zero::<N>(), One::one::<N>());
-        return $t::new( [ $( copy $value, )+ ] )
-      }
-    }
-  )
-)
-
-
-macro_rules! zero_impl(
-  ($t: ident, [ $($value: ident)|+ ] ) => (
-    impl<N: Copy + Zero> Zero for $t<N>
-    {
-      #[inline]
-      fn zero() -> $t<N>
-      {
-        let _0 = Zero::zero();
-        return $t::new( [ $( copy $value, )+ ] )
-      }
-
-     #[inline]
-     fn is_zero(&self) -> bool
-     { self.mij.iter().all(|e| e.is_zero()) }
-    }
-  )
-)
-
-macro_rules! dim_impl(
-  ($t: ident, $dim: expr) => (
-    impl<N> Dim for $t<N>
-    {
-      #[inline]
-      fn dim() -> uint
-      { $dim }
-    }
-  )
-)
-
-macro_rules! mat_indexing_impl(
-  ($t: ident, $dim: expr) => (
-    impl<N: Copy> $t<N>
-    {
-      #[inline]
-      pub fn offset(&self, i: uint, j: uint) -> uint
-      { i * $dim + j }
-    
-      #[inline]
-      pub fn set(&mut self, i: uint, j: uint, t: &N)
-      {
-        self.mij[self.offset(i, j)] = copy *t
-      }
-    
-      #[inline]
-      pub fn at(&self, i: uint, j: uint) -> N
-      {
-        copy self.mij[self.offset(i, j)]
-      }
-    }
-  )
-)
-
-macro_rules! mul_impl(
-  ($t: ident, $dim: expr) => (
-    impl<N: Copy + Ring>
-    Mul<$t<N>, $t<N>> for $t<N>
-    {
-      fn mul(&self, other: &$t<N>) -> $t<N>
-      {
-        let mut res: $t<N> = Zero::zero();
-    
-        for iterate(0u, $dim) |i|
-        {
-          for iterate(0u, $dim) |j|
-          {
-            let mut acc = Zero::zero::<N>();
-    
-            for iterate(0u, $dim) |k|
-            { acc = acc + self.at(i, k) * other.at(k, j); }
-    
-            res.set(i, j, &acc);
-          }
-        }
-    
-        res
-      }
-    }
-  )
-)
-
-macro_rules! rmul_impl(
-  ($t: ident, $v: ident, $dim: expr) => (
-    impl<N: Copy + Ring>
-    RMul<$v<N>> for $t<N>
-    {
-      fn rmul(&self, other: &$v<N>) -> $v<N>
-      {
-        let mut res : $v<N> = Zero::zero();
-    
-        for iterate(0u, $dim) |i|
-        {
-          for iterate(0u, $dim) |j|
-          { res.at[i] = res.at[i] + other.at[j] * self.at(i, j); }
-        }
-    
-        res
-      }
-    }
-  )
-)
-
-macro_rules! lmul_impl(
-  ($t: ident, $v: ident, $dim: expr) => (
-    impl<N: Copy + Ring>
-    LMul<$v<N>> for $t<N>
-    {
-      fn lmul(&self, other: &$v<N>) -> $v<N>
-      {
-    
-        let mut res : $v<N> = Zero::zero();
-    
-        for iterate(0u, $dim) |i|
-        {
-          for iterate(0u, $dim) |j|
-          { res.at[i] = res.at[i] + other.at[j] * self.at(j, i); }
-        }
-    
-        res
-      }
-    }
-  )
-)
-
-macro_rules! transform_impl(
-  ($t: ident, $v: ident) => (
-    impl<N: Copy + DivisionRing + Eq>
-    Transform<$v<N>> for $t<N>
-    {
-      #[inline]
-      fn transform_vec(&self, v: &$v<N>) -> $v<N>
-      { self.rmul(v) }
-    
-      #[inline]
-      fn inv_transform(&self, v: &$v<N>) -> $v<N>
-      { self.inverse().transform_vec(v) }
-    }
-  )
-)
-
-macro_rules! inv_impl(
-  ($t: ident, $dim: expr) => (
-    impl<N: Copy + Eq + DivisionRing>
-    Inv for $t<N>
-    {
-      #[inline]
-      fn inverse(&self) -> $t<N>
-      {
-        let mut res : $t<N> = copy *self;
-    
-        res.invert();
-    
-        res
-      }
-    
-      fn invert(&mut self)
-      {
-        let mut res: $t<N> = One::one();
-        let     _0N: N     = Zero::zero();
-    
-        // inversion using Gauss-Jordan elimination
-        for iterate(0u, $dim) |k|
-        {
-          // search a non-zero value on the k-th column
-          // FIXME: would it be worth it to spend some more time searching for the
-          // max instead?
-    
-          let mut n0 = k; // index of a non-zero entry
-    
-          while (n0 != $dim)
-          {
-            if self.at(n0, k) != _0N
-            { break; }
-    
-            n0 = n0 + 1;
-          }
-    
-          // swap pivot line
-          if n0 != k
-          {
-            for iterate(0u, $dim) |j|
-            {
-              let off_n0_j = self.offset(n0, j);
-              let off_k_j  = self.offset(k, j);
-    
-              swap(self.mij, off_n0_j, off_k_j);
-              swap(res.mij,  off_n0_j, off_k_j);
-            }
-          }
-    
-          let pivot = self.at(k, k);
-    
-          for iterate(k, $dim) |j|
-          {
-            let selfval = &(self.at(k, j) / pivot);
-            self.set(k, j, selfval);
-          }
-    
-          for iterate(0u, $dim) |j|
-          {
-            let resval  = &(res.at(k, j)   / pivot);
-            res.set(k, j, resval);
-          }
-    
-          for iterate(0u, $dim) |l|
-          {
-            if l != k
-            {
-              let normalizer = self.at(l, k);
-    
-              for iterate(k, $dim) |j|
-              {
-                let selfval = &(self.at(l, j) - self.at(k, j) * normalizer);
-                self.set(l, j, selfval);
-              }
-    
-              for iterate(0u, $dim) |j|
-              {
-                let resval  = &(res.at(l, j) - res.at(k, j) * normalizer);
-                res.set(l, j, resval);
-              }
-            }
-          }
-        }
-    
-        *self = res;
-      }
-    }
-  )
-)
-
-macro_rules! transpose_impl(
-  ($t: ident, $dim: expr) => (
-    impl<N: Copy> Transpose for $t<N>
-    {
-      #[inline]
-      fn transposed(&self) -> $t<N>
-      {
-        let mut res = copy *self;
-    
-        res.transpose();
-    
-        res
-      }
-    
-      fn transpose(&mut self)
-      {
-        for iterate(1u, $dim) |i|
-        {
-          for iterate(0u, $dim - 1) |j|
-          {
-            let off_i_j = self.offset(i, j);
-            let off_j_i = self.offset(j, i);
-    
-            swap(self.mij, off_i_j, off_j_i);
-          }
-        }
-      }
-    }
-  )
-)
-
-macro_rules! approx_eq_impl(
-  ($t: ident) => (
-    impl<N: ApproxEq<N>> ApproxEq<N> for $t<N>
-    {
-      #[inline]
-      fn approx_epsilon() -> N
-      { ApproxEq::approx_epsilon::<N, N>() }
-    
-      #[inline]
-      fn approx_eq(&self, other: &$t<N>) -> bool
-      {
-        let mut zip = self.mij.iter().zip(other.mij.iter());
-    
-        do zip.all |(a, b)| { a.approx_eq(b) }
-      }
-    
-      #[inline]
-      fn approx_eq_eps(&self, other: &$t<N>, epsilon: &N) -> bool
-      {
-        let mut zip = self.mij.iter().zip(other.mij.iter());
-    
-        do zip.all |(a, b)| { a.approx_eq_eps(b, epsilon) }
-      }
-    }
-  )
-)
-
-macro_rules! rand_impl(
-  ($t: ident, $param: ident, [ $($elem: ident)|+ ]) => (
-    impl<N: Rand> Rand for $t<N>
-    {
-      #[inline]
-      fn rand<R: Rng>($param: &mut R) -> $t<N>
-      { $t::new([ $( $elem.gen(), )+ ]) }
-    }
-  )
-)
+mod mat_impl;
 
 #[deriving(ToStr)]
 pub struct Mat1<N>
@@ -508,93 +185,3 @@ rand_impl!(Mat6, rng, [
            rng | rng | rng | rng | rng | rng |
            rng | rng | rng | rng | rng | rng
            ])
-
-// some specializations:
-impl<N: Copy + DivisionRing>
-Inv for Mat1<N>
-{
-  #[inline]
-  fn inverse(&self) -> Mat1<N>
-  {
-    let mut res : Mat1<N> = copy *self;
-
-    res.invert();
-
-    res
-  }
-
-  #[inline]
-  fn invert(&mut self)
-  {
-    assert!(!self.mij[0].is_zero());
-
-    self.mij[0] = One::one::<N>() / self.mij[0]
-  }
-}
-
-impl<N: Copy + DivisionRing>
-Inv for Mat2<N>
-{
-  #[inline]
-  fn inverse(&self) -> Mat2<N>
-  {
-    let mut res : Mat2<N> = copy *self;
-
-    res.invert();
-
-    res
-  }
-
-  #[inline]
-  fn invert(&mut self)
-  {
-    let det = self.mij[0 * 2 + 0] * self.mij[1 * 2 + 1] - self.mij[1 * 2 + 0] * self.mij[0 * 2 + 1];
-
-    assert!(!det.is_zero());
-
-    *self = Mat2::new([self.mij[1 * 2 + 1] / det , -self.mij[0 * 2 + 1] / det,
-                           -self.mij[1 * 2 + 0] / det, self.mij[0 * 2 + 0] / det])
-  }
-}
-
-impl<N: Copy + DivisionRing>
-Inv for Mat3<N>
-{
-  #[inline]
-  fn inverse(&self) -> Mat3<N>
-  {
-    let mut res = copy *self;
-
-    res.invert();
-
-    res
-  }
-
-  #[inline]
-  fn invert(&mut self)
-  {
-    let minor_m12_m23 = self.mij[1 * 3 + 1] * self.mij[2 * 3 + 2] - self.mij[2 * 3 + 1] * self.mij[1 * 3 + 2];
-    let minor_m11_m23 = self.mij[1 * 3 + 0] * self.mij[2 * 3 + 2] - self.mij[2 * 3 + 0] * self.mij[1 * 3 + 2];
-    let minor_m11_m22 = self.mij[1 * 3 + 0] * self.mij[2 * 3 + 1] - self.mij[2 * 3 + 0] * self.mij[1 * 3 + 1];
-
-    let det = self.mij[0 * 3 + 0] * minor_m12_m23
-              - self.mij[0 * 3 + 1] * minor_m11_m23
-              + self.mij[0 * 3 + 2] * minor_m11_m22;
-
-    assert!(!det.is_zero());
-
-    *self = Mat3::new( [
-      (minor_m12_m23  / det),
-      ((self.mij[0 * 3 + 2] * self.mij[2 * 3 + 1] - self.mij[2 * 3 + 2] * self.mij[0 * 3 + 1]) / det),
-      ((self.mij[0 * 3 + 1] * self.mij[1 * 3 + 2] - self.mij[1 * 3 + 1] * self.mij[0 * 3 + 2]) / det),
-
-      (-minor_m11_m23 / det),
-      ((self.mij[0 * 3 + 0] * self.mij[2 * 3 + 2] - self.mij[2 * 3 + 0] * self.mij[0 * 3 + 2]) / det),
-      ((self.mij[0 * 3 + 2] * self.mij[1 * 3 + 0] - self.mij[1 * 3 + 2] * self.mij[0 * 3 + 0]) / det),
-
-      (minor_m11_m22  / det),
-      ((self.mij[0 * 3 + 1] * self.mij[2 * 3 + 0] - self.mij[2 * 3 + 1] * self.mij[0 * 3 + 0]) / det),
-      ((self.mij[0 * 3 + 0] * self.mij[1 * 3 + 1] - self.mij[1 * 3 + 0] * self.mij[0 * 3 + 1]) / det)
-    ] )
-  }
-}
diff --git a/src/mat_impl.rs b/src/mat_impl.rs
new file mode 100644
index 00000000..f0bfdde6
--- /dev/null
+++ b/src/mat_impl.rs
@@ -0,0 +1,325 @@
+#[macro_escape];
+
+macro_rules! mat_impl(
+  ($t: ident, $dim: expr) => (
+    impl<N> $t<N>
+    {
+      #[inline]
+      pub fn new(mij: [N, ..$dim * $dim]) -> $t<N>
+      { $t { mij: mij } }
+    }
+  )
+)
+
+macro_rules! one_impl(
+  ($t: ident, [ $($value: ident)|+ ] ) => (
+    impl<N: Copy + One + Zero> One for $t<N>
+    {
+      #[inline]
+      fn one() -> $t<N>
+      {
+        let (_0, _1) = (Zero::zero::<N>(), One::one::<N>());
+        return $t::new( [ $( copy $value, )+ ] )
+      }
+    }
+  )
+)
+
+macro_rules! zero_impl(
+  ($t: ident, [ $($value: ident)|+ ] ) => (
+    impl<N: Copy + Zero> Zero for $t<N>
+    {
+      #[inline]
+      fn zero() -> $t<N>
+      {
+        let _0 = Zero::zero();
+        return $t::new( [ $( copy $value, )+ ] )
+      }
+
+     #[inline]
+     fn is_zero(&self) -> bool
+     { self.mij.iter().all(|e| e.is_zero()) }
+    }
+  )
+)
+
+macro_rules! dim_impl(
+  ($t: ident, $dim: expr) => (
+    impl<N> Dim for $t<N>
+    {
+      #[inline]
+      fn dim() -> uint
+      { $dim }
+    }
+  )
+)
+
+macro_rules! mat_indexing_impl(
+  ($t: ident, $dim: expr) => (
+    impl<N: Copy> $t<N>
+    {
+      #[inline]
+      pub fn offset(&self, i: uint, j: uint) -> uint
+      { i * $dim + j }
+    
+      #[inline]
+      pub fn set(&mut self, i: uint, j: uint, t: &N)
+      {
+        self.mij[self.offset(i, j)] = copy *t
+      }
+    
+      #[inline]
+      pub fn at(&self, i: uint, j: uint) -> N
+      {
+        copy self.mij[self.offset(i, j)]
+      }
+    }
+  )
+)
+
+macro_rules! mul_impl(
+  ($t: ident, $dim: expr) => (
+    impl<N: Copy + Ring>
+    Mul<$t<N>, $t<N>> for $t<N>
+    {
+      fn mul(&self, other: &$t<N>) -> $t<N>
+      {
+        let mut res: $t<N> = Zero::zero();
+    
+        for iterate(0u, $dim) |i|
+        {
+          for iterate(0u, $dim) |j|
+          {
+            let mut acc = Zero::zero::<N>();
+    
+            for iterate(0u, $dim) |k|
+            { acc = acc + self.at(i, k) * other.at(k, j); }
+    
+            res.set(i, j, &acc);
+          }
+        }
+    
+        res
+      }
+    }
+  )
+)
+
+macro_rules! rmul_impl(
+  ($t: ident, $v: ident, $dim: expr) => (
+    impl<N: Copy + Ring>
+    RMul<$v<N>> for $t<N>
+    {
+      fn rmul(&self, other: &$v<N>) -> $v<N>
+      {
+        let mut res : $v<N> = Zero::zero();
+    
+        for iterate(0u, $dim) |i|
+        {
+          for iterate(0u, $dim) |j|
+          { res.at[i] = res.at[i] + other.at[j] * self.at(i, j); }
+        }
+    
+        res
+      }
+    }
+  )
+)
+
+macro_rules! lmul_impl(
+  ($t: ident, $v: ident, $dim: expr) => (
+    impl<N: Copy + Ring>
+    LMul<$v<N>> for $t<N>
+    {
+      fn lmul(&self, other: &$v<N>) -> $v<N>
+      {
+    
+        let mut res : $v<N> = Zero::zero();
+    
+        for iterate(0u, $dim) |i|
+        {
+          for iterate(0u, $dim) |j|
+          { res.at[i] = res.at[i] + other.at[j] * self.at(j, i); }
+        }
+    
+        res
+      }
+    }
+  )
+)
+
+macro_rules! transform_impl(
+  ($t: ident, $v: ident) => (
+    impl<N: Copy + DivisionRing + Eq>
+    Transform<$v<N>> for $t<N>
+    {
+      #[inline]
+      fn transform_vec(&self, v: &$v<N>) -> $v<N>
+      { self.rmul(v) }
+    
+      #[inline]
+      fn inv_transform(&self, v: &$v<N>) -> $v<N>
+      { self.inverse().transform_vec(v) }
+    }
+  )
+)
+
+macro_rules! inv_impl(
+  ($t: ident, $dim: expr) => (
+    impl<N: Copy + Eq + DivisionRing>
+    Inv for $t<N>
+    {
+      #[inline]
+      fn inverse(&self) -> $t<N>
+      {
+        let mut res : $t<N> = copy *self;
+    
+        res.invert();
+    
+        res
+      }
+    
+      fn invert(&mut self)
+      {
+        let mut res: $t<N> = One::one();
+        let     _0N: N     = Zero::zero();
+    
+        // inversion using Gauss-Jordan elimination
+        for iterate(0u, $dim) |k|
+        {
+          // search a non-zero value on the k-th column
+          // FIXME: would it be worth it to spend some more time searching for the
+          // max instead?
+    
+          let mut n0 = k; // index of a non-zero entry
+    
+          while (n0 != $dim)
+          {
+            if self.at(n0, k) != _0N
+            { break; }
+    
+            n0 = n0 + 1;
+          }
+    
+          // swap pivot line
+          if n0 != k
+          {
+            for iterate(0u, $dim) |j|
+            {
+              let off_n0_j = self.offset(n0, j);
+              let off_k_j  = self.offset(k, j);
+    
+              swap(self.mij, off_n0_j, off_k_j);
+              swap(res.mij,  off_n0_j, off_k_j);
+            }
+          }
+    
+          let pivot = self.at(k, k);
+    
+          for iterate(k, $dim) |j|
+          {
+            let selfval = &(self.at(k, j) / pivot);
+            self.set(k, j, selfval);
+          }
+    
+          for iterate(0u, $dim) |j|
+          {
+            let resval  = &(res.at(k, j)   / pivot);
+            res.set(k, j, resval);
+          }
+    
+          for iterate(0u, $dim) |l|
+          {
+            if l != k
+            {
+              let normalizer = self.at(l, k);
+    
+              for iterate(k, $dim) |j|
+              {
+                let selfval = &(self.at(l, j) - self.at(k, j) * normalizer);
+                self.set(l, j, selfval);
+              }
+    
+              for iterate(0u, $dim) |j|
+              {
+                let resval  = &(res.at(l, j) - res.at(k, j) * normalizer);
+                res.set(l, j, resval);
+              }
+            }
+          }
+        }
+    
+        *self = res;
+      }
+    }
+  )
+)
+
+macro_rules! transpose_impl(
+  ($t: ident, $dim: expr) => (
+    impl<N: Copy> Transpose for $t<N>
+    {
+      #[inline]
+      fn transposed(&self) -> $t<N>
+      {
+        let mut res = copy *self;
+    
+        res.transpose();
+    
+        res
+      }
+    
+      fn transpose(&mut self)
+      {
+        for iterate(1u, $dim) |i|
+        {
+          for iterate(0u, $dim - 1) |j|
+          {
+            let off_i_j = self.offset(i, j);
+            let off_j_i = self.offset(j, i);
+    
+            swap(self.mij, off_i_j, off_j_i);
+          }
+        }
+      }
+    }
+  )
+)
+
+macro_rules! approx_eq_impl(
+  ($t: ident) => (
+    impl<N: ApproxEq<N>> ApproxEq<N> for $t<N>
+    {
+      #[inline]
+      fn approx_epsilon() -> N
+      { ApproxEq::approx_epsilon::<N, N>() }
+    
+      #[inline]
+      fn approx_eq(&self, other: &$t<N>) -> bool
+      {
+        let mut zip = self.mij.iter().zip(other.mij.iter());
+    
+        do zip.all |(a, b)| { a.approx_eq(b) }
+      }
+    
+      #[inline]
+      fn approx_eq_eps(&self, other: &$t<N>, epsilon: &N) -> bool
+      {
+        let mut zip = self.mij.iter().zip(other.mij.iter());
+    
+        do zip.all |(a, b)| { a.approx_eq_eps(b, epsilon) }
+      }
+    }
+  )
+)
+
+macro_rules! rand_impl(
+  ($t: ident, $param: ident, [ $($elem: ident)|+ ]) => (
+    impl<N: Rand> Rand for $t<N>
+    {
+      #[inline]
+      fn rand<R: Rng>($param: &mut R) -> $t<N>
+      { $t::new([ $( $elem.gen(), )+ ]) }
+    }
+  )
+)
diff --git a/src/mat_spec.rs b/src/mat_spec.rs
new file mode 100644
index 00000000..3f2613c0
--- /dev/null
+++ b/src/mat_spec.rs
@@ -0,0 +1,94 @@
+use std::num::{Zero, One};
+use mat::{Mat1, Mat2, Mat3};
+use traits::division_ring::DivisionRing;
+use traits::inv::Inv;
+
+// some specializations:
+impl<N: Copy + DivisionRing>
+Inv for Mat1<N>
+{
+  #[inline]
+  fn inverse(&self) -> Mat1<N>
+  {
+    let mut res : Mat1<N> = copy *self;
+
+    res.invert();
+
+    res
+  }
+
+  #[inline]
+  fn invert(&mut self)
+  {
+    assert!(!self.mij[0].is_zero());
+
+    self.mij[0] = One::one::<N>() / self.mij[0]
+  }
+}
+
+impl<N: Copy + DivisionRing>
+Inv for Mat2<N>
+{
+  #[inline]
+  fn inverse(&self) -> Mat2<N>
+  {
+    let mut res : Mat2<N> = copy *self;
+
+    res.invert();
+
+    res
+  }
+
+  #[inline]
+  fn invert(&mut self)
+  {
+    let det = self.mij[0 * 2 + 0] * self.mij[1 * 2 + 1] - self.mij[1 * 2 + 0] * self.mij[0 * 2 + 1];
+
+    assert!(!det.is_zero());
+
+    *self = Mat2::new([self.mij[1 * 2 + 1] / det , -self.mij[0 * 2 + 1] / det,
+                           -self.mij[1 * 2 + 0] / det, self.mij[0 * 2 + 0] / det])
+  }
+}
+
+impl<N: Copy + DivisionRing>
+Inv for Mat3<N>
+{
+  #[inline]
+  fn inverse(&self) -> Mat3<N>
+  {
+    let mut res = copy *self;
+
+    res.invert();
+
+    res
+  }
+
+  #[inline]
+  fn invert(&mut self)
+  {
+    let minor_m12_m23 = self.mij[1 * 3 + 1] * self.mij[2 * 3 + 2] - self.mij[2 * 3 + 1] * self.mij[1 * 3 + 2];
+    let minor_m11_m23 = self.mij[1 * 3 + 0] * self.mij[2 * 3 + 2] - self.mij[2 * 3 + 0] * self.mij[1 * 3 + 2];
+    let minor_m11_m22 = self.mij[1 * 3 + 0] * self.mij[2 * 3 + 1] - self.mij[2 * 3 + 0] * self.mij[1 * 3 + 1];
+
+    let det = self.mij[0 * 3 + 0] * minor_m12_m23
+              - self.mij[0 * 3 + 1] * minor_m11_m23
+              + self.mij[0 * 3 + 2] * minor_m11_m22;
+
+    assert!(!det.is_zero());
+
+    *self = Mat3::new( [
+      (minor_m12_m23  / det),
+      ((self.mij[0 * 3 + 2] * self.mij[2 * 3 + 1] - self.mij[2 * 3 + 2] * self.mij[0 * 3 + 1]) / det),
+      ((self.mij[0 * 3 + 1] * self.mij[1 * 3 + 2] - self.mij[1 * 3 + 1] * self.mij[0 * 3 + 2]) / det),
+
+      (-minor_m11_m23 / det),
+      ((self.mij[0 * 3 + 0] * self.mij[2 * 3 + 2] - self.mij[2 * 3 + 0] * self.mij[0 * 3 + 2]) / det),
+      ((self.mij[0 * 3 + 2] * self.mij[1 * 3 + 0] - self.mij[1 * 3 + 2] * self.mij[0 * 3 + 0]) / det),
+
+      (minor_m11_m22  / det),
+      ((self.mij[0 * 3 + 1] * self.mij[2 * 3 + 0] - self.mij[2 * 3 + 1] * self.mij[0 * 3 + 0]) / det),
+      ((self.mij[0 * 3 + 0] * self.mij[1 * 3 + 1] - self.mij[1 * 3 + 0] * self.mij[0 * 3 + 1]) / det)
+    ] )
+  }
+}
diff --git a/src/nalgebra.rc b/src/nalgebra.rc
index ffd783b6..afbcc79a 100644
--- a/src/nalgebra.rc
+++ b/src/nalgebra.rc
@@ -17,15 +17,12 @@ extern mod std;
 
 pub mod vec;
 pub mod mat;
+pub mod dmat;
+pub mod dvec;
 
-/// n-dimensional linear algebra (slower).
-pub mod ndim
-{
-  pub mod dmat;
-  pub mod dvec;
-  pub mod nvec;
-  pub mod nmat;
-}
+// specialization for some 1d, 2d and 3d operations
+pub mod mat_spec;
+pub mod vec_spec;
 
 /// Wrappers around raw matrices to restrict their behaviour.
 pub mod adaptors
diff --git a/src/ndim/nmat.rs b/src/ndim/nmat.rs
deleted file mode 100644
index 998fe510..00000000
--- a/src/ndim/nmat.rs
+++ /dev/null
@@ -1,184 +0,0 @@
-use std::uint::iterate;
-use std::num::{One, Zero};
-use std::rand::{Rand, Rng, RngUtil};
-use std::cmp::ApproxEq;
-use traits::dim::Dim;
-use traits::inv::Inv;
-use traits::division_ring::DivisionRing;
-use traits::transpose::Transpose;
-use traits::rlmul::{RMul, LMul};
-use ndim::dmat::{DMat, one_mat_with_dim, zero_mat_with_dim, is_zero_mat};
-use ndim::nvec::NVec;
-
-// D is a phantom type parameter, used only as a dimensional token.
-// Its allows use to encode the vector dimension at the type-level.
-// It can be anything implementing the Dim trait. However, to avoid confusion,
-// using d0, d1, d2, d3 and d4 tokens are prefered.
-#[deriving(Eq, ToStr)]
-pub struct NMat<D, N>
-{ mij: DMat<N> }
-
-impl<D: Dim, N: Copy> NMat<D, N>
-{
-  #[inline]
-  fn offset(i: uint, j: uint) -> uint
-  { i * Dim::dim::<D>() + j }
-
-  #[inline]
-  fn set(&mut self, i: uint, j: uint, t: &N)
-  { self.mij.set(i, j, t) }
-}
-
-impl<D: Dim, N> Dim for NMat<D, N>
-{
-  #[inline]
-  fn dim() -> uint
-  { Dim::dim::<D>() }
-}
-
-impl<D: Dim, N: Copy> Index<(uint, uint), N> for NMat<D, N>
-{
-  #[inline]
-  fn index(&self, &idx: &(uint, uint)) -> N
-  { self.mij[idx] }
-}
-
-impl<D: Dim, N: Copy + One + Zero> One for NMat<D, N>
-{
-  #[inline]
-  fn one() -> NMat<D, N>
-  { NMat { mij: one_mat_with_dim(Dim::dim::<D>()) } }
-}
-
-impl<D: Dim, N: Copy + Zero> Zero for NMat<D, N>
-{
-  #[inline]
-  fn zero() -> NMat<D, N>
-  { NMat { mij: zero_mat_with_dim(Dim::dim::<D>()) } }
-
-  #[inline]
-  fn is_zero(&self) -> bool
-  { is_zero_mat(&self.mij) }
-}
-
-impl<D: Dim, N: Copy + Mul<N, N> + Add<N, N> + Zero>
-Mul<NMat<D, N>, NMat<D, N>> for NMat<D, N>
-{
-  fn mul(&self, other: &NMat<D, N>) -> NMat<D, N>
-  {
-    let     dim = Dim::dim::<D>();
-    let mut res = Zero::zero::<NMat<D, N>>();
-
-    for iterate(0u, dim) |i|
-    {
-      for iterate(0u, dim) |j|
-      {
-        let mut acc: N = Zero::zero();
-
-        for iterate(0u, dim) |k|
-        { acc = acc + self[(i, k)] * other[(k, j)]; }
-
-        res.set(i, j, &acc);
-      }
-    }
-
-    res
-  }
-}
-
-impl<D: Dim, N: Copy + Add<N, N> + Mul<N, N> + Zero>
-RMul<NVec<D, N>> for NMat<D, N>
-{
-  fn rmul(&self, other: &NVec<D, N>) -> NVec<D, N>
-  {
-    let     dim              = Dim::dim::<D>();
-    let mut res : NVec<D, N> = Zero::zero();
-
-    for iterate(0u, dim) |i|
-    {
-      for iterate(0u, dim) |j|
-      { res.at.at[i] = res.at.at[i] + other.at.at[j] * self[(i, j)]; }
-    }
-
-    res
-  }
-}
-
-impl<D: Dim, N: Copy + Add<N, N> + Mul<N, N> + Zero>
-LMul<NVec<D, N>> for NMat<D, N>
-{
-  fn lmul(&self, other: &NVec<D, N>) -> NVec<D, N>
-  {
-    let     dim              = Dim::dim::<D>();
-    let mut res : NVec<D, N> = Zero::zero();
-
-    for iterate(0u, dim) |i|
-    {
-      for iterate(0u, dim) |j|
-      { res.at.at[i] = res.at.at[i] + other.at.at[j] * self[(j, i)]; }
-    }
-
-    res
-  }
-}
-
-impl<D: Dim, N: Clone + Copy + Eq + DivisionRing>
-Inv for NMat<D, N>
-{
-  #[inline]
-  fn inverse(&self) -> NMat<D, N>
-  { NMat { mij: self.mij.inverse() } }
-
-  #[inline]
-  fn invert(&mut self)
-  { self.mij.invert() }
-}
-
-impl<D: Dim, N:Copy> Transpose for NMat<D, N>
-{
-  #[inline]
-  fn transposed(&self) -> NMat<D, N>
-  {
-    let mut res = copy *self;
-
-    res.transpose();
-
-    res
-  }
-
-  #[inline]
-  fn transpose(&mut self)
-  { self.mij.transpose() }
-}
-
-impl<D, N: ApproxEq<N>> ApproxEq<N> for NMat<D, N>
-{
-  #[inline]
-  fn approx_epsilon() -> N
-  { ApproxEq::approx_epsilon::<N, N>() }
-
-  #[inline]
-  fn approx_eq(&self, other: &NMat<D, N>) -> bool
-  { self.mij.approx_eq(&other.mij) }
-
-  #[inline]
-  fn approx_eq_eps(&self, other: &NMat<D, N>, epsilon: &N) -> bool
-  { self.mij.approx_eq_eps(&other.mij, epsilon) }
-}
-
-impl<D: Dim, N: Rand + Zero + Copy> Rand for NMat<D, N>
-{
-  fn rand<R: Rng>(rng: &mut R) -> NMat<D, N>
-  {
-    let     dim = Dim::dim::<D>();
-    let mut res : NMat<D, N> = Zero::zero();
-
-    for iterate(0u, dim) |i|
-    {
-      for iterate(0u, dim) |j|
-      { res.set(i, j, &rng.gen()); }
-    }
-
-    res
-  }
-}
diff --git a/src/ndim/nvec.rs b/src/ndim/nvec.rs
deleted file mode 100644
index 57dbfbe3..00000000
--- a/src/ndim/nvec.rs
+++ /dev/null
@@ -1,276 +0,0 @@
-use std::uint::iterate;
-use std::num::{Zero, Algebraic, Bounded};
-use std::rand::{Rand, Rng, RngUtil};
-use std::vec::{VecIterator, VecMutIterator};
-use std::cmp::ApproxEq;
-use std::iterator::{IteratorUtil, FromIterator};
-use ndim::dvec::{DVec, zero_vec_with_dim, is_zero_vec};
-use traits::iterable::{Iterable, IterableMut};
-use traits::basis::Basis;
-use traits::ring::Ring;
-use traits::division_ring::DivisionRing;
-use traits::dim::Dim;
-use traits::dot::Dot;
-use traits::sub_dot::SubDot;
-use traits::norm::Norm;
-use traits::translation::{Translation, Translatable};
-use traits::scalar_op::{ScalarMul, ScalarDiv, ScalarAdd, ScalarSub};
-
-// D is a phantom parameter, used only as a dimensional token.
-// Its allows use to encode the vector dimension at the type-level.
-// It can be anything implementing the Dim trait. However, to avoid confusion,
-// using d0, d1, d2, d3, ..., d7 (or your own dn) are prefered.
-#[deriving(Eq, Ord, ToStr)]
-pub struct NVec<D, N>
-{ at: DVec<N> }
-
-
-impl<D: Dim, N> Dim for NVec<D, N>
-{
-  #[inline]
-  fn dim() -> uint
-  { Dim::dim::<D>() }
-}
-
-impl<D, N: Clone> Clone for NVec<D, N>
-{
-  #[inline]
-  fn clone(&self) -> NVec<D, N>
-  { NVec{ at: self.at.clone() } }
-}
-
-impl<D, N> Iterable<N> for NVec<D, N>
-{
-  fn iter<'l>(&'l self) -> VecIterator<'l, N>
-  { self.at.iter() }
-}
-
-impl<D, N> IterableMut<N> for NVec<D, N>
-{
-  fn mut_iter<'l>(&'l mut self) -> VecMutIterator<'l, N>
-  { self.at.mut_iter() }
-}
-
-impl<D: Dim, N: Zero + Copy, Iter: Iterator<N>> FromIterator<N, Iter> for NVec<D, N>
-{
-  fn from_iterator(param: &mut Iter) -> NVec<D, N>
-  {
-    let mut res: NVec<D, N> = Zero::zero();
-    let mut i = 0;
-
-    for param.advance |e|
-    {
-      res.at.at[i] = e;
-      i = i + 1;
-    }
-
-    res
-  }
-}
-
-impl<D, N: Copy + Add<N,N>> Add<NVec<D, N>, NVec<D, N>> for NVec<D, N>
-{
-  #[inline]
-  fn add(&self, other: &NVec<D, N>) -> NVec<D, N>
-  { NVec { at: self.at + other.at } }
-}
-
-impl<D, N: Copy + Sub<N,N>> Sub<NVec<D, N>, NVec<D, N>> for NVec<D, N>
-{
-  #[inline]
-  fn sub(&self, other: &NVec<D, N>) -> NVec<D, N>
-  { NVec { at: self.at - other.at } }
-}
-
-impl<D, N: Neg<N>> Neg<NVec<D, N>> for NVec<D, N>
-{
-  #[inline]
-  fn neg(&self) -> NVec<D, N>
-  { NVec { at: -self.at } }
-}
-
-impl<D: Dim, N: Ring>
-Dot<N> for NVec<D, N>
-{
-  #[inline]
-  fn dot(&self, other: &NVec<D, N>) -> N
-  { self.at.dot(&other.at) } 
-}
-
-impl<D: Dim, N: Ring> SubDot<N> for NVec<D, N>
-{
-  #[inline]
-  fn sub_dot(&self, a: &NVec<D, N>, b: &NVec<D, N>) -> N
-  { self.at.sub_dot(&a.at, &b.at) } 
-}
-
-impl<D: Dim, N: Mul<N, N>>
-ScalarMul<N> for NVec<D, N>
-{
-  #[inline]
-  fn scalar_mul(&self, s: &N) -> NVec<D, N>
-  { NVec { at: self.at.scalar_mul(s) } }
-
-  #[inline]
-  fn scalar_mul_inplace(&mut self, s: &N)
-  { self.at.scalar_mul_inplace(s) }
-}
-
-
-impl<D: Dim, N: Div<N, N>>
-ScalarDiv<N> for NVec<D, N>
-{
-  #[inline]
-  fn scalar_div(&self, s: &N) -> NVec<D, N>
-  { NVec { at: self.at.scalar_div(s) } }
-
-  #[inline]
-  fn scalar_div_inplace(&mut self, s: &N)
-  { self.at.scalar_div_inplace(s) }
-}
-
-impl<D: Dim, N: Add<N, N>>
-ScalarAdd<N> for NVec<D, N>
-{
-  #[inline]
-  fn scalar_add(&self, s: &N) -> NVec<D, N>
-  { NVec { at: self.at.scalar_add(s) } }
-
-  #[inline]
-  fn scalar_add_inplace(&mut self, s: &N)
-  { self.at.scalar_add_inplace(s) }
-}
-
-impl<D: Dim, N: Sub<N, N>>
-ScalarSub<N> for NVec<D, N>
-{
-  #[inline]
-  fn scalar_sub(&self, s: &N) -> NVec<D, N>
-  { NVec { at: self.at.scalar_sub(s) } }
-
-  #[inline]
-  fn scalar_sub_inplace(&mut self, s: &N)
-  { self.at.scalar_sub_inplace(s) }
-}
-
-impl<D: Dim, N: Clone + Copy + Add<N, N> + Neg<N>> Translation<NVec<D, N>> for NVec<D, N>
-{
-  #[inline]
-  fn translation(&self) -> NVec<D, N>
-  { self.clone() }
-
-  #[inline]
-  fn inv_translation(&self) -> NVec<D, N>
-  { -self.clone() }
-
-
-  #[inline]
-  fn translate_by(&mut self, t: &NVec<D, N>)
-  { *self = *self + *t; }
-}
-
-impl<D: Dim, N: Add<N, N> + Copy> Translatable<NVec<D, N>, NVec<D, N>> for NVec<D, N>
-{
-  #[inline]
-  fn translated(&self, t: &NVec<D, N>) -> NVec<D, N>
-  { self + *t }
-}
-
-impl<D: Dim, N: Copy + DivisionRing + Algebraic + Clone>
-Norm<N> for NVec<D, N>
-{
-  #[inline]
-  fn sqnorm(&self) -> N
-  { self.at.sqnorm() }
-
-  #[inline]
-  fn norm(&self) -> N
-  { self.at.norm() }
-
-  #[inline]
-  fn normalized(&self) -> NVec<D, N>
-  {
-    let mut res : NVec<D, N> = self.clone();
-
-    res.normalize();
-
-    res
-  }
-
-  #[inline]
-  fn normalize(&mut self) -> N
-  { self.at.normalize() }
-}
-
-impl<D: Dim,
-     N: Copy + DivisionRing + Algebraic + Clone + ApproxEq<N>>
-Basis for NVec<D, N>
-{
-  #[inline]
-  fn canonical_basis() -> ~[NVec<D, N>]
-  { DVec::canonical_basis_with_dim(Dim::dim::<D>()).map(|&e| NVec { at: e }) }
-
-  #[inline]
-  fn orthogonal_subspace_basis(&self) -> ~[NVec<D, N>]
-  { self.at.orthogonal_subspace_basis().map(|&e| NVec { at: e }) }
-}
-
-// FIXME: I dont really know how te generalize the cross product int
-// n-dimensions…
-// impl<N: Copy + Mul<N, N> + Sub<N, N>> Cross<N> for NVec<D, N>
-// {
-//   fn cross(&self, other: &NVec<D, N>) -> N
-//   { self.x * other.y - self.y * other.x }
-// }
-
-impl<D: Dim, N: Copy + Zero> Zero for NVec<D, N>
-{
-  #[inline]
-  fn zero() -> NVec<D, N>
-  { NVec { at: zero_vec_with_dim(Dim::dim::<D>()) } }
-
-  #[inline]
-  fn is_zero(&self) -> bool
-  { is_zero_vec(&self.at) }
-}
-
-impl<D, N: ApproxEq<N>> ApproxEq<N> for NVec<D, N>
-{
-  #[inline]
-  fn approx_epsilon() -> N
-  { ApproxEq::approx_epsilon::<N, N>() }
-
-  #[inline]
-  fn approx_eq(&self, other: &NVec<D, N>) -> bool
-  { self.at.approx_eq(&other.at) }
-
-  #[inline]
-  fn approx_eq_eps(&self, other: &NVec<D, N>, epsilon: &N) -> bool
-  { self.at.approx_eq_eps(&other.at, epsilon) }
-}
-
-impl<D: Dim, N: Rand + Zero + Copy> Rand for NVec<D, N>
-{
-  #[inline]
-  fn rand<R: Rng>(rng: &mut R) -> NVec<D, N>
-  {
-    let     dim = Dim::dim::<D>();
-    let mut res : NVec<D, N> = Zero::zero();
-
-    for iterate(0u, dim) |i|
-    { res.at.at[i] = rng.gen() }
-
-    res
-  }
-}
-
-impl<D: Dim, N: Bounded + Zero + Add<N, N> + Copy> Bounded for NVec<D, N>
-{
-  #[inline]
-  fn max_value() -> NVec<D, N>
-  { Zero::zero::<NVec<D, N>>().scalar_add(&Bounded::max_value()) }
-
-  #[inline]
-  fn min_value() -> NVec<D, N>
-  { Zero::zero::<NVec<D, N>>().scalar_add(&Bounded::min_value()) }
-}
diff --git a/src/tests/vec.rs b/src/tests/vec.rs
index f3602e27..5db1d5a5 100644
--- a/src/tests/vec.rs
+++ b/src/tests/vec.rs
@@ -7,11 +7,7 @@ use std::rand::{random};
 #[test]
 use std::cmp::ApproxEq;
 #[test]
-use vec::{Vec1, Vec2, Vec3};
-#[test]
-use ndim::nvec::NVec;
-#[test]
-use traits::dim::d7;
+use vec::{Vec1, Vec2, Vec3, Vec4, Vec5, Vec6};
 #[test]
 use traits::basis::Basis;
 #[test]
@@ -133,7 +129,7 @@ fn test_cross_vec3()
 
 #[test]
 fn test_commut_dot_nvec()
-{ test_commut_dot_impl!(NVec<d7, f64>); }
+{ test_commut_dot_impl!(Vec6<f64>); }
  
 #[test]
 fn test_commut_dot_vec3()
@@ -160,8 +156,16 @@ fn test_basis_vec3()
 { test_basis_impl!(Vec3<f64>); }
 
 #[test]
-fn test_basis_nvec()
-{ test_basis_impl!(NVec<d7, f64>); }
+fn test_basis_vec4()
+{ test_basis_impl!(Vec4<f64>); }
+
+#[test]
+fn test_basis_vec5()
+{ test_basis_impl!(Vec5<f64>); }
+
+#[test]
+fn test_basis_vec6()
+{ test_basis_impl!(Vec6<f64>); }
 
 #[test]
 fn test_subspace_basis_vec1()
@@ -176,8 +180,16 @@ fn test_subspace_basis_vec3()
 { test_subspace_basis_impl!(Vec3<f64>); }
 
 #[test]
-fn test_subspace_basis_nvec()
-{ test_subspace_basis_impl!(NVec<d7, f64>); }
+fn test_subspace_basis_vec4()
+{ test_subspace_basis_impl!(Vec4<f64>); }
+
+#[test]
+fn test_subspace_basis_vec5()
+{ test_subspace_basis_impl!(Vec5<f64>); }
+
+#[test]
+fn test_subspace_basis_vec6()
+{ test_subspace_basis_impl!(Vec6<f64>); }
 
 #[test]
 fn test_scalar_op_vec1()
@@ -190,10 +202,18 @@ fn test_scalar_op_vec2()
 #[test]
 fn test_scalar_op_vec3()
 { test_scalar_op_impl!(Vec3<f64>, f64); }
-
+ 
 #[test]
-fn test_scalar_op_nvec()
-{ test_scalar_op_impl!(NVec<d7, f64>, f64); }
+fn test_scalar_op_vec4()
+{ test_scalar_op_impl!(Vec4<f64>, f64); }
+ 
+#[test]
+fn test_scalar_op_vec5()
+{ test_scalar_op_impl!(Vec5<f64>, f64); }
+ 
+#[test]
+fn test_scalar_op_vec6()
+{ test_scalar_op_impl!(Vec6<f64>, f64); }
 
 #[test]
 fn test_iterator_vec1()
@@ -206,7 +226,15 @@ fn test_iterator_vec2()
 #[test]
 fn test_iterator_vec3()
 { test_iterator_impl!(Vec3<f64>, f64); }
-
+ 
 #[test]
-fn test_iterator_nvec()
-{ test_iterator_impl!(NVec<d7, f64>, f64); }
+fn test_iterator_vec4()
+{ test_iterator_impl!(Vec4<f64>, f64); }
+ 
+#[test]
+fn test_iterator_vec5()
+{ test_iterator_impl!(Vec5<f64>, f64); }
+ 
+#[test]
+fn test_iterator_vec6()
+{ test_iterator_impl!(Vec6<f64>, f64); }
diff --git a/src/vec.rs b/src/vec.rs
index 686cf409..9ba016b0 100644
--- a/src/vec.rs
+++ b/src/vec.rs
@@ -1,457 +1,22 @@
-use std::num::{abs, Zero, One, Algebraic, Bounded};
+use std::num::{Zero, One, Algebraic, Bounded};
 use std::rand::{Rand, Rng, RngUtil};
 use std::vec::{VecIterator, VecMutIterator};
-use std::iterator::{Iterator, FromIterator};
+use std::iterator::{Iterator, IteratorUtil, FromIterator};
 use std::cmp::ApproxEq;
+use std::uint::iterate;
 use traits::iterable::{Iterable, IterableMut, FromAnyIterator};
 use traits::basis::Basis;
-use traits::cross::Cross;
 use traits::dim::Dim;
 use traits::dot::Dot;
 use traits::sub_dot::SubDot;
 use traits::norm::Norm;
 use traits::translation::{Translation, Translatable};
 use traits::scalar_op::{ScalarMul, ScalarDiv, ScalarAdd, ScalarSub};
-use std::uint::iterate;
-
-use std::iterator::IteratorUtil;
 use traits::ring::Ring;
 use traits::division_ring::DivisionRing;
 
-// FIXME: is there a way to split this file:
-//  − one file for macros
-//  − another file for macro calls and specializations
-// ?
+mod vec_impl;
 
-macro_rules! new_impl(
-  ($t: ident, $dim: expr) => (
-    impl<N> $t<N>
-    {
-      #[inline]
-      pub fn new(at: [N, ..$dim]) -> $t<N>
-      { $t { at: at } }
-    }
-  )
-)
-
-macro_rules! new_repeat_impl(
-  ($t: ident, $param: ident, [ $($elem: ident)|+ ]) => (
-    impl<N: Copy> $t<N>
-    {
-      #[inline]
-      pub fn new_repeat($param: N) -> $t<N>
-      { $t{ at: [ $( copy $elem, )+ ] } }
-    }
-  )
-)
-
-macro_rules! iterable_impl(
-  ($t: ident) => (
-    impl<N> Iterable<N> for $t<N>
-    {
-      fn iter<'l>(&'l self) -> VecIterator<'l, N>
-      { self.at.iter() }
-    }
-  )
-)
-
-macro_rules! iterable_mut_impl(
-  ($t: ident) => (
-    impl<N> IterableMut<N> for $t<N>
-    {
-      fn mut_iter<'l>(&'l mut self) -> VecMutIterator<'l, N>
-      { self.at.mut_iter() }
-    }
-  )
-)
-
-macro_rules! eq_impl(
-  ($t: ident) => (
-    impl<N: Eq> Eq for $t<N>
-    {
-      #[inline]
-      fn eq(&self, other: &$t<N>) -> bool
-      { self.at.iter().zip(other.at.iter()).all(|(a, b)| a == b) }
-    
-      #[inline]
-      fn ne(&self, other: &$t<N>) -> bool
-      { self.at.iter().zip(other.at.iter()).all(|(a, b)| a != b) }
-    }
-  )
-)
-
-macro_rules! dim_impl(
-  ($t: ident, $dim: expr) => (
-    impl<N> Dim for $t<N>
-    {
-      #[inline]
-      fn dim() -> uint
-      { $dim }
-    }
-  )
-)
-
-// FIXME: add the possibility to specialize that
-macro_rules! basis_impl(
-  ($t: ident, $dim: expr) => (
-    impl<N: Copy + DivisionRing + Algebraic + ApproxEq<N>> Basis for $t<N>
-    {
-      pub fn canonical_basis() -> ~[$t<N>]
-      {
-        let mut res : ~[$t<N>] = ~[];
-    
-        for iterate(0u, $dim) |i|
-        {
-          let mut basis_element : $t<N> = Zero::zero();
-    
-          basis_element.at[i] = One::one();
-    
-          res.push(basis_element);
-        }
-    
-        res
-      }
-    
-      pub fn orthogonal_subspace_basis(&self) -> ~[$t<N>]
-      {
-        // compute the basis of the orthogonal subspace using Gram-Schmidt
-        // orthogonalization algorithm
-        let mut res : ~[$t<N>] = ~[];
-    
-        for iterate(0u, $dim) |i|
-        {
-          let mut basis_element : $t<N> = Zero::zero();
-    
-          basis_element.at[i] = One::one();
-    
-          if res.len() == $dim - 1
-          { break; }
-    
-          let mut elt = copy basis_element;
-    
-          elt = elt - self.scalar_mul(&basis_element.dot(self));
-    
-          for res.iter().advance |v|
-          { elt = elt - v.scalar_mul(&elt.dot(v)) };
-    
-          if !elt.sqnorm().approx_eq(&Zero::zero())
-          { res.push(elt.normalized()); }
-        }
-    
-        res
-      }
-    }
-  )
-)
-
-macro_rules! add_impl(
-  ($t: ident) => (
-    impl<N: Copy + Add<N,N>> Add<$t<N>, $t<N>> for $t<N>
-    {
-      #[inline]
-      fn add(&self, other: &$t<N>) -> $t<N>
-      {
-        self.at.iter()
-               .zip(other.at.iter())
-               .transform(|(a, b)| { *a + *b })
-               .collect()
-      }
-    }
-  )
-)
-
-macro_rules! sub_impl(
-  ($t: ident) => (
-    impl<N: Copy + Sub<N,N>> Sub<$t<N>, $t<N>> for $t<N>
-    {
-      #[inline]
-      fn sub(&self, other: &$t<N>) -> $t<N>
-      {
-        self.at.iter()
-               .zip(other.at.iter())
-               .transform(| (a, b) | { *a - *b })
-               .collect()
-      }
-    }
-  )
-)
-
-macro_rules! neg_impl(
-  ($t: ident) => (
-    impl<N: Neg<N>> Neg<$t<N>> for $t<N>
-    {
-      #[inline]
-      fn neg(&self) -> $t<N>
-      { self.at.iter().transform(|a| -a).collect() }
-    }
-  )
-)
-
-macro_rules! dot_impl(
-  ($t: ident, $dim: expr) => (
-    impl<N: Ring> Dot<N> for $t<N>
-    {
-      #[inline]
-      fn dot(&self, other: &$t<N>) -> N
-      {
-        let mut res = Zero::zero::<N>();
-    
-        for iterate(0u, $dim) |i|
-        { res = res + self.at[i] * other.at[i]; }
-    
-        res
-      } 
-    }
-  )
-)
-
-macro_rules! sub_dot_impl(
-  ($t: ident, $dim: expr) => (
-    impl<N: Ring> SubDot<N> for $t<N>
-    {
-      #[inline]
-      fn sub_dot(&self, a: &$t<N>, b: &$t<N>) -> N
-      {
-        let mut res = Zero::zero::<N>();
-    
-        for iterate(0u, $dim) |i|
-        { res = res + (self.at[i] - a.at[i]) * b.at[i]; }
-    
-        res
-      } 
-    }
-  )
-)
-
-macro_rules! scalar_mul_impl(
-  ($t: ident, $dim: expr) => (
-    impl<N: Mul<N, N>> ScalarMul<N> for $t<N>
-    {
-      #[inline]
-      fn scalar_mul(&self, s: &N) -> $t<N>
-      { self.at.iter().transform(|a| a * *s).collect() }
-    
-      #[inline]
-      fn scalar_mul_inplace(&mut self, s: &N)
-      {
-        for iterate(0u, $dim) |i|
-        { self.at[i] = self.at[i] * *s; }
-      }
-    }
-  )
-)
-
-
-macro_rules! scalar_div_impl(
-  ($t: ident, $dim: expr) => (
-    impl<N: Div<N, N>> ScalarDiv<N> for $t<N>
-    {
-      #[inline]
-      fn scalar_div(&self, s: &N) -> $t<N>
-      { self.at.iter().transform(|a| a / *s).collect() }
-    
-      #[inline]
-      fn scalar_div_inplace(&mut self, s: &N)
-      {
-        for iterate(0u, $dim) |i|
-        { self.at[i] = self.at[i] / *s; }
-      }
-    }
-  )
-)
-
-macro_rules! scalar_add_impl(
-  ($t: ident, $dim: expr) => (
-    impl<N: Add<N, N>> ScalarAdd<N> for $t<N>
-    {
-      #[inline]
-      fn scalar_add(&self, s: &N) -> $t<N>
-      { self.at.iter().transform(|a| a + *s).collect() }
-    
-      #[inline]
-      fn scalar_add_inplace(&mut self, s: &N)
-      {
-        for iterate(0u, $dim) |i|
-        { self.at[i] = self.at[i] + *s; }
-      }
-    }
-  )
-)
-
-macro_rules! scalar_sub_impl(
-  ($t: ident, $dim: expr) => (
-    impl<N: Sub<N, N>> ScalarSub<N> for $t<N>
-    {
-      #[inline]
-      fn scalar_sub(&self, s: &N) -> $t<N>
-      { self.at.iter().transform(|a| a - *s).collect() }
-    
-      #[inline]
-      fn scalar_sub_inplace(&mut self, s: &N)
-      {
-        for iterate(0u, $dim) |i|
-        { self.at[i] = self.at[i] - *s; }
-      }
-    }
-  )
-)
-
-macro_rules! translation_impl(
-  ($t: ident) => (
-    impl<N: Copy + Add<N, N> + Neg<N>> Translation<$t<N>> for $t<N>
-    {
-      #[inline]
-      fn translation(&self) -> $t<N>
-      { copy *self }
-
-      #[inline]
-      fn inv_translation(&self) -> $t<N>
-      { -self }
-    
-      #[inline]
-      fn translate_by(&mut self, t: &$t<N>)
-      { *self = *self + *t; }
-    }
-  )
-)
-
-macro_rules! translatable_impl(
-  ($t: ident) => (
-    impl<N: Add<N, N> + Copy> Translatable<$t<N>, $t<N>> for $t<N>
-    {
-      #[inline]
-      fn translated(&self, t: &$t<N>) -> $t<N>
-      { self + *t }
-    }
-  )
-)
-
-macro_rules! norm_impl(
-  ($t: ident, $dim: expr) => (
-    impl<N: Copy + DivisionRing + Algebraic> Norm<N> for $t<N>
-    {
-      #[inline]
-      fn sqnorm(&self) -> N
-      { self.dot(self) }
-    
-      #[inline]
-      fn norm(&self) -> N
-      { self.sqnorm().sqrt() }
-    
-      #[inline]
-      fn normalized(&self) -> $t<N>
-      {
-        let mut res : $t<N> = copy *self;
-    
-        res.normalize();
-    
-        res
-      }
-    
-      #[inline]
-      fn normalize(&mut self) -> N
-      {
-        let l = self.norm();
-    
-        for iterate(0u, $dim) |i|
-        { self.at[i] = self.at[i] / l; }
-    
-        l
-      }
-    }
-  )
-)
-
-macro_rules! approx_eq_impl(
-  ($t: ident) => (
-    impl<N: ApproxEq<N>> ApproxEq<N> for $t<N>
-    {
-      #[inline]
-      fn approx_epsilon() -> N
-      { ApproxEq::approx_epsilon::<N, N>() }
-    
-      #[inline]
-      fn approx_eq(&self, other: &$t<N>) -> bool
-      {
-        let mut zip = self.at.iter().zip(other.at.iter());
-    
-        do zip.all |(a, b)| { a.approx_eq(b) }
-      }
-    
-      #[inline]
-      fn approx_eq_eps(&self, other: &$t<N>, epsilon: &N) -> bool
-      {
-        let mut zip = self.at.iter().zip(other.at.iter());
-    
-        do zip.all |(a, b)| { a.approx_eq_eps(b, epsilon) }
-      }
-    }
-  )
-)
-
-macro_rules! zero_impl(
-  ($t: ident) => (
-    impl<N: Copy + Zero> Zero for $t<N>
-    {
-      #[inline]
-      fn zero() -> $t<N>
-      { $t::new_repeat(Zero::zero()) }
-    
-      #[inline]
-      fn is_zero(&self) -> bool
-      { self.at.iter().all(|e| e.is_zero()) }
-    }
-  )
-)
-
-macro_rules! rand_impl(
-  ($t: ident, $param: ident, [ $($elem: ident)|+ ]) => (
-    impl<N: Rand> Rand for $t<N>
-    {
-      #[inline]
-      fn rand<R: Rng>($param: &mut R) -> $t<N>
-      { $t::new([ $( $elem.gen(), )+ ]) }
-    }
-  )
-)
-
-macro_rules! from_any_iterator_impl(
-  ($t: ident, $param: ident, [ $($elem: ident)|+ ]) => (
-    impl<N: Copy> FromAnyIterator<N> for $t<N>
-    {
-      fn from_iterator<'l>($param: &mut VecIterator<'l, N>) -> $t<N>
-      { $t { at: [ $( copy *$elem.next().unwrap(), )+ ] } }
-
-      fn from_mut_iterator<'l>($param: &mut VecMutIterator<'l, N>) -> $t<N>
-      { $t { at: [ $( copy *$elem.next().unwrap(), )+ ] } }
-    }
-  )
-)
-
-macro_rules! from_iterator_impl(
-  ($t: ident, $param: ident, [ $($elem: ident)|+ ]) => (
-    impl<N, Iter: Iterator<N>> FromIterator<N, Iter> for $t<N>
-    {
-      fn from_iterator($param: &mut Iter) -> $t<N>
-      { $t { at: [ $( $elem.next().unwrap(), )+ ] } }
-    }
-  )
-)
-
-macro_rules! bounded_impl(
-  ($t: ident) => (
-    impl<N: Bounded + Copy> Bounded for $t<N>
-    {
-      #[inline]
-      fn max_value() -> $t<N>
-      { $t::new_repeat(Bounded::max_value()) }
-    
-      #[inline]
-      fn min_value() -> $t<N>
-      { $t::new_repeat(Bounded::min_value()) }
-    }
-  )
-)
 
 #[deriving(Ord, ToStr)]
 pub struct Vec1<N>
@@ -632,74 +197,3 @@ from_any_iterator_impl!(Vec6, iterator, [iterator | iterator | iterator | iterat
 bounded_impl!(Vec6)
 iterable_impl!(Vec6)
 iterable_mut_impl!(Vec6)
-
-impl<N: Mul<N, N> + Sub<N, N>> Cross<Vec1<N>> for Vec2<N>
-{
-  #[inline]
-  fn cross(&self, other : &Vec2<N>) -> Vec1<N>
-  { Vec1::new([self.at[0] * other.at[1] - self.at[1] * other.at[0]]) }
-}
-
-impl<N: Mul<N, N> + Sub<N, N>> Cross<Vec3<N>> for Vec3<N>
-{
-  #[inline]
-  fn cross(&self, other : &Vec3<N>) -> Vec3<N>
-  {
-    Vec3::new(
-      [self.at[1] * other.at[2] - self.at[2] * other.at[1],
-       self.at[2] * other.at[0] - self.at[0] * other.at[2],
-       self.at[0] * other.at[1] - self.at[1] * other.at[0]]
-    )
-  }
-}
-
-impl<N: One> Basis for Vec1<N>
-{
-  #[inline]
-  fn canonical_basis() -> ~[Vec1<N>]
-  { ~[ Vec1::new([One::one()]) ] } // FIXME: this should be static
-
-  #[inline]
-  fn orthogonal_subspace_basis(&self) -> ~[Vec1<N>]
-  { ~[] }
-}
-
-impl<N: Copy + One + Zero + Neg<N>> Basis for Vec2<N>
-{
-  #[inline]
-  fn canonical_basis()     -> ~[Vec2<N>]
-  {
-    // FIXME: this should be static
-    ~[ Vec2::new([One::one(), Zero::zero()]),
-       Vec2::new([Zero::zero(), One::one()]) ]
-  }
-
-  #[inline]
-  fn orthogonal_subspace_basis(&self) -> ~[Vec2<N>]
-  { ~[ Vec2::new([-self.at[1], copy self.at[0]]) ] }
-}
-
-impl<N: Copy + DivisionRing + Ord + Algebraic>
-Basis for Vec3<N>
-{
-  #[inline]
-  fn canonical_basis() -> ~[Vec3<N>]
-  {
-    // FIXME: this should be static
-    ~[ Vec3::new([One::one(), Zero::zero(), Zero::zero()]),
-       Vec3::new([Zero::zero(), One::one(), Zero::zero()]),
-       Vec3::new([Zero::zero(), Zero::zero(), One::one()]) ]
-  }
-
-  #[inline]
-  fn orthogonal_subspace_basis(&self) -> ~[Vec3<N>]
-  {
-      let a = 
-        if abs(copy self.at[0]) > abs(copy self.at[1])
-        { Vec3::new([copy self.at[2], Zero::zero(), -copy self.at[0]]).normalized() }
-        else
-        { Vec3::new([Zero::zero(), -self.at[2], copy self.at[1]]).normalized() };
-
-      ~[ a.cross(self), a ]
-  }
-}
diff --git a/src/vec_impl.rs b/src/vec_impl.rs
new file mode 100644
index 00000000..787dbe49
--- /dev/null
+++ b/src/vec_impl.rs
@@ -0,0 +1,431 @@
+#[macro_escape];
+
+macro_rules! new_impl(
+  ($t: ident, $dim: expr) => (
+    impl<N> $t<N>
+    {
+      #[inline]
+      pub fn new(at: [N, ..$dim]) -> $t<N>
+      { $t { at: at } }
+    }
+  )
+)
+
+macro_rules! new_repeat_impl(
+  ($t: ident, $param: ident, [ $($elem: ident)|+ ]) => (
+    impl<N: Copy> $t<N>
+    {
+      #[inline]
+      pub fn new_repeat($param: N) -> $t<N>
+      { $t{ at: [ $( copy $elem, )+ ] } }
+    }
+  )
+)
+
+macro_rules! iterable_impl(
+  ($t: ident) => (
+    impl<N> Iterable<N> for $t<N>
+    {
+      fn iter<'l>(&'l self) -> VecIterator<'l, N>
+      { self.at.iter() }
+    }
+  )
+)
+
+macro_rules! iterable_mut_impl(
+  ($t: ident) => (
+    impl<N> IterableMut<N> for $t<N>
+    {
+      fn mut_iter<'l>(&'l mut self) -> VecMutIterator<'l, N>
+      { self.at.mut_iter() }
+    }
+  )
+)
+
+macro_rules! eq_impl(
+  ($t: ident) => (
+    impl<N: Eq> Eq for $t<N>
+    {
+      #[inline]
+      fn eq(&self, other: &$t<N>) -> bool
+      { self.at.iter().zip(other.at.iter()).all(|(a, b)| a == b) }
+    
+      #[inline]
+      fn ne(&self, other: &$t<N>) -> bool
+      { self.at.iter().zip(other.at.iter()).all(|(a, b)| a != b) }
+    }
+  )
+)
+
+macro_rules! dim_impl(
+  ($t: ident, $dim: expr) => (
+    impl<N> Dim for $t<N>
+    {
+      #[inline]
+      fn dim() -> uint
+      { $dim }
+    }
+  )
+)
+
+// FIXME: add the possibility to specialize that
+macro_rules! basis_impl(
+  ($t: ident, $dim: expr) => (
+    impl<N: Copy + DivisionRing + Algebraic + ApproxEq<N>> Basis for $t<N>
+    {
+      pub fn canonical_basis() -> ~[$t<N>]
+      {
+        let mut res : ~[$t<N>] = ~[];
+    
+        for iterate(0u, $dim) |i|
+        {
+          let mut basis_element : $t<N> = Zero::zero();
+    
+          basis_element.at[i] = One::one();
+    
+          res.push(basis_element);
+        }
+    
+        res
+      }
+    
+      pub fn orthogonal_subspace_basis(&self) -> ~[$t<N>]
+      {
+        // compute the basis of the orthogonal subspace using Gram-Schmidt
+        // orthogonalization algorithm
+        let mut res : ~[$t<N>] = ~[];
+    
+        for iterate(0u, $dim) |i|
+        {
+          let mut basis_element : $t<N> = Zero::zero();
+    
+          basis_element.at[i] = One::one();
+    
+          if res.len() == $dim - 1
+          { break; }
+    
+          let mut elt = copy basis_element;
+    
+          elt = elt - self.scalar_mul(&basis_element.dot(self));
+    
+          for res.iter().advance |v|
+          { elt = elt - v.scalar_mul(&elt.dot(v)) };
+    
+          if !elt.sqnorm().approx_eq(&Zero::zero())
+          { res.push(elt.normalized()); }
+        }
+    
+        res
+      }
+    }
+  )
+)
+
+macro_rules! add_impl(
+  ($t: ident) => (
+    impl<N: Copy + Add<N,N>> Add<$t<N>, $t<N>> for $t<N>
+    {
+      #[inline]
+      fn add(&self, other: &$t<N>) -> $t<N>
+      {
+        self.at.iter()
+               .zip(other.at.iter())
+               .transform(|(a, b)| { *a + *b })
+               .collect()
+      }
+    }
+  )
+)
+
+macro_rules! sub_impl(
+  ($t: ident) => (
+    impl<N: Copy + Sub<N,N>> Sub<$t<N>, $t<N>> for $t<N>
+    {
+      #[inline]
+      fn sub(&self, other: &$t<N>) -> $t<N>
+      {
+        self.at.iter()
+               .zip(other.at.iter())
+               .transform(| (a, b) | { *a - *b })
+               .collect()
+      }
+    }
+  )
+)
+
+macro_rules! neg_impl(
+  ($t: ident) => (
+    impl<N: Neg<N>> Neg<$t<N>> for $t<N>
+    {
+      #[inline]
+      fn neg(&self) -> $t<N>
+      { self.at.iter().transform(|a| -a).collect() }
+    }
+  )
+)
+
+macro_rules! dot_impl(
+  ($t: ident, $dim: expr) => (
+    impl<N: Ring> Dot<N> for $t<N>
+    {
+      #[inline]
+      fn dot(&self, other: &$t<N>) -> N
+      {
+        let mut res = Zero::zero::<N>();
+    
+        for iterate(0u, $dim) |i|
+        { res = res + self.at[i] * other.at[i]; }
+    
+        res
+      } 
+    }
+  )
+)
+
+macro_rules! sub_dot_impl(
+  ($t: ident, $dim: expr) => (
+    impl<N: Ring> SubDot<N> for $t<N>
+    {
+      #[inline]
+      fn sub_dot(&self, a: &$t<N>, b: &$t<N>) -> N
+      {
+        let mut res = Zero::zero::<N>();
+    
+        for iterate(0u, $dim) |i|
+        { res = res + (self.at[i] - a.at[i]) * b.at[i]; }
+    
+        res
+      } 
+    }
+  )
+)
+
+macro_rules! scalar_mul_impl(
+  ($t: ident, $dim: expr) => (
+    impl<N: Mul<N, N>> ScalarMul<N> for $t<N>
+    {
+      #[inline]
+      fn scalar_mul(&self, s: &N) -> $t<N>
+      { self.at.iter().transform(|a| a * *s).collect() }
+    
+      #[inline]
+      fn scalar_mul_inplace(&mut self, s: &N)
+      {
+        for iterate(0u, $dim) |i|
+        { self.at[i] = self.at[i] * *s; }
+      }
+    }
+  )
+)
+
+
+macro_rules! scalar_div_impl(
+  ($t: ident, $dim: expr) => (
+    impl<N: Div<N, N>> ScalarDiv<N> for $t<N>
+    {
+      #[inline]
+      fn scalar_div(&self, s: &N) -> $t<N>
+      { self.at.iter().transform(|a| a / *s).collect() }
+    
+      #[inline]
+      fn scalar_div_inplace(&mut self, s: &N)
+      {
+        for iterate(0u, $dim) |i|
+        { self.at[i] = self.at[i] / *s; }
+      }
+    }
+  )
+)
+
+macro_rules! scalar_add_impl(
+  ($t: ident, $dim: expr) => (
+    impl<N: Add<N, N>> ScalarAdd<N> for $t<N>
+    {
+      #[inline]
+      fn scalar_add(&self, s: &N) -> $t<N>
+      { self.at.iter().transform(|a| a + *s).collect() }
+    
+      #[inline]
+      fn scalar_add_inplace(&mut self, s: &N)
+      {
+        for iterate(0u, $dim) |i|
+        { self.at[i] = self.at[i] + *s; }
+      }
+    }
+  )
+)
+
+macro_rules! scalar_sub_impl(
+  ($t: ident, $dim: expr) => (
+    impl<N: Sub<N, N>> ScalarSub<N> for $t<N>
+    {
+      #[inline]
+      fn scalar_sub(&self, s: &N) -> $t<N>
+      { self.at.iter().transform(|a| a - *s).collect() }
+    
+      #[inline]
+      fn scalar_sub_inplace(&mut self, s: &N)
+      {
+        for iterate(0u, $dim) |i|
+        { self.at[i] = self.at[i] - *s; }
+      }
+    }
+  )
+)
+
+macro_rules! translation_impl(
+  ($t: ident) => (
+    impl<N: Copy + Add<N, N> + Neg<N>> Translation<$t<N>> for $t<N>
+    {
+      #[inline]
+      fn translation(&self) -> $t<N>
+      { copy *self }
+
+      #[inline]
+      fn inv_translation(&self) -> $t<N>
+      { -self }
+    
+      #[inline]
+      fn translate_by(&mut self, t: &$t<N>)
+      { *self = *self + *t; }
+    }
+  )
+)
+
+macro_rules! translatable_impl(
+  ($t: ident) => (
+    impl<N: Add<N, N> + Copy> Translatable<$t<N>, $t<N>> for $t<N>
+    {
+      #[inline]
+      fn translated(&self, t: &$t<N>) -> $t<N>
+      { self + *t }
+    }
+  )
+)
+
+macro_rules! norm_impl(
+  ($t: ident, $dim: expr) => (
+    impl<N: Copy + DivisionRing + Algebraic> Norm<N> for $t<N>
+    {
+      #[inline]
+      fn sqnorm(&self) -> N
+      { self.dot(self) }
+    
+      #[inline]
+      fn norm(&self) -> N
+      { self.sqnorm().sqrt() }
+    
+      #[inline]
+      fn normalized(&self) -> $t<N>
+      {
+        let mut res : $t<N> = copy *self;
+    
+        res.normalize();
+    
+        res
+      }
+    
+      #[inline]
+      fn normalize(&mut self) -> N
+      {
+        let l = self.norm();
+    
+        for iterate(0u, $dim) |i|
+        { self.at[i] = self.at[i] / l; }
+    
+        l
+      }
+    }
+  )
+)
+
+macro_rules! approx_eq_impl(
+  ($t: ident) => (
+    impl<N: ApproxEq<N>> ApproxEq<N> for $t<N>
+    {
+      #[inline]
+      fn approx_epsilon() -> N
+      { ApproxEq::approx_epsilon::<N, N>() }
+    
+      #[inline]
+      fn approx_eq(&self, other: &$t<N>) -> bool
+      {
+        let mut zip = self.at.iter().zip(other.at.iter());
+    
+        do zip.all |(a, b)| { a.approx_eq(b) }
+      }
+    
+      #[inline]
+      fn approx_eq_eps(&self, other: &$t<N>, epsilon: &N) -> bool
+      {
+        let mut zip = self.at.iter().zip(other.at.iter());
+    
+        do zip.all |(a, b)| { a.approx_eq_eps(b, epsilon) }
+      }
+    }
+  )
+)
+
+macro_rules! zero_impl(
+  ($t: ident) => (
+    impl<N: Copy + Zero> Zero for $t<N>
+    {
+      #[inline]
+      fn zero() -> $t<N>
+      { $t::new_repeat(Zero::zero()) }
+    
+      #[inline]
+      fn is_zero(&self) -> bool
+      { self.at.iter().all(|e| e.is_zero()) }
+    }
+  )
+)
+
+macro_rules! rand_impl(
+  ($t: ident, $param: ident, [ $($elem: ident)|+ ]) => (
+    impl<N: Rand> Rand for $t<N>
+    {
+      #[inline]
+      fn rand<R: Rng>($param: &mut R) -> $t<N>
+      { $t::new([ $( $elem.gen(), )+ ]) }
+    }
+  )
+)
+
+macro_rules! from_any_iterator_impl(
+  ($t: ident, $param: ident, [ $($elem: ident)|+ ]) => (
+    impl<N: Copy> FromAnyIterator<N> for $t<N>
+    {
+      fn from_iterator<'l>($param: &mut VecIterator<'l, N>) -> $t<N>
+      { $t { at: [ $( copy *$elem.next().unwrap(), )+ ] } }
+
+      fn from_mut_iterator<'l>($param: &mut VecMutIterator<'l, N>) -> $t<N>
+      { $t { at: [ $( copy *$elem.next().unwrap(), )+ ] } }
+    }
+  )
+)
+
+macro_rules! from_iterator_impl(
+  ($t: ident, $param: ident, [ $($elem: ident)|+ ]) => (
+    impl<N, Iter: Iterator<N>> FromIterator<N, Iter> for $t<N>
+    {
+      fn from_iterator($param: &mut Iter) -> $t<N>
+      { $t { at: [ $( $elem.next().unwrap(), )+ ] } }
+    }
+  )
+)
+
+macro_rules! bounded_impl(
+  ($t: ident) => (
+    impl<N: Bounded + Copy> Bounded for $t<N>
+    {
+      #[inline]
+      fn max_value() -> $t<N>
+      { $t::new_repeat(Bounded::max_value()) }
+    
+      #[inline]
+      fn min_value() -> $t<N>
+      { $t::new_repeat(Bounded::min_value()) }
+    }
+  )
+)
diff --git a/src/vec_spec.rs b/src/vec_spec.rs
new file mode 100644
index 00000000..6d663443
--- /dev/null
+++ b/src/vec_spec.rs
@@ -0,0 +1,77 @@
+use std::num::{Zero, One, abs};
+use traits::basis::Basis;
+use traits::cross::Cross;
+use traits::division_ring::DivisionRing;
+use traits::norm::Norm;
+use vec::{Vec1, Vec2, Vec3};
+
+impl<N: Mul<N, N> + Sub<N, N>> Cross<Vec1<N>> for Vec2<N>
+{
+  #[inline]
+  fn cross(&self, other : &Vec2<N>) -> Vec1<N>
+  { Vec1::new([self.at[0] * other.at[1] - self.at[1] * other.at[0]]) }
+}
+
+impl<N: Mul<N, N> + Sub<N, N>> Cross<Vec3<N>> for Vec3<N>
+{
+  #[inline]
+  fn cross(&self, other : &Vec3<N>) -> Vec3<N>
+  {
+    Vec3::new(
+      [self.at[1] * other.at[2] - self.at[2] * other.at[1],
+       self.at[2] * other.at[0] - self.at[0] * other.at[2],
+       self.at[0] * other.at[1] - self.at[1] * other.at[0]]
+    )
+  }
+}
+
+impl<N: One> Basis for Vec1<N>
+{
+  #[inline]
+  fn canonical_basis() -> ~[Vec1<N>]
+  { ~[ Vec1::new([One::one()]) ] } // FIXME: this should be static
+
+  #[inline]
+  fn orthogonal_subspace_basis(&self) -> ~[Vec1<N>]
+  { ~[] }
+}
+
+impl<N: Copy + One + Zero + Neg<N>> Basis for Vec2<N>
+{
+  #[inline]
+  fn canonical_basis()     -> ~[Vec2<N>]
+  {
+    // FIXME: this should be static
+    ~[ Vec2::new([One::one(), Zero::zero()]),
+       Vec2::new([Zero::zero(), One::one()]) ]
+  }
+
+  #[inline]
+  fn orthogonal_subspace_basis(&self) -> ~[Vec2<N>]
+  { ~[ Vec2::new([-self.at[1], copy self.at[0]]) ] }
+}
+
+impl<N: Copy + DivisionRing + Ord + Algebraic>
+Basis for Vec3<N>
+{
+  #[inline]
+  fn canonical_basis() -> ~[Vec3<N>]
+  {
+    // FIXME: this should be static
+    ~[ Vec3::new([One::one(), Zero::zero(), Zero::zero()]),
+       Vec3::new([Zero::zero(), One::one(), Zero::zero()]),
+       Vec3::new([Zero::zero(), Zero::zero(), One::one()]) ]
+  }
+
+  #[inline]
+  fn orthogonal_subspace_basis(&self) -> ~[Vec3<N>]
+  {
+      let a = 
+        if abs(copy self.at[0]) > abs(copy self.at[1])
+        { Vec3::new([copy self.at[2], Zero::zero(), -copy self.at[0]]).normalized() }
+        else
+        { Vec3::new([Zero::zero(), -self.at[2], copy self.at[1]]).normalized() };
+
+      ~[ a.cross(self), a ]
+  }
+}