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