gonzui

    1: /****************************************************************
    2: 
    3: The author of this software is David M. Gay.
    4: 
    5: Copyright (C) 1998-2001 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: #ifndef NO_FENV_H
   34: #include <fenv.h>
   35: #endif
   36: 
   37: #ifdef USE_LOCALE
   38: #include "locale.h"
   39: #endif
   40: 
   41: #ifdef IEEE_Arith
   42: #ifndef NO_IEEE_Scale
   43: #define Avoid_Underflow
   44: #undef tinytens
   45: /* The factor of 2^106 in tinytens[4] helps us avoid setting the underflow */
   46: /* flag unnecessarily.  It leads to a song and dance at the end of strtod. */
   47: static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
   48:                 9007199254740992.*9007199254740992.e-256
   49:                 };
   50: #endif
   51: #endif
   52: 
   53: #ifdef Honor_FLT_ROUNDS
   54: #undef Check_FLT_ROUNDS
   55: #define Check_FLT_ROUNDS
   56: #else
   57: #define Rounding Flt_Rounds
   58: #endif
   59: 
   60: #ifdef Avoid_Underflow /*{*/
   61:  static double
   62: sulp
   63: #ifdef KR_headers
   64:         (x, scale) U *x; int scale;
   65: #else
   66:         (U *x, int scale)
   67: #endif
   68: {
   69:         U u;
   70:         double rv;
   71:         int i;
   72: 
   73:         rv = ulp(x);
   74:         if (!scale || (i = 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift)) <= 0)
   75:                 return rv; /* Is there an example where i <= 0 ? */
   76:         word0(&u) = Exp_1 + (i << Exp_shift);
   77:         word1(&u) = 0;
   78:         return rv * u.d;
   79:         }
   80: #endif /*}*/
   81: 
   82:  double
   83: strtod
   84: #ifdef KR_headers
   85:         (s00, se) CONST char *s00; char **se;
   86: #else
   87:         (CONST char *s00, char **se)
   88: #endif
   89: {
   90: #ifdef Avoid_Underflow
   91:         int scale;
   92: #endif
   93:         int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign,
   94:                  e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
   95:         CONST char *s, *s0, *s1;
   96:         double aadj;
   97:         Long L;
   98:         U adj, aadj1, rv, rv0;
   99:         ULong y, z;
  100:         Bigint *bb = NULL, *bb1, *bd = NULL, *bd0 = NULL, *bs = NULL, *delta = NULL;
  101: #ifdef Avoid_Underflow
  102:         ULong Lsb, Lsb1;
  103: #endif
  104: #ifdef SET_INEXACT
  105:         int inexact, oldinexact;
  106: #endif
  107: #ifdef USE_LOCALE /*{{*/
  108: #ifdef NO_LOCALE_CACHE
  109:         char *decimalpoint = localeconv()->decimal_point;
  110:         int dplen = strlen(decimalpoint);
  111: #else
  112:         char *decimalpoint;
  113:         static char *decimalpoint_cache;
  114:         static int dplen;
  115:         if (!(s0 = decimalpoint_cache)) {
  116:                 s0 = localeconv()->decimal_point;
  117:                 if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
  118:                         strlcpy(decimalpoint_cache, s0, strlen(s0) + 1);
  119:                         s0 = decimalpoint_cache;
  120:                         }
  121:                 dplen = strlen(s0);
  122:                 }
  123:         decimalpoint = (char*)s0;
  124: #endif /*NO_LOCALE_CACHE*/
  125: #else  /*USE_LOCALE}{*/
  126: #define dplen 1
  127: #endif /*USE_LOCALE}}*/
  128: 
  129: #ifdef Honor_FLT_ROUNDS /*{*/
  130:         int Rounding;
  131: #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
  132:         Rounding = Flt_Rounds;
  133: #else /*}{*/
  134:         Rounding = 1;
  135:         switch(fegetround()) {
  136:           case FE_TOWARDZERO:  Rounding = 0; break;
  137:           case FE_UPWARD:      Rounding = 2; break;
  138:           case FE_DOWNWARD:    Rounding = 3;
  139:           }
  140: #endif /*}}*/
  141: #endif /*}*/
  142: 
  143:         sign = nz0 = nz = decpt = 0;
  144:         dval(&rv) = 0.;
  145:         for(s = s00;;s++) switch(*s) {
  146:                 case '-':
  147:                         sign = 1;
  148:                         /* no break */
  149:                 case '+':
  150:                         if (*++s)
  151:                                 goto break2;
  152:                         /* no break */
  153:                 case 0:
  154:                         goto ret0;
  155:                 case '\t':
  156:                 case '\n':
  157:                 case '\v':
  158:                 case '\f':
  159:                 case '\r':
  160:                 case ' ':
  161:                         continue;
  162:                 default:
  163:                         goto break2;
  164:                 }
  165:  break2:
  166:         if (*s == '0') {
  167: #ifndef NO_HEX_FP /*{*/
  168:                 {
  169:                 static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
  170:                 Long exp;
  171:                 ULong bits[2];
  172:                 switch(s[1]) {
  173:                   case 'x':
  174:                   case 'X':
  175:                         {
  176: #ifdef Honor_FLT_ROUNDS
  177:                         FPI fpi1 = fpi;
  178:                         fpi1.rounding = Rounding;
  179: #else
  180: #define fpi1 fpi
  181: #endif
  182:                         switch((i = gethex(&s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) {
  183:                           case STRTOG_NoMemory:
  184:                                 goto ovfl;
  185:                           case STRTOG_NoNumber:
  186:                                 s = s00;
  187:                                 sign = 0;
  188:                           case STRTOG_Zero:
  189:                                 break;
  190:                           default:
  191:                                 if (bb) {
  192:                                         copybits(bits, fpi.nbits, bb);
  193:                                         Bfree(bb);
  194:                                         }
  195:                                 ULtod(((U*)&rv)->L, bits, exp, i);
  196:                           }}
  197:                         goto ret;
  198:                   }
  199:                 }
  200: #endif /*}*/
  201:                 nz0 = 1;
  202:                 while(*++s == '0') ;
  203:                 if (!*s)
  204:                         goto ret;
  205:                 }
  206:         s0 = s;
  207:         y = z = 0;
  208:         for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
  209:                 if (nd < 9)
  210:                         y = 10*y + c - '0';
  211:                 else if (nd < 16)
  212:                         z = 10*z + c - '0';
  213:         nd0 = nd;
  214: #ifdef USE_LOCALE
  215:         if (c == *decimalpoint) {
  216:                 for(i = 1; decimalpoint[i]; ++i)
  217:                         if (s[i] != decimalpoint[i])
  218:                                 goto dig_done;
  219:                 s += i;
  220:                 c = *s;
  221: #else
  222:         if (c == '.') {
  223:                 c = *++s;
  224: #endif
  225:                 decpt = 1;
  226:                 if (!nd) {
  227:                         for(; c == '0'; c = *++s)
  228:                                 nz++;
  229:                         if (c > '0' && c <= '9') {
  230:                                 s0 = s;
  231:                                 nf += nz;
  232:                                 nz = 0;
  233:                                 goto have_dig;
  234:                                 }
  235:                         goto dig_done;
  236:                         }
  237:                 for(; c >= '0' && c <= '9'; c = *++s) {
  238:  have_dig:
  239:                         nz++;
  240:                         if (c -= '0') {
  241:                                 nf += nz;
  242:                                 for(i = 1; i < nz; i++)
  243:                                         if (nd++ < 9)
  244:                                                 y *= 10;
  245:                                         else if (nd <= DBL_DIG + 1)
  246:                                                 z *= 10;
  247:                                 if (nd++ < 9)
  248:                                         y = 10*y + c;
  249:                                 else if (nd <= DBL_DIG + 1)
  250:                                         z = 10*z + c;
  251:                                 nz = 0;
  252:                                 }
  253:                         }
  254:                 }/*}*/
  255:  dig_done:
  256:         e = 0;
  257:         if (c == 'e' || c == 'E') {
  258:                 if (!nd && !nz && !nz0) {
  259:                         goto ret0;
  260:                         }
  261:                 s00 = s;
  262:                 esign = 0;
  263:                 switch(c = *++s) {
  264:                         case '-':
  265:                                 esign = 1;
  266:                         case '+':
  267:                                 c = *++s;
  268:                         }
  269:                 if (c >= '0' && c <= '9') {
  270:                         while(c == '0')
  271:                                 c = *++s;
  272:                         if (c > '0' && c <= '9') {
  273:                                 L = c - '0';
  274:                                 s1 = s;
  275:                                 while((c = *++s) >= '0' && c <= '9')
  276:                                         L = 10*L + c - '0';
  277:                                 if (s - s1 > 8 || L > 19999)
  278:                                         /* Avoid confusion from exponents
  279:                                          * so large that e might overflow.
  280:                                          */
  281:                                         e = 19999; /* safe for 16 bit ints */
  282:                                 else
  283:                                         e = (int)L;
  284:                                 if (esign)
  285:                                         e = -e;
  286:                                 }
  287:                         else
  288:                                 e = 0;
  289:                         }
  290:                 else
  291:                         s = s00;
  292:                 }
  293:         if (!nd) {
  294:                 if (!nz && !nz0) {
  295: #ifdef INFNAN_CHECK
  296:                         /* Check for Nan and Infinity */
  297:                         ULong bits[2];
  298:                         static FPI fpinan =  /* only 52 explicit bits */
  299:                                 { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI };
  300:                         if (!decpt)
  301:                          switch(c) {
  302:                           case 'i':
  303:                           case 'I':
  304:                                 if (match(&s,"nf")) {
  305:                                         --s;
  306:                                         if (!match(&s,"inity"))
  307:                                                 ++s;
  308:                                         word0(&rv) = 0x7ff00000;
  309:                                         word1(&rv) = 0;
  310:                                         goto ret;
  311:                                         }
  312:                                 break;
  313:                           case 'n':
  314:                           case 'N':
  315:                                 if (match(&s, "an")) {
  316: #ifndef No_Hex_NaN
  317:                                         if (*s == '(' /*)*/
  318:                                          && hexnan(&s, &fpinan, bits)
  319:                                                         == STRTOG_NaNbits) {
  320:                                                 word0(&rv) = 0x7ff00000 | bits[1];
  321:                                                 word1(&rv) = bits[0];
  322:                                                 }
  323:                                         else {
  324: #endif
  325:                                                 word0(&rv) = NAN_WORD0;
  326:                                                 word1(&rv) = NAN_WORD1;
  327: #ifndef No_Hex_NaN
  328:                                                 }
  329: #endif
  330:                                         goto ret;
  331:                                         }
  332:                           }
  333: #endif /* INFNAN_CHECK */
  334:  ret0:
  335:                         s = s00;
  336:                         sign = 0;
  337:                         }
  338:                 goto ret;
  339:                 }
  340:         e1 = e -= nf;
  341: 
  342:         /* Now we have nd0 digits, starting at s0, followed by a
  343:          * decimal point, followed by nd-nd0 digits.  The number we're
  344:          * after is the integer represented by those digits times
  345:          * 10**e */
  346: 
  347:         if (!nd0)
  348:                 nd0 = nd;
  349:         k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
  350:         dval(&rv) = y;
  351:         if (k > 9) {
  352: #ifdef SET_INEXACT
  353:                 if (k > DBL_DIG)
  354:                         oldinexact = get_inexact();
  355: #endif
  356:                 dval(&rv) = tens[k - 9] * dval(&rv) + z;
  357:                 }
  358:         if (nd <= DBL_DIG
  359: #ifndef RND_PRODQUOT
  360: #ifndef Honor_FLT_ROUNDS
  361:                 && Flt_Rounds == 1
  362: #endif
  363: #endif
  364:                         ) {
  365:                 if (!e)
  366:                         goto ret;
  367: #ifndef ROUND_BIASED_without_Round_Up
  368:                 if (e > 0) {
  369:                         if (e <= Ten_pmax) {
  370: #ifdef VAX
  371:                                 goto vax_ovfl_check;
  372: #else
  373: #ifdef Honor_FLT_ROUNDS
  374:                                 /* round correctly FLT_ROUNDS = 2 or 3 */
  375:                                 if (sign) {
  376:                                         rv.d = -rv.d;
  377:                                         sign = 0;
  378:                                         }
  379: #endif
  380:                                 /* rv = */ rounded_product(dval(&rv), tens[e]);
  381:                                 goto ret;
  382: #endif
  383:                                 }
  384:                         i = DBL_DIG - nd;
  385:                         if (e <= Ten_pmax + i) {
  386:                                 /* A fancier test would sometimes let us do
  387:                                  * this for larger i values.
  388:                                  */
  389: #ifdef Honor_FLT_ROUNDS
  390:                                 /* round correctly FLT_ROUNDS = 2 or 3 */
  391:                                 if (sign) {
  392:                                         rv.d = -rv.d;
  393:                                         sign = 0;
  394:                                         }
  395: #endif
  396:                                 e -= i;
  397:                                 dval(&rv) *= tens[i];
  398: #ifdef VAX
  399:                                 /* VAX exponent range is so narrow we must
  400:                                  * worry about overflow here...
  401:                                  */
  402:  vax_ovfl_check:
  403:                                 word0(&rv) -= P*Exp_msk1;
  404:                                 /* rv = */ rounded_product(dval(&rv), tens[e]);
  405:                                 if ((word0(&rv) & Exp_mask)
  406:                                  > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
  407:                                         goto ovfl;
  408:                                 word0(&rv) += P*Exp_msk1;
  409: #else
  410:                                 /* rv = */ rounded_product(dval(&rv), tens[e]);
  411: #endif
  412:                                 goto ret;
  413:                                 }
  414:                         }
  415: #ifndef Inaccurate_Divide
  416:                 else if (e >= -Ten_pmax) {
  417: #ifdef Honor_FLT_ROUNDS
  418:                         /* round correctly FLT_ROUNDS = 2 or 3 */
  419:                         if (sign) {
  420:                                 rv.d = -rv.d;
  421:                                 sign = 0;
  422:                                 }
  423: #endif
  424:                         /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
  425:                         goto ret;
  426:                         }
  427: #endif
  428: #endif /* ROUND_BIASED_without_Round_Up */
  429:                 }
  430:         e1 += nd - k;
  431: 
  432: #ifdef IEEE_Arith
  433: #ifdef SET_INEXACT
  434:         inexact = 1;
  435:         if (k <= DBL_DIG)
  436:                 oldinexact = get_inexact();
  437: #endif
  438: #ifdef Avoid_Underflow
  439:         scale = 0;
  440: #endif
  441: #ifdef Honor_FLT_ROUNDS
  442:         if (Rounding >= 2) {
  443:                 if (sign)
  444:                         Rounding = Rounding == 2 ? 0 : 2;
  445:                 else
  446:                         if (Rounding != 2)
  447:                                 Rounding = 0;
  448:                 }
  449: #endif
  450: #endif /*IEEE_Arith*/
  451: 
  452:         /* Get starting approximation = rv * 10**e1 */
  453: 
  454:         if (e1 > 0) {
  455:                 if ( (i = e1 & 15) !=0)
  456:                         dval(&rv) *= tens[i];
  457:                 if (e1 &= ~15) {
  458:                         if (e1 > DBL_MAX_10_EXP) {
  459:  ovfl:
  460:                                 /* Can't trust HUGE_VAL */
  461: #ifdef IEEE_Arith
  462: #ifdef Honor_FLT_ROUNDS
  463:                                 switch(Rounding) {
  464:                                   case 0: /* toward 0 */
  465:                                   case 3: /* toward -infinity */
  466:                                         word0(&rv) = Big0;
  467:                                         word1(&rv) = Big1;
  468:                                         break;
  469:                                   default:
  470:                                         word0(&rv) = Exp_mask;
  471:                                         word1(&rv) = 0;
  472:                                   }
  473: #else /*Honor_FLT_ROUNDS*/
  474:                                 word0(&rv) = Exp_mask;
  475:                                 word1(&rv) = 0;
  476: #endif /*Honor_FLT_ROUNDS*/
  477: #ifdef SET_INEXACT
  478:                                 /* set overflow bit */
  479:                                 dval(&rv0) = 1e300;
  480:                                 dval(&rv0) *= dval(&rv0);
  481: #endif
  482: #else /*IEEE_Arith*/
  483:                                 word0(&rv) = Big0;
  484:                                 word1(&rv) = Big1;
  485: #endif /*IEEE_Arith*/
  486:  range_err:
  487:                                 if (bd0) {
  488:                                         Bfree(bb);
  489:                                         Bfree(bd);
  490:                                         Bfree(bs);
  491:                                         Bfree(bd0);
  492:                                         Bfree(delta);
  493:                                         }
  494: #ifndef NO_ERRNO
  495:                                 /*errno = ERANGE*/;
  496: #endif
  497:                                 goto ret;
  498:                                 }
  499:                         e1 >>= 4;
  500:                         for(j = 0; e1 > 1; j++, e1 >>= 1)
  501:                                 if (e1 & 1)
  502:                                         dval(&rv) *= bigtens[j];
  503:                 /* The last multiplication could overflow. */
  504:                         word0(&rv) -= P*Exp_msk1;
  505:                         dval(&rv) *= bigtens[j];
  506:                         if ((z = word0(&rv) & Exp_mask)
  507:                          > Exp_msk1*(DBL_MAX_EXP+Bias-P))
  508:                                 goto ovfl;
  509:                         if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
  510:                                 /* set to largest number */
  511:                                 /* (Can't trust DBL_MAX) */
  512:                                 word0(&rv) = Big0;
  513:                                 word1(&rv) = Big1;
  514:                                 }
  515:                         else
  516:                                 word0(&rv) += P*Exp_msk1;
  517:                         }
  518:                 }
  519:         else if (e1 < 0) {
  520:                 e1 = -e1;
  521:                 if ( (i = e1 & 15) !=0)
  522:                         dval(&rv) /= tens[i];
  523:                 if (e1 >>= 4) {
  524:                         if (e1 >= 1 << n_bigtens)
  525:                                 goto undfl;
  526: #ifdef Avoid_Underflow
  527:                         if (e1 & Scale_Bit)
  528:                                 scale = 2*P;
  529:                         for(j = 0; e1 > 0; j++, e1 >>= 1)
  530:                                 if (e1 & 1)
  531:                                         dval(&rv) *= tinytens[j];
  532:                         if (scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
  533:                                                 >> Exp_shift)) > 0) {
  534:                                 /* scaled rv is denormal; zap j low bits */
  535:                                 if (j >= 32) {
  536:                                         word1(&rv) = 0;
  537:                                         if (j >= 53)
  538:                                          word0(&rv) = (P+2)*Exp_msk1;
  539:                                         else
  540:                                          word0(&rv) &= 0xffffffff << (j-32);
  541:                                         }
  542:                                 else
  543:                                         word1(&rv) &= 0xffffffff << j;
  544:                                 }
  545: #else
  546:                         for(j = 0; e1 > 1; j++, e1 >>= 1)
  547:                                 if (e1 & 1)
  548:                                         dval(&rv) *= tinytens[j];
  549:                         /* The last multiplication could underflow. */
  550:                         dval(&rv0) = dval(&rv);
  551:                         dval(&rv) *= tinytens[j];
  552:                         if (!dval(&rv)) {
  553:                                 dval(&rv) = 2.*dval(&rv0);
  554:                                 dval(&rv) *= tinytens[j];
  555: #endif
  556:                                 if (!dval(&rv)) {
  557:  undfl:
  558:                                         dval(&rv) = 0.;
  559:                                         goto range_err;
  560:                                         }
  561: #ifndef Avoid_Underflow
  562:                                 word0(&rv) = Tiny0;
  563:                                 word1(&rv) = Tiny1;
  564:                                 /* The refinement below will clean
  565:                                  * this approximation up.
  566:                                  */
  567:                                 }
  568: #endif
  569:                         }
  570:                 }
  571: 
  572:         /* Now the hard part -- adjusting rv to the correct value.*/
  573: 
  574:         /* Put digits into bd: true value = bd * 10^e */
  575: 
  576:         bd0 = s2b(s0, nd0, nd, y, dplen);
  577:         if (bd0 == NULL)
  578:                 goto ovfl;
  579: 
  580:         for(;;) {
  581:                 bd = Balloc(bd0->k);
  582:                 if (bd == NULL)
  583:                         goto ovfl;
  584:                 Bcopy(bd, bd0);
  585:                 bb = d2b(dval(&rv), &bbe, &bbbits);   /* rv = bb * 2^bbe */
  586:                 if (bb == NULL)
  587:                         goto ovfl;
  588:                 bs = i2b(1);
  589:                 if (bs == NULL)
  590:                         goto ovfl;
  591: 
  592:                 if (e >= 0) {
  593:                         bb2 = bb5 = 0;
  594:                         bd2 = bd5 = e;
  595:                         }
  596:                 else {
  597:                         bb2 = bb5 = -e;
  598:                         bd2 = bd5 = 0;
  599:                         }
  600:                 if (bbe >= 0)
  601:                         bb2 += bbe;
  602:                 else
  603:                         bd2 -= bbe;
  604:                 bs2 = bb2;
  605: #ifdef Honor_FLT_ROUNDS
  606:                 if (Rounding != 1)
  607:                         bs2++;
  608: #endif
  609: #ifdef Avoid_Underflow
  610:                 Lsb = LSB;
  611:                 Lsb1 = 0;
  612:                 j = bbe - scale;
  613:                 i = j + bbbits - 1;   /* logb(rv) */
  614:                 j = P + 1 - bbbits;
  615:                 if (i < Emin) {       /* denormal */
  616:                         i = Emin - i;
  617:                         j -= i;
  618:                         if (i < 32)
  619:                                 Lsb <<= i;
  620:                         else
  621:                                 Lsb1 = Lsb << (i-32);
  622:                         }
  623: #else /*Avoid_Underflow*/
  624: #ifdef Sudden_Underflow
  625: #ifdef IBM
  626:                 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
  627: #else
  628:                 j = P + 1 - bbbits;
  629: #endif
  630: #else /*Sudden_Underflow*/
  631:                 j = bbe;
  632:                 i = j + bbbits - 1;   /* logb(&rv) */
  633:                 if (i < Emin) /* denormal */
  634:                         j += P - Emin;
  635:                 else
  636:                         j = P + 1 - bbbits;
  637: #endif /*Sudden_Underflow*/
  638: #endif /*Avoid_Underflow*/
  639:                 bb2 += j;
  640:                 bd2 += j;
  641: #ifdef Avoid_Underflow
  642:                 bd2 += scale;
  643: #endif
  644:                 i = bb2 < bd2 ? bb2 : bd2;
  645:                 if (i > bs2)
  646:                         i = bs2;
  647:                 if (i > 0) {
  648:                         bb2 -= i;
  649:                         bd2 -= i;
  650:                         bs2 -= i;
  651:                         }
  652:                 if (bb5 > 0) {
  653:                         bs = pow5mult(bs, bb5);
  654:                         if (bs == NULL)
  655:                                 goto ovfl;
  656:                         bb1 = mult(bs, bb);
  657:                         if (bb1 == NULL)
  658:                                 goto ovfl;
  659:                         Bfree(bb);
  660:                         bb = bb1;
  661:                         }
  662:                 if (bb2 > 0) {
  663:                         bb = lshift(bb, bb2);
  664:                         if (bb == NULL)
  665:                                 goto ovfl;
  666:                         }
  667:                 if (bd5 > 0) {
  668:                         bd = pow5mult(bd, bd5);
  669:                         if (bd == NULL)
  670:                                 goto ovfl;
  671:                         }
  672:                 if (bd2 > 0) {
  673:                         bd = lshift(bd, bd2);
  674:                         if (bd == NULL)
  675:                                 goto ovfl;
  676:                         }
  677:                 if (bs2 > 0) {
  678:                         bs = lshift(bs, bs2);
  679:                         if (bs == NULL)
  680:                                 goto ovfl;
  681:                         }
  682:                 delta = diff(bb, bd);
  683:                 if (delta == NULL)
  684:                         goto ovfl;
  685:                 dsign = delta->sign;
  686:                 delta->sign = 0;
  687:                 i = cmp(delta, bs);
  688: #ifdef Honor_FLT_ROUNDS
  689:                 if (Rounding != 1) {
  690:                         if (i < 0) {
  691:                                 /* Error is less than an ulp */
  692:                                 if (!delta->x[0] && delta->wds <= 1) {
  693:                                         /* exact */
  694: #ifdef SET_INEXACT
  695:                                         inexact = 0;
  696: #endif
  697:                                         break;
  698:                                         }
  699:                                 if (Rounding) {
  700:                                         if (dsign) {
  701:                                                 dval(&adj) = 1.;
  702:                                                 goto apply_adj;
  703:                                                 }
  704:                                         }
  705:                                 else if (!dsign) {
  706:                                         dval(&adj) = -1.;
  707:                                         if (!word1(&rv)
  708:                                          && !(word0(&rv) & Frac_mask)) {
  709:                                                 y = word0(&rv) & Exp_mask;
  710: #ifdef Avoid_Underflow
  711:                                                 if (!scale || y > 2*P*Exp_msk1)
  712: #else
  713:                                                 if (y)
  714: #endif
  715:                                                   {
  716:                                                   delta = lshift(delta,Log2P);
  717:                                                   if (delta == NULL)
  718:                                                         goto ovfl;
  719:                                                   if (cmp(delta, bs) <= 0)
  720:                                                         dval(&adj) = -0.5;
  721:                                                   }
  722:                                                 }
  723:  apply_adj:
  724: #ifdef Avoid_Underflow
  725:                                         if (scale && (y = word0(&rv) & Exp_mask)
  726:                                                 <= 2*P*Exp_msk1)
  727:                                           word0(&adj) += (2*P+1)*Exp_msk1 - y;
  728: #else
  729: #ifdef Sudden_Underflow
  730:                                         if ((word0(&rv) & Exp_mask) <=
  731:                                                         P*Exp_msk1) {
  732:                                                 word0(&rv) += P*Exp_msk1;
  733:                                                 dval(&rv) += adj*ulp(&rv);
  734:                                                 word0(&rv) -= P*Exp_msk1;
  735:                                                 }
  736:                                         else
  737: #endif /*Sudden_Underflow*/
  738: #endif /*Avoid_Underflow*/
  739:                                         dval(&rv) += adj.d*ulp(&rv);
  740:                                         }
  741:                                 break;
  742:                                 }
  743:                         dval(&adj) = ratio(delta, bs);
  744:                         if (adj.d < 1.)
  745:                                 dval(&adj) = 1.;
  746:                         if (adj.d <= 0x7ffffffe) {
  747:                                 /* dval(&adj) = Rounding ? ceil(&adj) : floor(&adj); */
  748:                                 y = adj.d;
  749:                                 if (y != adj.d) {
  750:                                         if (!((Rounding>>1) ^ dsign))
  751:                                                 y++;
  752:                                         dval(&adj) = y;
  753:                                         }
  754:                                 }
  755: #ifdef Avoid_Underflow
  756:                         if (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
  757:                                 word0(&adj) += (2*P+1)*Exp_msk1 - y;
  758: #else
  759: #ifdef Sudden_Underflow
  760:                         if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
  761:                                 word0(&rv) += P*Exp_msk1;
  762:                                 dval(&adj) *= ulp(&rv);
  763:                                 if (dsign)
  764:                                         dval(&rv) += adj;
  765:                                 else
  766:                                         dval(&rv) -= adj;
  767:                                 word0(&rv) -= P*Exp_msk1;
  768:                                 goto cont;
  769:                                 }
  770: #endif /*Sudden_Underflow*/
  771: #endif /*Avoid_Underflow*/
  772:                         dval(&adj) *= ulp(&rv);
  773:                         if (dsign) {
  774:                                 if (word0(&rv) == Big0 && word1(&rv) == Big1)
  775:                                         goto ovfl;
  776:                                 dval(&rv) += adj.d;
  777:                                 }
  778:                         else
  779:                                 dval(&rv) -= adj.d;
  780:                         goto cont;
  781:                         }
  782: #endif /*Honor_FLT_ROUNDS*/
  783: 
  784:                 if (i < 0) {
  785:                         /* Error is less than half an ulp -- check for
  786:                          * special case of mantissa a power of two.
  787:                          */
  788:                         if (dsign || word1(&rv) || word0(&rv) & Bndry_mask
  789: #ifdef IEEE_Arith
  790: #ifdef Avoid_Underflow
  791:                          || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
  792: #else
  793:                          || (word0(&rv) & Exp_mask) <= Exp_msk1
  794: #endif
  795: #endif
  796:                                 ) {
  797: #ifdef SET_INEXACT
  798:                                 if (!delta->x[0] && delta->wds <= 1)
  799:                                         inexact = 0;
  800: #endif
  801:                                 break;
  802:                                 }
  803:                         if (!delta->x[0] && delta->wds <= 1) {
  804:                                 /* exact result */
  805: #ifdef SET_INEXACT
  806:                                 inexact = 0;
  807: #endif
  808:                                 break;
  809:                                 }
  810:                         delta = lshift(delta,Log2P);
  811:                         if (delta == NULL)
  812:                                 goto ovfl;
  813:                         if (cmp(delta, bs) > 0)
  814:                                 goto drop_down;
  815:                         break;
  816:                         }
  817:                 if (i == 0) {
  818:                         /* exactly half-way between */
  819:                         if (dsign) {
  820:                                 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
  821:                                  &&  word1(&rv) == (
  822: #ifdef Avoid_Underflow
  823:                         (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
  824:                 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
  825: #endif
  826:                                                    0xffffffff)) {
  827:                                         /*boundary case -- increment exponent*/
  828:                                         if (word0(&rv) == Big0 && word1(&rv) == Big1)
  829:                                                 goto ovfl;
  830:                                         word0(&rv) = (word0(&rv) & Exp_mask)
  831:                                                 + Exp_msk1
  832: #ifdef IBM
  833:                                                 | Exp_msk1 >> 4
  834: #endif
  835:                                                 ;
  836:                                         word1(&rv) = 0;
  837: #ifdef Avoid_Underflow
  838:                                         dsign = 0;
  839: #endif
  840:                                         break;
  841:                                         }
  842:                                 }
  843:                         else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
  844:  drop_down:
  845:                                 /* boundary case -- decrement exponent */
  846: #ifdef Sudden_Underflow /*{{*/
  847:                                 L = word0(&rv) & Exp_mask;
  848: #ifdef IBM
  849:                                 if (L <  Exp_msk1)
  850: #else
  851: #ifdef Avoid_Underflow
  852:                                 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
  853: #else
  854:                                 if (L <= Exp_msk1)
  855: #endif /*Avoid_Underflow*/
  856: #endif /*IBM*/
  857:                                         goto undfl;
  858:                                 L -= Exp_msk1;
  859: #else /*Sudden_Underflow}{*/
  860: #ifdef Avoid_Underflow
  861:                                 if (scale) {
  862:                                         L = word0(&rv) & Exp_mask;
  863:                                         if (L <= (2*P+1)*Exp_msk1) {
  864:                                                 if (L > (P+2)*Exp_msk1)
  865:                                                         /* round even ==> */
  866:                                                         /* accept rv */
  867:                                                         break;
  868:                                                 /* rv = smallest denormal */
  869:                                                 goto undfl;
  870:                                                 }
  871:                                         }
  872: #endif /*Avoid_Underflow*/
  873:                                 L = (word0(&rv) & Exp_mask) - Exp_msk1;
  874: #endif /*Sudden_Underflow}}*/
  875:                                 word0(&rv) = L | Bndry_mask1;
  876:                                 word1(&rv) = 0xffffffff;
  877: #ifdef IBM
  878:                                 goto cont;
  879: #else
  880:                                 break;
  881: #endif
  882:                                 }
  883: #ifndef ROUND_BIASED
  884: #ifdef Avoid_Underflow
  885:                         if (Lsb1) {
  886:                                 if (!(word0(&rv) & Lsb1))
  887:                                         break;
  888:                                 }
  889:                         else if (!(word1(&rv) & Lsb))
  890:                                 break;
  891: #else
  892:                         if (!(word1(&rv) & LSB))
  893:                                 break;
  894: #endif
  895: #endif
  896:                         if (dsign)
  897: #ifdef Avoid_Underflow
  898:                                 dval(&rv) += sulp(&rv, scale);
  899: #else
  900:                                 dval(&rv) += ulp(&rv);
  901: #endif
  902: #ifndef ROUND_BIASED
  903:                         else {
  904: #ifdef Avoid_Underflow
  905:                                 dval(&rv) -= sulp(&rv, scale);
  906: #else
  907:                                 dval(&rv) -= ulp(&rv);
  908: #endif
  909: #ifndef Sudden_Underflow
  910:                                 if (!dval(&rv))
  911:                                         goto undfl;
  912: #endif
  913:                                 }
  914: #ifdef Avoid_Underflow
  915:                         dsign = 1 - dsign;
  916: #endif
  917: #endif
  918:                         break;
  919:                         }
  920:                 if ((aadj = ratio(delta, bs)) <= 2.) {
  921:                         if (dsign)
  922:                                 aadj = dval(&aadj1) = 1.;
  923:                         else if (word1(&rv) || word0(&rv) & Bndry_mask) {
  924: #ifndef Sudden_Underflow
  925:                                 if (word1(&rv) == Tiny1 && !word0(&rv))
  926:                                         goto undfl;
  927: #endif
  928:                                 aadj = 1.;
  929:                                 dval(&aadj1) = -1.;
  930:                                 }
  931:                         else {
  932:                                 /* special case -- power of FLT_RADIX to be */
  933:                                 /* rounded down... */
  934: 
  935:                                 if (aadj < 2./FLT_RADIX)
  936:                                         aadj = 1./FLT_RADIX;
  937:                                 else
  938:                                         aadj *= 0.5;
  939:                                 dval(&aadj1) = -aadj;
  940:                                 }
  941:                         }
  942:                 else {
  943:                         aadj *= 0.5;
  944:                         dval(&aadj1) = dsign ? aadj : -aadj;
  945: #ifdef Check_FLT_ROUNDS
  946:                         switch(Rounding) {
  947:                                 case 2: /* towards +infinity */
  948:                                         dval(&aadj1) -= 0.5;
  949:                                         break;
  950:                                 case 0: /* towards 0 */
  951:                                 case 3: /* towards -infinity */
  952:                                         dval(&aadj1) += 0.5;
  953:                                 }
  954: #else
  955:                         if (Flt_Rounds == 0)
  956:                                 dval(&aadj1) += 0.5;
  957: #endif /*Check_FLT_ROUNDS*/
  958:                         }
  959:                 y = word0(&rv) & Exp_mask;
  960: 
  961:                 /* Check for overflow */
  962: 
  963:                 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
  964:                         dval(&rv0) = dval(&rv);
  965:                         word0(&rv) -= P*Exp_msk1;
  966:                         dval(&adj) = dval(&aadj1) * ulp(&rv);
  967:                         dval(&rv) += dval(&adj);
  968:                         if ((word0(&rv) & Exp_mask) >=
  969:                                         Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
  970:                                 if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
  971:                                         goto ovfl;
  972:                                 word0(&rv) = Big0;
  973:                                 word1(&rv) = Big1;
  974:                                 goto cont;
  975:                                 }
  976:                         else
  977:                                 word0(&rv) += P*Exp_msk1;
  978:                         }
  979:                 else {
  980: #ifdef Avoid_Underflow
  981:                         if (scale && y <= 2*P*Exp_msk1) {
  982:                                 if (aadj <= 0x7fffffff) {
  983:                                         if ((z = aadj) <= 0)
  984:                                                 z = 1;
  985:                                         aadj = z;
  986:                                         dval(&aadj1) = dsign ? aadj : -aadj;
  987:                                         }
  988:                                 word0(&aadj1) += (2*P+1)*Exp_msk1 - y;
  989:                                 }
  990:                         dval(&adj) = dval(&aadj1) * ulp(&rv);
  991:                         dval(&rv) += dval(&adj);
  992: #else
  993: #ifdef Sudden_Underflow
  994:                         if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
  995:                                 dval(&rv0) = dval(&rv);
  996:                                 word0(&rv) += P*Exp_msk1;
  997:                                 dval(&adj) = dval(&aadj1) * ulp(&rv);
  998:                                 dval(&rv) += dval(&adj);
  999: #ifdef IBM
 1000:                                 if ((word0(&rv) & Exp_mask) <  P*Exp_msk1)
 1001: #else
 1002:                                 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1)
 1003: #endif
 1004:                                         {
 1005:                                         if (word0(&rv0) == Tiny0
 1006:                                          && word1(&rv0) == Tiny1)
 1007:                                                 goto undfl;
 1008:                                         word0(&rv) = Tiny0;
 1009:                                         word1(&rv) = Tiny1;
 1010:                                         goto cont;
 1011:                                         }
 1012:                                 else
 1013:                                         word0(&rv) -= P*Exp_msk1;
 1014:                                 }
 1015:                         else {
 1016:                                 dval(&adj) = dval(&aadj1) * ulp(&rv);
 1017:                                 dval(&rv) += dval(&adj);
 1018:                                 }
 1019: #else /*Sudden_Underflow*/
 1020:                         /* Compute dval(&adj) so that the IEEE rounding rules will
 1021:                          * correctly round rv + dval(&adj) in some half-way cases.
 1022:                          * If rv * ulp(&rv) is denormalized (i.e.,
 1023:                          * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
 1024:                          * trouble from bits lost to denormalization;
 1025:                          * example: 1.2e-307 .
 1026:                          */
 1027:                         if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
 1028:                                 dval(&aadj1) = (double)(int)(aadj + 0.5);
 1029:                                 if (!dsign)
 1030:                                         dval(&aadj1) = -dval(&aadj1);
 1031:                                 }
 1032:                         dval(&adj) = dval(&aadj1) * ulp(&rv);
 1033:                         dval(&rv) += adj;
 1034: #endif /*Sudden_Underflow*/
 1035: #endif /*Avoid_Underflow*/
 1036:                         }
 1037:                 z = word0(&rv) & Exp_mask;
 1038: #ifndef SET_INEXACT
 1039: #ifdef Avoid_Underflow
 1040:                 if (!scale)
 1041: #endif
 1042:                 if (y == z) {
 1043:                         /* Can we stop now? */
 1044:                         L = (Long)aadj;
 1045:                         aadj -= L;
 1046:                         /* The tolerances below are conservative. */
 1047:                         if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
 1048:                                 if (aadj < .4999999 || aadj > .5000001)
 1049:                                         break;
 1050:                                 }
 1051:                         else if (aadj < .4999999/FLT_RADIX)
 1052:                                 break;
 1053:                         }
 1054: #endif
 1055:  cont:
 1056:                 Bfree(bb);
 1057:                 Bfree(bd);
 1058:                 Bfree(bs);
 1059:                 Bfree(delta);
 1060:                 }
 1061:         Bfree(bb);
 1062:         Bfree(bd);
 1063:         Bfree(bs);
 1064:         Bfree(bd0);
 1065:         Bfree(delta);
 1066: #ifdef SET_INEXACT
 1067:         if (inexact) {
 1068:                 if (!oldinexact) {
 1069:                         word0(&rv0) = Exp_1 + (70 << Exp_shift);
 1070:                         word1(&rv0) = 0;
 1071:                         dval(&rv0) += 1.;
 1072:                         }
 1073:                 }
 1074:         else if (!oldinexact)
 1075:                 clear_inexact();
 1076: #endif
 1077: #ifdef Avoid_Underflow
 1078:         if (scale) {
 1079:                 word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
 1080:                 word1(&rv0) = 0;
 1081:                 dval(&rv) *= dval(&rv0);
 1082: #ifndef NO_ERRNO
 1083:                 /* try to avoid the bug of testing an 8087 register value */
 1084: #ifdef IEEE_Arith
 1085:                 if (!(word0(&rv) & Exp_mask))
 1086: #else
 1087:                 if (word0(&rv) == 0 && word1(&rv) == 0)
 1088: #endif
 1089:                         /*errno = ERANGE*/;
 1090: #endif
 1091:                 }
 1092: #endif /* Avoid_Underflow */
 1093: #ifdef SET_INEXACT
 1094:         if (inexact && !(word0(&rv) & Exp_mask)) {
 1095:                 /* set underflow bit */
 1096:                 dval(&rv0) = 1e-300;
 1097:                 dval(&rv0) *= dval(&rv0);
 1098:                 }
 1099: #endif
 1100:  ret:
 1101:         if (se)
 1102:                 *se = (char *)s;
 1103:         return sign ? -dval(&rv) : dval(&rv);
 1104:         }
 1105: