Add n-dimensional vector and matrix.

README.md

@@ -6,6 +6,10 @@ programming language.
 It is mainly focused on features needed for real-time physics. It should be
 usable for graphics too.
 
+## Disclaimer
+
+As of today, nalgebra is largely untested.
+
 ## Licence
 
 nalgebra is provided "as is", under the BSD 3-Clause License.

src/dim2/mat2.rs

@@ -1,7 +1,10 @@
 use core::num::{One, Zero};
+use traits::dim::Dim;
 use traits::inv::Inv;
+use traits::transpose::Transpose;
 use dim2::vec2::Vec2;
 
+#[deriving(Eq)]
 pub struct Mat2<T>
 {
     m11: T, m12: T,
@@ -17,6 +20,12 @@ pub fn Mat2<T:Copy>(m11: T, m12: T, m21: T, m22: T) -> Mat2<T>
   }
 }
 
+impl<T> Dim for Mat2<T>
+{
+  fn dim() -> uint
+  { 2 }
+}
+
 impl<T:Copy + One + Zero> One for Mat2<T>
 {
   fn one() -> Mat2<T>
@@ -45,7 +54,7 @@ impl<T:Copy + Zero> Zero for Mat2<T>
 
 impl<T:Copy + Mul<T, T> + Add<T, T>> Mul<Mat2<T>, Mat2<T>> for Mat2<T>
 {
-  fn mul(&self, other : &Mat2<T>) -> Mat2<T>
+  fn mul(&self, other: &Mat2<T>) -> Mat2<T>
   {
     Mat2
     (self.m11 * other.m11 + self.m12 * other.m21,
@@ -60,18 +69,18 @@ impl<T:Copy + Mul<T, T> + Add<T, T>> Mul<Mat2<T>, Mat2<T>> for Mat2<T>
 //
 // impl<T:Copy + Mul<T, T> + Add<T, T>> Mul<Vec2<T>, Vec2<T>> for Mat2<T>
 // {
-//   fn mul(&self, v : &Vec2<T>) -> Vec2<T>
+//   fn mul(&self, v: &Vec2<T>) -> Vec2<T>
 //   { Vec2(self.m11 * v.x + self.m12 * v.y, self.m21 * v.x + self.m22 * v.y) }
 // }
 
 impl<T:Copy + Mul<T, T> + Add<T, T>> Mul<Mat2<T>, Vec2<T>> for Vec2<T>
 {
-  fn mul(&self, m : &Mat2<T>) -> Vec2<T>
+  fn mul(&self, m: &Mat2<T>) -> Vec2<T>
   { Vec2(self.x * m.m11 + self.y * m.m21, self.x * m.m12 + self.y * m.m22) }
 }
 
 impl<T:Copy + Mul<T, T> + Div<T, T> + Sub<T, T> + Neg<T> + Eq + Zero>
-Inv<T> for Mat2<T>
+Inv for Mat2<T>
 {
   fn inv(&self) -> Mat2<T>
   {
@@ -84,6 +93,18 @@ Inv<T> for Mat2<T>
   }
 }
 
+impl<T:Copy> Transpose for Mat2<T>
+{
+  fn transposed(&self) -> Mat2<T>
+  {
+    Mat2(self.m11, self.m21,
+         self.m12, self.m22)
+  }
+
+  fn transpose(&mut self)
+  { self.m21 <-> self.m12; }
+}
+
 impl<T:ToStr> ToStr for Mat2<T>
 {
   fn to_str(&self) -> ~str

src/dim2/rotmat2.rs

@@ -9,7 +9,7 @@ use dim2::mat2::Mat2;
 #[deriving(Eq)]
 pub struct Rotmat2<T>
 {
-  submat: Mat2<T>
+  priv submat: Mat2<T>
 }
 
 pub fn Rotmat2<T:Copy + Trigonometric + Neg<T>>(angle: T) -> Rotmat2<T>

src/dim2/vec2.rs

@@ -1,5 +1,8 @@
+use core::num::Zero;
 use traits::dot::Dot;
-use traits::sqrt::Sqrt;
+use traits::dim::Dim;
+use traits::cross::Cross;
+use traits::workarounds::sqrt::Sqrt;
 
 #[deriving(Eq)]
 pub struct Vec2<T>
@@ -11,6 +14,11 @@ pub struct Vec2<T>
 pub fn Vec2<T:Copy>(x: T, y: T) -> Vec2<T>
 { Vec2 {x: x, y: y} }
 
+impl<T> Dim for Vec2<T>
+{
+  fn dim() -> uint
+  { 2 }
+}
 
 impl<T:Copy + Add<T,T>> Add<Vec2<T>, Vec2<T>> for Vec2<T>
 {
@@ -36,6 +44,21 @@ impl<T:Copy + Mul<T, T> + Add<T, T> + Sqrt> Dot<T> for Vec2<T>
   { self.sqnorm().sqrt() }
 }
 
+impl<T:Copy + Mul<T, T> + Sub<T, T>> Cross<T> for Vec2<T>
+{
+  fn cross(&self, other : &Vec2<T>) -> T
+  { self.x * other.y - self.y * other.x }
+}
+
+impl<T:Copy + Zero> Zero for Vec2<T>
+{
+  fn zero() -> Vec2<T>
+  {
+    let _0 = Zero::zero();
+    Vec2(_0, _0)
+  }
+}
+
 impl<T:ToStr> ToStr for Vec2<T>
 {
   fn to_str(&self) -> ~str

src/dim3/mat3.rs

@@ -1,7 +1,11 @@
 use core::num::{One, Zero};
+// use core::rand::{Rand, Rng};
+use traits::dim::Dim;
 use traits::inv::Inv;
+use traits::transpose::Transpose;
 use dim3::vec3::Vec3;
 
+#[deriving(Eq)]
 pub struct Mat3<T>
 {
     m11: T, m12: T, m13: T,
@@ -21,6 +25,12 @@ pub fn Mat3<T:Copy>(m11: T, m12: T, m13: T,
   }
 }
 
+impl<T> Dim for Mat3<T>
+{
+  fn dim() -> uint
+  { 3 }
+}
+
 impl<T:Copy + One + Zero> One for Mat3<T>
 {
   fn one() -> Mat3<T>
@@ -95,7 +105,7 @@ impl<T:Copy + Mul<T, T> + Add<T, T>> Mul<Mat3<T>, Vec3<T>> for Vec3<T>
 
 impl<T:Copy + Mul<T, T> + Div<T, T> + Sub<T, T> + Add<T, T> + Neg<T>
      + Eq + Zero>
-Inv<T> for Mat3<T>
+Inv for Mat3<T>
 {
   fn inv(&self) -> Mat3<T>
   {
@@ -125,6 +135,36 @@ Inv<T> for Mat3<T>
   }
 }
 
+impl<T:Copy> Transpose for Mat3<T>
+{
+  fn transposed(&self) -> Mat3<T>
+  {
+    Mat3(self.m11, self.m21, self.m31,
+         self.m12, self.m22, self.m32,
+         self.m13, self.m23, self.m33)
+  }
+
+  fn transpose(&mut self)
+  {
+    self.m12 <-> self.m21;
+    self.m13 <-> self.m31;
+    self.m23 <-> self.m32;
+  }
+}
+
+// impl<T:Rand> Rand for Mat3<T>
+// {
+//   fn rand<R:Rng>(rng: &mut R) -> Mat3<T>
+//   {
+//     Mat3
+//     {
+//       m11: rng.next(), m12: rng.next(), m13: rng.next(),
+//       m21: rng.next(), m22: rng.next(), m23: rng.next(),
+//       m31: rng.next(), m32: rng.next(), m33: rng.next()
+//     }
+//   }
+// }
+
 impl<T:ToStr> ToStr for Mat3<T>
 {
   fn to_str(&self) -> ~str

src/dim3/vec3.rs

@@ -1,5 +1,8 @@
+use core::num::Zero;
+use traits::dim::Dim;
 use traits::dot::Dot;
-use traits::sqrt::Sqrt;
+use traits::cross::Cross;
+use traits::workarounds::sqrt::Sqrt;
 
 #[deriving(Eq)]
 pub struct Vec3<T>
@@ -12,6 +15,11 @@ pub struct Vec3<T>
 pub fn Vec3<T:Copy>(x: T, y: T, z: T) -> Vec3<T>
 { Vec3 {x: x, y: y, z: z} }
 
+impl<T> Dim for Vec3<T>
+{
+  fn dim() -> uint
+  { 3 }
+}
 
 impl<T:Copy + Add<T,T>> Add<Vec3<T>, Vec3<T>> for Vec3<T>
 {
@@ -25,18 +33,6 @@ impl<T:Copy + Sub<T,T>> Sub<Vec3<T>, Vec3<T>> for Vec3<T>
   { Vec3(self.x - other.x, self.y - other.y, self.z - other.z) }
 }
 
-impl<T:ToStr> ToStr for Vec3<T>
-{
-  fn to_str(&self) -> ~str
-  {
-    ~"Vec3 "
-    + "{ x : " + self.x.to_str()
-    + ", y : " + self.y.to_str()
-    + ", z : " + self.z.to_str()
-    + " }"
-  }
-}
-
 impl<T:Copy + Mul<T, T> + Add<T, T> + Sqrt> Dot<T> for Vec3<T>
 {
   fn dot(&self, other : &Vec3<T>) -> T
@@ -48,3 +44,36 @@ impl<T:Copy + Mul<T, T> + Add<T, T> + Sqrt> Dot<T> for Vec3<T>
   fn norm(&self) -> T
   { self.sqnorm().sqrt() }
 }
+
+impl<T:Copy + Mul<T, T> + Sub<T, T>> Cross<Vec3<T>> for Vec3<T>
+{
+  fn cross(&self, other : &Vec3<T>) -> Vec3<T>
+  {
+    Vec3(
+      self.y * other.z - self.z * other.y,
+      self.z * other.x - self.x * other.z,
+      self.x * other.y - self.y * other.x
+    )
+  }
+}
+
+impl<T:Copy + Zero> Zero for Vec3<T>
+{
+  fn zero() -> Vec3<T>
+  {
+    let _0 = Zero::zero();
+    Vec3(_0, _0, _0)
+  }
+}
+
+impl<T:ToStr> ToStr for Vec3<T>
+{
+  fn to_str(&self) -> ~str
+  {
+    ~"Vec3 "
+    + "{ x : " + self.x.to_str()
+    + ", y : " + self.y.to_str()
+    + ", z : " + self.z.to_str()
+    + " }"
+  }
+}

src/nalgebra.rc

@@ -7,11 +7,16 @@
 
 extern mod std;
 
-pub use dim2::vec2::Vec2;
 pub use dim3::vec3::Vec3;
+pub use dim3::mat3::Mat3;
+
+pub use dim2::vec2::Vec2;
 pub use dim2::mat2::Mat2;
 pub use dim2::rotmat2::Rotmat2;
 
+pub use ndim::nvec::NVec;
+pub use ndim::nmat::NMat;
+
 pub use traits::dot::Dot;
 pub use traits::cross::Cross;
 pub use traits::dim::Dim;
@@ -37,6 +42,8 @@ mod dim3
 
 mod ndim
 {
+  mod nvec;
+  mod nmat;
 }
 
 mod traits

src/ndim/nmat.rs

@@ -0,0 +1,187 @@
+use core::num::{One, Zero};
+use core::vec::{from_elem, swap};
+use traits::dim::Dim;
+use traits::inv::Inv;
+use traits::transpose::Transpose;
+// use ndim::nvec::NVec;
+
+// D is a phantom type parameter, used only as a dimensional token.
+// Its allows use to encode the vector dimension at the type-level.
+// It can be anything implementing the Dim trait. However, to avoid confusion,
+// using d0, d1, d2, d3 and d4 tokens are prefered.
+#[deriving(Eq)]
+pub struct NMat<D, T>
+{
+  mij: ~[T]
+}
+
+impl<D: Dim, T: Copy> NMat<D, T>
+{
+  fn offset(i: uint, j: uint) -> uint
+  { i * Dim::dim::<D>() + j }
+
+  fn set(&mut self, i: uint, j: uint, t: &T)
+  { self.mij[NMat::offset::<D, T>(i, j)] = *t }
+}
+
+impl<D: Dim, T> Dim for NMat<D, T>
+{
+  fn dim() -> uint
+  { Dim::dim::<D>() }
+}
+
+impl<D: Dim, T:Copy> Index<(uint, uint), T> for NMat<D, T>
+{
+  fn index(&self, &(i, j): &(uint, uint)) -> T
+  { self.mij[NMat::offset::<D, T>(i, j)] }
+}
+
+impl<D: Dim, T:Copy + One + Zero> One for NMat<D, T>
+{
+  fn one() -> NMat<D, T>
+  {
+    let     dim = Dim::dim::<D>();
+    let mut res = NMat{ mij: from_elem(dim * dim, Zero::zero()) };
+    let     _1  = One::one::<T>();
+
+    for uint::range(0u, dim) |i|
+    { res.set(i, i, &_1); }
+
+    res
+  }
+}
+
+impl<D: Dim, T:Copy + Zero> Zero for NMat<D, T>
+{
+  fn zero() -> NMat<D, T>
+  {
+    let dim = Dim::dim::<D>();
+
+    NMat{ mij: from_elem(dim * dim, Zero::zero()) }
+  }
+}
+
+impl<D: Dim, T:Copy + Mul<T, T> + Add<T, T> + Zero>
+Mul<NMat<D, T>, NMat<D, T>> for NMat<D, T>
+{
+  fn mul(&self, other: &NMat<D, T>) -> NMat<D, T>
+  {
+    let     dim = Dim::dim::<D>();
+    let mut res = Zero::zero::<NMat<D, T>>();
+
+    for uint::range(0u, dim) |i|
+    {
+      for uint::range(0u, dim) |j|
+      {
+        let mut acc: T = Zero::zero();
+
+        for uint::range(0u, dim) |k|
+        { acc += self[(i, k)] * other[(k, j)]; }
+
+        res.set(i, j, &acc);
+      }
+    }
+
+    res
+  }
+}
+
+impl<D: Dim,
+     T: Copy + Mul<T, T> + Div<T, T> + Sub<T, T> + Neg<T> + Eq + One + Zero>
+Inv for NMat<D, T>
+{
+  fn inv(&self) -> NMat<D, T>
+  {
+    let mut cpy = copy *self;
+    let     dim = Dim::dim::<D>();
+    let mut res = One::one::<NMat<D, T>>();
+    let     _0T = Zero::zero::<T>();
+
+    // inversion using Gauss-Jordan elimination
+    for uint::range(0u, dim) |k|
+    {
+      // search a non-zero value on the k-th column
+      // FIXME: is it worth it to spend some more time searching for the max
+      // instead?
+
+      // FIXME: this is kind of uggly…
+      // … but we cannot use position_betwee since we are iterating on one
+      // columns
+      let mut n0 = dim; // index of a non-zero entry
+
+      for uint::range(k, dim) |i|
+      {
+        n0 = k;
+
+        if (cpy[(i, k)] != _0T)
+        { break; }
+      }
+
+      assert!(n0 != dim); // non inversible matrix
+
+      // swap pivot line
+      for uint::range(0u, dim) |j|
+      {
+        swap(cpy.mij, NMat::offset::<D, T>(n0, j), NMat::offset::<D, T>(k, j));
+        swap(res.mij, NMat::offset::<D, T>(n0, j), NMat::offset::<D, T>(k, j));
+      }
+
+      let pivot = cpy[(k, k)];
+
+      for uint::range(k, dim) |j|
+      {
+        cpy.set(k, j, &(cpy[(k, j)] / pivot));
+        res.set(k, j, &(res[(k, j)] / pivot));
+      }
+
+      for uint::range(0u, dim) |l|
+      {
+        if (l != k)
+        {
+          let normalizer = cpy[(l, k)] / pivot;
+
+          for uint::range(k, dim) |j|
+          {
+            cpy.set(k, j, &(cpy[(l, j)] - cpy[(k, j)] * normalizer));
+            res.set(k, j, &(res[(l, j)] - res[(k, j)] * normalizer));
+          }
+        }
+      }
+    }
+
+    res
+  }
+}
+
+impl<D: Dim, T:Copy> Transpose for NMat<D, T>
+{
+  fn transposed(&self) -> NMat<D, T>
+  {
+    let mut res = copy *self;
+
+    res.transpose();
+
+    res
+  }
+
+  fn transpose(&mut self)
+  {
+    let dim = Dim::dim::<D>();
+
+    for uint::range(1u, dim) |i|
+    {
+      for uint::range(0u, dim - 1) |j|
+      {
+        swap(self.mij,
+             NMat::offset::<D, T>(i, j),
+             NMat::offset::<D, T>(j, i));
+      }
+    }
+  }
+}
+
+impl<D: Dim, T:ToStr> ToStr for NMat<D, T>
+{
+  fn to_str(&self) -> ~str
+  { ~"Mat" + Dim::dim::<D>().to_str() + " {" + self.mij.to_str() + " }" }
+}

src/ndim/nvec.rs

@@ -0,0 +1,79 @@
+use core::vec::{map2, from_elem};
+use core::num::Zero;
+use traits::dim::Dim;
+use traits::dot::Dot;
+use traits::workarounds::sqrt::Sqrt;
+
+// D is a phantom parameter, used only as a dimensional token.
+// Its allows use to encode the vector dimension at the type-level.
+// It can be anything implementing the Dim trait. However, to avoid confusion,
+// using d0, d1, d2, d3 and d4 tokens are prefered.
+// FIXME: it might be possible to implement type-level integers and use them
+// here?
+#[deriving(Eq)]
+pub struct NVec<D, T>
+{
+  at: ~[T]
+}
+
+
+impl<D: Dim, T> Dim for NVec<D, T>
+{
+  fn dim() -> uint
+  { Dim::dim::<D>() }
+}
+
+impl<D, T:Copy + Add<T,T>> Add<NVec<D, T>, NVec<D, T>> for NVec<D, T>
+{
+  fn add(&self, other: &NVec<D, T>) -> NVec<D, T>
+  { NVec { at: map2(self.at, other.at, | a, b | { *a + *b }) } }
+}
+
+impl<D, T:Copy + Sub<T,T>> Sub<NVec<D, T>, NVec<D, T>> for NVec<D, T>
+{
+  fn sub(&self, other: &NVec<D, T>) -> NVec<D, T>
+  { NVec { at: map2(self.at, other.at, | a, b | *a - *b) } }
+}
+
+impl<D: Dim, T:Copy + Mul<T, T> + Add<T, T> + Sqrt + Zero> Dot<T> for NVec<D, T>
+{
+  fn dot(&self, other: &NVec<D, T>) -> T
+  {
+    let mut res = Zero::zero::<T>();
+
+    for uint::range(0u, Dim::dim::<D>()) |i|
+    { res += self.at[i] * other.at[i]; }
+
+    res
+  } 
+
+  fn sqnorm(&self) -> T
+  { self.dot(self) }
+
+  fn norm(&self) -> T
+  { self.sqnorm().sqrt() }
+}
+
+// FIXME: I dont really know how te generalize the cross product int
+// n-dimensions…
+// impl<T:Copy + Mul<T, T> + Sub<T, T>> Cross<T> for NVec<D, T>
+// {
+//   fn cross(&self, other: &NVec<D, T>) -> T
+//   { self.x * other.y - self.y * other.x }
+// }
+
+impl<D: Dim, T:Copy + Zero> Zero for NVec<D, T>
+{
+  fn zero() -> NVec<D, T>
+  {
+    let _0 = Zero::zero();
+
+    NVec { at: from_elem(Dim::dim::<D>(), _0) }
+  }
+}
+
+impl<D: Dim, T:ToStr> ToStr for NVec<D, T>
+{
+  fn to_str(&self) -> ~str
+  { ~"Vec" + Dim::dim::<D>().to_str() + self.at.to_str() }
+}

src/traits/dim.rs

@@ -2,3 +2,26 @@ pub trait Dim
 {
   fn dim() -> uint;
 }
+
+// Some dimension token. Useful to restrict the dimension of n-dimensional
+// object at the type-level.
+pub struct d0;
+pub struct d1;
+pub struct d2;
+pub struct d3;
+pub struct d4;
+
+impl Dim for d0
+{ fn dim() -> uint { 0 } }
+
+impl Dim for d1
+{ fn dim() -> uint { 1 } }
+
+impl Dim for d2
+{ fn dim() -> uint { 2 } }
+
+impl Dim for d3
+{ fn dim() -> uint { 3 } }
+
+impl Dim for d4
+{ fn dim() -> uint { 4 } }

src/traits/inv.rs

@@ -1,4 +1,4 @@
-pub trait Inv<T>
+pub trait Inv
 {
   fn inv(&self) -> Self;
 }

src/traits/sqrt.rs

@@ -1,17 +0,0 @@
-// FIXME: this does not seem to exist already
-// but it will surely be added someday.
-
-pub trait Sqrt
-{
-  fn sqrt(&self) -> Self;
-}
-
-impl Sqrt for f64
-{
-  fn sqrt(&self) -> f64 { f64::sqrt(*self) }
-}
-
-impl Sqrt for f32
-{
-  fn sqrt(&self) -> f32 { f32::sqrt(*self) }
-}