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