Complex solutions of a cubic equation

[EmilianoG]

I was asked by my professor to code a program which solves the roots of a cubic equation following the discriminant approach found here: https://www.wikihow.com/Solve-a-Cubic-Equation#Using_a_.22Discriminant.22_Approach_sub

After following the steps this is what I have. I'm having trouble getting the solutions to display since I don't know which function/header to use in order to display a complex root, therefore, all of my solutions are showing up as NaN. Any suggestions on what I should change/add in order to make this work?

#include <iostream>
#include <cmath>
using namespace std;
//Solving a cubic equation using the discriminant approach
int main() {
	double a, b, c, d;
cout << "Please enter the values for the a, b, c, d coefficients in order: ";
cin >> a >> b >> c >> d;
double Delta0 = pow(b,2.0) - 3.0*a*c;
double Delta1 = (2.0*pow(b,3.0)- (9.0*a*b*c) + (27.0*pow(a,2.0)*d));
double Disc = (pow(Delta1,2.0) - 4.0*pow(Delta0,3.0))/(27.0*pow(a,2.0));
double C = cbrt(sqrt((pow(Delta1,2.0)- 4.0*pow(Delta1,3.0) + Delta1)/2.0));
double u = (-1.0 + sqrt(-3.0))/2.0;
if(Disc>0){      //Three real solutions
    double x1 = (b+pow(u,1.0)*C+Delta0/(pow(u,1.0)*C))/(3.0*a);
	double x2 = (b+pow(u,2.0)*C+Delta0/(pow(u,2.0)*C))/(3.0*a);
	double x3 = (b+pow(u,3.0)*C+Delta0/(pow(u,3.0)*C))/(3.0*a);
	cout << "D is positive and the roots are: "<< x1 << ", " << x2 << ", " << x3 << endl;	
} else if (Disc==0) {//one or two real solutions
double x1 = (b+pow(u,1.0)*C+Delta0/(pow(u,1.0)*C))/(3.0*a);
double x2 = (b+pow(u,2.0)*C+Delta0/(pow(u,2.0)*C))/(3.0*a);
cout << "D is zero and the roots are: "<< x1 << ", " << x2 << endl;
}else { //D is negative and there is only one solution

 double x1 = (b+pow(u,1.0)*C+Delta0/(pow(u,1.0)*C))/(3.0*a);
cout << "D is negative and the roots is: "<< x1 << endl;
}
	
return 0; 
}

[jonnin]

complex is a type in c++.

complex <double> x;


the lazy / easy fix is to make them all complex.

[icy1]

Note: I should also add that I discarded your WikiHow source as utter garbage after seeing Step6. This "step" randomly introduces a substitution variable casually set to something involving the square root of negative three, without mentioning imaginary, complex, or sourcing anything relevant to back up this craziness.

From what I read on this, the discriminant, as per the actual wiki at https://en.wikipedia.org/wiki/Cubic_function#The_discriminant , should be
double discwiki = pow(b,2)*pow(c,2) - 4*a*pow(c,3) - 4*pow(b,3)*d - 27*pow(a,2)*pow(d,2) + 18*a*b*c*d;

Unfortunately, the article skates around exactly how to get the roots; instead it prefers obscure language like "most convenient form for multiple roots". After some more digging I got to this answer https://stackoverflow.com/a/13328773/1831037 , which is really the pasted JavaScript of a working cubic equation solver at https://www.easycalculation.com/algebra/cubic-equation.php .

I translated it to C++ here, using some of the relatively simple examples from http://www.mash.dept.shef.ac.uk/Resources/web-cubicequations-john.pdf to see if it works. Before you ask, yeah, their value for discriminant is completely different, and if you have a math background you can probably better navigate what their values represent.

Can test/run online at https://repl.it/repls/WoozyValidTabs

#include <iostream>
#include <string>
#include <cmath>
using namespace std;

