-- ======================================================================
-- ENHANCED CONSTANTS (higher precision)
-- ======================================================================

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

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

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

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

-- IEEE 754 NaN check
function math.isnan(x)
	return x ~= x
end

-- Check if number is finite
function math.isfinite(x)
	return x > -math.infinity and x < math.infinity
end

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

-- Constrain value to range
function math.clamp(x, min, max)
	return x < min and min or (x > max and max or x)
end

-- Linear interpolation
function math.lerp(a, b, t)
	return a + (b - a) * t
end

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

-- Map value from input range to output range
function math.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
function math.round(x)
	return x >= 0 and math.floor(x + 0.5) or math.ceil(x - 0.5)
end

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

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

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

-- Factorial with bounds checking
function math.factorial(n)
	if n < 0 or n ~= math.floor(n) or n > 170 then
		return nil
	end
	local result = 1
	for i = 2, n do
		result = result * i
	end
	return result
end

-- Greatest common divisor using Euclidean algorithm
function math.gcd(a, b)
	a, b = math.floor(math.abs(a)), math.floor(math.abs(b))
	while b ~= 0 do
		a, b = b, a % b
	end
	return a
end

-- Least common multiple
function math.lcm(a, b)
	if a == 0 or b == 0 then return 0 end
	return math.abs(a * b) / math.gcd(a, b)
end

-- ======================================================================
-- ENHANCED RANDOM FUNCTIONS
-- ======================================================================

-- Random float in range
function math.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
function math.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 probability
function math.randboolean(p)
	p = p or 0.5
	return math.random() < p
end

-- ======================================================================
-- STATISTICS FUNCTIONS
-- ======================================================================

function math.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

function math.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

function math.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

function math.variance(t)
	if type(t) ~= "table" then return 0 end
	local count = 0
	local m = math.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

function math.stdev(t)
	return math.sqrt(math.variance(t))
end

function math.pvariance(t)
	if type(t) ~= "table" then return 0 end
	local count = 0
	local m = math.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

function math.pstdev(t)
	return math.sqrt(math.pvariance(t))
end

function math.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

function math.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

-- ======================================================================
-- 2D VECTOR OPERATIONS
-- ======================================================================

math.vec2 = {
	new = function(x, y)
		return {x = x or 0, y = y or 0}
	end,

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

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

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

	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,

	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 = function(a, b)
		return a.x * b.x + a.y * b.y
	end,

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

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

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

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

	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 = 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,

	angle = function(v)
		return math.atan2(v.y, v.x)
	end,

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

	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
-- ======================================================================

math.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 = 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.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
-- ======================================================================

math.mat2 = {
	new = function(a, b, c, d)
		return {
			{a or 1, b or 0},
			{c or 0, d or 1}
		}
	end,

	identity = function()
		return {{1, 0}, {0, 1}}
	end,

	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,

	det = function(m)
		return m[1][1] * m[2][2] - m[1][2] * m[2][1]
	end,

	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
		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,

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

	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,

	scale = function(sx, sy)
		sy = sy or sx
		return {
			{sx, 0},
			{0, sy}
		}
	end
}

math.mat3 = {
	identity = function()
		return {
			{1, 0, 0},
			{0, 1, 0},
			{0, 0, 1}
		}
	end,

	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,

	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_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 = function(x, y)
		return {
			{1, 0, x},
			{0, 1, y},
			{0, 0, 1}
		}
	end,

	rotation = function(angle)
		local cos, sin = math.cos(angle), math.sin(angle)
		return {
			{cos, -sin, 0},
			{sin, cos, 0},
			{0, 0, 1}
		}
	end,

	scale = function(sx, sy)
		sy = sy or sx
		return {
			{sx, 0, 0},
			{0, sy, 0},
			{0, 0, 1}
		}
	end,

	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
-- ======================================================================

math.geometry = {
	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.distance(px, py, x1, y1)
		end
		local t = ((px - x1) * dx + (py - y1) * dy) / len_sq
		t = math.clamp(t, 0, 1)
		local nearestX = x1 + t * dx
		local nearestY = y1 + t * dy
		return math.distance(px, py, nearestX, nearestY)
	end,

	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,

	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 = function(px, py, x1, y1, x2, y2, x3, y3)
		local area = math.geometry.triangle_area(x1, y1, x2, y2, x3, y3)
		local area1 = math.geometry.triangle_area(px, py, x2, y2, x3, y3)
		local area2 = math.geometry.triangle_area(x1, y1, px, py, x3, y3)
		local area3 = math.geometry.triangle_area(x1, y1, x2, y2, px, py)
		return math.abs(area - (area1 + area2 + area3)) < 1e-10
	end,

	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
		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,

	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.clamp(t, 0, 1)
		return x1 + t * dx, y1 + t * dy
	end
}

-- ======================================================================
-- INTERPOLATION FUNCTIONS
-- ======================================================================

math.interpolation = {
	bezier = function(t, p0, p1, p2, p3)
		t = math.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 = function(t, p0, p1, p2, p3)
		t = math.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 = function(t, p0, p1, m0, m1)
		t = math.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_bezier = function(t, p0, p1, p2)
		t = math.clamp(t, 0, 1)
		local mt = 1 - t
		return mt * mt * p0 + 2 * mt * t * p1 + t * t * p2
	end,

	step = function(t, edge, x)
		return t < edge and 0 or x
	end,

	smootherstep = function(edge0, edge1, x)
		local t = math.clamp((x - edge0) / (edge1 - edge0), 0, 1)
		return t * t * t * (t * (t * 6 - 15) + 10)
	end
}