355 lines · plain
1// (C) Copyright John Maddock 2005-2006.2// (C) Copyright Matt Borland 2024.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_MATH_TOOLS_FRACTION_INCLUDED8#define BOOST_MATH_TOOLS_FRACTION_INCLUDED9 10#ifdef _MSC_VER11#pragma once12#endif13 14#include <boost/math/tools/config.hpp>15#include <boost/math/tools/type_traits.hpp>16#include <boost/math/tools/numeric_limits.hpp>17#include <boost/math/tools/tuple.hpp>18#include <boost/math/tools/precision.hpp>19#include <boost/math/tools/complex.hpp>20#include <boost/math/tools/cstdint.hpp>21 22namespace boost{ namespace math{ namespace tools{23 24namespace detail25{26 27 template <typename T>28 struct is_pair : public boost::math::false_type{};29 30 template <typename T, typename U>31 struct is_pair<boost::math::pair<T,U>> : public boost::math::true_type{};32 33 template <typename Gen>34 struct fraction_traits_simple35 {36 using result_type = typename Gen::result_type;37 using value_type = typename Gen::result_type;38 39 BOOST_MATH_GPU_ENABLED static result_type a(const value_type&) BOOST_MATH_NOEXCEPT(value_type)40 {41 return 1;42 }43 BOOST_MATH_GPU_ENABLED static result_type b(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)44 {45 return v;46 }47 };48 49 template <typename Gen>50 struct fraction_traits_pair51 {52 using value_type = typename Gen::result_type;53 using result_type = typename value_type::first_type;54 55 BOOST_MATH_GPU_ENABLED static result_type a(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)56 {57 return v.first;58 }59 BOOST_MATH_GPU_ENABLED static result_type b(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)60 {61 return v.second;62 }63 };64 65 template <typename Gen>66 struct fraction_traits67 : public boost::math::conditional<68 is_pair<typename Gen::result_type>::value,69 fraction_traits_pair<Gen>,70 fraction_traits_simple<Gen>>::type71 {72 };73 74 template <typename T, bool = is_complex_type<T>::value>75 struct tiny_value76 {77 // For float, double, and long double, 1/min_value<T>() is finite.78 // But for mpfr_float and cpp_bin_float, 1/min_value<T>() is inf.79 // Multiply the min by 16 so that the reciprocal doesn't overflow.80 BOOST_MATH_GPU_ENABLED static T get() {81 return 16*tools::min_value<T>();82 }83 };84 template <typename T>85 struct tiny_value<T, true>86 {87 using value_type = typename T::value_type;88 BOOST_MATH_GPU_ENABLED static T get() {89 return 16*tools::min_value<value_type>();90 }91 };92 93} // namespace detail94 95namespace detail {96 97//98// continued_fraction_b99// Evaluates:100//101// b0 + a1102// ---------------103// b1 + a2104// ----------105// b2 + a3106// -----107// b3 + ...108//109// Note that the first a0 returned by generator Gen is discarded.110//111 112template <typename Gen, typename U>113BOOST_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)114 noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) 115 #ifndef BOOST_MATH_HAS_GPU_SUPPORT116 // SYCL can not handle this condition so we only check float on that platform117 && noexcept(std::declval<Gen>()())118 #endif119 )120{121 BOOST_MATH_STD_USING // ADL of std names122 123 using traits = detail::fraction_traits<Gen>;124 using result_type = typename traits::result_type;125 using value_type = typename traits::value_type;126 using integer_type = typename integer_scalar_type<result_type>::type;127 using scalar_type = typename scalar_type<result_type>::type;128 129 integer_type const zero(0), one(1);130 131 result_type tiny = detail::tiny_value<result_type>::get();132 scalar_type terminator = abs(factor);133 134 value_type v = g();135 136 result_type f, C, D, delta;137 f = traits::b(v);138 if(f == zero)139 f = tiny;140 C = f;141 D = 0;142 143 boost::math::uintmax_t counter(max_terms);144 do{145 v = g();146 D = traits::b(v) + traits::a(v) * D;147 if(D == result_type(0))148 D = tiny;149 C = traits::b(v) + traits::a(v) / C;150 if(C == zero)151 C = tiny;152 D = one/D;153 delta = C*D;154 f = f * delta;155 }while((abs(delta - one) > terminator) && --counter);156 157 max_terms = max_terms - counter;158 159 return f;160}161 162} // namespace detail163 164template <typename Gen, typename U>165BOOST_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)166 noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) 167 #ifndef BOOST_MATH_HAS_GPU_SUPPORT168 && noexcept(std::declval<Gen>()())169 #endif170 )171{172 return detail::continued_fraction_b_impl(g, factor, max_terms);173}174 175template <typename Gen, typename U>176BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor)177 noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) 178 #ifndef BOOST_MATH_HAS_GPU_SUPPORT179 && noexcept(std::declval<Gen>()())180 #endif181 )182{183 boost::math::uintmax_t max_terms = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();184 return detail::continued_fraction_b_impl(g, factor, max_terms);185}186 187template <typename Gen>188BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits)189 noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) 190 #ifndef BOOST_MATH_HAS_GPU_SUPPORT191 && noexcept(std::declval<Gen>()())192 #endif193 )194{195 BOOST_MATH_STD_USING // ADL of std names196 197 using traits = detail::fraction_traits<Gen>;198 using result_type = typename traits::result_type;199 200 result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);201 boost::math::uintmax_t max_terms = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();202 return detail::continued_fraction_b_impl(g, factor, max_terms);203}204 205template <typename Gen>206BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits, boost::math::uintmax_t& max_terms)207 noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) 208 #ifndef BOOST_MATH_HAS_GPU_SUPPORT209 && noexcept(std::declval<Gen>()())210 #endif211 )212{213 BOOST_MATH_STD_USING // ADL of std names214 215 using traits = detail::fraction_traits<Gen>;216 using result_type = typename traits::result_type;217 218 result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);219 return detail::continued_fraction_b_impl(g, factor, max_terms);220}221 222namespace detail {223 224//225// continued_fraction_a226// Evaluates:227//228// a1229// ---------------230// b1 + a2231// ----------232// b2 + a3233// -----234// b3 + ...235//236// Note that the first a1 and b1 returned by generator Gen are both used.237//238template <typename Gen, typename U>239BOOST_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)240 noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) 241 #ifndef BOOST_MATH_HAS_GPU_SUPPORT242 && noexcept(std::declval<Gen>()())243 #endif244 )245{246 BOOST_MATH_STD_USING // ADL of std names247 248 using traits = detail::fraction_traits<Gen>;249 using result_type = typename traits::result_type;250 using value_type = typename traits::value_type;251 using integer_type = typename integer_scalar_type<result_type>::type;252 using scalar_type = typename scalar_type<result_type>::type;253 254 integer_type const zero(0), one(1);255 256 result_type tiny = detail::tiny_value<result_type>::get();257 scalar_type terminator = abs(factor);258 259 value_type v = g();260 261 result_type f, C, D, delta, a0;262 f = traits::b(v);263 a0 = traits::a(v);264 if(f == zero)265 f = tiny;266 C = f;267 D = 0;268 269 boost::math::uintmax_t counter(max_terms);270 271 do{272 v = g();273 D = traits::b(v) + traits::a(v) * D;274 if(D == zero)275 D = tiny;276 C = traits::b(v) + traits::a(v) / C;277 if(C == zero)278 C = tiny;279 D = one/D;280 delta = C*D;281 f = f * delta;282 }while((abs(delta - one) > terminator) && --counter);283 284 max_terms = max_terms - counter;285 286 return a0/f;287}288 289} // namespace detail290 291template <typename Gen, typename U>292BOOST_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)293 noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) 294 #ifndef BOOST_MATH_HAS_GPU_SUPPORT295 && noexcept(std::declval<Gen>()())296 #endif297 )298{299 return detail::continued_fraction_a_impl(g, factor, max_terms);300}301 302template <typename Gen, typename U>303BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor)304 noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)305 #ifndef BOOST_MATH_HAS_GPU_SUPPORT306 && noexcept(std::declval<Gen>()())307 #endif308 )309{310 boost::math::uintmax_t max_iter = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();311 return detail::continued_fraction_a_impl(g, factor, max_iter);312}313 314template <typename Gen>315BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits)316 noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) 317 #ifndef BOOST_MATH_HAS_GPU_SUPPORT318 && noexcept(std::declval<Gen>()())319 #endif320 )321{322 BOOST_MATH_STD_USING // ADL of std names323 324 typedef detail::fraction_traits<Gen> traits;325 typedef typename traits::result_type result_type;326 327 result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);328 boost::math::uintmax_t max_iter = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();329 330 return detail::continued_fraction_a_impl(g, factor, max_iter);331}332 333template <typename Gen>334BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits, boost::math::uintmax_t& max_terms)335 noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) 336 #ifndef BOOST_MATH_HAS_GPU_SUPPORT337 && noexcept(std::declval<Gen>()())338 #endif339 )340{341 BOOST_MATH_STD_USING // ADL of std names342 343 using traits = detail::fraction_traits<Gen>;344 using result_type = typename traits::result_type;345 346 result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);347 return detail::continued_fraction_a_impl(g, factor, max_terms);348}349 350} // namespace tools351} // namespace math352} // namespace boost353 354#endif // BOOST_MATH_TOOLS_FRACTION_INCLUDED355