glam/f32/sse2/

mat4.rs

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

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

use crate::{
    EulerRot, Mat3, Mat3A, Quat, Vec3, Vec3A, Vec4,
    euler::{FromEuler, ToEuler},
    f32::math,
    sse2::*,
    swizzles::*,
};
use core::fmt;
use core::iter::{Product, Sum};
use core::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign};

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

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

/// Creates a 4x4 matrix from four column vectors.
#[inline(always)]
#[must_use]
pub const fn mat4(x_axis: Vec4, y_axis: Vec4, z_axis: Vec4, w_axis: Vec4) -> Mat4 {
    Mat4::from_cols(x_axis, y_axis, z_axis, w_axis)
}

/// A 4x4 column major matrix.
///
/// If you are primarily dealing with 3D affine transformations
/// considering using [`Affine3A`](crate::Affine3A) which is faster than a 4x4 matrix
/// for some affine operations.
///
/// Affine transformations including 3D translation, rotation and scale can be created
/// using methods such as [`Self::from_translation()`], [`Self::from_quat()`],
/// [`Self::from_scale()`] and [`Self::from_scale_rotation_translation()`].
///
/// The [`Self::transform_point3()`] and [`Self::transform_vector3()`] convenience methods
/// are provided for performing affine transformations on 3D vectors and points. These
/// multiply 3D inputs as 4D vectors with an implicit `w` value of `1` for points and `0`
/// for vectors respectively. These methods assume that `Self` contains a valid affine
/// transform.
///
/// 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(C)]
pub struct Mat4 {
    pub x_axis: Vec4,
    pub y_axis: Vec4,
    pub z_axis: Vec4,
    pub w_axis: Vec4,
}

impl Mat4 {
    /// A 4x4 matrix with all elements set to `0.0`.
    pub const ZERO: Self = Self::from_cols(Vec4::ZERO, Vec4::ZERO, Vec4::ZERO, Vec4::ZERO);

    /// A 4x4 identity matrix, where all diagonal elements are `1`, and all off-diagonal elements are `0`.
    pub const IDENTITY: Self = Self::from_cols(Vec4::X, Vec4::Y, Vec4::Z, Vec4::W);

    /// All NAN:s.
    pub const NAN: Self = Self::from_cols(Vec4::NAN, Vec4::NAN, Vec4::NAN, Vec4::NAN);

    #[allow(clippy::too_many_arguments)]
    #[inline(always)]
    #[must_use]
    const fn new(
        m00: f32,
        m01: f32,
        m02: f32,
        m03: f32,
        m10: f32,
        m11: f32,
        m12: f32,
        m13: f32,
        m20: f32,
        m21: f32,
        m22: f32,
        m23: f32,
        m30: f32,
        m31: f32,
        m32: f32,
        m33: f32,
    ) -> Self {
        Self {
            x_axis: Vec4::new(m00, m01, m02, m03),
            y_axis: Vec4::new(m10, m11, m12, m13),
            z_axis: Vec4::new(m20, m21, m22, m23),
            w_axis: Vec4::new(m30, m31, m32, m33),
        }
    }

    /// Creates a 4x4 matrix from four column vectors.
    ///
    /// See also [`Self::from_rows`] when the data is in row major order.
    #[inline(always)]
    #[must_use]
    pub const fn from_cols(x_axis: Vec4, y_axis: Vec4, z_axis: Vec4, w_axis: Vec4) -> Self {
        Self {
            x_axis,
            y_axis,
            z_axis,
            w_axis,
        }
    }

    /// Creates a 4x4 matrix from four row vectors.
    ///
    /// Matrices are stored in column major order, so the given rows are permuted into
    /// the matrix layout. Use [`Self::from_cols`] instead when the data is already in
    /// column major order.
    #[inline(always)]
    #[must_use]
    pub const fn from_rows(row0: Vec4, row1: Vec4, row2: Vec4, row3: Vec4) -> Self {
        let [m00, m01, m02, m03] = row0.to_array();
        let [m10, m11, m12, m13] = row1.to_array();
        let [m20, m21, m22, m23] = row2.to_array();
        let [m30, m31, m32, m33] = row3.to_array();
        Self::new(
            m00, m10, m20, m30, m01, m11, m21, m31, m02, m12, m22, m32, m03, m13, m23, m33,
        )
    }

    /// Creates a 4x4 matrix from a `[f32; 16]` array stored in column major order.
    ///
    /// If the data is in row major order use [`Self::from_rows_array`] instead.
    #[inline]
    #[must_use]
    pub const fn from_cols_array(m: &[f32; 16]) -> Self {
        Self::new(
            m[0], m[1], m[2], m[3], m[4], m[5], m[6], m[7], m[8], m[9], m[10], m[11], m[12], m[13],
            m[14], m[15],
        )
    }

    /// Creates a `[f32; 16]` array storing data in column major order.
    ///
    /// If you require the data in row major order use [`Self::to_rows_array`] instead.
    #[inline]
    #[must_use]
    pub const fn to_cols_array(&self) -> [f32; 16] {
        let [x_axis_x, x_axis_y, x_axis_z, x_axis_w] = self.x_axis.to_array();
        let [y_axis_x, y_axis_y, y_axis_z, y_axis_w] = self.y_axis.to_array();
        let [z_axis_x, z_axis_y, z_axis_z, z_axis_w] = self.z_axis.to_array();
        let [w_axis_x, w_axis_y, w_axis_z, w_axis_w] = self.w_axis.to_array();

        [
            x_axis_x, x_axis_y, x_axis_z, x_axis_w, y_axis_x, y_axis_y, y_axis_z, y_axis_w,
            z_axis_x, z_axis_y, z_axis_z, z_axis_w, w_axis_x, w_axis_y, w_axis_z, w_axis_w,
        ]
    }

    /// Creates a 4x4 matrix from a `[[f32; 4]; 4]` 4D array stored in column major order.
    ///
    /// If the data is in row major order `transpose` the returned matrix.
    #[inline]
    #[must_use]
    pub const fn from_cols_array_2d(m: &[[f32; 4]; 4]) -> Self {
        Self::from_cols(
            Vec4::from_array(m[0]),
            Vec4::from_array(m[1]),
            Vec4::from_array(m[2]),
            Vec4::from_array(m[3]),
        )
    }

    /// Creates a `[[f32; 4]; 4]` 4D array storing data in column major order.
    ///
    /// If you require row major order `transpose` the matrix first.
    #[inline]
    #[must_use]
    pub const fn to_cols_array_2d(&self) -> [[f32; 4]; 4] {
        [
            self.x_axis.to_array(),
            self.y_axis.to_array(),
            self.z_axis.to_array(),
            self.w_axis.to_array(),
        ]
    }

