Bit syntax float16 construction corrupts subnormals and misrounds normals in the software fallback (all Windows builds)

[a1x-an]

Describe the bug

Constructing 16-bit floats with bit syntax (<<X:16/float>>) produces wrong bit
patterns for nearly all values that are subnormal in binary16, and misrounds some
normal values, on builds that use the software conversion fallback in
erts/emulator/beam/erl_bits.c (i.e. when _Float16 is not available - notably
all official Windows/MSVC builds). Matching is correct, so the VM cannot
re-construct a value it has just successfully matched:

1> <<F:16/float>> = <<0,51>>.   %% 51 × 2⁻²⁴, exactly representable
2> <<F:16/float>>.
<<0,39>>                        %% construct(match(B)) ≠ B

To Reproduce

Erlang/OTP 29.0.3 (erts-17.0.3), Windows x86_64, official installer. Reproduces
identically in compiled code (JIT) and in the shell / erl_eval:

<<3.0517578125e-5:16/float>>.   %% = 2⁻¹⁵, exactly representable in binary16
%% expected <<2,0>> (0x0200), got <<0,1>> (0x0001) - off by a factor of 512

input (f64)	expected binary16	actual (Windows)
51 × 2⁻²⁴ = 3.039836883544922e-6	0x0033	0x0027
5 × 2⁻²⁴	0x0005	0x0003
2⁻¹⁵ = 3.0517578125e-5	0x0200	0x0001
2⁻²⁴ (min subnormal)	0x0001	0x0001 (correct by accident)
1023 × 2⁻²⁴ (max subnormal)	0x03FF	0x03FF (correct by accident)
1 + 4097 × 2⁻²³	0x3C01	0x3C00 (normal-range misrounding)

Every row except the last is exactly representable in binary16, so no rounding
mode can explain the results; the last row is strictly above the rounding
midpoint and must round up.

For scale: sweeping all finite float32 inputs against a compiler-provided
_Float16 oracle (gcc 13), the current fallback disagrees on 417,333,236
inputs, and fails construct∘match identity for 2026 of the 63,488 finite
binary16 values (all subnormals of both signs except twenty accidental
fixed points).

Expected behavior

IEEE 754 round-to-nearest-even, identical to the native _Float16 path, with
construct(match(B)) == B for every finite binary16 B.

Affected versions

Defects 1 and 2 (numbered as in Root cause under Additional context below)
are a regression introduced between OTP-27.0 and OTP-28.0. Up to and
including OTP 27 the fallback was a port of the FP16 library
(erl_bits_f16.h, https://github.com/Maratyszcza/FP16), whose subnormal
handling is correct; between OTP-27.0 and OTP-28.0 it was replaced by the
hand-written f32_to_f16(), which introduced the corruption. OTP 28.0
through 29.0.3 and current maint/master are affected. (The existing
bs_construct_SUITE:fp16/1 case only exercises the smallest and largest
subnormal - coincidentally the two values the broken code still gets right,
which is why the regression went unnoticed.) Defect 3 (double rounding
through float) dates back to the original implementation (OTP 24, PR #2890),
whose library API took a float.

Measured on shipped Windows installers (all of them fallback-path builds),
the latest patch of every line:

Shipped Windows installer	f16 segments	Defect 1: subnormals	Defect 2: sticky bits	construct∘match identity	Defect 3: double rounding
23.3.4.20	not supported (compiler warning "bad float bit size", badarg at run time)	-	-	-	-
24.3.4.17	yes	correct (51 × 2⁻²⁴ → 0x0033)	correct (→ 0x3C01)	0 failures	present (→ 0x3C00)
27.3.4.14	yes	correct	correct	0 failures	present
28.5.0.3 (latest 28 patch)	yes	broken (→ 39, 3, 1)	broken (→ 0x3C00)	2026 failures	present
29.0.3	yes	broken	broken	2026 failures	present

The 28 line still receives patches and remains affected in its latest
release, so the fix may be worth backporting there as well.

Only builds without a usable _Float16 are affected - which includes all
official Windows installers (MSVC has no _Float16), and any other
toolchain where the configure probe fails.

Additional context

Root cause (erts/emulator/beam/erl_bits.c, f32_to_f16()):

1. Subnormal branch. mantissa has already been reduced to the top 10 bits
   of the f32 mantissa ((u32 >> (23 - 10)) & 0x3ff), but the branch then ORs
   the implicit leading bit at position 23 (mantissa |= 1 << 23) - the position
   that belongs to the full 23-bit mantissa - and shifts by an exponent-derived
   amount. Hand-trace for 51×2⁻²⁴ (f32 biased exp 108, top-10 mantissa 608):
   significand = 0x800000+608 = 8389216, shift = 4 →
   (8389216 + 16) >> 4 = 524327; & 0x3ff = 39 - exactly the observed value.
   The same trace reproduces 2⁻¹⁵ → 1 and 5×2⁻²⁴ → 3, and shows min/max
   subnormals coming out right only by coincidence. The rounding increment
   round = 1 << shift is also a full ULP where half an ULP is required.

2. Normal branch. The round-to-nearest-even sticky test only inspects bits
   11 and 10 of the f32 word and ignores bits 9..0. A value with the guard bit
   set, sticky bits only in the low 10 bits, and an even kept LSB is treated as
   a tie and rounded down: <<1.00048840045928955078125:16/float>>
   (= 1 + 4097×2⁻²³) yields 0x3C00 instead of 0x3C01.

3. Double rounding. The fallback converts through float first
   (FP16_FROM_FP64(x) = f32_to_f16((float) x)), i.e. double→float→binary16
   with two roundings, while the native path casts double→_Float16 with one.
   <<1.000488282181322574615478515625:16/float>> (= 1 + 2⁻¹¹ + 2⁻³⁰) yields
   0x3C00 through the fallback but 0x3C01 on a native build - verified
   empirically on both sides (Windows 29.0.3 vs Linux/Alpine OTP 27,
   erts-15.2.7.8), so identical code gives platform-dependent results.

The native path (SIZEOF__FLOAT16 == 2 && FLOAT16_IS_CONVERTIBLE) is
unaffected: on a Linux build all of the values above construct correctly and
construct∘match holds for every finite binary16. The match direction
(f16_to_f32) is correct on all platforms.

I have a fix ready - replacing the fallback with a direct integer-exact
double→binary16 conversion (a single rounding, which also removes the platform
divergence) plus test cases for bs_construct_SUITE - validated exhaustively
against a compiler-provided _Float16 oracle: all finite float32 inputs,
identity over all finite binary16 values, and 10⁹ random finite doubles
produce identical results, with zero mismatches. PR to follow.