Teuchos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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>
47
52namespace Teuchos {
53
54template <int n>
56
57template <int n>
58BigUInt<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
64template <int n>
65BigUInt<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
72template <int n>
73BigUInt<n>::operator bool() const noexcept {
74 for (int i = 0; i < n; ++i) if (x[i]) return true;
75 return false;
76}
77
78template <int n>
79std::uint32_t& BigUInt<n>::operator[](int i) { return x[i]; }
80
81template <int n>
82std::uint32_t const& BigUInt<n>::operator[](int i) const { return x[i]; }
83
84template <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
96template <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
109template <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
126template <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
144template <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
156template <int n>
157BigUInt<n>& BigUInt<n>::operator<<=(std::uint32_t b) {
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
174template <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
187template <int n>
188std::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
209template <int n>
211 auto c = a;
212 c += b;
213 return c;
214}
215
216template <int n>
218 auto c = a;
219 c -= b;
220 return c;
221}
222
223template <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
236template <int n>
237BigUInt<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
244template <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
274template <int n>
275int 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
285template <int n>
286bool operator>=(BigUInt<n> const& a, BigUInt<n> const& b) {
287 return comp(a, b) > -1;
288}
289
290template <int n>
291bool operator<=(BigUInt<n> const& a, BigUInt<n> const& b) {
292 return comp(a, b) < 1;
293}
294
295template <int n>
296bool operator<(BigUInt<n> const& a, BigUInt<n> const& b) {
297 return comp(a, b) == -1;
298}
299
300template <int n>
301bool operator>(BigUInt<n> const& a, BigUInt<n> const& b) {
302 return comp(a, b) == 1;
303}
304
305template <int n>
306bool 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
311template <int n>
312void 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
Arbitrary-precision unsigned integer declaration.
Arbitrary-precision unsigned integer class.
BigUInt & operator-=(std::uint32_t b)
BigUInt & operator+=(std::uint32_t b)
BigUInt & operator>>=(std::uint32_t b)
std::uint32_t & operator[](int i)
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)
BigUInt< n > operator+(BigUInt< n > const &a, BigUInt< n > const &b)
bool operator>=(BigUInt< n > const &a, BigUInt< n > const &b)
std::ostream & operator<<(std::ostream &os, BigUInt< n > a)
int comp(BigUInt< n > const &a, BigUInt< n > const &b)
bool operator<(BigUInt< n > const &a, BigUInt< n > const &b)
BigUInt< n > operator-(BigUInt< n > const &a, BigUInt< n > const &b)
bool operator<=(BigUInt< n > const &a, BigUInt< n > const &b)
BigUInt< n > operator*(BigUInt< n > const &a, BigUInt< n > const &b)
BigUInt< n > operator/(BigUInt< n > const &a, std::uint32_t const &b)
bool operator==(BigUInt< n > const &a, BigUInt< n > const &b)
bool operator>(BigUInt< n > const &a, BigUInt< n > const &b)