OR-Tools  8.1
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 
39 #include "ortools/base/logging.h"
40 
41 #if defined(_MSC_VER)
42 static inline double isnan(double value) { return _isnan(value); }
43 static inline double round(double value) { return floor(value + 0.5); }
44 #elif defined(__APPLE__) || __GNUC__ >= 5
45 using std::isnan;
46 #endif
47 
48 namespace 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 
106 template <typename FloatType>
107 inline 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.
119 template <typename FloatType>
120 bool 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.
140 template <typename FloatType>
141 bool 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.
152 template <typename FloatType>
153 bool 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.
160 template <typename FloatType>
161 inline 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.
205 void GetBestScalingOfDoublesToInt64(const std::vector<double>& input,
206  int64 max_absolute_sum,
207  double* scaling_factor,
208  double* max_relative_coeff_error);
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]].
214 double 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.
226 void ComputeScalingErrors(const std::vector<double>& input,
227  const std::vector<double>& lb,
228  const std::vector<double>& ub,
229  const double scaling_factor,
230  double* max_relative_coeff_error,
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.
237 int64 ComputeGcdOfRoundedDoubles(const std::vector<double>& x,
238  double scaling_factor);
239 
240 // Returns alpha * x + (1 - alpha) * y.
241 template <typename FloatType>
242 inline 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_
min
int64 min
Definition: alldiff_cst.cc:138
max
int64 max
Definition: alldiff_cst.cc:139
operations_research::ComputeGcdOfRoundedDoubles
int64 ComputeGcdOfRoundedDoubles(const std::vector< double > &x, double scaling_factor)
Definition: fp_utils.cc:189
operations_research::ScopedFloatingPointEnv::ScopedFloatingPointEnv
ScopedFloatingPointEnv()
Definition: fp_utils.h:62
logging.h
operations_research::AreWithinAbsoluteTolerance
bool AreWithinAbsoluteTolerance(FloatType x, FloatType y, FloatType absolute_tolerance)
Definition: fp_utils.h:141
DCHECK_GT
#define DCHECK_GT(val1, val2)
Definition: base/logging.h:890
value
int64 value
Definition: demon_profiler.cc:43
operations_research
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
Definition: dense_doubly_linked_list.h:21
operations_research::IsIntegerWithinTolerance
bool IsIntegerWithinTolerance(FloatType x, FloatType tolerance)
Definition: fp_utils.h:161
int64
int64_t int64
Definition: integral_types.h:34
operations_research::ComputeScalingErrors
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
CHECK_EQ
#define CHECK_EQ(val1, val2)
Definition: base/logging.h:697
operations_research::ScopedFloatingPointEnv::~ScopedFloatingPointEnv
~ScopedFloatingPointEnv()
Definition: fp_utils.h:70
operations_research::Interpolate
FloatType Interpolate(FloatType x, FloatType y, FloatType alpha)
Definition: fp_utils.h:242
operations_research::GetBestScalingOfDoublesToInt64
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
input
static int input(yyscan_t yyscanner)
operations_research::ScopedFloatingPointEnv::EnableExceptions
void EnableExceptions(int excepts)
Definition: fp_utils.h:78
operations_research::IsPositiveOrNegativeInfinity
bool IsPositiveOrNegativeInfinity(FloatType x)
Definition: fp_utils.h:107
max_relative_coeff_error
double max_relative_coeff_error
Definition: sat/lp_utils.cc:490
operations_research::IsSmallerWithinTolerance
bool IsSmallerWithinTolerance(FloatType x, FloatType y, FloatType tolerance)
Definition: fp_utils.h:153
DCHECK_LE
#define DCHECK_LE(val1, val2)
Definition: base/logging.h:887
operations_research::AreWithinAbsoluteOrRelativeTolerances
bool AreWithinAbsoluteOrRelativeTolerances(FloatType x, FloatType y, FloatType relative_tolerance, FloatType absolute_tolerance)
Definition: fp_utils.h:120
operations_research::ScopedFloatingPointEnv
Definition: fp_utils.h:60