gdtoa/FreeBSD/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) {
 653                 rvb->wds = 0;
 654                 irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
 655 #ifndef NO_ERRNO
 656                 errno = ERANGE;
 657 #endif
 658  infnanexp:
 659                 *exp = fpi->emax + 1;
 660                 goto ret;
 661                 }
 662         rve1 = rve + rvbits - nbits;
 663         if (e2 < (emin = fpi->emin)) {
 664                 denorm = 1;
 665                 j = rve - emin;
 666                 if (j > 0) {
 667                         rvb = lshift(rvb, j);
 668                         rvbits += j;
 669                         }
 670                 else if (j < 0) {
 671                         rvbits += j;
 672                         if (rvbits <= 0) {
 673                                 if (rvbits < -1) {
 674  ufl:
 675                                         rvb->wds = 0;
 676                                         rvb->x[0] = 0;
 677                                         *exp = emin;
 678                                         irv = STRTOG_Underflow | STRTOG_Inexlo;
 679                                         goto ret;
 680                                         }
 681                                 rvb->x[0] = rvb->wds = rvbits = 1;
 682                                 }
 683                         else
 684                                 rshift(rvb, -j);
 685                         }
 686                 rve = rve1 = emin;
 687                 if (sudden_underflow && e2 + 1 < emin)
 688                         goto ufl;
 689                 }
 690
 691         /* Now the hard part -- adjusting rv to the correct value.*/
 692
 693         /* Put digits into bd: true value = bd * 10^e */
 694
 695         bd0 = s2b(s0, nd0, nd, y);
 696
 697         for(;;) {
 698                 bd = Balloc(bd0->k);
 699                 Bcopy(bd, bd0);
 700                 bb = Balloc(rvb->k);
 701                 Bcopy(bb, rvb);
 702                 bbbits = rvbits - bb0;
 703                 bbe = rve + bb0;
 704                 bs = i2b(1);
 705
 706                 if (e >= 0) {
 707                         bb2 = bb5 = 0;
 708                         bd2 = bd5 = e;
 709                         }
 710                 else {
 711                         bb2 = bb5 = -e;
 712                         bd2 = bd5 = 0;
 713                         }
 714                 if (bbe >= 0)
 715                         bb2 += bbe;
 716                 else
 717                         bd2 -= bbe;
 718                 bs2 = bb2;
 719                 j = nbits + 1 - bbbits;
 720                 i = bbe + bbbits - nbits;
 721                 if (i < emin)   /* denormal */
 722                         j += i - emin;
 723                 bb2 += j;
 724                 bd2 += j;
 725                 i = bb2 < bd2 ? bb2 : bd2;
 726                 if (i > bs2)
 727                         i = bs2;
 728                 if (i > 0) {
 729                         bb2 -= i;
 730                         bd2 -= i;
 731                         bs2 -= i;
 732                         }
 733                 if (bb5 > 0) {
 734                         bs = pow5mult(bs, bb5);
 735                         bb1 = mult(bs, bb);
 736                         Bfree(bb);
 737                         bb = bb1;
 738                         }
 739                 bb2 -= bb0;
 740                 if (bb2 > 0)
 741                         bb = lshift(bb, bb2);
 742                 else if (bb2 < 0)
 743                         rshift(bb, -bb2);
 744                 if (bd5 > 0)
 745                         bd = pow5mult(bd, bd5);
 746                 if (bd2 > 0)
 747                         bd = lshift(bd, bd2);
 748                 if (bs2 > 0)
 749                         bs = lshift(bs, bs2);
 750                 asub = 1;
 751                 inex = STRTOG_Inexhi;
 752                 delta = diff(bb, bd);
 753                 if (delta->wds <= 1 && !delta->x[0])
 754                         break;
 755                 dsign = delta->sign;
 756                 delta->sign = finished = 0;
 757                 L = 0;
 758                 i = cmp(delta, bs);
 759                 if (rd && i <= 0) {
 760                         irv = STRTOG_Normal;
 761                         if ( (finished = dsign ^ (rd&1)) !=0) {
 762                                 if (dsign != 0) {
 763                                         irv |= STRTOG_Inexhi;
 764                                         goto adj1;
 765                                         }
 766                                 irv |= STRTOG_Inexlo;
 767                                 if (rve1 == emin)
 768                                         goto adj1;
 769                                 for(i = 0, j = nbits; j >= ULbits;
 770                                                 i++, j -= ULbits) {
 771                                         if (rvb->x[i] & ALL_ON)
 772                                                 goto adj1;
 773                                         }
 774                                 if (j > 1 && lo0bits(rvb->x + i) < j - 1)
 775                                         goto adj1;
 776                                 rve = rve1 - 1;
 777                                 rvb = set_ones(rvb, rvbits = nbits);
 778                                 break;
 779                                 }
 780                         irv |= dsign ? STRTOG_Inexlo : STRTOG_Inexhi;
 781                         break;
 782                         }
 783                 if (i < 0) {
 784                         /* Error is less than half an ulp -- check for
 785                          * special case of mantissa a power of two.
 786                          */
 787                         irv = dsign
 788                                 ? STRTOG_Normal | STRTOG_Inexlo
 789                                 : STRTOG_Normal | STRTOG_Inexhi;
 790                         if (dsign || bbbits > 1 || denorm || rve1 == emin)
 791                                 break;
 792                         delta = lshift(delta,1);
 793                         if (cmp(delta, bs) > 0) {
 794                                 irv = STRTOG_Normal | STRTOG_Inexlo;
 795                                 goto drop_down;
 796                                 }
 797                         break;
 798                         }
 799                 if (i == 0) {
 800                         /* exactly half-way between */
 801                         if (dsign) {
 802                                 if (denorm && all_on(rvb, rvbits)) {
 803                                         /*boundary case -- increment exponent*/
 804                                         rvb->wds = 1;
 805                                         rvb->x[0] = 1;
 806                                         rve = emin + nbits - (rvbits = 1);
 807                                         irv = STRTOG_Normal | STRTOG_Inexhi;
 808                                         denorm = 0;
 809                                         break;
 810                                         }
 811                                 irv = STRTOG_Normal | STRTOG_Inexlo;
 812                                 }
 813                         else if (bbbits == 1) {
 814                                 irv = STRTOG_Normal;
 815  drop_down:
 816                                 /* boundary case -- decrement exponent */
 817                                 if (rve1 == emin) {
 818                                         irv = STRTOG_Normal | STRTOG_Inexhi;
 819                                         if (rvb->wds == 1 && rvb->x[0] == 1)
 820                                                 sudden_underflow = 1;
 821                                         break;
 822                                         }
 823                                 rve -= nbits;
 824                                 rvb = set_ones(rvb, rvbits = nbits);
 825                                 break;
 826                                 }
 827                         else
 828                                 irv = STRTOG_Normal | STRTOG_Inexhi;
 829                         if (bbbits < nbits && !denorm || !(rvb->x[0] & 1))
 830                                 break;
 831                         if (dsign) {
 832                                 rvb = increment(rvb);
 833                                 if ( (j = rvbits & kmask) !=0)
 834                                         j = ULbits - j;
 835                                 if (hi0bits(rvb->x[(rvb->wds - 1) >> kshift])
 836                                                 != j)
 837                                         rvbits++;
 838                                 irv = STRTOG_Normal | STRTOG_Inexhi;
 839                                 }
 840                         else {
 841                                 if (bbbits == 1)
 842                                         goto undfl;
 843                                 decrement(rvb);
 844                                 irv = STRTOG_Normal | STRTOG_Inexlo;
 845                                 }
 846                         break;
 847                         }
 848                 if ((dval(adj) = ratio(delta, bs)) <= 2.) {
 849  adj1:
 850                         inex = STRTOG_Inexlo;
 851                         if (dsign) {
 852                                 asub = 0;
 853                                 inex = STRTOG_Inexhi;
 854                                 }
 855                         else if (denorm && bbbits <= 1) {
 856  undfl:
 857                                 rvb->wds = 0;
 858                                 rve = emin;
 859                                 irv = STRTOG_Underflow | STRTOG_Inexlo;
 860                                 break;
 861                                 }
 862                         adj0 = dval(adj) = 1.;
 863                         }
 864                 else {
 865                         adj0 = dval(adj) *= 0.5;
 866                         if (dsign) {
 867                                 asub = 0;
 868                                 inex = STRTOG_Inexlo;
 869                                 }
 870                         if (dval(adj) < 2147483647.) {
 871                                 L = adj0;
 872                                 adj0 -= L;
 873                                 switch(rd) {
 874                                   case 0:
 875                                         if (adj0 >= .5)
 876                                                 goto inc_L;
 877                                         break;
 878                                   case 1:
 879                                         if (asub && adj0 > 0.)
 880                                                 goto inc_L;
 881                                         break;
 882                                   case 2:
 883                                         if (!asub && adj0 > 0.) {
 884  inc_L:
 885                                                 L++;
 886                                                 inex = STRTOG_Inexact - inex;
 887                                                 }
 888                                   }
 889                                 dval(adj) = L;
 890                                 }
 891                         }
 892                 y = rve + rvbits;
 893
 894                 /* adj *= ulp(dval(rv)); */
 895                 /* if (asub) rv -= adj; else rv += adj; */
 896
 897                 if (!denorm && rvbits < nbits) {
 898                         rvb = lshift(rvb, j = nbits - rvbits);
 899                         rve -= j;
 900                         rvbits = nbits;
 901                         }
 902                 ab = d2b(dval(adj), &abe, &abits);
 903                 if (abe < 0)
 904                         rshift(ab, -abe);
 905                 else if (abe > 0)
 906                         ab = lshift(ab, abe);
 907                 rvb0 = rvb;
 908                 if (asub) {
 909                         /* rv -= adj; */
 910                         j = hi0bits(rvb->x[rvb->wds-1]);
 911                         rvb = diff(rvb, ab);
 912                         k = rvb0->wds - 1;
 913                         if (denorm)
 914                                 /* do nothing */;
 915                         else if (rvb->wds <= k
 916                                 || hi0bits( rvb->x[k]) >
 917                                    hi0bits(rvb0->x[k])) {
 918                                 /* unlikely; can only have lost 1 high bit */
 919                                 if (rve1 == emin) {
 920                                         --rvbits;
 921                                         denorm = 1;
 922                                         }
 923                                 else {
 924                                         rvb = lshift(rvb, 1);
 925                                         --rve;
 926                                         --rve1;
 927                                         L = finished = 0;
 928                                         }
 929                                 }
 930                         }
 931                 else {
 932                         rvb = sum(rvb, ab);
 933                         k = rvb->wds - 1;
 934                         if (k >= rvb0->wds
 935                          || hi0bits(rvb->x[k]) < hi0bits(rvb0->x[k])) {
 936                                 if (denorm) {
 937                                         if (++rvbits == nbits)
 938                                                 denorm = 0;
 939                                         }
 940                                 else {
 941                                         rshift(rvb, 1);
 942                                         rve++;
 943                                         rve1++;
 944                                         L = 0;
 945                                         }
 946                                 }
 947                         }
 948                 Bfree(ab);
 949                 Bfree(rvb0);
 950                 if (finished)
 951                         break;
 952
 953                 z = rve + rvbits;
 954                 if (y == z && L) {
 955                         /* Can we stop now? */
 956                         tol = dval(adj) * 5e-16; /* > max rel error */
 957                         dval(adj) = adj0 - .5;
 958                         if (dval(adj) < -tol) {
 959                                 if (adj0 > tol) {
 960                                         irv |= inex;
 961                                         break;
 962                                         }
 963                                 }
 964                         else if (dval(adj) > tol && adj0 < 1. - tol) {
 965                                 irv |= inex;
 966                                 break;
 967                                 }
 968                         }
 969                 bb0 = denorm ? 0 : trailz(rvb);
 970                 Bfree(bb);
 971                 Bfree(bd);
 972                 Bfree(bs);
 973                 Bfree(delta);
 974                 }
 975         if (!denorm && (j = nbits - rvbits)) {
 976                 if (j > 0)
 977                         rvb = lshift(rvb, j);
 978                 else
 979                         rshift(rvb, -j);
 980                 rve -= j;
 981                 }
 982         *exp = rve;
 983         Bfree(bb);
 984         Bfree(bd);
 985         Bfree(bs);
 986         Bfree(bd0);
 987         Bfree(delta);
 988  ret:
 989         if (denorm) {
 990                 if (sudden_underflow) {
 991                         rvb->wds = 0;
 992                         irv = STRTOG_Underflow | STRTOG_Inexlo;
 993                         }
 994                 else  {
 995                         irv = (irv & ~STRTOG_Retmask) |
 996                                 (rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero);
 997                         if (irv & STRTOG_Inexact)
 998                                 irv |= STRTOG_Underflow;
 999                         }
1000                 }
1001         if (se)
1002                 *se = (char *)s;
1003         if (sign)
1004                 irv |= STRTOG_Neg;
1005         if (rvb) {
1006                 copybits(bits, nbits, rvb);
1007                 Bfree(rvb);
1008                 }
1009         return irv;
1010         }