| Summary: | __builtin_modfl returns incorrect value | ||
|---|---|---|---|
| Product: | gcc | Reporter: | Matthias Kretz (Vir) <mkretz> |
| Component: | target | Assignee: | Not yet assigned to anyone <unassigned> |
| Status: | RESOLVED MOVED | ||
| Severity: | normal | CC: | amodra, fw, jakub |
| Priority: | P3 | ||
| Version: | 10.1.1 | ||
| Target Milestone: | --- | ||
| URL: | https://sourceware.org/bugzilla/show_bug.cgi?id=18031 | ||
| See Also: | https://sourceware.org/bugzilla/show_bug.cgi?id=18031 | ||
| Host: | Target: | powerpc64le-*-* | |
| Build: | Known to work: | ||
| Known to fail: | Last reconfirmed: | ||
The testcase just assumes that modfl will be evaluated the same at compile time and at runtime, but for the broken by design double double format that is not the case. GCC internally treats the type as one with fixed 106 bit precision and that is why you get those results, while at runtime it actually is a variable bit precision (53 to 2000-ish mantissa bits) depending on the exact value. The runtime modf actually returns a large number. This is not about precision but about completely bogus values. You can adjust the testcase to:
int e = 69;
int main() {
long double a = -__builtin_ldexpl(
1.9446689187403240306919491832695730985733566864714824565497322973045558e+00l,
e);
long double tmp;
long double aa = __builtin_modfl(a, &tmp);
if (aa >= 1) {
__builtin_abort();
}
}
I should be more precise. Take this test case:
int e = 69;
int main() {
__ibm128 a = -__builtin_ldexpl(
1.9446689187403240306919491832695730985733566864714824565497322973045558e+00l, e);
__ibm128 tmp;
__ibm128 aa = __builtin_modfl(a, &tmp);
if (aa >= 1) {
__builtin_printf("modfl: %Lf\n", aa);
__builtin_printf("copysign(x-trunc(x), x): %Lf\n", __builtin_copysignl(a - __builtin_truncl(a), a));
return 1;
}
return 0;
}
it prints:
modfl: 43227.412695
copysign(x-trunc(x), x): -0.587305
Printing the numbers in detail (for all of them the long double and the two double portions after that)): a -0x1.f1d5d27f89914ab924b2cd099500000000p+69 -0x1.f1d5d27f89915000p+69 0x1.51b6d34cbd9ac000p+15 b -0x1.f1d5d27f89914ab924b2cd099500000000p+69 -0x1.f1d5d27f89915000p+69 0x1.51b6d34cbd9ac000p+15 aa 0x1.51b6d34cbd9ac000000000000000000000p+15 0x1.51b6d34cbd9ac000p+15 0x0.0000000000000000p+0 bb -0x1.2cb3426540000000000000000000000000p-1 -0x1.2cb3426540000000p-1 0x0.0000000000000000p+0 tmp1 -0x1.f1d5d27f89915000000000000000000000p+69 -0x1.f1d5d27f89915000p+69 0x0.5000000000000000p-1022 tmp2 -0x1.f1d5d27f89914ab9200000000000000000p+69 -0x1.f1d5d27f89915000p+69 0x1.51b8000000000000p+15 But e.g. __builtin_truncl works the same on both, I'll have a look in a debugger. Actually, I think this is a glibc bug rather than gcc, because the values computed at compile time look right, rather than those at runtime.
Consider:
int e = 69;
int main() {
long double a = -__builtin_ldexpl(1.9446689187403240306919491832695730985733566864714824565497322973045558e+00l, e);
long double b = -__builtin_ldexpl(1.9446689187403240306919491832695730985733566864714824565497322973045558e+00l, 69);
long double tmp1, tmp2;
long double aa = __builtin_modfl(a, &tmp1);
long double bb = __builtin_modfl(b, &tmp2);
long double c = __builtin_truncl(a);
long double d = __builtin_truncl(b);
__builtin_printf ("a %.34La %.16a %.16a\n", a, (double) a, (double) (a - (double) a));
__builtin_printf ("b %.34La %.16a %.16a\n", b, (double) b, (double) (b - (double) b));
__builtin_printf ("aa %.34La %.16a %.16a\n", aa, (double) aa, (double) (aa - (double) aa));
__builtin_printf ("bb %.34La %.16a %.16a\n", bb, (double) bb, (double) (bb - (double) bb));
__builtin_printf ("tmp1 %.34La %.16a %.16a\n", tmp1, (double) tmp1, (double) (tmp1 - (double) tmp1));
__builtin_printf ("tmp2 %.34La %.16a %.16a\n", tmp2, (double) tmp2, (double) (tmp2 - (double) tmp2));
__builtin_printf ("c %.34La %.16a %.16a\n", c, (double) c, (double) (c - (double) c));
__builtin_printf ("d %.34La %.16a %.16a\n", d, (double) d, (double) (d - (double) d));
return 0;
}
The tmp1, tmp2, c and d values should be the same, i.e. the integral part.
But only tmp2, c and d are the same, -0x1.f1d5d27f89914ab9200000000000000000p+69
while tmp1 is -0x1.f1d5d27f89914ab9200000000000000000p+69.
c is truncl computed at runtime, and tmp1 is the integral part from modfl computed at runtime, so they should be the same, but they are not.
Seems glibc for modfl in this case just uses the upper double mantissa bits and ignores the lower double mantissa bits.
Please file against glibc.
truncl has:
hi = trunc (xh);
if (hi != xh)
{
/* The high part is not an integer; the low part does not
affect the result. */
xh = hi;
xl = 0;
}
else
{
/* The high part is a nonzero integer. */
lo = xh > 0 ? floor (xl) : ceil (xl);
xh = hi;
xl = lo;
ldbl_canonicalize_int (&xh, &xl);
}
i.e. if the double trunc of the highpart double is equal to the highpart double (this case), it does another floor or ceil.
While the modfl code assumes that if the exponent is > 52 and smaller than 103, the fractional part is the second double, which is wrong.
The (i1&0x7fffffffffffffffLL) in there is weird because there is
i1 &= 0x000fffffffffffffLL; earlier, so it could just test i1, but more importantly, the } else { ... code assumes that the low double must have a particular exponent (it could have any smaller than the higher exponent minus 50something).
Even the j0 > 103 case is wrong. First of all, it is unclear where that 103 comes from, e.g. for the quad version it uses j0 > 111, i.e. two smaller than the 113 bit precision, but 103 is 3 smaller than 106. But more importantly, it assumes that such values have no fractional part, which is not given. E.g. the upper double could be 0x1.0p+110 and lower 0x1.0f0f0f0f0f0fp+30. There is a glibc bug about this already for years, just hasn't been fixed yet: https://sourceware.org/bugzilla/show_bug.cgi?id=18031 |
Test case: int e = 69; int main() { long double a = -__builtin_ldexpl( 1.9446689187403240306919491832695730985733566864714824565497322973045558e+00l, e); long double b = -__builtin_ldexpl( 1.9446689187403240306919491832695730985733566864714824565497322973045558e+00l, 69); long double tmp; long double aa = __builtin_modfl(a, &tmp); long double bb = __builtin_modfl(b, &tmp); if (aa != bb) { if (a == b) return 2; else return 1; } return 0; } This returns 2, i.e. the internal representation of a and b differ in such a way that modf gets broken. GCC was built with `--target=powerpc64le-linux-gnu --disable-nls --without-headers --with-long-double-128`