174 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_CENTERED_CONTINUED_FRACTION_HPP7#define BOOST_MATH_TOOLS_CENTERED_CONTINUED_FRACTION_HPP8 9#include <cmath>10#include <cstdint>11#include <vector>12#include <ostream>13#include <iomanip>14#include <limits>15#include <stdexcept>16#include <sstream>17#include <array>18#include <type_traits>19#include <boost/math/tools/is_standalone.hpp>20 21#ifndef BOOST_MATH_STANDALONE22#include <boost/config.hpp>23#ifdef BOOST_MATH_NO_CXX17_IF_CONSTEXPR24#error "The header <boost/math/norms.hpp> can only be used in C++17 and later."25#endif26#endif27 28#ifndef BOOST_MATH_STANDALONE29#include <boost/core/demangle.hpp>30#endif31 32namespace boost::math::tools {33 34template<typename Real, typename Z = int64_t>35class centered_continued_fraction {36public:37 centered_continued_fraction(Real x) : x_{x} {38 static_assert(std::is_integral_v<Z> && std::is_signed_v<Z>,39 "Centered continued fractions require signed integer types.");40 using std::round;41 using std::abs;42 using std::sqrt;43 using std::isfinite;44 if (!isfinite(x))45 {46 throw std::domain_error("Cannot convert non-finites into continued fractions."); 47 }48 b_.reserve(50);49 Real bj = round(x);50 b_.push_back(static_cast<Z>(bj));51 if (bj == x)52 {53 b_.shrink_to_fit();54 return;55 }56 x = 1/(x-bj);57 Real f = bj;58 if (bj == 0)59 {60 f = 16*(std::numeric_limits<Real>::min)();61 }62 Real C = f;63 Real D = 0;64 int i = 0;65 while (abs(f - x_) >= (1 + i++)*std::numeric_limits<Real>::epsilon()*abs(x_))66 {67 bj = round(x);68 b_.push_back(static_cast<Z>(bj));69 x = 1/(x-bj);70 D += bj;71 if (D == 0) {72 D = 16*(std::numeric_limits<Real>::min)();73 }74 C = bj + 1/C;75 if (C==0)76 {77 C = 16*(std::numeric_limits<Real>::min)();78 }79 D = 1/D;80 f *= (C*D);81 }82 // Deal with non-uniqueness of continued fractions: [a0; a1, ..., an, 1] = a0; a1, ..., an + 1].83 if (b_.size() > 2 && b_.back() == 1)84 {85 b_[b_.size() - 2] += 1;86 b_.resize(b_.size() - 1);87 }88 b_.shrink_to_fit();89 90 for (size_t i = 1; i < b_.size(); ++i)91 {92 if (b_[i] == 0) {93 std::ostringstream oss;94 oss << "Found a zero partial denominator: b[" << i << "] = " << b_[i] << "."95 #ifndef BOOST_MATH_STANDALONE96 << " This means the integer type '" << boost::core::demangle(typeid(Z).name())97 #else98 << " This means the integer type '" << typeid(Z).name()99 #endif100 << "' has overflowed and you need to use a wider type,"101 << " or there is a bug.";102 throw std::overflow_error(oss.str());103 }104 }105 }106 107 Real khinchin_geometric_mean() const {108 if (b_.size() == 1)109 { 110 return std::numeric_limits<Real>::quiet_NaN();111 }112 using std::log;113 using std::exp;114 using std::abs;115 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))};116 Real log_prod = 0;117 for (size_t i = 1; i < b_.size(); ++i)118 {119 if (abs(b_[i]) < static_cast<Z>(logs.size()))120 {121 log_prod += logs[abs(b_[i])];122 }123 else124 {125 log_prod += log(static_cast<Real>(abs(b_[i])));126 }127 }128 log_prod /= (b_.size()-1);129 return exp(log_prod);130 }131 132 const std::vector<Z>& partial_denominators() const {133 return b_;134 }135 136 template<typename T, typename Z2>137 friend std::ostream& operator<<(std::ostream& out, centered_continued_fraction<T, Z2>& ccf);138 139private:140 const Real x_;141 std::vector<Z> b_;142};143 144 145template<typename Real, typename Z2>146std::ostream& operator<<(std::ostream& out, centered_continued_fraction<Real, Z2>& scf) {147 constexpr const int p = std::numeric_limits<Real>::max_digits10;148 if constexpr (p == 2147483647)149 {150 out << std::setprecision(scf.x_.backend().precision());151 }152 else153 {154 out << std::setprecision(p);155 }156 157 out << "[" << scf.b_.front();158 if (scf.b_.size() > 1)159 {160 out << "; ";161 for (size_t i = 1; i < scf.b_.size() -1; ++i)162 {163 out << scf.b_[i] << ", ";164 }165 out << scf.b_.back();166 }167 out << "]";168 return out;169}170 171 172}173#endif174