//! # Mat4
//!
//! A 4x4 column-major matrix.
//!
//! # Conventions
//!
//! - Column-major
//! - Fields are named `m_COL_ROW` (e. g. `m01` is column `0`, row `1`)

use std::ops::Mul;

use crate::math::{Quat, Vec3, Vec4};

/// A 4x4 column-major matrix.
///
/// Fields are named `m_COL_ROW` (e.g. `m01` is column `0`, row `1`).
#[derive(Clone, Copy, Debug, PartialEq)]
#[repr(C)]
pub struct Mat4 {
    /// Column 0, row 0 of matrix.
    pub m00: f32,
    /// Column 0, row 1 of matrix.
    pub m01: f32,
    /// Column 0, row 2 of matrix.
    pub m02: f32,
    /// Column 0, row 3 of matrix.
    pub m03: f32,
    /// Column 1, row 0 of matrix.
    pub m10: f32,
    /// Column 1, row 1 of matrix.
    pub m11: f32,
    /// Column 1, row 2 of matrix.
    pub m12: f32,
    /// Column 1, row 3 of matrix.
    pub m13: f32,
    /// Column 2, row 0 of matrix.
    pub m20: f32,
    /// Column 2, row 1 of matrix.
    pub m21: f32,
    /// Column 2, row 2 of matrix.
    pub m22: f32,
    /// Column 2, row 3 of matrix.
    pub m23: f32,
    /// Column 3, row 0 of matrix.
    pub m30: f32,
    /// Column 3, row 1 of matrix.
    pub m31: f32,
    /// Column 3, row 2 of matrix.
    pub m32: f32,
    /// Column 3, row 3 of matrix.
    pub m33: f32,
}

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

    /// A 4x4 matrix with all elements `NaN`.
    pub const NAN: Self = Self::from_cols(Vec4::NAN, Vec4::NAN, Vec4::NAN, Vec4::NAN);

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

    /// A 4x4 matrix, where all diagonal elements are `1.0` and all others are `0.0`.
    #[rustfmt::skip]
    pub const IDENTITY: Self = Self::from_cols(Vec4::X, Vec4::Y, Vec4::Z, Vec4::W);

    /// Creates a `Mat4` from four column `Vec4`.
    #[inline]
    #[must_use]
    #[rustfmt::skip]
    pub const fn from_cols(col0: Vec4, col1: Vec4, col2: Vec4, col3: Vec4) -> Self {
        Self {
            m00: col0.x, m01: col0.y, m02: col0.z, m03: col0.w,
            m10: col1.x, m11: col1.y, m12: col1.z, m13: col1.w,
            m20: col2.x, m21: col2.y, m22: col2.z, m23: col2.w,
            m30: col3.x, m31: col3.y, m32: col3.z, m33: col3.w,
        }
    }

    /// Returns column `0` of `self` as a `Vec4`.
    #[inline]
    #[must_use]
    pub const fn col0(self) -> Vec4 {
        Vec4::new(self.m00, self.m01, self.m02, self.m03)
    }

    /// Returns column `1` of `self` as a `Vec4`.
    #[inline]
    #[must_use]
    pub const fn col1(self) -> Vec4 {
        Vec4::new(self.m10, self.m11, self.m12, self.m13)
    }

    /// Returns column `2` of `self` as a `Vec4`.
    #[inline]
    #[must_use]
    pub const fn col2(self) -> Vec4 {
        Vec4::new(self.m20, self.m21, self.m22, self.m23)
    }

    /// Returns column `3` of `self` as a `Vec4`.
    #[inline]
    #[must_use]
    pub const fn col3(self) -> Vec4 {
        Vec4::new(self.m30, self.m31, self.m32, self.m33)
    }

    /// Returns all columns of `self` as a `[Vec4; 4]`.
    #[inline]
    #[must_use]
    pub const fn cols(self) -> [Vec4; 4] {
        [self.col0(), self.col1(), self.col2(), self.col3()]
    }

    /// Returns row `0` of `self` as a `Vec4`.
    #[inline]
    #[must_use]
    pub const fn row0(self) -> Vec4 {
        Vec4::new(self.m00, self.m10, self.m20, self.m30)
    }

    /// Returns row `1` of `self` as a `Vec4`.
    #[inline]
    #[must_use]
    pub const fn row1(self) -> Vec4 {
        Vec4::new(self.m01, self.m11, self.m21, self.m31)
    }

    /// Returns row `2` of `self` as a `Vec4`.
    #[inline]
    #[must_use]
    pub const fn row2(self) -> Vec4 {
        Vec4::new(self.m02, self.m12, self.m22, self.m32)
    }

    /// Returns row `3` of `self` as a `Vec4`.
    #[inline]
    #[must_use]
    pub const fn row3(self) -> Vec4 {
        Vec4::new(self.m03, self.m13, self.m23, self.m33)
    }

    /// Returns all rows of `self` as a `[Vec4; 4]`.
    #[inline]
    #[must_use]
    pub const fn rows(self) -> [Vec4; 4] {
        [self.row0(), self.row1(), self.row2(), self.row3()]
    }

    /// Creates a `Mat4` that translates by `x`, `y` and `z` components of `v`.
    #[inline]
    #[must_use]
    pub const fn from_translation(v: Vec3) -> Self {
        Self::from_cols(Vec4::X, Vec4::Y, Vec4::Z, Vec4::new(v.x, v.y, v.z, 1.0))
    }

    /// Returns the translation component of `self` as a `Vec3`.
    #[inline]
    #[must_use]
    pub const fn translation(self) -> Vec3 {
        Vec3::new(self.m30, self.m31, self.m32)
    }

    /// Creates a `Mat4` that scales by `x`, `y` and `z` components of `v` along each axis.
    #[inline]
    #[must_use]
    pub const fn from_scale(v: Vec3) -> Self {
        Self::from_cols(
            Vec4::new(v.x, 0.0, 0.0, 0.0),
            Vec4::new(0.0, v.y, 0.0, 0.0),
            Vec4::new(0.0, 0.0, v.z, 0.0),
            Vec4::W,
        )
    }

    /// Creates a `Mat4` that rotates the X axis by `angle` in radians.
    #[inline]
    #[must_use]
    pub fn from_rotation_x(angle: f32) -> Self {
        let (s, c): (f32, f32) = angle.sin_cos();

        Self::from_cols(
            Vec4::X,
            Vec4::new(0.0, c, s, 0.0),
            Vec4::new(0.0, -s, c, 0.0),
            Vec4::W,
        )
    }

    /// Creates a `Mat4` that rotates the X axis by `angle` in degrees.
    #[inline]
    #[must_use]
    pub fn from_rotation_degrees_x(angle: f32) -> Self {
        Self::from_rotation_x(angle.to_radians())
    }

    /// Creates a `Mat4` that rotates the Y axis by `angle` in radians.
    #[inline]
    #[must_use]
    pub fn from_rotation_y(angle: f32) -> Self {
        let (s, c): (f32, f32) = angle.sin_cos();

        Self::from_cols(
            Vec4::new(c, 0.0, -s, 0.0),
            Vec4::Y,
            Vec4::new(s, 0.0, c, 0.0),
            Vec4::W,
        )
    }

    /// Creates a `Mat4` that rotates the Y axis by `angle` in degrees.
    #[inline]
    #[must_use]
    pub fn from_rotation_degrees_y(angle: f32) -> Self {
        Self::from_rotation_y(angle.to_radians())
    }

    /// Creates a `Mat4` that rotates the Z axis by `angle` in radians.
    #[inline]
    #[must_use]
    pub fn from_rotation_z(angle: f32) -> Self {
        let (s, c): (f32, f32) = angle.sin_cos();

        Self::from_cols(
            Vec4::new(c, s, 0.0, 0.0),
            Vec4::new(-s, c, 0.0, 0.0),
            Vec4::Z,
            Vec4::W,
        )
    }

    /// Creates a `Mat4` that rotates the Z axis by `angle` in degrees.
    #[inline]
    #[must_use]
    pub fn from_rotation_degrees_z(angle: f32) -> Self {
        Self::from_rotation_z(angle.to_radians())
    }

    /// Creates a `Mat4` from a `translation`, a `rotation` as euler angles in a `Vec3` in radians and a `scale`.
    ///
    /// Applied in order: scale -> rotation -> translate.
    #[inline]
    #[must_use]
    pub fn from_trs_euler(translation: Vec3, rotation: Vec3, scale: Vec3) -> Self {
        let t: Self = Self::from_translation(translation);
        let r: Self = Self::from_rotation_x(rotation.x)
            * Self::from_rotation_y(rotation.y)
            * Self::from_rotation_z(rotation.z);
        let s: Self = Self::from_scale(scale);

        t * r * s
    }

    /// Creates a `Mat4` from a `translation`, a `rotation` as euler angles in a `Vec3` in degrees and a `scale`.
    ///
    /// Applied in order: scale -> rotation -> translate.
    #[inline]
    #[must_use]
    pub fn from_trs_euler_degrees(translation: Vec3, rotation: Vec3, scale: Vec3) -> Self {
        let t: Self = Self::from_translation(translation);
        let r: Self = Self::from_rotation_degrees_x(rotation.x)
            * Self::from_rotation_degrees_y(rotation.y)
            * Self::from_rotation_degrees_z(rotation.z);
        let s: Self = Self::from_scale(scale);

        t * r * s
    }

    /// Creates a `Mat4` from a `translation`, a `rotation` as a `Quat` and a `scale`.
    ///
    /// Applied in order: scale -> rotate -> translate.
    #[inline]
    #[must_use]
    pub fn from_trs(translation: Vec3, rotation: Quat, scale: Vec3) -> Self {
        let t: Self = Self::from_translation(translation);
        let r: Self = Self::from_quat(rotation);
        let s: Self = Self::from_scale(scale);

        t * r * s
    }

    /// Creates a `Mat4` from a `Quat`.
    ///
    /// # Precondition
    /// Expects `q` to be normalized.
    #[inline]
    #[must_use]
    pub const fn from_quat(q: Quat) -> Self {
        let x2: f32 = q.x * q.x;
        let y2: f32 = q.y * q.y;
        let z2: f32 = q.z * q.z;
        let xy: f32 = q.x * q.y;
        let xz: f32 = q.x * q.z;
        let yz: f32 = q.y * q.z;
        let wx: f32 = q.w * q.x;
        let wy: f32 = q.w * q.y;
        let wz: f32 = q.w * q.z;

        Self {
            // Column 0
            m00: 1.0 - 2.0 * (y2 + z2),
            m01: 2.0 * (xy + wz),
            m02: 2.0 * (xz - wy),
            m03: 0.0,
            // Column 1
            m10: 2.0 * (xy - wz),
            m11: 1.0 - 2.0 * (x2 + z2),
            m12: 2.0 * (yz + wx),
            m13: 0.0,
            // Column 2
            m20: 2.0 * (xz + wy),
            m21: 2.0 * (yz - wx),
            m22: 1.0 - 2.0 * (x2 + y2),
            m23: 0.0,
            // Column 3
            m30: 0.0,
            m31: 0.0,
            m32: 0.0,
            m33: 1.0,
        }
    }

    /// Creates a `Mat4` that rotates around `axis` by `angle` in radians.
    ///
    /// # Precondition
    /// Expects `axis` to be normalized.
    #[inline]
    #[must_use]
    pub fn from_axis_angle(axis: Vec3, angle: f32) -> Self {
        Self::from_quat(Quat::from_axis_angle(axis, angle))
    }

    /// Creates a `Mat4` that rotates around `axis` by `angle` in degrees.
    ///
    /// # Precondition
    /// Expects `axis` to be normalized.
    #[inline]
    #[must_use]
    pub fn from_axis_angle_degrees(axis: Vec3, angle: f32) -> Self {
        Self::from_quat(Quat::from_axis_angle_degrees(axis, angle))
    }

    /// Creates a `Mat4` that is a view matrix
    /// positioned at `eye`, looking at `center`, with `up` as the up direction.
    ///
    /// # Convention
    /// Uses right-handed coordinate system.
    ///
    /// # Precondition
    /// Expects `up` to be normalized and not parallel to `(center - eye)`.
    ///
    /// # Edge cases
    /// If `eye` and `center` are the same point, returns [`IDENTITY`](Self::IDENTITY).
    #[inline]
    #[must_use]
    pub fn look_at_rh(eye: Vec3, center: Vec3, up: Vec3) -> Self {
        let dir: Vec3 = center - eye;

        if dir.length_squared() < Self::EPSILON_SQ {
            return Self::IDENTITY;
        }

        let f: Vec3 = dir.normalized();
        let r: Vec3 = f.cross(up).normalized();
        let u: Vec3 = r.cross(f);

        Self::from_cols(
            Vec4::new(r.x, u.x, -f.x, 0.0),
            Vec4::new(r.y, u.y, -f.y, 0.0),
            Vec4::new(r.z, u.z, -f.z, 0.0),
            Vec4::new(-r.dot(eye), -u.dot(eye), f.dot(eye), 1.0),
        )
    }

    /// Creates a `Mat4` that is a view matrix
    /// positioned at `eye`, looking at `center`, with `up` as the up direction.
    ///
    /// # Convention
    /// Uses left-handed coordinate system.
    ///
    /// # Precondition
    /// Expects `up` to be normalized and not parallel to `(center - eye)`.
    ///
    /// # Edge cases
    /// If `eye` and `center` are the same point, returns [`IDENTITY`](Self::IDENTITY).
    #[inline]
    #[must_use]
    pub fn look_at_lh(eye: Vec3, center: Vec3, up: Vec3) -> Self {
        let dir: Vec3 = center - eye;

        if dir.length_squared() < Self::EPSILON_SQ {
            return Self::IDENTITY;
        }

        let f: Vec3 = dir.normalized();
        let r: Vec3 = up.cross(f).normalized();
        let u: Vec3 = f.cross(r);

        Self::from_cols(
            Vec4::new(r.x, u.x, f.x, 0.0),
            Vec4::new(r.y, u.y, f.y, 0.0),
            Vec4::new(r.z, u.z, f.z, 0.0),
            Vec4::new(-r.dot(eye), -u.dot(eye), -f.dot(eye), 1.0),
        )
    }

    /// Creates a `Mat4` that is a perspective projection matrix with vertical field of view `fov_y` in radians,
    /// aspect ratio `aspect` (width / height) and near/far clip planes `z_near` and `z_far`.
    ///
    /// # Convention
    /// Uses right-handed coordinate system. Depth maps to `[0, 1]`.
    ///
    /// # Precondition
    /// Expects `fov_y > 0.0`, `aspect > 0.0` and `z_near < z_far`.
    #[inline]
    #[must_use]
    pub fn perspective_rh(fov_y: f32, aspect: f32, z_near: f32, z_far: f32) -> Self {
        let tan_half_fov: f32 = (fov_y * 0.5).tan();
        let sy: f32 = 1.0 / tan_half_fov;
        let sx: f32 = sy / aspect;
        let range: f32 = z_far - z_near;

        Self::from_cols(
            Vec4::new(sx, 0.0, 0.0, 0.0),
            Vec4::new(0.0, sy, 0.0, 0.0),
            Vec4::new(0.0, 0.0, -z_far / range, -1.0),
            Vec4::new(0.0, 0.0, -(z_near * z_far) / range, 0.0),
        )
    }

    /// Creates a `Mat4` that is a perspective projection matrix with vertical field of view `fov_y` in radians,
    /// aspect ratio `aspect` (width / height) and near/far clip planes `z_near` and `z_far`.
    ///
    /// # Convention
    /// Uses left-handed coordinate system. Depth maps to `[0, 1]`.
    ///
    /// # Precondition
    /// Expects `fov_y > 0.0`, `aspect > 0.0` and `z_near < z_far`.
    #[inline]
    #[must_use]
    pub fn perspective_lh(fov_y: f32, aspect: f32, z_near: f32, z_far: f32) -> Self {
        let tan_half_fov: f32 = (fov_y * 0.5).tan();
        let sy: f32 = 1.0 / tan_half_fov;
        let sx: f32 = sy / aspect;
        let range: f32 = z_far - z_near;

        Self::from_cols(
            Vec4::new(sx, 0.0, 0.0, 0.0),
            Vec4::new(0.0, sy, 0.0, 0.0),
            Vec4::new(0.0, 0.0, z_far / range, 1.0),
            Vec4::new(0.0, 0.0, -(z_near * z_far) / range, 0.0),
        )
    }

    /// Creates a `Mat4` that is a perspective projection matrix with vertical field of view `fov_y` in radians,
    /// aspect ratio `aspect` (width / height) and near/far clip planes `z_near` and `z_far`.
    ///
    /// # Convention
    /// Uses right-handed coordinate system. Depth maps to `[-1, 1]`.
    ///
    /// # Precondition
    /// Expects `fov_y > 0.0`, `aspect > 0.0` and `z_near < z_far`.
    #[inline]
    #[must_use]
    pub fn perspective_rh_gl(fov_y: f32, aspect: f32, z_near: f32, z_far: f32) -> Self {
        let tan_half_fov: f32 = (fov_y * 0.5).tan();
        let sy: f32 = 1.0 / tan_half_fov;
        let sx: f32 = sy / aspect;
        let range: f32 = z_far - z_near;

        Self::from_cols(
            Vec4::new(sx, 0.0, 0.0, 0.0),
            Vec4::new(0.0, sy, 0.0, 0.0),
            Vec4::new(0.0, 0.0, -(z_far + z_near) / range, -1.0),
            Vec4::new(0.0, 0.0, -(2.0 * z_far * z_near) / range, 0.0),
        )
    }

    /// Creates a `Mat4` that is a perspective projection matrix with vertical field of view `fov_y` in radians,
    /// aspect ratio `aspect` (width / height) and near/far clip planes `z_near` and `z_far`.
    ///
    /// # Convention
    /// Uses left-handed coordinate system. Depth maps to `[-1, 1]`.
    ///
    /// # Precondition
    /// Expects `fov_y > 0.0`, `aspect > 0.0` and `z_near < z_far`.
    #[inline]
    #[must_use]
    pub fn perspective_lh_gl(fov_y: f32, aspect: f32, z_near: f32, z_far: f32) -> Self {
        let tan_half_fov: f32 = (fov_y * 0.5).tan();
        let sy: f32 = 1.0 / tan_half_fov;
        let sx: f32 = sy / aspect;
        let range: f32 = z_far - z_near;

        Self::from_cols(
            Vec4::new(sx, 0.0, 0.0, 0.0),
            Vec4::new(0.0, sy, 0.0, 0.0),
            Vec4::new(0.0, 0.0, (z_far + z_near) / range, 1.0),
            Vec4::new(0.0, 0.0, -(2.0 * z_far * z_near) / range, 0.0),
        )
    }

    /// Creates a `Mat4` that is an orthographic projection matrix with the given
    /// left, right, bottom, top, near and far clip planes.
    ///
    /// # Convention
    /// Uses right-handed coordinate system. Depth maps to `[0, 1]`.
    ///
    /// # Precondition
    /// Expects `left != right`, `bottom != top` and `z_near != z_far`.
    #[inline]
    #[must_use]
    pub const fn orthographic_rh(
        left: f32,
        right: f32,
        bottom: f32,
        top: f32,
        z_near: f32,
        z_far: f32,
    ) -> Self {
        let rcp_width: f32 = 1.0 / (right - left);
        let rcp_height: f32 = 1.0 / (top - bottom);
        let rcp_range: f32 = 1.0 / (z_near - z_far);

        Self::from_cols(
            Vec4::new(2.0 * rcp_width, 0.0, 0.0, 0.0),
            Vec4::new(0.0, 2.0 * rcp_height, 0.0, 0.0),
            Vec4::new(0.0, 0.0, rcp_range, 0.0),
            Vec4::new(
                -(right + left) * rcp_width,
                -(top + bottom) * rcp_height,
                z_near * rcp_range,
                1.0,
            ),
        )
    }

    /// Creates a `Mat4` that is an orthographic projection matrix with the given
    /// left, right, bottom, top, near and far clip planes.
    ///
    /// # Convention
    /// Uses left-handed coordinate system. Depth maps to `[0, 1]`.
    ///
    /// # Precondition
    /// Expects `left != right`, `bottom != top` and `z_near != z_far`.
    #[inline]
    #[must_use]
    pub const fn orthographic_lh(
        left: f32,
        right: f32,
        bottom: f32,
        top: f32,
        z_near: f32,
        z_far: f32,
    ) -> Self {
        let rcp_width: f32 = 1.0 / (right - left);
        let rcp_height: f32 = 1.0 / (top - bottom);
        let rcp_range: f32 = 1.0 / (z_far - z_near);

        Self::from_cols(
            Vec4::new(2.0 * rcp_width, 0.0, 0.0, 0.0),
            Vec4::new(0.0, 2.0 * rcp_height, 0.0, 0.0),
            Vec4::new(0.0, 0.0, rcp_range, 0.0),
            Vec4::new(
                -(right + left) * rcp_width,
                -(top + bottom) * rcp_height,
                -z_near * rcp_range,
                1.0,
            ),
        )
    }

    /// Creates a `Mat4` that is an orthographic projection matrix with the given
    /// left, right, bottom, top, near and far clip planes.
    ///
    /// # Convention
    /// Uses right-handed coordinate system. Depth maps to `[-1, 1]`.
    ///
    /// # Precondition
    /// Expects `left != right`, `bottom != top` and `z_near != z_far`.
    #[inline]
    #[must_use]
    pub const fn orthographic_rh_gl(
        left: f32,
        right: f32,
        bottom: f32,
        top: f32,
        z_near: f32,
        z_far: f32,
    ) -> Self {
        let rcp_width: f32 = 1.0 / (right - left);
        let rcp_height: f32 = 1.0 / (top - bottom);
        let rcp_range: f32 = 1.0 / (z_far - z_near);

        Self::from_cols(
            Vec4::new(2.0 * rcp_width, 0.0, 0.0, 0.0),
            Vec4::new(0.0, 2.0 * rcp_height, 0.0, 0.0),
            Vec4::new(0.0, 0.0, -2.0 * rcp_range, 0.0),
            Vec4::new(
                -(right + left) * rcp_width,
                -(top + bottom) * rcp_height,
                -(z_far + z_near) * rcp_range,
                1.0,
            ),
        )
    }

    /// Creates a `Mat4` that is an orthographic projection matrix with the given
    /// left, right, bottom, top, near and far clip planes.
    ///
    /// # Convention
    /// Uses left-handed coordinate system. Depth maps to `[-1, 1]`.
    ///
    /// # Precondition
    /// Expects `left != right`, `bottom != top` and `z_near != z_far`.
    #[inline]
    #[must_use]
    pub const fn orthographic_lh_gl(
        left: f32,
        right: f32,
        bottom: f32,
        top: f32,
        z_near: f32,
        z_far: f32,
    ) -> Self {
        let rcp_width: f32 = 1.0 / (right - left);
        let rcp_height: f32 = 1.0 / (top - bottom);
        let rcp_range: f32 = 1.0 / (z_far - z_near);

        Self::from_cols(
            Vec4::new(2.0 * rcp_width, 0.0, 0.0, 0.0),
            Vec4::new(0.0, 2.0 * rcp_height, 0.0, 0.0),
            Vec4::new(0.0, 0.0, 2.0 * rcp_range, 0.0),
            Vec4::new(
                -(right + left) * rcp_width,
                -(top + bottom) * rcp_height,
                -(z_far + z_near) * rcp_range,
                1.0,
            ),
        )
    }

    /// Swaps rows and columns of `self`, turning column `n` into row `n`.
    #[inline]
    pub const fn transpose(&mut self) {
        *self = self.transposed();
    }

    /// Returns a `Mat4` that is the transpose of `self`.
    #[inline]
    #[must_use]
    #[rustfmt::skip]
    pub const fn transposed(self) -> Self {
        Self {
            m00: self.m00, m01: self.m10, m02: self.m20, m03: self.m30,
            m10: self.m01, m11: self.m11, m12: self.m21, m13: self.m31,
            m20: self.m02, m21: self.m12, m22: self.m22, m23: self.m32,
            m30: self.m03, m31: self.m13, m32: self.m23, m33: self.m33,
        }
    }

    /// Returns the determinant of `self`.
    #[must_use]
    pub const fn determinant(self) -> f32 {
        // 2x2 sub-determinants reused for cofactor expansion
        let a: f32 = self.m22 * self.m33 - self.m23 * self.m32;
        let b: f32 = self.m21 * self.m33 - self.m23 * self.m31;
        let c: f32 = self.m21 * self.m32 - self.m22 * self.m31;
        let d: f32 = self.m20 * self.m33 - self.m23 * self.m30;
        let e: f32 = self.m20 * self.m32 - self.m22 * self.m30;
        let f: f32 = self.m20 * self.m31 - self.m21 * self.m30;

        // Cofactors along row 0
        let c00: f32 = self.m11 * a - self.m12 * b + self.m13 * c;
        let c10: f32 = -(self.m10 * a - self.m12 * d + self.m13 * e);
        let c20: f32 = self.m10 * b - self.m11 * d + self.m13 * f;
        let c30: f32 = -(self.m10 * c - self.m11 * e + self.m12 * f);

        self.m00 * c00 + self.m01 * c10 + self.m02 * c20 + self.m03 * c30
    }

    /// Inverts `self`.
    ///
    /// # Edge cases
    /// If `self` is singular, `self` stays same.
    pub const fn inverse(&mut self) {
        if let Some(inv) = self.inversed() {
            *self = inv;
        }
    }

    /// Returns a `Mat4` that is the inverse of `self` or `None` if `self` is singular.
    ///
    /// # Edge cases
    /// If `self` is singular, returns `None`.
    #[must_use]
    pub const fn inversed(self) -> Option<Self> {
        // 2x2 minors from bottom-right 3x3 block
        let a: f32 = self.m22 * self.m33 - self.m23 * self.m32;
        let b: f32 = self.m21 * self.m33 - self.m23 * self.m31;
        let c: f32 = self.m21 * self.m32 - self.m22 * self.m31;
        let d: f32 = self.m20 * self.m33 - self.m23 * self.m30;
        let e: f32 = self.m20 * self.m32 - self.m22 * self.m30;
        let f: f32 = self.m20 * self.m31 - self.m21 * self.m30;

        // 2x2 minors from top-right 3x3 block
        let g: f32 = self.m12 * self.m33 - self.m13 * self.m32;
        let h: f32 = self.m11 * self.m33 - self.m13 * self.m31;
        let i: f32 = self.m11 * self.m32 - self.m12 * self.m31;
        let j: f32 = self.m10 * self.m33 - self.m13 * self.m30;
        let k: f32 = self.m10 * self.m32 - self.m12 * self.m30;
        let l: f32 = self.m10 * self.m31 - self.m11 * self.m30;

        // 2x2 minors from top-left 3x3 block
        let m: f32 = self.m12 * self.m23 - self.m13 * self.m22;
        let n: f32 = self.m11 * self.m23 - self.m13 * self.m21;
        let o: f32 = self.m11 * self.m22 - self.m12 * self.m21;
        let p: f32 = self.m10 * self.m23 - self.m13 * self.m20;
        let q: f32 = self.m10 * self.m22 - self.m12 * self.m20;
        let r: f32 = self.m10 * self.m21 - self.m11 * self.m20;

        // cofactors
        let c00: f32 = self.m11 * a - self.m12 * b + self.m13 * c;
        let c10: f32 = -(self.m10 * a - self.m12 * d + self.m13 * e);
        let c20: f32 = self.m10 * b - self.m11 * d + self.m13 * f;
        let c30: f32 = -(self.m10 * c - self.m11 * e + self.m12 * f);

        let det: f32 = self.m00 * c00 + self.m01 * c10 + self.m02 * c20 + self.m03 * c30;

        if det.abs() < Self::EPSILON {
            return None;
        }

        let inv_det: f32 = 1.0 / det;

        Some(Self {
            // Column 0
            m00: c00 * inv_det,
            m01: -(self.m01 * a - self.m02 * b + self.m03 * c) * inv_det,
            m02: (self.m01 * g - self.m02 * h + self.m03 * i) * inv_det,
            m03: -(self.m01 * m - self.m02 * n + self.m03 * o) * inv_det,

            // Column 1
            m10: c10 * inv_det,
            m11: (self.m00 * a - self.m02 * d + self.m03 * e) * inv_det,
            m12: -(self.m00 * g - self.m02 * j + self.m03 * k) * inv_det,
            m13: (self.m00 * m - self.m02 * p + self.m03 * q) * inv_det,

            // Column 2
            m20: c20 * inv_det,
            m21: -(self.m00 * b - self.m01 * d + self.m03 * f) * inv_det,
            m22: (self.m00 * h - self.m01 * j + self.m03 * l) * inv_det,
            m23: -(self.m00 * n - self.m01 * p + self.m03 * r) * inv_det,

            // Column 3
            m30: c30 * inv_det,
            m31: (self.m00 * c - self.m01 * e + self.m02 * f) * inv_det,
            m32: -(self.m00 * i - self.m01 * k + self.m02 * l) * inv_det,
            m33: (self.m00 * o - self.m01 * q + self.m02 * r) * inv_det,
        })
    }

    /// Returns `true` is `self` is invertible.
    ///
    /// Uses [`EPSILON`](Self::EPSILON) to determine whether the determinant
    /// is sufficiently far from zero.
    #[inline]
    #[must_use]
    pub const fn is_invertible(self) -> bool {
        self.determinant().abs() >= Self::EPSILON
    }

    /// Transform `v` as a point by `self`, applying rotation, scale and translation.
    #[inline]
    #[must_use]
    pub const fn transform_point(self, v: Vec3) -> Vec3 {
        Vec3::new(
            self.m00 * v.x + self.m10 * v.y + self.m20 * v.z + self.m30,
            self.m01 * v.x + self.m11 * v.y + self.m21 * v.z + self.m31,
            self.m02 * v.x + self.m12 * v.y + self.m22 * v.z + self.m32,
        )
    }

    /// Transform `v` as a direction by `self`, applying rotation and scale.
    #[inline]
    #[must_use]
    pub const fn transform_direction(self, v: Vec3) -> Vec3 {
        Vec3::new(
            self.m00 * v.x + self.m10 * v.y + self.m20 * v.z,
            self.m01 * v.x + self.m11 * v.y + self.m21 * v.z,
            self.m02 * v.x + self.m12 * v.y + self.m22 * v.z,
        )
    }

    /// Returns `true` if all components of `self` are approximately equal to
    /// the equivalent component of `other`.
    ///
    /// Uses [`EPSILON`](Self::EPSILON) as threshold.
    #[inline]
    #[must_use]
    pub fn approx_eq(self, other: Self) -> bool {
        let a: [f32; 16] = self.into();
        let b: [f32; 16] = other.into();
        a.iter()
            .zip(b.iter())
            .all(|(x, y)| (x - y).abs() < Self::EPSILON)
    }

    /// Returns `true` if all component of `self` are finite.
    #[inline]
    #[must_use]
    pub fn is_finite(self) -> bool {
        let m: [f32; 16] = self.into();
        m.iter().all(|x| x.is_finite())
    }

    /// Returns `true` if any component of `self` is `NaN`.
    #[inline]
    #[must_use]
    pub fn is_nan(self) -> bool {
        let m: [f32; 16] = self.into();
        m.iter().any(|x| x.is_nan())
    }
}

