Add dynamically sized vector.

src/nalgebra.rc

@@ -39,6 +39,8 @@ pub mod dim3
 /// n-dimensional linear algebra (slower).
 pub mod ndim
 {
+  // pub mod dmat;
+  pub mod dvec;
   pub mod nvec;
   pub mod nmat;
 }

src/ndim/dvec.rs

@@ -0,0 +1,233 @@
+use core::num::{Zero, One, Algebraic};
+use core::vec::{map_zip, map, all2, len, from_elem, all};
+use core::cmp::ApproxEq;
+use traits::ring::Ring;
+use traits::division_ring::DivisionRing;
+use traits::dot::Dot;
+use traits::sub_dot::SubDot;
+use traits::norm::Norm;
+use traits::translation::Translation;
+use traits::workarounds::scalar_op::{ScalarMul, ScalarDiv, ScalarAdd, ScalarSub};
+
+#[deriving(Eq, ToStr, Clone)]
+pub struct DVec<T>
+{
+  at: ~[T]
+}
+
+pub fn zero_with_dim<T: Zero + Copy>(dim: uint) -> DVec<T>
+{ DVec { at: from_elem(dim, Zero::zero::<T>()) } }
+
+pub fn is_zero<T: Zero>(vec: &DVec<T>) -> bool
+{ all(vec.at, |e| e.is_zero()) }
+
+// FIXME: is Clone needed?
+impl<T: Copy + DivisionRing + Algebraic + Clone + ApproxEq<T>> DVec<T>
+{
+  pub fn canonical_basis_with_dim(dim: uint) -> ~[DVec<T>]
+  {
+    let mut res : ~[DVec<T>] = ~[];
+
+    for uint::range(0u, dim) |i|
+    {
+      let mut basis_element : DVec<T> = zero_with_dim(dim);
+
+      basis_element.at[i] = One::one();
+
+      res.push(basis_element);
+    }
+
+    res
+  }
+
+  pub fn orthogonal_subspace_basis(&self) -> ~[DVec<T>]
+  {
+    // compute the basis of the orthogonal subspace using Gram-Schmidt
+    // orthogonalization algorithm
+    let     dim              = len(self.at);
+    let mut res : ~[DVec<T>] = ~[];
+
+    for uint::range(0u, dim) |i|
+    {
+      let mut basis_element : DVec<T> = zero_with_dim(len(self.at));
+
+      basis_element.at[i] = One::one();
+
+      if (res.len() == dim - 1)
+      { break; }
+
+      let mut elt = basis_element.clone();
+
+      elt -= self.scalar_mul(&basis_element.dot(self));
+
+      for res.each |v|
+      { elt -= v.scalar_mul(&elt.dot(v)) };
+
+      if (!elt.sqnorm().approx_eq(&Zero::zero()))
+      { res.push(elt.normalized()); }
+    }
+
+    assert!(res.len() == dim - 1);
+
+    res
+  }
+}
+
+impl<T: Copy + Add<T,T>> Add<DVec<T>, DVec<T>> for DVec<T>
+{
+  fn add(&self, other: &DVec<T>) -> DVec<T>
+  {
+    assert!(len(self.at) == len(other.at));
+    DVec { at: map_zip(self.at, other.at, | a, b | { *a + *b }) }
+  }
+}
+
+impl<T: Copy + Sub<T,T>> Sub<DVec<T>, DVec<T>> for DVec<T>
+{
+  fn sub(&self, other: &DVec<T>) -> DVec<T>
+  {
+    assert!(len(self.at) == len(other.at));
+    DVec { at: map_zip(self.at, other.at, | a, b | *a - *b) }
+  }
+}
+
+impl<T: Copy + Neg<T>> Neg<DVec<T>> for DVec<T>
+{
+  fn neg(&self) -> DVec<T>
+  { DVec { at: map(self.at, |a| -a) } }
+}
+
+impl<T: Copy + Ring>
+Dot<T> for DVec<T>
+{
+  fn dot(&self, other: &DVec<T>) -> T
+  {
+    assert!(len(self.at) == len(other.at));
+
+    let mut res = Zero::zero::<T>();
+
+    for uint::range(0u, len(self.at)) |i|
+    { res += self.at[i] * other.at[i]; }
+
+    res
+  } 
+}
+
+impl<T: Copy + Ring> SubDot<T> for DVec<T>
+{
+  fn sub_dot(&self, a: &DVec<T>, b: &DVec<T>) -> T
+  {
+    let mut res = Zero::zero::<T>();
+
+    for uint::range(0u, len(self.at)) |i|
+    { res += (self.at[i] - a.at[i]) * b.at[i]; }
+
+    res
+  } 
+}
+
+impl<T: Copy + Mul<T, T>>
+ScalarMul<T> for DVec<T>
+{
+  fn scalar_mul(&self, s: &T) -> DVec<T>
+  { DVec { at: map(self.at, |a| a * *s) } }
+
+  fn scalar_mul_inplace(&mut self, s: &T)
+  {
+    for uint::range(0u, len(self.at)) |i|
+    { self.at[i] *= *s; }
+  }
+}
+
+
+impl<T: Copy + Div<T, T>>
+ScalarDiv<T> for DVec<T>
+{
+  fn scalar_div(&self, s: &T) -> DVec<T>
+  { DVec { at: map(self.at, |a| a / *s) } }
+
+  fn scalar_div_inplace(&mut self, s: &T)
+  {
+    for uint::range(0u, len(self.at)) |i|
+    { self.at[i] /= *s; }
+  }
+}
+
+impl<T: Copy + Add<T, T>>
+ScalarAdd<T> for DVec<T>
+{
+  fn scalar_add(&self, s: &T) -> DVec<T>
+  { DVec { at: map(self.at, |a| a + *s) } }
+
+  fn scalar_add_inplace(&mut self, s: &T)
+  {
+    for uint::range(0u, len(self.at)) |i|
+    { self.at[i] += *s; }
+  }
+}
+
+impl<T: Copy + Sub<T, T>>
+ScalarSub<T> for DVec<T>
+{
+  fn scalar_sub(&self, s: &T) -> DVec<T>
+  { DVec { at: map(self.at, |a| a - *s) } }
+
+  fn scalar_sub_inplace(&mut self, s: &T)
+  {
+    for uint::range(0u, len(self.at)) |i|
+    { self.at[i] -= *s; }
+  }
+}
+
+impl<T: Clone + Copy + Add<T, T>> Translation<DVec<T>> for DVec<T>
+{
+  fn translation(&self) -> DVec<T>
+  { self.clone() }
+
+  fn translated(&self, t: &DVec<T>) -> DVec<T>
+  { self + *t }
+
+  fn translate(&mut self, t: &DVec<T>)
+  { *self = *self + *t; }
+}
+
+impl<T: Copy + DivisionRing + Algebraic + Clone>
+Norm<T> for DVec<T>
+{
+  fn sqnorm(&self) -> T
+  { self.dot(self) }
+
+  fn norm(&self) -> T
+  { self.sqnorm().sqrt() }
+
+  fn normalized(&self) -> DVec<T>
+  {
+    let mut res : DVec<T> = self.clone();
+
+    res.normalize();
+
+    res
+  }
+
+  fn normalize(&mut self) -> T
+  {
+    let l = self.norm();
+
+    for uint::range(0u, len(self.at)) |i|
+    { self.at[i] /= l; }
+
+    l
+  }
+}
+
+impl<T: ApproxEq<T>> ApproxEq<T> for DVec<T>
+{
+  fn approx_epsilon() -> T
+  { ApproxEq::approx_epsilon::<T, T>() }
+
+  fn approx_eq(&self, other: &DVec<T>) -> bool
+  { all2(self.at, other.at, |a, b| a.approx_eq(b)) }
+
+  fn approx_eq_eps(&self, other: &DVec<T>, epsilon: &T) -> bool
+  { all2(self.at, other.at, |a, b| a.approx_eq_eps(b, epsilon)) }
+}

