glam/f32/sse2/

quat.rs

// Generated from quat.rs.tera template. Edit the template, not the generated file.

use crate::{
    Mat3, Mat3A, Mat4, Vec2, Vec3, Vec3A, Vec4,
    euler::{EulerRot, FromEuler, ToEuler},
    f32::math,
    sse2::*,
};

#[cfg(feature = "f64")]
use crate::DQuat;

#[cfg(target_arch = "x86")]
use core::arch::x86::*;
#[cfg(target_arch = "x86_64")]
use core::arch::x86_64::*;

use core::fmt;
use core::iter::{Product, Sum};
use core::ops::{
    Add, AddAssign, Deref, DerefMut, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign,
};

#[cfg(feature = "zerocopy-08")]
use zerocopy_derive_08::*;

#[repr(C)]
union UnionCast {
    a: [f32; 4],
    v: Quat,
}

/// Creates a quaternion from `x`, `y`, `z` and `w` values.
///
/// This should generally not be called manually unless you know what you are doing. Use
/// one of the other constructors instead such as `identity` or `from_axis_angle`.
#[inline]
#[must_use]
pub const fn quat(x: f32, y: f32, z: f32, w: f32) -> Quat {
    Quat::from_xyzw(x, y, z, w)
}

/// A quaternion representing an orientation.
///
/// This quaternion is intended to be of unit length but may denormalize due to
/// floating point "error creep" which can occur when successive quaternion
/// operations are applied.
///
/// SIMD vector types are used for storage on supported platforms.
///
/// This type is 16 byte aligned.
#[derive(Clone, Copy)]
#[cfg_attr(feature = "bytemuck", derive(bytemuck::Pod, bytemuck::Zeroable))]
#[cfg_attr(
    feature = "zerocopy-08",
    derive(FromBytes, Immutable, IntoBytes, KnownLayout)
)]
#[repr(transparent)]
pub struct Quat(pub(crate) __m128);

impl Quat {
    /// All zeros.
    const ZERO: Self = Self::from_array([0.0; 4]);

    /// The identity quaternion. Corresponds to no rotation.
    pub const IDENTITY: Self = Self::from_xyzw(0.0, 0.0, 0.0, 1.0);

    /// All NANs.
    pub const NAN: Self = Self::from_array([f32::NAN; 4]);

    /// Creates a new rotation quaternion.
    ///
    /// This should generally not be called manually unless you know what you are doing.
    /// Use one of the other constructors instead such as `identity` or `from_axis_angle`.
    ///
    /// `from_xyzw` is mostly used by unit tests and `serde` deserialization.
    ///
    /// # Preconditions
    ///
    /// This function does not check if the input is normalized, it is up to the user to
    /// provide normalized input or to normalized the resulting quaternion.
    #[inline(always)]
    #[must_use]
    pub const fn from_xyzw(x: f32, y: f32, z: f32, w: f32) -> Self {
        unsafe { UnionCast { a: [x, y, z, w] }.v }
    }

    /// Creates a rotation quaternion from an array.
    ///
    /// # Preconditions
    ///
    /// This function does not check if the input is normalized, it is up to the user to
    /// provide normalized input or to normalized the resulting quaternion.
    #[inline]
    #[must_use]
    pub const fn from_array(a: [f32; 4]) -> Self {
        Self::from_xyzw(a[0], a[1], a[2], a[3])
    }

    /// Creates a new rotation quaternion from a 4D vector.
    ///
    /// # Preconditions
    ///
    /// This function does not check if the input is normalized, it is up to the user to
    /// provide normalized input or to normalized the resulting quaternion.
    #[inline]
    #[must_use]
    pub const fn from_vec4(v: Vec4) -> Self {
        Self(v.0)
    }

    /// Creates a rotation quaternion from a slice.
    ///
    /// # Preconditions
    ///
    /// This function does not check if the input is normalized, it is up to the user to
    /// provide normalized input or to normalized the resulting quaternion.
    ///
    /// # Panics
    ///
    /// Panics if `slice` length is less than 4.
    #[inline]
    #[must_use]
    #[track_caller]
    pub fn from_slice(slice: &[f32]) -> Self {
        assert!(slice.len() >= 4);
        Self(unsafe { _mm_loadu_ps(slice.as_ptr()) })
    }

    /// Writes the quaternion to an unaligned slice.
    ///
    /// # Panics
    ///
    /// Panics if `slice` length is less than 4.
    #[inline]
    #[track_caller]
    pub fn write_to_slice(self, slice: &mut [f32]) {
        assert!(slice.len() >= 4);
        unsafe { _mm_storeu_ps(slice.as_mut_ptr(), self.0) }
    }

