strtod.c

Newlib 嵌入式 C库 标准实现代码

			dval(rv) *= bigtens[j];
			if ((z = dword0(rv) & Exp_mask)
			 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
				goto ovfl;
			if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
				/* set to largest number */
				/* (Can't trust DBL_MAX) */
				dword0(rv) = Big0;
#ifndef _DOUBLE_IS_32BITS
				dword1(rv) = Big1;
#endif /*!_DOUBLE_IS_32BITS*/
				}
			else
				dword0(rv) += P*Exp_msk1;
			}
		}
	else if (e1 < 0) {
		e1 = -e1;
		if ( (i = e1 & 15) !=0)
			dval(rv) /= tens[i];
		if (e1 >>= 4) {
			if (e1 >= 1 << n_bigtens)
				goto undfl;
#ifdef Avoid_Underflow
			if (e1 & Scale_Bit)
				scale = 2*P;
			for(j = 0; e1 > 0; j++, e1 >>= 1)
				if (e1 & 1)
					dval(rv) *= tinytens[j];
			if (scale && (j = 2*P + 1 - ((dword0(rv) & Exp_mask)
						>> Exp_shift)) > 0) {
				/* scaled rv is denormal; zap j low bits */
				if (j >= 32) {
#ifndef _DOUBLE_IS_32BITS
					dword1(rv) = 0;
#endif /*!_DOUBLE_IS_32BITS*/
					if (j >= 53)
					 dword0(rv) = (P+2)*Exp_msk1;
					else
					 dword0(rv) &= 0xffffffff << (j-32);
					}
#ifndef _DOUBLE_IS_32BITS
				else
					dword1(rv) &= 0xffffffff << j;
#endif /*!_DOUBLE_IS_32BITS*/
				}
#else
			for(j = 0; e1 > 1; j++, e1 >>= 1)
				if (e1 & 1)
					dval(rv) *= tinytens[j];
			/* The last multiplication could underflow. */
			dval(rv0) = dval(rv);
			dval(rv) *= tinytens[j];
			if (!dval(rv)) {
				dval(rv) = 2.*dval(rv0);
				dval(rv) *= tinytens[j];
#endif
				if (!dval(rv)) {
 undfl:
					dval(rv) = 0.;
#ifndef NO_ERRNO
					ptr->_errno = ERANGE;
#endif
					if (bd0)
						goto retfree;
					goto ret;
					}
#ifndef Avoid_Underflow
#ifndef _DOUBLE_IS_32BITS
				dword0(rv) = Tiny0;
				dword1(rv) = Tiny1;
#else
				dword0(rv) = Tiny1;
#endif /*_DOUBLE_IS_32BITS*/
				/* The refinement below will clean
				 * this approximation up.
				 */
				}
#endif
			}
		}

	/* Now the hard part -- adjusting rv to the correct value.*/

	/* Put digits into bd: true value = bd * 10^e */

	bd0 = s2b(ptr, s0, nd0, nd, y);

	for(;;) {
		bd = Balloc(ptr,bd0->_k);
		Bcopy(bd, bd0);
		bb = d2b(ptr,dval(rv), &bbe, &bbbits);	/* rv = bb * 2^bbe */
		bs = i2b(ptr,1);

		if (e >= 0) {
			bb2 = bb5 = 0;
			bd2 = bd5 = e;
			}
		else {
			bb2 = bb5 = -e;
			bd2 = bd5 = 0;
			}
		if (bbe >= 0)
			bb2 += bbe;
		else
			bd2 -= bbe;
		bs2 = bb2;
#ifdef Honor_FLT_ROUNDS
		if (rounding != 1)
			bs2++;
#endif
#ifdef Avoid_Underflow
		j = bbe - scale;
		i = j + bbbits - 1;	/* logb(rv) */
		if (i < Emin)	/* denormal */
			j += P - Emin;
		else
			j = P + 1 - bbbits;
#else /*Avoid_Underflow*/
#ifdef Sudden_Underflow
#ifdef IBM
		j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
#else
		j = P + 1 - bbbits;
#endif
#else /*Sudden_Underflow*/
		j = bbe;
		i = j + bbbits - 1;	/* logb(rv) */
		if (i < Emin)	/* denormal */
			j += P - Emin;
		else
			j = P + 1 - bbbits;
#endif /*Sudden_Underflow*/
#endif /*Avoid_Underflow*/
		bb2 += j;
		bd2 += j;
#ifdef Avoid_Underflow
		bd2 += scale;
#endif
		i = bb2 < bd2 ? bb2 : bd2;
		if (i > bs2)
			i = bs2;
		if (i > 0) {
			bb2 -= i;
			bd2 -= i;
			bs2 -= i;
			}
		if (bb5 > 0) {
			bs = pow5mult(ptr, bs, bb5);
			bb1 = mult(ptr, bs, bb);
			Bfree(ptr, bb);
			bb = bb1;
			}
		if (bb2 > 0)
			bb = lshift(ptr, bb, bb2);
		if (bd5 > 0)
			bd = pow5mult(ptr, bd, bd5);
		if (bd2 > 0)
			bd = lshift(ptr, bd, bd2);
		if (bs2 > 0)
			bs = lshift(ptr, bs, bs2);
		delta = diff(ptr, bb, bd);
		dsign = delta->_sign;
		delta->_sign = 0;
		i = cmp(delta, bs);
#ifdef Honor_FLT_ROUNDS
		if (rounding != 1) {
			if (i < 0) {
				/* Error is less than an ulp */
				if (!delta->_x[0] && delta->_wds <= 1) {
					/* exact */
#ifdef SET_INEXACT
					inexact = 0;
#endif
					break;
					}
				if (rounding) {
					if (dsign) {
						adj = 1.;
						goto apply_adj;
						}
					}
				else if (!dsign) {
					adj = -1.;
					if (!dword1(rv)
					 && !(dword0(rv) & Frac_mask)) {
						y = dword0(rv) & Exp_mask;
#ifdef Avoid_Underflow
						if (!scale || y > 2*P*Exp_msk1)
#else
						if (y)
#endif
						  {
						  delta = lshift(ptr, delta,Log2P);
						  if (cmp(delta, bs) <= 0)
							adj = -0.5;
						  }
						}
 apply_adj:
#ifdef Avoid_Underflow
					if (scale && (y = dword0(rv) & Exp_mask)
						<= 2*P*Exp_msk1)
					  dword0(adj) += (2*P+1)*Exp_msk1 - y;
#else
#ifdef Sudden_Underflow
					if ((dword0(rv) & Exp_mask) <=
							P*Exp_msk1) {
						dword0(rv) += P*Exp_msk1;
						dval(rv) += adj*ulp(dval(rv));
						dword0(rv) -= P*Exp_msk1;
						}
					else
#endif /*Sudden_Underflow*/
#endif /*Avoid_Underflow*/
					dval(rv) += adj*ulp(dval(rv));
					}
				break;
				}
			adj = ratio(delta, bs);
			if (adj < 1.)
				adj = 1.;
			if (adj <= 0x7ffffffe) {
				/* adj = rounding ? ceil(adj) : floor(adj); */
				y = adj;
				if (y != adj) {
					if (!((rounding>>1) ^ dsign))
						y++;
					adj = y;
					}
				}
#ifdef Avoid_Underflow
			if (scale && (y = dword0(rv) & Exp_mask) <= 2*P*Exp_msk1)
				dword0(adj) += (2*P+1)*Exp_msk1 - y;
#else
#ifdef Sudden_Underflow
			if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) {
				dword0(rv) += P*Exp_msk1;
				adj *= ulp(dval(rv));
				if (dsign)
					dval(rv) += adj;
				else
					dval(rv) -= adj;
				dword0(rv) -= P*Exp_msk1;
				goto cont;
				}
#endif /*Sudden_Underflow*/
#endif /*Avoid_Underflow*/
			adj *= ulp(dval(rv));
			if (dsign)
				dval(rv) += adj;
			else
				dval(rv) -= adj;
			goto cont;
			}
