//! # Quat
//!
//! Quaternion type for representing 3D rotations.
//!
//! Quaternions avoid gimbal lock and interpolate smoothly,
//! making them the preferred rotation for game-engine use cases.
//!
//! # Conventions
//!
//! - Euler angles use ZYX order (yaw, pitch, roll)
//!
//! # Example
//!
//! ```rust
//! # use roofy::{Quat, Vec3};
//!
//! let rotation: Quat = Quat::from_axis_angle_degrees(Vec3::Y, 90.0);
//!
//! let forward: Vec3 = Vec3::Z;
//!
//! let rotated: Vec3 = rotation * forward;
//!
//! assert!((rotated.x - 1.0).abs() < 1e-5);
//! assert!(rotated.y.abs() < 1e-5);
//! assert!(rotated.z.abs() < 1e-5);
//! ```
use std::{
    fmt::Display,
    ops::{Index, IndexMut, Mul, Neg},
};

use crate::{Vec3, Vec4};

/// Quaternion representing a 3D rotation, with `f32` components.
/// Same memory layout as `[f32; 4]`.
///
/// Stored as `(x, y, z, w)`, where `(x, y, z)` is the vector part
/// and `w` is the scalar part.
#[derive(Clone, Copy, Debug, PartialEq)]
#[repr(C)]
pub struct Quat {
    /// X component of the quaternion.
    pub x: f32,
    /// Y component of the quaternion.
    pub y: f32,
    /// Z component of the quaternion.
    pub z: f32,
    /// W component of the quaternion.
    pub w: f32,
}

impl Quat {
    /// Threshold for treating a value as zero.
    pub const EPSILON: f32 = 1e-6;
    /// [`EPSILON`](Self::EPSILON) squared.
    pub const EPSILON_SQ: f32 = Self::EPSILON * Self::EPSILON;

    /// Shorthand for writing `Quat::new(f32::NAN, f32::NAN, f32::NAN, f32::NAN)`.
    pub const NAN: Self = Self::new(f32::NAN, f32::NAN, f32::NAN, f32::NAN);

    /// Shorthand for writing `Quat::new(0.0, 0.0, 0.0, 1.0)`.
    pub const IDENTITY: Self = Self::new(0.0, 0.0, 0.0, 1.0);

    /// Creates a `Quat` with `x`, `y`, `z` and `w` components.
    #[inline(always)]
    #[must_use]
    pub const fn new(x: f32, y: f32, z: f32, w: f32) -> Self {
        Self { x, y, z, w }
    }

    /// Creates a `Quat` from `[x, y, z, w]`
    #[inline]
    #[must_use]
    pub const fn from_array(a: [f32; 4]) -> Self {
        Self::new(a[0], a[1], a[2], a[3])
    }

    /// Creates `[x, y, z, w]` from a `Quat`.
    #[inline]
    #[must_use]
    pub const fn to_array(self) -> [f32; 4] {
        [self.x, self.y, self.z, self.w]
    }

    /// Creates a `Quat` from first 4 values of `slice`.
    ///
    /// # Panics
    /// Panics if `slice` is less than 4 elements long.
    #[inline]
    #[must_use]
    pub const fn from_slice(slice: &[f32]) -> Self {
        Self::new(slice[0], slice[1], slice[2], slice[3])
    }

    /// Writes the components of `self` to the first 4 elements of `slice`.
    ///
    /// # Panics
    /// Panics if `slice` is less than 4 elements long.
    #[inline]
    pub const fn write_to_slice(self, slice: &mut [f32]) {
        slice[0] = self.x;
        slice[1] = self.y;
        slice[2] = self.z;
        slice[3] = self.w;
    }

    /// Creates a `Quat` from a `Vec4`.
    #[inline]
    pub const fn from_vec4(v: Vec4) -> Self {
        Self::new(v.x, v.y, v.z, v.w)
    }

    /// Creates a `Quat` from the `angle` (in radians) around the x axis.
    #[inline]
    #[must_use]
    pub fn from_rotation_x(angle: f32) -> Self {
        let (s, c): (f32, f32) = (angle * 0.5).sin_cos();
        Self::new(s, 0.0, 0.0, c)
    }

