Home | History | Annotate | Line # | Download | only in gdtoa
strtod.c revision 1.18.6.1 
      1 /* $NetBSD: strtod.c,v 1.18.6.1 2024/07/20 15:03:06 martin 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 	int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign,
    100 		 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
    101 	CONST char *s, *s0, *s1;
    102 	double aadj;
    103 	Long L;
    104 	U adj, aadj1, rv, rv0;
    105 	ULong y, z;
    106 	Bigint *bb = NULL, *bb1, *bd0;
    107 	Bigint *bd = NULL, *bs = NULL, *delta = NULL; /* pacify gcc */
    108 #ifdef Avoid_Underflow
    109 	ULong Lsb, Lsb1;
    110 #endif
    111 #ifdef SET_INEXACT
    112 	int inexact, oldinexact;
    113 #endif
    114 #ifdef USE_LOCALE /*{{*/
    115 	char *decimalpoint = localeconv_l(loc)->decimal_point;
    116 	size_t dplen = strlen(decimalpoint);
    117 #endif /*USE_LOCALE}}*/
    118 
    119 #ifdef Honor_FLT_ROUNDS /*{*/
    120 	int Rounding;
    121 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
    122 	Rounding = Flt_Rounds;
    123 #else /*}{*/
    124 	Rounding = 1;
    125 	switch(fegetround()) {
    126 	  case FE_TOWARDZERO:	Rounding = 0; break;
    127 	  case FE_UPWARD:	Rounding = 2; break;
    128 	  case FE_DOWNWARD:	Rounding = 3;
    129 	  }
    130 #endif /*}}*/
    131 #endif /*}*/
    132 
    133 	sign = nz0 = nz = decpt = 0;
    134 	dval(&rv) = 0.;
    135 	for(s = s00;;s++) switch(*s) {
    136 		case '-':
    137 			sign = 1;
    138 			/* FALLTHROUGH */
    139 		case '+':
    140 			if (*++s)
    141 				goto break2;
    142 			/* FALLTHROUGH */
    143 		case 0:
    144 			goto ret0;
    145 		case '\t':
    146 		case '\n':
    147 		case '\v':
    148 		case '\f':
    149 		case '\r':
    150 		case ' ':
    151 			continue;
    152 		default:
    153 			goto break2;
    154 		}
    155  break2:
    156 	if (*s == '0') {
    157 #ifndef NO_HEX_FP /*{*/
    158 		{
    159 		static CONST FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
    160 		Long expt;
    161 		ULong bits[2];
    162 		switch(s[1]) {
    163 		  case 'x':
    164 		  case 'X':
    165 			{
    166 #ifdef Honor_FLT_ROUNDS
    167 			FPI fpi1 = fpi;
    168 			fpi1.rounding = Rounding;
    169 #else
    170 #define fpi1 fpi
    171 #endif
    172 			switch((i = gethex(&s, &fpi1, &expt, &bb, sign, loc)) & STRTOG_Retmask) {
    173 			  case STRTOG_NoNumber:
    174 				s = s00;
    175 				sign = 0;
    176 				/* FALLTHROUGH */
    177 			  case STRTOG_Zero:
    178 				break;
    179 			  default:
    180 				if (bb) {
    181 					copybits(bits, fpi.nbits, bb);
    182 					Bfree(bb);
    183 					}
    184 				ULtod((/* LINTED */(U*)&rv)->L, bits, expt, i);
    185 			  }}
    186 			goto ret;
    187 		  }
    188 		}
    189 #endif /*}*/
    190 		nz0 = 1;
    191 		while(*++s == '0') ;
    192 		if (!*s)
    193 			goto ret;
    194 		}
    195 	s0 = s;
    196 	y = z = 0;
    197 	for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
    198 		if (nd < 9)
    199 			y = 10*y + c - '0';
    200 		else if (nd < DBL_DIG + 2)
    201 			z = 10*z + c - '0';
    202 	nd0 = nd;
    203 #ifdef USE_LOCALE
    204 	if (c == *decimalpoint) {
    205 		for(i = 1; decimalpoint[i]; ++i)
    206 			if (s[i] != decimalpoint[i])
    207 				goto dig_done;
    208 		s += i;
    209 		c = *s;
    210 #else
    211 	if (c == '.') {
    212 		c = *++s;
    213 #endif
    214 		decpt = 1;
    215 		if (!nd) {
    216 			for(; c == '0'; c = *++s)
    217 				nz++;
    218 			if (c > '0' && c <= '9') {
    219 				s0 = s;
    220 				nf += nz;
    221 				nz = 0;
    222 				goto have_dig;
    223 				}
    224 			goto dig_done;
    225 			}
    226 		for(; c >= '0' && c <= '9'; c = *++s) {
    227  have_dig:
    228 			nz++;
    229 			if (c -= '0') {
    230 				nf += nz;
    231 				for(i = 1; i < nz; i++)
    232 					if (nd++ < 9)
    233 						y *= 10;
    234 					else if (nd <= DBL_DIG + 2)
    235 						z *= 10;
    236 				if (nd++ < 9)
    237 					y = 10*y + c;
    238 				else if (nd <= DBL_DIG + 2)
    239 					z = 10*z + c;
    240 				nz = 0;
    241 				}
    242 			}
    243 		}/*}*/
    244  dig_done:
    245 	e = 0;
    246 	if (c == 'e' || c == 'E') {
    247 		if (!nd && !nz && !nz0) {
    248 			goto ret0;
    249 			}
    250 		s00 = s;
    251 		esign = 0;
    252 		switch(c = *++s) {
    253 			case '-':
    254 				esign = 1;
    255 				/* FALLTHROUGH */
    256 			case '+':
    257 				c = *++s;
    258 			}
    259 		if (c >= '0' && c <= '9') {
    260 			while(c == '0')
    261 				c = *++s;
    262 			if (c > '0' && c <= '9') {
    263 				L = c - '0';
    264 				s1 = s;
    265 				while((c = *++s) >= '0' && c <= '9')
    266 					L = 10*L + c - '0';
    267 				if (s - s1 > 8 || L > 19999)
    268 					/* Avoid confusion from exponents
    269 					 * so large that e might overflow.
    270 					 */
    271 					e = 19999; /* safe for 16 bit ints */
    272 				else
    273 					e = (int)L;
    274 				if (esign)
    275 					e = -e;
    276 				}
    277 			else
    278 				e = 0;
    279 			}
    280 		else
    281 			s = s00;
    282 		}
    283 	if (!nd) {
    284 		if (!nz && !nz0) {
    285 #ifdef INFNAN_CHECK
    286 			/* Check for Nan and Infinity */
    287 			ULong bits[2];
    288 			static CONST FPI fpinan = /* only 52 explicit bits */
    289 				{ 52, 1-1023-53+1, 2046-1023-53+1, 1, SI };
    290 			if (!decpt)
    291 			 switch(c) {
    292 			  case 'i':
    293 			  case 'I':
    294 				if (match(&s,"nf")) {
    295 					--s;
    296 					if (!match(&s,"inity"))
    297 						++s;
    298 					word0(&rv) = 0x7ff00000;
    299 					word1(&rv) = 0;
    300 					goto ret;
    301 					}
    302 				break;
    303 			  case 'n':
    304 			  case 'N':
    305 				if (match(&s, "an")) {
    306 #ifndef No_Hex_NaN
    307 					if (*s == '(' /*)*/
    308 					 && hexnan(&s, &fpinan, bits)
    309 							== STRTOG_NaNbits) {
    310 						word0(&rv) = 0x7ff00000 | bits[1];
    311 						word1(&rv) = bits[0];
    312 						}
    313 					else {
    314 #endif
    315 						word0(&rv) = NAN_WORD0;
    316 						word1(&rv) = NAN_WORD1;
    317 #ifndef No_Hex_NaN
    318 						}
    319 #endif
    320 					goto ret;
    321 					}
    322 			  }
    323 #endif /* INFNAN_CHECK */
    324  ret0:
    325 			s = s00;
    326 			sign = 0;
    327 			}
    328 		goto ret;
    329 		}
    330 	e1 = e -= nf;
    331 
    332 	/* Now we have nd0 digits, starting at s0, followed by a
    333 	 * decimal point, followed by nd-nd0 digits.  The number we're
    334 	 * after is the integer represented by those digits times
    335 	 * 10**e */
    336 
    337 	if (!nd0)
    338 		nd0 = nd;
    339 	k = nd < DBL_DIG + 2 ? nd : DBL_DIG + 2;
    340 	dval(&rv) = y;
    341 	if (k > 9) {
    342 #ifdef SET_INEXACT
    343 		if (k > DBL_DIG)
    344 			oldinexact = get_inexact();
    345 #endif
    346 		dval(&rv) = tens[k - 9] * dval(&rv) + z;
    347 		}
    348 	bd0 = 0;
    349 	if (nd <= DBL_DIG
    350 #ifndef RND_PRODQUOT
    351 #ifndef Honor_FLT_ROUNDS
    352 		&& Flt_Rounds == 1
    353 #endif
    354 #endif
    355 			) {
    356 		if (!e)
    357 			goto ret;
    358 #ifndef ROUND_BIASED_without_Round_Up
    359 		if (e > 0) {
    360 			if (e <= Ten_pmax) {
    361 #ifdef VAX
    362 				goto vax_ovfl_check;
    363 #else
    364 #ifdef Honor_FLT_ROUNDS
    365 				/* round correctly FLT_ROUNDS = 2 or 3 */
    366 				if (sign) {
    367 					rv.d = -rv.d;
    368 					sign = 0;
    369 					}
    370 #endif
    371 				/* rv = */ rounded_product(dval(&rv), tens[e]);
    372 				goto ret;
    373 #endif
    374 				}
    375 			i = DBL_DIG - nd;
    376 			if (e <= Ten_pmax + i) {
    377 				/* A fancier test would sometimes let us do
    378 				 * this for larger i values.
    379 				 */
    380 #ifdef Honor_FLT_ROUNDS
    381 				/* round correctly FLT_ROUNDS = 2 or 3 */
    382 				if (sign) {
    383 					rv.d = -rv.d;
    384 					sign = 0;
    385 					}
    386 #endif
    387 				e -= i;
    388 				dval(&rv) *= tens[i];
    389 #ifdef VAX
    390 				/* VAX exponent range is so narrow we must
    391 				 * worry about overflow here...
    392 				 */
    393  vax_ovfl_check:
    394 				word0(&rv) -= P*Exp_msk1;
    395 				/* rv = */ rounded_product(dval(&rv), tens[e]);
    396 				if ((word0(&rv) & Exp_mask)
    397 				 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
    398 					goto ovfl;
    399 				word0(&rv) += P*Exp_msk1;
    400 #else
    401 				/* rv = */ rounded_product(dval(&rv), tens[e]);
    402 #endif
    403 				goto ret;
    404 				}
    405 			}
    406 #ifndef Inaccurate_Divide
    407 		else if (e >= -Ten_pmax) {
    408 #ifdef Honor_FLT_ROUNDS
    409 			/* round correctly FLT_ROUNDS = 2 or 3 */
    410 			if (sign) {
    411 				rv.d = -rv.d;
    412 				sign = 0;
    413 				}
    414 #endif
    415 			/* rv = */ rounded_quotient(dval(&rv), tens[-e]);
    416 			goto ret;
    417 			}
    418 #endif
    419 #endif /* ROUND_BIASED_without_Round_Up */
    420 		}
    421 	e1 += nd - k;
    422 
    423 #ifdef IEEE_Arith
    424 #ifdef SET_INEXACT
    425 	inexact = 1;
    426 	if (k <= DBL_DIG)
    427 		oldinexact = get_inexact();
    428 #endif
    429 #ifdef Avoid_Underflow
    430 	scale = 0;
    431 #endif
    432 #ifdef Honor_FLT_ROUNDS
    433 	if (Rounding >= 2) {
    434 		if (sign)
    435 			Rounding = Rounding == 2 ? 0 : 2;
    436 		else
    437 			if (Rounding != 2)
    438 				Rounding = 0;
    439 		}
    440 #endif
    441 #endif /*IEEE_Arith*/
    442 
    443 	/* Get starting approximation = rv * 10**e1 */
    444 
    445 	if (e1 > 0) {
    446 		if ( (i = e1 & 15) !=0)
    447 			dval(&rv) *= tens[i];
    448 		if (e1 &= ~15) {
    449 			if (e1 > DBL_MAX_10_EXP) {
    450  ovfl:
    451 				/* Can't trust HUGE_VAL */
    452 #ifdef IEEE_Arith
    453 #ifdef Honor_FLT_ROUNDS
    454 				switch(Rounding) {
    455 				  case 0: /* toward 0 */
    456 				  case 3: /* toward -infinity */
    457 					word0(&rv) = Big0;
    458 					word1(&rv) = Big1;
    459 					break;
    460 				  default:
    461 					word0(&rv) = Exp_mask;
    462 					word1(&rv) = 0;
    463 				  }
    464 #else /*Honor_FLT_ROUNDS*/
    465 				word0(&rv) = Exp_mask;
    466 				word1(&rv) = 0;
    467 #endif /*Honor_FLT_ROUNDS*/
    468 #ifdef SET_INEXACT
    469 				/* set overflow bit */
    470 				dval(&rv0) = 1e300;
    471 				dval(&rv0) *= dval(&rv0);
    472 #endif
    473 #else /*IEEE_Arith*/
    474 				word0(&rv) = Big0;
    475 				word1(&rv) = Big1;
    476 #endif /*IEEE_Arith*/
    477  range_err:
    478 				if (bd0) {
    479 					Bfree(bb);
    480 					Bfree(bd);
    481 					Bfree(bs);
    482 					Bfree(bd0);
    483 					Bfree(delta);
    484 					}
    485 #ifndef NO_ERRNO
    486 				errno = ERANGE;
    487 #endif
    488 				goto ret;
    489 				}
    490 			e1 = (unsigned int)e1 >> 4;
    491 			for(j = 0; e1 > 1; j++, e1 = (unsigned int)e1 >> 1)
    492 				if (e1 & 1)
    493 					dval(&rv) *= bigtens[j];
    494 		/* The last multiplication could overflow. */
    495 			word0(&rv) -= P*Exp_msk1;
    496 			dval(&rv) *= bigtens[j];
    497 			if ((z = word0(&rv) & Exp_mask)
    498 			 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
    499 				goto ovfl;
    500 			if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
    501 				/* set to largest number */
    502 				/* (Can't trust DBL_MAX) */
    503 				word0(&rv) = Big0;
    504 				word1(&rv) = Big1;
    505 				}
    506 			else
    507 				word0(&rv) += P*Exp_msk1;
    508 			}
    509 		}
    510 	else if (e1 < 0) {
    511 		e1 = -e1;
    512 		if ( (i = e1 & 15) !=0)
    513 			dval(&rv) /= tens[i];
    514 		if (e1 >>= 4) {
    515 			if (e1 >= 1 << n_bigtens)
    516 				goto undfl;
    517 #ifdef Avoid_Underflow
    518 			if (e1 & Scale_Bit)
    519 				scale = 2*P;
    520 			for(j = 0; e1 > 0; j++, e1 = (unsigned int)e1 >> 1)
    521 				if (e1 & 1)
    522 					dval(&rv) *= tinytens[j];
    523 			if (scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
    524 						>> Exp_shift)) > 0) {
    525 				/* scaled rv is denormal; zap j low bits */
    526 				if (j >= 32) {
    527 					word1(&rv) = 0;
    528 					if (j >= 53)
    529 					 word0(&rv) = (P+2)*Exp_msk1;
    530 					else
    531 					 word0(&rv) &= 0xffffffffU << (j-32);
    532 					}
    533 				else
    534 					word1(&rv) &= 0xffffffffU << j;
    535 				}
    536 #else
    537 			for(j = 0; e1 > 1; j++, e1 = (unsigned int)e1 >> 1)
    538 				if (e1 & 1)
    539 					dval(&rv) *= tinytens[j];
    540 			/* The last multiplication could underflow. */
    541 			dval(&rv0) = dval(&rv);
    542 			dval(&rv) *= tinytens[j];
    543 			if (!dval(&rv)) {
    544 				dval(&rv) = 2.*dval(&rv0);
    545 				dval(&rv) *= tinytens[j];
    546 #endif
    547 				if (!dval(&rv)) {
    548  undfl:
    549 					dval(&rv) = 0.;
    550 #ifdef Honor_FLT_ROUNDS
    551 					if (Rounding == 2)
    552 						word1(&rv) = 1;
    553 #endif
    554 					goto range_err;
    555 					}
    556 #ifndef Avoid_Underflow
    557 				word0(&rv) = Tiny0;
    558 				word1(&rv) = Tiny1;
    559 				/* The refinement below will clean
    560 				 * this approximation up.
    561 				 */
    562 				}
    563 #endif
    564 			}
    565 		}
    566 
    567 	/* Now the hard part -- adjusting rv to the correct value.*/
    568 
    569 	/* Put digits into bd: true value = bd * 10^e */
    570 
    571 	bd0 = s2b(s0, nd0, nd, y, dplen);
    572 	if (bd0 == NULL)
    573 		goto ovfl;
    574 
    575 	for(;;) {
    576 		bd = Balloc(bd0->k);
    577 		if (bd == NULL)
    578 			goto ovfl;
    579 		Bcopy(bd, bd0);
    580 		bb = d2b(dval(&rv), &bbe, &bbbits);	/* rv = bb * 2^bbe */
    581 		if (bb == NULL)
    582 			goto ovfl;
    583 		bs = i2b(1);
    584 		if (bs == NULL)
    585 			goto ovfl;
    586 
    587 		if (e >= 0) {
    588 			bb2 = bb5 = 0;
    589 			bd2 = bd5 = e;
    590 			}
    591 		else {
    592 			bb2 = bb5 = -e;
    593 			bd2 = bd5 = 0;
    594 			}
    595 		if (bbe >= 0)
    596 			bb2 += bbe;
    597 		else
    598 			bd2 -= bbe;
    599 		bs2 = bb2;
    600 #ifdef Honor_FLT_ROUNDS
    601 		if (Rounding != 1)
    602 			bs2++;
    603 #endif
    604 #ifdef Avoid_Underflow
    605 		Lsb = LSB;
    606 		Lsb1 = 0;
    607 		j = bbe - scale;
    608 		i = j + bbbits - 1;	/* logb(rv) */
    609 		j = P + 1 - bbbits;
    610 		if (i < Emin) {	/* denormal */
    611 			i = Emin - i;
    612 			j -= i;
    613 			if (i < 32)
    614 				Lsb <<= i;
    615 			else
    616 				Lsb1 = Lsb << (i-32);
    617 			}
    618 #else /*Avoid_Underflow*/
    619 #ifdef Sudden_Underflow
    620 #ifdef IBM
    621 		j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
    622 #else
    623 		j = P + 1 - bbbits;
    624 #endif
    625 #else /*Sudden_Underflow*/
    626 		j = bbe;
    627 		i = j + bbbits - 1;	/* logb(&rv) */
    628 		if (i < Emin)	/* denormal */
    629 			j += P - Emin;
    630 		else
    631 			j = P + 1 - bbbits;
    632 #endif /*Sudden_Underflow*/
    633 #endif /*Avoid_Underflow*/
    634 		bb2 += j;
    635 		bd2 += j;
    636 #ifdef Avoid_Underflow
    637 		bd2 += scale;
    638 #endif
    639 		i = bb2 < bd2 ? bb2 : bd2;
    640 		if (i > bs2)
    641 			i = bs2;
    642 		if (i > 0) {
    643 			bb2 -= i;
    644 			bd2 -= i;
    645 			bs2 -= i;
    646 			}
    647 		if (bb5 > 0) {
    648 			bs = pow5mult(bs, bb5);
    649 			if (bs == NULL)
    650 				goto ovfl;
    651 			bb1 = mult(bs, bb);
    652 			if (bb1 == NULL)
    653 				goto ovfl;
    654 			Bfree(bb);
    655 			bb = bb1;
    656 			}
    657 		if (bb2 > 0) {
    658 			bb = lshift(bb, bb2);
    659 			if (bb == NULL)
    660 				goto ovfl;
    661 			}
    662 		if (bd5 > 0) {
    663 			bd = pow5mult(bd, bd5);
    664 			if (bd == NULL)
    665 				goto ovfl;
    666 			}
    667 		if (bd2 > 0) {
    668 			bd = lshift(bd, bd2);
    669 			if (bd == NULL)
    670 				goto ovfl;
    671 			}
    672 		if (bs2 > 0) {
    673 			bs = lshift(bs, bs2);
    674 			if (bs == NULL)
    675 				goto ovfl;
    676 			}
    677 		delta = diff(bb, bd);
    678 		if (delta == NULL)
    679 			goto ovfl;
    680 		dsign = delta->sign;
    681 		delta->sign = 0;
    682 		i = cmp(delta, bs);
    683 #ifdef Honor_FLT_ROUNDS
    684 		if (Rounding != 1) {
    685 			if (i < 0) {
    686 				/* Error is less than an ulp */
    687 				if (!delta->x[0] && delta->wds <= 1) {
    688 					/* exact */
    689 #ifdef SET_INEXACT
    690 					inexact = 0;
    691 #endif
    692 					break;
    693 					}
    694 				if (Rounding) {
    695 					if (dsign) {
    696 						dval(&adj) = 1.;
    697 						goto apply_adj;
    698 						}
    699 					}
    700 				else if (!dsign) {
    701 					dval(&adj) = -1.;
    702 					if (!word1(&rv)
    703 					 && !(word0(&rv) & Frac_mask)) {
    704 						y = word0(&rv) & Exp_mask;
    705 #ifdef Avoid_Underflow
    706 						if (!scale || y > 2*P*Exp_msk1)
    707 #else
    708 						if (y)
    709 #endif
    710 						  {
    711 						  delta = lshift(delta,Log2P);
    712 						  if (delta == NULL)
    713 							goto ovfl;
    714 						  if (cmp(delta, bs) <= 0)
    715 							dval(&adj) = -0.5;
    716 						  }
    717 						}
    718  apply_adj:
    719 #ifdef Avoid_Underflow
    720 					if (scale && (y = word0(&rv) & Exp_mask)
    721 						<= 2*P*Exp_msk1)
    722 					  word0(&adj) += (2*P+1)*Exp_msk1 - y;
    723 #else
    724 #ifdef Sudden_Underflow
    725 					if ((word0(&rv) & Exp_mask) <=
    726 							P*Exp_msk1) {
    727 						word0(&rv) += P*Exp_msk1;
    728 						dval(&rv) += adj*ulp(&rv);
    729 						word0(&rv) -= P*Exp_msk1;
    730 						}
    731 					else
    732 #endif /*Sudden_Underflow*/
    733 #endif /*Avoid_Underflow*/
    734 					dval(&rv) += adj.d*ulp(&rv);
    735 					}
    736 				break;
    737 				}
    738 			dval(&adj) = ratio(delta, bs);
    739 			if (adj.d < 1.)
    740 				dval(&adj) = 1.;
    741 			if (adj.d <= 0x7ffffffe) {
    742 				/* dval(&adj) = Rounding ? ceil(&adj) : floor(&adj); */
    743 				y = adj.d;
    744 				if (y != adj.d) {
    745 					if (!(((unsigned int)Rounding>>1) ^ (unsigned int)dsign))
    746 						y++;
    747 					dval(&adj) = y;
    748 					}
    749 				}
    750 #ifdef Avoid_Underflow
    751 			if (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
    752 				word0(&adj) += (2*P+1)*Exp_msk1 - y;
    753 #else
    754 #ifdef Sudden_Underflow
    755 			if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
    756 				word0(&rv) += P*Exp_msk1;
    757 				dval(&adj) *= ulp(&rv);
    758 				if (dsign)
    759 					dval(&rv) += adj;
    760 				else
    761 					dval(&rv) -= adj;
    762 				word0(&rv) -= P*Exp_msk1;
    763 				goto cont;
    764 				}
    765 #endif /*Sudden_Underflow*/
    766 #endif /*Avoid_Underflow*/
    767 			dval(&adj) *= ulp(&rv);
    768 			if (dsign) {
    769 				if (word0(&rv) == Big0 && word1(&rv) == Big1)
    770 					goto ovfl;
    771 				dval(&rv) += adj.d;
    772 				}
    773 			else
    774 				dval(&rv) -= adj.d;
    775 			goto cont;
    776 			}
    777 #endif /*Honor_FLT_ROUNDS*/
    778 
    779 		if (i < 0) {
    780 			/* Error is less than half an ulp -- check for
    781 			 * special case of mantissa a power of two.
    782 			 */
    783 			if (dsign || word1(&rv) || word0(&rv) & Bndry_mask
    784 #ifdef IEEE_Arith
    785 #ifdef Avoid_Underflow
    786 			 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
    787 #else
    788 			 || (word0(&rv) & Exp_mask) <= Exp_msk1
    789 #endif
    790 #endif
    791 				) {
    792 #ifdef SET_INEXACT
    793 				if (!delta->x[0] && delta->wds <= 1)
    794 					inexact = 0;
    795 #endif
    796 				break;
    797 				}
    798 			if (!delta->x[0] && delta->wds <= 1) {
    799 				/* exact result */
    800 #ifdef SET_INEXACT
    801 				inexact = 0;
    802 #endif
    803 				break;
    804 				}
    805 			delta = lshift(delta,Log2P);
    806 			if (delta == NULL)
    807 				goto ovfl;
    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