src/ndim/nmat.rs

@@ -109,7 +109,7 @@ RMul<NVec<D, T>> for NMat<D, T>
     for uint::range(0u, dim) |i|
     {
       for uint::range(0u, dim) |j|
-      { res.at[i] = res.at[i] + other.at[j] * self[(i, j)]; }
+      { res.at.at[i] = res.at.at[i] + other.at.at[j] * self[(i, j)]; }
     }
 
     res
@@ -127,7 +127,7 @@ LMul<NVec<D, T>> for NMat<D, T>
     for uint::range(0u, dim) |i|
     {
       for uint::range(0u, dim) |j|
-      { res.at[i] = res.at[i] + other.at[j] * self[(j, i)]; }
+      { res.at.at[i] = res.at.at[i] + other.at.at[j] * self[(j, i)]; }
     }
 
     res
@@ -158,8 +158,8 @@ Inv for NMat<D, T>
     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
+      // FIXME: would it be worth it to spend some more time searching for the
-      // instead?
+      // max instead?
 
       // FIXME: this is kind of uggly…
       // … but we cannot use position_between since we are iterating on one

src/ndim/nvec.rs

@@ -1,7 +1,8 @@
-use core::num::{Zero, One, Algebraic};
+use core::num::{Zero, Algebraic};
 use core::rand::{Rand, Rng, RngUtil};
