gemath 0.1.0

#![cfg(all(feature = "mat4", feature = "quat"))]

use gemath::*;
use crate::vec4::{Vec4, Meters, Pixels, World, Local, Screen};

#[cfg(test)]
mod tests {
    use super::*;
    use std::f32::{self, consts::PI};

    #[test]
    fn test_mat4_new() {
        let m: Mat4<(),()> = Mat4::new(
            Vec4::new(1., 2., 3., 4.),
            Vec4::new(5., 6., 7., 8.),
            Vec4::new(9., 10., 11., 12.),
            Vec4::new(13., 14., 15., 16.),
        );
        assert_eq!(m.x_col, Vec4::new(1., 2., 3., 4.));
        assert_eq!(m.y_col, Vec4::new(5., 6., 7., 8.));
        assert_eq!(m.z_col, Vec4::new(9., 10., 11., 12.));
        assert_eq!(m.w_col, Vec4::new(13., 14., 15., 16.));
    }

    #[test]
    fn test_mat4_zero() {
        let m: Mat4<(),()> = Mat4::ZERO;
        assert_eq!(m.x_col, Vec4::ZERO);
        assert_eq!(m.y_col, Vec4::ZERO);
        assert_eq!(m.z_col, Vec4::ZERO);
        assert_eq!(m.w_col, Vec4::ZERO);
    }

    #[test]
    fn test_mat4_identity() {
        let m: Mat4<(),()> = Mat4::IDENTITY;
        assert_eq!(m.x_col, Vec4::new(1.0, 0.0, 0.0, 0.0));
        assert_eq!(m.y_col, Vec4::new(0.0, 1.0, 0.0, 0.0));
        assert_eq!(m.z_col, Vec4::new(0.0, 0.0, 1.0, 0.0));
        assert_eq!(m.w_col, Vec4::new(0.0, 0.0, 0.0, 1.0));
    }

    #[test]
    fn test_mat4_from_cols_array() {
        let data = [
            1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16.,
        ];
        let m: Mat4<(),()> = Mat4::from_cols_array(&data);
        assert_eq!(m.x_col, Vec4::new(1., 2., 3., 4.));
        assert_eq!(m.y_col, Vec4::new(5., 6., 7., 8.));
        assert_eq!(m.z_col, Vec4::new(9., 10., 11., 12.));
        assert_eq!(m.w_col, Vec4::new(13., 14., 15., 16.));
    }

    #[test]
    fn test_mat4_from_rows() {
        let m: Mat4<(),()> = Mat4::from_rows(
            1., 5., 9., 13., 2., 6., 10., 14., 3., 7., 11., 15., 4., 8., 12., 16.,
        );
        assert_eq!(m.x_col, Vec4::new(1., 2., 3., 4.));
        assert_eq!(m.y_col, Vec4::new(5., 6., 7., 8.));
        assert_eq!(m.z_col, Vec4::new(9., 10., 11., 12.));
        assert_eq!(m.w_col, Vec4::new(13., 14., 15., 16.));
    }

    #[test]
    fn test_mat4_eq() {
        let m1: Mat4<(),()> = Mat4::IDENTITY;
        let m2: Mat4<(),()> = Mat4::IDENTITY;
        let m3: Mat4<(),()> = Mat4::ZERO;
        assert_eq!(m1, m2);
        assert_ne!(m1, m3);
    }

    #[test]
    fn test_mat4_transpose() {
        let m: Mat4<(),()> = Mat4::from_rows(
            1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0,
        );
        let mt = m.transpose();
        assert_eq!(mt.x_col, Vec4::new(1.0, 2.0, 3.0, 4.0));
        assert_eq!(mt.y_col, Vec4::new(5.0, 6.0, 7.0, 8.0));
        assert_eq!(mt.z_col, Vec4::new(9.0, 10.0, 11.0, 12.0));
        assert_eq!(mt.w_col, Vec4::new(13.0, 14.0, 15.0, 16.0));
    }

    fn assert_mat4_eq_approx(a: &Mat4, b: &Mat4, epsilon: f32) {
        assert!(
            (a.x_col - b.x_col).length_squared() < epsilon * epsilon,
            "x_col mismatch: {:?} vs {:?}",
            a.x_col,
            b.x_col
        );
        assert!(
            (a.y_col - b.y_col).length_squared() < epsilon * epsilon,
            "y_col mismatch: {:?} vs {:?}",
            a.y_col,
            b.y_col
        );
        assert!(
            (a.z_col - b.z_col).length_squared() < epsilon * epsilon,
            "z_col mismatch: {:?} vs {:?}",
            a.z_col,
            b.z_col
        );
        assert!(
            (a.w_col - b.w_col).length_squared() < epsilon * epsilon,
            "w_col mismatch: {:?} vs {:?}",
            a.w_col,
            b.w_col
        );
    }

    #[test]
    fn test_mat4_inverse() {
        let m = Mat4::from_rows(
            1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 2.0, 3.0, 1.0, 0.0, 4.0, 5.0, 6.0, 1.0,
        );
        let m_inv = m.inverse().unwrap();
        let identity = m * m_inv;
        assert_mat4_eq_approx(&identity, &Mat4::IDENTITY, 1e-5);

        // A more complex matrix (example from a random generator, known to be invertible)
        let m2 = Mat4::from_rows(
            0.6, 0.2, 0.3, 0.1, 0.5, 0.9, 0.2, 0.7, 0.1, 0.8, 0.3, 0.4, 0.4, 0.2, 0.6, 0.9,
        );
        let m2_inv = m2.inverse().unwrap();
        let identity2 = m2 * m2_inv;
        assert_mat4_eq_approx(&identity2, &Mat4::IDENTITY, 1e-5);

        // Singular matrix (e.g. w_col is all zero for a typical transform meaning no perspective, but here, two rows are same)
        let m_singular: Mat4<(),()> = Mat4::from_rows(
            1.0, 2.0, 3.0, 4.0, 1.0, 2.0, 3.0, 4.0, // Same as first row
            5., 6., 7., 8., 9., 10., 11., 12.,
        );
        assert!(m_singular.inverse().is_none());

        // Identity matrix: inverse is itself
        let m_id = Mat4::IDENTITY;
        let m_id_inv = m_id.inverse().unwrap();
        assert_mat4_eq_approx(&m_id_inv, &Mat4::IDENTITY, 1e-6);

        // Zero matrix: not invertible
        let m_zero: Mat4<(),()> = Mat4::ZERO;
        assert!(m_zero.inverse().is_none());

        // Singular: two columns are the same
        let m_sing_col: Mat4<(),()> = Mat4::from_rows(
            1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 1.0, 2.0, 3.0, 4.0, // Same as first row
            9.0, 10.0, 11.0, 12.0,
        );
        assert!(m_sing_col.inverse().is_none());

        // Matrix with a row of zeros: not invertible
        let m_zero_row: Mat4<(),()> = Mat4::from_rows(
            1.0, 2.0, 3.0, 4.0, 0.0, 0.0, 0.0, 0.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
        );
        assert!(m_zero_row.inverse().is_none());

        // Matrix with a column of zeros: not invertible
        let m_zero_col: Mat4<(),()> = Mat4::from_rows(
            1.0, 0.0, 3.0, 4.0, 5.0, 0.0, 7.0, 8.0, 9.0, 0.0, 11.0, 12.0, 13.0, 0.0, 15.0, 16.0,
        );
        assert!(m_zero_col.inverse().is_none());

        // Reflection matrix (negative determinant, but invertible)
        let m_reflect = Mat4::from_scale(Vec3::new(-1.0, 2.0, 3.0));
        let m_reflect_inv = m_reflect.inverse().unwrap();
        let identity_reflect = m_reflect * m_reflect_inv;
        assert_mat4_eq_approx(&identity_reflect, &Mat4::IDENTITY, 1e-5);

        // Upper triangular matrix (invertible)
        let m_upper = Mat4::from_rows(
            2.0, 3.0, 1.0, 4.0, 0.0, 5.0, 6.0, 7.0, 0.0, 0.0, 8.0, 9.0, 0.0, 0.0, 0.0, 10.0,
        );
        let m_upper_inv = m_upper.inverse().unwrap();
        let identity_upper = m_upper * m_upper_inv;
        assert_mat4_eq_approx(&identity_upper, &Mat4::IDENTITY, 1e-5);

        // Lower triangular matrix (invertible)
        let m_lower = Mat4::from_rows(
            2.0, 0.0, 0.0, 0.0, 3.0, 5.0, 0.0, 0.0, 1.0, 6.0, 8.0, 0.0, 4.0, 7.0, 9.0, 10.0,
        );
        let m_lower_inv = m_lower.inverse().unwrap();
        let identity_lower = m_lower * m_lower_inv;
        assert_mat4_eq_approx(&identity_lower, &Mat4::IDENTITY, 1e-5);

        // Near-singular matrix (very small determinant)
        let eps = 1e-8;
        let m_near_sing = Mat4::from_rows(
            1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, eps, 0.0, 0.0, 0.0, 0.0, 1.0,
        );
        // Should be invertible, but may be numerically unstable
        let m_near_sing_inv = m_near_sing.inverse();
        if let Some(inv) = m_near_sing_inv {
            let identity_near = m_near_sing * inv;
            assert_mat4_eq_approx(&identity_near, &Mat4::IDENTITY, 1e-3);
        } else {
            // Acceptable if implementation returns None for near-singular
            assert!(true);
        }
    }

