BUG: Fix rounding of denormals in double and float to half casts … #12722
Conversation
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.
The test assumes that half to even is active, if this is ever changed or allowed to change, the test will fail.
numpy/core/src/npymath/halffloat.c
Outdated
| } | ||
| #else | ||
| d_sig += 0x0000020000000000ULL; | ||
| d_sig += 0x0000e00000000000ULL; |
There was a problem hiding this comment.
Is this right? I'm not sure I'm following the logic, but it seems you're only adding a single bit before, and now you're adding something with 3 bits set. It looks like you changed 2 to 2 * 70 = 0xe0, but shouldn't it be changed to 2 << 7 = 0x100 as above? I.e., 0x0001000000000000ULL
There was a problem hiding this comment.
wow, yes you are right :(, that is embaressing.
There was a problem hiding this comment.
Ah, that was some left-over in a code path that is never used (and thus cannot be tested), we could remove that completely in principle.
There was a problem hiding this comment.
Yes, I guess we always do proper round-to-even now.
numpy/core/tests/test_half.py
Outdated
| @@ -69,6 +69,58 @@ def test_half_conversions(self): | ||
| j = np.array(i_f16, dtype=int) | ||
| assert_equal(i_int, j) | ||
|
|
||
| @pytest.mark.parametrize("offset", [None, "up", "down"]) | ||
| @pytest.mark.parametrize("shift", [None, "up", "down"]) | ||
| @pytest.mark.parametrize("float_t", [np.float32,]) |
There was a problem hiding this comment.
Also test with np.float64? Otherwise why parametrize it?
There was a problem hiding this comment.
Yeah, forgot to read after trying, hopefully the float64s fail ;) because of that 3 bit thing above.
numpy/core/src/npymath/halffloat.c
Outdated
| */ | ||
| if ((f_sig&0x00003fffu) != 0x00001000u) { | ||
| if (((f_sig&0x00003fffu) != 0x00001000u) || (f&0x008ffu)) { |
There was a problem hiding this comment.
Where did 008ff come from? And why are you giving only 5 hex digits?
There was a problem hiding this comment.
Jeesh, I messed up, probably a typo, it should have been 7ff and of course full digits. 7ff for the last 11 bits, that would otherwise be gone.
There was a problem hiding this comment.
Dang, the tests probably do not catch this, because it is in the wrong bits... maybe should add a test for this part more specifically.
| /* Handle rounding by adding 1 to the bit beyond half precision */ | ||
| #if NPY_HALF_ROUND_TIES_TO_EVEN | ||
| /* | ||
| * If the last bit in the half significand is 0 (already even), and | ||
| * the remaining bit pattern is 1000...0, then we do not add one | ||
| * to the bit after the half significand. In all other cases, we do. | ||
| */ | ||
| if ((d_sig&0x000007ffffffffffULL) != 0x0000020000000000ULL) { | ||
| d_sig += 0x0000020000000000ULL; | ||
| if ((d_sig&0x003fffffffffffffULL) != 0x0010000000000000ULL) { |
There was a problem hiding this comment.
Might be more readable to express these as (1 << n) - 1 and (1 << n) for some n. Fine to leave as it too.
There was a problem hiding this comment.
Yeah, should change it everywhere then though probably.
There was a problem hiding this comment.
Given you got this wrong at least twice in this patch, I think it might be wise to put this patch on hold, and make this change first in a standalone MAINT PR. A simple npy_uint64 n_low_bits(int n) function would probably be helpful here.
There was a problem hiding this comment.
Probably a good idea, these bit patterns are a error prone :/, so I may do that but not super soon maybe. For now I had already fixed it up and added a test for the one actual problem.
There was a problem hiding this comment.
Without spaces is par for the course for this file - I don't think we need to block this PR on that.
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.
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.
|
Would having #10741 make testing easier here? |
|
I looked at the change made to the following code and it appears correct. |
| * significand. This changes the constant compared to normal branch. | ||
| */ | ||
| assert(d_exp - 998 >= 0); | ||
| d_sig <<= (d_exp - 998); |
There was a problem hiding this comment.
If you can shed any light in a comment on on where 998 comes from, that would be appreciated.
|
Actually, first had @eric-wieser I guess float to integer ratio is likely nice more generally. I am not sure how it would help here specifically (unless maybe to not rely on float16->float32/64 being correct). But we may want to reactivate that PR in any case, since python has long moved on and the maintenance burden seems not high? |
|
LGTM. Still needs release note |
|
Hmm, I thought/think release notes make sense, because it is a subtle change (we can still remove it anyway). But would do we want to backport it, and so that it should be already in 1.16.1 or both? |
|
Plus, in that case I should maybe look into giving it a nice backportable merge base. |
Needs to be decided for the LTS branch, we don't yet have a policy for that. I'm inclined to go for duplicate entries because the LTS branch is the end of Python 2.7 support, so a bit more separate than usual. |
|
Checking if reviewers all think this is ready. |
|
Sorry, got a bit side tracked, I don't quite remember myself if I thought it was. I think I was considering doing the simple fix in both versions still. I doubt it makes a noticable speed difference anyway. |
|
It seems to me this is ready - there was a sensible question about a handy helper macro, but probably best to put it in: at some level, the tests are what really counts here, and those now pass. |
|
OK, no more feedback, so putting this in. Thanks Sebastian. |
…umpy#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.
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:
make a difference (should happen rarely)
thus not lose the lowest orginal bits.
I am not sure if the double branch is not being a bit too smart and it is nicer to just also use the
||logic from above, since it uses different constants compared to the normal branch because of this.I do not know if there may be nicer ways with larger code changes, but I could not think of one.
This closes gh-12721