ceph/src/boost/boost/math/tools/condition_numbers.hpp

   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:
  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
  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
  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
  52     [[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
  62     [[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
  70     [[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>
  82 Real 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