git.saurik.com Git - apple/libc.git/blob - gdtoa/FreeBSD/gdtoa-gdtoa.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 "gdtoaimp.h"

34 static Bigint *

35 #ifdef KR_headers

 bitstob(bits, nbits, bbits) ULong *bits; int nbits; int *bbits;

37 #else

 bitstob(ULong *bits, int nbits, int *bbits)

39 #endif

40 {

41 int i, k;

42 Bigint *b;

43 ULong *be, *x, *x0;

45 i = ULbits;

46 k = 0;

47 while(i < nbits) {

48 i <<= 1;

49 k++;

50 }

51 #ifndef Pack_32

52 if (!k)

53 k = 1;

54 #endif

55 b = Balloc(k);

         be = bits + ((nbits - 1) >> kshift);

57 x = x0 = b->x;

58 do {

59 *x++ = *bits & ALL_ON;

60 #ifdef Pack_16

                 *x++ = (*bits >> 16) & ALL_ON;

62 #endif

                 } while(++bits <= be);

64 i = x - x0;

65 while(!x0[--i])

66 if (!i) {

67 b->wds = 0;

68 *bbits = 0;

69 goto ret;

70 }

         b->wds = i + 1;

         *bbits = i*ULbits + 32 - hi0bits(b->x[i]);

73 ret:

74 return b;

75 }

77 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.

78 *

79 * Inspired by "How to Print Floating-Point Numbers Accurately" by

80 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].

81 *

82 * Modifications:

83 * 1. Rather than iterating, we use a simple numeric overestimate

84 * to determine k = floor(log10(d)). We scale relevant

85 * quantities using O(log2(k)) rather than O(k) multiplications.

86 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't

87 * try to generate digits strictly left to right. Instead, we

88 * compute with fewer bits and propagate the carry if necessary

89 * when rounding the final digit up. This is often faster.

90 * 3. Under the assumption that input will be rounded nearest,

91 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.

92 * That is, we allow equality in stopping tests when the

93 * round-nearest rule will give the same floating-point value

94 * as would satisfaction of the stopping test with strict

95 * inequality.

96 * 4. We remove common factors of powers of 2 from relevant

97 * quantities.

98 * 5. When converting floating-point integers less than 1e16,

99 * we use floating-point arithmetic rather than resorting

100 * to multiple-precision integers.

101 * 6. When asked to produce fewer than 15 digits, we first try

102 * to get by with floating-point arithmetic; we resort to

103 * multiple-precision integer arithmetic only if we cannot

104 * guarantee that the floating-point calculation has given

105 * the correctly rounded result. For k requested digits and

106 * "uniformly" distributed input, the probability is

107 * something like 10^(k-15) that we must resort to the Long

108 * calculation.

109 */

110

111 char *

112 gdtoa

113 #ifdef KR_headers

         (fpi, be, bits, kindp, mode, ndigits, decpt, rve)

         FPI *fpi; int be; ULong *bits;

         int *kindp, mode, ndigits, *decpt; char **rve;

117 #else

         (FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, int *decpt, char **rve)

119 #endif

