Chris@16: //  (C) Copyright John Maddock 2005-2006.
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_MATH_TOOLS_FRACTION_INCLUDED
Chris@16: #define BOOST_MATH_TOOLS_FRACTION_INCLUDED
Chris@16: 
Chris@16: #ifdef _MSC_VER
Chris@16: #pragma once
Chris@16: #endif
Chris@16: 
Chris@16: #include <boost/config/no_tr1/cmath.hpp>
Chris@16: #include <boost/cstdint.hpp>
Chris@16: #include <boost/type_traits/integral_constant.hpp>
Chris@16: #include <boost/mpl/if.hpp>
Chris@16: #include <boost/math/tools/precision.hpp>
Chris@16: 
Chris@16: namespace boost{ namespace math{ namespace tools{
Chris@16: 
Chris@16: namespace detail
Chris@16: {
Chris@16: 
Chris@16:    template <class T>
Chris@16:    struct is_pair : public boost::false_type{};
Chris@16: 
Chris@16:    template <class T, class U>
Chris@16:    struct is_pair<std::pair<T,U> > : public boost::true_type{};
Chris@16: 
Chris@16:    template <class Gen>
Chris@16:    struct fraction_traits_simple
Chris@16:    {
Chris@16:        typedef typename Gen::result_type result_type;
Chris@16:        typedef typename Gen::result_type value_type;
Chris@16: 
Chris@16:        static result_type a(const value_type&)
Chris@16:        {
Chris@16:           return 1;
Chris@16:        }
Chris@16:        static result_type b(const value_type& v)
Chris@16:        {
Chris@16:           return v;
Chris@16:        }
Chris@16:    };
Chris@16: 
Chris@16:    template <class Gen>
Chris@16:    struct fraction_traits_pair
Chris@16:    {
Chris@16:        typedef typename Gen::result_type value_type;
Chris@16:        typedef typename value_type::first_type result_type;
Chris@16: 
Chris@16:        static result_type a(const value_type& v)
Chris@16:        {
Chris@16:           return v.first;
Chris@16:        }
Chris@16:        static result_type b(const value_type& v)
Chris@16:        {
Chris@16:           return v.second;
Chris@16:        }
Chris@16:    };
Chris@16: 
Chris@16:    template <class Gen>
Chris@16:    struct fraction_traits
Chris@16:        : public boost::mpl::if_c<
Chris@16:          is_pair<typename Gen::result_type>::value,
Chris@16:          fraction_traits_pair<Gen>,
Chris@16:          fraction_traits_simple<Gen> >::type
Chris@16:    {
Chris@16:    };
Chris@16: 
Chris@16: } // namespace detail
Chris@16: 
Chris@16: //
Chris@16: // continued_fraction_b
Chris@16: // Evaluates:
Chris@16: //
Chris@16: // b0 +       a1
Chris@16: //      ---------------
Chris@16: //      b1 +     a2
Chris@16: //           ----------
Chris@16: //           b2 +   a3
Chris@16: //                -----
Chris@16: //                b3 + ...
Chris@16: //
Chris@16: // Note that the first a0 returned by generator Gen is disarded.
Chris@16: //
Chris@16: template <class Gen, class U>
Chris@16: inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor, boost::uintmax_t& max_terms)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING // ADL of std names
Chris@16: 
Chris@16:    typedef detail::fraction_traits<Gen> traits;
Chris@16:    typedef typename traits::result_type result_type;
Chris@16:    typedef typename traits::value_type value_type;
Chris@16: 
Chris@16:    result_type tiny = tools::min_value<result_type>();
Chris@16: 
Chris@16:    value_type v = g();
Chris@16: 
Chris@16:    result_type f, C, D, delta;
Chris@16:    f = traits::b(v);
Chris@16:    if(f == 0)
Chris@16:       f = tiny;
Chris@16:    C = f;
Chris@16:    D = 0;
Chris@16: 
Chris@16:    boost::uintmax_t counter(max_terms);
Chris@16: 
Chris@16:    do{
Chris@16:       v = g();
Chris@16:       D = traits::b(v) + traits::a(v) * D;
Chris@16:       if(D == 0)
Chris@16:          D = tiny;
Chris@16:       C = traits::b(v) + traits::a(v) / C;
Chris@16:       if(C == 0)
Chris@16:          C = tiny;
Chris@16:       D = 1/D;
Chris@16:       delta = C*D;
Chris@16:       f = f * delta;
Chris@16:    }while((fabs(delta - 1) > factor) && --counter);
Chris@16: 
Chris@16:    max_terms = max_terms - counter;
Chris@16: 
Chris@16:    return f;
Chris@16: }
Chris@16: 
Chris@16: template <class Gen, class U>
Chris@16: inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor)
Chris@16: {
Chris@16:    boost::uintmax_t max_terms = (std::numeric_limits<boost::uintmax_t>::max)();
Chris@16:    return continued_fraction_b(g, factor, max_terms);
Chris@16: }
Chris@16: 
Chris@16: template <class Gen>
Chris@16: inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING // ADL of std names
Chris@16: 
Chris@16:    typedef detail::fraction_traits<Gen> traits;
Chris@16:    typedef typename traits::result_type result_type;
Chris@16: 
Chris@16:    result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
Chris@16:    boost::uintmax_t max_terms = (std::numeric_limits<boost::uintmax_t>::max)();
Chris@16:    return continued_fraction_b(g, factor, max_terms);
Chris@16: }
Chris@16: 
Chris@16: template <class Gen>
Chris@16: inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits, boost::uintmax_t& max_terms)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING // ADL of std names
Chris@16: 
Chris@16:    typedef detail::fraction_traits<Gen> traits;
Chris@16:    typedef typename traits::result_type result_type;
Chris@16: 
Chris@16:    result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
Chris@16:    return continued_fraction_b(g, factor, max_terms);
Chris@16: }
Chris@16: 
Chris@16: //
Chris@16: // continued_fraction_a
Chris@16: // Evaluates:
Chris@16: //
Chris@16: //            a1
Chris@16: //      ---------------
Chris@16: //      b1 +     a2
Chris@16: //           ----------
Chris@16: //           b2 +   a3
Chris@16: //                -----
Chris@16: //                b3 + ...
Chris@16: //
Chris@16: // Note that the first a1 and b1 returned by generator Gen are both used.
Chris@16: //
Chris@16: template <class Gen, class U>
Chris@16: inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor, boost::uintmax_t& max_terms)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING // ADL of std names
Chris@16: 
Chris@16:    typedef detail::fraction_traits<Gen> traits;
Chris@16:    typedef typename traits::result_type result_type;
Chris@16:    typedef typename traits::value_type value_type;
Chris@16: 
Chris@16:    result_type tiny = tools::min_value<result_type>();
Chris@16: 
Chris@16:    value_type v = g();
Chris@16: 
Chris@16:    result_type f, C, D, delta, a0;
Chris@16:    f = traits::b(v);
Chris@16:    a0 = traits::a(v);
Chris@16:    if(f == 0)
Chris@16:       f = tiny;
Chris@16:    C = f;
Chris@16:    D = 0;
Chris@16: 
Chris@16:    boost::uintmax_t counter(max_terms);
Chris@16: 
Chris@16:    do{
Chris@16:       v = g();
Chris@16:       D = traits::b(v) + traits::a(v) * D;
Chris@16:       if(D == 0)
Chris@16:          D = tiny;
Chris@16:       C = traits::b(v) + traits::a(v) / C;
Chris@16:       if(C == 0)
Chris@16:          C = tiny;
Chris@16:       D = 1/D;
Chris@16:       delta = C*D;
Chris@16:       f = f * delta;
Chris@16:    }while((fabs(delta - 1) > factor) && --counter);
Chris@16: 
Chris@16:    max_terms = max_terms - counter;
Chris@16: 
Chris@16:    return a0/f;
Chris@16: }
Chris@16: 
Chris@16: template <class Gen, class U>
Chris@16: inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor)
Chris@16: {
Chris@16:    boost::uintmax_t max_iter = (std::numeric_limits<boost::uintmax_t>::max)();
Chris@16:    return continued_fraction_a(g, factor, max_iter);
Chris@16: }
Chris@16: 
Chris@16: template <class Gen>
Chris@16: inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING // ADL of std names
Chris@16: 
Chris@16:    typedef detail::fraction_traits<Gen> traits;
Chris@16:    typedef typename traits::result_type result_type;
Chris@16: 
Chris@16:    result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
Chris@16:    boost::uintmax_t max_iter = (std::numeric_limits<boost::uintmax_t>::max)();
Chris@16: 
Chris@16:    return continued_fraction_a(g, factor, max_iter);
Chris@16: }
Chris@16: 
Chris@16: template <class Gen>
Chris@16: inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits, boost::uintmax_t& max_terms)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING // ADL of std names
Chris@16: 
Chris@16:    typedef detail::fraction_traits<Gen> traits;
Chris@16:    typedef typename traits::result_type result_type;
Chris@16: 
Chris@16:    result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
Chris@16:    return continued_fraction_a(g, factor, max_terms);
Chris@16: }
Chris@16: 
Chris@16: } // namespace tools
Chris@16: } // namespace math
Chris@16: } // namespace boost
Chris@16: 
Chris@16: #endif // BOOST_MATH_TOOLS_FRACTION_INCLUDED
Chris@16: