//! # Mat4
//!
//! A 4x4 column-major matrix.
//!
//! # Conventions
//!
//! - Column-major
//! - Stored as 4 columns `col0`..`col3`, each a [`Vec4`]
//!
//! # Examples
//!
//! ```rust
//! # use roofy::{Mat4, Vec3, Vec4};
//!
//! let m: Mat4 = Mat4::IDENTITY;
//! let transformed: Mat4 = m * Mat4::from_translation(Vec3::new(1.0, 2.0, 3.0));
//!
//! assert_eq!(transformed.col3, Vec4::new(1.0, 2.0, 3.0, 1.0));
//! ```

use std::{
    fmt::Display,
    ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign},
};

use crate::{Affine3, Mat3, Quat, Vec3, Vec4};

/// A 4x4 column-major matrix.
///
/// Stored as 4 columns `col0`..`col3`, each a [`Vec4`]
#[derive(Clone, Copy, Debug, PartialEq)]
#[repr(C)]
pub struct Mat4 {
    /// Column 0 of the matrix.
    pub col0: Vec4,
    /// Column 1 of the matrix.
    pub col1: Vec4,
    /// Column 2 of the matrix.
    pub col2: Vec4,
    /// Column 3 of the matrix.
    pub col3: Vec4,
}

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 components `NaN`.
    pub const NAN: Self = Self::from_cols(Vec4::NAN, Vec4::NAN, Vec4::NAN, Vec4::NAN);

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

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

    /// Creates a `Mat4` from 16 `f32`, in column-major order.
    ///
    /// # Convention
    /// Parameters are named `mCR` where `C` is the column index and `R` is the row index.
    #[allow(clippy::too_many_arguments)]
    #[inline(always)]
    #[must_use]
    pub 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::from_cols(
            Vec4::new(m00, m01, m02, m03),
            Vec4::new(m10, m11, m12, m13),
            Vec4::new(m20, m21, m22, m23),
            Vec4::new(m30, m31, m32, m33),
        )
    }

    /// Creates a `Mat4` from 4 column `Vec4`.
    #[inline(always)]
    #[must_use]
    pub const fn from_cols(col0: Vec4, col1: Vec4, col2: Vec4, col3: Vec4) -> Self {
        Self {
            col0,
            col1,
            col2,
            col3,
        }
    }

    /// Creates a `Mat4` with diagonal set to `v` and all other components `0.0`.
    #[inline]
    #[must_use]
    pub const fn from_diagonal(v: Vec4) -> Self {
        Self::new(
            v.x, 0.0, 0.0, 0.0, 0.0, v.y, 0.0, 0.0, 0.0, 0.0, v.z, 0.0, 0.0, 0.0, 0.0, v.w,
        )
    }

    /// Creates a `Mat4` from a `[f32; 16]`, in column-major order:
    /// `[col0.x, col0.y, col0.z, col0.w, col1.x, ..., col3.w]`.
    #[inline]
    #[must_use]
    pub const fn from_array(a: [f32; 16]) -> Self {
        Self::new(
            a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8], a[9], a[10], a[11], a[12], a[13],
            a[14], a[15],
        )
    }

    /// Returns a copy of the components as `[col0, col1, col2, col3]`,
    /// with each column as `[x, y, z, w]`, in column-major order.
    #[inline]
    #[must_use]
    pub const fn to_array(self) -> [f32; 16] {
        [
            self.col0.x,
            self.col0.y,
            self.col0.z,
            self.col0.w,
            self.col1.x,
            self.col1.y,
            self.col1.z,
            self.col1.w,
            self.col2.x,
            self.col2.y,
            self.col2.z,
            self.col2.w,
            self.col3.x,
            self.col3.y,
            self.col3.z,
            self.col3.w,
        ]
    }

    /// Creates a `Mat4` from first 16 values of `slice`, in column-major order.
    ///
    /// # Panics
    /// Panics if `slice` is less than 16 elements long.
    #[inline]
    #[must_use]
    pub const fn from_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 of `slice`.
    ///
    /// # Panics
    /// Panics if `slice` is less than 16 elements long.
    #[inline]
    pub const fn write_to_slice(self, slice: &mut [f32]) {
        slice[0] = self.col0.x;
        slice[1] = self.col0.y;
        slice[2] = self.col0.z;
        slice[3] = self.col0.w;
        slice[4] = self.col1.x;
        slice[5] = self.col1.y;
        slice[6] = self.col1.z;
        slice[7] = self.col1.w;
        slice[8] = self.col2.x;
        slice[9] = self.col2.y;
        slice[10] = self.col2.z;
        slice[11] = self.col2.w;
        slice[12] = self.col3.x;
        slice[13] = self.col3.y;
        slice[14] = self.col3.z;
        slice[15] = self.col3.w;
    }

    /// Creates a `Mat4` from a `Mat3`, setting `col3` and `row3` to identity.
    #[inline]
    #[must_use]
    pub const fn from_mat3(m: Mat3) -> Self {
        Self::from_cols(
            m.col0.to_vec4(),
            m.col1.to_vec4(),
            m.col2.to_vec4(),
            Vec4::W,
        )
    }

    /// Creates a `Mat4` from a `Mat3` and a `translation`,
    /// placing `matrix` in the upper-left 3x3 and `translation` in `col3`,
    /// with `row3` set to identity.
    #[inline]
    #[must_use]
    pub const fn from_mat3_translation(matrix: Mat3, translation: Vec3) -> Self {
        Self::from_cols(
            matrix.col0.to_vec4(),
            matrix.col1.to_vec4(),
            matrix.col2.to_vec4(),
            Vec4::new(translation.x, translation.y, translation.z, 1.0),
        )
    }

    /// 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))
    }

    /// 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 around 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 around the X axis by `angle` in degrees.
    #[inline]
    #[must_use]
    pub fn from_rotation_x_degrees(angle: f32) -> Self {
        Self::from_rotation_x(angle.to_radians())
    }

    /// Creates a `Mat4` that rotates around 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 around the Y axis by `angle` in degrees.
    #[inline]
    #[must_use]
    pub fn from_rotation_y_degrees(angle: f32) -> Self {
        Self::from_rotation_y(angle.to_radians())
    }

    /// Creates a `Mat4` that rotates around 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 around the Z axis by `angle` in degrees.
    #[inline]
    #[must_use]
    pub fn from_rotation_z_degrees(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_x_degrees(rotation.x)
            * Self::from_rotation_y_degrees(rotation.y)
            * Self::from_rotation_z_degrees(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` that rotates by `Quat`.
    ///
    /// # Preconditions
    /// 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::from_cols(
            Vec4::new(1.0 - 2.0 * (y2 + z2), 2.0 * (xy + wz), 2.0 * (xz - wy), 0.0),
            Vec4::new(2.0 * (xy - wz), 1.0 - 2.0 * (x2 + z2), 2.0 * (yz + wx), 0.0),
            Vec4::new(2.0 * (xz + wy), 2.0 * (yz - wx), 1.0 - 2.0 * (x2 + y2), 0.0),
            Vec4::new(0.0, 0.0, 0.0, 1.0),
        )
    }

    /// Creates a `Mat4` that rotates around `axis` by `angle` in radians.
    ///
    /// # Preconditions
    /// 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.
    ///
    /// # Preconditions
    /// 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 to `dir`, with `up` as the up direction.
    ///
    /// # Convention
    /// Uses right-handed coordinate system.
    ///
    /// # Preconditions
    /// Expects `dir` to be normalized.
    /// Expects the length of `dir.cross(up)` to not be approximately `0.0` and finite.
    #[inline]
    #[must_use]
    pub fn look_to_rh(eye: Vec3, dir: Vec3, up: Vec3) -> Self {
        let f: Vec3 = dir;
        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 to `dir`, with `up` as the up direction.
    ///
    /// # Convention
    /// Uses left-handed coordinate system.
    ///
    /// # Preconditions
    /// Expects `dir` to be normalized.
    /// Expects the length of `dir.cross(up)` to not be approximately `0.0` and finite.
    #[inline]
    #[must_use]
    pub fn look_to_lh(eye: Vec3, dir: Vec3, up: Vec3) -> Self {
        Self::look_to_rh(eye, -dir, up)
    }

    /// 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.
    ///
    /// # Preconditions
    /// Expects the length of `(center - eye)` to not be approximately `0.0` and finite.
    /// Expects the length of `(center - eye).normalized().cross(up)` to not be approximately `0.0` and finite.
    #[inline]
    #[must_use]
    pub fn look_at_rh(eye: Vec3, center: Vec3, up: Vec3) -> Self {
        Self::look_to_rh(eye, (center - eye).normalized(), up)
    }

    /// 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.
    ///
    /// # Preconditions
    /// Expects the length of `(center - eye)` to not be approximately `0.0` and finite.
    /// Expects the length of `(center - eye).normalized().cross(up)` to not be approximately `0.0` and finite.
    #[inline]
    #[must_use]
    pub fn look_at_lh(eye: Vec3, center: Vec3, up: Vec3) -> Self {
        Self::look_to_lh(eye, (center - eye).normalized(), up)
    }

    /// Creates a right-handed perspective projection matrix with `[0, 1]` depth range.
    ///
    /// # Preconditions
    /// Expects `fov_y_radians`, `aspect_ratio` and `z_near` to be positive.
    /// Expects `z_far > z_near`.
    #[inline]
    #[must_use]
    pub fn perspective_rh(fov_y_radians: f32, aspect_ratio: f32, z_near: f32, z_far: f32) -> Self {
        let f: f32 = 1.0 / (fov_y_radians * 0.5).tan();
        let inv: f32 = 1.0 / (z_near - z_far);

        Self::from_cols(
            Vec4::new(f / aspect_ratio, 0.0, 0.0, 0.0),
            Vec4::new(0.0, f, 0.0, 0.0),
            Vec4::new(0.0, 0.0, z_far * inv, -1.0),
            Vec4::new(0.0, 0.0, z_far * z_near * inv, 0.0),
        )
    }

    /// Creates a left-handed perspective projection matrix with `[0, 1]` depth range.
    ///
    /// # Preconditions
    /// Expects `fov_y_radians`, `aspect_ratio` and `z_near` to be positive.
    /// Expects `z_far > z_near`.
    #[inline]
    #[must_use]
    pub fn perspective_lh(fov_y_radians: f32, aspect_ratio: f32, z_near: f32, z_far: f32) -> Self {
        let f: f32 = 1.0 / (fov_y_radians * 0.5).tan();
        let inv: f32 = 1.0 / (z_far - z_near);

        Self::from_cols(
            Vec4::new(f / aspect_ratio, 0.0, 0.0, 0.0),
            Vec4::new(0.0, f, 0.0, 0.0),
            Vec4::new(0.0, 0.0, z_far * inv, 1.0),
            Vec4::new(0.0, 0.0, -z_far * z_near * inv, 0.0),
        )
    }

    /// Creates a right-handed perspective projection matrix with `[-1, 1]` depth range.
    ///
    /// # Preconditions
    /// Expects `fov_y_radians`, `aspect_ratio` and `z_near` to be positive.
    /// Expects `z_far > z_near`.
    #[inline]
    #[must_use]
    pub fn perspective_rh_no(
        fov_y_radians: f32,
        aspect_ratio: f32,
        z_near: f32,
        z_far: f32,
    ) -> Self {
        let f: f32 = (fov_y_radians * 0.5).tan();
        let inv: f32 = 1.0 / (z_near - z_far);

        Self::from_cols(
            Vec4::new(f / aspect_ratio, 0.0, 0.0, 0.0),
            Vec4::new(0.0, f, 0.0, 0.0),
            Vec4::new(0.0, 0.0, (z_far + z_near) * inv, -1.0),
            Vec4::new(0.0, 0.0, (2.0 * z_far * z_near) * inv, 0.0),
        )
    }

    /// Creates a left-handed perspective projection matrix with `[-1, 1]` depth range.
    ///
    /// # Preconditions
    /// Expects `fov_y_radians`, `aspect_ratio` and `z_near` to be positive.
    /// Expects `z_far > z_near`.
    #[inline]
    #[must_use]
    pub fn perspective_lh_no(
        fov_y_radians: f32,
        aspect_ratio: f32,
        z_near: f32,
        z_far: f32,
    ) -> Self {
        let f: f32 = (fov_y_radians * 0.5).tan();
        let inv: f32 = 1.0 / (z_far - z_near);

        Self::from_cols(
            Vec4::new(f / aspect_ratio, 0.0, 0.0, 0.0),
            Vec4::new(0.0, f, 0.0, 0.0),
            Vec4::new(0.0, 0.0, (z_far + z_near) * inv, 1.0),
            Vec4::new(0.0, 0.0, -(2.0 * z_far * z_near) * inv, 0.0),
        )
    }

    /// Creates an infinite right-handed perspective projection matrix with `[0, 1]` depth range.
    ///
    /// # Preconditions
    /// Expects `fov_y_radians`, `aspect_ratio` and `z_near` to be positive.
    #[inline]
    #[must_use]
    pub fn perspective_infinite_rh(fov_y_radians: f32, aspect_ratio: f32, z_near: f32) -> Self {
        let f: f32 = 1.0 / (fov_y_radians * 0.5).tan();

        Self::from_cols(
            Vec4::new(f / aspect_ratio, 0.0, 0.0, 0.0),
            Vec4::new(0.0, f, 0.0, 0.0),
            Vec4::new(0.0, 0.0, -1.0, -1.0),
            Vec4::new(0.0, 0.0, -z_near, 0.0),
        )
    }

    /// Creates an infinite reverse right-handed perspective projection matrix with `[0, 1]` depth range.
    ///
    /// # Preconditions
    /// Expects `fov_y_radians`, `aspect_ratio` and `z_near` to be positive.
    #[inline]
    #[must_use]
    pub fn perspective_infinite_reverse_rh(
        fov_y_radians: f32,
        aspect_ratio: f32,
        z_near: f32,
    ) -> Self {
        let f: f32 = 1.0 / (fov_y_radians * 0.5).tan();

        Self::from_cols(
            Vec4::new(f / aspect_ratio, 0.0, 0.0, 0.0),
            Vec4::new(0.0, f, 0.0, 0.0),
            Vec4::new(0.0, 0.0, 0.0, -1.0),
            Vec4::new(0.0, 0.0, z_near, 0.0),
        )
    }

    /// Creates an infinite left-handed perspective projection matrix with `[0, 1]` depth range.
    ///
    /// # Preconditions
    /// Expects `fov_y_radians`, `aspect_ratio` and `z_near` to be positive.
    #[inline]
    #[must_use]
    pub fn perspective_infinite_lh(fov_y_radians: f32, aspect_ratio: f32, z_near: f32) -> Self {
        let f: f32 = 1.0 / (fov_y_radians * 0.5).tan();

        Self::from_cols(
            Vec4::new(f / aspect_ratio, 0.0, 0.0, 0.0),
            Vec4::new(0.0, f, 0.0, 0.0),
            Vec4::new(0.0, 0.0, 1.0, 1.0),
            Vec4::new(0.0, 0.0, -z_near, 0.0),
        )
    }

    /// Creates an infinite reverse left-handed perspective projection matrix with `[0, 1]` depth range.
    ///
    /// # Preconditions
    /// Expects `fov_y_radians`, `aspect_ratio` and `z_near` to be positive.
    #[inline]
    #[must_use]
    pub fn perspective_infinite_reverse_lh(
        fov_y_radians: f32,
        aspect_ratio: f32,
        z_near: f32,
    ) -> Self {
        let f: f32 = 1.0 / (fov_y_radians * 0.5).tan();

        Self::from_cols(
            Vec4::new(f / aspect_ratio, 0.0, 0.0, 0.0),
            Vec4::new(0.0, f, 0.0, 0.0),
            Vec4::new(0.0, 0.0, 0.0, 1.0),
            Vec4::new(0.0, 0.0, z_near, 0.0),
        )
    }

    /// Creates a right-handed orthographic projection matrix with `[0, 1]` depth range.
    ///
    /// # Preconditions
    /// 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 rl: f32 = right - left;
        let tb: f32 = top - bottom;
        let fn_: f32 = z_far - z_near;

        Self::from_cols(
            Vec4::new(2.0 / rl, 0.0, 0.0, 0.0),
            Vec4::new(0.0, 2.0 / tb, 0.0, 0.0),
            Vec4::new(0.0, 0.0, -1.0 / fn_, 0.0),
            Vec4::new(
                -(right + left) / rl,
                -(top + bottom) / tb,
                -z_near / fn_,
                1.0,
            ),
        )
    }

    /// Creates a left-handed orthographic projection matrix with `[0, 1]` depth range.
    ///
    /// # Preconditions
    /// 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 rl: f32 = right - left;
        let tb: f32 = top - bottom;
        let fn_: f32 = z_far - z_near;

        Self::from_cols(
            Vec4::new(2.0 / rl, 0.0, 0.0, 0.0),
            Vec4::new(0.0, 2.0 / tb, 0.0, 0.0),
            Vec4::new(0.0, 0.0, 1.0 / fn_, 0.0),
            Vec4::new(
                -(right + left) / rl,
                -(top + bottom) / tb,
                -z_near / fn_,
                1.0,
            ),
        )
    }

    /// Creates a right-handed orthographic projection matrix with `[-1, 1]` depth range.
    ///
    /// # Preconditions
    /// Expects `left != right`, `bottom != top` and `z_near != z_far`.
    #[inline]
    #[must_use]
    pub const fn orthographic_rh_no(
        left: f32,
        right: f32,
        bottom: f32,
        top: f32,
        z_near: f32,
        z_far: f32,
    ) -> Self {
        let rl: f32 = right - left;
        let tb: f32 = top - bottom;
        let fn_: f32 = z_far - z_near;

        Self::from_cols(
            Vec4::new(2.0 / rl, 0.0, 0.0, 0.0),
            Vec4::new(0.0, 2.0 / tb, 0.0, 0.0),
            Vec4::new(0.0, 0.0, -2.0 / fn_, 0.0),
            Vec4::new(
                -(right + left) / rl,
                -(top + bottom) / tb,
                -(z_far + z_near) / fn_,
                1.0,
            ),
        )
    }

    /// Creates a left-handed orthographic projection matrix with `[-1, 1]` depth range.
    ///
    /// # Preconditions
    /// Expects `left != right`, `bottom != top` and `z_near != z_far`.
    #[inline]
    #[must_use]
    pub const fn orthographic_lh_no(
        left: f32,
        right: f32,
        bottom: f32,
        top: f32,
        z_near: f32,
        z_far: f32,
    ) -> Self {
        let rl: f32 = right - left;
        let tb: f32 = top - bottom;
        let fn_: f32 = z_far - z_near;

        Self::from_cols(
            Vec4::new(2.0 / rl, 0.0, 0.0, 0.0),
            Vec4::new(0.0, 2.0 / tb, 0.0, 0.0),
            Vec4::new(0.0, 0.0, 2.0 / fn_, 0.0),
            Vec4::new(
                -(right + left) / rl,
                -(top + bottom) / tb,
                -(z_far + z_near) / fn_,
                1.0,
            ),
        )
    }

    /// Returns the component-wise sum of `self` and `rhs`.
    #[inline]
    #[must_use]
    pub const fn add(self, rhs: Self) -> Self {
        Self::from_cols(
            self.col0.add(rhs.col0),
            self.col1.add(rhs.col1),
            self.col2.add(rhs.col2),
            self.col3.add(rhs.col3),
        )
    }

    /// Sets self to the component-wise sum of `self` and `rhs`.
    #[inline]
    pub const fn add_assign(&mut self, rhs: Self) {
        self.col0.add_assign(rhs.col0);
        self.col1.add_assign(rhs.col1);
        self.col2.add_assign(rhs.col2);
        self.col3.add_assign(rhs.col3);
    }

    /// Returns the component-wise difference of `self` and `rhs`.
    #[inline]
    #[must_use]
    pub const fn sub(self, rhs: Self) -> Self {
        Self::from_cols(
            self.col0.sub(rhs.col0),
            self.col1.sub(rhs.col1),
            self.col2.sub(rhs.col2),
            self.col3.sub(rhs.col3),
        )
    }

    /// Sets self to the component-wise difference of `self` and `rhs`.
    #[inline]
    pub const fn sub_assign(&mut self, rhs: Self) {
        self.col0.sub_assign(rhs.col0);
        self.col1.sub_assign(rhs.col1);
        self.col2.sub_assign(rhs.col2);
        self.col3.sub_assign(rhs.col3);
    }

    /// Returns the matrix product of `self` and `rhs`.
    #[inline]
    #[must_use]
    pub const fn mul(self, rhs: Self) -> Self {
        Self::from_cols(
            self.mul_vec4(rhs.col0),
            self.mul_vec4(rhs.col1),
            self.mul_vec4(rhs.col2),
            self.mul_vec4(rhs.col3),
        )
    }

    /// Sets `self` to the matrix product of `self` and `rhs`.
    #[inline]
    pub const fn mul_assign(&mut self, rhs: Self) {
        *self = self.mul(rhs)
    }

    /// Returns `rhs` transformed by `self`.
    #[inline]
    #[must_use]
    pub const fn mul_vec4(self, rhs: Vec4) -> Vec4 {
        self.col0
            .mul_scalar(rhs.x)
            .add(self.col1.mul_scalar(rhs.y))
            .add(self.col2.mul_scalar(rhs.z))
            .add(self.col3.mul_scalar(rhs.w))
    }

    /// Returns `self` with each component multiplied by `rhs`.
    #[inline]
    #[must_use]
    pub const fn mul_scalar(self, rhs: f32) -> Self {
        Self::from_cols(
            self.col0.mul_scalar(rhs),
            self.col1.mul_scalar(rhs),
            self.col2.mul_scalar(rhs),
            self.col3.mul_scalar(rhs),
        )
    }

    /// Sets `self` to `self` with each component multiplied by `rhs`.
    #[inline]
    pub const fn mul_assign_scalar(&mut self, rhs: f32) {
        *self = self.mul_scalar(rhs);
    }

    /// Returns `self` with each component divided by `rhs`.
    #[inline]
    #[must_use]
    pub const fn div_scalar(self, rhs: f32) -> Self {
        Self::from_cols(
            self.col0.div_scalar(rhs),
            self.col1.div_scalar(rhs),
            self.col2.div_scalar(rhs),
            self.col3.div_scalar(rhs),
        )
    }

    /// Sets `self` to `self` with each component divided by `rhs`.
    #[inline]
    pub const fn div_assign_scalar(&mut self, rhs: f32) {
        *self = self.div_scalar(rhs);
    }

    /// Returns `self` with negated components.
    #[inline]
    #[must_use]
    pub const fn neg(self) -> Self {
        Self::from_cols(
            self.col0.neg(),
            self.col1.neg(),
            self.col2.neg(),
            self.col3.neg(),
        )
    }

    /// Returns a reference to the components as `[col0, col1, col2, col3]`,
    /// with each column as `[x, y, z, w]`.
    ///
    /// # Safety
    /// Relies on `Mat4` being `#[repr(C)]`
    #[inline]
    #[must_use]
    pub const fn as_ref(&self) -> &[f32; 16] {
        unsafe { &*(self as *const Mat4 as *const [f32; 16]) }
    }

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

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

    /// Returns column of `self` at `index`.
    ///
    /// `0` maps to `col0`, `1` maps to `col1`, `2` maps to `col2`, `3` maps to `col3`.
    ///
    /// # Panics
    /// Panics if `index` is greater than `3`.
    #[inline]
    pub const fn col(self, index: usize) -> Vec4 {
        match index {
            0 => self.col0,
            1 => self.col1,
            2 => self.col2,
            3 => self.col3,
            _ => panic!("index out of bounds"),
        }
    }

    /// Returns row 0 of `self`.
    #[inline]
    #[must_use]
    pub const fn row0(self) -> Vec4 {
        Vec4::new(self.col0.x, self.col1.x, self.col2.x, self.col3.x)
    }

    /// Returns row 1 of `self`.
    #[inline]
    #[must_use]
    pub const fn row1(self) -> Vec4 {
        Vec4::new(self.col0.y, self.col1.y, self.col2.y, self.col3.y)
    }

    /// Returns row 2 of `self`.
    #[inline]
    #[must_use]
    pub const fn row2(self) -> Vec4 {
        Vec4::new(self.col0.z, self.col1.z, self.col2.z, self.col3.z)
    }

    /// Returns row 3 of `self`.
    #[inline]
    #[must_use]
    pub const fn row3(self) -> Vec4 {
        Vec4::new(self.col0.w, self.col1.w, self.col2.w, self.col3.w)
    }

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

    /// Returns row of `self` at `index`.
    ///
    /// `0` maps to `row0`, `1` maps to `row1`, `2` maps to `row2`, `3` maps to `row3`.
    ///
    /// # Panics
    /// Panics if `index` is greater than `3`.
    #[inline]
    #[must_use]
    pub const fn row(self, index: usize) -> Vec4 {
        match index {
            0 => self.row0(),
            1 => self.row1(),
            2 => self.row2(),
            3 => self.row3(),
            _ => panic!("index out of bounds"),
        }
    }

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

    /// Returns the transpose of `self`.
    #[inline]
    #[must_use]
    pub const fn transpose(self) -> Self {
        Self::from_cols(self.row0(), self.row1(), self.row2(), self.row3())
    }

    /// Returns the determinant of `self`.
    #[must_use]
    pub const fn determinant(self) -> f32 {
        let a: Vec4 = self.col0;
        let b: Vec4 = self.col1;
        let c: Vec4 = self.col2;
        let d: Vec4 = self.col3;

        // 2x2 minors from columns a/b.
        let s01: f32 = a.x * b.y - a.y * b.x;
        let s02: f32 = a.x * b.z - a.z * b.x;
        let s03: f32 = a.x * b.w - a.w * b.x;
        let s12: f32 = a.y * b.z - a.z * b.y;
        let s13: f32 = a.y * b.w - a.w * b.y;
        let s23: f32 = a.z * b.w - a.w * b.z;

        // 2x2 minors from columns c/d.
        let t01: f32 = c.x * d.y - c.y * d.x;
        let t02: f32 = c.x * d.z - c.z * d.x;
        let t03: f32 = c.x * d.w - c.w * d.x;
        let t12: f32 = c.y * d.z - c.z * d.y;
        let t13: f32 = c.y * d.w - c.w * d.y;
        let t23: f32 = c.z * d.w - c.w * d.z;

        s01 * t23 - s02 * t13 + s03 * t12 + s12 * t03 - s13 * t02 + s23 * t01
    }

    /// Returns the inverse of `self` or `None` if `self` is not invertible.
    ///
    /// # Edge cases
    /// If `self` is not invertible, returns `None`.
    #[must_use]
    pub const fn try_inverse(self) -> Option<Self> {
        let a: Vec4 = self.col0;
        let b: Vec4 = self.col1;
        let c: Vec4 = self.col2;
        let d: Vec4 = self.col3;

        let s01: f32 = a.x * b.y - a.y * b.x;
        let s02: f32 = a.x * b.z - a.z * b.x;
        let s03: f32 = a.x * b.w - a.w * b.x;
        let s12: f32 = a.y * b.z - a.z * b.y;
        let s13: f32 = a.y * b.w - a.w * b.y;
        let s23: f32 = a.z * b.w - a.w * b.z;

        let t01: f32 = c.x * d.y - c.y * d.x;
        let t02: f32 = c.x * d.z - c.z * d.x;
        let t03: f32 = c.x * d.w - c.w * d.x;
        let t12: f32 = c.y * d.z - c.z * d.y;
        let t13: f32 = c.y * d.w - c.w * d.y;
        let t23: f32 = c.z * d.w - c.w * d.z;

        let det: f32 = s01 * t23 - s02 * t13 + s03 * t12 + s12 * t03 - s13 * t02 + s23 * t01;

        if det.abs() < Self::EPSILON {
            return None;
        }
        let det_recip: f32 = det.recip();

        let col0: Vec4 = Vec4::new(
            (b.y * t23 - b.z * t13 + b.w * t12) * det_recip,
            (-a.y * t23 + a.z * t13 - a.w * t12) * det_recip,
            (d.y * s23 - d.z * s13 + d.w * s12) * det_recip,
            (-c.y * s23 + c.z * s13 - c.w * s12) * det_recip,
        );

        let col1: Vec4 = Vec4::new(
            (-b.x * t23 + b.z * t03 - b.w * t02) * det_recip,
            (a.x * t23 - a.z * t03 + a.w * t02) * det_recip,
            (-d.x * s23 + d.z * s03 - d.w * s02) * det_recip,
            (c.x * s23 - c.z * s03 + c.w * s02) * det_recip,
        );

        let col2: Vec4 = Vec4::new(
            (b.x * t13 - b.y * t03 + b.w * t01) * det_recip,
            (-a.x * t13 + a.y * t03 - a.w * t01) * det_recip,
            (d.x * s13 - d.y * s03 + d.w * s01) * det_recip,
            (-c.x * s13 + c.y * s03 - c.w * s01) * det_recip,
        );

        let col3: Vec4 = Vec4::new(
            (-b.x * t12 + b.y * t02 - b.z * t01) * det_recip,
            (a.x * t12 - a.y * t02 + a.z * t01) * det_recip,
            (-d.x * s12 + d.y * s02 - d.z * s01) * det_recip,
            (c.x * s12 - c.y * s02 + c.z * s01) * det_recip,
        );

        Some(Mat4::from_cols(col0, col1, col2, col3))
    }

    /// Returns the inverse of `self`.
    ///
    /// # Panics
    /// Panics if `self` is not invertible.
    #[must_use]
    pub const fn inverse(self) -> Self {
        self.try_inverse().expect("matrix not invertible")
    }

    /// Returns `true` if `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
    }

    /// Returns `v` transformed as a point by `self`.
    ///
    /// Applies the full matrix including translation.
    #[inline]
    #[must_use]
    pub const fn transform_point(self, v: Vec3) -> Vec3 {
        Vec3::new(
            self.col0.x * v.x + self.col1.x * v.y + self.col2.x * v.z + self.col3.x,
            self.col0.y * v.x + self.col1.y * v.y + self.col2.y * v.z + self.col3.y,
            self.col0.z * v.x + self.col1.z * v.y + self.col2.z * v.z + self.col3.z,
        )
    }

    /// Returns `v` transformed as a direction by `self`.
    ///
    /// Applies only the linear part of `self`, ignoring translation.
    /// Use this for directions, normals, and offsets.
    #[inline]
    #[must_use]
    pub const fn transform_vector(self, v: Vec3) -> Vec3 {
        Vec3::new(
            self.col0.x * v.x + self.col1.x * v.y + self.col2.x * v.z,
            self.col0.y * v.x + self.col1.y * v.y + self.col2.y * v.z,
            self.col0.z * v.x + self.col1.z * v.y + self.col2.z * v.z,
        )
    }

    /// Returns `self` with the absolute value of each component.
    #[inline]
    #[must_use]
    pub const fn abs(self) -> Self {
        Self::from_cols(
            self.col0.abs(),
            self.col1.abs(),
            self.col2.abs(),
            self.col3.abs(),
        )
    }

    /// Returns `self` with the reciprocal value of each component.
    #[inline]
    #[must_use]
    pub const fn recip(self) -> Self {
        Self::from_cols(
            self.col0.recip(),
            self.col1.recip(),
            self.col2.recip(),
            self.col3.recip(),
        )
    }

    /// Returns `true` if `self` and `other` approximately represent the same matrix.
    ///
    /// Uses [`EPSILON`](Self::EPSILON) as the threshold.
    #[inline]
    #[must_use]
    pub const fn approx_eq(self, other: Self) -> bool {
        let a: [f32; 16] = self.to_array();
        let b: [f32; 16] = other.to_array();

        let mut i: usize = 0;
        while i < 16 {
            if (a[i] - b[i]).abs() >= Self::EPSILON {
                return false;
            }
            i += 1;
        }
        true
    }

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

        let mut i: usize = 0;
        while i < 16 {
            if !m[i].is_finite() {
                return false;
            }
            i += 1;
        }
        true
    }

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

        let mut i: usize = 0;
        while i < 16 {
            if m[i].is_nan() {
                return true;
            }
            i += 1;
        }
        false
    }
}

