math.cpp

void solve_quartic_equation(std::complex<double> x[], std::complex<double> a, std::complex<double> b, std::complex<double> c, std::complex<double> d, std::complex<double> e)
{
    if (!(islessequal(std::abs(std::real(a)), std::numeric_limits<double>::epsilon()) && islessequal(std::abs(std::imag(a)), std::numeric_limits<double>::epsilon())))
    {
    	std::complex<double> alpha = std::complex<double>(-3.0) * std::pow(b, 2) / (std::complex<double>(8.0) * std::pow(a, 2)) + c / a;
    	std::complex<double> beta = std::pow(b, 3) / (std::complex<double>(8.0) * std::pow(a, 3)) - b * c / (std::complex<double>(2.0) * std::pow(a, 2)) + d / a;
    	std::complex<double> gamma = std::complex<double>(-3.0) * std::pow(b, 4) / (std::complex<double>(256.0) * std::pow(a, 4)) + c * std::pow(b, 2) / (std::complex<double>(16.0) * std::pow(a, 3)) - b * d / (std::complex<double>(4.0) * std::pow(a, 2)) + e / a;

    	std::complex<double> P = -std::pow(alpha, 2) / std::complex<double>(12.0) - gamma;
    	std::complex<double> Q = -std::pow(alpha, 3) / std::complex<double>(108.0) + alpha * gamma / std::complex<double>(3.0) - std::pow(beta, 2) / std::complex<double>(8.0);

    	std::complex<double> R = Q / std::complex<double>(2.0) + std::sqrt(std::pow(Q, 2) / std::complex<double>(4.0) + std::pow(P, 3) / std::complex<double>(27.0));
    	std::complex<double> U = std::pow(R, 1.0 / 3.0);

    	std::complex<double> y = std::complex<double>(-5.0 / 6.0) * alpha - U;
    	if (U != std::complex<double>(0.0)) {
    		y += P / (std::complex<double>(3.0) * U);
    	}

    	std::complex<double> W = std::sqrt(alpha + std::complex<double>(2.0) * y);

    	std::complex<double> minus_b_div_4a = -b / (std::complex<double>(4.0) * a);
    	std::complex<double> two_beta_div_w = std::complex<double>(2.0) * beta / W;
    	std::complex<double> three_alpha_plus_two_y = std::complex<double>(3.0) * alpha + std::complex<double>(2.0) * y;

    	x[0] = minus_b_div_4a + (W + std::sqrt(-(three_alpha_plus_two_y + two_beta_div_w))) / std::complex<double>(2.0);
    	x[1] = minus_b_div_4a + (W - std::sqrt(-(three_alpha_plus_two_y + two_beta_div_w))) / std::complex<double>(2.0);
	    x[2] = minus_b_div_4a + (-W + std::sqrt(-(three_alpha_plus_two_y - two_beta_div_w))) / std::complex<double>(2.0);
    	x[3] = minus_b_div_4a + (-W - std::sqrt(-(three_alpha_plus_two_y - two_beta_div_w))) / std::complex<double>(2.0);
    } else {
        solve_cubic(x, b, c, d, e);
        x[3] = x[0];
    }
}