Poly

--//Port of http://www.realtimerendering.com/resources/GraphicsGems/gems/Roots3And4.c from Graphics Gems

local Polynomial = {}

local math = math
local zero = 1e-30

local onethird = 1/3
local twotwsevs = 2/27
local threetwohund = 3/256
local onesixteenth = 1/16
local oneeighth = .125
local threeeighths = .375
local onehalf = 0.5
local onefourth = 0.25

local thirdpi = math.pi / 3

local function isZero(num)
	return num < zero and num > -zero
end

local function cubicRoot(num)
	return (num > 0 and math.pow(num, onethird)) or (num < 0 and -math.pow(-num, onethird)) or 0
end

function Polynomial.SolveQuadric(c0, c1, c2, reftable)
	local p, q, D

	p = c1 / (2 * c0)
	q = c2 / c0

	D = p * p - q

	if isZero(D) then
		reftable[1] = -p
		return 1
	elseif D < 0 then
		return 0
	else
		local sqr = math.sqrt(D)

		reftable[1] = sqr - p
		reftable[2] = -sqr - p

		return 2
	end
end

function Polynomial.SolveCubic(c0, c1, c2, c3, reftable)
	local num, sub, A, B, C, sqA, p, q, cbp, D

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

	sqA = A * A
	p = onethird * (-onethird * sqA + B)
	q = onehalf * (twotwsevs * A * sqA - onethird * A * B + C)

	cbp = p * p * p
	D = q * q + cbp

	sub = onethird * A

	if isZero(D) then
		if isZero(q) then
			reftable[1] = 0 - sub
			num = 1
		else
			local u = cubicRoot(-q)

			reftable[1] = 2 * u - sub
			reftable[2] = -u - sub

			num = 2
		end
	elseif D < 0 then
		local phi = onethird * math.acos(-q / math.sqrt(-cbp))
		local t = 2 * math.sqrt(-p)

		reftable[1] = t * math.cos(phi) - sub
		reftable[2] = -t * math.cos(phi + thirdpi) - sub
		reftable[3] = -t * math.cos(phi - thirdpi) - sub

		num = 3
	else
		local sqrD = math.sqrt(D)

		reftable[1] = cubicRoot(sqrD - q) - cubicRoot(sqrD + q) - sub
		num = 1
	end

	return num
end

function Polynomial.SolveQuartic(c0, c1, c2, c3, c4, reftable)
	local num, z, u, v, sub, A, B, C, D, sqA, p, q, r

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

	sqA = A * A
	p = -threeeighths * sqA + B
	q = oneeighth * sqA * A - onehalf * A * B + C
	r = -threetwohund * sqA * sqA + onesixteenth * sqA * B - onefourth * A * C + D

	sub = onefourth * A

	if isZero(r) then
		num = Polynomial.SolveCubic(1, 0, p, q, reftable)
	else
		local holder = {}
		local h, j, k

		Polynomial.SolveCubic(1, -onehalf * p, -r, onehalf * r * p - oneeighth * q * q, reftable)
		z = reftable[1]

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

		u = (isZero(u) and 0) or (u > 0 and math.sqrt(u)) or nil
		v = (isZero(v) and 0) or (v > 0 and math.sqrt(v)) or nil

		if u == nil or v == nil then return 0 end

		h = 1
		j = q < 0 and -v or v
		k = z - u

		num = Polynomial.SolveQuadric(h, j, k, reftable)

		j = -j
		k = z + u

		if num == 0 then
			num = num + Polynomial.SolveQuadric(h, j, k, reftable)
		elseif num == 1 then
			num = num + Polynomial.SolveQuadric(h, j, k, holder)

			reftable[2] = holder[1]
			reftable[3] = holder[2]
		elseif num == 2 then
			num = num + Polynomial.SolveQuadric(h, j, k, holder)

			reftable[3] = holder[1]
			reftable[4] = holder[2]
		end
	end

	if num > 0 then reftable[1] = reftable[1] - sub end
	if num > 1 then reftable[2] = reftable[2] - sub end
	if num > 2 then reftable[3] = reftable[3] - sub end
	if num > 3 then reftable[4] = reftable[4] - sub end

	return num
end

return Polynomial