ehmm... yeah about that now i have a rust crate

FILE ./angle/mod.rs ----

use std::f32::consts::{PI, TAU};

#[derive(Debug, Clone, Copy, PartialEq, PartialOrd)]
pub struct Radians(pub f32);

#[derive(Debug, Clone, Copy, PartialEq, PartialOrd)]
pub struct Degrees(pub f32);

#[derive(Debug, Clone, Copy, PartialEq, PartialOrd)]
pub struct Turns(pub f32);

#[derive(Debug, Clone, Copy, PartialEq, PartialOrd)]
pub struct Gradians(pub f32);


impl From<Degrees> for Radians {
    fn from(d: Degrees) -> Self { Radians(d.0.to_radians()) }
}

impl From<Radians> for Degrees {
    fn from(r: Radians) -> Self { Degrees(r.0.to_degrees()) }
}

impl From<Turns> for Radians {
    fn from(t: Turns) -> Self { Radians(t.0 * TAU) }
}

impl From<Radians> for Turns {
    fn from(r: Radians) -> Self { Turns(r.0 / TAU) }
}

impl From<Gradians> for Radians {
    fn from(g: Gradians) -> Self { Radians(g.0 * PI / 200.0) }
}

impl From<Radians> for Gradians {
    fn from(r: Radians) -> Self { Gradians(r.0 * 200.0 / PI) }
}


impl From<Degrees> for Turns {
    fn from(d: Degrees) -> Self { Radians::from(d).into() }
}

impl From<Turns> for Degrees {
    fn from(t: Turns) -> Self { Radians::from(t).into() }
}

impl From<Degrees> for Gradians {
    fn from(d: Degrees) -> Self { Radians::from(d).into() }
}

impl From<Gradians> for Degrees {
    fn from(g: Gradians) -> Self { Radians::from(g).into() }
}

impl From<Turns> for Gradians {
    fn from(t: Turns) -> Self { Radians::from(t).into() }
}

impl From<Gradians> for Turns {
    fn from(g: Gradians) -> Self { Radians::from(g).into() }
}


impl From<Radians> for f32 {
    fn from(r: Radians) -> f32 { r.0 }
}

impl From<f32> for Radians {
    fn from(f: f32) -> Radians { Radians(f) }
}
FILE ./lib.rs ----
pub mod vec;
pub mod matrix;
pub mod angle;
pub mod simd;
pub mod quat;
FILE ./matrix/mat2.rs ----
FILE ./matrix/mat3.rs ----
FILE ./matrix/mat4.rs ----
use std::ops::{
    Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign,
};

use crate::simd::*;
use crate::{angle::Radians, vec::Vec3};

#[repr(C)]
#[derive(Clone, Copy, Debug)]
pub struct Matrix4 {
    pub cols: [[f32; 4]; 4],
}

impl Matrix4 {
    fn as_flat(&self) -> &[f32] {
        unsafe { std::slice::from_raw_parts(self.cols.as_ptr() as *const f32, 16) }
    }

    fn as_flat_mut(&mut self) -> &mut [f32] {
        unsafe { std::slice::from_raw_parts_mut(self.cols.as_mut_ptr() as *mut f32, 16) }
    }

    pub const IDENTITY: Matrix4 = Matrix4 {
        cols: [
            [1.0, 0.0, 0.0, 0.0],
            [0.0, 1.0, 0.0, 0.0],
            [0.0, 0.0, 1.0, 0.0],
            [0.0, 0.0, 0.0, 1.0],
        ],
    };

    pub fn from_translation(t: Vec3) -> Self {
        let mut tmp = Self::IDENTITY;
        tmp[3] = [t.x, t.y, t.z, 1.0];
        tmp
    }

    pub fn from_scale(s: Vec3) -> Self {
        let mut tmp = Self::IDENTITY;
        tmp[0][0] = s.x;
        tmp[1][1] = s.y;
        tmp[2][2] = s.z;
        tmp
    }

    pub fn from_rotation_z(angle: impl Into<Radians>) -> Self {
        let Radians(rad) = angle.into();
        let mut tmp = Self::IDENTITY;

        let c = rad.cos();
        let s = rad.sin();
        tmp[0][0] = c;
        tmp[0][1] = s;
        tmp[1][0] = -s;
        tmp[1][1] = c;
        tmp
    }

    pub fn from_rotation_x(angle: impl Into<Radians>) -> Self {
        let Radians(rad) = angle.into();
        let mut tmp = Self::IDENTITY;

        let c = rad.cos();
        let s = rad.sin();
        tmp[1][1] = c;
        tmp[1][2] = s;
        tmp[2][1] = -s;
        tmp[2][2] = c;
        tmp
    }

    pub fn from_rotation_y(angle: impl Into<Radians>) -> Self {
        let Radians(rad) = angle.into();
        let mut tmp = Self::IDENTITY;

        let c = rad.cos();
        let s = rad.sin();
        tmp[0][0] = c;
        tmp[0][2] = -s;
        tmp[2][0] = s;
        tmp[2][2] = c;
        tmp
    }

    pub fn from_axis_angle(axis: Vec3, angle: impl Into<Radians>) -> Self {
        let Radians(rad) = angle.into();
        let mut tmp = Self::IDENTITY;
        let c = rad.cos();
        let s = rad.sin();
        let t = 1.0 - c;
        let x = axis.x;
        let y = axis.y;
        let z = axis.z;

        tmp[0][0] = t * x * x + c;
        tmp[0][1] = t * x * y + s * z;
        tmp[0][2] = t * x * z - s * y;

        tmp[1][0] = t * x * y - s * z;
        tmp[1][1] = t * y * y + c;
        tmp[1][2] = t * y * z + s * x;

        tmp[2][0] = t * x * z + s * y;
        tmp[2][1] = t * y * z - s * x;
        tmp[2][2] = t * z * z + c;

        tmp
    }

    fn mul_impl(a: &Matrix4, b: &Matrix4) -> Matrix4 {
        let mut out = Matrix4 {
            cols: [[0.0; 4]; 4],
        };
        for c in 0..4 {
            for r in 0..4 {
                let mut sum = 0.0;
                for k in 0..4 {
                    sum += a.cols[k][r] * b.cols[c][k];
                }
                out.cols[c][r] = sum;
            }
        }
        out
    }
}

macro_rules! impl_mat4_binop {
  ...
}

