Функция подсчёта корней кубического уравнения

int Cubic(double a, double b, double c,double d) {
    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 h = ((p*p*p) / 27.0);
    double j = (q*q) / 4.0;
    double s = (4 * (3 * a*c - b*b)*(3 * a*c - b*b)*(3 * a*c - b*b) + (2 * b*b*b - 9 * a*b*c + 27 * a*a*d)*(2 * b*b*b - 9 * a*b*c + 27 * a*a*d))/ (2916 * a*a*a*a*a*a);
    if (s < 0)
        {
            double f;
            if (q < 0)
            {
                f = atan(sqrt(-s) / (-q / 2));
            }
            else if (q > 0)
            {
                f = atan(sqrt(-s) / (-q / 2)) + M_PI;
            }
            else
            {
                f = (M_PI / 2);
            }
            double k;
            if (p == 0)
                k = 0;
            else
                k = sqrt(-p / 3.0);
            double x1 = (2.0*k*cos(f / 3.0) - b / 3.0*a);
            double x2 = (2.0*k*cos(f / 3.0 + (2.0*M_PI) / 3.0) - b / 3.0*a);
            double x3 = (2.0*k*cos(f / 3.0 + (4.0*M_PI) / 3.0) - b / 3.0*a);
            if (x1 != x2&&x1 != x3&&x2 != x3)
                return 1;
            else
                return 2;
        }
        else
            if (s > 0)
            {
                double mi;
                double pl;
                double mic;
                double plc;
                mi = -q / 2.0 - sqrt(s);
                pl = -q / 2.0 + sqrt(s);
                if (pl == 0)
                    plc = 0;
                else
                    plc = pow(abs(pl), 1 / 3.0)*(abs(pl) / pl);
                if (mi == 0)
                    mic = 0;
                else
                    mic = pow(abs(mi), 1 / 3.0)*(abs(mi) / mi);
                double x1 = plc + mic - (b / (3.0*a));
                double x2 = (-1 / 2.0)*(plc + mic) - (b / (3.0*a));
                double x2i = (sqrt(3.0) / 2.0)*(plc - mic);
                double x3 = (-1 / 2.0)*(plc + mic) - (b / (3.0*a));
                double x3i = -(sqrt(3.0) / 2.0)*(plc - mic);
                if (x3i != x2i&&x2 != x3)
                    return 2;
                else
                    return 1;
            }
            else
            {
                double k = -q / 2.0;
                double f;
                if (k == 0)
                    f = 0;
                else
                    f = pow(abs(k), 1 / 3.0)*(abs(k) / k);
                double x1 = (2.0*f - (b / (3.0*a)));
                double x2 = (-1.0*f - (b / (3.0*a)));
                double x3 = (-1.0*f - (b / (3.0*a)));
                if (x1 != x2&&x1 != x3&&x2 != x3)
                    return 1;
                else
                    return 2;
            }
}