1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100

use std::ops::Mul;

#[derive(Debug, Clone)]
pub struct Matrix4 {
    a: (f64, f64, f64, f64),
    b: (f64, f64, f64, f64),
    c: (f64, f64, f64, f64),
    d: (f64, f64, f64, f64),
}

#[derive(Debug, Clone, Copy)]
pub struct Matrix2 {
    a: f64,
    b: f64,
    c: f64,
    d: f64,
}
impl Matrix2 {
    pub fn new(r0: (f64, f64), r1: (f64, f64)) -> Self {
        Self {
            a: r0.0,
            b: r0.1,
            c: r1.0,
            d: r1.1,
        }
    }

    pub fn diagonal(a: f64, b: f64) -> Self {
        Self::new((a, 0.0), (0.0, b))
    }

    pub fn from_lower_triangular(x: (f64, f64, f64)) -> Self {
        Self::new((x.0, 0.0), (x.1, x.2))
    }

    pub fn transpose(&self) -> Self {
        Self {
            a: self.a,
            b: self.c,
            c: self.b,
            d: self.d,
        }
    }

    pub fn det(&self) -> f64 {
        self.a * self.d - self.b * self.c
    }

    pub fn inverse(&self) -> Self {
        let det = self.det();
        Self {
            a: self.d / det,
            b: -self.b / det,
            c: -self.c / det,
            d: self.a / det,
        }
    }

    pub fn lower_triangular(&self) -> impl Iterator<Item = f64> {
        use std::iter::once;
        once(self.a).chain(once(self.c)).chain(once(self.d))
    }

    pub fn lower_triangular_tuple(&self) -> (f64, f64, f64) {
        (self.a, self.c, self.d)
    }
}
impl Mul<(f64, f64)> for Matrix2 {
    type Output = (f64, f64);

    fn mul(self, rhs: (f64, f64)) -> Self::Output {
        let x = self.a * rhs.0 + self.b * rhs.1;
        let y = self.c * rhs.0 + self.d * rhs.1;
        (x, y)
    }
}
impl Mul<Matrix2> for Matrix2 {
    type Output = Matrix2;

    fn mul(self, rhs: Matrix2) -> Self::Output {
        let c0 = self * (rhs.a, rhs.c);
        let c1 = self * (rhs.b, rhs.d);
        Self {
            a: c0.0,
            b: c1.0,
            c: c0.1,
            d: c1.1,
        }
    }
}

#[derive(Debug)]
pub struct Transpose<T>(pub T);
impl Mul<(f64, f64)> for Transpose<(f64, f64)> {
    type Output = f64;

    fn mul(self, rhs: (f64, f64)) -> Self::Output {
        (self.0).0 * rhs.0 + (self.0).1 * rhs.1
    }
}