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