1 /*
2 * (C) Copyright 2005- ECMWF.
3 *
4 * This software is licensed under the terms of the Apache Licence Version 2.0
5 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
6 *
7 * In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
8 * virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
9 */
10
11 /***************************************************************************
12 * Enrico Fucile - 06.01.2009 *
13 * *
14 ***************************************************************************/
15 #include "grib_api_internal.h"
16
17 #if GRIB_PTHREADS
18 static pthread_once_t once = PTHREAD_ONCE_INIT;
19 static pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
20
init()21 static void init()
22 {
23 pthread_mutexattr_t attr;
24 pthread_mutexattr_init(&attr);
25 pthread_mutexattr_settype(&attr, PTHREAD_MUTEX_RECURSIVE);
26 pthread_mutex_init(&mutex, &attr);
27 pthread_mutexattr_destroy(&attr);
28 }
29 #elif GRIB_OMP_THREADS
30 static int once = 0;
31 static omp_nest_lock_t mutex;
32
init()33 static void init()
34 {
35 GRIB_OMP_CRITICAL(lock_grib_ieeefloat_c)
36 {
37 if (once == 0) {
38 omp_init_nest_lock(&mutex);
39 once = 1;
40 }
41 }
42 }
43 #endif
44
45 #if 1
46
47 typedef struct ieee_table_t ieee_table_t;
48
49 struct ieee_table_t
50 {
51 int inited;
52 double e[255];
53 double v[255];
54 double vmin;
55 double vmax;
56 };
57
58 static ieee_table_t ieee_table = { 0, {0,}, {0,}, 0, 0 };
59
init_ieee_table()60 static void init_ieee_table()
61 {
62 if (!ieee_table.inited) {
63 unsigned long i;
64 unsigned long mmin = 0x800000;
65 unsigned long mmax = 0xffffff;
66 double e = 1;
67 for (i = 1; i <= 104; i++) {
68 e *= 2;
69 ieee_table.e[i + 150] = e;
70 ieee_table.v[i + 150] = e * mmin;
71 }
72 ieee_table.e[150] = 1;
73 ieee_table.v[150] = mmin;
74 e = 1;
75 for (i = 1; i < 150; i++) {
76 e /= 2;
77 ieee_table.e[150 - i] = e;
78 ieee_table.v[150 - i] = e * mmin;
79 }
80 ieee_table.vmin = ieee_table.v[1];
81 ieee_table.vmax = ieee_table.e[254] * mmax;
82 ieee_table.inited = 1;
83 /*for (i=0;i<128;i++) printf("++++ ieee_table.v[%d]=%g\n",i,ieee_table.v[i]);*/
84 }
85 }
86
init_table_if_needed()87 static void init_table_if_needed()
88 {
89 GRIB_MUTEX_INIT_ONCE(&once, &init)
90 GRIB_MUTEX_LOCK(&mutex)
91
92 if (!ieee_table.inited)
93 init_ieee_table();
94
95 GRIB_MUTEX_UNLOCK(&mutex)
96 }
97
binary_search(double xx[],const unsigned long n,double x,unsigned long * j)98 static void binary_search(double xx[], const unsigned long n, double x, unsigned long* j)
99 {
100 /*These routine works only on ascending ordered arrays*/
101 unsigned long ju, jm, jl;
102 jl = 0;
103 ju = n;
104 while (ju - jl > 1) {
105 jm = (ju + jl) >> 1;
106 /* printf("jl=%lu jm=%lu ju=%lu\n",jl,jm,ju); */
107 /* printf("xx[jl]=%.10e xx[jm]=%.10e xx[ju]=%.10e\n",xx[jl],xx[jm],xx[ju]); */
108 if (x >= xx[jm])
109 jl = jm;
110 else
111 ju = jm;
112 }
113 *j = jl;
114 }
115
grib_ieee_table_e(unsigned long e)116 double grib_ieee_table_e(unsigned long e)
117 {
118 init_table_if_needed();
119 return ieee_table.e[e];
120 }
121
grib_ieee_table_v(unsigned long e)122 double grib_ieee_table_v(unsigned long e)
123 {
124 init_table_if_needed();
125 return ieee_table.v[e];
126 }
127
grib_ieee_to_long(double x)128 unsigned long grib_ieee_to_long(double x)
129 {
130 unsigned long s = 0;
131 unsigned long mmax = 0xffffff;
132 unsigned long mmin = 0x800000;
133 unsigned long m = 0;
134 unsigned long e = 0;
135 double rmmax = mmax + 0.5;
136
137 init_table_if_needed();
138
139 /* printf("\ngrib_ieee_to_long: x=%.20e\n",x); */
140 if (x < 0) {
141 s = 1;
142 x = -x;
143 }
144
145 /* Underflow */
146 if (x < ieee_table.vmin) {
147 /*printf("grib_ieee_to_long: (x < ieee_table.vmin) x=%.20e vmin=%.20e v=0x%lX\n",x,ieee_table.vmin,(s<<31));*/
148 return (s << 31);
149 }
150
151 /* Overflow */
152 if (x > ieee_table.vmax) {
153 fprintf(stderr, "grib_ieee_to_long: Number is too large: x=%.20e > xmax=%.20e\n", x, ieee_table.vmax);
154 Assert(0);
155 return 0;
156 }
157
158 binary_search(ieee_table.v, 254, x, &e);
159
160 /* printf("grib_ieee_to_long: e=%ld\n",e); */
161
162 x /= ieee_table.e[e];
163
164 /* printf("grib_ieee_to_long: x=%.20e\n",x); */
165
166 while (x < mmin) {
167 x *= 2;
168 e--;
169 /* printf("grib_ieee_to_long (e--): x=%.20e e=%ld \n",x,e); */
170 }
171
172 while (x > rmmax) {
173 x /= 2;
174 e++;
175 /* printf("grib_ieee_to_long (e++): x=%.20e e=%ld \n",x,e); */
176 }
177
178 m = x + 0.5;
179 /* printf("grib_ieee_to_long: m=0x%lX (%lu) x=%.10e \n",m,m,x ); */
180 if (m > mmax) {
181 e++;
182 m = 0x800000;
183 /* printf("grib_ieee_to_long: ( m > mmax ) m=0x%lX (%lu) x=%.10e \n",m,m,x ); */
184 }
185
186 /* printf("grib_ieee_to_long: s=%lu c=%lu (0x%lX) m=%lu (0x%lX)\n",s,e,e,m,m ); */
187
188 return (s << 31) | (e << 23) | (m & 0x7fffff);
189 }
190
grib_ieeefloat_error(double x)191 double grib_ieeefloat_error(double x)
192 {
193 unsigned long e = 0;
194
195 init_table_if_needed();
196
197 if (x < 0)
198 x = -x;
199
200 /* Underflow */
201 if (x < ieee_table.vmin)
202 return ieee_table.vmin;
203
204 /* Overflow */
205 if (x > ieee_table.vmax) {
206 fprintf(stderr, "grib_ieeefloat_error: Number is too large: x=%.20e > xmax=%.20e\n", x, ieee_table.vmax);
207 Assert(0);
208 return 0;
209 }
210
211 binary_search(ieee_table.v, 254, x, &e);
212
213 return ieee_table.e[e];
214 }
215
grib_long_to_ieee(unsigned long x)216 double grib_long_to_ieee(unsigned long x)
217 {
218 unsigned long s = x & 0x80000000;
219 unsigned long c = (x & 0x7f800000) >> 23;
220 unsigned long m = (x & 0x007fffff);
221
222 double val;
223
224 #ifdef DEBUG
225 if (x > 0 && x < 0x800000) {
226 fprintf(stderr, "grib_long_to_ieee: Invalid input %lu\n", x);
227 Assert(0);
228 }
229 #endif
230 init_table_if_needed();
231
232 if (c == 0 && m == 0)
233 return 0;
234
235 if (c == 0) {
236 m |= 0x800000;
237 c = 1;
238 }
239 else
240 m |= 0x800000;
241
242 val = m * ieee_table.e[c];
243 if (s)
244 val = -val;
245
246 return val;
247 }
248
grib_ieee_nearest_smaller_to_long(double x)249 unsigned long grib_ieee_nearest_smaller_to_long(double x)
250 {
251 unsigned long l;
252 unsigned long e;
253 unsigned long m;
254 unsigned long s;
255 unsigned long mmin = 0x800000;
256 double y, eps;
257
258 if (x == 0)
259 return 0;
260
261 init_table_if_needed();
262
263 l = grib_ieee_to_long(x);
264 y = grib_long_to_ieee(l);
265
266 if (x < y) {
267 if (x < 0 && -x < ieee_table.vmin) {
268 l = 0x80800000;
269 }
270 else {
271 e = (l & 0x7f800000) >> 23;
272 m = (l & 0x007fffff) | 0x800000;
273 s = l & 0x80000000;
274
275 if (m == mmin) {
276 /* printf("grib_ieee_nearest_smaller_to_long: m == mmin (0x%lX) e=%lu\n",m,e); */
277 e = s ? e : e - 1;
278 if (e < 1)
279 e = 1;
280 if (e > 254)
281 e = 254;
282 /* printf("grib_ieee_nearest_smaller_to_long: e=%lu \n",e); */
283 }
284
285 eps = ieee_table.e[e];
286
287 /* printf("grib_ieee_nearest_smaller_to_long: x<y\n"); */
288 l = grib_ieee_to_long(y - eps);
289 /* printf("grib_ieee_nearest_smaller_to_long: grib_ieee_to_long(y-eps)=0x%lX y=%.10e eps=%.10e x=%.10e\n",l,y,eps,x); */
290 }
291 }
292 else
293 return l;
294
295 if (x < grib_long_to_ieee(l)) {
296 printf("grib_ieee_nearest_smaller_to_long: x=%.20e grib_long_to_ieee(0x%lX)=%.20e\n", x, l, grib_long_to_ieee(l));
297 Assert(x >= grib_long_to_ieee(l));
298 }
299
300 return l;
301 }
302
grib_nearest_smaller_ieee_float(double a,double * ret)303 int grib_nearest_smaller_ieee_float(double a, double* ret)
304 {
305 unsigned long l = 0;
306
307 init_table_if_needed();
308
309 if (a > ieee_table.vmax)
310 return GRIB_INTERNAL_ERROR;
311
312 l = grib_ieee_nearest_smaller_to_long(a);
313 *ret = grib_long_to_ieee(l);
314 return GRIB_SUCCESS;
315 }
316
317 #else
318 /* old code to be deleted */
319
grib_ieeefloat_error(double x)320 double grib_ieeefloat_error(double x)
321 {
322 return 0;
323 }
324
grib_long_to_ieee(unsigned long x)325 double grib_long_to_ieee(unsigned long x)
326 {
327 unsigned long s = x & 0x80000000;
328 unsigned long c = (x & 0x7f800000) >> 23;
329 unsigned long m;
330 double val;
331 long e;
332
333 if (x == 0)
334 return 0;
335 Assert(c != 255);
336
337 if (c == 0) {
338 m = x & 0x007fffff;
339 e = -126 - 23;
340 }
341 else {
342 m = (x & 0x007fffff) | (1 << 23);
343 e = c - 127 - 23;
344 }
345
346 val = m;
347
348 while (e < 0) {
349 val /= 2.0;
350 e++;
351 }
352 while (e > 0) {
353 val *= 2.0;
354 e--;
355 }
356
357 if (s)
358 val = -val;
359
360 return val;
361 }
362
grib_nearest_smaller_ieee_float(double a,double * x)363 int grib_nearest_smaller_ieee_float(double a, double* x)
364 {
365 double e = grib_long_to_ieee(grib_ieee_to_long(a));
366 double b = a;
367
368 /* printf("----> a=%g e=%g e-a=%g\n",a,e,e-a); */
369
370 if (e > b) {
371 unsigned long ub = grib_ieee_to_long(b);
372 unsigned long ue;
373 while (e > b) {
374 /* printf("a=%g e=%g e-a=%g\n",a,e,e-a); */
375 a -= (e - a);
376 e = grib_long_to_ieee(grib_ieee_to_long(a));
377 }
378 ue = grib_ieee_to_long(e);
379 Assert((ue - ub) == 1);
380 }
381
382 Assert(b >= e);
383 *x = e;
384 return GRIB_SUCCESS;
385 }
386
grib_ieee_to_long(double x)387 unsigned long grib_ieee_to_long(double x)
388 {
389 /* double y = x; */
390 unsigned long s = 0;
391 unsigned long m;
392 long p = 0;
393 unsigned long e = 0;
394
395 if (x == 0)
396 return 0;
397
398 if (x < 0) {
399 s = 1;
400 x = -x;
401 }
402 while (x < 2) {
403 x *= 2;
404 p--;
405 }
406
407 while (x >= 2) {
408 x /= 2;
409 p++;
410 }
411
412 if (p > 127) {
413 /* Overflow */
414 e = 255;
415 m = 0;
416 }
417 else if (p < -126) {
418 /* int i; */
419 e = 0;
420 /* printf("p=%ld x=%g %ld\n",p,x,p+126+23); */
421 m = x * grib_power(p + 126 + 23, 2);
422 }
423 else {
424 e = p + 127;
425 m = x * (1 << 23);
426 m &= 0x007fffff;
427 }
428
429 m = (s << 31) | (e << 23) | m;
430
431 return m;
432 }
433
434 #endif
435
436 #ifdef IEEE
437
438
439 /*
440 * To make these two routines consistent to grib_ieee_to_long and grib_long_to_ieee,
441 * we should not do any byte swapping but rather perform a raw copy.
442 * Byte swapping is actually implemented in grib_decode_unsigned_long and
443 * grib_encode_unsigned_long.
444 */
grib_ieee64_to_long(double x)445 unsigned long grib_ieee64_to_long(double x)
446 {
447 unsigned long lval;
448 DebugAssert(sizeof(double) == sizeof(long));
449 memcpy(&lval, &x, sizeof(long));
450 return lval;
451 }
452
grib_long_to_ieee64(unsigned long x)453 double grib_long_to_ieee64(unsigned long x)
454 {
455 double dval;
456 DebugAssert(sizeof(double) == sizeof(long));
457 memcpy(&dval, &x, sizeof(long));
458 return dval;
459 }
460
grib_ieee_decode_array(grib_context * c,unsigned char * buf,size_t nvals,int bytes,double * val)461 int grib_ieee_decode_array(grib_context* c, unsigned char* buf, size_t nvals, int bytes, double* val)
462 {
463 int err = 0, i = 0, j = 0;
464 unsigned char s[8] = {0,};
465 float fval;
466 double* pval = val;
467
468 switch (bytes) {
469 case 4:
470 for (i = 0; i < nvals; i++) {
471 #if IEEE_LE
472 for (j = 3; j >= 0; j--)
473 s[j] = *(buf++);
474 memcpy(&fval, s, 4);
475 val[i] = (double)fval;
476 #elif IEEE_BE
477 memcpy(&fval, buf, 4);
478 val[i] = (double)fval;
479 buf += 4;
480 #endif
481 }
482 break;
483 case 8:
484 for (i = 0; i < nvals; i++) {
485 #if IEEE_LE
486 for (j = 7; j >= 0; j--)
487 s[j] = *(buf++);
488 memcpy(pval++, s, 8);
489 #elif IEEE_BE
490 memcpy(pval++, buf, 8);
491 buf += 8;
492 #endif
493 }
494 break;
495 default:
496 grib_context_log(c, GRIB_LOG_ERROR,
497 "grib_ieee_decode_array: %d bits not implemented", bytes * 8);
498 return GRIB_NOT_IMPLEMENTED;
499 }
500
501 return err;
502 }
503
504 #else
505
grib_ieee_decode_array(grib_context * c,unsigned char * buf,size_t nvals,int bytes,double * val)506 int grib_ieee_decode_array(grib_context* c, unsigned char* buf, size_t nvals, int bytes, double* val)
507 {
508 int err = 0, i = 0;
509 long bitr = 0;
510
511 for (i = 0; i < nvals; i++)
512 val[i] = grib_long_to_ieee(grib_decode_unsigned_long(buf, &bitr, bytes * 8));
513
514 return err;
515 }
516
517 #endif
518
519 #ifdef IEEE
520
grib_ieee_encode_array(grib_context * c,double * val,size_t nvals,int bytes,unsigned char * buf)521 int grib_ieee_encode_array(grib_context* c, double* val, size_t nvals, int bytes,
522 unsigned char* buf)
523 {
524 int err = 0, i = 0, j = 0;
525 unsigned char s4[4];
526 unsigned char s8[8];
527 float fval = 0;
528 double* pval = val;
529
530 switch (bytes) {
531 case 4:
532 for (i = 0; i < nvals; i++) {
533 fval = (float)val[i];
534
535 #if IEEE_LE
536 memcpy(s4, &(fval), 4);
537 for (j = 3; j >= 0; j--)
538 *(buf++) = s4[j];
539 #elif IEEE_BE
540 memcpy(buf, &(fval), 4);
541 buf += 4;
542 #endif
543 }
544 break;
545 case 8:
546 for (i = 0; i < nvals; i++) {
547 #if IEEE_LE
548 memcpy(s8, pval++, 8);
549 for (j = 7; j >= 0; j--)
550 *(buf++) = s8[j];
551 #elif IEEE_BE
552 memcpy(buf, pval++, 8);
553 buf += 8;
554 #endif
555 }
556 break;
557 default:
558 grib_context_log(c, GRIB_LOG_ERROR,
559 "grib_ieee_encode_array: %d bits not implemented", bytes * 8);
560 return GRIB_NOT_IMPLEMENTED;
561 }
562
563 return err;
564 }
565
566 #else
567
grib_ieee_encode_array(grib_context * c,double * val,size_t nvals,int bytes,unsigned char * buf)568 int grib_ieee_encode_array(grib_context* c, double* val, size_t nvals, int bytes, unsigned char* buf)
569 {
570 int err = 0, i = 0;
571 long bitr = 0;
572
573 for (i = 0; i < nvals; i++)
574 grib_encode_unsigned_long(buf, grib_ieee_to_long(val[i]), &bitr, bytes * 8);
575
576 return err;
577 }
578
579 #endif
580