StdLib/LibC/gdtoa/strtodg.c

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