git.saurik.com Git - apple/libc.git/blob - gdtoa/FreeBSD/gdtoa-strtodg.c

1 /****************************************************************

3 The author of this software is David M. Gay.

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.

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.

27 ****************************************************************/

29 /* Please send bug reports to David M. Gay (dmg at acm dot org,

30 * with " at " changed at "@" and " dot " changed to "."). */

32 #include "xlocale_private.h"

34 #include "gdtoaimp.h"

36 #ifdef USE_LOCALE

37 #include "locale.h"

38 #endif

40 #define fivesbits __fivesbits_D2A

41 #define all_on __all_on_D2A

42 #define set_ones __set_ones_D2A

43 #define rvOK __rvOK_D2A

44 #define mantbits __mantbits_D2A

46 #ifdef BUILDING_VARIANT

47 extern CONST int fivesbits[];

 int all_on(Bigint *b, int n);

 Bigint *set_ones(Bigint *b, int n);

 int rvOK(U *d, CONST FPI *fpi, Long *exp, ULong *bits, int exact, int rd, int *irv);

 int mantbits(U *d);

52 #else /* !BUILDING_VARIANT */

54 __private_extern__ CONST int

 fivesbits[] = {  0,  3,  5,  7, 10, 12, 14, 17, 19, 21,

                 24, 26, 28, 31, 33, 35, 38, 40, 42, 45,

                 47, 49, 52

58 #ifdef VAX

59 , 54, 56

60 #endif

61 };

63 Bigint *

64 #ifdef KR_headers

 increment(b) Bigint *b;

66 #else

67 increment(Bigint *b)

68 #endif

69 {

70 ULong *x, *xe;

71 Bigint *b1;

72 #ifdef Pack_16

73 ULong carry = 1, y;

74 #endif

76 x = b->x;

77 xe = x + b->wds;

78 #ifdef Pack_32

79 do {

                 if (*x < (ULong)0xffffffffL) {

81 ++*x;

82 return b;

83 }

84 *x++ = 0;

                 } while(x < xe);

86 #else

87 do {

88 y = *x + carry;

89 carry = y >> 16;

                 *x++ = y & 0xffff;

91 if (!carry)

92 return b;

                 } while(x < xe);

94 if (carry)

95 #endif

96 {

                 if (b->wds >= b->maxwds) {

                         b1 = Balloc(b->k+1);

99 Bcopy(b1,b);

100 Bfree(b);

101 b = b1;

102 }

                 b->x[b->wds++] = 1;

104 }

105 return b;

106 }

107

108 void

109 #ifdef KR_headers

 decrement(b) Bigint *b;

111 #else

112 decrement(Bigint *b)

113 #endif

114 {

115 ULong *x, *xe;

116 #ifdef Pack_16

117 ULong borrow = 1, y;

118 #endif

119

120 x = b->x;

121 xe = x + b->wds;

122 #ifdef Pack_32

123 do {

124 if (*x) {

125 --*x;

126 break;

127 }

128 *x++ = 0xffffffffL;

129 }

130 while(x < xe);

131 #else

132 do {

133 y = *x - borrow;

                 borrow = (y & 0x10000) >> 16;

                 *x++ = y & 0xffff;

                 } while(borrow && x < xe);

137 #endif

138 }

139

140 __private_extern__ int

141 #ifdef KR_headers

 all_on(b, n) Bigint *b; int n;

143 #else

 all_on(Bigint *b, int n)

145 #endif

146 {

147 ULong *x, *xe;

148

149 x = b->x;

150 xe = x + (n >> kshift);

151 while(x < xe)

                 if ((*x++ & ALL_ON) != ALL_ON)

153 return 0;

154 if (n &= kmask)

                 return ((*x | (ALL_ON << n)) & ALL_ON) == ALL_ON;

156 return 1;

157 }

158

159 Bigint *

160 #ifdef KR_headers

 set_ones(b, n) Bigint *b; int n;

162 #else

 set_ones(Bigint *b, int n)

164 #endif

165 {

166 int k;

167 ULong *x, *xe;

168

         k = (n + ((1 << kshift) - 1)) >> kshift;

         if (b->k < k) {

171 Bfree(b);

172 b = Balloc(k);

173 }

174 k = n >> kshift;

175 if (n &= kmask)

176 k++;

177 b->wds = k;

178 x = b->x;

179 xe = x + k;

180 while(x < xe)

181 *x++ = ALL_ON;

182 if (n)

                 x[-1] >>= ULbits - n;

184 return b;

185 }

