boost/math/special_functions/airy.hpp

   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 <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>
  14
  15 namespace boost{ namespace math{
  16
  17 namespace detail{
  18
  19 template <class T, class Policy>
  20 T airy_ai_imp(T x, const Policy& pol)
  21 {
  22    BOOST_MATH_STD_USING
  23
  24    if(x < 0)
  25    {
  26       T p = (-x * sqrt(-x) * 2) / 3;
  27       T v = T(1) / 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);
  32       return ai;
  33    }
  34    else if(fabs(x * x * x) / 6 < tools::epsilon<T>())
  35    {
  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);
  39       return ai;
  40    }
  41    else
  42    {
  43       T p = 2 * x * sqrt(x) / 3;
  44       T v = T(1) / 3;
  45       //T j1 = boost::math::cyl_bessel_i(-v, p, pol);
  46       //T j2 = boost::math::cyl_bessel_i(v, p, pol);
  47       //
  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:
  50       //
  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);
  53       return ai;
  54    }
  55 }
  56
  57 template <class T, class Policy>
  58 T airy_bi_imp(T x, const Policy& pol)
  59 {
  60    BOOST_MATH_STD_USING
  61
  62    if(x < 0)
  63    {
  64       T p = (-x * sqrt(-x) * 2) / 3;
  65       T v = T(1) / 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);
  70       return bi;
  71    }
  72    else if(fabs(x * x * x) / 6 < tools::epsilon<T>())
  73    {
  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);
  77       return bi;
  78    }
  79    else
  80    {
  81       T p = 2 * x * sqrt(x) / 3;
  82       T v = T(1) / 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);
  86       return bi;
  87    }
  88 }
  89
  90 template <class T, class Policy>
  91 T airy_ai_prime_imp(T x, const Policy& pol)
  92 {
  93    BOOST_MATH_STD_USING
  94
  95    if(x < 0)
  96    {
  97       T p = (-x * sqrt(-x) * 2) / 3;
  98       T v = T(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;
 102       return aip;
 103    }
 104    else if(fabs(x * x) / 2 < tools::epsilon<T>())
 105    {
 106       T tg = boost::math::tgamma(constants::third<T>(), pol);
 107       T aip = 1 / (boost::math::cbrt(T(3)) * tg);
 108       return -aip;
 109    }
 110    else
 111    {
 112       T p = 2 * x * sqrt(x) / 3;
 113       T v = T(2) / 3;
 114       //T j1 = boost::math::cyl_bessel_i(-v, p, pol);
 115       //T j2 = boost::math::cyl_bessel_i(v, p, pol);
 116       //
 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:
 119       //
 120       T aip = -cyl_bessel_k(v, p, pol) * x / (boost::math::constants::root_three<T>() * boost::math::constants::pi<T>());
 121       return aip;
 122    }
 123 }
 124
 125 template <class T, class Policy>
 126 T airy_bi_prime_imp(T x, const Policy& pol)
 127 {
 128    BOOST_MATH_STD_USING
 129
 130    if(x < 0)
 131    {
 132       T p = (-x * sqrt(-x) * 2) / 3;
 133       T v = T(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>();
 137       return aip;
 138    }
 139    else if(fabs(x * x) / 2 < tools::epsilon<T>())
 140    {
 141       T tg = boost::math::tgamma(constants::third<T>(), pol);
 142       T bip = sqrt(boost::math::cbrt(T(3))) / tg;
 143       return bip;
 144    }
 145    else
 146    {
 147       T p = 2 * x * sqrt(x) / 3;
 148       T v = T(2) / 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>();
 152       return aip;
 153    }
 154 }
 155
 156 template <class T, class Policy>
 157 T airy_ai_zero_imp(unsigned m, const Policy& pol)
 158 {
 159    BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt.
 160
 161    // Handle case when the zero'th zero is requested.
 162    if(m == 0U)
 163    {
 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);
 166    }
 167
 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);
 170
 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));
 173
 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))
 178                                  / log(2.0F)));
 179
 180    // Use a dynamic tolerance because the roots get closer the higher m gets.
 181    T tolerance;
 182
 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)); }
 187
 188    // Perform the root-finding using Newton-Raphson iteration from Boost.Math.
 189    const T am =
 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),
 192          guess_root,
 193          T(guess_root - tolerance),
 194          T(guess_root + tolerance),
 195          my_digits2,
 196          number_of_iterations);
 197
 198    static_cast<void>(number_of_iterations);
 199
 200    return am;
 201 }
 202
 203 template <class T, class Policy>
 204 T airy_bi_zero_imp(unsigned m, const Policy& pol)
 205 {
 206    BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt.
 207
 208    // Handle case when the zero'th zero is requested.
 209    if(m == 0U)
 210    {
 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);
 213    }
 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);
 216
 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));
 219
 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))
 224                                  / log(2.0F)));
 225
 226    // Use a dynamic tolerance because the roots get closer the higher m gets.
 227    T tolerance;
 228
 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)); }
 233
 234    // Perform the root-finding using Newton-Raphson iteration from Boost.Math.
 235    const T bm =
 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),
 238          guess_root,
 239          T(guess_root - tolerance),
 240          T(guess_root + tolerance),
 241          my_digits2,
 242          number_of_iterations);
 243
 244    static_cast<void>(number_of_iterations);
 245
 246    return bm;
 247 }
 248
 249 } // namespace detail
 250
 251 template <class T, class Policy>
 252 inline typename tools::promote_args<T>::type airy_ai(T x, const Policy&)
 253 {
 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<
 258       Policy,
 259       policies::promote_float<false>,
 260       policies::promote_double<false>,
 261       policies::discrete_quantile<>,
 262       policies::assert_undefined<> >::type forwarding_policy;
 263
 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%)");
 265 }
 266
 267 template <class T>
 268 inline typename tools::promote_args<T>::type airy_ai(T x)
 269 {
 270    return airy_ai(x, policies::policy<>());
 271 }
 272
 273 template <class T, class Policy>
 274 inline typename tools::promote_args<T>::type airy_bi(T x, const Policy&)
 275 {
 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<
 280       Policy,
 281       policies::promote_float<false>,
 282       policies::promote_double<false>,
 283       policies::discrete_quantile<>,
 284       policies::assert_undefined<> >::type forwarding_policy;
 285
 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%)");
 287 }
 288
 289 template <class T>
 290 inline typename tools::promote_args<T>::type airy_bi(T x)
 291 {
 292    return airy_bi(x, policies::policy<>());
 293 }
 294
 295 template <class T, class Policy>
 296 inline typename tools::promote_args<T>::type airy_ai_prime(T x, const Policy&)
 297 {
 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<
 302       Policy,
 303       policies::promote_float<false>,
 304       policies::promote_double<false>,
 305       policies::discrete_quantile<>,
 306       policies::assert_undefined<> >::type forwarding_policy;
 307
 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%)");
 309 }
 310
 311 template <class T>
 312 inline typename tools::promote_args<T>::type airy_ai_prime(T x)
 313 {
 314    return airy_ai_prime(x, policies::policy<>());
 315 }
 316
 317 template <class T, class Policy>
 318 inline typename tools::promote_args<T>::type airy_bi_prime(T x, const Policy&)
 319 {
 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<
 324       Policy,
 325       policies::promote_float<false>,
 326       policies::promote_double<false>,
 327       policies::discrete_quantile<>,
 328       policies::assert_undefined<> >::type forwarding_policy;
 329
 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%)");
 331 }
 332
 333 template <class T>
 334 inline typename tools::promote_args<T>::type airy_bi_prime(T x)
 335 {
 336    return airy_bi_prime(x, policies::policy<>());
 337 }
 338
 339 template <class T, class Policy>
 340 inline T airy_ai_zero(unsigned m, const Policy& /*pol*/)
 341 {
 342    BOOST_FPU_EXCEPTION_GUARD
 343    typedef typename policies::evaluation<T, Policy>::type value_type;
 344    typedef typename policies::normalise<
 345       Policy,
 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)");
 352 }
 353
 354 template <class T>
 355 inline T airy_ai_zero(unsigned m)
 356 {
 357    return airy_ai_zero<T>(m, policies::policy<>());
 358 }
 359
 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,
 365                          const Policy& pol)
 366 {
 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.");
 369
 370    for(unsigned i = 0; i < number_of_zeros; ++i)
 371    {
 372       *out_it = boost::math::airy_ai_zero<result_type>(start_index + i, pol);
 373       ++out_it;
 374    }
 375    return out_it;
 376 }
 377
 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)
 383 {
 384    return airy_ai_zero<T>(start_index, number_of_zeros, out_it, policies::policy<>());
 385 }
 386
 387 template <class T, class Policy>
 388 inline T airy_bi_zero(unsigned m, const Policy& /*pol*/)
 389 {
 390    BOOST_FPU_EXCEPTION_GUARD
 391    typedef typename policies::evaluation<T, Policy>::type value_type;
 392    typedef typename policies::normalise<
 393       Policy,
 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)");
 400 }
 401
 402 template <typename T>
 403 inline T airy_bi_zero(unsigned m)
 404 {
 405    return airy_bi_zero<T>(m, policies::policy<>());
 406 }
 407
 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,
 413                          const Policy& pol)
 414 {
 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.");
 417
 418    for(unsigned i = 0; i < number_of_zeros; ++i)
 419    {
 420       *out_it = boost::math::airy_bi_zero<result_type>(start_index + i, pol);
 421       ++out_it;
 422    }
 423    return out_it;
 424 }
 425
 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)
 431 {
 432    return airy_bi_zero<T>(start_index, number_of_zeros, out_it, policies::policy<>());
 433 }
 434
 435 }} // namespaces
 436
 437 #endif // BOOST_MATH_AIRY_HPP