dist/ntp/libparse/ieee754io.c

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