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    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          if(tools::max_value<T>() * p < gam)
0159          {
0160             return -policies::raise_overflow_error<T>(function, nullptr, pol);
0161          }
0162       }
0163       prefix = -gam / (constants::pi<T>() * p);
0164    }
0165    else
0166    {
0167       gam = boost::math::lgamma(v, pol);
0168       p = v * p;
0169       prefix = gam - log(constants::pi<T>()) - p;
0170       if(tools::log_max_value<T>() < prefix)
0171       {
0172          prefix -= log(tools::max_value<T>() / 4);
0173          scale /= (tools::max_value<T>() / 4);
0174          if(tools::log_max_value<T>() < prefix)
0175          {
0176             return -policies::raise_overflow_error<T>(function, nullptr, pol);
0177          }
0178       }
0179       prefix = -exp(prefix);
0180    }
0181    bessel_y_small_z_series_term_a<T, Policy> s(v, x);
0182    std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0183    *pscale = scale;
0184
0185    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0186
0187    policies::check_series_iterations<T>("boost::math::bessel_y_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
0188    result *= prefix;
0189
0190    if(!need_logs)
0191    {
0192       prefix = boost::math::tgamma(-v, pol) * boost::math::cos_pi(v, pol) * p / constants::pi<T>();
0193    }
0194    else
0195    {
0196       int sgn {};
0197       prefix = boost::math::lgamma(-v, &sgn, pol) + p;
0198       prefix = exp(prefix) * sgn / constants::pi<T>();
0199    }
0200    bessel_y_small_z_series_term_b<T, Policy> s2(v, x);
0201    max_iter = policies::get_max_series_iterations<Policy>();
0202
0203    T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0204
0205    result -= scale * prefix * b;
0206    return result;
0207 }
0208
0209 template <class T, class Policy>
0210 T bessel_yn_small_z(int n, T z, T* scale, const Policy& pol)
0211 {
0212    //
0213    // See http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/02/
0214    //
0215    // Note that when called we assume that x < epsilon and n is a positive integer.
0216    //
0217    BOOST_MATH_STD_USING
0218    BOOST_MATH_ASSERT(n >= 0);
0219    BOOST_MATH_ASSERT((z < policies::get_epsilon<T, Policy>()));
0220
0221    if(n == 0)
0222    {
0223       return (2 / constants::pi<T>()) * (log(z / 2) +  constants::euler<T>());
0224    }
0225    else if(n == 1)
0226    {
0227       return (z / constants::pi<T>()) * log(z / 2)
0228          - 2 / (constants::pi<T>() * z)
0229          - (z / (2 * constants::pi<T>())) * (1 - 2 * constants::euler<T>());
0230    }
0231    else if(n == 2)
0232    {
0233       return (z * z) / (4 * constants::pi<T>()) * log(z / 2)
0234          - (4 / (constants::pi<T>() * z * z))
0235          - ((z * z) / (8 * constants::pi<T>())) * (T(3)/2 - 2 * constants::euler<T>());
0236    }
0237    else
0238    {
0239       #if (defined(__GNUC__) && __GNUC__ == 13)
0240       auto p = static_cast<T>(pow(z / 2, T(n)));
0241       #else
0242       auto p = static_cast<T>(pow(z / 2, n));
0243       #endif
0244
0245       T result = -((boost::math::factorial<T>(n - 1, pol) / constants::pi<T>()));
0246       if(p * tools::max_value<T>() < result)
0247       {
0248          T div = tools::max_value<T>() / 8;
0249          result /= div;
0250          *scale /= div;
0251          if(p * tools::max_value<T>() < result)
0252          {
0253             return -policies::raise_overflow_error<T>("bessel_yn_small_z<%1%>(%1%,%1%)", nullptr, pol);
0254          }
0255       }
0256       return result / p;
0257    }
0258 }
0259
0260 }}} // namespaces
0261
0262 #endif // BOOST_MATH_BESSEL_JN_SERIES_HPP
0263