Subject: [NEWS] Re: fp-errors / compiling gcc 2.3.3 on i386.MACH

From: [email-masked] (Lars Koeller,D0-231/258,1065375,0526189607)
Date: Mon, 15 Feb 1993 09:25:10 GMT

---

Hello!

Here are the patches you need. We have had the same problems some time ago,
and spent a lot of time to solve it. The solution is to write a new 'atof'
and 'cvt' (see in doprnt), cause the original code makes a lot of errors
in the case of big and small numbers.
If you have access to the libc-sources, replace the original doprnt.c, atofs.s
and atof.c files and build a new library. In the other case compile the .c
files with 'cc -c -O <filename>' and the assembler code with (use as 1.36!)
'/lib/cpp -E -DASSEMBLER -DLIBC_SCCS atofs.s | as -o atofs.o'. Then replace
these objects in /lib/libc.a with 'ar'. When everything is o.k. enquire is
happy now. Good luck and have fun. Bye.

			Stefan Foerster & Lars Koeller 

###############################################################################

FILE: <source-tree>/lib/libc/stdio/doprnt.c

/* 
 * Mach Operating System
 * Copyright (c) 1990 Carnegie-Mellon University
 * All rights reserved.  The CMU software License Agreement specifies
 * the terms and conditions for use and redistribution.
 */
/*
 * HISTORY
 * 10-Dec-92  Stefan Foerster (sf) at Universitaet Bielefeld
 *      Improved precision to at least 15 digits for %lf etc.
 *      Removed braindead 'denormal' handling
 *
 * 14-Sep-90  Robert Baron (rvb) at Carnegie-Mellon University
 *	Handle ieee Infinity and Nan.
 *
 * $Log:	pmap.c,v $
 */

/*
 * Copyright (c) 1988 Regents of the University of California.
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms are permitted
 * provided that the above copyright notice and this paragraph are
 * duplicated in all such forms and that any documentation,
 * advertising materials, and other materials related to such
 * distribution and use acknowledge that the software was developed
 * by the University of California, Berkeley.  The name of the
 * University may not be used to endorse or promote products derived
 * from this software without specific prior written permission.
 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
 * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 */

#if defined(LIBC_SCCS) && !defined(lint)
static char sccsid[] = "@(#)doprnt.c	5.35 (Berkeley) 6/27/88";
#endif /* LIBC_SCCS and not lint */

#include <sys/types.h>
#include <varargs.h>
#include <stdio.h>
#include <ctype.h>

/* 11-bit exponent (VAX G floating point) is 308 decimal digits */
#define	MAXEXP		308
/* 128 bit fraction takes up 39 decimal digits; max reasonable precision */
#define	MAXFRACT	39

#define	DEFPREC		6

#define	BUF		(MAXEXP+MAXFRACT+1)	/* + decimal point */

#define	PUTC(ch)	(void) putc(ch, fp)

#define	ARG() \
	_ulong = flags&LONGINT ? va_arg(argp, long) : \
	    flags&SHORTINT ? va_arg(argp, short) : va_arg(argp, int);

#define	todigit(c)	((c) - '0')
#define	tochar(n)	((n) + '0')

/* have to deal with the negative buffer count kludge */
#define	NEGATIVE_COUNT_KLUDGE

#define	LONGINT		0x01		/* long integer */
#define	LONGDBL		0x02		/* long double; unimplemented */
#define	SHORTINT	0x04		/* short integer */
#define	ALT		0x08		/* alternate form */
#define	LADJUST		0x10		/* left adjustment */
#define	ZEROPAD		0x20		/* zero (as opposed to blank) pad */
#define	HEXPREFIX	0x40		/* add 0x or 0X prefix */

static int cvt();
static char *round();
extern void _f2bcd__();

