Cubic solver using Nickalls method

cubic_solver.c

void cubic_solver(double A, double B, double C, double D, double *roots) {
    double xN, yN;                // yN = f(xN)                        
    double del;
    double h;
    double theta;
    double alpha, beta, gamma;    // roots

    xN = -B/(3*A);
    yN = (A*xN*xN*xN) + (B*xN*xN) + (C*xN) + D;
    del = sqrt( (B*B - 3*A*C) / ((9*A*A)) );
    h = 2*A*del*del*del;

    // CASE 1: 1 real root
    if ((yN*yN) > (h*h)) {
        alpha = xN + cbrt((0.5*A) * (-yN + sqrt(yN*yN - h*h))) + 
                cbrt((0.5*A) * (-yN - sqrt(yN*yN - h*h)));

        roots[0] = alpha;
        
        return *roots;
    }
    // CASE 2: 3 real roots (two or three equal roots)
    else if ( ((yN*yN) == (h*h)) ) {
        // two equal roots
        if (h != 0) {
            //del = cbrt(yN/(2*A));
            alpha = beta = xN + del;
            gamma = xN - 2*del;

            roots[0] = alpha;
            roots[1] = beta;
            roots[2] = gamma;

            return *roots;
        }
        // three equal roots
        else {
            alpha = beta = gamma = xN;

            roots[0] = alpha;
            roots[1] = beta;
            roots[2] = gamma;

            return *roots;
        }
    }
    // CASE 3: 3 distinct real roots (yN^2 < h^2)
    else {
        theta = acos(-yN/h) / 3;
        
        alpha = xN + 2*del*cos(theta);
        beta  = xN + 2*del*cos(theta + 2*PI/3);
        gamma = xN + 2*del*cos(theta + 4*PI/3);

        roots[0] = alpha;
        roots[1] = beta;
        roots[2] = gamma;

        return *roots;
    }
    
    return *roots;
}