    /// Creates a `Quat` from the `angle` (in degrees) around the x axis.
    #[inline]
    #[must_use]
    pub fn from_rotation_x_degrees(angle: f32) -> Self {
        Self::from_rotation_x(angle.to_radians())
    }

    /// Creates a `Quat` from the `angle` (in radians) around the y axis.
    #[inline]
    #[must_use]
    pub fn from_rotation_y(angle: f32) -> Self {
        let (s, c): (f32, f32) = (angle * 0.5).sin_cos();
        Self::new(0.0, s, 0.0, c)
    }

    /// Creates a `Quat` from the `angle` (in degrees) around the y axis.
    #[inline]
    #[must_use]
    pub fn from_rotation_y_degrees(angle: f32) -> Self {
        Self::from_rotation_y(angle.to_radians())
    }

    /// Creates a `Quat` from the `angle` (in radians) around the z axis.
    #[inline]
    #[must_use]
    pub fn from_rotation_z(angle: f32) -> Self {
        let (s, c): (f32, f32) = (angle * 0.5).sin_cos();
        Self::new(0.0, 0.0, s, c)
    }

    /// Creates a `Quat` from the `angle` (in degrees) around the z axis.
    #[inline]
    #[must_use]
    pub fn from_rotation_z_degrees(angle: f32) -> Self {
        Self::from_rotation_z(angle.to_radians())
    }

    /// Creates a `Quat` representing a rotation of `angle` (in radians) around `axis`.
    ///
    /// # Preconditions
    /// Expects `axis` to be normalized.
    #[inline]
    #[must_use]
    pub fn from_axis_angle(axis: Vec3, angle: f32) -> Self {
        let (sin, cos): (f32, f32) = (angle * 0.5).sin_cos();
        Self::new(axis.x * sin, axis.y * sin, axis.z * sin, cos)
    }

    /// Returns the axis and angle (in radians) of `self`.
    ///
    /// # Preconditions
    /// Expects `self` to be normalized.
    ///
    /// # Edge cases
    /// If `self` is near-identity, returns `(Vec3::X, 0.0)`.
    #[inline]
    #[must_use]
    pub fn to_axis_angle(self) -> (Vec3, f32) {
        let angle: f32 = self.angle();
        let sin_half: f32 = (angle * 0.5).sin();

        if sin_half.abs() < Self::EPSILON {
            return (Vec3::X, 0.0);
        }

        let inv: f32 = 1.0 / sin_half;
        (Vec3::new(self.x * inv, self.y * inv, self.z * inv), angle)
    }

    /// Creates a `Quat` representing a rotation of `angle` (in degrees) around `axis`.
    ///
    /// # Preconditions
    /// Expects `axis` to be normalized.
    #[inline]
    #[must_use]
    pub fn from_axis_angle_degrees(axis: Vec3, angle: f32) -> Self {
        Self::from_axis_angle(axis, angle.to_radians())
    }

    /// Returns the axis and angle (in degrees) of `self`.
    ///
    /// # Preconditions
    /// Expects `self` to be normalized.
    ///
    /// # Edge cases
    /// If `self` is near-identity, returns `(Vec3::X, 0.0)`.
    #[inline]
    #[must_use]
    pub fn to_axis_angle_degrees(self) -> (Vec3, f32) {
        let (axis, angle): (Vec3, f32) = self.to_axis_angle();
        (axis, angle.to_degrees())
    }

    /// Creates a `Quat` from Euler angles in ZYX order (yaw, pitch, roll) in radians.
    ///
    /// # Convention
    /// Assumes right-handed coordinate system.
    #[inline]
    #[must_use]
    pub fn from_euler_zyx(yaw: f32, pitch: f32, roll: f32) -> Self {
        let (sy, cy): (f32, f32) = (yaw * 0.5).sin_cos();
        let (sp, cp): (f32, f32) = (pitch * 0.5).sin_cos();
        let (sr, cr): (f32, f32) = (roll * 0.5).sin_cos();

        Self::new(
            cy * cp * sr - sy * sp * cr,
            cy * sp * cr + sy * cp * sr,
            sy * cp * cr - cy * sp * sr,
            cy * cp * cr + sy * sp * sr,
        )
    }