    /// Create a quaternion for a normalized rotation `axis` and `angle` (in radians).
    ///
    /// The axis must be a unit vector.
    ///
    /// # Panics
    ///
    /// Will panic if `axis` is not normalized when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn from_axis_angle(axis: Vec3, angle: f32) -> Self {
        glam_assert!(axis.is_normalized());
        let (s, c) = math::sin_cos(angle * 0.5);
        let v = axis * s;
        Self::from_xyzw(v.x, v.y, v.z, c)
    }

    /// Create a quaternion that rotates `v.length()` radians around `v.normalize()`.
    ///
    /// `from_scaled_axis(Vec3::ZERO)` results in the identity quaternion.
    ///
    /// # Panics
    ///
    /// Will panic if `v` is not finite when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn from_scaled_axis(v: Vec3) -> Self {
        let length = v.length();
        if length == 0.0 {
            Self::IDENTITY
        } else {
            Self::from_axis_angle(v / length, length)
        }
    }

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

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

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

    /// Creates a quaternion from the given Euler rotation sequence and the angles (in radians).
    #[inline]
    #[must_use]
    pub fn from_euler(euler: EulerRot, a: f32, b: f32, c: f32) -> Self {
        Self::from_euler_angles(euler, a, b, c)
    }

    /// From the columns of a 3x3 rotation matrix.
    ///
    /// Note if the input axes contain scales, shears, or other non-rotation transformations then
    /// the output of this function is ill-defined.
    ///
    /// # Panics
    ///
    /// Will panic if any axis is not normalized when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn from_rotation_axes(x_axis: Vec3, y_axis: Vec3, z_axis: Vec3) -> Self {
        glam_assert!(x_axis.is_normalized() && y_axis.is_normalized() && z_axis.is_normalized());
        // Based on https://github.com/microsoft/DirectXMath `XMQuaternionRotationMatrix`
        let (m00, m01, m02) = x_axis.into();
        let (m10, m11, m12) = y_axis.into();
        let (m20, m21, m22) = z_axis.into();
        if m22 <= 0.0 {
            // x^2 + y^2 >= z^2 + w^2
            let dif10 = m11 - m00;
            let omm22 = 1.0 - m22;
            if dif10 <= 0.0 {
                // x^2 >= y^2
                let four_xsq = omm22 - dif10;
                let inv4x = 0.5 / math::sqrt(four_xsq);
                Self::from_xyzw(
                    four_xsq * inv4x,
                    (m01 + m10) * inv4x,
                    (m02 + m20) * inv4x,
                    (m12 - m21) * inv4x,
                )
            } else {
                // y^2 >= x^2
                let four_ysq = omm22 + dif10;
                let inv4y = 0.5 / math::sqrt(four_ysq);
                Self::from_xyzw(
                    (m01 + m10) * inv4y,
                    four_ysq * inv4y,
                    (m12 + m21) * inv4y,
                    (m20 - m02) * inv4y,
                )
            }
        } else {
            // z^2 + w^2 >= x^2 + y^2
            let sum10 = m11 + m00;
            let opm22 = 1.0 + m22;
            if sum10 <= 0.0 {
                // z^2 >= w^2
                let four_zsq = opm22 - sum10;
                let inv4z = 0.5 / math::sqrt(four_zsq);
                Self::from_xyzw(
                    (m02 + m20) * inv4z,
                    (m12 + m21) * inv4z,
                    four_zsq * inv4z,
                    (m01 - m10) * inv4z,
                )
            } else {
                // w^2 >= z^2
                let four_wsq = opm22 + sum10;
                let inv4w = 0.5 / math::sqrt(four_wsq);
                Self::from_xyzw(
                    (m12 - m21) * inv4w,
                    (m20 - m02) * inv4w,
                    (m01 - m10) * inv4w,
                    four_wsq * inv4w,
                )
            }
        }
    }

    /// Creates a quaternion from a 3x3 rotation matrix.
    ///
    /// Note if the input matrix contain scales, shears, or other non-rotation transformations then
    /// the resulting quaternion will be ill-defined.
    ///
    /// # Panics
    ///
    /// Will panic if any input matrix column is not normalized when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn from_mat3(mat: &Mat3) -> Self {
        Self::from_rotation_axes(mat.x_axis, mat.y_axis, mat.z_axis)
    }

    /// Creates a quaternion from a 3x3 SIMD aligned rotation matrix.
    ///
    /// Note if the input matrix contain scales, shears, or other non-rotation transformations then
    /// the resulting quaternion will be ill-defined.
    ///
    /// # Panics
    ///
    /// Will panic if any input matrix column is not normalized when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn from_mat3a(mat: &Mat3A) -> Self {
        Self::from_rotation_axes(mat.x_axis.into(), mat.y_axis.into(), mat.z_axis.into())
    }

    /// Creates a quaternion from the upper 3x3 rotation matrix inside a homogeneous 4x4 matrix.
    ///
    /// Note if the upper 3x3 matrix contain scales, shears, or other non-rotation transformations
    /// then the resulting quaternion will be ill-defined.
    ///
    /// # Panics
    ///
    /// Will panic if any column of the upper 3x3 rotation matrix is not normalized when
    /// `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn from_mat4(mat: &Mat4) -> Self {
        Self::from_rotation_axes(
            mat.x_axis.truncate(),
            mat.y_axis.truncate(),
            mat.z_axis.truncate(),
        )
    }

    /// Gets the minimal rotation for transforming `from` to `to`.  The rotation is in the
    /// plane spanned by the two vectors.  Will rotate at most 180 degrees.
    ///
    /// The inputs must be unit vectors.
    ///
    /// `from_rotation_arc(from, to) * from ≈ to`.
    ///
    /// For near-singular cases (from≈to and from≈-to) the current implementation
    /// is only accurate to about 0.001 (for `f32`).
    ///
    /// # Panics
    ///
    /// Will panic if `from` or `to` are not normalized when `glam_assert` is enabled.
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn from_rotation_arc(from: Vec3, to: Vec3) -> Self {
        glam_assert!(from.is_normalized());
        glam_assert!(to.is_normalized());

        const ONE_MINUS_EPS: f32 = 1.0 - 2.0 * f32::EPSILON;
        let dot = from.dot(to);
        if dot > ONE_MINUS_EPS {
            // 0° singularity: from ≈ to
            Self::IDENTITY
        } else if dot < -ONE_MINUS_EPS {
            // 180° singularity: from ≈ -to
            use core::f32::consts::PI; // half a turn = 𝛕/2 = 180°
            Self::from_axis_angle(from.any_orthonormal_vector(), PI)
        } else {
            let c = from.cross(to);
            Self::from_xyzw(c.x, c.y, c.z, 1.0 + dot).normalize()
        }
    }

    /// Gets the minimal rotation for transforming `from` to either `to` or `-to`.  This means
    /// that the resulting quaternion will rotate `from` so that it is colinear with `to`.
    ///
    /// The rotation is in the plane spanned by the two vectors.  Will rotate at most 90
    /// degrees.
    ///
    /// The inputs must be unit vectors.
    ///
    /// `to.dot(from_rotation_arc_colinear(from, to) * from).abs() ≈ 1`.
    ///
    /// # Panics
    ///
    /// Will panic if `from` or `to` are not normalized when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn from_rotation_arc_colinear(from: Vec3, to: Vec3) -> Self {
        if from.dot(to) < 0.0 {
            Self::from_rotation_arc(from, -to)
        } else {
            Self::from_rotation_arc(from, to)
        }
    }

    /// Gets the minimal rotation for transforming `from` to `to`.  The resulting rotation is
    /// around the z axis. Will rotate at most 180 degrees.
    ///
    /// The inputs must be unit vectors.
    ///
    /// `from_rotation_arc_2d(from, to) * from ≈ to`.
    ///
    /// For near-singular cases (from≈to and from≈-to) the current implementation
    /// is only accurate to about 0.001 (for `f32`).
    ///
    /// # Panics
    ///
    /// Will panic if `from` or `to` are not normalized when `glam_assert` is enabled.
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn from_rotation_arc_2d(from: Vec2, to: Vec2) -> Self {
        glam_assert!(from.is_normalized());
        glam_assert!(to.is_normalized());

        const ONE_MINUS_EPSILON: f32 = 1.0 - 2.0 * f32::EPSILON;
        let dot = from.dot(to);
        if dot > ONE_MINUS_EPSILON {
            // 0° singularity: from ≈ to
            Self::IDENTITY
        } else if dot < -ONE_MINUS_EPSILON {
            // 180° singularity: from ≈ -to
            const COS_FRAC_PI_2: f32 = 0.0;
            const SIN_FRAC_PI_2: f32 = 1.0;
            // rotation around z by PI radians
            Self::from_xyzw(0.0, 0.0, SIN_FRAC_PI_2, COS_FRAC_PI_2)
        } else {
            // vector3 cross where z=0
            let z = from.x * to.y - to.x * from.y;
            let w = 1.0 + dot;
            // calculate length with x=0 and y=0 to normalize
            let len_rcp = 1.0 / math::sqrt(z * z + w * w);
            Self::from_xyzw(0.0, 0.0, z * len_rcp, w * len_rcp)
        }
    }

    /// Returns the rotation axis (normalized) and angle (in radians) of `self`.
    #[inline]
    #[must_use]
    pub fn to_axis_angle(self) -> (Vec3, f32) {
        const EPSILON: f32 = 1.0e-8;
        let v = Vec3::new(self.x, self.y, self.z);
        let length = v.length();
        if length >= EPSILON {
            let angle = 2.0 * math::atan2(length, self.w);
            let axis = v / length;
            (axis, angle)
        } else {
            (Vec3::X, 0.0)
        }
    }

    /// Returns the rotation axis scaled by the rotation in radians.
    #[inline]
    #[must_use]
    pub fn to_scaled_axis(self) -> Vec3 {
        let (axis, angle) = self.to_axis_angle();
        axis * angle
    }

    /// Returns the rotation angles for the given euler rotation sequence.
    ///
    /// # Panics
    ///
    /// Will panic if `self` is not normalized when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn to_euler(self, order: EulerRot) -> (f32, f32, f32) {
        Mat3::from_quat(self).to_euler_angles(order)
    }

    /// Converts `self` to `[x, y, z, w]`
    #[inline]
    #[must_use]
    pub const fn to_array(&self) -> [f32; 4] {
        unsafe { *(self as *const Self as *const [f32; 4]) }
    }

    /// Returns the vector part of the quaternion.
    #[inline]
    #[must_use]
    pub fn xyz(self) -> Vec3 {
        Vec3::new(self.x, self.y, self.z)
    }

    /// Returns the quaternion conjugate of `self`. For a unit quaternion the
    /// conjugate is also the inverse.
    #[inline]
    #[must_use]
    pub fn conjugate(self) -> Self {
        const SIGN: __m128 = m128_from_f32x4([-0.0, -0.0, -0.0, 0.0]);
        Self(unsafe { _mm_xor_ps(self.0, SIGN) })
    }

    /// Returns the inverse of a normalized quaternion.
    ///
    /// Typically quaternion inverse returns the conjugate of a normalized quaternion.
    /// Because `self` is assumed to already be unit length this method *does not* normalize
    /// before returning the conjugate.
    ///
    /// # Panics
    ///
    /// Will panic if `self` is not normalized when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn inverse(self) -> Self {
        glam_assert!(self.is_normalized());
        self.conjugate()
    }

    /// Computes the dot product of `self` and `rhs`. The dot product is
    /// equal to the cosine of the angle between two quaternion rotations.
    #[inline]
    #[must_use]
    pub fn dot(self, rhs: Self) -> f32 {
        Vec4::from(self).dot(Vec4::from(rhs))
    }

    /// Computes the length of `self`.
    #[doc(alias = "magnitude")]
    #[inline]
    #[must_use]
    pub fn length(self) -> f32 {
        Vec4::from(self).length()
    }

    /// Computes the squared length of `self`.
    ///
    /// This is generally faster than `length()` as it avoids a square
    /// root operation.
    #[doc(alias = "magnitude2")]
    #[inline]
    #[must_use]
    pub fn length_squared(self) -> f32 {
        Vec4::from(self).length_squared()
    }

    /// Computes `1.0 / length()`.
    ///
    /// For valid results, `self` must _not_ be of length zero.
    #[inline]
    #[must_use]
    pub fn length_recip(self) -> f32 {
        Vec4::from(self).length_recip()
    }

    /// Returns `self` normalized to length 1.0.
    ///
    /// For valid results, `self` must _not_ be of length zero.
    ///
    /// # Panics
    ///
    /// Will panic if `self` is zero length when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn normalize(self) -> Self {
        Self::from_vec4(Vec4::from(self).normalize())
    }

    /// Returns `true` if, and only if, all elements are finite.
    /// If any element is either `NaN`, positive or negative infinity, this will return `false`.
    #[inline]
    #[must_use]
    pub fn is_finite(self) -> bool {
        Vec4::from(self).is_finite()
    }

    /// Returns `true` if any elements are `NAN`.
    #[inline]
    #[must_use]
    pub fn is_nan(self) -> bool {
        Vec4::from(self).is_nan()
    }

    /// Returns whether `self` of length `1.0` or not.
    ///
    /// Uses a precision threshold of `1e-6`.
    #[inline]
    #[must_use]
    pub fn is_normalized(self) -> bool {
        Vec4::from(self).is_normalized()
    }

    /// Returns `true` if `self` represents a rotation near the identity.
    #[inline]
    #[must_use]
    pub fn is_near_identity(self) -> bool {
        // Based on https://github.com/nfrechette/rtm `rtm::quat_near_identity`
        // The shortest rotation angle is `2 * acos(abs(w))`. Since `acos` is monotonically
        // decreasing, comparing `abs(w)` to the cosine threshold avoids calculating the angle.
        // Equivalent to an angular threshold of `(1.0 - 1e-6).acos() * 2.0`.
        const THRESHOLD: f32 = 1.0 - 1e-6;
        self.w.abs() > THRESHOLD
    }

    /// Returns the angle (in radians) for the minimal rotation between two quaternions
    /// in the range `[0, +π]`.
    ///
    /// Both quaternions must be normalized.
    ///
    /// # Panics
    ///
    /// Will panic if `self` or `rhs` are not normalized when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn angle_between(self, rhs: Self) -> f32 {
        glam_assert!(self.is_normalized() && rhs.is_normalized());
        math::acos_approx(self.dot(rhs).abs()) * 2.0
    }

    /// Rotates towards `rhs` up to `max_angle` (in radians).
    ///
    /// When `max_angle` is `0.0`, the result will be equal to `self`. When `max_angle` is equal to
    /// `self.angle_between(rhs)`, the result will be equal to `rhs`. If `max_angle` is negative,
    /// rotates towards the exact opposite of `rhs`. Will not go past the target.
    ///
    /// Both quaternions must be normalized.
    ///
    /// # Panics
    ///
    /// Will panic if `self` or `rhs` are not normalized when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn rotate_towards(self, rhs: Self, max_angle: f32) -> Self {
        glam_assert!(self.is_normalized() && rhs.is_normalized());
        let angle = self.angle_between(rhs);
        if angle <= 1e-4 {
            return rhs;
        }
        if max_angle < 0.0 {
            // The rotation opposite `rhs` is the component of `self` orthogonal to it.
            let target = (self - rhs * self.dot(rhs)).normalize();
            let angle = self.angle_between(target);
            if angle <= 1e-4 {
                return target;
            }
            return self.slerp(target, (-max_angle / angle).clamp(0.0, 1.0));
        }
        self.slerp(rhs, (max_angle / angle).clamp(0.0, 1.0))
    }

    /// Returns true if the absolute difference of all elements between `self` and `rhs`
    /// is less than or equal to `max_abs_diff`.
    ///
    /// This can be used to compare if two quaternions contain similar elements. It works
    /// best when comparing with a known value. The `max_abs_diff` that should be used used
    /// depends on the values being compared against.
    ///
    /// For more see
    /// [comparing floating point numbers](https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/).
    #[inline]
    #[must_use]
    pub fn abs_diff_eq(self, rhs: Self, max_abs_diff: f32) -> bool {
        Vec4::from(self).abs_diff_eq(Vec4::from(rhs), max_abs_diff)
    }

    #[inline(always)]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    fn lerp_impl(self, end: Self, s: f32) -> Self {
        (self * (1.0 - s) + end * s).normalize()
    }

    /// Performs a linear interpolation between `self` and `rhs` based on
    /// the value `s`, using the form `self * (1.0 - s) + end * s` before normalizing.
    ///
    /// When `s` is `0.0`, the result will be equal to `self`.  When `s`
    /// is `1.0`, the result will be equal to `rhs`.
    ///
    /// This interpolates linearly between the two rotations and does not rotate at a constant
    /// angular velocity; see [`slerp`](Self::slerp) if that is required.
    ///
    /// # Panics
    ///
    /// Will panic if `self` or `end` are not normalized when `glam_assert` is enabled.
    #[doc(alias = "mix")]
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn lerp(self, end: Self, s: f32) -> Self {
        glam_assert!(self.is_normalized());
        glam_assert!(end.is_normalized());

        const NEG_ZERO: __m128 = m128_from_f32x4([-0.0; 4]);
        unsafe {
            let dot = dot4_into_m128(self.0, end.0);
            // Calculate the bias, if the dot product is positive or zero, there is no bias
            // but if it is negative, we want to flip the 'end' rotation XYZW components
            let bias = _mm_and_ps(dot, NEG_ZERO);
            self.lerp_impl(Self(_mm_xor_ps(end.0, bias)), s)
        }
    }

    #[inline(always)]
    #[must_use]
    fn slerp_impl(self, end: Self, dot: f32, s: f32) -> Self {
        let theta = math::acos_approx(dot);

        let x = 1.0 - s;
        let y = s;
        let z = 1.0;

        unsafe {
            let tmp = _mm_mul_ps(_mm_set_ps1(theta), _mm_set_ps(0.0, z, y, x));
            let tmp = m128_sin(tmp);

            let scale1 = _mm_shuffle_ps(tmp, tmp, 0b00_00_00_00);
            let scale2 = _mm_shuffle_ps(tmp, tmp, 0b01_01_01_01);
            let theta_sin = _mm_shuffle_ps(tmp, tmp, 0b10_10_10_10);

            Self(_mm_div_ps(
                m128_mul_add(end.0, scale2, _mm_mul_ps(self.0, scale1)),
                theta_sin,
            ))
        }
    }

    /// Performs a spherical linear interpolation between `self` and `end`
    /// based on the value `s`.
    ///
    /// When `s` is `0.0`, the result will be equal to `self`.  When `s`
    /// is `1.0`, the result will be equal to `end`.
    ///
    /// # Panics
    ///
    /// Will panic if `self` or `end` are not normalized when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn slerp(self, mut end: Self, s: f32) -> Self {
        // http://number-none.com/product/Understanding%20Slerp,%20Then%20Not%20Using%20It/
        glam_assert!(self.is_normalized());
        glam_assert!(end.is_normalized());

        // Note that a rotation can be represented by two quaternions: `q` and
        // `-q`. The slerp path between `q` and `end` will be different from the
        // path between `-q` and `end`. One path will take the long way around and
        // one will take the short way. In order to correct for this, the `dot`
        // product between `self` and `end` should be positive. If the `dot`
        // product is negative, slerp between `self` and `-end`.
        let mut dot = self.dot(end);
        if dot < 0.0 {
            end = -end;
            dot = -dot;
        }

        const DOT_THRESHOLD: f32 = 1.0 - f32::EPSILON;
        if dot > DOT_THRESHOLD {
            // if above threshold perform linear interpolation to avoid divide by zero
            self.lerp_impl(end, s)
        } else {
            self.slerp_impl(end, dot, s)
        }
    }

    /// Performs a spherical linear interpolation between `self` and `end` based on the value `s`,
    /// preserving the rotation direction.
    ///
    /// When `s` is `0.0`, the result will be equal to `self`.  When `s` is `1.0`, the result will
    /// be equal to `end`.
    ///
    /// When the dot product of `self` and `end` is negative, the standard [`slerp`](Self::slerp)
    /// will flip the end quaternion to take the shortest path, while this method will take the
    /// longer arc. This is useful when the intended rotation direction must be preserved.
    ///
    /// # Panics
    ///
    /// Will panic if `self` or `end` are not normalized when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn slerp_long(self, end: Self, s: f32) -> Self {
        glam_assert!(self.is_normalized());
        glam_assert!(end.is_normalized());

        let dot = self.dot(end);

        const DOT_THRESHOLD: f32 = 1.0 - f32::EPSILON;
        if dot.abs() > DOT_THRESHOLD {
            // if above threshold perform linear interpolation to avoid divide by zero
            self.lerp_impl(end, s)
        } else {
            self.slerp_impl(end, dot, s)
        }
    }

    /// Multiplies a quaternion and a 3D vector, returning the rotated vector.
    ///
    /// # Panics
    ///
    /// Will panic if `self` is not normalized when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn mul_vec3(self, rhs: Vec3) -> Vec3 {
        glam_assert!(self.is_normalized());

        self.mul_vec3a(rhs.into()).into()
    }

    /// Multiplies two quaternions. If they each represent a rotation, the result will
    /// represent the combined rotation.
    ///
    /// Note that due to floating point rounding the result may not be perfectly normalized.
    #[inline]
    #[must_use]
    pub fn mul_quat(self, rhs: Self) -> Self {
        // Based on https://github.com/nfrechette/rtm `rtm::quat_mul`
        const CONTROL_WZYX: __m128 = m128_from_f32x4([1.0, -1.0, 1.0, -1.0]);
        const CONTROL_ZWXY: __m128 = m128_from_f32x4([1.0, 1.0, -1.0, -1.0]);
        const CONTROL_YXWZ: __m128 = m128_from_f32x4([-1.0, 1.0, 1.0, -1.0]);

        let lhs = self.0;
        let rhs = rhs.0;

        unsafe {
            let r_xxxx = _mm_shuffle_ps(lhs, lhs, 0b00_00_00_00);
            let r_yyyy = _mm_shuffle_ps(lhs, lhs, 0b01_01_01_01);
            let r_zzzz = _mm_shuffle_ps(lhs, lhs, 0b10_10_10_10);
            let r_wwww = _mm_shuffle_ps(lhs, lhs, 0b11_11_11_11);

            let lxrw_lyrw_lzrw_lwrw = _mm_mul_ps(r_wwww, rhs);
            let l_wzyx = _mm_shuffle_ps(rhs, rhs, 0b00_01_10_11);

            let lwrx_lzrx_lyrx_lxrx = _mm_mul_ps(r_xxxx, l_wzyx);
            let l_zwxy = _mm_shuffle_ps(l_wzyx, l_wzyx, 0b10_11_00_01);

            let lzry_lwry_lxry_lyry = _mm_mul_ps(r_yyyy, l_zwxy);
            let l_yxwz = _mm_shuffle_ps(l_zwxy, l_zwxy, 0b00_01_10_11);

            let lzry_lwry_nlxry_nlyry = _mm_mul_ps(lzry_lwry_lxry_lyry, CONTROL_ZWXY);

            let lyrz_lxrz_lwrz_lzrz = _mm_mul_ps(r_zzzz, l_yxwz);
            let result0 = m128_mul_add(lwrx_lzrx_lyrx_lxrx, CONTROL_WZYX, lxrw_lyrw_lzrw_lwrw);

            let result1 = m128_mul_add(lyrz_lxrz_lwrz_lzrz, CONTROL_YXWZ, lzry_lwry_nlxry_nlyry);

            Self(_mm_add_ps(result0, result1))
        }
    }

    /// Creates a quaternion from a 3x3 rotation matrix inside a 3D affine transform.
    ///
    /// Note if the input affine matrix contain scales, shears, or other non-rotation
    /// transformations then the resulting quaternion will be ill-defined.
    ///
    /// # Panics
    ///
    /// Will panic if any input affine matrix column is not normalized when `glam_assert` is
    /// enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn from_affine3(a: &crate::Affine3) -> Self {
        Self::from_rotation_axes(a.matrix3.x_axis, a.matrix3.y_axis, a.matrix3.z_axis)
    }

    /// Creates a quaternion from a 3x3 rotation matrix inside a 3D affine transform.
    ///
    /// Note if the input affine matrix contain scales, shears, or other non-rotation
    /// transformations then the resulting quaternion will be ill-defined.
    ///
    /// # Panics
    ///
    /// Will panic if any input affine matrix column is not normalized when `glam_assert` is
    /// enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn from_affine3a(a: &crate::Affine3A) -> Self {
        Self::from_rotation_axes(
            a.matrix3.x_axis.into(),
            a.matrix3.y_axis.into(),
            a.matrix3.z_axis.into(),
        )
    }

    /// Multiplies a quaternion and a 3D vector, returning the rotated vector.
    ///
    /// # Panics
    ///
    /// Will panic if `self` is not normalized when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn mul_vec3a(self, rhs: Vec3A) -> Vec3A {
        glam_assert!(self.is_normalized());

        unsafe {
            const TWO: __m128 = m128_from_f32x4([2.0; 4]);
            let w = _mm_shuffle_ps(self.0, self.0, 0b11_11_11_11);
            let b = self.0;
            let b2 = dot3_into_m128(b, b);
            Vec3A(_mm_add_ps(
                _mm_add_ps(
                    _mm_mul_ps(rhs.0, _mm_sub_ps(_mm_mul_ps(w, w), b2)),
                    _mm_mul_ps(b, _mm_mul_ps(dot3_into_m128(rhs.0, b), TWO)),
                ),
                _mm_mul_ps(Vec3A(b).cross(rhs).into(), _mm_mul_ps(w, TWO)),
            ))
        }
    }

    #[cfg(feature = "f64")]
    #[inline]
    #[must_use]
    pub fn as_dquat(self) -> DQuat {
        DQuat::from_xyzw(self.x as f64, self.y as f64, self.z as f64, self.w as f64)
    }
}

