core
/
eccodes
mirror of https://github.com/ecmwf/eccodes.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'feature/ECC-1467_grib_power_single_precision' into develop

This commit is contained in:
Shahram Najm 2023-06-22 19:54:11 +01:00
parent 03496f2fed 53a62ee885
commit 72e9d5dc78
36 changed files with 185 additions and 160 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.gitignore vendored
View File

@ -15,6 +15,8 @@ configure
grib_api.spec
grib_api.pc
grib_api_f90.pc
*.analyzerinfo
*.snalyzerinfo
# compiled source #
###################

4
src/deprecated/grib_accessor_class_apply_operators.cc
View File

@ -8,7 +8,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_scaling.h"
#include "grib_api_internal.h"
/*
This is used by make_class.pl
@ -423,7 +423,7 @@ Assert(0);
case 7:
if (Y) {
extraScale = Y;
referenceValueFactor = grib_power(Y, 10);
referenceValueFactor = grib_power<double>(Y, 10);
extraWidth = ((10 * Y) + 2) / 3;
}
else {

3
src/functions.cc
View File

@ -8,6 +8,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_scaling.h"
#include "grib_api_internal.h"
long grib_op_eq(long a, long b)
@ -69,7 +70,7 @@ double grib_op_neg_d(double a)
long grib_op_pow(long a, long b)
{
/* Note: This is actually 'a' to the power 'b' */
return grib_power(b, a);
return codes_power<double>(b, a);
}
long grib_op_add(long a, long b)

10
src/grib_accessor_class_bufr_data_array.cc
View File

@ -7,7 +7,7 @@
* In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_scaling.h"
#include "grib_api_internal.h"
/*
@ -803,7 +803,7 @@ static int encode_double_array(grib_context* c, grib_buffer* buff, long* pos, bu
modifiedReference = bd->reference;
modifiedFactor = bd->factor;
inverseFactor = grib_power(bd->scale, 10);
inverseFactor = codes_power<double>(bd->scale, 10);
modifiedWidth = bd->width;
err = descriptor_get_min_max(bd, modifiedWidth, modifiedReference, modifiedFactor, &minAllowed, &maxAllowed);
@ -945,10 +945,10 @@ static int encode_double_array(grib_context* c, grib_buffer* buff, long* pos, bu
localRange = (max - min) * inverseFactor + 1;
localWidth = ceil(log(localRange) / log(2.0));
lval = round(max * inverseFactor) - reference;
allone = grib_power(localWidth, 2) - 1;
allone = codes_power<double>(localWidth, 2) - 1;
while (allone <= lval) {
localWidth++;
allone = grib_power(localWidth, 2) - 1;
allone = codes_power<double>(localWidth, 2) - 1;
}
if (localWidth == 1)
localWidth++;
@ -3233,7 +3233,7 @@ static int process_elements(grib_accessor* a, int flag, long onlySubset, long st
return err;
}
bd = grib_bufr_descriptor_clone(self->expanded->v[index]);
bd->reference = -grib_power(bd->width, 2);
bd->reference = -codes_power<double>(bd->width, 2);
bd->width++;
err = codec_element(c, self, iss, buffer, data, &pos, index, bd, elementIndex, dval, sval);

3
src/grib_accessor_class_bufr_elements_table.cc
View File

@ -12,6 +12,7 @@
* Enrico Fucile
****************************************/
#include "grib_scaling.h"
#include "grib_api_internal.h"
#if GRIB_PTHREADS
@ -350,7 +351,7 @@ static int bufr_get_from_table(grib_accessor* a, bufr_descriptor* v)
/* ECC-985: Scale and reference are often 0 so we can reduce calls to atol */
v->scale = atol_fast(list[5]);
v->factor = grib_power(-v->scale, 10);
v->factor = codes_power<double>(-v->scale, 10);
v->reference = atol_fast(list[6]);
v->width = atol(list[7]);

9
src/grib_accessor_class_data_2order_packing.cc
View File

