math/tools/fraction.hpp

0001 //  (C) Copyright John Maddock 2005-2006.
0002 //  (C) Copyright Matt Borland 2024.
0003 //  Use, modification and distribution are subject to the
0004 //  Boost Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_MATH_TOOLS_FRACTION_INCLUDED
0008 #define BOOST_MATH_TOOLS_FRACTION_INCLUDED
0009
0010 #ifdef _MSC_VER
0011 #pragma once
0012 #endif
0013
0014 #include <boost/math/tools/config.hpp>
0015 #include <boost/math/tools/type_traits.hpp>
0016 #include <boost/math/tools/numeric_limits.hpp>
0017 #include <boost/math/tools/tuple.hpp>
0018 #include <boost/math/tools/precision.hpp>
0019 #include <boost/math/tools/complex.hpp>
0020 #include <boost/math/tools/cstdint.hpp>
0021
0022 namespace boost{ namespace math{ namespace tools{
0023
0024 namespace detail
0025 {
0026
0027    template <typename T>
0028    struct is_pair : public boost::math::false_type{};
0029
0030    template <typename T, typename U>
0031    struct is_pair<boost::math::pair<T,U>> : public boost::math::true_type{};
0032
0033    template <typename Gen>
0034    struct fraction_traits_simple
0035    {
0036       using result_type = typename Gen::result_type;
0037       using  value_type = typename Gen::result_type;
0038
0039       BOOST_MATH_GPU_ENABLED static result_type a(const value_type&) BOOST_MATH_NOEXCEPT(value_type)
0040       {
0041          return 1;
0042       }
0043       BOOST_MATH_GPU_ENABLED static result_type b(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)
0044       {
0045          return v;
0046       }
0047    };
0048
0049    template <typename Gen>
0050    struct fraction_traits_pair
0051    {
0052       using  value_type = typename Gen::result_type;
0053       using result_type = typename value_type::first_type;
0054
0055       BOOST_MATH_GPU_ENABLED static result_type a(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)
0056       {
0057          return v.first;
0058       }
0059       BOOST_MATH_GPU_ENABLED static result_type b(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)
0060       {
0061          return v.second;
0062       }
0063    };
0064
0065    template <typename Gen>
0066    struct fraction_traits
0067        : public boost::math::conditional<
0068          is_pair<typename Gen::result_type>::value,
0069          fraction_traits_pair<Gen>,
0070          fraction_traits_simple<Gen>>::type
0071    {
0072    };
0073
0074    template <typename T, bool = is_complex_type<T>::value>
0075    struct tiny_value
0076    {
0077       // For float, double, and long double, 1/min_value<T>() is finite.
0078       // But for mpfr_float and cpp_bin_float, 1/min_value<T>() is inf.
0079       // Multiply the min by 16 so that the reciprocal doesn't overflow.
0080       BOOST_MATH_GPU_ENABLED static T get() {
0081          return 16*tools::min_value<T>();
0082       }
0083    };
0084    template <typename T>
0085    struct tiny_value<T, true>
0086    {
0087       using value_type = typename T::value_type;
0088       BOOST_MATH_GPU_ENABLED static T get() {
0089          return 16*tools::min_value<value_type>();
0090       }
0091    };
0092
0093 } // namespace detail
0094
0095 namespace detail {
0096
0097 //
0098 // continued_fraction_b
0099 // Evaluates:
0100 //
0101 // b0 +       a1
0102 //      ---------------
0103 //      b1 +     a2
0104 //           ----------
0105 //           b2 +   a3
0106 //                -----
0107 //                b3 + ...
0108 //
0109 // Note that the first a0 returned by generator Gen is discarded.
0110 //
0111
0112 template <typename Gen, typename U>
0113 BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b_impl(Gen& g, const U& factor, boost::math::uintmax_t& max_terms)
0114       noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
0115       #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0116       // SYCL can not handle this condition so we only check float on that platform
0117       && noexcept(std::declval<Gen>()())
0118       #endif
0119       )
0120 {
0121    BOOST_MATH_STD_USING // ADL of std names
0122
0123    using traits = detail::fraction_traits<Gen>;
0124    using result_type = typename traits::result_type;
0125    using value_type = typename traits::value_type;
0126    using integer_type = typename integer_scalar_type<result_type>::type;
0127    using scalar_type = typename scalar_type<result_type>::type;
0128
0129    integer_type const zero(0), one(1);
0130
0131    result_type tiny = detail::tiny_value<result_type>::get();
0132    scalar_type terminator = abs(factor);
0133
0134    value_type v = g();
0135
0136    result_type f, C, D, delta;
0137    f = traits::b(v);
0138    if(f == zero)
0139       f = tiny;
0140    C = f;
0141    D = 0;
0142
0143    boost::math::uintmax_t counter(max_terms);
0144    do{
0145       v = g();
0146       D = traits::b(v) + traits::a(v) * D;
0147       if(D == result_type(0))
0148          D = tiny;
0149       C = traits::b(v) + traits::a(v) / C;
0150       if(C == zero)
0151          C = tiny;
0152       D = one/D;
0153       delta = C*D;
0154       f = f * delta;
0155    }while((abs(delta - one) > terminator) && --counter);
0156
0157    max_terms = max_terms - counter;
0158
0159    return f;
0160 }
0161
0162 } // namespace detail
0163
0164 template <typename Gen, typename U>
0165 BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor, boost::math::uintmax_t& max_terms)
0166    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
0167          #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0168          && noexcept(std::declval<Gen>()())
0169          #endif
0170          )
0171 {
0172    return detail::continued_fraction_b_impl(g, factor, max_terms);
0173 }
0174
0175 template <typename Gen, typename U>
0176 BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor)
0177    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
0178    #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0179    && noexcept(std::declval<Gen>()())
0180    #endif
0181    )
0182 {
0183    boost::math::uintmax_t max_terms = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();
0184    return detail::continued_fraction_b_impl(g, factor, max_terms);
0185 }
0186
0187 template <typename Gen>
0188 BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits)
0189    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
0190    #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0191    && noexcept(std::declval<Gen>()())
0192    #endif
0193    )
0194 {
0195    BOOST_MATH_STD_USING // ADL of std names
0196
0197    using traits = detail::fraction_traits<Gen>;
0198    using result_type = typename traits::result_type;
0199
0200    result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
0201    boost::math::uintmax_t max_terms = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();
0202    return detail::continued_fraction_b_impl(g, factor, max_terms);
0203 }
0204
0205 template <typename Gen>
0206 BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits, boost::math::uintmax_t& max_terms)
0207    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
0208    #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0209    && noexcept(std::declval<Gen>()())
0210    #endif
0211    )
0212 {
0213    BOOST_MATH_STD_USING // ADL of std names
0214
0215    using traits = detail::fraction_traits<Gen>;
0216    using result_type = typename traits::result_type;
0217
0218    result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
0219    return detail::continued_fraction_b_impl(g, factor, max_terms);
0220 }
0221
0222 namespace detail {
0223
0224 //
0225 // continued_fraction_a
0226 // Evaluates:
0227 //
0228 //            a1
0229 //      ---------------
0230 //      b1 +     a2
0231 //           ----------
0232 //           b2 +   a3
0233 //                -----
0234 //                b3 + ...
0235 //
0236 // Note that the first a1 and b1 returned by generator Gen are both used.
0237 //
0238 template <typename Gen, typename U>
0239 BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a_impl(Gen& g, const U& factor, boost::math::uintmax_t& max_terms)
0240    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
0241    #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0242    && noexcept(std::declval<Gen>()())
0243    #endif
0244    )
0245 {
0246    BOOST_MATH_STD_USING // ADL of std names
0247
0248    using traits = detail::fraction_traits<Gen>;
0249    using result_type = typename traits::result_type;
0250    using value_type = typename traits::value_type;
0251    using integer_type = typename integer_scalar_type<result_type>::type;
0252    using scalar_type = typename scalar_type<result_type>::type;
0253
0254    integer_type const zero(0), one(1);
0255
0256    result_type tiny = detail::tiny_value<result_type>::get();
0257    scalar_type terminator = abs(factor);
0258
0259    value_type v = g();
0260
0261    result_type f, C, D, delta, a0;
0262    f = traits::b(v);
0263    a0 = traits::a(v);
0264    if(f == zero)
0265       f = tiny;
0266    C = f;
0267    D = 0;
0268
0269    boost::math::uintmax_t counter(max_terms);
0270
0271    do{
0272       v = g();
0273       D = traits::b(v) + traits::a(v) * D;
0274       if(D == zero)
0275          D = tiny;
0276       C = traits::b(v) + traits::a(v) / C;
0277       if(C == zero)
0278          C = tiny;
0279       D = one/D;
0280       delta = C*D;
0281       f = f * delta;
0282    }while((abs(delta - one) > terminator) && --counter);
0283
0284    max_terms = max_terms - counter;
0285
0286    return a0/f;
0287 }
0288
0289 } // namespace detail
0290
0291 template <typename Gen, typename U>
0292 BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor, boost::math::uintmax_t& max_terms)
0293    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
0294    #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0295    && noexcept(std::declval<Gen>()())
0296    #endif
0297    )
0298 {
0299    return detail::continued_fraction_a_impl(g, factor, max_terms);
0300 }
0301
0302 template <typename Gen, typename U>
0303 BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor)
0304    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
0305    #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0306    && noexcept(std::declval<Gen>()())
0307    #endif
0308    )
0309 {
0310    boost::math::uintmax_t max_iter = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();
0311    return detail::continued_fraction_a_impl(g, factor, max_iter);
0312 }
0313
0314 template <typename Gen>
0315 BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits)
0316    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
0317    #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0318    && noexcept(std::declval<Gen>()())
0319    #endif
0320    )
0321 {
0322    BOOST_MATH_STD_USING // ADL of std names
0323
0324    typedef detail::fraction_traits<Gen> traits;
0325    typedef typename traits::result_type result_type;
0326
0327    result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
0328    boost::math::uintmax_t max_iter = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();
0329
0330    return detail::continued_fraction_a_impl(g, factor, max_iter);
0331 }
0332
0333 template <typename Gen>
0334 BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits, boost::math::uintmax_t& max_terms)
0335    noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
0336    #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0337    && noexcept(std::declval<Gen>()())
0338    #endif
0339    )
0340 {
0341    BOOST_MATH_STD_USING // ADL of std names
0342
0343    using traits = detail::fraction_traits<Gen>;
0344    using result_type = typename traits::result_type;
0345
0346    result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
0347    return detail::continued_fraction_a_impl(g, factor, max_terms);
0348 }
0349
0350 } // namespace tools
0351 } // namespace math
0352 } // namespace boost
0353
0354 #endif // BOOST_MATH_TOOLS_FRACTION_INCLUDED