Summary

Current implementation of REPR_C works by clearing the two lower bits of the mantissa to zero. As it happens after each floating point operation, this tends to bias floating point numbers towards zero, causing decimals like .9997 instead of rounded numbers. This is visible in test cases involving repeated computations, such as misc/rge_sm for instance, or when looking at the numerical results of perf_bench/bm_fft, which shows significantly higher errors than when using full 32 bit floats.

The standard way to fix this bias toward zero would be to add 2 to the mantissa before clearing the last two bits, effectively performing a round operation. This approach works quite well to fix the bias, but the repeated rounding operation still has some impact on accuracy on the longer term.

By experimenting with alternative approches, I found out that the most effective approach is to fix the bias not when saving the floating point number into an mp_obj_t, but when reloading the float value from the mp_obj_t, simply by copying the the two previous bits. Although less mathematically sound, this method appears empirically to work significantly better, although I don't have a mathematical proof to offer. The only explanation I have for this good behaviour is that statistically, it adds the required half-bit bias, while properly restoring the missing decimals for numbers with a mantissa ending in 0000 and 1111.

Here is a summary of the comparison between current baseline, and these two methods:

Given this improvement in REPR_C accuracy, and as discussed in PR 16953, I have updated the unix/longlong variant to use REPR_C, to ensure that it is included in CI builds. A few existing test cases had to be tweaked as the updated code did cause in a few cases the last digit to change, but this has been made explicit in the test code.

The only possibly problematic case found during testing of this code was found in the implementation of bm_fft, which computes log_2(n) as int(log(n) / log(2)). This expression is lacking a call to round() to make it safe, as nothing guarantees that the result of the division is not just a fraction below the expected integer value. This is actually what was happening when running the new code, causing a test failure. The fix to the test case was easy, but we cannot exclude that other pieces of code in the wild have been using the same expression. They would fail in the same way, unless the call to roundis added.

Testing

The new code has been tested as described above, on the latest development builds.

Trade-offs and Alternatives

Two methods to fix the bias have been described, and the most effective has been chosen.

Another approach would be to use an alternate REPR which favors floats over smallints, but discussions in MicroPython forum showed that this was not desirable.