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

source code of include/boost/math/special_functions/detail/bessel_y1.hpp