186

187 __private_extern__ int

188 rvOK

189 #ifdef KR_headers

  (d, fpi, exp, bits, exact, rd, irv)

  U *d; CONST FPI *fpi; Long *exp; ULong *bits; int exact, rd, *irv;

192 #else

  (U *d, CONST FPI *fpi, Long *exp, ULong *bits, int exact, int rd, int *irv)

194 #endif

195 {

196 Bigint *b;

197 ULong carry, inex, lostbits;

         int bdif, e, j, k, k1, nb, rv;

199

200 carry = rv = 0;

         b = d2b(dval(d), &e, &bdif);

202 bdif -= nb = fpi->nbits;

203 e += bdif;

         if (bdif <= 0) {

205 if (exact)

206 goto trunc;

207 goto ret;

208 }

209 if (P == nb) {

210 if (

211 #ifndef IMPRECISE_INEXACT

212 exact &&

213 #endif

214 fpi->rounding ==

215 #ifdef RND_PRODQUOT

216 FPI_Round_near

217 #else

218 Flt_Rounds

219 #endif

220 ) goto trunc;

221 goto ret;

222 }

223 switch(rd) {

224 case 1: /* round down (toward -Infinity) */

225 goto trunc;

226 case 2: /* round up (toward +Infinity) */

227 break;

228 default: /* round near */

229 k = bdif - 1;

                 if (k < 0)

231 goto trunc;

232 if (!k) {

233 if (!exact)

234 goto ret;

                         if (b->x[0] & 2)

236 break;

237 goto trunc;

238 }

                 if (b->x[k>>kshift] & ((ULong)1 << (k & kmask)))

240 break;

241 goto trunc;

242 }

243 /* "break" cases: round up 1 bit, then truncate; bdif > 0 */

244 carry = 1;

245 trunc:

246 inex = lostbits = 0;

         if (bdif > 0) {

                 if ( (lostbits = any_on(b, bdif)) !=0)

249 inex = STRTOG_Inexlo;

250 rshift(b, bdif);

251 if (carry) {

252 inex = STRTOG_Inexhi;

253 b = increment(b);

                         if ( (j = nb & kmask) !=0)

255 j = ULbits - j;

                         if (hi0bits(b->x[b->wds - 1]) != j) {

257 if (!lostbits)

                                         lostbits = b->x[0] & 1;

                                 rshift(b, 1);

260 e++;

261 }

262 }

263 }

         else if (bdif < 0)

                 b = lshift(b, -bdif);

         if (e < fpi->emin) {

267 k = fpi->emin - e;

268 e = fpi->emin;

                 if (k > nb || fpi->sudden_underflow) {

                         b->wds = inex = 0;

271 *irv = STRTOG_Underflow | STRTOG_Inexlo;

272 }

273 else {

274 k1 = k - 1;

                         if (k1 > 0 && !lostbits)

                                 lostbits = any_on(b, k1);

277 if (!lostbits && !exact)

278 goto ret;

279 lostbits |=

                           carry = b->x[k1>>kshift] & (1 << (k1 & kmask));

281 rshift(b, k);

282 *irv = STRTOG_Denormal;

283 if (carry) {

284 b = increment(b);

285 inex = STRTOG_Inexhi | STRTOG_Underflow;

286 }

287 else if (lostbits)

288 inex = STRTOG_Inexlo | STRTOG_Underflow;

289 }

290 }

         else if (e > fpi->emax) {

                 e = fpi->emax + 1;

                 *irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;

294 #ifndef NO_ERRNO

295 errno = ERANGE;

296 #endif

                 b->wds = inex = 0;

298 }

299 *exp = e;

         copybits(bits, nb, b);

301 *irv |= inex;

302 rv = 1;

303 ret:

304 Bfree(b);

305 return rv;

306 }

307

308 __private_extern__ int

309 #ifdef KR_headers

 mantbits(d) U *d;

311 #else

312 mantbits(U *d)

313 #endif