macro_rules! impl_mat4_assignop {
   ...
}

macro_rules! impl_mat4_scalar {
    ...
}

macro_rules! impl_mat4_scalar_assign {
    ($trait:ident, $method:ident, $simd_fn:ident) => {
        impl $trait<f32> for Matrix4 {
            fn $method(&mut self, rhs: f32) {
                $simd_fn(self.as_flat_mut(), rhs);
            }
        }
    };
}

macro_rules! impl_mat4_neg {
    () => {
        impl Neg for Matrix4 {
            type Output = Matrix4;
            fn neg(self) -> Matrix4 {
                let mut out = Matrix4 { cols: [[0.0; 4]; 4] };
                neg_f32_slice(self.as_flat(), out.as_flat_mut());
                out
            }
        }
        impl Neg for &Matrix4 {
            type Output = Matrix4;
            fn neg(self) -> Matrix4 {
                let mut out = Matrix4 { cols: [[0.0; 4]; 4] };
                neg_f32_slice(self.as_flat(), out.as_flat_mut());
                out
            }
        }
        impl Neg for &mut Matrix4 {
            type Output = Matrix4;
            fn neg(self) -> Matrix4 {
                let mut out = Matrix4 { cols: [[0.0; 4]; 4] };
                neg_f32_slice(self.as_flat(), out.as_flat_mut());
                out
            }
        }
    };
}

macro_rules! impl_mat4_mul {
    () => {
        impl Mul<Matrix4> for Matrix4 {
            type Output = Matrix4;
            fn mul(self, rhs: Matrix4) -> Matrix4 {
                let mut out = Matrix4 { cols: [[0.0; 4]; 4] };
                mat4_mul_f32_slices(self.as_flat(), rhs.as_flat(), out.as_flat_mut());
                out
            }
        }
        impl Mul<&Matrix4> for Matrix4 {
            type Output = Matrix4;
            fn mul(self, rhs: &Matrix4) -> Matrix4 {
                let mut out = Matrix4 { cols: [[0.0; 4]; 4] };
                mat4_mul_f32_slices(self.as_flat(), rhs.as_flat(), out.as_flat_mut());
                out
            }
        }
        impl Mul<&mut Matrix4> for Matrix4 {
            type Output = Matrix4;
            fn mul(self, rhs: &mut Matrix4) -> Matrix4 {
                let mut out = Matrix4 { cols: [[0.0; 4]; 4] };
                mat4_mul_f32_slices(self.as_flat(), rhs.as_flat(), out.as_flat_mut());
                out
            }
        }
        impl Mul<Matrix4> for &Matrix4 {
            type Output = Matrix4;
            fn mul(self, rhs: Matrix4) -> Matrix4 {
                let mut out = Matrix4 { cols: [[0.0; 4]; 4] };
                mat4_mul_f32_slices(self.as_flat(), rhs.as_flat(), out.as_flat_mut());
                out
            }
        }
        impl Mul<&Matrix4> for &Matrix4 {
            type Output = Matrix4;
            fn mul(self, rhs: &Matrix4) -> Matrix4 {
                let mut out = Matrix4 { cols: [[0.0; 4]; 4] };
                mat4_mul_f32_slices(self.as_flat(), rhs.as_flat(), out.as_flat_mut());
                out
            }
        }
        impl Mul<&mut Matrix4> for &Matrix4 {
            type Output = Matrix4;
            fn mul(self, rhs: &mut Matrix4) -> Matrix4 {
                let mut out = Matrix4 { cols: [[0.0; 4]; 4] };
                mat4_mul_f32_slices(self.as_flat(), rhs.as_flat(), out.as_flat_mut());
                out
            }
        }
        impl Mul<Matrix4> for &mut Matrix4 {
            type Output = Matrix4;
            fn mul(self, rhs: Matrix4) -> Matrix4 {
                let mut out = Matrix4 { cols: [[0.0; 4]; 4] };
                mat4_mul_f32_slices(self.as_flat(), rhs.as_flat(), out.as_flat_mut());
                out
            }
        }
        impl Mul<&Matrix4> for &mut Matrix4 {
            type Output = Matrix4;
            fn mul(self, rhs: &Matrix4) -> Matrix4 {
                let mut out = Matrix4 { cols: [[0.0; 4]; 4] };
                mat4_mul_f32_slices(self.as_flat(), rhs.as_flat(), out.as_flat_mut());
                out
            }
        }
        impl Mul<&mut Matrix4> for &mut Matrix4 {
            type Output = Matrix4;
            fn mul(self, rhs: &mut Matrix4) -> Matrix4 {
                let mut out = Matrix4 { cols: [[0.0; 4]; 4] };
                mat4_mul_f32_slices(self.as_flat(), rhs.as_flat(), out.as_flat_mut());
                out
            }
        }
    };
}

impl_mat4_binop!(Add, add, add_f32_slices);
impl_mat4_binop!(Sub, sub, sub_f32_slices);

impl_mat4_assignop!(AddAssign, add_assign, add_assign_f32_slices);
impl_mat4_assignop!(SubAssign, sub_assign, sub_assign_f32_slices);

impl_mat4_scalar!(Mul, mul, mul_f32_slice_scalar);
impl_mat4_scalar!(Div, div, div_f32_slice_scalar);

impl_mat4_scalar_assign!(MulAssign, mul_assign, mul_assign_f32_slice_scalar);
impl_mat4_scalar_assign!(DivAssign, div_assign, div_assign_f32_slice_scalar);

impl_mat4_neg!();
impl_mat4_mul!();

impl Index<usize> for Matrix4 {
    type Output = [f32; 4];
    fn index(&self, c: usize) -> &[f32; 4] {
        &self.cols[c]
    }
}

impl IndexMut<usize> for Matrix4 {
    fn index_mut(&mut self, c: usize) -> &mut [f32; 4] {
        &mut self.cols[c]
    }
}
FILE ./matrix/mod.rs ----
mod mat2;
mod mat3;
mod mat4;

pub use mat2::*;
pub use mat3::*;
pub use mat4::*;
FILE ./quat/mod.rs ----
use crate::{angle::Radians, vec::Vec3};
use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign};

#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Quat {
    w: f32,
    x: f32,
    y: f32,
    z: f32,
}

impl Quat {
    pub const IDENTITY: Quat = Quat {
        w: 1.0,
        x: 0.0,
        y: 0.0,
        z: 0.0,
    };