impl Add for Mat4 {
    type Output = Self;

    /// Returns the component-wise sum of `self` and `rhs`.
    ///
    /// Equivalent to [`add`](Self::add()).
    #[inline]
    fn add(self, rhs: Self) -> Self::Output {
        self.add(rhs)
    }
}

impl AddAssign for Mat4 {
    /// Sets self to the component-wise sum of `self` and `rhs`.
    ///
    /// Equivalent to [`add_assign`](Self::add_assign()).
    #[inline]
    fn add_assign(&mut self, rhs: Self) {
        self.add_assign(rhs);
    }
}

impl Sub for Mat4 {
    type Output = Self;

    /// Returns the component-wise difference of `self` and `rhs`.
    ///
    /// Equivalent to [`sub`](Self::sub()).
    #[inline]
    fn sub(self, rhs: Self) -> Self::Output {
        self.sub(rhs)
    }
}

impl SubAssign for Mat4 {
    /// Sets self to the component-wise difference of `self` and `rhs`.
    ///
    /// Equivalent to [`sub_assign`](Self::sub_assign()).
    #[inline]
    fn sub_assign(&mut self, rhs: Self) {
        self.sub_assign(rhs);
    }
}

impl Mul for Mat4 {
    type Output = Self;

    /// Returns the matrix product of `self` and `rhs`.
    ///
    /// Equivalent to [`mul`](Self::mul()).
    #[inline]
    fn mul(self, rhs: Self) -> Self::Output {
        self.mul(rhs)
    }
}

