Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

simple_continued_fraction(std::vector<Z>) #975

Closed
wants to merge 3 commits into from
Closed

simple_continued_fraction(std::vector<Z>) #975

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d83bd53
Typedef for Z `int_type`; replaced member variable `x_` with `precisi…
Apr 11, 2023
25ab7c5
Implemented `simple_continued_fraction(std::vector<Z>)`
Apr 11, 2023
c4c4592
Implemented templated free function `Rational to_rational(simple_cont…
Apr 11, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
134 changes: 96 additions & 38 deletions include/boost/math/tools/simple_continued_fraction.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@
#include <sstream>
#include <utility>
#include <cstdint>
#include <cassert>

#include <boost/math/tools/is_standalone.hpp>
#ifndef BOOST_MATH_STANDALONE
Expand All @@ -35,34 +36,59 @@ namespace boost::math::tools {
template<typename Real, typename Z = int64_t>
class simple_continued_fraction {
public:
simple_continued_fraction(Real x) : x_{x} {
typedef Z int_type;

simple_continued_fraction(std::vector<Z> data) : b_{std::move(data)} {
const size_t size_ = b_.size();
if (size_ == 0) {
throw std::length_error("Array of coefficients is empty.");
}

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] << ".";
throw std::domain_error(oss.str());
}
}

canonicalize();
}

simple_continued_fraction(Real x) : b_{} {
using std::floor;
using std::abs;
using std::sqrt;
using std::isfinite;
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.");
}

if constexpr (std_precision == 2147483647) {
precision_ = x.backend().precision();
}

b_.reserve(50);
Real bj = floor(x);
b_.push_back(static_cast<Z>(bj));
if (bj == x) {
b_.shrink_to_fit();
return;
}

const Real orig_x = x;
x = 1/(x-bj);
Real f = bj;
if (bj == 0) {
f = 16*(std::numeric_limits<Real>::min)();
}
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<Real>::epsilon()*abs(x_))
{
const Real eps_abs_orig_x = std::numeric_limits<Real>::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<Z>(bj));
x = 1/(x-bj);
Expand All @@ -77,16 +103,10 @@ class simple_continued_fraction {
D = 1/D;
f *= (C*D);
}
// 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);
}
b_.shrink_to_fit();

for (size_t i = 1; i < b_.size(); ++i) {
canonicalize();

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] << "."
Expand All @@ -101,9 +121,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<Real>::quiet_NaN();
}
using std::log;
Expand All @@ -113,7 +134,7 @@ class simple_continued_fraction {
// A random partial denominator has ~80% chance of being in this table:
const std::array<Real, 7> logs{std::numeric_limits<Real>::quiet_NaN(), static_cast<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))};
Real log_prod = 0;
for (size_t i = 1; i < b_.size(); ++i) {
for (size_t i = 1; i < size_; ++i) {
if (b_[i] < static_cast<Z>(logs.size())) {
log_prod += logs[b_[i]];
}
Expand All @@ -122,49 +143,57 @@ class simple_continued_fraction {
log_prod += log(static_cast<Real>(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<Real>::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<Real>(b_[i]);
}
return n/denom;
}


// Note that this also includes the integer part (i.e. all the coefficients)
const std::vector<Z>& partial_denominators() const {
return b_;
}

inline std::vector<Z>&& get_data() noexcept
{
inline std::vector<Z>&& get_data() noexcept {
return std::move(b_);
}

template<typename T, typename Z2>
friend std::ostream& operator<<(std::ostream& out, simple_continued_fraction<T, Z2>& scf);

private:
const Real x_;
static constexpr int std_precision = std::numeric_limits<Real>::max_digits10;

void canonicalize() {
// 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_.pop_back();
b_.back() += 1;
}
b_.shrink_to_fit();
}

std::vector<Z> b_;

int precision_{std_precision};
};


template<typename Real, typename Z2>
std::ostream& operator<<(std::ostream& out, simple_continued_fraction<Real, Z2>& scf) {
constexpr const int p = std::numeric_limits<Real>::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)
{
Expand All @@ -182,9 +211,38 @@ std::ostream& operator<<(std::ostream& out, simple_continued_fraction<Real, Z2>&
template<typename Real, typename Z = std::int64_t>
inline auto simple_continued_fraction_coefficients(Real x)
{
auto temp = simple_continued_fraction(x);
auto temp = simple_continued_fraction<Real, Z>(x);
return temp.get_data();
}

// Can be used with `boost::rational` from <boost/rational.hpp>
template <typename Rational, typename Real, typename Z = std::int64_t>
inline Rational to_rational(const simple_continued_fraction<Real, Z>& scf)
{
using int_t = typename Rational::int_type;

auto& coefs = scf.partial_denominators();
const size_t size_ = coefs.size();
assert(size_ >= 1);
if (size_ == 1) return static_cast<int_t>(coefs[0]);

// p0 = a0, p1 = a1.a0 + 1, pn = an.pn-1 + pn-2 for 2 <= n
// q0 = 1, q1 = a1, qn = an.qn-1 + qn-2 for 2 <= n

int_t p0 = coefs[0];
int_t p1 = p0*coefs[1] + 1;
int_t q0 = 1;
int_t q1 = coefs[1];
for (size_t i = 2; i < size_; ++i) {
const Z cn = coefs[i];
const int_t pn = cn*p1 + p0;
const int_t qn = cn*q1 + q0;
p0 = std::exchange(p1, pn);
q0 = std::exchange(q1, qn);
}

return {p1, q1};
}

}
#endif