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