From 20a92b7317ce271176d16fd8ffbab37dcfb4e7a3 Mon Sep 17 00:00:00 2001 From: Giulio Eulisse <10544+ktf@users.noreply.github.com> Date: Tue, 29 Sep 2026 11:40:54 +0200 Subject: [PATCH 1/2] 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 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. --- .../src/TrackParametrizationWithError.cxx | 12 +- GPU/Common/GPUCommonDouble.h | 49 ++ GPU/Common/GPUCommonDoubleBinary64.h | 589 ++++++++++++++++++ GPU/Common/GPUCommonMath.h | 4 + 4 files changed, 648 insertions(+), 6 deletions(-) create mode 100644 GPU/Common/GPUCommonDouble.h create mode 100644 GPU/Common/GPUCommonDoubleBinary64.h diff --git a/DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx b/DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx index 748cb47094d26..946492a4e2a54 100644 --- a/DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx +++ b/DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx @@ -180,10 +180,10 @@ GPUd() bool TrackParametrizationWithError::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::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::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 } diff --git a/GPU/Common/GPUCommonDouble.h b/GPU/Common/GPUCommonDouble.h new file mode 100644 index 0000000000000..2f7108b409a7e --- /dev/null +++ b/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 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(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 diff --git a/GPU/Common/GPUCommonDoubleBinary64.h b/GPU/Common/GPUCommonDoubleBinary64.h new file mode 100644 index 0000000000000..572e17797c0a2 --- /dev/null +++ b/GPU/Common/GPUCommonDoubleBinary64.h @@ -0,0 +1,589 @@ +// Copyright 2019-2026 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 GPUCommonDoubleBinary64.h +/// \brief IEEE-754 binary64 in software, for a device that has none + +#ifndef GPUCOMMONDOUBLEBINARY64_H +#define GPUCOMMONDOUBLEBINARY64_H + +#if !defined(__METAL__) && !defined(GPUCA_B64_HOST_REFERENCE) +#error "the emulated binary64 is for Metal, which has no double of its own; every other backend has a real one. Define GPUCA_B64_HOST_REFERENCE to build it on the host as a test reference." +#endif + +#include "GPUCommonDef.h" +#ifndef GPUCA_GPUCODE_DEVICE +#include +#endif + +// IEEE-754 binary64 in software, for a device that has no double at all. Round to +// nearest even only, with subnormals, infinities and NaNs; NaN propagation follows +// the ARM64 order, so an Apple host is a bit-exact reference down to the payload. +// There is no fused multiply-add, no square root and no other rounding mode. +// +// Addition, subtraction, multiplication, division and the conversions to and from +// float and the 32-bit integers are exact, so a Metal build reproduces the CPU +// result bit for bit wherever it uses only those. sin and cos come from fdlibm and +// land within 2 ulp of libm rather than matching it. +// +// 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. + +namespace o2::gpu +{ + +namespace binary64_detail +{ + +#ifdef __METAL__ +typedef ulong u64; +typedef long i64; +typedef uint u32; +#else +typedef uint64_t u64; +typedef int64_t i64; +typedef uint32_t u32; +#endif + +#define GPUCA_B64_ALWAYS inline __attribute__((always_inline)) +#ifdef __METAL__ +// The arithmetic has to stay out of line: inlined into a kernel it drops the +// occupancy (maxTotalThreadsPerThreadgroup 832 -> 384) and runs two to five +// times slower than the call. +#define GPUCA_B64_OP __attribute__((noinline)) +#else +#define GPUCA_B64_OP inline +#endif + +#ifdef __METAL__ +// metal::clz is not constant-evaluable, which would keep the whole class out of +// constant expressions; the builtin folds to the same ctlz at run time. +GPUCA_B64_ALWAYS constexpr int32_t clz64(u64 x) { return x ? __builtin_clzl(x) : 64; } +GPUCA_B64_ALWAYS constexpr int32_t clz32(u32 x) { return x ? __builtin_clz(x) : 32; } +GPUCA_B64_ALWAYS constexpr u32 asu32(float f) { return as_type(f); } +GPUCA_B64_ALWAYS constexpr float asf32(u32 u) { return as_type(u); } +#else +GPUCA_B64_ALWAYS constexpr int32_t clz64(u64 x) { return x ? __builtin_clzll(x) : 64; } +GPUCA_B64_ALWAYS constexpr int32_t clz32(u32 x) { return x ? __builtin_clz(x) : 32; } + +GPUCA_B64_ALWAYS constexpr u32 asu32(float f) { return __builtin_bit_cast(uint32_t, f); } +GPUCA_B64_ALWAYS constexpr float asf32(u32 u) { return __builtin_bit_cast(float, u); } +#endif + +// The high half of a 64x64 product, from 32-bit partial products. Neither +// metal::mulhi nor __int128 can appear in a constant expression, and the +// program-scope constants are constexpr. +GPUCA_B64_ALWAYS constexpr u64 mulhiu(u64 a, u64 b) +{ + const u64 al = a & 0xffffffffu, ah = a >> 32, bl = b & 0xffffffffu, bh = b >> 32; + const u64 ll = al * bl, lh = al * bh, hl = ah * bl, hh = ah * bh; + const u64 mid = (ll >> 32) + (lh & 0xffffffffu) + (hl & 0xffffffffu); + return hh + (lh >> 32) + (hl >> 32) + (mid >> 32); +} +GPUCA_B64_ALWAYS constexpr i64 mulhis(i64 a, i64 b) +{ + u64 h = mulhiu((u64)a, (u64)b); + h -= (a < 0) ? (u64)b : 0ULL; + h -= (b < 0) ? (u64)a : 0ULL; + return (i64)h; +} + +#define GPUCA_B64_SIGN 0x8000000000000000ULL +#define GPUCA_B64_FRAC 0x000fffffffffffffULL +#define GPUCA_B64_IMPL 0x0010000000000000ULL +#define GPUCA_B64_QUIET 0x0008000000000000ULL +#define GPUCA_B64_INF 0x7ff0000000000000ULL +#define GPUCA_B64_DNAN 0x7ff8000000000000ULL // ARM default NaN + +// right shift keeping a sticky bit; any count >= 0, no shift count ever reaches 64 +GPUCA_B64_ALWAYS constexpr u64 shrJam(u64 m, int32_t s) +{ + const int32_t sc = s > 63 ? 63 : s; + const u64 r = m >> sc; + const u64 lost = (m << (63 - sc)) << 1; + const bool big = s > 63; + const u64 rr = big ? 0ULL : r; + const u64 st = big ? m : lost; + return rr | (st != 0 ? 1ULL : 0ULL); +} + +// m: leading bit at 62 for a normal result, 10 guard bits below the 53-bit significand, +// value = m * 2^(e - 1085) with e the biased exponent. Packing (e-1)<<52 + m lets a +// rounding carry ripple into the exponent. +GPUCA_B64_ALWAYS constexpr u64 roundPack(u64 sign, int32_t e, u64 m) +{ + if (e >= 0x7ff) { + return sign | GPUCA_B64_INF; + } + if (e <= 0) { + m = shrJam(m, 1 - e); + e = 1; + } + const u64 r = m & 0x3ffULL; + m >>= 10; + if (r > 0x200ULL || (r == 0x200ULL && (m & 1ULL))) { + ++m; + } + return sign | (((u64)(e - 1) << 52) + m); +} + +// an sNaN operand wins over a qNaN one, among equals the first wins, the result is quietened +GPUCA_B64_ALWAYS constexpr u64 propNaN(u64 a, u64 b, bool an, bool bn) +{ + const bool as = an && !(a & GPUCA_B64_QUIET), bs = bn && !(b & GPUCA_B64_QUIET); + if (as) { + return a | GPUCA_B64_QUIET; + } + if (bs) { + return b | GPUCA_B64_QUIET; + } + return an ? a : b; +} + +GPUCA_B64_OP constexpr u64 addsub(u64 a, u64 b0, bool neg) +{ + const u64 b = neg ? (b0 ^ GPUCA_B64_SIGN) : b0; + const bool sw = (a & ~GPUCA_B64_SIGN) < (b & ~GPUCA_B64_SIGN); + const u64 x = sw ? b : a, y = sw ? a : b; // |x| >= |y| + int32_t ex = (int32_t)((x >> 52) & 0x7ff), ey = (int32_t)((y >> 52) & 0x7ff); + u64 mx = x & GPUCA_B64_FRAC, my = y & GPUCA_B64_FRAC; + if (ex == 0x7ff) { // y can only be inf/NaN if x is too + const bool an = ((a >> 52) & 0x7ff) == 0x7ff && (a & GPUCA_B64_FRAC) != 0, bn = ((b0 >> 52) & 0x7ff) == 0x7ff && (b0 & GPUCA_B64_FRAC) != 0; + if (an || bn) { + return propNaN(a, b0, an, bn); + } + if (ey == 0x7ff) { + return ((x ^ y) & GPUCA_B64_SIGN) ? GPUCA_B64_DNAN : x; + } + return x; + } + const u64 sx = x & GPUCA_B64_SIGN; + const bool sub = ((x ^ y) >> 63) != 0; + mx = (ex ? (mx | GPUCA_B64_IMPL) : mx) << 10; + ex = ex ? ex : 1; + my = (ey ? (my | GPUCA_B64_IMPL) : my) << 10; + ey = ey ? ey : 1; + my = shrJam(my, ex - ey); + const u64 s = sub ? mx - my : mx + my; + const bool carry = (s >> 63) != 0; + const int32_t lz = clz64(s) - 1; // -1 on carry, 63 on zero + const u64 sN = carry ? ((s >> 1) | (s & 1ULL)) : (s << (lz & 63)); + const int32_t eN = carry ? ex + 1 : ex - lz; + const bool zero = s == 0; + return roundPack((zero && sub) ? 0ULL : sx, zero ? -1 : eN, sN); +} + +GPUCA_B64_OP constexpr u64 mul(u64 a, u64 b) +{ + const u64 sign = (a ^ b) & GPUCA_B64_SIGN; + int32_t ea = (int32_t)((a >> 52) & 0x7ff), eb = (int32_t)((b >> 52) & 0x7ff); + u64 ma = a & GPUCA_B64_FRAC, mb = b & GPUCA_B64_FRAC; + if (ea == 0x7ff || eb == 0x7ff) { + const bool an = ea == 0x7ff && ma != 0, bn = eb == 0x7ff && mb != 0; + if (an || bn) { + return propNaN(a, b, an, bn); + } + if ((ea == 0 && ma == 0) || (eb == 0 && mb == 0)) { + return GPUCA_B64_DNAN; // inf * 0 + } + return sign | GPUCA_B64_INF; + } + if (ea == 0) { + if (ma == 0) { + return sign; + } + const int32_t lz = clz64(ma) - 11; + ma <<= lz; + ea = 1 - lz; + } else { + ma |= GPUCA_B64_IMPL; + } + if (eb == 0) { + if (mb == 0) { + return sign; + } + const int32_t lz = clz64(mb) - 11; + mb <<= lz; + eb = 1 - lz; + } else { + mb |= GPUCA_B64_IMPL; + } + const u64 lo = ma * mb, hi = mulhiu(ma, mb); // 106-bit product in [2^104, 2^106) + const bool top = (hi & (1ULL << 41)) != 0; + const u64 m1 = (hi << 21) | (lo >> 43), m0 = (hi << 22) | (lo >> 42); + const u64 st = top ? (lo & ((1ULL << 43) - 1)) : (lo & ((1ULL << 42) - 1)); + const u64 m = (top ? m1 : m0) | (st != 0 ? 1ULL : 0ULL); + return roundPack(sign, ea + eb - 1023 + (top ? 1 : 0), m); +} + +GPUCA_B64_OP constexpr u64 div(u64 a, u64 b) +{ + const u64 sign = (a ^ b) & GPUCA_B64_SIGN; + int32_t ea = (int32_t)((a >> 52) & 0x7ff), eb = (int32_t)((b >> 52) & 0x7ff); + u64 ma = a & GPUCA_B64_FRAC, mb = b & GPUCA_B64_FRAC; + if (ea == 0x7ff || eb == 0x7ff) { + const bool an = ea == 0x7ff && ma != 0, bn = eb == 0x7ff && mb != 0; + if (an || bn) { + return propNaN(a, b, an, bn); + } + if (ea == 0x7ff && eb == 0x7ff) { + return GPUCA_B64_DNAN; // inf / inf + } + return ea == 0x7ff ? (sign | GPUCA_B64_INF) : sign; // inf / x, x / inf + } + if (eb == 0 && mb == 0) { + return (ea == 0 && ma == 0) ? GPUCA_B64_DNAN : (sign | GPUCA_B64_INF); // x / 0 + } + if (ea == 0) { + if (ma == 0) { + return sign; + } + const int32_t lz = clz64(ma) - 11; + ma <<= lz; + ea = 1 - lz; + } else { + ma |= GPUCA_B64_IMPL; + } + if (eb == 0) { + const int32_t lz = clz64(mb) - 11; + mb <<= lz; + eb = 1 - lz; + } else { + mb |= GPUCA_B64_IMPL; + } + const bool lt = ma < mb; + const int32_t e = ea - eb + 1023 - (lt ? 1 : 0); + const u64 A2 = lt ? (ma << 1) : ma; // A2 / mb in [1, 2) + // reciprocal R ~ 2^114 / mb in (2^61, 2^62], seeded from a float division on the top 24 bits + const u32 rb = asu32(1.0f / (float)(u32)(mb >> 29)); + u64 R = (u64)((rb & 0x7fffffu) | 0x800000u) << ((int32_t)((rb >> 23) & 0xff) - 127 + 62); + const u64 Bn = mb << 11; + for (int32_t it = 0; it < 2; ++it) { + const i64 E = (i64)(1ULL << 61) - (i64)mulhiu(Bn, R); + R = (u64)((i64)R + mulhis((i64)R, E * 8)); // not E << 3: shifting a negative value is not a constant expression + } + R = R > (1ULL << 62) ? (1ULL << 62) : R; + // Q ~ A2 * 2^62 / mb, then the exact remainder, which fits in a signed 64-bit + // word because Q is within a few units + u64 Q = mulhiu(A2 << 10, (R << 2) - 1); + i64 rem = (i64)((A2 << 62) - Q * mb); + const i64 adj = mulhis(rem, (i64)R) >> 50; + Q = (u64)((i64)Q + adj); + rem -= adj * (i64)mb; + // after adj the remainder is within one divisor of [0, mb): one predicated step each way + const bool ng = rem < 0; + Q = ng ? Q - 1 : Q; + rem = ng ? rem + (i64)mb : rem; + const bool bg = rem >= (i64)mb; + Q = bg ? Q + 1 : Q; + rem = bg ? rem - (i64)mb : rem; + return roundPack(sign, e, Q | (rem != 0 ? 1ULL : 0ULL)); +} + +// an int32 or a uint32 always fits the 53-bit significand, so these are exact +GPUCA_B64_ALWAYS constexpr u64 fromU32(u32 x) +{ + if (x == 0) { + return 0ULL; + } + const int32_t lz = clz32(x); + return ((u64)(31 - lz + 1023) << 52) | (((u64)x << (21 + lz)) & GPUCA_B64_FRAC); +} + +GPUCA_B64_ALWAYS constexpr u64 fromI32(int32_t x) +{ + return fromU32(x < 0 ? (u32)(-(i64)x) : (u32)x) | (x < 0 ? GPUCA_B64_SIGN : 0ULL); +} + +GPUCA_B64_ALWAYS constexpr u64 fromFloat(float f) +{ + const u32 u = asu32(f); + const u64 sign = (u64)(u & 0x80000000u) << 32; + int32_t e = (int32_t)((u >> 23) & 0xff); + u32 m = u & 0x7fffffu; + if (e == 0xff) { + return sign | GPUCA_B64_INF | ((u64)m << 29) | (m ? GPUCA_B64_QUIET : 0ULL); + } + if (e == 0) { + if (m == 0) { + return sign; + } + const int32_t lz = clz32(m) - 8; + m <<= lz; + e = 1 - lz; + } + return sign | ((u64)(e - 127 + 1023) << 52) | ((u64)(m & 0x7fffffu) << 29); +} + +// binary64 -> binary32, round to nearest even, subnormals, inf, NaN (payload kept, quietened) +GPUCA_B64_OP constexpr float toFloat(u64 d) +{ + const u32 sign = (u32)(d >> 32) & 0x80000000u; + const int32_t be = (int32_t)((d >> 52) & 0x7ff); + const u32 man = (u32)((d & GPUCA_B64_FRAC) >> 29); + const u32 drop = (u32)d & 0x1fffffffu; + u32 bits = 0; + if (be == 0x7ff) { + bits = sign | 0x7f800000u | man | ((d & GPUCA_B64_FRAC) ? 0x400000u : 0u); + } else if (be == 0) { + bits = sign; + } else { + const int32_t e = be - 1023 + 127; + if (e >= 0xff) { + bits = sign | 0x7f800000u; + } else if (e > 0) { + bits = sign | ((u32)e << 23) | man; + if ((drop & 0x10000000u) && ((drop & 0x0fffffffu) || (man & 1u))) { + bits += 1u; + } + } else if (e > -24) { + const u32 full = man | 0x800000u; + const u32 sh = (u32)(1 - e); + const u32 lost = full & ((1u << sh) - 1u); + const u32 halfb = 1u << (sh - 1); + u32 sub = full >> sh; + if (lost > halfb || (lost == halfb && ((sub & 1u) || drop))) { + sub += 1u; + } + bits = sign | sub; + } else { + bits = sign; + } + } + return asf32(bits); +} + +} // namespace binary64_detail + +class GPUdoubleBinary64 +{ + public: + GPUdDefault() GPUdoubleBinary64() = default; + GPUdi() constexpr GPUdoubleBinary64(float v) : mBits(binary64_detail::fromFloat(v)) {} + GPUdi() constexpr operator float() const { return binary64_detail::toFloat(mBits); } + + GPUdi() static constexpr GPUdoubleBinary64 fromBits(binary64_detail::u64 b) { return GPUdoubleBinary64(b, FromBits{}); } + GPUdi() constexpr binary64_detail::u64 bits() const { return mBits; } + + GPUdi() constexpr GPUdoubleBinary64 operator-() const { return fromBits(mBits ^ GPUCA_B64_SIGN); } + GPUdi() constexpr GPUdoubleBinary64 operator+(GPUdoubleBinary64 b) const { return fromBits(binary64_detail::addsub(mBits, b.mBits, false)); } + GPUdi() constexpr GPUdoubleBinary64 operator-(GPUdoubleBinary64 b) const { return fromBits(binary64_detail::addsub(mBits, b.mBits, true)); } + GPUdi() constexpr GPUdoubleBinary64 operator*(GPUdoubleBinary64 b) const { return fromBits(binary64_detail::mul(mBits, b.mBits)); } + GPUdi() constexpr GPUdoubleBinary64 operator/(GPUdoubleBinary64 b) const { return fromBits(binary64_detail::div(mBits, b.mBits)); } + + GPUdi() constexpr GPUdoubleBinary64 operator+(float b) const { return *this + GPUdoubleBinary64(b); } + GPUdi() constexpr GPUdoubleBinary64 operator-(float b) const { return *this - GPUdoubleBinary64(b); } + GPUdi() constexpr GPUdoubleBinary64 operator*(float b) const { return *this * GPUdoubleBinary64(b); } + GPUdi() constexpr GPUdoubleBinary64 operator/(float b) const { return *this / GPUdoubleBinary64(b); } +#ifndef __METAL__ + GPUdi() constexpr GPUdoubleBinary64 operator+(double b) const { return *this + (float)b; } + GPUdi() constexpr GPUdoubleBinary64 operator-(double b) const { return *this - (float)b; } + GPUdi() constexpr GPUdoubleBinary64 operator*(double b) const { return *this * (float)b; } + GPUdi() constexpr GPUdoubleBinary64 operator/(double b) const { return *this / (float)b; } +#endif + + // integral operands: an exact match, so `2 * x` does not sit ambiguously between + // converting the int up and converting *this down + GPUdi() constexpr GPUdoubleBinary64 operator+(int32_t b) const { return *this + fromBits(binary64_detail::fromI32(b)); } + GPUdi() constexpr GPUdoubleBinary64 operator-(int32_t b) const { return *this - fromBits(binary64_detail::fromI32(b)); } + GPUdi() constexpr GPUdoubleBinary64 operator*(int32_t b) const { return *this * fromBits(binary64_detail::fromI32(b)); } + GPUdi() constexpr GPUdoubleBinary64 operator/(int32_t b) const { return *this / fromBits(binary64_detail::fromI32(b)); } + GPUdi() constexpr GPUdoubleBinary64 operator+(uint32_t b) const { return *this + fromBits(binary64_detail::fromU32(b)); } + GPUdi() constexpr GPUdoubleBinary64 operator-(uint32_t b) const { return *this - fromBits(binary64_detail::fromU32(b)); } + GPUdi() constexpr GPUdoubleBinary64 operator*(uint32_t b) const { return *this * fromBits(binary64_detail::fromU32(b)); } + GPUdi() constexpr GPUdoubleBinary64 operator/(uint32_t b) const { return *this / fromBits(binary64_detail::fromU32(b)); } +#ifdef __METAL__ + // the same surface for an object that lives in the constant address space, which + // a generic `this` does not reach + GPUdi() constexpr GPUdoubleBinary64(float v) constant : mBits(binary64_detail::fromFloat(v)) {} + GPUdi() constexpr operator float() constant { return binary64_detail::toFloat(mBits); } + GPUdi() constexpr GPUdoubleBinary64 operator+(GPUdoubleBinary64 b) constant { return fromBits(binary64_detail::addsub(mBits, b.mBits, false)); } + GPUdi() constexpr GPUdoubleBinary64 operator-(GPUdoubleBinary64 b) constant { return fromBits(binary64_detail::addsub(mBits, b.mBits, true)); } + GPUdi() constexpr GPUdoubleBinary64 operator*(GPUdoubleBinary64 b) constant { return fromBits(binary64_detail::mul(mBits, b.mBits)); } + GPUdi() constexpr GPUdoubleBinary64 operator/(GPUdoubleBinary64 b) constant { return fromBits(binary64_detail::div(mBits, b.mBits)); } + GPUdi() constexpr GPUdoubleBinary64 operator+(float b) constant { return *this + GPUdoubleBinary64(b); } + GPUdi() constexpr GPUdoubleBinary64 operator-(float b) constant { return *this - GPUdoubleBinary64(b); } + GPUdi() constexpr GPUdoubleBinary64 operator*(float b) constant { return *this * GPUdoubleBinary64(b); } + GPUdi() constexpr GPUdoubleBinary64 operator/(float b) constant { return *this / GPUdoubleBinary64(b); } + GPUdi() constexpr GPUdoubleBinary64 operator*(int32_t b) constant { return *this * fromBits(binary64_detail::fromI32(b)); } + GPUdi() constexpr GPUdoubleBinary64 operator/(int32_t b) constant { return *this / fromBits(binary64_detail::fromI32(b)); } +#endif + + GPUdi() GPUdoubleBinary64& operator+=(GPUdoubleBinary64 b) { return *this = *this + b; } + GPUdi() GPUdoubleBinary64& operator-=(GPUdoubleBinary64 b) { return *this = *this - b; } + GPUdi() GPUdoubleBinary64& operator*=(GPUdoubleBinary64 b) { return *this = *this * b; } + GPUdi() GPUdoubleBinary64& operator/=(GPUdoubleBinary64 b) { return *this = *this / b; } + + private: + struct FromBits { + }; + GPUdi() constexpr GPUdoubleBinary64(binary64_detail::u64 b, FromBits) : mBits(b) {} + + binary64_detail::u64 mBits; +}; + +GPUdi() constexpr GPUdoubleBinary64 operator+(int32_t a, GPUdoubleBinary64 b) { return GPUdoubleBinary64::fromBits(binary64_detail::fromI32(a)) + b; } +GPUdi() constexpr GPUdoubleBinary64 operator-(int32_t a, GPUdoubleBinary64 b) { return GPUdoubleBinary64::fromBits(binary64_detail::fromI32(a)) - b; } +GPUdi() constexpr GPUdoubleBinary64 operator*(int32_t a, GPUdoubleBinary64 b) { return GPUdoubleBinary64::fromBits(binary64_detail::fromI32(a)) * b; } +GPUdi() constexpr GPUdoubleBinary64 operator/(int32_t a, GPUdoubleBinary64 b) { return GPUdoubleBinary64::fromBits(binary64_detail::fromI32(a)) / b; } +GPUdi() constexpr GPUdoubleBinary64 operator+(uint32_t a, GPUdoubleBinary64 b) { return GPUdoubleBinary64::fromBits(binary64_detail::fromU32(a)) + b; } +GPUdi() constexpr GPUdoubleBinary64 operator-(uint32_t a, GPUdoubleBinary64 b) { return GPUdoubleBinary64::fromBits(binary64_detail::fromU32(a)) - b; } +GPUdi() constexpr GPUdoubleBinary64 operator*(uint32_t a, GPUdoubleBinary64 b) { return GPUdoubleBinary64::fromBits(binary64_detail::fromU32(a)) * b; } +GPUdi() constexpr GPUdoubleBinary64 operator/(uint32_t a, GPUdoubleBinary64 b) { return GPUdoubleBinary64::fromBits(binary64_detail::fromU32(a)) / b; } +GPUdi() constexpr GPUdoubleBinary64 operator+(float a, GPUdoubleBinary64 b) { return GPUdoubleBinary64(a) + b; } +GPUdi() constexpr GPUdoubleBinary64 operator-(float a, GPUdoubleBinary64 b) { return GPUdoubleBinary64(a) - b; } +GPUdi() constexpr GPUdoubleBinary64 operator*(float a, GPUdoubleBinary64 b) { return GPUdoubleBinary64(a) * b; } +GPUdi() constexpr GPUdoubleBinary64 operator/(float a, GPUdoubleBinary64 b) { return GPUdoubleBinary64(a) / b; } +#ifndef __METAL__ +GPUdi() constexpr GPUdoubleBinary64 operator+(double a, GPUdoubleBinary64 b) { return GPUdoubleBinary64((float)a) + b; } +GPUdi() constexpr GPUdoubleBinary64 operator-(double a, GPUdoubleBinary64 b) { return GPUdoubleBinary64((float)a) - b; } +GPUdi() constexpr GPUdoubleBinary64 operator*(double a, GPUdoubleBinary64 b) { return GPUdoubleBinary64((float)a) * b; } +GPUdi() constexpr GPUdoubleBinary64 operator/(double a, GPUdoubleBinary64 b) { return GPUdoubleBinary64((float)a) / b; } +#endif + +// rounds once, as `someFloat += someDouble` does on the host +#ifdef __METAL__ +GPUdi() thread float& operator+=(thread float& a, GPUdoubleBinary64 b) { return a = (float)(GPUdoubleBinary64(a) + b); } +GPUdi() device float& operator+=(device float& a, GPUdoubleBinary64 b) { return a = (float)(GPUdoubleBinary64(a) + b); } +GPUdi() threadgroup float& operator+=(threadgroup float& a, GPUdoubleBinary64 b) { return a = (float)(GPUdoubleBinary64(a) + b); } +#else +GPUdi() float& operator+=(float& a, GPUdoubleBinary64 b) { return a = (float)(GPUdoubleBinary64(a) + b); } +#endif + +namespace binary64_detail +{ +// sin and cos are fdlibm's __kernel_sin / __kernel_cos and the medium-range +// branch of __ieee754_rem_pio2. The coefficients are spelled as bit patterns +// because MSL has no double literals. The argument reduction is exact for +// |x| <= 2^19 * pi/2; beyond that the accuracy degrades gracefully. +GPUCA_B64_OP GPUdoubleBinary64 kernelSin(GPUdoubleBinary64 x, GPUdoubleBinary64 y, bool iy) +{ + typedef GPUdoubleBinary64 b64; + const b64 S1 = b64::fromBits(0xBFC5555555555549ULL); + const b64 S2 = b64::fromBits(0x3F8111111110F8A6ULL); + const b64 S3 = b64::fromBits(0xBF2A01A019C161D5ULL); + const b64 S4 = b64::fromBits(0x3EC71DE357B1FE7DULL); + const b64 S5 = b64::fromBits(0xBE5AE5E68A2B9CEBULL); + const b64 S6 = b64::fromBits(0x3DE5D93A5ACFD57CULL); + const b64 z = x * x; + const b64 v = z * x; + const b64 r = S2 + z * (S3 + z * (S4 + z * (S5 + z * S6))); + if (!iy) { + return x + v * (S1 + z * r); + } + return x - ((z * (b64::fromBits(0x3FE0000000000000ULL) * y - v * r) - y) - v * S1); +} + +GPUCA_B64_OP GPUdoubleBinary64 kernelCos(GPUdoubleBinary64 x, GPUdoubleBinary64 y) +{ + typedef GPUdoubleBinary64 b64; + const b64 C1 = b64::fromBits(0x3FA555555555554CULL); + const b64 C2 = b64::fromBits(0xBF56C16C16C15177ULL); + const b64 C3 = b64::fromBits(0x3EFA01A019CB1590ULL); + const b64 C4 = b64::fromBits(0xBE927E4F809C52ADULL); + const b64 C5 = b64::fromBits(0x3E21EE9EBDB4B1C4ULL); + const b64 C6 = b64::fromBits(0xBDA8FAE9BE8838D4ULL); + const b64 one = b64::fromBits(0x3FF0000000000000ULL); + const b64 oneHalf = b64::fromBits(0x3FE0000000000000ULL); + const u32 ix = (u32)(x.bits() >> 32) & 0x7fffffffu; + const b64 z = x * x; + const b64 r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * C6))))); + if (ix < 0x3FD33333u) { // |x| < 0.3 + return one - (oneHalf * z - (z * r - x * y)); + } + const b64 qx = (ix > 0x3FE90000u) ? b64::fromBits(0x3FD2000000000000ULL) : b64::fromBits((u64)(ix - 0x00200000u) << 32); // 0.28125, else |x| / 4 + return (one - qx) - ((oneHalf * z - qx) - (z * r - x * y)); +} + +GPUCA_B64_ALWAYS int32_t truncToInt32(GPUdoubleBinary64 x) +{ + const u64 b = x.bits(); + const int32_t e = (int32_t)((b >> 52) & 0x7ff) - 1023; + if (e < 0) { + return 0; + } + const u64 m = (b & GPUCA_B64_FRAC) | GPUCA_B64_IMPL; + const int32_t v = (int32_t)(e >= 52 ? (m << (e - 52)) : (m >> (52 - e))); + return (b & GPUCA_B64_SIGN) ? -v : v; +} + +struct SinCosPair { + GPUdoubleBinary64 s, c; +}; + +GPUCA_B64_OP SinCosPair sincos(GPUdoubleBinary64 x) +{ + typedef GPUdoubleBinary64 b64; + const u32 ix = (u32)(x.bits() >> 32) & 0x7fffffffu; + SinCosPair out; + if (ix >= 0x7FF00000u) { // inf or NaN + out.s = out.c = b64::fromBits(GPUCA_B64_DNAN); + return out; + } + + if (ix < 0x3E400000u) { // |x| < 2^-27, where the sign of a zero x has to survive + out.s = x; + out.c = b64::fromBits(0x3FF0000000000000ULL); + return out; + } + + b64 y0 = x, y1 = b64::fromBits(0ULL); + int32_t n = 0; + const bool reduced = ix > 0x3FE921FBu; // |x| > pi/4 + if (reduced) { + const b64 t = b64::fromBits(x.bits() & ~GPUCA_B64_SIGN); + n = truncToInt32(t * b64::fromBits(0x3FE45F306DC9C883ULL) + b64::fromBits(0x3FE0000000000000ULL)); + const b64 fn = b64((float)n); + // Cody-Waite with pi/2 split over three terms, good to 151 bits + b64 r = t - fn * b64::fromBits(0x3FF921FB54400000ULL); + b64 s = r; + b64 w = fn * b64::fromBits(0x3DD0B4611A600000ULL); + r = s - w; + w = fn * b64::fromBits(0x3BA3198A2E037073ULL) - ((s - r) - w); + s = r; + w = fn * b64::fromBits(0x3BA3198A2E000000ULL); + r = s - w; + w = fn * b64::fromBits(0x397B839A252049C1ULL) - ((s - r) - w); + y0 = r - w; + y1 = (r - y0) - w; + if (x.bits() & GPUCA_B64_SIGN) { + y0 = -y0; + y1 = -y1; + n = -n; + } + } + + switch (n & 3) { + case 0: + out.s = kernelSin(y0, y1, reduced); + out.c = kernelCos(y0, y1); + break; + case 1: + out.s = kernelCos(y0, y1); + out.c = -kernelSin(y0, y1, reduced); + break; + case 2: + out.s = -kernelSin(y0, y1, reduced); + out.c = -kernelCos(y0, y1); + break; + default: + out.s = -kernelCos(y0, y1); + out.c = kernelSin(y0, y1, reduced); + break; + } + return out; +} +} // namespace binary64_detail + +} // namespace o2::gpu + +#endif // GPUCOMMONDOUBLEBINARY64_H diff --git a/GPU/Common/GPUCommonMath.h b/GPU/Common/GPUCommonMath.h index 7a78a5881dcfa..7cbcbb3151b4d 100644 --- a/GPU/Common/GPUCommonMath.h +++ b/GPU/Common/GPUCommonMath.h @@ -308,6 +308,7 @@ GPUhdi() void GPUCommonMath::SinCos(float x, float& s, float& c) ) // clang-format on } +#ifndef __METAL__ GPUhdi() void GPUCommonMath::SinCosd(double x, double& s, double& c) { #if !defined(GPUCA_GPUCODE_DEVICE) && defined(__APPLE__) @@ -318,6 +319,7 @@ GPUhdi() void GPUCommonMath::SinCosd(double x, double& s, double& c) GPUCA_CHOICE((void)((s = sin(x)) + (c = cos(x))), sincos(x, &s, &c), s = sincos(x, &c)); #endif } +#endif GPUdi() constexpr uint32_t GPUCommonMath::Clz(uint32_t x) { @@ -444,11 +446,13 @@ GPUhdi() constexpr float GPUCommonMath::Abs(float x) return GPUCA_CHOICE(fabsf(x), fabsf(x), fabs(x)); } +#ifndef __METAL__ template <> GPUhdi() constexpr double GPUCommonMath::Abs(double x) { return GPUCA_CHOICE(fabs(x), fabs(x), fabs(x)); } +#endif template <> GPUhdi() constexpr int32_t GPUCommonMath::Abs(int32_t x) From d5a26f44a0af65d92e3e8f9394d1e97a105e482f Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Tue, 29 Sep 2026 13:24:08 +0000 Subject: [PATCH 2/2] Please consider the following formatting changes --- GPU/Common/GPUCommonDouble.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPU/Common/GPUCommonDouble.h b/GPU/Common/GPUCommonDouble.h index 2f7108b409a7e..88232d0b89915 100644 --- a/GPU/Common/GPUCommonDouble.h +++ b/GPU/Common/GPUCommonDouble.h @@ -34,7 +34,7 @@ static_assert(alignof(double) == 8, "the emulated double must match the alignmen template <> GPUhdi() constexpr double GPUCommonMath::Abs(double x) { - return double::fromBits(x.bits() & ~GPUCA_B64_SIGN); + return double ::fromBits(x.bits() & ~GPUCA_B64_SIGN); } // metal::fabs is not constant-evaluable, so this also fails to compile if the