Doc + slerp + conversions.

Cargo.toml

@@ -19,7 +19,6 @@ path = "src/lib.rs"
 arbitrary = [ "quickcheck" ]
 
 [dependencies]
-rustc-serialize = "0.3"
 typenum         = "1.4"
 generic-array   = "0.2"
 rand            = "0.3"
@@ -27,8 +26,13 @@ num-traits      = "0.1"
 num-complex     = "0.1"
 approx          = "0.1"
 alga            = "0.4"
+serde           = "0.9"
+serde_derive    = "0.9"
 # clippy = "*"
 
 [dependencies.quickcheck]
 optional = true
 version  = "0.3"
+
+[dev-dependencies]
+serde_json = "0.9"

Makefile

@@ -1,11 +1,11 @@
 all:
-	cargo build --features "arbitrary"
+	CARGO_INCREMENTAL=1 cargo build --features "arbitrary"
 
 doc:
-	cargo doc
+	CARGO_INCREMENTAL=1 cargo doc
 
 bench:
 	cargo bench
 
 test:
-	cargo test --features "arbitrary"
+	CARGO_INCREMENTAL=1 cargo test --features "arbitrary"

README.md

@@ -1,3 +1,6 @@
+<p align="center">
+  <img src="http://nalgebra.org/img/logo_nalgebra.svg" alt="crates.io">
+</p>
 <p align="center">
     <a href="https://crates.io/crates/nalgebra">
          <img src="http://meritbadge.herokuapp.com/nalgebra?style=flat-square" alt="crates.io">
@@ -8,75 +11,6 @@
 </p>
 <p align = "center">
     <strong>
-        <a href="http://nalgebra.org/doc/nalgebra">Documentation</a> | <a href="http://users.nphysics.org">Forum</a>
+        <a href="http://nalgebra.org">Users guide</a> | <a href="http://nalgebra.org/rustdoc/nalgebra/index.html">Documentation</a> | <a href="http://users.nphysics.org">Forum</a>
     </strong>
 </p>
