Thread overview
Accuracy problem with powl
Mar 14, 2004
Ronald Barrett
Mar 14, 2004
Walter
Mar 14, 2004
Ronald Barrett
Mar 19, 2004
Walter
Mar 19, 2004
Ronald Barrett
Mar 21, 2004
Walter
Mar 22, 2004
Ronald Barrett
Mar 22, 2004
Walter
March 14, 2004
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


March 14, 2004
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
>
>


March 14, 2004
#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 -.


March 19, 2004
"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.


March 19, 2004
"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.


March 21, 2004
"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?


March 22, 2004
"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 more information).


March 22, 2004
"Ronald Barrett" <ronaldb@sebcorinc.com> wrote in message news:c3n2c1$23v$1@digitaldaemon.com...
> 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 more information).

Would you like to take a look at \dm\src\core\pow.c for the powl()
implementation?