multiprecision/example/gauss_laguerre_quadrature.cpp

///////////////////////////////////////////////////////////////////////////////
//      Copyright Christopher Kormanyos 2012 - 2015, 2020.
// Distributed under 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)
//

// This example uses Boost.Multiprecision to implement
// a high-precision Gauss-Laguerre quadrature integration.
// The quadrature is used to calculate the airy_ai(x) function
// for real-valued x on the positive axis with x.ge.1.

// In this way, the integral representation could be seen
// as part of a scheme to calculate real-valued Airy functions
// on the positive axis for medium to large argument.
// A Taylor series or hypergeometric function (not part
// of this example) could be used for smaller arguments.

// This example has been tested with decimal digits counts
// ranging from 21...301, by adjusting the parameter
// local::my_digits10 at compile time.

// The quadrature integral representaion of airy_ai(x) used
// in this example can be found in:

// A. Gil, J. Segura, N.M. Temme, "Numerical Methods for Special
// Functions" (SIAM Press 2007), Sect. 5.3.3, in particular Eq. 5.110,
// page 145. Subsequently, Gil et al's book cites the another work:
// W. Gautschi, "Computation of Bessel and Airy functions and of
// related Gaussian quadrature formulae", BIT, 42 (2002), pp. 110-118.

#include <cmath>
#include <cstdint>
#include <functional>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <sstream>
#include <tuple>
#include <vector>

#include <boost/cstdfloat.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/cbrt.hpp>
#include <boost/math/special_functions/bessel.hpp>
#include <boost/math/special_functions/factorials.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/tools/roots.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/noncopyable.hpp>

namespace gauss { namespace laguerre {

namespace util {

void progress_bar(std::ostream& os, const float percent);

void progress_bar(std::ostream& os, const float percent)
{
  std::stringstream strstrm;

  strstrm.precision(1);

  strstrm << std::fixed << percent << "%";

  os << strstrm.str() << "\r";

  os.flush();
}

}

namespace detail
{
  template<typename T>
  class laguerre_l_object BOOST_FINAL
  {
  public:
    laguerre_l_object(const int n, const T a) BOOST_NOEXCEPT
      : order(n),
        alpha(a),
        p1   (0),
        d2   (0) { }

    laguerre_l_object& operator=(const laguerre_l_object& other)
    {
      if(this != other)
      {
        order = other.order;
        alpha = other.alpha;
        p1    = other.p1;
        d2    = other.d2;
      }

      return *this;
    }

    T operator()(const T& x) const BOOST_NOEXCEPT
    {
      // 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;

      const T my_two = 2;

      for(int 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 * order) - (j_plus_alpha * p1)) / x;

      // Return the value of the Laguerre function.
      return p2;
    }

    const T previous  () const BOOST_NOEXCEPT { return p1; }
    const T derivative() const BOOST_NOEXCEPT { return d2; }

    static bool root_tolerance(const T& a, const T& b) BOOST_NOEXCEPT
    {
      using std::fabs;

      // The relative tolerance here is: ((a - b) * 2) / (a + b).
      return ((fabs(a - b) * 2) < (fabs(a + b) * std::numeric_limits<T>::epsilon()));
    }

  private:
    const   int order;
    const   T   alpha;
    mutable T   p1;
    mutable T   d2;
  };

  template<typename T>
  class abscissas_and_weights : private boost::noncopyable
  {
  public:
    abscissas_and_weights(const int n, const T a) : order(n),
                                                    alpha(a),
                                                    xi   (),
                                                    wi   ()
    {
      if(alpha < -20.0F)
      {
        // Users can decide to perform different range checking.
        std::cout << "Range error: the order of the Laguerre function must exceed -20.0."
                  << std::endl;
      }
      else
      {
        calculate();
      }
    }

    const std::vector<T>& abscissa_n() const BOOST_NOEXCEPT { return xi; }
    const std::vector<T>& weight_n  () const BOOST_NOEXCEPT { return wi; }

  private:
    const int order;
    const T   alpha;

    std::vector<T> xi;
    std::vector<T> wi;

