/*      $NetBSD: ieee754io.c,v 1.6 2024/08/18 20:47:17 christos Exp $   */

/*
* /src/NTP/ntp4-dev/libntp/ieee754io.c,v 4.12 2005/04/16 17:32:10 kardel RELEASE_20050508_A
*
* ieee754io.c,v 4.12 2005/04/16 17:32:10 kardel RELEASE_20050508_A
*
* $Created: Sun Jul 13 09:12:02 1997 $
*
* Copyright (c) 1997-2005 by Frank Kardel <kardel <AT> ntp.org>
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
*    notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
*    notice, this list of conditions and the following disclaimer in the
*    documentation and/or other materials provided with the distribution.
* 3. Neither the name of the author nor the names of its contributors
*    may be used to endorse or promote products derived from this software
*    without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*
*/

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <stdio.h>
#include "l_stdlib.h"
#include "ntp_stdlib.h"
#include "ntp_fp.h"
#include "ieee754io.h"

static unsigned char get_byte (unsigned char *, offsets_t, int *);
#ifdef __not_yet__
static void put_byte (unsigned char *, offsets_t, int *, unsigned char);
#endif

#ifdef LIBDEBUG

static char *
fmt_blong(
         unsigned long val,
         int cnt
         )
{
 char *buf, *s;
 int i = cnt;

 val <<= 32 - cnt;
 LIB_GETBUF(buf);
 s = buf;

 while (i--)
   {
     if (val & 0x80000000)
       {
         *s++ = '1';
       }
     else
       {
         *s++ = '0';
       }
     val <<= 1;
   }
 *s = '\0';
 return buf;
}

static char *
fmt_flt(
       unsigned int sign,
       unsigned long mh,
       unsigned long ml,
       unsigned long ch
       )
{
       char *buf;

       LIB_GETBUF(buf);
       snprintf(buf, LIB_BUFLENGTH, "%c %s %s %s", sign ? '-' : '+',
                fmt_blong(ch, 11),
                fmt_blong(mh, 20),
                fmt_blong(ml, 32));

       return buf;
}

static char *
fmt_hex(
       unsigned char *bufp,
       int length
       )
{
       char *  buf;
       char    hex[4];
       int     i;

       LIB_GETBUF(buf);
       buf[0] = '\0';
       for (i = 0; i < length; i++) {
               snprintf(hex, sizeof(hex), "%02x", bufp[i]);
               strlcat(buf, hex, LIB_BUFLENGTH);
       }

       return buf;
}

#endif

static unsigned char
get_byte(
        unsigned char *bufp,
        offsets_t offset,
        int *fieldindex
        )
{
 unsigned char val;

 val     = *(bufp + offset[*fieldindex]);
#ifdef LIBDEBUG
 if (debug > 4)
   printf("fetchieee754: getbyte(0x%08x, %d) = 0x%02x\n", (unsigned int)(bufp)+offset[*fieldindex], *fieldindex, val);
#endif
 (*fieldindex)++;
 return val;
}

#ifdef __not_yet__
static void
put_byte(
        unsigned char *bufp,
        offsets_t offsets,
        int *fieldindex,
        unsigned char val
        )
{
 *(bufp + offsets[*fieldindex]) = val;
 (*fieldindex)++;
}
#endif

