math/ccmath/fma.hpp

0001 //  (C) Copyright Matt Borland 2022.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_CCMATH_FMA_HPP
0007 #define BOOST_MATH_CCMATH_FMA_HPP
0008
0009 #include <boost/math/ccmath/detail/config.hpp>
0010
0011 #ifdef BOOST_MATH_NO_CCMATH
0012 #error "The header <boost/math/fma.hpp> can only be used in C++17 and later."
0013 #endif
0014
0015 #include <boost/math/ccmath/isinf.hpp>
0016 #include <boost/math/ccmath/isnan.hpp>
0017
0018 namespace boost::math::ccmath {
0019
0020 namespace detail {
0021
0022 template <typename T>
0023 constexpr T fma_imp(const T x, const T y, const T z) noexcept
0024 {
0025     #if defined(__GNUC__) && !defined(__clang__) && !defined(__INTEL_COMPILER) && !defined(__INTEL_LLVM_COMPILER)
0026     if constexpr (std::is_same_v<T, float>)
0027     {
0028         return __builtin_fmaf(x, y, z);
0029     }
0030     else if constexpr (std::is_same_v<T, double>)
0031     {
0032         return __builtin_fma(x, y, z);
0033     }
0034     else if constexpr (std::is_same_v<T, long double>)
0035     {
0036         return __builtin_fmal(x, y, z);
0037     }
0038     #endif
0039
0040     // If we can't use compiler intrinsics hope that -fma flag optimizes this call to fma instruction
0041     return (x * y) + z;
0042 }
0043
0044 } // Namespace detail
0045
0046 template <typename Real, std::enable_if_t<!std::is_integral_v<Real>, bool> = true>
0047 constexpr Real fma(Real x, Real y, Real z) noexcept
0048 {
0049     if (BOOST_MATH_IS_CONSTANT_EVALUATED(x))
0050     {
0051         if (x == 0 && boost::math::ccmath::isinf(y))
0052         {
0053             return std::numeric_limits<Real>::quiet_NaN();
0054         }
0055         else if (y == 0 && boost::math::ccmath::isinf(x))
0056         {
0057             return std::numeric_limits<Real>::quiet_NaN();
0058         }
0059         else if (boost::math::ccmath::isnan(x))
0060         {
0061             return std::numeric_limits<Real>::quiet_NaN();
0062         }
0063         else if (boost::math::ccmath::isnan(y))
0064         {
0065             return std::numeric_limits<Real>::quiet_NaN();
0066         }
0067         else if (boost::math::ccmath::isnan(z))
0068         {
0069             return std::numeric_limits<Real>::quiet_NaN();
0070         }
0071
0072         return boost::math::ccmath::detail::fma_imp(x, y, z);
0073     }
0074     else
0075     {
0076         using std::fma;
0077         return fma(x, y, z);
0078     }
0079 }
0080
0081 template <typename T1, typename T2, typename T3>
0082 constexpr auto fma(T1 x, T2 y, T3 z) noexcept
0083 {
0084     if (BOOST_MATH_IS_CONSTANT_EVALUATED(x))
0085     {
0086         // If the type is an integer (e.g. epsilon == 0) then set the epsilon value to 1 so that type is at a minimum
0087         // cast to double
0088         constexpr auto T1p = std::numeric_limits<T1>::epsilon() > 0 ? std::numeric_limits<T1>::epsilon() : 1;
0089         constexpr auto T2p = std::numeric_limits<T2>::epsilon() > 0 ? std::numeric_limits<T2>::epsilon() : 1;
0090         constexpr auto T3p = std::numeric_limits<T3>::epsilon() > 0 ? std::numeric_limits<T3>::epsilon() : 1;
0091
0092         using promoted_type =
0093                               #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
0094                               std::conditional_t<T1p <= LDBL_EPSILON && T1p <= T2p, T1,
0095                               std::conditional_t<T2p <= LDBL_EPSILON && T2p <= T1p, T2,
0096                               std::conditional_t<T3p <= LDBL_EPSILON && T3p <= T2p, T3,
0097                               #endif
0098                               std::conditional_t<T1p <= DBL_EPSILON && T1p <= T2p, T1,
0099                               std::conditional_t<T2p <= DBL_EPSILON && T2p <= T1p, T2,
0100                               std::conditional_t<T3p <= DBL_EPSILON && T3p <= T2p, T3, double
0101                               #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
0102                               >>>>>>;
0103                               #else
0104                               >>>;
0105                               #endif
0106
0107         return boost::math::ccmath::fma(promoted_type(x), promoted_type(y), promoted_type(z));
0108     }
0109     else
0110     {
0111         using std::fma;
0112         return fma(x, y, z);
0113     }
0114 }
0115
0116 constexpr float fmaf(float x, float y, float z) noexcept
0117 {
0118     return boost::math::ccmath::fma(x, y, z);
0119 }
0120
0121 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
0122 constexpr long double fmal(long double x, long double y, long double z) noexcept
0123 {
0124     return boost::math::ccmath::fma(x, y, z);
0125 }
0126 #endif
0127
0128 } // Namespace boost::math::ccmath
0129
0130 #endif // BOOST_MATH_CCMATH_FMA_HPP