@ -8,6 +8,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_scaling.h"
#include "grib_api_internal.h"
/*
@ -621,8 +622,8 @@ static int unpack_double(grib_accessor* a, double* val, size_t* len)
if (boustrophedonic)
reverse_rows(sec_val, n_vals, Ni, bitmap, bitmap_len);
s = grib_power(binary_scale_factor, 2);
d = grib_power(-decimal_scale_factor, 10);
s = codes_power<double>(binary_scale_factor, 2);
d = codes_power<double>(-decimal_scale_factor, 10);
for (i = 0; i < n_vals; i++)
val[i] = (double)((((double)sec_val[i]) * s) + reference_value) * d;
@ -764,7 +765,7 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
/* calculation of integer array */
sec_val = (unsigned long*)grib_context_malloc(a->context, (n_vals) * sizeof(long));
d = grib_power(decimal_scale_factor, 10);
d = codes_power<double>(decimal_scale_factor, 10);
max = val[0];
min = max;
for (i = 0; i < n_vals; i++) {
@ -786,7 +787,7 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
/* the scale factor in Grib 1 is adjusted in gribex, for "normalization purpose" ... ?*/
binary_scale_factor = grib_get_binary_scale_fact(max, reference_value, bits_per_value, &err);
divisor = grib_power(-binary_scale_factor, 2);
divisor = codes_power<double>(-binary_scale_factor, 2);
for (i = 0; i < n_vals; i++)

2
src/grib_accessor_class_data_apply_bitmap.cc
View File

@ -8,7 +8,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_api_internal_cpp.h"
#include "grib_value.h"
/*
This is used by make_class.pl

19
src/grib_accessor_class_data_ccsds_packing.cc
View File

@ -8,7 +8,8 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_api_internal_cpp.h"
#include "grib_bits_any_endian_simple.h"
#include "grib_scaling.h"
#include <cstdint>
/*
@ -315,7 +316,7 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
return err;
if (bits_per_value == 0 || (binary_scale_factor == 0 && decimal_scale_factor != 0)) {
d = grib_power(decimal_scale_factor, 10);
d = codes_power<double>(decimal_scale_factor, 10);
min *= d;
max *= d;
@ -345,9 +346,9 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
range = (max - min);
unscaled_min = min;
unscaled_max = max;
f = (grib_power(bits_per_value, 2) - 1);
minrange = grib_power(-last, 2) * f;
maxrange = grib_power(last, 2) * f;
f = (codes_power<double>(bits_per_value, 2) - 1);
minrange = codes_power<double>(-last, 2) * f;
maxrange = codes_power<double>(last, 2) * f;
while (range < minrange) {
decimal_scale_factor += 1;
@ -369,11 +370,11 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
"%s %s: unable to find nearest_smaller_value of %g for %s", cclass_name, __func__, min, self->reference_value);
return GRIB_INTERNAL_ERROR;
}
d = grib_power(decimal_scale_factor, 10);
d = codes_power<double>(decimal_scale_factor, 10);
}
binary_scale_factor = grib_get_binary_scale_fact(max, reference_value, bits_per_value, &err);
divisor = grib_power(-binary_scale_factor, 2);
divisor = codes_power<double>(-binary_scale_factor, 2);
size_t nbytes = (bits_per_value + 7) / 8;
// ECC-1602: use native a data type (4 bytes for uint32_t) for values that require only 3 bytes
@ -565,8 +566,8 @@ static int unpack(grib_accessor* a, T* val, size_t* len)
return GRIB_SUCCESS;
}
bscale = grib_power(binary_scale_factor, 2);
dscale = grib_power(-decimal_scale_factor, 10);
bscale = codes_power<T>(binary_scale_factor, 2);
dscale = codes_power<T>(-decimal_scale_factor, 10);
buflen = grib_byte_count(a);
buf = (unsigned char*)hand->buffer->data;

33
src/grib_accessor_class_data_complex_packing.cc
View File

@ -8,7 +8,8 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_api_internal_cpp.h"
#include "grib_ieeefloat.h"
#include "grib_scaling.h"
#include <math.h>
#include <algorithm>
/*
@ -446,7 +447,7 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
if (pen_j == sub_j) {
double* values;
d = grib_power(decimal_scale_factor, 10);
d = codes_power<double>(decimal_scale_factor, 10);
if (d) {
values = (double*)grib_context_malloc_clear(a->context, sizeof(double) * n_vals);
for (i = 0; i < n_vals; i++)
@ -546,10 +547,10 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
"%s: unable to find nearest_smaller_value of %g for %s", cclass_name, min, self->reference_value);
return GRIB_INTERNAL_ERROR;
}
d = grib_power(+decimal_scale_factor, 10);
d = codes_power<double>(+decimal_scale_factor, 10);
}
else {
d = grib_power(+decimal_scale_factor, 10);
d = codes_power<double>(+decimal_scale_factor, 10);
if (grib_get_nearest_smaller_value(gh, self->reference_value, d * min, &reference_value) != GRIB_SUCCESS) {
grib_context_log(gh->context, GRIB_LOG_ERROR,
"%s: unable to find nearest_smaller_value of %g for %s", cclass_name, d * min, self->reference_value);
@ -569,7 +570,7 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
}
}
}
s = grib_power(-binary_scale_factor, 2);
s = codes_power<double>(-binary_scale_factor, 2);
i = 0;
@ -698,9 +699,9 @@ static int unpack(grib_accessor* a, T* val, size_t* len)
T* scals = NULL;
T* pscals = NULL, *pval = NULL;
double s = 0;
double d = 0;
double laplacianOperator = 0;
T s = 0;
T d = 0;
T laplacianOperator = 0;
unsigned char* buf = NULL;
unsigned char* hres = NULL;
unsigned char* lres = NULL;
@ -713,7 +714,7 @@ static int unpack(grib_accessor* a, T* val, size_t* len)
long offsetdata = 0;
long bits_per_value = 0;
double reference_value = 0;
T reference_value = 0;
long binary_scale_factor = 0;
long decimal_scale_factor = 0;
@ -727,6 +728,7 @@ static int unpack(grib_accessor* a, T* val, size_t* len)
T operat = 0;
int bytes;
int err = 0;
double tmp;
decode_float_proc decode_float = NULL;
@ -743,8 +745,9 @@ static int unpack(grib_accessor* a, T* val, size_t* len)
return ret;
if ((ret = grib_get_long_internal(gh, self->bits_per_value, &bits_per_value)) != GRIB_SUCCESS)
return ret;
if ((ret = grib_get_double_internal(gh, self->reference_value, &reference_value)) != GRIB_SUCCESS)
if ((ret = grib_get_double_internal(gh, self->reference_value, &tmp)) != GRIB_SUCCESS)
return ret;
reference_value = tmp;
if ((ret = grib_get_long_internal(gh, self->binary_scale_factor, &binary_scale_factor)) != GRIB_SUCCESS)
return ret;
@ -758,8 +761,10 @@ static int unpack(grib_accessor* a, T* val, size_t* len)
if ((ret = grib_get_long(gh, self->ieee_floats, &ieee_floats)) != GRIB_SUCCESS)
return ret;
if ((ret = grib_get_double_internal(gh, self->laplacianOperator, &laplacianOperator)) != GRIB_SUCCESS)
if ((ret = grib_get_double_internal(gh, self->laplacianOperator, &tmp)) != GRIB_SUCCESS)
return ret;
laplacianOperator = tmp;
if ((ret = grib_get_long_internal(gh, self->sub_j, &sub_j)) != GRIB_SUCCESS)
return ret;
if ((ret = grib_get_long_internal(gh, self->sub_k, &sub_k)) != GRIB_SUCCESS)
@ -807,7 +812,7 @@ static int unpack(grib_accessor* a, T* val, size_t* len)
if (pen_j == sub_j) {
n_vals = (pen_j + 1) * (pen_j + 2);
d = grib_power(-decimal_scale_factor, 10);
d = codes_power<T>(-decimal_scale_factor, 10);
grib_ieee_decode_array<T>(a->context, buf, n_vals, bytes, val);
if (d) {
@ -821,8 +826,8 @@ static int unpack(grib_accessor* a, T* val, size_t* len)
lpos = 8 * (packed_offset - offsetdata);
s = grib_power(binary_scale_factor, 2);
d = grib_power(-decimal_scale_factor, 10);
s = codes_power<T>(binary_scale_factor, 2);
d = codes_power<T>(-decimal_scale_factor, 10);
scals = (T*)grib_context_malloc(a->context, maxv * sizeof(T));
if (!scals) return GRIB_OUT_OF_MEMORY;

5
src/grib_accessor_class_data_g1second_order_constant_width_packing.cc
View File

@ -8,6 +8,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_scaling.h"
#include "grib_api_internal.h"
/*
@ -298,8 +299,8 @@ static int unpack_double(grib_accessor* a, double* values, size_t* len)
printf("XXXXXXX extrabits=%ld pos=%ld\n",extrabits,pos);
}*/
s = grib_power(binary_scale_factor, 2);
d = grib_power(-decimal_scale_factor, 10);
s = codes_power<double>(binary_scale_factor, 2);
d = codes_power<double>(-decimal_scale_factor, 10);
for (i = 0; i < numberOfSecondOrderPackedValues; i++) {
values[i] = (double)(((X[i] * s) + reference_value) * d);
}

42
src/grib_accessor_class_data_g1second_order_general_extended_packing.cc
View File

@ -8,6 +8,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_scaling.h"
#include "grib_api_internal.h"
#include "grib_optimize_decimal_factor.h"
@ -364,7 +365,6 @@ static int unpack(grib_accessor* a, double* dvalues, float* fvalues, size_t* len
double reference_value;
long binary_scale_factor;
long decimal_scale_factor;
double s, d;
long j, count = 0;
long *groupWidths = NULL, *groupLengths = NULL;
long orderOfSPD = 0;
@ -532,8 +532,8 @@ static int unpack(grib_accessor* a, double* dvalues, float* fvalues, size_t* len
self->dvalues = (double*)grib_context_malloc_clear(a->context, sizeof(double) * numberOfValues);
}
s = grib_power(binary_scale_factor, 2);
d = grib_power(-decimal_scale_factor, 10);
double s = codes_power<double>(binary_scale_factor, 2);
double d = codes_power<double>(-decimal_scale_factor, 10);
for (i = 0; i < numberOfValues; i++) {
dvalues[i] = (double)(((X[i] * s) + reference_value) * d);
self->dvalues[i] = dvalues[i];
@ -551,8 +551,8 @@ static int unpack(grib_accessor* a, double* dvalues, float* fvalues, size_t* len
self->fvalues = (float*)grib_context_malloc_clear(a->context, sizeof(float) * numberOfValues);
}
s = grib_power(binary_scale_factor, 2);
d = grib_power(-decimal_scale_factor, 10);
float s = codes_power<float>(binary_scale_factor, 2);
float d = codes_power<float>(-decimal_scale_factor, 10);
for (i = 0; i < numberOfValues; i++) {
fvalues[i] = (float)(((X[i] * s) + reference_value) * d);
self->fvalues[i] = fvalues[i];
@ -756,7 +756,7 @@ static int pack_double_old(grib_accessor* a, const double* val, size_t *len)
if((ret = grib_get_long_internal(handle,self->decimal_scale_factor, &decimal_scale_factor))
!= GRIB_SUCCESS)
return ret;
decimal = grib_power(decimal_scale_factor,10);
decimal = codes_power<double>(decimal_scale_factor,10);
max*=decimal;
min*=decimal;
@ -790,7 +790,7 @@ static int pack_double_old(grib_accessor* a, const double* val, size_t *len)
GRIB_SUCCESS)
return ret;
divisor = grib_power(-binary_scale_factor,2);
divisor = codes_power<double>(-binary_scale_factor,2);
X=(long*)grib_context_malloc_clear(a->context,sizeof(long)*numberOfValues);
for(i=0;i< numberOfValues;i++){
X[i] = (((val[i]*decimal)-reference_value)*divisor)+0.5;
@ -861,7 +861,7 @@ static int pack_double_old(grib_accessor* a, const double* val, size_t *len)
}
}
groupWidthA=number_of_bits(handle, maxA-minA);
range=(long)grib_power(groupWidthA,2)-1;
range=(long)codes_power<double>(groupWidthA,2)-1;
offsetC=count+groupLengthA;
if (offsetC==numberOfValues) {
@ -949,7 +949,7 @@ static int pack_double_old(grib_accessor* a, const double* val, size_t *len)
if (minA>X[count+i]) minA=X[count+i];
}
groupWidthA=number_of_bits(maxA-minA);
range=(long)grib_power(groupWidthA,2)-1;
range=(long)codes_power<double>(groupWidthA,2)-1;
groupLengths[numberOfGroups]=groupLengthA;
groupWidths[numberOfGroups]=groupWidthA;
firstOrderValues[numberOfGroups] = maxA-range > 0 ? maxA-range : 0;
@ -1013,7 +1013,7 @@ static int pack_double_old(grib_accessor* a, const double* val, size_t *len)
if (minA>X[count+i]) minA=X[count+i];
}
groupWidthA=number_of_bits(maxA-minA);
range=(long)grib_power(groupWidthA,2)-1;
range=(long)codes_power<double>(groupWidthA,2)-1;
groupLengths[numberOfGroups]=groupLengthA;
groupWidths[numberOfGroups]=groupWidthA;
firstOrderValues[numberOfGroups] = maxA-range > 0 ? maxA-range : 0;
@ -1055,7 +1055,7 @@ static int pack_double_old(grib_accessor* a, const double* val, size_t *len)
if (min>firstOrderValues[i]) min=firstOrderValues[i];
}
widthOfFirstOrderValues=number_of_bits(handle, max-min);
firstOrderValuesMax=(long)grib_power(widthOfFirstOrderValues,2)-1;
firstOrderValuesMax=(long)codes_power<double>(widthOfFirstOrderValues,2)-1;
if (numberOfGroups>2) {
/* loop through all the groups except the last in reverse order to
@ -1079,7 +1079,7 @@ static int pack_double_old(grib_accessor* a, const double* val, size_t *len)
if (min>X[offset+j]) min=X[offset+j];
}
groupWidth=number_of_bits(handle, max-min);
range=(long)grib_power(groupWidth,2)-1;
range=(long)codes_power<double>(groupWidth,2)-1;
/* width of first order values has to be unchanged.*/
for (j=groupWidth;j<groupWidths[i];j++) {
@ -1392,8 +1392,8 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
&decimal_scale_factor, &binary_scale_factor, &reference_value)) != GRIB_SUCCESS)
return ret;
decimal = grib_power(decimal_scale_factor, 10);
divisor = grib_power(-binary_scale_factor, 2);
decimal = codes_power<double>(decimal_scale_factor, 10);
divisor = codes_power<double>(-binary_scale_factor, 2);
/*min = min * decimal;*/
/*max = max * decimal;*/
@ -1410,7 +1410,7 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
if ((ret = grib_get_long_internal(handle, self->decimal_scale_factor, &decimal_scale_factor)) != GRIB_SUCCESS)
return ret;
decimal = grib_power(decimal_scale_factor, 10);
decimal = codes_power<double>(decimal_scale_factor, 10);
min = min * decimal;
max = max * decimal;
@ -1421,7 +1421,7 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
}
binary_scale_factor = grib_get_binary_scale_fact(max, reference_value, bits_per_value, &ret);
divisor = grib_power(-binary_scale_factor, 2);
divisor = codes_power<double>(-binary_scale_factor, 2);
}
if ((ret = grib_set_long_internal(handle, self->binary_scale_factor, binary_scale_factor)) !=
@ -1515,7 +1515,7 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
}
}
groupWidthA = number_of_bits(handle, maxA - minA);
range = (long)grib_power(groupWidthA, 2) - 1;
range = (long)codes_power<double>(groupWidthA, 2) - 1;
offsetC = count + groupLengthA;
if (offsetC == numberOfValues) {
@ -1608,7 +1608,7 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
if (minA>X[count+i]) minA=X[count+i];
}
groupWidthA=number_of_bits(maxA-minA);
range=(long)grib_power(groupWidthA,2)-1;
range=(long)codes_power<double>(groupWidthA,2)-1;
groupLengths[numberOfGroups]=groupLengthA;
groupWidths[numberOfGroups]=groupWidthA;
firstOrderValues[numberOfGroups] = maxA-range > 0 ? maxA-range : 0;
@ -1673,7 +1673,7 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
if (minA>X[count+i]) minA=X[count+i];
}
groupWidthA=number_of_bits(maxA-minA);
range=(long)grib_power(groupWidthA,2)-1;
range=(long)codes_power<double>(groupWidthA,2)-1;
groupLengths[numberOfGroups]=groupLengthA;
groupWidths[numberOfGroups]=groupWidthA;
firstOrderValues[numberOfGroups] = maxA-range > 0 ? maxA-range : 0;
@ -1717,7 +1717,7 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
min = firstOrderValues[i];
}
widthOfFirstOrderValues = number_of_bits(handle, max - min);
firstOrderValuesMax = (long)grib_power(widthOfFirstOrderValues, 2) - 1;
firstOrderValuesMax = (long)codes_power<double>(widthOfFirstOrderValues, 2) - 1;
if (numberOfGroups > 2) {
/* loop through all the groups except the last in reverse order to
@ -1745,7 +1745,7 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
min = X[offset + j];
}
groupWidth = number_of_bits(handle, max - min);
range = (long)grib_power(groupWidth, 2) - 1;
range = (long)codes_power<double>(groupWidth, 2) - 1;
/* width of first order values has to be unchanged.*/
for (j = groupWidth; j < groupWidths[i]; j++) {

5
src/grib_accessor_class_data_g1second_order_general_packing.cc
View File

@ -8,6 +8,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_scaling.h"
#include "grib_api_internal.h"
#include <type_traits>
@ -287,8 +288,8 @@ static int unpack(grib_accessor* a, T* values, size_t* len)
}
}
s = grib_power(binary_scale_factor, 2);
d = grib_power(-decimal_scale_factor, 10);
s = codes_power<T>(binary_scale_factor, 2);
d = codes_power<T>(-decimal_scale_factor, 10);
for (i = 0; i < numberOfSecondOrderPackedValues; i++) {
values[i] = (T)(((X[i] * s) + reference_value) * d);
}

