145 lines · plain
1// (C) Copyright Nick Thompson 2019.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_TOOLS_CONDITION_NUMBERS_HPP7#define BOOST_MATH_TOOLS_CONDITION_NUMBERS_HPP8#include <cmath>9#include <limits>10#include <boost/math/differentiation/finite_difference.hpp>11#include <boost/math/tools/config.hpp>12 13namespace boost { namespace math { namespace tools {14 15template<class Real, bool kahan=true>16class summation_condition_number {17public:18 summation_condition_number(Real const x = 0)19 {20 using std::abs;21 m_l1 = abs(x);22 m_sum = x;23 m_c = 0;24 }25 26 void operator+=(Real const & x)27 {28 using std::abs;29 // No need to Kahan the l1 calc; it's well conditioned:30 m_l1 += abs(x);31 BOOST_MATH_IF_CONSTEXPR (kahan)32 {33 Real y = x - m_c;34 Real t = m_sum + y;35 m_c = (t-m_sum) -y;36 m_sum = t;37 }38 else39 {40 m_sum += x;41 }42 }43 44 inline void operator-=(Real const & x)45 {46 this->operator+=(-x);47 }48 49 // Is operator*= relevant? Presumably everything gets rescaled,50 // (m_sum -> k*m_sum, m_l1->k*m_l1, m_c->k*m_c),51 // but is this sensible? More important is it useful?52 // In addition, it might change the condition number.53 54 Real operator()() const55 {56 using std::abs;57 if (m_sum == Real(0) && m_l1 != Real(0))58 {59 return std::numeric_limits<Real>::infinity();60 }61 return m_l1/abs(m_sum);62 }63 64 Real sum() const65 {66 // Higham, 1993, "The Accuracy of Floating Point Summation":67 // "In [17] and [18], Kahan describes a variation of compensated summation in which the final sum is also corrected68 // thus s=s+e is appended to the algorithm above)."69 return m_sum + m_c;70 }71 72 Real l1_norm() const73 {74 return m_l1;75 }76 77private:78 Real m_l1;79 Real m_sum;80 Real m_c;81};82 83template<class F, class Real>84Real evaluation_condition_number(F const & f, Real const & x)85{86 using std::abs;87 using std::isnan;88 using std::sqrt;89 using boost::math::differentiation::finite_difference_derivative;90 91 Real fx = f(x);92 if (isnan(fx))93 {94 return std::numeric_limits<Real>::quiet_NaN();95 }96 bool caught_exception = false;97 Real fp;98#ifndef BOOST_MATH_NO_EXCEPTIONS99 try100 {101#endif102 fp = finite_difference_derivative(f, x);103#ifndef BOOST_MATH_NO_EXCEPTIONS104 }105 catch(...)106 {107 caught_exception = true;108 }109#endif110 if (isnan(fp) || caught_exception)111 {112 // Check if the right derivative exists:113 fp = finite_difference_derivative<decltype(f), Real, 1>(f, x);114 if (isnan(fp))115 {116 // Check if a left derivative exists:117 const Real eps = (std::numeric_limits<Real>::epsilon)();118 Real h = - 2 * sqrt(eps);119 h = boost::math::differentiation::detail::make_xph_representable(x, h);120 Real yh = f(x + h);121 Real y0 = f(x);122 Real diff = yh - y0;123 fp = diff / h;124 if (isnan(fp))125 {126 return std::numeric_limits<Real>::quiet_NaN();127 }128 }129 }130 131 if (fx == 0)132 {133 if (x==0 || fp==0)134 {135 return std::numeric_limits<Real>::quiet_NaN();136 }137 return std::numeric_limits<Real>::infinity();138 }139 140 return abs(x*fp/fx);141}142 143}}} // Namespaces144#endif145