omelet 0.1.2

use rand::Rng;

use crate::mat4::Mat4;
use crate::utils;
use crate::utils::epsilon_eq;
use crate::utils::epsilon_eq_default;
use crate::utils::is_near_zero_default;
use crate::vec::Vec3;

use core::fmt;
use std::cmp::PartialEq;
use std::ops::{
    Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign,
};

/// A 4D vector with x, y, z, and w components.
#[derive(Debug, Clone, Copy)]
pub struct Vec4 {
    pub x: f32,
    pub y: f32,
    pub z: f32,
    pub w: f32,
}

impl Vec4 {
    // ============= Construction and Conversion =============

    /// Creates a 4D vector with the given `x`, `y`, `z`, and `w` components.
    ///
    /// # Parameters
    /// - `x`: The x component of the vector.
    /// - `y`: The y component of the vector.
    /// - `z`: The z component of the vector.
    /// - `w`: The w component of the vector.
    ///
    /// # Returns
    /// A new `Vec4` instance with the specified components.
    pub const fn new(x: f32, y: f32, z: f32, w: f32) -> Vec4 {
        Vec4 { x, y, z, w }
    }

    /// Converts the vector into an array of 4 `f32` components `[x, y, z, w]`.
    ///
    /// # Returns
    /// An array [f32; 4].
    pub fn to_array(self) -> [f32; 4] {
        [self.x, self.y, self.z, self.w]
    }

    /// Converts an array into a vector.
    ///
    /// # Parameters
    /// - `arr`: The array to convert.
    ///
    /// # Returns
    /// `Vec4` which is equivalent to `(arr[0], arr[1], arr[2], arr[3])`.
    pub fn from_array(arr: [f32; 4]) -> Vec4 {
        Vec4::new(arr[0], arr[1], arr[2], arr[3])
    }

    /// Converts the vector into a tuple `(f32, f32, f32, f32)`.
    ///
    /// # Returns
    /// `(f32, f32, f32, f32)` which is equivalent to (self.x, self.y, self.z, self.w).
    pub fn to_tuple(self) -> (f32, f32, f32, f32) {
        (self.x, self.y, self.z, self.w)
    }

    /// Converts a tuple to a vector.
    ///
    /// # Parameters
    /// - `t`: Tuple `(f32, f32, f32, f32)` to be converted.
    ///
    /// # Returns
    /// `Vec4` which is equivalent to `(t.0, t.1, t.2, t.3)`.
    pub fn from_tuple(tuple: (f32, f32, f32, f32)) -> Vec4 {
        Vec4::new(tuple.0, tuple.1, tuple.2, tuple.3)
    }

    /// Method to create a Vec4 from a Vec3.
    ///
    /// # Parameters
    /// - `v`: Vec2 to use for `x`, `y`, and `z` components.
    /// - `w`: f32 to use for the `w` component.
    ///
    /// # Returns
    /// `Vec4`.
    pub fn from_vec3_w(v: Vec3, w: f32) -> Vec4 {
        Vec4::new(v.x, v.y, v.z, w)
    }

    /// Returns a Vec3 that is the `x`, `y`, and `z` components of `self`.
    ///
    /// # Returns
    /// `Vec3`.
    pub fn xyz(self) -> Vec3 {
        Vec3::new(self.x, self.y, self.z)
    }

    // ============= Constants ==============
    /// Returns a `Vec4` where all components are zero.
    pub const ZERO: Self = Self {
        x: 0.0,
        y: 0.0,
        z: 0.0,
        w: 0.0,
    };

    /// Returns a `Vec4` where all components are NaN.
    pub const NAN: Self = Self {
        x: f32::NAN,
        y: f32::NAN,
        z: f32::NAN,
        w: f32::NAN,
    };

    /// Returns a `Vec4` where all components are infinity.
    pub const INFINITY: Self = Self {
        x: f32::INFINITY,
        y: f32::INFINITY,
        z: f32::INFINITY,
        w: f32::INFINITY,
    };

    /// Returns a `Vec4` equivalent to `(1.0, 0.0, 0.0, 0.0)`.
    pub const X: Self = Self {
        x: 1.0,
        y: 0.0,
        z: 0.0,
        w: 0.0,
    };

    /// Returns a `Vec4` equivalent to `(0.0, 1.0, 0.0, 0.0)`.
    pub const Y: Self = Self {
        x: 0.0,
        y: 1.0,
        z: 0.0,
        w: 0.0,
    };

    /// Returns a `Vec4` equivalent to `(0.0, 0.0, 1.0, 0.0)`.
    pub const Z: Self = Self {
        x: 0.0,
        y: 0.0,
        z: 1.0,
        w: 0.0,
    };

    /// Returns a `Vec4` equivalent to `(0.0, 0.0, 0.0, 1.0)`.
    pub const W: Self = Self {
        x: 0.0,
        y: 0.0,
        z: 0.0,
        w: 1.0,
    };

    // ============= Math Utilities =============

    /// Returns a new vector with each component replaced by its absolute value.
    ///
    /// This is useful for obtaining the magnitude of each axis regardless of sign.
    ///
    /// # Returns
    /// A `Vec4` where `x`, `y`, `z`, and `w` are the absolute vaklues of the original vectors components.
    pub fn abs(self) -> Vec4 {
        Vec4::new(self.x.abs(), self.y.abs(), self.z.abs(), self.w.abs())
    }

    /// Returns a new vector where each component is replaced by its sign.
    ///
    /// The sign of each component can be `-1.0`, `0.0`, or `1.0` depending on whether the
    /// component is negative, zero, or positive respectively.
    ///
    /// # Returns
    /// A `Vec4` representing the sign of each component.
    pub fn signum(self) -> Self {
        Vec4::new(
            if self.x > 0.0 {
                1.0
            } else if self.x < 0.0 {
                -1.0
            } else {
                0.0
            },
            if self.y > 0.0 {
                1.0
            } else if self.y < 0.0 {
                -1.0
            } else {
                0.0
            },
            if self.z > 0.0 {
                1.0
            } else if self.z < 0.0 {
                -1.0
            } else {
                0.0
            },
            if self.w > 0.0 {
                1.0
            } else if self.w < 0.0 {
                -1.0
            } else {
                0.0
            },
        )
    }

