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