/*
 *  quartic.cpp
 *  
 *
 *  Created by David Yu on 8/14/07.
 *  Copyright 2007 __MyCompanyName__. All rights reserved.
 *
 */

#define Pi 3.1415926
#include <math.h>
#include <complex>
#include <iostream>
using namespace std;

//gives a root of a cubic of form x^3+ax^2+bx+c
void cubic(double a, double b, double c, complex<double> & t0, complex<double> & t1, complex<double> & t2)
{
	//cout << "Passed to function cubic: a = " << a << ", b = " << b << ", c = " << c << "\n";
	complex<double> zeta = complex<double>(-0.5,sqrt(3)/2);
	double discriminant = 4*pow(a,3)*c-pow(a,2)*pow(b,2)+4*pow(b,3)-18*a*b*c+27*pow(c,2);
	if (discriminant < 0)
	{
		//cout << "Cubic discriminant < 0\n";
		double r = sqrt( pow((4.5*a*b-pow(a,3)-13.5*c),2)-6.75*discriminant );
		double theta = atan2(1.5*sqrt(-3*discriminant) , (4.5*a*b-pow(a,3)-13.5*c));
		t0 = (-a + 2*cbrt(r)*cos(theta/3))/3;
		t1 = (-a + 2*cbrt(r)*cos(theta/3-2*Pi/3))/3;
		t2 = (-a + 2*cbrt(r)*cos(theta/3+2*Pi/3))/3;
	} else if (discriminant > 0)
	{
		//cout << "Cubic discriminant > 0\n";
		double s0 = -a;
		double s1 = cbrt( 4.5*a*b-pow(a,3)-13.5*c+1.5*sqrt(3.0*discriminant) );
		double s2 = cbrt( 4.5*a*b-pow(a,3)-13.5*c-1.5*sqrt(3.0*discriminant) );
		//cout << "Lagrange resolvents are " << s0 << ", " << s1 << ", " << s2 << "\n";
		
		t0 = (s0 + s1 + s2)/3.0;
		t1 = (s0 + zeta*zeta*s1 + zeta*s2)/3.0;
		t2 = (s0 + zeta*s1 + zeta*zeta*s2)/3.0;
	}
}
	

//solves a quartic of form x^4+ax^3+bx^2+cx+d
void quartic(double a, double b, double c, double d, complex<double> & r1, complex<double> & r2, complex<double> & r3, complex<double> & r4)
{
	//cout << "quartic is x^4 + " << a << "x^3 + " << b << "x^2 + " << c << "x + " << d << "\n";
	double discriminant = (256*pow(d,3)-pow(d,2)*(27*pow(a,4)-144*pow(a,2)*b+192*a*c)-pow(c,2)*(27*pow(c,2)+4*pow(a,3)*c+4*pow(b,3)-pow(a*b,2))-2*d*(a*b*c*(9*pow(a,2)-40*b)-2*pow(b,3)*(pow(a,2)-4*b)-3*pow(c,2)*(pow(a,2)-24*b)));
	//cout << "discriminant = " << discriminant << "\n";
	complex<double> t0, t1, t2, r1plusr2, r3plusr4, r1r2, r3r4;
	cubic(-2*b, pow(b,2)+a*c-4*d, pow(c,2)+pow(a,2)*d-a*b*c, t0, t1, t2);
	//cout << "Solutions of resolvent cubic are " << t0 << ", " << t1 << ", " << t2 << "\n";
	r1plusr2 = (-a+sqrt(pow(a,2)-4.0*t0))/2.0;
	r3plusr4 = (-a-sqrt(pow(a,2)-4.0*t0))/2.0;
	r1r2 = (t1+t2-t0+sqrt(pow(t1+t2-t0,2)-16.0*d))/4.0;
	r3r4 = (t1+t2-t0-sqrt(pow(t1+t2-t0,2)-16.0*d))/4.0;
	r1 = (r1plusr2 + sqrt(pow(r1plusr2,2)-4.0*r1r2))/2.0;
	r2 = (r1plusr2 - sqrt(pow(r1plusr2,2)-4.0*r1r2))/2.0;
	r3 = (r3plusr4 + sqrt(pow(r3plusr4,2)-4.0*r3r4))/2.0;
	r4 = (r3plusr4 - sqrt(pow(r3plusr4,2)-4.0*r3r4))/2.0;
	//cout << "Check: " << pow(r1,4)+a*pow(r1,3)+b*pow(r1,2)+c*r1+d << ", " << pow(r2,4)+a*pow(r2,3)+b*pow(r2,2)+c*r2+d << ", " << pow(r3,4)+a*pow(r3,3)+b*pow(r3,2)+c*r3+d << ", " << pow(r4,4)+a*pow(r4,3)+b*pow(r4,2)+c*r4+d << "\n";
	//cout << "Quartic roots are " << r1 << ", " << r2 << ", " << r3 << ", " << r4 << "\n";
}