StdLib/LibC/gdtoa/strtodg.c

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