1/****************************************************************
2
3The author of this software is David M. Gay.
4
5Copyright (C) 1998, 1999 by Lucent Technologies
6All Rights Reserved
7
8Permission to use, copy, modify, and distribute this software and

--- 61 unchanged lines hidden (view full) ---

70#define Check_FLT_ROUNDS
71#else
72#define Rounding Flt_Rounds
73#endif
74
75 char *
76dtoa
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

--- 29 unchanged lines hidden (view full) ---

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;

--- 8 unchanged lines hidden (view full) ---

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 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 *sign = 0;
161
162#if defined(IEEE_Arith) + defined(VAX)
163#ifdef IEEE_Arith
164 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 if (!dval(&d)) {
182 *decpt = 1;
183 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 b = d2b(dval(&d), &be, &bbits);
201#ifdef Sudden_Underflow
202 i = (int)(word0(&d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
203#else
204 if (( i = (int)(word0(&d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)) )!=0) {
205#endif
206 dval(&d2) = dval(&d);
207 word0(&d2) &= Frac_mask1;
208 word0(&d2) |= Exp_11;
209#ifdef IBM
210 if (( j = 11 - hi0bits(word0(&d2) & Frac_mask) )!=0)
211 dval(&d2) /= 1 << j;
212#endif
213
214 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
215 * log10(x) = log(x) / log(10)
216 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
217 * log10(&d) = (i-Bias)*log(2)/log(10) + log10(&d2)
218 *
219 * This suggests computing an approximation k to log10(&d) by
220 *
221 * k = (i - Bias)*0.301029995663981
222 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
223 *
224 * We want k to be too large rather than too small.
225 * The error in the first-order Taylor series approximation
226 * is in our favor, so we just round up the constant enough
227 * to compensate for any error in the multiplication of

--- 12 unchanged lines hidden (view full) ---

240#endif
241#ifndef Sudden_Underflow
242 denorm = 0;
243 }
244 else {
245 /* d is denormalized */
246
247 i = bbits + be + (Bias + (P-1) - 1);
248 x = i > 32 ? word0(&d) << (64 - i) | word1(&d) >> (i - 32)
249 : word1(&d) << (32 - i);
250 dval(&d2) = x;
251 word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
252 i -= (Bias + (P-1) - 1) + 1;
253 denorm = 1;
254 }
255#endif
256 ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
257 k = (int)ds;
258 if (ds < 0. && ds != k)
259 k--; /* want k = floor(ds) */
260 k_check = 1;
261 if (k >= 0 && k <= Ten_pmax) {
262 if (dval(&d) < tens[k])
263 k--;
264 k_check = 0;
265 }
266 j = bbits - i - 1;
267 if (j >= 0) {
268 b2 = 0;
269 s2 = j;
270 }

--- 22 unchanged lines hidden (view full) ---