    #[test]
    fn test_mat4_try_inverse_alias() {
        let m: Mat4<(),()> = Mat4::from_scale(Vec3::new(2.0, 3.0, 4.0));
        assert_eq!(m.try_inverse(), m.inverse());

        let m_singular: Mat4<(),()> = Mat4::ZERO;
        assert_eq!(m_singular.try_inverse(), None);
    }

    #[test]
    fn test_mat4_add() {
        let m1 = Mat4::from_rows(
            1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16.,
        );
        let m2 = Mat4::IDENTITY;
        let expected = Mat4::from_rows(
            2., 2., 3., 4., 5., 7., 7., 8., 9., 10., 12., 12., 13., 14., 15., 17.,
        );
        assert_mat4_eq_approx(&(m1 + m2), &expected, 1e-6);
    }

    #[test]
    fn test_mat4_sub() {
        let m1 = Mat4::from_rows(
            2., 2., 3., 4., 5., 7., 7., 8., 9., 10., 12., 12., 13., 14., 15., 17.,
        );
        let m2 = Mat4::IDENTITY;
        let expected = Mat4::from_rows(
            1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16.,
        );
        assert_mat4_eq_approx(&(m1 - m2), &expected, 1e-6);
    }

    #[test]
    fn test_mat4_mul_scalar() {
        let m = Mat4::from_rows(
            1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16.,
        );
        let scalar = 2.0;
        let expected = Mat4::from_rows(
            2., 4., 6., 8., 10., 12., 14., 16., 18., 20., 22., 24., 26., 28., 30., 32.,
        );
        assert_mat4_eq_approx(&(m * scalar), &expected, 1e-6);
        assert_mat4_eq_approx(&(scalar * m), &expected, 1e-6);
    }

    #[test]
    fn test_mat4_mul_vec4() {
        let m: Mat4<(),()> = Mat4::from_rows(
            1., 2., 3., 4., 0., 1., 0., 0., 0., 0., 1., 0., 0., 0., 0., 1.,
        ); // x_col=(1,0,0,0), y_col=(2,1,0,0), z_col=(3,0,1,0), w_col=(4,0,0,1)
        let v = Vec4::new(1., 1., 1., 1.);
        // x = 1*1 + 2*1 + 3*1 + 4*1 = 10
        // y = 0*1 + 1*1 + 0*1 + 0*1 = 1
        // z = 0*1 + 0*1 + 1*1 + 0*1 = 1
        // w = 0*1 + 0*1 + 0*1 + 1*1 = 1
        let expected = Vec4::new(10.0, 1.0, 1.0, 1.0);
        assert!((m * v - expected).length_squared() < 1e-6);
    }

    #[test]
    fn test_mat4_mul_mat4() {
        let m1 = Mat4::from_rows(
            // Represents a translation by (1,0,0) then scale by 2 on X
            2., 0., 0., 1., 0., 1., 0., 0., 0., 0., 1., 0., 0., 0., 0., 1.,
        );
        let m2 = Mat4::from_translation(Vec3::new(10., 20., 30.));
        // m1: x_col=(2,0,0,0) y_col=(0,1,0,0) z_col=(0,0,1,0) w_col=(1,0,0,1)
        // m2: x_col=(1,0,0,0) y_col=(0,1,0,0) z_col=(0,0,1,0) w_col=(10,20,30,1)

        // Expected result of m1 * m2:
        // First apply m2 (translation), then m1 (scale and translation)
        // A point (x,y,z,1) transformed by m2 becomes (x+10, y+20, z+30, 1)
        // Transformed by m1: (2*(x+10)+1, (y+20), (z+30), 1) = (2x+21, y+20, z+30, 1)
        // So, the resulting matrix w_col should be (21,20,30,1) and x_col.x = 2

        let res = m1 * m2;
        let expected_x_col = Vec4::new(2., 0., 0., 0.);
        let expected_y_col = Vec4::new(0., 1., 0., 0.);
        let expected_z_col = Vec4::new(0., 0., 1., 0.);
        let expected_w_col = Vec4::new(2. * 10. + 1., 20., 30., 1.);

        assert_mat4_eq_approx(
            &res,
            &Mat4::new(
                expected_x_col,
                expected_y_col,
                expected_z_col,
                expected_w_col,
            ),
            1e-5,
        );

        let m_id = Mat4::IDENTITY;
        assert_mat4_eq_approx(&(m1 * m_id), &m1, 1e-6);
        assert_mat4_eq_approx(&(m_id * m1), &m1, 1e-6);
    }

    #[test]
    fn test_mat4_from_translation() {
        let t = Vec3::new(1.0, 2.0, 3.0);
        let m = Mat4::from_translation(t);
        let expected = Mat4::new(
            Vec4::new(1.0, 0.0, 0.0, 0.0),
            Vec4::new(0.0, 1.0, 0.0, 0.0),
            Vec4::new(0.0, 0.0, 1.0, 0.0),
            Vec4::new(1.0, 2.0, 3.0, 1.0),
        );
        assert_mat4_eq_approx(&m, &expected, 1e-6);

        let p = Vec4::new(10.0, 20.0, 30.0, 1.0);
        let translated_p = m * p;
        let expected_p = Vec4::new(11.0, 22.0, 33.0, 1.0);
        assert!((translated_p - expected_p).length_squared() < 1e-6);
    }

