libs/math/test/sinh_sinh_quadrature_test.cpp

   1 // Copyright Nick Thompson, 2017
   2 // Use, modification and distribution are subject to the
   3 // Boost Software License, Version 1.0.
   4 // (See accompanying file LICENSE_1_0.txt
   5 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 #define BOOST_TEST_MODULE sinh_sinh_quadrature_test
   8 #include <complex>
   9 #include <boost/multiprecision/cpp_complex.hpp>
  10 #include <boost/math/concepts/real_concept.hpp>
  11 #include <boost/test/included/unit_test.hpp>
  12 #include <boost/test/tools/floating_point_comparison.hpp>
  13 #include <boost/math/quadrature/sinh_sinh.hpp>
  14 #include <boost/math/special_functions/sinc.hpp>
  15 #include <boost/multiprecision/cpp_bin_float.hpp>
  16
  17 #if !BOOST_WORKAROUND(BOOST_MSVC, < 1900)
  18 // MSVC-12 has problems if we include 2 different multiprecision types in the same program,
  19 // basically random things stop compiling, even though they work fine in isolation :(
  20 #include <boost/multiprecision/cpp_dec_float.hpp>
  21 #endif
  22
  23 using std::expm1;
  24 using std::exp;
  25 using std::sin;
  26 using std::cos;
  27 using std::atan;
  28 using std::tan;
  29 using std::log;
  30 using std::log1p;
  31 using std::asinh;
  32 using std::atanh;
  33 using std::sqrt;
  34 using std::isnormal;
  35 using std::abs;
  36 using std::sinh;
  37 using std::tanh;
  38 using std::cosh;
  39 using std::pow;
  40 using std::string;
  41 using boost::multiprecision::cpp_bin_float_quad;
  42 using boost::math::quadrature::sinh_sinh;
  43 using boost::math::constants::pi;
  44 using boost::math::constants::pi_sqr;
  45 using boost::math::constants::half_pi;
  46 using boost::math::constants::two_div_pi;
  47 using boost::math::constants::half;
  48 using boost::math::constants::third;
  49 using boost::math::constants::half;
  50 using boost::math::constants::third;
  51 using boost::math::constants::catalan;
  52 using boost::math::constants::ln_two;
  53 using boost::math::constants::root_two;
  54 using boost::math::constants::root_two_pi;
  55 using boost::math::constants::root_pi;
  56 //
  57 // Code for generating the coefficients:
  58 //
  59 template <class T>
  60 void print_levels(const T& v, const char* suffix)
  61 {
  62    std::cout << "{\n";
  63    for (unsigned i = 0; i < v.size(); ++i)
  64    {
  65       std::cout << "      { ";
  66       for (unsigned j = 0; j < v[i].size(); ++j)
  67       {
  68          std::cout << v[i][j] << suffix << ", ";
  69       }
  70       std::cout << "},\n";
  71    }
  72    std::cout << "   };\n";
  73 }
  74
  75 template <class T>
  76 void print_levels(const std::pair<T, T>& p, const char* suffix = "")
  77 {
  78    std::cout << "   static const std::vector<std::vector<Real> > abscissa = ";
  79    print_levels(p.first, suffix);
  80    std::cout << "   static const std::vector<std::vector<Real> > weights = ";
  81    print_levels(p.second, suffix);
  82 }
  83
  84 template <class Real, class TargetType>
  85 std::pair<std::vector<std::vector<Real>>, std::vector<std::vector<Real>> > generate_constants(unsigned max_rows)
  86 {
  87    using boost::math::constants::half_pi;
  88    using boost::math::constants::two_div_pi;
  89    using boost::math::constants::pi;
  90    auto g = [](Real t) { return sinh(half_pi<Real>()*sinh(t)); };
  91    auto w = [](Real t) { return cosh(t)*half_pi<Real>()*cosh(half_pi<Real>()*sinh(t)); };
  92
  93    std::vector<std::vector<Real>> abscissa, weights;
  94
  95    std::vector<Real> temp;
  96
  97    Real t_max = log(2 * two_div_pi<Real>()*log(2 * two_div_pi<Real>()*sqrt(boost::math::tools::max_value<TargetType>())));
  98
  99    std::cout << "m_t_max = " << t_max << ";\n";
 100
 101    Real h = 1;
 102    for (Real t = 1; t < t_max; t += h)
 103    {
 104       temp.push_back(g(t));
 105    }
 106    abscissa.push_back(temp);
 107    temp.clear();
 108
 109    for (Real t = 1; t < t_max; t += h)
 110    {
 111       temp.push_back(w(t * h));
 112    }
 113    weights.push_back(temp);
 114    temp.clear();
 115
 116    for (unsigned row = 1; row < max_rows; ++row)
 117    {
 118       h /= 2;
 119       for (Real t = h; t < t_max; t += 2 * h)
 120          temp.push_back(g(t));
 121       abscissa.push_back(temp);
 122       temp.clear();
 123    }
 124    h = 1;
 125    for (unsigned row = 1; row < max_rows; ++row)
 126    {
 127       h /= 2;
 128       for (Real t = h; t < t_max; t += 2 * h)
 129          temp.push_back(w(t));
 130       weights.push_back(temp);
 131       temp.clear();
 132    }
 133
 134    return std::make_pair(abscissa, weights);
 135 }
 136
 137 template<class Real>
 138 void test_nr_examples()
 139 {
 140     std::cout << "Testing type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 141     Real integration_limit = sqrt(boost::math::tools::epsilon<Real>());
 142     Real tol = 10 * boost::math::tools::epsilon<Real>();
 143     std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
 144     Real Q;
 145     Real Q_expected;
 146     Real L1;
 147     Real error;
 148     sinh_sinh<Real> integrator(10);
 149
 150     auto f0 = [](Real)->Real { return (Real) 0; };
 151     Q = integrator.integrate(f0, integration_limit, &error, &L1);
 152     Q_expected = 0;
 153     BOOST_CHECK_SMALL(Q, tol);
 154     BOOST_CHECK_SMALL(L1, tol);
 155
 156     // In spite of the poles at \pm i, we still get a doubling of the correct digits at each level of refinement.
 157     auto f1 = [](const Real& t)->Real { return 1/(1+t*t); };
 158     Q = integrator.integrate(f1, integration_limit, &error, &L1);
 159     Q_expected = pi<Real>();
 160     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 161     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 162 #if defined(BOOST_MSVC) && (BOOST_MSVC < 1900)
 163     auto f2 = [](const Real& x)->Real { return fabs(x) > boost::math::tools::log_max_value<Real>() ? 0 : exp(-x*x); };
 164 #else
 165     auto f2 = [](const Real& x)->Real { return exp(-x*x); };
 166 #endif
 167     Q = integrator.integrate(f2, integration_limit, &error, &L1);
 168     Q_expected = root_pi<Real>();
 169     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 170     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 171
 172     auto f5 = [](const Real& t)->Real { return 1/cosh(t);};
 173     Q = integrator.integrate(f5, integration_limit, &error, &L1);
 174     Q_expected = pi<Real>();
 175     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 176     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 177
 178     // This oscillatory integral has rapid convergence because the oscillations get swamped by the exponential growth of the denominator,
 179     // none the less the error is slightly higher than for the other cases:
 180     tol *= 10;
 181     auto f8 = [](const Real& t)->Real { return cos(t)/cosh(t);};
 182     Q = integrator.integrate(f8, integration_limit, &error, &L1);
 183     Q_expected = pi<Real>()/cosh(half_pi<Real>());
 184     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 185     // Try again with progressively fewer arguments:
 186     Q = integrator.integrate(f8, integration_limit);
 187     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 188     Q = integrator.integrate(f8);
 189     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 190 }
 191
 192 // Test formulas for in the CRC Handbook of Mathematical functions, 32nd edition.
 193 template<class Real>
 194 void test_crc()
 195 {
 196     std::cout << "Testing CRC formulas on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 197     Real integration_limit = sqrt(boost::math::tools::epsilon<Real>());
 198     Real tol = 10 * boost::math::tools::epsilon<Real>();
 199     std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
 200     Real Q;
 201     Real Q_expected;
 202     Real L1;
 203     Real error;
 204     sinh_sinh<Real> integrator(10);
 205
 206     // CRC Definite integral 698:
 207     auto f0 = [](Real x)->Real {
 208       using std::sinh;
 209       if(x == 0) {
 210         return (Real) 1;
 211       }
 212       return x/sinh(x);
 213     };
 214     Q = integrator.integrate(f0, integration_limit, &error, &L1);
 215     Q_expected = pi_sqr<Real>()/2;
 216     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 217     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 218
 219
 220     // CRC Definite integral 695:
 221     auto f1 = [](Real x)->Real {
 222       using std::sin; using std::sinh;
 223       if(x == 0) {
 224         return (Real) 1;
 225       }
 226       return (Real) sin(x)/sinh(x);
 227     };
 228     Q = integrator.integrate(f1, integration_limit, &error, &L1);
 229     Q_expected = pi<Real>()*tanh(half_pi<Real>());
 230     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 231 }
 232
 233 template<class Complex>
 234 void test_dirichlet_eta()
 235 {
 236   typedef typename Complex::value_type Real;
 237   std::cout << "Testing Dirichlet eta function on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 238   Real tol = 10 * boost::math::tools::epsilon<Real>();
 239   Complex Q;
 240   sinh_sinh<Real> integrator(10);
 241
 242   //https://en.wikipedia.org/wiki/Dirichlet_eta_function, integral representations:
 243   Complex z = {1,1};
 244   auto eta = [&z](Real t)->Complex {
 245     using std::exp;
 246     using std::pow;
 247     using boost::math::constants::pi;
 248     Complex i = {0,1};
 249     Complex num = pow((Real)1/ (Real)2  + i*t, -z);
 250     Real denom = exp(pi<Real>()*t) + exp(-pi<Real>()*t);
 251     Complex res = num/denom;
 252     return res;
 253   };
 254   Q = integrator.integrate(eta);
 255   // N[DirichletEta[1 + I], 150]
 256   Complex Q_expected = {boost::lexical_cast<Real>("0.726559775062463263201495728547241386311129502735725787103568290594808442332084045617744978600192784188182345866652233650512117834307254514480657408096"),
 257                         boost::lexical_cast<Real>("0.158095863901207324355426285544321998253687969756843115763682522207208309489794631247865357375538028170751576870244296106203144195376645765556607038775")};
 258
 259   BOOST_CHECK_CLOSE_FRACTION(Q.real(), Q_expected.real(), tol);
 260   BOOST_CHECK_CLOSE_FRACTION(Q.imag(), Q_expected.imag(), tol);
 261 }
 262
 263
 264 BOOST_AUTO_TEST_CASE(sinh_sinh_quadrature_test)
 265 {
 266     //
 267     // Uncomment the following to print out the coefficients:
 268     //
 269     /*
 270     std::cout << std::scientific << std::setprecision(8);
 271     print_levels(generate_constants<cpp_bin_float_100, float>(8), "f");
 272     std::cout << std::setprecision(18);
 273     print_levels(generate_constants<cpp_bin_float_100, double>(8), "");
 274     std::cout << std::setprecision(35);
 275     print_levels(generate_constants<cpp_bin_float_100, cpp_bin_float_quad>(8), "L");
 276     */
 277     test_nr_examples<float>();
 278     test_nr_examples<double>();
 279 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 280     test_nr_examples<long double>();
 281 #endif
 282     test_nr_examples<cpp_bin_float_quad>();
 283 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 284     test_nr_examples<boost::math::concepts::real_concept>();
 285 #endif
 286 #if !BOOST_WORKAROUND(BOOST_MSVC, < 1900)
 287     test_nr_examples<boost::multiprecision::cpp_dec_float_50>();
 288 #endif
 289
 290     test_crc<float>();
 291     test_crc<double>();
 292     test_dirichlet_eta<std::complex<double>>();
 293 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 294     test_crc<long double>();
 295     test_dirichlet_eta<std::complex<long double>>();
 296 #endif
 297     test_crc<cpp_bin_float_quad>();
 298     test_dirichlet_eta<boost::multiprecision::cpp_complex_quad>();
 299 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 300     test_crc<boost::math::concepts::real_concept>();
 301 #endif
 302 #if !BOOST_WORKAROUND(BOOST_MSVC, < 1900)
 303     test_crc<boost::multiprecision::cpp_dec_float_50>();
 304 #endif
 305
 306
 307 }