fma.hpp 4.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136
  1. // (C) Copyright Matt Borland 2022.
  2. // Use, modification and distribution are subject to the
  3. // Boost Software License, Version 1.0. (See accompanying file
  4. // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  5. #ifndef BOOST_MATH_CCMATH_FMA_HPP
  6. #define BOOST_MATH_CCMATH_FMA_HPP
  7. #include <cmath>
  8. #include <limits>
  9. #include <type_traits>
  10. #include <boost/math/tools/is_constant_evaluated.hpp>
  11. #include <boost/math/ccmath/isinf.hpp>
  12. #include <boost/math/ccmath/isnan.hpp>
  13. #include <boost/math/tools/is_standalone.hpp>
  14. #ifndef BOOST_MATH_STANDALONE
  15. #include <boost/config.hpp>
  16. #ifdef BOOST_NO_CXX17_IF_CONSTEXPR
  17. #error "The header <boost/math/norms.hpp> can only be used in C++17 and later."
  18. #endif
  19. #endif
  20. namespace boost::math::ccmath {
  21. namespace detail {
  22. template <typename T>
  23. constexpr T fma_imp(const T x, const T y, const T z) noexcept
  24. {
  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. #endif
  39. // If we can't use compiler intrinsics hope that -fma flag optimizes this call to fma instruction
  40. return (x * y) + z;
  41. }
  42. } // Namespace detail
  43. template <typename Real, std::enable_if_t<!std::is_integral_v<Real>, bool> = true>
  44. constexpr Real fma(Real x, Real y, Real z) noexcept
  45. {
  46. if (BOOST_MATH_IS_CONSTANT_EVALUATED(x))
  47. {
  48. if (x == 0 && boost::math::ccmath::isinf(y))
  49. {
  50. return std::numeric_limits<Real>::quiet_NaN();
  51. }
  52. else if (y == 0 && boost::math::ccmath::isinf(x))
  53. {
  54. return std::numeric_limits<Real>::quiet_NaN();
  55. }
  56. else if (boost::math::ccmath::isnan(x))
  57. {
  58. return std::numeric_limits<Real>::quiet_NaN();
  59. }
  60. else if (boost::math::ccmath::isnan(y))
  61. {
  62. return std::numeric_limits<Real>::quiet_NaN();
  63. }
  64. else if (boost::math::ccmath::isnan(z))
  65. {
  66. return std::numeric_limits<Real>::quiet_NaN();
  67. }
  68. return boost::math::ccmath::detail::fma_imp(x, y, z);
  69. }
  70. else
  71. {
  72. using std::fma;
  73. return fma(x, y, z);
  74. }
  75. }
  76. template <typename T1, typename T2, typename T3>
  77. constexpr auto fma(T1 x, T2 y, T3 z) noexcept
  78. {
  79. if (BOOST_MATH_IS_CONSTANT_EVALUATED(x))
  80. {
  81. // If the type is an integer (e.g. epsilon == 0) then set the epsilon value to 1 so that type is at a minimum
  82. // cast to double
  83. constexpr auto T1p = std::numeric_limits<T1>::epsilon() > 0 ? std::numeric_limits<T1>::epsilon() : 1;
  84. constexpr auto T2p = std::numeric_limits<T2>::epsilon() > 0 ? std::numeric_limits<T2>::epsilon() : 1;
  85. constexpr auto T3p = std::numeric_limits<T3>::epsilon() > 0 ? std::numeric_limits<T3>::epsilon() : 1;
  86. using promoted_type =
  87. #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
  88. std::conditional_t<T1p <= LDBL_EPSILON && T1p <= T2p, T1,
  89. std::conditional_t<T2p <= LDBL_EPSILON && T2p <= T1p, T2,
  90. std::conditional_t<T3p <= LDBL_EPSILON && T3p <= T2p, T3,
  91. #endif
  92. std::conditional_t<T1p <= DBL_EPSILON && T1p <= T2p, T1,
  93. std::conditional_t<T2p <= DBL_EPSILON && T2p <= T1p, T2,
  94. std::conditional_t<T3p <= DBL_EPSILON && T3p <= T2p, T3, double
  95. #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
  96. >>>>>>;
  97. #else
  98. >>>;
  99. #endif
  100. return boost::math::ccmath::fma(promoted_type(x), promoted_type(y), promoted_type(z));
  101. }
  102. else
  103. {
  104. using std::fma;
  105. return fma(x, y, z);
  106. }
  107. }
  108. constexpr float fmaf(float x, float y, float z) noexcept
  109. {
  110. return boost::math::ccmath::fma(x, y, z);
  111. }
  112. #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
  113. constexpr long double fmal(long double x, long double y, long double z) noexcept
  114. {
  115. return boost::math::ccmath::fma(x, y, z);
  116. }
  117. #endif
  118. } // Namespace boost::math::ccmath
  119. #endif // BOOST_MATH_CCMATH_FMA_HPP