    /// Returns a new vector where each component is replaced by its IEEE 754 signum.
    ///
    /// Unlike [`signum()`](#method.signum), this method treats zero as positive:
    /// `0.0.signum()` is `1.0`, and `-0.0.signum()` is `-1.0`.
    ///
    /// # Returns
    /// A `Vec4` where each component is `-1.0` if negative, `1.0` otherwise.
    pub fn ieee_signum(self) -> Self {
        Vec4::new(
            self.x.signum(),
            self.y.signum(),
            self.z.signum(),
            self.w.signum(),
        )
    }

    /// Clamps each component of the vector between the specified `min` and `max` values.
    ///
    /// # Parameters
    /// - `min`: The minimum allowed value for each component.
    /// - `max`: The maximum allowed value for each component.
    ///
    /// # Returns
    /// A new `Vec4` where each component is limited to the range `[min, max]`.
    pub fn clamp(self, min: f32, max: f32) -> Vec4 {
        Vec4::new(
            utils::clamp(self.x, min, max),
            utils::clamp(self.y, min, max),
            utils::clamp(self.z, min, max),
            utils::clamp(self.w, min, max),
        )
    }

    /// Returns the component-wise minimum between `self` and another vector.
    ///
    /// # Parammeters
    /// - `other`: The other vector.
    ///
    /// # Returns
    /// A `Vec4` with each component being the minimum of the corresponding components.
    ///
    /// ```rust
    /// use omelet::vec4::Vec4;
    /// let a = Vec4::new(1.0, 4.0, 2.0, 1.0);
    /// let b = Vec4::new(3.0, 1.0, 1.0, 5.0);
    /// assert_eq!(a.min(b), Vec4::new(1.0, 1.0, 1.0, 1.0));
    /// ```
    pub fn min(self, other: Vec4) -> Vec4 {
        Vec4::new(
            self.x.min(other.x),
            self.y.min(other.y),
            self.z.min(other.z),
            self.w.min(other.w),
        )
    }

    /// Returns the component-wise maximum between `self` and another vector.
    ///
    /// # Parammeters
    /// - `other`: The other vector.
    ///
    /// # Returns
    /// A `Vec4` with each component being the maximum of the corresponding components.
    ///
    /// ```rust
    /// use omelet::vec4::Vec4;
    /// let a = Vec4::new(1.0, 4.0, 2.0, 1.0);
    /// let b = Vec4::new(3.0, 1.0, 1.0, 5.0);
    /// assert_eq!(a.max(b), Vec4::new(3.0, 4.0, 2.0, 5.0));
    /// ```
    pub fn max(self, other: Vec4) -> Vec4 {
        Vec4::new(
            self.x.max(other.x),
            self.y.max(other.y),
            self.z.max(other.z),
            self.w.max(other.w),
        )
    }

    /// Computes the **scalar 4D triple product** for four vectors:
    /// `a · (b × c × d)` - interpreted as the **determinant of the 4x4 matrix**
    /// whose columns are the vectors `a`, `b`, `c`, and `d`.
    ///
    /// This is used to compute **signed 3D hypervolume** (also called the 4-parallelotope volume).
    ///
    /// # Parameters
    /// - `a`, `b`, `c`, `d`: 4D vectors used to form a 4x4 matrix.
    ///
    /// # Returns
    /// Signed scalar volume (can be negative depending on handedness).
    pub fn triple_product_4d(a: Vec4, b: Vec4, c: Vec4, d: Vec4) -> f32 {
        let m = Mat4::new(
            Vec4::new(a.x, b.x, c.x, d.x),
            Vec4::new(a.y, b.y, c.y, d.y),
            Vec4::new(a.z, b.z, c.z, d.z),
            Vec4::new(a.w, b.w, c.w, d.w),
        );

        // Compute determinant using Laplace expansion or precomputed formula
        // This expansion is hardcoded for performance
        m[0][0]
            * (m[1][1] * (m[2][2] * m[3][3] - m[2][3] * m[3][2])
                - m[1][2] * (m[2][1] * m[3][3] - m[2][3] * m[3][1])
                + m[1][3] * (m[2][1] * m[3][2] - m[2][2] * m[3][1]))
            - m[0][1]
                * (m[1][0] * (m[2][2] * m[3][3] - m[2][3] * m[3][2])
                    - m[1][2] * (m[2][0] * m[3][3] - m[2][3] * m[3][0])
                    + m[1][3] * (m[2][0] * m[3][2] - m[2][2] * m[3][0]))
            + m[0][2]
                * (m[1][0] * (m[2][1] * m[3][3] - m[2][3] * m[3][1])
                    - m[1][1] * (m[2][0] * m[3][3] - m[2][3] * m[3][0])
                    + m[1][3] * (m[2][0] * m[3][1] - m[2][1] * m[3][0]))
            - m[0][3]
                * (m[1][0] * (m[2][1] * m[3][2] - m[2][2] * m[3][1])
                    - m[1][1] * (m[2][0] * m[3][2] - m[2][2] * m[3][0])
                    + m[1][2] * (m[2][0] * m[3][1] - m[2][1] * m[3][0]))
    }

    /// Returns the **unsigned 4D hypervolume** formed by the vectors `a`, `b`, `c`, and `d`,
    /// using the **absolute value** of the 4D scalar triple product.
    ///
    /// # parameters
    /// - `a`, `b`, `c`, `d`: 4D vectors forming the 4D volume.
    ///
    /// # Returns
    /// The absolute hypervolume (alsways non-negative).
    pub fn hypervolume_4d(a: Vec4, b: Vec4, c: Vec4, d: Vec4) -> f32 {
        Vec4::triple_product_4d(a, b, c, d).abs()
    }

    /// Returns a **unit-length vector perpendicular** to the current vector in 4D space.
    ///
    /// Chooses a perpendicular direction by nulling the smallest component and rotating others.
    /// If the vector is degenerate (zero vector), a defauly fallback of (1, 0, 0, 0) is returned.
    ///
    /// # Returns
    /// A normalized perpendicular Vec4.
    pub fn perpendicular(self) -> Vec4 {
        // Pick a perpendicular vector by zeroing the smallest component
        let mut v = if self.x.abs() <= self.y.abs()
            && self.x.abs() <= self.z.abs()
            && self.x.abs() <= self.w.abs()
        {
            Vec4::new(0.0, -self.z, self.y, -self.w)
        } else if self.y.abs() <= self.z.abs() && self.y.abs() <= self.w.abs() {
            Vec4::new(-self.z, 0.0, self.x, -self.w)
        } else if self.z.abs() <= self.w.abs() {
            Vec4::new(-self.y, self.x, 0.0, -self.w)
        } else {
            Vec4::new(-self.y, self.x, -self.z, 0.0)
        };
        v = v.normalize_or_zero();
        if v.is_zero() {
            Vec4::new(1.0, 0.0, 0.0, 0.0) // Fallback
        } else {
            v
        }
    }

