la-stack 0.5.0

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
109
110
111
112
113

#![forbid(unsafe_code)]

//! Independently checked fixtures shared by angle benchmarks and their tests.

use core::f64::consts::{FRAC_PI_4, PI};

use la_stack::{Vector, VectorAngle};

use crate::bench_utils::OrAbort;

/// The ordinary max-scaled Kahan control is valid for these vetted fixtures.
/// It is deliberately not a general-purpose API: pre-normalization can lose
/// subnormal differences. No input checks are included in this kernel control.
pub(crate) fn stable_control<const D: usize>(left: &[f64; D], right: &[f64; D]) -> f64 {
    let maximum = |values: &[f64; D]| values.iter().fold(0.0_f64, |s, x| s.max(x.abs()));
    let left_scale = maximum(left);
    let right_scale = maximum(right);
    let mut left_norm = 0.0_f64;
    let mut right_norm = 0.0_f64;
    for (&x, &y) in left.iter().zip(right) {
        left_norm = left_norm.hypot(x / left_scale);
        right_norm = right_norm.hypot(y / right_scale);
    }
    let mut difference = 0.0_f64;
    let mut sum = 0.0_f64;
    for (&x, &y) in left.iter().zip(right) {
        let x = (x / left_scale) * right_norm;
        let y = (y / right_scale) * left_norm;
        difference = difference.hypot(x - y);
        sum = sum.hypot(x + y);
    }
    (2.0 * difference * sum).atan2((sum - difference) * (sum + difference))
}

/// Private construction makes numerical validation a prerequisite to timing.
pub(crate) struct AngleInput<const D: usize> {
    name: &'static str,
    left: Vector<D>,
    right: Vector<D>,
}

impl<const D: usize> AngleInput<D> {
    fn new(name: &'static str, left: [f64; D], right: [f64; D], expected: f64) -> Self {
        let left = Vector::try_new(left).or_abort("finite benchmark fixture");
        let right = Vector::try_new(right).or_abort("finite benchmark fixture");
        let tolerance = expected.abs() * (8.0 * f64::EPSILON);
        for actual in [
            left.as_array()
                .as_slice()
                .angle(right.as_array())
                .or_abort("valid angle"),
            left.angle(&right).or_abort("valid vector angle"),
            stable_control(left.as_array(), right.as_array()),
        ] {
            assert!(
                (actual - expected).abs() <= tolerance,
                "{name}: {actual:e} vs {expected:e}"
            );
        }
        Self { name, left, right }
    }

    pub(crate) const fn name(&self) -> &'static str {
        self.name
    }
    pub(crate) const fn left(&self) -> &Vector<D> {
        &self.left
    }
    pub(crate) const fn right(&self) -> &Vector<D> {
        &self.right
    }
}

#[expect(
    clippy::cast_precision_loss,
    reason = "benchmark dimensions are 3 through 6"
)]
pub(crate) fn fixtures<const D: usize>() -> [AngleInput<D>; 5] {
    assert!((3..=6).contains(&D));
    let planar = |x, y| {
        let mut values = [0.0; D];
        values[0] = x;
        values[D - 1] = y;
        values
    };
    let mut dense_right = [1.0; D];
    dense_right[0] = -1.0;
    // Integer dot is D-2; D-1 nonzero exterior minors each have magnitude 2.
    let dense_reference = (4.0 * (D - 1) as f64).sqrt().atan2((D - 2) as f64);
    let t = 2.0_f64.powi(-30);
    [
        AngleInput::new("dense", [1.0; D], dense_right, dense_reference),
        AngleInput::new("near_parallel", planar(1.0, 0.0), planar(1.0, t), t.atan()),
        AngleInput::new(
            "near_antipodal",
            planar(1.0, 0.0),
            planar(-1.0, t),
            PI - t.atan(),
        ),
        AngleInput::new(
            "mixed_scale",
            planar(f64::MAX, 0.0),
            planar(f64::MIN_POSITIVE, f64::MIN_POSITIVE),
            FRAC_PI_4,
        ),
        AngleInput::new(
            "subnormal",
            planar(1.0, 0.0),
            planar(1.0, f64::from_bits(1)),
            f64::from_bits(1),
        ),
    ]
}