math/tools/centered_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_CENTERED_CONTINUED_FRACTION_HPP
0007 #define BOOST_MATH_TOOLS_CENTERED_CONTINUED_FRACTION_HPP
0008
0009 #include <cmath>
0010 #include <cstdint>
0011 #include <vector>
0012 #include <ostream>
0013 #include <iomanip>
0014 #include <limits>
0015 #include <stdexcept>
0016 #include <sstream>
0017 #include <array>
0018 #include <type_traits>
0019 #include <boost/math/tools/is_standalone.hpp>
0020
0021 #ifndef BOOST_MATH_STANDALONE
0022 #include <boost/config.hpp>
0023 #ifdef BOOST_NO_CXX17_IF_CONSTEXPR
0024 #error "The header <boost/math/norms.hpp> can only be used in C++17 and later."
0025 #endif
0026 #endif
0027
0028 #ifndef BOOST_MATH_STANDALONE
0029 #include <boost/core/demangle.hpp>
0030 #endif
0031
0032 namespace boost::math::tools {
0033
0034 template<typename Real, typename Z = int64_t>
0035 class centered_continued_fraction {
0036 public:
0037     centered_continued_fraction(Real x) : x_{x} {
0038         static_assert(std::is_integral_v<Z> && std::is_signed_v<Z>,
0039                       "Centered continued fractions require signed integer types.");
0040         using std::round;
0041         using std::abs;
0042         using std::sqrt;
0043         using std::isfinite;
0044         if (!isfinite(x))
0045         {
0046             throw std::domain_error("Cannot convert non-finites into continued fractions.");
0047         }
0048         b_.reserve(50);
0049         Real bj = round(x);
0050         b_.push_back(static_cast<Z>(bj));
0051         if (bj == x)
0052         {
0053             b_.shrink_to_fit();
0054             return;
0055         }
0056         x = 1/(x-bj);
0057         Real f = bj;
0058         if (bj == 0)
0059         {
0060             f = 16*(std::numeric_limits<Real>::min)();
0061         }
0062         Real C = f;
0063         Real D = 0;
0064         int i = 0;
0065         while (abs(f - x_) >= (1 + i++)*std::numeric_limits<Real>::epsilon()*abs(x_))
0066         {
0067             bj = round(x);
0068             b_.push_back(static_cast<Z>(bj));
0069             x = 1/(x-bj);
0070             D += bj;
0071             if (D == 0) {
0072                 D = 16*(std::numeric_limits<Real>::min)();
0073             }
0074             C = bj + 1/C;
0075             if (C==0)
0076             {
0077                 C = 16*(std::numeric_limits<Real>::min)();
0078             }
0079             D = 1/D;
0080             f *= (C*D);
0081         }
0082         // Deal with non-uniqueness of continued fractions: [a0; a1, ..., an, 1] = a0; a1, ..., an + 1].
0083         if (b_.size() > 2 && b_.back() == 1)
0084         {
0085             b_[b_.size() - 2] += 1;
0086             b_.resize(b_.size() - 1);
0087         }
0088         b_.shrink_to_fit();
0089
0090         for (size_t i = 1; i < b_.size(); ++i)
0091         {
0092             if (b_[i] == 0) {
0093                 std::ostringstream oss;
0094                 oss << "Found a zero partial denominator: b[" << i << "] = " << b_[i] << "."
0095                     #ifndef BOOST_MATH_STANDALONE
0096                     << " This means the integer type '" << boost::core::demangle(typeid(Z).name())
0097                     #else
0098                     << " This means the integer type '" << typeid(Z).name()
0099                     #endif
0100                     << "' has overflowed and you need to use a wider type,"
0101                     << " or there is a bug.";
0102                 throw std::overflow_error(oss.str());
0103             }
0104         }
0105     }
0106
0107     Real khinchin_geometric_mean() const {
0108         if (b_.size() == 1)
0109         {
0110             return std::numeric_limits<Real>::quiet_NaN();
0111         }
0112         using std::log;
0113         using std::exp;
0114         using std::abs;
0115         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))};
0116         Real log_prod = 0;
0117         for (size_t i = 1; i < b_.size(); ++i)
0118         {
0119             if (abs(b_[i]) < static_cast<Z>(logs.size()))
0120             {
0121                 log_prod += logs[abs(b_[i])];
0122             }
0123             else
0124             {
0125                 log_prod += log(static_cast<Real>(abs(b_[i])));
0126             }
0127         }
0128         log_prod /= (b_.size()-1);
0129         return exp(log_prod);
0130     }
0131
0132     const std::vector<Z>& partial_denominators() const {
0133         return b_;
0134     }
0135
0136     template<typename T, typename Z2>
0137     friend std::ostream& operator<<(std::ostream& out, centered_continued_fraction<T, Z2>& ccf);
0138
0139 private:
0140     const Real x_;
0141     std::vector<Z> b_;
0142 };
0143
0144
0145 template<typename Real, typename Z2>
0146 std::ostream& operator<<(std::ostream& out, centered_continued_fraction<Real, Z2>& scf) {
0147     constexpr const int p = std::numeric_limits<Real>::max_digits10;
0148     if constexpr (p == 2147483647)
0149     {
0150         out << std::setprecision(scf.x_.backend().precision());
0151     }
0152     else
0153     {
0154         out << std::setprecision(p);
0155     }
0156
0157     out << "[" << scf.b_.front();
0158     if (scf.b_.size() > 1)
0159     {
0160         out << "; ";
0161         for (size_t i = 1; i < scf.b_.size() -1; ++i)
0162         {
0163             out << scf.b_[i] << ", ";
0164         }
0165         out << scf.b_.back();
0166     }
0167     out << "]";
0168     return out;
0169 }
0170
0171
0172 }
0173 #endif