Issue 16026 - std.math.frexp!float() wrong for very small subnormal values
Summary: std.math.frexp!float() wrong for very small subnormal values
Status: RESOLVED FIXED
Alias: None
Product: D
Classification: Unclassified
Component: phobos (show other issues)
Version: D2
Hardware: x86_64 Linux
: P1 normal
Assignee: thomas.bockman
URL:
Keywords:
Depends on:
Blocks:
 
Reported: 2016-05-15 10:04 UTC by thomas.bockman
Modified: 2016-05-22 21:07 UTC (History)
1 user (show)

See Also:


Attachments

Note You need to log in before you can comment on or make changes to this issue.
Description thomas.bockman 2016-05-15 10:04:23 UTC
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 Iain Buclaw 2016-05-15 12:21:41 UTC
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 thomas.bockman 2016-05-18 00:06:53 UTC
> 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 Iain Buclaw 2016-05-18 09:31:33 UTC
Maybe you're bit twiddling the float differently then.  I never said it was wrong, just that it wasn't the compiler. :-)
Comment 4 thomas.bockman 2016-05-18 10:05:44 UTC
> 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 thomas.bockman 2016-05-18 10:50:34 UTC
I decided to go ahead and start by submitting the quick fix, since it takes so long to get anything substantial accepted:
    https://github.com/dlang/phobos/pull/4337
Comment 6 github-bugzilla 2016-05-22 21:07:39 UTC
Commits pushed to master at https://github.com/dlang/phobos

https://github.com/dlang/phobos/commit/a6a1957be301860b1fa5d9205f1edd7a39cc0a1a
Fix issue 16026: std.math.frexp!float() wrong for very small subnormal values

https://github.com/dlang/phobos/commit/4991b82304b421f56e5394ba849826ee3251c89f
Merge pull request #4337 from tsbockman/issue_16026

Fix issue 16026: std.math.frexp!float() wrong for very small subnormals