123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307 |
- // 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)
- #define BOOST_TEST_MODULE sinh_sinh_quadrature_test
- #include <complex>
- #include <boost/multiprecision/cpp_complex.hpp>
- #include <boost/math/concepts/real_concept.hpp>
- #include <boost/test/included/unit_test.hpp>
- #include <boost/test/tools/floating_point_comparison.hpp>
- #include <boost/math/quadrature/sinh_sinh.hpp>
- #include <boost/math/special_functions/sinc.hpp>
- #include <boost/multiprecision/cpp_bin_float.hpp>
- #if !BOOST_WORKAROUND(BOOST_MSVC, < 1900)
- // MSVC-12 has problems if we include 2 different multiprecision types in the same program,
- // basically random things stop compiling, even though they work fine in isolation :(
- #include <boost/multiprecision/cpp_dec_float.hpp>
- #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_quad;
- using boost::math::quadrature::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;
- //
- // Code for generating the coefficients:
- //
- template <class T>
- void print_levels(const T& v, const char* suffix)
- {
- std::cout << "{\n";
- for (unsigned i = 0; i < v.size(); ++i)
- {
- std::cout << " { ";
- for (unsigned j = 0; j < v[i].size(); ++j)
- {
- std::cout << v[i][j] << suffix << ", ";
- }
- std::cout << "},\n";
- }
- std::cout << " };\n";
- }
- template <class T>
- void print_levels(const std::pair<T, T>& p, const char* suffix = "")
- {
- std::cout << " static const std::vector<std::vector<Real> > abscissa = ";
- print_levels(p.first, suffix);
- std::cout << " static const std::vector<std::vector<Real> > weights = ";
- print_levels(p.second, suffix);
- }
- template <class Real, class TargetType>
- std::pair<std::vector<std::vector<Real>>, std::vector<std::vector<Real>> > generate_constants(unsigned max_rows)
- {
- using boost::math::constants::half_pi;
- using boost::math::constants::two_div_pi;
- using boost::math::constants::pi;
- auto g = [](Real t) { return sinh(half_pi<Real>()*sinh(t)); };
- auto w = [](Real t) { return cosh(t)*half_pi<Real>()*cosh(half_pi<Real>()*sinh(t)); };
- std::vector<std::vector<Real>> abscissa, weights;
- std::vector<Real> temp;
- Real t_max = log(2 * two_div_pi<Real>()*log(2 * two_div_pi<Real>()*sqrt(boost::math::tools::max_value<TargetType>())));
- std::cout << "m_t_max = " << t_max << ";\n";
- Real h = 1;
- for (Real t = 1; t < t_max; t += h)
- {
- temp.push_back(g(t));
- }
- abscissa.push_back(temp);
- temp.clear();
- for (Real t = 1; t < t_max; t += h)
- {
- temp.push_back(w(t * h));
- }
- weights.push_back(temp);
- temp.clear();
- for (unsigned row = 1; row < max_rows; ++row)
- {
- h /= 2;
- for (Real t = h; t < t_max; t += 2 * h)
- temp.push_back(g(t));
- abscissa.push_back(temp);
- temp.clear();
- }
- h = 1;
- for (unsigned row = 1; row < max_rows; ++row)
- {
- h /= 2;
- for (Real t = h; t < t_max; t += 2 * h)
- temp.push_back(w(t));
- weights.push_back(temp);
- temp.clear();
- }
- return std::make_pair(abscissa, weights);
- }
- template<class Real>
- void test_nr_examples()
- {
- std::cout << "Testing type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- Real integration_limit = sqrt(boost::math::tools::epsilon<Real>());
- Real tol = 10 * boost::math::tools::epsilon<Real>();
- std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
- Real Q;
- Real Q_expected;
- Real L1;
- Real error;
- sinh_sinh<Real> integrator(10);
- auto f0 = [](Real)->Real { return (Real) 0; };
- Q = integrator.integrate(f0, integration_limit, &error, &L1);
- Q_expected = 0;
- BOOST_CHECK_SMALL(Q, tol);
- BOOST_CHECK_SMALL(L1, 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 = [](const Real& t)->Real { return 1/(1+t*t); };
- Q = integrator.integrate(f1, integration_limit, &error, &L1);
- Q_expected = pi<Real>();
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
- #if defined(BOOST_MSVC) && (BOOST_MSVC < 1900)
- auto f2 = [](const Real& x)->Real { return fabs(x) > boost::math::tools::log_max_value<Real>() ? 0 : exp(-x*x); };
- #else
- auto f2 = [](const Real& x)->Real { return exp(-x*x); };
- #endif
- Q = integrator.integrate(f2, integration_limit, &error, &L1);
- Q_expected = root_pi<Real>();
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
- auto f5 = [](const Real& t)->Real { return 1/cosh(t);};
- Q = integrator.integrate(f5, integration_limit, &error, &L1);
- Q_expected = pi<Real>();
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
- // This oscillatory integral has rapid convergence because the oscillations get swamped by the exponential growth of the denominator,
- // none the less the error is slightly higher than for the other cases:
- tol *= 10;
- auto f8 = [](const Real& t)->Real { return cos(t)/cosh(t);};
- Q = integrator.integrate(f8, integration_limit, &error, &L1);
- Q_expected = pi<Real>()/cosh(half_pi<Real>());
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- // Try again with progressively fewer arguments:
- Q = integrator.integrate(f8, integration_limit);
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- Q = integrator.integrate(f8);
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- }
- // Test formulas for in the CRC Handbook of Mathematical functions, 32nd edition.
- template<class Real>
- void test_crc()
- {
- std::cout << "Testing CRC formulas on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- Real integration_limit = sqrt(boost::math::tools::epsilon<Real>());
- Real tol = 10 * boost::math::tools::epsilon<Real>();
- std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
- Real Q;
- Real Q_expected;
- Real L1;
- Real error;
- sinh_sinh<Real> integrator(10);
- // CRC Definite integral 698:
- auto f0 = [](Real x)->Real {
- using std::sinh;
- if(x == 0) {
- return (Real) 1;
- }
- return x/sinh(x);
- };
- Q = integrator.integrate(f0, integration_limit, &error, &L1);
- Q_expected = pi_sqr<Real>()/2;
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
- // CRC Definite integral 695:
- auto f1 = [](Real x)->Real {
- using std::sin; using std::sinh;
- if(x == 0) {
- return (Real) 1;
- }
- return (Real) sin(x)/sinh(x);
- };
- Q = integrator.integrate(f1, integration_limit, &error, &L1);
- Q_expected = pi<Real>()*tanh(half_pi<Real>());
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- }
- template<class Complex>
- void test_dirichlet_eta()
- {
- typedef typename Complex::value_type Real;
- std::cout << "Testing Dirichlet eta function on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- Real tol = 10 * boost::math::tools::epsilon<Real>();
- Complex Q;
- sinh_sinh<Real> integrator(10);
- //https://en.wikipedia.org/wiki/Dirichlet_eta_function, integral representations:
- Complex z = {1,1};
- auto eta = [&z](Real t)->Complex {
- using std::exp;
- using std::pow;
- using boost::math::constants::pi;
- Complex i = {0,1};
- Complex num = pow((Real)1/ (Real)2 + i*t, -z);
- Real denom = exp(pi<Real>()*t) + exp(-pi<Real>()*t);
- Complex res = num/denom;
- return res;
- };
- Q = integrator.integrate(eta);
- // N[DirichletEta[1 + I], 150]
- Complex Q_expected = {boost::lexical_cast<Real>("0.726559775062463263201495728547241386311129502735725787103568290594808442332084045617744978600192784188182345866652233650512117834307254514480657408096"),
- boost::lexical_cast<Real>("0.158095863901207324355426285544321998253687969756843115763682522207208309489794631247865357375538028170751576870244296106203144195376645765556607038775")};
- BOOST_CHECK_CLOSE_FRACTION(Q.real(), Q_expected.real(), tol);
- BOOST_CHECK_CLOSE_FRACTION(Q.imag(), Q_expected.imag(), tol);
- }
- BOOST_AUTO_TEST_CASE(sinh_sinh_quadrature_test)
- {
- //
- // Uncomment the following to print out the coefficients:
- //
- /*
- std::cout << std::scientific << std::setprecision(8);
- print_levels(generate_constants<cpp_bin_float_100, float>(8), "f");
- std::cout << std::setprecision(18);
- print_levels(generate_constants<cpp_bin_float_100, double>(8), "");
- std::cout << std::setprecision(35);
- print_levels(generate_constants<cpp_bin_float_100, cpp_bin_float_quad>(8), "L");
- */
- test_nr_examples<float>();
- test_nr_examples<double>();
- #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
- test_nr_examples<long double>();
- #endif
- test_nr_examples<cpp_bin_float_quad>();
- #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
- test_nr_examples<boost::math::concepts::real_concept>();
- #endif
- #if !BOOST_WORKAROUND(BOOST_MSVC, < 1900)
- test_nr_examples<boost::multiprecision::cpp_dec_float_50>();
- #endif
- test_crc<float>();
- test_crc<double>();
- test_dirichlet_eta<std::complex<double>>();
- #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
- test_crc<long double>();
- test_dirichlet_eta<std::complex<long double>>();
- #endif
- test_crc<cpp_bin_float_quad>();
- test_dirichlet_eta<boost::multiprecision::cpp_complex_quad>();
- #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
- test_crc<boost::math::concepts::real_concept>();
- #endif
- #if !BOOST_WORKAROUND(BOOST_MSVC, < 1900)
- test_crc<boost::multiprecision::cpp_dec_float_50>();
- #endif
- }
|