Chris@16: ///////////////////////////////////////////////////////////////
Chris@16: //  Copyright 2012 John Maddock. Distributed under the Boost
Chris@16: //  Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_
Chris@16: 
Chris@16: #ifndef BOOST_MP_INTEGER_HPP
Chris@16: #define BOOST_MP_INTEGER_HPP
Chris@16: 
Chris@16: #include <boost/multiprecision/cpp_int.hpp>
Chris@16: #include <boost/multiprecision/detail/bitscan.hpp>
Chris@16: 
Chris@16: namespace boost{
Chris@16: namespace multiprecision{
Chris@16: 
Chris@16: template <class Integer, class I2>
Chris@16: typename enable_if_c<is_integral<Integer>::value && is_integral<I2>::value, Integer&>::type
Chris@16:    multiply(Integer& result, const I2& a, const I2& b)
Chris@16: {
Chris@16:    return result = static_cast<Integer>(a) * static_cast<Integer>(b);
Chris@16: }
Chris@16: template <class Integer, class I2>
Chris@16: typename enable_if_c<is_integral<Integer>::value && is_integral<I2>::value, Integer&>::type
Chris@16:    add(Integer& result, const I2& a, const I2& b)
Chris@16: {
Chris@16:    return result = static_cast<Integer>(a) + static_cast<Integer>(b);
Chris@16: }
Chris@16: template <class Integer, class I2>
Chris@16: typename enable_if_c<is_integral<Integer>::value && is_integral<I2>::value, Integer&>::type
Chris@16:    subtract(Integer& result, const I2& a, const I2& b)
Chris@16: {
Chris@16:    return result = static_cast<Integer>(a) - static_cast<Integer>(b);
Chris@16: }
Chris@16: 
Chris@16: template <class Integer>
Chris@16: typename enable_if_c<is_integral<Integer>::value>::type divide_qr(const Integer& x, const Integer& y, Integer& q, Integer& r)
Chris@16: {
Chris@16:    q = x / y;
Chris@16:    r = x % y;
Chris@16: }
Chris@16: 
Chris@16: template <class I1, class I2>
Chris@16: typename enable_if_c<is_integral<I1>::value && is_integral<I2>::value, I2>::type integer_modulus(const I1& x, I2 val)
Chris@16: {
Chris@16:    return static_cast<I2>(x % val);
Chris@16: }
Chris@16: 
Chris@16: namespace detail{
Chris@16: //
Chris@16: // Figure out the kind of integer that has twice as many bits as some builtin
Chris@16: // integer type I.  Use a native type if we can (including types which may not
Chris@16: // be recognised by boost::int_t because they're larger than long long),
Chris@16: // otherwise synthesize a cpp_int to do the job.
Chris@16: //
Chris@16: template <class I>
Chris@16: struct double_integer
Chris@16: {
Chris@16:    static const unsigned int_t_digits =
Chris@16:       2 * sizeof(I) <= sizeof(long long) ? std::numeric_limits<I>::digits * 2 : 1;
Chris@16: 
Chris@16:    typedef typename mpl::if_c<
Chris@16:       2 * sizeof(I) <= sizeof(long long),
Chris@16:       typename mpl::if_c<
Chris@16:          is_signed<I>::value,
Chris@16:          typename boost::int_t<int_t_digits>::least,
Chris@16:          typename boost::uint_t<int_t_digits>::least
Chris@16:       >::type,
Chris@16:       typename mpl::if_c<
Chris@16:          2 * sizeof(I) <= sizeof(double_limb_type),
Chris@16:          typename mpl::if_c<
Chris@16:             is_signed<I>::value,
Chris@16:             signed_double_limb_type,
Chris@16:             double_limb_type
Chris@16:          >::type,
Chris@16:          number<cpp_int_backend<sizeof(I)*CHAR_BIT*2, sizeof(I)*CHAR_BIT*2, (is_signed<I>::value ? signed_magnitude : unsigned_magnitude), unchecked, void> >
Chris@16:       >::type
Chris@16:    >::type type;
Chris@16: };
Chris@16: 
Chris@16: }
Chris@16: 
Chris@16: template <class I1, class I2, class I3>
Chris@101: typename enable_if_c<is_integral<I1>::value && is_unsigned<I2>::value && is_integral<I3>::value, I1>::type
Chris@16:    powm(const I1& a, I2 b, I3 c)
Chris@16: {
Chris@16:    typedef typename detail::double_integer<I1>::type double_type;
Chris@16: 
Chris@16:    I1 x(1), y(a);
Chris@16:    double_type result;
Chris@16: 
Chris@16:    while(b > 0)
Chris@16:    {
Chris@16:       if(b & 1)
Chris@16:       {
Chris@16:          multiply(result, x, y);
Chris@16:          x = integer_modulus(result, c);
Chris@16:       }
Chris@16:       multiply(result, y, y);
Chris@16:       y = integer_modulus(result, c);
Chris@16:       b >>= 1;
Chris@16:    }
Chris@16:    return x % c;
Chris@16: }
Chris@16: 
Chris@101: template <class I1, class I2, class I3>
Chris@101: inline typename enable_if_c<is_integral<I1>::value && is_signed<I2>::value && is_integral<I3>::value, I1>::type
Chris@101:    powm(const I1& a, I2 b, I3 c)
Chris@101: {
Chris@101:    if(b < 0)
Chris@101:    {
Chris@101:       BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
Chris@101:    }
Chris@101:    return powm(a, static_cast<typename make_unsigned<I2>::type>(b), c);
Chris@101: }
Chris@101: 
Chris@16: template <class Integer>
Chris@16: typename enable_if_c<is_integral<Integer>::value, unsigned>::type lsb(const Integer& val)
Chris@16: {
Chris@16:    if(val <= 0)
Chris@16:    {
Chris@16:       if(val == 0)
Chris@16:       {
Chris@16:          BOOST_THROW_EXCEPTION(std::range_error("No bits were set in the operand."));
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          BOOST_THROW_EXCEPTION(std::range_error("Testing individual bits in negative values is not supported - results are undefined."));
Chris@16:       }
Chris@16:    }
Chris@16:    return detail::find_lsb(val);
Chris@16: }
Chris@16: 
Chris@16: template <class Integer>
Chris@16: typename enable_if_c<is_integral<Integer>::value, unsigned>::type msb(Integer val)
Chris@16: {
Chris@16:    if(val <= 0)
Chris@16:    {
Chris@16:       if(val == 0)
Chris@16:       {
Chris@16:          BOOST_THROW_EXCEPTION(std::range_error("No bits were set in the operand."));
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          BOOST_THROW_EXCEPTION(std::range_error("Testing individual bits in negative values is not supported - results are undefined."));
Chris@16:       }
Chris@16:    }
Chris@16:    return detail::find_msb(val);
Chris@16: }
Chris@16: 
Chris@16: template <class Integer>
Chris@16: typename enable_if_c<is_integral<Integer>::value, bool>::type bit_test(const Integer& val, unsigned index)
Chris@16: {
Chris@16:    Integer mask = 1;
Chris@16:    if(index >= sizeof(Integer) * CHAR_BIT)
Chris@16:       return 0;
Chris@16:    if(index)
Chris@16:       mask <<= index;
Chris@16:    return val & mask ? true : false;
Chris@16: }
Chris@16: 
Chris@16: template <class Integer>
Chris@16: typename enable_if_c<is_integral<Integer>::value, Integer&>::type bit_set(Integer& val, unsigned index)
Chris@16: {
Chris@16:    Integer mask = 1;
Chris@16:    if(index >= sizeof(Integer) * CHAR_BIT)
Chris@16:       return val;
Chris@16:    if(index)
Chris@16:       mask <<= index;
Chris@16:    val |= mask;
Chris@16:    return val;
Chris@16: }
Chris@16: 
Chris@16: template <class Integer>
Chris@16: typename enable_if_c<is_integral<Integer>::value, Integer&>::type bit_unset(Integer& val, unsigned index)
Chris@16: {
Chris@16:    Integer mask = 1;
Chris@16:    if(index >= sizeof(Integer) * CHAR_BIT)
Chris@16:       return val;
Chris@16:    if(index)
Chris@16:       mask <<= index;
Chris@16:    val &= ~mask;
Chris@16:    return val;
Chris@16: }
Chris@16: 
Chris@16: template <class Integer>
Chris@16: typename enable_if_c<is_integral<Integer>::value, Integer&>::type bit_flip(Integer& val, unsigned index)
Chris@16: {
Chris@16:    Integer mask = 1;
Chris@16:    if(index >= sizeof(Integer) * CHAR_BIT)
Chris@16:       return val;
Chris@16:    if(index)
Chris@16:       mask <<= index;
Chris@16:    val ^= mask;
Chris@16:    return val;
Chris@16: }
Chris@16: 
Chris@16: template <class Integer>
Chris@16: typename enable_if_c<is_integral<Integer>::value, Integer>::type sqrt(const Integer& x, Integer& r)
Chris@16: {
Chris@16:    //
Chris@16:    // This is slow bit-by-bit integer square root, see for example
Chris@16:    // http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_.28base_2.29
Chris@16:    // There are better methods such as http://hal.inria.fr/docs/00/07/28/54/PDF/RR-3805.pdf
Chris@16:    // and http://hal.inria.fr/docs/00/07/21/13/PDF/RR-4475.pdf which should be implemented
Chris@16:    // at some point.
Chris@16:    //
Chris@16:    Integer s = 0;
Chris@16:    if(x == 0)
Chris@16:    {
Chris@16:       r = 0;
Chris@16:       return s;
Chris@16:    }
Chris@16:    int g = msb(x);
Chris@16:    if(g == 0)
Chris@16:    {
Chris@16:       r = 1;
Chris@16:       return s;
Chris@16:    }
Chris@16:    
Chris@101:    Integer t = 0;
Chris@16:    r = x;
Chris@16:    g /= 2;
Chris@16:    bit_set(s, g);
Chris@16:    bit_set(t, 2 * g);
Chris@16:    r = x - t;
Chris@16:    --g;
Chris@16:    do
Chris@16:    {
Chris@16:       t = s;
Chris@16:       t <<= g + 1;
Chris@16:       bit_set(t, 2 * g);
Chris@16:       if(t <= r)
Chris@16:       {
Chris@16:          bit_set(s, g);
Chris@16:          r -= t;
Chris@16:       }
Chris@16:       --g;
Chris@16:    }
Chris@16:    while(g >= 0);
Chris@16:    return s;
Chris@16: }
Chris@16: 
Chris@16: template <class Integer>
Chris@16: typename enable_if_c<is_integral<Integer>::value, Integer>::type sqrt(const Integer& x)
Chris@16: {
Chris@16:    Integer r;
Chris@16:    return sqrt(x, r);
Chris@16: }
Chris@16: 
Chris@16: }} // namespaces
Chris@16: 
Chris@16: #endif