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
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148

use std::ops::MulAssign;

use num_traits::Float;

/// Computes one DQD transform in ping-pong form, with protection against
/// underflow and overflow.
///
/// # Panics
///
/// * The number of elements in `z` is not a multiple of 4.
/// * The length of `z` is not `(n0 + 1) * 4`.
/// * `n0` is not greater than `i0 + 1`.
#[allow(dead_code)]
pub(crate) fn lasq6<A, const PP: usize>(i0: usize, n0: usize, z: &mut [A]) -> (A, A, A, A, A, A)
where
    A: Float + MulAssign,
{
    assert!(PP <= 1, "`PP` must be either 0 or 1.");
    assert_eq!(z.len(), (n0 + 1) * 4);
    assert!(i0 + 1 < n0);

    let mut j4 = 4 * i0 + PP;
    let mut e_min = z[j4 + 4];
    let mut d = z[j4];
    let mut d_min = d;

    if n0 >= 3 {
        for j4 in (4 * i0 + 3..=4 * (n0 - 3) + 3).step_by(4) {
            z[j4 - PP - 2] = d + z[j4 + PP - 1];
            if z[j4 - PP - 2] == A::zero() {
                z[j4 - PP] = A::zero();
                d = z[j4 + PP + 1];
                d_min = d;
                e_min = A::zero();
            } else if A::min_positive_value() * z[j4 + PP + 1] < z[j4 - PP - 2]
                && A::min_positive_value() * z[j4 - PP - 2] < z[j4 + PP + 1]
            {
                let tmp = z[j4 + PP + 1] / z[j4 - PP - 2];
                z[j4 - PP] = z[j4 + PP - 1] * tmp;
                d *= tmp;
            } else {
                z[j4 - PP] = z[j4 + PP + 1] * (z[j4 + PP - 1] / z[j4 - PP - 2]);
                d = z[j4 + PP + 1] * (d / z[j4 - PP - 2]);
            }
            if d < d_min {
                d_min = d;
            }
            let e = z[j4 - PP];
            if e < e_min {
                e_min = e;
            }
        }
    }

    let d_nm2 = d;
    let d_min_2 = d_min;
    j4 = 4 * (n0 - 2) + 3 - PP;
    let mut j4_p2 = j4 + 2 * PP - 1;
    z[j4 - 2] = d_nm2 + z[j4_p2];
    let d_nm1;
    if z[j4 - 2] == A::zero() {
        z[j4] = A::zero();
        d_nm1 = z[j4_p2 + 2];
        d_min = d_nm1;
        e_min = A::zero();
    } else if A::min_positive_value() * z[j4_p2 + 2] < z[j4 - 2]
        && A::min_positive_value() * z[j4 - 2] < z[j4_p2 + 2]
    {
        let tmp = z[j4_p2 + 2] / z[j4 - 2];
        z[j4] = z[j4_p2] * tmp;
        d_nm1 = d_nm2 * tmp;
    } else {
        z[j4] = z[j4_p2 + 2] * (z[j4_p2] / z[j4 - 2]);
        d_nm1 = z[j4_p2 + 2] * (d_nm2 / z[j4 - 2]);
    }
    if d_nm1 < d_min {
        d_min = d_nm1;
    }

    let d_min_1 = d_min;
    j4 += 4;
    j4_p2 = j4 + 2 * PP - 1;
    z[j4 - 2] = d_nm1 + z[j4_p2];
    let d_n;
    if z[j4 - 2] == A::zero() {
        z[j4] = A::zero();
        d_n = z[j4_p2 + 2];
        d_min = d_n;
        e_min = A::zero();
    } else if A::min_positive_value() * z[j4_p2 + 2] < z[j4 - 2]
        && A::min_positive_value() * z[j4 - 2] < z[j4_p2 + 2]
    {
        let tmp = z[j4_p2 + 2] / z[j4 - 2];
        z[j4] = z[j4_p2] * tmp;
        d_n = d_nm1 * tmp;
    } else {
        z[j4] = z[j4_p2 + 2] * (z[j4_p2] / z[j4 - 2]);
        d_n = z[j4_p2 + 2] * (d_nm1 / z[j4 - 2]);
    }
    if d_n < d_min {
        d_min = d_n;
    }

    z[j4 + 2] = d_n;
    z[4 * n0 + 3 - PP] = e_min;
    (d_min, d_min_1, d_min_2, d_n, d_nm1, d_nm2)
}

#[cfg(test)]
mod tests {
    use approx::assert_ulps_eq;

    #[test]
    fn lasq6_ping() {
        let mut z = [
            2.0, -1.0, 3.0, 1.0, -1.0, 3.0, 2.0, 2.0, 1.0, 1.0, -1.0, -1.0, 2.0, -3.0, -1.0, 2.0,
        ];
        let res = super::lasq6::<_, 0>(0, 3, &mut z);
        assert_eq!(
            z,
            [
                2.0, 5.0, 3.0, -0.6, -1.0, 1.6, 2.0, 1.25, 1.0, -1.25, -1.0, 1.6, 2.0, 0.4, -1.0,
                -1.0,
            ]
        );
        assert_eq!(res, (-0.4, -0.4, -0.4, 0.4, -0.25, -0.4));
    }

    #[test]
    fn lasq6_pong() {
        let mut z = [
            2.0, -1.0, 3.0, 1.0, -1.0, 3.0, 2.0, 2.0, 1.0, 1.0, -1.0, -1.0, 2.0, -3.0, -1.0, 2.0,
        ];
        let z_after = [
            0.0, -1.0, 0.0, 1.0, 5.0, 3.0, 0.4, 2.0, -0.4, 1.0, -7.5, -1.0, 4.5, -3.0, 0.0, 2.0,
        ];
        let res = super::lasq6::<_, 1>(0, 3, &mut z);
        for (&actual, &expected) in z.iter().zip(z_after.iter()) {
            assert_ulps_eq!(actual, expected);
        }
        assert_ulps_eq!(res.0, 0.6);
        assert_ulps_eq!(res.1, 0.6);
        assert_ulps_eq!(res.2, 3.0);
        assert_ulps_eq!(res.3, 4.5);
        assert_ulps_eq!(res.4, 0.6);
        assert_ulps_eq!(res.5, 3.0);
    }
}