]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 |