?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Cygwin: mmap: allow remapping part of an existing anonymous mapping
[newlib-cygwin.git] / newlib / libc / stdlib / gdtoa-gdtoa.c
blobda2338c48a4c190e07326deeee14ede43a5527ae
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 <newlib.h>
33 #include <sys/config.h>
34
35 #ifdef _USE_GDTOA
36 #include "gdtoaimp.h"
37
38 static Bigint *
39 #ifdef KR_headers
40 bitstob(ptr, bits, nbits, bbits) ULong *bits;
41 struct _reent ptr, int nbits; int *bbits;
42 #else
43 bitstob(struct _reent *ptr, ULong *bits, int nbits, int *bbits)
44 #endif
45 {
46 int i, k;
47 Bigint *b;
48 ULong *be, *x, *x0;
49
50 i = ULbits;
51 k = 0;
52 while(i < nbits) {
53 i <<= 1;
54 k++;
55 }
56 #ifndef Pack_32
57 if (!k)
58 k = 1;
59 #endif
60 b = Balloc(ptr, k);
61 if (b == NULL)
62 return (NULL);
63 be = bits + ((nbits - 1) >> kshift);
64 x = x0 = b->_x;
65 do {
66 *x++ = *bits & ALL_ON;
67 #ifdef Pack_16
68 *x++ = (*bits >> 16) & ALL_ON;
69 #endif
70 } while(++bits <= be);
71 i = x - x0;
72 while(!x0[--i])
73 if (!i) {
74 b->_wds = 0;
75 *bbits = 0;
76 goto ret;
77 }
78 b->_wds = i + 1;
79 *bbits = i*ULbits + 32 - hi0bits(b->_x[i]);
80 ret:
81 return b;
82 }
83
84 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
85 *
86 * Inspired by "How to Print Floating-Point Numbers Accurately" by
87 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
88 *
89 * Modifications:
90 * 1. Rather than iterating, we use a simple numeric overestimate
91 * to determine k = floor(log10(d)). We scale relevant
92 * quantities using O(log2(k)) rather than O(k) multiplications.
93 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
94 * try to generate digits strictly left to right. Instead, we
95 * compute with fewer bits and propagate the carry if necessary
96 * when rounding the final digit up. This is often faster.
97 * 3. Under the assumption that input will be rounded nearest,
98 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
99 * That is, we allow equality in stopping tests when the
100 * round-nearest rule will give the same floating-point value
101 * as would satisfaction of the stopping test with strict
102 * inequality.
103 * 4. We remove common factors of powers of 2 from relevant
104 * quantities.
105 * 5. When converting floating-point integers less than 1e16,
106 * we use floating-point arithmetic rather than resorting
107 * to multiple-precision integers.
108 * 6. When asked to produce fewer than 15 digits, we first try
109 * to get by with floating-point arithmetic; we resort to
110 * multiple-precision integer arithmetic only if we cannot
111 * guarantee that the floating-point calculation has given
112 * the correctly rounded result. For k requested digits and
113 * "uniformly" distributed input, the probability is
114 * something like 10^(k-15) that we must resort to the Long
115 * calculation.
116 */
117
118 char *
119 gdtoa
120 #ifdef KR_headers
121 (ptr, fpi, be, bits, kindp, mode, ndigits, decpt, rve)
122 struct _reent *ptr, FPI *fpi; int be; ULong *bits;
123 int *kindp, mode, ndigits, *decpt; char **rve;
124 #else
125 (struct _reent *ptr, FPI *fpi, int be, ULong *bits, int *kindp,
126 int mode, int ndigits, int *decpt, char **rve)
127 #endif
128 {
129 /* Arguments ndigits and decpt are similar to the second and third
130 arguments of ecvt and fcvt; trailing zeros are suppressed from
131 the returned string. If not null, *rve is set to point
132 to the end of the return value. If d is +-Infinity or NaN,
133 then *decpt is set to 9999.
134 be = exponent: value = (integer represented by bits) * (2 to the power of be).
135
136 mode:
137 0 ==> shortest string that yields d when read in
138 and rounded to nearest.
139 1 ==> like 0, but with Steele & White stopping rule;
140 e.g. with IEEE P754 arithmetic , mode 0 gives
141 1e23 whereas mode 1 gives 9.999999999999999e22.
142 2 ==> max(1,ndigits) significant digits. This gives a
143 return value similar to that of ecvt, except
144 that trailing zeros are suppressed.
145 3 ==> through ndigits past the decimal point. This
146 gives a return value similar to that from fcvt,
147 except that trailing zeros are suppressed, and
148 ndigits can be negative.
149 4-9 should give the same return values as 2-3, i.e.,
150 4 <= mode <= 9 ==> same return as mode
151 2 + (mode & 1). These modes are mainly for
152 debugging; often they run slower but sometimes
153 faster than modes 2-3.
154 4,5,8,9 ==> left-to-right digit generation.
155 6-9 ==> don't try fast floating-point estimate
156 (if applicable).
157
158 Values of mode other than 0-9 are treated as mode 0.
159
160 Sufficient space is allocated to the return value
161 to hold the suppressed trailing zeros.
162 */
163
164 int bbits, b2, b5, be0, dig, i, ieps, ilim, ilim0, ilim1, inex;
165 int j, j1, k, k0, k_check, kind, leftright, m2, m5, nbits;
166 int rdir, s2, s5, spec_case, try_quick;
167 Long L;
168 Bigint *b, *b1, *delta, *mlo, *mhi, *mhi1, *S;
169 double d2, ds;
170 char *s, *s0;
171 U d, eps;
172
173 #ifndef MULTIPLE_THREADS
174 if (dtoa_result) {
175 freedtoa(ptr, dtoa_result);
176 dtoa_result = 0;
177 }
178 #endif
179 inex = 0;
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(ptr, "Infinity", rve, 8);
190 case STRTOG_NaN:
191 *decpt = -32768;
192 return nrv_alloc(ptr, "NaN", rve, 3);
193 default:
194 return 0;
195 }
196 b = bitstob(ptr, 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(ptr, b);
207 ret_zero:
208 *decpt = 1;
209 return nrv_alloc(ptr, "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 ( (j1 = j & 3) !=0)
261 dval(d) *= 1 << j1;
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 /* no break */
311 case 4:
312 if (ndigits <= 0)
313 ndigits = 1;
314 ilim = ilim1 = i = ndigits;
315 break;
316 case 3:
317 leftright = 0;
318 /* no break */
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(ptr, 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 = 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 >>= 1, i++)
366 if (j & 1) {
367 ieps++;
368 ds *= bigtens[i];
369 }
370 }
371 else {
372 ds = 1.;
373 if ( (j1 = -k) !=0) {
374 dval(d) *= tens[j1 & 0xf];
375 for(j = j1 >> 4; j; 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(ptr, 1);
543 if (mhi == NULL)
544 return (NULL);
545 }
546 if (m2 > 0 && s2 > 0) {
547 i = m2 < s2 ? m2 : s2;
548 b2 -= i;
549 m2 -= i;
550 s2 -= i;
551 }
552 if (b5 > 0) {
553 if (leftright) {
554 if (m5 > 0) {
555 mhi = pow5mult(ptr, mhi, m5);
556 if (mhi == NULL)
557 return (NULL);
558 b1 = mult(ptr, mhi, b);
559 if (b1 == NULL)
560 return (NULL);
561 Bfree(ptr, b);
562 b = b1;
563 }
564 if ( (j = b5 - m5) !=0) {
565 b = pow5mult(ptr, b, j);
566 if (b == NULL)
567 return (NULL);
568 }
569 }
570 else {
571 b = pow5mult(ptr, b, b5);
572 if (b == NULL)
573 return (NULL);
574 }
575 }
576 S = i2b(ptr, 1);
577 if (S == NULL)
578 return (NULL);
579 if (s5 > 0) {
580 S = pow5mult(ptr, S, s5);
581 if (S == NULL)
582 return (NULL);
583 }
584
585 /* Check for special case that d is a normalized power of 2. */
586
587 spec_case = 0;
588 if (mode < 2) {
589 if (bbits == 1 && be0 > fpi->emin + 1) {
590 /* The special case */
591 b2++;
592 s2++;
593 spec_case = 1;
594 }
595 }
596
597 /* Arrange for convenient computation of quotients:
598 * shift left if necessary so divisor has 4 leading 0 bits.
599 *
600 * Perhaps we should just compute leading 28 bits of S once
601 * and for all and pass them and a shift to quorem, so it
602 * can do shifts and ors to compute the numerator for q.
603 */
604 i = ((s5 ? hi0bits(S->_x[S->_wds-1]) : ULbits - 1) - s2 - 4) & kmask;
605 m2 += i;
606 if ((b2 += i) > 0) {
607 b = lshift(ptr, b, b2);
608 if (b == NULL)
609 return (NULL);
610 }
611 if ((s2 += i) > 0) {
612 S = lshift(ptr, S, s2);
613 if (S == NULL)
614 return (NULL);
615 }
616 if (k_check) {
617 if (cmp(b,S) < 0) {
618 k--;
619 b = multadd(ptr, b, 10, 0); /* we botched the k estimate */
620 if (b == NULL)
621 return (NULL);
622 if (leftright) {
623 mhi = multadd(ptr, mhi, 10, 0);
624 if (mhi == NULL)
625 return (NULL);
626 }
627 ilim = ilim1;
628 }
629 }
630 if (ilim <= 0 && mode > 2) {
631 S = multadd(ptr, S,5,0);
632 if (S == NULL)
633 return (NULL);
634 if (ilim < 0 || cmp(b,S) <= 0) {
635 /* no digits, fcvt style */
636 no_digits:
637 k = -1 - ndigits;
638 inex = STRTOG_Inexlo;
639 goto ret;
640 }
641 one_digit:
642 inex = STRTOG_Inexhi;
643 *s++ = '1';
644 k++;
645 goto ret;
646 }
647 if (leftright) {
648 if (m2 > 0) {
649 mhi = lshift(ptr, mhi, m2);
650 if (mhi == NULL)
651 return (NULL);
652 }
653
654 /* Compute mlo -- check for special case
655 * that d is a normalized power of 2.
656 */
657
658 mlo = mhi;
659 if (spec_case) {
660 mhi = Balloc(ptr, mhi->_k);
661 if (mhi == NULL)
662 return (NULL);
663 Bcopy(mhi, mlo);
664 mhi = lshift(ptr, mhi, 1);
665 if (mhi == NULL)
666 return (NULL);
667 }
668
669 for(i = 1;;i++) {
670 dig = quorem(b,S) + '0';
671 /* Do we yet have the shortest decimal string
672 * that will round to d?
673 */
674 j = cmp(b, mlo);
675 delta = diff(ptr, S, mhi);
676 if (delta == NULL)
677 return (NULL);
678 j1 = delta->_sign ? 1 : cmp(b, delta);
679 Bfree(ptr, delta);
680 #ifndef ROUND_BIASED
681 if (j1 == 0 && !mode && !(bits[0] & 1) && !rdir) {
682 if (dig == '9')
683 goto round_9_up;
684 if (j <= 0) {
685 if (b->_wds > 1 || b->_x[0])
686 inex = STRTOG_Inexlo;
687 }
688 else {
689 dig++;
690 inex = STRTOG_Inexhi;
691 }
692 *s++ = dig;
693 goto ret;
694 }
695 #endif
696 if (j < 0 || (j == 0 && !mode
697 #ifndef ROUND_BIASED
698 && !(bits[0] & 1)
699 #endif
700 )) {
701 if (rdir && (b->_wds > 1 || b->_x[0])) {
702 if (rdir == 2) {
703 inex = STRTOG_Inexlo;
704 goto accept;
705 }
706 while (cmp(S,mhi) > 0) {
707 *s++ = dig;
708 mhi1 = multadd(ptr, mhi, 10, 0);
709 if (mhi1 == NULL)
710 return (NULL);
711 if (mlo == mhi)
712 mlo = mhi1;
713 mhi = mhi1;
714 b = multadd(ptr, b, 10, 0);
715 if (b == NULL)
716 return (NULL);
717 dig = quorem(b,S) + '0';
718 }
719 if (dig++ == '9')
720 goto round_9_up;
721 inex = STRTOG_Inexhi;
722 goto accept;
723 }
724 if (j1 > 0) {
725 b = lshift(ptr, b, 1);
726 if (b == NULL)
727 return (NULL);
728 j1 = cmp(b, S);
729 #ifdef ROUND_BIASED
730 if (j1 >= 0 /*)*/
731 #else
732 if ((j1 > 0 || (j1 == 0 && dig & 1))
733 #endif
734 && dig++ == '9')
735 goto round_9_up;
736 inex = STRTOG_Inexhi;
737 }
738 if (b->_wds > 1 || b->_x[0])
739 inex = STRTOG_Inexlo;
740 accept:
741 *s++ = dig;
742 goto ret;
743 }
744 if (j1 > 0 && rdir != 2) {
745 if (dig == '9') { /* possible if i == 1 */
746 round_9_up:
747 *s++ = '9';
748 inex = STRTOG_Inexhi;
749 goto roundoff;
750 }
751 inex = STRTOG_Inexhi;
752 *s++ = dig + 1;
753 goto ret;
754 }
755 *s++ = dig;
756 if (i == ilim)
757 break;
758 b = multadd(ptr, b, 10, 0);
759 if (b == NULL)
760 return (NULL);
761 if (mlo == mhi) {
762 mlo = mhi = multadd(ptr, mhi, 10, 0);
763 if (mlo == NULL)
764 return (NULL);
765 }
766 else {
767 mlo = multadd(ptr, mlo, 10, 0);
768 if (mlo == NULL)
769 return (NULL);
770 mhi = multadd(ptr, mhi, 10, 0);
771 if (mhi == NULL)
772 return (NULL);
773 }
774 }
775 }
776 else
777 for(i = 1;; i++) {
778 *s++ = dig = quorem(b,S) + '0';
779 if (i >= ilim)
780 break;
781 b = multadd(ptr, b, 10, 0);
782 if (b == NULL)
783 return (NULL);
784 }
785
786 /* Round off last digit */
787
788 if (rdir) {
789 if (rdir == 2 || (b->_wds <= 1 && !b->_x[0]))
790 goto chopzeros;
791 goto roundoff;
792 }
793 b = lshift(ptr, b, 1);
794 if (b == NULL)
795 return (NULL);
796 j = cmp(b, S);
797 #ifdef ROUND_BIASED
798 if (j >= 0)
799 #else
800 if (j > 0 || (j == 0 && dig & 1))
801 #endif
802 {
803 roundoff:
804 inex = STRTOG_Inexhi;
805 while(*--s == '9')
806 if (s == s0) {
807 k++;
808 *s++ = '1';
809 goto ret;
810 }
811 ++*s++;
812 }
813 else {
814 chopzeros:
815 if (b->_wds > 1 || b->_x[0])
816 inex = STRTOG_Inexlo;
817 while(*--s == '0'){}
818 ++s;
819 }
820 ret:
821 Bfree(ptr, S);
822 if (mhi) {
823 if (mlo && mlo != mhi)
824 Bfree(ptr, mlo);
825 Bfree(ptr, mhi);
826 }
827 ret1:
828 Bfree(ptr, b);
829 *s = 0;
830 *decpt = k + 1;
831 if (rve)
832 *rve = s;
833 *kindp |= inex;
834 return s0;
835 }
836 DEF_STRONG(gdtoa);
837 #endif /* _USE_GDTOA */