/*
 * 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.
 *
 */

#include <spar/sl_poly.h>
#include <spar/sl_math.h>
#include <spar/sl_math_library.h>
#include <spar/sl_complex.h>

/* 
 * finds the real roots of a 
 *
 * 		a z^2 + b z + c = 0 
 */
int
sl_poly_zsolve_quadratic (double a, double b, double c, sl_complex * z0,
			  sl_complex * z1)
{
  double disc = sl_pow_2 (b) - 4 * a * c;

  if (SL_IS_POSITIVE (disc))
    {
      if (b == 0)
	{
	  double s = sl_abs (0.5 * sl_sqrt (disc) / a);
	  SL_COMPLEX_RE (*z0) = -s;
	  SL_COMPLEX_IM (*z0) = 0;
	  SL_COMPLEX_RE (*z1) = s;
	  SL_COMPLEX_IM (*z1) = 0;
	}
      else
	{
	  double sgnb = (b > 0 ? 1 : -1);
	  double temp = -0.5 * (b + sgnb * sl_sqrt (disc));
	  double r1 = temp / a;
	  double r2 = c / temp;

	  if (r1 < r2)
	    {
	      SL_COMPLEX_RE (*z0) = r1;
	      SL_COMPLEX_IM (*z0) = 0;
	      SL_COMPLEX_RE (*z1) = r2;
	      SL_COMPLEX_IM (*z1) = 0;
	    }
	  else
	    {
	      SL_COMPLEX_RE (*z0) = r2;
	      SL_COMPLEX_IM (*z0) = 0;
	      SL_COMPLEX_RE (*z1) = r1;
	      SL_COMPLEX_IM (*z1) = 0;
	    }
	}
      return 2;
    }
  else if (SL_IS_ZERO (disc))
    {
      SL_COMPLEX_RE (*z0) = -0.5 * b / a;
      SL_COMPLEX_IM (*z0) = 0;

      SL_COMPLEX_RE (*z1) = -0.5 * b / a;
      SL_COMPLEX_IM (*z1) = 0;

      return 2;
    }
  else
    {
      double s = sl_abs (0.5 * sl_sqrt (-disc) / a);

      SL_COMPLEX_RE (*z0) = -0.5 * b / a;
      SL_COMPLEX_IM (*z0) = -s;

      SL_COMPLEX_RE (*z1) = -0.5 * b / a;
      SL_COMPLEX_IM (*z1) = s;

      return 2;
    }
}