1: /*
2: * /src/NTP/ntp4-dev/libntp/ieee754io.c,v 4.12 2005/04/16 17:32:10 kardel RELEASE_20050508_A
3: *
4: * ieee754io.c,v 4.12 2005/04/16 17:32:10 kardel RELEASE_20050508_A
5: *
6: * $Created: Sun Jul 13 09:12:02 1997 $
7: *
8: * Copyright (c) 1997-2005 by Frank Kardel <kardel <AT> ntp.org>
9: *
10: * Redistribution and use in source and binary forms, with or without
11: * modification, are permitted provided that the following conditions
12: * are met:
13: * 1. Redistributions of source code must retain the above copyright
14: * notice, this list of conditions and the following disclaimer.
15: * 2. Redistributions in binary form must reproduce the above copyright
16: * notice, this list of conditions and the following disclaimer in the
17: * documentation and/or other materials provided with the distribution.
18: * 3. Neither the name of the author nor the names of its contributors
19: * may be used to endorse or promote products derived from this software
20: * without specific prior written permission.
21: *
22: * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
23: * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
24: * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
25: * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
26: * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
27: * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
28: * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
29: * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
30: * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
31: * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
32: * SUCH DAMAGE.
33: *
34: */
35:
36: #ifdef HAVE_CONFIG_H
37: #include "config.h"
38: #endif
39:
40: #include <stdio.h>
41: #include "l_stdlib.h"
42: #include "ntp_stdlib.h"
43: #include "ntp_fp.h"
44: #include "ieee754io.h"
45:
46: static unsigned char get_byte (unsigned char *, offsets_t, int *);
47: #ifdef __not_yet__
48: static void put_byte (unsigned char *, offsets_t, int *, unsigned char);
49: #endif
50:
51: #ifdef LIBDEBUG
52:
53: #include "lib_strbuf.h"
54:
55: static char *
56: fmt_blong(
57: unsigned long val,
58: int cnt
59: )
60: {
61: char *buf, *s;
62: int i = cnt;
63:
64: val <<= 32 - cnt;
65: LIB_GETBUF(buf);
66: s = buf;
67:
68: while (i--)
69: {
70: if (val & 0x80000000)
71: {
72: *s++ = '1';
73: }
74: else
75: {
76: *s++ = '0';
77: }
78: val <<= 1;
79: }
80: *s = '\0';
81: return buf;
82: }
83:
84: static char *
85: fmt_flt(
86: unsigned int sign,
87: unsigned long mh,
88: unsigned long ml,
89: unsigned long ch
90: )
91: {
92: char *buf;
93:
94: LIB_GETBUF(buf);
95: sprintf(buf, "%c %s %s %s", sign ? '-' : '+',
96: fmt_blong(ch, 11),
97: fmt_blong(mh, 20),
98: fmt_blong(ml, 32));
99: return buf;
100: }
101:
102: static char *
103: fmt_hex(
104: unsigned char *bufp,
105: int length
106: )
107: {
108: char *buf;
109: int i;
110:
111: LIB_GETBUF(buf);
112: for (i = 0; i < length; i++)
113: {
114: sprintf(buf+i*2, "%02x", bufp[i]);
115: }
116: return buf;
117: }
118:
119: #endif
120:
121: static unsigned char
122: get_byte(
123: unsigned char *bufp,
124: offsets_t offset,
125: int *fieldindex
126: )
127: {
128: unsigned char val;
129:
130: val = *(bufp + offset[*fieldindex]);
131: #ifdef LIBDEBUG
132: if (debug > 4)
133: printf("fetchieee754: getbyte(0x%08x, %d) = 0x%02x\n", (unsigned int)(bufp)+offset[*fieldindex], *fieldindex, val);
134: #endif
135: (*fieldindex)++;
136: return val;
137: }
138:
139: #ifdef __not_yet__
140: static void
141: put_byte(
142: unsigned char *bufp,
143: offsets_t offsets,
144: int *fieldindex,
145: unsigned char val
146: )
147: {
148: *(bufp + offsets[*fieldindex]) = val;
149: (*fieldindex)++;
150: }
151: #endif
152:
153: /*
154: * make conversions to and from external IEEE754 formats and internal
155: * NTP FP format.
156: */
157: int
158: fetch_ieee754(
159: unsigned char **buffpp,
160: int size,
161: l_fp *lfpp,
162: offsets_t offsets
163: )
164: {
165: unsigned char *bufp = *buffpp;
166: unsigned int sign;
167: unsigned int bias;
168: unsigned int maxexp;
169: int mbits;
170: u_long mantissa_low;
171: u_long mantissa_high;
172: u_long characteristic;
173: long exponent;
174: #ifdef LIBDEBUG
175: int length;
176: #endif
177: unsigned char val;
178: int fieldindex = 0;
179:
180: switch (size)
181: {
182: case IEEE_DOUBLE:
183: #ifdef LIBDEBUG
184: length = 8;
185: #endif
186: mbits = 52;
187: bias = 1023;
188: maxexp = 2047;
189: break;
190:
191: case IEEE_SINGLE:
192: #ifdef LIBDEBUG
193: length = 4;
194: #endif
195: mbits = 23;
196: bias = 127;
197: maxexp = 255;
198: break;
199:
200: default:
201: return IEEE_BADCALL;
202: }
203:
204: val = get_byte(bufp, offsets, &fieldindex); /* fetch sign byte & first part of characteristic */
205:
206: sign = (val & 0x80) != 0;
207: characteristic = (val & 0x7F);
208:
209: val = get_byte(bufp, offsets, &fieldindex); /* fetch rest of characteristic and start of mantissa */
210:
211: switch (size)
212: {
213: case IEEE_SINGLE:
214: characteristic <<= 1;
215: characteristic |= (val & 0x80) != 0; /* grab last characteristic bit */
216:
217: mantissa_high = 0;
218:
219: mantissa_low = (val &0x7F) << 16;
220: mantissa_low |= get_byte(bufp, offsets, &fieldindex) << 8;
221: mantissa_low |= get_byte(bufp, offsets, &fieldindex);
222: break;
223:
224: case IEEE_DOUBLE:
225: characteristic <<= 4;
226: characteristic |= (val & 0xF0) >> 4; /* grab lower characteristic bits */
227:
228: mantissa_high = (val & 0x0F) << 16;
229: mantissa_high |= get_byte(bufp, offsets, &fieldindex) << 8;
230: mantissa_high |= get_byte(bufp, offsets, &fieldindex);
231:
232: mantissa_low = get_byte(bufp, offsets, &fieldindex) << 24;
233: mantissa_low |= get_byte(bufp, offsets, &fieldindex) << 16;
234: mantissa_low |= get_byte(bufp, offsets, &fieldindex) << 8;
235: mantissa_low |= get_byte(bufp, offsets, &fieldindex);
236: break;
237:
238: default:
239: return IEEE_BADCALL;
240: }
241: #ifdef LIBDEBUG
242: if (debug > 4)
243: {
244: double d;
245: float f;
246:
247: if (size == IEEE_SINGLE)
248: {
249: int i;
250:
251: for (i = 0; i < length; i++)
252: {
253: *((unsigned char *)(&f)+i) = *(*buffpp + offsets[i]);
254: }
255: d = f;
256: }
257: else
258: {
259: int i;
260:
261: for (i = 0; i < length; i++)
262: {
263: *((unsigned char *)(&d)+i) = *(*buffpp + offsets[i]);
264: }
265: }
266:
267: printf("fetchieee754: FP: %s -> %s -> %e(=%s)\n", fmt_hex(*buffpp, length),
268: fmt_flt(sign, mantissa_high, mantissa_low, characteristic),
269: d, fmt_hex((unsigned char *)&d, length));
270: }
271: #endif
272:
273: *buffpp += fieldindex;
274:
275: /*
276: * detect funny numbers
277: */
278: if (characteristic == maxexp)
279: {
280: /*
281: * NaN or Infinity
282: */
283: if (mantissa_low || mantissa_high)
284: {
285: /*
286: * NaN
287: */
288: return IEEE_NAN;
289: }
290: else
291: {
292: /*
293: * +Inf or -Inf
294: */
295: return sign ? IEEE_NEGINFINITY : IEEE_POSINFINITY;
296: }
297: }
298: else
299: {
300: /*
301: * collect real numbers
302: */
303:
304: L_CLR(lfpp);
305:
306: /*
307: * check for overflows
308: */
309: exponent = characteristic - bias;
310:
311: if (exponent > 31) /* sorry - hardcoded */
312: {
313: /*
314: * overflow only in respect to NTP-FP representation
315: */
316: return sign ? IEEE_NEGOVERFLOW : IEEE_POSOVERFLOW;
317: }
318: else
319: {
320: int frac_offset; /* where the fraction starts */
321:
322: frac_offset = mbits - exponent;
323:
324: if (characteristic == 0)
325: {
326: /*
327: * de-normalized or tiny number - fits only as 0
328: */
329: return IEEE_OK;
330: }
331: else
332: {
333: /*
334: * adjust for implied 1
335: */
336: if (mbits > 31)
337: mantissa_high |= 1 << (mbits - 32);
338: else
339: mantissa_low |= 1 << mbits;
340:
341: /*
342: * take mantissa apart - if only all machine would support
343: * 64 bit operations 8-(
344: */
345: if (frac_offset > mbits)
346: {
347: lfpp->l_ui = 0; /* only fractional number */
348: frac_offset -= mbits + 1; /* will now contain right shift count - 1*/
349: if (mbits > 31)
350: {
351: lfpp->l_uf = mantissa_high << (63 - mbits);
352: lfpp->l_uf |= mantissa_low >> (mbits - 33);
353: lfpp->l_uf >>= frac_offset;
354: }
355: else
356: {
357: lfpp->l_uf = mantissa_low >> frac_offset;
358: }
359: }
360: else
361: {
362: if (frac_offset > 32)
363: {
364: /*
365: * must split in high word
366: */
367: lfpp->l_ui = mantissa_high >> (frac_offset - 32);
368: lfpp->l_uf = (mantissa_high & ((1 << (frac_offset - 32)) - 1)) << (64 - frac_offset);
369: lfpp->l_uf |= mantissa_low >> (frac_offset - 32);
370: }
371: else
372: {
373: /*
374: * must split in low word
375: */
376: lfpp->l_ui = mantissa_high << (32 - frac_offset);
377: lfpp->l_ui |= (mantissa_low >> frac_offset) & ((1 << (32 - frac_offset)) - 1);
378: lfpp->l_uf = (mantissa_low & ((1 << frac_offset) - 1)) << (32 - frac_offset);
379: }
380: }
381:
382: /*
383: * adjust for sign
384: */
385: if (sign)
386: {
387: L_NEG(lfpp);
388: }
389:
390: return IEEE_OK;
391: }
392: }
393: }
394: }
395:
396: int
397: put_ieee754(
398: unsigned char **bufpp,
399: int size,
400: l_fp *lfpp,
401: offsets_t offsets
402: )
403: {
404: l_fp outlfp;
405: #ifdef LIBDEBUG
406: unsigned int sign;
407: unsigned int bias;
408: #endif
409: /*unsigned int maxexp;*/
410: int mbits;
411: int msb;
412: u_long mantissa_low = 0;
413: u_long mantissa_high = 0;
414: #ifdef LIBDEBUG
415: u_long characteristic = 0;
416: long exponent;
417: #endif
418: /*int length;*/
419: unsigned long mask;
420:
421: outlfp = *lfpp;
422:
423: switch (size)
424: {
425: case IEEE_DOUBLE:
426: /*length = 8;*/
427: mbits = 52;
428: #ifdef LIBDEBUG
429: bias = 1023;
430: #endif
431: /*maxexp = 2047;*/
432: break;
433:
434: case IEEE_SINGLE:
435: /*length = 4;*/
436: mbits = 23;
437: #ifdef LIBDEBUG
438: bias = 127;
439: #endif
440: /*maxexp = 255;*/
441: break;
442:
443: default:
444: return IEEE_BADCALL;
445: }
446:
447: /*
448: * find sign
449: */
450: if (L_ISNEG(&outlfp))
451: {
452: L_NEG(&outlfp);
453: #ifdef LIBDEBUG
454: sign = 1;
455: #endif
456: }
457: else
458: {
459: #ifdef LIBDEBUG
460: sign = 0;
461: #endif
462: }
463:
464: if (L_ISZERO(&outlfp))
465: {
466: #ifdef LIBDEBUG
467: exponent = mantissa_high = mantissa_low = 0; /* true zero */
468: #endif
469: }
470: else
471: {
472: /*
473: * find number of significant integer bits
474: */
475: mask = 0x80000000;
476: if (outlfp.l_ui)
477: {
478: msb = 63;
479: while (mask && ((outlfp.l_ui & mask) == 0))
480: {
481: mask >>= 1;
482: msb--;
483: }
484: }
485: else
486: {
487: msb = 31;
488: while (mask && ((outlfp.l_uf & mask) == 0))
489: {
490: mask >>= 1;
491: msb--;
492: }
493: }
494:
495: switch (size)
496: {
497: case IEEE_SINGLE:
498: mantissa_high = 0;
499: if (msb >= 32)
500: {
501: mantissa_low = (outlfp.l_ui & ((1 << (msb - 32)) - 1)) << (mbits - (msb - 32));
502: mantissa_low |= outlfp.l_uf >> (mbits - (msb - 32));
503: }
504: else
505: {
506: mantissa_low = (outlfp.l_uf << (mbits - msb)) & ((1 << mbits) - 1);
507: }
508: break;
509:
510: case IEEE_DOUBLE:
511: if (msb >= 32)
512: {
513: mantissa_high = (outlfp.l_ui << (mbits - msb)) & ((1 << (mbits - 32)) - 1);
514: mantissa_high |= outlfp.l_uf >> (32 - (mbits - msb));
515: mantissa_low = (outlfp.l_ui & ((1 << (msb - mbits)) - 1)) << (32 - (msb - mbits));
516: mantissa_low |= outlfp.l_uf >> (msb - mbits);
517: }
518: else
519: {
520: mantissa_high = outlfp.l_uf << (mbits - 32 - msb);
521: mantissa_low = outlfp.l_uf << (mbits - 32);
522: }
523: }
524:
525: #ifdef LIBDEBUG
526: exponent = msb - 32;
527: characteristic = exponent + bias;
528:
529: if (debug > 4)
530: printf("FP: %s\n", fmt_flt(sign, mantissa_high, mantissa_low, characteristic));
531: #endif
532: }
533: return IEEE_OK;
534: }
535:
536:
537: #if defined(DEBUG) && defined(LIBDEBUG)
538: int main(
539: int argc,
540: char **argv
541: )
542: {
543: static offsets_t native_off = { 0, 1, 2, 3, 4, 5, 6, 7 };
544: double f = 1.0;
545: double *f_p = &f;
546: l_fp fp;
547:
548: if (argc == 2)
549: {
550: if (sscanf(argv[1], "%lf", &f) != 1)
551: {
552: printf("cannot convert %s to a float\n", argv[1]);
553: return 1;
554: }
555: }
556:
557: printf("double: %s %s\n", fmt_blong(*(unsigned long *)&f, 32), fmt_blong(*(unsigned long *)((char *)(&f)+4), 32));
558: printf("fetch from %f = %d\n", f, fetch_ieee754((void *)&f_p, IEEE_DOUBLE, &fp, native_off));
559: printf("fp [%s %s] = %s\n", fmt_blong(fp.l_ui, 32), fmt_blong(fp.l_uf, 32), mfptoa(fp.l_ui, fp.l_uf, 15));
560: f_p = &f;
561: put_ieee754((void *)&f_p, IEEE_DOUBLE, &fp, native_off);
562:
563: return 0;
564: }
565:
566: #endif
567: /*
568: * History:
569: *
570: * ieee754io.c,v
571: * Revision 4.12 2005/04/16 17:32:10 kardel
572: * update copyright
573: *
574: * Revision 4.11 2004/11/14 15:29:41 kardel
575: * support PPSAPI, upgrade Copyright to Berkeley style
576: *
577: * Revision 4.8 1999/02/21 12:17:36 kardel
578: * 4.91f reconcilation
579: *
580: * Revision 4.7 1999/02/21 11:26:03 kardel
581: * renamed index to fieldindex to avoid index() name clash
582: *
583: * Revision 4.6 1998/11/15 20:27:52 kardel
584: * Release 4.0.73e13 reconcilation
585: *
586: * Revision 4.5 1998/08/16 19:01:51 kardel
587: * debug information only compile for LIBDEBUG case
588: *
589: * Revision 4.4 1998/08/09 09:39:28 kardel
590: * Release 4.0.73e2 reconcilation
591: *
592: * Revision 4.3 1998/06/13 11:56:19 kardel
593: * disabled putbute() for the time being
594: *
595: * Revision 4.2 1998/06/12 15:16:58 kardel
596: * ansi2knr compatibility
597: *
598: * Revision 4.1 1998/05/24 07:59:56 kardel
599: * conditional debug support
600: *
601: * Revision 4.0 1998/04/10 19:46:29 kardel
602: * Start 4.0 release version numbering
603: *
604: * Revision 1.1 1998/04/10 19:27:46 kardel
605: * initial NTP VERSION 4 integration of PARSE with GPS166 binary support
606: *
607: * Revision 1.1 1997/10/06 21:05:45 kardel
608: * new parse structure
609: *
610: */
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>