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 static Bigint * 35112158Sdas#ifdef KR_headers 36112158Sdasbitstob(bits, nbits, bbits) ULong *bits; int nbits; int *bbits; 37112158Sdas#else 38112158Sdasbitstob(ULong *bits, int nbits, int *bbits) 39112158Sdas#endif 40112158Sdas{ 41112158Sdas int i, k; 42112158Sdas Bigint *b; 43112158Sdas ULong *be, *x, *x0; 44112158Sdas 45112158Sdas i = ULbits; 46112158Sdas k = 0; 47112158Sdas while(i < nbits) { 48112158Sdas i <<= 1; 49112158Sdas k++; 50112158Sdas } 51112158Sdas#ifndef Pack_32 52112158Sdas if (!k) 53112158Sdas k = 1; 54112158Sdas#endif 55112158Sdas b = Balloc(k); 56112158Sdas be = bits + ((nbits - 1) >> kshift); 57112158Sdas x = x0 = b->x; 58112158Sdas do { 59112158Sdas *x++ = *bits & ALL_ON; 60112158Sdas#ifdef Pack_16 61112158Sdas *x++ = (*bits >> 16) & ALL_ON; 62112158Sdas#endif 63112158Sdas } while(++bits <= be); 64112158Sdas i = x - x0; 65112158Sdas while(!x0[--i]) 66112158Sdas if (!i) { 67112158Sdas b->wds = 0; 68112158Sdas *bbits = 0; 69112158Sdas goto ret; 70112158Sdas } 71112158Sdas b->wds = i + 1; 72112158Sdas *bbits = i*ULbits + 32 - hi0bits(b->x[i]); 73112158Sdas ret: 74112158Sdas return b; 75112158Sdas } 76112158Sdas 77112158Sdas/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. 78112158Sdas * 79112158Sdas * Inspired by "How to Print Floating-Point Numbers Accurately" by 80165743Sdas * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126]. 81112158Sdas * 82112158Sdas * Modifications: 83112158Sdas * 1. Rather than iterating, we use a simple numeric overestimate 84112158Sdas * to determine k = floor(log10(d)). We scale relevant 85112158Sdas * quantities using O(log2(k)) rather than O(k) multiplications. 86112158Sdas * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't 87112158Sdas * try to generate digits strictly left to right. Instead, we 88112158Sdas * compute with fewer bits and propagate the carry if necessary 89112158Sdas * when rounding the final digit up. This is often faster. 90112158Sdas * 3. Under the assumption that input will be rounded nearest, 91112158Sdas * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. 92112158Sdas * That is, we allow equality in stopping tests when the 93112158Sdas * round-nearest rule will give the same floating-point value 94112158Sdas * as would satisfaction of the stopping test with strict 95112158Sdas * inequality. 96112158Sdas * 4. We remove common factors of powers of 2 from relevant 97112158Sdas * quantities. 98112158Sdas * 5. When converting floating-point integers less than 1e16, 99112158Sdas * we use floating-point arithmetic rather than resorting 100112158Sdas * to multiple-precision integers. 101112158Sdas * 6. When asked to produce fewer than 15 digits, we first try 102112158Sdas * to get by with floating-point arithmetic; we resort to 103112158Sdas * multiple-precision integer arithmetic only if we cannot 104112158Sdas * guarantee that the floating-point calculation has given 105112158Sdas * the correctly rounded result. For k requested digits and 106112158Sdas * "uniformly" distributed input, the probability is 107112158Sdas * something like 10^(k-15) that we must resort to the Long 108112158Sdas * calculation. 109112158Sdas */ 110112158Sdas 111112158Sdas char * 112112158Sdasgdtoa 113112158Sdas#ifdef KR_headers 114112158Sdas (fpi, be, bits, kindp, mode, ndigits, decpt, rve) 115112158Sdas FPI *fpi; int be; ULong *bits; 116112158Sdas int *kindp, mode, ndigits, *decpt; char **rve; 117112158Sdas#else 118112158Sdas (FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, int *decpt, char **rve) 119112158Sdas#endif 120112158Sdas{ 121112158Sdas /* Arguments ndigits and decpt are similar to the second and third 122112158Sdas arguments of ecvt and fcvt; trailing zeros are suppressed from 123112158Sdas the returned string. If not null, *rve is set to point 124112158Sdas to the end of the return value. If d is +-Infinity or NaN, 125112158Sdas then *decpt is set to 9999. 126219557Sdas be = exponent: value = (integer represented by bits) * (2 to the power of be). 127112158Sdas 128112158Sdas mode: 129112158Sdas 0 ==> shortest string that yields d when read in 130112158Sdas and rounded to nearest. 131112158Sdas 1 ==> like 0, but with Steele & White stopping rule; 132112158Sdas e.g. with IEEE P754 arithmetic , mode 0 gives 133112158Sdas 1e23 whereas mode 1 gives 9.999999999999999e22. 134112158Sdas 2 ==> max(1,ndigits) significant digits. This gives a 135112158Sdas return value similar to that of ecvt, except 136112158Sdas that trailing zeros are suppressed. 137112158Sdas 3 ==> through ndigits past the decimal point. This 138112158Sdas gives a return value similar to that from fcvt, 139112158Sdas except that trailing zeros are suppressed, and 140112158Sdas ndigits can be negative. 141112158Sdas 4-9 should give the same return values as 2-3, i.e., 142112158Sdas 4 <= mode <= 9 ==> same return as mode 143112158Sdas 2 + (mode & 1). These modes are mainly for 144112158Sdas debugging; often they run slower but sometimes 145112158Sdas faster than modes 2-3. 146112158Sdas 4,5,8,9 ==> left-to-right digit generation. 147112158Sdas 6-9 ==> don't try fast floating-point estimate 148112158Sdas (if applicable). 149112158Sdas 150112158Sdas Values of mode other than 0-9 are treated as mode 0. 151112158Sdas 152112158Sdas Sufficient space is allocated to the return value 153112158Sdas to hold the suppressed trailing zeros. 154112158Sdas */ 155112158Sdas 156112158Sdas int bbits, b2, b5, be0, dig, i, ieps, ilim, ilim0, ilim1, inex; 157112158Sdas int j, j1, k, k0, k_check, kind, leftright, m2, m5, nbits; 158112158Sdas int rdir, s2, s5, spec_case, try_quick; 159112158Sdas Long L; 160112158Sdas Bigint *b, *b1, *delta, *mlo, *mhi, *mhi1, *S; 161219557Sdas double d2, ds; 162112158Sdas char *s, *s0; 163219557Sdas U d, eps; 164112158Sdas 165112158Sdas#ifndef MULTIPLE_THREADS 166112158Sdas if (dtoa_result) { 167112158Sdas freedtoa(dtoa_result); 168112158Sdas dtoa_result = 0; 169112158Sdas } 170112158Sdas#endif 171112158Sdas inex = 0; 172112158Sdas kind = *kindp &= ~STRTOG_Inexact; 173112158Sdas switch(kind & STRTOG_Retmask) { 174112158Sdas case STRTOG_Zero: 175112158Sdas goto ret_zero; 176112158Sdas case STRTOG_Normal: 177112158Sdas case STRTOG_Denormal: 178112158Sdas break; 179112158Sdas case STRTOG_Infinite: 180112158Sdas *decpt = -32768; 181112158Sdas return nrv_alloc("Infinity", rve, 8); 182112158Sdas case STRTOG_NaN: 183112158Sdas *decpt = -32768; 184112158Sdas return nrv_alloc("NaN", rve, 3); 185112158Sdas default: 186112158Sdas return 0; 187112158Sdas } 188112158Sdas b = bitstob(bits, nbits = fpi->nbits, &bbits); 189112158Sdas be0 = be; 190112158Sdas if ( (i = trailz(b)) !=0) { 191112158Sdas rshift(b, i); 192112158Sdas be += i; 193112158Sdas bbits -= i; 194112158Sdas } 195112158Sdas if (!b->wds) { 196112158Sdas Bfree(b); 197112158Sdas ret_zero: 198112158Sdas *decpt = 1; 199112158Sdas return nrv_alloc("0", rve, 1); 200112158Sdas } 201112158Sdas 202219557Sdas dval(&d) = b2d(b, &i); 203112158Sdas i = be + bbits - 1; 204219557Sdas word0(&d) &= Frac_mask1; 205219557Sdas word0(&d) |= Exp_11; 206112158Sdas#ifdef IBM 207219557Sdas if ( (j = 11 - hi0bits(word0(&d) & Frac_mask)) !=0) 208219557Sdas dval(&d) /= 1 << j; 209112158Sdas#endif 210112158Sdas 211112158Sdas /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 212112158Sdas * log10(x) = log(x) / log(10) 213112158Sdas * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) 214219557Sdas * log10(&d) = (i-Bias)*log(2)/log(10) + log10(d2) 215112158Sdas * 216219557Sdas * This suggests computing an approximation k to log10(&d) by 217112158Sdas * 218112158Sdas * k = (i - Bias)*0.301029995663981 219112158Sdas * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); 220112158Sdas * 221112158Sdas * We want k to be too large rather than too small. 222112158Sdas * The error in the first-order Taylor series approximation 223112158Sdas * is in our favor, so we just round up the constant enough 224112158Sdas * to compensate for any error in the multiplication of 225112158Sdas * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, 226112158Sdas * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, 227112158Sdas * adding 1e-13 to the constant term more than suffices. 228112158Sdas * Hence we adjust the constant term to 0.1760912590558. 229112158Sdas * (We could get a more accurate k by invoking log10, 230112158Sdas * but this is probably not worthwhile.) 231112158Sdas */ 232112158Sdas#ifdef IBM 233112158Sdas i <<= 2; 234112158Sdas i += j; 235112158Sdas#endif 236219557Sdas ds = (dval(&d)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981; 237112158Sdas 238112158Sdas /* correct assumption about exponent range */ 239112158Sdas if ((j = i) < 0) 240112158Sdas j = -j; 241112158Sdas if ((j -= 1077) > 0) 242112158Sdas ds += j * 7e-17; 243112158Sdas 244112158Sdas k = (int)ds; 245112158Sdas if (ds < 0. && ds != k) 246112158Sdas k--; /* want k = floor(ds) */ 247112158Sdas k_check = 1; 248112158Sdas#ifdef IBM 249112158Sdas j = be + bbits - 1; 250112158Sdas if ( (j1 = j & 3) !=0) 251219557Sdas dval(&d) *= 1 << j1; 252219557Sdas word0(&d) += j << Exp_shift - 2 & Exp_mask; 253112158Sdas#else 254219557Sdas word0(&d) += (be + bbits - 1) << Exp_shift; 255112158Sdas#endif 256112158Sdas if (k >= 0 && k <= Ten_pmax) { 257219557Sdas if (dval(&d) < tens[k]) 258112158Sdas k--; 259112158Sdas k_check = 0; 260112158Sdas } 261112158Sdas j = bbits - i - 1; 262112158Sdas if (j >= 0) { 263112158Sdas b2 = 0; 264112158Sdas s2 = j; 265112158Sdas } 266112158Sdas else { 267112158Sdas b2 = -j; 268112158Sdas s2 = 0; 269112158Sdas } 270112158Sdas if (k >= 0) { 271112158Sdas b5 = 0; 272112158Sdas s5 = k; 273112158Sdas s2 += k; 274112158Sdas } 275112158Sdas else { 276112158Sdas b2 -= k; 277112158Sdas b5 = -k; 278112158Sdas s5 = 0; 279112158Sdas } 280112158Sdas if (mode < 0 || mode > 9) 281112158Sdas mode = 0; 282112158Sdas try_quick = 1; 283112158Sdas if (mode > 5) { 284112158Sdas mode -= 4; 285112158Sdas try_quick = 0; 286112158Sdas } 287219557Sdas else if (i >= -4 - Emin || i < Emin) 288219557Sdas try_quick = 0; 289112158Sdas leftright = 1; 290219557Sdas ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */ 291219557Sdas /* silence erroneous "gcc -Wall" warning. */ 292112158Sdas switch(mode) { 293112158Sdas case 0: 294112158Sdas case 1: 295112158Sdas i = (int)(nbits * .30103) + 3; 296112158Sdas ndigits = 0; 297112158Sdas break; 298112158Sdas case 2: 299112158Sdas leftright = 0; 300112158Sdas /* no break */ 301112158Sdas case 4: 302112158Sdas if (ndigits <= 0) 303112158Sdas ndigits = 1; 304112158Sdas ilim = ilim1 = i = ndigits; 305112158Sdas break; 306112158Sdas case 3: 307112158Sdas leftright = 0; 308112158Sdas /* no break */ 309112158Sdas case 5: 310112158Sdas i = ndigits + k + 1; 311112158Sdas ilim = i; 312112158Sdas ilim1 = i - 1; 313112158Sdas if (i <= 0) 314112158Sdas i = 1; 315112158Sdas } 316112158Sdas s = s0 = rv_alloc(i); 317112158Sdas 318112158Sdas if ( (rdir = fpi->rounding - 1) !=0) { 319112158Sdas if (rdir < 0) 320112158Sdas rdir = 2; 321112158Sdas if (kind & STRTOG_Neg) 322112158Sdas rdir = 3 - rdir; 323112158Sdas } 324112158Sdas 325112158Sdas /* Now rdir = 0 ==> round near, 1 ==> round up, 2 ==> round down. */ 326112158Sdas 327112158Sdas if (ilim >= 0 && ilim <= Quick_max && try_quick && !rdir 328112158Sdas#ifndef IMPRECISE_INEXACT 329112158Sdas && k == 0 330112158Sdas#endif 331112158Sdas ) { 332112158Sdas 333112158Sdas /* Try to get by with floating-point arithmetic. */ 334112158Sdas 335112158Sdas i = 0; 336219557Sdas d2 = dval(&d); 337112158Sdas#ifdef IBM 338219557Sdas if ( (j = 11 - hi0bits(word0(&d) & Frac_mask)) !=0) 339219557Sdas dval(&d) /= 1 << j; 340112158Sdas#endif 341112158Sdas k0 = k; 342112158Sdas ilim0 = ilim; 343112158Sdas ieps = 2; /* conservative */ 344112158Sdas if (k > 0) { 345112158Sdas ds = tens[k&0xf]; 346112158Sdas j = k >> 4; 347112158Sdas if (j & Bletch) { 348112158Sdas /* prevent overflows */ 349112158Sdas j &= Bletch - 1; 350219557Sdas dval(&d) /= bigtens[n_bigtens-1]; 351112158Sdas ieps++; 352112158Sdas } 353112158Sdas for(; j; j >>= 1, i++) 354112158Sdas if (j & 1) { 355112158Sdas ieps++; 356112158Sdas ds *= bigtens[i]; 357112158Sdas } 358112158Sdas } 359112158Sdas else { 360112158Sdas ds = 1.; 361112158Sdas if ( (j1 = -k) !=0) { 362219557Sdas dval(&d) *= tens[j1 & 0xf]; 363112158Sdas for(j = j1 >> 4; j; j >>= 1, i++) 364112158Sdas if (j & 1) { 365112158Sdas ieps++; 366219557Sdas dval(&d) *= bigtens[i]; 367112158Sdas } 368112158Sdas } 369112158Sdas } 370219557Sdas if (k_check && dval(&d) < 1. && ilim > 0) { 371112158Sdas if (ilim1 <= 0) 372112158Sdas goto fast_failed; 373112158Sdas ilim = ilim1; 374112158Sdas k--; 375219557Sdas dval(&d) *= 10.; 376112158Sdas ieps++; 377112158Sdas } 378219557Sdas dval(&eps) = ieps*dval(&d) + 7.; 379219557Sdas word0(&eps) -= (P-1)*Exp_msk1; 380112158Sdas if (ilim == 0) { 381112158Sdas S = mhi = 0; 382219557Sdas dval(&d) -= 5.; 383219557Sdas if (dval(&d) > dval(&eps)) 384112158Sdas goto one_digit; 385219557Sdas if (dval(&d) < -dval(&eps)) 386112158Sdas goto no_digits; 387112158Sdas goto fast_failed; 388112158Sdas } 389112158Sdas#ifndef No_leftright 390112158Sdas if (leftright) { 391112158Sdas /* Use Steele & White method of only 392112158Sdas * generating digits needed. 393112158Sdas */ 394219557Sdas dval(&eps) = ds*0.5/tens[ilim-1] - dval(&eps); 395112158Sdas for(i = 0;;) { 396219557Sdas L = (Long)(dval(&d)/ds); 397219557Sdas dval(&d) -= L*ds; 398112158Sdas *s++ = '0' + (int)L; 399219557Sdas if (dval(&d) < dval(&eps)) { 400219557Sdas if (dval(&d)) 401112158Sdas inex = STRTOG_Inexlo; 402112158Sdas goto ret1; 403112158Sdas } 404219557Sdas if (ds - dval(&d) < dval(&eps)) 405112158Sdas goto bump_up; 406112158Sdas if (++i >= ilim) 407112158Sdas break; 408219557Sdas dval(&eps) *= 10.; 409219557Sdas dval(&d) *= 10.; 410112158Sdas } 411112158Sdas } 412112158Sdas else { 413112158Sdas#endif 414112158Sdas /* Generate ilim digits, then fix them up. */ 415219557Sdas dval(&eps) *= tens[ilim-1]; 416219557Sdas for(i = 1;; i++, dval(&d) *= 10.) { 417219557Sdas if ( (L = (Long)(dval(&d)/ds)) !=0) 418219557Sdas dval(&d) -= L*ds; 419112158Sdas *s++ = '0' + (int)L; 420112158Sdas if (i == ilim) { 421112158Sdas ds *= 0.5; 422219557Sdas if (dval(&d) > ds + dval(&eps)) 423112158Sdas goto bump_up; 424219557Sdas else if (dval(&d) < ds - dval(&eps)) { 425219557Sdas if (dval(&d)) 426112158Sdas inex = STRTOG_Inexlo; 427187808Sdas goto clear_trailing0; 428112158Sdas } 429112158Sdas break; 430112158Sdas } 431112158Sdas } 432112158Sdas#ifndef No_leftright 433112158Sdas } 434112158Sdas#endif 435112158Sdas fast_failed: 436112158Sdas s = s0; 437219557Sdas dval(&d) = d2; 438112158Sdas k = k0; 439112158Sdas ilim = ilim0; 440112158Sdas } 441112158Sdas 442112158Sdas /* Do we have a "small" integer? */ 443112158Sdas 444112158Sdas if (be >= 0 && k <= Int_max) { 445112158Sdas /* Yes. */ 446112158Sdas ds = tens[k]; 447112158Sdas if (ndigits < 0 && ilim <= 0) { 448112158Sdas S = mhi = 0; 449219557Sdas if (ilim < 0 || dval(&d) <= 5*ds) 450112158Sdas goto no_digits; 451112158Sdas goto one_digit; 452112158Sdas } 453219557Sdas for(i = 1;; i++, dval(&d) *= 10.) { 454219557Sdas L = dval(&d) / ds; 455219557Sdas dval(&d) -= L*ds; 456112158Sdas#ifdef Check_FLT_ROUNDS 457112158Sdas /* If FLT_ROUNDS == 2, L will usually be high by 1 */ 458219557Sdas if (dval(&d) < 0) { 459112158Sdas L--; 460219557Sdas dval(&d) += ds; 461112158Sdas } 462112158Sdas#endif 463112158Sdas *s++ = '0' + (int)L; 464219557Sdas if (dval(&d) == 0.) 465112158Sdas break; 466112158Sdas if (i == ilim) { 467112158Sdas if (rdir) { 468112158Sdas if (rdir == 1) 469112158Sdas goto bump_up; 470112158Sdas inex = STRTOG_Inexlo; 471112158Sdas goto ret1; 472112158Sdas } 473219557Sdas dval(&d) += dval(&d); 474219557Sdas#ifdef ROUND_BIASED 475219557Sdas if (dval(&d) >= ds) 476219557Sdas#else 477219557Sdas if (dval(&d) > ds || (dval(&d) == ds && L & 1)) 478219557Sdas#endif 479219557Sdas { 480112158Sdas bump_up: 481112158Sdas inex = STRTOG_Inexhi; 482112158Sdas while(*--s == '9') 483112158Sdas if (s == s0) { 484112158Sdas k++; 485112158Sdas *s = '0'; 486112158Sdas break; 487112158Sdas } 488112158Sdas ++*s++; 489112158Sdas } 490187808Sdas else { 491112158Sdas inex = STRTOG_Inexlo; 492187808Sdas clear_trailing0: 493187808Sdas while(*--s == '0'){} 494187808Sdas ++s; 495187808Sdas } 496112158Sdas break; 497112158Sdas } 498112158Sdas } 499112158Sdas goto ret1; 500112158Sdas } 501112158Sdas 502112158Sdas m2 = b2; 503112158Sdas m5 = b5; 504112158Sdas mhi = mlo = 0; 505112158Sdas if (leftright) { 506219557Sdas i = nbits - bbits; 507219557Sdas if (be - i++ < fpi->emin && mode != 3 && mode != 5) { 508219557Sdas /* denormal */ 509219557Sdas i = be - fpi->emin + 1; 510219557Sdas if (mode >= 2 && ilim > 0 && ilim < i) 511219557Sdas goto small_ilim; 512112158Sdas } 513219557Sdas else if (mode >= 2) { 514219557Sdas small_ilim: 515112158Sdas j = ilim - 1; 516112158Sdas if (m5 >= j) 517112158Sdas m5 -= j; 518112158Sdas else { 519112158Sdas s5 += j -= m5; 520112158Sdas b5 += j; 521112158Sdas m5 = 0; 522112158Sdas } 523112158Sdas if ((i = ilim) < 0) { 524112158Sdas m2 -= i; 525112158Sdas i = 0; 526112158Sdas } 527112158Sdas } 528112158Sdas b2 += i; 529112158Sdas s2 += i; 530112158Sdas mhi = i2b(1); 531112158Sdas } 532112158Sdas if (m2 > 0 && s2 > 0) { 533112158Sdas i = m2 < s2 ? m2 : s2; 534112158Sdas b2 -= i; 535112158Sdas m2 -= i; 536112158Sdas s2 -= i; 537112158Sdas } 538112158Sdas if (b5 > 0) { 539112158Sdas if (leftright) { 540112158Sdas if (m5 > 0) { 541112158Sdas mhi = pow5mult(mhi, m5); 542112158Sdas b1 = mult(mhi, b); 543112158Sdas Bfree(b); 544112158Sdas b = b1; 545112158Sdas } 546112158Sdas if ( (j = b5 - m5) !=0) 547112158Sdas b = pow5mult(b, j); 548112158Sdas } 549112158Sdas else 550112158Sdas b = pow5mult(b, b5); 551112158Sdas } 552112158Sdas S = i2b(1); 553112158Sdas if (s5 > 0) 554112158Sdas S = pow5mult(S, s5); 555112158Sdas 556112158Sdas /* Check for special case that d is a normalized power of 2. */ 557112158Sdas 558112158Sdas spec_case = 0; 559112158Sdas if (mode < 2) { 560112158Sdas if (bbits == 1 && be0 > fpi->emin + 1) { 561112158Sdas /* The special case */ 562112158Sdas b2++; 563112158Sdas s2++; 564112158Sdas spec_case = 1; 565112158Sdas } 566112158Sdas } 567112158Sdas 568112158Sdas /* Arrange for convenient computation of quotients: 569112158Sdas * shift left if necessary so divisor has 4 leading 0 bits. 570112158Sdas * 571112158Sdas * Perhaps we should just compute leading 28 bits of S once 572112158Sdas * and for all and pass them and a shift to quorem, so it 573112158Sdas * can do shifts and ors to compute the numerator for q. 574112158Sdas */ 575219557Sdas i = ((s5 ? hi0bits(S->x[S->wds-1]) : ULbits - 1) - s2 - 4) & kmask; 576219557Sdas m2 += i; 577219557Sdas if ((b2 += i) > 0) 578112158Sdas b = lshift(b, b2); 579219557Sdas if ((s2 += i) > 0) 580112158Sdas S = lshift(S, s2); 581112158Sdas if (k_check) { 582112158Sdas if (cmp(b,S) < 0) { 583112158Sdas k--; 584112158Sdas b = multadd(b, 10, 0); /* we botched the k estimate */ 585112158Sdas if (leftright) 586112158Sdas mhi = multadd(mhi, 10, 0); 587112158Sdas ilim = ilim1; 588112158Sdas } 589112158Sdas } 590112158Sdas if (ilim <= 0 && mode > 2) { 591112158Sdas if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) { 592112158Sdas /* no digits, fcvt style */ 593112158Sdas no_digits: 594112158Sdas k = -1 - ndigits; 595112158Sdas inex = STRTOG_Inexlo; 596112158Sdas goto ret; 597112158Sdas } 598112158Sdas one_digit: 599112158Sdas inex = STRTOG_Inexhi; 600112158Sdas *s++ = '1'; 601112158Sdas k++; 602112158Sdas goto ret; 603112158Sdas } 604112158Sdas if (leftright) { 605112158Sdas if (m2 > 0) 606112158Sdas mhi = lshift(mhi, m2); 607112158Sdas 608112158Sdas /* Compute mlo -- check for special case 609112158Sdas * that d is a normalized power of 2. 610112158Sdas */ 611112158Sdas 612112158Sdas mlo = mhi; 613112158Sdas if (spec_case) { 614112158Sdas mhi = Balloc(mhi->k); 615112158Sdas Bcopy(mhi, mlo); 616112158Sdas mhi = lshift(mhi, 1); 617112158Sdas } 618112158Sdas 619112158Sdas for(i = 1;;i++) { 620112158Sdas dig = quorem(b,S) + '0'; 621112158Sdas /* Do we yet have the shortest decimal string 622112158Sdas * that will round to d? 623112158Sdas */ 624112158Sdas j = cmp(b, mlo); 625112158Sdas delta = diff(S, mhi); 626112158Sdas j1 = delta->sign ? 1 : cmp(b, delta); 627112158Sdas Bfree(delta); 628112158Sdas#ifndef ROUND_BIASED 629112158Sdas if (j1 == 0 && !mode && !(bits[0] & 1) && !rdir) { 630112158Sdas if (dig == '9') 631112158Sdas goto round_9_up; 632112158Sdas if (j <= 0) { 633112158Sdas if (b->wds > 1 || b->x[0]) 634112158Sdas inex = STRTOG_Inexlo; 635112158Sdas } 636112158Sdas else { 637112158Sdas dig++; 638112158Sdas inex = STRTOG_Inexhi; 639112158Sdas } 640112158Sdas *s++ = dig; 641112158Sdas goto ret; 642112158Sdas } 643112158Sdas#endif 644219557Sdas if (j < 0 || (j == 0 && !mode 645112158Sdas#ifndef ROUND_BIASED 646112158Sdas && !(bits[0] & 1) 647112158Sdas#endif 648219557Sdas )) { 649112158Sdas if (rdir && (b->wds > 1 || b->x[0])) { 650112158Sdas if (rdir == 2) { 651112158Sdas inex = STRTOG_Inexlo; 652112158Sdas goto accept; 653112158Sdas } 654112158Sdas while (cmp(S,mhi) > 0) { 655112158Sdas *s++ = dig; 656112158Sdas mhi1 = multadd(mhi, 10, 0); 657112158Sdas if (mlo == mhi) 658112158Sdas mlo = mhi1; 659112158Sdas mhi = mhi1; 660112158Sdas b = multadd(b, 10, 0); 661112158Sdas dig = quorem(b,S) + '0'; 662112158Sdas } 663112158Sdas if (dig++ == '9') 664112158Sdas goto round_9_up; 665112158Sdas inex = STRTOG_Inexhi; 666112158Sdas goto accept; 667112158Sdas } 668112158Sdas if (j1 > 0) { 669112158Sdas b = lshift(b, 1); 670112158Sdas j1 = cmp(b, S); 671219557Sdas#ifdef ROUND_BIASED 672219557Sdas if (j1 >= 0 /*)*/ 673219557Sdas#else 674219557Sdas if ((j1 > 0 || (j1 == 0 && dig & 1)) 675219557Sdas#endif 676112158Sdas && dig++ == '9') 677112158Sdas goto round_9_up; 678112158Sdas inex = STRTOG_Inexhi; 679112158Sdas } 680112158Sdas if (b->wds > 1 || b->x[0]) 681112158Sdas inex = STRTOG_Inexlo; 682112158Sdas accept: 683112158Sdas *s++ = dig; 684112158Sdas goto ret; 685112158Sdas } 686112158Sdas if (j1 > 0 && rdir != 2) { 687112158Sdas if (dig == '9') { /* possible if i == 1 */ 688112158Sdas round_9_up: 689112158Sdas *s++ = '9'; 690112158Sdas inex = STRTOG_Inexhi; 691112158Sdas goto roundoff; 692112158Sdas } 693112158Sdas inex = STRTOG_Inexhi; 694112158Sdas *s++ = dig + 1; 695112158Sdas goto ret; 696112158Sdas } 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 (i >= ilim) 713112158Sdas break; 714112158Sdas b = multadd(b, 10, 0); 715112158Sdas } 716112158Sdas 717112158Sdas /* Round off last digit */ 718112158Sdas 719112158Sdas if (rdir) { 720219557Sdas if (rdir == 2 || (b->wds <= 1 && !b->x[0])) 721112158Sdas goto chopzeros; 722112158Sdas goto roundoff; 723112158Sdas } 724112158Sdas b = lshift(b, 1); 725112158Sdas j = cmp(b, S); 726219557Sdas#ifdef ROUND_BIASED 727219557Sdas if (j >= 0) 728219557Sdas#else 729219557Sdas if (j > 0 || (j == 0 && dig & 1)) 730219557Sdas#endif 731219557Sdas { 732112158Sdas roundoff: 733112158Sdas inex = STRTOG_Inexhi; 734112158Sdas while(*--s == '9') 735112158Sdas if (s == s0) { 736112158Sdas k++; 737112158Sdas *s++ = '1'; 738112158Sdas goto ret; 739112158Sdas } 740112158Sdas ++*s++; 741112158Sdas } 742112158Sdas else { 743112158Sdas chopzeros: 744112158Sdas if (b->wds > 1 || b->x[0]) 745112158Sdas inex = STRTOG_Inexlo; 746112158Sdas while(*--s == '0'){} 747187808Sdas ++s; 748112158Sdas } 749112158Sdas ret: 750112158Sdas Bfree(S); 751112158Sdas if (mhi) { 752112158Sdas if (mlo && mlo != mhi) 753112158Sdas Bfree(mlo); 754112158Sdas Bfree(mhi); 755112158Sdas } 756112158Sdas ret1: 757112158Sdas Bfree(b); 758112158Sdas *s = 0; 759112158Sdas *decpt = k + 1; 760112158Sdas if (rve) 761112158Sdas *rve = s; 762112158Sdas *kindp |= inex; 763112158Sdas return s0; 764112158Sdas } 765