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

New issue

Open

opened 2023-11-29 18:46:26 -06:00 by dpgeorge · 2 comments

dpgeorge commented

2023-11-29 18:46:26 -06:00

(Migrated from github.com)

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.

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: ```c #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 commented

2023-11-30 03:49:05 -06:00

(Migrated from github.com)

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);
}

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 commented

2023-12-03 19:04:27 -06:00

(Migrated from github.com)

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

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