From ee2cd5d5e59733004fef3b7bf49ca9e29ed384cb Mon Sep 17 00:00:00 2001 From: Nick Date: Sat, 18 Jul 2020 09:28:39 -0400 Subject: [PATCH] Luroth expansions (#401) --- doc/equations/luroth_expansion.svg | 2 + doc/internals/luroth_expansion.qbk | 62 ++++++++ doc/math.qbk | 1 + example/luroth.cpp | 18 +++ include/boost/math/tools/luroth_expansion.hpp | 138 ++++++++++++++++++ test/Jamfile.v2 | 1 + test/luroth_expansion_test.cpp | 95 ++++++++++++ 7 files changed, 317 insertions(+) create mode 100644 doc/equations/luroth_expansion.svg create mode 100644 doc/internals/luroth_expansion.qbk create mode 100644 example/luroth.cpp create mode 100644 include/boost/math/tools/luroth_expansion.hpp create mode 100644 test/luroth_expansion_test.cpp diff --git a/doc/equations/luroth_expansion.svg b/doc/equations/luroth_expansion.svg new file mode 100644 index 000000000..1f7b4af5c --- /dev/null +++ b/doc/equations/luroth_expansion.svg @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/doc/internals/luroth_expansion.qbk b/doc/internals/luroth_expansion.qbk new file mode 100644 index 000000000..70002dac4 --- /dev/null +++ b/doc/internals/luroth_expansion.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:luroth_expansion Luroth Expansions] + + #include + namespace boost::math::tools { + + template + class luroth_expansion { + public: + luroth_expansion(Real x); + + std::vector const & digits() const; + + Real digit_geometric_mean() const; + + template + friend std::ostream& operator<<(std::ostream& out, luroth_expansion& luroth); + }; + } + + +The `luroth_expansion` class provided by Boost expands a floating point number into a Lüroth representation, i.e., + +[$../equations/luroth_expansion.svg] + +The numbers /d/[sub i] are called digits or denominators; we use the terminology digits, since technically in our notation /d/[sub 0] is not a denominator. + +Here's a minimal working example: + + using boost::math::constants::pi; + using boost::math::tools::luroth_expansion; + auto luroth = luroth_expansion(pi()); + std::cout << "π ≈ " << luroth << "\n"; + // Prints: + // π ≈ ((3; 7, 1, 1, 1, 2, 1, 4, 23, 4, 1, 1, 1, 1, 80, 1, 1, 5)) + +The class computes denominators while simultaneously computing convergents. +Once a convergent is within a few ulps of the input value, the computation stops. + +/Nota bene:/ There is an alternative definition of the Lüroth representation where every digit is shifted by 1. +We follow the definition given in Kalpazidou; with the modification that we do not constrain the input to be in the interval [0,1] +and let the first digit be the floor of the input. + +For almost all real numbers, the geometric mean of the digits converges to a constant which is approximately 2.2001610580. +This is "Khinchin's constant" for the Lüroth representation. + + + +[heading References] + +* Kalpazidou, Sofia. "Khintchine's constant for Lüroth representation." Journal of Number Theory 29.2 (1988): 196-205. + +* Finch, Steven R. Mathematical constants. Cambridge university press, 2003. + + +[endsect] diff --git a/doc/math.qbk b/doc/math.qbk index df0770f6f..dd97e7228 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -753,6 +753,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22. [include internals/fraction.qbk] [include internals/simple_continued_fraction.qbk] [include internals/centered_continued_fraction.qbk] +[include internals/luroth_expansion.qbk] [include internals/recurrence.qbk] [/include internals/rational.qbk] [/moved to tools] [include internals/tuple.qbk] diff --git a/example/luroth.cpp b/example/luroth.cpp new file mode 100644 index 000000000..cfaa55387 --- /dev/null +++ b/example/luroth.cpp @@ -0,0 +1,18 @@ +// (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 + +using boost::math::constants::pi; +using boost::math::tools::luroth_expansion; +using boost::multiprecision::mpfr_float; + +int main() { + using Real = mpfr_float; + mpfr_float::default_precision(1024); + auto luroth = luroth_expansion(pi()); + std::cout << luroth << "\n"; +} diff --git a/include/boost/math/tools/luroth_expansion.hpp b/include/boost/math/tools/luroth_expansion.hpp new file mode 100644 index 000000000..385a02e04 --- /dev/null +++ b/include/boost/math/tools/luroth_expansion.hpp @@ -0,0 +1,138 @@ +// (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_LUROTH_EXPANSION_HPP +#define BOOST_MATH_TOOLS_LUROTH_EXPANSION_HPP + +#include +#include +#include +#include +#include +#include + +namespace boost::math::tools { + +template +class luroth_expansion { +public: + luroth_expansion(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 a Lüroth representation."); + } + d_.reserve(50); + Real dn = floor(x); + d_.push_back(static_cast(dn)); + if (dn == x) + { + d_.shrink_to_fit(); + return; + } + // This attempts to follow the notation of: + // "Khinchine's constant for Luroth Representation", by Sophia Kalpazidou. + x = x - dn; + Real computed = dn; + Real prod = 1; + // Let the error bound grow by 1 ULP/iteration. + // I haven't done the error analysis to show that this is an expected rate of error growth, + // but if you don't do this, you can easily get into an infinite loop. + Real i = 1; + Real scale = std::numeric_limits::epsilon()*abs(x_)/2; + while (abs(x_ - computed) > (i++)*scale) + { + Real recip = 1/x; + Real dn = floor(recip); + // x = n + 1/k => lur(x) = ((n; k - 1)) + // Note that this is a bit different than Kalpazidou (examine the half-open interval of definition carefully). + // One way to examine this definition is better for rationals (it never happens for irrationals) + // is to consider i + 1/3. If you follow Kalpazidou, then you get ((i, 3, 0)); a zero digit! + // That's bad since it destroys uniqueness and also breaks the computation of the geometric mean. + if (recip == dn) { + d_.push_back(static_cast(dn - 1)); + break; + } + d_.push_back(static_cast(dn)); + Real tmp = 1/(dn+1); + computed += prod*tmp; + prod *= tmp/dn; + x = dn*(dn+1)*(x - tmp); + } + + for (size_t i = 1; i < d_.size(); ++i) + { + // Sanity check: + if (d_[i] <= 0) + { + throw std::domain_error("Found a digit <= 0; this is an error."); + } + } + d_.shrink_to_fit(); + } + + + const std::vector& digits() const { + return d_; + } + + // Under the assumption of 'randomness', this mean converges to 2.2001610580. + // See Finch, Mathematical Constants, section 1.8.1. + Real digit_geometric_mean() const { + if (d_.size() == 1) { + return std::numeric_limits::quiet_NaN(); + } + using std::log; + using std::exp; + Real g = 0; + for (size_t i = 1; i < d_.size(); ++i) { + g += log(static_cast(d_[i])); + } + return exp(g/(d_.size() - 1)); + } + + template + friend std::ostream& operator<<(std::ostream& out, luroth_expansion& scf); + +private: + const Real x_; + std::vector d_; +}; + + +template +std::ostream& operator<<(std::ostream& out, luroth_expansion& luroth) +{ + constexpr const int p = std::numeric_limits::max_digits10; + if constexpr (p == 2147483647) + { + out << std::setprecision(luroth.x_.backend().precision()); + } + else + { + out << std::setprecision(p); + } + + out << "((" << luroth.d_.front(); + if (luroth.d_.size() > 1) + { + out << "; "; + for (size_t i = 1; i < luroth.d_.size() -1; ++i) + { + out << luroth.d_[i] << ", "; + } + out << luroth.d_.back(); + } + out << "))"; + return out; +} + + +} +#endif diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index c629dc290..a98f8dc31 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -896,6 +896,7 @@ test-suite misc : [ 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 centered_continued_fraction_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] + [ run luroth_expansion_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/luroth_expansion_test.cpp b/test/luroth_expansion_test.cpp new file mode 100644 index 000000000..f9664f108 --- /dev/null +++ b/test/luroth_expansion_test.cpp @@ -0,0 +1,95 @@ +/* + * 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::luroth_expansion; +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 luroth = luroth_expansion(ii); + auto const & a = luroth.digits(); + CHECK_EQUAL(size_t(1), a.size()); + CHECK_EQUAL(i, a.front()); + } +} + + +template +void test_halves() +{ + // x = n + 1/k => lur(x) = ((n; k - 1)) + // Note that this is a bit different that Kalpazidou (examine the half-open interval of definition carefully). + // One way to examine this definition is correct for rationals (it never happens for irrationals) + // is to consider i + 1/3. If you follow Kalpazidou, then you get ((i, 3, 0)); a zero digit! + // That's bad since it destroys uniqueness and also breaks the computation of the geometric mean. + for (int64_t i = -20; i < 20; ++i) { + Real x = i + Real(1)/Real(2); + auto luroth = luroth_expansion(x); + auto const & a = luroth.digits(); + CHECK_EQUAL(size_t(2), a.size()); + CHECK_EQUAL(i, a.front()); + CHECK_EQUAL(int64_t(1), a.back()); + } + + for (int64_t i = -20; i < 20; ++i) { + Real x = i + Real(1)/Real(4); + auto luroth = luroth_expansion(x); + auto const & a = luroth.digits(); + CHECK_EQUAL(size_t(2), a.size()); + CHECK_EQUAL(i, a.front()); + CHECK_EQUAL(int64_t(3), a.back()); + } + + for (int64_t i = -20; i < 20; ++i) { + Real x = i + Real(1)/Real(8); + auto luroth = luroth_expansion(x); + auto const & a = luroth.digits(); + CHECK_EQUAL(size_t(2), a.size()); + CHECK_EQUAL(i, a.front()); + CHECK_EQUAL(int64_t(7), a.back()); + } + // 1/3 is a pain because it's not representable: + Real x = Real(1)/Real(3); + auto luroth = luroth_expansion(x); + auto const & a = luroth.digits(); + CHECK_EQUAL(size_t(2), a.size()); + CHECK_EQUAL(int64_t(0), a.front()); + CHECK_EQUAL(int64_t(2), a.back()); +} + + +int main() +{ + test_integral(); + test_integral(); + test_integral(); + test_integral(); + + test_halves(); + test_halves(); + test_halves(); + test_halves(); + + #ifdef BOOST_HAS_FLOAT128 + test_integral(); + test_halves(); + #endif + return boost::math::test::report_errors(); +}