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 P((unsigned char *, offsets_t, int *));
47 #ifdef __not_yet__
48 static void put_byte P((unsigned char *, offsets_t, int *, unsigned char));
49 #endif
50
51 #ifdef LIBDEBUG
52
53 #include "lib_strbuf.h"
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 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 *
fmt_hex(unsigned char * bufp,int length)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
get_byte(unsigned char * bufp,offsets_t offset,int * fieldindex)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
put_byte(unsigned char * bufp,offsets_t offsets,int * fieldindex,unsigned char val)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
fetch_ieee754(unsigned char ** buffpp,int size,l_fp * lfpp,offsets_t offsets)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
put_ieee754(unsigned char ** bufpp,int size,l_fp * lfpp,offsets_t offsets)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)
main(int argc,char ** argv)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 */
611