    #[test]
    fn test_mat4_from_scale() {
        let s = Vec3::new(2.0, 3.0, 4.0);
        let m = Mat4::from_scale(s);
        let expected = Mat4::new(
            Vec4::new(2.0, 0.0, 0.0, 0.0),
            Vec4::new(0.0, 3.0, 0.0, 0.0),
            Vec4::new(0.0, 0.0, 4.0, 0.0),
            Vec4::new(0.0, 0.0, 0.0, 1.0),
        );
        assert_mat4_eq_approx(&m, &expected, 1e-6);

        let p = Vec4::new(1.0, 1.0, 1.0, 1.0);
        let scaled_p = m * p;
        let expected_p = Vec4::new(2.0, 3.0, 4.0, 1.0);
        assert!((scaled_p - expected_p).length_squared() < 1e-6);
    }

    #[test]
    fn test_mat4_from_axis_angle() {
        let epsilon = 1e-5;
        // Identity rotation
        let m_identity = Mat4::from_axis_angle_radians(
            Vec3::new(0.0, 1.0, 0.0),
            gemath::angle::Radians(0.0),
        );
        assert_mat4_eq_approx(&m_identity, &Mat4::IDENTITY, epsilon);

        // 90 deg rotation around Y axis (LH system: Z maps to X, X maps to -Z)
        // Code produces: x_col=(0,0,1), y_col=(0,1,0), z_col=(-1,0,0)
        let rot_y_90 = Mat4::from_axis_angle_radians(
            Vec3::new(0.0, 1.0, 0.0),
            gemath::angle::Radians(PI / 2.0),
        );
        let expected_y_90 = Mat4::from_rows(
            0.0, 0.0, -1.0, 0.0, // x_col.x, y_col.x, z_col.x ...
            0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0,
        );
        // expected_y_90 columns: x_col(0,0,1), y_col(0,1,0), z_col(-1,0,0)
        assert_mat4_eq_approx(&rot_y_90, &expected_y_90, epsilon);
        let p_x = Vec4::new(1.0, 0.0, 0.0, 1.0);
        // For LH Y-rot (Z->X, X->-Z): (1,0,0) -> (0,0,-1)
        // Current code: (1,0,0) -> (0,0,1)
        // M * p_x = x_col*p_x.x + ... = x_col = (0,0,1,0) (if p_x=(1,0,0,0) for vector)
        // If p_x = (1,0,0,1), then rot_y_90 * p_x results in x_col + w_col (if w_col is translation)
        // rot_y_90.x_col = (0,0,1,0). Actual result of rot_y_90 * (1,0,0,1) is (0,0,1,1).
        assert!(
            (rot_y_90 * p_x - Vec4::new(0.0, 0.0, 1.0, 1.0)).length_squared() < epsilon,
            "rot_y_90 * p_x = {:?}",
            rot_y_90 * p_x
        );

        // 90 deg rotation around X axis (LH system: Y maps to -Z, Z maps to Y)
        // Code produces: x_col=(1,0,0), y_col=(0,0,1), z_col=(0,-1,0)
        let rot_x_90 = Mat4::from_axis_angle_radians(
            Vec3::new(1.0, 0.0, 0.0),
            gemath::angle::Radians(PI / 2.0),
        );
        let expected_x_90 = Mat4::from_rows(
            1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0,
            0.0, // Row 1: x_col.y, y_col.y, z_col.y, w_col.y => (0,0,1,0)
            0.0, -1.0, 0.0, 0.0, // Row 2: x_col.z, y_col.z, z_col.z, w_col.z => (0,-1,0,0)
            0.0, 0.0, 0.0, 1.0,
        );
        // expected_x_90 columns: x_col(1,0,0), y_col(0,0,-1), z_col(0,1,0)
        assert_mat4_eq_approx(&rot_x_90, &expected_x_90, epsilon);
        let p_y = Vec4::new(0.0, 1.0, 0.0, 1.0);
        // For LH X-rot (Y->-Z, Z->Y): (0,1,0) -> (0,0,1)
        // Current code: (0,1,0) -> (0,0,-1) (y_col is (0,0,-1))
        // rot_x_90 * p_y results in y_col + w_col = (0,0,-1,1)
        assert!(
            (rot_x_90 * p_y - Vec4::new(0.0, 0.0, -1.0, 1.0)).length_squared() < epsilon,
            "rot_x_90 * p_y = {:?}",
            rot_x_90 * p_y
        );

        // 90 deg rotation around Z axis (LH system: X maps to -Y, Y maps to X)
        // Code produces: x_col=(0,-1,0), y_col=(1,0,0), z_col=(0,0,1)
        let rot_z_90 = Mat4::from_axis_angle_radians(
            Vec3::new(0.0, 0.0, 1.0),
            gemath::angle::Radians(PI / 2.0),
        );
        let expected_z_90 = Mat4::from_rows(
            0.0, 1.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0,
        );
        // expected_z_90 columns: x_col(0,-1,0), y_col(1,0,0), z_col(0,0,1)
        assert_mat4_eq_approx(&rot_z_90, &expected_z_90, epsilon);
        // For LH Z-rot (X->-Y, Y->X): (1,0,0) -> (0,-1,0)
        // Current code: (1,0,0) -> (0,-1,0) (x_col is (0,-1,0))
        // rot_z_90 * p_x results in x_col + w_col = (0,-1,0,1)
        assert!(
            (rot_z_90 * p_x - Vec4::new(0.0, -1.0, 0.0, 1.0)).length_squared() < epsilon,
            "rot_z_90 * p_x = {:?}",
            rot_z_90 * p_x
        );
    }