-use core::vec::{map_zip, from_elem, map, all, all2};
+use core::vec::{map};
 use core::cmp::ApproxEq;
+use ndim::dvec::{DVec, zero_with_dim, is_zero};
 use traits::basis::Basis;
 use traits::ring::Ring;
 use traits::division_ring::DivisionRing;
@@ -15,14 +16,12 @@ use traits::workarounds::scalar_op::{ScalarMul, ScalarDiv, ScalarAdd, ScalarSub}
 // D is a phantom parameter, used only as a dimensional token.
 // Its allows use to encode the vector dimension at the type-level.
 // It can be anything implementing the Dim trait. However, to avoid confusion,
-// using d0, d1, d2, d3 and d4 tokens are prefered.
+// using d0, d1, d2, d3, ..., d7 (or your own dn) are prefered.
 // FIXME: it might be possible to implement type-level integers and use them
 // here?
 #[deriving(Eq, ToStr)]
 pub struct NVec<D, T>
-{
+{ at: DVec<T> }
-  at: ~[T]
-}
 
 
 impl<D: Dim, T> Dim for NVec<D, T>
@@ -40,59 +39,42 @@ impl<D, T: Clone> Clone for NVec<D, T>
 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: map_zip(self.at, other.at, | a, b | { *a + *b }) } }
+  { NVec { at: self.at + other.at } }
 }
 
 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: map_zip(self.at, other.at, | a, b | *a - *b) } }
+  { NVec { at: self.at - other.at } }
 }
 
 impl<D, T: Copy + Neg<T>> Neg<NVec<D, T>> for NVec<D, T>
 {
   fn neg(&self) -> NVec<D, T>
-  { NVec { at: map(self.at, |a| -a) } }
+  { NVec { at: -self.at } }
 }
 
 impl<D: Dim, T: Copy + Ring>
 Dot<T> for NVec<D, T>
 {
   fn dot(&self, other: &NVec<D, T>) -> T
-  {
+  { self.at.dot(&other.at) } 
-    let mut res = Zero::zero::<T>();
-
-    for uint::range(0u, Dim::dim::<D>()) |i|
-    { res += self.at[i] * other.at[i]; }
-
-    res
-  } 
 }
 
 impl<D: Dim, T: Copy + Ring> SubDot<T> for NVec<D, T>
 {
   fn sub_dot(&self, a: &NVec<D, T>, b: &NVec<D, T>) -> T
-  {
+  { self.at.sub_dot(&a.at, &b.at) } 
-    let mut res = Zero::zero::<T>();
-
-    for uint::range(0u, Dim::dim::<D>()) |i|
-    { res += (self.at[i] - a.at[i]) * b.at[i]; }
-
-    res
-  } 
 }
 
 impl<D: Dim, T: Copy + Mul<T, T>>
 ScalarMul<T> for NVec<D, T>
 {
   fn scalar_mul(&self, s: &T) -> NVec<D, T>
-  { NVec { at: map(self.at, |a| a * *s) } }
+  { NVec { at: self.at.scalar_mul(s) } }
 
   fn scalar_mul_inplace(&mut self, s: &T)
-  {
+  { self.at.scalar_mul_inplace(s) }
-    for uint::range(0u, Dim::dim::<D>()) |i|
-    { self.at[i] *= *s; }
-  }
 }
 
 
@@ -100,39 +82,30 @@ impl<D: Dim, T: Copy + Div<T, T>>
 ScalarDiv<T> for NVec<D, T>
 {
   fn scalar_div(&self, s: &T) -> NVec<D, T>
-  { NVec { at: map(self.at, |a| a / *s) } }
+  { NVec { at: self.at.scalar_div(s) } }
 
   fn scalar_div_inplace(&mut self, s: &T)
-  {
+  { self.at.scalar_div_inplace(s) }
-    for uint::range(0u, Dim::dim::<D>()) |i|
-    { self.at[i] /= *s; }
-  }
 }
 
 impl<D: Dim, T: Copy + Add<T, T>>
 ScalarAdd<T> for NVec<D, T>
 {
   fn scalar_add(&self, s: &T) -> NVec<D, T>
-  { NVec { at: map(self.at, |a| a + *s) } }
+  { NVec { at: self.at.scalar_add(s) } }
 
   fn scalar_add_inplace(&mut self, s: &T)
-  {
+  { self.at.scalar_add_inplace(s) }
-    for uint::range(0u, Dim::dim::<D>()) |i|
-    { self.at[i] += *s; }
-  }
 }
 
 impl<D: Dim, T: Copy + Sub<T, T>>
 ScalarSub<T> for NVec<D, T>
 {
   fn scalar_sub(&self, s: &T) -> NVec<D, T>
-  { NVec { at: map(self.at, |a| a - *s) } }
+  { NVec { at: self.at.scalar_sub(s) } }
 
   fn scalar_sub_inplace(&mut self, s: &T)
-  {
+  { self.scalar_sub_inplace(s) }
-    for uint::range(0u, Dim::dim::<D>()) |i|
-    { self.at[i] -= *s; }
-  }
 }
 
 impl<D: Dim, T: Clone + Copy + Add<T, T>> Translation<NVec<D, T>> for NVec<D, T>