/*
* make conversions to and from external IEEE754 formats and internal
* NTP FP format.
*/
int
fetch_ieee754(
             unsigned char **buffpp,
             int size,
             l_fp *lfpp,
             offsets_t offsets
             )
{
 unsigned char *bufp = *buffpp;
 unsigned int sign;
 unsigned int bias;
 unsigned int maxexp;
 int mbits;
 u_long mantissa_low;
 u_long mantissa_high;
 u_long characteristic;
 long exponent;
#ifdef LIBDEBUG
 int length;
#endif
 unsigned char val;
 int fieldindex = 0;

 switch (size)
   {
   case IEEE_DOUBLE:
#ifdef LIBDEBUG
     length = 8;
#endif
     mbits  = 52;
     bias   = 1023;
     maxexp = 2047;
     break;

   case IEEE_SINGLE:
#ifdef LIBDEBUG
     length = 4;
#endif
     mbits  = 23;
     bias   = 127;
     maxexp = 255;
     break;

   default:
     return IEEE_BADCALL;
   }

 val = get_byte(bufp, offsets, &fieldindex); /* fetch sign byte & first part of characteristic */

 sign     = (val & 0x80) != 0;
 characteristic = (val & 0x7F);

 val = get_byte(bufp, offsets, &fieldindex); /* fetch rest of characteristic and start of mantissa */

 switch (size)
   {
   case IEEE_SINGLE:
     characteristic <<= 1;
     characteristic  |= (val & 0x80) != 0; /* grab last characteristic bit */

     mantissa_high  = 0;

     mantissa_low   = (val &0x7F) << 16;
     mantissa_low  |= (u_long)get_byte(bufp, offsets, &fieldindex) << 8;
     mantissa_low  |= get_byte(bufp, offsets, &fieldindex);
     break;

   case IEEE_DOUBLE:
     characteristic <<= 4;
     characteristic  |= (val & 0xF0) >> 4; /* grab lower characteristic bits */

     mantissa_high  = (val & 0x0F) << 16;
     mantissa_high |= (u_long)get_byte(bufp, offsets, &fieldindex) << 8;
     mantissa_high |= get_byte(bufp, offsets, &fieldindex);

     mantissa_low   = (u_long)get_byte(bufp, offsets, &fieldindex) << 24;
     mantissa_low  |= (u_long)get_byte(bufp, offsets, &fieldindex) << 16;
     mantissa_low  |= (u_long)get_byte(bufp, offsets, &fieldindex) << 8;
     mantissa_low  |= get_byte(bufp, offsets, &fieldindex);
     break;

   default:
     return IEEE_BADCALL;
   }
#ifdef LIBDEBUG
 if (debug > 4)
 {
   double d;
   float f;

   if (size == IEEE_SINGLE)
     {
       int i;

       for (i = 0; i < length; i++)
         {
           *((unsigned char *)(&f)+i) = *(*buffpp + offsets[i]);
         }
       d = f;
     }
   else
     {
       int i;

       for (i = 0; i < length; i++)
         {
           *((unsigned char *)(&d)+i) = *(*buffpp + offsets[i]);
         }
     }

   printf("fetchieee754: FP: %s -> %s -> %e(=%s)\n", fmt_hex(*buffpp, length),
          fmt_flt(sign, mantissa_high, mantissa_low, characteristic),
          d, fmt_hex((unsigned char *)&d, length));
 }
#endif

 *buffpp += fieldindex;

 /*
  * detect funny numbers
  */
 if (characteristic == maxexp)
   {
     /*
      * NaN or Infinity
      */
     if (mantissa_low || mantissa_high)
       {
         /*
          * NaN
          */
         return IEEE_NAN;
       }
     else
       {
         /*
          * +Inf or -Inf
          */
         return sign ? IEEE_NEGINFINITY : IEEE_POSINFINITY;
       }
   }
 else
   {
     /*
      * collect real numbers
      */

     L_CLR(lfpp);

     /*
      * check for overflows
      */
     exponent = characteristic - bias;

     if (exponent > 31)        /* sorry - hardcoded */
       {
         /*
          * overflow only in respect to NTP-FP representation
          */
         return sign ? IEEE_NEGOVERFLOW : IEEE_POSOVERFLOW;
       }
     else
       {
         int frac_offset;      /* where the fraction starts */

         frac_offset = mbits - exponent;

         if (characteristic == 0)
           {
             /*
              * de-normalized or tiny number - fits only as 0
              */
             return IEEE_OK;
           }
         else
           {
             /*
              * adjust for implied 1
              */
             if (mbits > 31)
               mantissa_high |= 1 << (mbits - 32);
             else
               mantissa_low  |= 1 << mbits;

             /*
              * take mantissa apart - if only all machine would support
              * 64 bit operations 8-(
              */
             if (frac_offset > mbits)
               {
                 lfpp->l_ui = 0; /* only fractional number */
                 frac_offset -= mbits + 1; /* will now contain right shift count - 1*/
                 if (mbits > 31)
                   {
                     lfpp->l_uf   = mantissa_high << (63 - mbits);
                     lfpp->l_uf  |= mantissa_low  >> (mbits - 33);
                     lfpp->l_uf >>= frac_offset;
                   }
                 else
                   {
                     lfpp->l_uf = mantissa_low >> frac_offset;
                   }
               }
             else
               {
                 if (frac_offset > 32)
                   {
                     /*
                      * must split in high word
                      */
                     lfpp->l_ui  =  mantissa_high >> (frac_offset - 32);
                     lfpp->l_uf  = (mantissa_high & ((1 << (frac_offset - 32)) - 1)) << (64 - frac_offset);
                     lfpp->l_uf |=  mantissa_low  >> (frac_offset - 32);
                   }
                 else
                   {
                     /*
                      * must split in low word
                      */
                     lfpp->l_ui  =  mantissa_high << (32 - frac_offset);
                     lfpp->l_ui |= (mantissa_low >> frac_offset) & ((1 << (32 - frac_offset)) - 1);
                     lfpp->l_uf  = (mantissa_low & ((1 << frac_offset) - 1)) << (32 - frac_offset);
                   }
               }

             /*
              * adjust for sign
              */
             if (sign)
               {
                 L_NEG(lfpp);
               }

             return IEEE_OK;
           }
       }
   }
}

