OR-Tools  8.2
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
15
16#include <cmath>
17
18#include "ortools/util/bitset.h"
19
20namespace operations_research {
21
22namespace {
23
24void 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
30template <bool use_bounds>
31void 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
61template <bool use_bounds>
62void 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
159void ComputeScalingErrors(const std::vector<double>& input,
160 const std::vector<double>& lb,
161 const std::vector<double>& ub, double scaling_factor,
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
168double 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
178void 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
189int64 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
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 DCHECK_GE(val1, val2)
Definition: base/logging.h:889
int64 value
int64_t int64
uint64_t uint64
const double kInfinity
Definition: lp_types.h:83
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
int64 ComputeGcdOfRoundedDoubles(const std::vector< double > &x, double scaling_factor)
Definition: fp_utils.cc:189
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
int MostSignificantBitPosition64(uint64 n)
Definition: bitset.h:231
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