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)
7 #ifndef BOOST_MATH_AIRY_HPP
8 #define BOOST_MATH_AIRY_HPP
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>
17 namespace boost{ namespace math{
21 template <class T, class Policy>
22 T airy_ai_imp(T x, const Policy& pol)
28 T p = (-x * sqrt(-x) * 2) / 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);
36 else if(fabs(x * x * x) / 6 < tools::epsilon<T>())
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);
45 T p = 2 * x * sqrt(x) / 3;
47 //T j1 = boost::math::cyl_bessel_i(-v, p, pol);
48 //T j2 = boost::math::cyl_bessel_i(v, p, pol);
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:
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);
59 template <class T, class Policy>
60 T airy_bi_imp(T x, const Policy& pol)
66 T p = (-x * sqrt(-x) * 2) / 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);
74 else if(fabs(x * x * x) / 6 < tools::epsilon<T>())
76 T tg = boost::math::tgamma(constants::twothirds<T>(), pol);
77 //T ai = 1 / (pow(T(3), constants::twothirds<T>()) * tg);
78 T bi = 1 / (sqrt(boost::math::cbrt(T(3))) * tg);
83 T p = 2 * x * sqrt(x) / 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);
92 template <class T, class Policy>
93 T airy_ai_prime_imp(T x, const Policy& pol)
99 T p = (-x * sqrt(-x) * 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;
106 else if(fabs(x * x) / 2 < tools::epsilon<T>())
108 T tg = boost::math::tgamma(constants::third<T>(), pol);
109 T aip = 1 / (boost::math::cbrt(T(3)) * tg);
114 T p = 2 * x * sqrt(x) / 3;
116 //T j1 = boost::math::cyl_bessel_i(-v, p, pol);
117 //T j2 = boost::math::cyl_bessel_i(v, p, pol);
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:
122 T aip = -cyl_bessel_k(v, p, pol) * x / (boost::math::constants::root_three<T>() * boost::math::constants::pi<T>());
127 template <class T, class Policy>
128 T airy_bi_prime_imp(T x, const Policy& pol)
134 T p = (-x * sqrt(-x) * 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>();
141 else if(fabs(x * x) / 2 < tools::epsilon<T>())
143 T tg = boost::math::tgamma(constants::third<T>(), pol);
144 T bip = sqrt(boost::math::cbrt(T(3))) / tg;
149 T p = 2 * x * sqrt(x) / 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>();
158 template <class T, class Policy>
159 T airy_ai_zero_imp(int m, const Policy& pol)
161 BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt.
163 // Handle cases when a negative zero (negative rank) is requested.
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);
170 // Handle case when the zero'th zero is requested.
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);
177 // Set up the initial guess for the upcoming root-finding.
178 const T guess_root = boost::math::detail::airy_zero::airy_ai_zero_detail::initial_guess<T>(m);
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));
184 const boost::uintmax_t iterations_allowed = static_cast<boost::uintmax_t>((std::max)(12, my_digits10 * 2));
186 boost::uintmax_t iterations_used = iterations_allowed;
188 // Use a dynamic tolerance because the roots get closer the higher m gets.
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)); }
196 // Perform the root-finding using Newton-Raphson iteration from Boost.Math.
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),
201 T(guess_root - tolerance),
202 T(guess_root + tolerance),
203 policies::digits<T, Policy>(),
206 static_cast<void>(iterations_used);
211 template <class T, class Policy>
212 T airy_bi_zero_imp(int m, const Policy& pol)
214 BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt.
216 // Handle cases when a negative zero (negative rank) is requested.
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);
223 // Handle case when the zero'th zero is requested.
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);
229 // Set up the initial guess for the upcoming root-finding.
230 const T guess_root = boost::math::detail::airy_zero::airy_bi_zero_detail::initial_guess<T>(m);
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));
236 const boost::uintmax_t iterations_allowed = static_cast<boost::uintmax_t>((std::max)(12, my_digits10 * 2));
238 boost::uintmax_t iterations_used = iterations_allowed;
240 // Use a dynamic tolerance because the roots get closer the higher m gets.
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)); }
248 // Perform the root-finding using Newton-Raphson iteration from Boost.Math.
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),
253 T(guess_root - tolerance),
254 T(guess_root + tolerance),
255 policies::digits<T, Policy>(),
258 static_cast<void>(iterations_used);
263 } // namespace detail
265 template <class T, class Policy>
266 inline typename tools::promote_args<T>::type airy_ai(T x, const Policy&)
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<
273 policies::promote_float<false>,
274 policies::promote_double<false>,
275 policies::discrete_quantile<>,
276 policies::assert_undefined<> >::type forwarding_policy;
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%)");
282 inline typename tools::promote_args<T>::type airy_ai(T x)
284 return airy_ai(x, policies::policy<>());
287 template <class T, class Policy>
288 inline typename tools::promote_args<T>::type airy_bi(T x, const Policy&)
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<
295 policies::promote_float<false>,
296 policies::promote_double<false>,
297 policies::discrete_quantile<>,
298 policies::assert_undefined<> >::type forwarding_policy;
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%)");
304 inline typename tools::promote_args<T>::type airy_bi(T x)
306 return airy_bi(x, policies::policy<>());
309 template <class T, class Policy>
310 inline typename tools::promote_args<T>::type airy_ai_prime(T x, const Policy&)
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<
317 policies::promote_float<false>,
318 policies::promote_double<false>,
319 policies::discrete_quantile<>,
320 policies::assert_undefined<> >::type forwarding_policy;
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%)");
326 inline typename tools::promote_args<T>::type airy_ai_prime(T x)
328 return airy_ai_prime(x, policies::policy<>());
331 template <class T, class Policy>
332 inline typename tools::promote_args<T>::type airy_bi_prime(T x, const Policy&)
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<
339 policies::promote_float<false>,
340 policies::promote_double<false>,
341 policies::discrete_quantile<>,
342 policies::assert_undefined<> >::type forwarding_policy;
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%)");
348 inline typename tools::promote_args<T>::type airy_bi_prime(T x)
350 return airy_bi_prime(x, policies::policy<>());
353 template <class T, class Policy>
354 inline T airy_ai_zero(int m, const Policy& /*pol*/)
356 BOOST_FPU_EXCEPTION_GUARD
357 typedef typename policies::evaluation<T, Policy>::type value_type;
358 typedef typename policies::normalise<
360 policies::promote_float<false>,
361 policies::promote_double<false>,
362 policies::discrete_quantile<>,
363 policies::assert_undefined<> >::type forwarding_policy;
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.");
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)");
374 inline T airy_ai_zero(int m)
376 return airy_ai_zero<T>(m, policies::policy<>());
379 template <class T, class OutputIterator, class Policy>
380 inline OutputIterator airy_ai_zero(
382 unsigned number_of_zeros,
383 OutputIterator out_it,
386 typedef T result_type;
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.");
393 for(unsigned i = 0; i < number_of_zeros; ++i)
395 *out_it = boost::math::airy_ai_zero<result_type>(start_index + i, pol);
401 template <class T, class OutputIterator>
402 inline OutputIterator airy_ai_zero(
404 unsigned number_of_zeros,
405 OutputIterator out_it)
407 return airy_ai_zero<T>(start_index, number_of_zeros, out_it, policies::policy<>());
410 template <class T, class Policy>
411 inline T airy_bi_zero(int m, const Policy& /*pol*/)
413 BOOST_FPU_EXCEPTION_GUARD
414 typedef typename policies::evaluation<T, Policy>::type value_type;
415 typedef typename policies::normalise<
417 policies::promote_float<false>,
418 policies::promote_double<false>,
419 policies::discrete_quantile<>,
420 policies::assert_undefined<> >::type forwarding_policy;
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.");
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)");
430 template <typename T>
431 inline T airy_bi_zero(int m)
433 return airy_bi_zero<T>(m, policies::policy<>());
436 template <class T, class OutputIterator, class Policy>
437 inline OutputIterator airy_bi_zero(
439 unsigned number_of_zeros,
440 OutputIterator out_it,
443 typedef T result_type;
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.");
450 for(unsigned i = 0; i < number_of_zeros; ++i)
452 *out_it = boost::math::airy_bi_zero<result_type>(start_index + i, pol);
458 template <class T, class OutputIterator>
459 inline OutputIterator airy_bi_zero(
461 unsigned number_of_zeros,
462 OutputIterator out_it)
464 return airy_bi_zero<T>(start_index, number_of_zeros, out_it, policies::policy<>());
469 #endif // BOOST_MATH_AIRY_HPP