Solving a cubic equation

Question

As part of a program I'm writing, I need to solve a cubic equation exactly (rather than using a numerical root finder):

a*x**3 + b*x**2 + c*x + d = 0.

I'm trying to use the equations from here. However, consider the following code (this is Python but it's pretty generic code):

a =  1.0
b =  0.0
c =  0.2 - 1.0
d = -0.7 * 0.2

q = (3*a*c - b**2) / (9 * a**2)
r = (9*a*b*c - 27*a**2*d - 2*b**3) / (54*a**3)

print "q = ",q
print "r = ",r

delta = q**3 + r**2

print "delta = ",delta

# here delta is less than zero so we use the second set of equations from the article:

rho = (-q**3)**0.5

# For x1 the imaginary part is unimportant since it cancels out
s_real = rho**(1./3.)
t_real = rho**(1./3.)

print "s [real] = ",s_real
print "t [real] = ",t_real

x1 = s_real + t_real - b / (3. * a)

print "x1 = ", x1

print "should be zero: ",a*x1**3+b*x1**2+c*x1+d

But the output is:

q =  -0.266666666667
r =  0.07
delta =  -0.014062962963
s [real] =  0.516397779494
t [real] =  0.516397779494
x1 =  1.03279555899
should be zero:  0.135412149064

so the output is not zero, and so x1 isn't actually a solution. Is there a mistake in the Wikipedia article?

ps: I know that numpy.roots will solve this kind of equation but I need to do this for millions of equations and so I need to implement this to work on arrays of coefficients.

I had some working code that solved for the real roots in VB.NET and I know it worked... I'll try to see if I can find it amongst my folders and post it tomarrow (no net access at home). — J. Polfer, Dec 01 '09 at 22:23
For the x1 solution in the wikipedia article, you add s and t. s has imaginary part theta/3, and t has imaginary part -theta/3, so for that solution the imaginary parts should cancel exactly. — astrofrog, Dec 01 '09 at 22:24
Can you be sure that, for those coefficients, the imaginary part is unimportant? After all `delta` is not zero. — pavium, Dec 01 '09 at 22:25
I'm not saying that the imaginary part is unimportant, I'm saying that s + t is real whatever the imaginary parts for s and t (in the delta < 0 case). — astrofrog, Dec 01 '09 at 22:30
One comment is confusing/wrong : # here delta is zero <- should be "less than zero" — Mark Byers, Dec 01 '09 at 22:32

score 24 · Accepted Answer · answered Dec 01 '09 at 22:44

24

Wikipedia's notation (rho^(1/3), theta/3) does not mean that rho^(1/3) is the real part and theta/3 is the imaginary part. Rather, this is in polar coordinates. Thus, if you want the real part, you would take rho^(1/3) * cos(theta/3).

I made these changes to your code and it worked for me:

theta = arccos(r/rho)
s_real = rho**(1./3.) * cos( theta/3)
t_real = rho**(1./3.) * cos(-theta/3)

(Of course, s_real = t_real here because cos is even.)

answered Dec 01 '09 at 22:44

A. Rex

31,633
21
89
96

There it is. Can't believe I read through that several times. – John Dec 01 '09 at 22:52
Thanks so much for spotting this! I was going crazy trying to understand why it wouldn't work. – astrofrog Dec 01 '09 at 23:01

score 2 · Answer 2 · answered Dec 01 '09 at 22:37

2

I've looked at the Wikipedia article and your program.

I also solved the equation using Wolfram Alpha and the results there don't match what you get.

I'd just go through your program at each step, use a lot of print statements, and get each intermediate result. Then go through with a calculator and do it yourself.

I can't find what's happening, but where your hand calculations and the program diverge is a good place to look.

answered Dec 01 '09 at 22:37

John

15,990
10
70
110

I did go through very carefully, and also tried using a calculator, and I'm pretty confident I'm not doing anything wrong. Could this be a problem with the wikipedia article? – astrofrog Dec 01 '09 at 22:45
See A. Rex's answer. As they say, "Well, THERE'S your problem!" – John Dec 01 '09 at 22:52

