Solve Quadratic Equation in C++

P

6

8

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!

Pyrexia answered 22/5, 2009 at 14:12 Comment(3)

The other guys provided good answers so no reason for me to try and outshine them ;) However, if you want a more general solution to the equation ax^2+bx+c=0, remember that a==0 should be a valid value. This would result in a division of zero, so you must take care of this case separately. In this case it would mean that you are left with a linear equation with one root. Cheers ! – Accursed 22/5, 2009 at 14:38

You imply you are concerned about complex roots, but what about complex coefficients in the original equation? – Silence 22/5, 2009 at 14:47

Complex coefficients requires another approach completely. So that will be the next question :) scurrs off and prepares an answer before hand – Accursed 22/5, 2009 at 16:20

A

7

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!

Ammoniate answered 22/5, 2009 at 14:26 Comment(8)

If the sqrt succeeds, the result is >= 0. And if the argument is negative, your program crashes. You should test the sign first, and calculate the sqrt afterwards. If the sign would be negative, you'd set result.first.i = +sqrt(4*ac-bb)/2/a. (Why define your own type if there's a perfectly fine std::pair<std::complex> ?) – Manysided 22/5, 2009 at 14:56

should be delta=bb-4*ac, and take sqrt only when delta>=0. delta=0 or a=0 are valid cases when we have one root. What if a=b=0 and c=1? – Kerbstone 22/5, 2009 at 14:59

What if? That's not a quadratic function, and the /2/a part will fail. This works reasonably well if delta=0, except that tou'll return the same result twice. – Manysided 22/5, 2009 at 15:7

Why are you defining a pair and a complex type? Both already exist! – Notify 22/5, 2009 at 19:7

Portable? Idk.. I wrote it all in the text editor after like 24 hours of no sleep, I'm surprised it even compiles <.< – Ammoniate 22/5, 2009 at 23:53

I would have gone for an intermediate variable to store sqrt(delta): the engineer in me cries whenever he sees unnecessary duplication of a floating point operation... – Cloven 23/5, 2009 at 18:4

Sorry about the late comment but a new user is using this code and they are getting errors. Why are you using sqrt(delta) when delta < 0? Most implementations will return NaN in that case. – Allopathy 5/9, 2017 at 15:53

Must have been a typo, it's supposed to be sqrt(-delta), ie the square root of the absolute value of delta. – Ammoniate 12/9, 2017 at 17:59

N

24

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.

Nadda answered 22/5, 2009 at 21:44 Comment(2)

+1 this is significantly computationally more robust than (-b +/- sqrt(b*b - 4*a*c))/(2a). BTW: since temp could be 0.0, the usually pre-division check is needed. (e. g. a,b,c = 1,0,0). – Waldos 6/6, 2013 at 22:53

sign(b) can be defined as +1 for b >= 0 and -1 for b < 0. – Galop 24/10, 2020 at 15:20

A

7

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;
}