GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: lib/libc/gdtoa/dtoa.c Lines: 127 310 41.0 %
Date: 2017-11-07 Branches: 78 277 28.2 %

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
/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
35
 *
36
 * Inspired by "How to Print Floating-Point Numbers Accurately" by
37
 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
38
 *
39
 * Modifications:
40
 *	1. Rather than iterating, we use a simple numeric overestimate
41
 *	   to determine k = floor(log10(d)).  We scale relevant
42
 *	   quantities using O(log2(k)) rather than O(k) multiplications.
43
 *	2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
44
 *	   try to generate digits strictly left to right.  Instead, we
45
 *	   compute with fewer bits and propagate the carry if necessary
46
 *	   when rounding the final digit up.  This is often faster.
47
 *	3. Under the assumption that input will be rounded nearest,
48
 *	   mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
49
 *	   That is, we allow equality in stopping tests when the
50
 *	   round-nearest rule will give the same floating-point value
51
 *	   as would satisfaction of the stopping test with strict
52
 *	   inequality.
53
 *	4. We remove common factors of powers of 2 from relevant
54
 *	   quantities.
55
 *	5. When converting floating-point integers less than 1e16,
56
 *	   we use floating-point arithmetic rather than resorting
57
 *	   to multiple-precision integers.
58
 *	6. When asked to produce fewer than 15 digits, we first try
59
 *	   to get by with floating-point arithmetic; we resort to
60
 *	   multiple-precision integer arithmetic only if we cannot
61
 *	   guarantee that the floating-point calculation has given
62
 *	   the correctly rounded result.  For k requested digits and
63
 *	   "uniformly" distributed input, the probability is
64
 *	   something like 10^(k-15) that we must resort to the Long
65
 *	   calculation.
66
 */
67
68
#ifdef Honor_FLT_ROUNDS
69
#undef Check_FLT_ROUNDS
70
#define Check_FLT_ROUNDS
71
#else
72
#define Rounding Flt_Rounds
73
#endif
74
75
 char *
76
dtoa
77
#ifdef KR_headers
78
	(d0, mode, ndigits, decpt, sign, rve)
79
	double d0; int mode, ndigits, *decpt, *sign; char **rve;
80
#else
81
	(double d0, int mode, int ndigits, int *decpt, int *sign, char **rve)