    /// Returns a `Vec3` with Euler angles from `self` in ZYX order (yaw, pitch, roll) in radians.
    ///
    /// # Convention
    /// Assumes right-handed coordinate system.
    ///
    /// # Gimbal lock
    /// When `pitch` ≈ `PI / 2`, yaw and roll become indistinguishable,
    /// then yaw is set to `0.0` and all rotation is expressed as roll.
    #[inline]
    #[must_use]
    pub fn to_euler_zyx(self) -> Vec3 {
        let Self { x, y, z, w } = self;

        let sin_pitch: f32 = 2.0 * (w * y - z * x);
        let pitch: f32 = sin_pitch.clamp(-1.0, 1.0).asin();

        if sin_pitch.abs() - 1.0 > Self::EPSILON {
            let yaw: f32 = 0.0;
            let roll: f32 = x.atan2(w) * 2.0;
            return Vec3::new(yaw, pitch, roll);
        }

        let yaw: f32 = (2.0 * (w * z + x * y)).atan2(1.0 - 2.0 * (y * y + z * z));
        let roll: f32 = (2.0 * (w * x + y * z)).atan2(1.0 - 2.0 * (x * x + y * y));

        Vec3::new(yaw, pitch, roll)
    }

    /// Creates a `Quat` from Euler angles in ZYX order (yaw, pitch, roll) in degrees.
    ///
    /// # Convention
    /// Assumes right-handed coordinate system.
    #[inline]
    #[must_use]
    pub fn from_euler_zyx_degrees(yaw: f32, pitch: f32, roll: f32) -> Self {
        Self::from_euler_zyx(yaw.to_radians(), pitch.to_radians(), roll.to_radians())
    }

    /// Returns a `Vec3` with Euler angles from `self` in ZYX order (yaw, pitch, roll) in degrees.
    ///
    /// # Convention
    /// Assumes right-handed coordinate system.
    ///
    /// # Gimbal lock
    /// When `pitch` ≈ `90.0`, yaw and roll become indistinguishable,
    /// then yaw is set to `0.0` and all rotation is expressed as roll.
    #[inline]
    #[must_use]
    pub fn to_euler_degrees_zyx(self) -> Vec3 {
        let e: Vec3 = self.to_euler_zyx();
        Vec3::new(e.x.to_degrees(), e.y.to_degrees(), e.z.to_degrees())
    }

    /// Creates a `Quat` that rotates from `+Z` world-forward to `dir`.
    ///
    /// # Convention
    /// Uses right-handed coordinate system.
    ///
    /// # Preconditions
    /// Expects `dir` to be normalized.
    /// Expects the length of `dir.cross(up)` to not be approximately `0.0` and finite.
    #[must_use]
    pub fn look_to_rh(dir: Vec3, up: Vec3) -> Self {
        let f: Vec3 = dir;
        let r: Vec3 = f.cross(up).normalized();
        let u: Vec3 = f.cross(r);

        let trace: f32 = r.x + u.y + f.z;

        if trace > 0.0 {
            let s: f32 = 0.5 / (trace + 1.0).sqrt();
            Self::new((u.z - f.y) * s, (f.x - r.z) * s, (r.y - u.x) * s, 0.25 / s)
        } else if r.x > u.y && r.x > f.z {
            let s: f32 = 2.0 * (1.0 + r.x - u.y - f.z).sqrt();
            Self::new(0.25 * s, (u.x + r.y) / s, (f.x + r.z) / s, (u.z - f.y) / s)
        } else if u.y > f.z {
            let s: f32 = 2.0 * (1.0 + u.y - r.x - f.z).sqrt();
            Self::new((u.x + r.y) / s, 0.25 * s, (f.y + u.z) / s, (f.x - r.z) / s)
        } else {
            let s: f32 = 2.0 * (1.0 + f.z - r.x - u.y).sqrt();
            Self::new((f.x + r.z) / s, (f.y + u.z) / s, 0.25 * s, (r.y - u.x) / s)
        }
    }

    /// Creates a `Quat` that rotates from `-Z` world-forward to `dir`.
    ///
    /// # Convention
    /// Uses left-handed coordinate system.
    ///
    /// # Preconditions
    /// Expects `dir` to be normalized.
    /// Expects the length of `dir.cross(up)` to not be approximately `0.0` and finite.
    #[inline]
    #[must_use]
    pub fn look_to_lh(dir: Vec3, up: Vec3) -> Self {
        Self::look_to_rh(-dir, up)
    }