_doprnt(fmt0, argp, fp)
	u_char *fmt0;
	va_list argp;
	register FILE *fp;
{
	register u_char *fmt;	/* format string */
	register int ch;	/* character from fmt */
	register int cnt;	/* return value accumulator */
	register int n;		/* random handy integer */
	register char *t;	/* buffer pointer */
	double _double;		/* double precision arguments %[eEfgG] */
	u_long _ulong;		/* integer arguments %[diouxX] */
	int base;		/* base for [diouxX] conversion */
	int dprec;		/* decimal precision in [diouxX] */
	int fieldsz;		/* field size expanded by sign, etc */
	int flags;		/* flags as above */
	int fpprec;		/* `extra' floating precision in [eEfgG] */
	int prec;		/* precision from format (%.3d), or -1 */
	int realsz;		/* field size expanded by decimal precision */
	int size;		/* size of converted field or string */
	int width;		/* width from format (%8d), or 0 */
	char sign;		/* sign prefix (' ', '+', '-', or \0) */
	char softsign;		/* temporary negative sign for floats */
	char *digs;		/* digits for [diouxX] conversion */
	char buf[BUF];		/* space for %c, %[diouxX], %[eEfgG] */

	if (fp->_flag & _IORW) {
		fp->_flag |= _IOWRT;
		fp->_flag &= ~(_IOEOF|_IOREAD);
	}
	if ((fp->_flag & _IOWRT) == 0)
		return (EOF);

	fmt = fmt0;
	digs = "0123456789abcdef";
	for (cnt = 0;; ++fmt) {
		n = fp->_cnt;
		for (t = (char *)fp->_ptr; (ch = *fmt) && ch != '%';
		     ++cnt, ++fmt)
			if (--n < 0
#ifdef NEGATIVE_COUNT_KLUDGE
			    && (!(fp->_flag & _IOLBF) || -n >= fp->_bufsiz)
#endif
			    || ch == '\n' && fp->_flag & _IOLBF) {
				fp->_cnt = n;
				fp->_ptr = t;
				(void) _flsbuf((u_char)ch, fp);
				n = fp->_cnt;
				t = (char *)fp->_ptr;
			} else
				*t++ = ch;
		fp->_cnt = n;
		fp->_ptr = t;
		if (!ch)
			return (cnt);

		flags = 0; dprec = 0; fpprec = 0; width = 0;
		prec = -1;
		sign = '\0';

rflag:		switch (*++fmt) {
		case ' ':
			/*
			 * ``If the space and + flags both appear, the space
			 * flag will be ignored.''
			 *	-- ANSI X3J11
			 */
			if (!sign)
				sign = ' ';
			goto rflag;
		case '#':
			flags |= ALT;
			goto rflag;
		case '*':
			/*
			 * ``A negative field width argument is taken as a
			 * - flag followed by a  positive field width.''
			 *	-- ANSI X3J11
			 * They don't exclude field widths read from args.
			 */
			if ((width = va_arg(argp, int)) >= 0)
				goto rflag;
			width = -width;
			/* FALLTHROUGH */
		case '-':
			flags |= LADJUST;
			goto rflag;
		case '+':
			sign = '+';
			goto rflag;
		case '.':
			if (*++fmt == '*')
				n = va_arg(argp, int);
			else {
				n = 0;
				while (isascii(*fmt) && isdigit(*fmt))
					n = 10 * n + todigit(*fmt++);
				--fmt;
			}
			prec = n < 0 ? -1 : n;
			goto rflag;
		case '0':
			/*
			 * ``Note that 0 is taken as a flag, not as the
			 * beginning of a field width.''
			 *	-- ANSI X3J11
			 */
			flags |= ZEROPAD;
			goto rflag;
		case '1': case '2': case '3': case '4':
		case '5': case '6': case '7': case '8': case '9':
			n = 0;
			do {
				n = 10 * n + todigit(*fmt);
			} while (isascii(*++fmt) && isdigit(*fmt));
			width = n;
			--fmt;
			goto rflag;
		case 'L':
			flags |= LONGDBL;
			goto rflag;
		case 'h':
			flags |= SHORTINT;
			goto rflag;
		case 'l':
			flags |= LONGINT;
			goto rflag;
		case 'c':
			*(t = buf) = va_arg(argp, int);
			size = 1;
			sign = '\0';
			goto pforw;
		case 'D':
			flags |= LONGINT;
			/*FALLTHROUGH*/
		case 'd':
		case 'i':
			ARG();
			if ((long)_ulong < 0) {
				_ulong = -_ulong;
				sign = '-';
			}
			base = 10;
			goto number;
		case 'e':
		case 'E':
		case 'f':
		case 'g':
		case 'G':
			_double = va_arg(argp, double);
			/*
			 * don't do unrealistic precision; just pad it with
			 * zeroes later, so buffer size stays rational.
			 */
			if (prec > MAXFRACT) {
				if (*fmt != 'g' && *fmt != 'G' || (flags&ALT))
					fpprec = prec - MAXFRACT;
				prec = MAXFRACT;
			}
			else if (prec == -1)
				prec = DEFPREC;
			/*
			 * softsign avoids negative 0 if _double is < 0 and
			 * no significant digits will be shown
			 */
#ifdef	i386
			if (!(size = special_floats(_double, buf, &softsign))) {
#endif	i386
			if (_double < 0) {
				softsign = '-';
				_double = -_double;
			}
			else
				softsign = 0;
			/*
			 * cvt may have to round up past the "start" of the
			 * buffer, i.e. ``intf("%.2f", (double)9.999);'';
			 * if the first char isn't NULL, it did.
			 */
			*buf = NULL;
			size = cvt(_double, prec, flags, &softsign, *fmt, buf,
			    buf + sizeof(buf));
#ifdef	i386
			}
#endif	i386
			if (softsign)
				sign = '-';
			t = *buf ? buf : buf + 1;
			goto pforw;
		case 'n':
			if (flags & LONGINT)
				*va_arg(argp, long *) = cnt;
			else if (flags & SHORTINT)
				*va_arg(argp, short *) = cnt;
			else
				*va_arg(argp, int *) = cnt;
			break;
		case 'O':
			flags |= LONGINT;
			/*FALLTHROUGH*/
		case 'o':
			ARG();
			base = 8;
			goto nosign;
		case 'p':
			/*
			 * ``The argument shall be a pointer to void.  The
			 * value of the pointer is converted to a sequence
			 * of printable characters, in an implementation-
			 * defined manner.''
			 *	-- ANSI X3J11
			 */
			/* NOSTRICT */
			_ulong = (u_long)va_arg(argp, void *);
			base = 16;
			goto nosign;
		case 's':
			if (!(t = va_arg(argp, char *)))
				t = "(null)";
			if (prec >= 0) {
				/*
				 * can't use strlen; can only look for the
				 * NUL in the first `prec' characters, and
				 * strlen() will go further.
				 */
				char *p, *memchr();

				if (p = memchr(t, 0, prec)) {
					size = p - t;
					if (size > prec)
						size = prec;
				} else
					size = prec;
			} else
				size = strlen(t);
			sign = '\0';
			goto pforw;
		case 'U':
			flags |= LONGINT;
			/*FALLTHROUGH*/
		case 'u':
			ARG();
			base = 10;
			goto nosign;
		case 'X':
			digs = "0123456789ABCDEF";
			/* FALLTHROUGH */
		case 'x':
			ARG();
			base = 16;
			/* leading 0x/X only if non-zero */
			if (flags & ALT && _ulong != 0)
				flags |= HEXPREFIX;

			/* unsigned conversions */
nosign:			sign = '\0';
			/*
			 * ``... diouXx conversions ... if a precision is
			 * specified, the 0 flag will be ignored.''
			 *	-- ANSI X3J11
			 */
number:			if ((dprec = prec) >= 0)
				flags &= ~ZEROPAD;

			/*
			 * ``The result of converting a zero value with an
			 * explicit precision of zero is no characters.''
			 *	-- ANSI X3J11
			 */
			t = buf + BUF;
			if (_ulong != 0 || prec != 0) {
				do {
					*--t = digs[_ulong % base];
					_ulong /= base;
				} while (_ulong);
				digs = "0123456789abcdef";
				if (flags & ALT && base == 8 && *t != '0')
					*--t = '0'; /* octal leading 0 */
			}
			size = buf + BUF - t;

pforw:
			/*
			 * All reasonable formats wind up here.  At this point,
			 * `t' points to a string which (if not flags&LADJUST)
			 * should be padded out to `width' places.  If
			 * flags&ZEROPAD, it should first be prefixed by any
			 * sign or other prefix; otherwise, it should be blank
			 * padded before the prefix is emitted.  After any
			 * left-hand padding and prefixing, emit zeroes
			 * required by a decimal [diouxX] precision, then print
			 * the string proper, then emit zeroes required by any
			 * leftover floating precision; finally, if LADJUST,
			 * pad with blanks.
			 */

			/*
			 * compute actual size, so we know how much to pad
			 * fieldsz excludes decimal prec; realsz includes it
			 */
			fieldsz = size + fpprec;
			if (sign)
				fieldsz++;
			if (flags & HEXPREFIX)
				fieldsz += 2;
			realsz = dprec > fieldsz ? dprec : fieldsz;

			/* right-adjusting blank padding */
			if ((flags & (LADJUST|ZEROPAD)) == 0 && width)
				for (n = realsz; n < width; n++)
					PUTC(' ');
			/* prefix */
			if (sign)
				PUTC(sign);
			if (flags & HEXPREFIX) {
				PUTC('0');
				PUTC((char)*fmt);
			}
			/* right-adjusting zero padding */
			if ((flags & (LADJUST|ZEROPAD)) == ZEROPAD)
				for (n = realsz; n < width; n++)
					PUTC('0');
			/* leading zeroes from decimal precision */
			for (n = fieldsz; n < dprec; n++)
				PUTC('0');

			/* the string or number proper */
			n = size;

			/* code below don't work !!
			if (fp->_cnt - (n = size) >= 0 &&
			    (fp->_flag & _IOLBF) == 0) {
				fp->_cnt -= n;
				bcopy(t, (char *)fp->_ptr, n);
				fp->_ptr += n;
			} else
				while (--n >= 0)
					PUTC(*t++);
			  end of shit !! */

			while ((*t != *fmt) && (--n >= 0)) PUTC(*t++);
			/* trailing f.p. zeroes */
			while (--fpprec >= 0) PUTC('0');
			while (--n >= 0) PUTC(*t++);

			/* left-adjusting padding (always blank) */
			if (flags & LADJUST)
				for (n = realsz; n < width; n++)
					PUTC(' ');
			/* finally, adjust cnt */
			cnt += width > realsz ? width : realsz;
			break;
		case '\0':	/* "%?" prints ?, unless ? is NULL */
			return (cnt);
		default:
			PUTC((char)*fmt);
			cnt++;
		}
	}
	/* NOTREACHED */
}

