diff --git a/failtest/CMakeLists.txt b/failtest/CMakeLists.txt index aac11efc8..a9f946731 100644 --- a/failtest/CMakeLists.txt +++ b/failtest/CMakeLists.txt @@ -69,6 +69,9 @@ ei_add_failtest("eigensolver_cplx") ei_add_failtest("initializer_list_1") ei_add_failtest("initializer_list_2") +ei_add_failtest("erf_no_scalar_overload") +ei_add_failtest("erfc_no_scalar_overload") + ei_add_failtest("structured_bindings_dynamic_matrix") ei_add_failtest("structured_bindings_dynamic_array") ei_add_failtest("structured_bindings_rowmajor") diff --git a/failtest/erf_no_scalar_overload.cpp b/failtest/erf_no_scalar_overload.cpp new file mode 100644 index 000000000..9255dc70f --- /dev/null +++ b/failtest/erf_no_scalar_overload.cpp @@ -0,0 +1,24 @@ +#include "../unsupported/Eigen/SpecialFunctions" + +// A custom scalar that does not provide erf. Eigen's approximations are tuned +// for float and double and are not valid for it, so numext::erf must be +// rejected at compile time instead of computing a wrong result. +namespace custom_scalar { +struct CustomReal { + double value; +}; +#ifndef EIGEN_SHOULD_FAIL_TO_BUILD +inline CustomReal erf(const CustomReal& x) { return CustomReal{std::erf(x.value)}; } +#endif +} // namespace custom_scalar + +namespace Eigen { +template <> +struct NumTraits : NumTraits { + using Real = custom_scalar::CustomReal; + using NonInteger = custom_scalar::CustomReal; + using Nested = custom_scalar::CustomReal; +}; +} // namespace Eigen + +int main() { return Eigen::numext::erf(custom_scalar::CustomReal{0.5}).value > 0.0 ? 0 : 1; } diff --git a/failtest/erfc_no_scalar_overload.cpp b/failtest/erfc_no_scalar_overload.cpp new file mode 100644 index 000000000..57e9af291 --- /dev/null +++ b/failtest/erfc_no_scalar_overload.cpp @@ -0,0 +1,24 @@ +#include "../unsupported/Eigen/SpecialFunctions" + +// A custom scalar that does not provide erfc. Eigen's approximations are tuned +// for float and double and are not valid for it, so numext::erfc must be +// rejected at compile time instead of computing a wrong result. +namespace custom_scalar { +struct CustomReal { + double value; +}; +#ifndef EIGEN_SHOULD_FAIL_TO_BUILD +inline CustomReal erfc(const CustomReal& x) { return CustomReal{std::erfc(x.value)}; } +#endif +} // namespace custom_scalar + +namespace Eigen { +template <> +struct NumTraits : NumTraits { + using Real = custom_scalar::CustomReal; + using NonInteger = custom_scalar::CustomReal; + using Nested = custom_scalar::CustomReal; +}; +} // namespace Eigen + +int main() { return Eigen::numext::erfc(custom_scalar::CustomReal{0.5}).value > 0.0 ? 0 : 1; } diff --git a/unsupported/Eigen/SpecialFunctions b/unsupported/Eigen/SpecialFunctions index c8a209f57..70520549d 100644 --- a/unsupported/Eigen/SpecialFunctions +++ b/unsupported/Eigen/SpecialFunctions @@ -53,6 +53,12 @@ namespace Eigen { * - bessel_y0 * - bessel_y1 * + * The implementations are tuned for \c float and \c double, so they do not carry + * over to an arbitrary scalar. \c erf and \c erfc accept any other scalar type + * that declares its own \c erf / \c erfc in the namespace of that type, where + * argument-dependent lookup finds it; a scalar type that declares neither is + * rejected at compile time. + * * \code * #include * \endcode diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h index fd19f6ee4..49724a52f 100644 --- a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h @@ -259,6 +259,27 @@ struct digamma_impl { } }; +// Does unqualified lookup of erf/erfc succeed for T? The lookup below mirrors +// the one at the call sites in erf_impl/erfc_impl, so it finds std::erf/std::erfc +// and any overload visible through argument-dependent lookup. +namespace unqualified_erf { +EIGEN_USING_STD(erf) +EIGEN_USING_STD(erfc) +template +auto test_erf(int) -> decltype(void(erf(std::declval())), std::true_type{}); +template +std::false_type test_erf(...); +template +auto test_erfc(int) -> decltype(void(erfc(std::declval())), std::true_type{}); +template +std::false_type test_erfc(...); +} // namespace unqualified_erf + +template +struct has_erf : decltype(unqualified_erf::test_erf(0)) {}; +template +struct has_erfc : decltype(unqualified_erf::test_erfc(0)) {}; + /*************************************************************************** * Implementation of erfc. ****************************************************************************/ @@ -416,7 +437,32 @@ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T generic_fast_erfc::run(const T& template struct erfc_impl { typedef typename unpacket_traits::type Scalar; - EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T& x) { return generic_fast_erfc::run(x); } + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T& x) { return run_impl(x, std::is_same()); } + + private: + // Packets of float/double: vectorized rational approximation. + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run_impl(const T& x, std::false_type) { + return generic_fast_erfc::run(x); + } + // Any other scalar type: defer to an erfc found by argument-dependent lookup + // (or std::erfc), keeping custom scalars on their own implementation instead + // of the float/double-tuned polynomials. + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run_impl(const T& x, std::true_type) { + EIGEN_STATIC_ASSERT_NON_INTEGER(T) + return run_scalar(x, has_erfc()); + } + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run_scalar(const T& x, std::true_type) { + EIGEN_USING_STD(erfc); + return erfc(x); + } + // Reject the type here instead of letting overload resolution fail inside the + // call above: the approximations in this file are tuned for float and double + // and are not valid for an arbitrary scalar, so a scalar type that wants erfc + // has to supply it. + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run_scalar(const T& x, std::false_type) { + EIGEN_STATIC_ASSERT(has_erfc::value, SCALAR_TYPE_MUST_PROVIDE_AN_ERFC_OVERLOAD_FOUND_BY_ADL_OR_IN_NAMESPACE_STD) + return x; + } }; template <> @@ -511,7 +557,32 @@ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T generic_fast_erf::run(const T& x template struct erf_impl { typedef typename unpacket_traits::type Scalar; - EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T& x) { return generic_fast_erf::run(x); } + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T& x) { return run_impl(x, std::is_same()); } + + private: + // Packets of float/double: vectorized rational approximation. + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run_impl(const T& x, std::false_type) { + return generic_fast_erf::run(x); + } + // Any other scalar type: defer to an erf found by argument-dependent lookup + // (or std::erf), keeping custom scalars on their own implementation instead + // of the float/double-tuned polynomials. + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run_impl(const T& x, std::true_type) { + EIGEN_STATIC_ASSERT_NON_INTEGER(T) + return run_scalar(x, has_erf()); + } + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run_scalar(const T& x, std::true_type) { + EIGEN_USING_STD(erf); + return erf(x); + } + // Reject the type here instead of letting overload resolution fail inside the + // call above: the approximations in this file are tuned for float and double + // and are not valid for an arbitrary scalar, so a scalar type that wants erf + // has to supply it. + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run_scalar(const T& x, std::false_type) { + EIGEN_STATIC_ASSERT(has_erf::value, SCALAR_TYPE_MUST_PROVIDE_AN_ERF_OVERLOAD_FOUND_BY_ADL_OR_IN_NAMESPACE_STD) + return x; + } }; template <> diff --git a/unsupported/test/special_functions.cpp b/unsupported/test/special_functions.cpp index ab9721b9d..c1702a0e5 100644 --- a/unsupported/test/special_functions.cpp +++ b/unsupported/test/special_functions.cpp @@ -441,10 +441,81 @@ void scalar_ndtri() { VERIFY(positive_inf > Scalar(0)); } +// A scalar type without a vectorized erf/erfc implementation must defer to +// std::erf/std::erfc instead of the float/double-tuned polynomials. +template +void scalar_erf_erfc() { + EIGEN_USING_STD(erf); + EIGEN_USING_STD(erfc); + // Approximate comparison: the compiler may constant-fold one side to the + // correctly rounded value while the other runs the libm implementation. + for (Scalar x : {Scalar(-2.5L), Scalar(-0.5L), Scalar(0.0L), Scalar(0.25L), Scalar(3.0L)}) { + VERIFY_IS_APPROX(numext::erf(x), Scalar(erf(x))); + VERIFY_IS_APPROX(numext::erfc(x), Scalar(erfc(x))); + } +} + +namespace custom_scalar { +// Minimal custom real scalar providing erf/erfc via argument-dependent lookup, +// the way multiprecision types do. numext::erf/erfc must route to these instead +// of the float/double-tuned polynomials (issue #3023). +struct CustomReal { + double value; + CustomReal() : value(0) {} + explicit CustomReal(double v) : value(v) {} +}; +inline CustomReal erf(const CustomReal& x) { return CustomReal(std::erf(x.value)); } +inline CustomReal erfc(const CustomReal& x) { return CustomReal(std::erfc(x.value)); } + +// The shape of boost::numeric::interval: neither erf nor erfc. Eigen's +// approximations are not valid for such a type, so it must not be routed into +// them; see failtest/erf_no_scalar_overload.cpp for the resulting error. +struct NoErfReal { + double value; +}; +} // namespace custom_scalar + +namespace Eigen { +template <> +struct NumTraits : NumTraits { + using Real = custom_scalar::CustomReal; + using NonInteger = custom_scalar::CustomReal; + using Nested = custom_scalar::CustomReal; +}; +} // namespace Eigen + +void custom_scalar_erf_erfc() { + using custom_scalar::CustomReal; + using custom_scalar::NoErfReal; + // Pin the routing decision in both directions. + STATIC_CHECK((internal::has_erf::value)); + STATIC_CHECK((internal::has_erfc::value)); + STATIC_CHECK((internal::has_erf::value)); + STATIC_CHECK((internal::has_erfc::value)); + STATIC_CHECK((!internal::has_erf::value)); + STATIC_CHECK((!internal::has_erfc::value)); + + Eigen::Array x; + x(0) = CustomReal(-2.5); + x(1) = CustomReal(-0.5); + x(2) = CustomReal(0.25); + x(3) = CustomReal(3.0); + Eigen::Array e = x.erf(); + Eigen::Array c = x.erfc(); + // The polynomial paths cannot even be instantiated for CustomReal, so + // compiling proves the routing; the value check guards the plumbing. + for (int i = 0; i < 4; ++i) { + VERIFY_IS_APPROX(e(i).value, std::erf(x(i).value)); + VERIFY_IS_APPROX(c(i).value, std::erfc(x(i).value)); + } +} + EIGEN_DECLARE_TEST(special_functions) { CALL_SUBTEST_1(array_special_functions()); CALL_SUBTEST_2(array_special_functions()); CALL_SUBTEST_3(scalar_ndtri()); + CALL_SUBTEST_4(scalar_erf_erfc()); + CALL_SUBTEST_5(custom_scalar_erf_erfc()); // TODO(cantonios): half/bfloat16 don't have enough precision to reproduce results above. // CALL_SUBTEST_4(array_special_functions>()); // CALL_SUBTEST_5(array_special_functions>());