gemath 0.1.0

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

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


#[cfg(test)]
mod tests {
    use super::*;
    const EPSILON: f32 = 1e-6;
    const PI: f32 = std::f32::consts::PI;

    #[test]
    fn test_quat_new() {
        let q: Quat<(),()> = Quat::new(1.0, 2.0, 3.0, 4.0);
        assert!((q.x - 1.0).abs() < EPSILON);
        assert!((q.y - 2.0).abs() < EPSILON);
        assert!((q.z - 3.0).abs() < EPSILON);
        assert!((q.w - 4.0).abs() < EPSILON);
    }

    #[test]
    fn test_quat_identity() {
        let q: Quat<(),()> = Quat::IDENTITY;
        assert!((q.x - 0.0).abs() < EPSILON);
        assert!((q.y - 0.0).abs() < EPSILON);
        assert!((q.z - 0.0).abs() < EPSILON);
        assert!((q.w - 1.0).abs() < EPSILON);
    }

    #[test]
    fn test_quat_zero() {
        let q: Quat<(),()> = Quat::ZERO;
        assert!((q.x - 0.0).abs() < EPSILON);
        assert!((q.y - 0.0).abs() < EPSILON);
        assert!((q.z - 0.0).abs() < EPSILON);
        assert!((q.w - 0.0).abs() < EPSILON);
    }

    #[test]
    fn test_quat_default() {
        let q: Quat = Default::default(); // Should be 0,0,0,0 if Default is derived.
        // If w=1.0 is desired for Default, Default trait needs custom impl.
        assert_eq!(q, Quat::ZERO);
    }

    #[test]
    fn test_quat_from_vec4() {
        let v: Vec4<(),()> = Vec4::new(1.0, 2.0, 3.0, 4.0);
        let q = Quat::from_vec4(v);
        assert!((q.x - 1.0).abs() < EPSILON);
        assert!((q.y - 2.0).abs() < EPSILON);
        assert!((q.z - 3.0).abs() < EPSILON);
        assert!((q.w - 4.0).abs() < EPSILON);
    }

    #[test]
    fn test_quat_to_vec4() {
        let q: Quat<(),()> = Quat::new(1.0, 2.0, 3.0, 4.0);
        let v = q.to_vec4();
        assert!((v.x - 1.0).abs() < EPSILON);
        assert!((v.y - 2.0).abs() < EPSILON);
        assert!((v.z - 3.0).abs() < EPSILON);
        assert!((v.w - 4.0).abs() < EPSILON);
    }

    #[test]
    fn test_quat_xyz() {
        let q: Quat<(),()> = Quat::new(1.0, 2.0, 3.0, 4.0);
        let xyz = q.xyz();
        assert_eq!(xyz, Vec3::new(1.0, 2.0, 3.0));
    }

    #[test]
    fn test_quat_eq() {
        let q1: Quat<(),()> = Quat::new(1.0, 2.0, 3.0, 4.0);
        let q2: Quat<(),()> = Quat::new(1.0, 2.0, 3.0, 4.0);
        let q3: Quat<(),()> = Quat::new(4.0, 3.0, 2.0, 1.0);
        assert_eq!(q1, q2);
        assert_ne!(q1, q3);
    }

    #[test]
    fn test_quat_neg() {
        let q: Quat<(),()> = Quat::new(1.0, -2.0, 3.0, -4.0);
        let neg_q = -q;
        assert!((neg_q.x - (-1.0)).abs() < EPSILON);
        assert!((neg_q.y - 2.0).abs() < EPSILON);
        assert!((neg_q.z - (-3.0)).abs() < EPSILON);
        assert!((neg_q.w - 4.0).abs() < EPSILON);
    }

    #[test]
    fn test_quat_add() {
        let q1: Quat<(),()> = Quat::new(1.0, 2.0, 3.0, 4.0);
        let q2: Quat<(),()> = Quat::new(5.0, 6.0, 7.0, 8.0);
        let res = q1 + q2;
        assert_eq!(res, Quat::new(6.0, 8.0, 10.0, 12.0));
    }

    #[test]
    fn test_quat_add_assign() {
        let mut q1: Quat<(),()> = Quat::new(1.0, 2.0, 3.0, 4.0);
        let q2: Quat<(),()> = Quat::new(5.0, 6.0, 7.0, 8.0);
        q1 += q2;
        assert_eq!(q1, Quat::new(6.0, 8.0, 10.0, 12.0));
    }

    #[test]
    fn test_quat_sub() {
        let q1: Quat<(),()> = Quat::new(5.0, 8.0, 10.0, 12.0);
        let q2: Quat<(),()> = Quat::new(1.0, 2.0, 3.0, 4.0);
        let res = q1 - q2;
        assert_eq!(res, Quat::new(4.0, 6.0, 7.0, 8.0));
    }

    #[test]
    fn test_quat_sub_assign() {
        let mut q1: Quat<(),()> = Quat::new(5.0, 8.0, 10.0, 12.0);
        let q2: Quat<(),()> = Quat::new(1.0, 2.0, 3.0, 4.0);
        q1 -= q2;
        assert_eq!(q1, Quat::new(4.0, 6.0, 7.0, 8.0));
    }

    #[test]
    fn test_quat_mul_quat() {
        // Example from https://www.mathworks.com/help/nav/ref/quaternion.multiply.html
        // q1 = [1 0 1 0], q2 = [1 0.5 0.3 0.1]
        // In (w,x,y,z) order. Our struct is (x,y,z,w)
        // So q1_xyzw = (0,1,0,1), q2_xyzw = (0.5,0.3,0.1,1)
        // w_res = w1w2 - dot(v1,v2) = 1*1 - (0*0.5 + 1*0.3 + 0*0.1) = 1 - 0.3 = 0.7
        // v_res = w1v2 + w2v1 + cross(v1,v2)
        // w1v2 = 1*(0.5,0.3,0.1) = (0.5,0.3,0.1)
        // w2v1 = 1*(0,1,0) = (0,1,0)
        // cross(v1,v2) = cross((0,1,0), (0.5,0.3,0.1))
        //   cx = y1z2 - z1y2 = 1*0.1 - 0*0.3 = 0.1
        //   cy = z1x2 - x1z2 = 0*0.5 - 0*0.1 = 0
        //   cz = x1y2 - y1x2 = 0*0.3 - 1*0.5 = -0.5
        //   cross = (0.1, 0, -0.5)
        // v_res = (0.5,0.3,0.1) + (0,1,0) + (0.1,0,-0.5) = (0.6, 1.3, -0.4)
        // Result: (w,x,y,z) = (0.7, 0.6, 1.3, -0.4)
        // Result in (x,y,z,w): (0.6, 1.3, -0.4, 0.7)
        let q1: Quat<(),()> = Quat::new(0.0, 1.0, 0.0, 1.0); // w, x, y, z -> x,y,z,w
        let q2: Quat<(),()> = Quat::new(0.5, 0.3, 0.1, 1.0);

        let res = q1 * q2;
        assert!((res.x - 0.6).abs() < EPSILON);
        assert!((res.y - 1.3).abs() < EPSILON);
        assert!((res.z + 0.4).abs() < EPSILON);
        assert!((res.w - 0.7).abs() < EPSILON);

        let q_id: Quat<(),()> = Quat::IDENTITY;
        let q_val: Quat<(),()> = Quat::new(1.0, 2.0, 3.0, 4.0);
        assert_eq!(q_id * q_val, q_val);
        assert_eq!(q_val * q_id, q_val);
    }

