[RESOLVED] Solve a cubic polynomial

[Kallog]

Sounds easy, but I can't quite do it reliably on the computer.

I'm using what I think is the Cardini algorithm. It first converts to the depressed form, then chooses one of 3 options depending on if some expression is <0, =0 or >0.

Problem is it's hard to do floating point comparisons safely, so some of my equations are ending up with the wrong number of solutions.

Is there a robust way to do it?

Oh, also I only need real solutions. Complex solutions should become 0. This also worries me a bit. I have a feeling rounding errors might cause a solution to have a very tiny imaginary component, then when I set it to zero I lose the entire real component.

[jemidiah]

I think you mean Cardano instead of Cardini. Without more thought than it's worth, I can't tell you if that method is sensitive to roundoff errors. I can say that the analytic solution to the cubic equation is perhaps of more theoretical interest than practical interest, considering the large number of numerical root finding algorithms.

If I were you, I'd use Newton's Method (a search or the Wikipedia article would give details) to find the real root of your cubic (one is guaranteed to exist since the polynomial is of odd degree). Newton's Method is easy to implement, fast in this case, and should converge to the limit of your machine arithmetic. You could then use polynomial division (equivalently, synthetic division) to get a quadratic which is easily solved using the standard quadratic formula. The benefit to this setup is that your only floating point comparison is testing the discriminant of the final quadratic which is successfully done every day by some programmer. It's pretty robust to rounding errors unless your real root and the discriminant of the resulting quadratic are of wildly different orders of magnitude. You're also guaranteed one "good" solution no matter what.

If you want more details, I'm happy to provide them, though the methods I mentioned are well-known, well-studied, and well-documented. Analytic solutions are fine for theory, but in reality huge swaths of applied math is done using numerical (i.e. approximate) methods.

Edit: In response to your worry about losing valid solutions to rounding errors, if you find it's a problem you could test a/(sqrt(a^2 + b^2)) and see if it's close to 1 (0.9999, say) for complex solution a+bi. This corresponds to b much smaller than a. A similar test, |a/b| > 10000, say, may perform better since the Sqrt(_^2 + _^2) is pretty susceptible to floating point rounding errors in the addition.

[jemidiah]

I got to thinking it would take me much less time to implement the above than for someone who's not familiar with those methods, so I quick implemented it in VB6. Of course it could be cleaned up a fair amount.

Const NearZero = 0.00000000000001

Private Sub Form_Load()
Dim a As Double, b As Double, c As Double, d As Double
a = 1
b = -25
c = 200
d = -500

Dim x As Double
Dim xn As Double
Dim diff As Double
diff = 1
x = 0
Do Until diff < NearZero
    xn = x - f(x, a, b, c, d) / fp(x, a, b, c)
    diff = Abs(xn - x)
    x = xn
Loop

'Found a real root
Dim r As Double
r = x

MsgBox "Root " & r

'Polynomial division gives quadratic
'a x^2 + (b + r a ) x + (c + r b + r^2 a)
Dim u As Double, v As Double, w As Double
u = a
v = b + r * a
w = c + r * b + r * r * a

Dim discriminant As Double
discriminant = v * v - 4 * u * w
If v = 0 Then
    If Abs(discriminant) < NearZero Then
        'Basically zero; suppose it is exactly zero
        MsgBox "Double root: " & -v / (2 * u)
        End
        Exit Sub
    End If
ElseIf Abs(Sqr(Abs(discriminant)) / v) < NearZero Then
    'Basically zero; suppose it is exactly zero
    MsgBox "Double root: " & -v / (2 * u)
    End
    Exit Sub
End If

If discriminant < 0 Then
    MsgBox "Complex roots."
Else
    MsgBox "Two real roots: " & _
           (-v + Sqr(discriminant)) / (2 * u) & ", " & _
           (-v - Sqr(discriminant)) / (2 * u)
End If
End
End Sub

