123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521 |
- // Copyright Nick Thompson, 2017
- // Copyright John Maddock 2017
- // Use, modification and distribution are subject to the
- // Boost Software License, Version 1.0.
- // (See accompanying file LICENSE_1_0.txt
- // or copy at http://www.boost.org/LICENSE_1_0.txt)
- #include <cmath>
- #include <cstdint>
- #include <functional>
- #include <iomanip>
- #include <iostream>
- #include <numeric>
- #include <boost/math/constants/constants.hpp>
- #include <boost/math/special_functions/cbrt.hpp>
- #include <boost/math/special_functions/factorials.hpp>
- #include <boost/math/special_functions/gamma.hpp>
- #include <boost/math/tools/roots.hpp>
- #include <boost/noncopyable.hpp>
- #define CPP_BIN_FLOAT 1
- #define CPP_DEC_FLOAT 2
- #define CPP_MPFR_FLOAT 3
- //#define MP_TYPE CPP_BIN_FLOAT
- #define MP_TYPE CPP_DEC_FLOAT
- //#define MP_TYPE CPP_MPFR_FLOAT
- namespace
- {
- struct digits_characteristics
- {
- static const int digits10 = 300;
- static const int guard_digits = 6;
- };
- }
- #if (MP_TYPE == CPP_BIN_FLOAT)
- #include <boost/multiprecision/cpp_bin_float.hpp>
- namespace mp = boost::multiprecision;
- typedef mp::number<mp::cpp_bin_float<digits_characteristics::digits10 + digits_characteristics::guard_digits>, mp::et_off> mp_type;
- #elif (MP_TYPE == CPP_DEC_FLOAT)
- #include <boost/multiprecision/cpp_dec_float.hpp>
- namespace mp = boost::multiprecision;
- typedef mp::number<mp::cpp_dec_float<digits_characteristics::digits10 + digits_characteristics::guard_digits>, mp::et_off> mp_type;
- #elif (MP_TYPE == CPP_MPFR_FLOAT)
- #include <boost/multiprecision/mpfr.hpp>
- namespace mp = boost::multiprecision;
- typedef mp::number<mp::mpfr_float_backend<digits_characteristics::digits10 + digits_characteristics::guard_digits>, mp::et_off> mp_type;
- #else
- #error MP_TYPE is undefined
- #endif
- template<typename T>
- class laguerre_function_object
- {
- public:
- laguerre_function_object(const int n, const T a) : order(n),
- alpha(a),
- p1 (0),
- d2 (0) { }
- laguerre_function_object(const laguerre_function_object& other) : order(other.order),
- alpha(other.alpha),
- p1 (other.p1),
- d2 (other.d2) { }
- ~laguerre_function_object() { }
- T operator()(const T& x) const
- {
- // Calculate (via forward recursion):
- // * the value of the Laguerre function L(n, alpha, x), called (p2),
- // * the value of the derivative of the Laguerre function (d2),
- // * and the value of the corresponding Laguerre function of
- // previous order (p1).
- // Return the value of the function (p2) in order to be used as a
- // function object with Boost.Math root-finding. Store the values
- // of the Laguerre function derivative (d2) and the Laguerre function
- // of previous order (p1) in class members for later use.
- p1 = T(0);
- T p2 = T(1);
- d2 = T(0);
- T j_plus_alpha(alpha);
- T two_j_plus_one_plus_alpha_minus_x(1 + alpha - x);
- int j;
- const T my_two(2);
- for(j = 0; j < order; ++j)
- {
- const T p0(p1);
- // Set the value of the previous Laguerre function.
- p1 = p2;
- // Use a recurrence relation to compute the value of the Laguerre function.
- p2 = ((two_j_plus_one_plus_alpha_minus_x * p1) - (j_plus_alpha * p0)) / (j + 1);
- ++j_plus_alpha;
- two_j_plus_one_plus_alpha_minus_x += my_two;
- }
- // Set the value of the derivative of the Laguerre function.
- d2 = ((p2 * j) - (j_plus_alpha * p1)) / x;
- // Return the value of the Laguerre function.
- return p2;
- }
- const T& previous () const { return p1; }
- const T& derivative() const { return d2; }
- static bool root_tolerance(const T& a, const T& b)
- {
- using std::abs;
- // The relative tolerance here is: ((a - b) * 2) / (a + b).
- return (abs((a - b) * 2) < ((a + b) * boost::math::tools::epsilon<T>()));
- }
- private:
- const int order;
- const T alpha;
- mutable T p1;
- mutable T d2;
- laguerre_function_object();
- const laguerre_function_object& operator=(const laguerre_function_object&);
- };
- template<typename T>
- class guass_laguerre_abscissas_and_weights : private boost::noncopyable
- {
- public:
- guass_laguerre_abscissas_and_weights(const int n, const T a) : order(n),
- alpha(a),
- valid(true),
- xi (),
- wi ()
- {
- if(alpha < -20.0F)
- {
- // TBD: If we ever boostify this, throw a range error here.
- // If so, then also document it in the docs.
- std::cout << "Range error: the order of the Laguerre function must exceed -20.0." << std::endl;
- }
- else
- {
- calculate();
- }
- }
- virtual ~guass_laguerre_abscissas_and_weights() { }
- const std::vector<T>& abscissas() const { return xi; }
- const std::vector<T>& weights () const { return wi; }
- bool get_valid() const { return valid; }
- private:
- const int order;
- const T alpha;
- bool valid;
- std::vector<T> xi;
- std::vector<T> wi;
- void calculate()
- {
- using std::abs;
- std::cout << "finding approximate roots..." << std::endl;
- std::vector<boost::math::tuple<T, T> > root_estimates;
- root_estimates.reserve(static_cast<typename std::vector<boost::math::tuple<T, T> >::size_type>(order));
- const laguerre_function_object<T> laguerre_object(order, alpha);
- // Set the initial values of the step size and the running step
- // to be used for finding the estimate of the first root.
- T step_size = 0.01F;
- T step = step_size;
- T first_laguerre_root = 0.0F;
- bool first_laguerre_root_has_been_found = true;
- if(alpha < -1.0F)
- {
- // Iteratively step through the Laguerre function using a
- // small step-size in order to find a rough estimate of
- // the first zero.
- bool this_laguerre_value_is_negative = (laguerre_object(mp_type(0)) < 0);
- static const int j_max = 10000;
- int j;
- for(j = 0; (j < j_max) && (this_laguerre_value_is_negative != (laguerre_object(step) < 0)); ++j)
- {
- // Increment the step size until the sign of the Laguerre function
- // switches. This indicates a zero-crossing, signalling the next root.
- step += step_size;
- }
- if(j >= j_max)
- {
- first_laguerre_root_has_been_found = false;
- }
- else
- {
- // We have found the first zero-crossing. Put a loose bracket around
- // the root using a window. Here, we know that the first root lies
- // between (x - step_size) < root < x.
- // Before storing the approximate root, perform a couple of
- // bisection steps in order to tighten up the root bracket.
- boost::uintmax_t a_couple_of_iterations = 3U;
- const std::pair<T, T>
- first_laguerre_root = boost::math::tools::bisect(laguerre_object,
- step - step_size,
- step,
- laguerre_function_object<T>::root_tolerance,
- a_couple_of_iterations);
- static_cast<void>(a_couple_of_iterations);
- }
- }
- else
- {
- // Calculate an estimate of the 1st root of a generalized Laguerre
- // function using either a Taylor series or an expansion in Bessel
- // function zeros. The Bessel function zeros expansion is from Tricomi.
- // Here, we obtain an estimate of the first zero of J_alpha(x).
- T j_alpha_m1;
- if(alpha < 1.4F)
- {
- // For small alpha, use a short series obtained from Mathematica(R).
- // Series[BesselJZero[v, 1], {v, 0, 3}]
- // N[%, 12]
- j_alpha_m1 = ((( 0.09748661784476F
- * alpha - 0.17549359276115F)
- * alpha + 1.54288974259931F)
- * alpha + 2.40482555769577F);
- }
- else
- {
- // For larger alpha, use the first line of Eqs. 10.21.40 in the NIST Handbook.
- const T alpha_pow_third(boost::math::cbrt(alpha));
- const T alpha_pow_minus_two_thirds(T(1) / (alpha_pow_third * alpha_pow_third));
- j_alpha_m1 = alpha * ((((( + 0.043F
- * alpha_pow_minus_two_thirds - 0.0908F)
- * alpha_pow_minus_two_thirds - 0.00397F)
- * alpha_pow_minus_two_thirds + 1.033150F)
- * alpha_pow_minus_two_thirds + 1.8557571F)
- * alpha_pow_minus_two_thirds + 1.0F);
- }
- const T vf = ((order * 4.0F) + (alpha * 2.0F) + 2.0F);
- const T vf2 = vf * vf;
- const T j_alpha_m1_sqr = j_alpha_m1 * j_alpha_m1;
- first_laguerre_root = (j_alpha_m1_sqr * (-0.6666666666667F + ((0.6666666666667F * alpha) * alpha) + (0.3333333333333F * j_alpha_m1_sqr) + vf2)) / (vf2 * vf);
- }
- if(first_laguerre_root_has_been_found)
- {
- bool this_laguerre_value_is_negative = (laguerre_object(mp_type(0)) < 0);
- // Re-set the initial value of the step-size based on the
- // estimate of the first root.
- step_size = first_laguerre_root / 2;
- step = step_size;
- // Step through the Laguerre function using a step-size
- // of dynamic width in order to find the zero crossings
- // of the Laguerre function, providing rough estimates
- // of the roots. Refine the brackets with a few bisection
- // steps, and store the results as bracketed root estimates.
- while(static_cast<int>(root_estimates.size()) < order)
- {
- // Increment the step size until the sign of the Laguerre function
- // switches. This indicates a zero-crossing, signalling the next root.
- step += step_size;
- if(this_laguerre_value_is_negative != (laguerre_object(step) < 0))
- {
- // We have found the next zero-crossing.
- // Change the running sign of the Laguerre function.
- this_laguerre_value_is_negative = (!this_laguerre_value_is_negative);
- // We have found the first zero-crossing. Put a loose bracket around
- // the root using a window. Here, we know that the first root lies
- // between (x - step_size) < root < x.
- // Before storing the approximate root, perform a couple of
- // bisection steps in order to tighten up the root bracket.
- boost::uintmax_t a_couple_of_iterations = 3U;
- const std::pair<T, T>
- root_estimate_bracket = boost::math::tools::bisect(laguerre_object,
- step - step_size,
- step,
- laguerre_function_object<T>::root_tolerance,
- a_couple_of_iterations);
- static_cast<void>(a_couple_of_iterations);
- // Store the refined root estimate as a bracketed range in a tuple.
- root_estimates.push_back(boost::math::tuple<T, T>(root_estimate_bracket.first,
- root_estimate_bracket.second));
- if(root_estimates.size() >= static_cast<std::size_t>(2U))
- {
- // Determine the next step size. This is based on the distance between
- // the previous two roots, whereby the estimates of the previous roots
- // are computed by taking the average of the lower and upper range of
- // the root-estimate bracket.
- const T r0 = ( boost::math::get<0>(*(root_estimates.rbegin() + 1U))
- + boost::math::get<1>(*(root_estimates.rbegin() + 1U))) / 2;
- const T r1 = ( boost::math::get<0>(*root_estimates.rbegin())
- + boost::math::get<1>(*root_estimates.rbegin())) / 2;
- const T distance_between_previous_roots = r1 - r0;
- step_size = distance_between_previous_roots / 3;
- }
- }
- }
- const T norm_g =
- ((alpha == 0) ? T(-1)
- : -boost::math::tgamma(alpha + order) / boost::math::factorial<T>(order - 1));
- xi.reserve(root_estimates.size());
- wi.reserve(root_estimates.size());
- // Calculate the abscissas and weights to full precision.
- for(std::size_t i = static_cast<std::size_t>(0U); i < root_estimates.size(); ++i)
- {
- std::cout << "calculating abscissa and weight for index: " << i << std::endl;
- // Calculate the abscissas using iterative root-finding.
- // Select the maximum allowed iterations, being at least 20.
- // The determination of the maximum allowed iterations is
- // based on the number of decimal digits in the numerical
- // type T.
- const int my_digits10 = static_cast<int>(static_cast<float>(boost::math::tools::digits<T>()) * 0.301F);
- const boost::uintmax_t number_of_iterations_allowed = (std::max)(20, my_digits10 / 2);
- boost::uintmax_t number_of_iterations_used = number_of_iterations_allowed;
- // Perform the root-finding using ACM TOMS 748 from Boost.Math.
- const std::pair<T, T>
- laguerre_root_bracket = boost::math::tools::toms748_solve(laguerre_object,
- boost::math::get<0>(root_estimates[i]),
- boost::math::get<1>(root_estimates[i]),
- laguerre_function_object<T>::root_tolerance,
- number_of_iterations_used);
- // Based on the result of *each* root-finding operation, re-assess
- // the validity of the Guass-Laguerre abscissas and weights object.
- valid &= (number_of_iterations_used < number_of_iterations_allowed);
- // Compute the Laguerre root as the average of the values from
- // the solved root bracket.
- const T laguerre_root = ( laguerre_root_bracket.first
- + laguerre_root_bracket.second) / 2;
- // Calculate the weight for this Laguerre root. Here, we calculate
- // the derivative of the Laguerre function and the value of the
- // previous Laguerre function on the x-axis at the value of this
- // Laguerre root.
- static_cast<void>(laguerre_object(laguerre_root));
- // Store the abscissa and weight for this index.
- xi.push_back(laguerre_root);
- wi.push_back(norm_g / ((laguerre_object.derivative() * order) * laguerre_object.previous()));
- }
- }
- }
- };
- namespace
- {
- template<typename T>
- struct gauss_laguerre_ai
- {
- public:
- gauss_laguerre_ai(const T X) : x(X)
- {
- using std::exp;
- using std::sqrt;
- zeta = ((sqrt(x) * x) * 2) / 3;
- const T zeta_times_48_pow_sixth = sqrt(boost::math::cbrt(zeta * 48));
- factor = 1 / ((sqrt(boost::math::constants::pi<T>()) * zeta_times_48_pow_sixth) * (exp(zeta) * gamma_of_five_sixths()));
- }
- gauss_laguerre_ai(const gauss_laguerre_ai& other) : x (other.x),
- zeta (other.zeta),
- factor(other.factor) { }
- T operator()(const T& t) const
- {
- using std::sqrt;
- return factor / sqrt(boost::math::cbrt(2 + (t / zeta)));
- }
- private:
- const T x;
- T zeta;
- T factor;
- static const T& gamma_of_five_sixths()
- {
- static const T value = boost::math::tgamma(T(5) / 6);
- return value;
- }
- const gauss_laguerre_ai& operator=(const gauss_laguerre_ai&);
- };
- template<typename T>
- T gauss_laguerre_airy_ai(const T x)
- {
- static const float digits_factor = static_cast<float>(std::numeric_limits<mp_type>::digits10) / 300.0F;
- static const int laguerre_order = static_cast<int>(600.0F * digits_factor);
- static const guass_laguerre_abscissas_and_weights<T> abscissas_and_weights(laguerre_order, -T(1) / 6);
- T airy_ai_result;
- if(abscissas_and_weights.get_valid())
- {
- const gauss_laguerre_ai<T> this_gauss_laguerre_ai(x);
- airy_ai_result =
- std::inner_product(abscissas_and_weights.abscissas().begin(),
- abscissas_and_weights.abscissas().end(),
- abscissas_and_weights.weights().begin(),
- T(0),
- std::plus<T>(),
- [&this_gauss_laguerre_ai](const T& this_abscissa, const T& this_weight) -> T
- {
- return this_gauss_laguerre_ai(this_abscissa) * this_weight;
- });
- }
- else
- {
- // TBD: Consider an error message.
- airy_ai_result = T(0);
- }
- return airy_ai_result;
- }
- }
- int main()
- {
- // Use Gauss-Laguerre integration to compute airy_ai(120 / 7).
- // 9 digits
- // 3.89904210e-22
- // 10 digits
- // 3.899042098e-22
- // 50 digits.
- // 3.8990420982303275013276114626640705170145070824318e-22
- // 100 digits.
- // 3.899042098230327501327611462664070517014507082431797677146153303523108862015228
- // 864136051942933142648e-22
- // 200 digits.
- // 3.899042098230327501327611462664070517014507082431797677146153303523108862015228
- // 86413605194293314264788265460938200890998546786740097437064263800719644346113699
- // 77010905030516409847054404055843899790277e-22
- // 300 digits.
- // 3.899042098230327501327611462664070517014507082431797677146153303523108862015228
- // 86413605194293314264788265460938200890998546786740097437064263800719644346113699
- // 77010905030516409847054404055843899790277083960877617919088116211775232728792242
- // 9346416823281460245814808276654088201413901972239996130752528e-22
- // 500 digits.
- // 3.899042098230327501327611462664070517014507082431797677146153303523108862015228
- // 86413605194293314264788265460938200890998546786740097437064263800719644346113699
- // 77010905030516409847054404055843899790277083960877617919088116211775232728792242
- // 93464168232814602458148082766540882014139019722399961307525276722937464859521685
- // 42826483602153339361960948844649799257455597165900957281659632186012043089610827
- // 78871305322190941528281744734605934497977375094921646511687434038062987482900167
- // 45127557400365419545e-22
- // Mathematica(R) or Wolfram's Alpha:
- // N[AiryAi[120 / 7], 300]
- std::cout << std::setprecision(digits_characteristics::digits10)
- << gauss_laguerre_airy_ai(mp_type(120) / 7)
- << std::endl;
- }
|