5
src/grib_accessor_class_data_g1second_order_row_by_row_packing.cc
View File

@ -8,6 +8,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_scaling.h"
#include "grib_api_internal.h"
/*
@ -411,8 +412,8 @@ static int unpack(grib_accessor* a, T* values, size_t* len)
k++;
}
s = grib_power(binary_scale_factor, 2);
d = grib_power(-decimal_scale_factor, 10);
s = codes_power<T>(binary_scale_factor, 2);
d = codes_power<T>(-decimal_scale_factor, 10);
for (i = 0; i < n; i++) {
values[i] = (T)(((X[i] * s) + reference_value) * d);
}

5
src/grib_accessor_class_data_g1simple_packing.cc
View File

@ -8,6 +8,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_scaling.h"
#include "grib_api_internal.h"
/*
@ -290,8 +291,8 @@ static int pack_double(grib_accessor* a, const double* cval, size_t* len)
if ((ret = grib_get_long_internal(grib_handle_of_accessor(a), self->offsetsection, &offsetsection)) != GRIB_SUCCESS)
return ret;
decimal = grib_power(decimal_scale_factor, 10);
divisor = grib_power(-binary_scale_factor, 2);
decimal = codes_power<double>(decimal_scale_factor, 10);
divisor = codes_power<double>(-binary_scale_factor, 2);
buflen = (((bits_per_value * n_vals) + 7) / 8) * sizeof(unsigned char);
if ((buflen + (offsetdata - offsetsection)) % 2) {

5
src/grib_accessor_class_data_g22order_packing.cc
View File

@ -9,6 +9,7 @@
*/
#include "grib_scaling.h"
#include "grib_api_internal.h"
#include <type_traits>
@ -749,8 +750,8 @@ static int unpack(grib_accessor* a, T* val, const size_t* len)
// de_spatial_difference (a->context, sec_val, n_vals, orderOfSpatialDifferencing, bias);
}
binary_s = (T)grib_power(binary_scale_factor, 2);
decimal_s = (T)grib_power(-decimal_scale_factor, 10);
binary_s = (T)codes_power<T>(binary_scale_factor, 2);
decimal_s = (T)codes_power<T>(-decimal_scale_factor, 10);
for (i = 0; i < n_vals; i++) {
if (sec_val[i] == LONG_MAX) {

9
src/grib_accessor_class_data_g2bifourier_packing.cc
View File

@ -12,6 +12,7 @@
* philippe.marguinaud@meteo.fr
*******************************/
#include "grib_scaling.h"
#include "grib_api_internal.h"
#include "grib_optimize_decimal_factor.h"
#include <math.h>
@ -665,8 +666,8 @@ static int unpack_double(grib_accessor* a, double* val, size_t* len)
buf = (unsigned char*)gh->buffer->data;
buf += grib_byte_offset(a);
s = grib_power(bt->binary_scale_factor, 2);
d = grib_power(-bt->decimal_scale_factor, 10);
s = codes_power<double>(bt->binary_scale_factor, 2);
d = codes_power<double>(-bt->decimal_scale_factor, 10);
/*
* Decode data
@ -815,8 +816,8 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
if (ret != GRIB_SUCCESS)
goto cleanup;
s = grib_power(-bt->binary_scale_factor, 2);
d = grib_power(+bt->decimal_scale_factor, 10);
s = codes_power<double>(-bt->binary_scale_factor, 2);
d = codes_power<double>(+bt->decimal_scale_factor, 10);
}
else {
bt->decimal_scale_factor = 0;

5
src/grib_accessor_class_data_g2simple_packing.cc
View File

@ -8,6 +8,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_scaling.h"
#include "grib_api_internal.h"
/*
@ -244,8 +245,8 @@ static int pack_double(grib_accessor* a, const double* cval, size_t* len)
if ((ret = grib_get_long_internal(grib_handle_of_accessor(a), self->decimal_scale_factor, &decimal_scale_factor)) != GRIB_SUCCESS)
return ret;
decimal = grib_power(decimal_scale_factor, 10);
divisor = grib_power(-binary_scale_factor, 2);
decimal = codes_power<double>(decimal_scale_factor, 10);
divisor = codes_power<double>(-binary_scale_factor, 2);
buflen = (((bits_per_value * n_vals) + 7) / 8) * sizeof(unsigned char);
buf = (unsigned char*)grib_context_buffer_malloc_clear(a->context, buflen);

9
src/grib_accessor_class_data_jpeg2000_packing.cc
View File

@ -8,6 +8,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_scaling.h"
#include "grib_api_internal.h"
/*
@ -271,8 +272,8 @@ static int unpack_double(grib_accessor* a, double* val, size_t* len)
self->dirty = 0;
bscale = grib_power(binary_scale_factor, 2);
dscale = grib_power(-decimal_scale_factor, 10);
bscale = codes_power<double>(binary_scale_factor, 2);
dscale = codes_power<double>(-decimal_scale_factor, 10);
/* TODO: This should be called upstream */
if (*len < n_vals)
@ -411,8 +412,8 @@ static int pack_double(grib_accessor* a, const double* cval, size_t* len)
if ((ret = grib_get_long_internal(grib_handle_of_accessor(a), self->decimal_scale_factor, &decimal_scale_factor)) != GRIB_SUCCESS)
return ret;
decimal = grib_power(decimal_scale_factor, 10);
divisor = grib_power(-binary_scale_factor, 2);
decimal = codes_power<double>(decimal_scale_factor, 10);
divisor = codes_power<double>(-binary_scale_factor, 2);
simple_packing_size = (((bits_per_value * n_vals) + 7) / 8) * sizeof(unsigned char);
buf = (unsigned char*)grib_context_malloc_clear(a->context, simple_packing_size + EXTRA_BUFFER_SIZE);

9
src/grib_accessor_class_data_png_packing.cc
View File

@ -8,6 +8,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_scaling.h"
#include "grib_api_internal.h"
#define PNG_ANYBITS
@ -248,8 +249,8 @@ static int unpack_double(grib_accessor* a, double* val, size_t* len)
if ((err = grib_get_long_internal(grib_handle_of_accessor(a), self->decimal_scale_factor, &decimal_scale_factor)) != GRIB_SUCCESS)
return err;
bscale = grib_power(binary_scale_factor, 2);
dscale = grib_power(-decimal_scale_factor, 10);
bscale = codes_power<double>(binary_scale_factor, 2);
dscale = codes_power<double>(-decimal_scale_factor, 10);
/* TODO: This should be called upstream */
if (*len < n_vals)
@ -502,7 +503,7 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
return GRIB_SUCCESS;
}
d = grib_power(decimal_scale_factor, 10);
d = codes_power<double>(decimal_scale_factor, 10);
max = val[0];
min = max;
@ -526,7 +527,7 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
}
binary_scale_factor = grib_get_binary_scale_fact(max, reference_value, bits_per_value, &err);
divisor = grib_power(-binary_scale_factor, 2);
divisor = codes_power<double>(-binary_scale_factor, 2);
#ifndef PNG_ANYBITS
Assert(bits_per_value % 8 == 0);

