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