std.mathspecial: for 64-bit reals, provide MAXGAMMA, MAXLOG, and MINLOG #3045
Conversation
|
I've added you to the autotester whitelist: https://auto-tester.puremagic.com/pulls.ghtml?projectid=1 |
|
Ah, great @yebblies , thanks! |
| @@ -231,7 +238,12 @@ real igammaTemmeLarge(real a, real x) | |||
|
|
|||
| public: | |||
| /// The maximum value of x for which gamma(x) < real.infinity. | |||
| enum real MAXGAMMA = 1755.5483429L; | |||
| static if (real.mant_dig == 64) // 80-bit reals | |||
There was a problem hiding this comment.
Could you use floatTraits from std.math ?
See also #3011
There was a problem hiding this comment.
Indeed I tried, but floatTraits is still private. I agree it'd be better to use floatTraits.
When is #3011 going to be merged?
There was a problem hiding this comment.
I don't know about time, but DMD 2.068.
This PR can be good argument for #3011.
You may want to add dependency note in PR description: somebody from core team will add dependency label.
There was a problem hiding this comment.
I think now I screwed up this branch/PR while trying to update it to upstream master :'(
Edit: OK fixed it.
JohanEngelen
force-pushed
the
gamma64bit
branch
from
7c62dfc
to
2ac2581
Compare
|
Updated the PR using floatTraits from std.math. Depends on PR #3011. |
|
The DMC build for DMD uses 80 bit reals, even while supporting linking with MSVC libraries. Frankly, it's a bug in the MSVC generated build of DMD that 80 bit reals are not supported. That doesn't affect this PR, though. |
|
This PR requires PR #3011 |
|
PR #3011 was merged. Please remove trailing white spaces and rebase. |
JohanEngelen
force-pushed
the
gamma64bit
branch
2 times, most recently
from
f6596f4
to
00fda8c
Compare
| @@ -111,7 +111,14 @@ real gammaStirling(real x) | |||
| } | |||
| else { | |||
| w = 1.0L + w * poly( w, SmallStirlingCoeffs); | |||
| y = pow( x, x - 0.5L ) / y; | |||
| if ( (floatTraits!(real).realFormat == RealFormat.ieeeDouble) && (x > 143.0L) ) { | |||
There was a problem hiding this comment.
Would the first condition be optimised at compile time?
There was a problem hiding this comment.
You think it is better to do a static if on the first condition?
I think I chose against it, because of having to duplicate the y = pow( x, x - 0.5L ) / y;, but I have no preference. The alternative could be:
else {
w = 1.0L + w * poly( w, SmallStirlingCoeffs);
static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble) {
// Avoid overflow in pow() for 64-bit reals
if (x > 143.0L) ) {
real v = pow( x, 0.5L * x - 0.25L );
y = v * (v / y);
}
else {
y = pow( x, x - 0.5L ) / y;
}
}
else {
y = pow( x, x - 0.5L ) / y;
}
}
There was a problem hiding this comment.
Yes. This code looks better for math. Furthermore we would have this functions as template in the future, so please replace 0.5L for 64bit float with 0.5.
Inside gammaStirling(), fix a call to pow() that would overflow for 64-bit reals.
|
@9il I've made some additional fixes and rebased again. Thank for the reviews. At this point I have no further motivation to work on this (certainly not motivated to port CEPHES's |
| } | ||
| else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble) { | ||
| enum real MAXLOG = 0x1.62e42fefa39efp+9L; // log(real.max) | ||
| enum real MINLOG = -0x1.74385446d71c3p+9L; // log(real.min_normal*real.epsilon) = log(smallest denormal) |
There was a problem hiding this comment.
Please replace with this lines. The MINLOG was changed, the MAXLOG remains the same.
enum real MAXLOG = 0x1.62e42fefa39efp+9; /// log(2^+1024) at Wolfram Alpha, 21.07.2015
enum real MINLOG = -0x1.74910d52d3052p+9; /// log(2^-1075) at Wolfram Alpha, 21.07.2015There was a problem hiding this comment.
@9il Are you sure your MINLOG is correct? I believe the smallest denormal is 2^-1074, according to Wikipedia. In that case, I think my value (the one calculated by DMD) is correct and is the same as calculated by Wolfram Alpha and then converted to hex (do you know how to have Wolfram Alpha do that conversion to 64-bit float hex?).
There was a problem hiding this comment.
Smallest denormal is 2^-1074, however CEPHES uses 2^-1075 (probably because rounding rule?). Lets it be 2^-1074 ;)
CEPHES github fork fo reference (origin has the same constant).
Wolfram:
- Example
- D code example for transformation
writefln("%a", 0x2e9.221aa5a60a3p+0);
There was a problem hiding this comment.
Never believe what you read on Wikipedia. ;)
Just, kidding, this is fine.
|
LGTM |
|
Looks fine - anything else to review @ibuclaw? |
|
Auto-merge toggled on |
std.mathspecial: for 64-bit reals, provide MAXGAMMA, MAXLOG, and MINLOG
Because MAXGAMMA, MAXLOG, and MINLOG are only defined for 80-bit reals, gammafunction.d unittests fail for targets that have 64-bit reals (notably MSVC).
This PR conditionally defines the 64-bit equivalent constants. Also, it contains a fix for gammaStirling(x) in which a call to pow() would overflow for 64-bit reals for x > 143. This fixes several failing unittests. The changes should be a noop for 80-bit reals.
floatTraits from std.math is used to detect whether 80- or 64-bit reals are used; therefore there is a dependency on PR #3011.
The changes are tested using DMD 2.066.1 (Windows) and LDC2 (Windows, version 192583 based on DMD v2.066.1 and LLVM 3.7.0svn-r230699).
Comments are very welcome.