## c++ - Accuracy problem with powl

• Ronald Barrett (23/23) Mar 13 2004 From the following source:
• Walter (4/27) Mar 13 2004 The result could actually be the same - gcc may be rounding the result. ...
• Ronald Barrett (55/57) Mar 14 2004 #include
• Walter (4/6) Mar 18 2004 Ah, there's the problem. MinGW is dropping the last few bits of the resu...
"Ronald Barrett" <ronaldb sebcorinc.com> writes:
```From the following source:
#include <stdio.h>
#include <math.h>

int main(void){
long double a, b, x;

x = 0x1.6000022202b1076ap-58l; //4.77049e-18l;
a = powl(1.0l - x, 1.0l - x);
b = powl(x, 1.0l - powl(x, 1.0l - x));

printf("x = %.20Lf\n", x);
printf("a = %.20Lf\n", a);
printf("b = %.20Lf\n", b);

printf("result = %Lg\n", a + b - 1.0l);

return 0;
}

the output is:
x = 0.00000000000000000477
a = 0.99999999999999999514
b = 0.00000000000000000477
result = -5.42101e-20

Mathematically the result should be >= 0 for x in the range [0;1].
For example gcc v3.3.1 (mingw special 20030804-1) outputs result = 0 for
similar cases (including this).

Ronald
```
Mar 13 2004
"Walter" <walter digitalmars.com> writes:
```The result could actually be the same - gcc may be rounding the result. Try
printing the result out in binary and see.

"Ronald Barrett" <ronaldb sebcorinc.com> wrote in message
news:c30bsj\$1enq\$1 digitaldaemon.com...
From the following source:
#include <stdio.h>
#include <math.h>

int main(void){
long double a, b, x;

x = 0x1.6000022202b1076ap-58l; //4.77049e-18l;
a = powl(1.0l - x, 1.0l - x);
b = powl(x, 1.0l - powl(x, 1.0l - x));

printf("x = %.20Lf\n", x);
printf("a = %.20Lf\n", a);
printf("b = %.20Lf\n", b);

printf("result = %Lg\n", a + b - 1.0l);

return 0;
}

the output is:
x = 0.00000000000000000477
a = 0.99999999999999999514
b = 0.00000000000000000477
result = -5.42101e-20

Mathematically the result should be >= 0 for x in the range [0;1].
For example gcc v3.3.1 (mingw special 20030804-1) outputs result = 0 for
similar cases (including this).

Ronald

```
Mar 13 2004
"Ronald Barrett" <ronaldb sebcorinc.com> writes:
```#include <stdio.h>
#include <math.h>
#include <stdint.h>
#include <inttypes.h>

void printfp(long double value)
{
for (int counter = 0; counter < 10; counter++)
printf("%02" PRIX8, *(((uint8_t*) &value) + counter));

printf("\n\n");
}

int main(void){
long double a, b, x, result;

x = 0x1.6000022202b1076ap-58l; //4.77049e-18l;
a = powl(1.0l - x, 1.0l - x);
b = powl(x, 1.0l - powl(x, 1.0l - x));

printf("x = %.20Lf\n", x);
printfp(x);

printf("a = %.20Lf\n", a);
printfp(a);

printf("b = %.20Lf\n", b);
printfp(b);

result = a + b - 1.0l;

printf("result = %Lg\n", result);
printfp(result);

return 0;
}

DMC v8.40 output:
x = 0.00000000000000000477
B5835801110100B0C53F

a = 0.99999999999999999514
A7FFFFFFFFFFFFFFFE3F

b = 0.00000000000000000477
2A8D5801110100B0C53F

result = -5.42101e-20
0000000000000080BFBF

gcc v3.3.1 (mingw special 20030804-1) output:
x = -0.00000000000000000000
B5835801110100B0C53F

a = -1.#QNAN000000000000000
A8FFFFFFFFFFFFFFFE3F

b = -0.00000000000000000000
228D5801110100B0C53F

result = 0
00000000000000000000

The C runtime library of MinGW gcc is MSVC based and therefore printf
recognizes 80 bit long double as 64 bit one.

"Walter" <walter digitalmars.com> wrote in message
news:c312v6\$2jhk\$2 digitaldaemon.com...
The result could actually be the same - gcc may be rounding the result.

Try
printing the result out in binary and see.

There is a difference between the outputs:
for DMC: A7FFFFFFFFFFFFFFFE3F + 2A8D5801110100B0C53F - 1.0l =
0000000000000080BFBF < 0.0l
for gcc: A8FFFFFFFFFFFFFFFE3F + 228D5801110100B0C53F - 1.0l = 0.0l

Also gcc and DMC have the same default rounding direction (FE_TONEAREST) for
operations + and -.
```
Mar 14 2004
"Walter" <walter digitalmars.com> writes:
```"Ronald Barrett" <ronaldb sebcorinc.com> wrote in message
news:c31r78\$o5h\$1 digitaldaemon.com...
The C runtime library of MinGW gcc is MSVC based and therefore printf
recognizes 80 bit long double as 64 bit one.

Ah, there's the problem. MinGW is dropping the last few bits of the result,
which is a bug in MinGW's long double support.
```
Mar 18 2004
"Ronald Barrett" <ronaldb sebcorinc.com> writes:
```"Walter" <walter digitalmars.com> wrote in message
news:c3du2o\$4tn\$1 digitaldaemon.com...
"Ronald Barrett" <ronaldb sebcorinc.com> wrote in message
news:c31r78\$o5h\$1 digitaldaemon.com...
The C runtime library of MinGW gcc is MSVC based and therefore printf
recognizes 80 bit long double as 64 bit one.

Ah, there's the problem. MinGW is dropping the last few bits of the

result,
which is a bug in MinGW's long double support.

With Lcc-win32 the result is 0.0l (lcc haven't the MinGW gcc printf problem
because it uses different runtime library). Also with MinGW gcc (with the
Dinkum library) the result is also 0.0l: A8FFFFFFFFFFFFFFFE3F +
238D5801110100B0C53F - 1.0l = 0.0l.

The MinGW gcc printf problem is unrelated to the powl accuracy.
```
Mar 19 2004
"Walter" <walter digitalmars.com> writes:
```"Ronald Barrett" <ronaldb sebcorinc.com> wrote in message
news:c3fhfd\$2t3c\$1 digitaldaemon.com...
"Walter" <walter digitalmars.com> wrote in message
news:c3du2o\$4tn\$1 digitaldaemon.com...
"Ronald Barrett" <ronaldb sebcorinc.com> wrote in message
news:c31r78\$o5h\$1 digitaldaemon.com...
The C runtime library of MinGW gcc is MSVC based and therefore printf
recognizes 80 bit long double as 64 bit one.

Ah, there's the problem. MinGW is dropping the last few bits of the

result,
which is a bug in MinGW's long double support.

With Lcc-win32 the result is 0.0l (lcc haven't the MinGW gcc printf

problem
because it uses different runtime library). Also with MinGW gcc (with the
Dinkum library) the result is also 0.0l: A8FFFFFFFFFFFFFFFE3F +
238D5801110100B0C53F - 1.0l = 0.0l.

The MinGW gcc printf problem is unrelated to the powl accuracy.

Why is 'a' being set to a QNAN?
```
Mar 21 2004
"Ronald Barrett" <ronaldb sebcorinc.com> writes:
```"Walter" <walter digitalmars.com> wrote in message
news:c3jjgt\$irg\$1 digitaldaemon.com...
Why is 'a' being set to a QNAN?

#include <stdio.h>
#include <stdint.h>
#include <inttypes.h>

void printfp(long double value)
{
for (int counter = 0; counter < 10; counter++)
printf("%02" PRIX8, *(((uint8_t*) &value) + counter));

printf("\n");
}

int main(void){
long double a = 0x1.ffffffffffffff5p-1l;
printfp(a);
printf("a = %.20lf\n", *((double*) &a));

return 0;
}

dmc test.c & test
A8FFFFFFFFFFFFFFFE3F
a = -nan

As I said the printf problem of the standard MinGW gcc C library is
unrelated to the powl accuracy (you can ask the MinGW gcc developers for
```"Ronald Barrett" <ronaldb sebcorinc.com> wrote in message