score 1 · Answer 3 · answered Oct 25 '22 at 18:03

from cmath import *

solutions = set()


def cbrt(polynomial):
    solution = set()
    root1 = polynomial ** (1 / 3)
    root2 = (polynomial ** (1 / 3)) * (-1 / 2 + (sqrt(3) * 1j) / 2)
    root3 = (polynomial ** (1 / 3)) * (-1 / 2 - (sqrt(3) * 1j) / 2)
    solution.update({root1, root2, root3})
    return solution


def cardano(a, b, c, d):
    p = (3 * a * c - b ** 2) / (3 * a ** 2)
    q = (2 * b ** 3 - 9 * a * b * c + 27 * a ** 2 * d) / (27 * a ** 3)
    alpha = cbrt(-q / 2 + sqrt((q / 2) ** 2 + (p / 3) ** 3))
    beta = cbrt(-q / 2 - sqrt((q / 2) ** 2 + (p / 3) ** 3))
    for i in alpha:
        for j in beta:
            if abs((i * j) + p / 3) <= 0.00001:
                x = i + j - b / (3 * a)
                solutions.add(x)


def quadratic(a, b, c):
    D = b ** 2 - 4 * a * c
    x1 = (-b + sqrt(D)) / 2 * a
    x2 = (-b - sqrt(D)) / 2 * a
    solutions.update({x1, x2})


def linear(a, b):
    if a == 0 and b == 0:
        solutions.add("True")

    if a == 0 and b != 0:
        solutions.add("False")

    if a != 0:
        solutions.add(-b / a)


print('ax^3+bx^2+cx+d=0')
a, b, c, d = map(complex, input().split())

if a != 0:
    cardano(a, b, c, d)

elif b != 0:
    quadratic(b, c, d)

else:
    linear(c, d)

print(solutions)

score 0 · Answer 4 · answered Jun 01 '13 at 00:55

Here's A. Rex's solution in JavaScript:

a =  1.0;
b =  0.0;
c =  0.2 - 1.0;
d = -0.7 * 0.2;

q = (3*a*c - Math.pow(b, 2)) / (9 * Math.pow(a, 2));
r = (9*a*b*c - 27*Math.pow(a, 2)*d - 2*Math.pow(b, 3)) / (54*Math.pow(a, 3));
console.log("q = "+q);
console.log("r = "+r);

delta = Math.pow(q, 3) + Math.pow(r, 2);
console.log("delta = "+delta);

// here delta is less than zero so we use the second set of equations from the article:
rho = Math.pow((-Math.pow(q, 3)), 0.5);
theta = Math.acos(r/rho);

// For x1 the imaginary part is unimportant since it cancels out
s_real = Math.pow(rho, (1./3.)) * Math.cos( theta/3);
t_real = Math.pow(rho, (1./3.)) * Math.cos(-theta/3);

console.log("s [real] = "+s_real);
console.log("t [real] = "+t_real);

x1 = s_real + t_real - b / (3. * a);

console.log("x1 = "+x1);
console.log("should be zero: "+(a*Math.pow(x1, 3)+b*Math.pow(x1, 2)+c*x1+d));

Seeing as you use Math.pow a lot, you might consider this extension: `Number.prototype.pow = function(a) { return Math.pow(this, a); }; 2..pow(3) // 8` — P i, Jul 20 '15 at 13:09

Seonghoon Jang · Answer 5 · 2016-11-12T07:46:03.473

Here, I put a cubic equation (with complex coefficients) solver.

#include <string>
#include <fstream>
#include <iostream>
#include <cstdlib>

using namespace std;

#define PI 3.141592

long double complex_multiply_r(long double xr, long double xi, long double yr, long double yi) {
    return (xr * yr - xi * yi);
}

long double complex_multiply_i(long double xr, long double xi, long double yr, long double yi) {
    return (xr * yi + xi * yr);
}

long double complex_triple_multiply_r(long double xr, long double xi, long double yr, long double yi, long double zr, long double zi) {
    return (xr * yr * zr - xi * yi * zr - xr * yi * zi - xi * yr * zi);
}

