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