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