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 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 *
  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   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 *
 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
 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
 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
 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
 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)
 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  */