cubicRoots.lua

local pi    = math.pi
local abs   = math.abs
local atan2 = math.atan2
local cos   = math.cos
local sqrt  = math.sqrt
local sign  = math.sign

local EPS = 1e-9

local function cbrt(x)
	return sign(x)*abs(x)^(1/3)
end

-- solves for x where ax^3 + bx^2 + cx + d = 0
local function cubicRoots(a, b, c, d)
	local r0, r1, r2

	if abs(a) > EPS then
		b = b/a
		c = c/a
		d = d/a
		
		local u = b*b*(b*d - c*c/4)/27 + c*(c*c/27 - b*d/6) + d*d/4
		local v = b*(c/6 - b*b/27) - d/2

		local off = b/3

		if u > 0 then
			local su = sqrt(u)
			local i = cbrt(v + su)
			local j = cbrt(v - su)

			r0 = i + j - off

		elseif u < 0 then
			local scale = 2/3*sqrt(b*b - 3*c)
			local angle = atan2(sqrt(abs(u)), v)/3

			r0 = scale*cos(angle) - off
			r1 = scale*cos(angle + 2/3*pi) - off
			r2 = scale*cos(angle + 4/3*pi) - off

		else
			local cv = cbrt(v)
			
			r0 = 2*cv - off
			r1 = cv - off
		end

	elseif abs(b) > EPS then
		local det = c*c - 4*b*d
		if det >= 0 then
			local s = c + sqrt(det)*(c >= 0 and 1 or -1)
			r0 = s/(-2*b)
			r1 = (-2*d)/s
		end

	elseif abs(c) > EPS then
		r0 = -d/c
	end

	return r0, r1, r2
end