math/tools/fraction.hpp

0001 //  (C) Copyright John Maddock 2005-2006.
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_FRACTION_INCLUDED
0007 #define BOOST_MATH_TOOLS_FRACTION_INCLUDED
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #endif
0012
0013 #include <boost/math/tools/precision.hpp>
0014 #include <boost/math/tools/complex.hpp>
0015 #include <type_traits>
0016 #include <cstdint>
0017 #include <cmath>
0018
0019 namespace boost{ namespace math{ namespace tools{
0020
0021 namespace detail
0022 {
0023
0024    template <typename T>
0025    struct is_pair : public std::false_type{};
0026
0027    template <typename T, typename U>
0028    struct is_pair<std::pair<T,U>> : public std::true_type{};
0029
0030    template <typename Gen>
0031    struct fraction_traits_simple
0032    {
0033       using result_type = typename Gen::result_type;
0034       using  value_type = typename Gen::result_type;
0035
0036       static result_type a(const value_type&) BOOST_MATH_NOEXCEPT(value_type)
0037       {
0038          return 1;
0039       }
0040       static result_type b(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)
0041       {
0042          return v;
0043       }
0044    };
0045
0046    template <typename Gen>
0047    struct fraction_traits_pair
0048    {
0049       using  value_type = typename Gen::result_type;
0050       using result_type = typename value_type::first_type;
0051
0052       static result_type a(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)
0053       {
0054          return v.first;
0055       }
0056       static result_type b(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)
0057       {
0058          return v.second;
0059       }
0060    };
0061
0062    template <typename Gen>
0063    struct fraction_traits
0064        : public std::conditional<
0065          is_pair<typename Gen::result_type>::value,
0066          fraction_traits_pair<Gen>,
0067          fraction_traits_simple<Gen>>::type
0068    {
0069    };
0070
0071    template <typename T, bool = is_complex_type<T>::value>
0072    struct tiny_value
0073    {
0074       // For float, double, and long double, 1/min_value<T>() is finite.
0075       // But for mpfr_float and cpp_bin_float, 1/min_value<T>() is inf.
0076       // Multiply the min by 16 so that the reciprocal doesn't overflow.
0077       static T get() {
0078          return 16*tools::min_value<T>();
0079       }
0080    };
0081    template <typename T>
0082    struct tiny_value<T, true>
0083    {
0084       using value_type = typename T::value_type;
0085       static T get() {
0086          return 16*tools::min_value<value_type>();
0087       }
0088    };
0089
0090 } // namespace detail
0091
0092 //
0093 // continued_fraction_b
0094 // Evaluates:
0095 //
0096 // b0 +       a1
0097 //      ---------------
0098 //      b1 +     a2
0099 //           ----------
0100 //           b2 +   a3
0101 //                -----
0102 //                b3 + ...
0103 //
0104 // Note that the first a0 returned by generator Gen is discarded.
0105 //
0106 template <typename Gen, typename U>
0107 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor, std::uintmax_t& max_terms)
0108       noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
0109 {
0110    BOOST_MATH_STD_USING // ADL of std names
0111
0112    using traits = detail::fraction_traits<Gen>;
0113    using result_type = typename traits::result_type;
0114    using value_type = typename traits::value_type;
0115    using integer_type = typename integer_scalar_type<result_type>::type;
0116    using scalar_type = typename scalar_type<result_type>::type;
0117
0118    integer_type const zero(0), one(1);
0119
0120    result_type tiny = detail::tiny_value<result_type>::get();
0121    scalar_type terminator = abs(factor);
0122
0123    value_type v = g();
0124
0125    result_type f, C, D, delta;
0126    f = traits::b(v);
0127    if(f == zero)
0128       f = tiny;
0129    C = f;
0130    D = 0;
0131
0132    std::uintmax_t counter(max_terms);
0133    do{
0134       v = g();
0135       D = traits::b(v) + traits::a(v) * D;
0136       if(D == result_type(0))
0137          D = tiny;
0138       C = traits::b(v) + traits::a(v) / C;
0139       if(C == zero)
0140          C = tiny;
0141       D = one/D;
0142       delta = C*D;
0143       f = f * delta;
0144    }while((abs(delta - one) > terminator) && --counter);
0145
0146    max_terms = max_terms - counter;
0147
0148    return f;
0149 }
0150
0151 template <typename Gen, typename U>
0152 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor)
0153    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
0154 {
0155    std::uintmax_t max_terms = (std::numeric_limits<std::uintmax_t>::max)();
0156    return continued_fraction_b(g, factor, max_terms);
0157 }
0158
0159 template <typename Gen>
0160 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits)
0161    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
0162 {
0163    BOOST_MATH_STD_USING // ADL of std names
0164
0165    using traits = detail::fraction_traits<Gen>;
0166    using result_type = typename traits::result_type;
0167
0168    result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
0169    std::uintmax_t max_terms = (std::numeric_limits<std::uintmax_t>::max)();
0170    return continued_fraction_b(g, factor, max_terms);
0171 }
0172
0173 template <typename Gen>
0174 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits, std::uintmax_t& max_terms)
0175    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
0176 {
0177    BOOST_MATH_STD_USING // ADL of std names
0178
0179    using traits = detail::fraction_traits<Gen>;
0180    using result_type = typename traits::result_type;
0181
0182    result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
0183    return continued_fraction_b(g, factor, max_terms);
0184 }
0185
0186 //
0187 // continued_fraction_a
0188 // Evaluates:
0189 //
0190 //            a1
0191 //      ---------------
0192 //      b1 +     a2
0193 //           ----------
0194 //           b2 +   a3
0195 //                -----
0196 //                b3 + ...
0197 //
0198 // Note that the first a1 and b1 returned by generator Gen are both used.
0199 //
0200 template <typename Gen, typename U>
0201 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor, std::uintmax_t& max_terms)
0202    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
0203 {
0204    BOOST_MATH_STD_USING // ADL of std names
0205
0206    using traits = detail::fraction_traits<Gen>;
0207    using result_type = typename traits::result_type;
0208    using value_type = typename traits::value_type;
0209    using integer_type = typename integer_scalar_type<result_type>::type;
0210    using scalar_type = typename scalar_type<result_type>::type;
0211
0212    integer_type const zero(0), one(1);
0213
0214    result_type tiny = detail::tiny_value<result_type>::get();
0215    scalar_type terminator = abs(factor);
0216
0217    value_type v = g();
0218
0219    result_type f, C, D, delta, a0;
0220    f = traits::b(v);
0221    a0 = traits::a(v);
0222    if(f == zero)
0223       f = tiny;
0224    C = f;
0225    D = 0;
0226
0227    std::uintmax_t counter(max_terms);
0228
0229    do{
0230       v = g();
0231       D = traits::b(v) + traits::a(v) * D;
0232       if(D == zero)
0233          D = tiny;
0234       C = traits::b(v) + traits::a(v) / C;
0235       if(C == zero)
0236          C = tiny;
0237       D = one/D;
0238       delta = C*D;
0239       f = f * delta;
0240    }while((abs(delta - one) > terminator) && --counter);
0241
0242    max_terms = max_terms - counter;
0243
0244    return a0/f;
0245 }
0246
0247 template <typename Gen, typename U>
0248 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor)
0249    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
0250 {
0251    std::uintmax_t max_iter = (std::numeric_limits<std::uintmax_t>::max)();
0252    return continued_fraction_a(g, factor, max_iter);
0253 }
0254
0255 template <typename Gen>
0256 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits)
0257    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
0258 {
0259    BOOST_MATH_STD_USING // ADL of std names
0260
0261    typedef detail::fraction_traits<Gen> traits;
0262    typedef typename traits::result_type result_type;
0263
0264    result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
0265    std::uintmax_t max_iter = (std::numeric_limits<std::uintmax_t>::max)();
0266
0267    return continued_fraction_a(g, factor, max_iter);
0268 }
0269
0270 template <typename Gen>
0271 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits, std::uintmax_t& max_terms)
0272    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
0273 {
0274    BOOST_MATH_STD_USING // ADL of std names
0275
0276    using traits = detail::fraction_traits<Gen>;
0277    using result_type = typename traits::result_type;
0278
0279    result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
0280    return continued_fraction_a(g, factor, max_terms);
0281 }
0282
0283 } // namespace tools
0284 } // namespace math
0285 } // namespace boost
0286
0287 #endif // BOOST_MATH_TOOLS_FRACTION_INCLUDED