diff --git a/doc/internals/barycentric_rational_interpolation.qbk b/doc/internals/barycentric_rational_interpolation.qbk new file mode 100644 index 0000000000..892d6983b4 --- /dev/null +++ b/doc/internals/barycentric_rational_interpolation.qbk @@ -0,0 +1,50 @@ +[section:barycentric Barycentric Rational Interpolation] + +[h4 Synopsis] + +`` +#include +`` + + +[h4 Description] + +Barycentric rational interpolation is a high-precision interpolation class for non-uniformly spaced samples. +It requires [bigo](N) time for construction, and [bigo](N) time for each evaluation. +Linear time evaluation is not optimal; for instance the cubic b spline can be evaluated in constant time. +However, using the cubic b spline requires uniformly spaced samples, which are not always available. + +Use of the class requires a vector of independent variables x[0], x[1], .... x[n-1] where x[i+1] > x[i], +and a vector of dependent variables y[0], y[1], ... , y[n-1]. +The call is trivial: + + __boost::math::tools::barycentric_rational interpolant(x.data(), y.data(), y.size()); + +This implicitly calls the constructor with approximation order 3, and hence the accuracy is [bigo](h^4). +In general, if you require an approximation order d, then the error is [bigo](h^d+1). +A call to the constructor with an explicit approximation order could be + + __boost::math::tools::barycentric_rational(x.data(), y.data(), y.size(), 5); + +To evaluate the interpolant, + +Although this algorithm is robust, it can surprise you. The main way this occurs is if there is some i such that x[i] - x[i-1] is many orders of magnitude smaller than x[i+1] - x[i]. +This can cause a large "swing" where the interpolant makes a fast dive to hit y[i], but then flies very high since it is unconstrained until y[i+1], which is far away. + +The precise quantification of this idea is encoded in the "local mesh ratio", discussed in Floater's paper on the barycentric rational interpolation. + +See + +* [@http://www.mn.uio.no/math/english/people/aca/michaelf/papers/rational.pdf Barycentric rational interpolation with no poles and a high rate of interpolation ] +* Floater, M.S. & Hormann, K. Numer. Math. (2007) 107: 315. doi:10.1007/s00211-007-0093-y + + +[endsect] [/section:barycentric Barycentric Rational Interpolation] + +[/ + Copyright 2017 Nick Thompson + + 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). +] diff --git a/doc/math.qbk b/doc/math.qbk index 2f4ceab3b0..50b6439cf6 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -658,6 +658,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22. [include background/error.qbk] [/relative error NOT handling] [include background/lanczos.qbk] [include background/remez.qbk] +[include quadrature/double_exponential.qbk] [include background/references.qbk] [section:logs_and_tables Error logs and tables] diff --git a/doc/quadrature/double_exponential.qbk b/doc/quadrature/double_exponential.qbk new file mode 100644 index 0000000000..170cc7336c --- /dev/null +++ b/doc/quadrature/double_exponential.qbk @@ -0,0 +1,206 @@ +[/ +Copyright (c) 2017 Nick Thompson +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) +] + +[section:double_exponential Double-exponential quadrature] + +[heading Synopsis] + +`` + #include + #include + #include + + namespace boost{ namespace math{ + + template + class tanh_sinh + { + public: + tanh_sinh(Real tol = sqrt(std::numeric_limits::epsilon()), size_t max_refinements = 15); + + template + Real integrate(const F f, Real a, Real b, Real* error = nullptr, Real* L1 = nullptr) const; + } + + template + class exp_sinh + { + public: + exp_sinh(Real tol = sqrt(std::numeric_limits::epsilon()), size_t max_refinements = 9); + + template + Real integrate(const F f, Real a = 0, Real b = std::numeric_limits::infinity(), Real* error = nullptr, Real* L1 = nullptr) const; + } + + template + class sinh_sinh + { + public: + sinh_sinh(Real tol = sqrt(std::numeric_limits::epsilon()), size_t max_refinements = 9); + + template + Real integrate(const F f, Real* error = nullptr, Real* L1 = nullptr) const; + } + +}} +`` + +[heading Overview] + +The tanh-sinh quadrature routine provided by boost is a rapidly convergent numerical integration scheme for holomorphic integrands. +By this we mean that the integrand is the restriction to the real line of a complex-differentiable function which is bounded on the interior of the unit disk /|z| < 1/. +More precisely, the integrands must lie within a Hardy space. +If your integrand obeys these conditions, it can be shown that tanh-sinh integration is optimal, +meaning that it requires the fewest function evaluations for a given accuracy of any quadrature algorithm for a random element from the Hardy space. +A basic example of how to use the tanh-sinh quadrature is shown below: + + tanh_sinh integrator; + auto f = [](double x) { return 5*x + 7; }; + double Q = integrator.integrate(f, (double) 0, (double) 1); + +The basic idea of tanh-sinh quadrature is that a variable transformation can cause the endpoint derivatives to decay rapidly. +When the derivatives at the endpoints decay much faster than the Bernoulli numbers grow, +the Euler-Maclaurin summation formula tells us that simple trapezoidal quadrature converges exponentially fast. + + +One very nice property of tanh-sinh quadrature is that it can handle singularities at the endpoints of the integration domain. +For instance, the following integrand, singular at both endpoints, can be efficiently evaluated to 100 binary digits: + + auto f = [](Real x) { return log(x)*log(1-x); }; + Real Q = integrator.integrate(f, (Real) 0, (Real) 1); + +Although boost's implementation of tanh-sinh quadrature can achieve very high precision, it is not truly arbitrary precision. +The reason is that vast speed improvements can be obtained by caching the weights and abscissas of the variable transformation + + tanh(half_pi()*sinh(t)); + +These numbers have been precomputed up to 100 binary digits, so this is roughly the most precision that can be expected from the quadrature. +We think that this is a reasonable compromise between efficiency, accuracy, and compile time that will satisfy the vast majority of users. + + +Now onto the caveats: As stated before, the integrands must lie in a Hardy space to ensure rapid convergence. +Attempting to integrate a function which is not bounded on the unit disk by tanh-sinh can lead to very slow convergence. +For example, take the Runge function: + + auto f1 = [](double t) { return 1/(1+25*t*t); }; + Q = integrator.integrate(f1, (double) -1, (double) 1); + +This function has poles at \u00B1 \u2148/5, and as such it is not bounded on the unit disk. +However, the related function + + auto f2 = [](double t) { return 1/(1+0.04*t*t); }; + Q = integrator.integrate(f2, (double) -1, (double) 1); + +has poles outside the unit disk (at \u00B1 5\u2148), and is therefore in the Hardy space. +The integration is performed 22x faster for the second function! +If you do not understand the structure of your integrand in the complex plane, do performance testing before deployment. + +Like the trapezoidal quadrature, the tanh-sinh quadrature produces an estimate of the L[sub 1] norm of the integral along with the requested integral. +This is to establish a scale against which to measure the tolerance, and provides an estimate of the condition number of the summation. +This can be queries as follows: + + + tanh_sinh integrator; + auto f = [](double x) { return 5*x + 7; }; + double error; + double L1; + double Q = integrator.integrate(f, (double) 0, (double) 1, &error, &L1); + double condition_number = L1/std::abs(Q); + +If the condition number is large, the computed integral is worthless. +Occasionally, the condition number will be so large that a `boost::math::evaluation_error` will be thrown. +However, the class is very conservative; +throwing only when the condition number is so large that the integral is essentially guaranteed to be worthless (when it exceeds the reciprocal of the unit roundoff). +Often, the computation is very poor well before the condition number exceeds this threshold, +so it is the responsibility of the user to check the L[sub 1] norm computation and evaluate whether it is acceptable. + +Although the tanh-sinh quadrature can compute integral over infinite domains by variable transformations, these transformations can create a very poorly behaved integrand. +For this reason, double-exponential variable transformations have been provided that allow stable computation over infinite domains; these being the exp-sinh and sinh-sinh quadrature. + +The sinh-sinh quadrature allows computation over the entire real line, and is called as follows: + + sinh_sinh integrator; + auto f = [](double x) { return exp(-x*x); }; + double error; + double L1; + double Q = integrator.integrate(f, &error, &L1); + +Note that the limits of integration are understood to be (-\u221E, \u221E). + +For half-infinite intervals, the exp-sinh quadrature is provided: + + exp_sinh integrator; + auto f = [](double x) { return exp(-3*x); }; + double error; + double L1; + double Q = integrator.integrate(f, (Real) 0, std::numeric_limits::infinity(), &error, &L1); + +Either the left limit or the right limit must be infinite; not both. +Endpoint singularities are supported by sinh-sinh and exp-sinh. +These quadrature algorithms are truly arbitrary precision, but this requires the abscissas and weights to be computed at runtime. +This means that the call to the constructor is expensive, so make sure to reuse the object for multiple integrations (if, say you are assembling a stiffness matrix). +Subsequent calls the `integrator.integrate` are threadsafe, so assembling a stiffness matrix can be done in parallel. + +[heading Tolerance and Max-interval Halvings] + +The constructor for all three double-exponential quadratures supports two additional arguments: A tolerance and max-interval halvings. +The tolerance is met when two subsequent estimates of the integral have absolute error less than `tol*L1`. +It is highly recommended that the tolerance be left at the default value of [radic][epsilon]. +Since double exponential quadrature converges exponentially fast for functions in Hardy spaces, then once the routine has *proved* that the error is ~[radic][epsilon], +then the error in fact is ~[epsilon]. +If you request that the error be ~[epsilon], this tolerance might never be achieved (as the summation is not stabilized ala Kahan), and the routine will simply flounder, +dividing the interval in half in order to increase the precision of the integrand, only to be thwarted by floating point roundoff. + +The max interval halvings is the maximum number number of times the interval can be cut in half before giving up. +If the accuracy is not met at that time, the routine returns its best estimate, along with the `error` and `L1`, +which allows the user to decide if another quadrature routine should be employed. + +An example of this is + + double tol = std::sqrt(std::numeric_limits::epsilon()); + size_t max_halvings = 12; + tanh_sinh integrator(tol, max_halvings); + auto f = [](double x) { return 5*x + 7; }; + double error, L1; + double Q = integrator.integrate(f, (double) 0, (double) 1, &error, &L1); + if (error*L1 > 0.01) + { + Q = some_other_quadrature_method(f, (double) 0, (double) 1); + } + +[heading Caveats] + +A few things to keep in mind while using the tanh-sinh, exp-sinh, and sinh-sinh quadratures: + +These routines are *very* aggressive about approaching the endpoint singularities. +This allows lots of significant digits to be extracted, but also has another problem: Roundoff error can cause the function to be evaluated at the endpoints. +A few ways to avoid this: Narrow up the bounds of integration to say, [a + [epsilon], b - [epsilon]], make sure (a+b)/2 and (b-a)/2 are representable, and finally, if you think the compromise between accuracy an usability has gone too far in the direction of accuracy, file a ticket. + +Both exp-sinh and sinh-sinh quadratures evaluate the functions they are passed at *very* large argument. +You might understand that x[super 12]exp(-x) is should be zero when x[super 12] overflows, but IEEE floating point arithmetic does not. +Hence `std::pow(x, 12)*std::exp(-x)` is an indeterminate form whenever `std::pow(x, 12)` overflows. +So make sure your functions have the correct limiting behavior; for example + + auto f = [](double x) { + double t = exp(-x); + if(t == 0) + { + return 0; + } + return t*pow(x, 12); + }; + +has the correct behavior for large /x/, but `auto f = [](double x) { return exp(-x)*pow(x, 12); };` does not. + +Oscillatory integrals, such as the sinc integral, are poorly approximated by double-exponential quadrature. +Fortunately the error estimates and L1 norm are massive for these integrals, but nonetheless, oscillatory integrals require different techniques. + +References: + +Tanaka, Ken’ichiro, et al. ['Function classes for double exponential integration formulas.] Numerische Mathematik 111.4 (2009): 631-655. + +[endsect] diff --git a/example/barycentric_interpolation_example.cpp b/example/barycentric_interpolation_example.cpp new file mode 100644 index 0000000000..06ac3267d9 --- /dev/null +++ b/example/barycentric_interpolation_example.cpp @@ -0,0 +1,87 @@ + +// Copyright Nick Thompson, 2017 + +// 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). + +#include +#include +#include + +//[barycentric_rational_example + +/*`This example shows how to use barycentric rational interpolation, using Walter Kohn's classic paper + "Solution of the Schrodinger Equation in Periodic Lattices with an Application to Metallic Lithium" + In this paper, Kohn needs to repeatedly solve an ODE (the radial Schrodinger equation) given a potential + which is only known at non-equally samples data. + + If he'd only had the barycentric rational interpolant of boost::math! + References: + Kohn, W., and N. Rostoker. "Solution of the Schrödinger equation in periodic lattices with an application to metallic lithium." Physical Review 94.5 (1954): 1111. +*/ + +#include +//] [/barycentric_rational_example] + +int main() +{ + // The lithium potential is given in Kohn's paper, Table I: + std::vector r(45); + std::vector mrV(45); + + r[0] = 0.02; mrV[0] = 5.727; + r[1] = 0.04, mrV[1] = 5.544; + r[2] = 0.06, mrV[2] = 5.450; + r[3] = 0.08, mrV[3] = 5.351; + r[4] = 0.10, mrV[4] = 5.253; + r[5] = 0.12, mrV[5] = 5.157; + r[6] = 0.14, mrV[6] = 5.058; + r[7] = 0.16, mrV[7] = 4.960; + r[8] = 0.18, mrV[8] = 4.862; + r[9] = 0.20, mrV[9] = 4.762; + r[10] = 0.24, mrV[10] = 4.563; + r[11] = 0.28, mrV[11] = 4.360; + r[12] = 0.32, mrV[12] = 4.1584; + r[13] = 0.36, mrV[13] = 3.9463; + r[14] = 0.40, mrV[14] = 3.7360; + r[15] = 0.44, mrV[15] = 3.5429; + r[16] = 0.48, mrV[16] = 3.3797; + r[17] = 0.52, mrV[17] = 3.2417; + r[18] = 0.56, mrV[18] = 3.1209; + r[19] = 0.60, mrV[19] = 3.0138; + r[20] = 0.68, mrV[20] = 2.8342; + r[21] = 0.76, mrV[21] = 2.6881; + r[22] = 0.84, mrV[22] = 2.5662; + r[23] = 0.92, mrV[23] = 2.4242; + r[24] = 1.00, mrV[24] = 2.3766; + r[25] = 1.08, mrV[25] = 2.3058; + r[26] = 1.16, mrV[26] = 2.2458; + r[27] = 1.24, mrV[27] = 2.2035; + r[28] = 1.32, mrV[28] = 2.1661; + r[29] = 1.40, mrV[29] = 2.1350; + r[30] = 1.48, mrV[30] = 2.1090; + r[31] = 1.64, mrV[31] = 2.0697; + r[32] = 1.80, mrV[32] = 2.0466; + r[33] = 1.96, mrV[33] = 2.0325; + r[34] = 2.12, mrV[34] = 2.0288; + r[35] = 2.28, mrV[35] = 2.0292; + r[36] = 2.44, mrV[36] = 2.0228; + r[37] = 2.60, mrV[37] = 2.0124; + r[38] = 2.76, mrV[38] = 2.0065; + r[39] = 2.92, mrV[39] = 2.0031; + r[40] = 3.08, mrV[40] = 2.0015; + r[41] = 3.24, mrV[41] = 2.0008; + r[42] = 3.40, mrV[42] = 2.0004; + r[43] = 3.56, mrV[43] = 2.0002; + r[44] = 3.72, mrV[44] = 2.0001; + + // Now we want to interpolate this potential at any r: + boost::math::barycentric_rational b(r.data(), mrV.data(), r.size()); + + for (size_t i = 1; i < 8; ++i) + { + double r = i*0.5; + std::cout << "(r, V) = (" << r << ", " << -b(r)/r << ")\n"; + } +} diff --git a/include/boost/math/interpolators/barycentric_rational.hpp b/include/boost/math/interpolators/barycentric_rational.hpp new file mode 100644 index 0000000000..83116cd960 --- /dev/null +++ b/include/boost/math/interpolators/barycentric_rational.hpp @@ -0,0 +1,60 @@ +/* + * Copyright Nick Thompson, 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) + * + * Given N samples (t_i, y_i) which are irregularly spaced, this routine constructs an + * interpolant s which is constructed in O(N) time, occupies O(N) space, and can be evaluated in O(N) time. + * The interpolation is stable, unless one point is incredibly close to another, and the next point is incredibly far. + * The measure of this stability is the "local mesh ratio", which can be queried from the routine. + * Pictorially, the following t_i spacing is bad (has a high local mesh ratio) + * || | | | | + * and this t_i spacing is good (has a low local mesh ratio) + * | | | | | | | | | | + * + * + * If f is C^{d+2}, then the interpolant is O(h^(d+1)) accurate, where d is the interpolation order. + * A disadvantage of this interpolant is that it does not reproduce rational functions; for example, 1/(1+x^2) is not interpolated exactly. + * + * References: + * Floater, Michael S., and Kai Hormann. "Barycentric rational interpolation with no poles and high rates of approximation." Numerische Mathematik 107.2 (2007): 315-331. + * Press, William H., et al. "Numerical recipes third edition: the art of scientific computing." Cambridge University Press 32 (2007): 10013-2473. + */ + +#ifndef BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_HPP +#define BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_HPP + +#include +#include + +namespace boost{ namespace math{ + +template +class barycentric_rational +{ +public: + barycentric_rational(const Real* const x, const Real* const y, size_t n, size_t approximation_order = 3); + + Real operator()(Real x) const; + +private: + std::shared_ptr> m_imp; +}; + +template +barycentric_rational::barycentric_rational(const Real* const x, const Real* const y, size_t n, size_t approximation_order): + m_imp(std::make_shared>(x, y, n, approximation_order)) +{ + return; +} + +template +Real barycentric_rational::operator()(Real x) const +{ + return m_imp->operator()(x); +} + + +}} +#endif diff --git a/include/boost/math/interpolators/detail/barycentric_rational_detail.hpp b/include/boost/math/interpolators/detail/barycentric_rational_detail.hpp new file mode 100644 index 0000000000..b0a202f0f0 --- /dev/null +++ b/include/boost/math/interpolators/detail/barycentric_rational_detail.hpp @@ -0,0 +1,154 @@ +/* + * Copyright Nick Thompson, 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) + */ + +#ifndef BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_DETAIL_HPP +#define BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_DETAIL_HPP + +#include +#include +#include +#include + +namespace boost{ namespace math{ namespace detail{ + +template +class barycentric_rational_imp +{ +public: + barycentric_rational_imp(const Real* const x, const Real* const y, size_t n, size_t approximation_order = 3); + + Real operator()(Real x) const; + + // The barycentric weights are not really that interesting; except to the unit tests! + Real weight(size_t i) const { return m_w[i]; } + +private: + // Technically, we do not need to copy over y to m_y, or x to m_x. + // We could simply store a pointer. However, in doing so, + // we'd need to make sure the user managed the lifetime of m_y, + // and didn't alter its data. Since we are unlikely to run out of + // memory for a linearly scaling algorithm, it seems best just to make a copy. + std::vector m_y; + std::vector m_x; + std::vector m_w; +}; + +template +barycentric_rational_imp::barycentric_rational_imp(const Real* const x, const Real* const y, size_t n, size_t approximation_order) +{ + using std::abs; + + if (approximation_order >= n) + { + throw std::domain_error("Approximation order must be < data length."); + } + + if (x == nullptr) + { + throw std::domain_error("Independent variable passed to barycentric_rational is a nullptr"); + } + + if (y == nullptr) + { + throw std::domain_error("Dependent variable passed to barycentric_rational is a nullptr"); + } + + // Big sad memcpy to make sure the object is easy to use. + m_x.resize(n); + m_y.resize(n); + for(size_t i = 0; i < n; ++i) + { + // But if we're going to do a memcpy, we can do some error checking which is inexpensive relative to the copy: + if(boost::math::isnan(x[i])) + { + boost::format fmtr = boost::format("x[%1%] is a NAN") % i; + throw std::domain_error(fmtr.str()); + } + + if(boost::math::isnan(y[i])) + { + boost::format fmtr = boost::format("y[%1%] is a NAN") % i; + throw std::domain_error(fmtr.str()); + } + + m_x[i] = x[i]; + m_y[i] = y[i]; + } + + m_w.resize(n, 0); + for(int64_t k = 0; k < n; ++k) + { + int64_t i_min = std::max(k - (int64_t) approximation_order, (int64_t) 0); + int64_t i_max = k; + if (k >= n - approximation_order) + { + i_max = n - approximation_order - 1; + } + + for(int64_t i = i_min; i <= i_max; ++i) + { + Real inv_product = 1; + int64_t j_max = std::min(i + approximation_order, n - 1); + for(int64_t j = i; j <= j_max; ++j) + { + if (j == k) + { + continue; + } + + Real diff = m_x[k] - m_x[j]; + if (abs(diff) < std::numeric_limits::epsilon()) + { + boost::format fmtr = boost::format("Spacing between x[%1%] and x[%2%] is %3%, which is smaller than the epsilon of %4%") % k % i % diff % boost::typeindex::type_id().pretty_name(); + std::cout << typeid(fmtr).name() << std::endl; + throw std::logic_error(fmtr.str()); + } + inv_product *= diff; + } + if (i % 2 == 0) + { + m_w[k] += 1/inv_product; + } + else + { + m_w[k] -= 1/inv_product; + } + } + } +} + +template +Real barycentric_rational_imp::operator()(Real x) const +{ + Real numerator = 0; + Real denominator = 0; + for(size_t i = 0; i < m_x.size(); ++i) + { + // Presumably we should see if the accuracy is improved by using ULP distance of say, 5 here, instead of testing for floating point equality. + // However, it has been shown that if x approx x_i, but x != x_i, then inaccuracy in the numerator cancels the inaccuracy in the denominator, + // and the result is fairly accurate. See: http://epubs.siam.org/doi/pdf/10.1137/S0036144502417715 + if (x == m_x[i]) + { + return m_y[i]; + } + Real t = m_w[i]/(x - m_x[i]); + numerator += t*m_y[i]; + denominator += t; + } + return numerator/denominator; +} + +/* + * A formula for computing the derivative of the barycentric representation is given in + * "Some New Aspects of Rational Interpolation", by Claus Schneider and Wilhelm Werner, + * Mathematics of Computation, v47, number 175, 1986. + * http://www.ams.org/journals/mcom/1986-47-175/S0025-5718-1986-0842136-8/S0025-5718-1986-0842136-8.pdf + * However, this requires a lot of machinery which is not built into the library at present. + * So we wait until there is a requirement to interpolate the derivative. + */ +}}} +#endif diff --git a/include/boost/math/quadrature/detail/exp_sinh_detail.hpp b/include/boost/math/quadrature/detail/exp_sinh_detail.hpp new file mode 100644 index 0000000000..820a43b7c3 --- /dev/null +++ b/include/boost/math/quadrature/detail/exp_sinh_detail.hpp @@ -0,0 +1,205 @@ +#ifndef BOOST_MATH_QUADRATURE_DETAIL_EXP_SINH_DETAIL_HPP +#define BOOST_MATH_QUADRATURE_DETAIL_EXP_SINH_DETAIL_HPP + +#include +#include +#include +#include +#include +#include + +namespace boost{ namespace math{ namespace detail{ + + +// Returns the exp-sinh quadrature of a function f over the open interval (0, infinity) + +template +class exp_sinh_detail +{ +public: + exp_sinh_detail(Real tol, size_t max_refinements); + + template + Real integrate(const F f, Real* error, Real* L1) const; + +private: + Real m_tol; + + std::vector> m_abscissas; + std::vector> m_weights; +}; + +template +exp_sinh_detail::exp_sinh_detail(Real tol, size_t max_refinements) +{ + using std::exp; + using std::log; + using std::sqrt; + using std::asinh; + using std::cosh; + using std::sinh; + using std::ceil; + using boost::math::constants::two_div_pi; + using boost::math::constants::half_pi; + + m_tol = tol; + std::cout << std::setprecision(std::numeric_limits::digits10); + // t_min is chosen such that x = exp(pi/2 sinh(-t_max)) = eps/10000. + // This is a compromise; using some denormals so as to pick up information about singularities at zero, + // while not risking hitting the singularity through roundoff error. + const Real t_min = asinh(two_div_pi()*log(std::numeric_limits::epsilon()/10000)); + + // t_max is chosen to make g'(t_max) ~ sqrt(max) (g' grows faster than g). + // This will allow some flexibility on the users part; they can at least square a number function without overflow. + // But there is no unique choice; the further out we can evaluate the function, the better we can do on slowly decaying integrands. + const Real t_max = log(2*two_div_pi()*log(2*two_div_pi()*sqrt(std::numeric_limits::max()))); + m_abscissas.resize(max_refinements); + m_weights.resize(max_refinements); + + for (size_t i = 0; i < max_refinements; ++i) + { + Real h = (Real) 1/ (Real) (1<()*sinh(arg)); + m_abscissas[i].emplace_back(x); + Real w = cosh(arg)*half_pi()*x; + m_weights[i].emplace_back(w); + + // TODO: Repeatedly adding h is the incorrect method; it accumalates roundoff error. + if (i != 0) + { + arg = arg + 2*h; + } + else + { + arg = arg + h; + } + } + } +} + +template +template +Real exp_sinh_detail::integrate(const F f, Real* error, Real* L1) const +{ + using std::abs; + using std::floor; + using std::tanh; + using std::sinh; + using std::sqrt; + using boost::math::constants::half; + using boost::math::constants::half_pi; + + Real y_max = f(std::numeric_limits::max()); + if(abs(y_max) > std::numeric_limits::epsilon() || !isfinite(y_max)) + { + boost::basic_format err_msg = boost::format("\nThe function you are trying to integrate does not go to zero at infinity, and instead evaluates to %1% at x=%2%.\n") % y_max % std::numeric_limits::max(); + throw std::domain_error(err_msg.str()); + } + + // std::cout << std::setprecision(std::numeric_limits::digits10); + + // Get the party started with two estimates of the integral: + Real I0 = 0; + Real L1_I0 = 0; + for(size_t i = 0; i < m_abscissas[0].size(); ++i) + { + Real y = f(m_abscissas[0][i]); + I0 += y*m_weights[0][i]; + L1_I0 += abs(y)*m_weights[0][i]; + } + + // std::cout << "First estimate : " << I0 << std::endl; + Real I1 = I0; + Real L1_I1 = L1_I0; + for (size_t i = 0; i < m_abscissas[1].size(); ++i) + { + Real y = f(m_abscissas[1][i]); + I1 += y*m_weights[1][i]; + L1_I1 += abs(y)*m_weights[1][i]; + } + + I1 *= half(); + L1_I1 *= half(); + Real err = abs(I0 - I1); + // std::cout << "Second estimate: " << I1 << " Error estimate at level " << 1 << " = " << err << std::endl; + + for(size_t i = 2; i < m_abscissas.size(); ++i) + { + I0 = I1; + L1_I0 = L1_I1; + + I1 = half()*I0; + L1_I1 = half()*L1_I0; + Real h = (Real) 1/ (Real) (1 << i); + Real sum = 0; + Real absum = 0; + + Real abterm1 = 1; + Real eps = std::numeric_limits::epsilon()*L1_I1; + for(size_t j = 0; j < m_weights[i].size(); ++j) + { + Real x = m_abscissas[i][j]; + Real y = f(x); + sum += y*m_weights[i][j]; + Real abterm0 = abs(y)*m_weights[i][j]; + absum += abterm0; + + // We require two consecutive terms to be < eps in case we hit a zero of f. + // Numerical experiments indicate that we should start this check at ~30 for floats, + // ~60 for doubles, and ~100 for long doubles. + // However, starting the check at x = 10 rather than x = 100 will only save two function evaluations. + if (x > (Real) 100 && abterm0 < eps && abterm1 < eps) + { + break; + } + abterm1 = abterm0; + } + + I1 += sum*h; + L1_I1 += absum*h; + err = abs(I0 - I1); + // std::cout << "Estimate: " << I1 << " Error estimate at level " << i << " = " << err << std::endl; + if (!isfinite(I1)) + { + throw std::domain_error("The exp_sinh quadrature evaluated your function at a singular point. Please ensure your function evaluates to a finite number of its entire domain.\n"); + } + if (err <= m_tol*L1_I1 || err <= std::numeric_limits::epsilon() /* catch the 0 case. */) + { + break; + } + } + + if (error) + { + *error = err; + } + + if(L1) + { + *L1 = L1_I1; + } + + if(L1_I1 != (Real) 0 && L1_I1/abs(I1) > 1/std::numeric_limits::epsilon()) + { + Real cond = abs(I1)/L1_I1; + Real inv_prec = 1.0/std::numeric_limits::epsilon(); + boost::basic_format err_msg = boost::format("\nThe condition number of the quadrature sum is %1%, which exceeds the inverse of the precision of the type %2%.\nNo correct digits can be expected for the integral.\n") % cond % inv_prec; + throw boost::math::evaluation_error(err_msg.str()); + } + + return I1; +} + +}}} +#endif diff --git a/include/boost/math/quadrature/detail/sinh_sinh_detail.hpp b/include/boost/math/quadrature/detail/sinh_sinh_detail.hpp new file mode 100644 index 0000000000..59d4c2cd9f --- /dev/null +++ b/include/boost/math/quadrature/detail/sinh_sinh_detail.hpp @@ -0,0 +1,207 @@ +#ifndef BOOST_MATH_QUADRATURE_DETAIL_SINH_SINH_DETAIL_HPP +#define BOOST_MATH_QUADRATURE_DETAIL_SINH_SINH_DETAIL_HPP + +#include +#include +#include +#include +#include +#include +#include + +namespace boost{ namespace math{ namespace detail{ + + +// Returns the sinh-sinh quadrature of a function f over the entire real line + +template +class sinh_sinh_detail +{ +public: + sinh_sinh_detail(Real tol, size_t max_refinements); + + template + Real integrate(const F f, Real* error, Real* L1) const; + +private: + Real m_tol; + + std::vector> m_abscissas; + std::vector> m_weights; +}; + +template +sinh_sinh_detail::sinh_sinh_detail(Real tol, size_t max_refinements) +{ + using std::log; + using std::sqrt; + using std::cosh; + using std::sinh; + using boost::math::constants::two_div_pi; + using boost::math::constants::half_pi; + + m_tol = tol; + + // t_max is chosen to make g'(t_max) ~ sqrt(max) (g' grows faster than g). + // This will allow some flexibility on the users part; they can at least square a number function without overflow. + // But there is no unique choice; the further out we can evaluate the function, the better we can do on slowly decaying integrands. + const Real t_max = log(2*two_div_pi()*log(2*two_div_pi()*sqrt(std::numeric_limits::max()))); + m_abscissas.resize(max_refinements); + m_weights.resize(max_refinements); + + for (size_t i = 0; i < max_refinements; ++i) + { + Real h = (Real) 1/ (Real) (1<()*sinh(arg); + Real x = sinh(tmp); + m_abscissas[i].emplace_back(x); + Real w = cosh(arg)*half_pi()*cosh(tmp); + m_weights[i].emplace_back(w); + + // TODO: Repeated addition of h accumulates roundoff error. + if (i != 0) + { + arg = arg + 2*h; + } + else + { + arg = arg + h; + } + + } + + } +} + +template +template +Real sinh_sinh_detail::integrate(const F f, Real* error, Real* L1) const +{ + using std::abs; + using std::sqrt; + using std::isfinite; + using boost::math::constants::half; + using boost::math::constants::half_pi; + + Real y_max = f(std::numeric_limits::max()); + if(abs(y_max) > std::numeric_limits::epsilon() || !isfinite(y_max)) + { + boost::basic_format err_msg = boost::format("\nThe function you are trying to integrate does not go to zero at infinity, and instead evaluates to %1% at x=%2%.\n") % y_max % std::numeric_limits::max(); + throw std::domain_error(err_msg.str()); + } + + Real y_min = f(std::numeric_limits::lowest()); + if(abs(y_min) > std::numeric_limits::epsilon() || !isfinite(y_min)) + { + boost::basic_format err_msg = boost::format("\nThe function you are trying to integrate does not go to zero at infinity, and instead evaluates to %1% at x=%2%.\n") % y_min % std::numeric_limits::lowest(); + throw std::domain_error(err_msg.str()); + } + + // Get the party started with two estimates of the integral: + Real I0 = f(0)*half_pi(); + Real L1_I0 = abs(I0); + for(size_t i = 0; i < m_abscissas[0].size(); ++i) + { + Real x = m_abscissas[0][i]; + Real yp = f(x); + Real ym = f(-x); + I0 += (yp + ym)*m_weights[0][i]; + L1_I0 += (abs(yp)+abs(ym))*m_weights[0][i]; + } + + // Uncomment the estimates to work the convergence on the command line. + // std::cout << std::setprecision(std::numeric_limits::digits10); + // std::cout << "First estimate : " << I0 << std::endl; + Real I1 = I0; + Real L1_I1 = L1_I0; + for (size_t i = 0; i < m_abscissas[1].size(); ++i) + { + Real x= m_abscissas[1][i]; + Real yp = f(x); + Real ym = f(-x); + I1 += (yp + ym)*m_weights[1][i]; + L1_I1 += (abs(yp) + abs(ym))*m_weights[1][i]; + } + + I1 *= half(); + L1_I1 *= half(); + Real err = abs(I0 - I1); + // std::cout << "Second estimate: " << I1 << " Error estimate at level " << 1 << " = " << err << std::endl; + + for(size_t i = 2; i < m_abscissas.size(); ++i) + { + I0 = I1; + L1_I0 = L1_I1; + + I1 = half()*I0; + L1_I1 = half()*L1_I0; + Real h = (Real) 1/ (Real) (1 << i); + Real sum = 0; + Real absum = 0; + + Real abterm1 = 1; + Real eps = std::numeric_limits::epsilon()*L1_I1; + for(size_t j = 0; j < m_weights[i].size(); ++j) + { + Real x = m_abscissas[i][j]; + Real yp = f(x); + Real ym = f(-x); + sum += (yp + ym)*m_weights[i][j]; + Real abterm0 = (abs(yp) + abs(ym))*m_weights[i][j]; + absum += abterm0; + + // We require two consecutive terms to be < eps in case we hit a zero of f. + if (x > (Real) 100 && abterm0 < eps && abterm1 < eps) + { + break; + } + abterm1 = abterm0; + } + + I1 += sum*h; + L1_I1 += absum*h; + err = abs(I0 - I1); + // std::cout << "Estimate: " << I1 << " Error estimate at level " << i << " = " << err << std::endl; + if (!isfinite(I1)) + { + std::string err_msg = "\nThe sinh_sinh quadrature evaluated your function at a singular point.\n"; + err_msg += "sinh_sinh quadrature cannot handle singularities in the domain.\n"; + err_msg += "If you are sure your function has no singularities, please submit a bug against boost.math\n"; + throw std::domain_error(err_msg); + } + if (err <= m_tol*L1_I1) + { + break; + } + } + + if (error) + { + *error = err; + } + + if (L1) + { + *L1 = L1_I1; + } + + Real cond = abs(I1)/L1_I1; + Real inv_prec = 1.0/std::numeric_limits::epsilon(); + if (L1_I1 != (Real) 0 && cond > inv_prec) + { + boost::basic_format err_msg = boost::format("\nThe condition number of the quadrature sum is %1%, which exceeds the inverse of the precision of the type %2%.\nNo correct digits can be expected for the integral.\n") % cond % inv_prec; + throw boost::math::evaluation_error(err_msg.str()); + } + + return I1; +} + +}}} +#endif diff --git a/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp b/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp new file mode 100644 index 0000000000..ec9749a9db --- /dev/null +++ b/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp @@ -0,0 +1,911 @@ +#ifndef BOOST_MATH_QUADRATURE_DETAIL_TANH_SINH_DETAIL_HPP +#define BOOST_MATH_QUADRATURE_DETAIL_TANH_SINH_DETAIL_HPP + +#include +#include +#include +#include + + +namespace boost{ namespace math{ namespace detail{ + + +// Returns the tanh-sinh quadrature of a function f over the open interval (-1, 1) + +template +class tanh_sinh_detail +{ +public: + tanh_sinh_detail(Real tol, size_t max_refinements); + + template + Real integrate(const F f, Real* error, Real* L1) const; + +private: + Real m_tol; + Real m_t_max; + size_t m_max_refinements; + + const std::vector> m_abscissas{ + { + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"), // g(0) + lexical_cast("0.951367964072746945727055362904639667492765811307212865380079106867050650113429723656597692697291568999499"), // g(1) + lexical_cast("0.999977477192461592864899636308688849285982470957489530009950811164291603926066803141755297920571692976244"), // g(2) + lexical_cast("0.999999999999957058389441217592223051901253805502353310119906858210171731776098247943305316472169355401371"), // g(3) + lexical_cast("0.999999999999999999999999999999999999883235110249013906725663510362671044720752808415603434458608954302982"), // g(4) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999988520"), // g(5) + }, + { + lexical_cast("0.674271492248435826080420090632052144267244477154740136576044169121421222924677035505101849981126759183585"), // g(0.5) + lexical_cast("0.997514856457224386832717419238820368149231215381809295391585321457677448277585202664213469426402672227688"), // g(1.5) + lexical_cast("0.999999988875664881984668015033322737014902900245095922058323073481945768209599289672119775131473502267885"), // g(2.5) + lexical_cast("0.999999999999999999999946215084086063112254432391666747077319911504923816981286361121293026490456057993796"), // g(3.5) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999920569786807778838966034206923747918174840316"), // g(4.5) + }, + { + lexical_cast("0.377209738164034173791478637625934394046524589091544494109155461140878077301078976168419776169200094544589"), // g(0.25) + lexical_cast("0.859569058689896635174660197262069516075083168252338276643044595162944527564048480673943822401281000252926"), // g(0.75) + lexical_cast("0.987040560507376891687385718614135001229151758131034693314584168837680972959341356928470215492975267090176"), // g(1.25) + lexical_cast("0.999688264028353209050521178688807537891398170976819589774088477780370538141423443919794848343175486617410"), // g(1.75) + lexical_cast("0.999999204737114712664465542942453822161328064663254514326490697871837585328183799811639137651868933587404"), // g(2.25) + lexical_cast("0.999999999952856448176777424944457862129587495731371252310623438621153306790574628302457884905941423030536"), // g(2.75) + lexical_cast("0.999999999999999994584777176192769154034548477151229790685620860861730773100739128006210236909975145697421"), // g(3.25) + lexical_cast("0.999999999999999999999999999979596996056475056705929184501160957954417209773314067152744883862197472178031"), // g(3.75) + lexical_cast("0.999999999999999999999999999999999999999999999996940309864635549970180556579324161834657856283984500867338"), // g(4.25) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999999999713780158074124883878031702"), // g(4.75) + }, + { + lexical_cast("0.194357003324935431614643585437365635031489393758436888182689674537462819812964064346338646610576297581586"), // g(0.125) + lexical_cast("0.539146705387967769049755260308830589945204337905779109884127027471561734390039416581961384016815411085580"), // g(0.375) + lexical_cast("0.780607438983200299254796988339627743834947936136904016926299134081823852355991340267593451506957297252642"), // g(0.625) + lexical_cast("0.914879263264574610907375523744910820577829749873996555598632174862675662265560700507766434443713921307207"), // g(0.875) + lexical_cast("0.973966868195677448563302789223612658807481226082861590717047117507988640718384015095866496549988100356307"), // g(1.125) + lexical_cast("0.994055506631402143292689401897413806059678470641395104619869077085166056854518863874214331797351767289690"), // g(1.375) + lexical_cast("0.999065196455785846424759456709674994467769868654311697576247713203495811920284255853592032063572034922686"), // g(1.625) + lexical_cast("0.999909384695143999838580661336134716111211722349817780096644956892300806327379163871993823738779424348069"), // g(1.875) + lexical_cast("0.999995316041220528430078670445852706137191337226395937272536393951221809216824715464045171619432016510042"), // g(2.125) + lexical_cast("0.999999892781612418381905659520380205133858442023594319170176418926992040688171632229230452149685824165825"), // g(2.375) + lexical_cast("0.999999999142705092178322670809525899422862401124370994655071602795868824046322936142150835978433636746740"), // g(2.625) + lexical_cast("0.999999999998232165307161279630254920657948830921426033212848744385677415250369788914260210901037491192034"), // g(2.875) + lexical_cast("0.999999999999999362512150495556035217718543739839598370606597786212238422008494757151263655913568329795137"), // g(3.125) + lexical_cast("0.999999999999999999975570627209147826822539587252630093662718420988465083426398009049228737368581275877248"), // g(3.375) + lexical_cast("0.999999999999999999999999947484535269804251437086861669491172732379804801794867463751892542463157039399705"), // g(3.625) + lexical_cast("0.999999999999999999999999999999997210532837710582255794317318253110654027520929260603115511315387451260374"), // g(3.875) + lexical_cast("0.999999999999999999999999999999999999999998721889101906195502025128665366339855321543654319959192027632657"), // g(4.125) + lexical_cast("0.999999999999999999999999999999999999999999999999999998691727663146819195361825683537382068319470132598104"), // g(4.375) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999472002187738974562351500885507843247"), // g(4.625) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999990938089594581899"), // g(4.875) + }, + { + lexical_cast("0.097923885287832333262426257841800739277991416860553014192169233925851692222761737403400646946996810714243"), // g(0.0625) + lexical_cast("0.287879932742715914564047412640584527472150498931093229176206925153530639212626250201384282745863057140995"), // g(0.1875) + lexical_cast("0.461253543939585704403089605478914762863209230623056694674993009583620156127192848514703129689772208211723"), // g(0.3125) + lexical_cast("0.610273657500638944882295196142322929272126071895168709332784368858457580338160458175008100036181213129260"), // g(0.4375) + lexical_cast("0.731018034792561511486243174126473035522958183942327321173501458391658882821561878668831641953955876427689"), // g(0.5625) + lexical_cast("0.823317005506402370062529821336585278778686354622394018907952364068212256641462115457246024984689378719324"), // g(0.6875) + lexical_cast("0.889891402784260198076869185915921905547677505107904531082251186493539771393653372732219535266764055746948"), // g(0.8125) + lexical_cast("0.935160857521984683228604353505559414643233174259410614443095874852547155290772317379778701463465547117041"), // g(0.9375) + lexical_cast("0.964112164223547291931468134197655554746442479999771045384745514462535509114997136284568448771488182152344"), // g(1.0625) + lexical_cast("0.981454826677335170026708563678276787119098621035665101299653042509293471745034185419176689483762517329856"), // g(1.1875) + lexical_cast("0.991126992441698802231143679454846914217586259006443485654561472284595775494251990563370501434560309592333"), // g(1.3125) + lexical_cast("0.996108665437508542535621637146298489121149585549952407779989838186439679719334321312223452465221335939924"), // g(1.4375) + lexical_cast("0.998454208767697737510844453998095652865528721174750512677554577296976878473984948660132937507588562719585"), // g(1.5625) + lexical_cast("0.999451434435274605841915072392823257052451949418516049867244221879506069754651427125490290081959377494427"), // g(1.6875) + lexical_cast("0.999828822072874941659766880829150059268532266322457928470788602237570528282362735805474461821547990270187"), // g(1.8125) + lexical_cast("0.999953871005627960747340800441686620980855509673364243556079565213417238874480960295627275099729241808849"), // g(1.9375) + lexical_cast("0.999989482014818503611471895899007720602902492343952737673740948074896439547231461927799345322653909942883"), // g(2.0625) + lexical_cast("0.999998017140595432079371981517432411433781680968439373654166553518514233346092449748225001874140137160282"), // g(2.1875) + lexical_cast("0.999999698894152611220345020698910138201776866792200349281689337716739578683277465315125136688228484355588"), // g(2.3125) + lexical_cast("0.999999964239080915342285626382935410646072872602826135584664520982694595347321344813284653649501363550641"), // g(2.4375) + lexical_cast("0.999999996787199098304239500034340852599841652502343477742523446571681728437217933803884887492857151527030"), // g(2.5625) + lexical_cast("0.999999999789732862235885118956133030249877860325502789814259697102439593977121468853074906349395599492798"), // g(2.6875) + lexical_cast("0.999999999990393933521491534996761230126495752035058750371207104450643553909962268167080693233470464022371"), // g(2.8125) + lexical_cast("0.999999999999708097335898248606616341262414631780830812553586730476909924106339053824128432576795745904030"), // g(2.9375) + lexical_cast("0.999999999999994413883361161772478258894903291272354533648485685173893761566697648044831140121719610701918"), // g(3.0625) + lexical_cast("0.999999999999999936717928653168317632857497076013739688631421486566501253923588291792312872556197791673211"), // g(3.1875) + lexical_cast("0.999999999999999999604353866053235217151851825200411035151346544106038020511604101850368937090131102545215"), // g(3.3125) + lexical_cast("0.999999999999999999998739024155283451527384988865470705871154114370361554458335844183866997811666616798976"), // g(3.4375) + lexical_cast("0.999999999999999999999998127507825472391716958201647120368693798321655943245592754535450849909074747681960"), // g(3.5625) + lexical_cast("0.999999999999999999999999998829969706086880918982632682814675565650518550599127288905129036425692254822915"), // g(3.6875) + lexical_cast("0.999999999999999999999999999999725901382164875904751900132567664126424404651014456483941344862882672111885"), // g(3.8125) + lexical_cast("0.999999999999999999999999999999999978877172785523880230046108100497343478408080122636562274664462699186454"), // g(3.9375) + lexical_cast("0.999999999999999999999999999999999999999538277605174049828569616310243872135620815984308962291209405831177"), // g(4.0625) + lexical_cast("0.999999999999999999999999999999999999999999997579637802373759923392302353674547296197097139608305096662812"), // g(4.1875) + lexical_cast("0.999999999999999999999999999999999999999999999999997484438136198054797417578855545676097374818698933178207"), // g(4.3125) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999582136093029780596306092017648999137976557361095"), // g(4.4375) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999991310181991788327697443705660229702848592"), // g(4.5625) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999982845642234997521097056816153005"), // g(4.6875) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999999999999997650696758780466193299"), // g(4.8125) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999984354567611"), // g(4.9375) + }, + { + lexical_cast("0.049055967305077886314518733812558024186190594522577929133874258403716720136373526724058314689500719879108"), // g(0.03125) + lexical_cast("0.146417984290587940527215566547532640091887719157378510760862597160477250510941162865658587391118599123435"), // g(0.09375) + lexical_cast("0.241566319538883658379307858792962251533817550555192887394117855391242655014383231530195532341198725923960"), // g(0.15625) + lexical_cast("0.333142264577638092438479117881900415393269597145667388229881859612778271640892458711999039673113782322759"), // g(0.21875) + lexical_cast("0.419952111278447158491989828248472390262100664047003718717660961186725980035992633816314430865219256697879"), // g(0.28125) + lexical_cast("0.501013389379309101516506181245109658303016252588117809314943148482857730037964859991354653353614312466427"), // g(0.34375) + lexical_cast("0.575584490635151659953583015206098981937114190718398659114692582600533608360562365795951326255435451324348"), // g(0.40625) + lexical_cast("0.643176758985204701281418660658073036744286341686133235184095903629200246787480818209575639010204336375753"), // g(0.46875) + lexical_cast("0.703550005147142015655495685290319504959682522849123919858984809033031483935095634544309955308190291479053"), // g(0.53125) + lexical_cast("0.756693908633729949407188822489216751981287087446937619321916526900373368852010370466818995015863222862122"), // g(0.59375) + lexical_cast("0.802798741343241265763884181796827847819974819581943422154209276681599958254661449550227532122613671982280"), // g(0.65625) + lexical_cast("0.842219246350756863819973516488706219419512640440301688556941320851442821117057183898770374058443380334815"), // g(0.71875) + lexical_cast("0.875435397630408678372299833411464372063507828193182139899686484625077152372266959771742362259418939591530"), // g(0.78125) + lexical_cast("0.903013281513573870636377959712083386985779993156177959755760035383594055905622749089460693896406222693699"), // g(0.84375) + lexical_cast("0.925568634068612666452141988294418409321506208703172424613590831530263626380800105774028929974373962947218"), // g(0.90625) + lexical_cast("0.943734786052757156854901123618049687883653784322084165689218388450236355243461010301711218796114310673224"), // g(0.96875) + lexical_cast("0.958136022710213690118176450267088471838797224811226203107796332098153155562176068894372980605720885214409"), // g(1.03125) + lexical_cast("0.969366732896917335170154421885106935309845588596038873516433608917130176242697220090959634923155658205546"), // g(1.09375) + lexical_cast("0.977976235186664972979677612058492152933862943354101932731649113297325348962487336855509900187272370184059"), // g(1.15625) + lexical_cast("0.984458831167430830868210653713099901261389062270985127324761800738960911996834155126856635359428165185540"), // g(1.21875) + lexical_cast("0.989248431090133896006986599116918357283502325600362842139835217624907528641466219844620505200692818625881"), // g(1.28125) + lexical_cast("0.992716997196827285379385026669772451664282781030200524124855111594391527718191193755969172969469363705360"), // g(1.34375) + lexical_cast("0.995176026155327354256389972821064571954433337153651487314979185050603642051679069145248180162524268289931"), // g(1.40625) + lexical_cast("0.996880318128191873720647564584261150768550013433900442031828095968019716148998166040133539693795556295480"), // g(1.46875) + lexical_cast("0.998033336315433754022139892013135001514338132834733651914009933714605897127619531730183627753259379942317"), // g(1.53125) + lexical_cast("0.998793534298805899288730607295838650390005311078193298102275453258885228928615060342059178900126999186133"), // g(1.59375) + lexical_cast("0.999281111921791955408551352736963049816503284164678230635189418868597949348817544441358325059187746499100"), // g(1.65625) + lexical_cast("0.999584750351517587316853184319625485667328434389759918903046804784972056092828510019746746746099755782323"), // g(1.71875) + lexical_cast("0.999767971599560835061275314865898944608406416204156618176675346892260692302166389619017568238095126086076"), // g(1.78125) + lexical_cast("0.999874865048780346484697959952056776461887112585682055798709546559427779995450728973085421270772880505664"), // g(1.84375) + lexical_cast("0.999935019925082423686205613430833558155717950000316332442819539848274817113103996480193599096124843935465"), // g(1.90625) + lexical_cast("0.999967593067943459764303089027105328744793599367129242153214575272063224646639331393439819088793577030274"), // g(1.96875) + lexical_cast("0.999984519902270824421965821319650857202000814040840055434710368185052033517924904839634456671486761306484"), // g(2.03125) + lexical_cast("0.999992937876662885647744209318649390050089014936069507171850876567467797652614071870719947353272182519293"), // g(2.09375) + lexical_cast("0.999996932449190357505753493376348767025420602209514502822612289757719467839648455475581472679251593842557"), // g(2.15625) + lexical_cast("0.999998735471865909541730060108096603068366020733848393916677540220757131626396028122330000940209475162201"), // g(2.21875) + lexical_cast("0.999999507005719436888839381971237783832643160123064068160251469835451176020601472504727788247634444794293"), // g(2.28125) + lexical_cast("0.999999818893712767008630838021784833984656997500195308224445639235892462209935570086496911533071977087455"), // g(2.34375) + lexical_cast("0.999999937554078373775024439238869485338620595967968725964389270831692664575370940573486036255773721533551"), // g(2.40625) + lexical_cast("0.999999979874503201754204332281602613820838916879252472683936726944720799255042026941497070536352059158539"), // g(2.46875) + lexical_cast("0.999999993964134201647723624042684414011786944695370213696965260099841383366326401134942210197432768340451"), // g(2.53125) + lexical_cast("0.999999998323361948257805712376050755873763138258297382903914749049017157703250731548370050786872441191487"), // g(2.59375) + lexical_cast("0.999999999570787772613611936445673357498715599182514952506697906806381724508207891577721637577534975420184"), // g(2.65625) + lexical_cast("0.999999999899277723256269151505430020254489776931030883431607669632248209663928689307481753766353500226849"), // g(2.71875) + lexical_cast("0.999999999978455337410957976346300673706685619977882519466522213739703034694765144082155706700617983978517"), // g(2.78125) + lexical_cast("0.999999999995824606875803486698674153541612557781649716847358732403462135338986320491033101425017651823568"), // g(2.84375) + lexical_cast("0.999999999999271526273571226595753644459303012290016807228590124651719266142241023582760851504364325764823"), // g(2.90625) + lexical_cast("0.999999999999886361301456078084894772996885909329652818311789809343152357798533275803234141997050908281812"), // g(2.96875) + lexical_cast("0.999999999999984264464902073469815143022235852317886335270798347860649303075014866860982924784797280179393"), // g(3.03125) + lexical_cast("0.999999999999998080781090289237860478050517213025811746604534278515428437494789627187697928612861722915267"), // g(3.09375) + lexical_cast("0.999999999999999795504045915326443649480437826477882085241327238988393113454004029399049044413682594307796"), // g(3.15625) + lexical_cast("0.999999999999999981130416932706590317890609617287796581196117972272916774467682104728121565361533718859654"), // g(3.21875) + lexical_cast("0.999999999999999998506128351164277103355072159575292493569012746356125826505649217249611405324750888246499"), // g(3.28125) + lexical_cast("0.999999999999999999899530618999624235194896915855549164249906264207980669354314670982085968695886664710980"), // g(3.34375) + lexical_cast("0.999999999999999999994320070821926803299871482831640759432671719063485660287245955331606927623793629089590"), // g(3.40625) + lexical_cast("0.999999999999999999999733089604738648485948153081402766331866202125433024422683657580593309236642932034925"), // g(3.46875) + lexical_cast("0.999999999999999999999989698218618105885316519484716600003641424604944834081762386786826575300824223422458"), // g(3.53125) + lexical_cast("0.999999999999999999999999677551512360875379497532361882754873173523099213379021741376314894561339098393438"), // g(3.59375) + lexical_cast("0.999999999999999999999999991925217099578311319488639767077882000766163827207953449215216503897022509887221"), // g(3.65625) + lexical_cast("0.999999999999999999999999999840534539800671445292714948079656583742156376567760349074494710054769476096983"), // g(3.71875) + lexical_cast("0.999999999999999999999999999997554276024724717412084127804639383434763009476144753753992477175065716930923"), // g(3.78125) + lexical_cast("0.999999999999999999999999999999971340802280921344318626319578389964342657715312145390178385524991742710395"), // g(3.84375) + lexical_cast("0.999999999999999999999999999999999747831747481710059017376411688776331073718895234952452078955161356666502"), // g(3.90625) + lexical_cast("0.999999999999999999999999999999999998364485559689927255326065859921103684150090149812857583677798540583655"), // g(3.96875) + lexical_cast("0.999999999999999999999999999999999999992333340900932164277332575751590554219016841681674832012390635791863"), // g(4.03125) + lexical_cast("0.999999999999999999999999999999999999999974564228711862306525042342323570536582349887001249681816795247227"), // g(4.09375) + lexical_cast("0.999999999999999999999999999999999999999999941590527719283517683183762803793531153548442221910149304532843"), // g(4.15625) + lexical_cast("0.999999999999999999999999999999999999999999999909341304352060297447888980553820078655558718947594338854856"), // g(4.21875) + lexical_cast("0.999999999999999999999999999999999999999999999999907264301839411474095172332986675686118842325985300288496"), // g(4.28125) + lexical_cast("0.999999999999999999999999999999999999999999999999999939142851957947722395964533012630701610319864700124513"), // g(4.34375) + lexical_cast("0.999999999999999999999999999999999999999999999999999999975101910263277853393767914602374615667211494447848"), // g(4.40625) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999993840137913682945570996670024098102535164354421"), // g(4.46875) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999107858729277732880967814578246516778492447"), // g(4.53125) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999926927756739959548610221950520326960748"), // g(4.59375) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999996737439137448614151954888046024572"), // g(4.65625) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999999923642148550163775039117320280"), // g(4.71875) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999999999999101448253097118029843742"), // g(4.78125) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999999999999999994914128148451271000"), // g(4.84375) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999986792603998015"), // g(4.90625) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999985035135"), // g(4.96875) + }, + { + lexical_cast("0.024539763574649160378815204133417875445863215064464608956895044609613165807743160719464451280015177366047"), // g(0.015625) + lexical_cast("0.073525122985671294475493956399705178658389065355171533915990397969116950529498708276835699551602576315925"), // g(0.046875) + lexical_cast("0.122229122201557642351355436474843080727525559168851961307326117592014044180619783289732445222070567026224"), // g(0.078125) + lexical_cast("0.170467972382010518106974587954629422247819259678386906532369539462580024173199368300907551727925894826348"), // g(0.109375) + lexical_cast("0.218063473469712004630197728122759489993591721911968349772635335605013554251005227751602273788509068827768"), // g(0.140625) + lexical_cast("0.264845076583447950461212665118686193377689138924463915604546680494327903082237086492686049952455516196338"), // g(0.171875) + lexical_cast("0.310651780552845960831223570228276117790900830291233480329114422969341275733912045874813905619598144755182"), // g(0.203125) + lexical_cast("0.355333825165074533298754211585183943731716978542867342975713539314813074990332709426902766329500035280122"), // g(0.234375) + lexical_cast("0.398754150467237756442581972104564984183196589603101833773161974728621743001028869891654015488955372569429"), // g(0.265625) + lexical_cast("0.440789599033900866267280244194151053047579390745283782307335001757946437605373983919088061197784856983562"), // g(0.296875) + lexical_cast("0.481331846116905044218492487066532263691894294558410034310919914283949473453346449800853248305213115832126"), // g(0.328125) + lexical_cast("0.520288050691230159575623043642628348233377550276259011535573556712044275760300024941164324417881655260915"), // g(0.359375) + lexical_cast("0.557581228260778230800143467991427986140509811221120710295524345563561959837957800378428132350456008959903"), // g(0.390625) + lexical_cast("0.593150353591953158798252310382019052720472512461466593316789052867105150353197901884985062826170286137806"), // g(0.421875) + lexical_cast("0.626950208051042879496190026689355520279110883033792384384658002934684283337790046078152714136007556881568"), // g(0.453125) + lexical_cast("0.658950991743350124383142853856164010931375645572951466745863593066519424002725927426397714518354148617487"), // g(0.484375) + lexical_cast("0.689137725061667671760526624252008601106750026097972108206033999197424406440192350815404478184554705206257"), // g(0.515625) + lexical_cast("0.717509467487324127212121983621397955277712107841604834988436858336961794561995678665405386180061291186569"), // g(0.546875) + lexical_cast("0.744078383547347399125495563474695309784414124523934238477669187901457818018589261965516293218009220524845"), // g(0.578125) + lexical_cast("0.768868686768246584594556881359924616505401561118656571885272932843285345273991987917365780176625871156636"), // g(0.609375) + lexical_cast("0.791915492376142114474766097125173136071887123795151582468762215511693180036774136547014378362649358506442"), // g(0.640625) + lexical_cast("0.813263608502973851675287583639738347628286388868501267508240391684309992119024653284842237788083864933246"), // g(0.671875) + lexical_cast("0.832966293919410875636012972278592969146439177629268660471223900662097982688700327775233447671477333195295"), // g(0.703125) + lexical_cast("0.851084007987848732614307971738962906523440935719961053106347532970333379131566141191856770896987167106777"), // g(0.734375) + lexical_cast("0.867683175775645986694538915759513734859817770240354328433422368786978718497346146495097167087337333124758"), // g(0.765625) + lexical_cast("0.882834988244668955133093021086816913091662889815942272008298349735524564954895061176615551347722890044403"), // g(0.796875) + lexical_cast("0.896614254280076025787104770742455766232755883620157918788379141962473985621702317217391357385305459551381"), // g(0.828125) + lexical_cast("0.909098318163020435111602266942887955960342926100384359781951030650288122540773234320029408541951937094189"), // g(0.859375) + lexical_cast("0.920366053031952802346959831658231750041733952809073468247631516399787578082498442883770474797966683886924"), // g(0.890625) + lexical_cast("0.930496937997153406311385871717394461677024716743194227098506460563969756852293181618838727220737819009115"), // g(0.921875) + lexical_cast("0.939570223933274755385079922262738480317128204000199829408208173148871526215930675137729398813770704691999"), // g(0.953125) + lexical_cast("0.947664190615153097336932443870647421984719723532535015103119532834778799958165286361562789931722580329484"), // g(0.984375) + lexical_cast("0.954855495805022685406798767343986404221087726922756308613895534964635510152245670365044944835963277500194"), // g(1.015625) + lexical_cast("0.961218615151116407531539792833464814239348141424012723447122892828778812261555210082509575266985875043571"), // g(1.046875) + lexical_cast("0.966825370312355852840021670595606120902775862752839967214604078688367279975132732137152201052503164323018"), // g(1.078125) + lexical_cast("0.971744541565487308923313067648354430256221182743941226733519714630576013967971193573588079157358687735840"), // g(1.109375) + lexical_cast("0.976041560256576739334153026648223010102795289497865077752379571913515848644516440075766293165165274037803"), // g(1.140625) + lexical_cast("0.979778275800615762646502362632280397298673568492270916950415712300131176717110184126189958086261658149965"), // g(1.171875) + lexical_cast("0.983012791481101105576939675455149880687921814247707579509332424119188076725848739069430625018819367499986"), // g(1.203125) + lexical_cast("0.985799363025283435967696415480740550667241416811113828831880235871676397583466503239273324760169219777944"), // g(1.234375) + lexical_cast("0.988188353800742642430716734601052207317333451871274167575722174936569321060220884820380390317458363942270"), // g(1.265625) + lexical_cast("0.990226240467527746936097513319723273560804970370770863320681507664240065241689296263176408382151118320659"), // g(1.296875) + lexical_cast("0.991955663002677615620130768103503643638358675579998145245637497789115454511727184261416776236597257312777"), // g(1.328125) + lexical_cast("0.993415513169264039001291394207886186471136662643955148522315076542909342031246524239981262598429713013047"), // g(1.359375) + lexical_cast("0.994641055712511196722084841575214872266652321301838125470334023086563238065417463278588259848778772721434"), // g(1.390625) + lexical_cast("0.995664076816953169647824835643592605604255578367091325574986552147312999719032619493425800507480569723560"), // g(1.421875) + lexical_cast("0.996513054640253773172879441381629316940644119753128535778483168093423939763922554287566472555643564905551"), // g(1.453125) + lexical_cast("0.997213347043468702236798438923979253284680487723457727227735693623081416988063542707594620320337843765442"), // g(1.484375) + lexical_cast("0.997787391958906530826502855641535679172961288898190898885885789727418242099600873694429770204632167412480"), // g(1.515625) + lexical_cast("0.998254916171996293438409660013903137634975936104914528844970997166876470406307475408661060280450091324018"), // g(1.546875) + lexical_cast("0.998633148640677477622875993802159463019262067333395110056756095635094759421279931934684219243871401865800"), // g(1.578125) + lexical_cast("0.998937034833512173726489175890272462064922637031812443679145922041778828159670660697034088108187551576730"), // g(1.609375) + lexical_cast("0.999179448934885917157557276520126160472957962669292385535106281773260643318307883121231522925081031532969"), // g(1.640625) + lexical_cast("0.999371401140937686896420594457834348577291750383900316575762289527695901511569897912559434414893319746842"), // g(1.671875) + lexical_cast("0.999522237651217204224111367085105825522012295508224046633384623939059543388641920420401673403377094830154"), // g(1.703125) + lexical_cast("0.999639831345600365190217316559820536015055047356781384040847519699899857413145027408233548259988542137061"), // g(1.734375) + lexical_cast("0.999730761519808482630401199111790997353137359325798926522030129197880231229184096116317466868886230470463"), // g(1.765625) + lexical_cast("0.999800481431138386301967810550116384764787615213991040322156393795144469664898240627371239799899054160244"), // g(1.796875) + lexical_cast("0.999853472773111411714825359948522723059543925767081332723268305717111157129995937235104347810321857624981"), // g(1.828125) + lexical_cast("0.999893386547592564256491398103345365738857546255238856501342870780996234580137339944467621762519392320071"), // g(1.859375) + lexical_cast("0.999923170129289328694327906747469096657836526675953732330689707054096435980605112725360972833623996484720"), // g(1.890625) + lexical_cast("0.999945180614458693089605653812338969204804898674074691045932524591547444523053449214931929043688991926743"), // g(1.921875) + lexical_cast("0.999961284807856666131592875140322654964386741993096829332438202532772537452712562987164115971176607273760"), // g(1.953125) + lexical_cast("0.999972946425232236564296551405861928488393447001512979747749895668450972654972928812395117466408027282956"), // g(1.984375) + lexical_cast("0.999981301270120726791718785112613848672595112549858429096586200250084136463600498980061577127305570650369"), // g(2.015625) + lexical_cast("0.999987221282000628110403125088942275790835744576020157472242570908368305740261044660399904176923807917287"), // g(2.046875) + lexical_cast("0.999991368448344873444076514459820294025024291982800024550356196252113148423079976267422555486726944695860"), // g(2.078125) + lexical_cast("0.999994239627616634782013358791490234898949013166028190314689279751288516674766732138410087580437926055702"), // g(2.109375) + lexical_cast("0.999996203347166176753078412505608961497047570078602026635723094411347050066445365758029401697256825964352"), // g(2.140625) + lexical_cast("0.999997529623805167930325832978177062792864667554571428047348632997366248366870090253270772186902244457523"), // g(2.171875) + lexical_cast("0.999998413810964735424197651428792024270859125490162777774589619201616842080978402241545971123088802465345"), // g(2.203125) + lexical_cast("0.999998995410689969623491570252751514540790446083209277324974472698210419125987174633572269691681778748388"), // g(2.234375) + lexical_cast("0.999999372707335369470408065328269510622000820195014514782586629806822913123694905791529274958727999575234"), // g(2.265625) + lexical_cast("0.999999613988550242748977111585844080493329229887891295148587379622015596112918102999514710365309084931588"), // g(2.296875) + lexical_cast("0.999999766023332433120478826363329286107295963668722188789077419686879966508320322753796458956840046022307"), // g(2.328125) + lexical_cast("0.999999860371214599407711186401370807923696229197684810839582879648163434310478325293997656363711618300901"), // g(2.359375) + lexical_cast("0.999999918004794710562251438113801568778259940436067788001646433851504120721491427662620358315670934848881"), // g(2.390625) + lexical_cast("0.999999952642664461848026000044314367945034834950818691991040100849072036057628998371722435136215664109513"), // g(2.421875) + lexical_cast("0.999999973113235943623083169854869397745128002613407767216515340339776540129298040295912524593924112195344"), // g(2.453125) + lexical_cast("0.999999985003076311733367548176404664439498573729865861609501261728514520871338200270754248548994865324468"), // g(2.484375) + lexical_cast("0.999999991786456099071980956791533766280975098545545439443600528014542930211812970251850407150921831151341"), // g(2.515625) + lexical_cast("0.999999995585633615838753599126475596488193895929739538678756986286051775589140326875249634732087708344684"), // g(2.546875) + lexical_cast("0.999999997673236737899538770171294615029525379682022727961052725153761673551869149081088514655528112615723"), // g(2.578125) + lexical_cast("0.999999998797983500398859233110536764723370506014230505733946581769681293073500199066363731441981179698668"), // g(2.609375) + lexical_cast("0.999999999391776875832213848946106251733222623406697155203862827599511530548696684554690294395044460402188"), // g(2.640625) + lexical_cast("0.999999999698754369251676774371962641972303016749904290546063818997066776513061336189682166498719398944505"), // g(2.671875) + lexical_cast("0.999999999854056115499871428560854000803460484952606218446207511740848374699327458672845817325452074300416"), // g(2.703125) + lexical_cast("0.999999999930888395014768921333032751771561462357776403267387698410910157284722533983890902013975177726611"), // g(2.734375) + lexical_cast("0.999999999968033216735554630277149570764721028155009147038916015563025328071758823521563015242651537903406"), // g(2.765625) + lexical_cast("0.999999999985568790075957347505492654540535325180712834768216408753204137307652756418229700527473627968818"), // g(2.796875) + lexical_cast("0.999999999993646323871508818572328398040876268802630587696310772380050296314289284088827434267491345072485"), // g(2.828125) + lexical_cast("0.999999999997274049481127480279890956867731032250306826061026076773018983692894231844753659546719950003982"), // g(2.859375) + lexical_cast("0.999999999998861265427759215034028892835419035064005497119819632748241441455045548565733019611100832082288"), // g(2.890625) + lexical_cast("0.999999999999537226537823068657094388275825668825768255803239377404791245063232344852624997490140016897718"), // g(2.921875) + lexical_cast("0.999999999999817200977515688951599009732342304725120189233537910012424069109148211340137931293203647917330"), // g(2.953125) + lexical_cast("0.999999999999929879532615094313912814301197441877086628318390846930022940168414376904920853375432682687685"), // g(2.984375) + lexical_cast("0.999999999999973903946338620453625381179011373065415854711612340117275559406451722886433991855180609703337"), // g(3.015625) + lexical_cast("0.999999999999990586638664711144423931453473543946375173168784428229755021554384267101399389484316374784633"), // g(3.046875) + lexical_cast("0.999999999999996712074305063506579609420494394694814271462713124900358985684585644928436634392284769770366"), // g(3.078125) + lexical_cast("0.999999999999998889137059817810156874902680103568441605068555530814349981256386951911742679432410358109363"), // g(3.109375) + lexical_cast("0.999999999999999637339773579678779802821695809371718632063128391240509853494297468740705558582565789306447"), // g(3.140625) + lexical_cast("0.999999999999999885721234749473301367201222927198470102374965213685001124129999378829702104389172344477206"), // g(3.171875) + lexical_cast("0.999999999999999965280977307505063204075886852395058039071502137488669959645062303252728871472589128115277"), // g(3.203125) + lexical_cast("0.999999999999999989842180931694875412550702445624143893055015768061161465215549282742511675317409286414989"), // g(3.234375) + lexical_cast("0.999999999999999997141467175491851826072519244255165504448953132277060372277893588237238102438615305399735"), // g(3.265625) + lexical_cast("0.999999999999999999227216521261074474315562149281863083927192433309338485624907040900104468036863405543457"), // g(3.296875) + lexical_cast("0.999999999999999999799557600778248277950206850471623227374381187035560338059265635734697070049409665738103"), // g(3.328125) + lexical_cast("0.999999999999999999950184316235300885068961773562422713606138181060098612187655205829612361585819533948919"), // g(3.359375) + lexical_cast("0.999999999999999999988153315075776597632355126174005520227700785943127093556947602824089332125864400129929"), // g(3.390625) + lexical_cast("0.999999999999999999997308015096389609429298796861591756550312366158425103234132330040727377287784163866254"), // g(3.421875) + lexical_cast("0.999999999999999999999416333255806605858701058702611105194050020387210497348612151663233096746303243479768"), // g(3.453125) + lexical_cast("0.999999999999999999999879433841134512764866456281397142486130414885448482190806795767003151711453109127778"), // g(3.484375) + lexical_cast("0.999999999999999999999976308906491921932349406437681865603051313874505442746511506957918206402657708420626"), // g(3.515625) + lexical_cast("0.999999999999999999999995578660538271998224864063183274180710327093731489309148308151208754306245398458519"), // g(3.546875) + lexical_cast("0.999999999999999999999999217616499589131418178103885675337588998727790164441016998178184320629717468797985"), // g(3.578125) + lexical_cast("0.999999999999999999999999868946685587452184514479879444808787116938931736195336569126670158288153262213351"), // g(3.609375) + lexical_cast("0.999999999999999999999999979256549869264849260037256635499721719365364384464776488808336965445974843567748"), // g(3.640625) + lexical_cast("0.999999999999999999999999996903031674449250417163036747168003461618781879063247701790152002357385598623179"), // g(3.671875) + lexical_cast("0.999999999999999999999999999564679947180324402422599346898721478673555226731002774524230012051141383347060"), // g(3.703125) + lexical_cast("0.999999999999999999999999999942500443319296077099306861248145445450795420368313946552211166150789384747398"), // g(3.734375) + lexical_cast("0.999999999999999999999999999992877284073361827258793136861565431252575172626145737925177835263625060324817"), // g(3.765625) + lexical_cast("0.999999999999999999999999999999174216457964419542319955486789754785939078150605472111982162098188134509617"), // g(3.796875) + lexical_cast("0.999999999999999999999999999999910584589933750444562352062330033311295140225751179436546909270146329553775"), // g(3.828125) + lexical_cast("0.999999999999999999999999999999990977203349969133025734556703837225049292006409995763226680819494755149947"), // g(3.859375) + lexical_cast("0.999999999999999999999999999999999153396988046840282473113845798648089047071875648546359194078420443211209"), // g(3.890625) + lexical_cast("0.999999999999999999999999999999999926307271934354701831382163381875333775364498721513748568754428840374083"), // g(3.921875) + lexical_cast("0.999999999999999999999999999999999994063367643991125022372144215977924493444982440970666750261770833501016"), // g(3.953125) + lexical_cast("0.999999999999999999999999999999999999558472216130530404546187628565445984008004315791155352111445810861783"), // g(3.984375) + lexical_cast("0.999999999999999999999999999999999999969760398641364685083035283418365832057343804122335170816451269504560"), // g(4.015625) + lexical_cast("0.999999999999999999999999999999999999998097796208342295120705633213220124161225158854624041662628043139974"), // g(4.046875) + lexical_cast("0.999999999999999999999999999999999999999890395667139654221174510368013567712749108748450937018733493364962"), // g(4.078125) + lexical_cast("0.999999999999999999999999999999999999999994231304401995930517267920264471971517126429270010022117510177093"), // g(4.109375) + lexical_cast("0.999999999999999999999999999999999999999999723460061826228977881120567397993972071695951028748597577842147"), // g(4.140625) + lexical_cast("0.999999999999999999999999999999999999999999987961290471437425948891541215692215373303586235782453845147660"), // g(4.171875) + lexical_cast("0.999999999999999999999999999999999999999999999525523767466450128221678490112046052186534432663346484517441"), // g(4.203125) + lexical_cast("0.999999999999999999999999999999999999999999999983123208718946690858956148015221090678809747386260398948125"), // g(4.234375) + lexical_cast("0.999999999999999999999999999999999999999999999999460003300170027151912938166349753544319699832203992897436"), // g(4.265625) + lexical_cast("0.999999999999999999999999999999999999999999999999984509759505519173856534930303705769734720891642360192687"), // g(4.296875) + lexical_cast("0.999999999999999999999999999999999999999999999999999603005274474386502345348727988708714041059894095275250"), // g(4.328125) + lexical_cast("0.999999999999999999999999999999999999999999999999999990942325073476064144744046894900843568862749908403913"), // g(4.359375) + lexical_cast("0.999999999999999999999999999999999999999999999999999999816704956651630686192426467157780890405204216080278"), // g(4.390625) + lexical_cast("0.999999999999999999999999999999999999999999999999999999996722577667602635610322309320474579393404987417070"), // g(4.421875) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999948423128469160636359110852375729788440977714938"), // g(4.453125) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999288528477859180200000678913998253560344994272"), // g(4.484375) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999991433110982878975067646676928663510925470880"), // g(4.515625) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999910344414457886381411766229763114044559702"), // g(4.546875) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999188130620309452561947632987885842177083"), // g(4.578125) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999993667806114092266084417443887109377185"), // g(4.609375) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999957662736575396719127741239143326530"), // g(4.640625) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999758527395356232098031548653832407"), // g(4.671875) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999998831044408905178547422145308590"), // g(4.703125) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999999995221933101447355392021268151"), // g(4.734375) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999999999983597964283047841009690344"), // g(4.765625) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999999999999952975108512357894918527"), // g(4.796875) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999999999999999888040953299220192760"), // g(4.828125) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999779945664926699750"), // g(4.859375) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999645107597593674"), // g(4.890625) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999533305906312"), // g(4.921875) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999502809333"), // g(4.953125) + lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999573749"), // g(4.984375) + }, + }; + + const std::vector> m_weights{ + { + lexical_cast("1.570796326794896619231321691639751442098584699687552910487472296153908203143104499314017412671058533991074"), // g'(0) + lexical_cast("0.230022394514788685000412470422321665303851255802659059205889049267344079034811721955914622224110769925389"), // g'(1) + lexical_cast("0.000266200513752716908657010159372233158103339260303474890448151422406465941700219597184051859505048039379"), // g'(2) + lexical_cast("0.000000000001358178427453909083422196787474500211182703205221379182701148467473091863958082265061102415190"), // g'(3) + lexical_cast("0.000000000000000000000000000000000010017416784066252963809895613167040078319571113599666377944404151233916"), // g'(4) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002676308"), // g'(5) + }, + { + lexical_cast("0.965976579412301148012086924538029475282925173953839319280640651228612016942374156051084481340637789686773"), // g'(0.5) + lexical_cast("0.018343166989927842087331266912053799182780844891859123704139702537001454134135727608312038892655885289502"), // g'(1.5) + lexical_cast("0.000000214312045569430393576972333072321177878392994404158296872167216420137367626015660887389125958297359"), // g'(2.5) + lexical_cast("0.000000000000000000002800315101977588958258001699217015336310581249269449114860803391121477177123095973212"), // g'(3.5) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000011232705345486918789827474356787339538750684404"), // g'(4.5) + }, + { + lexical_cast("1.389614759247256322860819129532051320578963754203760794373378106104353527766753015814332241420248231414085"), // g'(0.25) + lexical_cast("0.531078275428053974761132317615314081316383438040798655367002821918297730288693432971141471528424730268920"), // g'(0.75) + lexical_cast("0.076385743570832304188340574831642800088057531968167335207418652655026443958989111001740359820912513366775"), // g'(1.25) + lexical_cast("0.002902517747901313593593294890458021493283931520621269849938702028801100032945565619776694568712365440204"), // g'(1.75) + lexical_cast("0.000011983701363170720046901264217342609850414053186881688297911977061939129002941852260994641452968753149"), // g'(2.25) + lexical_cast("0.000000001163116581425578276559715526238292550912119341858913103096243311657573123490443124099153951208264"), // g'(2.75) + lexical_cast("0.000000000000000219707923629797991740920411247835398364032243317196139216771669502867834292155305036890622"), // g'(3.25) + lexical_cast("0.000000000000000000000000001363510330763761541372474765508158531409135540695703430420767286724584829502841"), // g'(3.75) + lexical_cast("0.000000000000000000000000000000000000000000000337005685404192649899341739285505662357053600685059902409156"), // g'(4.25) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000000051969783800898552138641449339"), // g'(4.75) + }, + { + lexical_cast("1.523283718634705213194962790158845334535887058581160414232292195083923441711089834872202018481714068265874"), // g'(0.125) + lexical_cast("1.193463025849156963909371648222971962983135675758906039639808708295899934988800331391641633505274800733739"), // g'(0.375) + lexical_cast("0.737437848361547841364500858485266812076333434271679688817118172760237363857747855352433500512107413836274"), // g'(0.625) + lexical_cast("0.360461418469343674165419405305111915999562236843279741777216862896838317598331490355112158453382264499935"), // g'(0.875) + lexical_cast("0.137422107733167723411102816000757634955546811302944952911944638642075130409924435566222978577348630239795"), // g'(1.125) + lexical_cast("0.039175005493600779071814125724013543926474361440430503895786595281886103593125778916562474039476045546479"), // g'(1.375) + lexical_cast("0.007742601026064240712330911164066815715739004138914557972278805220602461623801176256818815341175689849120"), // g'(1.625) + lexical_cast("0.000949946804283468716905391803588290645599496466362504375579729890166137558846764281909082013792857588283"), // g'(1.875) + lexical_cast("0.000062482559240744082890784437584871091195765427233728160580340356687494550913554122518585235582888183362"), // g'(2.125) + lexical_cast("0.000001826332059371065969910928097449472730027834534513400058702156695265420575321758936306719563351875792"), // g'(2.375) + lexical_cast("0.000000018687282268736410131523743935312466546469465751985818253470776559331427674118404092705647447708478"), // g'(2.625) + lexical_cast("0.000000000049378538776631926963708240981686385688538089122379464179611658093414201820783209158947548099290"), // g'(2.875) + lexical_cast("0.000000000000022834926702613953995564216678996857823700190922893992062567474697649971246018412201710541049"), // g'(3.125) + lexical_cast("0.000000000000000001122753142818155150094255452340294487178754520669590337027693465434158001687151621567161"), // g'(3.375) + lexical_cast("0.000000000000000000000003097653970117354371545851432763107826514605577040924324913965609976812540667954525"), // g'(3.625) + lexical_cast("0.000000000000000000000000000000211212334353722559135268119776231624097073988695606211819781446821221407215"), // g'(3.875) + lexical_cast("0.000000000000000000000000000000000000000124241475706160523670077003963444149150044959415404297798776679701"), // g'(4.125) + lexical_cast("0.000000000000000000000000000000000000000000000000000163277073317994932378129476269327467680577634204011231"), // g'(4.375) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000084606887310962137960222762714813370163"), // g'(4.625) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000001864449207958865027"), // g'(4.875) + }, + { + lexical_cast("1.558773355533330145062422255303911650913653419716123258575823844242311218191742707506332260947433822741044"), // g'(0.0625) + lexical_cast("1.466014426716965781027541193665889531415409392722064361086248282756748739496420298831321519828378556002658"), // g'(0.1875) + lexical_cast("1.297475750424977997988538308528435261157132579154910704420759249609048525111522825715836900012713756740793"), // g'(0.3125) + lexical_cast("1.081634985490070407444853274933647081975577566901795789661775590774791605470760834259212179529338772899439"), // g'(0.4375) + lexical_cast("0.850172856456620068952731098052211300592762473235164319319730375930401185061654589540728127181453457995871"), // g'(0.5625) + lexical_cast("0.630405135164743691060150580239202409283620999547554574405345432687011289480732546270043535835546239806306"), // g'(0.6875) + lexical_cast("0.440833236273858237067712706109932220090079370355383513418787238650680549127633879301802645739940335717360"), // g'(0.8125) + lexical_cast("0.290240679312454185000612314799668780962191429526944375904765893963488780983625043853285447755293573461581"), // g'(0.9375) + lexical_cast("0.179324412110728292963459783971217654539874988705586245112221317641182103571339286966223492929524044754476"), // g'(1.0625) + lexical_cast("0.103432154223332900624823850529514176586578503073172859130601674265212005856359181729852654906098803578827"), // g'(1.1875) + lexical_cast("0.055289683742240583845301977440481556933408044706616680391861008544193557784606931443468908103837024579093"), // g'(1.3125) + lexical_cast("0.027133510013712003218708018921405436444642216545174487001029628100858795149505971331120118430786805813066"), // g'(1.4375) + lexical_cast("0.012083543599157953493134951286413130934561539333133237043169739476307291054367023200000385825615221945447"), // g'(1.5625) + lexical_cast("0.004816298143928463017275766007138771453339922288894138500249107155688575955404714471338997856655589497183"), // g'(1.6875) + lexical_cast("0.001690873998142639647215541751024903372105472295812356286827736078638936018893984073130997792093832518523"), // g'(1.8125) + lexical_cast("0.000513393824067903360165889064488588831319697335602551298104262323633326714857215022112175476113997205237"), // g'(1.9375) + lexical_cast("0.000132052341256099748786804020743407578728552770515294022030238142601687785776521371211170470527518902633"), // g'(2.0625) + lexical_cast("0.000028110164327940134748736157546988307382789043257543279202085484174573606187127697055691202877555401179"), // g'(2.1875) + lexical_cast("0.000004823718203261550212402544034394316885029070995428950091523875242090745127519612129898123057387580549"), // g'(2.3125) + lexical_cast("0.000000647775660359297199077339874176974315780234569374097629681497095863536185282276464511564351373607156"), // g'(2.4375) + lexical_cast("0.000000065835185127183396672340995882066948606156286957895280018929061991227529612096505273354107037836464"), // g'(2.5625) + lexical_cast("0.000000004876006097424062586890460642680934720698091346275265980808135812651900878569853770505821434068830"), // g'(2.6875) + lexical_cast("0.000000000252163479185301485718266564918543983391535784401652466437582530827387158493030119120918329405106"), // g'(2.8125) + lexical_cast("0.000000000008675931414979604650195607762477884261062369187782542267591220091953822676501175854014087120773"), // g'(2.9375) + lexical_cast("0.000000000000188020717307506498094762558437719749824538089416626388343837727165608674244286899191331882231"), // g'(3.0625) + lexical_cast("0.000000000000002412423038430878639389930773055029518471728488024078428020750943576498858962639362755793751"), // g'(3.1875) + lexical_cast("0.000000000000000017084532772405701711664118263431986062467708215229851301797632080501043865563725375991858"), // g'(3.3125) + lexical_cast("0.000000000000000000061682568490762382593952567143955272099148473965377186569434997700282162435358180953024"), // g'(3.4375) + lexical_cast("0.000000000000000000000103767972385287061606108632709158267387585366439286051516655090912960812409101940501"), // g'(3.5625) + lexical_cast("0.000000000000000000000000073459841032226935608549262812034877184478626860884665235801772548317888502733967"), // g'(3.6875) + lexical_cast("0.000000000000000000000000000019497833624335174811763249596143101399876957245735428968235842941245604813388"), // g'(3.8125) + lexical_cast("0.000000000000000000000000000000001702438776125754721867532441086720904090474496331484536273502315569284797"), // g'(3.9375) + lexical_cast("0.000000000000000000000000000000000000042164863709484278882204462399998678775113208285972487939873276630037"), // g'(4.0625) + lexical_cast("0.000000000000000000000000000000000000000000250442771162752843178118475337847815354603491905071780281742971"), // g'(4.1875) + lexical_cast("0.000000000000000000000000000000000000000000000000294936014574619331371951595670994262584278418532490648159"), // g'(4.3125) + lexical_cast("0.000000000000000000000000000000000000000000000000000000055513323569653710605229337238942258327417785511979"), // g'(4.4375) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000001308116512021082164368611090789921129262011"), // g'(4.5625) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000002926082446382397532862965945236761"), // g'(4.6875) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000000000000454077346699978241489410"), // g'(4.8125) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003426563569209"), // g'(4.9375) + }, + { + lexical_cast("1.567781431307221857184578395608129985508731044605633255795643429120138259524502701742482715390419743940777"), // g'(0.03125) + lexical_cast("1.543881116176959220412019565230455395681644041043861723595162102954998162906440561996844264719288523984838"), // g'(0.09375) + lexical_cast("1.497226222541036289617512110625370605944473084396164048384849286249507314636821432486108538395745995071173"), // g'(0.15625) + lexical_cast("1.430008354872299667629414512169884537246354149720184463003734364143612515271259946631802092783120723129071"), // g'(0.21875) + lexical_cast("1.345278884766251661463188142158828373672838450871367155499486867249420628703086880438041422396860114207358"), // g'(0.28125) + lexical_cast("1.246701207451857704817137316675678337329938931463972848778333137286136044911173767019357851543515360084940"), // g'(0.34375) + lexical_cast("1.138272243376305373371832830171750926260504280858791638179927315314246259275503532583433377703199393165722"), // g'(0.40625) + lexical_cast("1.024044933111811448259402245494254120254049912197634153284183238679245275199953891255656527115306979942845"), // g'(0.46875) + lexical_cast("0.907879379154895316930979720597322569577291079355158209255092870544770061064288836788563808366846645777866"), // g'(0.53125) + lexical_cast("0.793242700820516717873852598629955893388585169481474291350088672560007276198658983401795105363814513427149"), // g'(0.59375) + lexical_cast("0.683068516344263754641188931872023686010045438275124159828436415019300544929096961591360171014695325357182"), // g'(0.65625) + lexical_cast("0.579678103087787647081104524297701196102396798987896883287943431946792649799001507662982753379788446917319"), // g'(0.71875) + lexical_cast("0.484758091214755392865907754198688861835957473837081540217565397499578160440294356369986977943420575999271"), // g'(0.78125) + lexical_cast("0.399384741525717135146966196552751829485322444839752200305051054663621386574293421127304951583313147231335"), // g'(0.84375) + lexical_cast("0.324082539611528904017760150000608521496610285955450235939434098847767708174500119693298381850542941636110"), // g'(0.90625) + lexical_cast("0.258904639514053516002513872589102127104634909933912066661562881185787066303053216189795980213644036906794"), // g'(0.96875) + lexical_cast("0.203523998858601745186478849132323237332201554075909184023685755157419223830048644915744598116785548292294"), // g'(1.03125) + lexical_cast("0.157326203484366150267385633762597340513018425856811416700022493599965664726952329950890984836925162081776"), // g'(1.09375) + lexical_cast("0.119497411288695924275942759058221440732591462102393865923809745453652225048676424614704499759122418120112"), // g'(1.15625) + lexical_cast("0.089103139240941462840959442033942473840983177593666127837963626878433665232163310822511715363358368012614"), // g'(1.21875) + lexical_cast("0.065155533432536205042255277931864103027470024503861546227085609363824257186720490477503981303197700490281"), // g'(1.28125) + lexical_cast("0.046668208054846613643791300289316396039224003031511647004468288367153705920923668212902762934313474116113"), // g'(1.34375) + lexical_cast("0.032698732726609031112522702100248948990202653053699579743688734952338451372781337540824178631609888378938"), // g'(1.40625) + lexical_cast("0.022379471063648476483258477281403216114463225472686280254823920842131679381849488698099157182778552585696"), // g'(1.46875) + lexical_cast("0.014937835096050129695520628452448084817491418397098138574408056027824742558359036353042202729259689475779"), // g'(1.53125) + lexical_cast("0.009707223739391689269235578630758993667430834815092132317202977202728846588535795538939786744891523508391"), // g'(1.59375) + lexical_cast("0.006130037632083030125244577194087324639280075762202651535474189288959362910416267797977171536000338405928"), // g'(1.65625) + lexical_cast("0.003754250977431834302296714460279130616059681869177316891945506584252205866750408952324976329360048041744"), // g'(1.71875) + lexical_cast("0.002225082706478642702158462041170389618927480048039676593550739331725261004361394313795759779229985688013"), // g'(1.78125) + lexical_cast("0.001273327944708238202674090353557735272857893449667410082070223063619923848468084914029029029452352122423"), // g'(1.84375) + lexical_cast("0.000701859515684242270804743047185673315267810259446867393407417805858312131136536526522872430086830242775"), // g'(1.90625) + lexical_cast("0.000371666936216777603012957602189683547361029702146538431198647394039953879809033420775946520426760365369"), // g'(1.96875) + lexical_cast("0.000188564429767003185715299221155754737487923589299711135953997833183257028113286668727279094268195580509"), // g'(2.03125) + lexical_cast("0.000091390817490710122732277133049597671605200249047028801062145237116184330202586168290278123433996992006"), // g'(2.09375) + lexical_cast("0.000042183183841757600604161049520839394732745162386101627744778040590027419635042422269594227475926856493"), // g'(2.15625) + lexical_cast("0.000018481813599879217116302847163026167343229526622133814885444866635226594107310616201879681803979268598"), // g'(2.21875) + lexical_cast("0.000007659575852520316256156860619950629602221322327668123619654350335655900260321461248258401686383450872"), // g'(2.28125) + lexical_cast("0.000002991661587813878709443954037223499025021052024899534528690163884768913124192316308752952221239364875"), // g'(2.34375) + lexical_cast("0.000001096883512590126473196793140106175059329465735584851612728554325930415320046299752063818350148758993"), // g'(2.40625) + lexical_cast("0.000000375954118623606300911838176681463408070589556116345501748040734638886414965561255409187797947974925"), // g'(2.46875) + lexical_cast("0.000000119924427829027702186799922060706889373331936612294691392478456284269112386616767757113062219199407"), // g'(2.53125) + lexical_cast("0.000000035434777171421953042822362946016498779589683183675332697044799658134169915821826432359139146781942"), // g'(2.59375) + lexical_cast("0.000000009649888896108963360920618196763059172060396135408883397117980894733376845656283202132084956193859"), // g'(2.65625) + lexical_cast("0.000000002409177325647594077859608843933217848949198452356641362987145917749352997079714644193317406807080"), // g'(2.71875) + lexical_cast("0.000000000548283577970949775501245692628858095972575807076726362350725571330454450349783337533439866790293"), // g'(2.78125) + lexical_cast("0.000000000113060553474946805357898794970458309140084494916872577635069922418330209552574406913559269814134"), // g'(2.84375) + lexical_cast("0.000000000020989335404511469108589724179284344077119940052820532945558095739022557242864542615740574013845"), // g'(2.90625) + lexical_cast("0.000000000003484193767026105968529891782697496782392224937035299906089661427166315608771327082526680975547"), // g'(2.96875) + lexical_cast("0.000000000000513412752450142074742769387618083006551339099252922922260923080821580257441827928430048962416"), // g'(3.03125) + lexical_cast("0.000000000000066639922833087653243643099189129361080639834075520364446277112540163227315881566635173210918"), // g'(3.09375) + lexical_cast("0.000000000000007556721775780565189445526183865755658208623494563436730029973698620661370460716469199723406"), // g'(3.15625) + lexical_cast("0.000000000000000742099323099221675885386620710989719376118530280158051428109311268002669669806056042011289"), // g'(3.21875) + lexical_cast("0.000000000000000062528048446104553564782731960667282431028411480030870948856657091252911728874557758278268"), // g'(3.28125) + lexical_cast("0.000000000000000004475759506669096971393956684263666647035805905451937875326762461236603085923164471607374"), // g'(3.34375) + lexical_cast("0.000000000000000000269312066148696950577145258104178242537085822059480713001141057884061799266451750094679"), // g'(3.40625) + lexical_cast("0.000000000000000000013469941569542286092134858963015561957068862655850623081136682314255544860877000996881"), // g'(3.46875) + lexical_cast("0.000000000000000000000553358349941557114556332284896446443290601834217024573566655194255561252181383565972"), // g'(3.53125) + lexical_cast("0.000000000000000000000018435469747181493780961793761587470494572758329229028729042387873382333895611282461"), // g'(3.59375) + lexical_cast("0.000000000000000000000000491393687126490400833082009577485742508307494739377169222065565536809285392995916"), // g'(3.65625) + lexical_cast("0.000000000000000000000000010329391306928575387668383646132905267564979978946790362051244483324690166898666"), // g'(3.71875) + lexical_cast("0.000000000000000000000000000168627700384926065277190131694716102552828133425615722292602284904843076675660"), // g'(3.78125) + lexical_cast("0.000000000000000000000000000002103305749001808952384333210000493449079478161554745190912387376533447476822"), // g'(3.84375) + lexical_cast("0.000000000000000000000000000000019699209796232343253581550830906288944882338601002066557038963464059035686"), // g'(3.90625) + lexical_cast("0.000000000000000000000000000000000135998946163037956979409055971971623477353232403235467218525345012385250"), // g'(3.96875) + lexical_cast("0.000000000000000000000000000000000000678597883755924790855405829688203467120274542406668117110849232223909"), // g'(4.03125) + lexical_cast("0.000000000000000000000000000000000000002396506369944321740606723199153688848649681614787602593437029678159"), // g'(4.09375) + lexical_cast("0.000000000000000000000000000000000000000005857956948308421084638634596113137471738261348961749219470996058"), // g'(4.15625) + lexical_cast("0.000000000000000000000000000000000000000000009678392775571709557136698357096443369145187187297977630530985"), // g'(4.21875) + lexical_cast("0.000000000000000000000000000000000000000000000010538361132564208838364959652125833459207883163470694926441"), // g'(4.28125) + lexical_cast("0.000000000000000000000000000000000000000000000000007361585830978764921078757415438121748497657856592622296"), // g'(4.34375) + lexical_cast("0.000000000000000000000000000000000000000000000000000003205978535283386613396816070688529440452628300680685"), // g'(4.40625) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000844308926618642560822341111358361213153156398923"), // g'(4.46875) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000130166948744281729919957292529774299273116940"), // g'(4.53125) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000011348985048372046635260385887339234534102"), // g'(4.59375) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000539388136951158084102448622955692896"), // g'(4.65625) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000013438019476233711028911671531262"), // g'(4.71875) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000000000168330959366332402630225083"), // g'(4.78125) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000000000000001014205907446815853819"), // g'(4.84375) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002803613619344987"), // g'(4.90625) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003381538789"), // g'(4.96875) + }, + { + lexical_cast("1.570042029279593146749298643773506406814086664704211409620743165557095475459206175917795800288827901535159"), // g'(0.015625) + lexical_cast("1.564021403773232099877114326609256956911418138079007203782656883101047408188541309713005504608886652171494"), // g'(0.046875) + lexical_cast("1.552053169845412119220857160738363356586316176993686065295094432448602562558992253418868577735906380650163"), // g'(0.078125) + lexical_cast("1.534281738154303431635377214537564601260501752273321101306796522290938238587552763019996938086320242027937"), // g'(0.109375) + lexical_cast("1.510919723074169712706179187363428784413953706025910740462077188975653010432297801005140712708726920545807"), // g'(0.140625) + lexical_cast("1.482243297885538069991851555750554483010306477303941134907131825367652175447657674446878021580864448650637"), // g'(0.171875) + lexical_cast("1.448586254961322591610626268497558230170223858747360936806953445653102211642913858887350319280208271657703"), // g'(0.203125) + lexical_cast("1.410332971446259012949044799400294954301898627186633342241168192591777573884752147008141550835945611294496"), // g'(0.234375) + lexical_cast("1.367910511680896488078889666873193284438102768431083960444814279979159678575452637837455834146916421212506"), // g'(0.265625) + lexical_cast("1.321780117443772857894583479184842537571591408745703009579739899003716528422133681645270092231900420827150"), // g'(0.296875) + lexical_cast("1.272428345537862708234621301871357262147717256984072891079400459864079308600214969023660326846987095558686"), // g'(0.328125) + lexical_cast("1.220358109579358220738638953299856684492511943586747141621126612743181016956293456941328702642750163756344"), // g'(0.359375) + lexical_cast("1.166079869932434576632765865091072864305670126415549676173134966612509443875539501185930661534132661935453"), // g'(0.390625) + lexical_cast("1.110103193965340379556876518609851780910441897724068922205543285834036684652901324628135282821996211834341"), // g'(0.421875) + lexical_cast("1.052928879955266655551187700034453673940330202986808761458284686753221730146087158244559719140390280328647"), // g'(0.453125) + lexical_cast("0.995041804046132715136568973654376926257305447497770119916811783573051625756641694631375051739552451554297"), // g'(0.484375) + lexical_cast("0.936904612745667933656663094921637477964830274329132749108110307970948926132400847952022519804849967561534"), // g'(0.515625) + lexical_cast("0.878952345552782120390564821964559257153771691121306682522234582360965715540208272322640867795296879091994"), // g'(0.546875) + lexical_cast("0.821588035266964703341846909903376178823710177814231982112742030164074686054430487677703743009966290959465"), // g'(0.578125) + lexical_cast("0.765179298908956136701604036106714533065059966771117614841446656450553668987266113665955582929621488700271"), // g'(0.609375) + lexical_cast("0.710055901205468983853355833069699896946037929979636279429463515583289789059832935375165659248459758404114"), // g'(0.640625) + lexical_cast("0.656508246131627530756371874069751788892809075707880616945228311304900136436163527616325116228825847364109"), // g'(0.671875) + lexical_cast("0.604786730578403621576636951410117080727209915783816537759975928661680491031798428382325277636413350972804"), // g'(0.703125) + lexical_cast("0.555101878003633509592664715171720141800942797619421185423997564370970754413541571257543660152161513494345"), // g'(0.734375) + lexical_cast("0.507625158831908099697441295086072121614833574539416807394886434591712490606916388594466498892454273251129"), // g'(0.765625) + lexical_cast("0.462490398055367761295640568255341077403273127548291267196344494930532799529885501995939841083288190265935"), // g'(0.796875) + lexical_cast("0.419795668445015480655895454933802487550723118496189647500581075312799046213883645569407253613440101670130"), // g'(0.828125) + lexical_cast("0.379605569386651609992271273702887643026599498568472638870347635328473379544173684130896724340925704191336"), // g'(0.859375) + lexical_cast("0.341953795923016832303703636637002678298132433609185927721103216849472031446233191415425656206133763452558"), // g'(0.890625) + lexical_cast("0.306845909417916949320770025172500171932757625129106531127035036866109901738360194280724087448787883491006"), // g'(0.921875) + lexical_cast("0.274262229689068106370909601214370028469444994893463333836460575013533444158848354755160598383368753799504"), // g'(0.953125) + lexical_cast("0.244160777869839908675206128139678921488760784834565537915351968670913878613710851086025790670715529690192"), // g'(0.984375) + lexical_cast("0.216480209117296170381315245287803516998395334333622786452468494056186127392094125586255533494464103824687"), // g'(1.015625) + lexical_cast("0.191142684133427495322250957465087808288759909097051830915963580310338313269429897603922679940286942202189"), // g'(1.046875) + lexical_cast("0.168056637948269162330314135966437939639159876882688238522888791212549229742151667284607236231343861055819"), // g'(1.078125) + lexical_cast("0.147119413257856932477963110961851225376139575535567887401811129731399585464271565902323262216224330113093"), // g'(1.109375) + lexical_cast("0.128219733631200986752782738812952782200889020485884590258057184849730400243209345821714481205668251977169"), // g'(1.140625) + lexical_cast("0.111239998988744530348740542251160908315985045286433642870454488668697320429550207947265692776313325385783"), // g'(1.171875) + lexical_cast("0.096058391865189467849178662072631058481111936168344706567642860168891311467280932468611840536857385286659"), // g'(1.203125) + lexical_cast("0.082550788110701737653537523573008139539626250354113812806122012801880295654963209039406446997847977585216"), // g'(1.234375) + lexical_cast("0.070592469906866999351840896875301285930871379942130429339219840984873549638194867403335627317258620428011"), // g'(1.265625) + lexical_cast("0.060059642358636300319399035772689020859372782390819237694891224415362476529538938910277541629653752239857"), // g'(1.296875) + lexical_cast("0.050830757572570471070050742610585546538777855107871728732686017515674115769714999939010821344254060136281"), // g'(1.328125) + lexical_cast("0.042787652157725676034469840927028069225873654519021738405746193482980595320530938533835208879166923498550"), // g'(1.359375) + lexical_cast("0.035816505604196436523119991586426561468497089504860912906314914221922063255629240973776624566780655171724"), // g'(1.390625) + lexical_cast("0.029808628117310126968601751537697736254089545807637046082708655932803041328577131366526634532564703021917"), // g'(1.421875) + lexical_cast("0.024661087314753282510573737140694345104475111106211794879844709982238091486098551757845149347618610032885"), // g'(1.453125) + lexical_cast("0.020277183817500123925718353350832363606274271379990517119797309236732993394075485991675565378835117399990"), // g'(1.484375) + lexical_cast("0.016566786254247575375293749008565116603738813921237655148202097271118978259277010195933424056662734396069"), // g'(1.515625) + lexical_cast("0.013446536605285730674148393318910798784847889463939012778503099339814210725517618275874439253140906795321"), // g'(1.546875) + lexical_cast("0.010839937168255907210869762269731543536117453194208729751318975711261680570139811325908435421637044586857"), // g'(1.578125) + lexical_cast("0.008677330749539181585355913836929722539275966612469987097585830633746799618215565932575864813323592569064"), // g'(1.609375) + lexical_cast("0.006895785969066003532912041569588146385648482133964273925008846005414392356354762540479686458456915182422"), // g'(1.640625) + lexical_cast("0.005438899797623998433124926177868203505129237029268854683709340550916273934963934680338382849644109991203"), // g'(1.671875) + lexical_cast("0.004256529599017858016492207077537233128063198584173996714275176466493951517013230345961773807279245818242"), // g'(1.703125) + lexical_cast("0.003304466994034830236341392229394134530658656092624668672443607085027530811446824346183323643773146431650"), // g'(1.734375) + lexical_cast("0.002544065767529172967773767418371787185744918421169971121717939040485439469176241721346003713126884990789"), // g'(1.765625) + lexical_cast("0.001941835775984367581433723487743118747896161816112075901249984974807012598751430317814488690531071486413"), // g'(1.796875) + lexical_cast("0.001469014359942979105844078986521256862651478804967239809867901326844812545018570707402922167337070318165"), // g'(1.828125) + lexical_cast("0.001101126113451938386192529068663844073379771262193065469815230950219555247693085164131775351486650875579"), // g'(1.859375) + lexical_cast("0.000817541013324694931147962586297788606685750454786825608886470281060088745573971194437715594337819829362"), // g'(1.890625) + lexical_cast("0.000601039879911474225734015276657317979859035267662209367133195654905210933494993937372767466898730470954"), // g'(1.921875) + lexical_cast("0.000437394956159116877860869411609646253161082903561375724962103397200094976253287110042483591754739259377"), // g'(1.953125) + lexical_cast("0.000314972091860212002738163289809563177930090611263425563998770969626606901739196068292862372824103237676"), // g'(1.984375) + lexical_cast("0.000224359652050085491041071474105375858702315883318488669316880414932511699763514637604056020585206729702"), // g'(2.015625) + lexical_cast("0.000158027884007011919492306175750514536346161529152522385060799084390555018606145439756610554625333671622"), // g'(2.046875) + lexical_cast("0.000110021128466666972244550229942641180564730883151356076430662293982979610721079431010333297684046009961"), // g'(2.078125) + lexical_cast("0.000075683996586201477788165374446595019081401614817640599970019960133137090438903927564511393112632437065"), // g'(2.109375) + lexical_cast("0.000051421497447658802091816641948792100558148902462530083957875539317483754154795205671660781697877342351"), // g'(2.140625) + lexical_cast("0.000034492124759343197699875355951697823383607724526995377193557806925119821757676060929080600688708285749"), // g'(2.171875) + lexical_cast("0.000022832118109036146591351441225788428615349407345705869110374116111811995689284875925667954577184155616"), // g'(2.203125) + lexical_cast("0.000014908514031870608449073273504082824708512965090788681136021554541816186977506405534445749983533425679"), // g'(2.234375) + lexical_cast("0.000009598194128378471077646466502228564631552196031060212700520428669493372565393532967025917421159451055"), // g'(2.265625) + lexical_cast("0.000006089910032094903925589071885768365033552154835852742380372200679140144575813694932014500333245992150"), // g'(2.296875) + lexical_cast("0.000003806198326464489904501471179191111027481090154526936313549675432967861039485776341397959769354546390"), // g'(2.328125) + lexical_cast("0.000002342166720852809684271714796661027895635740590680844209665411750528687006930156784955258417831088291"), // g'(2.359375) + lexical_cast("0.000001418306715549391752314987533629911529497333350741978774845090291216724311125257284507102899422956865"), // g'(2.390625) + lexical_cast("0.000000844737563848598634692437745232670534830111918817650772401536131863058168276807373058727322485757840"), // g'(2.421875) + lexical_cast("0.000000494582887027541985082961041588843194493377682057522062572717738874999396226021708036282926453353941"), // g'(2.453125) + lexical_cast("0.000000284499236591598063392547619068660729818049847005360315769911532204922664017734687362775492164759895"), // g'(2.484375) + lexical_cast("0.000000160693945790762249108543722438150594237346027537882772733624318355157175540333342507476026614814130"), // g'(2.515625) + lexical_cast("0.000000089071395140242387123752443718466410476313538889832255161675156512656219521550834665275976938868356"), // g'(2.546875) + lexical_cast("0.000000048420950198072369668651027805560415741991640925058727737520004599224608796883917135571227944705408"), // g'(2.578125) + lexical_cast("0.000000025799568229535892380414249079124994790308394653761259451441767353646230715445642196863431273168636"), // g'(2.609375) + lexical_cast("0.000000013464645522302038795727293982948898016459367623322889579168050607638027227263299371956638777679387"), // g'(2.640625) + lexical_cast("0.000000006878461095589900111136392082034183014083825816558701395279223167039349900136169062290105069997768"), // g'(2.671875) + lexical_cast("0.000000003437185674465009051135961298221448530081865035385056331928932641931922771730843393793643760327334"), // g'(2.703125) + lexical_cast("0.000000001678889768216190680690345047883480132698093706961415511874822067072901599214236274008898918308398"), // g'(2.734375) + lexical_cast("0.000000000800997844797296653557042389274831630361435553706219916861697657604610128769179525418316481404345"), // g'(2.765625) + lexical_cast("0.000000000372995018430527900380251336584019890026152442892339581968995175499821780284935610164053236169322"), // g'(2.796875) + lexical_cast("0.000000000169394577894116468756512995325939417364175960421778711053720888629092040431505248908626943368199"), // g'(2.828125) + lexical_cast("0.000000000074967397573818224522463983518852212900281332337849881083903113441220398338524319131856022840614"), // g'(2.859375) + lexical_cast("0.000000000032304464333252365759775564205312505307622216733879458398560977659687434066187535677925097333693"), // g'(2.890625) + lexical_cast("0.000000000013542512912336274431500361753892705431312441925822439518907856136082193371315555649521173566303"), // g'(2.921875) + lexical_cast("0.000000000005518236946817488582064057310269826608690191333604484745956639099226874591346995610017895744578"), // g'(2.953125) + lexical_cast("0.000000000002183592209923360905209996960682785475770239830035281640763341169350903702318435418313485484506"), // g'(2.984375) + lexical_cast("0.000000000000838312896050266709354740899668484911221353018140659322756581004639495581497849888708173298815"), // g'(3.015625) + lexical_cast("0.000000000000311949772868480812347780516506912409007323288594064564370343874294076599123929297097943701034"), // g'(3.046875) + lexical_cast("0.000000000000112402089599228614755295762794129173217506414466830176829682374326571276860325607154373658184"), // g'(3.078125) + lexical_cast("0.000000000000039176794506016467939451671035436525222416381468156163273536964030667848655076274979841889932"), // g'(3.109375) + lexical_cast("0.000000000000013194342231967989407842390664342775745117191635312195908676033481558696665209428463486208527"), // g'(3.140625) + lexical_cast("0.000000000000004289196222067908001815994863079952340273897838658392450299677661912128035665377808982343013"), // g'(3.171875) + lexical_cast("0.000000000000001344322287539522146221080617673048211975604193697619423606843740606538997605021696052730313"), // g'(3.203125) + lexical_cast("0.000000000000000405755770226285762137849766040194966492490426242400304963816867248696901614193771132632268"), // g'(3.234375) + lexical_cast("0.000000000000000117798121272483481430795429265227823389875207214758138237599792271851893488029564688863643"), // g'(3.265625) + lexical_cast("0.000000000000000032853861628847006575461912014493979358530709067619584531982236429803843165396338610640044"), // g'(3.296875) + lexical_cast("0.000000000000000008791316558919890092526902998792709028047982544408870690081761661321472588740075392508761"), // g'(3.328125) + lexical_cast("0.000000000000000002254074830436881944615810338867844822175525846534259182547196778944865636295869684455645"), // g'(3.359375) + lexical_cast("0.000000000000000000553017691284033758949053178662806068054806326469125516762919577442779068366903379217460"), // g'(3.390625) + lexical_cast("0.000000000000000000129645271406893694281896506281765229596780479816278314054216272866619778633667630574782"), // g'(3.421875) + lexical_cast("0.000000000000000000028999645564315719864785175397007915761860516014741538325797768478450821415402826624430"), // g'(3.453125) + lexical_cast("0.000000000000000000006180143249339988447269052329417001570073091758594717163491894057763074807165678977332"), // g'(3.484375) + lexical_cast("0.000000000000000000001252867643227321043362563007101318081461568219865046511052023371704463428404680889318"), // g'(3.515625) + lexical_cast("0.000000000000000000000241225054683610100547624673480448147701505727875355844542994650779545976378140495633"), // g'(3.546875) + lexical_cast("0.000000000000000000000044039066999398680900169681674054684798630835786940549653827399216605257470202939021"), // g'(3.578125) + lexical_cast("0.000000000000000000000007610577807582062575573263413677229331772124424281336240218035661098314864102221811"), // g'(3.609375) + lexical_cast("0.000000000000000000000001242805165212316494523676552887136305778441127251023744635059768345553662775558813"), // g'(3.640625) + lexical_cast("0.000000000000000000000000191431069022397606764880347277712106636642704343748749955269034686788637609751431"), // g'(3.671875) + lexical_cast("0.000000000000000000000000027761251025850583169467065718383927058685109711523535827817309393218519985998495"), // g'(3.703125) + lexical_cast("0.000000000000000000000000003783124072813779705228356271438163581045319609772458825344530460745782474772479"), // g'(3.734375) + lexical_cast("0.000000000000000000000000000483491015481884766524481138636610326553791961644557566572126932357916549658622"), // g'(3.765625) + lexical_cast("0.000000000000000000000000000057831786972290529603663947121744896347612469302867477188325934498277736898412"), // g'(3.796875) + lexical_cast("0.000000000000000000000000000006460575703441721989631276323582594708943380055313130523026515535423593747853"), // g'(3.828125) + lexical_cast("0.000000000000000000000000000000672603738958794053564859771086340933383167706149735466581641424943518853899"), // g'(3.859375) + lexical_cast("0.000000000000000000000000000000065111534511374516546657070065569591183533449520255894748314043378512284243"), // g'(3.890625) + lexical_cast("0.000000000000000000000000000000005847409074548101988461965425962519818528933382988527264024098242130376043"), // g'(3.921875) + lexical_cast("0.000000000000000000000000000000000486004605514227333525776098099024879586315726517256290832491329641099103"), // g'(3.953125) + lexical_cast("0.000000000000000000000000000000000037292395298436024198774848925631641863619212674246073207802818468360567"), // g'(3.984375) + lexical_cast("0.000000000000000000000000000000000002635123061680366694961156623031568922293337299526836931288645171639104"), // g'(4.015625) + lexical_cast("0.000000000000000000000000000000000000171019264014906972413652933319853622109200887319960194386873195193761"), // g'(4.046875) + lexical_cast("0.000000000000000000000000000000000000010166685309552127228580017693019009338225337573956469339645525466356"), // g'(4.078125) + lexical_cast("0.000000000000000000000000000000000000000552069094434276256759290252786473269399999179845742278717819055721"), // g'(4.109375) + lexical_cast("0.000000000000000000000000000000000000000027304755116605008208950089540576931808633121962393012091151683420"), // g'(4.140625) + lexical_cast("0.000000000000000000000000000000000000000001226380966652512847885874020657184963132052922860038712169812377"), // g'(4.171875) + lexical_cast("0.000000000000000000000000000000000000000000049868392983536125725869145787511356632015632683793328037522360"), // g'(4.203125) + lexical_cast("0.000000000000000000000000000000000000000000001830065415771163301567568873064584609169884696317729419132589"), // g'(4.234375) + lexical_cast("0.000000000000000000000000000000000000000000000060413503225082640435207520715838787982618017939346839877514"), // g'(4.265625) + lexical_cast("0.000000000000000000000000000000000000000000000001788000308525712677504883767028310272269396132176347017318"), // g'(4.296875) + lexical_cast("0.000000000000000000000000000000000000000000000000047278206635290633557485891499143010933492804256323487467"), // g'(4.328125) + lexical_cast("0.000000000000000000000000000000000000000000000000001112910171804977195020337584091692655306606526691894796"), // g'(4.359375) + lexical_cast("0.000000000000000000000000000000000000000000000000000023236007287936054634490409808175661127701919340803619"), // g'(4.390625) + lexical_cast("0.000000000000000000000000000000000000000000000000000000428657921354384008188648326796452926850537472514601"), // g'(4.421875) + lexical_cast("0.000000000000000000000000000000000000000000000000000000006959873522553751553524717208478990181867788537220"), // g'(4.453125) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000099053998128674644167844413280198940468591365341"), // g'(4.484375) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000001230569032237017264765448373626102406066245667"), // g'(4.515625) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000013287055463682438112745418077294775519617148"), // g'(4.546875) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000124138446090017620919198186310425312846149"), // g'(4.578125) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000998948898309891630051620362306290253591"), // g'(4.609375) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000006890979289020980623954223971870828237"), // g'(4.640625) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000040550412837679391530484571616360602"), // g'(4.671875) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000202532547805183379185534866559734"), // g'(4.703125) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000000854119386081712900563561139533"), // g'(4.734375) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000000003025058495847807518839405848"), // g'(4.765625) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000000000008948157498373784728600920"), // g'(4.796875) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000000000000021980365957192520850760"), // g'(4.828125) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000044573390842518213190"), // g'(4.859375) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000074167316712100725"), // g'(4.890625) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000100627894464704"), // g'(4.921875) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000110606170249"), // g'(4.953125) + lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000097834582"), // g'(4.984375) + }, + }; +}; + +template +tanh_sinh_detail::tanh_sinh_detail(Real tol, size_t max_refinements) +{ + m_tol = tol; + m_max_refinements = max_refinements; + /* + * Our goal is to calculate t_max such that tanh(pi/2 sinh(t_max)) < 1 in the requested precision. + * What follows is a good estimate for t_max, but in fact we can get closer by starting with this estimate + * and then calculating tanh(pi/2 sinh(t_max + eps)) until it = 1 (the second to last is t_max). + * However, this is computationally expensive, so we can't do it. + * An alternative is to cache the value of t_max for various types (float, double, long double, float128, cpp_bin_float_50, cpp_bin_float_100) + * and then simply use them, but unfortunately the values will be platform dependent. + * As such we are then susceptible to the catastrophe where we evaluate the function at x = 1, when we have promised we wouldn't do that. + * Other authors solve this problem by computing the abscissas in double the requested precision, and then returning the result at the request precision. + * This seems to be overkill to me, but presumably it's an option if we find integrals on which this method struggles. + */ + + using std::tanh; + using std::sinh; + using std::asinh; + using std::atanh; + using boost::math::constants::half_pi; + using boost::math::constants::two_div_pi; + + auto g = [](Real t) { return tanh(half_pi()*sinh(t)); }; + auto g_inv = [](Real x) { return asinh(two_div_pi()*atanh(x)); }; + + Real x = float_prior((Real) 1); + m_t_max = g_inv(x); + while(!isnormal(m_t_max)) + { + // Although you might suspect that we could be in this loop essentially for ever, in point of fact it is only called once + // even for 100 digit accuracy, and isn't called at all up to float128. + x = float_prior(x); + m_t_max = g_inv(x); + } + // This occurs once on 100 digit arithmetic: + while(!(g(m_t_max) < (Real) 1)) + { + x = float_prior(x); + m_t_max = g_inv(x); + } +} + + +template +template +Real tanh_sinh_detail::integrate(const F f, Real* error, Real* L1) const +{ + using std::abs; + using std::floor; + using std::tanh; + using std::sinh; + using std::sqrt; + using boost::math::constants::half; + using boost::math::constants::half_pi; + Real h = 1; + Real I0 = half_pi()*f(0); + Real L1_I0 = abs(I0); + for(size_t i = 1; i <= m_t_max; ++i) + { + Real x = m_abscissas[0][i]; + Real w = m_weights[0][i]; + Real yp = f(x); + Real ym = f(-x); + I0 += (yp + ym)*w; + L1_I0 += (abs(yp) + abs(ym))*w; + } + // Uncomment to understand the convergence rate: + // std::cout << std::setprecision(std::numeric_limits::digits10); + // std::cout << "First estimate : " << I0 << std::endl; + Real I1 = half()*I0; + Real L1_I1 = half()*L1_I0; + h /= 2; + Real sum = 0; + Real absum = 0; + for(size_t i = 0; i < m_weights[1].size() && h + i <= m_t_max; ++i) + { + Real x = m_abscissas[1][i]; + Real w = m_weights[1][i]; + + Real yp = f(x); + Real ym = f(-x); + sum += (yp + ym)*w; + absum += (abs(yp) + abs(ym))*w; + } + I1 += sum*h; + L1_I1 += absum*h; + Real err = abs(I0 - I1); + + // std::cout << "Second estimate: " << I1 << " Error estimate at level " << 1 << " = " << err << std::endl; + + size_t k = 2; + + while (k < 4 || (k < m_weights.size() && k < m_max_refinements) ) + { + I0 = I1; + L1_I0 = L1_I1; + + I1 = half()*I0; + L1_I1 = half()*L1_I0; + h *= half(); + Real sum = 0; + Real absum = 0; + auto const& abscissa_row = m_abscissas[k]; + auto const& weight_row = m_weights[k]; + + // We have no guarantee that round-off error won't cause the function to be evaluated strictly at the endpoints. + // However, making sure x < 1 - eps is a reasonable compromise between accuracy and usability.. + for(size_t j = 0; j < weight_row.size() && abscissa_row[j] < (Real) 1 - std::numeric_limits::epsilon(); ++j) + { + Real x = abscissa_row[j]; + Real w = weight_row[j]; + + Real yp = f(x); + Real ym = f(-x); + Real term = (yp + ym)*w; + sum += term; + + // A question arises as to how accurately we actually need to estimate the L1 integral. + // For simple integrands, computing the L1 norm makes the integration 20% slower, + // but for more complicated integrands, this calculation is not noticeable. + Real abterm = (abs(yp) + abs(ym))*w; + absum += abterm; + } + + I1 += sum*h; + L1_I1 += absum*h; + ++k; + err = abs(I0 - I1); + // std::cout << "Estimate: " << I1 << " Error estimate at level " << k << " = " << err << std::endl; + + if (!isfinite(I1)) + { + throw std::domain_error("The tanh_sinh quadrature evaluated your function at a singular point. Please narrow the bounds of integration or check your function for singularities.\n"); + } + + if (err <= m_tol*L1_I1) + { + break; + } + + } + + // Since we are out of precomputed abscissas and weights, switch to computing abscissas and weights on the fly. + while (k < m_max_refinements && err > m_tol*L1_I1) + { + I0 = I1; + L1_I0 = L1_I1; + + I1 = half()*I0; + L1_I1 = half()*L1_I0; + h *= half(); + Real sum = 0; + Real absum = 0; + for(Real t = h; t < m_t_max - std::numeric_limits::epsilon(); t += 2*h) + { + Real s = sinh(t); + Real c = sqrt(1+s*s); + Real x = tanh(half_pi()*s); + Real w = half_pi()*c*(1-x*x); + + Real yp = f(x); + Real ym = f(-x); + Real term = (yp + ym)*w; + sum += term; + Real abterm = (abs(yp) + abs(ym))*w; + absum += abterm; + // There are claims that this test improves performance, + // however my benchmarks show that it's slower! + // However, I leave this comment here because it totally stands to reason that this *should* help: + //if (abterm < std::numeric_limits::epsilon()) { break; } + } + + I1 += sum*h; + L1_I1 += absum*h; + ++k; + err = abs(I0 - I1); + + // std::cout << "Estimate: " << I1 << " Error estimate at level " << k << " = " << err << std::endl; + if (!isfinite(I1)) + { + throw std::domain_error("The tanh_sinh quadrature evaluated your function at a singular point. Please narrow the bounds of integration or check your function for singularities.\n"); + } + + if (err <= m_tol*L1_I1) + { + break; + } + } + + if (error) + { + *error = err; + } + + if (L1) + { + *L1 = L1_I1; + } + + return I1; +} + +}}} +#endif diff --git a/include/boost/math/quadrature/exp_sinh.hpp b/include/boost/math/quadrature/exp_sinh.hpp new file mode 100644 index 0000000000..0dc023c5ae --- /dev/null +++ b/include/boost/math/quadrature/exp_sinh.hpp @@ -0,0 +1,88 @@ +// Copyright Nick Thompson, 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) + +/* + * This class performs exp-sinh quadrature on half infinite intervals. + * + * References: + * + * 1) Tanaka, Ken’ichiro, et al. "Function classes for double exponential integration formulas." Numerische Mathematik 111.4 (2009): 631-655. + */ + +#ifndef BOOST_MATH_QUADRATURE_EXP_SINH_HPP +#define BOOST_MATH_QUADRATURE_EXP_SINH_HPP + +#include +#include +#include +#include + +namespace boost{ namespace math{ + + +template +class exp_sinh +{ +public: + exp_sinh(Real tol = sqrt(std::numeric_limits::epsilon()), size_t max_refinements = 9); + + template + Real integrate(const F f, Real a = 0, Real b = std::numeric_limits::infinity(), Real* error = nullptr, Real* L1 = nullptr) const; + +private: + std::shared_ptr> m_imp; +}; + +template +exp_sinh::exp_sinh(Real tol, size_t max_refinements) : m_imp(std::make_shared>(tol, max_refinements)) +{ + return; +} + +template +template +Real exp_sinh::integrate(const F f, Real a, Real b, Real* error, Real* L1) const +{ + using std::isfinite; + using std::abs; + using boost::math::constants::half; + using boost::math::detail::exp_sinh_detail; + + // Right limit is infinite: + if (isfinite(a) && b >= std::numeric_limits::max()) + { + // If a = 0, don't use an additional level of indirection: + if (a == (Real) 0) + { + return m_imp->integrate(f, error, L1); + } + const auto u = [&](Real t) { return f(t + a); }; + return m_imp->integrate(u, error, L1); + } + + if (isfinite(b) && a <= std::numeric_limits::lowest()) + { + const auto u = [&](Real t) { return f(b-t);}; + return m_imp->integrate(u, error, L1); + } + + if (isfinite(a) && isfinite(b)) + { + throw std::domain_error("Use tanh_sinh quadrature for integration over finite domains; exp_sinh is for half infinite integrals.\n"); + } + + // Infinite limits: + if (a <= std::numeric_limits::lowest() && b >= std::numeric_limits::max()) + { + throw std::domain_error("Use sinh_sinh quadrature for integration over the whole real line; exp_sinh is for half infinite integrals.\n"); + } + + throw std::logic_error("The domain of integration is not sensible; please check the bounds.\n"); +} + + +}} +#endif diff --git a/include/boost/math/quadrature/sinh_sinh.hpp b/include/boost/math/quadrature/sinh_sinh.hpp new file mode 100644 index 0000000000..18516e2eb9 --- /dev/null +++ b/include/boost/math/quadrature/sinh_sinh.hpp @@ -0,0 +1,53 @@ +// Copyright Nick Thompson, 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) + +/* + * This class performs sinh-sinh quadrature over the entire real line. + * + * References: + * + * 1) Tanaka, Ken’ichiro, et al. "Function classes for double exponential integration formulas." Numerische Mathematik 111.4 (2009): 631-655. + */ + +#ifndef BOOST_MATH_QUADRATURE_SINH_SINH_HPP +#define BOOST_MATH_QUADRATURE_SINH_SINH_HPP + +#include +#include +#include +#include + +namespace boost{ namespace math{ + + +template +class sinh_sinh +{ +public: + sinh_sinh(Real tol = sqrt(std::numeric_limits::epsilon()), size_t max_refinements = 9); + + template + Real integrate(const F f, Real* error = nullptr, Real* L1 = nullptr) const; + +private: + std::shared_ptr> m_imp; +}; + +template +sinh_sinh::sinh_sinh(Real tol, size_t max_refinements) : m_imp(std::make_shared>(tol, max_refinements)) +{ + return; +} + +template +template +Real sinh_sinh::integrate(const F f, Real* error, Real* L1) const +{ + return m_imp->integrate(f, error, L1); +} + +}} +#endif diff --git a/include/boost/math/quadrature/tanh_sinh.hpp b/include/boost/math/quadrature/tanh_sinh.hpp new file mode 100644 index 0000000000..9ca1e886e6 --- /dev/null +++ b/include/boost/math/quadrature/tanh_sinh.hpp @@ -0,0 +1,123 @@ +// Copyright Nick Thompson, 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) + +/* + * This class performs tanh-sinh quadrature on the real line. + * Tanh-sinh quadrature is exponentially convergent for integrands in Hardy spaces, + * (see https://en.wikipedia.org/wiki/Hardy_space for a formal definition), and is optimal for a random function from that class. + * + * The tanh-sinh quadrature is one of a class of so called "double exponential quadratures"-there is a large family of them, + * but this one seems to be the most commonly used. + * + * As always, there are caveats: For instance, if the function you want to integrate is not holomorphic on the unit disk, + * then the rapid convergence will be spoiled. In this case, a more appropriate quadrature is (say) Romberg, which does not + * require the function to be holomorphic, only differentiable up to some order. + * + * In addition, if you are integrating a periodic function over a period, the trapezoidal rule is better. + * + * References: + * + * 1) Mori, Masatake. "Quadrature formulas obtained by variable transformation and the DE-rule." Journal of Computational and Applied Mathematics 12 (1985): 119-130. + * 2) Bailey, David H., Karthik Jeyabalan, and Xiaoye S. Li. "A comparison of three high-precision quadrature schemes." Experimental Mathematics 14.3 (2005): 317-329. + * 3) Press, William H., et al. "Numerical recipes third edition: the art of scientific computing." Cambridge University Press 32 (2007): 10013-2473. + * + */ + +#ifndef BOOST_MATH_QUADRATURE_TANH_SINH_HPP +#define BOOST_MATH_QUADRATURE_TANH_SINH_HPP + +#include +#include +#include +#include + +namespace boost{ namespace math{ + + +template +class tanh_sinh +{ +public: + tanh_sinh(Real tol = sqrt(std::numeric_limits::epsilon()), size_t max_refinements = 15); + + template + Real integrate(const F f, Real a, Real b, Real* error = nullptr, Real* L1 = nullptr) const; + +private: + std::shared_ptr> m_imp; +}; + +template +tanh_sinh::tanh_sinh(Real tol, size_t max_refinements) : m_imp(std::make_shared>(tol, max_refinements)) +{ + return; +} + +template +template +Real tanh_sinh::integrate(const F f, Real a, Real b, Real* error, Real* L1) const +{ + using std::isfinite; + using boost::math::constants::half; + using boost::math::detail::tanh_sinh_detail; + + if (isfinite(a) && isfinite(b)) + { + if (b <= a) + { + throw std::domain_error("Arguments to integrate are in wrong order; integration over [a,b] must have b > a.\n"); + } + Real avg = (a+b)*half(); + Real diff = (b-a)*half(); + auto u = [&](Real z) { return f(avg + diff*z); }; + Real Q = diff*m_imp->integrate(u, error, L1); + + if(L1) + { + *L1 *= diff; + } + return Q; + } + + // Infinite limits: + if (a <= std::numeric_limits::lowest() && b >= std::numeric_limits::max()) + { + auto u = [&](Real t) { auto t_sq = t*t; auto inv = 1/(1 - t_sq); return f(t*inv)*(1+t_sq)*inv*inv; }; + return m_imp->integrate(u, error, L1); + } + + // Right limit is infinite: + if (isfinite(a) && b >= std::numeric_limits::max()) + { + auto u = [&](Real t) { auto z = 1/(t+1); auto arg = 2*z + a - 1; return f(arg)*z*z; }; + Real Q = 2*m_imp->integrate(u, error, L1); + if(L1) + { + *L1 *= 2; + } + + return Q; + } + + if (isfinite(b) && a <= std::numeric_limits::lowest()) + { + auto u = [&](Real t) { return f(b-t);}; + auto v = [&](Real t) { auto z = 1/(t+1); auto arg = 2*z - 1; return u(arg)*z*z; }; + + Real Q = 2*m_imp->integrate(v, error, L1); + if (L1) + { + *L1 *= 2; + } + return Q; + } + + throw std::logic_error("The domain of integration is not sensible; please check the bounds.\n"); +} + + +}} +#endif diff --git a/test/exp_sinh_quadrature_test.cpp b/test/exp_sinh_quadrature_test.cpp new file mode 100644 index 0000000000..cf9e4d6893 --- /dev/null +++ b/test/exp_sinh_quadrature_test.cpp @@ -0,0 +1,327 @@ +#define BOOST_TEST_MODULE exp_sinh_quadrature_test + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#ifdef __GNUC__ +#ifndef __clang__ +#include +using boost::multiprecision::float128; +#endif +#endif + +using std::expm1; +using std::exp; +using std::cos; +using std::atan; +using std::tan; +using std::log; +using std::log1p; +using std::asinh; +using std::atanh; +using std::sqrt; +using std::isnormal; +using std::abs; +using std::sinh; +using std::tanh; +using std::cosh; +using std::pow; +using std::string; +using boost::multiprecision::cpp_bin_float_50; +using boost::multiprecision::cpp_bin_float_100; +using boost::multiprecision::cpp_dec_float_50; +using boost::multiprecision::cpp_dec_float_100; +using boost::multiprecision::cpp_bin_float; +using boost::math::exp_sinh; +using boost::math::constants::pi; +using boost::math::constants::half_pi; +using boost::math::constants::two_div_pi; +using boost::math::constants::half; +using boost::math::constants::third; +using boost::math::constants::half; +using boost::math::constants::third; +using boost::math::constants::catalan; +using boost::math::constants::ln_two; +using boost::math::constants::root_two; +using boost::math::constants::root_two_pi; +using boost::math::constants::root_pi; + + +template +void test_right_limit_infinite() +{ + std::cout << "Testing right limit infinite for tanh_sinh in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = sqrt(std::numeric_limits::epsilon()); + Real Q; + Real Q_expected; + Real error; + Real L1; + exp_sinh integrator(tol, 12); + + // Example 12 + const auto f2 = [](Real t) { return exp(-t)/sqrt(t); }; + Q = integrator.integrate(f2, 0, std::numeric_limits::infinity(), &error, &L1); + Q_expected = root_pi(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + // The integrand is strictly positive, so it coincides with the value of the integral: + BOOST_CHECK_CLOSE(L1, Q_expected, 100*tol); + + auto f3 = [](Real t) { return exp(-t)*cos(t); }; + Q = integrator.integrate(f3, 0, std::numeric_limits::infinity(), &error, &L1); + Q_expected = half(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + auto f4 = [](Real t) { return 1/(1+t*t); }; + Q = integrator.integrate(f4, 1, std::numeric_limits::infinity(), &error, &L1); + Q_expected = pi()/4; + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + BOOST_CHECK_CLOSE(L1, Q_expected, 100*tol); +} + +template +void test_left_limit_infinite() +{ + std::cout << "Testing left limit infinite for 1/(1+t^2) on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = sqrt(std::numeric_limits::epsilon()); + Real Q; + Real Q_expected; + Real error; + Real L1; + exp_sinh integrator(tol, 8); + + // Example 11: + auto f1 = [](Real t) { return 1/(1+t*t);}; + Q = integrator.integrate(f1, -std::numeric_limits::infinity(), 0, &error, &L1); + Q_expected = half_pi(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + BOOST_CHECK_CLOSE(L1, Q_expected, 100*tol); +} + + +// Some examples of tough integrals from NR, section 4.5.4: +template +void test_nr_examples() +{ + using std::sin; + using std::pow; + using std::exp; + using std::sqrt; + std::cout << "Testing type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = sqrt(std::numeric_limits::epsilon()); + std::cout << std::setprecision(std::numeric_limits::digits10); + Real Q; + Real Q_expected; + Real L1; + Real error; + exp_sinh integrator(tol, 12); + + auto f0 = [](Real x) { return (Real) 0; }; + Q = integrator.integrate(f0, 0, std::numeric_limits::infinity(), &error, &L1); + Q_expected = 0; + BOOST_CHECK_CLOSE(Q, 0, 100*tol); + BOOST_CHECK_CLOSE(L1, 0, 100*tol); + + auto f = [](Real x) { return 1/(1+x*x); }; + Q = integrator.integrate(f, 0, std::numeric_limits::infinity(), &error, &L1); + Q_expected = half_pi(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + BOOST_CHECK_CLOSE(L1, Q_expected, 100*tol); + + auto f1 = [](Real x) { + Real z1 = exp(-x); + if (z1 == 0) + { + return (Real) 0; + } + Real z2 = pow(x, -3*half())*z1; + if (z2 == 0) + { + return (Real) 0; + } + return sin(x*half())*z2; + }; + + Q = integrator.integrate(f1, 0, std::numeric_limits::infinity(), &error, &L1); + Q_expected = sqrt(pi()*(sqrt((Real) 5) - 2)); + + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + auto f2 = [](Real x) { return pow(x, -(Real) 2/(Real) 7)*exp(-x*x); }; + Q = integrator.integrate(f2, 0, std::numeric_limits::infinity(), &error, &L1); + Q_expected = half()*boost::math::tgamma((Real) 5/ (Real) 14); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + BOOST_CHECK_CLOSE(L1, Q_expected, 100*tol); + + auto f3 = [](Real x) { return (Real) 1/ (sqrt(x)*(1+x)); }; + Q = integrator.integrate(f3, 0, std::numeric_limits::infinity(), &error, &L1); + Q_expected = pi(); + + BOOST_CHECK_CLOSE(Q, Q_expected, 1000*std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(L1, Q_expected, 1000*std::numeric_limits::epsilon()); + + auto f4 = [](Real t) { return exp(-t*t*half()); }; + Q = integrator.integrate(f4, 0, std::numeric_limits::infinity(), &error, &L1); + Q_expected = root_two_pi()/2; + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + BOOST_CHECK_CLOSE(L1, Q_expected, 100*tol); + + // This test shows how oscillatory integrals with 1/t decay are approximated very poorly by this method. + // In fact |f(x)| <= C/(1+x^(1+eps)) for large x for this method to converge. + //Q = integrator.integrate(boost::math::sinc_pi, 0, std::numeric_limits::infinity()); + //Q_expected = half_pi(); + //BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + auto f5 = [](Real t) { return 1/cosh(t);}; + Q = integrator.integrate(f5, 0, std::numeric_limits::infinity(), &error, &L1); + Q_expected = half_pi(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + BOOST_CHECK_CLOSE(L1, Q_expected, 100*tol); +} + +// Definite integrals found in the CRC Handbook of Mathematical Formulas +template +void test_crc() +{ + using std::sin; + using std::pow; + using std::exp; + using std::sqrt; + using std::log; + std::cout << "Testing integral from CRC handbook on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = sqrt(std::numeric_limits::epsilon()); + std::cout << std::setprecision(std::numeric_limits::digits10); + Real Q; + Real Q_expected; + Real L1; + Real error; + exp_sinh integrator(tol, 14); + + auto f0 = [](Real x) { return log(x)*exp(-x); }; + Q = integrator.integrate(f0, 0, std::numeric_limits::infinity(), &error, &L1); + Q_expected = -boost::math::constants::euler(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + + // Test the integral representation of the gamma function: + auto f1 = [](Real t) { Real x = exp(-t); + if(x == 0) + { + return (Real) 0; + } + return pow(t, (Real) 12 - 1)*x; + }; + + Q = integrator.integrate(f1, 0, std::numeric_limits::infinity(), &error, &L1); + Q_expected = boost::math::tgamma(12); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + + // Integral representation of the modified bessel function: + // K_5(12) + auto f2 = [](Real t) { + Real x = exp(-12*cosh(t)); + if (x == 0) + { + return (Real) 0; + } + return x*cosh(5*t); + }; + Q = integrator.integrate(f2, 0, std::numeric_limits::infinity(), &error, &L1); + Q_expected = boost::math::cyl_bessel_k(5, 12); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + // Laplace transform of cos(at) + Real a = 20; + Real s = 1; + auto f3 = [&](Real t) { + Real x = exp(-s*t); + if (x == 0) + { + return (Real) 0; + } + return cos(a*t)*x; + }; + + // For high oscillation frequency, the quadrature sum is ill-conditioned. + Q = integrator.integrate(f3, 0, std::numeric_limits::infinity(), &error, &L1); + Q_expected = s/(a*a+s*s); + // Since the integrand is oscillatory, we increase the tolerance: + BOOST_CHECK_CLOSE(Q, Q_expected, 10000*tol); + + // Laplace transform of J_0(t): + auto f4 = [&](Real t) { + Real x = exp(-s*t); + if (x == 0) + { + return (Real) 0; + } + return boost::math::cyl_bessel_j(0, t)*x; + }; + + Q = integrator.integrate(f4, 0, std::numeric_limits::infinity(), &error, &L1); + Q_expected = 1/sqrt(1+s*s); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + // The summation condition number is good, but the integral is oscillatory. + // Convergence is slow, but the final answer is reasonable. + auto f5 = [](Real t) { Real x = boost::math::sinc_pi(t); return x*x; }; + Q = integrator.integrate(f5, 0, std::numeric_limits::infinity(), &error, &L1); + Q_expected = boost::math::constants::half_pi(); + BOOST_CHECK_CLOSE(Q, Q_expected, 0.1); + + auto f6 = [](Real t) { return exp(-t*t)*log(t);}; + Q = integrator.integrate(f6, 0, std::numeric_limits::infinity(), &error, &L1); + Q_expected = -boost::math::constants::root_pi()*(boost::math::constants::euler() + 2*ln_two())/4; + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + +} + +BOOST_AUTO_TEST_CASE(exp_sinh_quadrature_test) +{ + test_left_limit_infinite(); + + test_right_limit_infinite(); + test_right_limit_infinite(); + test_right_limit_infinite(); + #ifdef __GNUC__ + #ifndef __clang__ + test_right_limit_infinite(); + #endif + #endif + + // Unfortunately this fails with message + // boost/multiprecision/detail/functions/trig.hpp(69): fatal error: in "void boost::multiprecision::default_ops::hyp0F1(T&, const T&, const T&) + // [with T = boost::multiprecision::backends::cpp_bin_float<50u>]": boost::exception_detail::clone_impl >: + // H0F1 Failed to Converge + // test_nr_examples(); + test_nr_examples(); + test_nr_examples(); + test_nr_examples(); + #ifdef __GNUC__ + #ifndef __clang__ + test_nr_examples(); + #endif + #endif + + test_crc(); + test_crc(); + test_crc(); + #ifdef __GNUC__ + #ifndef __clang__ + test_crc(); + #endif + #endif + +} diff --git a/test/sinh_sinh_quadrature_test.cpp b/test/sinh_sinh_quadrature_test.cpp new file mode 100644 index 0000000000..e1b97d3e20 --- /dev/null +++ b/test/sinh_sinh_quadrature_test.cpp @@ -0,0 +1,185 @@ +#define BOOST_TEST_MODULE sinh_sinh_quadrature_test + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#ifdef __GNUC__ +#ifndef __clang__ +#include +using boost::multiprecision::float128; +#endif +#endif + +using std::expm1; +using std::exp; +using std::sin; +using std::cos; +using std::atan; +using std::tan; +using std::log; +using std::log1p; +using std::asinh; +using std::atanh; +using std::sqrt; +using std::isnormal; +using std::abs; +using std::sinh; +using std::tanh; +using std::cosh; +using std::pow; +using std::string; +using boost::multiprecision::cpp_bin_float_50; +using boost::multiprecision::cpp_bin_float_100; +using boost::multiprecision::cpp_dec_float_50; +using boost::multiprecision::cpp_dec_float_100; +using boost::multiprecision::cpp_bin_float; +using boost::math::sinh_sinh; +using boost::math::constants::pi; +using boost::math::constants::pi_sqr; +using boost::math::constants::half_pi; +using boost::math::constants::two_div_pi; +using boost::math::constants::half; +using boost::math::constants::third; +using boost::math::constants::half; +using boost::math::constants::third; +using boost::math::constants::catalan; +using boost::math::constants::ln_two; +using boost::math::constants::root_two; +using boost::math::constants::root_two_pi; +using boost::math::constants::root_pi; + + +template +void test_nr_examples() +{ + std::cout << "Testing type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = sqrt(std::numeric_limits::epsilon()); + std::cout << std::setprecision(std::numeric_limits::digits10); + Real Q; + Real Q_expected; + Real L1; + Real error; + sinh_sinh integrator(tol, 10); + + auto f0 = [](Real x) { return (Real) 0; }; + Q = integrator.integrate(f0, &error, &L1); + Q_expected = 0; + BOOST_CHECK_CLOSE(Q, 0, 100*tol); + BOOST_CHECK_CLOSE(L1, 0, 100*tol); + + // In spite of the poles at \pm i, we still get a doubling of the correct digits at each level of refinement. + auto f1 = [](Real t) { return 1/(1+t*t); }; + Q = integrator.integrate(f1, &error, &L1); + Q_expected = pi(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + BOOST_CHECK_CLOSE(L1, Q_expected, 100*tol); + + auto f2 = [](Real x) { return exp(-x*x); }; + Q = integrator.integrate(f2, &error, &L1); + Q_expected = root_pi(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + BOOST_CHECK_CLOSE(L1, Q_expected, 100*tol); + + + // This test shows how oscillatory integrals with 1/t decay are approximated very poorly by this method. + // In fact |f(x)| <= C/(1+x^(1+eps)) for large x for this method to converge. + //Q = integrator.integrate(boost::math::sinc_pi); + //Q_expected = pi(); + //BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + auto f5 = [](Real t) { return 1/cosh(t);}; + Q = integrator.integrate(f5, &error, &L1); + Q_expected = pi(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + BOOST_CHECK_CLOSE(L1, Q_expected, 100*tol); + + // This demonstrates that oscillatory integrals which are nonetheless L1 integrable converge, but so slowly as to be useless. + //std::cout << "cos(t)/(1+t*t)\n"; + //auto f6 = [](Real t) { return cos(t)/(1+t*t);}; + //Q = integrator.integrate(f6, &error, &L1); + //Q_expected = pi()/boost::math::constants::e(); + //BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + //std::cout << "Exact " << Q_expected << std::endl; + //std::cout << "L1 " << L1 << std::endl; + + // This oscillatory integral has rapid convergence because the oscillations get swamped by the exponential growth of the denominator. + auto f8 = [](Real t) { return cos(t)/cosh(t);}; + Q = integrator.integrate(f8, &error, &L1); + Q_expected = pi()/cosh(half_pi()); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + +} + +// Test formulas for in the CRC Handbook of Mathematical functions, 32nd edition. +template +void test_crc() +{ + std::cout << "Testing CRC formulas on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = sqrt(std::numeric_limits::epsilon()); + std::cout << std::setprecision(std::numeric_limits::digits10); + Real Q; + Real Q_expected; + Real L1; + Real error; + sinh_sinh integrator(tol, 10); + + // CRC Definite integral 698: + auto f0 = [](Real x) { + if(x == 0) { + return (Real) 1; + } + return x/sinh(x); + }; + Q = integrator.integrate(f0, &error, &L1); + Q_expected = pi_sqr()/2; + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + BOOST_CHECK_CLOSE(L1, Q_expected, 100*tol); + + + // CRC Definite integral 695: + auto f1 = [](Real x) { + if(x == 0) { + return (Real) 1; + } + return (Real) sin(x)/sinh(x); + }; + Q = integrator.integrate(f1, &error, &L1); + Q_expected = pi()*tanh(half_pi()); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); +} + + +BOOST_AUTO_TEST_CASE(sinh_sinh_quadrature_test) +{ + test_nr_examples(); + test_nr_examples(); + test_nr_examples(); + #ifdef __GNUC__ + #ifndef __clang__ + test_nr_examples(); + #endif + #endif + + test_crc(); + test_crc(); + test_crc(); + #ifdef __GNUC__ + #ifndef __clang__ + test_crc(); + #endif + #endif + + + // Unfortunately this fails with message + // boost/multiprecision/detail/functions/trig.hpp(69): fatal error: in "void boost::multiprecision::default_ops::hyp0F1(T&, const T&, const T&) + // [with T = boost::multiprecision::backends::cpp_bin_float<50u>]": boost::exception_detail::clone_impl >: + // H0F1 Failed to Converge + //test_nr_examples(); +} diff --git a/test/tanh_sinh_quadrature_test.cpp b/test/tanh_sinh_quadrature_test.cpp new file mode 100644 index 0000000000..50440eb77f --- /dev/null +++ b/test/tanh_sinh_quadrature_test.cpp @@ -0,0 +1,598 @@ +#define BOOST_TEST_MODULE tanh_sinh_quadrature_test + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#ifdef __GNUC__ +#ifndef __clang__ +#include +using boost::multiprecision::float128; +#endif +#endif + +using std::expm1; +using std::atan; +using std::tan; +using std::log; +using std::log1p; +using std::asinh; +using std::atanh; +using std::sqrt; +using std::isnormal; +using std::abs; +using std::sinh; +using std::tanh; +using std::cosh; +using std::pow; +using std::string; +using boost::multiprecision::cpp_bin_float_50; +using boost::multiprecision::cpp_bin_float_100; +using boost::multiprecision::cpp_dec_float_50; +using boost::multiprecision::cpp_dec_float_100; +using boost::multiprecision::cpp_bin_float; +using boost::math::tanh_sinh; +using boost::math::detail::tanh_sinh_detail; +using boost::math::constants::pi; +using boost::math::constants::half_pi; +using boost::math::constants::two_div_pi; +using boost::math::constants::two_pi; +using boost::math::constants::half; +using boost::math::constants::third; +using boost::math::constants::half; +using boost::math::constants::third; +using boost::math::constants::catalan; +using boost::math::constants::ln_two; +using boost::math::constants::root_two; +using boost::math::constants::root_two_pi; +using boost::math::constants::root_pi; + +template +Real g(Real t) +{ + return tanh(half_pi()*sinh(t)); +} + +template +Real g_prime(Real t) +{ + + Real denom = cosh(half_pi()*sinh(t)); + return half_pi()*cosh(t)/(denom*denom); +} + +void print_xt(cpp_dec_float_100 x, cpp_dec_float_100 t) +{ + std::cout << std::setprecision(std::numeric_limits::digits10 + 5) << std::fixed; + std::cout << std::fixed << " lexical_cast(\"" << x << "\"), // g(" << std::defaultfloat << t << ")\n"; +} + +void print_xt_prime(cpp_dec_float_100 x, cpp_dec_float_100 t) +{ + std::cout << std::setprecision(std::numeric_limits::digits10 + 5) << std::fixed; + std::cout << std::fixed << " lexical_cast(\"" << x << "\"), // g'(" << std::defaultfloat << t << ")\n"; +} + + +void generate_constants() +{ + + std::cout << " static constexpr const std::vector> abscissas{\n"; + std::cout << " {\n"; + print_xt(g(0), 0); + print_xt(g(1), 1); + print_xt(g(2), 2); + print_xt(g(3), 3); + print_xt(g(4), 4); + print_xt(g(5), 5); + std::cout << " },\n"; + size_t k = 1; + while(k < 9) + { + cpp_dec_float_100 h = (cpp_dec_float_100) 1/ (cpp_dec_float_100) (1 << k); + + std::cout << " {\n"; + for(cpp_dec_float_100 t = h; t <= 5; t += 2*h) + { + auto x = g(t); + print_xt(x, t); + } + std::cout << " },\n"; + ++k; + } + std::cout << "};\n"; + + + std::cout << " static constexpr const std::vector> weights{\n"; + std::cout << " {\n"; + print_xt_prime(g_prime(0), 0); + print_xt_prime(g_prime(1), 1); + print_xt_prime(g_prime(2), 2); + print_xt_prime(g_prime(3), 3); + print_xt_prime(g_prime(4), 4); + print_xt_prime(g_prime(5), 5); + std::cout << " },\n"; + k = 1; + while(k < 9) + { + cpp_dec_float_100 h = (cpp_dec_float_100) 1/ (cpp_dec_float_100) (1 << k); + + std::cout << " {\n"; + for(cpp_dec_float_100 t = h; t <= 5; t += 2*h) + { + auto x = g_prime(t); + print_xt_prime(x, t); + } + std::cout << " },\n"; + ++k; + } + std::cout << " };\n"; + +} + + +template +void test_detail() +{ + std::cout << "Testing tanh_sinh_detail on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = sqrt(std::numeric_limits::epsilon()); + tanh_sinh_detail integrator(tol, 20); + auto f = [](Real x) { return x*x; }; + Real err; + Real L1; + Real Q = integrator.integrate(f, &err, &L1); + BOOST_CHECK_CLOSE(Q, 2*third(), 100*tol); + BOOST_CHECK_CLOSE(L1, 2*third(), 100*tol); +} + + +template +void test_linear() +{ + std::cout << "Testing linear functions are integrated properly by tanh_sinh on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = 100*std::numeric_limits::epsilon(); + tanh_sinh integrator(tol, 20); + auto f = [](Real x) { return 5*x + 7; }; + Real error; + Real L1; + Real Q = integrator.integrate(f, (Real) 0, (Real) 1, &error, &L1); + BOOST_CHECK_CLOSE(Q, 9.5, 100*tol); + BOOST_CHECK_CLOSE(L1, 9.5, 100*tol); +} + + +template +void test_quadratic() +{ + std::cout << "Testing quadratic functions are integrated properly by tanh_sinh on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = 100*std::numeric_limits::epsilon(); + tanh_sinh integrator(tol, 20); + auto f = [](Real x) { return 5*x*x + 7*x + 12; }; + Real error; + Real L1; + Real Q = integrator.integrate(f, (Real) 0, (Real) 1, &error, &L1); + BOOST_CHECK_CLOSE(Q, (Real) 17 + half()*third(), 100*tol); + BOOST_CHECK_CLOSE(L1, (Real) 17 + half()*third(), 100*tol); +} + + +template +void test_singular() +{ + std::cout << "Testing integration of endpoint singularities on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = 100*std::numeric_limits::epsilon(); + Real error; + Real L1; + tanh_sinh integrator(tol, 20); + auto f = [](Real x) { return log(x)*log(1-x); }; + Real Q = integrator.integrate(f, (Real) 0, (Real) 1, &error, &L1); + Real Q_expected = 2 - pi()*pi()*half()*third(); + + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + BOOST_CHECK_CLOSE(L1, Q_expected, 100*tol); +} + + +// Examples taken from +//http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf +template +void test_ca() +{ + std::cout << "Testing integration of C(a) on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = sqrt(std::numeric_limits::epsilon()); + Real error; + Real L1; + + tanh_sinh integrator(tol, 20); + auto f1 = [](Real x) { return atan(x)/(x*(x*x + 1)) ; }; + Real Q = integrator.integrate(f1, (Real) 0, (Real) 1, &error, &L1); + Real Q_expected = pi()*ln_two()/8 + catalan()*half(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + BOOST_CHECK_CLOSE(L1, Q_expected, 100*tol); + + auto f2 = [](Real x) { Real t0 = x*x + 1; Real t1 = sqrt(t0); return atan(t1)/(t0*t1); }; + Q = integrator.integrate(f2, (Real) 0 , (Real) 1, &error, &L1); + Q_expected = pi()/4 - pi()/root_two() + 3*atan(root_two())/root_two(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + BOOST_CHECK_CLOSE(L1, Q_expected, 100*tol); + + auto f5 = [](Real t) { return t*t*log(t)/((t*t - 1)*(t*t*t*t + 1)); }; + Q = integrator.integrate(f5, (Real) 0 , (Real) 1); + Q_expected = pi()*pi()*(2 - root_two())/32; + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + + // Oh it suffers on this one. + auto f6 = [](Real t) { Real x = log(t); return x*x; }; + Q = integrator.integrate(f6, (Real) 0 , (Real) 1, &error, &L1); + Q_expected = 2; + + BOOST_CHECK_CLOSE(Q, Q_expected, 5000*tol); + BOOST_CHECK_CLOSE(L1, Q_expected, 5000*tol); + + + // Although it doesn't get to the requested tolerance on this integral, the error bounds can be queried and are reasonable: + auto f7 = [](Real t) { return sqrt(tan(t)); }; + Q = integrator.integrate(f7, (Real) 0 , (Real) half_pi(), &error, &L1); + Q_expected = pi()/root_two(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + BOOST_CHECK_CLOSE(L1, Q_expected, 100*tol); + + auto f8 = [](Real t) { return log(cos(t)); }; + Q = integrator.integrate(f8, (Real) 0 , half_pi(), &error, &L1); + Q_expected = -pi()*ln_two()*half(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + BOOST_CHECK_CLOSE(L1, -Q_expected, 100*tol); +} + + +template +void test_three_quadrature_schemes_examples() +{ + std::cout << "Testing integral in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = sqrt(std::numeric_limits::epsilon()); + Real Q; + Real Q_expected; + + tanh_sinh integrator(tol, 20); + // Example 1: + auto f1 = [](Real t) { return t*log1p(t); }; + Q = integrator.integrate(f1, (Real) 0 , (Real) 1); + Q_expected = half()*half(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + + // Example 2: + auto f2 = [](Real t) { return t*t*atan(t); }; + Q = integrator.integrate(f2, (Real) 0 , (Real) 1); + Q_expected = (pi() -2 + 2*ln_two())/12; + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + // Example 3: + auto f3 = [](Real t) { return exp(t)*cos(t); }; + Q = integrator.integrate(f3, (Real) 0, half_pi()); + Q_expected = expm1(half_pi())*half(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + // Example 4: + auto f4 = [](Real x) { Real t0 = sqrt(x*x + 2); return atan(t0)/(t0*(x*x+1)); }; + Q = integrator.integrate(f4, (Real) 0 , (Real) 1); + Q_expected = 5*pi()*pi()/96; + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + // Example 5: + auto f5 = [](Real t) { return sqrt(t)*log(t); }; + Q = integrator.integrate(f5, (Real) 0 , (Real) 1); + Q_expected = -4/ (Real) 9; + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + // Example 6: + auto f6 = [](Real t) { return sqrt(1 - t*t); }; + Q = integrator.integrate(f6, (Real) 0 , (Real) 1); + Q_expected = pi()/4; + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); +} + + +template +void test_integration_over_real_line() +{ + std::cout << "Testing integrals over entire real line in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = sqrt(std::numeric_limits::epsilon()); + Real Q; + Real Q_expected; + Real error; + Real L1; + tanh_sinh integrator(tol, 20); + + auto f1 = [](Real t) { return 1/(1+t*t);}; + Q = integrator.integrate(f1, -std::numeric_limits::infinity(), std::numeric_limits::infinity(), &error, &L1); + Q_expected = pi(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + BOOST_CHECK_CLOSE(L1, Q_expected, 100*tol); + + auto f2 = [](Real t) { return exp(-t*t*half()); }; + Q = integrator.integrate(f2, -std::numeric_limits::infinity(), std::numeric_limits::infinity(), &error, &L1); + Q_expected = root_two_pi(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + BOOST_CHECK_CLOSE(L1, Q_expected, 100*tol); + + // This test shows how oscillatory integrals are approximated very poorly by this method: + //std::cout << "Testing sinc integral: \n"; + //Q = integrator.integrate(boost::math::sinc_pi, -std::numeric_limits::infinity(), std::numeric_limits::infinity(), &error, &L1); + //std::cout << "Error estimate of sinc integral: " << error << std::endl; + //std::cout << "L1 norm of sinc integral " << L1 << std::endl; + //Q_expected = pi(); + //BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + auto f4 = [](Real t) { return 1/cosh(t);}; + Q = integrator.integrate(f4, -std::numeric_limits::infinity(), std::numeric_limits::infinity(), &error, &L1); + Q_expected = pi(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + BOOST_CHECK_CLOSE(L1, Q_expected, 100*tol); + +} + +template +void test_right_limit_infinite() +{ + std::cout << "Testing right limit infinite for tanh_sinh in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = sqrt(std::numeric_limits::epsilon()); + Real Q; + Real Q_expected; + Real error; + Real L1; + tanh_sinh integrator; + + // Example 11: + auto f1 = [](Real t) { return 1/(1+t*t);}; + Q = integrator.integrate(f1, 0, std::numeric_limits::infinity(), &error, &L1); + Q_expected = half_pi(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + // Example 12 + auto f2 = [](Real t) { return exp(-t)/sqrt(t); }; + Q = integrator.integrate(f2, 0, std::numeric_limits::infinity(), &error, &L1); + Q_expected = root_pi(); + BOOST_CHECK_CLOSE(Q, Q_expected, 1000*tol); + + auto f3 = [](Real t) { return exp(-t)*cos(t); }; + Q = integrator.integrate(f3, 0, std::numeric_limits::infinity(), &error, &L1); + Q_expected = half(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + auto f4 = [](Real t) { return 1/(1+t*t); }; + Q = integrator.integrate(f4, 1, std::numeric_limits::infinity(), &error, &L1); + Q_expected = pi()/4; + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); +} + +template +void test_left_limit_infinite() +{ + std::cout << "Testing left limit infinite for tanh_sinh in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = sqrt(std::numeric_limits::epsilon()); + Real Q; + Real Q_expected; + tanh_sinh integrator; + + // Example 11: + auto f1 = [](Real t) { return 1/(1+t*t);}; + Q = integrator.integrate(f1, -std::numeric_limits::infinity(), 0); + Q_expected = half_pi(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); +} + + +// A horrible function taken from +// http://www.chebfun.org/examples/quad/GaussClenCurt.html +template +void test_horrible() +{ + std::cout << "Testing Trefenthen's horrible integral on type " << boost::typeindex::type_id().pretty_name() << "\n"; + // We only know the integral to double precision, so requesting a higher tolerance doesn't make sense. + Real tol = sqrt(std::numeric_limits::epsilon()); + Real Q; + Real Q_expected; + Real error; + Real L1; + tanh_sinh integrator; + + auto f = [](Real x) { return x*sin(2*exp(2*sin(2*exp(2*x) ) ) ); }; + Q = integrator.integrate(f, (Real) -1, (Real) 1, &error, &L1); + // The example does not have more digits than this. + // If this integration routine is correct, then all the digits here are correct. + Q_expected = 0.336732834781728; + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); +} + +// Some examples of tough integrals from NR, section 4.5.4: +template +void test_nr_examples() +{ + std::cout << "Testing singular integrals from NR 4.5.4 on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = sqrt(std::numeric_limits::epsilon()); + Real Q; + Real Q_expected; + Real error; + Real L1; + tanh_sinh integrator(tol, 15); + + auto f1 = [](Real x) { return sin(x*half())*pow(x, -3*half())*exp(-x); }; + Q = integrator.integrate(f1, 0, std::numeric_limits::infinity(), &error, &L1); + Q_expected = sqrt(pi()*(sqrt((Real) 5) - 2)); + BOOST_CHECK_CLOSE(Q, Q_expected, 1000*tol); + + auto f2 = [](Real x) { return pow(x, -(Real) 2/(Real) 7)*exp(-x*x); }; + Q = integrator.integrate(f2, 0, std::numeric_limits::infinity()); + Q_expected = half()*boost::math::tgamma((Real) 5/ (Real) 14); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + +} + +// Test integrand known to fool some termination schemes: +template +void test_early_termination() +{ + std::cout << "Testing Clenshaw & Curtis's example of integrand which fools termination schemes on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = sqrt(std::numeric_limits::epsilon()); + Real Q; + Real Q_expected; + Real error; + Real L1; + tanh_sinh integrator; + + auto f1 = [](Real x) { return 23*cosh(x)/ (Real) 25 - cos(x) ; }; + Q = integrator.integrate(f1, (Real) -1, (Real) 1, &error, &L1); + Q_expected = 46*sinh((Real) 1)/(Real) 25 - 2*sin((Real) 1); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + +} + +// Test some definite integrals from the CRC handbook, 32nd edition: +template +void test_crc() +{ + std::cout << "Testing CRC formulas on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = sqrt(std::numeric_limits::epsilon()); + Real Q; + Real Q_expected; + Real error; + Real L1; + tanh_sinh integrator; + + // CRC Definite integral 585 + auto f1 = [](Real x) { Real t = log(1/x); return x*x*t*t*t; }; + Q = integrator.integrate(f1, (Real) 0, (Real) 1, &error, &L1); + Q_expected = (Real) 2/ (Real) 27; + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + // CRC 636: + //auto f2 = [](Real x) { return sqrt(cos(x)); }; + //Q = integrator.integrate(f2, (Real) 0, (Real) pi(), &error, &L1); + //Q_expected = pow(two_pi(), 3*half())/pow(tgamma(0.25), 2); + //BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); +} + +BOOST_AUTO_TEST_CASE(tanh_sinh_quadrature_test) +{ + //generate_constants(); + test_detail(); + test_detail(); + test_detail(); + + test_right_limit_infinite(); + test_right_limit_infinite(); + test_right_limit_infinite(); + #ifdef __GNUC__ + #ifndef __clang__ + test_right_limit_infinite(); + #endif + #endif + + test_left_limit_infinite(); + test_left_limit_infinite(); + test_left_limit_infinite(); + + #ifdef __GNUC__ + #ifndef __clang__ + test_left_limit_infinite(); + #endif + #endif + + + test_linear(); + test_linear(); + test_linear(); + + #ifdef __GNUC__ + #ifndef __clang__ + test_linear(); + #endif + #endif + + test_quadratic(); + test_quadratic(); + test_quadratic(); + + #ifdef __GNUC__ + #ifndef __clang__ + test_quadratic(); + #endif + #endif + + + test_singular(); + test_singular(); + test_singular(); + + + test_ca(); + test_ca(); + test_ca(); + + #ifdef __GNUC__ + #ifndef __clang__ + test_ca(); + #endif + #endif + + + test_three_quadrature_schemes_examples(); + test_three_quadrature_schemes_examples(); + test_three_quadrature_schemes_examples(); + + #ifdef __GNUC__ + #ifndef __clang__ + test_three_quadrature_schemes_examples(); + #endif + #endif + + test_integration_over_real_line(); + test_integration_over_real_line(); + test_integration_over_real_line(); + + #ifdef __GNUC__ + #ifndef __clang__ + test_integration_over_real_line(); + #endif + #endif + + + test_horrible(); + test_horrible(); + test_horrible(); + + #ifdef __GNUC__ + #ifndef __clang__ + test_horrible(); + #endif + #endif + + + test_singular(); + test_singular(); + + test_nr_examples(); + test_nr_examples(); + test_nr_examples(); + + test_early_termination(); + test_early_termination(); + test_early_termination(); + + test_crc(); + test_crc(); + test_crc(); + +} diff --git a/test/test_barycentric_rational.cpp b/test/test_barycentric_rational.cpp new file mode 100644 index 0000000000..17b6038154 --- /dev/null +++ b/test/test_barycentric_rational.cpp @@ -0,0 +1,250 @@ +#define BOOST_TEST_MODULE barycentric_rational + +#include +#include +#include +#include +#include +#include +#include +#ifdef __GNUC__ +#ifndef __clang__ +#include +#endif +#endif + +using boost::multiprecision::cpp_bin_float_50; + +template +void test_interpolation_condition() +{ + std::cout << "Testing interpolation condition for barycentric interpolation on type " << boost::typeindex::type_id().pretty_name() << "\n"; + std::random_device rd; + std::mt19937 gen(rd()); + boost::random::uniform_real_distribution dis(0.1, 1); + std::vector x(500); + std::vector y(500); + x[0] = dis(gen); + y[0] = dis(gen); + for (size_t i = 1; i < x.size(); ++i) + { + x[i] = x[i-1] + dis(gen); + y[i] = dis(gen); + } + + boost::math::barycentric_rational interpolator(x.data(), y.data(), y.size()); + + for (size_t i = 0; i < x.size(); ++i) + { + auto z = interpolator(x[i]); + BOOST_CHECK_CLOSE(z, y[i], 100*std::numeric_limits::epsilon()); + } +} + +template +void test_interpolation_condition_high_order() +{ + std::cout << "Testing interpolation condition in high order for barycentric interpolation on type " << boost::typeindex::type_id().pretty_name() << "\n"; + std::random_device rd; + std::mt19937 gen(rd()); + boost::random::uniform_real_distribution dis(0.1, 1); + std::vector x(500); + std::vector y(500); + x[0] = dis(gen); + y[0] = dis(gen); + for (size_t i = 1; i < x.size(); ++i) + { + x[i] = x[i-1] + dis(gen); + y[i] = dis(gen); + } + + // Order 5 approximation: + boost::math::barycentric_rational interpolator(x.data(), y.data(), y.size(), 5); + + for (size_t i = 0; i < x.size(); ++i) + { + auto z = interpolator(x[i]); + BOOST_CHECK_CLOSE(z, y[i], 100*std::numeric_limits::epsilon()); + } +} + + +template +void test_constant() +{ + std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on type " << boost::typeindex::type_id().pretty_name() << "\n"; + + std::random_device rd; + std::mt19937 gen(rd()); + boost::random::uniform_real_distribution dis(0.1, 1); + std::vector x(500); + std::vector y(500); + Real constant = -8; + x[0] = dis(gen); + y[0] = constant; + for (size_t i = 1; i < x.size(); ++i) + { + x[i] = x[i-1] + dis(gen); + y[i] = y[0]; + } + + boost::math::barycentric_rational interpolator(x.data(), y.data(), y.size()); + + for (size_t i = 0; i < x.size(); ++i) + { + // Don't evaluate the constant at x[i]; that's already tested in the interpolation condition test. + auto z = interpolator(x[i] + dis(gen)); + BOOST_CHECK_CLOSE(z, constant, 100*sqrt(std::numeric_limits::epsilon())); + } +} + +template +void test_constant_high_order() +{ + std::cout << "Testing that constants are interpolated correctly in high order using barycentric interpolation on type " << boost::typeindex::type_id().pretty_name() << "\n"; + + std::random_device rd; + std::mt19937 gen(rd()); + boost::random::uniform_real_distribution dis(0.1, 1); + std::vector x(500); + std::vector y(500); + Real constant = 5; + x[0] = dis(gen); + y[0] = constant; + for (size_t i = 1; i < x.size(); ++i) + { + x[i] = x[i-1] + dis(gen); + y[i] = y[0]; + } + + // Set interpolation order to 7: + boost::math::barycentric_rational interpolator(x.data(), y.data(), y.size(), 7); + + for (size_t i = 0; i < x.size(); ++i) + { + auto z = interpolator(x[i] + dis(gen)); + BOOST_CHECK_CLOSE(z, constant, 1000*sqrt(std::numeric_limits::epsilon())); + } +} + + +template +void test_runge() +{ + std::cout << "Testing interpolation of Runge's 1/(1+25x^2) function using barycentric interpolation on type " << boost::typeindex::type_id().pretty_name() << "\n"; + + std::random_device rd; + std::mt19937 gen(rd()); + boost::random::uniform_real_distribution dis(0.005, 0.01); + std::vector x(500); + std::vector y(500); + x[0] = -2; + y[0] = 1/(1+25*x[0]*x[0]); + for (size_t i = 1; i < x.size(); ++i) + { + x[i] = x[i-1] + dis(gen); + y[i] = 1/(1+25*x[i]*x[i]); + } + + boost::math::barycentric_rational interpolator(x.data(), y.data(), y.size(), 5); + + for (size_t i = 0; i < x.size(); ++i) + { + auto t = x[i] + dis(gen); + auto z = interpolator(t); + BOOST_CHECK_CLOSE(z, 1/(1+25*t*t), 0.01); + } +} + +template +void test_weights() +{ + std::cout << "Testing weights are calculated correctly using barycentric interpolation on type " << boost::typeindex::type_id().pretty_name() << "\n"; + + std::random_device rd; + std::mt19937 gen(rd()); + boost::random::uniform_real_distribution dis(0.005, 0.01); + std::vector x(500); + std::vector y(500); + x[0] = -2; + y[0] = 1/(1+25*x[0]*x[0]); + for (size_t i = 1; i < x.size(); ++i) + { + x[i] = x[i-1] + dis(gen); + y[i] = 1/(1+25*x[i]*x[i]); + } + + boost::math::detail::barycentric_rational_imp interpolator(x.data(), y.data(), y.size(), 0); + + for (size_t i = 0; i < x.size(); ++i) + { + auto w = interpolator.weight(i); + if (i % 2 == 0) + { + BOOST_CHECK_CLOSE(w, 1, 0.00001); + } + else + { + BOOST_CHECK_CLOSE(w, -1, 0.00001); + } + } + + // d = 1: + interpolator = boost::math::detail::barycentric_rational_imp(x.data(), y.data(), y.size(), 1); + + for (size_t i = 1; i < x.size() -1; ++i) + { + auto w = interpolator.weight(i); + auto w_expect = 1/(x[i] - x[i - 1]) + 1/(x[i+1] - x[i]); + if (i % 2 == 0) + { + BOOST_CHECK_CLOSE(w, -w_expect, 0.00001); + } + else + { + BOOST_CHECK_CLOSE(w, w_expect, 0.00001); + } + } + +} + + +BOOST_AUTO_TEST_CASE(barycentric_rational) +{ + test_weights(); + test_constant(); + test_constant(); + test_constant(); + test_constant(); + + test_constant_high_order(); + test_constant_high_order(); + test_constant_high_order(); + test_constant_high_order(); + + test_interpolation_condition(); + test_interpolation_condition(); + test_interpolation_condition(); + test_interpolation_condition(); + + test_interpolation_condition_high_order(); + test_interpolation_condition_high_order(); + test_interpolation_condition_high_order(); + test_interpolation_condition_high_order(); + + test_runge(); + test_runge(); + test_runge(); + test_runge(); + + #ifdef __GNUC__ + #ifndef __clang__ + test_interpolation_condition(); + test_constant(); + test_constant_high_order(); + test_interpolation_condition_high_order(); + test_runge(); + #endif + #endif + +}