long double complex_triple_multiply_i(long double xr, long double xi, long double yr, long double yi, long double zr, long double zi) {
    return (xr * yr * zi - xi * yi * zi + xr * yi * zr + xi * yr * zr);
}

long double complex_quadraple_multiply_r(long double xr, long double xi, long double yr, long double yi, long double zr, long double zi, long double wr, long double wi) {
    long double z1r, z1i, z2r, z2i;    
    z1r = complex_multiply_r(xr, xi, yr, yi);
    z1i = complex_multiply_i(xr, xi, yr, yi);
    z2r = complex_multiply_r(zr, zi, wr, wi);
    z2i = complex_multiply_i(zr, zi, wr, wi);
    return (complex_multiply_r(z1r, z1i, z2r, z2i));
}

long double complex_quadraple_multiply_i(long double xr, long double xi, long double yr, long double yi, long double zr, long double zi, long double wr, long double wi) {
    long double z1r, z1i, z2r, z2i;
    z1r = complex_multiply_r(xr, xi, yr, yi);
    z1i = complex_multiply_i(xr, xi, yr, yi);
    z2r = complex_multiply_r(zr, zi, wr, wi);
    z2i = complex_multiply_i(zr, zi, wr, wi);
    return (complex_multiply_i(z1r, z1i, z2r, z2i));
}

long double complex_divide_r(long double xr, long double xi, long double yr, long double yi) {
    return ((xr * yr + xi * yi) / (yr * yr + yi * yi));
}

long double complex_divide_i(long double xr, long double xi, long double yr, long double yi) {
    return ((-xr * yi + xi * yr) / (yr * yr + yi * yi));
}

long double complex_root_r(long double xr, long double xi) {
    long double r, theta;
    r = sqrt(xr*xr + xi*xi);
    if (r != 0.0) {
        if (xr >= 0 && xi >= 0) {
            theta = atan(xi / xr);
        }
        else if (xr < 0 && xi >= 0) {
            theta = PI - abs(atan(xi / xr));
        }
        else if (xr < 0 && xi < 0) {
            theta = PI + abs(atan(xi / xr));
        }
        else {
            theta = 2.0 * PI + atan(xi / xr);
        }
        return (sqrt(r) * cos(theta / 2.0));
    }
    else {
        return 0.0;
    }

}    

long double complex_root_i(long double xr, long double xi) {
    long double r, theta;
    r = sqrt(xr*xr + xi*xi);
    if (r != 0.0) {
        if (xr >= 0 && xi >= 0) {
            theta = atan(xi / xr);
        }
        else if (xr < 0 && xi >= 0) {
            theta = PI - abs(atan(xi / xr));
        }
        else if (xr < 0 && xi < 0) {
            theta = PI + abs(atan(xi / xr));
        }
        else {
            theta = 2.0 * PI + atan(xi / xr);
        }
        return (sqrt(r) * sin(theta / 2.0));
    }
    else {
        return 0.0;
    }
}    

long double complex_cuberoot_r(long double xr, long double xi) {
    long double r, theta;
    r = sqrt(xr*xr + xi*xi);
    if (r != 0.0) {
        if (xr >= 0 && xi >= 0) {
            theta = atan(xi / xr);
        }
        else if (xr < 0 && xi >= 0) {
            theta = PI - abs(atan(xi / xr));
        }
        else if (xr < 0 && xi < 0) {
            theta = PI + abs(atan(xi / xr));
        }
        else {
            theta = 2.0 * PI + atan(xi / xr);
        }
        return (pow(r, 1.0 / 3.0) * cos(theta / 3.0));
    }
    else {
        return 0.0;
    }
}    