    pub fn new(w: f32, x: f32, y: f32, z: f32) -> Self {
        Self { w, x, y, z }
    }

    pub fn from_axis_angle(axis: Vec3, angle: impl Into<Radians>) -> Self {
        let Radians(rad) = angle.into();
        let axis = axis.normalize();
        Self {
            w: (rad / 2.0).cos(),
            x: axis.x * (rad / 2.0).sin(),
            y: axis.y * (rad / 2.0).sin(),
            z: axis.z * (rad / 2.0).sin(),
        }
    }

    pub fn length_square(&self) -> f32 {
        self.w * self.w + self.x * self.x + self.y * self.y + self.z * self.z
    }

    pub fn length(&self) -> f32 {
        (self.w * self.w + self.x * self.x + self.y * self.y + self.z * self.z).sqrt()
    }

    pub fn normalize(&self) -> Self {
        let len = self.length();
        if len < 1e-6 {
            return Quat::IDENTITY;
        };
        Self {
            w: self.w / len,
            x: self.x / len,
            y: self.y / len,
            z: self.z / len,
        }
    }

    pub fn normalize_in_place(&mut self) {
        let len = self.length();
        if len < 1e-6 {
            self.w = 1.0;
            self.x = 0.0;
            self.y = 0.0;
            self.z = 0.0;
            return;
        }
        self.w /= len;
        self.x /= len;
        self.y /= len;
        self.z /= len;
    }

    pub fn conjugate(&self) -> Self {
        Self {
            w: self.w,
            x: -self.x,
            y: -self.y,
            z: -self.z,
        }
    }

    pub fn conjugate_in_place(&mut self) {
        self.x = -self.x;
        self.y = -self.y;
        self.z = -self.z;
    }

    pub fn inverse(&self) -> Self {
        let len_sq = self.length_square();
        if len_sq < 1e-6 {
            return Quat::IDENTITY;
        }
        let con = self.conjugate();
        Self {
            w: con.w / len_sq,
            x: con.x / len_sq,
            y: con.y / len_sq,
            z: con.z / len_sq,
        }
    }

    pub fn inverse_in_place(&mut self) {
        let len_sq = self.length_square();
        if len_sq < 1e-6 {
            self.w = 1.0;
            self.x = 0.0;
            self.y = 0.0;
            self.z = 0.0;
            return;
        }
        self.w /= len_sq;
        self.x = -self.x / len_sq;
        self.y = -self.y / len_sq;
        self.z = -self.z / len_sq;
    }
}


macro_rules! impl_quat_binop {
    ($trait:ident, $method:ident, $op:tt) => {
        impl $trait<Quat> for Quat {
            type Output = Quat;
            fn $method(self, rhs: Quat) -> Quat {
                Quat { w: self.w $op rhs.w, x: self.x $op rhs.x,
                       y: self.y $op rhs.y, z: self.z $op rhs.z }
            }
        }
        impl $trait<&Quat> for Quat {
            type Output = Quat;
            fn $method(self, rhs: &Quat) -> Quat {
                Quat { w: self.w $op rhs.w, x: self.x $op rhs.x,
                       y: self.y $op rhs.y, z: self.z $op rhs.z }
            }
        }
        impl $trait<&mut Quat> for Quat {
            type Output = Quat;
            fn $method(self, rhs: &mut Quat) -> Quat {
                Quat { w: self.w $op rhs.w, x: self.x $op rhs.x,
                       y: self.y $op rhs.y, z: self.z $op rhs.z }
            }
        }
        impl $trait<Quat> for &Quat {
            type Output = Quat;
            fn $method(self, rhs: Quat) -> Quat {
                Quat { w: self.w $op rhs.w, x: self.x $op rhs.x,
                       y: self.y $op rhs.y, z: self.z $op rhs.z }
            }
        }
        impl $trait<&Quat> for &Quat {
            type Output = Quat;
            fn $method(self, rhs: &Quat) -> Quat {
                Quat { w: self.w $op rhs.w, x: self.x $op rhs.x,
                       y: self.y $op rhs.y, z: self.z $op rhs.z }
            }
        }
        impl $trait<&mut Quat> for &Quat {
            type Output = Quat;
            fn $method(self, rhs: &mut Quat) -> Quat {
                Quat { w: self.w $op rhs.w, x: self.x $op rhs.x,
                       y: self.y $op rhs.y, z: self.z $op rhs.z }
            }
        }
        impl $trait<Quat> for &mut Quat {
            type Output = Quat;
            fn $method(self, rhs: Quat) -> Quat {
                Quat { w: self.w $op rhs.w, x: self.x $op rhs.x,
                       y: self.y $op rhs.y, z: self.z $op rhs.z }
            }
        }
        impl $trait<&Quat> for &mut Quat {
            type Output = Quat;
            fn $method(self, rhs: &Quat) -> Quat {
                Quat { w: self.w $op rhs.w, x: self.x $op rhs.x,
                       y: self.y $op rhs.y, z: self.z $op rhs.z }
            }
        }
        impl $trait<&mut Quat> for &mut Quat {
            type Output = Quat;
            fn $method(self, rhs: &mut Quat) -> Quat {
                Quat { w: self.w $op rhs.w, x: self.x $op rhs.x,
                       y: self.y $op rhs.y, z: self.z $op rhs.z }
            }
        }
    };
}

macro_rules! impl_quat_assignop {
    ($trait:ident, $method:ident, $op:tt) => {
        impl $trait<Quat> for Quat {
            fn $method(&mut self, rhs: Quat) {
                self.w $op rhs.w; self.x $op rhs.x;
                self.y $op rhs.y; self.z $op rhs.z;
            }
        }
        impl $trait<&Quat> for Quat {
            fn $method(&mut self, rhs: &Quat) {
                self.w $op rhs.w; self.x $op rhs.x;
                self.y $op rhs.y; self.z $op rhs.z;
            }
        }
        impl $trait<&mut Quat> for Quat {
            fn $method(&mut self, rhs: &mut Quat) {
                self.w $op rhs.w; self.x $op rhs.x;
                self.y $op rhs.y; self.z $op rhs.z;
            }
        }
    };
}