314 {

315 ULong L;

316 #ifdef VAX

         L = word1(d) << 16 | word1(d) >> 16;

318 if (L)

319 #else

         if ( (L = word1(d)) !=0)

321 #endif

                 return P - lo0bits(&L);

323 #ifdef VAX

         L = word0(d) << 16 | word0(d) >> 16 | Exp_msk11;

325 #else

         L = word0(d) | Exp_msk1;

327 #endif

         return P - 32 - lo0bits(&L);

329 }

330

331 #endif /* BUILDING_VARIANT */

332

333 int

334 strtodg

335 #ifdef KR_headers

         (s00, se, fpi, exp, bits, loc)

         CONST char *s00; char **se; CONST FPI *fpi; Long *exp; ULong *bits; locale_t loc;

338 #else

         (CONST char *s00, char **se, CONST FPI *fpi, Long *exp, ULong *bits, locale_t loc)

340 #endif

341 {

342 int abe, abits, asub;

         int bb0, bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, denorm;

         int dsign, e, e1, e2, emin, esign, finished, i, inex, irv;

         int j, k, nbits, nd, nd0, nf, nz, nz0, rd, rvbits, rve, rve1, sign;

346 int sudden_underflow;

         CONST char *s, *s0, *s1;

348 char *strunc = NULL;

349 double adj0, tol;

350 Long L;

351 U adj, rv;

         ULong *b, *be, y, z;

         Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0;

354 #ifdef USE_LOCALE /*{{*/

355 NORMALIZE_LOCALE(loc);

356 #ifdef NO_LOCALE_CACHE

         char *decimalpoint = localeconv_l(loc)->decimal_point;

         int dplen = strlen(decimalpoint);

359 #else

360 char *decimalpoint;

361 static char *decimalpoint_cache;

362 static int dplen;

363 if (!(s0 = decimalpoint_cache)) {

                 s0 = localeconv_l(loc)->decimal_point;

                 if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {

366 strcpy(decimalpoint_cache, s0);

367 s0 = decimalpoint_cache;

368 }

369 dplen = strlen(s0);

370 }

371 decimalpoint = (char*)s0;

372 #endif /*NO_LOCALE_CACHE*/

373 #else /*USE_LOCALE}{*/

374 #define dplen 1

375 #endif /*USE_LOCALE}}*/

376

377 irv = STRTOG_Zero;

         denorm = sign = nz0 = nz = 0;

         dval(&rv) = 0.;

380 rvb = 0;

381 nbits = fpi->nbits;

         for(s = s00;;s++) switch(*s) {

383 case '-':

384 sign = 1;

385 /* no break */

386 case '+':

387 if (*++s)

388 goto break2;

389 /* no break */

390 case 0:

391 sign = 0;

392 irv = STRTOG_NoNumber;

393 s = s00;

394 goto ret;

                 case '\t':

                 case '\n':

                 case '\v':

                 case '\f':

                 case '\r':

400 case ' ':

401 continue;

402 default:

403 goto break2;

404 }

405 break2:

         if (*s == '0') {

407 #ifndef NO_HEX_FP

                 switch(s[1]) {

409 case 'x':

410 case 'X':

                         irv = gethex(&s, fpi, exp, &rvb, sign, loc);

412 if (irv == STRTOG_NoNumber) {

413 s = s00;

414 sign = 0;

415 }

416 goto ret;

417 }

418 #endif

419 nz0 = 1;

                 while(*++s == '0') ;

421 if (!*s)

422 goto ret;

423 }

424 sudden_underflow = fpi->sudden_underflow;

425 s0 = s;

426 y = z = 0;

         for(decpt = nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)

                 if (nd < 9)

                         y = 10*y + c - '0';

                 else if (nd < 16)

                         z = 10*z + c - '0';

432 nd0 = nd;

433 #ifdef USE_LOCALE

434 if (c == *decimalpoint) {

                 for(i = 1; decimalpoint[i]; ++i)

                         if (s[i] != decimalpoint[i])

437 goto dig_done;

438 s += i;

439 c = *s;

440 #else

         if (c == '.') {

442 c = *++s;

443 #endif

444 decpt = 1;

445 if (!nd) {

                         for(; c == '0'; c = *++s)

447 nz++;

                         if (c > '0' && c <= '9') {

449 s0 = s;

450 nf += nz;

451 nz = 0;

452 goto have_dig;

453 }

454 goto dig_done;

455 }

                 for(; c >= '0' && c <= '9'; c = *++s) {

457 have_dig:

458 nz++;

                         if (c -= '0') {

460 nf += nz;

                                 for(i = 1; i < nz; i++)

                                         if (nd++ < 9)

463 y *= 10;

                                         else if (nd <= DBL_DIG + 1)

465 z *= 10;

                                 if (nd++ < 9)

                                         y = 10*y + c;

                                 else if (nd <= DBL_DIG + 1)

                                         z = 10*z + c;

470 nz = 0;

471 }

472 }

473 }/*}*/

474 dig_done:

475 e = 0;

         if (c == 'e' || c == 'E') {

                 if (!nd && !nz && !nz0) {

478 irv = STRTOG_NoNumber;

479 s = s00;

480 goto ret;

481 }

482 s00 = s;

483 esign = 0;

484 switch(c = *++s) {

485 case '-':

486 esign = 1;

487 case '+':

488 c = *++s;

489 }

                 if (c >= '0' && c <= '9') {

                         while(c == '0')

492 c = *++s;

                         if (c > '0' && c <= '9') {

494 L = c - '0';

495 s1 = s;

                                 while((c = *++s) >= '0' && c <= '9')

                                         L = 10*L + c - '0';

                                 if (s - s1 > 8 || L > 19999)

499 /* Avoid confusion from exponents

500 * so large that e might overflow.

501 */

502 e = 19999; /* safe for 16 bit ints */

503 else

504 e = (int)L;

505 if (esign)

506 e = -e;

507 }

508 else

509 e = 0;

510 }

511 else

512 s = s00;

513 }

514 if (!nd) {

515 if (!nz && !nz0) {

516 #ifdef INFNAN_CHECK

517 /* Check for Nan and Infinity */

518 if (!decpt)

519 switch(c) {

520 case 'i':

521 case 'I':

                                 if (match(&s,"nf")) {

523 --s;

                                         if (!match(&s,"inity"))

525 ++s;

526 irv = STRTOG_Infinite;

527 goto infnanexp;

528 }

529 break;

530 case 'n':

531 case 'N':

                                 if (match(&s, "an")) {

533 irv = STRTOG_NaN;

                                         *exp = fpi->emax + 1;

535 #ifndef No_Hex_NaN

                                         if (*s == '(') /*)*/

                                                 irv = hexnan(&s, fpi, bits);

538 #endif

539 goto infnanexp;

540 }

541 }

542 #endif /* INFNAN_CHECK */

543 irv = STRTOG_NoNumber;

544 s = s00;

545 }

546 goto ret;

547 }

         TRUNCATE_DIGITS(s0, strunc, nd, nd0, nf, fpi->nbits, fpi->emin, dplen);

549

550 irv = STRTOG_Normal;

551 e1 = e -= nf;

552 rd = 0;

         switch(fpi->rounding & 3) {

554 case FPI_Round_up:

555 rd = 2 - sign;

556 break;

557 case FPI_Round_zero:

558 rd = 1;

559 break;

560 case FPI_Round_down:

561 rd = 1 + sign;

562 }

563

564 /* Now we have nd0 digits, starting at s0, followed by a

565 * decimal point, followed by nd-nd0 digits. The number we're

566 * after is the integer represented by those digits times

567 * 10**e */

568

569 if (!nd0)

570 nd0 = nd;

         k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;

572 dval(&rv) = y;

         if (k > 9)

                 dval(&rv) = tens[k - 9] * dval(&rv) + z;

575 bd0 = 0;

         if (nbits <= P && nd <= DBL_DIG) {

577 if (!e) {

                         if (rvOK(&rv, fpi, exp, bits, 1, rd, &irv))

579 goto ret;

580 }

                 else if (e > 0) {

582 if (e <= Ten_pmax) {

583 #ifdef VAX

584 goto vax_ovfl_check;

585 #else

                                 i = fivesbits[e] + mantbits(&rv) <= P;

                                 /* rv = */ rounded_product(dval(&rv), tens[e]);

                                 if (rvOK(&rv, fpi, exp, bits, i, rd, &irv))

589 goto ret;

590 e1 -= e;

591 goto rv_notOK;

592 #endif

593 }

594 i = DBL_DIG - nd;

                         if (e <= Ten_pmax + i) {

596 /* A fancier test would sometimes let us do

597 * this for larger i values.

598 */

599 e2 = e - i;

600 e1 -= i;

                                 dval(&rv) *= tens[i];

602 #ifdef VAX

603 /* VAX exponent range is so narrow we must

604 * worry about overflow here...

605 */

606 vax_ovfl_check:

                                 dval(&adj) = dval(&rv);

                                 word0(&adj) -= P*Exp_msk1;

                                 /* adj = */ rounded_product(dval(&adj), tens[e2]);

                                 if ((word0(&adj) & Exp_mask)

                                  > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))

612 goto rv_notOK;

                                 word0(&adj) += P*Exp_msk1;

                                 dval(&rv) = dval(&adj);

615 #else

                                 /* rv = */ rounded_product(dval(&rv), tens[e2]);

617 #endif

                                 if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv))

619 goto ret;

620 e1 -= e2;

621 }

622 }

