lib/libc/gdtoa/gdtoa.c

   1 /* $NetBSD: gdtoa.c,v 1.4 2008/03/21 23:13:48 christos Exp $ */
   2
   3 /****************************************************************
   4
   5 The author of this software is David M. Gay.
   6
   7 Copyright (C) 1998, 1999 by Lucent Technologies
   8 All Rights Reserved
   9
  10 Permission to use, copy, modify, and distribute this software and
  11 its documentation for any purpose and without fee is hereby
  12 granted, provided that the above copyright notice appear in all
  13 copies and that both that the copyright notice and this
  14 permission notice and warranty disclaimer appear in supporting
  15 documentation, and that the name of Lucent or any of its entities
  16 not be used in advertising or publicity pertaining to
  17 distribution of the software without specific, written prior
  18 permission.
  19
  20 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
  21 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
  22 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
  23 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  24 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
  25 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
  26 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
  27 THIS SOFTWARE.
  28
  29 ****************************************************************/
  30
  31 /* Please send bug reports to David M. Gay (dmg at acm dot org,
  32  * with " at " changed at "@" and " dot " changed to ".").      */
  33
  34 #include "gdtoaimp.h"
  35
  36  static Bigint *
  37 #ifdef KR_headers
  38 bitstob(bits, nbits, bbits) ULong *bits; int nbits; int *bbits;
  39 #else
  40 bitstob(ULong *bits, int nbits, int *bbits)
  41 #endif
  42 {
  43         int i, k;
  44         Bigint *b;
  45         ULong *be, *x, *x0;
  46
  47         i = ULbits;
  48         k = 0;
  49         while(i < nbits) {
  50                 i <<= 1;
  51                 k++;
  52                 }
  53 #ifndef Pack_32
  54         if (!k)
  55                 k = 1;
  56 #endif
  57         b = Balloc(k);
  58         if (b == NULL)
  59                 return NULL;
  60         be = bits + (((unsigned int)nbits - 1) >> kshift);
  61         x = x0 = b->x;
  62         do {
  63                 *x++ = *bits & ALL_ON;
  64 #ifdef Pack_16
  65                 *x++ = (*bits >> 16) & ALL_ON;
  66 #endif
  67                 } while(++bits <= be);
  68         i = x - x0;
  69         while(!x0[--i])
  70                 if (!i) {
  71                         b->wds = 0;
  72                         *bbits = 0;
  73                         goto ret;
  74                         }
  75         b->wds = i + 1;
  76         *bbits = i*ULbits + 32 - hi0bits(b->x[i]);
  77  ret:
  78         return b;
  79         }
  80
  81 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
  82  *
  83  * Inspired by "How to Print Floating-Point Numbers Accurately" by
  84  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
  85  *
  86  * Modifications:
  87  *      1. Rather than iterating, we use a simple numeric overestimate
  88  *         to determine k = floor(log10(d)).  We scale relevant
  89  *         quantities using O(log2(k)) rather than O(k) multiplications.
  90  *      2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
  91  *         try to generate digits strictly left to right.  Instead, we
  92  *         compute with fewer bits and propagate the carry if necessary
  93  *         when rounding the final digit up.  This is often faster.
  94  *      3. Under the assumption that input will be rounded nearest,
  95  *         mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
  96  *         That is, we allow equality in stopping tests when the
  97  *         round-nearest rule will give the same floating-point value
  98  *         as would satisfaction of the stopping test with strict
  99  *         inequality.
 100  *      4. We remove common factors of powers of 2 from relevant
 101  *         quantities.
 102  *      5. When converting floating-point integers less than 1e16,
 103  *         we use floating-point arithmetic rather than resorting
 104  *         to multiple-precision integers.
 105  *      6. When asked to produce fewer than 15 digits, we first try
 106  *         to get by with floating-point arithmetic; we resort to
 107  *         multiple-precision integer arithmetic only if we cannot
 108  *         guarantee that the floating-point calculation has given
 109  *         the correctly rounded result.  For k requested digits and
 110  *         "uniformly" distributed input, the probability is
 111  *         something like 10^(k-15) that we must resort to the Long
 112  *         calculation.
 113  */
 114
 115  char *
 116 gdtoa
 117 #ifdef KR_headers
 118         (fpi, be, bits, kindp, mode, ndigits, decpt, rve)
 119         FPI *fpi; int be; ULong *bits;
 120         int *kindp, mode, ndigits, *decpt; char **rve;
 121 #else
 122         (FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, int *decpt, char **rve)
 123 #endif
 124 {
 125  /*     Arguments ndigits and decpt are similar to the second and third
 126         arguments of ecvt and fcvt; trailing zeros are suppressed from
 127         the returned string.  If not null, *rve is set to point
 128         to the end of the return value.  If d is +-Infinity or NaN,
 129         then *decpt is set to 9999.
 130
 131         mode:
 132                 0 ==> shortest string that yields d when read in
 133                         and rounded to nearest.
 134                 1 ==> like 0, but with Steele & White stopping rule;
 135                         e.g. with IEEE P754 arithmetic , mode 0 gives
 136                         1e23 whereas mode 1 gives 9.999999999999999e22.
 137                 2 ==> max(1,ndigits) significant digits.  This gives a
 138                         return value similar to that of ecvt, except
 139                         that trailing zeros are suppressed.
 140                 3 ==> through ndigits past the decimal point.  This
 141                         gives a return value similar to that from fcvt,
 142                         except that trailing zeros are suppressed, and
 143                         ndigits can be negative.
 144                 4-9 should give the same return values as 2-3, i.e.,
 145                         4 <= mode <= 9 ==> same return as mode
 146                         2 + (mode & 1).  These modes are mainly for
 147                         debugging; often they run slower but sometimes
 148                         faster than modes 2-3.
 149                 4,5,8,9 ==> left-to-right digit generation.
 150                 6-9 ==> don't try fast floating-point estimate
 151                         (if applicable).
 152
 153                 Values of mode other than 0-9 are treated as mode 0.
 154
 155                 Sufficient space is allocated to the return value
 156                 to hold the suppressed trailing zeros.
 157         */
 158
 159         int bbits, b2, b5, be0, dig, i, ieps, ilim = 0, ilim0, ilim1 = 0, inex;
 160         int j, jj1, k, k0, k_check, kind, leftright, m2, m5, nbits;
 161         int rdir, s2, s5, spec_case, try_quick;
 162         Long L;
 163         Bigint *b, *b1, *delta, *mlo, *mhi, *mhi1, *S;
 164         double d, d2, ds, eps;
 165         char *s, *s0;
 166
 167 #ifndef MULTIPLE_THREADS
 168         if (dtoa_result) {
 169                 freedtoa(dtoa_result);
 170                 dtoa_result = 0;
 171                 }
 172 #endif
 173         inex = 0;
 174         if (*kindp & STRTOG_NoMemory)
 175                 return NULL;
 176         kind = *kindp &= ~STRTOG_Inexact;
 177         switch(kind & STRTOG_Retmask) {
 178           case STRTOG_Zero:
 179                 goto ret_zero;
 180           case STRTOG_Normal:
 181           case STRTOG_Denormal:
 182                 break;
 183           case STRTOG_Infinite:
 184                 *decpt = -32768;
 185                 return nrv_alloc("Infinity", rve, 8);
 186           case STRTOG_NaN:
 187                 *decpt = -32768;
 188                 return nrv_alloc("NaN", rve, 3);
 189           default:
 190                 return 0;
 191           }
 192         b = bitstob(bits, nbits = fpi->nbits, &bbits);
 193         if (b == NULL)
 194                 return NULL;
 195         be0 = be;
 196         if ( (i = trailz(b)) !=0) {
 197                 rshift(b, i);
 198                 be += i;
 199                 bbits -= i;
 200                 }
 201         if (!b->wds) {
 202                 Bfree(b);
 203  ret_zero:
 204                 *decpt = 1;
 205                 return nrv_alloc("0", rve, 1);
 206                 }
 207
 208         dval(d) = b2d(b, &i);
 209         i = be + bbits - 1;
 210         word0(d) &= Frac_mask1;
 211         word0(d) |= Exp_11;
 212 #ifdef IBM
 213         if ( (j = 11 - hi0bits(word0(d) & Frac_mask)) !=0)
 214                 dval(d) /= 1 << j;
 215 #endif
 216
 217         /* log(x)       ~=~ log(1.5) + (x-1.5)/1.5
 218          * log10(x)      =  log(x) / log(10)
 219          *              ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
 220          * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
 221          *
 222          * This suggests computing an approximation k to log10(d) by
 223          *
 224          * k = (i - Bias)*0.301029995663981
 225          *      + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
 226          *
 227          * We want k to be too large rather than too small.
 228          * The error in the first-order Taylor series approximation
 229          * is in our favor, so we just round up the constant enough
 230          * to compensate for any error in the multiplication of
 231          * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
 232          * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
 233          * adding 1e-13 to the constant term more than suffices.
 234          * Hence we adjust the constant term to 0.1760912590558.
 235          * (We could get a more accurate k by invoking log10,
 236          *  but this is probably not worthwhile.)
 237          */
 238 #ifdef IBM
 239         i <<= 2;
 240         i += j;
 241 #endif
 242         ds = (dval(d)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
 243
 244         /* correct assumption about exponent range */
 245         if ((j = i) < 0)
 246                 j = -j;
 247         if ((j -= 1077) > 0)
 248                 ds += j * 7e-17;
 249
 250         k = (int)ds;
 251         if (ds < 0. && ds != k)
 252                 k--;    /* want k = floor(ds) */
 253         k_check = 1;
 254 #ifdef IBM
 255         j = be + bbits - 1;
 256         if ( (jj1 = j & 3) !=0)
 257                 dval(d) *= 1 << jj1;
 258         word0(d) += j << Exp_shift - 2 & Exp_mask;
 259 #else
 260         word0(d) += (be + bbits - 1) << Exp_shift;
 261 #endif
 262         if (k >= 0 && k <= Ten_pmax) {
 263                 if (dval(d) < tens[k])
 264                         k--;
 265                 k_check = 0;
 266                 }
 267         j = bbits - i - 1;
 268         if (j >= 0) {
 269                 b2 = 0;
 270                 s2 = j;
 271                 }
 272         else {
 273                 b2 = -j;
 274                 s2 = 0;
 275                 }
 276         if (k >= 0) {
 277                 b5 = 0;
 278                 s5 = k;
 279                 s2 += k;
 280                 }
 281         else {
 282                 b2 -= k;
 283                 b5 = -k;
 284                 s5 = 0;
 285                 }
 286         if (mode < 0 || mode > 9)
 287                 mode = 0;
 288         try_quick = 1;
 289         if (mode > 5) {
 290                 mode -= 4;
 291                 try_quick = 0;
 292                 }
 293         leftright = 1;
 294         switch(mode) {
 295                 case 0:
 296                 case 1:
 297                         ilim = ilim1 = -1;
 298                         i = (int)(nbits * .30103) + 3;
 299                         ndigits = 0;
 300                         break;
 301                 case 2:
 302                         leftright = 0;
 303                         /*FALLTHROUGH*/
 304                 case 4:
 305                         if (ndigits <= 0)
 306                                 ndigits = 1;
 307                         ilim = ilim1 = i = ndigits;
 308                         break;
 309                 case 3:
 310                         leftright = 0;
 311                         /*FALLTHROUGH*/
 312                 case 5:
 313                         i = ndigits + k + 1;
 314                         ilim = i;
 315                         ilim1 = i - 1;
 316                         if (i <= 0)
 317                                 i = 1;
 318                 }
 319         s = s0 = rv_alloc((size_t)i);
 320         if (s == NULL)
 321                 return NULL;
 322
 323         if ( (rdir = fpi->rounding - 1) !=0) {
 324                 if (rdir < 0)
 325                         rdir = 2;
 326                 if (kind & STRTOG_Neg)
 327                         rdir = 3 - rdir;
 328                 }
 329
 330         /* Now rdir = 0 ==> round near, 1 ==> round up, 2 ==> round down. */
 331
 332         if (ilim >= 0 && ilim <= Quick_max && try_quick && !rdir
 333 #ifndef IMPRECISE_INEXACT
 334                 && k == 0
 335 #endif
 336                                                                 ) {
 337
 338                 /* Try to get by with floating-point arithmetic. */
 339
 340                 i = 0;
 341                 d2 = dval(d);
 342 #ifdef IBM
 343                 if ( (j = 11 - hi0bits(word0(d) & Frac_mask)) !=0)
 344                         dval(d) /= 1 << j;
 345 #endif
 346                 k0 = k;
 347                 ilim0 = ilim;
 348                 ieps = 2; /* conservative */
 349                 if (k > 0) {
 350                         ds = tens[k&0xf];
 351                         j = (unsigned int)k >> 4;
 352                         if (j & Bletch) {
 353                                 /* prevent overflows */
 354                                 j &= Bletch - 1;
 355                                 dval(d) /= bigtens[n_bigtens-1];
 356                                 ieps++;
 357                                 }
 358                         for(; j; j /= 2, i++)
 359                                 if (j & 1) {
 360                                         ieps++;
 361                                         ds *= bigtens[i];
 362                                         }
 363                         }
 364                 else  {
 365                         ds = 1.;
 366                         if ( (jj1 = -k) !=0) {
 367                                 dval(d) *= tens[jj1 & 0xf];
 368                                 for(j = jj1 >> 4; j; j >>= 1, i++)
 369                                         if (j & 1) {
 370                                                 ieps++;
 371                                                 dval(d) *= bigtens[i];
 372                                                 }
 373                                 }
 374                         }
 375                 if (k_check && dval(d) < 1. && ilim > 0) {
 376                         if (ilim1 <= 0)
 377                                 goto fast_failed;
 378                         ilim = ilim1;
 379                         k--;
 380                         dval(d) *= 10.;
 381                         ieps++;
 382                         }
 383                 dval(eps) = ieps*dval(d) + 7.;
 384                 word0(eps) -= (P-1)*Exp_msk1;
 385                 if (ilim == 0) {
 386                         S = mhi = 0;
 387                         dval(d) -= 5.;
 388                         if (dval(d) > dval(eps))
 389                                 goto one_digit;
 390                         if (dval(d) < -dval(eps))
 391                                 goto no_digits;
 392                         goto fast_failed;
 393                         }
 394 #ifndef No_leftright
 395                 if (leftright) {
 396                         /* Use Steele & White method of only
 397                          * generating digits needed.
 398                          */
 399                         dval(eps) = ds*0.5/tens[ilim-1] - dval(eps);
 400                         for(i = 0;;) {
 401                                 L = (Long)(dval(d)/ds);
 402                                 dval(d) -= L*ds;
 403                                 *s++ = '0' + (int)L;
 404                                 if (dval(d) < dval(eps)) {
 405                                         if (dval(d))
 406                                                 inex = STRTOG_Inexlo;
 407                                         goto ret1;
 408                                         }
 409                                 if (ds - dval(d) < dval(eps))
 410                                         goto bump_up;
 411                                 if (++i >= ilim)
 412                                         break;
 413                                 dval(eps) *= 10.;
 414                                 dval(d) *= 10.;
 415                                 }
 416                         }
 417                 else {
 418 #endif
 419                         /* Generate ilim digits, then fix them up. */
 420                         dval(eps) *= tens[ilim-1];
 421                         for(i = 1;; i++, dval(d) *= 10.) {
 422                                 if ( (L = (Long)(dval(d)/ds)) !=0)
 423                                         dval(d) -= L*ds;
 424                                 *s++ = '0' + (int)L;
 425                                 if (i == ilim) {
 426                                         ds *= 0.5;
 427                                         if (dval(d) > ds + dval(eps))
 428                                                 goto bump_up;
 429                                         else if (dval(d) < ds - dval(eps)) {
 430                                                 while(*--s == '0'){}
 431                                                 s++;
 432                                                 if (dval(d))
 433                                                         inex = STRTOG_Inexlo;
 434                                                 goto ret1;
 435                                                 }
 436                                         break;
 437                                         }
 438                                 }
 439 #ifndef No_leftright
 440                         }
 441 #endif
 442  fast_failed:
 443                 s = s0;
 444                 dval(d) = d2;
 445                 k = k0;
 446                 ilim = ilim0;
 447                 }
 448
 449         /* Do we have a "small" integer? */
 450
 451         if (be >= 0 && k <= Int_max) {
 452                 /* Yes. */
 453                 ds = tens[k];
 454                 if (ndigits < 0 && ilim <= 0) {
 455                         S = mhi = 0;
 456                         if (ilim < 0 || dval(d) <= 5*ds)
 457                                 goto no_digits;
 458                         goto one_digit;
 459                         }
 460                 for(i = 1;; i++, dval(d) *= 10.) {
 461                         L = dval(d) / ds;
 462                         dval(d) -= L*ds;
 463 #ifdef Check_FLT_ROUNDS
 464                         /* If FLT_ROUNDS == 2, L will usually be high by 1 */
 465                         if (dval(d) < 0) {
 466                                 L--;
 467                                 dval(d) += ds;
 468                                 }
 469 #endif
 470                         *s++ = '0' + (int)L;
 471                         if (dval(d) == 0.)
 472                                 break;
 473                         if (i == ilim) {
 474                                 if (rdir) {
 475                                         if (rdir == 1)
 476                                                 goto bump_up;
 477                                         inex = STRTOG_Inexlo;
 478                                         goto ret1;
 479                                         }
 480                                 dval(d) += dval(d);
 481                                 if (dval(d) > ds || (dval(d) == ds && L & 1)) {
 482  bump_up:
 483                                         inex = STRTOG_Inexhi;
 484                                         while(*--s == '9')
 485                                                 if (s == s0) {
 486                                                         k++;
 487                                                         *s = '0';
 488                                                         break;
 489                                                         }
 490                                         ++*s++;
 491                                         }
 492                                 else
 493                                         inex = STRTOG_Inexlo;
 494                                 break;
 495                                 }
 496                         }
 497                 goto ret1;
 498                 }
 499
 500         m2 = b2;
 501         m5 = b5;
 502         mhi = mlo = 0;
 503         if (leftright) {
 504                 if (mode < 2) {
 505                         i = nbits - bbits;
 506                         if (be - i++ < fpi->emin)
 507                                 /* denormal */
 508                                 i = be - fpi->emin + 1;
 509                         }
 510                 else {
 511                         j = ilim - 1;
 512                         if (m5 >= j)
 513                                 m5 -= j;
 514                         else {
 515                                 s5 += j -= m5;
 516                                 b5 += j;
 517                                 m5 = 0;
 518                                 }
 519                         if ((i = ilim) < 0) {
 520                                 m2 -= i;
 521                                 i = 0;
 522                                 }
 523                         }
 524                 b2 += i;
 525                 s2 += i;
 526                 mhi = i2b(1);
 527                 }
 528         if (m2 > 0 && s2 > 0) {
 529                 i = m2 < s2 ? m2 : s2;
 530                 b2 -= i;
 531                 m2 -= i;
 532                 s2 -= i;
 533                 }
 534         if (b5 > 0) {
 535                 if (leftright) {
 536                         if (m5 > 0) {
 537                                 mhi = pow5mult(mhi, m5);
 538                                 if (mhi == NULL)
 539                                         return NULL;
 540                                 b1 = mult(mhi, b);
 541                                 if (b1 == NULL)
 542                                         return NULL;
 543                                 Bfree(b);
 544                                 b = b1;
 545                                 }
 546                         if ( (j = b5 - m5) !=0) {
 547                                 b = pow5mult(b, j);
 548                                 if (b == NULL)
 549                                         return NULL;
 550                                 }
 551                         }
 552                 else {
 553                         b = pow5mult(b, b5);
 554                         if (b == NULL)
 555                                 return NULL;
 556                         }
 557                 }
 558         S = i2b(1);
 559         if (S == NULL)
 560                 return NULL;
 561         if (s5 > 0) {
 562                 S = pow5mult(S, s5);
 563                 if (S == NULL)
 564                         return NULL;
 565                 }
 566
 567         /* Check for special case that d is a normalized power of 2. */
 568
 569         spec_case = 0;
 570         if (mode < 2) {
 571                 if (bbits == 1 && be0 > fpi->emin + 1) {
 572                         /* The special case */
 573                         b2++;
 574                         s2++;
 575                         spec_case = 1;
 576                         }
 577                 }
 578
 579         /* Arrange for convenient computation of quotients:
 580          * shift left if necessary so divisor has 4 leading 0 bits.
 581          *
 582          * Perhaps we should just compute leading 28 bits of S once
 583          * and for all and pass them and a shift to quorem, so it
 584          * can do shifts and ors to compute the numerator for q.
 585          */
 586 #ifdef Pack_32
 587         if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) !=0)
 588                 i = 32 - i;
 589 #else
 590         if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) !=0)
 591                 i = 16 - i;
 592 #endif
 593         if (i > 4) {
 594                 i -= 4;
 595                 b2 += i;
 596                 m2 += i;
 597                 s2 += i;
 598                 }
 599         else if (i < 4) {
 600                 i += 28;
 601                 b2 += i;
 602                 m2 += i;
 603                 s2 += i;
 604                 }
 605         if (b2 > 0)
 606                 b = lshift(b, b2);
 607         if (s2 > 0)
 608                 S = lshift(S, s2);
 609         if (k_check) {
 610                 if (cmp(b,S) < 0) {
 611                         k--;
 612                         b = multadd(b, 10, 0);  /* we botched the k estimate */
 613                         if (b == NULL)
 614                                 return NULL;
 615                         if (leftright) {
 616                                 mhi = multadd(mhi, 10, 0);
 617                                 if (mhi == NULL)
 618                                         return NULL;
 619                                 }
 620                         ilim = ilim1;
 621                         }
 622                 }
 623         if (ilim <= 0 && mode > 2) {
 624                 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
 625                         /* no digits, fcvt style */
 626  no_digits:
 627                         k = -1 - ndigits;
 628                         inex = STRTOG_Inexlo;
 629                         goto ret;
 630                         }
 631  one_digit:
 632                 inex = STRTOG_Inexhi;
 633                 *s++ = '1';
 634                 k++;
 635                 goto ret;
 636                 }
 637         if (leftright) {
 638                 if (m2 > 0) {
 639                         mhi = lshift(mhi, m2);
 640                         if (mhi == NULL)
 641                                 return NULL;
 642                         }
 643
 644                 /* Compute mlo -- check for special case
 645                  * that d is a normalized power of 2.
 646                  */
 647
 648                 mlo = mhi;
 649                 if (spec_case) {
 650                         mhi = Balloc(mhi->k);
 651                         if (mhi == NULL)
 652                                 return NULL;
 653                         Bcopy(mhi, mlo);
 654                         mhi = lshift(mhi, 1);
 655                         if (mhi == NULL)
 656                                 return NULL;
 657                         }
 658
 659                 for(i = 1;;i++) {
 660                         dig = quorem(b,S) + '0';
 661                         /* Do we yet have the shortest decimal string
 662                          * that will round to d?
 663                          */
 664                         j = cmp(b, mlo);
 665                         delta = diff(S, mhi);
 666                         if (delta == NULL)
 667                                 return NULL;
 668                         jj1 = delta->sign ? 1 : cmp(b, delta);
 669                         Bfree(delta);
 670 #ifndef ROUND_BIASED
 671                         if (jj1 == 0 && !mode && !(bits[0] & 1) && !rdir) {
 672                                 if (dig == '9')
 673                                         goto round_9_up;
 674                                 if (j <= 0) {
 675                                         if (b->wds > 1 || b->x[0])
 676                                                 inex = STRTOG_Inexlo;
 677                                         }
 678                                 else {
 679                                         dig++;
 680                                         inex = STRTOG_Inexhi;
 681                                         }
 682                                 *s++ = dig;
 683                                 goto ret;
 684                                 }
 685 #endif
 686                         if (j < 0 || (j == 0 && !mode
 687 #ifndef ROUND_BIASED
 688                                                         && !(bits[0] & 1)
 689 #endif
 690                                         )) {
 691                                 if (rdir && (b->wds > 1 || b->x[0])) {
 692                                         if (rdir == 2) {
 693                                                 inex = STRTOG_Inexlo;
 694                                                 goto accept;
 695                                                 }
 696                                         while (cmp(S,mhi) > 0) {
 697                                                 *s++ = dig;
 698                                                 mhi1 = multadd(mhi, 10, 0);
 699                                                 if (mhi1 == NULL)
 700                                                         return NULL;
 701                                                 if (mlo == mhi)
 702                                                         mlo = mhi1;
 703                                                 mhi = mhi1;
 704                                                 b = multadd(b, 10, 0);
 705                                                 if (b == NULL)
 706                                                         return NULL;
 707                                                 dig = quorem(b,S) + '0';
 708                                                 }
 709                                         if (dig++ == '9')
 710                                                 goto round_9_up;
 711                                         inex = STRTOG_Inexhi;
 712                                         goto accept;
 713                                         }
 714                                 if (jj1 > 0) {
 715                                         b = lshift(b, 1);
 716                                         if (b == NULL)
 717                                                 return NULL;
 718                                         jj1 = cmp(b, S);
 719                                         if ((jj1 > 0 || (jj1 == 0 && dig & 1))
 720                                         && dig++ == '9')
 721                                                 goto round_9_up;
 722                                         inex = STRTOG_Inexhi;
 723                                         }
 724                                 if (b->wds > 1 || b->x[0])
 725                                         inex = STRTOG_Inexlo;
 726  accept:
 727                                 *s++ = dig;
 728                                 goto ret;
 729                                 }
 730                         if (jj1 > 0 && rdir != 2) {
 731                                 if (dig == '9') { /* possible if i == 1 */
 732  round_9_up:
 733                                         *s++ = '9';
 734                                         inex = STRTOG_Inexhi;
 735                                         goto roundoff;
 736                                         }
 737                                 inex = STRTOG_Inexhi;
 738                                 *s++ = dig + 1;
 739                                 goto ret;
 740                                 }
 741                         *s++ = dig;
 742                         if (i == ilim)
 743                                 break;
 744                         b = multadd(b, 10, 0);
 745                         if (b == NULL)
 746                                 return NULL;
 747                         if (mlo == mhi) {
 748                                 mlo = mhi = multadd(mhi, 10, 0);
 749                                 if (mlo == NULL)
 750                                         return NULL;
 751                                 }
 752                         else {
 753                                 mlo = multadd(mlo, 10, 0);
 754                                 if (mlo == NULL)
 755                                         return NULL;
 756                                 mhi = multadd(mhi, 10, 0);
 757                                 if (mhi == NULL)
 758                                         return NULL;
 759                                 }
 760                         }
 761                 }
 762         else
 763                 for(i = 1;; i++) {
 764                         *s++ = dig = quorem(b,S) + '0';
 765                         if (i >= ilim)
 766                                 break;
 767                         b = multadd(b, 10, 0);
 768                         if (b == NULL)
 769                                 return NULL;
 770                         }
 771
 772         /* Round off last digit */
 773
 774         if (rdir) {
 775                 if (rdir == 2 || (b->wds <= 1 && !b->x[0]))
 776                         goto chopzeros;
 777                 goto roundoff;
 778                 }
 779         b = lshift(b, 1);
 780         if (b == NULL)
 781                 return NULL;
 782         j = cmp(b, S);
 783         if (j > 0 || (j == 0 && dig & 1)) {
 784  roundoff:
 785                 inex = STRTOG_Inexhi;
 786                 while(*--s == '9')
 787                         if (s == s0) {
 788                                 k++;
 789                                 *s++ = '1';
 790                                 goto ret;
 791                                 }
 792                 ++*s++;
 793                 }
 794         else {
 795  chopzeros:
 796                 if (b->wds > 1 || b->x[0])
 797                         inex = STRTOG_Inexlo;
 798                 while(*--s == '0'){}
 799                 s++;
 800                 }
 801  ret:
 802         Bfree(S);
 803         if (mhi) {
 804                 if (mlo && mlo != mhi)
 805                         Bfree(mlo);
 806                 Bfree(mhi);
 807                 }
 808  ret1:
 809         Bfree(b);
 810         *s = 0;
 811         *decpt = k + 1;
 812         if (rve)
 813                 *rve = s;
 814         *kindp |= inex;
 815         return s0;
 816         }