CubeEquation

#include <stdio.h>
#include <stdlib.h>
#define _USE_MATH_DEFINE
#include <math.h>
#include <float.h>
typedef enum { REAL, IMAGE } complex_t;
int main(void) {
  double a, b, c, d, p, q, D, x1, x2, x3, y1, y2, y3, Im, u, v, d_u, d_v, phi;
  complex_t kind;
  printf("3次方程式(ax^3+bx^2+cx+d=0)の各係数:\n");
  scanf("%lf%lf%lf%lf", &a, &b, &c, &d);
  if (a == 0.) exit(1);
  printf("3次方程式 (%g x^3 %+g x^2 %+g x %+g = 0) の解は\n", a, b, c, d);
  p = (3 * a * c - b * b) / (9 * a * a);
  q = (2 * b * b * b - 9 * a * b * c + 27 * a * a * d) / (54 * a * a * a);
  D = q * q + p * p * p;
  if (fabs(D) <= DBL_MIN) {
    if (fabs(q) <= DBL_MIN) y1 = y2 = y3 = 0;
    else { y1 = 2 * pow(-q, 1.0 / 3.0); y2 = y3 = -y1 / 2; }
    kind = REAL;
  }
  else if (D < 0) {
    phi = acos(-q / pow(-p * p * p, 1.0 / 3.0));
    y1 = 2 * sqrt(-p) * cos(phi / 3);
    y2 = -2 * sqrt(-p) * cos(phi / 3 + M_PI / 3);
    y3 = -2 * sqrt(-p) * cos(phi / 3 - M_PI / 3);
    kind = REAL;
  }
  else {
    d_u = -q + sqrt(D); d_v = -q - sqrt(D);
    if (d_u >= 0) u = pow(d_u, 1.0 / 3.0); else u = -pow(-d_u, 1.0 / 3.0);
    if (d_v >= 0) v = pow(d_v, 1.0 / 3.0); else v = -pow(-d_v, 1.0 / 3.0);
    y1 = u + v; y2 = y3 = -y1 / 2; Im = sqrt(3) / 2 * (u - v); kind = IMAGE;
  }
  x1 = y1 - b / (3 * a); x2 = y2 - b / (3 * a); x3 = y3 - b / (3 * a);
  if (kind == REAL) printf("x1 = %g, x2 = %g, x3 = %g\n", x1, x2, x3);
  else printf("x1 = %g, x2 = %g + %g i, x3 = %g - %g i\n", x1, x2, Im, x3, Im);
  return 0;
}