125 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_ENGEL_EXPANSION_HPP7#define BOOST_MATH_TOOLS_ENGEL_EXPANSION_HPP8 9#include <cmath>10#include <cstdint>11#include <vector>12#include <ostream>13#include <iomanip>14#include <limits>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 engel_expansion {29public:30 engel_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 an Engel expansion.");39 }40 41 if(x==0)42 {43 throw std::domain_error("Zero does not have an Engel expansion.");44 }45 a_.reserve(64);46 // Let the error bound grow by 1 ULP/iteration.47 // I haven't done the error analysis to show that this is an expected rate of error growth,48 // but if you don't do this, you can easily get into an infinite loop.49 Real i = 1;50 Real computed = 0;51 Real term = 1;52 Real scale = std::numeric_limits<Real>::epsilon()*abs(x_)/2;53 Real u = x;54 while (abs(x_ - computed) > (i++)*scale)55 {56 Real recip = 1/u;57 Real ak = ceil(recip);58 a_.push_back(static_cast<Z>(ak));59 u = u*ak - 1;60 if (u==0)61 {62 break;63 }64 term /= ak;65 computed += term;66 }67 68 for (size_t j = 1; j < a_.size(); ++j)69 {70 // Sanity check: This should only happen when wraparound occurs:71 if (a_[j] < a_[j-1])72 {73 throw std::domain_error("The digits of an Engel expansion must form a non-decreasing sequence; consider increasing the wide of the integer type.");74 }75 // Watch out for saturating behavior:76 if (a_[j] == (std::numeric_limits<Z>::max)())77 {78 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.");79 }80 }81 a_.shrink_to_fit();82 }83 84 85 const std::vector<Z>& digits() const86 {87 return a_;88 }89 90 template<typename T, typename Z2>91 friend std::ostream& operator<<(std::ostream& out, engel_expansion<T, Z2>& eng);92 93private:94 Real x_;95 std::vector<Z> a_;96};97 98 99template<typename Real, typename Z2>100std::ostream& operator<<(std::ostream& out, engel_expansion<Real, Z2>& engel)101{102 constexpr const int p = std::numeric_limits<Real>::max_digits10;103 if constexpr (p == 2147483647)104 {105 out << std::setprecision(engel.x_.backend().precision());106 }107 else108 {109 out << std::setprecision(p);110 }111 112 out << "{";113 for (size_t i = 0; i < engel.a_.size() - 1; ++i)114 {115 out << engel.a_[i] << ", ";116 }117 out << engel.a_.back();118 out << "}";119 return out;120}121 122 123}124#endif125