Bug 16026 – std.math.frexp!float() wrong for very small subnormal values
Status
RESOLVED
Resolution
FIXED
Severity
normal
Priority
P1
Component
phobos
Product
D
Version
D2
Platform
x86_64
OS
Linux
Creation time
2016-05-15T10:04:00Z
Last change time
2016-05-22T21:07:39Z
Assigned to
thomas.bockman
Creator
thomas.bockman
Comments
Comment #0 by thomas.bockman — 2016-05-15T10:04:23Z
import std.math, std.traits, std.stdio;
T test(T)()
if (isFloatingPoint!T)
{
int exp;
return frexp(3 * (T.epsilon * T.min_normal), exp);
}
void main() {
writeln(test!real()); // 0.75
writeln(test!double()); // 0.5?? (But 0.75 on GDC.)
writeln(test!float()); // 0.5?? (But 0.75 on GDC.)
}
---
I tried tracking this down myself in the std.math code, but I couldn't find anything wrong. Maybe it's a codegen bug?
Comment #1 by ibuclaw — 2016-05-15T12:21:41Z
This is actually a result of turning frexp into a template.
Before it did the operation at real precision, then down casted. Now, it calls frexp with the precision of the type T.
It's not a codegen bug, I can reproduce the same using GDC and Phobos 2.067.
Comment #2 by thomas.bockman — 2016-05-18T00:06:53Z
> Before it did the operation at real precision, then down casted.
The function is bugged. I did a straightforward rewrite of frexp() just based on the docs, and my version produces the expected answer without any need for extra precision.
However, my rewrite depends upon a revised version of my floating-point decomposition union type, so fixing #16026 will have to wait until after the union itself gets accepted into std.math or std.bitmanip.
(And that, in turn, may have to wait until I get a newly discovered DMD codegen bug fixed.)
Comment #3 by ibuclaw — 2016-05-18T09:31:33Z
Maybe you're bit twiddling the float differently then. I never said it was wrong, just that it wasn't the compiler. :-)
Comment #4 by thomas.bockman — 2016-05-18T10:05:44Z
> Maybe you're bit twiddling the float differently then.
Ah! I found it:
// float subnormal
vf *= F.RECIP_EPSILON;
ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
exp = ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1;
vu[F.EXPPOS_SHORT] =
cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) | 0x3F00);
Above, the 0x8000 should really be ~F.EXPMASK, because the fractional part overlaps with the highest 16 bits, where the exponent and sign are. The double version has the same problem. I bet it got copy-pasta-ed from the 80-bit and 128-bit code, where 0x8000 == ~F.EXPMASK.
I could submit a quick two-line fix, but I'd rather actually clean up all those "magic numbers" and the copy-pasta: https://github.com/dlang/phobos/pull/4336
If PR #4336 gets accepted, I can just replace the whole frexp() implementation with something more maintainable. (I already have it working; but need to add ibmExtended support.)
Comment #5 by thomas.bockman — 2016-05-18T10:50:34Z