    /// Creates a `Quat` that rotates from `+Z` world-forward to face `center` from `eye`.
    ///
    /// # Convention
    /// Uses right-handed coordinate system.
    ///
    /// # Preconditions
    /// Expects the length of `(center - eye)` to not be approximately `0.0` and finite.
    /// Expects the length of `(center - eye).normalized().cross(up)` to not be approximately `0.0` and finite.
    #[inline]
    #[must_use]
    pub fn look_at_rh(eye: Vec3, center: Vec3, up: Vec3) -> Self {
        Self::look_to_rh((center - eye).normalized(), up)
    }

    /// Creates a `Quat` that rotates from `-Z` world-forward to face `center` from `eye`.
    ///
    /// # Convention
    /// Uses left-handed coordinate system.
    ///
    /// # Preconditions
    /// Expects the length of `(center - eye)` to not be approximately `0.0` and finite.
    /// Expects the length of `(center - eye).normalized().cross(up)` to not be approximately `0.0` and finite.
    #[inline]
    #[must_use]
    pub fn look_at_lh(eye: Vec3, center: Vec3, up: Vec3) -> Self {
        Self::look_to_lh((center - eye).normalized(), up)
    }

    /// Creates a `Quat` that rotates `from` to `to` along the shortest arc.
    ///
    /// # Preconditions
    /// Expects `from` and `to` to be normalized.
    #[inline]
    #[must_use]
    pub fn from_to_rotation(from: Vec3, to: Vec3) -> Self {
        let dot: f32 = from.dot(to).clamp(-1.0, 1.0);

        if dot < -1.0 + Self::EPSILON {
            let axis: Vec3 = Vec3::X.cross(from);
            let axis: Vec3 = if axis.length_squared() < Self::EPSILON_SQ {
                Vec3::Y.cross(from)
            } else {
                axis
            };
            return Self::from_axis_angle(axis.normalized(), std::f32::consts::PI);
        }

        let axis: Vec3 = from.cross(to);
        let s: f32 = ((1.0 + dot) * 2.0).sqrt();

        Self::new(axis.x / s, axis.y / s, axis.z / s, s * 0.5).normalized()
    }

    /// Returns the quaternion product of `self` and `rhs`, applying `rhs` first.
    ///
    /// # Preconditions
    /// Expects `self` and `rhs` to be normalized.
    #[inline]
    #[must_use]
    pub const fn mul(self, rhs: Self) -> Self {
        Self::new(
            self.w * rhs.x + self.x * rhs.w + self.y * rhs.z - self.z * rhs.y,
            self.w * rhs.y - self.x * rhs.z + self.y * rhs.w + self.z * rhs.x,
            self.w * rhs.z + self.x * rhs.y - self.y * rhs.x + self.z * rhs.w,
            self.w * rhs.w - self.x * rhs.x - self.y * rhs.y - self.z * rhs.z,
        )
    }

    /// Returns `rhs` rotated by `self`.
    ///
    /// # Preconditions
    /// Expects `self` to be a normalized.
    #[inline]
    #[must_use]
    pub const fn mul_vec3(self, rhs: Vec3) -> Vec3 {
        let u: Vec3 = Vec3::new(self.x, self.y, self.z);
        let t: Vec3 = u.cross(rhs).mul_scalar(2.0);

        rhs.add(t.mul_scalar(self.w)).add(u.cross(t))
    }

    /// Returns `self` with negated components.
    #[inline]
    #[must_use]
    pub const fn neg(self) -> Self {
        Self::new(-self.x, -self.y, -self.z, -self.w)
    }

    /// Returns a reference to the component at `index`.
    ///
    /// `0` maps to `x`, `1` maps to `y`, `2` maps to `z`, `3` maps to `w`.
    ///
    /// # Panics
    /// Panics if `index` is greater than `3`.
    #[inline]
    #[must_use]
    pub const fn index(&self, index: usize) -> &f32 {
        match index {
            0 => &self.x,
            1 => &self.y,
            2 => &self.z,
            3 => &self.w,
            _ => panic!("index out of bounds"),
        }
    }

