`powf(10, -11)` is not accurate

[dpgeorge]

Background: trying to solve a floating-point accuracy issue in MicroPython running on RP2040: https://github.com/micropython/micropython/issues/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.

[peterharperuk]

This is what I get running on ubuntu

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

Spraying a few doubles around gets you this on Pico...

static double fpow_int2(double x,int y) {
    double u;
    if(y==1) return x;
    u=fpow_int2(x,y/2);
    u*=u;
    if(y&1) u*=x;
    return u;
}

// for the case where x not zero or infinity, y small and not zero
static inline float fpowint_1(double x,int y) {
    if(y<0) x=1.0f/x,y=-y;
    return fpow_int2(x,y);
}

[dpgeorge]

See the fix suggested here: https://github.com/micropython/micropython/issues/13094#issuecomment-1837652290 That computes the reciprocal first, then the power.