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