--2nd and 3rd degree equation solvers --Written by Cosmin Apreutesei. Public domain. local abs, sqrt, acos, cos, pi = math.abs, math.sqrt, math.acos, math.cos, math.pi --2nd degree equation solver, see http://en.wikipedia.org/wiki/Quadratic_equation#Floating-point_implementation. --epsilon controls the precision of the solver. local function solve2(a, b, c, epsilon) epsilon = epsilon or 1e-30 if abs(a) < epsilon then if abs(b) < epsilon then return end --contradiction return -c / b --1st degree end local D = b^2 - 4*a*c --discriminant if D > epsilon then --D > 0, real root local q = -1/2 * (b + (b >= 0 and 1 or -1) * sqrt(D)) return q / a, c / q elseif D < -epsilon then --D < 0, complex root return else --D == 0, double root return -b / (2 * a) end end --3rd degree equation solver using the trigonometric method, see http://en.wikipedia.org/wiki/Cubic_function. --epsilon controls the precision of the solver. --TODO: review this code for catastrophic cancelations and alike. local function solve3(a, b, c, d, epsilon) epsilon = epsilon or 1e-16 if abs(a) < epsilon then --2nd degree return solve2(b, c, d) end local z = a a = b / z b = c / z c = d / z local p = b - a^2 / 3 local q = a * (2 * a^2 - 9 * b) / 27 + c local p3 = p^3 local D = q^2 + 4 * p3 / 27 local offset = -a / 3 if D > epsilon then z = sqrt(D) local u = ( -q + z) / 2 local v = ( -q - z) / 2 u = u >= 0 and u^(1 / 3) or -((-u)^(1 / 3)) v = v >= 0 and v^(1 / 3) or -((-v)^(1 / 3)) return u + v + offset elseif D < -epsilon then local u = 2 * sqrt( -p / 3) local v = acos( -sqrt( -27 / p3) * q / 2) / 3 return u * cos(v) + offset, u * cos(v + 2 * pi / 3) + offset, u * cos(v + 4 * pi / 3) + offset else --D == 0 local u = q < 0 and (-q / 2)^(1 / 3) or -((q / 2)^(1 / 3)) if abs(u) < epsilon then return offset end return 2*u + offset, -u + offset end end if not ... then require'eq_test' end return { solve2 = solve2, solve3 = solve3, }