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
24namespace boost { namespace math { namespace detail{
25
26template <typename T, typename Policy>
27T 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, pol);
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_MATH_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 value = (2 * k * current / x) - prev;
79 prev = current;
80 current = value;
81 }
82 }
83 else if((x < 1) || (n > x * x / 4) || (x < 5))
84 {
85 return factor * bessel_j_small_z_series(T(n), x, pol);
86 }
87 else // backward recurrence
88 {
89 T fn; int s; // fn = J_(n+1) / J_n
90 // |x| <= n, fast convergence for continued fraction CF1
91 boost::math::detail::CF1_jy(static_cast<T>(n), x, &fn, &s, pol);
92 prev = fn;
93 current = 1;
94 // Check recursion won't go on too far:
95 policies::check_series_iterations<T>("boost::math::bessel_j_n<%1%>(%1%,%1%)", n, pol);
96 for (int k = n; k > 0; k--)
97 {
98 T fact = 2 * k / x;
99 if((fabs(fact) > 1) && ((tools::max_value<T>() - fabs(prev)) / fabs(fact) < fabs(current)))
100 {
101 prev /= current;
102 scale /= current;
103 current = 1;
104 }
105 next = fact * current - prev;
106 prev = current;
107 current = next;
108 }
109 value = bessel_j0(x) / current; // normalization
110 scale = 1 / scale;
111 }
112 value *= factor;
113
114 if(tools::max_value<T>() * scale < fabs(value))
115 return policies::raise_overflow_error<T>("boost::math::bessel_jn<%1%>(%1%,%1%)", nullptr, pol); // LCOV_EXCL_LINE we should never get here!
116
117 return value / scale;
118}
119
120}}} // namespaces
121
122#endif // BOOST_MATH_BESSEL_JN_HPP
123
124

source code of boost/libs/math/include/boost/math/special_functions/detail/bessel_jn.hpp