    /// Creates a 4x4 matrix from a `[f32; 16]` array stored in row major order.
    ///
    /// Matrices are stored in column major order, so the array is permuted into the
    /// matrix layout. Use [`Self::from_cols_array`] instead when the data is already in
    /// column major order.
    #[inline]
    #[must_use]
    pub const fn from_rows_array(m: &[f32; 16]) -> Self {
        Self::new(
            m[0], m[4], m[8], m[12], m[1], m[5], m[9], m[13], m[2], m[6], m[10], m[14], m[3], m[7],
            m[11], m[15],
        )
    }

    /// Creates a `[f32; 16]` array storing data in row major order.
    ///
    /// Matrices are stored in column major order, so the array is permuted out of the
    /// column major storage. Use [`Self::to_cols_array`] instead when you want data in
    /// column major order.
    #[inline]
    #[must_use]
    pub const fn to_rows_array(&self) -> [f32; 16] {
        let m = self.to_cols_array();
        [
            m[0], m[4], m[8], m[12], m[1], m[5], m[9], m[13], m[2], m[6], m[10], m[14], m[3], m[7],
            m[11], m[15],
        ]
    }

    /// Creates a 4x4 matrix with its diagonal set to `diagonal` and all other entries set to 0.
    #[doc(alias = "scale")]
    #[inline]
    #[must_use]
    pub const fn from_diagonal(diagonal: Vec4) -> Self {
        // diagonal.x, diagonal.y etc can't be done in a const-context
        let [x, y, z, w] = diagonal.to_array();
        Self::new(
            x, 0.0, 0.0, 0.0, 0.0, y, 0.0, 0.0, 0.0, 0.0, z, 0.0, 0.0, 0.0, 0.0, w,
        )
    }

    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    fn quat_to_axes(rotation: Quat) -> (Vec4, Vec4, Vec4) {
        glam_assert!(rotation.is_normalized());

        let (x, y, z, w) = rotation.into();
        let x2 = x + x;
        let y2 = y + y;
        let z2 = z + z;
        let xx = x * x2;
        let xy = x * y2;
        let xz = x * z2;
        let yy = y * y2;
        let yz = y * z2;
        let zz = z * z2;
        let wx = w * x2;
        let wy = w * y2;
        let wz = w * z2;

        let x_axis = Vec4::new(1.0 - (yy + zz), xy + wz, xz - wy, 0.0);
        let y_axis = Vec4::new(xy - wz, 1.0 - (xx + zz), yz + wx, 0.0);
        let z_axis = Vec4::new(xz + wy, yz - wx, 1.0 - (xx + yy), 0.0);
        (x_axis, y_axis, z_axis)
    }