impl MulAssign for Mat4 {
    /// Sets `self` to the matrix product of `self` and `rhs`.
    ///
    /// Equivalent to [`mul_assign`](Self::mul_assign()).
    #[inline]
    fn mul_assign(&mut self, rhs: Self) {
        self.mul_assign(rhs);
    }
}

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

    /// Returns `rhs` transformed by `self`.
    ///
    /// Equivalent to [`mul_vec4`](Self::mul_vec4()).
    #[inline]
    fn mul(self, rhs: Vec4) -> Self::Output {
        self.mul_vec4(rhs)
    }
}

impl Mul<f32> for Mat4 {
    type Output = Self;

    /// Returns `self` with each component multiplied by `rhs`.
    ///
    /// Equivalent to [`mul_scalar`](Self::mul_scalar()).
    #[inline]
    fn mul(self, rhs: f32) -> Self::Output {
        self.mul_scalar(rhs)
    }
}

impl MulAssign<f32> for Mat4 {
    /// Sets `self` to `self` with each component multiplied by `rhs`.
    ///
    /// Equivalent to [`mul_assign_scalar`](Self::mul_assign_scalar()).
    #[inline]
    fn mul_assign(&mut self, rhs: f32) {
        self.mul_assign_scalar(rhs);
    }
}

impl Div<f32> for Mat4 {
    type Output = Self;

