GPU: give Metal an IEEE-754 binary64 in software

Metal Shading Language has no double at all, and the tracking needs one in a
handful of places: the track parametrisation, the propagator and the material
budget all carry double members and do double arithmetic.

GPUCommonDoubleBinary64.h implements binary64 on top of a 64-bit integer. Round
to nearest even only, with subnormals, infinities and NaNs, and NaN propagation
in the ARM64 order, so an Apple host is a bit-exact reference down to the
payload. Addition, subtraction, multiplication, division and the conversions to
and from float and the 32-bit integers are exact; sin and cos are fdlibm's and
land within 2 ulp of libm. There is no fused multiply-add and no square root.

GPUCommonMath's own double helpers cannot serve the class: SinCosd's OpenCL arm
passes the cosine by pointer where MSL takes a thread reference, so it does not
compile, and Abs<double> is worse because it does -- fabs resolves to
metal::fabs(float) through the implicit conversion, with no error and no warning,
rounding the mantissa away at 17 call sites including the SMatrixGPU pivot
selection and the Kalman guards against Almost0 and Almost1. Both are defined for
Metal in GPUCommonDouble.h instead, and a static_assert there turns the silent
case into a compile error, since metal::fabs is not constant-evaluable.

The Metal entry point aliases the double keyword to the class, so the shared
headers go on saying double and no call site changes. The header refuses to
build anywhere else, since every other backend has a real double.

It costs of the order of a hundred times plain float on an M-series GPU, which
the tracking can afford because double is a small fraction of its floating point
work. In exchange a Metal build reproduces the CPU result bit for bit wherever
it uses only the exact operations.

Author: ktf

DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx

@@ -180,10 +180,10 @@ GPUd() bool TrackParametrizationWithError<value_T>::propagateTo(value_t xk, Trac
   }
   value_t kb = bz * constants::math::B2C;
   // evaluate in double prec.
-  double snpRef0 = linRef0.getSnp(), cspRef0 = gpu::CAMath::Sqrt((1 - snpRef0) * (1 + snpRef0));
-  double snpRef1 = linRef1.getSnp(), cspRef1 = gpu::CAMath::Sqrt((1 - snpRef1) * (1 + snpRef1));
-  double cspRef0Inv = 1 / cspRef0, cspRef1Inv = 1 / cspRef1, cc = cspRef0 + cspRef1, ccInv = 1 / cc, dy2dx = (snpRef0 + snpRef1) * ccInv;
-  double dxccInv = dx * ccInv, hh = dxccInv * cspRef1Inv * (1 + cspRef0 * cspRef1 + snpRef0 * snpRef1), jj = dx * (dy2dx - snpRef1 * cspRef1Inv);
+  double snpRef0 = linRef0.getSnp(), cspRef0 = gpu::CAMath::Sqrt((1.f - snpRef0) * (1.f + snpRef0));
+  double snpRef1 = linRef1.getSnp(), cspRef1 = gpu::CAMath::Sqrt((1.f - snpRef1) * (1.f + snpRef1));
+  double cspRef0Inv = 1.f / cspRef0, cspRef1Inv = 1.f / cspRef1, cc = cspRef0 + cspRef1, ccInv = 1.f / cc, dy2dx = (snpRef0 + snpRef1) * ccInv;
+  double dxccInv = dx * ccInv, hh = dxccInv * cspRef1Inv * (1.f + cspRef0 * cspRef1 + snpRef0 * snpRef1), jj = dx * (dy2dx - snpRef1 * cspRef1Inv);
 
   double f02 = hh * cspRef0Inv;
   double f04 = hh * dxccInv * kb;
@@ -638,7 +638,7 @@ GPUd() bool TrackParametrizationWithError<value_T>::propagateTo(value_t xk, cons
   }
   double r1pr2Inv = 1. / (r1 + r2), r2inv = 1. / r2, r1inv = 1. / r1;
   double dy2dx = (f1 + f2) * r1pr2Inv, dx2r1pr2 = dx * r1pr2Inv;
-  value_t step = (gpu::CAMath::Abs(x2r) < 0.05f) ? dx * gpu::CAMath::Abs(r2 + f2 * dy2dx)                                                   // chord
+  value_t step = (gpu::CAMath::Abs(x2r) < 0.05f) ? value_t(dx * gpu::CAMath::Abs(r2 + f2 * dy2dx))                                          // chord
                                                  : 2.f * gpu::CAMath::ASin(0.5f * dx * gpu::CAMath::Sqrt(1.f + dy2dx * dy2dx) * crv) / crv; // arc
   step *= gpu::CAMath::Sqrt(1.f + this->getTgl() * this->getTgl());
   //
@@ -1086,7 +1086,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

@@ -0,0 +1,49 @@
+// Copyright 2019-2025 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+/// \file GPUCommonDouble.h
+/// \brief Storage for double precision members shared with a device that has none
+
+/// \brief What MSL needs once the double keyword names the emulated type
+
+#ifndef GPUCOMMONDOUBLE_H
+#define GPUCOMMONDOUBLE_H
+
+#include "GPUCommonDef.h"
+#include "GPUCommonMath.h"
+
+#ifdef __METAL__
+
+namespace o2::gpu
+{
+
+static_assert(sizeof(double) == 8, "the emulated double must match the size of a real one");
+static_assert(alignof(double) == 8, "the emulated double must match the alignment of a real one");
+
+// CAMath::Abs deduces its parameter rather than taking a float, so a call on a
+// double picks the primary template, which has no definition. The rest of CAMath
+// takes float and is reached through the implicit conversion.
+template <>
+GPUhdi() constexpr double GPUCommonMath::Abs<double>(double x)
+{
+  return double::fromBits(x.bits() & ~GPUCA_B64_SIGN);
+}
+
+// metal::fabs is not constant-evaluable, so this also fails to compile if the
+// specialisation above is ever dropped and the call falls back to it in float
+static_assert(GPUCommonMath::Abs<double>(GPUdoubleBinary64::fromBits(0xBFF0000000000001ULL)).bits() == 0x3FF0000000000001ULL,
+              "Abs on the emulated double must clear the sign bit and keep every other one");
+
+} // namespace o2::gpu
+
+#endif // __METAL__
+
+#endif // GPUCOMMONDOUBLE_H