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