-- math.lua - Extended math library with advanced functions and utilities

local math_ext = {}

-- Import standard math functions to maintain compatibility
for name, func in pairs(_G.math) do
	math_ext[name] = func
end

-- ======================================================================
-- CONSTANTS (higher precision than standard library)
-- ======================================================================

math_ext.pi = 3.14159265358979323846
math_ext.tau = 6.28318530717958647693 -- 2*pi, useful for full rotations
math_ext.e = 2.71828182845904523536
math_ext.phi = 1.61803398874989484820 -- Golden ratio (1 + sqrt(5)) / 2
math_ext.sqrt2 = 1.41421356237309504880
math_ext.sqrt3 = 1.73205080756887729353
math_ext.ln2 = 0.69314718055994530942 -- Natural log of 2
math_ext.ln10 = 2.30258509299404568402 -- Natural log of 10
math_ext.infinity = 1/0
math_ext.nan = 0/0

-- ======================================================================
-- EXTENDED FUNCTIONS
-- ======================================================================

-- Cube root that handles negative numbers correctly (unlike x^(1/3))
function math_ext.cbrt(x)
    return x < 0 and -(-x)^(1/3) or x^(1/3)
end

-- Euclidean distance (hypotenuse) - more accurate than sqrt(x*x + y*y) for edge cases
function math_ext.hypot(x, y)
    return math.sqrt(x * x + y * y)
end

-- IEEE 754 NaN check - only reliable way to test for NaN
function math_ext.isnan(x)
    return x ~= x
end

-- Check if number is finite (not infinity or NaN)
function math_ext.isfinite(x)
    return x > -math_ext.infinity and x < math_ext.infinity
end

-- Mathematical sign function returning -1, 0, or 1
function math_ext.sign(x)
    return x > 0 and 1 or (x < 0 and -1 or 0)
end

-- Constrain value to range [min, max] - essential for safe calculations
function math_ext.clamp(x, min, max)
    return x < min and min or (x > max and max or x)
end

-- Linear interpolation - fundamental for animation and gradients
function math_ext.lerp(a, b, t)
    return a + (b - a) * t
end

-- Smooth interpolation using Hermite polynomial - gives eased motion
function math_ext.smoothstep(a, b, t)
    t = math_ext.clamp((t - a) / (b - a), 0, 1)
    return t * t * (3 - 2 * t)
end

-- Map value from input range to output range - useful for scaling
function math_ext.map(x, in_min, in_max, out_min, out_max)
    return (x - in_min) * (out_max - out_min) / (in_max - in_min) + out_min
end

-- Round to nearest integer (more predictable than math.floor(x + 0.5))
function math_ext.round(x)
    return x >= 0 and math.floor(x + 0.5) or math.ceil(x - 0.5)
end

-- Round to specified decimal places using multiplication/division
function math_ext.roundto(x, decimals)
    local mult = 10 ^ (decimals or 0)
    return math.floor(x * mult + 0.5) / mult
end

-- Normalize angle to [-π, π] range for consistent angle calculations
function math_ext.normalize_angle(angle)
    return angle - 2 * math_ext.pi * math.floor((angle + math_ext.pi) / (2 * math_ext.pi))
end

-- 2D Euclidean distance between two points
function math_ext.distance(x1, y1, x2, y2)
    local dx, dy = x2 - x1, y2 - y1
    return math.sqrt(dx * dx + dy * dy)
end

-- ======================================================================
-- RANDOM NUMBER FUNCTIONS - Enhanced random utilities
-- ======================================================================

-- Random float in specified range [min, max) - more flexible than math.random()
function math_ext.randomf(min, max)
    if not min and not max then
        return math.random()
    elseif not max then
        max = min
        min = 0
    end
    return min + math.random() * (max - min)
end

-- Random integer in range [min, max] inclusive - handles single argument case
function math_ext.randint(min, max)
    if not max then
        max = min
        min = 1
    end
    return math.floor(math.random() * (max - min + 1) + min)
end

