Решение кубического уравнения методом Кардано

#include <iostream>
#include <cmath>
#include <algorithm>
#include <vector>
#include <limits>
#include <queue>
#include <string>
#include <set>
#include <stack>
#include <tuple>
 
using namespace std;
 
//Cube root
double croot(double x){
    if (x < 0)
        return -pow(-x, 1.0/3.0);
    return pow(x, 1.0/3.0);
}
 
int main(){         
    freopen("input.txt", "rt", stdin);
    freopen("output.txt", "wt", stdout);
    //Ax + By + Cz + D = 0; (A != 0)
    const double M_PI = 3.141592653589793238462;
    //-----
    double A, B, C, D;
    cin >> A >> B >> C >> D;
    if (A == 0){
        cout << "X is any";
        return 0;
    }
    //-----
    double p = (3.0*A*C-B*B)/(3.0*A*A);
    double q = (2.0*B*B*B-9.0*A*B*C+27.0*A*A*D)/(27.0*A*A*A);
    double S = (q*q/4.0) + (p*p*p/27.0);
    //-----
    double F;
    if (q == 0)
        F = M_PI/2.0;
    if (q < 0)
        F = atan(-2.0*sqrt(-S)/q);
    if (q > 0)
        F = atan(-2.0*sqrt(-S)/q) + M_PI;
    //-----
    pair<double, double> x[3];
    for (int i = 0; i < 3; i++)
        x[i].first = x[i].second = 0;
    if (S < 0){
        x[0].first = 2.0*sqrt(-p/3.0)*cos(F/3.0)-B/(3.0*A);
        x[1].first = 2.0*sqrt(-p/3.0)*cos((F/3.0)+2.0*M_PI/3.0)-B/(3.0*A);
        x[2].first = 2.0*sqrt(-p/3.0)*cos((F/3.0)+4.0*M_PI/3.0)-B/(3.0*A);
    }
    if (S == 0){
        x[0].first = 2.0*croot(-q/2.0)-B/(3.0*A);
        x[1].first = -croot(-q/2.0)-B/(3.0*A);
        x[2].first = -croot(-q/2.0)-B/(3.0*A);
    }
    if (S > 0){
        double temp1 = croot((-q/2.0)+sqrt(S)) + croot((-q/2.0)-sqrt(S));
        double temp2 = croot((-q/2.0)+sqrt(S)) - croot((-q/2.0)-sqrt(S));
        x[0].first = temp1 - B/(3.0*A);
        x[1].first = -temp1/2.0 - B/(3.0*A); x[1].second = sqrt(3)*temp2/2.0;
        x[2].first = -temp1/2.0 - B/(3.0*A); x[2].second = -sqrt(3)*temp2/2.0;
    }
    //-----
    printf("x1 = %.15lf + i * %.15lf\n", x[0].first, x[0].second);
    printf("x2 = %.15lf + i * %.15lf\n", x[1].first, x[1].second);
    printf("x3 = %.15lf + i * %.15lf\n", x[2].first, x[2].second);
    return 0;
}