    abscissas_and_weights() : order(),
                              alpha(),
                              xi   (),
                              wi   () { }

    void calculate()
    {
      using std::abs;

      std::cout << "Finding the approximate roots..." << std::endl;

      std::vector<std::tuple<T, T>> root_estimates;

      root_estimates.reserve(static_cast<typename std::vector<std::tuple<T, T>>::size_type>(order));

      const laguerre_l_object<T> laguerre_root_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;

      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.

        const bool this_laguerre_value_is_negative = (laguerre_root_object(T(0)) < 0);

        BOOST_CONSTEXPR_OR_CONST int j_max = 10000;

        int j = 0;

        while((j < j_max) && (this_laguerre_value_is_negative != (laguerre_root_object(step) < 0)))
        {
          // Increment the step size until the sign of the Laguerre function
          // switches. This indicates a zero-crossing, signalling the next root.
          step += step_size;

          ++j;
        }
      }
      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 cyl_bessel_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 = (((          T(0.09748661784476F)
                          * alpha - T(0.17549359276115F))
                          * alpha + T(1.54288974259931F))
                          * alpha + T(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 * (((((                             + T(0.043F)
                                    * alpha_pow_minus_two_thirds - T(0.0908F))
                                    * alpha_pow_minus_two_thirds - T(0.00397F))
                                    * alpha_pow_minus_two_thirds + T(1.033150F))
                                    * alpha_pow_minus_two_thirds + T(1.8557571F))
                                    * alpha_pow_minus_two_thirds + T(1.0F));
        }

        const T vf             = (  T(order * 4U)
                                  + T(alpha * 2U)
                                  + T(2U));
        const T vf2            = vf * vf;
        const T j_alpha_m1_sqr = j_alpha_m1 * j_alpha_m1;

        first_laguerre_root = (j_alpha_m1_sqr * (   -T(0.6666666666667F)
                                                 + ((T(0.6666666666667F) * alpha) * alpha)
                                                 +  (T(0.3333333333333F) * j_alpha_m1_sqr) + vf2)) / (vf2 * vf);
      }

