[RFC PATCH] gas model: use virial cubic polynomial form

[Linus Torvalds]

From: Linus Torvalds <torvalds at linux-foundation.org>
Date: Thu, 3 Mar 2016 21:30:28 -0800
Subject: [PATCH] gas model: use virial cubic polynomial form

The "virial" form of the Z compression factor is of the form

   Z = 1.0 + A*p + B*p^2 + C*p^3 + ..

and it's considered the "right" polynomial form to use.  It happens to
also make for one constant less per gas (since the 1.0 can be added
later), and can be used to simplify the expression and avoid a few
floating point operations.

However, in order for that kind of expression simplification to make
sense, we need to make sure that we don't calculate the powers of the
pressure multiple times either, and that means we have to inline all the
actual calculations.

Our compiler options still mean that the generated code isn't optimal,
but that's a separate issue. And it is a lot better than it used to be.

Being clever about this does potentially make the code a tiny bit less
legible, but maybe that's not too bad for something that we'd expect to
not ever touch once we get it right.

Signed-off-by: Linus Torvalds <torvalds at linux-foundation.org>
---

NOTE THE "RFC" part. I'm not sure this is worth it. It came from the G+ 
discussion: the virial form is apparently the preferred form for the Z 
calculation and the difference to the pure cubic least-square is minimal.

And it does generate better code, just because there are fewer constants 
and fewer FP operations and we can delay adding the 1.0 term until after 
the gas mix, rather than doing it for each gas.

At the same time, I'm not exactly claiming that this really helps make the 
code obvious.  Maybe it's premature optimization. And really, considering 
that we don't even inline the "get_o2()" trivial function etc, the code 
generation really isn't very optional.

So for this to make sense, we'd need to feel that performance matters. We 
*do* end up doing the gas use for each plot point, so it's not like this 
function might not get called a lot.

But compiler options will change the code generation even more than this 
ever will, and aggressive compiler optimizations would have made the old 
code better too (ie just inlining and CSE on the pressure power series 
would do good things to the thing without this patch too).

In other words: I think this patch is correct, and it *does* make the code 
faster. But do we even care at this point?

 subsurface-core/gas-model.c | 98 +++++++++++++++++++--------------------------
 1 file changed, 42 insertions(+), 56 deletions(-)

diff --git a/subsurface-core/gas-model.c b/subsurface-core/gas-model.c
index 420233612cb9..81765e003002 100644
--- a/subsurface-core/gas-model.c
+++ b/subsurface-core/gas-model.c
@@ -4,72 +4,58 @@
 #include <stdlib.h>
 #include "dive.h"
 
-/*
- * Generic cubic polynomial
- */
-static double cubic(double bar, const double coefficient[])
-{
-	double	x0 = 1.0,
-		x1 = bar,
-		x2 = x1*x1,
-		x3 = x2*x1;
-
-	return	x0 * coefficient[0] +
-		x1 * coefficient[1] +
-		x2 * coefficient[2] +
-		x3 * coefficient[3];
-}
+/* "Virial minus one" - the virial cubic form without the initial 1.0 */
+#define virial_m1(C, x1, x2, x3) (C[0]*x1+C[1]*x2+C[2]*x3)
 
 /*
- * Cubic least-square coefficients for O2/N2/He based on data from
+ * Cubic virial least-square coefficients for O2/N2/He based on data from
  *
  *    PERRY’S CHEMICAL ENGINEERS’ HANDBOOK SEVENTH EDITION
  *
  * with the lookup and curve fitting by Lubomir.
  *
+ * The "virial" form of the compression factor polynomial is
+ *
+ *      Z = 1.0 + C[0]*P + C[1]*P^2 + C[2]*P^3 ...
+ *
+ * and these tables do not contain the initial 1.0 term.
+ *
  * NOTE! Helium coefficients are a linear mix operation between the
  * 323K and one for 273K isotherms, to make everything be at 300K.
  */
-static const double o2_coefficients[4] = {
-	+1.00117935180448264158,
-	-0.00074149079841471884,
-	+0.00000291901111247509,
-	-0.00000000162092217461
-};
-
-static const double n2_coefficients[4] = {
-	+1.00030344355797817778,
-	-0.00022528077251905598,
-	+0.00000295430303276288,
-	-0.00000000210649996337
-};
-
-static const double he_coefficients[4] = {
-	+1.00000137322788319261,
-	+0.000488393024887620665,
-	-0.000000054210166760015,
-	+0.000000000010908069275
-};
-
-static double o2_compressibility_factor(double bar) { return cubic(bar, o2_coefficients); }
-static double n2_compressibility_factor(double bar) { return cubic(bar, n2_coefficients); }
-static double he_compressibility_factor(double bar) { return cubic(bar, he_coefficients); }
-
-/*
- * We end up using a simple linear mix of the gas-specific functions.
- */
 double gas_compressibility_factor(struct gasmix *gas, double bar)
 {
-	double o2, n2, he;	// Z factors
-	double of, nf, hf;	// gas fractions ("partial pressures")
-
-	o2 = o2_compressibility_factor(bar);
-	n2 = n2_compressibility_factor(bar);
-	he = he_compressibility_factor(bar);
-
-	of = get_o2(gas) / 1000.0;
-	hf = get_he(gas) / 1000.0;
-	nf = 1.0 - of - nf;
-
-	return o2*of + n2*nf + he*hf;
+	static const double o2_coefficients[3] = {
+		-0.00071809207370164567,
+		+0.00000281852572807643,
+		-0.00000000150290620491
+	};
+	static const double n2_coefficients[3] = {
+		-0.00021926035329221337,
+		+0.00000292844845531647,
+		-0.00000000207613482075
+	};
+	static const double he_coefficients[3] = {
+		+0.00047961098687979363,
+		-0.00000004077670019935,
+		+0.00000000000077707035
+	};
+	double o2, he;
+	double x1, x2, x3;
+	double Z;
+
+	o2 = get_o2(gas) / 1000.0;
+	he = get_he(gas) / 1000.0;
+
+	x1 = bar; x2 = x1*x1; x3 = x2*x1;
+
+	Z = virial_m1(o2_coefficients, x1, x2, x3) * o2 +
+	    virial_m1(he_coefficients, x1, x2, x3) * he +
+	    virial_m1(n2_coefficients, x1, x2, x3) * (1.0 - o2 - he);
+
+	/*
+	 * We add the 1.0 at the very end - the linear mixing of the
+	 * three 1.0 terms is still 1.0 regardless of the gas mix.
+	 */
+	return Z + 1.0;
 }
-- 
2.7.2.334.g7c0da37