623 #ifndef Inaccurate_Divide

624 else if (e >= -Ten_pmax) {

                         /* rv = */ rounded_quotient(dval(&rv), tens[-e]);

                         if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv))

627 goto ret;

628 e1 -= e;

629 }

630 #endif

631 }

632 rv_notOK:

633 e1 += nd - k;

634

635 /* Get starting approximation = rv * 10**e1 */

636

637 e2 = 0;

         if (e1 > 0) {

                 if ( (i = e1 & 15) !=0)

                         dval(&rv) *= tens[i];

                 if (e1 &= ~15) {

642 e1 >>= 4;

                         while(e1 >= (1 << (n_bigtens-1))) {

                                 e2 += ((word0(&rv) & Exp_mask)

645 >> Exp_shift1) - Bias;

646 word0(&rv) &= ~Exp_mask;

                                 word0(&rv) |= Bias << Exp_shift1;

                                 dval(&rv) *= bigtens[n_bigtens-1];

                                 e1 -= 1 << (n_bigtens-1);

650 }

                         e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;

652 word0(&rv) &= ~Exp_mask;

                         word0(&rv) |= Bias << Exp_shift1;

                         for(j = 0; e1 > 0; j++, e1 >>= 1)

                                 if (e1 & 1)

                                         dval(&rv) *= bigtens[j];

657 }

658 }

         else if (e1 < 0) {

660 e1 = -e1;

                 if ( (i = e1 & 15) !=0)

                         dval(&rv) /= tens[i];

                 if (e1 &= ~15) {

664 e1 >>= 4;

                         while(e1 >= (1 << (n_bigtens-1))) {

                                 e2 += ((word0(&rv) & Exp_mask)

667 >> Exp_shift1) - Bias;

668 word0(&rv) &= ~Exp_mask;

                                 word0(&rv) |= Bias << Exp_shift1;

                                 dval(&rv) *= tinytens[n_bigtens-1];

                                 e1 -= 1 << (n_bigtens-1);

672 }

                         e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;

674 word0(&rv) &= ~Exp_mask;

                         word0(&rv) |= Bias << Exp_shift1;

                         for(j = 0; e1 > 0; j++, e1 >>= 1)

                                 if (e1 & 1)

                                         dval(&rv) *= tinytens[j];

679 }

680 }

