BUG: Fix rounding of denormals in double and float to half casts … (#12722)

* BUG: Fix rounding of denormals in double and float to half casts

Previously the significand was shifted right to align denormals of
different magnitude. This loses some bits that can make a difference
for rounding. This is fixed:
 1. For floats, by inspecting the original last bits when this may
    make a difference (should happen rarely)
 2. For doubles by shifting the bits left to align the denromals and
    thus not lose the lowest orginal bits.

* TST: Test half denormal rounding

The test assumes that half to even is active, if this is ever
changed or allowed to change, the test will fail.

* Fixup: Fixup for halffloat.c

The one code path cannot be used. The other must have been a typo, but
is a valid bug, a new test for it in the next commit.

* TST: half casting lower bits are not lost for denormal results

The first test only tested the off by one, this one specifically tests
that all bits are used to decide if "round to nearest even" should
be used, in the example of rounding towards 0.

* TST: Test not just denormals but all positive finite float16s

* DOC: Add lots of comments and add a short release note.

1 files changed, 7 insertions, 0 deletions

diff --git a/doc/release/1.17.0-notes.rst b/doc/release/1.17.0-notes.rst
index 51e195437..d15936cb6 100644
--- a/doc/release/1.17.0-notes.rst
+++ b/doc/release/1.17.0-notes.rst
@@ -30,6 +30,13 @@ Expired deprecations
 Compatibility notes
 ===================
 
+float16 subnormal rounding
+--------------------------
+Casting from a different floating point precision to float16 used incorrect
+rounding in some edge cases. This means in rare cases, subnormal results will
+now be rounded up instead of down, changing the last bit (ULP) of the result.
+
+
 
 C API changes
 =============