math/tools/engel_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_ENGEL_EXPANSION_HPP
0007 #define BOOST_MATH_TOOLS_ENGEL_EXPANSION_HPP
0008
0009 #include <cmath>
0010 #include <cstdint>
0011 #include <vector>
0012 #include <ostream>
0013 #include <iomanip>
0014 #include <limits>
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_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 engel_expansion {
0029 public:
0030     engel_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 an Engel expansion.");
0039         }
0040
0041         if(x==0)
0042         {
0043             throw std::domain_error("Zero does not have an Engel expansion.");
0044         }
0045         a_.reserve(64);
0046         // Let the error bound grow by 1 ULP/iteration.
0047         // I haven't done the error analysis to show that this is an expected rate of error growth,
0048         // but if you don't do this, you can easily get into an infinite loop.
0049         Real i = 1;
0050         Real computed = 0;
0051         Real term = 1;
0052         Real scale = std::numeric_limits<Real>::epsilon()*abs(x_)/2;
0053         Real u = x;
0054         while (abs(x_ - computed) > (i++)*scale)
0055         {
0056             Real recip = 1/u;
0057             Real ak = ceil(recip);
0058             a_.push_back(static_cast<Z>(ak));
0059             u = u*ak - 1;
0060             if (u==0)
0061             {
0062                 break;
0063             }
0064             term /= ak;
0065             computed += term;
0066         }
0067
0068         for (size_t j = 1; j < a_.size(); ++j)
0069         {
0070             // Sanity check: This should only happen when wraparound occurs:
0071             if (a_[j] < a_[j-1])
0072             {
0073                 throw std::domain_error("The digits of an Engel expansion must form a non-decreasing sequence; consider increasing the wide of the integer type.");
0074             }
0075             // Watch out for saturating behavior:
0076             if (a_[j] == (std::numeric_limits<Z>::max)())
0077             {
0078                 throw std::domain_error("The integer type Z does not have enough width to hold the terms of the Engel expansion; please widen the type.");
0079             }
0080         }
0081         a_.shrink_to_fit();
0082     }
0083
0084
0085     const std::vector<Z>& digits() const
0086     {
0087         return a_;
0088     }
0089
0090     template<typename T, typename Z2>
0091     friend std::ostream& operator<<(std::ostream& out, engel_expansion<T, Z2>& eng);
0092
0093 private:
0094     Real x_;
0095     std::vector<Z> a_;
0096 };
0097
0098
0099 template<typename Real, typename Z2>
0100 std::ostream& operator<<(std::ostream& out, engel_expansion<Real, Z2>& engel)
0101 {
0102     constexpr const int p = std::numeric_limits<Real>::max_digits10;
0103     if constexpr (p == 2147483647)
0104     {
0105         out << std::setprecision(engel.x_.backend().precision());
0106     }
0107     else
0108     {
0109         out << std::setprecision(p);
0110     }
0111
0112     out << "{";
0113     for (size_t i = 0; i < engel.a_.size() - 1; ++i)
0114     {
0115         out << engel.a_[i] << ", ";
0116     }
0117     out << engel.a_.back();
0118     out << "}";
0119     return out;
0120 }
0121
0122
0123 }
0124 #endif