#include <complex>
#include <cstdio>
#include <cstdlib>
#include <iostream>

#ifndef printf_s
#define printf_s printf
#endif

#ifndef scanf_s
#define scanf_s scanf
#endif

#if defined(_WIN32) || defined(WIN32) || defined(_WIN64) || defined(WIN64)
#define clear() system("cls")
#elif defined(__unix__) || defined(linux)
#define clear() system("clear")
#endif

#define eps (double)0.00000001
#define PI (double)3.141592

using namespace std;

int main(void) {
  double fa, fb, fc, fd, fnp, fq, fD;

  clear();

  printf_s("Podaj liczbe a: ");
  cin >> fa;
  if (abs(fa) < eps) {
    printf_s("a == 0\n");
    exit(-1);
  }

  printf_s("Podaj liczbe b: ");
  cin >> fb;

  printf_s("Podaj liczbe c: ");
  cin >> fc;

  printf_s("Podaj liczbe d: ");
  cin >> fd;

  /*
   *hornera
   */
  // complex<double> w1 = ((a * z1 + b) * z1 + c) * z1 + d;

  /*
   * klasyczna
   */
  // complex<double> w1p = a * pow(z1, 3) + b * pow(z1, 2) + c * z1 + d;

  fnp = (-(fb * fb * fb) / (27 * fa * fa * fa)) + ((fb * fc) / (6 * fa * fa)) -
        (fd / (2 * fa));
  fq = (fc / (3 * fa)) - ((fb * fb) / (9 * fa * fa));

  fD = (fnp * fnp) + (fq * fq * fq);

  if (fD > 0) {
    /*
     * jedyn rzeczywisty, dwa zespolōne
     * src: https://pl.wikipedia.org/wiki/R%C3%B3wnanie_sze%C5%9Bcienne
     */
    double fu, fv, fx, fre_z, fim_z;
    double fe = sqrt(3.0) / 2.0;

    fu = cbrt(fnp + sqrt(fD));
    fv = cbrt(fnp - sqrt(fD));
    fx = fu + fv - (fb / (3 * fa));
    fre_z = -((fu / 2.0) + (fv / 2.0) + (fb / (3.0 * fa)));

    fim_z = (fu * fe) - (fv * fe);

    complex<double> z1(fx, 0);
    complex<double> z2(fre_z, fim_z);
    complex<double> z3(fre_z, -fim_z);

    cout << z1 << endl << z2 << endl << z3 << endl;
  } else if (fD < 0) {
    /*
     * trzi rzeczywiste
     * src: http://www.1728.org/cubic2.htm
     */
    double ff = (((3 * fc) / fa) - ((fb * fb) / (fa * fa))) / 3.0;
    double fg = (((2 * fb * fb * fb) / (fa * fa * fa)) -
                 ((9 * fb * fc) / (fa * fa)) + ((27 * fd) / fa)) /
                27.0;
    double fh = ((fg * fg) / 4.0) + ((ff * ff * ff) / 27.0);
    double fi = sqrt(((fg * fg) / 4) - fh);
    double fj = cbrt(fi);
    double fk = acos(-(fg / (2 * fi)));
    double fL = -fj;
    double fM = cos(fk / 3);
    double fN = sqrt(3.0) * sin(fk / 3);
    double fP = -(fb / (3 * fa));
    double fdelta_lt_0_x1 = (2 * fj * fM) + fP;
    double fdelta_lt_0_x2 = fL * (fM + fN) + fP;
    double fdelta_lt_0_x3 = fL * (fM - fN) + fP;
    cout << "x_1 = " << fdelta_lt_0_x1 << endl;
    cout << "x_2 = " << fdelta_lt_0_x2 << endl;
    cout << "x_3 = " << fdelta_lt_0_x3 << endl;
  } else if (abs(fD) < eps) {
    /*
     * max. dwa rzeczywiste
     * src: https://pl.wikipedia.org/wiki/R%C3%B3wnanie_sze%C5%9Bcienne
     */
    double fdelta0_x1 = cbrt(fq / 2.0) - (fb / (3.0 * fa));
    double fdelta0_x2 = ((-2) * cbrt(fq / 2.0)) - (fb / (3.0 * fa));
    if (abs(abs(fdelta0_x1) - abs(fdelta0_x2)) < eps) {
      cout << "x_0 = " << fdelta0_x1 << endl;
    } else {
      cout << "x_1 = " << fdelta0_x1 << endl;
      cout << "x_2 = " << fdelta0_x2 << endl;
    }
  } else {
    cout << "Cos je niy on" << endl;
  }

  getchar();
  getchar();

  exit(0);
}