def cubic(a, b, c, d=None):
if d: <EOL><INDENT>a, b, c = b / float(a), c / float(a), d / float(a)<EOL><DEDENT>t = a / <NUM_LIT><EOL>p, q = b - <NUM_LIT:3> * t**<NUM_LIT:2>, c - b * t + <NUM_LIT:2> * t**<NUM_LIT:3><EOL>u, v = quadratic(q, -(p/<NUM_LIT>)**<NUM_LIT:3>)<EOL>if type(u) == type(<NUM_LIT>): <EOL><INDENT>r, w = polar(u.real, u.imag)<EOL>y1 = <NUM_LIT:2> * cbrt(r) * cos(w / <NUM_LIT>)<EOL><DEDENT>else: <EOL><INDENT>y1 = cbrt(u) + cbrt(v)<EOL><DEDENT>y2, y3 = quadratic(y1, p + y1**<NUM_LIT:2>)<EOL>return y1 - t, y2 - t, y3 - t<EOL>
x^3 + ax^2 + bx + c = 0  (or ax^3 + bx^2 + cx + d = 0)
    With substitution x = y-t and t = a/3, the cubic equation reduces to
        y^3 + py + q = 0,
    where p = b-3t^2 and q = c-bt+2t^3.  Then, one real root y1 = u+v can
    be determined by solving
        w^2 + qw - (p/3)^3 = 0
    where w = u^3, v^3.  From Vieta's theorem,
        y1 + y2 + y3 = 0
        y1 y2 + y1 y3 + y2 y3 = p
        y1 y2 y3 = -q,
    the other two (real or complex) roots can be obtained by solving
        y^2 + (y1)y + (p+y1^2) = 0