82
#endif
83
{
84
 /*	Arguments ndigits, decpt, sign are similar to those
85
	of ecvt and fcvt; trailing zeros are suppressed from
86
	the returned string.  If not null, *rve is set to point
87
	to the end of the return value.  If d is +-Infinity or NaN,
88
	then *decpt is set to 9999.
89
90
	mode:
91
		0 ==> shortest string that yields d when read in
92
			and rounded to nearest.
93
		1 ==> like 0, but with Steele & White stopping rule;
94
			e.g. with IEEE P754 arithmetic , mode 0 gives
95
			1e23 whereas mode 1 gives 9.999999999999999e22.
96
		2 ==> max(1,ndigits) significant digits.  This gives a
97
			return value similar to that of ecvt, except
98
			that trailing zeros are suppressed.
99
		3 ==> through ndigits past the decimal point.  This
100
			gives a return value similar to that from fcvt,
101
			except that trailing zeros are suppressed, and
102
			ndigits can be negative.
103
		4,5 ==> similar to 2 and 3, respectively, but (in
104
			round-nearest mode) with the tests of mode 0 to
105
			possibly return a shorter string that rounds to d.
106
			With IEEE arithmetic and compilation with
107
			-DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
108
			as modes 2 and 3 when FLT_ROUNDS != 1.
109
		6-9 ==> Debugging modes similar to mode - 4:  don't try
110
			fast floating-point estimate (if applicable).
111
112
		Values of mode other than 0-9 are treated as mode 0.
113
114
		Sufficient space is allocated to the return value
115
		to hold the suppressed trailing zeros.
116
	*/
117
118
844
	int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
119
		j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
120
		spec_case, try_quick;
121
	Long L;
122
#ifndef Sudden_Underflow
123
	int denorm;
124
	ULong x;
125
#endif
126
	Bigint *b, *b1, *delta, *mlo, *mhi, *S;
127
	U d, d2, eps;
128
	double ds;
129
	char *s, *s0;
130
#ifdef SET_INEXACT
131
	int inexact, oldinexact;
132
#endif
133
#ifdef Honor_FLT_ROUNDS /*{*/
134
	int Rounding;
135
#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
136
	Rounding = Flt_Rounds;
137
#else /*}{*/
138
	Rounding = 1;
139
	switch(fegetround()) {
140
	  case FE_TOWARDZERO:	Rounding = 0; break;
141
	  case FE_UPWARD:	Rounding = 2; break;
142
	  case FE_DOWNWARD:	Rounding = 3;
143
	  }
144
#endif /*}}*/
145
#endif /*}*/
146
147
#ifndef MULTIPLE_THREADS
148
	if (dtoa_result) {
149
		freedtoa(dtoa_result);
150
		dtoa_result = 0;
151
		}
152
#endif
153
	d.d = d0;
154
422
	if (word0(&d) & Sign_bit) {
155
		/* set sign for everything, including 0's and NaNs */
156
		*sign = 1;
157
		word0(&d) &= ~Sign_bit;	/* clear sign bit */
158
		}
159
	else
160
422
		*sign = 0;
161
162
#if defined(IEEE_Arith) + defined(VAX)
163
#ifdef IEEE_Arith
164
422
	if ((word0(&d) & Exp_mask) == Exp_mask)
165
#else
166
	if (word0(&d)  == 0x8000)
167
#endif
168
		{
169
		/* Infinity or NaN */
170
		*decpt = 9999;
171
#ifdef IEEE_Arith
172
		if (!word1(&d) && !(word0(&d) & 0xfffff))
173
			return nrv_alloc("Infinity", rve, 8);
174
#endif
175
		return nrv_alloc("NaN", rve, 3);
176
		}
177
#endif
178
#ifdef IBM
179
	dval(&d) += 0; /* normalize */
180
#endif
181
422
	if (!dval(&d)) {
182
200
		*decpt = 1;
183
200
		return nrv_alloc("0", rve, 1);
184
		}
185
186
#ifdef SET_INEXACT
187
	try_quick = oldinexact = get_inexact();
188
	inexact = 1;
189
#endif
190
#ifdef Honor_FLT_ROUNDS
191
	if (Rounding >= 2) {
192
		if (*sign)
193
			Rounding = Rounding == 2 ? 0 : 2;
194
		else
195
			if (Rounding != 2)
196
				Rounding = 0;
197
		}
198
#endif
199
200
222
	b = d2b(dval(&d), &be, &bbits);
201
222
	if (b == NULL)
202
		return (NULL);
203
#ifdef Sudden_Underflow
204
	i = (int)(word0(&d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
205
#else
206
222
	if (( i = (int)(word0(&d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)) )!=0) {
207
#endif
208
		dval(&d2) = dval(&d);
209
222
		word0(&d2) &= Frac_mask1;
210
222
		word0(&d2) |= Exp_11;
211
#ifdef IBM
212
		if (( j = 11 - hi0bits(word0(&d2) & Frac_mask) )!=0)
213
			dval(&d2) /= 1 << j;
214
#endif
215
216
		/* log(x)	~=~ log(1.5) + (x-1.5)/1.5
217
		 * log10(x)	 =  log(x) / log(10)
218
		 *		~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
219
		 * log10(&d) = (i-Bias)*log(2)/log(10) + log10(&d2)
220
		 *
221
		 * This suggests computing an approximation k to log10(&d) by
222
		 *
223
		 * k = (i - Bias)*0.301029995663981
224
		 *	+ ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
225
		 *
226
		 * We want k to be too large rather than too small.
227
		 * The error in the first-order Taylor series approximation
228
		 * is in our favor, so we just round up the constant enough
229
		 * to compensate for any error in the multiplication of
230
		 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
231
		 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
232
		 * adding 1e-13 to the constant term more than suffices.
233
		 * Hence we adjust the constant term to 0.1760912590558.
234
		 * (We could get a more accurate k by invoking log10,
235
		 *  but this is probably not worthwhile.)
236
		 */
237
238
222
		i -= Bias;
239
#ifdef IBM
240
		i <<= 2;
241
		i += j;
242
#endif
243
#ifndef Sudden_Underflow
244
		denorm = 0;
245
222
		}
246
	else {
247
		/* d is denormalized */
248
249
		i = bbits + be + (Bias + (P-1) - 1);
250
		x = i > 32  ? word0(&d) << (64 - i) | word1(&d) >> (i - 32)
251
			    : word1(&d) << (32 - i);
252
		dval(&d2) = x;
253
		word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
254
		i -= (Bias + (P-1) - 1) + 1;
255
		denorm = 1;
256
		}
257
#endif
258
222
	ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
259
222
	k = (int)ds;
260

412
	if (ds < 0. && ds != k)
261
190
		k--;	/* want k = floor(ds) */
262
	k_check = 1;
263
222
	if (k >= 0 && k <= Ten_pmax) {
264
32
		if (dval(&d) < tens[k])
265
			k--;
266
		k_check = 0;
267
32
		}
268
222
	j = bbits - i - 1;
269
222
	if (j >= 0) {
270
		b2 = 0;
271
		s2 = j;
272
222
		}
273
	else {
274
		b2 = -j;
275
		s2 = 0;
276
		}
277
222
	if (k >= 0) {
278
		b5 = 0;
279
		s5 = k;
280
32
		s2 += k;
281
32
		}
282
	else {
283
190
		b2 -= k;
284
190
		b5 = -k;
285
		s5 = 0;
286
		}
287
222
	if (mode < 0 || mode > 9)
288
		mode = 0;
289
290
#ifndef SET_INEXACT
291
#ifdef Check_FLT_ROUNDS
292
	try_quick = Rounding == 1;
293
#else
294
	try_quick = 1;
295
#endif
296
#endif /*SET_INEXACT*/
297
298
222
	if (mode > 5) {
299
		mode -= 4;
300
		try_quick = 0;
301
		}
302
	leftright = 1;
303
	ilim = ilim1 = -1;	/* Values for cases 0 and 1; done here to */
304
				/* silence erroneous "gcc -Wall" warning. */
305

666
	switch(mode) {
306
		case 0:
307
		case 1:
308
			i = 18;
309
			ndigits = 0;
310
			break;
311
		case 2:
312
			leftright = 0;
313
			/* no break */
314
		case 4:
315
			if (ndigits <= 0)
316
				ndigits = 1;
317
			ilim = ilim1 = i = ndigits;
318
			break;
319
		case 3:
320
222
			leftright = 0;
321
			/* no break */
322
		case 5:
323
222
			i = ndigits + k + 1;
324
			ilim = i;
325
			ilim1 = i - 1;
326
222
			if (i <= 0)
327
38
				i = 1;
328
		}
329
222
	s = s0 = rv_alloc(i);
330
222
	if (s == NULL)
331
		return (NULL);
332
333
#ifdef Honor_FLT_ROUNDS
334
	if (mode > 1 && Rounding != 1)
335
		leftright = 0;
336
#endif
337
338
222
	if (ilim >= 0 && ilim <= Quick_max && try_quick) {
339
340
		/* Try to get by with floating-point arithmetic. */
341
342
		i = 0;
343
		dval(&d2) = dval(&d);
344
		k0 = k;
345
		ilim0 = ilim;
346
		ieps = 2; /* conservative */
347
204
		if (k > 0) {
348
12
			ds = tens[k&0xf];
349
12
			j = k >> 4;
350
12
			if (j & Bletch) {
351
				/* prevent overflows */
352
				j &= Bletch - 1;
353
				dval(&d) /= bigtens[n_bigtens-1];
354
				ieps++;
355
				}
356
12
			for(; j; j >>= 1, i++)
357
				if (j & 1) {
358
					ieps++;
359
					ds *= bigtens[i];
360
					}
361
12
			dval(&d) /= ds;
362
12
			}
363
192
		else if (( j1 = -k )!=0) {
364
172
			dval(&d) *= tens[j1 & 0xf];
365
344
			for(j = j1 >> 4; j; j >>= 1, i++)
366
				if (j & 1) {
367
					ieps++;
368
					dval(&d) *= bigtens[i];
369
					}
370
			}
371

376
		if (k_check && dval(&d) < 1. && ilim > 0) {
372
			if (ilim1 <= 0)
373
				goto fast_failed;
374
			ilim = ilim1;
375
			k--;
376
			dval(&d) *= 10.;
377
			ieps++;
378
			}
379
204
		dval(&eps) = ieps*dval(&d) + 7.;
380
204
		word0(&eps) -= (P-1)*Exp_msk1;
381
204
		if (ilim == 0) {
382
			S = mhi = 0;
383
20
			dval(&d) -= 5.;
384
20
			if (dval(&d) > dval(&eps))
385
				goto one_digit;
386
			if (dval(&d) < -dval(&eps))
387
				goto no_digits;
388
			goto fast_failed;
389
			}
390
#ifndef No_leftright
391
184
		if (leftright) {
392
			/* Use Steele & White method of only
393
			 * generating digits needed.
394
			 */
395
			dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
396
			for(i = 0;;) {
397
				L = dval(&d);
398
				dval(&d) -= L;
399
				*s++ = '0' + (int)L;
400
				if (dval(&d) < dval(&eps))
401
					goto ret1;
402
				if (1. - dval(&d) < dval(&eps))
403
					goto bump_up;
404
				if (++i >= ilim)
405
					break;
406
				dval(&eps) *= 10.;
407
				dval(&d) *= 10.;
408
				}
409
			}
410
		else {
411
#endif
412
			/* Generate ilim digits, then fix them up. */
413
184
			dval(&eps) *= tens[ilim-1];
414
262
			for(i = 1;; i++, dval(&d) *= 10.) {
415
262
				L = (Long)(dval(&d));
416
262
				if (!(dval(&d) -= L))
417
					ilim = i;
418
262
				*s++ = '0' + (int)L;
419
262
				if (i == ilim) {
420
184
					if (dval(&d) > 0.5 + dval(&eps))
421
						goto bump_up;
422
147
					else if (dval(&d) < 0.5 - dval(&eps)) {
423
153
						while(*--s == '0');
424
147
						s++;
425
147
						goto ret1;
426
						}
427
					break;
428
					}
429
				}
430
#ifndef No_leftright
431
			}
432
#endif
433
 fast_failed:
434
		s = s0;
435
		dval(&d) = dval(&d2);
436
		k = k0;
437
		ilim = ilim0;
438
		}
439
440
	/* Do we have a "small" integer? */
441
442
18
	if (be >= 0 && k <= Int_max) {
443
		/* Yes. */
444
		ds = tens[k];
445
		if (ndigits < 0 && ilim <= 0) {
446
			S = mhi = 0;
447
			if (ilim < 0 || dval(&d) <= 5*ds)
448
				goto no_digits;
449
			goto one_digit;
450
			}
451
		for(i = 1;; i++, dval(&d) *= 10.) {
452
			L = (Long)(dval(&d) / ds);
453
			dval(&d) -= L*ds;
454
#ifdef Check_FLT_ROUNDS
455
			/* If FLT_ROUNDS == 2, L will usually be high by 1 */
456
			if (dval(&d) < 0) {
457
				L--;
458
				dval(&d) += ds;
459
				}
460
#endif
461
			*s++ = '0' + (int)L;
462
			if (!dval(&d)) {
463
#ifdef SET_INEXACT
464
				inexact = 0;
465
#endif
466
				break;
467
				}
468
			if (i == ilim) {
469
#ifdef Honor_FLT_ROUNDS
470
				if (mode > 1)
471
				switch(Rounding) {
472
				  case 0: goto ret1;
473
				  case 2: goto bump_up;
474
				  }
475
#endif
476
				dval(&d) += dval(&d);
477
#ifdef ROUND_BIASED
478
				if (dval(&d) >= ds)
479
#else
480
				if (dval(&d) > ds || (dval(&d) == ds && L & 1))
481
#endif
482
					{
483
 bump_up:
484
38
					while(*--s == '9')
485
1
						if (s == s0) {
486
							k++;
487
							*s = '0';
488
							break;
489
							}
490
37
					++*s++;
491
37
					}
492
				break;
493
				}
494
			}
495
		goto ret1;
496
		}
497
498
	m2 = b2;
499
	m5 = b5;
500
	mhi = mlo = 0;
501
18
	if (leftright) {
502
		i =
503
#ifndef Sudden_Underflow
504
			denorm ? be + (Bias + (P-1) - 1 + 1) :
505
#endif
506
#ifdef IBM
507
			1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
508
#else
509
			1 + P - bbits;
510
#endif
511
		b2 += i;
512
		s2 += i;
513
		mhi = i2b(1);
514
		if (mhi == NULL)
515
			return (NULL);
516
		}
517
18
	if (m2 > 0 && s2 > 0) {
518
18
		i = m2 < s2 ? m2 : s2;
519
18
		b2 -= i;
520
18
		m2 -= i;
521
18
		s2 -= i;
522
18
		}
523
18
	if (b5 > 0) {
524
18
		if (leftright) {
525
			if (m5 > 0) {
526
				mhi = pow5mult(mhi, m5);
527
				if (mhi == NULL)
528
					return (NULL);
529
				b1 = mult(mhi, b);
530
				if (b1 == NULL)
531
					return (NULL);
532
				Bfree(b);
533
				b = b1;
534
				}
535
			if (( j = b5 - m5 )!=0) {
536
				b = pow5mult(b, j);
537
				if (b == NULL)
538
					return (NULL);
539
				}
540
			}
541
		else {
542
18
			b = pow5mult(b, b5);
543
18
			if (b == NULL)
544
				return (NULL);
545
			}
546
		}
547
18
	S = i2b(1);
548
18
	if (S == NULL)
549
		return (NULL);
550
18
	if (s5 > 0) {
551
		S = pow5mult(S, s5);
552
		if (S == NULL)
553
			return (NULL);
554
		}
555
556
	/* Check for special case that d is a normalized power of 2. */
557
558
	spec_case = 0;
559
18
	if ((mode < 2 || leftright)
560
#ifdef Honor_FLT_ROUNDS
561
			&& Rounding == 1
562
#endif
563
				) {
564
		if (!word1(&d) && !(word0(&d) & Bndry_mask)
565
#ifndef Sudden_Underflow
566
		 && word0(&d) & (Exp_mask & ~Exp_msk1)
567
#endif
568
				) {
569
			/* The special case */
570
			b2 += Log2P;
571
			s2 += Log2P;
572
			spec_case = 1;
573
			}
574
		}
575
576
	/* Arrange for convenient computation of quotients:
577
	 * shift left if necessary so divisor has 4 leading 0 bits.
578
	 *
579
	 * Perhaps we should just compute leading 28 bits of S once
580
	 * and for all and pass them and a shift to quorem, so it
581
	 * can do shifts and ors to compute the numerator for q.
582
	 */
583
#ifdef Pack_32
584

36
	if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f )!=0)
585
18
		i = 32 - i;
586
#else
587
	if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf )!=0)