#ifdef	i386
special_floats(dbl, str, signp)
double dbl;
char *str;
char *signp;
{
	union {double d; int i[2];} dd;

	dd.d = dbl;

	if (!dd.i[0] && !dd.i[1])
		return 0;
	else if ((dd.i[1] & 0x7ff00000) == 0x7ff00000) {
		*signp = (dd.i[1] & 0x80000000) ? '-' : 0 ;
		if (!dd.i[0] && !(dd.i[1] & 0xfffff)) {
			strcpy(str, "Infinity");
			return 8;
		} else {
			strcpy(str, "NaN");
			return 3;
		}
/* braindamaged code !! */
/*	} else 	if (!(dd.i[1] & 0x7ff00000)) {
		*signp = (dd.i[1] & 0x80000000) ? '-' : 0 ;
		strcpy(str, "Denormal");
		return 8;
*/
	} else
		return 0;
}
#endif	i386

#define BCDTODIGIT(x) ((bcd[ (x) >> 1] >> 4*( (x) & 0x01)) & 0x0f) 

static
cvt(number, prec, flags, signp, fmtch, startp, endp)
     double number;
     register int prec;
     int flags;
     u_char fmtch;
     char *signp, *startp, *endp;
{
  char bcd[10];
  char *t, c, *p, gformat = 0;
  int  exp, digit;
  int  cnt, bcddigit = 16;

  t = ++startp;

  _f2bcd__( number, bcd, &exp );

  switch(fmtch) {

  case 'f':

    /* Vorkommastellen ausgeben */
    cnt = (exp >= 0) ? exp + 1 : 0;
    if (!cnt) *t++ = '0';
    while (cnt--)
      if (bcddigit >= 0) {
	*t++ = BCDTODIGIT(bcddigit) + '0';
	bcddigit--;
      } else *t++ = '0';

    /* Dezimalpunkt */
    if (prec || flags&ALT) *t++ = '.';

    /* Nachkommastellen */
    cnt = (exp < 0) ? -exp - 1 : 0;
    for (; cnt-- && prec ; *t++ = '0', prec--);
    for (; prec && (bcddigit >= 0); bcddigit--, prec--)
      *t++ = BCDTODIGIT(bcddigit) + '0';

    /* Runden wenn noch Ziffern uebrig */
    if ((bcddigit >= 0) && (cnt < 0))
      digit = BCDTODIGIT(bcddigit);
    else
      digit = 0;
    startp = round( digit, startp, t, NULL, signp );

    /* Nachkommastellen auffuellen */
    while(prec--) *t++ = '0';

    break;

  case 'e':
  case 'E':

    /* Vorkommastelle ausgeben */
ef: *t++ = BCDTODIGIT(bcddigit) + '0'; 
    bcddigit--;

    /* Dezimalpunkt */
    if (prec || flags&ALT) *t++ = '.';

    /* Nachkommastellen */
    for (; prec && (bcddigit >= 0) ; prec--, bcddigit--)
      *t++ = BCDTODIGIT(bcddigit) + '0'; 

    /* Runden wenn noch Ziffern uebrig */
    if (bcddigit >= 0) {
      digit = BCDTODIGIT(bcddigit);
      startp = round( digit, startp, t, &exp, signp );
    }

    /* Nachkommastellen auffuellen */
    for ( ; prec--; *t++ = '0' );

    if (gformat) {
      /* Nachfolgende Nullen beseitigen */
      p = t-1;
      for ( ; (*p == '0') && (p != startp) ; p--, t--);
      if (*p == '.') t--;
    }

    /* Exponent ausgeben */
    *t++ = fmtch;
    if (exp < 0) {
      exp = -exp;
      *t++ = '-';
    } else *t++ = '+';
    if ( !((c = tochar( exp / 100 )) == '0') ) *t++ = c;
    exp %= 100;
    *t++ = tochar( exp / 10 );
    exp %= 10;
    *t++ = tochar( exp );

    break;

  case 'g':
  case 'G':

    if (!prec) prec++;

    /* Exponentialformat wenn exp > prec oder exp < -4 */
    if (exp > prec || (exp < -4)) {
      --prec;
      fmtch -= 2;
      gformat = 1;
      goto ef;
    }

    /* %f-Format, aber prec gibt Anzahl signifikanter Stellen an */
    /* Vorkommastellen */
    cnt = (exp >= 0) ? exp + 1 : 0;
    if (cnt) {
      for ( ; prec && (bcddigit >= 0) && cnt; prec--, cnt--, bcddigit--)
	*t++ = BCDTODIGIT(bcddigit) + '0';
      if (!prec && (bcddigit >= 0)) {
	digit = BCDTODIGIT(bcddigit);
	startp = round( digit, startp, t, &exp, signp );
      }
      for ( ; cnt; cnt-- ) *t++ = '0';
    }

    /* Dezimalpunkt */
    if (prec || flags&ALT) *t++ = '.';

    /* Nachkommastellen */
    cnt = (exp < 0) ? -exp - 1 : 0;             /* LK!!! patch cause leading */
    for (; cnt-- && prec ; *t++ = '0', prec--); /* '0' are missing !!!!      */
    if (prec) {
      for ( ; prec && (bcddigit >= 0); prec--, bcddigit--)
	*t++ = BCDTODIGIT(bcddigit) + '0';

      /* Runden wenn noch Ziffern uebrig */
      if (bcddigit >= 0) {
	digit = BCDTODIGIT(bcddigit);
	startp = round( digit, startp, t, NULL, signp );
      }
    
      /* Nachfolgende Nullen beseitigen */
      p = t-1;
      for ( ; (*p == '0') && (p != startp) ; p--, t--);
      if (*p == '.') t--;
    }

    break;

  }
  return( t - startp );
}

