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