588
		i = 16 - i;
589
#endif
590
18
	if (i > 4) {
591
15
		i -= 4;
592
15
		b2 += i;
593
15
		m2 += i;
594
15
		s2 += i;
595
15
		}
596
3
	else if (i < 4) {
597
1
		i += 28;
598
1
		b2 += i;
599
1
		m2 += i;
600
1
		s2 += i;
601
1
		}
602
18
	if (b2 > 0) {
603
16
		b = lshift(b, b2);
604
16
		if (b == NULL)
605
			return (NULL);
606
		}
607
18
	if (s2 > 0) {
608
18
		S = lshift(S, s2);
609
18
		if (S == NULL)
610
			return (NULL);
611
		}
612
18
	if (k_check) {
613
18
		if (cmp(b,S) < 0) {
614
			k--;
615
			b = multadd(b, 10, 0);	/* we botched the k estimate */
616
			if (b == NULL)
617
				return (NULL);
618
			if (leftright) {
619
				mhi = multadd(mhi, 10, 0);
620
				if (mhi == NULL)
621
					return (NULL);
622
				}
623
			ilim = ilim1;
624
			}
625
		}
626

36
	if (ilim <= 0 && (mode == 3 || mode == 5)) {
627
18
		S = multadd(S,5,0);
628
18
		if (S == NULL)
629
			return (NULL);
630

18
		if (ilim < 0 || cmp(b,S) <= 0) {
631
			/* no digits, fcvt style */
632
 no_digits:
633
18
			k = -1 - ndigits;
634
18
			goto ret;
635
			}
636
 one_digit:
637
20
		*s++ = '1';
638
20
		k++;
639
20
		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, Log2P);
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 != 1 && !(word1(&d) & 1)
676
#ifdef Honor_FLT_ROUNDS
677
				&& Rounding >= 1
