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