la_stack/vector.rs
1#![forbid(unsafe_code)]
2
3//! Fixed-size, stack-allocated vectors.
4
5use core::hint::cold_path;
6
7use crate::{ArithmeticOperation, LaError};
8
9/// Finite fixed-size vector of length `D`, stored inline.
10///
11/// Public construction rejects NaN and infinity through [`try_new`](Self::try_new),
12/// and the storage field is private, so a `Vector` value carries the invariant
13/// that every stored entry is finite. Algorithms therefore do not re-scan stored
14/// entries at every use; user-visible non-finite errors come from construction
15/// boundaries or from values computed during arithmetic, such as overflowed
16/// accumulators.
17///
18/// Direct field construction is intentionally unavailable to downstream callers:
19///
20/// ```compile_fail
21/// use la_stack::Vector;
22///
23/// let _ = Vector::<2> {
24/// data: [1.0, f64::NAN],
25/// };
26/// ```
27#[must_use]
28#[derive(Clone, Copy, Debug, PartialEq)]
29pub struct Vector<const D: usize> {
30 data: [f64; D],
31}
32
33impl<const D: usize> Vector<D> {
34 /// Test-only infallible constructor for finite literal fixtures.
35 #[cfg(test)]
36 #[inline]
37 pub(crate) const fn new(data: [f64; D]) -> Self {
38 match Self::try_new(data) {
39 Ok(vector) => vector,
40 Err(_) => panic!("Vector::new requires finite entries"),
41 }
42 }
43
44 /// Try to create a finite vector from a backing array.
45 ///
46 /// This is the public raw-storage boundary for vectors. Successful
47 /// construction makes the returned [`Vector`] a finite-storage proof.
48 ///
49 /// # Examples
50 /// ```
51 /// use la_stack::prelude::*;
52 ///
53 /// # fn main() -> Result<(), LaError> {
54 /// let v = Vector::<3>::try_new([1.0, 2.0, 3.0])?;
55 /// assert_eq!(v.into_array(), [1.0, 2.0, 3.0]);
56 /// # Ok(())
57 /// # }
58 /// ```
59 ///
60 /// # Errors
61 /// Returns [`LaError::NonFinite`] with the first offending entry index when
62 /// `data` contains NaN or infinity.
63 #[inline]
64 pub const fn try_new(data: [f64; D]) -> Result<Self, LaError> {
65 if let Some(index) = Self::first_non_finite_entry(&data) {
66 Err(LaError::non_finite_input_vector(index))
67 } else {
68 Ok(Self { data })
69 }
70 }
71
72 /// Finalize vector storage produced by an arithmetic operation.
73 ///
74 /// Keeping this validation in the type that owns the finite-storage
75 /// invariant prevents a new computation path from accidentally turning raw
76 /// non-finite storage into a [`Vector`].
77 #[inline]
78 pub(crate) const fn from_computation(
79 data: [f64; D],
80 operation: ArithmeticOperation,
81 ) -> Result<Self, LaError> {
82 if let Some(index) = Self::first_non_finite_entry(&data) {
83 Err(LaError::non_finite_computation_step(operation, index))
84 } else {
85 Ok(Self { data })
86 }
87 }
88
89 /// Return the first non-finite stored entry in index order.
90 ///
91 /// Used by the public raw-storage boundary to report the first offending
92 /// index with [`LaError::NonFinite`].
93 const fn first_non_finite_entry(data: &[f64; D]) -> Option<usize> {
94 let mut i = 0;
95 while i < D {
96 if !data[i].is_finite() {
97 return Some(i);
98 }
99 i += 1;
100 }
101 None
102 }
103
104 /// All-zeros finite vector.
105 ///
106 /// # Examples
107 /// ```
108 /// use la_stack::prelude::*;
109 ///
110 /// let z = Vector::<2>::zero();
111 /// assert_eq!(z.into_array(), [0.0, 0.0]);
112 /// ```
113 #[inline]
114 pub const fn zero() -> Self {
115 Self { data: [0.0; D] }
116 }
117
118 /// Borrow the finite backing array.
119 ///
120 /// # Examples
121 /// ```
122 /// use la_stack::prelude::*;
123 ///
124 /// # fn main() -> Result<(), LaError> {
125 /// let v = Vector::<2>::try_new([1.0, -2.0])?;
126 /// assert_eq!(v.as_array(), &[1.0, -2.0]);
127 /// # Ok(())
128 /// # }
129 /// ```
130 #[inline]
131 #[must_use]
132 pub const fn as_array(&self) -> &[f64; D] {
133 &self.data
134 }
135
136 /// Consume and return the finite backing array.
137 ///
138 /// # Examples
139 /// ```
140 /// use la_stack::prelude::*;
141 ///
142 /// # fn main() -> Result<(), LaError> {
143 /// let v = Vector::<2>::try_new([1.0, 2.0])?;
144 /// let a = v.into_array();
145 /// assert_eq!(a, [1.0, 2.0]);
146 /// # Ok(())
147 /// # }
148 /// ```
149 #[inline]
150 #[must_use]
151 pub const fn into_array(self) -> [f64; D] {
152 self.data
153 }
154
155 /// Dot product.
156 ///
157 /// Terms are accumulated in `f64` using [`f64::mul_add`] at each index.
158 /// Intermediate rounding occurs, and this method does not provide a
159 /// certified absolute rounding bound for the returned dot product. Raw
160 /// `Vector` values are finite by construction, so this method only checks
161 /// whether the accumulation overflows to NaN or infinity.
162 ///
163 /// # Examples
164 /// ```
165 /// use la_stack::prelude::*;
166 ///
167 /// # fn main() -> Result<(), LaError> {
168 /// let a = Vector::<3>::try_new([1.0, 2.0, 3.0])?;
169 /// let b = Vector::<3>::try_new([-2.0, 0.5, 4.0])?;
170 /// assert!((a.dot(&b)? - 11.0).abs() <= 1e-12);
171 /// # Ok(())
172 /// # }
173 /// ```
174 ///
175 /// # Errors
176 /// Returns [`LaError::NonFinite`] when the accumulated dot product overflows
177 /// to NaN or infinity.
178 #[inline]
179 pub const fn dot(&self, other: &Self) -> Result<f64, LaError> {
180 self.dot_with_operation(other, ArithmeticOperation::VectorDotProduct)
181 }
182
183 /// Accumulate a dot product while retaining the public operation that owns it.
184 const fn dot_with_operation(
185 &self,
186 other: &Self,
187 operation: ArithmeticOperation,
188 ) -> Result<f64, LaError> {
189 let lhs = self.as_array();
190 let rhs = other.as_array();
191 let mut acc = 0.0;
192 let mut i = 0;
193 while i < D {
194 acc = lhs[i].mul_add(rhs[i], acc);
195 i += 1;
196 }
197 if acc.is_finite() {
198 Ok(acc)
199 } else {
200 cold_path();
201 Err(Self::dot_non_finite_error(lhs, rhs, operation))
202 }
203 }
204
205 /// Replay a non-finite dot product to locate the first failing step.
206 ///
207 /// This runs only after the success-path traversal has produced a non-finite
208 /// final accumulator. Stored entries are finite, so once a fused multiply-add
209 /// produces a non-finite accumulator, later steps cannot make it finite again.
210 /// Replaying the same left-to-right operations must therefore find the first
211 /// failing index.
212 #[cold]
213 const fn dot_non_finite_error(
214 lhs: &[f64; D],
215 rhs: &[f64; D],
216 operation: ArithmeticOperation,
217 ) -> LaError {
218 let mut acc = 0.0;
219 let mut i = 0;
220 let last = D.saturating_sub(1);
221 while i < last {
222 acc = lhs[i].mul_add(rhs[i], acc);
223 if !acc.is_finite() {
224 return LaError::non_finite_computation_step(operation, i);
225 }
226 i += 1;
227 }
228
229 LaError::non_finite_computation_step(operation, last)
230 }
231
232 /// Squared Euclidean norm.
233 ///
234 /// This is computed as `dot(self, self)`, so `norm2_sq` has the same
235 /// `f64` [`mul_add`](f64::mul_add) accumulation behavior as [`dot`](Self::dot).
236 /// Intermediate rounding occurs, and this method does not provide a
237 /// certified absolute rounding bound for the returned squared norm.
238 /// `Vector` values are finite by construction, so this method only checks
239 /// whether the accumulation overflows to NaN or infinity.
240 ///
241 /// # Examples
242 /// ```
243 /// use la_stack::prelude::*;
244 ///
245 /// # fn main() -> Result<(), LaError> {
246 /// let v = Vector::<3>::try_new([1.0, 2.0, 3.0])?;
247 /// assert!((v.norm2_sq()? - 14.0).abs() <= 1e-12);
248 /// # Ok(())
249 /// # }
250 /// ```
251 ///
252 /// # Errors
253 /// Returns [`LaError::NonFinite`] when the accumulated norm overflows to NaN
254 /// or infinity.
255 #[inline]
256 pub const fn norm2_sq(&self) -> Result<f64, LaError> {
257 self.dot_with_operation(self, ArithmeticOperation::VectorSquaredNorm)
258 }
259}
260
261impl<const D: usize> Default for Vector<D> {
262 #[inline]
263 fn default() -> Self {
264 Self::zero()
265 }
266}
267
268#[cfg(test)]
269mod tests {
270 use core::hint::black_box;
271
272 use approx::assert_abs_diff_eq;
273 use pastey::paste;
274
275 use super::*;
276
277 macro_rules! gen_vector_tests {
278 ($d:literal) => {
279 paste! {
280 #[test]
281 fn [<vector_new_as_array_into_array_ $d d>]() {
282 let arr = {
283 let mut arr = [0.0f64; $d];
284 let values = [1.0f64, 2.0, 3.0, 4.0, 5.0];
285 for (dst, src) in arr.iter_mut().zip(values.iter()) {
286 *dst = *src;
287 }
288 arr
289 };
290
291 let v = Vector::<$d>::new(arr);
292
293 for i in 0..$d {
294 assert_abs_diff_eq!(v.as_array()[i], arr[i], epsilon = 0.0);
295 }
296
297 let out = v.into_array();
298 for i in 0..$d {
299 assert_abs_diff_eq!(out[i], arr[i], epsilon = 0.0);
300 }
301 }
302
303 #[test]
304 fn [<vector_zero_as_array_into_array_default_ $d d>]() {
305 let z = Vector::<$d>::zero();
306 for &x in z.as_array() {
307 assert_abs_diff_eq!(x, 0.0, epsilon = 0.0);
308 }
309 for x in z.into_array() {
310 assert_abs_diff_eq!(x, 0.0, epsilon = 0.0);
311 }
312
313 let d = Vector::<$d>::default();
314 for x in d.into_array() {
315 assert_abs_diff_eq!(x, 0.0, epsilon = 0.0);
316 }
317 }
318
319 #[test]
320 fn [<vector_dot_and_norm2_sq_ $d d>]() {
321 // Use black_box to avoid constant-folding/inlining eliminating the actual dot loop,
322 // which can make coverage tools report the mul_add line as uncovered.
323
324 let a_arr = {
325 let mut arr = [0.0f64; $d];
326 let values = [1.0f64, 2.0, 3.0, 4.0, 5.0];
327 for (dst, src) in arr.iter_mut().zip(values.iter()) {
328 *dst = black_box(*src);
329 }
330 arr
331 };
332 let b_arr = {
333 let mut arr = [0.0f64; $d];
334 let values = [-2.0f64, 0.5, 4.0, -1.0, 2.0];
335 for (dst, src) in arr.iter_mut().zip(values.iter()) {
336 *dst = black_box(*src);
337 }
338 arr
339 };
340
341 let expected_dot = {
342 let mut acc = 0.0;
343 let mut i = 0;
344 while i < $d {
345 acc = a_arr[i].mul_add(b_arr[i], acc);
346 i += 1;
347 }
348 acc
349 };
350 let expected_norm2_sq = {
351 let mut acc = 0.0;
352 let mut i = 0;
353 while i < $d {
354 acc = a_arr[i].mul_add(a_arr[i], acc);
355 i += 1;
356 }
357 acc
358 };
359
360 let a = Vector::<$d>::new(black_box(a_arr));
361 let b = Vector::<$d>::new(black_box(b_arr));
362
363 // Call via (black_boxed) fn pointers to discourage inlining, improving line-level coverage
364 // attribution for the loop body.
365 let dot_fn: fn(&Vector<$d>, &Vector<$d>) -> Result<f64, LaError> =
366 black_box(Vector::<$d>::dot);
367 let norm2_sq_fn: fn(&Vector<$d>) -> Result<f64, LaError> =
368 black_box(Vector::<$d>::norm2_sq);
369
370 assert_abs_diff_eq!(
371 dot_fn(black_box(&a), black_box(&b)).unwrap(),
372 expected_dot,
373 epsilon = 1e-14
374 );
375 assert_abs_diff_eq!(
376 norm2_sq_fn(black_box(&a)).unwrap(),
377 expected_norm2_sq,
378 epsilon = 1e-14
379 );
380 }
381
382 #[test]
383 fn [<vector_try_new_rejects_non_finite_ $d d>]() {
384 for value in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
385 let mut data = [1.0f64; $d];
386 data[$d - 1] = value;
387 assert_eq!(
388 Vector::<$d>::try_new(data),
389 Err(LaError::non_finite_input_vector($d - 1))
390 );
391 }
392
393 let mut data = [1.0f64; $d];
394 data[0] = f64::INFINITY;
395 data[$d - 1] = f64::NAN;
396 assert_eq!(
397 Vector::<$d>::try_new(data),
398 Err(LaError::non_finite_input_vector(0))
399 );
400 }
401
402 #[test]
403 fn [<vector_from_computation_preserves_failure_provenance_ $d d>]() {
404 let mut data = [1.0f64; $d];
405 data[$d - 1] = f64::INFINITY;
406
407 assert_eq!(
408 Vector::<$d>::from_computation(
409 data,
410 ArithmeticOperation::LuSolve,
411 ),
412 Err(LaError::non_finite_computation_step(
413 ArithmeticOperation::LuSolve,
414 $d - 1,
415 ))
416 );
417 }
418
419 #[test]
420 fn [<vector_dot_and_norm2_sq_reject_overflow_ $d d>]() {
421 let mut a_arr = [1.0f64; $d];
422 a_arr[0] = f64::MAX;
423 let a = Vector::<$d>::new(a_arr);
424
425 let mut b_arr = [1.0f64; $d];
426 b_arr[0] = 2.0;
427 let b = Vector::<$d>::new(b_arr);
428
429 assert_eq!(
430 a.dot(&b),
431 Err(LaError::non_finite_computation_step(
432 ArithmeticOperation::VectorDotProduct,
433 0,
434 ))
435 );
436 assert_eq!(
437 a.norm2_sq(),
438 Err(LaError::non_finite_computation_step(
439 ArithmeticOperation::VectorSquaredNorm,
440 0,
441 ))
442 );
443 }
444
445 }
446 };
447 }
448
449 // Mirror delaunay-style multi-dimension tests.
450 gen_vector_tests!(1);
451 gen_vector_tests!(2);
452 gen_vector_tests!(3);
453 gen_vector_tests!(4);
454 gen_vector_tests!(5);
455
456 macro_rules! gen_vector_replay_tests {
457 ($d:literal) => {
458 paste! {
459 #[test]
460 fn [<vector_dot_and_norm2_sq_report_last_overflowing_step_ $d d>]() {
461 let mut dot_lhs = [1.0f64; $d];
462 dot_lhs[$d - 1] = f64::MAX;
463 let mut dot_rhs = [1.0f64; $d];
464 dot_rhs[$d - 1] = 2.0;
465 let dot_lhs = Vector::<$d>::new(dot_lhs);
466 let dot_rhs = Vector::<$d>::new(dot_rhs);
467
468 assert_eq!(
469 dot_lhs.dot(&dot_rhs),
470 Err(LaError::non_finite_computation_step(
471 ArithmeticOperation::VectorDotProduct,
472 $d - 1,
473 ))
474 );
475
476 let mut norm_data = [1.0f64; $d];
477 norm_data[$d - 1] = f64::MAX;
478 let vector = Vector::<$d>::new(norm_data);
479
480 assert_eq!(
481 vector.norm2_sq(),
482 Err(LaError::non_finite_computation_step(
483 ArithmeticOperation::VectorSquaredNorm,
484 $d - 1,
485 ))
486 );
487 }
488 }
489 };
490 }
491
492 gen_vector_replay_tests!(2);
493 gen_vector_replay_tests!(3);
494 gen_vector_replay_tests!(4);
495 gen_vector_replay_tests!(5);
496
497 macro_rules! gen_vector_const_eval_tests {
498 ($d:literal, $dot:literal, $norm2_sq:literal) => {
499 paste! {
500 #[test]
501 fn [<vector_dot_and_norm2_sq_const_eval_ $d d>]() {
502 const DOT: Result<f64, LaError> = Vector::<$d>::new([1.0; $d])
503 .dot(&Vector::<$d>::new([2.0; $d]));
504 const NORM2_SQ: Result<f64, LaError> =
505 Vector::<$d>::new([1.0; $d]).norm2_sq();
506
507 assert_eq!(DOT, Ok($dot));
508 assert_eq!(NORM2_SQ, Ok($norm2_sq));
509 }
510 }
511 };
512 }
513
514 gen_vector_const_eval_tests!(2, 4.0, 2.0);
515 gen_vector_const_eval_tests!(3, 6.0, 3.0);
516 gen_vector_const_eval_tests!(4, 8.0, 4.0);
517 gen_vector_const_eval_tests!(5, 10.0, 5.0);
518
519 #[test]
520 fn vector_dot_and_norm2_sq_overflow_const_eval() {
521 const DOT: Result<f64, LaError> =
522 Vector::<2>::new([f64::MAX; 2]).dot(&Vector::<2>::new([1.0; 2]));
523 const NORM2_SQ: Result<f64, LaError> = Vector::<2>::new([f64::MAX; 2]).norm2_sq();
524
525 assert_eq!(
526 DOT,
527 Err(LaError::non_finite_computation_step(
528 ArithmeticOperation::VectorDotProduct,
529 1,
530 ))
531 );
532 assert_eq!(
533 NORM2_SQ,
534 Err(LaError::non_finite_computation_step(
535 ArithmeticOperation::VectorSquaredNorm,
536 0,
537 ))
538 );
539 }
540
541 #[test]
542 fn vector_dot_and_norm2_sq_preserve_fma_and_left_to_right_order() {
543 let dot_large = 9_007_199_254_740_992.0;
544 let dot_lhs = Vector::<4>::new([dot_large, 1.0, 1.0, 1.0]);
545 let dot_rhs = Vector::<4>::new([1.0; 4]);
546 assert_eq!(dot_lhs.dot(&dot_rhs), Ok(dot_large));
547
548 let fused_lhs = Vector::<2>::new([f64::MAX, f64::MAX]);
549 let fused_rhs = Vector::<2>::new([-1.0, 2.0]);
550 assert_eq!(fused_lhs.dot(&fused_rhs), Ok(f64::MAX));
551
552 let norm_large = 134_217_728.0;
553 let vector = Vector::<4>::new([norm_large, 1.0, 1.0, 1.0]);
554 assert_eq!(vector.norm2_sq(), Ok(norm_large * norm_large));
555 }
556
557 #[test]
558 fn vector_dot_and_norm2_sq_report_first_middle_overflowing_step() {
559 let dot_lhs = Vector::<3>::new([f64::MAX, f64::MAX, 1.0]);
560 let dot_rhs = Vector::<3>::new([1.0; 3]);
561 assert_eq!(
562 dot_lhs.dot(&dot_rhs),
563 Err(LaError::non_finite_computation_step(
564 ArithmeticOperation::VectorDotProduct,
565 1,
566 ))
567 );
568
569 let norm_large = 1.0e154;
570 let vector = Vector::<3>::new([norm_large, norm_large, 1.0]);
571 assert_eq!(
572 vector.norm2_sq(),
573 Err(LaError::non_finite_computation_step(
574 ArithmeticOperation::VectorSquaredNorm,
575 1,
576 ))
577 );
578 }
579
580 #[test]
581 fn zero_dimension_vector_has_zero_dot_and_norm() {
582 let vector = Vector::<0>::try_new([]).unwrap();
583
584 assert!(vector.as_array().is_empty());
585 assert!(vector.into_array().is_empty());
586 assert_eq!(vector.dot(&Vector::zero()), Ok(0.0));
587 assert_eq!(vector.norm2_sq(), Ok(0.0));
588 }
589}