static char *
round( digit, start, t, exp, signp )
     char *start, *t, *signp;
     int  digit, *exp;
{
  char *p;
  
  if (digit > 4) {
    for ( p = t-1; ; p-- ) {
      if (*p == '.') p--;
      if (++*p <= '9') break;
      *p = '0';
      if (p == start) {
	if (exp) {
	  *p = '1';
	  ++*exp;
	} else {
	  *--p = '1';
	  --start;
	}
	break;
      }
    }
  } else if (*signp == '-') {
    for ( p = t-1; ; p-- ) {
      if (*p == '.') --p;
      if (*p != '0') break;
      if (p == start) *signp = 0;
    }
  }
  return(start);
} 

###############################################################################

FILE: <source-tree>/lib/libc/gen/gen.i386/atof.c

#if defined(LIBC_SCCS) && !defined(lint)
static char sccsid[] = "@(#)atof.c      1.0 (Forester) 8/12/92";
#endif /* LIBC_SCCS and not lint */

#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>

extern double _bcd2f__();

double
atof(s)
     register const char *s;
{
  char bcd[10];
  char bcddigit = 17;
  char sign = 0, esign = 0;
  int  exp = 0;
  int  digits = 0, fracdigits = 0;

  bzero(bcd,10);

  while (isspace(*s)) s++;
  if ((sign=(*s == '-')) || (*s == '+')) s++;
  while (*s == '0') s++;

  while (isdigit(*s)) {
    digits++;
    if ( bcddigit >= 0) {
      bcd[bcddigit >> 1] |= (*s++ - '0') << 4*(bcddigit & 0x01);
      bcddigit--;
    } else s++;
  }

  if (*s == '.') {
    s++;
    if (!digits) 
      while (*s == '0') { 
	fracdigits++; 
	s++;
      }
    while (isdigit(*s)) {
      if ( bcddigit >= 0) {
	bcd[bcddigit >> 1] |= (*s++ - '0') << 4*(bcddigit & 0x01);
	bcddigit--;
      } else s++;
    }
  }

  if (toupper(*s) == 'E') {
    if ((esign=(*++s == '-')) || (*s == '+')) s++;
    while (isdigit(*s)) exp = exp*10 + (*s++ - '0');
    if (esign) exp = -exp;
  }

  if (digits) {
    exp += digits-1;
  } else if (fracdigits) {
    exp -= fracdigits+1;
  } else exp -= 1;

  bcd[9] = (sign) ? 0x80 : 0x00;

  return( _bcd2f__( bcd, exp ));
}