    /// Returns `self` with each component divided by `rhs`.
    ///
    /// Equivalent to [`div_scalar`](Self::div_scalar()).
    #[inline]
    fn div(self, rhs: f32) -> Self::Output {
        self.div_scalar(rhs)
    }
}

impl DivAssign<f32> for Mat4 {
    /// Sets `self` to `self` with each component divided by `rhs`.
    ///
    /// Equivalent to [`div_assign_scalar`](Self::div_assign_scalar()).
    #[inline]
    fn div_assign(&mut self, rhs: f32) {
        self.div_assign_scalar(rhs);
    }
}

impl Neg for Mat4 {
    type Output = Self;

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

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

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

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

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

impl From<[f32; 16]> for Mat4 {
    /// Creates a `Mat4` from a `[f32; 16]`, in column-major order:
    /// `[col0.x, col0.y, col0.z, col0.w, col1.x, ..., col3.w]`.
    ///
    /// Equivalent to [`from_array`](Self::from_array()).
    #[inline]
    fn from(v: [f32; 16]) -> Self {
        Self::from_array(v)
    }
}

impl From<Mat4> for [f32; 16] {
    /// Returns a copy of the components as `[col0, col1, col2, col3]`,
    /// with each column as `[x, y, z, w]`.
    ///
    /// Equivalent to [`to_array`](Mat4::to_array()).
    #[inline]
    fn from(v: Mat4) -> Self {
        v.to_array()
    }
}

impl From<Mat3> for Mat4 {
    /// Creates a `Mat4` from a `Mat3`, setting `col3` and `row3` to identity.
    ///
    /// Equivalent to [`from_mat3`](Self::from_mat3()).
    #[inline]
    fn from(v: Mat3) -> Self {
        Self::from_mat3(v)
    }
}

impl From<Affine3> for Mat4 {
    /// Returns the equivalent `Mat4` of `v`.
    ///
    /// Equivalent to [`to_mat4`](Affine3::to_mat4()).
    #[inline]
    fn from(v: Affine3) -> Self {
        v.to_mat4()
    }
}

#[cfg(feature = "serde")]
impl serde::Serialize for Mat4 {
    /// Serializes as a `[f32; 16]` in column-major order.
    ///
    /// # Format
    /// Encoded as a flat 16-element array. Component order matches [`to_array`](Self::to_array()).
    fn serialize<S: serde::Serializer>(&self, serializer: S) -> Result<S::Ok, S::Error> {
        self.to_array().serialize(serializer)
    }
}

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

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

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