lair 0.8.0

Linear algebra in Rust
Documentation 
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
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125

//! Matrix functions and special matrices.

use ndarray::Array2;

use crate::{InvalidInput, Scalar};

/// Constructs a circulant matrix.
///
/// # Examples
///
/// ```
/// use lair::matrix::circulant;
///
/// let a = vec![1., 2., 3.];
/// let c = circulant(&a);
/// assert_eq!(c, ndarray::aview2(&[[1., 3., 2.], [2., 1., 3.], [3., 2., 1.]]));
/// ```
pub fn circulant<A>(a: &[A]) -> Array2<A>
where
    A: Copy,
{
    let mut x = Array2::<A>::uninit((a.len(), a.len()));
    unsafe {
        for (i, a_elem) in a.iter().enumerate() {
            for j in 0..a.len() {
                *(x[[(i + j) % a.len(), j]].as_mut_ptr()) = *a_elem;
            }
        }
        x.assume_init()
    }
}

/// Constructs a companion matrix.
///
/// # Errors
///
/// * [`InvalidInput::Shape`] if `a` contains less than two coefficients.
/// * [`InvalidInput::Value`] if `a[0]` is zero.
///
/// [`InvalidInput::Shape`]: ../enum.InvalidInput.html#variant.Shape
/// [`InvalidInput::Value`]: ../enum.InvalidInput.html#variant.Value
///
/// # Examples
///
/// ```
/// use lair::matrix::companion;
///
/// let a = vec![1., -10., 31., -30.];
/// let c = companion(&a).expect("valid input");
/// assert_eq!(c, ndarray::aview2(&[[10., -31., 30.], [1., 0., 0.], [0., 1., 0.]]));
/// ```
pub fn companion<A>(a: &[A]) -> Result<Array2<A>, InvalidInput>
where
    A: Scalar,
{
    if a.len() < 2 {
        return Err(InvalidInput::Shape(format!(
            "input polynomial has {} coefficient; expected at least two",
            a.len()
        )));
    }
    if a[0] == A::zero() {
        return Err(InvalidInput::Value(format!(
            "invalid first coefficient {}, expected a non-zero value",
            a[0]
        )));
    }

    let mut matrix = Array2::<A>::zeros((a.len() - 1, a.len() - 1));
    for (mv, av) in matrix.row_mut(0).into_iter().zip(a.iter().skip(1)) {
        *mv = -*av / a[0];
    }
    for i in 1..a.len() - 1 {
        matrix[[i, i - 1]] = A::one();
    }
    Ok(matrix)
}

#[cfg(test)]
mod tests {
    use ndarray::aview2;
    use num_complex::{Complex32, Complex64};

    #[test]
    fn circulant() {
        let a = Vec::<f32>::new();
        let c = super::circulant(&a);
        assert_eq!(c.shape(), [0, 0]);

        let a = vec![Complex32::new(1., 2.)];
        let c = super::circulant(&a);
        assert_eq!(c, aview2(&[[Complex32::new(1., 2.)]]));

        let a = vec![1., 2., 4.];
        let c = super::circulant(&a);
        assert_eq!(c, aview2(&[[1., 4., 2.], [2., 1., 4.], [4., 2., 1.]]));
    }

    #[test]
    fn companion() {
        let a = Vec::<f32>::new();
        assert!(super::companion(&a).is_err());

        let a = vec![Complex64::new(1., 2.)];
        assert!(super::companion(&a).is_err());

        let a = vec![
            Complex32::new(0., 0.),
            Complex32::new(1., 1.),
            Complex32::new(2., 2.),
        ];
        assert!(super::companion(&a).is_err());

        let a = vec![Complex32::new(1., 2.), Complex32::new(3., 4.)];
        let c = super::companion(&a).expect("valid input");
        assert_eq!(c, aview2(&[[Complex32::new(-2.2, 0.4)]]));

        let a = vec![2., -4., 8., -10.];
        let c = super::companion(&a).expect("valid input");
        assert_eq!(
            c,
            aview2(&[[2.0, -4.0, 5.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]])
        );
    }
}