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