diff --git a/src/eccodes_prototypes.h b/src/eccodes_prototypes.h index 278fdcf7e..3b405a05a 100644 --- a/src/eccodes_prototypes.h +++ b/src/eccodes_prototypes.h @@ -765,8 +765,6 @@ unsigned long grib_ibm_nearest_smaller_to_long(double x); int grib_nearest_smaller_ibm_float(double a, double* ret); /* grib_ieeefloat.cc*/ -double grib_ieee_table_e(unsigned long e); -double grib_ieee_table_v(unsigned long e); unsigned long grib_ieee_to_long(double x); double grib_ieeefloat_error(double x); double grib_long_to_ieee(unsigned long x); diff --git a/src/grib_ibmfloat.cc b/src/grib_ibmfloat.cc index a4b54b358..47f0eda35 100644 --- a/src/grib_ibmfloat.cc +++ b/src/grib_ibmfloat.cc @@ -9,91 +9,9 @@ */ #include "grib_api_internal.h" +#include "grib_ibmfloat.h" -#if GRIB_PTHREADS -static pthread_once_t once = PTHREAD_ONCE_INIT; -static pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER; - -static void init() -{ - pthread_mutexattr_t attr; - pthread_mutexattr_init(&attr); - pthread_mutexattr_settype(&attr, PTHREAD_MUTEX_RECURSIVE); - pthread_mutex_init(&mutex, &attr); - pthread_mutexattr_destroy(&attr); -} -#elif GRIB_OMP_THREADS -static int once = 0; -static omp_nest_lock_t mutex; - -static void init() -{ - GRIB_OMP_CRITICAL(lock_grib_ibmfloat_c) - { - if (once == 0) { - omp_init_nest_lock(&mutex); - once = 1; - } - } -} -#endif - -typedef struct ibm_table_t ibm_table_t; - -struct ibm_table_t -{ - int inited; - double e[128]; - double v[128]; - double vmin; - double vmax; -}; - -static ibm_table_t ibm_table = { 0, {0,}, {0,}, 0, 0 }; - - -/** -.. seealso:: - Documentation for ``init_ieee_float()`` in ``grib_ieeefloat.cc`` -*/ - -static void init_ibm_table() -{ - if (!ibm_table.inited) { - unsigned long i; - unsigned long mmin = 0x100000; - unsigned long mmax = 0xffffff; - double e = 1; - for (i = 1; i <= 57; i++) { - e *= 16; - ibm_table.e[i + 70] = e; - ibm_table.v[i + 70] = e * mmin; - } - ibm_table.e[70] = 1; - ibm_table.v[70] = mmin; - e = 1; - for (i = 1; i <= 70; i++) { - e /= 16; - ibm_table.e[70 - i] = e; - ibm_table.v[70 - i] = e * mmin; - } - ibm_table.vmin = ibm_table.v[0]; - ibm_table.vmax = ibm_table.e[127] * mmax; - ibm_table.inited = 1; - /*for (i=0;i<128;i++) printf("++++ ibm_table.v[%d]=%g\n",i,ibm_table.v[i]);*/ - } -} - -static void init_table_if_needed() -{ - GRIB_MUTEX_INIT_ONCE(&once, &init); - GRIB_MUTEX_LOCK(&mutex); - - if (!ibm_table.inited) - init_ibm_table(); - - GRIB_MUTEX_UNLOCK(&mutex); -} +constexpr auto ibm_table = IbmTable{}; static void binary_search(const double xx[], const unsigned long n, double x, unsigned long* j) { @@ -122,8 +40,6 @@ unsigned long grib_ibm_to_long(double x) unsigned long e = 0; double rmmax = mmax + 0.5; - init_table_if_needed(); - /* printf("\ngrib_ibm_to_long: x=%.20e\n",x); */ if (x < 0) { s = 1; @@ -143,7 +59,7 @@ unsigned long grib_ibm_to_long(double x) return 0; } - binary_search(ibm_table.v, 127, x, &e); + binary_search(ibm_table.v.data(), 127, x, &e); /* printf("grib_ibm_to_long: e=%ld\n",e); */ @@ -180,8 +96,6 @@ double grib_ibmfloat_error(double x) { unsigned long e = 0; - init_table_if_needed(); - if (x < 0) x = -x; @@ -196,7 +110,7 @@ double grib_ibmfloat_error(double x) return 0; } - binary_search(ibm_table.v, 127, x, &e); + binary_search(ibm_table.v.data(), 127, x, &e); return ibm_table.e[e]; } @@ -209,8 +123,6 @@ double grib_long_to_ibm(unsigned long x) double val = m; - init_table_if_needed(); - /*if(x == 0) return 0;*/ if (c == 0 && m <= 1) return 0; @@ -225,13 +137,11 @@ double grib_long_to_ibm(unsigned long x) double grib_ibm_table_e(unsigned long e) { - init_table_if_needed(); return ibm_table.e[e]; } double grib_ibm_table_v(unsigned long e) { - init_table_if_needed(); return ibm_table.v[e]; } @@ -247,8 +157,6 @@ unsigned long grib_ibm_nearest_smaller_to_long(double x) if (x == 0) return 0; - init_table_if_needed(); - l = grib_ibm_to_long(x); y = grib_long_to_ibm(l); @@ -293,8 +201,6 @@ int grib_nearest_smaller_ibm_float(double a, double* ret) { unsigned long l = 0; - init_table_if_needed(); - if (a > ibm_table.vmax) return GRIB_INTERNAL_ERROR; diff --git a/src/grib_ibmfloat.h b/src/grib_ibmfloat.h new file mode 100644 index 000000000..e521076f5 --- /dev/null +++ b/src/grib_ibmfloat.h @@ -0,0 +1,74 @@ +/* + * (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_scaling.h" + +#include +#include +#include + +/** +.. _init_ieee_table: + +Init IBM Floats Table +===================== + +Initializes the ibm_table with IBM Float values. Nearest smaller values (e.g., reference values for grid_simple) are taken from this table. + +Details +------- + +The table layout is as follows: + ++-------+----------------+------------------------+ +| idx (i) | multiplier (e) | value (v = mmin * e) | ++-------+----------------+------------------------+ +| 0 | 16^(-70) | 0x100000 * 2^(-70) | +| 1 | 16^(-69) | 0x100000 * 2^(-69) | +| ... | ... | ... | +| 126 | 16^56 | 0x100000 * 2^56 | +| 127 | 16^57 | 0x100000 * 2^57 | ++-------+----------------+------------------------+ + +The vmin and vmax boundaries are defined as: + +- vmin = 0x100000 * 2^(-70) +- vmax = 0xffffff * 2^57 +*/ + +struct IbmTable { +private: + using ValueType = double; + static constexpr int TABLESIZE = 128; + static constexpr uint32_t mantissa_min = 0x100000; + static constexpr uint32_t mantissa_max = 0xffffff; + +public: + static constexpr std::array e = []() { + std::array multiplier{}; + multiplier[0] = 0; + for (int i = 0; i < TABLESIZE; ++i) { + multiplier[i] = codes_power(i - 70, 16); + } + return multiplier; + }(); + static constexpr std::array v = []() { + std::array values{}; + values[0] = 0; + for (int i = 0; i < TABLESIZE; ++i) { + values[i] = e[i] * mantissa_min; + } + return values; + }(); + static constexpr ValueType vmin = e[0] * mantissa_min; + static constexpr ValueType vmax = e[127] * mantissa_max; +}; diff --git a/src/grib_ieeefloat.cc b/src/grib_ieeefloat.cc index e638ab4e9..5c8ffbb18 100644 --- a/src/grib_ieeefloat.cc +++ b/src/grib_ieeefloat.cc @@ -8,121 +8,11 @@ * virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. */ -/*************************************************************************** - * Enrico Fucile - 06.01.2009 * - ***************************************************************************/ #include "grib_ieeefloat.h" -#if GRIB_PTHREADS -static pthread_once_t once = PTHREAD_ONCE_INIT; -static pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER; - -static void init() -{ - pthread_mutexattr_t attr; - pthread_mutexattr_init(&attr); - pthread_mutexattr_settype(&attr, PTHREAD_MUTEX_RECURSIVE); - pthread_mutex_init(&mutex, &attr); - pthread_mutexattr_destroy(&attr); -} -#elif GRIB_OMP_THREADS -static int once = 0; -static omp_nest_lock_t mutex; - -static void init() -{ - GRIB_OMP_CRITICAL(lock_grib_ieeefloat_c) - { - if (once == 0) { - omp_init_nest_lock(&mutex); - once = 1; - } - } -} -#endif - /* See old implementation in src/deprecated/grib_ieeefloat.c */ -typedef struct ieee_table_t ieee_table_t; - -struct ieee_table_t -{ - int inited; - double e[255]; - double v[255]; - double vmin; - double vmax; -}; - -static ieee_table_t ieee_table = { 0, {0,}, {0,}, 0, 0 }; - - -/** -.. _init_ieee_table: - -Init IEEE Table -=============== - -Initializes the ieee_table with IEEE754 single precision (32-bit) values. Nearest smaller values (e.g., reference values for grid_simple and grid_ccsds) are taken from this table. - -Details -------- - -The table layout is as follows: - -+-------+----------------+----------------------+ -| idx (i) | multiplier (e) | value (v = mmin * e) | -+-------+----------------+----------------------+ -| 1 | 2^(-149) | 0x800000 * 2^(-149) | -| 2 | 2^(-148) | 0x800000 * 2^(-148) | -| ... | ... | ... | -| 253 | 2^103 | 0x800000 * 2^103 | -| 254 | 2^104 | 0x800000 * 2^104 | -+-------+----------------+----------------------+ - -The vmin and vmax boundaries are defined as: - -- vmin = 0x800000 * 2^(-149) -- vmax = 0xffffff * 2^104 -*/ - -static void init_ieee_table() -{ - if (!ieee_table.inited) { - unsigned long i; - unsigned long mmin = 0x800000; // minimum mantissa - unsigned long mmax = 0xffffff; // maximum mantissa - double e = 1; - for (i = 1; i <= 104; i++) { - e *= 2; - ieee_table.e[i + 150] = e; - ieee_table.v[i + 150] = e * mmin; - } - ieee_table.e[150] = 1; - ieee_table.v[150] = mmin; - e = 1; - for (i = 1; i < 150; i++) { - e /= 2; - ieee_table.e[150 - i] = e; - ieee_table.v[150 - i] = e * mmin; - } - ieee_table.vmin = ieee_table.v[1]; - ieee_table.vmax = ieee_table.e[254] * mmax; - ieee_table.inited = 1; - /*for (i=0;i<128;i++) printf("++++ ieee_table.v[%d]=%g\n",i,ieee_table.v[i]);*/ - } -} - -static void init_table_if_needed() -{ - GRIB_MUTEX_INIT_ONCE(&once, &init); - GRIB_MUTEX_LOCK(&mutex); - - if (!ieee_table.inited) - init_ieee_table(); - - GRIB_MUTEX_UNLOCK(&mutex); -} +constexpr auto ieee_table = IeeeTable(); static void binary_search(const double xx[], const unsigned long n, double x, unsigned long* j) { @@ -142,18 +32,6 @@ static void binary_search(const double xx[], const unsigned long n, double x, un *j = jl; } -double grib_ieee_table_e(unsigned long e) -{ - init_table_if_needed(); - return ieee_table.e[e]; -} - -double grib_ieee_table_v(unsigned long e) -{ - init_table_if_needed(); - return ieee_table.v[e]; -} - unsigned long grib_ieee_to_long(double x) { unsigned long s = 0; @@ -163,8 +41,6 @@ unsigned long grib_ieee_to_long(double x) unsigned long e = 0; double rmmax = mmax + 0.5; - init_table_if_needed(); - /* printf("\ngrib_ieee_to_long: x=%.20e\n",x); */ if (x < 0) { s = 1; @@ -184,7 +60,7 @@ unsigned long grib_ieee_to_long(double x) return 0; } - binary_search(ieee_table.v, 254, x, &e); + binary_search(ieee_table.v.data(), 254, x, &e); /* printf("grib_ieee_to_long: e=%ld\n",e); */ @@ -221,8 +97,6 @@ double grib_ieeefloat_error(double x) { unsigned long e = 0; - init_table_if_needed(); - if (x < 0) x = -x; @@ -237,7 +111,7 @@ double grib_ieeefloat_error(double x) return 0; } - binary_search(ieee_table.v, 254, x, &e); + binary_search(ieee_table.v.data(), 254, x, &e); return ieee_table.e[e]; } @@ -256,7 +130,6 @@ double grib_long_to_ieee(unsigned long x) Assert(0); } #endif - init_table_if_needed(); if (c == 0 && m == 0) return 0; @@ -288,8 +161,6 @@ unsigned long grib_ieee_nearest_smaller_to_long(double x) if (x == 0) return 0; - init_table_if_needed(); - l = grib_ieee_to_long(x); y = grib_long_to_ieee(l); @@ -334,10 +205,8 @@ int grib_nearest_smaller_ieee_float(double a, double* ret) { unsigned long l = 0; - init_table_if_needed(); - if (a > ieee_table.vmax) { - grib_context* c = grib_context_get_default(); + const grib_context* c = grib_context_get_default(); grib_context_log(c, GRIB_LOG_ERROR, "Number is too large: x=%e > xmax=%e (IEEE float)", a, ieee_table.vmax); return GRIB_INTERNAL_ERROR; diff --git a/src/grib_ieeefloat.h b/src/grib_ieeefloat.h index 74944f2ac..fb198a3cb 100644 --- a/src/grib_ieeefloat.h +++ b/src/grib_ieeefloat.h @@ -11,5 +11,69 @@ #pragma once #include "grib_api_internal.h" +#include "grib_scaling.h" + +#include +#include +#include template int grib_ieee_decode_array(grib_context* c, unsigned char* buf, size_t nvals, int bytes, T* val); + +/** +.. _init_ieee_table: + +Init IEEE Table +=============== + +Initializes the ieee_table with IEEE754 single precision (32-bit) values. Nearest smaller values (e.g., reference values for grid_simple and grid_ccsds) are taken from this table. + +Details +------- + +The table layout is as follows: + ++-------+----------------+------------------------+ +| idx (i) | multiplier (e) | value (v = mmin * e) | ++-------+----------------+------------------------+ +| 0 | 0 | 0 | +| 1 | 2^(-149) | 0x800000 * 2^(-149) | +| 2 | 2^(-148) | 0x800000 * 2^(-148) | +| ... | ... | ... | +| 253 | 2^103 | 0x800000 * 2^103 | +| 254 | 2^104 | 0x800000 * 2^104 | ++-------+----------------+------------------------+ + +The vmin and vmax boundaries are defined as: + +- vmin = 0x800000 * 2^(-149) +- vmax = 0xffffff * 2^104 +*/ + +template +struct IeeeTable { +private: + static_assert(std::is_floating_point::value, "ValueType must be a floating point type"); + static constexpr int TABLESIZE = 255; + static constexpr uint32_t mantissa_min = 0x800000; + static constexpr uint32_t mantissa_max = 0xffffff; + +public: + static constexpr std::array e = []() { + std::array multiplier{}; + multiplier[0] = 0; + for (int i = 1; i < TABLESIZE; ++i) { + multiplier[i] = codes_power(i - 150, 2); + } + return multiplier; + }(); + static constexpr std::array v = []() { + std::array values{}; + values[0] = 0; + for (int i = 1; i < TABLESIZE; ++i) { + values[i] = e[i] * mantissa_min; + } + return values; + }(); + static constexpr ValueType vmin = e[1] * mantissa_min; + static constexpr ValueType vmax = e[254] * mantissa_max; +}; diff --git a/src/grib_scaling.h b/src/grib_scaling.h index 53e3b1d4d..c48b3aa2c 100644 --- a/src/grib_scaling.h +++ b/src/grib_scaling.h @@ -1,10 +1,8 @@ #pragma once -template T codes_power(long s, long n); - /* Return n to the power of s */ template -T codes_power(long s, long n) +constexpr T codes_power(long s, long n) { T divisor = 1.0; if (s == 0)