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
46static unsigned char get_byte P((unsigned char *, offsets_t, int *));
47#ifdef __not_yet__
48static void put_byte P((unsigned char *, offsets_t, int *, unsigned char));
49#endif
50
51#ifdef LIBDEBUG
52
53#include "lib_strbuf.h"
54
55static char *
56fmt_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
84static char *
85fmt_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
102static char *
103fmt_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
121static unsigned char
122get_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__
140static void
141put_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 */
157int
158fetch_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
396int
397put_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)
538int 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 */
611