Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:42:37

0001 ///////////////////////////////////////////////////////////////
0002 //  Copyright 2020 John Maddock. Distributed under the Boost
0003 //  Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
0005 
0006 #ifndef BOOST_MP_RATIONAL_ADAPTOR_HPP
0007 #define BOOST_MP_RATIONAL_ADAPTOR_HPP
0008 
0009 #include <boost/multiprecision/number.hpp>
0010 #include <boost/multiprecision/detail/hash.hpp>
0011 #include <boost/multiprecision/detail/float128_functions.hpp>
0012 #include <boost/multiprecision/detail/no_exceptions_support.hpp>
0013 
0014 namespace boost {
0015 namespace multiprecision {
0016 namespace backends {
0017 
0018 template <class Backend>
0019 struct rational_adaptor
0020 {
0021    //
0022    // Each backend need to declare 3 type lists which declare the types
0023    // with which this can interoperate.  These lists must at least contain
0024    // the widest type in each category - so "long long" must be the final
0025    // type in the signed_types list for example.  Any narrower types if not
0026    // present in the list will get promoted to the next wider type that is
0027    // in the list whenever mixed arithmetic involving that type is encountered.
0028    //
0029    typedef typename Backend::signed_types    signed_types;
0030    typedef typename Backend::unsigned_types  unsigned_types;
0031    typedef typename Backend::float_types     float_types;
0032 
0033    typedef typename std::tuple_element<0, unsigned_types>::type ui_type;
0034 
0035    static Backend get_one()
0036    {
0037       Backend t;
0038       t = static_cast<ui_type>(1);
0039       return t;
0040    }
0041    static Backend get_zero()
0042    {
0043       Backend t;
0044       t = static_cast<ui_type>(0);
0045       return t;
0046    }
0047 
0048    static const Backend& one()
0049    {
0050       static const Backend result(get_one());
0051       return result;
0052    }
0053    static const Backend& zero()
0054    {
0055       static const Backend result(get_zero());
0056       return result;
0057    }
0058 
0059    void normalize()
0060    {
0061       using default_ops::eval_gcd;
0062       using default_ops::eval_eq;
0063       using default_ops::eval_divide;
0064       using default_ops::eval_get_sign;
0065 
0066       int s = eval_get_sign(m_denom);
0067 
0068       if(s == 0)
0069       {
0070          BOOST_MP_THROW_EXCEPTION(std::overflow_error("Integer division by zero"));
0071       }
0072       else if (s < 0)
0073       {
0074          m_num.negate();
0075          m_denom.negate();
0076       }
0077 
0078       Backend g, t;
0079       eval_gcd(g, m_num, m_denom);
0080       if (!eval_eq(g, one()))
0081       {
0082          eval_divide(t, m_num, g);
0083          m_num.swap(t);
0084          eval_divide(t, m_denom, g);
0085          m_denom = std::move(t);
0086       }
0087    }
0088 
0089    // We must have a default constructor:
0090    rational_adaptor()
0091       : m_num(zero()), m_denom(one()) {}
0092 
0093    rational_adaptor(const rational_adaptor& o) : m_num(o.m_num), m_denom(o.m_denom) {}
0094    rational_adaptor(rational_adaptor&& o) = default;
0095 
0096    // Optional constructors, we can make this type slightly more efficient
0097    // by providing constructors from any type we can handle natively.
0098    // These will also cause number<> to be implicitly constructible
0099    // from these types unless we make such constructors explicit.
0100    //
0101    template <class Arithmetic>
0102    rational_adaptor(const Arithmetic& val, typename std::enable_if<std::is_constructible<Backend, Arithmetic>::value && !std::is_floating_point<Arithmetic>::value>::type const* = nullptr)
0103       : m_num(val), m_denom(one()) {}
0104 
0105    //
0106    // Pass-through 2-arg construction of components:
0107    //
0108    template <class T, class U>
0109    rational_adaptor(const T& a, const U& b, typename std::enable_if<std::is_constructible<Backend, T const&>::value && std::is_constructible<Backend, U const&>::value>::type const* = nullptr)
0110       : m_num(a), m_denom(b) 
0111    {
0112       normalize();
0113    }
0114    template <class T, class U>
0115    rational_adaptor(T&& a, const U& b, typename std::enable_if<std::is_constructible<Backend, T>::value && std::is_constructible<Backend, U>::value>::type const* = nullptr)
0116       : m_num(static_cast<T&&>(a)), m_denom(b) 
0117    {
0118       normalize();
0119    }
0120    template <class T, class U>
0121    rational_adaptor(T&& a, U&& b, typename std::enable_if<std::is_constructible<Backend, T>::value && std::is_constructible<Backend, U>::value>::type const* = nullptr)
0122       : m_num(static_cast<T&&>(a)), m_denom(static_cast<U&&>(b)) 
0123    {
0124       normalize();
0125    }
0126    template <class T, class U>
0127    rational_adaptor(const T& a, U&& b, typename std::enable_if<std::is_constructible<Backend, T>::value && std::is_constructible<Backend, U>::value>::type const* = nullptr)
0128       : m_num(a), m_denom(static_cast<U&&>(b)) 
0129    {
0130       normalize();
0131    }
0132    //
0133    // In the absense of converting constructors, operator= takes the strain.
0134    // In addition to the usual suspects, there must be one operator= for each type
0135    // listed in signed_types, unsigned_types, and float_types plus a string constructor.
0136    //
0137    rational_adaptor& operator=(const rational_adaptor& o) = default;
0138    rational_adaptor& operator=(rational_adaptor&& o) = default;
0139    template <class Arithmetic>
0140    inline typename std::enable_if<!std::is_floating_point<Arithmetic>::value, rational_adaptor&>::type operator=(const Arithmetic& i)
0141    {
0142       m_num = i;
0143       m_denom = one();
0144       return *this;
0145    }
0146    rational_adaptor& operator=(const char* s)
0147    {
0148       using default_ops::eval_eq;
0149 
0150       std::string                        s1;
0151       multiprecision::number<Backend>    v1, v2;
0152       char                               c;
0153       bool                               have_hex = false;
0154       const char* p = s; // saved for later
0155 
0156       while ((0 != (c = *s)) && (c == 'x' || c == 'X' || c == '-' || c == '+' || (c >= '0' && c <= '9') || (have_hex && (c >= 'a' && c <= 'f')) || (have_hex && (c >= 'A' && c <= 'F'))))
0157       {
0158          if (c == 'x' || c == 'X')
0159             have_hex = true;
0160          s1.append(1, c);
0161          ++s;
0162       }
0163       v1.assign(s1);
0164       s1.erase();
0165       if (c == '/')
0166       {
0167          ++s;
0168          while ((0 != (c = *s)) && (c == 'x' || c == 'X' || c == '-' || c == '+' || (c >= '0' && c <= '9') || (have_hex && (c >= 'a' && c <= 'f')) || (have_hex && (c >= 'A' && c <= 'F'))))
0169          {
0170             if (c == 'x' || c == 'X')
0171                have_hex = true;
0172             s1.append(1, c);
0173             ++s;
0174          }
0175          v2.assign(s1);
0176       }
0177       else
0178          v2 = 1;
0179       if (*s)
0180       {
0181          BOOST_MP_THROW_EXCEPTION(std::runtime_error(std::string("Could not parse the string \"") + p + std::string("\" as a valid rational number.")));
0182       }
0183       multiprecision::number<Backend> gcd;
0184       eval_gcd(gcd.backend(), v1.backend(), v2.backend());
0185       if (!eval_eq(gcd.backend(), one()))
0186       {
0187          v1 /= gcd;
0188          v2 /= gcd;
0189       }
0190       num() = std::move(std::move(v1).backend());
0191       denom() = std::move(std::move(v2).backend());
0192       return *this;
0193    }
0194    template <class Float>
0195    typename std::enable_if<std::is_floating_point<Float>::value, rational_adaptor&>::type operator=(Float i)
0196    {
0197       using default_ops::eval_eq;
0198       BOOST_MP_FLOAT128_USING using std::floor; using std::frexp; using std::ldexp;
0199 
0200       int   e;
0201       Float f = frexp(i, &e);
0202 #ifdef BOOST_HAS_FLOAT128
0203       f = ldexp(f, std::is_same<float128_type, Float>::value ? 113 : std::numeric_limits<Float>::digits);
0204       e -= std::is_same<float128_type, Float>::value ? 113 : std::numeric_limits<Float>::digits;
0205 #else
0206       f = ldexp(f, std::numeric_limits<Float>::digits);
0207       e -= std::numeric_limits<Float>::digits;
0208 #endif
0209       number<Backend> num(f);
0210       number<Backend> denom(1u);
0211       if (e > 0)
0212       {
0213          num <<= e;
0214       }
0215       else if (e < 0)
0216       {
0217          denom <<= -e;
0218       }
0219       number<Backend> gcd;
0220       eval_gcd(gcd.backend(), num.backend(), denom.backend());
0221       if (!eval_eq(gcd.backend(), one()))
0222       {
0223          num /= gcd;
0224          denom /= gcd;
0225       }
0226       this->num() = std::move(std::move(num).backend());
0227       this->denom() = std::move(std::move(denom).backend());
0228       return *this;
0229    }
0230 
0231    void swap(rational_adaptor& o)
0232    {
0233       m_num.swap(o.m_num);
0234       m_denom.swap(o.m_denom);
0235    }
0236    std::string str(std::streamsize digits, std::ios_base::fmtflags f) const
0237    {
0238       using default_ops::eval_eq;
0239       //
0240       // We format the string ourselves so we can match what GMP's mpq type does:
0241       //
0242       std::string result = num().str(digits, f);
0243       if (!eval_eq(denom(), one()))
0244       {
0245          result.append(1, '/');
0246          result.append(denom().str(digits, f));
0247       }
0248       return result;
0249    }
0250    void negate()
0251    {
0252       m_num.negate();
0253    }
0254    int compare(const rational_adaptor& o) const
0255    {
0256       std::ptrdiff_t s1 = eval_get_sign(*this);
0257       std::ptrdiff_t s2 = eval_get_sign(o);
0258       if (s1 != s2)
0259       {
0260          return s1 < s2 ? -1 : 1;
0261       }
0262       else if (s1 == 0)
0263          return 0; // both zero.
0264 
0265       bool neg = false;
0266       if (s1 >= 0)
0267       {
0268          s1 = eval_msb(num()) + eval_msb(o.denom());
0269          s2 = eval_msb(o.num()) + eval_msb(denom());
0270       }
0271       else
0272       {
0273          Backend t(num());
0274          t.negate();
0275          s1 = eval_msb(t) + eval_msb(o.denom());
0276          t = o.num();
0277          t.negate();
0278          s2 = eval_msb(t) + eval_msb(denom());
0279          neg = true;
0280       }
0281       s1 -= s2;
0282       if (s1 < -1)
0283          return neg ? 1 : -1;
0284       else if (s1 > 1)
0285          return neg ? -1 : 1;
0286 
0287       Backend t1, t2;
0288       eval_multiply(t1, num(), o.denom());
0289       eval_multiply(t2, o.num(), denom());
0290       return t1.compare(t2);
0291    }
0292    //
0293    // Comparison with arithmetic types, default just constructs a temporary:
0294    //
0295    template <class A>
0296    typename std::enable_if<boost::multiprecision::detail::is_arithmetic<A>::value, int>::type compare(A i) const
0297    {
0298       rational_adaptor t;
0299       t = i;  //  Note: construct directly from i if supported.
0300       return compare(t);
0301    }
0302 
0303    Backend& num() { return m_num; }
0304    const Backend& num()const { return m_num; }
0305    Backend& denom() { return m_denom; }
0306    const Backend& denom()const { return m_denom; }
0307 
0308    #ifndef BOOST_MP_STANDALONE
0309    template <class Archive>
0310    void serialize(Archive& ar, const std::integral_constant<bool, true>&)
0311    {
0312       // Saving
0313       number<Backend> n(num()), d(denom());
0314       ar& boost::make_nvp("numerator", n);
0315       ar& boost::make_nvp("denominator", d);
0316    }
0317    template <class Archive>
0318    void serialize(Archive& ar, const std::integral_constant<bool, false>&)
0319    {
0320       // Loading
0321       number<Backend> n, d;
0322       ar& boost::make_nvp("numerator", n);
0323       ar& boost::make_nvp("denominator", d);
0324       num() = n.backend();
0325       denom() = d.backend();
0326    }
0327    template <class Archive>
0328    void serialize(Archive& ar, const unsigned int /*version*/)
0329    {
0330       using tag = typename Archive::is_saving;
0331       using saving_tag = std::integral_constant<bool, tag::value>;
0332       serialize(ar, saving_tag());
0333    }
0334    #endif // BOOST_MP_STANDALONE
0335    
0336  private:
0337    Backend m_num, m_denom;
0338 };
0339 
0340 //
0341 // Helpers:
0342 //
0343 template <class T>
0344 inline constexpr typename std::enable_if<std::numeric_limits<T>::is_specialized && !std::numeric_limits<T>::is_signed, bool>::type
0345 is_minus_one(const T&)
0346 {
0347    return false;
0348 }
0349 template <class T>
0350 inline constexpr typename std::enable_if<!std::numeric_limits<T>::is_specialized || std::numeric_limits<T>::is_signed, bool>::type
0351 is_minus_one(const T& val)
0352 {
0353    return val == -1;
0354 }
0355 
0356 //
0357 // Required non-members:
0358 //
0359 template <class Backend> 
0360 inline void eval_add(rational_adaptor<Backend>& a, const rational_adaptor<Backend>& b)
0361 {
0362    eval_add_subtract_imp(a, a, b, true);
0363 }
0364 template <class Backend> 
0365 inline void eval_subtract(rational_adaptor<Backend>& a, const rational_adaptor<Backend>& b)
0366 {
0367    eval_add_subtract_imp(a, a, b, false);
0368 }
0369 
0370 template <class Backend> 
0371 inline void eval_multiply(rational_adaptor<Backend>& a, const rational_adaptor<Backend>& b)
0372 {
0373    eval_multiply_imp(a, a, b.num(), b.denom());
0374 }
0375 
0376 template <class Backend> 
0377 void eval_divide(rational_adaptor<Backend>& a, const rational_adaptor<Backend>& b)
0378 {
0379    using default_ops::eval_divide;
0380    rational_adaptor<Backend> t;
0381    eval_divide(t, a, b);
0382    a = std::move(t);
0383 }
0384 //
0385 // Conversions:
0386 //
0387 template <class R, class IntBackend>
0388 inline typename std::enable_if<number_category<R>::value == number_kind_floating_point>::type eval_convert_to(R* result, const rational_adaptor<IntBackend>& backend)
0389 {
0390    //
0391    // The generic conversion is as good as anything we can write here:
0392    //
0393    ::boost::multiprecision::detail::generic_convert_rational_to_float(*result, backend);
0394 }
0395 
0396 template <class R, class IntBackend>
0397 inline typename std::enable_if<(number_category<R>::value != number_kind_integer) && (number_category<R>::value != number_kind_floating_point) && !std::is_enum<R>::value>::type eval_convert_to(R* result, const rational_adaptor<IntBackend>& backend)
0398 {
0399    using default_ops::eval_convert_to;
0400    R d;
0401    eval_convert_to(result, backend.num());
0402    eval_convert_to(&d, backend.denom());
0403    *result /= d;
0404 }
0405 
0406 template <class R, class Backend>
0407 inline typename std::enable_if<number_category<R>::value == number_kind_integer>::type eval_convert_to(R* result, const rational_adaptor<Backend>& backend)
0408 {
0409    using default_ops::eval_divide;
0410    using default_ops::eval_convert_to;
0411    Backend t;
0412    eval_divide(t, backend.num(), backend.denom());
0413    eval_convert_to(result, t);
0414 }
0415 
0416 //
0417 // Hashing support, not strictly required, but it is used in our tests:
0418 //
0419 template <class Backend>
0420 inline std::size_t hash_value(const rational_adaptor<Backend>& arg)
0421 {
0422    std::size_t result = hash_value(arg.num());
0423    std::size_t result2 = hash_value(arg.denom());
0424    boost::multiprecision::detail::hash_combine(result, result2);
0425    return result;
0426 }
0427 //
0428 // assign_components:
0429 //
0430 template <class Backend>
0431 void assign_components(rational_adaptor<Backend>& result, Backend const& a, Backend const& b)
0432 {
0433    using default_ops::eval_gcd;
0434    using default_ops::eval_divide;
0435    using default_ops::eval_eq;
0436    using default_ops::eval_is_zero;
0437    using default_ops::eval_get_sign;
0438 
0439    if (eval_is_zero(b))
0440    {
0441       BOOST_MP_THROW_EXCEPTION(std::overflow_error("Integer division by zero"));
0442    }
0443    Backend g;
0444    eval_gcd(g, a, b);
0445    if (eval_eq(g, rational_adaptor<Backend>::one()))
0446    {
0447       result.num() = a;
0448       result.denom() = b;
0449    }
0450    else
0451    {
0452       eval_divide(result.num(), a, g);
0453       eval_divide(result.denom(), b, g);
0454    }
0455    if (eval_get_sign(result.denom()) < 0)
0456    {
0457       result.num().negate();
0458       result.denom().negate();
0459    }
0460 }
0461 //
0462 // Again for arithmetic types, overload for whatever arithmetic types are directly supported:
0463 //
0464 template <class Backend, class Arithmetic1, class Arithmetic2>
0465 inline void assign_components(rational_adaptor<Backend>& result, const Arithmetic1& a, typename std::enable_if<std::is_arithmetic<Arithmetic1>::value && std::is_arithmetic<Arithmetic2>::value, const Arithmetic2&>::type b)
0466 {
0467    using default_ops::eval_gcd;
0468    using default_ops::eval_divide;
0469    using default_ops::eval_eq;
0470 
0471    if (b == 0)
0472    {
0473       BOOST_MP_THROW_EXCEPTION(std::overflow_error("Integer division by zero"));
0474    }
0475 
0476    Backend g;
0477    result.num()   = a;
0478    eval_gcd(g, result.num(), b);
0479    if (eval_eq(g, rational_adaptor<Backend>::one()))
0480    {
0481       result.denom() = b;
0482    }
0483    else
0484    {
0485       eval_divide(result.num(), g);
0486       eval_divide(result.denom(), b, g);
0487    }
0488    if (eval_get_sign(result.denom()) < 0)
0489    {
0490       result.num().negate();
0491       result.denom().negate();
0492    }
0493 }
0494 template <class Backend, class Arithmetic1, class Arithmetic2>
0495 inline void assign_components(rational_adaptor<Backend>& result, const Arithmetic1& a, typename std::enable_if<!std::is_arithmetic<Arithmetic1>::value || !std::is_arithmetic<Arithmetic2>::value, const Arithmetic2&>::type b)
0496 {
0497    using default_ops::eval_gcd;
0498    using default_ops::eval_divide;
0499    using default_ops::eval_eq;
0500 
0501    Backend g;
0502    result.num()   = a;
0503    result.denom() = b;
0504 
0505    if (eval_get_sign(result.denom()) == 0)
0506    {
0507       BOOST_MP_THROW_EXCEPTION(std::overflow_error("Integer division by zero"));
0508    }
0509 
0510    eval_gcd(g, result.num(), result.denom());
0511    if (!eval_eq(g, rational_adaptor<Backend>::one()))
0512    {
0513       eval_divide(result.num(), g);
0514       eval_divide(result.denom(), g);
0515    }
0516    if (eval_get_sign(result.denom()) < 0)
0517    {
0518       result.num().negate();
0519       result.denom().negate();
0520    }
0521 }
0522 //
0523 // Optional comparison operators:
0524 //
0525 template <class Backend>
0526 inline bool eval_is_zero(const rational_adaptor<Backend>& arg)
0527 {
0528    using default_ops::eval_is_zero;
0529    return eval_is_zero(arg.num());
0530 }
0531 
0532 template <class Backend>
0533 inline int eval_get_sign(const rational_adaptor<Backend>& arg)
0534 {
0535    using default_ops::eval_get_sign;
0536    return eval_get_sign(arg.num());
0537 }
0538 
0539 template <class Backend>
0540 inline bool eval_eq(const rational_adaptor<Backend>& a, const rational_adaptor<Backend>& b)
0541 {
0542    using default_ops::eval_eq;
0543    return eval_eq(a.num(), b.num()) && eval_eq(a.denom(), b.denom());
0544 }
0545 
0546 template <class Backend, class Arithmetic>
0547 inline typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value&& std::is_integral<Arithmetic>::value, bool>::type 
0548    eval_eq(const rational_adaptor<Backend>& a, Arithmetic b)
0549 {
0550    using default_ops::eval_eq;
0551    return eval_eq(a.denom(), rational_adaptor<Backend>::one()) && eval_eq(a.num(), b);
0552 }
0553 
0554 template <class Backend, class Arithmetic>
0555 inline typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value&& std::is_integral<Arithmetic>::value, bool>::type 
0556    eval_eq(Arithmetic b, const rational_adaptor<Backend>& a)
0557 {
0558    using default_ops::eval_eq;
0559    return eval_eq(a.denom(), rational_adaptor<Backend>::one()) && eval_eq(a.num(), b);
0560 }
0561 
0562 //
0563 // Arithmetic operations, starting with addition:
0564 //
0565 template <class Backend, class Arithmetic> 
0566 void eval_add_subtract_imp(rational_adaptor<Backend>& result, const Arithmetic& arg, bool isaddition)
0567 {
0568    using default_ops::eval_multiply;
0569    using default_ops::eval_divide;
0570    using default_ops::eval_add;
0571    using default_ops::eval_gcd;
0572    Backend t;
0573    eval_multiply(t, result.denom(), arg);
0574    if (isaddition)
0575       eval_add(result.num(), t);
0576    else
0577       eval_subtract(result.num(), t);
0578    //
0579    // There is no need to re-normalize here, we have 
0580    // (a + bm) / b
0581    // and gcd(a + bm, b) = gcd(a, b) = 1
0582    //
0583    /*
0584    eval_gcd(t, result.num(), result.denom());
0585    if (!eval_eq(t, rational_adaptor<Backend>::one()) != 0)
0586    {
0587       Backend t2;
0588       eval_divide(t2, result.num(), t);
0589       t2.swap(result.num());
0590       eval_divide(t2, result.denom(), t);
0591       t2.swap(result.denom());
0592    }
0593    */
0594 }
0595 
0596 template <class Backend, class Arithmetic> 
0597 inline typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && (std::is_integral<Arithmetic>::value || std::is_same<Arithmetic, Backend>::value)>::type
0598    eval_add(rational_adaptor<Backend>& result, const Arithmetic& arg)
0599 {
0600    eval_add_subtract_imp(result, arg, true);
0601 }
0602 
0603 template <class Backend, class Arithmetic> 
0604 inline typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && (std::is_integral<Arithmetic>::value || std::is_same<Arithmetic, Backend>::value)>::type
0605    eval_subtract(rational_adaptor<Backend>& result, const Arithmetic& arg)
0606 {
0607    eval_add_subtract_imp(result, arg, false);
0608 }
0609 
0610 template <class Backend>
0611 void eval_add_subtract_imp(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, const rational_adaptor<Backend>& b, bool isaddition)
0612 {
0613    using default_ops::eval_eq;
0614    using default_ops::eval_multiply;
0615    using default_ops::eval_divide;
0616    using default_ops::eval_add;
0617    using default_ops::eval_subtract;
0618    //
0619    // Let  a = an/ad
0620    //      b = bn/bd
0621    //      g = gcd(ad, bd)
0622    // result = rn/rd
0623    //
0624    // Then:
0625    // rn = an * (bd/g) + bn * (ad/g)
0626    // rd = ad * (bd/g)
0627    //    = (ad/g) * (bd/g) * g
0628    //
0629    // And the whole thing can then be rescaled by
0630    //      gcd(rn, g)
0631    //
0632    Backend gcd, t1, t2, t3, t4;
0633    //
0634    // Begin by getting the gcd of the 2 denominators:
0635    //
0636    eval_gcd(gcd, a.denom(), b.denom());
0637    //
0638    // Do we have gcd > 1:
0639    //
0640    if (!eval_eq(gcd, rational_adaptor<Backend>::one()))
0641    {
0642       //
0643       // Scale the denominators by gcd, and put the results in t1 and t2:
0644       //
0645       eval_divide(t1, b.denom(), gcd);
0646       eval_divide(t2, a.denom(), gcd);
0647       //
0648       // multiply the numerators by the scale denominators and put the results in t3, t4:
0649       //
0650       eval_multiply(t3, a.num(), t1);
0651       eval_multiply(t4, b.num(), t2);
0652       //
0653       // Add them up:
0654       //
0655       if (isaddition)
0656          eval_add(t3, t4);
0657       else
0658          eval_subtract(t3, t4);
0659       //
0660       // Get the gcd of gcd and our numerator (t3):
0661       //
0662       eval_gcd(t4, t3, gcd);
0663       if (eval_eq(t4, rational_adaptor<Backend>::one()))
0664       {
0665          result.num() = t3;
0666          eval_multiply(result.denom(), t1, a.denom());
0667       }
0668       else
0669       {
0670          //
0671          // Uncommon case where gcd is not 1, divide the numerator
0672          // and the denominator terms by the new gcd.  Note we perform division
0673          // on the existing gcd value as this is the smallest of the 3 denominator
0674          // terms we'll be multiplying together, so there's a good chance it's a
0675          // single limb value already:
0676          //
0677          eval_divide(result.num(), t3, t4);
0678          eval_divide(t3, gcd, t4);
0679          eval_multiply(t4, t1, t2);
0680          eval_multiply(result.denom(), t4, t3);
0681       }
0682    }
0683    else
0684    {
0685       //
0686       // Most common case (approx 60%) where gcd is one:
0687       //
0688       eval_multiply(t1, a.num(), b.denom());
0689       eval_multiply(t2, a.denom(), b.num());
0690       if (isaddition)
0691          eval_add(result.num(), t1, t2);
0692       else
0693          eval_subtract(result.num(), t1, t2);
0694       eval_multiply(result.denom(), a.denom(), b.denom());
0695    }
0696 }
0697 
0698 
0699 template <class Backend>
0700 inline void eval_add(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, const rational_adaptor<Backend>& b)
0701 {
0702    eval_add_subtract_imp(result, a, b, true);
0703 }
0704 template <class Backend>
0705 inline void eval_subtract(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, const rational_adaptor<Backend>& b)
0706 {
0707    eval_add_subtract_imp(result, a, b, false);
0708 }
0709 
0710 template <class Backend, class Arithmetic>
0711 void eval_add_subtract_imp(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, const Arithmetic& b, bool isaddition)
0712 {
0713    using default_ops::eval_add;
0714    using default_ops::eval_subtract;
0715    using default_ops::eval_multiply;
0716 
0717    if (&result == &a)
0718       return eval_add_subtract_imp(result, b, isaddition);
0719 
0720    eval_multiply(result.num(), a.denom(), b);
0721    if (isaddition)
0722       eval_add(result.num(), a.num());
0723    else
0724       BOOST_IF_CONSTEXPR(std::numeric_limits<Backend>::is_signed == false)
0725    {
0726       Backend t;
0727       eval_subtract(t, a.num(), result.num());
0728       result.num() = std::move(t);
0729    }
0730    else
0731    {
0732       eval_subtract(result.num(), a.num());
0733       result.negate();
0734    }
0735    result.denom() = a.denom();
0736    //
0737    // There is no need to re-normalize here, we have 
0738    // (a + bm) / b
0739    // and gcd(a + bm, b) = gcd(a, b) = 1
0740    //
0741 }
0742 template <class Backend, class Arithmetic>
0743 inline typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && (std::is_integral<Arithmetic>::value || std::is_same<Arithmetic, Backend>::value)>::type
0744    eval_add(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, const Arithmetic& b)
0745 {
0746    eval_add_subtract_imp(result, a, b, true);
0747 }
0748 template <class Backend, class Arithmetic>
0749 inline typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && (std::is_integral<Arithmetic>::value || std::is_same<Arithmetic, Backend>::value)>::type
0750    eval_subtract(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, const Arithmetic& b)
0751 {
0752    eval_add_subtract_imp(result, a, b, false);
0753 }
0754 
0755 //
0756 // Multiplication:
0757 //
0758 template <class Backend> 
0759 void eval_multiply_imp(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, const Backend& b_num, const Backend& b_denom)
0760 {
0761    using default_ops::eval_multiply;
0762    using default_ops::eval_divide;
0763    using default_ops::eval_gcd;
0764    using default_ops::eval_get_sign;
0765    using default_ops::eval_eq;
0766 
0767    Backend gcd_left, gcd_right, t1, t2;
0768    eval_gcd(gcd_left, a.num(), b_denom);
0769    eval_gcd(gcd_right, b_num, a.denom());
0770    //
0771    // Unit gcd's are the most likely case:
0772    //
0773    bool b_left = eval_eq(gcd_left, rational_adaptor<Backend>::one());
0774    bool b_right = eval_eq(gcd_right, rational_adaptor<Backend>::one());
0775 
0776    if (b_left && b_right)
0777    {
0778       eval_multiply(result.num(), a.num(), b_num);
0779       eval_multiply(result.denom(), a.denom(), b_denom);
0780    }
0781    else if (b_left)
0782    {
0783       eval_divide(t2, b_num, gcd_right);
0784       eval_multiply(result.num(), a.num(), t2);
0785       eval_divide(t1, a.denom(), gcd_right);
0786       eval_multiply(result.denom(), t1, b_denom);
0787    }
0788    else if (b_right)
0789    {
0790       eval_divide(t1, a.num(), gcd_left);
0791       eval_multiply(result.num(), t1, b_num);
0792       eval_divide(t2, b_denom, gcd_left);
0793       eval_multiply(result.denom(), a.denom(), t2);
0794    }
0795    else
0796    {
0797       eval_divide(t1, a.num(), gcd_left);
0798       eval_divide(t2, b_num, gcd_right);
0799       eval_multiply(result.num(), t1, t2);
0800       eval_divide(t1, a.denom(), gcd_right);
0801       eval_divide(t2, b_denom, gcd_left);
0802       eval_multiply(result.denom(), t1, t2);
0803    }
0804    //
0805    // We may have b_denom negative if this is actually division, if so just correct things now:
0806    //
0807    if (eval_get_sign(b_denom) < 0)
0808    {
0809       result.num().negate();
0810       result.denom().negate();
0811    }
0812 }
0813 
0814 template <class Backend> 
0815 void eval_multiply(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, const rational_adaptor<Backend>& b)
0816 {
0817    using default_ops::eval_multiply;
0818 
0819    if (&a == &b)
0820    {
0821       // squaring, gcd's are 1:
0822       eval_multiply(result.num(), a.num(), b.num());
0823       eval_multiply(result.denom(), a.denom(), b.denom());
0824       return;
0825    }
0826    eval_multiply_imp(result, a, b.num(), b.denom());
0827 }
0828 
0829 template <class Backend, class Arithmetic> 
0830 void eval_multiply_imp(Backend& result_num, Backend& result_denom, Arithmetic arg)
0831 {
0832    if (arg == 0)
0833    {
0834       result_num = rational_adaptor<Backend>::zero();
0835       result_denom = rational_adaptor<Backend>::one();
0836       return;
0837    }
0838    else if (arg == 1)
0839       return;
0840 
0841    using default_ops::eval_multiply;
0842    using default_ops::eval_divide;
0843    using default_ops::eval_gcd;
0844    using default_ops::eval_convert_to;
0845 
0846    Backend gcd, t;
0847    Arithmetic integer_gcd;
0848    eval_gcd(gcd, result_denom, arg);
0849    eval_convert_to(&integer_gcd, gcd);
0850    arg /= integer_gcd;
0851    if (boost::multiprecision::detail::unsigned_abs(arg) > 1)
0852    {
0853       eval_multiply(t, result_num, arg);
0854       result_num = std::move(t);
0855    }
0856    else if (is_minus_one(arg))
0857       result_num.negate();
0858    if (integer_gcd > 1)
0859    {
0860       eval_divide(t, result_denom, integer_gcd);
0861       result_denom = std::move(t);
0862    }
0863 }
0864 template <class Backend> 
0865 void eval_multiply_imp(Backend& result_num, Backend& result_denom, Backend arg)
0866 {
0867    using default_ops::eval_multiply;
0868    using default_ops::eval_divide;
0869    using default_ops::eval_gcd;
0870    using default_ops::eval_convert_to;
0871    using default_ops::eval_is_zero;
0872    using default_ops::eval_eq;
0873    using default_ops::eval_get_sign;
0874 
0875    if (eval_is_zero(arg))
0876    {
0877       result_num = rational_adaptor<Backend>::zero();
0878       result_denom = rational_adaptor<Backend>::one();
0879       return;
0880    }
0881    else if (eval_eq(arg, rational_adaptor<Backend>::one()))
0882       return;
0883 
0884    Backend gcd, t;
0885    eval_gcd(gcd, result_denom, arg);
0886    if (!eval_eq(gcd, rational_adaptor<Backend>::one()))
0887    {
0888       eval_divide(t, arg, gcd);
0889       arg = t;
0890    }
0891    else
0892       t = arg;
0893    if (eval_get_sign(arg) < 0)
0894       t.negate();
0895 
0896    if (!eval_eq(t, rational_adaptor<Backend>::one()))
0897    {
0898       eval_multiply(t, result_num, arg);
0899       result_num = std::move(t);
0900    }
0901    else if (eval_get_sign(arg) < 0)
0902       result_num.negate();
0903    if (!eval_eq(gcd, rational_adaptor<Backend>::one()))
0904    {
0905       eval_divide(t, result_denom, gcd);
0906       result_denom = std::move(t);
0907    }
0908 }
0909 
0910 template <class Backend, class Arithmetic> 
0911 inline typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && (std::is_integral<Arithmetic>::value || std::is_same<Arithmetic, Backend>::value)>::type
0912    eval_multiply(rational_adaptor<Backend>& result, const Arithmetic& arg)
0913 {
0914    eval_multiply_imp(result.num(), result.denom(), arg);
0915 }
0916 
0917 template <class Backend, class Arithmetic> 
0918 typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && std::is_integral<Arithmetic>::value>::type
0919    eval_multiply_imp(rational_adaptor<Backend>& result, const Backend& a_num, const Backend& a_denom, Arithmetic b)
0920 {
0921    if (b == 0)
0922    {
0923       result.num() = rational_adaptor<Backend>::zero();
0924       result.denom() = rational_adaptor<Backend>::one();
0925       return;
0926    }
0927    else if (b == 1)
0928    {
0929       result.num() = a_num;
0930       result.denom() = a_denom;
0931       return;
0932    }
0933 
0934    using default_ops::eval_multiply;
0935    using default_ops::eval_divide;
0936    using default_ops::eval_gcd;
0937    using default_ops::eval_convert_to;
0938 
0939    Backend gcd;
0940    Arithmetic integer_gcd;
0941    eval_gcd(gcd, a_denom, b);
0942    eval_convert_to(&integer_gcd, gcd);
0943    b /= integer_gcd;
0944    if (boost::multiprecision::detail::unsigned_abs(b) > 1)
0945       eval_multiply(result.num(), a_num, b);
0946    else if (is_minus_one(b))
0947    {
0948       result.num() = a_num;
0949       result.num().negate();
0950    }
0951    else
0952       result.num() = a_num;
0953    if (integer_gcd > 1)
0954       eval_divide(result.denom(), a_denom, integer_gcd);
0955    else
0956       result.denom() = a_denom;
0957 }
0958 template <class Backend> 
0959 inline void eval_multiply_imp(rational_adaptor<Backend>& result, const Backend& a_num, const Backend& a_denom, const Backend& b)
0960 {
0961    result.num() = a_num;
0962    result.denom() = a_denom;
0963    eval_multiply_imp(result.num(), result.denom(), b);
0964 }
0965 
0966 template <class Backend, class Arithmetic> 
0967 inline typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && (std::is_integral<Arithmetic>::value || std::is_same<Arithmetic, Backend>::value)>::type
0968    eval_multiply(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, const Arithmetic& b)
0969 {
0970    if (&result == &a)
0971       return eval_multiply(result, b);
0972 
0973    eval_multiply_imp(result, a.num(), a.denom(), b);
0974 }
0975 
0976 template <class Backend, class Arithmetic> 
0977 inline typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && (std::is_integral<Arithmetic>::value || std::is_same<Arithmetic, Backend>::value)>::type
0978    eval_multiply(rational_adaptor<Backend>& result, const Arithmetic& b, const rational_adaptor<Backend>& a)
0979 {
0980    return eval_multiply(result, a, b);
0981 }
0982 
0983 //
0984 // Division:
0985 //
0986 template <class Backend>
0987 inline void eval_divide(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, const rational_adaptor<Backend>& b)
0988 {
0989    using default_ops::eval_multiply;
0990    using default_ops::eval_get_sign;
0991 
0992    if (eval_get_sign(b.num()) == 0)
0993    {
0994       BOOST_MP_THROW_EXCEPTION(std::overflow_error("Integer division by zero"));
0995       return;
0996    }
0997    if (&a == &b)
0998    {
0999       // Huh? Really?
1000       result.num() = result.denom() = rational_adaptor<Backend>::one();
1001       return;
1002    }
1003    if (&result == &b)
1004    {
1005       rational_adaptor<Backend> t(b);
1006       return eval_divide(result, a, t);
1007    }
1008    eval_multiply_imp(result, a, b.denom(), b.num());
1009 }
1010 
1011 template <class Backend, class Arithmetic> 
1012 inline typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && (std::is_integral<Arithmetic>::value || std::is_same<Arithmetic, Backend>::value)>::type
1013    eval_divide(rational_adaptor<Backend>& result, const Arithmetic& b, const rational_adaptor<Backend>& a)
1014 {
1015    using default_ops::eval_get_sign;
1016 
1017    if (eval_get_sign(a.num()) == 0)
1018    {
1019       BOOST_MP_THROW_EXCEPTION(std::overflow_error("Integer division by zero"));
1020       return;
1021    }
1022    if (&a == &result)
1023    {
1024       eval_multiply_imp(result.denom(), result.num(), b);
1025       result.num().swap(result.denom());
1026    }
1027    else
1028       eval_multiply_imp(result, a.denom(), a.num(), b);
1029 
1030    if (eval_get_sign(result.denom()) < 0)
1031    {
1032       result.num().negate();
1033       result.denom().negate();
1034    }
1035 }
1036 
1037 template <class Backend, class Arithmetic>
1038 typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && std::is_integral<Arithmetic>::value>::type
1039 eval_divide(rational_adaptor<Backend>& result, Arithmetic arg)
1040 {
1041    if (arg == 0)
1042    {
1043       BOOST_MP_THROW_EXCEPTION(std::overflow_error("Integer division by zero"));
1044       return;
1045    }
1046    else if (arg == 1)
1047       return;
1048    else if (is_minus_one(arg))
1049    {
1050       result.negate();
1051       return;
1052    }
1053    if (eval_get_sign(result) == 0)
1054    {
1055       return;
1056    }
1057 
1058 
1059    using default_ops::eval_multiply;
1060    using default_ops::eval_gcd;
1061    using default_ops::eval_convert_to;
1062    using default_ops::eval_divide;
1063 
1064    Backend gcd, t;
1065    Arithmetic integer_gcd;
1066    eval_gcd(gcd, result.num(), arg);
1067    eval_convert_to(&integer_gcd, gcd);
1068    arg /= integer_gcd;
1069 
1070    eval_multiply(t, result.denom(), boost::multiprecision::detail::unsigned_abs(arg));
1071    result.denom() = std::move(t);
1072    if (arg < 0)
1073    {
1074       result.num().negate();
1075    }
1076    if (integer_gcd > 1)
1077    {
1078       eval_divide(t, result.num(), integer_gcd);
1079       result.num() = std::move(t);
1080    }
1081 }
1082 template <class Backend>
1083 void eval_divide(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, Backend arg)
1084 {
1085    using default_ops::eval_multiply;
1086    using default_ops::eval_gcd;
1087    using default_ops::eval_convert_to;
1088    using default_ops::eval_divide;
1089    using default_ops::eval_is_zero;
1090    using default_ops::eval_eq;
1091    using default_ops::eval_get_sign;
1092 
1093    if (eval_is_zero(arg))
1094    {
1095       BOOST_MP_THROW_EXCEPTION(std::overflow_error("Integer division by zero"));
1096       return;
1097    }
1098    else if (eval_eq(a, rational_adaptor<Backend>::one()) || (eval_get_sign(a) == 0))
1099    {
1100       if (&result != &a)
1101          result = a;
1102       result.denom() = arg;
1103       return;
1104    }
1105 
1106    Backend gcd, u_arg, t;
1107    eval_gcd(gcd, a.num(), arg);
1108    bool has_unit_gcd = eval_eq(gcd, rational_adaptor<Backend>::one());
1109    if (!has_unit_gcd)
1110    {
1111       eval_divide(u_arg, arg, gcd);
1112       arg = u_arg;
1113    }
1114    else
1115       u_arg = arg;
1116    if (eval_get_sign(u_arg) < 0)
1117       u_arg.negate();
1118 
1119    eval_multiply(t, a.denom(), u_arg);
1120    result.denom() = std::move(t);
1121    
1122    if (!has_unit_gcd)
1123    {
1124       eval_divide(t, a.num(), gcd);
1125       result.num() = std::move(t);
1126    }
1127    else if (&result != &a)
1128       result.num() = a.num();
1129 
1130    if (eval_get_sign(arg) < 0)
1131    {
1132       result.num().negate();
1133    }
1134 }
1135 template <class Backend>
1136 void eval_divide(rational_adaptor<Backend>& result, const Backend& arg)
1137 {
1138    eval_divide(result, result, arg);
1139 }
1140 
1141 template <class Backend, class Arithmetic>
1142 typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && std::is_integral<Arithmetic>::value>::type
1143    eval_divide(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, Arithmetic arg)
1144 {
1145    if (&result == &a)
1146       return eval_divide(result, arg);
1147    if (arg == 0)
1148    {
1149       BOOST_MP_THROW_EXCEPTION(std::overflow_error("Integer division by zero"));
1150       return;
1151    }
1152    else if (arg == 1)
1153    {
1154       result = a;
1155       return;
1156    }
1157    else if (is_minus_one(arg))
1158    {
1159       result = a;
1160       result.num().negate();
1161       return;
1162    }
1163 
1164    if (eval_get_sign(a) == 0)
1165    {
1166       result = a;
1167       return;
1168    }
1169 
1170    using default_ops::eval_multiply;
1171    using default_ops::eval_divide;
1172    using default_ops::eval_gcd;
1173    using default_ops::eval_convert_to;
1174 
1175    Backend gcd;
1176    Arithmetic integer_gcd;
1177    eval_gcd(gcd, a.num(), arg);
1178    eval_convert_to(&integer_gcd, gcd);
1179    arg /= integer_gcd;
1180    eval_multiply(result.denom(), a.denom(), boost::multiprecision::detail::unsigned_abs(arg));
1181 
1182    if (integer_gcd > 1)
1183    {
1184       eval_divide(result.num(), a.num(), integer_gcd);
1185    }
1186    else
1187       result.num() = a.num();
1188    if (arg < 0)
1189    {
1190       result.num().negate();
1191    }
1192 }
1193 
1194 //
1195 // Increment and decrement:
1196 //
1197 template <class Backend> 
1198 inline void eval_increment(rational_adaptor<Backend>& arg)
1199 {
1200    using default_ops::eval_add;
1201    eval_add(arg.num(), arg.denom());
1202 }
1203 template <class Backend> 
1204 inline void eval_decrement(rational_adaptor<Backend>& arg)
1205 {
1206    using default_ops::eval_subtract;
1207    eval_subtract(arg.num(), arg.denom());
1208 }
1209 
1210 //
1211 // abs:
1212 //
1213 template <class Backend> 
1214 inline void eval_abs(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& arg)
1215 {
1216    using default_ops::eval_abs;
1217    eval_abs(result.num(), arg.num());
1218    result.denom() = arg.denom();
1219 }
1220 
1221 } // namespace backends
1222 
1223 //
1224 // Define a category for this number type, one of:
1225 // 
1226 //    number_kind_integer
1227 //    number_kind_floating_point
1228 //    number_kind_rational
1229 //    number_kind_fixed_point
1230 //    number_kind_complex
1231 //
1232 template<class Backend>
1233 struct number_category<rational_adaptor<Backend> > : public std::integral_constant<int, number_kind_rational>
1234 {};
1235 
1236 template <class Backend, expression_template_option ExpressionTemplates>
1237 struct component_type<number<rational_adaptor<Backend>, ExpressionTemplates> >
1238 {
1239    typedef number<Backend, ExpressionTemplates> type;
1240 };
1241 
1242 template <class IntBackend, expression_template_option ET>
1243 inline number<IntBackend, ET> numerator(const number<rational_adaptor<IntBackend>, ET>& val)
1244 {
1245    return val.backend().num();
1246 }
1247 template <class IntBackend, expression_template_option ET>
1248 inline number<IntBackend, ET> denominator(const number<rational_adaptor<IntBackend>, ET>& val)
1249 {
1250    return val.backend().denom();
1251 }
1252 
1253 template <class Backend>
1254 struct is_unsigned_number<rational_adaptor<Backend> > : public is_unsigned_number<Backend>
1255 {};
1256 
1257 
1258 }} // namespace boost::multiprecision
1259 
1260 namespace std {
1261 
1262    template <class IntBackend, boost::multiprecision::expression_template_option ExpressionTemplates>
1263    class numeric_limits<boost::multiprecision::number<boost::multiprecision::rational_adaptor<IntBackend>, ExpressionTemplates> > : public std::numeric_limits<boost::multiprecision::number<IntBackend, ExpressionTemplates> >
1264    {
1265       using base_type = std::numeric_limits<boost::multiprecision::number<IntBackend> >;
1266       using number_type = boost::multiprecision::number<boost::multiprecision::rational_adaptor<IntBackend> >;
1267 
1268    public:
1269       static constexpr bool is_integer = false;
1270       static constexpr bool is_exact = true;
1271       static constexpr      number_type(min)() { return (base_type::min)(); }
1272       static constexpr      number_type(max)() { return (base_type::max)(); }
1273       static constexpr number_type lowest() { return -(max)(); }
1274       static constexpr number_type epsilon() { return base_type::epsilon(); }
1275       static constexpr number_type round_error() { return epsilon() / 2; }
1276       static constexpr number_type infinity() { return base_type::infinity(); }
1277       static constexpr number_type quiet_NaN() { return base_type::quiet_NaN(); }
1278       static constexpr number_type signaling_NaN() { return base_type::signaling_NaN(); }
1279       static constexpr number_type denorm_min() { return base_type::denorm_min(); }
1280    };
1281 
1282    template <class IntBackend, boost::multiprecision::expression_template_option ExpressionTemplates>
1283    constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::rational_adaptor<IntBackend>, ExpressionTemplates> >::is_integer;
1284    template <class IntBackend, boost::multiprecision::expression_template_option ExpressionTemplates>
1285    constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::rational_adaptor<IntBackend>, ExpressionTemplates> >::is_exact;
1286 
1287 } // namespace std
1288 
1289 #endif