]>
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_Y1_HPP | |
7 | #define BOOST_MATH_BESSEL_Y1_HPP | |
8 | ||
9 | #ifdef _MSC_VER | |
10 | #pragma once | |
11 | #pragma warning(push) | |
12 | #pragma warning(disable:4702) // Unreachable code (release mode only warning) | |
13 | #endif | |
14 | ||
15 | #include <boost/math/special_functions/detail/bessel_j1.hpp> | |
16 | #include <boost/math/constants/constants.hpp> | |
17 | #include <boost/math/tools/rational.hpp> | |
18 | #include <boost/math/tools/big_constant.hpp> | |
19 | #include <boost/math/policies/error_handling.hpp> | |
20 | #include <boost/assert.hpp> | |
21 | ||
22 | // Bessel function of the second kind of order one | |
23 | // x <= 8, minimax rational approximations on root-bracketing intervals | |
24 | // x > 8, Hankel asymptotic expansion in Hart, Computer Approximations, 1968 | |
25 | ||
26 | namespace boost { namespace math { namespace detail{ | |
27 | ||
28 | template <typename T, typename Policy> | |
29 | T bessel_y1(T x, const Policy&); | |
30 | ||
31 | template <class T, class Policy> | |
32 | struct bessel_y1_initializer | |
33 | { | |
34 | struct init | |
35 | { | |
36 | init() | |
37 | { | |
38 | do_init(); | |
39 | } | |
40 | static void do_init() | |
41 | { | |
42 | bessel_y1(T(1), Policy()); | |
43 | } | |
44 | void force_instantiate()const{} | |
45 | }; | |
46 | static const init initializer; | |
47 | static void force_instantiate() | |
48 | { | |
49 | initializer.force_instantiate(); | |
50 | } | |
51 | }; | |
52 | ||
53 | template <class T, class Policy> | |
54 | const typename bessel_y1_initializer<T, Policy>::init bessel_y1_initializer<T, Policy>::initializer; | |
55 | ||
56 | template <typename T, typename Policy> | |
57 | T bessel_y1(T x, const Policy& pol) | |
58 | { | |
59 | bessel_y1_initializer<T, Policy>::force_instantiate(); | |
60 | ||
61 | static const T P1[] = { | |
62 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0535726612579544093e+13)), | |
63 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.4708611716525426053e+12)), | |
64 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.7595974497819597599e+11)), | |
65 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.2144548214502560419e+09)), | |
66 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.9157479997408395984e+07)), | |
67 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.2157953222280260820e+05)), | |
68 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.1714424660046133456e+02)), | |
69 | }; | |
70 | static const T Q1[] = { | |
71 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.0737873921079286084e+14)), | |
72 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.1272286200406461981e+12)), | |
73 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.7800352738690585613e+10)), | |
74 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.2250435122182963220e+08)), | |
75 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.8136470753052572164e+05)), | |
76 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.2079908168393867438e+02)), | |
77 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)), | |
78 | }; | |
79 | static const T P2[] = { | |
80 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1514276357909013326e+19)), | |
81 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.6808094574724204577e+18)), | |
82 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.3638408497043134724e+16)), | |
83 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0686275289804744814e+15)), | |
84 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.9530713129741981618e+13)), | |
85 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.7453673962438488783e+11)), | |
86 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.1957961912070617006e+09)), | |
87 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.9153806858264202986e+06)), | |
88 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.2337180442012953128e+03)), | |
89 | }; | |
90 | static const T Q2[] = { | |
91 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.3321844313316185697e+20)), | |
92 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.6968198822857178911e+18)), | |
93 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.0837179548112881950e+16)), | |
94 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1187010065856971027e+14)), | |
95 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.0221766852960403645e+11)), | |
96 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.3550318087088919566e+08)), | |
97 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0453748201934079734e+06)), | |
98 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.2855164849321609336e+03)), | |
99 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)), | |
100 | }; | |
101 | static const T PC[] = { | |
102 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4357578167941278571e+06)), | |
103 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.9422465050776411957e+06)), | |
104 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.6033732483649391093e+06)), | |
105 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5235293511811373833e+06)), | |
106 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0982405543459346727e+05)), | |
107 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.6116166443246101165e+03)), | |
108 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0)), | |
109 | }; | |
110 | static const T QC[] = { | |
111 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4357578167941278568e+06)), | |
112 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.9341243899345856590e+06)), | |
113 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.5853394797230870728e+06)), | |
114 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5118095066341608816e+06)), | |
115 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0726385991103820119e+05)), | |
116 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4550094401904961825e+03)), | |
117 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)), | |
118 | }; | |
119 | static const T PS[] = { | |
120 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.3220913409857223519e+04)), | |
121 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.5145160675335701966e+04)), | |
122 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6178836581270835179e+04)), | |
123 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8494262873223866797e+04)), | |
124 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7063754290207680021e+03)), | |
125 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.5265133846636032186e+01)), | |
126 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0)), | |
127 | }; | |
128 | static const T QS[] = { | |
129 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.0871281941028743574e+05)), | |
130 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8194580422439972989e+06)), | |
131 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.4194606696037208929e+06)), | |
132 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0029443582266975117e+05)), | |
133 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.7890229745772202641e+04)), | |
134 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.6383677696049909675e+02)), | |
135 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)), | |
136 | }; | |
137 | static const T x1 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.1971413260310170351e+00)), | |
138 | x2 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.4296810407941351328e+00)), | |
139 | x11 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.620e+02)), | |
140 | x12 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8288260310170351490e-03)), | |
141 | x21 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3900e+03)), | |
142 | x22 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.4592058648672279948e-06)) | |
143 | ; | |
144 | T value, factor, r, rc, rs; | |
145 | ||
146 | BOOST_MATH_STD_USING | |
147 | using namespace boost::math::tools; | |
148 | using namespace boost::math::constants; | |
149 | ||
150 | if (x <= 0) | |
151 | { | |
152 | return policies::raise_domain_error<T>("bost::math::bessel_y1<%1%>(%1%,%1%)", | |
153 | "Got x == %1%, but x must be > 0, complex result not supported.", x, pol); | |
154 | } | |
155 | if (x <= 4) // x in (0, 4] | |
156 | { | |
157 | T y = x * x; | |
158 | T z = 2 * log(x/x1) * bessel_j1(x) / pi<T>(); | |
159 | r = evaluate_rational(P1, Q1, y); | |
160 | factor = (x + x1) * ((x - x11/256) - x12) / x; | |
161 | value = z + factor * r; | |
162 | } | |
163 | else if (x <= 8) // x in (4, 8] | |
164 | { | |
165 | T y = x * x; | |
166 | T z = 2 * log(x/x2) * bessel_j1(x) / pi<T>(); | |
167 | r = evaluate_rational(P2, Q2, y); | |
168 | factor = (x + x2) * ((x - x21/256) - x22) / x; | |
169 | value = z + factor * r; | |
170 | } | |
171 | else // x in (8, \infty) | |
172 | { | |
173 | T y = 8 / x; | |
174 | T y2 = y * y; | |
175 | rc = evaluate_rational(PC, QC, y2); | |
176 | rs = evaluate_rational(PS, QS, y2); | |
177 | factor = 1 / (sqrt(x) * root_pi<T>()); | |
178 | // | |
179 | // This code is really just: | |
180 | // | |
181 | // T z = x - 0.75f * pi<T>(); | |
182 | // value = factor * (rc * sin(z) + y * rs * cos(z)); | |
183 | // | |
184 | // But using the sin/cos addition rules, plus constants for sin/cos of 3PI/4 | |
185 | // which then cancel out with corresponding terms in "factor". | |
186 | // | |
187 | T sx = sin(x); | |
188 | T cx = cos(x); | |
189 | value = factor * (y * rs * (sx - cx) - rc * (sx + cx)); | |
190 | } | |
191 | ||
192 | return value; | |
193 | } | |
194 | ||
195 | }}} // namespaces | |
196 | ||
197 | #ifdef _MSC_VER | |
198 | #pragma warning(pop) | |
199 | #endif | |
200 | ||
201 | #endif // BOOST_MATH_BESSEL_Y1_HPP | |
202 |