    /// Returns a mutable reference to the component at `index`.
    ///
    /// `0` maps to `x`, `1` maps to `y`, `2` maps to `z`, `3` maps to `w`.
    ///
    /// # Panics
    /// Panics if `index` is greater than `3`.
    #[inline]
    #[must_use]
    pub const fn index_mut(&mut self, index: usize) -> &mut f32 {
        match index {
            0 => &mut self.x,
            1 => &mut self.y,
            2 => &mut self.z,
            3 => &mut self.w,
            _ => panic!("index out of bounds"),
        }
    }

    /// Returns a reference to the components as `[x, y, z, w]`.
    ///
    /// # Safety
    /// Relies on `Quat` being `#[repr(C)]`
    #[inline]
    #[must_use]
    pub const fn as_ref(&self) -> &[f32; 4] {
        unsafe { &*(self as *const Quat as *const [f32; 4]) }
    }

    /// Returns a mutable reference to the components as `[x, y, z, w]`.
    ///
    /// # Safety
    /// Relies on `Quat` being `#[repr(C)]`
    #[inline]
    #[must_use]
    pub const fn as_mut(&mut self) -> &mut [f32; 4] {
        unsafe { &mut *(self as *mut Quat as *mut [f32; 4]) }
    }

    /// Returns the dot product of `self` and `other`.
    #[inline]
    #[must_use]
    pub const fn dot(self, other: Self) -> f32 {
        self.x * other.x + self.y * other.y + self.z * other.z + self.w * other.w
    }

    /// Returns the length (magnitude) of `self`.
    #[inline]
    #[must_use]
    pub fn length(self) -> f32 {
        self.dot(self).sqrt()
    }

    /// Returns the squared length (magnitude) of `self`.
    #[inline]
    #[must_use]
    pub const fn length_squared(self) -> f32 {
        self.dot(self)
    }

    /// Normalizes `self`.
    ///
    /// # Preconditions
    /// Expects the length of `self` to not be approximately `0.0` and finite.
    #[inline]
    pub fn normalize(&mut self) {
        let len_recip: f32 = self.length().recip();
        self.x *= len_recip;
        self.y *= len_recip;
        self.z *= len_recip;
        self.w *= len_recip;
    }

    /// Returns `self` normalized.
    ///
    /// # Preconditions
    /// Expects the length of `self` to not be approximately `0.0` and finite.
    #[inline]
    #[must_use]
    pub fn normalized(self) -> Self {
        let len_recip: f32 = self.length().recip();
        Self::new(
            self.x * len_recip,
            self.y * len_recip,
            self.z * len_recip,
            self.w * len_recip,
        )
    }

    /// Returns `true` and normalizes `self` if the length of `self`
    /// is not approximately `0.0` and finite, otherwise leaves `self` unchanged
    /// and returns `false`.
    ///
    /// Uses [`EPSILON`](Self::EPSILON) as a threshold.
    #[inline]
    #[must_use]
    pub fn try_normalize(&mut self) -> bool {
        let len_sq: f32 = self.length_squared();

        if len_sq < Self::EPSILON_SQ || !len_sq.is_finite() {
            return false;
        }

        self.normalize();
        true
    }

    /// Returns `self` normalized if the length of `self` is not
    /// approximately `0.0` and finite, otherwise returns `None`.
    ///
    /// Uses [`EPSILON`](Self::EPSILON) as a threshold.
    #[inline]
    #[must_use]
    pub fn try_normalized(self) -> Option<Self> {
        let len_sq: f32 = self.length_squared();

        if len_sq < Self::EPSILON_SQ || !len_sq.is_finite() {
            return None;
        }

        Some(self.normalized())
    }

    /// Returns `self` normalized if the length of `self` is not
    /// approximately `0.0` and finite, otherwise returns `fallback`.
    ///
    /// Uses [`EPSILON`](Self::EPSILON) as a threshold.
    #[inline]
    #[must_use]
    pub fn normalized_or(self, fallback: Self) -> Self {
        self.try_normalized().unwrap_or(fallback)
    }

