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