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 } |