...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
#include <boost/math/tools/condition_numbers.hpp> namespace boost::math::tools { template<class Real, bool kahan=true> class summation_condition_number { public: summation_condition_number(Real const x = 0); void operator+=(Real const & x); inline void operator-=(Real const & x); [[nodiscard]] Real operator()() const; [[nodiscard]] Real sum() const; [[nodiscard]] Real l1_norm() const; }; template<class F, class Real> auto evaluation_condition_number(F const & f, Real const & x); } // namespaces
Here we compute the condition number of the alternating harmonic sum:
using boost::math::tools::summation_condition_number; auto cond = summation_condition_number<float, /* kahan = */ false>(); float max_n = 10000000; for (float n = 1; n < max_n; n += 2) { cond += 1/n; cond -= 1/(n+1); } std::cout << std::setprecision(std::numeric_limits<float>::digits10); std::cout << "ln(2) = " << boost::math::constants::ln_two<float>() << "\n"; std::cout << "Sum = " << cond.sum() << "\n"; std::cout << "Condition number = " << cond() << "\n";
Output:
ln(2) = 0.693147 Sum = 0.693137 Condition number = 22.22282
We see that we have lost roughly two digits of accuracy, consistent with the heuristic that if the condition number is 10^{k}, then we lose k significant digits in the sum.
Our guess it that if you're worried about whether your sum is ill-conditioned, the last thing you want is for your condition number estimate to be inaccurate as well. Since the condition number estimate relies on computing the (perhaps ill-conditioned) sum, we have defaulted the accumulation to use Kahan summation:
auto cond = boost::math::tools::summation_condition_number<float>(); // will use Kahan summation. // ...
Output:
ln(2) = 0.693147 Kahan sum = 0.693147 Condition number = 22.2228
If you are interested, the L_{1} norm is also generated by this computation, so you may query it if you like:
float l1 = cond.l1_norm(); // l1 = 15.4
The condition number of function evaluation is defined as the absolute value of xf'(x)f(x)^{-1}. It is computed as follows:
using boost::math::tools::evaluation_condition_number; auto f = [](double x)->double { return std::log(x); }; double x = 1.125; double cond = evaluation_condition_number(f, 1.125); // cond = 1/log(x)
The condition number of function evaluation gives us a measure of how sensitive our function is to roundoff error. Unfortunately, evaluating the condition number requires evaluating the function and its derivative, and this calculation is itself inaccurate whenever the condition number of function evaluation is large. Sadly, this is also the regime when you are most interested in evaluating a condition number!
This seems to be a fundamental problem. However, it should not be necessary to evaluate the condition number to high accuracy, valuable insights can be obtained simply by looking at the change in condition number as the function evolves over its domain.