/contrib/ntp/libparse/ieee754io.c
C | 610 lines | 415 code | 61 blank | 134 comment | 44 complexity | 91e0acfe0e3487a6c6350f2c1a39c33c MD5 | raw file
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 46static unsigned char get_byte P((unsigned char *, offsets_t, int *)); 47#ifdef __not_yet__ 48static void put_byte P((unsigned char *, offsets_t, int *, unsigned char)); 49#endif 50 51#ifdef LIBDEBUG 52 53#include "lib_strbuf.h" 54 55static char * 56fmt_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 84static char * 85fmt_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 102static char * 103fmt_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 121static unsigned char 122get_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__ 140static void 141put_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 */ 157int 158fetch_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 396int 397put_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) 538int 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 */