phynexis-deps/eigen
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

SpecialFunctions: Restore erf/erfc for custom scalars and long double

Browse Source 
libeigen/eigen!2840

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
This commit is contained in:
Rasmus Munk Larsen
2026-08-17 19:32:40 -07:00
co-authored by Rasmus Munk Larsen
parent b0d40c93fc
commit e70059f86a
6 changed files with 201 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files
failtest/CMakeLists.txt
+3
View File
@@ -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")
failtest/erf_no_scalar_overload.cpp
+24
View File
@@ -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<custom_scalar::CustomReal> : NumTraits<double> {
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; }
failtest/erfc_no_scalar_overload.cpp
+24
View File
@@ -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<custom_scalar::CustomReal> : NumTraits<double> {
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; }
unsupported/Eigen/SpecialFunctions
+6
View File
@@ -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 <unsupported/Eigen/SpecialFunctions>
* \endcode
unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
+73 -2
View File
@@ -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 <typename T>
auto test_erf(int) -> decltype(void(erf(std::declval<const T&>())), std::true_type{});
template <typename T>
std::false_type test_erf(...);
template <typename T>
auto test_erfc(int) -> decltype(void(erfc(std::declval<const T&>())), std::true_type{});
template <typename T>
std::false_type test_erfc(...);
} // namespace unqualified_erf
template <typename T>
struct has_erf : decltype(unqualified_erf::test_erf<T>(0)) {};
template <typename T>
struct has_erfc : decltype(unqualified_erf::test_erfc<T>(0)) {};
/***************************************************************************
* Implementation of erfc.
****************************************************************************/
@@ -416,7 +437,32 @@ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T generic_fast_erfc<double>::run(const T&
template <typename T>
struct erfc_impl {
typedef typename unpacket_traits<T>::type Scalar;
EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T& x) { return generic_fast_erfc<Scalar>::run(x); }
EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T& x) { return run_impl(x, std::is_same<T, Scalar>()); }
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<Scalar>::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<T>());
}
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<T>::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<double>::run(const T& x
template <typename T>
struct erf_impl {
typedef typename unpacket_traits<T>::type Scalar;
EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T& x) { return generic_fast_erf<Scalar>::run(x); }
EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T& x) { return run_impl(x, std::is_same<T, Scalar>()); }
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<Scalar>::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<T>());
}
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<T>::value, SCALAR_TYPE_MUST_PROVIDE_AN_ERF_OVERLOAD_FOUND_BY_ADL_OR_IN_NAMESPACE_STD)
return x;
}
};
template <>
unsupported/test/special_functions.cpp
+71
View File
@@ -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 <typename Scalar>
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<custom_scalar::CustomReal> : NumTraits<double> {
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<CustomReal>::value));
STATIC_CHECK((internal::has_erfc<CustomReal>::value));
STATIC_CHECK((internal::has_erf<long double>::value));
STATIC_CHECK((internal::has_erfc<long double>::value));
STATIC_CHECK((!internal::has_erf<NoErfReal>::value));
STATIC_CHECK((!internal::has_erfc<NoErfReal>::value));
Eigen::Array<CustomReal, 4, 1> x;
x(0) = CustomReal(-2.5);
x(1) = CustomReal(-0.5);
x(2) = CustomReal(0.25);
x(3) = CustomReal(3.0);
Eigen::Array<CustomReal, 4, 1> e = x.erf();
Eigen::Array<CustomReal, 4, 1> 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<ArrayXf>());
CALL_SUBTEST_2(array_special_functions<ArrayXd>());
CALL_SUBTEST_3(scalar_ndtri<long double>());
CALL_SUBTEST_4(scalar_erf_erfc<long double>());
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<ArrayX<Eigen::half>>());
// CALL_SUBTEST_5(array_special_functions<ArrayX<Eigen::bfloat16>>());