diff --git a/include/boost/math/tools/simple_continued_fraction.hpp b/include/boost/math/tools/simple_continued_fraction.hpp index 0fe17fc59d..d4679f04aa 100644 --- a/include/boost/math/tools/simple_continued_fraction.hpp +++ b/include/boost/math/tools/simple_continued_fraction.hpp @@ -32,17 +32,26 @@ namespace boost::math::tools { template class simple_continued_fraction { public: - simple_continued_fraction(Real x) : x_{x} { + using int_type = Z; + + simple_continued_fraction() = default; + ~simple_continued_fraction() = default; + simple_continued_fraction(const simple_continued_fraction&) = default; + simple_continued_fraction& operator =(const simple_continued_fraction&) = default; + simple_continued_fraction(simple_continued_fraction&&) = default; + simple_continued_fraction& operator =(simple_continued_fraction&&) = default; + 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."); } - b_.reserve(50); + reserve(50); Real bj = floor(x); - b_.push_back(static_cast(bj)); + push_back(bj); if (bj == x) { b_.shrink_to_fit(); return; @@ -54,14 +63,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)); + push_back(bj); x = 1/(x-bj); D += bj; if (D == 0) { @@ -77,13 +85,14 @@ class simple_continued_fraction { // Deal with non-uniqueness of continued fractions: [a0; a1, ..., an, 1] = a0; a1, ..., an + 1]. // 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); + if (size() > 2 && back() == 1) { + pop_back(); + back() += 1; } b_.shrink_to_fit(); - - for (size_t i = 1; i < b_.size(); ++i) { + + const int size_ = size(); + for (int i = 1; i < size_; ++i) { if (b_[i] <= 0) { std::ostringstream oss; oss << "Found a negative partial denominator: b[" << i << "] = " << b_[i] << "." @@ -98,9 +107,10 @@ class simple_continued_fraction { } } } - + Real khinchin_geometric_mean() const { - if (b_.size() == 1) { + const int size_ = size(); + if (size_ == 1) { return std::numeric_limits::quiet_NaN(); } using std::log; @@ -110,7 +120,7 @@ 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) { + for (int i = 1; i < size_; ++i) { if (b_[i] < static_cast(logs.size())) { log_prod += logs[b_[i]]; } @@ -119,53 +129,90 @@ 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 int size_ = 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 (int i = 1; i < size_; ++i) { denom += 1/static_cast(b_[i]); } return n/denom; } - - const std::vector& partial_denominators() const { - return b_; + + size_t size() const noexcept { + return b_.size(); + } + int_type operator [](int idx) const { + return b_[idx]; + } + int_type& operator [](int idx) { + return b_[idx]; + } + int_type front() const { + return b_.front(); + } + int_type& front() { + return b_.front(); + } + int_type back() const { + return b_.back(); + } + int_type& back() { + return b_.back(); } - + auto begin() const noexcept { + return b_.begin(); + } + auto begin() noexcept { + return b_.begin(); + } + auto end() const noexcept { + return b_.end(); + } + auto end() noexcept { + return b_.end(); + } + + void reserve(size_t size_) { + b_.reserve(size_); + } + void push_back(int_type c) { + b_.push_back(c); + } + void pop_back() { + b_.pop_back(); + } + template friend std::ostream& operator<<(std::ostream& out, simple_continued_fraction& scf); private: - const Real x_; - std::vector b_; + std::vector b_{}; }; 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 << "[" << scf.b_.front(); - if (scf.b_.size() > 1) + static_assert(p != 2147483647, "numeric_limits::max_digits10 == 2147483647"); + out << std::setprecision(p); + + out << "[" << scf.front(); + if (scf.size() > 1) { out << "; "; - for (size_t i = 1; i < scf.b_.size() -1; ++i) + for (size_t i = 1; i < scf.size() -1; ++i) { out << scf.b_[i] << ", "; } - out << scf.b_.back(); + out << scf.back(); } out << "]"; return out;