macro_rules! impl_quat_scalar {
    ($trait:ident, $method:ident, $op:tt) => {
        impl $trait<f32> for Quat {
            type Output = Quat;
            fn $method(self, rhs: f32) -> Quat {
                Quat { w: self.w $op rhs, x: self.x $op rhs,
                       y: self.y $op rhs, z: self.z $op rhs }
            }
        }
        impl $trait<f32> for &Quat {
            type Output = Quat;
            fn $method(self, rhs: f32) -> Quat {
                Quat { w: self.w $op rhs, x: self.x $op rhs,
                       y: self.y $op rhs, z: self.z $op rhs }
            }
        }
        impl $trait<f32> for &mut Quat {
            type Output = Quat;
            fn $method(self, rhs: f32) -> Quat {
                Quat { w: self.w $op rhs, x: self.x $op rhs,
                       y: self.y $op rhs, z: self.z $op rhs }
            }
        }
    };
}

macro_rules! impl_quat_scalar_assign {
    ($trait:ident, $method:ident, $op:tt) => {
        impl $trait<f32> for Quat {
            fn $method(&mut self, rhs: f32) {
                self.w $op rhs; self.x $op rhs;
                self.y $op rhs; self.z $op rhs;
            }
        }
    };
}

macro_rules! impl_quat_mul {
    ($lhs:ty, $rhs:ty) => {
        impl Mul<$rhs> for $lhs {
            type Output = Quat;
            fn mul(self, rhs: $rhs) -> Quat {
                Quat {
                    w: self.w * rhs.w - self.x * rhs.x - self.y * rhs.y - self.z * rhs.z,
                    x: self.w * rhs.x + self.x * rhs.w + self.y * rhs.z - self.z * rhs.y,
                    y: self.w * rhs.y - self.x * rhs.z + self.y * rhs.w + self.z * rhs.x,
                    z: self.w * rhs.z + self.x * rhs.y - self.y * rhs.x + self.z * rhs.w,
                }
            }
        }
    };
}

impl_quat_mul!(Quat, Quat);
impl_quat_mul!(Quat, &Quat);
impl_quat_mul!(Quat, &mut Quat);
impl_quat_mul!(&Quat, Quat);
impl_quat_mul!(&Quat, &Quat);
impl_quat_mul!(&Quat, &mut Quat);
impl_quat_mul!(&mut Quat, Quat);
impl_quat_mul!(&mut Quat, &Quat);
impl_quat_mul!(&mut Quat, &mut Quat);

impl MulAssign<Quat> for Quat {
    fn mul_assign(&mut self, rhs: Quat) {
        *self = &*self * rhs;
    }
}
impl MulAssign<&Quat> for Quat {
    fn mul_assign(&mut self, rhs: &Quat) {
        *self = &*self * rhs;
    }
}
impl MulAssign<&mut Quat> for Quat {
    fn mul_assign(&mut self, rhs: &mut Quat) {
        *self = &*self * &*rhs;
    }
}

impl Neg for Quat {
    type Output = Quat;
    fn neg(self) -> Quat {
        Quat {
            w: -self.w,
            x: -self.x,
            y: -self.y,
            z: -self.z,
        }
    }
}
impl Neg for &Quat {
    type Output = Quat;
    fn neg(self) -> Quat {
        Quat {
            w: -self.w,
            x: -self.x,
            y: -self.y,
            z: -self.z,
        }
    }
}
impl Neg for &mut Quat {
    type Output = Quat;
    fn neg(self) -> Quat {
        Quat {
            w: -self.w,
            x: -self.x,
            y: -self.y,
            z: -self.z,
        }
    }
}

impl_quat_binop!(Add, add, +);
impl_quat_binop!(Sub, sub, -);

impl_quat_assignop!(AddAssign, add_assign, +=);
impl_quat_assignop!(SubAssign, sub_assign, -=);

impl_quat_scalar!(Mul, mul, *);
impl_quat_scalar!(Div, div, /);

impl_quat_scalar_assign!(MulAssign, mul_assign, *=);
impl_quat_scalar_assign!(DivAssign, div_assign, /=);

impl PartialOrd for Quat {
    fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
        self.length_square().partial_cmp(&other.length_square())
    }
    fn lt(&self, other: &Self) -> bool {
        self.length_square() < other.length_square()
    }
    fn gt(&self, other: &Self) -> bool {
        self.length_square() > other.length_square()
    }
    fn le(&self, other: &Self) -> bool {
        self.length_square() <= other.length_square()
    }
    fn ge(&self, other: &Self) -> bool {
        self.length_square() >= other.length_square()
    }
}
FILE ./simd/mod.rs ----
use std::arch::x86_64::*;

pub fn add_f32_slices(a: &[f32], b: &[f32], out: &mut [f32]) {
    assert_eq!(a.len(), b.len());
    assert_eq!(a.len(), out.len());
    #[cfg(target_arch = "x86_64")]
    {
        if is_x86_feature_detected!("sse2") {
            unsafe { add_f32_sse2(a, b, out) };
            return;
        }
    }
    add_f32_scalar(a, b, out);
}

#[inline]
fn add_f32_scalar(a: &[f32], b: &[f32], out: &mut [f32]) {
    for i in 0..a.len() {
        out[i] = a[i] + b[i];
    }
}

#[target_feature(enable = "sse2")]
unsafe fn add_f32_sse2(a: &[f32], b: &[f32], out: &mut [f32]) {
    let chunks = a.len() / 4;
    for i in 0..chunks {
        let off = i * 4;
        let va = _mm_loadu_ps(a.as_ptr().add(off));
        let vb = _mm_loadu_ps(b.as_ptr().add(off));
        _mm_storeu_ps(out.as_mut_ptr().add(off), _mm_add_ps(va, vb));
    }
    let rem = chunks * 4;
    add_f32_scalar(&a[rem..], &b[rem..], &mut out[rem..]);
}

pub fn sub_f32_slices(a: &[f32], b: &[f32], out: &mut [f32]) {
    assert_eq!(a.len(), b.len());
    assert_eq!(a.len(), out.len());
    #[cfg(target_arch = "x86_64")]
    {
        if is_x86_feature_detected!("sse2") {
            unsafe { sub_f32_sse2(a, b, out) };
            return;
        }
    }
    sub_f32_scalar(a, b, out);
}

#[inline]
fn sub_f32_scalar(a: &[f32], b: &[f32], out: &mut [f32]) {
    for i in 0..a.len() {
        out[i] = a[i] - b[i];
    }
}

