boost/multiprecision/integer.hpp

0001 ///////////////////////////////////////////////////////////////
0002 //  Copyright 2012-21 John Maddock.
0003 //  Copyright 2021 Iskandarov Lev. Distributed under the Boost
0004 //  Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
0006
0007 #ifndef BOOST_MP_INTEGER_HPP
0008 #define BOOST_MP_INTEGER_HPP
0009
0010 #include <type_traits>
0011 #include <boost/multiprecision/cpp_int.hpp>
0012 #include <boost/multiprecision/detail/bitscan.hpp>
0013 #include <boost/multiprecision/detail/no_exceptions_support.hpp>
0014 #include <boost/multiprecision/detail/standalone_config.hpp>
0015
0016 namespace boost {
0017 namespace multiprecision {
0018
0019 template <class Integer, class I2>
0020 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value && boost::multiprecision::detail::is_integral<I2>::value, Integer&>::type
0021 multiply(Integer& result, const I2& a, const I2& b)
0022 {
0023    return result = static_cast<Integer>(a) * static_cast<Integer>(b);
0024 }
0025 template <class Integer, class I2>
0026 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value && boost::multiprecision::detail::is_integral<I2>::value, Integer&>::type
0027 add(Integer& result, const I2& a, const I2& b)
0028 {
0029    return result = static_cast<Integer>(a) + static_cast<Integer>(b);
0030 }
0031 template <class Integer, class I2>
0032 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value && boost::multiprecision::detail::is_integral<I2>::value, Integer&>::type
0033 subtract(Integer& result, const I2& a, const I2& b)
0034 {
0035    return result = static_cast<Integer>(a) - static_cast<Integer>(b);
0036 }
0037
0038 template <class Integer>
0039 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value>::type divide_qr(const Integer& x, const Integer& y, Integer& q, Integer& r)
0040 {
0041    q = x / y;
0042    r = x % y;
0043 }
0044
0045 template <class I1, class I2>
0046 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<I1>::value && boost::multiprecision::detail::is_integral<I2>::value, I2>::type integer_modulus(const I1& x, I2 val)
0047 {
0048    return static_cast<I2>(x % val);
0049 }
0050
0051 namespace detail {
0052 //
0053 // Figure out the kind of integer that has twice as many bits as some builtin
0054 // integer type I.  Use a native type if we can (including types which may not
0055 // be recognised by boost::int_t because they're larger than long long),
0056 // otherwise synthesize a cpp_int to do the job.
0057 //
0058 template <class I>
0059 struct double_integer
0060 {
0061    static constexpr const unsigned int_t_digits =
0062        2 * sizeof(I) <= sizeof(long long) ? std::numeric_limits<I>::digits * 2 : 1;
0063
0064    using type = typename std::conditional<
0065        2 * sizeof(I) <= sizeof(long long),
0066        typename std::conditional<
0067            boost::multiprecision::detail::is_signed<I>::value && boost::multiprecision::detail::is_integral<I>::value,
0068            typename boost::multiprecision::detail::int_t<int_t_digits>::least,
0069            typename boost::multiprecision::detail::uint_t<int_t_digits>::least>::type,
0070        typename std::conditional<
0071            2 * sizeof(I) <= sizeof(double_limb_type),
0072            typename std::conditional<
0073                boost::multiprecision::detail::is_signed<I>::value && boost::multiprecision::detail::is_integral<I>::value,
0074                signed_double_limb_type,
0075                double_limb_type>::type,
0076            number<cpp_int_backend<sizeof(I) * CHAR_BIT * 2, sizeof(I) * CHAR_BIT * 2, (boost::multiprecision::detail::is_signed<I>::value ? signed_magnitude : unsigned_magnitude), unchecked, void> > >::type>::type;
0077 };
0078
0079 } // namespace detail
0080
0081 template <class I1, class I2, class I3>
0082 BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<I1>::value && boost::multiprecision::detail::is_unsigned<I2>::value && boost::multiprecision::detail::is_integral<I3>::value, I1>::type
0083 powm(const I1& a, I2 b, I3 c)
0084 {
0085    using double_type = typename detail::double_integer<I1>::type;
0086
0087    I1          x(1), y(a);
0088    double_type result(0);
0089
0090    while (b > 0)
0091    {
0092       if (b & 1)
0093       {
0094          multiply(result, x, y);
0095          x = integer_modulus(result, c);
0096       }
0097       multiply(result, y, y);
0098       y = integer_modulus(result, c);
0099       b >>= 1;
0100    }
0101    return x % c;
0102 }
0103
0104 template <class I1, class I2, class I3>
0105 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<I1>::value && boost::multiprecision::detail::is_signed<I2>::value && boost::multiprecision::detail::is_integral<I2>::value && boost::multiprecision::detail::is_integral<I3>::value, I1>::type
0106 powm(const I1& a, I2 b, I3 c)
0107 {
0108    if (b < 0)
0109    {
0110       BOOST_MP_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
0111    }
0112    return powm(a, static_cast<typename boost::multiprecision::detail::make_unsigned<I2>::type>(b), c);
0113 }
0114
0115 template <class Integer>
0116 BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, std::size_t>::type lsb(const Integer& val)
0117 {
0118    if (val <= 0)
0119    {
0120       if (val == 0)
0121       {
0122          BOOST_MP_THROW_EXCEPTION(std::domain_error("No bits were set in the operand."));
0123       }
0124       else
0125       {
0126          BOOST_MP_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined."));
0127       }
0128    }
0129    return detail::find_lsb(val);
0130 }
0131
0132 template <class Integer>
0133 BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, std::size_t>::type msb(Integer val)
0134 {
0135    if (val <= 0)
0136    {
0137       if (val == 0)
0138       {
0139          BOOST_MP_THROW_EXCEPTION(std::domain_error("No bits were set in the operand."));
0140       }
0141       else
0142       {
0143          BOOST_MP_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined."));
0144       }
0145    }
0146    return detail::find_msb(val);
0147 }
0148
0149 template <class Integer>
0150 BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, bool>::type bit_test(const Integer& val, std::size_t index)
0151 {
0152    Integer mask = 1;
0153    if (index >= sizeof(Integer) * CHAR_BIT)
0154       return 0;
0155    if (index)
0156       mask <<= index;
0157    return val & mask ? true : false;
0158 }
0159
0160 template <class Integer>
0161 BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, Integer&>::type bit_set(Integer& val, std::size_t index)
0162 {
0163    Integer mask = 1;
0164    if (index >= sizeof(Integer) * CHAR_BIT)
0165       return val;
0166    if (index)
0167       mask <<= index;
0168    val |= mask;
0169    return val;
0170 }
0171
0172 template <class Integer>
0173 BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, Integer&>::type bit_unset(Integer& val, std::size_t index)
0174 {
0175    Integer mask = 1;
0176    if (index >= sizeof(Integer) * CHAR_BIT)
0177       return val;
0178    if (index)
0179       mask <<= index;
0180    val &= ~mask;
0181    return val;
0182 }
0183
0184 template <class Integer>
0185 BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, Integer&>::type bit_flip(Integer& val, std::size_t index)
0186 {
0187    Integer mask = 1;
0188    if (index >= sizeof(Integer) * CHAR_BIT)
0189       return val;
0190    if (index)
0191       mask <<= index;
0192    val ^= mask;
0193    return val;
0194 }
0195
0196 namespace detail {
0197
0198 template <class Integer>
0199 BOOST_MP_CXX14_CONSTEXPR Integer karatsuba_sqrt(const Integer& x, Integer& r, size_t bits)
0200 {
0201    //
0202    // Define the floating point type used for std::sqrt, in our tests, sqrt(double) and sqrt(long double) take
0203    // about the same amount of time as long as long double is not an emulated 128-bit type (ie the same type
0204    // as __float128 from libquadmath).  So only use long double if it's an 80-bit type:
0205    //
0206 #ifndef __clang__
0207    typedef typename std::conditional<(std::numeric_limits<long double>::digits == 64), long double, double>::type real_cast_type;
0208 #else
0209    // clang has buggy __int128 -> long double conversion:
0210    typedef double real_cast_type;
0211 #endif
0212    //
0213    // As per the Karatsuba sqrt algorithm, the low order bits/4 bits pay no part in the result, only in the remainder,
0214    // so define the number of bits our argument must have before passing to std::sqrt is safe, even if doing so
0215    // looses a few bits:
0216    //
0217    constexpr std::size_t cutoff = (std::numeric_limits<real_cast_type>::digits * 4) / 3;
0218    //
0219    // Type which can hold at least "cutoff" bits:
0220    //
0221 #ifdef BOOST_HAS_INT128
0222    using cutoff_t = typename std::conditional<(cutoff > 64), uint128_type, std::uint64_t>::type;
0223 #else
0224    using cutoff_t = std::uint64_t;
0225 #endif
0226    //
0227    // See if we can take the fast path:
0228    //
0229    if (bits <= cutoff)
0230    {
0231       constexpr cutoff_t half_bits = (cutoff_t(1u) << ((sizeof(cutoff_t) * CHAR_BIT) / 2)) - 1;
0232       cutoff_t       val = static_cast<cutoff_t>(x);
0233       real_cast_type real_val = static_cast<real_cast_type>(val);
0234       cutoff_t       s64 = static_cast<cutoff_t>(std::sqrt(real_val));
0235       // converting to long double can loose some precision, and `sqrt` can give eps error, so we'll fix this
0236       // this is needed
0237       while ((s64 > half_bits) || (s64 * s64 > val))
0238          s64--;
0239       // in my tests this never fired, but theoretically this might be needed
0240       while ((s64 < half_bits) && ((s64 + 1) * (s64 + 1) <= val))
0241          s64++;
0242       r = static_cast<Integer>(val - s64 * s64);
0243       return static_cast<Integer>(s64);
0244    }
0245    // https://hal.inria.fr/file/index/docid/72854/filename/RR-3805.pdf
0246    std::size_t b = bits / 4;
0247    Integer q = x;
0248    q >>= b * 2;
0249    Integer s = karatsuba_sqrt(q, r, bits - b * 2);
0250    Integer t = 0u;
0251    bit_set(t, static_cast<unsigned>(b * 2));
0252    r <<= b;
0253    t--;
0254    t &= x;
0255    t >>= b;
0256    t += r;
0257    s <<= 1;
0258    divide_qr(t, s, q, r);
0259    r <<= b;
0260    t = 0u;
0261    bit_set(t, static_cast<unsigned>(b));
0262    t--;
0263    t &= x;
0264    r += t;
0265    s <<= (b - 1); // we already <<1 it before
0266    s += q;
0267    q *= q;
0268    // we substract after, so it works for unsigned integers too
0269    if (r < q)
0270    {
0271       t = s;
0272       t <<= 1;
0273       t--;
0274       r += t;
0275       s--;
0276    }
0277    r -= q;
0278    return s;
0279 }
0280
0281 template <class Integer>
0282 BOOST_MP_CXX14_CONSTEXPR Integer bitwise_sqrt(const Integer& x, Integer& r)
0283 {
0284    //
0285    // This is slow bit-by-bit integer square root, see for example
0286    // http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_.28base_2.29
0287    // There are better methods such as http://hal.inria.fr/docs/00/07/28/54/PDF/RR-3805.pdf
0288    // and http://hal.inria.fr/docs/00/07/21/13/PDF/RR-4475.pdf which should be implemented
0289    // at some point.
0290    //
0291    Integer s = 0;
0292    switch (x)
0293    {
0294    case 0:
0295       r = 0;
0296       return s;
0297    case 1:
0298       r = 0;
0299       return 1;
0300    case 2:
0301       r = 1;
0302       return 1;
0303    case 3:
0304       r = 2;
0305       return 1;
0306    default:
0307       break;
0308       // fall through:
0309    }
0310    std::ptrdiff_t g = msb(x);
0311
0312    Integer t = 0;
0313    r = x;
0314    g /= 2;
0315    bit_set(s, g);
0316    bit_set(t, 2 * g);
0317    r = x - t;
0318    --g;
0319    do
0320    {
0321       t = s;
0322       t <<= g + 1;
0323       bit_set(t, 2 * g);
0324       if (t <= r)
0325       {
0326          bit_set(s, g);
0327          r -= t;
0328       }
0329       --g;
0330    } while (g >= 0);
0331    return s;
0332 }
0333
0334 } // namespace detail
0335
0336 template <class Integer>
0337 BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, Integer>::type sqrt(const Integer& x, Integer& r)
0338 {
0339 #ifndef BOOST_MP_NO_CONSTEXPR_DETECTION
0340    // recursive Karatsuba sqrt can cause issues in constexpr context:
0341    if (BOOST_MP_IS_CONST_EVALUATED(x))
0342    {
0343       return detail::bitwise_sqrt(x, r);
0344    }
0345 #endif
0346    if (x == 0u) {
0347       r = 0u;
0348       return 0u;
0349    }
0350
0351    return detail::karatsuba_sqrt(x, r, msb(x) + 1);
0352 }
0353
0354 template <class Integer>
0355 BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, Integer>::type sqrt(const Integer& x)
0356 {
0357    Integer r(0);
0358    return sqrt(x, r);
0359 }
0360
0361 }} // namespace boost::multiprecision
0362
0363 #endif