impl Mul for Mat4 {
    type Output = Self;

    /// Multiplies each column of `rhs` by `self`, returns a new `Mat4`.
    #[inline]
    fn mul(self, rhs: Self) -> Self::Output {
        Self::from_cols(
            self * rhs.col0(),
            self * rhs.col1(),
            self * rhs.col2(),
            self * rhs.col3(),
        )
    }
}

impl Mul<Vec4> for Mat4 {
    type Output = Vec4;

    /// Multiplies `self` by `rhs`, returns a transformed `Vec4`.
    #[inline]
    fn mul(self, rhs: Vec4) -> Self::Output {
        Vec4::new(
            self.row0().dot(rhs),
            self.row1().dot(rhs),
            self.row2().dot(rhs),
            self.row3().dot(rhs),
        )
    }
}

impl From<[f32; 16]> for Mat4 {
    /// Creates a `Mat4` from a `[f32; 16]`.
    /// Elements map to columns in order, each column maps to `[x, y, z, w]`.
    ///
    /// # Convention
    /// Converts in column-major order.
    #[inline]
    fn from(v: [f32; 16]) -> Self {
        Self::from_cols(
            Vec4::new(v[0], v[1], v[2], v[3]),
            Vec4::new(v[4], v[5], v[6], v[7]),
            Vec4::new(v[8], v[9], v[10], v[11]),
            Vec4::new(v[12], v[13], v[14], v[15]),
        )
    }
}

