boost/random/detail/large_arithmetic.hpp

   1 /* boost random/detail/large_arithmetic.hpp header file
   2  *
   3  * Copyright Steven Watanabe 2011
   4  * Distributed under the Boost Software License, Version 1.0. (See
   5  * accompanying file LICENSE_1_0.txt or copy at
   6  * http://www.boost.org/LICENSE_1_0.txt)
   7  *
   8  * See http://www.boost.org for most recent version including documentation.
   9  *
  10  * $Id: large_arithmetic.hpp 71018 2011-04-05 21:27:52Z steven_watanabe $
  11  */
  12
  13 #ifndef BOOST_RANDOM_DETAIL_LARGE_ARITHMETIC_HPP
  14 #define BOOST_RANDOM_DETAIL_LARGE_ARITHMETIC_HPP
  15
  16 #include <boost/cstdint.hpp>
  17 #include <boost/integer.hpp>
  18 #include <boost/limits.hpp>
  19 #include <boost/random/detail/integer_log2.hpp>
  20
  21 #include <boost/random/detail/disable_warnings.hpp>
  22
  23 namespace boost {
  24 namespace random {
  25 namespace detail {
  26
  27 struct div_t {
  28     boost::uintmax_t quotient;
  29     boost::uintmax_t remainder;
  30 };
  31
  32 inline div_t muldivmod(boost::uintmax_t a, boost::uintmax_t b, boost::uintmax_t m)
  33 {
  34     static const int bits =
  35         ::std::numeric_limits< ::boost::uintmax_t>::digits / 2;
  36     static const ::boost::uintmax_t mask = (::boost::uintmax_t(1) << bits) - 1;
  37     typedef ::boost::uint_t<bits>::fast digit_t;
  38
  39     int shift = std::numeric_limits< ::boost::uintmax_t>::digits - 1
  40         - detail::integer_log2(m);
  41
  42     a <<= shift;
  43     m <<= shift;
  44
  45     digit_t product[4] = { 0, 0, 0, 0 };
  46     digit_t a_[2] = { digit_t(a & mask), digit_t((a >> bits) & mask) };
  47     digit_t b_[2] = { digit_t(b & mask), digit_t((b >> bits) & mask) };
  48     digit_t m_[2] = { digit_t(m & mask), digit_t((m >> bits) & mask) };
  49
  50     // multiply a * b
  51     for(int i = 0; i < 2; ++i) {
  52         digit_t carry = 0;
  53         for(int j = 0; j < 2; ++j) {
  54             ::boost::uint64_t temp = ::boost::uintmax_t(a_[i]) * b_[j] +
  55                 carry + product[i + j];
  56             product[i + j] = digit_t(temp & mask);
  57             carry = digit_t(temp >> bits);
  58         }
  59         if(carry != 0) {
  60             product[i + 2] += carry;
  61         }
  62     }
  63
  64     digit_t quotient[2];
  65
  66     if(m == 0) {
  67         div_t result = {
  68             ((::boost::uintmax_t(product[3]) << bits) | product[2]),
  69             ((::boost::uintmax_t(product[1]) << bits) | product[0]) >> shift,
  70         };
  71         return result;
  72     }
  73
  74     // divide product / m
  75     for(int i = 3; i >= 2; --i) {
  76         ::boost::uintmax_t temp =
  77             ::boost::uintmax_t(product[i]) << bits | product[i - 1];
  78
  79         digit_t q = digit_t((product[i] == m_[1]) ? mask : temp / m_[1]);
  80
  81         ::boost::uintmax_t rem =
  82             ((temp - ::boost::uintmax_t(q) * m_[1]) << bits) + product[i - 2];
  83
  84         ::boost::uintmax_t diff = m_[0] * ::boost::uintmax_t(q);
  85
  86         int error = 0;
  87         if(diff > rem) {
  88             if(diff - rem > m) {
  89                 error = 2;
  90             } else {
  91                 error = 1;
  92             }
  93         }
  94         q -= error;
  95         rem = rem + error * m - diff;
  96
  97         quotient[i - 2] = q;
  98         product[i] = 0;
  99         product[i-1] = (rem >> bits) & mask;
 100         product[i-2] = rem & mask;
 101     }
 102
 103     div_t result = {
 104         ((::boost::uintmax_t(quotient[1]) << bits) | quotient[0]),
 105         ((::boost::uintmax_t(product[1]) << bits) | product[0]) >> shift,
 106     };
 107     return result;
 108 }
 109
 110 inline boost::uintmax_t muldiv(boost::uintmax_t a, boost::uintmax_t b, boost::uintmax_t m)
 111 { return detail::muldivmod(a, b, m).quotient; }
 112
 113 inline boost::uintmax_t mulmod(boost::uintmax_t a, boost::uintmax_t b, boost::uintmax_t m)
 114 { return detail::muldivmod(a, b, m).remainder; }
 115
 116 } // namespace detail
 117 } // namespace random
 118 } // namespace boost
 119
 120 #include <boost/random/detail/enable_warnings.hpp>
 121
 122 #endif // BOOST_RANDOM_DETAIL_LARGE_ARITHMETIC_HPP