Home | History | Annotate | Line # | Download | only in libparse
      1 /*	$NetBSD: ieee754io.c,v 1.6 2024/08/18 20:47:17 christos Exp $	*/
      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 (unsigned char *, offsets_t, int *);
     49 #ifdef __not_yet__
     50 static void put_byte (unsigned char *, offsets_t, int *, unsigned char);
     51 #endif
     52 
     53 #ifdef LIBDEBUG
     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 /*
    401  * DLH: This function is currently unused in ntpd.  If you think about
    402  * using it, be sure it does what you intend.  I notice the bufpp arg
    403  * is never referenced, and the calculated mantissa_high & mantissa_low
    404  * are only referenced in debug output.  It seems they're supposed to
    405  * be composed into an ieee754-format float and stored at *bufpp or
    406  * possibly **bufpp.  Brought to my attention by this:
    407  *
    408  * ieee754io.c:414:10: warning: variable 'mantissa_low' set but not used
    409  * [-Wunused-but-set-variable]
    410  *
    411  * To quiet it I'm #ifdef'ing the function away for now, here and below
    412  * the call to it in main().
    413  */
    414 #ifdef PUT_IEEE754_UNUSED_FUNC
    415 int
    416 put_ieee754(
    417 	    unsigned char **bufpp,
    418 	    int size,
    419 	    l_fp *lfpp,
    420 	    offsets_t offsets
    421 	    )
    422 {
    423   l_fp outlfp;
    424 #ifdef LIBDEBUG
    425   unsigned int sign;
    426   unsigned int bias;
    427 #endif
    428 /*unsigned int maxexp;*/
    429   int mbits;
    430   int msb;
    431   u_long mantissa_low = 0;
    432   u_long mantissa_high = 0;
    433 #ifdef LIBDEBUG
    434   u_long characteristic = 0;
    435   long exponent;
    436 #endif
    437 /*int length;*/
    438   unsigned long mask;
    439 
    440   outlfp = *lfpp;
    441 
    442   switch (size)
    443     {
    444     case IEEE_DOUBLE:
    445     /*length = 8;*/
    446       mbits  = 52;
    447 #ifdef LIBDEBUG
    448       bias   = 1023;
    449 #endif
    450     /*maxexp = 2047;*/
    451       break;
    452 
    453     case IEEE_SINGLE:
    454     /*length = 4;*/
    455       mbits  = 23;
    456 #ifdef LIBDEBUG
    457       bias   = 127;
    458 #endif
    459     /*maxexp = 255;*/
    460       break;
    461 
    462     default:
    463       return IEEE_BADCALL;
    464     }
    465 
    466   /*
    467    * find sign
    468    */
    469   if (L_ISNEG(&outlfp))
    470     {
    471       L_NEG(&outlfp);
    472 #ifdef LIBDEBUG
    473       sign = 1;
    474 #endif
    475     }
    476   else
    477     {
    478 #ifdef LIBDEBUG
    479       sign = 0;
    480 #endif
    481     }
    482 
    483   if (L_ISZERO(&outlfp))
    484     {
    485 #ifdef LIBDEBUG
    486       exponent = mantissa_high = mantissa_low = 0; /* true zero */
    487 #endif
    488     }
    489   else
    490     {
    491       /*
    492        * find number of significant integer bits
    493        */
    494       mask = 0x80000000;
    495       if (outlfp.l_ui)
    496 	{
    497 	  msb = 63;
    498 	  while (mask && ((outlfp.l_ui & mask) == 0))
    499 	    {
    500 	      mask >>= 1;
    501 	      msb--;
    502 	    }
    503 	}
    504       else
    505 	{
    506 	  msb = 31;
    507 	  while (mask && ((outlfp.l_uf & mask) == 0))
    508 	    {
    509 	      mask >>= 1;
    510 	      msb--;
    511 	    }
    512 	}
    513 
    514       switch (size)
    515 	{
    516 	case IEEE_SINGLE:
    517 	  mantissa_high = 0;
    518 	  if (msb >= 32)
    519 	    {
    520 	      mantissa_low  = (outlfp.l_ui & ((1 << (msb - 32)) - 1)) << (mbits - (msb - 32));
    521 	      mantissa_low |=  outlfp.l_uf >> (mbits - (msb - 32));
    522 	    }
    523 	  else
    524 	    {
    525 	      mantissa_low  = (outlfp.l_uf << (mbits - msb)) & ((1 << mbits) - 1);
    526 	    }
    527 	  break;
    528 
    529 	case IEEE_DOUBLE:
    530 	  if (msb >= 32)
    531 	    {
    532 	      mantissa_high  = (outlfp.l_ui << (mbits - msb)) & ((1 << (mbits - 32)) - 1);
    533 	      mantissa_high |=  outlfp.l_uf >> (32 - (mbits - msb));
    534 	      mantissa_low   = (outlfp.l_ui & ((1 << (msb - mbits)) - 1)) << (32 - (msb - mbits));
    535 	      mantissa_low  |=  outlfp.l_uf >> (msb - mbits);
    536 	    }
    537 	  else
    538 	    {
    539 	      mantissa_high  = outlfp.l_uf << (mbits - 32 - msb);
    540 	      mantissa_low   = outlfp.l_uf << (mbits - 32);
    541 	    }
    542 	}
    543 
    544 #ifdef LIBDEBUG
    545       exponent = msb - 32;
    546       characteristic = exponent + bias;
    547 
    548       if (debug > 4)
    549 	printf("FP: %s\n", fmt_flt(sign, mantissa_high, mantissa_low, characteristic));
    550 #endif
    551     }
    552   return IEEE_OK;
    553 }
    554 #endif	/* PUT_IEEE754_UNUSED_FUNC */
    555 
    556 
    557 #if defined(DEBUG) && defined(LIBDEBUG)
    558 int main(
    559 	 int argc,
    560 	 char **argv
    561 	 )
    562 {
    563   static offsets_t native_off = { 0, 1, 2, 3, 4, 5, 6, 7 };
    564   double f = 1.0;
    565   double *f_p = &f;
    566   l_fp fp;
    567 
    568   if (argc == 2)
    569     {
    570       if (sscanf(argv[1], "%lf", &f) != 1)
    571 	{
    572 	  printf("cannot convert %s to a float\n", argv[1]);
    573 	  return 1;
    574 	}
    575     }
    576 
    577   printf("double: %s %s\n", fmt_blong(*(unsigned long *)&f, 32), fmt_blong(*(unsigned long *)((char *)(&f)+4), 32));
    578   printf("fetch from %f = %d\n", f, fetch_ieee754((void *)&f_p, IEEE_DOUBLE, &fp, native_off));
    579   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));
    580   f_p = &f;
    581 #ifdef PUT_IEEE754_UNUSED_FUNC
    582   put_ieee754((void *)&f_p, IEEE_DOUBLE, &fp, native_off);
    583 /* there should be a check on *f_p (f) having the expected result here */
    584 #endif	/* PUT_IEEE754_UNUSED_FUNC */
    585 
    586   return 0;
    587 }
    588 
    589 #endif
    590 /*
    591  * History:
    592  *
    593  * ieee754io.c,v
    594  * Revision 4.12  2005/04/16 17:32:10  kardel
    595  * update copyright
    596  *
    597  * Revision 4.11  2004/11/14 15:29:41  kardel
    598  * support PPSAPI, upgrade Copyright to Berkeley style
    599  *
    600  * Revision 4.8  1999/02/21 12:17:36  kardel
    601  * 4.91f reconcilation
    602  *
    603  * Revision 4.7  1999/02/21 11:26:03  kardel
    604  * renamed index to fieldindex to avoid index() name clash
    605  *
    606  * Revision 4.6  1998/11/15 20:27:52  kardel
    607  * Release 4.0.73e13 reconcilation
    608  *
    609  * Revision 4.5  1998/08/16 19:01:51  kardel
    610  * debug information only compile for LIBDEBUG case
    611  *
    612  * Revision 4.4  1998/08/09 09:39:28  kardel
    613  * Release 4.0.73e2 reconcilation
    614  *
    615  * Revision 4.3  1998/06/13 11:56:19  kardel
    616  * disabled putbute() for the time being
    617  *
    618  * Revision 4.2  1998/06/12 15:16:58  kardel
    619  * ansi2knr compatibility
    620  *
    621  * Revision 4.1  1998/05/24 07:59:56  kardel
    622  * conditional debug support
    623  *
    624  * Revision 4.0  1998/04/10 19:46:29  kardel
    625  * Start 4.0 release version numbering
    626  *
    627  * Revision 1.1  1998/04/10 19:27:46  kardel
    628  * initial NTP VERSION 4 integration of PARSE with GPS166 binary support
    629  *
    630  * Revision 1.1  1997/10/06 21:05:45  kardel
    631  * new parse structure
    632  *
    633  */
    634