Skip to content

std.mathspecial: for 64-bit reals, provide MAXGAMMA, MAXLOG, and MINLOG #3045

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ibuclaw merged 1 commit into from
Jul 30, 2015
Merged

std.mathspecial: for 64-bit reals, provide MAXGAMMA, MAXLOG, and MINLOG #3045

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 30 additions & 4 deletions std/internal/math/gammafunction.d
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,7 +111,19 @@ real gammaStirling(real x)
}
else {
w = 1.0L + w * poly( w, SmallStirlingCoeffs);
y = pow( x, x - 0.5L ) / y;
static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble) {
// Avoid overflow in pow() for 64-bit reals
if (x > 143.0) {
real v = pow( x, 0.5 * x - 0.25 );
y = v * (v / y);
}
else {
y = pow( x, x - 0.5 ) / y;
}
}
else {
y = pow( x, x - 0.5L ) / y;
}
}
y = SQRT2PI * y * w;
return y;
Expand Down Expand Up @@ -231,7 +243,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 (floatTraits!(real).realFormat == RealFormat.ieeeExtended)
enum real MAXGAMMA = 1755.5483429L;
else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
enum real MAXGAMMA = 171.6243769L;
else
static assert(0, "missing MAXGAMMA for other real types");


/*****************************************************
Expand Down Expand Up @@ -531,8 +548,17 @@ unittest {


private {
enum real MAXLOG = 0x1.62e42fefa39ef358p+13L; // log(real.max)
enum real MINLOG = -0x1.6436716d5406e6d8p+13L; // log(real.min*real.epsilon) = log(smallest denormal)
static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended) {
enum real MAXLOG = 0x1.62e42fefa39ef358p+13L; // log(real.max)
enum real MINLOG = -0x1.6436716d5406e6d8p+13L; // log(real.min_normal*real.epsilon) = log(smallest denormal)
}
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)
Copy link
Member

@9il 9il Jul 20, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.2015

Copy link
Contributor Author

@JohanEngelen JohanEngelen Jul 22, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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?).

Copy link
Member

@9il 9il Jul 22, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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:

  1. Example
  2. D code example for transformation writefln("%a", 0x2e9.221aa5a60a3p+0);

Copy link
Member

@9il 9il Jul 22, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Current constant is OK

Copy link
Member

@ibuclaw ibuclaw Jul 30, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Never believe what you read on Wikipedia. ;)

Just, kidding, this is fine.

}
else
static assert(0, "missing MAXLOG and MINLOG for other real types");

enum real BETA_BIG = 9.223372036854775808e18L;
enum real BETA_BIGINV = 1.084202172485504434007e-19L;
}
Expand Down