StdLib/LibC/gdtoa/strtodg.c

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