187 lines · plain
1// Copyright John Maddock 2012.2// Copyright Matt Borland 2024.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0.5// (See accompanying file LICENSE_1_0.txt6// or copy at http://www.boost.org/LICENSE_1_0.txt)7 8#ifndef BOOST_MATH_HANKEL_HPP9#define BOOST_MATH_HANKEL_HPP10 11#include <boost/math/tools/config.hpp>12#include <boost/math/tools/complex.hpp>13#include <boost/math/special_functions/math_fwd.hpp>14#include <boost/math/special_functions/bessel.hpp>15#include <boost/math/special_functions/detail/iconv.hpp>16#include <boost/math/constants/constants.hpp>17#include <boost/math/policies/error_handling.hpp>18 19namespace boost{ namespace math{20 21namespace detail{22 23template <class T, class Policy>24BOOST_MATH_GPU_ENABLED boost::math::complex<T> hankel_imp(T v, T x, const bessel_no_int_tag&, const Policy& pol, int sign)25{26 BOOST_MATH_STD_USING27 constexpr auto function = "boost::math::cyl_hankel_1<%1%>(%1%,%1%)";28 29 if(x < 0)30 {31 bool isint_v = floor(v) == v;32 T j, y;33 bessel_jy(v, -x, &j, &y, need_j | need_y, pol);34 boost::math::complex<T> cx(x), cv(v);35 boost::math::complex<T> j_result, y_result;36 if(isint_v)37 {38 int s = (iround(v) & 1) ? -1 : 1;39 j_result = j * s;40 y_result = T(s) * (y - (2 / constants::pi<T>()) * (log(-x) - log(cx)) * j);41 }42 else43 {44 j_result = pow(cx, v) * pow(-cx, -v) * j;45 T p1 = pow(-x, v);46 boost::math::complex<T> p2 = pow(cx, v);47 y_result = p1 * y / p248 + (p2 / p1 - p1 / p2) * j / tan(constants::pi<T>() * v);49 }50 // multiply y_result by i:51 y_result = boost::math::complex<T>(-sign * y_result.imag(), sign * y_result.real());52 return j_result + y_result;53 }54 55 if(x == 0)56 {57 if(v == 0)58 {59 // J is 1, Y is -INF60 return boost::math::complex<T>(1, sign * -policies::raise_overflow_error<T>(function, nullptr, pol));61 }62 else63 {64 // At least one of J and Y is complex infinity:65 return boost::math::complex<T>(policies::raise_overflow_error<T>(function, nullptr, pol), sign * policies::raise_overflow_error<T>(function, nullptr, pol));66 }67 }68 69 T j, y;70 bessel_jy(v, x, &j, &y, need_j | need_y, pol);71 return boost::math::complex<T>(j, sign * y);72}73 74template <class T, class Policy>75BOOST_MATH_GPU_ENABLED boost::math::complex<T> hankel_imp(int v, T x, const bessel_int_tag&, const Policy& pol, int sign);76 77template <class T, class Policy>78BOOST_MATH_GPU_ENABLED inline boost::math::complex<T> hankel_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol, int sign)79{80 BOOST_MATH_STD_USING // ADL of std names.81 int ival = detail::iconv(v, pol);82 if(0 == v - ival)83 {84 return hankel_imp(ival, x, bessel_int_tag(), pol, sign);85 }86 return hankel_imp(v, x, bessel_no_int_tag(), pol, sign);87}88 89template <class T, class Policy>90BOOST_MATH_GPU_ENABLED inline boost::math::complex<T> hankel_imp(int v, T x, const bessel_int_tag&, const Policy& pol, int sign)91{92 BOOST_MATH_STD_USING93 if((abs(v) < 200) && (x > 0))94 return boost::math::complex<T>(bessel_jn(v, x, pol), sign * bessel_yn(v, x, pol));95 return hankel_imp(static_cast<T>(v), x, bessel_no_int_tag(), pol, sign);96}97 98template <class T, class Policy>99BOOST_MATH_GPU_ENABLED inline boost::math::complex<T> sph_hankel_imp(T v, T x, const Policy& pol, int sign)100{101 BOOST_MATH_STD_USING102 return constants::root_half_pi<T>() * hankel_imp(v + 0.5f, x, bessel_no_int_tag(), pol, sign) / sqrt(boost::math::complex<T>(x));103}104 105} // namespace detail106 107template <class T1, class T2, class Policy>108BOOST_MATH_GPU_ENABLED inline boost::math::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> cyl_hankel_1(T1 v, T2 x, const Policy& pol)109{110 BOOST_FPU_EXCEPTION_GUARD111 typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;112 typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;113 typedef typename policies::evaluation<result_type, Policy>::type value_type;114 return policies::checked_narrowing_cast<boost::math::complex<result_type>, Policy>(detail::hankel_imp<value_type>(v, static_cast<value_type>(x), tag_type(), pol, 1), "boost::math::cyl_hankel_1<%1%>(%1%,%1%)");115}116 117template <class T1, class T2>118BOOST_MATH_GPU_ENABLED inline boost::math::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> cyl_hankel_1(T1 v, T2 x)119{120 return cyl_hankel_1(v, x, policies::policy<>());121}122 123template <class T1, class T2, class Policy>124BOOST_MATH_GPU_ENABLED inline boost::math::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> cyl_hankel_2(T1 v, T2 x, const Policy& pol)125{126 BOOST_FPU_EXCEPTION_GUARD127 typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;128 typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;129 typedef typename policies::evaluation<result_type, Policy>::type value_type;130 return policies::checked_narrowing_cast<boost::math::complex<result_type>, Policy>(detail::hankel_imp<value_type>(v, static_cast<value_type>(x), tag_type(), pol, -1), "boost::math::cyl_hankel_1<%1%>(%1%,%1%)");131}132 133template <class T1, class T2>134BOOST_MATH_GPU_ENABLED inline boost::math::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> cyl_hankel_2(T1 v, T2 x)135{136 return cyl_hankel_2(v, x, policies::policy<>());137}138 139template <class T1, class T2, class Policy>140BOOST_MATH_GPU_ENABLED inline boost::math::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> sph_hankel_1(T1 v, T2 x, const Policy&)141{142 BOOST_FPU_EXCEPTION_GUARD143 typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;144 typedef typename policies::evaluation<result_type, Policy>::type value_type;145 typedef typename policies::normalise<146 Policy,147 policies::promote_float<false>,148 policies::promote_double<false>,149 policies::discrete_quantile<>,150 policies::assert_undefined<> >::type forwarding_policy;151 152 return policies::checked_narrowing_cast<boost::math::complex<result_type>, Policy>(detail::sph_hankel_imp<value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy(), 1), "boost::math::sph_hankel_1<%1%>(%1%,%1%)");153}154 155template <class T1, class T2>156BOOST_MATH_GPU_ENABLED inline boost::math::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> sph_hankel_1(T1 v, T2 x)157{158 return sph_hankel_1(v, x, policies::policy<>());159}160 161template <class T1, class T2, class Policy>162BOOST_MATH_GPU_ENABLED inline boost::math::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> sph_hankel_2(T1 v, T2 x, const Policy&)163{164 BOOST_FPU_EXCEPTION_GUARD165 typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;166 typedef typename policies::evaluation<result_type, Policy>::type value_type;167 typedef typename policies::normalise<168 Policy,169 policies::promote_float<false>,170 policies::promote_double<false>,171 policies::discrete_quantile<>,172 policies::assert_undefined<> >::type forwarding_policy;173 174 return policies::checked_narrowing_cast<boost::math::complex<result_type>, Policy>(detail::sph_hankel_imp<value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy(), -1), "boost::math::sph_hankel_1<%1%>(%1%,%1%)");175}176 177template <class T1, class T2>178BOOST_MATH_GPU_ENABLED inline boost::math::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> sph_hankel_2(T1 v, T2 x)179{180 return sph_hankel_2(v, x, policies::policy<>());181}182 183}} // namespaces184 185#endif // BOOST_MATH_HANKEL_HPP186 187