Google OR-Tools v9.14
a fast and portable software suite for combinatorial optimization
Loading...
Searching...
No Matches
binary_search.h
Go to the documentation of this file.
1// Copyright 2010-2025 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_ALGORITHMS_BINARY_SEARCH_H_
15#define OR_TOOLS_ALGORITHMS_BINARY_SEARCH_H_
16
17#include <cmath>
18#include <cstdint>
19#include <functional>
20#include <utility>
21
22#include "absl/base/log_severity.h"
23#include "absl/functional/function_ref.h"
24#include "absl/log/check.h"
25#include "absl/numeric/int128.h"
26#include "absl/types/span.h"
28
29namespace operations_research {
30// Finds a point in [x_true, x_false) where f changes from true to false.
31// If check_bounds is true, it will CHECK that f(x_true) = true and
32// f(x_false) = false.
33//
34// EXAMPLE:
35// // Finds the value x in [0,Pi/2] such that cos(x)=2sin(x):
36// BinarySearch<double>(/*x_true=*/0.0, /*x_false=*/M_PI/2,
37// [](double x) { return cos(x) >= 2*sin(x); });
38//
39// MONOTONIC FUNCTIONS: Suppose f is a monotonic boolean function. See below for
40// the NON-MONOTONIC case.
41//
42// If x_true < x_false, this returns X such that:
43// - x_true < X < x_false,
44// - f((x_true, X]) = true (i.e. for all x in (x_true, X], f(x) = true),
45// - f((X, x_false)) = false (i.e. for all x in (X, x_false), f(x) = false)
46// or this returns x_true if such an X does not exist.
47//
48// If x_true > x_false, this function returns X such that:
49// - x_false < X < x_true
50// - f((x_false, X)) = false
51// - f([X, x_true)) = true
52// or this return x_true if such an X does not exist.
53//
54// Also note that 'Point' may be floating-point types: the function will still
55// converge when the midpoint can't be distinguished from one of the limits,
56// which will always happen.
57// You can use other types than numerical types, too. absl::Duration is
58// naturally supported.
59//
60// OVERFLOWS and NON-NUMERICAL TYPES: If your points may be subject to overflow
61// (e.g. ((kint32max-1) + (kint32max))/2 will overflow an int32_t), or they
62// don't support doing (x+y)/2, you can specialize BinarySearchMidpoint() to
63// suit your needs. See the examples in the unit test.
64//
65// NON-MONOTONIC FUNCTIONS: If f isn't monotonic, the binary search will still
66// run with it typical complexity, and finish. The Point X that it returns
67// will be a "local" inflexion point, meaning that the smallest possible move
68// of that point X to a point X' (in the x_true->x_false direction deduced from
69// the arguments) would make f(X') return false. EXAMPLES:
70// - If Point=int32_t, then the returned X verifies f(X) = true and f(X') =
71// false
72// with X' = X + (x_true < x_false ? 1 : -1).
73// - If Point=double, ditto with X' = X * (1 + (x_true < x_false ? ε : -ε)),
74// where ε = std::numeric_limits<double>::epsilon().
75//
76// Note also that even if f() is non-deterministic, i.e. f(X) can sometimes
77// return true and sometimes false for the same X, then the binary search will
78// still finish, but it's hard to say anything about the returned X.
79template <class Point, bool check_bounds = DEBUG_MODE>
80Point BinarySearch(Point x_true, Point x_false, std::function<bool(Point)> f);
81
82// Used by BinarySearch(). This is just (x+y)/2, with some DCHECKs to catch
83// overflows. You should override (i.e. specialize) it if you risk getting
84// overflows or need something more complicated than (x+y)/2.
85// See examples in the unit test.
86template <class Point>
87Point BinarySearchMidpoint(Point x, Point y);
88
89// Returns the minimum of a convex function on a discrete set of sorted points.
90// It is an error to call this with an empty set of points.
91//
92// We assume the function is "unimodal" with potentially more than one minimum.
93// That is strictly decreasing, then have a minimum that can have many points,
94// and then is strictly increasing. In this case if we have two points with
95// exactly the same value, one of the minimum is always between the two. We will
96// only return one of the minimum.
97//
98// Note that if we allow for non strictly decreasing/increasing, then you have
99// corner cases where one need to check all points to find the minimum. For
100// instance, if the function is constant except at one point where it is lower.
101//
102// The usual algorithm to optimize such a function is a ternary search. However
103// here we assume calls to f() are expensive, and we try to minimize those. So
104// we use a slightly different algorithm than:
105// https://en.wikipedia.org/wiki/Ternary_search
106//
107// TODO(user): Some relevant optimizations:
108// - Abort early if we know a lower bound on the min.
109// - Seed with a starting point if we know one.
110// - We technically do not need the points to be sorted and can use
111// linear-time median computation to speed this up.
112//
113// TODO(user): replace std::function by absl::FunctionRef here and in
114// BinarySearch().
115template <class Point, class Value>
116std::pair<Point, Value> ConvexMinimum(absl::Span<const Point> sorted_points,
117 std::function<Value(Point)> f);
118
119// Internal part of ConvexMinimum() that can also be used directly in some
120// situation when we already know some value of f(). This assumes that we
121// already have a current_min candidate that is either before or after all the
122// points in sorted_points.
123template <class Point, class Value>
124std::pair<Point, Value> ConvexMinimum(bool is_to_the_right,
125 std::pair<Point, Value> current_min,
126 absl::Span<const Point> sorted_points,
127 std::function<Value(Point)> f);
128
129// Searches in the range [begin, end), where Point supports basic arithmetic.
130template <class Point, class Value>
131std::pair<Point, Value> RangeConvexMinimum(Point begin, Point end,
132 absl::FunctionRef<Value(Point)> f);
133template <class Point, class Value>
134std::pair<Point, Value> RangeConvexMinimum(std::pair<Point, Value> current_min,
135 Point begin, Point end,
136 absl::FunctionRef<Value(Point)> f);
137
138// _____________________________________________________________________________
139// Implementation.
140
141namespace internal {
142template <class T>
143bool IsNanGeneric(const T&) {
144 return false;
145}
146
147template <>
148inline bool IsNanGeneric(const float& x) {
149 return std::isnan(x);
150}
151template <>
152inline bool IsNanGeneric(const double& x) {
153 return std::isnan(x);
154}
155template <>
156inline bool IsNanGeneric(const long double& x) {
157 return std::isnan(x);
158}
159
160template <typename T>
161bool AreNumbersOfSameSign(const T& x, const T& y) {
162 if constexpr (std::is_same_v<T, absl::int128>) {
163 return (x < 0) == (y < 0);
164 } else if constexpr (std::is_same_v<T, absl::uint128> ||
165 std::is_unsigned_v<T>) {
166 return true;
167 } else if constexpr (std::is_signed_v<T>) {
168 return std::signbit(x) == std::signbit(y);
169 }
170 return false;
171}
172
173template <typename IntT>
174constexpr bool UsesTwosComplement() {
175 if constexpr (std::is_integral_v<IntT>) {
176 if constexpr (std::is_unsigned_v<IntT>) {
177 return true;
178 } else {
179 return (IntT{-1} & IntT{3}) == IntT{3};
180 }
181 }
182 return false;
183}
184
185} // namespace internal
186
187template <class Point>
188Point BinarySearchMidpoint(Point x, Point y) {
189 Point midpoint;
190 if constexpr (internal::UsesTwosComplement<Point>()) {
191 // For integers using two's complement (all compilers in practice up to
192 // c++17, and all compilers in theory starting from c++20), we can use a
193 // trick from Hacker's delight. See e.g.
194 // https://lemire.me/blog/2022/12/06/fast-midpoint-between-two-integers-without-overflow/
195 midpoint = (x | y) - ((x ^ y) >> 1);
196 } else {
197 midpoint = !internal::AreNumbersOfSameSign(x, y)
198 ? (x + y) / 2
199 :
200 // For integers of the same sign, avoid overflows with a
201 // simple trick.
202 x < y ? x + (y - x) / 2
203 : y + (x - y) / 2;
204 }
205 if (DEBUG_MODE &&
206 !internal::IsNanGeneric(midpoint)) { // Check for overflows.
207 if (x < y) {
208 DCHECK_LE(midpoint, y) << "Overflow?";
209 DCHECK_GE(midpoint, x) << "Overflow?";
210 } else {
211 DCHECK_GE(midpoint, y) << "Overflow?";
212 DCHECK_LE(midpoint, x) << "Overflow?";
213 }
214 }
215 return midpoint;
216}
217
218template <class Point, bool check_bounds>
219Point BinarySearch(Point x_true, Point x_false, std::function<bool(Point)> f) {
220 if constexpr (check_bounds) {
221 CHECK(f(x_true)) << x_true;
222 CHECK(!f(x_false)) << x_false;
223 }
224 int num_iterations = 0;
225 constexpr int kMaxNumIterations = 1000000;
226 while (true) {
227 // NOTE(user): If your "Point" type doesn't support the + or the /2
228 // operations, we could imagine using the same trick as IsNanGeneric() to
229 // make BinarySearch() work for your use case.
230 const Point x = BinarySearchMidpoint(x_true, x_false);
231 if (internal::IsNanGeneric(x) || x == x_true || x == x_false) return x_true;
232 if (++num_iterations > kMaxNumIterations) {
233 LOG(DFATAL)
234 << "The binary search seems to loop forever! This indicates that your"
235 " input types don't converge when repeatedly taking the midpoint";
236 return x_true;
237 }
238 if (f(x)) {
239 x_true = x;
240 } else {
241 x_false = x;
242 }
243 }
244}
245
246template <class Point, class Value>
247std::pair<Point, Value> RangeConvexMinimum(Point begin, Point end,
248 absl::FunctionRef<Value(Point)> f) {
249 DCHECK_LT(begin, end);
250 const Value size = end - begin;
251 if (size == 1) {
252 return {begin, f(begin)};
253 }
254
255 // Starts by splitting interval in two with two queries and getting some info.
256 // Note the current min will be outside the interval.
257 std::pair<Point, Value> current_min;
258 {
259 DCHECK_GE(size, 2);
260 const Point mid = begin + (end - begin) / 2;
261 DCHECK_GT(mid, begin);
262 const Value v = f(mid);
263 const Point before_mid = mid - 1;
264 const Value before_v = f(before_mid);
265 if (before_v == v) return {before_mid, before_v};
266 if (before_v < v) {
267 // Note that we exclude before_mid from the range.
268 current_min = {before_mid, before_v};
269 end = before_mid;
270 } else {
271 current_min = {mid, v};
272 begin = mid + 1;
273 }
274 }
275 if (begin >= end) return current_min;
276 return RangeConvexMinimum<Point, Value>(current_min, begin, end, f);
277}
278
279template <class Point, class Value>
280std::pair<Point, Value> RangeConvexMinimum(std::pair<Point, Value> current_min,
281 Point begin, Point end,
282 absl::FunctionRef<Value(Point)> f) {
283 DCHECK_LT(begin, end);
284 while ((end - begin) > 1) {
285 DCHECK(current_min.first < begin || current_min.first >= end);
286 bool current_is_after_end = current_min.first >= end;
287 const Point mid = begin + (end - begin) / 2;
288 const Value v = f(mid);
289 if (v >= current_min.second) {
290 // If the midpoint is no better than our current minimum, then the
291 // global min must lie between our midpoint and our current min.
292 if (current_is_after_end) {
293 begin = mid + 1;
294 } else {
295 end = mid;
296 }
297 } else {
298 // v < current_min.second, we cannot decide, so we use a second value
299 // close to v like in the initial step.
300 DCHECK_GT(mid, begin);
301 const Point before_mid = mid - 1;
302 const Value before_v = f(before_mid);
303 if (before_v == v) return {before_mid, before_v};
304 if (before_v < v) {
305 current_min = {before_mid, before_v};
306 current_is_after_end = true;
307 end = before_mid;
308 } else {
309 current_is_after_end = false;
310 current_min = {mid, v};
311 begin = mid + 1;
312 }
313 }
314 }
315
316 if (end - begin == 1) {
317 const Value v = f(begin);
318 if (v <= current_min.second) return {begin, v};
319 }
320 return current_min;
321}
322
323template <class Point, class Value>
324std::pair<Point, Value> ConvexMinimum(absl::Span<const Point> sorted_points,
325 std::function<Value(Point)> f) {
326 auto index_f = [&](int index) -> Value { return f(sorted_points[index]); };
327 const auto& [index, v] =
328 RangeConvexMinimum<int64_t, Value>(0, sorted_points.size(), index_f);
329 return {sorted_points[index], v};
330}
331
332template <class Point, class Value>
333std::pair<Point, Value> ConvexMinimum(bool is_to_the_right,
334 std::pair<Point, Value> current_min,
335 absl::Span<const Point> sorted_points,
336 std::function<Value(Point)> f) {
337 auto index_f = [&](int index) -> Value { return f(sorted_points[index]); };
338 std::pair<int, Value> index_current_min = std::make_pair(
339 is_to_the_right ? sorted_points.size() : -1, current_min.second);
340 const auto& [index, v] = RangeConvexMinimum<int64_t, Value>(
341 index_current_min, 0, sorted_points.size(), index_f);
342 if (index == index_current_min.first) return current_min;
343 return {sorted_points[index], v};
344}
345} // namespace operations_research
346
347#endif // OR_TOOLS_ALGORITHMS_BINARY_SEARCH_H_
bool AreNumbersOfSameSign(const T &x, const T &y)
constexpr bool UsesTwosComplement()
std::function< int64_t(const Model &)> Value(IntegerVariable v)
This checks that the variable is fixed.
Definition integer.h:1601
In SWIG mode, we don't want anything besides these top-level includes.
Point BinarySearchMidpoint(Point x, Point y)
ClosedInterval::Iterator end(ClosedInterval interval)
const bool DEBUG_MODE
Definition radix_sort.h:266
std::pair< Point, Value > RangeConvexMinimum(Point begin, Point end, absl::FunctionRef< Value(Point)> f)
Searches in the range [begin, end), where Point supports basic arithmetic.
Point BinarySearch(Point x_true, Point x_false, std::function< bool(Point)> f)
std::pair< Point, Value > ConvexMinimum(absl::Span< const Point > sorted_points, std::function< Value(Point)> f)
ClosedInterval::Iterator begin(ClosedInterval interval)