From 7f09d3487de97882585a72819ec7e0f9100c7121 Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Thu, 18 Feb 2021 20:49:18 +0000 Subject: Use the Cephes double subtraction trick in pexp even when FMA is available. Otherwise the accuracy drops from 1 ulp to 3 ulp. --- Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h | 15 ++++----------- 1 file changed, 4 insertions(+), 11 deletions(-) (limited to 'Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h') diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h index abdbcdbb9..e4c7e6223 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h @@ -1,4 +1,3 @@ - // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // @@ -468,16 +467,10 @@ Packet pexp_float(const Packet _x) // Get r = x - m*ln(2). If no FMA instructions are available, m*ln(2) is // subtracted out in two parts, m*C1+m*C2 = m*ln(2), to avoid accumulating // truncation errors. - Packet r; -#ifdef EIGEN_HAS_SINGLE_INSTRUCTION_MADD - const Packet cst_nln2 = pset1(-0.6931471805599453f); - r = pmadd(m, cst_nln2, x); -#else - const Packet cst_cephes_exp_C1 = pset1(0.693359375f); - const Packet cst_cephes_exp_C2 = pset1(-2.12194440e-4f); - r = psub(x, pmul(m, cst_cephes_exp_C1)); - r = psub(r, pmul(m, cst_cephes_exp_C2)); -#endif + const Packet cst_cephes_exp_C1 = pset1(-0.693359375f); + const Packet cst_cephes_exp_C2 = pset1(2.12194440e-4f); + Packet r = pmadd(m, cst_cephes_exp_C1, x); + r = pmadd(m, cst_cephes_exp_C2, r); Packet r2 = pmul(r, r); Packet r3 = pmul(r2, r); -- cgit v1.2.3