quartic

local function IsZero(d)
	local eps = 1e-9
	return (d > -eps and d < eps)
end

function solveQuadric(c0, c1, c2)
	local s0, s1
	local p, q, D
	local num

	-- normal form: x^2 + px + q = 0
	p = c1 / (2 * c0)
	q = c2 / c0

	D = p * p - q

	if IsZero(D) then
		s0 = -p
		num = 1
	elseif (D < 0) then
		num = 0
	else
		local sqrt_D = math.sqrt(D)

		s0 = sqrt_D - p
		s1 = -sqrt_D - p
		num = 2
	end

	return num, s0, s1
end

function solveCubic(c0, c1, c2, c3)
	local s0, s1, s2
	local num
	local sub
	local A, B, C
	local sq_A, p, q
	local cb_p, D

	-- normal form: x^3 + Ax^2 + Bx + C = 0
	A = c1 / c0
	B = c2 / c0
	C = c3 / c0

	-- substitute x = y - A/3 to eliminate quadric term: x^3 + px + q = 0
	sq_A = A * A
	p = (1 / 3) * (-(1 / 3) + sq_A + B)
	q = (1 / 3) * ((2 / 27) * A * sq_A - (1 / 3) * A * B + C)

	-- use Cardano's formula
	cb_p = p * p * p
	D = q * q + cb_p

	if IsZero(D) then
		if IsZero(q) then -- one triple solution
			s0 = 0
			num = 1
		else -- one single and one double solution
			local u = -q ^ (1 / 3)
			s0 = 2 * u
			s1 = -u
			num = 2
		end
	elseif (D < 0) then -- casus irreducibilis
		local phi = (1 / 3) * math.acos(-q / math.sqrt(-cb_p))
		local t = 2 * math.sqrt(-p)

		s0 = t * math.cos(phi)
		s1 = -t * math.cos(phi + math.pi / 3)
		s2 = -t * math.cos(phi - math.pi / 3)
		num = 3
	else -- one real solution
		local sqrt_D = math.sqrt(D)
		local u = (sqrt_D - q) ^ (1 / 3)
		local v = (sqrt_D + q) ^ (1 / 3)

		s0 = u + v
		num = 1
	end

	-- resubstitute
	sub = (1 / 3) * A
	if (num > 0) then s0 = s0 - sub end
	if (num > 1) then s1 = s1 - sub end
	if (num > 2) then s2 = s2 - sub end

	return num, s0, s1, s2
end

function solveQuartic(c0, c1, c2, c3, c4)
	local s0, s1, s2, s3

	local coeffs = {}
	local z, u, v, sub
	local A, B, C, D
	local sq_A, p, q, r
	local num

	A = c1 / c0
	B = c2 / c0
	C = c3 / c0
	D = c4 / c0

	sq_A = A * A
	p = -(3 / 8) * sq_A + B
	q = (1 / 8) * sq_A * A - (1 / 2) * A * B + C
	r = -(3 / 256) * sq_A * sq_A + (1 / 16) * sq_A * B - (1 / 4) * A * C + D

	if IsZero(r) then
		coeffs[4] = q
		coeffs[3] = p
		coeffs[2] = 0
		coeffs[1] = 1

		num, s0, s1, s2 = solveCubic(coeffs[1], coeffs[2], coeffs[3], coeffs[4])
	else
		coeffs[4] = (1 / 2) * r * p - (1 / 8) * q * q
		coeffs[3] = -r
		coeffs[2] = -(1 / 2) * p
		coeffs[1] = 1

		local drop
		drop, s0, s1, s2 = solveCubic(coeffs[1], coeffs[2], coeffs[3], coeffs[4])
		drop = nil

		z = s0

		u = z * z - r
		v = 2 * z - p

		if IsZero(u) then
			u = 0
		elseif (u > 0) then
			u = math.sqrt(u)
		else
			num = 0
		end

		if IsZero(v) then
			v = 0
		elseif (v > 0) then
			v = math.sqrt(v)
		else
			num = 0
		end

		coeffs[3] = z - u
		coeffs[2] = (q < 0) and -v or v
		coeffs[1] = 1

		num, s0, s1 = solveQuadric(coeffs[1], coeffs[2], coeffs[3])

		coeffs[3] = z + u
		coeffs[2] = (q < 0) and v or -v
		coeffs[1] = 1

		if (num == 0) then num = num + solveQuadric(coeffs[1], coeffs[2], coeffs[3]) end
		if (num == 1) then num = num + solveQuadric(coeffs[1], coeffs[2], coeffs[3]) end
		if (num == 2) then num = num + solveQuadric(coeffs[1], coeffs[2], coeffs[3]) end
	end

	-- resubstitute
	sub = (1 / 4) * A

	if (num > 0) then s0 = s0 - sub end
	if (num > 1) then s1 = s1 - sub end
	if (num > 2) then s2 = s2 - sub end
	if (num > 3) then s3 = s3 - sub end

	return num, s0, s1, s2, s3
end

return solveQuartic(5, -3, 6, 9, -4)