From c54eb562ecc5fda22ac83180c59cfddef8e036f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Fri, 28 Jun 2013 22:55:09 +0000 Subject: [PATCH] Refactor code for matrices. --- src/adaptors/rotmat.rs | 29 +- src/dim1/mat1.rs | 137 ---------- src/dim2/mat2.rs | 190 ------------- src/dim3/mat3.rs | 244 ----------------- src/mat.rs | 600 +++++++++++++++++++++++++++++++++++++++++ src/nalgebra.rc | 19 +- src/tests/mat.rs | 6 +- 7 files changed, 615 insertions(+), 610 deletions(-) delete mode 100644 src/dim1/mat1.rs delete mode 100644 src/dim2/mat2.rs delete mode 100644 src/dim3/mat3.rs create mode 100644 src/mat.rs diff --git a/src/adaptors/rotmat.rs b/src/adaptors/rotmat.rs index 59354e9f..87a0ff26 100644 --- a/src/adaptors/rotmat.rs +++ b/src/adaptors/rotmat.rs @@ -1,6 +1,8 @@ use std::num::{One, Zero}; use std::rand::{Rand, Rng, RngUtil}; use std::cmp::ApproxEq; +use traits::ring::Ring; +use traits::division_ring::DivisionRing; use traits::rlmul::{RMul, LMul}; use traits::dim::Dim; use traits::inv::Inv; @@ -8,8 +10,7 @@ use traits::transpose::Transpose; use traits::rotation::{Rotation, Rotate, Rotatable}; use traits::transformation::{Transform}; // FIXME: implement Transformation and Transformable use vec::Vec1; -use dim2::mat2::Mat2; -use dim3::mat3::Mat3; +use mat::{Mat2, Mat3}; use vec::Vec3; #[deriving(Eq, ToStr)] @@ -22,11 +23,10 @@ pub fn rotmat2>(angle: N) -> Rotmat> let sia = angle.sin(); Rotmat - { submat: Mat2::new(copy coa, -sia, copy sia, copy coa) } + { submat: Mat2::new( [ copy coa, -sia, copy sia, copy coa ] ) } } -pub fn rotmat3 + One + Sub + Add + - Mul> +pub fn rotmat3 (axis: &Vec3, angle: N) -> Rotmat> { let _1 = One::one::(); @@ -41,7 +41,7 @@ pub fn rotmat3 + One + Sub + Add + let sin = angle.sin(); Rotmat { - submat: Mat3::new( + submat: Mat3::new( [ (sqx + (_1 - sqx) * cos), (ux * uy * one_m_cos - uz * sin), (ux * uz * one_m_cos + uy * sin), @@ -52,16 +52,16 @@ pub fn rotmat3 + One + Sub + Add + (ux * uz * one_m_cos - uy * sin), (uy * uz * one_m_cos + ux * sin), - (sqz + (_1 - sqz) * cos)) + (sqz + (_1 - sqz) * cos) ] ) } } -impl + Trigonometric + Neg + Mul + Add + Copy> +impl Rotation> for Rotmat> { #[inline] fn rotation(&self) -> Vec1 - { Vec1::new([-(self.submat.m12 / self.submat.m11).atan()]) } + { Vec1::new([ -(self.submat.at(0, 1) / self.submat.at(0, 0)).atan() ]) } #[inline] fn inv_rotation(&self) -> Vec1 @@ -72,7 +72,7 @@ Rotation> for Rotmat> { *self = self.rotated(rot) } } -impl + Trigonometric + Neg + Mul + Add + Copy> +impl Rotatable, Rotmat>> for Rotmat> { #[inline] @@ -80,8 +80,7 @@ Rotatable, Rotmat>> for Rotmat> { rotmat2(copy rot.at[0]) * *self } } -impl + Trigonometric + Neg + Mul + Add + Copy + - One + Sub> +impl Rotation<(Vec3, N)> for Rotmat> { #[inline] @@ -98,8 +97,7 @@ Rotation<(Vec3, N)> for Rotmat> { *self = self.rotated(rot) } } -impl + Trigonometric + Neg + Mul + Add + Copy + - One + Sub> +impl Rotatable<(Vec3, N), Rotmat>> for Rotmat> { #[inline] @@ -136,8 +134,7 @@ impl + LMul, V> Transform for Rotmat { self.inv_rotate(v) } } -impl + One + Sub + Add + - Mul> +impl Rand for Rotmat> { #[inline] diff --git a/src/dim1/mat1.rs b/src/dim1/mat1.rs deleted file mode 100644 index 540a3947..00000000 --- a/src/dim1/mat1.rs +++ /dev/null @@ -1,137 +0,0 @@ -use std::num::{One, Zero}; -use std::rand::{Rand, Rng, RngUtil}; -use std::cmp::ApproxEq; -use traits::division_ring::DivisionRing; -use traits::dim::Dim; -use traits::inv::Inv; -use traits::transpose::Transpose; -use traits::transformation::Transform; // FIXME: implement Transformable, Transformation -use traits::rlmul::{RMul, LMul}; -use vec::Vec1; - -#[deriving(Eq, ToStr)] -pub struct Mat1 -{ m11: N } - -impl Mat1 -{ - #[inline] - pub fn new(m11: N) -> Mat1 - { - Mat1 - { m11: m11 } - } -} - -impl Dim for Mat1 -{ - #[inline] - fn dim() -> uint - { 1 } -} - -impl One for Mat1 -{ - #[inline] - fn one() -> Mat1 - { return Mat1::new(One::one()) } -} - -impl Zero for Mat1 -{ - #[inline] - fn zero() -> Mat1 - { Mat1::new(Zero::zero()) } - - #[inline] - fn is_zero(&self) -> bool - { self.m11.is_zero() } -} - -impl + Add> Mul, Mat1> for Mat1 -{ - #[inline] - fn mul(&self, other: &Mat1) -> Mat1 - { Mat1::new(self.m11 * other.m11) } -} - -impl -Transform> for Mat1 -{ - #[inline] - fn transform_vec(&self, v: &Vec1) -> Vec1 - { self.rmul(v) } - - #[inline] - fn inv_transform(&self, v: &Vec1) -> Vec1 - { self.inverse().transform_vec(v) } -} - -impl + Mul> RMul> for Mat1 -{ - #[inline] - fn rmul(&self, other: &Vec1) -> Vec1 - { Vec1::new([self.m11 * other.at[0]]) } -} - -impl + Mul> LMul> for Mat1 -{ - #[inline] - fn lmul(&self, other: &Vec1) -> Vec1 - { Vec1::new([self.m11 * other.at[0]]) } -} - -impl -Inv for Mat1 -{ - #[inline] - fn inverse(&self) -> Mat1 - { - let mut res : Mat1 = copy *self; - - res.invert(); - - res - } - - #[inline] - fn invert(&mut self) - { - assert!(!self.m11.is_zero()); - - self.m11 = One::one::() / self.m11 - } -} - -impl Transpose for Mat1 -{ - #[inline] - fn transposed(&self) -> Mat1 - { copy *self } - - #[inline] - fn transpose(&mut self) - { } -} - -impl> ApproxEq for Mat1 -{ - #[inline] - fn approx_epsilon() -> N - { ApproxEq::approx_epsilon::() } - - #[inline] - fn approx_eq(&self, other: &Mat1) -> bool - { self.m11.approx_eq(&other.m11) } - - #[inline] - fn approx_eq_eps(&self, other: &Mat1, epsilon: &N) -> bool - { self.m11.approx_eq_eps(&other.m11, epsilon) } -} - -impl Rand for Mat1 -{ - #[inline] - fn rand(rng: &mut R) -> Mat1 - { Mat1::new(rng.gen()) } -} diff --git a/src/dim2/mat2.rs b/src/dim2/mat2.rs deleted file mode 100644 index e24eca28..00000000 --- a/src/dim2/mat2.rs +++ /dev/null @@ -1,190 +0,0 @@ -use std::num::{One, Zero}; -use std::rand::{Rand, Rng, RngUtil}; -use std::cmp::ApproxEq; -use std::util::swap; -use traits::transformation::Transform; -use traits::division_ring::DivisionRing; -use traits::dim::Dim; -use traits::inv::Inv; -use traits::transpose::Transpose; -use traits::rlmul::{RMul, LMul}; -use vec::Vec2; - -#[deriving(Eq, ToStr)] -pub struct Mat2 -{ - m11: N, m12: N, - m21: N, m22: N -} - -impl Mat2 -{ - #[inline] - pub fn new(m11: N, m12: N, m21: N, m22: N) -> Mat2 - { - Mat2 - { - m11: m11, m12: m12, - m21: m21, m22: m22, - } - } -} - -impl Dim for Mat2 -{ - #[inline] - fn dim() -> uint - { 2 } -} - -impl One for Mat2 -{ - #[inline] - fn one() -> Mat2 - { - let (_0, _1) = (Zero::zero(), One::one()); - return Mat2::new(copy _1, copy _0, - _0, _1) - } -} - -impl Zero for Mat2 -{ - #[inline] - fn zero() -> Mat2 - { - let _0 = Zero::zero(); - return Mat2::new(copy _0, copy _0, - copy _0, _0) - } - - #[inline] - fn is_zero(&self) -> bool - { - self.m11.is_zero() && self.m12.is_zero() && - self.m21.is_zero() && self.m22.is_zero() - } -} - -impl + Add> Mul, Mat2> for Mat2 -{ - #[inline] - fn mul(&self, other: &Mat2) -> Mat2 - { - Mat2::new( - self.m11 * other.m11 + self.m12 * other.m21, - self.m11 * other.m12 + self.m12 * other.m22, - self.m21 * other.m11 + self.m22 * other.m21, - self.m21 * other.m12 + self.m22 * other.m22 - ) - } -} - -impl -Transform> for Mat2 -{ - #[inline] - fn transform_vec(&self, v: &Vec2) -> Vec2 - { self.rmul(v) } - - #[inline] - fn inv_transform(&self, v: &Vec2) -> Vec2 - { self.inverse().transform_vec(v) } -} - -impl + Mul> RMul> for Mat2 -{ - #[inline] - fn rmul(&self, other: &Vec2) -> Vec2 - { - Vec2::new( - [ self.m11 * other.at[0] + self.m12 * other.at[1], - self.m21 * other.at[0] + self.m22 * other.at[1] ] - ) - } -} - -impl + Mul> LMul> for Mat2 -{ - #[inline] - fn lmul(&self, other: &Vec2) -> Vec2 - { - Vec2::new( - [ self.m11 * other.at[0] + self.m21 * other.at[1], - self.m12 * other.at[0] + self.m22 * other.at[1] ] - ) - } -} - -impl -Inv for Mat2 -{ - #[inline] - fn inverse(&self) -> Mat2 - { - let mut res : Mat2 = copy *self; - - res.invert(); - - res - } - - #[inline] - fn invert(&mut self) - { - let det = self.m11 * self.m22 - self.m21 * self.m12; - - assert!(!det.is_zero()); - - *self = Mat2::new(self.m22 / det , -self.m12 / det, - -self.m21 / det, self.m11 / det) - } -} - -impl Transpose for Mat2 -{ - #[inline] - fn transposed(&self) -> Mat2 - { - Mat2::new(copy self.m11, copy self.m21, - copy self.m12, copy self.m22) - } - - #[inline] - fn transpose(&mut self) - { swap(&mut self.m21, &mut self.m12); } -} - -impl> ApproxEq for Mat2 -{ - #[inline] - fn approx_epsilon() -> N - { ApproxEq::approx_epsilon::() } - - #[inline] - fn approx_eq(&self, other: &Mat2) -> bool - { - self.m11.approx_eq(&other.m11) && - self.m12.approx_eq(&other.m12) && - - self.m21.approx_eq(&other.m21) && - self.m22.approx_eq(&other.m22) - } - - #[inline] - fn approx_eq_eps(&self, other: &Mat2, epsilon: &N) -> bool - { - self.m11.approx_eq_eps(&other.m11, epsilon) && - self.m12.approx_eq_eps(&other.m12, epsilon) && - - self.m21.approx_eq_eps(&other.m21, epsilon) && - self.m22.approx_eq_eps(&other.m22, epsilon) - } -} - -impl Rand for Mat2 -{ - #[inline] - fn rand(rng: &mut R) -> Mat2 - { Mat2::new(rng.gen(), rng.gen(), rng.gen(), rng.gen()) } -} diff --git a/src/dim3/mat3.rs b/src/dim3/mat3.rs deleted file mode 100644 index 098d0f57..00000000 --- a/src/dim3/mat3.rs +++ /dev/null @@ -1,244 +0,0 @@ -use std::num::{One, Zero}; -use std::rand::{Rand, Rng, RngUtil}; -use std::cmp::ApproxEq; -use std::util::swap; -use traits::dim::Dim; -use traits::inv::Inv; -use traits::transformation::Transform; -use traits::division_ring::DivisionRing; -use traits::transpose::Transpose; -use traits::rlmul::{RMul, LMul}; -use vec::Vec3; - -#[deriving(Eq, ToStr)] -pub struct Mat3 -{ - m11: N, m12: N, m13: N, - m21: N, m22: N, m23: N, - m31: N, m32: N, m33: N -} - -impl Mat3 -{ - #[inline] - pub fn new(m11: N, m12: N, m13: N, - m21: N, m22: N, m23: N, - m31: N, m32: N, m33: N) -> Mat3 - { - Mat3 - { - m11: m11, m12: m12, m13: m13, - m21: m21, m22: m22, m23: m23, - m31: m31, m32: m32, m33: m33 - } - } -} - -impl Dim for Mat3 -{ - #[inline] - fn dim() -> uint - { 3 } -} - -impl One for Mat3 -{ - #[inline] - fn one() -> Mat3 - { - let (_0, _1) = (Zero::zero(), One::one()); - return Mat3::new(copy _1, copy _0, copy _0, - copy _0, copy _1, copy _0, - copy _0, _0, _1) - } -} - -impl Zero for Mat3 -{ - #[inline] - fn zero() -> Mat3 - { - let _0 = Zero::zero(); - return Mat3::new(copy _0, copy _0, copy _0, - copy _0, copy _0, copy _0, - copy _0, copy _0, _0) - } - - #[inline] - fn is_zero(&self) -> bool - { - self.m11.is_zero() && self.m12.is_zero() && self.m13.is_zero() && - self.m21.is_zero() && self.m22.is_zero() && self.m23.is_zero() && - self.m31.is_zero() && self.m32.is_zero() && self.m33.is_zero() - } -} - -impl + Add> Mul, Mat3> for Mat3 -{ - #[inline] - fn mul(&self, other: &Mat3) -> Mat3 - { - Mat3::new( - self.m11 * other.m11 + self.m12 * other.m21 + self.m13 * other.m31, - self.m11 * other.m12 + self.m12 * other.m22 + self.m13 * other.m32, - self.m11 * other.m13 + self.m12 * other.m23 + self.m13 * other.m33, - - self.m21 * other.m11 + self.m22 * other.m21 + self.m23 * other.m31, - self.m21 * other.m12 + self.m22 * other.m22 + self.m23 * other.m32, - self.m21 * other.m13 + self.m22 * other.m23 + self.m23 * other.m33, - - self.m31 * other.m11 + self.m32 * other.m21 + self.m33 * other.m31, - self.m31 * other.m12 + self.m32 * other.m22 + self.m33 * other.m32, - self.m31 * other.m13 + self.m32 * other.m23 + self.m33 * other.m33 - ) - } -} - -impl -Transform> for Mat3 -{ - #[inline] - fn transform_vec(&self, v: &Vec3) -> Vec3 - { self.rmul(v) } - - #[inline] - fn inv_transform(&self, v: &Vec3) -> Vec3 - { self.inverse().transform_vec(v) } -} - -impl + Mul> RMul> for Mat3 -{ - #[inline] - fn rmul(&self, other: &Vec3) -> Vec3 - { - Vec3::new( - [self.m11 * other.at[0] + self.m12 * other.at[1] + self.m13 * other.at[2], - self.m21 * other.at[0] + self.m22 * other.at[1] + self.m33 * other.at[2], - self.m31 * other.at[0] + self.m32 * other.at[1] + self.m33 * other.at[2]] - ) - } -} - -impl + Mul> LMul> for Mat3 -{ - #[inline] - fn lmul(&self, other: &Vec3) -> Vec3 - { - Vec3::new( - [self.m11 * other.at[0] + self.m21 * other.at[1] + self.m31 * other.at[2], - self.m12 * other.at[0] + self.m22 * other.at[1] + self.m32 * other.at[2], - self.m13 * other.at[0] + self.m23 * other.at[1] + self.m33 * other.at[2]] - ) - } -} - -impl -Inv for Mat3 -{ - #[inline] - fn inverse(&self) -> Mat3 - { - let mut res = copy *self; - - res.invert(); - - res - } - - #[inline] - fn invert(&mut self) - { - let minor_m22_m33 = self.m22 * self.m33 - self.m32 * self.m23; - let minor_m21_m33 = self.m21 * self.m33 - self.m31 * self.m23; - let minor_m21_m32 = self.m21 * self.m32 - self.m31 * self.m22; - - let det = self.m11 * minor_m22_m33 - - self.m12 * minor_m21_m33 - + self.m13 * minor_m21_m32; - - assert!(!det.is_zero()); - - *self = Mat3::new( - (minor_m22_m33 / det), - ((self.m13 * self.m32 - self.m33 * self.m12) / det), - ((self.m12 * self.m23 - self.m22 * self.m13) / det), - - (-minor_m21_m33 / det), - ((self.m11 * self.m33 - self.m31 * self.m13) / det), - ((self.m13 * self.m21 - self.m23 * self.m11) / det), - - (minor_m21_m32 / det), - ((self.m12 * self.m31 - self.m32 * self.m11) / det), - ((self.m11 * self.m22 - self.m21 * self.m12) / det) - ) - } -} - -impl Transpose for Mat3 -{ - #[inline] - fn transposed(&self) -> Mat3 - { - Mat3::new(copy self.m11, copy self.m21, copy self.m31, - copy self.m12, copy self.m22, copy self.m32, - copy self.m13, copy self.m23, copy self.m33) - } - - #[inline] - fn transpose(&mut self) - { - swap(&mut self.m12, &mut self.m21); - swap(&mut self.m13, &mut self.m31); - swap(&mut self.m23, &mut self.m32); - } -} - -impl> ApproxEq for Mat3 -{ - #[inline] - fn approx_epsilon() -> N - { ApproxEq::approx_epsilon::() } - - #[inline] - fn approx_eq(&self, other: &Mat3) -> bool - { - self.m11.approx_eq(&other.m11) && - self.m12.approx_eq(&other.m12) && - self.m13.approx_eq(&other.m13) && - - self.m21.approx_eq(&other.m21) && - self.m22.approx_eq(&other.m22) && - self.m23.approx_eq(&other.m23) && - - self.m31.approx_eq(&other.m31) && - self.m32.approx_eq(&other.m32) && - self.m33.approx_eq(&other.m33) - } - - #[inline] - fn approx_eq_eps(&self, other: &Mat3, epsilon: &N) -> bool - { - self.m11.approx_eq_eps(&other.m11, epsilon) && - self.m12.approx_eq_eps(&other.m12, epsilon) && - self.m13.approx_eq_eps(&other.m13, epsilon) && - - self.m21.approx_eq_eps(&other.m21, epsilon) && - self.m22.approx_eq_eps(&other.m22, epsilon) && - self.m23.approx_eq_eps(&other.m23, epsilon) && - - self.m31.approx_eq_eps(&other.m31, epsilon) && - self.m32.approx_eq_eps(&other.m32, epsilon) && - self.m33.approx_eq_eps(&other.m33, epsilon) - } -} - -impl Rand for Mat3 -{ - #[inline] - fn rand(rng: &mut R) -> Mat3 - { - Mat3::new(rng.gen(), rng.gen(), rng.gen(), - rng.gen(), rng.gen(), rng.gen(), - rng.gen(), rng.gen(), rng.gen()) - } -} diff --git a/src/mat.rs b/src/mat.rs new file mode 100644 index 00000000..d1e0995d --- /dev/null +++ b/src/mat.rs @@ -0,0 +1,600 @@ +use std::uint::iterate; +use std::num::{One, Zero}; +use std::vec::swap; +use std::cmp::ApproxEq; +use std::rand::{Rand, Rng, RngUtil}; +use std::iterator::IteratorUtil; +use vec::{Vec1, Vec2, Vec3, Vec4, Vec5, Vec6}; +use traits::dim::Dim; +use traits::ring::Ring; +use traits::inv::Inv; +use traits::division_ring::DivisionRing; +use traits::transpose::Transpose; +use traits::rlmul::{RMul, LMul}; +use traits::transformation::Transform; + +macro_rules! mat_impl( + ($t: ident, $dim: expr) => ( + impl $t + { + #[inline] + pub fn new(mij: [N, ..$dim * $dim]) -> $t + { $t { mij: mij } } + } + ) +) + +macro_rules! one_impl( + ($t: ident, [ $($value: ident)|+ ] ) => ( + impl One for $t + { + #[inline] + fn one() -> $t + { + let (_0, _1) = (Zero::zero::(), One::one::()); + return $t::new( [ $( copy $value, )+ ] ) + } + } + ) +) + + +macro_rules! zero_impl( + ($t: ident, [ $($value: ident)|+ ] ) => ( + impl Zero for $t + { + #[inline] + fn zero() -> $t + { + let _0 = Zero::zero(); + return $t::new( [ $( copy $value, )+ ] ) + } + + #[inline] + fn is_zero(&self) -> bool + { self.mij.iter().all(|e| e.is_zero()) } + } + ) +) + +macro_rules! dim_impl( + ($t: ident, $dim: expr) => ( + impl Dim for $t + { + #[inline] + fn dim() -> uint + { $dim } + } + ) +) + +macro_rules! mat_indexing_impl( + ($t: ident, $dim: expr) => ( + impl $t + { + #[inline] + pub fn offset(&self, i: uint, j: uint) -> uint + { i * $dim + j } + + #[inline] + pub fn set(&mut self, i: uint, j: uint, t: &N) + { + self.mij[self.offset(i, j)] = copy *t + } + + #[inline] + pub fn at(&self, i: uint, j: uint) -> N + { + copy self.mij[self.offset(i, j)] + } + } + ) +) + +macro_rules! mul_impl( + ($t: ident, $dim: expr) => ( + impl + Mul<$t, $t> for $t + { + fn mul(&self, other: &$t) -> $t + { + let mut res: $t = Zero::zero(); + + for iterate(0u, $dim) |i| + { + for iterate(0u, $dim) |j| + { + let mut acc = Zero::zero::(); + + for iterate(0u, $dim) |k| + { acc = acc + self.at(i, k) * other.at(k, j); } + + res.set(i, j, &acc); + } + } + + res + } + } + ) +) + +macro_rules! rmul_impl( + ($t: ident, $v: ident, $dim: expr) => ( + impl + RMul<$v> for $t + { + fn rmul(&self, other: &$v) -> $v + { + let mut res : $v = Zero::zero(); + + for iterate(0u, $dim) |i| + { + for iterate(0u, $dim) |j| + { res.at[i] = res.at[i] + other.at[j] * self.at(i, j); } + } + + res + } + } + ) +) + +macro_rules! lmul_impl( + ($t: ident, $v: ident, $dim: expr) => ( + impl + LMul<$v> for $t + { + fn lmul(&self, other: &$v) -> $v + { + + let mut res : $v = Zero::zero(); + + for iterate(0u, $dim) |i| + { + for iterate(0u, $dim) |j| + { res.at[i] = res.at[i] + other.at[j] * self.at(j, i); } + } + + res + } + } + ) +) + +macro_rules! transform_impl( + ($t: ident, $v: ident) => ( + impl + Transform<$v> for $t + { + #[inline] + fn transform_vec(&self, v: &$v) -> $v + { self.rmul(v) } + + #[inline] + fn inv_transform(&self, v: &$v) -> $v + { self.inverse().transform_vec(v) } + } + ) +) + +macro_rules! inv_impl( + ($t: ident, $dim: expr) => ( + impl + Inv for $t + { + #[inline] + fn inverse(&self) -> $t + { + let mut res : $t = copy *self; + + res.invert(); + + res + } + + fn invert(&mut self) + { + let mut res: $t = One::one(); + let _0N: N = Zero::zero(); + + // inversion using Gauss-Jordan elimination + for iterate(0u, $dim) |k| + { + // search a non-zero value on the k-th column + // FIXME: would it be worth it to spend some more time searching for the + // max instead? + + let mut n0 = k; // index of a non-zero entry + + while (n0 != $dim) + { + if self.at(n0, k) != _0N + { break; } + + n0 = n0 + 1; + } + + // swap pivot line + if n0 != k + { + for iterate(0u, $dim) |j| + { + let off_n0_j = self.offset(n0, j); + let off_k_j = self.offset(k, j); + + swap(self.mij, off_n0_j, off_k_j); + swap(res.mij, off_n0_j, off_k_j); + } + } + + let pivot = self.at(k, k); + + for iterate(k, $dim) |j| + { + let selfval = &(self.at(k, j) / pivot); + self.set(k, j, selfval); + } + + for iterate(0u, $dim) |j| + { + let resval = &(res.at(k, j) / pivot); + res.set(k, j, resval); + } + + for iterate(0u, $dim) |l| + { + if l != k + { + let normalizer = self.at(l, k); + + for iterate(k, $dim) |j| + { + let selfval = &(self.at(l, j) - self.at(k, j) * normalizer); + self.set(l, j, selfval); + } + + for iterate(0u, $dim) |j| + { + let resval = &(res.at(l, j) - res.at(k, j) * normalizer); + res.set(l, j, resval); + } + } + } + } + + *self = res; + } + } + ) +) + +macro_rules! transpose_impl( + ($t: ident, $dim: expr) => ( + impl Transpose for $t + { + #[inline] + fn transposed(&self) -> $t + { + let mut res = copy *self; + + res.transpose(); + + res + } + + fn transpose(&mut self) + { + for iterate(1u, $dim) |i| + { + for iterate(0u, $dim - 1) |j| + { + let off_i_j = self.offset(i, j); + let off_j_i = self.offset(j, i); + + swap(self.mij, off_i_j, off_j_i); + } + } + } + } + ) +) + +macro_rules! approx_eq_impl( + ($t: ident) => ( + impl> ApproxEq for $t + { + #[inline] + fn approx_epsilon() -> N + { ApproxEq::approx_epsilon::() } + + #[inline] + fn approx_eq(&self, other: &$t) -> bool + { + let mut zip = self.mij.iter().zip(other.mij.iter()); + + do zip.all |(a, b)| { a.approx_eq(b) } + } + + #[inline] + fn approx_eq_eps(&self, other: &$t, epsilon: &N) -> bool + { + let mut zip = self.mij.iter().zip(other.mij.iter()); + + do zip.all |(a, b)| { a.approx_eq_eps(b, epsilon) } + } + } + ) +) + +macro_rules! rand_impl( + ($t: ident, $param: ident, [ $($elem: ident)|+ ]) => ( + impl Rand for $t + { + #[inline] + fn rand($param: &mut R) -> $t + { $t::new([ $( $elem.gen(), )+ ]) } + } + ) +) + +#[deriving(ToStr)] +pub struct Mat1 +{ mij: [N, ..1 * 1] } + +mat_impl!(Mat1, 1) +one_impl!(Mat1, [ _1 ]) +zero_impl!(Mat1, [ _0 ]) +dim_impl!(Mat1, 1) +mat_indexing_impl!(Mat1, 1) +mul_impl!(Mat1, 1) +rmul_impl!(Mat1, Vec1, 1) +lmul_impl!(Mat1, Vec1, 1) +transform_impl!(Mat1, Vec1) +// inv_impl!(Mat1, 1) +transpose_impl!(Mat1, 1) +approx_eq_impl!(Mat1) +rand_impl!(Mat1, rng, [ rng ]) + +#[deriving(ToStr)] +pub struct Mat2 +{ mij: [N, ..2 * 2] } + +mat_impl!(Mat2, 2) +one_impl!(Mat2, [ _1 | _0 | + _0 | _1 ]) +zero_impl!(Mat2, [ _0 | _0 | + _0 | _0 ]) +dim_impl!(Mat2, 2) +mat_indexing_impl!(Mat2, 2) +mul_impl!(Mat2, 2) +rmul_impl!(Mat2, Vec2, 2) +lmul_impl!(Mat2, Vec2, 2) +transform_impl!(Mat2, Vec2) +// inv_impl!(Mat2, 2) +transpose_impl!(Mat2, 2) +approx_eq_impl!(Mat2) +rand_impl!(Mat2, rng, [ rng | rng | + rng | rng ]) + +#[deriving(ToStr)] +pub struct Mat3 +{ mij: [N, ..3 * 3] } + +mat_impl!(Mat3, 3) +one_impl!(Mat3, [ _1 | _0 | _0 | + _0 | _1 | _0 | + _0 | _0 | _1 ]) +zero_impl!(Mat3, [ _0 | _0 | _0 | + _0 | _0 | _0 | + _0 | _0 | _0 ]) +dim_impl!(Mat3, 3) +mat_indexing_impl!(Mat3, 3) +mul_impl!(Mat3, 3) +rmul_impl!(Mat3, Vec3, 3) +lmul_impl!(Mat3, Vec3, 3) +transform_impl!(Mat3, Vec3) +// inv_impl!(Mat3, 3) +transpose_impl!(Mat3, 3) +approx_eq_impl!(Mat3) +rand_impl!(Mat3, rng, [ rng | rng | rng | + rng | rng | rng | + rng | rng | rng]) + +#[deriving(ToStr)] +pub struct Mat4 +{ mij: [N, ..4 * 4] } + +mat_impl!(Mat4, 4) +one_impl!(Mat4, [ + _1 | _0 | _0 | _0 | + _0 | _1 | _0 | _0 | + _0 | _0 | _1 | _0 | + _0 | _0 | _0 | _1 + ]) +zero_impl!(Mat4, [ + _0 | _0 | _0 | _0 | + _0 | _0 | _0 | _0 | + _0 | _0 | _0 | _0 | + _0 | _0 | _0 | _0 + ]) +dim_impl!(Mat4, 4) +mat_indexing_impl!(Mat4, 4) +mul_impl!(Mat4, 4) +rmul_impl!(Mat4, Vec4, 4) +lmul_impl!(Mat4, Vec4, 4) +transform_impl!(Mat4, Vec4) +inv_impl!(Mat4, 4) +transpose_impl!(Mat4, 4) +approx_eq_impl!(Mat4) +rand_impl!(Mat4, rng, [ + rng | rng | rng | rng | + rng | rng | rng | rng | + rng | rng | rng | rng | + rng | rng | rng | rng + ]) + +#[deriving(ToStr)] +pub struct Mat5 +{ mij: [N, ..5 * 5] } + +mat_impl!(Mat5, 5) +one_impl!(Mat5, [ + _1 | _0 | _0 | _0 | _0 | + _0 | _1 | _0 | _0 | _0 | + _0 | _0 | _1 | _0 | _0 | + _0 | _0 | _0 | _1 | _0 | + _0 | _0 | _0 | _0 | _1 + ]) +zero_impl!(Mat5, [ + _0 | _0 | _0 | _0 | _0 | + _0 | _0 | _0 | _0 | _0 | + _0 | _0 | _0 | _0 | _0 | + _0 | _0 | _0 | _0 | _0 | + _0 | _0 | _0 | _0 | _0 + ]) +dim_impl!(Mat5, 5) +mat_indexing_impl!(Mat5, 5) +mul_impl!(Mat5, 5) +rmul_impl!(Mat5, Vec5, 5) +lmul_impl!(Mat5, Vec5, 5) +transform_impl!(Mat5, Vec5) +inv_impl!(Mat5, 5) +transpose_impl!(Mat5, 5) +approx_eq_impl!(Mat5) +rand_impl!(Mat5, rng, [ + rng | rng | rng | rng | rng | + rng | rng | rng | rng | rng | + rng | rng | rng | rng | rng | + rng | rng | rng | rng | rng | + rng | rng | rng | rng | rng + ]) + +#[deriving(ToStr)] +pub struct Mat6 +{ mij: [N, ..6 * 6] } + +mat_impl!(Mat6, 6) +one_impl!(Mat6, [ + _1 | _0 | _0 | _0 | _0 | _0 | + _0 | _1 | _0 | _0 | _0 | _0 | + _0 | _0 | _1 | _0 | _0 | _0 | + _0 | _0 | _0 | _1 | _0 | _0 | + _0 | _0 | _0 | _0 | _1 | _0 | + _0 | _0 | _0 | _0 | _0 | _1 + ]) +zero_impl!(Mat6, [ + _0 | _0 | _0 | _0 | _0 | _0 | + _0 | _0 | _0 | _0 | _0 | _0 | + _0 | _0 | _0 | _0 | _0 | _0 | + _0 | _0 | _0 | _0 | _0 | _0 | + _0 | _0 | _0 | _0 | _0 | _0 | + _0 | _0 | _0 | _0 | _0 | _0 + ]) +dim_impl!(Mat6, 6) +mat_indexing_impl!(Mat6, 6) +mul_impl!(Mat6, 6) +rmul_impl!(Mat6, Vec6, 6) +lmul_impl!(Mat6, Vec6, 6) +transform_impl!(Mat6, Vec6) +inv_impl!(Mat6, 6) +transpose_impl!(Mat6, 6) +approx_eq_impl!(Mat6) +rand_impl!(Mat6, rng, [ + rng | rng | rng | rng | rng | rng | + rng | rng | rng | rng | rng | rng | + rng | rng | rng | rng | rng | rng | + rng | rng | rng | rng | rng | rng | + rng | rng | rng | rng | rng | rng | + rng | rng | rng | rng | rng | rng + ]) + +// some specializations: +impl +Inv for Mat1 +{ + #[inline] + fn inverse(&self) -> Mat1 + { + let mut res : Mat1 = copy *self; + + res.invert(); + + res + } + + #[inline] + fn invert(&mut self) + { + assert!(!self.mij[0].is_zero()); + + self.mij[0] = One::one::() / self.mij[0] + } +} + +impl +Inv for Mat2 +{ + #[inline] + fn inverse(&self) -> Mat2 + { + let mut res : Mat2 = copy *self; + + res.invert(); + + res + } + + #[inline] + fn invert(&mut self) + { + let det = self.mij[0 * 2 + 0] * self.mij[1 * 2 + 1] - self.mij[1 * 2 + 0] * self.mij[0 * 2 + 1]; + + assert!(!det.is_zero()); + + *self = Mat2::new([self.mij[1 * 2 + 1] / det , -self.mij[0 * 2 + 1] / det, + -self.mij[1 * 2 + 0] / det, self.mij[0 * 2 + 0] / det]) + } +} + +impl +Inv for Mat3 +{ + #[inline] + fn inverse(&self) -> Mat3 + { + let mut res = copy *self; + + res.invert(); + + res + } + + #[inline] + fn invert(&mut self) + { + let minor_m12_m23 = self.mij[1 * 3 + 1] * self.mij[2 * 3 + 2] - self.mij[2 * 3 + 1] * self.mij[1 * 3 + 2]; + let minor_m11_m23 = self.mij[1 * 3 + 0] * self.mij[2 * 3 + 2] - self.mij[2 * 3 + 0] * self.mij[1 * 3 + 2]; + let minor_m11_m22 = self.mij[1 * 3 + 0] * self.mij[2 * 3 + 1] - self.mij[2 * 3 + 0] * self.mij[1 * 3 + 1]; + + let det = self.mij[0 * 3 + 0] * minor_m12_m23 + - self.mij[0 * 3 + 1] * minor_m11_m23 + + self.mij[0 * 3 + 2] * minor_m11_m22; + + assert!(!det.is_zero()); + + *self = Mat3::new( [ + (minor_m12_m23 / det), + ((self.mij[0 * 3 + 2] * self.mij[2 * 3 + 1] - self.mij[2 * 3 + 2] * self.mij[0 * 3 + 1]) / det), + ((self.mij[0 * 3 + 1] * self.mij[1 * 3 + 2] - self.mij[1 * 3 + 1] * self.mij[0 * 3 + 2]) / det), + + (-minor_m11_m23 / det), + ((self.mij[0 * 3 + 0] * self.mij[2 * 3 + 2] - self.mij[2 * 3 + 0] * self.mij[0 * 3 + 2]) / det), + ((self.mij[0 * 3 + 2] * self.mij[1 * 3 + 0] - self.mij[1 * 3 + 2] * self.mij[0 * 3 + 0]) / det), + + (minor_m11_m22 / det), + ((self.mij[0 * 3 + 1] * self.mij[2 * 3 + 0] - self.mij[2 * 3 + 1] * self.mij[0 * 3 + 0]) / det), + ((self.mij[0 * 3 + 0] * self.mij[1 * 3 + 1] - self.mij[1 * 3 + 0] * self.mij[0 * 3 + 1]) / det) + ] ) + } +} diff --git a/src/nalgebra.rc b/src/nalgebra.rc index dcaef33b..ffd783b6 100644 --- a/src/nalgebra.rc +++ b/src/nalgebra.rc @@ -16,24 +16,7 @@ extern mod std; pub mod vec; - -/// 1-dimensional linear algebra. -pub mod dim1 -{ - pub mod mat1; -} - -/// 2-dimensional linear algebra. -pub mod dim2 -{ - pub mod mat2; -} - -/// 3-dimensional linear algebra. -pub mod dim3 -{ - pub mod mat3; -} +pub mod mat; /// n-dimensional linear algebra (slower). pub mod ndim diff --git a/src/tests/mat.rs b/src/tests/mat.rs index b6fd00bf..aa34c711 100644 --- a/src/tests/mat.rs +++ b/src/tests/mat.rs @@ -11,11 +11,7 @@ use traits::rotation::{Rotation, Rotatable}; #[test] use vec::Vec1; #[test] -use dim1::mat1::Mat1; -#[test] -use dim2::mat2::Mat2; -#[test] -use dim3::mat3::Mat3; +use mat::{Mat1, Mat2, Mat3}; #[test] use adaptors::rotmat::Rotmat;