Talk:Cubic function/Archive 4

Archive 1

Archive 2

Archive 3

Archive 4

Archive 5

Archive 6

→

Archive 8

Needs clear algorithm

I'm writing software that needs to do a cubic solve. In Python, numpy.roots works great, but is slow for solving thousands of them in parallel, so I came here in search of the algorithm... but it isn't at all clear what the full algorithm is. Other web pages fail me too.

Am I just blind, or are none of the algorithms listed completely? In the "General formula of roots" section, it says if $-27a^{2}\Delta >0$, then the solutions are ___. But then it gives caveats and the version with Q and C, but that also has caveats. Is that Q and C version complete? If so, it doesn't read like an algorithm... —Ben FrantzDale (talk) 13:40, 11 October 2011 (UTC)

Do you really need a cubic solver whose output is thousands formulas with cubic and square roots inside? This seems incredible. I guess that you want numerical output. For this it is clearly better to use a root finding algorithm, because computing a numerical approximation of the roots of a cubic equation with a root finding algorithm is usually much faster than evaluating numerically the formula. In fact evaluating the formula implies to compute a square root, a cubic root and several other operations. Cubic root extraction is not really faster than a general root finding algorithm applied to a general cubic equation. D.Lazard (talk) 21:08, 12 October 2011 (UTC)

You may be right about performance. I had sort of assumed that since cubics can be solved in closed form, that it would be faster than an iterative approach. Still, the "algorithms" on the page aren't at all clear. I should be able to use the information on this page to find the roots of a cubic function, but attempts to follow the steps failed, usually in ambiguity. —Ben FrantzDale (talk) 13:08, 13 October 2011 (UTC)

I do not agree with you about the clarity of the "algorithms": The formula involving $C$ and $Q$ looks perfectly like an algorithm. It specifies the auxiliary functions (square and cubic roots). Its main difference with a well written algorithm is that the simple cases are at the end while a program is easier to read if the simple cases are at the beginning. Another difference with a well written algorithm is that in a program, names will be given to the expressions which appear several times, to avoid to compute them several times. For better readability, names have been given only to the quantities whose sign or equality to 0 have to be tested. I do not see any way to improve the algorithmic aspect of this page, while keeping a good readability, especially for beginners. Note that the first general formula is not useful from a logical point of view, and could thus be suppressed. But it has to appear, as it is the formula which is the most popular. D.Lazard (talk) 10:58, 14 October 2011 (UTC)

Some times ago I saw that sympy implement Cardano's formulas. However, I would have to consult the source again to find out if it is in any way stabilized. It is not to much effort to first implement the shift and then use the simple steps of the derivation.--LutzL (talk) 13:35, 13 October 2011 (UTC)

I must have been mistaken. I just re-tried coding the Q, C algorithm and it seems to work up to rounding error. For reference. This implementation solves many 4xn vectors at the same time (hence the lack of if/else and instead the use of masked arrays):

def solveCubicRoots(p):
    p = asarray(p)
    if p.ndim == 1:
        return solveCubicRoots(p[:,None])[:,0]
    p = cast[complex](p)
    a, b, c, d = p
    b2m3ac = b**2 - 3*a*c
    Q = sqrt((2*b**3 - 9*a*b*c+27*a**2*d)**2 - 4 * (b**2-3*a*c)**3)
    def CofQ(Q, m): 
        a, b, c, d = p[:,m]
        return (0.5*(Q+2*b**3-9*a*b*c+27*a**2*d))**(1./3.)

    result = zeros(p[:3].shape, dtype=complex)
    
    # Four cases: zero or not for Q and for b2m3ac:
    ZZ = logical_and(Q == 0, b2m3ac == 0)
    ZN = logical_and(Q == 0, b2m3ac != 0)
    NZ = logical_and(Q != 0, b2m3ac == 0)
    NN = logical_and(Q != 0, b2m3ac != 0)

    a, b, c, d = p[:,ZZ]
    result[:,ZZ] = -b/3*a
        
    a, b, c, d = p[:,ZN]
    result[:2,ZN] = (b*c - 9*a*d)/(2*(3*a*c-b**2))
    result[2,ZN] = (9*a**2*d-4*a*b*c+b**3)/(a*(3*a*c-b**2))
    
    QNZ = (Q != 0)
    C = CofQ(Q[QNZ], QNZ)
    mask = QNZ.copy()
    mask[QNZ] = logical_and(b2m3ac[QNZ] == 0, C == 0)
    C[mask] = CofQ(-Q[mask], mask)
    a, b, c, d = p[:,QNZ]
    i = complex(0,1)
    result[0,QNZ] = (-b - C - (b**2-3*a*c)/C)/(3*a)
    result[1,QNZ] = -b/(3*a) + C*(1+i*sqrt(3))/(6*a) + (1-i*sqrt(3))*(b**2-3*a*c)/(6*a*C)
    result[2,QNZ] = -b/(3*a) + C*(1-i*sqrt(3))/(6*a) + (1+i*sqrt(3))*(b**2-3*a*c)/(6*a*C)
    
    a, b, c, d = p[:,ZZ]
    result[:,ZZ] = -b/(3*a)

    return result

I may rephrase the description of this algorithm to make the description more imperative. —Ben FrantzDale (talk) 14:53, 17 October 2011 (UTC)