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