1 /****************************************************************
2 
3 The author of this software is David M. Gay.
4 
5 Copyright (C) 1998-2001 by Lucent Technologies
6 All Rights Reserved
7 
8 Permission to use, copy, modify, and distribute this software and
9 its documentation for any purpose and without fee is hereby
10 granted, provided that the above copyright notice appear in all
11 copies and that both that the copyright notice and this
12 permission notice and warranty disclaimer appear in supporting
13 documentation, and that the name of Lucent or any of its entities
14 not be used in advertising or publicity pertaining to
15 distribution of the software without specific, written prior
16 permission.
17 
18 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
19 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
20 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
21 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
22 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
23 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
24 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
25 THIS 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 <_ansi.h>
33 #include <errno.h>
34 #include <stdlib.h>
35 #include <string.h>
36 #include "mprec.h"
37 #include "gdtoa.h"
38 #include "gd_qnan.h"
39 
40 #include "locale.h"
41 
42 #if defined (_HAVE_LONG_DOUBLE) && !defined (_LDBL_EQ_DBL)
43 
44 #define USE_LOCALE
45 
46  static const int
47 fivesbits[] = {	 0,  3,  5,  7, 10, 12, 14, 17, 19, 21,
48 		24, 26, 28, 31, 33, 35, 38, 40, 42, 45,
49 		47, 49, 52
50 #ifdef VAX
51 		, 54, 56
52 #endif
53 		};
54 
55 static _Bigint *
56 #ifdef KR_headers
sum(p,a,b)57 sum(p, a, b) struct _reent *p; _Bigint *a; _Bigint *b;
58 #else
59 sum(struct _reent *p, _Bigint *a, _Bigint *b)
60 #endif
61 {
62 	_Bigint *c;
63 	__ULong carry, *xc, *xa, *xb, *xe, y;
64 #ifdef Pack_32
65 	__ULong z;
66 #endif
67 
68 	if (a->_wds < b->_wds) {
69 		c = b; b = a; a = c;
70 		}
71 	c = Balloc(p, a->_k);
72 	c->_wds = a->_wds;
73 	carry = 0;
74 	xa = a->_x;
75 	xb = b->_x;
76 	xc = c->_x;
77 	xe = xc + b->_wds;
78 #ifdef Pack_32
79 	do {
80 		y = (*xa & 0xffff) + (*xb & 0xffff) + carry;
81 		carry = (y & 0x10000) >> 16;
82 		z = (*xa++ >> 16) + (*xb++ >> 16) + carry;
83 		carry = (z & 0x10000) >> 16;
84 		Storeinc(xc, z, y);
85 		}
86 		while(xc < xe);
87 	xe += a->_wds - b->_wds;
88 	while(xc < xe) {
89 		y = (*xa & 0xffff) + carry;
90 		carry = (y & 0x10000) >> 16;
91 		z = (*xa++ >> 16) + carry;
92 		carry = (z & 0x10000) >> 16;
93 		Storeinc(xc, z, y);
94 		}
95 #else
96 	do {
97 		y = *xa++ + *xb++ + carry;
98 		carry = (y & 0x10000) >> 16;
99 		*xc++ = y & 0xffff;
100 		}
101 		while(xc < xe);
102 	xe += a->_wds - b->_wds;
103 	while(xc < xe) {
104 		y = *xa++ + carry;
105 		carry = (y & 0x10000) >> 16;
106 		*xc++ = y & 0xffff;
107 		}
108 #endif
109 	if (carry) {
110 		if (c->_wds == c->_maxwds) {
111 			b = Balloc(p, c->_k + 1);
112 			Bcopy(b, c);
113 			Bfree(p, c);
114 			c = b;
115 			}
116 		c->_x[c->_wds++] = 1;
117 		}
118 	return c;
119 	}
120 
121 static void
122 #ifdef KR_headers
rshift(b,k)123 rshift(b, k) _Bigint *b; int k;
124 #else
125 rshift(_Bigint *b, int k)
126 #endif
127 {
128 	__ULong *x, *x1, *xe, y;
129 	int n;
130 
131 	x = x1 = b->_x;
132 	n = k >> kshift;
133 	if (n < b->_wds) {
134 		xe = x + b->_wds;
135 		x += n;
136 		if (k &= kmask) {
137 			n = ULbits - k;
138 			y = *x++ >> k;
139 			while(x < xe) {
140 				*x1++ = (y | (*x << n)) & ALL_ON;
141 				y = *x++ >> k;
142 				}
143 			if ((*x1 = y) !=0)
144 				x1++;
145 			}
146 		else
147 			while(x < xe)
148 				*x1++ = *x++;
149 		}
150 	if ((b->_wds = x1 - b->_x) == 0)
151 		b->_x[0] = 0;
152 	}
153 
154 static int
155 #ifdef KR_headers
trailz(b)156 trailz(b) _Bigint *b;
157 #else
158 trailz(_Bigint *b)
159 #endif
160 {
161 	__ULong L, *x, *xe;
162 	int n = 0;
163 
164 	x = b->_x;
165 	xe = x + b->_wds;
166 	for(n = 0; x < xe && !*x; x++)
167 		n += ULbits;
168 	if (x < xe) {
169 		L = *x;
170 		n += lo0bits(&L);
171 		}
172 	return n;
173 	}
174 
175  _Bigint *
176 #ifdef KR_headers
increment(p,b)177 increment(p, b) struct _reent *p; _Bigint *b;
178 #else
179 increment(struct _reent *p, _Bigint *b)
180 #endif
181 {
182 	__ULong *x, *xe;
183 	_Bigint *b1;
184 #ifdef Pack_16
185 	__ULong carry = 1, y;
186 #endif
187 
188 	x = b->_x;
189 	xe = x + b->_wds;
190 #ifdef Pack_32
191 	do {
192 		if (*x < (__ULong)0xffffffffL) {
193 			++*x;
194 			return b;
195 			}
196 		*x++ = 0;
197 		} while(x < xe);
198 #else
199 	do {
200 		y = *x + carry;
201 		carry = y >> 16;
202 		*x++ = y & 0xffff;
203 		if (!carry)
204 			return b;
205 		} while(x < xe);
206 	if (carry)
207 #endif
208 	{
209 		if (b->_wds >= b->_maxwds) {
210 			b1 = Balloc(p,b->_k+1);
211 			Bcopy(b1,b);
212 			Bfree(p,b);
213 			b = b1;
214 			}
215 		b->_x[b->_wds++] = 1;
216 		}
217 	return b;
218 	}
219 
220  int
221 #ifdef KR_headers
decrement(b)222 decrement(b) _Bigint *b;
223 #else
224 decrement(_Bigint *b)
225 #endif
226 {
227 	__ULong *x, *xe;
228 #ifdef Pack_16
229 	__ULong borrow = 1, y;
230 #endif
231 
232 	x = b->_x;
233 	xe = x + b->_wds;
234 #ifdef Pack_32
235 	do {
236 		if (*x) {
237 			--*x;
238 			break;
239 			}
240 		*x++ = 0xffffffffL;
241 		}
242 		while(x < xe);
243 #else
244 	do {
245 		y = *x - borrow;
246 		borrow = (y & 0x10000) >> 16;
247 		*x++ = y & 0xffff;
248 		} while(borrow && x < xe);
249 #endif
250 	return STRTOG_Inexlo;
251 	}
252 
253  static int
254 #ifdef KR_headers
all_on(b,n)255 all_on(b, n) _Bigint *b; int n;
256 #else
257 all_on(_Bigint *b, int n)
258 #endif
259 {
260 	__ULong *x, *xe;
261 
262 	x = b->_x;
263 	xe = x + (n >> kshift);
264 	while(x < xe)
265 		if ((*x++ & ALL_ON) != ALL_ON)
266 			return 0;
267 	if (n &= kmask)
268 		return ((*x | (ALL_ON << n)) & ALL_ON) == ALL_ON;
269 	return 1;
270 	}
271 
272  _Bigint *
273 #ifdef KR_headers
set_ones(p,b,n)274 set_ones(p, b, n) struct _reent *p; _Bigint *b; int n;
275 #else
276 set_ones(struct _reent *p, _Bigint *b, int n)
277 #endif
278 {
279 	int k;
280 	__ULong *x, *xe;
281 
282 	k = (n + ((1 << kshift) - 1)) >> kshift;
283 	if (b->_k < k) {
284 		Bfree(p,b);
285 		b = Balloc(p,k);
286 		}
287 	k = n >> kshift;
288 	if (n &= kmask)
289 		k++;
290 	b->_wds = k;
291 	x = b->_x;
292 	xe = x + k;
293 	while(x < xe)
294 		*x++ = ALL_ON;
295 	if (n)
296 		x[-1] >>= ULbits - n;
297 	return b;
298 	}
299 
300  static int
rvOK(p,d,fpi,exp,bits,exact,rd,irv)301 rvOK
302 #ifdef KR_headers
303  (p, d, fpi, exp, bits, exact, rd, irv)
304  struct _reent *p; double d; FPI *fpi; Long *exp; __ULong *bits; int exact, rd, *irv;
305 #else
306  (struct _reent *p, double d, FPI *fpi, Long *exp, __ULong *bits, int exact, int rd, int *irv)
307 #endif
308 {
309 	_Bigint *b;
310 	__ULong carry, inex, lostbits;
311 	int bdif, e, j, k, k1, nb, rv;
312 
313 	carry = rv = 0;
314 	b = d2b(p, d, &e, &bdif);
315 	bdif -= nb = fpi->nbits;
316 	e += bdif;
317 	if (bdif <= 0) {
318 		if (exact)
319 			goto trunc;
320 		goto ret;
321 		}
322 	if (P == nb) {
323 		if (
324 #ifndef IMPRECISE_INEXACT
325 			exact &&
326 #endif
327 			fpi->rounding ==
328 #ifdef RND_PRODQUOT
329 					FPI_Round_near
330 #else
331 					Flt_Rounds
332 #endif
333 			) goto trunc;
334 		goto ret;
335 		}
336 	switch(rd) {
337 	  case 1:
338 		goto trunc;
339 	  case 2:
340 		break;
341 	  default: /* round near */
342 		k = bdif - 1;
343 		if (k < 0)
344 			goto trunc;
345 		if (!k) {
346 			if (!exact)
347 				goto ret;
348 			if (b->_x[0] & 2)
349 				break;
350 			goto trunc;
351 			}
352 		if (b->_x[k>>kshift] & ((__ULong)1 << (k & kmask)))
353 			break;
354 		goto trunc;
355 	  }
356 	/* "break" cases: round up 1 bit, then truncate; bdif > 0 */
357 	carry = 1;
358  trunc:
359 	inex = lostbits = 0;
360 	if (bdif > 0) {
361 		if ( (lostbits = any_on(b, bdif)) !=0)
362 			inex = STRTOG_Inexlo;
363 		rshift(b, bdif);
364 		if (carry) {
365 			inex = STRTOG_Inexhi;
366 			b = increment(p, b);
367 			if ( (j = nb & kmask) !=0)
368 				j = ULbits - j;
369 			if (hi0bits(b->_x[b->_wds - 1]) != j) {
370 				if (!lostbits)
371 					lostbits = b->_x[0] & 1;
372 				rshift(b, 1);
373 				e++;
374 				}
375 			}
376 		}
377 	else if (bdif < 0)
378 		b = lshift(p, b, -bdif);
379 	if (e < fpi->emin) {
380 		k = fpi->emin - e;
381 		e = fpi->emin;
382 		if (k > nb || fpi->sudden_underflow) {
383 			b->_wds = inex = 0;
384 			*irv = STRTOG_Underflow | STRTOG_Inexlo;
385 			}
386 		else {
387 			k1 = k - 1;
388 			if (k1 > 0 && !lostbits)
389 				lostbits = any_on(b, k1);
390 			if (!lostbits && !exact)
391 				goto ret;
392 			lostbits |=
393 			  carry = b->_x[k1>>kshift] & (1 << (k1 & kmask));
394 			rshift(b, k);
395 			*irv = STRTOG_Denormal;
396 			if (carry) {
397 				b = increment(p, b);
398 				inex = STRTOG_Inexhi | STRTOG_Underflow;
399 				}
400 			else if (lostbits)
401 				inex = STRTOG_Inexlo | STRTOG_Underflow;
402 			}
403 		}
404 	else if (e > fpi->emax) {
405 		e = fpi->emax + 1;
406 		*irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
407 #ifndef NO_ERRNO
408 		errno = ERANGE;
409 #endif
410 		b->_wds = inex = 0;
411 		}
412 	*exp = e;
413 	copybits(bits, nb, b);
414 	*irv |= inex;
415 	rv = 1;
416  ret:
417 	Bfree(p,b);
418 	return rv;
419 	}
420 
421  static int
422 #ifdef KR_headers
mantbits(d)423 mantbits(d) double d;
424 #else
425 mantbits(U d)
426 #endif
427 {
428 	__ULong L;
429 #ifdef VAX
430 	L = word1(d) << 16 | word1(d) >> 16;
431 	if (L)
432 #else
433 	if ( (L = word1(d)) !=0)
434 #endif
435 		return P - lo0bits(&L);
436 #ifdef VAX
437 	L = word0(d) << 16 | word0(d) >> 16 | Exp_msk11;
438 #else
439 	L = word0(d) | Exp_msk1;
440 #endif
441 	return P - 32 - lo0bits(&L);
442 	}
443 
444  int
_strtodg_r(p,s00,se,fpi,exp,bits)445 _strtodg_r
446 #ifdef KR_headers
447 	(p, s00, se, fpi, exp, bits)
448 	struct _reent *p; const char *s00; char **se; FPI *fpi; Long *exp; __ULong *bits;
449 #else
450 	(struct _reent *p, const char *s00, char **se, FPI *fpi, Long *exp, __ULong *bits)
451 #endif
452 {
453 	int abe, abits, asub;
454 	int bb0, bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, denorm;
455 	int dsign, e, e1, e2, emin, esign, finished, i, inex, irv;
456 	int j, k, nbits, nd, nd0, nf, nz, nz0, rd, rvbits, rve, rve1, sign;
457 	int sudden_underflow;
458 	const char *s, *s0, *s1;
459 	//double adj, adj0, rv, tol;
460 	double adj0, tol;
461 	U adj, rv;
462 	Long L;
463 	__ULong y, z;
464 	_Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0;
465 
466 	irv = STRTOG_Zero;
467 	denorm = sign = nz0 = nz = 0;
468 	dval(rv) = 0.;
469 	rvb = 0;
470 	nbits = fpi->nbits;
471 	for(s = s00;;s++) switch(*s) {
472 		case '-':
473 			sign = 1;
474 			/* no break */
475 		case '+':
476 			if (*++s)
477 				goto break2;
478 			/* no break */
479 		case 0:
480 			sign = 0;
481 			irv = STRTOG_NoNumber;
482 			s = s00;
483 			goto ret;
484 		case '\t':
485 		case '\n':
486 		case '\v':
487 		case '\f':
488 		case '\r':
489 		case ' ':
490 			continue;
491 		default:
492 			goto break2;
493 		}
494  break2:
495 	if (*s == '0') {
496 #ifndef NO_HEX_FP
497 		switch(s[1]) {
498 		  case 'x':
499 		  case 'X':
500 			irv = gethex(p, &s, fpi, exp, &rvb, sign);
501 			if (irv == STRTOG_NoNumber) {
502 				s = s00;
503 				sign = 0;
504 				}
505 			goto ret;
506 		  }
507 #endif
508 		nz0 = 1;
509 		while(*++s == '0') ;
510 		if (!*s)
511 			goto ret;
512 		}
513 	sudden_underflow = fpi->sudden_underflow;
514 	s0 = s;
515 	y = z = 0;
516 	for(decpt = nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
517 		if (nd < 9)
518 			y = 10*y + c - '0';
519 		else if (nd < 16)
520 			z = 10*z + c - '0';
521 	nd0 = nd;
522 #ifdef USE_LOCALE
523 	if (strncmp (s, _localeconv_r (p)->decimal_point,
524 		     strlen (_localeconv_r (p)->decimal_point)) == 0)
525 #else
526 	if (c == '.')
527 #endif
528 		{
529 		decpt = 1;
530 #ifdef USE_LOCALE
531 		c = *(s += strlen (_localeconv_r (p)->decimal_point));
532 #else
533 		c = *++s;
534 #endif
535 		if (!nd) {
536 			for(; c == '0'; c = *++s)
537 				nz++;
538 			if (c > '0' && c <= '9') {
539 				s0 = s;
540 				nf += nz;
541 				nz = 0;
542 				goto have_dig;
543 				}
544 			goto dig_done;
545 			}
546 		for(; c >= '0' && c <= '9'; c = *++s) {
547  have_dig:
548 			nz++;
549 			if (c -= '0') {
550 				nf += nz;
551 				for(i = 1; i < nz; i++)
552 					if (nd++ < 9)
553 						y *= 10;
554 					else if (nd <= DBL_DIG + 1)
555 						z *= 10;
556 				if (nd++ < 9)
557 					y = 10*y + c;
558 				else if (nd <= DBL_DIG + 1)
559 					z = 10*z + c;
560 				nz = 0;
561 				}
562 			}
563 		}
564  dig_done:
565 	e = 0;
566 	if (c == 'e' || c == 'E') {
567 		if (!nd && !nz && !nz0) {
568 			irv = STRTOG_NoNumber;
569 			s = s00;
570 			goto ret;
571 			}
572 		s00 = s;
573 		esign = 0;
574 		switch(c = *++s) {
575 			case '-':
576 				esign = 1;
577 			case '+':
578 				c = *++s;
579 			}
580 		if (c >= '0' && c <= '9') {
581 			while(c == '0')
582 				c = *++s;
583 			if (c > '0' && c <= '9') {
584 				L = c - '0';
585 				s1 = s;
586 				while((c = *++s) >= '0' && c <= '9')
587 					L = 10*L + c - '0';
588 				if (s - s1 > 8 || L > 19999)
589 					/* Avoid confusion from exponents
590 					 * so large that e might overflow.
591 					 */
592 					e = 19999; /* safe for 16 bit ints */
593 				else
594 					e = (int)L;
595 				if (esign)
596 					e = -e;
597 				}
598 			else
599 				e = 0;
600 			}
601 		else
602 			s = s00;
603 		}
604 	if (!nd) {
605 		if (!nz && !nz0) {
606 #ifdef INFNAN_CHECK
607 			/* Check for Nan and Infinity */
608 			if (!decpt)
609 			 switch(c) {
610 			  case 'i':
611 			  case 'I':
612 				if (match(&s,"nf")) {
613 					--s;
614 					if (!match(&s,"inity"))
615 						++s;
616 					irv = STRTOG_Infinite;
617 					goto infnanexp;
618 					}
619 				break;
620 			  case 'n':
621 			  case 'N':
622 				if (match(&s, "an")) {
623 					irv = STRTOG_NaN;
624 					*exp = fpi->emax + 1;
625 #ifndef No_Hex_NaN
626 					if (*s == '(') /*)*/
627 						irv = hexnan(&s, fpi, bits);
628 #endif
629 					goto infnanexp;
630 					}
631 			  }
632 #endif /* INFNAN_CHECK */
633 			irv = STRTOG_NoNumber;
634 			s = s00;
635 			}
636 		goto ret;
637 		}
638 
639 	irv = STRTOG_Normal;
640 	e1 = e -= nf;
641 	rd = 0;
642 	switch(fpi->rounding & 3) {
643 	  case FPI_Round_up:
644 		rd = 2 - sign;
645 		break;
646 	  case FPI_Round_zero:
647 		rd = 1;
648 		break;
649 	  case FPI_Round_down:
650 		rd = 1 + sign;
651 	  }
652 
653 	/* Now we have nd0 digits, starting at s0, followed by a
654 	 * decimal point, followed by nd-nd0 digits.  The number we're
655 	 * after is the integer represented by those digits times
656 	 * 10**e */
657 
658 	if (!nd0)
659 		nd0 = nd;
660 	k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
661 	dval(rv) = y;
662 	if (k > 9)
663 		dval(rv) = tens[k - 9] * dval(rv) + z;
664 	bd0 = 0;
665 	if (nbits <= P && nd <= DBL_DIG) {
666 		if (!e) {
667 			if (rvOK(p, dval(rv), fpi, exp, bits, 1, rd, &irv))
668 				goto ret;
669 			}
670 		else if (e > 0) {
671 			if (e <= Ten_pmax) {
672 #ifdef VAX
673 				goto vax_ovfl_check;
674 #else
675 				i = fivesbits[e] + mantbits(rv) <= P;
676 				/* rv = */ rounded_product(dval(rv), tens[e]);
677 				if (rvOK(p, dval(rv), fpi, exp, bits, i, rd, &irv))
678 					goto ret;
679 				e1 -= e;
680 				goto rv_notOK;
681 #endif
682 				}
683 			i = DBL_DIG - nd;
684 			if (e <= Ten_pmax + i) {
685 				/* A fancier test would sometimes let us do
686 				 * this for larger i values.
687 				 */
688 				e2 = e - i;
689 				e1 -= i;
690 				dval(rv) *= tens[i];
691 #ifdef VAX
692 				/* VAX exponent range is so narrow we must
693 				 * worry about overflow here...
694 				 */
695  vax_ovfl_check:
696 				dval(adj) = dval(rv);
697 				word0(adj) -= P*Exp_msk1;
698 				/* adj = */ rounded_product(dval(adj), tens[e2]);
699 				if ((word0(adj) & Exp_mask)
700 				 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
701 					goto rv_notOK;
702 				word0(adj) += P*Exp_msk1;
703 				dval(rv) = dval(adj);
704 #else
705 				/* rv = */ rounded_product(dval(rv), tens[e2]);
706 #endif
707 				if (rvOK(p, dval(rv), fpi, exp, bits, 0, rd, &irv))
708 					goto ret;
709 				e1 -= e2;
710 				}
711 			}
712 #ifndef Inaccurate_Divide
713 		else if (e >= -Ten_pmax) {
714 			/* rv = */ rounded_quotient(dval(rv), tens[-e]);
715 			if (rvOK(p, dval(rv), fpi, exp, bits, 0, rd, &irv))
716 				goto ret;
717 			e1 -= e;
718 			}
719 #endif
720 		}
721  rv_notOK:
722 	e1 += nd - k;
723 
724 	/* Get starting approximation = rv * 10**e1 */
725 
726 	e2 = 0;
727 	if (e1 > 0) {
728 		if ( (i = e1 & 15) !=0)
729 			dval(rv) *= tens[i];
730 		if (e1 &= ~15) {
731 			e1 >>= 4;
732 			while(e1 >= (1 << n_bigtens-1)) {
733 				e2 += ((word0(rv) & Exp_mask)
734 					>> Exp_shift1) - Bias;
735 				word0(rv) &= ~Exp_mask;
736 				word0(rv) |= Bias << Exp_shift1;
737 				dval(rv) *= bigtens[n_bigtens-1];
738 				e1 -= 1 << n_bigtens-1;
739 				}
740 			e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias;
741 			word0(rv) &= ~Exp_mask;
742 			word0(rv) |= Bias << Exp_shift1;
743 			for(j = 0; e1 > 0; j++, e1 >>= 1)
744 				if (e1 & 1)
745 					dval(rv) *= bigtens[j];
746 			}
747 		}
748 	else if (e1 < 0) {
749 		e1 = -e1;
750 		if ( (i = e1 & 15) !=0)
751 			dval(rv) /= tens[i];
752 		if (e1 &= ~15) {
753 			e1 >>= 4;
754 			while(e1 >= (1 << n_bigtens-1)) {
755 				e2 += ((word0(rv) & Exp_mask)
756 					>> Exp_shift1) - Bias;
757 				word0(rv) &= ~Exp_mask;
758 				word0(rv) |= Bias << Exp_shift1;
759 				dval(rv) *= tinytens[n_bigtens-1];
760 				e1 -= 1 << n_bigtens-1;
761 				}
762 			e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias;
763 			word0(rv) &= ~Exp_mask;
764 			word0(rv) |= Bias << Exp_shift1;
765 			for(j = 0; e1 > 0; j++, e1 >>= 1)
766 				if (e1 & 1)
767 					dval(rv) *= tinytens[j];
768 			}
769 		}
770 #ifdef IBM
771 	/* e2 is a correction to the (base 2) exponent of the return
772 	 * value, reflecting adjustments above to avoid overflow in the
773 	 * native arithmetic.  For native IBM (base 16) arithmetic, we
774 	 * must multiply e2 by 4 to change from base 16 to 2.
775 	 */
776 	e2 <<= 2;
777 #endif
778 	rvb = d2b(p, dval(rv), &rve, &rvbits);	/* rv = rvb * 2^rve */
779 	rve += e2;
780 	if ((j = rvbits - nbits) > 0) {
781 		rshift(rvb, j);
782 		rvbits = nbits;
783 		rve += j;
784 		}
785 	bb0 = 0;	/* trailing zero bits in rvb */
786 	e2 = rve + rvbits - nbits;
787 	if (e2 > fpi->emax + 1)
788 		goto huge;
789 	rve1 = rve + rvbits - nbits;
790 	if (e2 < (emin = fpi->emin)) {
791 		denorm = 1;
792 		j = rve - emin;
793 		if (j > 0) {
794 			rvb = lshift(p, rvb, j);
795 			rvbits += j;
796 			}
797 		else if (j < 0) {
798 			rvbits += j;
799 			if (rvbits <= 0) {
800 				if (rvbits < -1) {
801  ufl:
802 					rvb->_wds = 0;
803 					rvb->_x[0] = 0;
804 					*exp = emin;
805 					irv = STRTOG_Underflow | STRTOG_Inexlo;
806 					goto ret;
807 					}
808 				rvb->_x[0] = rvb->_wds = rvbits = 1;
809 				}
810 			else
811 				rshift(rvb, -j);
812 			}
813 		rve = rve1 = emin;
814 		if (sudden_underflow && e2 + 1 < emin)
815 			goto ufl;
816 		}
817 
818 	/* Now the hard part -- adjusting rv to the correct value.*/
819 
820 	/* Put digits into bd: true value = bd * 10^e */
821 
822 	bd0 = s2b(p, s0, nd0, nd, y);
823 
824 	for(;;) {
825 		bd = Balloc(p,bd0->_k);
826 		Bcopy(bd, bd0);
827 		bb = Balloc(p,rvb->_k);
828 		Bcopy(bb, rvb);
829 		bbbits = rvbits - bb0;
830 		bbe = rve + bb0;
831 		bs = i2b(p, 1);
832 
833 		if (e >= 0) {
834 			bb2 = bb5 = 0;
835 			bd2 = bd5 = e;
836 			}
837 		else {
838 			bb2 = bb5 = -e;
839 			bd2 = bd5 = 0;
840 			}
841 		if (bbe >= 0)
842 			bb2 += bbe;
843 		else
844 			bd2 -= bbe;
845 		bs2 = bb2;
846 		j = nbits + 1 - bbbits;
847 		i = bbe + bbbits - nbits;
848 		if (i < emin)	/* denormal */
849 			j += i - emin;
850 		bb2 += j;
851 		bd2 += j;
852 		i = bb2 < bd2 ? bb2 : bd2;
853 		if (i > bs2)
854 			i = bs2;
855 		if (i > 0) {
856 			bb2 -= i;
857 			bd2 -= i;
858 			bs2 -= i;
859 			}
860 		if (bb5 > 0) {
861 			bs = pow5mult(p, bs, bb5);
862 			bb1 = mult(p, bs, bb);
863 			Bfree(p,bb);
864 			bb = bb1;
865 			}
866 		bb2 -= bb0;
867 		if (bb2 > 0)
868 			bb = lshift(p, bb, bb2);
869 		else if (bb2 < 0)
870 			rshift(bb, -bb2);
871 		if (bd5 > 0)
872 			bd = pow5mult(p, bd, bd5);
873 		if (bd2 > 0)
874 			bd = lshift(p, bd, bd2);
875 		if (bs2 > 0)
876 			bs = lshift(p, bs, bs2);
877 		asub = 1;
878 		inex = STRTOG_Inexhi;
879 		delta = diff(p, bb, bd);
880 		if (delta->_wds <= 1 && !delta->_x[0])
881 			break;
882 		dsign = delta->_sign;
883 		delta->_sign = finished = 0;
884 		L = 0;
885 		i = cmp(delta, bs);
886 		if (rd && i <= 0) {
887 			irv = STRTOG_Normal;
888 			if ( (finished = dsign ^ (rd&1)) !=0) {
889 				if (dsign != 0) {
890 					irv |= STRTOG_Inexhi;
891 					goto adj1;
892 					}
893 				irv |= STRTOG_Inexlo;
894 				if (rve1 == emin)
895 					goto adj1;
896 				for(i = 0, j = nbits; j >= ULbits;
897 						i++, j -= ULbits) {
898 					if (rvb->_x[i] & ALL_ON)
899 						goto adj1;
900 					}
901 				if (j > 1 && lo0bits(rvb->_x + i) < j - 1)
902 					goto adj1;
903 				rve = rve1 - 1;
904 				rvb = set_ones(p, rvb, rvbits = nbits);
905 				break;
906 				}
907 			irv |= dsign ? STRTOG_Inexlo : STRTOG_Inexhi;
908 			break;
909 			}
910 		if (i < 0) {
911 			/* Error is less than half an ulp -- check for
912 			 * special case of mantissa a power of two.
913 			 */
914 			irv = dsign
915 				? STRTOG_Normal | STRTOG_Inexlo
916 				: STRTOG_Normal | STRTOG_Inexhi;
917 			if (dsign || bbbits > 1 || denorm || rve1 == emin)
918 				break;
919 			delta = lshift(p, delta,1);
920 			if (cmp(delta, bs) > 0) {
921 				irv = STRTOG_Normal | STRTOG_Inexlo;
922 				goto drop_down;
923 				}
924 			break;
925 			}
926 		if (i == 0) {
927 			/* exactly half-way between */
928 			if (dsign) {
929 				if (denorm && all_on(rvb, rvbits)) {
930 					/*boundary case -- increment exponent*/
931 					rvb->_wds = 1;
932 					rvb->_x[0] = 1;
933 					rve = emin + nbits - (rvbits = 1);
934 					irv = STRTOG_Normal | STRTOG_Inexhi;
935 					denorm = 0;
936 					break;
937 					}
938 				irv = STRTOG_Normal | STRTOG_Inexlo;
939 				}
940 			else if (bbbits == 1) {
941 				irv = STRTOG_Normal;
942  drop_down:
943 				/* boundary case -- decrement exponent */
944 				if (rve1 == emin) {
945 					irv = STRTOG_Normal | STRTOG_Inexhi;
946 					if (rvb->_wds == 1 && rvb->_x[0] == 1)
947 						sudden_underflow = 1;
948 					break;
949 					}
950 				rve -= nbits;
951 				rvb = set_ones(p, rvb, rvbits = nbits);
952 				break;
953 				}
954 			else
955 				irv = STRTOG_Normal | STRTOG_Inexhi;
956 			if (bbbits < nbits && !denorm || !(rvb->_x[0] & 1))
957 				break;
958 			if (dsign) {
959 				rvb = increment(p, rvb);
960 				j = kmask & (ULbits - (rvbits & kmask));
961 				if (hi0bits(rvb->_x[rvb->_wds - 1]) != j)
962 					rvbits++;
963 				irv = STRTOG_Normal | STRTOG_Inexhi;
964 				}
965 			else {
966 				if (bbbits == 1)
967 					goto undfl;
968 				decrement(rvb);
969 				irv = STRTOG_Normal | STRTOG_Inexlo;
970 				}
971 			break;
972 			}
973 		if ((dval(adj) = ratio(delta, bs)) <= 2.) {
974  adj1:
975 			inex = STRTOG_Inexlo;
976 			if (dsign) {
977 				asub = 0;
978 				inex = STRTOG_Inexhi;
979 				}
980 			else if (denorm && bbbits <= 1) {
981  undfl:
982 				rvb->_wds = 0;
983 				rve = emin;
984 				irv = STRTOG_Underflow | STRTOG_Inexlo;
985 				break;
986 				}
987 			adj0 = dval(adj) = 1.;
988 			}
989 		else {
990 			adj0 = dval(adj) *= 0.5;
991 			if (dsign) {
992 				asub = 0;
993 				inex = STRTOG_Inexlo;
994 				}
995 			if (dval(adj) < 2147483647.) {
996 				L = adj0;
997 				adj0 -= L;
998 				switch(rd) {
999 				  case 0:
1000 					if (adj0 >= .5)
1001 						goto inc_L;
1002 					break;
1003 				  case 1:
1004 					if (asub && adj0 > 0.)
1005 						goto inc_L;
1006 					break;
1007 				  case 2:
1008 					if (!asub && adj0 > 0.) {
1009  inc_L:
1010 						L++;
1011 						inex = STRTOG_Inexact - inex;
1012 						}
1013 				  }
1014 				dval(adj) = L;
1015 				}
1016 			}
1017 		y = rve + rvbits;
1018 
1019 		/* adj *= ulp(dval(rv)); */
1020 		/* if (asub) rv -= adj; else rv += adj; */
1021 
1022 		if (!denorm && rvbits < nbits) {
1023 			rvb = lshift(p, rvb, j = nbits - rvbits);
1024 			rve -= j;
1025 			rvbits = nbits;
1026 			}
1027 		ab = d2b(p, dval(adj), &abe, &abits);
1028 		if (abe < 0)
1029 			rshift(ab, -abe);
1030 		else if (abe > 0)
1031 			ab = lshift(p, ab, abe);
1032 		rvb0 = rvb;
1033 		if (asub) {
1034 			/* rv -= adj; */
1035 			j = hi0bits(rvb->_x[rvb->_wds-1]);
1036 			rvb = diff(p, rvb, ab);
1037 			k = rvb0->_wds - 1;
1038 			if (denorm)
1039 				/* do nothing */;
1040 			else if (rvb->_wds <= k
1041 				|| hi0bits( rvb->_x[k]) >
1042 				   hi0bits(rvb0->_x[k])) {
1043 				/* unlikely; can only have lost 1 high bit */
1044 				if (rve1 == emin) {
1045 					--rvbits;
1046 					denorm = 1;
1047 					}
1048 				else {
1049 					rvb = lshift(p, rvb, 1);
1050 					--rve;
1051 					--rve1;
1052 					L = finished = 0;
1053 					}
1054 				}
1055 			}
1056 		else {
1057 			rvb = sum(p, rvb, ab);
1058 			k = rvb->_wds - 1;
1059 			if (k >= rvb0->_wds
1060 			 || hi0bits(rvb->_x[k]) < hi0bits(rvb0->_x[k])) {
1061 				if (denorm) {
1062 					if (++rvbits == nbits)
1063 						denorm = 0;
1064 					}
1065 				else {
1066 					rshift(rvb, 1);
1067 					rve++;
1068 					rve1++;
1069 					L = 0;
1070 					}
1071 				}
1072 			}
1073 		Bfree(p,ab);
1074 		Bfree(p,rvb0);
1075 		if (finished)
1076 			break;
1077 
1078 		z = rve + rvbits;
1079 		if (y == z && L) {
1080 			/* Can we stop now? */
1081 			tol = dval(adj) * 5e-16; /* > max rel error */
1082 			dval(adj) = adj0 - .5;
1083 			if (dval(adj) < -tol) {
1084 				if (adj0 > tol) {
1085 					irv |= inex;
1086 					break;
1087 					}
1088 				}
1089 			else if (dval(adj) > tol && adj0 < 1. - tol) {
1090 				irv |= inex;
1091 				break;
1092 				}
1093 			}
1094 		bb0 = denorm ? 0 : trailz(rvb);
1095 		Bfree(p,bb);
1096 		Bfree(p,bd);
1097 		Bfree(p,bs);
1098 		Bfree(p,delta);
1099 		}
1100 	if (!denorm && (j = nbits - rvbits)) {
1101 		if (j > 0)
1102 			rvb = lshift(p, rvb, j);
1103 		else
1104 			rshift(rvb, -j);
1105 		rve -= j;
1106 		}
1107 	*exp = rve;
1108 	Bfree(p,bb);
1109 	Bfree(p,bd);
1110 	Bfree(p,bs);
1111 	Bfree(p,bd0);
1112 	Bfree(p,delta);
1113 	if (rve > fpi->emax) {
1114  huge:
1115 		rvb->_wds = 0;
1116 		irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
1117 #ifndef NO_ERRNO
1118 		errno = ERANGE;
1119 #endif
1120  infnanexp:
1121 		*exp = fpi->emax + 1;
1122 		}
1123  ret:
1124 	if (denorm) {
1125 		if (sudden_underflow) {
1126 			rvb->_wds = 0;
1127 			irv = STRTOG_Underflow | STRTOG_Inexlo;
1128 			}
1129 		else  {
1130 			irv = (irv & ~STRTOG_Retmask) |
1131 				(rvb->_wds > 0 ? STRTOG_Denormal : STRTOG_Zero);
1132 			if (irv & STRTOG_Inexact)
1133 				irv |= STRTOG_Underflow;
1134 			}
1135 		}
1136 	if (se)
1137 		*se = (char *)s;
1138 	if (sign)
1139 		irv |= STRTOG_Neg;
1140 	if (rvb) {
1141 		copybits(bits, nbits, rvb);
1142 		Bfree(p,rvb);
1143 		}
1144 	return irv;
1145 	}
1146 
1147 #endif /* _HAVE_LONG_DOUBLE && !_LDBL_EQ_DBL */
1148