contrib/ntp/libparse/ieee754io.c

   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 #include "lib_strbuf.h"
  54
  55 static char *
  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 *
  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 *
 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
 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
 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
 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 int
 401 put_ieee754(
 402             unsigned char **bufpp,
 403             int size,
 404             l_fp *lfpp,
 405             offsets_t offsets
 406             )
 407 {
 408   l_fp outlfp;
 409 #ifdef LIBDEBUG
 410   unsigned int sign;
 411   unsigned int bias;
 412 #endif
 413 /*unsigned int maxexp;*/
 414   int mbits;
 415   int msb;
 416   u_long mantissa_low = 0;
 417   u_long mantissa_high = 0;
 418 #ifdef LIBDEBUG
 419   u_long characteristic = 0;
 420   long exponent;
 421 #endif
 422 /*int length;*/
 423   unsigned long mask;
 424
 425   outlfp = *lfpp;
 426
 427   switch (size)
 428     {
 429     case IEEE_DOUBLE:
 430     /*length = 8;*/
 431       mbits  = 52;
 432 #ifdef LIBDEBUG
 433       bias   = 1023;
 434 #endif
 435     /*maxexp = 2047;*/
 436       break;
 437
 438     case IEEE_SINGLE:
 439     /*length = 4;*/
 440       mbits  = 23;
 441 #ifdef LIBDEBUG
 442       bias   = 127;
 443 #endif
 444     /*maxexp = 255;*/
 445       break;
 446
 447     default:
 448       return IEEE_BADCALL;
 449     }
 450
 451   /*
 452    * find sign
 453    */
 454   if (L_ISNEG(&outlfp))
 455     {
 456       L_NEG(&outlfp);
 457 #ifdef LIBDEBUG
 458       sign = 1;
 459 #endif
 460     }
 461   else
 462     {
 463 #ifdef LIBDEBUG
 464       sign = 0;
 465 #endif
 466     }
 467
 468   if (L_ISZERO(&outlfp))
 469     {
 470 #ifdef LIBDEBUG
 471       exponent = mantissa_high = mantissa_low = 0; /* true zero */
 472 #endif
 473     }
 474   else
 475     {
 476       /*
 477        * find number of significant integer bits
 478        */
 479       mask = 0x80000000;
 480       if (outlfp.l_ui)
 481         {
 482           msb = 63;
 483           while (mask && ((outlfp.l_ui & mask) == 0))
 484             {
 485               mask >>= 1;
 486               msb--;
 487             }
 488         }
 489       else
 490         {
 491           msb = 31;
 492           while (mask && ((outlfp.l_uf & mask) == 0))
 493             {
 494               mask >>= 1;
 495               msb--;
 496             }
 497         }
 498
 499       switch (size)
 500         {
 501         case IEEE_SINGLE:
 502           mantissa_high = 0;
 503           if (msb >= 32)
 504             {
 505               mantissa_low  = (outlfp.l_ui & ((1 << (msb - 32)) - 1)) << (mbits - (msb - 32));
 506               mantissa_low |=  outlfp.l_uf >> (mbits - (msb - 32));
 507             }
 508           else
 509             {
 510               mantissa_low  = (outlfp.l_uf << (mbits - msb)) & ((1 << mbits) - 1);
 511             }
 512           break;
 513
 514         case IEEE_DOUBLE:
 515           if (msb >= 32)
 516             {
 517               mantissa_high  = (outlfp.l_ui << (mbits - msb)) & ((1 << (mbits - 32)) - 1);
 518               mantissa_high |=  outlfp.l_uf >> (32 - (mbits - msb));
 519               mantissa_low   = (outlfp.l_ui & ((1 << (msb - mbits)) - 1)) << (32 - (msb - mbits));
 520               mantissa_low  |=  outlfp.l_uf >> (msb - mbits);
 521             }
 522           else
 523             {
 524               mantissa_high  = outlfp.l_uf << (mbits - 32 - msb);
 525               mantissa_low   = outlfp.l_uf << (mbits - 32);
 526             }
 527         }
 528
 529 #ifdef LIBDEBUG
 530       exponent = msb - 32;
 531       characteristic = exponent + bias;
 532
 533       if (debug > 4)
 534         printf("FP: %s\n", fmt_flt(sign, mantissa_high, mantissa_low, characteristic));
 535 #endif
 536     }
 537   return IEEE_OK;
 538 }
 539
 540
 541 #if defined(DEBUG) && defined(LIBDEBUG)
 542 int main(
 543          int argc,
 544          char **argv
 545          )
 546 {
 547   static offsets_t native_off = { 0, 1, 2, 3, 4, 5, 6, 7 };
 548   double f = 1.0;
 549   double *f_p = &f;
 550   l_fp fp;
 551
 552   if (argc == 2)
 553     {
 554       if (sscanf(argv[1], "%lf", &f) != 1)
 555         {
 556           printf("cannot convert %s to a float\n", argv[1]);
 557           return 1;
 558         }
 559     }
 560
 561   printf("double: %s %s\n", fmt_blong(*(unsigned long *)&f, 32), fmt_blong(*(unsigned long *)((char *)(&f)+4), 32));
 562   printf("fetch from %f = %d\n", f, fetch_ieee754((void *)&f_p, IEEE_DOUBLE, &fp, native_off));
 563   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));
 564   f_p = &f;
 565   put_ieee754((void *)&f_p, IEEE_DOUBLE, &fp, native_off);
 566
 567   return 0;
 568 }
 569
 570 #endif
 571 /*
 572  * History:
 573  *
 574  * ieee754io.c,v
 575  * Revision 4.12  2005/04/16 17:32:10  kardel
 576  * update copyright
 577  *
 578  * Revision 4.11  2004/11/14 15:29:41  kardel
 579  * support PPSAPI, upgrade Copyright to Berkeley style
 580  *
 581  * Revision 4.8  1999/02/21 12:17:36  kardel
 582  * 4.91f reconcilation
 583  *
 584  * Revision 4.7  1999/02/21 11:26:03  kardel
 585  * renamed index to fieldindex to avoid index() name clash
 586  *
 587  * Revision 4.6  1998/11/15 20:27:52  kardel
 588  * Release 4.0.73e13 reconcilation
 589  *
 590  * Revision 4.5  1998/08/16 19:01:51  kardel
 591  * debug information only compile for LIBDEBUG case
 592  *
 593  * Revision 4.4  1998/08/09 09:39:28  kardel
 594  * Release 4.0.73e2 reconcilation
 595  *
 596  * Revision 4.3  1998/06/13 11:56:19  kardel
 597  * disabled putbute() for the time being
 598  *
 599  * Revision 4.2  1998/06/12 15:16:58  kardel
 600  * ansi2knr compatibility
 601  *
 602  * Revision 4.1  1998/05/24 07:59:56  kardel
 603  * conditional debug support
 604  *
 605  * Revision 4.0  1998/04/10 19:46:29  kardel
 606  * Start 4.0 release version numbering
 607  *
 608  * Revision 1.1  1998/04/10 19:27:46  kardel
 609  * initial NTP VERSION 4 integration of PARSE with GPS166 binary support
 610  *
 611  * Revision 1.1  1997/10/06 21:05:45  kardel
 612  * new parse structure
 613  *
 614  */