    /// Returns `self` normalized if the length of `self` is not
    /// approximately `0.0` and finite, otherwise returns [`IDENTITY`](Self::IDENTITY).
    ///
    /// Uses [`EPSILON`](Self::EPSILON) as a threshold.
    #[inline]
    #[must_use]
    pub fn normalized_or_identity(self) -> Self {
        self.normalized_or(Self::IDENTITY)
    }

    /// Returns `true` if `self` approximately has a length of `1.0`.
    ///
    /// Uses [`EPSILON`](Self::EPSILON) as threshold.
    #[inline]
    #[must_use]
    pub const fn is_normalized(self) -> bool {
        (self.length_squared() - 1.0).abs() < Self::EPSILON
    }

    /// Returns a `Quat` that is the conjugate of `self`.
    #[inline]
    #[must_use]
    pub const fn conjugate(self) -> Self {
        Self::new(-self.x, -self.y, -self.z, self.w)
    }

    /// Returns a `Quat` that is the inverse of `self`.
    ///
    /// # Edge cases
    /// If `self` is near-zero, returns [`IDENTITY`](Self::IDENTITY).
    #[inline]
    #[must_use]
    pub const fn inverse(self) -> Self {
        let len_sq: f32 = self.length_squared();

        if len_sq < Self::EPSILON_SQ {
            return Self::IDENTITY;
        }

        let len_sq_recip: f32 = len_sq.recip();
        Self::new(
            -self.x * len_sq_recip,
            -self.y * len_sq_recip,
            -self.z * len_sq_recip,
            self.w * len_sq_recip,
        )
    }

    /// Returns the rotation angle of `self` in radians.
    ///
    /// # Preconditions
    /// Expects `self` to be normalized.
    #[inline]
    #[must_use]
    pub fn angle(self) -> f32 {
        self.w.clamp(-1.0, 1.0).acos() * 2.0
    }

    /// Returns the rotation angle of `self` in degrees.
    ///
    /// # Preconditions
    /// Expects `self` to be normalized.
    #[inline]
    #[must_use]
    pub fn angle_degrees(self) -> f32 {
        self.angle().to_degrees()
    }

    /// Returns the angle  between `self` and `rhs` in radians.
    ///
    /// # Preconditions
    /// Expects `self` and `rhs` to be normalized.
    #[inline]
    #[must_use]
    pub fn angle_between(self, rhs: Self) -> f32 {
        self.dot(rhs).abs().clamp(0.0, 1.0).acos() * 2.0
    }

    /// Returns the angle  between `self` and `rhs` in degrees.
    ///
    /// # Preconditions
    /// Expects `self` and `rhs` to be normalized.
    #[inline]
    #[must_use]
    pub fn angle_between_degrees(self, rhs: Self) -> f32 {
        self.angle_between(rhs).to_degrees()
    }

    /// Returns a `Quat` linearly interpolated between `self` and `to`, with factor `t`, normalized.
    ///
    /// # Preconditions
    /// Expects `self` and `to` to be normalized.
    #[inline]
    #[must_use]
    pub fn lerp(self, mut to: Self, t: f32) -> Self {
        if self.dot(to) < 0.0 {
            to = -to;
        }

        Self::new(
            self.x + (to.x - self.x) * t,
            self.y + (to.y - self.y) * t,
            self.z + (to.z - self.z) * t,
            self.w + (to.w - self.w) * t,
        )
        .normalized()
    }

    /// Returns a `Quat` linearly interpolated between `self` and `to`, with factor `t`, normalized.
    ///
    /// Clamps `t` between `0.0` and `1.0`.
    ///
    /// # Preconditions
    /// Expects `self` and `to` to be normalized.
    #[inline]
    #[must_use]
    pub fn lerp_clamped(self, to: Self, t: f32) -> Self {
        self.lerp(to, t.clamp(0.0, 1.0))
    }

    /// Returns a `Quat` spherically interpolated between `self` and `to`, with factor `t`.
    ///
    /// Falls back to [`lerp`](Self::lerp), when `self` and `to` are nearly identical.
    ///
    /// # Preconditions
    /// Expects `self` and `to` to be normalized.
    #[inline]
    #[must_use]
    pub fn slerp(self, mut to: Self, t: f32) -> Self {
        let mut dot: f32 = self.dot(to).clamp(-1.0, 1.0);

        if dot < 0.0 {
            to = -to;
            dot = -dot;
        }

        if dot > 1.0 - Self::EPSILON {
            return self.lerp(to, t);
        }

        let theta: f32 = dot.acos();
        let sin_theta: f32 = theta.sin();

        let a: f32 = ((1.0 - t) * theta).sin() / sin_theta;
        let b: f32 = (t * theta).sin() / sin_theta;

        Self::new(
            a * self.x + b * to.x,
            a * self.y + b * to.y,
            a * self.z + b * to.z,
            a * self.w + b * to.w,
        )
    }

