1 files changed, 139 insertions, 0 deletions
diff --git a/boost/math/tools/condition_numbers.hpp b/boost/math/tools/condition_numbers.hpp
new file mode 100644
index 0000000000..afbea75cb5
--- /dev/null
+++ b/boost/math/tools/condition_numbers.hpp
@@ -0,0 +1,139 @@
+//  (C) Copyright Nick Thompson 2019.
+//  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)
+
+#ifndef BOOST_MATH_TOOLS_CONDITION_NUMBERS_HPP
+#define BOOST_MATH_TOOLS_CONDITION_NUMBERS_HPP
+#include <cmath>
+#include <boost/math/differentiation/finite_difference.hpp>
+
+namespace boost::math::tools {
+
+template<class Real, bool kahan=true>
+class summation_condition_number {
+public:
+    summation_condition_number(Real const x = 0)
+    {
+        using std::abs;
+        m_l1 = abs(x);
+        m_sum = x;
+        m_c = 0;
+    }
+
+    void operator+=(Real const & x)
+    {
+        using std::abs;
+        // No need to Kahan the l1 calc; it's well conditioned:
+        m_l1 += abs(x);
+        if constexpr(kahan)
+        {
+            Real y = x - m_c;
+            Real t = m_sum + y;
+            m_c = (t-m_sum) -y;
+            m_sum = t;
+        }
+        else
+        {
+            m_sum += x;
+        }
+    }
+
+    inline void operator-=(Real const & x)
+    {
+        this->operator+=(-x);
+    }
+
+    // Is operator*= relevant? Presumably everything gets rescaled,
+    // (m_sum -> k*m_sum, m_l1->k*m_l1, m_c->k*m_c),
+    // but is this sensible? More important is it useful?
+    // In addition, it might change the condition number.
+
+    [[nodiscard]] Real operator()() const
+    {
+        using std::abs;
+        if (m_sum == Real(0) && m_l1 != Real(0))
+        {
+            return std::numeric_limits<Real>::infinity();
+        }
+        return m_l1/abs(m_sum);
+    }
+
+    [[nodiscard]] Real sum() const
+    {
+        // Higham, 1993, "The Accuracy of Floating Point Summation":
+        // "In [17] and [18], Kahan describes a variation of compensated summation in which the final sum is also corrected
+        // thus s=s+e is appended to the algorithm above)."
+        return m_sum + m_c;
+    }
+
+    [[nodiscard]] Real l1_norm() const
+    {
+        return m_l1;
+    }
+
+private:
+    Real m_l1;
+    Real m_sum;
+    Real m_c;
+};
+
+template<class F, class Real>
+auto evaluation_condition_number(F const & f, Real const & x)
+{
+    using std::abs;
+    using std::isnan;
+    using std::sqrt;
+    using boost::math::differentiation::finite_difference_derivative;
+
+    Real fx = f(x);
+    if (isnan(fx))
+    {
+        return std::numeric_limits<Real>::quiet_NaN();
+    }
+    bool caught_exception = false;
+    Real fp;
+    try
+    {
+        fp = finite_difference_derivative(f, x);
+    }
+    catch(...)
+    {
+        caught_exception = true;
+    }
+
+    if (isnan(fp) || caught_exception)
+    {
+        // Check if the right derivative exists:
+        fp = finite_difference_derivative<decltype(f), Real, 1>(f, x);
+        if (isnan(fp))
+        {
+            // Check if a left derivative exists:
+            const Real eps = (std::numeric_limits<Real>::epsilon)();
+            Real h = - 2 * sqrt(eps);
+            h = boost::math::differentiation::detail::make_xph_representable(x, h);
+            Real yh = f(x + h);
+            Real y0 = f(x);
+            Real diff = yh - y0;
+            fp = diff / h;
+            if (isnan(fp))
+            {
+                return std::numeric_limits<Real>::quiet_NaN();
+            }
+        }
+    }
+
+    if (fx == 0)
+    {
+        if (x==0 || fp==0)
+        {
+            return std::numeric_limits<Real>::quiet_NaN();
+        }
+        return std::numeric_limits<Real>::infinity();
+    }
+
+    return abs(x*fp/fx);
+}
+
+}
+#endif