GPU: make Metal's cluster widths round as on the CPU

Needs review. On a 30-collision Pb-Pb time frame, 7 of 4,416,506 clusters
differed between the CPU and the Metal backend, each by one step of a
packed width, sigmaTime or sigmaPad. Those widths were all within one
float ulp below a quantisation step: sigma * 32 was 11.4999994,
14.4999992, 29.4999996 or 31.4999988, where the packing rounds at n + 0.5,
so any difference in the last bit decides between n and n + 1.

The CPU build compiles this with -ffast-math. Its arm64 code for
finalize() computes 1 / Q once, the means as products with it, the
variance as one fused multiply-add, fma(SS, 1 / Q, -mean * mean), and an
IEEE sqrt. Metal's fast-math division and sqrt are not correctly rounded.
Do the same operations on Metal explicitly, with precise::divide, fma and
precise::sqrt: the pbpb clusters then match the CPU bit for bit, with no
measurable cost; the rest of the Metal code keeps fast math.

This matches the CPU build that the compiler produces on arm64 with
clang. An x86 build, or another compiler, may choose a different sequence
under -ffast-math. Writing the sequence explicitly for all backends would
make the widths independent of those choices, but can change the CPU
results on some platforms.

Author: ktf

GPU/GPUTracking/TPCClusterFinder/ClusterAccumulator.cxx

@@ -62,13 +62,26 @@ GPUd() void ClusterAccumulator::finalize(const CfChargePos& pos, const Charge q,
 {
   mQtot += q;
 
+#ifdef __METAL__
+  // Metal's fast-math division and sqrt are not correctly rounded, which moves
+  // widths that sit within an ulp of a packing step. Do the operations the
+  // arm64 -ffast-math CPU build does here: reciprocal, products, fused variance.
+  const float invQ = metal::precise::divide(1.f, mQtot);
+  mPadMean *= invQ;
+  mTimeMean *= invQ;
+  const float padMean2 = mPadMean * mPadMean;
+  const float timeMean2 = mTimeMean * mTimeMean;
+  mPadSigma = metal::precise::sqrt(metal::fma(mPadSigma, invQ, -padMean2));
+  mTimeSigma = metal::precise::sqrt(metal::fma(mTimeSigma, invQ, -timeMean2));
+#else
   mPadMean /= mQtot;
   mTimeMean /= mQtot;
   mPadSigma /= mQtot;
   mTimeSigma /= mQtot;
 
   mPadSigma = CAMath::Sqrt(mPadSigma - mPadMean * mPadMean);
   mTimeSigma = CAMath::Sqrt(mTimeSigma - mTimeMean * mTimeMean);
+#endif
 
   Pad pad = pos.pad();
   mPadMean += pad;