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/config/no_tr1/cmath.hpp>
  14 #include <boost/cstdint.hpp>
  15 #include <boost/type_traits/integral_constant.hpp>
  16 #include <boost/mpl/if.hpp>
  17 #include <boost/math/tools/precision.hpp>
  18
  19 namespace boost{ namespace math{ namespace tools{
  20
  21 namespace detail
  22 {
  23
  24    template <class T>
  25    struct is_pair : public boost::false_type{};
  26
  27    template <class T, class U>
  28    struct is_pair<std::pair<T,U> > : public boost::true_type{};
  29
  30    template <class Gen>
  31    struct fraction_traits_simple
  32    {
  33        typedef typename Gen::result_type result_type;
  34        typedef typename Gen::result_type value_type;
  35
  36        static result_type a(const value_type&)
  37        {
  38           return 1;
  39        }
  40        static result_type b(const value_type& v)
  41        {
  42           return v;
  43        }
  44    };
  45
  46    template <class Gen>
  47    struct fraction_traits_pair
  48    {
  49        typedef typename Gen::result_type value_type;
  50        typedef typename value_type::first_type result_type;
  51
  52        static result_type a(const value_type& v)
  53        {
  54           return v.first;
  55        }
  56        static result_type b(const value_type& v)
  57        {
  58           return v.second;
  59        }
  60    };
  61
  62    template <class Gen>
  63    struct fraction_traits
  64        : public boost::mpl::if_c<
  65          is_pair<typename Gen::result_type>::value,
  66          fraction_traits_pair<Gen>,
  67          fraction_traits_simple<Gen> >::type
  68    {
  69    };
  70
  71 } // namespace detail
  72
  73 //
  74 // continued_fraction_b
  75 // Evaluates:
  76 //
  77 // b0 +       a1
  78 //      ---------------
  79 //      b1 +     a2
  80 //           ----------
  81 //           b2 +   a3
  82 //                -----
  83 //                b3 + ...
  84 //
  85 // Note that the first a0 returned by generator Gen is disarded.
  86 //
  87 template <class Gen, class U>
  88 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor, boost::uintmax_t& max_terms)
  89 {
  90    BOOST_MATH_STD_USING // ADL of std names
  91
  92    typedef detail::fraction_traits<Gen> traits;
  93    typedef typename traits::result_type result_type;
  94    typedef typename traits::value_type value_type;
  95
  96    result_type tiny = tools::min_value<result_type>();
  97
  98    value_type v = g();
  99
 100    result_type f, C, D, delta;
 101    f = traits::b(v);
 102    if(f == 0)
 103       f = tiny;
 104    C = f;
 105    D = 0;
 106
 107    boost::uintmax_t counter(max_terms);
 108
 109    do{
 110       v = g();
 111       D = traits::b(v) + traits::a(v) * D;
 112       if(D == 0)
 113          D = tiny;
 114       C = traits::b(v) + traits::a(v) / C;
 115       if(C == 0)
 116          C = tiny;
 117       D = 1/D;
 118       delta = C*D;
 119       f = f * delta;
 120    }while((fabs(delta - 1) > factor) && --counter);
 121
 122    max_terms = max_terms - counter;
 123
 124    return f;
 125 }
 126
 127 template <class Gen, class U>
 128 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor)
 129 {
 130    boost::uintmax_t max_terms = (std::numeric_limits<boost::uintmax_t>::max)();
 131    return continued_fraction_b(g, factor, max_terms);
 132 }
 133
 134 template <class Gen>
 135 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits)
 136 {
 137    BOOST_MATH_STD_USING // ADL of std names
 138
 139    typedef detail::fraction_traits<Gen> traits;
 140    typedef typename traits::result_type result_type;
 141
 142    result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
 143    boost::uintmax_t max_terms = (std::numeric_limits<boost::uintmax_t>::max)();
 144    return continued_fraction_b(g, factor, max_terms);
 145 }
 146
 147 template <class Gen>
 148 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits, boost::uintmax_t& max_terms)
 149 {
 150    BOOST_MATH_STD_USING // ADL of std names
 151
 152    typedef detail::fraction_traits<Gen> traits;
 153    typedef typename traits::result_type result_type;
 154
 155    result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
 156    return continued_fraction_b(g, factor, max_terms);
 157 }
 158
 159 //
 160 // continued_fraction_a
 161 // Evaluates:
 162 //
 163 //            a1
 164 //      ---------------
 165 //      b1 +     a2
 166 //           ----------
 167 //           b2 +   a3
 168 //                -----
 169 //                b3 + ...
 170 //
 171 // Note that the first a1 and b1 returned by generator Gen are both used.
 172 //
 173 template <class Gen, class U>
 174 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor, boost::uintmax_t& max_terms)
 175 {
 176    BOOST_MATH_STD_USING // ADL of std names
 177
 178    typedef detail::fraction_traits<Gen> traits;
 179    typedef typename traits::result_type result_type;
 180    typedef typename traits::value_type value_type;
 181
 182    result_type tiny = tools::min_value<result_type>();
 183
 184    value_type v = g();
 185
 186    result_type f, C, D, delta, a0;
 187    f = traits::b(v);
 188    a0 = traits::a(v);
 189    if(f == 0)
 190       f = tiny;
 191    C = f;
 192    D = 0;
 193
 194    boost::uintmax_t counter(max_terms);
 195
 196    do{
 197       v = g();
 198       D = traits::b(v) + traits::a(v) * D;
 199       if(D == 0)
 200          D = tiny;
 201       C = traits::b(v) + traits::a(v) / C;
 202       if(C == 0)
 203          C = tiny;
 204       D = 1/D;
 205       delta = C*D;
 206       f = f * delta;
 207    }while((fabs(delta - 1) > factor) && --counter);
 208
 209    max_terms = max_terms - counter;
 210
 211    return a0/f;
 212 }
 213
 214 template <class Gen, class U>
 215 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor)
 216 {
 217    boost::uintmax_t max_iter = (std::numeric_limits<boost::uintmax_t>::max)();
 218    return continued_fraction_a(g, factor, max_iter);
 219 }
 220
 221 template <class Gen>
 222 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits)
 223 {
 224    BOOST_MATH_STD_USING // ADL of std names
 225
 226    typedef detail::fraction_traits<Gen> traits;
 227    typedef typename traits::result_type result_type;
 228
 229    result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
 230    boost::uintmax_t max_iter = (std::numeric_limits<boost::uintmax_t>::max)();
 231
 232    return continued_fraction_a(g, factor, max_iter);
 233 }
 234
 235 template <class Gen>
 236 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits, boost::uintmax_t& max_terms)
 237 {
 238    BOOST_MATH_STD_USING // ADL of std names
 239
 240    typedef detail::fraction_traits<Gen> traits;
 241    typedef typename traits::result_type result_type;
 242
 243    result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
 244    return continued_fraction_a(g, factor, max_terms);
 245 }
 246
 247 } // namespace tools
 248 } // namespace math
 249 } // namespace boost
 250
 251 #endif // BOOST_MATH_TOOLS_FRACTION_INCLUDED
 252