/*
* DLH: This function is currently unused in ntpd.  If you think about
* using it, be sure it does what you intend.  I notice the bufpp arg
* is never referenced, and the calculated mantissa_high & mantissa_low
* are only referenced in debug output.  It seems they're supposed to
* be composed into an ieee754-format float and stored at *bufpp or
* possibly **bufpp.  Brought to my attention by this:
*
* ieee754io.c:414:10: warning: variable 'mantissa_low' set but not used
* [-Wunused-but-set-variable]
*
* To quiet it I'm #ifdef'ing the function away for now, here and below
* the call to it in main().
*/
#ifdef PUT_IEEE754_UNUSED_FUNC
int
put_ieee754(
           unsigned char **bufpp,
           int size,
           l_fp *lfpp,
           offsets_t offsets
           )
{
 l_fp outlfp;
#ifdef LIBDEBUG
 unsigned int sign;
 unsigned int bias;
#endif
/*unsigned int maxexp;*/
 int mbits;
 int msb;
 u_long mantissa_low = 0;
 u_long mantissa_high = 0;
#ifdef LIBDEBUG
 u_long characteristic = 0;
 long exponent;
#endif
/*int length;*/
 unsigned long mask;

 outlfp = *lfpp;

 switch (size)
   {
   case IEEE_DOUBLE:
   /*length = 8;*/
     mbits  = 52;
#ifdef LIBDEBUG
     bias   = 1023;
#endif
   /*maxexp = 2047;*/
     break;

   case IEEE_SINGLE:
   /*length = 4;*/
     mbits  = 23;
#ifdef LIBDEBUG
     bias   = 127;
#endif
   /*maxexp = 255;*/
     break;

   default:
     return IEEE_BADCALL;
   }

 /*
  * find sign
  */
 if (L_ISNEG(&outlfp))
   {
     L_NEG(&outlfp);
#ifdef LIBDEBUG
     sign = 1;
#endif
   }
 else
   {
#ifdef LIBDEBUG
     sign = 0;
#endif
   }

 if (L_ISZERO(&outlfp))
   {
#ifdef LIBDEBUG
     exponent = mantissa_high = mantissa_low = 0; /* true zero */
#endif
   }
 else
   {
     /*
      * find number of significant integer bits
      */
     mask = 0x80000000;
     if (outlfp.l_ui)
       {
         msb = 63;
         while (mask && ((outlfp.l_ui & mask) == 0))
           {
             mask >>= 1;
             msb--;
           }
       }
     else
       {
         msb = 31;
         while (mask && ((outlfp.l_uf & mask) == 0))
           {
             mask >>= 1;
             msb--;
           }
       }

     switch (size)
       {
       case IEEE_SINGLE:
         mantissa_high = 0;
         if (msb >= 32)
           {
             mantissa_low  = (outlfp.l_ui & ((1 << (msb - 32)) - 1)) << (mbits - (msb - 32));
             mantissa_low |=  outlfp.l_uf >> (mbits - (msb - 32));
           }
         else
           {
             mantissa_low  = (outlfp.l_uf << (mbits - msb)) & ((1 << mbits) - 1);
           }
         break;

       case IEEE_DOUBLE:
         if (msb >= 32)
           {
             mantissa_high  = (outlfp.l_ui << (mbits - msb)) & ((1 << (mbits - 32)) - 1);
             mantissa_high |=  outlfp.l_uf >> (32 - (mbits - msb));
             mantissa_low   = (outlfp.l_ui & ((1 << (msb - mbits)) - 1)) << (32 - (msb - mbits));
             mantissa_low  |=  outlfp.l_uf >> (msb - mbits);
           }
         else
           {
             mantissa_high  = outlfp.l_uf << (mbits - 32 - msb);
             mantissa_low   = outlfp.l_uf << (mbits - 32);
           }
       }

#ifdef LIBDEBUG
     exponent = msb - 32;
     characteristic = exponent + bias;

     if (debug > 4)
       printf("FP: %s\n", fmt_flt(sign, mantissa_high, mantissa_low, characteristic));
#endif
   }
 return IEEE_OK;
}
#endif  /* PUT_IEEE754_UNUSED_FUNC */


