Teuchos Package Browser (Single Doxygen Collection)  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 //
4 // Teuchos: Common Tools Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #ifndef TEUCHOS_BIG_UINT_HPP
43 #define TEUCHOS_BIG_UINT_HPP
44 
45 #include <iostream>
46 #include <Teuchos_BigUIntDecl.hpp>
47 
52 namespace Teuchos {
53 
54 template <int n>
56 
57 template <int n>
58 BigUInt<n>::BigUInt(std::uint64_t v) {
59  for (int i = 2; i < n; ++i) x[i] = 0;
60  if (n > 1) x[1] = std::uint32_t(v >> 32);
61  x[0] = std::uint32_t(v);
62 }
63 
64 template <int n>
65 BigUInt<n>::BigUInt(std::string const& s) : BigUInt(std::uint32_t(0)) {
66  for (auto c : s) {
67  operator*=(10);
68  operator+=(c - '0');
69  }
70 }
71 
72 template <int n>
73 BigUInt<n>::operator bool() const noexcept {
74  for (int i = 0; i < n; ++i) if (x[i]) return true;
75  return false;
76 }
77 
78 template <int n>
79 std::uint32_t& BigUInt<n>::operator[](int i) { return x[i]; }
80 
81 template <int n>
82 std::uint32_t const& BigUInt<n>::operator[](int i) const { return x[i]; }
83 
84 template <int n>
86  std::uint32_t carry = b;
87  for (int i = 0; i < n; ++i) {
88  std::uint64_t ax = x[i];
89  auto cx = ax + std::uint64_t(carry);
90  carry = std::uint32_t(cx >> 32);
91  x[i] = std::uint32_t(cx);
92  }
93  return *this;
94 }
95 
96 template <int n>
98  std::uint32_t carry = 0;
99  for (int i = 0; i < n; ++i) {
100  std::uint64_t ax = x[i];
101  std::uint64_t bx = b[i];
102  auto cx = ax + bx + std::uint64_t(carry);
103  carry = std::uint32_t(cx >> 32);
104  x[i] = std::uint32_t(cx);
105  }
106  return *this;
107 }
108 
109 template <int n>
111  std::int64_t carry = b;
112  for (int i = 0; i < n; ++i) {
113  std::int64_t ax = x[i];
114  auto cx = ax - carry;
115  if (cx < 0) {
116  carry = 1;
117  cx += std::int64_t(1) << 32;
118  } else {
119  carry = 0;
120  }
121  x[i] = std::uint32_t(cx);
122  }
123  return *this;
124 }
125 
126 template <int n>
128  std::int64_t carry = 0;
129  for (int i = 0; i < n; ++i) {
130  std::int64_t ax = x[i];
131  std::int64_t bx = b[i];
132  auto cx = ax - bx - carry;
133  if (cx < 0) {
134  carry = 1;
135  cx += std::int64_t(1) << 32;
136  } else {
137  carry = 0;
138  }
139  x[i] = std::uint32_t(cx);
140  }
141  return *this;
142 }
143 
144 template <int n>
146  std::uint32_t carry = 0;
147  for (int i = 0; i < n; ++i) {
148  std::uint64_t ax = x[i];
149  auto cx = (ax * std::uint64_t(b)) + std::uint64_t(carry);
150  carry = std::uint32_t(cx >> 32);
151  x[i] = std::uint32_t(cx);
152  }
153  return *this;
154 }
155 
156 template <int n>
158  auto ndigits = b / 32;
159  auto nbits = b - (ndigits * 32);
160  for (int i = n - 1; i >= 0; --i) {
161  std::uint32_t xi = 0;
162  if (i >= int(ndigits)) {
163  xi = x[i - ndigits] << nbits;
164  }
165  // nbits &&, because apparently shifting a 32-bit value by 32 is not allowed
166  if (nbits && (i > int(ndigits))) {
167  xi |= x[i - ndigits - 1] >> (32 - nbits);
168  }
169  x[i] = xi;
170  }
171  return *this;
172 }
173 
174 template <int n>
176  auto ndigits = b / 32;
177  auto nbits = b - (ndigits * 32);
178  for (int i = 0; i < n; ++i) {
179  std::uint32_t xi = 0;
180  if (i + ndigits < n) xi = x[i + ndigits] >> nbits;
181  if (nbits && i + ndigits + 1 < n) xi |= x[i + ndigits + 1] << (32 - nbits);
182  x[i] = xi;
183  }
184  return *this;
185 }
186 
187 template <int n>
188 std::ostream& operator<<(std::ostream& os, BigUInt<n> a) {
189  char buf[n * 20];
190  int i = 0;
191  while (a) {
192  BigUInt<n> quotient;
193  divmod(quotient, a, 10);
194  auto remainder = a[0];
195  a = quotient;
196  buf[i++] = char(remainder) + '0';
197  }
198  for (int j = 0; j < i / 2; ++j) {
199  auto tmp = buf[i - j - 1];
200  buf[i - j - 1] = buf[j];
201  buf[j] = tmp;
202  }
203  if (i == 0) buf[i++] = '0';
204  buf[i] = '\0';
205  os << buf;
206  return os;
207 }
208 
209 template <int n>
211  auto c = a;
212  c += b;
213  return c;
214 }
215 
216 template <int n>
218  auto c = a;
219  c -= b;
220  return c;
221 }
222 
223 template <int n>
225  BigUInt<n> a_times_b_i;
226  BigUInt<n> c{0};
227  for (int i = n - 1; i >= 0; --i) {
228  a_times_b_i = a;
229  a_times_b_i *= b[i];
230  c <<= 32;
231  c += a_times_b_i;
232  }
233  return c;
234 }
235 
236 template <int n>
237 BigUInt<n> operator/(BigUInt<n> const& a, std::uint32_t const& b) {
238  BigUInt<n> quotient;
239  auto x = a;
240  divmod(quotient, x, b);
241  return quotient;
242 }
243 
244 template <int n>
246  if (b > a) return BigUInt<n>(0);
247  BigUInt<n> quotient(1);
248  auto c = b;
249  while (c < a) {
250  c <<= 1;
251  quotient <<= 1;
252  }
253  auto factor = quotient;
254  factor >>= 1;
255  while (factor) {
256  int d = comp(a, c);
257  if (d == 0) break;
258  if (d == -1) {
259  c -= b * factor;
260  quotient -= factor;
261  } else {
262  c += b * factor;
263  quotient += factor;
264  }
265  factor >>= 1;
266  }
267  if (c > a) {
268  c -= b;
269  quotient -= 1;
270  }
271  return quotient;
272 }
273 
274 template <int n>
275 int comp(BigUInt<n> const& a, BigUInt<n> const& b) {
276  for (int i = n - 1; i >= 0; --i) {
277  if (a[i] != b[i]) {
278  if (a[i] > b[i]) return 1;
279  else return -1;
280  }
281  }
282  return 0;
283 }
284 
285 template <int n>
286 bool operator>=(BigUInt<n> const& a, BigUInt<n> const& b) {
287  return comp(a, b) > -1;
288 }
289 
290 template <int n>
291 bool operator<=(BigUInt<n> const& a, BigUInt<n> const& b) {
292  return comp(a, b) < 1;
293 }
294 
295 template <int n>
296 bool operator<(BigUInt<n> const& a, BigUInt<n> const& b) {
297  return comp(a, b) == -1;
298 }
299 
300 template <int n>
301 bool operator>(BigUInt<n> const& a, BigUInt<n> const& b) {
302  return comp(a, b) == 1;
303 }
304 
305 template <int n>
306 bool operator==(BigUInt<n> const& a, BigUInt<n> const& b) {
307  for (int i = 0; i < n; ++i) if (a[i] != b[i]) return false;
308  return true;
309 }
310 
311 template <int n>
312 void divmod(BigUInt<n>& quotient, BigUInt<n>& x, std::uint32_t const& b) {
313  quotient = BigUInt<n>(std::uint32_t(0));
314  for (int i = n - 1; i >= 0;) {
315  if (x[i]) {
316  if (x[i] >= b) {
317  auto dividend = x[i];
318  auto quotient2 = dividend / b;
319  auto remainder = dividend - (quotient2 * b);
320  quotient[i] = quotient2;
321  x[i] = remainder;
322  } else if (i > 0) {
323  auto dividend = std::uint64_t(x[i]);
324  dividend <<= 32;
325  dividend |= x[i - 1];
326  auto quotient2 = dividend / std::uint64_t(b);
327  auto remainder = dividend - (quotient2 * std::uint64_t(b));
328  quotient[i - 1] = std::uint32_t(quotient2);
329  x[i] = 0;
330  x[i - 1] = std::uint32_t(remainder);
331  i = i - 1;
332  } else {
333  break;
334  }
335  } else {
336  i = i - 1;
337  }
338  }
339 }
340 
341 }
342 
343 #endif
344 
BigUInt< n > operator+(BigUInt< n > const &a, BigUInt< n > const &b)
BigUInt< n > operator/(BigUInt< n > const &a, std::uint32_t const &b)
BigUInt< n > operator-(BigUInt< n > const &a, BigUInt< n > const &b)
Arbitrary-precision unsigned integer class.
BigUInt & operator*=(std::uint32_t b)
std::uint32_t & operator[](int i)
BigUInt & operator<<=(std::uint32_t b)
BigUInt< n > operator*(BigUInt< n > const &a, BigUInt< n > const &b)
bool operator>(BigUInt< n > const &a, BigUInt< n > const &b)
bool operator>=(BigUInt< n > const &a, BigUInt< n > const &b)
BigUInt & operator-=(std::uint32_t b)
BigUInt & operator>>=(std::uint32_t b)
void divmod(BigUInt< n > &quotient, BigUInt< n > &x, std::uint32_t const &b)
Arbitrary-precision unsigned integer declaration.
BigUInt & operator+=(std::uint32_t b)
bool operator==(BigUInt< n > const &a, BigUInt< n > const &b)
int comp(BigUInt< n > const &a, BigUInt< n > const &b)