contrib/gdtoa/strtodg.c

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