bessel_jy_asym.hpp

来自「Boost provides free peer-reviewed portab」· HPP 代码 · 共 303 行

//  Copyright (c) 2007 John Maddock
//  Use, modification and distribution are subject to the
//  Boost Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

//
// This is a partial header, do not include on it's own!!!
//
// Contains asymptotic expansions for Bessel J(v,x) and Y(v,x)
// functions, as x -> INF.
//
#ifndef BOOST_MATH_SF_DETAIL_BESSEL_JY_ASYM_HPP
#define BOOST_MATH_SF_DETAIL_BESSEL_JY_ASYM_HPP

#ifdef _MSC_VER
#pragma once
#endif

#include <boost/math/special_functions/factorials.hpp>

namespace boost{ namespace math{ namespace detail{

template <class T>
inline T asymptotic_bessel_j_large_x_P(T v, T x)
{
   // A&S 9.2.9
   T s = 1;
   T mu = 4 * v * v;
   T ez2 = 8 * x;
   ez2 *= ez2;
   s -= (mu-1) * (mu-9) / (2 * ez2);
   s += (mu-1) * (mu-9) * (mu-25) * (mu - 49) / (24 * ez2 * ez2);
   return s;
}

template <class T>
inline T asymptotic_bessel_j_large_x_Q(T v, T x)
{
   // A&S 9.2.10
   T s = 0;
   T mu = 4 * v * v;
   T ez = 8*x;
   s += (mu-1) / ez;
   s -= (mu-1) * (mu-9) * (mu-25) / (6 * ez*ez*ez);
   return s;
}

template <class T>
inline T asymptotic_bessel_j_large_x(T v, T x)
{
   // 
   // See http://functions.wolfram.com/BesselAiryStruveFunctions/BesselJ/06/02/02/0001/
   //
   // Also A&S 9.2.5
   //
   BOOST_MATH_STD_USING // ADL of std names
   T chi = fabs(x) - constants::pi<T>() * (2 * v + 1) / 4;
   return sqrt(2 / (constants::pi<T>() * x))
      * (asymptotic_bessel_j_large_x_P(v, x) * cos(chi) 
         - asymptotic_bessel_j_large_x_Q(v, x) * sin(chi));
}

template <class T>
inline T asymptotic_bessel_y_large_x(T v, T x)
{
   // 
   // See http://functions.wolfram.com/BesselAiryStruveFunctions/BesselJ/06/02/02/0001/
   //
   // Also A&S 9.2.5
   //
   BOOST_MATH_STD_USING // ADL of std names
   T chi = fabs(x) - constants::pi<T>() * (2 * v + 1) / 4;
   return sqrt(2 / (constants::pi<T>() * x))
      * (asymptotic_bessel_j_large_x_P(v, x) * sin(chi) 
         - asymptotic_bessel_j_large_x_Q(v, x) * cos(chi));
}

template <class T>
inline T asymptotic_bessel_amplitude(T v, T x)
{
   // Calculate the amplitude of J(v, x) and Y(v, x) for large
   // x: see A&S 9.2.28.
   BOOST_MATH_STD_USING
   T s = 1;
   T mu = 4 * v * v;
   T txq = 2 * x;
   txq *= txq;

   s += (mu - 1) / (2 * txq);
   s += 3 * (mu - 1) * (mu - 9) / (txq * txq * 8);
   s += 15 * (mu - 1) * (mu - 9) * (mu - 25) / (txq * txq * txq * 8 * 6);

   return sqrt(s * 2 / (constants::pi<T>() * x));
}

template <class T>
T asymptotic_bessel_phase_mx(T v, T x)
{
   //
   // Calculate the phase of J(v, x) and Y(v, x) for large x.
   // See A&S 9.2.29.
   // Note that the result returned is the phase less x.
   //
   T mu = 4 * v * v;
   T denom = 4 * x;
   T denom_mult = denom * denom;

   T s = -constants::pi<T>() * (v / 2 + 0.25f);
   s += (mu - 1) / (2 * denom);
   denom *= denom_mult;
   s += (mu - 1) * (mu - 25) / (6 * denom);
   denom *= denom_mult;
   s += (mu - 1) * (mu * mu - 114 * mu + 1073) / (5 * denom);
   denom *= denom_mult;
   s += (mu - 1) * (5 * mu * mu * mu - 1535 * mu * mu + 54703 * mu - 375733) / (14 * denom);
   return s;
}

template <class T>
inline T asymptotic_bessel_y_large_x_2(T v, T x)
{
   // See A&S 9.2.19.
   BOOST_MATH_STD_USING
   // Get the phase and amplitude:
   T ampl = asymptotic_bessel_amplitude(v, x);
   T phase = asymptotic_bessel_phase_mx(v, x);
   //
   // Calculate the sine of the phase, using:
   // sin(x+p) = sin(x)cos(p) + cos(x)sin(p)
   //
   T sin_phase = sin(phase) * cos(x) + cos(phase) * sin(x);
   return sin_phase * ampl;
}

template <class T>
inline T asymptotic_bessel_j_large_x_2(T v, T x)
{
   // See A&S 9.2.19.
   BOOST_MATH_STD_USING
   // Get the phase and amplitude:
   T ampl = asymptotic_bessel_amplitude(v, x);
   T phase = asymptotic_bessel_phase_mx(v, x);
   //
   // Calculate the sine of the phase, using:
   // cos(x+p) = cos(x)cos(p) - sin(x)sin(p)
   //
   T sin_phase = cos(phase) * cos(x) - sin(phase) * sin(x);
   return sin_phase * ampl;
}

//
// Various limits for the J and Y asymptotics
// (the asympotic expansions are safe to use if
// x is less than the limit given).
// We assume that if we don't use these expansions then the
// error will likely be >100eps, so the limits given are chosen
// to lead to < 100eps truncation error.
//
template <class T>
inline T asymptotic_bessel_y_limit(const mpl::int_<0>&)
{
   // default case:
   BOOST_MATH_STD_USING
   return 2.25 / pow(100 * tools::epsilon<T>() / T(0.001f), T(0.2f));
}
template <class T>
inline T asymptotic_bessel_y_limit(const mpl::int_<53>&)
{
   // double case:
   return 304 /*780*/;
}
template <class T>
inline T asymptotic_bessel_y_limit(const mpl::int_<64>&)
{
   // 80-bit extended-double case:
   return 1552 /*3500*/;
}
template <class T>
inline T asymptotic_bessel_y_limit(const mpl::int_<113>&)
{
   // 128-bit long double case:
   return 1245243 /*3128000*/;
}

template <class T, class Policy>
struct bessel_asymptotic_tag
{
   typedef typename policies::precision<T, Policy>::type precision_type;
   typedef typename mpl::if_<
      mpl::or_<
         mpl::equal_to<precision_type, mpl::int_<0> >,
         mpl::greater<precision_type, mpl::int_<113> > >,
      mpl::int_<0>,
      typename mpl::if_<
         mpl::greater<precision_type, mpl::int_<64> >,
         mpl::int_<113>,
         typename mpl::if_<
            mpl::greater<precision_type, mpl::int_<53> >,
            mpl::int_<64>,
            mpl::int_<53>
         >::type
      >::type
   >::type type;
};

template <class T>
inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<0>&)
{
   // default case:
   BOOST_MATH_STD_USING
   T v2 = (std::max)(T(3), v * v);
   return v2 / pow(100 * tools::epsilon<T>() / T(2e-5f), T(0.17f));
}
template <class T>
inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<53>&)
{
   // double case:
   T v2 = (std::max)(T(3), v * v);
   return v2 * 33 /*73*/;
}
template <class T>
inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<64>&)
{
   // 80-bit extended-double case:
   T v2 = (std::max)(T(3), v * v);
   return v2 * 121 /*266*/;
}
template <class T>
inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<113>&)
{
   // 128-bit long double case:
   T v2 = (std::max)(T(3), v * v);
   return v2 * 39154 /*85700*/;
}

