contrib/ntp/libntp/ieee754io.c

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