ieee754io.c revision 290001 
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 (unsigned char *, offsets_t, int *);
47#ifdef __not_yet__
48static void put_byte (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	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
103static char *
104fmt_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
125static unsigned char
126get_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__
144static void
145put_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 */
161int
162fetch_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
400int
401put_ieee754(
402	    unsigned char **bufpp,
403	    int size,
404	    l_fp *lfpp,
405	    offsets_t offsets
406	    )
407{
408  l_fp outlfp;
409#ifdef LIBDEBUG
410  unsigned int sign;
411  unsigned int bias;
412#endif
413/*unsigned int maxexp;*/
414  int mbits;
415  int msb;
416  u_long mantissa_low = 0;
417  u_long mantissa_high = 0;
418#ifdef LIBDEBUG
419  u_long characteristic = 0;
420  long exponent;
421#endif
422/*int length;*/
423  unsigned long mask;
424
425  outlfp = *lfpp;
426
427  switch (size)
428    {
429    case IEEE_DOUBLE:
430    /*length = 8;*/
431      mbits  = 52;
432#ifdef LIBDEBUG
433      bias   = 1023;
434#endif
435    /*maxexp = 2047;*/
436      break;
437
438    case IEEE_SINGLE:
439    /*length = 4;*/
440      mbits  = 23;
441#ifdef LIBDEBUG
442      bias   = 127;
443#endif
444    /*maxexp = 255;*/
445      break;
446
447    default:
448      return IEEE_BADCALL;
449    }
450
451  /*
452   * find sign
453   */
454  if (L_ISNEG(&outlfp))
455    {
456      L_NEG(&outlfp);
457#ifdef LIBDEBUG
458      sign = 1;
459#endif
460    }
461  else
462    {
463#ifdef LIBDEBUG
464      sign = 0;
465#endif
466    }
467
468  if (L_ISZERO(&outlfp))
469    {
470#ifdef LIBDEBUG
471      exponent = mantissa_high = mantissa_low = 0; /* true zero */
472#endif
473    }
474  else
475    {
476      /*
477       * find number of significant integer bits
478       */
479      mask = 0x80000000;
480      if (outlfp.l_ui)
481	{
482	  msb = 63;
483	  while (mask && ((outlfp.l_ui & mask) == 0))
484	    {
485	      mask >>= 1;
486	      msb--;
487	    }
488	}
489      else
490	{
491	  msb = 31;
492	  while (mask && ((outlfp.l_uf & mask) == 0))
493	    {
494	      mask >>= 1;
495	      msb--;
496	    }
497	}
498
499      switch (size)
500	{
501	case IEEE_SINGLE:
502	  mantissa_high = 0;
503	  if (msb >= 32)
504	    {
505	      mantissa_low  = (outlfp.l_ui & ((1 << (msb - 32)) - 1)) << (mbits - (msb - 32));
506	      mantissa_low |=  outlfp.l_uf >> (mbits - (msb - 32));
507	    }
508	  else
509	    {
510	      mantissa_low  = (outlfp.l_uf << (mbits - msb)) & ((1 << mbits) - 1);
511	    }
512	  break;
513
514	case IEEE_DOUBLE:
515	  if (msb >= 32)
516	    {
517	      mantissa_high  = (outlfp.l_ui << (mbits - msb)) & ((1 << (mbits - 32)) - 1);
518	      mantissa_high |=  outlfp.l_uf >> (32 - (mbits - msb));
519	      mantissa_low   = (outlfp.l_ui & ((1 << (msb - mbits)) - 1)) << (32 - (msb - mbits));
520	      mantissa_low  |=  outlfp.l_uf >> (msb - mbits);
521	    }
522	  else
523	    {
524	      mantissa_high  = outlfp.l_uf << (mbits - 32 - msb);
525	      mantissa_low   = outlfp.l_uf << (mbits - 32);
526	    }
527	}
528
529#ifdef LIBDEBUG
530      exponent = msb - 32;
531      characteristic = exponent + bias;
532
533      if (debug > 4)
534	printf("FP: %s\n", fmt_flt(sign, mantissa_high, mantissa_low, characteristic));
535#endif
536    }
537  return IEEE_OK;
538}
539
540
541#if defined(DEBUG) && defined(LIBDEBUG)
542int main(
543	 int argc,
544	 char **argv
545	 )
546{
547  static offsets_t native_off = { 0, 1, 2, 3, 4, 5, 6, 7 };
548  double f = 1.0;
549  double *f_p = &f;
550  l_fp fp;
551
552  if (argc == 2)
553    {
554      if (sscanf(argv[1], "%lf", &f) != 1)
555	{
556	  printf("cannot convert %s to a float\n", argv[1]);
557	  return 1;
558	}
559    }
560
561  printf("double: %s %s\n", fmt_blong(*(unsigned long *)&f, 32), fmt_blong(*(unsigned long *)((char *)(&f)+4), 32));
562  printf("fetch from %f = %d\n", f, fetch_ieee754((void *)&f_p, IEEE_DOUBLE, &fp, native_off));
563  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));
564  f_p = &f;
565  put_ieee754((void *)&f_p, IEEE_DOUBLE, &fp, native_off);
566
567  return 0;
568}
569
570#endif
571/*
572 * History:
573 *
574 * ieee754io.c,v
575 * Revision 4.12  2005/04/16 17:32:10  kardel
576 * update copyright
577 *
578 * Revision 4.11  2004/11/14 15:29:41  kardel
579 * support PPSAPI, upgrade Copyright to Berkeley style
580 *
581 * Revision 4.8  1999/02/21 12:17:36  kardel
582 * 4.91f reconcilation
583 *
584 * Revision 4.7  1999/02/21 11:26:03  kardel
585 * renamed index to fieldindex to avoid index() name clash
586 *
587 * Revision 4.6  1998/11/15 20:27:52  kardel
588 * Release 4.0.73e13 reconcilation
589 *
590 * Revision 4.5  1998/08/16 19:01:51  kardel
591 * debug information only compile for LIBDEBUG case
592 *
593 * Revision 4.4  1998/08/09 09:39:28  kardel
594 * Release 4.0.73e2 reconcilation
595 *
596 * Revision 4.3  1998/06/13 11:56:19  kardel
597 * disabled putbute() for the time being
598 *
599 * Revision 4.2  1998/06/12 15:16:58  kardel
600 * ansi2knr compatibility
601 *
602 * Revision 4.1  1998/05/24 07:59:56  kardel
603 * conditional debug support
604 *
605 * Revision 4.0  1998/04/10 19:46:29  kardel
606 * Start 4.0 release version numbering
607 *
608 * Revision 1.1  1998/04/10 19:27:46  kardel
609 * initial NTP VERSION 4 integration of PARSE with GPS166 binary support
610 *
611 * Revision 1.1  1997/10/06 21:05:45  kardel
612 * new parse structure
613 *
614 */
615