powm1.hpp 3.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102
  1. // (C) Copyright John Maddock 2006.
  2. // (C) Copyright Matt Borland 2024.
  3. // Use, modification and distribution are subject to the
  4. // Boost Software License, Version 1.0. (See accompanying file
  5. // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  6. #ifndef BOOST_MATH_POWM1
  7. #define BOOST_MATH_POWM1
  8. #ifdef _MSC_VER
  9. #pragma once
  10. #pragma warning(push)
  11. #pragma warning(disable:4702) // Unreachable code (release mode only warning)
  12. #endif
  13. #include <boost/math/tools/config.hpp>
  14. #include <boost/math/special_functions/math_fwd.hpp>
  15. #include <boost/math/special_functions/log1p.hpp>
  16. #include <boost/math/special_functions/expm1.hpp>
  17. #include <boost/math/special_functions/trunc.hpp>
  18. #include <boost/math/special_functions/sign.hpp>
  19. #include <boost/math/tools/assert.hpp>
  20. namespace boost{ namespace math{ namespace detail{
  21. template <class T, class Policy>
  22. BOOST_MATH_GPU_ENABLED inline T powm1_imp(const T x, const T y, const Policy& pol)
  23. {
  24. BOOST_MATH_STD_USING
  25. constexpr auto function = "boost::math::powm1<%1%>(%1%, %1%)";
  26. if ((fabs(y * (x - 1)) < T(0.5)) || (fabs(y) < T(0.2)))
  27. {
  28. // We don't have any good/quick approximation for log(x) * y
  29. // so just try it and see:
  30. T l = y * log(x);
  31. if (l < T(0.5))
  32. return boost::math::expm1(l, pol);
  33. if (l > boost::math::tools::log_max_value<T>())
  34. return boost::math::policies::raise_overflow_error<T>(function, nullptr, pol);
  35. // fall through....
  36. }
  37. T result = pow(x, y) - 1;
  38. if((boost::math::isinf)(result))
  39. return result < 0 ? -boost::math::policies::raise_overflow_error<T>(function, nullptr, pol) : boost::math::policies::raise_overflow_error<T>(function, nullptr, pol);
  40. if((boost::math::isnan)(result))
  41. return boost::math::policies::raise_domain_error<T>(function, "Result of pow is complex or undefined", x, pol);
  42. return result;
  43. }
  44. template <class T, class Policy>
  45. BOOST_MATH_GPU_ENABLED inline T powm1_imp_dispatch(const T x, const T y, const Policy& pol)
  46. {
  47. BOOST_MATH_STD_USING
  48. if ((boost::math::signbit)(x)) // Need to error check -0 here as well
  49. {
  50. constexpr auto function = "boost::math::powm1<%1%>(%1%, %1%)";
  51. // y had better be an integer:
  52. if (boost::math::trunc(y) != y)
  53. return boost::math::policies::raise_domain_error<T>(function, "For non-integral exponent, expected base > 0 but got %1%", x, pol);
  54. if (boost::math::trunc(y / 2) == y / 2)
  55. return powm1_imp(T(-x), T(y), pol);
  56. }
  57. return powm1_imp(T(x), T(y), pol);
  58. }
  59. } // detail
  60. template <class T1, class T2>
  61. BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2>::type
  62. powm1(const T1 a, const T2 z)
  63. {
  64. typedef typename tools::promote_args<T1, T2>::type result_type;
  65. return detail::powm1_imp_dispatch(static_cast<result_type>(a), static_cast<result_type>(z), policies::policy<>());
  66. }
  67. template <class T1, class T2, class Policy>
  68. BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2>::type
  69. powm1(const T1 a, const T2 z, const Policy& pol)
  70. {
  71. typedef typename tools::promote_args<T1, T2>::type result_type;
  72. return detail::powm1_imp_dispatch(static_cast<result_type>(a), static_cast<result_type>(z), pol);
  73. }
  74. } // namespace math
  75. } // namespace boost
  76. #ifdef _MSC_VER
  77. #pragma warning(pop)
  78. #endif
  79. #endif // BOOST_MATH_POWM1