Still Faster, After 30 Years

Back in the thermodynamic code I keep writing lately, one tiny problem shows up in every single iteration of a flash calculation. You need to solve a cubic equation in the compressibility factor Z, the heart of any cubic equation of state (Peng-Robinson, Soave-Redlich-Kwong, you name it).

A cubic has three roots. Depending on the conditions, you take one or another depending on your goal, if you are looking at the liquid, the vapor or the root minimizing the energy. A simulator calls this solver millions of times per full plant simulation, even more in reservoir simulations. Every nanosecond you save matters.

So, this weekend, I got curious. The thermolib code I have been using for years ships a tight little Newton solver for this cubic, written in Fortran in the 1990s. Is it still the fastest thing in the room, or has an analytical approach caught up in thirty years of compilers and CPUs?

The problem

The cubic we need to solve has a known form:

$Z³ - X₂·Z² + X₁·Z - X₀ = 0$

where $X₀$, $X₁$, $X₂$ come from the reduced EOS parameters $A$, $B$ and the EOS-specific constants (different for PR and SRK, the same for both solvers in my comparison). The physical root is always a bit above $B$ for the liquid phase, or a bit above $B + 1$ for the vapor phase. That is the key domain fact.

Two approaches

The old Fortran way: Newton with Halley (cubic) corrections.

Given the polynomial $f(Z) = Z^3 - X_2 Z^2 + X_1 Z - X_0$ and its first two derivatives $f'$ and $f''$, pick a starting value $Z_0$ by checking the sign of $f$ at the inflection point $Z = X_2 / 3$. Take $Z_0 = B$ on the liquid side ($f > 0$), or $Z_0 = B + 1$ on the vapor side.

Then iterate, with a third-order correction:

$$\delta_n = \frac{f(Z_n)}{f'(Z_n)} \left(1 + \frac{f(Z_n) \, f''(Z_n)}{f'(Z_n)^2}\right), \qquad Z_{n+1} = Z_n - \delta_n$$

Stop when $|\delta_n| < 10^{-7}$.

Two or three iterations and it is done. If a second root is needed, factor the first one out and apply one more Newton step on the remaining quadratic.

The old Cardano/mathematical way: closed form, no iteration. The formula is from 1545, that is not a typo. If you are interested with history, the Wikipedia page about the cubic equation is a rabbit hole you can spent a couple of hours into and then come back here.

Rewrite the cubic as $Z^3 + p Z^2 + q Z + r = 0$ with $p = -X_2$, $q = X_1$, $r = -X_0$, then substitute $Z = t - p/3$ to kill the quadratic term:

$$t^3 + P\,t + Q = 0, \qquad P = q - \frac{p^2}{3}, \qquad Q = \frac{2 p^3}{27} - \frac{p q}{3} + r$$

The discriminant decides which branch to take:

$$\Delta = \left(\frac{Q}{2}\right)^2 + \left(\frac{P}{3}\right)^3$$

If $\Delta > 0$, one real root, given by a sum of two real cube roots:

$$t = \sqrt[3]{-\tfrac{Q}{2} + \sqrt{\Delta}} + \sqrt[3]{-\tfrac{Q}{2} - \sqrt{\Delta}}$$

If $\Delta \leq 0$, three real roots, given by the trigonometric form:

$$t_k = 2 \sqrt{-\tfrac{P}{3}} \, \cos!\left(\frac{1}{3} \arccos!\left(\frac{-Q/2}{\sqrt{-(P/3)^3}}\right) - \frac{2 \pi k}{3}\right), \qquad k = 0, 1, 2$$

Recover $Z = t - p/3$ and pick the physical root.

No loop and no tolerance, a single pass through a 480-year-old formula. On paper, this should be hard to beat.

The benchmark

I reimplemented both solvers in Fortran, in a small standalone file, and tested five cases: single-root, two-phase and dense-liquid regimes. Both agree to machine precision, max error around 2·10⁻¹⁶.

Then I timed them with my benchmarking tools:

The old Fortran Newton wins every case. Between 1.22× and 3.16× faster. Cardano stays flat at around 46 ns per call, Newton varies from 14 to 38 ns depending on how close the starting guess was.

Why the old code still wins

Cardano pays for a cbrtf, a sqrt, an acos and a cos on every call. Those are not cheap on a modern CPU. Each takes tens of nanoseconds of hardware work, and the analytical formula needs all of them. There is no shortcut, because the formula is general and does not know anything about where the root should be.

Newton pays for nothing exotic, just a handful of multiplies and adds per iteration. And it converges in two or three iterations, because the starting guess is already close. The author of the Fortran code knew that the physical root is always near $B$ or near $B + 1$. That one piece of domain knowledge is worth a factor of three.

This is the whole point. Cardano is general and cannot skip anything. The Fortran solver is specialised and skips a lot because the author knew where the root lives.

Not a coincidence

The author of that Fortran Newton solver is Michael L. Michelsen, the one from Thermodynamic Models: Fundamentals & Computational Aspects, and the one whose department I was lucky to do my PhD in, in Denmark.

Thirty years later, with a modern Fortran compiler and a modern 2025 (AMD) CPU, I would have expected the general and analytical code to have caught up (maybe through vectorization, better tables, whatever).

It has not. The old code still wins, because it embbeds the domain knowledge of the expert in it.

PS: By the way, a good start Newton/Halley's method implementation can be found on the CD accompanying the book of Michelsen and Mollerup.