GPU: add an emulated binary64 and a selector for GPUdoubleCalc

GPUdoubleCalc now names one of four representations, chosen by GPUCA_DOUBLECALC:
hardware double, plain float, the compensated two-float type, or a full IEEE-754
binary64 in software. The default is unchanged, double everywhere except Metal
and the two-float type there, so this is a no-op for every existing build. It
replaces GPUCA_FORCE_DOUBLECALC and GPUCA_FORCE_FLOATCALC, which could only
express two of the four.

The binary64 emulation is round to nearest even with subnormals, infinities and
NaNs, and correctly rounded division; NaN propagation follows the ARM64 order so
an Apple host is a reference down to the payload. It is not for production: over
the whole propagateTo kernel on an M1 Max it costs of the order of a hundred
times plain float, against roughly two for the two-float type. It is here because it
reproduces the CPU result exactly -- every covariance element bit-identical --
which turns a disagreement between Metal and the CPU into a search over code
rather than over numerics.

Measured on that kernel, 524288 tracks with strongly correlated covariances,
median of ten runs, error as |dC_ij| / sqrt(C_ii C_jj) against the CPU double
result: plain float 0.108 ns and 1.2e-5, two-float 0.256 ns and 5.6e-7, binary64
12.2 ns and zero. Only the binary64 kernel is unstable run to run, between 9.1
and 12.3 ns; the others vary by a few percent. For scale, the CPU double result
is itself 1.1e-6 from a cancellation-free reference, so the two-float type is
already at the noise floor of the formula.

For scale in the other direction, the same propagateTo on this machine's CPU,
where double is native, costs 23.5 ns in float and 27.1 ns in double, so real
hardware double is a 1.15x proposition. None of this arithmetic would be needed
if Metal had it.

Validated against hardware double on 9M operand pairs per operation, over uniform
bit patterns, near-equal exponents, the subnormal band, sparse mantissas for the
exact and halfway cases, the special values and track-like magnitudes, plus 4M
conversions each way: no mismatches.

The arithmetic must stay out of line on Metal. Inlined into a kernel it drops the
occupancy and runs two to five times slower than the call.

Author: ktf

DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx

@@ -1087,7 +1087,7 @@ GPUd() auto TrackParametrizationWithError<value_T>::getPredictedChi2(const value
   auto chi2 = (d * (szz * d - sdz * z) + z * (sdd * z - d * sdz)) / det;
   if (chi2 < 0.) {
 #ifndef GPUCA_ALIGPUCODE
-    LOGP(warning, "Negative chi2={}, Cluster: {} {} {} Dy:{} Dz:{} | sdd:{} sdz:{} szz:{} det:{}", chi2, cov[0], cov[1], cov[2], d, z, sdd, sdz, szz, det);
+    LOGP(warning, "Negative chi2={}, Cluster: {} {} {} Dy:{} Dz:{} | sdd:{} sdz:{} szz:{} det:{}", double(chi2), cov[0], cov[1], cov[2], d, z, double(sdd), double(sdz), double(szz), double(det));
     LOGP(warning, "Track: {}", asString());
 #endif
   }

GPU/Common/GPUCommonDouble.h

@@ -17,6 +17,7 @@
 
 #include "GPUCommonDef.h"
 #include "GPUCommonMath.h"
+#include "GPUCommonDoubleBinary64.h"
 
 #ifndef GPUCA_GPUCODE_DEVICE
 #include <cstdint>
@@ -194,27 +195,56 @@ GPUdi() threadgroup float& operator+=(threadgroup float& a, GPUdoubleCalcImpl b)
 GPUdi() float& operator+=(float& a, GPUdoubleCalcImpl b) { return a = (float)(GPUdoubleCalcImpl(a) + b); }
 #endif
 
-// CAMath::Abs is a template that deduces its parameter, so a call on this type
-// picks the primary template, which has no definition. The rest of CAMath takes
-// float and is reached through the implicit conversion.
+// CAMath::Abs is a template that deduces its parameter, so a call on one of these
+// types picks the primary template, which has no definition. The rest of CAMath
+// takes float and is reached through the implicit conversion.
 template <>
 GPUhdi() GPUdoubleCalcImpl GPUCommonMath::Abs<GPUdoubleCalcImpl>(GPUdoubleCalcImpl x)
 {
   return (float)x < 0.f ? -x : x;
 }
 
-// GPUCA_FORCE_DOUBLECALC lets a host test exercise the Metal representation and
-// compare it against the double one.
-#if defined(__METAL__) && defined(__FAST_MATH__)
-// Fast math reassociates the compensation terms away: the two-float type would
-// then cost 1.5x for the accuracy of a plain float.
+template <>
+GPUhdi() GPUdoubleBinary64 GPUCommonMath::Abs<GPUdoubleBinary64>(GPUdoubleBinary64 x)
+{
+  return GPUdoubleBinary64::fromBits(x.bits() & ~GPUCA_B64_SIGN);
+}
+
+// What GPUdoubleCalc is. Define GPUCA_DOUBLECALC to one of these to override the
+// default, which is hardware double everywhere except Metal, and the compensated
+// two-float type there. A host build can be set to any of them to compare the
+// representations against each other; a Metal build can be set to BINARY64 to
+// reproduce the CPU result bit for bit, at two orders of magnitude the cost.
+#define GPUCA_DOUBLECALC_DOUBLE 1
+#define GPUCA_DOUBLECALC_FLOAT 2
+#define GPUCA_DOUBLECALC_TWOFLOAT 3
+#define GPUCA_DOUBLECALC_BINARY64 4
+
+#ifndef GPUCA_DOUBLECALC
+  #if defined(__METAL__) && defined(__FAST_MATH__)
+    // Fast math reassociates the compensation terms away, so the two-float type
+    // would cost 1.5x for the accuracy of a plain float.
+    #define GPUCA_DOUBLECALC GPUCA_DOUBLECALC_FLOAT
+  #elif defined(__METAL__)
+    #define GPUCA_DOUBLECALC GPUCA_DOUBLECALC_TWOFLOAT
+  #else
+    #define GPUCA_DOUBLECALC GPUCA_DOUBLECALC_DOUBLE
+  #endif
+#endif
+
+#if GPUCA_DOUBLECALC == GPUCA_DOUBLECALC_DOUBLE
+  #ifdef __METAL__
+    #error "MSL has no double; GPUCA_DOUBLECALC_DOUBLE cannot be selected for Metal"
+  #endif
+typedef double GPUdoubleCalc;
+#elif GPUCA_DOUBLECALC == GPUCA_DOUBLECALC_FLOAT
 typedef float GPUdoubleCalc;
-#elif defined(__METAL__) || defined(GPUCA_FORCE_DOUBLECALC)
+#elif GPUCA_DOUBLECALC == GPUCA_DOUBLECALC_TWOFLOAT
 typedef GPUdoubleCalcImpl GPUdoubleCalc;
-#elif defined(GPUCA_FORCE_FLOATCALC) // for the host test only, to show what plain float would cost
-typedef float GPUdoubleCalc;
+#elif GPUCA_DOUBLECALC == GPUCA_DOUBLECALC_BINARY64
+typedef GPUdoubleBinary64 GPUdoubleCalc;
 #else
-typedef double GPUdoubleCalc;
+  #error "Invalid setting for GPUCA_DOUBLECALC"
 #endif
 
 } // namespace o2::gpu