681 #ifdef IBM

682 /* e2 is a correction to the (base 2) exponent of the return

683 * value, reflecting adjustments above to avoid overflow in the

684 * native arithmetic. For native IBM (base 16) arithmetic, we

685 * must multiply e2 by 4 to change from base 16 to 2.

686 */

687 e2 <<= 2;

688 #endif

         rvb = d2b(dval(&rv), &rve, &rvbits);    /* rv = rvb * 2^rve */

690 rve += e2;

         if ((j = rvbits - nbits) > 0) {

692 rshift(rvb, j);

693 rvbits = nbits;

694 rve += j;

695 }

696 bb0 = 0; /* trailing zero bits in rvb */

697 e2 = rve + rvbits - nbits;

         if (e2 > fpi->emax + 1)

699 goto huge;

700 rve1 = rve + rvbits - nbits;

         if (e2 < (emin = fpi->emin)) {

702 denorm = 1;

703 j = rve - emin;

                 if (j > 0) {

                         rvb = lshift(rvb, j);

706 rvbits += j;

707 }

                 else if (j < 0) {

709 rvbits += j;

                         if (rvbits <= 0) {

                                 if (rvbits < -1) {

712 ufl:

713 rvb->wds = 0;

                                         rvb->x[0] = 0;

715 *exp = emin;

716 irv = STRTOG_Underflow | STRTOG_Inexlo;

717 /* When __DARWIN_UNIX03 is set, we don't need this (errno is set later) */

718 #if !defined(NO_ERRNO) && !__DARWIN_UNIX03

719 errno = ERANGE;

720 #endif

721 goto ret;

722 }

                                 rvb->x[0] = rvb->wds = rvbits = 1;

724 }

725 else

726 rshift(rvb, -j);

727 }

728 rve = rve1 = emin;

                 if (sudden_underflow && e2 + 1 < emin)

730 goto ufl;

731 }

