int256.cc 6.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204
  1. /*
  2. * Copyright 2018 Google LLC.
  3. * Licensed under the Apache License, Version 2.0 (the "License");
  4. * you may not use this file except in compliance with the License.
  5. * You may obtain a copy of the License at
  6. *
  7. * https://www.apache.org/licenses/LICENSE-2.0
  8. *
  9. * Unless required by applicable law or agreed to in writing, software
  10. * distributed under the License is distributed on an "AS IS" BASIS,
  11. * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  12. * See the License for the specific language governing permissions and
  13. * limitations under the License.
  14. */
  15. #include "int256.h"
  16. #include <iomanip>
  17. #include <iostream>
  18. #include <sstream>
  19. #include <glog/logging.h>
  20. #include "absl/numeric/int128.h"
  21. namespace rlwe {
  22. // Returns the 0-based position of the last set bit (i.e., most significant bit)
  23. // in the given uint64. The argument may not be 0.
  24. //
  25. // For example:
  26. // Given: 5 (decimal) == 101 (binary)
  27. // Returns: 2
  28. #define STEP(T, n, pos, sh) \
  29. do { \
  30. if ((n) >= (static_cast<T>(1) << (sh))) { \
  31. (n) = (n) >> (sh); \
  32. (pos) |= (sh); \
  33. } \
  34. } while (0)
  35. static inline int Fls64(Uint64 n) {
  36. //DCHECK_NE(0, n);
  37. int pos = 0;
  38. STEP(Uint64, n, pos, 0x20);
  39. Uint32 n32 = n;
  40. STEP(Uint32, n32, pos, 0x10);
  41. STEP(Uint32, n32, pos, 0x08);
  42. STEP(Uint32, n32, pos, 0x04);
  43. return pos + ((static_cast<Uint64>(0x3333333322221100) >> (n32 << 2)) & 0x3);
  44. }
  45. #undef STEP
  46. // Like Fls64() above, but returns the 0-based position of the last set bit
  47. // (i.e., most significant bit) in the given uint128. The argument may not be 0.
  48. static inline int Fls128(absl::uint128 n) {
  49. if (Uint64 hi = absl::Uint128High64(n)) {
  50. return Fls64(hi) + 64;
  51. }
  52. return Fls64(absl::Uint128Low64(n));
  53. }
  54. // Like Fls128() above, but returns the 0-based position of the last set bit
  55. // (i.e., most significant bit) in the given uint256. The argument may not be 0.
  56. static inline int Fls256(uint256 n) {
  57. absl::uint128 hi = Uint256High128(n);
  58. if (hi != 0) {
  59. return Fls128(hi) + 128;
  60. }
  61. return Fls128(Uint256Low128(n));
  62. }
  63. // Long division/modulo for uint256 implemented using the shift-subtract
  64. // division algorithm adapted from:
  65. // http://stackoverflow.com/questions/5386377/division-without-using
  66. void uint256::DivModImpl(uint256 dividend, uint256 divisor,
  67. uint256* quotient_ret, uint256* remainder_ret) {
  68. if (divisor == static_cast<uint256>(0)) {
  69. LOG(FATAL) << "Division or mod by zero: dividend.hi=" << dividend.hi_
  70. << ", lo=" << dividend.lo_;
  71. }
  72. if (divisor > dividend) {
  73. *quotient_ret = 0;
  74. *remainder_ret = dividend;
  75. return;
  76. }
  77. if (divisor == dividend) {
  78. *quotient_ret = 1;
  79. *remainder_ret = 0;
  80. return;
  81. }
  82. uint256 denominator = divisor;
  83. uint256 quotient = 0;
  84. // Left aligns the MSB of the denominator and the dividend.
  85. const int shift = Fls256(dividend) - Fls256(denominator);
  86. denominator <<= shift;
  87. // Uses shift-subtract algorithm to divide dividend by denominator. The
  88. // remainder will be left in dividend.
  89. for (int i = 0; i <= shift; ++i) {
  90. quotient <<= 1;
  91. if (dividend >= denominator) {
  92. dividend -= denominator;
  93. quotient |= 1;
  94. }
  95. denominator >>= 1;
  96. }
  97. *quotient_ret = quotient;
  98. *remainder_ret = dividend;
  99. }
  100. uint256& uint256::operator/=(const uint256& divisor) {
  101. uint256 quotient = 0;
  102. uint256 remainder = 0;
  103. DivModImpl(*this, divisor, &quotient, &remainder);
  104. *this = quotient;
  105. return *this;
  106. }
  107. uint256& uint256::operator%=(const uint256& divisor) {
  108. uint256 quotient = 0;
  109. uint256 remainder = 0;
  110. DivModImpl(*this, divisor, &quotient, &remainder);
  111. *this = remainder;
  112. return *this;
  113. }
  114. std::ostream& operator<<(std::ostream& o, const uint256& b) {
  115. std::ios_base::fmtflags flags = o.flags();
  116. // Select a divisor which is the largest power of the base < 2^64.
  117. uint256 div;
  118. std::streamsize div_base_log;
  119. switch (flags & std::ios::basefield) {
  120. case std::ios::hex:
  121. div = static_cast<Uint64>(0x1000000000000000); // 16^15
  122. div_base_log = 15;
  123. break;
  124. case std::ios::oct:
  125. div = static_cast<Uint64>(01000000000000000000000); // 8^21
  126. div_base_log = 21;
  127. break;
  128. default: // std::ios::dec
  129. div = static_cast<Uint64>(10000000000000000000ull); // 10^19
  130. div_base_log = 19;
  131. break;
  132. }
  133. // Now piece together the uint256 representation from five chunks of
  134. // the original value, each less than "div" and therefore representable
  135. // as a uint64.
  136. std::ostringstream os;
  137. std::ios_base::fmtflags copy_mask =
  138. std::ios::basefield | std::ios::showbase | std::ios::uppercase;
  139. os.setf(flags & copy_mask, copy_mask);
  140. uint256 high = b;
  141. uint256 low;
  142. uint256::DivModImpl(high, div, &high, &low);
  143. uint256 mid_low;
  144. uint256::DivModImpl(high, div, &high, &mid_low);
  145. uint256 mid_mid;
  146. uint256::DivModImpl(high, div, &high, &mid_mid);
  147. uint256 mid_high;
  148. uint256::DivModImpl(high, div, &high, &mid_high);
  149. bool print = high.lo_ != 0;
  150. if (print) {
  151. os << high.lo_;
  152. os << std::noshowbase << std::setfill('0') << std::setw(div_base_log);
  153. }
  154. print |= mid_high.lo_ != 0;
  155. if (print) {
  156. os << mid_high.lo_;
  157. os << std::noshowbase << std::setfill('0') << std::setw(div_base_log);
  158. }
  159. print |= mid_mid.lo_ != 0;
  160. if (print) {
  161. os << mid_mid.lo_;
  162. os << std::noshowbase << std::setfill('0') << std::setw(div_base_log);
  163. }
  164. print |= mid_low.lo_ != 0;
  165. if (print) {
  166. os << mid_low.lo_;
  167. os << std::noshowbase << std::setfill('0') << std::setw(div_base_log);
  168. }
  169. os << low.lo_;
  170. std::string rep = os.str();
  171. // Add the requisite padding.
  172. std::streamsize width = o.width(0);
  173. if (width > rep.size()) {
  174. if ((flags & std::ios::adjustfield) == std::ios::left) {
  175. rep.append(width - rep.size(), o.fill());
  176. } else {
  177. rep.insert(0, width - rep.size(), o.fill());
  178. }
  179. }
  180. // Stream the final representation in a single "<<" call.
  181. return o << rep;
  182. }
  183. } // namespace rlwe