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