    /// Returns a `Quat` spherically interpolated between `self` and `to`, with factor `t`.
    ///
    /// Falls back to [`lerp`](Self::lerp), when `self` and `to` are nearly identical.
    ///
    /// Clamps `t` between `0.0` and `1.0`.
    ///
    /// # Preconditions
    /// Expects `self` and `to` to be normalized.
    #[inline]
    #[must_use]
    pub fn slerp_clamped(self, to: Self, t: f32) -> Self {
        self.slerp(to, t.clamp(0.0, 1.0))
    }

    /// Returns a `Quat` rotated from `self` towards `to`, by at most `max_angle` radians.
    ///
    /// # Preconditions
    /// Expects `self` and `to` to be normalized and `max_angle` to be non-negative.
    ///
    /// # Edge cases
    /// If the angle between `self` and `to` is less than or equal to `max_angle`,
    /// returns `to`.
    #[inline]
    #[must_use]
    pub fn rotate_towards(self, to: Self, max_angle: f32) -> Self {
        let angle: f32 = self.angle_between(to);

        if angle < Self::EPSILON || angle <= max_angle {
            return to;
        }

        let t: f32 = max_angle / angle;
        self.slerp(to, t)
    }

    /// Returns `true` if `self` and `other` approximately represent the same rotation.
    ///
    /// Uses [`EPSILON`](Self::EPSILON) as the threshold.
    ///
    /// # Preconditions
    /// Expects `self` and `other` to be normalized.
    #[inline]
    #[must_use]
    pub const fn approx_eq(self, other: Self) -> bool {
        self.dot(other).abs() > 1.0 - Self::EPSILON
    }

    /// Returns `true` if all components of `self` are finite.
    #[inline]
    #[must_use]
    pub const fn is_finite(self) -> bool {
        self.x.is_finite() && self.y.is_finite() && self.z.is_finite() && self.w.is_finite()
    }

    /// Returns `true` if any components of `self` is `NaN`.
    #[inline]
    #[must_use]
    pub const fn is_nan(self) -> bool {
        self.x.is_nan() || self.y.is_nan() || self.z.is_nan() || self.w.is_nan()
    }
}

impl Mul for Quat {
    type Output = Self;

    /// Returns the quaternion product of `self` and `rhs`, applying `rhs` first.
    ///
    /// Equivalent to [`mul`](Self::mul()).
    ///
    /// # Preconditions
    /// Expects `self` and `rhs` to be normalized.
    #[inline]
    fn mul(self, rhs: Self) -> Self::Output {
        self.mul(rhs)
    }
}

impl Mul<Vec3> for Quat {
    type Output = Vec3;

    /// Returns `rhs` rotated by `self`.
    ///
    /// Equivalent to [`mul_vec3`](Self::mul_vec3()).
    ///
    /// # Preconditions
    /// Expects `self` to be a normalized.
    #[inline]
    fn mul(self, rhs: Vec3) -> Self::Output {
        self.mul_vec3(rhs)
    }
}

impl Neg for Quat {
    type Output = Self;

    /// Returns `self` with negated components.
    ///
    /// Equivalent to [`neg`](Self::neg()).
    #[inline]
    fn neg(self) -> Self::Output {
        self.neg()
    }
}

impl Index<usize> for Quat {
    type Output = f32;

    /// Returns a reference to the component at `index`.
    ///
    /// `0` maps to `x`, `1` maps to `y`, `2` maps to `z`, `3` maps to `w`.
    ///
    /// Equivalent to [`index`](Self::index()).
    ///
    /// # Panics
    /// Panics if `index` is greater than `3`.
    #[inline]
    fn index(&self, index: usize) -> &Self::Output {
        self.index(index)
    }
}

