C++: Solving Cubic Equations

Question

This site solves cubic equations, and has the equations it uses

I wrote this function to get the same results but it's not working

void cubicFormula(float a, float b, float c, float d, float *results)
{
    float a2 = a*a;
    float b2 = b*b;
    float b3 = b2*b;

    float q = (3*c- b2)/9;
    float q2 = q*q;
    float q3 = q2*q;
    float r = -27*d + b*(9*c-2*b2);
    float r2 = r*r;

    float discriminant = q3 + r2;

    float s = r + sqrtf(discriminant);
    float t = r - sqrtf(discriminant);

    float term1 = powf((-t + s)/2, 1/3.);

    float r13= 2 * sqrtf(q);

    results[0] = (- term1 + r13 * cosf(q3/3));
    results[1] = (- term1 + r13 * cosf(q3 + (2*M_PI)/3));
    results[2] = (- term1 + r13 * cosf(q3 + (4*M_PI)/3));
}

What could be the problem, is there something I am missing?

Update:

For instance, using these values

a = -1.000000;
b = 36.719475;
c = -334.239960;
d = 629.877808;

The discriminant is less then zero, so r13 is nan

But the site mentioned above and a couple of online solvers I've tested find the three roots

23.775687816485593
2.548516232734247
10.395270950780164

All of them real numbers

Answer 1

The site you link to provide a Cubic Equation Calculator. You could start by looking for differences between your code and theirs:

function cubicsolve(dataForm)
{
    var a = parseFloat(dataForm.aIn.value);
    var b = parseFloat(dataForm.bIn.value);
    var c = parseFloat(dataForm.cIn.value);
    var d = parseFloat(dataForm.dIn.value);
    if (a == 0)
    {
        alert("The coefficient of the cube of x is 0. Please use the utility for a SECOND degree quadratic. No further action taken.");
        return;
    } //End if a == 0

    if (d == 0)
    {
        alert("One root is 0. Now divide through by x and use the utility for a SECOND degree quadratic to solve the resulting equation for the other two roots. No further action taken.");
        return;
    } //End if d == 0
    b /= a;
    c /= a;
    d /= a;
    var disc, q, r, dum1, s, t, term1, r13;
    q = (3.0*c - (b*b))/9.0;
    r = -(27.0*d) + b*(9.0*c - 2.0*(b*b));
    r /= 54.0;
    disc = q*q*q + r*r;
    dataForm.x1Im.value = 0; //The first root is always real.
    term1 = (b/3.0);
    if (disc > 0) { // one root real, two are complex
        s = r + Math.sqrt(disc);
        s = ((s < 0) ? -Math.pow(-s, (1.0/3.0)) : Math.pow(s, (1.0/3.0)));
        t = r - Math.sqrt(disc);
        t = ((t < 0) ? -Math.pow(-t, (1.0/3.0)) : Math.pow(t, (1.0/3.0)));
        dataForm.x1Re.value = -term1 + s + t;
        term1 += (s + t)/2.0;
        dataForm.x3Re.value = dataForm.x2Re.value = -term1;
        term1 = Math.sqrt(3.0)*(-t + s)/2;
        dataForm.x2Im.value = term1;
        dataForm.x3Im.value = -term1;
        return;
    } 
    // End if (disc > 0)
    // The remaining options are all real
    dataForm.x3Im.value = dataForm.x2Im.value = 0;
    if (disc == 0){ // All roots real, at least two are equal.
        r13 = ((r < 0) ? -Math.pow(-r,(1.0/3.0)) : Math.pow(r,(1.0/3.0)));
        dataForm.x1Re.value = -term1 + 2.0*r13;
        dataForm.x3Re.value = dataForm.x2Re.value = -(r13 + term1);
        return;
    } // End if (disc == 0)
    // Only option left is that all roots are real and unequal (to get here, q < 0)
    q = -q;
    dum1 = q*q*q;
    dum1 = Math.acos(r/Math.sqrt(dum1));
    r13 = 2.0*Math.sqrt(q);
    dataForm.x1Re.value = -term1 + r13*Math.cos(dum1/3.0);
    dataForm.x2Re.value = -term1 + r13*Math.cos((dum1 + 2.0*Math.PI)/3.0);
    dataForm.x3Re.value = -term1 + r13*Math.cos((dum1 + 4.0*Math.PI)/3.0);
    return;
}  //End of cubicSolve