Using libquadmath with eigen

Question

I would like to use __float128 with eigen's erf() function, but found out that currently it only supports floats and doubles:

This function supports only float and double scalar types in c++11 mode. To support other scalar types, or float/double in non c++11 mode, the user has to provide implementations of erf(T) for any scalar type T to be supported.

As I want to use __float128, I want to rely on libquadmaths erfq implementation if that is possible. But how to do that? The only (ugly?) way I can currently think of is using eigens unaryExpr(). Are there other possibilities?

Answer 1

You can specialize Eigen::internal::erf_impl (similar for any other function of course):

#include <quadmath.h>
#include <iostream>
#include <unsupported/Eigen/SpecialFunctions>

namespace Eigen { namespace internal {
template<>
struct erf_impl<__float128> {
  EIGEN_DEVICE_FUNC
  static EIGEN_STRONG_INLINE __float128 run(__float128 x) { return ::erfq(x); }
};
}}

int main()
{
    typedef Eigen::Array<__float128, Eigen::Dynamic, 1> ArrayXF;
    ArrayXF a(4); a << 0, 0.25, 0.5, 0.75;
    ArrayXF b = a.erf();

    for(int i=0; i<4; ++i){
        char buf[100];
        quadmath_snprintf(buf, 100, "%.50Qe", b[i]); std::cout << buf << '\n';
    }
}

Output:

0.00000000000000000000000000000000000000000000000000e+00 2.76326390168236932985068267764815703534647315720851e-01 5.20499877813046537682746653891964513119913394193564e-01 7.11155633653515131598937834591410814324096358715387e-01