X-Git-Url: https://git.donarmstrong.com/?p=rsem.git;a=blobdiff_plain;f=boost%2Fmath%2Fspecial_functions%2Fdetail%2Fbessel_jy_asym.hpp;fp=boost%2Fmath%2Fspecial_functions%2Fdetail%2Fbessel_jy_asym.hpp;h=81f6238e586eb428d378574c42905147a6e4b795;hp=0000000000000000000000000000000000000000;hb=2d71eb92104693ca9baa5a2e1c23eeca776d8fd3;hpb=da57529b92adbb7ae74a89861cb39fb35ac7c62d diff --git a/boost/math/special_functions/detail/bessel_jy_asym.hpp b/boost/math/special_functions/detail/bessel_jy_asym.hpp new file mode 100644 index 0000000..81f6238 --- /dev/null +++ b/boost/math/special_functions/detail/bessel_jy_asym.hpp @@ -0,0 +1,205 @@ +// 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 + +namespace boost{ namespace math{ namespace detail{ + +template +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() * x)); +} + +template +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 - PI(v/2 + 1/4)) + // which we'll factor in later when we calculate the sines/cosines of the result: + // + T mu = 4 * v * v; + T denom = 4 * x; + T denom_mult = denom * denom; + + T s = 0; + 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 +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); + BOOST_MATH_INSTRUMENT_VARIABLE(ampl); + BOOST_MATH_INSTRUMENT_VARIABLE(phase); + // + // Calculate the sine of the phase, using + // sine/cosine addition rules to factor in + // the x - PI(v/2 + 1/4) term not added to the + // phase when we calculated it. + // + T cx = cos(x); + T sx = sin(x); + T ci = cos_pi(v / 2 + 0.25f); + T si = sin_pi(v / 2 + 0.25f); + T sin_phase = sin(phase) * (cx * ci + sx * si) + cos(phase) * (sx * ci - cx * si); + BOOST_MATH_INSTRUMENT_CODE(sin(phase)); + BOOST_MATH_INSTRUMENT_CODE(cos(x)); + BOOST_MATH_INSTRUMENT_CODE(cos(phase)); + BOOST_MATH_INSTRUMENT_CODE(sin(x)); + return sin_phase * ampl; +} + +template +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); + BOOST_MATH_INSTRUMENT_VARIABLE(ampl); + BOOST_MATH_INSTRUMENT_VARIABLE(phase); + // + // Calculate the sine of the phase, using + // sine/cosine addition rules to factor in + // the x - PI(v/2 + 1/4) term not added to the + // phase when we calculated it. + // + BOOST_MATH_INSTRUMENT_CODE(cos(phase)); + BOOST_MATH_INSTRUMENT_CODE(cos(x)); + BOOST_MATH_INSTRUMENT_CODE(sin(phase)); + BOOST_MATH_INSTRUMENT_CODE(sin(x)); + T cx = cos(x); + T sx = sin(x); + T ci = cos_pi(v / 2 + 0.25f); + T si = sin_pi(v / 2 + 0.25f); + T sin_phase = cos(phase) * (cx * ci + sx * si) - sin(phase) * (sx * ci - cx * si); + BOOST_MATH_INSTRUMENT_VARIABLE(sin_phase); + return sin_phase * ampl; +} + +template +inline bool asymptotic_bessel_large_x_limit(const T& v, const T& x) +{ + BOOST_MATH_STD_USING + // + // Determines if x is large enough compared to v to take the asymptotic + // forms above. From A&S 9.2.28 we require: + // v < x * eps^1/8 + // and from A&S 9.2.29 we require: + // v^12/10 < 1.5 * x * eps^1/10 + // using the former seems to work OK in practice with broadly similar + // error rates either side of the divide for v < 10000. + // At double precision eps^1/8 ~= 0.01. + // + return (std::max)(T(fabs(v)), T(1)) < x * sqrt(tools::forth_root_epsilon()); +} + +template +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(); ++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(k, pol); + g = f + g_prefix * q; + h = -k * g + p; + y += c * g; + y1 += c * h; + if(c * g / tools::epsilon() < y) + break; + } + + *Y = -y; + *Y1 = (-2 / x) * y1; +} + +template +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())); + + return (boost::math::isfinite)(s) ? + s : policies::raise_overflow_error("boost::math::asymptotic_bessel_i_large_x<%1%>(%1%,%1%)", 0, pol); +} + +}}} // namespaces + +#endif +