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