From 1ac89b2b02169d4c19c04771df01ff4ed77bb7c1 Mon Sep 17 00:00:00 2001 From: Nick Date: Fri, 26 Jun 2020 14:50:04 -0400 Subject: [PATCH] Simple continued fraction (#377) * Simple continued fraction [CI SKIP] * Comments on error analysis [CI SKIP] * Simple continued fraction [CI SKIP] * Clarify comment and kick off build. --- doc/equations/khinchin_geometric.svg | 2 + doc/equations/khinchin_harmonic.svg | 2 + doc/internals/simple_continued_fraction.qbk | 62 +++++++ doc/math.qbk | 1 + example/to_continued_fraction.cpp | 36 ++++ .../math/tools/simple_continued_fraction.hpp | 159 +++++++++++++++++ test/Jamfile.v2 | 1 + test/math_unit_test.hpp | 12 +- test/simple_continued_fraction_test.cpp | 160 ++++++++++++++++++ 9 files changed, 432 insertions(+), 3 deletions(-) create mode 100644 doc/equations/khinchin_geometric.svg create mode 100644 doc/equations/khinchin_harmonic.svg create mode 100644 doc/internals/simple_continued_fraction.qbk create mode 100644 example/to_continued_fraction.cpp create mode 100644 include/boost/math/tools/simple_continued_fraction.hpp create mode 100644 test/simple_continued_fraction_test.cpp 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(); +}