diff options
Diffstat (limited to 'mathfuncs_log.h')
| -rw-r--r-- | mathfuncs_log.h | 137 |
1 files changed, 63 insertions, 74 deletions
diff --git a/mathfuncs_log.h b/mathfuncs_log.h index cd71eb3..fa517ba 100644 --- a/mathfuncs_log.h +++ b/mathfuncs_log.h @@ -7,93 +7,82 @@ #include <cmath> +namespace vecmathlib { +template <typename realvec_t> +realvec_t mathfuncs<realvec_t>::vml_log2(realvec_t x) { + // Algorithm inspired by SLEEF 2.80 -namespace vecmathlib { - - template<typename realvec_t> - realvec_t mathfuncs<realvec_t>::vml_log2(realvec_t x) - { - // Algorithm inspired by SLEEF 2.80 - - // Rescale - intvec_t ilogb_x = ilogb(x * RV(M_SQRT2)); - x = ldexp(x, -ilogb_x); - VML_ASSERT(all(x >= RV(M_SQRT1_2) && x <= RV(M_SQRT2))); - - realvec_t y = (x - RV(1.0)) / (x + RV(1.0)); - realvec_t y2 = y*y; - - realvec_t r; - switch (sizeof(real_t)) { - case 4: - // float, error=7.09807175879142775648452461821e-8 - r = RV(0.59723611417135718739797302426); - r = mad(r, y2, RV(0.961524413175528426101613434)); - r = mad(r, y2, RV(2.88539097665498228703236701)); - break; - case 8: + // Rescale + intvec_t ilogb_x = ilogb(x * RV(M_SQRT2)); + x = ldexp(x, -ilogb_x); + VML_ASSERT(all(x >= RV(M_SQRT1_2) && x <= RV(M_SQRT2))); + + realvec_t y = (x - RV(1.0)) / (x + RV(1.0)); + realvec_t y2 = y * y; + + realvec_t r; + switch (sizeof(real_t)) { + case 4: + // float, error=7.09807175879142775648452461821e-8 + r = RV(0.59723611417135718739797302426); + r = mad(r, y2, RV(0.961524413175528426101613434)); + r = mad(r, y2, RV(2.88539097665498228703236701)); + break; + case 8: #ifdef VML_HAVE_FP_CONTRACT - // double, error=1.48294180185938512675770096324e-16 - r = RV(0.243683403415639178527756320773); - r = mad(r, y2, RV(0.26136626803870009948502658)); - r = mad(r, y2, RV(0.320619429891299265439389)); - r = mad(r, y2, RV(0.4121983452028499242926)); - r = mad(r, y2, RV(0.577078017761894161436)); - r = mad(r, y2, RV(0.96179669392233355927)); - r = mad(r, y2, RV(2.8853900817779295236)); + // double, error=1.48294180185938512675770096324e-16 + r = RV(0.243683403415639178527756320773); + r = mad(r, y2, RV(0.26136626803870009948502658)); + r = mad(r, y2, RV(0.320619429891299265439389)); + r = mad(r, y2, RV(0.4121983452028499242926)); + r = mad(r, y2, RV(0.577078017761894161436)); + r = mad(r, y2, RV(0.96179669392233355927)); + r = mad(r, y2, RV(2.8853900817779295236)); #else - // double, error=2.1410114030383689267772704676e-14 - r = RV(0.283751646449323373643963474845); - r = mad(r, y2, RV(0.31983138095551191299118812)); - r = mad(r, y2, RV(0.412211603844146279666022)); - r = mad(r, y2, RV(0.5770779098948940070516)); - r = mad(r, y2, RV(0.961796694295973716912)); - r = mad(r, y2, RV(2.885390081777562819196)); + // double, error=2.1410114030383689267772704676e-14 + r = RV(0.283751646449323373643963474845); + r = mad(r, y2, RV(0.31983138095551191299118812)); + r = mad(r, y2, RV(0.412211603844146279666022)); + r = mad(r, y2, RV(0.5770779098948940070516)); + r = mad(r, y2, RV(0.961796694295973716912)); + r = mad(r, y2, RV(2.885390081777562819196)); #endif - break; - default: - __builtin_unreachable(); - } - r *= y; - - // Undo rescaling - r += convert_float(ilogb_x); - - return r; - } - - - - template<typename realvec_t> - inline - realvec_t mathfuncs<realvec_t>::vml_log(realvec_t x) - { - return log2(x) * RV(M_LN2); + break; + default: + __builtin_unreachable(); } + r *= y; - template<typename realvec_t> - inline - realvec_t mathfuncs<realvec_t>::vml_log10(realvec_t x) - { - return log(x) * RV(M_LOG10E); - } + // Undo rescaling + r += convert_float(ilogb_x); - template<typename realvec_t> - inline - realvec_t mathfuncs<realvec_t>::vml_log1p(realvec_t x) - { - // TODO: Check SLEEF 2.80 algorithm - - return log(RV(1.0) + x); + return r; +} + +template <typename realvec_t> +inline realvec_t mathfuncs<realvec_t>::vml_log(realvec_t x) { + return log2(x) * RV(M_LN2); +} + +template <typename realvec_t> +inline realvec_t mathfuncs<realvec_t>::vml_log10(realvec_t x) { + return log(x) * RV(M_LOG10E); +} + +template <typename realvec_t> +inline realvec_t mathfuncs<realvec_t>::vml_log1p(realvec_t x) { + // TODO: Check SLEEF 2.80 algorithm + + return log(RV(1.0) + x); #if 0 // Goldberg, theorem 4 realvec_t x1 = RV(1.0) + x; x1.barrier(); return ifthen(x1 == x, x, x * log(x1) / (x1 - RV(1.0))); #endif - } - +} + }; // namespace vecmathlib -#endif // #ifndef MATHFUNCS_LOG_H +#endif // #ifndef MATHFUNCS_LOG_H |
