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 9its documentation for any purpose and without fee is hereby 10granted, provided that the above copyright notice appear in all 11copies and that both that the copyright notice and this 12permission notice and warranty disclaimer appear in supporting 13documentation, and that the name of Lucent or any of its entities 14not be used in advertising or publicity pertaining to 15distribution of the software without specific, written prior 16permission. 17 18LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, 19INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. 20IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY 21SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 22WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER 23IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, 24ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF 25THIS 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 * 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 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 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 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 228 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, 229 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, 230 * adding 1e-13 to the constant term more than suffices. 231 * Hence we adjust the constant term to 0.1760912590558. 232 * (We could get a more accurate k by invoking log10, 233 * but this is probably not worthwhile.) 234 */ 235 236 i -= Bias; 237#ifdef IBM 238 i <<= 2; 239 i += j; 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 } 271 else { 272 b2 = -j; 273 s2 = 0; 274 } 275 if (k >= 0) { 276 b5 = 0; 277 s5 = k; 278 s2 += k; 279 } 280 else { 281 b2 -= k; 282 b5 = -k; 283 s5 = 0; 284 } 285 if (mode < 0 || mode > 9) 286 mode = 0; 287 288#ifndef SET_INEXACT 289#ifdef Check_FLT_ROUNDS 290 try_quick = Rounding == 1; 291#else 292 try_quick = 1; 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) 314 ndigits = 1; 315 ilim = ilim1 = i = ndigits; 316 break; 317 case 3: 318 leftright = 0; 319 /* no break */ 320 case 5: 321 i = ndigits + k + 1; 322 ilim = i; 323 ilim1 = i - 1; 324 if (i <= 0) 325 i = 1; 326 } 327 s = s0 = rv_alloc(i); 328 329#ifdef Honor_FLT_ROUNDS 330 if (mode > 1 && Rounding != 1) 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 if (dval(&d) > ds || (dval(&d) == ds && L & 1)) { 474 bump_up: 475 while(*--s == '9') 476 if (s == s0) { 477 k++; 478 *s = '0'; 479 break; 480 } 481 ++*s++; 482 } 483 break; 484 } 485 } 486 goto ret1; 487 } 488 489 m2 = b2; 490 m5 = b5; 491 mhi = mlo = 0; 492 if (leftright) { 493 i = 494#ifndef Sudden_Underflow 495 denorm ? be + (Bias + (P-1) - 1 + 1) : 496#endif 497#ifdef IBM 498 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3); 499#else 500 1 + P - bbits; 501#endif 502 b2 += i; 503 s2 += i; 504 mhi = i2b(1); 505 } 506 if (m2 > 0 && s2 > 0) { 507 i = m2 < s2 ? m2 : s2; 508 b2 -= i; 509 m2 -= i; 510 s2 -= i; 511 } 512 if (b5 > 0) { 513 if (leftright) { 514 if (m5 > 0) { 515 mhi = pow5mult(mhi, m5); 516 b1 = mult(mhi, b); 517 Bfree(b); 518 b = b1; 519 } 520 if (( j = b5 - m5 )!=0) 521 b = pow5mult(b, j); 522 } 523 else 524 b = pow5mult(b, b5); 525 } 526 S = i2b(1); 527 if (s5 > 0) 528 S = pow5mult(S, s5); 529 530 /* Check for special case that d is a normalized power of 2. */ 531 532 spec_case = 0; 533 if ((mode < 2 || leftright) 534#ifdef Honor_FLT_ROUNDS 535 && Rounding == 1 536#endif 537 ) { 538 if (!word1(&d) && !(word0(&d) & Bndry_mask) 539#ifndef Sudden_Underflow 540 && word0(&d) & (Exp_mask & ~Exp_msk1) 541#endif 542 ) { 543 /* The special case */ 544 b2 += Log2P; 545 s2 += Log2P; 546 spec_case = 1; 547 } 548 } 549 550 /* Arrange for convenient computation of quotients: 551 * shift left if necessary so divisor has 4 leading 0 bits. 552 * 553 * Perhaps we should just compute leading 28 bits of S once 554 * and for all and pass them and a shift to quorem, so it 555 * can do shifts and ors to compute the numerator for q. 556 */ 557#ifdef Pack_32 558 if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f )!=0) 559 i = 32 - i; 560#else 561 if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf )!=0) 562 i = 16 - i; 563#endif 564 if (i > 4) { 565 i -= 4; 566 b2 += i; 567 m2 += i; 568 s2 += i; 569 } 570 else if (i < 4) { 571 i += 28; 572 b2 += i; 573 m2 += i; 574 s2 += i; 575 } 576 if (b2 > 0) 577 b = lshift(b, b2); 578 if (s2 > 0) 579 S = lshift(S, s2); 580 if (k_check) { 581 if (cmp(b,S) < 0) { 582 k--; 583 b = multadd(b, 10, 0); /* we botched the k estimate */ 584 if (leftright) 585 mhi = multadd(mhi, 10, 0); 586 ilim = ilim1; 587 } 588 } 589 if (ilim <= 0 && (mode == 3 || mode == 5)) { 590 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) { 591 /* no digits, fcvt style */ 592 no_digits: 593 k = -1 - ndigits; 594 goto ret; 595 } 596 one_digit: 597 *s++ = '1'; 598 k++; 599 goto ret; 600 } 601 if (leftright) { 602 if (m2 > 0) 603 mhi = lshift(mhi, m2); 604 605 /* Compute mlo -- check for special case 606 * that d is a normalized power of 2. 607 */ 608 609 mlo = mhi; 610 if (spec_case) { 611 mhi = Balloc(mhi->k); 612 Bcopy(mhi, mlo); 613 mhi = lshift(mhi, Log2P); 614 } 615 616 for(i = 1;;i++) { 617 dig = quorem(b,S) + '0'; 618 /* Do we yet have the shortest decimal string 619 * that will round to d? 620 */ 621 j = cmp(b, mlo); 622 delta = diff(S, mhi); 623 j1 = delta->sign ? 1 : cmp(b, delta); 624 Bfree(delta); 625#ifndef ROUND_BIASED 626 if (j1 == 0 && mode != 1 && !(word1(&d) & 1) 627#ifdef Honor_FLT_ROUNDS 628 && Rounding >= 1 629#endif 630 ) { 631 if (dig == '9') 632 goto round_9_up; 633 if (j > 0) 634 dig++; 635#ifdef SET_INEXACT 636 else if (!b->x[0] && b->wds <= 1) 637 inexact = 0; 638#endif 639 *s++ = dig; 640 goto ret; 641 } 642#endif 643 if (j < 0 || (j == 0 && mode != 1 644#ifndef ROUND_BIASED 645 && !(word1(&d) & 1) 646#endif 647 )) { 648 if (!b->x[0] && b->wds <= 1) { 649#ifdef SET_INEXACT 650 inexact = 0; 651#endif 652 goto accept_dig; 653 } 654#ifdef Honor_FLT_ROUNDS 655 if (mode > 1) 656 switch(Rounding) { 657 case 0: goto accept_dig; 658 case 2: goto keep_dig; 659 } 660#endif /*Honor_FLT_ROUNDS*/ 661 if (j1 > 0) { 662 b = lshift(b, 1); 663 j1 = cmp(b, S); 664 if ((j1 > 0 || (j1 == 0 && dig & 1)) 665 && dig++ == '9') 666 goto round_9_up; 667 } 668 accept_dig: 669 *s++ = dig; 670 goto ret; 671 } 672 if (j1 > 0) { 673#ifdef Honor_FLT_ROUNDS 674 if (!Rounding) 675 goto accept_dig; 676#endif 677 if (dig == '9') { /* possible if i == 1 */ 678 round_9_up: 679 *s++ = '9'; 680 goto roundoff; 681 } 682 *s++ = dig + 1; 683 goto ret; 684 } 685#ifdef Honor_FLT_ROUNDS 686 keep_dig: 687#endif 688 *s++ = dig; 689 if (i == ilim) 690 break; 691 b = multadd(b, 10, 0); 692 if (mlo == mhi) 693 mlo = mhi = multadd(mhi, 10, 0); 694 else { 695 mlo = multadd(mlo, 10, 0); 696 mhi = multadd(mhi, 10, 0); 697 } 698 } 699 } 700 else 701 for(i = 1;; i++) { 702 *s++ = dig = quorem(b,S) + '0'; 703 if (!b->x[0] && b->wds <= 1) { 704#ifdef SET_INEXACT 705 inexact = 0; 706#endif 707 goto ret; 708 } 709 if (i >= ilim) 710 break; 711 b = multadd(b, 10, 0); 712 } 713 714 /* Round off last digit */ 715 716#ifdef Honor_FLT_ROUNDS 717 switch(Rounding) { 718 case 0: goto trimzeros; 719 case 2: goto roundoff; 720 } 721#endif 722 b = lshift(b, 1); 723 j = cmp(b, S); 724 if (j > 0 || (j == 0 && dig & 1)) { 725 roundoff: 726 while(*--s == '9') 727 if (s == s0) { 728 k++; 729 *s++ = '1'; 730 goto ret; 731 } 732 ++*s++; 733 } 734 else { 735#ifdef Honor_FLT_ROUNDS 736 trimzeros: 737#endif 738 while(*--s == '0'); 739 s++; 740 } 741 ret: 742 Bfree(S); 743 if (mhi) { 744 if (mlo && mlo != mhi) 745 Bfree(mlo); 746 Bfree(mhi); 747 } 748 ret1: 749#ifdef SET_INEXACT 750 if (inexact) { 751 if (!oldinexact) { 752 word0(&d) = Exp_1 + (70 << Exp_shift); 753 word1(&d) = 0; 754 dval(&d) += 1.; 755 } 756 } 757 else if (!oldinexact) 758 clear_inexact(); 759#endif 760 Bfree(b); 761 *s = 0; 762 *decpt = k + 1; 763 if (rve) 764 *rve = s; 765 return s0; 766 } 767