Bug 5613 – std.mathspecial.betaIncomplete makes excessively stringent assumptions about FP Precision
Status
RESOLVED
Resolution
WORKSFORME
Severity
normal
Priority
P2
Component
phobos
Product
D
Version
D2
Platform
Other
OS
Windows
Creation time
2011-02-19T08:00:00Z
Last change time
2011-12-28T06:44:20Z
Assigned to
nobody
Creator
dsimcha
Comments
Comment #0 by dsimcha — 2011-02-19T08:00:49Z
std.mathspecial.betaIncomplete fails horribly anytime floating point precision is slightly degraded. For example, in GDC, exp2() is implemented (don't ask me why) in terms of powl(). This isn't terribly accurate, but the performance of std.mathspecial should degrade more gracefully in the face of this.
The relevant bug reports filed against GDC are:
https://bitbucket.org/goshawk/gdc/issue/140/stdmathspecialbetaincomplete-broken-for-64https://bitbucket.org/goshawk/gdc/issue/153/exp-not-computed-to-full-80-bit-precision
A test case that completely blows up if precision is somewhat reduced is:
import std.mathspecial, std.stdio;
void main() {
writeln(betaIncomplete(950, 51, 0.96));
}
From tracing through the code, it appears that betaIncomplete takes the gamma() instead of the logGamma() branch anytime the result can fit in an 80-bit real, including in my test case. This is silly. Instead it should leave some margin for safety and make sure the high order bits are right at the expense of some fuzz in the low order bits.
Comment #1 by clugdbug — 2011-02-19T11:54:02Z
If exp is only 64 bits, there is ZERO value in supporting 80-bit reals. None at all. If tests are only failing because of low exp() precision, they deserve to fail. I would say that DMD's implementation of 80-bit exp and transcendentals is as inaccurate as you could possibly allow. Probably too inaccurate. If anything does worse than it, it's really, really bad.
OTOH we will need 64-bit function results, for machines where real == double.
Comment #2 by dsimcha — 2011-02-19T12:13:06Z
Ok, so I guess the fix is to put in a version statement to redefine all those constants that determine what branch is used, so that it's still correct in 64-bit precision. Then, to get GDC working until its exp() function stops sucking, maybe the GDC patches should set the flag to assume 64-bit precision on GDC x64, since as you note, for all practical purposes it is only 64-bit.
Comment #3 by dsimcha — 2011-12-28T06:44:20Z
I'm closing this because it now works fine on GDC, LDC and 64-bit DMD, so it wasn't as much of a problem as I thought it was.