diff --git a/src/funs/boolv.rs b/src/funs/boolv.rs
index 7578199..f5a48a7 100644
--- a/src/funs/boolv.rs
+++ b/src/funs/boolv.rs
@@ -1,4 +1,4 @@
-use vector::{Vector, Vec2, Vec3, Vec4};
+use vec::{Vector, Vec2, Vec3, Vec4};
 
 pub trait BooleanVector: Vector<bool> {
     pure fn any() -> bool;
diff --git a/src/funs/common.rs b/src/funs/common.rs
index bf292f9..f9fde37 100644
--- a/src/funs/common.rs
+++ b/src/funs/common.rs
@@ -1,7 +1,7 @@
 use core::cmp::{Eq, Ord};
 
 use num::cast::*;
-use vector::{Vec2, Vec3, Vec4};
+use vec::{Vec2, Vec3, Vec4};
 
 /**
  * Common Functions for all numeric types
diff --git a/src/funs/exp.rs b/src/funs/exp.rs
index cddb288..29c7003 100644
--- a/src/funs/exp.rs
+++ b/src/funs/exp.rs
@@ -2,7 +2,7 @@
  * Exponential Functions
  */
 use num::cast::*;
-use vector::{Vec2, Vec3, Vec4};
+use vec::{Vec2, Vec3, Vec4};
 
 pub trait Exp {
     pure fn pow<N:NumCast>(n: &N) -> self;
diff --git a/src/funs/projection.rs b/src/funs/projection.rs
index 6fd804d..3200987 100644
--- a/src/funs/projection.rs
+++ b/src/funs/projection.rs
@@ -1,8 +1,8 @@
-use core::float::consts::pi;
-use core::float::tan;
-
-use matrix::Mat4;
-use num::cast::*;
+use funs::trig::*;
+use mat::Mat4;
+use num::cast::cast;
+use num::consts::pi;
+use num::ext::FloatExt;
 
 //
 //  Create a perspective projection matrix
@@ -11,8 +11,8 @@ use num::cast::*;
 //  http://www.opengl.org/wiki/GluPerspective_code
 //
 #[inline(always)]
-pure fn perspective<T:Copy NumCast Neg<T>>(fovy: float, aspectRatio: float, near: float, far: float) -> Mat4<T> {
-    let ymax = near * tan(fovy * pi / 360f);
+pure fn perspective<T:Copy FloatExt Trig>(fovy: T, aspectRatio: T, near: T, far: T) -> Mat4<T> {
+    let ymax = near * tan(&(fovy * pi() / cast(360)));
     let xmax = ymax * aspectRatio;
     return frustum(-xmax, xmax, -ymax, ymax, near, far);
 }
@@ -25,24 +25,26 @@ pure fn perspective<T:Copy NumCast Neg<T>>(fovy: float, aspectRatio: float, near
 //  TODO: double check algorithm
 //
 #[inline(always)]
-pure fn frustum<T:Copy NumCast Neg<T>>(left: float, right: float, bottom: float, top: float, near: float, far: float) -> Mat4<T> {
-    let _0 = cast(0);
+pure fn frustum<T:Copy FloatExt>(left: T, right: T, bottom: T, top: T, near: T, far: T) -> Mat4<T> {
+    let _0: T = cast(0);
+    let _2: T = cast(2);
+    let neg_1 = cast(-1);
     
-    let c0r0 = cast((2f * near) / (right - left));
+    let c0r0 = (_2 * near) / (right - left);
     let c0r1 = _0;
     let c0r2 = _0;
     let c0r3 = _0;
     let c1r0 = _0;
-    let c1r1 = cast((2f * near) / (top - bottom));
+    let c1r1 = (_2 * near) / (top - bottom);
     let c1r2 = _0;
     let c1r3 = _0;
-    let c2r0 = cast((right + left) / (right - left));
-    let c2r1 = cast((top + bottom) / (top - bottom));
-    let c2r2 = cast(-(far + near) / (far - near));
-    let c2r3 = cast(-1);
+    let c2r0 = (right + left) / (right - left);
+    let c2r1 = (top + bottom) / (top - bottom);
+    let c2r2 = -(far + near) / (far - near);
+    let c2r3 = neg_1;
     let c3r0 = _0;
     let c3r1 = _0;
-    let c3r2 = cast(-(2f * far * near) / (far - near));
+    let c3r2 = -(_2 * far * near) / (far - near);
     let c3r3 = _0;
     
     return Mat4::new(c0r0, c0r1, c0r2, c0r3,
diff --git a/src/funs/relv.rs b/src/funs/relv.rs
index 6bf3b14..0ed2587 100644
--- a/src/funs/relv.rs
+++ b/src/funs/relv.rs
@@ -4,7 +4,7 @@
 
 use core::cmp::{Eq, Ord};
 
-use vector::{Vec2, Vec3, Vec4};
+use vec::{Vec2, Vec3, Vec4};
 
 pub trait RelVector<BVec> {
     pure fn less_than(y: &self) -> BVec;
diff --git a/src/funs/test/test_boolv.rs b/src/funs/test/test_boolv.rs
index 5497207..dce9ede 100644
--- a/src/funs/test/test_boolv.rs
+++ b/src/funs/test/test_boolv.rs
@@ -1,4 +1,4 @@
-use vector::*;
+use vec::*;
 use boolv::*;
 
 #[test]
diff --git a/src/funs/test/test_transform.rs b/src/funs/test/test_transform.rs
index 63fdb89..c25cbfb 100644
--- a/src/funs/test/test_transform.rs
+++ b/src/funs/test/test_transform.rs
@@ -1,6 +1,6 @@
 use funs::transform::*;
-use matrix::Mat4;
-use vector::{Vec3, Vec4};
+use mat::Mat4;
+use vec::{Vec3, Vec4};
 
 #[test]
 fn test_mat4_from_rotation() {
diff --git a/src/funs/transform.rs b/src/funs/transform.rs
index fca1acb..1efaa93 100644
--- a/src/funs/transform.rs
+++ b/src/funs/transform.rs
@@ -1,5 +1,5 @@
 use funs::trig::*;
-use matrix::{Mat3, Mat4};
+use mat::{Mat3, Mat4};
 use num::cast::*;
 
 pub pure fn mat3_from_rotation<T:Copy Num NumCast AngleConv Trig>(theta: T, axis: Vec3<T>) -> Mat3<T> {
diff --git a/src/funs/trig.rs b/src/funs/trig.rs
index 5fc7e8b..b5619bf 100644
--- a/src/funs/trig.rs
+++ b/src/funs/trig.rs
@@ -1,5 +1,5 @@
 use num::cast::*;
-use vector::{Vec3, Vec2, Vec4};
+use vec::{Vec3, Vec2, Vec4};
 
 pub trait AngleConv {
     pure fn to_degrees() -> self;
diff --git a/src/lmath.rc b/src/lmath.rc
index c436922..5f35a5a 100644
--- a/src/lmath.rc
+++ b/src/lmath.rc
@@ -10,15 +10,16 @@
 
 extern mod std;
 
-pub mod matrix;
-pub mod quaternion;
-pub mod vector;
+pub mod dim;
+pub mod mat;
+pub mod quat;
+pub mod vec;
 
 #[test]
 mod test {
-    mod test_matrix;
-    mod test_quaternion;
-    mod test_vector;
+    mod test_mat;
+    mod test_quat;
+    mod test_vec;
 }
 
 use common::*;
@@ -29,7 +30,7 @@ pub mod common {
 pub mod num {
     pub mod cast;
     pub mod consts;
-    pub mod traits;
+    pub mod ext;
 }
 
 pub mod funs {
diff --git a/src/mat.rs b/src/mat.rs
new file mode 100644
index 0000000..b7a552a
--- /dev/null
+++ b/src/mat.rs
@@ -0,0 +1,945 @@
+use core::cast::transmute;
+use core::cmp::Eq;
+use core::ptr::to_unsafe_ptr;
+use core::vec::raw::buf_as_slice;
+
+use std::cmp::FuzzyEq;
+
+use dim::Dimensional;
+use funs::common::*;
+use funs::exp::*;
+use math::*;
+use num::cast::*;
+use quat::{Quat, ToQuat};
+use vec::{Vec2, Vec3, Vec4};
+
+
+// GLSL equivalent type aliases
+
+pub type mat2 = Mat2<f32>;              /// a 2×2 single-precision floating-point matrix
+pub type mat3 = Mat3<f32>;              /// a 3×3 single-precision floating-point matrix
+pub type mat4 = Mat4<f32>;              /// a 4×4 single-precision floating-point matrix
+pub type mat2x2 = Mat2<f32>;            /// same as a `mat2`
+// pub type mat2x3 =                    /// a single-precision floating-point matrix with 2 columns and 3 rows
+// pub type mat2x4 =                    /// a single-precision floating-point matrix with 2 columns and 4 rows
+// pub type mat3x2 =                    /// a single-precision floating-point matrix with 3 columns and 2 rows
+pub type mat3x3 = Mat3<f32>;            /// same as a `mat3`
+// pub type mat3x4 =                    /// a single-precision floating-point matrix with 3 columns and 4 rows
+// pub type mat4x2 =                    /// a single-precision floating-point matrix with 4 columns and 2 rows
+// pub type mat4x3 =                    /// a single-precision floating-point matrix with 4 columns and 3 rows
+pub type mat4x4 = Mat4<f32>;            /// same as a `mat4`
+
+pub type dmat2 = Mat2<f64>;             /// a 2×2 double-precision floating-point matrix
+pub type dmat3 = Mat3<f64>;             /// a 3×3 double-precision floating-point matrix
+pub type dmat4 = Mat4<f64>;             /// a 4×4 double-precision floating-point matrix
+pub type dmat2x2 = Mat2<f64>;           /// same as a `dmat2`
+// pub type dmat2x3 =                   /// a double-precision floating-point matrix with 2 columns and 3 rows
+// pub type dmat2x4 =                   /// a double-precision floating-point matrix with 2 columns and 4 rows
+// pub type dmat3x2 =                   /// a double-precision floating-point matrix with 3 columns and 2 rows
+pub type dmat3x3 = Mat3<f64>;           /// same as a `dmat3`
+// pub type dmat3x4 =                   /// a double-precision floating-point matrix with 3 columns and 4 rows
+// pub type dmat4x2 =                   /// a double-precision floating-point matrix with 4 columns and 2 rows
+// pub type dmat4x3 =                   /// a double-precision floating-point matrix with 4 columns and 3 rows
+pub type dmat4x4 = Mat4<f64>;           /// same as a `dmat4`
+
+
+pub trait Matrix<T, Col, Row>: Dimensional<T>, Eq {
+    pure fn rows() -> uint;
+    pure fn cols() -> uint;
+    pure fn is_col_major() -> bool;
+    pure fn is_square() -> bool;
+    
+    pure fn col(i: uint) -> Col;
+    pure fn row(i: uint) -> Row;
+}
+
+pub trait NumericMatrix<T, Col, Row>: Matrix<T, Col, Row>, Neg<self> {
+    pure fn mul_t(value: T) -> self;
+    pure fn mul_v(other: &Col) -> Col;
+    pure fn add_m(other: &self) -> self;
+    pure fn sub_m(other: &self) -> self;
+}
+
+pub trait NumericMatrix_NxN<T, ColRow>: NumericMatrix<T, ColRow, ColRow> {
+    pure fn mul_m(other: &self) -> self;
+    
+    pure fn det() -> T;
+    
+    pure fn invert() -> Option<self>;
+    pure fn transpose() -> self;
+    
+    pure fn is_identity() -> bool;
+    pure fn is_symmetric() -> bool;
+    pure fn is_diagonal() -> bool;
+    pure fn is_rotated() -> bool;
+    pure fn is_invertible() -> bool;
+}
+
+pub trait Matrix2<T>: Matrix<T, Mat2<T>, Mat2<T>> {
+    pure fn to_Mat3() -> Mat3<T>;
+    pure fn to_Mat4() -> Mat4<T>;
+}
+
+pub trait Matrix3<T>: Matrix<T, Mat3<T>, Mat3<T>> {
+    pure fn to_Mat4() -> Mat4<T>;
+}
+
+pub trait Matrix4<T>: Matrix<T, Mat4<T>, Mat4<T>> {
+    
+}
+
+
+
+
+
+
+//
+//  Mat2: A 2x2, column major matrix
+//
+pub struct Mat2<T> { x: Vec2<T>, y: Vec2<T> }
+
+pub mod Mat2 {
+    
+    #[inline(always)]
+    pub pure fn new<T>(c0r0: T, c0r1: T,
+                       c1r0: T, c1r1: T) -> Mat2<T> {
+        Mat2::from_cols(Vec2::new(move c0r0, move c0r1),
+                        Vec2::new(move c1r0, move c1r1))
+    }
+
+    #[inline(always)]
+    pub pure fn from_cols<T>(c0: Vec2<T>, c1: Vec2<T>) -> Mat2<T> {
+        Mat2 { x: move c0,
+               y: move c1 }
+    }
+    
+    #[inline(always)]
+    pub pure fn from_value<T:Copy NumCast>(value: T) -> Mat2<T> {
+        let _0 = cast(0);
+        Mat2::new(value,    _0,
+                     _0, value)
+    }
+    
+    #[inline(always)]
+    pub pure fn zero<T:Copy NumCast>() -> Mat2<T> {
+        let _0 = cast(0);
+        Mat2::new(_0, _0,
+                  _0, _0)
+    }
+    
+    #[inline(always)]
+    pub pure fn identity<T:Copy NumCast>() -> Mat2<T> {
+        let _0 = cast(0);
+        let _1 = cast(1);
+        Mat2::new(_1, _0,
+                  _0, _1)
+    }
+}
+
+pub impl<T:Copy> Mat2<T>: Matrix<T, Vec2<T>, Vec2<T>> {
+    #[inline(always)]
+    static pure fn dim() -> uint { 2 }
+    
+    #[inline(always)]
+    pure fn rows() -> uint { 2 }
+    
+    #[inline(always)]
+    pure fn cols() -> uint { 2 }
+    
+    #[inline(always)]
+    pure fn is_col_major() -> bool { true }
+    
+    #[inline(always)]
+    pure fn is_square() -> bool { true }
+    
+    #[inline(always)]
+    pure fn col(i: uint) -> Vec2<T> { self[i] }
+    
+    #[inline(always)]
+    pure fn row(i: uint) -> Vec2<T> {
+        Vec2::new(self[0][i],
+                  self[1][i])
+    }
+    
+    #[inline(always)]
+    pure fn index(i: uint) -> Vec2<T> {
+        unsafe { do buf_as_slice(
+            transmute::<*Mat2<T>, *Vec2<T>>(
+                to_unsafe_ptr(&self)), 2) |slice| { slice[i] }
+        }
+    }
+    
+    #[inline(always)]
+    pure fn to_ptr() -> *T {
+        self[0].to_ptr()
+    }
+}
+
+pub impl<T:Copy Num> Mat2<T>: NumericMatrix<T, Vec2<T>, Vec2<T>> {
+    #[inline(always)]
+    pure fn neg() -> Mat2<T> {
+        Mat2::from_cols(-self[0], -self[1])
+    }
+    
+    #[inline(always)]
+    pure fn mul_t(value: T) -> Mat2<T> {
+        Mat2::from_cols(self[0].mul_t(value),
+                        self[1].mul_t(value))
+    }
+    
+    #[inline(always)]
+    pure fn mul_v(other: &Vec2<T>) -> Vec2<T> {
+        Vec2::new(self.row(0).dot(other),
+                  self.row(1).dot(other))
+    }
+    
+    #[inline(always)]
+    pure fn add_m(other: &Mat2<T>) -> Mat2<T> {
+        Mat2::from_cols(self[0].add_v(&other[0]),
+                        self[1].add_v(&other[1]))
+    }
+    
+    #[inline(always)]
+    pure fn sub_m(other: &Mat2<T>) -> Mat2<T> {
+        Mat2::from_cols(self[0].sub_v(&other[0]),
+                        self[1].sub_v(&other[1]))
+    }
+}
+
+pub impl<T:Copy Num NumCast FuzzyEq> Mat2<T>: NumericMatrix_NxN<T, Vec2<T>> {
+    #[inline(always)]
+    pure fn mul_m(other: &Mat2<T>) -> Mat2<T> {
+        Mat2::new(self.row(0).dot(&other.col(0)), self.row(1).dot(&other.col(0)),
+                  self.row(0).dot(&other.col(1)), self.row(1).dot(&other.col(1)))
+    }
+    
+    pure fn det() -> T {
+       self[0][0]*self[1][1] - self[1][0]*self[0][1]
+    }
+
+    #[inline(always)]
+    pure fn invert() -> Option<Mat2<T>> {
+        let _0 = cast(0);
+        let d = self.det();
+        if d.fuzzy_eq(&_0) {
+            None
+        } else {
+            Some(Mat2::new(self[1][1]/d, -self[0][1]/d,
+                           -self[1][0]/d, self[0][0]/d))
+        }
+    }
+    
+    #[inline(always)]
+    pure fn transpose() -> Mat2<T> {
+        Mat2::new(self[0][0], self[1][0],
+                  self[0][1], self[1][1])
+    }
+    
+    #[inline(always)]
+    pure fn is_identity() -> bool {
+        self.fuzzy_eq(&Mat2::identity())
+    }
+    
+    #[inline(always)]
+    pure fn is_symmetric() -> bool {
+        self[0][1].fuzzy_eq(&self[1][0]) &&
+        self[1][0].fuzzy_eq(&self[0][1])
+    }
+    
+    #[inline(always)]
+    pure fn is_diagonal() -> bool {
+        let _0 = cast(0);
+        self[0][1].fuzzy_eq(&_0) &&
+        self[1][0].fuzzy_eq(&_0)
+    }
+    
+    #[inline(always)]
+    pure fn is_rotated() -> bool {
+        !self.fuzzy_eq(&Mat2::identity())
+    }
+
+    #[inline(always)]
+    pure fn is_invertible() -> bool {
+        let _0 = cast(0);
+        !self.det().fuzzy_eq(&_0)
+    }
+}
+
+pub impl<T:Copy NumCast> Mat2<T>: Matrix2<T> {
+    #[inline(always)]
+    pure fn to_Mat3() -> Mat3<T> {
+        Mat3::from_Mat2(&self)
+    }
+    
+    #[inline(always)]
+    pure fn to_Mat4() -> Mat4<T> {
+        Mat4::from_Mat2(&self)
+    }
+}
+
+// TODO: make work for T:Integer
+pub impl<T:Copy FuzzyEq> Mat2<T>: Eq {
+    #[inline(always)]
+    pure fn eq(other: &Mat2<T>) -> bool {
+        self.fuzzy_eq(other)
+    }
+    
+    #[inline(always)]
+    pure fn ne(other: &Mat2<T>) -> bool {
+        !(self == *other)
+    }
+}
+
+impl<T:Copy Eq> Mat2<T>: ExactEq {
+    #[inline(always)]
+    pure fn exact_eq(other: &Mat2<T>) -> bool {
+        self[0].exact_eq(&other[0]) &&
+        self[1].exact_eq(&other[1])
+    }
+}
+
+pub impl<T:Copy FuzzyEq> Mat2<T>: FuzzyEq {
+    #[inline(always)]
+    pure fn fuzzy_eq(other: &Mat2<T>) -> bool {
+        self[0].fuzzy_eq(&other[0]) &&
+        self[1].fuzzy_eq(&other[1])
+    }
+}
+
+
+
+
+
+
+//
+//  Mat3: A 3x3, column major matrix
+//
+pub struct Mat3<T> { x: Vec3<T>, y: Vec3<T>, z: Vec3<T> }
+
+pub mod Mat3 {
+    
+    #[inline(always)]
+    pub pure fn new<T>(c0r0:T, c0r1:T, c0r2:T,
+                       c1r0:T, c1r1:T, c1r2:T,
+                       c2r0:T, c2r1:T, c2r2:T) -> Mat3<T> {
+        Mat3::from_cols(Vec3::new(move c0r0, move c0r1, move c0r2),
+                        Vec3::new(move c1r0, move c1r1, move c1r2),
+                        Vec3::new(move c2r0, move c2r1, move c2r2))
+    }
+    
+    #[inline(always)]
+    pub pure fn from_cols<T>(c0: Vec3<T>, c1: Vec3<T>, c2: Vec3<T>) -> Mat3<T> {
+        Mat3 { x: move c0,
+               y: move c1,
+               z: move c2 }
+    }
+    
+    #[inline(always)]
+    pub pure fn from_value<T:Copy NumCast>(value: T) -> Mat3<T> {
+        let _0 = cast(0);
+        Mat3::new(value,    _0,    _0,
+                     _0, value,    _0,
+                     _0,    _0, value)
+    }
+    
+    #[inline(always)]
+    pub pure fn from_Mat2<T:Copy NumCast>(m: &Mat2<T>) -> Mat3<T> {
+        let _0 = cast(0);
+        let _1 = cast(1);
+        Mat3::new(m[0][0], m[0][1], _0,
+                  m[1][0], m[1][1], _0,
+                       _0,      _0, _1)
+    }
+    
+    #[inline(always)]
+    pub pure fn zero<T:Copy NumCast>() -> Mat3<T> {
+        let _0 = cast(0);
+        Mat3::new(_0, _0, _0,
+                  _0, _0, _0,
+                  _0, _0, _0)
+    }
+    
+    #[inline(always)]
+    pub pure fn identity<T:Copy NumCast>() -> Mat3<T> {
+        let _0 = cast(0);
+        let _1 = cast(1);
+        Mat3::new(_1, _0, _0,
+                  _0, _1, _0,
+                  _0, _0, _1)
+    }
+}
+
+pub impl<T:Copy> Mat3<T>: Matrix<T, Vec3<T>, Vec3<T>> {
+    #[inline(always)]
+    static pure fn dim() -> uint { 3 }
+    
+    #[inline(always)]
+    pure fn rows() -> uint { 3 }
+    
+    #[inline(always)]
+    pure fn cols() -> uint { 3 }
+    
+    #[inline(always)]
+    pure fn is_col_major() -> bool { true }
+    
+    #[inline(always)]
+    pure fn is_square() -> bool { true }
+    
+    #[inline(always)]
+    pure fn col(i: uint) -> Vec3<T> { self[i] }
+    
+    #[inline(always)]
+    pure fn row(i: uint) -> Vec3<T> {
+        Vec3::new(self[0][i],
+                  self[1][i],
+                  self[2][i])
+    }
+    
+    #[inline(always)]
+    pure fn index(i: uint) -> Vec3<T> {
+        unsafe { do buf_as_slice(
+            transmute::<*Mat3<T>, *Vec3<T>>(
+                to_unsafe_ptr(&self)), 3) |slice| { slice[i] }
+        }
+    }
+    
+    #[inline(always)]
+    pure fn to_ptr() -> *T {
+        self[0].to_ptr()
+    }
+}
+
+pub impl<T:Copy Num> Mat3<T>: NumericMatrix<T, Vec3<T>, Vec3<T>> {
+    #[inline(always)]
+    pure fn neg() -> Mat3<T> {
+        Mat3::from_cols(-self[0], -self[1], -self[2])
+    }
+    
+    #[inline(always)]
+    pure fn mul_t(value: T) -> Mat3<T> {
+        Mat3::from_cols(self[0].mul_t(value),
+                        self[1].mul_t(value),
+                        self[2].mul_t(value))
+    }
+    
+    #[inline(always)]
+    pure fn mul_v(other: &Vec3<T>) -> Vec3<T> {
+        Vec3::new(self.row(0).dot(other),
+                  self.row(1).dot(other),
+                  self.row(2).dot(other))
+    }
+    
+    #[inline(always)]
+    pure fn add_m(other: &Mat3<T>) -> Mat3<T> {
+        Mat3::from_cols(self[0].add_v(&other[0]),
+                        self[1].add_v(&other[1]),
+                        self[2].add_v(&other[2]))
+    }
+    
+    #[inline(always)]
+    pure fn sub_m(other: &Mat3<T>) -> Mat3<T> {
+        Mat3::from_cols(self[0].sub_v(&other[0]),
+                        self[1].sub_v(&other[1]),
+                        self[2].sub_v(&other[2]))
+    }
+}
+
+pub impl<T:Copy Num NumCast FuzzyEq> Mat3<T>: NumericMatrix_NxN<T, Vec3<T>> {
+    #[inline(always)]
+    pure fn mul_m(other: &Mat3<T>) -> Mat3<T> {
+        Mat3::new(self.row(0).dot(&other.col(0)), self.row(1).dot(&other.col(0)), self.row(2).dot(&other.col(0)),
+                  self.row(0).dot(&other.col(1)), self.row(1).dot(&other.col(1)), self.row(2).dot(&other.col(1)),
+                  self.row(0).dot(&other.col(2)), self.row(1).dot(&other.col(2)), self.row(2).dot(&other.col(2)))
+    }
+    
+    pure fn det() -> T {
+        self.col(0).dot(&self.col(1).cross(&self.col(2)))
+    }
+
+    // #[inline(always)]
+    pure fn invert() -> Option<Mat3<T>> {
+        let d = self.det();
+        let _0 = cast(0);
+        if d.fuzzy_eq(&_0) {
+            None
+        } else {
+            Some(Mat3::from_cols(self[1].cross(&self[2]).div_t(d),
+                                 self[2].cross(&self[0]).div_t(d),
+                                 self[0].cross(&self[1]).div_t(d))
+            .transpose())
+        }
+    }
+    
+    #[inline(always)]
+    pure fn transpose() -> Mat3<T> {
+        Mat3::new(self[0][0], self[1][0], self[2][0],
+                  self[0][1], self[1][1], self[2][1],
+                  self[0][2], self[1][2], self[2][2])
+    }
+    
+    #[inline(always)]
+    pure fn is_identity() -> bool {
+        self.fuzzy_eq(&Mat3::identity())
+    }
+    
+    #[inline(always)]
+    pure fn is_symmetric() -> bool {
+        self[0][1].fuzzy_eq(&self[1][0]) &&
+        self[0][2].fuzzy_eq(&self[2][0]) &&
+        
+        self[1][0].fuzzy_eq(&self[0][1]) &&
+        self[1][2].fuzzy_eq(&self[2][1]) &&
+        
+        self[2][0].fuzzy_eq(&self[0][2]) &&
+        self[2][1].fuzzy_eq(&self[1][2])
+    }
+    
+    #[inline(always)]
+    pure fn is_diagonal() -> bool {
+        let _0 = cast(0);
+        self[0][1].fuzzy_eq(&_0) &&
+        self[0][2].fuzzy_eq(&_0) &&
+        
+        self[1][0].fuzzy_eq(&_0) &&
+        self[1][2].fuzzy_eq(&_0) &&
+        
+        self[2][0].fuzzy_eq(&_0) &&
+        self[2][1].fuzzy_eq(&_0)
+    }
+    
+    #[inline(always)]
+    pure fn is_rotated() -> bool {
+        !self.fuzzy_eq(&Mat3::identity())
+    }
+
+    #[inline(always)]
+    pure fn is_invertible() -> bool {
+        let _0 = cast(0);
+        !self.det().fuzzy_eq(&_0)
+    }
+}
+
+pub impl<T:Copy NumCast> Mat3<T>: Matrix3<T> {
+    #[inline(always)]
+    pure fn to_Mat4() -> Mat4<T> {
+        Mat4::from_Mat3(&self)
+    }
+}
+
+pub impl<T:Copy Num NumCast Ord> Mat3<T>: ToQuat<T> {
+    pure fn to_Quat() -> Quat<T> {
+        // Implemented using a mix of ideas from jMonkeyEngine and Ken Shoemake's
+        // paper on Quaternions: http://www.cs.ucr.edu/~vbz/resources/Quatut.pdf
+        
+        let mut s: float;
+        let w: float, x: float, y: float, z: float;
+        let trace: float = cast(self[0][0] + self[1][1] + self[2][2]);
+        
+        if trace >= cast(0) {
+            s = (trace + 1f).sqrt();
+            w = 0.5 * s;
+            s = 0.5 / s;
+            x = (self[1][2] - self[2][1]).cast::<float>() * s;
+            y = (self[2][0] - self[0][2]).cast::<float>() * s;
+            z = (self[0][1] - self[1][0]).cast::<float>() * s;
+        } else if (self[0][0] > self[1][1]) && (self[0][0] > self[2][2]) {
+            s = (1f + (self[0][0] - self[1][1] - self[2][2]).cast::<float>()).sqrt();
+            w = 0.5 * s;
+            s = 0.5 / s;
+            x = (self[0][1] - self[1][0]).cast::<float>() * s;
+            y = (self[2][0] - self[0][2]).cast::<float>() * s;
+            z = (self[1][2] - self[2][1]).cast::<float>() * s;
+        } else if self[1][1] > self[2][2] {
+            s = (1f + (self[1][1] - self[0][0] - self[2][2]).cast::<float>()).sqrt();
+            w = 0.5 * s;
+            s = 0.5 / s;
+            x = (self[0][1] - self[1][0]).cast::<float>() * s;
+            y = (self[1][2] - self[2][1]).cast::<float>() * s;
+            z = (self[2][0] - self[0][2]).cast::<float>() * s;
+        } else {
+            s = (1f + (self[2][2] - self[0][0] - self[1][1]).cast::<float>()).sqrt();
+            w = 0.5 * s;
+            s = 0.5 / s;
+            x = (self[2][0] - self[0][2]).cast::<float>() * s;
+            y = (self[1][2] - self[2][1]).cast::<float>() * s;
+            z = (self[0][1] - self[1][0]).cast::<float>() * s;
+        }
+        
+        Quat::new(cast(w), cast(x), cast(y), cast(z))
+    }
+}
+
+// TODO: make work for T:Integer
+pub impl<T:Copy FuzzyEq> Mat3<T>: Eq {
+    #[inline(always)]
+    pure fn eq(other: &Mat3<T>) -> bool {
+        self.fuzzy_eq(other)
+    }
+    
+    #[inline(always)]
+    pure fn ne(other: &Mat3<T>) -> bool {
+        !(self == *other)
+    }
+}
+
+pub impl<T:Copy Eq> Mat3<T>: ExactEq {
+    #[inline(always)]
+    pure fn exact_eq(other: &Mat3<T>) -> bool {
+        self[0].exact_eq(&other[0]) &&
+        self[1].exact_eq(&other[1]) &&
+        self[2].exact_eq(&other[2])
+    }
+}
+
+pub impl<T:Copy FuzzyEq> Mat3<T>: FuzzyEq {
+    #[inline(always)]
+    pure fn fuzzy_eq(other: &Mat3<T>) -> bool {
+        self[0].fuzzy_eq(&other[0]) &&
+        self[1].fuzzy_eq(&other[1]) &&
+        self[2].fuzzy_eq(&other[2])
+    }
+}
+
+
+
+
+
+
+//
+//  Mat4: A 4x4, column major matrix
+//
+pub struct Mat4<T> { x: Vec4<T>, y: Vec4<T>, z: Vec4<T>, w: Vec4<T> }
+
+pub mod Mat4 {
+    
+    #[inline(always)]
+    pub pure fn new<T>(c0r0: T, c0r1: T, c0r2: T, c0r3: T,
+                       c1r0: T, c1r1: T, c1r2: T, c1r3: T,
+                       c2r0: T, c2r1: T, c2r2: T, c2r3: T,
+                       c3r0: T, c3r1: T, c3r2: T, c3r3: T) -> Mat4<T>  {
+        Mat4::from_cols(Vec4::new(move c0r0, move c0r1, move c0r2, move c0r3),
+                        Vec4::new(move c1r0, move c1r1, move c1r2, move c1r3),
+                        Vec4::new(move c2r0, move c2r1, move c2r2, move c2r3),
+                        Vec4::new(move c3r0, move c3r1, move c3r2, move c3r3))
+    }
+    
+    #[inline(always)]
+    pub pure fn from_cols<T>(c0: Vec4<T>, c1: Vec4<T>, c2: Vec4<T>, c3: Vec4<T>) -> Mat4<T> {
+        Mat4 { x: move c0,
+               y: move c1,
+               z: move c2,
+               w: move c3 }
+    }
+    
+    #[inline(always)]
+    pub pure fn from_value<T:Copy NumCast>(value: T) -> Mat4<T> {
+        let _0 = cast(0);
+        Mat4::new(value,    _0,    _0,    _0,
+                     _0, value,    _0,    _0,
+                     _0,    _0, value,    _0,
+                     _0,    _0,    _0, value)
+    }
+    
+    #[inline(always)]
+    pub pure fn from_Mat2<T:Copy NumCast>(m: &Mat2<T>) -> Mat4<T> {
+        let _0 = cast(0);
+        let _1 = cast(1);
+        Mat4::new(m[0][0], m[0][1], _0, _0,
+                  m[1][0], m[1][1], _0, _0,
+                       _0,      _0, _1, _0,
+                       _0,      _0, _0, _1)
+    }
+    
+    #[inline(always)]
+    pub pure fn from_Mat3<T:Copy NumCast>(m: &Mat3<T>) -> Mat4<T> {
+        let _0 = cast(0);
+        let _1 = cast(1);
+        Mat4::new(m[0][0], m[0][1], m[0][2], _0,
+                  m[1][0], m[1][1], m[1][2], _0,
+                  m[2][0], m[2][1], m[2][2], _0,
+                       _0,      _0,      _0, _1)
+    }
+    
+    #[inline(always)]
+    pub pure fn zero<T:Copy NumCast>() -> Mat4<T> {
+        let _0 = cast(0);
+        Mat4::new(_0, _0, _0, _0,
+                  _0, _0, _0, _0,
+                  _0, _0, _0, _0,
+                  _0, _0, _0, _0)
+    }
+    
+    #[inline(always)]
+    pub pure fn identity<T:Copy NumCast>() -> Mat4<T> {
+        let _0 = cast(0);
+        let _1 = cast(1);
+        Mat4::new(_1, _0, _0, _0,
+                  _0, _1, _0, _0,
+                  _0, _0, _1, _0,
+                  _0, _0, _0, _1)
+    }
+}
+
+pub impl<T:Copy> Mat4<T>: Matrix<T, Vec4<T>, Vec4<T>> {
+    #[inline(always)]
+    static pure fn dim() -> uint { 4 }
+    
+    #[inline(always)]
+    pure fn rows() -> uint { 4 }
+    
+    #[inline(always)]
+    pure fn cols() -> uint { 4 }
+    
+    #[inline(always)]
+    pure fn is_col_major() -> bool { true }
+    
+    #[inline(always)]
+    pure fn is_square() -> bool { true }
+    
+    #[inline(always)]
+    pure fn col(i: uint) -> Vec4<T> { self[i] }
+    
+    #[inline(always)]
+    pure fn row(i: uint) -> Vec4<T> {
+        Vec4::new(self[0][i],
+                  self[1][i],
+                  self[2][i],
+                  self[3][i])
+    }
+    
+    #[inline(always)]
+    pure fn index(i: uint) -> Vec4<T> {
+        unsafe { do buf_as_slice(
+            transmute::<*Mat4<T>, *Vec4<T>>(
+                to_unsafe_ptr(&self)), 4) |slice| { slice[i] }
+        }
+    }
+    
+    #[inline(always)]
+    pure fn to_ptr() -> *T {
+        self[0].to_ptr()
+    }
+}
+
+pub impl<T:Copy Num> Mat4<T>: NumericMatrix<T, Vec4<T>, Vec4<T>> {
+    #[inline(always)]
+    pure fn neg() -> Mat4<T> {
+        Mat4::from_cols(-self[0], -self[1], -self[2], -self[3])
+    }
+    
+    #[inline(always)]
+    pure fn mul_t(value: T) -> Mat4<T> {
+        Mat4::from_cols(self[0].mul_t(value),
+                        self[1].mul_t(value),
+                        self[2].mul_t(value),
+                        self[3].mul_t(value))
+    }
+    
+    #[inline(always)]
+    pure fn mul_v(other: &Vec4<T>) -> Vec4<T> {
+        Vec4::new(self.row(0).dot(other),
+                  self.row(1).dot(other),
+                  self.row(2).dot(other),
+                  self.row(3).dot(other))
+    }
+    
+    #[inline(always)]
+    pure fn add_m(other: &Mat4<T>) -> Mat4<T> {
+        Mat4::from_cols(self[0].add_v(&other[0]),
+                        self[1].add_v(&other[1]),
+                        self[2].add_v(&other[2]),
+                        self[3].add_v(&other[3]))
+    }
+    
+    #[inline(always)]
+    pure fn sub_m(other: &Mat4<T>) -> Mat4<T> {
+        Mat4::from_cols(self[0].sub_v(&other[0]),
+                        self[1].sub_v(&other[1]),
+                        self[2].sub_v(&other[2]),
+                        self[3].sub_v(&other[3]))
+    }
+}
+
+pub impl<T:Copy Num NumCast FuzzyEq Signed Ord> Mat4<T>: NumericMatrix_NxN<T, Vec4<T>> {
+    #[inline(always)]
+    pure fn mul_m(other: &Mat4<T>) -> Mat4<T> {
+        // Surprisingly when building with optimisation turned on this is actually
+        // faster than writing out the matrix multiplication in expanded form.
+        // If you don't believe me, see ./test/performance/matrix_mul.rs
+        Mat4::new(self.row(0).dot(&other.col(0)), self.row(1).dot(&other.col(0)), self.row(2).dot(&other.col(0)), self.row(3).dot(&other.col(0)),
+                  self.row(0).dot(&other.col(1)), self.row(1).dot(&other.col(1)), self.row(2).dot(&other.col(1)), self.row(3).dot(&other.col(1)),
+                  self.row(0).dot(&other.col(2)), self.row(1).dot(&other.col(2)), self.row(2).dot(&other.col(2)), self.row(3).dot(&other.col(2)),
+                  self.row(0).dot(&other.col(3)), self.row(1).dot(&other.col(3)), self.row(2).dot(&other.col(3)), self.row(3).dot(&other.col(3)))
+    }
+    
+    pure fn det() -> T {
+        self[0][0]*Mat3::new(self[1][1], self[2][1], self[3][1],
+                             self[1][2], self[2][2], self[3][2],
+                             self[1][3], self[2][3], self[3][3]).det() -
+        self[1][0]*Mat3::new(self[0][1], self[2][1], self[3][1],
+                             self[0][2], self[2][2], self[3][2],
+                             self[0][3], self[2][3], self[3][3]).det() +
+        self[2][0]*Mat3::new(self[0][1], self[1][1], self[3][1],
+                             self[0][2], self[1][2], self[3][2],
+                             self[0][3], self[1][3], self[3][3]).det() -
+        self[3][0]*Mat3::new(self[0][1], self[1][1], self[2][1],
+                             self[0][2], self[1][2], self[2][2],
+                             self[0][3], self[1][3], self[2][3]).det()
+    }
+
+    pure fn invert() -> Option<Mat4<T>> {
+        let d = self.det();
+        let _0 = cast(0);
+        if d.fuzzy_eq(&_0) {
+            None
+        } else {
+
+            // Gauss Jordan Elimination with partial pivoting
+
+            let mut a = self.transpose();
+            let mut inv = Mat4::identity::<T>();
+
+            // Find largest pivot column j among rows j..3
+            for uint::range(0, 4) |j| {
+                let mut i1 = j;
+                for uint::range(j + 1, 4) |i| {
+                    if abs(&a[i][j]) > abs(&a[i1][j]) {
+                        i1 = i;
+                    }
+                }
+
+                // Swap rows i1 and j in a and inv to
+                // put pivot on diagonal
+                let c = [mut a.x, a.y, a.z, a.w];
+                c[i1] <-> c[j];
+                a = Mat4::from_cols(c[0], c[1], c[2], c[3]);
+                let c = [mut inv.x, inv.y, inv.z, inv.w];
+                c[i1] <-> c[j];
+                inv = Mat4::from_cols(c[0], c[1], c[2], c[3]);
+
+                // Scale row j to have a unit diagonal
+                let c = [mut inv.x, inv.y, inv.z, inv.w];
+                c[j] = c[j].div_t(a[j][j]);
+                inv = Mat4::from_cols(c[0], c[1], c[2], c[3]);
+                let c = [mut a.x, a.y, a.z, a.w];
+                c[j] = c[j].div_t(a[j][j]);
+                a = Mat4::from_cols(c[0], c[1], c[2], c[3]);
+
+                // Eliminate off-diagonal elems in col j of a,
+                // doing identical ops to inv
+                for uint::range(0, 4) |i| {
+                    if i != j {
+                        let c = [mut inv.x, inv.y, inv.z, inv.w];
+                        c[i] = c[i].sub_v(&c[j].mul_t(a[i][j]));
+                        inv = Mat4::from_cols(c[0], c[1], c[2], c[3]);
+
+                        let c = [mut a.x, a.y, a.z, a.w];
+                        c[i] = c[i].sub_v(&c[j].mul_t(a[i][j]));
+                        a = Mat4::from_cols(c[0], c[1], c[2], c[3]); 
+                    }
+                }
+            }
+            Some(inv.transpose())
+        }
+    }
+    
+    #[inline(always)]
+    pure fn transpose() -> Mat4<T> {
+        Mat4::new(self[0][0], self[1][0], self[2][0], self[3][0],
+                  self[0][1], self[1][1], self[2][1], self[3][1],
+                  self[0][2], self[1][2], self[2][2], self[3][2],
+                  self[0][3], self[1][3], self[2][3], self[3][3])
+    }
+    
+    #[inline(always)]
+    pure fn is_identity() -> bool {
+        self.fuzzy_eq(&Mat4::identity())
+    }
+    
+    #[inline(always)]
+    pure fn is_symmetric() -> bool {
+        self[0][1].fuzzy_eq(&self[1][0]) &&
+        self[0][2].fuzzy_eq(&self[2][0]) &&
+        self[0][3].fuzzy_eq(&self[3][0]) &&
+        
+        self[1][0].fuzzy_eq(&self[0][1]) &&
+        self[1][2].fuzzy_eq(&self[2][1]) &&
+        self[1][3].fuzzy_eq(&self[3][1]) &&
+        
+        self[2][0].fuzzy_eq(&self[0][2]) &&
+        self[2][1].fuzzy_eq(&self[1][2]) &&
+        self[2][3].fuzzy_eq(&self[3][2]) &&
+        
+        self[3][0].fuzzy_eq(&self[0][3]) &&
+        self[3][1].fuzzy_eq(&self[1][3]) &&
+        self[3][2].fuzzy_eq(&self[2][3])
+    }
+    
+    #[inline(always)]
+    pure fn is_diagonal() -> bool {
+        let _0 = cast(0);
+        self[0][1].fuzzy_eq(&_0) &&
+        self[0][2].fuzzy_eq(&_0) &&
+        self[0][3].fuzzy_eq(&_0) &&
+        
+        self[1][0].fuzzy_eq(&_0) &&
+        self[1][2].fuzzy_eq(&_0) &&
+        self[1][3].fuzzy_eq(&_0) &&
+        
+        self[2][0].fuzzy_eq(&_0) &&
+        self[2][1].fuzzy_eq(&_0) &&
+        self[2][3].fuzzy_eq(&_0) &&
+        
+        self[3][0].fuzzy_eq(&_0) &&
+        self[3][1].fuzzy_eq(&_0) &&
+        self[3][2].fuzzy_eq(&_0)
+    }
+    
+    #[inline(always)]
+    pure fn is_rotated() -> bool {
+        !self.fuzzy_eq(&Mat4::identity())
+    }
+
+    #[inline(always)]
+    pure fn is_invertible() -> bool {
+        let _0 = cast(0);
+        !self.det().fuzzy_eq(&_0)
+    }
+}
+
+pub impl<T> Mat4<T>: Matrix4<T> {
+    
+}
+
+// TODO: make work for T:Integer
+pub impl<T:Copy FuzzyEq> Mat4<T>: Eq {
+    #[inline(always)]
+    pure fn eq(other: &Mat4<T>) -> bool {
+        self.fuzzy_eq(other)
+    }
+    
+    #[inline(always)]
+    pure fn ne(other: &Mat4<T>) -> bool {
+        !(self == *other)
+    }
+}
+
+pub impl<T:Copy Eq> Mat4<T>: ExactEq {
+    #[inline(always)]
+    pure fn exact_eq(other: &Mat4<T>) -> bool {
+        self[0].exact_eq(&other[0]) &&
+        self[1].exact_eq(&other[1]) &&
+        self[2].exact_eq(&other[2]) &&
+        self[3].exact_eq(&other[3])
+    }
+}
+
+pub impl<T:Copy FuzzyEq> Mat4<T>: FuzzyEq {
+    #[inline(always)]
+    pure fn fuzzy_eq(other: &Mat4<T>) -> bool {
+        self[0].fuzzy_eq(&other[0]) &&
+        self[1].fuzzy_eq(&other[1]) &&
+        self[2].fuzzy_eq(&other[2]) &&
+        self[3].fuzzy_eq(&other[3])
+    }
+}
diff --git a/src/num/ext.rs b/src/num/ext.rs
new file mode 100644
index 0000000..8232ba8
--- /dev/null
+++ b/src/num/ext.rs
@@ -0,0 +1,59 @@
+/**
+ * Various traits intended to be used with the built in numeric types. These
+ * allow one to be more specific with trait bounds when using generics.
+ *
+ * Note: These traits won't be able to be used to their full potential until
+ * trait inheritence is implemented.
+ */
+
+use core::cmp::{Eq, Ord};
+
+use std::cmp::FuzzyEq;
+
+use num::cast::*;
+use num::consts::*;
+
+trait NumExt: Copy, Eq, Num, NumCast, Ord {}
+
+trait UnSignedExt: NumExt {}
+
+pub impl u8:   UnSignedExt {}
+pub impl u16:  UnSignedExt {}
+pub impl u32:  UnSignedExt {}
+pub impl u64:  UnSignedExt {}
+pub impl uint: UnSignedExt {}
+
+
+trait SignedExt: NumExt {}
+
+pub impl i8:    SignedExt {}
+pub impl i16:   SignedExt {}
+pub impl i32:   SignedExt {}
+pub impl i64:   SignedExt {}
+pub impl int:   SignedExt {}
+
+pub impl f32:   SignedExt {}
+pub impl f64:   SignedExt {}
+pub impl float: SignedExt {}
+
+
+trait IntegerExt: NumExt, IntConsts {}
+
+pub impl u8:   IntegerExt {}
+pub impl u16:  IntegerExt {}
+pub impl u32:  IntegerExt {}
+pub impl u64:  IntegerExt {}
+pub impl uint: IntegerExt {}
+
+pub impl i8:   IntegerExt {}
+pub impl i16:  IntegerExt {}
+pub impl i32:  IntegerExt {}
+pub impl i64:  IntegerExt {}
+pub impl int:  IntegerExt {}
+
+
+trait FloatExt: NumExt, FloatConsts, FuzzyEq {}
+
+pub impl f32:   FloatExt {}
+pub impl f64:   FloatExt {}
+pub impl float: FloatExt {}
\ No newline at end of file
diff --git a/src/num/traits.rs b/src/num/traits.rs
deleted file mode 100644
index 6aabfcd..0000000
--- a/src/num/traits.rs
+++ /dev/null
@@ -1,59 +0,0 @@
-/**
- * Various traits intended to be used with the built in numeric types. These
- * allow one to be more specific with trait bounds when using generics.
- *
- * Note: These traits won't be able to be used to their full potential until
- * trait inheritence is implemented.
- */
-
-use core::cmp::{Eq, Ord};
-
-use std::cmp::FuzzyEq;
-
-use num::cast::*;
-use num::consts::*;
-
-trait NumExt: Copy, Eq, Num, NumCast, Ord {}
-
-trait UnSignedNum: NumExt {}
-
-pub impl u8:   UnSignedNum {}
-pub impl u16:  UnSignedNum {}
-pub impl u32:  UnSignedNum {}
-pub impl u64:  UnSignedNum {}
-pub impl uint: UnSignedNum {}
-
-
-trait SignedNum: NumExt {}
-
-pub impl i8:    SignedNum {}
-pub impl i16:   SignedNum {}
-pub impl i32:   SignedNum {}
-pub impl i64:   SignedNum {}
-pub impl int:   SignedNum {}
-
-pub impl f32:   SignedNum {}
-pub impl f64:   SignedNum {}
-pub impl float: SignedNum {}
-
-
-trait IntegerNum: NumExt, IntConsts {}
-
-pub impl u8:   IntegerNum {}
-pub impl u16:  IntegerNum {}
-pub impl u32:  IntegerNum {}
-pub impl u64:  IntegerNum {}
-pub impl uint: IntegerNum {}
-
-pub impl i8:   IntegerNum {}
-pub impl i16:  IntegerNum {}
-pub impl i32:  IntegerNum {}
-pub impl i64:  IntegerNum {}
-pub impl int:  IntegerNum {}
-
-
-trait FloatNum: NumExt, FloatConsts, FuzzyEq {}
-
-pub impl f32:   FloatNum {}
-pub impl f64:   FloatNum {}
-pub impl float: FloatNum {}
\ No newline at end of file
diff --git a/src/quat.rs b/src/quat.rs
new file mode 100644
index 0000000..e8ecfae
--- /dev/null
+++ b/src/quat.rs
@@ -0,0 +1,306 @@
+use core::cast::transmute;
+use core::cmp::Eq;
+use core::ptr::{addr_of, to_unsafe_ptr};
+use core::vec::raw::buf_as_slice;
+
+use std::cmp::FuzzyEq;
+
+use dim::Dimensional;
+use funs::exp::*;
+use funs::trig::*;
+use funs::common::*;
+use math::*;
+use mat::{Mat3, Mat4};
+use num::cast::*;
+use vec::Vec3;
+
+
+// These quaternion type aliases are not actually specified in the GLSL spec
+// but they follow the same nomenclature
+
+pub type quat4  = Quat<f32>;         /// a single-precision floating-point quaternion
+pub type dquat4 = Quat<f64>;         /// a double-precision floating-point quaternion
+
+//
+//  Quaternion
+//
+pub trait Quaternion<T>: Dimensional<T>, Eq, Neg<self> {
+    pure fn mul_t(value: T) -> self;
+    pure fn div_t(value: T) -> self;
+    
+    pure fn mul_v(vec: &Vec3<T>) -> Vec3<T>;
+    
+    pure fn add_q(other: &self) -> self;
+    pure fn sub_q(other: &self) -> self;
+    pure fn mul_q(other: &self) -> self;
+    
+    pure fn dot(other: &self) -> T;
+    
+    pure fn conjugate() -> self;
+    pure fn inverse() -> self;
+    pure fn length2() -> T;
+    pure fn length() -> T;
+    pure fn normalize() -> self;
+    
+    pure fn nlerp(other: &self, amount: T) -> self;
+    pure fn slerp(other: &self, amount: T) -> self;
+    
+    pure fn to_Mat3() -> Mat3<T>;
+    pure fn to_Mat4() -> Mat4<T>;
+}
+
+pub trait ToQuat<T> {
+    pure fn to_Quat() -> Quat<T>;
+}
+
+
+
+
+
+
+pub struct Quat<T> { w: T, x: T, y: T, z: T }
+
+pub mod Quat {
+    #[inline(always)]
+    pub pure fn new<T>(w: T, x: T, y: T, z: T) -> Quat<T> {
+        Quat { w: move w, x: move x, y: move y, z: move z }
+    }
+    
+    #[inline(always)]
+    pub pure fn from_sv<T:Copy>(s: T, v: Vec3<T>) -> Quat<T> {
+        Quat::new(s, v.x, v.y, v.z)
+    }
+
+    #[inline(always)]
+    pub pure fn from_axis_angle<T:Copy Num NumCast Trig AngleConv>(axis: Vec3<T>, theta: T) -> Quat<T> {
+        let half = radians(&theta) / cast(2);
+        from_sv(cos(&half), axis.mul_t(sin(&half)))
+    }
+    
+    #[inline(always)]
+    pub pure fn zero<T:Copy NumCast>() -> Quat<T> {
+        let _0 = cast(0);
+        Quat::new(_0, _0, _0, _0)
+    }
+    
+    #[inline(always)]
+    pub pure fn identity<T:Copy NumCast>() -> Quat<T> {
+        let _0 = cast(0);
+        Quat::new(cast(1), _0, _0, _0)
+    }
+}
+
+pub impl<T:Copy Num NumCast Trig Exp Extent Ord FuzzyEq> Quat<T>: Quaternion<T> {
+    #[inline(always)]
+    static pure fn dim() -> uint { 4 }
+    
+    #[inline(always)]
+    pure fn to_ptr() -> *T {
+        addr_of(&self[0])
+    }
+    
+    #[inline(always)]
+    pure fn neg() -> Quat<T> {
+        Quat::new(-self[0], -self[1], -self[2], -self[3])
+    }
+    
+    #[inline(always)]
+    pure fn mul_t(value: T) -> Quat<T> {
+        Quat::new(self[0] * value,
+                  self[1] * value,
+                  self[2] * value,
+                  self[3] * value)
+    }
+    
+    #[inline(always)]
+    pure fn div_t(value: T) -> Quat<T> {
+        Quat::new(self[0] / value,
+                  self[1] / value,
+                  self[2] / value,
+                  self[3] / value)
+    }
+
+    #[inline(always)]
+    pure fn mul_v(vec: &Vec3<T>) -> Vec3<T>  {
+        let base = Vec3{ x:self.x, y:self.y, z:self.z };
+        let tmp = base.cross(vec).add_v(&vec.mul_t(self.w));
+        base.cross(&tmp).mul_t(cast(2)).add_v(vec)
+    }
+    
+    #[inline(always)]
+    pure fn add_q(other: &Quat<T>) -> Quat<T> {
+        Quat::new(self[0] + other[0],
+                  self[1] + other[1],
+                  self[2] + other[2],
+                  self[3] + other[3])
+    }
+    
+    #[inline(always)]
+    pure fn sub_q(other: &Quat<T>) -> Quat<T> {
+        Quat::new(self[0] - other[0],
+                  self[1] - other[1],
+                  self[2] - other[2],
+                  self[3] - other[3])
+    }
+    
+    #[inline(always)]
+    pure fn mul_q(other: &Quat<T>) -> Quat<T> {
+        Quat::new(self.w * other.w - self.x * other.x - self.y * other.y - self.z * other.z,
+                  self.w * other.x + self.x * other.w + self.y * other.z - self.z * other.y, 
+                  self.w * other.y + self.y * other.w + self.z * other.x - self.x * other.z, 
+                  self.w * other.z + self.z * other.w + self.x * other.y - self.y * other.x) 
+    }
+    
+    #[inline(always)]
+    pure fn dot(other: &Quat<T>) -> T {
+        self.w * other.w +
+        self.x * other.x +
+        self.y * other.y +
+        self.z * other.z
+    }
+    
+    #[inline(always)]
+    pure fn conjugate() -> Quat<T> {
+        Quat::new(self.w, -self.x, -self.y, -self.z)
+    }
+    
+    #[inline(always)]
+    pure fn inverse() -> Quat<T> {
+        let mut n: T = cast(1);
+        n /= self.length2();
+        self.conjugate().mul_t(n)
+    }
+    
+    #[inline(always)]
+    pure fn length2() -> T {
+        self.w * self.w +
+        self.x * self.x +
+        self.y * self.y +
+        self.z * self.z
+    }
+    
+    #[inline(always)]
+    pure fn length() -> T {
+        self.length2().sqrt()
+    }
+    
+    #[inline(always)]
+    pure fn normalize() -> Quat<T> {
+        let mut n: T = cast(1);
+        n /= self.length();
+        return self.mul_t(n);
+    }
+    
+    #[inline(always)]
+    pure fn nlerp(other: &Quat<T>, amount: T) -> Quat<T> {
+        let _1: T = cast(1);
+        self.mul_t(_1 - amount).add_q(&other.mul_t(amount)).normalize()
+    }
+    
+    /**
+     * Spherical Linear Intoperlation
+     *
+     * Both quaternions should be normalized first, or else strange things will
+     * will happen...
+     *
+     * Note: The `acos` used in `slerp` is an expensive operation, so unless your
+     * quarternions a far away from each other it's generally more advisable to
+     * use nlerp when you know your rotations are going to be small.
+     *
+     * See *[Understanding Slerp, Then Not Using It]
+     * (http://number-none.com/product/Understanding%20Slerp,%20Then%20Not%20Using%20It/)*
+     * for more information. The [Arcsynthesis OpenGL tutorial]
+     * (http://www.arcsynthesis.org/gltut/Positioning/Tut08%20Interpolation.html)
+     * also provides a good explanation.
+     */
+    #[inline(always)]
+    pure fn slerp(other: &Quat<T>, amount: T) -> Quat<T> {
+        let dot: T = cast(self.dot(other));
+        
+        // if quaternions are close together use `nlerp`
+        let dot_threshold = cast(0.9995);
+        if dot > dot_threshold { return self.nlerp(other, amount) }
+        
+        let robust_dot = dot.clamp(&-cast(1), &cast(1));    // stay within the domain of acos()
+        let theta_0 = acos(&robust_dot);                    // the angle between the quaternions
+        let theta = theta_0 * amount;                       // the fraction of theta specified by `amount`
+        
+        let q = other.sub_q(&self.mul_t(robust_dot))
+                     .normalize();
+        
+        self.mul_t(cos(&theta)).add_q(&q.mul_t(sin(&theta)))
+    }
+    
+    #[inline(always)]
+    pure fn to_Mat3() -> Mat3<T> {
+        let x2 = self.x + self.x;
+        let y2 = self.y + self.y;
+        let z2 = self.z + self.z;
+        
+        let xx2 = x2 * self.x;
+        let xy2 = x2 * self.y;
+        let xz2 = x2 * self.z;
+        
+        let yy2 = y2 * self.y;
+        let yz2 = y2 * self.z;
+        let zz2 = z2 * self.z;
+        
+        let wy2 = y2 * self.w;
+        let wz2 = z2 * self.w;
+        let wx2 = x2 * self.w;
+        
+        let _1: T = cast(1);
+        
+        Mat3::new(_1 - yy2 - zz2,      xy2 - wz2,      xz2 + wy2,
+                       xy2 + wz2, _1 - xx2 - zz2,      yz2 - wx2,
+                       xz2 - wy2,      yz2 + wx2, _1 - xx2 - yy2)
+    }
+    
+    #[inline(always)]
+    pure fn to_Mat4() -> Mat4<T> {
+        self.to_Mat3().to_Mat4()
+    }
+}
+
+pub impl<T:Copy> Quat<T>: Index<uint, T> {
+    #[inline(always)]
+    pure fn index(i: uint) -> T {
+        unsafe { do buf_as_slice(
+            transmute::<*Quat<T>, *T>(
+                to_unsafe_ptr(&self)), 4) |slice| { slice[i] }
+        }
+    }
+}
+
+// TODO: make work for T:Integer
+pub impl<T:Copy FuzzyEq> Quat<T>: Eq {
+    #[inline(always)]
+    pure fn eq(other: &Quat<T>) -> bool {
+        self.fuzzy_eq(other)
+    }
+    
+    #[inline(always)]
+    pure fn ne(other: &Quat<T>) -> bool {
+        !(self == *other)
+    }
+}
+
+pub impl<T:Copy Eq> Quat<T>: ExactEq {
+    #[inline(always)]
+    pure fn exact_eq(other: &Quat<T>) -> bool {
+        self[0] == other[0] &&
+        self[1] == other[1] &&
+        self[2] == other[2] &&
+        self[3] == other[3]
+    }
+}
+
+pub impl<T:Copy FuzzyEq> Quat<T>: FuzzyEq {
+    #[inline(always)]
+    pure fn fuzzy_eq(other: &Quat<T>) -> bool {
+        self[0].fuzzy_eq(&other[0]) &&
+        self[1].fuzzy_eq(&other[1]) &&
+        self[2].fuzzy_eq(&other[2]) &&
+        self[3].fuzzy_eq(&other[3])
+    }
+}
\ No newline at end of file
diff --git a/src/test/test_mat.rs b/src/test/test_mat.rs
new file mode 100644
index 0000000..8f08fec
--- /dev/null
+++ b/src/test/test_mat.rs
@@ -0,0 +1,335 @@
+use mat::*;
+use vec::*;
+
+// TODO
+
+#[test]
+fn test_Mat2() {
+    let a = Mat2 { x: Vec2 { x: 1f, y: 3f },
+                   y: Vec2 { x: 2f, y: 4f } };
+    let b = Mat2 { x: Vec2 { x: 2f, y: 4f },
+                   y: Vec2 { x: 3f, y: 5f } };
+    let v1 = Vec2::new(1f, 2f);
+    let f1 = 0.5f;
+    
+    assert a == Mat2::new(1f, 3f,
+                          2f, 4f);
+    
+    assert a == Mat2::from_cols(Vec2::new(1f, 3f),
+                                Vec2::new(2f, 4f));
+    
+    assert Mat2::from_value(4f64) == Mat2::new(4f64, 0f64,
+                                               0f64, 4f64);
+    
+    assert a[0] == Vec2::new(1f, 3f);
+    assert a[1] == Vec2::new(2f, 4f);
+    
+    assert a.row(0) == Vec2::new(1f, 2f);
+    assert a.row(1) == Vec2::new(3f, 4f);
+    
+    assert a.col(0) == Vec2::new(1f, 3f);
+    assert a.col(1) == Vec2::new(2f, 4f);
+
+    assert a.det() == -2f;
+    
+    assert a.neg() == Mat2::new(-1f, -3f,
+                                -2f, -4f);
+    assert -a == a.neg();
+    
+    assert a.mul_t(f1) == Mat2::new(0.5f, 1.5f,
+                                    1.0f, 2.0f);
+    assert a.mul_v(&v1) == Vec2::new(5f, 11f);
+    
+    assert a.add_m(&b) == Mat2::new(3f, 7f,
+                                    5f, 9f);
+    assert a.sub_m(&b) == Mat2::new(-1f, -1f,
+                                    -1f, -1f);
+    assert a.mul_m(&b) == Mat2::new(10.0, 22.0,
+                                    13.0, 29.0);
+    
+    assert a.transpose() == Mat2::new(1f, 2f,
+                                      3f, 4f);
+
+    assert option::unwrap(a.invert()) == Mat2::new(-2f,  1.5f,
+                                                    1f, -0.5f);
+
+    assert Mat2::new(0f, 2f,
+                     0f, 5f).invert().is_none();
+    
+    // exact_eq
+    // fuzzy_eq
+    // eq
+    
+    assert Mat2::identity::<float>().is_identity();
+    assert Mat2::identity::<float>().is_symmetric();
+    assert Mat2::identity::<float>().is_diagonal();
+    assert !Mat2::identity::<float>().is_rotated();
+    assert Mat2::identity::<float>().is_invertible();
+    
+    assert !a.is_identity();
+    assert !a.is_symmetric();
+    assert !a.is_diagonal();
+    assert a.is_rotated();
+    assert a.is_invertible();
+    
+    let c = Mat2::new(2f, 1f,
+                      1f, 2f);
+    assert !c.is_identity();
+    assert c.is_symmetric();
+    assert !c.is_diagonal();
+    assert c.is_rotated();
+    assert c.is_invertible();
+    
+    assert Mat2::from_value(6f).is_diagonal();
+    
+    assert a.to_Mat3() == Mat3::new(1f, 3f, 0f,
+                                    2f, 4f, 0f,
+                                    0f, 0f, 1f);
+    
+    assert a.to_Mat4() == Mat4::new(1f, 3f, 0f, 0f,
+                                    2f, 4f, 0f, 0f,
+                                    0f, 0f, 1f, 0f,
+                                    0f, 0f, 0f, 1f);
+}
+
+#[test]
+fn test_Mat3() {
+    let a = Mat3 { x: Vec3 { x: 1f, y: 4f, z:  7f },
+                   y: Vec3 { x: 2f, y: 5f, z:  8f },
+                   z: Vec3 { x: 3f, y: 6f, z:  9f } };
+    let b = Mat3 { x: Vec3 { x: 2f, y: 5f, z:  8f },
+                   y: Vec3 { x: 3f, y: 6f, z:  9f },
+                   z: Vec3 { x: 4f, y: 7f, z: 10f } };
+    let v1 = Vec3::new(1f, 2f, 3f);
+    let f1 = 0.5f;
+    
+    assert a == Mat3::new(1f, 4f, 7f,
+                          2f, 5f, 8f,
+                          3f, 6f, 9f);
+    
+    assert a == Mat3::from_cols(Vec3::new(1f, 4f, 7f),
+                                Vec3::new(2f, 5f, 8f),
+                                Vec3::new(3f, 6f, 9f));
+    
+    assert Mat3::from_value(4f64) == Mat3::new(4f64, 0f64, 0f64,
+                                               0f64, 4f64, 0f64,
+                                               0f64, 0f64, 4f64);
+    
+    assert Mat3::from_Mat2(&Mat2::new(1f32, 3f32,
+                                      2f32, 4f32)) == Mat3::new(1f32, 3f32, 0f32,
+                                                                2f32, 4f32, 0f32,
+                                                                0f32, 0f32, 1f32);
+    
+    assert a[0] == Vec3::new(1f, 4f, 7f);
+    assert a[1] == Vec3::new(2f, 5f, 8f);
+    assert a[2] == Vec3::new(3f, 6f, 9f);
+    
+    assert a.row(0) == Vec3::new(1f, 2f, 3f);
+    assert a.row(1) == Vec3::new(4f, 5f, 6f);
+    assert a.row(2) == Vec3::new(7f, 8f, 9f);
+    
+    assert a.col(0) == Vec3::new(1f, 4f, 7f);
+    assert a.col(1) == Vec3::new(2f, 5f, 8f);
+    assert a.col(2) == Vec3::new(3f, 6f, 9f);
+
+    assert a.det() == 0f;
+    
+    assert a.neg() == Mat3::new(-1f, -4f, -7f,
+                                -2f, -5f, -8f,
+                                -3f, -6f, -9f);
+    assert -a == a.neg();
+    
+    assert a.mul_t(f1) == Mat3::new(0.5f, 2.0f, 3.5f,
+                                    1.0f, 2.5f, 4.0f,
+                                    1.5f, 3.0f, 4.5f);
+    assert a.mul_v(&v1) == Vec3::new(14f, 32f, 50f);
+    
+    assert a.add_m(&b) == Mat3::new(3f,  9f, 15f,
+                                    5f, 11f, 17f,
+                                    7f, 13f, 19f);
+    assert a.sub_m(&b) == Mat3::new(-1f, -1f, -1f,
+                                    -1f, -1f, -1f,
+                                    -1f, -1f, -1f);
+    assert a.mul_m(&b) == Mat3::new(36f,  81f, 126f,
+                                    42f,  96f, 150f,
+                                    48f, 111f, 174f);
+    
+    assert a.transpose() == Mat3::new(1f, 2f, 3f,
+                                      4f, 5f, 6f,
+                                      7f, 8f, 9f);
+
+    assert a.invert().is_none();
+
+    assert option::unwrap(Mat3::identity::<float>().invert())
+        == Mat3::identity::<float>();
+    
+    assert option::unwrap(Mat3::new(2f, 4f, 6f,
+                                    0f, 2f, 4f,
+                                    0f, 0f, 1f).invert())
+        == Mat3::new(0.5f,  -1f,  1f,
+                       0f, 0.5f, -2f,
+                       0f,   0f,  1f);
+    // exact_eq
+    // fuzzy_eq
+    // eq
+    
+    assert Mat3::identity::<float>().is_identity();
+    assert Mat3::identity::<float>().is_symmetric();
+    assert Mat3::identity::<float>().is_diagonal();
+    assert !Mat3::identity::<float>().is_rotated();
+    assert Mat3::identity::<float>().is_invertible();
+    
+    assert !a.is_identity();
+    assert !a.is_symmetric();
+    assert !a.is_diagonal();
+    assert a.is_rotated();
+    assert !a.is_invertible();
+    
+    let c = Mat3::new(3f, 2f, 1f,
+                      2f, 3f, 2f,
+                      1f, 2f, 3f);
+    assert !c.is_identity();
+    assert c.is_symmetric();
+    assert !c.is_diagonal();
+    assert c.is_rotated();
+    assert c.is_invertible();
+    
+    assert Mat3::from_value(6f).is_diagonal();
+    
+    assert a.to_Mat4() == Mat4::new(1f, 4f, 7f, 0f,
+                                    2f, 5f, 8f, 0f,
+                                    3f, 6f, 9f, 0f,
+                                    0f, 0f, 0f, 1f);
+    
+    // to_Quaternion
+}
+
+#[test]
+fn test_Mat4() {
+    let a = Mat4 { x: Vec4 { x: 1f, y: 5f, z:  9f, w: 13f },
+                   y: Vec4 { x: 2f, y: 6f, z: 10f, w: 14f },
+                   z: Vec4 { x: 3f, y: 7f, z: 11f, w: 15f },
+                   w: Vec4 { x: 4f, y: 8f, z: 12f, w: 16f } };
+    let b = Mat4 { x: Vec4 { x: 2f, y: 6f, z: 10f, w: 14f },
+                   y: Vec4 { x: 3f, y: 7f, z: 11f, w: 15f },
+                   z: Vec4 { x: 4f, y: 8f, z: 12f, w: 16f },
+                   w: Vec4 { x: 5f, y: 9f, z: 13f, w: 17f } };
+    let c = Mat4 { x: Vec4 { x: 3f, y: 2f, z:  1f, w:  1f },
+                   y: Vec4 { x: 2f, y: 3f, z:  2f, w:  2f },
+                   z: Vec4 { x: 1f, y: 2f, z:  3f, w:  3f },
+                   w: Vec4 { x: 0f, y: 1f, z:  1f, w:  0f } };
+    let v1 = Vec4::new(1f, 2f, 3f, 4f);
+    let f1 = 0.5f;
+    
+    assert a == Mat4::new(1f, 5f,  9f, 13f,
+                          2f, 6f, 10f, 14f,
+                          3f, 7f, 11f, 15f,
+                          4f, 8f, 12f, 16f);
+    
+    assert a == Mat4::from_cols(Vec4::new(1f, 5f,  9f, 13f),
+                                Vec4::new(2f, 6f, 10f, 14f),
+                                Vec4::new(3f, 7f, 11f, 15f),
+                                Vec4::new(4f, 8f, 12f, 16f));
+    
+    assert Mat4::from_value(4f64) == Mat4::new(4f64, 0f64, 0f64, 0f64,
+                                               0f64, 4f64, 0f64, 0f64,
+                                               0f64, 0f64, 4f64, 0f64,
+                                               0f64, 0f64, 0f64, 4f64);
+    
+    assert Mat4::from_Mat2(&Mat2::new(1f, 3f,
+                                      2f, 4f)) == Mat4::new(1f, 3f, 0f, 0f,
+                                                            2f, 4f, 0f, 0f,
+                                                            0f, 0f, 1f, 0f,
+                                                            0f, 0f, 0f, 1f);
+    
+    assert Mat4::from_Mat3(&Mat3::new(1f32, 4f32, 7f32,
+                                      2f32, 5f32, 8f32,
+                                      3f32, 6f32, 9f32)) == Mat4::new(1f32, 4f32, 7f32, 0f32,
+                                                                      2f32, 5f32, 8f32, 0f32,
+                                                                      3f32, 6f32, 9f32, 0f32,
+                                                                      0f32, 0f32, 0f32, 1f32);
+    
+    assert a[0] == Vec4::new(1f, 5f,  9f, 13f);
+    assert a[1] == Vec4::new(2f, 6f, 10f, 14f);
+    assert a[2] == Vec4::new(3f, 7f, 11f, 15f);
+    assert a[3] == Vec4::new(4f, 8f, 12f, 16f);
+    
+    assert a.row(0) == Vec4::new( 1f,  2f,  3f,  4f);
+    assert a.row(1) == Vec4::new( 5f,  6f,  7f,  8f);
+    assert a.row(2) == Vec4::new( 9f, 10f, 11f, 12f);
+    assert a.row(3) == Vec4::new(13f, 14f, 15f, 16f);
+    
+    assert a.col(0) == Vec4::new(1f, 5f,  9f, 13f);
+    assert a.col(1) == Vec4::new(2f, 6f, 10f, 14f);
+    assert a.col(2) == Vec4::new(3f, 7f, 11f, 15f);
+    assert a.col(3) == Vec4::new(4f, 8f, 12f, 16f);
+
+    assert a.det() == 0f;
+    
+    assert a.neg() == Mat4::new(-1f, -5f,  -9f, -13f,
+                                -2f, -6f, -10f, -14f,
+                                -3f, -7f, -11f, -15f,
+                                -4f, -8f, -12f, -16f);
+    assert -a == a.neg();
+    
+    assert a.mul_t(f1) == Mat4::new(0.5f, 2.5f, 4.5f, 6.5f,
+                                    1.0f, 3.0f, 5.0f, 7.0f,
+                                    1.5f, 3.5f, 5.5f, 7.5f,
+                                    2.0f, 4.0f, 6.0f, 8.0f);
+    assert a.mul_v(&v1) == Vec4::new(30.0, 70.0, 110.0, 150.0);
+    
+    assert a.add_m(&b) == Mat4::new(3f, 11f, 19f, 27f,
+                                    5f, 13f, 21f, 29f,
+                                    7f, 15f, 23f, 31f,
+                                    9f, 17f, 25f, 33f);
+    assert a.sub_m(&b) == Mat4::new(-1f, -1f, -1f, -1f,
+                                    -1f, -1f, -1f, -1f,
+                                    -1f, -1f, -1f, -1f,
+                                    -1f, -1f, -1f, -1f);
+    assert a.mul_m(&b) == Mat4::new(100f, 228f, 356f, 484f,
+                                    110f, 254f, 398f, 542f,
+                                    120f, 280f, 440f, 600f,
+                                    130f, 306f, 482f, 658f);
+    
+    assert a.transpose() == Mat4::new( 1f,  2f,  3f,  4f,
+                                       5f,  6f,  7f,  8f,
+                                       9f, 10f, 11f, 12f,
+                                      13f, 14f, 15f, 16f);
+
+    assert option::unwrap(c.invert())
+        == Mat4::new( 5f, -4f,  1f,  0f,
+                     -4f,  8f, -4f,  0f,
+                      4f, -8f,  4f,  8f,
+                     -3f,  4f,  1f, -8f).mul_t(0.125f);
+
+    assert option::unwrap(Mat4::identity::<float>().invert())
+        == Mat4::identity::<float>();
+    
+    // exact_eq
+    // fuzzy_eq
+    // eq
+    
+    assert Mat4::identity::<float>().is_identity();
+    assert Mat4::identity::<float>().is_symmetric();
+    assert Mat4::identity::<float>().is_diagonal();
+    assert !Mat4::identity::<float>().is_rotated();
+    assert Mat4::identity::<float>().is_invertible();
+    
+    assert !a.is_identity();
+    assert !a.is_symmetric();
+    assert !a.is_diagonal();
+    assert a.is_rotated();
+    assert !a.is_invertible();
+    
+    let c = Mat4::new(4f, 3f, 2f, 1f,
+                      3f, 4f, 3f, 2f,
+                      2f, 3f, 4f, 3f,
+                      1f, 2f, 3f, 4f);
+    assert !c.is_identity();
+    assert c.is_symmetric();
+    assert !c.is_diagonal();
+    assert c.is_rotated();
+    assert c.is_invertible();
+    
+    assert Mat4::from_value(6f).is_diagonal();
+}
diff --git a/src/test/test_quat.rs b/src/test/test_quat.rs
new file mode 100644
index 0000000..f859d68
--- /dev/null
+++ b/src/test/test_quat.rs
@@ -0,0 +1,29 @@
+use mat::*;
+use quat::*;
+use vec::*;
+
+// TODO
+
+#[test]
+fn test_Quat() {
+    let a = Quat { w: 1f, x: 2f, y: 3f, z: 4f };
+    // let b = Quat { data: [ 5f, 6f, 7f, 8f ] };
+    // let f1 = 1.5f;
+    // let f2 = 0.5f;
+    
+    assert a == Quat::new(1f, 2f, 3f, 4f);
+    
+    assert Quat::zero()     == Quat::new(0f, 0f, 0f, 0f);
+    assert Quat::identity() == Quat::new(1f, 0f, 0f, 0f);
+    
+    assert a.w == 1f;
+    assert a.x == 2f;
+    assert a.y == 3f;
+    assert a.z == 4f;
+    assert a[0] == 1f;
+    assert a[1] == 2f;
+    assert a[2] == 3f;
+    assert a[3] == 4f;
+    
+    // TODO
+}
\ No newline at end of file
diff --git a/src/test/test_vec.rs b/src/test/test_vec.rs
new file mode 100644
index 0000000..db7950f
--- /dev/null
+++ b/src/test/test_vec.rs
@@ -0,0 +1,200 @@
+use std::cmp::FuzzyEq;
+use vec::*;
+
+// TODO
+
+#[test]
+fn test_Vec2() {
+    // assert Vec2::dim == 2;
+    
+    let a = Vec2 { x: 1f, y: 2f };
+    let b = Vec2 { x: 3f, y: 4f };
+    let f1 = 1.5f;
+    let f2 = 0.5f;
+    
+    assert Vec2::new(1f, 2f) == a;
+    assert Vec2::from_value(1f32) == Vec2::new(1f32, 1f32);
+    
+    assert Vec2::zero()     == Vec2::new(0f, 0f);
+    assert Vec2::unit_x()   == Vec2::new(1f, 0f);
+    assert Vec2::unit_y()   == Vec2::new(0f, 1f);
+    assert Vec2::identity() == Vec2::new(1f, 1f);
+    
+    assert a.x == 1f;
+    assert a.y == 2f;
+    assert a[0] == 1f;
+    assert a[1] == 2f;
+    
+    assert -a      == Vec2::new(-1f, -2f);
+    assert a.neg() == Vec2::new(-1f, -2f);
+    
+    assert a.mul_t(f1) == Vec2::new( 1.5f, 3.0f);
+    assert a.div_t(f2) == Vec2::new( 2.0f, 4.0f);
+    
+    assert a.add_v(&b) == Vec2::new( 4f,  6f);
+    assert a.sub_v(&b) == Vec2::new(-2f, -2f);
+    
+    // exact_eq
+    // fuzzy_eq
+    // eq
+    
+    // assert c.abs()       == Vec2::new( 2.0f,  1.0f);
+    // assert c.min(&d)      == Vec2::new(-2.0f, -1.0f);
+    // assert c.max(&d)      == Vec2::new( 1.0f,  0.0f);
+}
+
+#[test]
+fn test_Vec2_geometric() {
+    let a = Vec2::new(5f, 12f); // (5, 12, 13) Pythagorean triple
+    let b0 = Vec2::new(3f, 4f); // (3, 4, 5) Pythagorean triple
+    let b = a.add_v(&b0);
+    
+    assert a.length() == 13f;
+    assert a.length2() == 13f * 13f;
+    
+    assert b0.length() == 5f;
+    assert b0.length2() == 5f * 5f;
+    
+    assert a.distance(&b) == 5f;
+    assert a.distance2(&b) == 5f * 5f;
+    
+    let c = Vec2::new(-2.0f, -1.0f);
+    let d = Vec2::new( 1.0f,  0.0f);
+    
+    assert c.lerp(&d, 0.75f) == Vec2::new(0.250f, -0.250f);
+}
+
+#[test]
+fn test_Vec3() {
+    // assert Vec3::dim == 3;
+    
+    let a = Vec3 { x: 1f, y: 2f, z: 3f };
+    let b = Vec3 { x: 4f, y: 5f, z: 6f };
+    let f1 = 1.5f;
+    let f2 = 0.5f;
+    
+    assert Vec3::new(1f, 2f, 3f) == a;
+    assert Vec3::from_value(1f32) == Vec3::new(1f32, 1f32, 1f32);
+    
+    assert Vec3::zero()     == Vec3::new(0f, 0f, 0f);
+    assert Vec3::unit_x()   == Vec3::new(1f, 0f, 0f);
+    assert Vec3::unit_y()   == Vec3::new(0f, 1f, 0f);
+    assert Vec3::unit_z()   == Vec3::new(0f, 0f, 1f);
+    assert Vec3::identity() == Vec3::new(1f, 1f, 1f);
+    
+    assert a.x == 1f;
+    assert a.y == 2f;
+    assert a.z == 3f;
+    assert a[0] == 1f;
+    assert a[1] == 2f;
+    assert a[2] == 3f;
+    
+    assert a.cross(&b) == Vec3::new(-3f, 6f, -3f);
+    
+    assert -a      == Vec3::new(-1f, -2f, -3f);
+    assert a.neg() == Vec3::new(-1f, -2f, -3f);
+    
+    assert a.mul_t(f1) == Vec3::new( 1.5f, 3.0f, 4.5f);
+    assert a.div_t(f2) == Vec3::new( 2.0f, 4.0f, 6.0f);
+    
+    assert a.add_v(&b) == Vec3::new( 5f,  7f,  9f);
+    assert a.sub_v(&b) == Vec3::new(-3f, -3f, -3f);
+    
+    // exact_eq
+    // fuzzy_eq
+    // eq
+    
+    // assert c.abs()        == Vec3::new( 2.0f,  1.0f, 1.0f);
+    // assert c.min(&d)      == Vec3::new(-2.0f, -1.0f, 0.5f);
+    // assert c.max(&d)      == Vec3::new( 1.0f,  0.0f, 1.0f);
+}
+
+#[test]
+fn test_Vec3_geometric() {
+    let a = Vec3::new(2f, 3f, 6f); // (2, 3, 6, 7) Pythagorean quadruple
+    let b0 = Vec3::new(1f, 4f, 8f); // (1, 4, 8, 9) Pythagorean quadruple
+    let b = a.add_v(&b0);
+    
+    assert a.length() == 7f;
+    assert a.length2() == 7f * 7f;
+    
+    assert b0.length() == 9f;
+    assert b0.length2() == 9f * 9f;
+    
+    assert a.distance(&b) == 9f;
+    assert a.distance2(&b) == 9f * 9f;
+    
+    let c = Vec3::new(-2.0f, -1.0f, 1.0f);
+    let d = Vec3::new( 1.0f,  0.0f, 0.5f);
+    
+    assert c.lerp(&d, 0.75f) == Vec3::new(0.250f, -0.250f, 0.625f);
+}
+
+#[test]
+fn test_Vec4() {
+    // assert Vec4::dim == 4;
+    
+    let a = Vec4 { x: 1f, y: 2f, z: 3f, w: 4f };
+    let b = Vec4 { x: 5f, y: 6f, z: 7f, w: 8f };
+    let f1 = 1.5f;
+    let f2 = 0.5f;
+    
+    assert Vec4::new(1f, 2f, 3f, 4f) == a;
+    assert Vec4::from_value(1f32) == Vec4::new(1f32, 1f32, 1f32, 1f32);
+    
+    assert Vec4::zero()     == Vec4::new(0f, 0f, 0f, 0f);
+    assert Vec4::unit_x()   == Vec4::new(1f, 0f, 0f, 0f);
+    assert Vec4::unit_y()   == Vec4::new(0f, 1f, 0f, 0f);
+    assert Vec4::unit_z()   == Vec4::new(0f, 0f, 1f, 0f);
+    assert Vec4::unit_w()   == Vec4::new(0f, 0f, 0f, 1f);
+    assert Vec4::identity() == Vec4::new(1f, 1f, 1f, 1f);
+    
+    assert a.x == 1f;
+    assert a.y == 2f;
+    assert a.z == 3f;
+    assert a.w == 4f;
+    assert a[0] == 1f;
+    assert a[1] == 2f;
+    assert a[2] == 3f;
+    assert a[3] == 4f;
+    
+    assert -a      == Vec4::new(-1f, -2f, -3f, -4f);
+    assert a.neg() == Vec4::new(-1f, -2f, -3f, -4f);
+    
+    assert a.mul_t(f1) == Vec4::new( 1.5f, 3.0f, 4.5f, 6.0f);
+    assert a.div_t(f2) == Vec4::new( 2.0f, 4.0f, 6.0f, 8.0f);
+    
+    assert a.add_v(&b) == Vec4::new( 6f,  8f, 10f, 12f);
+    assert a.sub_v(&b) == Vec4::new(-4f, -4f, -4f, -4f);
+    
+    assert a.dot(&b) == 70f;
+    
+    // exact_eq
+    // fuzzy_eq
+    // eq
+    
+    // assert c.abs()        == Vec4::new( 2.0f,  1.0f, 1.0f, 2.0f);
+    // assert c.min(&d)      == Vec4::new(-2.0f, -1.0f, 0.5f, 1.0f);
+    // assert c.max(&d)      == Vec4::new( 1.0f,  0.0f, 1.0f, 2.0f);
+}
+
+#[test]
+fn test_Vec4_geometric() {
+    let a = Vec4::new(1f, 2f, 4f, 10f); // (1, 2, 4, 10, 11) Pythagorean quintuple
+    let b0 = Vec4::new(1f, 2f, 8f, 10f); // (1, 2, 8, 10, 13) Pythagorean quintuple
+    let b = a.add_v(&b0);
+    
+    assert a.length() == 11f;
+    assert a.length2() == 11f * 11f;
+    
+    assert b0.length() == 13f;
+    assert b0.length2() == 13f * 13f;
+    
+    assert a.distance(&b) == 13f;
+    assert a.distance2(&b) == 13f * 13f;
+
+    let c = Vec4::new(-2.0f, -1.0f, 1.0f, 2.0f);
+    let d = Vec4::new( 1.0f,  0.0f, 0.5f, 1.0f);
+    
+    assert c.lerp(&d, 0.75f) == Vec4::new(0.250f, -0.250f, 0.625f, 1.250f);
+}
\ No newline at end of file
diff --git a/src/vec.rs b/src/vec.rs
new file mode 100644
index 0000000..8e640f2
--- /dev/null
+++ b/src/vec.rs
@@ -0,0 +1,622 @@
+use core::cast::transmute;
+use core::cmp::Eq;
+use core::ptr::{addr_of, to_unsafe_ptr};
+use core::vec::raw::buf_as_slice;
+
+use std::cmp::FuzzyEq;
+
+use dim::Dimensional;
+use funs::exp::Exp;
+use math::*;
+use num::cast::*;
+
+// GLSL equivalent type aliases
+
+pub type vec2  = Vec2<f32>;         /// a two-component single-precision floating-point vector
+pub type vec3  = Vec3<f32>;         /// a three-component single-precision floating-point vector
+pub type vec4  = Vec4<f32>;         /// a four-component single-precision floating-point vector
+
+pub type dvec2 = Vec2<f64>;         /// a two-component double-precision floating-point vector
+pub type dvec3 = Vec3<f64>;         /// a three-component double-precision floating-point vector
+pub type dvec4 = Vec4<f64>;         /// a four-component double-precision floating-point vector
+
+pub type bvec2 = Vec2<bool>;        /// a two-component Boolean vector
+pub type bvec3 = Vec3<bool>;        /// a three-component Boolean vector
+pub type bvec4 = Vec4<bool>;        /// a four-component Boolean vector
+
+pub type ivec2 = Vec2<i32>;         /// a two-component signed integer vector
+pub type ivec3 = Vec3<i32>;         /// a three-component signed integer vector
+pub type ivec4 = Vec4<i32>;         /// a four-component signed integer vector
+
+pub type uvec2 = Vec2<u32>;         /// a two-component unsigned integer vector
+pub type uvec3 = Vec3<u32>;         /// a three-component unsigned integer vector
+pub type uvec4 = Vec4<u32>;         /// a four-component unsigned integer vector
+
+
+pub trait Vector<T>: Dimensional<T>, Eq {}
+
+pub trait NumericVector<T>: Vector<T>, Neg<self>{
+    pure fn mul_t(value: T) -> self;
+    pure fn div_t(value: T) -> self;
+    
+    pure fn add_v(other: &self) -> self;
+    pure fn sub_v(other: &self) -> self;
+    
+    pure fn dot(other: &self) -> T;
+}
+
+pub trait GeometricVector<T>: NumericVector<T> {
+    pure fn length2() -> T;
+    pure fn length() -> T;
+    pure fn distance2(other: &self) -> T;
+    pure fn distance(other: &self) -> T;
+    pure fn normalize() -> self;
+    pure fn lerp(other: &self, amount: T) -> self;
+}
+
+pub trait Vector2<T>: Vector<T> {
+    // static pure fn new(x: T, y: T) -> self;
+    // static pure fn from_value(value: T) -> self;
+}
+
+pub trait Vector3<T>: Vector<T> {
+    // static pure fn new(x: T, y: T, z: T) -> self;
+    // static pure fn from_value(value: T) -> self;
+    
+    pure fn cross(other: &self) -> self;
+}
+
+pub trait Vector4<T>: Vector<T> {
+    // pub static pure fn new(x: T, y: T, z: T, w: T) -> self;
+    // pub static pure fn from_value(value: T) -> self;
+}
+
+
+
+
+
+//
+//  Vec2
+//
+pub struct Vec2<T> { x: T, y: T }
+
+pub mod Vec2 {
+    
+    #[inline(always)]
+    pub pure fn new<T>(x: T, y: T) -> Vec2<T> {
+        Vec2 { x: move x, y: move y }
+    }
+    
+    #[inline(always)]
+    pub pure fn from_value<T:Copy>(value: T) -> Vec2<T> {
+        Vec2::new(value, value)
+    }
+    
+    #[inline(always)]
+    pub pure fn zero<T:Copy NumCast>() -> Vec2<T> {
+        let _0 = cast(0);
+        Vec2::new(_0, _0)
+    }
+    
+    #[inline(always)]
+    pub pure fn unit_x<T:Copy NumCast>() -> Vec2<T> {
+        let _0 = cast(0);
+        let _1 = cast(1);
+        Vec2::new(_1, _0)
+    }
+    
+    #[inline(always)]
+    pub pure fn unit_y<T:Copy NumCast>() -> Vec2<T> {
+        let _0 = cast(0);
+        let _1 = cast(1);
+        Vec2::new(_0, _1)
+    }
+    
+    #[inline(always)]
+    pub pure fn identity<T:Copy NumCast>() -> Vec2<T> {
+        let _1 = cast(1);
+        Vec2::new(_1, _1)
+    }
+}
+
+pub impl<T:Copy> Vec2<T>: Vector<T> {
+    #[inline(always)]
+    static pure fn dim() -> uint { 2 }
+    
+    #[inline(always)]
+    pure fn index(i: uint) -> T {
+        unsafe { do buf_as_slice(
+            transmute::<*Vec2<T>, *T>(
+                to_unsafe_ptr(&self)), 2) |slice| { slice[i] }
+        }
+    }
+    
+    #[inline(always)]
+    pure fn to_ptr() -> *T {
+        ptr::addr_of(&self[0])
+    }
+}
+    
+pub impl<T:Copy Num> Vec2<T>: NumericVector<T> {
+    #[inline(always)]
+    pure fn neg() -> Vec2<T> {
+        Vec2::new(-self[0], -self[1])
+    }
+    
+    #[inline(always)]
+    pure fn mul_t(value: T) -> Vec2<T> {
+        Vec2::new(self[0] * value,
+                  self[1] * value)
+    }
+    
+    #[inline(always)]
+    pure fn div_t(value: T) -> Vec2<T> {
+        Vec2::new(self[0] / value,
+                  self[1] / value)
+    }
+    
+    #[inline(always)]
+    pure fn add_v(other: &Vec2<T>) -> Vec2<T> {
+        Vec2::new(self[0] + other[0],
+                  self[1] + other[1])
+    }
+    
+    #[inline(always)]
+    pure fn sub_v(other: &Vec2<T>) -> Vec2<T> {
+        Vec2::new(self[0] - other[0],
+                  self[1] - other[1])
+    }
+    
+    #[inline(always)]
+    pure fn dot(other: &Vec2<T>) -> T {
+        self[0] * other[0] +
+        self[1] * other[1]
+    }
+}
+    
+pub impl<T:Copy Num NumCast Exp> Vec2<T>: GeometricVector<T> {
+    #[inline(always)]
+    pure fn length2() -> T {
+        self.dot(&self)
+    }
+    
+    #[inline(always)]
+    pure fn length() -> T {
+        self.length2().sqrt()
+    }
+    
+    #[inline(always)]
+    pure fn distance2(other: &Vec2<T>) -> T {
+        other.sub_v(&self).length2()
+    }
+    
+    #[inline(always)]
+    pure fn distance(other: &Vec2<T>) -> T {
+        other.distance2(&self).sqrt()
+    }
+    
+    #[inline(always)]
+    pure fn normalize() -> Vec2<T> {
+        let mut n: T = cast(1); 
+        n /= self.length();
+        return self.mul_t(n);
+    }
+    
+    #[inline(always)]
+    pure fn lerp(other: &Vec2<T>, amount: T) -> Vec2<T> {
+        self.add_v(&other.sub_v(&self).mul_t(amount))
+    }
+}
+
+// TODO: make work for T:Integer
+pub impl<T:Copy FuzzyEq> Vec2<T>: Eq {
+    #[inline(always)]
+    pure fn eq(other: &Vec2<T>) -> bool {
+        self.fuzzy_eq(other)
+    }
+    
+    #[inline(always)]
+    pure fn ne(other: &Vec2<T>) -> bool {
+        !(self == *other)
+    }
+}
+
+impl<T:Copy Eq> Vec2<T>: ExactEq {
+    #[inline(always)]
+    pure fn exact_eq(other: &Vec2<T>) -> bool {
+        self[0] == other[0] &&
+        self[1] == other[1]
+    }
+}
+
+pub impl<T:Copy FuzzyEq> Vec2<T>: FuzzyEq {
+    #[inline(always)]
+    pure fn fuzzy_eq(other: &Vec2<T>) -> bool {
+        self[0].fuzzy_eq(&other[0]) &&
+        self[1].fuzzy_eq(&other[1])
+    }
+}
+
+
+
+
+
+
+//
+//  Vec3
+//
+pub struct Vec3<T> { x: T, y: T, z: T }
+
+pub mod Vec3 {
+    
+    #[inline(always)]
+    pub pure fn new<T>(x: T, y: T, z: T) -> Vec3<T> {
+        Vec3 { x: move x, y: move y, z: move z }
+    }
+    
+    #[inline(always)]
+    pub pure fn from_value<T:Copy>(value: T) -> Vec3<T> {
+        Vec3::new(value, value, value)
+    }
+    
+    #[inline(always)]
+    pub pure fn zero<T:Copy NumCast>() -> Vec3<T> {
+        let _0 = cast(0);
+        Vec3::new(_0, _0, _0)
+    }
+    
+    #[inline(always)]
+    pub pure fn unit_x<T:Copy NumCast>() -> Vec3<T> {
+        let _0 = cast(0);
+        let _1 = cast(1);
+        Vec3::new(_1, _0, _0)
+    }
+    
+    #[inline(always)]
+    pub pure fn unit_y<T:Copy NumCast>() -> Vec3<T> {
+        let _0 = cast(0);
+        let _1 = cast(1);
+        Vec3::new(_0, _1, _0)
+    }
+    
+    #[inline(always)]
+    pub pure fn unit_z<T:Copy NumCast>() -> Vec3<T> {
+        let _0 = cast(0);
+        let _1 = cast(1);
+        Vec3::new(_0, _0, _1)
+    }
+    
+    #[inline(always)]
+    pub pure fn identity<T:Copy NumCast>() -> Vec3<T> {
+        let _1 = cast(1);
+        Vec3::new(_1, _1, _1)
+    }
+}
+
+pub impl<T:Copy Num> Vec3<T>: Vector3<T> {
+    #[inline(always)]
+    pure fn cross(other: &Vec3<T>) -> Vec3<T> {
+        Vec3::new((self[1] * other[2]) - (self[2] * other[1]),
+                  (self[2] * other[0]) - (self[0] * other[2]),
+                  (self[0] * other[1]) - (self[1] * other[0]))
+    }
+}
+
+pub impl<T:Copy> Vec3<T>: Vector<T> {
+    #[inline(always)]
+    static pure fn dim() -> uint { 3 }
+    
+    #[inline(always)]
+    pure fn index(i: uint) -> T {
+        unsafe { do buf_as_slice(
+            transmute::<*Vec3<T>, *T>(
+                to_unsafe_ptr(&self)), 3) |slice| { slice[i] }
+        }
+    }
+    
+    #[inline(always)]
+    pure fn to_ptr() -> *T {
+        addr_of(&self[0])
+    }
+}
+
+pub impl<T:Copy Num> Vec3<T>: NumericVector<T> {
+    #[inline(always)]
+    pure fn neg() -> Vec3<T> {
+        Vec3::new(-self[0], -self[1], -self[2])
+    }
+    
+    #[inline(always)]
+    pure fn mul_t(value: T) -> Vec3<T> {
+        Vec3::new(self[0] * value,
+                  self[1] * value,
+                  self[2] * value)
+    }
+    
+    #[inline(always)]
+    pure fn div_t(value: T) -> Vec3<T> {
+        Vec3::new(self[0] / value,
+                  self[1] / value,
+                  self[2] / value)
+    }
+    
+    #[inline(always)]
+    pure fn add_v(other: &Vec3<T>) -> Vec3<T>{
+        Vec3::new(self[0] + other[0],
+                  self[1] + other[1],
+                  self[2] + other[2])
+    }
+    
+    #[inline(always)]
+    pure fn sub_v(other: &Vec3<T>) -> Vec3<T>{
+        Vec3::new(self[0] - other[0],
+                  self[1] - other[1],
+                  self[2] - other[2])
+    }
+    
+    #[inline(always)]
+    pure fn dot(other: &Vec3<T>) -> T {
+        self[0] * other[0] +
+        self[1] * other[1] +
+        self[2] * other[2]
+    }
+}
+
+pub impl<T:Copy Num NumCast Exp> Vec3<T>: GeometricVector<T> {
+    #[inline(always)]
+    pure fn length2() -> T {
+        self.dot(&self)
+    }
+    
+    #[inline(always)]
+    pure fn length() -> T {
+        self.length2().sqrt()
+    }
+    
+    #[inline(always)]
+    pure fn distance2(other: &Vec3<T>) -> T {
+        other.sub_v(&self).length2()
+    }
+    
+    #[inline(always)]
+    pure fn distance(other: &Vec3<T>) -> T {
+        other.distance2(&self).sqrt()
+    }
+    
+    #[inline(always)]
+    pure fn normalize() -> Vec3<T> {
+        let mut n: T = cast(1);
+        n /= self.length();
+        return self.mul_t(n);
+    }
+    
+    #[inline(always)]
+    pure fn lerp(other: &Vec3<T>, amount: T) -> Vec3<T> {
+        self.add_v(&other.sub_v(&self).mul_t(amount))
+    }
+}
+
+// TODO: make work for T:Integer
+pub impl<T:Copy FuzzyEq> Vec3<T>: Eq {
+    #[inline(always)]
+    pure fn eq(other: &Vec3<T>) -> bool {
+        self.fuzzy_eq(other)
+    }
+    
+    #[inline(always)]
+    pure fn ne(other: &Vec3<T>) -> bool {
+        !(self == *other)
+    }
+}
+
+pub impl<T:Copy Eq> Vec3<T>: ExactEq {
+    #[inline(always)]
+    pure fn exact_eq(other: &Vec3<T>) -> bool {
+        self[0] == other[0] &&
+        self[1] == other[1] &&
+        self[2] == other[2]
+    }
+}
+
+pub impl<T:Copy FuzzyEq> Vec3<T>: FuzzyEq {
+    #[inline(always)]
+    pure fn fuzzy_eq(other: &Vec3<T>) -> bool {
+        self[0].fuzzy_eq(&other[0]) &&
+        self[1].fuzzy_eq(&other[1]) &&
+        self[2].fuzzy_eq(&other[2])
+    }
+}
+
+
+
+
+
+
+//
+//  Vec4
+//
+pub struct Vec4<T> { x: T, y: T, z: T, w: T }
+
+pub mod Vec4 {
+    #[inline(always)]
+    pub pure fn new<T>(x: T, y: T, z: T, w: T) -> Vec4<T> {
+        Vec4 { x: move x, y: move y, z: move z, w: move w }
+    }
+    
+    #[inline(always)]
+    pub pure fn from_value<T:Copy>(value: T) -> Vec4<T> {
+        Vec4::new(value, value, value, value)
+    }
+    
+    #[inline(always)]
+    pub pure fn zero<T:Copy NumCast>() -> Vec4<T> {
+        let _0 = cast(0);
+        Vec4::new(_0, _0, _0, _0)
+    }
+    
+    #[inline(always)]
+    pub pure fn unit_x<T:Copy NumCast>() -> Vec4<T> {
+        let _0 = cast(0);
+        let _1 = cast(1);
+        Vec4::new(_1, _0, _0, _0)
+    }
+    
+    #[inline(always)]
+    pub pure fn unit_y<T:Copy NumCast>() -> Vec4<T> {
+        let _0 = cast(0);
+        let _1 = cast(1);
+        Vec4::new(_0, _1, _0, _0)
+    }
+    
+    #[inline(always)]
+    pub pure fn unit_z<T:Copy NumCast>() -> Vec4<T> {
+        let _0 = cast(0);
+        let _1 = cast(1);
+        Vec4::new(_0, _0, _1, _0)
+    }
+    
+    #[inline(always)]
+    pub pure fn unit_w<T:Copy NumCast>() -> Vec4<T> {
+        let _0 = cast(0);
+        let _1 = cast(1);
+        Vec4::new(_0, _0, _0, _1)
+    }
+    
+    #[inline(always)]
+    pub pure fn identity<T:Copy NumCast>() -> Vec4<T> {
+        let _1 = cast(1);
+        Vec4::new(_1, _1, _1, _1)
+    }
+}
+
+pub impl<T:Copy> Vec4<T>: Vector<T> {
+    #[inline(always)]
+    static pure fn dim() -> uint { 4 }
+    
+    #[inline(always)]
+    pure fn index(i: uint) -> T {
+        unsafe { do buf_as_slice(
+            transmute::<*Vec4<T>, *T>(
+                to_unsafe_ptr(&self)), 4) |slice| { slice[i] }
+        }
+    }
+    
+    #[inline(always)]
+    pure fn to_ptr() -> *T {
+        addr_of(&self[0])
+    }
+}
+
+pub impl<T:Copy Num> Vec4<T>: NumericVector<T> {
+    #[inline(always)]
+    pure fn neg() -> Vec4<T> {
+        Vec4::new(-self[0], -self[1], -self[2], -self[3])
+    }
+    
+    #[inline(always)]
+    pure fn mul_t(value: T) -> Vec4<T> {
+        Vec4::new(self[0] * value,
+                  self[1] * value,
+                  self[2] * value,
+                  self[3] * value)
+    }
+    
+    #[inline(always)]
+    pure fn div_t(value: T) -> Vec4<T> {
+        Vec4::new(self[0] / value,
+                  self[1] / value,
+                  self[2] / value,
+                  self[3] / value)
+    }
+    
+    #[inline(always)]
+    pure fn add_v(other: &Vec4<T>) -> Vec4<T> {
+        Vec4::new(self[0] + other[0],
+                  self[1] + other[1],
+                  self[2] + other[2],
+                  self[3] + other[3])
+    }
+    
+    #[inline(always)]
+    pure fn sub_v(other: &Vec4<T>) -> Vec4<T> {
+        Vec4::new(self[0] - other[0],
+                  self[1] - other[1],
+                  self[2] - other[2],
+                  self[3] - other[3])
+    }
+    
+    #[inline(always)]
+    pure fn dot(other: &Vec4<T>) -> T {
+        self[0] * other[0] +
+        self[1] * other[1] +
+        self[2] * other[2] +
+        self[3] * other[3]
+    }
+}
+
+pub impl<T:Copy Num NumCast Exp> Vec4<T>: GeometricVector<T> {
+    #[inline(always)]
+    pure fn length2() -> T {
+        self.dot(&self)
+    }
+    
+    #[inline(always)]
+    pure fn length() -> T {
+        self.length2().sqrt()
+    }
+    
+    #[inline(always)]
+    pure fn distance2(other: &Vec4<T>) -> T {
+        other.sub_v(&self).length2()
+    }
+    
+    #[inline(always)]
+    pure fn distance(other: &Vec4<T>) -> T {
+        other.distance2(&self).sqrt()
+    }
+    
+    #[inline(always)]
+    pure fn normalize() -> Vec4<T> {
+        let mut n: T = cast(1);
+        n /= self.length();
+        return self.mul_t(n);
+    }
+    
+    #[inline(always)]
+    pure fn lerp(other: &Vec4<T>, amount: T) -> Vec4<T> {
+        self.add_v(&other.sub_v(&self).mul_t(amount))
+    }
+}
+
+pub impl<T:Copy FuzzyEq> Vec4<T>: Eq {
+    #[inline(always)]
+    pure fn eq(other: &Vec4<T>) -> bool {
+        self.fuzzy_eq(other)
+    }
+    
+    #[inline(always)]
+    pure fn ne(other: &Vec4<T>) -> bool {
+        !(self == *other)
+    }
+}
+
+// TODO: make work for T:Integer
+pub impl<T:Copy Eq> Vec4<T>: ExactEq {
+    #[inline(always)]
+    pure fn exact_eq(other: &Vec4<T>) -> bool {
+        self[0] == other[0] &&
+        self[1] == other[1] &&
+        self[2] == other[2] &&
+        self[3] == other[3]
+    }
+}
+
+pub impl<T:Copy FuzzyEq> Vec4<T>: FuzzyEq {
+    #[inline(always)]
+    pure fn fuzzy_eq(other: &Vec4<T>) -> bool {
+        self[0].fuzzy_eq(&other[0]) &&
+        self[1].fuzzy_eq(&other[1]) &&
+        self[2].fuzzy_eq(&other[2]) &&
+        self[3].fuzzy_eq(&other[3])
+    }
+}
\ No newline at end of file