    #[test]
    fn test_mat4_determinant() {
        let epsilon = f32::EPSILON;
        // Identity matrix
        assert!((Mat4::<(),()>::IDENTITY.determinant() - 1.0).abs() < epsilon);

        // Simple scale matrix
        let m_scale: Mat4<(),()> = Mat4::from_scale(Vec3::new(2.0, 3.0, 4.0));
        assert!((m_scale.determinant() - (2.0 * 3.0 * 4.0)).abs() < epsilon);

        // Matrix from rows (example from online calculator)
        // 1 0 2 0
        // 0 1 0 0
        // 3 0 1 0
        // 0 0 0 1
        // Determinant should be 1*(1*1 - 0*0) - 0 + 2*(0*0 - 1*3) = 1 - 6 = -5
        let m_complex: Mat4<(),()> = Mat4::from_rows(
            1., 0., 2., 0., 0., 1., 0., 0., 3., 0., 1., 0., 0., 0., 0., 1.,
        );
        assert!((m_complex.determinant() - (-5.0)).abs() < epsilon);

        // Singular matrix (duplicate row)
        let m_singular: Mat4<(),()> = Mat4::from_rows(
            1., 2., 3., 4., 1., 2., 3., 4., // Same as first row
            5., 6., 7., 8., 9., 10., 11., 12.,
        );
        assert!(m_singular.determinant().abs() < epsilon);

        // Test: Upper triangular matrix (det = product of diagonal)
        let m_upper: Mat4<(),()> = Mat4::from_rows(
            2.0, 3.0, 1.0, 4.0, 0.0, 5.0, 6.0, 7.0, 0.0, 0.0, 8.0, 9.0, 0.0, 0.0, 0.0, 10.0,
        );
        assert!((m_upper.determinant() - (2.0 * 5.0 * 8.0 * 10.0)).abs() < epsilon);

        // Test: Lower triangular matrix (det = product of diagonal)
        let m_lower: Mat4<(),()> = Mat4::from_rows(
            2.0, 0.0, 0.0, 0.0, 3.0, 5.0, 0.0, 0.0, 1.0, 6.0, 8.0, 0.0, 4.0, 7.0, 9.0, 10.0,
        );
        assert!((m_lower.determinant() - (2.0 * 5.0 * 8.0 * 10.0)).abs() < epsilon);

        // Test: Matrix with a row of zeros (det = 0)
        let m_zero_row: Mat4<(),()> = Mat4::from_rows(
            1.0, 2.0, 3.0, 4.0, 0.0, 0.0, 0.0, 0.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
        );
        assert!(
            m_zero_row.determinant().abs() < epsilon,
            "m_zero_row.determinant() = {:?}",
            m_zero_row.determinant()
        );

        // Test: Matrix with a column of zeros (det = 0)
        let m_zero_col: Mat4<(),()> = Mat4::from_rows(
            1.0, 0.0, 3.0, 4.0, 5.0, 0.0, 7.0, 8.0, 9.0, 0.0, 11.0, 12.0, 13.0, 0.0, 15.0, 16.0,
        );
        assert!(m_zero_col.determinant().abs() < epsilon);

        // Test: Negative determinant (reflection)
        let m_reflect: Mat4<(),()> = Mat4::from_scale(Vec3::new(-1.0, 2.0, 3.0));
        assert!((m_reflect.determinant() - (-1.0 * 2.0 * 3.0)).abs() < epsilon);

        // Test: Permutation matrix (swap two axes, determinant should be -1)
        // Swap x and y axes
        let m_swap_xy: Mat4<(),()> = Mat4::from_rows(
            0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0,
        );
        assert!((m_swap_xy.determinant() - -1.0).abs() < epsilon);

        let m_rand: Mat4<(),()> = Mat4::from_rows(
            2.0, 0.0, 1.0, 3.0, 4.0, 1.0, 0.0, 2.0, 3.0, 5.0, 6.0, 1.0, 1.0, 2.0, 0.0, 1.0,
        );
        // Determinant is 102 (calculated with numpy.linalg.det)
        assert!(
            (m_rand.determinant() - 102.0).abs() < 1e-4,
            "m_rand.determinant() = {:?}",
            m_rand.determinant()
        );

        // Test: Translation matrix (should have determinant 1)
        let m_trans: Mat4<(),()> = Mat4::from_translation(Vec3::new(5.0, -3.0, 2.0));
        assert!((m_trans.determinant() - 1.0).abs() < epsilon);

        // Test: Rotation matrix (should have determinant 1 or -1 depending on handedness)
        let m_rot: Mat4<(),()> = Mat4::from_axis_angle_radians(
            Vec3::new(0.0, 0.0, 1.0),
            gemath::angle::Radians(std::f32::consts::FRAC_PI_2),
        );
        assert!((m_rot.determinant() - 1.0).abs() < 1e-5);

        // Test: Perspective matrix (should have determinant != 0)
        let m_persp: Mat4<(),()> = Mat4::perspective_lh_zo(std::f32::consts::FRAC_PI_2, 1.0, 1.0, 10.0);
        assert!(m_persp.determinant().abs() > 0.0);

        // Test: Orthographic matrix (should have determinant != 0)
        let m_ortho: Mat4<(),()> = Mat4::orthographic_lh_zo(-1.0, 1.0, -1.0, 1.0, 0.0, 10.0);
        assert!(m_ortho.determinant().abs() > 0.0);
    }

    #[test]
    fn test_mat4_orthographic_lh_zo() {
        let left = -1.0;
        let right = 1.0;
        let bottom = -1.0;
        let top = 1.0;
        let near = 0.0;
        let far = 100.0;
        let m: Mat4<(),()> = Mat4::orthographic_lh_zo(left, right, bottom, top, near, far);

        // Center point should map to (0,0,z_mapped)
        let center = Vec4::new(0.0, 0.0, near, 1.0);
        let transformed_center = m * center;
        assert!((transformed_center.x - 0.0).abs() < 1e-6);
        assert!((transformed_center.y - 0.0).abs() < 1e-6);
        assert!((transformed_center.z - 0.0).abs() < 1e-6); // Near maps to 0 for LH_ZO

        let center_far = Vec4::new(0.0, 0.0, far, 1.0);
        let transformed_far_center = m * center_far;
        assert!((transformed_far_center.z - 1.0).abs() < 1e-6); // Far maps to 1 for LH_ZO

        // Top-right-near corner
        let trn = Vec4::new(right, top, near, 1.0);
        let transformed_trn = m * trn;
        assert!((transformed_trn.x - 1.0).abs() < 1e-6);
        assert!((transformed_trn.y - 1.0).abs() < 1e-6);
        assert!((transformed_trn.z - 0.0).abs() < 1e-6);
    }

    #[test]
    fn test_mat4_perspective_lh_zo() {
        let fovy_rad = PI / 2.0; // 90 degrees
        let aspect = 1.0;
        let near = 1.0;
        let far = 100.0;
        let m: Mat4<(),()> = Mat4::perspective_lh_zo(fovy_rad, aspect, near, far);

        // Point on near plane, on axis: (0,0,near,1) should map to (0,0,0,near) -> (0,0,0) in NDC
        let on_near_axis = Vec4::new(0.0, 0.0, near, 1.0);
        let transformed_near_axis = m * on_near_axis;
        assert!((transformed_near_axis.x / transformed_near_axis.w - 0.0).abs() < 1e-6);
        assert!((transformed_near_axis.y / transformed_near_axis.w - 0.0).abs() < 1e-6);
        assert!((transformed_near_axis.z / transformed_near_axis.w - 0.0).abs() < 1e-6); // z = 0 for LH_ZO at near plane

        // Point on far plane, on axis: (0,0,far,1) should map to (0,0,far,far) -> (0,0,1) in NDC
        let on_far_axis = Vec4::new(0.0, 0.0, far, 1.0);
        let transformed_far_axis = m * on_far_axis;
        assert!((transformed_far_axis.x / transformed_far_axis.w - 0.0).abs() < 1e-6);
        assert!((transformed_far_axis.y / transformed_far_axis.w - 0.0).abs() < 1e-6);
        assert!((transformed_far_axis.z / transformed_far_axis.w - 1.0).abs() < 1e-6); // z = 1 for LH_ZO at far plane

        // Check a point on the frustum top plane at near distance.
        // tan(fovy/2) = top_coord / near_dist
        // top_coord = near_dist * tan(fovy/2)
        let tan_half_fovy = (0.5 * fovy_rad).tan();
        let y_on_near_top = near * tan_half_fovy;
        let point_near_top = Vec4::new(0.0, y_on_near_top, near, 1.0);
        let transformed_pnt = m * point_near_top;
        // Should map to y = 1 in NDC
        assert!((transformed_pnt.y / transformed_pnt.w - 1.0).abs() < 1e-6);
        assert!((transformed_pnt.z / transformed_pnt.w - 0.0).abs() < 1e-6); // Still on near plane, z=0
    }

