LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
BigInteger.cc
Go to the documentation of this file.
1 /*
2  * LSST Data Management System
3  * Copyright 2014-2015 AURA/LSST.
4  *
5  * This product includes software developed by the
6  * LSST Project (http://www.lsst.org/).
7  *
8  * This program is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the LSST License Statement and
19  * the GNU General Public License along with this program. If not,
20  * see <https://www.lsstcorp.org/LegalNotices/>.
21  */
22 
25 
27 
28 #include <algorithm>
29 
30 
31 namespace lsst {
32 namespace sphgeom {
33 
34 namespace {
35 
36 // `_add` computes the sum of the digits in m1 and m2 and stores the result in
37 // m, which may equal either m1 or m2. It returns the number of digits in the
38 // sum.
39 unsigned _add(uint32_t * const m,
40  uint32_t const * const m1,
41  uint32_t const * const m2,
42  unsigned const size1,
43  unsigned const size2)
44 {
45  uint64_t sum = 0;
46  unsigned i = 0;
47  unsigned const size = std::min(size1, size2);
48  for (; i < size; ++i) {
49  sum = static_cast<uint64_t>(m1[i]) + (sum >> 32) +
50  static_cast<uint64_t>(m2[i]);
51  m[i] = static_cast<uint32_t>(sum);
52  }
53  for (; i < size1; ++i) {
54  sum = static_cast<uint64_t>(m1[i]) + (sum >> 32);
55  m[i] = static_cast<uint32_t>(sum);
56  }
57  for (; i < size2; ++i) {
58  sum = static_cast<uint64_t>(m2[i]) + (sum >> 32);
59  m[i] = static_cast<uint32_t>(sum);
60  }
61  uint32_t carry = static_cast<uint32_t>(sum >> 32);
62  if (carry != 0) {
63  m[i] = carry;
64  ++i;
65  }
66  return i;
67 }
68 
69 // `_sub` computes the difference of the digits in m1 and m2 and stores the
70 // result in m, which may equal either m1 or m2. It assumes that the difference
71 // is non-negative, and returns the the number of digits in the difference.
72 unsigned _sub(uint32_t * const m,
73  uint32_t const * const m1,
74  uint32_t const * const m2,
75  unsigned const size1,
76  unsigned const size2)
77 {
78  int64_t diff = 0;
79  unsigned i = 0;
80  // Note that right shifting a negative integer is implementation defined
81  // (not undefined) behavior. Here we assume arithmetic right shifts for
82  // negative integers.
83  //
84  // TODO(smm): check for this at build time, and provide an alternate
85  // implementation when this assumption does not hold.
86  for (; i < size2; ++i) {
87  diff = static_cast<int64_t>(m1[i]) + (diff >> 32) -
88  static_cast<int64_t>(m2[i]);
89  m[i] = static_cast<uint32_t>(diff);
90  }
91  // Borrow from the remaining digits in m1 while necessary.
92  for (; diff < 0; ++i) {
93  diff = static_cast<int64_t>(m1[i]) - 1;
94  m[i] = static_cast<uint32_t>(diff);
95  }
96  // If we subtracted anything from the most significant digit in m1, strip
97  // any leading zeroes that may have been introduced.
98  if (i == size1) {
99  while (m1[i - 1] == 0) { --i; }
100  return i;
101  }
102  // Otherwise, copy the remaining digits from m1 to m.
103  for (; i < size1; ++i) {
104  m[i] = m1[i];
105  }
106  return i;
107 }
108 
109 // `_mul` computes the product of m1 and m2 and stores the result in m,
110 // which may equal either m1 or m2. It assumes that m1 has at least as many
111 // digits as m2, and returns the the number of digits in the product.
112 unsigned _mul(uint32_t * const m,
113  uint32_t const * const m1,
114  uint32_t const * const m2,
115  unsigned const size1,
116  unsigned const size2)
117 {
118  // The algorithm implemented is long multiplication, and is good for small
119  // digit counts. A requirement for anything fancier involves enough
120  // complexity that switching to GMP is definitely warranted.
121  //
122  // The only small subtlety here is that because m is allowed to alias m1
123  // or m2, we compute products of m2 from most to least significant digit of
124  // m1. This way no as yet unprocessed input digits are overwritten.
125  unsigned const size = size1 + size2;
126  for (unsigned i = size1; i < size; ++i) {
127  m[i] = 0;
128  }
129  for (unsigned i = size1; i > 0; --i) {
130  uint64_t digit = m1[i - 1];
131  uint64_t carry = static_cast<uint64_t>(m2[0]) * digit;
132  unsigned j = 1;
133  m[i - 1] = static_cast<uint32_t>(carry);
134  carry >>= 32;
135  for (; j < size2; ++j) {
136  carry = static_cast<uint64_t>(m2[j]) * digit +
137  static_cast<uint64_t>(m[i + j - 1]) +
138  carry;
139  m[i + j - 1] = static_cast<uint32_t>(carry);
140  carry >>= 32;
141  }
142  for (; carry != 0; ++j) {
143  carry += static_cast<uint64_t>(m[i + j - 1]);
144  m[i + j - 1] = static_cast<uint32_t>(carry);
145  carry >>= 32;
146  }
147  }
148  // The product of m1 and m2 contains size1 + size2 or size1 + size2 - 1
149  // digits.
150  return m[size - 1] == 0 ? size - 1 : size;
151 }
152 
153 } // unnamed namespace
154 
155 
157  if (b._sign == 0) {
158  return *this;
159  }
160  if (_sign == 0) {
161  *this = b;
162  return *this;
163  }
164  if (this == &b) {
165  return multiplyPow2(1);
166  }
167  // When adding two magnitudes, the maximum number of bits in the result is
168  // one greater than the number of bits B in the larger input. When
169  // subtracting them, the maximum result size is B.
170  _checkCapacity(std::max(_size, b._size) + 1);
171  if (_sign == b._sign) {
172  // If the signs of both integers match, add their magnitudes.
173  _size = _add(_digits, _digits, b._digits, _size, b._size);
174  return *this;
175  }
176  // Otherwise, subtract the smaller magnitude from the larger. The sign of
177  // the result is the sign of the input with the larger magnitude.
178  if (_size > b._size) {
179  _size = _sub(_digits, _digits, b._digits, _size, b._size);
180  } else if (_size < b._size) {
181  _size = _sub(_digits, b._digits, _digits, b._size, _size);
182  _sign = b._sign;
183  } else {
184  // Both magnitudes have the same number of digits. Compare and discard
185  // leading common digits until we find a digit position where the
186  // magnitudes differ, or we determine that they are identical.
187  int i = _size;
188  for (; i > 0 && _digits[i - 1] == b._digits[i - 1]; --i) {}
189  if (i == 0) {
190  setToZero();
191  } else if (_digits[i - 1] > b._digits[i - 1]) {
192  _size = _sub(_digits, _digits, b._digits, i, i);
193  } else {
194  _size = _sub(_digits, b._digits, _digits, i, i);
195  _sign = b._sign;
196  }
197  }
198  return *this;
199 }
200 
202  if (this != &b) {
203  // Avoid code duplication by computing a - b = -(-a + b).
204  // This only works if a and b are distinct.
205  negate();
206  add(b);
207  negate();
208  } else {
209  setToZero();
210  }
211  return *this;
212 }
213 
215  if (_sign == 0 || n == 0) {
216  return *this;
217  }
218  // Decompose n into (s, z), where s is a shift by less than 32 bits, and
219  // z is the number of trailing zero digits introduced by the shift.
220  unsigned const z = (n >> 5);
221  unsigned const s = (n & 0x1f);
222  unsigned const size = _size + z;
223  _checkCapacity(size + 1);
224  if (s == 0) {
225  // Right-shifting an unsigned 32 bit integer by 32 bits is undefined
226  // behavior. Avoid that using this special case code.
227  for (unsigned i = _size; i != 0; --i) {
228  _digits[i + z - 1] = _digits[i - 1];
229  }
230  for (unsigned i = z; i != 0; --i) {
231  _digits[i - 1] = 0;
232  }
233  _size = size;
234  } else {
235  uint32_t low, high = 0;
236  for (unsigned i = _size; i != 0; --i, high = low) {
237  low = _digits[i - 1];
238  _digits[i + z] = (high << s) | (low >> (32 - s));
239  }
240  _digits[z] = high << s;
241  for (unsigned i = z; i != 0; --i) {
242  _digits[i] = 0;
243  }
244  _size = (_digits[size] == 0) ? size: size + 1;
245  }
246  return *this;
247 }
248 
250  _sign *= b._sign;
251  if (_sign == 0) {
252  _size = 0;
253  return *this;
254  }
255  _checkCapacity(_size + b._size);
256  if (_size >= b._size) {
257  _size = _mul(_digits, _digits, b._digits, _size, b._size);
258  } else {
259  _size = _mul(_digits, b._digits, _digits, b._size, _size);
260  }
261  return *this;
262 }
263 
264 }} // namespace lsst::sphgeom
This file declares a class for arbitrary precision integers.
double z
Definition: Match.cc:44
int m
Definition: SpanSet.cc:49
table::Key< int > b
BigInteger is an arbitrary precision signed integer class.
Definition: BigInteger.h:45
BigInteger & subtract(BigInteger const &b)
subtract subtracts b from this integer.
Definition: BigInteger.cc:201
BigInteger & multiply(BigInteger const &b)
multiply multiplies this integer by b.
Definition: BigInteger.cc:249
BigInteger & multiplyPow2(unsigned n)
multiplyPow2 multiplies this integer by 2ⁿ.
Definition: BigInteger.cc:214
void setToZero()
setToZero sets this integer to zero.
Definition: BigInteger.h:81
void negate()
negate multiplies this integer by -1.
Definition: BigInteger.h:103
BigInteger & add(BigInteger const &b)
add adds b to this integer.
Definition: BigInteger.cc:156
T max(T... args)
T min(T... args)
A base class for image defects.