special_functions/detail/bessel_jy_series.hpp

0001 //  Copyright (c) 2011 John Maddock
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_BESSEL_JN_SERIES_HPP
0007 #define BOOST_MATH_BESSEL_JN_SERIES_HPP
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #endif
0012
0013 #include <boost/math/tools/config.hpp>
0014 #include <boost/math/tools/assert.hpp>
0015 #include <boost/math/tools/cstdint.hpp>
0016
0017 namespace boost { namespace math { namespace detail{
0018
0019 template <class T, class Policy>
0020 struct bessel_j_small_z_series_term
0021 {
0022    typedef T result_type;
0023
0024    BOOST_MATH_GPU_ENABLED bessel_j_small_z_series_term(T v_, T x)
0025       : N(0), v(v_)
0026    {
0027       BOOST_MATH_STD_USING
0028       mult = x / 2;
0029       mult *= -mult;
0030       term = 1;
0031    }
0032    BOOST_MATH_GPU_ENABLED T operator()()
0033    {
0034       T r = term;
0035       ++N;
0036       term *= mult / (N * (N + v));
0037       return r;
0038    }
0039 private:
0040    unsigned N;
0041    T v;
0042    T mult;
0043    T term;
0044 };
0045 //
0046 // Series evaluation for BesselJ(v, z) as z -> 0.
0047 // See http://functions.wolfram.com/Bessel-TypeFunctions/BesselJ/06/01/04/01/01/0003/
0048 // Converges rapidly for all z << v.
0049 //
0050 template <class T, class Policy>
0051 BOOST_MATH_GPU_ENABLED inline T bessel_j_small_z_series(T v, T x, const Policy& pol)
0052 {
0053    BOOST_MATH_STD_USING
0054    T prefix;
0055    if(v < max_factorial<T>::value)
0056    {
0057       prefix = pow(x / 2, v) / boost::math::tgamma(v+1, pol);
0058    }
0059    else
0060    {
0061       prefix = v * log(x / 2) - boost::math::lgamma(v+1, pol);
0062       prefix = exp(prefix);
0063    }
0064    if(0 == prefix)
0065       return prefix;
0066
0067    bessel_j_small_z_series_term<T, Policy> s(v, x);
0068    boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0069
0070    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0071
0072    policies::check_series_iterations<T>("boost::math::bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
0073    return prefix * result;
0074 }
0075
0076 template <class T, class Policy>
0077 struct bessel_y_small_z_series_term_a
0078 {
0079    typedef T result_type;
0080
0081    BOOST_MATH_GPU_ENABLED bessel_y_small_z_series_term_a(T v_, T x)
0082       : N(0), v(v_)
0083    {
0084       BOOST_MATH_STD_USING
0085       mult = x / 2;
0086       mult *= -mult;
0087       term = 1;
0088    }
0089    BOOST_MATH_GPU_ENABLED T operator()()
0090    {
0091       BOOST_MATH_STD_USING
0092       T r = term;
0093       ++N;
0094       term *= mult / (N * (N - v));
0095       return r;
0096    }
0097 private:
0098    unsigned N;
0099    T v;
0100    T mult;
0101    T term;
0102 };
0103
0104 template <class T, class Policy>
0105 struct bessel_y_small_z_series_term_b
0106 {
0107    typedef T result_type;
0108
0109    BOOST_MATH_GPU_ENABLED bessel_y_small_z_series_term_b(T v_, T x)
0110       : N(0), v(v_)
0111    {
0112       BOOST_MATH_STD_USING
0113       mult = x / 2;
0114       mult *= -mult;
0115       term = 1;
0116    }
0117    BOOST_MATH_GPU_ENABLED T operator()()
0118    {
0119       T r = term;
0120       ++N;
0121       term *= mult / (N * (N + v));
0122       return r;
0123    }
0124 private:
0125    unsigned N;
0126    T v;
0127    T mult;
0128    T term;
0129 };
0130 //
0131 // Series form for BesselY as z -> 0,
0132 // see: http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/01/0003/
0133 // This series is only useful when the second term is small compared to the first
0134 // otherwise we get catastrophic cancellation errors.
0135 //
0136 // Approximating tgamma(v) by v^v, and assuming |tgamma(-z)| < eps we end up requiring:
0137 // eps/2 * v^v(x/2)^-v > (x/2)^v or log(eps/2) > v log((x/2)^2/v)
0138 //
0139 template <class T, class Policy>
0140 BOOST_MATH_GPU_ENABLED inline T bessel_y_small_z_series(T v, T x, T* pscale, const Policy& pol)
0141 {
0142    BOOST_MATH_STD_USING
0143    constexpr auto function = "bessel_y_small_z_series<%1%>(%1%,%1%)";
0144    T prefix;
0145    T gam;
0146    T p = log(x / 2);
0147    T scale = 1;
0148    bool need_logs = (v >= max_factorial<T>::value) || (tools::log_max_value<T>() / v < fabs(p));
0149
0150    if(!need_logs)
0151    {
0152       gam = boost::math::tgamma(v, pol);
0153       p = pow(x / 2, v);
0154       if(tools::max_value<T>() * p < gam)
0155       {
0156          scale /= gam;
0157          gam = 1;
0158          /*
0159          * We can never get here, it would require p < 1/max_value.
0160          if(tools::max_value<T>() * p < gam)
0161          {
0162             return -policies::raise_overflow_error<T>(function, nullptr, pol);
0163          }
0164          */
0165       }
0166       prefix = -gam / (constants::pi<T>() * p);
0167    }
0168    else
0169    {
0170       gam = boost::math::lgamma(v, pol);
0171       p = v * p;
0172       prefix = gam - log(constants::pi<T>()) - p;
0173       if(tools::log_max_value<T>() < prefix)
0174       {
0175          prefix -= log(tools::max_value<T>() / 4);
0176          scale /= (tools::max_value<T>() / 4);
0177          if(tools::log_max_value<T>() < prefix)
0178          {
0179             return -policies::raise_overflow_error<T>(function, nullptr, pol);
0180          }
0181       }
0182       prefix = -exp(prefix);
0183    }
0184    bessel_y_small_z_series_term_a<T, Policy> s(v, x);
0185    boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0186    *pscale = scale;
0187
0188    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0189
0190    policies::check_series_iterations<T>("boost::math::bessel_y_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
0191    result *= prefix;
0192
0193    if(!need_logs)
0194    {
0195       prefix = boost::math::tgamma(-v, pol) * boost::math::cos_pi(v, pol) * p / constants::pi<T>();
0196    }
0197    else
0198    {
0199       int sgn {};
0200       prefix = boost::math::lgamma(-v, &sgn, pol) + p;
0201       prefix = exp(prefix) * sgn / constants::pi<T>();
0202    }
0203    bessel_y_small_z_series_term_b<T, Policy> s2(v, x);
0204    max_iter = policies::get_max_series_iterations<Policy>();
0205
0206    T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0207
0208    result -= scale * prefix * b;
0209    return result;
0210 }
0211
0212 template <class T, class Policy>
0213 BOOST_MATH_GPU_ENABLED T bessel_yn_small_z(int n, T z, T* scale, const Policy& pol)
0214 {
0215    //
0216    // See http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/02/
0217    //
0218    // Note that when called we assume that x < epsilon and n is a positive integer.
0219    //
0220    BOOST_MATH_STD_USING
0221    BOOST_MATH_ASSERT(n >= 0);
0222    BOOST_MATH_ASSERT((z < policies::get_epsilon<T, Policy>()));
0223
0224    if(n == 0)
0225    {
0226       return (2 / constants::pi<T>()) * (log(z / 2) +  constants::euler<T>());
0227    }
0228    else if(n == 1)
0229    {
0230       return (z / constants::pi<T>()) * log(z / 2)
0231          - 2 / (constants::pi<T>() * z)
0232          - (z / (2 * constants::pi<T>())) * (1 - 2 * constants::euler<T>());
0233    }
0234    else if(n == 2)
0235    {
0236       return (z * z) / (4 * constants::pi<T>()) * log(z / 2)
0237          - (4 / (constants::pi<T>() * z * z))
0238          - ((z * z) / (8 * constants::pi<T>())) * (T(3)/2 - 2 * constants::euler<T>());
0239    }
0240    else
0241    {
0242       #if (defined(__GNUC__) && __GNUC__ == 13)
0243       auto p = static_cast<T>(pow(z / 2, T(n)));
0244       #else
0245       auto p = static_cast<T>(pow(z / 2, n));
0246       #endif
0247
0248       T result = -((boost::math::factorial<T>(n - 1, pol) / constants::pi<T>()));
0249       if(p * tools::max_value<T>() < fabs(result))
0250       {
0251          T div = tools::max_value<T>() / 8;
0252          result /= div;
0253          *scale /= div;
0254          if(p * tools::max_value<T>() < result)
0255          {
0256             // Impossible to get here??
0257             return -policies::raise_overflow_error<T>("bessel_yn_small_z<%1%>(%1%,%1%)", nullptr, pol); // LCOV_EXCL_LINE
0258          }
0259       }
0260       return result / p;
0261    }
0262 }
0263
0264 }}} // namespaces
0265
0266 #endif // BOOST_MATH_BESSEL_JN_SERIES_HPP
0267