[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 at acm dot org,
 * with " at " changed at "@" and " dot " changed to ".").	*/

#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, 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

	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':
			switch((i = gethex(&s, &fpi, &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(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;
					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 {
#endif
						word0(rv) = NAN_WORD0;
						word1(rv) = NAN_WORD1;
#ifndef No_Hex_NaN
						}
#endif
					goto ret;
					}
			  }
#endif /* INFNAN_CHECK */
 ret0:
			s = s00;
			sign = 0;
			}
		goto ret;
		}
	e1 = e -= nf;

	/* Now we have nd0 digits, starting at s0, followed by a
	 * decimal point, followed by nd-nd0 digits.  The number we're
	 * after is the integer represented by those digits times
	 * 10**e */

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