1 | // Copyright (c) 2006 Xiaogang Zhang |
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 | #ifndef BOOST_MATH_BESSEL_JN_HPP |
7 | #define BOOST_MATH_BESSEL_JN_HPP |
8 | |
9 | #ifdef _MSC_VER |
10 | #pragma once |
11 | #endif |
12 | |
13 | #include <boost/math/special_functions/detail/bessel_j0.hpp> |
14 | #include <boost/math/special_functions/detail/bessel_j1.hpp> |
15 | #include <boost/math/special_functions/detail/bessel_jy.hpp> |
16 | #include <boost/math/special_functions/detail/bessel_jy_asym.hpp> |
17 | #include <boost/math/special_functions/detail/bessel_jy_series.hpp> |
18 | |
19 | // Bessel function of the first kind of integer order |
20 | // J_n(z) is the minimal solution |
21 | // n < abs(z), forward recurrence stable and usable |
22 | // n >= abs(z), forward recurrence unstable, use Miller's algorithm |
23 | |
24 | namespace boost { namespace math { namespace detail{ |
25 | |
26 | template <typename T, typename Policy> |
27 | T bessel_jn(int n, T x, const Policy& pol) |
28 | { |
29 | T value(0), factor, current, prev, next; |
30 | |
31 | BOOST_MATH_STD_USING |
32 | |
33 | // |
34 | // Reflection has to come first: |
35 | // |
36 | if (n < 0) |
37 | { |
38 | factor = static_cast<T>((n & 0x1) ? -1 : 1); // J_{-n}(z) = (-1)^n J_n(z) |
39 | n = -n; |
40 | } |
41 | else |
42 | { |
43 | factor = 1; |
44 | } |
45 | if(x < 0) |
46 | { |
47 | factor *= (n & 0x1) ? -1 : 1; // J_{n}(-z) = (-1)^n J_n(z) |
48 | x = -x; |
49 | } |
50 | // |
51 | // Special cases: |
52 | // |
53 | if(asymptotic_bessel_large_x_limit(T(n), x)) |
54 | return factor * asymptotic_bessel_j_large_x_2<T>(T(n), x); |
55 | if (n == 0) |
56 | { |
57 | return factor * bessel_j0(x); |
58 | } |
59 | if (n == 1) |
60 | { |
61 | return factor * bessel_j1(x); |
62 | } |
63 | |
64 | if (x == 0) // n >= 2 |
65 | { |
66 | return static_cast<T>(0); |
67 | } |
68 | |
69 | BOOST_ASSERT(n > 1); |
70 | T scale = 1; |
71 | if (n < abs(x)) // forward recurrence |
72 | { |
73 | prev = bessel_j0(x); |
74 | current = bessel_j1(x); |
75 | policies::check_series_iterations<T>("boost::math::bessel_j_n<%1%>(%1%,%1%)" , n, pol); |
76 | for (int k = 1; k < n; k++) |
77 | { |
78 | T fact = 2 * k / x; |
79 | // |
80 | // rescale if we would overflow or underflow: |
81 | // |
82 | if((fabs(fact) > 1) && ((tools::max_value<T>() - fabs(prev)) / fabs(fact) < fabs(current))) |
83 | { |
84 | scale /= current; |
85 | prev /= current; |
86 | current = 1; |
87 | } |
88 | value = fact * current - prev; |
89 | prev = current; |
90 | current = value; |
91 | } |
92 | } |
93 | else if((x < 1) || (n > x * x / 4) || (x < 5)) |
94 | { |
95 | return factor * bessel_j_small_z_series(T(n), x, pol); |
96 | } |
97 | else // backward recurrence |
98 | { |
99 | T fn; int s; // fn = J_(n+1) / J_n |
100 | // |x| <= n, fast convergence for continued fraction CF1 |
101 | boost::math::detail::CF1_jy(static_cast<T>(n), x, &fn, &s, pol); |
102 | prev = fn; |
103 | current = 1; |
104 | // Check recursion won't go on too far: |
105 | policies::check_series_iterations<T>("boost::math::bessel_j_n<%1%>(%1%,%1%)" , n, pol); |
106 | for (int k = n; k > 0; k--) |
107 | { |
108 | T fact = 2 * k / x; |
109 | if((fabs(fact) > 1) && ((tools::max_value<T>() - fabs(prev)) / fabs(fact) < fabs(current))) |
110 | { |
111 | prev /= current; |
112 | scale /= current; |
113 | current = 1; |
114 | } |
115 | next = fact * current - prev; |
116 | prev = current; |
117 | current = next; |
118 | } |
119 | value = bessel_j0(x) / current; // normalization |
120 | scale = 1 / scale; |
121 | } |
122 | value *= factor; |
123 | |
124 | if(tools::max_value<T>() * scale < fabs(value)) |
125 | return policies::raise_overflow_error<T>("boost::math::bessel_jn<%1%>(%1%,%1%)" , 0, pol); |
126 | |
127 | return value / scale; |
128 | } |
129 | |
130 | }}} // namespaces |
131 | |
132 | #endif // BOOST_MATH_BESSEL_JN_HPP |
133 | |
134 | |