impl fmt::Debug for Quat {
    fn fmt(&self, fmt: &mut fmt::Formatter<'_>) -> fmt::Result {
        fmt.debug_tuple(stringify!(Quat))
            .field(&self.x)
            .field(&self.y)
            .field(&self.z)
            .field(&self.w)
            .finish()
    }
}

impl fmt::Display for Quat {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> 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 Add for Quat {
    type Output = Self;
    /// Adds two quaternions.
    ///
    /// The sum is not guaranteed to be normalized.
    ///
    /// Note that addition is not the same as combining the rotations represented by the
    /// two quaternions! That corresponds to multiplication.
    #[inline]
    fn add(self, rhs: Self) -> Self {
        Self::from_vec4(Vec4::from(self) + Vec4::from(rhs))
    }
}

impl Add<&Self> for Quat {
    type Output = Self;
    #[inline]
    fn add(self, rhs: &Self) -> Self {
        self.add(*rhs)
    }
}

impl Add<&Quat> for &Quat {
    type Output = Quat;
    #[inline]
    fn add(self, rhs: &Quat) -> Quat {
        (*self).add(*rhs)
    }
}

impl Add<Quat> for &Quat {
    type Output = Quat;
    #[inline]
    fn add(self, rhs: Quat) -> Quat {
        (*self).add(rhs)
    }
}

impl AddAssign for Quat {
    #[inline]
    fn add_assign(&mut self, rhs: Self) {
        *self = self.add(rhs);
    }
}

impl AddAssign<&Self> for Quat {
    #[inline]
    fn add_assign(&mut self, rhs: &Self) {
        self.add_assign(*rhs);
    }
}

impl Sub for Quat {
    type Output = Self;
    /// Subtracts the `rhs` quaternion from `self`.
    ///
    /// The difference is not guaranteed to be normalized.
    #[inline]
    fn sub(self, rhs: Self) -> Self {
        Self::from_vec4(Vec4::from(self) - Vec4::from(rhs))
    }
}

impl Sub<&Self> for Quat {
    type Output = Self;
    #[inline]
    fn sub(self, rhs: &Self) -> Self {
        self.sub(*rhs)
    }
}

impl Sub<&Quat> for &Quat {
    type Output = Quat;
    #[inline]
    fn sub(self, rhs: &Quat) -> Quat {
        (*self).sub(*rhs)
    }
}

impl Sub<Quat> for &Quat {
    type Output = Quat;
    #[inline]
    fn sub(self, rhs: Quat) -> Quat {
        (*self).sub(rhs)
    }
}

impl SubAssign for Quat {
    #[inline]
    fn sub_assign(&mut self, rhs: Self) {
        *self = self.sub(rhs);
    }
}

impl SubAssign<&Self> for Quat {
    #[inline]
    fn sub_assign(&mut self, rhs: &Self) {
        self.sub_assign(*rhs);
    }
}

impl Mul<f32> for Quat {
    type Output = Self;
    /// Multiplies a quaternion by a scalar value.
    ///
    /// The product is not guaranteed to be normalized.
    #[inline]
    fn mul(self, rhs: f32) -> Self {
        Self::from_vec4(Vec4::from(self) * rhs)
    }
}

impl Mul<&f32> for Quat {
    type Output = Self;
    #[inline]
    fn mul(self, rhs: &f32) -> Self {
        self.mul(*rhs)
    }
}

impl Mul<&f32> for &Quat {
    type Output = Quat;
    #[inline]
    fn mul(self, rhs: &f32) -> Quat {
        (*self).mul(*rhs)
    }
}

impl Mul<f32> for &Quat {
    type Output = Quat;
    #[inline]
    fn mul(self, rhs: f32) -> Quat {
        (*self).mul(rhs)
    }
}

impl MulAssign<f32> for Quat {
    #[inline]
    fn mul_assign(&mut self, rhs: f32) {
        *self = self.mul(rhs);
    }
}

impl MulAssign<&f32> for Quat {
    #[inline]
    fn mul_assign(&mut self, rhs: &f32) {
        self.mul_assign(*rhs);
    }
}

impl Div<f32> for Quat {
    type Output = Self;
    /// Divides a quaternion by a scalar value.
    /// The quotient is not guaranteed to be normalized.
    #[inline]
    fn div(self, rhs: f32) -> Self {
        Self::from_vec4(Vec4::from(self) / rhs)
    }
}

impl Div<&f32> for Quat {
    type Output = Self;
    #[inline]
    fn div(self, rhs: &f32) -> Self {
        self.div(*rhs)
    }
}

impl Div<&f32> for &Quat {
    type Output = Quat;
    #[inline]
    fn div(self, rhs: &f32) -> Quat {
        (*self).div(*rhs)
    }
}

impl Div<f32> for &Quat {
    type Output = Quat;
    #[inline]
    fn div(self, rhs: f32) -> Quat {
        (*self).div(rhs)
    }
}

impl DivAssign<f32> for Quat {
    #[inline]
    fn div_assign(&mut self, rhs: f32) {
        *self = self.div(rhs);
    }
}

impl DivAssign<&f32> for Quat {
    #[inline]
    fn div_assign(&mut self, rhs: &f32) {
        self.div_assign(*rhs);
    }
}

impl Mul for Quat {
    type Output = Self;
    /// Multiplies two quaternions. If they each represent a rotation, the result will
    /// represent the combined rotation.
    ///
    /// Note that due to floating point rounding the result may not be perfectly
    /// normalized.
    #[inline]
    fn mul(self, rhs: Self) -> Self {
        self.mul_quat(rhs)
    }
}

impl Mul<&Self> for Quat {
    type Output = Self;
    #[inline]
    fn mul(self, rhs: &Self) -> Self {
        self.mul(*rhs)
    }
}

impl Mul<&Quat> for &Quat {
    type Output = Quat;
    #[inline]
    fn mul(self, rhs: &Quat) -> Quat {
        (*self).mul(*rhs)
    }
}

impl Mul<Quat> for &Quat {
    type Output = Quat;
    #[inline]
    fn mul(self, rhs: Quat) -> Quat {
        (*self).mul(rhs)
    }
}

impl MulAssign for Quat {
    #[inline]
    fn mul_assign(&mut self, rhs: Self) {
        *self = self.mul(rhs);
    }
}

impl MulAssign<&Self> for Quat {
    #[inline]
    fn mul_assign(&mut self, rhs: &Self) {
        self.mul_assign(*rhs);
    }
}

impl Mul<Vec3> for Quat {
    type Output = Vec3;
    /// Multiplies a quaternion and a 3D vector, returning the rotated vector.
    ///
    /// # Panics
    ///
    /// Will panic if `self` is not normalized when `glam_assert` is enabled.
    #[inline]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    fn mul(self, rhs: Vec3) -> Self::Output {
        self.mul_vec3(rhs)
    }
}

