From 010105113de72a9f4a76cee806feef56f093154b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Kol=C3=A1rik?= Date: Fri, 31 Mar 2023 15:34:28 +0200 Subject: [PATCH] simple_continued_fraction: added protected non-const accessor to partial denominators (#970) --- .../math/tools/simple_continued_fraction.hpp | 85 +++++++++++-------- 1 file changed, 50 insertions(+), 35 deletions(-) diff --git a/include/boost/math/tools/simple_continued_fraction.hpp b/include/boost/math/tools/simple_continued_fraction.hpp index 0fe17fc59d..94d2f3155a 100644 --- a/include/boost/math/tools/simple_continued_fraction.hpp +++ b/include/boost/math/tools/simple_continued_fraction.hpp @@ -32,17 +32,28 @@ namespace boost::math::tools { template class simple_continued_fraction { public: - simple_continued_fraction(Real x) : x_{x} { + typedef Z int_type; + + simple_continued_fraction(Real x) { using std::floor; using std::abs; using std::sqrt; using std::isfinite; + const Real orig_x = x; if (!isfinite(x)) { - throw std::domain_error("Cannot convert non-finites into continued fractions."); + throw std::domain_error("Cannot convert non-finites into continued fractions."); + } + + constexpr int p = std::numeric_limits::max_digits10; + if constexpr (p == 2147483647) { + precision_ = orig_x.backend().precision(); + } else { + precision_ = p; } + b_.reserve(50); Real bj = floor(x); - b_.push_back(static_cast(bj)); + b_.push_back(static_cast(bj)); if (bj == x) { b_.shrink_to_fit(); return; @@ -54,14 +65,13 @@ class simple_continued_fraction { } Real C = f; Real D = 0; - int i = 0; - // the "1 + i++" lets the error bound grow slowly with the number of convergents. + // the "1 + i" lets the error bound grow slowly with the number of convergents. // I have not worked out the error propagation of the Modified Lentz's method to see if it does indeed grow at this rate. // Numerical Recipes claims that no one has worked out the error analysis of the modified Lentz's method. - while (abs(f - x_) >= (1 + i++)*std::numeric_limits::epsilon()*abs(x_)) - { + const Real eps_abs_orig_x = std::numeric_limits::epsilon()*abs(orig_x); + for (int i = 0; abs(f - orig_x) >= (1 + i)*eps_abs_orig_x; ++i) { bj = floor(x); - b_.push_back(static_cast(bj)); + b_.push_back(static_cast(bj)); x = 1/(x-bj); D += bj; if (D == 0) { @@ -78,19 +88,20 @@ class simple_continued_fraction { // The shorter representation is considered the canonical representation, // so if we compute a non-canonical representation, change it to canonical: if (b_.size() > 2 && b_.back() == 1) { - b_[b_.size() - 2] += 1; - b_.resize(b_.size() - 1); + b_.pop_back(); + b_.back() += 1; } b_.shrink_to_fit(); - - for (size_t i = 1; i < b_.size(); ++i) { + + const size_t size_ = b_.size(); + for (size_t i = 1; i < size_; ++i) { if (b_[i] <= 0) { std::ostringstream oss; oss << "Found a negative partial denominator: b[" << i << "] = " << b_[i] << "." #ifndef BOOST_MATH_STANDALONE - << " This means the integer type '" << boost::core::demangle(typeid(Z).name()) + << " This means the integer type '" << boost::core::demangle(typeid(int_type).name()) #else - << " This means the integer type '" << typeid(Z).name() + << " This means the integer type '" << typeid(int_type).name() #endif << "' has overflowed and you need to use a wider type," << " or there is a bug."; @@ -98,9 +109,10 @@ class simple_continued_fraction { } } } - + Real khinchin_geometric_mean() const { - if (b_.size() == 1) { + const size_t size_ = b_.size(); + if (size_ == 1) { return std::numeric_limits::quiet_NaN(); } using std::log; @@ -110,8 +122,8 @@ class simple_continued_fraction { // A random partial denominator has ~80% chance of being in this table: const std::array logs{std::numeric_limits::quiet_NaN(), Real(0), log(static_cast(2)), log(static_cast(3)), log(static_cast(4)), log(static_cast(5)), log(static_cast(6))}; Real log_prod = 0; - for (size_t i = 1; i < b_.size(); ++i) { - if (b_[i] < static_cast(logs.size())) { + for (size_t i = 1; i < size_; ++i) { + if (b_[i] < static_cast(logs.size())) { log_prod += logs[b_[i]]; } else @@ -119,44 +131,47 @@ class simple_continued_fraction { log_prod += log(static_cast(b_[i])); } } - log_prod /= (b_.size()-1); + log_prod /= (size_-1); return exp(log_prod); } - + Real khinchin_harmonic_mean() const { - if (b_.size() == 1) { + const size_t size_ = b_.size(); + if (size_ == 1) { return std::numeric_limits::quiet_NaN(); } - Real n = b_.size() - 1; + Real n = size_ - 1; Real denom = 0; - for (size_t i = 1; i < b_.size(); ++i) { + for (size_t i = 1; i < size_; ++i) { denom += 1/static_cast(b_[i]); } return n/denom; } - - const std::vector& partial_denominators() const { + + // Note that this also includes the integer part (i.e. all the coefficients) + const std::vector& partial_denominators() const { return b_; } - + template friend std::ostream& operator<<(std::ostream& out, simple_continued_fraction& scf); +protected: + simple_continued_fraction() = default; + // Note that non-const operations may result in invalid simple continued fraction + std::vector& partial_denominators() { + return b_; + } private: - const Real x_; - std::vector b_; + std::vector b_{}; + + int precision_{}; }; template std::ostream& operator<<(std::ostream& out, simple_continued_fraction& scf) { - constexpr const int p = std::numeric_limits::max_digits10; - if constexpr (p == 2147483647) { - out << std::setprecision(scf.x_.backend().precision()); - } else { - out << std::setprecision(p); - } - + out << std::setprecision(scf.precision_); out << "[" << scf.b_.front(); if (scf.b_.size() > 1) {