Teuchos - Trilinos Tools Package  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Teuchos_BigUInt.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Teuchos: Common Tools Package
4 //
5 // Copyright 2004 NTESS and the Teuchos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef TEUCHOS_BIG_UINT_HPP
11 #define TEUCHOS_BIG_UINT_HPP
12 
13 #include <iostream>
14 #include <Teuchos_BigUIntDecl.hpp>
15 
20 namespace Teuchos {
21 
22 template <int n>
23 BigUInt<n>::BigUInt() {}
24 
25 template <int n>
26 BigUInt<n>::BigUInt(std::uint64_t v) {
27  for (int i = 2; i < n; ++i) x[i] = 0;
28  if (n > 1) x[1] = std::uint32_t(v >> 32);
29  x[0] = std::uint32_t(v);
30 }
31 
32 template <int n>
33 BigUInt<n>::BigUInt(std::string const& s) : BigUInt(std::uint32_t(0)) {
34  for (auto c : s) {
35  operator*=(10);
36  operator+=(c - '0');
37  }
38 }
39 
40 template <int n>
41 BigUInt<n>::operator bool() const noexcept {
42  for (int i = 0; i < n; ++i) if (x[i]) return true;
43  return false;
44 }
45 
46 template <int n>
47 std::uint32_t& BigUInt<n>::operator[](int i) { return x[i]; }
48 
49 template <int n>
50 std::uint32_t const& BigUInt<n>::operator[](int i) const { return x[i]; }
51 
52 template <int n>
53 BigUInt<n>& BigUInt<n>::operator+=(std::uint32_t b) {
54  std::uint32_t carry = b;
55  for (int i = 0; i < n; ++i) {
56  std::uint64_t ax = x[i];
57  auto cx = ax + std::uint64_t(carry);
58  carry = std::uint32_t(cx >> 32);
59  x[i] = std::uint32_t(cx);
60  }
61  return *this;
62 }
63 
64 template <int n>
65 BigUInt<n>& BigUInt<n>::operator+=(BigUInt<n> const& b) {
66  std::uint32_t carry = 0;
67  for (int i = 0; i < n; ++i) {
68  std::uint64_t ax = x[i];
69  std::uint64_t bx = b[i];
70  auto cx = ax + bx + std::uint64_t(carry);
71  carry = std::uint32_t(cx >> 32);
72  x[i] = std::uint32_t(cx);
73  }
74  return *this;
75 }
76 
77 template <int n>
78 BigUInt<n>& BigUInt<n>::operator-=(std::uint32_t b) {
79  std::int64_t carry = b;
80  for (int i = 0; i < n; ++i) {
81  std::int64_t ax = x[i];
82  auto cx = ax - carry;
83  if (cx < 0) {
84  carry = 1;
85  cx += std::int64_t(1) << 32;
86  } else {
87  carry = 0;
88  }
89  x[i] = std::uint32_t(cx);
90  }
91  return *this;
92 }
93 
94 template <int n>
95 BigUInt<n>& BigUInt<n>::operator-=(BigUInt<n> const& b) {
96  std::int64_t carry = 0;
97  for (int i = 0; i < n; ++i) {
98  std::int64_t ax = x[i];
99  std::int64_t bx = b[i];
100  auto cx = ax - bx - carry;
101  if (cx < 0) {
102  carry = 1;
103  cx += std::int64_t(1) << 32;
104  } else {
105  carry = 0;
106  }
107  x[i] = std::uint32_t(cx);
108  }
109  return *this;
110 }
111 
112 template <int n>
113 BigUInt<n>& BigUInt<n>::operator*=(std::uint32_t b) {
114  std::uint32_t carry = 0;
115  for (int i = 0; i < n; ++i) {
116  std::uint64_t ax = x[i];
117  auto cx = (ax * std::uint64_t(b)) + std::uint64_t(carry);
118  carry = std::uint32_t(cx >> 32);
119  x[i] = std::uint32_t(cx);
120  }
121  return *this;
122 }
123 
124 template <int n>
125 BigUInt<n>& BigUInt<n>::operator<<=(std::uint32_t b) {
126  auto ndigits = b / 32;
127  auto nbits = b - (ndigits * 32);
128  for (int i = n - 1; i >= 0; --i) {
129  std::uint32_t xi = 0;
130  if (i >= int(ndigits)) {
131  xi = x[i - ndigits] << nbits;
132  }
133  // nbits &&, because apparently shifting a 32-bit value by 32 is not allowed
134  if (nbits && (i > int(ndigits))) {
135  xi |= x[i - ndigits - 1] >> (32 - nbits);
136  }
137  x[i] = xi;
138  }
139  return *this;
140 }
141 
142 template <int n>
143 BigUInt<n>& BigUInt<n>::operator>>=(std::uint32_t b) {
144  auto ndigits = b / 32;
145  auto nbits = b - (ndigits * 32);
146  for (int i = 0; i < n; ++i) {
147  std::uint32_t xi = 0;
148  if (i + ndigits < n) xi = x[i + ndigits] >> nbits;
149  if (nbits && i + ndigits + 1 < n) xi |= x[i + ndigits + 1] << (32 - nbits);
150  x[i] = xi;
151  }
152  return *this;
153 }
154 
155 template <int n>
156 std::ostream& operator<<(std::ostream& os, BigUInt<n> a) {
157  char buf[n * 20];
158  int i = 0;
159  while (a) {
160  BigUInt<n> quotient;
161  divmod(quotient, a, 10);
162  auto remainder = a[0];
163  a = quotient;
164  buf[i++] = char(remainder) + '0';
165  }
166  for (int j = 0; j < i / 2; ++j) {
167  auto tmp = buf[i - j - 1];
168  buf[i - j - 1] = buf[j];
169  buf[j] = tmp;
170  }
171  if (i == 0) buf[i++] = '0';
172  buf[i] = '\0';
173  os << buf;
174  return os;
175 }
176 
177 template <int n>
178 BigUInt<n> operator+(BigUInt<n> const& a, BigUInt<n> const& b) {
179  auto c = a;
180  c += b;
181  return c;
182 }
183 
184 template <int n>
185 BigUInt<n> operator-(BigUInt<n> const& a, BigUInt<n> const& b) {
186  auto c = a;
187  c -= b;
188  return c;
189 }
190 
191 template <int n>
192 BigUInt<n> operator*(BigUInt<n> const& a, BigUInt<n> const& b) {
193  BigUInt<n> a_times_b_i;
194  BigUInt<n> c{0};
195  for (int i = n - 1; i >= 0; --i) {
196  a_times_b_i = a;
197  a_times_b_i *= b[i];
198  c <<= 32;
199  c += a_times_b_i;
200  }
201  return c;
202 }
203 
204 template <int n>
205 BigUInt<n> operator/(BigUInt<n> const& a, std::uint32_t const& b) {
206  BigUInt<n> quotient;
207  auto x = a;
208  divmod(quotient, x, b);
209  return quotient;
210 }
211 
212 template <int n>
213 BigUInt<n> operator/(BigUInt<n> const& a, BigUInt<n> const& b) {
214  if (b > a) return BigUInt<n>(0);
215  BigUInt<n> quotient(1);
216  auto c = b;
217  while (c < a) {
218  c <<= 1;
219  quotient <<= 1;
220  }
221  auto factor = quotient;
222  factor >>= 1;
223  while (factor) {
224  int d = comp(a, c);
225  if (d == 0) break;
226  if (d == -1) {
227  c -= b * factor;
228  quotient -= factor;
229  } else {
230  c += b * factor;
231  quotient += factor;
232  }
233  factor >>= 1;
234  }
235  if (c > a) {
236  c -= b;
237  quotient -= 1;
238  }
239  return quotient;
240 }
241 
242 template <int n>
243 int comp(BigUInt<n> const& a, BigUInt<n> const& b) {
244  for (int i = n - 1; i >= 0; --i) {
245  if (a[i] != b[i]) {
246  if (a[i] > b[i]) return 1;
247  else return -1;
248  }
249  }
250  return 0;
251 }
252 
253 template <int n>
254 bool operator>=(BigUInt<n> const& a, BigUInt<n> const& b) {
255  return comp(a, b) > -1;
256 }
257 
258 template <int n>
259 bool operator<=(BigUInt<n> const& a, BigUInt<n> const& b) {
260  return comp(a, b) < 1;
261 }
262 
263 template <int n>
264 bool operator<(BigUInt<n> const& a, BigUInt<n> const& b) {
265  return comp(a, b) == -1;
266 }
267 
268 template <int n>
269 bool operator>(BigUInt<n> const& a, BigUInt<n> const& b) {
270  return comp(a, b) == 1;
271 }
272 
273 template <int n>
274 bool operator==(BigUInt<n> const& a, BigUInt<n> const& b) {
275  for (int i = 0; i < n; ++i) if (a[i] != b[i]) return false;
276  return true;
277 }
278 
279 template <int n>
280 void divmod(BigUInt<n>& quotient, BigUInt<n>& x, std::uint32_t const& b) {
281  quotient = BigUInt<n>(std::uint32_t(0));
282  for (int i = n - 1; i >= 0;) {
283  if (x[i]) {
284  if (x[i] >= b) {
285  auto dividend = x[i];
286  auto quotient2 = dividend / b;
287  auto remainder = dividend - (quotient2 * b);
288  quotient[i] = quotient2;
289  x[i] = remainder;
290  } else if (i > 0) {
291  auto dividend = std::uint64_t(x[i]);
292  dividend <<= 32;
293  dividend |= x[i - 1];
294  auto quotient2 = dividend / std::uint64_t(b);
295  auto remainder = dividend - (quotient2 * std::uint64_t(b));
296  quotient[i - 1] = std::uint32_t(quotient2);
297  x[i] = 0;
298  x[i - 1] = std::uint32_t(remainder);
299  i = i - 1;
300  } else {
301  break;
302  }
303  } else {
304  i = i - 1;
305  }
306  }
307 }
308 
309 }
310 
311 #endif
312 
Arbitrary-precision unsigned integer declaration.