template <class T, class Policy>
void temme_asyptotic_y_small_x(T v, T x, T* Y, T* Y1, const Policy& pol)
{
   T c = 1;
   T p = (v / boost::math::sin_pi(v, pol)) * pow(x / 2, -v) / boost::math::tgamma(1 - v, pol);
   T q = (v / boost::math::sin_pi(v, pol)) * pow(x / 2, v) / boost::math::tgamma(1 + v, pol);
   T f = (p - q) / v;
   T g_prefix = boost::math::sin_pi(v / 2, pol);
   g_prefix *= g_prefix * 2 / v;
   T g = f + g_prefix * q;
   T h = p;
   T c_mult = -x * x / 4;

   T y(c * g), y1(c * h);

   for(int k = 1; k < policies::get_max_series_iterations<Policy>(); ++k)
   {
      f = (k * f + p + q) / (k*k - v*v);
      p /= k - v;
      q /= k + v;
      c *= c_mult / k;
      T c1 = pow(-x * x / 4, k) / factorial<T>(k, pol);
      g = f + g_prefix * q;
      h = -k * g + p;
      y += c * g;
      y1 += c * h;
      if(c * g / tools::epsilon<T>() < y)
         break;
   }

   *Y = -y;
   *Y1 = (-2 / x) * y1;
}

template <class T, class Policy>
T asymptotic_bessel_i_large_x(T v, T x, const Policy& pol)
{
   BOOST_MATH_STD_USING  // ADL of std names
   T s = 1;
   T mu = 4 * v * v;
   T ex = 8 * x;
   T num = mu - 1;
   T denom = ex;

   s -= num / denom;

   num *= mu - 9;
   denom *= ex * 2;
   s += num / denom;

   num *= mu - 25;
   denom *= ex * 3;
   s -= num / denom;

   // Try and avoid overflow to the last minute:
   T e = exp(x/2);

   s = e * (e * s / sqrt(2 * x * constants::pi<T>()));

   return (boost::math::isfinite)(s) ? 
      s : policies::raise_overflow_error<T>("boost::math::asymptotic_bessel_i_large_x<%1%>(%1%,%1%)", 0, pol);
}

}}} // namespaces

#endif