//Solving a cubic equation using the discriminant approach
int main() 
{
    double a, b, c, d;
    cout << "Please enter space-separated a, b, c, d coefficients:\n";
    cin >> a >> b >> c >> d;
    
    b /= a;
    c /= a;
    d /= a;
    
    double 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;
    term1 = (b/3.0);
    
    double x1_real, x2_real, x3_real;
    double x2_imag, x3_imag;
    string x2_imag_s, x3_imag_s;
    if (disc > 0)   // One root real, two are complex
    {
        s = r + sqrt(disc);
        s = s<0 ? -cbrt(-s) : cbrt(s);
        t = r - sqrt(disc);
        t = t<0 ? -cbrt(-t) : cbrt(t);
        x1_real = -term1 + s + t;
        term1 += (s + t)/2.0;
        x3_real = x2_real = -term1;
        term1 = sqrt(3.0)*(-t + s)/2;
        x2_imag = term1;
        x3_imag = -term1;
        x2_imag_s =  " + "+ to_string(x2_imag) + "i";
        x3_imag_s =  " - "+ to_string(x2_imag) + "i";
    } 
    // The remaining options are all real
    else if (disc == 0)  // All roots real, at least two are equal.
    { 
        x3_imag = x2_imag = 0;
        r13 = r<0 ? -cbrt(-r) : cbrt(r);
        x1_real = -term1 + 2.0*r13;
        x3_real = x2_real = -(r13 + term1);
    }
    // Only option left is that all roots are real and unequal (to get here, q < 0)
    else
    {
        x3_imag = x2_imag = 0;
        q = -q;
        dum1 = q*q*q;
        dum1 = acos(r/sqrt(dum1));
        r13 = 2.0*sqrt(q);
        x1_real = -term1 + r13*cos(dum1/3.0);
        x2_real = -term1 + r13*cos((dum1 + 2.0*M_PI)/3.0);
        x3_real = -term1 + r13*cos((dum1 + 4.0*M_PI)/3.0);
    }
    
    cout << "\nRoots:" << endl <<
            "  x = " << x1_real << endl <<
            "  x = " << x2_real << x2_imag_s << endl <<
            "  x = " << x3_real << x3_imag_s << endl;
    
    return 0; 
}

Example output:

Please enter space-separated a, b, c, d coefficients:
 1 1 1 -3

Roots:
  x = 1
  x = -1 + 1.414214i
  x = -1 - 1.414214i

Please enter space-separated a, b, c, d coefficients:
 1 -6 11 -6

Roots:
  x = 3
  x = 1
  x = 2

[lastchance]

Agree with @icy1 that the OP's original link was seriously flawed (especially steps 4 and 5). For real coefficients you can skirt round the fact that some of the roots can be complex.

#include <iostream>
#include <cmath>
using namespace std;


int main()
{
   const double PI = 4.0 * atan( 1.0 );

   double a, b, c, d;
   cout << "Solve ax^3 + bx^2 + cx + d = 0    ( a/=0, b, c, d all real)\n";
   cout << "\nInput a b c d (separated by spaces): ";
   cin >> a >> b >> c >> d;

   // Reduced equation: X^3 - 3pX - 2q = 0, where X = x-b/(3a)
   double p = ( b * b - 3.0 * a * c ) / ( 9.0 * a * a );
   double q = ( 9.0 * a * b * c - 27.0 * a * a * d - 2.0 * b * b * b ) / ( 54.0 * a * a * a );
   double offset = b / ( 3.0 * a );

   // Discriminant
   double discriminant = p * p * p - q * q;

   cout << "\nRoots:\n";
   if ( discriminant > 0 )           // set X = 2 sqrt(p) cos(theta) and compare 4 cos^3(theta)-3 cos(theta) = cos(3 theta)
   {
      double theta = acos( q / ( p * sqrt( p ) ) );
      double r = 2.0 * sqrt( p );
      for ( int n = 0; n < 3; n++ ) cout << r * cos( ( theta + 2.0 * n * PI ) / 3.0 ) - offset << '\n';
   }
   else 
   {
      double gamma1 = cbrt( q + sqrt( -discriminant ) );
      double gamma2 = cbrt( q - sqrt( -discriminant ) );

      cout << gamma1 + gamma2 - offset << '\n';

      double re = -0.5 * ( gamma1 + gamma2 ) - offset;
      double im = ( gamma1 - gamma2 ) * sqrt( 3.0 ) / 2.0;
      if ( discriminant == 0.0 )                // Equal roots (hmmm, floating point ...)
      {
         cout << re << '\n';
         cout << re << '\n';
      }
      else
      {
         cout << re << " + " << im << " i\n";
         cout << re << " - " << im << " i\n";
      }
   }
}

[EmilianoG]

Thanks guys, very helpful!