131 lines · plain
1// (C) Copyright Matt Borland 2022.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_CCMATH_FMA_HPP7#define BOOST_MATH_CCMATH_FMA_HPP8 9#include <boost/math/ccmath/detail/config.hpp>10 11#ifdef BOOST_MATH_NO_CCMATH12#error "The header <boost/math/fma.hpp> can only be used in C++17 and later."13#endif14 15#include <boost/math/ccmath/isinf.hpp>16#include <boost/math/ccmath/isnan.hpp>17 18namespace boost::math::ccmath {19 20namespace detail {21 22template <typename T>23constexpr T fma_imp(const T x, const T y, const T z) noexcept24{25 #if defined(__GNUC__) && !defined(__clang__) && !defined(__INTEL_COMPILER) && !defined(__INTEL_LLVM_COMPILER)26 if constexpr (std::is_same_v<T, float>)27 {28 return __builtin_fmaf(x, y, z);29 }30 else if constexpr (std::is_same_v<T, double>)31 {32 return __builtin_fma(x, y, z);33 }34 else if constexpr (std::is_same_v<T, long double>)35 {36 return __builtin_fmal(x, y, z);37 }38 #endif39 40 // If we can't use compiler intrinsics hope that -fma flag optimizes this call to fma instruction41 return (x * y) + z;42}43 44} // Namespace detail45 46template <typename Real, std::enable_if_t<!std::is_integral_v<Real>, bool> = true>47constexpr Real fma(Real x, Real y, Real z) noexcept48{49 if (BOOST_MATH_IS_CONSTANT_EVALUATED(x))50 {51 if (x == 0 && boost::math::ccmath::isinf(y))52 {53 return std::numeric_limits<Real>::quiet_NaN();54 }55 else if (y == 0 && boost::math::ccmath::isinf(x))56 {57 return std::numeric_limits<Real>::quiet_NaN();58 }59 else if (boost::math::ccmath::isnan(x))60 {61 return std::numeric_limits<Real>::quiet_NaN();62 }63 else if (boost::math::ccmath::isnan(y))64 {65 return std::numeric_limits<Real>::quiet_NaN();66 }67 else if (boost::math::ccmath::isnan(z))68 {69 return std::numeric_limits<Real>::quiet_NaN();70 }71 72 return boost::math::ccmath::detail::fma_imp(x, y, z);73 }74 else75 {76 using std::fma;77 return fma(x, y, z);78 }79}80 81template <typename T1, typename T2, typename T3>82constexpr auto fma(T1 x, T2 y, T3 z) noexcept83{84 if (BOOST_MATH_IS_CONSTANT_EVALUATED(x))85 {86 // If the type is an integer (e.g. epsilon == 0) then set the epsilon value to 1 so that type is at a minimum 87 // cast to double88 constexpr auto T1p = std::numeric_limits<T1>::epsilon() > 0 ? std::numeric_limits<T1>::epsilon() : 1;89 constexpr auto T2p = std::numeric_limits<T2>::epsilon() > 0 ? std::numeric_limits<T2>::epsilon() : 1;90 constexpr auto T3p = std::numeric_limits<T3>::epsilon() > 0 ? std::numeric_limits<T3>::epsilon() : 1;91 92 using promoted_type = 93 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS94 std::conditional_t<T1p <= LDBL_EPSILON && T1p <= T2p, T1,95 std::conditional_t<T2p <= LDBL_EPSILON && T2p <= T1p, T2,96 std::conditional_t<T3p <= LDBL_EPSILON && T3p <= T2p, T3,97 #endif98 std::conditional_t<T1p <= DBL_EPSILON && T1p <= T2p, T1,99 std::conditional_t<T2p <= DBL_EPSILON && T2p <= T1p, T2, 100 std::conditional_t<T3p <= DBL_EPSILON && T3p <= T2p, T3, double101 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS102 >>>>>>;103 #else104 >>>;105 #endif106 107 return boost::math::ccmath::fma(promoted_type(x), promoted_type(y), promoted_type(z));108 }109 else110 {111 using std::fma;112 return fma(x, y, z);113 }114}115 116constexpr float fmaf(float x, float y, float z) noexcept117{118 return boost::math::ccmath::fma(x, y, z);119}120 121#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS122constexpr long double fmal(long double x, long double y, long double z) noexcept123{124 return boost::math::ccmath::fma(x, y, z);125}126#endif127 128} // Namespace boost::math::ccmath129 130#endif // BOOST_MATH_CCMATH_FMA_HPP131