120 {

121 /* Arguments ndigits and decpt are similar to the second and third

122 arguments of ecvt and fcvt; trailing zeros are suppressed from

123 the returned string. If not null, *rve is set to point

124 to the end of the return value. If d is +-Infinity or NaN,

125 then *decpt is set to 9999.

126

127 mode:

128 0 ==> shortest string that yields d when read in

129 and rounded to nearest.

130 1 ==> like 0, but with Steele & White stopping rule;

131 e.g. with IEEE P754 arithmetic , mode 0 gives

132 1e23 whereas mode 1 gives 9.999999999999999e22.

133 2 ==> max(1,ndigits) significant digits. This gives a

134 return value similar to that of ecvt, except

135 that trailing zeros are suppressed.

136 3 ==> through ndigits past the decimal point. This

137 gives a return value similar to that from fcvt,

138 except that trailing zeros are suppressed, and

139 ndigits can be negative.

140 4-9 should give the same return values as 2-3, i.e.,

141 4 <= mode <= 9 ==> same return as mode

142 2 + (mode & 1). These modes are mainly for

143 debugging; often they run slower but sometimes

144 faster than modes 2-3.

145 4,5,8,9 ==> left-to-right digit generation.

146 6-9 ==> don't try fast floating-point estimate

147 (if applicable).

148

149 Values of mode other than 0-9 are treated as mode 0.

150

151 Sufficient space is allocated to the return value

152 to hold the suppressed trailing zeros.

153 */

154

         int bbits, b2, b5, be0, dig, i, ieps, ilim, ilim0, ilim1, inex;

         int j, j1, k, k0, k_check, kind, leftright, m2, m5, nbits;

         int rdir, s2, s5, spec_case, try_quick;

158 Long L;

         Bigint *b, *b1, *delta, *mlo, *mhi, *mhi1, *S;

160 double d2, ds;

161 char *s, *s0;

162 U d, eps;

163

164 #ifndef MULTIPLE_THREADS

165 if (dtoa_result) {

166 freedtoa(dtoa_result);

167 dtoa_result = 0;

168 }

169 #endif

170 inex = 0;

171 kind = *kindp &= ~STRTOG_Inexact;

172 switch(kind & STRTOG_Retmask) {

173 case STRTOG_Zero:

174 goto ret_zero;

175 case STRTOG_Normal:

176 case STRTOG_Denormal:

177 break;

178 case STRTOG_Infinite:

179 *decpt = -32768;

                 return nrv_alloc("Infinity", rve, 8);

181 case STRTOG_NaN:

182 *decpt = -32768;

                 return nrv_alloc("NaN", rve, 3);

184 default:

185 return 0;

186 }

         b = bitstob(bits, nbits = fpi->nbits, &bbits);

188 be0 = be;

         if ( (i = trailz(b)) !=0) {

190 rshift(b, i);

191 be += i;

192 bbits -= i;

193 }

194 if (!b->wds) {

195 Bfree(b);

196 ret_zero:

197 *decpt = 1;

                 return nrv_alloc("0", rve, 1);

199 }

200

         dval(&d) = b2d(b, &i);

         i = be + bbits - 1;

203 word0(&d) &= Frac_mask1;

204 word0(&d) |= Exp_11;

205 #ifdef IBM

         if ( (j = 11 - hi0bits(word0(&d) & Frac_mask)) !=0)

                 dval(&d) /= 1 << j;

208 #endif

209

210 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5

211 * log10(x) = log(x) / log(10)

212 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))

213 * log10(&d) = (i-Bias)*log(2)/log(10) + log10(d2)

214 *

215 * This suggests computing an approximation k to log10(&d) by

216 *

217 * k = (i - Bias)*0.301029995663981

218 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );

219 *

220 * We want k to be too large rather than too small.

221 * The error in the first-order Taylor series approximation

222 * is in our favor, so we just round up the constant enough

223 * to compensate for any error in the multiplication of

224 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,

225 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,

226 * adding 1e-13 to the constant term more than suffices.

227 * Hence we adjust the constant term to 0.1760912590558.

228 * (We could get a more accurate k by invoking log10,

229 * but this is probably not worthwhile.)

230 */

231 #ifdef IBM

232 i <<= 2;

233 i += j;

