devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

betainc edge case checks at start of calculation

Browse Source 
libeigen/eigen!2123

Closes #2359
This commit is contained in:
Blake
2026-02-08 13:05:06 -05:00
committed by Rasmus Munk Larsen
parent afb4380534
commit 752911927f
3 changed files with 34 additions and 35 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/CoeffwiseMathFunctionsTable.dox
View File

@@ -572,7 +572,7 @@ This also means that, unless specified, if the function \c std::foo is available
\anchor cwisetable_betainc
betainc(a,b,x);
</td>
<td><a href="https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function">Incomplete beta function</a></td>
<td><a href="https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function">regularized incomplete beta function</a></td>
<td>
built-in for float and double,\n but requires \cpp11
</td>

36
unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
View File

@@ -1844,14 +1844,15 @@ struct betainc_impl<float> {
const float nan = NumTraits<float>::quiet_NaN();
float ans, t;
if (a <= 0.0f) return nan;
if (b <= 0.0f) return nan;
if ((x <= 0.0f) || (x >= 1.0f)) {
if (x == 0.0f) return 0.0f;
if (x == 1.0f) return 1.0f;
// mtherr("betaincf", DOMAIN);
return nan;
}
if (a == 0.0f && b == 0.0f) return nan;
if (x < 0.0f || x > 1.0f) return nan;
if (a < 0.0f) return nan;
if (b < 0.0f) return nan;
if (a == 0.0f) return 1.0f;
if (b == 0.0f) return 0.0f;
if (x == 0.0f) return 0.0f;
if (x == 1.0f) return 1.0f;
// mtherr("betaincf", DOMAIN);
/* transformation for small aa */
if (a <= 1.0f) {
@@ -1914,16 +1915,15 @@ struct betainc_impl<double> {
double a, b, t, x, xc, w, y;
bool reversed_a_b = false;
if (aa <= 0.0 || bb <= 0.0) {
return nan; // goto domerr;
}
if ((xx <= 0.0) || (xx >= 1.0)) {
if (xx == 0.0) return (0.0);
if (xx == 1.0) return (1.0);
// mtherr("incbet", DOMAIN);
return nan;
}
if (aa == 0.0 && bb == 0.0) return nan;
if (xx < 0.0 || xx > 1.0) return nan;
if (aa < 0.0) return nan;
if (bb < 0.0) return nan;
if (aa == 0.0) return 1.0;
if (bb == 0.0) return 0.0;
if (xx == 0.0) return 0.0;
if (xx == 1.0) return 1.0;
// mtherr("incbet", DOMAIN);
if ((bb * xx) <= 1.0 && xx <= 0.95) {
return betainc_helper<double>::incbps(aa, bb, xx);

31
unsupported/test/special_functions.cpp
View File

@@ -223,10 +223,11 @@ void array_special_functions() {
// b = np.logspace(-3, 3, 5) - 1e-3
// x = np.linspace(-0.1, 1.1, 5)
// (full_a, full_b, full_x) = np.vectorize(lambda a, b, x: (a, b, x))(*np.ix_(a, b, x))
// full_a = full_a.flatten().tolist() # same for full_b, full_x
// full_a = full_a.flatten().tolist()
// full_b = full_b.flatten().tolist()
// full_x = full_x.flatten().tolist()
// v = scipy.special.betainc(full_a, full_b, full_x).flatten().tolist()
//
// Note in Eigen, we call betainc with arguments in the order (x, a, b).
// Note in Eigen, we call betainc with arguments in the order (a, b, x);
ArrayType a(125);
ArrayType b(125);
ArrayType x(125);
@@ -272,19 +273,17 @@ void array_special_functions() {
1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1,
-0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1;
v << nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
nan, nan, nan, nan, nan, nan, nan, nan, nan, 0.47972119876364683, 0.5, 0.5202788012363533, nan, nan,
0.9518683957740043, 0.9789663010413743, 0.9931729188073435, nan, nan, 0.999995949033062, 0.9999999999993698,
0.9999999999999999, nan, nan, 0.9999999999999999, 0.9999999999999999, 0.9999999999999999, nan, nan, nan, nan,
nan, nan, nan, 0.006827081192655869, 0.0210336989586256, 0.04813160422599567, nan, nan, 0.20014344256217678,
0.5000000000000001, 0.7998565574378232, nan, nan, 0.9991401428435834, 0.999999999698403, 0.9999999999999999,
nan, nan, 0.9999999999999999, 0.9999999999999999, 0.9999999999999999, nan, nan, nan, nan, nan, nan, nan,
1.0646600232370887e-25, 6.301722877826246e-13, 4.050966937974938e-06, nan, nan, 7.864342668429763e-23,
3.015969667594166e-10, 0.0008598571564165444, nan, nan, 6.031987710123844e-08, 0.5000000000000007,
0.9999999396801229, nan, nan, 0.9999999999999999, 0.9999999999999999, 0.9999999999999999, nan, nan, nan, nan,
nan, nan, nan, 0.0, 7.029920380986636e-306, 2.2450728208591345e-101, nan, nan, 0.0, 9.275871147869727e-302,
1.2232913026152827e-97, nan, nan, 0.0, 3.0891393081932924e-252, 2.9303043666183996e-60, nan, nan,
2.248913486879199e-196, 0.5000000000004947, 0.9999999999999999, nan;
v << nan, nan, nan, nan, nan, nan, 1.0, 1.0, 1.0, nan, nan, 1.0, 1.0, 1.0, nan, nan, 1.0, 1.0, 1.0, nan, nan, 1.0,
1.0, 1.0, nan, nan, 0.0, 0.0, 0.0, nan, nan, 0.47972119876364683, 0.4999999999999997, 0.5202788012363533, nan,
nan, 0.9518683957740043, 0.9789663010413745, 0.993172918807344, nan, nan, 0.999995949033062, 0.99999999999937,
1.0, nan, nan, 1.0, 1.0, 1.0, nan, nan, 0.0, 0.0, 0.0, nan, nan, 0.0068270811926558735, 0.0210336989586256,
0.048131604225995696, nan, nan, 0.20014344256217687, 0.5000000000000002, 0.7998565574378237, nan, nan,
0.9991401428435834, 0.9999999996984031, 1.0, nan, nan, 1.0, 1.0, 1.0, nan, nan, 0.0, 0.0, 0.0, nan, nan,
1.0646600232370862e-25, 6.30172287782625e-13, 4.050966937974916e-06, nan, nan, 7.864342668429712e-23,
3.015969667594142e-10, 0.0008598571564165379, nan, nan, 6.031987710123845e-08, 0.5000000000000001,
0.9999999396801229, nan, nan, 1.0, 1.0, 1.0, nan, nan, 0.0, 0.0, 0.0, nan, nan, 0.0, 7.029920380993552e-306,
2.245072820861537e-101, nan, nan, 0.0, 9.275871147875753e-302, 1.223291302616039e-97, nan, nan, 0.0,
3.089139308197642e-252, 2.9303043666230076e-60, nan, nan, 2.248913486881819e-196, 0.5000000000000036, 1.0, nan;
CALL_SUBTEST(res = betainc(a, b, x); verify_component_wise(res, v););
}