OR-Tools  8.1
fp_utils.cc
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 #include "ortools/util/fp_utils.h"
15 
16 #include <cmath>
17 
18 #include "ortools/util/bitset.h"
19 
20 namespace operations_research {
21 
22 namespace {
23 
24 void ReorderAndCapTerms(double* min, double* max) {
25  if (*min > *max) std::swap(*min, *max);
26  if (*min > 0.0) *min = 0.0;
27  if (*max < 0.0) *max = 0.0;
28 }
29 
30 template <bool use_bounds>
31 void ComputeScalingErrors(const std::vector<double>& input,
32  const std::vector<double>& lb,
33  const std::vector<double>& ub, double scaling_factor,
35  double* max_scaled_sum_error) {
36  double max_error = 0.0;
37  double min_error = 0.0;
39  const int size = input.size();
40  for (int i = 0; i < size; ++i) {
41  const double x = input[i];
42  if (x == 0.0) continue;
43  const double scaled = x * scaling_factor;
44 
45  if (scaled == 0.0) {
46  *max_relative_coeff_error = std::numeric_limits<double>::infinity();
47  } else {
49  *max_relative_coeff_error, std::abs(std::round(scaled) / scaled - 1));
50  }
51 
52  const double error = std::round(scaled) - scaled;
53  const double error_lb = (use_bounds ? error * lb[i] : -error);
54  const double error_ub = (use_bounds ? error * ub[i] : error);
55  max_error += std::max(error_lb, error_ub);
56  min_error += std::min(error_lb, error_ub);
57  }
58  *max_scaled_sum_error = std::max(std::abs(max_error), std::abs(min_error));
59 }
60 
61 template <bool use_bounds>
62 void GetBestScalingOfDoublesToInt64(const std::vector<double>& input,
63  const std::vector<double>& lb,
64  const std::vector<double>& ub,
65  int64 max_absolute_sum,
66  double* scaling_factor) {
67  const double kInfinity = std::numeric_limits<double>::infinity();
68 
69  // We start by initializing the returns value to the "error" state.
70  *scaling_factor = 0;
71 
72  // Abort in the "error" state if max_absolute_sum doesn't make sense.
73  if (max_absolute_sum < 0) return;
74 
75  // Our scaling scaling_factor will be 2^factor_exponent.
76  //
77  // TODO(user): Consider using a non-power of two factor if the error can't be
78  // zero? Note however that using power of two has the extra advantage that
79  // subsequent int64 -> double -> scaled back to int64 will loose no extra
80  // information.
81  int factor_exponent = 0;
82  uint64 sum_min = 0; // negated.
83  uint64 sum_max = 0;
84  bool recompute_sum = false;
85  bool is_first_value = true;
86  const int msb = MostSignificantBitPosition64(max_absolute_sum);
87  const int size = input.size();
88  for (int i = 0; i < size; ++i) {
89  const double x = input[i];
90  double min_term = use_bounds ? x * lb[i] : -x;
91  double max_term = use_bounds ? x * ub[i] : x;
92  ReorderAndCapTerms(&min_term, &max_term);
93 
94  // If min_term or max_term is not finite, then abort in the "error" state.
95  if (!(min_term > -kInfinity && max_term < kInfinity)) return;
96 
97  // A value of zero can just be skipped (and needs to because the code below
98  // doesn't handle it correctly).
99  if (min_term == 0.0 && max_term == 0.0) continue;
100 
101  // Compute the greatest candidate such that
102  // round(fabs(c).2^candidate) <= max_absolute_sum.
103  const double c = std::max(-min_term, max_term);
104  int candidate = msb - ilogb(c);
105  if (std::round(ldexp(std::abs(c), candidate)) > max_absolute_sum) {
106  --candidate;
107  }
108  DCHECK_LE(std::abs(static_cast<int64>(round(ldexp(c, candidate)))),
109  max_absolute_sum);
110 
111  // Update factor_exponent which is the min of all the candidates.
112  if (is_first_value || candidate < factor_exponent) {
113  is_first_value = false;
114  factor_exponent = candidate;
115  recompute_sum = true;
116  } else {
117  // Update the sum of absolute values of the numbers seen so far.
118  sum_min -=
119  static_cast<int64>(std::round(ldexp(min_term, factor_exponent)));
120  sum_max +=
121  static_cast<int64>(std::round(ldexp(max_term, factor_exponent)));
122  if (sum_min > max_absolute_sum || sum_max > max_absolute_sum) {
123  factor_exponent--;
124  recompute_sum = true;
125  }
126  }
127 
128  // TODO(user): This is not super efficient, but note that in practice we
129  // will only scan the vector of values about log(size) times. It is possible
130  // to maintain an upper bound on the abs_sum in linear time instead, but the
131  // code and corner cases are a lot more involved. Also, we currently only
132  // use this code in situations where its run-time is negligeable compared to
133  // the rest.
134  while (recompute_sum) {
135  sum_min = 0;
136  sum_max = 0;
137  for (int j = 0; j <= i; ++j) {
138  const double x = input[j];
139  double min_term = use_bounds ? x * lb[j] : -x;
140  double max_term = use_bounds ? x * ub[j] : x;
141  ReorderAndCapTerms(&min_term, &max_term);
142  sum_min -=
143  static_cast<int64>(std::round(ldexp(min_term, factor_exponent)));
144  sum_max +=
145  static_cast<int64>(std::round(ldexp(max_term, factor_exponent)));
146  }
147  if (sum_min > max_absolute_sum || sum_max > max_absolute_sum) {
148  factor_exponent--;
149  } else {
150  recompute_sum = false;
151  }
152  }
153  }
154  *scaling_factor = ldexp(1.0, factor_exponent);
155 }
156 
157 } // namespace
158 
159 void ComputeScalingErrors(const std::vector<double>& input,
160  const std::vector<double>& lb,
161  const std::vector<double>& ub, double scaling_factor,
162  double* max_relative_coeff_error,
163  double* max_scaled_sum_error) {
164  ComputeScalingErrors<true>(input, lb, ub, scaling_factor,
165  max_relative_coeff_error, max_scaled_sum_error);
166 }
167 
168 double GetBestScalingOfDoublesToInt64(const std::vector<double>& input,
169  const std::vector<double>& lb,
170  const std::vector<double>& ub,
171  int64 max_absolute_sum) {
172  double scaling_factor;
173  GetBestScalingOfDoublesToInt64<true>(input, lb, ub, max_absolute_sum,
174  &scaling_factor);
175  return scaling_factor;
176 }
177 
178 void GetBestScalingOfDoublesToInt64(const std::vector<double>& input,
179  int64 max_absolute_sum,
180  double* scaling_factor,
181  double* max_relative_coeff_error) {
182  double max_scaled_sum_error;
183  GetBestScalingOfDoublesToInt64<false>(input, {}, {}, max_absolute_sum,
184  scaling_factor);
185  ComputeScalingErrors<false>(input, {}, {}, *scaling_factor,
186  max_relative_coeff_error, &max_scaled_sum_error);
187 }
188 
189 int64 ComputeGcdOfRoundedDoubles(const std::vector<double>& x,
190  double scaling_factor) {
191  int64 gcd = 0;
192  for (int i = 0; i < x.size() && gcd != 1; ++i) {
193  int64 value = std::abs(std::round(x[i] * scaling_factor));
194  DCHECK_GE(value, 0);
195  if (value == 0) continue;
196  if (gcd == 0) {
197  gcd = value;
198  continue;
199  }
200  // GCD(gcd, value) = GCD(value, gcd % value);
201  while (value != 0) {
202  const int64 r = gcd % value;
203  gcd = value;
204  value = r;
205  }
206  }
207  DCHECK_GE(gcd, 0);
208  return gcd > 0 ? gcd : 1;
209 }
210 
211 } // namespace operations_research
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
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
int64
int64_t int64
Definition: integral_types.h:34
operations_research::glop::kInfinity
const double kInfinity
Definition: lp_types.h:83
fp_utils.h
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
uint64
uint64_t uint64
Definition: integral_types.h:39
operations_research::MostSignificantBitPosition64
int MostSignificantBitPosition64(uint64 n)
Definition: bitset.h:231
DCHECK_GE
#define DCHECK_GE(val1, val2)
Definition: base/logging.h:889
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)
max_relative_coeff_error
double max_relative_coeff_error
Definition: sat/lp_utils.cc:490
DCHECK_LE
#define DCHECK_LE(val1, val2)
Definition: base/logging.h:887
bitset.h