mirror of https://github.com/ecmwf/eccodes.git
Merge pull request #166 from joobog/feature/compile_time_ieee_ibm_tables
Compile time initialisation of IBM and IEEE tables
This commit is contained in:
shahramn
GitHub
commit
9250c83cb6
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
||||
|
|
|
|||
|
|
@ -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_api_internal.h"
|
||||
#include "grib_scaling.h"
|
||||
|
||||
#include <array>
|
||||
#include <type_traits>
|
||||
#include <cstdint>
|
||||
|
||||
|
||||
/**
|
||||
.. _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 uint8_t TABLESIZE = 128;
|
||||
static constexpr uint32_t mantissa_min = 0x100000;
|
||||
static constexpr uint32_t mantissa_max = 0xffffff;
|
||||
|
||||
public:
|
||||
static constexpr std::array<ValueType, TABLESIZE> e = []() {
|
||||
std::array<ValueType, TABLESIZE> multiplier{};
|
||||
for (uint8_t i = 1; i < TABLESIZE; ++i) {
|
||||
multiplier[i] = codes_power<ValueType>(i - 70, 16);
|
||||
}
|
||||
return multiplier;
|
||||
}();
|
||||
static constexpr std::array<ValueType, TABLESIZE> v = []() {
|
||||
std::array<ValueType, TABLESIZE> values{};
|
||||
for (uint8_t 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;
|
||||
};
|
||||
|
|
@ -13,116 +13,9 @@
|
|||
***************************************************************************/
|
||||
#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<double>();
|
||||
|
||||
static void binary_search(const double xx[], const unsigned long n, double x, unsigned long* j)
|
||||
{
|
||||
|
|
@ -142,18 +35,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 +44,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 +63,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,7 +100,7 @@ double grib_ieeefloat_error(double x)
|
|||
{
|
||||
unsigned long e = 0;
|
||||
|
||||
init_table_if_needed();
|
||||
//init_table_if_needed();
|
||||
|
||||
if (x < 0)
|
||||
x = -x;
|
||||
|
|
@ -237,7 +116,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 +135,7 @@ double grib_long_to_ieee(unsigned long x)
|
|||
Assert(0);
|
||||
}
|
||||
#endif
|
||||
init_table_if_needed();
|
||||
//init_table_if_needed();
|
||||
|
||||
if (c == 0 && m == 0)
|
||||
return 0;
|
||||
|
|
@ -288,7 +167,7 @@ unsigned long grib_ieee_nearest_smaller_to_long(double x)
|
|||
if (x == 0)
|
||||
return 0;
|
||||
|
||||
init_table_if_needed();
|
||||
//init_table_if_needed();
|
||||
|
||||
l = grib_ieee_to_long(x);
|
||||
y = grib_long_to_ieee(l);
|
||||
|
|
@ -334,7 +213,7 @@ int grib_nearest_smaller_ieee_float(double a, double* ret)
|
|||
{
|
||||
unsigned long l = 0;
|
||||
|
||||
init_table_if_needed();
|
||||
//init_table_if_needed();
|
||||
|
||||
if (a > ieee_table.vmax) {
|
||||
grib_context* c = grib_context_get_default();
|
||||
|
|
|
|||
|
|
@ -11,5 +11,66 @@
|
|||
#pragma once
|
||||
|
||||
#include "grib_api_internal.h"
|
||||
#include "grib_scaling.h"
|
||||
|
||||
#include <array>
|
||||
#include <type_traits>
|
||||
#include <cstdint>
|
||||
|
||||
template <typename T> 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) |
|
||||
+-------+----------------+----------------------+
|
||||
| 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 <typename ValueType>
|
||||
struct IeeeTable {
|
||||
private:
|
||||
static_assert(std::is_floating_point<ValueType>::value, "ValueType must be a floating point type");
|
||||
static constexpr uint8_t TABLESIZE = 255;
|
||||
static constexpr uint32_t mantissa_min = 0x800000;
|
||||
static constexpr uint32_t mantissa_max = 0xffffff;
|
||||
|
||||
public:
|
||||
static constexpr std::array<ValueType, TABLESIZE> e = []() {
|
||||
std::array<ValueType, TABLESIZE> multiplier{};
|
||||
for (uint8_t i = 1; i < TABLESIZE; ++i) {
|
||||
multiplier[i] = codes_power<ValueType>(i - 150, 2);
|
||||
}
|
||||
return multiplier;
|
||||
}();
|
||||
static constexpr std::array<ValueType, TABLESIZE> v = []() {
|
||||
std::array<ValueType, TABLESIZE> values{};
|
||||
for (uint8_t 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;
|
||||
};
|
||||
|
|
|
|||
|
|
@ -1,10 +1,8 @@
|
|||
#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)
|
||||
constexpr T codes_power(long s, long n)
|
||||
{
|
||||
T divisor = 1.0;
|
||||
if (s == 0)
|
||||
|
|
|
|||
Loading…
Reference in New Issue