manifold_rust/
linalg.rs

1// linalg.rs — Linear algebra types for manifold-rust
2// Ports: include/manifold/linalg.h and the type aliases in include/manifold/common.h
3//
4// C++ naming:  mat<T, M, N> = M rows, N cols, column-major (N columns of vec<T,M>)
5// Rust types use the same naming as the C++ `using` aliases in common.h.
6//
7// All f64 types match C++ `double`; f32 appears only at the MeshGL boundary.
8
9use std::hash::{Hash, Hasher};
10use std::ops::{
11    Add, AddAssign, BitAnd, BitAndAssign, BitOr, BitOrAssign, BitXor, BitXorAssign, Div, DivAssign,
12    Index, IndexMut, Mul, MulAssign, Neg, Not, Rem, RemAssign, Shl, ShlAssign, Shr, ShrAssign,
13    Sub, SubAssign,
14};
15
16use crate::math;
17
18// ─── Vec2 ────────────────────────────────────────────────────────────────────
19
20/// `la::vec<double, 2>` — 2-component f64 column vector
21#[derive(Clone, Copy, Debug, Default, PartialEq)]
22#[repr(C)]
23pub struct Vec2 {
24    pub x: f64,
25    pub y: f64,
26}
27
28impl Vec2 {
29    #[inline]
30    pub const fn new(x: f64, y: f64) -> Self {
31        Self { x, y }
32    }
33    #[inline]
34    pub const fn splat(v: f64) -> Self {
35        Self { x: v, y: v }
36    }
37    #[inline]
38    pub fn xy(self) -> Vec2 {
39        self
40    }
41}
42
43impl Index<usize> for Vec2 {
44    type Output = f64;
45    fn index(&self, i: usize) -> &f64 {
46        match i {
47            0 => &self.x,
48            1 => &self.y,
49            _ => panic!("Vec2 index out of range: {i}"),
50        }
51    }
52}
53impl IndexMut<usize> for Vec2 {
54    fn index_mut(&mut self, i: usize) -> &mut f64 {
55        match i {
56            0 => &mut self.x,
57            1 => &mut self.y,
58            _ => panic!("Vec2 index out of range: {i}"),
59        }
60    }
61}
62
63// Lexicographic ordering — matches C++ compare(vec<T,2>, vec<T,2>)
64impl PartialOrd for Vec2 {
65    fn partial_cmp(&self, o: &Self) -> Option<std::cmp::Ordering> {
66        if self.x != o.x {
67            self.x.partial_cmp(&o.x)
68        } else {
69            self.y.partial_cmp(&o.y)
70        }
71    }
72}
73
74impl From<[f64; 2]> for Vec2 {
75    fn from(a: [f64; 2]) -> Self {
76        Self::new(a[0], a[1])
77    }
78}
79impl From<Vec2> for [f64; 2] {
80    fn from(v: Vec2) -> [f64; 2] {
81        [v.x, v.y]
82    }
83}
84impl From<(f64, f64)> for Vec2 {
85    fn from((x, y): (f64, f64)) -> Self {
86        Self::new(x, y)
87    }
88}
89
90macro_rules! impl_vec2_ops {
91    () => {
92        impl Neg for Vec2 {
93            type Output = Vec2;
94            fn neg(self) -> Vec2 {
95                Vec2::new(-self.x, -self.y)
96            }
97        }
98        impl Add for Vec2 {
99            type Output = Vec2;
100            fn add(self, b: Vec2) -> Vec2 {
101                Vec2::new(self.x + b.x, self.y + b.y)
102            }
103        }
104        impl Sub for Vec2 {
105            type Output = Vec2;
106            fn sub(self, b: Vec2) -> Vec2 {
107                Vec2::new(self.x - b.x, self.y - b.y)
108            }
109        }
110        // vec * vec is element-wise (cmul in C++)
111        impl Mul for Vec2 {
112            type Output = Vec2;
113            fn mul(self, b: Vec2) -> Vec2 {
114                Vec2::new(self.x * b.x, self.y * b.y)
115            }
116        }
117        impl Div for Vec2 {
118            type Output = Vec2;
119            fn div(self, b: Vec2) -> Vec2 {
120                Vec2::new(self.x / b.x, self.y / b.y)
121            }
122        }
123        impl Mul<f64> for Vec2 {
124            type Output = Vec2;
125            fn mul(self, s: f64) -> Vec2 {
126                Vec2::new(self.x * s, self.y * s)
127            }
128        }
129        impl Mul<Vec2> for f64 {
130            type Output = Vec2;
131            fn mul(self, v: Vec2) -> Vec2 {
132                Vec2::new(self * v.x, self * v.y)
133            }
134        }
135        impl Div<f64> for Vec2 {
136            type Output = Vec2;
137            fn div(self, s: f64) -> Vec2 {
138                Vec2::new(self.x / s, self.y / s)
139            }
140        }
141        impl Add<f64> for Vec2 {
142            type Output = Vec2;
143            fn add(self, s: f64) -> Vec2 {
144                Vec2::new(self.x + s, self.y + s)
145            }
146        }
147        impl Sub<f64> for Vec2 {
148            type Output = Vec2;
149            fn sub(self, s: f64) -> Vec2 {
150                Vec2::new(self.x - s, self.y - s)
151            }
152        }
153        impl AddAssign for Vec2 {
154            fn add_assign(&mut self, b: Vec2) {
155                *self = *self + b;
156            }
157        }
158        impl SubAssign for Vec2 {
159            fn sub_assign(&mut self, b: Vec2) {
160                *self = *self - b;
161            }
162        }
163        impl MulAssign for Vec2 {
164            fn mul_assign(&mut self, b: Vec2) {
165                *self = *self * b;
166            }
167        }
168        impl MulAssign<f64> for Vec2 {
169            fn mul_assign(&mut self, s: f64) {
170                *self = *self * s;
171            }
172        }
173        impl DivAssign for Vec2 {
174            fn div_assign(&mut self, b: Vec2) {
175                *self = *self / b;
176            }
177        }
178        impl DivAssign<f64> for Vec2 {
179            fn div_assign(&mut self, s: f64) {
180                *self = *self / s;
181            }
182        }
183    };
184}
185impl_vec2_ops!();
186
187// ─── Vec3 ────────────────────────────────────────────────────────────────────
188
189/// `la::vec<double, 3>`
190#[derive(Clone, Copy, Debug, Default, PartialEq)]
191#[repr(C)]
192pub struct Vec3 {
193    pub x: f64,
194    pub y: f64,
195    pub z: f64,
196}
197
198impl Vec3 {
199    #[inline]
200    pub const fn new(x: f64, y: f64, z: f64) -> Self {
201        Self { x, y, z }
202    }
203    #[inline]
204    pub const fn splat(v: f64) -> Self {
205        Self { x: v, y: v, z: v }
206    }
207    #[inline]
208    pub fn xy(self) -> Vec2 {
209        Vec2::new(self.x, self.y)
210    }
211    #[inline]
212    pub fn yz(self) -> Vec2 {
213        Vec2::new(self.y, self.z)
214    }
215}
216
217impl Index<usize> for Vec3 {
218    type Output = f64;
219    fn index(&self, i: usize) -> &f64 {
220        match i {
221            0 => &self.x,
222            1 => &self.y,
223            2 => &self.z,
224            _ => panic!("Vec3 index out of range: {i}"),
225        }
226    }
227}
228impl IndexMut<usize> for Vec3 {
229    fn index_mut(&mut self, i: usize) -> &mut f64 {
230        match i {
231            0 => &mut self.x,
232            1 => &mut self.y,
233            2 => &mut self.z,
234            _ => panic!("Vec3 index out of range: {i}"),
235        }
236    }
237}
238
239impl PartialOrd for Vec3 {
240    fn partial_cmp(&self, o: &Self) -> Option<std::cmp::Ordering> {
241        if self.x != o.x {
242            self.x.partial_cmp(&o.x)
243        } else if self.y != o.y {
244            self.y.partial_cmp(&o.y)
245        } else {
246            self.z.partial_cmp(&o.z)
247        }
248    }
249}
250
251impl From<[f64; 3]> for Vec3 {
252    fn from(a: [f64; 3]) -> Self {
253        Self::new(a[0], a[1], a[2])
254    }
255}
256impl From<Vec3> for [f64; 3] {
257    fn from(v: Vec3) -> [f64; 3] {
258        [v.x, v.y, v.z]
259    }
260}
261impl From<(f64, f64, f64)> for Vec3 {
262    fn from((x, y, z): (f64, f64, f64)) -> Self {
263        Self::new(x, y, z)
264    }
265}
266impl From<(Vec2, f64)> for Vec3 {
267    fn from((xy, z): (Vec2, f64)) -> Self {
268        Self::new(xy.x, xy.y, z)
269    }
270}
271
272macro_rules! impl_vec3_ops {
273    () => {
274        impl Neg for Vec3 {
275            type Output = Vec3;
276            fn neg(self) -> Vec3 {
277                Vec3::new(-self.x, -self.y, -self.z)
278            }
279        }
280        impl Add for Vec3 {
281            type Output = Vec3;
282            fn add(self, b: Vec3) -> Vec3 {
283                Vec3::new(self.x + b.x, self.y + b.y, self.z + b.z)
284            }
285        }
286        impl Sub for Vec3 {
287            type Output = Vec3;
288            fn sub(self, b: Vec3) -> Vec3 {
289                Vec3::new(self.x - b.x, self.y - b.y, self.z - b.z)
290            }
291        }
292        impl Mul for Vec3 {
293            type Output = Vec3;
294            fn mul(self, b: Vec3) -> Vec3 {
295                Vec3::new(self.x * b.x, self.y * b.y, self.z * b.z)
296            }
297        }
298        impl Div for Vec3 {
299            type Output = Vec3;
300            fn div(self, b: Vec3) -> Vec3 {
301                Vec3::new(self.x / b.x, self.y / b.y, self.z / b.z)
302            }
303        }
304        impl Mul<f64> for Vec3 {
305            type Output = Vec3;
306            fn mul(self, s: f64) -> Vec3 {
307                Vec3::new(self.x * s, self.y * s, self.z * s)
308            }
309        }
310        impl Mul<Vec3> for f64 {
311            type Output = Vec3;
312            fn mul(self, v: Vec3) -> Vec3 {
313                Vec3::new(self * v.x, self * v.y, self * v.z)
314            }
315        }
316        impl Div<f64> for Vec3 {
317            type Output = Vec3;
318            fn div(self, s: f64) -> Vec3 {
319                Vec3::new(self.x / s, self.y / s, self.z / s)
320            }
321        }
322        impl Add<f64> for Vec3 {
323            type Output = Vec3;
324            fn add(self, s: f64) -> Vec3 {
325                Vec3::new(self.x + s, self.y + s, self.z + s)
326            }
327        }
328        impl Sub<f64> for Vec3 {
329            type Output = Vec3;
330            fn sub(self, s: f64) -> Vec3 {
331                Vec3::new(self.x - s, self.y - s, self.z - s)
332            }
333        }
334        impl AddAssign for Vec3 {
335            fn add_assign(&mut self, b: Vec3) {
336                *self = *self + b;
337            }
338        }
339        impl SubAssign for Vec3 {
340            fn sub_assign(&mut self, b: Vec3) {
341                *self = *self - b;
342            }
343        }
344        impl MulAssign for Vec3 {
345            fn mul_assign(&mut self, b: Vec3) {
346                *self = *self * b;
347            }
348        }
349        impl MulAssign<f64> for Vec3 {
350            fn mul_assign(&mut self, s: f64) {
351                *self = *self * s;
352            }
353        }
354        impl DivAssign for Vec3 {
355            fn div_assign(&mut self, b: Vec3) {
356                *self = *self / b;
357            }
358        }
359        impl DivAssign<f64> for Vec3 {
360            fn div_assign(&mut self, s: f64) {
361                *self = *self / s;
362            }
363        }
364    };
365}
366impl_vec3_ops!();
367
368// ─── Vec4 ────────────────────────────────────────────────────────────────────
369
370/// `la::vec<double, 4>` / `quat`
371#[derive(Clone, Copy, Debug, Default, PartialEq)]
372#[repr(C)]
373pub struct Vec4 {
374    pub x: f64,
375    pub y: f64,
376    pub z: f64,
377    pub w: f64,
378}
379
380impl Vec4 {
381    #[inline]
382    pub const fn new(x: f64, y: f64, z: f64, w: f64) -> Self {
383        Self { x, y, z, w }
384    }
385    #[inline]
386    pub const fn splat(v: f64) -> Self {
387        Self {
388            x: v,
389            y: v,
390            z: v,
391            w: v,
392        }
393    }
394    #[inline]
395    pub fn xy(self) -> Vec2 {
396        Vec2::new(self.x, self.y)
397    }
398    #[inline]
399    pub fn xyz(self) -> Vec3 {
400        Vec3::new(self.x, self.y, self.z)
401    }
402}
403
404impl Index<usize> for Vec4 {
405    type Output = f64;
406    fn index(&self, i: usize) -> &f64 {
407        match i {
408            0 => &self.x,
409            1 => &self.y,
410            2 => &self.z,
411            3 => &self.w,
412            _ => panic!("Vec4 index out of range: {i}"),
413        }
414    }
415}
416impl IndexMut<usize> for Vec4 {
417    fn index_mut(&mut self, i: usize) -> &mut f64 {
418        match i {
419            0 => &mut self.x,
420            1 => &mut self.y,
421            2 => &mut self.z,
422            3 => &mut self.w,
423            _ => panic!("Vec4 index out of range: {i}"),
424        }
425    }
426}
427
428impl PartialOrd for Vec4 {
429    fn partial_cmp(&self, o: &Self) -> Option<std::cmp::Ordering> {
430        if self.x != o.x {
431            self.x.partial_cmp(&o.x)
432        } else if self.y != o.y {
433            self.y.partial_cmp(&o.y)
434        } else if self.z != o.z {
435            self.z.partial_cmp(&o.z)
436        } else {
437            self.w.partial_cmp(&o.w)
438        }
439    }
440}
441
442impl From<[f64; 4]> for Vec4 {
443    fn from(a: [f64; 4]) -> Self {
444        Self::new(a[0], a[1], a[2], a[3])
445    }
446}
447impl From<Vec4> for [f64; 4] {
448    fn from(v: Vec4) -> [f64; 4] {
449        [v.x, v.y, v.z, v.w]
450    }
451}
452impl From<(Vec3, f64)> for Vec4 {
453    fn from((xyz, w): (Vec3, f64)) -> Self {
454        Self::new(xyz.x, xyz.y, xyz.z, w)
455    }
456}
457impl From<(Vec2, f64, f64)> for Vec4 {
458    fn from((xy, z, w): (Vec2, f64, f64)) -> Self {
459        Self::new(xy.x, xy.y, z, w)
460    }
461}
462
463macro_rules! impl_vec4_ops {
464    () => {
465        impl Neg for Vec4 {
466            type Output = Vec4;
467            fn neg(self) -> Vec4 {
468                Vec4::new(-self.x, -self.y, -self.z, -self.w)
469            }
470        }
471        impl Add for Vec4 {
472            type Output = Vec4;
473            fn add(self, b: Vec4) -> Vec4 {
474                Vec4::new(self.x + b.x, self.y + b.y, self.z + b.z, self.w + b.w)
475            }
476        }
477        impl Sub for Vec4 {
478            type Output = Vec4;
479            fn sub(self, b: Vec4) -> Vec4 {
480                Vec4::new(self.x - b.x, self.y - b.y, self.z - b.z, self.w - b.w)
481            }
482        }
483        impl Mul for Vec4 {
484            type Output = Vec4;
485            fn mul(self, b: Vec4) -> Vec4 {
486                Vec4::new(self.x * b.x, self.y * b.y, self.z * b.z, self.w * b.w)
487            }
488        }
489        impl Div for Vec4 {
490            type Output = Vec4;
491            fn div(self, b: Vec4) -> Vec4 {
492                Vec4::new(self.x / b.x, self.y / b.y, self.z / b.z, self.w / b.w)
493            }
494        }
495        impl Mul<f64> for Vec4 {
496            type Output = Vec4;
497            fn mul(self, s: f64) -> Vec4 {
498                Vec4::new(self.x * s, self.y * s, self.z * s, self.w * s)
499            }
500        }
501        impl Mul<Vec4> for f64 {
502            type Output = Vec4;
503            fn mul(self, v: Vec4) -> Vec4 {
504                Vec4::new(self * v.x, self * v.y, self * v.z, self * v.w)
505            }
506        }
507        impl Div<f64> for Vec4 {
508            type Output = Vec4;
509            fn div(self, s: f64) -> Vec4 {
510                Vec4::new(self.x / s, self.y / s, self.z / s, self.w / s)
511            }
512        }
513        impl Add<f64> for Vec4 {
514            type Output = Vec4;
515            fn add(self, s: f64) -> Vec4 {
516                Vec4::new(self.x + s, self.y + s, self.z + s, self.w + s)
517            }
518        }
519        impl Sub<f64> for Vec4 {
520            type Output = Vec4;
521            fn sub(self, s: f64) -> Vec4 {
522                Vec4::new(self.x - s, self.y - s, self.z - s, self.w - s)
523            }
524        }
525        impl AddAssign for Vec4 {
526            fn add_assign(&mut self, b: Vec4) {
527                *self = *self + b;
528            }
529        }
530        impl SubAssign for Vec4 {
531            fn sub_assign(&mut self, b: Vec4) {
532                *self = *self - b;
533            }
534        }
535        impl MulAssign for Vec4 {
536            fn mul_assign(&mut self, b: Vec4) {
537                *self = *self * b;
538            }
539        }
540        impl MulAssign<f64> for Vec4 {
541            fn mul_assign(&mut self, s: f64) {
542                *self = *self * s;
543            }
544        }
545        impl DivAssign for Vec4 {
546            fn div_assign(&mut self, b: Vec4) {
547                *self = *self / b;
548            }
549        }
550        impl DivAssign<f64> for Vec4 {
551            fn div_assign(&mut self, s: f64) {
552                *self = *self / s;
553            }
554        }
555    };
556}
557impl_vec4_ops!();
558
559/// Quaternion — same bits as Vec4, used where the value represents xi+yj+zk+w
560pub type Quat = Vec4;
561
562// ─── IVec2 / IVec3 / IVec4 ───────────────────────────────────────────────────
563
564/// `la::vec<int, 2>`
565#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
566#[repr(C)]
567pub struct IVec2 {
568    pub x: i32,
569    pub y: i32,
570}
571
572impl IVec2 {
573    #[inline]
574    pub const fn new(x: i32, y: i32) -> Self {
575        Self { x, y }
576    }
577    #[inline]
578    pub const fn splat(v: i32) -> Self {
579        Self { x: v, y: v }
580    }
581}
582
583impl Index<usize> for IVec2 {
584    type Output = i32;
585    fn index(&self, i: usize) -> &i32 {
586        match i {
587            0 => &self.x,
588            1 => &self.y,
589            _ => panic!("IVec2 index out of range: {i}"),
590        }
591    }
592}
593impl IndexMut<usize> for IVec2 {
594    fn index_mut(&mut self, i: usize) -> &mut i32 {
595        match i {
596            0 => &mut self.x,
597            1 => &mut self.y,
598            _ => panic!("IVec2 index out of range: {i}"),
599        }
600    }
601}
602
603impl PartialOrd for IVec2 {
604    fn partial_cmp(&self, o: &Self) -> Option<std::cmp::Ordering> {
605        Some(self.cmp(o))
606    }
607}
608impl Ord for IVec2 {
609    fn cmp(&self, o: &Self) -> std::cmp::Ordering {
610        self.x.cmp(&o.x).then(self.y.cmp(&o.y))
611    }
612}
613
614impl Hash for IVec2 {
615    fn hash<H: Hasher>(&self, state: &mut H) {
616        self.x.hash(state);
617        self.y.hash(state);
618    }
619}
620
621impl From<[i32; 2]> for IVec2 {
622    fn from(a: [i32; 2]) -> Self {
623        Self::new(a[0], a[1])
624    }
625}
626
627impl Neg for IVec2 {
628    type Output = IVec2;
629    fn neg(self) -> IVec2 {
630        IVec2::new(-self.x, -self.y)
631    }
632}
633impl Add for IVec2 {
634    type Output = IVec2;
635    fn add(self, b: IVec2) -> IVec2 {
636        IVec2::new(self.x + b.x, self.y + b.y)
637    }
638}
639impl Sub for IVec2 {
640    type Output = IVec2;
641    fn sub(self, b: IVec2) -> IVec2 {
642        IVec2::new(self.x - b.x, self.y - b.y)
643    }
644}
645impl Mul for IVec2 {
646    type Output = IVec2;
647    fn mul(self, b: IVec2) -> IVec2 {
648        IVec2::new(self.x * b.x, self.y * b.y)
649    }
650}
651impl Mul<i32> for IVec2 {
652    type Output = IVec2;
653    fn mul(self, s: i32) -> IVec2 {
654        IVec2::new(self.x * s, self.y * s)
655    }
656}
657impl Mul<IVec2> for i32 {
658    type Output = IVec2;
659    fn mul(self, v: IVec2) -> IVec2 {
660        IVec2::new(self * v.x, self * v.y)
661    }
662}
663impl AddAssign for IVec2 {
664    fn add_assign(&mut self, b: IVec2) {
665        *self = *self + b;
666    }
667}
668
669/// `la::vec<int, 3>`
670#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
671#[repr(C)]
672pub struct IVec3 {
673    pub x: i32,
674    pub y: i32,
675    pub z: i32,
676}
677
678impl IVec3 {
679    #[inline]
680    pub const fn new(x: i32, y: i32, z: i32) -> Self {
681        Self { x, y, z }
682    }
683    #[inline]
684    pub const fn splat(v: i32) -> Self {
685        Self { x: v, y: v, z: v }
686    }
687    #[inline]
688    pub fn xy(self) -> IVec2 {
689        IVec2::new(self.x, self.y)
690    }
691}
692
693impl Index<usize> for IVec3 {
694    type Output = i32;
695    fn index(&self, i: usize) -> &i32 {
696        match i {
697            0 => &self.x,
698            1 => &self.y,
699            2 => &self.z,
700            _ => panic!("IVec3 index out of range: {i}"),
701        }
702    }
703}
704impl IndexMut<usize> for IVec3 {
705    fn index_mut(&mut self, i: usize) -> &mut i32 {
706        match i {
707            0 => &mut self.x,
708            1 => &mut self.y,
709            2 => &mut self.z,
710            _ => panic!("IVec3 index out of range: {i}"),
711        }
712    }
713}
714
715impl PartialOrd for IVec3 {
716    fn partial_cmp(&self, o: &Self) -> Option<std::cmp::Ordering> {
717        Some(self.cmp(o))
718    }
719}
720impl Ord for IVec3 {
721    fn cmp(&self, o: &Self) -> std::cmp::Ordering {
722        self.x.cmp(&o.x).then(self.y.cmp(&o.y)).then(self.z.cmp(&o.z))
723    }
724}
725
726impl Hash for IVec3 {
727    fn hash<H: Hasher>(&self, state: &mut H) {
728        self.x.hash(state);
729        self.y.hash(state);
730        self.z.hash(state);
731    }
732}
733
734impl From<[i32; 3]> for IVec3 {
735    fn from(a: [i32; 3]) -> Self {
736        Self::new(a[0], a[1], a[2])
737    }
738}
739impl From<IVec3> for [i32; 3] {
740    fn from(v: IVec3) -> [i32; 3] {
741        [v.x, v.y, v.z]
742    }
743}
744
745impl Neg for IVec3 {
746    type Output = IVec3;
747    fn neg(self) -> IVec3 {
748        IVec3::new(-self.x, -self.y, -self.z)
749    }
750}
751impl Add for IVec3 {
752    type Output = IVec3;
753    fn add(self, b: IVec3) -> IVec3 {
754        IVec3::new(self.x + b.x, self.y + b.y, self.z + b.z)
755    }
756}
757impl Sub for IVec3 {
758    type Output = IVec3;
759    fn sub(self, b: IVec3) -> IVec3 {
760        IVec3::new(self.x - b.x, self.y - b.y, self.z - b.z)
761    }
762}
763impl Mul for IVec3 {
764    type Output = IVec3;
765    fn mul(self, b: IVec3) -> IVec3 {
766        IVec3::new(self.x * b.x, self.y * b.y, self.z * b.z)
767    }
768}
769impl Mul<i32> for IVec3 {
770    type Output = IVec3;
771    fn mul(self, s: i32) -> IVec3 {
772        IVec3::new(self.x * s, self.y * s, self.z * s)
773    }
774}
775impl AddAssign for IVec3 {
776    fn add_assign(&mut self, b: IVec3) {
777        *self = *self + b;
778    }
779}
780
781/// `la::vec<int, 4>`
782#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
783#[repr(C)]
784pub struct IVec4 {
785    pub x: i32,
786    pub y: i32,
787    pub z: i32,
788    pub w: i32,
789}
790
791impl IVec4 {
792    #[inline]
793    pub const fn new(x: i32, y: i32, z: i32, w: i32) -> Self {
794        Self { x, y, z, w }
795    }
796    #[inline]
797    pub const fn splat(v: i32) -> Self {
798        Self {
799            x: v,
800            y: v,
801            z: v,
802            w: v,
803        }
804    }
805    #[inline]
806    pub fn xyz(self) -> IVec3 {
807        IVec3::new(self.x, self.y, self.z)
808    }
809}
810
811impl Index<usize> for IVec4 {
812    type Output = i32;
813    fn index(&self, i: usize) -> &i32 {
814        match i {
815            0 => &self.x,
816            1 => &self.y,
817            2 => &self.z,
818            3 => &self.w,
819            _ => panic!("IVec4 index out of range: {i}"),
820        }
821    }
822}
823impl IndexMut<usize> for IVec4 {
824    fn index_mut(&mut self, i: usize) -> &mut i32 {
825        match i {
826            0 => &mut self.x,
827            1 => &mut self.y,
828            2 => &mut self.z,
829            3 => &mut self.w,
830            _ => panic!("IVec4 index out of range: {i}"),
831        }
832    }
833}
834
835impl PartialOrd for IVec4 {
836    fn partial_cmp(&self, o: &Self) -> Option<std::cmp::Ordering> {
837        Some(self.cmp(o))
838    }
839}
840impl Ord for IVec4 {
841    fn cmp(&self, o: &Self) -> std::cmp::Ordering {
842        self.x
843            .cmp(&o.x)
844            .then(self.y.cmp(&o.y))
845            .then(self.z.cmp(&o.z))
846            .then(self.w.cmp(&o.w))
847    }
848}
849
850impl Hash for IVec4 {
851    fn hash<H: Hasher>(&self, state: &mut H) {
852        self.x.hash(state);
853        self.y.hash(state);
854        self.z.hash(state);
855        self.w.hash(state);
856    }
857}
858
859impl From<[i32; 4]> for IVec4 {
860    fn from(a: [i32; 4]) -> Self {
861        Self::new(a[0], a[1], a[2], a[3])
862    }
863}
864
865impl Neg for IVec4 {
866    type Output = IVec4;
867    fn neg(self) -> IVec4 {
868        IVec4::new(-self.x, -self.y, -self.z, -self.w)
869    }
870}
871impl Add for IVec4 {
872    type Output = IVec4;
873    fn add(self, b: IVec4) -> IVec4 {
874        IVec4::new(self.x + b.x, self.y + b.y, self.z + b.z, self.w + b.w)
875    }
876}
877impl Sub for IVec4 {
878    type Output = IVec4;
879    fn sub(self, b: IVec4) -> IVec4 {
880        IVec4::new(self.x - b.x, self.y - b.y, self.z - b.z, self.w - b.w)
881    }
882}
883impl Mul<i32> for IVec4 {
884    type Output = IVec4;
885    fn mul(self, s: i32) -> IVec4 {
886        IVec4::new(self.x * s, self.y * s, self.z * s, self.w * s)
887    }
888}
889
890// ─── BVec4 ───────────────────────────────────────────────────────────────────
891
892/// `la::vec<bool, 4>`
893#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
894pub struct BVec4 {
895    pub x: bool,
896    pub y: bool,
897    pub z: bool,
898    pub w: bool,
899}
900
901impl BVec4 {
902    #[inline]
903    pub const fn new(x: bool, y: bool, z: bool, w: bool) -> Self {
904        Self { x, y, z, w }
905    }
906    #[inline]
907    pub const fn splat(v: bool) -> Self {
908        Self {
909            x: v,
910            y: v,
911            z: v,
912            w: v,
913        }
914    }
915}
916
917impl Index<usize> for BVec4 {
918    type Output = bool;
919    fn index(&self, i: usize) -> &bool {
920        match i {
921            0 => &self.x,
922            1 => &self.y,
923            2 => &self.z,
924            3 => &self.w,
925            _ => panic!("BVec4 index out of range: {i}"),
926        }
927    }
928}
929impl IndexMut<usize> for BVec4 {
930    fn index_mut(&mut self, i: usize) -> &mut bool {
931        match i {
932            0 => &mut self.x,
933            1 => &mut self.y,
934            2 => &mut self.z,
935            3 => &mut self.w,
936            _ => panic!("BVec4 index out of range: {i}"),
937        }
938    }
939}
940
941impl Not for BVec4 {
942    type Output = BVec4;
943    fn not(self) -> BVec4 {
944        BVec4::new(!self.x, !self.y, !self.z, !self.w)
945    }
946}
947impl BitAnd for BVec4 {
948    type Output = BVec4;
949    fn bitand(self, b: BVec4) -> BVec4 {
950        BVec4::new(self.x & b.x, self.y & b.y, self.z & b.z, self.w & b.w)
951    }
952}
953impl BitOr for BVec4 {
954    type Output = BVec4;
955    fn bitor(self, b: BVec4) -> BVec4 {
956        BVec4::new(self.x | b.x, self.y | b.y, self.z | b.z, self.w | b.w)
957    }
958}
959
960// ─── UVec3 ───────────────────────────────────────────────────────────────────
961
962/// u32 index vector — used for triangle index triples in Morton sorting
963#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
964#[repr(C)]
965pub struct UVec3 {
966    pub x: u32,
967    pub y: u32,
968    pub z: u32,
969}
970
971impl UVec3 {
972    #[inline]
973    pub const fn new(x: u32, y: u32, z: u32) -> Self {
974        Self { x, y, z }
975    }
976}
977
978impl Index<usize> for UVec3 {
979    type Output = u32;
980    fn index(&self, i: usize) -> &u32 {
981        match i {
982            0 => &self.x,
983            1 => &self.y,
984            2 => &self.z,
985            _ => panic!("UVec3 index out of range: {i}"),
986        }
987    }
988}
989
990// ─── Matrix types ─────────────────────────────────────────────────────────────
991// mat<T, M, N>: M rows, N cols, stored as N column vectors each of length M.
992// Named following C++ common.h aliases (mat3x4 = 3 rows, 4 cols).
993
994/// `la::mat<double, 2, 2>` — 2×2 matrix stored as 2 columns of Vec2
995#[derive(Clone, Copy, Debug, Default, PartialEq)]
996pub struct Mat2 {
997    pub x: Vec2, // column 0
998    pub y: Vec2, // column 1
999}
1000
1001impl Mat2 {
1002    #[inline]
1003    pub const fn from_cols(x: Vec2, y: Vec2) -> Self {
1004        Self { x, y }
1005    }
1006    /// Identity matrix
1007    pub fn identity() -> Self {
1008        Self {
1009            x: Vec2::new(1.0, 0.0),
1010            y: Vec2::new(0.0, 1.0),
1011        }
1012    }
1013    pub fn row(&self, i: usize) -> Vec2 {
1014        Vec2::new(self.x[i], self.y[i])
1015    }
1016    pub fn transpose(&self) -> Self {
1017        Self {
1018            x: self.row(0),
1019            y: self.row(1),
1020        }
1021    }
1022    pub fn determinant(&self) -> f64 {
1023        self.x.x * self.y.y - self.x.y * self.y.x
1024    }
1025}
1026
1027impl Index<usize> for Mat2 {
1028    type Output = Vec2;
1029    fn index(&self, j: usize) -> &Vec2 {
1030        match j {
1031            0 => &self.x,
1032            1 => &self.y,
1033            _ => panic!("Mat2 column index out of range: {j}"),
1034        }
1035    }
1036}
1037impl IndexMut<usize> for Mat2 {
1038    fn index_mut(&mut self, j: usize) -> &mut Vec2 {
1039        match j {
1040            0 => &mut self.x,
1041            1 => &mut self.y,
1042            _ => panic!("Mat2 column index out of range: {j}"),
1043        }
1044    }
1045}
1046
1047impl Mul<Vec2> for Mat2 {
1048    type Output = Vec2;
1049    fn mul(self, b: Vec2) -> Vec2 {
1050        self.x * b.x + self.y * b.y
1051    }
1052}
1053impl Mul for Mat2 {
1054    type Output = Mat2;
1055    fn mul(self, b: Mat2) -> Mat2 {
1056        Mat2::from_cols(self * b.x, self * b.y)
1057    }
1058}
1059impl Mul<f64> for Mat2 {
1060    type Output = Mat2;
1061    fn mul(self, s: f64) -> Mat2 {
1062        Mat2::from_cols(self.x * s, self.y * s)
1063    }
1064}
1065impl Add for Mat2 {
1066    type Output = Mat2;
1067    fn add(self, b: Mat2) -> Mat2 {
1068        Mat2::from_cols(self.x + b.x, self.y + b.y)
1069    }
1070}
1071impl Sub for Mat2 {
1072    type Output = Mat2;
1073    fn sub(self, b: Mat2) -> Mat2 {
1074        Mat2::from_cols(self.x - b.x, self.y - b.y)
1075    }
1076}
1077impl Neg for Mat2 {
1078    type Output = Mat2;
1079    fn neg(self) -> Mat2 {
1080        Mat2::from_cols(-self.x, -self.y)
1081    }
1082}
1083
1084/// `la::mat<double, 3, 3>` — 3×3 matrix
1085#[derive(Clone, Copy, Debug, Default, PartialEq)]
1086pub struct Mat3 {
1087    pub x: Vec3, // column 0
1088    pub y: Vec3, // column 1
1089    pub z: Vec3, // column 2
1090}
1091
1092impl Mat3 {
1093    #[inline]
1094    pub const fn from_cols(x: Vec3, y: Vec3, z: Vec3) -> Self {
1095        Self { x, y, z }
1096    }
1097    pub fn identity() -> Self {
1098        Self {
1099            x: Vec3::new(1.0, 0.0, 0.0),
1100            y: Vec3::new(0.0, 1.0, 0.0),
1101            z: Vec3::new(0.0, 0.0, 1.0),
1102        }
1103    }
1104    pub fn row(&self, i: usize) -> Vec3 {
1105        Vec3::new(self.x[i], self.y[i], self.z[i])
1106    }
1107    pub fn transpose(&self) -> Self {
1108        Self {
1109            x: self.row(0),
1110            y: self.row(1),
1111            z: self.row(2),
1112        }
1113    }
1114    /// C++ adjugate(mat<T,3,3>)
1115    pub fn adjugate(&self) -> Self {
1116        let a = self;
1117        Self {
1118            x: Vec3::new(
1119                a.y.y * a.z.z - a.z.y * a.y.z,
1120                a.z.y * a.x.z - a.x.y * a.z.z,
1121                a.x.y * a.y.z - a.y.y * a.x.z,
1122            ),
1123            y: Vec3::new(
1124                a.y.z * a.z.x - a.z.z * a.y.x,
1125                a.z.z * a.x.x - a.x.z * a.z.x,
1126                a.x.z * a.y.x - a.y.z * a.x.x,
1127            ),
1128            z: Vec3::new(
1129                a.y.x * a.z.y - a.z.x * a.y.y,
1130                a.z.x * a.x.y - a.x.x * a.z.y,
1131                a.x.x * a.y.y - a.y.x * a.x.y,
1132            ),
1133        }
1134    }
1135    pub fn determinant(&self) -> f64 {
1136        let a = self;
1137        a.x.x * (a.y.y * a.z.z - a.z.y * a.y.z)
1138            + a.x.y * (a.y.z * a.z.x - a.z.z * a.y.x)
1139            + a.x.z * (a.y.x * a.z.y - a.z.x * a.y.y)
1140    }
1141    pub fn inverse(&self) -> Self {
1142        self.adjugate() * (1.0 / self.determinant())
1143    }
1144    pub fn diagonal(&self) -> Vec3 {
1145        Vec3::new(self.x.x, self.y.y, self.z.z)
1146    }
1147    pub fn trace(&self) -> f64 {
1148        self.x.x + self.y.y + self.z.z
1149    }
1150}
1151
1152impl Index<usize> for Mat3 {
1153    type Output = Vec3;
1154    fn index(&self, j: usize) -> &Vec3 {
1155        match j {
1156            0 => &self.x,
1157            1 => &self.y,
1158            2 => &self.z,
1159            _ => panic!("Mat3 column index out of range: {j}"),
1160        }
1161    }
1162}
1163impl IndexMut<usize> for Mat3 {
1164    fn index_mut(&mut self, j: usize) -> &mut Vec3 {
1165        match j {
1166            0 => &mut self.x,
1167            1 => &mut self.y,
1168            2 => &mut self.z,
1169            _ => panic!("Mat3 column index out of range: {j}"),
1170        }
1171    }
1172}
1173
1174impl Mul<Vec3> for Mat3 {
1175    type Output = Vec3;
1176    fn mul(self, b: Vec3) -> Vec3 {
1177        self.x * b.x + self.y * b.y + self.z * b.z
1178    }
1179}
1180impl Mul for Mat3 {
1181    type Output = Mat3;
1182    fn mul(self, b: Mat3) -> Mat3 {
1183        Mat3::from_cols(self * b.x, self * b.y, self * b.z)
1184    }
1185}
1186impl Mul<f64> for Mat3 {
1187    type Output = Mat3;
1188    fn mul(self, s: f64) -> Mat3 {
1189        Mat3::from_cols(self.x * s, self.y * s, self.z * s)
1190    }
1191}
1192impl Add for Mat3 {
1193    type Output = Mat3;
1194    fn add(self, b: Mat3) -> Mat3 {
1195        Mat3::from_cols(self.x + b.x, self.y + b.y, self.z + b.z)
1196    }
1197}
1198impl Sub for Mat3 {
1199    type Output = Mat3;
1200    fn sub(self, b: Mat3) -> Mat3 {
1201        Mat3::from_cols(self.x - b.x, self.y - b.y, self.z - b.z)
1202    }
1203}
1204impl Neg for Mat3 {
1205    type Output = Mat3;
1206    fn neg(self) -> Mat3 {
1207        Mat3::from_cols(-self.x, -self.y, -self.z)
1208    }
1209}
1210
1211/// `la::mat<double, 4, 4>` — 4×4 matrix
1212#[derive(Clone, Copy, Debug, Default, PartialEq)]
1213pub struct Mat4 {
1214    pub x: Vec4, // column 0
1215    pub y: Vec4, // column 1
1216    pub z: Vec4, // column 2
1217    pub w: Vec4, // column 3
1218}
1219
1220impl Mat4 {
1221    #[inline]
1222    pub const fn from_cols(x: Vec4, y: Vec4, z: Vec4, w: Vec4) -> Self {
1223        Self { x, y, z, w }
1224    }
1225    pub fn identity() -> Self {
1226        Self {
1227            x: Vec4::new(1.0, 0.0, 0.0, 0.0),
1228            y: Vec4::new(0.0, 1.0, 0.0, 0.0),
1229            z: Vec4::new(0.0, 0.0, 1.0, 0.0),
1230            w: Vec4::new(0.0, 0.0, 0.0, 1.0),
1231        }
1232    }
1233    pub fn row(&self, i: usize) -> Vec4 {
1234        Vec4::new(self.x[i], self.y[i], self.z[i], self.w[i])
1235    }
1236    pub fn transpose(&self) -> Self {
1237        Self {
1238            x: self.row(0),
1239            y: self.row(1),
1240            z: self.row(2),
1241            w: self.row(3),
1242        }
1243    }
1244    pub fn determinant(&self) -> f64 {
1245        let a = self;
1246        a.x.x
1247            * (a.y.y * a.z.z * a.w.w
1248                + a.w.y * a.y.z * a.z.w
1249                + a.z.y * a.w.z * a.y.w
1250                - a.y.y * a.w.z * a.z.w
1251                - a.z.y * a.y.z * a.w.w
1252                - a.w.y * a.z.z * a.y.w)
1253            + a.x.y
1254                * (a.y.z * a.w.w * a.z.x
1255                    + a.z.z * a.y.w * a.w.x
1256                    + a.w.z * a.z.w * a.y.x
1257                    - a.y.z * a.z.w * a.w.x
1258                    - a.w.z * a.y.w * a.z.x
1259                    - a.z.z * a.w.w * a.y.x)
1260            + a.x.z
1261                * (a.y.w * a.z.x * a.w.y
1262                    + a.w.w * a.y.x * a.z.y
1263                    + a.z.w * a.w.x * a.y.y
1264                    - a.y.w * a.w.x * a.z.y
1265                    - a.z.w * a.y.x * a.w.y
1266                    - a.w.w * a.z.x * a.y.y)
1267            + a.x.w
1268                * (a.y.x * a.w.y * a.z.z
1269                    + a.z.x * a.y.y * a.w.z
1270                    + a.w.x * a.z.y * a.y.z
1271                    - a.y.x * a.z.y * a.w.z
1272                    - a.w.x * a.y.y * a.z.z
1273                    - a.z.x * a.w.y * a.y.z)
1274    }
1275    pub fn adjugate(&self) -> Self {
1276        let a = self;
1277        Self {
1278            x: Vec4::new(
1279                a.y.y * a.z.z * a.w.w
1280                    + a.w.y * a.y.z * a.z.w
1281                    + a.z.y * a.w.z * a.y.w
1282                    - a.y.y * a.w.z * a.z.w
1283                    - a.z.y * a.y.z * a.w.w
1284                    - a.w.y * a.z.z * a.y.w,
1285                a.x.y * a.w.z * a.z.w
1286                    + a.z.y * a.x.z * a.w.w
1287                    + a.w.y * a.z.z * a.x.w
1288                    - a.w.y * a.x.z * a.z.w
1289                    - a.z.y * a.w.z * a.x.w
1290                    - a.x.y * a.z.z * a.w.w,
1291                a.x.y * a.y.z * a.w.w
1292                    + a.w.y * a.x.z * a.y.w
1293                    + a.y.y * a.w.z * a.x.w
1294                    - a.x.y * a.w.z * a.y.w
1295                    - a.y.y * a.x.z * a.w.w
1296                    - a.w.y * a.y.z * a.x.w,
1297                a.x.y * a.z.z * a.y.w
1298                    + a.y.y * a.x.z * a.z.w
1299                    + a.z.y * a.y.z * a.x.w
1300                    - a.x.y * a.y.z * a.z.w
1301                    - a.z.y * a.x.z * a.y.w
1302                    - a.y.y * a.z.z * a.x.w,
1303            ),
1304            y: Vec4::new(
1305                a.y.z * a.w.w * a.z.x
1306                    + a.z.z * a.y.w * a.w.x
1307                    + a.w.z * a.z.w * a.y.x
1308                    - a.y.z * a.z.w * a.w.x
1309                    - a.w.z * a.y.w * a.z.x
1310                    - a.z.z * a.w.w * a.y.x,
1311                a.x.z * a.z.w * a.w.x
1312                    + a.w.z * a.x.w * a.z.x
1313                    + a.z.z * a.w.w * a.x.x
1314                    - a.x.z * a.w.w * a.z.x
1315                    - a.z.z * a.x.w * a.w.x
1316                    - a.w.z * a.z.w * a.x.x,
1317                a.x.z * a.w.w * a.y.x
1318                    + a.y.z * a.x.w * a.w.x
1319                    + a.w.z * a.y.w * a.x.x
1320                    - a.x.z * a.y.w * a.w.x
1321                    - a.w.z * a.x.w * a.y.x
1322                    - a.y.z * a.w.w * a.x.x,
1323                a.x.z * a.y.w * a.z.x
1324                    + a.z.z * a.x.w * a.y.x
1325                    + a.y.z * a.z.w * a.x.x
1326                    - a.x.z * a.z.w * a.y.x
1327                    - a.y.z * a.x.w * a.z.x
1328                    - a.z.z * a.y.w * a.x.x,
1329            ),
1330            z: Vec4::new(
1331                a.y.w * a.z.x * a.w.y
1332                    + a.w.w * a.y.x * a.z.y
1333                    + a.z.w * a.w.x * a.y.y
1334                    - a.y.w * a.w.x * a.z.y
1335                    - a.z.w * a.y.x * a.w.y
1336                    - a.w.w * a.z.x * a.y.y,
1337                a.x.w * a.w.x * a.z.y
1338                    + a.z.w * a.x.x * a.w.y
1339                    + a.w.w * a.z.x * a.x.y
1340                    - a.x.w * a.z.x * a.w.y
1341                    - a.w.w * a.x.x * a.z.y
1342                    - a.z.w * a.w.x * a.x.y,
1343                a.x.w * a.y.x * a.w.y
1344                    + a.w.w * a.x.x * a.y.y
1345                    + a.y.w * a.w.x * a.x.y
1346                    - a.x.w * a.w.x * a.y.y
1347                    - a.y.w * a.x.x * a.w.y
1348                    - a.w.w * a.y.x * a.x.y,
1349                a.x.w * a.z.x * a.y.y
1350                    + a.y.w * a.x.x * a.z.y
1351                    + a.z.w * a.y.x * a.x.y
1352                    - a.x.w * a.y.x * a.z.y
1353                    - a.z.w * a.x.x * a.y.y
1354                    - a.y.w * a.z.x * a.x.y,
1355            ),
1356            w: Vec4::new(
1357                a.y.x * a.w.y * a.z.z
1358                    + a.z.x * a.y.y * a.w.z
1359                    + a.w.x * a.z.y * a.y.z
1360                    - a.y.x * a.z.y * a.w.z
1361                    - a.w.x * a.y.y * a.z.z
1362                    - a.z.x * a.w.y * a.y.z,
1363                a.x.x * a.z.y * a.w.z
1364                    + a.w.x * a.x.y * a.z.z
1365                    + a.z.x * a.w.y * a.x.z
1366                    - a.x.x * a.w.y * a.z.z
1367                    - a.z.x * a.x.y * a.w.z
1368                    - a.w.x * a.z.y * a.x.z,
1369                a.x.x * a.w.y * a.y.z
1370                    + a.y.x * a.x.y * a.w.z
1371                    + a.w.x * a.y.y * a.x.z
1372                    - a.x.x * a.y.y * a.w.z
1373                    - a.w.x * a.x.y * a.y.z
1374                    - a.y.x * a.w.y * a.x.z,
1375                a.x.x * a.y.y * a.z.z
1376                    + a.z.x * a.x.y * a.y.z
1377                    + a.y.x * a.z.y * a.x.z
1378                    - a.x.x * a.z.y * a.y.z
1379                    - a.y.x * a.x.y * a.z.z
1380                    - a.z.x * a.y.y * a.x.z,
1381            ),
1382        }
1383    }
1384    pub fn inverse(&self) -> Self {
1385        self.adjugate() * (1.0 / self.determinant())
1386    }
1387}
1388
1389impl Index<usize> for Mat4 {
1390    type Output = Vec4;
1391    fn index(&self, j: usize) -> &Vec4 {
1392        match j {
1393            0 => &self.x,
1394            1 => &self.y,
1395            2 => &self.z,
1396            3 => &self.w,
1397            _ => panic!("Mat4 column index out of range: {j}"),
1398        }
1399    }
1400}
1401impl IndexMut<usize> for Mat4 {
1402    fn index_mut(&mut self, j: usize) -> &mut Vec4 {
1403        match j {
1404            0 => &mut self.x,
1405            1 => &mut self.y,
1406            2 => &mut self.z,
1407            3 => &mut self.w,
1408            _ => panic!("Mat4 column index out of range: {j}"),
1409        }
1410    }
1411}
1412
1413impl Mul<Vec4> for Mat4 {
1414    type Output = Vec4;
1415    fn mul(self, b: Vec4) -> Vec4 {
1416        self.x * b.x + self.y * b.y + self.z * b.z + self.w * b.w
1417    }
1418}
1419impl Mul for Mat4 {
1420    type Output = Mat4;
1421    fn mul(self, b: Mat4) -> Mat4 {
1422        Mat4::from_cols(self * b.x, self * b.y, self * b.z, self * b.w)
1423    }
1424}
1425impl Mul<f64> for Mat4 {
1426    type Output = Mat4;
1427    fn mul(self, s: f64) -> Mat4 {
1428        Mat4::from_cols(self.x * s, self.y * s, self.z * s, self.w * s)
1429    }
1430}
1431impl Add for Mat4 {
1432    type Output = Mat4;
1433    fn add(self, b: Mat4) -> Mat4 {
1434        Mat4::from_cols(self.x + b.x, self.y + b.y, self.z + b.z, self.w + b.w)
1435    }
1436}
1437impl Sub for Mat4 {
1438    type Output = Mat4;
1439    fn sub(self, b: Mat4) -> Mat4 {
1440        Mat4::from_cols(self.x - b.x, self.y - b.y, self.z - b.z, self.w - b.w)
1441    }
1442}
1443impl Neg for Mat4 {
1444    type Output = Mat4;
1445    fn neg(self) -> Mat4 {
1446        Mat4::from_cols(-self.x, -self.y, -self.z, -self.w)
1447    }
1448}
1449
1450/// `la::mat<double, 3, 4>` — 3 rows, 4 columns (affine transform type)
1451/// Stored as 4 columns each of Vec3: [x|y|z|translation]
1452#[derive(Clone, Copy, Debug, Default, PartialEq)]
1453pub struct Mat3x4 {
1454    pub x: Vec3, // column 0
1455    pub y: Vec3, // column 1
1456    pub z: Vec3, // column 2
1457    pub w: Vec3, // column 3 (translation)
1458}
1459
1460impl Mat3x4 {
1461    #[inline]
1462    pub const fn from_cols(x: Vec3, y: Vec3, z: Vec3, w: Vec3) -> Self {
1463        Self { x, y, z, w }
1464    }
1465    /// Identity: rotation=I, translation=0
1466    pub fn identity() -> Self {
1467        Self {
1468            x: Vec3::new(1.0, 0.0, 0.0),
1469            y: Vec3::new(0.0, 1.0, 0.0),
1470            z: Vec3::new(0.0, 0.0, 1.0),
1471            w: Vec3::new(0.0, 0.0, 0.0),
1472        }
1473    }
1474    pub fn row(&self, i: usize) -> Vec4 {
1475        Vec4::new(self.x[i], self.y[i], self.z[i], self.w[i])
1476    }
1477    /// Extract the 3×3 rotation/scale submatrix (first 3 columns)
1478    pub fn rotation(&self) -> Mat3 {
1479        Mat3::from_cols(self.x, self.y, self.z)
1480    }
1481    /// Extract the translation column
1482    pub fn translation(&self) -> Vec3 {
1483        self.w
1484    }
1485}
1486
1487impl Index<usize> for Mat3x4 {
1488    type Output = Vec3;
1489    fn index(&self, j: usize) -> &Vec3 {
1490        match j {
1491            0 => &self.x,
1492            1 => &self.y,
1493            2 => &self.z,
1494            3 => &self.w,
1495            _ => panic!("Mat3x4 column index out of range: {j}"),
1496        }
1497    }
1498}
1499impl IndexMut<usize> for Mat3x4 {
1500    fn index_mut(&mut self, j: usize) -> &mut Vec3 {
1501        match j {
1502            0 => &mut self.x,
1503            1 => &mut self.y,
1504            2 => &mut self.z,
1505            3 => &mut self.w,
1506            _ => panic!("Mat3x4 column index out of range: {j}"),
1507        }
1508    }
1509}
1510
1511/// mat3x4 * vec4 → vec3  (matrix-vector multiplication)
1512impl Mul<Vec4> for Mat3x4 {
1513    type Output = Vec3;
1514    fn mul(self, b: Vec4) -> Vec3 {
1515        self.x * b.x + self.y * b.y + self.z * b.z + self.w * b.w
1516    }
1517}
1518
1519/// mat3x4 * mat4x4 → mat3x4  (chain transforms)
1520impl Mul<Mat4> for Mat3x4 {
1521    type Output = Mat3x4;
1522    fn mul(self, b: Mat4) -> Mat3x4 {
1523        Mat3x4::from_cols(self * b.x, self * b.y, self * b.z, self * b.w)
1524    }
1525}
1526
1527/// mat4x4 * mat3x4 — promotes mat3x4 to mat4x4 then multiplies; returns Mat4
1528impl Mul<Mat3x4> for Mat4 {
1529    type Output = Mat4;
1530    fn mul(self, b: Mat3x4) -> Mat4 {
1531        // Treat mat3x4 as mat4x4 with bottom row [0,0,0,1]
1532        let b4 = mat3x4_to_mat4(b);
1533        self * b4
1534    }
1535}
1536
1537impl Mul<f64> for Mat3x4 {
1538    type Output = Mat3x4;
1539    fn mul(self, s: f64) -> Mat3x4 {
1540        Mat3x4::from_cols(self.x * s, self.y * s, self.z * s, self.w * s)
1541    }
1542}
1543impl Add for Mat3x4 {
1544    type Output = Mat3x4;
1545    fn add(self, b: Mat3x4) -> Mat3x4 {
1546        Mat3x4::from_cols(self.x + b.x, self.y + b.y, self.z + b.z, self.w + b.w)
1547    }
1548}
1549impl Sub for Mat3x4 {
1550    type Output = Mat3x4;
1551    fn sub(self, b: Mat3x4) -> Mat3x4 {
1552        Mat3x4::from_cols(self.x - b.x, self.y - b.y, self.z - b.z, self.w - b.w)
1553    }
1554}
1555impl Neg for Mat3x4 {
1556    type Output = Mat3x4;
1557    fn neg(self) -> Mat3x4 {
1558        Mat3x4::from_cols(-self.x, -self.y, -self.z, -self.w)
1559    }
1560}
1561
1562/// `la::mat<double, 4, 3>` — 4 rows, 3 cols
1563#[derive(Clone, Copy, Debug, Default, PartialEq)]
1564pub struct Mat4x3 {
1565    pub x: Vec4,
1566    pub y: Vec4,
1567    pub z: Vec4,
1568}
1569
1570impl Mat4x3 {
1571    #[inline]
1572    pub const fn from_cols(x: Vec4, y: Vec4, z: Vec4) -> Self {
1573        Self { x, y, z }
1574    }
1575    pub fn row(&self, i: usize) -> Vec3 {
1576        Vec3::new(self.x[i], self.y[i], self.z[i])
1577    }
1578}
1579
1580impl Mul<Vec3> for Mat4x3 {
1581    type Output = Vec4;
1582    fn mul(self, b: Vec3) -> Vec4 {
1583        self.x * b.x + self.y * b.y + self.z * b.z
1584    }
1585}
1586
1587/// `la::mat<double, 3, 2>` — 3 rows, 2 cols
1588#[derive(Clone, Copy, Debug, Default, PartialEq)]
1589pub struct Mat3x2 {
1590    pub x: Vec3,
1591    pub y: Vec3,
1592}
1593
1594impl Mat3x2 {
1595    #[inline]
1596    pub const fn from_cols(x: Vec3, y: Vec3) -> Self {
1597        Self { x, y }
1598    }
1599    pub fn row(&self, i: usize) -> Vec2 {
1600        Vec2::new(self.x[i], self.y[i])
1601    }
1602}
1603
1604impl Mul<Vec2> for Mat3x2 {
1605    type Output = Vec3;
1606    fn mul(self, b: Vec2) -> Vec3 {
1607        self.x * b.x + self.y * b.y
1608    }
1609}
1610
1611/// `la::mat<double, 2, 3>` — 2 rows, 3 cols
1612#[derive(Clone, Copy, Debug, Default, PartialEq)]
1613pub struct Mat2x3 {
1614    pub x: Vec2,
1615    pub y: Vec2,
1616    pub z: Vec2,
1617}
1618
1619impl Mat2x3 {
1620    #[inline]
1621    pub const fn from_cols(x: Vec2, y: Vec2, z: Vec2) -> Self {
1622        Self { x, y, z }
1623    }
1624}
1625
1626impl Mul<Vec3> for Mat2x3 {
1627    type Output = Vec2;
1628    fn mul(self, b: Vec3) -> Vec2 {
1629        self.x * b.x + self.y * b.y + self.z * b.z
1630    }
1631}
1632
1633// ─── Helper conversions ───────────────────────────────────────────────────────
1634
1635/// Embed Mat3x4 into Mat4x4 by appending bottom row [0,0,0,1]
1636pub fn mat3x4_to_mat4(m: Mat3x4) -> Mat4 {
1637    Mat4::from_cols(
1638        Vec4::new(m.x.x, m.x.y, m.x.z, 0.0),
1639        Vec4::new(m.y.x, m.y.y, m.y.z, 0.0),
1640        Vec4::new(m.z.x, m.z.y, m.z.z, 0.0),
1641        Vec4::new(m.w.x, m.w.y, m.w.z, 1.0),
1642    )
1643}
1644
1645/// Extract upper-left 3 rows from a Mat4x4 (drops 4th row)
1646pub fn mat4_to_mat3x4(m: Mat4) -> Mat3x4 {
1647    Mat3x4::from_cols(m.x.xyz(), m.y.xyz(), m.z.xyz(), m.w.xyz())
1648}
1649
1650// ─── Vector algebra functions ─────────────────────────────────────────────────
1651
1652/// 2D cross product: `a.x*b.y - a.y*b.x`
1653#[inline]
1654pub fn cross2(a: Vec2, b: Vec2) -> f64 {
1655    a.x * b.y - a.y * b.x
1656}
1657
1658/// 2D: rotate vector 90° CCW by scalar `a` (scalar × vec2 cross)
1659#[inline]
1660pub fn cross_sv(a: f64, b: Vec2) -> Vec2 {
1661    Vec2::new(-a * b.y, a * b.x)
1662}
1663
1664/// 2D: vec2 × scalar
1665#[inline]
1666pub fn cross_vs(a: Vec2, b: f64) -> Vec2 {
1667    Vec2::new(a.y * b, -a.x * b)
1668}
1669
1670/// 3D cross product
1671#[inline]
1672pub fn cross(a: Vec3, b: Vec3) -> Vec3 {
1673    Vec3::new(
1674        a.y * b.z - a.z * b.y,
1675        a.z * b.x - a.x * b.z,
1676        a.x * b.y - a.y * b.x,
1677    )
1678}
1679
1680/// Dot product
1681#[inline]
1682pub fn dot2(a: Vec2, b: Vec2) -> f64 {
1683    a.x * b.x + a.y * b.y
1684}
1685#[inline]
1686pub fn dot(a: Vec3, b: Vec3) -> f64 {
1687    a.x * b.x + a.y * b.y + a.z * b.z
1688}
1689#[inline]
1690pub fn dot4(a: Vec4, b: Vec4) -> f64 {
1691    a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w
1692}
1693
1694#[inline]
1695pub fn length2_2(a: Vec2) -> f64 {
1696    dot2(a, a)
1697}
1698#[inline]
1699pub fn length2(a: Vec3) -> f64 {
1700    dot(a, a)
1701}
1702#[inline]
1703pub fn length2_4(a: Vec4) -> f64 {
1704    dot4(a, a)
1705}
1706
1707#[inline]
1708pub fn length_2(a: Vec2) -> f64 {
1709    length2_2(a).sqrt()
1710}
1711#[inline]
1712pub fn length(a: Vec3) -> f64 {
1713    length2(a).sqrt()
1714}
1715#[inline]
1716pub fn length_4(a: Vec4) -> f64 {
1717    length2_4(a).sqrt()
1718}
1719
1720#[inline]
1721pub fn normalize2(a: Vec2) -> Vec2 {
1722    a / length_2(a)
1723}
1724#[inline]
1725pub fn normalize(a: Vec3) -> Vec3 {
1726    a / length(a)
1727}
1728#[inline]
1729pub fn normalize4(a: Vec4) -> Vec4 {
1730    a / length_4(a)
1731}
1732
1733#[inline]
1734pub fn distance2_2(a: Vec2, b: Vec2) -> f64 {
1735    length2_2(b - a)
1736}
1737#[inline]
1738pub fn distance2(a: Vec3, b: Vec3) -> f64 {
1739    length2(b - a)
1740}
1741#[inline]
1742pub fn distance_2(a: Vec2, b: Vec2) -> f64 {
1743    length_2(b - a)
1744}
1745#[inline]
1746pub fn distance(a: Vec3, b: Vec3) -> f64 {
1747    length(b - a)
1748}
1749
1750/// Angle between two unit vectors (clamped to [0, π])
1751#[inline]
1752pub fn uangle(a: Vec3, b: Vec3) -> f64 {
1753    let d = dot(a, b);
1754    if d > 1.0 {
1755        0.0
1756    } else {
1757        math::acos(if d < -1.0 { -1.0 } else { d })
1758    }
1759}
1760
1761/// Angle between two non-unit vectors
1762#[inline]
1763pub fn angle(a: Vec3, b: Vec3) -> f64 {
1764    uangle(normalize(a), normalize(b))
1765}
1766
1767/// 2D rotation: rotate `v` CCW by angle `a` (radians)
1768#[inline]
1769pub fn rot2(a: f64, v: Vec2) -> Vec2 {
1770    let (s, c) = (math::sin(a), math::cos(a));
1771    Vec2::new(v.x * c - v.y * s, v.x * s + v.y * c)
1772}
1773
1774/// Rotate `v` CCW around X axis by `a` radians
1775#[inline]
1776pub fn rotx(a: f64, v: Vec3) -> Vec3 {
1777    let (s, c) = (math::sin(a), math::cos(a));
1778    Vec3::new(v.x, v.y * c - v.z * s, v.y * s + v.z * c)
1779}
1780
1781/// Rotate `v` CCW around Y axis by `a` radians
1782#[inline]
1783pub fn roty(a: f64, v: Vec3) -> Vec3 {
1784    let (s, c) = (math::sin(a), math::cos(a));
1785    Vec3::new(v.x * c + v.z * s, v.y, -v.x * s + v.z * c)
1786}
1787
1788/// Rotate `v` CCW around Z axis by `a` radians
1789#[inline]
1790pub fn rotz(a: f64, v: Vec3) -> Vec3 {
1791    let (s, c) = (math::sin(a), math::cos(a));
1792    Vec3::new(v.x * c - v.y * s, v.x * s + v.y * c, v.z)
1793}
1794
1795/// Normalized linear interpolation
1796#[inline]
1797pub fn nlerp(a: Vec3, b: Vec3, t: f64) -> Vec3 {
1798    normalize(lerp3(a, b, t))
1799}
1800
1801/// Spherical linear interpolation between unit vectors
1802#[inline]
1803pub fn slerp(a: Vec3, b: Vec3, t: f64) -> Vec3 {
1804    let th = uangle(a, b);
1805    if th == 0.0 {
1806        a
1807    } else {
1808        a * math::sin(th * (1.0 - t)) / math::sin(th) + b * math::sin(th * t) / math::sin(th)
1809    }
1810}
1811
1812// ─── Component-wise math functions ───────────────────────────────────────────
1813
1814#[inline]
1815pub fn abs2(a: Vec2) -> Vec2 {
1816    Vec2::new(a.x.abs(), a.y.abs())
1817}
1818#[inline]
1819pub fn abs3(a: Vec3) -> Vec3 {
1820    Vec3::new(a.x.abs(), a.y.abs(), a.z.abs())
1821}
1822#[inline]
1823pub fn abs4(a: Vec4) -> Vec4 {
1824    Vec4::new(a.x.abs(), a.y.abs(), a.z.abs(), a.w.abs())
1825}
1826
1827#[inline]
1828pub fn floor2(a: Vec2) -> Vec2 {
1829    Vec2::new(a.x.floor(), a.y.floor())
1830}
1831#[inline]
1832pub fn floor3(a: Vec3) -> Vec3 {
1833    Vec3::new(a.x.floor(), a.y.floor(), a.z.floor())
1834}
1835#[inline]
1836pub fn floor4(a: Vec4) -> Vec4 {
1837    Vec4::new(a.x.floor(), a.y.floor(), a.z.floor(), a.w.floor())
1838}
1839
1840#[inline]
1841pub fn ceil2(a: Vec2) -> Vec2 {
1842    Vec2::new(a.x.ceil(), a.y.ceil())
1843}
1844#[inline]
1845pub fn ceil3(a: Vec3) -> Vec3 {
1846    Vec3::new(a.x.ceil(), a.y.ceil(), a.z.ceil())
1847}
1848
1849#[inline]
1850pub fn round3(a: Vec3) -> Vec3 {
1851    Vec3::new(a.x.round(), a.y.round(), a.z.round())
1852}
1853#[inline]
1854pub fn round4(a: Vec4) -> Vec4 {
1855    Vec4::new(a.x.round(), a.y.round(), a.z.round(), a.w.round())
1856}
1857
1858#[inline]
1859pub fn sqrt3(a: Vec3) -> Vec3 {
1860    Vec3::new(a.x.sqrt(), a.y.sqrt(), a.z.sqrt())
1861}
1862#[inline]
1863pub fn sqrt4(a: Vec4) -> Vec4 {
1864    Vec4::new(a.x.sqrt(), a.y.sqrt(), a.z.sqrt(), a.w.sqrt())
1865}
1866
1867#[inline]
1868pub fn isfinite3(a: Vec3) -> bool {
1869    a.x.is_finite() && a.y.is_finite() && a.z.is_finite()
1870}
1871#[inline]
1872pub fn isfinite4(a: Vec4) -> bool {
1873    a.x.is_finite() && a.y.is_finite() && a.z.is_finite() && a.w.is_finite()
1874}
1875
1876// ─── Component-wise min/max/clamp ─────────────────────────────────────────────
1877
1878#[inline]
1879pub fn min2(a: Vec2, b: Vec2) -> Vec2 {
1880    Vec2::new(a.x.min(b.x), a.y.min(b.y))
1881}
1882#[inline]
1883pub fn min3(a: Vec3, b: Vec3) -> Vec3 {
1884    Vec3::new(a.x.min(b.x), a.y.min(b.y), a.z.min(b.z))
1885}
1886#[inline]
1887pub fn min4(a: Vec4, b: Vec4) -> Vec4 {
1888    Vec4::new(a.x.min(b.x), a.y.min(b.y), a.z.min(b.z), a.w.min(b.w))
1889}
1890
1891#[inline]
1892pub fn max2(a: Vec2, b: Vec2) -> Vec2 {
1893    Vec2::new(a.x.max(b.x), a.y.max(b.y))
1894}
1895#[inline]
1896pub fn max3(a: Vec3, b: Vec3) -> Vec3 {
1897    Vec3::new(a.x.max(b.x), a.y.max(b.y), a.z.max(b.z))
1898}
1899#[inline]
1900pub fn max4(a: Vec4, b: Vec4) -> Vec4 {
1901    Vec4::new(a.x.max(b.x), a.y.max(b.y), a.z.max(b.z), a.w.max(b.w))
1902}
1903
1904/// Clamp `x` component-wise between `lo` and `hi`
1905#[inline]
1906pub fn clamp3(x: Vec3, lo: Vec3, hi: Vec3) -> Vec3 {
1907    Vec3::new(
1908        x.x.max(lo.x).min(hi.x),
1909        x.y.max(lo.y).min(hi.y),
1910        x.z.max(lo.z).min(hi.z),
1911    )
1912}
1913#[inline]
1914pub fn clamp4(x: Vec4, lo: Vec4, hi: Vec4) -> Vec4 {
1915    Vec4::new(
1916        x.x.max(lo.x).min(hi.x),
1917        x.y.max(lo.y).min(hi.y),
1918        x.z.max(lo.z).min(hi.z),
1919        x.w.max(lo.w).min(hi.w),
1920    )
1921}
1922
1923/// Scalar clamp
1924#[inline]
1925pub fn clamp_s(x: f64, lo: f64, hi: f64) -> f64 {
1926    x.max(lo).min(hi)
1927}
1928
1929/// Linear interpolation: `a*(1-t) + b*t`
1930#[inline]
1931pub fn lerp2(a: Vec2, b: Vec2, t: f64) -> Vec2 {
1932    a * (1.0 - t) + b * t
1933}
1934#[inline]
1935pub fn lerp3(a: Vec3, b: Vec3, t: f64) -> Vec3 {
1936    a * (1.0 - t) + b * t
1937}
1938#[inline]
1939pub fn lerp4(a: Vec4, b: Vec4, t: f64) -> Vec4 {
1940    a * (1.0 - t) + b * t
1941}
1942
1943/// Reductions
1944#[inline]
1945pub fn minelem2(a: Vec2) -> f64 {
1946    a.x.min(a.y)
1947}
1948#[inline]
1949pub fn minelem3(a: Vec3) -> f64 {
1950    a.x.min(a.y).min(a.z)
1951}
1952#[inline]
1953pub fn minelem4(a: Vec4) -> f64 {
1954    a.x.min(a.y).min(a.z).min(a.w)
1955}
1956
1957#[inline]
1958pub fn maxelem2(a: Vec2) -> f64 {
1959    a.x.max(a.y)
1960}
1961#[inline]
1962pub fn maxelem3(a: Vec3) -> f64 {
1963    a.x.max(a.y).max(a.z)
1964}
1965#[inline]
1966pub fn maxelem4(a: Vec4) -> f64 {
1967    a.x.max(a.y).max(a.z).max(a.w)
1968}
1969
1970#[inline]
1971pub fn sum3(a: Vec3) -> f64 {
1972    a.x + a.y + a.z
1973}
1974#[inline]
1975pub fn sum4(a: Vec4) -> f64 {
1976    a.x + a.y + a.z + a.w
1977}
1978
1979/// Index of minimum element (argmin)
1980#[inline]
1981pub fn argmin3(a: Vec3) -> usize {
1982    if a.x <= a.y && a.x <= a.z {
1983        0
1984    } else if a.y <= a.z {
1985        1
1986    } else {
1987        2
1988    }
1989}
1990/// Index of maximum element (argmax)
1991#[inline]
1992pub fn argmax3(a: Vec3) -> usize {
1993    if a.x >= a.y && a.x >= a.z {
1994        0
1995    } else if a.y >= a.z {
1996        1
1997    } else {
1998        2
1999    }
2000}
2001#[inline]
2002pub fn argmax4(a: Vec4) -> usize {
2003    let mut j = 0usize;
2004    for i in 1..4 {
2005        if a[i] > a[j] {
2006            j = i;
2007        }
2008    }
2009    j
2010}
2011
2012// ─── Quaternion functions ─────────────────────────────────────────────────────
2013
2014/// Quaternion conjugate: `{-x, -y, -z, w}`
2015#[inline]
2016pub fn qconj(q: Quat) -> Quat {
2017    Vec4::new(-q.x, -q.y, -q.z, q.w)
2018}
2019
2020/// Quaternion inverse
2021#[inline]
2022pub fn qinv(q: Quat) -> Quat {
2023    qconj(q) / length2_4(q)
2024}
2025
2026/// Quaternion Hamilton product
2027#[inline]
2028pub fn qmul(a: Quat, b: Quat) -> Quat {
2029    Quat::new(
2030        a.x * b.w + a.w * b.x + a.y * b.z - a.z * b.y,
2031        a.y * b.w + a.w * b.y + a.z * b.x - a.x * b.z,
2032        a.z * b.w + a.w * b.z + a.x * b.y - a.y * b.x,
2033        a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z,
2034    )
2035}
2036
2037/// X-axis direction from quaternion: `qrot(q, {1,0,0})`
2038#[inline]
2039pub fn qxdir(q: Quat) -> Vec3 {
2040    Vec3::new(
2041        q.w * q.w + q.x * q.x - q.y * q.y - q.z * q.z,
2042        (q.x * q.y + q.z * q.w) * 2.0,
2043        (q.z * q.x - q.y * q.w) * 2.0,
2044    )
2045}
2046
2047/// Y-axis direction from quaternion: `qrot(q, {0,1,0})`
2048#[inline]
2049pub fn qydir(q: Quat) -> Vec3 {
2050    Vec3::new(
2051        (q.x * q.y - q.z * q.w) * 2.0,
2052        q.w * q.w - q.x * q.x + q.y * q.y - q.z * q.z,
2053        (q.y * q.z + q.x * q.w) * 2.0,
2054    )
2055}
2056
2057/// Z-axis direction from quaternion: `qrot(q, {0,0,1})`
2058#[inline]
2059pub fn qzdir(q: Quat) -> Vec3 {
2060    Vec3::new(
2061        (q.z * q.x + q.y * q.w) * 2.0,
2062        (q.y * q.z - q.x * q.w) * 2.0,
2063        q.w * q.w - q.x * q.x - q.y * q.y + q.z * q.z,
2064    )
2065}
2066
2067/// Rotation matrix from quaternion
2068#[inline]
2069pub fn qmat(q: Quat) -> Mat3 {
2070    Mat3::from_cols(qxdir(q), qydir(q), qzdir(q))
2071}
2072
2073/// Rotate vector `v` by quaternion `q`
2074#[inline]
2075pub fn qrot(q: Quat, v: Vec3) -> Vec3 {
2076    qxdir(q) * v.x + qydir(q) * v.y + qzdir(q) * v.z
2077}
2078
2079/// Rotation angle of a unit quaternion
2080#[inline]
2081pub fn qangle(q: Quat) -> f64 {
2082    math::atan2(length(q.xyz()), q.w) * 2.0
2083}
2084
2085/// Rotation axis of a unit quaternion
2086#[inline]
2087pub fn qaxis(q: Quat) -> Vec3 {
2088    normalize(q.xyz())
2089}
2090
2091/// Quaternion nlerp — shortest path
2092#[inline]
2093pub fn qnlerp(a: Quat, b: Quat, t: f64) -> Quat {
2094    let b2 = if dot4(a, b) < 0.0 { -b } else { b };
2095    normalize4(lerp4(a, b2, t))
2096}
2097
2098/// Quaternion slerp — shortest path
2099#[inline]
2100pub fn qslerp(a: Quat, b: Quat, t: f64) -> Quat {
2101    let b2 = if dot4(a, b) < 0.0 { -b } else { b };
2102    // slerp on Vec4 (unit quaternion treated as unit 4D vector)
2103    let th = {
2104        let d = dot4(a, b2).max(-1.0).min(1.0);
2105        if d > 1.0 {
2106            0.0
2107        } else {
2108            math::acos(d)
2109        }
2110    };
2111    if th == 0.0 {
2112        a
2113    } else {
2114        a * math::sin(th * (1.0 - t)) / math::sin(th) + b2 * math::sin(th * t) / math::sin(th)
2115    }
2116}
2117
2118/// Unit quaternion from axis + angle
2119#[inline]
2120pub fn rotation_quat_axis_angle(axis: Vec3, angle: f64) -> Quat {
2121    Quat::new(
2122        axis.x * math::sin(angle / 2.0),
2123        axis.y * math::sin(angle / 2.0),
2124        axis.z * math::sin(angle / 2.0),
2125        math::cos(angle / 2.0),
2126    )
2127}
2128
2129/// Unit quaternion representing shortest rotation from `orig` to `dest`
2130pub fn rotation_quat_vec(orig: Vec3, dest: Vec3) -> Quat {
2131    let cos_theta = dot(orig, dest);
2132    let eps = f64::EPSILON;
2133    if cos_theta >= 1.0 - eps {
2134        return Quat::new(0.0, 0.0, 0.0, 1.0);
2135    }
2136    if cos_theta < -1.0 + eps {
2137        let mut axis = cross(Vec3::new(0.0, 0.0, 1.0), orig);
2138        if length2(axis) < eps {
2139            axis = cross(Vec3::new(1.0, 0.0, 0.0), orig);
2140        }
2141        return rotation_quat_axis_angle(normalize(axis), std::f64::consts::PI);
2142    }
2143    let axis = cross(orig, dest);
2144    let s = ((1.0 + cos_theta) * 2.0).sqrt();
2145    Quat::new(axis.x / s, axis.y / s, axis.z / s, s * 0.5)
2146}
2147
2148/// Unit quaternion from a rotation matrix
2149pub fn rotation_quat_mat(m: Mat3) -> Quat {
2150    let q = Vec4::new(
2151        m.x.x - m.y.y - m.z.z,
2152        m.y.y - m.x.x - m.z.z,
2153        m.z.z - m.x.x - m.y.y,
2154        m.x.x + m.y.y + m.z.z,
2155    );
2156    // s[argmax(q)] gives the sign correction
2157    let s = [
2158        Vec4::new(1.0, m.x.y + m.y.x, m.z.x + m.x.z, m.y.z - m.z.y),
2159        Vec4::new(m.x.y + m.y.x, 1.0, m.y.z + m.z.y, m.z.x - m.x.z),
2160        Vec4::new(m.x.z + m.z.x, m.y.z + m.z.y, 1.0, m.x.y - m.y.x),
2161        Vec4::new(m.y.z - m.z.y, m.z.x - m.x.z, m.x.y - m.y.x, 1.0),
2162    ];
2163    let idx = argmax4(q);
2164    // copysign(normalize(sqrt(max(0, 1+q))), s[idx])
2165    let sq = Vec4::new(
2166        (0.0f64).max(1.0 + q.x).sqrt(),
2167        (0.0f64).max(1.0 + q.y).sqrt(),
2168        (0.0f64).max(1.0 + q.z).sqrt(),
2169        (0.0f64).max(1.0 + q.w).sqrt(),
2170    );
2171    let n = normalize4(sq);
2172    let si = s[idx];
2173    Vec4::new(
2174        n.x.copysign(si.x),
2175        n.y.copysign(si.y),
2176        n.z.copysign(si.z),
2177        n.w.copysign(si.w),
2178    )
2179}
2180
2181// ─── Matrix factory functions ─────────────────────────────────────────────────
2182
2183pub fn translation_matrix(t: Vec3) -> Mat4 {
2184    Mat4::from_cols(
2185        Vec4::new(1.0, 0.0, 0.0, 0.0),
2186        Vec4::new(0.0, 1.0, 0.0, 0.0),
2187        Vec4::new(0.0, 0.0, 1.0, 0.0),
2188        Vec4::new(t.x, t.y, t.z, 1.0),
2189    )
2190}
2191
2192pub fn rotation_matrix(q: Quat) -> Mat4 {
2193    Mat4::from_cols(
2194        Vec4::from((qxdir(q), 0.0)),
2195        Vec4::from((qydir(q), 0.0)),
2196        Vec4::from((qzdir(q), 0.0)),
2197        Vec4::new(0.0, 0.0, 0.0, 1.0),
2198    )
2199}
2200
2201pub fn scaling_matrix(s: Vec3) -> Mat4 {
2202    Mat4::from_cols(
2203        Vec4::new(s.x, 0.0, 0.0, 0.0),
2204        Vec4::new(0.0, s.y, 0.0, 0.0),
2205        Vec4::new(0.0, 0.0, s.z, 0.0),
2206        Vec4::new(0.0, 0.0, 0.0, 1.0),
2207    )
2208}
2209
2210pub fn pose_matrix(q: Quat, p: Vec3) -> Mat4 {
2211    Mat4::from_cols(
2212        Vec4::from((qxdir(q), 0.0)),
2213        Vec4::from((qydir(q), 0.0)),
2214        Vec4::from((qzdir(q), 0.0)),
2215        Vec4::new(p.x, p.y, p.z, 1.0),
2216    )
2217}
2218
2219/// Outer product: vec3 ⊗ vec3 → Mat3
2220pub fn outerprod(a: Vec3, b: Vec3) -> Mat3 {
2221    Mat3::from_cols(a * b.x, a * b.y, a * b.z)
2222}
2223
2224// ─── Hash impls for f64 vectors (using bit-cast) ──────────────────────────────
2225// Matches C++ hash: h(v.x) ^ (h(v.y) << 1) ^ ...
2226
2227fn hash_f64<H: Hasher>(v: f64, state: &mut H) {
2228    // Use bit representation for hashing (NaN will hash consistently)
2229    v.to_bits().hash(state);
2230}
2231
2232impl Hash for Vec2 {
2233    fn hash<H: Hasher>(&self, state: &mut H) {
2234        hash_f64(self.x, state);
2235        // XOR with shift — mirrors C++ std::hash specialization
2236        let mut h2 = std::collections::hash_map::DefaultHasher::new();
2237        hash_f64(self.y, &mut h2);
2238        state.write_u64(std::hash::BuildHasher::build_hasher(&std::collections::hash_map::RandomState::new()).finish() ^ (std::hash::Hasher::finish(&h2) << 1));
2239    }
2240}
2241
2242// Simpler hash that is consistent: just hash all fields sequentially
2243impl Hash for Vec3 {
2244    fn hash<H: Hasher>(&self, state: &mut H) {
2245        self.x.to_bits().hash(state);
2246        self.y.to_bits().hash(state);
2247        self.z.to_bits().hash(state);
2248    }
2249}
2250
2251impl Hash for Vec4 {
2252    fn hash<H: Hasher>(&self, state: &mut H) {
2253        self.x.to_bits().hash(state);
2254        self.y.to_bits().hash(state);
2255        self.z.to_bits().hash(state);
2256        self.w.to_bits().hash(state);
2257    }
2258}
2259
2260// ─── Tests ────────────────────────────────────────────────────────────────────
2261
2262
2263#[cfg(test)]
2264#[path = "linalg_tests.rs"]
2265mod tests;
manifold_rust/linalg.rs

manifold_rust/
linalg.rs