    #[test]
    fn test_quat_mul_assign_quat() {
        let mut q1: Quat<(),()> = Quat::new(0.0, 1.0, 0.0, 1.0);
        let q2 = Quat::new(0.5, 0.3, 0.1, 1.0);
        q1 *= q2;
        assert!((q1.x - 0.6).abs() < EPSILON);
        assert!((q1.y - 1.3).abs() < EPSILON);
        assert!((q1.z + 0.4).abs() < EPSILON);
        assert!((q1.w - 0.7).abs() < EPSILON);
    }

    #[test]
    fn test_quat_mul_f32() {
        let q: Quat<(),()> = Quat::new(1.0, 2.0, 3.0, 4.0);
        let s = 2.0;
        let res_qs = q * s;
        let res_sq = s * q;
        assert_eq!(res_qs, Quat::new(2.0, 4.0, 6.0, 8.0));
        assert_eq!(res_sq, Quat::new(2.0, 4.0, 6.0, 8.0));
    }

    #[test]
    fn test_quat_mul_assign_f32() {
        let mut q: Quat<(),()> = Quat::new(1.0, 2.0, 3.0, 4.0);
        q *= 2.0;
        assert_eq!(q, Quat::new(2.0, 4.0, 6.0, 8.0));
    }

    #[test]
    fn test_quat_div_f32() {
        let q: Quat<(),()> = Quat::new(2.0, 4.0, 6.0, 8.0);
        let s = 2.0;
        let res = q / s;
        assert_eq!(res, Quat::new(1.0, 2.0, 3.0, 4.0));
    }

    #[test]
    fn test_quat_checked_div_scalar() {
        let q: Quat<(),()> = Quat::new(2.0, 4.0, 6.0, 8.0);
        assert_eq!(q.checked_div_scalar(2.0), Some(Quat::new(1.0, 2.0, 3.0, 4.0)));
        assert_eq!(q.checked_div_scalar(0.0), None);
        assert_eq!(q.checked_div_scalar(-0.0), None);
        assert_eq!(q.checked_div_scalar(f32::NAN), None);
        assert_eq!(q.checked_div_scalar(f32::INFINITY), None);
    }

    #[test]
    fn test_quat_div_assign_f32() {
        let mut q: Quat<(),()> = Quat::new(2.0, 4.0, 6.0, 8.0);
        q /= 2.0;
        assert_eq!(q, Quat::new(1.0, 2.0, 3.0, 4.0));
    }

    #[test]
    fn test_quat_mul_vec3_rotation() {
        // Rotate (1,0,0) by 90 degrees (PI/2) around Z-axis (0,0,1)
        // Expected: (0,1,0)
        let axis: Vec3<(),()> = Vec3::new(0.0, 0.0, 1.0);
        let angle = PI / 2.0;
        // let q = Quat::from_axis_angle(axis, angle); // Use this when from_axis_angle is confirmed in main code
        let q = Quat::from_axis_angle_radians(
            axis.normalize(),
            gemath::angle::Radians(angle),
        );

        let v = Vec3::new(1.0, 0.0, 0.0);
        let rotated_v = q * v;

        assert!(
            (rotated_v.x - 0.0).abs() < EPSILON,
            "rotated_v.x: {}",
            rotated_v.x
        );
        assert!(
            (rotated_v.y - 1.0).abs() < EPSILON,
            "rotated_v.y: {}",
            rotated_v.y
        );
        assert!(
            (rotated_v.z - 0.0).abs() < EPSILON,
            "rotated_v.z: {}",
            rotated_v.z
        );

        // Rotate (1,0,0) by 180 degrees (PI) around Y-axis (0,1,0)
        // Expected: (-1,0,0)
        let axis2: Vec3<(),()> = Vec3::new(0.0, 1.0, 0.0);
        let angle2 = PI;
        // let q2 = Quat::from_axis_angle(axis2, angle2);
        let q2 = Quat::from_axis_angle_radians(
            axis2.normalize(),
            gemath::angle::Radians(angle2),
        );
        let v2 = Vec3::new(1.0, 0.0, 0.0);
        let rotated_v2 = q2 * v2;

        assert!(
            (rotated_v2.x - (-1.0)).abs() < EPSILON,
            "rotated_v2.x: {}",
            rotated_v2.x
        );
        assert!(
            (rotated_v2.y - 0.0).abs() < EPSILON,
            "rotated_v2.y: {}",
            rotated_v2.y
        );
        assert!(
            (rotated_v2.z - 0.0).abs() < EPSILON,
            "rotated_v2.z: {}",
            rotated_v2.z
        );

        // Test with identity quaternion
        let q_identity: Quat<(),()> = Quat::IDENTITY;
        let v_orig = Vec3::new(1.2, 3.4, 5.6);
        let rotated_v_id = q_identity * v_orig;
        assert!((rotated_v_id.x - v_orig.x).abs() < EPSILON);
        assert!((rotated_v_id.y - v_orig.y).abs() < EPSILON);
        assert!((rotated_v_id.z - v_orig.z).abs() < EPSILON);
    }

