|
diff --git a/doc/html/math_toolkit/navigation.html b/doc/html/math_toolkit/navigation.html
index 6d4285f59..54b9d334b 100644
--- a/doc/html/math_toolkit/navigation.html
+++ b/doc/html/math_toolkit/navigation.html
@@ -27,7 +27,7 @@
Navigation
-
+
Boost.Math documentation is provided in both HTML and PDF formats.
diff --git a/doc/quadrature/double_exponential.qbk b/doc/quadrature/double_exponential.qbk
index 920565718..fdd4a1088 100644
--- a/doc/quadrature/double_exponential.qbk
+++ b/doc/quadrature/double_exponential.qbk
@@ -213,6 +213,24 @@ to the square root of epsilon, and all tests were conducted in double precision:
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.
+[h4 Complex integrals]
+
+The tanh_sinh integrator supports integration of functions which return complex results, for example the sine-integral `Si(z)` has the integral representation:
+
+[equation sine_integral]
+
+Which we can code up directly as:
+
+ template
+ Complex Si(Complex z)
+ {
+ typedef typename Complex::value_type value_type;
+ using std::sin; using std::cos; using std::exp;
+ auto f = [&z](value_type t) { return -exp(-z * cos(t)) * cos(z * sin(t)); };
+ boost::math::quadrature::tanh_sinh integrator;
+ return integrator.integrate(f, 0, boost::math::constants::half_pi()) + boost::math::constants::half_pi();
+ }
+
[endsect]
[section:de_tanh_sinh_2_arg Handling functions with large features near an endpoint with tanh-sinh quadrature]
@@ -273,7 +291,27 @@ The sinh-sinh quadrature allows computation over the entire real line, and is ca
double L1;
double Q = integrator.integrate(f, &error, &L1);
-Note that the limits of integration are understood to be (-\u221E, \u221E). Complex valued integrands are supported.
+Note that the limits of integration are understood to be (-\u221E, \u221E).
+
+Complex valued integrands are supported as well, for example the [@https://en.wikipedia.org/wiki/Dirichlet_eta_function Dirichlet Eta function]
+can be represented via:
+
+[equation complex_eta_integral]
+
+which we can directly code up as:
+
+ template
+ Complex eta(Complex s)
+ {
+ typedef typename Complex::value_type value_type;
+ using std::pow; using std::exp;
+ Complex i(0, 1);
+ value_type pi = boost::math::constants::pi();
+ auto f = [&, s, i](value_type t) { return pow(0.5 + i * t, -s) / (exp(pi * t) + exp(-pi * t)); };
+ boost::math::quadrature::sinh_sinh integrator;
+ return integrator.integrate(f);
+ }
+
[endsect]
@@ -313,6 +351,32 @@ The native integration range of this integrator is (0, [infin]), but we also sup
Endpoint singularities and complex-valued integrands are supported by exp-sinh.
+For example, the modified Bessel function K can be represented via:
+
+[equation complex_bessel_k_integral]
+
+Which we can code up as:
+
+ template
+ Complex bessel_K(Complex alpha, Complex z)
+ {
+ typedef typename Complex::value_type value_type;
+ using std::cosh; using std::exp;
+ auto f = [&, alpha, z](value_type t)
+ {
+ value_type ct = cosh(t);
+ if (ct > log(std::numeric_limits::max()))
+ return Complex(0);
+ return exp(-z * ct) * cosh(alpha * t);
+ };
+ boost::math::quadrature::exp_sinh integrator;
+ return integrator.integrate(f);
+ }
+
+The only wrinkle in the above code is the need to check for large `cosh(t)` in which case we assume that
+`exp(-x cosh(t))` tends to zero faster than `cosh(alpha x)` tends to infinity and return `0`. Without that
+check we end up with `0 * Infinity` as the result (a NaN).
+
[endsect]
[section:de_tol Setting the Termination Condition for Integration]
diff --git a/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp b/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp
index e5f2bc200..b881e6b4b 100644
--- a/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp
+++ b/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp
@@ -45,7 +45,7 @@ public:
}
template
- Real integrate(const F f, Real* error, Real* L1, const char* function, Real left_min_complement, Real right_min_complement, Real tolerance, std::size_t* levels) const;
+ decltype(std::declval()(std::declval(), std::declval())) integrate(const F f, Real* error, Real* L1, const char* function, Real left_min_complement, Real right_min_complement, Real tolerance, std::size_t* levels) const;
private:
const std::vector& get_abscissa_row(std::size_t n)const
@@ -181,7 +181,7 @@ private:
template
template
-Real tanh_sinh_detail::integrate(const F f, Real* error, Real* L1, const char* function, Real left_min_complement, Real right_min_complement, Real tolerance, std::size_t* levels) const
+decltype(std::declval()(std::declval(), std::declval())) tanh_sinh_detail::integrate(const F f, Real* error, Real* L1, const char* function, Real left_min_complement, Real right_min_complement, Real tolerance, std::size_t* levels) const
{
using std::abs;
using std::fabs;
@@ -224,10 +224,12 @@ Real tanh_sinh_detail::integrate(const F f, Real* error, Real* L1,
//
BOOST_ASSERT(m_abscissas[0][max_left_position] < 0);
BOOST_ASSERT(m_abscissas[0][max_right_position] < 0);
-
+ //
+ // The type of the result:
+ typedef decltype(std::declval()(std::declval(), std::declval())) result_type;
Real h = m_t_max / m_inital_row_length;
- Real I0 = half_pi()*f(0, 1);
+ result_type I0 = half_pi()*f(0, 1);
Real L1_I0 = abs(I0);
for(size_t i = 1; i < m_abscissas[0].size(); ++i)
{
@@ -243,7 +245,7 @@ Real tanh_sinh_detail::integrate(const F f, Real* error, Real* L1,
}
else
xc = x - 1;
- Real yp, ym;
+ result_type yp, ym;
yp = i <= max_right_position ? f(x, -xc) : 0;
ym = i <= max_left_position ? f(-x, xc) : 0;
I0 += (yp + ym)*w;
@@ -257,7 +259,7 @@ Real tanh_sinh_detail::integrate(const F f, Real* error, Real* L1,
// L1_I0 and L1_I1 are the absolute integral values.
//
size_t k = 1;
- Real I1 = I0;
+ result_type I1 = I0;
Real L1_I1 = L1_I0;
Real err = 0;
//
@@ -274,7 +276,7 @@ Real tanh_sinh_detail::integrate(const F f, Real* error, Real* L1,
I1 = half()*I0;
L1_I1 = half()*L1_I0;
h *= half();
- Real sum = 0;
+ result_type sum = 0;
Real absum = 0;
auto const& abscissa_row = this->get_abscissa_row(k);
auto const& weight_row = this->get_weight_row(k);
@@ -327,9 +329,9 @@ Real tanh_sinh_detail::integrate(const F f, Real* error, Real* L1,
xc = x - 1;
}
- Real yp = j > max_right_index ? 0 : f(x, -xc);
- Real ym = j > max_left_index ? 0 : f(-x, xc);
- Real term = (yp + ym)*w;
+ result_type yp = j > max_right_index ? 0 : f(x, -xc);
+ result_type ym = j > max_left_index ? 0 : f(-x, xc);
+ result_type term = (yp + ym)*w;
sum += term;
// A question arises as to how accurately we actually need to estimate the L1 integral.
@@ -372,7 +374,7 @@ Real tanh_sinh_detail::integrate(const F f, Real* error, Real* L1,
// parameters. We could keep hunting until we find something, but that would handicap
// integrals which really are zero.... so a compromise then!
//
- if (err <= tolerance*L1_I1)
+ if (err <= abs(tolerance*L1_I1))
{
break;
}
diff --git a/include/boost/math/quadrature/tanh_sinh.hpp b/include/boost/math/quadrature/tanh_sinh.hpp
index 7305d93aa..d08b3305a 100644
--- a/include/boost/math/quadrature/tanh_sinh.hpp
+++ b/include/boost/math/quadrature/tanh_sinh.hpp
@@ -44,14 +44,14 @@ public:
: m_imp(std::make_shared>(max_refinements, min_complement)) {}
template
- auto integrate(const F f, Real a, Real b, Real tolerance = tools::root_epsilon(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) ->decltype(Real(std::declval()(std::declval()))) const;
+ auto integrate(const F f, Real a, Real b, Real tolerance = tools::root_epsilon(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) ->decltype(std::declval()(std::declval())) const;
template
- auto integrate(const F f, Real a, Real b, Real tolerance = tools::root_epsilon(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) ->decltype(Real(std::declval()(std::declval(), std::declval()))) const;
+ auto integrate(const F f, Real a, Real b, Real tolerance = tools::root_epsilon(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) ->decltype(std::declval()(std::declval(), std::declval())) const;
template
- auto integrate(const F f, Real tolerance = tools::root_epsilon(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) ->decltype(Real(std::declval()(std::declval()))) const;
+ auto integrate(const F f, Real tolerance = tools::root_epsilon(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) ->decltype(std::declval()(std::declval())) const;
template
- auto integrate(const F f, Real tolerance = tools::root_epsilon(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) ->decltype(Real(std::declval()(std::declval(), std::declval()))) const;
+ auto integrate(const F f, Real tolerance = tools::root_epsilon(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) ->decltype(std::declval()(std::declval(), std::declval())) const;
private:
std::shared_ptr> m_imp;
@@ -59,7 +59,7 @@ private:
template
template
-auto tanh_sinh::integrate(const F f, Real a, Real b, Real tolerance, Real* error, Real* L1, std::size_t* levels) ->decltype(Real(std::declval()(std::declval()))) const
+auto tanh_sinh::integrate(const F f, Real a, Real b, Real tolerance, Real* error, Real* L1, std::size_t* levels) ->decltype(std::declval()(std::declval())) const
{
BOOST_MATH_STD_USING
using boost::math::constants::half;
@@ -67,13 +67,15 @@ auto tanh_sinh::integrate(const F f, Real a, Real b, Real toleranc
static const char* function = "tanh_sinh<%1%>::integrate";
+ typedef decltype(std::declval()(std::declval())) result_type;
+
if (!(boost::math::isnan)(a) && !(boost::math::isnan)(b))
{
// Infinite limits:
if ((a <= -tools::max_value()) && (b >= tools::max_value()))
{
- auto u = [&](const Real& t, const Real& tc)->Real
+ auto u = [&](const Real& t, const Real& tc)->result_type
{
Real t_sq = t*t;
Real inv;
@@ -92,7 +94,7 @@ auto tanh_sinh::integrate(const F f, Real a, Real b, Real toleranc
// Right limit is infinite:
if ((boost::math::isfinite)(a) && (b >= tools::max_value()))
{
- auto u = [&](const Real& t, const Real& tc)->Real
+ auto u = [&](const Real& t, const Real& tc)->result_type
{
Real z, arg;
if (t > -0.5f)
@@ -106,7 +108,7 @@ auto tanh_sinh::integrate(const F f, Real a, Real b, Real toleranc
return f(arg)*z*z;
};
Real left_limit = sqrt(tools::min_value()) * 4;
- Real Q = 2 * m_imp->integrate(u, error, L1, function, left_limit, tools::min_value(), tolerance, levels);
+ result_type Q = Real(2) * m_imp->integrate(u, error, L1, function, left_limit, tools::min_value(), tolerance, levels);
if (L1)
{
*L1 *= 2;
@@ -117,7 +119,7 @@ auto tanh_sinh::integrate(const F f, Real a, Real b, Real toleranc
if ((boost::math::isfinite)(b) && (a <= -tools::max_value()))
{
- auto v = [&](const Real& t, const Real& tc)->Real
+ auto v = [&](const Real& t, const Real& tc)->result_type
{
Real z;
if (t > -0.5)
@@ -133,7 +135,7 @@ auto tanh_sinh::integrate(const F f, Real a, Real b, Real toleranc
};
Real left_limit = sqrt(tools::min_value()) * 4;
- Real Q = 2 * m_imp->integrate(v, error, L1, function, left_limit, tools::min_value(), tolerance, levels);
+ result_type Q = Real(2) * m_imp->integrate(v, error, L1, function, left_limit, tools::min_value(), tolerance, levels);
if (L1)
{
*L1 *= 2;
@@ -169,7 +171,7 @@ auto tanh_sinh::integrate(const F f, Real a, Real b, Real toleranc
//
BOOST_ASSERT((left_min_complement * diff + a) > a);
BOOST_ASSERT((b - right_min_complement * diff) < b);
- auto u = [&](Real z, Real zc)->Real
+ auto u = [&](Real z, Real zc)->result_type
{
Real position;
if (z < -0.5)
@@ -190,7 +192,7 @@ auto tanh_sinh::integrate(const F f, Real a, Real b, Real toleranc
BOOST_ASSERT(position != b);
return f(position);
};
- Real Q = diff*m_imp->integrate(u, error, L1, function, left_min_complement, right_min_complement, tolerance, levels);
+ result_type Q = diff*m_imp->integrate(u, error, L1, function, left_min_complement, right_min_complement, tolerance, levels);
if (L1)
{
@@ -204,7 +206,7 @@ auto tanh_sinh::integrate(const F f, Real a, Real b, Real toleranc
template
template
-auto tanh_sinh::integrate(const F f, Real a, Real b, Real tolerance, Real* error, Real* L1, std::size_t* levels) ->decltype(Real(std::declval()(std::declval(), std::declval()))) const
+auto tanh_sinh::integrate(const F f, Real a, Real b, Real tolerance, Real* error, Real* L1, std::size_t* levels) ->decltype(std::declval()(std::declval(), std::declval())) const
{
BOOST_MATH_STD_USING
using boost::math::constants::half;
@@ -241,7 +243,7 @@ auto tanh_sinh::integrate(const F f, Real a, Real b, Real toleranc
template
template
-auto tanh_sinh::integrate(const F f, Real tolerance, Real* error, Real* L1, std::size_t* levels) ->decltype(Real(std::declval()(std::declval()))) const
+auto tanh_sinh::integrate(const F f, Real tolerance, Real* error, Real* L1, std::size_t* levels) ->decltype(std::declval()(std::declval())) const
{
using boost::math::quadrature::detail::tanh_sinh_detail;
static const char* function = "tanh_sinh<%1%>::integrate";
@@ -251,7 +253,7 @@ auto tanh_sinh::integrate(const F f, Real tolerance, Real* error,
template
template
-auto tanh_sinh::integrate(const F f, Real tolerance, Real* error, Real* L1, std::size_t* levels) ->decltype(Real(std::declval()(std::declval(), std::declval()))) const
+auto tanh_sinh::integrate(const F f, Real tolerance, Real* error, Real* L1, std::size_t* levels) ->decltype(std::declval()(std::declval(), std::declval())) const
{
using boost::math::quadrature::detail::tanh_sinh_detail;
static const char* function = "tanh_sinh<%1%>::integrate";
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index dfe12f163..eb101eb75 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -1057,6 +1057,10 @@ test-suite misc :
: : : release msvc:/bigobj TEST8 [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ]
[ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr cxx11_unified_initialization_syntax sfinae_expr ] :
tanh_sinh_quadrature_test_8 ]
+ [ run tanh_sinh_quadrature_test.cpp ../../test/build//boost_unit_test_framework
+ : : : release msvc:/bigobj TEST9
+ [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr cxx11_unified_initialization_syntax sfinae_expr ] :
+ tanh_sinh_quadrature_test_9 ]
[ run sinh_sinh_quadrature_test.cpp ../../test/build//boost_unit_test_framework
: : : release [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr cxx11_unified_initialization_syntax ] ]
diff --git a/test/tanh_sinh_quadrature_test.cpp b/test/tanh_sinh_quadrature_test.cpp
index 8fbc62618..beaa06526 100644
--- a/test/tanh_sinh_quadrature_test.cpp
+++ b/test/tanh_sinh_quadrature_test.cpp
@@ -33,7 +33,7 @@
#endif
#if !defined(TEST1) && !defined(TEST2) && !defined(TEST3) && !defined(TEST4) && !defined(TEST5) && !defined(TEST6) && !defined(TEST7) && !defined(TEST8)\
- && !defined(TEST1A) && !defined(TEST2A) && !defined(TEST3A) && !defined(TEST6A)
+ && !defined(TEST1A) && !defined(TEST2A) && !defined(TEST3A) && !defined(TEST6A) && !defined(TEST9)
# define TEST1
# define TEST2
# define TEST3
@@ -46,6 +46,7 @@
# define TEST2A
# define TEST3A
# define TEST6A
+# define TEST9
#endif
using std::expm1;
@@ -563,7 +564,7 @@ void test_crc()
// We also use a 2 argument functor so that 1-x is evaluated accurately:
if (std::numeric_limits::max_exponent > std::numeric_limits::max_exponent)
{
- for (Real p = -0.99; p < 1; p += 0.1) {
+ for (Real p = Real (-0.99); p < 1; p += Real(0.1)) {
auto f = [&](Real x, Real cx)->Real
{
//return pow(x, p) / pow(1 - x, p);
@@ -579,7 +580,7 @@ void test_crc()
// domain (0, INF). Internally we need to expand out the terms and evaluate using logs to avoid spurous overflow,
// this gives us
// for p > 0:
- for (Real p = 0.99; p > 0; p -= 0.1) {
+ for (Real p = Real(0.99); p > 0; p -= Real(0.1)) {
auto f = [&](Real x)->Real
{
return exp(-x * (1 - p) + p * log(-boost::math::expm1(-x)));
@@ -589,7 +590,7 @@ void test_crc()
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 10 * tol);
}
// and for p < 1:
- for (Real p = -0.99; p < 0; p += 0.1) {
+ for (Real p = Real (-0.99); p < 0; p += Real(0.1)) {
auto f = [&](Real x)->Real
{
return exp(-p * log(-boost::math::expm1(-x)) - (1 + p) * x);
@@ -616,7 +617,7 @@ void test_crc()
//
Real limit = std::numeric_limits::max_exponent > std::numeric_limits::max_exponent
? .95f : .85f;
- for (Real h = 0.01; h < limit; h += 0.1) {
+ for (Real h = Real(0.01); h < limit; h += Real(0.1)) {
auto f = [&](Real x, Real xc)->Real { return xc > 0 ? pow(1/tan(xc), h) : pow(tan(x), h); };
Q = integrator.integrate(f, (Real) 0, half_pi(), get_convergence_tolerance(), &error, &L1);
Q_expected = half_pi()/cos(h*half_pi());
@@ -706,7 +707,7 @@ void test_sf()
Real x = 2, y = 3, z = 0.5, p = 0.25;
// Elliptic integral RC:
- Real Q = integrator.integrate([&](const Real& t) { return 1 / (sqrt(t + x) * (t + y)); }, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value()) / 2;
+ Real Q = integrator.integrate([&](const Real& t)->Real { return 1 / (sqrt(t + x) * (t + y)); }, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value()) / 2;
BOOST_CHECK_CLOSE_FRACTION(Q, boost::math::ellint_rc(x, y), tol);
// Elliptic Integral RJ:
BOOST_CHECK_CLOSE_FRACTION(Real((Real(3) / 2) * integrator.integrate([&](Real t)->Real { return 1 / (sqrt((t + x) * (t + y) * (t + z)) * (t + p)); }, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value())), boost::math::ellint_rj(x, y, z, p), tol);
@@ -719,7 +720,7 @@ void test_sf()
// tgamma expressed as an integral:
BOOST_CHECK_CLOSE_FRACTION(integrator.integrate([&](Real t)->Real { using std::pow; using std::exp; return t > 10000 ? Real(0) : Real(pow(t, z - 1) * exp(-t)); }, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value()),
boost::math::tgamma(z), tol);
- BOOST_CHECK_CLOSE_FRACTION(integrator.integrate([](const Real& t) { using std::exp; return exp(-t*t); }, std::numeric_limits::has_infinity ? -std::numeric_limits::infinity() : -boost::math::tools::max_value(), std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value()),
+ BOOST_CHECK_CLOSE_FRACTION(integrator.integrate([](const Real& t)->Real { using std::exp; return exp(-t*t); }, std::numeric_limits::has_infinity ? -std::numeric_limits::infinity() : -boost::math::tools::max_value(), std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value()),
boost::math::constants::root_pi(), tol);
}
@@ -788,6 +789,50 @@ void test_2_arg()
BOOST_CHECK_CLOSE_FRACTION(L1, -Q_expected, tol);
}
+template
+void test_complex()
+{
+ typedef typename Complex::value_type value_type;
+ //
+ // Integral version of the confluent hypergeometric function:
+ // https://dlmf.nist.gov/13.4#E1
+ //
+ value_type tol = 10 * boost::math::tools::epsilon();
+ Complex a(2, 3), b(3, 4), z(0.5, -2);
+
+ auto f = [&](value_type t)
+ {
+ return exp(z * t) * pow(t, a - value_type(1)) * pow(value_type(1) - t, b - a - value_type(1));
+ };
+
+ auto integrator = get_integrator();
+ auto Q = integrator.integrate(f, value_type(0), value_type(1), get_convergence_tolerance());
+ //
+ // Expected result computed from http://www.wolframalpha.com/input/?i=1F1%5B(2%2B3i),+(3%2B4i);+(0.5-2i)%5D+*+gamma(2%2B3i)+*+gamma(1%2Bi)+%2F+gamma(3%2B4i)
+ //
+ Complex expected(boost::lexical_cast("-0.2911081612888249710582867318081776512805281815037891183828405999609246645054069649838607112484426042883371996"),
+ boost::lexical_cast("0.4507983563969959578849120188097153649211346293694903758252662015991543519595834937475296809912196906074655385"));
+
+ value_type error = abs(expected - Q);
+ BOOST_CHECK_LE(error, tol);
+
+ //
+ // Sin Integral https://dlmf.nist.gov/6.2#E9
+ //
+ auto f2 = [z](value_type t)
+ {
+ return -exp(-z * cos(t)) * cos(z * sin(t));
+ };
+ Q = integrator.integrate(f2, value_type(0), boost::math::constants::half_pi(), get_convergence_tolerance());
+
+ expected = Complex(boost::lexical_cast("0.8893822921008980697856313681734926564752476188106405688951257340480164694708337246829840859633322683740376134733"),
+ -boost::lexical_cast("2.381380802906111364088958767973164614925936185337231718483495612539455538280372745733208000514737758457795502168"));
+ expected -= boost::math::constants::half_pi();
+
+ error = abs(expected - Q);
+ BOOST_CHECK_LE(error, tol);
+}
+
BOOST_AUTO_TEST_CASE(tanh_sinh_quadrature_test)
{
@@ -927,6 +972,12 @@ BOOST_AUTO_TEST_CASE(tanh_sinh_quadrature_test)
test_2_arg();
#endif
+#ifdef TEST9
+ test_complex >();
+ test_complex >();
+#endif
+
+
#endif
}
+ 
