sycl: fix f32 to q2_0/q4_0/q5_0 rounding in CPY (#31120)

ggml/src/ggml-metal/kernels/quantize.h

@@ -20,6 +20,7 @@ void quantize_q1_0(device const float * src, device block_q1_0 & dst) {
 }
 
 void quantize_q2_0(device const float * src, device block_q2_0 & dst) {
+#pragma METAL fp math_mode(safe)
     float amax = 0.0f;
     for (int j = 0; j < QK2_0; j++) {
         float a = fabs(src[j]);

ggml/src/ggml-sycl/cpy.hpp

@@ -70,8 +70,22 @@ inline void cpy_blck_f32_q1_0(const char * cxi, char * cdsti) {
     }
 }
 
-inline int round_nearest_int(float x) {
-    return (int)(x >= 0.0f ? x + 0.5f : x - 0.5f);
+// exact roundf (half away from zero); x + 0.5f carries for x = 0.5 - 2^-25, and sycl::round is not exact under fast fp
+inline int round_half_away_int(float x) {
+    const float t = sycl::trunc(x);
+    return (int) (sycl::fabs(x - t) >= 0.5f ? t + sycl::copysign(1.0f, x) : t);
+}
+
+// correctly rounded a / b to match the CPU: fast fp gives an approximate reciprocal (~12% off by an ulp on B580),
+// so pick the approximation or a neighbour by smallest exact FMA residual
+inline float div_correctly_rounded(float a, float b) {
+    const float y  = a / b;
+    const float lo = sycl::nextafter(y, -INFINITY);
+    const float hi = sycl::nextafter(y, INFINITY);
+    const float e  = sycl::fabs(sycl::fma(-b, y, a));
+    const float el = sycl::fabs(sycl::fma(-b, lo, a));
+    const float eh = sycl::fabs(sycl::fma(-b, hi, a));
+    return el < e && el <= eh ? lo : (eh < e ? hi : y);
 }
 
 inline void cpy_blck_f32_q2_0(const char * cxi, char * cdsti) {
@@ -85,7 +99,7 @@ inline void cpy_blck_f32_q2_0(const char * cxi, char * cdsti) {
     }
 
     const float d  = amax;
-    const float id = d > 0.0f ? 1.0f / d : 0.0f;
+    const float id = d > 0.0f ? div_correctly_rounded(1.0f, d) : 0.0f;
 
     dsti->d = d;
 
@@ -94,7 +108,7 @@ inline void cpy_blck_f32_q2_0(const char * cxi, char * cdsti) {
     }
 
     for (int j = 0; j < QK2_0; ++j) {
-        int q = round_nearest_int(xi[j] * id) + 1;
+        int q = round_half_away_int(xi[j] * id) + 1;
         q = dpct::max(0, dpct::min(3, q));
 
         const int byte_index = j / 4;
@@ -378,7 +392,7 @@ inline void cpy_blck_f32_q4_0(const char * cxi, char * cdsti) {
     }
 
     const float d  = vmax / -8;
-    const float id = d ? 1.0f / d : 0.0f;
+    const float id = d ? div_correctly_rounded(1.0f, d) : 0.0f;
 
     dsti->d = d;
 
@@ -386,8 +400,9 @@ inline void cpy_blck_f32_q4_0(const char * cxi, char * cdsti) {
         const float x0 = xi[0 + j] * id;
         const float x1 = xi[QK4_0 / 2 + j] * id;
 
-        const uint8_t xi0 = dpct::min(15, (int8_t) (x0 + 8.5f));
-        const uint8_t xi1 = dpct::min(15, (int8_t) (x1 + 8.5f));
+        // no-op fmin stops x * id + 8.5f being fused into an FMA, so x * id is rounded as on the CPU
+        const uint8_t xi0 = dpct::min(15, (int8_t) (sycl::fmin(x0, 6.5f) + 8.5f));
+        const uint8_t xi1 = dpct::min(15, (int8_t) (sycl::fmin(x1, 6.5f) + 8.5f));
 
         dsti->qs[j] = xi0;
         dsti->qs[j] |= xi1 << 4;
@@ -442,7 +457,7 @@ inline void cpy_blck_f32_q5_0(const char * cxi, char * cdsti) {
     }
 
     const float d  = vmax / -16;
-    const float id = d ? 1.0f / d : 0.0f;
+    const float id = d ? div_correctly_rounded(1.0f, d) : 0.0f;
 
     dsti->d = d;
 
@@ -451,8 +466,9 @@ inline void cpy_blck_f32_q5_0(const char * cxi, char * cdsti) {
         const float x0 = xi[0 + j] * id;
         const float x1 = xi[QK5_0 / 2 + j] * id;
 
-        const uint8_t xi0 = dpct::min(31, (int8_t) (x0 + 16.5f));
-        const uint8_t xi1 = dpct::min(31, (int8_t) (x1 + 16.5f));
+        // fmin: no FMA, see cpy_blck_f32_q4_0
+        const uint8_t xi0 = dpct::min(31, (int8_t) (sycl::fmin(x0, 14.5f) + 16.5f));
+        const uint8_t xi1 = dpct::min(31, (int8_t) (sycl::fmin(x1, 14.5f) + 16.5f));
 
         dsti->qs[j] = (xi0 & 0xf) | ((xi1 & 0xf) << 4);
         qh |= ((xi0 & 0x10u) >> 4) << (j + 0);