ceph/src/boost/boost/math/tools/luroth_expansion.hpp

   1 //  (C) Copyright Nick Thompson 2020.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_TOOLS_LUROTH_EXPANSION_HPP
   7 #define BOOST_MATH_TOOLS_LUROTH_EXPANSION_HPP
   8
   9 #include <vector>
  10 #include <ostream>
  11 #include <iomanip>
  12 #include <cmath>
  13 #include <limits>
  14 #include <stdexcept>
  15
  16 namespace boost::math::tools {
  17
  18 template<typename Real, typename Z = int64_t>
  19 class luroth_expansion {
  20 public:
  21     luroth_expansion(Real x) : x_{x}
  22     {
  23         using std::floor;
  24         using std::abs;
  25         using std::sqrt;
  26         using std::isfinite;
  27         if (!isfinite(x))
  28         {
  29             throw std::domain_error("Cannot convert non-finites into a Luroth representation.");
  30         }
  31         d_.reserve(50);
  32         Real dn1 = floor(x);
  33         d_.push_back(static_cast<Z>(dn1));
  34         if (dn1 == x)
  35         {
  36            d_.shrink_to_fit();
  37            return;
  38         }
  39         // This attempts to follow the notation of:
  40         // "Khinchine's constant for Luroth Representation", by Sophia Kalpazidou.
  41         x = x - dn1;
  42         Real computed = dn1;
  43         Real prod = 1;
  44         // Let the error bound grow by 1 ULP/iteration.
  45         // I haven't done the error analysis to show that this is an expected rate of error growth,
  46         // but if you don't do this, you can easily get into an infinite loop.
  47         Real i = 1;
  48         Real scale = std::numeric_limits<Real>::epsilon()*abs(x_)/2;
  49         while (abs(x_ - computed) > (i++)*scale)
  50         {
  51            Real recip = 1/x;
  52            Real dn = floor(recip);
  53            // x = n + 1/k => lur(x) = ((n; k - 1))
  54            // Note that this is a bit different than Kalpazidou (examine the half-open interval of definition carefully).
  55            // One way to examine this definition is better for rationals (it never happens for irrationals)
  56            // is to consider i + 1/3. If you follow Kalpazidou, then you get ((i, 3, 0)); a zero digit!
  57            // That's bad since it destroys uniqueness and also breaks the computation of the geometric mean.
  58            if (recip == dn) {
  59               d_.push_back(static_cast<Z>(dn - 1));
  60               break;
  61            }
  62            d_.push_back(static_cast<Z>(dn));
  63            Real tmp = 1/(dn+1);
  64            computed += prod*tmp;
  65            prod *= tmp/dn;
  66            x = dn*(dn+1)*(x - tmp);
  67         }
  68
  69         for (size_t i = 1; i < d_.size(); ++i)
  70         {
  71             // Sanity check:
  72             if (d_[i] <= 0)
  73             {
  74                 throw std::domain_error("Found a digit <= 0; this is an error.");
  75             }
  76         }
  77         d_.shrink_to_fit();
  78     }
  79
  80
  81     const std::vector<Z>& digits() const {
  82       return d_;
  83     }
  84
  85     // Under the assumption of 'randomness', this mean converges to 2.2001610580.
  86     // See Finch, Mathematical Constants, section 1.8.1.
  87     Real digit_geometric_mean() const {
  88         if (d_.size() == 1) {
  89             return std::numeric_limits<Real>::quiet_NaN();
  90         }
  91         using std::log;
  92         using std::exp;
  93         Real g = 0;
  94         for (size_t i = 1; i < d_.size(); ++i) {
  95             g += log(static_cast<Real>(d_[i]));
  96         }
  97         return exp(g/(d_.size() - 1));
  98     }
  99
 100     template<typename T, typename Z2>
 101     friend std::ostream& operator<<(std::ostream& out, luroth_expansion<T, Z2>& scf);
 102
 103 private:
 104     const Real x_;
 105     std::vector<Z> d_;
 106 };
 107
 108
 109 template<typename Real, typename Z2>
 110 std::ostream& operator<<(std::ostream& out, luroth_expansion<Real, Z2>& luroth)
 111 {
 112    constexpr const int p = std::numeric_limits<Real>::max_digits10;
 113    if constexpr (p == 2147483647)
 114    {
 115       out << std::setprecision(luroth.x_.backend().precision());
 116    }
 117    else
 118    {
 119       out << std::setprecision(p);
 120    }
 121
 122    out << "((" << luroth.d_.front();
 123    if (luroth.d_.size() > 1)
 124    {
 125       out << "; ";
 126       for (size_t i = 1; i < luroth.d_.size() -1; ++i)
 127       {
 128          out << luroth.d_[i] << ", ";
 129       }
 130       out << luroth.d_.back();
 131    }
 132    out << "))";
 133    return out;
 134 }
 135
 136
 137 }
 138 #endif