    /// Creates an affine transformation matrix from the given 3D `scale`, `rotation` and
    /// `translation`.
    ///
    /// The resulting matrix can be used to transform 3D points and vectors. See
    /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
    ///
    /// # Panics
    ///
    /// Will panic if `rotation` is not normalized when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn from_scale_rotation_translation(scale: Vec3, rotation: Quat, translation: Vec3) -> Self {
        let (x_axis, y_axis, z_axis) = Self::quat_to_axes(rotation);
        Self::from_cols(
            x_axis.mul(scale.x),
            y_axis.mul(scale.y),
            z_axis.mul(scale.z),
            Vec4::from((translation, 1.0)),
        )
    }

    /// Creates an affine transformation matrix from the given 3D `translation`.
    ///
    /// The resulting matrix can be used to transform 3D points and vectors. See
    /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
    ///
    /// # Panics
    ///
    /// Will panic if `rotation` 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_translation(rotation: Quat, translation: Vec3) -> Self {
        let (x_axis, y_axis, z_axis) = Self::quat_to_axes(rotation);
        Self::from_cols(x_axis, y_axis, z_axis, Vec4::from((translation, 1.0)))
    }

    /// Extracts `scale`, `rotation` and `translation` from `self`. The input matrix is
    /// expected to be a 3D affine transformation matrix otherwise the output will be invalid.
    ///
    /// # Panics
    /// Will panic if `self` is not a valid affine transformation matrix, if the determinant of the
    /// 3x3 linear part (the rotation and scale part of the transform) is zero, when `glam_assert`
    /// is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn to_scale_rotation_translation(&self) -> (Vec3, Quat, Vec3) {
        glam_assert!(self.row(3).abs_diff_eq(Vec4::W, 1e-6));

        let r = Mat3A::from_mat4(*self);

        let det = r.determinant();

        glam_assert!(det != 0.0);

        let scale = Vec3::new(
            r.x_axis.length() * det.signum(),
            r.y_axis.length(),
            r.z_axis.length(),
        );

        glam_assert!(scale.cmpne(Vec3::ZERO).all());

        let inv_scale = scale.recip();

        let rotation = Quat::from_mat3a(&Mat3A::from_cols(
            r.x_axis.mul(inv_scale.x),
            r.y_axis.mul(inv_scale.y),
            r.z_axis.mul(inv_scale.z),
        ));

        let translation = self.w_axis.xyz();

        (scale, rotation, translation)
    }

    /// Creates an affine transformation matrix from the given `rotation` quaternion.
    ///
    /// The resulting matrix can be used to transform 3D points and vectors. See
    /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
    ///
    /// # Panics
    ///
    /// Will panic if `rotation` is not normalized when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn from_quat(rotation: Quat) -> Self {
        let (x_axis, y_axis, z_axis) = Self::quat_to_axes(rotation);
        Self::from_cols(x_axis, y_axis, z_axis, Vec4::W)
    }

    /// Creates an affine transformation matrix from the given 3x3 linear transformation
    /// matrix.
    ///
    /// The resulting matrix can be used to transform 3D points and vectors. See
    /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
    #[inline]
    #[must_use]
    pub fn from_mat3(m: Mat3) -> Self {
        Self::from_cols(
            Vec4::from((m.x_axis, 0.0)),
            Vec4::from((m.y_axis, 0.0)),
            Vec4::from((m.z_axis, 0.0)),
            Vec4::W,
        )
    }

    /// Creates an affine transformation matrics from a 3x3 matrix (expressing scale, shear and
    /// rotation) and a translation vector.
    ///
    /// Equivalent to `Mat4::from_translation(translation) * Mat4::from_mat3(mat3)`
    #[inline]
    #[must_use]
    pub fn from_mat3_translation(mat3: Mat3, translation: Vec3) -> Self {
        Self::from_cols(
            Vec4::from((mat3.x_axis, 0.0)),
            Vec4::from((mat3.y_axis, 0.0)),
            Vec4::from((mat3.z_axis, 0.0)),
            Vec4::from((translation, 1.0)),
        )
    }

    /// Creates an affine transformation matrix from the given 3x3 linear transformation
    /// matrix.
    ///
    /// The resulting matrix can be used to transform 3D points and vectors. See
    /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
    #[inline]
    #[must_use]
    pub fn from_mat3a(m: Mat3A) -> Self {
        Self::from_cols(
            Vec4::from((m.x_axis, 0.0)),
            Vec4::from((m.y_axis, 0.0)),
            Vec4::from((m.z_axis, 0.0)),
            Vec4::W,
        )
    }

    /// Creates an affine transformation matrix from the given 3D `translation`.
    ///
    /// The resulting matrix can be used to transform 3D points and vectors. See
    /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
    #[inline]
    #[must_use]
    pub fn from_translation(translation: Vec3) -> Self {
        Self::from_cols(
            Vec4::X,
            Vec4::Y,
            Vec4::Z,
            Vec4::new(translation.x, translation.y, translation.z, 1.0),
        )
    }

    /// Creates an affine transformation matrix containing a 3D rotation around a normalized
    /// rotation `axis` of `angle` (in radians).
    ///
    /// The resulting matrix can be used to transform 3D points and vectors. See
    /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
    ///
    /// # 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 (sin, cos) = math::sin_cos(angle);
        let axis_sin = axis.mul(sin);
        let axis_sq = axis.mul(axis);
        let omc = 1.0 - cos;
        let xyomc = axis.x * axis.y * omc;
        let xzomc = axis.x * axis.z * omc;
        let yzomc = axis.y * axis.z * omc;
        Self::from_cols(
            Vec4::new(
                axis_sq.x * omc + cos,
                xyomc + axis_sin.z,
                xzomc - axis_sin.y,
                0.0,
            ),
            Vec4::new(
                xyomc - axis_sin.z,
                axis_sq.y * omc + cos,
                yzomc + axis_sin.x,
                0.0,
            ),
            Vec4::new(
                xzomc + axis_sin.y,
                yzomc - axis_sin.x,
                axis_sq.z * omc + cos,
                0.0,
            ),
            Vec4::W,
        )
    }

    /// Creates a affine transformation matrix containing a rotation from the given euler
    /// rotation sequence and angles (in radians).
    ///
    /// The resulting matrix can be used to transform 3D points and vectors. See
    /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
    #[inline]
    #[must_use]
    pub fn from_euler(order: EulerRot, a: f32, b: f32, c: f32) -> Self {
        Self::from_euler_angles(order, a, b, c)
    }

    /// Extract Euler angles with the given Euler rotation order.
    ///
    /// Note if the upper 3x3 matrix contain scales, shears, or other non-rotation transformations
    /// then the resulting Euler angles 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 to_euler(&self, order: EulerRot) -> (f32, f32, f32) {
        glam_assert!(
            self.x_axis.xyz().is_normalized()
                && self.y_axis.xyz().is_normalized()
                && self.z_axis.xyz().is_normalized()
        );
        Mat3::from_mat4(*self).to_euler_angles(order)
    }

    /// Creates an affine transformation matrix containing a 3D rotation around the x axis of
    /// `angle` (in radians).
    ///
    /// The resulting matrix can be used to transform 3D points and vectors. See
    /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
    #[inline]
    #[must_use]
    pub fn from_rotation_x(angle: f32) -> Self {
        let (sina, cosa) = math::sin_cos(angle);
        Self::from_cols(
            Vec4::X,
            Vec4::new(0.0, cosa, sina, 0.0),
            Vec4::new(0.0, -sina, cosa, 0.0),
            Vec4::W,
        )
    }

    /// Creates an affine transformation matrix containing a 3D rotation around the y axis of
    /// `angle` (in radians).
    ///
    /// The resulting matrix can be used to transform 3D points and vectors. See
    /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
    #[inline]
    #[must_use]
    pub fn from_rotation_y(angle: f32) -> Self {
        let (sina, cosa) = math::sin_cos(angle);
        Self::from_cols(
            Vec4::new(cosa, 0.0, -sina, 0.0),
            Vec4::Y,
            Vec4::new(sina, 0.0, cosa, 0.0),
            Vec4::W,
        )
    }

    /// Creates an affine transformation matrix containing a 3D rotation around the z axis of
    /// `angle` (in radians).
    ///
    /// The resulting matrix can be used to transform 3D points and vectors. See
    /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
    #[inline]
    #[must_use]
    pub fn from_rotation_z(angle: f32) -> Self {
        let (sina, cosa) = math::sin_cos(angle);
        Self::from_cols(
            Vec4::new(cosa, sina, 0.0, 0.0),
            Vec4::new(-sina, cosa, 0.0, 0.0),
            Vec4::Z,
            Vec4::W,
        )
    }

    /// Creates an affine transformation matrix containing the given 3D non-uniform `scale`.
    ///
    /// The resulting matrix can be used to transform 3D points and vectors. See
    /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
    ///
    /// # Panics
    ///
    /// Will panic if all elements of `scale` are zero when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn from_scale(scale: Vec3) -> Self {
        // Do not panic as long as any component is non-zero
        glam_assert!(scale.cmpne(Vec3::ZERO).any());

        Self::from_cols(
            Vec4::new(scale.x, 0.0, 0.0, 0.0),
            Vec4::new(0.0, scale.y, 0.0, 0.0),
            Vec4::new(0.0, 0.0, scale.z, 0.0),
            Vec4::W,
        )
    }

    /// Creates a 4x4 matrix from the first 16 values in `slice`.
    ///
    /// See also [`Self::from_rows_slice`] when the slice is in row major order.
    ///
    /// # Panics
    ///
    /// Panics if `slice` is less than 16 elements long.
    #[inline]
    #[must_use]
    #[track_caller]
    pub const fn from_cols_slice(slice: &[f32]) -> Self {
        Self::new(
            slice[0], slice[1], slice[2], slice[3], slice[4], slice[5], slice[6], slice[7],
            slice[8], slice[9], slice[10], slice[11], slice[12], slice[13], slice[14], slice[15],
        )
    }

    /// Writes the columns of `self` to the first 16 elements in `slice`.
    ///
    /// # Panics
    ///
    /// Panics if `slice` is less than 16 elements long.
    #[inline]
    #[track_caller]
    pub fn write_cols_to_slice(&self, slice: &mut [f32]) {
        slice[0] = self.x_axis.x;
        slice[1] = self.x_axis.y;
        slice[2] = self.x_axis.z;
        slice[3] = self.x_axis.w;
        slice[4] = self.y_axis.x;
        slice[5] = self.y_axis.y;
        slice[6] = self.y_axis.z;
        slice[7] = self.y_axis.w;
        slice[8] = self.z_axis.x;
        slice[9] = self.z_axis.y;
        slice[10] = self.z_axis.z;
        slice[11] = self.z_axis.w;
        slice[12] = self.w_axis.x;
        slice[13] = self.w_axis.y;
        slice[14] = self.w_axis.z;
        slice[15] = self.w_axis.w;
    }

    /// Creates a 4x4 matrix from the first 16 values in `slice`, stored in row
    /// major order.
    ///
    /// Matrices are stored in column major order, so the slice is permuted into the
    /// matrix layout. Use [`Self::from_cols_slice`] instead when the slice is already in
    /// column major order.
    ///
    /// # Panics
    ///
    /// Panics if `slice` is less than 16 elements long.
    #[inline]
    #[must_use]
    #[track_caller]
    pub const fn from_rows_slice(slice: &[f32]) -> Self {
        Self::new(
            slice[0], slice[4], slice[8], slice[12], slice[1], slice[5], slice[9], slice[13],
            slice[2], slice[6], slice[10], slice[14], slice[3], slice[7], slice[11], slice[15],
        )
    }

    /// Returns the matrix column for the given `index`.
    ///
    /// # Panics
    ///
    /// Panics if `index` is greater than 3.
    #[inline]
    #[must_use]
    #[track_caller]
    pub fn col(&self, index: usize) -> Vec4 {
        match index {
            0 => self.x_axis,
            1 => self.y_axis,
            2 => self.z_axis,
            3 => self.w_axis,
            _ => panic!("index out of bounds"),
        }
    }

    /// Returns a mutable reference to the matrix column for the given `index`.
    ///
    /// # Panics
    ///
    /// Panics if `index` is greater than 3.
    #[inline]
    #[track_caller]
    pub fn col_mut(&mut self, index: usize) -> &mut Vec4 {
        match index {
            0 => &mut self.x_axis,
            1 => &mut self.y_axis,
            2 => &mut self.z_axis,
            3 => &mut self.w_axis,
            _ => panic!("index out of bounds"),
        }
    }

    /// Returns the matrix row for the given `index`.
    ///
    /// See also [`Self::set_row`] when you need to change the row.
    ///
    /// # Panics
    ///
    /// Panics if `index` is greater than 3.
    #[inline]
    #[must_use]
    #[track_caller]
    pub fn row(&self, index: usize) -> Vec4 {
        match index {
            0 => Vec4::new(self.x_axis.x, self.y_axis.x, self.z_axis.x, self.w_axis.x),
            1 => Vec4::new(self.x_axis.y, self.y_axis.y, self.z_axis.y, self.w_axis.y),
            2 => Vec4::new(self.x_axis.z, self.y_axis.z, self.z_axis.z, self.w_axis.z),
            3 => Vec4::new(self.x_axis.w, self.y_axis.w, self.z_axis.w, self.w_axis.w),
            _ => panic!("index out of bounds"),
        }
    }

    /// Sets the matrix row for the given `index`.
    ///
    /// Matrices are stored in column major order, so the row is spread across all
    /// 4 columns and writing it touches every column. Use [`Self::col_mut`]
    /// instead when you can work with columns. See also [`Self::row`].
    ///
    /// # Panics
    ///
    /// Panics if `index` is greater than 3.
    #[inline]
    #[track_caller]
    pub fn set_row(&mut self, index: usize, row: Vec4) {
        match index {
            0 => {
                self.x_axis.x = row.x;
                self.y_axis.x = row.y;
                self.z_axis.x = row.z;
                self.w_axis.x = row.w;
            }
            1 => {
                self.x_axis.y = row.x;
                self.y_axis.y = row.y;
                self.z_axis.y = row.z;
                self.w_axis.y = row.w;
            }
            2 => {
                self.x_axis.z = row.x;
                self.y_axis.z = row.y;
                self.z_axis.z = row.z;
                self.w_axis.z = row.w;
            }
            3 => {
                self.x_axis.w = row.x;
                self.y_axis.w = row.y;
                self.z_axis.w = row.z;
                self.w_axis.w = row.w;
            }
            _ => panic!("index out of bounds"),
        }
    }

    /// 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 {
        self.x_axis.is_finite()
            && self.y_axis.is_finite()
            && self.z_axis.is_finite()
            && self.w_axis.is_finite()
    }

    /// Returns `true` if any elements are `NaN`.
    #[inline]
    #[must_use]
    pub fn is_nan(&self) -> bool {
        self.x_axis.is_nan() || self.y_axis.is_nan() || self.z_axis.is_nan() || self.w_axis.is_nan()
    }

    /// Returns the transpose of `self`.
    #[inline]
    #[must_use]
    pub fn transpose(&self) -> Self {
        unsafe {
            // Based on https://github.com/microsoft/DirectXMath `XMMatrixTranspose`
            let tmp0 = _mm_shuffle_ps(self.x_axis.0, self.y_axis.0, 0b01_00_01_00);
            let tmp1 = _mm_shuffle_ps(self.x_axis.0, self.y_axis.0, 0b11_10_11_10);
            let tmp2 = _mm_shuffle_ps(self.z_axis.0, self.w_axis.0, 0b01_00_01_00);
            let tmp3 = _mm_shuffle_ps(self.z_axis.0, self.w_axis.0, 0b11_10_11_10);

            Self {
                x_axis: Vec4(_mm_shuffle_ps(tmp0, tmp2, 0b10_00_10_00)),
                y_axis: Vec4(_mm_shuffle_ps(tmp0, tmp2, 0b11_01_11_01)),
                z_axis: Vec4(_mm_shuffle_ps(tmp1, tmp3, 0b10_00_10_00)),
                w_axis: Vec4(_mm_shuffle_ps(tmp1, tmp3, 0b11_01_11_01)),
            }
        }
    }

    /// Returns the diagonal of `self`.
    #[inline]
    #[must_use]
    pub fn diagonal(&self) -> Vec4 {
        Vec4::new(self.x_axis.x, self.y_axis.y, self.z_axis.z, self.w_axis.w)
    }

    /// Returns the determinant of `self`.
    #[must_use]
    pub fn determinant(&self) -> f32 {
        unsafe {
            // Based on https://github.com/g-truc/glm `glm_mat4_determinant_lowp`
            let swp2a = _mm_shuffle_ps(self.z_axis.0, self.z_axis.0, 0b00_01_01_10);
            let swp3a = _mm_shuffle_ps(self.w_axis.0, self.w_axis.0, 0b11_10_11_11);
            let swp2b = _mm_shuffle_ps(self.z_axis.0, self.z_axis.0, 0b11_10_11_11);
            let swp3b = _mm_shuffle_ps(self.w_axis.0, self.w_axis.0, 0b00_01_01_10);
            let swp2c = _mm_shuffle_ps(self.z_axis.0, self.z_axis.0, 0b00_00_01_10);
            let swp3c = _mm_shuffle_ps(self.w_axis.0, self.w_axis.0, 0b01_10_00_00);

            let mula = _mm_mul_ps(swp2a, swp3a);
            let mulb = _mm_mul_ps(swp2b, swp3b);
            let mulc = _mm_mul_ps(swp2c, swp3c);
            let sube = _mm_sub_ps(mula, mulb);
            let subf = _mm_sub_ps(_mm_movehl_ps(mulc, mulc), mulc);

            let subfaca = _mm_shuffle_ps(sube, sube, 0b10_01_00_00);
            let swpfaca = _mm_shuffle_ps(self.y_axis.0, self.y_axis.0, 0b00_00_00_01);
            let mulfaca = _mm_mul_ps(swpfaca, subfaca);

            let subtmpb = _mm_shuffle_ps(sube, subf, 0b00_00_11_01);
            let subfacb = _mm_shuffle_ps(subtmpb, subtmpb, 0b11_01_01_00);
            let swpfacb = _mm_shuffle_ps(self.y_axis.0, self.y_axis.0, 0b01_01_10_10);
            let mulfacb = _mm_mul_ps(swpfacb, subfacb);

            let subres = _mm_sub_ps(mulfaca, mulfacb);
            let subtmpc = _mm_shuffle_ps(sube, subf, 0b01_00_10_10);
            let subfacc = _mm_shuffle_ps(subtmpc, subtmpc, 0b11_11_10_00);
            let swpfacc = _mm_shuffle_ps(self.y_axis.0, self.y_axis.0, 0b10_11_11_11);
            let mulfacc = _mm_mul_ps(swpfacc, subfacc);

            let addres = _mm_add_ps(subres, mulfacc);
            let detcof = _mm_mul_ps(addres, _mm_setr_ps(1.0, -1.0, 1.0, -1.0));

            dot4(self.x_axis.0, detcof)
        }
    }

    /// If `CHECKED` is true then if the determinant is zero this function will return a tuple
    /// containing a zero matrix and false. If the determinant is non zero a tuple containing the
    /// inverted matrix and true is returned.
    ///
    /// If `CHECKED` is false then the determinant is not checked and if it is zero the resulting
    /// inverted matrix will be invalid. Will panic if the resulting inverted matrix is not finite
    /// when `glam_assert` is enabled.
    ///
    /// A tuple containing the inverted matrix and a bool is used instead of an option here as
    /// regular Rust enums put the discriminant first which can result in a lot of padding if the
    /// matrix is aligned.
    #[inline(always)]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    fn inverse_checked<const CHECKED: bool>(&self) -> (Self, bool) {
        unsafe {
            // Based on https://github.com/g-truc/glm `glm_mat4_inverse`
            let fac0 = {
                let swp0a = _mm_shuffle_ps(self.w_axis.0, self.z_axis.0, 0b11_11_11_11);
                let swp0b = _mm_shuffle_ps(self.w_axis.0, self.z_axis.0, 0b10_10_10_10);

                let swp00 = _mm_shuffle_ps(self.z_axis.0, self.y_axis.0, 0b10_10_10_10);
                let swp01 = _mm_shuffle_ps(swp0a, swp0a, 0b10_00_00_00);
                let swp02 = _mm_shuffle_ps(swp0b, swp0b, 0b10_00_00_00);
                let swp03 = _mm_shuffle_ps(self.z_axis.0, self.y_axis.0, 0b11_11_11_11);

                let mul00 = _mm_mul_ps(swp00, swp01);
                let mul01 = _mm_mul_ps(swp02, swp03);
                _mm_sub_ps(mul00, mul01)
            };
            let fac1 = {
                let swp0a = _mm_shuffle_ps(self.w_axis.0, self.z_axis.0, 0b11_11_11_11);
                let swp0b = _mm_shuffle_ps(self.w_axis.0, self.z_axis.0, 0b01_01_01_01);

                let swp00 = _mm_shuffle_ps(self.z_axis.0, self.y_axis.0, 0b01_01_01_01);
                let swp01 = _mm_shuffle_ps(swp0a, swp0a, 0b10_00_00_00);
                let swp02 = _mm_shuffle_ps(swp0b, swp0b, 0b10_00_00_00);
                let swp03 = _mm_shuffle_ps(self.z_axis.0, self.y_axis.0, 0b11_11_11_11);

                let mul00 = _mm_mul_ps(swp00, swp01);
                let mul01 = _mm_mul_ps(swp02, swp03);
                _mm_sub_ps(mul00, mul01)
            };
            let fac2 = {
                let swp0a = _mm_shuffle_ps(self.w_axis.0, self.z_axis.0, 0b10_10_10_10);
                let swp0b = _mm_shuffle_ps(self.w_axis.0, self.z_axis.0, 0b01_01_01_01);

                let swp00 = _mm_shuffle_ps(self.z_axis.0, self.y_axis.0, 0b01_01_01_01);
                let swp01 = _mm_shuffle_ps(swp0a, swp0a, 0b10_00_00_00);
                let swp02 = _mm_shuffle_ps(swp0b, swp0b, 0b10_00_00_00);
                let swp03 = _mm_shuffle_ps(self.z_axis.0, self.y_axis.0, 0b10_10_10_10);

                let mul00 = _mm_mul_ps(swp00, swp01);
                let mul01 = _mm_mul_ps(swp02, swp03);
                _mm_sub_ps(mul00, mul01)
            };
            let fac3 = {
                let swp0a = _mm_shuffle_ps(self.w_axis.0, self.z_axis.0, 0b11_11_11_11);
                let swp0b = _mm_shuffle_ps(self.w_axis.0, self.z_axis.0, 0b00_00_00_00);

                let swp00 = _mm_shuffle_ps(self.z_axis.0, self.y_axis.0, 0b00_00_00_00);
                let swp01 = _mm_shuffle_ps(swp0a, swp0a, 0b10_00_00_00);
                let swp02 = _mm_shuffle_ps(swp0b, swp0b, 0b10_00_00_00);
                let swp03 = _mm_shuffle_ps(self.z_axis.0, self.y_axis.0, 0b11_11_11_11);

                let mul00 = _mm_mul_ps(swp00, swp01);
                let mul01 = _mm_mul_ps(swp02, swp03);
                _mm_sub_ps(mul00, mul01)
            };
            let fac4 = {
                let swp0a = _mm_shuffle_ps(self.w_axis.0, self.z_axis.0, 0b10_10_10_10);
                let swp0b = _mm_shuffle_ps(self.w_axis.0, self.z_axis.0, 0b00_00_00_00);

                let swp00 = _mm_shuffle_ps(self.z_axis.0, self.y_axis.0, 0b00_00_00_00);
                let swp01 = _mm_shuffle_ps(swp0a, swp0a, 0b10_00_00_00);
                let swp02 = _mm_shuffle_ps(swp0b, swp0b, 0b10_00_00_00);
                let swp03 = _mm_shuffle_ps(self.z_axis.0, self.y_axis.0, 0b10_10_10_10);

                let mul00 = _mm_mul_ps(swp00, swp01);
                let mul01 = _mm_mul_ps(swp02, swp03);
                _mm_sub_ps(mul00, mul01)
            };
            let fac5 = {
                let swp0a = _mm_shuffle_ps(self.w_axis.0, self.z_axis.0, 0b01_01_01_01);
                let swp0b = _mm_shuffle_ps(self.w_axis.0, self.z_axis.0, 0b00_00_00_00);

                let swp00 = _mm_shuffle_ps(self.z_axis.0, self.y_axis.0, 0b00_00_00_00);
                let swp01 = _mm_shuffle_ps(swp0a, swp0a, 0b10_00_00_00);
                let swp02 = _mm_shuffle_ps(swp0b, swp0b, 0b10_00_00_00);
                let swp03 = _mm_shuffle_ps(self.z_axis.0, self.y_axis.0, 0b01_01_01_01);

                let mul00 = _mm_mul_ps(swp00, swp01);
                let mul01 = _mm_mul_ps(swp02, swp03);
                _mm_sub_ps(mul00, mul01)
            };
            let sign_a = _mm_set_ps(1.0, -1.0, 1.0, -1.0);
            let sign_b = _mm_set_ps(-1.0, 1.0, -1.0, 1.0);

            let temp0 = _mm_shuffle_ps(self.y_axis.0, self.x_axis.0, 0b00_00_00_00);
            let vec0 = _mm_shuffle_ps(temp0, temp0, 0b10_10_10_00);

            let temp1 = _mm_shuffle_ps(self.y_axis.0, self.x_axis.0, 0b01_01_01_01);
            let vec1 = _mm_shuffle_ps(temp1, temp1, 0b10_10_10_00);

            let temp2 = _mm_shuffle_ps(self.y_axis.0, self.x_axis.0, 0b10_10_10_10);
            let vec2 = _mm_shuffle_ps(temp2, temp2, 0b10_10_10_00);

            let temp3 = _mm_shuffle_ps(self.y_axis.0, self.x_axis.0, 0b11_11_11_11);
            let vec3 = _mm_shuffle_ps(temp3, temp3, 0b10_10_10_00);

            let mul00 = _mm_mul_ps(vec1, fac0);
            let mul01 = _mm_mul_ps(vec2, fac1);
            let mul02 = _mm_mul_ps(vec3, fac2);
            let sub00 = _mm_sub_ps(mul00, mul01);
            let add00 = _mm_add_ps(sub00, mul02);
            let inv0 = _mm_mul_ps(sign_b, add00);

            let mul03 = _mm_mul_ps(vec0, fac0);
            let mul04 = _mm_mul_ps(vec2, fac3);
            let mul05 = _mm_mul_ps(vec3, fac4);
            let sub01 = _mm_sub_ps(mul03, mul04);
            let add01 = _mm_add_ps(sub01, mul05);
            let inv1 = _mm_mul_ps(sign_a, add01);

            let mul06 = _mm_mul_ps(vec0, fac1);
            let mul07 = _mm_mul_ps(vec1, fac3);
            let mul08 = _mm_mul_ps(vec3, fac5);
            let sub02 = _mm_sub_ps(mul06, mul07);
            let add02 = _mm_add_ps(sub02, mul08);
            let inv2 = _mm_mul_ps(sign_b, add02);

            let mul09 = _mm_mul_ps(vec0, fac2);
            let mul10 = _mm_mul_ps(vec1, fac4);
            let mul11 = _mm_mul_ps(vec2, fac5);
            let sub03 = _mm_sub_ps(mul09, mul10);
            let add03 = _mm_add_ps(sub03, mul11);
            let inv3 = _mm_mul_ps(sign_a, add03);

            let row0 = _mm_shuffle_ps(inv0, inv1, 0b00_00_00_00);
            let row1 = _mm_shuffle_ps(inv2, inv3, 0b00_00_00_00);
            let row2 = _mm_shuffle_ps(row0, row1, 0b10_00_10_00);

            let dot0 = dot4(self.x_axis.0, row2);

            let rcp0 = _mm_set1_ps(dot0.recip());
            let m = Self {
                x_axis: Vec4(_mm_mul_ps(inv0, rcp0)),
                y_axis: Vec4(_mm_mul_ps(inv1, rcp0)),
                z_axis: Vec4(_mm_mul_ps(inv2, rcp0)),
                w_axis: Vec4(_mm_mul_ps(inv3, rcp0)),
            };
            if CHECKED {
                if !m.is_finite() {
                    return (Self::ZERO, false);
                }
            } else {
                glam_assert!(m.is_finite());
            }
            (m, true)
        }
    }

    /// Returns the inverse of `self`.
    ///
    /// If the matrix is not invertible the returned matrix will be invalid. The
    /// returned matrix will also be invalid if the inverse is not finite, which can
    /// happen when `self` contains very large or very small values. Use
    /// [`Self::try_inverse`] or [`Self::inverse_or_zero`] to detect these cases.
    ///
    /// # Panics
    ///
    /// Will panic if the resulting inverted matrix is not finite when `glam_assert`
    /// is enabled.
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn inverse(&self) -> Self {
        self.inverse_checked::<false>().0
    }

    /// Returns the inverse of `self` or `None` if the matrix is not invertible, or if
    /// the inverse is not finite.
    #[must_use]
    pub fn try_inverse(&self) -> Option<Self> {
        let (m, is_valid) = self.inverse_checked::<true>();
        if is_valid { Some(m) } else { None }
    }

    /// Returns the inverse of `self` or `Mat4::ZERO` if the matrix is not
    /// invertible, or if the inverse is not finite.
    #[must_use]
    pub fn inverse_or_zero(&self) -> Self {
        self.inverse_checked::<true>().0
    }

    /// Transforms the given 3D vector as a point, applying perspective correction.
    ///
    /// This is the equivalent of multiplying the 3D vector as a 4D vector where `w` is `1.0`.
    /// The perspective divide is performed meaning the resulting 3D vector is divided by `w`.
    ///
    /// This method assumes that `self` contains a projective transform.
    #[inline]
    #[must_use]
    pub fn project_point3(&self, rhs: Vec3) -> Vec3 {
        let mut res = self.x_axis.mul(rhs.x);
        res = self.y_axis.mul(rhs.y).add(res);
        res = self.z_axis.mul(rhs.z).add(res);
        res = self.w_axis.add(res);
        res = res.div(res.w);
        res.xyz()
    }

    /// Transforms the given 3D vector as a point.
    ///
    /// This is the equivalent of multiplying the 3D vector as a 4D vector where `w` is
    /// `1.0`.
    ///
    /// This method assumes that `self` contains a valid affine transform. It does not perform
    /// a perspective divide, if `self` contains a perspective transform, or if you are unsure,
    /// the [`Self::project_point3()`] method should be used instead.
    ///
    /// # Panics
    ///
    /// Will panic if the 3rd row of `self` is not `(0, 0, 0, 1)` when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn transform_point3(&self, rhs: Vec3) -> Vec3 {
        glam_assert!(self.row(3).abs_diff_eq(Vec4::W, 1e-6));
        let mut res = self.x_axis.mul(rhs.x);
        res = self.y_axis.mul(rhs.y).add(res);
        res = self.z_axis.mul(rhs.z).add(res);
        res = self.w_axis.add(res);
        res.xyz()
    }

    /// Transforms the given 3D vector as a direction.
    ///
    /// This is the equivalent of multiplying the 3D vector as a 4D vector where `w` is
    /// `0.0`.
    ///
    /// This method assumes that `self` contains a valid affine transform.
    ///
    /// # Panics
    ///
    /// Will panic if the 3rd row of `self` is not `(0, 0, 0, 1)` when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn transform_vector3(&self, rhs: Vec3) -> Vec3 {
        glam_assert!(self.row(3).abs_diff_eq(Vec4::W, 1e-6));
        let mut res = self.x_axis.mul(rhs.x);
        res = self.y_axis.mul(rhs.y).add(res);
        res = self.z_axis.mul(rhs.z).add(res);
        res.xyz()
    }

    /// Transforms the given [`Vec3A`] as a 3D point, applying perspective correction.
    ///
    /// This is the equivalent of multiplying the [`Vec3A`] as a 4D vector where `w` is `1.0`.
    /// The perspective divide is performed meaning the resulting 3D vector is divided by `w`.
    ///
    /// This method assumes that `self` contains a projective transform.
    #[inline]
    #[must_use]
    pub fn project_point3a(&self, rhs: Vec3A) -> Vec3A {
        let mut res = self.x_axis.mul(rhs.xxxx());
        res = self.y_axis.mul(rhs.yyyy()).add(res);
        res = self.z_axis.mul(rhs.zzzz()).add(res);
        res = self.w_axis.add(res);
        res = res.div(res.wwww());
        Vec3A::from_vec4(res)
    }

    /// Transforms the given [`Vec3A`] as 3D point.
    ///
    /// This is the equivalent of multiplying the [`Vec3A`] as a 4D vector where `w` is `1.0`.
    ///
    /// This method assumes that `self` contains a valid affine transform. If `self` contains a
    /// perspective transform, or if you are unsure, the [`Self::project_point3a()`] method should
    /// be used instead.
    ///
    /// # Panics
    ///
    /// Will panic if the 3rd row of `self` is not `(0, 0, 0, 1)` when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn transform_point3a(&self, rhs: Vec3A) -> Vec3A {
        glam_assert!(self.row(3).abs_diff_eq(Vec4::W, 1e-6));
        let mut res = self.x_axis.mul(rhs.xxxx());
        res = self.y_axis.mul(rhs.yyyy()).add(res);
        res = self.z_axis.mul(rhs.zzzz()).add(res);
        res = self.w_axis.add(res);
        Vec3A::from_vec4(res)
    }

    /// Transforms the given [`Vec3A`] as 3D vector.
    ///
    /// This is the equivalent of multiplying the [`Vec3A`] as a 4D vector where `w` is `0.0`.
    ///
    /// This method assumes that `self` contains a valid affine transform.
    ///
    /// # Panics
    ///
    /// Will panic if the 3rd row of `self` is not `(0, 0, 0, 1)` when `glam_assert` is enabled.
    #[inline]
    #[must_use]
    #[cfg_attr(any(debug_assertions, feature = "glam-assert"), track_caller)]
    pub fn transform_vector3a(&self, rhs: Vec3A) -> Vec3A {
        glam_assert!(self.row(3).abs_diff_eq(Vec4::W, 1e-6));
        let mut res = self.x_axis.mul(rhs.xxxx());
        res = self.y_axis.mul(rhs.yyyy()).add(res);
        res = self.z_axis.mul(rhs.zzzz()).add(res);
        Vec3A::from_vec4(res)
    }

    /// Transforms a 4D vector.
    #[inline]
    #[must_use]
    pub fn mul_vec4(&self, rhs: Vec4) -> Vec4 {
        let mut res = self.x_axis.mul(rhs.xxxx());
        res = res.add(self.y_axis.mul(rhs.yyyy()));
        res = res.add(self.z_axis.mul(rhs.zzzz()));
        res = res.add(self.w_axis.mul(rhs.wwww()));
        res
    }

    /// Transforms a 4D vector by the transpose of `self`.
    #[inline]
    #[must_use]
    pub fn mul_transpose_vec4(&self, rhs: Vec4) -> Vec4 {
        Vec4::new(
            self.x_axis.dot(rhs),
            self.y_axis.dot(rhs),
            self.z_axis.dot(rhs),
            self.w_axis.dot(rhs),
        )
    }

    /// Multiplies two 4x4 matrices.
    #[inline]
    #[must_use]
    pub fn mul_mat4(&self, rhs: &Self) -> Self {
        self.mul(rhs)
    }

    /// Adds two 4x4 matrices.
    #[inline]
    #[must_use]
    pub fn add_mat4(&self, rhs: &Self) -> Self {
        self.add(rhs)
    }

    /// Subtracts two 4x4 matrices.
    #[inline]
    #[must_use]
    pub fn sub_mat4(&self, rhs: &Self) -> Self {
        self.sub(rhs)
    }

    /// Multiplies a 4x4 matrix by a scalar.
    #[inline]
    #[must_use]
    pub fn mul_scalar(&self, rhs: f32) -> Self {
        Self::from_cols(
            self.x_axis.mul(rhs),
            self.y_axis.mul(rhs),
            self.z_axis.mul(rhs),
            self.w_axis.mul(rhs),
        )
    }

    /// Multiply `self` by a scaling vector `scale`.
    /// This is faster than creating a whole diagonal scaling matrix and then multiplying that.
    /// This operation is commutative.
    #[inline]
    #[must_use]
    pub fn mul_diagonal_scale(&self, scale: Vec4) -> Self {
        Self::from_cols(
            self.x_axis * scale.x,
            self.y_axis * scale.y,
            self.z_axis * scale.z,
            self.w_axis * scale.w,
        )
    }

    /// Divides a 4x4 matrix by a scalar.
    #[inline]
    #[must_use]
    pub fn div_scalar(&self, rhs: f32) -> Self {
        let rhs = Vec4::splat(rhs);
        Self::from_cols(
            self.x_axis.div(rhs),
            self.y_axis.div(rhs),
            self.z_axis.div(rhs),
            self.w_axis.div(rhs),
        )
    }

    /// Returns a matrix containing the reciprocal `1.0/n` of each element of `self`.
    #[inline]
    #[must_use]
    pub fn recip(&self) -> Self {
        Self::from_cols(
            self.x_axis.recip(),
            self.y_axis.recip(),
            self.z_axis.recip(),
            self.w_axis.recip(),
        )
    }

    /// 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 matrices 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 {
        self.x_axis.abs_diff_eq(rhs.x_axis, max_abs_diff)
            && self.y_axis.abs_diff_eq(rhs.y_axis, max_abs_diff)
            && self.z_axis.abs_diff_eq(rhs.z_axis, max_abs_diff)
            && self.w_axis.abs_diff_eq(rhs.w_axis, max_abs_diff)
    }

    /// Takes the absolute value of each element in `self`
    #[inline]
    #[must_use]
    pub fn abs(&self) -> Self {
        Self::from_cols(
            self.x_axis.abs(),
            self.y_axis.abs(),
            self.z_axis.abs(),
            self.w_axis.abs(),
        )
    }

    #[cfg(feature = "f64")]
    #[inline]
    #[must_use]
    pub fn as_dmat4(&self) -> DMat4 {
        DMat4::from_cols(
            self.x_axis.as_dvec4(),
            self.y_axis.as_dvec4(),
            self.z_axis.as_dvec4(),
            self.w_axis.as_dvec4(),
        )
    }
}

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