impl IndexMut<usize> for Quat {
    /// Returns a mutable reference to the component at `index`.
    ///
    /// `0` maps to `x`, `1` maps to `y`, `2` maps to `z`, `3` maps to `w`.
    ///
    /// Equivalent to [`index_mut`](Self::index_mut()).
    ///
    /// # Panics
    /// Panics if `index` is greater than `3`.
    #[inline]
    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
        self.index_mut(index)
    }
}

impl AsRef<[f32; 4]> for Quat {
    /// Returns a reference to the components as `[x, y, z, w]`.
    ///
    /// Equivalent to [`as_ref`](Self::as_ref()).
    ///
    /// # Safety
    /// Relies on `Quat` being `#[repr(C)]`
    #[inline]
    fn as_ref(&self) -> &[f32; 4] {
        self.as_ref()
    }
}

impl AsMut<[f32; 4]> for Quat {
    /// Returns a mutable reference to the components as `[x, y, z, w]`.
    ///
    /// Equivalent to [`as_mut`](Self::as_mut()).
    ///
    /// # Safety
    /// Relies on `Quat` being `#[repr(C)]`
    #[inline]
    fn as_mut(&mut self) -> &mut [f32; 4] {
        self.as_mut()
    }
}

impl Display for Quat {
    /// Formats `self` as `(x, y, z, w)`.
    ///
    /// Respects the precision of `f`.
    #[inline]
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        if let Some(p) = f.precision() {
            write!(
                f,
                "({:.*}, {:.*}, {:.*}, {:.*})",
                p, self.x, p, self.y, p, self.z, p, self.w
            )
        } else {
            write!(f, "({}, {}, {}, {})", self.x, self.y, self.z, self.w)
        }
    }
}

impl Default for Quat {
    /// Default is [`IDENTITY`](Self::IDENTITY).
    #[inline]
    fn default() -> Self {
        Self::IDENTITY
    }
}

impl From<[f32; 4]> for Quat {
    /// Creates a `Quat` from `[x, y, z, w]`.
    ///
    /// Equivalent to [`from_array`](Self::from_array()).
    #[inline]
    fn from(v: [f32; 4]) -> Self {
        Self::from_array(v)
    }
}

impl From<Quat> for [f32; 4] {
    /// Creates `[x, y, z, w]` from a `Quat`.
    ///
    /// Equivalent to [`to_array`](Quat::to_array()).
    #[inline]
    fn from(v: Quat) -> Self {
        v.to_array()
    }
}

impl From<&[f32]> for Quat {
    /// Creates a `Quat` from first 4 values of `v`.
    ///
    /// Equivalent to [`from_slice`](Self::from_slice()).
    ///
    /// # Panics
    /// Panics if `v` is less than 4 elements long.
    #[inline]
    fn from(v: &[f32]) -> Self {
        Self::from_slice(v)
    }
}

impl From<Vec4> for Quat {
    /// Creates a `Quat` from a `Vec4`.
    #[inline]
    fn from(v: Vec4) -> Self {
        Self::new(v.x, v.y, v.z, v.w)
    }
}

#[cfg(feature = "serde")]
impl serde::Serialize for Quat {
    /// Serializes as `[x, y, z, w]`.
    ///
    /// # Format
    /// Encoded as a flat 4-element array. Component order matches [`to_array`](Self::to_array()).
    fn serialize<S: serde::Serializer>(&self, serializer: S) -> Result<S::Ok, S::Error> {
        self.to_array().serialize(serializer)
    }
}

#[cfg(feature = "serde")]
impl<'de> serde::Deserialize<'de> for Quat {
    /// Deserializes from `[x, y, z, w]`.
    ///
    /// # Format
    /// Decoded as a flat 4-element array. Component order matches [`from_array`](Self::from_array()).
    ///
    /// # Errors
    /// Returns an error if array is not exactly 4 elements.
    fn deserialize<D: serde::Deserializer<'de>>(deserializer: D) -> Result<Self, D::Error> {
        <[f32; 4]>::deserialize(deserializer).map(Self::from_array)
    }
}

#[cfg(feature = "bytemuck")]
unsafe impl bytemuck::Zeroable for Quat {}

#[cfg(feature = "bytemuck")]
unsafe impl bytemuck::Pod for Quat {}