293#endif
294#endif /*SET_INEXACT*/
295
296 if (mode > 5) {
297 mode -= 4;
298 try_quick = 0;
299 }
300 leftright = 1;
301 ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
302 /* silence erroneous "gcc -Wall" warning. */
303 switch(mode) {
304 case 0:
305 case 1:
306 i = 18;
307 ndigits = 0;
308 break;
309 case 2:
310 leftright = 0;
311 /* no break */
312 case 4:
313 if (ndigits <= 0)

--- 17 unchanged lines hidden (view full) ---

331 leftright = 0;
332#endif
333
334 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
335
336 /* Try to get by with floating-point arithmetic. */
337
338 i = 0;
339 dval(&d2) = dval(&d);
340 k0 = k;
341 ilim0 = ilim;
342 ieps = 2; /* conservative */
343 if (k > 0) {
344 ds = tens[k&0xf];
345 j = k >> 4;
346 if (j & Bletch) {
347 /* prevent overflows */
348 j &= Bletch - 1;
349 dval(&d) /= bigtens[n_bigtens-1];
350 ieps++;
351 }
352 for(; j; j >>= 1, i++)
353 if (j & 1) {
354 ieps++;
355 ds *= bigtens[i];
356 }
357 dval(&d) /= ds;
358 }
359 else if (( j1 = -k )!=0) {
360 dval(&d) *= tens[j1 & 0xf];
361 for(j = j1 >> 4; j; j >>= 1, i++)
362 if (j & 1) {
363 ieps++;
364 dval(&d) *= bigtens[i];
365 }
366 }
367 if (k_check && dval(&d) < 1. && ilim > 0) {
368 if (ilim1 <= 0)
369 goto fast_failed;
370 ilim = ilim1;
371 k--;
372 dval(&d) *= 10.;
373 ieps++;
374 }
375 dval(&eps) = ieps*dval(&d) + 7.;
376 word0(&eps) -= (P-1)*Exp_msk1;
377 if (ilim == 0) {
378 S = mhi = 0;
379 dval(&d) -= 5.;
380 if (dval(&d) > dval(&eps))
381 goto one_digit;
382 if (dval(&d) < -dval(&eps))
383 goto no_digits;
384 goto fast_failed;
385 }
386#ifndef No_leftright
387 if (leftright) {
388 /* Use Steele & White method of only
389 * generating digits needed.
390 */
391 dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
392 for(i = 0;;) {
393 L = dval(&d);
394 dval(&d) -= L;
395 *s++ = '0' + (int)L;
396 if (dval(&d) < dval(&eps))
397 goto ret1;
398 if (1. - dval(&d) < dval(&eps))
399 goto bump_up;
400 if (++i >= ilim)
401 break;
402 dval(&eps) *= 10.;
403 dval(&d) *= 10.;
404 }
405 }
406 else {
407#endif
408 /* Generate ilim digits, then fix them up. */
409 dval(&eps) *= tens[ilim-1];
410 for(i = 1;; i++, dval(&d) *= 10.) {
411 L = (Long)(dval(&d));
412 if (!(dval(&d) -= L))
413 ilim = i;
414 *s++ = '0' + (int)L;
415 if (i == ilim) {
416 if (dval(&d) > 0.5 + dval(&eps))
417 goto bump_up;
418 else if (dval(&d) < 0.5 - dval(&eps)) {
419 while(*--s == '0');
420 s++;
421 goto ret1;
422 }
423 break;
424 }
425 }
426#ifndef No_leftright
427 }
428#endif
429 fast_failed:
430 s = s0;
431 dval(&d) = dval(&d2);
432 k = k0;
433 ilim = ilim0;
434 }
435
436 /* Do we have a "small" integer? */
437
438 if (be >= 0 && k <= Int_max) {
439 /* Yes. */
440 ds = tens[k];
441 if (ndigits < 0 && ilim <= 0) {
442 S = mhi = 0;
443 if (ilim < 0 || dval(&d) <= 5*ds)
444 goto no_digits;
445 goto one_digit;
446 }
447 for(i = 1;; i++, dval(&d) *= 10.) {
448 L = (Long)(dval(&d) / ds);
449 dval(&d) -= L*ds;
450#ifdef Check_FLT_ROUNDS
451 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
452 if (dval(&d) < 0) {
453 L--;
454 dval(&d) += ds;
455 }
456#endif
457 *s++ = '0' + (int)L;
458 if (!dval(&d)) {
459#ifdef SET_INEXACT
460 inexact = 0;
461#endif
462 break;
463 }
464 if (i == ilim) {
465#ifdef Honor_FLT_ROUNDS
466 if (mode > 1)
467 switch(Rounding) {
468 case 0: goto ret1;
469 case 2: goto bump_up;
470 }
471#endif
472 dval(&d) += dval(&d);
473#ifdef ROUND_BIASED
474 if (dval(&d) >= ds)
475#else
476 if (dval(&d) > ds || (dval(&d) == ds && L & 1))
477#endif
478 {
479 bump_up:
480 while(*--s == '9')
481 if (s == s0) {
482 k++;
483 *s = '0';
484 break;
485 }
486 ++*s++;

--- 48 unchanged lines hidden (view full) ---

535 /* Check for special case that d is a normalized power of 2. */
536
537 spec_case = 0;
538 if ((mode < 2 || leftright)
539#ifdef Honor_FLT_ROUNDS
540 && Rounding == 1
541#endif
542 ) {
543 if (!word1(&d) && !(word0(&d) & Bndry_mask)
544#ifndef Sudden_Underflow
545 && word0(&d) & (Exp_mask & ~Exp_msk1)
546#endif
547 ) {
548 /* The special case */
549 b2 += Log2P;
550 s2 += Log2P;
551 spec_case = 1;
552 }
553 }

--- 69 unchanged lines hidden (view full) ---

623 /* Do we yet have the shortest decimal string
624 * that will round to d?
625 */
626 j = cmp(b, mlo);
627 delta = diff(S, mhi);
628 j1 = delta->sign ? 1 : cmp(b, delta);
629 Bfree(delta);
630#ifndef ROUND_BIASED
631 if (j1 == 0 && mode != 1 && !(word1(&d) & 1)
632#ifdef Honor_FLT_ROUNDS
633 && Rounding >= 1
634#endif
635 ) {
636 if (dig == '9')
637 goto round_9_up;
638 if (j > 0)
639 dig++;
640#ifdef SET_INEXACT
641 else if (!b->x[0] && b->wds <= 1)
642 inexact = 0;
643#endif
644 *s++ = dig;
645 goto ret;
646 }
647#endif
648 if (j < 0 || (j == 0 && mode != 1
649#ifndef ROUND_BIASED
650 && !(word1(&d) & 1)
651#endif
652 )) {
653 if (!b->x[0] && b->wds <= 1) {
654#ifdef SET_INEXACT
655 inexact = 0;
656#endif
657 goto accept_dig;
658 }
659#ifdef Honor_FLT_ROUNDS
660 if (mode > 1)
661 switch(Rounding) {
662 case 0: goto accept_dig;
663 case 2: goto keep_dig;
664 }
665#endif /*Honor_FLT_ROUNDS*/
666 if (j1 > 0) {
667 b = lshift(b, 1);
668 j1 = cmp(b, S);
669#ifdef ROUND_BIASED
670 if (j1 >= 0 /*)*/
671#else
672 if ((j1 > 0 || (j1 == 0 && dig & 1))
673#endif
674 && dig++ == '9')
675 goto round_9_up;
676 }
677 accept_dig:
678 *s++ = dig;
679 goto ret;
680 }
681 if (j1 > 0) {

--- 43 unchanged lines hidden (view full) ---

725#ifdef Honor_FLT_ROUNDS
726 switch(Rounding) {
727 case 0: goto trimzeros;
728 case 2: goto roundoff;
729 }
730#endif
731 b = lshift(b, 1);
732 j = cmp(b, S);
733#ifdef ROUND_BIASED
734 if (j >= 0)
735#else
736 if (j > 0 || (j == 0 && dig & 1))
737#endif
738 {
739 roundoff:
740 while(*--s == '9')
741 if (s == s0) {
742 k++;
743 *s++ = '1';
744 goto ret;
745 }
746 ++*s++;
747 }
748 else {
749#ifdef Honor_FLT_ROUNDS
750 trimzeros:
751#endif
752 while(*--s == '0');
753 s++;
754 }
755 ret:
756 Bfree(S);
757 if (mhi) {
758 if (mlo && mlo != mhi)
759 Bfree(mlo);
760 Bfree(mhi);
761 }
762 ret1:
763#ifdef SET_INEXACT
764 if (inexact) {
765 if (!oldinexact) {
766 word0(&d) = Exp_1 + (70 << Exp_shift);
767 word1(&d) = 0;
768 dval(&d) += 1.;
769 }
770 }
771 else if (!oldinexact)
772 clear_inexact();
773#endif
774 Bfree(b);
775 *s = 0;
776 *decpt = k + 1;
777 if (rve)
778 *rve = s;
779 return s0;
780 }