]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Copyright John Maddock 2012. |
2 | // Use, modification and distribution are subject to the | |
3 | // Boost Software License, Version 1.0. | |
4 | // (See accompanying file LICENSE_1_0.txt | |
5 | // or copy at http://www.boost.org/LICENSE_1_0.txt) | |
6 | ||
7 | #ifndef BOOST_MATH_AIRY_HPP | |
8 | #define BOOST_MATH_AIRY_HPP | |
9 | ||
10 | #include <limits> | |
11 | #include <boost/math/special_functions/math_fwd.hpp> | |
12 | #include <boost/math/special_functions/bessel.hpp> | |
13 | #include <boost/math/special_functions/cbrt.hpp> | |
14 | #include <boost/math/special_functions/detail/airy_ai_bi_zero.hpp> | |
15 | #include <boost/math/tools/roots.hpp> | |
16 | ||
17 | namespace boost{ namespace math{ | |
18 | ||
19 | namespace detail{ | |
20 | ||
21 | template <class T, class Policy> | |
22 | T airy_ai_imp(T x, const Policy& pol) | |
23 | { | |
24 | BOOST_MATH_STD_USING | |
25 | ||
26 | if(x < 0) | |
27 | { | |
28 | T p = (-x * sqrt(-x) * 2) / 3; | |
29 | T v = T(1) / 3; | |
30 | T j1 = boost::math::cyl_bessel_j(v, p, pol); | |
31 | T j2 = boost::math::cyl_bessel_j(-v, p, pol); | |
32 | T ai = sqrt(-x) * (j1 + j2) / 3; | |
33 | //T bi = sqrt(-x / 3) * (j2 - j1); | |
34 | return ai; | |
35 | } | |
36 | else if(fabs(x * x * x) / 6 < tools::epsilon<T>()) | |
37 | { | |
38 | T tg = boost::math::tgamma(constants::twothirds<T>(), pol); | |
39 | T ai = 1 / (pow(T(3), constants::twothirds<T>()) * tg); | |
40 | //T bi = 1 / (sqrt(boost::math::cbrt(T(3))) * tg); | |
41 | return ai; | |
42 | } | |
43 | else | |
44 | { | |
45 | T p = 2 * x * sqrt(x) / 3; | |
46 | T v = T(1) / 3; | |
47 | //T j1 = boost::math::cyl_bessel_i(-v, p, pol); | |
48 | //T j2 = boost::math::cyl_bessel_i(v, p, pol); | |
49 | // | |
50 | // Note that although we can calculate ai from j1 and j2, the accuracy is horrible | |
51 | // as we're subtracting two very large values, so use the Bessel K relation instead: | |
52 | // | |
53 | T ai = cyl_bessel_k(v, p, pol) * sqrt(x / 3) / boost::math::constants::pi<T>(); //sqrt(x) * (j1 - j2) / 3; | |
54 | //T bi = sqrt(x / 3) * (j1 + j2); | |
55 | return ai; | |
56 | } | |
57 | } | |
58 | ||
59 | template <class T, class Policy> | |
60 | T airy_bi_imp(T x, const Policy& pol) | |
61 | { | |
62 | BOOST_MATH_STD_USING | |
63 | ||
64 | if(x < 0) | |
65 | { | |
66 | T p = (-x * sqrt(-x) * 2) / 3; | |
67 | T v = T(1) / 3; | |
68 | T j1 = boost::math::cyl_bessel_j(v, p, pol); | |
69 | T j2 = boost::math::cyl_bessel_j(-v, p, pol); | |
70 | //T ai = sqrt(-x) * (j1 + j2) / 3; | |
71 | T bi = sqrt(-x / 3) * (j2 - j1); | |
72 | return bi; | |
73 | } | |
74 | else if(fabs(x * x * x) / 6 < tools::epsilon<T>()) | |
75 | { | |
76 | T tg = boost::math::tgamma(constants::twothirds<T>(), pol); | |
77 | //T ai = 1 / (pow(T(3), constants::twothirds<T>()) * tg); | |
20effc67 | 78 | T bi = 1 / (sqrt(boost::math::cbrt(T(3), pol)) * tg); |
7c673cae FG |
79 | return bi; |
80 | } | |
81 | else | |
82 | { | |
83 | T p = 2 * x * sqrt(x) / 3; | |
84 | T v = T(1) / 3; | |
85 | T j1 = boost::math::cyl_bessel_i(-v, p, pol); | |
86 | T j2 = boost::math::cyl_bessel_i(v, p, pol); | |
87 | T bi = sqrt(x / 3) * (j1 + j2); | |
88 | return bi; | |
89 | } | |
90 | } | |
91 | ||
92 | template <class T, class Policy> | |
93 | T airy_ai_prime_imp(T x, const Policy& pol) | |
94 | { | |
95 | BOOST_MATH_STD_USING | |
96 | ||
97 | if(x < 0) | |
98 | { | |
99 | T p = (-x * sqrt(-x) * 2) / 3; | |
100 | T v = T(2) / 3; | |
101 | T j1 = boost::math::cyl_bessel_j(v, p, pol); | |
102 | T j2 = boost::math::cyl_bessel_j(-v, p, pol); | |
103 | T aip = -x * (j1 - j2) / 3; | |
104 | return aip; | |
105 | } | |
106 | else if(fabs(x * x) / 2 < tools::epsilon<T>()) | |
107 | { | |
108 | T tg = boost::math::tgamma(constants::third<T>(), pol); | |
20effc67 | 109 | T aip = 1 / (boost::math::cbrt(T(3), pol) * tg); |
7c673cae FG |
110 | return -aip; |
111 | } | |
112 | else | |
113 | { | |
114 | T p = 2 * x * sqrt(x) / 3; | |
115 | T v = T(2) / 3; | |
116 | //T j1 = boost::math::cyl_bessel_i(-v, p, pol); | |
117 | //T j2 = boost::math::cyl_bessel_i(v, p, pol); | |
118 | // | |
119 | // Note that although we can calculate ai from j1 and j2, the accuracy is horrible | |
120 | // as we're subtracting two very large values, so use the Bessel K relation instead: | |
121 | // | |
122 | T aip = -cyl_bessel_k(v, p, pol) * x / (boost::math::constants::root_three<T>() * boost::math::constants::pi<T>()); | |
123 | return aip; | |
124 | } | |
125 | } | |
126 | ||
127 | template <class T, class Policy> | |
128 | T airy_bi_prime_imp(T x, const Policy& pol) | |
129 | { | |
130 | BOOST_MATH_STD_USING | |
131 | ||
132 | if(x < 0) | |
133 | { | |
134 | T p = (-x * sqrt(-x) * 2) / 3; | |
135 | T v = T(2) / 3; | |
136 | T j1 = boost::math::cyl_bessel_j(v, p, pol); | |
137 | T j2 = boost::math::cyl_bessel_j(-v, p, pol); | |
138 | T aip = -x * (j1 + j2) / constants::root_three<T>(); | |
139 | return aip; | |
140 | } | |
141 | else if(fabs(x * x) / 2 < tools::epsilon<T>()) | |
142 | { | |
143 | T tg = boost::math::tgamma(constants::third<T>(), pol); | |
20effc67 | 144 | T bip = sqrt(boost::math::cbrt(T(3), pol)) / tg; |
7c673cae FG |
145 | return bip; |
146 | } | |
147 | else | |
148 | { | |
149 | T p = 2 * x * sqrt(x) / 3; | |
150 | T v = T(2) / 3; | |
151 | T j1 = boost::math::cyl_bessel_i(-v, p, pol); | |
152 | T j2 = boost::math::cyl_bessel_i(v, p, pol); | |
153 | T aip = x * (j1 + j2) / boost::math::constants::root_three<T>(); | |
154 | return aip; | |
155 | } | |
156 | } | |
157 | ||
158 | template <class T, class Policy> | |
159 | T airy_ai_zero_imp(int m, const Policy& pol) | |
160 | { | |
161 | BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt. | |
162 | ||
163 | // Handle cases when a negative zero (negative rank) is requested. | |
164 | if(m < 0) | |
165 | { | |
166 | return policies::raise_domain_error<T>("boost::math::airy_ai_zero<%1%>(%1%, int)", | |
167 | "Requested the %1%'th zero, but the rank must be 1 or more !", static_cast<T>(m), pol); | |
168 | } | |
169 | ||
170 | // Handle case when the zero'th zero is requested. | |
171 | if(m == 0U) | |
172 | { | |
173 | return policies::raise_domain_error<T>("boost::math::airy_ai_zero<%1%>(%1%,%1%)", | |
174 | "The requested rank of the zero is %1%, but must be 1 or more !", static_cast<T>(m), pol); | |
175 | } | |
176 | ||
177 | // Set up the initial guess for the upcoming root-finding. | |
20effc67 | 178 | const T guess_root = boost::math::detail::airy_zero::airy_ai_zero_detail::initial_guess<T>(m, pol); |
7c673cae FG |
179 | |
180 | // Select the maximum allowed iterations based on the number | |
181 | // of decimal digits in the numeric type T, being at least 12. | |
182 | const int my_digits10 = static_cast<int>(static_cast<float>(policies::digits<T, Policy>() * 0.301F)); | |
183 | ||
184 | const boost::uintmax_t iterations_allowed = static_cast<boost::uintmax_t>((std::max)(12, my_digits10 * 2)); | |
185 | ||
186 | boost::uintmax_t iterations_used = iterations_allowed; | |
187 | ||
188 | // Use a dynamic tolerance because the roots get closer the higher m gets. | |
189 | T tolerance; | |
190 | ||
191 | if (m <= 10) { tolerance = T(0.3F); } | |
192 | else if(m <= 100) { tolerance = T(0.1F); } | |
193 | else if(m <= 1000) { tolerance = T(0.05F); } | |
194 | else { tolerance = T(1) / sqrt(T(m)); } | |
195 | ||
196 | // Perform the root-finding using Newton-Raphson iteration from Boost.Math. | |
197 | const T am = | |
198 | boost::math::tools::newton_raphson_iterate( | |
199 | boost::math::detail::airy_zero::airy_ai_zero_detail::function_object_ai_and_ai_prime<T, Policy>(pol), | |
200 | guess_root, | |
201 | T(guess_root - tolerance), | |
202 | T(guess_root + tolerance), | |
203 | policies::digits<T, Policy>(), | |
204 | iterations_used); | |
205 | ||
206 | static_cast<void>(iterations_used); | |
207 | ||
208 | return am; | |
209 | } | |
210 | ||
211 | template <class T, class Policy> | |
212 | T airy_bi_zero_imp(int m, const Policy& pol) | |
213 | { | |
214 | BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt. | |
215 | ||
216 | // Handle cases when a negative zero (negative rank) is requested. | |
217 | if(m < 0) | |
218 | { | |
219 | return policies::raise_domain_error<T>("boost::math::airy_bi_zero<%1%>(%1%, int)", | |
220 | "Requested the %1%'th zero, but the rank must 1 or more !", static_cast<T>(m), pol); | |
221 | } | |
222 | ||
223 | // Handle case when the zero'th zero is requested. | |
224 | if(m == 0U) | |
225 | { | |
226 | return policies::raise_domain_error<T>("boost::math::airy_bi_zero<%1%>(%1%,%1%)", | |
227 | "The requested rank of the zero is %1%, but must be 1 or more !", static_cast<T>(m), pol); | |
228 | } | |
229 | // Set up the initial guess for the upcoming root-finding. | |
20effc67 | 230 | const T guess_root = boost::math::detail::airy_zero::airy_bi_zero_detail::initial_guess<T>(m, pol); |
7c673cae FG |
231 | |
232 | // Select the maximum allowed iterations based on the number | |
233 | // of decimal digits in the numeric type T, being at least 12. | |
234 | const int my_digits10 = static_cast<int>(static_cast<float>(policies::digits<T, Policy>() * 0.301F)); | |
235 | ||
236 | const boost::uintmax_t iterations_allowed = static_cast<boost::uintmax_t>((std::max)(12, my_digits10 * 2)); | |
237 | ||
238 | boost::uintmax_t iterations_used = iterations_allowed; | |
239 | ||
240 | // Use a dynamic tolerance because the roots get closer the higher m gets. | |
241 | T tolerance; | |
242 | ||
243 | if (m <= 10) { tolerance = T(0.3F); } | |
244 | else if(m <= 100) { tolerance = T(0.1F); } | |
245 | else if(m <= 1000) { tolerance = T(0.05F); } | |
246 | else { tolerance = T(1) / sqrt(T(m)); } | |
247 | ||
248 | // Perform the root-finding using Newton-Raphson iteration from Boost.Math. | |
249 | const T bm = | |
250 | boost::math::tools::newton_raphson_iterate( | |
251 | boost::math::detail::airy_zero::airy_bi_zero_detail::function_object_bi_and_bi_prime<T, Policy>(pol), | |
252 | guess_root, | |
253 | T(guess_root - tolerance), | |
254 | T(guess_root + tolerance), | |
255 | policies::digits<T, Policy>(), | |
256 | iterations_used); | |
257 | ||
258 | static_cast<void>(iterations_used); | |
259 | ||
260 | return bm; | |
261 | } | |
262 | ||
263 | } // namespace detail | |
264 | ||
265 | template <class T, class Policy> | |
266 | inline typename tools::promote_args<T>::type airy_ai(T x, const Policy&) | |
267 | { | |
268 | BOOST_FPU_EXCEPTION_GUARD | |
269 | typedef typename tools::promote_args<T>::type result_type; | |
270 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
271 | typedef typename policies::normalise< | |
272 | Policy, | |
273 | policies::promote_float<false>, | |
274 | policies::promote_double<false>, | |
275 | policies::discrete_quantile<>, | |
276 | policies::assert_undefined<> >::type forwarding_policy; | |
277 | ||
278 | return policies::checked_narrowing_cast<result_type, Policy>(detail::airy_ai_imp<value_type>(static_cast<value_type>(x), forwarding_policy()), "boost::math::airy<%1%>(%1%)"); | |
279 | } | |
280 | ||
281 | template <class T> | |
282 | inline typename tools::promote_args<T>::type airy_ai(T x) | |
283 | { | |
284 | return airy_ai(x, policies::policy<>()); | |
285 | } | |
286 | ||
287 | template <class T, class Policy> | |
288 | inline typename tools::promote_args<T>::type airy_bi(T x, const Policy&) | |
289 | { | |
290 | BOOST_FPU_EXCEPTION_GUARD | |
291 | typedef typename tools::promote_args<T>::type result_type; | |
292 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
293 | typedef typename policies::normalise< | |
294 | Policy, | |
295 | policies::promote_float<false>, | |
296 | policies::promote_double<false>, | |
297 | policies::discrete_quantile<>, | |
298 | policies::assert_undefined<> >::type forwarding_policy; | |
299 | ||
300 | return policies::checked_narrowing_cast<result_type, Policy>(detail::airy_bi_imp<value_type>(static_cast<value_type>(x), forwarding_policy()), "boost::math::airy<%1%>(%1%)"); | |
301 | } | |
302 | ||
303 | template <class T> | |
304 | inline typename tools::promote_args<T>::type airy_bi(T x) | |
305 | { | |
306 | return airy_bi(x, policies::policy<>()); | |
307 | } | |
308 | ||
309 | template <class T, class Policy> | |
310 | inline typename tools::promote_args<T>::type airy_ai_prime(T x, const Policy&) | |
311 | { | |
312 | BOOST_FPU_EXCEPTION_GUARD | |
313 | typedef typename tools::promote_args<T>::type result_type; | |
314 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
315 | typedef typename policies::normalise< | |
316 | Policy, | |
317 | policies::promote_float<false>, | |
318 | policies::promote_double<false>, | |
319 | policies::discrete_quantile<>, | |
320 | policies::assert_undefined<> >::type forwarding_policy; | |
321 | ||
322 | return policies::checked_narrowing_cast<result_type, Policy>(detail::airy_ai_prime_imp<value_type>(static_cast<value_type>(x), forwarding_policy()), "boost::math::airy<%1%>(%1%)"); | |
323 | } | |
324 | ||
325 | template <class T> | |
326 | inline typename tools::promote_args<T>::type airy_ai_prime(T x) | |
327 | { | |
328 | return airy_ai_prime(x, policies::policy<>()); | |
329 | } | |
330 | ||
331 | template <class T, class Policy> | |
332 | inline typename tools::promote_args<T>::type airy_bi_prime(T x, const Policy&) | |
333 | { | |
334 | BOOST_FPU_EXCEPTION_GUARD | |
335 | typedef typename tools::promote_args<T>::type result_type; | |
336 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
337 | typedef typename policies::normalise< | |
338 | Policy, | |
339 | policies::promote_float<false>, | |
340 | policies::promote_double<false>, | |
341 | policies::discrete_quantile<>, | |
342 | policies::assert_undefined<> >::type forwarding_policy; | |
343 | ||
344 | return policies::checked_narrowing_cast<result_type, Policy>(detail::airy_bi_prime_imp<value_type>(static_cast<value_type>(x), forwarding_policy()), "boost::math::airy<%1%>(%1%)"); | |
345 | } | |
346 | ||
347 | template <class T> | |
348 | inline typename tools::promote_args<T>::type airy_bi_prime(T x) | |
349 | { | |
350 | return airy_bi_prime(x, policies::policy<>()); | |
351 | } | |
352 | ||
353 | template <class T, class Policy> | |
354 | inline T airy_ai_zero(int m, const Policy& /*pol*/) | |
355 | { | |
356 | BOOST_FPU_EXCEPTION_GUARD | |
357 | typedef typename policies::evaluation<T, Policy>::type value_type; | |
358 | typedef typename policies::normalise< | |
359 | Policy, | |
360 | policies::promote_float<false>, | |
361 | policies::promote_double<false>, | |
362 | policies::discrete_quantile<>, | |
363 | policies::assert_undefined<> >::type forwarding_policy; | |
364 | ||
365 | BOOST_STATIC_ASSERT_MSG( false == std::numeric_limits<T>::is_specialized | |
366 | || ( true == std::numeric_limits<T>::is_specialized | |
367 | && false == std::numeric_limits<T>::is_integer), | |
368 | "Airy value type must be a floating-point type."); | |
369 | ||
370 | return policies::checked_narrowing_cast<T, Policy>(detail::airy_ai_zero_imp<value_type>(m, forwarding_policy()), "boost::math::airy_ai_zero<%1%>(unsigned)"); | |
371 | } | |
372 | ||
373 | template <class T> | |
374 | inline T airy_ai_zero(int m) | |
375 | { | |
376 | return airy_ai_zero<T>(m, policies::policy<>()); | |
377 | } | |
378 | ||
379 | template <class T, class OutputIterator, class Policy> | |
380 | inline OutputIterator airy_ai_zero( | |
381 | int start_index, | |
382 | unsigned number_of_zeros, | |
383 | OutputIterator out_it, | |
384 | const Policy& pol) | |
385 | { | |
386 | typedef T result_type; | |
387 | ||
388 | BOOST_STATIC_ASSERT_MSG( false == std::numeric_limits<T>::is_specialized | |
389 | || ( true == std::numeric_limits<T>::is_specialized | |
390 | && false == std::numeric_limits<T>::is_integer), | |
391 | "Airy value type must be a floating-point type."); | |
392 | ||
393 | for(unsigned i = 0; i < number_of_zeros; ++i) | |
394 | { | |
395 | *out_it = boost::math::airy_ai_zero<result_type>(start_index + i, pol); | |
396 | ++out_it; | |
397 | } | |
398 | return out_it; | |
399 | } | |
400 | ||
401 | template <class T, class OutputIterator> | |
402 | inline OutputIterator airy_ai_zero( | |
403 | int start_index, | |
404 | unsigned number_of_zeros, | |
405 | OutputIterator out_it) | |
406 | { | |
407 | return airy_ai_zero<T>(start_index, number_of_zeros, out_it, policies::policy<>()); | |
408 | } | |
409 | ||
410 | template <class T, class Policy> | |
411 | inline T airy_bi_zero(int m, const Policy& /*pol*/) | |
412 | { | |
413 | BOOST_FPU_EXCEPTION_GUARD | |
414 | typedef typename policies::evaluation<T, Policy>::type value_type; | |
415 | typedef typename policies::normalise< | |
416 | Policy, | |
417 | policies::promote_float<false>, | |
418 | policies::promote_double<false>, | |
419 | policies::discrete_quantile<>, | |
420 | policies::assert_undefined<> >::type forwarding_policy; | |
421 | ||
422 | BOOST_STATIC_ASSERT_MSG( false == std::numeric_limits<T>::is_specialized | |
423 | || ( true == std::numeric_limits<T>::is_specialized | |
424 | && false == std::numeric_limits<T>::is_integer), | |
425 | "Airy value type must be a floating-point type."); | |
426 | ||
427 | return policies::checked_narrowing_cast<T, Policy>(detail::airy_bi_zero_imp<value_type>(m, forwarding_policy()), "boost::math::airy_bi_zero<%1%>(unsigned)"); | |
428 | } | |
429 | ||
430 | template <typename T> | |
431 | inline T airy_bi_zero(int m) | |
432 | { | |
433 | return airy_bi_zero<T>(m, policies::policy<>()); | |
434 | } | |
435 | ||
436 | template <class T, class OutputIterator, class Policy> | |
437 | inline OutputIterator airy_bi_zero( | |
438 | int start_index, | |
439 | unsigned number_of_zeros, | |
440 | OutputIterator out_it, | |
441 | const Policy& pol) | |
442 | { | |
443 | typedef T result_type; | |
444 | ||
445 | BOOST_STATIC_ASSERT_MSG( false == std::numeric_limits<T>::is_specialized | |
446 | || ( true == std::numeric_limits<T>::is_specialized | |
447 | && false == std::numeric_limits<T>::is_integer), | |
448 | "Airy value type must be a floating-point type."); | |
449 | ||
450 | for(unsigned i = 0; i < number_of_zeros; ++i) | |
451 | { | |
452 | *out_it = boost::math::airy_bi_zero<result_type>(start_index + i, pol); | |
453 | ++out_it; | |
454 | } | |
455 | return out_it; | |
456 | } | |
457 | ||
458 | template <class T, class OutputIterator> | |
459 | inline OutputIterator airy_bi_zero( | |
460 | int start_index, | |
461 | unsigned number_of_zeros, | |
462 | OutputIterator out_it) | |
463 | { | |
464 | return airy_bi_zero<T>(start_index, number_of_zeros, out_it, policies::policy<>()); | |
465 | } | |
466 | ||
467 | }} // namespaces | |
468 | ||
469 | #endif // BOOST_MATH_AIRY_HPP |