@@ -166,14 +139,7 @@ Norm<T> for NVec<D, T>
   }
 
   fn normalize(&mut self) -> T
-  {
+  { self.at.normalize() }
-    let l = self.norm();
-
-    for uint::range(0u, Dim::dim::<D>()) |i|
-    { self.at[i] /= l; }
-
-    l
-  }
 }
 
 impl<D: Dim,
@@ -181,53 +147,10 @@ impl<D: Dim,
 Basis for NVec<D, T>
 {
   fn canonical_basis() -> ~[NVec<D, T>]
-  {
+  { map(DVec::canonical_basis_with_dim(Dim::dim::<D>()), |&e| NVec { at: e }) }
-    let     dim = Dim::dim::<D>();
-    let mut res : ~[NVec<D, T>] = ~[];
-
-    for uint::range(0u, dim) |i|
-    {
-      let mut basis_element : NVec<D, T> = Zero::zero();
-
-      basis_element.at[i] = One::one();
-
-      res.push(basis_element);
-    }
-
-    res
-  }
 
   fn orthogonal_subspace_basis(&self) -> ~[NVec<D, T>]
-  {
+  { map(self.at.orthogonal_subspace_basis(), |&e| NVec { at: e }) }
-    // compute the basis of the orthogonal subspace using Gram-Schmidt
-    // orthogonalization algorithm
-    let     dim = Dim::dim::<D>();
-    let mut res : ~[NVec<D, T>] = ~[];
-
-    for uint::range(0u, dim) |i|
-    {
-      let mut basis_element : NVec<D, T> = Zero::zero();
-
-      basis_element.at[i] = One::one();
-
-      if (res.len() == dim - 1)
-      { break; }
-
-      let mut elt = basis_element.clone();
-
-      elt -= self.scalar_mul(&basis_element.dot(self));
-
-      for res.each |v|
-      { elt -= v.scalar_mul(&elt.dot(v)) };
-
-      if (!elt.sqnorm().approx_eq(&Zero::zero()))
-      { res.push(elt.normalized()); }
-    }
-
-    assert!(res.len() == dim - 1);
-
-    res
-  }
 }
 
 // FIXME: I dont really know how te generalize the cross product int
@@ -241,16 +164,10 @@ Basis for NVec<D, T>
 impl<D: Dim, T: Copy + Zero> Zero for NVec<D, T>
 {
   fn zero() -> NVec<D, T>
-  {
+  { NVec { at: zero_with_dim(Dim::dim::<D>()) } }
-    let _0 = Zero::zero();
-
-    NVec { at: from_elem(Dim::dim::<D>(), _0) }
-  }
 
   fn is_zero(&self) -> bool
-  {
+  { is_zero(&self.at) }
-    all(self.at, |e| e.is_zero())
-  }
 }
 
 impl<D, T: ApproxEq<T>> ApproxEq<T> for NVec<D, T>
@@ -259,10 +176,10 @@ impl<D, T: ApproxEq<T>> ApproxEq<T> for NVec<D, T>
   { ApproxEq::approx_epsilon::<T, T>() }
 
   fn approx_eq(&self, other: &NVec<D, T>) -> bool
-  { all2(self.at, other.at, |a, b| a.approx_eq(b)) }
+  { self.at.approx_eq(&other.at) }
 
   fn approx_eq_eps(&self, other: &NVec<D, T>, epsilon: &T) -> bool
-  { all2(self.at, other.at, |a, b| a.approx_eq_eps(b, epsilon)) }
+  { self.at.approx_eq_eps(&other.at, epsilon) }
 }
 
 impl<D: Dim, T: Rand + Zero + Copy> Rand for NVec<D, T>
@@ -273,7 +190,7 @@ impl<D: Dim, T: Rand + Zero + Copy> Rand for NVec<D, T>
     let mut res : NVec<D, T> = Zero::zero();
 
     for uint::range(0u, dim) |i|
-    { res.at[i] = rng.gen() }
+    { res.at.at[i] = rng.gen() }
 
     res
   }