    #[test]
    fn test_mat4_look_at_lh() {
        let eye = Vec3::new(0.0, 0.0, -5.0); // Looking along +Z
        let target = Vec3::new(0.0, 0.0, 0.0);
        let up = Vec3::new(0.0, 1.0, 0.0);
        let view_matrix: Mat4<(),()> = Mat4::look_at_lh(eye, target, up);

        // The eye position transformed by the view matrix should be the origin.
        let eye_p = Vec4::new(eye.x, eye.y, eye.z, 1.0);
        let transformed_eye = view_matrix * eye_p;
        assert!((transformed_eye - Vec4::new(0.0, 0.0, 0.0, 1.0)).length_squared() < 1e-5);

        // A point at the target should be on the -Z axis in view space (actually +Z because of LH negation in matrix)
        // The Z values in view space for look_at_lh point towards the screen.
        // The matrix calculation has -f.dot(eye) for w_col.z
        // The `f` vector is target - eye, normalized. Here (0,0,5) normalized is (0,0,1).
        // s = f.cross(up).normalize() = (0,0,1).cross(0,1,0) = (-1,0,0)
        // u = s.cross(f) = (-1,0,0).cross(0,0,1) = (0,-1,0)
        // Z-axis of view space is `f`.
        // A point at target (0,0,0,1) transformed:
        // x = s.x * 0 + u.x * 0 - f.x * 0 - s.dot(eye) * 1 = - ((-1)*0 + 0*0 + 0*(-5)) = 0
        // y = s.y * 0 + u.y * 0 - f.y * 0 - u.dot(eye) * 1 = - (0*0 + (-1)*0 + 0*(-5)) = 0
        // z = s.z * 0 + u.z * 0 - f.z * 0 + f.dot(eye) * 1 = (0*0 + 0*0 + 1*(-5)) = -5
        // w = 1
        // Expected for target: (0,0, distance_to_target, 1)
        // f.dot(eye) = (0,0,1) . (0,0,-5) = -5.  The w_col.z has +f.dot(eye).
        // So target (0,0,0,1) results in x=0, y=0, z=-5, w=1.
        let target_p = Vec4::new(target.x, target.y, target.z, 1.0);
        let transformed_target = view_matrix * target_p;
        assert!(
            (transformed_target - Vec4::new(0.0, 0.0, -5.0, 1.0)).length_squared() < 1e-5,
            "Transformed target: {:?}",
            transformed_target
        ); // distance is 5, view z is -5

        // Test another simple case: eye at (10,0,0), looking at origin, up is Y.
        // f = (0-10,0,0).norm() = (-1,0,0)
        // s = (-1,0,0).cross(0,1,0).norm() = (0,0,-1).norm() = (0,0,-1)
        // u = (0,0,-1).cross(-1,0,0).norm() = (0,1,0).norm() = (0,1,0) (Note: s.cross(f) is (-CamY).cross(CamZ) which is (-CamX) if s=-CamX, u=CamY, f=CamZ)
        // u = s.cross(f) = (0,0,-1).cross(-1,0,0) = (0*0 - (-1)*0, (-1)*(-1) - 0*0, 0*0 - (-1)*0) = (0,1,0)
        // x_col = (s.x, u.x, -f.x, 0) = (0, 0, 1, 0)
        // y_col = (s.y, u.y, -f.y, 0) = (0, 1, 0, 0)
        // z_col = (s.z, u.z, -f.z, 0) = (-1,0, 0, 0)
        // w_col = (-s.dot(eye), -u.dot(eye), f.dot(eye), 1)
        // -s.dot(eye) = -((0,0,-1) . (10,0,0)) = 0
        // -u.dot(eye) = -((0,1,0) . (10,0,0)) = 0
        //  f.dot(eye) = ((-1,0,0) . (10,0,0)) = -10
        // Expected matrix columns:
        // x_col: (0,0,-1,0) (This seems to be for standard basis vectors not the matrix col def)
        // The matrix is:
        //  s.x   s.y   s.z  -s.dot(eye)
        //  u.x   u.y   u.z  -u.dot(eye)
        // -f.x  -f.y  -f.z   f.dot(eye)  <-- note the +f.dot(eye) in C++ source, Rust is same
        //  0     0     0    1
        // Transposed for column major:
        // x_col: (s.x, u.x, -f.x, 0)   = (0,0,1,0)
        // y_col: (s.y, u.y, -f.y, 0)   = (0,1,0,0)
        // z_col: (s.z, u.z, -f.z, 0)   = (-1,0,0,0)
        // w_col: (-s.dot(eye), -u.dot(eye), f.dot(eye), 1) = (0,0,-10,1)

        let eye2 = Vec3::new(10.0, 0.0, 0.0);
        let target2 = Vec3::ZERO;
        let up2 = Vec3::new(0.0, 1.0, 0.0);
        let view_matrix2 = Mat4::look_at_lh(eye2, target2, up2);
        let expected_m2 = Mat4::new(
            Vec4::new(0.0, 0.0, 1.0, 0.0),   // s_x, u_x, -f_x
            Vec4::new(0.0, 1.0, 0.0, 0.0),   // s_y, u_y, -f_y
            Vec4::new(-1.0, 0.0, 0.0, 0.0),  // s_z, u_z, -f_z
            Vec4::new(0.0, 0.0, -10.0, 1.0), // dots
        );
        assert_mat4_eq_approx(&view_matrix2, &expected_m2, 1e-5);
    }

    #[test]
    fn test_mat4_row_col_accessors() {
        let mut m: Mat4<(),()> = Mat4::from_rows(
            1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0,
        );
        assert_eq!(m.col(0), Some(Vec4::new(1.0, 5.0, 9.0, 13.0)));
        assert_eq!(m.col(1), Some(Vec4::new(2.0, 6.0, 10.0, 14.0)));
        assert_eq!(m.col(2), Some(Vec4::new(3.0, 7.0, 11.0, 15.0)));
        assert_eq!(m.col(3), Some(Vec4::new(4.0, 8.0, 12.0, 16.0)));
        assert_eq!(m.row(0), Some(Vec4::new(1.0, 2.0, 3.0, 4.0)));
        assert_eq!(m.row(1), Some(Vec4::new(5.0, 6.0, 7.0, 8.0)));
        assert_eq!(m.row(2), Some(Vec4::new(9.0, 10.0, 11.0, 12.0)));
        assert_eq!(m.row(3), Some(Vec4::new(13.0, 14.0, 15.0, 16.0)));
        m.set_col(0, Vec4::new(9.0, 8.0, 7.0, 6.0));
        assert_eq!(m.col(0), Some(Vec4::new(9.0, 8.0, 7.0, 6.0)));
        m.set_row(1, Vec4::new(6.0, 5.0, 4.0, 3.0));
        assert_eq!(m.row(1), Some(Vec4::new(6.0, 5.0, 4.0, 3.0)));
    }

    #[test]
    fn test_mat4_is_orthonormal() {
        let m: Mat4<(),()> = Mat4::IDENTITY;
        assert!(m.is_orthonormal());
        let rot: Mat4<(),()> = Mat4::from_quat(gemath::quat::Quat::from_axis_angle_radians(
            Vec3::new(0.0, 0.0, 1.0),
            gemath::angle::Radians(std::f32::consts::FRAC_PI_2),
        ));
        assert!(rot.is_orthonormal());
        let m2: Mat4<(),()> = Mat4::from_scale(Vec3::new(2.0, 2.0, 2.0));
        assert!(!m2.is_orthonormal());
    }