#[target_feature(enable = "sse2")]
unsafe fn sub_f32_sse2(a: &[f32], b: &[f32], out: &mut [f32]) {
    let chunks = a.len() / 4;
    for i in 0..chunks {
        let off = i * 4;
        let va = _mm_loadu_ps(a.as_ptr().add(off));
        let vb = _mm_loadu_ps(b.as_ptr().add(off));
        _mm_storeu_ps(out.as_mut_ptr().add(off), _mm_sub_ps(va, vb));
    }
    let rem = chunks * 4;
    sub_f32_scalar(&a[rem..], &b[rem..], &mut out[rem..]);
}

pub fn mul_f32_slices(a: &[f32], b: &[f32], out: &mut [f32]) {
    assert_eq!(a.len(), b.len());
    assert_eq!(a.len(), out.len());
    #[cfg(target_arch = "x86_64")]
    {
        if is_x86_feature_detected!("sse2") {
            unsafe { mul_f32_sse2(a, b, out) };
            return;
        }
    }
    mul_f32_scalar(a, b, out);
}

#[inline]
fn mul_f32_scalar(a: &[f32], b: &[f32], out: &mut [f32]) {
    for i in 0..a.len() {
        out[i] = a[i] * b[i];
    }
}

#[target_feature(enable = "sse2")]
unsafe fn mul_f32_sse2(a: &[f32], b: &[f32], out: &mut [f32]) {
    let chunks = a.len() / 4;
    for i in 0..chunks {
        let off = i * 4;
        let va = _mm_loadu_ps(a.as_ptr().add(off));
        let vb = _mm_loadu_ps(b.as_ptr().add(off));
        _mm_storeu_ps(out.as_mut_ptr().add(off), _mm_mul_ps(va, vb));
    }
    let rem = chunks * 4;
    mul_f32_scalar(&a[rem..], &b[rem..], &mut out[rem..]);
}

pub fn div_f32_slices(a: &[f32], b: &[f32], out: &mut [f32]) {
    assert_eq!(a.len(), b.len());
    assert_eq!(a.len(), out.len());
    #[cfg(target_arch = "x86_64")]
    {
        if is_x86_feature_detected!("sse2") {
            unsafe { div_f32_sse2(a, b, out) };
            return;
        }
    }
    div_f32_scalar(a, b, out);
}

#[inline]
fn div_f32_scalar(a: &[f32], b: &[f32], out: &mut [f32]) {
    for i in 0..a.len() {
        out[i] = a[i] / b[i];
    }
}

#[target_feature(enable = "sse2")]
unsafe fn div_f32_sse2(a: &[f32], b: &[f32], out: &mut [f32]) {
    let chunks = a.len() / 4;
    for i in 0..chunks {
        let off = i * 4;
        let va = _mm_loadu_ps(a.as_ptr().add(off));
        let vb = _mm_loadu_ps(b.as_ptr().add(off));
        _mm_storeu_ps(out.as_mut_ptr().add(off), _mm_div_ps(va, vb));
    }
    let rem = chunks * 4;
    div_f32_scalar(&a[rem..], &b[rem..], &mut out[rem..]);
}

pub fn add_assign_f32_slices(a: &mut [f32], b: &[f32]) {
    assert_eq!(a.len(), b.len());
    #[cfg(target_arch = "x86_64")]
    {
        if is_x86_feature_detected!("sse2") {
            unsafe { add_assign_f32_sse2(a, b) };
            return;
        }
    }
    add_assign_f32_scalar(a, b);
}

#[inline]
fn add_assign_f32_scalar(a: &mut [f32], b: &[f32]) {
    for i in 0..a.len() {
        a[i] += b[i];
    }
}

#[target_feature(enable = "sse2")]
unsafe fn add_assign_f32_sse2(a: &mut [f32], b: &[f32]) {
    let chunks = a.len() / 4;
    for i in 0..chunks {
        let off = i * 4;
        let va = _mm_loadu_ps(a.as_ptr().add(off));
        let vb = _mm_loadu_ps(b.as_ptr().add(off));
        _mm_storeu_ps(a.as_mut_ptr().add(off), _mm_add_ps(va, vb));
    }
    let rem = chunks * 4;
    add_assign_f32_scalar(&mut a[rem..], &b[rem..]);
}

pub fn sub_assign_f32_slices(a: &mut [f32], b: &[f32]) {
    assert_eq!(a.len(), b.len());
    #[cfg(target_arch = "x86_64")]
    {
        if is_x86_feature_detected!("sse2") {
            unsafe { sub_assign_f32_sse2(a, b) };
            return;
        }
    }
    sub_assign_f32_scalar(a, b);
}

#[inline]
fn sub_assign_f32_scalar(a: &mut [f32], b: &[f32]) {
    for i in 0..a.len() {
        a[i] -= b[i];
    }
}

#[target_feature(enable = "sse2")]
unsafe fn sub_assign_f32_sse2(a: &mut [f32], b: &[f32]) {
    let chunks = a.len() / 4;
    for i in 0..chunks {
        let off = i * 4;
        let va = _mm_loadu_ps(a.as_ptr().add(off));
        let vb = _mm_loadu_ps(b.as_ptr().add(off));
        _mm_storeu_ps(a.as_mut_ptr().add(off), _mm_sub_ps(va, vb));
    }
    let rem = chunks * 4;
    sub_assign_f32_scalar(&mut a[rem..], &b[rem..]);
}

pub fn mul_assign_f32_slices(a: &mut [f32], b: &[f32]) {
    assert_eq!(a.len(), b.len());
    #[cfg(target_arch = "x86_64")]
    {
        if is_x86_feature_detected!("sse2") {
            unsafe { mul_assign_f32_sse2(a, b) };
            return;
        }
    }
    mul_assign_f32_scalar(a, b);
}

#[inline]
fn mul_assign_f32_scalar(a: &mut [f32], b: &[f32]) {
    for i in 0..a.len() {
        a[i] *= b[i];
    }
}

#[target_feature(enable = "sse2")]
unsafe fn mul_assign_f32_sse2(a: &mut [f32], b: &[f32]) {
    let chunks = a.len() / 4;
    for i in 0..chunks {
        let off = i * 4;
        let va = _mm_loadu_ps(a.as_ptr().add(off));
        let vb = _mm_loadu_ps(b.as_ptr().add(off));
        _mm_storeu_ps(a.as_mut_ptr().add(off), _mm_mul_ps(va, vb));
    }
    let rem = chunks * 4;
    mul_assign_f32_scalar(&mut a[rem..], &b[rem..]);
}

