summaryrefslogtreecommitdiffstats
path: root/newlib/libc/stdlib/strtod.c
diff options
context:
space:
mode:
Diffstat (limited to 'newlib/libc/stdlib/strtod.c')
-rw-r--r--newlib/libc/stdlib/strtod.c1191
1 files changed, 0 insertions, 1191 deletions
diff --git a/newlib/libc/stdlib/strtod.c b/newlib/libc/stdlib/strtod.c
deleted file mode 100644
index 703fbce1d..000000000
--- a/newlib/libc/stdlib/strtod.c
+++ /dev/null
@@ -1,1191 +0,0 @@
-/*
-FUNCTION
- <<strtod>>, <<strtof>>---string to double or float
-
-INDEX
- strtod
-INDEX
- _strtod_r
-INDEX
- strtof
-
-ANSI_SYNOPSIS
- #include <stdlib.h>
- double strtod(const char *<[str]>, char **<[tail]>);
- float strtof(const char *<[str]>, char **<[tail]>);
-
- double _strtod_r(void *<[reent]>,
- const char *<[str]>, char **<[tail]>);
-
-TRAD_SYNOPSIS
- #include <stdlib.h>
- double strtod(<[str]>,<[tail]>)
- char *<[str]>;
- char **<[tail]>;
-
- float strtof(<[str]>,<[tail]>)
- char *<[str]>;
- char **<[tail]>;
-
- double _strtod_r(<[reent]>,<[str]>,<[tail]>)
- char *<[reent]>;
- char *<[str]>;
- char **<[tail]>;
-
-DESCRIPTION
- The function <<strtod>> parses the character string <[str]>,
- producing a substring which can be converted to a double
- value. The substring converted is the longest initial
- subsequence of <[str]>, beginning with the first
- non-whitespace character, that has one of these formats:
- .[+|-]<[digits]>[.[<[digits]>]][(e|E)[+|-]<[digits]>]
- .[+|-].<[digits]>[(e|E)[+|-]<[digits]>]
- .[+|-](i|I)(n|N)(f|F)[(i|I)(n|N)(i|I)(t|T)(y|Y)]
- .[+|-](n|N)(a|A)(n|N)[<(>[<[hexdigits]>]<)>]
- .[+|-]0(x|X)<[hexdigits]>[.[<[hexdigits]>]][(p|P)[+|-]<[digits]>]
- .[+|-]0(x|X).<[hexdigits]>[(p|P)[+|-]<[digits]>]
- The substring contains no characters if <[str]> is empty, consists
- entirely of whitespace, or if the first non-whitespace
- character is something other than <<+>>, <<->>, <<.>>, or a
- digit, and cannot be parsed as infinity or NaN. If the platform
- does not support NaN, then NaN is treated as an empty substring.
- If the substring is empty, no conversion is done, and
- the value of <[str]> is stored in <<*<[tail]>>>. Otherwise,
- the substring is converted, and a pointer to the final string
- (which will contain at least the terminating null character of
- <[str]>) is stored in <<*<[tail]>>>. If you want no
- assignment to <<*<[tail]>>>, pass a null pointer as <[tail]>.
- <<strtof>> is identical to <<strtod>> except for its return type.
-
- This implementation returns the nearest machine number to the
- input decimal string. Ties are broken by using the IEEE
- round-even rule. However, <<strtof>> is currently subject to
- double rounding errors.
-
- The alternate function <<_strtod_r>> is a reentrant version.
- The extra argument <[reent]> is a pointer to a reentrancy structure.
-
-RETURNS
- <<strtod>> returns the converted substring value, if any. If
- no conversion could be performed, 0 is returned. If the
- correct value is out of the range of representable values,
- plus or minus <<HUGE_VAL>> is returned, and <<ERANGE>> is
- stored in errno. If the correct value would cause underflow, 0
- is returned and <<ERANGE>> is stored in errno.
-
-Supporting OS subroutines required: <<close>>, <<fstat>>, <<isatty>>,
-<<lseek>>, <<read>>, <<sbrk>>, <<write>>.
-*/
-
-/****************************************************************
-
-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 <errno.h>
-#include <string.h>
-#include "mprec.h"
-#include "gdtoa.h"
-#include "gd_qnan.h"
-
-/* #ifndef NO_FENV_H */
-/* #include <fenv.h> */
-/* #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
-
-#ifndef NO_HEX_FP
-
-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;
-}
-#endif /* !NO_HEX_FP */
-
-#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),
- struct _reent *ptr _AND
- _CONST char *s00 _AND
- char **se)
-{
-#ifdef Avoid_Underflow
- int scale;
-#endif
- 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, adj;
- U aadj1, 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 the number is not hex, then the parse of
- 0 is still valid. */
- s00 = s + 1;
- {
-#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
-#define fpi1 fpi
-#endif
- switch((i = gethex(ptr, &s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) {
- case STRTOG_NoNumber:
- s = s00;
- case STRTOG_Zero:
- break;
- default:
- if (bb) {
- copybits(bits, fpi.nbits, bb);
- Bfree(ptr,bb);
- }
- ULtod(rv.i, 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;
-#ifndef _DOUBLE_IS_32BITS
- dword1(rv) = 0;
-#endif /*!_DOUBLE_IS_32BITS*/
- 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];
-#ifndef _DOUBLE_IS_32BITS
- dword1(rv) = bits[0];
-#endif /*!_DOUBLE_IS_32BITS*/
- }
- else {
-#endif
- dword0(rv) = NAN_WORD0;
-#ifndef _DOUBLE_IS_32BITS
- dword1(rv) = NAN_WORD1;
-#endif /*!_DOUBLE_IS_32BITS*/
-#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) {
- dval(rv) = -dval(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.
- */
-#ifdef Honor_FLT_ROUNDS
- /* round correctly FLT_ROUNDS = 2 or 3 */
- if (sign) {
- dval(rv) = -dval(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:
- 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 = */ rounded_product(dval(rv), tens[e]);
-#endif
- goto ret;
- }
- }
-#ifndef Inaccurate_Divide
- else if (e >= -Ten_pmax) {
-#ifdef Honor_FLT_ROUNDS
- /* round correctly FLT_ROUNDS = 2 or 3 */
- if (sign) {
- dval(rv) = -dval(rv);
- sign = 0;
- }
-#endif
- /* rv = */ rounded_quotient(dval(rv), tens[-e]);
- goto ret;
- }
-#endif
- }
- e1 += nd - k;
-
-#ifdef IEEE_Arith
-#ifdef SET_INEXACT
- inexact = 1;
- if (k <= DBL_DIG)
- oldinexact = get_inexact();
-#endif
-#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;
- }
-#endif
-#endif /*IEEE_Arith*/
-
- /* Get starting approximation = rv * 10**e1 */
-
- 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;
-#ifndef _DOUBLE_IS_32BITS
- dword1(rv) = Big1;
-#endif /*!_DOUBLE_IS_32BITS*/
- break;
- default:
- dword0(rv) = Exp_mask;
-#ifndef _DOUBLE_IS_32BITS
- dword1(rv) = 0;
-#endif /*!_DOUBLE_IS_32BITS*/
- }
-#else /*Honor_FLT_ROUNDS*/
- dword0(rv) = Exp_mask;
-#ifndef _DOUBLE_IS_32BITS
- dword1(rv) = 0;
-#endif /*!_DOUBLE_IS_32BITS*/
-#endif /*Honor_FLT_ROUNDS*/
-#ifdef SET_INEXACT
- /* set overflow bit */
- dval(rv0) = 1e300;
- dval(rv0) *= dval(rv0);
-#endif
-#else /*IEEE_Arith*/
- dword0(rv) = Big0;
-#ifndef _DOUBLE_IS_32BITS
- dword1(rv) = Big1;
-#endif /*!_DOUBLE_IS_32BITS*/
-#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;
-#ifndef _DOUBLE_IS_32BITS
- dword1(rv) = Big1;
-#endif /*!_DOUBLE_IS_32BITS*/
- }
- 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) {
-#ifndef _DOUBLE_IS_32BITS
- dword1(rv) = 0;
-#endif /*!_DOUBLE_IS_32BITS*/
- if (j >= 53)
- dword0(rv) = (P+2)*Exp_msk1;
- else
- dword0(rv) &= 0xffffffff << (j-32);
- }
-#ifndef _DOUBLE_IS_32BITS
- else
- dword1(rv) &= 0xffffffff << j;
-#endif /*!_DOUBLE_IS_32BITS*/
- }
-#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
-#ifndef _DOUBLE_IS_32BITS
- dword0(rv) = Tiny0;
- dword1(rv) = Tiny1;
-#else
- dword0(rv) = Tiny1;
-#endif /*_DOUBLE_IS_32BITS*/
- /* The refinement below will clean
- * this approximation up.
- */
- }
-#endif
- }
- }
-
- /* Now the hard part -- adjusting rv to the correct value.*/
-
- /* Put digits into bd: true value = bd * 10^e */
-
- 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);
-#else
- 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
- if (y)
-#endif
- {
- 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
- 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
-#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 (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
- ;
-#ifndef _DOUBLE_IS_32BITS
- dword1(rv) = 0;
-#endif /*!_DOUBLE_IS_32BITS*/
-#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;
-#ifndef _DOUBLE_IS_32BITS
- dword1(rv) = 0xffffffff;
-#endif /*!_DOUBLE_IS_32BITS*/
-#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 = dval(aadj1) = 1.;
- else if (dword1(rv) || dword0(rv) & Bndry_mask) {
-#ifndef Sudden_Underflow
- if (dword1(rv) == Tiny1 && !dword0(rv))
- goto undfl;
-#endif
- aadj = 1.;
- dval(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;
- dval(aadj1) = -aadj;
- }
- }
- else {
- aadj *= 0.5;
- dval(aadj1) = dsign ? aadj : -aadj;
-#ifdef Check_FLT_ROUNDS
- switch(Rounding) {
- case 2: /* towards +infinity */
- dval(aadj1) -= 0.5;
- break;
- case 0: /* towards 0 */
- case 3: /* towards -infinity */
- dval(aadj1) += 0.5;
- }
-#else
- if (Flt_Rounds == 0)
- dval(aadj1) += 0.5;
-#endif /*Check_FLT_ROUNDS*/
- }
- y = dword0(rv) & Exp_mask;
-
- /* Check for overflow */
-
- if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
- dval(rv0) = dval(rv);
- dword0(rv) -= P*Exp_msk1;
- adj = dval(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;
-#ifndef _DOUBLE_IS_32BITS
- dword1(rv) = Big1;
-#endif /*!_DOUBLE_IS_32BITS*/
- 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;
- dval(aadj1) = dsign ? aadj : -aadj;
- }
- dword0(aadj1) += (2*P+1)*Exp_msk1 - y;
- }
- adj = dval(aadj1) * ulp(dval(rv));
- dval(rv) += adj;
-#else
-#ifdef Sudden_Underflow
- if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) {
- dval(rv0) = dval(rv);
- dword0(rv) += P*Exp_msk1;
- adj = dval(aadj1) * ulp(dval(rv));
- dval(rv) += adj;
-#ifdef IBM
- if ((dword0(rv) & Exp_mask) < P*Exp_msk1)
-#else
- if ((dword0(rv) & Exp_mask) <= P*Exp_msk1)
-#endif
- {
- if (dword0(rv0) == Tiny0
- && dword1(rv0) == Tiny1)
- goto undfl;
-#ifndef _DOUBLE_IS_32BITS
- dword0(rv) = Tiny0;
- dword1(rv) = Tiny1;
-#else
- dword0(rv) = Tiny1;
-#endif /*_DOUBLE_IS_32BITS*/
- goto cont;
- }
- else
- dword0(rv) -= P*Exp_msk1;
- }
- else {
- adj = dval(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.) {
- dval(aadj1) = (double)(int)(aadj + 0.5);
- if (!dsign)
- dval(aadj1) = -dval(aadj1);
- }
- adj = dval(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);
- }
-#ifdef SET_INEXACT
- if (inexact) {
- if (!oldinexact) {
- dword0(rv0) = Exp_1 + (70 << Exp_shift);
-#ifndef _DOUBLE_IS_32BITS
- dword1(rv0) = 0;
-#endif /*!_DOUBLE_IS_32BITS*/
- dval(rv0) += 1.;
- }
- }
- else if (!oldinexact)
- clear_inexact();
-#endif
-#ifdef Avoid_Underflow
- if (scale) {
- dword0(rv0) = Exp_1 - 2*P*Exp_msk1;
-#ifndef _DOUBLE_IS_32BITS
- dword1(rv0) = 0;
-#endif /*!_DOUBLE_IS_32BITS*/
- 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
-
-double
-_DEFUN (strtod, (s00, se),
- _CONST char *s00 _AND char **se)
-{
- return _strtod_r (_REENT, s00, se);
-}
-
-float
-_DEFUN (strtof, (s00, se),
- _CONST char *s00 _AND
- char **se)
-{
- double retval = _strtod_r (_REENT, s00, se);
- if (isnan (retval))
- return nanf (NULL);
- return (float)retval;
-}
-
-#endif