[apple/libc.git] / gdtoa / FreeBSD / gdtoa-strtod.c

/****************************************************************

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@acm.org
 */

#include "gdtoaimp.h"

#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

 double
strtod
#ifdef KR_headers
	(s00, se) CONST char *s00; char **se;
#else
	(CONST char *s00, char **se)
#endif
{
#ifdef Avoid_Underflow
	int scale;
#endif
	int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, 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

	sign = nz0 = nz = 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':
			switch(i = gethex(&s, &fpi, &exp, &bb, sign)) {
			  case STRTOG_NoNumber:
				s = s00;
				sign = 0;
			  case STRTOG_Zero:
				break;
			  default:
				copybits(bits, fpi.nbits, bb);
				Bfree(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
		{
		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 };
			switch(c) {
			  case 'i':
			  case 'I':
				if (match(&s,"nf")) {
					--s;
					if (!match(&s,"inity"))
						++s;
					word0(rv) = 0x7ff00000;
					word1(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) {
						word0(rv) = 0x7ff00000 | bits[1];
						word1(rv) = bits[0];
						}
					else {
						word0(rv) = NAN_WORD0;
						word1(rv) = NAN_WORD1;
						}
#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.
				 */
#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 = */ rounded_product(dval(rv), tens[e]);
				if ((word0(rv) & Exp_mask)
				 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
					goto ovfl;
				word0(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) {
				rv = -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
				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 */
					word0(rv) = Big0;
					word1(rv) = Big1;
					break;
				  default:
					word0(rv) = Exp_mask;
					word1(rv) = 0;
				  }
#else /*Honor_FLT_ROUNDS*/
				word0(rv) = Exp_mask;
				word1(rv) = 0;
#endif /*Honor_FLT_ROUNDS*/
#ifdef SET_INEXACT
				/* set overflow bit */
				dval(rv0) = 1e300;
				dval(rv0) *= dval(rv0);
#endif
#else /*IEEE_Arith*/
				word0(rv) = Big0;
				word1(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. */
			word0(rv) -= P*Exp_msk1;
			dval(rv) *= 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;
				word1(rv) = Big1;
				}
			else
				word0(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 - ((word0(rv) & Exp_mask)
						>> Exp_shift)) > 0) {
				/* scaled rv is denormal; zap j low bits */
				if (j >= 32) {
					word1(rv) = 0;
					if (j >= 53)
					 word0(rv) = (P+2)*Exp_msk1;
					else
					 word0(rv) &= 0xffffffff << j-32;
					}
				else
					word1(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
					errno = ERANGE;
#endif
					if (bd0)
						goto retfree;
					goto ret;
					}
#ifndef Avoid_Underflow
				word0(rv) = Tiny0;
				word1(rv) = Tiny1;
				/* 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(s0, nd0, nd, y);

	for(;;) {
		bd = Balloc(bd0->k);
		Bcopy(bd, bd0);
		bb = d2b(dval(rv), &bbe, &bbbits);	/* rv = bb * 2^bbe */
		bs = i2b(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(bs, bb5);
			bb1 = mult(bs, bb);
			Bfree(bb);
			bb = bb1;
			}
		if (bb2 > 0)
			bb = lshift(bb, bb2);
		if (bd5 > 0)
			bd = pow5mult(bd, bd5);
		if (bd2 > 0)
			bd = lshift(bd, bd2);
		if (bs2 > 0)
			bs = lshift(bs, bs2);
		delta = diff(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 (!word1(rv)
					 && !(word0(rv) & Frac_mask)) {
						y = word0(rv) & Exp_mask;
#ifdef Avoid_Underflow
						if (!scale || y > 2*P*Exp_msk1)
#else
						if (y)
#endif
						  {
						  delta = lshift(delta,Log2P);
						  if (cmp(delta, bs) <= 0)
							adj = -0.5;
						  }
						}
 apply_adj:
#ifdef Avoid_Underflow
					if (scale && (y = word0(rv) & Exp_mask)
						<= 2*P*Exp_msk1)
					  word0(adj) += (2*P+1)*Exp_msk1 - y;
#else
#ifdef Sudden_Underflow
					if ((word0(rv) & Exp_mask) <=
							P*Exp_msk1) {
						word0(rv) += P*Exp_msk1;
						dval(rv) += adj*ulp(dval(rv));
						word0(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 = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
				word0(adj) += (2*P+1)*Exp_msk1 - y;
#else
#ifdef Sudden_Underflow
			if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
				word0(rv) += P*Exp_msk1;
				adj *= ulp(dval(rv));
				if (dsign)
					dval(rv) += adj;
				else
					dval(rv) -= adj;
				word0(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 || word1(rv) || word0(rv) & Bndry_mask
#ifdef IEEE_Arith
#ifdef Avoid_Underflow
			 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
#else
			 || (word0(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(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) == (
#ifdef Avoid_Underflow
			(scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
		? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
#endif
						   0xffffffff)) {
					/*boundary case -- increment exponent*/
					word0(rv) = (word0(rv) & Exp_mask)
						+ Exp_msk1
#ifdef IBM
						| Exp_msk1 >> 4
#endif
						;
					word1(rv) = 0;
#ifdef Avoid_Underflow
					dsign = 0;
#endif
					break;
					}
				}
			else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
 drop_down:
				/* boundary case -- decrement exponent */
#ifdef Sudden_Underflow /*{{*/
				L = word0(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 = word0(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 = (word0(rv) & Exp_mask) - Exp_msk1;
#endif /*Sudden_Underflow}*/
				word0(rv) = L | Bndry_mask1;
				word1(rv) = 0xffffffff;
#ifdef IBM
				goto cont;
#else
				break;
#endif
				}
#ifndef ROUND_BIASED
			if (!(word1(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 (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... */

				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(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 /*Check_FLT_ROUNDS*/
			}
		y = word0(rv) & Exp_mask;

		/* Check for overflow */

		if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
			dval(rv0) = dval(rv);
			word0(rv) -= P*Exp_msk1;
			adj = aadj1 * ulp(dval(rv));
			dval(rv) += adj;
			if ((word0(rv) & Exp_mask) >=
					Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
				if (word0(rv0) == Big0 && word1(rv0) == Big1)
					goto ovfl;
				word0(rv) = Big0;
				word1(rv) = Big1;
				goto cont;
				}
			else
				word0(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;
					}
				word0(aadj1) += (2*P+1)*Exp_msk1 - y;
				}
			adj = aadj1 * ulp(dval(rv));
			dval(rv) += adj;
#else
#ifdef Sudden_Underflow
			if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
				dval(rv0) = dval(rv);
				word0(rv) += P*Exp_msk1;
				adj = aadj1 * ulp(dval(rv));
				dval(rv) += adj;
#ifdef IBM
				if ((word0(rv) & Exp_mask) <  P*Exp_msk1)
#else
				if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
#endif
					{
					if (word0(rv0) == Tiny0
					 && word1(rv0) == Tiny1)
						goto undfl;
					word0(rv) = Tiny0;
					word1(rv) = Tiny1;
					goto cont;
					}
				else
					word0(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 = word0(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 || word1(rv) || word0(rv) & Bndry_mask) {
				if (aadj < .4999999 || aadj > .5000001)
					break;
				}
			else if (aadj < .4999999/FLT_RADIX)
				break;
			}
#endif
 cont:
		Bfree(bb);
		Bfree(bd);
		Bfree(bs);
		Bfree(delta);
		}
#ifdef SET_INEXACT
	if (inexact) {
		if (!oldinexact) {
			word0(rv0) = Exp_1 + (70 << Exp_shift);
			word1(rv0) = 0;
			dval(rv0) += 1.;
			}
		}
	else if (!oldinexact)
		clear_inexact();
#endif
#ifdef Avoid_Underflow
	if (scale) {
		word0(rv0) = Exp_1 - 2*P*Exp_msk1;
		word1(rv0) = 0;
		dval(rv) *= dval(rv0);
#ifndef NO_ERRNO
		/* try to avoid the bug of testing an 8087 register value */
		if (word0(rv) == 0 && word1(rv) == 0)
			errno = ERANGE;
#endif
		}
#endif /* Avoid_Underflow */
#ifdef SET_INEXACT
	if (inexact && !(word0(rv) & Exp_mask)) {
		/* set underflow bit */
		dval(rv0) = 1e-300;
		dval(rv0) *= dval(rv0);
		}
#endif
 retfree:
	Bfree(bb);
	Bfree(bd);
	Bfree(bs);
	Bfree(bd0);
	Bfree(delta);
 ret:
	if (se)
		*se = (char *)s;
	return sign ? -dval(rv) : dval(rv);
	}
Commit	Line	Data
9385eb3d A	1	/****************************************************************
	2
	3	The author of this software is David M. Gay.
	4
	5	Copyright (C) 1998-2001 by Lucent Technologies
	6	All Rights Reserved
	7
	8	Permission to use, copy, modify, and distribute this software and
	9	its documentation for any purpose and without fee is hereby
	10	granted, provided that the above copyright notice appear in all
	11	copies and that both that the copyright notice and this
	12	permission notice and warranty disclaimer appear in supporting
	13	documentation, and that the name of Lucent or any of its entities
	14	not be used in advertising or publicity pertaining to
	15	distribution of the software without specific, written prior
	16	permission.
	17
	18	LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
	19	INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
	20	IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
	21	SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
	22	WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
	23	IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
	24	ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
	25	THIS SOFTWARE.
	26
	27	****************************************************************/
	28
	29	/* Please send bug reports to
	30	David M. Gay
	31	dmg@acm.org
	32	*/
	33
	34	#include "gdtoaimp.h"
	35
	36	#ifdef USE_LOCALE
	37	#include "locale.h"
	38	#endif
	39
	40	#ifdef IEEE_Arith
	41	#ifndef NO_IEEE_Scale
	42	#define Avoid_Underflow
	43	#undef tinytens
	44	/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
	45	/* flag unnecessarily. It leads to a song and dance at the end of strtod. */
	46	static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
	47	9007199254740992.e-256
	48	};
	49	#endif
	50	#endif
	51
	52	#ifdef Honor_FLT_ROUNDS
	53	#define Rounding rounding
	54	#undef Check_FLT_ROUNDS
	55	#define Check_FLT_ROUNDS
	56	#else
	57	#define Rounding Flt_Rounds
	58	#endif
	59
	60	double
	61	strtod
	62	#ifdef KR_headers
	63	(s00, se) CONST char s00; char *se;
	64	#else
65	(CONST char s00, char *se)
66	#endif
67	{
68	#ifdef Avoid_Underflow
69	int scale;
70	#endif
71	int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
72	e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
73	CONST char s, s0, *s1;
74	double aadj, aadj1, adj, rv, rv0;
75	Long L;
76	ULong y, z;
77	Bigint bb, bb1, bd, bd0, bs, delta;
78	#ifdef SET_INEXACT
79	int inexact, oldinexact;
80	#endif
81	#ifdef Honor_FLT_ROUNDS
82	int rounding;
83	#endif
84
85	sign = nz0 = nz = 0;
86	dval(rv) = 0.;
87	for(s = s00;;s++) switch(*s) {
88	case '-':
89	sign = 1;
90	/* no break */
91	case '+':
92	if (*++s)
93	goto break2;
94	/* no break */
95	case 0:
96	goto ret0;
97	case '\t':
98	case '\n':
99	case '\v':
100	case '\f':
101	case '\r':
102	case ' ':
103	continue;
104	default:
105	goto break2;
106	}
107	break2:
108	if (*s == '0') {
109	#ifndef NO_HEX_FP
110	{
111	static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
112	Long exp;
113	ULong bits[2];
114	switch(s[1]) {
115	case 'x':
116	case 'X':
117	switch(i = gethex(&s, &fpi, &exp, &bb, sign)) {
118	case STRTOG_NoNumber:
119	s = s00;
120	sign = 0;
121	case STRTOG_Zero:
122	break;
123	default:
124	copybits(bits, fpi.nbits, bb);
125	Bfree(bb);
126	ULtod(((U*)&rv)->L, bits, exp, i);
127	}
128	goto ret;
129	}
130	}
131	#endif
132	nz0 = 1;
133	while(*++s == '0') ;
134	if (!*s)
135	goto ret;
136	}
137	s0 = s;
138	y = z = 0;
139	for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
140	if (nd < 9)
141	y = 10*y + c - '0';
142	else if (nd < 16)
143	z = 10*z + c - '0';
144	nd0 = nd;
145	#ifdef USE_LOCALE
146	if (c == *localeconv()->decimal_point)
147	#else
148	if (c == '.')
149	#endif
150	{
151	c = *++s;
152	if (!nd) {
153	for(; c == '0'; c = *++s)
154	nz++;
155	if (c > '0' && c <= '9') {
156	s0 = s;
157	nf += nz;
158	nz = 0;
159	goto have_dig;
160	}
161	goto dig_done;
162	}
163	for(; c >= '0' && c <= '9'; c = *++s) {
164	have_dig:
165	nz++;
166	if (c -= '0') {
167	nf += nz;
168	for(i = 1; i < nz; i++)
169	if (nd++ < 9)
170	y *= 10;
171	else if (nd <= DBL_DIG + 1)
172	z *= 10;
173	if (nd++ < 9)
174	y = 10*y + c;
175	else if (nd <= DBL_DIG + 1)
176	z = 10*z + c;
177	nz = 0;
178	}
179	}
180	}
181	dig_done:
182	e = 0;
183	if (c == 'e' \|\| c == 'E') {
184	if (!nd && !nz && !nz0) {
185	goto ret0;
186	}
187	s00 = s;
188	esign = 0;
189	switch(c = *++s) {
190	case '-':
191	esign = 1;
192	case '+':
193	c = *++s;
194	}
195	if (c >= '0' && c <= '9') {
196	while(c == '0')
197	c = *++s;
198	if (c > '0' && c <= '9') {
199	L = c - '0';
200	s1 = s;
201	while((c = *++s) >= '0' && c <= '9')
202	L = 10*L + c - '0';
203	if (s - s1 > 8 \|\| L > 19999)
204	/* Avoid confusion from exponents
205	* so large that e might overflow.
206	*/
207	e = 19999; /* safe for 16 bit ints */
208	else
209	e = (int)L;
210	if (esign)
211	e = -e;
212	}
213	else
214	e = 0;
215	}
216	else
217	s = s00;
218	}
219	if (!nd) {
220	if (!nz && !nz0) {
221	#ifdef INFNAN_CHECK
222	/* Check for Nan and Infinity */
223	ULong bits[2];
224	static FPI fpinan = /* only 52 explicit bits */
225	{ 52, 1-1023-53+1, 2046-1023-53+1, 1, SI };
226	switch(c) {
227	case 'i':
228	case 'I':
229	if (match(&s,"nf")) {
230	--s;
231	if (!match(&s,"inity"))
232	++s;
233	word0(rv) = 0x7ff00000;
234	word1(rv) = 0;
235	goto ret;
236	}
237	break;
238	case 'n':
239	case 'N':
240	if (match(&s, "an")) {
241	#ifndef No_Hex_NaN
242	if (s == '(' /)*/
243	&& hexnan(&s, &fpinan, bits)
244	== STRTOG_NaNbits) {
245	word0(rv) = 0x7ff00000 \| bits[1];
246	word1(rv) = bits[0];
247	}
248	else {
249	word0(rv) = NAN_WORD0;
250	word1(rv) = NAN_WORD1;
251	}
252	#endif
253	goto ret;
254	}
255	}
256	#endif /* INFNAN_CHECK */
257	ret0:
258	s = s00;
259	sign = 0;
260	}
261	goto ret;
262	}
263	e1 = e -= nf;
264
265	/* Now we have nd0 digits, starting at s0, followed by a
266	* decimal point, followed by nd-nd0 digits. The number we're
267	* after is the integer represented by those digits times
268	* 10*e /
269
270	if (!nd0)
271	nd0 = nd;
272	k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
273	dval(rv) = y;
274	if (k > 9) {
275	#ifdef SET_INEXACT
276	if (k > DBL_DIG)
277	oldinexact = get_inexact();
278	#endif
279	dval(rv) = tens[k - 9] * dval(rv) + z;
280	}
281	bd0 = 0;
282	if (nd <= DBL_DIG
283	#ifndef RND_PRODQUOT
284	#ifndef Honor_FLT_ROUNDS
285	&& Flt_Rounds == 1
286	#endif
287	#endif
288	) {
289	if (!e)
290	goto ret;
291	if (e > 0) {
292	if (e <= Ten_pmax) {
293	#ifdef VAX
294	goto vax_ovfl_check;
295	#else
296	#ifdef Honor_FLT_ROUNDS
297	/* round correctly FLT_ROUNDS = 2 or 3 */
298	if (sign) {
299	rv = -rv;
300	sign = 0;
301	}
302	#endif
303	/* rv = */ rounded_product(dval(rv), tens[e]);
304	goto ret;
305	#endif
306	}
307	i = DBL_DIG - nd;
308	if (e <= Ten_pmax + i) {
309	/* A fancier test would sometimes let us do
310	* this for larger i values.
311	*/
312	#ifdef Honor_FLT_ROUNDS
313	/* round correctly FLT_ROUNDS = 2 or 3 */
314	if (sign) {
315	rv = -rv;
316	sign = 0;
317	}
318	#endif
319	e -= i;
320	dval(rv) *= tens[i];
321	#ifdef VAX
322	/* VAX exponent range is so narrow we must
323	* worry about overflow here...
324	*/
325	vax_ovfl_check:
326	word0(rv) -= P*Exp_msk1;
327	/* rv = */ rounded_product(dval(rv), tens[e]);
328	if ((word0(rv) & Exp_mask)
329	> Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
330	goto ovfl;
331	word0(rv) += P*Exp_msk1;
332	#else
333	/* rv = */ rounded_product(dval(rv), tens[e]);
334	#endif
335	goto ret;
336	}
337	}
338	#ifndef Inaccurate_Divide
339	else if (e >= -Ten_pmax) {
340	#ifdef Honor_FLT_ROUNDS
341	/* round correctly FLT_ROUNDS = 2 or 3 */
342	if (sign) {
343	rv = -rv;
344	sign = 0;
345	}
346	#endif
347	/* rv = */ rounded_quotient(dval(rv), tens[-e]);
348	goto ret;
349	}
350	#endif
351	}
352	e1 += nd - k;
353
354	#ifdef IEEE_Arith
355	#ifdef SET_INEXACT
356	inexact = 1;
357	if (k <= DBL_DIG)
358	oldinexact = get_inexact();
359	#endif
360	#ifdef Avoid_Underflow
361	scale = 0;
362	#endif
363	#ifdef Honor_FLT_ROUNDS
364	if ((rounding = Flt_Rounds) >= 2) {
365	if (sign)
366	rounding = rounding == 2 ? 0 : 2;
367	else
368	if (rounding != 2)
369	rounding = 0;
370	}
371	#endif
372	#endif /IEEE_Arith/
373
374	/* Get starting approximation = rv * 10*e1 /
375
376	if (e1 > 0) {
377	if ( (i = e1 & 15) !=0)
378	dval(rv) *= tens[i];
379	if (e1 &= ~15) {
380	if (e1 > DBL_MAX_10_EXP) {
381	ovfl:
382	#ifndef NO_ERRNO
383	errno = ERANGE;
384	#endif
385	/* Can't trust HUGE_VAL */
386	#ifdef IEEE_Arith
387	#ifdef Honor_FLT_ROUNDS
388	switch(rounding) {
389	case 0: /* toward 0 */
390	case 3: /* toward -infinity */
391	word0(rv) = Big0;
392	word1(rv) = Big1;
393	break;
394	default:
395	word0(rv) = Exp_mask;
396	word1(rv) = 0;
397	}
398	#else /Honor_FLT_ROUNDS/
399	word0(rv) = Exp_mask;
400	word1(rv) = 0;
401	#endif /Honor_FLT_ROUNDS/
402	#ifdef SET_INEXACT
403	/* set overflow bit */
404	dval(rv0) = 1e300;
405	dval(rv0) *= dval(rv0);
406	#endif
407	#else /IEEE_Arith/
408	word0(rv) = Big0;
409	word1(rv) = Big1;
410	#endif /IEEE_Arith/
411	if (bd0)
412	goto retfree;
413	goto ret;
414	}
415	e1 >>= 4;
416	for(j = 0; e1 > 1; j++, e1 >>= 1)
417	if (e1 & 1)
418	dval(rv) *= bigtens[j];
419	/* The last multiplication could overflow. */
420	word0(rv) -= P*Exp_msk1;
421	dval(rv) *= bigtens[j];
422	if ((z = word0(rv) & Exp_mask)
423	> Exp_msk1*(DBL_MAX_EXP+Bias-P))
424	goto ovfl;
425	if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
426	/* set to largest number */
427	/* (Can't trust DBL_MAX) */
428	word0(rv) = Big0;
429	word1(rv) = Big1;
430	}
431	else
432	word0(rv) += P*Exp_msk1;
433	}
434	}
435	else if (e1 < 0) {
436	e1 = -e1;
437	if ( (i = e1 & 15) !=0)
438	dval(rv) /= tens[i];
439	if (e1 >>= 4) {
440	if (e1 >= 1 << n_bigtens)
441	goto undfl;
442	#ifdef Avoid_Underflow
443	if (e1 & Scale_Bit)
444	scale = 2*P;
445	for(j = 0; e1 > 0; j++, e1 >>= 1)
446	if (e1 & 1)
447	dval(rv) *= tinytens[j];
448	if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
449	>> Exp_shift)) > 0) {
450	/* scaled rv is denormal; zap j low bits */
451	if (j >= 32) {
452	word1(rv) = 0;
453	if (j >= 53)
454	word0(rv) = (P+2)*Exp_msk1;
455	else
456	word0(rv) &= 0xffffffff << j-32;
457	}
458	else
459	word1(rv) &= 0xffffffff << j;
460	}
461	#else
462	for(j = 0; e1 > 1; j++, e1 >>= 1)
463	if (e1 & 1)
464	dval(rv) *= tinytens[j];
465	/* The last multiplication could underflow. */
466	dval(rv0) = dval(rv);
467	dval(rv) *= tinytens[j];
468	if (!dval(rv)) {
469	dval(rv) = 2.*dval(rv0);
470	dval(rv) *= tinytens[j];
471	#endif
472	if (!dval(rv)) {
473	undfl:
474	dval(rv) = 0.;
475	#ifndef NO_ERRNO
476	errno = ERANGE;
477	#endif
478	if (bd0)
479	goto retfree;
480	goto ret;
481	}
482	#ifndef Avoid_Underflow
483	word0(rv) = Tiny0;
484	word1(rv) = Tiny1;
485	/* The refinement below will clean
486	* this approximation up.
487	*/
488	}
489	#endif
490	}
491	}
492
493	/* Now the hard part -- adjusting rv to the correct value.*/
494
495	/* Put digits into bd: true value = bd * 10^e */
496
497	bd0 = s2b(s0, nd0, nd, y);
498
499	for(;;) {
500	bd = Balloc(bd0->k);
501	Bcopy(bd, bd0);
502	bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
503	bs = i2b(1);
504
505	if (e >= 0) {
506	bb2 = bb5 = 0;
507	bd2 = bd5 = e;
508	}
509	else {
510	bb2 = bb5 = -e;
511	bd2 = bd5 = 0;
512	}
513	if (bbe >= 0)
514	bb2 += bbe;
515	else
516	bd2 -= bbe;
517	bs2 = bb2;
518	#ifdef Honor_FLT_ROUNDS
519	if (rounding != 1)
520	bs2++;
521	#endif
522	#ifdef Avoid_Underflow
523	j = bbe - scale;
524	i = j + bbbits - 1; /* logb(rv) */
525	if (i < Emin) /* denormal */
526	j += P - Emin;
527	else
528	j = P + 1 - bbbits;
529	#else /Avoid_Underflow/
530	#ifdef Sudden_Underflow
531	#ifdef IBM
532	j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
533	#else
534	j = P + 1 - bbbits;
535	#endif
536	#else /Sudden_Underflow/
537	j = bbe;
538	i = j + bbbits - 1; /* logb(rv) */
539	if (i < Emin) /* denormal */
540	j += P - Emin;
541	else
542	j = P + 1 - bbbits;
543	#endif /Sudden_Underflow/
544	#endif /Avoid_Underflow/
545	bb2 += j;
546	bd2 += j;
547	#ifdef Avoid_Underflow
548	bd2 += scale;
549	#endif
550	i = bb2 < bd2 ? bb2 : bd2;
551	if (i > bs2)
552	i = bs2;
553	if (i > 0) {
554	bb2 -= i;
555	bd2 -= i;
556	bs2 -= i;
557	}
558	if (bb5 > 0) {
559	bs = pow5mult(bs, bb5);
560	bb1 = mult(bs, bb);
561	Bfree(bb);
562	bb = bb1;
563	}
564	if (bb2 > 0)
565	bb = lshift(bb, bb2);
566	if (bd5 > 0)
567	bd = pow5mult(bd, bd5);
568	if (bd2 > 0)
569	bd = lshift(bd, bd2);
570	if (bs2 > 0)
571	bs = lshift(bs, bs2);
572	delta = diff(bb, bd);
573	dsign = delta->sign;
574	delta->sign = 0;
575	i = cmp(delta, bs);
576	#ifdef Honor_FLT_ROUNDS
577	if (rounding != 1) {
578	if (i < 0) {
579	/* Error is less than an ulp */
580	if (!delta->x[0] && delta->wds <= 1) {
581	/* exact */
582	#ifdef SET_INEXACT
583	inexact = 0;
584	#endif
585	break;
586	}
587	if (rounding) {
588	if (dsign) {
589	adj = 1.;
590	goto apply_adj;
591	}
592	}
593	else if (!dsign) {
594	adj = -1.;
595	if (!word1(rv)
596	&& !(word0(rv) & Frac_mask)) {
597	y = word0(rv) & Exp_mask;
598	#ifdef Avoid_Underflow
599	if (!scale \|\| y > 2PExp_msk1)
600	#else
601	if (y)
602	#endif
603	{
604	delta = lshift(delta,Log2P);
605	if (cmp(delta, bs) <= 0)
606	adj = -0.5;
607	}
608	}
609	apply_adj:
610	#ifdef Avoid_Underflow
611	if (scale && (y = word0(rv) & Exp_mask)
612	<= 2PExp_msk1)
613	word0(adj) += (2P+1)Exp_msk1 - y;
614	#else
615	#ifdef Sudden_Underflow
616	if ((word0(rv) & Exp_mask) <=
617	P*Exp_msk1) {
618	word0(rv) += P*Exp_msk1;
619	dval(rv) += adj*ulp(dval(rv));
620	word0(rv) -= P*Exp_msk1;
621	}
622	else
623	#endif /Sudden_Underflow/
624	#endif /Avoid_Underflow/
625	dval(rv) += adj*ulp(dval(rv));
626	}
627	break;
628	}
629	adj = ratio(delta, bs);
630	if (adj < 1.)
631	adj = 1.;
632	if (adj <= 0x7ffffffe) {
633	/* adj = rounding ? ceil(adj) : floor(adj); */
634	y = adj;
635	if (y != adj) {
636	if (!((rounding>>1) ^ dsign))
637	y++;
638	adj = y;
639	}
640	}
641	#ifdef Avoid_Underflow
642	if (scale && (y = word0(rv) & Exp_mask) <= 2PExp_msk1)
643	word0(adj) += (2P+1)Exp_msk1 - y;
644	#else
645	#ifdef Sudden_Underflow
646	if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
647	word0(rv) += P*Exp_msk1;
648	adj *= ulp(dval(rv));
649	if (dsign)
650	dval(rv) += adj;
651	else
652	dval(rv) -= adj;
653	word0(rv) -= P*Exp_msk1;
654	goto cont;
655	}
656	#endif /Sudden_Underflow/
657	#endif /Avoid_Underflow/
658	adj *= ulp(dval(rv));
659	if (dsign)
660	dval(rv) += adj;
661	else
662	dval(rv) -= adj;
663	goto cont;
664	}
665	#endif /Honor_FLT_ROUNDS/
666
667	if (i < 0) {
668	/* Error is less than half an ulp -- check for
669	* special case of mantissa a power of two.
670	*/
671	if (dsign \|\| word1(rv) \|\| word0(rv) & Bndry_mask
672	#ifdef IEEE_Arith
673	#ifdef Avoid_Underflow
674	\|\| (word0(rv) & Exp_mask) <= (2P+1)Exp_msk1
675	#else
676	\|\| (word0(rv) & Exp_mask) <= Exp_msk1
677	#endif
678	#endif
679	) {
680	#ifdef SET_INEXACT
681	if (!delta->x[0] && delta->wds <= 1)
682	inexact = 0;
683	#endif
684	break;
685	}
686	if (!delta->x[0] && delta->wds <= 1) {
687	/* exact result */
688	#ifdef SET_INEXACT
689	inexact = 0;
690	#endif
691	break;
692	}
693	delta = lshift(delta,Log2P);
694	if (cmp(delta, bs) > 0)
695	goto drop_down;
696	break;
697	}
698	if (i == 0) {
699	/* exactly half-way between */
700	if (dsign) {
701	if ((word0(rv) & Bndry_mask1) == Bndry_mask1
702	&& word1(rv) == (
703	#ifdef Avoid_Underflow
704	(scale && (y = word0(rv) & Exp_mask) <= 2PExp_msk1)
705	? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
706	#endif
707	0xffffffff)) {
708	/boundary case -- increment exponent/
709	word0(rv) = (word0(rv) & Exp_mask)
710	+ Exp_msk1
711	#ifdef IBM
712	\| Exp_msk1 >> 4
713	#endif
714	;
715	word1(rv) = 0;
716	#ifdef Avoid_Underflow
717	dsign = 0;
718	#endif
719	break;
720	}
721	}
722	else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
723	drop_down:
724	/* boundary case -- decrement exponent */
725	#ifdef Sudden_Underflow /{{/
726	L = word0(rv) & Exp_mask;
727	#ifdef IBM
728	if (L < Exp_msk1)
729	#else
730	#ifdef Avoid_Underflow
731	if (L <= (scale ? (2P+1)Exp_msk1 : Exp_msk1))
732	#else
733	if (L <= Exp_msk1)
734	#endif /Avoid_Underflow/
735	#endif /IBM/
736	goto undfl;
737	L -= Exp_msk1;
738	#else /Sudden_Underflow}{/
739	#ifdef Avoid_Underflow
740	if (scale) {
741	L = word0(rv) & Exp_mask;
742	if (L <= (2P+1)Exp_msk1) {
743	if (L > (P+2)*Exp_msk1)
744	/* round even ==> */
745	/* accept rv */
746	break;
747	/* rv = smallest denormal */
748	goto undfl;
749	}
750	}
751	#endif /Avoid_Underflow/
752	L = (word0(rv) & Exp_mask) - Exp_msk1;
753	#endif /Sudden_Underflow}/
754	word0(rv) = L \| Bndry_mask1;
755	word1(rv) = 0xffffffff;
756	#ifdef IBM
757	goto cont;
758	#else
759	break;
760	#endif
761	}
762	#ifndef ROUND_BIASED
763	if (!(word1(rv) & LSB))
764	break;
765	#endif
766	if (dsign)
767	dval(rv) += ulp(dval(rv));
768	#ifndef ROUND_BIASED
769	else {
770	dval(rv) -= ulp(dval(rv));
771	#ifndef Sudden_Underflow
772	if (!dval(rv))
773	goto undfl;
774	#endif
775	}
776	#ifdef Avoid_Underflow
777	dsign = 1 - dsign;
778	#endif
779	#endif
780	break;
781	}
782	if ((aadj = ratio(delta, bs)) <= 2.) {
783	if (dsign)
784	aadj = aadj1 = 1.;
785	else if (word1(rv) \|\| word0(rv) & Bndry_mask) {
786	#ifndef Sudden_Underflow
787	if (word1(rv) == Tiny1 && !word0(rv))
788	goto undfl;
789	#endif
790	aadj = 1.;
791	aadj1 = -1.;
792	}
793	else {
794	/* special case -- power of FLT_RADIX to be */
795	/* rounded down... */
796
797	if (aadj < 2./FLT_RADIX)
798	aadj = 1./FLT_RADIX;
799	else
800	aadj *= 0.5;
801	aadj1 = -aadj;
802	}
803	}
804	else {
805	aadj *= 0.5;
806	aadj1 = dsign ? aadj : -aadj;
807	#ifdef Check_FLT_ROUNDS
808	switch(Rounding) {
809	case 2: /* towards +infinity */
810	aadj1 -= 0.5;
811	break;
812	case 0: /* towards 0 */
813	case 3: /* towards -infinity */
814	aadj1 += 0.5;
815	}
816	#else
817	if (Flt_Rounds == 0)
818	aadj1 += 0.5;
819	#endif /Check_FLT_ROUNDS/
820	}
821	y = word0(rv) & Exp_mask;
822
823	/* Check for overflow */
824
825	if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
826	dval(rv0) = dval(rv);
827	word0(rv) -= P*Exp_msk1;
828	adj = aadj1 * ulp(dval(rv));
829	dval(rv) += adj;
830	if ((word0(rv) & Exp_mask) >=
831	Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
832	if (word0(rv0) == Big0 && word1(rv0) == Big1)
833	goto ovfl;
834	word0(rv) = Big0;
835	word1(rv) = Big1;
836	goto cont;
837	}
838	else
839	word0(rv) += P*Exp_msk1;
840	}
841	else {
842	#ifdef Avoid_Underflow
843	if (scale && y <= 2PExp_msk1) {
844	if (aadj <= 0x7fffffff) {
845	if ((z = aadj) <= 0)
846	z = 1;
847	aadj = z;
848	aadj1 = dsign ? aadj : -aadj;
849	}
850	word0(aadj1) += (2P+1)Exp_msk1 - y;
851	}
852	adj = aadj1 * ulp(dval(rv));
853	dval(rv) += adj;
854	#else
855	#ifdef Sudden_Underflow
856	if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
857	dval(rv0) = dval(rv);
858	word0(rv) += P*Exp_msk1;
859	adj = aadj1 * ulp(dval(rv));
860	dval(rv) += adj;
861	#ifdef IBM
862	if ((word0(rv) & Exp_mask) < P*Exp_msk1)
863	#else
864	if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
865	#endif
866	{
867	if (word0(rv0) == Tiny0
868	&& word1(rv0) == Tiny1)
869	goto undfl;
870	word0(rv) = Tiny0;
871	word1(rv) = Tiny1;
872	goto cont;
873	}
874	else
875	word0(rv) -= P*Exp_msk1;
876	}
877	else {
878	adj = aadj1 * ulp(dval(rv));
879	dval(rv) += adj;
880	}
881	#else /Sudden_Underflow/
882	/* Compute adj so that the IEEE rounding rules will
883	* correctly round rv + adj in some half-way cases.
884	* If rv * ulp(rv) is denormalized (i.e.,
885	* y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
886	* trouble from bits lost to denormalization;
887	* example: 1.2e-307 .
888	*/
889	if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
890	aadj1 = (double)(int)(aadj + 0.5);
891	if (!dsign)
892	aadj1 = -aadj1;
893	}
894	adj = aadj1 * ulp(dval(rv));
895	dval(rv) += adj;
896	#endif /Sudden_Underflow/
897	#endif /Avoid_Underflow/
898	}
899	z = word0(rv) & Exp_mask;
900	#ifndef SET_INEXACT
901	#ifdef Avoid_Underflow
902	if (!scale)
903	#endif
904	if (y == z) {
905	/* Can we stop now? */
906	L = (Long)aadj;
907	aadj -= L;
908	/* The tolerances below are conservative. */
909	if (dsign \|\| word1(rv) \|\| word0(rv) & Bndry_mask) {
910	if (aadj < .4999999 \|\| aadj > .5000001)
911	break;
912	}
913	else if (aadj < .4999999/FLT_RADIX)
914	break;
915	}
916	#endif
917	cont:
918	Bfree(bb);
919	Bfree(bd);
920	Bfree(bs);
921	Bfree(delta);
922	}
923	#ifdef SET_INEXACT
924	if (inexact) {
925	if (!oldinexact) {
926	word0(rv0) = Exp_1 + (70 << Exp_shift);
927	word1(rv0) = 0;
928	dval(rv0) += 1.;
929	}
930	}
931	else if (!oldinexact)
932	clear_inexact();
933	#endif
934	#ifdef Avoid_Underflow
935	if (scale) {
936	word0(rv0) = Exp_1 - 2PExp_msk1;
937	word1(rv0) = 0;
938	dval(rv) *= dval(rv0);
939	#ifndef NO_ERRNO
940	/* try to avoid the bug of testing an 8087 register value */
941	if (word0(rv) == 0 && word1(rv) == 0)
942	errno = ERANGE;
943	#endif
944	}
945	#endif /* Avoid_Underflow */
946	#ifdef SET_INEXACT
947	if (inexact && !(word0(rv) & Exp_mask)) {
948	/* set underflow bit */
949	dval(rv0) = 1e-300;
950	dval(rv0) *= dval(rv0);
951	}
952	#endif
953	retfree:
954	Bfree(bb);
955	Bfree(bd);
956	Bfree(bs);
957	Bfree(bd0);
958	Bfree(delta);
959	ret:
960	if (se)
961	se = (char )s;
962	return sign ? -dval(rv) : dval(rv);
963	}
964