###############################################################################

FILE: <source-tree>/lib/libc/gen/gen.i386/atofs.s

#include <machine/asm.h>
/* double _bcd2f__( BCD *bcd, int exp ); */

i17:	.long	17

ENTRY(_bcd2f__)
	FRAME
	pushl	%eax
	movl	8(%ebp),%eax	
	  subl	$4,%esp
	fstcw	(%esp)
	  pushl (%esp)
	  andl  $0xf0ff,(%esp)
	  orl   $0x0300,(%esp)
	fldcw	(%esp)		
	fbld	(%eax)		
	fildl	12(%ebp)	
	fisubl	i17		
	fldl2t			
	fmulp	%st,%st(1)	
	fld1			
	fscale			
	fxch	%st(1)		
	fld	%st(0)		
	  orl	$0x0c00,(%esp)	
	fldcw	(%esp)		
	frndint
	  andl  $0xf3ff,(%esp)
	fldcw	(%esp)		
	fsubrp	%st,%st(1)	
	f2xm1			
	fld1			
	faddp	%st,%st(1)	
	fmulp	%st,%st(1)	
	fmulp	%st,%st(1)	
	
	  addl	$4,%esp
	fldcw	(%esp)
	  addl  $4,%esp
	popl	%eax
	EMARF
	ret

