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