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 static char *
fmt_blong(unsigned long val,int cnt)54 fmt_blong(
55 unsigned long val,
56 int cnt
57 )
58 {
59 char *buf, *s;
60 int i = cnt;
61
62 val <<= 32 - cnt;
63 LIB_GETBUF(buf);
64 s = buf;
65
66 while (i--)
67 {
68 if (val & 0x80000000)
69 {
70 *s++ = '1';
71 }
72 else
73 {
74 *s++ = '0';
75 }
76 val <<= 1;
77 }
78 *s = '\0';
79 return buf;
80 }
81
82 static char *
fmt_flt(unsigned int sign,unsigned long mh,unsigned long ml,unsigned long ch)83 fmt_flt(
84 unsigned int sign,
85 unsigned long mh,
86 unsigned long ml,
87 unsigned long ch
88 )
89 {
90 char *buf;
91
92 LIB_GETBUF(buf);
93 snprintf(buf, LIB_BUFLENGTH, "%c %s %s %s", sign ? '-' : '+',
94 fmt_blong(ch, 11),
95 fmt_blong(mh, 20),
96 fmt_blong(ml, 32));
97
98 return buf;
99 }
100
101 static char *
fmt_hex(unsigned char * bufp,int length)102 fmt_hex(
103 unsigned char *bufp,
104 int length
105 )
106 {
107 char * buf;
108 char hex[4];
109 int i;
110
111 LIB_GETBUF(buf);
112 buf[0] = '\0';
113 for (i = 0; i < length; i++) {
114 snprintf(hex, sizeof(hex), "%02x", bufp[i]);
115 strlcat(buf, hex, LIB_BUFLENGTH);
116 }
117
118 return buf;
119 }
120
121 #endif
122
123 static unsigned char
get_byte(unsigned char * bufp,offsets_t offset,int * fieldindex)124 get_byte(
125 unsigned char *bufp,
126 offsets_t offset,
127 int *fieldindex
128 )
129 {
130 unsigned char val;
131
132 val = *(bufp + offset[*fieldindex]);
133 #ifdef LIBDEBUG
134 if (debug > 4)
135 printf("fetchieee754: getbyte(0x%08x, %d) = 0x%02x\n", (unsigned int)(bufp)+offset[*fieldindex], *fieldindex, val);
136 #endif
137 (*fieldindex)++;
138 return val;
139 }
140
141 #ifdef __not_yet__
142 static void
put_byte(unsigned char * bufp,offsets_t offsets,int * fieldindex,unsigned char val)143 put_byte(
144 unsigned char *bufp,
145 offsets_t offsets,
146 int *fieldindex,
147 unsigned char val
148 )
149 {
150 *(bufp + offsets[*fieldindex]) = val;
151 (*fieldindex)++;
152 }
153 #endif
154
155 /*
156 * make conversions to and from external IEEE754 formats and internal
157 * NTP FP format.
158 */
159 int
fetch_ieee754(unsigned char ** buffpp,int size,l_fp * lfpp,offsets_t offsets)160 fetch_ieee754(
161 unsigned char **buffpp,
162 int size,
163 l_fp *lfpp,
164 offsets_t offsets
165 )
166 {
167 unsigned char *bufp = *buffpp;
168 unsigned int sign;
169 unsigned int bias;
170 unsigned int maxexp;
171 int mbits;
172 u_long mantissa_low;
173 u_long mantissa_high;
174 u_long characteristic;
175 long exponent;
176 #ifdef LIBDEBUG
177 int length;
178 #endif
179 unsigned char val;
180 int fieldindex = 0;
181
182 switch (size)
183 {
184 case IEEE_DOUBLE:
185 #ifdef LIBDEBUG
186 length = 8;
187 #endif
188 mbits = 52;
189 bias = 1023;
190 maxexp = 2047;
191 break;
192
193 case IEEE_SINGLE:
194 #ifdef LIBDEBUG
195 length = 4;
196 #endif
197 mbits = 23;
198 bias = 127;
199 maxexp = 255;
200 break;
201
202 default:
203 return IEEE_BADCALL;
204 }
205
206 val = get_byte(bufp, offsets, &fieldindex); /* fetch sign byte & first part of characteristic */
207
208 sign = (val & 0x80) != 0;
209 characteristic = (val & 0x7F);
210
211 val = get_byte(bufp, offsets, &fieldindex); /* fetch rest of characteristic and start of mantissa */
212
213 switch (size)
214 {
215 case IEEE_SINGLE:
216 characteristic <<= 1;
217 characteristic |= (val & 0x80) != 0; /* grab last characteristic bit */
218
219 mantissa_high = 0;
220
221 mantissa_low = (val &0x7F) << 16;
222 mantissa_low |= (u_long)get_byte(bufp, offsets, &fieldindex) << 8;
223 mantissa_low |= get_byte(bufp, offsets, &fieldindex);
224 break;
225
226 case IEEE_DOUBLE:
227 characteristic <<= 4;
228 characteristic |= (val & 0xF0) >> 4; /* grab lower characteristic bits */
229
230 mantissa_high = (val & 0x0F) << 16;
231 mantissa_high |= (u_long)get_byte(bufp, offsets, &fieldindex) << 8;
232 mantissa_high |= get_byte(bufp, offsets, &fieldindex);
233
234 mantissa_low = (u_long)get_byte(bufp, offsets, &fieldindex) << 24;
235 mantissa_low |= (u_long)get_byte(bufp, offsets, &fieldindex) << 16;
236 mantissa_low |= (u_long)get_byte(bufp, offsets, &fieldindex) << 8;
237 mantissa_low |= get_byte(bufp, offsets, &fieldindex);
238 break;
239
240 default:
241 return IEEE_BADCALL;
242 }
243 #ifdef LIBDEBUG
244 if (debug > 4)
245 {
246 double d;
247 float f;
248
249 if (size == IEEE_SINGLE)
250 {
251 int i;
252
253 for (i = 0; i < length; i++)
254 {
255 *((unsigned char *)(&f)+i) = *(*buffpp + offsets[i]);
256 }
257 d = f;
258 }
259 else
260 {
261 int i;
262
263 for (i = 0; i < length; i++)
264 {
265 *((unsigned char *)(&d)+i) = *(*buffpp + offsets[i]);
266 }
267 }
268
269 printf("fetchieee754: FP: %s -> %s -> %e(=%s)\n", fmt_hex(*buffpp, length),
270 fmt_flt(sign, mantissa_high, mantissa_low, characteristic),
271 d, fmt_hex((unsigned char *)&d, length));
272 }
273 #endif
274
275 *buffpp += fieldindex;
276
277 /*
278 * detect funny numbers
279 */
280 if (characteristic == maxexp)
281 {
282 /*
283 * NaN or Infinity
284 */
285 if (mantissa_low || mantissa_high)
286 {
287 /*
288 * NaN
289 */
290 return IEEE_NAN;
291 }
292 else
293 {
294 /*
295 * +Inf or -Inf
296 */
297 return sign ? IEEE_NEGINFINITY : IEEE_POSINFINITY;
298 }
299 }
300 else
301 {
302 /*
303 * collect real numbers
304 */
305
306 L_CLR(lfpp);
307
308 /*
309 * check for overflows
310 */
311 exponent = characteristic - bias;
312
313 if (exponent > 31) /* sorry - hardcoded */
314 {
315 /*
316 * overflow only in respect to NTP-FP representation
317 */
318 return sign ? IEEE_NEGOVERFLOW : IEEE_POSOVERFLOW;
319 }
320 else
321 {
322 int frac_offset; /* where the fraction starts */
323
324 frac_offset = mbits - exponent;
325
326 if (characteristic == 0)
327 {
328 /*
329 * de-normalized or tiny number - fits only as 0
330 */
331 return IEEE_OK;
332 }
333 else
334 {
335 /*
336 * adjust for implied 1
337 */
338 if (mbits > 31)
339 mantissa_high |= 1 << (mbits - 32);
340 else
341 mantissa_low |= 1 << mbits;
342
343 /*
344 * take mantissa apart - if only all machine would support
345 * 64 bit operations 8-(
346 */
347 if (frac_offset > mbits)
348 {
349 lfpp->l_ui = 0; /* only fractional number */
350 frac_offset -= mbits + 1; /* will now contain right shift count - 1*/
351 if (mbits > 31)
352 {
353 lfpp->l_uf = mantissa_high << (63 - mbits);
354 lfpp->l_uf |= mantissa_low >> (mbits - 33);
355 lfpp->l_uf >>= frac_offset;
356 }
357 else
358 {
359 lfpp->l_uf = mantissa_low >> frac_offset;
360 }
361 }
362 else
363 {
364 if (frac_offset > 32)
365 {
366 /*
367 * must split in high word
368 */
369 lfpp->l_ui = mantissa_high >> (frac_offset - 32);
370 lfpp->l_uf = (mantissa_high & ((1 << (frac_offset - 32)) - 1)) << (64 - frac_offset);
371 lfpp->l_uf |= mantissa_low >> (frac_offset - 32);
372 }
373 else
374 {
375 /*
376 * must split in low word
377 */
378 lfpp->l_ui = mantissa_high << (32 - frac_offset);
379 lfpp->l_ui |= (mantissa_low >> frac_offset) & ((1 << (32 - frac_offset)) - 1);
380 lfpp->l_uf = (mantissa_low & ((1 << frac_offset) - 1)) << (32 - frac_offset);
381 }
382 }
383
384 /*
385 * adjust for sign
386 */
387 if (sign)
388 {
389 L_NEG(lfpp);
390 }
391
392 return IEEE_OK;
393 }
394 }
395 }
396 }
397
398 /*
399 * DLH: This function is currently unused in ntpd. If you think about
400 * using it, be sure it does what you intend. I notice the bufpp arg
401 * is never referenced, and the calculated mantissa_high & mantissa_low
402 * are only referenced in debug output. It seems they're supposed to
403 * be composed into an ieee754-format float and stored at *bufpp or
404 * possibly **bufpp. Brought to my attention by this:
405 *
406 * ieee754io.c:414:10: warning: variable 'mantissa_low' set but not used
407 * [-Wunused-but-set-variable]
408 *
409 * To quiet it I'm #ifdef'ing the function away for now, here and below
410 * the call to it in main().
411 */
412 #ifdef PUT_IEEE754_UNUSED_FUNC
413 int
put_ieee754(unsigned char ** bufpp,int size,l_fp * lfpp,offsets_t offsets)414 put_ieee754(
415 unsigned char **bufpp,
416 int size,
417 l_fp *lfpp,
418 offsets_t offsets
419 )
420 {
421 l_fp outlfp;
422 #ifdef LIBDEBUG
423 unsigned int sign;
424 unsigned int bias;
425 #endif
426 /*unsigned int maxexp;*/
427 int mbits;
428 int msb;
429 u_long mantissa_low = 0;
430 u_long mantissa_high = 0;
431 #ifdef LIBDEBUG
432 u_long characteristic = 0;
433 long exponent;
434 #endif
435 /*int length;*/
436 unsigned long mask;
437
438 outlfp = *lfpp;
439
440 switch (size)
441 {
442 case IEEE_DOUBLE:
443 /*length = 8;*/
444 mbits = 52;
445 #ifdef LIBDEBUG
446 bias = 1023;
447 #endif
448 /*maxexp = 2047;*/
449 break;
450
451 case IEEE_SINGLE:
452 /*length = 4;*/
453 mbits = 23;
454 #ifdef LIBDEBUG
455 bias = 127;
456 #endif
457 /*maxexp = 255;*/
458 break;
459
460 default:
461 return IEEE_BADCALL;
462 }
463
464 /*
465 * find sign
466 */
467 if (L_ISNEG(&outlfp))
468 {
469 L_NEG(&outlfp);
470 #ifdef LIBDEBUG
471 sign = 1;
472 #endif
473 }
474 else
475 {
476 #ifdef LIBDEBUG
477 sign = 0;
478 #endif
479 }
480
481 if (L_ISZERO(&outlfp))
482 {
483 #ifdef LIBDEBUG
484 exponent = mantissa_high = mantissa_low = 0; /* true zero */
485 #endif
486 }
487 else
488 {
489 /*
490 * find number of significant integer bits
491 */
492 mask = 0x80000000;
493 if (outlfp.l_ui)
494 {
495 msb = 63;
496 while (mask && ((outlfp.l_ui & mask) == 0))
497 {
498 mask >>= 1;
499 msb--;
500 }
501 }
502 else
503 {
504 msb = 31;
505 while (mask && ((outlfp.l_uf & mask) == 0))
506 {
507 mask >>= 1;
508 msb--;
509 }
510 }
511
512 switch (size)
513 {
514 case IEEE_SINGLE:
515 mantissa_high = 0;
516 if (msb >= 32)
517 {
518 mantissa_low = (outlfp.l_ui & ((1 << (msb - 32)) - 1)) << (mbits - (msb - 32));
519 mantissa_low |= outlfp.l_uf >> (mbits - (msb - 32));
520 }
521 else
522 {
523 mantissa_low = (outlfp.l_uf << (mbits - msb)) & ((1 << mbits) - 1);
524 }
525 break;
526
527 case IEEE_DOUBLE:
528 if (msb >= 32)
529 {
530 mantissa_high = (outlfp.l_ui << (mbits - msb)) & ((1 << (mbits - 32)) - 1);
531 mantissa_high |= outlfp.l_uf >> (32 - (mbits - msb));
532 mantissa_low = (outlfp.l_ui & ((1 << (msb - mbits)) - 1)) << (32 - (msb - mbits));
533 mantissa_low |= outlfp.l_uf >> (msb - mbits);
534 }
535 else
536 {
537 mantissa_high = outlfp.l_uf << (mbits - 32 - msb);
538 mantissa_low = outlfp.l_uf << (mbits - 32);
539 }
540 }
541
542 #ifdef LIBDEBUG
543 exponent = msb - 32;
544 characteristic = exponent + bias;
545
546 if (debug > 4)
547 printf("FP: %s\n", fmt_flt(sign, mantissa_high, mantissa_low, characteristic));
548 #endif
549 }
550 return IEEE_OK;
551 }
552 #endif /* PUT_IEEE754_UNUSED_FUNC */
553
554
555 #if defined(DEBUG) && defined(LIBDEBUG)
main(int argc,char ** argv)556 int main(
557 int argc,
558 char **argv
559 )
560 {
561 static offsets_t native_off = { 0, 1, 2, 3, 4, 5, 6, 7 };
562 double f = 1.0;
563 double *f_p = &f;
564 l_fp fp;
565
566 if (argc == 2)
567 {
568 if (sscanf(argv[1], "%lf", &f) != 1)
569 {
570 printf("cannot convert %s to a float\n", argv[1]);
571 return 1;
572 }
573 }
574
575 printf("double: %s %s\n", fmt_blong(*(unsigned long *)&f, 32), fmt_blong(*(unsigned long *)((char *)(&f)+4), 32));
576 printf("fetch from %f = %d\n", f, fetch_ieee754((void *)&f_p, IEEE_DOUBLE, &fp, native_off));
577 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));
578 f_p = &f;
579 #ifdef PUT_IEEE754_UNUSED_FUNC
580 put_ieee754((void *)&f_p, IEEE_DOUBLE, &fp, native_off);
581 /* there should be a check on *f_p (f) having the expected result here */
582 #endif /* PUT_IEEE754_UNUSED_FUNC */
583
584 return 0;
585 }
586
587 #endif
588 /*
589 * History:
590 *
591 * ieee754io.c,v
592 * Revision 4.12 2005/04/16 17:32:10 kardel
593 * update copyright
594 *
595 * Revision 4.11 2004/11/14 15:29:41 kardel
596 * support PPSAPI, upgrade Copyright to Berkeley style
597 *
598 * Revision 4.8 1999/02/21 12:17:36 kardel
599 * 4.91f reconcilation
600 *
601 * Revision 4.7 1999/02/21 11:26:03 kardel
602 * renamed index to fieldindex to avoid index() name clash
603 *
604 * Revision 4.6 1998/11/15 20:27:52 kardel
605 * Release 4.0.73e13 reconcilation
606 *
607 * Revision 4.5 1998/08/16 19:01:51 kardel
608 * debug information only compile for LIBDEBUG case
609 *
610 * Revision 4.4 1998/08/09 09:39:28 kardel
611 * Release 4.0.73e2 reconcilation
612 *
613 * Revision 4.3 1998/06/13 11:56:19 kardel
614 * disabled putbute() for the time being
615 *
616 * Revision 4.2 1998/06/12 15:16:58 kardel
617 * ansi2knr compatibility
618 *
619 * Revision 4.1 1998/05/24 07:59:56 kardel
620 * conditional debug support
621 *
622 * Revision 4.0 1998/04/10 19:46:29 kardel
623 * Start 4.0 release version numbering
624 *
625 * Revision 1.1 1998/04/10 19:27:46 kardel
626 * initial NTP VERSION 4 integration of PARSE with GPS166 binary support
627 *
628 * Revision 1.1 1997/10/06 21:05:45 kardel
629 * new parse structure
630 *
631 */
632