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

   1 //  (C) Copyright John Maddock 2005-2006.
   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_FRACTION_INCLUDED
   7 #define BOOST_MATH_TOOLS_FRACTION_INCLUDED
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 #include <boost/math/tools/precision.hpp>
  14 #include <boost/math/tools/complex.hpp>
  15 #include <type_traits>
  16 #include <cstdint>
  17 #include <cmath>
  18
  19 namespace boost{ namespace math{ namespace tools{
  20
  21 namespace detail
  22 {
  23
  24    template <typename T>
  25    struct is_pair : public std::false_type{};
  26
  27    template <typename T, typename U>
  28    struct is_pair<std::pair<T,U>> : public std::true_type{};
  29
  30    template <typename Gen>
  31    struct fraction_traits_simple
  32    {
  33       using result_type = typename Gen::result_type;
  34       using  value_type = typename Gen::result_type;
  35
  36       static result_type a(const value_type&) BOOST_MATH_NOEXCEPT(value_type)
  37       {
  38          return 1;
  39       }
  40       static result_type b(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)
  41       {
  42          return v;
  43       }
  44    };
  45
  46    template <typename Gen>
  47    struct fraction_traits_pair
  48    {
  49       using  value_type = typename Gen::result_type;
  50       using result_type = typename value_type::first_type;
  51
  52       static result_type a(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)
  53       {
  54          return v.first;
  55       }
  56       static result_type b(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)
  57       {
  58          return v.second;
  59       }
  60    };
  61
  62    template <typename Gen>
  63    struct fraction_traits
  64        : public std::conditional<
  65          is_pair<typename Gen::result_type>::value,
  66          fraction_traits_pair<Gen>,
  67          fraction_traits_simple<Gen>>::type
  68    {
  69    };
  70
  71    template <typename T, bool = is_complex_type<T>::value>
  72    struct tiny_value
  73    {
  74       // For float, double, and long double, 1/min_value<T>() is finite.
  75       // But for mpfr_float and cpp_bin_float, 1/min_value<T>() is inf.
  76       // Multiply the min by 16 so that the reciprocal doesn't overflow.
  77       static T get() {
  78          return 16*tools::min_value<T>();
  79       }
  80    };
  81    template <typename T>
  82    struct tiny_value<T, true>
  83    {
  84       using value_type = typename T::value_type;
  85       static T get() {
  86          return 16*tools::min_value<value_type>();
  87       }
  88    };
  89
  90 } // namespace detail
  91
  92 //
  93 // continued_fraction_b
  94 // Evaluates:
  95 //
  96 // b0 +       a1
  97 //      ---------------
  98 //      b1 +     a2
  99 //           ----------
 100 //           b2 +   a3
 101 //                -----
 102 //                b3 + ...
 103 //
 104 // Note that the first a0 returned by generator Gen is discarded.
 105 //
 106 template <typename Gen, typename U>
 107 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor, std::uintmax_t& max_terms)
 108       noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
 109 {
 110    BOOST_MATH_STD_USING // ADL of std names
 111
 112    using traits = detail::fraction_traits<Gen>;
 113    using result_type = typename traits::result_type;
 114    using value_type = typename traits::value_type;
 115    using integer_type = typename integer_scalar_type<result_type>::type;
 116    using scalar_type = typename scalar_type<result_type>::type;
 117
 118    integer_type const zero(0), one(1);
 119
 120    result_type tiny = detail::tiny_value<result_type>::get();
 121    scalar_type terminator = abs(factor);
 122
 123    value_type v = g();
 124
 125    result_type f, C, D, delta;
 126    f = traits::b(v);
 127    if(f == zero)
 128       f = tiny;
 129    C = f;
 130    D = 0;
 131
 132    std::uintmax_t counter(max_terms);
 133    do{
 134       v = g();
 135       D = traits::b(v) + traits::a(v) * D;
 136       if(D == result_type(0))
 137          D = tiny;
 138       C = traits::b(v) + traits::a(v) / C;
 139       if(C == zero)
 140          C = tiny;
 141       D = one/D;
 142       delta = C*D;
 143       f = f * delta;
 144    }while((abs(delta - one) > terminator) && --counter);
 145
 146    max_terms = max_terms - counter;
 147
 148    return f;
 149 }
 150
 151 template <typename Gen, typename U>
 152 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor)
 153    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
 154 {
 155    std::uintmax_t max_terms = (std::numeric_limits<std::uintmax_t>::max)();
 156    return continued_fraction_b(g, factor, max_terms);
 157 }
 158
 159 template <typename Gen>
 160 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits)
 161    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
 162 {
 163    BOOST_MATH_STD_USING // ADL of std names
 164
 165    using traits = detail::fraction_traits<Gen>;
 166    using result_type = typename traits::result_type;
 167
 168    result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
 169    std::uintmax_t max_terms = (std::numeric_limits<std::uintmax_t>::max)();
 170    return continued_fraction_b(g, factor, max_terms);
 171 }
 172
 173 template <typename Gen>
 174 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits, std::uintmax_t& max_terms)
 175    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
 176 {
 177    BOOST_MATH_STD_USING // ADL of std names
 178
 179    using traits = detail::fraction_traits<Gen>;
 180    using result_type = typename traits::result_type;
 181
 182    result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
 183    return continued_fraction_b(g, factor, max_terms);
 184 }
 185
 186 //
 187 // continued_fraction_a
 188 // Evaluates:
 189 //
 190 //            a1
 191 //      ---------------
 192 //      b1 +     a2
 193 //           ----------
 194 //           b2 +   a3
 195 //                -----
 196 //                b3 + ...
 197 //
 198 // Note that the first a1 and b1 returned by generator Gen are both used.
 199 //
 200 template <typename Gen, typename U>
 201 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor, std::uintmax_t& max_terms)
 202    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
 203 {
 204    BOOST_MATH_STD_USING // ADL of std names
 205
 206    using traits = detail::fraction_traits<Gen>;
 207    using result_type = typename traits::result_type;
 208    using value_type = typename traits::value_type;
 209    using integer_type = typename integer_scalar_type<result_type>::type;
 210    using scalar_type = typename scalar_type<result_type>::type;
 211
 212    integer_type const zero(0), one(1);
 213
 214    result_type tiny = detail::tiny_value<result_type>::get();
 215    scalar_type terminator = abs(factor);
 216
 217    value_type v = g();
 218
 219    result_type f, C, D, delta, a0;
 220    f = traits::b(v);
 221    a0 = traits::a(v);
 222    if(f == zero)
 223       f = tiny;
 224    C = f;
 225    D = 0;
 226
 227    std::uintmax_t counter(max_terms);
 228
 229    do{
 230       v = g();
 231       D = traits::b(v) + traits::a(v) * D;
 232       if(D == zero)
 233          D = tiny;
 234       C = traits::b(v) + traits::a(v) / C;
 235       if(C == zero)
 236          C = tiny;
 237       D = one/D;
 238       delta = C*D;
 239       f = f * delta;
 240    }while((abs(delta - one) > terminator) && --counter);
 241
 242    max_terms = max_terms - counter;
 243
 244    return a0/f;
 245 }
 246
 247 template <typename Gen, typename U>
 248 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor)
 249    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
 250 {
 251    std::uintmax_t max_iter = (std::numeric_limits<std::uintmax_t>::max)();
 252    return continued_fraction_a(g, factor, max_iter);
 253 }
 254
 255 template <typename Gen>
 256 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits)
 257    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
 258 {
 259    BOOST_MATH_STD_USING // ADL of std names
 260
 261    typedef detail::fraction_traits<Gen> traits;
 262    typedef typename traits::result_type result_type;
 263
 264    result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
 265    std::uintmax_t max_iter = (std::numeric_limits<std::uintmax_t>::max)();
 266
 267    return continued_fraction_a(g, factor, max_iter);
 268 }
 269
 270 template <typename Gen>
 271 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits, std::uintmax_t& max_terms)
 272    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
 273 {
 274    BOOST_MATH_STD_USING // ADL of std names
 275
 276    using traits = detail::fraction_traits<Gen>;
 277    using result_type = typename traits::result_type;
 278
 279    result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
 280    return continued_fraction_a(g, factor, max_terms);
 281 }
 282
 283 } // namespace tools
 284 } // namespace math
 285 } // namespace boost
 286
 287 #endif // BOOST_MATH_TOOLS_FRACTION_INCLUDED