678
#endif
679
								   ) {
680
				if (dig == '9')
681
					goto round_9_up;
682
				if (j > 0)
683
					dig++;
684
#ifdef SET_INEXACT
685
				else if (!b->x[0] && b->wds <= 1)
686
					inexact = 0;
687
#endif
688
				*s++ = dig;
689
				goto ret;
690
				}
691
#endif
692
			if (j < 0 || (j == 0 && mode != 1
693
#ifndef ROUND_BIASED
694
							&& !(word1(&d) & 1)
695
#endif
696
					)) {
697
				if (!b->x[0] && b->wds <= 1) {
698
#ifdef SET_INEXACT
699
					inexact = 0;
700
#endif
701
					goto accept_dig;
702
					}
703
#ifdef Honor_FLT_ROUNDS
704
				if (mode > 1)
705
				 switch(Rounding) {
706
				  case 0: goto accept_dig;
707
				  case 2: goto keep_dig;
708
				  }
709
#endif /*Honor_FLT_ROUNDS*/
710
				if (j1 > 0) {
711
					b = lshift(b, 1);
712
					if (b == NULL)
713
						return (NULL);
714
					j1 = cmp(b, S);
715
#ifdef ROUND_BIASED
716
					if (j1 >= 0 /*)*/
717
#else
718
					if ((j1 > 0 || (j1 == 0 && dig & 1))
719
#endif
720
					&& dig++ == '9')
721
						goto round_9_up;
722
					}
723
 accept_dig:
724
				*s++ = dig;
725
				goto ret;
726
				}
727
			if (j1 > 0) {
728
#ifdef Honor_FLT_ROUNDS
729
				if (!Rounding)
730
					goto accept_dig;
731
#endif
732
				if (dig == '9') { /* possible if i == 1 */
733
 round_9_up:
734
					*s++ = '9';
735
					goto roundoff;
736
					}
737
				*s++ = dig + 1;
738
				goto ret;
739
				}
740
#ifdef Honor_FLT_ROUNDS
741
 keep_dig:
742
#endif
743
			*s++ = dig;
744
			if (i == ilim)
745
				break;
746
			b = multadd(b, 10, 0);
747
			if (b == NULL)
748
				return (NULL);
749
			if (mlo == mhi) {
750
				mlo = mhi = multadd(mhi, 10, 0);
751
				if (mlo == NULL)
752
					return (NULL);
753
				}
754
			else {
755
				mlo = multadd(mlo, 10, 0);
756
				if (mlo == NULL)
757
					return (NULL);
758
				mhi = multadd(mhi, 10, 0);
759
				if (mhi == NULL)
760
					return (NULL);
761
				}
762
			}
