File 2321-erts-Replace-FP16-implementation-with-a-new-one.patch of Package erlang

From f7c3e3614b05de73f57ee19a6ccfd53456ce2bc5 Mon Sep 17 00:00:00 2001
From: lucioleKi <isabell@erlang.org>
Date: Thu, 20 Mar 2025 10:36:14 +0100
Subject: [PATCH] erts: Replace FP16 implementation with a new one

The new implementation has the same behavior as the previous one.
The reason we are writing this new implementation is that the
previous version is external and has an MIT license. If it gets a newer
version, we need to upgrade it and keep it up to date. Because FP16 has
a clear specification and is stand-alone, it is easier to just write
our own implementation.

The newer compilers already have native support for FP16, so this
implementation is only relevant for older compilers. For this reason,
the new implementation has not been tested for speed.
---
 .github/scripts/otp-compliance.es          |   2 +-
 erts/emulator/beam/erl_bits.c              | 114 +++++++++++-
 erts/emulator/beam/erl_bits_f16.h          | 205 ---------------------
 erts/emulator/beam/vendor.info             |  13 --
 erts/emulator/test/bs_construct_SUITE.erl  |   2 +
 erts/emulator/test/bs_match_misc_SUITE.erl |   4 +
 6 files changed, 118 insertions(+), 222 deletions(-)
 delete mode 100644 erts/emulator/beam/erl_bits_f16.h

diff --git a/erts/emulator/beam/erl_bits.c b/erts/emulator/beam/erl_bits.c
index a1843715dd..f23648af2c 100644
--- a/erts/emulator/beam/erl_bits.c
+++ b/erts/emulator/beam/erl_bits.c
@@ -48,9 +48,8 @@ typedef _Float16 erlfp16;
 #define FP16_TO_FP64(x) ((double) x)
 #else
 typedef Uint16 erlfp16;