impl Mul<&Vec3> for Quat {
    type Output = Vec3;
    #[inline]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    fn mul(self, rhs: &Vec3) -> Vec3 {
        self.mul(*rhs)
    }
}

impl Mul<&Vec3> for &Quat {
    type Output = Vec3;
    #[inline]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    fn mul(self, rhs: &Vec3) -> Vec3 {
        (*self).mul(*rhs)
    }
}

impl Mul<Vec3> for &Quat {
    type Output = Vec3;
    #[inline]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    fn mul(self, rhs: Vec3) -> Vec3 {
        (*self).mul(rhs)
    }
}

impl Mul<Vec3A> for Quat {
    type Output = Vec3A;
    /// Multiplies a quaternion and a 3D vector, returning the rotated vector.
    ///
    /// # Panics
    ///
    /// Will panic if `self` is not normalized when `glam_assert` is enabled.
    #[inline]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    fn mul(self, rhs: Vec3A) -> Self::Output {
        self.mul_vec3a(rhs)
    }
}

impl Mul<&Vec3A> for Quat {
    type Output = Vec3A;
    #[inline]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    fn mul(self, rhs: &Vec3A) -> Vec3A {
        self.mul(*rhs)
    }
}

impl Mul<&Vec3A> for &Quat {
    type Output = Vec3A;
    #[inline]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    fn mul(self, rhs: &Vec3A) -> Vec3A {
        (*self).mul(*rhs)
    }
}

