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