1112158Sdas/**************************************************************** 2112158Sdas 3112158SdasThe author of this software is David M. Gay. 4112158Sdas 5112158SdasCopyright (C) 1998, 1999 by Lucent Technologies 6112158SdasAll Rights Reserved 7112158Sdas 8112158SdasPermission to use, copy, modify, and distribute this software and 9112158Sdasits documentation for any purpose and without fee is hereby 10112158Sdasgranted, provided that the above copyright notice appear in all 11112158Sdascopies and that both that the copyright notice and this 12112158Sdaspermission notice and warranty disclaimer appear in supporting 13112158Sdasdocumentation, and that the name of Lucent or any of its entities 14112158Sdasnot be used in advertising or publicity pertaining to 15112158Sdasdistribution of the software without specific, written prior 16112158Sdaspermission. 17112158Sdas 18112158SdasLUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, 19112158SdasINCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. 20112158SdasIN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY 21112158SdasSPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 22112158SdasWHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER 23112158SdasIN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, 24112158SdasARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF 25112158SdasTHIS SOFTWARE. 26112158Sdas 27112158Sdas****************************************************************/ 28112158Sdas 29165743Sdas/* Please send bug reports to David M. Gay (dmg at acm dot org, 30165743Sdas * with " at " changed at "@" and " dot " changed to "."). */ 31112158Sdas 32112158Sdas#include "gdtoaimp.h" 33112158Sdas 34112158Sdas/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. 35112158Sdas * 36112158Sdas * Inspired by "How to Print Floating-Point Numbers Accurately" by 37165743Sdas * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126]. 38112158Sdas * 39112158Sdas * Modifications: 40112158Sdas * 1. Rather than iterating, we use a simple numeric overestimate 41112158Sdas * to determine k = floor(log10(d)). We scale relevant 42112158Sdas * quantities using O(log2(k)) rather than O(k) multiplications. 43112158Sdas * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't 44112158Sdas * try to generate digits strictly left to right. Instead, we 45112158Sdas * compute with fewer bits and propagate the carry if necessary 46112158Sdas * when rounding the final digit up. This is often faster. 47112158Sdas * 3. Under the assumption that input will be rounded nearest, 48112158Sdas * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. 49112158Sdas * That is, we allow equality in stopping tests when the 50112158Sdas * round-nearest rule will give the same floating-point value 51112158Sdas * as would satisfaction of the stopping test with strict 52112158Sdas * inequality. 53112158Sdas * 4. We remove common factors of powers of 2 from relevant 54112158Sdas * quantities. 55112158Sdas * 5. When converting floating-point integers less than 1e16, 56112158Sdas * we use floating-point arithmetic rather than resorting 57112158Sdas * to multiple-precision integers. 58112158Sdas * 6. When asked to produce fewer than 15 digits, we first try 59112158Sdas * to get by with floating-point arithmetic; we resort to 60112158Sdas * multiple-precision integer arithmetic only if we cannot 61112158Sdas * guarantee that the floating-point calculation has given 62112158Sdas * the correctly rounded result. For k requested digits and 63112158Sdas * "uniformly" distributed input, the probability is 64112158Sdas * something like 10^(k-15) that we must resort to the Long 65112158Sdas * calculation. 66112158Sdas */ 67112158Sdas 68112158Sdas#ifdef Honor_FLT_ROUNDS 69112158Sdas#undef Check_FLT_ROUNDS 70112158Sdas#define Check_FLT_ROUNDS 71112158Sdas#else 72112158Sdas#define Rounding Flt_Rounds 73112158Sdas#endif 74112158Sdas 75112158Sdas char * 76112158Sdasdtoa 77112158Sdas#ifdef KR_headers 78219557Sdas (d0, mode, ndigits, decpt, sign, rve) 79219557Sdas double d0; int mode, ndigits, *decpt, *sign; char **rve; 80112158Sdas#else 81219557Sdas (double d0, int mode, int ndigits, int *decpt, int *sign, char **rve) 82112158Sdas#endif 83112158Sdas{ 84112158Sdas /* Arguments ndigits, decpt, sign are similar to those 85112158Sdas of ecvt and fcvt; trailing zeros are suppressed from 86112158Sdas the returned string. If not null, *rve is set to point 87112158Sdas to the end of the return value. If d is +-Infinity or NaN, 88112158Sdas then *decpt is set to 9999. 89112158Sdas 90112158Sdas mode: 91112158Sdas 0 ==> shortest string that yields d when read in 92112158Sdas and rounded to nearest. 93112158Sdas 1 ==> like 0, but with Steele & White stopping rule; 94112158Sdas e.g. with IEEE P754 arithmetic , mode 0 gives 95112158Sdas 1e23 whereas mode 1 gives 9.999999999999999e22. 96112158Sdas 2 ==> max(1,ndigits) significant digits. This gives a 97112158Sdas return value similar to that of ecvt, except 98112158Sdas that trailing zeros are suppressed. 99112158Sdas 3 ==> through ndigits past the decimal point. This 100112158Sdas gives a return value similar to that from fcvt, 101112158Sdas except that trailing zeros are suppressed, and 102112158Sdas ndigits can be negative. 103112158Sdas 4,5 ==> similar to 2 and 3, respectively, but (in 104112158Sdas round-nearest mode) with the tests of mode 0 to 105112158Sdas possibly return a shorter string that rounds to d. 106112158Sdas With IEEE arithmetic and compilation with 107112158Sdas -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same 108112158Sdas as modes 2 and 3 when FLT_ROUNDS != 1. 109112158Sdas 6-9 ==> Debugging modes similar to mode - 4: don't try 110112158Sdas fast floating-point estimate (if applicable). 111112158Sdas 112112158Sdas Values of mode other than 0-9 are treated as mode 0. 113112158Sdas 114112158Sdas Sufficient space is allocated to the return value 115112158Sdas to hold the suppressed trailing zeros. 116112158Sdas */ 117112158Sdas 118112158Sdas int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1, 119112158Sdas j, j1, k, k0, k_check, leftright, m2, m5, s2, s5, 120112158Sdas spec_case, try_quick; 121112158Sdas Long L; 122112158Sdas#ifndef Sudden_Underflow 123112158Sdas int denorm; 124112158Sdas ULong x; 125112158Sdas#endif 126112158Sdas Bigint *b, *b1, *delta, *mlo, *mhi, *S; 127219557Sdas U d, d2, eps; 128219557Sdas double ds; 129112158Sdas char *s, *s0; 130112158Sdas#ifdef SET_INEXACT 131112158Sdas int inexact, oldinexact; 132112158Sdas#endif 133182709Sdas#ifdef Honor_FLT_ROUNDS /*{*/ 134182709Sdas int Rounding; 135182709Sdas#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */ 136182709Sdas Rounding = Flt_Rounds; 137182709Sdas#else /*}{*/ 138182709Sdas Rounding = 1; 139182709Sdas switch(fegetround()) { 140182709Sdas case FE_TOWARDZERO: Rounding = 0; break; 141182709Sdas case FE_UPWARD: Rounding = 2; break; 142182709Sdas case FE_DOWNWARD: Rounding = 3; 143182709Sdas } 144182709Sdas#endif /*}}*/ 145182709Sdas#endif /*}*/ 146112158Sdas 147112158Sdas#ifndef MULTIPLE_THREADS 148112158Sdas if (dtoa_result) { 149112158Sdas freedtoa(dtoa_result); 150112158Sdas dtoa_result = 0; 151112158Sdas } 152112158Sdas#endif 153219557Sdas d.d = d0; 154219557Sdas if (word0(&d) & Sign_bit) { 155112158Sdas /* set sign for everything, including 0's and NaNs */ 156112158Sdas *sign = 1; 157219557Sdas word0(&d) &= ~Sign_bit; /* clear sign bit */ 158112158Sdas } 159112158Sdas else 160112158Sdas *sign = 0; 161112158Sdas 162112158Sdas#if defined(IEEE_Arith) + defined(VAX) 163112158Sdas#ifdef IEEE_Arith 164219557Sdas if ((word0(&d) & Exp_mask) == Exp_mask) 165112158Sdas#else 166219557Sdas if (word0(&d) == 0x8000) 167112158Sdas#endif 168112158Sdas { 169112158Sdas /* Infinity or NaN */ 170112158Sdas *decpt = 9999; 171112158Sdas#ifdef IEEE_Arith 172219557Sdas if (!word1(&d) && !(word0(&d) & 0xfffff)) 173112158Sdas return nrv_alloc("Infinity", rve, 8); 174112158Sdas#endif 175112158Sdas return nrv_alloc("NaN", rve, 3); 176112158Sdas } 177112158Sdas#endif 178112158Sdas#ifdef IBM 179219557Sdas dval(&d) += 0; /* normalize */ 180112158Sdas#endif 181219557Sdas if (!dval(&d)) { 182112158Sdas *decpt = 1; 183112158Sdas return nrv_alloc("0", rve, 1); 184112158Sdas } 185112158Sdas 186112158Sdas#ifdef SET_INEXACT 187112158Sdas try_quick = oldinexact = get_inexact(); 188112158Sdas inexact = 1; 189112158Sdas#endif 190112158Sdas#ifdef Honor_FLT_ROUNDS 191182709Sdas if (Rounding >= 2) { 192112158Sdas if (*sign) 193182709Sdas Rounding = Rounding == 2 ? 0 : 2; 194112158Sdas else 195182709Sdas if (Rounding != 2) 196182709Sdas Rounding = 0; 197112158Sdas } 198112158Sdas#endif 199112158Sdas 200219557Sdas b = d2b(dval(&d), &be, &bbits); 201112158Sdas#ifdef Sudden_Underflow 202219557Sdas i = (int)(word0(&d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)); 203112158Sdas#else 204219557Sdas if (( i = (int)(word0(&d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)) )!=0) { 205112158Sdas#endif 206219557Sdas dval(&d2) = dval(&d); 207219557Sdas word0(&d2) &= Frac_mask1; 208219557Sdas word0(&d2) |= Exp_11; 209112158Sdas#ifdef IBM 210219557Sdas if (( j = 11 - hi0bits(word0(&d2) & Frac_mask) )!=0) 211219557Sdas dval(&d2) /= 1 << j; 212112158Sdas#endif 213112158Sdas 214112158Sdas /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 215112158Sdas * log10(x) = log(x) / log(10) 216112158Sdas * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) 217219557Sdas * log10(&d) = (i-Bias)*log(2)/log(10) + log10(&d2) 218112158Sdas * 219219557Sdas * This suggests computing an approximation k to log10(&d) by 220112158Sdas * 221112158Sdas * k = (i - Bias)*0.301029995663981 222112158Sdas * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); 223112158Sdas * 224112158Sdas * We want k to be too large rather than too small. 225112158Sdas * The error in the first-order Taylor series approximation 226112158Sdas * is in our favor, so we just round up the constant enough 227112158Sdas * to compensate for any error in the multiplication of 228112158Sdas * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, 229112158Sdas * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, 230112158Sdas * adding 1e-13 to the constant term more than suffices. 231112158Sdas * Hence we adjust the constant term to 0.1760912590558. 232112158Sdas * (We could get a more accurate k by invoking log10, 233112158Sdas * but this is probably not worthwhile.) 234112158Sdas */ 235112158Sdas 236112158Sdas i -= Bias; 237112158Sdas#ifdef IBM 238112158Sdas i <<= 2; 239112158Sdas i += j; 240112158Sdas#endif 241112158Sdas#ifndef Sudden_Underflow 242112158Sdas denorm = 0; 243112158Sdas } 244112158Sdas else { 245112158Sdas /* d is denormalized */ 246112158Sdas 247112158Sdas i = bbits + be + (Bias + (P-1) - 1); 248219557Sdas x = i > 32 ? word0(&d) << (64 - i) | word1(&d) >> (i - 32) 249219557Sdas : word1(&d) << (32 - i); 250219557Sdas dval(&d2) = x; 251219557Sdas word0(&d2) -= 31*Exp_msk1; /* adjust exponent */ 252112158Sdas i -= (Bias + (P-1) - 1) + 1; 253112158Sdas denorm = 1; 254112158Sdas } 255112158Sdas#endif 256219557Sdas ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981; 257112158Sdas k = (int)ds; 258112158Sdas if (ds < 0. && ds != k) 259112158Sdas k--; /* want k = floor(ds) */ 260112158Sdas k_check = 1; 261112158Sdas if (k >= 0 && k <= Ten_pmax) { 262219557Sdas if (dval(&d) < tens[k]) 263112158Sdas k--; 264112158Sdas k_check = 0; 265112158Sdas } 266112158Sdas j = bbits - i - 1; 267112158Sdas if (j >= 0) { 268112158Sdas b2 = 0; 269112158Sdas s2 = j; 270112158Sdas } 271112158Sdas else { 272112158Sdas b2 = -j; 273112158Sdas s2 = 0; 274112158Sdas } 275112158Sdas if (k >= 0) { 276112158Sdas b5 = 0; 277112158Sdas s5 = k; 278112158Sdas s2 += k; 279112158Sdas } 280112158Sdas else { 281112158Sdas b2 -= k; 282112158Sdas b5 = -k; 283112158Sdas s5 = 0; 284112158Sdas } 285112158Sdas if (mode < 0 || mode > 9) 286112158Sdas mode = 0; 287112158Sdas 288112158Sdas#ifndef SET_INEXACT 289112158Sdas#ifdef Check_FLT_ROUNDS 290112158Sdas try_quick = Rounding == 1; 291112158Sdas#else 292112158Sdas try_quick = 1; 293112158Sdas#endif 294112158Sdas#endif /*SET_INEXACT*/ 295112158Sdas 296112158Sdas if (mode > 5) { 297112158Sdas mode -= 4; 298112158Sdas try_quick = 0; 299112158Sdas } 300112158Sdas leftright = 1; 301219557Sdas ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */ 302219557Sdas /* silence erroneous "gcc -Wall" warning. */ 303112158Sdas switch(mode) { 304112158Sdas case 0: 305112158Sdas case 1: 306112158Sdas i = 18; 307112158Sdas ndigits = 0; 308112158Sdas break; 309112158Sdas case 2: 310112158Sdas leftright = 0; 311112158Sdas /* no break */ 312112158Sdas case 4: 313112158Sdas if (ndigits <= 0) 314112158Sdas ndigits = 1; 315112158Sdas ilim = ilim1 = i = ndigits; 316112158Sdas break; 317112158Sdas case 3: 318112158Sdas leftright = 0; 319112158Sdas /* no break */ 320112158Sdas case 5: 321112158Sdas i = ndigits + k + 1; 322112158Sdas ilim = i; 323112158Sdas ilim1 = i - 1; 324112158Sdas if (i <= 0) 325112158Sdas i = 1; 326112158Sdas } 327112158Sdas s = s0 = rv_alloc(i); 328112158Sdas 329112158Sdas#ifdef Honor_FLT_ROUNDS 330182709Sdas if (mode > 1 && Rounding != 1) 331112158Sdas leftright = 0; 332112158Sdas#endif 333112158Sdas 334112158Sdas if (ilim >= 0 && ilim <= Quick_max && try_quick) { 335112158Sdas 336112158Sdas /* Try to get by with floating-point arithmetic. */ 337112158Sdas 338112158Sdas i = 0; 339219557Sdas dval(&d2) = dval(&d); 340112158Sdas k0 = k; 341112158Sdas ilim0 = ilim; 342112158Sdas ieps = 2; /* conservative */ 343112158Sdas if (k > 0) { 344112158Sdas ds = tens[k&0xf]; 345112158Sdas j = k >> 4; 346112158Sdas if (j & Bletch) { 347112158Sdas /* prevent overflows */ 348112158Sdas j &= Bletch - 1; 349219557Sdas dval(&d) /= bigtens[n_bigtens-1]; 350112158Sdas ieps++; 351112158Sdas } 352112158Sdas for(; j; j >>= 1, i++) 353112158Sdas if (j & 1) { 354112158Sdas ieps++; 355112158Sdas ds *= bigtens[i]; 356112158Sdas } 357219557Sdas dval(&d) /= ds; 358112158Sdas } 359112158Sdas else if (( j1 = -k )!=0) { 360219557Sdas dval(&d) *= tens[j1 & 0xf]; 361112158Sdas for(j = j1 >> 4; j; j >>= 1, i++) 362112158Sdas if (j & 1) { 363112158Sdas ieps++; 364219557Sdas dval(&d) *= bigtens[i]; 365112158Sdas } 366112158Sdas } 367219557Sdas if (k_check && dval(&d) < 1. && ilim > 0) { 368112158Sdas if (ilim1 <= 0) 369112158Sdas goto fast_failed; 370112158Sdas ilim = ilim1; 371112158Sdas k--; 372219557Sdas dval(&d) *= 10.; 373112158Sdas ieps++; 374112158Sdas } 375219557Sdas dval(&eps) = ieps*dval(&d) + 7.; 376219557Sdas word0(&eps) -= (P-1)*Exp_msk1; 377112158Sdas if (ilim == 0) { 378112158Sdas S = mhi = 0; 379219557Sdas dval(&d) -= 5.; 380219557Sdas if (dval(&d) > dval(&eps)) 381112158Sdas goto one_digit; 382219557Sdas if (dval(&d) < -dval(&eps)) 383112158Sdas goto no_digits; 384112158Sdas goto fast_failed; 385112158Sdas } 386112158Sdas#ifndef No_leftright 387112158Sdas if (leftright) { 388112158Sdas /* Use Steele & White method of only 389112158Sdas * generating digits needed. 390112158Sdas */ 391219557Sdas dval(&eps) = 0.5/tens[ilim-1] - dval(&eps); 392112158Sdas for(i = 0;;) { 393219557Sdas L = dval(&d); 394219557Sdas dval(&d) -= L; 395112158Sdas *s++ = '0' + (int)L; 396219557Sdas if (dval(&d) < dval(&eps)) 397112158Sdas goto ret1; 398219557Sdas if (1. - dval(&d) < dval(&eps)) 399112158Sdas goto bump_up; 400112158Sdas if (++i >= ilim) 401112158Sdas break; 402219557Sdas dval(&eps) *= 10.; 403219557Sdas dval(&d) *= 10.; 404112158Sdas } 405112158Sdas } 406112158Sdas else { 407112158Sdas#endif 408112158Sdas /* Generate ilim digits, then fix them up. */ 409219557Sdas dval(&eps) *= tens[ilim-1]; 410219557Sdas for(i = 1;; i++, dval(&d) *= 10.) { 411219557Sdas L = (Long)(dval(&d)); 412219557Sdas if (!(dval(&d) -= L)) 413112158Sdas ilim = i; 414112158Sdas *s++ = '0' + (int)L; 415112158Sdas if (i == ilim) { 416219557Sdas if (dval(&d) > 0.5 + dval(&eps)) 417112158Sdas goto bump_up; 418219557Sdas else if (dval(&d) < 0.5 - dval(&eps)) { 419112158Sdas while(*--s == '0'); 420112158Sdas s++; 421112158Sdas goto ret1; 422112158Sdas } 423112158Sdas break; 424112158Sdas } 425112158Sdas } 426112158Sdas#ifndef No_leftright 427112158Sdas } 428112158Sdas#endif 429112158Sdas fast_failed: 430112158Sdas s = s0; 431219557Sdas dval(&d) = dval(&d2); 432112158Sdas k = k0; 433112158Sdas ilim = ilim0; 434112158Sdas } 435112158Sdas 436112158Sdas /* Do we have a "small" integer? */ 437112158Sdas 438112158Sdas if (be >= 0 && k <= Int_max) { 439112158Sdas /* Yes. */ 440112158Sdas ds = tens[k]; 441112158Sdas if (ndigits < 0 && ilim <= 0) { 442112158Sdas S = mhi = 0; 443219557Sdas if (ilim < 0 || dval(&d) <= 5*ds) 444112158Sdas goto no_digits; 445112158Sdas goto one_digit; 446112158Sdas } 447219557Sdas for(i = 1;; i++, dval(&d) *= 10.) { 448219557Sdas L = (Long)(dval(&d) / ds); 449219557Sdas dval(&d) -= L*ds; 450112158Sdas#ifdef Check_FLT_ROUNDS 451112158Sdas /* If FLT_ROUNDS == 2, L will usually be high by 1 */ 452219557Sdas if (dval(&d) < 0) { 453112158Sdas L--; 454219557Sdas dval(&d) += ds; 455112158Sdas } 456112158Sdas#endif 457112158Sdas *s++ = '0' + (int)L; 458219557Sdas if (!dval(&d)) { 459112158Sdas#ifdef SET_INEXACT 460112158Sdas inexact = 0; 461112158Sdas#endif 462112158Sdas break; 463112158Sdas } 464112158Sdas if (i == ilim) { 465112158Sdas#ifdef Honor_FLT_ROUNDS 466112158Sdas if (mode > 1) 467182709Sdas switch(Rounding) { 468112158Sdas case 0: goto ret1; 469112158Sdas case 2: goto bump_up; 470112158Sdas } 471112158Sdas#endif 472219557Sdas dval(&d) += dval(&d); 473219557Sdas#ifdef ROUND_BIASED 474219557Sdas if (dval(&d) >= ds) 475219557Sdas#else 476219557Sdas if (dval(&d) > ds || (dval(&d) == ds && L & 1)) 477219557Sdas#endif 478219557Sdas { 479112158Sdas bump_up: 480112158Sdas while(*--s == '9') 481112158Sdas if (s == s0) { 482112158Sdas k++; 483112158Sdas *s = '0'; 484112158Sdas break; 485112158Sdas } 486112158Sdas ++*s++; 487112158Sdas } 488112158Sdas break; 489112158Sdas } 490112158Sdas } 491112158Sdas goto ret1; 492112158Sdas } 493112158Sdas 494112158Sdas m2 = b2; 495112158Sdas m5 = b5; 496112158Sdas mhi = mlo = 0; 497112158Sdas if (leftright) { 498112158Sdas i = 499112158Sdas#ifndef Sudden_Underflow 500112158Sdas denorm ? be + (Bias + (P-1) - 1 + 1) : 501112158Sdas#endif 502112158Sdas#ifdef IBM 503112158Sdas 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3); 504112158Sdas#else 505112158Sdas 1 + P - bbits; 506112158Sdas#endif 507112158Sdas b2 += i; 508112158Sdas s2 += i; 509112158Sdas mhi = i2b(1); 510112158Sdas } 511112158Sdas if (m2 > 0 && s2 > 0) { 512112158Sdas i = m2 < s2 ? m2 : s2; 513112158Sdas b2 -= i; 514112158Sdas m2 -= i; 515112158Sdas s2 -= i; 516112158Sdas } 517112158Sdas if (b5 > 0) { 518112158Sdas if (leftright) { 519112158Sdas if (m5 > 0) { 520112158Sdas mhi = pow5mult(mhi, m5); 521112158Sdas b1 = mult(mhi, b); 522112158Sdas Bfree(b); 523112158Sdas b = b1; 524112158Sdas } 525112158Sdas if (( j = b5 - m5 )!=0) 526112158Sdas b = pow5mult(b, j); 527112158Sdas } 528112158Sdas else 529112158Sdas b = pow5mult(b, b5); 530112158Sdas } 531112158Sdas S = i2b(1); 532112158Sdas if (s5 > 0) 533112158Sdas S = pow5mult(S, s5); 534112158Sdas 535112158Sdas /* Check for special case that d is a normalized power of 2. */ 536112158Sdas 537112158Sdas spec_case = 0; 538112158Sdas if ((mode < 2 || leftright) 539112158Sdas#ifdef Honor_FLT_ROUNDS 540182709Sdas && Rounding == 1 541112158Sdas#endif 542112158Sdas ) { 543219557Sdas if (!word1(&d) && !(word0(&d) & Bndry_mask) 544112158Sdas#ifndef Sudden_Underflow 545219557Sdas && word0(&d) & (Exp_mask & ~Exp_msk1) 546112158Sdas#endif 547112158Sdas ) { 548112158Sdas /* The special case */ 549112158Sdas b2 += Log2P; 550112158Sdas s2 += Log2P; 551112158Sdas spec_case = 1; 552112158Sdas } 553112158Sdas } 554112158Sdas 555112158Sdas /* Arrange for convenient computation of quotients: 556112158Sdas * shift left if necessary so divisor has 4 leading 0 bits. 557112158Sdas * 558112158Sdas * Perhaps we should just compute leading 28 bits of S once 559112158Sdas * and for all and pass them and a shift to quorem, so it 560112158Sdas * can do shifts and ors to compute the numerator for q. 561112158Sdas */ 562112158Sdas#ifdef Pack_32 563112158Sdas if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f )!=0) 564112158Sdas i = 32 - i; 565112158Sdas#else 566112158Sdas if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf )!=0) 567112158Sdas i = 16 - i; 568112158Sdas#endif 569112158Sdas if (i > 4) { 570112158Sdas i -= 4; 571112158Sdas b2 += i; 572112158Sdas m2 += i; 573112158Sdas s2 += i; 574112158Sdas } 575112158Sdas else if (i < 4) { 576112158Sdas i += 28; 577112158Sdas b2 += i; 578112158Sdas m2 += i; 579112158Sdas s2 += i; 580112158Sdas } 581112158Sdas if (b2 > 0) 582112158Sdas b = lshift(b, b2); 583112158Sdas if (s2 > 0) 584112158Sdas S = lshift(S, s2); 585112158Sdas if (k_check) { 586112158Sdas if (cmp(b,S) < 0) { 587112158Sdas k--; 588112158Sdas b = multadd(b, 10, 0); /* we botched the k estimate */ 589112158Sdas if (leftright) 590112158Sdas mhi = multadd(mhi, 10, 0); 591112158Sdas ilim = ilim1; 592112158Sdas } 593112158Sdas } 594112158Sdas if (ilim <= 0 && (mode == 3 || mode == 5)) { 595112158Sdas if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) { 596112158Sdas /* no digits, fcvt style */ 597112158Sdas no_digits: 598112158Sdas k = -1 - ndigits; 599112158Sdas goto ret; 600112158Sdas } 601112158Sdas one_digit: 602112158Sdas *s++ = '1'; 603112158Sdas k++; 604112158Sdas goto ret; 605112158Sdas } 606112158Sdas if (leftright) { 607112158Sdas if (m2 > 0) 608112158Sdas mhi = lshift(mhi, m2); 609112158Sdas 610112158Sdas /* Compute mlo -- check for special case 611112158Sdas * that d is a normalized power of 2. 612112158Sdas */ 613112158Sdas 614112158Sdas mlo = mhi; 615112158Sdas if (spec_case) { 616112158Sdas mhi = Balloc(mhi->k); 617112158Sdas Bcopy(mhi, mlo); 618112158Sdas mhi = lshift(mhi, Log2P); 619112158Sdas } 620112158Sdas 621112158Sdas for(i = 1;;i++) { 622112158Sdas dig = quorem(b,S) + '0'; 623112158Sdas /* Do we yet have the shortest decimal string 624112158Sdas * that will round to d? 625112158Sdas */ 626112158Sdas j = cmp(b, mlo); 627112158Sdas delta = diff(S, mhi); 628112158Sdas j1 = delta->sign ? 1 : cmp(b, delta); 629112158Sdas Bfree(delta); 630112158Sdas#ifndef ROUND_BIASED 631219557Sdas if (j1 == 0 && mode != 1 && !(word1(&d) & 1) 632112158Sdas#ifdef Honor_FLT_ROUNDS 633182709Sdas && Rounding >= 1 634112158Sdas#endif 635112158Sdas ) { 636112158Sdas if (dig == '9') 637112158Sdas goto round_9_up; 638112158Sdas if (j > 0) 639112158Sdas dig++; 640112158Sdas#ifdef SET_INEXACT 641112158Sdas else if (!b->x[0] && b->wds <= 1) 642112158Sdas inexact = 0; 643112158Sdas#endif 644112158Sdas *s++ = dig; 645112158Sdas goto ret; 646112158Sdas } 647112158Sdas#endif 648219557Sdas if (j < 0 || (j == 0 && mode != 1 649112158Sdas#ifndef ROUND_BIASED 650219557Sdas && !(word1(&d) & 1) 651112158Sdas#endif 652219557Sdas )) { 653112158Sdas if (!b->x[0] && b->wds <= 1) { 654112158Sdas#ifdef SET_INEXACT 655112158Sdas inexact = 0; 656112158Sdas#endif 657112158Sdas goto accept_dig; 658112158Sdas } 659112158Sdas#ifdef Honor_FLT_ROUNDS 660112158Sdas if (mode > 1) 661182709Sdas switch(Rounding) { 662112158Sdas case 0: goto accept_dig; 663112158Sdas case 2: goto keep_dig; 664112158Sdas } 665112158Sdas#endif /*Honor_FLT_ROUNDS*/ 666112158Sdas if (j1 > 0) { 667112158Sdas b = lshift(b, 1); 668112158Sdas j1 = cmp(b, S); 669219557Sdas#ifdef ROUND_BIASED 670219557Sdas if (j1 >= 0 /*)*/ 671219557Sdas#else 672219557Sdas if ((j1 > 0 || (j1 == 0 && dig & 1)) 673219557Sdas#endif 674112158Sdas && dig++ == '9') 675112158Sdas goto round_9_up; 676112158Sdas } 677112158Sdas accept_dig: 678112158Sdas *s++ = dig; 679112158Sdas goto ret; 680112158Sdas } 681112158Sdas if (j1 > 0) { 682112158Sdas#ifdef Honor_FLT_ROUNDS 683182709Sdas if (!Rounding) 684112158Sdas goto accept_dig; 685112158Sdas#endif 686112158Sdas if (dig == '9') { /* possible if i == 1 */ 687112158Sdas round_9_up: 688112158Sdas *s++ = '9'; 689112158Sdas goto roundoff; 690112158Sdas } 691112158Sdas *s++ = dig + 1; 692112158Sdas goto ret; 693112158Sdas } 694112158Sdas#ifdef Honor_FLT_ROUNDS 695112158Sdas keep_dig: 696112158Sdas#endif 697112158Sdas *s++ = dig; 698112158Sdas if (i == ilim) 699112158Sdas break; 700112158Sdas b = multadd(b, 10, 0); 701112158Sdas if (mlo == mhi) 702112158Sdas mlo = mhi = multadd(mhi, 10, 0); 703112158Sdas else { 704112158Sdas mlo = multadd(mlo, 10, 0); 705112158Sdas mhi = multadd(mhi, 10, 0); 706112158Sdas } 707112158Sdas } 708112158Sdas } 709112158Sdas else 710112158Sdas for(i = 1;; i++) { 711112158Sdas *s++ = dig = quorem(b,S) + '0'; 712112158Sdas if (!b->x[0] && b->wds <= 1) { 713112158Sdas#ifdef SET_INEXACT 714112158Sdas inexact = 0; 715112158Sdas#endif 716112158Sdas goto ret; 717112158Sdas } 718112158Sdas if (i >= ilim) 719112158Sdas break; 720112158Sdas b = multadd(b, 10, 0); 721112158Sdas } 722112158Sdas 723112158Sdas /* Round off last digit */ 724112158Sdas 725112158Sdas#ifdef Honor_FLT_ROUNDS 726182709Sdas switch(Rounding) { 727112158Sdas case 0: goto trimzeros; 728112158Sdas case 2: goto roundoff; 729112158Sdas } 730112158Sdas#endif 731112158Sdas b = lshift(b, 1); 732112158Sdas j = cmp(b, S); 733219557Sdas#ifdef ROUND_BIASED 734219557Sdas if (j >= 0) 735219557Sdas#else 736219557Sdas if (j > 0 || (j == 0 && dig & 1)) 737219557Sdas#endif 738219557Sdas { 739112158Sdas roundoff: 740112158Sdas while(*--s == '9') 741112158Sdas if (s == s0) { 742112158Sdas k++; 743112158Sdas *s++ = '1'; 744112158Sdas goto ret; 745112158Sdas } 746112158Sdas ++*s++; 747112158Sdas } 748112158Sdas else { 749219557Sdas#ifdef Honor_FLT_ROUNDS 750112158Sdas trimzeros: 751219557Sdas#endif 752112158Sdas while(*--s == '0'); 753112158Sdas s++; 754112158Sdas } 755112158Sdas ret: 756112158Sdas Bfree(S); 757112158Sdas if (mhi) { 758112158Sdas if (mlo && mlo != mhi) 759112158Sdas Bfree(mlo); 760112158Sdas Bfree(mhi); 761112158Sdas } 762112158Sdas ret1: 763112158Sdas#ifdef SET_INEXACT 764112158Sdas if (inexact) { 765112158Sdas if (!oldinexact) { 766219557Sdas word0(&d) = Exp_1 + (70 << Exp_shift); 767219557Sdas word1(&d) = 0; 768219557Sdas dval(&d) += 1.; 769112158Sdas } 770112158Sdas } 771112158Sdas else if (!oldinexact) 772112158Sdas clear_inexact(); 773112158Sdas#endif 774112158Sdas Bfree(b); 775112158Sdas *s = 0; 776112158Sdas *decpt = k + 1; 777112158Sdas if (rve) 778112158Sdas *rve = s; 779112158Sdas return s0; 780112158Sdas } 781