1/****************************************************************
2
3The author of this software is David M. Gay.
4
5Copyright (C) 1998, 1999 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 static Bigint *
35#ifdef KR_headers
36bitstob(bits, nbits, bbits) ULong *bits; int nbits; int *bbits;
37#else
38bitstob(ULong *bits, int nbits, int *bbits)
39#endif
40{
41	int i, k;
42	Bigint *b;
43	ULong *be, *x, *x0;
44
45	i = ULbits;
46	k = 0;
47	while(i < nbits) {
48		i <<= 1;
49		k++;
50		}
51#ifndef Pack_32
52	if (!k)
53		k = 1;
54#endif
55	b = Balloc(k);
56	if (b == NULL)
57		return (NULL);
58	be = bits + ((nbits - 1) >> kshift);
59	x = x0 = b->x;
60	do {
61		*x++ = *bits & ALL_ON;
62#ifdef Pack_16
63		*x++ = (*bits >> 16) & ALL_ON;
64#endif
65		} while(++bits <= be);
66	i = x - x0;
67	while(!x0[--i])
68		if (!i) {
69			b->wds = 0;
70			*bbits = 0;
71			goto ret;
72			}
73	b->wds = i + 1;
74	*bbits = i*ULbits + 32 - hi0bits(b->x[i]);
75 ret:
76	return b;
77	}
78
79/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
80 *
81 * Inspired by "How to Print Floating-Point Numbers Accurately" by
82 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
83 *
84 * Modifications:
85 *	1. Rather than iterating, we use a simple numeric overestimate
86 *	   to determine k = floor(log10(d)).  We scale relevant
87 *	   quantities using O(log2(k)) rather than O(k) multiplications.
88 *	2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
89 *	   try to generate digits strictly left to right.  Instead, we
90 *	   compute with fewer bits and propagate the carry if necessary
91 *	   when rounding the final digit up.  This is often faster.
92 *	3. Under the assumption that input will be rounded nearest,
93 *	   mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
94 *	   That is, we allow equality in stopping tests when the
95 *	   round-nearest rule will give the same floating-point value
96 *	   as would satisfaction of the stopping test with strict
97 *	   inequality.
98 *	4. We remove common factors of powers of 2 from relevant
99 *	   quantities.
100 *	5. When converting floating-point integers less than 1e16,
101 *	   we use floating-point arithmetic rather than resorting
102 *	   to multiple-precision integers.
103 *	6. When asked to produce fewer than 15 digits, we first try
104 *	   to get by with floating-point arithmetic; we resort to
105 *	   multiple-precision integer arithmetic only if we cannot
106 *	   guarantee that the floating-point calculation has given
107 *	   the correctly rounded result.  For k requested digits and
108 *	   "uniformly" distributed input, the probability is
109 *	   something like 10^(k-15) that we must resort to the Long
110 *	   calculation.
111 */
112
113 char *
114gdtoa
115#ifdef KR_headers
116	(fpi, be, bits, kindp, mode, ndigits, decpt, rve)
117	FPI *fpi; int be; ULong *bits;
118	int *kindp, mode, ndigits, *decpt; char **rve;
119#else
120	(FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, int *decpt, char **rve)
121#endif
122{
123 /*	Arguments ndigits and decpt are similar to the second and third
124	arguments of ecvt and fcvt; trailing zeros are suppressed from
125	the returned string.  If not null, *rve is set to point
126	to the end of the return value.  If d is +-Infinity or NaN,
127	then *decpt is set to 9999.
128	be = exponent: value = (integer represented by bits) * (2 to the power of be).
129
130	mode:
131		0 ==> shortest string that yields d when read in
132			and rounded to nearest.
133		1 ==> like 0, but with Steele & White stopping rule;
134			e.g. with IEEE P754 arithmetic , mode 0 gives
135			1e23 whereas mode 1 gives 9.999999999999999e22.
136		2 ==> max(1,ndigits) significant digits.  This gives a
137			return value similar to that of ecvt, except
138			that trailing zeros are suppressed.
139		3 ==> through ndigits past the decimal point.  This
140			gives a return value similar to that from fcvt,
141			except that trailing zeros are suppressed, and
142			ndigits can be negative.
143		4-9 should give the same return values as 2-3, i.e.,
144			4 <= mode <= 9 ==> same return as mode
145			2 + (mode & 1).  These modes are mainly for
146			debugging; often they run slower but sometimes
147			faster than modes 2-3.
148		4,5,8,9 ==> left-to-right digit generation.
149		6-9 ==> don't try fast floating-point estimate
150			(if applicable).
151
152		Values of mode other than 0-9 are treated as mode 0.
153
154		Sufficient space is allocated to the return value
155		to hold the suppressed trailing zeros.
156	*/
157
158	int bbits, b2, b5, be0, dig, i, ieps, ilim, ilim0, ilim1, inex;
159	int j, j1, k, k0, k_check, kind, leftright, m2, m5, nbits;
160	int rdir, s2, s5, spec_case, try_quick;
161	Long L;
162	Bigint *b, *b1, *delta, *mlo, *mhi, *mhi1, *S;
163	double d2, ds;
164	char *s, *s0;
165	U d, eps;
166
167#ifndef MULTIPLE_THREADS
168	if (dtoa_result) {
169		freedtoa(dtoa_result);
170		dtoa_result = 0;
171		}
172#endif
173	inex = 0;
174	kind = *kindp &= ~STRTOG_Inexact;
175	switch(kind & STRTOG_Retmask) {
176	  case STRTOG_Zero:
177		goto ret_zero;
178	  case STRTOG_Normal:
179	  case STRTOG_Denormal:
180		break;
181	  case STRTOG_Infinite:
182		*decpt = -32768;
183		return nrv_alloc("Infinity", rve, 8);
184	  case STRTOG_NaN:
185		*decpt = -32768;
186		return nrv_alloc("NaN", rve, 3);
187	  default:
188		return 0;
189	  }
190	b = bitstob(bits, nbits = fpi->nbits, &bbits);
191	if (b == NULL)
192		return (NULL);
193	be0 = be;
194	if ( (i = trailz(b)) !=0) {
195		rshift(b, i);
196		be += i;
197		bbits -= i;
198		}
199	if (!b->wds) {
200		Bfree(b);
201 ret_zero:
202		*decpt = 1;
203		return nrv_alloc("0", rve, 1);
204		}
205
206	dval(&d) = b2d(b, &i);
207	i = be + bbits - 1;
208	word0(&d) &= Frac_mask1;
209	word0(&d) |= Exp_11;
210#ifdef IBM
211	if ( (j = 11 - hi0bits(word0(&d) & Frac_mask)) !=0)
212		dval(&d) /= 1 << j;
213#endif
214
215	/* log(x)	~=~ log(1.5) + (x-1.5)/1.5
216	 * log10(x)	 =  log(x) / log(10)
217	 *		~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
218	 * log10(&d) = (i-Bias)*log(2)/log(10) + log10(d2)
219	 *
220	 * This suggests computing an approximation k to log10(&d) by
221	 *
222	 * k = (i - Bias)*0.301029995663981
223	 *	+ ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
224	 *
225	 * We want k to be too large rather than too small.
226	 * The error in the first-order Taylor series approximation
227	 * is in our favor, so we just round up the constant enough
228	 * to compensate for any error in the multiplication of
229	 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
230	 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
231	 * adding 1e-13 to the constant term more than suffices.
232	 * Hence we adjust the constant term to 0.1760912590558.
233	 * (We could get a more accurate k by invoking log10,
234	 *  but this is probably not worthwhile.)
235	 */
236#ifdef IBM
237	i <<= 2;
238	i += j;
239#endif
240	ds = (dval(&d)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
241
242	/* correct assumption about exponent range */
243	if ((j = i) < 0)
244		j = -j;
245	if ((j -= 1077) > 0)
246		ds += j * 7e-17;
247
248	k = (int)ds;
249	if (ds < 0. && ds != k)
250		k--;	/* want k = floor(ds) */
251	k_check = 1;
252#ifdef IBM
253	j = be + bbits - 1;
254	if ( (j1 = j & 3) !=0)
255		dval(&d) *= 1 << j1;
256	word0(&d) += j << Exp_shift - 2 & Exp_mask;
257#else
258	word0(&d) += (be + bbits - 1) << Exp_shift;
259#endif
260	if (k >= 0 && k <= Ten_pmax) {
261		if (dval(&d) < tens[k])
262			k--;
263		k_check = 0;
264		}
265	j = bbits - i - 1;
266	if (j >= 0) {
267		b2 = 0;
268		s2 = j;
269		}
270	else {
271		b2 = -j;
272		s2 = 0;
273		}
274	if (k >= 0) {
275		b5 = 0;
276		s5 = k;
277		s2 += k;
278		}
279	else {
280		b2 -= k;
281		b5 = -k;
282		s5 = 0;
283		}
284	if (mode < 0 || mode > 9)
285		mode = 0;
286	try_quick = 1;
287	if (mode > 5) {
288		mode -= 4;
289		try_quick = 0;
290		}
291	else if (i >= -4 - Emin || i < Emin)
292		try_quick = 0;
293	leftright = 1;
294	ilim = ilim1 = -1;	/* Values for cases 0 and 1; done here to */
295				/* silence erroneous "gcc -Wall" warning. */
296	switch(mode) {
297		case 0:
298		case 1:
299			i = (int)(nbits * .30103) + 3;
300			ndigits = 0;
301			break;
302		case 2:
303			leftright = 0;
304			/* no break */
305		case 4:
306			if (ndigits <= 0)
307				ndigits = 1;
308			ilim = ilim1 = i = ndigits;
309			break;
310		case 3:
311			leftright = 0;
312			/* no break */
313		case 5:
314			i = ndigits + k + 1;
315			ilim = i;
316			ilim1 = i - 1;
317			if (i <= 0)
318				i = 1;
319		}
320	s = s0 = rv_alloc(i);
321	if (s == NULL)
322		return (NULL);
323
324	if ( (rdir = fpi->rounding - 1) !=0) {
325		if (rdir < 0)
326			rdir = 2;
327		if (kind & STRTOG_Neg)
328			rdir = 3 - rdir;
329		}
330
331	/* Now rdir = 0 ==> round near, 1 ==> round up, 2 ==> round down. */
332
333	if (ilim >= 0 && ilim <= Quick_max && try_quick && !rdir
334#ifndef IMPRECISE_INEXACT
335		&& k == 0
336#endif
337								) {
338
339		/* Try to get by with floating-point arithmetic. */
340
341		i = 0;
342		d2 = dval(&d);
343#ifdef IBM
344		if ( (j = 11 - hi0bits(word0(&d) & Frac_mask)) !=0)
345			dval(&d) /= 1 << j;
346#endif
347		k0 = k;
348		ilim0 = ilim;
349		ieps = 2; /* conservative */
350		if (k > 0) {
351			ds = tens[k&0xf];
352			j = k >> 4;
353			if (j & Bletch) {
354				/* prevent overflows */
355				j &= Bletch - 1;
356				dval(&d) /= bigtens[n_bigtens-1];
357				ieps++;
358				}
359			for(; j; j >>= 1, i++)
360				if (j & 1) {
361					ieps++;
362					ds *= bigtens[i];
363					}
364			}
365		else  {
366			ds = 1.;
367			if ( (j1 = -k) !=0) {
368				dval(&d) *= tens[j1 & 0xf];
369				for(j = j1 >> 4; j; j >>= 1, i++)
370					if (j & 1) {
371						ieps++;
372						dval(&d) *= bigtens[i];
373						}
374				}
375			}
376		if (k_check && dval(&d) < 1. && ilim > 0) {
377			if (ilim1 <= 0)
378				goto fast_failed;
379			ilim = ilim1;
380			k--;
381			dval(&d) *= 10.;
382			ieps++;
383			}
384		dval(&eps) = ieps*dval(&d) + 7.;
385		word0(&eps) -= (P-1)*Exp_msk1;
386		if (ilim == 0) {
387			S = mhi = 0;
388			dval(&d) -= 5.;
389			if (dval(&d) > dval(&eps))
390				goto one_digit;
391			if (dval(&d) < -dval(&eps))
392				goto no_digits;
393			goto fast_failed;
394			}
395#ifndef No_leftright
396		if (leftright) {
397			/* Use Steele & White method of only
398			 * generating digits needed.
399			 */
400			dval(&eps) = ds*0.5/tens[ilim-1] - dval(&eps);
401			for(i = 0;;) {
402				L = (Long)(dval(&d)/ds);
403				dval(&d) -= L*ds;
404				*s++ = '0' + (int)L;
405				if (dval(&d) < dval(&eps)) {
406					if (dval(&d))
407						inex = STRTOG_Inexlo;
408					goto ret1;
409					}
410				if (ds - dval(&d) < dval(&eps))
411					goto bump_up;
412				if (++i >= ilim)
413					break;
414				dval(&eps) *= 10.;
415				dval(&d) *= 10.;
416				}
417			}
418		else {
419#endif
420			/* Generate ilim digits, then fix them up. */
421			dval(&eps) *= tens[ilim-1];
422			for(i = 1;; i++, dval(&d) *= 10.) {
423				if ( (L = (Long)(dval(&d)/ds)) !=0)
424					dval(&d) -= L*ds;
425				*s++ = '0' + (int)L;
426				if (i == ilim) {
427					ds *= 0.5;
428					if (dval(&d) > ds + dval(&eps))
429						goto bump_up;
430					else if (dval(&d) < ds - dval(&eps)) {
431						if (dval(&d))
432							inex = STRTOG_Inexlo;
433						goto clear_trailing0;
434						}
435					break;
436					}
437				}
438#ifndef No_leftright
439			}
440#endif
441 fast_failed:
442		s = s0;
443		dval(&d) = d2;
444		k = k0;
445		ilim = ilim0;
446		}
447
448	/* Do we have a "small" integer? */
449
450	if (be >= 0 && k <= Int_max) {
451		/* Yes. */
452		ds = tens[k];
453		if (ndigits < 0 && ilim <= 0) {
454			S = mhi = 0;
455			if (ilim < 0 || dval(&d) <= 5*ds)
456				goto no_digits;
457			goto one_digit;
458			}
459		for(i = 1;; i++, dval(&d) *= 10.) {
460			L = dval(&d) / ds;
461			dval(&d) -= L*ds;
462#ifdef Check_FLT_ROUNDS
463			/* If FLT_ROUNDS == 2, L will usually be high by 1 */
464			if (dval(&d) < 0) {
465				L--;
466				dval(&d) += ds;
467				}
468#endif
469			*s++ = '0' + (int)L;
470			if (dval(&d) == 0.)
471				break;
472			if (i == ilim) {
473				if (rdir) {
474					if (rdir == 1)
475						goto bump_up;
476					inex = STRTOG_Inexlo;
477					goto ret1;
478					}
479				dval(&d) += dval(&d);
480#ifdef ROUND_BIASED
481				if (dval(&d) >= ds)
482#else
483				if (dval(&d) > ds || (dval(&d) == ds && L & 1))
484#endif
485					{
486 bump_up:
487					inex = STRTOG_Inexhi;
488					while(*--s == '9')
489						if (s == s0) {
490							k++;
491							*s = '0';
492							break;
493							}
494					++*s++;
495					}
496				else {
497					inex = STRTOG_Inexlo;
498 clear_trailing0:
499					while(*--s == '0'){}
500					++s;
501					}
502				break;
503				}
504			}
505		goto ret1;
506		}
507
508	m2 = b2;
509	m5 = b5;
510	mhi = mlo = 0;
511	if (leftright) {
512		i = nbits - bbits;
513		if (be - i++ < fpi->emin && mode != 3 && mode != 5) {
514			/* denormal */
515			i = be - fpi->emin + 1;
516			if (mode >= 2 && ilim > 0 && ilim < i)
517				goto small_ilim;
518			}
519		else if (mode >= 2) {
520 small_ilim:
521			j = ilim - 1;
522			if (m5 >= j)
523				m5 -= j;
524			else {
525				s5 += j -= m5;
526				b5 += j;
527				m5 = 0;
528				}
529			if ((i = ilim) < 0) {
530				m2 -= i;
531				i = 0;
532				}
533			}
534		b2 += i;
535		s2 += i;
536		mhi = i2b(1);
537		if (mhi == NULL)
538			return (NULL);
539		}
540	if (m2 > 0 && s2 > 0) {
541		i = m2 < s2 ? m2 : s2;
542		b2 -= i;
543		m2 -= i;
544		s2 -= i;
545		}
546	if (b5 > 0) {
547		if (leftright) {
548			if (m5 > 0) {
549				mhi = pow5mult(mhi, m5);
550				if (mhi == NULL)
551					return (NULL);
552				b1 = mult(mhi, b);
553				if (b1 == NULL)
554					return (NULL);
555				Bfree(b);
556				b = b1;
557				}
558			if ( (j = b5 - m5) !=0) {
559				b = pow5mult(b, j);
560				if (b == NULL)
561					return (NULL);
562				}
563			}
564		else {
565			b = pow5mult(b, b5);
566			if (b == NULL)
567				return (NULL);
568			}
569		}
570	S = i2b(1);
571	if (S == NULL)
572		return (NULL);
573	if (s5 > 0) {
574		S = pow5mult(S, s5);
575		if (S == NULL)
576			return (NULL);
577		}
578
579	/* Check for special case that d is a normalized power of 2. */
580
581	spec_case = 0;
582	if (mode < 2) {
583		if (bbits == 1 && be0 > fpi->emin + 1) {
584			/* The special case */
585			b2++;
586			s2++;
587			spec_case = 1;
588			}
589		}
590
591	/* Arrange for convenient computation of quotients:
592	 * shift left if necessary so divisor has 4 leading 0 bits.
593	 *
594	 * Perhaps we should just compute leading 28 bits of S once
595	 * and for all and pass them and a shift to quorem, so it
596	 * can do shifts and ors to compute the numerator for q.
597	 */
598	i = ((s5 ? hi0bits(S->x[S->wds-1]) : ULbits - 1) - s2 - 4) & kmask;
599	m2 += i;
600	if ((b2 += i) > 0) {
601		b = lshift(b, b2);
602		if (b == NULL)
603			return (NULL);
604		}
605	if ((s2 += i) > 0) {
606		S = lshift(S, s2);
607		if (S == NULL)
608			return (NULL);
609		}
610	if (k_check) {
611		if (cmp(b,S) < 0) {
612			k--;
613			b = multadd(b, 10, 0);	/* we botched the k estimate */
614			if (b == NULL)
615				return (NULL);
616			if (leftright) {
617				mhi = multadd(mhi, 10, 0);
618				if (mhi == NULL)
619					return (NULL);
620				}
621			ilim = ilim1;
622			}
623		}
624	if (ilim <= 0 && mode > 2) {
625		S = multadd(S,5,0);
626		if (S == NULL)
627			return (NULL);
628		if (ilim < 0 || cmp(b,S) <= 0) {
629			/* no digits, fcvt style */
630 no_digits:
631			k = -1 - ndigits;
632			inex = STRTOG_Inexlo;
633			goto ret;
634			}
635 one_digit:
636		inex = STRTOG_Inexhi;
637		*s++ = '1';
638		k++;
639		goto ret;
640		}
641	if (leftright) {
642		if (m2 > 0) {
643			mhi = lshift(mhi, m2);
644			if (mhi == NULL)
645				return (NULL);
646			}
647
648		/* Compute mlo -- check for special case
649		 * that d is a normalized power of 2.
650		 */
651
652		mlo = mhi;
653		if (spec_case) {
654			mhi = Balloc(mhi->k);
655			if (mhi == NULL)
656				return (NULL);
657			Bcopy(mhi, mlo);
658			mhi = lshift(mhi, 1);
659			if (mhi == NULL)
660				return (NULL);
661			}
662
663		for(i = 1;;i++) {
664			dig = quorem(b,S) + '0';
665			/* Do we yet have the shortest decimal string
666			 * that will round to d?
667			 */
668			j = cmp(b, mlo);
669			delta = diff(S, mhi);
670			if (delta == NULL)
671				return (NULL);
672			j1 = delta->sign ? 1 : cmp(b, delta);
673			Bfree(delta);
674#ifndef ROUND_BIASED
675			if (j1 == 0 && !mode && !(bits[0] & 1) && !rdir) {
676				if (dig == '9')
677					goto round_9_up;
678				if (j <= 0) {
679					if (b->wds > 1 || b->x[0])
680						inex = STRTOG_Inexlo;
681					}
682				else {
683					dig++;
684					inex = STRTOG_Inexhi;
685					}
686				*s++ = dig;
687				goto ret;
688				}
689#endif
690			if (j < 0 || (j == 0 && !mode
691#ifndef ROUND_BIASED
692							&& !(bits[0] & 1)
693#endif
694					)) {
695				if (rdir && (b->wds > 1 || b->x[0])) {
696					if (rdir == 2) {
697						inex = STRTOG_Inexlo;
698						goto accept;
699						}
700					while (cmp(S,mhi) > 0) {
701						*s++ = dig;
702						mhi1 = multadd(mhi, 10, 0);
703						if (mhi1 == NULL)
704							return (NULL);
705						if (mlo == mhi)
706							mlo = mhi1;
707						mhi = mhi1;
708						b = multadd(b, 10, 0);
709						if (b == NULL)
710							return (NULL);
711						dig = quorem(b,S) + '0';
712						}
713					if (dig++ == '9')
714						goto round_9_up;
715					inex = STRTOG_Inexhi;
716					goto accept;
717					}
718				if (j1 > 0) {
719					b = lshift(b, 1);
720					if (b == NULL)
721						return (NULL);
722					j1 = cmp(b, S);
723#ifdef ROUND_BIASED
724					if (j1 >= 0 /*)*/
725#else
726					if ((j1 > 0 || (j1 == 0 && dig & 1))
727#endif
728					&& dig++ == '9')
729						goto round_9_up;
730					inex = STRTOG_Inexhi;
731					}
732				if (b->wds > 1 || b->x[0])
733					inex = STRTOG_Inexlo;
734 accept:
735				*s++ = dig;
736				goto ret;
737				}
738			if (j1 > 0 && rdir != 2) {
739				if (dig == '9') { /* possible if i == 1 */
740 round_9_up:
741					*s++ = '9';
742					inex = STRTOG_Inexhi;
743					goto roundoff;
744					}
745				inex = STRTOG_Inexhi;
746				*s++ = dig + 1;
747				goto ret;
748				}
749			*s++ = dig;
750			if (i == ilim)
751				break;
752			b = multadd(b, 10, 0);
753			if (b == NULL)
754				return (NULL);
755			if (mlo == mhi) {
756				mlo = mhi = multadd(mhi, 10, 0);
757				if (mlo == NULL)
758					return (NULL);
759				}
760			else {
761				mlo = multadd(mlo, 10, 0);
762				if (mlo == NULL)
763					return (NULL);
764				mhi = multadd(mhi, 10, 0);
765				if (mhi == NULL)
766					return (NULL);
767				}
768			}
769		}
770	else
771		for(i = 1;; i++) {
772			*s++ = dig = quorem(b,S) + '0';
773			if (i >= ilim)
774				break;
775			b = multadd(b, 10, 0);
776			if (b == NULL)
777				return (NULL);
778			}
779
780	/* Round off last digit */
781
782	if (rdir) {
783		if (rdir == 2 || (b->wds <= 1 && !b->x[0]))
784			goto chopzeros;
785		goto roundoff;
786		}
787	b = lshift(b, 1);
788	if (b == NULL)
789		return (NULL);
790	j = cmp(b, S);
791#ifdef ROUND_BIASED
792	if (j >= 0)
793#else
794	if (j > 0 || (j == 0 && dig & 1))
795#endif
796		{
797 roundoff:
798		inex = STRTOG_Inexhi;
799		while(*--s == '9')
800			if (s == s0) {
801				k++;
802				*s++ = '1';
803				goto ret;
804				}
805		++*s++;
806		}
807	else {
808 chopzeros:
809		if (b->wds > 1 || b->x[0])
810			inex = STRTOG_Inexlo;
811		while(*--s == '0'){}
812		++s;
813		}
814 ret:
815	Bfree(S);
816	if (mhi) {
817		if (mlo && mlo != mhi)
818			Bfree(mlo);
819		Bfree(mhi);
820		}
821 ret1:
822	Bfree(b);
823	*s = 0;
824	*decpt = k + 1;
825	if (rve)
826		*rve = s;
827	*kindp |= inex;
828	return s0;
829	}
830