    // ============= Magnitude and Normalization =============

    /// Calculates and returns the magnitude (length) of the vector.
    ///
    /// The length is computed as the squre root of the sum of sqaures of `x`, `y`, `z`, and `w`.
    ///
    /// # Returns
    /// A `f32` representing the Euclidean length of the vector.
    ///
    /// Note: This calculation uses sqrt(), which can be taxing on the system.
    pub fn length(&self) -> f32 {
        (self.x * self.x + self.y * self.y + self.z * self.z + self.w * self.w).sqrt()
    }

    /// Computes and returns the squared magnitude (length) of the vector.
    ///
    /// This avoids the expensive square root calculation used in `length()`,
    /// and it useful when only relative length or comparisons are needed.
    ///
    /// # Returns
    /// A `f32` representing the squared length `(x² + y² + z² + w²)`.
    pub fn squared_length(&self) -> f32 {
        self.x * self.x + self.y * self.y + self.z * self.z + self.w * self.w
    }

    /// Returns a normalized (unit length) vector in the same direction as `self`.
    ///
    /// If the length is near 0 within the default tolerance, then it will return a zero vector.
    ///
    /// # Returns
    /// `Vec4` either normalized or zero.
    pub fn normalize_or_zero(&self) -> Vec4 {
        let len = self.length();
        if len == 0.0 {
            Vec4::new(0.0, 0.0, 0.0, 0.0)
        } else {
            Vec4::new(self.x / len, self.y / len, self.z / len, self.w / len)
        }
    }

    /// Attempts to return a normalized version of the vector, or `None` if the length
    /// is zero or nearly zero.
    ///
    /// # Returns
    /// - `Some(Vec4)` containing the normalized vector if length > epsilon.
    /// - `None` if vector length is zero or near zero (within `1e-6`).
    pub fn try_normalize(&self) -> Option<Vec4> {
        let len = self.length();
        if is_near_zero_default(len) {
            None
        } else {
            Some(Vec4::new(
                self.x / len,
                self.y / len,
                self.z / len,
                self.w / len,
            ))
        }
    }

    /// Returns a normalized (unit length) vector in the same direction as `self`.
    ///
    /// # Panics
    /// This function panics if the vector length is zero to prevent division by zero.
    ///
    /// # Returns
    /// A `Vec4` with length 1 pointing in the same direction.
    pub fn normalize(&self) -> Vec4 {
        let len = self.length();
        assert!(
            !is_near_zero_default(len),
            "Cannot normalize zero length vector"
        );
        Vec4::new(self.x / len, self.y / len, self.z / len, self.w / len)
    }

    /// Returns `true` if the vector is normalized within the small epsilon tolerance.
    ///
    /// This checks if the length is approx 1.0 within 1e-6.
    ///
    /// # Returns
    /// Boolean indicating normalization status.
    pub fn is_normalized(self) -> bool {
        epsilon_eq(self.length(), 1.0, 1e-6)
    }

    /// Faster approx check if the squared length is approx 1 within epsilon.
    ///
    /// Avoids computing a square root, useful for performance-critical code.
    ///
    /// # Returns
    /// Boolean indicating if vector is approx normalized.
    pub fn is_normalized_fast(self) -> bool {
        epsilon_eq(self.squared_length(), 1.0, 1e-6)
    }

    // ============= Dot, Cross, and Angles =============

    /// Computes the dot product between `self` and another vector.
    ///
    /// The dot product is defined as `x1 * x2 + y1 * y2 + z1 * z2 + w1 * w2` and measures the similarity
    /// of the two vectors directions.
    ///
    /// # Parameters
    /// - `other`: The other `Vec4` to perfom the dot with.
    ///
    /// # Returns
    /// A scalar `f32` value representing the dot product.
    pub fn dot(&self, other: Vec4) -> f32 {
        self.x * other.x + self.y * other.y + self.z * other.z + self.w * other.w
    }

    /// Computes the cross product of two vectors.
    ///
    /// The cross product is defined as:
    /// ```text
    /// x = v1.y * v2.z - v1.z * v2.y
    /// y = v1.z * v2.x - v1.x * v2.z
    /// z = v1.x * v2.y - v1.y * v2.x
    /// ```
    ///
    /// # Parameters
    /// - `other`: The other vector in the equation.
    ///
    /// # Returns
    /// A new `Vec4` which is the result of the cross product between `self` and `other`.
    ///
    /// Note: This calculation ignores `w`.
    pub fn cross_xyz(&self, other: Vec4) -> Vec4 {
        Vec4::new(
            self.y * other.z - self.z * other.y, //x
            self.z * other.x - self.x * other.z, //y
            self.x * other.y - self.y * other.x, //z
            self.w,
        )
    }

    /// Computes the angle (in radians) between this vector and another.
    ///
    /// # Parameters
    /// - `other`: A 4D vector to compare against.
    ///
    /// # Returns
    /// Angle in radians between the vectors. Returns 0.0 if either vector is zero-length.
    pub fn angle_to(self, other: Vec4) -> f32 {
        let dot = self.dot(other);
        let len_product = self.length() * other.length();
        if epsilon_eq_default(len_product, 0.0) {
            return 0.0;
        }
        (dot / len_product).clamp(-1.0, 1.0).acos()
    }

    /// Calculates the angle in radians between `self` and another vector.
    ///
    /// The angle is in the range `[0, π]`.
    ///
    /// # Parameters
    /// - `other`: The other vector.
    ///
    /// # Returns
    /// The angle `f32` between the two vectors in radians.
    pub fn angle_between_radians(a: Vec4, b: Vec4) -> f32 {
        a.dot(b).clamp(-1.0, 1.0).acos()
    }

