Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: CPL
4 : * Purpose: Convert between VAX and IEEE floating point formats
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2000, Avenza Systems Inc, http://www.avenza.com/
9 : *
10 : * Permission is hereby granted, free of charge, to any person obtaining a
11 : * copy of this software and associated documentation files (the "Software"),
12 : * to deal in the Software without restriction, including without limitation
13 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
14 : * and/or sell copies of the Software, and to permit persons to whom the
15 : * Software is furnished to do so, subject to the following conditions:
16 : *
17 : * The above copyright notice and this permission notice shall be included
18 : * in all copies or substantial portions of the Software.
19 : *
20 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26 : * DEALINGS IN THE SOFTWARE.
27 : ****************************************************************************/
28 :
29 : #include "cpl_port.h"
30 : #include "cpl_vax.h"
31 :
32 : namespace
33 : {
34 : typedef struct dbl
35 : {
36 : // cppcheck-suppress unusedStructMember
37 : GUInt32 hi;
38 : // cppcheck-suppress unusedStructMember
39 : GUInt32 lo;
40 : } double64_t;
41 : } // namespace
42 :
43 : /************************************************************************/
44 : /* CPLVaxToIEEEDouble() */
45 : /************************************************************************/
46 :
47 12107 : void CPLVaxToIEEEDouble(void *dbl)
48 :
49 : {
50 : double64_t dt;
51 : GUInt32 sign;
52 : int exponent;
53 : GUInt32 rndbits;
54 :
55 : /* -------------------------------------------------------------------- */
56 : /* Arrange the VAX double so that it may be accessed by a */
57 : /* double64_t structure, (two GUInt32s). */
58 : /* -------------------------------------------------------------------- */
59 : {
60 12107 : const unsigned char *src = static_cast<const unsigned char *>(dbl);
61 : unsigned char dest[8];
62 : #ifdef CPL_LSB
63 12107 : dest[2] = src[0];
64 12107 : dest[3] = src[1];
65 12107 : dest[0] = src[2];
66 12107 : dest[1] = src[3];
67 12107 : dest[6] = src[4];
68 12107 : dest[7] = src[5];
69 12107 : dest[4] = src[6];
70 12107 : dest[5] = src[7];
71 : #else
72 : dest[1] = src[0];
73 : dest[0] = src[1];
74 : dest[3] = src[2];
75 : dest[2] = src[3];
76 : dest[5] = src[4];
77 : dest[4] = src[5];
78 : dest[7] = src[6];
79 : dest[6] = src[7];
80 : #endif
81 12107 : memcpy(&dt, dest, 8);
82 : }
83 :
84 : /* -------------------------------------------------------------------- */
85 : /* Save the sign of the double */
86 : /* -------------------------------------------------------------------- */
87 12107 : sign = dt.hi & 0x80000000;
88 :
89 : /* -------------------------------------------------------------------- */
90 : /* Adjust the exponent so that we may work with it */
91 : /* -------------------------------------------------------------------- */
92 12107 : exponent = (dt.hi >> 23) & 0x000000ff;
93 :
94 12107 : if (exponent)
95 3894 : exponent = exponent - 129 + 1023;
96 :
97 : /* -------------------------------------------------------------------- */
98 : /* Save the bits that we are discarding so we can round properly */
99 : /* -------------------------------------------------------------------- */
100 12107 : rndbits = dt.lo & 0x00000007;
101 :
102 12107 : dt.lo = dt.lo >> 3;
103 12107 : dt.lo = (dt.lo & 0x1fffffff) | (dt.hi << 29);
104 :
105 12107 : if (rndbits)
106 0 : dt.lo = dt.lo | 0x00000001;
107 :
108 : /* -------------------------------------------------------------------- */
109 : /* Shift the hi-order int over 3 and insert the exponent and sign */
110 : /* -------------------------------------------------------------------- */
111 12107 : dt.hi = dt.hi >> 3;
112 12107 : dt.hi = dt.hi & 0x000fffff;
113 12107 : dt.hi = dt.hi | (static_cast<GUInt32>(exponent) << 20) | sign;
114 :
115 : #ifdef CPL_LSB
116 : /* -------------------------------------------------------------------- */
117 : /* Change the number to a byte swapped format */
118 : /* -------------------------------------------------------------------- */
119 12107 : const unsigned char *src = reinterpret_cast<const unsigned char *>(&dt);
120 12107 : unsigned char *dest = static_cast<unsigned char *>(dbl);
121 :
122 12107 : memcpy(dest + 0, src + 4, 4);
123 12107 : memcpy(dest + 4, src + 0, 4);
124 : #else
125 : memcpy(dbl, &dt, 8);
126 : #endif
127 12107 : }
128 :
129 : /************************************************************************/
130 : /* CPLIEEEToVaxDouble() */
131 : /************************************************************************/
132 :
133 111 : void CPLIEEEToVaxDouble(void *dbl)
134 :
135 : {
136 : double64_t dt;
137 :
138 : #ifdef CPL_LSB
139 : {
140 111 : const GByte *src = static_cast<const GByte *>(dbl);
141 : GByte dest[8];
142 :
143 111 : dest[0] = src[4];
144 111 : dest[1] = src[5];
145 111 : dest[2] = src[6];
146 111 : dest[3] = src[7];
147 111 : dest[4] = src[0];
148 111 : dest[5] = src[1];
149 111 : dest[6] = src[2];
150 111 : dest[7] = src[3];
151 111 : memcpy(&dt, dest, 8);
152 : }
153 : #else
154 : memcpy(&dt, dbl, 8);
155 : #endif
156 :
157 111 : GInt32 sign = dt.hi & 0x80000000;
158 111 : GInt32 exponent = dt.hi >> 20;
159 111 : exponent = exponent & 0x000007ff;
160 :
161 : /* -------------------------------------------------------------------- */
162 : /* An exponent of zero means a zero value. */
163 : /* -------------------------------------------------------------------- */
164 111 : if (exponent)
165 110 : exponent = exponent - 1023 + 129;
166 :
167 : /* -------------------------------------------------------------------- */
168 : /* In the case of overflow, return the largest number we can */
169 : /* -------------------------------------------------------------------- */
170 111 : if (exponent > 255)
171 : {
172 : GByte dest[8];
173 :
174 0 : if (sign)
175 0 : dest[1] = 0xff;
176 : else
177 0 : dest[1] = 0x7f;
178 :
179 0 : dest[0] = 0xff;
180 0 : dest[2] = 0xff;
181 0 : dest[3] = 0xff;
182 0 : dest[4] = 0xff;
183 0 : dest[5] = 0xff;
184 0 : dest[6] = 0xff;
185 0 : dest[7] = 0xff;
186 0 : memcpy(dbl, dest, 8);
187 :
188 0 : return;
189 : }
190 :
191 : /* -------------------------------------------------------------------- */
192 : /* In the case of of underflow return zero */
193 : /* -------------------------------------------------------------------- */
194 111 : else if ((exponent < 0) || (exponent == 0 && sign == 0))
195 : {
196 1 : memset(dbl, 0, 8);
197 :
198 1 : return;
199 : }
200 : else
201 : {
202 : /* --------------------------------------------------------------------
203 : */
204 : /* Shift the fraction 3 bits left and set the exponent and
205 : * sign*/
206 : /* --------------------------------------------------------------------
207 : */
208 110 : dt.hi = dt.hi << 3;
209 110 : dt.hi = dt.hi | (dt.lo >> 29);
210 110 : dt.hi = dt.hi & 0x007fffff;
211 110 : dt.hi = dt.hi | (exponent << 23) | sign;
212 :
213 110 : dt.lo = dt.lo << 3;
214 : }
215 :
216 : /* -------------------------------------------------------------------- */
217 : /* Convert the double back to VAX format */
218 : /* -------------------------------------------------------------------- */
219 110 : const GByte *src = reinterpret_cast<GByte *>(&dt);
220 :
221 : #ifdef CPL_LSB
222 110 : GByte *dest = static_cast<GByte *>(dbl);
223 110 : memcpy(dest + 2, src + 0, 2);
224 110 : memcpy(dest + 0, src + 2, 2);
225 110 : memcpy(dest + 6, src + 4, 2);
226 110 : memcpy(dest + 4, src + 6, 2);
227 : #else
228 : GByte dest[8];
229 : dest[1] = src[0];
230 : dest[0] = src[1];
231 : dest[3] = src[2];
232 : dest[2] = src[3];
233 : dest[5] = src[4];
234 : dest[4] = src[5];
235 : dest[7] = src[6];
236 : dest[6] = src[7];
237 : memcpy(dbl, dest, 8);
238 : #endif
239 : }
240 :
241 : //////////////////////////////////////////////////////////////////////////
242 : /// Below code is adapted from Public Domain VICAR project
243 : /// https://github.com/nasa/VICAR/blob/master/vos/rtl/source/conv_vax_ieee_r.c
244 : //////////////////////////////////////////////////////////////////////////
245 :
246 192 : static void real_byte_swap(const unsigned char from[4], unsigned char to[4])
247 : {
248 192 : to[0] = from[1];
249 192 : to[1] = from[0];
250 192 : to[2] = from[3];
251 192 : to[3] = from[2];
252 192 : }
253 :
254 : /* Shift x[1]..x[3] right one bit by bytes, don't bother with x[0] */
255 : #define SHIFT_RIGHT(x) \
256 : { \
257 : x[3] = ((x[3] >> 1) & 0x7F) | ((x[2] << 7) & 0x80); \
258 : x[2] = ((x[2] >> 1) & 0x7F) | ((x[1] << 7) & 0x80); \
259 : x[1] = (x[1] >> 1) & 0x7F; \
260 : }
261 :
262 : /* Shift x[1]..x[3] left one bit by bytes, don't bother with x[0] */
263 : #define SHIFT_LEFT(x) \
264 : { \
265 : x[1] = ((x[1] << 1) & 0xFE) | ((x[2] >> 7) & 0x01); \
266 : x[2] = ((x[2] << 1) & 0xFE) | ((x[3] >> 7) & 0x01); \
267 : x[3] = (x[3] << 1) & 0xFE; \
268 : }
269 :
270 : /************************************************************************/
271 : /* Convert between IEEE and Vax single-precision floating point. */
272 : /* Both formats are represented as: */
273 : /* (-1)^s * f * 2^(e-bias) */
274 : /* where s is the sign bit, f is the mantissa (see below), e is the */
275 : /* exponent, and bias is the exponent bias (see below). */
276 : /* There is an assumed leading 1 on the mantissa (except for IEEE */
277 : /* denormalized numbers), but the placement of the binary point varies. */
278 : /* */
279 : /* IEEE format: seeeeeee efffffff 8*f 8*f */
280 : /* where e is exponent with bias of 127 and f is of the */
281 : /* form 1.fffff... */
282 : /* Special cases: */
283 : /* e=255, f!=0: NaN (Not a Number) */
284 : /* e=255, f=0: Infinity (+/- depending on s) */
285 : /* e=0, f!=0: Denormalized numbers, of the form */
286 : /* (-1)^s * (0.ffff) * 2^(-126) */
287 : /* e=0, f=0: Zero (can be +/-) */
288 : /* */
289 : /* VAX format: seeeeeee efffffff 8*f 8*f */
290 : /* where e is exponent with bias of 128 and f is of the */
291 : /* form .1fffff... */
292 : /* Byte swapping: Note that the above format is the logical format, */
293 : /* which can be represented as bytes SE1 E2F1 F2 F3. */
294 : /* The actual order in memory is E2F1 SE1 F3 F2 (which is */
295 : /* two half-word swaps, NOT a full-word swap). */
296 : /* Special cases: */
297 : /* e=0, s=0: Zero (no +/-) */
298 : /* e=0, s=1: Invalid, causes Reserved Operand error */
299 : /* */
300 : /* The same code works on all byte-order machines because only byte */
301 : /* operations are performed. It could perhaps be done more efficiently */
302 : /* on a longword basis, but then the code would be byte-order dependent.*/
303 : /* MAKE SURE any mods will work on either byte order!!! */
304 : /************************************************************************/
305 :
306 : /************************************************************************/
307 : /* This routine will convert VAX F floating point values to IEEE */
308 : /* single precision floating point. */
309 : /************************************************************************/
310 :
311 156 : static void vax_ieee_r(const unsigned char *from, unsigned char *ieee)
312 : {
313 : unsigned char vaxf[4];
314 : unsigned char exp;
315 :
316 156 : real_byte_swap(from, vaxf); /* Put bytes in rational order */
317 156 : memcpy(ieee, vaxf, 4); /* Since most bits are the same */
318 :
319 156 : exp = ((vaxf[0] << 1) & 0xFE) | ((vaxf[1] >> 7) & 0x01);
320 :
321 156 : if (exp == 0)
322 : { /* Zero or invalid pattern */
323 0 : if (vaxf[0] & 0x80)
324 : { /* Sign bit set, which is illegal for VAX */
325 0 : ieee[0] = 0x7F; /* IEEE NaN */
326 0 : ieee[1] = 0xFF;
327 0 : ieee[2] = 0xFF;
328 0 : ieee[3] = 0xFF;
329 : }
330 : else
331 : { /* Zero */
332 0 : ieee[0] = ieee[1] = ieee[2] = ieee[3] = 0;
333 : }
334 : }
335 :
336 156 : else if (exp >= 3)
337 : { /* Normal case */
338 156 : exp -= 2;
339 156 : ieee[0] =
340 156 : (vaxf[0] & 0x80) | ((exp >> 1) & 0x7F); /* remake sign + exponent */
341 : } /* Low bit of exp can't change, so don't bother w/it */
342 :
343 0 : else if (exp == 2)
344 : { /* Denormalize the number */
345 0 : SHIFT_RIGHT(ieee); /* Which means shift right 1, */
346 0 : ieee[1] = (ieee[1] & 0x3F) | 0x40; /* Add suppressed most signif bit, */
347 0 : ieee[0] = vaxf[0] & 0x80; /* and set exponent to 0 (preserving sign) */
348 : }
349 :
350 : else
351 : { /* Exp==1, denormalize again */
352 0 : SHIFT_RIGHT(ieee); /* Like above but shift by 2 */
353 0 : SHIFT_RIGHT(ieee);
354 0 : ieee[1] = (ieee[1] & 0x1F) | 0x20;
355 0 : ieee[0] = vaxf[0] & 0x80;
356 : }
357 :
358 : #ifdef CPL_LSB
359 156 : CPL_SWAP32PTR(ieee);
360 : #endif
361 156 : }
362 :
363 : /************************************************************************/
364 : /* This routine will convert IEEE single precision floating point */
365 : /* values to VAX F floating point. */
366 : /************************************************************************/
367 :
368 36 : static void ieee_vax_r(unsigned char *ieee, unsigned char *to)
369 : {
370 : unsigned char vaxf[4];
371 : unsigned char exp;
372 :
373 : #ifdef CPL_LSB
374 36 : CPL_SWAP32PTR(ieee);
375 : #endif
376 :
377 36 : memcpy(vaxf, ieee, 4); /* Since most bits are the same */
378 :
379 36 : exp = ((ieee[0] << 1) & 0xFE) | ((ieee[1] >> 7) & 0x01);
380 :
381 : /* Exponent 255 means NaN or Infinity, exponent 254 is too large for */
382 : /* VAX notation. In either case, set to sign * highest possible number */
383 :
384 36 : if (exp == 255 || exp == 254)
385 : { /* Infinity or NaN or too big */
386 0 : vaxf[0] = 0x7F | (ieee[0] & 0x80);
387 0 : vaxf[1] = 0xFF;
388 0 : vaxf[2] = 0xFF;
389 0 : vaxf[3] = 0xFF;
390 : }
391 :
392 36 : else if (exp != 0)
393 : { /* Normal case */
394 36 : exp += 2;
395 36 : vaxf[0] =
396 36 : (ieee[0] & 0x80) | ((exp >> 1) & 0x7F); /* remake sign + exponent */
397 : } /* Low bit of exp can't change, so don't bother w/it */
398 :
399 : else
400 : { /* exp == 0, zero or denormalized number */
401 0 : if (ieee[1] == 0 && ieee[2] == 0 && ieee[3] == 0)
402 : { /* +/- 0 */
403 0 : vaxf[0] = vaxf[1] = vaxf[2] = vaxf[3] = 0;
404 : }
405 : else
406 : { /* denormalized number */
407 0 : if (ieee[1] & 0x40)
408 : { /* hi bit set (0.1ffff) */
409 0 : SHIFT_LEFT(vaxf); /* Renormalize */
410 0 : vaxf[1] = vaxf[1] & 0x7F; /* Set vax exponent to 2 */
411 0 : vaxf[0] = (ieee[0] & 0x80) | 0x01; /* sign, exponent==2 */
412 : }
413 0 : else if (ieee[1] & 0x20)
414 : { /* next bit set (0.01ffff) */
415 0 : SHIFT_LEFT(vaxf); /* Renormalize */
416 0 : SHIFT_LEFT(vaxf);
417 0 : vaxf[1] = vaxf[1] | 0x80; /* Set vax exponent to 1 */
418 0 : vaxf[0] = ieee[0] & 0x80; /* sign, exponent==1 */
419 : }
420 : else
421 : { /* Number too small for VAX */
422 0 : vaxf[0] = vaxf[1] = vaxf[2] = vaxf[3] = 0; /* so set to 0 */
423 : }
424 : }
425 : }
426 :
427 36 : real_byte_swap(vaxf, to); /* Put bytes in weird VAX order */
428 36 : }
429 :
430 156 : void CPLVaxToIEEEFloat(void *f)
431 : {
432 : unsigned char res[4];
433 156 : vax_ieee_r(static_cast<const unsigned char *>(f), res);
434 156 : memcpy(f, res, 4);
435 156 : }
436 :
437 36 : void CPLIEEEToVaxFloat(void *f)
438 : {
439 : unsigned char res[4];
440 36 : ieee_vax_r(static_cast<unsigned char *>(f), res);
441 36 : memcpy(f, res, 4);
442 36 : }
|