impl From<Mat4> for [f32; 16] {
    /// Creates a `[f32; 16]` from a `Mat4`.
    /// Elements map to `[col0, col1, col2, col3]`, each column maps to `[x, y, z, w]`.
    ///
    /// # Convention
    /// Converts in column-major order.
    #[inline]
    fn from(v: Mat4) -> Self {
        [
            v.m00, v.m01, v.m02, v.m03, v.m10, v.m11, v.m12, v.m13, v.m20, v.m21, v.m22, v.m23,
            v.m30, v.m31, v.m32, v.m33,
        ]
    }
}

impl From<&[f32]> for Mat4 {
    /// Creates a `Mat4` from a `&[f32]`.
    /// First 16 elements map to `[col0, col1, col2, col3], each column maps to `[x, y, z, w]`.
    ///
    /// # Panics
    /// Panics if `v` is less than 16 elements long.
    ///
    /// # Convention
    /// Converts in column-major order.
    #[inline]
    fn from(v: &[f32]) -> Self {
        Self::from_cols(
            Vec4::new(v[0], v[1], v[2], v[3]),
            Vec4::new(v[4], v[5], v[6], v[7]),
            Vec4::new(v[8], v[9], v[10], v[11]),
            Vec4::new(v[12], v[13], v[14], v[15]),
        )
    }
}

