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
  30         David M. Gay
  31         dmg@acm.org
  32  */
  33
  34 #include "gdtoaimp.h"
  35
  36 #ifdef USE_LOCALE
  37 #include "locale.h"
  38 #endif
  39
  40  static CONST int
  41 fivesbits[] = {  0,  3,  5,  7, 10, 12, 14, 17, 19, 21,
  42                 24, 26, 28, 31, 33, 35, 38, 40, 42, 45,
  43                 47, 49, 52
  44 #ifdef VAX
  45                 , 54, 56
  46 #endif
  47                 };
  48
  49  Bigint *
  50 #ifdef KR_headers
  51 increment(b) Bigint *b;
  52 #else
  53 increment(Bigint *b)
  54 #endif
  55 {
  56         ULong *x, *xe;
  57         Bigint *b1;
  58 #ifdef Pack_16
  59         ULong carry = 1, y;
  60 #endif
  61
  62         x = b->x;
  63         xe = x + b->wds;
  64 #ifdef Pack_32
  65         do {
  66                 if (*x < (ULong)0xffffffffL) {
  67                         ++*x;
  68                         return b;
  69                         }
  70                 *x++ = 0;
  71                 } while(x < xe);
  72 #else
  73         do {
  74                 y = *x + carry;
  75                 carry = y >> 16;
  76                 *x++ = y & 0xffff;
  77                 if (!carry)
  78                         return b;
  79                 } while(x < xe);
  80         if (carry)
  81 #endif
  82         {
  83                 if (b->wds >= b->maxwds) {
  84                         b1 = Balloc(b->k+1);
  85                         Bcopy(b1,b);
  86                         Bfree(b);
  87                         b = b1;
  88                         }
  89                 b->x[b->wds++] = 1;
  90                 }
  91         return b;
  92         }
  93
  94  int
  95 #ifdef KR_headers
  96 decrement(b) Bigint *b;
  97 #else
  98 decrement(Bigint *b)
  99 #endif
 100 {
 101         ULong *x, *xe;
 102 #ifdef Pack_16
 103         ULong borrow = 1, y;
 104 #endif
 105
 106         x = b->x;
 107         xe = x + b->wds;
 108 #ifdef Pack_32
 109         do {
 110                 if (*x) {
 111                         --*x;
 112                         break;
 113                         }
 114                 *x++ = 0xffffffffL;
 115                 }
 116                 while(x < xe);
 117 #else
 118         do {
 119                 y = *x - borrow;
 120                 borrow = (y & 0x10000) >> 16;
 121                 *x++ = y & 0xffff;
 122                 } while(borrow && x < xe);
 123 #endif
 124         return STRTOG_Inexlo;
 125         }
 126
 127  static int
 128 #ifdef KR_headers
 129 all_on(b, n) Bigint *b; int n;
 130 #else
 131 all_on(Bigint *b, int n)
 132 #endif
 133 {
 134         ULong *x, *xe;
 135
 136         x = b->x;
 137         xe = x + (n >> kshift);
 138         while(x < xe)
 139                 if ((*x++ & ALL_ON) != ALL_ON)
 140                         return 0;
 141         if (n &= kmask)
 142                 return ((*x | (ALL_ON << n)) & ALL_ON) == ALL_ON;
 143         return 1;
 144         }
 145
 146  Bigint *
 147 #ifdef KR_headers
 148 set_ones(b, n) Bigint *b; int n;
 149 #else
 150 set_ones(Bigint *b, int n)
 151 #endif
 152 {
 153         int k;
 154         ULong *x, *xe;
 155
 156         k = (n + ((1 << kshift) - 1)) >> kshift;
 157         if (b->k < k) {
 158                 Bfree(b);
 159                 b = Balloc(k);
 160                 }
 161         k = n >> kshift;
 162         if (n &= kmask)
 163                 k++;
 164         b->wds = k;
 165         x = b->x;
 166         xe = x + k;
 167         while(x < xe)
 168                 *x++ = ALL_ON;
 169         if (n)
 170                 x[-1] >>= ULbits - n;
 171         return b;
 172         }
 173
 174  static int
 175 rvOK
 176 #ifdef KR_headers
 177  (d, fpi, exp, bits, exact, rd, irv)
 178  double d; FPI *fpi; Long *exp; ULong *bits; int exact, rd, *irv;
 179 #else
 180  (double d, FPI *fpi, Long *exp, ULong *bits, int exact, int rd, int *irv)
 181 #endif
 182 {
 183         Bigint *b;
 184         ULong carry, inex, lostbits;
 185         int bdif, e, j, k, k1, nb, rv;
 186
 187         carry = rv = 0;
 188         b = d2b(d, &e, &bdif);
 189         bdif -= nb = fpi->nbits;
 190         e += bdif;
 191         if (bdif <= 0) {
 192                 if (exact)
 193                         goto trunc;
 194                 goto ret;
 195                 }
 196         if (P == nb) {
 197                 if (
 198 #ifndef IMPRECISE_INEXACT
 199                         exact &&
 200 #endif
 201                         fpi->rounding ==
 202 #ifdef RND_PRODQUOT
 203                                         FPI_Round_near
 204 #else
 205                                         Flt_Rounds
 206 #endif
 207                         ) goto trunc;
 208                 goto ret;
 209                 }
 210         switch(rd) {
 211           case 1:
 212                 goto trunc;
 213           case 2:
 214                 break;
 215           default: /* round near */
 216                 k = bdif - 1;
 217                 if (k < 0)
 218                         goto trunc;
 219                 if (!k) {
 220                         if (!exact)
 221                                 goto ret;
 222                         if (b->x[0] & 2)
 223                                 break;
 224                         goto trunc;
 225                         }
 226                 if (b->x[k>>kshift] & ((ULong)1 << (k & kmask)))
 227                         break;
 228                 goto trunc;
 229           }
 230         /* "break" cases: round up 1 bit, then truncate; bdif > 0 */
 231         carry = 1;
 232  trunc:
 233         inex = lostbits = 0;
 234         if (bdif > 0) {
 235                 if ( (lostbits = any_on(b, bdif)) !=0)
 236                         inex = STRTOG_Inexlo;
 237                 rshift(b, bdif);
 238                 if (carry) {
 239                         inex = STRTOG_Inexhi;
 240                         b = increment(b);
 241                         if ( (j = nb & kmask) !=0)
 242                                 j = 32 - j;
 243                         if (hi0bits(b->x[b->wds - 1]) != j) {
 244                                 if (!lostbits)
 245                                         lostbits = b->x[0] & 1;
 246                                 rshift(b, 1);
 247                                 e++;
 248                                 }
 249                         }
 250                 }
 251         else if (bdif < 0)
 252                 b = lshift(b, -bdif);
 253         if (e < fpi->emin) {
 254                 k = fpi->emin - e;
 255                 e = fpi->emin;
 256                 if (k > nb || fpi->sudden_underflow) {
 257                         b->wds = inex = 0;
 258                         *irv = STRTOG_Underflow | STRTOG_Inexlo;
 259                         }
 260                 else {
 261                         k1 = k - 1;
 262                         if (k1 > 0 && !lostbits)
 263                                 lostbits = any_on(b, k1);
 264                         if (!lostbits && !exact)
 265                                 goto ret;
 266                         lostbits |=
 267                           carry = b->x[k1>>kshift] & (1 << (k1 & kmask));
 268                         rshift(b, k);
 269                         *irv = STRTOG_Denormal;
 270                         if (carry) {
 271                                 b = increment(b);
 272                                 inex = STRTOG_Inexhi | STRTOG_Underflow;
 273                                 }
 274                         else if (lostbits)
 275                                 inex = STRTOG_Inexlo | STRTOG_Underflow;
 276                         }
 277                 }
 278         else if (e > fpi->emax) {
 279                 e = fpi->emax + 1;
 280                 *irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
 281 #ifndef NO_ERRNO
 282                 errno = ERANGE;
 283 #endif
 284                 b->wds = inex = 0;
 285                 }
 286         *exp = e;
 287         copybits(bits, nb, b);
 288         *irv |= inex;
 289         rv = 1;
 290  ret:
 291         Bfree(b);
 292         return rv;
 293         }
 294
 295  static int
 296 #ifdef KR_headers
 297 mantbits(d) double d;
 298 #else
 299 mantbits(double d)
 300 #endif
 301 {
 302         ULong L;
 303 #ifdef VAX
 304         L = word1(d) << 16 | word1(d) >> 16;
 305         if (L)
 306 #else
 307         if ( (L = word1(d)) !=0)
 308 #endif
 309                 return P - lo0bits(&L);
 310 #ifdef VAX
 311         L = word0(d) << 16 | word0(d) >> 16 | Exp_msk11;
 312 #else
 313         L = word0(d) | Exp_msk1;
 314 #endif
 315         return P - 32 - lo0bits(&L);
 316         }
 317
 318  int
 319 strtodg
 320 #ifdef KR_headers
 321         (s00, se, fpi, exp, bits)
 322         CONST char *s00; char **se; FPI *fpi; Long *exp; ULong *bits;
 323 #else
 324         (CONST char *s00, char **se, FPI *fpi, Long *exp, ULong *bits)
 325 #endif
 326 {
 327         int abe, abits, asub;
 328         int bb0, bb2, bb5, bbe, bd2, bd5, bbbits, bs2;
 329         int c, denorm, dsign, e, e1, e2, emin, esign, finished, i, inex, irv;
 330         int j, k, nbits, nd, nd0, nf, nz, nz0, rd, rvbits, rve, rve1, sign;
 331         int sudden_underflow;
 332         CONST char *s, *s0, *s1;
 333         double adj, adj0, rv, tol;
 334         Long L;
 335         ULong y, z;
 336         Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0;
 337
 338         irv = STRTOG_Zero;
 339         denorm = sign = nz0 = nz = 0;
 340         dval(rv) = 0.;
 341         rvb = 0;
 342         nbits = fpi->nbits;
 343         for(s = s00;;s++) switch(*s) {
 344                 case '-':
 345                         sign = 1;
 346                         /* no break */
 347                 case '+':
 348                         if (*++s)
 349                                 goto break2;
 350                         /* no break */
 351                 case 0:
 352                         sign = 0;
 353                         irv = STRTOG_NoNumber;
 354                         s = s00;
 355                         goto ret;
 356                 case '\t':
 357                 case '\n':
 358                 case '\v':
 359                 case '\f':
 360                 case '\r':
 361                 case ' ':
 362                         continue;
 363                 default:
 364                         goto break2;
 365                 }
 366  break2:
 367         if (*s == '0') {
 368 #ifndef NO_HEX_FP
 369                 switch(s[1]) {
 370                   case 'x':
 371                   case 'X':
 372                         irv = gethex(&s, fpi, exp, &rvb, sign);
 373                         if (irv == STRTOG_NoNumber) {
 374                                 s = s00;
 375                                 sign = 0;
 376                                 }
 377                         goto ret;
 378                   }
 379 #endif
 380                 nz0 = 1;
 381                 while(*++s == '0') ;
 382                 if (!*s)
 383                         goto ret;
 384                 }
 385         sudden_underflow = fpi->sudden_underflow;
 386         s0 = s;
 387         y = z = 0;
 388         for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
 389                 if (nd < 9)
 390                         y = 10*y + c - '0';
 391                 else if (nd < 16)
 392                         z = 10*z + c - '0';
 393         nd0 = nd;
 394 #ifdef USE_LOCALE
 395         if (c == *localeconv()->decimal_point)
 396 #else
 397         if (c == '.')
 398 #endif
 399                 {
 400                 c = *++s;
 401                 if (!nd) {
 402                         for(; c == '0'; c = *++s)
 403                                 nz++;
 404                         if (c > '0' && c <= '9') {
 405                                 s0 = s;
 406                                 nf += nz;
 407                                 nz = 0;
 408                                 goto have_dig;
 409                                 }
 410                         goto dig_done;
 411                         }
 412                 for(; c >= '0' && c <= '9'; c = *++s) {
 413  have_dig:
 414                         nz++;
 415                         if (c -= '0') {
 416                                 nf += nz;
 417                                 for(i = 1; i < nz; i++)
 418                                         if (nd++ < 9)
 419                                                 y *= 10;
 420                                         else if (nd <= DBL_DIG + 1)
 421                                                 z *= 10;
 422                                 if (nd++ < 9)
 423                                         y = 10*y + c;
 424                                 else if (nd <= DBL_DIG + 1)
 425                                         z = 10*z + c;
 426                                 nz = 0;
 427                                 }
 428                         }
 429                 }
 430  dig_done:
 431         e = 0;
 432         if (c == 'e' || c == 'E') {
 433                 if (!nd && !nz && !nz0) {
 434                         irv = STRTOG_NoNumber;
 435                         s = s00;
 436                         goto ret;
 437                         }
 438                 s00 = s;
 439                 esign = 0;
 440                 switch(c = *++s) {
 441                         case '-':
 442                                 esign = 1;
 443                         case '+':
 444                                 c = *++s;
 445                         }
 446                 if (c >= '0' && c <= '9') {
 447                         while(c == '0')
 448                                 c = *++s;
 449                         if (c > '0' && c <= '9') {
 450                                 L = c - '0';
 451                                 s1 = s;
 452                                 while((c = *++s) >= '0' && c <= '9')
 453                                         L = 10*L + c - '0';
 454                                 if (s - s1 > 8 || L > 19999)
 455                                         /* Avoid confusion from exponents
 456                                          * so large that e might overflow.
 457                                          */
 458                                         e = 19999; /* safe for 16 bit ints */
 459                                 else
 460                                         e = (int)L;
 461                                 if (esign)
 462                                         e = -e;
 463                                 }
 464                         else
 465                                 e = 0;
 466                         }
 467                 else
 468                         s = s00;
 469                 }
 470         if (!nd) {
 471                 if (!nz && !nz0) {
 472 #ifdef INFNAN_CHECK
 473                         /* Check for Nan and Infinity */
 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
 636         rvb = d2b(dval(rv), &rve, &rvbits);     /* rv = rvb * 2^rve */
 637         rve += e2;
 638         if ((j = rvbits - nbits) > 0) {
 639                 rshift(rvb, j);
 640                 rvbits = nbits;
 641                 rve += j;
 642                 }
 643         bb0 = 0;        /* trailing zero bits in rvb */
 644         e2 = rve + rvbits - nbits;
 645         if (e2 > fpi->emax) {
 646                 rvb->wds = 0;
 647                 irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
 648 #ifndef NO_ERRNO
 649                 errno = ERANGE;
 650 #endif
 651  infnanexp:
 652                 *exp = fpi->emax + 1;
 653                 goto ret;
 654                 }
 655         rve1 = rve + rvbits - nbits;
 656         if (e2 < (emin = fpi->emin)) {
 657                 denorm = 1;
 658                 j = rve - emin;
 659                 if (j > 0) {
 660                         rvb = lshift(rvb, j);
 661                         rvbits += j;
 662                         }
 663                 else if (j < 0) {
 664                         rvbits += j;
 665                         if (rvbits <= 0) {
 666                                 if (rvbits < -1) {
 667  ufl:
 668                                         rvb->wds = 0;
 669                                         rvb->x[0] = 0;
 670                                         *exp = emin;
 671                                         irv = STRTOG_Underflow | STRTOG_Inexlo;
 672                                         goto ret;
 673                                         }
 674                                 rvb->x[0] = rvb->wds = rvbits = 1;
 675                                 }
 676                         else
 677                                 rshift(rvb, -j);
 678                         }
 679                 rve = rve1 = emin;
 680                 if (sudden_underflow && e2 + 1 < emin)
 681                         goto ufl;
 682                 }
 683
 684         /* Now the hard part -- adjusting rv to the correct value.*/
 685
 686         /* Put digits into bd: true value = bd * 10^e */
 687
 688         bd0 = s2b(s0, nd0, nd, y);
 689
 690         for(;;) {
 691                 bd = Balloc(bd0->k);
 692                 Bcopy(bd, bd0);
 693                 bb = Balloc(rvb->k);
 694                 Bcopy(bb, rvb);
 695                 bbbits = rvbits - bb0;
 696                 bbe = rve + bb0;
 697                 bs = i2b(1);
 698
 699                 if (e >= 0) {
 700                         bb2 = bb5 = 0;
 701                         bd2 = bd5 = e;
 702                         }
 703                 else {
 704                         bb2 = bb5 = -e;
 705                         bd2 = bd5 = 0;
 706                         }
 707                 if (bbe >= 0)
 708                         bb2 += bbe;
 709                 else
 710                         bd2 -= bbe;
 711                 bs2 = bb2;
 712                 j = nbits + 1 - bbbits;
 713                 i = bbe + bbbits - nbits;
 714                 if (i < emin)   /* denormal */
 715                         j += i - emin;
 716                 bb2 += j;
 717                 bd2 += j;
 718                 i = bb2 < bd2 ? bb2 : bd2;
 719                 if (i > bs2)
 720                         i = bs2;
 721                 if (i > 0) {
 722                         bb2 -= i;
 723                         bd2 -= i;
 724                         bs2 -= i;
 725                         }
 726                 if (bb5 > 0) {
 727                         bs = pow5mult(bs, bb5);
 728                         bb1 = mult(bs, bb);
 729                         Bfree(bb);
 730                         bb = bb1;
 731                         }
 732                 bb2 -= bb0;
 733                 if (bb2 > 0)
 734                         bb = lshift(bb, bb2);
 735                 else if (bb2 < 0)
 736                         rshift(bb, -bb2);
 737                 if (bd5 > 0)
 738                         bd = pow5mult(bd, bd5);
 739                 if (bd2 > 0)
 740                         bd = lshift(bd, bd2);
 741                 if (bs2 > 0)
 742                         bs = lshift(bs, bs2);
 743                 asub = 1;
 744                 inex = STRTOG_Inexhi;
 745                 delta = diff(bb, bd);
 746                 if (delta->wds <= 1 && !delta->x[0])
 747                         break;
 748                 dsign = delta->sign;
 749                 delta->sign = finished = 0;
 750                 L = 0;
 751                 i = cmp(delta, bs);
 752                 if (rd && i <= 0) {
 753                         irv = STRTOG_Normal;
 754                         if ( (finished = dsign ^ (rd&1)) !=0) {
 755                                 if (dsign != 0) {
 756                                         irv |= STRTOG_Inexhi;
 757                                         goto adj1;
 758                                         }
 759                                 irv |= STRTOG_Inexlo;
 760                                 if (rve1 == emin)
 761                                         goto adj1;
 762                                 for(i = 0, j = nbits; j >= ULbits;
 763                                                 i++, j -= ULbits) {
 764                                         if (rvb->x[i] & ALL_ON)
 765                                                 goto adj1;
 766                                         }
 767                                 if (j > 1 && lo0bits(rvb->x + i) < j - 1)
 768                                         goto adj1;
 769                                 rve = rve1 - 1;
 770                                 rvb = set_ones(rvb, rvbits = nbits);
 771                                 break;
 772                                 }
 773                         irv |= dsign ? STRTOG_Inexlo : STRTOG_Inexhi;
 774                         break;
 775                         }
 776                 if (i < 0) {
 777                         /* Error is less than half an ulp -- check for
 778                          * special case of mantissa a power of two.
 779                          */
 780                         irv = dsign
 781                                 ? STRTOG_Normal | STRTOG_Inexlo
 782                                 : STRTOG_Normal | STRTOG_Inexhi;
 783                         if (dsign || bbbits > 1 || denorm || rve1 == emin)
 784                                 break;
 785                         delta = lshift(delta,1);
 786                         if (cmp(delta, bs) > 0) {
 787                                 irv = STRTOG_Normal | STRTOG_Inexlo;
 788                                 goto drop_down;
 789                                 }
 790                         break;
 791                         }
 792                 if (i == 0) {
 793                         /* exactly half-way between */
 794                         if (dsign) {
 795                                 if (denorm && all_on(rvb, rvbits)) {
 796                                         /*boundary case -- increment exponent*/
 797                                         rvb->wds = 1;
 798                                         rvb->x[0] = 1;
 799                                         rve = emin + nbits - (rvbits = 1);
 800                                         irv = STRTOG_Normal | STRTOG_Inexhi;
 801                                         denorm = 0;
 802                                         break;
 803                                         }
 804                                 irv = STRTOG_Normal | STRTOG_Inexlo;
 805                                 }
 806                         else if (bbbits == 1) {
 807                                 irv = STRTOG_Normal;
 808  drop_down:
 809                                 /* boundary case -- decrement exponent */
 810                                 if (rve1 == emin) {
 811                                         irv = STRTOG_Normal | STRTOG_Inexhi;
 812                                         if (rvb->wds == 1 && rvb->x[0] == 1)
 813                                                 sudden_underflow = 1;
 814                                         break;
 815                                         }
 816                                 rve -= nbits;
 817                                 rvb = set_ones(rvb, rvbits = nbits);
 818                                 break;
 819                                 }
 820                         else
 821                                 irv = STRTOG_Normal | STRTOG_Inexhi;
 822                         if (bbbits < nbits && !denorm || !(rvb->x[0] & 1))
 823                                 break;
 824                         if (dsign) {
 825                                 rvb = increment(rvb);
 826                                 if ( (j = rvbits >> kshift) !=0)
 827                                         j = 32 - j;
 828                                 if (hi0bits(rvb->x[(rvb->wds - 1) >> kshift])
 829                                                 != j)
 830                                         rvbits++;
 831                                 irv = STRTOG_Normal | STRTOG_Inexhi;
 832                                 }
 833                         else {
 834                                 if (bbbits == 1)
 835                                         goto undfl;
 836                                 decrement(rvb);
 837                                 irv = STRTOG_Normal | STRTOG_Inexlo;
 838                                 }
 839                         break;
 840                         }
 841                 if ((dval(adj) = ratio(delta, bs)) <= 2.) {
 842  adj1:
 843                         inex = STRTOG_Inexlo;
 844                         if (dsign) {
 845                                 asub = 0;
 846                                 inex = STRTOG_Inexhi;
 847                                 }
 848                         else if (denorm && bbbits <= 1) {
 849  undfl:
 850                                 rvb->wds = 0;
 851                                 rve = emin;
 852                                 irv = STRTOG_Underflow | STRTOG_Inexlo;
 853                                 break;
 854                                 }
 855                         adj0 = dval(adj) = 1.;
 856                         }
 857                 else {
 858                         adj0 = dval(adj) *= 0.5;
 859                         if (dsign) {
 860                                 asub = 0;
 861                                 inex = STRTOG_Inexlo;
 862                                 }
 863                         if (dval(adj) < 2147483647.) {
 864                                 L = adj0;
 865                                 adj0 -= L;
 866                                 switch(rd) {
 867                                   case 0:
 868                                         if (adj0 >= .5)
 869                                                 goto inc_L;
 870                                         break;
 871                                   case 1:
 872                                         if (asub && adj0 > 0.)
 873                                                 goto inc_L;
 874                                         break;
 875                                   case 2:
 876                                         if (!asub && adj0 > 0.) {
 877  inc_L:
 878                                                 L++;
 879                                                 inex = STRTOG_Inexact - inex;
 880                                                 }
 881                                   }
 882                                 dval(adj) = L;
 883                                 }
 884                         }
 885                 y = rve + rvbits;
 886
 887                 /* adj *= ulp(dval(rv)); */
 888                 /* if (asub) rv -= adj; else rv += adj; */
 889
 890                 if (!denorm && rvbits < nbits) {
 891                         rvb = lshift(rvb, j = nbits - rvbits);
 892                         rve -= j;
 893                         rvbits = nbits;
 894                         }
 895                 ab = d2b(dval(adj), &abe, &abits);
 896                 if (abe < 0)
 897                         rshift(ab, -abe);
 898                 else if (abe > 0)
 899                         ab = lshift(ab, abe);
 900                 rvb0 = rvb;
 901                 if (asub) {
 902                         /* rv -= adj; */
 903                         j = hi0bits(rvb->x[rvb->wds-1]);
 904                         rvb = diff(rvb, ab);
 905                         k = rvb0->wds - 1;
 906                         if (denorm)
 907                                 /* do nothing */;
 908                         else if (rvb->wds <= k
 909                                 || hi0bits( rvb->x[k]) >
 910                                    hi0bits(rvb0->x[k])) {
 911                                 /* unlikely; can only have lost 1 high bit */
 912                                 if (rve1 == emin) {
 913                                         --rvbits;
 914                                         denorm = 1;
 915                                         }
 916                                 else {
 917                                         rvb = lshift(rvb, 1);
 918                                         --rve;
 919                                         --rve1;
 920                                         L = finished = 0;
 921                                         }
 922                                 }
 923                         }
 924                 else {
 925                         rvb = sum(rvb, ab);
 926                         k = rvb->wds - 1;
 927                         if (k >= rvb0->wds
 928                          || hi0bits(rvb->x[k]) < hi0bits(rvb0->x[k])) {
 929                                 if (denorm) {
 930                                         if (++rvbits == nbits)
 931                                                 denorm = 0;
 932                                         }
 933                                 else {
 934                                         rshift(rvb, 1);
 935                                         rve++;
 936                                         rve1++;
 937                                         L = 0;
 938                                         }
 939                                 }
 940                         }
 941                 Bfree(ab);
 942                 Bfree(rvb0);
 943                 if (finished)
 944                         break;
 945
 946                 z = rve + rvbits;
 947                 if (y == z && L) {
 948                         /* Can we stop now? */
 949                         tol = dval(adj) * 5e-16; /* > max rel error */
 950                         dval(adj) = adj0 - .5;
 951                         if (dval(adj) < -tol) {
 952                                 if (adj0 > tol) {
 953                                         irv |= inex;
 954                                         break;
 955                                         }
 956                                 }
 957                         else if (dval(adj) > tol && adj0 < 1. - tol) {
 958                                 irv |= inex;
 959                                 break;
 960                                 }
 961                         }
 962                 bb0 = denorm ? 0 : trailz(rvb);
 963                 Bfree(bb);
 964                 Bfree(bd);
 965                 Bfree(bs);
 966                 Bfree(delta);
 967                 }
 968         if (!denorm && rvbits < nbits) {
 969                 j = nbits - rvbits;
 970                 rvb = lshift(rvb, j);
 971                 rve -= j;
 972                 }
 973         *exp = rve;
 974         Bfree(bb);
 975         Bfree(bd);
 976         Bfree(bs);
 977         Bfree(bd0);
 978         Bfree(delta);
 979  ret:
 980         if (denorm) {
 981                 if (sudden_underflow) {
 982                         rvb->wds = 0;
 983                         irv = STRTOG_Underflow | STRTOG_Inexlo;
 984                         }
 985                 else  {
 986                         irv = (irv & ~STRTOG_Retmask) |
 987                                 (rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero);
 988                         if (irv & STRTOG_Inexact)
 989                                 irv |= STRTOG_Underflow;
 990                         }
 991                 }
 992         if (se)
 993                 *se = (char *)s;
 994         if (sign)
 995                 irv |= STRTOG_Neg;
 996         if (rvb) {
 997                 copybits(bits, nbits, rvb);
 998                 Bfree(rvb);
 999                 }
1000         return irv;
1001         }