math/tools/simple_continued_fraction.hpp

0001 //  (C) Copyright Nick Thompson 2020.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_TOOLS_SIMPLE_CONTINUED_FRACTION_HPP
0007 #define BOOST_MATH_TOOLS_SIMPLE_CONTINUED_FRACTION_HPP
0008
0009 #include <array>
0010 #include <vector>
0011 #include <ostream>
0012 #include <iomanip>
0013 #include <cmath>
0014 #include <cstdint>
0015 #include <limits>
0016 #include <stdexcept>
0017 #include <sstream>
0018
0019 #include <boost/math/tools/is_standalone.hpp>
0020 #ifndef BOOST_MATH_STANDALONE
0021 #include <boost/config.hpp>
0022 #ifdef BOOST_NO_CXX17_IF_CONSTEXPR
0023 #error "The header <boost/math/norms.hpp> can only be used in C++17 and later."
0024 #endif
0025 #endif
0026
0027 #ifndef BOOST_MATH_STANDALONE
0028 #include <boost/core/demangle.hpp>
0029 #endif
0030
0031 namespace boost::math::tools {
0032
0033 template<typename Real, typename Z = int64_t>
0034 class simple_continued_fraction {
0035 public:
0036     simple_continued_fraction(Real x) : x_{x} {
0037         using std::floor;
0038         using std::abs;
0039         using std::sqrt;
0040         using std::isfinite;
0041         if (!isfinite(x)) {
0042             throw std::domain_error("Cannot convert non-finites into continued fractions.");
0043         }
0044         b_.reserve(50);
0045         Real bj = floor(x);
0046         b_.push_back(static_cast<Z>(bj));
0047         if (bj == x) {
0048            b_.shrink_to_fit();
0049            return;
0050         }
0051         x = 1/(x-bj);
0052         Real f = bj;
0053         if (bj == 0) {
0054            f = 16*(std::numeric_limits<Real>::min)();
0055         }
0056         Real C = f;
0057         Real D = 0;
0058         int i = 0;
0059         // the "1 + i++" lets the error bound grow slowly with the number of convergents.
0060         // I have not worked out the error propagation of the Modified Lentz's method to see if it does indeed grow at this rate.
0061         // Numerical Recipes claims that no one has worked out the error analysis of the modified Lentz's method.
0062         while (abs(f - x_) >= (1 + i++)*std::numeric_limits<Real>::epsilon()*abs(x_))
0063         {
0064           bj = floor(x);
0065           b_.push_back(static_cast<Z>(bj));
0066           x = 1/(x-bj);
0067           D += bj;
0068           if (D == 0) {
0069              D = 16*(std::numeric_limits<Real>::min)();
0070           }
0071           C = bj + 1/C;
0072           if (C==0) {
0073              C = 16*(std::numeric_limits<Real>::min)();
0074           }
0075           D = 1/D;
0076           f *= (C*D);
0077        }
0078        // Deal with non-uniqueness of continued fractions: [a0; a1, ..., an, 1] = a0; a1, ..., an + 1].
0079        // The shorter representation is considered the canonical representation,
0080        // so if we compute a non-canonical representation, change it to canonical:
0081        if (b_.size() > 2 && b_.back() == 1) {
0082           b_[b_.size() - 2] += 1;
0083           b_.resize(b_.size() - 1);
0084        }
0085        b_.shrink_to_fit();
0086
0087        for (size_t i = 1; i < b_.size(); ++i) {
0088          if (b_[i] <= 0) {
0089             std::ostringstream oss;
0090             oss << "Found a negative partial denominator: b[" << i << "] = " << b_[i] << "."
0091                 #ifndef BOOST_MATH_STANDALONE
0092                 << " This means the integer type '" << boost::core::demangle(typeid(Z).name())
0093                 #else
0094                 << " This means the integer type '" << typeid(Z).name()
0095                 #endif
0096                 << "' has overflowed and you need to use a wider type,"
0097                 << " or there is a bug.";
0098             throw std::overflow_error(oss.str());
0099          }
0100        }
0101     }
0102
0103     Real khinchin_geometric_mean() const {
0104         if (b_.size() == 1) {
0105          return std::numeric_limits<Real>::quiet_NaN();
0106         }
0107          using std::log;
0108          using std::exp;
0109          // Precompute the most probable logarithms. See the Gauss-Kuzmin distribution for details.
0110          // Example: b_i = 1 has probability -log_2(3/4) ~ .415:
0111          // A random partial denominator has ~80% chance of being in this table:
0112          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))};
0113          Real log_prod = 0;
0114          for (size_t i = 1; i < b_.size(); ++i) {
0115             if (b_[i] < static_cast<Z>(logs.size())) {
0116                log_prod += logs[b_[i]];
0117             }
0118             else
0119             {
0120                log_prod += log(static_cast<Real>(b_[i]));
0121             }
0122          }
0123          log_prod /= (b_.size()-1);
0124          return exp(log_prod);
0125     }
0126
0127     Real khinchin_harmonic_mean() const {
0128         if (b_.size() == 1) {
0129           return std::numeric_limits<Real>::quiet_NaN();
0130         }
0131         Real n = b_.size() - 1;
0132         Real denom = 0;
0133         for (size_t i = 1; i < b_.size(); ++i) {
0134             denom += 1/static_cast<Real>(b_[i]);
0135         }
0136         return n/denom;
0137     }
0138
0139     const std::vector<Z>& partial_denominators() const {
0140       return b_;
0141     }
0142
0143     template<typename T, typename Z2>
0144     friend std::ostream& operator<<(std::ostream& out, simple_continued_fraction<T, Z2>& scf);
0145
0146 private:
0147     const Real x_;
0148     std::vector<Z> b_;
0149 };
0150
0151
0152 template<typename Real, typename Z2>
0153 std::ostream& operator<<(std::ostream& out, simple_continued_fraction<Real, Z2>& scf) {
0154    constexpr const int p = std::numeric_limits<Real>::max_digits10;
0155    if constexpr (p == 2147483647) {
0156       out << std::setprecision(scf.x_.backend().precision());
0157    } else {
0158       out << std::setprecision(p);
0159    }
0160
0161    out << "[" << scf.b_.front();
0162    if (scf.b_.size() > 1)
0163    {
0164       out << "; ";
0165       for (size_t i = 1; i < scf.b_.size() -1; ++i)
0166       {
0167          out << scf.b_[i] << ", ";
0168       }
0169       out << scf.b_.back();
0170    }
0171    out << "]";
0172    return out;
0173 }
0174
0175
0176 }
0177 #endif