732

733 /* Now the hard part -- adjusting rv to the correct value.*/

734

735 /* Put digits into bd: true value = bd * 10^e */

736

         bd0 = s2b(s0, nd0, nd, y, dplen);

738

739 for(;;) {

                 bd = Balloc(bd0->k);

741 Bcopy(bd, bd0);

                 bb = Balloc(rvb->k);

743 Bcopy(bb, rvb);

744 bbbits = rvbits - bb0;

745 bbe = rve + bb0;

                 bs = i2b(1);

747

                 if (e >= 0) {

749 bb2 = bb5 = 0;

750 bd2 = bd5 = e;

751 }

752 else {

753 bb2 = bb5 = -e;

754 bd2 = bd5 = 0;

755 }

                 if (bbe >= 0)

757 bb2 += bbe;

758 else

759 bd2 -= bbe;

760 bs2 = bb2;

                 j = nbits + 1 - bbbits;

762 i = bbe + bbbits - nbits;

                 if (i < emin)   /* denormal */

764 j += i - emin;

765 bb2 += j;

766 bd2 += j;

                 i = bb2 < bd2 ? bb2 : bd2;

768 if (i > bs2)

769 i = bs2;

                 if (i > 0) {

771 bb2 -= i;

772 bd2 -= i;

773 bs2 -= i;

774 }

                 if (bb5 > 0) {

                         bs = pow5mult(bs, bb5);

                         bb1 = mult(bs, bb);

778 Bfree(bb);

779 bb = bb1;

780 }

781 bb2 -= bb0;

                 if (bb2 > 0)

                         bb = lshift(bb, bb2);

                 else if (bb2 < 0)

785 rshift(bb, -bb2);

                 if (bd5 > 0)

                         bd = pow5mult(bd, bd5);

                 if (bd2 > 0)

                         bd = lshift(bd, bd2);

                 if (bs2 > 0)

                         bs = lshift(bs, bs2);

792 asub = 1;

793 inex = STRTOG_Inexhi;

                 delta = diff(bb, bd);

                 if (delta->wds <= 1 && !delta->x[0])

796 break;

797 dsign = delta->sign;

                 delta->sign = finished = 0;

799 L = 0;

                 i = cmp(delta, bs);

                 if (rd && i <= 0) {

802 irv = STRTOG_Normal;

                         if ( (finished = dsign ^ (rd&1)) !=0) {

                                 if (dsign != 0) {

805 irv |= STRTOG_Inexhi;

806 goto adj1;

807 }

808 irv |= STRTOG_Inexlo;

809 if (rve1 == emin)

810 goto adj1;

                                 for(i = 0, j = nbits; j >= ULbits;

812 i++, j -= ULbits) {

                                         if (rvb->x[i] & ALL_ON)

814 goto adj1;

815 }

                                 if (j > 1 && lo0bits(rvb->x + i) < j - 1)

817 goto adj1;

818 rve = rve1 - 1;

                                 rvb = set_ones(rvb, rvbits = nbits);

820 break;

821 }

822 irv |= dsign ? STRTOG_Inexlo : STRTOG_Inexhi;

823 break;

824 }

                 if (i < 0) {

826 /* Error is less than half an ulp -- check for

827 * special case of mantissa a power of two.

828 */

829 irv = dsign

830 ? STRTOG_Normal | STRTOG_Inexlo

831 : STRTOG_Normal | STRTOG_Inexhi;

                         if (dsign || bbbits > 1 || denorm || rve1 == emin)

833 break;

                         delta = lshift(delta,1);

                         if (cmp(delta, bs) > 0) {

836 irv = STRTOG_Normal | STRTOG_Inexlo;

837 goto drop_down;

838 }

839 break;

840 }

                 if (i == 0) {

842 /* exactly half-way between */

843 if (dsign) {

                                 if (denorm && all_on(rvb, rvbits)) {

845 /*boundary case -- increment exponent*/

846 rvb->wds = 1;

                                         rvb->x[0] = 1;

                                         rve = emin + nbits - (rvbits = 1);

849 irv = STRTOG_Normal | STRTOG_Inexhi;

850 denorm = 0;

851 break;

852 }

853 irv = STRTOG_Normal | STRTOG_Inexlo;

854 }

                         else if (bbbits == 1) {

856 irv = STRTOG_Normal;

857 drop_down:

858 /* boundary case -- decrement exponent */

859 if (rve1 == emin) {

860 irv = STRTOG_Normal | STRTOG_Inexhi;

                                         if (rvb->wds == 1 && rvb->x[0] == 1)

862 sudden_underflow = 1;

863 break;

864 }

865 rve -= nbits;

                                 rvb = set_ones(rvb, rvbits = nbits);

867 break;

868 }

869 else

870 irv = STRTOG_Normal | STRTOG_Inexhi;

                         if ((bbbits < nbits && !denorm) || !(rvb->x[0] & 1))

872 break;

873 if (dsign) {

874 rvb = increment(rvb);

                                 j = kmask & (ULbits - (rvbits & kmask));

                                 if (hi0bits(rvb->x[rvb->wds - 1]) != j)

877 rvbits++;

878 irv = STRTOG_Normal | STRTOG_Inexhi;

879 }

880 else {

                                 if (bbbits == 1)

882 goto undfl;

883 decrement(rvb);

884 irv = STRTOG_Normal | STRTOG_Inexlo;

885 }

886 break;

887 }

                 if ((dval(&adj) = ratio(delta, bs)) <= 2.) {

889 adj1:

890 inex = STRTOG_Inexlo;

891 if (dsign) {

892 asub = 0;

893 inex = STRTOG_Inexhi;

894 }

                         else if (denorm && bbbits <= 1) {

896 undfl:

897 rvb->wds = 0;

898 rve = emin;

899 irv = STRTOG_Underflow | STRTOG_Inexlo;

900 break;

901 }

                         adj0 = dval(&adj) = 1.;

903 }

904 else {

                         adj0 = dval(&adj) *= 0.5;

906 if (dsign) {

907 asub = 0;

908 inex = STRTOG_Inexlo;

909 }

                         if (dval(&adj) < 2147483647.) {

911 L = adj0;

912 adj0 -= L;

913 switch(rd) {

914 case 0:

                                         if (adj0 >= .5)

916 goto inc_L;

917 break;

918 case 1:

                                         if (asub && adj0 > 0.)

920 goto inc_L;

921 break;

922 case 2:

                                         if (!asub && adj0 > 0.) {

924 inc_L:

925 L++;

926 inex = STRTOG_Inexact - inex;

927 }

928 }

929 dval(&adj) = L;

930 }

931 }

932 y = rve + rvbits;

933

934 /* adj *= ulp(dval(&rv)); */

935 /* if (asub) rv -= adj; else rv += adj; */

936

                 if (!denorm && rvbits < nbits) {

                         rvb = lshift(rvb, j = nbits - rvbits);

939 rve -= j;

940 rvbits = nbits;

941 }

                 ab = d2b(dval(&adj), &abe, &abits);

                 if (abe < 0)

944 rshift(ab, -abe);

                 else if (abe > 0)

                         ab = lshift(ab, abe);

947 rvb0 = rvb;

948 if (asub) {

949 /* rv -= adj; */

                         j = hi0bits(rvb->x[rvb->wds-1]);

                         rvb = diff(rvb, ab);

                         k = rvb0->wds - 1;

953 if (denorm)

954 /* do nothing */;

955 else if (rvb->wds <= k

                                 || hi0bits( rvb->x[k]) >

                                    hi0bits(rvb0->x[k])) {

958 /* unlikely; can only have lost 1 high bit */

959 if (rve1 == emin) {

960 --rvbits;

961 denorm = 1;

962 }

963 else {

                                         rvb = lshift(rvb, 1);

965 --rve;

966 --rve1;

967 L = finished = 0;

968 }

969 }

970 }

971 else {

                         rvb = sum(rvb, ab);

                         k = rvb->wds - 1;

974 if (k >= rvb0->wds

                          || hi0bits(rvb->x[k]) < hi0bits(rvb0->x[k])) {

976 if (denorm) {

977 if (++rvbits == nbits)

978 denorm = 0;

979 }

980 else {

                                         rshift(rvb, 1);

982 rve++;

983 rve1++;

984 L = 0;

985 }

986 }

987 }

988 Bfree(ab);

989 Bfree(rvb0);

990 if (finished)

991 break;

992

993 z = rve + rvbits;

                 if (y == z && L) {

995 /* Can we stop now? */

                         tol = dval(&adj) * 5e-16; /* > max rel error */

                         dval(&adj) = adj0 - .5;

                         if (dval(&adj) < -tol) {

999 if (adj0 > tol) {

1000 irv |= inex;

1001 break;

1002 }

1003 }

                         else if (dval(&adj) > tol && adj0 < 1. - tol) {

1005 irv |= inex;

1006 break;

1007 }

1008 }

                 bb0 = denorm ? 0 : trailz(rvb);

1010 Bfree(bb);

1011 Bfree(bd);

1012 Bfree(bs);

1013 Bfree(delta);

1014 }

         if (!denorm && (j = nbits - rvbits)) {

                 if (j > 0)

                         rvb = lshift(rvb, j);

1018 else

1019 rshift(rvb, -j);

1020 rve -= j;

1021 }

