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