234 #endif

         ds = (dval(&d)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;

236

237 /* correct assumption about exponent range */

         if ((j = i) < 0)

239 j = -j;

         if ((j -= 1077) > 0)

241 ds += j * 7e-17;

242

243 k = (int)ds;

         if (ds < 0. && ds != k)

245 k--; /* want k = floor(ds) */

246 k_check = 1;

247 #ifdef IBM

         j = be + bbits - 1;

         if ( (j1 = j & 3) !=0)

                 dval(&d) *= 1 << j1;

         word0(&d) += j << Exp_shift - 2 & Exp_mask;

252 #else

         word0(&d) += (be + bbits - 1) << Exp_shift;

254 #endif

         if (k >= 0 && k <= Ten_pmax) {

                 if (dval(&d) < tens[k])

257 k--;

258 k_check = 0;

259 }

         j = bbits - i - 1;

         if (j >= 0) {

262 b2 = 0;

263 s2 = j;

264 }

265 else {

266 b2 = -j;

267 s2 = 0;

268 }

         if (k >= 0) {

270 b5 = 0;

271 s5 = k;

272 s2 += k;

273 }

274 else {

275 b2 -= k;

276 b5 = -k;

277 s5 = 0;

278 }

         if (mode < 0 || mode > 9)

280 mode = 0;

281 try_quick = 1;

         if (mode > 5) {

283 mode -= 4;

284 try_quick = 0;

285 }

286 leftright = 1;

         ilim = ilim1 = -1;      /* Values for cases 0 and 1; done here to */

288 /* silence erroneous "gcc -Wall" warning. */

289 switch(mode) {

290 case 0:

291 case 1:

                         i = (int)(nbits * .30103) + 3;

293 ndigits = 0;

294 break;

295 case 2:

296 leftright = 0;

297 /* no break */

298 case 4:

                         if (ndigits <= 0)

300 ndigits = 1;

301 ilim = ilim1 = i = ndigits;

302 break;

303 case 3:

304 leftright = 0;

305 /* no break */

306 case 5:

                         i = ndigits + k + 1;

308 ilim = i;

309 ilim1 = i - 1;

                         if (i <= 0)

311 i = 1;

312 }

         s = s0 = rv_alloc(i);

314

         if ( (rdir = fpi->rounding - 1) !=0) {

                 if (rdir < 0)

317 rdir = 2;

318 if (kind & STRTOG_Neg)

319 rdir = 3 - rdir;

320 }

321

322 /* Now rdir = 0 ==> round near, 1 ==> round up, 2 ==> round down. */

323

         if (ilim >= 0 && ilim <= Quick_max && try_quick && !rdir

325 #ifndef IMPRECISE_INEXACT

326 && k == 0

327 #endif

328 ) {

329

330 /* Try to get by with floating-point arithmetic. */

331

332 i = 0;

333 d2 = dval(&d);

334 #ifdef IBM

                 if ( (j = 11 - hi0bits(word0(&d) & Frac_mask)) !=0)

                         dval(&d) /= 1 << j;

337 #endif

338 k0 = k;

339 ilim0 = ilim;

340 ieps = 2; /* conservative */

                 if (k > 0) {

                         ds = tens[k&0xf];

343 j = k >> 4;

344 if (j & Bletch) {

345 /* prevent overflows */

346 j &= Bletch - 1;

                                 dval(&d) /= bigtens[n_bigtens-1];

348 ieps++;

349 }

                         for(; j; j >>= 1, i++)

                                 if (j & 1) {

352 ieps++;

353 ds *= bigtens[i];

354 }

355 }

356 else {

357 ds = 1.;

                         if ( (j1 = -k) !=0) {

                                 dval(&d) *= tens[j1 & 0xf];

                                 for(j = j1 >> 4; j; j >>= 1, i++)

                                         if (j & 1) {

362 ieps++;

                                                 dval(&d) *= bigtens[i];

364 }

365 }

366 }

                 if (k_check && dval(&d) < 1. && ilim > 0) {

                         if (ilim1 <= 0)

369 goto fast_failed;

370 ilim = ilim1;

371 k--;

                         dval(&d) *= 10.;

373 ieps++;

374 }

                 dval(&eps) = ieps*dval(&d) + 7.;

                 word0(&eps) -= (P-1)*Exp_msk1;

                 if (ilim == 0) {

378 S = mhi = 0;

                         dval(&d) -= 5.;

                         if (dval(&d) > dval(&eps))

381 goto one_digit;

                         if (dval(&d) < -dval(&eps))

383 goto no_digits;

384 goto fast_failed;

385 }

386 #ifndef No_leftright

387 if (leftright) {

388 /* Use Steele & White method of only

389 * generating digits needed.

390 */

                         dval(&eps) = ds*0.5/tens[ilim-1] - dval(&eps);

                         for(i = 0;;) {

                                 L = (Long)(dval(&d)/ds);

                                 dval(&d) -= L*ds;

                                 *s++ = '0' + (int)L;

                                 if (dval(&d) < dval(&eps)) {

                                         if (dval(&d))

398 inex = STRTOG_Inexlo;

399 goto ret1;

400 }

                                 if (ds - dval(&d) < dval(&eps))

402 goto bump_up;

403 if (++i >= ilim)

404 break;

                                 dval(&eps) *= 10.;

                                 dval(&d) *= 10.;

407 }

408 }

409 else {

410 #endif

411 /* Generate ilim digits, then fix them up. */

                         dval(&eps) *= tens[ilim-1];

                         for(i = 1;; i++, dval(&d) *= 10.) {

                                 if ( (L = (Long)(dval(&d)/ds)) !=0)

                                         dval(&d) -= L*ds;

                                 *s++ = '0' + (int)L;

417 if (i == ilim) {

418 ds *= 0.5;

                                         if (dval(&d) > ds + dval(&eps))

420 goto bump_up;

                                         else if (dval(&d) < ds - dval(&eps)) {

                                                 if (dval(&d))

423 inex = STRTOG_Inexlo;

424 goto clear_trailing0;

425 }

426 break;

427 }

428 }

429 #ifndef No_leftright

430 }

431 #endif

432 fast_failed:

433 s = s0;

434 dval(&d) = d2;

435 k = k0;

436 ilim = ilim0;

437 }

438

439 /* Do we have a "small" integer? */

440

         if (be >= 0 && k <= Int_max) {

442 /* Yes. */

443 ds = tens[k];

                 if (ndigits < 0 && ilim <= 0) {

445 S = mhi = 0;

                         if (ilim < 0 || dval(&d) <= 5*ds)

447 goto no_digits;

448 goto one_digit;

449 }

                 for(i = 1;; i++, dval(&d) *= 10.) {

                         L = dval(&d) / ds;

                         dval(&d) -= L*ds;

453 #ifdef Check_FLT_ROUNDS

454 /* If FLT_ROUNDS == 2, L will usually be high by 1 */

                         if (dval(&d) < 0) {

456 L--;

457 dval(&d) += ds;

458 }

459 #endif

                         *s++ = '0' + (int)L;

                         if (dval(&d) == 0.)

462 break;

463 if (i == ilim) {

464 if (rdir) {

                                         if (rdir == 1)

466 goto bump_up;

467 inex = STRTOG_Inexlo;

468 goto ret1;

469 }

                                 dval(&d) += dval(&d);

                                 if (dval(&d) > ds || (dval(&d) == ds && L & 1)) {

472 bump_up:

473 inex = STRTOG_Inexhi;

                                         while(*--s == '9')

475 if (s == s0) {

476 k++;

477 *s = '0';

478 break;

479 }

480 ++*s++;

481 }

482 else {

483 inex = STRTOG_Inexlo;

484 clear_trailing0:

                                         while(*--s == '0'){}

486 ++s;

487 }

488 break;

489 }

490 }

491 goto ret1;

492 }

493

494 m2 = b2;

495 m5 = b5;

496 mhi = mlo = 0;

497 if (leftright) {

                 if (mode < 2) {

499 i = nbits - bbits;

                         if (be - i++ < fpi->emin)

501 /* denormal */

                                 i = be - fpi->emin + 1;

503 }

504 else {

505 j = ilim - 1;

506 if (m5 >= j)

507 m5 -= j;

508 else {

509 s5 += j -= m5;

510 b5 += j;

511 m5 = 0;

512 }

                         if ((i = ilim) < 0) {

514 m2 -= i;

515 i = 0;

516 }

517 }

518 b2 += i;

519 s2 += i;

                 mhi = i2b(1);

521 }

         if (m2 > 0 && s2 > 0) {

                 i = m2 < s2 ? m2 : s2;

524 b2 -= i;

525 m2 -= i;

526 s2 -= i;

527 }

         if (b5 > 0) {

529 if (leftright) {

                         if (m5 > 0) {

                                 mhi = pow5mult(mhi, m5);

                                 b1 = mult(mhi, b);

533 Bfree(b);

534 b = b1;

535 }

                         if ( (j = b5 - m5) !=0)

                                 b = pow5mult(b, j);

538 }

539 else

                         b = pow5mult(b, b5);

541 }

         S = i2b(1);

         if (s5 > 0)

                 S = pow5mult(S, s5);

545

546 /* Check for special case that d is a normalized power of 2. */

547

548 spec_case = 0;

         if (mode < 2) {

                 if (bbits == 1 && be0 > fpi->emin + 1) {

551 /* The special case */

552 b2++;

553 s2++;

554 spec_case = 1;

555 }

556 }

557

558 /* Arrange for convenient computation of quotients:

559 * shift left if necessary so divisor has 4 leading 0 bits.

560 *

561 * Perhaps we should just compute leading 28 bits of S once

562 * and for all and pass them and a shift to quorem, so it

563 * can do shifts and ors to compute the numerator for q.

564 */

         i = ((s5 ? hi0bits(S->x[S->wds-1]) : ULbits - 1) - s2 - 4) & kmask;

566 m2 += i;

         if ((b2 += i) > 0)

                 b = lshift(b, b2);

         if ((s2 += i) > 0)

                 S = lshift(S, s2);

571 if (k_check) {

                 if (cmp(b,S) < 0) {

573 k--;

                         b = multadd(b, 10, 0);  /* we botched the k estimate */

575 if (leftright)

                                 mhi = multadd(mhi, 10, 0);

577 ilim = ilim1;

578 }

579 }

         if (ilim <= 0 && mode > 2) {

                 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {

582 /* no digits, fcvt style */

583 no_digits:

584 k = -1 - ndigits;

585 inex = STRTOG_Inexlo;

586 goto ret;

587 }

588 one_digit:

589 inex = STRTOG_Inexhi;

590 *s++ = '1';

591 k++;

592 goto ret;

593 }

594 if (leftright) {

                 if (m2 > 0)

                         mhi = lshift(mhi, m2);

597

598 /* Compute mlo -- check for special case

599 * that d is a normalized power of 2.

600 */

601

602 mlo = mhi;

603 if (spec_case) {

                         mhi = Balloc(mhi->k);

605 Bcopy(mhi, mlo);

                         mhi = lshift(mhi, 1);

607 }

608

                 for(i = 1;;i++) {

                         dig = quorem(b,S) + '0';

611 /* Do we yet have the shortest decimal string

612 * that will round to d?

613 */

                         j = cmp(b, mlo);

                         delta = diff(S, mhi);

                         j1 = delta->sign ? 1 : cmp(b, delta);

617 Bfree(delta);

618 #ifndef ROUND_BIASED

                         if (j1 == 0 && !mode && !(bits[0] & 1) && !rdir) {

                                 if (dig == '9')

621 goto round_9_up;

                                 if (j <= 0) {

                                         if (b->wds > 1 || b->x[0])

624 inex = STRTOG_Inexlo;

625 }

626 else {

627 dig++;

628 inex = STRTOG_Inexhi;

629 }

630 *s++ = dig;

631 goto ret;

632 }

633 #endif

                         if (j < 0 || (j == 0 && !mode

635 #ifndef ROUND_BIASED

                                                         && !(bits[0] & 1)

637 #endif

638 )) {

                                 if (rdir && (b->wds > 1 || b->x[0])) {

                                         if (rdir == 2) {

641 inex = STRTOG_Inexlo;

642 goto accept;

643 }

                                         while (cmp(S,mhi) > 0) {

645 *s++ = dig;

                                                 mhi1 = multadd(mhi, 10, 0);

647 if (mlo == mhi)

648 mlo = mhi1;

649 mhi = mhi1;

                                                 b = multadd(b, 10, 0);

                                                 dig = quorem(b,S) + '0';

652 }

                                         if (dig++ == '9')

654 goto round_9_up;

655 inex = STRTOG_Inexhi;

656 goto accept;

657 }

                                 if (j1 > 0) {

                                         b = lshift(b, 1);

                                         j1 = cmp(b, S);

                                         if ((j1 > 0 || (j1 == 0 && dig & 1))

662 && dig++ == '9')

663 goto round_9_up;

664 inex = STRTOG_Inexhi;

665 }

                                 if (b->wds > 1 || b->x[0])

667 inex = STRTOG_Inexlo;

668 accept:

669 *s++ = dig;

670 goto ret;

671 }

                         if (j1 > 0 && rdir != 2) {

                                 if (dig == '9') { /* possible if i == 1 */

674 round_9_up:

675 *s++ = '9';

676 inex = STRTOG_Inexhi;

677 goto roundoff;

678 }

679 inex = STRTOG_Inexhi;

                                 *s++ = dig + 1;

681 goto ret;

682 }

683 *s++ = dig;

684 if (i == ilim)

685 break;

                         b = multadd(b, 10, 0);

687 if (mlo == mhi)

                                 mlo = mhi = multadd(mhi, 10, 0);

689 else {

                                 mlo = multadd(mlo, 10, 0);

                                 mhi = multadd(mhi, 10, 0);

692 }

693 }

694 }

695 else

                 for(i = 1;; i++) {

                         *s++ = dig = quorem(b,S) + '0';

698 if (i >= ilim)

699 break;

                         b = multadd(b, 10, 0);

701 }

702

703 /* Round off last digit */

704

705 if (rdir) {

                 if (rdir == 2 || (b->wds <= 1 && !b->x[0]))

707 goto chopzeros;

708 goto roundoff;

709 }

         b = lshift(b, 1);

         j = cmp(b, S);

         if (j > 0 || (j == 0 && dig & 1)) {

713 roundoff:

714 inex = STRTOG_Inexhi;

                 while(*--s == '9')

716 if (s == s0) {

717 k++;

718 *s++ = '1';

719 goto ret;

720 }

721 ++*s++;

722 }

723 else {

724 chopzeros:

                 if (b->wds > 1 || b->x[0])

726 inex = STRTOG_Inexlo;

                 while(*--s == '0'){}

728 ++s;

729 }

730 ret:

731 Bfree(S);

732 if (mhi) {

                 if (mlo && mlo != mhi)

734 Bfree(mlo);

735 Bfree(mhi);

736 }

737 ret1:

738 Bfree(b);

739 *s = 0;

         *decpt = k + 1;

741 if (rve)

742 *rve = s;

743 *kindp |= inex;

744 return s0;

745 }