-- Random boolean with configurable probability - useful for procedural generation
function math_ext.randboolean(p)
    p = p or 0.5
    return math.random() < p
end

-- ======================================================================
-- STATISTICS FUNCTIONS - Robust statistical calculations
-- ======================================================================

-- Sum of numeric values in table - filters out non-numbers automatically
function math_ext.sum(t)
    if type(t) ~= "table" then return 0 end

    local sum = 0
    for i=1, #t do
        if type(t[i]) == "number" then
            sum = sum + t[i]
        end
    end
    return sum
end

-- Arithmetic mean (average) - handles mixed-type arrays safely
function math_ext.mean(t)
    if type(t) ~= "table" or #t == 0 then return 0 end

    local sum = 0
    local count = 0
    for i=1, #t do
        if type(t[i]) == "number" then
            sum = sum + t[i]
            count = count + 1
        end
    end
    return count > 0 and sum / count or 0
end

-- Median value - sorts data to find middle value (handles even/odd lengths)
function math_ext.median(t)
    if type(t) ~= "table" or #t == 0 then return 0 end

    local nums = {}
    local count = 0
    for i=1, #t do
        if type(t[i]) == "number" then
            count = count + 1
            nums[count] = t[i]
        end
    end

    if count == 0 then return 0 end

    table.sort(nums)

    if count % 2 == 0 then
        return (nums[count/2] + nums[count/2 + 1]) / 2
    else
        return nums[math.ceil(count/2)]
    end
end

