LSST Applications g042eb84c57+730a74494b,g04e9c324dd+8c5ae1fdc5,g134cb467dc+1f1e3e7524,g199a45376c+0ba108daf9,g1fd858c14a+fa7d31856b,g210f2d0738+f66ac109ec,g262e1987ae+83a3acc0e5,g29ae962dfc+d856a2cb1f,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+a1e0c9f713,g47891489e3+0d594cb711,g4d44eb3520+c57ec8f3ed,g4d7b6aa1c5+f66ac109ec,g53246c7159+8c5ae1fdc5,g56a1a4eaf3+fd7ad03fde,g64539dfbff+f66ac109ec,g67b6fd64d1+0d594cb711,g67fd3c3899+f66ac109ec,g6985122a63+0d594cb711,g74acd417e5+3098891321,g786e29fd12+668abc6043,g81db2e9a8d+98e2ab9f28,g87389fa792+8856018cbb,g89139ef638+0d594cb711,g8d7436a09f+80fda9ce03,g8ea07a8fe4+760ca7c3fc,g90f42f885a+033b1d468d,g97be763408+a8a29bda4b,g99822b682c+e3ec3c61f9,g9d5c6a246b+0d5dac0c3d,ga41d0fce20+9243b26dd2,gbf99507273+8c5ae1fdc5,gd7ef33dd92+0d594cb711,gdab6d2f7ff+3098891321,ge410e46f29+0d594cb711,geaed405ab2+c4bbc419c6,gf9a733ac38+8c5ae1fdc5,w.2025.38
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)