--- /dev/null
+// 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_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 - 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 <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);
+ 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 <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);
+ 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 <class T>
+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<T>());
+}
+
+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
+
projects / rsem.git / blobdiff