Private Function f(x As Double, a As Double, b As Double, c As Double, d As Double) As Double
    f = a * x * x * x + b * x * x + c * x + d
End Function

Private Function fp(x As Double, a As Double, b As Double, c As Double) As Double
    fp = 3 * a * x * x + 2 * b * x + c
End Function

[Kallog]

Cheers, that's pretty awesome!

But I think I'll do a kind of hybrid method for now, unless I find that doesn't work:

1. Use Cardano to find the 1st root. This includes a > test but values just either side of 0 seem to end up with the same result so it should be OK.

2. Use long division to get a quadratic as you suggested.

3. Use the quadratic equation to find the other 2 roots, as you suggested.


Problem is in step 2. I don't understand the maths. I can do the long division. Then I end up with a quadratic quotient. Do I have to set this to 0 and solve? What happens to the remainder? I'm just not sure how this makes sense

I could copy and paste your code but I'd rather know what's going on. I got into this problem in the first place because of somebody else's code.

[Kallog]

Oh it's OK.

I've realized the remainder has to be zero because we divided the cubic by one of its factors. Also the remainder has the same form as the cubic but with r substituted for x, so it has to be zero again.

I don't understand how you got your quadratic:
'a x^2 + (b + r a ) x + (c + r b + r^2 a)

But the one I got matched Wikipedia so I'll go with that.

[jemidiah]

I got my quadratic through standard polynomial division of ax^3 + bx^2 + cx + d by (x-r). I also checked it for several values and had Mathematica multiply out (a x^2 + (b + r a ) x + (c + r b + r^2 a))*(x-r) to make sure I didn't make a mistake (I got back ax^3 + bx^2 + cx + d - ar^3 - br^2 - cr - d = original, as you found out, since the last four terms are 0). What's the Wikipedia one you found?

[Kallog]

I got my quadratic through standard polynomial division of ax^3 + bx^2 + cx + d by (x-r). I also checked it for several values

Ah oops. I didn't notice you're using different coefficient labels to Wikipedia. They divided them all by a and renamed what was left to a,b,c. I'm doing it with only 2 coefficients because that's the form Cardino has it in.

[jemidiah]

Oh okie. Glad it seems to have worked out .

[Kallog]

Oh okie. Glad it seems to have worked out .

Actually it hasn't anymore :P There's something wrong with my Cardano, but I'll sort that out - might have to use Newton's method for that afterall.

I found out that the equation I'm solving is guaranteed to have only real roots. When they become imaginary it's only because of floating point errors. So rather than checking for small imaginary components, I just force the discriminant to be zero if it happens to come out negative.

[jemidiah]

Newton's Method iteratively finds roots of a curve P(x). It starts at an arbitrary x value, draws the line tangent to P(x) at the current x value, and calculates where that line intersects the x-axis. The intersection point is then used as the new x value, where the process repeats. Drawing it out for, say, a parabola shows you how successive iterations get arbitrarily close to the root.

If you've never had Calculus, the calculation of the slope of the tangent curve using the derivative of P(x) will have to remain magical, but the rest is quite understandable IMO. (The derivative of ax^3 + bx^2 + cx + d is 3ax^2 + 2bx + c, as is coded in my fp routine.)

Say you start at x0. If you do the math, the line tangent to P(x) at x0 happensto insersect the x-axis at

x1 = x0 - P(x0) / P'(x0)

where P'(x0) is the slope of P(x) at x0. This leads to an easy iteration scheme where you keep iterating until you get no appreciable change in the above. (At a root, P(x0) = 0, so x1 = x0, and the value doesn't change during the iteration.)

[Kallog]

Yea I know Newton's method.

I think I've got it fine now using another formula off Wikipedia which finds the/a real root.

I hadn't realised you can't do the cube root of a negative number tho!!! I can't believe I've been programming for like 15 years and never worked that out. So I wrote a little CubeRoot() function that returns the negative of the positive square root if the input is negative.

