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