763
		}
764
	else
765
		for(i = 1;; i++) {
766
			*s++ = dig = quorem(b,S) + '0';
767
			if (!b->x[0] && b->wds <= 1) {
768
#ifdef SET_INEXACT
769
				inexact = 0;
770
#endif
771
				goto ret;
772
				}
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
#ifdef Honor_FLT_ROUNDS
783
	switch(Rounding) {
784
	  case 0: goto trimzeros;
785
	  case 2: goto roundoff;
786
	  }
787
#endif
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
		while(*--s == '9')
800
			if (s == s0) {
801
				k++;
802
				*s++ = '1';
803
				goto ret;
804
				}
805
		++*s++;
806
		}
807
	else {
808
#ifdef Honor_FLT_ROUNDS
809
 trimzeros:
810
#endif
811
		while(*--s == '0');
812
		s++;
813
		}
814
 ret:
815
38
	Bfree(S);
816
38
	if (mhi) {
817
		if (mlo && mlo != mhi)
818
			Bfree(mlo);
819
		Bfree(mhi);
820
		}
821
 ret1:
822
#ifdef SET_INEXACT
823
	if (inexact) {
824
		if (!oldinexact) {
825
			word0(&d) = Exp_1 + (70 << Exp_shift);
826
			word1(&d) = 0;
827
			dval(&d) += 1.;
828
			}
829
		}
830
	else if (!oldinexact)
831
		clear_inexact();
832
#endif
833
222
	Bfree(b);
834
222
	*s = 0;
835
222
	*decpt = k + 1;
836
222
	if (rve)
837
222
		*rve = s;
838
222
	return s0;
839
422
	}
840
DEF_STRONG(dtoa);