Here's the final result. Only works where all roots are real, or close to. And it has no arbitrary numbers, which I like.

        'a,b,c are the coefficients of x^2 and x and the constant.

        'Find one root using the general formula
        Dim m As Double
        Dim k As Double
        m = 2 * a ^ 3 - 9 * a * b + 27 * c
        k = a ^ 2 - 3 * b
        discriminant = Math.Max(m ^ 2 - 4 * k ^ 3, 0)
        x1 = -a / 3 _
           - 1 / 3 * CubeRoot(0.5 * (m + Math.Sqrt(discriminant))) _
           - 1 / 3 * CubeRoot(0.5 * (m - Math.Sqrt(discriminant)))

        'The remaining 2 roots are found by solving with the quadratic formula.
        discriminant = Math.Max(-3 * x1 ^ 2 - 2 * a * x1 + a ^ 2 - 4 * b, 0)
        x2 = (-a - x1 + Math.Sqrt(discriminant)) / 2
        x3 = (-a - x1 - Math.Sqrt(discriminant)) / 2

[jemidiah]

You can actually take the cube root of a negative number. (For instance, (-2)^3 = -8, so CubeRoot(-8) = -2.) However, on a computer, you might have trouble since any slight rounding error in the exponent of an expression like (-x)^(1/3) will cause the result to be imaginary. Maybe that's what you meant.

[Kallog]

computer, you might have trouble since any slight rounding error in the exponent of an expression like (-x)^(1/3) will cause the result to be imaginary. Maybe that's what you meant.

Yep exactly that. (-5)^(1/3) results in NaN. At least in VB.Net, but I think most languages are the same.

[Kallog]

Wow, this problem's getting harder and harder. Turns out that code I posted doens't actually work except in some special cases. Now I've gone back to Cardano, which is more right, but is sometimes a bit unstable.

Newton's method also has plenty of problems. I suppose it'd be fine if you had clever convergence checks and got it to try again with different start points when it failed.

I can't believe there isn't a safe way to solve a cubic with real coefficients and real solutions.

	'Solve x^3 + a x^2 + b x + c = 0

        'Find one root using Cardano's method
        Dim p As Double
        Dim q As Double
        p = b - a ^ 2 / 3
        q = c + (2 * a ^ 3 - 9 * a * b) / 27
        discriminant = q ^ 2 / 4 + p ^ 3 / 27
       
        'Solve   u = CubeRoot(-q / 2 + Math.Sqrt(discriminant))
        'and     x1 = u - p / (3 * u) - a / 3
        'with complex u
        'Any of the 3 solutions will do, they're all real.
        '
        'Find one of the 3 complex cube roots, u, using De Moivre:
        Dim Re_u As Double
        Dim Im_u As Double
        Dim Modu As Double
        Dim Argu As Double
        Modu = CubeRoot(Math.Sqrt((-q / 2) ^ 2 + Math.Abs(discriminant)))
        Argu = Math.Atan2(Math.Sqrt(Math.Abs(discriminant)), -q / 2) / 3
        Re_u = Modu * Math.Cos(Argu)
        Im_u = Modu * Math.Sin(Argu)
        '
        'Evaluate x1. Ignore the imaginary component because all 3 solutions must be real
        If Re_u = 0 And Im_u = 0 Then
            'Avoid division by zero. This seems to be the limit as Re_u and Im_u tend to zero.
            x1 = -a / 3
        Else
            x1 = Re_u - p * Re_u / (3 * Re_u ^ 2 + 3 * Im_u ^ 2) - a / 3
        End If


        'The remaining 2 roots are found by solving with the quadratic formula.
        discriminant = Math.Max(-3 * x1 ^ 2 - 2 * a * x1 + a ^ 2 - 4 * b, 0)
        x2 = (-a - x1 + Math.Sqrt(discriminant)) / 2
        x3 = (-a - x1 - Math.Sqrt(discriminant)) / 2

