diff --git a/newlib/ChangeLog b/newlib/ChangeLog index f62c61946..dd4a3363a 100644 --- a/newlib/ChangeLog +++ b/newlib/ChangeLog @@ -1,3 +1,19 @@ +2006-06-22 Jeff Johnston + + * libc/stdlib/Makefile.am: Add new gdtoa routines. + * libc/stdlib/Makefile.in: Regenerated. + * libc/stdlib/gd_qnan.h: New file. + * libc/stdlib/gdtoa-gethex.c: Ditto. + * libc/stdlib/gdtoa-hexnan.c: Ditto. + * libc/stdlib/gdtoa.h: Ditto. + * libc/stdlib/mprec.c: Add new helper routines needed by + the new gdtoa code. + * libc/stdlib/mprec.h: Integrate some defines and prototypes + used by gdtoa routines here. + * libc/stdlib/strtod.c: Rebased on David M. Gay's gdtoa-strtod.c + which adds C99 support such as nan, inf, and hexadecimal input + format. + 2006-06-15 Corinna Vinschen * libc/include/stdio.h (__sgetc_r): Fix typo. diff --git a/newlib/libc/stdlib/Makefile.am b/newlib/libc/stdlib/Makefile.am index fdc52e04e..0874729c5 100644 --- a/newlib/libc/stdlib/Makefile.am +++ b/newlib/libc/stdlib/Makefile.am @@ -27,6 +27,8 @@ GENERAL_SOURCES = \ envlock.c \ eprintf.c \ exit.c \ + gdtoa-gethex.c \ + gdtoa-hexnan.c \ getenv.c \ getenv_r.c \ labs.c \ @@ -256,6 +258,8 @@ $(lpfx)mbtowc_r.$(oext): mbtowc_r.c mbctype.h $(lpfx)mprec.$(oext): mprec.c mprec.h $(lpfx)strtod.$(oext): strtod.c mprec.h +$(lpfx)gdtoa-gethex.$(oext): gdtoa-gethex.c mprec.h +$(lpfx)gdtoa-hexnan.$(oext): gdtoa-hexnan.c mprec.h $(lpfx)wctomb_r.$(oext): wctomb_r.c mbctype.h $(lpfx)drand48.$(oext): drand48.c rand48.h $(lpfx)erand48.$(oext): erand48.c rand48.h diff --git a/newlib/libc/stdlib/Makefile.in b/newlib/libc/stdlib/Makefile.in index 56f4a365f..0c6c5e75a 100644 --- a/newlib/libc/stdlib/Makefile.in +++ b/newlib/libc/stdlib/Makefile.in @@ -74,6 +74,7 @@ am__objects_1 = lib_a-__adjust.$(OBJEXT) lib_a-__atexit.$(OBJEXT) \ lib_a-dtoa.$(OBJEXT) lib_a-dtoastub.$(OBJEXT) \ lib_a-environ.$(OBJEXT) lib_a-envlock.$(OBJEXT) \ lib_a-eprintf.$(OBJEXT) lib_a-exit.$(OBJEXT) \ + lib_a-gdtoa-gethex.$(OBJEXT) lib_a-gdtoa-hexnan.$(OBJEXT) \ lib_a-getenv.$(OBJEXT) lib_a-getenv_r.$(OBJEXT) \ lib_a-labs.$(OBJEXT) lib_a-ldiv.$(OBJEXT) \ lib_a-ldtoa.$(OBJEXT) lib_a-malloc.$(OBJEXT) \ @@ -124,12 +125,12 @@ LTLIBRARIES = $(noinst_LTLIBRARIES) am__objects_7 = __adjust.lo __atexit.lo __call_atexit.lo __exp10.lo \ __ten_mu.lo _Exit.lo abort.lo abs.lo assert.lo atexit.lo \ atof.lo atoff.lo atoi.lo atol.lo calloc.lo div.lo dtoa.lo \ - dtoastub.lo environ.lo envlock.lo eprintf.lo exit.lo getenv.lo \ - getenv_r.lo labs.lo ldiv.lo ldtoa.lo malloc.lo mblen.lo \ - mblen_r.lo mbstowcs.lo mbstowcs_r.lo mbtowc.lo mbtowc_r.lo \ - mlock.lo mprec.lo mstats.lo rand.lo rand_r.lo realloc.lo \ - strtod.lo strtol.lo strtoul.lo wcstombs.lo wcstombs_r.lo \ - wctomb.lo wctomb_r.lo + dtoastub.lo environ.lo envlock.lo eprintf.lo exit.lo \ + gdtoa-gethex.lo gdtoa-hexnan.lo getenv.lo getenv_r.lo labs.lo \ + ldiv.lo ldtoa.lo malloc.lo mblen.lo mblen_r.lo mbstowcs.lo \ + mbstowcs_r.lo mbtowc.lo mbtowc_r.lo mlock.lo mprec.lo \ + mstats.lo rand.lo rand_r.lo realloc.lo strtod.lo strtol.lo \ + strtoul.lo wcstombs.lo wcstombs_r.lo wctomb.lo wctomb_r.lo am__objects_8 = cxa_atexit.lo cxa_finalize.lo drand48.lo ecvtbuf.lo \ efgcvt.lo erand48.lo jrand48.lo lcong48.lo lrand48.lo \ mrand48.lo msize.lo mtrim.lo nrand48.lo rand48.lo seed48.lo \ @@ -249,6 +250,7 @@ PACKAGE_TARNAME = @PACKAGE_TARNAME@ PACKAGE_VERSION = @PACKAGE_VERSION@ PATH_SEPARATOR = @PATH_SEPARATOR@ RANLIB = @RANLIB@ +READELF = @READELF@ SET_MAKE = @SET_MAKE@ SHELL = @SHELL@ STRIP = @STRIP@ @@ -259,6 +261,7 @@ ac_ct_AR = @ac_ct_AR@ ac_ct_AS = @ac_ct_AS@ ac_ct_CC = @ac_ct_CC@ ac_ct_RANLIB = @ac_ct_RANLIB@ +ac_ct_READELF = @ac_ct_READELF@ ac_ct_STRIP = @ac_ct_STRIP@ aext = @aext@ am__fastdepCC_FALSE = @am__fastdepCC_FALSE@ @@ -329,6 +332,8 @@ GENERAL_SOURCES = \ envlock.c \ eprintf.c \ exit.c \ + gdtoa-gethex.c \ + gdtoa-hexnan.c \ getenv.c \ getenv_r.c \ labs.c \ @@ -685,6 +690,18 @@ lib_a-exit.o: exit.c lib_a-exit.obj: exit.c $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-exit.obj `if test -f 'exit.c'; then $(CYGPATH_W) 'exit.c'; else $(CYGPATH_W) '$(srcdir)/exit.c'; fi` +lib_a-gdtoa-gethex.o: gdtoa-gethex.c + $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-gdtoa-gethex.o `test -f 'gdtoa-gethex.c' || echo '$(srcdir)/'`gdtoa-gethex.c + +lib_a-gdtoa-gethex.obj: gdtoa-gethex.c + $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-gdtoa-gethex.obj `if test -f 'gdtoa-gethex.c'; then $(CYGPATH_W) 'gdtoa-gethex.c'; else $(CYGPATH_W) '$(srcdir)/gdtoa-gethex.c'; fi` + +lib_a-gdtoa-hexnan.o: gdtoa-hexnan.c + $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-gdtoa-hexnan.o `test -f 'gdtoa-hexnan.c' || echo '$(srcdir)/'`gdtoa-hexnan.c + +lib_a-gdtoa-hexnan.obj: gdtoa-hexnan.c + $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-gdtoa-hexnan.obj `if test -f 'gdtoa-hexnan.c'; then $(CYGPATH_W) 'gdtoa-hexnan.c'; else $(CYGPATH_W) '$(srcdir)/gdtoa-hexnan.c'; fi` + lib_a-getenv.o: getenv.c $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-getenv.o `test -f 'getenv.c' || echo '$(srcdir)/'`getenv.c @@ -1298,6 +1315,8 @@ $(lpfx)mbtowc_r.$(oext): mbtowc_r.c mbctype.h $(lpfx)mprec.$(oext): mprec.c mprec.h $(lpfx)strtod.$(oext): strtod.c mprec.h +$(lpfx)gdtoa-gethex.$(oext): gdtoa-gethex.c mprec.h +$(lpfx)gdtoa-hexnan.$(oext): gdtoa-hexnan.c mprec.h $(lpfx)wctomb_r.$(oext): wctomb_r.c mbctype.h $(lpfx)drand48.$(oext): drand48.c rand48.h $(lpfx)erand48.$(oext): erand48.c rand48.h diff --git a/newlib/libc/stdlib/gd_qnan.h b/newlib/libc/stdlib/gd_qnan.h new file mode 100644 index 000000000..68f16cd12 --- /dev/null +++ b/newlib/libc/stdlib/gd_qnan.h @@ -0,0 +1,33 @@ +#ifdef __IEEE_BIG_ENDIAN + +#define f_QNAN 0x7fc00000 +#define d_QNAN0 0x7ff80000 +#define d_QNAN1 0x0 +#define ld_QNAN0 0x7ff80000 +#define ld_QNAN1 0x0 +#define ld_QNAN2 0x0 +#define ld_QNAN3 0x0 +#define ldus_QNAN0 0x7ff8 +#define ldus_QNAN1 0x0 +#define ldus_QNAN2 0x0 +#define ldus_QNAN3 0x0 +#define ldus_QNAN4 0x0 + +#elif defined(__IEEE_LITTLE_ENDIAN) + +#define f_QNAN 0xffc00000 +#define d_QNAN0 0x0 +#define d_QNAN1 0xfff80000 +#define ld_QNAN0 0x0 +#define ld_QNAN1 0xc0000000 +#define ld_QNAN2 0xffff +#define ld_QNAN3 0x0 +#define ldus_QNAN0 0x0 +#define ldus_QNAN1 0x0 +#define ldus_QNAN2 0x0 +#define ldus_QNAN3 0xc000 +#define ldus_QNAN4 0xffff + +#else +#error IEEE endian not defined +#endif diff --git a/newlib/libc/stdlib/gdtoa-gethex.c b/newlib/libc/stdlib/gdtoa-gethex.c new file mode 100644 index 000000000..92f30fca0 --- /dev/null +++ b/newlib/libc/stdlib/gdtoa-gethex.c @@ -0,0 +1,351 @@ +/**************************************************************** + +The author of this software is David M. Gay. + +Copyright (C) 1998 by Lucent Technologies +All Rights Reserved + +Permission to use, copy, modify, and distribute this software and +its documentation for any purpose and without fee is hereby +granted, provided that the above copyright notice appear in all +copies and that both that the copyright notice and this +permission notice and warranty disclaimer appear in supporting +documentation, and that the name of Lucent or any of its entities +not be used in advertising or publicity pertaining to +distribution of the software without specific, written prior +permission. + +LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, +INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. +IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY +SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES +WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER +IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, +ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF +THIS SOFTWARE. + +****************************************************************/ + +/* Please send bug reports to David M. Gay (dmg at acm dot org, + * with " at " changed at "@" and " dot " changed to "."). */ + +#include <_ansi.h> +#include +#include +#include "mprec.h" +#include "gdtoa.h" +#include "gd_qnan.h" + +#ifdef USE_LOCALE +#include "locale.h" +#endif + +unsigned char hexdig[256]; + +static void +_DEFUN (htinit, (h, s, inc), + unsigned char *h _AND + unsigned char *s _AND + int inc) +{ + int i, j; + for(i = 0; (j = s[i]) !=0; i++) + h[j] = i + inc; +} + +void +_DEFUN_VOID (hexdig_init) +{ +#define USC (unsigned char *) + htinit(hexdig, USC "0123456789", 0x10); + htinit(hexdig, USC "abcdef", 0x10 + 10); + htinit(hexdig, USC "ABCDEF", 0x10 + 10); +} + +static void +_DEFUN(rshift, (b, k), + _Bigint *b _AND + int k) +{ + __ULong *x, *x1, *xe, y; + int n; + + x = x1 = b->_x; + n = k >> kshift; + if (n < b->_wds) { + xe = x + b->_wds; + x += n; + if (k &= kmask) { + n = ULbits - k; + y = *x++ >> k; + while(x < xe) { + *x1++ = (y | (*x << n)) & ALL_ON; + y = *x++ >> k; + } + if ((*x1 = y) !=0) + x1++; + } + else + while(x < xe) + *x1++ = *x++; + } + if ((b->_wds = x1 - b->_x) == 0) + b->_x[0] = 0; +} + +static _Bigint * +_DEFUN (increment, (ptr, b), + struct _reent *ptr _AND + _Bigint *b) +{ + __ULong *x, *xe; + _Bigint *b1; +#ifdef Pack_16 + __ULong carry = 1, y; +#endif + + x = b->_x; + xe = x + b->_wds; +#ifdef Pack_32 + do { + if (*x < (__ULong)0xffffffffL) { + ++*x; + return b; + } + *x++ = 0; + } while(x < xe); +#else + do { + y = *x + carry; + carry = y >> 16; + *x++ = y & 0xffff; + if (!carry) + return b; + } while(x < xe); + if (carry) +#endif + { + if (b->_wds >= b->_maxwds) { + b1 = Balloc(ptr, b->_k+1); + Bcopy(b1, b); + Bfree(ptr, b); + b = b1; + } + b->_x[b->_wds++] = 1; + } + return b; +} + + +int +_DEFUN(gethex, (ptr, sp, fpi, exp, bp, sign), + struct _reent *ptr _AND + _CONST char **sp _AND + FPI *fpi _AND + Long *exp _AND + _Bigint **bp _AND + int sign) +{ + _Bigint *b; + _CONST unsigned char *decpt, *s0, *s, *s1; + int esign, havedig, irv, k, n, nbits, up, zret; + __ULong L, lostbits, *x; + Long e, e1; +#ifdef USE_LOCALE + unsigned char decimalpoint = *localeconv()->decimal_point; +#else +#define decimalpoint '.' +#endif + + if (!hexdig['0']) + hexdig_init(); + havedig = 0; + s0 = *(_CONST unsigned char **)sp + 2; + while(s0[havedig] == '0') + havedig++; + s0 += havedig; + s = s0; + decpt = 0; + zret = 0; + e = 0; + if (!hexdig[*s]) { + zret = 1; + if (*s != decimalpoint) + goto pcheck; + decpt = ++s; + if (!hexdig[*s]) + goto pcheck; + while(*s == '0') + s++; + if (hexdig[*s]) + zret = 0; + havedig = 1; + s0 = s; + } + while(hexdig[*s]) + s++; + if (*s == decimalpoint && !decpt) { + decpt = ++s; + while(hexdig[*s]) + s++; + } + if (decpt) + e = -(((Long)(s-decpt)) << 2); + pcheck: + s1 = s; + switch(*s) { + case 'p': + case 'P': + esign = 0; + switch(*++s) { + case '-': + esign = 1; + /* no break */ + case '+': + s++; + } + if ((n = hexdig[*s]) == 0 || n > 0x19) { + s = s1; + break; + } + e1 = n - 0x10; + while((n = hexdig[*++s]) !=0 && n <= 0x19) + e1 = 10*e1 + n - 0x10; + if (esign) + e1 = -e1; + e += e1; + } + *sp = (char*)s; + if (zret) + return havedig ? STRTOG_Zero : STRTOG_NoNumber; + n = s1 - s0 - 1; + for(k = 0; n > 7; n >>= 1) + k++; + b = Balloc(ptr, k); + x = b->_x; + n = 0; + L = 0; + while(s1 > s0) { + if (*--s1 == decimalpoint) + continue; + if (n == 32) { + *x++ = L; + L = 0; + n = 0; + } + L |= (hexdig[*s1] & 0x0f) << n; + n += 4; + } + *x++ = L; + b->_wds = n = x - b->_x; + n = 32*n - hi0bits(L); + nbits = fpi->nbits; + lostbits = 0; + x = b->_x; + if (n > nbits) { + n -= nbits; + if (any_on(b,n)) { + lostbits = 1; + k = n - 1; + if (x[k>>kshift] & 1 << (k & kmask)) { + lostbits = 2; + if (k > 1 && any_on(b,k-1)) + lostbits = 3; + } + } + rshift(b, n); + e += n; + } + else if (n < nbits) { + n = nbits - n; + b = lshift(ptr, b, n); + e -= n; + x = b->_x; + } + if (e > fpi->emax) { + ovfl: + Bfree(ptr, b); + *bp = 0; + return STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi; + } + irv = STRTOG_Normal; + if (e < fpi->emin) { + irv = STRTOG_Denormal; + n = fpi->emin - e; + if (n >= nbits) { + switch (fpi->rounding) { + case FPI_Round_near: + if (n == nbits && (n < 2 || any_on(b,n-1))) + goto one_bit; + break; + case FPI_Round_up: + if (!sign) + goto one_bit; + break; + case FPI_Round_down: + if (sign) { + one_bit: + *exp = fpi->emin; + x[0] = b->_wds = 1; + *bp = b; + return STRTOG_Denormal | STRTOG_Inexhi + | STRTOG_Underflow; + } + } + Bfree(ptr, b); + *bp = 0; + return STRTOG_Zero | STRTOG_Inexlo | STRTOG_Underflow; + } + k = n - 1; + if (lostbits) + lostbits = 1; + else if (k > 0) + lostbits = any_on(b,k); + if (x[k>>kshift] & 1 << (k & kmask)) + lostbits |= 2; + nbits -= n; + rshift(b,n); + e = fpi->emin; + } + if (lostbits) { + up = 0; + switch(fpi->rounding) { + case FPI_Round_zero: + break; + case FPI_Round_near: + if ((lostbits & 2) + && ((lostbits & 1) | (x[0] & 1))) + up = 1; + break; + case FPI_Round_up: + up = 1 - sign; + break; + case FPI_Round_down: + up = sign; + } + if (up) { + k = b->_wds; + b = increment(ptr, b); + x = b->_x; + if (irv == STRTOG_Denormal) { + if (nbits == fpi->nbits - 1 + && x[nbits >> kshift] & 1 << (nbits & kmask)) + irv = STRTOG_Normal; + } + else if ((b->_wds > k) + || ((n = nbits & kmask) !=0 + && (hi0bits(x[k-1]) < 32-n))) { + rshift(b,1); + if (++e > fpi->emax) + goto ovfl; + } + irv |= STRTOG_Inexhi; + } + else + irv |= STRTOG_Inexlo; + } + *bp = b; + *exp = e; + return irv; +} + diff --git a/newlib/libc/stdlib/gdtoa-hexnan.c b/newlib/libc/stdlib/gdtoa-hexnan.c new file mode 100644 index 000000000..058bbd17d --- /dev/null +++ b/newlib/libc/stdlib/gdtoa-hexnan.c @@ -0,0 +1,142 @@ +/**************************************************************** + +The author of this software is David M. Gay. + +Copyright (C) 2000 by Lucent Technologies +All Rights Reserved + +Permission to use, copy, modify, and distribute this software and +its documentation for any purpose and without fee is hereby +granted, provided that the above copyright notice appear in all +copies and that both that the copyright notice and this +permission notice and warranty disclaimer appear in supporting +documentation, and that the name of Lucent or any of its entities +not be used in advertising or publicity pertaining to +distribution of the software without specific, written prior +permission. + +LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, +INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. +IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY +SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES +WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER +IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, +ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF +THIS SOFTWARE. + +****************************************************************/ + +/* Please send bug reports to + David M. Gay + Bell Laboratories, Room 2C-463 + 600 Mountain Avenue + Murray Hill, NJ 07974-0636 + U.S.A. + dmg@bell-labs.com + */ + +/* Modified 06-21-2006 by Jeff Johnston to work with newlib. */ + +#include <_ansi.h> +#include +#include +#include "mprec.h" +#include "gdtoa.h" + +#ifdef INFNAN_CHECK +static void +_DEFUN (L_shift, (x, x1, i), + __ULong *x _AND + __ULong *x1 _AND + int i) +{ + int j; + + i = 8 - i; + i <<= 2; + j = ULbits - i; + do { + *x |= x[1] << j; + x[1] >>= i; + } while(++x < x1); +} + +int +_DEFUN (hexnan, (sp, fpi, x0), + _CONST char **sp _AND + FPI *fpi _AND + __ULong *x0) +{ + __ULong c, h, *x, *x1, *xe; + _CONST char *s; + int havedig, hd0, i, nbits; + + if (!hexdig['0']) + hexdig_init(); + nbits = fpi->nbits; + x = x0 + (nbits >> kshift); + if (nbits & kmask) + x++; + *--x = 0; + x1 = xe = x; + havedig = hd0 = i = 0; + s = *sp; + while((c = *(_CONST unsigned char*)++s)) { + if (!(h = hexdig[c])) { + if (c <= ' ') { + if (hd0 < havedig) { + if (x < x1 && i < 8) + L_shift(x, x1, i); + if (x <= x0) { + i = 8; + continue; + } + hd0 = havedig; + *--x = 0; + x1 = x; + i = 0; + } + continue; + } + if (/*(*/ c == ')' && havedig) { + *sp = s + 1; + break; + } + return STRTOG_NaN; + } + havedig++; + if (++i > 8) { + if (x <= x0) + continue; + i = 1; + *--x = 0; + } + *x = ((*x << 4) | (h & 0xf)); + } + if (!havedig) + return STRTOG_NaN; + if (x < x1 && i < 8) + L_shift(x, x1, i); + if (x > x0) { + x1 = x0; + do *x1++ = *x++; + while(x <= xe); + do *x1++ = 0; + while(x1 <= xe); + } + else { + /* truncate high-order word if necessary */ + if ( (i = nbits & (ULbits-1)) !=0) + *xe &= ((__ULong)0xffffffff) >> (ULbits - i); + } + for(x1 = xe;; --x1) { + if (*x1 != 0) + break; + if (x1 == x0) { + *x1 = 1; + break; + } + } + return STRTOG_NaNbits; +} +#endif /* INFNAN_CHECK */ diff --git a/newlib/libc/stdlib/gdtoa.h b/newlib/libc/stdlib/gdtoa.h new file mode 100644 index 000000000..5029e58de --- /dev/null +++ b/newlib/libc/stdlib/gdtoa.h @@ -0,0 +1,72 @@ +/**************************************************************** + +The author of this software is David M. Gay. + +Copyright (C) 1998 by Lucent Technologies +All Rights Reserved + +Permission to use, copy, modify, and distribute this software and +its documentation for any purpose and without fee is hereby +granted, provided that the above copyright notice appear in all +copies and that both that the copyright notice and this +permission notice and warranty disclaimer appear in supporting +documentation, and that the name of Lucent or any of its entities +not be used in advertising or publicity pertaining to +distribution of the software without specific, written prior +permission. + +LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, +INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. +IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY +SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES +WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER +IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, +ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF +THIS SOFTWARE. + +****************************************************************/ + +/* Please send bug reports to David M. Gay (dmg at acm dot org, + * with " at " changed at "@" and " dot " changed to "."). */ + +#ifndef GDTOA_H_INCLUDED +#define GDTOA_H_INCLUDED + + + enum { /* return values from strtodg */ + STRTOG_Zero = 0, + STRTOG_Normal = 1, + STRTOG_Denormal = 2, + STRTOG_Infinite = 3, + STRTOG_NaN = 4, + STRTOG_NaNbits = 5, + STRTOG_NoNumber = 6, + STRTOG_Retmask = 7, + + /* The following may be or-ed into one of the above values. */ + + STRTOG_Neg = 0x08, + STRTOG_Inexlo = 0x10, + STRTOG_Inexhi = 0x20, + STRTOG_Inexact = 0x30, + STRTOG_Underflow= 0x40, + STRTOG_Overflow = 0x80 + }; + + typedef struct +FPI { + int nbits; + int emin; + int emax; + int rounding; + int sudden_underflow; + } FPI; + +enum { /* FPI.rounding values: same as FLT_ROUNDS */ + FPI_Round_zero = 0, + FPI_Round_near = 1, + FPI_Round_up = 2, + FPI_Round_down = 3 + }; + +#endif /* GDTOA_H_INCLUDED */ diff --git a/newlib/libc/stdlib/mprec.c b/newlib/libc/stdlib/mprec.c index 0ef28c745..6e84ece5b 100644 --- a/newlib/libc/stdlib/mprec.c +++ b/newlib/libc/stdlib/mprec.c @@ -985,3 +985,61 @@ _DEFUN (_mprec_log10, (dig), } return v; } + +void +_DEFUN (copybits, (c, n, b), + __ULong *c _AND + int n _AND + _Bigint *b) +{ + __ULong *ce, *x, *xe; +#ifdef Pack_16 + int nw, nw1; +#endif + + ce = c + ((n-1) >> kshift) + 1; + x = b->_x; +#ifdef Pack_32 + xe = x + b->_wds; + while(x < xe) + *c++ = *x++; +#else + nw = b->_wds; + nw1 = nw & 1; + for(xe = x + (nw - nw1); x < xe; x += 2) + Storeinc(c, x[1], x[0]); + if (nw1) + *c++ = *x; +#endif + while(c < ce) + *c++ = 0; +} + +__ULong +_DEFUN (any_on, (b, k), + _Bigint *b _AND + int k) +{ + int n, nwds; + __ULong *x, *x0, x1, x2; + + x = b->_x; + nwds = b->_wds; + n = k >> kshift; + if (n > nwds) + n = nwds; + else if (n < nwds && (k &= kmask)) { + x1 = x2 = x[n]; + x1 >>= k; + x1 <<= k; + if (x1 != x2) + return 1; + } + x0 = x; + x += n; + while(x > x0) + if (*--x) + return 1; + return 0; +} + diff --git a/newlib/libc/stdlib/mprec.h b/newlib/libc/stdlib/mprec.h index 4ca48f22f..4fd496831 100644 --- a/newlib/libc/stdlib/mprec.h +++ b/newlib/libc/stdlib/mprec.h @@ -78,6 +78,39 @@ union double_union #define word1(x) (x.i[1]) #endif + +/* The following is taken from gdtoaimp.h for use with new strtod. */ +typedef __int32_t Long; +typedef union { double d; __ULong L[2]; } U; + +#ifdef YES_ALIAS +#define dval(x) x +#ifdef IEEE_8087 +#define dword0(x) ((__ULong *)&x)[1] +#define dword1(x) ((__ULong *)&x)[0] +#else +#define dword0(x) ((__ULong *)&x)[0] +#define dword1(x) ((__ULong *)&x)[1] +#endif +#else /* !YES_ALIAS */ +#ifdef IEEE_8087 +#define dword0(x) ((U*)&x)->L[1] +#define dword1(x) ((U*)&x)->L[0] +#else +#define dword0(x) ((U*)&x)->L[0] +#define dword1(x) ((U*)&x)->L[1] +#endif +#define dval(x) ((U*)&x)->d +#endif /* YES_ALIAS */ + + +#undef SI +#ifdef Sudden_Underflow +#define SI 1 +#else +#define SI 0 +#endif + /* The following definition of Storeinc is appropriate for MIPS processors. * An alternative that might be better on some machines is * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) @@ -161,12 +194,22 @@ union double_union #define Quick_max 14 #define Int_max 14 #define Infinite(x) (word0(x) == ((__uint32_t)0x7ff00000L)) /* sufficient test for here */ + +#ifndef Flt_Rounds +#ifdef FLT_ROUNDS +#define Flt_Rounds FLT_ROUNDS +#else +#define Flt_Rounds 1 +#endif +#endif /*Flt_Rounds*/ + #endif #else #undef Sudden_Underflow #define Sudden_Underflow #ifdef IBM +#define Flt_Rounds 0 #define Exp_shift 24 #define Exp_shift1 24 #define Exp_msk1 ((__uint32_t)0x1000000L) @@ -191,6 +234,7 @@ union double_union #define Quick_max 14 #define Int_max 15 #else /* VAX */ +#define Flt_Rounds 1 #define Exp_shift 23 #define Exp_shift1 7 #define Exp_msk1 0x80 @@ -219,6 +263,58 @@ union double_union #ifndef IEEE_Arith #define ROUND_BIASED +#else +#define Scale_Bit 0x10 +#if defined(_DOUBLE_IS_32BITS) && defined(__v800) +#define n_bigtens 2 +#else +#define n_bigtens 5 +#endif +#endif + +#ifdef IBM +#define n_bigtens 3 +#endif + +#ifdef VAX +#define n_bigtens 2 +#endif + +#ifndef __NO_INFNAN_CHECK +#define INFNAN_CHECK +#endif + +/* + * NAN_WORD0 and NAN_WORD1 are only referenced in strtod.c. Prior to + * 20050115, they used to be hard-wired here (to 0x7ff80000 and 0, + * respectively), but now are determined by compiling and running + * qnan.c to generate gd_qnan.h, which specifies d_QNAN0 and d_QNAN1. + * Formerly gdtoaimp.h recommended supplying suitable -DNAN_WORD0=... + * and -DNAN_WORD1=... values if necessary. This should still work. + * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.) + */ +#ifdef IEEE_Arith +#ifdef IEEE_MC68k +#define _0 0 +#define _1 1 +#ifndef NAN_WORD0 +#define NAN_WORD0 d_QNAN0 +#endif +#ifndef NAN_WORD1 +#define NAN_WORD1 d_QNAN1 +#endif +#else +#define _0 1 +#define _1 0 +#ifndef NAN_WORD0 +#define NAN_WORD0 d_QNAN1 +#endif +#ifndef NAN_WORD1 +#define NAN_WORD1 d_QNAN0 +#endif +#endif +#else +#undef INFNAN_CHECK #endif #ifdef RND_PRODQUOT @@ -249,6 +345,17 @@ extern double rnd_prod(double, double), rnd_quot(double, double); #endif #endif +#ifdef Pack_32 +#define ULbits 32 +#define kshift 5 +#define kmask 31 +#define ALL_ON 0xffffffff +#else +#define ULbits 16 +#define kshift 4 +#define kmask 15 +#define ALL_ON 0xffff +#endif #ifdef __cplusplus extern "C" double strtod(const char *s00, char **se); @@ -261,26 +368,34 @@ typedef struct _Bigint _Bigint; #define Balloc _Balloc #define Bfree _Bfree -#define multadd _multadd -#define s2b _s2b -#define lo0bits _lo0bits -#define hi0bits _hi0bits -#define i2b _i2b -#define mult _multiply -#define pow5mult _pow5mult -#define lshift _lshift +#define multadd __multadd +#define s2b __s2b +#define lo0bits __lo0bits +#define hi0bits __hi0bits +#define i2b __i2b +#define mult __multiply +#define pow5mult __pow5mult +#define lshift __lshift #define cmp __mcmp #define diff __mdiff -#define ulp _ulp -#define b2d _b2d -#define d2b _d2b -#define ratio _ratio +#define ulp __ulp +#define b2d __b2d +#define d2b __d2b +#define ratio __ratio +#define any_on __any_on +#define gethex __gethex +#define copybits __copybits +#define hexnan __hexnan +#define hexdig_init __hexdig_init + +#define hexdig __hexdig #define tens __mprec_tens #define bigtens __mprec_bigtens #define tinytens __mprec_tinytens struct _reent ; +struct FPI; double _EXFUN(ulp,(double x)); double _EXFUN(b2d,(_Bigint *a , int *e)); _Bigint * _EXFUN(Balloc,(struct _reent *p, int k)); @@ -292,23 +407,25 @@ _Bigint * _EXFUN(mult, (struct _reent *, _Bigint *, _Bigint *)); _Bigint * _EXFUN(pow5mult, (struct _reent *, _Bigint *, int k)); int _EXFUN(hi0bits,(__ULong)); int _EXFUN(lo0bits,(__ULong *)); -_Bigint * _EXFUN(d2b,(struct _reent *p, double d, int *e, int *bits)); -_Bigint * _EXFUN(lshift,(struct _reent *p, _Bigint *b, int k)); -_Bigint * _EXFUN(diff,(struct _reent *p, _Bigint *a, _Bigint *b)); -int _EXFUN(cmp,(_Bigint *a, _Bigint *b)); - +_Bigint * _EXFUN(d2b,(struct _reent *p, double d, int *e, int *bits)); +_Bigint * _EXFUN(lshift,(struct _reent *p, _Bigint *b, int k)); +_Bigint * _EXFUN(diff,(struct _reent *p, _Bigint *a, _Bigint *b)); +int _EXFUN(cmp,(_Bigint *a, _Bigint *b)); +int _EXFUN(gethex,(struct _reent *p, _CONST char **sp, struct FPI *fpi, Long *exp, _Bigint **bp, int sign)); double _EXFUN(ratio,(_Bigint *a, _Bigint *b)); -#define Bcopy(x,y) memcpy((char *)&x->_sign, (char *)&y->_sign, y->_wds*sizeof(__Long) + 2*sizeof(int)) - -#if defined(_DOUBLE_IS_32BITS) && defined(__v800) -#define n_bigtens 2 -#else -#define n_bigtens 5 +__ULong _EXFUN(any_on,(_Bigint *b, int k)); +void _EXFUN(copybits,(__ULong *c, int n, _Bigint *b)); +void _EXFUN(hexdig_init,(void)); +#ifdef INFNAN_CHECK +int _EXFUN(hexnan,(_CONST char **sp, struct FPI *fpi, __ULong *x0)); #endif +#define Bcopy(x,y) memcpy((char *)&x->_sign, (char *)&y->_sign, y->_wds*sizeof(__Long) + 2*sizeof(int)) + extern _CONST double tinytens[]; extern _CONST double bigtens[]; extern _CONST double tens[]; +extern unsigned char hexdig[]; double _EXFUN(_mprec_log10,(int)); diff --git a/newlib/libc/stdlib/strtod.c b/newlib/libc/stdlib/strtod.c index 9b70dfc3c..57297e05b 100644 --- a/newlib/libc/stdlib/strtod.c +++ b/newlib/libc/stdlib/strtod.c @@ -70,37 +70,130 @@ Supporting OS subroutines required: <>, <>, <>, */ /**************************************************************** - * - * The author of this software is David M. Gay. - * - * Copyright (c) 1991 by AT&T. - * - * Permission to use, copy, modify, and distribute this software for any - * purpose without fee is hereby granted, provided that this entire notice - * is included in all copies of any software which is or includes a copy - * or modification of this software and in all copies of the supporting - * documentation for such software. - * - * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED - * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY - * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY - * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. - * - ***************************************************************/ -/* Please send bug reports to - David M. Gay - AT&T Bell Laboratories, Room 2C-463 - 600 Mountain Avenue - Murray Hill, NJ 07974-2070 - U.S.A. - dmg@research.att.com or research!dmg - */ +The author of this software is David M. Gay. + +Copyright (C) 1998-2001 by Lucent Technologies +All Rights Reserved + +Permission to use, copy, modify, and distribute this software and +its documentation for any purpose and without fee is hereby +granted, provided that the above copyright notice appear in all +copies and that both that the copyright notice and this +permission notice and warranty disclaimer appear in supporting +documentation, and that the name of Lucent or any of its entities +not be used in advertising or publicity pertaining to +distribution of the software without specific, written prior +permission. + +LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, +INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. +IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY +SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES +WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER +IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, +ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF +THIS SOFTWARE. + +****************************************************************/ + +/* Please send bug reports to David M. Gay (dmg at acm dot org, + * with " at " changed at "@" and " dot " changed to "."). */ + +/* Original file gdtoa-strtod.c Modified 06-21-2006 by Jeff Johnston to work within newlib. */ #include <_ansi.h> -#include +#include #include #include "mprec.h" +#include "gdtoa.h" +#include "gd_qnan.h" + +/* #ifndef NO_FENV_H */ +/* #include */ +/* #endif */ + +#ifdef USE_LOCALE +#include "locale.h" +#endif + +#ifdef IEEE_Arith +#ifndef NO_IEEE_Scale +#define Avoid_Underflow +#undef tinytens +/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */ +/* flag unnecessarily. It leads to a song and dance at the end of strtod. */ +static _CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, + 9007199254740992.e-256 + }; +#endif +#endif + +#ifdef Honor_FLT_ROUNDS +#define Rounding rounding +#undef Check_FLT_ROUNDS +#define Check_FLT_ROUNDS +#else +#define Rounding Flt_Rounds +#endif + +static void +_DEFUN (ULtod, (L, bits, exp, k), + __ULong *L _AND + __ULong *bits _AND + Long exp _AND + int k) +{ + switch(k & STRTOG_Retmask) { + case STRTOG_NoNumber: + case STRTOG_Zero: + L[0] = L[1] = 0; + break; + + case STRTOG_Denormal: + L[_1] = bits[0]; + L[_0] = bits[1]; + break; + + case STRTOG_Normal: + case STRTOG_NaNbits: + L[_1] = bits[0]; + L[_0] = (bits[1] & ~0x100000) | ((exp + 0x3ff + 52) << 20); + break; + + case STRTOG_Infinite: + L[_0] = 0x7ff00000; + L[_1] = 0; + break; + + case STRTOG_NaN: + L[_0] = 0x7fffffff; + L[_1] = (__ULong)-1; + } + if (k & STRTOG_Neg) + L[_0] |= 0x80000000L; +} + +#ifdef INFNAN_CHECK +static int +_DEFUN (match, (sp, t), + _CONST char **sp _AND + char *t) +{ + int c, d; + _CONST char *s = *sp; + + while( (d = *t++) !=0) { + if ((c = *++s) >= 'A' && c <= 'Z') + c += 'a' - 'A'; + if (c != d) + return 0; + } + *sp = s + 1; + return 1; +} +#endif /* INFNAN_CHECK */ + double _DEFUN (_strtod_r, (ptr, s00, se), @@ -108,627 +201,919 @@ _DEFUN (_strtod_r, (ptr, s00, se), _CONST char *s00 _AND char **se) { - int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, e1, esign, i, j, - k, nd, nd0, nf, nz, nz0, sign; - long e; - _CONST char *s, *s0, *s1, *s2; - double aadj, aadj1, adj; - long L; - unsigned long z; - __ULong y; - union double_union rv, rv0; - int nanflag; - - _Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; - sign = nz0 = nz = nanflag = 0; - rv.d = 0.; - for (s = s00;; s++) - switch (*s) - { - case '-': - sign = 1; - /* no break */ - case '+': - if (*++s) - goto break2; - /* no break */ - case 0: - s = s00; - goto ret; - case '\t': - case '\n': - case '\v': - case '\f': - case '\r': - case ' ': - continue; - default: - goto break2; - } -break2: - if (*s == 'n' || *s == 'N') - { - ++s; - if (*s == 'a' || *s == 'A') - { - ++s; - if (*s == 'n' || *s == 'N') - { - nanflag = 1; - ++s; - goto ret; - } - } - s = s00; - goto ret; - } - else if (*s == '0') - { - nz0 = 1; - while (*++s == '0'); - if (!*s) - goto ret; - } - s0 = s; - y = z = 0; - for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) - if (nd < 9) - y = 10 * y + c - '0'; - else if (nd < 16) - z = 10 * z + c - '0'; - nd0 = nd; - if (c == '.') - { - c = *++s; - if (!nd) - { - for (; c == '0'; c = *++s) - nz++; - if (c > '0' && c <= '9') - { - s0 = s; - nf += nz; - nz = 0; - goto have_dig; - } - goto dig_done; - } - for (; c >= '0' && c <= '9'; c = *++s) - { - have_dig: - nz++; - if (c -= '0') - { - nf += nz; - for (i = 1; i < nz; i++) - if (nd++ < 9) - y *= 10; - else if (nd <= DBL_DIG + 1) - z *= 10; - if (nd++ < 9) - y = 10 * y + c; - else if (nd <= DBL_DIG + 1) - z = 10 * z + c; - nz = 0; - } - } - } -dig_done: - e = 0; - if (c == 'e' || c == 'E') - { - if (!nd && !nz && !nz0) - { - s = s00; - goto ret; - } - s2 = s; - esign = 0; - switch (c = *++s) - { - case '-': - esign = 1; - case '+': - c = *++s; - } - if (c >= '0' && c <= '9') - { - while (c == '0') - c = *++s; - if (c > '0' && c <= '9') - { - e = c - '0'; - s1 = s; - while ((c = *++s) >= '0' && c <= '9') - e = 10 * e + c - '0'; - if (s - s1 > 8) - /* Avoid confusion from exponents - * so large that e might overflow. - */ - e = 9999999L; - if (esign) - e = -e; - } - else - e = 0; - } - else - s = s2; - } - if (!nd) - { - if (!nz && !nz0) - s = s00; - goto ret; - } - e1 = e -= nf; - - /* Now we have nd0 digits, starting at s0, followed by a - * decimal point, followed by nd-nd0 digits. The number we're - * after is the integer represented by those digits times - * 10**e */ - - if (!nd0) - nd0 = nd; - k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; - rv.d = y; - if (k > 9) - rv.d = tens[k - 9] * rv.d + z; - bd0 = 0; - if (nd <= DBL_DIG -#ifndef RND_PRODQUOT - && FLT_ROUNDS == 1 +#ifdef Avoid_Underflow + int scale; #endif - ) - { - if (!e) - goto ret; - if (e > 0) - { - if (e <= Ten_pmax) - { -#ifdef VAX - goto vax_ovfl_check; + int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign, + e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; + _CONST char *s, *s0, *s1; + double aadj, aadj1, adj, rv, rv0; + Long L; + __ULong y, z; + _Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; +#ifdef SET_INEXACT + int inexact, oldinexact; +#endif +#ifdef Honor_FLT_ROUNDS + int rounding; +#endif + + delta = bs = bd = NULL; + sign = nz0 = nz = decpt = 0; + dval(rv) = 0.; + for(s = s00;;s++) switch(*s) { + case '-': + sign = 1; + /* no break */ + case '+': + if (*++s) + goto break2; + /* no break */ + case 0: + goto ret0; + case '\t': + case '\n': + case '\v': + case '\f': + case '\r': + case ' ': + continue; + default: + goto break2; + } + break2: + if (*s == '0') { +#ifndef NO_HEX_FP + { + static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; + Long exp; + __ULong bits[2]; + switch(s[1]) { + case 'x': + case 'X': + { +#if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD) + FPI fpi1 = fpi; + switch(fegetround()) { + case FE_TOWARDZERO: fpi1.rounding = 0; break; + case FE_UPWARD: fpi1.rounding = 2; break; + case FE_DOWNWARD: fpi1.rounding = 3; + } #else - /* rv.d = */ rounded_product (rv.d, tens[e]); - goto ret; +#define fpi1 fpi #endif - } - i = DBL_DIG - nd; - if (e <= Ten_pmax + i) - { - /* A fancier test would sometimes let us do + switch((i = gethex(ptr, &s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) { + case STRTOG_NoNumber: + s = s00; + sign = 0; + case STRTOG_Zero: + break; + default: + if (bb) { + copybits(bits, fpi.nbits, bb); + Bfree(ptr,bb); + } + ULtod(((U*)&rv)->L, bits, exp, i); + }} + goto ret; + } + } +#endif + nz0 = 1; + while(*++s == '0') ; + if (!*s) + goto ret; + } + s0 = s; + y = z = 0; + for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) + if (nd < 9) + y = 10*y + c - '0'; + else if (nd < 16) + z = 10*z + c - '0'; + nd0 = nd; +#ifdef USE_LOCALE + if (c == *localeconv()->decimal_point) +#else + if (c == '.') +#endif + { + decpt = 1; + c = *++s; + if (!nd) { + for(; c == '0'; c = *++s) + nz++; + if (c > '0' && c <= '9') { + s0 = s; + nf += nz; + nz = 0; + goto have_dig; + } + goto dig_done; + } + for(; c >= '0' && c <= '9'; c = *++s) { + have_dig: + nz++; + if (c -= '0') { + nf += nz; + for(i = 1; i < nz; i++) + if (nd++ < 9) + y *= 10; + else if (nd <= DBL_DIG + 1) + z *= 10; + if (nd++ < 9) + y = 10*y + c; + else if (nd <= DBL_DIG + 1) + z = 10*z + c; + nz = 0; + } + } + } + dig_done: + e = 0; + if (c == 'e' || c == 'E') { + if (!nd && !nz && !nz0) { + goto ret0; + } + s00 = s; + esign = 0; + switch(c = *++s) { + case '-': + esign = 1; + case '+': + c = *++s; + } + if (c >= '0' && c <= '9') { + while(c == '0') + c = *++s; + if (c > '0' && c <= '9') { + L = c - '0'; + s1 = s; + while((c = *++s) >= '0' && c <= '9') + L = 10*L + c - '0'; + if (s - s1 > 8 || L > 19999) + /* Avoid confusion from exponents + * so large that e might overflow. + */ + e = 19999; /* safe for 16 bit ints */ + else + e = (int)L; + if (esign) + e = -e; + } + else + e = 0; + } + else + s = s00; + } + if (!nd) { + if (!nz && !nz0) { +#ifdef INFNAN_CHECK + /* Check for Nan and Infinity */ + __ULong bits[2]; + static FPI fpinan = /* only 52 explicit bits */ + { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI }; + if (!decpt) + switch(c) { + case 'i': + case 'I': + if (match(&s,"nf")) { + --s; + if (!match(&s,"inity")) + ++s; + dword0(rv) = 0x7ff00000; + dword1(rv) = 0; + goto ret; + } + break; + case 'n': + case 'N': + if (match(&s, "an")) { +#ifndef No_Hex_NaN + if (*s == '(' /*)*/ + && hexnan(&s, &fpinan, bits) + == STRTOG_NaNbits) { + dword0(rv) = 0x7ff00000 | bits[1]; + dword1(rv) = bits[0]; + } + else { +#endif + dword0(rv) = NAN_WORD0; + dword1(rv) = NAN_WORD1; +#ifndef No_Hex_NaN + } +#endif + goto ret; + } + } +#endif /* INFNAN_CHECK */ + ret0: + s = s00; + sign = 0; + } + goto ret; + } + e1 = e -= nf; + + /* Now we have nd0 digits, starting at s0, followed by a + * decimal point, followed by nd-nd0 digits. The number we're + * after is the integer represented by those digits times + * 10**e */ + + if (!nd0) + nd0 = nd; + k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; + dval(rv) = y; + if (k > 9) { +#ifdef SET_INEXACT + if (k > DBL_DIG) + oldinexact = get_inexact(); +#endif + dval(rv) = tens[k - 9] * dval(rv) + z; + } + bd0 = 0; + if (nd <= DBL_DIG +#ifndef RND_PRODQUOT +#ifndef Honor_FLT_ROUNDS + && Flt_Rounds == 1 +#endif +#endif + ) { + if (!e) + goto ret; + if (e > 0) { + if (e <= Ten_pmax) { +#ifdef VAX + goto vax_ovfl_check; +#else +#ifdef Honor_FLT_ROUNDS + /* round correctly FLT_ROUNDS = 2 or 3 */ + if (sign) { + rv = -rv; + sign = 0; + } +#endif + /* rv = */ rounded_product(dval(rv), tens[e]); + goto ret; +#endif + } + i = DBL_DIG - nd; + if (e <= Ten_pmax + i) { + /* A fancier test would sometimes let us do * this for larger i values. */ - e -= i; - rv.d *= tens[i]; +#ifdef Honor_FLT_ROUNDS + /* round correctly FLT_ROUNDS = 2 or 3 */ + if (sign) { + rv = -rv; + sign = 0; + } +#endif + e -= i; + dval(rv) *= tens[i]; #ifdef VAX - /* VAX exponent range is so narrow we must - * worry about overflow here... - */ - vax_ovfl_check: - word0 (rv) -= P * Exp_msk1; - /* rv.d = */ rounded_product (rv.d, tens[e]); - if ((word0 (rv) & Exp_mask) - > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P)) - goto ovfl; - word0 (rv) += P * Exp_msk1; + /* VAX exponent range is so narrow we must + * worry about overflow here... + */ + vax_ovfl_check: + dword0(rv) -= P*Exp_msk1; + /* rv = */ rounded_product(dval(rv), tens[e]); + if ((dword0(rv) & Exp_mask) + > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) + goto ovfl; + dword0(rv) += P*Exp_msk1; #else - /* rv.d = */ rounded_product (rv.d, tens[e]); + /* rv = */ rounded_product(dval(rv), tens[e]); #endif - goto ret; - } - } + goto ret; + } + } #ifndef Inaccurate_Divide - else if (e >= -Ten_pmax) - { - /* rv.d = */ rounded_quotient (rv.d, tens[-e]); - goto ret; - } + else if (e >= -Ten_pmax) { +#ifdef Honor_FLT_ROUNDS + /* round correctly FLT_ROUNDS = 2 or 3 */ + if (sign) { + rv = -rv; + sign = 0; + } #endif - } - e1 += nd - k; + /* rv = */ rounded_quotient(dval(rv), tens[-e]); + goto ret; + } +#endif + } + e1 += nd - k; - /* Get starting approximation = rv.d * 10**e1 */ - - if (e1 > 0) - { - if ((i = e1 & 15) != 0) - rv.d *= tens[i]; - if (e1 &= ~15) - { - if (e1 > DBL_MAX_10_EXP) - { - ovfl: - ptr->_errno = ERANGE; -#ifdef _HAVE_STDC - rv.d = HUGE_VAL; -#else - /* Can't trust HUGE_VAL */ #ifdef IEEE_Arith - word0 (rv) = Exp_mask; -#ifndef _DOUBLE_IS_32BITS - word1 (rv) = 0; +#ifdef SET_INEXACT + inexact = 1; + if (k <= DBL_DIG) + oldinexact = get_inexact(); #endif -#else - word0 (rv) = Big0; -#ifndef _DOUBLE_IS_32BITS - word1 (rv) = Big1; -#endif -#endif -#endif - if (bd0) - goto retfree; - goto ret; - } - if (e1 >>= 4) - { - for (j = 0; e1 > 1; j++, e1 >>= 1) - if (e1 & 1) - rv.d *= bigtens[j]; - /* The last multiplication could overflow. */ - word0 (rv) -= P * Exp_msk1; - rv.d *= bigtens[j]; - if ((z = word0 (rv) & Exp_mask) - > Exp_msk1 * (DBL_MAX_EXP + Bias - P)) - goto ovfl; - if (z > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P)) - { - /* set to largest number */ - /* (Can't trust DBL_MAX) */ - word0 (rv) = Big0; -#ifndef _DOUBLE_IS_32BITS - word1 (rv) = Big1; +#ifdef Avoid_Underflow + scale = 0; #endif +#ifdef Honor_FLT_ROUNDS + if ((rounding = Flt_Rounds) >= 2) { + if (sign) + rounding = rounding == 2 ? 0 : 2; + else + if (rounding != 2) + rounding = 0; } - else - word0 (rv) += P * Exp_msk1; - } - - } - } - else if (e1 < 0) - { - e1 = -e1; - if ((i = e1 & 15) != 0) - rv.d /= tens[i]; - if (e1 &= ~15) - { - e1 >>= 4; - if (e1 >= 1 << n_bigtens) - goto undfl; - for (j = 0; e1 > 1; j++, e1 >>= 1) - if (e1 & 1) - rv.d *= tinytens[j]; - /* The last multiplication could underflow. */ - rv0.d = rv.d; - rv.d *= tinytens[j]; - if (!rv.d) - { - rv.d = 2. * rv0.d; - rv.d *= tinytens[j]; - if (!rv.d) - { - undfl: - rv.d = 0.; - ptr->_errno = ERANGE; - if (bd0) - goto retfree; - goto ret; - } -#ifndef _DOUBLE_IS_32BITS - word0 (rv) = Tiny0; - word1 (rv) = Tiny1; -#else - word0 (rv) = Tiny1; #endif - /* The refinement below will clean - * this approximation up. - */ - } - } - } +#endif /*IEEE_Arith*/ - /* Now the hard part -- adjusting rv to the correct value.*/ + /* Get starting approximation = rv * 10**e1 */ - /* Put digits into bd: true value = bd * 10^e */ + if (e1 > 0) { + if ( (i = e1 & 15) !=0) + dval(rv) *= tens[i]; + if (e1 &= ~15) { + if (e1 > DBL_MAX_10_EXP) { + ovfl: +#ifndef NO_ERRNO + ptr->_errno = ERANGE; +#endif + /* Can't trust HUGE_VAL */ +#ifdef IEEE_Arith +#ifdef Honor_FLT_ROUNDS + switch(rounding) { + case 0: /* toward 0 */ + case 3: /* toward -infinity */ + dword0(rv) = Big0; + dword1(rv) = Big1; + break; + default: + dword0(rv) = Exp_mask; + dword1(rv) = 0; + } +#else /*Honor_FLT_ROUNDS*/ + dword0(rv) = Exp_mask; + dword1(rv) = 0; +#endif /*Honor_FLT_ROUNDS*/ +#ifdef SET_INEXACT + /* set overflow bit */ + dval(rv0) = 1e300; + dval(rv0) *= dval(rv0); +#endif +#else /*IEEE_Arith*/ + dword0(rv) = Big0; + dword1(rv) = Big1; +#endif /*IEEE_Arith*/ + if (bd0) + goto retfree; + goto ret; + } + e1 >>= 4; + for(j = 0; e1 > 1; j++, e1 >>= 1) + if (e1 & 1) + dval(rv) *= bigtens[j]; + /* The last multiplication could overflow. */ + dword0(rv) -= P*Exp_msk1; + dval(rv) *= bigtens[j]; + if ((z = dword0(rv) & Exp_mask) + > Exp_msk1*(DBL_MAX_EXP+Bias-P)) + goto ovfl; + if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { + /* set to largest number */ + /* (Can't trust DBL_MAX) */ + dword0(rv) = Big0; + dword1(rv) = Big1; + } + else + dword0(rv) += P*Exp_msk1; + } + } + else if (e1 < 0) { + e1 = -e1; + if ( (i = e1 & 15) !=0) + dval(rv) /= tens[i]; + if (e1 >>= 4) { + if (e1 >= 1 << n_bigtens) + goto undfl; +#ifdef Avoid_Underflow + if (e1 & Scale_Bit) + scale = 2*P; + for(j = 0; e1 > 0; j++, e1 >>= 1) + if (e1 & 1) + dval(rv) *= tinytens[j]; + if (scale && (j = 2*P + 1 - ((dword0(rv) & Exp_mask) + >> Exp_shift)) > 0) { + /* scaled rv is denormal; zap j low bits */ + if (j >= 32) { + dword1(rv) = 0; + if (j >= 53) + dword0(rv) = (P+2)*Exp_msk1; + else + dword0(rv) &= 0xffffffff << (j-32); + } + else + dword1(rv) &= 0xffffffff << j; + } +#else + for(j = 0; e1 > 1; j++, e1 >>= 1) + if (e1 & 1) + dval(rv) *= tinytens[j]; + /* The last multiplication could underflow. */ + dval(rv0) = dval(rv); + dval(rv) *= tinytens[j]; + if (!dval(rv)) { + dval(rv) = 2.*dval(rv0); + dval(rv) *= tinytens[j]; +#endif + if (!dval(rv)) { + undfl: + dval(rv) = 0.; +#ifndef NO_ERRNO + ptr->_errno = ERANGE; +#endif + if (bd0) + goto retfree; + goto ret; + } +#ifndef Avoid_Underflow + dword0(rv) = Tiny0; + dword1(rv) = Tiny1; + /* The refinement below will clean + * this approximation up. + */ + } +#endif + } + } - bd0 = s2b (ptr, s0, nd0, nd, y); + /* Now the hard part -- adjusting rv to the correct value.*/ - for (;;) - { - bd = Balloc (ptr, bd0->_k); - Bcopy (bd, bd0); - bb = d2b (ptr, rv.d, &bbe, &bbbits); /* rv.d = bb * 2^bbe */ - bs = i2b (ptr, 1); + /* Put digits into bd: true value = bd * 10^e */ - if (e >= 0) - { - bb2 = bb5 = 0; - bd2 = bd5 = e; - } - else - { - bb2 = bb5 = -e; - bd2 = bd5 = 0; - } - if (bbe >= 0) - bb2 += bbe; - else - bd2 -= bbe; - bs2 = bb2; + bd0 = s2b(ptr, s0, nd0, nd, y); + + for(;;) { + bd = Balloc(ptr,bd0->_k); + Bcopy(bd, bd0); + bb = d2b(ptr,dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */ + bs = i2b(ptr,1); + + if (e >= 0) { + bb2 = bb5 = 0; + bd2 = bd5 = e; + } + else { + bb2 = bb5 = -e; + bd2 = bd5 = 0; + } + if (bbe >= 0) + bb2 += bbe; + else + bd2 -= bbe; + bs2 = bb2; +#ifdef Honor_FLT_ROUNDS + if (rounding != 1) + bs2++; +#endif +#ifdef Avoid_Underflow + j = bbe - scale; + i = j + bbbits - 1; /* logb(rv) */ + if (i < Emin) /* denormal */ + j += P - Emin; + else + j = P + 1 - bbbits; +#else /*Avoid_Underflow*/ #ifdef Sudden_Underflow #ifdef IBM - j = 1 + 4 * P - 3 - bbbits + ((bbe + bbbits - 1) & 3); + j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); #else - j = P + 1 - bbbits; + j = P + 1 - bbbits; #endif +#else /*Sudden_Underflow*/ + j = bbe; + i = j + bbbits - 1; /* logb(rv) */ + if (i < Emin) /* denormal */ + j += P - Emin; + else + j = P + 1 - bbbits; +#endif /*Sudden_Underflow*/ +#endif /*Avoid_Underflow*/ + bb2 += j; + bd2 += j; +#ifdef Avoid_Underflow + bd2 += scale; +#endif + i = bb2 < bd2 ? bb2 : bd2; + if (i > bs2) + i = bs2; + if (i > 0) { + bb2 -= i; + bd2 -= i; + bs2 -= i; + } + if (bb5 > 0) { + bs = pow5mult(ptr, bs, bb5); + bb1 = mult(ptr, bs, bb); + Bfree(ptr, bb); + bb = bb1; + } + if (bb2 > 0) + bb = lshift(ptr, bb, bb2); + if (bd5 > 0) + bd = pow5mult(ptr, bd, bd5); + if (bd2 > 0) + bd = lshift(ptr, bd, bd2); + if (bs2 > 0) + bs = lshift(ptr, bs, bs2); + delta = diff(ptr, bb, bd); + dsign = delta->_sign; + delta->_sign = 0; + i = cmp(delta, bs); +#ifdef Honor_FLT_ROUNDS + if (rounding != 1) { + if (i < 0) { + /* Error is less than an ulp */ + if (!delta->_x[0] && delta->_wds <= 1) { + /* exact */ +#ifdef SET_INEXACT + inexact = 0; +#endif + break; + } + if (rounding) { + if (dsign) { + adj = 1.; + goto apply_adj; + } + } + else if (!dsign) { + adj = -1.; + if (!dword1(rv) + && !(dword0(rv) & Frac_mask)) { + y = dword0(rv) & Exp_mask; +#ifdef Avoid_Underflow + if (!scale || y > 2*P*Exp_msk1) #else - i = bbe + bbbits - 1; /* logb(rv.d) */ - if (i < Emin) /* denormal */ - j = bbe + (P - Emin); - else - j = P + 1 - bbbits; + if (y) #endif - bb2 += j; - bd2 += j; - i = bb2 < bd2 ? bb2 : bd2; - if (i > bs2) - i = bs2; - if (i > 0) - { - bb2 -= i; - bd2 -= i; - bs2 -= i; - } - if (bb5 > 0) - { - bs = pow5mult (ptr, bs, bb5); - bb1 = mult (ptr, bs, bb); - Bfree (ptr, bb); - bb = bb1; - } - if (bb2 > 0) - bb = lshift (ptr, bb, bb2); - if (bd5 > 0) - bd = pow5mult (ptr, bd, bd5); - if (bd2 > 0) - bd = lshift (ptr, bd, bd2); - if (bs2 > 0) - bs = lshift (ptr, bs, bs2); - delta = diff (ptr, bb, bd); - dsign = delta->_sign; - delta->_sign = 0; - i = cmp (delta, bs); - if (i < 0) - { - /* Error is less than half an ulp -- check for - * special case of mantissa a power of two. - */ - if (dsign || word1 (rv) || word0 (rv) & Bndry_mask) - break; - delta = lshift (ptr, delta, Log2P); - if (cmp (delta, bs) > 0) - goto drop_down; - break; - } - if (i == 0) - { - /* exactly half-way between */ - if (dsign) - { - if ((word0 (rv) & Bndry_mask1) == Bndry_mask1 - && word1 (rv) == 0xffffffff) - { - /*boundary case -- increment exponent*/ - word0 (rv) = (word0 (rv) & Exp_mask) - + Exp_msk1 -#ifdef IBM - | Exp_msk1 >> 4 -#endif - ; -#ifndef _DOUBLE_IS_32BITS - word1 (rv) = 0; -#endif - break; - } - } - else if (!(word0 (rv) & Bndry_mask) && !word1 (rv)) - { - drop_down: - /* boundary case -- decrement exponent */ + { + delta = lshift(ptr, delta,Log2P); + if (cmp(delta, bs) <= 0) + adj = -0.5; + } + } + apply_adj: +#ifdef Avoid_Underflow + if (scale && (y = dword0(rv) & Exp_mask) + <= 2*P*Exp_msk1) + dword0(adj) += (2*P+1)*Exp_msk1 - y; +#else #ifdef Sudden_Underflow - L = word0 (rv) & Exp_mask; -#ifdef IBM - if (L < Exp_msk1) + if ((dword0(rv) & Exp_mask) <= + P*Exp_msk1) { + dword0(rv) += P*Exp_msk1; + dval(rv) += adj*ulp(dval(rv)); + dword0(rv) -= P*Exp_msk1; + } + else +#endif /*Sudden_Underflow*/ +#endif /*Avoid_Underflow*/ + dval(rv) += adj*ulp(dval(rv)); + } + break; + } + adj = ratio(delta, bs); + if (adj < 1.) + adj = 1.; + if (adj <= 0x7ffffffe) { + /* adj = rounding ? ceil(adj) : floor(adj); */ + y = adj; + if (y != adj) { + if (!((rounding>>1) ^ dsign)) + y++; + adj = y; + } + } +#ifdef Avoid_Underflow + if (scale && (y = dword0(rv) & Exp_mask) <= 2*P*Exp_msk1) + dword0(adj) += (2*P+1)*Exp_msk1 - y; #else - if (L <= Exp_msk1) -#endif - goto undfl; - L -= Exp_msk1; -#else - L = (word0 (rv) & Exp_mask) - Exp_msk1; -#endif - word0 (rv) = L | Bndry_mask1; -#ifndef _DOUBLE_IS_32BITS - word1 (rv) = 0xffffffff; -#endif -#ifdef IBM - goto cont; -#else - break; -#endif - } -#ifndef ROUND_BIASED - if (!(word1 (rv) & LSB)) - break; -#endif - if (dsign) - rv.d += ulp (rv.d); -#ifndef ROUND_BIASED - else - { - rv.d -= ulp (rv.d); -#ifndef Sudden_Underflow - if (!rv.d) - goto undfl; -#endif - } -#endif - break; - } - if ((aadj = ratio (delta, bs)) <= 2.) - { - if (dsign) - aadj = aadj1 = 1.; - else if (word1 (rv) || word0 (rv) & Bndry_mask) - { -#ifndef Sudden_Underflow - if (word1 (rv) == Tiny1 && !word0 (rv)) - goto undfl; -#endif - aadj = 1.; - aadj1 = -1.; - } - else - { - /* special case -- power of FLT_RADIX to be */ - /* rounded down... */ +#ifdef Sudden_Underflow + if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) { + dword0(rv) += P*Exp_msk1; + adj *= ulp(dval(rv)); + if (dsign) + dval(rv) += adj; + else + dval(rv) -= adj; + dword0(rv) -= P*Exp_msk1; + goto cont; + } +#endif /*Sudden_Underflow*/ +#endif /*Avoid_Underflow*/ + adj *= ulp(dval(rv)); + if (dsign) + dval(rv) += adj; + else + dval(rv) -= adj; + goto cont; + } +#endif /*Honor_FLT_ROUNDS*/ - if (aadj < 2. / FLT_RADIX) - aadj = 1. / FLT_RADIX; - else - aadj *= 0.5; - aadj1 = -aadj; - } - } - else - { - aadj *= 0.5; - aadj1 = dsign ? aadj : -aadj; + if (i < 0) { + /* Error is less than half an ulp -- check for + * special case of mantissa a power of two. + */ + if (dsign || dword1(rv) || dword0(rv) & Bndry_mask +#ifdef IEEE_Arith +#ifdef Avoid_Underflow + || (dword0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1 +#else + || (dword0(rv) & Exp_mask) <= Exp_msk1 +#endif +#endif + ) { +#ifdef SET_INEXACT + if (!delta->x[0] && delta->wds <= 1) + inexact = 0; +#endif + break; + } + if (!delta->_x[0] && delta->_wds <= 1) { + /* exact result */ +#ifdef SET_INEXACT + inexact = 0; +#endif + break; + } + delta = lshift(ptr,delta,Log2P); + if (cmp(delta, bs) > 0) + goto drop_down; + break; + } + if (i == 0) { + /* exactly half-way between */ + if (dsign) { + if ((dword0(rv) & Bndry_mask1) == Bndry_mask1 + && dword1(rv) == ( +#ifdef Avoid_Underflow + (scale && (y = dword0(rv) & Exp_mask) <= 2*P*Exp_msk1) + ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : +#endif + 0xffffffff)) { + /*boundary case -- increment exponent*/ + dword0(rv) = (dword0(rv) & Exp_mask) + + Exp_msk1 +#ifdef IBM + | Exp_msk1 >> 4 +#endif + ; + dword1(rv) = 0; +#ifdef Avoid_Underflow + dsign = 0; +#endif + break; + } + } + else if (!(dword0(rv) & Bndry_mask) && !dword1(rv)) { + drop_down: + /* boundary case -- decrement exponent */ +#ifdef Sudden_Underflow /*{{*/ + L = dword0(rv) & Exp_mask; +#ifdef IBM + if (L < Exp_msk1) +#else +#ifdef Avoid_Underflow + if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1)) +#else + if (L <= Exp_msk1) +#endif /*Avoid_Underflow*/ +#endif /*IBM*/ + goto undfl; + L -= Exp_msk1; +#else /*Sudden_Underflow}{*/ +#ifdef Avoid_Underflow + if (scale) { + L = dword0(rv) & Exp_mask; + if (L <= (2*P+1)*Exp_msk1) { + if (L > (P+2)*Exp_msk1) + /* round even ==> */ + /* accept rv */ + break; + /* rv = smallest denormal */ + goto undfl; + } + } +#endif /*Avoid_Underflow*/ + L = (dword0(rv) & Exp_mask) - Exp_msk1; +#endif /*Sudden_Underflow}*/ + dword0(rv) = L | Bndry_mask1; + dword1(rv) = 0xffffffff; +#ifdef IBM + goto cont; +#else + break; +#endif + } +#ifndef ROUND_BIASED + if (!(dword1(rv) & LSB)) + break; +#endif + if (dsign) + dval(rv) += ulp(dval(rv)); +#ifndef ROUND_BIASED + else { + dval(rv) -= ulp(dval(rv)); +#ifndef Sudden_Underflow + if (!dval(rv)) + goto undfl; +#endif + } +#ifdef Avoid_Underflow + dsign = 1 - dsign; +#endif +#endif + break; + } + if ((aadj = ratio(delta, bs)) <= 2.) { + if (dsign) + aadj = aadj1 = 1.; + else if (dword1(rv) || dword0(rv) & Bndry_mask) { +#ifndef Sudden_Underflow + if (dword1(rv) == Tiny1 && !dword0(rv)) + goto undfl; +#endif + aadj = 1.; + aadj1 = -1.; + } + else { + /* special case -- power of FLT_RADIX to be */ + /* rounded down... */ + + if (aadj < 2./FLT_RADIX) + aadj = 1./FLT_RADIX; + else + aadj *= 0.5; + aadj1 = -aadj; + } + } + else { + aadj *= 0.5; + aadj1 = dsign ? aadj : -aadj; #ifdef Check_FLT_ROUNDS - switch (FLT_ROUNDS) - { - case 2: /* towards +infinity */ - aadj1 -= 0.5; - break; - case 0: /* towards 0 */ - case 3: /* towards -infinity */ - aadj1 += 0.5; - } + switch(Rounding) { + case 2: /* towards +infinity */ + aadj1 -= 0.5; + break; + case 0: /* towards 0 */ + case 3: /* towards -infinity */ + aadj1 += 0.5; + } #else - if (FLT_ROUNDS == 0) - aadj1 += 0.5; -#endif - } - y = word0 (rv) & Exp_mask; + if (Flt_Rounds == 0) + aadj1 += 0.5; +#endif /*Check_FLT_ROUNDS*/ + } + y = dword0(rv) & Exp_mask; - /* Check for overflow */ + /* Check for overflow */ - if (y == Exp_msk1 * (DBL_MAX_EXP + Bias - 1)) - { - rv0.d = rv.d; - word0 (rv) -= P * Exp_msk1; - adj = aadj1 * ulp (rv.d); - rv.d += adj; - if ((word0 (rv) & Exp_mask) >= - Exp_msk1 * (DBL_MAX_EXP + Bias - P)) - { - if (word0 (rv0) == Big0 && word1 (rv0) == Big1) - goto ovfl; -#ifdef _DOUBLE_IS_32BITS - word0 (rv) = Big1; + if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { + dval(rv0) = dval(rv); + dword0(rv) -= P*Exp_msk1; + adj = aadj1 * ulp(dval(rv)); + dval(rv) += adj; + if ((dword0(rv) & Exp_mask) >= + Exp_msk1*(DBL_MAX_EXP+Bias-P)) { + if (dword0(rv0) == Big0 && dword1(rv0) == Big1) + goto ovfl; + dword0(rv) = Big0; + dword1(rv) = Big1; + goto cont; + } + else + dword0(rv) += P*Exp_msk1; + } + else { +#ifdef Avoid_Underflow + if (scale && y <= 2*P*Exp_msk1) { + if (aadj <= 0x7fffffff) { + if ((z = aadj) <= 0) + z = 1; + aadj = z; + aadj1 = dsign ? aadj : -aadj; + } + dword0(aadj1) += (2*P+1)*Exp_msk1 - y; + } + adj = aadj1 * ulp(dval(rv)); + dval(rv) += adj; #else - word0 (rv) = Big0; - word1 (rv) = Big1; -#endif - goto cont; - } - else - word0 (rv) += P * Exp_msk1; - } - else - { #ifdef Sudden_Underflow - if ((word0 (rv) & Exp_mask) <= P * Exp_msk1) - { - rv0.d = rv.d; - word0 (rv) += P * Exp_msk1; - adj = aadj1 * ulp (rv.d); - rv.d += adj; + if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) { + dval(rv0) = dval(rv); + dword0(rv) += P*Exp_msk1; + adj = aadj1 * ulp(dval(rv)); + dval(rv) += adj; #ifdef IBM - if ((word0 (rv) & Exp_mask) < P * Exp_msk1) + if ((dword0(rv) & Exp_mask) < P*Exp_msk1) #else - if ((word0 (rv) & Exp_mask) <= P * Exp_msk1) + if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) #endif - { - if (word0 (rv0) == Tiny0 - && word1 (rv0) == Tiny1) - goto undfl; - word0 (rv) = Tiny0; - word1 (rv) = Tiny1; - goto cont; + { + if (dword0(rv0) == Tiny0 + && dword1(rv0) == Tiny1) + goto undfl; + dword0(rv) = Tiny0; + dword1(rv) = Tiny1; + goto cont; + } + else + dword0(rv) -= P*Exp_msk1; + } + else { + adj = aadj1 * ulp(dval(rv)); + dval(rv) += adj; + } +#else /*Sudden_Underflow*/ + /* Compute adj so that the IEEE rounding rules will + * correctly round rv + adj in some half-way cases. + * If rv * ulp(rv) is denormalized (i.e., + * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid + * trouble from bits lost to denormalization; + * example: 1.2e-307 . + */ + if (y <= (P-1)*Exp_msk1 && aadj > 1.) { + aadj1 = (double)(int)(aadj + 0.5); + if (!dsign) + aadj1 = -aadj1; + } + adj = aadj1 * ulp(dval(rv)); + dval(rv) += adj; +#endif /*Sudden_Underflow*/ +#endif /*Avoid_Underflow*/ + } + z = dword0(rv) & Exp_mask; +#ifndef SET_INEXACT +#ifdef Avoid_Underflow + if (!scale) +#endif + if (y == z) { + /* Can we stop now? */ + L = (Long)aadj; + aadj -= L; + /* The tolerances below are conservative. */ + if (dsign || dword1(rv) || dword0(rv) & Bndry_mask) { + if (aadj < .4999999 || aadj > .5000001) + break; + } + else if (aadj < .4999999/FLT_RADIX) + break; + } +#endif + cont: + Bfree(ptr,bb); + Bfree(ptr,bd); + Bfree(ptr,bs); + Bfree(ptr,delta); } - else - word0 (rv) -= P * Exp_msk1; - } - else - { - adj = aadj1 * ulp (rv.d); - rv.d += adj; - } -#else - /* Compute adj so that the IEEE rounding rules will - * correctly round rv.d + adj in some half-way cases. - * If rv.d * ulp(rv.d) is denormalized (i.e., - * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid - * trouble from bits lost to denormalization; - * example: 1.2e-307 . - */ - if (y <= (P - 1) * Exp_msk1 && aadj >= 1.) - { - aadj1 = (double) (int) (aadj + 0.5); - if (!dsign) - aadj1 = -aadj1; - } - adj = aadj1 * ulp (rv.d); - rv.d += adj; +#ifdef SET_INEXACT + if (inexact) { + if (!oldinexact) { + dword0(rv0) = Exp_1 + (70 << Exp_shift); + dword1(rv0) = 0; + dval(rv0) += 1.; + } + } + else if (!oldinexact) + clear_inexact(); #endif - } - z = word0 (rv) & Exp_mask; - if (y == z) - { - /* Can we stop now? */ - L = aadj; - aadj -= L; - /* The tolerances below are conservative. */ - if (dsign || word1 (rv) || word0 (rv) & Bndry_mask) - { - if (aadj < .4999999 || aadj > .5000001) - break; - } - else if (aadj < .4999999 / FLT_RADIX) - break; - } - cont: - Bfree (ptr, bb); - Bfree (ptr, bd); - Bfree (ptr, bs); - Bfree (ptr, delta); - } -retfree: - Bfree (ptr, bb); - Bfree (ptr, bd); - Bfree (ptr, bs); - Bfree (ptr, bd0); - Bfree (ptr, delta); -ret: - if (se) - *se = (char *) s; - - if (nanflag) - return nan (NULL); - return (sign && (s != s00)) ? -rv.d : rv.d; +#ifdef Avoid_Underflow + if (scale) { + dword0(rv0) = Exp_1 - 2*P*Exp_msk1; + dword1(rv0) = 0; + dval(rv) *= dval(rv0); +#ifndef NO_ERRNO + /* try to avoid the bug of testing an 8087 register value */ + if (dword0(rv) == 0 && dword1(rv) == 0) + ptr->_errno = ERANGE; +#endif + } +#endif /* Avoid_Underflow */ +#ifdef SET_INEXACT + if (inexact && !(dword0(rv) & Exp_mask)) { + /* set underflow bit */ + dval(rv0) = 1e-300; + dval(rv0) *= dval(rv0); + } +#endif + retfree: + Bfree(ptr,bb); + Bfree(ptr,bd); + Bfree(ptr,bs); + Bfree(ptr,bd0); + Bfree(ptr,delta); + ret: + if (se) + *se = (char *)s; + return sign ? -dval(rv) : dval(rv); } #ifndef NO_REENT