1// Copyright (c) 2007 John Maddock
2// Use, modification and distribution are subject to the
3// Boost Software License, Version 1.0. (See accompanying file
4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5
6//
7// This is a partial header, do not include on it's own!!!
8//
9// Contains asymptotic expansions for Bessel J(v,x) and Y(v,x)
10// functions, as x -> INF.
11//
12#ifndef BOOST_MATH_SF_DETAIL_BESSEL_JY_ASYM_HPP
13#define BOOST_MATH_SF_DETAIL_BESSEL_JY_ASYM_HPP
14
15#ifdef _MSC_VER
16#pragma once
17#endif
18
19#include <boost/math/special_functions/factorials.hpp>
20
21namespace boost{ namespace math{ namespace detail{
22
23template <class T>
24inline T asymptotic_bessel_amplitude(T v, T x)
25{
26 // Calculate the amplitude of J(v, x) and Y(v, x) for large
27 // x: see A&S 9.2.28.
28 BOOST_MATH_STD_USING
29 T s = 1;
30 T mu = 4 * v * v;
31 T txq = 2 * x;
32 txq *= txq;
33
34 s += (mu - 1) / (2 * txq);
35 s += 3 * (mu - 1) * (mu - 9) / (txq * txq * 8);
36 s += 15 * (mu - 1) * (mu - 9) * (mu - 25) / (txq * txq * txq * 8 * 6);
37
38 return sqrt(s * 2 / (constants::pi<T>() * x));
39}
40
41template <class T>
42T asymptotic_bessel_phase_mx(T v, T x)
43{
44 //
45 // Calculate the phase of J(v, x) and Y(v, x) for large x.
46 // See A&S 9.2.29.
47 // Note that the result returned is the phase less (x - PI(v/2 + 1/4))
48 // which we'll factor in later when we calculate the sines/cosines of the result:
49 //
50 T mu = 4 * v * v;
51 T denom = 4 * x;
52 T denom_mult = denom * denom;
53
54 T s = 0;
55 s += (mu - 1) / (2 * denom);
56 denom *= denom_mult;
57 s += (mu - 1) * (mu - 25) / (6 * denom);
58 denom *= denom_mult;
59 s += (mu - 1) * (mu * mu - 114 * mu + 1073) / (5 * denom);
60 denom *= denom_mult;
61 s += (mu - 1) * (5 * mu * mu * mu - 1535 * mu * mu + 54703 * mu - 375733) / (14 * denom);
62 return s;
63}
64
65template <class T>
66inline T asymptotic_bessel_y_large_x_2(T v, T x)
67{
68 // See A&S 9.2.19.
69 BOOST_MATH_STD_USING
70 // Get the phase and amplitude:
71 T ampl = asymptotic_bessel_amplitude(v, x);
72 T phase = asymptotic_bessel_phase_mx(v, x);
73 BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
74 BOOST_MATH_INSTRUMENT_VARIABLE(phase);
75 //
76 // Calculate the sine of the phase, using
77 // sine/cosine addition rules to factor in
78 // the x - PI(v/2 + 1/4) term not added to the
79 // phase when we calculated it.
80 //
81 T cx = cos(x);
82 T sx = sin(x);
83 T ci = cos_pi(v / 2 + 0.25f);
84 T si = sin_pi(v / 2 + 0.25f);
85 T sin_phase = sin(phase) * (cx * ci + sx * si) + cos(phase) * (sx * ci - cx * si);
86 BOOST_MATH_INSTRUMENT_CODE(sin(phase));
87 BOOST_MATH_INSTRUMENT_CODE(cos(x));
88 BOOST_MATH_INSTRUMENT_CODE(cos(phase));
89 BOOST_MATH_INSTRUMENT_CODE(sin(x));
90 return sin_phase * ampl;
91}
92
93template <class T>
94inline T asymptotic_bessel_j_large_x_2(T v, T x)
95{
96 // See A&S 9.2.19.
97 BOOST_MATH_STD_USING
98 // Get the phase and amplitude:
99 T ampl = asymptotic_bessel_amplitude(v, x);
100 T phase = asymptotic_bessel_phase_mx(v, x);
101 BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
102 BOOST_MATH_INSTRUMENT_VARIABLE(phase);
103 //
104 // Calculate the sine of the phase, using
105 // sine/cosine addition rules to factor in
106 // the x - PI(v/2 + 1/4) term not added to the
107 // phase when we calculated it.
108 //
109 BOOST_MATH_INSTRUMENT_CODE(cos(phase));
110 BOOST_MATH_INSTRUMENT_CODE(cos(x));
111 BOOST_MATH_INSTRUMENT_CODE(sin(phase));
112 BOOST_MATH_INSTRUMENT_CODE(sin(x));
113 T cx = cos(x);
114 T sx = sin(x);
115 T ci = cos_pi(v / 2 + 0.25f);
116 T si = sin_pi(v / 2 + 0.25f);
117 T sin_phase = cos(phase) * (cx * ci + sx * si) - sin(phase) * (sx * ci - cx * si);
118 BOOST_MATH_INSTRUMENT_VARIABLE(sin_phase);
119 return sin_phase * ampl;
120}
121
122template <class T>
123inline bool asymptotic_bessel_large_x_limit(int v, const T& x)
124{
125 BOOST_MATH_STD_USING
126 //
127 // Determines if x is large enough compared to v to take the asymptotic
128 // forms above. From A&S 9.2.28 we require:
129 // v < x * eps^1/8
130 // and from A&S 9.2.29 we require:
131 // v^12/10 < 1.5 * x * eps^1/10
132 // using the former seems to work OK in practice with broadly similar
133 // error rates either side of the divide for v < 10000.
134 // At double precision eps^1/8 ~= 0.01.
135 //
136 BOOST_ASSERT(v >= 0);
137 return (v ? v : 1) < x * 0.004f;
138}
139
140template <class T>
141inline bool asymptotic_bessel_large_x_limit(const T& v, const T& x)
142{
143 BOOST_MATH_STD_USING
144 //
145 // Determines if x is large enough compared to v to take the asymptotic
146 // forms above. From A&S 9.2.28 we require:
147 // v < x * eps^1/8
148 // and from A&S 9.2.29 we require:
149 // v^12/10 < 1.5 * x * eps^1/10
150 // using the former seems to work OK in practice with broadly similar
151 // error rates either side of the divide for v < 10000.
152 // At double precision eps^1/8 ~= 0.01.
153 //
154 return (std::max)(T(fabs(v)), T(1)) < x * sqrt(tools::forth_root_epsilon<T>());
155}
156
157template <class T, class Policy>
158void temme_asyptotic_y_small_x(T v, T x, T* Y, T* Y1, const Policy& pol)
159{
160 T c = 1;
161 T p = (v / boost::math::sin_pi(v, pol)) * pow(x / 2, -v) / boost::math::tgamma(1 - v, pol);
162 T q = (v / boost::math::sin_pi(v, pol)) * pow(x / 2, v) / boost::math::tgamma(1 + v, pol);
163 T f = (p - q) / v;
164 T g_prefix = boost::math::sin_pi(v / 2, pol);
165 g_prefix *= g_prefix * 2 / v;
166 T g = f + g_prefix * q;
167 T h = p;
168 T c_mult = -x * x / 4;
169
170 T y(c * g), y1(c * h);
171
172 for(int k = 1; k < policies::get_max_series_iterations<Policy>(); ++k)
173 {
174 f = (k * f + p + q) / (k*k - v*v);
175 p /= k - v;
176 q /= k + v;
177 c *= c_mult / k;
178 T c1 = pow(-x * x / 4, k) / factorial<T>(k, pol);
179 g = f + g_prefix * q;
180 h = -k * g + p;
181 y += c * g;
182 y1 += c * h;
183 if(c * g / tools::epsilon<T>() < y)
184 break;
185 }
186
187 *Y = -y;
188 *Y1 = (-2 / x) * y1;
189}
190
191template <class T, class Policy>
192T asymptotic_bessel_i_large_x(T v, T x, const Policy& pol)
193{
194 BOOST_MATH_STD_USING // ADL of std names
195 T s = 1;
196 T mu = 4 * v * v;
197 T ex = 8 * x;
198 T num = mu - 1;
199 T denom = ex;
200
201 s -= num / denom;
202
203 num *= mu - 9;
204 denom *= ex * 2;
205 s += num / denom;
206
207 num *= mu - 25;
208 denom *= ex * 3;
209 s -= num / denom;
210
211 // Try and avoid overflow to the last minute:
212 T e = exp(x/2);
213
214 s = e * (e * s / sqrt(2 * x * constants::pi<T>()));
215
216 return (boost::math::isfinite)(s) ?
217 s : policies::raise_overflow_error<T>("boost::math::asymptotic_bessel_i_large_x<%1%>(%1%,%1%)", 0, pol);
218}
219
220}}} // namespaces
221
222#endif
223
224

source code of include/boost/math/special_functions/detail/bessel_jy_asym.hpp