StdLib/LibC/gdtoa/gdtoa.c

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