    /// Calculates the angle in degrees between `self` and another vector.
    ///
    /// The angle is in the range `[0°, 180°]`.
    ///
    /// # Parameters
    /// - `other`: The other vector.
    ///
    /// # Returns
    /// The angle `f32` between vectors in degrees.
    pub fn angle_between_degrees(a: Vec4, b: Vec4) -> f32 {
        utils::radians_to_degrees(a.dot(b).clamp(-1.0, 1.0).acos())
    }

    // ============= Interpolation =============

    /// Linearly interpolates between `self` and vector `b` by factor `t`.
    ///
    /// This performs a weighted average:
    /// `result = self * (1.0 - t) + b * t`
    ///
    /// # Parameters
    /// - `b`: The target vector to interpolate towards.
    /// - `t`: The interpolation factor in the range `[0.0, 1.0]`.
    ///
    /// # Returns
    /// A `Vec4` representing the point `t` fraction between `self` and `b`.
    pub fn lerp(self, target: Vec4, t: f32) -> Vec4 {
        self * (1.0 - t) + target * t
    }

    /// Clamped version of `lerp` that restricts `t` to the range `[0.0, 1.0]`.
    ///
    /// Prevents overshooting during interpolation.
    ///
    /// # Parameters
    /// - `b`: The target vector.
    /// - `t`: The interpolation factor, automatically clamped between `0.0` and `1.0`.
    ///
    /// # Returns
    /// A clamped interpolated `Vec4`.
    pub fn lerp_clamped(self, target: Vec4, t: f32) -> Vec4 {
        let t = utils::clamp(t, 0.0, 1.0);
        self.lerp(target, t)
    }

    /// Linearly interpolates between two vectors `a` and `b` by factor `t`.
    ///
    /// Unlike `lerp` and `lerp_clamped`, this is a static method that does not require `self`.
    ///
    /// # Parameters
    /// - `a`: The starting vector.
    /// - `b`: The target vector.
    /// - `t`: The interpolation factor.
    ///
    /// # Returns
    /// A vector interpolated between `a` and `b`.
    pub fn lerp_between(a: Vec4, b: Vec4, t: f32) -> Vec4 {
        a * (1.0 - t) + b * t
    }

    /// Same as `lerp_between`, but clamps `t` to `[0.0, 1.0]`.
    ///
    /// Ensures the result remains between `a` and `b`.
    ///
    /// # Parameters
    /// - `a`: The starting vector.
    /// - `b`: The target vector.
    /// - `t`: The interpolation factor, clamped.
    ///
    /// # Returns
    /// A clamped vector interpolated between `a` and `b`.
    pub fn lerp_between_clamped(a: Vec4, b: Vec4, t: f32) -> Vec4 {
        let t = utils::clamp(t, 0.0, 1.0);
        Vec4::lerp_between(a, b, t)
    }

    /// Performs spherical linear interpolation (SLERP) between two 4D vectors.
    ///
    /// SLERP smoothly interpolates between two directions over the surface of a sphere,
    /// preserving constant angular velocity and the length of the interpolated vector.
    /// This version handles vectors of different magnitudes and falls back to
    /// linear interpolation if the angle between them is too small.
    ///
    /// # Parameters
    /// - `a`: Starting vector.
    /// - `b`: Ending vector.
    /// - `t`: Interpolation factor in `[0.0, 1.0]`.
    ///
    /// # Returns
    /// A vector that represents the spherical interpolation between `a` and `b` at `t`.
    ///
    /// # Behavior
    /// - If either `a` or `b` is near-zero length, it blends based on magnitude.
    /// - If vectors are nearly parallel (dot > 1 - epsilon), falls back to clamped lerp.
    /// - Interpolates both direction and magnitude.
    ///
    /// # Example
    /// ```rust
    /// use omelet::vec4::Vec4;
    /// let a = Vec4::new(1.0, 0.0, 0.0, 0.0);
    /// let b = Vec4::new(0.0, 1.0, 0.0, 0.0);
    /// let halfway = Vec4::slerp(a, b, 0.5);
    /// ```
    pub fn slerp(a: Vec4, b: Vec4, t: f32) -> Vec4 {
        // Handle zero vectors
        if a.length() < 1e-6 {
            return b * t;
        }
        if b.length() < 1e-6 {
            return a * (1.0 - t);
        }

        let a_norm = a.normalize_or_zero();
        let b_norm = b.normalize_or_zero();
        let dot = a_norm.dot(b_norm).clamp(-1.0, 1.0);

        if dot > 1.0 - 1e-6 {
            return Vec4::lerp(a, b, t);
        }

        let theta = dot.acos();
        let sin_theta = theta.sin();
        let wa = ((1.0 - t) * theta).sin() / sin_theta;
        let wb = (t * theta).sin() / sin_theta;

        (a_norm * wa + b_norm * wb) * (a.length() * (1.0 - t) + b.length() * t)
    }

    /// Performs spherical linear interpolation (SLERP) between two 4D vectors
    /// using angle interpolation with shortest path logic.
    ///
    /// This function is conceptually similar to [`slerp`], but is optimised for
    /// interpolating orientations with the shortest arc. It assumes the vectors represent directions,
    /// and normalizes them before interpolation.
    ///
    /// # Parameters
    /// - `a`: Starting vector.
    /// - `b`: Ending vector.
    /// - `t`: Interpolation factor in `[0.0, 1.0]`.
    ///
    /// # Returns
    /// A vectopr representing the interpolated direction, scaled to smoothly blend
    /// the lengths of `a` and `b`.
    ///
    /// # Behavior
    /// - Normalizes both vectors before interpolation.
    /// - If either vector is near-zero, falls back to linear interpolation.
    /// - If the angle between vectors is too small, falls back to clamped lerp.
    ///
    /// # Example
    /// ```rust
    /// use omelet::vec4::Vec4;
    /// let start = Vec4::new(1.0, 0.0, 0.0, 0.0);
    /// let end = Vec4::new(0.0, 1.0, 0.0, 0.0);
    /// let mid = Vec4::slerp_angle(start, end, 0.5);
    /// ```
    pub fn slerp_angle(a: Vec4, b: Vec4, t: f32) -> Vec4 {
        // Handle zero vectors
        if a.length() < 1e-6 || b.length() < 1e-6 {
            return Vec4::lerp_between_clamped(a, b, t);
        }

        let a_norm = a.normalize();
        let b_norm = b.normalize();
        let dot = a_norm.dot(b_norm).clamp(-1.0, 1.0);

        // If the vectors are nearly parallel, use linear interpolation
        if dot > 1.0 - 1e-6 {
            return Vec4::lerp_between_clamped(a, b, t);
        }

        // Calculate the angle between the vectors
        let theta = dot.acos();

        // Interpolate using the slerp formula
        let sin_theta = theta.sin();
        let wa = ((1.0 - t) * theta).sin() / sin_theta;
        let wb = (t * theta).sin() / sin_theta;

        // Interpolate both direction and magnitude
        (a_norm * wa + b_norm * wb) * (a.length() * (1.0 - t) + b.length() * t)
    }