    #[test]
    fn test_quat_from_axis_angle() {
        // Rotate 90 degrees around Y axis
        let axis: Vec3<(),()> = Vec3::new(0.0, 1.0, 0.0);
        let angle = PI / 2.0;
        let q = Quat::from_axis_angle_radians(axis, gemath::angle::Radians(angle));
        assert!((q.x - 0.0).abs() < EPSILON);
        assert!((q.y - (PI / 4.0).sin()).abs() < EPSILON);
        assert!((q.z - 0.0).abs() < EPSILON);
        assert!((q.w - (PI / 4.0).cos()).abs() < EPSILON);

        // Rotate 180 degrees around X axis
        let axis_x: Vec3<(),()> = Vec3::new(1.0, 0.0, 0.0);
        let angle_pi = PI;
        let q_x = Quat::from_axis_angle_radians(axis_x, gemath::angle::Radians(angle_pi));
        assert!((q_x.x - 1.0).abs() < EPSILON); // sin(pi/2) = 1
        assert!((q_x.y - 0.0).abs() < EPSILON);
        assert!((q_x.z - 0.0).abs() < EPSILON);
        assert!((q_x.w - 0.0).abs() < EPSILON); // cos(pi/2) = 0

        // Zero angle
        let q_zero_angle: Quat<(),()> = Quat::from_axis_angle_radians(
            Vec3::new(1.0, 2.0, 3.0),
            gemath::angle::Radians(0.0),
        );
        assert_eq!(q_zero_angle, Quat::IDENTITY);

        // Zero axis (should result in identity or a zero quaternion based on normalization)
        // gb_vec3_norm of zero vector results in zero vector. sin(angle/2)*0 = 0.
        // cos(angle/2) remains.
        let q_zero_axis: Quat<(),()> =
            Quat::from_axis_angle_radians(Vec3::ZERO, gemath::angle::Radians(PI / 3.0));
        assert!((q_zero_axis.x).abs() < EPSILON);
        assert!((q_zero_axis.y).abs() < EPSILON);
        assert!((q_zero_axis.z).abs() < EPSILON);
        assert!((q_zero_axis.w - (PI / 6.0).cos()).abs() < EPSILON); // Not identity, but C++ behavior for zero axis norm
        // The C++ gb_quat_axis_angle normalizes axis.
        // If axis is zero, normalized is zero.
        // So xyz components become 0. w becomes cos(angle/2).
    }

    #[test]
    fn test_quat_from_euler_angles() {
        // Pitch 90 deg
        let q_pitch: Quat<(),()> = Quat::from_euler_angles_radians(
            gemath::angle::Radians(PI / 2.0),
            gemath::angle::Radians(0.0),
            gemath::angle::Radians(0.0),
        );
        let expected_pitch: Quat<(),()> = Quat::from_axis_angle_radians(
            Vec3::new(1.0, 0.0, 0.0),
            gemath::angle::Radians(PI / 2.0),
        );
        assert!(
            (q_pitch.x - expected_pitch.x).abs() < EPSILON,
            "Pitch X: {} vs {}",
            q_pitch.x,
            expected_pitch.x
        );
        assert!(
            (q_pitch.y - expected_pitch.y).abs() < EPSILON,
            "Pitch Y: {} vs {}",
            q_pitch.y,
            expected_pitch.y
        );
        assert!(
            (q_pitch.z - expected_pitch.z).abs() < EPSILON,
            "Pitch Z: {} vs {}",
            q_pitch.z,
            expected_pitch.z
        );
        assert!(
            (q_pitch.w - expected_pitch.w).abs() < EPSILON,
            "Pitch W: {} vs {}",
            q_pitch.w,
            expected_pitch.w
        );

        // Yaw 90 deg
        let q_yaw: Quat<(),()> = Quat::from_euler_angles_radians(
            gemath::angle::Radians(0.0),
            gemath::angle::Radians(PI / 2.0),
            gemath::angle::Radians(0.0),
        );
        let expected_yaw: Quat<(),()> = Quat::from_axis_angle_radians(
            Vec3::new(0.0, 1.0, 0.0),
            gemath::angle::Radians(PI / 2.0),
        );
        assert!(
            (q_yaw.x - expected_yaw.x).abs() < EPSILON,
            "Yaw X: {} vs {}",
            q_yaw.x,
            expected_yaw.x
        );
        assert!(
            (q_yaw.y - expected_yaw.y).abs() < EPSILON,
            "Yaw Y: {} vs {}",
            q_yaw.y,
            expected_yaw.y
        );
        assert!(
            (q_yaw.z - expected_yaw.z).abs() < EPSILON,
            "Yaw Z: {} vs {}",
            q_yaw.z,
            expected_yaw.z
        );
        assert!(
            (q_yaw.w - expected_yaw.w).abs() < EPSILON,
            "Yaw W: {} vs {}",
            q_yaw.w,
            expected_yaw.w
        );

        // Roll 90 deg
        let q_roll: Quat<(),()> = Quat::from_euler_angles_radians(
            gemath::angle::Radians(0.0),
            gemath::angle::Radians(0.0),
            gemath::angle::Radians(PI / 2.0),
        );
        let expected_roll: Quat<(),()> = Quat::from_axis_angle_radians(
            Vec3::new(0.0, 0.0, 1.0),
            gemath::angle::Radians(PI / 2.0),
        );
        assert!(
            (q_roll.x - expected_roll.x).abs() < EPSILON,
            "Roll X: {} vs {}",
            q_roll.x,
            expected_roll.x
        );
        assert!(
            (q_roll.y - expected_roll.y).abs() < EPSILON,
            "Roll Y: {} vs {}",
            q_roll.y,
            expected_roll.y
        );
        assert!(
            (q_roll.z - expected_roll.z).abs() < EPSILON,
            "Roll Z: {} vs {}",
            q_roll.z,
            expected_roll.z
        );
        assert!(
            (q_roll.w - expected_roll.w).abs() < EPSILON,
            "Roll W: {} vs {}",
            q_roll.w,
            expected_roll.w
        );

        let q_pyr: Quat<(),()> = Quat::from_euler_angles_radians(
            gemath::angle::Radians(0.1),
            gemath::angle::Radians(0.2),
            gemath::angle::Radians(0.3),
        );

        assert!(
            (q_pyr.x - 0.064071).abs() < EPSILON,
            "Euler X: {} vs {}",
            q_pyr.x,
            0.064071
        );
        assert!(
            (q_pyr.y - 0.091157).abs() < EPSILON,
            "Euler Y: {} vs {}",
            q_pyr.y,
            0.091157
        );
        assert!(
            (q_pyr.z - 0.143572).abs() < EPSILON,
            "Euler Z: {} vs {}",
            q_pyr.z,
            0.143572
        );
        assert!(
            (q_pyr.w - 0.983347).abs() < EPSILON,
            "Euler W: {} vs {}",
            q_pyr.w,
            0.983347
        );
    }