impl AsRef<[f32; 16]> for Mat4 {
    /// Returns a `&[f32; 16]` reference to the matrix 4x4 components.
    /// Elements map to `[col0, col1, col2, col3]`, each column maps to `[x, y, z, w]`.
    ///
    /// # Convention
    /// Array is in column-major order.
    #[inline]
    fn as_ref(&self) -> &[f32; 16] {
        unsafe { &*(self as *const Mat4 as *const [f32; 16]) }
    }
}

impl AsRef<[f32]> for Mat4 {
    /// Returns a `&[f32]` reference to the `Mat4` components.
    /// Elements map to `[col0, col1, col2, col3]`, each column maps to `[x, y, z, w]`.
    #[inline]
    fn as_ref(&self) -> &[f32] {
        let arr: &[f32; 16] = self.as_ref();
        &arr[..]
    }
}

impl AsMut<[f32; 16]> for Mat4 {
    /// Returns a `&mut [f32; 16]` mutable reference to the matrix 4x4 components.
    /// Elements map to `[col0, col1, col2, col3]`, each column maps to `[x, y, z, w]`.
    ///
    /// # Convention
    /// Array is in column-major order.
    #[inline]
    fn as_mut(&mut self) -> &mut [f32; 16] {
        unsafe { &mut *(self as *mut Mat4 as *mut [f32; 16]) }
    }
}

impl AsMut<[f32]> for Mat4 {
    /// Returns a `&[f32]` mutable reference to the `Mat4` components.
    /// Elements map to `[col0, col1, col2, col3]`, each column maps to `[x, y, z, w]`.
    #[inline]
    fn as_mut(&mut self) -> &mut [f32] {
        let arr: &mut [f32; 16] = self.as_mut();
        &mut arr[..]
    }
}

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