impl Mul<Vec3A> for &Quat {
    type Output = Vec3A;
    #[inline]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    fn mul(self, rhs: Vec3A) -> Vec3A {
        (*self).mul(rhs)
    }
}

impl Neg for Quat {
    type Output = Self;
    #[inline]
    fn neg(self) -> Self {
        self * -1.0
    }
}

impl Neg for &Quat {
    type Output = Quat;
    #[inline]
    fn neg(self) -> Quat {
        (*self).neg()
    }
}

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

impl PartialEq for Quat {
    #[inline]
    fn eq(&self, rhs: &Self) -> bool {
        Vec4::from(*self).eq(&Vec4::from(*rhs))
    }
}

impl AsRef<[f32; 4]> for Quat {
    #[inline]
    fn as_ref(&self) -> &[f32; 4] {
        unsafe { &*(self as *const Self as *const [f32; 4]) }
    }
}

impl Sum<Self> for Quat {
    fn sum<I>(iter: I) -> Self
    where
        I: Iterator<Item = Self>,
    {
        iter.fold(Self::ZERO, Self::add)
    }
}

impl<'a> Sum<&'a Self> for Quat {
    fn sum<I>(iter: I) -> Self
    where
        I: Iterator<Item = &'a Self>,
    {
        iter.fold(Self::ZERO, |a, &b| Self::add(a, b))
    }
}

