la-stack 0.4.3

Fast, stack-allocated linear algebra for fixed dimensions
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

//! Regression tests for bugs caught in public API behavior.

use la_stack::prelude::*;

#[test]
#[cfg(feature = "exact")]
fn det_exact_f64_preserves_min_positive_subnormal() -> Result<(), LaError> {
    let tiny = f64::from_bits(1);
    let m = Matrix::<1>::try_from_rows([[tiny]])?;

    assert_eq!(m.det_exact_f64().unwrap().to_bits(), tiny.to_bits());
    Ok(())
}

#[test]
#[cfg(feature = "exact")]
fn det_exact_f64_strict_vs_rounded_inexact_det() -> Result<(), LaError> {
    let m = Matrix::<2>::try_from_rows([[1.0 + f64::EPSILON, 0.0], [0.0, 1.0 - f64::EPSILON]])?;

    assert_eq!(
        m.det_exact_f64(),
        Err(LaError::Unrepresentable {
            index: None,
            reason: UnrepresentableReason::RequiresRounding,
        })
    );
    assert_eq!(
        m.det_exact_rounded_f64().unwrap().to_bits(),
        1.0f64.to_bits()
    );
    Ok(())
}

#[test]
#[cfg(feature = "exact")]
fn solve_exact_f64_strict_vs_rounded_non_dyadic() -> Result<(), LaError> {
    let a = Matrix::<1>::try_from_rows([[3.0]])?;
    let b = Vector::<1>::try_new([1.0])?;

    assert_eq!(
        a.solve_exact_f64(b),
        Err(LaError::Unrepresentable {
            index: Some(0),
            reason: UnrepresentableReason::RequiresRounding,
        })
    );
    assert_eq!(
        a.solve_exact_rounded_f64(b).unwrap().into_array()[0].to_bits(),
        (1.0f64 / 3.0).to_bits()
    );
    Ok(())
}

#[test]
#[cfg(feature = "exact")]
fn requires_rounding_error_can_fall_back_to_rounded_solve() -> Result<(), LaError> {
    let a = Matrix::<1>::try_from_rows([[3.0]])?;
    let b = Vector::<1>::try_new([1.0])?;

    let rounded = match a.solve_exact_f64(b) {
        Ok(x) => x,
        Err(err) if err.requires_rounding() => a.solve_exact_rounded_f64(b)?,
        Err(err) => return Err(err),
    };

    assert_eq!(rounded.into_array()[0].to_bits(), (1.0f64 / 3.0).to_bits());
    Ok(())
}

#[test]
fn det_direct_skips_zero_coefficient_terms_that_would_overflow() -> Result<(), LaError> {
    let d3 = Matrix::<3>::try_from_rows([
        [1.0, 0.0, 0.0],
        [1.0e300, 1.0, 1.0e300],
        [1.0e300, 0.0, 1.0e300],
    ])?;
    assert_eq!(d3.det_direct(), Ok(Some(1.0e300)));

    let d4 = Matrix::<4>::try_from_rows([
        [0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0],
        [1.0e300, 0.0, 1.0e300, 1.0e300],
        [1.0e300, 0.0, 1.0e300, -1.0e300],
    ])?;
    assert_eq!(d4.det_direct(), Ok(Some(0.0)));

    Ok(())
}

#[test]
fn det_errbound_skips_zero_coefficient_terms_that_would_overflow() -> Result<(), LaError> {
    let d3 = Matrix::<3>::try_from_rows([
        [1.0, 0.0, 0.0],
        [1.0e300, 1.0, 1.0e300],
        [1.0e300, 0.0, 1.0e300],
    ])?;
    assert_eq!(d3.det_errbound(), Ok(Some(ERR_COEFF_3 * 1.0e300)));

    let d4 = Matrix::<4>::try_from_rows([
        [0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0],
        [1.0e300, 0.0, 1.0e300, 1.0e300],
        [1.0e300, 0.0, 1.0e300, -1.0e300],
    ])?;
    assert_eq!(d4.det_errbound(), Ok(Some(0.0)));

    Ok(())
}