1 // Copyright 2010-2021 Google LLC
2 // Licensed under the Apache License, Version 2.0 (the "License");
3 // you may not use this file except in compliance with the License.
4 // You may obtain a copy of the License at
5 //
6 // http://www.apache.org/licenses/LICENSE-2.0
7 //
8 // Unless required by applicable law or agreed to in writing, software
9 // distributed under the License is distributed on an "AS IS" BASIS,
10 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11 // See the License for the specific language governing permissions and
12 // limitations under the License.
13
14 #ifndef OR_TOOLS_UTIL_SATURATED_ARITHMETIC_H_
15 #define OR_TOOLS_UTIL_SATURATED_ARITHMETIC_H_
16
17 #include "absl/base/casts.h"
18 #include "ortools/base/integral_types.h"
19 #include "ortools/util/bitset.h"
20
21 namespace operations_research {
22 // ---------- Overflow utility functions ----------
23
24 // Implement two's complement addition and subtraction on int64s.
25 //
26 // The C and C++ standards specify that the overflow of signed integers is
27 // undefined. This is because of the different possible representations that may
28 // be used for signed integers (one's complement, two's complement, sign and
29 // magnitude). Such overflows are detected by Address Sanitizer with
30 // -fsanitize=signed-integer-overflow.
31 //
32 // Simple, portable overflow detection on current machines relies on
33 // these two functions. For example, if the sign of the sum of two positive
34 // integers is negative, there has been an overflow.
35 //
36 // Note that the static assert will break if the code is compiled on machines
37 // which do not use two's complement.
TwosComplementAddition(int64_t x,int64_t y)38 inline int64_t TwosComplementAddition(int64_t x, int64_t y) {
39 static_assert(static_cast<uint64_t>(-1LL) == ~0ULL,
40 "The target architecture does not use two's complement.");
41 return absl::bit_cast<int64_t>(static_cast<uint64_t>(x) +
42 static_cast<uint64_t>(y));
43 }
44
TwosComplementSubtraction(int64_t x,int64_t y)45 inline int64_t TwosComplementSubtraction(int64_t x, int64_t y) {
46 static_assert(static_cast<uint64_t>(-1LL) == ~0ULL,
47 "The target architecture does not use two's complement.");
48 return absl::bit_cast<int64_t>(static_cast<uint64_t>(x) -
49 static_cast<uint64_t>(y));
50 }
51
52 // Helper function that returns true if an overflow has occurred in computing
53 // sum = x + y. sum is expected to be computed elsewhere.
AddHadOverflow(int64_t x,int64_t y,int64_t sum)54 inline bool AddHadOverflow(int64_t x, int64_t y, int64_t sum) {
55 // Overflow cannot occur if operands have different signs.
56 // It can only occur if sign(x) == sign(y) and sign(sum) != sign(x),
57 // which is equivalent to: sign(x) != sign(sum) && sign(y) != sign(sum).
58 // This is captured when the expression below is negative.
59 DCHECK_EQ(sum, TwosComplementAddition(x, y));
60 return ((x ^ sum) & (y ^ sum)) < 0;
61 }
62
SubHadOverflow(int64_t x,int64_t y,int64_t diff)63 inline bool SubHadOverflow(int64_t x, int64_t y, int64_t diff) {
64 // This is the same reasoning as for AddHadOverflow. We have x = diff + y.
65 // The formula is the same, with 'x' and diff exchanged.
66 DCHECK_EQ(diff, TwosComplementSubtraction(x, y));
67 return AddHadOverflow(diff, y, x);
68 }
69
70 // A note on overflow treatment.
71 // kint64min and kint64max are treated as infinity.
72 // Thus if the computation overflows, the result is always kint64m(ax/in).
73 //
74 // Note(user): this is actually wrong: when computing A-B, if A is kint64max
75 // and B is finite, then A-B won't be kint64max: overflows aren't sticky.
76 // TODO(user): consider making some operations overflow-sticky, some others
77 // not, but make an explicit choice throughout.
AddOverflows(int64_t x,int64_t y)78 inline bool AddOverflows(int64_t x, int64_t y) {
79 return AddHadOverflow(x, y, TwosComplementAddition(x, y));
80 }
81
SubOverflows(int64_t x,int64_t y)82 inline int64_t SubOverflows(int64_t x, int64_t y) {
83 return SubHadOverflow(x, y, TwosComplementSubtraction(x, y));
84 }
85
86 // Performs *b += a and returns false iff the addition overflow or underflow.
87 // This function only works for typed integer type (IntType<>).
88 template <typename IntegerType>
SafeAddInto(IntegerType a,IntegerType * b)89 bool SafeAddInto(IntegerType a, IntegerType* b) {
90 const int64_t x = a.value();
91 const int64_t y = b->value();
92 const int64_t sum = TwosComplementAddition(x, y);
93 if (AddHadOverflow(x, y, sum)) return false;
94 *b = sum;
95 return true;
96 }
97
98 // Returns kint64max if x >= 0 and kint64min if x < 0.
CapWithSignOf(int64_t x)99 inline int64_t CapWithSignOf(int64_t x) {
100 // return kint64max if x >= 0 or kint64max + 1 (== kint64min) if x < 0.
101 return TwosComplementAddition(kint64max, static_cast<int64_t>(x < 0));
102 }
103
CapAddGeneric(int64_t x,int64_t y)104 inline int64_t CapAddGeneric(int64_t x, int64_t y) {
105 const int64_t result = TwosComplementAddition(x, y);
106 return AddHadOverflow(x, y, result) ? CapWithSignOf(x) : result;
107 }
108
109 #if defined(__GNUC__) && defined(__x86_64__)
110 // TODO(user): port this to other architectures.
CapAddFast(int64_t x,int64_t y)111 inline int64_t CapAddFast(int64_t x, int64_t y) {
112 const int64_t cap = CapWithSignOf(x);
113 int64_t result = x;
114 // clang-format off
115 asm volatile( // 'volatile': ask compiler optimizer "keep as is".
116 "\t" "addq %[y],%[result]"
117 "\n\t" "cmovoq %[cap],%[result]" // Conditional move if overflow.
118 : [result] "=r"(result) // Output
119 : "[result]" (result), [y] "r"(y), [cap] "r"(cap) // Input.
120 : "cc" /* Clobbered registers */ );
121 // clang-format on
122 return result;
123 }
124 #endif
125
CapAdd(int64_t x,int64_t y)126 inline int64_t CapAdd(int64_t x, int64_t y) {
127 #if defined(__GNUC__) && defined(__x86_64__)
128 return CapAddFast(x, y);
129 #else
130 return CapAddGeneric(x, y);
131 #endif
132 }
133
CapAddTo(int64_t x,int64_t * y)134 inline void CapAddTo(int64_t x, int64_t* y) { *y = CapAdd(*y, x); }
135
CapSubGeneric(int64_t x,int64_t y)136 inline int64_t CapSubGeneric(int64_t x, int64_t y) {
137 const int64_t result = TwosComplementSubtraction(x, y);
138 return SubHadOverflow(x, y, result) ? CapWithSignOf(x) : result;
139 }
140
141 #if defined(__GNUC__) && defined(__x86_64__)
142 // TODO(user): port this to other architectures.
CapSubFast(int64_t x,int64_t y)143 inline int64_t CapSubFast(int64_t x, int64_t y) {
144 const int64_t cap = CapWithSignOf(x);
145 int64_t result = x;
146 // clang-format off
147 asm volatile( // 'volatile': ask compiler optimizer "keep as is".
148 "\t" "subq %[y],%[result]"
149 "\n\t" "cmovoq %[cap],%[result]" // Conditional move if overflow.
150 : [result] "=r"(result) // Output
151 : "[result]" (result), [y] "r"(y), [cap] "r"(cap) // Input.
152 : "cc" /* Clobbered registers */ );
153 // clang-format on
154 return result;
155 }
156 #endif
157
CapSub(int64_t x,int64_t y)158 inline int64_t CapSub(int64_t x, int64_t y) {
159 #if defined(__GNUC__) && defined(__x86_64__)
160 return CapSubFast(x, y);
161 #else
162 return CapSubGeneric(x, y);
163 #endif
164 }
165
166 // Note(user): -kint64min != kint64max, but kint64max == ~kint64min.
CapOpp(int64_t v)167 inline int64_t CapOpp(int64_t v) { return v == kint64min ? ~v : -v; }
168
169 namespace cap_prod_util {
170 // Returns an unsigned int equal to the absolute value of n, in a way that
171 // will not produce overflows.
uint_abs(int64_t n)172 inline uint64_t uint_abs(int64_t n) {
173 return n < 0 ? ~static_cast<uint64_t>(n) + 1 : static_cast<uint64_t>(n);
174 }
175 } // namespace cap_prod_util
176
177 // The generic algorithm computes a bound on the number of bits necessary to
178 // store the result. For this it uses the position of the most significant bits
179 // of each of the arguments.
180 // If the result needs at least 64 bits, then return a capped value.
181 // If the result needs at most 63 bits, then return the product.
182 // Otherwise, the result may use 63 or 64 bits: compute the product
183 // as a uint64_t, and cap it if necessary.
CapProdGeneric(int64_t x,int64_t y)184 inline int64_t CapProdGeneric(int64_t x, int64_t y) {
185 const uint64_t a = cap_prod_util::uint_abs(x);
186 const uint64_t b = cap_prod_util::uint_abs(y);
187 // Let MSB(x) denote the most significant bit of x. We have:
188 // MSB(x) + MSB(y) <= MSB(x * y) <= MSB(x) + MSB(y) + 1
189 const int msb_sum =
190 MostSignificantBitPosition64(a) + MostSignificantBitPosition64(b);
191 const int kMaxBitIndexInInt64 = 63;
192 if (msb_sum <= kMaxBitIndexInInt64 - 2) return x * y;
193 // Catch a == 0 or b == 0 now, as MostSignificantBitPosition64(0) == 0.
194 // TODO(user): avoid this by writing function Log2(a) with Log2(0) == -1.
195 if (a == 0 || b == 0) return 0;
196 const int64_t cap = CapWithSignOf(x ^ y);
197 if (msb_sum >= kMaxBitIndexInInt64) return cap;
198 // The corner case is when msb_sum == 62, i.e. at least 63 bits will be
199 // needed to store the product. The following product will never overflow
200 // on uint64_t, since msb_sum == 62.
201 const uint64_t u_prod = a * b;
202 // The overflow cases are captured by one of the following conditions:
203 // (cap >= 0 && u_prod >= static_cast<uint64_t>(kint64max) or
204 // (cap < 0 && u_prod >= static_cast<uint64_t>(kint64min)).
205 // These can be optimized as follows (and if the condition is false, it is
206 // safe to compute x * y.
207 if (u_prod >= static_cast<uint64_t>(cap)) return cap;
208 const int64_t abs_result = absl::bit_cast<int64_t>(u_prod);
209 return cap < 0 ? -abs_result : abs_result;
210 }
211
212 #if defined(__GNUC__) && defined(__x86_64__)
213 // TODO(user): port this to other architectures.
CapProdFast(int64_t x,int64_t y)214 inline int64_t CapProdFast(int64_t x, int64_t y) {
215 // cap = kint64max if x and y have the same sign, cap = kint64min
216 // otherwise.
217 const int64_t cap = CapWithSignOf(x ^ y);
218 int64_t result = x;
219 // Here, we use the fact that imul of two signed 64-integers returns a 128-bit
220 // result -- we care about the lower 64 bits. More importantly, imul also sets
221 // the carry flag if 64 bits were not enough.
222 // We therefore use cmovc to return cap if the carry was set.
223 // clang-format off
224 asm volatile( // 'volatile': ask compiler optimizer "keep as is".
225 "\n\t" "imulq %[y],%[result]"
226 "\n\t" "cmovcq %[cap],%[result]" // Conditional move if carry.
227 : [result] "=r"(result) // Output
228 : "[result]" (result), [y] "r"(y), [cap] "r"(cap) // Input.
229 : "cc" /* Clobbered registers */);
230 // clang-format on
231 return result;
232 }
233 #endif
234
CapProd(int64_t x,int64_t y)235 inline int64_t CapProd(int64_t x, int64_t y) {
236 #if defined(__GNUC__) && defined(__x86_64__)
237 return CapProdFast(x, y);
238 #else
239 return CapProdGeneric(x, y);
240 #endif
241 }
242 } // namespace operations_research
243
244 #endif // OR_TOOLS_UTIL_SATURATED_ARITHMETIC_H_
245