impl Product for Quat {
    fn product<I>(iter: I) -> Self
    where
        I: Iterator<Item = Self>,
    {
        iter.fold(Self::IDENTITY, Self::mul)
    }
}

impl<'a> Product<&'a Self> for Quat {
    fn product<I>(iter: I) -> Self
    where
        I: Iterator<Item = &'a Self>,
    {
        iter.fold(Self::IDENTITY, |a, &b| Self::mul(a, b))
    }
}

impl From<Quat> for Vec4 {
    #[inline]
    fn from(q: Quat) -> Self {
        Self(q.0)
    }
}

impl From<Quat> for (f32, f32, f32, f32) {
    #[inline]
    fn from(q: Quat) -> Self {
        Vec4::from(q).into()
    }
}

impl From<Quat> for [f32; 4] {
    #[inline]
    fn from(q: Quat) -> Self {
        Vec4::from(q).into()
    }
}

impl From<Quat> for __m128 {
    #[inline]
    fn from(q: Quat) -> Self {
        q.0
    }
}

impl Deref for Quat {
    type Target = crate::deref::Vec4<f32>;
    #[inline]
    fn deref(&self) -> &Self::Target {
        unsafe { &*(self as *const Self).cast() }
    }
}

impl DerefMut for Quat {
    #[inline]
    fn deref_mut(&mut self) -> &mut Self::Target {
        unsafe { &mut *(self as *mut Self).cast() }
    }
}