contrib/gdtoa/strtodg.c

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