    // ============= Projection and Reflection =============

    /// Project `self` onto another vector `onto`.
    ///
    /// This returns the vector component of `self` that lies in the direction of `onto`.
    /// If `onto` is a zero vector, returns the zero vector to avoid division by zero.
    ///
    /// # Parameters
    /// - `onto`: The vector to project onto.
    ///
    /// # Returns
    /// The projection vector of `self` onto `onto`.
    pub fn project(&self, onto: Vec4) -> Vec4 {
        let denominator = onto.dot(onto);
        if denominator == 0.0 {
            Vec4::new(0.0, 0.0, 0.0, 0.0)
        } else {
            onto * (self.dot(onto) / denominator)
        }
    }

    /// Returns the component of `self` orthogonal to the vector `other`.
    ///
    /// The rejection is defined as `self - projection of self onto other`.
    ///
    /// # Parameters
    /// - `other`: The vector to reject from.
    ///
    /// # Returns
    /// A `Vec4` representing the component of `self` perpendicular to `other`.
    pub fn reject(&self, other: Vec4) -> Vec4 {
        *self - self.project(other)
    }

    /// Reflects the vector `self` over a given normal vector.
    ///
    /// This performs a mirror reflection of the vector abou the plane defined by the normal.
    /// The normal vector is assumed to be normalized.
    ///
    /// # Parameters
    /// - `normal`: The normalized normal vector to reflect about.
    ///
    /// # Returns
    /// A `Vec4` reflected.
    pub fn reflect(&self, normal: Vec4) -> Vec4 {
        *self - normal * (2.0 * self.dot(normal))
    }

    /// Returns the mirrored vector over a given normal.
    ///
    /// # Parameters
    /// - `normal`: The normal to mirror the vector over.
    ///
    /// # Returns
    /// `Vec4` which is mirrored over `normal`.
    pub fn mirror(&self, normal: Vec4) -> Vec4 {
        *self - normal * (2.0 * self.dot(normal))
    }

    // ============= Distance =============

    /// Returns the Euclidean distance between `self` and another vector.
    ///
    /// Calculated as the length of the difference vector.
    ///
    /// # Parameters
    /// - `other`: The target vector.
    ///
    /// # Returns
    /// A `f32` distance value.
    pub fn distance(&self, other: Vec4) -> f32 {
        (*self - other).length()
    }

    /// Returns the squared Euclidean distance between `self` and another vector.
    ///
    /// This avoids the cost of the square root and is useful when only relative comparisons
    /// of distance are required.
    ///
    /// # Parameters
    /// - `other`: The target vector.
    ///
    /// # Returns
    /// A `f32` squared distance value.
    pub fn squared_distance(&self, other: Vec4) -> f32 {
        (*self - other).squared_length()
    }

    /// Returns the normalized direction vector pointing from `self` towards `other`.
    ///
    /// # Parameters
    /// - `other`: The target vector.
    ///
    /// # Returns
    /// A unit length `Vec4` pointing from `self` to `other`.
    ///
    /// # Panics
    /// Panics if the direction vector length is zero (i.e. both vectors equal).
    pub fn direction_to(&self, other: Vec4) -> Vec4 {
        let delta = other - *self;
        delta.normalize_or_zero()
    }

    /// Returns the raw direction vector from `self` to `other` without normalization.
    ///
    /// # Parameters
    /// - `other`: The target vector.
    ///
    /// # Returns
    /// The difference vector `other - self`.
    pub fn direction_to_raw(&self, other: Vec4) -> Vec4 {
        other - *self
    }

    /// Moves the `current` vector towards the `target` vector by at most `max_delta`.
    ///
    /// If the distance between `current` and `target` is less than or equal to `max_delta`,
    /// returns `target`. Otherwise, returns a vector closer to `target` by `max_delta`.
    ///
    /// # Parameters
    /// - `current`: The starting position vector.
    /// - `target`: The target position vector.
    /// - `max_delta`: The maximum distance to move towards the target.
    ///
    /// # Returns
    /// A new `Vec4` moved towards the target by `max_delta` or exactly `target` if close enough.
    ///
    /// # Example
    /// ```rust
    /// use omelet::vec4::Vec4;
    /// let current = Vec4::new(0.0, 0.0, 0.0, 0.0);
    /// let target = Vec4::new(10.0, 0.0, 0.0, 0.0);
    /// let moved = Vec4::move_towards(current, target, 3.0);
    /// assert_eq!(moved, Vec4::new(3.0, 0.0, 0.0, 0.0));
    /// ```
    pub fn move_towards(current: Vec4, target: Vec4, max_delta: f32) -> Vec4 {
        let delta = target - current;
        let distance = delta.length();

        if distance <= max_delta || distance < f32::EPSILON {
            target
        } else {
            current + delta / distance * max_delta
        }
    }

    // ============= Geometry =============

    /// Returns two orthonormal vectors.
    ///
    /// Given trwo non-parallel vectors, returns two orthonormal vectors:
    /// the first is a normalized version of the input `a` and the second is
    /// perpendicular to `a` and lies in the same plane as `a` and `b`.
    ///
    /// # Parameters
    /// - `a`: First input vector.
    /// - `b`: Second input vector to orthonormalize relative to `a`.
    ///
    /// # Returns
    /// Tuple `(a_orthonormal, b_orthonormal)` where both vectors are unit length
    /// and perpendicular to each other.
    ///
    /// # Panics
    /// Panics if input vectors are zero or nearly parallel.
    pub fn orthonormalize(a: Vec4, b: Vec4) -> (Vec4, Vec4) {
        let a_norm = a.normalize();
        let projection = b.project(a_norm);
        let b_orthogonal = b - projection;
        let b_norm = b_orthogonal.normalize();
        (a_norm, b_norm)
    }