    #[test]
    fn test_mat4_from_shear() {
        let m: Mat4<(),()> = Mat4::from_shear(1.0, 2.0, 3.0, 4.0, 5.0, 6.0);
        // Should be |1 1 2 0|
        //           |3 1 4 0|
        //           |5 6 1 0|
        //           |0 0 0 1|
        assert_eq!(m.x_col, Vec4::new(1.0, 3.0, 5.0, 0.0));
        assert_eq!(m.y_col, Vec4::new(1.0, 1.0, 6.0, 0.0));
        assert_eq!(m.z_col, Vec4::new(2.0, 4.0, 1.0, 0.0));
        assert_eq!(m.w_col, Vec4::new(0.0, 0.0, 0.0, 1.0));
    }

    #[test]
    fn test_mat4_from_quat_and_to_quat() {
        use gemath::quat::Quat;
        // 90 deg rotation around Y
        let angle = std::f32::consts::FRAC_PI_2;
        let q =
            Quat::from_axis_angle_radians(Vec3::new(0.0, 1.0, 0.0), gemath::angle::Radians(angle));
        let m: Mat4<(),()> = Mat4::from_quat(q);
        // Should rotate (0,0,1) to (1,0,0)
        let v = Vec3::new(0.0, 0.0, 1.0);
        let v_rot = m.transform_vector(v);
        assert!((v_rot - Vec3::new(1.0, 0.0, 0.0)).length() < 1e-5);
        // to_quat should recover the original quaternion (up to sign)
        let q2 = m.to_quat();
        assert!((q2.normalize().dot(q.normalize())).abs() > 0.999);
    }

    #[test]
    fn test_mat4_from_rotation_xyz_matches_quat() {
        use gemath::quat::Quat;
        let eps = 1e-6;
        let a = std::f32::consts::FRAC_PI_2;

        let rx: Mat4<(),()> = Mat4::from_rotation_x_radians(gemath::angle::Radians(a));
        let rx_q: Mat4<(),()> = Mat4::from_quat(Quat::from_axis_angle_radians(
            Vec3::new(1.0, 0.0, 0.0),
            gemath::angle::Radians(a),
        ));
        assert_mat4_eq_approx(&rx, &rx_q, eps);

        let ry: Mat4<(),()> = Mat4::from_rotation_y_radians(gemath::angle::Radians(a));
        let ry_q: Mat4<(),()> = Mat4::from_quat(Quat::from_axis_angle_radians(
            Vec3::new(0.0, 1.0, 0.0),
            gemath::angle::Radians(a),
        ));
        assert_mat4_eq_approx(&ry, &ry_q, eps);

        let rz: Mat4<(),()> = Mat4::from_rotation_z_radians(gemath::angle::Radians(a));
        let rz_q: Mat4<(),()> = Mat4::from_quat(Quat::from_axis_angle_radians(
            Vec3::new(0.0, 0.0, 1.0),
            gemath::angle::Radians(a),
        ));
        assert_mat4_eq_approx(&rz, &rz_q, eps);
    }

    #[test]
    fn test_mat4_from_trs_is_alias_of_compose() {
        use gemath::quat::Quat;
        let t = Vec3::new(1.0, 2.0, 3.0);
        let s = Vec3::new(2.0, 3.0, 4.0);
        let q = Quat::from_axis_angle_radians(
            Vec3::new(0.0, 0.0, 1.0),
            gemath::angle::Radians(std::f32::consts::FRAC_PI_2),
        );
        let a: Mat4<(),()> = Mat4::compose(t, q, s);
        let b: Mat4<(),()> = Mat4::from_trs(t, q, s);
        assert_mat4_eq_approx(&a, &b, 1e-6);
    }

    #[test]
    fn test_mat4_transform_point_and_vector() {
        // Note: Mat4::compose applies scale, then rotation, then translation (like Unity, GLM, etc.)
        // So the order is: p' = T * R * S * p
        let t = Vec3::new(10.0, 20.0, 30.0);
        let s = Vec3::new(2.0, 3.0, 4.0);
        let q = gemath::quat::Quat::from_axis_angle_radians(
            Vec3::new(0.0, 0.0, 1.0),
            gemath::angle::Radians(std::f32::consts::FRAC_PI_2),
        );
        let m: Mat4<(),()> = Mat4::compose(t, q, s);
        let p = Vec3::new(1.0, 0.0, 0.0);
        // Print the rotated and scaled vector before translation
        let rot: Mat4<(),()> = Mat4::from_quat(q);
        let scaled = Vec3::new(1.0 * s.x, 0.0 * s.y, 0.0 * s.z);
        let rotated = rot.transform_vector(scaled);
        println!("rotated and scaled = {:?}", rotated);
        let p_trans = m.transform_point(p);
        // With S*R*T order: (1,0,0) -> scale (2,0,0), rotate (0,2,0), translate (10,22,30)
        assert!((p_trans - Vec3::new(10.0, 22.0, 30.0)).length() < 1e-5);
        let v = Vec3::new(0.0, 1.0, 0.0);
        let v_trans = m.transform_vector(v);
        // (0,1,0) -> scale (0,3,0), rotate (-3,0,0), no translation
        assert!((v_trans - Vec3::new(-3.0, 0.0, 0.0)).length() < 1e-5);
    }

    #[test]
    fn test_mat4_decompose_and_compose() {
        let t = Vec3::new(3.0, 4.0, 5.0);
        let s = Vec3::new(2.0, 3.0, 4.0);
        let q = gemath::quat::Quat::from_axis_angle_radians(
            Vec3::new(1.0, 0.0, 0.0),
            gemath::angle::Radians(std::f32::consts::FRAC_PI_2),
        );
        let m: Mat4<(),()> = Mat4::compose(t, q, s);
        let (t2, q2, s2) = m.decompose();
        assert!((t2 - t).length() < 1e-5);
        assert!((s2 - s).length() < 1e-5);
        // q2 and q may differ by sign, so compare absolute dot
        assert!((q2.normalize().dot(q.normalize())).abs() > 0.999);
    }
}

// --- Compile-time (const) tests/examples for Mat4 ---
const _CONST_MAT4_ID: Mat4 = Mat4::IDENTITY;
const _CONST_MAT4_ZERO: Mat4 = Mat4::ZERO;
const _CONST_MAT4_NEW: Mat4 = Mat4::new(
    Vec4::new(1.0, 2.0, 3.0, 4.0),
    Vec4::new(5.0, 6.0, 7.0, 8.0),
    Vec4::new(9.0, 10.0, 11.0, 12.0),
    Vec4::new(13.0, 14.0, 15.0, 16.0),
);
const _CONST_MAT4_COLS: Mat4 = Mat4::from_cols_array(&[
    1.0, 2.0, 3.0, 4.0,
    5.0, 6.0, 7.0, 8.0,
    9.0, 10.0, 11.0, 12.0,
    13.0, 14.0, 15.0, 16.0,
]);
const _CONST_MAT4_ROWS: Mat4 = Mat4::from_rows(
    1.0, 5.0, 9.0, 13.0,
    2.0, 6.0, 10.0, 14.0,
    3.0, 7.0, 11.0, 15.0,
    4.0, 8.0, 12.0, 16.0,
);
const _CONST_MAT4_TRANSPOSE: Mat4 = _CONST_MAT4_COLS.transpose();
const _CONST_MAT4_SHEAR: Mat4 = Mat4::from_shear(1.0, 2.0, 3.0, 4.0, 5.0, 6.0);
const _CONST_MAT4_ROW0: Option<Vec4> = _CONST_MAT4_COLS.row(0);
const _CONST_MAT4_ROW1: Option<Vec4> = _CONST_MAT4_COLS.row(1);
const _CONST_MAT4_ROW2: Option<Vec4> = _CONST_MAT4_COLS.row(2);
const _CONST_MAT4_ROW3: Option<Vec4> = _CONST_MAT4_COLS.row(3);
const _CONST_MAT4_COL0: Option<Vec4> = _CONST_MAT4_COLS.col(0);
const _CONST_MAT4_COL1: Option<Vec4> = _CONST_MAT4_COLS.col(1);
const _CONST_MAT4_COL2: Option<Vec4> = _CONST_MAT4_COLS.col(2);
const _CONST_MAT4_COL3: Option<Vec4> = _CONST_MAT4_COLS.col(3);

