26#include "absl/base/casts.h"
27#include "absl/log/check.h"
28#include "absl/types/span.h"
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;
41template <
bool use_bounds>
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) {
53 if (x == 0.0)
continue;
54 const double scaled =
x * scaling_factor;
57 *max_relative_coeff_error = std::numeric_limits<double>::infinity();
59 *max_relative_coeff_error = std::max(
60 *max_relative_coeff_error, std::abs(std::round(scaled) / scaled - 1));
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);
69 *max_scaled_sum_error = std::max(std::abs(max_error), std::abs(min_error));
72template <
bool use_bounds>
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();
84 if (max_absolute_sum < 0)
return;
92 int factor_exponent = 0;
95 bool recompute_sum =
false;
96 bool is_first_value =
true;
98 const int size =
input.size();
99 for (
int i = 0;
i < size; ++
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);
110 if (min_term == 0.0 && max_term == 0.0)
continue;
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;
119 if (std::round(ldexp(std::abs(c), candidate)) > max_absolute_sum) {
122 DCHECK_LE(std::abs(
static_cast<int64_t
>(round(ldexp(c, candidate)))),
126 if (is_first_value || candidate < factor_exponent) {
127 is_first_value =
false;
128 factor_exponent = candidate;
129 recompute_sum =
true;
133 static_cast<int64_t
>(std::round(ldexp(min_term, factor_exponent)));
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)) {
139 recompute_sum =
true;
149 while (recompute_sum) {
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);
158 static_cast<int64_t
>(std::round(ldexp(min_term, factor_exponent)));
160 static_cast<int64_t
>(std::round(ldexp(max_term, factor_exponent)));
162 if (sum_min >
static_cast<uint64_t
>(max_absolute_sum) ||
163 sum_max >
static_cast<uint64_t
>(max_absolute_sum)) {
166 recompute_sum =
false;
170 *scaling_factor = ldexp(1.0, factor_exponent);
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) {
181 max_relative_coeff_error, max_scaled_sum_error);
185 absl::Span<const double> lb,
186 absl::Span<const double> ub,
187 int64_t max_absolute_sum) {
188 double scaling_factor;
191 DCHECK(std::isfinite(scaling_factor));
192 return scaling_factor;
196 int64_t max_absolute_sum,
197 double* scaling_factor,
198 double* max_relative_coeff_error) {
199 double max_scaled_sum_error;
203 max_relative_coeff_error, &max_scaled_sum_error);
204 DCHECK(std::isfinite(*scaling_factor));
208 double scaling_factor) {
209 DCHECK(std::isfinite(scaling_factor));
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));
215 if (value == 0)
continue;
222 const int64_t r = gcd % value;
228 return gcd > 0 ? gcd : 1;
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)
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)
int64_t ComputeGcdOfRoundedDoubles(absl::Span< const double > x, double scaling_factor)
int MostSignificantBitPosition64(uint64_t n)
static int input(yyscan_t yyscanner)