www.digitalmars.com         C & C++   DMDScript  

c++ - Accuracy problem with powl

reply "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
parent reply "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
parent reply "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
parent reply "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
parent reply "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
parent reply "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
parent reply "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
more information).
Mar 22 2004
parent "Walter" <walter digitalmars.com> writes:
"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?
Mar 22 2004