diff --git a/doc/equations/khinchin_geometric.svg b/doc/equations/khinchin_geometric.svg
new file mode 100644
index 000000000..c8c033943
--- /dev/null
+++ b/doc/equations/khinchin_geometric.svg
@@ -0,0 +1,2 @@
+
+
\ No newline at end of file
diff --git a/doc/equations/khinchin_harmonic.svg b/doc/equations/khinchin_harmonic.svg
new file mode 100644
index 000000000..be22f1fa1
--- /dev/null
+++ b/doc/equations/khinchin_harmonic.svg
@@ -0,0 +1,2 @@
+
+
\ No newline at end of file
diff --git a/doc/internals/simple_continued_fraction.qbk b/doc/internals/simple_continued_fraction.qbk
new file mode 100644
index 000000000..ae4b52a5d
--- /dev/null
+++ b/doc/internals/simple_continued_fraction.qbk
@@ -0,0 +1,62 @@
+[/
+ Copyright Nick Thompson, 2020
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or copy at
+ http://www.boost.org/LICENSE_1_0.txt).
+]
+
+[section:simple_continued_fraction Simple Continued Fractions]
+
+ #include
+ namespace boost::math::tools {
+
+ template
+ class simple_continued_fraction {
+ public:
+ simple_continued_fraction(Real x);
+
+ Real khinchin_geometric_mean() const;
+
+ Real khinchin_harmonic_mean() const;
+
+ template
+ friend std::ostream& operator<<(std::ostream& out, simple_continued_fraction& scf);
+ };
+ }
+
+
+The `simple_continued_fraction` class provided by Boost affords the ability to convert a floating point number into a simple continued fraction.
+In addition, we can answer a few questions about the number in question using this representation.
+
+Here's a minimal working example:
+
+ using boost::math::constants::pi;
+ using boost::math::tools::simple_continued_fraction;
+ auto cfrac = simple_continued_fraction(pi());
+ std::cout << "π ≈ " << cfrac << "\n";
+ // Prints:
+ // π ≈ 3.14159265358979323851 ≈ [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2]
+
+The class computes partial denominators which simultaneously computing convergents with the modified Lentz's algorithm.
+Once a convergent is within a few ulps of the input value, the computation stops.
+
+Note that every floating point number is a rational number, and this exact rational can be exactly converted to a finite continued fraction.
+This is perfectly sensible behavior, but we do not do it here.
+This is because when examining known values like π, it creates a large number of incorrect partial denominators, even if every bit of the binary representation is correct.
+
+It may be the case the a few incorrect partial convergents is harmless, but we compute continued fractions because we would like to do something with them.
+One sensible thing to do it to ask whether the number is in some sense "random"; a question that can be partially answered by computing the Khinchin geometric mean
+
+[$../equations/khinchin_geometric.svg]
+
+and Khinchin harmonic mean
+
+[$../equations/khinchin_harmonic.svg]
+
+If these approach Khinchin's constant /K/[sub 0] and /K/[sub -1] as the number of partial denominators goes to infinity, then our number is "uninteresting" with respect to the characterization.
+These violations are washed out if too many incorrect partial denominators are included in the expansion.
+
+Note: The convergence of these means to the Khinchin limit is exceedingly slow; we've used 30,000 decimal digits of π and only found two digits of agreement with /K/[sub 0].
+However, clear violations of are obvious, such as the continued fraction expansion of √2, whose Khinchin geometric mean is precisely 2.
+
+[endsect]
diff --git a/doc/math.qbk b/doc/math.qbk
index 6cedb58c8..09d8ad9b8 100644
--- a/doc/math.qbk
+++ b/doc/math.qbk
@@ -751,6 +751,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22.
[include internals/series.qbk]
[include internals/agm.qbk]
[include internals/fraction.qbk]
+[include internals/simple_continued_fraction.qbk]
[include internals/recurrence.qbk]
[/include internals/rational.qbk] [/moved to tools]
[include internals/tuple.qbk]
diff --git a/example/to_continued_fraction.cpp b/example/to_continued_fraction.cpp
new file mode 100644
index 000000000..7c7db7295
--- /dev/null
+++ b/example/to_continued_fraction.cpp
@@ -0,0 +1,36 @@
+// (C) Copyright Nick Thompson 2020.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#include
+#include
+#include
+#include
+
+using boost::math::constants::root_two;
+using boost::math::constants::phi;
+using boost::math::constants::pi;
+using boost::math::constants::e;
+using boost::math::constants::zeta_three;
+using boost::math::tools::simple_continued_fraction;
+
+int main()
+{
+ using Real = boost::multiprecision::cpp_bin_float_100;
+ auto phi_cfrac = simple_continued_fraction(phi());
+ std::cout << "φ ≈ " << phi_cfrac << "\n\n";
+
+ auto pi_cfrac = simple_continued_fraction(pi());
+ std::cout << "π ≈ " << pi_cfrac << "\n";
+ std::cout << "Known: [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, 1, 84, 2, 1, 1, 15, 3, 13, 1, 4, 2, 6, 6, 99, 1, 2, 2, 6, 3, 5, 1, 1, 6, 8, 1, 7, 1, 2, 3, 7, 1, 2, 1, 1, 12, 1, 1, 1, 3, 1, 1, 8, 1, 1, 2, 1, 6, 1, 1, 5, 2, 2, 3, 1, 2, 4, 4, 16, 1, 161, 45, 1, 22, 1, 2, 2, 1, 4, 1, 2, 24, 1, 2, 1, 3, 1, 2, 1, ...]\n\n";
+
+ auto rt_cfrac = simple_continued_fraction(root_two());
+ std::cout << "√2 ≈ " << rt_cfrac << "\n\n";
+
+ auto e_cfrac = simple_continued_fraction(e());
+ std::cout << "e ≈ " << e_cfrac << "\n";
+
+ // Correctness can be checked in Mathematica via: ContinuedFraction[Zeta[3], 500]
+ auto z_cfrac = simple_continued_fraction(zeta_three());
+ std::cout << "ζ(3) ≈ " << z_cfrac << "\n";
+}
diff --git a/include/boost/math/tools/simple_continued_fraction.hpp b/include/boost/math/tools/simple_continued_fraction.hpp
new file mode 100644
index 000000000..4663e5952
--- /dev/null
+++ b/include/boost/math/tools/simple_continued_fraction.hpp
@@ -0,0 +1,159 @@
+// (C) Copyright Nick Thompson 2020.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_TOOLS_SIMPLE_CONTINUED_FRACTION_HPP
+#define BOOST_MATH_TOOLS_SIMPLE_CONTINUED_FRACTION_HPP
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+namespace boost::math::tools {
+
+template
+class simple_continued_fraction {
+public:
+ simple_continued_fraction(Real x) : x_{x} {
+ 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.");
+ }
+ b_.reserve(50);
+ Real bj = floor(x);
+ b_.push_back(static_cast(bj));
+ if (bj == x) {
+ b_.shrink_to_fit();
+ return;
+ }
+ x = 1/(x-bj);
+ Real f = bj;
+ if (bj == 0) {
+ f = 16*std::numeric_limits::min();
+ }
+ Real C = f;
+ Real D = 0;
+ int i = 0;
+ // 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_))
+ {
+ bj = floor(x);
+ b_.push_back(static_cast(bj));
+ x = 1/(x-bj);
+ D += bj;
+ if (D == 0) {
+ D = 16*std::numeric_limits::min();
+ }
+ C = bj + 1/C;
+ if (C==0) {
+ C = 16*std::numeric_limits::min();
+ }
+ 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) {
+ if (b_[i] <= 0) {
+ std::ostringstream oss;
+ oss << "Found a negative partial denominator: b[" << i << "] = " << b_[i] << "."
+ << " This means the integer type '" << boost::core::demangle(typeid(Z).name())
+ << "' has overflowed and you need to use a wider type,"
+ << " or there is a bug.";
+ throw std::overflow_error(oss.str());
+ }
+ }
+ }
+
+ Real khinchin_geometric_mean() const {
+ if (b_.size() == 1) {
+ return std::numeric_limits::quiet_NaN();
+ }
+ using std::log;
+ using std::exp;
+ // Precompute the most probable logarithms. See the Gauss-Kuzmin distribution for details.
+ // Example: b_i = 1 has probability -log_2(3/4) ≈ .415:
+ // 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) {
+ if (b_[i] < static_cast(logs.size())) {
+ log_prod += logs[b_[i]];
+ }
+ else
+ {
+ log_prod += log(static_cast(b_[i]));
+ }
+ }
+ log_prod /= (b_.size()-1);
+ return exp(log_prod);
+ }
+
+ Real khinchin_harmonic_mean() const {
+ if (b_.size() == 1) {
+ return std::numeric_limits::quiet_NaN();
+ }
+ Real n = b_.size() - 1;
+ Real denom = 0;
+ for (size_t i = 1; i < b_.size(); ++i) {
+ denom += 1/static_cast(b_[i]);
+ }
+ return n/denom;
+ }
+
+ const std::vector& partial_denominators() const {
+ return b_;
+ }
+
+ template
+ friend std::ostream& operator<<(std::ostream& out, simple_continued_fraction& scf);
+
+private:
+ const Real x_;
+ 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.x_ << " ≈ [" << scf.b_.front();
+ if (scf.b_.size() > 1)
+ {
+ out << "; ";
+ for (size_t i = 1; i < scf.b_.size() -1; ++i)
+ {
+ out << scf.b_[i] << ", ";
+ }
+ out << scf.b_.back();
+ }
+ out << "]\n";
+ return out;
+}
+
+
+}
+#endif
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index aa56d14c7..03c2fc7c5 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -894,6 +894,7 @@ test-suite misc :
test_tr1_c_long_double
]
[ run test_constants.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ]
+ [ run simple_continued_fraction_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ]
[ run test_classify.cpp pch ../../test/build//boost_unit_test_framework ]
[ run test_error_handling.cpp ../../test/build//boost_unit_test_framework ]
[ run legendre_stieltjes_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_auto_declarations cxx11_range_based_for ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ]
diff --git a/test/math_unit_test.hpp b/test/math_unit_test.hpp
index e6d5e241e..718e061af 100644
--- a/test/math_unit_test.hpp
+++ b/test/math_unit_test.hpp
@@ -323,9 +323,15 @@ bool check_equal(Real x, Real y, std::string const & filename, std::string const
if (x != y) {
std::ios_base::fmtflags f( std::cerr.flags() );
std::cerr << "\033[0;31mError at " << filename << ":" << function << ":" << line << "\n";
- std::cerr << "\033[0m Expected " << x << " == " << y << " but:\n";
- std::cerr << " Expected = " << std::defaultfloat << std::fixed << x << " = " << std::scientific << x << std::hexfloat << " = " << x << "\n";
- std::cerr << " Computed = " << std::defaultfloat << std::fixed << y << " = " << std::scientific << y << std::hexfloat << " = " << y << "\n";
+ std::cerr << "\033[0m Condition '" << x << " == " << y << "' is not satisfied:\n";
+ if (std::is_floating_point::value) {
+ std::cerr << " Expected = " << std::defaultfloat << std::fixed << x << " = " << std::scientific << x << std::hexfloat << " = " << x << "\n";
+ std::cerr << " Computed = " << std::defaultfloat << std::fixed << y << " = " << std::scientific << y << std::hexfloat << " = " << y << "\n";
+ } else {
+ std::cerr << " Expected: " << x << " = " << "0x" << std::hex << x << "\n";
+ std::cerr << std::dec;
+ std::cerr << " Computed: " << y << " = " << "0x" << std::hex << y << "\n";
+ }
std::cerr.flags(f);
++detail::global_error_count;
return false;
diff --git a/test/simple_continued_fraction_test.cpp b/test/simple_continued_fraction_test.cpp
new file mode 100644
index 000000000..637916496
--- /dev/null
+++ b/test/simple_continued_fraction_test.cpp
@@ -0,0 +1,160 @@
+/*
+ * Copyright Nick Thompson, 2020
+ * Use, modification and distribution are subject to the
+ * Boost Software License, Version 1.0. (See accompanying file
+ * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#include "math_unit_test.hpp"
+#include
+#include
+#ifdef BOOST_HAS_FLOAT128
+#include
+using boost::multiprecision::float128;
+#endif
+#include
+
+using boost::math::tools::simple_continued_fraction;
+using boost::multiprecision::cpp_bin_float_100;
+using boost::math::constants::pi;
+
+template
+void test_integral()
+{
+ for (int64_t i = -20; i < 20; ++i) {
+ Real ii = i;
+ auto cfrac = simple_continued_fraction(ii);
+ auto const & a = cfrac.partial_denominators();
+ CHECK_EQUAL(size_t(1), a.size());
+ CHECK_EQUAL(i, a.front());
+ }
+}
+
+template
+void test_halves()
+{
+ for (int64_t i = -20; i < 20; ++i) {
+ Real x = i + Real(1)/Real(2);
+ auto cfrac = simple_continued_fraction(x);
+ auto const & a = cfrac.partial_denominators();
+ CHECK_EQUAL(size_t(2), a.size());
+ CHECK_EQUAL(i, a.front());
+ CHECK_EQUAL(int64_t(2), a.back());
+ }
+
+ // We'll also test quarters; why not?
+ for (int64_t i = -20; i < 20; ++i) {
+ Real x = i + Real(1)/Real(4);
+ auto cfrac = simple_continued_fraction(x);
+ auto const & a = cfrac.partial_denominators();
+ CHECK_EQUAL(size_t(2), a.size());
+ CHECK_EQUAL(i, a.front());
+ CHECK_EQUAL(int64_t(4), a.back());
+ }
+
+ for (int64_t i = -20; i < 20; ++i) {
+ Real x = i + Real(1)/Real(8);
+ auto cfrac = simple_continued_fraction(x);
+ auto const & a = cfrac.partial_denominators();
+ CHECK_EQUAL(size_t(2), a.size());
+ CHECK_EQUAL(i, a.front());
+ CHECK_EQUAL(int64_t(8), a.back());
+ }
+
+ for (int64_t i = -20; i < 20; ++i) {
+ Real x = i + Real(3)/Real(4);
+ auto cfrac = simple_continued_fraction(x);
+ auto const & a = cfrac.partial_denominators();
+ CHECK_EQUAL(size_t(3), a.size());
+ CHECK_EQUAL(i, a.front());
+ CHECK_EQUAL(int64_t(1), a[1]);
+ CHECK_EQUAL(int64_t(3), a.back());
+ }
+
+ for (int64_t i = -20; i < 20; ++i) {
+ Real x = i + Real(7)/Real(8);
+ auto cfrac = simple_continued_fraction(x);
+ auto const & a = cfrac.partial_denominators();
+ CHECK_EQUAL(size_t(3), a.size());
+ CHECK_EQUAL(i, a.front());
+ CHECK_EQUAL(int64_t(1), a[1]);
+ CHECK_EQUAL(int64_t(7), a.back());
+ }
+}
+
+template
+void test_simple()
+{
+ std::cout << "Testing rational numbers on type " << boost::core::demangle(typeid(Real).name()) << "\n";
+ {
+ Real x = Real(649)/200;
+ // ContinuedFraction[649/200] = [3; 4, 12, 4]
+ auto cfrac = simple_continued_fraction(x);
+ auto const & a = cfrac.partial_denominators();
+ CHECK_EQUAL(size_t(4), a.size());
+ CHECK_EQUAL(int64_t(3), a[0]);
+ CHECK_EQUAL(int64_t(4), a[1]);
+ CHECK_EQUAL(int64_t(12), a[2]);
+ CHECK_EQUAL(int64_t(4), a[3]);
+ }
+
+ {
+ Real x = Real(415)/Real(93);
+ // [4; 2, 6, 7]:
+ auto cfrac = simple_continued_fraction(x);
+ auto const & a = cfrac.partial_denominators();
+ CHECK_EQUAL(size_t(4), a.size());
+ CHECK_EQUAL(int64_t(4), a[0]);
+ CHECK_EQUAL(int64_t(2), a[1]);
+ CHECK_EQUAL(int64_t(6), a[2]);
+ CHECK_EQUAL(int64_t(7), a[3]);
+ }
+
+}
+
+template
+void test_khinchin()
+{
+ // These are simply sanity checks; the convergence is too slow otherwise:
+ auto cfrac = simple_continued_fraction(pi());
+ auto K0 = cfrac.khinchin_geometric_mean();
+ CHECK_MOLLIFIED_CLOSE(Real(2.6854520010), K0, 0.1);
+ auto Km1 = cfrac.khinchin_harmonic_mean();
+ CHECK_MOLLIFIED_CLOSE(Real(1.74540566240), Km1, 0.1);
+
+ using std::sqrt;
+ auto rt_cfrac = simple_continued_fraction(sqrt(static_cast(2)));
+ K0 = rt_cfrac.khinchin_geometric_mean();
+ CHECK_ULP_CLOSE(Real(2), K0, 10);
+ Km1 = rt_cfrac.khinchin_harmonic_mean();
+ CHECK_ULP_CLOSE(Real(2), Km1, 10);
+}
+
+
+int main()
+{
+ test_integral();
+ test_integral();
+ test_integral();
+ test_integral();
+
+ test_halves();
+ test_halves();
+ test_halves();
+ test_halves();
+
+ test_simple();
+ test_simple();
+ test_simple();
+ test_simple();
+
+ test_khinchin();
+
+ #ifdef BOOST_HAS_FLOAT128
+ test_integral();
+ test_halves();
+ test_simple();
+ test_khinchin();
+ #endif
+ return boost::math::test::report_errors();
+}