    /// Returns a vector rotated around an axis by a specific angle.
    ///
    /// # Parameters
    /// - `axis`: Axis to rotate around.
    /// - `angle`: Angle to rotate by.
    ///
    /// # Returns
    /// A `Vec4` rotated around the `axis` by `angle` using Rodrigues rotation formula.
    pub fn rotate_around_axis(&self, axis: Vec4, angle: f32) -> Vec4 {
        let axis = axis.normalize();
        let cos = angle.cos();
        let sin = angle.sin();
        *self * cos + axis.cross_xyz(*self) * sin + axis * axis.dot(*self) * (1.0 - cos)
    }

    /// Rotates this vector **within a 2D plane** formed by two axis indices in 4D space.
    ///
    /// Applies a rotation by the given angle in radians in the plane defined by axis `a` and `b`.
    ///
    /// # Parameters
    /// - `a`, `b`: Axis indices (0 = x, 1 = y, 2 = z, 3 = w); must be distinct.
    /// - `angle`: Angle of rotation in radians.
    ///
    /// # Panics
    /// If `a` or `b` are out of bounds or equal.
    ///
    /// # Returns
    /// The rotated Vec4.
    pub fn rotate_in_plane(self, a: usize, b: usize, angle: f32) -> Vec4 {
        assert!(
            a < 4 && b < 4 && a != b,
            "Invalid acis indices for 4D plane rotation"
        );

        let sin = angle.sin();
        let cos = angle.cos();

        let mut components = [self.x, self.y, self.z, self.w];
        let (va, vb) = (components[a], components[b]);

        components[a] = va * cos - vb * sin;
        components[b] = va * sin + vb * cos;

        Vec4::new(components[0], components[1], components[2], components[3])
    }

    // ============= Random =============

    /// Returns a randomly generated unit vector using rand crate.
    ///
    /// # Returns
    /// `Vec4` which each component is randomly generated in range `[0.0, 1.0]`.
    pub fn random_unit_vector() -> Vec4 {
        let mut rng = rand::rng();
        loop {
            let x = rng.random_range(-1.0..=1.0);
            let y = rng.random_range(-1.0..=1.0);
            let z = rng.random_range(-1.0..=1.0);
            let w = rng.random_range(-1.0..=1.0);

            let v = Vec4::new(x, y, z, w);
            let len_sq = v.squared_length();
            if len_sq > 0.0 && len_sq <= 1.0 {
                return v / len_sq.sqrt();
            }
        }
    }

    /// Generates a random point inside the 4D unit hypersphere (radius ≤ 1).
    ///
    /// This uses rejection sampling within a [-1, 1]^4 hypercube.
    /// The returned vector lies inside or on the boundary of the unit 4D sphere.
    ///
    /// # Returns
    /// A randomly generated Vec4 within the unit 4D sphere.
    pub fn random_in_unit_sphere() -> Vec4 {
        let mut rng = rand::rng();
        loop {
            let x = rng.random_range(-1.0..=1.0);
            let y = rng.random_range(-1.0..=1.0);
            let z = rng.random_range(-1.0..=1.0);
            let w = rng.random_range(-1.0..=1.0);
            let v = Vec4::new(x, y, z, w);
            if v.squared_length() <= 1.0 {
                return v;
            }
        }
    }

    // ============= Barycentric and Triangles =============

    /// Computes the barycentric coordinates of point `p` relative to triangle (`a`, `b`, `c`).
    ///
    /// Useful for interpolation inside triangles or determining whether a point lies in a triangle.
    ///
    /// # Parameters
    /// - `p`: The point to evaluate.
    /// - `a`, `b`, `c`: The triangles vertices.
    ///
    /// # Returns
    /// A tupel `(u, v, w)` of barycentric weight where `u + v + w = 1`.
    ///
    /// Note: Works properly if a, b, c, p all lie in the same plane.
    ///
    /// /// # Example
    /// ```rust
    /// use omelet::vec4::Vec4;
    /// let p = Vec4::new(2.0, 1.0, 0.0, 0.0);
    /// let a = Vec4::new(0.0, 0.0, 0.0, 0.0);
    /// let b = Vec4::new(4.0, 0.0, 0.0, 0.0);
    /// let c = Vec4::new(0.0, 4.0, 0.0, 0.0);
    /// let (u, v, w) = Vec4::barycentric(p, a, b, c);
    /// assert!((u + v + w - 1.0).abs() < 1e-6);
    /// ```
    pub fn barycentric(p: Vec4, a: Vec4, b: Vec4, c: Vec4) -> (f32, f32, f32) {
        let v0 = b - a;
        let v1 = c - a;
        let v2 = p - a;

        let d00 = v0.dot(v0);
        let d01 = v0.dot(v1);
        let d11 = v1.dot(v1);
        let d20 = v2.dot(v0);
        let d21 = v2.dot(v1);

        let denom = d00 * d11 - d01 * d01;
        if is_near_zero_default(denom) {
            return (0.0, 0.0, 0.0); // degenerate triangle
        }

        let v = (d11 * d20 - d01 * d21) / denom;
        let w = (d00 * d21 - d01 * d20) / denom;
        let u = 1.0 - v - w;

        (u, v, w)
    }

    /// Determines if a point lies within a triangle defined by three vertices.
    ///
    /// Uses barycentric coordinates to evaluate the points location relative to the triangle.
    ///
    /// # Parameters
    /// - `p`: The point to test.
    /// - `a`, `b`, `c`: The triangles corners.
    ///
    /// # Returns
    /// `true` if `p` lies inside or on the triangle, `false` otherwise.
    pub fn in_triangle(p: Vec4, a: Vec4, b: Vec4, c: Vec4) -> bool {
        let (u, v, w) = Vec4::barycentric(p, a, b, c);
        u >= 0.0 && v >= 0.0 && w >= 0.0 && u <= 1.0 && v <= 1.0 && w <= 1.0
    }

    // ============= Comparison and Validity =============

