LSST Applications g04e9c324dd+8c5ae1fdc5,g134cb467dc+1b3060144d,g18429d2f64+f642bf4753,g199a45376c+0ba108daf9,g1fd858c14a+2dcf163641,g262e1987ae+7b8c96d2ca,g29ae962dfc+3bd6ecb08a,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+53e1a9e7c5,g4595892280+fef73a337f,g47891489e3+2efcf17695,g4d44eb3520+642b70b07e,g53246c7159+8c5ae1fdc5,g67b6fd64d1+2efcf17695,g67fd3c3899+b70e05ef52,g74acd417e5+317eb4c7d4,g786e29fd12+668abc6043,g87389fa792+8856018cbb,g89139ef638+2efcf17695,g8d7436a09f+3be3c13596,g8ea07a8fe4+9f5ccc88ac,g90f42f885a+a4e7b16d9b,g97be763408+ad77d7208f,g9dd6db0277+b70e05ef52,ga681d05dcb+a3f46e7fff,gabf8522325+735880ea63,gac2eed3f23+2efcf17695,gb89ab40317+2efcf17695,gbf99507273+8c5ae1fdc5,gd8ff7fe66e+b70e05ef52,gdab6d2f7ff+317eb4c7d4,gdc713202bf+b70e05ef52,gdfd2d52018+b10e285e0f,ge365c994fd+310e8507c4,ge410e46f29+2efcf17695,geaed405ab2+562b3308c0,gffca2db377+8c5ae1fdc5,w.2025.35
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(std::uint32_t * const m,
47 std::uint32_t const * const m1,
48 std::uint32_t const * const m2,
49 unsigned const size1,
50 unsigned const size2)
51{
52 std::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<std::uint64_t>(m1[i]) + (sum >> 32) +
57 static_cast<std::uint64_t>(m2[i]);
58 m[i] = static_cast<std::uint32_t>(sum);
59 }
60 for (; i < size1; ++i) {
61 sum = static_cast<std::uint64_t>(m1[i]) + (sum >> 32);
62 m[i] = static_cast<std::uint32_t>(sum);
63 }
64 for (; i < size2; ++i) {
65 sum = static_cast<std::uint64_t>(m2[i]) + (sum >> 32);
66 m[i] = static_cast<std::uint32_t>(sum);
67 }
68 std::uint32_t carry = static_cast<std::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(std::uint32_t * const m,
80 std::uint32_t const * const m1,
81 std::uint32_t const * const m2,
82 unsigned const size1,
83 unsigned const size2)
84{
85 std::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<std::int64_t>(m1[i]) + (diff >> 32) -
95 static_cast<std::int64_t>(m2[i]);
96 m[i] = static_cast<std::uint32_t>(diff);
97 }
98 // Borrow from the remaining digits in m1 while necessary.
99 for (; diff < 0; ++i) {
100 diff = static_cast<std::int64_t>(m1[i]) - 1;
101 m[i] = static_cast<std::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(std::uint32_t * const m,
120 std::uint32_t const * const m1,
121 std::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 std::uint64_t digit = m1[i - 1];
138 std::uint64_t carry = static_cast<std::uint64_t>(m2[0]) * digit;
139 unsigned j = 1;
140 m[i - 1] = static_cast<std::uint32_t>(carry);
141 carry >>= 32;
142 for (; j < size2; ++j) {
143 carry = static_cast<std::uint64_t>(m2[j]) * digit +
144 static_cast<std::uint64_t>(m[i + j - 1]) +
145 carry;
146 m[i + j - 1] = static_cast<std::uint32_t>(carry);
147 carry >>= 32;
148 }
149 for (; carry != 0; ++j) {
150 carry += static_cast<std::uint64_t>(m[i + j - 1]);
151 m[i + j - 1] = static_cast<std::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 std::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.
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:87
void negate()
negate multiplies this integer by -1.
Definition BigInteger.h:109
BigInteger & add(BigInteger const &b)
add adds b to this integer.
BigInteger(std::uint32_t *digits, unsigned capacity)
This constructor creates a zero-valued integer with the given digit array.
Definition BigInteger.h:55
T max(T... args)
T min(T... args)