#if defined(DEBUG) && defined(LIBDEBUG)
int main(
        int argc,
        char **argv
        )
{
 static offsets_t native_off = { 0, 1, 2, 3, 4, 5, 6, 7 };
 double f = 1.0;
 double *f_p = &f;
 l_fp fp;

 if (argc == 2)
   {
     if (sscanf(argv[1], "%lf", &f) != 1)
       {
         printf("cannot convert %s to a float\n", argv[1]);
         return 1;
       }
   }

 printf("double: %s %s\n", fmt_blong(*(unsigned long *)&f, 32), fmt_blong(*(unsigned long *)((char *)(&f)+4), 32));
 printf("fetch from %f = %d\n", f, fetch_ieee754((void *)&f_p, IEEE_DOUBLE, &fp, native_off));
 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));
 f_p = &f;
#ifdef PUT_IEEE754_UNUSED_FUNC
 put_ieee754((void *)&f_p, IEEE_DOUBLE, &fp, native_off);
/* there should be a check on *f_p (f) having the expected result here */
#endif  /* PUT_IEEE754_UNUSED_FUNC */

 return 0;
}

#endif
/*
* History:
*
* ieee754io.c,v
* Revision 4.12  2005/04/16 17:32:10  kardel
* update copyright
*
* Revision 4.11  2004/11/14 15:29:41  kardel
* support PPSAPI, upgrade Copyright to Berkeley style
*
* Revision 4.8  1999/02/21 12:17:36  kardel
* 4.91f reconcilation
*
* Revision 4.7  1999/02/21 11:26:03  kardel
* renamed index to fieldindex to avoid index() name clash
*
* Revision 4.6  1998/11/15 20:27:52  kardel
* Release 4.0.73e13 reconcilation
*
* Revision 4.5  1998/08/16 19:01:51  kardel
* debug information only compile for LIBDEBUG case
*
* Revision 4.4  1998/08/09 09:39:28  kardel
* Release 4.0.73e2 reconcilation
*
* Revision 4.3  1998/06/13 11:56:19  kardel
* disabled putbute() for the time being
*
* Revision 4.2  1998/06/12 15:16:58  kardel
* ansi2knr compatibility
*
* Revision 4.1  1998/05/24 07:59:56  kardel
* conditional debug support
*
* Revision 4.0  1998/04/10 19:46:29  kardel
* Start 4.0 release version numbering
*
* Revision 1.1  1998/04/10 19:27:46  kardel
* initial NTP VERSION 4 integration of PARSE with GPS166 binary support
*
* Revision 1.1  1997/10/06 21:05:45  kardel
* new parse structure
*
*/