-
-nalgebra
-========
-
-**nalgebra** is a low-dimensional linear algebra library written for Rust targeting:
-
-* General-purpose linear algebra (still lacks a lot of features…)
-* Real time computer graphics.
-* Real time computer physics.
-
-## Using **nalgebra**
-You will need the last stable build of the [rust compiler](http://www.rust-lang.org)
-and the official package manager: [cargo](https://github.com/rust-lang/cargo).
-
-Simply add the following to your `Cargo.toml` file:
-
-```.ignore
-[dependencies]
-nalgebra = "0.10.*"
-```
-
-
-All the functionality of **nalgebra** is grouped in one place: the root module `nalgebra::`.  This
-module re-exports everything and includes free functions for all traits methods performing
-out-of-place operations.
-
-Thus, you can import the whole prelude using:
-
-```.ignore
-use nalgebra::*;
-```
-
-However, the recommended way to use **nalgebra** is to import types and traits
-explicitly, and call free-functions using the `na::` prefix:
-
-```.rust
-extern crate nalgebra as na;
-use na::{Vector3, Rotation3, Rotation};
-
-fn main() {
-    let     a = Vector3::new(1.0f64, 1.0, 1.0);
-    let mut b = Rotation3::new(na::zero());
-
-    b.append_rotation_mut(&a);
-
-    assert!(na::approx_eq(&na::rotation(&b), &a));
-}
-```
-
-
-## Features
-**nalgebra** is meant to be a general-purpose, low-dimensional, linear algebra library, with
-an optimized set of tools for computer graphics and physics. Those features include:
-
-* Vectors with predefined static sizes: `Vector1`, `Vector2`, `Vector3`, `Vector4`, `Vector5`, `Vector6`.
-* Vector with a user-defined static size: `VectorN` (available only with the `generic_sizes` feature).
-* Points with static sizes: `Point1`, `Point2`, `Point3`, `Point4`, `Point5`, `Point6`.
-* Square matrices with static sizes: `Matrix1`, `Matrix2`, `Matrix3`, `Matrix4`, `Matrix5`, `Matrix6 `.
-* Rotation matrices: `Rotation2`, `Rotation3`
-* Quaternions: `Quaternion`, `Unit<Quaternion>`.
-* Unit-sized values (unit vectors, unit quaternions, etc.): `Unit<T>`, e.g., `Unit<Vector3<f32>>`.
-* Isometries (translation ⨯ rotation): `Isometry2`, `Isometry3`
-* Similarity transformations (translation ⨯ rotation ⨯ uniform scale): `Similarity2`, `Similarity3`.
-* 3D projections for computer graphics: `Persp3`, `PerspMatrix3`, `Ortho3`, `OrthoMatrix3`.
-* Dynamically sized heap-allocated vector: `DVector`.
-* Dynamically sized stack-allocated vectors with a maximum size: `DVector1` to `DVector6`.
-* Dynamically sized heap-allocated (square or rectangular) matrix: `DMatrix`.
-* Linear algebra and data analysis operators: `Covariance`, `Mean`, `qr`, `cholesky`.
-* Almost one trait per functionality: useful for generic programming.

examples/dimensional_genericity.rs

@@ -0,0 +1,66 @@
+extern crate alga;
+extern crate nalgebra as na;
+
+use alga::general::Real;
+use alga::linear::FiniteDimInnerSpace;
+use na::{Unit, ColumnVector, OwnedColumnVector, Vector2, Vector3};
+use na::storage::Storage;
+use na::dimension::{DimName, U1};
+
+/// Reflects a vector wrt. the hyperplane with normal `plane_normal`.
+fn reflect_wrt_hyperplane_with_algebraic_genericity<V>(plane_normal: &Unit<V>, vector: &V) -> V
+    where V: FiniteDimInnerSpace + Copy {
+    let n = plane_normal.as_ref(); // Get the underlying vector of type `V`.
+    *vector - *n * (n.dot(vector) * na::convert(2.0))
+}
+
+
+/// Reflects a vector wrt. the hyperplane with normal `plane_normal`.
+fn reflect_wrt_hyperplane_with_structural_genericity<N, D, S>(plane_normal: &Unit<ColumnVector<N, D, S>>,
+                                                               vector:       &ColumnVector<N, D, S>)
+                                                               -> OwnedColumnVector<N, D, S::Alloc>
+    where N: Real,
+          D: DimName,
+          S: Storage<N, D, U1> {
+    let n = plane_normal.as_ref(); // Get the underlying V.
+    vector - n * (n.dot(vector) * na::convert(2.0))
+}
+
+/// Reflects a 2D vector wrt. the 2D line with normal `plane_normal`.
+fn reflect_wrt_hyperplane2<N>(plane_normal: &Unit<Vector2<N>>,
+                              vector:       &Vector2<N>)
+                              -> Vector2<N>
+    where N: Real {
+    let n = plane_normal.as_ref(); // Get the underlying Vector2
+    vector - n * (n.dot(vector) * na::convert(2.0))
+}
+
+/// Reflects a 3D vector wrt. the 3D plane with normal `plane_normal`.
+/// /!\ This is an exact replicate of `reflect_wrt_hyperplane2, but for 3D.
+fn reflect_wrt_hyperplane3<N>(plane_normal: &Unit<Vector3<N>>,
+                              vector:       &Vector3<N>)
+                              -> Vector3<N>
+    where N: Real {
+    let n = plane_normal.as_ref(); // Get the underlying Vector3
+    vector - n * (n.dot(vector) * na::convert(2.0))
+}
+
+
+fn main() {
+    let plane2 = Vector2::y_axis(); // 2D plane normal.
+    let plane3 = Vector3::y_axis(); // 3D plane normal.
+
+    let v2 = Vector2::new(1.0, 2.0);      // 2D vector to be reflected.
+    let v3 = Vector3::new(1.0, 2.0, 3.0); // 3D vector to be reflected.
+
+    // We can call the same function for 2D and 3D.
+    assert_eq!(reflect_wrt_hyperplane_with_algebraic_genericity(&plane2, &v2).y, -2.0);
+    assert_eq!(reflect_wrt_hyperplane_with_algebraic_genericity(&plane3, &v3).y, -2.0);
+
+    assert_eq!(reflect_wrt_hyperplane_with_structural_genericity(&plane2, &v2).y, -2.0);
+    assert_eq!(reflect_wrt_hyperplane_with_structural_genericity(&plane3, &v3).y, -2.0);
+
+    // Call each specific implementation depending on the dimension.
+    assert_eq!(reflect_wrt_hyperplane2(&plane2, &v2).y, -2.0);
+    assert_eq!(reflect_wrt_hyperplane3(&plane3, &v3).y, -2.0);
+}

examples/homogeneous_coordinates.rs

@@ -0,0 +1,45 @@
+#[macro_use]
+extern crate approx;
+extern crate nalgebra as na;
+
+use std::f32;
+use na::{Vector2, Point2, Isometry2};
+
+
+fn use_dedicated_types() {
+    let iso = Isometry2::new(Vector2::new(1.0, 1.0), f32::consts::PI);
+    let pt  = Point2::new(1.0, 0.0);
+    let vec = Vector2::x();
+
+    let transformed_pt  = iso * pt;
+    let transformed_vec = iso * vec;
+
+    assert_relative_eq!(transformed_pt, Point2::new(0.0, 1.0));
+    assert_relative_eq!(transformed_vec, Vector2::new(-1.0, 0.0));
+}
+
+fn use_homogeneous_coordinates() {
+    let iso = Isometry2::new(Vector2::new(1.0, 1.0), f32::consts::PI);
+    let pt  = Point2::new(1.0, 0.0);
+    let vec = Vector2::x();
+
+    // Compute using homogeneous coordinates.
+    let hom_iso = iso.to_homogeneous();
+    let hom_pt  = pt.to_homogeneous();
+    let hom_vec = vec.to_homogeneous();
+
+    let hom_transformed_pt  = hom_iso * hom_pt;
+    let hom_transformed_vec = hom_iso * hom_vec;
+
+    // Convert back to the cartesian coordinates.
+    let transformed_pt  = Point2::from_homogeneous(hom_transformed_pt).unwrap();
+    let transformed_vec = Vector2::from_homogeneous(hom_transformed_vec).unwrap();
+
+    assert_relative_eq!(transformed_pt, Point2::new(0.0, 1.0));
+    assert_relative_eq!(transformed_vec, Vector2::new(-1.0, 0.0));
+}
+
+fn main() {
+    use_dedicated_types();
+    use_homogeneous_coordinates();
+}

examples/identity.rs

@@ -0,0 +1,42 @@
+extern crate alga;
+extern crate nalgebra as na;
+
+
+use alga::linear::Transformation;
+use na::{Id, Vector3, Point3, Isometry3};
+
+/*
+ * Applies `n` times the transformation `t` to the vector `v` and sum each
+ * intermediate value.
+ */
+fn complicated_algorithm<T>(v: &Vector3<f32>, t: &T, n: usize) -> Vector3<f32>
+  where T: Transformation<Point3<f32>> {
+
+  let mut result = *v;
+
+  // Do lots of operations involving t.
+  for _ in 0 .. n {
+    result = v + t.transform_vector(&result);
+  }
+
+  result
+}
+
+
+/*
+ * The two following calls are equivalent in term of result.
+ */
+fn main() {
+  let v = Vector3::new(1.0, 2.0, 3.0);
+
+  // The specialization generated by the compiler will do vector additions only.
+  let result1 = complicated_algorithm(&v, &Id::new(), 100000);
+
+  // The specialization generated by the compiler will also include matrix multiplications.
+  let iso     = Isometry3::identity();
+  let result2 = complicated_algorithm(&v, &iso, 100000);
+
+  // They both return the same result.
+  assert!(result1 == Vector3::new(100001.0, 200002.0, 300003.0));
+  assert!(result2 == Vector3::new(100001.0, 200002.0, 300003.0));
+}

examples/matrix_construction.rs

@@ -0,0 +1,62 @@
+extern crate nalgebra as na;
+
+use na::{Vector2, RowVector3, Matrix2x3, DMatrix};
+
+
+fn main() {
+    // All the following matrices are equal but constructed in different ways.
+    let m = Matrix2x3::new(1.1, 1.2, 1.3,
+                           2.1, 2.2, 2.3);
+
+    let m1 = Matrix2x3::from_rows(&[
+        RowVector3::new(1.1, 1.2, 1.3),
+        RowVector3::new(2.1, 2.2, 2.3)
+    ]);
+
+    let m2 = Matrix2x3::from_columns(&[
+        Vector2::new(1.1, 2.1),
+        Vector2::new(1.2, 2.2),
+        Vector2::new(1.3, 2.3)
+    ]);
+
+    let m3 = Matrix2x3::from_row_slice(&[
+        1.1, 1.2, 1.3,
+        2.1, 2.2, 2.3
+    ]);
+
+    let m4 = Matrix2x3::from_column_slice(&[
+        1.1, 2.1,
+        1.2, 2.2,
+        1.3, 2.3
+    ]);
+
+    let m5 = Matrix2x3::from_fn(|r, c| (r + 1) as f32 + (c + 1) as f32 / 10.0);
+
+    let m6 = Matrix2x3::from_iterator([ 1.1f32, 2.1, 1.2, 2.2, 1.3, 2.3 ].iter().cloned());
+
+    assert_eq!(m, m1); assert_eq!(m, m2); assert_eq!(m, m3);
+    assert_eq!(m, m4); assert_eq!(m, m5); assert_eq!(m, m6);
+
+    // All the following matrices are equal but constructed in different ways.
+    // This time, we used a dynamically-sized matrix to show the extra arguments
+    // for the matrix shape.
+    let dm = DMatrix::from_row_slice(4, 3, &[
+        1.0, 0.0, 0.0,
+        0.0, 1.0, 0.0,
+        0.0, 0.0, 1.0,
+        0.0, 0.0, 0.0
+    ]);
+
+    let dm1 = DMatrix::from_diagonal_element(4, 3, 1.0);
+    let dm2 = DMatrix::identity(4, 3);
+    let dm3 = DMatrix::from_fn(4, 3, |r, c| if r == c { 1.0 } else { 0.0 });
+    let dm4 = DMatrix::from_iterator(4, 3, [
+        // Components listed column-by-column.
+        1.0, 0.0, 0.0, 0.0,
+        0.0, 1.0, 0.0, 0.0,
+        0.0, 0.0, 1.0, 0.0
+    ].iter().cloned());
+
+    assert_eq!(dm, dm1); assert_eq!(dm, dm2);
+    assert_eq!(dm, dm3); assert_eq!(dm, dm4);
+}

examples/mvp.rs

@@ -0,0 +1,28 @@
+#![allow(unused_variables)]
+
+extern crate nalgebra as na;
+
+use na::{Vector3, Point3, Isometry3, Perspective3};
+
+fn main() {
+    // Our object is translated along the x axis.
+    let model = Isometry3::new(Vector3::x(), na::zero());
+
+    // Our camera looks toward the point (1.0, 0.0, 0.0).
+    // It is located at (0.0, 0.0, 1.0).
+    let eye    = Point3::new(0.0, 0.0, 1.0);
+    let target = Point3::new(1.0, 0.0, 0.0);
+    let view   = Isometry3::look_at_rh(&eye, &target, &Vector3::y());
+
+    // A perspective projection.
+    let projection = Perspective3::new(16.0 / 9.0, 3.14 / 2.0, 1.0, 1000.0);
+
+    // The combination of the model with the view is still an isometry.
+    let model_view = model * view;
+
+    // Convert everything to a `Matrix4` so that they can be combined.
+    let mat_model_view = model_view.to_homogeneous();
+
+    // Combine everything.
+    let model_view_projection = projection.as_matrix() * mat_model_view;
+}

examples/raw_pointer.rs

@@ -0,0 +1,35 @@
+extern crate nalgebra as na;
+
+use na::{Vector3, Point3, Matrix3};
+
+fn main() {
+    let v = Vector3::new(1.0f32, 0.0, 1.0);
+    let p = Point3::new(1.0f32, 0.0, 1.0);
+    let m = na::one::<Matrix3<f32>>();
+
+    // Convert to arrays.
+    let v_array = v.as_slice();
+    let p_array = p.coords.as_slice();
+    let m_array = m.as_slice();
+
+    // Get data pointers.
+    let v_pointer = v_array.as_ptr();
+    let p_pointer = p_array.as_ptr();
+    let m_pointer = m_array.as_ptr();
+
+    /* Then pass the raw pointers to some graphics API. */
+
+    unsafe {
+        assert_eq!(*v_pointer, 1.0);
+        assert_eq!(*v_pointer.offset(1), 0.0);
+        assert_eq!(*v_pointer.offset(2), 1.0);
+
+        assert_eq!(*p_pointer, 1.0);
+        assert_eq!(*p_pointer.offset(1), 0.0);
+        assert_eq!(*p_pointer.offset(2), 1.0);
+
+        assert_eq!(*m_pointer, 1.0);
+        assert_eq!(*m_pointer.offset(4), 1.0);
+        assert_eq!(*m_pointer.offset(8), 1.0);
+    }
+}

examples/scalar_genericity.rs

@@ -0,0 +1,34 @@
+extern crate alga;
+extern crate nalgebra as na;
+
+use alga::general::{RingCommutative, Real};
+use na::{Vector3, Scalar};
+
+fn print_vector<N: Scalar>(m: &Vector3<N>) {
+  println!("{:?}", m)
+}
+
+fn print_squared_norm<N: Scalar + RingCommutative>(v: &Vector3<N>) {
+  // NOTE: alternatively, nalgebra already defines `v.squared_norm()`.
+  let sqnorm = v.dot(v);
+
+  println!("{:?}", sqnorm);
+}
+
+fn print_norm<N: Real>(v: &Vector3<N>) {
+  // NOTE: alternatively, nalgebra already defines `v.norm()`.
+  let norm = v.dot(v).sqrt();
+
+  // The Real bound implies that N is Display so we can
+  // use "{}" instead of "{:?}" for the format string.
+  println!("{}", norm)
+}
+
+fn main() {
+  let v1 = Vector3::new(1, 2, 3);
+  let v2 = Vector3::new(1.0, 2.0, 3.0);
+
+  print_vector(&v1);
+  print_squared_norm(&v1);
+  print_norm(&v2);
+}

examples/screen_to_view_coords.rs

@@ -0,0 +1,24 @@
+#![allow(unused_variables)]
+
+extern crate nalgebra as na;
+
+use na::{Point2, Point3, Perspective3, Unit};
+
+
+fn main() {
+    let projection   = Perspective3::new(800.0 / 600.0, 3.14 / 2.0, 1.0, 1000.0);
+    let screen_point = Point2::new(10.0f32, 20.0);
+
+    // Compute two points in clip-space.
+    // "ndc" = normalized device coordinates.
+    let near_ndc_point = Point3::new(screen_point.x / 800.0, screen_point.y / 600.0, -1.0);
+    let far_ndc_point  = Point3::new(screen_point.x / 800.0, screen_point.y / 600.0,  1.0);
+
+    // Unproject them to view-space.
+    let near_view_point = projection.unproject_point(&near_ndc_point);
+    let far_view_point  = projection.unproject_point(&far_ndc_point);
+
+    // Compute the view-space line parameters.
+    let line_location  = near_view_point;
+    let line_direction = Unit::new_normalize(far_view_point - near_view_point);
+}

examples/transform_conversion.rs

@@ -0,0 +1,23 @@
+extern crate nalgebra as na;
+
+use na::{Vector2, Isometry2, Similarity2};
+
+fn main() {
+    // Isometry -> Similarity conversion always succeeds.
+    let iso = Isometry2::new(Vector2::new(1.0f32, 2.0), na::zero());
+    let _: Similarity2<f32> = na::convert(iso);
+
+    // Similarity -> Isometry conversion fails if the scaling factor is not 1.0.
+    let sim_without_scaling = Similarity2::new(Vector2::new(1.0f32, 2.0), 3.14, 1.0);
+    let sim_with_scaling    = Similarity2::new(Vector2::new(1.0f32, 2.0), 3.14, 2.0);
+
+    let iso_success: Option<Isometry2<f32>> = na::try_convert(sim_without_scaling);
+    let iso_fail:    Option<Isometry2<f32>> = na::try_convert(sim_with_scaling);
+
+    assert!(iso_success.is_some());
+    assert!(iso_fail.is_none());
+
+    // Similarity -> Isometry conversion can be forced at your own risks.
+    let iso_forced: Isometry2<f32> = unsafe { na::convert_unchecked(sim_with_scaling) };
+    assert_eq!(iso_success.unwrap(), iso_forced);
+}

examples/transform_matrix4.rs

@@ -0,0 +1,39 @@
+#[macro_use]
+extern crate approx;
+extern crate alga;
+extern crate nalgebra as na;
+
+use alga::linear::Transformation;
+use na::{Vector3, Point3, Matrix4};
+
+
+fn main() {
+    // Create a uniform scaling matrix with scaling factor 2.
+    let mut m =  Matrix4::new_scaling(2.0);
+
+    assert_eq!(m.transform_vector(&Vector3::x()), Vector3::x() * 2.0);
+    assert_eq!(m.transform_vector(&Vector3::y()), Vector3::y() * 2.0);
+    assert_eq!(m.transform_vector(&Vector3::z()), Vector3::z() * 2.0);
+
+    // Append a nonuniform scaling in-place.
+    m.append_nonuniform_scaling_mut(&Vector3::new(1.0, 2.0, 3.0));
+
+    assert_eq!(m.transform_vector(&Vector3::x()), Vector3::x() * 2.0);
+    assert_eq!(m.transform_vector(&Vector3::y()), Vector3::y() * 4.0);
+    assert_eq!(m.transform_vector(&Vector3::z()), Vector3::z() * 6.0);
+
+    // Append a translation out-of-place.
+    let m2 = m.append_translation(&Vector3::new(42.0, 0.0, 0.0));
+
+    assert_eq!(m2.transform_point(&Point3::new(1.0, 1.0, 1.0)), Point3::new(42.0 + 2.0, 4.0, 6.0));
+
+    // Create rotation.
+    let rot        = Matrix4::from_scaled_axis(&Vector3::x() * 3.14);
+    let rot_then_m = m * rot; // Right-multiplication is equivalent to prepending `rot` to `m`.
+    let m_then_rot = rot * m; // Left-multiplication is equivalent to appending `rot` to `m`.
+
+    let pt = Point3::new(1.0, 2.0, 3.0);
+
+    assert_relative_eq!(m.transform_point(&rot.transform_point(&pt)), rot_then_m.transform_point(&pt));
+    assert_relative_eq!(rot.transform_point(&m.transform_point(&pt)), m_then_rot.transform_point(&pt));
+}

examples/transform_vector_point.rs

@@ -0,0 +1,20 @@
+#[macro_use]
+extern crate approx;
+extern crate nalgebra as na;
+
+use std::f32;
+use na::{Vector2, Point2, Isometry2};
+
+fn main() {
+    let t = Isometry2::new(Vector2::new(1.0, 1.0), f32::consts::PI);
+    let p = Point2::new(1.0, 0.0); // Will be affected by te rotation and the translation.
+    let v = Vector2::x();          // Will *not* be affected by the translation.
+
+    assert_relative_eq!(t * p, Point2::new(-1.0 + 1.0, 1.0));
+    //                                     ^^^^ │ ^^^^^^^^
+    //                                  rotated │ translated
+
+    assert_relative_eq!(t * v, Vector2::new(-1.0, 0.0));
+    //                                      ^^^^^
+    //                                   rotated only
+}

examples/transform_vector_point3.rs

@@ -0,0 +1,27 @@
+extern crate alga;
+extern crate nalgebra as na;
+
+use alga::linear::Transformation;
+use na::{Vector3, Vector4, Point3, Matrix4};
+
+
+fn main() {
+    let mut m = Matrix4::new_rotation_wrt_point(Vector3::x() * 1.57, Point3::new(1.0, 2.0, 1.0));
+    m.append_scaling_mut(2.0);
+
+    let point1             = Point3::new(2.0, 3.0, 4.0);
+    let homogeneous_point2 = Vector4::new(2.0, 3.0, 4.0, 1.0);
+
+    // First option: use the dedicated `.transform_point(...)` method.
+    let transformed_point1 = m.transform_point(&point1);
+    // Second option: use the homogeneous coordinates of the point.
+    let transformed_homogeneous_point2 = m * homogeneous_point2;
+
+    // Recover the 3D point from its 4D homogeneous coordinates.
+    let transformed_point2 = Point3::from_homogeneous(transformed_homogeneous_point2);
+
+    // Check that transforming the 3D point with the `.transform_point` method is
+    // indeed equivalent to multiplying its 4D homogeneous coordinates by the 4x4
+    // matrix.
+    assert_eq!(transformed_point1, transformed_point2.unwrap());
+}

examples/transformation_pointer.rs

@@ -0,0 +1,25 @@
+extern crate nalgebra as na;
+
+use na::{Vector3, Isometry3};
+
+fn main() {
+    let iso = Isometry3::new(Vector3::new(1.0f32, 0.0, 1.0), na::zero());
+
+    // Compute the homogeneous coordinates first.
+    let iso_matrix  = iso.to_homogeneous();
+    let iso_array   = iso_matrix.as_slice();
+    let iso_pointer = iso_array.as_ptr();
+
+    /* Then pass the raw pointer to some graphics API. */
+
+    unsafe {
+        assert_eq!(*iso_pointer,            1.0);
+        assert_eq!(*iso_pointer.offset(5),  1.0);
+        assert_eq!(*iso_pointer.offset(10), 1.0);
+        assert_eq!(*iso_pointer.offset(15), 1.0);
+
+        assert_eq!(*iso_pointer.offset(12), 1.0);
+        assert_eq!(*iso_pointer.offset(13), 0.0);
+        assert_eq!(*iso_pointer.offset(14), 1.0);
+    }
+}

examples/unit_wrapper.rs

@@ -0,0 +1,30 @@
+extern crate nalgebra as na;
+
+use na::{Unit, Vector3};
+
+fn length_on_direction_with_unit(v: &Vector3<f32>, dir: &Unit<Vector3<f32>>) -> f32 {
+  // No need to normalize `dir`: we know that it is non-zero and normalized.
+  v.dot(dir.as_ref())
+}
+
+
+
+fn length_on_direction_without_unit(v: &Vector3<f32>, dir: &Vector3<f32>) -> f32 {
+  // Obligatory normalization of the direction vector (and test, for robustness).
+  if let Some(unit_dir) = dir.try_normalize(1.0e-6) {
+    v.dot(&unit_dir)
+  }
+  else {
+    // Normalization failed because the norm was too small.
+    panic!("Invalid input direction.")
+  }
+}
+
+fn main() {
+    let v = Vector3::new(1.0, 2.0, 3.0);
+
+    let l1 = length_on_direction_with_unit(&v, &Vector3::y_axis());
+    let l2 = length_on_direction_without_unit(&v, &Vector3::y());
+
+    assert_eq!(l1, l2)
+}

src/core/alias.rs

@@ -19,52 +19,88 @@ pub type MatrixNM<N, R, C> = Matrix<N, R, C, MatrixArray<N, R, C>>;
 /// A staticaly sized column-major square matrix with `D` rows and columns.
 pub type MatrixN<N, D> = MatrixNM<N, D, D>;
 
+/// A stack-allocated, column-major, 1x1 square matrix.
 pub type Matrix1<N> = MatrixN<N, U1>;
+/// A stack-allocated, column-major, 2x2 square matrix.
 pub type Matrix2<N> = MatrixN<N, U2>;
+/// A stack-allocated, column-major, 3x3 square matrix.
 pub type Matrix3<N> = MatrixN<N, U3>;
+/// A stack-allocated, column-major, 4x4 square matrix.
 pub type Matrix4<N> = MatrixN<N, U4>;
+/// A stack-allocated, column-major, 5x5 square matrix.
 pub type Matrix5<N> = MatrixN<N, U5>;
+/// A stack-allocated, column-major, 6x6 square matrix.
 pub type Matrix6<N> = MatrixN<N, U6>;
 
+/// A stack-allocated, column-major, 1x2 square matrix.
 pub type Matrix1x2<N> = MatrixNM<N, U1, U2>;
+/// A stack-allocated, column-major, 1x3 square matrix.
 pub type Matrix1x3<N> = MatrixNM<N, U1, U3>;
+/// A stack-allocated, column-major, 1x4 square matrix.
 pub type Matrix1x4<N> = MatrixNM<N, U1, U4>;
+/// A stack-allocated, column-major, 1x5 square matrix.
 pub type Matrix1x5<N> = MatrixNM<N, U1, U5>;
+/// A stack-allocated, column-major, 1x6 square matrix.
 pub type Matrix1x6<N> = MatrixNM<N, U1, U6>;
 
+/// A stack-allocated, column-major, 2x3 square matrix.
 pub type Matrix2x3<N> = MatrixNM<N, U2, U3>;
+/// A stack-allocated, column-major, 2x4 square matrix.
 pub type Matrix2x4<N> = MatrixNM<N, U2, U4>;
+/// A stack-allocated, column-major, 2x5 square matrix.
 pub type Matrix2x5<N> = MatrixNM<N, U2, U5>;
+/// A stack-allocated, column-major, 2x6 square matrix.
 pub type Matrix2x6<N> = MatrixNM<N, U2, U6>;
 
+/// A stack-allocated, column-major, 3x4 square matrix.
 pub type Matrix3x4<N> = MatrixNM<N, U3, U4>;
+/// A stack-allocated, column-major, 3x5 square matrix.
 pub type Matrix3x5<N> = MatrixNM<N, U3, U5>;
+/// A stack-allocated, column-major, 3x6 square matrix.
 pub type Matrix3x6<N> = MatrixNM<N, U3, U6>;
 
+/// A stack-allocated, column-major, 4x5 square matrix.
 pub type Matrix4x5<N> = MatrixNM<N, U4, U5>;
+/// A stack-allocated, column-major, 4x6 square matrix.
 pub type Matrix4x6<N> = MatrixNM<N, U4, U6>;
 
+/// A stack-allocated, column-major, 5x6 square matrix.
 pub type Matrix5x6<N> = MatrixNM<N, U5, U6>;
 
 
+/// A stack-allocated, column-major, 2x1 square matrix.
 pub type Matrix2x1<N> = MatrixNM<N, U2, U1>;
+/// A stack-allocated, column-major, 3x1 square matrix.
 pub type Matrix3x1<N> = MatrixNM<N, U3, U1>;
+/// A stack-allocated, column-major, 4x1 square matrix.
 pub type Matrix4x1<N> = MatrixNM<N, U4, U1>;
+/// A stack-allocated, column-major, 5x1 square matrix.
 pub type Matrix5x1<N> = MatrixNM<N, U5, U1>;
+/// A stack-allocated, column-major, 6x1 square matrix.
 pub type Matrix6x1<N> = MatrixNM<N, U6, U1>;
 
+/// A stack-allocated, column-major, 3x2 square matrix.
 pub type Matrix3x2<N> = MatrixNM<N, U3, U2>;
+/// A stack-allocated, column-major, 4x2 square matrix.
 pub type Matrix4x2<N> = MatrixNM<N, U4, U2>;
+/// A stack-allocated, column-major, 5x2 square matrix.
 pub type Matrix5x2<N> = MatrixNM<N, U5, U2>;
+/// A stack-allocated, column-major, 6x2 square matrix.
 pub type Matrix6x2<N> = MatrixNM<N, U6, U2>;
 
+/// A stack-allocated, column-major, 4x3 square matrix.
 pub type Matrix4x3<N> = MatrixNM<N, U4, U3>;
+/// A stack-allocated, column-major, 5x3 square matrix.
 pub type Matrix5x3<N> = MatrixNM<N, U5, U3>;
+/// A stack-allocated, column-major, 6x3 square matrix.
 pub type Matrix6x3<N> = MatrixNM<N, U6, U3>;
 
+/// A stack-allocated, column-major, 5x4 square matrix.
 pub type Matrix5x4<N> = MatrixNM<N, U5, U4>;
+/// A stack-allocated, column-major, 6x4 square matrix.
 pub type Matrix6x4<N> = MatrixNM<N, U6, U4>;
 
+/// A stack-allocated, column-major, 6x5 square matrix.
 pub type Matrix6x5<N> = MatrixNM<N, U6, U5>;
 
 
@@ -81,11 +117,17 @@ pub type DVector<N> = Matrix<N, Dynamic, U1, MatrixVec<N, Dynamic, U1>>;
 /// A statically sized D-dimensional column vector.
 pub type VectorN<N, D> = MatrixNM<N, D, U1>;
 
+/// A stack-allocated, 1-dimensional column vector.
 pub type Vector1<N> = VectorN<N, U1>;
+/// A stack-allocated, 2-dimensional column vector.
 pub type Vector2<N> = VectorN<N, U2>;
+/// A stack-allocated, 3-dimensional column vector.
 pub type Vector3<N> = VectorN<N, U3>;
+/// A stack-allocated, 4-dimensional column vector.
 pub type Vector4<N> = VectorN<N, U4>;
+/// A stack-allocated, 5-dimensional column vector.
 pub type Vector5<N> = VectorN<N, U5>;
+/// A stack-allocated, 6-dimensional column vector.
 pub type Vector6<N> = VectorN<N, U6>;
 
 
@@ -102,9 +144,15 @@ pub type RowDVector<N> = Matrix<N, U1, Dynamic, MatrixVec<N, U1, Dynamic>>;
 /// A statically sized D-dimensional row vector.
 pub type RowVectorN<N, D> = MatrixNM<N, U1, D>;
 
+/// A stack-allocated, 1-dimensional row vector.
 pub type RowVector1<N> = RowVectorN<N, U1>;
+/// A stack-allocated, 2-dimensional row vector.
 pub type RowVector2<N> = RowVectorN<N, U2>;
+/// A stack-allocated, 3-dimensional row vector.
 pub type RowVector3<N> = RowVectorN<N, U3>;
+/// A stack-allocated, 4-dimensional row vector.
 pub type RowVector4<N> = RowVectorN<N, U4>;
+/// A stack-allocated, 5-dimensional row vector.
 pub type RowVector5<N> = RowVectorN<N, U5>;
+/// A stack-allocated, 6-dimensional row vector.
 pub type RowVector6<N> = RowVectorN<N, U6>;

src/core/allocator.rs

@@ -1,3 +1,5 @@
+//! Abstract definition of a matrix data storage allocator.
+
 use std::any::Any;
 
 use core::Scalar;
@@ -15,6 +17,7 @@ use core::storage::{Storage, OwnedStorage};
 /// Every allocator must be both static and dynamic. Though not all implementations may share the
 /// same `Buffer` type.
 pub trait Allocator<N: Scalar, R: Dim, C: Dim>: Any + Sized {
+    /// The type of buffer this allocator can instanciate.
     type Buffer: OwnedStorage<N, R, C, Alloc = Self>;
 
     /// Allocates a buffer with the given number of rows and columns without initializing its content.

src/core/cg.rs

@@ -25,7 +25,7 @@ impl<N, D: DimName, S> SquareMatrix<N, D, S>
     #[inline]
     pub fn new_scaling(scaling: N) -> Self {
         let mut res = Self::from_diagonal_element(scaling);
-        res[(D::dim(), D::dim())] = N::one();
+        res[(D::dim() - 1, D::dim() - 1)] = N::one();
 
         res
     }
@@ -50,7 +50,7 @@ impl<N, D: DimName, S> SquareMatrix<N, D, S>
               SB: Storage<N, DimNameDiff<D, U1>, U1>,
               S::Alloc: Allocator<N, DimNameDiff<D, U1>, U1> {
         let mut res = Self::one();
-        res.fixed_slice_mut::<DimNameDiff<D, U1>, U1>(0, D::dim()).copy_from(translation);
+        res.fixed_slice_mut::<DimNameDiff<D, U1>, U1>(0, D::dim() - 1).copy_from(translation);
 
         res
     }
@@ -311,7 +311,7 @@ impl<N, D: DimName, S> SquareMatrix<N, D, S>
               S::Alloc: Allocator<N, DimNameDiff<D, U1>, U1> {
         for i in 0 .. D::dim() {
             for j in 0 .. D::dim() - 1 {
-                self[(j, i)] += shift[i] * self[(D::dim(), j)];
+                self[(j, i)] += shift[j] * self[(D::dim() - 1, i)];
             }
         }
     }
@@ -324,12 +324,12 @@ impl<N, D: DimName, S> SquareMatrix<N, D, S>
               S::Alloc: Allocator<N, DimNameDiff<D, U1>, U1> +
                         Allocator<N, DimNameDiff<D, U1>, DimNameDiff<D, U1>> +
                         Allocator<N, U1, DimNameDiff<D, U1>> {
-        let scale = self.fixed_slice::<U1, DimNameDiff<D, U1>>(D::dim(), 0).tr_dot(&shift);
+        let scale = self.fixed_slice::<U1, DimNameDiff<D, U1>>(D::dim() - 1, 0).tr_dot(&shift);
         let post_translation = self.fixed_slice::<DimNameDiff<D, U1>, DimNameDiff<D, U1>>(0, 0) * shift;
 
-        self[(D::dim(), D::dim())] += scale;
+        self[(D::dim() - 1, D::dim() - 1)] += scale;
 
-        let mut translation = self.fixed_slice_mut::<DimNameDiff<D, U1>, U1>(0, D::dim());
+        let mut translation = self.fixed_slice_mut::<DimNameDiff<D, U1>, U1>(0, D::dim() - 1);
         translation += post_translation;
     }
 }
@@ -349,7 +349,7 @@ impl<N, D, SA, SB> Transformation<PointBase<N, DimNameDiff<D, U1>, SB>> for Squa
     fn transform_vector(&self, v: &ColumnVector<N, DimNameDiff<D, U1>, SB>)
         -> ColumnVector<N, DimNameDiff<D, U1>, SB> {
         let transform  = self.fixed_slice::<DimNameDiff<D, U1>, DimNameDiff<D, U1>>(0, 0);
-        let normalizer = self.fixed_slice::<U1, DimNameDiff<D, U1>>(D::dim(), 0);
+        let normalizer = self.fixed_slice::<U1, DimNameDiff<D, U1>>(D::dim() - 1, 0);
         let n = normalizer.tr_dot(&v);
 
         if !n.is_zero() {
@@ -363,9 +363,9 @@ impl<N, D, SA, SB> Transformation<PointBase<N, DimNameDiff<D, U1>, SB>> for Squa
     fn transform_point(&self, pt: &PointBase<N, DimNameDiff<D, U1>, SB>)
         -> PointBase<N, DimNameDiff<D, U1>, SB> {
         let transform   = self.fixed_slice::<DimNameDiff<D, U1>, DimNameDiff<D, U1>>(0, 0);
-        let translation = self.fixed_slice::<DimNameDiff<D, U1>, U1>(0, D::dim());
-        let normalizer  = self.fixed_slice::<U1, DimNameDiff<D, U1>>(D::dim(), 0);
-        let n = normalizer.tr_dot(&pt.coords) + unsafe { *self.get_unchecked(D::dim(), D::dim()) };
+        let translation = self.fixed_slice::<DimNameDiff<D, U1>, U1>(0, D::dim() - 1);
+        let normalizer  = self.fixed_slice::<U1, DimNameDiff<D, U1>>(D::dim() - 1, 0);
+        let n = normalizer.tr_dot(&pt.coords) + unsafe { *self.get_unchecked(D::dim() - 1, D::dim() - 1) };
 
         if !n.is_zero() {
             return transform * (pt / n) + translation;

src/core/constraint.rs

@@ -1,6 +1,8 @@
+//! Compatibility constraints between matrix shapes, e.g., for addition or multiplication.
+
 use core::dimension::{Dim, DimName, Dynamic};
 
-/// A type for enforcing constraints.
+/// A type used in `where` clauses for enforcing constraints.
 pub struct ShapeConstraint;
 
 /// Constraints `C1` and `R2` to be equivalent.
@@ -14,9 +16,12 @@ where ShapeConstraint: DimEq<C1, R2> {
 }
 
 macro_rules! equality_trait_decl(
-    ($($Trait: ident),* $(,)*) => {$(
+    ($($doc: expr, $Trait: ident),* $(,)*) => {$(
         // XXX: we can't do something like `DimEq<D1> for D2` because we would require a blancket impl…
+        #[doc = $doc]
         pub trait $Trait<D1: Dim, D2: Dim> {
+            /// This is either equal to `D1` or `D2`, always choosing the one (if any) which is a type-level
+            /// constant.
             type Representative: Dim;
         }
 
@@ -34,11 +39,26 @@ macro_rules! equality_trait_decl(
     )*}
 );
 
-equality_trait_decl!(DimEq, SameNumberOfRows, SameNumberOfColumns);
+equality_trait_decl!(
+    "Constraints `D1` and `D2` to be equivalent.",
+    DimEq,
 
-/// Constraints D1 and D2 to be equivalent, where the both designates dimensions of algebraic
+    "Constraints `D1` and `D2` to be equivalent. \
+     They are both assumed to be the number of \
+     rows of a matrix.",
+    SameNumberOfRows,
+
+    "Constraints `D1` and `D2` to be equivalent. \
+     They are both assumed to be the number of \
+     columns of a matrix.",
+    SameNumberOfColumns
+);
+
+/// Constraints D1 and D2 to be equivalent, where they both designate dimensions of algebraic
 /// entities (e.g. square matrices).
 pub trait SameDimension<D1: Dim, D2: Dim>: SameNumberOfRows<D1, D2> + SameNumberOfColumns<D1, D2> {
+    /// This is either equal to `D1` or `D2`, always choosing the one (if any) which is a type-level
+    /// constant.
     type Representative: Dim;
 }
 

src/core/construction.rs

@@ -116,6 +116,10 @@ impl<N: Scalar, R: Dim, C: Dim, S: OwnedStorage<N, R, C>> Matrix<N, R, C, S>
         res
     }
 
+    /// Builds a new matrix from its rows.
+    ///
+    /// Panics if not enough rows are provided (for statically-sized matrices), or if all rows do
+    /// not have the same dimensions.
     #[inline]
     pub fn from_rows<SB>(rows: &[Matrix<N, U1, C, SB>]) -> Matrix<N, R, C, S>
         where SB: Storage<N, U1, C> {
@@ -127,13 +131,18 @@ impl<N: Scalar, R: Dim, C: Dim, S: OwnedStorage<N, R, C>> Matrix<N, R, C, S>
 
         if C::try_to_usize().is_none() {
             assert!(rows.iter().all(|r| r.len() == ncols),
-            "The rows provided must all have the same dimension.");
+            "The provided rows must all have the same dimension.");
         }
 
         // FIXME: optimize that.
         Self::from_fn_generic(R::from_usize(nrows), C::from_usize(ncols), |i, j| rows[i][(0, j)])
     }
 
+
+    /// Builds a new matrix from its columns.
+    ///
+    /// Panics if not enough columns are provided (for statically-sized matrices), or if all
+    /// columns do not have the same dimensions.
     #[inline]
     pub fn from_columns<SB>(columns: &[ColumnVector<N, R, SB>]) -> Matrix<N, R, C, S>
         where SB: Storage<N, R, U1> {
@@ -239,12 +248,17 @@ macro_rules! impl_constructors(
                 Self::from_fn_generic($($gargs, )* f)
             }
 
+            /// Creates an identity matrix. If the matrix is not square, the largest square
+            /// submatrix (starting at the first row and column) is set to the identity while all
+            /// other entries are set to zero.
             #[inline]
             pub fn identity($($args: usize,)*) -> Matrix<N $(, $Dims)*, S>
                 where N: Zero + One {
                 Self::identity_generic($($gargs),* )
             }
 
+            /// Creates a matrix filled with its diagonal filled with `elt` and all other
+            /// components set to zero.
             #[inline]
             pub fn from_diagonal_element($($args: usize,)* elt: N) -> Matrix<N $(, $Dims)*, S>
                 where N: Zero + One {

src/core/coordinates.rs

@@ -1,3 +1,9 @@
+#![allow(missing_docs)]
+
+//! Structures to which matrices and vector can be auto-dereferenced (through `Deref`) to access
+//! components using their names. For example, if `v` is a 3D vector, one can write `v.z` instead
+//! of `v[2]`.
+
 use std::mem;
 use std::ops::{Deref, DerefMut};
 
@@ -14,8 +20,10 @@ use core::allocator::OwnedAllocator;
 
 macro_rules! coords_impl(
     ($T: ident; $($comps: ident),*) => {
+        /// Data structure used to provide access to matrix and vector coordinates with the dot
+        /// notation, e.g., `v.x` is the same as `v[0]` for a vector.
         #[repr(C)]
-        #[derive(Eq, PartialEq, RustcEncodable, RustcDecodable, Clone, Hash, Debug, Copy)]
+        #[derive(Eq, PartialEq, Clone, Hash, Debug, Copy, Serialize, Deserialize)]
         pub struct $T<N: Scalar> {
             $(pub $comps: N),*
         }

src/core/decompositions.rs

@@ -272,14 +272,7 @@ impl<N: Real, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
     ///
     /// Matrix symmetricness is not checked. Returns `None` if `self` is not definite positive.
     #[inline]
-    pub fn cholesky(&self) -> Option<OwnedSquareMatrix<N, D, S::Alloc>> {
-        let out = self.transpose();
-        self.do_cholesky(out).ok()
-    }
-
-    /// Cholesky decomposition G of a square symmetric positive definite matrix A, such that A = G * G^T
-    #[inline]
-    pub fn cholesky_unchecked(&self) -> Result<OwnedSquareMatrix<N, D, S::Alloc>, &'static str> {
+    pub fn cholesky(&self) -> Result<OwnedSquareMatrix<N, D, S::Alloc>, &'static str> {
         let out = self.transpose();
 
         if !out.relative_eq(self, N::default_epsilon(), N::default_max_relative()) {
@@ -289,6 +282,13 @@ impl<N: Real, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
         self.do_cholesky(out)
     }
 
+    /// Cholesky decomposition G of a square symmetric positive definite matrix A, such that A = G * G^T
+    #[inline]
+    pub fn cholesky_unchecked(&self) -> Result<OwnedSquareMatrix<N, D, S::Alloc>, &'static str> {
+        let out = self.transpose();
+        self.do_cholesky(out)
+    }
+
     #[inline(always)]
     fn do_cholesky(&self, mut out: OwnedSquareMatrix<N, D, S::Alloc>)
                    -> Result<OwnedSquareMatrix<N, D, S::Alloc>, &'static str> {

src/core/default_allocator.rs

@@ -1,3 +1,8 @@
+//! The default matrix data storage allocator.
+//!
+//! This will use stack-allocated buffers for matrices with dimensions known at compile-time, and
+//! heap-allocated buffers for matrices with at least one dimension unknown at compile-time.
+
 use std::mem;
 use std::ops::Mul;
 

src/core/dimension.rs

@@ -1,10 +1,14 @@
+#![allow(missing_docs)]
+
+//! Traits and tags for identifying the dimension of all algebraic entities.
+
 use std::fmt::Debug;
 use std::any::Any;
 use std::ops::{Add, Sub, Mul, Div};
 use typenum::{self, Unsigned, UInt, B1, Bit, UTerm, Sum, Prod, Diff, Quot};
 
 /// Dim of dynamically-sized algebraic entities.
-#[derive(Clone, Copy, Eq, PartialEq, Debug)]
+#[derive(Clone, Copy, Eq, PartialEq, Debug, Serialize, Deserialize)]
 pub struct Dynamic {
     value: usize
 }
@@ -19,17 +23,27 @@ impl Dynamic {
     }
 }
 
-/// Trait implemented by dynamic dimensions.
+/// Trait implemented by `Dynamic`.
 pub trait IsDynamic { }
-/// Trait implemented by dimensions that are not equal to U1.
+/// Trait implemented by `Dynamic` and type-level integers different from `U1`.
 pub trait IsNotStaticOne { }
 
 impl IsDynamic for Dynamic { }
 impl IsNotStaticOne for Dynamic { }
 
+/// Trait implemented by any type that can be used as a dimension. This includes type-level
+/// integers and `Dynamic` (for dimensions not known at compile-time).
 pub trait Dim: Any + Debug + Copy + PartialEq + Send {
+    /// Gets the compile-time value of `Self`. Returns `None` if it is not known, i.e., if `Self =
+    /// Dynamic`.
     fn try_to_usize() -> Option<usize>;
+
+    /// Gets the run-time value of `self`. For type-level integers, this is the same as
+    /// `Self::try_to_usize().unwrap()`.
     fn value(&self) -> usize;
+
+    /// Builds an instance of `Self` from a run-time value. Panics if `Self` is a type-level
+    /// integer and `dim != Self::try_to_usize().unwrap()`.
     fn from_usize(dim: usize) -> Self;
 }
 
@@ -127,6 +141,7 @@ dim_ops!(
 );
 
 
+/// Trait implemented exclusively by type-level integers.
 pub trait DimName: Dim {
     type Value: NamedDim<Name = Self>;
 
@@ -146,7 +161,7 @@ pub trait NamedDim: Sized + Any + Unsigned {
     type Name: DimName<Value = Self>;
 }
 
-#[derive(Debug, Copy, Clone, Hash, PartialEq, Eq)]
+#[derive(Debug, Copy, Clone, Hash, PartialEq, Eq, Serialize, Deserialize)]
 pub struct U1;
 
 impl Dim for U1 {
@@ -182,7 +197,7 @@ impl NamedDim for typenum::U1{
 
 macro_rules! named_dimension(
     ($($D: ident),* $(,)*) => {$(
-        #[derive(Debug, Copy, Clone, Hash, PartialEq, Eq)]
+        #[derive(Debug, Copy, Clone, Hash, PartialEq, Eq, Serialize, Deserialize)]
         pub struct $D;
 
         impl Dim for $D {

src/core/iter.rs

@@ -1,3 +1,5 @@
+//! Matrix iterators.
+
 use std::marker::PhantomData;
 use std::mem;
 

src/core/matrix.rs

@@ -7,9 +7,9 @@ use std::any::TypeId;
 use std::mem;
 use approx::ApproxEq;
 
-use alga::general::{Field, Real};
+use alga::general::{Ring, Real};
 
-use core::Scalar;
+use core::{Scalar, Unit};
 use core::dimension::{Dim, DimAdd, DimSum, U1, U2};
 use core::constraint::{ShapeConstraint, SameNumberOfRows, SameNumberOfColumns};
 use core::iter::{MatrixIter, MatrixIterMut};
@@ -55,10 +55,37 @@ pub type MatrixTrMul<N, R1, C1, C2, SA> = Matrix<N, C1, C2, TrMulStorage<N, R1,
 pub type MatrixWithScalar<NOld, NNew, R, C, S> =
 Matrix<NNew, R, C, <<S as Storage<NOld, R, C>>::Alloc as Allocator<NNew, R, C>>::Buffer>;
 
+/// The most generic column-major matrix (and vector) type.
+///
+/// It combines four type parameters:
+/// - `N`: for the matrix components scalar type.
+/// - `R`: for the matrix number of rows.
+/// - `C`: for the matrix number of columns.
+/// - `S`: for the matrix data storage, i.e., the buffer that actually contains the matrix
+/// components.
+///
+/// The matrix dimensions parameters `R` and `C` can either be:
+/// - type-level unsigned integer contants (e.g. `U1`, `U124`) from the `nalgebra::` root module.
+/// All numbers from 0 to 127 are defined that way.
+/// - type-level unsigned integer constants (e.g. `U1024`, `U10000`) from the `typenum::` crate.
+/// Using those, you will not get error messages as nice as for numbers smaller than 128 defined on
+/// the `nalgebra::` module.
+/// - the special value `Dynamic` from the `nalgebra::` root module. This indicates that the
+/// specified dimension is not known at compile-time. Note that this will generally imply that the
+/// matrix data storage `S` performs a dynamic allocation and contains extra metadata for the
+/// matrix shape.
+///
+/// Note that mixing `Dynamic` with type-level unsigned integers is allowed. Actually, a
+/// dynamically-sized column vector should be represented as a `Matrix<N, Dynamic, U1, S>` (given
+/// some concrete types for `N` and a compatible data storage type `S`).
 #[repr(C)]
-#[derive(RustcEncodable, RustcDecodable, Hash, Debug, Clone, Copy)]
+#[derive(Serialize, Deserialize, Hash, Debug, Clone, Copy)]
 pub struct Matrix<N: Scalar, R: Dim, C: Dim, S> {
+    /// The data storage that contains all the matrix components and informations about its number
+    /// of rows and column (if needed).
     pub data:  S,
+
+    #[serde(skip_serializing, skip_deserializing)]
     _phantoms: PhantomData<(N, R, C)>
 }
 
@@ -149,7 +176,7 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
         nrows * ncols
     }
 
-    /// The shape (number of rows, number of columns) of this matrix.
+    /// The shape of this matrix returned as the tuple (number of rows, number of columns).
     #[inline]
     pub fn shape(&self) -> (usize, usize) {
         let (nrows, ncols) = self.data.shape();
@@ -405,6 +432,20 @@ impl<N, D, S> ColumnVector<N, D, S>
 
         res
     }
+
+    /// Constructs a vector from coordinates in projective space, i.e., removes a `0` at the end of
+    /// `self`. Returns `None` if this last component is not zero.
+    #[inline]
+    pub fn from_homogeneous<SB>(v: ColumnVector<N, DimSum<D, U1>, SB>) -> Option<OwnedColumnVector<N, D, S::Alloc>>
+        where SB: Storage<N, DimSum<D, U1>, U1, Alloc = S::Alloc> {
+        if v[v.len() - 1].is_zero() {
+            let nrows = D::from_usize(v.len() - 1);
+            Some(v.generic_slice((0, 0), (nrows, U1)).into_owned())
+        }
+        else {
+            None
+        }
+    }
 }
 
 
@@ -677,7 +718,7 @@ impl<N, R: Dim, C: Dim, S> PartialOrd for Matrix<N, R, C, S>
 
         if let Some(mut first_ord) = first_ord {
             let mut it = self.iter().zip(other.iter());
-            it.next(); // Drop the first elements (we already tested it).
+            let _ = it.next(); // Drop the first elements (we already tested it).
 
             for (left, right) in it {
                 if let Some(ord) = left.partial_cmp(right) {
@@ -806,7 +847,7 @@ impl<N, R: Dim, C: Dim, S> fmt::Display for Matrix<N, R, C, S>
 
 
 impl<N, R: Dim, C: Dim, S> Matrix<N, R, C, S>
-    where N: Scalar + Field,
+    where N: Scalar + Ring,
           S: Storage<N, R, C> {
     /// The dot product between two matrices (seen as vectors).
     #[inline]
@@ -1011,3 +1052,35 @@ impl<N, R: Dim, C: Dim, S> Matrix<N, R, C, S>
         }
     }
 }
+
+impl<N, R: Dim, C: Dim, S> ApproxEq for Unit<Matrix<N, R, C, S>>
+    where N: Scalar + ApproxEq,
+          S: Storage<N, R, C>,
+          N::Epsilon: Copy {
+    type Epsilon = N::Epsilon;
+
+    #[inline]
+    fn default_epsilon() -> Self::Epsilon {
+        N::default_epsilon()
+    }
+
+    #[inline]
+    fn default_max_relative() -> Self::Epsilon {
+        N::default_max_relative()
+    }
+
+    #[inline]
+    fn default_max_ulps() -> u32 {
+        N::default_max_ulps()
+    }
+
+    #[inline]
+    fn relative_eq(&self, other: &Self, epsilon: Self::Epsilon, max_relative: Self::Epsilon) -> bool {
+        self.as_ref().relative_eq(other.as_ref(), epsilon, max_relative)
+    }
+
+    #[inline]
+    fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool {
+        self.as_ref().ulps_eq(other.as_ref(), epsilon, max_ulps)
+    }
+}

src/core/matrix_alga.rs

@@ -263,7 +263,7 @@ impl<N, R: DimName, C: DimName, S> FiniteDimInnerSpace for Matrix<N, R, C, S>
                         a = Self::from_column_slice(&[N::zero(), -v[2], v[1]]);
                     };
 
-                    a.normalize_mut();
+                    let _ = a.normalize_mut();
 
                     if f(&a.cross(v)) {
                         f(&a);

src/core/matrix_array.rs

@@ -1,6 +1,11 @@
+use std::mem;
+use std::marker::PhantomData;
 use std::ops::{Deref, DerefMut, Mul};
-use std::fmt::{Debug, Formatter, Result};
+use std::fmt::{self, Debug, Formatter};
 use std::hash::{Hash, Hasher};
+use serde::{Serialize, Serializer, Deserialize, Deserializer};
+use serde::ser::SerializeSeq;
+use serde::de::{SeqVisitor, Visitor};
 
 use typenum::Prod;
 use generic_array::{ArrayLength, GenericArray};
@@ -28,6 +33,7 @@ where R: DimName,
     data: GenericArray<N, Prod<R::Value, C::Value>>
 }
 
+
 impl<N, R, C> Hash for MatrixArray<N, R, C>
 where N: Hash,
       R: DimName,
@@ -70,7 +76,7 @@ where N: Debug,
       R::Value: Mul<C::Value>,
       Prod<R::Value, C::Value>: ArrayLength<N> {
     #[inline]
-    fn fmt(&self, fmt: &mut Formatter) -> Result {
+    fn fmt(&self, fmt: &mut Formatter) -> fmt::Result {
         self.data.fmt(fmt)
     }
 }
@@ -185,3 +191,97 @@ unsafe impl<N, R, C> OwnedStorage<N, R, C> for MatrixArray<N, R, C>
         &mut self[..]
     }
 }
+
+
+/*
+ *
+ * Allocation-less serde impls.
+ *
+ */
+// XXX: open an issue for GenericArray so that it implements serde traits?
+impl<N, R, C> Serialize for MatrixArray<N, R, C>
+where N: Scalar + Serialize,
+      R: DimName,
+      C: DimName,
+      R::Value: Mul<C::Value>,
+      Prod<R::Value, C::Value>: ArrayLength<N> {
+
+
+    fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
+        where S: Serializer {
+            let mut serializer = serializer.serialize_seq_fixed_size(R::dim() * C::dim())?;
+
+            for e in self.iter() {
+                serializer.serialize_element(e)?;
+            }
+
+            serializer.end()
+        }
+}
+
+
+impl<N, R, C> Deserialize for MatrixArray<N, R, C>
+where N: Scalar + Deserialize,
+      R: DimName,
+      C: DimName,
+      R::Value: Mul<C::Value>,
+      Prod<R::Value, C::Value>: ArrayLength<N> {
+
+
+    fn deserialize<D>(deserializer: D) -> Result<Self, D::Error>
+        where D: Deserializer {
+
+            let len = R::dim() * C::dim();
+            deserializer.deserialize_seq_fixed_size(len, MatrixArrayVisitor::new())
+        }
+}
+
+
+/// A visitor that produces a matrix array.
+struct MatrixArrayVisitor<N, R, C> {
+    marker: PhantomData<(N, R, C)>
+}
+
+impl<N, R, C> MatrixArrayVisitor<N, R, C>
+where N: Scalar,
+      R: DimName,
+      C: DimName,
+      R::Value: Mul<C::Value>,
+      Prod<R::Value, C::Value>: ArrayLength<N> {
+
+    /// Construct a new sequence visitor.
+    pub fn new() -> Self {
+        MatrixArrayVisitor {
+            marker: PhantomData,
+        }
+    }
+}
+
+impl<N, R, C> Visitor for MatrixArrayVisitor<N, R, C>
+where N: Scalar + Deserialize,
+      R: DimName,
+      C: DimName,
+      R::Value: Mul<C::Value>,
+      Prod<R::Value, C::Value>: ArrayLength<N> {
+
+    type Value = MatrixArray<N, R, C>;
+
+    fn expecting(&self, formatter: &mut Formatter) -> fmt::Result {
+        formatter.write_str("a matrix array")
+    }
+
+    #[inline]
+    fn visit_seq<V>(self, mut visitor: V) -> Result<MatrixArray<N, R, C>, V::Error>
+        where V: SeqVisitor {
+
+        let mut out: Self::Value = unsafe { mem::uninitialized() };
+        let mut curr = 0;
+
+        while let Some(value) = try!(visitor.visit()) {
+            out[curr] = value;
+            curr += 1;
+        }
+
+        Ok(out)
+    }
+}

src/core/matrix_slice.rs

@@ -7,7 +7,8 @@ use core::storage::{Storage, StorageMut, Owned};
 use core::allocator::Allocator;
 
 macro_rules! slice_storage_impl(
-    ($Storage: ident as $SRef: ty; $T: ident.$get_addr: ident ($Ptr: ty as $Ref: ty)) => {
+    ($doc: expr; $Storage: ident as $SRef: ty; $T: ident.$get_addr: ident ($Ptr: ty as $Ref: ty)) => {
+        #[doc = $doc]
         pub struct $T<'a, N: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim, Alloc> {
             ptr:       $Ptr,
             shape:     (R, C),
@@ -54,8 +55,14 @@ macro_rules! slice_storage_impl(
     }
 );
 
-slice_storage_impl!(Storage as &'a S; SliceStorage.get_address_unchecked(*const N as &'a N));
-slice_storage_impl!(StorageMut as &'a mut S; SliceStorageMut.get_address_unchecked_mut(*mut N as &'a mut N));
+slice_storage_impl!("A matrix data storage for a matrix slice. Only contains an internal reference \
+                     to another matrix data storage.";
+    Storage as &'a S; SliceStorage.get_address_unchecked(*const N as &'a N));
+
+slice_storage_impl!("A mutable matrix data storage for mutable matrix slice. Only contains an \
+                     internal mutable reference to another matrix data storage.";
+    StorageMut as &'a mut S; SliceStorageMut.get_address_unchecked_mut(*mut N as &'a mut N)
+);
 
 
 impl<'a, N: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim, Alloc> Copy
@@ -345,6 +352,11 @@ macro_rules! matrix_slice_impl(
                 }
             }
 
+
+            /// Slices this matrix starting at its component `(start.0, start.1)` and with
+            /// `(shape.0, shape.1)` components. Each row (resp. column) of the sliced matrix is
+            /// separated by `steps.0` (resp. `steps.1`) ignored rows (resp. columns) of the
+            /// original matrix.
             #[inline]
             pub fn $slice_with_steps($me: $Me, start: (usize, usize), shape: (usize, usize), steps: (usize, usize))
                 -> $MatrixSlice<N, Dynamic, Dynamic, Dynamic, Dynamic, S::Alloc> {
@@ -371,6 +383,10 @@ macro_rules! matrix_slice_impl(
                 }
             }
 
+            /// Slices this matrix starting at its component `(start.0, start.1)` and with
+            /// `(R::dim(), CSlice::dim())` components. Each row (resp. column) of the sliced
+            /// matrix is separated by `steps.0` (resp. `steps.1`) ignored rows (resp. columns) of
+            /// the original matrix.
             #[inline]
             pub fn $fixed_slice_with_steps<RSlice, CSlice>($me: $Me, start: (usize, usize), steps: (usize, usize))
                 -> $MatrixSlice<N, RSlice, CSlice, Dynamic, Dynamic, S::Alloc>

src/core/matrix_vec.rs

@@ -5,7 +5,6 @@ use core::dimension::{Dim, DimName, Dynamic, U1};
 use core::storage::{Storage, StorageMut, Owned, OwnedStorage};
 use core::default_allocator::DefaultAllocator;
 
-
 /*
  *
  * Storage.
@@ -13,7 +12,7 @@ use core::default_allocator::DefaultAllocator;
  */
 /// A Vec-based matrix data storage. It may be dynamically-sized.
 #[repr(C)]
-#[derive(Eq, Debug, Clone, PartialEq)]
+#[derive(Eq, Debug, Clone, PartialEq, Serialize, Deserialize)]
 pub struct MatrixVec<N, R: Dim, C: Dim> {
     data:   Vec<N>,
     nrows:  R,
@@ -21,8 +20,10 @@ pub struct MatrixVec<N, R: Dim, C: Dim> {
 }
 
 impl<N, R: Dim, C: Dim> MatrixVec<N, R, C> {
+    /// Creates a new dynamic matrix data storage from the given vector and shape.
     #[inline]
     pub fn new(nrows: R, ncols: C, data: Vec<N>) -> MatrixVec<N, R, C> {
+        assert!(nrows.value() * ncols.value() == data.len(), "Data storage buffer dimension mismatch.");
         MatrixVec {
             data:   data,
             nrows:  nrows,

src/core/mod.rs

@@ -1,9 +1,11 @@
+//! Data structures for vector and matrix computations.
+
 pub mod dimension;
 pub mod constraint;
 pub mod allocator;
 pub mod storage;
 pub mod coordinates;
-pub mod ops;
+mod ops;
 pub mod iter;
 pub mod default_allocator;
 

src/core/ops.rs

@@ -33,7 +33,9 @@ impl<N, R: Dim, C: Dim, S> Index<(usize, usize)> for Matrix<N, R, C, S>
 
     #[inline]
     fn index(&self, ij: (usize, usize)) -> &N {
-        assert!(ij < self.shape(), "Matrix index out of bounds.");
+        let shape = self.shape();
+        assert!(ij.0 < shape.0 && ij.1 < shape.1, "Matrix index out of bounds.");
+
         unsafe { self.get_unchecked(ij.0, ij.1) }
     }
 }
@@ -53,7 +55,9 @@ impl<N, R: Dim, C: Dim, S> IndexMut<(usize, usize)> for Matrix<N, R, C, S>
 
     #[inline]
     fn index_mut(&mut self, ij: (usize, usize)) -> &mut N {
-        assert!(ij < self.shape(), "Matrix index out of bounds.");
+        let shape = self.shape();
+        assert!(ij.0 < shape.0 && ij.1 < shape.1, "Matrix index out of bounds.");
+
         unsafe { self.get_unchecked_mut(ij.0, ij.1) }
     }
 }
@@ -121,7 +125,13 @@ macro_rules! componentwise_binop_impl(
                 assert!(self.shape() == right.shape(), "Matrix addition/subtraction dimensions mismatch.");
                 let mut res = self.into_owned_sum::<R2, C2>();
 
-                for (left, right) in res.iter_mut().zip(right.iter()) {
+                // XXX: optimize our iterator!
+                //
+                // Using our own iterator prevents loop unrolling, wich breaks some optimization
+                // (like SIMD). On the other hand, using the slice iterator is 4x faster.
+
+                // for (left, right) in res.iter_mut().zip(right.iter()) {
+                for (left, right) in res.as_mut_slice().iter_mut().zip(right.iter()) {
                     *left = left.$method(*right)
                 }
 
@@ -143,7 +153,13 @@ macro_rules! componentwise_binop_impl(
                 assert!(self.shape() == right.shape(), "Matrix addition/subtraction dimensions mismatch.");
                 let mut res = right.into_owned_sum::<R1, C1>();
 
-                for (left, right) in self.iter().zip(res.iter_mut()) {
+                // XXX: optimize our iterator!
+                //
+                // Using our own iterator prevents loop unrolling, wich breaks some optimization
+                // (like SIMD). On the other hand, using the slice iterator is 4x faster.
+
+                // for (left, right) in self.iter().zip(res.iter_mut()) {
+                for (left, right) in self.iter().zip(res.as_mut_slice().iter_mut()) {
                     *right = left.$method(*right)
                 }
 
@@ -237,7 +253,13 @@ macro_rules! componentwise_scalarop_impl(
             fn $method(self, rhs: N) -> Self::Output {
                 let mut res = self.into_owned();
 
-                for left in res.iter_mut() {
+                // XXX: optimize our iterator!
+                //
+                // Using our own iterator prevents loop unrolling, wich breaks some optimization
+                // (like SIMD). On the other hand, using the slice iterator is 4x faster.
+
+                // for left in res.iter_mut() {
+                for left in res.as_mut_slice().iter_mut() {
                     *left = left.$method(rhs)
                 }
 
@@ -282,7 +304,13 @@ macro_rules! left_scalar_mul_impl(
             fn mul(self, right: Matrix<$T, R, C, S>) -> Self::Output {
                 let mut res = right.into_owned();
 
-                for right in res.iter_mut() {
+                // XXX: optimize our iterator!
+                //
+                // Using our own iterator prevents loop unrolling, wich breaks some optimization
+                // (like SIMD). On the other hand, using the slice iterator is 4x faster.
+
+                // for right in res.iter_mut() {
+                for right in res.as_mut_slice().iter_mut() {
                     *right = self * *right
                 }
 

src/core/storage.rs

@@ -1,3 +1,5 @@
+//! Abstract definition of a matrix data storage.
+
 use std::mem;
 use std::any::Any;
 
@@ -110,6 +112,12 @@ pub unsafe trait Storage<N: Scalar, R: Dim, C: Dim>: Sized {
     }
 }
 
+
+/// Trait implemented by matrix data storage that can provide a mutable access to its elements.
+///
+/// Note that a mutable access does not mean that the matrix owns its data. For example, a mutable
+/// matrix slice can provide mutable access to its elements even if it does not own its data (it
+/// contains only an internal reference to them).
 pub unsafe trait StorageMut<N: Scalar, R: Dim, C: Dim>: Storage<N, R, C> {
     /// The matrix mutable data pointer.
     fn ptr_mut(&mut self) -> *mut N;

src/core/unit.rs

@@ -1,4 +1,5 @@
-use std::ops::Neg;
+use std::mem;
+use std::ops::{Neg, Deref};
 use approx::ApproxEq;
 
 use alga::general::SubsetOf;
@@ -9,7 +10,7 @@ use alga::linear::NormedSpace;
 ///
 /// Use `.as_ref()` or `.unwrap()` to obtain the undelying value by-reference or by-move.
 #[repr(C)]
-#[derive(Eq, PartialEq, Clone, Hash, Debug, Copy)]
+#[derive(Eq, PartialEq, Clone, Hash, Debug, Copy, Serialize, Deserialize)]
 pub struct Unit<T> {
     value: T
 }
@@ -114,34 +115,34 @@ where T::Field: ApproxEq {
     }
 }
 
-impl<T: ApproxEq> ApproxEq for Unit<T> {
-    type Epsilon = T::Epsilon;
-
-    #[inline]
-    fn default_epsilon() -> Self::Epsilon {
-        T::default_epsilon()
-    }
-
-    #[inline]
-    fn default_max_relative() -> Self::Epsilon {
-        T::default_max_relative()
-    }
-
-    #[inline]
-    fn default_max_ulps() -> u32 {
-        T::default_max_ulps()
-    }
-
-    #[inline]
-    fn relative_eq(&self, other: &Self, epsilon: Self::Epsilon, max_relative: Self::Epsilon) -> bool {
-        self.value.relative_eq(&other.value, epsilon, max_relative)
-    }
-
-    #[inline]
-    fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool {
-        self.value.ulps_eq(&other.value, epsilon, max_ulps)
-    }
-}
+// impl<T: ApproxEq> ApproxEq for Unit<T> {
+//     type Epsilon = T::Epsilon;
+// 
+//     #[inline]
+//     fn default_epsilon() -> Self::Epsilon {
+//         T::default_epsilon()
+//     }
+// 
+//     #[inline]
+//     fn default_max_relative() -> Self::Epsilon {
+//         T::default_max_relative()
+//     }
+// 
+//     #[inline]
+//     fn default_max_ulps() -> u32 {
+//         T::default_max_ulps()
+//     }
+// 
+//     #[inline]
+//     fn relative_eq(&self, other: &Self, epsilon: Self::Epsilon, max_relative: Self::Epsilon) -> bool {
+//         self.value.relative_eq(&other.value, epsilon, max_relative)
+//     }
+// 
+//     #[inline]
+//     fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool {
+//         self.value.ulps_eq(&other.value, epsilon, max_ulps)
+//     }
+// }
 
 
 // FIXME:re-enable this impl when spacialization is possible.
@@ -163,3 +164,12 @@ impl<T: Neg> Neg for Unit<T> {
         Unit::new_unchecked(-self.value)
     }
 }
+
+impl<T> Deref for Unit<T> {
+    type Target = T;
+
+    #[inline]
+    fn deref(&self) -> &T {
+        unsafe { mem::transmute(self) }
+    }
+}

src/geometry/isometry.rs

@@ -11,13 +11,23 @@ use core::storage::{Storage, OwnedStorage};
 use core::allocator::{Allocator, OwnedAllocator};
 use geometry::{TranslationBase, PointBase};
 
+
+/// An isometry that uses a data storage deduced from the allocator `A`.
+pub type OwnedIsometryBase<N, D, A, R> =
+    IsometryBase<N, D, <A as Allocator<N, D, U1>>::Buffer, R>;
+
 /// A direct isometry, i.e., a rotation followed by a translation.
 #[repr(C)]
-#[derive(Hash, Debug, Clone, Copy)]
+#[derive(Hash, Debug, Clone, Copy, Serialize, Deserialize)]
 pub struct IsometryBase<N: Scalar, D: DimName, S, R> {
+    /// The pure rotational part of this isometry.
     pub rotation:    R,
+    /// The pure translational part of this isometry.
     pub translation: TranslationBase<N, D, S>,
-    // One private field just to prevent explicit construction.
+
+
+    // One dummy private field just to prevent explicit construction.
+    #[serde(skip_serializing, skip_deserializing)]
     _noconstruct: PhantomData<N>
 }
 

src/geometry/isometry_ops.rs

@@ -214,8 +214,8 @@ isometry_binop_assign_impl_all!(
 );
 
 
-// IsometryBase × RotationBase
-// IsometryBase ÷ RotationBase
+// IsometryBase × R
+// IsometryBase ÷ R
 isometry_binop_impl_all!(
     Mul, mul;
     self: IsometryBase<N, D, S, R>, rhs: R, Output = IsometryBase<N, D, S, R>;
@@ -371,9 +371,9 @@ isometry_from_composition_impl_all!(
     self: UnitQuaternionBase<N, SA>, right: TranslationBase<N, U3, SB>,
     Output = IsometryBase<N, U3, SB, UnitQuaternionBase<N, SA>>;
     [val val] => IsometryBase::from_parts(TranslationBase::from_vector(&self *  right.vector), self);
-    [ref val] => IsometryBase::from_parts(TranslationBase::from_vector(self * right.vector), self.clone());
+    [ref val] => IsometryBase::from_parts(TranslationBase::from_vector( self *  right.vector), self.clone());
     [val ref] => IsometryBase::from_parts(TranslationBase::from_vector(&self * &right.vector), self);
-    [ref ref] => IsometryBase::from_parts(TranslationBase::from_vector(self * &right.vector), self.clone());
+    [ref ref] => IsometryBase::from_parts(TranslationBase::from_vector( self * &right.vector), self.clone());
 );
 
 // RotationBase × IsometryBase

src/geometry/mod.rs

@@ -1,3 +1,5 @@
+//! Data structures for points and usual transformations (rotations, isometries, etc.)
+
 mod op_macros;
 
 mod point;
@@ -28,6 +30,7 @@ mod unit_complex;
 mod unit_complex_construction;
 mod unit_complex_ops;
 mod unit_complex_alga;
+mod unit_complex_conversion;
 
 mod translation;
 mod translation_construction;

src/geometry/orthographic.rs

@@ -13,7 +13,7 @@ use core::helper;
 use geometry::{PointBase, OwnedPoint};
 
 /// A 3D orthographic projection stored as an homogeneous 4x4 matrix.
-#[derive(Debug, Clone, Copy)] // FIXME: Hash
+#[derive(Debug, Clone, Copy, Serialize, Deserialize)] // FIXME: Hash
 pub struct OrthographicBase<N: Scalar, S: Storage<N, U4, U4>> {
     matrix: SquareMatrix<N, U4, S>
 }
@@ -25,9 +25,7 @@ impl<N, S> Eq for OrthographicBase<N, S>
     where N: Scalar + Eq,
           S: Storage<N, U4, U4> { }
 
-impl<N, S> PartialEq for OrthographicBase<N, S>
-    where N: Scalar,
-          S: Storage<N, U4, U4> {
+impl<N: Scalar, S: Storage<N, U4, U4>> PartialEq for OrthographicBase<N, S> {
     #[inline]
     fn eq(&self, right: &Self) -> bool {
         self.matrix == right.matrix
@@ -80,10 +78,7 @@ impl<N, S> OrthographicBase<N, S>
     }
 }
 
-impl<N, S> OrthographicBase<N, S>
-    where N: Real,
-          S: Storage<N, U4, U4> {
-
+impl<N: Real, S: Storage<N, U4, U4>> OrthographicBase<N, S> {
     /// A reference to the underlying homogeneous transformation matrix.
     #[inline]
     pub fn as_matrix(&self) -> &SquareMatrix<N, U4, S> {
@@ -96,6 +91,26 @@ impl<N, S> OrthographicBase<N, S>
         self.matrix
     }
 
+    /// Retrieves the inverse of the underlying homogeneous matrix.
+    #[inline]
+    pub fn inverse(&self) -> OwnedSquareMatrix<N, U4, S::Alloc> {
+        let mut res = self.to_homogeneous();
+
+        let inv_m11 = N::one() / self.matrix[(0, 0)];
+        let inv_m22 = N::one() / self.matrix[(1, 1)];
+        let inv_m33 = N::one() / self.matrix[(2, 2)];
+
+        res[(0, 0)] = inv_m11;
+        res[(1, 1)] = inv_m22;
+        res[(2, 2)] = inv_m33;
+
+        res[(0, 3)] = -self.matrix[(0, 3)] * inv_m11; 
+        res[(1, 3)] = -self.matrix[(1, 3)] * inv_m22; 
+        res[(2, 3)] = -self.matrix[(2, 3)] * inv_m33; 
+
+        res
+    }
+
     /// Computes the corresponding homogeneous matrix.
     #[inline]
     pub fn to_homogeneous(&self) -> OwnedSquareMatrix<N, U4, S::Alloc> {
@@ -151,6 +166,18 @@ impl<N, S> OrthographicBase<N, S>
         )
     }
 
+    /// Un-projects a point. Faster than multiplication by the underlying matrix inverse.
+    #[inline]
+    pub fn unproject_point<SB>(&self, p: &PointBase<N, U3, SB>) -> OwnedPoint<N, U3, SB::Alloc>
+        where SB: Storage<N, U3, U1> {
+
+        OwnedPoint::<N, U3, SB::Alloc>::new(
+            (p[0] - self.matrix[(0, 3)]) / self.matrix[(0, 0)],
+            (p[1] - self.matrix[(1, 3)]) / self.matrix[(1, 1)],
+            (p[2] - self.matrix[(2, 3)]) / self.matrix[(2, 2)]
+        )
+    }
+
     // FIXME: when we get specialization, specialize the Mul impl instead.
     /// Projects a vector. Faster than matrix multiplication.
     #[inline]
@@ -165,9 +192,7 @@ impl<N, S> OrthographicBase<N, S>
     }
 }
 
-impl<N, S> OrthographicBase<N, S>
-    where N: Real,
-          S: StorageMut<N, U4, U4> {
+impl<N: Real, S: StorageMut<N, U4, U4>> OrthographicBase<N, S> {
     /// Sets the smallest x-coordinate of the view cuboid.
     #[inline]
     pub fn set_left(&mut self, left: N) {

src/geometry/perspective.rs

@@ -13,7 +13,7 @@ use core::helper;
 use geometry::{PointBase, OwnedPoint};
 
 /// A 3D perspective projection stored as an homogeneous 4x4 matrix.
-#[derive(Debug, Clone, Copy)] // FIXME: Hash
+#[derive(Debug, Clone, Copy, Serialize, Deserialize)] // FIXME: Hash
 pub struct PerspectiveBase<N: Scalar, S: Storage<N, U4, U4>> {
     matrix: SquareMatrix<N, U4, S>
 }
@@ -85,6 +85,25 @@ impl<N, S> PerspectiveBase<N, S>
         self.matrix
     }
 
+    /// Retrieves the inverse of the underlying homogeneous matrix.
+    #[inline]
+    pub fn inverse(&self) -> OwnedSquareMatrix<N, U4, S::Alloc> {
+        let mut res = self.to_homogeneous();
+
+        res[(0, 0)] = N::one() / self.matrix[(0, 0)];
+        res[(1, 1)] = N::one() / self.matrix[(1, 1)];
+        res[(2, 2)] = N::zero();
+
+        let m23 = self.matrix[(2, 3)];
+        let m32 = self.matrix[(3, 2)];
+
+        res[(2, 3)] = N::one() / m32;
+        res[(3, 2)] = N::one() / m23;
+        res[(3, 3)] = -self.matrix[(2, 2)] / (m23 * m32);
+
+        res
+    }
+
     /// Computes the corresponding homogeneous matrix.
     #[inline]
     pub fn to_homogeneous(&self) -> OwnedSquareMatrix<N, U4, S::Alloc> {
@@ -137,6 +156,20 @@ impl<N, S> PerspectiveBase<N, S>
         )
     }
 
+    /// Un-projects a point. Faster than multiplication by the matrix inverse.
+    #[inline]
+    pub fn unproject_point<SB>(&self, p: &PointBase<N, U3, SB>) -> OwnedPoint<N, U3, SB::Alloc>
+        where SB: Storage<N, U3, U1> {
+
+        let inverse_denom = self.matrix[(2, 3)] / (p[2] + self.matrix[(2, 2)]);
+
+        OwnedPoint::<N, U3, SB::Alloc>::new(
+            p[0] * inverse_denom / self.matrix[(0, 0)],
+            p[1] * inverse_denom / self.matrix[(1, 1)],
+            -inverse_denom
+        )
+    }
+
     // FIXME: when we get specialization, specialize the Mul impl instead.
     /// Projects a vector. Faster than matrix multiplication.
     #[inline]

src/geometry/point.rs

@@ -9,13 +9,11 @@ use core::dimension::{DimName, DimNameSum, DimNameAdd, U1};
 use core::storage::{Storage, StorageMut, MulStorage};
 use core::allocator::{Allocator, SameShapeR};
 
-pub type PointSumStorage<N, R1, R2, SA> =
-    <<SA as Storage<N, R1, U1>>::Alloc as Allocator<N, SameShapeR<R1, R2>, U1>>::Buffer;
-
 // XXX Bad name: we can't even add points…
 /// The type of the result of the sum of a point with a vector.
 pub type PointSum<N, D1, D2, SA> =
-    PointBase<N, SameShapeR<D1, D2>, PointSumStorage<N, D1, D2, SA>>;
+    PointBase<N, SameShapeR<D1, D2>,
+    <<SA as Storage<N, D1, U1>>::Alloc as Allocator<N, SameShapeR<D1, D2>, U1>>::Buffer>;
 
 /// The type of the result of the multiplication of a point by a matrix.
 pub type PointMul<N, R1, C1, SA> = PointBase<N, R1, MulStorage<N, R1, C1, U1, SA>>;
@@ -23,10 +21,11 @@ pub type PointMul<N, R1, C1, SA> = PointBase<N, R1, MulStorage<N, R1, C1, U1, SA
 /// A point with an owned storage.
 pub type OwnedPoint<N, D, A> = PointBase<N, D, <A as Allocator<N, D, U1>>::Buffer>;
 
-/// A point in n-dimensional euclidean space.
+/// A point in a n-dimensional euclidean space.
 #[repr(C)]
-#[derive(Hash, Debug)]
+#[derive(Hash, Debug, Serialize, Deserialize)]
 pub struct PointBase<N: Scalar, D: DimName, S: Storage<N, D, U1>> {
+    /// The coordinates of this point, i.e., the shift from the origin.
     pub coords: ColumnVector<N, D, S>
 }
 

src/geometry/quaternion.rs

@@ -21,8 +21,9 @@ pub type OwnedUnitQuaternionBase<N, A> = UnitQuaternionBase<N, <A as Allocator<N
 /// A quaternion. See the type alias `UnitQuaternionBase = Unit<QuaternionBase>` for a quaternion
 /// that may be used as a rotation.
 #[repr(C)]
-#[derive(Hash, Debug, Copy, Clone)]
+#[derive(Hash, Debug, Copy, Clone, Serialize, Deserialize)]
 pub struct QuaternionBase<N: Real, S: Storage<N, U4, U1>> {
+    /// This quaternion as a 4D vector of coordinates in the `[ x, y, z, w ]` storage order.
     pub coords: ColumnVector<N, U4, S>
 }
 
@@ -86,6 +87,12 @@ impl<N, S> QuaternionBase<N, S>
         self.coords.norm_squared()
     }
 
+    /// Normalizes this quaternion.
+    #[inline]
+    pub fn normalize(&self) -> OwnedQuaternionBase<N, S::Alloc> {
+        QuaternionBase::from_vector(self.coords.normalize())
+    }
+
     /// Compute the conjugate of this quaternion.
     #[inline]
     pub fn conjugate(&self) -> OwnedQuaternionBase<N, S::Alloc> {
@@ -108,6 +115,13 @@ impl<N, S> QuaternionBase<N, S>
             None
         }
     }
+
+    /// Linear interpolation between two quaternion.
+    #[inline]
+    pub fn lerp<S2>(&self, other: &QuaternionBase<N, S2>, t: N) -> OwnedQuaternionBase<N, S::Alloc>
+        where S2: Storage<N, U4, U1> {
+        self * (N::one() - t) + other * t
+    }
 }
 
 
@@ -208,6 +222,12 @@ impl<N, S> QuaternionBase<N, S>
             true
         }
     }
+
+    /// Normalizes this quaternion.
+    #[inline]
+    pub fn normalize_mut(&mut self) -> N {
+        self.coords.normalize_mut()
+    }
 }
 
 impl<N, S> ApproxEq for QuaternionBase<N, S>
@@ -262,6 +282,18 @@ pub type UnitQuaternionBase<N, S> = Unit<QuaternionBase<N, S>>;
 impl<N, S> UnitQuaternionBase<N, S>
     where N: Real,
           S: Storage<N, U4, U1> {
+    /// Moves this unit quaternion into one that owns its data.
+    #[inline]
+    pub fn into_owned(self) -> OwnedUnitQuaternionBase<N, S::Alloc> {
+        UnitQuaternionBase::new_unchecked(self.unwrap().into_owned())
+    }
+
+    /// Clones this unit quaternion into one that owns its data.
+    #[inline]
+    pub fn clone_owned(&self) -> OwnedUnitQuaternionBase<N, S::Alloc> {
+        UnitQuaternionBase::new_unchecked(self.as_ref().clone_owned())
+    }
+
     /// The rotation angle in [0; pi] of this unit quaternion.
     #[inline]
     pub fn angle(&self) -> N {
@@ -312,6 +344,75 @@ impl<N, S> UnitQuaternionBase<N, S>
         where S2: Storage<N, U4, U1> {
         other / self
     }
+
+    /// Linear interpolation between two unit quaternions.
+    ///
+    /// The result is not normalized.
+    #[inline]
+    pub fn lerp<S2>(&self, other: &UnitQuaternionBase<N, S2>, t: N) -> OwnedQuaternionBase<N, S::Alloc>
+        where S2: Storage<N, U4, U1> {
+        self.as_ref().lerp(other.as_ref(), t)
+    }
+
+    /// Normalized linear interpolation between two unit quaternions.
+    #[inline]
+    pub fn nlerp<S2>(&self, other: &UnitQuaternionBase<N, S2>, t: N) -> OwnedUnitQuaternionBase<N, S::Alloc>
+        where S2: Storage<N, U4, U1> {
+        let mut res = self.lerp(other, t);
+        let _ = res.normalize_mut();
+
+        UnitQuaternionBase::new_unchecked(res)
+    }
+
+    /// Spherical linear interpolation between two unit quaternions.
+    ///
+    /// Panics if the angle between both quaternion is 180 degrees (in which case the interpolation
+    /// is not well-defined).
+    #[inline]
+    pub fn slerp<S2>(&self, other: &UnitQuaternionBase<N, S2>, t: N) -> OwnedUnitQuaternionBase<N, S::Alloc>
+        where S2: Storage<N, U4, U1, Alloc = S::Alloc> {
+        self.try_slerp(other, t, N::zero()).expect(
+            "Unable to perform a spherical quaternion interpolation when they \
+             are 180 degree apart (the result is not unique).")
+    }
+
+    /// Computes the spherical linear interpolation between two unit quaternions or returns `None`
+    /// if both quaternions are approximately 180 degrees apart (in which case the interpolation is
+    /// not well-defined).
+    ///
+    /// # Arguments
+    /// * `self`: the first quaternion to interpolate from.
+    /// * `other`: the second quaternion to interpolate toward.
+    /// * `t`: the interpolation parameter. Should be between 0 and 1.
+    /// * `epsilon`: the value bellow which the sinus of the angle separating both quaternion
+    /// must be to return `None`.
+    #[inline]
+    pub fn try_slerp<S2>(&self, other: &UnitQuaternionBase<N, S2>, t: N, epsilon: N)
+                         -> Option<OwnedUnitQuaternionBase<N, S::Alloc>>
+        where S2: Storage<N, U4, U1, Alloc = S::Alloc> {
+
+        let c_hang = self.coords.dot(&other.coords);
+
+        // self == other
+        if c_hang.abs() >= N::one() {
+            return Some(self.clone_owned())
+        }
+
+        let hang   = c_hang.acos();
+        let s_hang = (N::one() - c_hang * c_hang).sqrt();
+
+        // FIXME: what if s_hang is 0.0 ? The result is not well-defined.
+        if relative_eq!(s_hang, N::zero(), epsilon = epsilon) {
+            None
+        }
+        else {
+            let ta = ((N::one() - t) * hang).sin() / s_hang;
+            let tb = (t * hang).sin() / s_hang; 
+            let res = self.as_ref() * ta + other.as_ref() * tb;
+
+            Some(UnitQuaternionBase::new_unchecked(res))
+        }
+    }
 }
 
 impl<N, S> UnitQuaternionBase<N, S>
@@ -453,3 +554,34 @@ impl<N, S> fmt::Display for UnitQuaternionBase<N, S>
         }
     }
 }
+
+impl<N, S> ApproxEq for UnitQuaternionBase<N, S>
+    where N: Real + ApproxEq<Epsilon = N>,
+          S: Storage<N, U4, U1> {
+    type Epsilon = N;
+
+    #[inline]
+    fn default_epsilon() -> Self::Epsilon {
+        N::default_epsilon()
+    }
+
+    #[inline]
+    fn default_max_relative() -> Self::Epsilon {
+        N::default_max_relative()
+    }
+
+    #[inline]
+    fn default_max_ulps() -> u32 {
+        N::default_max_ulps()
+    }
+
+    #[inline]
+    fn relative_eq(&self, other: &Self, epsilon: Self::Epsilon, max_relative: Self::Epsilon) -> bool {
+        self.as_ref().relative_eq(other.as_ref(), epsilon, max_relative)
+    }
+
+    #[inline]
+    fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool {
+        self.as_ref().ulps_eq(other.as_ref(), epsilon, max_ulps)
+    }
+}

src/geometry/quaternion_construction.rs

@@ -33,7 +33,7 @@ impl<N, S> QuaternionBase<N, S>
     /// Creates a new quaternion from its individual components. Note that the arguments order does
     /// **not** follow the storage order.
     ///
-    /// The storage order is [ x, y, z, w ].
+    /// The storage order is `[ x, y, z, w ]`.
     #[inline]
     pub fn new(w: N, x: N, y: N, z: N) -> Self {
         let v = ColumnVector::<N, U4, S>::new(x, y, z, w);
@@ -169,19 +169,44 @@ impl<N, S> UnitQuaternionBase<N, S>
     pub fn from_rotation_matrix<SB>(rotmat: &RotationBase<N, U3, SB>) -> Self
         where SB: Storage<N, U3, U3>,
               SB::Alloc: Allocator<N, U3, U1> {
-        let angle = rotmat.angle();
 
-        if let Some(axis) = rotmat.axis() {
-            Self::from_axis_angle(&axis, angle)
+        // Robust matrix to quaternion transformation.
+        // See http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion
+        let tr = rotmat[(0, 0)] + rotmat[(1, 1)] + rotmat[(2, 2)];
+        let res;
+
+        let _0_25: N = ::convert(0.25);
+
+        if tr > N::zero() {
+          let denom = (tr + N::one()).sqrt() * ::convert(2.0);
+          res = QuaternionBase::new(_0_25 * denom,
+                                    (rotmat[(2, 1)] - rotmat[(1, 2)]) / denom,
+                                    (rotmat[(0, 2)] - rotmat[(2, 0)]) / denom,
+                                    (rotmat[(1, 0)] - rotmat[(0, 1)]) / denom);
         }
-        else if angle > ::convert(1.0f64) {
-            // The angle is 3.14.
-            -Self::identity()
+        else if rotmat[(0, 0)] > rotmat[(1, 1)] && rotmat[(0, 0)] > rotmat[(2, 2)] {
+          let denom = (N::one() + rotmat[(0, 0)] - rotmat[(1, 1)] - rotmat[(2, 2)]).sqrt() * ::convert(2.0);
+          res = QuaternionBase::new((rotmat[(2, 1)] - rotmat[(1, 2)]) / denom,
+                                    _0_25 * denom,
+                                    (rotmat[(0, 1)] + rotmat[(1, 0)]) / denom,
+                                    (rotmat[(0, 2)] + rotmat[(2, 0)]) / denom);
+        }
+        else if rotmat[(1, 1)] > rotmat[(2, 2)] {
+          let denom = (N::one() + rotmat[(1, 1)] - rotmat[(0, 0)] - rotmat[(2, 2)]).sqrt() * ::convert(2.0);
+          res = QuaternionBase::new((rotmat[(0, 2)] - rotmat[(2, 0)]) / denom,
+                                    (rotmat[(0, 1)] + rotmat[(1, 0)]) / denom,
+                                    _0_25 * denom,
+                                    (rotmat[(1, 2)] + rotmat[(2, 1)]) / denom);
         }
         else {
-            // The angle is 0.
-            Self::identity()
+          let denom = (N::one() + rotmat[(2, 2)] - rotmat[(0, 0)] - rotmat[(1, 1)]).sqrt() * ::convert(2.0);
+          res = QuaternionBase::new((rotmat[(1, 0)] - rotmat[(0, 1)]) / denom,
+                                    (rotmat[(0, 2)] + rotmat[(2, 0)]) / denom,
+                                    (rotmat[(1, 2)] + rotmat[(2, 1)]) / denom,
+                                    _0_25 * denom);
         }
+
+        Self::new_unchecked(res)
     }
 
     /// The unit quaternion needed to make `a` and `b` be collinear and point toward the same

src/geometry/quaternion_coordinates.rs

@@ -1,5 +1,5 @@
-use std::ops::{Deref, DerefMut};
 use std::mem;
+use std::ops::{Deref, DerefMut};
 
 use alga::general::Real;
 

src/geometry/rotation.rs

@@ -15,8 +15,8 @@ pub type OwnedRotation<N, D, A> = RotationBase<N, D, <A as Allocator<N, D, D>>::
 
 /// A rotation matrix.
 #[repr(C)]
-#[derive(Hash, Debug, Clone, Copy)]
-pub struct RotationBase<N: Scalar, D: DimName, S: Storage<N, D, D>> {
+#[derive(Hash, Debug, Clone, Copy, Serialize, Deserialize)]
+pub struct RotationBase<N: Scalar, D: DimName, S> {
     matrix: SquareMatrix<N, D, S>
 }
 

src/geometry/rotation_alga.rs

@@ -201,20 +201,23 @@ impl<N, D: DimName, SA, SB> Rotation<PointBase<N, D, SB>> for RotationBase<N, D,
           SA::Alloc: OwnedAllocator<N, D, D, SA>,
           SB::Alloc: OwnedAllocator<N, D, U1, SB> {
     #[inline]
-    fn powf(&self, n: N) -> Option<Self> {
-        // XXX: add the general case.
+    fn powf(&self, _: N) -> Option<Self> {
+        // XXX: Add the general case.
+        // XXX: Use specialization for 2D and 3D.
         unimplemented!()
     }
 
     #[inline]
-    fn rotation_between(a: &ColumnVector<N, D, SB>, b: &ColumnVector<N, D, SB>) -> Option<Self> {
-        // XXX: add the general case.
+    fn rotation_between(_: &ColumnVector<N, D, SB>, _: &ColumnVector<N, D, SB>) -> Option<Self> {
+        // XXX: Add the general case.
+        // XXX: Use specialization for 2D and 3D.
         unimplemented!()
     }
 
     #[inline]
-    fn scaled_rotation_between(a: &ColumnVector<N, D, SB>, b: &ColumnVector<N, D, SB>, n: N) -> Option<Self> {
-        // XXX: add the general case.
+    fn scaled_rotation_between(_: &ColumnVector<N, D, SB>, _: &ColumnVector<N, D, SB>, _: N) -> Option<Self> {
+        // XXX: Add the general case.
+        // XXX: Use specialization for 2D and 3D.
         unimplemented!()
     }
 }

src/geometry/rotation_ops.rs

@@ -118,7 +118,7 @@ md_impl_all!(
 );
 
 
-// RotationBase *= RotationBase
+// RotationBase ×= RotationBase
 // FIXME: try not to call `inverse()` explicitly.
 
 md_assign_impl_all!(

src/geometry/rotation_specialization.rs

@@ -44,7 +44,7 @@ where N: Real,
     pub fn rotation_between<SB, SC>(a: &ColumnVector<N, U2, SB>, b: &ColumnVector<N, U2, SC>) -> Self
         where SB: Storage<N, U2, U1>,
               SC: Storage<N, U2, U1> {
-        UnitComplex::rotation_between(a, b).to_rotation_matrix()
+        ::convert(UnitComplex::rotation_between(a, b).to_rotation_matrix())
     }
 
     /// The smallest rotation needed to make `a` and `b` collinear and point toward the same
@@ -53,7 +53,7 @@ where N: Real,
     pub fn scaled_rotation_between<SB, SC>(a: &ColumnVector<N, U2, SB>, b: &ColumnVector<N, U2, SC>, s: N) -> Self
         where SB: Storage<N, U2, U1>,
               SC: Storage<N, U2, U1> {
-        UnitComplex::scaled_rotation_between(a, b, s).to_rotation_matrix()
+        ::convert(UnitComplex::scaled_rotation_between(a, b, s).to_rotation_matrix())
     }
 }
 
@@ -81,6 +81,8 @@ where N: Real,
         other * self.inverse()
     }
 
+    /// Raise the quaternion to a given floating power, i.e., returns the rotation with the angle
+    /// of `self` multiplied by `n`.
     #[inline]
     pub fn powf(&self, n: N) -> OwnedRotation<N, U2, S::Alloc> {
         OwnedRotation::<_, _, S::Alloc>::new(self.angle() * n)
@@ -338,6 +340,8 @@ where N: Real,
         other * self.inverse()
     }
 
+    /// Raise the quaternion to a given floating power, i.e., returns the rotation with the same
+    /// axis as `self` and an angle equal to `self.angle()` multiplied by `n`.
     #[inline]
     pub fn powf(&self, n: N) -> OwnedRotation<N, U3, S::Alloc> {
         if let Some(axis) = self.axis() {

src/geometry/similarity.rs

@@ -10,10 +10,16 @@ use core::storage::{Storage, OwnedStorage};
 use core::allocator::{Allocator, OwnedAllocator};
 use geometry::{PointBase, TranslationBase, IsometryBase};
 
+
+/// A similarity that uses a data storage deduced from the allocator `A`.
+pub type OwnedSimilarityBase<N, D, A, R> =
+    SimilarityBase<N, D, <A as Allocator<N, D, U1>>::Buffer, R>;
+
 /// A similarity, i.e., an uniform scaling, followed by a rotation, followed by a translation.
 #[repr(C)]
-#[derive(Hash, Debug, Clone, Copy)]
+#[derive(Hash, Debug, Clone, Copy, Serialize, Deserialize)]
 pub struct SimilarityBase<N: Scalar, D: DimName, S, R> {
+    /// The part of this similarity that does not include the scaling factor.
     pub isometry: IsometryBase<N, D, S, R>,
     scaling:      N
 }
@@ -62,12 +68,6 @@ impl<N, D: DimName, S, R> SimilarityBase<N, D, S, R>
         self.isometry.translation.vector *= self.scaling;
     }
 
-    /// The scaling factor of this similarity transformation.
-    #[inline]
-    pub fn scaling(&self) -> N {
-        self.scaling
-    }
-
     /// The scaling factor of this similarity transformation.
     #[inline]
     pub fn set_scaling(&mut self, scaling: N) {
@@ -162,6 +162,12 @@ impl<N, D: DimName, S, R> SimilarityBase<N, D, S, R>
 
         res
     }
+
+    /// The scaling factor of this similarity transformation.
+    #[inline]
+    pub fn scaling(&self) -> N {
+        self.scaling
+    }
 }
 
 

src/geometry/similarity_alias.rs

@@ -12,7 +12,7 @@ pub type Similarity2<N> = SimilarityBase<N, U2, MatrixArray<N, U2, U1>, UnitComp
 /// A 3-dimensional similarity.
 pub type Similarity3<N> = SimilarityBase<N, U3, MatrixArray<N, U3, U1>, UnitQuaternion<N>>;
 
-/// A 3-dimensional similarity using a rotation matrix for its rotation part.
+/// A 2-dimensional similarity using a rotation matrix for its rotation part.
 pub type SimilarityMatrix2<N> = Similarity<N, U2>;
 
 /// A 3-dimensional similarity using a rotation matrix for its rotation part.

src/geometry/similarity_construction.rs

@@ -12,7 +12,8 @@ use core::dimension::{DimName, U1, U2, U3, U4};
 use core::allocator::{OwnedAllocator, Allocator};
 use core::storage::OwnedStorage;
 
-use geometry::{PointBase, TranslationBase, RotationBase, SimilarityBase, UnitQuaternionBase, IsometryBase};
+use geometry::{PointBase, TranslationBase, RotationBase, SimilarityBase,
+               UnitComplex, UnitQuaternionBase, IsometryBase};
 
 
 impl<N, D: DimName, S, R> SimilarityBase<N, D, S, R>
@@ -106,6 +107,17 @@ impl<N, S, SR> SimilarityBase<N, U2, S, RotationBase<N, U2, SR>>
     }
 }
 
+impl<N, S> SimilarityBase<N, U2, S, UnitComplex<N>>
+    where N: Real,
+          S: OwnedStorage<N, U2, U1>,
+          S::Alloc: OwnedAllocator<N, U2, U1, S> {
+    /// Creates a new similarity from a translation and a rotation angle.
+    #[inline]
+    pub fn new(translation: ColumnVector<N, U2, S>, angle: N, scaling: N) -> Self {
+        Self::from_parts(TranslationBase::from_vector(translation), UnitComplex::new(angle), scaling)
+    }
+}
+
 // 3D rotation.
 macro_rules! similarity_construction_impl(
     ($Rot: ty, $RotId: ident, $RRDim: ty, $RCDim: ty) => {

src/geometry/transform.rs

@@ -14,6 +14,8 @@ use core::allocator::Allocator;
 ///
 /// NOTE: this trait is not intended to be implementable outside of the `nalgebra` crate.
 pub trait TCategory: Any + Debug + Copy + PartialEq + Send {
+    /// Indicates whether a `Transform` with the category `Self` has a bottom-row different from
+    /// `0 0 .. 1`.
     #[inline]
     fn has_normalizer() -> bool {
         true
@@ -28,16 +30,19 @@ pub trait TCategory: Any + Debug + Copy + PartialEq + Send {
               N::Epsilon: Copy;
 }
 
-/// Traits that gives the transformation category that is compatible with the result of the
+/// Traits that gives the `Transform` category that is compatible with the result of the
 /// multiplication of transformations with categories `Self` and `Other`.
 pub trait TCategoryMul<Other: TCategory>: TCategory {
+    /// The transform category that results from the multiplication of a `Transform<Self>` to a
+    /// `Transform<Other>`. This is usually equal to `Self` or `Other`, whichever is the most
+    /// general category.
     type Representative: TCategory;
 }
 
-/// Indicates that `Self` is a more general transformation category than `Other`.
+/// Indicates that `Self` is a more general `Transform` category than `Other`.
 pub trait SuperTCategoryOf<Other: TCategory>: TCategory { }
 
-/// Indicates that `Self` is a more specific transformation category than `Other`.
+/// Indicates that `Self` is a more specific `Transform` category than `Other`.
 ///
 /// Automatically implemented based on `SuperTCategoryOf`.
 pub trait SubTCategoryOf<Other: TCategory>: TCategory { }
@@ -46,17 +51,17 @@ where T1: TCategory,
       T2: SuperTCategoryOf<T1> {
 }
 
-/// Tag representing the most general (not necessarily inversible) transformation type.
-#[derive(Debug, Copy, Clone, Hash, PartialEq, Eq)]
-pub struct TGeneral;
+/// Tag representing the most general (not necessarily inversible) `Transform` type.
+#[derive(Debug, Copy, Clone, Hash, PartialEq, Eq, Serialize, Deserialize)]
+pub enum TGeneral { }
 
-/// Tag representing the most general inversible transformation type.
-#[derive(Debug, Copy, Clone, Hash, PartialEq, Eq)]
-pub struct TProjective;
+/// Tag representing the most general inversible `Transform` type.
+#[derive(Debug, Copy, Clone, Hash, PartialEq, Eq, Serialize, Deserialize)]
+pub enum TProjective { }
 
-/// Tag representing an affine transformation. Its bottom-row is equal to `(0, 0 ... 0, 1)`.
-#[derive(Debug, Copy, Clone, Hash, PartialEq, Eq)]
-pub struct TAffine;
+/// Tag representing an affine `Transform`. Its bottom-row is equal to `(0, 0 ... 0, 1)`.
+#[derive(Debug, Copy, Clone, Hash, PartialEq, Eq, Serialize, Deserialize)]
+pub enum TAffine { }
 
 impl TCategory for TGeneral {
     #[inline]
@@ -150,9 +155,11 @@ pub type OwnedTransform<N, D, A, C>
 /// It is stored as a matrix with dimensions `(D + 1, D + 1)`, e.g., it stores a 4x4 matrix for a
 /// 3D transformation.
 #[repr(C)]
-#[derive(Debug, Clone, Copy)] // FIXME: Hash
+#[derive(Debug, Clone, Copy, Serialize, Deserialize)] // FIXME: Hash
 pub struct TransformBase<N: Scalar, D: DimNameAdd<U1>, S, C: TCategory> {
     matrix:   SquareMatrix<N, DimNameSum<D, U1>, S>,
+
+    #[serde(skip_serializing, skip_deserializing)]
     _phantom: PhantomData<C>
 }
 

src/geometry/transform_alias.rs

@@ -3,14 +3,27 @@ use core::dimension::{U1, U2, U3, DimNameSum};
 
 use geometry::{TransformBase, TGeneral, TProjective, TAffine};
 
+/// A `D`-dimensional general transformation that may not be inversible. Stored as an homogeneous
+/// `(D + 1) × (D + 1)` matrix.
 pub type Transform<N, D> = TransformBase<N, D, MatrixArray<N, DimNameSum<D, U1>, DimNameSum<D, U1>>, TGeneral>;
+
+/// An inversible `D`-dimensional general transformation. Stored as an homogeneous
+/// `(D + 1) × (D + 1)` matrix.
 pub type Projective<N, D> = TransformBase<N, D, MatrixArray<N, DimNameSum<D, U1>, DimNameSum<D, U1>>, TProjective>;
+
+/// A `D`-dimensional affine transformation. Stored as an homogeneous `(D + 1) × (D + 1)` matrix.
 pub type Affine<N, D> = TransformBase<N, D, MatrixArray<N, DimNameSum<D, U1>, DimNameSum<D, U1>>, TAffine>;
 
+/// A 2D general transformation that may not be inversible. Stored as an homogeneous 3x3 matrix.
 pub type Transform2<N> = Transform<N, U2>;
+/// An inversible 2D general transformation. Stored as an homogeneous 3x3 matrix.
 pub type Projective2<N> = Projective<N, U2>;
+/// A 2D affine transformation. Stored as an homogeneous 3x3 matrix.
 pub type Affine2<N> = Affine<N, U2>;
 
+/// A 3D general transformation that may not be inversible. Stored as an homogeneous 4x4 matrix.
 pub type Transform3<N> = Transform<N, U3>;
+/// An inversible 3D general transformation. Stored as an homogeneous 4x4 matrix.
 pub type Projective3<N> = Projective<N, U3>;
+/// A 3D affine transformation. Stored as an homogeneous 4x4 matrix.
 pub type Affine3<N> = Affine<N, U3>;

src/geometry/transform_ops.rs

@@ -13,6 +13,68 @@ use geometry::{PointBase, OwnedPoint, TransformBase, OwnedTransform, TCategory,
                SubTCategoryOf, SuperTCategoryOf, TGeneral, TProjective, TAffine, RotationBase,
                UnitQuaternionBase, IsometryBase, SimilarityBase, TranslationBase};
 
+/*
+ *
+ * In the following, we provide:
+ * =========================
+ *
+ * Index<(usize, usize)>
+ * IndexMut<(usize, usize)> (where TCategory == TGeneral)
+ *
+ * (Operators)
+ *
+ * TransformBase × IsometryBase
+ * TransformBase × RotationBase
+ * TransformBase × SimilarityBase
+ * TransformBase × TransformBase
+ * TransformBase × UnitQuaternion
+ * FIXME: TransformBase × UnitComplex
+ * TransformBase × TranslationBase
+ * TransformBase × ColumnVector
+ * TransformBase × PointBase
+ *
+ * IsometryBase       × TransformBase
+ * RotationBase       × TransformBase
+ * SimilarityBase     × TransformBase
+ * TranslationBase    × TransformBase
+ * UnitQuaternionBase × TransformBase
+ * FIXME: UnitComplex × TransformBase
+ *
+ * FIXME: TransformBase ÷ IsometryBase
+ * TransformBase ÷ RotationBase
+ * FIXME: TransformBase ÷ SimilarityBase
+ * TransformBase ÷ TransformBase
+ * TransformBase ÷ UnitQuaternion
+ * TransformBase ÷ TranslationBase
+ *
+ * FIXME: IsometryBase       ÷ TransformBase
+ * RotationBase       ÷ TransformBase
+ * FIXME: SimilarityBase     ÷ TransformBase
+ * TranslationBase    ÷ TransformBase
+ * UnitQuaternionBase ÷ TransformBase
+ * FIXME: UnitComplex ÷ TransformBase
+ *
+ *
+ * (Assignment Operators)
+ *
+ *
+ * TransformBase ×= TransformBase
+ * TransformBase ×= SimilarityBase
+ * TransformBase ×= IsometryBase
+ * TransformBase ×= RotationBase
+ * TransformBase ×= UnitQuaternionBase
+ * FIXME: TransformBase ×= UnitComplex
+ * TransformBase ×= TranslationBase
+ *
+ * TransformBase ÷= TransformBase
+ * FIXME: TransformBase ÷= SimilarityBase
+ * FIXME: TransformBase ÷= IsometryBase
+ * TransformBase ÷= RotationBase
+ * TransformBase ÷= UnitQuaternionBase
+ * FIXME: TransformBase ÷= UnitComplex
+ *
+ */
+
 /*
  *
  * Indexing.
@@ -41,61 +103,6 @@ impl<N, D, S> IndexMut<(usize, usize)> for TransformBase<N, D, S, TGeneral>
     }
 }
 
-/*
- *
- * In the following, we provide:
- * =========================
- *
- *
- * (Operators)
- *
- * [OK] TransformBase × IsometryBase
- * [OK] TransformBase × RotationBase
- * [OK] TransformBase × SimilarityBase
- * [OK] TransformBase × TransformBase
- * [OK] TransformBase × UnitQuaternion
- * [OK] TransformBase × TranslationBase
- * [OK] TransformBase × ColumnVector
- * [OK] TransformBase × PointBase
- *
- * [OK] IsometryBase       × TransformBase
- * [OK] RotationBase       × TransformBase
- * [OK] SimilarityBase     × TransformBase
- * [OK] TranslationBase    × TransformBase
- * [OK] UnitQuaternionBase × TransformBase
- *
- * TransformBase ÷ IsometryBase
- * [OK] TransformBase ÷ RotationBase
- * TransformBase ÷ SimilarityBase
- * [OK] TransformBase ÷ TransformBase
- * [OK] TransformBase ÷ UnitQuaternion
- * [OK] TransformBase ÷ TranslationBase
- *
- * IsometryBase       ÷ TransformBase
- * [OK] RotationBase       ÷ TransformBase
- * SimilarityBase     ÷ TransformBase
- * [OK] TranslationBase    ÷ TransformBase
- * [OK] UnitQuaternionBase ÷ TransformBase
- *
- *
- * (Assignment Operators)
- *
- *
- * [OK] TransformBase ×= TransformBase
- * [OK] TransformBase ×= SimilarityBase
- * [OK] TransformBase ×= IsometryBase
- * [OK] TransformBase ×= RotationBase
- * [OK] TransformBase ×= UnitQuaternionBase
- * [OK] TransformBase ×= TranslationBase
- *
- * [OK] TransformBase ÷= TransformBase
- * TransformBase ÷= SimilarityBase
- * TransformBase ÷= IsometryBase
- * [OK] TransformBase ÷= RotationBase
- * [OK] TransformBase ÷= UnitQuaternionBase
- *
- */
-
 
 // TransformBase × ColumnVector
 md_impl_all!(

src/geometry/translation.rs

@@ -14,8 +14,10 @@ pub type OwnedTranslation<N, D, S> = TranslationBase<N, D, Owned<N, D, U1, <S as
 
 /// A translation.
 #[repr(C)]
-#[derive(Hash, Debug, Clone, Copy)]
+#[derive(Hash, Debug, Clone, Copy, Serialize, Deserialize)]
 pub struct TranslationBase<N: Scalar, D: DimName, S/*: Storage<N, D, U1>*/> {
+    /// The translation coordinates, i.e., how much is added to a point's coordinates when it is
+    /// translated.
     pub vector: ColumnVector<N, D, S>
 }
 

src/geometry/unit_complex.rs

@@ -1,12 +1,11 @@
 use std::fmt;
+use approx::ApproxEq;
 use num_complex::Complex;
 
 use alga::general::Real;
-use core::{Unit, SquareMatrix, Vector1};
-use core::dimension::{U2, U3};
-use core::allocator::{OwnedAllocator, Allocator};
-use core::storage::OwnedStorage;
-use geometry::{RotationBase, OwnedRotation};
+use core::{Unit, SquareMatrix, Vector1, Matrix3};
+use core::dimension::U2;
+use geometry::Rotation2;
 
 /// A complex number with a norm equal to 1.
 pub type UnitComplex<N> = Unit<Complex<N>>;
@@ -15,7 +14,7 @@ impl<N: Real> UnitComplex<N> {
     /// The rotation angle in `]-pi; pi]` of this unit complex number.
     #[inline]
     pub fn angle(&self) -> N {
-        self.complex().im.atan2(self.complex().re)
+        self.im.atan2(self.re)
     }
 
     /// The rotation angle returned as a 1-dimensional vector.
@@ -35,7 +34,7 @@ impl<N: Real> UnitComplex<N> {
     /// Compute the conjugate of this unit complex number.
     #[inline]
     pub fn conjugate(&self) -> Self {
-        UnitComplex::new_unchecked(self.as_ref().conj())
+        UnitComplex::new_unchecked(self.conj())
     }
 
     /// Inverts this complex number if it is not zero.
@@ -83,13 +82,11 @@ impl<N: Real> UnitComplex<N> {
 
     /// Builds the rotation matrix corresponding to this unit complex number.
     #[inline]
-    pub fn to_rotation_matrix<S>(&self) -> RotationBase<N, U2, S>
-        where S: OwnedStorage<N, U2, U2>,
-              S::Alloc: OwnedAllocator<N, U2, U2, S> {
-        let r = self.complex().re;
-        let i = self.complex().im;
+    pub fn to_rotation_matrix(&self) -> Rotation2<N> {
+        let r = self.re;
+        let i = self.im;
 
-        RotationBase::from_matrix_unchecked(
+        Rotation2::from_matrix_unchecked(
             SquareMatrix::<_, U2, _>::new(
                 r, -i,
                 i,  r
@@ -99,12 +96,8 @@ impl<N: Real> UnitComplex<N> {
 
     /// Converts this unit complex number into its equivalent homogeneous transformation matrix.
     #[inline]
-    pub fn to_homogeneous<S>(&self) -> SquareMatrix<N, U3, S>
-        where S: OwnedStorage<N, U3, U3>,
-              S::Alloc: OwnedAllocator<N, U3, U3, S> +
-                        Allocator<N, U2, U2> {
-        let r: OwnedRotation<N, U2, S::Alloc> = self.to_rotation_matrix();
-        r.to_homogeneous()
+    pub fn to_homogeneous(&self) -> Matrix3<N> {
+        self.to_rotation_matrix().to_homogeneous()
     }
 }
 
@@ -113,3 +106,34 @@ impl<N: Real + fmt::Display> fmt::Display for UnitComplex<N> {
         write!(f, "UnitComplex angle: {}", self.angle())
     }
 }
+
+impl<N: Real> ApproxEq for UnitComplex<N> {
+    type Epsilon = N;
+
+    #[inline]
+    fn default_epsilon() -> Self::Epsilon {
+        N::default_epsilon()
+    }
+
+    #[inline]
+    fn default_max_relative() -> Self::Epsilon {
+        N::default_max_relative()
+    }
+
+    #[inline]
+    fn default_max_ulps() -> u32 {
+        N::default_max_ulps()
+    }
+
+    #[inline]
+    fn relative_eq(&self, other: &Self, epsilon: Self::Epsilon, max_relative: Self::Epsilon) -> bool {
+        self.re.relative_eq(&other.re, epsilon, max_relative) &&
+        self.im.relative_eq(&other.im, epsilon, max_relative)
+    }
+
+    #[inline]
+    fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool {
+        self.re.ulps_eq(&other.re, epsilon, max_ulps) &&
+        self.im.ulps_eq(&other.im, epsilon, max_ulps)
+    }
+}

src/geometry/unit_complex_conversion.rs

@@ -0,0 +1,173 @@
+use num::Zero;
+use num_complex::Complex;
+
+use alga::general::{SubsetOf, SupersetOf, Real};
+use alga::linear::Rotation as AlgaRotation;
+
+use core::SquareMatrix;
+use core::dimension::{U1, U2, U3};
+use core::storage::OwnedStorage;
+use core::allocator::{Allocator, OwnedAllocator};
+use geometry::{PointBase, UnitComplex, RotationBase, OwnedRotation, IsometryBase,
+               SimilarityBase, TransformBase, SuperTCategoryOf, TAffine, TranslationBase};
+
+/*
+ * This file provides the following conversions:
+ * =============================================
+ *
+ * UnitComplex -> UnitComplex
+ * UnitComplex -> RotationBase<U1>
+ * UnitComplex -> IsometryBase<U2>
+ * UnitComplex -> SimilarityBase<U2>
+ * UnitComplex -> TransformBase<U2>
+ * UnitComplex -> Matrix<U3> (homogeneous)
+ *
+ * NOTE:
+ * UnitComplex -> Complex is already provided by: Unit<T> -> T
+ */
+
+impl<N1, N2> SubsetOf<UnitComplex<N2>> for UnitComplex<N1>
+    where N1: Real,
+          N2: Real + SupersetOf<N1> {
+    #[inline]
+    fn to_superset(&self) -> UnitComplex<N2> {
+        UnitComplex::new_unchecked(self.as_ref().to_superset())
+    }
+
+    #[inline]
+    fn is_in_subset(uq: &UnitComplex<N2>) -> bool {
+        ::is_convertible::<_, Complex<N1>>(uq.as_ref())
+    }
+
+    #[inline]
+    unsafe fn from_superset_unchecked(uq: &UnitComplex<N2>) -> Self {
+        Self::new_unchecked(::convert_ref_unchecked(uq.as_ref()))
+    }
+}
+
+impl<N1, N2, S> SubsetOf<RotationBase<N2, U2, S>> for UnitComplex<N1>
+    where N1: Real,
+          N2: Real + SupersetOf<N1>,
+          S: OwnedStorage<N2, U2, U2>,
+          S::Alloc: OwnedAllocator<N2, U2, U2, S> +
+                    Allocator<N2, U3, U1> +
+                    Allocator<N2, U2, U1> +
+                    Allocator<N1, U2, U2> {
+    #[inline]
+    fn to_superset(&self) -> RotationBase<N2, U2, S> {
+        let q: UnitComplex<N2> = self.to_superset();
+        q.to_rotation_matrix().to_superset()
+    }
+
+    #[inline]
+    fn is_in_subset(rot: &RotationBase<N2, U2, S>) -> bool {
+        ::is_convertible::<_, OwnedRotation<N1, U2, S::Alloc>>(rot)
+    }
+
+    #[inline]
+    unsafe fn from_superset_unchecked(rot: &RotationBase<N2, U2, S>) -> Self {
+        let q = UnitComplex::<N2>::from_rotation_matrix(rot);
+        ::convert_unchecked(q)
+    }
+}
+
+
+impl<N1, N2, S, R> SubsetOf<IsometryBase<N2, U2, S, R>> for UnitComplex<N1>
+    where N1: Real,
+          N2: Real + SupersetOf<N1>,
+          S: OwnedStorage<N2, U2, U1>,
+          R: AlgaRotation<PointBase<N2, U2, S>> + SupersetOf<UnitComplex<N1>>,
+          S::Alloc: OwnedAllocator<N2, U2, U1, S> {
+    #[inline]
+    fn to_superset(&self) -> IsometryBase<N2, U2, S, R> {
+        IsometryBase::from_parts(TranslationBase::identity(), ::convert_ref(self))
+    }
+
+    #[inline]
+    fn is_in_subset(iso: &IsometryBase<N2, U2, S, R>) -> bool {
+        iso.translation.vector.is_zero()
+    }
+
+    #[inline]
+    unsafe fn from_superset_unchecked(iso: &IsometryBase<N2, U2, S, R>) -> Self {
+        ::convert_ref_unchecked(&iso.rotation)
+    }
+}
+
+
+impl<N1, N2, S, R> SubsetOf<SimilarityBase<N2, U2, S, R>> for UnitComplex<N1>
+    where N1: Real,
+          N2: Real + SupersetOf<N1>,
+          S: OwnedStorage<N2, U2, U1>,
+          R: AlgaRotation<PointBase<N2, U2, S>> + SupersetOf<UnitComplex<N1>>,
+          S::Alloc: OwnedAllocator<N2, U2, U1, S> {
+    #[inline]
+    fn to_superset(&self) -> SimilarityBase<N2, U2, S, R> {
+        SimilarityBase::from_isometry(::convert_ref(self), N2::one())
+    }
+
+    #[inline]
+    fn is_in_subset(sim: &SimilarityBase<N2, U2, S, R>) -> bool {
+        sim.isometry.translation.vector.is_zero() &&
+        sim.scaling() == N2::one()
+    }
+
+    #[inline]
+    unsafe fn from_superset_unchecked(sim: &SimilarityBase<N2, U2, S, R>) -> Self {
+        ::convert_ref_unchecked(&sim.isometry)
+    }
+}
+
+
+impl<N1, N2, S, C> SubsetOf<TransformBase<N2, U2, S, C>> for UnitComplex<N1>
+    where N1: Real,
+          N2: Real + SupersetOf<N1>,
+          S: OwnedStorage<N2, U3, U3>,
+          C: SuperTCategoryOf<TAffine>,
+          S::Alloc: OwnedAllocator<N2, U3, U3, S> +
+                    Allocator<N2, U2, U2>         +
+                    Allocator<N2, U1, U2>         +
+                    Allocator<N1, U2, U2>         +
+                    Allocator<N1, U3, U3> {
+    #[inline]
+    fn to_superset(&self) -> TransformBase<N2, U2, S, C> {
+        TransformBase::from_matrix_unchecked(self.to_homogeneous().to_superset())
+    }
+
+    #[inline]
+    fn is_in_subset(t: &TransformBase<N2, U2, S, C>) -> bool {
+        <Self as SubsetOf<_>>::is_in_subset(t.matrix())
+    }
+
+    #[inline]
+    unsafe fn from_superset_unchecked(t: &TransformBase<N2, U2, S, C>) -> Self {
+        Self::from_superset_unchecked(t.matrix())
+    }
+}
+
+
+impl<N1, N2, S> SubsetOf<SquareMatrix<N2, U3, S>> for UnitComplex<N1>
+    where N1: Real,
+          N2: Real + SupersetOf<N1>,
+          S: OwnedStorage<N2, U3, U3>,
+          S::Alloc: OwnedAllocator<N2, U3, U3, S> +
+                    Allocator<N2, U2, U2>         +
+                    Allocator<N2, U1, U2>         +
+                    Allocator<N1, U2, U2>         +
+                    Allocator<N1, U3, U3> {
+    #[inline]
+    fn to_superset(&self) -> SquareMatrix<N2, U3, S> {
+        self.to_homogeneous().to_superset()
+    }
+
+    #[inline]
+    fn is_in_subset(m: &SquareMatrix<N2, U3, S>) -> bool {
+        ::is_convertible::<_, OwnedRotation<N1, U2, S::Alloc>>(m)
+    }
+
+    #[inline]
+    unsafe fn from_superset_unchecked(m: &SquareMatrix<N2, U3, S>) -> Self {
+        let rot: OwnedRotation<N1, U2, S::Alloc> = ::convert_ref_unchecked(m);
+        Self::from_rotation_matrix(&rot)
+    }
+}

src/geometry/unit_complex_ops.rs

@@ -3,8 +3,13 @@ use std::ops::{Mul, MulAssign, Div, DivAssign};
 use alga::general::Real;
 use core::{Unit, ColumnVector, OwnedColumnVector};
 use core::dimension::{U1, U2};
-use core::storage::Storage;
-use geometry::{UnitComplex, RotationBase, PointBase, OwnedPoint};
+use core::storage::{Storage, OwnedStorage};
+use core::allocator::OwnedAllocator;
+use geometry::{UnitComplex, RotationBase,
+               PointBase, OwnedPoint,
+               IsometryBase, OwnedIsometryBase,
+               SimilarityBase, OwnedSimilarityBase,
+               TranslationBase};
 
 /*
  * This file provides:
@@ -23,8 +28,11 @@ use geometry::{UnitComplex, RotationBase, PointBase, OwnedPoint};
  * UnitComplex × ColumnVector
  * UnitComplex × Unit<T>
  *
- * NOTE: -UnitComplex is already provided by `Unit<T>`.
+ * UnitComplex × IsometryBase<UnitComplex>
+ * UnitComplex × SimilarityBase<UnitComplex>
+ * UnitComplex × TranslationBase -> IsometryBase<UnitComplex>
  *
+ * NOTE: -UnitComplex is already provided by `Unit<T>`.
  *
  * (Assignment Operators)
  *
@@ -34,8 +42,8 @@ use geometry::{UnitComplex, RotationBase, PointBase, OwnedPoint};
  * UnitComplex  ÷= UnitComplex
  * UnitComplex  ÷= RotationBase
  *
- * FIXME: RotationBase ×= UnitComplex
- * FIXME: RotationBase ÷= UnitComplex
+ * RotationBase ×= UnitComplex
+ * RotationBase ÷= UnitComplex
  *
  */
 
@@ -232,7 +240,7 @@ complex_op_impl_all!(
     [ref ref] => {
         let i = self.as_ref().im;
         let r = self.as_ref().re;
-        OwnedColumnVector::<_, U2, S::Alloc>::new(r * rhs[0] - i * rhs[0], i * rhs[1] + r * rhs[1])
+        OwnedColumnVector::<_, U2, S::Alloc>::new(r * rhs[0] - i * rhs[1], i * rhs[0] + r * rhs[1])
     };
 );
 
@@ -247,6 +255,45 @@ complex_op_impl_all!(
     [ref ref] => Unit::new_unchecked(self * rhs.as_ref());
 );
 
+// UnitComplex × IsometryBase<UnitComplex>
+complex_op_impl_all!(
+    Mul, mul;
+    (U2, U1);
+    self: UnitComplex<N>, rhs: IsometryBase<N, U2, S, UnitComplex<N>>,
+    Output = OwnedIsometryBase<N, U2, S::Alloc, UnitComplex<N>>;
+    [val val] => &self * &rhs;
+    [ref val] =>  self * &rhs;
+    [val ref] => &self *  rhs;
+    [ref ref] => {
+        let shift = self * &rhs.translation.vector;
+        IsometryBase::from_parts(TranslationBase::from_vector(shift), self * &rhs.rotation)
+    };
+);
+
+// UnitComplex × SimilarityBase<UnitComplex>
+complex_op_impl_all!(
+    Mul, mul;
+    (U2, U1);
+    self: UnitComplex<N>, rhs: SimilarityBase<N, U2, S, UnitComplex<N>>,
+    Output = OwnedSimilarityBase<N, U2, S::Alloc, UnitComplex<N>>;
+    [val val] => &self * &rhs;
+    [ref val] =>  self * &rhs;
+    [val ref] => &self *  rhs;
+    [ref ref] => SimilarityBase::from_isometry(self * &rhs.isometry, rhs.scaling());
+);
+
+// UnitComplex × TranslationBase
+complex_op_impl_all!(
+    Mul, mul;
+    (U2, U1);
+    self: UnitComplex<N>, rhs: TranslationBase<N, U2, S>,
+    Output = OwnedIsometryBase<N, U2, S::Alloc, UnitComplex<N>>;
+    [val val] => IsometryBase::from_parts(TranslationBase::from_vector(&self *  rhs.vector), self);
+    [ref val] => IsometryBase::from_parts(TranslationBase::from_vector( self *  rhs.vector), self.clone());
+    [val ref] => IsometryBase::from_parts(TranslationBase::from_vector(&self * &rhs.vector), self);
+    [ref ref] => IsometryBase::from_parts(TranslationBase::from_vector( self * &rhs.vector), self.clone());
+);
+
 // UnitComplex ×= UnitComplex
 impl<N: Real> MulAssign<UnitComplex<N>> for UnitComplex<N> {
     #[inline]
@@ -307,3 +354,42 @@ impl<'b, N: Real, S: Storage<N, U2, U2>> DivAssign<&'b RotationBase<N, U2, S>> f
         *self = &*self / rhs
     }
 }
+
+
+// RotationBase ×= UnitComplex
+impl<N: Real, S> MulAssign<UnitComplex<N>> for RotationBase<N, U2, S>
+    where S: OwnedStorage<N, U2, U2>,
+          S::Alloc: OwnedAllocator<N, U2, U2, S> {
+    #[inline]
+    fn mul_assign(&mut self, rhs: UnitComplex<N>) {
+        self.mul_assign(rhs.to_rotation_matrix())
+    }
+}
+
+impl<'b, N: Real, S: Storage<N, U2, U2>> MulAssign<&'b UnitComplex<N>> for RotationBase<N, U2, S>
+    where S: OwnedStorage<N, U2, U2>,
+          S::Alloc: OwnedAllocator<N, U2, U2, S> {
+    #[inline]
+    fn mul_assign(&mut self, rhs: &'b UnitComplex<N>) {
+        self.mul_assign(rhs.to_rotation_matrix())
+    }
+}
+
+// RotationBase ÷= UnitComplex
+impl<N: Real, S> DivAssign<UnitComplex<N>> for RotationBase<N, U2, S>
+    where S: OwnedStorage<N, U2, U2>,
+          S::Alloc: OwnedAllocator<N, U2, U2, S> {
+    #[inline]
+    fn div_assign(&mut self, rhs: UnitComplex<N>) {
+        self.div_assign(rhs.to_rotation_matrix())
+    }
+}
+
+impl<'b, N: Real, S: Storage<N, U2, U2>> DivAssign<&'b UnitComplex<N>> for RotationBase<N, U2, S>
+    where S: OwnedStorage<N, U2, U2>,
+          S::Alloc: OwnedAllocator<N, U2, U2, S> {
+    #[inline]
+    fn div_assign(&mut self, rhs: &'b UnitComplex<N>) {
+        self.div_assign(rhs.to_rotation_matrix())
+    }
+}

src/lib.rs

@@ -1,14 +1,97 @@
+/*!
+# nalgebra
+
+**nalgebra** is a linear algebra library written for Rust targeting:
+
+* General-purpose linear algebra (still lacks a lot of features…)
+* Real time computer graphics.
+* Real time computer physics.
+
+## Using **nalgebra**
+You will need the last stable build of the [rust compiler](http://www.rust-lang.org)
+and the official package manager: [cargo](https://github.com/rust-lang/cargo).
+
+Simply add the following to your `Cargo.toml` file:
+
+```.ignore
+[dependencies]
+nalgebra = "0.11"
+```
+
+
+Most useful functionalities of **nalgebra** are grouped in the root module `nalgebra::`.
+
+However, the recommended way to use **nalgebra** is to import types and traits
+explicitly, and call free-functions using the `na::` prefix:
+
+```.rust
+#[macro_use]
+extern crate approx; // For the macro relative_eq!
+extern crate nalgebra as na;
+use na::{Vector3, Rotation3};
+
+fn main() {
+    let axis  = Vector3::x_axis();
+    let angle = 1.57;
+    let b     = Rotation3::from_axis_angle(&axis, angle);
+
+    relative_eq!(b.axis().unwrap(), axis);
+    relative_eq!(b.angle(), angle);
+}
+```
+
+
+## Features
+**nalgebra** is meant to be a general-purpose, low-dimensional, linear algebra library, with
+an optimized set of tools for computer graphics and physics. Those features include:
+
+* A single parametrizable type `Matrix` for vectors, (square or rectangular) matrices, and slices
+  with dimensions known either at compile-time (using type-level integers) or at runtime.
+* Matrices and vectors with compile-time sizes are statically allocated while dynamic ones are
+  allocated on the heap.
+* Convenient aliases for low-dimensional matrices and vectors: `Vector1` to `Vector6` and
+  `Matrix1x1` to `Matrix6x6` (including rectangular matrices like `Matrix2x5`.
+* Points sizes known at compile time, and convenience aliases: `Point1` to `Point6`.
+* Translation (seen as a transformation that composes by multiplication): `Translation2`,
+  `Translation3`.
+* Rotation matrices: `Rotation2`, `Rotation3`.
+* Quaternions: `Quaternion`, `UnitQuaternion` (for 3D rotation).
+* Unit complex numbers can be used for 2D rotation: `UnitComplex`.
+* Algebraic entities with a norm equal to one: `Unit<T>`, e.g., `Unit<Vector3<f32>>`.
+* Isometries (translation ⨯ rotation): `Isometry2`, `Isometry3`
+* Similarity transformations (translation ⨯ rotation ⨯ uniform scale): `Similarity2`, `Similarity3`.
+* Affine transformations stored as an homogeneous matrix: `Affine2`, `Affine3`.
+* Projective (i.e. invertible) transformations stored as an homogeneous matrix: `Projective2`,
+  `Projective3`.
+* General transformations that does not have to be invertible, stored as an homogeneous matrix:
+  `Transform2`, `Transform3`.
+* 3D projections for computer graphics: `Perspective3`, `Orthographic3`.
+* Linear algebra and data analysis operators: QR decomposition, eigen-decomposition.
+* Implements all meaningful traits from the [alga](https://crates.io/crates/alga) crate for
+  generic programming.
+*/
+
+
+
+
+
 // #![feature(plugin)]
 //
 // #![plugin(clippy)]
-// #![allow(too_many_arguments)]
-// #![allow(derive_hash_xor_eq)]
-// #![allow(len_without_is_empty)]
-// #![allow(transmute_ptr_to_ref)]
+
+#![deny(non_camel_case_types)]
+#![deny(unused_parens)]
+#![deny(non_upper_case_globals)]
+#![deny(unused_qualifications)]
+#![deny(unused_results)]
+#![warn(missing_docs)]
+#![doc(html_root_url = "http://nalgebra.org/doc")]
 
 #[cfg(feature = "arbitrary")]
 extern crate quickcheck;
-extern crate rustc_serialize;
+extern crate serde;
+#[macro_use]
+extern crate serde_derive;
 extern crate num_traits as num;
 extern crate num_complex;
 extern crate rand;
@@ -32,11 +115,14 @@ pub use traits::*;
 use std::cmp::{self, PartialOrd, Ordering};
 
 use num::Signed;
-use alga::general::{Id, Identity, SupersetOf, MeetSemilattice, JoinSemilattice, Lattice, Inverse,
+use alga::general::{Identity, SupersetOf, MeetSemilattice, JoinSemilattice, Lattice, Inverse,
                     Multiplicative, Additive, AdditiveGroup};
 use alga::linear::SquareMatrix as AlgaSquareMatrix;
 use alga::linear::{InnerSpace, NormedSpace, FiniteDimVectorSpace, EuclideanSpace};
 
+pub use alga::general::Id;
+
+
 /*
  *
  * Multiplicative identity.
@@ -279,7 +365,7 @@ pub fn try_inverse<M: AlgaSquareMatrix>(m: &M) -> Option<M> {
     m.try_inverse()
 }
 
-/// Computes the the multiplicative inverse (transformation, unit quaternion, etc.)
+/// Computes the multiplicative inverse of an (always invertible) algebraic entity.
 #[inline]
 pub fn inverse<M: Inverse<Multiplicative>>(m: &M) -> M {
     m.inverse()
@@ -295,7 +381,7 @@ pub fn dot<V: FiniteDimVectorSpace>(a: &V, b: &V) -> V::Field {
     a.dot(b)
 }
 
-/// Computes the angle between two vectors.
+/// Computes the smallest angle between two vectors.
 #[inline]
 pub fn angle<V: InnerSpace>(a: &V, b: &V) -> V::Real {
     a.angle(b)
@@ -305,25 +391,25 @@ pub fn angle<V: InnerSpace>(a: &V, b: &V) -> V::Real {
  * Normed space
  */
 
-/// Computes the L2 norm of a vector.
+/// Computes the L2 (euclidean) norm of a vector.
 #[inline]
 pub fn norm<V: NormedSpace>(v: &V) -> V::Field {
     v.norm()
 }
 
-/// Computes the squared L2 norm of a vector.
+/// Computes the squared L2 (euclidean) norm of the vector `v`.
 #[inline]
 pub fn norm_squared<V: NormedSpace>(v: &V) -> V::Field {
     v.norm_squared()
 }
 
-/// Gets the normalized version of a vector.
+/// Computes the normalized version of the vector `v`.
 #[inline]
 pub fn normalize<V: NormedSpace>(v: &V) -> V {
     v.normalize()
 }
 
-/// Gets the normalized version of a vector or `None` if its norm is smaller than `min_norm`.
+/// Computes the normalized version of the vector `v` or returns `None` if its norm is smaller than `min_norm`.
 #[inline]
 pub fn try_normalize<V: NormedSpace>(v: &V, min_norm: V::Field) -> Option<V> {
     v.try_normalize(min_norm)

tests/isometry.rs

@@ -10,7 +10,10 @@ extern crate alga;
 extern crate nalgebra as na;
 
 use alga::linear::{Transformation, ProjectiveTransformation};
-use na::{Vector3, Point3, Rotation3, Isometry3, Translation3, UnitQuaternion};
+use na::{
+    Vector3, Point3, Rotation3, Isometry3, Translation3, UnitQuaternion,
+    Vector2, Point2, Rotation2, Isometry2, Translation2, UnitComplex
+};
 
 quickcheck!(
     fn append_rotation_wrt_point_to_id(r: UnitQuaternion<f64>, p: Point3<f64>) -> bool {
@@ -65,7 +68,38 @@ quickcheck!(
         relative_eq!(i.inverse() * p, i.inverse_transform_point(&p), epsilon = 1.0e-7)
     }
 
-    fn composition(i: Isometry3<f64>, uq: UnitQuaternion<f64>, r: Rotation3<f64>,
+    fn composition2(i: Isometry2<f64>,    uc: UnitComplex<f64>, r: Rotation2<f64>,
+                    t: Translation2<f64>, v:  Vector2<f64>,     p: Point2<f64>) -> bool {
+        // (rotation × translation) * point = rotation × (translation * point)
+        relative_eq!((uc * t) * v, uc * v, epsilon = 1.0e-7)       &&
+        relative_eq!((r  * t) * v, r  * v, epsilon = 1.0e-7)       &&
+        relative_eq!((uc * t) * p, uc * (t * p), epsilon = 1.0e-7) &&
+        relative_eq!((r * t)  * p, r  * (t * p), epsilon = 1.0e-7) &&
+
+        // (translation × rotation) * point = translation × (rotation * point)
+        (t * uc) * v == uc * v       &&
+        (t * r)  * v == r  * v       &&
+        (t * uc) * p == t * (uc * p) &&
+        (t * r)  * p == t * (r  * p) &&
+
+        // (rotation × isometry) * point = rotation × (isometry * point)
+        relative_eq!((uc * i) * v, uc * (i * v), epsilon = 1.0e-7) &&
+        relative_eq!((uc * i) * p, uc * (i * p), epsilon = 1.0e-7) &&
+
+        // (isometry × rotation) * point = isometry × (rotation * point)
+        relative_eq!((i * uc) * v, i * (uc * v), epsilon = 1.0e-7) &&
+        relative_eq!((i * uc) * p, i * (uc * p), epsilon = 1.0e-7) &&
+
+        // (translation × isometry) * point = translation × (isometry * point)
+        relative_eq!((t * i) * v,     (i * v), epsilon = 1.0e-7) &&
+        relative_eq!((t * i) * p, t * (i * p), epsilon = 1.0e-7) &&
+
+        // (isometry × translation) * point = isometry × (translation * point)
+        relative_eq!((i * t) * v, i * v,       epsilon = 1.0e-7) &&
+        relative_eq!((i * t) * p, i * (t * p), epsilon = 1.0e-7)
+    }
+
+    fn composition3(i: Isometry3<f64>,    uq: UnitQuaternion<f64>, r: Rotation3<f64>,
                     t: Translation3<f64>, v:  Vector3<f64>,        p: Point3<f64>) -> bool {
         // (rotation × translation) * point = rotation × (translation * point)
         relative_eq!((uq * t) * v, uq * v, epsilon = 1.0e-7)       &&

tests/matrix_slice.rs

@@ -1,8 +1,3 @@
-#[cfg(feature = "arbitrary")]
-#[macro_use]
-extern crate quickcheck;
-#[macro_use]
-extern crate approx;
 extern crate num_traits as num;
 extern crate nalgebra as na;
 

tests/projection.rs

@@ -0,0 +1,51 @@
+#[cfg(feature = "arbitrary")]
+#[macro_use]
+extern crate quickcheck;
+#[macro_use]
+extern crate approx;
+
+extern crate nalgebra as na;
+
+use na::{Point3, Perspective3, Orthographic3};
+
+#[test]
+fn perspective_inverse() {
+    let proj = Perspective3::new(800.0 / 600.0, 3.14 / 2.0, 1.0, 1000.0);
+    let inv  = proj.inverse();
+
+    let id = inv * proj.unwrap();
+
+    assert!(id.is_identity(1.0e-7));
+}
+
+#[test]
+fn orthographic_inverse() {
+    let proj = Orthographic3::new(1.0, 2.0, -3.0, -2.5, 10.0, 900.0);
+    let inv  = proj.inverse();
+
+    let id = inv * proj.unwrap();
+
+    assert!(id.is_identity(1.0e-7));
+}
+
+
+#[cfg(feature = "arbitrary")]
+quickcheck!{
+    fn perspective_project_unproject(pt: Point3<f64>) -> bool {
+        let proj = Perspective3::new(800.0 / 600.0, 3.14 / 2.0, 1.0, 1000.0);
+
+        let projected   = proj.project_point(&pt);
+        let unprojected = proj.unproject_point(&projected);
+
+        relative_eq!(pt, unprojected, epsilon = 1.0e-7)
+    }
+
+    fn orthographic_project_unproject(pt: Point3<f64>) -> bool {
+        let proj = Orthographic3::new(1.0, 2.0, -3.0, -2.5, 10.0, 900.0);
+
+        let projected   = proj.project_point(&pt);
+        let unprojected = proj.unproject_point(&projected);
+
+        relative_eq!(pt, unprojected, epsilon = 1.0e-7)
+    }
+}

tests/rotation.rs

@@ -7,6 +7,7 @@ extern crate num_traits as num;
 extern crate alga;
 extern crate nalgebra as na;
 
+use std::f64;
 use alga::general::Real;
 use na::{Vector2, Vector3, Rotation2, Rotation3, Unit};
 
@@ -67,7 +68,6 @@ quickcheck!(
         a.angle(&b) == b.angle(&a)
     }
 
-
     /*
      *
      * RotationBase matrix between vectors.

tests/serde.rs

@@ -0,0 +1,47 @@
+extern crate rand;
+extern crate nalgebra as na;
+
+extern crate serde_json;
+
+use na::{
+    DMatrix,
+    Matrix3x4,
+    Point3,
+    Translation3,
+    Rotation3,
+    Isometry3,
+    IsometryMatrix3,
+    Similarity3,
+    SimilarityMatrix3,
+    Quaternion
+};
+
+macro_rules! test_serde(
+    ($($test: ident, $ty: ident);* $(;)*) => {$(
+        #[test]
+        fn $test() {
+            let v: $ty<f32> = rand::random();
+            let serialized = serde_json::to_string(&v).unwrap();
+            assert_eq!(v, serde_json::from_str(&serialized).unwrap());
+        }
+    )*}
+);
+
+#[test]
+fn serde_dmatrix() {
+    let v: DMatrix<f32> = DMatrix::new_random(3, 4);
+    let serialized = serde_json::to_string(&v).unwrap();
+    assert_eq!(v, serde_json::from_str(&serialized).unwrap());
+}
+
+test_serde!(
+    serde_matrix3x4,          Matrix3x4;
+    serde_point3,             Point3;
+    serde_translation3,       Translation3;
+    serde_rotation3,          Rotation3;
+    serde_isometry3,          Isometry3;
+    serde_isometry_matrix3,   IsometryMatrix3;
+    serde_similarity3,        Similarity3;
+    serde_similarity_matrix3, SimilarityMatrix3;
+    serde_quaternion,         Quaternion;
+);

tests/unit_complex.rs

@@ -0,0 +1,171 @@
+#![allow(non_snake_case)]
+
+#[cfg(feature = "arbitrary")]
+#[macro_use]
+extern crate quickcheck;
+#[macro_use]
+extern crate approx;
+extern crate num_traits as num;
+extern crate alga;
+extern crate nalgebra as na;
+
+use na::{Unit, UnitComplex, Vector2, Point2, Rotation2};
+
+
+quickcheck!(
+
+    /*
+     *
+     * From/to rotation matrix.
+     *
+     */
+    fn unit_complex_rotation_conversion(c: UnitComplex<f64>) -> bool {
+        let r  = c.to_rotation_matrix();
+        let cc = UnitComplex::from_rotation_matrix(&r);
+        let rr = cc.to_rotation_matrix();
+
+        relative_eq!(c, cc, epsilon = 1.0e-7) &&
+        relative_eq!(r, rr, epsilon = 1.0e-7)
+    }
+
+    /*
+     *
+     * PointBase/Vector transformation.
+     *
+     */
+    fn unit_complex_transformation(c: UnitComplex<f64>, v: Vector2<f64>, p: Point2<f64>) -> bool {
+        let r  = c.to_rotation_matrix();
+        let rv = r * v;
+        let rp = r * p;
+
+        relative_eq!( c *  v, rv, epsilon = 1.0e-7) &&
+        relative_eq!( c * &v, rv, epsilon = 1.0e-7) &&
+        relative_eq!(&c *  v, rv, epsilon = 1.0e-7) &&
+        relative_eq!(&c * &v, rv, epsilon = 1.0e-7) &&
+
+        relative_eq!( c *  p, rp, epsilon = 1.0e-7) &&
+        relative_eq!( c * &p, rp, epsilon = 1.0e-7) &&
+        relative_eq!(&c *  p, rp, epsilon = 1.0e-7) &&
+        relative_eq!(&c * &p, rp, epsilon = 1.0e-7)
+    }
+
+    /*
+     *
+     * Inversion.
+     *
+     */
+    fn unit_complex_inv(c: UnitComplex<f64>) -> bool {
+        let iq = c.inverse();
+        relative_eq!(&iq * &c, UnitComplex::identity(), epsilon = 1.0e-7) &&
+        relative_eq!( iq * &c, UnitComplex::identity(), epsilon = 1.0e-7) &&
+        relative_eq!(&iq *  c, UnitComplex::identity(), epsilon = 1.0e-7) &&
+        relative_eq!( iq *  c, UnitComplex::identity(), epsilon = 1.0e-7) &&
+
+        relative_eq!(&c * &iq, UnitComplex::identity(), epsilon = 1.0e-7) &&
+        relative_eq!( c * &iq, UnitComplex::identity(), epsilon = 1.0e-7) &&
+        relative_eq!(&c *  iq, UnitComplex::identity(), epsilon = 1.0e-7) &&
+        relative_eq!( c *  iq, UnitComplex::identity(), epsilon = 1.0e-7)
+    }
+
+    /*
+     *
+     * Quaterion * Vector == RotationBase * Vector
+     *
+     */
+    fn unit_complex_mul_vector(c: UnitComplex<f64>, v: Vector2<f64>, p: Point2<f64>) -> bool {
+        let r = c.to_rotation_matrix();
+
+        relative_eq!(c * v, r * v, epsilon = 1.0e-7) &&
+        relative_eq!(c * p, r * p, epsilon = 1.0e-7)
+    }
+
+    // Test that all operators (incl. all combinations of references) work.
+    // See the top comment on `geometry/quaternion_ops.rs` for details on which operations are
+    // supported.
+    fn all_op_exist(uc: UnitComplex<f64>, v: Vector2<f64>, p: Point2<f64>, r: Rotation2<f64>) -> bool {
+        let uv = Unit::new_normalize(v);
+
+        let ucMuc = uc * uc;
+        let ucMr  = uc * r;
+        let rMuc  = r  * uc;
+        let ucDuc = uc / uc;
+        let ucDr  = uc / r;
+        let rDuc  = r / uc;
+
+        let ucMp  = uc * p;
+        let ucMv  = uc * v;
+        let ucMuv = uc * uv;
+
+        let mut ucMuc1 = uc;
+        let mut ucMuc2 = uc;
+
+        let mut ucMr1 = uc;
+        let mut ucMr2 = uc;
+
+        let mut ucDuc1 = uc;
+        let mut ucDuc2 = uc;
+
+        let mut ucDr1 = uc;
+        let mut ucDr2 = uc;
+
+        ucMuc1 *= uc;
+        ucMuc2 *= &uc;
+
+        ucMr1 *= r;
+        ucMr2 *= &r;
+
+        ucDuc1 /= uc;
+        ucDuc2 /= &uc;
+
+        ucDr1 /= r;
+        ucDr2 /= &r;
+
+        ucMuc1 == ucMuc  &&
+        ucMuc1 == ucMuc2 &&
+
+        ucMr1  == ucMr   &&
+        ucMr1  == ucMr2  &&
+
+        ucDuc1 == ucDuc  &&
+        ucDuc1 == ucDuc2 &&
+
+        ucDr1  == ucDr   &&
+        ucDr1  == ucDr2  &&
+
+        ucMuc == &uc * &uc &&
+        ucMuc ==  uc * &uc &&
+        ucMuc == &uc *  uc &&
+
+        ucMr == &uc * &r &&
+        ucMr ==  uc * &r &&
+        ucMr == &uc *  r &&
+
+        rMuc == &r * &uc &&
+        rMuc ==  r * &uc &&
+        rMuc == &r *  uc &&
+
+        ucDuc == &uc / &uc &&
+        ucDuc ==  uc / &uc &&
+        ucDuc == &uc /  uc &&
+
+        ucDr == &uc / &r &&
+        ucDr ==  uc / &r &&
+        ucDr == &uc /  r &&
+
+        rDuc == &r / &uc &&
+        rDuc ==  r / &uc &&
+        rDuc == &r /  uc &&
+
+        ucMp == &uc * &p &&
+        ucMp ==  uc * &p &&
+        ucMp == &uc *  p &&
+
+        ucMv == &uc * &v &&
+        ucMv ==  uc * &v &&
+        ucMv == &uc *  v &&
+
+        ucMuv == &uc * &uv &&
+        ucMuv ==  uc * &uv &&
+        ucMuv == &uc *  uv
+    }
+);