/*
 * The Spar Library - modular math parser
 * Copyright (C) 2000,2001 Davide Angelocola <davide178@inwind.it>
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA  02111-1307, USA.
 *
 */

/* solve_quadratic.c - finds the real roots of a x^2 + b x + c = 0 */
#include <spar/sl_poly.h>
#include <spar/sl_math.h>
#include <spar/sl_math_library.h>


int
sl_poly_solve_quadratic (double a, double b, double c, double *x1, double *x2)
{
  double disc = b * b - 4 * a * c;

  if (disc > 0)
    {
      if (b == 0)
	{
	  double r = SL_ABS (0.5 * sl_sqrn (disc, 2) / a);
	  *x1 = -r;
	  *x2 = r;
	}
      else
	{
	  double sgnb = SL_SGN (b);
	  double temp = -0.5 * (b + sgnb * sl_sqrn (disc, 2));
	  double r1 = temp / a;
	  double r2 = c / temp;

	  if (r1 < r2)
	    {
	      *x1 = r1;
	      *x2 = r2;
	    }
	  else
	    {
	      *x1 = r2;
	      *x2 = r1;
	    }
	}
      return 2;
    }
  else if (disc == 0)
    {
      *x1 = -0.5 * b / a;
      *x2 = -0.5 * b / a;
      return 2;
    }
  else
    {
      return 0;
    }
}