実装

三次方程式を関数として以下のように実装する。

/* 三次方程式 */
int WH_Algebra__solveCubicEquation
(double tolerance, 
 double a, double b, double c,
 double *OUT__x_0, double *OUT__x_1, double *OUT__x_2)
/*
  入力引数:
    tolerance は判別式に対するトレランス値
    a, b, c はスカラーで、方程式の係数
  出力引数:
    OUT__x_0, OUT__x_1, OUT__x_2 はスカラーで、解
      もし、戻り値が3なら、3つとも有効
      もし、戻り値が1なら、OUT__x_0のみ有効
  戻り値:
    実根の数を返す。
*/
{
  double Q = (a * a - 3 * b) / 9;
  double R = ( 2 * a * a * a - 9 * a * b + 27 * c) / 54;
  double d = Q * Q * Q - R * R;

  if (fabs (d) < tolerance) {
    d = 0.0;
  }
  if (d < 0.0) {
    double e = pow (sqrt (-d) + fabs (R), 1.0 / 3);

    if (0.0 < R) {
      *OUT__x_0 = -(e + Q / e) - a / 3;
    } else {
      *OUT__x_0 = (e + Q / e) - a / 3;
    }
    return 1;
  } else {
    double theta = acos (R / sqrt (Q * Q * Q));

    *OUT__x_0 
      = -2 * sqrt (Q) * cos (theta / 3) - a / 3;
    *OUT__x_1 
      = -2 * sqrt (Q) * cos ((theta + 2 * M_PI) / 3) - a / 3;
    *OUT__x_2 
      = -2 * sqrt (Q) * cos ((theta + 4 * M_PI) / 3) - a / 3;
    return 3;
  }
}