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