pub fn div_assign_f32_slices(a: &mut [f32], b: &[f32]) {
    assert_eq!(a.len(), b.len());
    #[cfg(target_arch = "x86_64")]
    {
        if is_x86_feature_detected!("sse2") {
            unsafe { div_assign_f32_sse2(a, b) };
            return;
        }
    }
    div_assign_f32_scalar(a, b);
}

#[inline]
fn div_assign_f32_scalar(a: &mut [f32], b: &[f32]) {
    for i in 0..a.len() {
        a[i] /= b[i];
    }
}

#[target_feature(enable = "sse2")]
unsafe fn div_assign_f32_sse2(a: &mut [f32], b: &[f32]) {
    let chunks = a.len() / 4;
    for i in 0..chunks {
        let off = i * 4;
        let va = _mm_loadu_ps(a.as_ptr().add(off));
        let vb = _mm_loadu_ps(b.as_ptr().add(off));
        _mm_storeu_ps(a.as_mut_ptr().add(off), _mm_div_ps(va, vb));
    }
    let rem = chunks * 4;
    div_assign_f32_scalar(&mut a[rem..], &b[rem..]);
}

pub fn neg_f32_slice(a: &[f32], out: &mut [f32]) {
    assert_eq!(a.len(), out.len());
    #[cfg(target_arch = "x86_64")]
    {
        if is_x86_feature_detected!("sse2") {
            unsafe { neg_f32_sse2(a, out) };
            return;
        }
    }
    neg_f32_scalar(a, out);
}

#[inline]
fn neg_f32_scalar(a: &[f32], out: &mut [f32]) {
    for i in 0..a.len() {
        out[i] = -a[i];
    }
}

#[target_feature(enable = "sse2")]
unsafe fn neg_f32_sse2(a: &[f32], out: &mut [f32]) {
    let sign_mask = _mm_set1_ps(-0.0f32);
    let chunks = a.len() / 4;
    for i in 0..chunks {
        let off = i * 4;
        let va = _mm_loadu_ps(a.as_ptr().add(off));
        _mm_storeu_ps(out.as_mut_ptr().add(off), _mm_xor_ps(va, sign_mask));
    }
    let rem = chunks * 4;
    neg_f32_scalar(&a[rem..], &mut out[rem..]);
}

pub fn mul_f32_slice_scalar(a: &[f32], b: f32, out: &mut [f32]) {
    assert_eq!(a.len(), out.len());
    #[cfg(target_arch = "x86_64")]
    {
        if is_x86_feature_detected!("sse2") {
            unsafe { mul_f32_scalar_sse2(a, b, out) };
            return;
        }
    }
    for i in 0..a.len() { out[i] = a[i] * b; }
}

#[target_feature(enable = "sse2")]
unsafe fn mul_f32_scalar_sse2(a: &[f32], b: f32, out: &mut [f32]) {
    let vb = _mm_set1_ps(b);
    let chunks = a.len() / 4;
    for i in 0..chunks {
        let off = i * 4;
        let va = _mm_loadu_ps(a.as_ptr().add(off));
        _mm_storeu_ps(out.as_mut_ptr().add(off), _mm_mul_ps(va, vb));
    }
    let rem = chunks * 4;
    for i in rem..a.len() { out[i] = a[i] * b; }
}

pub fn div_f32_slice_scalar(a: &[f32], b: f32, out: &mut [f32]) {
    assert_eq!(a.len(), out.len());
    #[cfg(target_arch = "x86_64")]
    {
        if is_x86_feature_detected!("sse2") {
            unsafe { div_f32_scalar_sse2(a, b, out) };
            return;
        }
    }
    for i in 0..a.len() { out[i] = a[i] / b; }
}

#[target_feature(enable = "sse2")]
unsafe fn div_f32_scalar_sse2(a: &[f32], b: f32, out: &mut [f32]) {
    let vb = _mm_set1_ps(b);
    let chunks = a.len() / 4;
    for i in 0..chunks {
        let off = i * 4;
        let va = _mm_loadu_ps(a.as_ptr().add(off));
        _mm_storeu_ps(out.as_mut_ptr().add(off), _mm_div_ps(va, vb));
    }
    let rem = chunks * 4;
    for i in rem..a.len() { out[i] = a[i] / b; }
}

pub fn mul_assign_f32_slice_scalar(a: &mut [f32], b: f32) {
    #[cfg(target_arch = "x86_64")]
    {
        if is_x86_feature_detected!("sse2") {
            unsafe { mul_assign_f32_scalar_sse2(a, b) };
            return;
        }
    }
    for x in a.iter_mut() { *x *= b; }
}

#[target_feature(enable = "sse2")]
unsafe fn mul_assign_f32_scalar_sse2(a: &mut [f32], b: f32) {
    let vb = _mm_set1_ps(b);
    let chunks = a.len() / 4;
    for i in 0..chunks {
        let off = i * 4;
        let va = _mm_loadu_ps(a.as_ptr().add(off));
        _mm_storeu_ps(a.as_mut_ptr().add(off), _mm_mul_ps(va, vb));
    }
    let rem = chunks * 4;
    for x in a[rem..].iter_mut() { *x *= b; }
}

pub fn div_assign_f32_slice_scalar(a: &mut [f32], b: f32) {
    #[cfg(target_arch = "x86_64")]
    {
        if is_x86_feature_detected!("sse2") {
            unsafe { div_assign_f32_scalar_sse2(a, b) };
            return;
        }
    }
    for x in a.iter_mut() { *x /= b; }
}

#[target_feature(enable = "sse2")]
unsafe fn div_assign_f32_scalar_sse2(a: &mut [f32], b: f32) {
    let vb = _mm_set1_ps(b);
    let chunks = a.len() / 4;
    for i in 0..chunks {
        let off = i * 4;
        let va = _mm_loadu_ps(a.as_ptr().add(off));
        _mm_storeu_ps(a.as_mut_ptr().add(off), _mm_div_ps(va, vb));
    }
    let rem = chunks * 4;
    for x in a[rem..].iter_mut() { *x /= b; }
}

