Background: trying to solve a floating-point accuracy issue in MicroPython running on RP2040: micropython/micropython#13094

There seems to be an accuracy issue with the computation done by powf for certain arguments. In particular we see issues with powf(10, -11).

Code to reproduce:

#include <math.h>

static void pr(float f_in) {
    volatile float f = f_in;
    const uint8_t *b = (void *)&f;
    printf("%02x:%02x:%02x:%02x\n", b[0], b[1], b[2], b[3]);
}

int main(int argc, char **argv) {
    stdio_init_all();
    {
        int exp = -11;
        pr(powf(10.0, exp));
        pr(powf(10.0, -exp));
    }
    {
        volatile int exp = -11;
        pr(powf(10.0, exp));
        pr(powf(10.0, -exp));
    }

    return 0;
}

The second set of prints shows the issue when volatile int exp is used, which forces the evaluation to go through the full powf function (otherwise the compiler may optimise it out).

The results I get on a standard Pico with pico-sdk 1.5.1:

ff:eb:2f:2d
b7:43:ba:51
03:ec:2f:2d
b7:43:ba:51

The first two lines are correct (actually the second one is out by an LSB, it should be b8:43:ba:51). The third line is quite inaccurate, out by 4 LSBs.

The bug seems to be in the implementation of fpowint_1() and fpow_int2(), it looks like the latter is accumulating an error in its recursion.