long double complex_cuberoot_i(long double xr, long double xi) {
    long double r, theta;
    r = sqrt(xr*xr + xi*xi);
    if (r != 0.0) {
        if (xr >= 0 && xi >= 0) {
            theta = atan(xi / xr);
        }
        else if (xr < 0 && xi >= 0) {
            theta = PI - abs(atan(xi / xr));
        }
        else if (xr < 0 && xi < 0) {
            theta = PI + abs(atan(xi / xr));
        }
        else {
            theta = 2.0 * PI + atan(xi / xr);
        }
        return (pow(r, 1.0 / 3.0) * sin(theta / 3.0));
    }
    else {
        return 0.0;
    }
}    

void main() {
    long double a[2], b[2], c[2], d[2], minusd[2];
    long double r, theta;
    cout << "ar?";
    cin >> a[0];
    cout << "ai?";
    cin >> a[1];
    cout << "br?";
    cin >> b[0];
    cout << "bi?";
    cin >> b[1];
    cout << "cr?";
    cin >> c[0];
    cout << "ci?";
    cin >> c[1];
    cout << "dr?";
    cin >> d[0];
    cout << "di?";
    cin >> d[1];

    if (b[0] == 0.0 && b[1] == 0.0 && c[0] == 0.0 && c[1] == 0.0) {
        if (d[0] == 0.0 && d[1] == 0.0) {
            cout << "x1r: 0.0 \n";
            cout << "x1i: 0.0 \n";
            cout << "x2r: 0.0 \n";
            cout << "x2i: 0.0 \n";
            cout << "x3r: 0.0 \n";
            cout << "x3i: 0.0 \n";
        }
        else {
                minusd[0] = -d[0];
                minusd[1] = -d[1];
                r = sqrt(minusd[0]*minusd[0] + minusd[1]*minusd[1]);
                if (minusd[0] >= 0 && minusd[1] >= 0) {
                    theta = atan(minusd[1] / minusd[0]);
                }
                else if (minusd[0] < 0 && minusd[1] >= 0) {
                    theta = PI - abs(atan(minusd[1] / minusd[0]));
                }
                else if (minusd[0] < 0 && minusd[1] < 0) {
                    theta = PI + abs(atan(minusd[1] / minusd[0]));
                }
                else {
                    theta = 2.0 * PI + atan(minusd[1] / minusd[0]);
                }
                cout << "x1r: " << pow(r, 1.0 / 3.0) * cos(theta / 3.0) << "\n";
                cout << "x1i: " << pow(r, 1.0 / 3.0) * sin(theta / 3.0) << "\n";
                cout << "x2r: " << pow(r, 1.0 / 3.0) * cos((theta + 2.0 * PI) / 3.0) << "\n";
                cout << "x2i: " << pow(r, 1.0 / 3.0) * sin((theta + 2.0 * PI) / 3.0) << "\n";
                cout << "x3r: " << pow(r, 1.0 / 3.0) * cos((theta + 4.0 * PI) / 3.0) << "\n";
                cout << "x3i: " << pow(r, 1.0 / 3.0) * sin((theta + 4.0 * PI) / 3.0) << "\n";
            }
        }
        else {
        // find eigenvalues
        long double term0[2], term1[2], term2[2], term3[2], term3buf[2];
        long double first[2], second[2], second2[2], third[2];
        term0[0] = -4.0 * complex_quadraple_multiply_r(a[0], a[1], c[0], c[1], c[0], c[1], c[0], c[1]);
        term0[1] = -4.0 * complex_quadraple_multiply_i(a[0], a[1], c[0], c[1], c[0], c[1], c[0], c[1]);
        term0[0] += complex_quadraple_multiply_r(b[0], b[1], b[0], b[1], c[0], c[1], c[0], c[1]);
        term0[1] += complex_quadraple_multiply_i(b[0], b[1], b[0], b[1], c[0], c[1], c[0], c[1]);
        term0[0] += -4.0 * complex_quadraple_multiply_r(b[0], b[1], b[0], b[1], b[0], b[1], d[0], d[1]);
        term0[1] += -4.0 * complex_quadraple_multiply_i(b[0], b[1], b[0], b[1], b[0], b[1], d[0], d[1]);
        term0[0] += 18.0 * complex_quadraple_multiply_r(a[0], a[1], b[0], b[1], c[0], c[1], d[0], d[1]);
        term0[1] += 18.0 * complex_quadraple_multiply_i(a[0], a[1], b[0], b[1], c[0], c[1], d[0], d[1]);
        term0[0] += -27.0 * complex_quadraple_multiply_r(a[0], a[1], a[0], a[1], d[0], d[1], d[0], d[1]);
        term0[1] += -27.0 * complex_quadraple_multiply_i(a[0], a[1], a[0], a[1], d[0], d[1], d[0], d[1]);
        term1[0] = -27.0 * complex_triple_multiply_r(a[0], a[1], a[0], a[1], d[0], d[1]);
        term1[1] = -27.0 * complex_triple_multiply_i(a[0], a[1], a[0], a[1], d[0], d[1]);
        term1[0] += 9.0 * complex_triple_multiply_r(a[0], a[1], b[0], b[1], c[0], c[1]);
        term1[1] += 9.0 * complex_triple_multiply_i(a[0], a[1], b[0], b[1], c[0], c[1]);
        term1[0] -= 2.0 * complex_triple_multiply_r(b[0], b[1], b[0], b[1], b[0], b[1]);
        term1[1] -= 2.0 * complex_triple_multiply_i(b[0], b[1], b[0], b[1], b[0], b[1]);
        term2[0] = 3.0 * complex_multiply_r(a[0], a[1], c[0], c[1]);
        term2[1] = 3.0 * complex_multiply_i(a[0], a[1], c[0], c[1]);
        term2[0] -= complex_multiply_r(b[0], b[1], b[0], b[1]);
        term2[1] -= complex_multiply_i(b[0], b[1], b[0], b[1]);
        term3[0] = complex_multiply_r(term1[0], term1[1], term1[0], term1[1]);
        term3[1] = complex_multiply_i(term1[0], term1[1], term1[0], term1[1]);
        term3[0] += 4.0 * complex_triple_multiply_r(term2[0], term2[1], term2[0], term2[1], term2[0], term2[1]);
        term3[1] += 4.0 * complex_triple_multiply_i(term2[0], term2[1], term2[0], term2[1], term2[0], term2[1]);
        term3buf[0] = term3[0];
        term3buf[1] = term3[1];
        term3[0] = complex_root_r(term3buf[0], term3buf[1]);
        term3[1] = complex_root_i(term3buf[0], term3buf[1]);

        if (term0[0] == 0.0 && term0[1] == 0.0 && term1[0] == 0.0 && term1[1] == 0.0) {
            cout << "x1r: " << -pow(d[0], 1.0 / 3.0) << "\n";
            cout << "x1i: " << 0.0 << "\n";
            cout << "x2r: " << -pow(d[0], 1.0 / 3.0) << "\n";
            cout << "x2i: " << 0.0 << "\n";
            cout << "x3r: " << -pow(d[0], 1.0 / 3.0) << "\n";
            cout << "x3i: " << 0.0 << "\n";
        }
        else {
            // eigenvalue1
            first[0] = complex_divide_r(complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1]), 3.0 * pow(2.0, 1.0 / 3.0) * a[0], 3.0 * pow(2.0, 1.0 / 3.0) * a[1]);
            first[1] = complex_divide_i(complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1]), 3.0 * pow(2.0, 1.0 / 3.0) * a[0], 3.0 * pow(2.0, 1.0 / 3.0) * a[1]);
            second[0] = complex_divide_r(pow(2.0, 1.0 / 3.0) * term2[0], pow(2.0, 1.0 / 3.0) * term2[1], 3.0 * complex_multiply_r(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 3.0 * complex_multiply_i(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])));
            second[1] = complex_divide_i(pow(2.0, 1.0 / 3.0) * term2[0], pow(2.0, 1.0 / 3.0) * term2[1], 3.0 * complex_multiply_r(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 3.0 * complex_multiply_i(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])));
            third[0] = complex_divide_r(b[0], b[1], 3.0 * a[0], 3.0 * a[1]);
            third[1] = complex_divide_i(b[0], b[1], 3.0 * a[0], 3.0 * a[1]);
            cout << "x1r: " << first[0] - second[0] - third[0] << "\n";
            cout << "x1i: " << first[1] - second[1] - third[1] << "\n";

            // eigenvalue2
            first[0] = complex_divide_r(complex_multiply_r(1.0, -sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), complex_multiply_i(1.0, -sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 6.0 * pow(2.0, 1.0 / 3.0) * a[0], 6.0 * pow(2.0, 1.0 / 3.0) * a[1]);
            first[1] = complex_divide_i(complex_multiply_r(1.0, -sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), complex_multiply_i(1.0, -sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 6.0 * pow(2.0, 1.0 / 3.0) * a[0], 6.0 * pow(2.0, 1.0 / 3.0) * a[1]);
            second[0] = complex_divide_r(complex_multiply_r(1.0, sqrt(3.0), term2[0], term2[1]), complex_multiply_i(1.0, sqrt(3.0), term2[0], term2[1]), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_r(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_i(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])));
            second[1] = complex_divide_i(complex_multiply_r(1.0, sqrt(3.0), term2[0], term2[1]), complex_multiply_i(1.0, sqrt(3.0), term2[0], term2[1]), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_r(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_i(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])));
            third[0] = complex_divide_r(b[0], b[1], 3.0 * a[0], 3.0 * a[1]);
            third[1] = complex_divide_i(b[0], b[1], 3.0 * a[0], 3.0 * a[1]);
            cout << "x2r: " << -first[0] + second[0] - third[0] << "\n";
            cout << "x2i: " << -first[1] + second[1] - third[1] << "\n";

            // eigenvalue3
            first[0] = complex_divide_r(complex_multiply_r(1.0, sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), complex_multiply_i(1.0, sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 6.0 * pow(2.0, 1.0 / 3.0) * a[0], 6.0 * pow(2.0, 1.0 / 3.0) * a[1]);
            first[1] = complex_divide_i(complex_multiply_r(1.0, sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), complex_multiply_i(1.0, sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 6.0 * pow(2.0, 1.0 / 3.0) * a[0], 6.0 * pow(2.0, 1.0 / 3.0) * a[1]);
            second[0] = complex_divide_r(complex_multiply_r(1.0, -sqrt(3.0), term2[0], term2[1]), complex_multiply_i(1.0, -sqrt(3.0), term2[0], term2[1]), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_r(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_i(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])));
            second[1] = complex_divide_i(complex_multiply_r(1.0, -sqrt(3.0), term2[0], term2[1]), complex_multiply_i(1.0, -sqrt(3.0), term2[0], term2[1]), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_r(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_i(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])));
            third[0] = complex_divide_r(b[0], b[1], 3.0 * a[0], 3.0 * a[1]);
            third[1] = complex_divide_i(b[0], b[1], 3.0 * a[0], 3.0 * a[1]);
            cout << "x3r: " << -first[0] + second[0] - third[0] << "\n";
            cout << "x3i: " << -first[1] + second[1] - third[1] << "\n";
        }
    }

    int end;
    cin >> end;
}

Welcome to Stack Overflow! While this code snippet may solve the question, including an explanation [really helps](//meta.stackexchange.com/q/114762) to improve the quality of your post. Remember that you are answering the question for readers in the future, not just the person asking now! Please [edit] your answer to add explanation, and give an indication of what limitations and assumptions apply. — Toby Speight, Nov 11 '16 at 11:13

score -2 · Answer 6 · answered Jul 11 '18 at 18:28

-2

In case someone needs C++ code, you can use this piece of OpenCV:

https://github.com/opencv/opencv/blob/master/modules/calib3d/src/polynom_solver.cpp

answered Jul 11 '18 at 18:28

Stepan Yakovenko

8,670
28
113
206

No! any computer language can answer with correct algorithmic. – Cletrix Apr 10 '19 at 00:03
solve_deg3 for a=3.0, b=3.0, c=1.0, d=1.0 gives an incorrect result = -nan (it should be -1.0) – SergeyT Feb 14 '21 at 01:06

Solving a cubic equation

6 Answers6

Linked