    /// Returns `true` if `self` and `other` are approx equal within a given epsilon tolerance.
    ///
    /// This compares each component individually and returns `true` only if both are within epsilon.
    ///
    /// # Parameters
    /// - `other`: The other vector to compare.
    /// - `epsilon`: The maximum allowed difference for each component.
    ///
    /// # Returns
    /// Boolean indicating approx equality.
    pub fn approx_eq(&self, other: Vec4) -> bool {
        epsilon_eq(self.x, other.x, 1e-6)
            && epsilon_eq(self.y, other.y, 1e-6)
            && epsilon_eq(self.z, other.z, 1e-6)
            && epsilon_eq(self.w, other.w, 1e-6)
    }

    /// Returns `true` if `self` and `other` are approx equal within `epsilon`.
    ///
    /// This compares each component individually and returns `true` only if both are within epsilon.
    ///
    /// # Parameters
    /// - `other`: The other vector to compare.
    /// - `epsilon`: The maximum allowed difference for each component.
    ///
    /// # Returns
    /// Boolean indicating approx equality.
    pub fn approx_eq_eps(&self, other: Vec4, epsilon: f32) -> bool {
        epsilon_eq(self.x, other.x, epsilon)
            && epsilon_eq(self.y, other.y, epsilon)
            && epsilon_eq(self.z, other.z, epsilon)
            && epsilon_eq(self.w, other.w, epsilon)
    }

    /// Checks whether all components of the vector are finite values.
    ///
    /// A component is considered finite if it is not `NaN` or infinite.
    ///
    /// # Returns
    /// `true` if both `x`, `y`, `z`, and `w` are finite numbers.
    pub fn is_finite(self) -> bool {
        self.x.is_finite() && self.y.is_finite() && self.z.is_finite() && self.w.is_finite()
    }

    /// Checks whether either component of the vector is `NaN`.
    ///
    /// # Returns
    /// `true` if `x` or `y` or `z` or `w` is `NaN`.
    pub fn is_nan(self) -> bool {
        self.x.is_nan() || self.y.is_nan() || self.z.is_nan() || self.w.is_nan()
    }

    /// Checks whether all components of the vector are approx zero.
    ///
    /// # Returns
    /// `true` if `x`, `y`, `z`, and `w` are approx equal to 0.0.
    pub fn is_zero(self) -> bool {
        epsilon_eq_default(self.x, 0.0)
            && epsilon_eq_default(self.y, 0.0)
            && epsilon_eq_default(self.z, 0.0)
            && epsilon_eq_default(self.w, 0.0)
    }

    /// Checks whether both components of the vector are approx equal to zero within a user-specified tolerance.
    ///
    /// Useful to check for effective zero while avoiding floating point precision issues.
    ///
    /// # Parameters
    /// - `epsilon`: The allowed margin of error.
    ///
    /// # Returns
    /// `true` if `x`, `y`, `z`, and `w` are approx equal to zero within epsilon, `false` otherwise.
    pub fn is_zero_eps(self, epsilon: f32) -> bool {
        epsilon_eq(self.x, 0.0, epsilon)
            && epsilon_eq(self.y, 0.0, epsilon)
            && epsilon_eq(self.z, 0.0, epsilon)
            && epsilon_eq(self.w, 0.0, epsilon)
    }
}

// ============= Operator Overloads =============

/// Adds two vectors together component-wise.
impl Add for Vec4 {
    type Output = Self;
    #[inline]
    fn add(self, rhs: Self) -> Self::Output {
        Self::new(
            self.x + rhs.x,
            self.y + rhs.y,
            self.z + rhs.z,
            self.w + rhs.w,
        )
    }
}

/// Adds a scalar to each component of the vector.
impl Add<f32> for Vec4 {
    type Output = Self;
    #[inline]
    fn add(self, rhs: f32) -> Self::Output {
        Self::new(self.x + rhs, self.y + rhs, self.z + rhs, self.w + rhs)
    }
}

/// Adds each component of a vector to a scalar.
impl Add<Vec4> for f32 {
    type Output = Vec4;
    #[inline]
    fn add(self, rhs: Vec4) -> Self::Output {
        Vec4::new(self + rhs.x, self + rhs.y, self + rhs.z, self + rhs.w)
    }
}

/// Subtracts `rhs` from `self` component-wise.
impl Sub for Vec4 {
    type Output = Self;
    #[inline]
    fn sub(self, rhs: Self) -> Self::Output {
        Self::new(
            self.x - rhs.x,
            self.y - rhs.y,
            self.z - rhs.z,
            self.w - rhs.w,
        )
    }
}

/// Subtracts a scalar from each component of the vector.
impl Sub<f32> for Vec4 {
    type Output = Self;
    #[inline]
    fn sub(self, rhs: f32) -> Self::Output {
        Self::new(self.x - rhs, self.y - rhs, self.z - rhs, self.w - rhs)
    }
}

/// Subtracts each component of a vector from a scalar.
impl Sub<Vec4> for f32 {
    type Output = Vec4;
    #[inline]
    fn sub(self, rhs: Vec4) -> Self::Output {
        Vec4::new(self - rhs.x, self - rhs.y, self - rhs.z, self - rhs.w)
    }
}

/// Multiplies two vectors together component-wise.
impl Mul for Vec4 {
    type Output = Self;
    #[inline]
    fn mul(self, rhs: Self) -> Self::Output {
        Self::new(
            self.x * rhs.x,
            self.y * rhs.y,
            self.z * rhs.z,
            self.w * rhs.w,
        )
    }
}

/// Multiplies each component of a vector by a scalar.
impl Mul<f32> for Vec4 {
    type Output = Self;
    #[inline]
    fn mul(self, rhs: f32) -> Self::Output {
        Self::new(self.x * rhs, self.y * rhs, self.z * rhs, self.w * rhs)
    }
}

/// Multiplies a scalar by each component of a vector.
impl Mul<Vec4> for f32 {
    type Output = Vec4;
    #[inline]
    fn mul(self, rhs: Vec4) -> Self::Output {
        Vec4::new(self * rhs.x, self * rhs.y, self * rhs.z, self * rhs.w)
    }
}

/// Divides `self` by `rhs` component-wise.
impl Div for Vec4 {
    type Output = Self;
    #[inline]
    fn div(self, rhs: Self) -> Self::Output {
        Self::new(
            self.x / rhs.x,
            self.y / rhs.y,
            self.z / rhs.z,
            self.w / rhs.w,
        )
    }
}

