diff --git a/doc/CoeffwiseMathFunctionsTable.dox b/doc/CoeffwiseMathFunctionsTable.dox index 3b9ba7bbf..eb0f36c1b 100644 --- a/doc/CoeffwiseMathFunctionsTable.dox +++ b/doc/CoeffwiseMathFunctionsTable.dox @@ -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); - Incomplete beta function + regularized incomplete beta function built-in for float and double,\n but requires \cpp11 diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h index 387836b72..49c74d526 100644 --- a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h @@ -1844,14 +1844,15 @@ struct betainc_impl { const float nan = NumTraits::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 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::incbps(aa, bb, xx); diff --git a/unsupported/test/special_functions.cpp b/unsupported/test/special_functions.cpp index b16a7bb82..9239b3327 100644 --- a/unsupported/test/special_functions.cpp +++ b/unsupported/test/special_functions.cpp @@ -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);); }