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