2
src/grib_accessor_class_data_raw_packing.cc
View File

@ -11,7 +11,7 @@
* Enrico Fucile
****************************/
#include "grib_api_internal_cpp.h"
#include "grib_ieeefloat.h"
#define PRE_PROCESSING_NONE 0
#define PRE_PROCESSING_DIFFERENCE 1

3
src/grib_accessor_class_data_run_length_packing.cc
View File

@ -8,6 +8,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_scaling.h"
#include "grib_api_internal.h"
/*
@ -203,7 +204,7 @@ static int unpack_double(grib_accessor* a, double* val, size_t* len)
if (decimal_scale_factor > 127) {
decimal_scale_factor = -(decimal_scale_factor - 128);
}
level_scale_factor = grib_power(-decimal_scale_factor, 10.0);
level_scale_factor = codes_power<double>(-decimal_scale_factor, 10.0);
levels = (double*)grib_context_malloc_clear(a->context, sizeof(double) * (number_of_level_values + 1));
levels[0] = missingValue;
for (i = 0; i < number_of_level_values; i++) {

5
src/grib_accessor_class_data_sh_packed.cc
View File

@ -8,6 +8,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_scaling.h"
#include "grib_api_internal.h"
#include <math.h>
/*
@ -328,8 +329,8 @@ static int unpack_double(grib_accessor* a, double* val, size_t* len)
lpos = 8 * (packed_offset - offsetdata);
s = grib_power(binary_scale_factor, 2);
d = grib_power(-decimal_scale_factor, 10);
s = codes_power<double>(binary_scale_factor, 2);
d = codes_power<double>(-decimal_scale_factor, 10);
scals = (double*)grib_context_malloc(a->context, maxv * sizeof(double));
Assert(scals);

5
src/grib_accessor_class_data_sh_unpacked.cc
View File

@ -8,6 +8,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_scaling.h"
#include "grib_api_internal.h"
#include <math.h>
/*
@ -307,8 +308,8 @@ static int unpack_double(grib_accessor* a, double* val, size_t* len)
lpos = 8 * (packed_offset - offsetdata);
s = grib_power(binary_scale_factor, 2);
d = grib_power(-decimal_scale_factor, 10);
s = codes_power<double>(binary_scale_factor, 2);
d = codes_power<double>(-decimal_scale_factor, 10);
scals = (double*)grib_context_malloc(a->context, maxv * sizeof(double));
Assert(scals);

28
src/grib_accessor_class_data_simple_packing.cc
View File

@ -8,9 +8,11 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_api_internal_cpp.h"
#include "grib_scaling.h"
#include "grib_bits_any_endian_simple.h"
#include "grib_optimize_decimal_factor.h"
#include <float.h>
#include <type_traits>
/*
This is used by make_class.pl
@ -235,8 +237,8 @@ static int unpack_double_element(grib_accessor* a, size_t idx, double* val)
}
Assert(idx < n_vals);
s = grib_power(binary_scale_factor, 2);
d = grib_power(-decimal_scale_factor, 10);
s = codes_power<double>(binary_scale_factor, 2);
d = codes_power<double>(-decimal_scale_factor, 10);
grib_context_log(a->context, GRIB_LOG_DEBUG,
"%s: %s: creating %s, %ld values (idx=%zu)",
@ -368,8 +370,8 @@ static int unpack(grib_accessor* a, T* val, size_t* len)
return GRIB_SUCCESS;
}
s = grib_power(binary_scale_factor, 2);
d = grib_power(-decimal_scale_factor, 10);
s = codes_power<T>(binary_scale_factor, 2);
d = codes_power<T>(-decimal_scale_factor, 10);
grib_context_log(a->context, GRIB_LOG_DEBUG,
"%s %s: Creating %s, %zu values", cclass_name, __func__, a->name, n_vals);
@ -506,8 +508,8 @@ static int _unpack_double(grib_accessor* a, double* val, size_t* len, unsigned c
return GRIB_SUCCESS;
}
s = grib_power(binary_scale_factor, 2);
d = grib_power(-decimal_scale_factor, 10);
s = codes_power<double>(binary_scale_factor, 2);
d = codes_power<double>(-decimal_scale_factor, 10);
grib_context_log(a->context, GRIB_LOG_DEBUG,
"%s %s: Creating %s, %zu values", cclass_name, __func__, a->name, n_vals);
@ -743,7 +745,7 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
/* decimal_scale_factor is given, binary_scale_factor=0 and bits_per_value is computed */
binary_scale_factor = 0;
decimal_scale_factor = decimal_scale_factor_get;
decimal = grib_power(decimal_scale_factor, 10);
decimal = codes_power<double>(decimal_scale_factor, 10);
min *= decimal;
max *= decimal;
@ -789,14 +791,14 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
return err;
}
else {
/* printf("max=%g reference_value=%g grib_power(-last,2)=%g decimal_scale_factor=%ld bits_per_value=%ld\n",
max,reference_value,grib_power(-last,2),decimal_scale_factor,bits_per_value);*/
/* printf("max=%g reference_value=%g codes_power<double>(-last,2)=%g decimal_scale_factor=%ld bits_per_value=%ld\n",
max,reference_value,codes_power<double>(-last,2),decimal_scale_factor,bits_per_value);*/
range = (max - min);
unscaled_min = min;
unscaled_max = max;
f = (grib_power(bits_per_value, 2) - 1);
minrange = grib_power(-last, 2) * f;
maxrange = grib_power(last, 2) * f;
f = (codes_power<double>(bits_per_value, 2) - 1);
minrange = codes_power<double>(-last, 2) * f;
maxrange = codes_power<double>(last, 2) * f;
while (range < minrange) {
decimal_scale_factor += 1;

3
src/grib_accessor_class_expanded_descriptors.cc
View File

@ -15,6 +15,7 @@
can appear
*/
#include "grib_scaling.h"
#include "grib_api_internal.h"
/*
This is used by make_class.pl
@ -483,7 +484,7 @@ static void __expand(grib_accessor* a, bufr_descriptors_array* unexpanded, bufr_
case 7:
if (us->Y) {
ccp->extraScale = us->Y;
ccp->referenceFactor = grib_power(us->Y, 10);
ccp->referenceFactor = codes_power<double>(us->Y, 10);
ccp->extraWidth = ((10 * us->Y) + 2) / 3;
}
else {

2
src/grib_accessor_class_gen.cc
View File

@ -14,9 +14,9 @@
* Shahram Najm *
***************************************************************************/
#include "grib_api_internal.h"
#include "grib_value.h"
#include <limits>
#include <cassert>
#include "grib_api_internal_cpp.h"
/*
This is used by make_class.pl

5
src/grib_accessor_class_second_order_bits_per_value.cc
View File

@ -12,6 +12,7 @@
* Enrico Fucile
**********************************/
#include "grib_scaling.h"
#include "grib_api_internal.h"
/*
This is used by make_class.pl
@ -225,8 +226,8 @@ static int unpack_long(grib_accessor* a, long* val, size_t* len)
min = values[i];
}
d = grib_power(decimalScaleFactor, 10);
b = grib_power(-binaryScaleFactor, 2);
d = codes_power<double>(decimalScaleFactor, 10);
b = codes_power<double>(-binaryScaleFactor, 2);
/* self->bitsPerValue=(long)ceil(log((double)((max-min)*d+1))/log(2.0))-binaryScaleFactor; */
/* See GRIB-540 for why we use ceil */

3
src/grib_accessor_class_simple_packing_error.cc
View File

@ -13,6 +13,7 @@
**************************************/
#include "grib_scaling.h"
#include "grib_api_internal.h"
/*
This is used by make_class.pl
@ -166,7 +167,7 @@ static int unpack_double(grib_accessor* a, double* val, size_t* len)
Assert(1 == 0);
if (bitsPerValue != 0)
*val = (*val + grib_power(binaryScaleFactor, 2)) * grib_power(-decimalScaleFactor, 10) * 0.5;
*val = (*val + codes_power<double>(binaryScaleFactor, 2)) * codes_power<double>(-decimalScaleFactor, 10) * 0.5;
if (ret == GRIB_SUCCESS)
*len = 1;

2
src/grib_api_internal.h
View File

@ -1563,8 +1563,6 @@ typedef struct j2k_encode_helper
#endif
/* This part is automatically generated by ./errors.pl, do not edit */
#ifndef grib_errors_internal_H
#define grib_errors_internal_H

16
src/grib_api_internal_cpp.h
View File

@ -1,16 +0,0 @@
/*
* (C) Copyright 2005- ECMWF.
*
* This software is licensed under the terms of the Apache Licence Version 2.0
* which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
*
* In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#pragma once
#include "grib_accessor.h"
#include "grib_value.h"
#include "grib_bits_any_endian_simple.h"
#include "grib_ieeefloat.h"

6
src/grib_bits_any_endian.cc
View File

@ -17,6 +17,8 @@
#include <stdint.h>
#endif
#include "grib_scaling.h"
#if GRIB_PTHREADS
static pthread_once_t once = PTHREAD_ONCE_INIT;
static pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
@ -392,7 +394,7 @@ int grib_encode_unsigned_longb(unsigned char* p, unsigned long val, long* bitp,
}
#ifdef DEBUG
{
unsigned long maxV = grib_power(nb, 2);
unsigned long maxV = codes_power<double>(nb, 2);
if (val > maxV) {
fprintf(stderr, "grib_encode_unsigned_longb: Value=%lu, but number of bits=%ld!\n", val, nb);
Assert(0);
@ -421,7 +423,7 @@ int grib_encode_size_tb(unsigned char* p, size_t val, long* bitp, long nb)
}
#ifdef DEBUG
{
size_t maxV = grib_power(nb, 2);
size_t maxV = codes_power<double>(nb, 2);
if (val > maxV) {
fprintf(stderr, "grib_encode_size_tb: Value=%lu, but number of bits=%ld!\n", val, nb);
Assert(0);

3
src/grib_bufr_descriptor.cc
View File

@ -9,6 +9,7 @@
*/
#include "grib_scaling.h"
#include "grib_api_internal.h"
bufr_descriptor* grib_bufr_descriptor_new(grib_accessor* tables_accessor, int code, int silent, int* err)
@ -109,7 +110,7 @@ void grib_bufr_descriptor_set_scale(bufr_descriptor* v, long scale)
v->scale = scale;
if (scale != 0)
v->type = BUFR_DESCRIPTOR_TYPE_DOUBLE;
v->factor = grib_power(-scale, 10);
v->factor = codes_power<double>(-scale, 10);
}
int grib_bufr_descriptor_can_be_missing(bufr_descriptor* v)

2
src/grib_header_compute.cc
View File

@ -48,7 +48,7 @@ static long op_neg(long a) {return -a;}
static double op_neg_d(double a) {return -a;}
static long op_pow(long a, long b) {return grib_power(a,b);}
static long op_pow(long a, long b) {return codes_power<double>(a,b);}
static long op_add(long a, long b) {return a+b;}
static long op_sub(long a, long b) {return a-b;}
static long op_div(long a, long b) {return a/b;}

21
src/grib_optimize_decimal_factor.cc
View File

@ -8,6 +8,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_scaling.h"
#include "grib_api_internal.h"
#include "grib_optimize_decimal_factor.h"
#include <math.h>
@ -50,9 +51,9 @@ static void factec(int* krep, const double pa, const int knbit, const long kdec,
}
/* Binary scale factor associated to kdec */
*ke = floor(log2((pa * grib_power(kdec, 10)) / (grib_power(knbit, 2) - 0.5))) + 1;
*ke = floor(log2((pa * codes_power<double>(kdec, 10)) / (codes_power<double>(knbit, 2) - 0.5))) + 1;
/* Encoded value for pa = max - min */
*knutil = floor(0.5 + pa * grib_power(kdec, 10) * grib_power(-*ke, 2));
*knutil = floor(0.5 + pa * codes_power<double>(kdec, 10) * codes_power<double>(-*ke, 2));
end:
return;
@ -98,14 +99,14 @@ int grib_optimize_decimal_factor(grib_accessor* a, const char* reference_value,
xtinyr4 = FLT_MIN;
xhuger4 = FLT_MAX;
inbint = grib_power(knbit, 2) - 1;
inbint = codes_power<double>(knbit, 2) - 1;
xnbint = (double)inbint;
/* Test decimal scale factors; keep the most suitable */
for (jdec = idecmin; jdec <= idecmax; jdec++) {
/* Fix a problem in GRIBEX */
if (compat_gribex)
if (pa * grib_power(jdec, 10) <= 1.E-12)
if (pa * codes_power<double>(jdec, 10) <= 1.E-12)
continue;
/* Check it will be possible to decode reference value with 32bit floats */
@ -125,7 +126,7 @@ int grib_optimize_decimal_factor(grib_accessor* a, const char* reference_value,
/* Check it will be possible to decode the maximum value of the fields using 32bit floats */
if (compat_32bit)
if (pmin * grib_power(jdec, 10) + xnbint * grib_power(ie, 2) >= xhuger4)
if (pmin * codes_power<double>(jdec, 10) + xnbint * codes_power<double>(ie, 2) >= xhuger4)
continue;
/* GRIB1 demands that the binary scale factor be encoded in a single byte */
@ -141,8 +142,8 @@ int grib_optimize_decimal_factor(grib_accessor* a, const char* reference_value,
}
if (inumax > 0) {
double decimal = grib_power(+*kdec, 10);
double divisor = grib_power(-*kbin, 2);
double decimal = codes_power<double>(+*kdec, 10);
double divisor = codes_power<double>(-*kbin, 2);
double min = pmin * decimal;
long vmin, vmax;
if (grib_get_nearest_smaller_value(gh, reference_value, min, ref) != GRIB_SUCCESS) {
@ -164,9 +165,9 @@ int grib_optimize_decimal_factor(grib_accessor* a, const char* reference_value,
int last = compat_gribex ? 99 : 127;
double min = pmin, max = pmax;
double range = max - min;
double f = grib_power(knbit, 2) - 1;
double minrange = grib_power(-last, 2) * f;
double maxrange = grib_power(+last, 2) * f;
double f = codes_power<double>(knbit, 2) - 1;
double minrange = codes_power<double>(-last, 2) * f;
double maxrange = codes_power<double>(+last, 2) * f;
double decimal = 1;
int err;

32
src/grib_scaling.cc
View File

@ -12,25 +12,13 @@
* Enrico Fucile
**************************************/
#include "grib_scaling.h"
#include "grib_api_internal.h"
/* Return n to the power of s */
double grib_power(long s, long n)
{
double divisor = 1.0;
if (s == 0)
return 1.0;
if (s == 1)
return n;
while (s < 0) {
divisor /= n;
s++;
}
while (s > 0) {
divisor *= n;
s--;
}
return divisor;
// Unfortunately, metkit uses grib_power() (illegal usage of private API)
// As soon as it is fixed, the wrapper below can be deleted
double grib_power(long s, long n) {
return codes_power<double>(s, n);
}
long grib_get_binary_scale_fact(double max, double min, long bpval, int* ret)
@ -43,14 +31,14 @@ long grib_get_binary_scale_fact(double max, double min, long bpval, int* ret)
const size_t ulong_size = sizeof(maxint) * 8;
/* See ECC-246
unsigned long maxint = grib_power(bpval,2) - 1;
unsigned long maxint = codes_power<double>(bpval,2) - 1;
double dmaxint=(double)maxint;
*/
if (bpval >= ulong_size) {
*ret = GRIB_OUT_OF_RANGE; /*overflow*/
return 0;
}
const double dmaxint = grib_power(bpval, 2) - 1;
const double dmaxint = codes_power<double>(bpval, 2) - 1;
maxint = (unsigned long)dmaxint; /* Now it's safe to cast */
*ret = 0;
@ -101,7 +89,7 @@ long grib_get_bits_per_value(double max, double min, long binary_scale_factor)
long scale = 0;
const long last = 127; /* Depends on edition, should be parameter */
unsigned long maxint = grib_power(binary_scale_factor, 2) - 1;
unsigned long maxint = codes_power<double>(binary_scale_factor, 2) - 1;
double dmaxint = (double)maxint;
if (maxint == 0)
@ -145,10 +133,10 @@ long grib_get_decimal_scale_fact(double max, double min, long bpval, long binary
long scale = 0;
const long last = 127; /* Depends on edition, should be parameter */
unsigned long maxint = grib_power(bpval, 2) - 1;
unsigned long maxint = codes_power<double>(bpval, 2) - 1;
double dmaxint = (double)maxint;
range *= grib_power(-binary_scale, 2);
range *= codes_power<double>(-binary_scale, 2);
Assert(bpval >= 1);
if (range == 0)

23
src/grib_scaling.h Normal file
View File

@ -0,0 +1,23 @@
#pragma once
template <typename T> T codes_power(long s, long n);
/* Return n to the power of s */
template <typename T>
T codes_power(long s, long n)
{
T divisor = 1.0;
if (s == 0)
return 1.0;
if (s == 1)
return n;
while (s < 0) {
divisor /= n;
s++;
}
while (s > 0) {
divisor *= n;
s--;
}
return divisor;
}