      bool this_laguerre_value_is_negative = (laguerre_root_object(T(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(root_estimates.size() < static_cast<std::size_t>(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_root_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 = 4U;

          const std::pair<T, T>
            root_estimate_bracket = boost::math::tools::bisect(laguerre_root_object,
                                                               step - step_size,
                                                               step,
                                                               laguerre_l_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(std::tuple<T, T>(std::get<0>(root_estimate_bracket),
                                                    std::get<1>(root_estimate_bracket)));

          if(    (root_estimates.size() == 1U)
             || ((root_estimates.size() % 8U) == 0U)
             ||  (root_estimates.size() == static_cast<std::size_t>(order)))
          {
            const float progress = (100.0F * static_cast<float>(root_estimates.size())) / static_cast<float>(order);

            std::cout << "root_estimates.size(): "
                      << root_estimates.size()
                      << ", "
                      ;

            util::progress_bar(std::cout, progress);
          }

          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 = (  std::get<0>(*(root_estimates.crbegin() + 1U))
                          + std::get<1>(*(root_estimates.crbegin() + 1U))) / 2;

            const T r1 = (  std::get<0>(*root_estimates.crbegin())
                          + std::get<1>(*root_estimates.crbegin())) / 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());

      std::cout << std::endl;

      // Calculate the abscissas and weights to full precision.
      for(std::size_t i = static_cast<std::size_t>(0U); i < root_estimates.size(); ++i)
      {
        if(   ((i % 8U) == 0U)
           || ( i == root_estimates.size() - 1U))
        {
          const float progress = (100.0F * static_cast<float>(i + 1U)) / static_cast<float>(root_estimates.size());

          std::cout << "Calculating abscissas and weights. Processed "
                    << (i + 1U)
                    << ", "
                    ;

          util::progress_bar(std::cout, progress);
        }

        // 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.
        BOOST_CONSTEXPR_OR_CONST int local_math_tools_digits10 =
          static_cast<int>(static_cast<boost::float_least32_t>(boost::math::tools::digits<T>()) * BOOST_FLOAT32_C(0.301));

        const boost::uintmax_t number_of_iterations_allowed = (std::max)(20, local_math_tools_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_root_object,
                                                                    std::get<0>(root_estimates.at(i)),
                                                                    std::get<1>(root_estimates.at(i)),
                                                                    laguerre_l_object<T>::root_tolerance,
                                                                    number_of_iterations_used);

        static_cast<void>(number_of_iterations_used);

        // Compute the Laguerre root as the average of the values from
        // the solved root bracket.
        const T laguerre_root = (  std::get<0>(laguerre_root_bracket)
                                 + std::get<1>(laguerre_root_bracket)) / 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_root_object(laguerre_root));

        // Store the abscissa and weight for this index.
        xi.push_back(laguerre_root);
        wi.push_back(norm_g / ((laguerre_root_object.derivative() * order) * laguerre_root_object.previous()));
      }

      std::cout << std::endl;
    }
  };

  template<typename T>
  struct airy_ai_object BOOST_FINAL
  {
  public:
    airy_ai_object(const T& x)
      : my_x     (x),
        my_zeta  (((sqrt(x) * x) * 2) / 3),
        my_factor(make_factor(my_zeta)) { }

    T operator()(const T& t) const
    {
      using std::sqrt;

      return my_factor / sqrt(boost::math::cbrt(2 + (t / my_zeta)));
    }

  private:
    const T my_x;
    const T my_zeta;
    const T my_factor;

    airy_ai_object() : my_x     (),
                       my_zeta  (),
                       my_factor() { }

    static T make_factor(const T& z)
    {
      using std::exp;
      using std::sqrt;

      return 1 / ((sqrt(boost::math::constants::pi<T>()) * sqrt(boost::math::cbrt(z * 48))) * (exp(z) * boost::math::tgamma(T(5) / 6)));
    }
  };
} // namespace detail

} } // namespace gauss::laguerre


// A float_type is created to handle the desired number of decimal digits from `cpp_dec_float` without using __expression_templates.
struct local
{
  BOOST_STATIC_CONSTEXPR unsigned int my_digits10 = 101U;