/// Divides each component of a vector by a scalar.
impl Div<f32> for Vec4 {
    type Output = Self;
    #[inline]
    fn div(self, rhs: f32) -> Self::Output {
        Self::new(self.x / rhs, self.y / rhs, self.z / rhs, self.w / rhs)
    }
}

/// Divides a scalar by each component of a vector.
impl Div<Vec4> for f32 {
    type Output = Vec4;
    #[inline]
    fn div(self, rhs: Vec4) -> Self::Output {
        Vec4::new(self / rhs.x, self / rhs.y, self / rhs.z, self / rhs.w)
    }
}

/// Negates each component of the vector.
impl Neg for Vec4 {
    type Output = Self;
    #[inline]
    fn neg(self) -> Self::Output {
        Self::new(-self.x, -self.y, -self.z, -self.w)
    }
}

// ============= Assignment Operator Overloads =============

impl AddAssign for Vec4 {
    #[inline]
    fn add_assign(&mut self, rhs: Self) {
        self.x += rhs.x;
        self.y += rhs.y;
        self.z += rhs.z;
        self.w += rhs.w;
    }
}

impl AddAssign<f32> for Vec4 {
    #[inline]
    fn add_assign(&mut self, rhs: f32) {
        self.x += rhs;
        self.y += rhs;
        self.z += rhs;
        self.w += rhs;
    }
}

impl SubAssign for Vec4 {
    #[inline]
    fn sub_assign(&mut self, rhs: Self) {
        self.x -= rhs.x;
        self.y -= rhs.y;
        self.z -= rhs.z;
        self.w -= rhs.w;
    }
}

impl SubAssign<f32> for Vec4 {
    #[inline]
    fn sub_assign(&mut self, rhs: f32) {
        self.x -= rhs;
        self.y -= rhs;
        self.z -= rhs;
        self.w -= rhs;
    }
}

impl MulAssign for Vec4 {
    #[inline]
    fn mul_assign(&mut self, rhs: Self) {
        self.x *= rhs.x;
        self.y *= rhs.y;
        self.z *= rhs.z;
        self.w *= rhs.w;
    }
}

impl MulAssign<f32> for Vec4 {
    #[inline]
    fn mul_assign(&mut self, rhs: f32) {
        self.x *= rhs;
        self.y *= rhs;
        self.z *= rhs;
        self.w *= rhs;
    }
}

impl DivAssign for Vec4 {
    #[inline]
    fn div_assign(&mut self, rhs: Self) {
        self.x /= rhs.x;
        self.y /= rhs.y;
        self.z /= rhs.z;
        self.w /= rhs.w;
    }
}

impl DivAssign<f32> for Vec4 {
    #[inline]
    fn div_assign(&mut self, rhs: f32) {
        self.x /= rhs;
        self.y /= rhs;
        self.z /= rhs;
        self.w /= rhs;
    }
}

// ============= Trait Implementations =============

impl Default for Vec4 {
    /// Returns a `Vec4` with all components set to zero.
    #[inline]
    fn default() -> Self {
        Self::ZERO // Assumes Vec4::ZERO constant exists
    }
}

/// Checks whether two vectors are exactly equal.
impl PartialEq for Vec4 {
    #[inline]
    fn eq(&self, other: &Self) -> bool {
        self.x == other.x && self.y == other.y && self.z == other.z && self.w == other.w
    }
}

/// Enables `v[index]` access. Panics if `index` is out of bounds.
impl Index<usize> for Vec4 {
    type Output = f32;
    #[inline]
    fn index(&self, index: usize) -> &Self::Output {
        match index {
            0 => &self.x,
            1 => &self.y,
            2 => &self.z,
            3 => &self.w,
            _ => panic!("Vec4 index out of bounds: {}", index),
        }
    }
}

/// Enables mutable `v[index]` access. Panics if `index` is out of bounds.
impl IndexMut<usize> for Vec4 {
    #[inline]
    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
        match index {
            0 => &mut self.x,
            1 => &mut self.y,
            2 => &mut self.z,
            3 => &mut self.w,
            _ => panic!("Vec4 index out of bounds: {}", index),
        }
    }
}

/// Implements the `Display` trait for pretty-printing.
impl fmt::Display for Vec4 {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(
            f,
            "Vec4({:.3}, {:.3}, {:.3}, {:.3})",
            self.x, self.y, self.z, self.w
        )
    }
}

// ============= Conversion Traits =============

/// Creates a `Vec4` from a tuple `(f32, f32, f32, f32)`.
impl From<(f32, f32, f32, f32)> for Vec4 {
    #[inline]
    fn from(t: (f32, f32, f32, f32)) -> Self {
        Self::new(t.0, t.1, t.2, t.3)
    }
}

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

/// Creates a `Vec4` from an array `[f32; 4]`.
impl From<[f32; 4]> for Vec4 {
    #[inline]
    fn from(arr: [f32; 4]) -> Self {
        Self::new(arr[0], arr[1], arr[2], arr[3])
    }
}

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

// ============= Approx Crate Implementations =============

/// Implements absolute difference equality comparison for `Vec4`.
impl approx::AbsDiffEq for Vec4 {
    type Epsilon = f32;

    #[inline]
    fn default_epsilon() -> f32 {
        f32::EPSILON
    }

    #[inline]
    fn abs_diff_eq(&self, other: &Self, epsilon: f32) -> bool {
        f32::abs_diff_eq(&self.x, &other.x, epsilon)
            && f32::abs_diff_eq(&self.y, &other.y, epsilon)
            && f32::abs_diff_eq(&self.z, &other.z, epsilon)
            && f32::abs_diff_eq(&self.w, &other.w, epsilon)
    }
}

/// Implements relative equality comparison for `Vec4`.
impl approx::RelativeEq for Vec4 {
    #[inline]
    fn default_max_relative() -> f32 {
        f32::EPSILON
    }

    #[inline]
    fn relative_eq(&self, other: &Self, epsilon: f32, max_relative: f32) -> bool {
        f32::relative_eq(&self.x, &other.x, epsilon, max_relative)
            && f32::relative_eq(&self.y, &other.y, epsilon, max_relative)
            && f32::relative_eq(&self.z, &other.z, epsilon, max_relative)
            && f32::relative_eq(&self.w, &other.w, epsilon, max_relative)
    }
}