Re: [avr-libc-dev] Pow Function in avr8

[Georg-Johann Lay]

Amine Najahi schrieb:

> The floating point standard leaves their accuracy implementation
> defined (for now, this will surely change in the future) which means
> that as soon as you use one of these functions, your program has
> probably no chance to be reproducible on a different machine.

Well, this will be the case if epsilon is greater than zero.

See for example the "most often reported non-bug in GCC".

http://gcc.gnu.org/bugs/#nonbugs

> In your case, all you can say, is that on your examples, your host
> machine has a better implementation that avr-libc. A small error
> occurs while log is computed, gets amplified once multiplied by the
> exponent i, which causes a bigger error while computing exp (knowing
> that exp itself introduces its proper error too).
>
> To see the error appear since log's computation, take an int pointer
> on the float value, and print it as an hexadecimal integer.

Computing 2^5, once by AVR Libc's pow and once by GCC, will yield:

1: pow (2.000000, 5.000000) = 32.000000 = 42000000, eps = 0.000000e+00
2: pow (2.000000, 5.000000) = 31.999989 = 41fffffa, eps = 3.576279e-07
3: pow (2.000000, 5) = 32.000000 = 42000000

To verify that the second call is folded by the compiler, skim the 
generated assembler.

This means the deviation from the correct result 32 is 3 LSBs,
i.e. 3 * 2^{-23} under the assumption that eps computation is not lossy.

Notice avr-gcc has built-in macro

#define __DBL_EPSILON__ ((double)1.19209290e-7L)


You may also take into account that the first derivative of 2^x at
5 is ~22, i.e. an epsilon in the input transforms to a 22-fold epsilon
of the output provided x is around 5.

Johann

--

Reference code used with avr-gcc 4.7.2 + AVR Libc 1.8.0:

#include <stdio.h>
#include <math.h>

static unsigned int long tox (double x)
{
    long l;
    asm ("" : "=r" (l) : "0" (x));
    return l;
}

static double eps (double x, double x0)
{
    return (x0 - x) / x0;
}

void testit (void)
{
    double f1, f2, f3, x1 = 5.0, x2 = 5.0, b = 2.0;

    asm ("" : "+r" (x2));

    f1 = pow (b, x1);
    f2 = pow (b, x2);
    f3 = pow (b, 5);
    printf ("1: pow (%f, %f) = %f = %lx, eps = %e\n",
            b, x1, f1, tox (f1), eps (f1, 32));
    printf ("2: pow (%f, %f) = %f = %lx, eps = %e\n",
            b, x2, f2, tox (f2), eps (f2, 32));
    printf ("3: pow (%f, 5) = %f = %lx\n", b, f3, tox (f3));
}