LSST Applications g2079a07aa2+86d27d4dc4,g2305ad1205+a659bff248,g2bbee38e9b+3c60f8fe34,g337abbeb29+3c60f8fe34,g33d1c0ed96+3c60f8fe34,g3502564af9+d77d6d1350,g3a166c0a6a+3c60f8fe34,g487adcacf7+25d9892218,g4be5004598+d77d6d1350,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+4d81263f9a,g5cd07815a0+980d2b1c3b,g607f77f49a+d77d6d1350,g858d7b2824+d77d6d1350,g88963caddf+83e433e629,g99cad8db69+a4d3c48eeb,g9ddcbc5298+9a081db1e4,ga1e77700b3+bcf1af89ad,ga57fefb910+9a39d7b2d7,gae0086650b+585e252eca,gb065fddaf9+4f9fd82a2c,gb0e22166c9+60f28cb32d,gb363559e06+d84b1d3d07,gb3b7280ab2+4563d032e1,gb4b16eec92+babe958938,gba4ed39666+c2a2e4ac27,gbb8dafda3b+ed6854b564,gc120e1dc64+b72d212f87,gc28159a63d+3c60f8fe34,gc3e9b769f7+921dbcd359,gcf0d15dbbd+9a39d7b2d7,gdaeeff99f8+f9a426f77a,gddc38dedce+585e252eca,ge79ae78c31+3c60f8fe34,w.2024.21
LSST Data Management Base Package
Loading...
Searching...
No Matches
BigInteger.cc
Go to the documentation of this file.
1/*
2 * This file is part of sphgeom.
3 *
4 * Developed for the LSST Data Management System.
5 * This product includes software developed by the LSST Project
6 * (http://www.lsst.org).
7 * See the COPYRIGHT file at the top-level directory of this distribution
8 * for details of code ownership.
9 *
10 * This software is dual licensed under the GNU General Public License and also
11 * under a 3-clause BSD license. Recipients may choose which of these licenses
12 * to use; please see the files gpl-3.0.txt and/or bsd_license.txt,
13 * respectively. If you choose the GPL option then the following text applies
14 * (but note that there is still no warranty even if you opt for BSD instead):
15 *
16 * This program is free software: you can redistribute it and/or modify
17 * it under the terms of the GNU General Public License as published by
18 * the Free Software Foundation, either version 3 of the License, or
19 * (at your option) any later version.
20 *
21 * This program is distributed in the hope that it will be useful,
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 * GNU General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License
27 * along with this program. If not, see <http://www.gnu.org/licenses/>.
28 */
29
32
34
35#include <algorithm>
36
37
38namespace lsst {
39namespace sphgeom {
40
41namespace {
42
43// `_add` computes the sum of the digits in m1 and m2 and stores the result in
44// m, which may equal either m1 or m2. It returns the number of digits in the
45// sum.
46unsigned _add(uint32_t * const m,
47 uint32_t const * const m1,
48 uint32_t const * const m2,
49 unsigned const size1,
50 unsigned const size2)
51{
52 uint64_t sum = 0;
53 unsigned i = 0;
54 unsigned const size = std::min(size1, size2);
55 for (; i < size; ++i) {
56 sum = static_cast<uint64_t>(m1[i]) + (sum >> 32) +
57 static_cast<uint64_t>(m2[i]);
58 m[i] = static_cast<uint32_t>(sum);
59 }
60 for (; i < size1; ++i) {
61 sum = static_cast<uint64_t>(m1[i]) + (sum >> 32);
62 m[i] = static_cast<uint32_t>(sum);
63 }
64 for (; i < size2; ++i) {
65 sum = static_cast<uint64_t>(m2[i]) + (sum >> 32);
66 m[i] = static_cast<uint32_t>(sum);
67 }
68 uint32_t carry = static_cast<uint32_t>(sum >> 32);
69 if (carry != 0) {
70 m[i] = carry;
71 ++i;
72 }
73 return i;
74}
75
76// `_sub` computes the difference of the digits in m1 and m2 and stores the
77// result in m, which may equal either m1 or m2. It assumes that the difference
78// is non-negative, and returns the the number of digits in the difference.
79unsigned _sub(uint32_t * const m,
80 uint32_t const * const m1,
81 uint32_t const * const m2,
82 unsigned const size1,
83 unsigned const size2)
84{
85 int64_t diff = 0;
86 unsigned i = 0;
87 // Note that right shifting a negative integer is implementation defined
88 // (not undefined) behavior. Here we assume arithmetic right shifts for
89 // negative integers.
90 //
91 // TODO(smm): check for this at build time, and provide an alternate
92 // implementation when this assumption does not hold.
93 for (; i < size2; ++i) {
94 diff = static_cast<int64_t>(m1[i]) + (diff >> 32) -
95 static_cast<int64_t>(m2[i]);
96 m[i] = static_cast<uint32_t>(diff);
97 }
98 // Borrow from the remaining digits in m1 while necessary.
99 for (; diff < 0; ++i) {
100 diff = static_cast<int64_t>(m1[i]) - 1;
101 m[i] = static_cast<uint32_t>(diff);
102 }
103 // If we subtracted anything from the most significant digit in m1, strip
104 // any leading zeroes that may have been introduced.
105 if (i == size1) {
106 while (m1[i - 1] == 0) { --i; }
107 return i;
108 }
109 // Otherwise, copy the remaining digits from m1 to m.
110 for (; i < size1; ++i) {
111 m[i] = m1[i];
112 }
113 return i;
114}
115
116// `_mul` computes the product of m1 and m2 and stores the result in m,
117// which may equal either m1 or m2. It assumes that m1 has at least as many
118// digits as m2, and returns the the number of digits in the product.
119unsigned _mul(uint32_t * const m,
120 uint32_t const * const m1,
121 uint32_t const * const m2,
122 unsigned const size1,
123 unsigned const size2)
124{
125 // The algorithm implemented is long multiplication, and is good for small
126 // digit counts. A requirement for anything fancier involves enough
127 // complexity that switching to GMP is definitely warranted.
128 //
129 // The only small subtlety here is that because m is allowed to alias m1
130 // or m2, we compute products of m2 from most to least significant digit of
131 // m1. This way no as yet unprocessed input digits are overwritten.
132 unsigned const size = size1 + size2;
133 for (unsigned i = size1; i < size; ++i) {
134 m[i] = 0;
135 }
136 for (unsigned i = size1; i > 0; --i) {
137 uint64_t digit = m1[i - 1];
138 uint64_t carry = static_cast<uint64_t>(m2[0]) * digit;
139 unsigned j = 1;
140 m[i - 1] = static_cast<uint32_t>(carry);
141 carry >>= 32;
142 for (; j < size2; ++j) {
143 carry = static_cast<uint64_t>(m2[j]) * digit +
144 static_cast<uint64_t>(m[i + j - 1]) +
145 carry;
146 m[i + j - 1] = static_cast<uint32_t>(carry);
147 carry >>= 32;
148 }
149 for (; carry != 0; ++j) {
150 carry += static_cast<uint64_t>(m[i + j - 1]);
151 m[i + j - 1] = static_cast<uint32_t>(carry);
152 carry >>= 32;
153 }
154 }
155 // The product of m1 and m2 contains size1 + size2 or size1 + size2 - 1
156 // digits.
157 return m[size - 1] == 0 ? size - 1 : size;
158}
159
160} // unnamed namespace
161
162
164 if (b._sign == 0) {
165 return *this;
166 }
167 if (_sign == 0) {
168 *this = b;
169 return *this;
170 }
171 if (this == &b) {
172 return multiplyPow2(1);
173 }
174 // When adding two magnitudes, the maximum number of bits in the result is
175 // one greater than the number of bits B in the larger input. When
176 // subtracting them, the maximum result size is B.
177 _checkCapacity(std::max(_size, b._size) + 1);
178 if (_sign == b._sign) {
179 // If the signs of both integers match, add their magnitudes.
180 _size = _add(_digits, _digits, b._digits, _size, b._size);
181 return *this;
182 }
183 // Otherwise, subtract the smaller magnitude from the larger. The sign of
184 // the result is the sign of the input with the larger magnitude.
185 if (_size > b._size) {
186 _size = _sub(_digits, _digits, b._digits, _size, b._size);
187 } else if (_size < b._size) {
188 _size = _sub(_digits, b._digits, _digits, b._size, _size);
189 _sign = b._sign;
190 } else {
191 // Both magnitudes have the same number of digits. Compare and discard
192 // leading common digits until we find a digit position where the
193 // magnitudes differ, or we determine that they are identical.
194 int i = _size;
195 for (; i > 0 && _digits[i - 1] == b._digits[i - 1]; --i) {}
196 if (i == 0) {
197 setToZero();
198 } else if (_digits[i - 1] > b._digits[i - 1]) {
199 _size = _sub(_digits, _digits, b._digits, i, i);
200 } else {
201 _size = _sub(_digits, b._digits, _digits, i, i);
202 _sign = b._sign;
203 }
204 }
205 return *this;
206}
207
209 if (this != &b) {
210 // Avoid code duplication by computing a - b = -(-a + b).
211 // This only works if a and b are distinct.
212 negate();
213 add(b);
214 negate();
215 } else {
216 setToZero();
217 }
218 return *this;
219}
220
222 if (_sign == 0 || n == 0) {
223 return *this;
224 }
225 // Decompose n into (s, z), where s is a shift by less than 32 bits, and
226 // z is the number of trailing zero digits introduced by the shift.
227 unsigned const z = (n >> 5);
228 unsigned const s = (n & 0x1f);
229 unsigned const size = _size + z;
230 _checkCapacity(size + 1);
231 if (s == 0) {
232 // Right-shifting an unsigned 32 bit integer by 32 bits is undefined
233 // behavior. Avoid that using this special case code.
234 for (unsigned i = _size; i != 0; --i) {
235 _digits[i + z - 1] = _digits[i - 1];
236 }
237 for (unsigned i = z; i != 0; --i) {
238 _digits[i - 1] = 0;
239 }
240 _size = size;
241 } else {
242 uint32_t low, high = 0;
243 for (unsigned i = _size; i != 0; --i, high = low) {
244 low = _digits[i - 1];
245 _digits[i + z] = (high << s) | (low >> (32 - s));
246 }
247 _digits[z] = high << s;
248 for (unsigned i = z; i != 0; --i) {
249 _digits[i] = 0;
250 }
251 _size = (_digits[size] == 0) ? size: size + 1;
252 }
253 return *this;
254}
255
257 _sign *= b._sign;
258 if (_sign == 0) {
259 _size = 0;
260 return *this;
261 }
262 _checkCapacity(_size + b._size);
263 if (_size >= b._size) {
264 _size = _mul(_digits, _digits, b._digits, _size, b._size);
265 } else {
266 _size = _mul(_digits, b._digits, _digits, b._size, _size);
267 }
268 return *this;
269}
270
271}} // namespace lsst::sphgeom
This file declares a class for arbitrary precision integers.
double z
Definition Match.cc:44
int m
Definition SpanSet.cc:48
table::Key< int > b
BigInteger is an arbitrary precision signed integer class.
Definition BigInteger.h:52
BigInteger & subtract(BigInteger const &b)
subtract subtracts b from this integer.
BigInteger & multiply(BigInteger const &b)
multiply multiplies this integer by b.
BigInteger & multiplyPow2(unsigned n)
multiplyPow2 multiplies this integer by 2ⁿ.
void setToZero()
setToZero sets this integer to zero.
Definition BigInteger.h:88
void negate()
negate multiplies this integer by -1.
Definition BigInteger.h:110
BigInteger & add(BigInteger const &b)
add adds b to this integer.
T max(T... args)
T min(T... args)