gonzui

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