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