-- Sample variance - measures spread using n-1 denominator (Bessel's correction)
function math_ext.variance(t)
    if type(t) ~= "table" then return 0 end

    local count = 0
    local m = math_ext.mean(t)
    local sum = 0

    for i=1, #t do
        if type(t[i]) == "number" then
            local dev = t[i] - m
            sum = sum + dev * dev
            count = count + 1
        end
    end

    return count > 1 and sum / count or 0
end

-- Sample standard deviation - square root of variance
function math_ext.stdev(t)
    return math.sqrt(math_ext.variance(t))
end

-- Population variance - uses n denominator instead of n-1
function math_ext.pvariance(t)
    if type(t) ~= "table" then return 0 end

    local count = 0
    local m = math_ext.mean(t)
    local sum = 0

    for i=1, #t do
        if type(t[i]) == "number" then
            local dev = t[i] - m
            sum = sum + dev * dev
            count = count + 1
        end
    end

    return count > 0 and sum / count or 0
end

-- Population standard deviation
function math_ext.pstdev(t)
    return math.sqrt(math_ext.pvariance(t))
end

-- Mode - most frequently occurring value (first encountered if tie)
function math_ext.mode(t)
    if type(t) ~= "table" or #t == 0 then return nil end

    local counts = {}
    local most_frequent = nil
    local max_count = 0

    for i=1, #t do
        local v = t[i]
        counts[v] = (counts[v] or 0) + 1
        if counts[v] > max_count then
            max_count = counts[v]
            most_frequent = v
        end
    end

    return most_frequent
end

-- Simultaneous min/max - more efficient than separate calls
function math_ext.minmax(t)
    if type(t) ~= "table" or #t == 0 then return nil, nil end

    local min, max
    for i=1, #t do
        if type(t[i]) == "number" then
            min = t[i]
            max = t[i]
            break
        end
    end

    if min == nil then return nil, nil end

    for i=1, #t do
        if type(t[i]) == "number" then
            if t[i] < min then min = t[i] end
            if t[i] > max then max = t[i] end
        end
    end

    return min, max
end

-- ======================================================================
-- VECTOR OPERATIONS - 2D/3D vector math for graphics and physics
-- ======================================================================

-- 2D Vector operations - fundamental for 2D graphics, physics, and UI
math_ext.vec2 = {
    -- Create new 2D vector with default zero values
    new = function(x, y)
        return {x = x or 0, y = y or 0}
    end,

    -- Create independent copy to avoid reference issues
    copy = function(v)
        return {x = v.x, y = v.y}
    end,

    -- Vector addition - combines two displacement vectors
    add = function(a, b)
        return {x = a.x + b.x, y = a.y + b.y}
    end,

    -- Vector subtraction - difference between two points/vectors
    sub = function(a, b)
        return {x = a.x - b.x, y = a.y - b.y}
    end,

    -- Scalar or component-wise multiplication
    mul = function(a, b)
        if type(b) == "number" then
            return {x = a.x * b, y = a.y * b}
        end
        return {x = a.x * b.x, y = a.y * b.y}
    end,

    -- Scalar or component-wise division (uses multiplication for efficiency)
    div = function(a, b)
        if type(b) == "number" then
            local inv = 1 / b
            return {x = a.x * inv, y = a.y * inv}
        end
        return {x = a.x / b.x, y = a.y / b.y}
    end,

    -- Dot product - measures vector similarity/projection
    dot = function(a, b)
        return a.x * b.x + a.y * b.y
    end,

    -- Euclidean length/magnitude of vector
    length = function(v)
        return math.sqrt(v.x * v.x + v.y * v.y)
    end,

    -- Squared length - avoids sqrt for performance when comparing lengths
    length_squared = function(v)
        return v.x * v.x + v.y * v.y
    end,

    -- Distance between two points
    distance = function(a, b)
        local dx, dy = b.x - a.x, b.y - a.y
        return math.sqrt(dx * dx + dy * dy)
    end,

    -- Squared distance - avoids sqrt for performance
    distance_squared = function(a, b)
        local dx, dy = b.x - a.x, b.y - a.y
        return dx * dx + dy * dy
    end,

    -- Convert to unit vector (length 1) - preserves direction
    normalize = function(v)
        local len = math.sqrt(v.x * v.x + v.y * v.y)
        if len > 1e-10 then
            local inv_len = 1 / len
            return {x = v.x * inv_len, y = v.y * inv_len}
        end
        return {x = 0, y = 0}
    end,

    -- Rotate vector by angle (counterclockwise)
    rotate = function(v, angle)
        local c, s = math.cos(angle), math.sin(angle)
        return {
            x = v.x * c - v.y * s,
            y = v.x * s + v.y * c
        }
    end,

    -- Get angle of vector from positive x-axis
    angle = function(v)
        return math.atan2(v.y, v.x)
    end,

    -- Linear interpolation between two vectors
    lerp = function(a, b, t)
        t = math_ext.clamp(t, 0, 1)
        return {
            x = a.x + (b.x - a.x) * t,
            y = a.y + (b.y - a.y) * t
        }
    end,

    -- Reflect vector across normal (like light bouncing off surface)
    reflect = function(v, normal)
        local dot = v.x * normal.x + v.y * normal.y
        return {
            x = v.x - 2 * dot * normal.x,
            y = v.y - 2 * dot * normal.y
        }
    end
}

-- 3D Vector operations - essential for 3D graphics and spatial calculations
math_ext.vec3 = {
    new = function(x, y, z)
        return {x = x or 0, y = y or 0, z = z or 0}
    end,

    copy = function(v)
        return {x = v.x, y = v.y, z = v.z}
    end,

    add = function(a, b)
        return {x = a.x + b.x, y = a.y + b.y, z = a.z + b.z}
    end,

    sub = function(a, b)
        return {x = a.x - b.x, y = a.y - b.y, z = a.z - b.z}
    end,

    mul = function(a, b)
        if type(b) == "number" then
            return {x = a.x * b, y = a.y * b, z = a.z * b}
        end
        return {x = a.x * b.x, y = a.y * b.y, z = a.z * b.z}
    end,

    div = function(a, b)
        if type(b) == "number" then
            local inv = 1 / b
            return {x = a.x * inv, y = a.y * inv, z = a.z * inv}
        end
        return {x = a.x / b.x, y = a.y / b.y, z = a.z / b.z}
    end,

    dot = function(a, b)
        return a.x * b.x + a.y * b.y + a.z * b.z
    end,

    -- Cross product - creates perpendicular vector (right-hand rule)
    cross = function(a, b)
        return {
            x = a.y * b.z - a.z * b.y,
            y = a.z * b.x - a.x * b.z,
            z = a.x * b.y - a.y * b.x
        }
    end,

    length = function(v)
        return math.sqrt(v.x * v.x + v.y * v.y + v.z * v.z)
    end,

    length_squared = function(v)
        return v.x * v.x + v.y * v.y + v.z * v.z
    end,

    distance = function(a, b)
        local dx, dy, dz = b.x - a.x, b.y - a.y, b.z - a.z
        return math.sqrt(dx * dx + dy * dy + dz * dz)
    end,

    distance_squared = function(a, b)
        local dx, dy, dz = b.x - a.x, b.y - a.y, b.z - a.z
        return dx * dx + dy * dy + dz * dz
    end,

    normalize = function(v)
        local len = math.sqrt(v.x * v.x + v.y * v.y + v.z * v.z)
        if len > 1e-10 then
            local inv_len = 1 / len
            return {x = v.x * inv_len, y = v.y * inv_len, z = v.z * inv_len}
        end
        return {x = 0, y = 0, z = 0}
    end,

    lerp = function(a, b, t)
        t = math_ext.clamp(t, 0, 1)
        return {
            x = a.x + (b.x - a.x) * t,
            y = a.y + (b.y - a.y) * t,
            z = a.z + (b.z - a.z) * t
        }
    end,

    reflect = function(v, normal)
        local dot = v.x * normal.x + v.y * normal.y + v.z * normal.z
        return {
            x = v.x - 2 * dot * normal.x,
            y = v.y - 2 * dot * normal.y,
            z = v.z - 2 * dot * normal.z
        }
    end
}

-- ======================================================================
-- MATRIX OPERATIONS - Linear transformations for graphics and math
-- ======================================================================

math_ext.mat2 = {
    -- Create 2x2 matrix with specified values (defaults to identity)
    new = function(a, b, c, d)
        return {
            {a or 1, b or 0},
            {c or 0, d or 1}
        }
    end,

    -- 2x2 identity matrix - no transformation
    identity = function()
        return {{1, 0}, {0, 1}}
    end,

    -- Matrix multiplication - combines transformations (order matters)
    mul = function(a, b)
        return {
            {
                a[1][1] * b[1][1] + a[1][2] * b[2][1],
                a[1][1] * b[1][2] + a[1][2] * b[2][2]
            },
            {
                a[2][1] * b[1][1] + a[2][2] * b[2][1],
                a[2][1] * b[1][2] + a[2][2] * b[2][2]
            }
        }
    end,

    -- Determinant - measures area scaling factor
    det = function(m)
        return m[1][1] * m[2][2] - m[1][2] * m[2][1]
    end,

    -- Inverse matrix - reverses transformation (if possible)
    inverse = function(m)
        local det = m[1][1] * m[2][2] - m[1][2] * m[2][1]
        if math.abs(det) < 1e-10 then
            return nil  -- Matrix is not invertible
        end

        local inv_det = 1 / det
        return {
            {m[2][2] * inv_det, -m[1][2] * inv_det},
            {-m[2][1] * inv_det, m[1][1] * inv_det}
        }
    end,

    -- Create rotation matrix for given angle
    rotation = function(angle)
        local cos, sin = math.cos(angle), math.sin(angle)
        return {
            {cos, -sin},
            {sin, cos}
        }
    end,

    -- Apply matrix transformation to 2D vector
    transform = function(m, v)
        return {
            x = m[1][1] * v.x + m[1][2] * v.y,
            y = m[2][1] * v.x + m[2][2] * v.y
        }
    end,

    -- Create scaling matrix (uniform if sy omitted)
    scale = function(sx, sy)
        sy = sy or sx
        return {
            {sx, 0},
            {0, sy}
        }
    end
}

math_ext.mat3 = {
    -- 3x3 identity matrix - useful for 2D transformations with translation
    identity = function()
        return {
            {1, 0, 0},
            {0, 1, 0},
            {0, 0, 1}
        }
    end,

    -- Complete 2D transformation matrix (translate, rotate, scale)
    transform = function(x, y, angle, sx, sy)
        sx = sx or 1
        sy = sy or sx
        local cos, sin = math.cos(angle), math.sin(angle)
        return {
            {cos * sx, -sin * sy, x},
            {sin * sx, cos * sy, y},
            {0, 0, 1}
        }
    end,

    -- 3x3 matrix multiplication - combines transformations
    mul = function(a, b)
        local result = {
            {0, 0, 0},
            {0, 0, 0},
            {0, 0, 0}
        }

        for i = 1, 3 do
            for j = 1, 3 do
                for k = 1, 3 do
                    result[i][j] = result[i][j] + a[i][k] * b[k][j]
                end
            end
        end

        return result
    end,

    -- Transform 2D point using homogeneous coordinates
    transform_point = function(m, v)
        local x = m[1][1] * v.x + m[1][2] * v.y + m[1][3]
        local y = m[2][1] * v.x + m[2][2] * v.y + m[2][3]
        local w = m[3][1] * v.x + m[3][2] * v.y + m[3][3]

        if math.abs(w) < 1e-10 then
            return {x = 0, y = 0}
        end

        return {x = x / w, y = y / w}
    end,

    -- Translation matrix
    translation = function(x, y)
        return {
            {1, 0, x},
            {0, 1, y},
            {0, 0, 1}
        }
    end,

    -- 2D rotation matrix in 3x3 form
    rotation = function(angle)
        local cos, sin = math.cos(angle), math.sin(angle)
        return {
            {cos, -sin, 0},
            {sin, cos, 0},
            {0, 0, 1}
        }
    end,

    -- Scaling matrix in 3x3 form
    scale = function(sx, sy)
        sy = sy or sx
        return {
            {sx, 0, 0},
            {0, sy, 0},
            {0, 0, 1}
        }
    end,

    -- 3x3 determinant - measures volume scaling
    det = function(m)
        return 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])
    end
}