-#define FP16_FROM_FP64(x) (fp16_ieee_from_fp32_value((float) x))
-#define FP16_TO_FP64(x) ((double) fp16_ieee_to_fp32_value(x))
-#include "erl_bits_f16.h"
+#define FP16_FROM_FP64(x) (f32_to_f16((float) x))
+#define FP16_TO_FP64(x) ((double) f16_to_f32(x))
 #endif
 
 /*
@@ -424,6 +423,115 @@ erts_bs_get_binary_2(Process *p, Uint num_bits, ErlSubBits *sb)
     return result;
 }
 
+static ERTS_INLINE Uint16 f32_to_f16(float fp)
+{
+    union {
+        float f32;
+        Uint32 u32;
+    } u;
+    Uint32 u32;
+    Uint32 sign, exp, mantissa;
+    int signed_exp;
+    Uint16 res;
+
+    u.f32 = fp;
+    u32 = u.u32;
+
+    sign = (u32 >> 31) & 0x1;
+    exp = (u32 >> 23) & 0xff;
+    mantissa = (u32 >> (23 - 10)) & 0x3ff;
+
+    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;
+        }
+    } 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)) {
+            /* 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;
+                }
+            }
+        }
+        res = sign << 15 | exp << 10 | (mantissa & 0x3ff);
+    }
+    return res;
+}
+
+static ERTS_INLINE float f16_to_f32(Uint16 fp)
+{
+    union {
+        float f32;
+        Uint32 u32;
+    }u;
+    Uint32 sign, exp, mantissa;
+    Uint32 res;
+
+    sign = (fp >> 15) & 0x1;
+    exp = (fp >> 10) & 0x1f;
+    mantissa = fp & 0x3ff;
+    res = sign << 31;
+    if (exp == 0) {
+        if (mantissa != 0) {
+            /* Normalize subnormals in f16 to f32 */
+            while ((mantissa & 0x400) == 0) {
+                mantissa <<= 1;
+                exp -= 1;
+            }
+            mantissa = mantissa & 0x3ff;
+            exp = 127 - 15 + 1 + exp;
+            res |= exp << 23;
+            res |= mantissa << 13;
+        }
+    } else if (exp == 0x1f) {
+        /* Convert infinity and NaN to f32 */
+        res |= 0xff << 23;
+        if (mantissa != 0) {
+            res |= mantissa << 13;
+        }
+    } else {
+        /* Convert normal values to f32. */
+        exp += (127 - 15);
+        res |= exp << 23;
+        res |= mantissa << 13;
+    }
+    u.u32 = res;
+    return u.f32;
+}
+
 Eterm
 erts_bs_get_float_2(Process *p, Uint num_bits, unsigned flags, ErlSubBits *sb)
 {
diff --git a/erts/emulator/beam/erl_bits_f16.h b/erts/emulator/beam/erl_bits_f16.h
deleted file mode 100644
index cda52ab719..0000000000
--- a/erts/emulator/beam/erl_bits_f16.h
+++ /dev/null
@@ -1,205 +0,0 @@
-/*
- * C-implementation of Float16 - https://github.com/Maratyszcza/FP16
- *
- * Version: 0a92994d729ff76a58f692d3028ca1b64b145d91
- *
- * The MIT License (MIT)
- *
- * Copyright (c) 2017 Facebook Inc.
- * Copyright (c) 2017 Georgia Institute of Technology
- * Copyright 2019 Google LLC
- *
- * Permission is hereby granted, free of charge, to any person obtaining a copy of this
- * software and associated documentation files (the "Software"), to deal in the Software
- * without restriction, including without limitation the rights to use, copy, modify,
- * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to
- * permit persons to whom the Software is furnished to do so, subject to the following
- * conditions:
- *
- * The above copyright notice and this permission notice shall be included in all copies
- * or substantial portions of the Software.
- *
- * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
- * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
- * PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE
- * FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
- * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
- * DEALINGS IN THE SOFTWARE.
- */
-
-/*
- * The following changes were done to adapt this to Erlang/OTP conventions:
- *
- * 1. uint16_t and uint32_t have been rewritten to Uint16 and Uint32
- * 2. Mixed declarations were moved to the top to avoid warnings
- * 3. inline was rewritten as ERTS_INLINE
- * 4. UINT16_C(x) and UINT32_C(x) were rewritten to xU as we don't support 16bits platform
- */
-
-static ERTS_INLINE float fp32_from_bits(Uint32 w) {
-    union {
-	Uint32 as_bits;
-	float as_value;
-    } fp32 = { w };
-    return fp32.as_value;
-}
-
-static ERTS_INLINE Uint32 fp32_to_bits(float f) {
-    union {
-	float as_value;
-	Uint32 as_bits;
-    } fp32 = { f };
-    return fp32.as_bits;
-}
-
-/*
- * Convert a 32-bit floating-point number in IEEE single-precision format to a 16-bit floating-point number in
- * IEEE half-precision format, in bit representation.
- *
- * @note The implementation relies on IEEE-like (no assumption about rounding mode and no operations on denormals)
- * floating-point operations and bitcasts between integer and floating-point variables.
- */
-static ERTS_INLINE Uint16 fp16_ieee_from_fp32_value(float f) {
-#if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L) || defined(__GNUC__) && !defined(__STRICT_ANSI__)
-    const float scale_to_inf = 0x1.0p+112f;
-    const float scale_to_zero = 0x1.0p-110f;
-#else
-    const float scale_to_inf = fp32_from_bits(0x77800000U);
-    const float scale_to_zero = fp32_from_bits(0x08800000U);
-#endif
-    float base = (fabsf(f) * scale_to_inf) * scale_to_zero;
-
-    const Uint32 w = fp32_to_bits(f);
-    const Uint32 shl1_w = w + w;
-    const Uint32 sign = w & 0x80000000U;
-
-    Uint32 bits, exp_bits, mantissa_bits, nonsign, bias;
-
-    bias = shl1_w & 0xFF000000U;
-
-    if (bias < 0x71000000U) {
-	bias = 0x71000000U;
-    }
-
-    base = fp32_from_bits((bias >> 1) + 0x07800000U) + base;
-    bits = fp32_to_bits(base);
-    exp_bits = (bits >> 13) & 0x00007C00U;
-    mantissa_bits = bits & 0x00000FFFU;
-    nonsign = exp_bits + mantissa_bits;
-    return (sign >> 16) | (shl1_w > 0xFF000000U ? 0x7E00U : nonsign);
-}
-
-/*
- * Convert a 16-bit floating-point number in IEEE half-precision format, in bit representation, to
- * a 32-bit floating-point number in IEEE single-precision format.
- *
- * @note The implementation relies on IEEE-like (no assumption about rounding mode and no operations on denormals)
- * floating-point operations and bitcasts between integer and floating-point variables.
- */
-static ERTS_INLINE float fp16_ieee_to_fp32_value(Uint16 h) {
-    /*
-     * Extend the half-precision floating-point number to 32 bits and shift to the upper part of the 32-bit word:
-     *      +---+-----+------------+-------------------+
-     *      | S |EEEEE|MM MMMM MMMM|0000 0000 0000 0000|
-     *      +---+-----+------------+-------------------+
-     * Bits  31  26-30    16-25            0-15
-     *
-     * S - sign bit, E - bits of the biased exponent, M - bits of the mantissa, 0 - zero bits.
-     */
-    const Uint32 w = (Uint32) h << 16;
-    /*
-     * Extract the sign of the input number into the high bit of the 32-bit word:
-     *
-     *      +---+----------------------------------+
-     *      | S |0000000 00000000 00000000 00000000|
-     *      +---+----------------------------------+
-     * Bits  31                 0-31
-     */
-    const Uint32 sign = w & 0x80000000U;
-    /*
-     * Extract mantissa and biased exponent of the input number into the high bits of the 32-bit word:
-     *
-     *      +-----+------------+---------------------+
-     *      |EEEEE|MM MMMM MMMM|0 0000 0000 0000 0000|
-     *      +-----+------------+---------------------+
-     * Bits  27-31    17-26            0-16
-     */
-    const Uint32 two_w = w + w;
-
-    /*
-     * Shift mantissa and exponent into bits 23-28 and bits 13-22 so they become mantissa and exponent
-     * of a single-precision floating-point number:
-     *
-     *       S|Exponent |          Mantissa
-     *      +-+---+-----+------------+----------------+
-     *      |0|000|EEEEE|MM MMMM MMMM|0 0000 0000 0000|
-     *      +-+---+-----+------------+----------------+
-     * Bits   | 23-31   |           0-22
-     *
-     * Next, there are some adjustments to the exponent:
-     * - The exponent needs to be corrected by the difference in exponent bias between single-precision and half-precision
-     *   formats (0x7F - 0xF = 0x70)
-     * - Inf and NaN values in the inputs should become Inf and NaN values after conversion to the single-precision number.
-     *   Therefore, if the biased exponent of the half-precision input was 0x1F (max possible value), the biased exponent
-     *   of the single-precision output must be 0xFF (max possible value). We do this correction in two steps:
-     *   - First, we adjust the exponent by (0xFF - 0x1F) = 0xE0 (see exp_offset below) rather than by 0x70 suggested
-     *     by the difference in the exponent bias (see above).
-     *   - Then we multiply the single-precision result of exponent adjustment by 2**(-112) to reverse the effect of
-     *     exponent adjustment by 0xE0 less the necessary exponent adjustment by 0x70 due to difference in exponent bias.
-     *     The floating-point multiplication hardware would ensure than Inf and NaN would retain their value on at least
-     *     partially IEEE754-compliant implementations.
-     *
-     * Note that the above operations do not handle denormal inputs (where biased exponent == 0). However, they also do not
-     * operate on denormal inputs, and do not produce denormal results.
-     */
-    const Uint32 exp_offset = 0xE0U << 23;
-#if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L) || defined(__GNUC__) && !defined(__STRICT_ANSI__)
-    const float exp_scale = 0x1.0p-112f;
-#else
-    const float exp_scale = fp32_from_bits(0x7800000U);
-#endif
-    const float normalized_value = fp32_from_bits((two_w >> 4) + exp_offset) * exp_scale;
-
-    /*
-     * Convert denormalized half-precision inputs into single-precision results (always normalized).
-     * Zero inputs are also handled here.
-     *
-     * In a denormalized number the biased exponent is zero, and mantissa has on-zero bits.
-     * First, we shift mantissa into bits 0-9 of the 32-bit word.
-     *
-     *                  zeros           |  mantissa
-     *      +---------------------------+------------+
-     *      |0000 0000 0000 0000 0000 00|MM MMMM MMMM|
-     *      +---------------------------+------------+
-     * Bits             10-31                0-9
-     *
-     * Now, remember that denormalized half-precision numbers are represented as:
-     *    FP16 = mantissa * 2**(-24).
-     * The trick is to construct a normalized single-precision number with the same mantissa and thehalf-precision input
-     * and with an exponent which would scale the corresponding mantissa bits to 2**(-24).
-     * A normalized single-precision floating-point number is represented as:
-     *    FP32 = (1 + mantissa * 2**(-23)) * 2**(exponent - 127)
-     * Therefore, when the biased exponent is 126, a unit change in the mantissa of the input denormalized half-precision
-     * number causes a change of the constructud single-precision number by 2**(-24), i.e. the same amount.
-     *
-     * The last step is to adjust the bias of the constructed single-precision number. When the input half-precision number
-     * is zero, the constructed single-precision number has the value of
-     *    FP32 = 1 * 2**(126 - 127) = 2**(-1) = 0.5
-     * Therefore, we need to subtract 0.5 from the constructed single-precision number to get the numerical equivalent of
-     * the input half-precision number.
-     */
-    const Uint32 magic_mask = 126U << 23;
-    const float magic_bias = 0.5f;
-    const float denormalized_value = fp32_from_bits((two_w >> 17) | magic_mask) - magic_bias;
-
-    /*
-     * - Choose either results of conversion of input as a normalized number, or as a denormalized number, depending on the
-     *   input exponent. The variable two_w contains input exponent in bits 27-31, therefore if its smaller than 2**27, the
-     *   input is either a denormal number, or zero.
-     * - Combine the result of conversion of exponent and mantissa with the sign of the input number.
-     */
-    const Uint32 denormalized_cutoff = 1U << 27;
-    const Uint32 result = sign |
-	(two_w < denormalized_cutoff ? fp32_to_bits(denormalized_value) : fp32_to_bits(normalized_value));
-    return fp32_from_bits(result);
-}
diff --git a/erts/emulator/test/bs_construct_SUITE.erl b/erts/emulator/test/bs_construct_SUITE.erl
index 9baad051d6..aceb28b800 100644
--- a/erts/emulator/test/bs_construct_SUITE.erl
+++ b/erts/emulator/test/bs_construct_SUITE.erl
@@ -1121,6 +1121,8 @@ fp16(_Config) ->
     ?FP16(16#3c01, 1.00097656),
     %% rounding of 1/3 to nearest
     ?FP16(16#3555, 0.33325195),
+    %% infinity
+    ?FP16(16#7c00, 1000000000),
     %% others
     ?FP16(16#4000, 2),
     ?FP16(16#4000, 2.0),
diff --git a/erts/emulator/test/bs_match_misc_SUITE.erl b/erts/emulator/test/bs_match_misc_SUITE.erl
index aaab728a4a..48ada29b12 100644
--- a/erts/emulator/test/bs_match_misc_SUITE.erl
+++ b/erts/emulator/test/bs_match_misc_SUITE.erl
@@ -82,6 +82,10 @@ t_float(Config) when is_list(Config) ->
 
     {'EXIT',{{badmatch,_},_}} = (catch match_float(<<0,0>>, 8, 0)),
     {'EXIT',{{badmatch,_},_}} = (catch match_float(<<0,0>>, 16#7fffffff, 0)),
+    NaN = <<0:1,31:5,42:10>>,
+    {'EXIT',{{badmatch,_},_}} = (catch match_float(NaN, 16, 0)),
+    Inf = <<0:1,31:5,0:10>>,
+    {'EXIT',{{badmatch,_},_}} = (catch match_float(Inf, 16, 0)),
 
     ok.
 
-- 
2.43.0

openSUSE Build Service is sponsored by