pub fn mat4_mul_f32_slices(a: &[f32], b: &[f32], out: &mut [f32]) {
    assert_eq!(a.len(), 16);
    assert_eq!(b.len(), 16);
    assert_eq!(out.len(), 16);
    for c in 0..4 {
        for r in 0..4 {
            let mut sum = 0.0f32;
            for k in 0..4 {
                sum += a[k * 4 + r] * b[c * 4 + k];
            }
            out[c * 4 + r] = sum;
        }
    }
}
FILE ./vec/mod.rs ----
mod vec2;
mod vec3;

pub use vec2::*;
pub use vec3::*;
FILE ./vec/vec2.rs ----

use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign};

use crate::angle::Radians;

#[derive(Clone, Copy, Debug)]
pub struct Vec2 {
    pub x: f32,
    pub y: f32,
}

impl Vec2 {
    pub const ZERO: Self = Self { x: 0.0, y: 0.0 };

    pub const ONE: Self = Self { x: 1.0, y: 1.0 };

    pub fn new(x: f32, y: f32) -> Self {
        Self { x, y }
    }

    pub fn length_square(&self) -> f32 {
        self.x * self.x + self.y * self.y
    }

    pub fn length(&self) -> f32 {
        (self.x * self.x + self.y * self.y).sqrt()
    }

    pub fn normalize(&self) -> Self {
        let len = self.length();
        if len <= 1e-6 {
            return Self::ZERO;
        }
        Self {
            x: self.x / len,
            y: self.y / len,
        }
    }

    pub fn normalize_in_place(&mut self) {
        let len = self.length();
        if len <= 1e-6 {
            self.x = 0f32;
            self.y = 0f32;
            return;
        }
        *self /= len;
    }

    pub fn set_mag(&self, mag: f32) -> Self {
        let len = self.length();
        if len < 1e-6 {
            return Self::ZERO;
        }
        Self {
            x: (self.x / len) * mag,
            y: (self.y / len) * mag,
        }
    }

    pub fn set_mag_in_place(&mut self, mag: f32) {
        let len = self.length();
        if len < 1e-6 {
            self.x = 0f32;
            self.y = 0.0;
            return;
        }
        self.x = (self.x / len) * mag;
        self.y = (self.y / len) * mag;
    }

    pub fn to_tuple(&self) -> (f32, f32) {
        (self.x, self.y)
    }

    pub fn from_tuple(t: (f32, f32)) -> Self {
        Self { x: t.0, y: t.1 }
    }

    pub fn dot(&self, rhs: &Vec2) -> f32 {
        self.x * rhs.x + self.y * rhs.y
    }

    pub fn distance(&self, rhs: &Vec2) -> f32 {
        (*self - *rhs).length()
    }

    pub fn distance_squared(&self, rhs: &Vec2) -> f32 {
        (*self - *rhs).length_square()
    }

    pub fn cross(&self, rhs: &Vec2) -> f32 {
        self.x * rhs.y - self.y * rhs.x
    }

    pub fn angle(&self) -> Radians {
        Radians(self.y.atan2(self.x))
    }

    pub fn from_angle(angle: impl Into<Radians>) -> Self {
        let Radians(r) = angle.into();
        Self {
            x: r.cos(),
            y: r.sin(),
        }
    }

    pub fn rotate(&self, angle: impl Into<Radians>) -> Self {
        let Radians(r) = angle.into();
        let (s, c) = r.sin_cos();
        Self {
            x: self.x * c - self.y * s,
            y: self.x * s + self.y * c,
        }
    }

    pub fn perp(&self) -> Self {
        Self {
            x: -self.y,
            y: self.x,
        }
    }

    pub fn lerp(&self, rhs: &Vec2, t: f32) -> Self {
        *self + (*rhs - *self) * t
    }

    pub fn clamp_length(&self, max: f32) -> Self {
        let len_sq = self.length_square();
        if len_sq > max * max {
            *self * (max / len_sq.sqrt())
        } else {
            *self
        }
    }

    pub fn abs(&self) -> Self {
        Self {
            x: self.x.abs(),
            y: self.y.abs(),
        }
    }
}


macro_rules! impl_vec2_binop {
    ...
}

macro_rules! impl_vec2_assignop {
    ...
}

macro_rules! impl_vec2_scalar {
    ...
}

macro_rules! impl_vec2_scalar_assign {
    ...
}

macro_rules! impl_vec2_neg {
    ...
}

impl_vec2_binop!(Add, add, +);
impl_vec2_binop!(Sub, sub, -);
impl_vec2_binop!(Mul, mul, *);
impl_vec2_binop!(Div, div, /);

impl_vec2_assignop!(AddAssign, add_assign, +=);
impl_vec2_assignop!(SubAssign, sub_assign, -=);
impl_vec2_assignop!(MulAssign, mul_assign, *=);
impl_vec2_assignop!(DivAssign, div_assign, /=);

impl_vec2_scalar!(Mul, mul, *);
impl_vec2_scalar!(Div, div, /);

impl_vec2_scalar_assign!(MulAssign, mul_assign, *=);
impl_vec2_scalar_assign!(DivAssign, div_assign, /=);

impl_vec2_neg!();

impl PartialEq for Vec2 {
    fn eq(&self, other: &Self) -> bool {
        self.x == other.x && self.y == other.y
    }

    fn ne(&self, other: &Self) -> bool {
        self.x != other.x || self.y != other.y
    }
}

impl PartialOrd for Vec2 {
    fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
        self.length_square().partial_cmp(&other.length_square())
    }

    fn lt(&self, other: &Self) -> bool {
        self.length_square() < other.length_square()
    }

    fn gt(&self, other: &Self) -> bool {
        self.length_square() > other.length_square()
    }

    fn le(&self, other: &Self) -> bool {
        self.length_square() <= other.length_square()
    }

    fn ge(&self, other: &Self) -> bool {
        self.length_square() >= other.length_square()
    }
}

FILE ./vec/vec3.rs ----

use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign};

use crate::angle::Radians;

#[derive(Clone, Copy, Debug)]
pub struct Vec3 {
    pub x: f32,
    pub y: f32,
    pub z: f32,
}

impl Vec3 {
    pub const ZERO: Self = Self { x: 0.0, y: 0.0, z: 0.0 };

    pub const ONE: Self = Self { x: 1.0, y: 1.0, z: 1.0 };

    pub fn new(x: f32, y: f32, z: f32) -> Self {
        Self { x, y, z }
    }

    pub fn length_square(&self) -> f32 {
        self.x * self.x + self.y * self.y + self.z * self.z
    }