/* void _f2bcd__( double d, BCD *bcd, int *exp, int prec ); */
	
i1:	.long 1
i10:	.long 10
i1e16:	.quad 1000000000000000
i1e17:	.quad 10000000000000000
i1e18:	.quad 100000000000000000

ENTRY(_f2bcd__)
	FRAME
	pushl	%eax		
	  subl  $4,%esp		
	fstcw (%esp)		
	  pushl (%esp)		
	  andl	$0xf0ff,(%esp)	
	  orl	$0x0300,(%esp)	
	fldcw	(%esp)		
	fldl	8(%ebp)		
	fabs
	ftst
	fnstsw	%ax
	  sahf
	  je	bnull		
	fxtract			
	fxch	%st(1)		
	fldlg2			
	fmulp	%st,%st(1)	
	fld	%st(0)		
	  orl	$0x0c00,(%esp)	
	fldcw	(%esp)		
	frndint			
	  andl	$0xf3ff,(%esp)	
	fldcw	(%esp)		
	fld	%st(0)		
	fsubrp	%st,%st(2)	
	fxch	%st(1)		
	fldl2t			
	fmulp	%st,%st(1)	
	fld	%st(2)		
	fscale			
	fxch	%st(1)		
	fld	%st(0)		
	  orl	$0x0c00,(%esp)	
	fldcw	(%esp)		
	frndint			
	  andl	$0xf3ff,(%esp)
	fldcw	(%esp)		
	fsubrp	%st,%st(1)	
	f2xm1			
	fld1			
	faddp	%st,%st(1)	
	fmulp	%st,%st(1)	
	ficoml	i1
	fnstsw	%ax
	  sahf
	  jb    emi
	ficoml	i10
	fnstsw	%ax
	  sahf
	  jb	enu
	fildq	i1e16
	fmulp	%st,%st(1)
	fld1
	  jmp	eend
emi:	fildq	i1e18
	fmulp	%st,%st(1)
	fld1
	fchs
	  jmp	eend
enu:	fildq	i1e17
	fmulp	%st,%st(1)
	fldz
eend:	faddp	%st,%st(2)
	fxch	%st(1)
	jmp	bout

bnull:	fldz	
	fldz
bout:	movl	20(%ebp),%eax	
	fistpl	(%eax)		
	frndint
	movl	16(%ebp),%eax	
	fbstp	(%eax)
	fstp	%st(0)
	  addl  $4,%esp
	fldcw	(%esp)
	  addl	$4,%esp
	popl	%eax
	EMARF
	ret

###############################################################################

-------------------------------------------------------------------------------

Lars Koeller			    E-Mail: [email-masked]
Universitaet Bielefeld (Germany)	    [email-masked]
Fakultaet fuer Physik / D0-231	    Phone:  +49 521-106-5375
Universitaetsstr. 25		    Fax:    +49 521-106-5244
4800 Bielefeld 1		    Telex:  932 362