-- ======================================================================
-- GEOMETRY FUNCTIONS - Computational geometry utilities
-- ======================================================================

math_ext.geometry = {
    -- Shortest distance from point to line segment
    point_line_distance = function(px, py, x1, y1, x2, y2)
        local dx, dy = x2 - x1, y2 - y1
        local len_sq = dx * dx + dy * dy

        if len_sq < 1e-10 then
            return math_ext.distance(px, py, x1, y1)
        end

        local t = ((px - x1) * dx + (py - y1) * dy) / len_sq
        t = math_ext.clamp(t, 0, 1)

        local nearestX = x1 + t * dx
        local nearestY = y1 + t * dy

        return math_ext.distance(px, py, nearestX, nearestY)
    end,

    -- Point-in-polygon test using ray casting algorithm
    point_in_polygon = function(px, py, vertices)
        local inside = false
        local n = #vertices / 2

        for i = 1, n do
            local x1, y1 = vertices[i*2-1], vertices[i*2]
            local x2, y2

            if i == n then
                x2, y2 = vertices[1], vertices[2]
            else
                x2, y2 = vertices[i*2+1], vertices[i*2+2]
            end

            if ((y1 > py) ~= (y2 > py)) and
               (px < (x2 - x1) * (py - y1) / (y2 - y1) + x1) then
                inside = not inside
            end
        end

        return inside
    end,

    -- Calculate triangle area using cross product method
    triangle_area = function(x1, y1, x2, y2, x3, y3)
        return math.abs((x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)) / 2)
    end,

    -- Point-in-triangle test using barycentric coordinates
    point_in_triangle = function(px, py, x1, y1, x2, y2, x3, y3)
        local area = math_ext.geometry.triangle_area(x1, y1, x2, y2, x3, y3)
        local area1 = math_ext.geometry.triangle_area(px, py, x2, y2, x3, y3)
        local area2 = math_ext.geometry.triangle_area(x1, y1, px, py, x3, y3)
        local area3 = math_ext.geometry.triangle_area(x1, y1, x2, y2, px, py)

        return math.abs(area - (area1 + area2 + area3)) < 1e-10
    end,

    -- Line segment intersection using parametric equations
    line_intersect = function(x1, y1, x2, y2, x3, y3, x4, y4)
        local d = (y4 - y3) * (x2 - x1) - (x4 - x3) * (y2 - y1)

        if math.abs(d) < 1e-10 then
            return false, nil, nil -- Lines are parallel
        end

        local ua = ((x4 - x3) * (y1 - y3) - (y4 - y3) * (x1 - x3)) / d
        local ub = ((x2 - x1) * (y1 - y3) - (y2 - y1) * (x1 - x3)) / d

        if ua >= 0 and ua <= 1 and ub >= 0 and ub <= 1 then
            local x = x1 + ua * (x2 - x1)
            local y = y1 + ua * (y2 - y1)
            return true, x, y
        end

        return false, nil, nil
    end,

    -- Find closest point on line segment to given point
    closest_point_on_segment = function(px, py, x1, y1, x2, y2)
        local dx, dy = x2 - x1, y2 - y1
        local len_sq = dx * dx + dy * dy

        if len_sq < 1e-10 then
            return x1, y1
        end

        local t = ((px - x1) * dx + (py - y1) * dy) / len_sq
        t = math_ext.clamp(t, 0, 1)

        return x1 + t * dx, y1 + t * dy
    end
}

