xref: /plan9/sys/src/ape/lib/ap/stdio/strtod.c (revision 43aadf5e5598b9159a98f64e309c3ae860328a56)
13e12c5d1SDavid du Colombier #include "fconv.h"
23e12c5d1SDavid du Colombier 
33e12c5d1SDavid du Colombier /* strtod for IEEE-, VAX-, and IBM-arithmetic machines (dmg).
43e12c5d1SDavid du Colombier  *
53e12c5d1SDavid du Colombier  * This strtod returns a nearest machine number to the input decimal
63e12c5d1SDavid du Colombier  * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
73e12c5d1SDavid du Colombier  * broken by the IEEE round-even rule.  Otherwise ties are broken by
83e12c5d1SDavid du Colombier  * biased rounding (add half and chop).
93e12c5d1SDavid du Colombier  *
103e12c5d1SDavid du Colombier  * Inspired loosely by William D. Clinger's paper "How to Read Floating
113e12c5d1SDavid du Colombier  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
123e12c5d1SDavid du Colombier  *
133e12c5d1SDavid du Colombier  * Modifications:
143e12c5d1SDavid du Colombier  *
153e12c5d1SDavid du Colombier  *	1. We only require IEEE, IBM, or VAX double-precision
163e12c5d1SDavid du Colombier  *		arithmetic (not IEEE double-extended).
173e12c5d1SDavid du Colombier  *	2. We get by with floating-point arithmetic in a case that
183e12c5d1SDavid du Colombier  *		Clinger missed -- when we're computing d * 10^n
193e12c5d1SDavid du Colombier  *		for a small integer d and the integer n is not too
203e12c5d1SDavid du Colombier  *		much larger than 22 (the maximum integer k for which
213e12c5d1SDavid du Colombier  *		we can represent 10^k exactly), we may be able to
223e12c5d1SDavid du Colombier  *		compute (d*10^k) * 10^(e-k) with just one roundoff.
233e12c5d1SDavid du Colombier  *	3. Rather than a bit-at-a-time adjustment of the binary
243e12c5d1SDavid du Colombier  *		result in the hard case, we use floating-point
253e12c5d1SDavid du Colombier  *		arithmetic to determine the adjustment to within
263e12c5d1SDavid du Colombier  *		one bit; only in really hard cases do we need to
273e12c5d1SDavid du Colombier  *		compute a second residual.
283e12c5d1SDavid du Colombier  *	4. Because of 3., we don't need a large table of powers of 10
293e12c5d1SDavid du Colombier  *		for ten-to-e (just some small tables, e.g. of 10^k
303e12c5d1SDavid du Colombier  *		for 0 <= k <= 22).
313e12c5d1SDavid du Colombier  */
323e12c5d1SDavid du Colombier 
333e12c5d1SDavid du Colombier #ifdef RND_PRODQUOT
343e12c5d1SDavid du Colombier #define rounded_product(a,b) a = rnd_prod(a, b)
353e12c5d1SDavid du Colombier #define rounded_quotient(a,b) a = rnd_quot(a, b)
363e12c5d1SDavid du Colombier extern double rnd_prod(double, double), rnd_quot(double, double);
373e12c5d1SDavid du Colombier #else
383e12c5d1SDavid du Colombier #define rounded_product(a,b) a *= b
393e12c5d1SDavid du Colombier #define rounded_quotient(a,b) a /= b
403e12c5d1SDavid du Colombier #endif
413e12c5d1SDavid du Colombier 
423e12c5d1SDavid du Colombier  static double
ulp(double xarg)433e12c5d1SDavid du Colombier ulp(double xarg)
443e12c5d1SDavid du Colombier {
453e12c5d1SDavid du Colombier 	register long L;
463e12c5d1SDavid du Colombier 	Dul a;
473e12c5d1SDavid du Colombier 	Dul x;
483e12c5d1SDavid du Colombier 
493e12c5d1SDavid du Colombier 	x.d = xarg;
503e12c5d1SDavid du Colombier 	L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
513e12c5d1SDavid du Colombier #ifndef Sudden_Underflow
523e12c5d1SDavid du Colombier 	if (L > 0) {
533e12c5d1SDavid du Colombier #endif
543e12c5d1SDavid du Colombier #ifdef IBM
553e12c5d1SDavid du Colombier 		L |= Exp_msk1 >> 4;
563e12c5d1SDavid du Colombier #endif
573e12c5d1SDavid du Colombier 		word0(a) = L;
583e12c5d1SDavid du Colombier 		word1(a) = 0;
593e12c5d1SDavid du Colombier #ifndef Sudden_Underflow
603e12c5d1SDavid du Colombier 		}
613e12c5d1SDavid du Colombier 	else {
623e12c5d1SDavid du Colombier 		L = -L >> Exp_shift;
633e12c5d1SDavid du Colombier 		if (L < Exp_shift) {
643e12c5d1SDavid du Colombier 			word0(a) = 0x80000 >> L;
653e12c5d1SDavid du Colombier 			word1(a) = 0;
663e12c5d1SDavid du Colombier 			}
673e12c5d1SDavid du Colombier 		else {
683e12c5d1SDavid du Colombier 			word0(a) = 0;
693e12c5d1SDavid du Colombier 			L -= Exp_shift;
703e12c5d1SDavid du Colombier 			word1(a) = L >= 31 ? 1 : 1 << 31 - L;
713e12c5d1SDavid du Colombier 			}
723e12c5d1SDavid du Colombier 		}
733e12c5d1SDavid du Colombier #endif
743e12c5d1SDavid du Colombier 	return a.d;
753e12c5d1SDavid du Colombier 	}
763e12c5d1SDavid du Colombier 
773e12c5d1SDavid du Colombier  static Bigint *
s2b(CONST char * s,int nd0,int nd,unsigned long y9)783e12c5d1SDavid du Colombier s2b(CONST char *s, int nd0, int nd, unsigned long y9)
793e12c5d1SDavid du Colombier {
803e12c5d1SDavid du Colombier 	Bigint *b;
813e12c5d1SDavid du Colombier 	int i, k;
823e12c5d1SDavid du Colombier 	long x, y;
833e12c5d1SDavid du Colombier 
843e12c5d1SDavid du Colombier 	x = (nd + 8) / 9;
853e12c5d1SDavid du Colombier 	for(k = 0, y = 1; x > y; y <<= 1, k++) ;
863e12c5d1SDavid du Colombier #ifdef Pack_32
873e12c5d1SDavid du Colombier 	b = Balloc(k);
883e12c5d1SDavid du Colombier 	b->x[0] = y9;
893e12c5d1SDavid du Colombier 	b->wds = 1;
903e12c5d1SDavid du Colombier #else
913e12c5d1SDavid du Colombier 	b = Balloc(k+1);
923e12c5d1SDavid du Colombier 	b->x[0] = y9 & 0xffff;
933e12c5d1SDavid du Colombier 	b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
943e12c5d1SDavid du Colombier #endif
953e12c5d1SDavid du Colombier 
963e12c5d1SDavid du Colombier 	i = 9;
973e12c5d1SDavid du Colombier 	if (9 < nd0) {
983e12c5d1SDavid du Colombier 		s += 9;
993e12c5d1SDavid du Colombier 		do b = multadd(b, 10, *s++ - '0');
1003e12c5d1SDavid du Colombier 			while(++i < nd0);
1013e12c5d1SDavid du Colombier 		s++;
1023e12c5d1SDavid du Colombier 		}
1033e12c5d1SDavid du Colombier 	else
1043e12c5d1SDavid du Colombier 		s += 10;
1053e12c5d1SDavid du Colombier 	for(; i < nd; i++)
1063e12c5d1SDavid du Colombier 		b = multadd(b, 10, *s++ - '0');
1073e12c5d1SDavid du Colombier 	return b;
1083e12c5d1SDavid du Colombier 	}
1093e12c5d1SDavid du Colombier 
1103e12c5d1SDavid du Colombier  static double
b2d(Bigint * a,int * e)1113e12c5d1SDavid du Colombier b2d(Bigint *a, int *e)
1123e12c5d1SDavid du Colombier {
1133e12c5d1SDavid du Colombier 	unsigned long *xa, *xa0, w, y, z;
1143e12c5d1SDavid du Colombier 	int k;
1153e12c5d1SDavid du Colombier 	Dul d;
1163e12c5d1SDavid du Colombier #ifdef VAX
1173e12c5d1SDavid du Colombier 	unsigned long d0, d1;
1183e12c5d1SDavid du Colombier #else
1193e12c5d1SDavid du Colombier #define d0 word0(d)
1203e12c5d1SDavid du Colombier #define d1 word1(d)
1213e12c5d1SDavid du Colombier #endif
1223e12c5d1SDavid du Colombier 
1233e12c5d1SDavid du Colombier 	xa0 = a->x;
1243e12c5d1SDavid du Colombier 	xa = xa0 + a->wds;
1253e12c5d1SDavid du Colombier 	y = *--xa;
1263e12c5d1SDavid du Colombier #ifdef DEBUG
1273e12c5d1SDavid du Colombier 	if (!y) Bug("zero y in b2d");
1283e12c5d1SDavid du Colombier #endif
1293e12c5d1SDavid du Colombier 	k = hi0bits(y);
1303e12c5d1SDavid du Colombier 	*e = 32 - k;
1313e12c5d1SDavid du Colombier #ifdef Pack_32
1323e12c5d1SDavid du Colombier 	if (k < Ebits) {
1333e12c5d1SDavid du Colombier 		d0 = Exp_1 | y >> Ebits - k;
1343e12c5d1SDavid du Colombier 		w = xa > xa0 ? *--xa : 0;
1353e12c5d1SDavid du Colombier 		d1 = y << (32-Ebits) + k | w >> Ebits - k;
1363e12c5d1SDavid du Colombier 		goto ret_d;
1373e12c5d1SDavid du Colombier 		}
1383e12c5d1SDavid du Colombier 	z = xa > xa0 ? *--xa : 0;
1393e12c5d1SDavid du Colombier 	if (k -= Ebits) {
1403e12c5d1SDavid du Colombier 		d0 = Exp_1 | y << k | z >> 32 - k;
1413e12c5d1SDavid du Colombier 		y = xa > xa0 ? *--xa : 0;
1423e12c5d1SDavid du Colombier 		d1 = z << k | y >> 32 - k;
1433e12c5d1SDavid du Colombier 		}
1443e12c5d1SDavid du Colombier 	else {
1453e12c5d1SDavid du Colombier 		d0 = Exp_1 | y;
1463e12c5d1SDavid du Colombier 		d1 = z;
1473e12c5d1SDavid du Colombier 		}
1483e12c5d1SDavid du Colombier #else
1493e12c5d1SDavid du Colombier 	if (k < Ebits + 16) {
1503e12c5d1SDavid du Colombier 		z = xa > xa0 ? *--xa : 0;
1513e12c5d1SDavid du Colombier 		d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1523e12c5d1SDavid du Colombier 		w = xa > xa0 ? *--xa : 0;
1533e12c5d1SDavid du Colombier 		y = xa > xa0 ? *--xa : 0;
1543e12c5d1SDavid du Colombier 		d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1553e12c5d1SDavid du Colombier 		goto ret_d;
1563e12c5d1SDavid du Colombier 		}
1573e12c5d1SDavid du Colombier 	z = xa > xa0 ? *--xa : 0;
1583e12c5d1SDavid du Colombier 	w = xa > xa0 ? *--xa : 0;
1593e12c5d1SDavid du Colombier 	k -= Ebits + 16;
1603e12c5d1SDavid du Colombier 	d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1613e12c5d1SDavid du Colombier 	y = xa > xa0 ? *--xa : 0;
1623e12c5d1SDavid du Colombier 	d1 = w << k + 16 | y << k;
1633e12c5d1SDavid du Colombier #endif
1643e12c5d1SDavid du Colombier  ret_d:
1653e12c5d1SDavid du Colombier #ifdef VAX
1663e12c5d1SDavid du Colombier 	word0(d) = d0 >> 16 | d0 << 16;
1673e12c5d1SDavid du Colombier 	word1(d) = d1 >> 16 | d1 << 16;
1683e12c5d1SDavid du Colombier #else
1693e12c5d1SDavid du Colombier #undef d0
1703e12c5d1SDavid du Colombier #undef d1
1713e12c5d1SDavid du Colombier #endif
1723e12c5d1SDavid du Colombier 	return d.d;
1733e12c5d1SDavid du Colombier 	}
1743e12c5d1SDavid du Colombier 
1753e12c5d1SDavid du Colombier  static double
ratio(Bigint * a,Bigint * b)1763e12c5d1SDavid du Colombier ratio(Bigint *a, Bigint *b)
1773e12c5d1SDavid du Colombier {
1783e12c5d1SDavid du Colombier 	Dul da, db;
1793e12c5d1SDavid du Colombier 	int k, ka, kb;
1803e12c5d1SDavid du Colombier 
1813e12c5d1SDavid du Colombier 	da.d = b2d(a, &ka);
1823e12c5d1SDavid du Colombier 	db.d = b2d(b, &kb);
1833e12c5d1SDavid du Colombier #ifdef Pack_32
1843e12c5d1SDavid du Colombier 	k = ka - kb + 32*(a->wds - b->wds);
1853e12c5d1SDavid du Colombier #else
1863e12c5d1SDavid du Colombier 	k = ka - kb + 16*(a->wds - b->wds);
1873e12c5d1SDavid du Colombier #endif
1883e12c5d1SDavid du Colombier #ifdef IBM
1893e12c5d1SDavid du Colombier 	if (k > 0) {
1903e12c5d1SDavid du Colombier 		word0(da) += (k >> 2)*Exp_msk1;
1913e12c5d1SDavid du Colombier 		if (k &= 3)
1923e12c5d1SDavid du Colombier 			da *= 1 << k;
1933e12c5d1SDavid du Colombier 		}
1943e12c5d1SDavid du Colombier 	else {
1953e12c5d1SDavid du Colombier 		k = -k;
1963e12c5d1SDavid du Colombier 		word0(db) += (k >> 2)*Exp_msk1;
1973e12c5d1SDavid du Colombier 		if (k &= 3)
1983e12c5d1SDavid du Colombier 			db *= 1 << k;
1993e12c5d1SDavid du Colombier 		}
2003e12c5d1SDavid du Colombier #else
2013e12c5d1SDavid du Colombier 	if (k > 0)
2023e12c5d1SDavid du Colombier 		word0(da) += k*Exp_msk1;
2033e12c5d1SDavid du Colombier 	else {
2043e12c5d1SDavid du Colombier 		k = -k;
2053e12c5d1SDavid du Colombier 		word0(db) += k*Exp_msk1;
2063e12c5d1SDavid du Colombier 		}
2073e12c5d1SDavid du Colombier #endif
2083e12c5d1SDavid du Colombier 	return da.d / db.d;
2093e12c5d1SDavid du Colombier 	}
2103e12c5d1SDavid du Colombier 
2113e12c5d1SDavid du Colombier  double
strtod(CONST char * s00,char ** se)2123e12c5d1SDavid du Colombier strtod(CONST char *s00, char **se)
2133e12c5d1SDavid du Colombier {
2143e12c5d1SDavid du Colombier 	int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
2153e12c5d1SDavid du Colombier 		 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
2163e12c5d1SDavid du Colombier 	CONST char *s, *s0, *s1;
2173e12c5d1SDavid du Colombier 	double aadj, aadj1, adj;
2183e12c5d1SDavid du Colombier 	Dul rv, rv0;
2193e12c5d1SDavid du Colombier 	long L;
2203e12c5d1SDavid du Colombier 	unsigned long y, z;
2213e12c5d1SDavid du Colombier 	Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
2223e12c5d1SDavid du Colombier 	sign = nz0 = nz = 0;
2233e12c5d1SDavid du Colombier 	rv.d = 0.;
2243e12c5d1SDavid du Colombier 	for(s = s00;;s++) switch(*s) {
2253e12c5d1SDavid du Colombier 		case '-':
2263e12c5d1SDavid du Colombier 			sign = 1;
2273e12c5d1SDavid du Colombier 			/* no break */
2283e12c5d1SDavid du Colombier 		case '+':
2293e12c5d1SDavid du Colombier 			if (*++s)
2303e12c5d1SDavid du Colombier 				goto break2;
2313e12c5d1SDavid du Colombier 			/* no break */
2323e12c5d1SDavid du Colombier 		case 0:
2333e12c5d1SDavid du Colombier 			s = s00;
2343e12c5d1SDavid du Colombier 			goto ret;
2353e12c5d1SDavid du Colombier 		case '\t':
2363e12c5d1SDavid du Colombier 		case '\n':
2373e12c5d1SDavid du Colombier 		case '\v':
2383e12c5d1SDavid du Colombier 		case '\f':
2393e12c5d1SDavid du Colombier 		case '\r':
2403e12c5d1SDavid du Colombier 		case ' ':
2413e12c5d1SDavid du Colombier 			continue;
2423e12c5d1SDavid du Colombier 		default:
2433e12c5d1SDavid du Colombier 			goto break2;
2443e12c5d1SDavid du Colombier 		}
2453e12c5d1SDavid du Colombier  break2:
2463e12c5d1SDavid du Colombier 	if (*s == '0') {
2473e12c5d1SDavid du Colombier 		nz0 = 1;
2483e12c5d1SDavid du Colombier 		while(*++s == '0') ;
2493e12c5d1SDavid du Colombier 		if (!*s)
2503e12c5d1SDavid du Colombier 			goto ret;
2513e12c5d1SDavid du Colombier 		}
2523e12c5d1SDavid du Colombier 	s0 = s;
2533e12c5d1SDavid du Colombier 	y = z = 0;
2543e12c5d1SDavid du Colombier 	for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
2553e12c5d1SDavid du Colombier 		if (nd < 9)
2563e12c5d1SDavid du Colombier 			y = 10*y + c - '0';
2573e12c5d1SDavid du Colombier 		else if (nd < 16)
2583e12c5d1SDavid du Colombier 			z = 10*z + c - '0';
2593e12c5d1SDavid du Colombier 	nd0 = nd;
2603e12c5d1SDavid du Colombier 	if (c == '.') {
2613e12c5d1SDavid du Colombier 		c = *++s;
2623e12c5d1SDavid du Colombier 		if (!nd) {
2633e12c5d1SDavid du Colombier 			for(; c == '0'; c = *++s)
2643e12c5d1SDavid du Colombier 				nz++;
2653e12c5d1SDavid du Colombier 			if (c > '0' && c <= '9') {
2663e12c5d1SDavid du Colombier 				s0 = s;
2673e12c5d1SDavid du Colombier 				nf += nz;
2683e12c5d1SDavid du Colombier 				nz = 0;
2693e12c5d1SDavid du Colombier 				goto have_dig;
2703e12c5d1SDavid du Colombier 				}
2713e12c5d1SDavid du Colombier 			goto dig_done;
2723e12c5d1SDavid du Colombier 			}
2733e12c5d1SDavid du Colombier 		for(; c >= '0' && c <= '9'; c = *++s) {
2743e12c5d1SDavid du Colombier  have_dig:
2753e12c5d1SDavid du Colombier 			nz++;
2763e12c5d1SDavid du Colombier 			if (c -= '0') {
2773e12c5d1SDavid du Colombier 				nf += nz;
2783e12c5d1SDavid du Colombier 				for(i = 1; i < nz; i++)
2793e12c5d1SDavid du Colombier 					if (nd++ < 9)
2803e12c5d1SDavid du Colombier 						y *= 10;
2813e12c5d1SDavid du Colombier 					else if (nd <= DBL_DIG + 1)
2823e12c5d1SDavid du Colombier 						z *= 10;
2833e12c5d1SDavid du Colombier 				if (nd++ < 9)
2843e12c5d1SDavid du Colombier 					y = 10*y + c;
2853e12c5d1SDavid du Colombier 				else if (nd <= DBL_DIG + 1)
2863e12c5d1SDavid du Colombier 					z = 10*z + c;
2873e12c5d1SDavid du Colombier 				nz = 0;
2883e12c5d1SDavid du Colombier 				}
2893e12c5d1SDavid du Colombier 			}
2903e12c5d1SDavid du Colombier 		}
2913e12c5d1SDavid du Colombier  dig_done:
2923e12c5d1SDavid du Colombier 	e = 0;
2933e12c5d1SDavid du Colombier 	if (c == 'e' || c == 'E') {
2943e12c5d1SDavid du Colombier 		if (!nd && !nz && !nz0) {
2953e12c5d1SDavid du Colombier 			s = s00;
2963e12c5d1SDavid du Colombier 			goto ret;
2973e12c5d1SDavid du Colombier 			}
2983e12c5d1SDavid du Colombier 		s00 = s;
2993e12c5d1SDavid du Colombier 		esign = 0;
3003e12c5d1SDavid du Colombier 		switch(c = *++s) {
3013e12c5d1SDavid du Colombier 			case '-':
3023e12c5d1SDavid du Colombier 				esign = 1;
3033e12c5d1SDavid du Colombier 			case '+':
3043e12c5d1SDavid du Colombier 				c = *++s;
3053e12c5d1SDavid du Colombier 			}
3063e12c5d1SDavid du Colombier 		if (c >= '0' && c <= '9') {
3073e12c5d1SDavid du Colombier 			while(c == '0')
3083e12c5d1SDavid du Colombier 				c = *++s;
3093e12c5d1SDavid du Colombier 			if (c > '0' && c <= '9') {
3103e12c5d1SDavid du Colombier 				e = c - '0';
3113e12c5d1SDavid du Colombier 				s1 = s;
3123e12c5d1SDavid du Colombier 				while((c = *++s) >= '0' && c <= '9')
3133e12c5d1SDavid du Colombier 					e = 10*e + c - '0';
3143e12c5d1SDavid du Colombier 				if (s - s1 > 8)
3153e12c5d1SDavid du Colombier 					/* Avoid confusion from exponents
3163e12c5d1SDavid du Colombier 					 * so large that e might overflow.
3173e12c5d1SDavid du Colombier 					 */
3183e12c5d1SDavid du Colombier 					e = 9999999;
3193e12c5d1SDavid du Colombier 				if (esign)
3203e12c5d1SDavid du Colombier 					e = -e;
3213e12c5d1SDavid du Colombier 				}
3223e12c5d1SDavid du Colombier 			else
3233e12c5d1SDavid du Colombier 				e = 0;
3243e12c5d1SDavid du Colombier 			}
3253e12c5d1SDavid du Colombier 		else
3263e12c5d1SDavid du Colombier 			s = s00;
3273e12c5d1SDavid du Colombier 		}
3283e12c5d1SDavid du Colombier 	if (!nd) {
3293e12c5d1SDavid du Colombier 		if (!nz && !nz0)
3303e12c5d1SDavid du Colombier 			s = s00;
3313e12c5d1SDavid du Colombier 		goto ret;
3323e12c5d1SDavid du Colombier 		}
3333e12c5d1SDavid du Colombier 	e1 = e -= nf;
3343e12c5d1SDavid du Colombier 
3353e12c5d1SDavid du Colombier 	/* Now we have nd0 digits, starting at s0, followed by a
3363e12c5d1SDavid du Colombier 	 * decimal point, followed by nd-nd0 digits.  The number we're
3373e12c5d1SDavid du Colombier 	 * after is the integer represented by those digits times
3383e12c5d1SDavid du Colombier 	 * 10**e */
3393e12c5d1SDavid du Colombier 
3403e12c5d1SDavid du Colombier 	if (!nd0)
3413e12c5d1SDavid du Colombier 		nd0 = nd;
3423e12c5d1SDavid du Colombier 	k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
3433e12c5d1SDavid du Colombier 	rv.d = y;
3443e12c5d1SDavid du Colombier 	if (k > 9)
3453e12c5d1SDavid du Colombier 		rv.d = tens[k - 9] * rv.d + z;
346219b2ee8SDavid du Colombier 	bd0 = 0;
3473e12c5d1SDavid du Colombier 	if (nd <= DBL_DIG
3483e12c5d1SDavid du Colombier #ifndef RND_PRODQUOT
3493e12c5d1SDavid du Colombier 		&& FLT_ROUNDS == 1
3503e12c5d1SDavid du Colombier #endif
3513e12c5d1SDavid du Colombier 			) {
3523e12c5d1SDavid du Colombier 		if (!e)
3533e12c5d1SDavid du Colombier 			goto ret;
3543e12c5d1SDavid du Colombier 		if (e > 0) {
3553e12c5d1SDavid du Colombier 			if (e <= Ten_pmax) {
3563e12c5d1SDavid du Colombier #ifdef VAX
3573e12c5d1SDavid du Colombier 				goto vax_ovfl_check;
3583e12c5d1SDavid du Colombier #else
3593e12c5d1SDavid du Colombier 				/* rv = */ rounded_product(rv.d, tens[e]);
3603e12c5d1SDavid du Colombier 				goto ret;
3613e12c5d1SDavid du Colombier #endif
3623e12c5d1SDavid du Colombier 				}
3633e12c5d1SDavid du Colombier 			i = DBL_DIG - nd;
3643e12c5d1SDavid du Colombier 			if (e <= Ten_pmax + i) {
3653e12c5d1SDavid du Colombier 				/* A fancier test would sometimes let us do
3663e12c5d1SDavid du Colombier 				 * this for larger i values.
3673e12c5d1SDavid du Colombier 				 */
3683e12c5d1SDavid du Colombier 				e -= i;
3693e12c5d1SDavid du Colombier 				rv.d *= tens[i];
3703e12c5d1SDavid du Colombier #ifdef VAX
3713e12c5d1SDavid du Colombier 				/* VAX exponent range is so narrow we must
3723e12c5d1SDavid du Colombier 				 * worry about overflow here...
3733e12c5d1SDavid du Colombier 				 */
3743e12c5d1SDavid du Colombier  vax_ovfl_check:
3753e12c5d1SDavid du Colombier 				word0(rv) -= P*Exp_msk1;
3763e12c5d1SDavid du Colombier 				/* rv = */ rounded_product(rv.d, tens[e]);
3773e12c5d1SDavid du Colombier 				if ((word0(rv) & Exp_mask)
3783e12c5d1SDavid du Colombier 				 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
3793e12c5d1SDavid du Colombier 					goto ovfl;
3803e12c5d1SDavid du Colombier 				word0(rv) += P*Exp_msk1;
3813e12c5d1SDavid du Colombier #else
3823e12c5d1SDavid du Colombier 				/* rv = */ rounded_product(rv.d, tens[e]);
3833e12c5d1SDavid du Colombier #endif
3843e12c5d1SDavid du Colombier 				goto ret;
3853e12c5d1SDavid du Colombier 				}
3863e12c5d1SDavid du Colombier 			}
3873e12c5d1SDavid du Colombier 		else if (e >= -Ten_pmax) {
3883e12c5d1SDavid du Colombier 			/* rv = */ rounded_quotient(rv.d, tens[-e]);
3893e12c5d1SDavid du Colombier 			goto ret;
3903e12c5d1SDavid du Colombier 			}
3913e12c5d1SDavid du Colombier 		}
3923e12c5d1SDavid du Colombier 	e1 += nd - k;
3933e12c5d1SDavid du Colombier 
3943e12c5d1SDavid du Colombier 	/* Get starting approximation = rv * 10**e1 */
3953e12c5d1SDavid du Colombier 
3963e12c5d1SDavid du Colombier 	if (e1 > 0) {
397*43aadf5eSDavid du Colombier 		if (nd0 + e1 - 1 > DBL_MAX_10_EXP)
398*43aadf5eSDavid du Colombier 			goto ovfl;
3993e12c5d1SDavid du Colombier 		if (i = e1 & 15)
4003e12c5d1SDavid du Colombier 			rv.d *= tens[i];
4013e12c5d1SDavid du Colombier 		if (e1 &= ~15) {
4023e12c5d1SDavid du Colombier 			if (e1 > DBL_MAX_10_EXP) {
4033e12c5d1SDavid du Colombier  ovfl:
4043e12c5d1SDavid du Colombier 				errno = ERANGE;
4053e12c5d1SDavid du Colombier 				rv.d = HUGE_VAL;
406219b2ee8SDavid du Colombier 				if (bd0)
407219b2ee8SDavid du Colombier 					goto retfree;
4083e12c5d1SDavid du Colombier 				goto ret;
4093e12c5d1SDavid du Colombier 				}
4103e12c5d1SDavid du Colombier 			if (e1 >>= 4) {
4113e12c5d1SDavid du Colombier 				for(j = 0; e1 > 1; j++, e1 >>= 1)
4123e12c5d1SDavid du Colombier 					if (e1 & 1)
4133e12c5d1SDavid du Colombier 						rv.d *= bigtens[j];
4143e12c5d1SDavid du Colombier 			/* The last multiplication could overflow. */
4153e12c5d1SDavid du Colombier 				word0(rv) -= P*Exp_msk1;
4163e12c5d1SDavid du Colombier 				rv.d *= bigtens[j];
4173e12c5d1SDavid du Colombier 				if ((z = word0(rv) & Exp_mask)
4183e12c5d1SDavid du Colombier 				 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
4193e12c5d1SDavid du Colombier 					goto ovfl;
4203e12c5d1SDavid du Colombier 				if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
4213e12c5d1SDavid du Colombier 					/* set to largest number */
4223e12c5d1SDavid du Colombier 					/* (Can't trust DBL_MAX) */
4233e12c5d1SDavid du Colombier 					word0(rv) = Big0;
4243e12c5d1SDavid du Colombier 					word1(rv) = Big1;
4253e12c5d1SDavid du Colombier 					}
4263e12c5d1SDavid du Colombier 				else
4273e12c5d1SDavid du Colombier 					word0(rv) += P*Exp_msk1;
4283e12c5d1SDavid du Colombier 				}
4293e12c5d1SDavid du Colombier 
4303e12c5d1SDavid du Colombier 			}
4313e12c5d1SDavid du Colombier 		}
4323e12c5d1SDavid du Colombier 	else if (e1 < 0) {
4333e12c5d1SDavid du Colombier 		e1 = -e1;
4343e12c5d1SDavid du Colombier 		if (i = e1 & 15)
4353e12c5d1SDavid du Colombier 			rv.d /= tens[i];
4363e12c5d1SDavid du Colombier 		if (e1 &= ~15) {
4373e12c5d1SDavid du Colombier 			e1 >>= 4;
438219b2ee8SDavid du Colombier 			if (e1 >= 1 << n_bigtens)
439219b2ee8SDavid du Colombier 				goto undfl;
4403e12c5d1SDavid du Colombier 			for(j = 0; e1 > 1; j++, e1 >>= 1)
4413e12c5d1SDavid du Colombier 				if (e1 & 1)
4423e12c5d1SDavid du Colombier 					rv.d *= tinytens[j];
4433e12c5d1SDavid du Colombier 			/* The last multiplication could underflow. */
4443e12c5d1SDavid du Colombier 			rv0.d = rv.d;
4453e12c5d1SDavid du Colombier 			rv.d *= tinytens[j];
4463e12c5d1SDavid du Colombier 			if (rv.d == 0) {
4473e12c5d1SDavid du Colombier 				rv.d = 2.*rv0.d;
4483e12c5d1SDavid du Colombier 				rv.d *= tinytens[j];
4493e12c5d1SDavid du Colombier 				if (rv.d == 0) {
4503e12c5d1SDavid du Colombier  undfl:
4513e12c5d1SDavid du Colombier 					rv.d = 0.;
4523e12c5d1SDavid du Colombier 					errno = ERANGE;
453219b2ee8SDavid du Colombier 					if (bd0)
454219b2ee8SDavid du Colombier 						goto retfree;
4553e12c5d1SDavid du Colombier 					goto ret;
4563e12c5d1SDavid du Colombier 					}
4573e12c5d1SDavid du Colombier 				word0(rv) = Tiny0;
4583e12c5d1SDavid du Colombier 				word1(rv) = Tiny1;
4593e12c5d1SDavid du Colombier 				/* The refinement below will clean
4603e12c5d1SDavid du Colombier 				 * this approximation up.
4613e12c5d1SDavid du Colombier 				 */
4623e12c5d1SDavid du Colombier 				}
4633e12c5d1SDavid du Colombier 			}
4643e12c5d1SDavid du Colombier 		}
4653e12c5d1SDavid du Colombier 
4663e12c5d1SDavid du Colombier 	/* Now the hard part -- adjusting rv to the correct value.*/
4673e12c5d1SDavid du Colombier 
4683e12c5d1SDavid du Colombier 	/* Put digits into bd: true value = bd * 10^e */
4693e12c5d1SDavid du Colombier 
4703e12c5d1SDavid du Colombier 	bd0 = s2b(s0, nd0, nd, y);
4713e12c5d1SDavid du Colombier 
4723e12c5d1SDavid du Colombier 	for(;;) {
4733e12c5d1SDavid du Colombier 		bd = Balloc(bd0->k);
4743e12c5d1SDavid du Colombier 		Bcopy(bd, bd0);
4753e12c5d1SDavid du Colombier 		bb = d2b(rv.d, &bbe, &bbbits);	/* rv = bb * 2^bbe */
4763e12c5d1SDavid du Colombier 		bs = i2b(1);
4773e12c5d1SDavid du Colombier 
4783e12c5d1SDavid du Colombier 		if (e >= 0) {
4793e12c5d1SDavid du Colombier 			bb2 = bb5 = 0;
4803e12c5d1SDavid du Colombier 			bd2 = bd5 = e;
4813e12c5d1SDavid du Colombier 			}
4823e12c5d1SDavid du Colombier 		else {
4833e12c5d1SDavid du Colombier 			bb2 = bb5 = -e;
4843e12c5d1SDavid du Colombier 			bd2 = bd5 = 0;
4853e12c5d1SDavid du Colombier 			}
4863e12c5d1SDavid du Colombier 		if (bbe >= 0)
4873e12c5d1SDavid du Colombier 			bb2 += bbe;
4883e12c5d1SDavid du Colombier 		else
4893e12c5d1SDavid du Colombier 			bd2 -= bbe;
4903e12c5d1SDavid du Colombier 		bs2 = bb2;
4913e12c5d1SDavid du Colombier #ifdef Sudden_Underflow
4923e12c5d1SDavid du Colombier #ifdef IBM
4933e12c5d1SDavid du Colombier 		j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
4943e12c5d1SDavid du Colombier #else
4953e12c5d1SDavid du Colombier 		j = P + 1 - bbbits;
4963e12c5d1SDavid du Colombier #endif
4973e12c5d1SDavid du Colombier #else
4983e12c5d1SDavid du Colombier 		i = bbe + bbbits - 1;	/* logb(rv) */
4993e12c5d1SDavid du Colombier 		if (i < Emin)	/* denormal */
5003e12c5d1SDavid du Colombier 			j = bbe + (P-Emin);
5013e12c5d1SDavid du Colombier 		else
5023e12c5d1SDavid du Colombier 			j = P + 1 - bbbits;
5033e12c5d1SDavid du Colombier #endif
5043e12c5d1SDavid du Colombier 		bb2 += j;
5053e12c5d1SDavid du Colombier 		bd2 += j;
5063e12c5d1SDavid du Colombier 		i = bb2 < bd2 ? bb2 : bd2;
5073e12c5d1SDavid du Colombier 		if (i > bs2)
5083e12c5d1SDavid du Colombier 			i = bs2;
5093e12c5d1SDavid du Colombier 		if (i > 0) {
5103e12c5d1SDavid du Colombier 			bb2 -= i;
5113e12c5d1SDavid du Colombier 			bd2 -= i;
5123e12c5d1SDavid du Colombier 			bs2 -= i;
5133e12c5d1SDavid du Colombier 			}
5143e12c5d1SDavid du Colombier 		if (bb5 > 0) {
5153e12c5d1SDavid du Colombier 			bs = pow5mult(bs, bb5);
5163e12c5d1SDavid du Colombier 			bb1 = mult(bs, bb);
5173e12c5d1SDavid du Colombier 			Bfree(bb);
5183e12c5d1SDavid du Colombier 			bb = bb1;
5193e12c5d1SDavid du Colombier 			}
5203e12c5d1SDavid du Colombier 		if (bb2 > 0)
5213e12c5d1SDavid du Colombier 			bb = lshift(bb, bb2);
5223e12c5d1SDavid du Colombier 		if (bd5 > 0)
5233e12c5d1SDavid du Colombier 			bd = pow5mult(bd, bd5);
5243e12c5d1SDavid du Colombier 		if (bd2 > 0)
5253e12c5d1SDavid du Colombier 			bd = lshift(bd, bd2);
5263e12c5d1SDavid du Colombier 		if (bs2 > 0)
5273e12c5d1SDavid du Colombier 			bs = lshift(bs, bs2);
5283e12c5d1SDavid du Colombier 		delta = diff(bb, bd);
5293e12c5d1SDavid du Colombier 		dsign = delta->sign;
5303e12c5d1SDavid du Colombier 		delta->sign = 0;
5313e12c5d1SDavid du Colombier 		i = cmp(delta, bs);
5323e12c5d1SDavid du Colombier 		if (i < 0) {
5333e12c5d1SDavid du Colombier 			/* Error is less than half an ulp -- check for
5343e12c5d1SDavid du Colombier 			 * special case of mantissa a power of two.
5353e12c5d1SDavid du Colombier 			 */
5363e12c5d1SDavid du Colombier 			if (dsign || word1(rv) || word0(rv) & Bndry_mask)
5373e12c5d1SDavid du Colombier 				break;
5383e12c5d1SDavid du Colombier 			delta = lshift(delta,Log2P);
5393e12c5d1SDavid du Colombier 			if (cmp(delta, bs) > 0)
5403e12c5d1SDavid du Colombier 				goto drop_down;
5413e12c5d1SDavid du Colombier 			break;
5423e12c5d1SDavid du Colombier 			}
5433e12c5d1SDavid du Colombier 		if (i == 0) {
5443e12c5d1SDavid du Colombier 			/* exactly half-way between */
5453e12c5d1SDavid du Colombier 			if (dsign) {
5463e12c5d1SDavid du Colombier 				if ((word0(rv) & Bndry_mask1) == Bndry_mask1
5473e12c5d1SDavid du Colombier 				 &&  word1(rv) == 0xffffffff) {
5483e12c5d1SDavid du Colombier 					/*boundary case -- increment exponent*/
5493e12c5d1SDavid du Colombier 					word0(rv) = (word0(rv) & Exp_mask)
5503e12c5d1SDavid du Colombier 						+ Exp_msk1
5513e12c5d1SDavid du Colombier #ifdef IBM
5523e12c5d1SDavid du Colombier 						| Exp_msk1 >> 4
5533e12c5d1SDavid du Colombier #endif
5543e12c5d1SDavid du Colombier 						;
5553e12c5d1SDavid du Colombier 					word1(rv) = 0;
5563e12c5d1SDavid du Colombier 					break;
5573e12c5d1SDavid du Colombier 					}
5583e12c5d1SDavid du Colombier 				}
5593e12c5d1SDavid du Colombier 			else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
5603e12c5d1SDavid du Colombier  drop_down:
5613e12c5d1SDavid du Colombier 				/* boundary case -- decrement exponent */
5623e12c5d1SDavid du Colombier #ifdef Sudden_Underflow
5633e12c5d1SDavid du Colombier 				L = word0(rv) & Exp_mask;
5643e12c5d1SDavid du Colombier #ifdef IBM
5653e12c5d1SDavid du Colombier 				if (L <  Exp_msk1)
5663e12c5d1SDavid du Colombier #else
5673e12c5d1SDavid du Colombier 				if (L <= Exp_msk1)
5683e12c5d1SDavid du Colombier #endif
5693e12c5d1SDavid du Colombier 					goto undfl;
5703e12c5d1SDavid du Colombier 				L -= Exp_msk1;
5713e12c5d1SDavid du Colombier #else
5723e12c5d1SDavid du Colombier 				L = (word0(rv) & Exp_mask) - Exp_msk1;
5733e12c5d1SDavid du Colombier #endif
5743e12c5d1SDavid du Colombier 				word0(rv) = L | Bndry_mask1;
5753e12c5d1SDavid du Colombier 				word1(rv) = 0xffffffff;
5763e12c5d1SDavid du Colombier #ifdef IBM
5773e12c5d1SDavid du Colombier 				goto cont;
5783e12c5d1SDavid du Colombier #else
5793e12c5d1SDavid du Colombier 				break;
5803e12c5d1SDavid du Colombier #endif
5813e12c5d1SDavid du Colombier 				}
5823e12c5d1SDavid du Colombier #ifndef ROUND_BIASED
5833e12c5d1SDavid du Colombier 			if (!(word1(rv) & LSB))
5843e12c5d1SDavid du Colombier 				break;
5853e12c5d1SDavid du Colombier #endif
5863e12c5d1SDavid du Colombier 			if (dsign)
5873e12c5d1SDavid du Colombier 				rv.d += ulp(rv.d);
5883e12c5d1SDavid du Colombier #ifndef ROUND_BIASED
5893e12c5d1SDavid du Colombier 			else {
5903e12c5d1SDavid du Colombier 				rv.d -= ulp(rv.d);
5913e12c5d1SDavid du Colombier #ifndef Sudden_Underflow
5923e12c5d1SDavid du Colombier 				if (rv.d == 0)
5933e12c5d1SDavid du Colombier 					goto undfl;
5943e12c5d1SDavid du Colombier #endif
5953e12c5d1SDavid du Colombier 				}
5963e12c5d1SDavid du Colombier #endif
5973e12c5d1SDavid du Colombier 			break;
5983e12c5d1SDavid du Colombier 			}
5993e12c5d1SDavid du Colombier 		if ((aadj = ratio(delta, bs)) <= 2.) {
6003e12c5d1SDavid du Colombier 			if (dsign)
6013e12c5d1SDavid du Colombier 				aadj = aadj1 = 1.;
6023e12c5d1SDavid du Colombier 			else if (word1(rv) || word0(rv) & Bndry_mask) {
6033e12c5d1SDavid du Colombier #ifndef Sudden_Underflow
6043e12c5d1SDavid du Colombier 				if (word1(rv) == Tiny1 && !word0(rv))
6053e12c5d1SDavid du Colombier 					goto undfl;
6063e12c5d1SDavid du Colombier #endif
6073e12c5d1SDavid du Colombier 				aadj = 1.;
6083e12c5d1SDavid du Colombier 				aadj1 = -1.;
6093e12c5d1SDavid du Colombier 				}
6103e12c5d1SDavid du Colombier 			else {
6113e12c5d1SDavid du Colombier 				/* special case -- power of FLT_RADIX to be */
6123e12c5d1SDavid du Colombier 				/* rounded down... */
6133e12c5d1SDavid du Colombier 
6143e12c5d1SDavid du Colombier 				if (aadj < 2./FLT_RADIX)
6153e12c5d1SDavid du Colombier 					aadj = 1./FLT_RADIX;
6163e12c5d1SDavid du Colombier 				else
6173e12c5d1SDavid du Colombier 					aadj *= 0.5;
6183e12c5d1SDavid du Colombier 				aadj1 = -aadj;
6193e12c5d1SDavid du Colombier 				}
6203e12c5d1SDavid du Colombier 			}
6213e12c5d1SDavid du Colombier 		else {
6223e12c5d1SDavid du Colombier 			aadj *= 0.5;
6233e12c5d1SDavid du Colombier 			aadj1 = dsign ? aadj : -aadj;
6243e12c5d1SDavid du Colombier #ifdef Check_FLT_ROUNDS
6253e12c5d1SDavid du Colombier 			switch(FLT_ROUNDS) {
6263e12c5d1SDavid du Colombier 				case 2: /* towards +infinity */
6273e12c5d1SDavid du Colombier 					aadj1 -= 0.5;
6283e12c5d1SDavid du Colombier 					break;
6293e12c5d1SDavid du Colombier 				case 0: /* towards 0 */
6303e12c5d1SDavid du Colombier 				case 3: /* towards -infinity */
6313e12c5d1SDavid du Colombier 					aadj1 += 0.5;
6323e12c5d1SDavid du Colombier 				}
6333e12c5d1SDavid du Colombier #else
6343e12c5d1SDavid du Colombier 			if (FLT_ROUNDS == 0)
6353e12c5d1SDavid du Colombier 				aadj1 += 0.5;
6363e12c5d1SDavid du Colombier #endif
6373e12c5d1SDavid du Colombier 			}
6383e12c5d1SDavid du Colombier 		y = word0(rv) & Exp_mask;
6393e12c5d1SDavid du Colombier 
6403e12c5d1SDavid du Colombier 		/* Check for overflow */
6413e12c5d1SDavid du Colombier 
6423e12c5d1SDavid du Colombier 		if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
6433e12c5d1SDavid du Colombier 			rv0.d = rv.d;
6443e12c5d1SDavid du Colombier 			word0(rv) -= P*Exp_msk1;
6453e12c5d1SDavid du Colombier 			adj = aadj1 * ulp(rv.d);
6463e12c5d1SDavid du Colombier 			rv.d += adj;
6473e12c5d1SDavid du Colombier 			if ((word0(rv) & Exp_mask) >=
6483e12c5d1SDavid du Colombier 					Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
6493e12c5d1SDavid du Colombier 				if (word0(rv0) == Big0 && word1(rv0) == Big1)
6503e12c5d1SDavid du Colombier 					goto ovfl;
6513e12c5d1SDavid du Colombier 				word0(rv) = Big0;
6523e12c5d1SDavid du Colombier 				word1(rv) = Big1;
6533e12c5d1SDavid du Colombier 				goto cont;
6543e12c5d1SDavid du Colombier 				}
6553e12c5d1SDavid du Colombier 			else
6563e12c5d1SDavid du Colombier 				word0(rv) += P*Exp_msk1;
6573e12c5d1SDavid du Colombier 			}
6583e12c5d1SDavid du Colombier 		else {
6593e12c5d1SDavid du Colombier #ifdef Sudden_Underflow
6603e12c5d1SDavid du Colombier 			if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
6613e12c5d1SDavid du Colombier 				rv0.d = rv.d;
6623e12c5d1SDavid du Colombier 				word0(rv) += P*Exp_msk1;
6633e12c5d1SDavid du Colombier 				adj = aadj1 * ulp(rv.d);
6643e12c5d1SDavid du Colombier 				rv.d += adj;
6653e12c5d1SDavid du Colombier #ifdef IBM
6663e12c5d1SDavid du Colombier 				if ((word0(rv) & Exp_mask) <  P*Exp_msk1)
6673e12c5d1SDavid du Colombier #else
6683e12c5d1SDavid du Colombier 				if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
6693e12c5d1SDavid du Colombier #endif
6703e12c5d1SDavid du Colombier 					{
6713e12c5d1SDavid du Colombier 					if (word0(rv0) == Tiny0
6723e12c5d1SDavid du Colombier 					 && word1(rv0) == Tiny1)
6733e12c5d1SDavid du Colombier 						goto undfl;
6743e12c5d1SDavid du Colombier 					word0(rv) = Tiny0;
6753e12c5d1SDavid du Colombier 					word1(rv) = Tiny1;
6763e12c5d1SDavid du Colombier 					goto cont;
6773e12c5d1SDavid du Colombier 					}
6783e12c5d1SDavid du Colombier 				else
6793e12c5d1SDavid du Colombier 					word0(rv) -= P*Exp_msk1;
6803e12c5d1SDavid du Colombier 				}
6813e12c5d1SDavid du Colombier 			else {
6823e12c5d1SDavid du Colombier 				adj = aadj1 * ulp(rv.d);
6833e12c5d1SDavid du Colombier 				rv.d += adj;
6843e12c5d1SDavid du Colombier 				}
6853e12c5d1SDavid du Colombier #else
6863e12c5d1SDavid du Colombier 			/* Compute adj so that the IEEE rounding rules will
6873e12c5d1SDavid du Colombier 			 * correctly round rv + adj in some half-way cases.
6883e12c5d1SDavid du Colombier 			 * If rv * ulp(rv) is denormalized (i.e.,
6893e12c5d1SDavid du Colombier 			 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
6903e12c5d1SDavid du Colombier 			 * trouble from bits lost to denormalization;
6913e12c5d1SDavid du Colombier 			 * example: 1.2e-307 .
6923e12c5d1SDavid du Colombier 			 */
6933e12c5d1SDavid du Colombier 			if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
6943e12c5d1SDavid du Colombier 				aadj1 = (double)(int)(aadj + 0.5);
6953e12c5d1SDavid du Colombier 				if (!dsign)
6963e12c5d1SDavid du Colombier 					aadj1 = -aadj1;
6973e12c5d1SDavid du Colombier 				}
6983e12c5d1SDavid du Colombier 			adj = aadj1 * ulp(rv.d);
6993e12c5d1SDavid du Colombier 			rv.d += adj;
7003e12c5d1SDavid du Colombier #endif
7013e12c5d1SDavid du Colombier 			}
7023e12c5d1SDavid du Colombier 		z = word0(rv) & Exp_mask;
7033e12c5d1SDavid du Colombier 		if (y == z) {
7043e12c5d1SDavid du Colombier 			/* Can we stop now? */
7053e12c5d1SDavid du Colombier 			L = aadj;
7063e12c5d1SDavid du Colombier 			aadj -= L;
7073e12c5d1SDavid du Colombier 			/* The tolerances below are conservative. */
7083e12c5d1SDavid du Colombier 			if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
7093e12c5d1SDavid du Colombier 				if (aadj < .4999999 || aadj > .5000001)
7103e12c5d1SDavid du Colombier 					break;
7113e12c5d1SDavid du Colombier 				}
7123e12c5d1SDavid du Colombier 			else if (aadj < .4999999/FLT_RADIX)
7133e12c5d1SDavid du Colombier 				break;
7143e12c5d1SDavid du Colombier 			}
7153e12c5d1SDavid du Colombier  cont:
7163e12c5d1SDavid du Colombier 		Bfree(bb);
7173e12c5d1SDavid du Colombier 		Bfree(bd);
7183e12c5d1SDavid du Colombier 		Bfree(bs);
7193e12c5d1SDavid du Colombier 		Bfree(delta);
7203e12c5d1SDavid du Colombier 		}
721219b2ee8SDavid du Colombier  retfree:
7223e12c5d1SDavid du Colombier 	Bfree(bb);
7233e12c5d1SDavid du Colombier 	Bfree(bd);
7243e12c5d1SDavid du Colombier 	Bfree(bs);
7253e12c5d1SDavid du Colombier 	Bfree(bd0);
7263e12c5d1SDavid du Colombier 	Bfree(delta);
7273e12c5d1SDavid du Colombier  ret:
7283e12c5d1SDavid du Colombier 	if (se)
7293e12c5d1SDavid du Colombier 		*se = (char *)s;
7303e12c5d1SDavid du Colombier 	return sign ? -rv.d : rv.d;
7313e12c5d1SDavid du Colombier 	}
732