    pub fn length(&self) -> f32 {
        (self.x * self.x + self.y * self.y + self.z * self.z).sqrt()
    }

    pub fn normalize(&self) -> Self {
        let len = self.length();
        if len <= 1e-6 {
            return Self::ZERO;
        }
        Self {
            x: self.x / len,
            y: self.y / len,
            z: self.z / len,
        }
    }

    pub fn normalize_in_place(&mut self) {
        let len = self.length();
        if len <= 1e-6 {
            self.x = 0f32;
            self.y = 0f32;
            self.z = 0f32;
            return;
        }
        *self /= len;
    }

    pub fn set_mag(&self, mag: f32) -> Self {
        let len = self.length();
        if len < 1e-6 {
            return Self::ZERO;
        }
        Self {
            x: (self.x / len) * mag,
            y: (self.y / len) * mag,
            z: (self.z / len) * mag,
        }
    }

    pub fn set_mag_in_place(&mut self, mag: f32) {
        let len = self.length();
        if len < 1e-6 {
            self.x = 0f32;
            self.y = 0.0;
            self.z = 0f32;
            return;
        }
        self.x = (self.x / len) * mag;
        self.y = (self.y / len) * mag;
        self.z = (self.z / len) * mag;
    }

    pub fn to_tuple(&self) -> (f32, f32, f32) {
        (self.x, self.y, self.z)
    }

    pub fn from_tuple(t: (f32, f32, f32)) -> Self {
        Self { x: t.0, y: t.1, z: t.2 }
    }

    pub fn dot(&self, rhs: &Vec3) -> f32 {
        self.x * rhs.x + self.y * rhs.y + self.z * rhs.z
    }

    pub fn distance(&self, rhs: &Vec3) -> f32 {
        (*self - *rhs).length()
    }

    pub fn distance_squared(&self, rhs: &Vec3) -> f32 {
        (*self - *rhs).length_square()
    }

    pub fn cross(&self, rhs: &Vec3) -> Vec3 {
        Vec3::new(
            self.y * rhs.z - self.z * rhs.y,
            self.z * rhs.x - self.x * rhs.z,
            self.x * rhs.y - self.y * rhs.x,
        )
    }

    pub fn angle_xy(&self) -> Radians {
        Radians(self.y.atan2(self.x))
    }

    pub fn angle_xz(&self) -> Radians {
        Radians(self.z.atan2(self.x))
    }

    pub fn angle_yz(&self) -> Radians {
        Radians(self.z.atan2(self.y))
    }

    pub fn from_yaw_pitch(yaw: impl Into<Radians>, pitch: impl Into<Radians>) -> Self {
        let Radians(y) = yaw.into();
        let Radians(p) = pitch.into();
        Vec3 {
            x: y.cos() * p.cos(),
            y: p.sin(),
            z: y.sin() * p.cos(),
        }
    }

    pub fn rotate(&self, axis: Vec3, angle: impl Into<Radians>) -> Self {
        let axis_len_sq = axis.length_square();
        if axis_len_sq < 1e-6 {
            return *self;
        }
        let axis = axis / axis_len_sq.sqrt();
        let Radians(theta) = angle.into();
        let (s, c) = theta.sin_cos();
        let v = self;
        v * c + axis.cross(v) * s + axis * (axis.dot(v) * (1.0 - c))
    }

    pub fn right_hand_perp(&self, up: Vec3) -> Vec3 {
        self.cross(&up)
    }

    pub fn lerp(&self, rhs: &Vec3, t: f32) -> Self {
        *self + (*rhs - *self) * t
    }

    pub fn clamp_length(&self, max: f32) -> Self {
        let len_sq = self.length_square();
        if len_sq > max * max {
            *self * (max / len_sq.sqrt())
        } else {
            *self
        }
    }
}



macro_rules! impl_vec3_binop {
    ...
}

macro_rules! impl_vec3_assignop {
    ...
}

macro_rules! impl_vec3_scalar {
    ...
}

macro_rules! impl_vec3_scalar_assign {
    ...
}

macro_rules! impl_vec3_neg {
    () => {
        impl Neg for Vec3 {
            type Output = Vec3;
            fn neg(self) -> Vec3 {
                Vec3 {
                    x: -self.x,
                    y: -self.y,
                    z: -self.z
                }
            }
        }
        impl Neg for &Vec3 {
            type Output = Vec3;
            fn neg(self) -> Vec3 {
                Vec3 {
                    x: -self.x,
                    y: -self.y,
                    z: -self.z
                }
            }
        }
        impl Neg for &mut Vec3 {
            type Output = Vec3;
            fn neg(self) -> Vec3 {
                Vec3 {
                    x: -self.x,
                    y: -self.y,
                    z: -self.z
                }
            }
        }
    };
}

impl_vec3_binop!(Add, add, +);
impl_vec3_binop!(Sub, sub, -);
impl_vec3_binop!(Mul, mul, *);
impl_vec3_binop!(Div, div, /);

impl_vec3_assignop!(AddAssign, add_assign, +=);
impl_vec3_assignop!(SubAssign, sub_assign, -=);
impl_vec3_assignop!(MulAssign, mul_assign, *=);
impl_vec3_assignop!(DivAssign, div_assign, /=);

impl_vec3_scalar!(Mul, mul, *);
impl_vec3_scalar!(Div, div, /);

impl_vec3_scalar_assign!(MulAssign, mul_assign, *=);
impl_vec3_scalar_assign!(DivAssign, div_assign, /=);

impl_vec3_neg!();

impl PartialEq for Vec3 {
    fn eq(&self, other: &Self) -> bool {
        self.x == other.x && self.y == other.y && self.z == other.z
    }

    fn ne(&self, other: &Self) -> bool {
        self.x != other.x || self.y != other.y || self.z != other.z
    }
}

impl PartialOrd for Vec3 {
    fn lt(&self, other: &Self) -> bool {
        self.length_square() < other.length_square()
    }

    fn gt(&self, other: &Self) -> bool {
        self.length_square() > other.length_square()
    }

    fn le(&self, other: &Self) -> bool {
        self.length_square() <= other.length_square()
    }

    fn ge(&self, other: &Self) -> bool {
        self.length_square() >= other.length_square()
    }

    fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
        self.length_square().partial_cmp(&other.length_square())
    }
}

this is mayth :)