#endif /*Honor_FLT_ROUNDS*/

		if (i < 0) {
			/* Error is less than half an ulp -- check for
			 * special case of mantissa a power of two.
			 */
			if (dsign || dword1(rv) || dword0(rv) & Bndry_mask
#ifdef IEEE_Arith
#ifdef Avoid_Underflow
			 || (dword0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
#else
			 || (dword0(rv) & Exp_mask) <= Exp_msk1
#endif
#endif
				) {
#ifdef SET_INEXACT
				if (!delta->x[0] && delta->wds <= 1)
					inexact = 0;
#endif
				break;
				}
			if (!delta->_x[0] && delta->_wds <= 1) {
				/* exact result */
#ifdef SET_INEXACT
				inexact = 0;
#endif
				break;
				}
			delta = lshift(ptr,delta,Log2P);
			if (cmp(delta, bs) > 0)
				goto drop_down;
			break;
			}
		if (i == 0) {
			/* exactly half-way between */
			if (dsign) {
				if ((dword0(rv) & Bndry_mask1) == Bndry_mask1
				 &&  dword1(rv) == (
#ifdef Avoid_Underflow
			(scale && (y = dword0(rv) & Exp_mask) <= 2*P*Exp_msk1)
		? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
#endif
						   0xffffffff)) {
					/*boundary case -- increment exponent*/
					dword0(rv) = (dword0(rv) & Exp_mask)
						+ Exp_msk1
#ifdef IBM
						| Exp_msk1 >> 4
#endif
						;
#ifndef _DOUBLE_IS_32BITS
					dword1(rv) = 0;
#endif /*!_DOUBLE_IS_32BITS*/
#ifdef Avoid_Underflow
					dsign = 0;
#endif
					break;
					}
				}
			else if (!(dword0(rv) & Bndry_mask) && !dword1(rv)) {
 drop_down:
				/* boundary case -- decrement exponent */
#ifdef Sudden_Underflow /*{{*/
				L = dword0(rv) & Exp_mask;
#ifdef IBM
				if (L <  Exp_msk1)
#else
#ifdef Avoid_Underflow
				if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
#else
				if (L <= Exp_msk1)
#endif /*Avoid_Underflow*/
#endif /*IBM*/
					goto undfl;
				L -= Exp_msk1;
#else /*Sudden_Underflow}{*/
#ifdef Avoid_Underflow
				if (scale) {
					L = dword0(rv) & Exp_mask;
					if (L <= (2*P+1)*Exp_msk1) {
						if (L > (P+2)*Exp_msk1)
							/* round even ==> */
							/* accept rv */
							break;
						/* rv = smallest denormal */
						goto undfl;
						}
					}
#endif /*Avoid_Underflow*/
				L = (dword0(rv) & Exp_mask) - Exp_msk1;
#endif /*Sudden_Underflow}*/
				dword0(rv) = L | Bndry_mask1;
#ifndef _DOUBLE_IS_32BITS
				dword1(rv) = 0xffffffff;
#endif /*!_DOUBLE_IS_32BITS*/
#ifdef IBM
				goto cont;
#else
				break;
#endif
				}
#ifndef ROUND_BIASED
			if (!(dword1(rv) & LSB))
				break;
#endif
			if (dsign)
				dval(rv) += ulp(dval(rv));
#ifndef ROUND_BIASED
			else {
				dval(rv) -= ulp(dval(rv));
#ifndef Sudden_Underflow
				if (!dval(rv))
					goto undfl;
#endif
				}
#ifdef Avoid_Underflow
			dsign = 1 - dsign;
#endif
#endif
			break;
			}
		if ((aadj = ratio(delta, bs)) <= 2.) {
			if (dsign)
				aadj = aadj1 = 1.;
			else if (dword1(rv) || dword0(rv) & Bndry_mask) {
#ifndef Sudden_Underflow
				if (dword1(rv) == Tiny1 && !dword0(rv))
					goto undfl;
#endif
				aadj = 1.;
				aadj1 = -1.;
				}
			else {
				/* special case -- power of FLT_RADIX to be */
				/* rounded down... */

				if (aadj < 2./FLT_RADIX)
					aadj = 1./FLT_RADIX;
				else
					aadj *= 0.5;
				aadj1 = -aadj;
				}
			}
		else {
			aadj *= 0.5;
			aadj1 = dsign ? aadj : -aadj;
#ifdef Check_FLT_ROUNDS
			switch(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 = 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;
#ifndef _DOUBLE_IS_32BITS
				dword1(rv) = Big1;
#endif /*!_DOUBLE_IS_32BITS*/
				goto cont;
				}
			else
				dword0(rv) += P*Exp_msk1;
			}
		else {
#ifdef Avoid_Underflow
			if (scale && y <= 2*P*Exp_msk1) {
				if (aadj <= 0x7fffffff) {
					if ((z = aadj) <= 0)
						z = 1;
					aadj = z;
					aadj1 = dsign ? aadj : -aadj;
					}
				dword0(aadj1) += (2*P+1)*Exp_msk1 - y;
				}
			adj = aadj1 * ulp(dval(rv));
			dval(rv) += adj;
#else
#ifdef Sudden_Underflow
			if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) {
				dval(rv0) = dval(rv);
				dword0(rv) += P*Exp_msk1;
				adj = aadj1 * ulp(dval(rv));
				dval(rv) += adj;
#ifdef IBM
				if ((dword0(rv) & Exp_mask) <  P*Exp_msk1)
#else
				if ((dword0(rv) & Exp_mask) <= P*Exp_msk1)
#endif
					{
					if (dword0(rv0) == Tiny0
					 && dword1(rv0) == Tiny1)
						goto undfl;
#ifndef _DOUBLE_IS_32BITS
					dword0(rv) = Tiny0;
					dword1(rv) = Tiny1;
#else
					dword0(rv) = Tiny1;
#endif /*_DOUBLE_IS_32BITS*/
					goto cont;
					}
				else
					dword0(rv) -= P*Exp_msk1;
				}
			else {
				adj = 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);
		}
