GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: lib/libc/gdtoa/gdtoa.c Lines: 0 361 0.0 %
Date: 2017-11-07 Branches: 0 333 0.0 %

Line Branch Exec Source
1
/****************************************************************
2
3
The author of this software is David M. Gay.
4
5
Copyright (C) 1998, 1999 by Lucent Technologies
6
All Rights Reserved
7
8
Permission to use, copy, modify, and distribute this software and
9
its documentation for any purpose and without fee is hereby
10
granted, provided that the above copyright notice appear in all
11
copies and that both that the copyright notice and this
12
permission notice and warranty disclaimer appear in supporting
13
documentation, and that the name of Lucent or any of its entities
14
not be used in advertising or publicity pertaining to
15
distribution of the software without specific, written prior
16
permission.
17
18
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
19
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
20
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
21
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
22
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
23
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
24
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
25
THIS 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
36
bitstob(bits, nbits, bbits) ULong *bits; int nbits; int *bbits;
37
#else
38
bitstob(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 *
114
gdtoa
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
DEF_STRONG(gdtoa);