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