    #[test]
    fn test_quat_dot() {
        let q1: Quat<(),()> = Quat::new(1.0, 0.0, 0.0, 0.0);
        let q2: Quat<(),()> = Quat::new(0.0, 1.0, 0.0, 0.0);
        let q3: Quat<(),()> = Quat::new(1.0, 2.0, 3.0, 4.0);
        let q4: Quat<(),()> = Quat::new(4.0, 3.0, 2.0, 1.0);
        assert!((q1.dot(q1) - 1.0).abs() < EPSILON);
        assert!((q1.dot(q2) - 0.0).abs() < EPSILON);
        assert!((q3.dot(q4) - (4.0 + 6.0 + 6.0 + 4.0)).abs() < EPSILON); // 20
    }

    #[test]
    fn test_quat_length() {
        let q1: Quat<(),()> = Quat::new(1.0, 0.0, 0.0, 0.0);
        let q2: Quat<(),()> = Quat::new(0.0, 0.0, 0.0, 1.0);
        let q3: Quat<(),()> = Quat::new(1.0, 2.0, 3.0, 4.0); // len_sq = 1+4+9+16 = 30
        assert!((q1.length_squared() - 1.0).abs() < EPSILON);
        assert!((q1.length() - 1.0).abs() < EPSILON);
        assert!((q2.length_squared() - 1.0).abs() < EPSILON);
        assert!((q2.length() - 1.0).abs() < EPSILON);
        assert!((q3.length_squared() - 30.0).abs() < EPSILON);
        assert!((q3.length() - 30.0_f32.sqrt()).abs() < EPSILON);
        assert!((Quat::<(),()>::ZERO.length_squared()).abs() < EPSILON);
        assert!((Quat::<(),()>::ZERO.length()).abs() < EPSILON);
    }

    #[test]
    fn test_quat_normalize() {
        let q1: Quat<(),()> = Quat::new(1.0, 2.0, 3.0, 4.0);
        let n1 = q1.normalize();
        assert!(
            (n1.length() - 1.0).abs() < EPSILON,
            "Normalized length: {}",
            n1.length()
        );
        let len = q1.length();
        assert!((n1.x - q1.x / len).abs() < EPSILON);
        assert!((n1.y - q1.y / len).abs() < EPSILON);
        assert!((n1.z - q1.z / len).abs() < EPSILON);
        assert!((n1.w - q1.w / len).abs() < EPSILON);

        let q_zero: Quat<(),()> = Quat::ZERO;
        let n_zero = q_zero.normalize(); // Should be zero
        assert_eq!(n_zero, Quat::ZERO);

        let q_ident: Quat<(),()> = Quat::IDENTITY;
        let n_ident = q_ident.normalize();
        assert_eq!(n_ident, Quat::IDENTITY);

        let tn1 = q1.try_normalize().unwrap();
        assert!((tn1.length() - 1.0).abs() < EPSILON);
        assert!(q_zero.try_normalize().is_none());
    }

    #[test]
    fn test_quat_conjugate() {
        let q: Quat<(),()> = Quat::new(1.0, 2.0, 3.0, 4.0);
        let conj = q.conjugate();
        assert_eq!(conj, Quat::new(-1.0, -2.0, -3.0, 4.0));
        assert_eq!(Quat::<(),()>::IDENTITY.conjugate(), Quat::<(),()>::IDENTITY);
    }

    #[test]
    fn test_quat_inverse() {
        let q: Quat<(),()> = Quat::new(1.0, 2.0, 3.0, 4.0); // len_sq = 30
        let inv_q = q.inverse().unwrap();
        let expected_inv = q.conjugate() / q.length_squared();
        assert_eq!(inv_q, expected_inv);

        // q * q_inv should be identity
        let prod = q * inv_q;
        assert!(
            (prod.x - Quat::<(),()>::IDENTITY.x).abs() < EPSILON,
            "prod.x: {}",
            prod.x
        );
        assert!(
            (prod.y - Quat::<(),()>::IDENTITY.y).abs() < EPSILON,
            "prod.y: {}",
            prod.y
        );
        assert!(
            (prod.z - Quat::<(),()>::IDENTITY.z).abs() < EPSILON,
            "prod.z: {}",
            prod.z
        );
        assert!(
            (prod.w - Quat::<(),()>::IDENTITY.w).abs() < EPSILON,
            "prod.w: {}",
            prod.w
        );

        assert!(Quat::<(),()>::ZERO.inverse().is_none());
        let id_inv = Quat::<(),()>::IDENTITY.inverse().unwrap();
        assert_eq!(id_inv, Quat::<(),()>::IDENTITY);
    }

