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