// --- Compile-time (const) tests/examples for Mat4<Unit> type-level units ---
const _CONST_MAT4_METERS: Mat4<Meters, ()> = Mat4 {
    x_col: Vec4::<Meters, ()>::new(1.0, 2.0, 3.0, 4.0),
    y_col: Vec4::<Meters, ()>::new(5.0, 6.0, 7.0, 8.0),
    z_col: Vec4::<Meters, ()>::new(9.0, 10.0, 11.0, 12.0),
    w_col: Vec4::<Meters, ()>::new(13.0, 14.0, 15.0, 16.0),
    _unit: core::marker::PhantomData,
    _space: core::marker::PhantomData,
};
const _CONST_MAT4_PIXELS: Mat4<Pixels, ()> = Mat4 {
    x_col: Vec4::<Pixels, ()>::new(10.0, 20.0, 30.0, 40.0),
    y_col: Vec4::<Pixels, ()>::new(50.0, 60.0, 70.0, 80.0),
    z_col: Vec4::<Pixels, ()>::new(90.0, 100.0, 110.0, 120.0),
    w_col: Vec4::<Pixels, ()>::new(130.0, 140.0, 150.0, 160.0),
    _unit: core::marker::PhantomData,
    _space: core::marker::PhantomData,
};

const fn _make_mat4_meters() -> Mat4<Meters, ()> {
    Mat4 {
        x_col: Vec4::<Meters, ()>::new(1.0, 2.0, 3.0, 4.0),
        y_col: Vec4::<Meters, ()>::new(5.0, 6.0, 7.0, 8.0),
        z_col: Vec4::<Meters, ()>::new(9.0, 10.0, 11.0, 12.0),
        w_col: Vec4::<Meters, ()>::new(13.0, 14.0, 15.0, 16.0),
        _unit: core::marker::PhantomData,
        _space: core::marker::PhantomData,
    }
}
const fn _make_mat4_pixels() -> Mat4<Pixels, ()> {
    Mat4 {
        x_col: Vec4::<Pixels, ()>::new(10.0, 20.0, 30.0, 40.0),
        y_col: Vec4::<Pixels, ()>::new(50.0, 60.0, 70.0, 80.0),
        z_col: Vec4::<Pixels, ()>::new(90.0, 100.0, 110.0, 120.0),
        w_col: Vec4::<Pixels, ()>::new(130.0, 140.0, 150.0, 160.0),
        _unit: core::marker::PhantomData,
        _space: core::marker::PhantomData,
    }
}
const _CONST_MAT4_METERS2: Mat4<Meters, ()> = _make_mat4_meters();
const _CONST_MAT4_PIXELS2: Mat4<Pixels, ()> = _make_mat4_pixels();

// --- Compile-time (const) tests/examples for Mat4<Unit, Space> phantom coordinate spaces ---
const _CONST_MAT4_WORLD: Mat4<(), World> = Mat4 {
    x_col: Vec4::<(), World>::new(1.0, 2.0, 3.0, 4.0),
    y_col: Vec4::<(), World>::new(5.0, 6.0, 7.0, 8.0),
    z_col: Vec4::<(), World>::new(9.0, 10.0, 11.0, 12.0),
    w_col: Vec4::<(), World>::new(13.0, 14.0, 15.0, 16.0),
    _unit: core::marker::PhantomData,
    _space: core::marker::PhantomData,
};
const _CONST_MAT4_LOCAL: Mat4<(), Local> = Mat4 {
    x_col: Vec4::<(), Local>::new(11.0, 12.0, 13.0, 14.0),
    y_col: Vec4::<(), Local>::new(15.0, 16.0, 17.0, 18.0),
    z_col: Vec4::<(), Local>::new(19.0, 20.0, 21.0, 22.0),
    w_col: Vec4::<(), Local>::new(23.0, 24.0, 25.0, 26.0),
    _unit: core::marker::PhantomData,
    _space: core::marker::PhantomData,
};
const _CONST_MAT4_SCREEN: Mat4<(), Screen> = Mat4 {
    x_col: Vec4::<(), Screen>::new(21.0, 22.0, 23.0, 24.0),
    y_col: Vec4::<(), Screen>::new(25.0, 26.0, 27.0, 28.0),
    z_col: Vec4::<(), Screen>::new(29.0, 30.0, 31.0, 32.0),
    w_col: Vec4::<(), Screen>::new(33.0, 34.0, 35.0, 36.0),
    _unit: core::marker::PhantomData,
    _space: core::marker::PhantomData,
};