    #[test]
    fn test_quat_try_inverse_alias() {
        let q: Quat<(),()> = Quat::new(0.1, 0.2, 0.3, 0.4);
        assert_eq!(q.try_inverse(), q.inverse());

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

    #[test]
    fn test_quat_to_axis_angle() {
        // 90 deg around Y
        let q_y_90: Quat<(),()> = Quat::from_axis_angle_radians(
            Vec3::new(0.0, 1.0, 0.0),
            gemath::angle::Radians(PI / 2.0),
        );
        let (axis, angle) = q_y_90.to_axis_angle();
        assert!((axis.x - 0.0).abs() < EPSILON, "Axis X: {}", axis.x);
        assert!((axis.y - 1.0).abs() < EPSILON, "Axis Y: {}", axis.y);
        assert!((axis.z - 0.0).abs() < EPSILON, "Axis Z: {}", axis.z);
        assert!((angle - PI / 2.0).abs() < EPSILON, "Angle: {}", angle);

        // Identity quaternion
        let (id_axis, id_angle) = Quat::<(),()>::IDENTITY.to_axis_angle();
        assert!(
            (id_angle - 0.0).abs() < EPSILON,
            "Identity angle: {}",
            id_angle
        );
        // Axis can be arbitrary for 0 angle, gb_math C returns (1,0,0) effectively.
        // Our implementation returns (1,0,0) for identity due to the s < EPSILON check.
        assert!(
            (id_axis.x - 1.0).abs() < EPSILON,
            "Identity axis X: {}",
            id_axis.x
        );

        // 180 deg around X
        let q_x_180: Quat<(),()> =
            Quat::from_axis_angle_radians(Vec3::new(1.0, 0.0, 0.0), gemath::angle::Radians(PI)); // x=1, y=0, z=0, w=0
        let (axis_x180, angle_x180) = q_x_180.to_axis_angle();
        assert!(
            (axis_x180.x - 1.0).abs() < EPSILON,
            "Axis X180 X: {}",
            axis_x180.x
        );
        assert!(
            (axis_x180.y - 0.0).abs() < EPSILON,
            "Axis X180 Y: {}",
            axis_x180.y
        );
        assert!(
            (axis_x180.z - 0.0).abs() < EPSILON,
            "Axis X180 Z: {}",
            axis_x180.z
        );
        assert!(
            (angle_x180 - PI).abs() < EPSILON,
            "Angle X180: {}",
            angle_x180
        );

        // A more complex quaternion approx 180 deg around (1,2,3)
        let q_complex: Quat<(),()> = Quat::new(0.267261, 0.534522, 0.801784, 0.0).normalize();
        let (axis_c, angle_c) = q_complex.to_axis_angle();
        assert!(
            (axis_c.x - 0.267261).abs() < EPSILON,
            "Axis C X: {}",
            axis_c.x
        );
        assert!(
            (axis_c.y - 0.534522).abs() < EPSILON,
            "Axis C Y: {}",
            axis_c.y
        );
        assert!(
            (axis_c.z - 0.801784).abs() < EPSILON,
            "Axis C Z: {}",
            axis_c.z
        );
        assert!((angle_c - 3.141593).abs() < EPSILON, "Angle C: {}", angle_c);

        // A more complex quaternion approx 120 deg around (1,2,3)
        // Expected axis (1,2,3) normalized = (0.267261, 0.534522, 0.801784)
        // Expected angle approx 120 deg = 2.094395 rad
        // If w=0 for normalized quat, angle = PI. Here w is not 0 after normalize.
        // q_complex.w is not 0 after normalize, it's q_complex.w / length

        /* Explanation:
           To construct a quaternion representing a 120° rotation (in radians: θ = 120° = 2.094395 radians) around the axis (1, 2, 3), using the `Quat::new(x, y, z, w)` API, we need to:
           1.	Normalize the axis (which we already did: (0.267261, 0.534522, 0.801784))
           2.	Compute quaternion components using the axis-angle formula:
           • x = axis.x * sin(θ/2)
           • y = axis.y * sin(θ/2)
           • z = axis.z * sin(θ/2)
           • w = cos(θ/2)

           Compute the values:
           1. Angle in radians and half-angle
               •	θ = 120° = 2.094395 radians
               •	θ/2 = 1.0471975 radians
           2. Sine and cosine of half-angle
               •	sin(θ/2) ≈ sin(1.0471975) ≈ 0.8660254
               •	cos(θ/2) ≈ cos(1.0471975) ≈ 0.5
           3. Quaternion components
               •	x = 0.267261 × 0.8660254 ≈ 0.231455
               •	y = 0.534522 × 0.8660254 ≈ 0.462910
               •	z = 0.801784 × 0.8660254 ≈ 0.694365
               •	w = 0.5
           4. Final construction
           Quat::new(0.231455, 0.462910, 0.694365, 0.5)
        */
        let q_complex: Quat<(),()> = Quat::new(0.231455, 0.462910, 0.694365, 0.5).normalize();
        let (axis_c, angle_c) = q_complex.to_axis_angle();
        assert!(
            (axis_c.x - 0.267261).abs() < EPSILON,
            "Axis C X: {}",
            axis_c.x
        );
        assert!(
            (axis_c.y - 0.534522).abs() < EPSILON,
            "Axis C Y: {}",
            axis_c.y
        );
        assert!(
            (axis_c.z - 0.801784).abs() < EPSILON,
            "Axis C Z: {}",
            axis_c.z
        );
        assert!((angle_c - 2.094395).abs() < EPSILON, "Angle C: {}", angle_c);

        /*
            New test:
            Purpose:
            •	Primary Goal: Ensure the conversion algorithm correctly implements the mathematical formulas
            •	Secondary Benefit: Catches errors in floating-point handling and normalization
            •	Smart Approach: Avoids external reference values that could contain their own calculation errors
        */
        let q_norm: Quat<(),()> = Quat::new(1.0, 2.0, 3.0, 0.5).normalize(); // Example
        let (axis_n, angle_n) = q_norm.to_axis_angle();
        let expected_angle_n = 2.0 * q_norm.w.acos();
        let s_n = (1.0 - q_norm.w * q_norm.w).sqrt();
        let expected_axis_n = q_norm.xyz() / s_n;

        assert!(
            (angle_n - expected_angle_n).abs() < EPSILON,
            "Angle N: {} vs {}",
            angle_n,
            expected_angle_n
        );
        assert!(
            (axis_n.x - expected_axis_n.x).abs() < EPSILON,
            "Axis N X: {} vs {}",
            axis_n.x,
            expected_axis_n.x
        );
        assert!(
            (axis_n.y - expected_axis_n.y).abs() < EPSILON,
            "Axis N Y: {} vs {}",
            axis_n.y,
            expected_axis_n.y
        );
        assert!(
            (axis_n.z - expected_axis_n.z).abs() < EPSILON,
            "Axis N Z: {} vs {}",
            axis_n.z,
            expected_axis_n.z
        );
    }

    #[test]
    fn test_quat_pitch_yaw_roll() {
        // Pure pitch 90 deg
        let q_pitch: Quat<(),()> = Quat::from_axis_angle_radians(
            Vec3::new(1.0, 0.0, 0.0),
            gemath::angle::Radians(PI / 2.0),
        );
        assert!(
            (q_pitch.pitch() - PI / 2.0).abs() < EPSILON,
            "Pitch: {}",
            q_pitch.pitch()
        );
        assert!(
            (q_pitch.yaw() - 0.0).abs() < EPSILON,
            "Yaw for pitch: {}",
            q_pitch.yaw()
        );
        assert!(
            (q_pitch.roll() - 0.0).abs() < EPSILON,
            "Roll for pitch: {}",
            q_pitch.roll()
        );

        // Pure yaw 90 deg
        let q_yaw: Quat<(),()> = Quat::from_axis_angle_radians(
            Vec3::new(0.0, 1.0, 0.0),
            gemath::angle::Radians(PI / 2.0),
        );
        assert!(
            (q_yaw.pitch() - 0.0).abs() < EPSILON,
            "Pitch for yaw: {}",
            q_yaw.pitch()
        );
        assert!(
            // TODO: why is this not 1e-6 at least?
            (q_yaw.yaw() - PI / 2.0).abs() < 1e-3,
            "Yaw: {}",
            q_yaw.yaw()
        );
        assert!(
            (q_yaw.roll() - 0.0).abs() < EPSILON,
            "Roll for yaw: {}",
            q_yaw.roll()
        );

        // Pure roll 90 deg
        let q_roll: Quat<(),()> = Quat::from_axis_angle_radians(
            Vec3::new(0.0, 0.0, 1.0),
            gemath::angle::Radians(PI / 2.0),
        );
        assert!(
            (q_roll.pitch() - 0.0).abs() < EPSILON,
            "Pitch for roll: {}",
            q_roll.pitch()
        );
        assert!(
            (q_roll.yaw() - 0.0).abs() < EPSILON,
            "Yaw for roll: {}",
            q_roll.yaw()
        );
        assert!(
            (q_roll.roll() - PI / 2.0).abs() < EPSILON,
            "Roll: {}",
            q_roll.roll()
        );

        // Combined: R(PI/6) Y(PI/4) P(PI/3)
        // Order of application in gb_quat_euler_angles: Y then P then R.
        let q_combi: Quat<(),()> = Quat::from_euler_angles_radians(
            gemath::angle::Radians(PI / 3.0),
            gemath::angle::Radians(PI / 4.0),
            gemath::angle::Radians(PI / 6.0),
        );
        // Note: recovering Euler angles can be tricky due to gimbal lock and multiple representations.
        // The C++ functions gb_quat_pitch/yaw/roll match specific formulas.
        // We test that our Rust functions match those formulas by re-deriving from the quaternion.
        let expected_pitch = (2.0 * (q_combi.y * q_combi.z + q_combi.w * q_combi.x)).atan2(
            q_combi.w * q_combi.w - q_combi.x * q_combi.x - q_combi.y * q_combi.y
                + q_combi.z * q_combi.z,
        );
        let expected_yaw = (-2.0 * (q_combi.x * q_combi.z - q_combi.w * q_combi.y)).asin();
        let expected_roll = (2.0 * (q_combi.x * q_combi.y + q_combi.w * q_combi.z)).atan2(
            q_combi.w * q_combi.w + q_combi.x * q_combi.x
                - q_combi.y * q_combi.y
                - q_combi.z * q_combi.z,
        );

        assert!(
            (q_combi.pitch() - expected_pitch).abs() < EPSILON,
            "Combi Pitch: {} vs {}",
            q_combi.pitch(),
            expected_pitch
        );
        assert!(
            (q_combi.yaw() - expected_yaw).abs() < EPSILON,
            "Combi Yaw: {} vs {}",
            q_combi.yaw(),
            expected_yaw
        );
        assert!(
            (q_combi.roll() - expected_roll).abs() < EPSILON,
            "Combi Roll: {} vs {}",
            q_combi.roll(),
            expected_roll
        );
    }

    // Helper for from_mat4 tests
    fn mat4_from_quat_manual(q: Quat) -> Mat4 {
        // Matches gb_mat4_from_quat
        let mut m = Mat4::IDENTITY;
        let nq = q.normalize(); // C++ gb_mat4_from_quat normalizes input quat

        let xx = nq.x * nq.x;
        let yy = nq.y * nq.y;
        let zz = nq.z * nq.z;
        let xy = nq.x * nq.y;
        let xz = nq.x * nq.z;
        let yz = nq.y * nq.z;
        let wx = nq.w * nq.x;
        let wy = nq.w * nq.y;
        let wz = nq.w * nq.z;

        m.x_col.x = 1.0 - 2.0 * (yy + zz);
        m.y_col.x = 2.0 * (xy - wz); // m[0][1] in C++ row-major, m.y_col.x in Rust col-major
        m.z_col.x = 2.0 * (xz + wy); // m[0][2] in C++ row-major

        m.x_col.y = 2.0 * (xy + wz); // m[1][0] in C++ row-major
        m.y_col.y = 1.0 - 2.0 * (xx + zz);
        m.z_col.y = 2.0 * (yz - wx); // m[1][2] in C++ row-major

        m.x_col.z = 2.0 * (xz - wy); // m[2][0] in C++ row-major
        m.y_col.z = 2.0 * (yz + wx); // m[2][1] in C++ row-major
        m.z_col.z = 1.0 - 2.0 * (xx + yy);
        m
    }

    fn assert_quat_eq_approx(q1: Quat, q2: Quat, epsilon: f32, msg: &str) {
        // Check if q1 is approximately q2 or -q2
        let direct_match = (q1.x - q2.x).abs() < epsilon
            && (q1.y - q2.y).abs() < epsilon
            && (q1.z - q2.z).abs() < epsilon
            && (q1.w - q2.w).abs() < epsilon;

        let negated_match = (q1.x + q2.x).abs() < epsilon
            && (q1.y + q2.y).abs() < epsilon
            && (q1.z + q2.z).abs() < epsilon
            && (q1.w + q2.w).abs() < epsilon;

        if direct_match || negated_match {
            // It's a match
        } else {
            // Log details if no match
            assert!(
                (q1.x - q2.x).abs() < epsilon || (q1.x + q2.x).abs() < epsilon,
                "{}: X mismatch: q1.x={}, q2.x={}",
                msg,
                q1.x,
                q2.x
            );
            assert!(
                (q1.y - q2.y).abs() < epsilon || (q1.y + q2.y).abs() < epsilon,
                "{}: Y mismatch: q1.y={}, q2.y={}",
                msg,
                q1.y,
                q2.y
            );
            assert!(
                (q1.z - q2.z).abs() < epsilon || (q1.z + q2.z).abs() < epsilon,
                "{}: Z mismatch: q1.z={}, q2.z={}",
                msg,
                q1.z,
                q2.z
            );
            assert!(
                (q1.w - q2.w).abs() < epsilon || (q1.w + q2.w).abs() < epsilon,
                "{}: W mismatch: q1.w={}, q2.w={}",
                msg,
                q1.w,
                q2.w
            );
            // If any of the above assertions fail, the test fails with a specific component message.
            // If we reach here, it means the combination failed but individual component checks would also fail.
            // For a clearer overall message if the logic itself is the point of test:
            panic!("{}: Quaternion mismatch. q1: {:?}, q2: {:?}", msg, q1, q2);
        }
    }

    #[test]
    fn test_quat_from_mat4() {
        // Identity matrix
        let m_ident = Mat4::IDENTITY;
        let q_ident_from_m = Quat::from_mat4(&m_ident);
        // Check if q_ident_from_m is (0,0,0,1) or (0,0,0,-1)
        assert_quat_eq_approx(
            q_ident_from_m,
            Quat::IDENTITY,
            EPSILON,
            "Identity Mat4 -> Quat",
        );

        // Rotation 90 deg around Y axis
        let q_y_90 = Quat::from_axis_angle_radians(
            Vec3::new(0.0, 1.0, 0.0),
            gemath::angle::Radians(PI / 2.0),
        );
        let m_y_90 = mat4_from_quat_manual(q_y_90); // Use our Mat4::from_quat if available and tested
        // For now, mat4_from_quat_manual matches gb_mat4_from_quat
        let q_y_90_from_m = Quat::from_mat4(&m_y_90);
        assert_quat_eq_approx(q_y_90_from_m, q_y_90, EPSILON, "Y-90 Mat4 -> Quat");

        // Rotation 180 deg around X axis
        let q_x_180 =
            Quat::from_axis_angle_radians(Vec3::new(1.0, 0.0, 0.0), gemath::angle::Radians(PI));
        let m_x_180 = mat4_from_quat_manual(q_x_180);
        let q_x_180_from_m = Quat::from_mat4(&m_x_180);
        assert_quat_eq_approx(q_x_180_from_m, q_x_180, EPSILON, "X-180 Mat4 -> Quat");

        // More complex rotation
        let q_complex = Quat::from_euler_angles_radians(
            gemath::angle::Radians(PI / 6.0),
            gemath::angle::Radians(PI / 4.0),
            gemath::angle::Radians(PI / 3.0),
        )
        .normalize();
        let m_complex = mat4_from_quat_manual(q_complex);
        let q_complex_from_m = Quat::from_mat4(&m_complex);
        assert_quat_eq_approx(q_complex_from_m, q_complex, EPSILON, "Complex Mat4 -> Quat");

        // Test case from gb_math.cpp comments (if any specific matrix values for quat conversion)
        // For now, use known quaternion, convert to matrix, convert back.
    }

    fn quat_approx_eq(a: Quat, b: Quat, eps: f32) -> bool {
        (a.x - b.x).abs() < eps
            && (a.y - b.y).abs() < eps
            && (a.z - b.z).abs() < eps
            && (a.w - b.w).abs() < eps
    }

    #[test]
    fn test_quat_interpolations() {
        // Basic quaternions
        let q1 = Quat::IDENTITY;
        let q2 = Quat::from_axis_angle_radians(
            Vec3::new(0.0, 1.0, 0.0),
            gemath::angle::Radians(PI),
        ); // 180 deg around Y
        let q3 = Quat::from_axis_angle_radians(
            Vec3::new(1.0, 0.0, 0.0),
            gemath::angle::Radians(PI / 2.0),
        ); // 90 deg around X
        let q4 = Quat::from_axis_angle_radians(
            Vec3::new(0.0, 0.0, 1.0),
            gemath::angle::Radians(PI / 2.0),
        ); // 90 deg around Z

        // LERP: t=0 returns a, t=1 returns b, t=0.5 is midpoint (not normalized)
        let lerp_0 = q1.lerp(q2, 0.0);
        let lerp_1 = q1.lerp(q2, 1.0);
        let lerp_half = q1.lerp(q2, 0.5);
        assert!(quat_approx_eq(lerp_0, q1, EPSILON));
        assert!(quat_approx_eq(lerp_1, q2, EPSILON));
        // LERP midpoint is not normalized, but should be between q1 and q2
        assert!((lerp_half.length() - 1.0).abs() > 1e-3); // Not normalized

        // NLERP: t=0 returns a, t=1 returns b, t=0.5 is normalized midpoint
        let nlerp_0 = q1.nlerp(q2, 0.0);
        let nlerp_1 = q1.nlerp(q2, 1.0);
        let nlerp_half = q1.nlerp(q2, 0.5);
        assert!(quat_approx_eq(nlerp_0, q1.normalize(), EPSILON));
        assert!(quat_approx_eq(nlerp_1, q2.normalize(), EPSILON));
        assert!((nlerp_half.length() - 1.0).abs() < EPSILON);

        // SLERP: t=0 returns a, t=1 returns b, t=0.5 is spherical midpoint
        let slerp_0 = q1.slerp(q2, 0.0);
        let slerp_1 = q1.slerp(q2, 1.0);
        let slerp_half = q1.slerp(q2, 0.5);
        assert!((slerp_0.x - q1.x).abs() < EPSILON);
        assert!((slerp_1.x - q2.x).abs() < EPSILON);
        assert!((slerp_half.length() - 1.0).abs() < EPSILON);
        // SLERP between identical quats is just the quat
        let slerp_same = q1.slerp(q1, 0.5);
        assert!((slerp_same.x - q1.x).abs() < EPSILON);
        // SLERP between opposite quats (should take shortest path)
        let slerp_opposite = q1.slerp(-q1, 0.5);
        assert!((slerp_opposite.length() - 1.0).abs() < EPSILON);

        // SLERP_APPROX: should be close to slerp for small angles
        let slerp_approx_half = q1.slerp_approx(q3, 0.5);
        let slerp_exact_half = q1.slerp(q3, 0.5);
        assert!((slerp_approx_half.length() - 1.0).abs() < EPSILON);
        assert!((slerp_approx_half.x - slerp_exact_half.x).abs() < 1e-2);

        // NLERP and SLERP should be close for small angles
        let nlerp_small = q1.nlerp(q3, 0.1);
        let slerp_small = q1.slerp(q3, 0.1);
        assert!((nlerp_small.x - slerp_small.x).abs() < 1e-2);

        // NQUAD: midpoint between two nlerps
        let nquad = Quat::nquad(q1, q3, q4, q2, 0.5);
        assert!((nquad.length() - 1.0).abs() < EPSILON);
        // SQUAD: midpoint between two slerps
        let squad = Quat::squad(q1, q3, q4, q2, 0.5);
        assert!((squad.length() - 1.0).abs() < EPSILON);
        // SQUAD_APPROX: midpoint between two slerp_approx
        let squad_approx = Quat::squad_approx(q1, q3, q4, q2, 0.5);
        assert!((squad_approx.length() - 1.0).abs() < EPSILON);

        // Edge: non-normalized input
        let q1_scaled = q1 * 2.0;
        let q2_scaled = q2 * 3.0;
        let nlerp_scaled = q1_scaled.nlerp(q2_scaled, 0.5);
        assert!((nlerp_scaled.length() - 1.0).abs() < EPSILON);
        let slerp_scaled = q1_scaled.slerp(q2_scaled, 0.5);

        // Slerp implementation does not normalize the result.
        // If the inputs are not normalized, the result will not be normalized either.
        // This is standard for most quaternion libraries: slerp expects normalized input.
        assert!(
            (slerp_scaled.normalize().length() - 1.0).abs() < EPSILON,
            "Slerp scaled length: {}",
            slerp_scaled.length()
        );
    }

    #[test]
    fn test_quat_from_mat3_and_to_mat3() {
        // 90 deg rotation around Z
        let angle = PI / 2.0;
        let q =
            Quat::from_axis_angle_radians(Vec3::new(0.0, 0.0, 1.0), gemath::angle::Radians(angle));
        let m = q.to_mat3();
        // Should rotate (1,0,0) to (0,1,0)
        let v = Vec3::new(1.0, 0.0, 0.0);
        let v_rot = m * v;
        assert!((v_rot - Vec3::new(0.0, 1.0, 0.0)).length() < EPSILON);
        // from_mat3 should recover the original quaternion (up to sign)
        let q2 = Quat::from_mat3(&m);
        assert!((q2.normalize().dot(q.normalize())).abs() > 0.999);
    }

    #[test]
    fn test_quat_is_normalized() {
        let q: Quat<(),()> = Quat::from_axis_angle_radians(
            Vec3::new(1.0, 0.0, 0.0),
            gemath::angle::Radians(PI / 2.0),
        );
        assert!(q.is_normalized());
        let q2 = q * 2.0;
        assert!(!q2.is_normalized());
        let q3 = q2.normalize();
        assert!(q3.is_normalized());
    }

    #[test]
    fn test_quat_angle_between() {
        let q1: Quat<(),()> = Quat::IDENTITY;
        let q2 =
            Quat::from_axis_angle_radians(Vec3::new(0.0, 1.0, 0.0), gemath::angle::Radians(PI / 2.0));
        let angle = q1.angle_between(q2);
        assert!((angle - (PI / 2.0)).abs() < EPSILON);
        let angle2 = q2.angle_between(q1);
        assert!((angle2 - (PI / 2.0)).abs() < EPSILON);
        let angle_same = q1.angle_between(q1);
        assert!(angle_same.abs() < EPSILON);
    }

    #[test]
    fn test_quat_rotate_vec3_in_place() {
        let mut v: Vec3<(),()> = Vec3::new(1.0, 0.0, 0.0);
        let q = Quat::from_axis_angle_radians(
            Vec3::new(0.0, 0.0, 1.0),
            gemath::angle::Radians(PI / 2.0),
        );
        q.rotate_vec3(&mut v);
        assert!((v - Vec3::new(0.0, 1.0, 0.0)).length() < EPSILON);
    }

    #[test]
    fn test_quat_ln_and_exp() {
        let q: Quat<(),()> = Quat::from_axis_angle_radians(
            Vec3::new(0.0, 1.0, 0.0),
            gemath::angle::Radians(PI / 3.0),
        );
        let ln_q = q.ln();
        let exp_ln_q = ln_q.exp();
        // exp(ln(q)) should be proportional to q (may differ by sign/scale)
        let qn = q.normalize();
        let expn = exp_ln_q.normalize();
        assert!((qn.dot(expn)).abs() > 0.999);
    }
}

// --- Compile-time (const) tests/examples for Quat<Unit, Space> type-level units and phantom coordinate spaces ---
const _CONST_Q_METERS_WORLD: Quat<Meters, World> = Quat::new(1.0, 2.0, 3.0, 4.0);
const _CONST_Q_PIXELS_SCREEN: Quat<Pixels, Screen> = Quat::new(5.0, 6.0, 7.0, 8.0);
const _CONST_Q_RADIANS_LOCAL: Quat<Radians, Local> = Quat::new(9.0, 10.0, 11.0, 12.0);

const fn _make_quat_meters_world() -> Quat<Meters, World> {
    Quat::new(3.0, 4.0, 5.0, 6.0)
}
const fn _make_quat_pixels_screen() -> Quat<Pixels, Screen> {
    Quat::new(30.0, 40.0, 50.0, 60.0)
}

const _CONST_Q_METERS_WORLD2: Quat<Meters, World> = _make_quat_meters_world();
const _CONST_Q_PIXELS_SCREEN2: Quat<Pixels, Screen> = _make_quat_pixels_screen();

const _: () = {
    // Compile-time assertions for const-everything (units and spaces)
    assert!(_CONST_Q_METERS_WORLD.x == 1.0 && _CONST_Q_METERS_WORLD.y == 2.0 && _CONST_Q_METERS_WORLD.z == 3.0 && _CONST_Q_METERS_WORLD.w == 4.0);
    assert!(_CONST_Q_PIXELS_SCREEN.x == 5.0 && _CONST_Q_PIXELS_SCREEN.y == 6.0 && _CONST_Q_PIXELS_SCREEN.z == 7.0 && _CONST_Q_PIXELS_SCREEN.w == 8.0);
    assert!(_CONST_Q_RADIANS_LOCAL.x == 9.0 && _CONST_Q_RADIANS_LOCAL.y == 10.0 && _CONST_Q_RADIANS_LOCAL.z == 11.0 && _CONST_Q_RADIANS_LOCAL.w == 12.0);
    assert!(_CONST_Q_METERS_WORLD2.x == 3.0 && _CONST_Q_METERS_WORLD2.y == 4.0 && _CONST_Q_METERS_WORLD2.z == 5.0 && _CONST_Q_METERS_WORLD2.w == 6.0);
    assert!(_CONST_Q_PIXELS_SCREEN2.x == 30.0 && _CONST_Q_PIXELS_SCREEN2.y == 40.0 && _CONST_Q_PIXELS_SCREEN2.z == 50.0 && _CONST_Q_PIXELS_SCREEN2.w == 60.0);
    // Compile-time type safety: the following line would fail to compile if uncommented
    // const _FAIL: Quat<Meters, World> = Quat::<Pixels, World>::new(1.0, 2.0, 3.0, 4.0); // error: mismatched types
};