ECC-1416: GRIB: Setting computed key firstSize to 3e-8 sets wrong coded keys

This commit is contained in:
Shahram Najm 2022-06-27 23:09:09 +01:00
parent dce03425ed
commit e6a42b7958
2 changed files with 101 additions and 46 deletions

View File

@ -180,6 +180,79 @@ static double eval_value_factor(int64_t value, int64_t factor)
return (double)value * pow(10.0, -factor);
}
/* Return 0 if no error and set the two outputs 'value' and 'factor'
* value cannot exceed maximum_value and factor cannot exceed maximum_factor
*/
static int get_scaled_value_and_scale_factor_algorithm1(
double input, int64_t maximum_value, int64_t maximum_factor,
int64_t* ret_value, int64_t* ret_factor)
{
int64_t factor = 0;
int64_t value = 0;
factor = floor(log10(maximum_value)) - floor(log10(input < 0 ? -input : input));
if (factor > maximum_factor) {
return GRIB_INTERNAL_ERROR;
}
value = (int64_t)round(input * pow(10, factor));
while ((value % 10 == 0) && (factor > 0)) {
value /= 10;
factor--;
}
if (value >= maximum_value)
return GRIB_INTERNAL_ERROR;
if (factor >= maximum_factor)
return GRIB_INTERNAL_ERROR;
*ret_factor = factor;
*ret_value = value;
return GRIB_SUCCESS;
}
/* Return 0 if no error and set the two outputs 'value' and 'factor'
* value cannot exceed maximum_value and factor cannot exceed maximum_factor
*/
static int get_scaled_value_and_scale_factor_algorithm2(
double input, int64_t maximum_value, int64_t maximum_factor,
int64_t* ret_value, int64_t* ret_factor)
{
int64_t factor = 0, prev_factor = 0;
int64_t value = 0, prev_value = 0;
double exact = input;
const float epsilon = float_epsilon();
int is_negative = 0;
/* Loop until we find a close enough approximation. Keep the last good values */
if (exact < 0) {
is_negative = 1;
exact *= -1;
}
factor = prev_factor = 0;
value = prev_value = round(exact);
while (!is_approximately_equal(exact, eval_value_factor(value, factor), epsilon) &&
value < maximum_value &&
factor < maximum_factor) {
value = round(exact * pow(10., ++factor));
if (value > maximum_value || factor > maximum_factor) {
/* One or more maxima exceeded. So stop and use the previous values */
value = prev_value;
factor = prev_factor;
break;
}
prev_factor = factor;
prev_value = value;
}
if (is_negative) {
value *= -1;
}
*ret_factor = factor;
*ret_value = value;
return GRIB_SUCCESS;
}
static int pack_double(grib_accessor* a, const double* val, size_t* len)
{
/* See ECC-979 */
@ -187,12 +260,10 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
grib_accessor_from_scale_factor_scaled_value* self = (grib_accessor_from_scale_factor_scaled_value*)a;
grib_handle* hand = grib_handle_of_accessor(a);
int ret = 0;
int64_t factor = 0, prev_factor = 0;
int64_t value = 0, prev_value = 0;
int64_t factor = 0;
int64_t value = 0;
double exact = *val; /*the input*/
const float epsilon = float_epsilon();
int is_negative = 0;
unsigned long maxval_value, maxval_factor; /*maximum allowable values*/
int64_t maxval_value, maxval_factor; /*maximum allowable values*/
grib_accessor *accessor_factor, *accessor_value;
if (exact == 0) {
@ -220,29 +291,13 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
maxval_value = (1UL << (accessor_value->length * 8)) - 2; /* exclude missing */
maxval_factor = (1UL << (accessor_factor->length * 8)) - 2; /* exclude missing */
/* Loop until we find a close enough approximation. Keep the last good values */
if (exact < 0) {
is_negative = 1;
exact *= -1;
}
factor = prev_factor = 0;
value = prev_value = round(exact);
while (!is_approximately_equal(exact, eval_value_factor(value, factor), epsilon) &&
value < maxval_value &&
factor < maxval_factor) {
value = round(exact * pow(10., ++factor));
if (value > maxval_value || factor > maxval_factor) {
/* One or more maxima exceeded. So stop and use the previous values */
value = prev_value;
factor = prev_factor;
break;
ret = get_scaled_value_and_scale_factor_algorithm1(exact, maxval_value, maxval_factor, &value, &factor);
if (ret) {
ret = get_scaled_value_and_scale_factor_algorithm2(exact, maxval_value, maxval_factor, &value, &factor);
if (ret) {
grib_context_log(a->context, GRIB_LOG_ERROR, "Failed to compute %s and %s from %g", self->scaleFactor, self->scaledValue, exact);
return ret;
}
prev_factor = factor;
prev_value = value;
}
if (is_negative) {
value *= -1;
}
if ((ret = grib_set_long_internal(hand, self->scaleFactor, factor)) != GRIB_SUCCESS)

View File

@ -88,25 +88,25 @@ ${tools_dir}/grib_set -s lowerLimit=-6.6,upperLimit=-1.02 $tempGrib $temp2
grib_check_key_equals $temp2 scaleFactorOfLowerLimit,scaledValueOfLowerLimit,lowerLimit "1 -66 -6.6"
grib_check_key_equals $temp2 scaleFactorOfUpperLimit,scaledValueOfUpperLimit,upperLimit "2 -102 -1.02"
# input factor value
set_lowerLimit_and_check 550 0 550
set_lowerLimit_and_check -99 0 -99
set_lowerLimit_and_check 6.77 2 677
set_lowerLimit_and_check 0.001 3 1
set_lowerLimit_and_check -6.6 1 -66
set_lowerLimit_and_check -1.02 2 -102
set_lowerLimit_and_check 3e-05 5 3
set_lowerLimit_and_check -3.9e-05 6 -39
#set_lowerLimit_and_check 3.14e-06 6 314
#set_lowerLimit_and_check 3.14e-07 7 314
#set_lowerLimit_and_check 3e-08 8 3
#set_lowerLimit_and_check 1.0e-10 10 1
#set_lowerLimit_and_check 0.03e-06 8 3
#set_lowerLimit_and_check 3.14e-09 11 314
#set_lowerLimit_and_check -3.1456e-09 13 -31456
#set_lowerLimit_and_check 0.0000123 7 123
# input factor value
set_lowerLimit_and_check 0 0 0
set_lowerLimit_and_check 1 0 1
set_lowerLimit_and_check 550 0 550
set_lowerLimit_and_check -99 0 -99
set_lowerLimit_and_check 6.77 2 677
set_lowerLimit_and_check 0.001 3 1
set_lowerLimit_and_check -6.6 1 -66
set_lowerLimit_and_check -1.02 2 -102
set_lowerLimit_and_check 3e-05 5 3
set_lowerLimit_and_check -3.9e-05 6 -39
set_lowerLimit_and_check 3.14e-06 8 314
set_lowerLimit_and_check 3.14e-07 9 314
set_lowerLimit_and_check 3e-08 8 3
set_lowerLimit_and_check -3e-08 8 -3
set_lowerLimit_and_check 1e-10 10 1
set_lowerLimit_and_check 3.14e-09 11 314
set_lowerLimit_and_check 1.23e-05 7 123
set_lowerLimit_and_check -3.1456e-09 13 -31456
# Clean up
rm -f $tempGrib $tempFilt $temp2