[jemidiah]

I'd love to see what cases you're feeding in to Newton's Method that's killing it on this problem. There are a few edge cases to watch out for (you accidentally iterate to a point with zero derivative, or very small derivative, for instance; you start at 1 when the root is at 10^1234, etc.) but I don't see them as difficult to overcome most of the time in the special case of a cubic. Perhaps I'm being naive. Similarly, what cases are causing your Cardano version to fail?

The first condition in the following worries me slightly:

        If Re_u = 0 And Im_u = 0 Then
            'Avoid division by zero. This seems to be the limit as Re_u and Im_u tend to zero.
            x1 = -a / 3
        Else
            x1 = Re_u - p * Re_u / (3 * Re_u ^ 2 + 3 * Im_u ^ 2) - a / 3
        End If

It's not the limit as Re_u and Im_u tend to zero [that limit is infinity]. I'm not sure it would give correct results ever, actually, but I suppose you tested it?

[Kallog]

I'd love to see what cases you're feeding in to Newton's Method that's killing it on this problem. There are a few edge cases to watch out for (you accidentally iterate to a point with zero derivative, or

I havn't implemented it yet. But yea, I don't want it to ever fail. The whole reason I'm doing this is because the original code was fine most of the time, and broke in a few special cases.

It's not the limit as Re_u and Im_u tend to zero [that limit is infinity]. I'm not sure it would give correct results ever, actually, but I suppose you tested it?

That's turned out to be the main problem now. I know it looks like infinity but it can't be because polynomials can't have infinite roots. I sort of got it by trial and error so I can't be sure if it's always correct. However it is correct for the input:

a=-3, b=3, c=-1

which is three repeated real roots at a point of inflection. What I'm worried about is other cases where the limit perhaps isn't -a/3. Also about having near-zero values that aren't picked up by the zero-test and cause division by zero.

I suppose basically I want to find failing cases because I don't trust this code.

[jemidiah]

I got to thinking that quadrature should work in this case, since you just want the single real root of the cubic (and you seem happy with the polynomial division->solve quadratic solution for the other roots). Quadrature takes an interval, (a, b) [representing real numbers x with a < x < b] where you know your root is in the interval. If your function is monotonic (always increasing or always decreasing) on the interval, you can cut the interval in half by just checking the value of the function at the midpoint and seeing if the root must lie between a and the midpoint or between the midpoint and b.

In this case, the major complication is that a cubic equation is often not monotonic. However, it *is* monotonic on certain intervals. If the equation has two local extrema, it is monotonic between them. Also, it is always monotonic between a local extrema and either +inf or -inf (depending on if there is another extrema between the one in question and the relevant infinity). If there are no extremas, the function is monotonic everywhere.

So, there are either 1 or 3 monotonic intervals for any cubic: (-inf, inf) [if the cubic is entirely monotonic]; or (-inf, d1) [if d1 is the leftmost local extrema], (d1, d2) [if d2 is the rightmost extrema], and (d2, inf). For quadrature to work, whatever monotonic interval we pick must contain a root. In the case of the interval (-inf, inf), we can't plug in "infinity" to a computer, so we need to search for an interval with a root. Any interval will work, though, since the polynomial is always monotonic. In this case, you can start "expanding" any arbitrary interval until the signs of the function at the endpoints differ, at which point you have a monotonic interval containing a root.

In the case of the other intervals, a root may lie in more than one interval. If the signs of the function at the endpoints of (d1, d2) differ, a root lies in the interval and, since the interval is monotonic, we're ready for quadrature. Otherwise, if both d1 and d2 are positive, a root cannot lie between them. However, we do know a root lies between (-inf, d1), since the polynomial goes to -inf as x-> -inf. We can find an actual left endpoint by "expanding left" pretty much like we did before. If both d1 and d2 are negative, the root lies between (d2, +inf), and a right endpoint can be found by "expanding right". In either case, we again have an interval on which the polynomial is monotonic, which contains a root.

Say the interval obtained from the above is (l, r). Quadrature just cuts this interval in half at each iteration, narrowing the size of the range in which the root lies exponentially quickly, so it converges very quickly. I implemented this in VB6 and got the following (correct) test results using cubic equation x^3 + ax^2 + bx + c:

a=96, b=15*10^15, c=14
=> r = -9.33...3 * 10^-16 (correct according to Mathematica's NSolve)

a=10, b=0, c=5
=> r = -10.0495085668579 (correct)

a=-10, b=0, c=-5
=> r = 10.0495085668579 (correct)

The nice thing about quadrature here is that it's pretty robust--no division by variables, anywhere, just >, <, = 0 tests.

The code is below. One thing I thought of but didn't have a convenient mechanism for was that each time an interval is expanded, the endpoints should be checked to see if they're zero, at which point the root should be returned straight off. For polynomials "in the wild" this should almost never be useful, but it is a valid edge case.

Option Explicit

Const NearZero = 0.00000000000001

Dim a As Double, b As Double, c As Double

Private Sub Form_Load()
'Find real root of
'   f(x) = x^3 + a x^2 + b x + d
'using quadrature. Since the leading
'term is x^3, the polynomial goes from
'positive to negative as it goes from
'-inf to +inf.
a = -10
b = 0
c = -5

'Find extremums, d1 and d2,
'so we can find a monotonic interval
'(l, r) containing a root.
Dim d1 As Double, d2 As Double
Dim fd1 As Double, fd2 As Double
Dim l As Double, r As Double
l = -1
r = 1

Dim discriminant As Double
discriminant = a * a - 3 * b
If discriminant <= 0 Then
    'd1 and d2 imaginary =>
    '    f(x) monotonic for real x
    'd1 = d2 =>
    '    f'(x) has double root =>
    '    f(x) monotonic
    'so f(x) monotonic here
    ExpandRootInterval l, r, True, True
Else
    d1 = -a / 3 - 1 / 3 * Sqr(discriminant)
    d2 = -a / 3 + 1 / 3 * Sqr(discriminant)
    fd1 = f(d1)
    fd2 = f(d2)
    
    If fd1 = 0 Then
        MsgBox "Got lucky; root = " & d1
        Exit Sub
    ElseIf fd2 = 0 Then
        MsgBox "Got lucky; root = " & d2
        Exit Sub
    End If
    
    If SignDiffers(fd1, fd2) Then
       'Root is between d1 and d2, since
       'sign changes. This interval
       'is also conveniently monotonic.
        l = fd1
        r = fd2
    ElseIf fd1 > 0 And fd2 > 0 Then
        'Root is "left" of d1, since
        'd1 is positive but function
        'eventually gets negative as we go left
        l = d1 - 1
        r = d1
        ExpandRootInterval l, r, True, False
    ElseIf fd1 < 0 And fd2 < 0 Then
        'Root is "right" of d2, since
        'd2 is negative but function
        'eventually gets positive as we go right
        l = d2
        r = d2 + 1
        ExpandRootInterval l, r, False, True
    End If
End If

'At this point, (l, r) contains
'a real root, and f is monotonic on this interval.
'So, quadrature can be applied to find the root.
Dim fl As Double, fr As Double
Dim m As Double 'midpoint of interval
Dim fm As Double 'f at midpoint of interval

Dim intLen As Double, dLen As Double
intLen = r - l
dLen = 1
While dLen > 0
    'At each step the interval length
    '    is halved, so this loop ends quickly.
    'Stop when the interval stops changing in length
    '    (as far as the machine can tell)

    m = (r + l) / 2
    fm = f(m)
    
    fl = f(l)
    fr = f(r)
    
    If SignDiffers(fl, fm) Then
        'Root between l and m
        r = m
    Else
        'Root between m and r,
        'since fl and fr always differ in sign
        l = m
    End If
    
    dLen = intLen - (r - l)
    intLen = r - l
Wend

MsgBox "Root is in interval [" & l & ", " & r & "]"
End
End Sub

Private Function f(x As Double) As Double
f = x * x * x + a * x * x + b * x + c
End Function

Private Function SignDiffers(a As Double, b As Double) As Boolean
SignDiffers = (a < 0 And b > 0) Or (a > 0 And b < 0)
End Function

Private Sub ExpandRootInterval(ByRef l As Double, ByRef r As Double, _
                               ExpandLeft As Boolean, ExpandRight As Double)
'Takes in interval (l, r), l < r.
'At each step, doubles the distance between
'   l and midpoint (if ExpandLeft) and between
'   midpoint and r (if ExpandRight)
'Stops when f(l) and f(r)'s signs differ
'Reaches double size maximum relatively quickly.
'Returns values in input variables l and r, which
'   are passed ByRef
Dim m As Double
Do Until SignDiffers(f(l), f(r))
    m = (r + l) / 2
    If ExpandLeft Then
        l = m - 2 * (m - l)
    End If
    If ExpandRight Then
        r = m + 2 * (r - m)
    End If
Loop
End Sub

[Kallog]

Getting a bit more confidence in that limit.

From the expression for u, it looks like u can only be zero if p is zero.
If p is zero then the expression for x1 could potentially have a limit of -a/3, depending on how fast p goes to zero.

[Kallog]

I got to thinking that quadrature should work in this case, since you just want the single real root of the cubic (and you seem

Hmm, that's quite an intuitive way. And you reckon it's 100% bulletproof? With all those different cases it considers, I wonder if you could make Newton bulletproof by treating each trouble case separately too.

By the way, speed is important, tho I see this seems pretty fast, I have to do half a million of these solves, ideally without the user waiting around. I have a suspicion Newton's method would take less iterations to get to double precision accuracy.

I might just stick with the last code I posted, and lots of test cases. Until a few weeks later when I discover another horrible fault again and come back, maybe then I'll have to bite the bullet and do it iteratively.

[jemidiah]

The method I posted would require a maximum of... log_2(1.79769E+308) = 1024 iterations to find the range (l, r) for use in quadrature. The quadrature itself would take a maximum of twice that many iterations to narrow the range down to 1 unit, followed by another ~50 to get 15 digits of precision. These are wildly extreme cases--your root would need to be extravagantly huge to take that many iterations. For instance, in the large test case (a=96, b=15*10^15, c=14) this method took 0 and then 104 iterations finding the range and then doing the quadrature, respectively. The case a=10, b=0, c=5 took 4 and 52 iterations.

Repeating the a=10, b=0, c=5 case 500,000 times took 19 seconds on my machine running VB6 code. It's not instant, but solving 500,000 cubic equations in well under a minute seems very acceptable to me without knowing your application. .NET would almost certainly be significantly faster. Applying an optimization...

Dim intLen As Double, dLen As Double
intLen = r - l
dLen = 1
fl = f(l)
fr = f(r)
While dLen > 0
    'At each step the interval length
    '    is halved, so this loop ends quickly.
    'Stop when the interval stops changing in length
    '    (as far as the machine can tell)

    m = (r + l) / 2
    fm = f(m)
    If SignDiffers(fl, fm) Then
        'Root between l and m
        r = m
        fr = fm
    Else
        'Root between m and r,
        'since fl and fr always differ in sign
        l = m
        fl = fm
    End If
    
    dLen = intLen - (r - l)
    intLen = r - l
Wend

...cut the run time to under 12 seconds, so I'm certain you could squeeze some more out if you really wanted, not to mention the method is completely parallelizable. If translation to .NET doubled the speed and all 4 of my cores were used, the run time would be about 1.5 seconds. Another optimization, effectively inlining SignDiffers, cut my runtime down to about 9.5 seconds--so we'd be down to 1.2 seconds in theoretical .NET on 4 cores.

And you reckon it's 100&#37; bulletproof? With all those different cases it considers, I wonder if you could make Newton bulletproof by treating each trouble case separately too.

This is far more bulletproof than Newton's Method. In retrospect, I hadn't considered that complex roots could attract Newton's Method if those roots are much closer to your guess than the real root is. These cases are hard to find ahead of time without computing the roots in the first place. Since all of your roots are real this might not be a problem, but it depends--if some rounding somewhere has made your roots slip into being barely complex, Newton's Method might well not converge. This will always converge, every time, within the above time limits.

Compared to Cardano's method, this again has the advantage that there's no division by possibly small (or 0) numbers. Everything is done with multiplication and is quite well-behaved. While there are a number of cases, they're all relatively simple, requiring pretty simple operations instead of numerically unstable algebra [the discriminant test at the start requires no division].

I might just stick with the last code I posted, and lots of test cases. Until a few weeks later when I discover another horrible fault again and come back, maybe then I'll have to bite the bullet and do it iteratively.

That sounds like a dangerous plan considering the number of times you've revised your strategy already. Quadrature seems like the silver bullet you've been looking for.

From the expression for u, it looks like u can only be zero if p is zero.

Yup. If p is zero, the depressed cubic has the simple form

t^3 + q = 0

so one solution is t = CubeRoot(-q). Substituting back for x, since x = t - a/3, x = CubeRoot(-q) - a/3, which is quite close to what you had. In the case you gave, a=-3, b=3, c=-1, q is given by

q = c + (2a^3 - 9ab)/27
= -1 + (-2*27 + 27*3) / 27
= -1 + (3-2) = 0

so CubeRoot(-q) = 0 and x reduces to -a/3 as you said. However, this is pure dumb luck. Considering limits isn't valid since p is assumed non-zero after that substitution (presumably with the assumption that p=0 can be easily special-cased like I did above).

[Kallog]

Yup. If p is zero, the depressed cubic has the simple form

t^3 + q = 0

so one solution is t = CubeRoot(-q). Substituting back for x, since x = t - a/3, x = CubeRoot(-q) - a/3, which is quite close to what you had.

OK I can rephrase the comments in my code to explain this case as a triple root, rather than a limit.

Although assuming p=0 leads to x=-a/3 might have been dumb luck, I think it's still valid luck -

When p=0:
If q is nonzero then there are 2 complex values of t=AllCubeRoots(-q)
If there are any complex values of t, then there are corresponding complex values of x=t-a/3
But I know my problem can't have any complex values of x, so q can't be non-zero and x=-a/3
If there's a slight rounding error that causes complex x's this must be in the form of q being slightly non-zero.
If q is slightly non-zero then x=CubeRoot(-q) - a/3 can be replaced with x=-a/3 to correct the error
Bazinga!


That only leaves the risk of numerical instability when a not-quite-zero denominator isn't caught by the = test. But I've tried it with cases that are so close to zero it sometimes gets caught in the test and sometimes not - either way the final results come out pretty much the same. Not amazing accuracy but good enough. If need be I can replace the = test with a < 1e-13 or so test.

[jemidiah]

...Bazinga!

Good point, you're right that in your case dropping the CubeRoot(-q) is in fact justified. I meant that the conclusion was dumb luck, not backed up by any solid reasoning, and was potentially susceptible to errors. But now it is backed up by solid reasoning, so yay .

[Kallog]

Good point, you're right that in your case dropping the CubeRoot(-q) is in fact justified. I meant that the conclusion was dumb luck, not backed up by any solid reasoning, and was potentially susceptible to errors. But now it is backed up by solid reasoning, so yay .

Cheers for your help. Let's hope I don't have to come crawling back haha.