Solve Quadratic Equation in C++

Question

Hello all,

I am trying to write a function in C++ that solves for X using the quadratic equation. This is what I have written initially, which seems to work as long as there are no complex numbers for an answer:

float solution1 = (float)(-1.0 * b) + (sqrt((b * b) - (4 * a * c)));
solution1 = solution1 / (2*a);

cout << "Solution 1: " << solution1 << endl;

float solution2 = (float)(-b) - (sqrt((b*b) - (4 * a * c)));
solution2 = solution2 / (2*a);
cout << "Solution 2: " << solution2;

If, for example, I use the equation: x^2 - x - 6, I get the solution 3, -2 correctly.

My question is how would I account for complex numbers....for example, given the equation:

x^2 + 2x + 5

Solving by hand, I would get -1 + 2i, -1 - 2i.

Well, I guess two question, can I write the above better and also make it account for the complex number?

Thanks for any help!

Answer 1

You more or less have it, just check to see if the part that's inside the square root is negative and then keep track of that separately in your reductions.

Answer 2

Something like this would work:

struct complex { double r,i; }
struct pair<T> { T p1, p2; }

pair<complex> GetResults(double a, double b, double c)
{
  pair<complex> result={0};

  if(a<0.000001)    // ==0
  {
    if(b>0.000001)  // !=0
      result.p1.r=result.p2.r=-c/b;
    else
      if(c>0.00001) throw exception("no solutions");
    return result;
  }

  double delta=b*b-4*a*c;
  if(delta>=0)
  {
    result.p1.r=(-b-sqrt(delta))/2/a;
    result.p2.r=(-b+sqrt(delta))/2/a;
  }
  else
  {
    result.p1.r=result.p2.r=-b/2/a;
    result.p1.i=sqrt(delta)/2/a;
    result.p2.i=-sqrt(delta)/2/a;
  }

  return result;
}

That way you get the results in a similar way for both real and complex results (the real results just have the imaginary part set to 0). Would look even prettier with boost!

edit: fixed for the delta thing and added a check for degenerate cases like a=0. Sleepless night ftl!

Answer 3

As a sidenote: When dividing, always check if the denominator is not zero. And remember for floating point numbers to use something like:

#inlcude <float.h>
if (fabs(a) < FLT_EPSILON)
    then a is considered 0

Answer 4

You could basically just use std::complex<float> instead of float to get support for complex numbers.

Answer 5

Nicking the idea from Blindy:

typedef std::complex<double> complex;
using std::pair;
pair<complex> GetResults(double a, double b, double c)
{
  double delta=(b*b-4*a*c);
  double inv_2a = 1/2/a;
  if(delta >= 0) {
    double root = sqrt(delta);
    return std::make_pair(
        complex((-b-root)*inv_2a),
        complex((-b+root)*inv_2a);
  } else {
    double root = sqrt(-delta);
    return std::make_pair(
        complex(-b*inv_2a, -root*inv_2a)),
        complex(-b*inv_2a, +root*inv_2a)));
  }
}

Answer 6

An important note to all of this. The solutions shown in these responses and in the original question are not robust.

The well known solution (-b +- sqrt(b^2 - 4ac)) / 2a is known to be non-robust in computation when ac is very small compered to b^2, because one is subtracting two very similar values. It is better to use the lesser known solution 2c / (-b -+ sqrt(b^2 -4ac)) for the other root.

A robust solution can be calculated as:

temp = -0.5 * (b + sign(b) * sqrt(b*b - 4*a*c);
x1 = temp / a;
x2 = c / temp;

The use of sign(b) ensures that we are not subtracting two similar values.

For the OP, modify this for complex numbers as shown by other posters.