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