/*
 * quad_eq.c
 * This program reads the coefficients of a quadratic equation and solves
 * the equation.
 */

#include <stdio.h>
#include <math.h>

/* Prototypes: */
int solve_quad (double a, double b, double c,
                double *x1, double *x2);

int main ()
{
  double    a, b, c;          /* The equation coefficients. */
  double    x1, x2;           /* Its solutions. */
  int       n_sols;           /* The number of real solutions. */
  
  /* Read the coeffcients. */
  printf ("Please enter the three coefficients of the quadratic equation: ");
  scanf ("%lf %lf %lf", &a, &b, &c);

  /* Solve the equation. */
  n_sols = solve_quad (a, b, c,
                       &x1, &x2);

  /* Print the solutions. */
  switch (n_sols)
  {
  case 0:
    printf ("The equation has no real solutions.\n");
    break;

  case 1:
    printf ("The equation has a single solution: (x = %g).\n", x1);
    break;

  case 2:
    printf ("The equation has two solutions: (x1 = %g) (x2 = %g).\n", x1, x2);
    break;

  default:
    printf ("You must be kidding!\n");
    return (1);
  }

  return (0);
}

/* ------------------------------------------------------------------------
 * Function: solve_quad
 * Purpose : Solve the quadratic equation: (a*x^2 + bx + c = 0).
 * Input   : a, b, c - The equation coefficients.
 * Output  : x1, x2  - The solutions.
 * Returns : The number of real solutions to the equation (0, 1 or 2).
 */
int solve_quad (double a, double b, double c,
                double *x1, double *x2)
{
  double   disc;
  double   sqrt_disc;

  /* Check if this is really a linear equation. */
  if (a == 0)
  {
    /* The equation is: (b*x + c = 0). */
    if (b == 0)
      return (0);

    *x1 = -c / b;
    return (1);
  }

  /* Compute the discriminant and act according to its sign. */
  disc = b*b - 4*a*c;

  if (disc < 0)
  {
    return (0);
  }
  else if (disc == 0)
  {
    *x1 = -b / (2*a);
    return (1);
  }

  /*                          -b +/- sqrt(b^2 - 4ac)
   * Use the formula: x1,2 = ------------------------
   *                                   2a
   */
  sqrt_disc = sqrt (disc);
  
  *x1 = (-b + sqrt_disc) / (2*a);
  *x2 = (-b - sqrt_disc) / (2*a);
  return (2);
}