boost/math/special_functions/detail/airy_ai_bi_zero.hpp

   1 //  Copyright (c) 2013 Christopher Kormanyos
   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 // This work is based on an earlier work:
   7 // "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
   8 // in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
   9 //
  10 // This header contains implementation details for estimating the zeros
  11 // of the Airy functions airy_ai and airy_bi on the negative real axis.
  12 //
  13 #ifndef _AIRY_AI_BI_ZERO_2013_01_20_HPP_
  14   #define _AIRY_AI_BI_ZERO_2013_01_20_HPP_
  15
  16   #include <boost/math/constants/constants.hpp>
  17   #include <boost/math/special_functions/cbrt.hpp>
  18
  19   namespace boost { namespace math {
  20   namespace detail
  21   {
  22     // Forward declarations of the needed Airy function implementations.
  23     template <class T, class Policy>
  24     T airy_ai_imp(T x, const Policy& pol);
  25     template <class T, class Policy>
  26     T airy_bi_imp(T x, const Policy& pol);
  27     template <class T, class Policy>
  28     T airy_ai_prime_imp(T x, const Policy& pol);
  29     template <class T, class Policy>
  30     T airy_bi_prime_imp(T x, const Policy& pol);
  31
  32     namespace airy_zero
  33     {
  34       template<class T>
  35       T equation_as_10_4_105(const T& z)
  36       {
  37         const T one_over_z        (T(1) / z);
  38         const T one_over_z_squared(one_over_z * one_over_z);
  39
  40         const T z_pow_third     (boost::math::cbrt(z));
  41         const T z_pow_two_thirds(z_pow_third * z_pow_third);
  42
  43         // Implement the top line of Eq. 10.4.105.
  44         const T fz(z_pow_two_thirds * (((((                     + (T(162375596875.0) / 334430208UL)
  45                                            * one_over_z_squared - (   T(108056875.0) /   6967296UL))
  46                                            * one_over_z_squared + (       T(77125UL) /     82944UL))
  47                                            * one_over_z_squared - (           T(5U)  /        36U))
  48                                            * one_over_z_squared + (           T(5U)  /        48U))
  49                                            * one_over_z_squared + (1)));
  50
  51         return fz;
  52       }
  53
  54       namespace airy_ai_zero_detail
  55       {
  56         template<class T>
  57         T initial_guess(const unsigned m)
  58         {
  59           T guess;
  60
  61           switch(m)
  62           {
  63             case 0U: { guess = T(0);                       break; }
  64             case 1U: { guess = T(-2.33810741045976703849); break; }
  65             case 2U: { guess = T(-4.08794944413097061664); break; }
  66             case 3U: { guess = T(-5.52055982809555105913); break; }
  67             case 4U: { guess = T(-6.78670809007175899878); break; }
  68             case 5U: { guess = T(-7.94413358712085312314); break; }
  69             case 6U: { guess = T(-9.02265085334098038016); break; }
  70             case 7U: { guess = T(-10.0401743415580859306); break; }
  71             case 8U: { guess = T(-11.0085243037332628932); break; }
  72             case 9U: { guess = T(-11.9360155632362625170); break; }
  73             case 10U:{ guess = T(-12.8287767528657572004); break; }
  74             default:
  75             {
  76               const T t(((boost::math::constants::pi<T>() * 3U) * ((T(m) * 4U) - 1)) / 8U);
  77               guess = -boost::math::detail::airy_zero::equation_as_10_4_105(t);
  78               break;
  79             }
  80           }
  81
  82           return guess;
  83         }
  84
  85         template<class T, class Policy>
  86         class function_object_ai_and_ai_prime
  87         {
  88         public:
  89           function_object_ai_and_ai_prime(const Policy pol) : my_pol(pol) { }
  90
  91           boost::math::tuple<T, T> operator()(const T& x) const
  92           {
  93             // Return a tuple containing both Ai(x) and Ai'(x).
  94             return boost::math::make_tuple(
  95               boost::math::detail::airy_ai_imp      (x, my_pol),
  96               boost::math::detail::airy_ai_prime_imp(x, my_pol));
  97           }
  98
  99         private:
 100           const Policy& my_pol;
 101           const function_object_ai_and_ai_prime& operator=(const function_object_ai_and_ai_prime&);
 102         };
 103       } // namespace airy_ai_zero_detail
 104
 105       namespace airy_bi_zero_detail
 106       {
 107         template<class T>
 108         T initial_guess(const unsigned m)
 109         {
 110           T guess;
 111
 112           switch(m)
 113           {
 114             case 0U: { guess = T(0);                       break; }
 115             case 1U: { guess = T(-1.17371322270912792492); break; }
 116             case 2U: { guess = T(-3.27109330283635271568); break; }
 117             case 3U: { guess = T(-4.83073784166201593267); break; }
 118             case 4U: { guess = T(-6.16985212831025125983); break; }
 119             case 5U: { guess = T(-7.37676207936776371360); break; }
 120             case 6U: { guess = T(-8.49194884650938801345); break; }
 121             case 7U: { guess = T(-9.53819437934623888663); break; }
 122             case 8U: { guess = T(-10.5299135067053579244); break; }
 123             case 9U: { guess = T(-11.4769535512787794379); break; }
 124             case 10U:{ guess = T(-12.3864171385827387456); break; }
 125             default:
 126             {
 127               const T t(((boost::math::constants::pi<T>() * 3U) * ((T(m) * 4U) - 3)) / 8U);
 128               guess = -boost::math::detail::airy_zero::equation_as_10_4_105(t);
 129               break;
 130             }
 131           }
 132
 133           return guess;
 134         }
 135
 136         template<class T, class Policy>
 137         class function_object_bi_and_bi_prime
 138         {
 139         public:
 140           function_object_bi_and_bi_prime(const Policy pol) : my_pol(pol) { }
 141
 142           boost::math::tuple<T, T> operator()(const T& x) const
 143           {
 144             // Return a tuple containing both Bi(x) and Bi'(x).
 145             return boost::math::make_tuple(
 146               boost::math::detail::airy_bi_imp      (x, my_pol),
 147               boost::math::detail::airy_bi_prime_imp(x, my_pol));
 148           }
 149
 150         private:
 151           const Policy& my_pol;
 152           const function_object_bi_and_bi_prime& operator=(const function_object_bi_and_bi_prime&);
 153         };
 154       } // namespace airy_bi_zero_detail
 155     } // namespace airy_zero
 156   } // namespace detail
 157   } // namespace math
 158   } // namespaces boost
 159
 160 #endif // _AIRY_AI_BI_ZERO_2013_01_20_HPP_