1022 *exp = rve;

1023 Bfree(bb);

1024 Bfree(bd);

1025 Bfree(bs);

1026 Bfree(bd0);

1027 Bfree(delta);

         if (rve > fpi->emax) {

                 switch(fpi->rounding & 3) {

1030 case FPI_Round_near:

1031 goto huge;

1032 case FPI_Round_up:

1033 if (!sign)

1034 goto huge;

1035 break;

1036 case FPI_Round_down:

1037 if (sign)

1038 goto huge;

1039 }

1040 /* Round to largest representable magnitude */

1041 Bfree(rvb);

1042 rvb = 0;

1043 irv = STRTOG_Normal | STRTOG_Inexlo;

1044 *exp = fpi->emax;

1045 b = bits;

                 be = b + ((fpi->nbits + 31) >> 5);

1047 while(b < be)

1048 *b++ = -1;

                 if ((j = fpi->nbits & 0x1f))

                         *--be >>= (32 - j);

1051 goto ret;

1052 huge:

1053 rvb->wds = 0;

1054 irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;

1055 #ifndef NO_ERRNO

1056 errno = ERANGE;

1057 #endif

1058 infnanexp:

                 *exp = fpi->emax + 1;

1060 }

1061 ret:

1062 if (denorm) {

1063 if (sudden_underflow) {

1064 rvb->wds = 0;

1065 irv = STRTOG_Underflow | STRTOG_Inexlo;

1066 #if !defined(NO_ERRNO) && __DARWIN_UNIX03

1067 errno = ERANGE;

1068 #endif

1069 }

1070 else {

1071 irv = (irv & ~STRTOG_Retmask) |

                                 (rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero);

1073 if (irv & STRTOG_Inexact) {

1074 irv |= STRTOG_Underflow;

1075 #if !defined(NO_ERRNO) && __DARWIN_UNIX03

1076 errno = ERANGE;

1077 #endif

1078 }

1079 }

1080 }

1081 if (se)

                 *se = (char *)s;

1083 if (sign)

1084 irv |= STRTOG_Neg;

1085 if (rvb) {

                 copybits(bits, nbits, rvb);

1087 Bfree(rvb);

1088 }

1089 if (strunc)

1090 #ifdef FREE

1091 FREE(strunc);

1092 #else

1093 free(strunc);

1094 #endif

1095 return irv;

1096 }