OR-Tools  8.2
fp_utils.h
Go to the documentation of this file.
1// Copyright 2010-2018 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// Utility functions on IEEE floating-point numbers.
15// Implemented on float, double, and long double.
16//
17// Also a placeholder for tools controlling and checking FPU rounding modes.
18//
19// IMPORTANT NOTICE: you need to compile your binary with -frounding-math if
20// you want to use rounding modes.
21
22#ifndef OR_TOOLS_UTIL_FP_UTILS_H_
23#define OR_TOOLS_UTIL_FP_UTILS_H_
24
25#if defined(_MSC_VER)
26#pragma fenv_access(on) // NOLINT
27#else
28#include <fenv.h> // NOLINT
29#endif
30
31#ifdef __SSE__
32#include <xmmintrin.h>
33#endif
34
35#include <algorithm>
36#include <cmath>
37#include <limits>
38
40
41#if defined(_MSC_VER)
42static inline double isnan(double value) { return _isnan(value); }
43static inline double round(double value) { return floor(value + 0.5); }
44#elif defined(__APPLE__) || __GNUC__ >= 5
45using std::isnan;
46#endif
47
48namespace operations_research {
49
50// ScopedFloatingPointEnv is used to easily enable Floating-point exceptions.
51// The initial state is automatically restored when the object is deleted.
52//
53// Note(user): For some reason, this causes an FPE exception to be triggered for
54// unknown reasons when compiled in 32 bits. Because of this, we do not turn
55// on FPE exception if __x86_64__ is not defined.
56//
57// TODO(user): Make it work on 32 bits.
58// TODO(user): Make it work on msvc, currently calls to _controlfp crash.
59
61 public:
63#if defined(_MSC_VER)
64 // saved_control_ = _controlfp(0, 0);
65#elif (defined(__GNUC__) || defined(__llvm__)) && defined(__x86_64__)
66 CHECK_EQ(0, fegetenv(&saved_fenv_));
67#endif
68 }
69
71#if defined(_MSC_VER)
72 // CHECK_EQ(saved_control_, _controlfp(saved_control_, 0xFFFFFFFF));
73#elif defined(__x86_64__) && defined(__GLIBC__)
74 CHECK_EQ(0, fesetenv(&saved_fenv_));
75#endif
76 }
77
78 void EnableExceptions(int excepts) {
79#if defined(_MSC_VER)
80 // _controlfp(static_cast<unsigned int>(excepts), _MCW_EM);
81#elif (defined(__GNUC__) || defined(__llvm__)) && defined(__x86_64__) && \
82 !defined(__ANDROID__)
83 CHECK_EQ(0, fegetenv(&fenv_));
84 excepts &= FE_ALL_EXCEPT;
85#if defined(__APPLE__)
86 fenv_.__control &= ~excepts;
87#elif defined(__FreeBSD__)
88 fenv_.__x87.__control &= ~excepts;
89#else // Linux
90 fenv_.__control_word &= ~excepts;
91#endif
92 fenv_.__mxcsr &= ~(excepts << 7);
93 CHECK_EQ(0, fesetenv(&fenv_));
94#endif
95 }
96
97 private:
98#if defined(_MSC_VER)
99 // unsigned int saved_control_;
100#elif (defined(__GNUC__) || defined(__llvm__)) && defined(__x86_64__)
101 fenv_t fenv_;
102 mutable fenv_t saved_fenv_;
103#endif
104};
105
106template <typename FloatType>
107inline bool IsPositiveOrNegativeInfinity(FloatType x) {
108 return x == std::numeric_limits<FloatType>::infinity() ||
109 x == -std::numeric_limits<FloatType>::infinity();
110}
111
112// Tests whether x and y are close to one another using absolute and relative
113// tolerances.
114// Returns true if |x - y| <= a (with a being the absolute_tolerance).
115// The above case is useful for values that are close to zero.
116// Returns true if |x - y| <= max(|x|, |y|) * r. (with r being the relative
117// tolerance.)
118// The cases for infinities are treated separately to avoid generating NaNs.
119template <typename FloatType>
120bool AreWithinAbsoluteOrRelativeTolerances(FloatType x, FloatType y,
121 FloatType relative_tolerance,
122 FloatType absolute_tolerance) {
123 DCHECK_LE(0.0, relative_tolerance);
124 DCHECK_LE(0.0, absolute_tolerance);
125 DCHECK_GT(1.0, relative_tolerance);
127 return x == y;
128 }
129 const FloatType difference = fabs(x - y);
130 if (difference <= absolute_tolerance) {
131 return true;
132 }
133 const FloatType largest_magnitude = std::max(fabs(x), fabs(y));
134 return difference <= largest_magnitude * relative_tolerance;
135}
136
137// Tests whether x and y are close to one another using an absolute tolerance.
138// Returns true if |x - y| <= a (with a being the absolute_tolerance).
139// The cases for infinities are treated separately to avoid generating NaNs.
140template <typename FloatType>
141bool AreWithinAbsoluteTolerance(FloatType x, FloatType y,
142 FloatType absolute_tolerance) {
143 DCHECK_LE(0.0, absolute_tolerance);
145 return x == y;
146 }
147 return fabs(x - y) <= absolute_tolerance;
148}
149
150// Returns true if x is less than y or slighlty greater than y with the given
151// absolute or relative tolerance.
152template <typename FloatType>
153bool IsSmallerWithinTolerance(FloatType x, FloatType y, FloatType tolerance) {
154 if (IsPositiveOrNegativeInfinity(y)) return x <= y;
155 return x <= y + tolerance * std::max(1.0, std::min(std::abs(x), std::abs(y)));
156}
157
158// Returns true if x is within tolerance of any integer. Always returns
159// false for x equal to +/- infinity.
160template <typename FloatType>
161inline bool IsIntegerWithinTolerance(FloatType x, FloatType tolerance) {
162 DCHECK_LE(0.0, tolerance);
163 if (IsPositiveOrNegativeInfinity(x)) return false;
164 return std::abs(x - std::round(x)) <= tolerance;
165}
166
167// Handy alternatives to EXPECT_NEAR(), using relative and absolute tolerance
168// instead of relative tolerance only, and with a proper support for infinity.
169// TODO(user): investigate moving this to ortools/base/ or some other place.
170#define EXPECT_COMPARABLE(expected, obtained, epsilon) \
171 EXPECT_TRUE(operations_research::AreWithinAbsoluteOrRelativeTolerances( \
172 expected, obtained, epsilon, epsilon)) \
173 << obtained << " != expected value " << expected \
174 << " within epsilon = " << epsilon;
175
176#define EXPECT_NOTCOMPARABLE(expected, obtained, epsilon) \
177 EXPECT_FALSE(operations_research::AreWithinAbsoluteOrRelativeTolerances( \
178 expected, obtained, epsilon, epsilon)) \
179 << obtained << " == expected value " << expected \
180 << " within epsilon = " << epsilon;
181
182// Given an array of doubles, this computes a positive scaling factor such that
183// the scaled doubles can then be rounded to integers with little or no loss of
184// precision, and so that the L1 norm of these integers is <= max_sum. More
185// precisely, the following formulas will hold (x[i] is input[i], for brevity):
186// - For all i, |round(factor * x[i]) / factor - x[i]| <= error * |x[i]|
187// - The sum over i of |round(factor * x[i])| <= max_sum.
188//
189// The algorithm tries to minimize "error" (which is the relative error for one
190// coefficient). Note however than in really broken cases, the error might be
191// infinity and the factor zero.
192//
193// Note on the algorithm:
194// - It only uses factors of the form 2^n (i.e. ldexp(1.0, n)) for simplicity.
195// - The error will be zero in many practical instances. For example, if x
196// contains only integers with low magnitude; or if x contains doubles whose
197// exponents cover a small range.
198// - It chooses the factor as high as possible under the given constraints, as
199// a result the numbers produced may be large. To balance this, we recommend
200// to divide the scaled integers by their gcd() which will result in no loss
201// of precision and will help in many practical cases.
202//
203// TODO(user): incorporate the gcd computation here? The issue is that I am
204// not sure if I just do factor /= gcd that round(x * factor) will be the same.
205void GetBestScalingOfDoublesToInt64(const std::vector<double>& input,
206 int64 max_absolute_sum,
207 double* scaling_factor,
209
210// Returns the scaling factor like above with the extra conditions:
211// - The sum over i of min(0, round(factor * x[i])) >= -max_sum.
212// - The sum over i of max(0, round(factor * x[i])) <= max_sum.
213// For any possible values of the x[i] such that x[i] is in [lb[i], ub[i]].
214double GetBestScalingOfDoublesToInt64(const std::vector<double>& input,
215 const std::vector<double>& lb,
216 const std::vector<double>& ub,
217 int64 max_absolute_sum);
218// This computes:
219//
220// The max_relative_coeff_error, which is the maximum over all coeff of
221// |round(factor * x[i]) / (factor * x[i]) - 1|.
222//
223// The max_scaled_sum_error which is a bound on the maximum difference between
224// the exact scaled sum and the rounded one. One needs to divide this by
225// scaling_factor to have the maximum absolute error on the original sum.
226void ComputeScalingErrors(const std::vector<double>& input,
227 const std::vector<double>& lb,
228 const std::vector<double>& ub,
229 const double scaling_factor,
231 double* max_scaled_sum_error);
232
233// Returns the Greatest Common Divisor of the numbers
234// round(fabs(x[i] * scaling_factor)). The numbers 0 are ignored and if they are
235// all zero then the result is 1. Note that round(fabs()) is the same as
236// fabs(round()) since the numbers are rounded away from zero.
237int64 ComputeGcdOfRoundedDoubles(const std::vector<double>& x,
238 double scaling_factor);
239
240// Returns alpha * x + (1 - alpha) * y.
241template <typename FloatType>
242inline FloatType Interpolate(FloatType x, FloatType y, FloatType alpha) {
243 return alpha * x + (1 - alpha) * y;
244}
245
246} // namespace operations_research
247
248#endif // OR_TOOLS_UTIL_FP_UTILS_H_
int64 min
Definition: alldiff_cst.cc:138
int64 max
Definition: alldiff_cst.cc:139
#define DCHECK_LE(val1, val2)
Definition: base/logging.h:887
#define CHECK_EQ(val1, val2)
Definition: base/logging.h:697
#define DCHECK_GT(val1, val2)
Definition: base/logging.h:890
int64 value
int64_t int64
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
bool AreWithinAbsoluteOrRelativeTolerances(FloatType x, FloatType y, FloatType relative_tolerance, FloatType absolute_tolerance)
Definition: fp_utils.h:120
bool IsPositiveOrNegativeInfinity(FloatType x)
Definition: fp_utils.h:107
int64 ComputeGcdOfRoundedDoubles(const std::vector< double > &x, double scaling_factor)
Definition: fp_utils.cc:189
FloatType Interpolate(FloatType x, FloatType y, FloatType alpha)
Definition: fp_utils.h:242
bool AreWithinAbsoluteTolerance(FloatType x, FloatType y, FloatType absolute_tolerance)
Definition: fp_utils.h:141
void ComputeScalingErrors(const std::vector< double > &input, const std::vector< double > &lb, const std::vector< double > &ub, double scaling_factor, double *max_relative_coeff_error, double *max_scaled_sum_error)
Definition: fp_utils.cc:159
bool IsSmallerWithinTolerance(FloatType x, FloatType y, FloatType tolerance)
Definition: fp_utils.h:153
bool IsIntegerWithinTolerance(FloatType x, FloatType tolerance)
Definition: fp_utils.h:161
double GetBestScalingOfDoublesToInt64(const std::vector< double > &input, const std::vector< double > &lb, const std::vector< double > &ub, int64 max_absolute_sum)
Definition: fp_utils.cc:168
static int input(yyscan_t yyscanner)
double max_relative_coeff_error