const _: () = {
    // Compile-time assertions for const-everything
    assert!(_CONST_MAT4_ID.x_col.x == 1.0 && _CONST_MAT4_ID.y_col.y == 1.0 && _CONST_MAT4_ID.z_col.z == 1.0 && _CONST_MAT4_ID.w_col.w 
    == 1.0);
    assert!(_CONST_MAT4_ZERO.x_col.x == 0.0 && _CONST_MAT4_ZERO.y_col.y == 0.0 && _CONST_MAT4_ZERO.z_col.z == 0.0 && _CONST_MAT4_ZERO.
    w_col.w == 0.0);
    assert!(_CONST_MAT4_NEW.x_col.x == 1.0 && _CONST_MAT4_NEW.x_col.y == 2.0 && _CONST_MAT4_NEW.x_col.z == 3.0 && _CONST_MAT4_NEW.
    x_col.w == 4.0);
    assert!(_CONST_MAT4_NEW.y_col.x == 5.0 && _CONST_MAT4_NEW.y_col.y == 6.0 && _CONST_MAT4_NEW.y_col.z == 7.0 && _CONST_MAT4_NEW.
    y_col.w == 8.0);
    assert!(_CONST_MAT4_NEW.z_col.x == 9.0 && _CONST_MAT4_NEW.z_col.y == 10.0 && _CONST_MAT4_NEW.z_col.z == 11.0 && _CONST_MAT4_NEW.
    z_col.w == 12.0);
    assert!(_CONST_MAT4_NEW.w_col.x == 13.0 && _CONST_MAT4_NEW.w_col.y == 14.0 && _CONST_MAT4_NEW.w_col.z == 15.0 && _CONST_MAT4_NEW.
    w_col.w == 16.0);
    assert!(_CONST_MAT4_COLS.x_col.x == 1.0 && _CONST_MAT4_COLS.x_col.y == 2.0 && _CONST_MAT4_COLS.x_col.z == 3.0 && _CONST_MAT4_COLS.
    x_col.w == 4.0);
    assert!(_CONST_MAT4_COLS.y_col.x == 5.0 && _CONST_MAT4_COLS.y_col.y == 6.0 && _CONST_MAT4_COLS.y_col.z == 7.0 && _CONST_MAT4_COLS.
    y_col.w == 8.0);
    assert!(_CONST_MAT4_COLS.z_col.x == 9.0 && _CONST_MAT4_COLS.z_col.y == 10.0 && _CONST_MAT4_COLS.z_col.z == 11.0 && 
    _CONST_MAT4_COLS.z_col.w == 12.0);
    assert!(_CONST_MAT4_COLS.w_col.x == 13.0 && _CONST_MAT4_COLS.w_col.y == 14.0 && _CONST_MAT4_COLS.w_col.z == 15.0 && 
    _CONST_MAT4_COLS.w_col.w == 16.0);
    assert!(_CONST_MAT4_ROWS.x_col.x == 1.0 && _CONST_MAT4_ROWS.x_col.y == 2.0 && _CONST_MAT4_ROWS.x_col.z == 3.0 && _CONST_MAT4_ROWS.
    x_col.w == 4.0);
    assert!(_CONST_MAT4_ROWS.y_col.x == 5.0 && _CONST_MAT4_ROWS.y_col.y == 6.0 && _CONST_MAT4_ROWS.y_col.z == 7.0 && _CONST_MAT4_ROWS.
    y_col.w == 8.0);
    assert!(_CONST_MAT4_ROWS.z_col.x == 9.0 && _CONST_MAT4_ROWS.z_col.y == 10.0 && _CONST_MAT4_ROWS.z_col.z == 11.0 && 
    _CONST_MAT4_ROWS.z_col.w == 12.0);
    assert!(_CONST_MAT4_ROWS.w_col.x == 13.0 && _CONST_MAT4_ROWS.w_col.y == 14.0 && _CONST_MAT4_ROWS.w_col.z == 15.0 && 
    _CONST_MAT4_ROWS.w_col.w == 16.0);
    assert!(_CONST_MAT4_TRANSPOSE.x_col.x == 1.0 && _CONST_MAT4_TRANSPOSE.x_col.y == 5.0 && _CONST_MAT4_TRANSPOSE.x_col.z == 9.0 && 
    _CONST_MAT4_TRANSPOSE.x_col.w == 13.0);
    assert!(_CONST_MAT4_TRANSPOSE.y_col.x == 2.0 && _CONST_MAT4_TRANSPOSE.y_col.y == 6.0 && _CONST_MAT4_TRANSPOSE.y_col.z == 10.0 && 
    _CONST_MAT4_TRANSPOSE.y_col.w == 14.0);
    assert!(_CONST_MAT4_TRANSPOSE.z_col.x == 3.0 && _CONST_MAT4_TRANSPOSE.z_col.y == 7.0 && _CONST_MAT4_TRANSPOSE.z_col.z == 11.0 && 
    _CONST_MAT4_TRANSPOSE.z_col.w == 15.0);
    assert!(_CONST_MAT4_TRANSPOSE.w_col.x == 4.0 && _CONST_MAT4_TRANSPOSE.w_col.y == 8.0 && _CONST_MAT4_TRANSPOSE.w_col.z == 12.0 && 
    _CONST_MAT4_TRANSPOSE.w_col.w == 16.0);
    assert!(_CONST_MAT4_SHEAR.x_col.x == 1.0 && _CONST_MAT4_SHEAR.x_col.y == 3.0 && _CONST_MAT4_SHEAR.x_col.z == 5.0 && 
    _CONST_MAT4_SHEAR.x_col.w == 0.0);
    assert!(_CONST_MAT4_SHEAR.y_col.x == 1.0 && _CONST_MAT4_SHEAR.y_col.y == 1.0 && _CONST_MAT4_SHEAR.y_col.z == 6.0 && 
    _CONST_MAT4_SHEAR.y_col.w == 0.0);
    assert!(_CONST_MAT4_SHEAR.z_col.x == 2.0 && _CONST_MAT4_SHEAR.z_col.y == 4.0 && _CONST_MAT4_SHEAR.z_col.z == 1.0 && 
    _CONST_MAT4_SHEAR.z_col.w == 0.0);
    assert!(_CONST_MAT4_SHEAR.w_col.x == 0.0 && _CONST_MAT4_SHEAR.w_col.y == 0.0 && _CONST_MAT4_SHEAR.w_col.z == 0.0 && 
    _CONST_MAT4_SHEAR.w_col.w == 1.0);
    // Option values for row/col
    match _CONST_MAT4_ROW0 { Some(v) => assert!(v.x == 1.0 && v.y == 5.0 && v.z == 9.0 && v.w == 13.0), None => panic!("row0") }
    match _CONST_MAT4_ROW1 { Some(v) => assert!(v.x == 2.0 && v.y == 6.0 && v.z == 10.0 && v.w == 14.0), None => panic!("row1") }
    match _CONST_MAT4_ROW2 { Some(v) => assert!(v.x == 3.0 && v.y == 7.0 && v.z == 11.0 && v.w == 15.0), None => panic!("row2") }
    match _CONST_MAT4_ROW3 { Some(v) => assert!(v.x == 4.0 && v.y == 8.0 && v.z == 12.0 && v.w == 16.0), None => panic!("row3") }
    match _CONST_MAT4_COL0 { Some(v) => assert!(v.x == 1.0 && v.y == 2.0 && v.z == 3.0 && v.w == 4.0), None => panic!("col0") }
    match _CONST_MAT4_COL1 { Some(v) => assert!(v.x == 5.0 && v.y == 6.0 && v.z == 7.0 && v.w == 8.0), None => panic!("col1") }
    match _CONST_MAT4_COL2 { Some(v) => assert!(v.x == 9.0 && v.y == 10.0 && v.z == 11.0 && v.w == 12.0), None => panic!("col2") }
    match _CONST_MAT4_COL3 { Some(v) => assert!(v.x == 13.0 && v.y == 14.0 && v.z == 15.0 && v.w == 16.0), None => panic!("col3") }
    // Compile-time assertions for type-level units
    assert!(_CONST_MAT4_METERS.x_col.x == 1.0 && _CONST_MAT4_METERS.w_col.w == 16.0);
    assert!(_CONST_MAT4_PIXELS.x_col.x == 10.0 && _CONST_MAT4_PIXELS.w_col.w == 160.0);
    assert!(_CONST_MAT4_METERS2.x_col.x == 1.0 && _CONST_MAT4_METERS2.w_col.w == 16.0);
    assert!(_CONST_MAT4_PIXELS2.x_col.x == 10.0 && _CONST_MAT4_PIXELS2.w_col.w == 160.0);
    // Compile-time assertions for phantom coordinate spaces
    assert!(_CONST_MAT4_WORLD.x_col.x == 1.0 && _CONST_MAT4_WORLD.w_col.w == 16.0);
    assert!(_CONST_MAT4_LOCAL.x_col.x == 11.0 && _CONST_MAT4_LOCAL.w_col.w == 26.0);
    assert!(_CONST_MAT4_SCREEN.x_col.x == 21.0 && _CONST_MAT4_SCREEN.w_col.w == 36.0);
};