*/
/****************************************************************
- *
- * 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 <reent.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
+
+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),
_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
+#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, 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
+#define fpi1 fpi
+#endif
+ 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
- && FLT_ROUNDS == 1
-#endif
- )
- {
- if (!e)
- goto ret;
- if (e > 0)
- {
- if (e <= Ten_pmax)
- {
+#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;
+ goto vax_ovfl_check;
#else
- /* rv.d = */ rounded_product (rv.d, tens[e]);
- goto ret;
-#endif
- }
- i = DBL_DIG - nd;
- if (e <= Ten_pmax + i)
- {
- /* A fancier test would sometimes let us do
+#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;
- }
-#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 */
+ else if (e >= -Ten_pmax) {
+#ifdef Honor_FLT_ROUNDS
+ /* round correctly FLT_ROUNDS = 2 or 3 */
+ if (sign) {
+ rv = -rv;
+ sign = 0;
+ }
+#endif
+ /* rv = */ rounded_quotient(dval(rv), tens[-e]);
+ goto ret;
+ }
+#endif
+ }
+ e1 += nd - k;
+
#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;
+#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;
+ 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;
+ }
}
-#ifndef _DOUBLE_IS_32BITS
- word0 (rv) = Tiny0;
- word1 (rv) = Tiny1;
-#else
- word0 (rv) = Tiny1;
-#endif
- /* The refinement below will clean
- * this approximation up.
- */
- }
- }
- }
-
- /* 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, rv.d, &bbe, &bbbits); /* rv.d = 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 Sudden_Underflow
-#ifdef IBM
- j = 1 + 4 * P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
+ 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
- j = P + 1 - bbbits;
+ 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
-#else
- i = bbe + bbbits - 1; /* logb(rv.d) */
- if (i < Emin) /* denormal */
- j = bbe + (P - Emin);
- else
- j = P + 1 - bbbits;
-#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
+ if (!dval(rv)) {
+ undfl:
+ dval(rv) = 0.;
+#ifndef NO_ERRNO
+ ptr->_errno = ERANGE;
#endif
- ;
-#ifndef _DOUBLE_IS_32BITS
- word1 (rv) = 0;
+ if (bd0)
+ goto retfree;
+ goto ret;
+ }
+#ifndef Avoid_Underflow
+ dword0(rv) = Tiny0;
+ dword1(rv) = Tiny1;
+ /* The refinement below will clean
+ * this approximation up.
+ */
+ }
#endif
- break;
+ }
}
- }
- else if (!(word0 (rv) & Bndry_mask) && !word1 (rv))
- {
- drop_down:
- /* boundary case -- decrement exponent */
+
+ /* 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
- L = word0 (rv) & Exp_mask;
#ifdef IBM
- if (L < Exp_msk1)
+ j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
#else
- if (L <= Exp_msk1)
+ 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
- goto undfl;
- L -= Exp_msk1;
+ 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
- L = (word0 (rv) & Exp_mask) - Exp_msk1;
+ if (y)
#endif
- word0 (rv) = L | Bndry_mask1;
-#ifndef _DOUBLE_IS_32BITS
- word1 (rv) = 0xffffffff;
+ {
+ 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
+ ;
+ 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;
+ goto cont;
#else
- break;
+ break;
#endif
- }
+ }
#ifndef ROUND_BIASED
- if (!(word1 (rv) & LSB))
- break;
+ if (!(dword1(rv) & LSB))
+ break;
#endif
- if (dsign)
- rv.d += ulp (rv.d);
+ if (dsign)
+ dval(rv) += ulp(dval(rv));
#ifndef ROUND_BIASED
- else
- {
- rv.d -= ulp (rv.d);
+ else {
+ dval(rv) -= ulp(dval(rv));
#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)
- {
+ 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 (word1 (rv) == Tiny1 && !word0 (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;
+ 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;
-
- /* 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 (Flt_Rounds == 0)
+ 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 = 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;
-#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 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
+#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
git://sourceware.org / newlib-cygwin.git / commitdiff