/*
	Author      : Mark Steeves, mlstvs
	Date        : 17 October 1998
	Task        : Solution to Assignment 2, Problem 2
	Description : This program will solve cubic equations of the form:
			1x^3 + Ax^2 + Bx + c = 0
	Filename    : cubic.cpp
*/

#include <math.h>
#include <iostream.h>

int main()
{
	double a, b, c, q, r, t, u, v, sgn, x1, x2, x3;
	const double PI = 3.141592654;

	//Get the value of a, b, c.  Also describe the program.
	cout<<"This program will solve a cubic equation of the form: \n"
		 <<"      1X^3 + aX^2 + bX + c = 0 \n\n"
		 <<"Please enter the value of a:  ";
	cin >> a;

	cout<<"Please enter the value of b:  ";
	cin >> b;

	cout<<"Please enter the value of c:  ";
	cin >> c;

	//Calculation of q and r
	q = (pow(a, 2.0) - 3.0 * b) / 9.0;
	r = (2.0*(pow(a, 3.0)) - (9.0*a*b) + (27.0 * c))/54.0;

	//Is r^2 less than or equal to q^3?
	if( r*r <= q*q*q )
	{//There are three solutions...
		//Find the value of t
		t = acos(r/sqrt(q*q*q));

		//Find the roots:
		x1 = -2.0 * sqrt(q) * cos(t/3.0) - a/3.0;
		x2 = -2.0 * sqrt(q) * cos((t + 2.0 * PI)/3.0) - a/3.0;
		x3 = -2.0 * sqrt(q) * cos((t - 2.0 * PI)/3.0) - a/3.0;

		//Output the roots:
		cout.setf(ios::fixed);
		cout.setf(ios::showpoint);
		cout.precision(5);
		cout<<"\nThere are three solutions: \n"
			 <<"\tX1 = "<<x1<<endl
			 <<"\tX2 = "<<x2<<endl
			 <<"\tX3 = "<<x3<<endl;
	}
	else
	{//There is one solution...
		//Find what sgn(r) is:
		if(r >= 0)
			sgn = 1.0;
		else
			sgn = -1.0;

		//calculate u
		u = -1*sgn*(pow(fabs(r) + sqrt(r*r - q*q*q), 1.0/3.0));
				/*NOTES: It is very important to make sure that you are not doing
				integer division.  The second parameter of the pow function takes a
				double value.  If you pass 1/3 as the second parameter, this is
				evaluated to be 0.  This is because 1/3 is actually 0.33333, and if
				this is an integer, it is truncated, so the number that is actually
				passed is 0.  If ever in doubt, put the decimal point in anyway.*/

		//Find the value of v:
		if (u != 0)
			v = q/u;
		else
			v = 0;

		//and finally, find the value of x1.
		x1 = u + v - a/3.0;

		//Output the results
		cout.setf(ios::fixed);
		cout.setf(ios::showpoint);
		cout.precision(5);

		cout<<"\nThere is one solution: \n"
			 <<"\tX1 = " << x1<<endl;
	}
	return(0);
}