On Tue, 14 Jun 2022, LIU Hao wrote:
在 2022/6/14 10:04, sisyphus 写道:
I think the different results can also be expressed as:
(1/3) ^ 5 produces 3f70db20a88f4695
1 / (3 ^ 5) produces 3f70db20a88f4696
The latter is a far less error prone approach, and therefore yields the
result that should be provided for (1 / 3) ^ 5.
It's not easy to tell which one is closer to the exact result, like this:
```
#include <stdio.h>
#include <quadmath.h>
volatile double x = 190000000000.6;
int
main(void)
{
char temp[64];
printf("x = %a\n", x);
printf("x^2 = %a\n", x * x);
printf("1/x = %a\n", 1 / x);
printf("1/(x^2) = %a\n", 1 / (x * x));
printf("(1/x)^2 = %a\n", (1 / x) * (1 / x));
quadmath_snprintf(temp, sizeof temp, "%Qa", powq(x, -2));
printf("exact = %s\n", temp);
}
```
gives (which is the same result as on Windows)
```
lh_mouse@lhmouse-xps ~/Desktop $ gcc test.c -Wall -Wextra -lquadmath
lh_mouse@lhmouse-xps ~/Desktop $ ./a.out
x = 0x1.61e70f6004ccdp+37
x^2 = 0x1.e93f08f38d71ep+74
1/x = 0x1.725c9fad04f66p-38
1/(x^2) = 0x1.0be7ef89a5262p-75
(1/x)^2 = 0x1.0be7ef89a5261p-75
exact = 0x1.0be7ef89a5261607a8abf8991a41p-75
```
It's actually the square of reciprocal that is more close to the quadruple
value, than the reciprocal of square.
In this case, both of those results are within 1 ULP from the desired
result, so it's all up to whichever rounding direction is better for this
particular value. But for the values involved in the original bug report,
the current algorithm produced an error more than a couple ULP, right?
// Martin
_______________________________________________
Mingw-w64-public mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public