-- ======================================================================
-- INTERPOLATION FUNCTIONS - Smooth transitions and curve generation
-- ======================================================================

math_ext.interpolation = {
    -- Cubic Bézier curve - standard for smooth animations
    bezier = function(t, p0, p1, p2, p3)
        t = math_ext.clamp(t, 0, 1)
        local t2 = t * t
        local t3 = t2 * t
        local mt = 1 - t
        local mt2 = mt * mt
        local mt3 = mt2 * mt

        return p0 * mt3 + 3 * p1 * mt2 * t + 3 * p2 * mt * t2 + p3 * t3
    end,

    -- Catmull-Rom spline - smooth interpolation through points
    catmull_rom = function(t, p0, p1, p2, p3)
        t = math_ext.clamp(t, 0, 1)
        local t2 = t * t
        local t3 = t2 * t

        return 0.5 * (
            (2 * p1) +
            (-p0 + p2) * t +
            (2 * p0 - 5 * p1 + 4 * p2 - p3) * t2 +
            (-p0 + 3 * p1 - 3 * p2 + p3) * t3
        )
    end,

    -- Hermite interpolation with tangent control
    hermite = function(t, p0, p1, m0, m1)
        t = math_ext.clamp(t, 0, 1)
        local t2 = t * t
        local t3 = t2 * t
        local h00 = 2 * t3 - 3 * t2 + 1
        local h10 = t3 - 2 * t2 + t
        local h01 = -2 * t3 + 3 * t2
        local h11 = t3 - t2

        return h00 * p0 + h10 * m0 + h01 * p1 + h11 * m1
    end,

    -- Quadratic Bézier - simpler curve with one control point
    quadratic_bezier = function(t, p0, p1, p2)
        t = math_ext.clamp(t, 0, 1)
        local mt = 1 - t
        return mt * mt * p0 + 2 * mt * t * p1 + t * t * p2
    end,

    -- Step function - instant transition at threshold
    step = function(t, edge, x)
        return t < edge and 0 or x
    end,

    -- Smoothstep - S-curve interpolation (3rd order)
    smoothstep = function(edge0, edge1, x)
        local t = math_ext.clamp((x - edge0) / (edge1 - edge0), 0, 1)
        return t * t * (3 - 2 * t)
    end,

    -- Smootherstep - Even smoother S-curve (5th order, Ken Perlin)
    smootherstep = function(edge0, edge1, x)
        local t = math_ext.clamp((x - edge0) / (edge1 - edge0), 0, 1)
        return t * t * t * (t * (t * 6 - 15) + 10)
    end
}

return math_ext