Fix intermediate overflow in bessel_i0 and bessel_i1 for large finite arguments libeigen/eigen!2920
diff --git a/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsImpl.h b/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsImpl.h index 1b75542..05690d3 100644 --- a/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsImpl.h +++ b/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsImpl.h
@@ -180,7 +180,21 @@ template <typename T, typename ScalarType = typename unpacket_traits<T>::type> struct generic_i0 { EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T& x) { - return pmul(pexp(pabs(x)), generic_i0e<T, ScalarType>::run(x)); + // Evaluating i0(x) = exp(|x|) * i0e(x) can prematurely cause intermediate overflow for large |x| + // i.e. 88.7228 < |x| <= 91.9008 for float + // 709.7827 < |x| <= 713.9869 for double + // Instead, use i0(x) = exp(|x|/2) * (exp(|x|/2) * i0e(x)). Cutting |x| in half keeps the intermediate result + // finite for all finite results. + const T ax = pabs(x); + const T i0e = generic_i0e<T, ScalarType>::run(x); + const T half_exp = pexp(pmul(pset1<T>(ScalarType(0.5)), ax)); + T scaled = pmul(half_exp, i0e); +#if defined(__FAST_MATH__) || defined(__ASSOCIATIVE_MATH__) || EIGEN_COMP_NVHPC + // Unfortunately fast-math reassociates (exp(|x|/2) * exp(|x|/2)) * i0e(x) which reintroduces the overflow, + // so we need to tell it not to. + EIGEN_OPTIMIZATION_BARRIER(scaled) +#endif + return pmul(half_exp, scaled); } }; @@ -326,7 +340,20 @@ template <typename T, typename ScalarType = typename unpacket_traits<T>::type> struct generic_i1 { EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T& x) { - return pmul(pexp(pabs(x)), generic_i1e<T, ScalarType>::run(x)); + // Evaluating i1(x) = exp(|x|) * i1e(x) can prematurely cause intermediate overflow for large |x| + // i.e. 88.7228 < |x| <= 91.9063 for float + // 709.7827 < |x| <= 713.9876 for double + // Instead, use i1(x) = exp(|x|/2) * (exp(|x|/2) * i1e(x)). Cutting |x| in half keeps the intermediate result + // finite for all finite results. + const T ax = pabs(x); + const T i1e = generic_i1e<T, ScalarType>::run(x); + const T half_exp = pexp(pmul(pset1<T>(ScalarType(0.5)), ax)); + T scaled = pmul(half_exp, i1e); +#if defined(__FAST_MATH__) || defined(__ASSOCIATIVE_MATH__) || EIGEN_COMP_NVHPC + // Reassociating back to (exp(|x|/2) * exp(|x|/2)) * i1e(x) would reintroduce the overflow. + EIGEN_OPTIMIZATION_BARRIER(scaled) +#endif + return pmul(half_exp, scaled); } };
diff --git a/unsupported/test/bessel_functions.cpp b/unsupported/test/bessel_functions.cpp index f8d4833..b316f11 100644 --- a/unsupported/test/bessel_functions.cpp +++ b/unsupported/test/bessel_functions.cpp
@@ -270,7 +270,52 @@ } } +// exp(|x|) overflows above |x| ~ 88.7228 (float) / ~709.7827 (double), while i0/i1 stay finite up to +// |x| ~ 91.9008 / ~713.9869, so these arguments exercise the intermediate-overflow path. +// Reference values computed with mpmath at 40 digits. +template <typename ArrayType> +void test_bessel_i_large_finite(const typename ArrayType::Scalar (&x_val)[2], + const typename ArrayType::Scalar (&i0_truth)[2], + const typename ArrayType::Scalar (&i1_truth)[2]) { + ArrayType x(4); + x << x_val[0], -x_val[0], x_val[1], -x_val[1]; + ArrayType i0_res = bessel_i0(x); + ArrayType i1_res = bessel_i1(x); + for (Index i = 0; i < x.size(); ++i) { + VERIFY((numext::isfinite)(i0_res(i))); + VERIFY((numext::isfinite)(i1_res(i))); + // The vectorized and scalar paths must agree. + VERIFY_IS_APPROX(numext::bessel_i0(x(i)), i0_res(i)); + VERIFY_IS_APPROX(numext::bessel_i1(x(i)), i1_res(i)); + } + // i0 is even, i1 is odd. + VERIFY_IS_APPROX(i0_res(0), i0_truth[0]); + VERIFY_IS_APPROX(i0_res(1), i0_truth[0]); + VERIFY_IS_APPROX(i0_res(2), i0_truth[1]); + VERIFY_IS_APPROX(i0_res(3), i0_truth[1]); + VERIFY_IS_APPROX(i1_res(0), i1_truth[0]); + VERIFY_IS_APPROX(i1_res(1), -i1_truth[0]); + VERIFY_IS_APPROX(i1_res(2), i1_truth[1]); + VERIFY_IS_APPROX(i1_res(3), -i1_truth[1]); +} + +void test_bessel_i_large_finite_float() { + const float x[2] = {90.0f, 91.0f}; + const float i0_truth[2] = {5.1392383455086638e+37f, 1.3892714060989622e+38f}; + const float i1_truth[2] = {5.1106068152565982e+37f, 1.3816168414593216e+38f}; + test_bessel_i_large_finite<ArrayXf>(x, i0_truth, i1_truth); +} + +void test_bessel_i_large_finite_double() { + const double x[2] = {712.0, 713.0}; + const double i0_truth[2] = {2.4684110577627524e+307, 6.7051282636709964e+307}; + const double i1_truth[2] = {2.4666770135246152e+307, 6.7004245591864022e+307}; + test_bessel_i_large_finite<ArrayXd>(x, i0_truth, i1_truth); +} + EIGEN_DECLARE_TEST(bessel_functions) { CALL_SUBTEST_1(array_bessel_functions<ArrayXf>()); + CALL_SUBTEST_1(test_bessel_i_large_finite_float()); CALL_SUBTEST_2(array_bessel_functions<ArrayXd>()); + CALL_SUBTEST_2(test_bessel_i_large_finite_double()); }
diff --git a/unsupported/test/special_packetmath.cpp b/unsupported/test/special_packetmath.cpp index 8bf1a40..9b500af 100644 --- a/unsupported/test/special_packetmath.cpp +++ b/unsupported/test/special_packetmath.cpp
@@ -143,6 +143,21 @@ CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i0, internal::pbessel_i0); CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i1, internal::pbessel_i1); + // Boundary values for which a naive i0(x) = exp(|x|) * i0e(x) overflows even though the result is + // finite. CHECK_CWISE1_IF cannot detect that on its own: its reference is the scalar path through + // the same generic_i0, so both sides would be +inf and compare equal. Check finiteness explicitly. + if (PacketTraits::HasBessel && + (internal::is_same<Scalar, float>::value || internal::is_same<Scalar, double>::value)) { + Scalar boundary = internal::is_same<Scalar, float>::value ? Scalar(90) : Scalar(713); + for (int i = 0; i < size; ++i) { + data1[i] = (i % 2 == 0) ? boundary : -boundary; + } + CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i0, internal::pbessel_i0); + for (int i = 0; i < PacketSize; ++i) VERIFY((numext::isfinite)(data2[i])); + CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i1, internal::pbessel_i1); + for (int i = 0; i < PacketSize; ++i) VERIFY((numext::isfinite)(data2[i])); + } + // y_i, and k_i are valid for x > 0. { const int max_exponent = numext::mini(std::numeric_limits<Scalar>::max_exponent10 - 1, 5);