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));
}