From 6530318cfa91f3e5bdc4d1027c81cfe67bdb3114 Mon Sep 17 00:00:00 2001 From: a1x-an <138564121+a1x-an@users.noreply.github.com> Date: Fri, 3 Jul 2026 14:36:37 +0300 Subject: [PATCH] erts: Fix float16 construction in the software fallback The software conversion used when _Float16 is not available (notably all official Windows/MSVC builds) produced wrong bit patterns for nearly all values that are subnormal in binary16, and misrounded some normal values: 1> <> = <<0,51>>. %% 51 * 2^-24, exact 2> <>. <<0,39>> %% construct(match(B)) =/= B The subnormal branch of f32_to_f16() OR:ed the implicit leading bit into position 23 after the mantissa had already been reduced to its top 10 bits, and used a whole ULP as the rounding increment. The normal branch ignored bits 9..0 in its sticky test, so values just above a rounding midpoint could round down as if they were ties. Additionally, converting through float (double -> float -> binary16) rounds twice, so results could differ from the native _Float16 path for the same input. Replace the fallback with a direct integer-exact double -> binary16 conversion with a single round-to-nearest-even step, giving the same results as the native path. The match direction (f16_to_f32) was correct and is unchanged. Verified exhaustively against a compiler-provided _Float16 oracle: all 2^32 float32 inputs, construct/match identity over all finite binary16 values, and 10^9 random doubles produce identical results. Also verified in a patched emulator built from maint on x86_64 Linux with the fallback both disabled (native path, unchanged) and forced via FLOAT16_IS_CONVERTIBLE (fixed path, now identical to native). --- erts/emulator/beam/erl_bits.c | 124 ++++++++++++---------- erts/emulator/test/bs_construct_SUITE.erl | 24 +++++ 2 files changed, 93 insertions(+), 55 deletions(-) diff --git a/erts/emulator/beam/erl_bits.c b/erts/emulator/beam/erl_bits.c index f23648af2c73..6b70a2e6625c 100644 --- a/erts/emulator/beam/erl_bits.c +++ b/erts/emulator/beam/erl_bits.c @@ -48,7 +48,7 @@ typedef _Float16 erlfp16; #define FP16_TO_FP64(x) ((double) x) #else typedef Uint16 erlfp16; -#define FP16_FROM_FP64(x) (f32_to_f16((float) x)) +#define FP16_FROM_FP64(x) (f64_to_f16(x)) #define FP16_TO_FP64(x) ((double) f16_to_f32(x)) #endif @@ -423,72 +423,86 @@ erts_bs_get_binary_2(Process *p, Uint num_bits, ErlSubBits *sb) return result; } -static ERTS_INLINE Uint16 f32_to_f16(float fp) +/* Software fallback for FP16_FROM_FP64: convert a double directly to a + * binary16 with a single round-to-nearest-even step, mirroring the + * semantics of the native (_Float16) cast used when _Float16 is + * available. Converting through float first would round twice and could + * give results that differ from the native path. */ +static ERTS_INLINE Uint16 f64_to_f16(double fp) { union { - float f32; - Uint32 u32; + double f64; + Uint64 u64; } u; - Uint32 u32; - Uint32 sign, exp, mantissa; - int signed_exp; - Uint16 res; + Uint64 man, keep, low; + Uint16 sign; + int exp, e, shift, guard, sticky; - u.f32 = fp; - u32 = u.u32; + u.f64 = fp; - sign = (u32 >> 31) & 0x1; - exp = (u32 >> 23) & 0xff; - mantissa = (u32 >> (23 - 10)) & 0x3ff; + sign = (Uint16) ((u.u64 >> 48) & 0x8000); + exp = (int) ((u.u64 >> 52) & 0x7ff); + man = u.u64 & ((((Uint64) 1) << 52) - 1); + + /* Erlang floats are always finite, small integers always fit in a + * finite double, and oversized bignums are rejected by + * big_to_double() before we are called, so exp == 0x7ff (Inf/NaN) + * cannot occur here. If it ever did, it would fall through to the + * overflow branch below and yield infinity. */ + ASSERT(exp != 0x7ff); if (exp == 0) { - /* Convert zero to f16. */ - res = sign << 15; - return res; - } - - signed_exp = exp - (127 - 15); - exp -= (127 - 15); /* Convert exponent from f32 bias to f16 bias. */ - if (signed_exp <= 0) { - Uint32 shift; - mantissa |= 1 << 23; - shift = -signed_exp; - if (shift <= 24 ) { - /* Subnormal value in f16. */ - Uint32 round = 1 << shift; - mantissa = (mantissa + round + ((mantissa >> shift) & 1)) >> shift; - res = sign << 15 | (mantissa & 0x3ff); - } else { - /* Underflow to 0. */ - res = sign << 15; + /* Subnormal doubles are far below the binary16 subnormal range; + * they round to zero. */ + return sign; + } + + e = exp - 1023; + man |= ((Uint64) 1) << 52; /* Full 53-bit significand. */ + + if (e >= -14) { + /* Normal value in binary16: keep 11 significand bits and round + * the 42 bits below them. */ + Uint64 k; + if (e > 15) { + /* Overflow becomes infinity. */ + return (Uint16) (sign | 0x7C00); } - } else if (exp > 0x1f) { - /* Overflow becomes infinity. */ - res = (sign << 15) | (0x1f << 10); - return res; - } else { - /* Normal value in f16. Apply rounding. */ - Uint32 bit_11 = u32 & (0x1 << 10); - Uint32 bit_12 = u32 & (0x1 << 11); - Uint32 bit_13 = u32 & (0x1 << 12); - Uint32 bit_14 = u32 & (0x1 << 13); - if (bit_13 && (bit_11 || bit_12 || bit_14)) { + k = man >> 42; /* 11 bits: 1024..2047 */ + guard = (int) ((man >> 41) & 1); + sticky = (man & ((((Uint64) 1) << 41) - 1)) != 0; + if (guard && (sticky || (k & 1))) { /* Round to nearest, ties to even. */ - mantissa += 1; - if (mantissa == 1 << 10) { - /* Mantissa overflow: carry into exponent. */ - exp += 1; - mantissa = 0; - if (exp > 0x1f) { - /* Overflow becomes infinity. */ - res = (sign << 15) | (0x1f << 10); - return res; + k++; + if (k == 2048) { + /* Significand overflow: carry into the exponent. */ + k = 1024; + e++; + if (e > 15) { + return (Uint16) (sign | 0x7C00); } } } - res = sign << 15 | exp << 10 | (mantissa & 0x3ff); - } - return res; + return (Uint16) (sign | ((Uint16) (e + 15) << 10) | (Uint16) (k & 0x3ff)); + } + + /* Subnormal value in binary16: the stored significand is + * RNE(|fp| * 2^24). A rounded result of 1024 naturally lands on + * 0x0400, the smallest normal number. */ + shift = 28 - e; /* e <= -15, so shift >= 43. */ + if (shift > 63) { + /* Underflow to zero. */ + return sign; + } + keep = man >> shift; + guard = (int) ((man >> (shift - 1)) & 1); + low = man & ((((Uint64) 1) << (shift - 1)) - 1); + sticky = low != 0; + if (guard && (sticky || (keep & 1))) { + /* Round to nearest, ties to even. */ + keep++; + } + return (Uint16) (sign | (Uint16) keep); } static ERTS_INLINE float f16_to_f32(Uint16 fp) diff --git a/erts/emulator/test/bs_construct_SUITE.erl b/erts/emulator/test/bs_construct_SUITE.erl index aceb28b800b8..9ca7e6c90cc9 100644 --- a/erts/emulator/test/bs_construct_SUITE.erl +++ b/erts/emulator/test/bs_construct_SUITE.erl @@ -1126,6 +1126,30 @@ fp16(_Config) -> %% others ?FP16(16#4000, 2), ?FP16(16#4000, 2.0), + + %% Subnormals used to be corrupted by the software conversion + %% fallback. All of these are exactly representable in binary16, + %% so no rounding is involved. + ?FP16(16#0033, 3.039836883544921875e-6), %% 51 * 2^-24 + ?FP16(16#0005, 2.98023223876953125e-7), %% 5 * 2^-24 + ?FP16(16#0200, 3.0517578125e-5), %% 2^-15 + ?FP16(16#0163, 2.1159648895263671875e-5), %% 355 * 2^-24 + %% Sticky bits below bit 10 used to be ignored when rounding + %% normal values; this is strictly above the rounding midpoint + %% and must round up. + ?FP16(16#3c01, 1.00048840045928955078125), %% 1 + 4097 * 2^-23 + %% An exact tie must round to even. + ?FP16(16#3c00, 1.00048828125), %% 1 + 2^-11 + %% double -> binary16 must round once; rounding through float + %% first lands on the float tie point and rounds to even + %% instead of up. + ?FP16(16#3c01, 1.000488282181322574615478515625), %% 1 + 2^-11 + 2^-30 + %% Constructing a matched binary16 must reproduce it exactly. + _ = [begin + <> = <>, + <> = <> + end || H <- lists:seq(0, 16#ffff), + (H bsr 10) band 16#1f =/= 16#1f], ok. zero_init(_Config) ->