Skip to content

Instantly share code, notes, and snippets.

@vrbadev

vrbadev/quartic_roots.cpp

Last active October 10, 2025 13:29
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save vrbadev/ef8d561f692b431e858d82fdc7c094e2 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save vrbadev/ef8d561f692b431e858d82fdc7c094e2 to your computer and use it in GitHub Desktop.
Download ZIP
Fast and simple quartic polynomial roots calculation implemented in plain C++
Raw
quartic_roots.cpp
template<typename T=double>
std::array<std::complex<T>, 4> solve_quartic_roots(std::complex<T> c4, std::complex<T> c3, std::complex<T> c2, std::complex<T> c1, std::complex<T> c0) const {
c3 /= c4; c2 /= c4; c1 /= c4; c0 /= c4;
const T eps = std::numeric_limits<T>::epsilon();
const T inv2 = 0.5, inv3 = 1.0/3.0, inv27 = 1.0/27.0;
auto a4q = 0.25 * c3;
auto a4q2 = a4q * a4q;
auto a4q4 = a4q2 * a4q2;
auto p = 3.0 * a4q2 - inv2 * c2;
auto q = c3 * a4q2 - c2 * a4q + inv2 * c1;
auto r = 3.0 * a4q4 - c2 * a4q2 + c1 * a4q - c0;
auto o = p * r - inv2 * q * q;
auto p_sq = p * p;
auto p_cu = p_sq * p;
auto p2 = (3.0 * r - p_sq) * inv3;
auto q2 = (2.0 * p_cu - 9.0 * p * r + 27.0 * o) * inv27;
auto q2h = -inv2 * q2;
auto p2_third = p2 * inv3;
auto d = std::sqrt(q2h * q2h + p2_third * p2_third * p2_third);
auto w = q2h + d;
auto w_mag = std::abs(w);
auto w_ang = std::arg(w);
auto u = std::polar(std::pow(w_mag, inv3), w_ang * inv3);
if (std::abs(u) < eps) {
auto w2 = q2h - d;
u = std::polar(std::pow(std::abs(w2), inv3), std::arg(w2) * inv3);
}
auto v = -p2 / (3.0 * u);
auto z0 = u + v - p * inv3;
auto s = std::sqrt(2.0 * (p + z0));
auto t = (std::abs(s) < eps) ? (z0 * z0 + r) : (-q / s);
auto sh = 0.5 * s;
auto base = sh * sh - z0;
auto d1 = std::sqrt(base - t);
auto d2 = std::sqrt(base + t);
return {-sh - d1 - a4q, -sh + d1 - a4q, sh - d2 - a4q, sh + d2 - a4q};
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment