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