## C Language Functions for the Cubic Formula

```
void CubicFormula(FILE *f, double a, double b, double c, double d)
{
double x1,x2,x3;
double p,q,r,pAbs;
double discriminant,sqrtD,phi,u,v,im;

p=(3.0*a*c - b*b) / (9.0*a*a);
q=(9.0*a*b*c - 27.0*a*a*d - 2.0*b*b*b) / (54.0*a*a*a);
r=-b/(3.0*a);
discriminant=p*p*p+q*q;

if(fabs(discriminant) < 1.0e-14) {
fprintf(f,"  two real roots (one of the roots is a double root)\n");
u=CubeRoot(q);
v=CubeRoot(q);
x1 = u+v+r;
x2 = -(u+v)/2.0+r;
fprintf(f,"  x1 = %g\n",x1);
fprintf(f,"  x2 = %g  (the double root)\n",x2);
} else if(discriminant>0) {
sqrtD=sqrt(discriminant);
u=CubeRoot(q+sqrtD);
v=CubeRoot(q-sqrtD);
x1 = u+v+r;
fprintf(f,"  the one real root = %g\n",x1);
im = (u-v)*sqrt(3.0)/2.0;
x2 = -(u+v)/2.0+r;
fprintf(f,"  the imaginary roots = %g + and - i * %g\n",x2,im);
} else {
fprintf(f,"  three real roots\n");
pAbs=fabs(p);
phi=acos(q/sqrt(pAbs*pAbs*pAbs));
x1 = 2.0*sqrt(pAbs)*cos(phi/3.0)+r;
x2 = -2.0*sqrt(pAbs)*cos((phi+M_PI)/3.0)+r;
x3 = -2.0*sqrt(pAbs)*cos((phi-M_PI)/3.0)+r;
fprintf(f,"  x1 = %g\n",x1);
fprintf(f,"  x2 = %g\n",x2);
fprintf(f,"  x3 = %g\n",x3);
}
fprintf(f,"\n");
}

double CubeRoot(double x)
{
if(x==0.0) {
return 0.0;
}
if(x>0.0) {
return pow(x,(1.0/3.0));
}
return -pow(-x,(1.0/3.0));
}
```

(email: ``` Jeff Stetekluh```)