#ifdef SET_INEXACT
	if (inexact) {
		if (!oldinexact) {
			dword0(rv0) = Exp_1 + (70 << Exp_shift);
#ifndef _DOUBLE_IS_32BITS
			dword1(rv0) = 0;
#endif /*!_DOUBLE_IS_32BITS*/
			dval(rv0) += 1.;
			}
		}
	else if (!oldinexact)
		clear_inexact();
#endif
#ifdef Avoid_Underflow
	if (scale) {
		dword0(rv0) = Exp_1 - 2*P*Exp_msk1;
#ifndef _DOUBLE_IS_32BITS
		dword1(rv0) = 0;
#endif /*!_DOUBLE_IS_32BITS*/
		dval(rv) *= dval(rv0);
#ifndef NO_ERRNO
		/* try to avoid the bug of testing an 8087 register value */
		if (dword0(rv) == 0 && dword1(rv) == 0)
			ptr->_errno = ERANGE;
#endif
		}
#endif /* Avoid_Underflow */
#ifdef SET_INEXACT
	if (inexact && !(dword0(rv) & Exp_mask)) {
		/* set underflow bit */
		dval(rv0) = 1e-300;
		dval(rv0) *= dval(rv0);
		}
#endif
 retfree:
	Bfree(ptr,bb);
	Bfree(ptr,bd);
	Bfree(ptr,bs);
	Bfree(ptr,bd0);
	Bfree(ptr,delta);
 ret:
	if (se)
		*se = (char *)s;
	return sign ? -dval(rv) : dval(rv);
}

#ifndef NO_REENT

double
_DEFUN (strtod, (s00, se),
	_CONST char *s00 _AND char **se)
{
  return _strtod_r (_REENT, s00, se);
}

float
_DEFUN (strtof, (s00, se),
	_CONST char *s00 _AND
	char **se)
{
  double retval = _strtod_r (_REENT, s00, se);
  if (isnan (retval))
    return nanf (NULL);
  return (float)retval;
}

#endif