special_functions/detail/bessel_jy_derivatives_series.hpp

0001 //  Copyright (c) 2013 Anton Bikineev
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_JY_DERIVATIVES_SERIES_HPP
0007 #define BOOST_MATH_BESSEL_JY_DERIVATIVES_SERIES_HPP
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #endif
0012
0013 #include <cmath>
0014 #include <cstdint>
0015
0016 namespace boost{ namespace math{ namespace detail{
0017
0018 template <class T, class Policy>
0019 struct bessel_j_derivative_small_z_series_term
0020 {
0021    typedef T result_type;
0022
0023    bessel_j_derivative_small_z_series_term(T v_, T x)
0024       : N(0), v(v_), term(1), mult(x / 2)
0025    {
0026       mult *= -mult;
0027       // iterate if v == 0; otherwise result of
0028       // first term is 0 and tools::sum_series stops
0029       if (v == 0)
0030          iterate();
0031    }
0032    T operator()()
0033    {
0034       T r = term * (v + 2 * N);
0035       iterate();
0036       return r;
0037    }
0038 private:
0039    void iterate()
0040    {
0041       ++N;
0042       term *= mult / (N * (N + v));
0043    }
0044    unsigned N;
0045    T v;
0046    T term;
0047    T mult;
0048 };
0049 //
0050 // Series evaluation for BesselJ'(v, z) as z -> 0.
0051 // It's derivative of http://functions.wolfram.com/Bessel-TypeFunctions/BesselJ/06/01/04/01/01/0003/
0052 // Converges rapidly for all z << v.
0053 //
0054 template <class T, class Policy>
0055 inline T bessel_j_derivative_small_z_series(T v, T x, const Policy& pol)
0056 {
0057    BOOST_MATH_STD_USING
0058    T prefix;
0059    if (v < boost::math::max_factorial<T>::value)
0060    {
0061       prefix = pow(x / 2, v - 1) / 2 / boost::math::tgamma(v + 1, pol);
0062    }
0063    else
0064    {
0065       prefix = (v - 1) * log(x / 2) - constants::ln_two<T>() - boost::math::lgamma(v + 1, pol);
0066       prefix = exp(prefix);
0067    }
0068    if (0 == prefix)
0069       return prefix;
0070
0071    bessel_j_derivative_small_z_series_term<T, Policy> s(v, x);
0072    std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
0073
0074    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0075
0076    boost::math::policies::check_series_iterations<T>("boost::math::bessel_j_derivative_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
0077    return prefix * result;
0078 }
0079
0080 template <class T, class Policy>
0081 struct bessel_y_derivative_small_z_series_term_a
0082 {
0083    typedef T result_type;
0084
0085    bessel_y_derivative_small_z_series_term_a(T v_, T x)
0086       : N(0), v(v_)
0087    {
0088       mult = x / 2;
0089       mult *= -mult;
0090       term = 1;
0091    }
0092    T operator()()
0093    {
0094       T r = term * (-v + 2 * N);
0095       ++N;
0096       term *= mult / (N * (N - v));
0097       return r;
0098    }
0099 private:
0100    unsigned N;
0101    T v;
0102    T mult;
0103    T term;
0104 };
0105
0106 template <class T, class Policy>
0107 struct bessel_y_derivative_small_z_series_term_b
0108 {
0109    typedef T result_type;
0110
0111    bessel_y_derivative_small_z_series_term_b(T v_, T x)
0112       : N(0), v(v_)
0113    {
0114       mult = x / 2;
0115       mult *= -mult;
0116       term = 1;
0117    }
0118    T operator()()
0119    {
0120       T r = term * (v + 2 * N);
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 // It's derivative of 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_derivative_small_z_series(T v, T x, const Policy& pol)
0142 {
0143    BOOST_MATH_STD_USING
0144    static const char* function = "bessel_y_derivative_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 >= boost::math::max_factorial<T>::value) || (boost::math::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 + 1) * 2;
0155       if (boost::math::tools::max_value<T>() * p < gam)
0156       {
0157          scale /= gam;
0158          gam = 1;
0159          if (boost::math::tools::max_value<T>() * p < gam)
0160          {
0161             // This term will overflow to -INF, when combined with the series below it becomes +INF:
0162             return boost::math::policies::raise_overflow_error<T>(function, nullptr, pol);
0163          }
0164       }
0165       prefix = -gam / (boost::math::constants::pi<T>() * p);
0166    }
0167    else
0168    {
0169       gam = boost::math::lgamma(v, pol);
0170       p = (v + 1) * p + constants::ln_two<T>();
0171       prefix = gam - log(boost::math::constants::pi<T>()) - p;
0172       if (boost::math::tools::log_max_value<T>() < prefix)
0173       {
0174          prefix -= log(boost::math::tools::max_value<T>() / 4);
0175          scale /= (boost::math::tools::max_value<T>() / 4);
0176          if (boost::math::tools::log_max_value<T>() < prefix)
0177          {
0178             return boost::math::policies::raise_overflow_error<T>(function, nullptr, pol);
0179          }
0180       }
0181       prefix = -exp(prefix);
0182    }
0183    bessel_y_derivative_small_z_series_term_a<T, Policy> s(v, x);
0184    std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
0185
0186    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0187
0188    boost::math::policies::check_series_iterations<T>("boost::math::bessel_y_derivative_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
0189    result *= prefix;
0190
0191    p = pow(x / 2, v - 1) / 2;
0192    if (!need_logs)
0193    {
0194       prefix = boost::math::tgamma(-v, pol) * boost::math::cos_pi(v, pol) * p / boost::math::constants::pi<T>();
0195    }
0196    else
0197    {
0198       int sgn {};
0199       prefix = boost::math::lgamma(-v, &sgn, pol) + (v - 1) * log(x / 2) - constants::ln_two<T>();
0200       prefix = exp(prefix) * sgn / boost::math::constants::pi<T>();
0201    }
0202    bessel_y_derivative_small_z_series_term_b<T, Policy> s2(v, x);
0203    max_iter = boost::math::policies::get_max_series_iterations<Policy>();
0204
0205    T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0206
0207    result += scale * prefix * b;
0208    if(scale * tools::max_value<T>() < result)
0209       return boost::math::policies::raise_overflow_error<T>(function, nullptr, pol);
0210    return result / scale;
0211 }
0212
0213 // Calculating of BesselY'(v,x) with small x (x < epsilon) and integer x using derivatives
0214 // of formulas in http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/02/
0215 // seems to lose precision. Instead using linear combination of regular Bessel is preferred.
0216
0217 }}} // namespaces
0218
0219 #endif // BOOST_MATH_BESSEL_JY_DERIVATIVES_SERIES_HPP