Code:
#include <stdio.h>
#include <math.h>
#define PI (3.141592653589793238462643)
#define a 0.448
#define b 4.29e-5
#define R 8.318
#define n 1.0
#define P 4.0e7
#define T 400.0
float Q, S, a1, a2, a3, A11, x1, x2, x3, theta, test1, test2, test3, Qcubed, Ssquared;
//Sign function prototype
float sign(float S);
int main()
{
//Code that handles coefficents
a1 = -((P*n*b + n*R*T)/P);
a2 = (a*n*n/P);
a3 = -(a*b*n*n*n/P);
//Doing some calculations to simplify our life
Q = (a1*a1 - 3.*a2)/(9.0);
S = (2.*a1*a1*a1 - 9.*a1*a2 + 27.*a3)/(54.0);
Qcubed = Q*Q*Q;
Ssquared = S*S;
//Printing out all the variables used so far to make sure no inf,NaN or other errors have occured
printf("%f = Q, %f = S, %f = a1, %f = a2 %f = a3\n", Q, S, a1, a2, a3);
printf("%f = sign(S) %f = A11\n", sign(S), A11);
printf("%f = Qcubed - Ssquared\n", Qcubed - Ssquared);
//The first case, if our Q^3 - S^2 > 0
//The cubic equation has 3 real, easy to find roots, found as below
if (Qcubed - Ssquared > 0)
{
theta = acos(S/sqrt(Qcubed));
x1 = -2*sqrt(Q)*cos(theta/3.) - b/3.;
x2 = -2*sqrt(Q)*cos((theta + 2*PI)/3.) - a1/3.;
x3 = -2*sqrt(Q)*cos((theta + 4*PI)/3.) - a1/3.;
printf("The cubic has three real roots, %f = x1, %f = x2, %f = x3\n", x1, x2, x3);
}
//The second case, if Q^3 - S^2 <= 0, first we calculate A11
//After A11 has been calculated we can decide which equation to use to solve the cubic
//The cubic has only one real root as expressed as below
if (Qcubed - Ssquared <= 0)
{
//Calculating A11 to find out which equation to use to solve the cubic
A11 = -sign(S)*pow((sqrt(Ssquared - Qcubed) + fabs(S)),1/3.);
//If A11 is zero, then we use this to find the root
if (A11 = 0)
{
x1 = A11 - a1/3.;
printf("The cubic has only one real root, %f = x1\n", x1);
}
//If A11 is anything APART from 0 then we use this to find the root
else
{
x1 = A11 + Q/A11 - a1/3.;
printf("The cubic has only one real root, %f = x1\n", x1);
}
}
//Accuracy testing, subbing x1, x2, x3 values into equation
//If these equal 0, then we know our solutions are genuine
test1 = x1*x1*x1 + a1*x1*x1 + a2*x1 + a3;
test2 = x2*x2*x2 + a1*x2*x2 + a2*x2 + a3;
test3 = x3*x3*x3 + a1*x3*x3 + a2*x3 + a3;
printf("Subbing x1, x2, and x3 into the equation yields %f = test1, %f = test2, %f = test3\n", test1, test2, test3);
printf("If these are equal to 0, our solutions are genuine\n");
return 0;
}
//The sign function
float sign(float S)
{
if (S > 0)
return 1;
else if (S < 0)
return -1;
else
return 0;
}