impl Add for Mat4 {
    type Output = Self;
    #[inline]
    fn add(self, rhs: Self) -> Self {
        Self::from_cols(
            self.x_axis.add(rhs.x_axis),
            self.y_axis.add(rhs.y_axis),
            self.z_axis.add(rhs.z_axis),
            self.w_axis.add(rhs.w_axis),
        )
    }
}

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

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

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

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

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

impl Sub for Mat4 {
    type Output = Self;
    #[inline]
    fn sub(self, rhs: Self) -> Self {
        Self::from_cols(
            self.x_axis.sub(rhs.x_axis),
            self.y_axis.sub(rhs.y_axis),
            self.z_axis.sub(rhs.z_axis),
            self.w_axis.sub(rhs.w_axis),
        )
    }
}

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

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

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

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

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

impl Neg for Mat4 {
    type Output = Self;
    #[inline]
    fn neg(self) -> Self::Output {
        Self::from_cols(
            self.x_axis.neg(),
            self.y_axis.neg(),
            self.z_axis.neg(),
            self.w_axis.neg(),
        )
    }
}

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

impl Mul for Mat4 {
    type Output = Self;
    #[inline]
    fn mul(self, rhs: Self) -> Self {
        Self::from_cols(
            self.mul(rhs.x_axis),
            self.mul(rhs.y_axis),
            self.mul(rhs.z_axis),
            self.mul(rhs.w_axis),
        )
    }
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

impl Div<Mat4> for f32 {
    type Output = Mat4;
    #[inline]
    fn div(self, rhs: Mat4) -> Self::Output {
        Mat4::from_cols(
            self.div(rhs.x_axis),
            self.div(rhs.y_axis),
            self.div(rhs.z_axis),
            self.div(rhs.w_axis),
        )
    }
}

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

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

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

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

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

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

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

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

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

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

impl<'a> Sum<&'a Self> for Mat4 {
    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 Mat4 {
    fn product<I>(iter: I) -> Self
    where
        I: Iterator<Item = Self>,
    {
        iter.fold(Self::IDENTITY, Self::mul)
    }
}

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

impl PartialEq for Mat4 {
    #[inline]
    fn eq(&self, rhs: &Self) -> bool {
        self.x_axis.eq(&rhs.x_axis)
            && self.y_axis.eq(&rhs.y_axis)
            && self.z_axis.eq(&rhs.z_axis)
            && self.w_axis.eq(&rhs.w_axis)
    }
}

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

impl AsMut<[f32; 16]> for Mat4 {
    #[inline]
    fn as_mut(&mut self) -> &mut [f32; 16] {
        unsafe { &mut *(self as *mut Self as *mut [f32; 16]) }
    }
}

impl fmt::Debug for Mat4 {
    fn fmt(&self, fmt: &mut fmt::Formatter<'_>) -> fmt::Result {
        fmt.debug_struct(stringify!(Mat4))
            .field("x_axis", &self.x_axis)
            .field("y_axis", &self.y_axis)
            .field("z_axis", &self.z_axis)
            .field("w_axis", &self.w_axis)
            .finish()
    }
}

impl fmt::Display for Mat4 {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        if let Some(p) = f.precision() {
            write!(
                f,
                "[{:.*}, {:.*}, {:.*}, {:.*}]",
                p, self.x_axis, p, self.y_axis, p, self.z_axis, p, self.w_axis
            )
        } else {
            write!(
                f,
                "[{}, {}, {}, {}]",
                self.x_axis, self.y_axis, self.z_axis, self.w_axis
            )
        }
    }
}