File: C:\NOAA\NEMS_11731\src\chem\gocart\src\GMAO_Shared\MAPL_cfio\ShaveMantissa.c
1
2 #include <stdio.h>
3 #include <math.h>
4 #include "ShaveMantissa.h"
5
6 #define MAXBITS 20
7
8 #define EXPMSK 0x7f800000
9
10 #define isZNN(E) ((E)==0x00000000)
11 #define isNAN(E) ((E)==EXPMSK )
12
13
14 #if 1
15 # define isUndef(A) ((A) == undef)
16 #else
17 # define isUndef(A) ((float32)fabs(undef-(A))<tol)
18 #endif
19
20 #define SKIP(I,A) (isZNN(I) || isUndef(A) || isNAN(I))
21
22
23
24 (float32 minv, float32 maxv)
25 {
26 float32 midv, mnabs, range;
27
28 range = (maxv-minv);
29 midv = (maxv+minv)*0.5;
30 mnabs = abs(maxv)>abs(minv) ? fabs(minv) : fabs(maxv);
31
32 return (range<mnabs) ? midv : midv*(mnabs/range);
33 }
34
35
36
37
38
39
40
41
42
43
44
45 int ShaveMantissa32 ( a, ain, len, xbits, has_undef, undef, chunksize )
46
47
48
49
50
51 ;
52 int xbits;
53 int has_undef;
54 ;
55 ;
56 [];
57
58
59
60
61
62 [];
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101 {
102 float32 maxv, minv, offset, *b, *c, *begnxt, *last, tol;
103 uint32 round, mask, e;
104
105 union{
106 float32 x;
107 uint32 i;
108 } aa;
109
110
111
112 if ( len < 1 || xbits < 1 ) {
113 fprintf(stderr,
114 "ShaveMantissa32: Bad length of mask bits: len= %d, xbits= %d\n",
115 len, xbits );
116 return 1;
117 }
118
119 if ( xbits > MAXBITS ) {
120 fprintf(stderr,
121 "ShaveMantissa32: Shaving too many bits: %d; maximum allowed is %d\n",
122 xbits, MAXBITS );
123 return 2;
124 }
125
126
127
128 if ( !has_undef ) undef = (float32) HUGE_VAL;
129
130
131
132 = 0.0001*undef;
133 mask = 0xFFFFFFFF<<xbits--;
134 round = 1 <<xbits ;
135 last = &a[len-1];
136
137
138
139
140
141 = a;
142 if(ain!=a) {
143 if(abs(ain-a)<len) {
144 fprintf(stderr,"ShaveMantissa32: Overlapping arrays");
145 return 3;
146 }
147 while(a<=last) *a++=*ain++;
148 }
149
150
151
152 while(b<=last) {
153
154
155
156 = b + chunksize;
157 if(begnxt>last) begnxt = last+1;
158
159
160
161 = b-1;
162 while(++a < begnxt) {
163 aa.x = *a;
164 e = aa.i & EXPMSK;
165 if(!SKIP(e,aa.x)) {maxv=aa.x; minv=aa.x; c=a; break;}
166 }
167
168
169
170 if(a==begnxt) {b=begnxt; continue;}
171
172
173
174 while(++a<begnxt) {
175 aa.x = *a;
176 e = aa.i & EXPMSK;
177 if(!SKIP(e,aa.x)) {
178 if(aa.x<minv) minv=aa.x;
179 if(aa.x>maxv) maxv=aa.x;
180 }
181 }
182
183
184
185 if(minv==maxv) {b=begnxt; continue;}
186
187
188
189 = SetOffset(minv,maxv);
190
191
192
193 = c-1;
194 while(++a<begnxt) {
195 aa.x = *a;
196 e = aa.i & EXPMSK;
197 if(!SKIP(e,aa.x)) {
198 aa.x -= offset;
199 aa.i = ((aa.i + round) & mask);
200 aa.x += offset;
201 if(aa.x>maxv) aa.x=maxv;
202 if(aa.x<minv) aa.x=minv;
203 *a = aa.x;
204 }
205 }
206
207
208
209 = begnxt;
210
211 }
212
213 return 0;
214 }
215
216
217
218
219
220 int SHAVEMANTISSA32 (float32 *a, float32 *ain, int32 *len, int *xbits,
221 int *has_undef, float32 *undef, int32 *chunksize)
222 {return (int)ShaveMantissa32(a,ain,*len,*xbits,*has_undef,*undef,*chunksize);}
223
224 int SHAVEMANTISSA32_ (float32 *a, float32 *ain, int32 *len, int *xbits,
225 int *has_undef, float32 *undef, int32 *chunksize)
226 {return (int)ShaveMantissa32(a,ain,*len,*xbits,*has_undef,*undef,*chunksize);}
227
228 int shavemantissa32 (float32 *a, float32 *ain, int32 *len, int *xbits,
229 int *has_undef, float32 *undef, int32 *chunksize)
230 {return (int)ShaveMantissa32(a,ain,*len,*xbits,*has_undef,*undef,*chunksize);}
231
232 int shavemantissa32_ (float32 *a, float32 *ain, int32 *len, int *xbits,
233 int *has_undef, float32 *undef, int32 *chunksize)
234 {return (int)ShaveMantissa32(a,ain,*len,*xbits,*has_undef,*undef,*chunksize);}
235
236
237
238
239
240