diff --git a/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp b/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp index 7cef821aae..b36ccf54ad 100644 --- a/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp +++ b/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp @@ -195,6 +195,14 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin using std::sqrt; using boost::math::constants::half; using boost::math::constants::half_pi; + + // + // The type of the result: + typedef decltype(std::declval()(std::declval(), std::declval())) result_type; + + Real h = m_t_max / m_inital_row_length; + result_type I0 = half_pi() * f(0, 1); + Real L1_I0 = abs(I0); // // We maintain 4 integer values: // max_left_position is the logical index of the abscissa value closest to the @@ -222,22 +230,42 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin while (max_right_position && fabs(m_abscissas[0][max_right_position]) < right_min_complement) --max_right_position; // + // Check for non-finite values at the end points: + // + result_type yp, ym; + do + { + yp = f(-1 - m_abscissas[0][max_left_position], m_abscissas[0][max_left_position]); + if ((boost::math::isfinite)(yp)) + break; + --max_left_position; + } while (m_abscissas[0][max_left_position] < 0); + do + { + ym = f(1 + m_abscissas[0][max_right_position], -m_abscissas[0][max_right_position]); + if ((boost::math::isfinite)(ym)) + break; + --max_right_position; + } while (m_abscissas[0][max_right_position] < 0); + + if ((max_left_position == 0) && (max_right_position == 0)) + { + return policies::raise_evaluation_error(function, "The tanh_sinh quadrature found your function to be non-finite everywhere! Please check your function for singularities.", ym, Policy()); + } + + I0 += yp * m_weights[0][max_left_position] + ym * m_weights[0][max_right_position]; + L1_I0 += abs(yp * m_weights[0][max_left_position]) + abs(ym * m_weights[0][max_right_position]); + // // Assumption: left_min_complement/right_min_complement are sufficiently small that we only // ever decrement through the stored values that are complements (the negative ones), and // never ever hit the true abscissa values (positive stored values). // BOOST_MATH_ASSERT(m_abscissas[0][max_left_position] < 0); BOOST_MATH_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; - result_type I0 = half_pi()*f(0, 1); - Real L1_I0 = abs(I0); for(size_t i = 1; i < m_abscissas[0].size(); ++i) { - if ((i > max_right_position) && (i > max_left_position)) + if ((i >= max_right_position) && (i >= max_left_position)) break; Real x = m_abscissas[0][i]; Real xc = x; @@ -249,9 +277,8 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin } else xc = x - 1; - result_type yp, ym; - yp = i <= max_right_position ? f(x, -xc) : 0; - ym = i <= max_left_position ? f(-x, xc) : 0; + yp = i < max_right_position ? f(x, -xc) : 0; + ym = i < max_left_position ? f(-x, xc) : 0; I0 += (yp + ym)*w; L1_I0 += (abs(yp) + abs(ym))*w; } @@ -282,6 +309,7 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin h *= half(); result_type sum = 0; Real absum = 0; + Real endpoint_error = 0; auto const& abscissa_row = this->get_abscissa_row(k); auto const& weight_row = this->get_weight_row(k); std::size_t first_complement_index = this->get_first_complement_index(k); @@ -311,12 +339,45 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin ++max_right_position; ++max_right_index; } + // + // We also check that our endpoints don't hit singularities: + // + do + { + yp = f(-1 - abscissa_row[max_left_index], abscissa_row[max_left_index]); + if ((boost::math::isfinite)(yp)) + break; + max_left_position -= 2; + --max_left_index; + } while (abscissa_row[max_left_index] < 0); + do + { + ym = f(1 + abscissa_row[max_right_index], -abscissa_row[max_right_index]); + if ((boost::math::isfinite)(ym)) + break; + --max_right_index; + max_right_position -= 2; + } while (abscissa_row[max_right_index] < 0); + + sum += yp * weight_row[max_left_index] + ym * weight_row[max_right_index]; + absum += abs(yp * weight_row[max_left_index]) + abs(ym * weight_row[max_right_index]); + // + // We estimate the error due to truncation as the value contributed by the two most extreme points. + // In most cases this is tiny and can be ignored, if it is significant then either the area of the + // integral is so far our in the tails that our exponent range can't reach it (example x^-p at double + // precision and p ~ 1), or our function is truncated near epsilon, and we have had to narrow our endpoints. + // In this latter case we may over-estimate the error, but this is the best we can do. + // In any event, we do not add endpoint_error to the error estimate until we terminate the main loop, + // otherwise it can make things appear to be non-converged, when in reality, they are as converged as they + // will ever be. + // + endpoint_error = absum; for(size_t j = 0; j < weight_row.size(); ++j) { // If both left and right abscissa values are out of bounds at this step // we can just stop this loop right now: - if ((j > max_left_index) && (j > max_right_index)) + if ((j >= max_left_index) && (j >= max_right_index)) break; Real x = abscissa_row[j]; Real xc = x; @@ -333,8 +394,8 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin xc = x - 1; } - result_type yp = j > max_right_index ? 0 : f(x, -xc); - result_type ym = j > max_left_index ? 0 : f(-x, xc); + yp = j >= max_right_index ? 0 : f(x, -xc); + ym = j >= max_left_index ? 0 : f(-x, xc); result_type term = (yp + ym)*w; sum += term; @@ -359,15 +420,38 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin // // If the error is increasing, and we're past level 4, something bad is very likely happening: // - if ((err > last_err) && (k > 4) && (++thrash_count > 1)) + if ((err * 1.5 > last_err) && (k > 4)) { - // We could raise an evaluation_error, but since we likely have some sort of result, just return the last one - // (ie before the error started going up) - I1 = I0; - L1_I1 = L1_I0; - --k; - err = last_err; - break; + bool terminate = false; + if ((++thrash_count > 1) && (last_err < 1e-3)) + // Probably just thrashing, abort: + terminate = true; + else if(thrash_count > 2) + // OK, terrible error, but giving up anyway! + terminate = true; + else if (last_err < boost::math::tools::root_epsilon()) + // Trying to squeeze precision that probably isn't there, abort: + terminate = true; + else + { + // Take a look at the end points, if there's significant new area being + // discovered, then we're not able to get close enough to the endpoints + // to ever find the integral: + if (abs(endpoint_error / sum) > err) + terminate = true; + } + + if (terminate) + { + // We could raise an evaluation_error, but since we likely have some sort of result, just return the last one + // (ie before the error started going up) + I1 = I0; + L1_I1 = L1_I0; + --k; + err = last_err + endpoint_error; + break; + } + // Fall through and keep going, assume we've discovered a new feature of f(x).... } // // Termination condition: @@ -380,6 +464,31 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin // if (err <= abs(tolerance*L1_I1)) { + // + // A quick sanity check: have we at some point narrowed our boundaries as a result + // of non-finite values? If so let's check that the area isn't on an increasing + // trajectory at our new end point, and increase our error estimate by the last + // good value as an estimate for what we may have discarded. + // + if ((max_left_index < abscissa_row.size() - 1) && (abs(abscissa_row[max_left_index + 1]) > left_min_complement)) + { + yp = f(-1 - abscissa_row[max_left_index], abscissa_row[max_left_index]) * weight_row[max_left_index]; + ym = f(-1 - abscissa_row[max_left_index - 1], abscissa_row[max_left_index - 1]) * weight_row[max_left_index - 1]; + if (abs(yp) > abs(ym)) + { + return policies::raise_evaluation_error(function, "The tanh_sinh quadrature evaluated your function at a singular point and got %1%. Inetgration bounds were automatically narrowed, but the integral was found to be increasing at the new endpoint. Please check your function, and consider providing a 2-argument functor.", I1, Policy()); + } + } + if ((max_right_index < abscissa_row.size() - 1) && (abs(abscissa_row[max_right_index + 1]) > right_min_complement)) + { + yp = f(-1 - abscissa_row[max_right_index], abscissa_row[max_right_index]) * weight_row[max_right_index]; + ym = f(-1 - abscissa_row[max_right_index - 1], abscissa_row[max_right_index - 1]) * weight_row[max_right_index - 1]; + if (abs(yp) > abs(ym)) + { + return policies::raise_evaluation_error(function, "The tanh_sinh quadrature evaluated your function at a singular point and got %1%. Inetgration bounds were automatically narrowed, but the integral was found to be increasing at the new endpoint. Please check your function, and consider providing a 2-argument functor.", I1, Policy()); + } + } + err += endpoint_error; break; }