# GL Profile Suite / boost_1_51_0 / boost / math / special_functions / detail / bessel_jn.hpp

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