math/tools/luroth_expansion.hpp

0001 //  (C) Copyright Nick Thompson 2020.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_TOOLS_LUROTH_EXPANSION_HPP
0007 #define BOOST_MATH_TOOLS_LUROTH_EXPANSION_HPP
0008
0009 #include <vector>
0010 #include <ostream>
0011 #include <iomanip>
0012 #include <cmath>
0013 #include <limits>
0014 #include <cstdint>
0015 #include <stdexcept>
0016
0017 #include <boost/math/tools/is_standalone.hpp>
0018 #ifndef BOOST_MATH_STANDALONE
0019 #include <boost/config.hpp>
0020 #ifdef BOOST_MATH_NO_CXX17_IF_CONSTEXPR
0021 #error "The header <boost/math/norms.hpp> can only be used in C++17 and later."
0022 #endif
0023 #endif
0024
0025 namespace boost::math::tools {
0026
0027 template<typename Real, typename Z = int64_t>
0028 class luroth_expansion {
0029 public:
0030     luroth_expansion(Real x) : x_{x}
0031     {
0032         using std::floor;
0033         using std::abs;
0034         using std::sqrt;
0035         using std::isfinite;
0036         if (!isfinite(x))
0037         {
0038             throw std::domain_error("Cannot convert non-finites into a Luroth representation.");
0039         }
0040         d_.reserve(50);
0041         Real dn1 = floor(x);
0042         d_.push_back(static_cast<Z>(dn1));
0043         if (dn1 == x)
0044         {
0045            d_.shrink_to_fit();
0046            return;
0047         }
0048         // This attempts to follow the notation of:
0049         // "Khinchine's constant for Luroth Representation", by Sophia Kalpazidou.
0050         x = x - dn1;
0051         Real computed = dn1;
0052         Real prod = 1;
0053         // Let the error bound grow by 1 ULP/iteration.
0054         // I haven't done the error analysis to show that this is an expected rate of error growth,
0055         // but if you don't do this, you can easily get into an infinite loop.
0056         Real i = 1;
0057         Real scale = std::numeric_limits<Real>::epsilon()*abs(x_)/2;
0058         while (abs(x_ - computed) > (i++)*scale)
0059         {
0060            Real recip = 1/x;
0061            Real dn = floor(recip);
0062            // x = n + 1/k => lur(x) = ((n; k - 1))
0063            // Note that this is a bit different than Kalpazidou (examine the half-open interval of definition carefully).
0064            // One way to examine this definition is better for rationals (it never happens for irrationals)
0065            // is to consider i + 1/3. If you follow Kalpazidou, then you get ((i, 3, 0)); a zero digit!
0066            // That's bad since it destroys uniqueness and also breaks the computation of the geometric mean.
0067            if (recip == dn) {
0068               d_.push_back(static_cast<Z>(dn - 1));
0069               break;
0070            }
0071            d_.push_back(static_cast<Z>(dn));
0072            Real tmp = 1/(dn+1);
0073            computed += prod*tmp;
0074            prod *= tmp/dn;
0075            x = dn*(dn+1)*(x - tmp);
0076         }
0077
0078         for (size_t i = 1; i < d_.size(); ++i)
0079         {
0080             // Sanity check:
0081             if (d_[i] <= 0)
0082             {
0083                 throw std::domain_error("Found a digit <= 0; this is an error.");
0084             }
0085         }
0086         d_.shrink_to_fit();
0087     }
0088
0089
0090     const std::vector<Z>& digits() const {
0091       return d_;
0092     }
0093
0094     // Under the assumption of 'randomness', this mean converges to 2.2001610580.
0095     // See Finch, Mathematical Constants, section 1.8.1.
0096     Real digit_geometric_mean() const {
0097         if (d_.size() == 1) {
0098             return std::numeric_limits<Real>::quiet_NaN();
0099         }
0100         using std::log;
0101         using std::exp;
0102         Real g = 0;
0103         for (size_t i = 1; i < d_.size(); ++i) {
0104             g += log(static_cast<Real>(d_[i]));
0105         }
0106         return exp(g/(d_.size() - 1));
0107     }
0108
0109     template<typename T, typename Z2>
0110     friend std::ostream& operator<<(std::ostream& out, luroth_expansion<T, Z2>& scf);
0111
0112 private:
0113     const Real x_;
0114     std::vector<Z> d_;
0115 };
0116
0117
0118 template<typename Real, typename Z2>
0119 std::ostream& operator<<(std::ostream& out, luroth_expansion<Real, Z2>& luroth)
0120 {
0121    constexpr const int p = std::numeric_limits<Real>::max_digits10;
0122    if constexpr (p == 2147483647)
0123    {
0124       out << std::setprecision(luroth.x_.backend().precision());
0125    }
0126    else
0127    {
0128       out << std::setprecision(p);
0129    }
0130
0131    out << "((" << luroth.d_.front();
0132    if (luroth.d_.size() > 1)
0133    {
0134       out << "; ";
0135       for (size_t i = 1; i < luroth.d_.size() -1; ++i)
0136       {
0137          out << luroth.d_[i] << ", ";
0138       }
0139       out << luroth.d_.back();
0140    }
0141    out << "))";
0142    return out;
0143 }
0144
0145
0146 }
0147 #endif