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  double d; FPI *fpi; Long *exp; ULong *bits; int exact, rd, *irv;
 176 #else
 177  (double 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(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) double d;
 295 #else
 296 mantbits(double 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
 317 #ifdef KR_headers
 318         (s00, se, fpi, exp, bits)
 319         CONST char *s00; char **se; FPI *fpi; Long *exp; ULong *bits;
 320 #else
 321         (CONST char *s00, char **se, FPI *fpi, Long *exp, ULong *bits)
 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 adj, adj0, rv, tol;
 331         Long L;
 332         ULong *b, *be, y, z;
 333         Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0;
 334
 335         irv = STRTOG_Zero;
 336         denorm = sign = nz0 = nz = 0;
 337         dval(rv) = 0.;
 338         rvb = 0;
 339         nbits = fpi->nbits;
 340         for(s = s00;;s++) switch(*s) {
 341                 case '-':
 342                         sign = 1;
 343                         /* no break */
 344                 case '+':
 345                         if (*++s)
 346                                 goto break2;
 347                         /* no break */
 348                 case 0:
 349                         sign = 0;
 350                         irv = STRTOG_NoNumber;
 351                         s = s00;
 352                         goto ret;
 353                 case '\t':
 354                 case '\n':
 355                 case '\v':
 356                 case '\f':
 357                 case '\r':
 358                 case ' ':
 359                         continue;
 360                 default:
 361                         goto break2;
 362                 }
 363  break2:
 364         if (*s == '0') {
 365 #ifndef NO_HEX_FP
 366                 switch(s[1]) {
 367                   case 'x':
 368                   case 'X':
 369                         irv = gethex(&s, fpi, exp, &rvb, sign);
 370                         if (irv == STRTOG_NoNumber) {
 371                                 s = s00;
 372                                 sign = 0;
 373                                 }
 374                         goto ret;
 375                   }
 376 #endif
 377                 nz0 = 1;
 378                 while(*++s == '0') ;
 379                 if (!*s)
 380                         goto ret;
 381                 }
 382         sudden_underflow = fpi->sudden_underflow;
 383         s0 = s;
 384         y = z = 0;
 385         for(decpt = nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
 386                 if (nd < 9)
 387                         y = 10*y + c - '0';
 388                 else if (nd < 16)
 389                         z = 10*z + c - '0';
 390         nd0 = nd;
 391 #ifdef USE_LOCALE
 392         if (c == *localeconv()->decimal_point)
 393 #else
 394         if (c == '.')
 395 #endif
 396                 {
 397                 decpt = 1;
 398                 c = *++s;
 399                 if (!nd) {
 400                         for(; c == '0'; c = *++s)
 401                                 nz++;
 402                         if (c > '0' && c <= '9') {
 403                                 s0 = s;
 404                                 nf += nz;
 405                                 nz = 0;
 406                                 goto have_dig;
 407                                 }
 408                         goto dig_done;
 409                         }
 410                 for(; c >= '0' && c <= '9'; c = *++s) {
 411  have_dig:
 412                         nz++;
 413                         if (c -= '0') {
 414                                 nf += nz;
 415                                 for(i = 1; i < nz; i++)
 416                                         if (nd++ < 9)
 417                                                 y *= 10;
 418                                         else if (nd <= DBL_DIG + 1)
 419                                                 z *= 10;
 420                                 if (nd++ < 9)
 421                                         y = 10*y + c;
 422                                 else if (nd <= DBL_DIG + 1)
 423                                         z = 10*z + c;
 424                                 nz = 0;
 425                                 }
 426                         }
 427                 }
 428  dig_done:
 429         e = 0;
 430         if (c == 'e' || c == 'E') {
 431                 if (!nd && !nz && !nz0) {
 432                         irv = STRTOG_NoNumber;
 433                         s = s00;
 434                         goto ret;
 435                         }
 436                 s00 = s;
 437                 esign = 0;
 438                 switch(c = *++s) {
 439                         case '-':
 440                                 esign = 1;
 441                         case '+':
 442                                 c = *++s;
 443                         }
 444                 if (c >= '0' && c <= '9') {
 445                         while(c == '0')
 446                                 c = *++s;
 447                         if (c > '0' && c <= '9') {
 448                                 L = c - '0';
 449                                 s1 = s;
 450                                 while((c = *++s) >= '0' && c <= '9')
 451                                         L = 10*L + c - '0';
 452                                 if (s - s1 > 8 || L > 19999)
 453                                         /* Avoid confusion from exponents
 454                                          * so large that e might overflow.
 455                                          */
 456                                         e = 19999; /* safe for 16 bit ints */
 457                                 else
 458                                         e = (int)L;
 459                                 if (esign)
 460                                         e = -e;
 461                                 }
 462                         else
 463                                 e = 0;
 464                         }
 465                 else
 466                         s = s00;
 467                 }
 468         if (!nd) {
 469                 if (!nz && !nz0) {
 470 #ifdef INFNAN_CHECK
 471                         /* Check for Nan and Infinity */
 472                         if (!decpt)
 473                          switch(c) {
 474                           case 'i':
 475                           case 'I':
 476                                 if (match(&s,"nf")) {
 477                                         --s;
 478                                         if (!match(&s,"inity"))
 479                                                 ++s;
 480                                         irv = STRTOG_Infinite;
 481                                         goto infnanexp;
 482                                         }
 483                                 break;
 484                           case 'n':
 485                           case 'N':
 486                                 if (match(&s, "an")) {
 487                                         irv = STRTOG_NaN;
 488                                         *exp = fpi->emax + 1;
 489 #ifndef No_Hex_NaN
 490                                         if (*s == '(') /*)*/
 491                                                 irv = hexnan(&s, fpi, bits);
 492 #endif
 493                                         goto infnanexp;
 494                                         }
 495                           }
 496 #endif /* INFNAN_CHECK */
 497                         irv = STRTOG_NoNumber;
 498                         s = s00;
 499                         }
 500                 goto ret;
 501                 }
 502
 503         irv = STRTOG_Normal;
 504         e1 = e -= nf;
 505         rd = 0;
 506         switch(fpi->rounding & 3) {
 507           case FPI_Round_up:
 508                 rd = 2 - sign;
 509                 break;
 510           case FPI_Round_zero:
 511                 rd = 1;
 512                 break;
 513           case FPI_Round_down:
 514                 rd = 1 + sign;
 515           }
 516
 517         /* Now we have nd0 digits, starting at s0, followed by a
 518          * decimal point, followed by nd-nd0 digits.  The number we're
 519          * after is the integer represented by those digits times
 520          * 10**e */
 521
 522         if (!nd0)
 523                 nd0 = nd;
 524         k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
 525         dval(rv) = y;
 526         if (k > 9)
 527                 dval(rv) = tens[k - 9] * dval(rv) + z;
 528         bd0 = 0;
 529         if (nbits <= P && nd <= DBL_DIG) {
 530                 if (!e) {
 531                         if (rvOK(dval(rv), fpi, exp, bits, 1, rd, &irv))
 532                                 goto ret;
 533                         }
 534                 else if (e > 0) {
 535                         if (e <= Ten_pmax) {
 536 #ifdef VAX
 537                                 goto vax_ovfl_check;
 538 #else
 539                                 i = fivesbits[e] + mantbits(dval(rv)) <= P;
 540                                 /* rv = */ rounded_product(dval(rv), tens[e]);
 541                                 if (rvOK(dval(rv), fpi, exp, bits, i, rd, &irv))
 542                                         goto ret;
 543                                 e1 -= e;
 544                                 goto rv_notOK;
 545 #endif
 546                                 }
 547                         i = DBL_DIG - nd;
 548                         if (e <= Ten_pmax + i) {
 549                                 /* A fancier test would sometimes let us do
 550                                  * this for larger i values.
 551                                  */
 552                                 e2 = e - i;
 553                                 e1 -= i;
 554                                 dval(rv) *= tens[i];
 555 #ifdef VAX
 556                                 /* VAX exponent range is so narrow we must
 557                                  * worry about overflow here...
 558                                  */
 559  vax_ovfl_check:
 560                                 dval(adj) = dval(rv);
 561                                 word0(adj) -= P*Exp_msk1;
 562                                 /* adj = */ rounded_product(dval(adj), tens[e2]);
 563                                 if ((word0(adj) & Exp_mask)
 564                                  > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
 565                                         goto rv_notOK;
 566                                 word0(adj) += P*Exp_msk1;
 567                                 dval(rv) = dval(adj);
 568 #else
 569                                 /* rv = */ rounded_product(dval(rv), tens[e2]);
 570 #endif
 571                                 if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv))
 572                                         goto ret;
 573                                 e1 -= e2;
 574                                 }
 575                         }
 576 #ifndef Inaccurate_Divide
 577                 else if (e >= -Ten_pmax) {
 578                         /* rv = */ rounded_quotient(dval(rv), tens[-e]);
 579                         if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv))
 580                                 goto ret;
 581                         e1 -= e;
 582                         }
 583 #endif
 584                 }
 585  rv_notOK:
 586         e1 += nd - k;
 587
 588         /* Get starting approximation = rv * 10**e1 */
 589
 590         e2 = 0;
 591         if (e1 > 0) {
 592                 if ( (i = e1 & 15) !=0)
 593                         dval(rv) *= tens[i];
 594                 if (e1 &= ~15) {
 595                         e1 >>= 4;
 596                         while(e1 >= (1 << n_bigtens-1)) {
 597                                 e2 += ((word0(rv) & Exp_mask)
 598                                         >> Exp_shift1) - Bias;
 599                                 word0(rv) &= ~Exp_mask;
 600                                 word0(rv) |= Bias << Exp_shift1;
 601                                 dval(rv) *= bigtens[n_bigtens-1];
 602                                 e1 -= 1 << n_bigtens-1;
 603                                 }
 604                         e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias;
 605                         word0(rv) &= ~Exp_mask;
 606                         word0(rv) |= Bias << Exp_shift1;
 607                         for(j = 0; e1 > 0; j++, e1 >>= 1)
 608                                 if (e1 & 1)
 609                                         dval(rv) *= bigtens[j];
 610                         }
 611                 }
 612         else if (e1 < 0) {
 613                 e1 = -e1;
 614                 if ( (i = e1 & 15) !=0)
 615                         dval(rv) /= tens[i];
 616                 if (e1 &= ~15) {
 617                         e1 >>= 4;
 618                         while(e1 >= (1 << n_bigtens-1)) {
 619                                 e2 += ((word0(rv) & Exp_mask)
 620                                         >> Exp_shift1) - Bias;
 621                                 word0(rv) &= ~Exp_mask;
 622                                 word0(rv) |= Bias << Exp_shift1;
 623                                 dval(rv) *= tinytens[n_bigtens-1];
 624                                 e1 -= 1 << n_bigtens-1;
 625                                 }
 626                         e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias;
 627                         word0(rv) &= ~Exp_mask;
 628                         word0(rv) |= Bias << Exp_shift1;
 629                         for(j = 0; e1 > 0; j++, e1 >>= 1)
 630                                 if (e1 & 1)
 631                                         dval(rv) *= tinytens[j];
 632                         }
 633                 }
 634 #ifdef IBM
 635         /* e2 is a correction to the (base 2) exponent of the return
 636          * value, reflecting adjustments above to avoid overflow in the
 637          * native arithmetic.  For native IBM (base 16) arithmetic, we
 638          * must multiply e2 by 4 to change from base 16 to 2.
 639          */
 640         e2 <<= 2;
 641 #endif
 642         rvb = d2b(dval(rv), &rve, &rvbits);     /* rv = rvb * 2^rve */
 643         rve += e2;
 644         if ((j = rvbits - nbits) > 0) {
 645                 rshift(rvb, j);
 646                 rvbits = nbits;
 647                 rve += j;
 648                 }
 649         bb0 = 0;        /* trailing zero bits in rvb */
 650         e2 = rve + rvbits - nbits;
 651         if (e2 > fpi->emax + 1)
 652                 goto huge;
 653         rve1 = rve + rvbits - nbits;
 654         if (e2 < (emin = fpi->emin)) {
 655                 denorm = 1;
 656                 j = rve - emin;
 657                 if (j > 0) {
 658                         rvb = lshift(rvb, j);
 659                         rvbits += j;
 660                         }
 661                 else if (j < 0) {
 662                         rvbits += j;
 663                         if (rvbits <= 0) {
 664                                 if (rvbits < -1) {
 665  ufl:
 666                                         rvb->wds = 0;
 667                                         rvb->x[0] = 0;
 668                                         *exp = emin;
 669                                         irv = STRTOG_Underflow | STRTOG_Inexlo;
 670                                         goto ret;
 671                                         }
 672                                 rvb->x[0] = rvb->wds = rvbits = 1;
 673                                 }
 674                         else
 675                                 rshift(rvb, -j);
 676                         }
 677                 rve = rve1 = emin;
 678                 if (sudden_underflow && e2 + 1 < emin)
 679                         goto ufl;
 680                 }
 681
 682         /* Now the hard part -- adjusting rv to the correct value.*/
 683
 684         /* Put digits into bd: true value = bd * 10^e */
 685
 686         bd0 = s2b(s0, nd0, nd, y);
 687
 688         for(;;) {
 689                 bd = Balloc(bd0->k);
 690                 Bcopy(bd, bd0);
 691                 bb = Balloc(rvb->k);
 692                 Bcopy(bb, rvb);
 693                 bbbits = rvbits - bb0;
 694                 bbe = rve + bb0;
 695                 bs = i2b(1);
 696
 697                 if (e >= 0) {
 698                         bb2 = bb5 = 0;
 699                         bd2 = bd5 = e;
 700                         }
 701                 else {
 702                         bb2 = bb5 = -e;
 703                         bd2 = bd5 = 0;
 704                         }
 705                 if (bbe >= 0)
 706                         bb2 += bbe;
 707                 else
 708                         bd2 -= bbe;
 709                 bs2 = bb2;
 710                 j = nbits + 1 - bbbits;
 711                 i = bbe + bbbits - nbits;
 712                 if (i < emin)   /* denormal */
 713                         j += i - emin;
 714                 bb2 += j;
 715                 bd2 += j;
 716                 i = bb2 < bd2 ? bb2 : bd2;
 717                 if (i > bs2)
 718                         i = bs2;
 719                 if (i > 0) {
 720                         bb2 -= i;
 721                         bd2 -= i;
 722                         bs2 -= i;
 723                         }
 724                 if (bb5 > 0) {
 725                         bs = pow5mult(bs, bb5);
 726                         bb1 = mult(bs, bb);
 727                         Bfree(bb);
 728                         bb = bb1;
 729                         }
 730                 bb2 -= bb0;
 731                 if (bb2 > 0)
 732                         bb = lshift(bb, bb2);
 733                 else if (bb2 < 0)
 734                         rshift(bb, -bb2);
 735                 if (bd5 > 0)
 736                         bd = pow5mult(bd, bd5);
 737                 if (bd2 > 0)
 738                         bd = lshift(bd, bd2);
 739                 if (bs2 > 0)
 740                         bs = lshift(bs, bs2);
 741                 asub = 1;
 742                 inex = STRTOG_Inexhi;
 743                 delta = diff(bb, bd);
 744                 if (delta->wds <= 1 && !delta->x[0])
 745                         break;
 746                 dsign = delta->sign;
 747                 delta->sign = finished = 0;
 748                 L = 0;
 749                 i = cmp(delta, bs);
 750                 if (rd && i <= 0) {
 751                         irv = STRTOG_Normal;
 752                         if ( (finished = dsign ^ (rd&1)) !=0) {
 753                                 if (dsign != 0) {
 754                                         irv |= STRTOG_Inexhi;
 755                                         goto adj1;
 756                                         }
 757                                 irv |= STRTOG_Inexlo;
 758                                 if (rve1 == emin)
 759                                         goto adj1;
 760                                 for(i = 0, j = nbits; j >= ULbits;
 761                                                 i++, j -= ULbits) {
 762                                         if (rvb->x[i] & ALL_ON)
 763                                                 goto adj1;
 764                                         }
 765                                 if (j > 1 && lo0bits(rvb->x + i) < j - 1)
 766                                         goto adj1;
 767                                 rve = rve1 - 1;
 768                                 rvb = set_ones(rvb, rvbits = nbits);
 769                                 break;
 770                                 }
 771                         irv |= dsign ? STRTOG_Inexlo : STRTOG_Inexhi;
 772                         break;
 773                         }
 774                 if (i < 0) {
 775                         /* Error is less than half an ulp -- check for
 776                          * special case of mantissa a power of two.
 777                          */
 778                         irv = dsign
 779                                 ? STRTOG_Normal | STRTOG_Inexlo
 780                                 : STRTOG_Normal | STRTOG_Inexhi;
 781                         if (dsign || bbbits > 1 || denorm || rve1 == emin)
 782                                 break;
 783                         delta = lshift(delta,1);
 784                         if (cmp(delta, bs) > 0) {
 785                                 irv = STRTOG_Normal | STRTOG_Inexlo;
 786                                 goto drop_down;
 787                                 }
 788                         break;
 789                         }
 790                 if (i == 0) {
 791                         /* exactly half-way between */
 792                         if (dsign) {
 793                                 if (denorm && all_on(rvb, rvbits)) {
 794                                         /*boundary case -- increment exponent*/
 795                                         rvb->wds = 1;
 796                                         rvb->x[0] = 1;
 797                                         rve = emin + nbits - (rvbits = 1);
 798                                         irv = STRTOG_Normal | STRTOG_Inexhi;
 799                                         denorm = 0;
 800                                         break;
 801                                         }
 802                                 irv = STRTOG_Normal | STRTOG_Inexlo;
 803                                 }
 804                         else if (bbbits == 1) {
 805                                 irv = STRTOG_Normal;
 806  drop_down:
 807                                 /* boundary case -- decrement exponent */
 808                                 if (rve1 == emin) {
 809                                         irv = STRTOG_Normal | STRTOG_Inexhi;
 810                                         if (rvb->wds == 1 && rvb->x[0] == 1)
 811                                                 sudden_underflow = 1;
 812                                         break;
 813                                         }
 814                                 rve -= nbits;
 815                                 rvb = set_ones(rvb, rvbits = nbits);
 816                                 break;
 817                                 }
 818                         else
 819                                 irv = STRTOG_Normal | STRTOG_Inexhi;
 820                         if (bbbits < nbits && !denorm || !(rvb->x[0] & 1))
 821                                 break;
 822                         if (dsign) {
 823                                 rvb = increment(rvb);
 824                                 j = kmask & (ULbits - (rvbits & kmask));
 825                                 if (hi0bits(rvb->x[rvb->wds - 1]) != j)
 826                                         rvbits++;
 827                                 irv = STRTOG_Normal | STRTOG_Inexhi;
 828                                 }
 829                         else {
 830                                 if (bbbits == 1)
 831                                         goto undfl;
 832                                 decrement(rvb);
 833                                 irv = STRTOG_Normal | STRTOG_Inexlo;
 834                                 }
 835                         break;
 836                         }
 837                 if ((dval(adj) = ratio(delta, bs)) <= 2.) {
 838  adj1:
 839                         inex = STRTOG_Inexlo;
 840                         if (dsign) {
 841                                 asub = 0;
 842                                 inex = STRTOG_Inexhi;
 843                                 }
 844                         else if (denorm && bbbits <= 1) {
 845  undfl:
 846                                 rvb->wds = 0;
 847                                 rve = emin;
 848                                 irv = STRTOG_Underflow | STRTOG_Inexlo;
 849                                 break;
 850                                 }
 851                         adj0 = dval(adj) = 1.;
 852                         }
 853                 else {
 854                         adj0 = dval(adj) *= 0.5;
 855                         if (dsign) {
 856                                 asub = 0;
 857                                 inex = STRTOG_Inexlo;
 858                                 }
 859                         if (dval(adj) < 2147483647.) {
 860                                 L = adj0;
 861                                 adj0 -= L;
 862                                 switch(rd) {
 863                                   case 0:
 864                                         if (adj0 >= .5)
 865                                                 goto inc_L;
 866                                         break;
 867                                   case 1:
 868                                         if (asub && adj0 > 0.)
 869                                                 goto inc_L;
 870                                         break;
 871                                   case 2:
 872                                         if (!asub && adj0 > 0.) {
 873  inc_L:
 874                                                 L++;
 875                                                 inex = STRTOG_Inexact - inex;
 876                                                 }
 877                                   }
 878                                 dval(adj) = L;
 879                                 }
 880                         }
 881                 y = rve + rvbits;
 882
 883                 /* adj *= ulp(dval(rv)); */
 884                 /* if (asub) rv -= adj; else rv += adj; */
 885
 886                 if (!denorm && rvbits < nbits) {
 887                         rvb = lshift(rvb, j = nbits - rvbits);
 888                         rve -= j;
 889                         rvbits = nbits;
 890                         }
 891                 ab = d2b(dval(adj), &abe, &abits);
 892                 if (abe < 0)
 893                         rshift(ab, -abe);
 894                 else if (abe > 0)
 895                         ab = lshift(ab, abe);
 896                 rvb0 = rvb;
 897                 if (asub) {
 898                         /* rv -= adj; */
 899                         j = hi0bits(rvb->x[rvb->wds-1]);
 900                         rvb = diff(rvb, ab);
 901                         k = rvb0->wds - 1;
 902                         if (denorm)
 903                                 /* do nothing */;
 904                         else if (rvb->wds <= k
 905                                 || hi0bits( rvb->x[k]) >
 906                                    hi0bits(rvb0->x[k])) {
 907                                 /* unlikely; can only have lost 1 high bit */
 908                                 if (rve1 == emin) {
 909                                         --rvbits;
 910                                         denorm = 1;
 911                                         }
 912                                 else {
 913                                         rvb = lshift(rvb, 1);
 914                                         --rve;
 915                                         --rve1;
 916                                         L = finished = 0;
 917                                         }
 918                                 }
 919                         }
 920                 else {
 921                         rvb = sum(rvb, ab);
 922                         k = rvb->wds - 1;
 923                         if (k >= rvb0->wds
 924                          || hi0bits(rvb->x[k]) < hi0bits(rvb0->x[k])) {
 925                                 if (denorm) {
 926                                         if (++rvbits == nbits)
 927                                                 denorm = 0;
 928                                         }
 929                                 else {
 930                                         rshift(rvb, 1);
 931                                         rve++;
 932                                         rve1++;
 933                                         L = 0;
 934                                         }
 935                                 }
 936                         }
 937                 Bfree(ab);
 938                 Bfree(rvb0);
 939                 if (finished)
 940                         break;
 941
 942                 z = rve + rvbits;
 943                 if (y == z && L) {
 944                         /* Can we stop now? */
 945                         tol = dval(adj) * 5e-16; /* > max rel error */
 946                         dval(adj) = adj0 - .5;
 947                         if (dval(adj) < -tol) {
 948                                 if (adj0 > tol) {
 949                                         irv |= inex;
 950                                         break;
 951                                         }
 952                                 }
 953                         else if (dval(adj) > tol && adj0 < 1. - tol) {
 954                                 irv |= inex;
 955                                 break;
 956                                 }
 957                         }
 958                 bb0 = denorm ? 0 : trailz(rvb);
 959                 Bfree(bb);
 960                 Bfree(bd);
 961                 Bfree(bs);
 962                 Bfree(delta);
 963                 }
 964         if (!denorm && (j = nbits - rvbits)) {
 965                 if (j > 0)
 966                         rvb = lshift(rvb, j);
 967                 else
 968                         rshift(rvb, -j);
 969                 rve -= j;
 970                 }
 971         *exp = rve;
 972         Bfree(bb);
 973         Bfree(bd);
 974         Bfree(bs);
 975         Bfree(bd0);
 976         Bfree(delta);
 977         if (rve > fpi->emax) {
 978                 switch(fpi->rounding & 3) {
 979                   case FPI_Round_near:
 980                         goto huge;
 981                   case FPI_Round_up:
 982                         if (!sign)
 983                                 goto huge;
 984                         break;
 985                   case FPI_Round_down:
 986                         if (sign)
 987                                 goto huge;
 988                   }
 989                 /* Round to largest representable magnitude */
 990                 Bfree(rvb);
 991                 rvb = 0;
 992                 irv = STRTOG_Normal | STRTOG_Inexlo;
 993                 *exp = fpi->emax;
 994                 b = bits;
 995                 be = b + (fpi->nbits >> 5) + 1;
 996                 while(b < be)
 997                         *b++ = -1;
 998                 if ((j = fpi->nbits & 0x1f))
 999                         *--be >>= (32 - j);
1000                 goto ret;
1001  huge:
1002                 rvb->wds = 0;
1003                 irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
1004 #ifndef NO_ERRNO
1005                 errno = ERANGE;
1006 #endif
1007  infnanexp:
1008                 *exp = fpi->emax + 1;
1009                 }
1010  ret:
1011         if (denorm) {
1012                 if (sudden_underflow) {
1013                         rvb->wds = 0;
1014                         irv = STRTOG_Underflow | STRTOG_Inexlo;
1015 #ifndef NO_ERRNO
1016                         errno = ERANGE;
1017 #endif
1018                         }
1019                 else  {
1020                         irv = (irv & ~STRTOG_Retmask) |
1021                                 (rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero);
1022                         if (irv & STRTOG_Inexact) {
1023                                 irv |= STRTOG_Underflow;
1024 #ifndef NO_ERRNO
1025                                 errno = ERANGE;
1026 #endif
1027                                 }
1028                         }
1029                 }
1030         if (se)
1031                 *se = (char *)s;
1032         if (sign)
1033                 irv |= STRTOG_Neg;
1034         if (rvb) {
1035                 copybits(bits, nbits, rvb);
1036                 Bfree(rvb);
1037                 }
1038         return irv;
1039         }