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
10 #include <boost/math/special_functions/bessel.hpp>
11 #include <boost/math/special_functions/cbrt.hpp>
12 #include <boost/math/special_functions/detail/airy_ai_bi_zero.hpp>
13 #include <boost/math/tools/roots.hpp>
15 namespace boost{ namespace math{
19 template <class T, class Policy>
20 T airy_ai_imp(T x, const Policy& pol)
26 T p = (-x * sqrt(-x) * 2) / 3;
28 T j1 = boost::math::cyl_bessel_j(v, p, pol);
29 T j2 = boost::math::cyl_bessel_j(-v, p, pol);
30 T ai = sqrt(-x) * (j1 + j2) / 3;
31 //T bi = sqrt(-x / 3) * (j2 - j1);
34 else if(fabs(x * x * x) / 6 < tools::epsilon<T>())
36 T tg = boost::math::tgamma(constants::twothirds<T>(), pol);
37 T ai = 1 / (pow(T(3), constants::twothirds<T>()) * tg);
38 //T bi = 1 / (sqrt(boost::math::cbrt(T(3))) * tg);
43 T p = 2 * x * sqrt(x) / 3;
45 //T j1 = boost::math::cyl_bessel_i(-v, p, pol);
46 //T j2 = boost::math::cyl_bessel_i(v, p, pol);
48 // Note that although we can calculate ai from j1 and j2, the accuracy is horrible
49 // as we're subtracting two very large values, so use the Bessel K relation instead:
51 T ai = cyl_bessel_k(v, p, pol) * sqrt(x / 3) / boost::math::constants::pi<T>(); //sqrt(x) * (j1 - j2) / 3;
52 //T bi = sqrt(x / 3) * (j1 + j2);
57 template <class T, class Policy>
58 T airy_bi_imp(T x, const Policy& pol)
64 T p = (-x * sqrt(-x) * 2) / 3;
66 T j1 = boost::math::cyl_bessel_j(v, p, pol);
67 T j2 = boost::math::cyl_bessel_j(-v, p, pol);
68 //T ai = sqrt(-x) * (j1 + j2) / 3;
69 T bi = sqrt(-x / 3) * (j2 - j1);
72 else if(fabs(x * x * x) / 6 < tools::epsilon<T>())
74 T tg = boost::math::tgamma(constants::twothirds<T>(), pol);
75 //T ai = 1 / (pow(T(3), constants::twothirds<T>()) * tg);
76 T bi = 1 / (sqrt(boost::math::cbrt(T(3))) * tg);
81 T p = 2 * x * sqrt(x) / 3;
83 T j1 = boost::math::cyl_bessel_i(-v, p, pol);
84 T j2 = boost::math::cyl_bessel_i(v, p, pol);
85 T bi = sqrt(x / 3) * (j1 + j2);
90 template <class T, class Policy>
91 T airy_ai_prime_imp(T x, const Policy& pol)
97 T p = (-x * sqrt(-x) * 2) / 3;
99 T j1 = boost::math::cyl_bessel_j(v, p, pol);
100 T j2 = boost::math::cyl_bessel_j(-v, p, pol);
101 T aip = -x * (j1 - j2) / 3;
104 else if(fabs(x * x) / 2 < tools::epsilon<T>())
106 T tg = boost::math::tgamma(constants::third<T>(), pol);
107 T aip = 1 / (boost::math::cbrt(T(3)) * tg);
112 T p = 2 * x * sqrt(x) / 3;
114 //T j1 = boost::math::cyl_bessel_i(-v, p, pol);
115 //T j2 = boost::math::cyl_bessel_i(v, p, pol);
117 // Note that although we can calculate ai from j1 and j2, the accuracy is horrible
118 // as we're subtracting two very large values, so use the Bessel K relation instead:
120 T aip = -cyl_bessel_k(v, p, pol) * x / (boost::math::constants::root_three<T>() * boost::math::constants::pi<T>());
125 template <class T, class Policy>
126 T airy_bi_prime_imp(T x, const Policy& pol)
132 T p = (-x * sqrt(-x) * 2) / 3;
134 T j1 = boost::math::cyl_bessel_j(v, p, pol);
135 T j2 = boost::math::cyl_bessel_j(-v, p, pol);
136 T aip = -x * (j1 + j2) / constants::root_three<T>();
139 else if(fabs(x * x) / 2 < tools::epsilon<T>())
141 T tg = boost::math::tgamma(constants::third<T>(), pol);
142 T bip = sqrt(boost::math::cbrt(T(3))) / tg;
147 T p = 2 * x * sqrt(x) / 3;
149 T j1 = boost::math::cyl_bessel_i(-v, p, pol);
150 T j2 = boost::math::cyl_bessel_i(v, p, pol);
151 T aip = x * (j1 + j2) / boost::math::constants::root_three<T>();
156 template <class T, class Policy>
157 T airy_ai_zero_imp(unsigned m, const Policy& pol)
159 BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt.
161 // Handle case when the zero'th zero is requested.
164 return policies::raise_domain_error<T>("boost::math::airy_ai_zero<%1%>(%1%,%1%)",
165 "The requested rank of the zero is %1%, but must be 1 or more !", static_cast<T>(m), pol);
168 // Set up the initial guess for the upcoming root-finding.
169 const T guess_root = boost::math::detail::airy_zero::airy_ai_zero_detail::initial_guess<T>(m);
171 // Select the maximum allowed iterations, being at least 24.
172 boost::uintmax_t number_of_iterations = (std::max)(24, int(std::numeric_limits<T>::digits10));
174 // Select the desired number of binary digits of precision.
175 // Account for the radix of number representations having non-two radix!
176 const int my_digits2 = int(float(std::numeric_limits<T>::digits)
177 * ( log(float(std::numeric_limits<T>::radix))
180 // Use a dynamic tolerance because the roots get closer the higher m gets.
183 if (m <= 10U) { tolerance = T(0.3F); }
184 else if(m <= 100U) { tolerance = T(0.1F); }
185 else if(m <= 1000U) { tolerance = T(0.05F); }
186 else { tolerance = T(1) / sqrt(T(m)); }
188 // Perform the root-finding using Newton-Raphson iteration from Boost.Math.
190 boost::math::tools::newton_raphson_iterate(
191 boost::math::detail::airy_zero::airy_ai_zero_detail::function_object_ai_and_ai_prime<T, Policy>(pol),
193 T(guess_root - tolerance),
194 T(guess_root + tolerance),
196 number_of_iterations);
198 static_cast<void>(number_of_iterations);
203 template <class T, class Policy>
204 T airy_bi_zero_imp(unsigned m, const Policy& pol)
206 BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt.
208 // Handle case when the zero'th zero is requested.
211 return policies::raise_domain_error<T>("boost::math::airy_bi_zero<%1%>(%1%,%1%)",
212 "The requested rank of the zero is %1%, but must be 1 or more !", static_cast<T>(m), pol);
214 // Set up the initial guess for the upcoming root-finding.
215 const T guess_root = boost::math::detail::airy_zero::airy_bi_zero_detail::initial_guess<T>(m);
217 // Select the maximum allowed iterations, being at least 24.
218 boost::uintmax_t number_of_iterations = (std::max)(24, int(std::numeric_limits<T>::digits10));
220 // Select the desired number of binary digits of precision.
221 // Account for the radix of number representations having non-two radix!
222 const int my_digits2 = int(float(std::numeric_limits<T>::digits)
223 * ( log(float(std::numeric_limits<T>::radix))
226 // Use a dynamic tolerance because the roots get closer the higher m gets.
229 if (m <= 10U) { tolerance = T(0.3F); }
230 else if(m <= 100U) { tolerance = T(0.1F); }
231 else if(m <= 1000U) { tolerance = T(0.05F); }
232 else { tolerance = T(1) / sqrt(T(m)); }
234 // Perform the root-finding using Newton-Raphson iteration from Boost.Math.
236 boost::math::tools::newton_raphson_iterate(
237 boost::math::detail::airy_zero::airy_bi_zero_detail::function_object_bi_and_bi_prime<T, Policy>(pol),
239 T(guess_root - tolerance),
240 T(guess_root + tolerance),
242 number_of_iterations);
244 static_cast<void>(number_of_iterations);
249 } // namespace detail
251 template <class T, class Policy>
252 inline typename tools::promote_args<T>::type airy_ai(T x, const Policy&)
254 BOOST_FPU_EXCEPTION_GUARD
255 typedef typename tools::promote_args<T>::type result_type;
256 typedef typename policies::evaluation<result_type, Policy>::type value_type;
257 typedef typename policies::normalise<
259 policies::promote_float<false>,
260 policies::promote_double<false>,
261 policies::discrete_quantile<>,
262 policies::assert_undefined<> >::type forwarding_policy;
264 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%)");
268 inline typename tools::promote_args<T>::type airy_ai(T x)
270 return airy_ai(x, policies::policy<>());
273 template <class T, class Policy>
274 inline typename tools::promote_args<T>::type airy_bi(T x, const Policy&)
276 BOOST_FPU_EXCEPTION_GUARD
277 typedef typename tools::promote_args<T>::type result_type;
278 typedef typename policies::evaluation<result_type, Policy>::type value_type;
279 typedef typename policies::normalise<
281 policies::promote_float<false>,
282 policies::promote_double<false>,
283 policies::discrete_quantile<>,
284 policies::assert_undefined<> >::type forwarding_policy;
286 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%)");
290 inline typename tools::promote_args<T>::type airy_bi(T x)
292 return airy_bi(x, policies::policy<>());
295 template <class T, class Policy>
296 inline typename tools::promote_args<T>::type airy_ai_prime(T x, const Policy&)
298 BOOST_FPU_EXCEPTION_GUARD
299 typedef typename tools::promote_args<T>::type result_type;
300 typedef typename policies::evaluation<result_type, Policy>::type value_type;
301 typedef typename policies::normalise<
303 policies::promote_float<false>,
304 policies::promote_double<false>,
305 policies::discrete_quantile<>,
306 policies::assert_undefined<> >::type forwarding_policy;
308 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%)");
312 inline typename tools::promote_args<T>::type airy_ai_prime(T x)
314 return airy_ai_prime(x, policies::policy<>());
317 template <class T, class Policy>
318 inline typename tools::promote_args<T>::type airy_bi_prime(T x, const Policy&)
320 BOOST_FPU_EXCEPTION_GUARD
321 typedef typename tools::promote_args<T>::type result_type;
322 typedef typename policies::evaluation<result_type, Policy>::type value_type;
323 typedef typename policies::normalise<
325 policies::promote_float<false>,
326 policies::promote_double<false>,
327 policies::discrete_quantile<>,
328 policies::assert_undefined<> >::type forwarding_policy;
330 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%)");
334 inline typename tools::promote_args<T>::type airy_bi_prime(T x)
336 return airy_bi_prime(x, policies::policy<>());
339 template <class T, class Policy>
340 inline T airy_ai_zero(unsigned m, const Policy& /*pol*/)
342 BOOST_FPU_EXCEPTION_GUARD
343 typedef typename policies::evaluation<T, Policy>::type value_type;
344 typedef typename policies::normalise<
346 policies::promote_float<false>,
347 policies::promote_double<false>,
348 policies::discrete_quantile<>,
349 policies::assert_undefined<> >::type forwarding_policy;
350 BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<T>::is_integer, "Airy return type must be a floating-point type.");
351 return policies::checked_narrowing_cast<T, Policy>(detail::airy_ai_zero_imp<value_type>(m, forwarding_policy()), "boost::math::airy_ai_zero<%1%>(unsigned)");
355 inline T airy_ai_zero(unsigned m)
357 return airy_ai_zero<T>(m, policies::policy<>());
360 template <class T, class OutputIterator, class Policy>
361 inline OutputIterator airy_ai_zero(
362 unsigned start_index,
363 unsigned number_of_zeros,
364 OutputIterator out_it,
367 typedef T result_type;
368 BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<result_type>::is_integer, "Airy return type must be a floating-point type.");
370 for(unsigned i = 0; i < number_of_zeros; ++i)
372 *out_it = boost::math::airy_ai_zero<result_type>(start_index + i, pol);
378 template <class T, class OutputIterator>
379 inline OutputIterator airy_ai_zero(
380 unsigned start_index,
381 unsigned number_of_zeros,
382 OutputIterator out_it)
384 return airy_ai_zero<T>(start_index, number_of_zeros, out_it, policies::policy<>());
387 template <class T, class Policy>
388 inline T airy_bi_zero(unsigned m, const Policy& /*pol*/)
390 BOOST_FPU_EXCEPTION_GUARD
391 typedef typename policies::evaluation<T, Policy>::type value_type;
392 typedef typename policies::normalise<
394 policies::promote_float<false>,
395 policies::promote_double<false>,
396 policies::discrete_quantile<>,
397 policies::assert_undefined<> >::type forwarding_policy;
398 BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<T>::is_integer, "Airy return type must be a floating-point type.");
399 return policies::checked_narrowing_cast<T, Policy>(detail::airy_bi_zero_imp<value_type>(m, forwarding_policy()), "boost::math::airy_bi_zero<%1%>(unsigned)");
402 template <typename T>
403 inline T airy_bi_zero(unsigned m)
405 return airy_bi_zero<T>(m, policies::policy<>());
408 template <class T, class OutputIterator, class Policy>
409 inline OutputIterator airy_bi_zero(
410 unsigned start_index,
411 unsigned number_of_zeros,
412 OutputIterator out_it,
415 typedef T result_type;
416 BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<result_type>::is_integer, "Airy return type must be a floating-point type.");
418 for(unsigned i = 0; i < number_of_zeros; ++i)
420 *out_it = boost::math::airy_bi_zero<result_type>(start_index + i, pol);
426 template <class T, class OutputIterator>
427 inline OutputIterator airy_bi_zero(
428 unsigned start_index,
429 unsigned number_of_zeros,
430 OutputIterator out_it)
432 return airy_bi_zero<T>(start_index, number_of_zeros, out_it, policies::policy<>());
437 #endif // BOOST_MATH_AIRY_HPP