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

   1 //  (C) Copyright Nick Thompson 2020.
   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_SIMPLE_CONTINUED_FRACTION_HPP
   7 #define BOOST_MATH_TOOLS_SIMPLE_CONTINUED_FRACTION_HPP
   8
   9 #include <vector>
  10 #include <ostream>
  11 #include <iomanip>
  12 #include <cmath>
  13 #include <limits>
  14 #include <stdexcept>
  15 #include <boost/core/demangle.hpp>
  16
  17 namespace boost::math::tools {
  18
  19 template<typename Real, typename Z = int64_t>
  20 class simple_continued_fraction {
  21 public:
  22     simple_continued_fraction(Real x) : x_{x} {
  23         using std::floor;
  24         using std::abs;
  25         using std::sqrt;
  26         using std::isfinite;
  27         if (!isfinite(x)) {
  28             throw std::domain_error("Cannot convert non-finites into continued fractions.");
  29         }
  30         b_.reserve(50);
  31         Real bj = floor(x);
  32         b_.push_back(static_cast<Z>(bj));
  33         if (bj == x) {
  34            b_.shrink_to_fit();
  35            return;
  36         }
  37         x = 1/(x-bj);
  38         Real f = bj;
  39         if (bj == 0) {
  40            f = 16*std::numeric_limits<Real>::min();
  41         }
  42         Real C = f;
  43         Real D = 0;
  44         int i = 0;
  45         // the "1 + i++" lets the error bound grow slowly with the number of convergents.
  46         // I have not worked out the error propagation of the Modified Lentz's method to see if it does indeed grow at this rate.
  47         // Numerical Recipes claims that no one has worked out the error analysis of the modified Lentz's method.
  48         while (abs(f - x_) >= (1 + i++)*std::numeric_limits<Real>::epsilon()*abs(x_))
  49         {
  50           bj = floor(x);
  51           b_.push_back(static_cast<Z>(bj));
  52           x = 1/(x-bj);
  53           D += bj;
  54           if (D == 0) {
  55              D = 16*std::numeric_limits<Real>::min();
  56           }
  57           C = bj + 1/C;
  58           if (C==0) {
  59              C = 16*std::numeric_limits<Real>::min();
  60           }
  61           D = 1/D;
  62           f *= (C*D);
  63        }
  64        // Deal with non-uniqueness of continued fractions: [a0; a1, ..., an, 1] = a0; a1, ..., an + 1].
  65        // The shorter representation is considered the canonical representation,
  66        // so if we compute a non-canonical representation, change it to canonical:
  67        if (b_.size() > 2 && b_.back() == 1) {
  68           b_[b_.size() - 2] += 1;
  69           b_.resize(b_.size() - 1);
  70        }
  71        b_.shrink_to_fit();
  72
  73        for (size_t i = 1; i < b_.size(); ++i) {
  74          if (b_[i] <= 0) {
  75             std::ostringstream oss;
  76             oss << "Found a negative partial denominator: b[" << i << "] = " << b_[i] << "."
  77                 << " This means the integer type '" << boost::core::demangle(typeid(Z).name())
  78                 << "' has overflowed and you need to use a wider type,"
  79                 << " or there is a bug.";
  80             throw std::overflow_error(oss.str());
  81          }
  82        }
  83     }
  84
  85     Real khinchin_geometric_mean() const {
  86         if (b_.size() == 1) {
  87          return std::numeric_limits<Real>::quiet_NaN();
  88         }
  89          using std::log;
  90          using std::exp;
  91          // Precompute the most probable logarithms. See the Gauss-Kuzmin distribution for details.
  92          // Example: b_i = 1 has probability -log_2(3/4) ≈ .415:
  93          // A random partial denominator has ~80% chance of being in this table:
  94          const std::array<Real, 7> logs{std::numeric_limits<Real>::quiet_NaN(), Real(0), log(static_cast<Real>(2)), log(static_cast<Real>(3)), log(static_cast<Real>(4)), log(static_cast<Real>(5)), log(static_cast<Real>(6))};
  95          Real log_prod = 0;
  96          for (size_t i = 1; i < b_.size(); ++i) {
  97             if (b_[i] < static_cast<Z>(logs.size())) {
  98                log_prod += logs[b_[i]];
  99             }
 100             else
 101             {
 102                log_prod += log(static_cast<Real>(b_[i]));
 103             }
 104          }
 105          log_prod /= (b_.size()-1);
 106          return exp(log_prod);
 107     }
 108
 109     Real khinchin_harmonic_mean() const {
 110         if (b_.size() == 1) {
 111           return std::numeric_limits<Real>::quiet_NaN();
 112         }
 113         Real n = b_.size() - 1;
 114         Real denom = 0;
 115         for (size_t i = 1; i < b_.size(); ++i) {
 116             denom += 1/static_cast<Real>(b_[i]);
 117         }
 118         return n/denom;
 119     }
 120
 121     const std::vector<Z>& partial_denominators() const {
 122       return b_;
 123     }
 124
 125     template<typename T, typename Z2>
 126     friend std::ostream& operator<<(std::ostream& out, simple_continued_fraction<T, Z2>& scf);
 127
 128 private:
 129     const Real x_;
 130     std::vector<Z> b_;
 131 };
 132
 133
 134 template<typename Real, typename Z2>
 135 std::ostream& operator<<(std::ostream& out, simple_continued_fraction<Real, Z2>& scf) {
 136    constexpr const int p = std::numeric_limits<Real>::max_digits10;
 137    if constexpr (p == 2147483647) {
 138       out << std::setprecision(scf.x_.backend().precision());
 139    } else {
 140       out << std::setprecision(p);
 141    }
 142
 143    out << "[" << scf.b_.front();
 144    if (scf.b_.size() > 1)
 145    {
 146       out << "; ";
 147       for (size_t i = 1; i < scf.b_.size() -1; ++i)
 148       {
 149          out << scf.b_[i] << ", ";
 150       }
 151       out << scf.b_.back();
 152    }
 153    out << "]";
 154    return out;
 155 }
 156
 157
 158 }
 159 #endif