  typedef boost::multiprecision::number<boost::multiprecision::cpp_dec_float<my_digits10>,
                                        boost::multiprecision::et_off>
  float_type;
};

BOOST_STATIC_ASSERT_MSG(local::my_digits10 > 20U,
                        "Error: This example is intended to have more than 20 decimal digits");

int main()
{
  // Use Gauss-Laguerre quadrature integration to compute airy_ai(x / 7)
  // with 7 <= x <= 120 and where x is incremented in steps of 1.

  // During development of this example, we have empirically found
  // the numbers of Gauss-Laguerre coefficients needed for convergence
  // when using various counts of base-10 digits.

  // Let's calibrate, for instance, the number of coefficients needed
  // at the point x = 1.

  // Empirical data lead to:
  // Fit[{{21.0, 3.5}, {51.0, 11.1}, {101.0, 22.5}, {201.0, 46.8}}, {1, d, d^2}, d]
  // FullSimplify[%]
  // -1.28301 + (0.235487 + 0.0000178915 d) d

  // We need significantly more coefficients at smaller real values than are needed
  // at larger real values because the slope derivative of airy_ai(x) gets more
  // steep as x approaches zero.

  // This Gauss-Laguerre quadrature is designed for airy_ai(x) with real-valued x >= 1.

  BOOST_CONSTEXPR_OR_CONST boost::float_least32_t d = static_cast<boost::float_least32_t>(std::numeric_limits<local::float_type>::digits10);

  BOOST_CONSTEXPR_OR_CONST boost::float_least32_t laguerre_order_factor = -1.28301F + ((0.235487F + (0.0000178915F * d)) * d);

  BOOST_CONSTEXPR_OR_CONST int laguerre_order = static_cast<int>(laguerre_order_factor * d);

  std::cout << "std::numeric_limits<local::float_type>::digits10: " << std::numeric_limits<local::float_type>::digits10 << std::endl;
  std::cout << "laguerre_order: " << laguerre_order << std::endl;

  typedef gauss::laguerre::detail::abscissas_and_weights<local::float_type> abscissas_and_weights_type;

  const abscissas_and_weights_type the_abscissas_and_weights(laguerre_order, local::float_type(-1) / 6);

  bool result_is_ok = true;

  for(std::uint32_t u = 7U; u < 121U; ++u)
  {
    const local::float_type x = local::float_type(u) / 7;

    typedef gauss::laguerre::detail::airy_ai_object<local::float_type> airy_ai_object_type;

    const airy_ai_object_type the_airy_ai_object(x);

    const local::float_type airy_ai_value =
      std::inner_product(the_abscissas_and_weights.abscissa_n().cbegin(),
                         the_abscissas_and_weights.abscissa_n().cend(),
                         the_abscissas_and_weights.weight_n().cbegin(),
                         local::float_type(0U),
                         std::plus<local::float_type>(),
                         [&the_airy_ai_object](const local::float_type& this_abscissa,
                                               const local::float_type& this_weight) -> local::float_type
                         {
                           return the_airy_ai_object(this_abscissa) * this_weight;
                         });

    static const local::float_type one_third = 1.0F / local::float_type(3U);

    static const local::float_type one_over_pi_times_one_over_sqrt_three =
      sqrt(one_third) * boost::math::constants::one_div_pi<local::float_type>();

    const local::float_type sqrt_x = sqrt(x);

    const local::float_type airy_ai_control =
       (sqrt_x * one_over_pi_times_one_over_sqrt_three)
      * boost::math::cyl_bessel_k(one_third, ((2.0F * x) * sqrt_x) * one_third);

    std::cout << std::setprecision(std::numeric_limits<local::float_type>::digits10)
              << "airy_ai_value  : "
              << airy_ai_value
              << std::endl;

    std::cout << std::setprecision(std::numeric_limits<local::float_type>::digits10)
              << "airy_ai_control: "
              << airy_ai_control
              << std::endl;

    const local::float_type delta = fabs(1.0F - (airy_ai_control / airy_ai_value));

    static const local::float_type tol("1E-" + boost::lexical_cast<std::string>(std::numeric_limits<local::float_type>::digits10 - 7U));

    result_is_ok &= (delta < tol);
  }

  std::cout << std::endl
            << "Total... result_is_ok: "
            << std::boolalpha
            << result_is_ok
            << std::endl;
} // int main()

/*


Partial output:

//[gauss_laguerre_quadrature_output_1

std::numeric_limits<local::float_type>::digits10: 101
laguerre_order: 2291

Finding the approximate roots...
root_estimates.size(): 1, 0.0%
root_estimates.size(): 8, 0.3%
root_estimates.size(): 16, 0.7%
...
root_estimates.size(): 2288, 99.9%
root_estimates.size(): 2291, 100.0%


Calculating abscissas and weights. Processed 1, 0.0%
Calculating abscissas and weights. Processed 9, 0.4%
...
Calculating abscissas and weights. Processed 2289, 99.9%
Calculating abscissas and weights. Processed 2291, 100.0%
//] [/gauss_laguerre_quadrature_output_1]

//[gauss_laguerre_quadrature_output_2

airy_ai_value  : 0.13529241631288141552414742351546630617494414298833070600910205475763353480226572366348710990874867334
airy_ai_control: 0.13529241631288141552414742351546630617494414298833070600910205475763353480226572366348710990874868323
airy_ai_value  : 0.11392302126009621102904231059693500086750049240884734708541630001378825889924647699516200868335286103
airy_ai_control: 0.1139230212600962110290423105969350008675004924088473470854163000137882588992464769951620086833528582
...
airy_ai_value  : 3.8990420982303275013276114626640705170145070824317976771461533035231088620152288641360519429331427451e-22
airy_ai_control: 3.8990420982303275013276114626640705170145070824317976771461533035231088620152288641360519429331426481e-22

Total... result_is_ok: true

//] [/gauss_laguerre_quadrature_output_2]


*/