src/third_party/dtoa/dtoa.c - Issue 5195003: Remove Gay's dtoa from sources.

Side by Side Diff: src/third_party/dtoa/dtoa.c

Issue 5195003: Remove Gay's dtoa from sources. (Closed) Base URL: https://v8.googlecode.com/svn/branches/bleeding_edge

Patch Set: Added files to ARM and x64 VS project files. Created 10 years ago

Use n/p to move between diff chunks; N/P to move between comments. Draft comments are only viewable by you.

Jump to:

View unified diff | Download patch | Annotate | Revision Log

OLD	NEW
	(Empty)
1 /****************************************************************

2 *

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

4 *

5 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.

6 *

7 * Permission to use, copy, modify, and distribute this software for any

8 * purpose without fee is hereby granted, provided that this entire notice

9 * is included in all copies of any software which is or includes a copy

10 * or modification of this software and in all copies of the supporting

11 * documentation for such software.

12 *

13 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED

14 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY

15 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY

16 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.

17 *

18 ***************************************************************/

19

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

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

22

23 /* On a machine with IEEE extended-precision registers, it is

24 * necessary to specify double-precision (53-bit) rounding precision

25 * before invoking strtod or dtoa. If the machine uses (the equivalent

26 * of) Intel 80x87 arithmetic, the call

27 * _control87(PC_53, MCW_PC);

28 * does this with many compilers. Whether this or another call is

29 * appropriate depends on the compiler; for this to work, it may be

30 * necessary to #include "float.h" or another system-dependent header

31 * file.

32 */

33

34 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.

35 *

36 * This strtod returns a nearest machine number to the input decimal

37 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are

38 * broken by the IEEE round-even rule. Otherwise ties are broken by

39 * biased rounding (add half and chop).

40 *

41 * Inspired loosely by William D. Clinger's paper "How to Read Floating

42 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].

43 *

44 * Modifications:

45 *

46 * 1. We only require IEEE, IBM, or VAX double-precision

47 * arithmetic (not IEEE double-extended).

48 * 2. We get by with floating-point arithmetic in a case that

49 * Clinger missed -- when we're computing d * 10^n

50 * for a small integer d and the integer n is not too

51 * much larger than 22 (the maximum integer k for which

52 * we can represent 10^k exactly), we may be able to

53 * compute (d10^k) 10^(e-k) with just one roundoff.

54 * 3. Rather than a bit-at-a-time adjustment of the binary

55 * result in the hard case, we use floating-point

56 * arithmetic to determine the adjustment to within

57 * one bit; only in really hard cases do we need to

58 * compute a second residual.

59 * 4. Because of 3., we don't need a large table of powers of 10

60 * for ten-to-e (just some small tables, e.g. of 10^k

61 * for 0 <= k <= 22).

62 */

63

64 /*

65 * #define IEEE_8087 for IEEE-arithmetic machines where the least

66 * significant byte has the lowest address.

67 * #define IEEE_MC68k for IEEE-arithmetic machines where the most

68 * significant byte has the lowest address.

69 * #define Long int on machines with 32-bit ints and 64-bit longs.

70 * #define IBM for IBM mainframe-style floating-point arithmetic.

71 * #define VAX for VAX-style floating-point arithmetic (D_floating).

72 * #define No_leftright to omit left-right logic in fast floating-point

73 * computation of dtoa.

74 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3

75 * and strtod and dtoa should round accordingly.

76 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3

77 * and Honor_FLT_ROUNDS is not #defined.

78 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines

79 * that use extended-precision instructions to compute rounded

80 * products and quotients) with IBM.

81 * #define ROUND_BIASED for IEEE-format with biased rounding.

82 * #define Inaccurate_Divide for IEEE-format with correctly rounded

83 * products but inaccurate quotients, e.g., for Intel i860.

84 * #define NO_LONG_LONG on machines that do not have a "long long"

85 * integer type (of >= 64 bits). On such machines, you can

86 * #define Just_16 to store 16 bits per 32-bit Long when doing

87 * high-precision integer arithmetic. Whether this speeds things

88 * up or slows things down depends on the machine and the number

89 * being converted. If long long is available and the name is

90 * something other than "long long", #define Llong to be the name,

91 * and if "unsigned Llong" does not work as an unsigned version of

92 * Llong, #define #ULLong to be the corresponding unsigned type.

93 * #define KR_headers for old-style C function headers.

94 * #define Bad_float_h if your system lacks a float.h or if it does not

95 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,

96 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.

97 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)

98 * if memory is available and otherwise does something you deem

99 * appropriate. If MALLOC is undefined, malloc will be invoked

100 * directly -- and assumed always to succeed.

101 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making

102 * memory allocations from a private pool of memory when possible.

103 * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes,

104 * unless #defined to be a different length. This default length

105 * suffices to get rid of MALLOC calls except for unusual cases,

106 * such as decimal-to-binary conversion of a very long string of

107 * digits. The longest string dtoa can return is about 751 bytes

108 * long. For conversions by strtod of strings of 800 digits and

109 * all dtoa conversions in single-threaded executions with 8-byte

110 * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte

111 * pointers, PRIVATE_MEM >= 7112 appears adequate.

112 * #define INFNAN_CHECK on IEEE systems to cause strtod to check for

113 * Infinity and NaN (case insensitively). On some systems (e.g.,

114 * some HP systems), it may be necessary to #define NAN_WORD0

115 * appropriately -- to the most significant word of a quiet NaN.

116 * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)

117 * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,

118 * strtod also accepts (case insensitively) strings of the form

119 * NaN(x), where x is a string of hexadecimal digits and spaces;

120 * if there is only one string of hexadecimal digits, it is taken

121 * for the 52 fraction bits of the resulting NaN; if there are two

122 * or more strings of hex digits, the first is for the high 20 bits,

123 * the second and subsequent for the low 32 bits, with intervening

124 * white space ignored; but if this results in none of the 52

125 * fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0

126 * and NAN_WORD1 are used instead.

127 * #define MULTIPLE_THREADS if the system offers preemptively scheduled

128 * multiple threads. In this case, you must provide (or suitably

129 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed

130 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed

131 * in pow5mult, ensures lazy evaluation of only one copy of high

132 * powers of 5; omitting this lock would introduce a small

133 * probability of wasting memory, but would otherwise be harmless.)

134 * You must also invoke freedtoa(s) to free the value s returned by

135 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.

136 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that

137 * avoids underflows on inputs whose result does not underflow.

138 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format

139 * floating-point numbers and flushes underflows to zero rather

140 * than implementing gradual underflow, then you must also #define

141 * Sudden_Underflow.

142 * #define YES_ALIAS to permit aliasing certain double values with

143 * arrays of ULongs. This leads to slightly better code with

144 * some compilers and was always used prior to 19990916, but it

145 * is not strictly legal and can cause trouble with aggressively

146 * optimizing compilers (e.g., gcc 2.95.1 under -O2).

147 * #define USE_LOCALE to use the current locale's decimal_point value.

148 * #define SET_INEXACT if IEEE arithmetic is being used and extra

149 * computation should be done to set the inexact flag when the

150 * result is inexact and avoid setting inexact when the result

151 * is exact. In this case, dtoa.c must be compiled in

152 * an environment, perhaps provided by #include "dtoa.c" in a

153 * suitable wrapper, that defines two functions,

154 * int get_inexact(void);

155 * void clear_inexact(void);

156 * such that get_inexact() returns a nonzero value if the

157 * inexact bit is already set, and clear_inexact() sets the

158 * inexact bit to 0. When SET_INEXACT is #defined, strtod

159 * also does extra computations to set the underflow and overflow

160 * flags when appropriate (i.e., when the result is tiny and

161 * inexact or when it is a numeric value rounded to +-infinity).

162 * #define NO_ERRNO if strtod should not assign errno = ERANGE when

163 * the result overflows to +-Infinity or underflows to 0.

164 */

165

166 #ifndef Long

167 #if __LP64__

168 #define Long int

169 #else

170 #define Long long

171 #endif

172 #endif

173 #ifndef ULong

174 typedef unsigned Long ULong;

175 #endif

176

177 #ifdef DEBUG

178 #include "stdio.h"

179 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}

180 #endif

181

182 #include "stdlib.h"

183 #include "string.h"

184

185 #ifdef USE_LOCALE

186 #include "locale.h"

187 #endif

188

189 #ifdef MALLOC

190 #ifdef KR_headers

191 extern char *MALLOC();

192 #else

193 extern void *MALLOC(size_t);

194 #endif

195 #else

196 #define MALLOC malloc

197 #endif

198

199 #ifndef Omit_Private_Memory

200 #ifndef PRIVATE_MEM

201 #define PRIVATE_MEM 2304

202 #endif

203 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))

204 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;

205 #endif

206

207 #undef IEEE_Arith

208 #undef Avoid_Underflow

209 #ifdef IEEE_MC68k

210 #define IEEE_Arith

211 #endif

212 #ifdef IEEE_8087

213 #define IEEE_Arith

214 #endif

215

216 #include "errno.h"

217

218 #ifdef Bad_float_h

219

220 #ifdef IEEE_Arith

221 #define DBL_DIG 15

222 #define DBL_MAX_10_EXP 308

223 #define DBL_MAX_EXP 1024

224 #define FLT_RADIX 2

225 #endif /IEEE_Arith/

226

227 #ifdef IBM

228 #define DBL_DIG 16

229 #define DBL_MAX_10_EXP 75

230 #define DBL_MAX_EXP 63

231 #define FLT_RADIX 16

232 #define DBL_MAX 7.2370055773322621e+75

233 #endif

234

235 #ifdef VAX

236 #define DBL_DIG 16

237 #define DBL_MAX_10_EXP 38

238 #define DBL_MAX_EXP 127

239 #define FLT_RADIX 2

240 #define DBL_MAX 1.7014118346046923e+38

241 #endif

242

243 #ifndef LONG_MAX

244 #define LONG_MAX 2147483647

245 #endif

246

247 #else /* ifndef Bad_float_h */

248 #include "float.h"

249 #endif /* Bad_float_h */

250

251 #ifndef __MATH_H__

252 #include "math.h"

253 #endif

254

255 #ifdef __cplusplus

256 extern "C" {

257 #endif

258

259 #ifndef CONST

260 #ifdef KR_headers

261 #define CONST /* blank */

262 #else

263 #define CONST const

264 #endif

265 #endif

266

267 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1

268 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.

269 #endif

270

271 typedef union { double d; ULong L[2]; } U;

272

273 #ifdef IEEE_8087

274 #define word0(x) (x).L[1]

275 #define word1(x) (x).L[0]

276 #else

277 #define word0(x) (x).L[0]

278 #define word1(x) (x).L[1]

279 #endif

280 #define dval(x) (x).d

281

282 /* The following definition of Storeinc is appropriate for MIPS processors.

283 * An alternative that might be better on some machines is

284 * #define Storeinc(a,b,c) (*a++ = b << 16 \| c & 0xffff)

285 */

286 #if defined(IEEE_8087) + defined(VAX)

287 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \

288 ((unsigned short *)a)[0] = (unsigned short)c, a++)

289 #else

290 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \

291 ((unsigned short *)a)[1] = (unsigned short)c, a++)

292 #endif

293

294 /* #define P DBL_MANT_DIG */

295 /* Ten_pmax = floor(Plog(2)/log(5)) /

296 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */

297 /* Quick_max = floor((P-1)log(FLT_RADIX)/log(10) - 1) /

298 /* Int_max = floor(Plog(FLT_RADIX)/log(10) - 1) /

299

300 #ifdef IEEE_Arith

301 #define Exp_shift 20

302 #define Exp_shift1 20

303 #define Exp_msk1 0x100000

304 #define Exp_msk11 0x100000

305 #define Exp_mask 0x7ff00000

306 #define P 53

307 #define Bias 1023

308 #define Emin (-1022)

309 #define Exp_1 0x3ff00000

310 #define Exp_11 0x3ff00000

311 #define Ebits 11

312 #define Frac_mask 0xfffff

313 #define Frac_mask1 0xfffff

314 #define Ten_pmax 22

315 #define Bletch 0x10

316 #define Bndry_mask 0xfffff

317 #define Bndry_mask1 0xfffff

318 #define LSB 1

319 #define Sign_bit 0x80000000

320 #define Log2P 1

321 #define Tiny0 0

322 #define Tiny1 1

323 #define Quick_max 14

324 #define Int_max 14

325 #ifndef NO_IEEE_Scale

326 #define Avoid_Underflow

327 #ifdef Flush_Denorm /* debugging option */

328 #undef Sudden_Underflow

329 #endif

330 #endif

331

332 #ifndef Flt_Rounds

333 #ifdef FLT_ROUNDS

334 #define Flt_Rounds FLT_ROUNDS

335 #else

336 #define Flt_Rounds 1

337 #endif

338 #endif /Flt_Rounds/

339

340 #ifdef Honor_FLT_ROUNDS

341 #define Rounding rounding

342 #undef Check_FLT_ROUNDS

343 #define Check_FLT_ROUNDS

344 #else

345 #define Rounding Flt_Rounds

346 #endif

347

348 #else /* ifndef IEEE_Arith */

349 #undef Check_FLT_ROUNDS

350 #undef Honor_FLT_ROUNDS

351 #undef SET_INEXACT

352 #undef Sudden_Underflow

353 #define Sudden_Underflow

354 #ifdef IBM

355 #undef Flt_Rounds

356 #define Flt_Rounds 0

357 #define Exp_shift 24

358 #define Exp_shift1 24

359 #define Exp_msk1 0x1000000

360 #define Exp_msk11 0x1000000

361 #define Exp_mask 0x7f000000

362 #define P 14

363 #define Bias 65

364 #define Exp_1 0x41000000

365 #define Exp_11 0x41000000

366 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */

367 #define Frac_mask 0xffffff

368 #define Frac_mask1 0xffffff

369 #define Bletch 4

370 #define Ten_pmax 22

371 #define Bndry_mask 0xefffff

372 #define Bndry_mask1 0xffffff

373 #define LSB 1

374 #define Sign_bit 0x80000000

375 #define Log2P 4

376 #define Tiny0 0x100000

377 #define Tiny1 0

378 #define Quick_max 14

379 #define Int_max 15

380 #else /* VAX */

381 #undef Flt_Rounds

382 #define Flt_Rounds 1

383 #define Exp_shift 23

384 #define Exp_shift1 7

385 #define Exp_msk1 0x80

386 #define Exp_msk11 0x800000

387 #define Exp_mask 0x7f80

388 #define P 56

389 #define Bias 129

390 #define Exp_1 0x40800000

391 #define Exp_11 0x4080

392 #define Ebits 8

393 #define Frac_mask 0x7fffff

394 #define Frac_mask1 0xffff007f

395 #define Ten_pmax 24

396 #define Bletch 2

397 #define Bndry_mask 0xffff007f

398 #define Bndry_mask1 0xffff007f

399 #define LSB 0x10000

400 #define Sign_bit 0x8000

401 #define Log2P 1

402 #define Tiny0 0x80

403 #define Tiny1 0

404 #define Quick_max 15

405 #define Int_max 15

406 #endif /* IBM, VAX */

407 #endif /* IEEE_Arith */

408

409 #ifndef IEEE_Arith

410 #define ROUND_BIASED

411 #endif

412

413 #ifdef RND_PRODQUOT

414 #define rounded_product(a,b) a = rnd_prod(a, b)

415 #define rounded_quotient(a,b) a = rnd_quot(a, b)

416 #ifdef KR_headers

417 extern double rnd_prod(), rnd_quot();

418 #else

419 extern double rnd_prod(double, double), rnd_quot(double, double);

420 #endif

421 #else

422 #define rounded_product(a,b) a *= b

423 #define rounded_quotient(a,b) a /= b

424 #endif

425

426 #define Big0 (Frac_mask1 \| Exp_msk1*(DBL_MAX_EXP+Bias-1))

427 #define Big1 0xffffffff

428

429 #ifndef Pack_32

430 #define Pack_32

431 #endif

432

433 #ifdef KR_headers

434 #define FFFFFFFF ((((unsigned long)0xffff)<<16)\|(unsigned long)0xffff)

435 #else

436 #define FFFFFFFF 0xffffffffUL

437 #endif

438

439 #ifdef NO_LONG_LONG

440 #undef ULLong

441 #ifdef Just_16

442 #undef Pack_32

443 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.

444 * This makes some inner loops simpler and sometimes saves work

445 * during multiplications, but it often seems to make things slightly

446 * slower. Hence the default is now to store 32 bits per Long.

447 */

448 #endif

449 #else /* long long available */

450 #ifndef Llong

451 #define Llong long long

452 #endif

453 #ifndef ULLong

454 #define ULLong unsigned Llong

455 #endif

456 #endif /* NO_LONG_LONG */

457

458 #ifndef MULTIPLE_THREADS

459 #define ACQUIRE_DTOA_LOCK(n) /nothing/

460 #define FREE_DTOA_LOCK(n) /nothing/

461 #endif

462

463 #define Kmax 15

464

465 #ifdef __cplusplus

466 extern "C" double strtod(const char s00, char *se);

467 extern "C" char *dtoa(double d, int mode, int ndigits,

468 int decpt, int sign, char **rve);

469 #endif

470

471 struct

472 Bigint {

473 struct Bigint *next;

474 int k, maxwds, sign, wds;

475 ULong x[1];

476 };

477

478 typedef struct Bigint Bigint;

479

480 static Bigint *freelist[Kmax+1];

481

482 static Bigint *

483 Balloc

484 #ifdef KR_headers

485 (k) int k;

486 #else

487 (int k)

488 #endif

489 {

490 int x;

491 Bigint *rv;

492 #ifndef Omit_Private_Memory

493 unsigned int len;

494 #endif

495

496 ACQUIRE_DTOA_LOCK(0);

497 /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0). */

498 /* but this case seems very unlikely. */

499 if (k <= Kmax && (rv = freelist[k])) {

500 freelist[k] = rv->next;

501 }

502 else {

503 x = 1 << k;

504 #ifdef Omit_Private_Memory

505 rv = (Bigint )MALLOC(sizeof(Bigint) + (x-1)sizeof(ULong));

506 #else

507 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1 )

508 /sizeof(double);

509 if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {

510 rv = (Bigint*)pmem_next;

511 pmem_next += len;

512 }

513 else

514 rv = (Bigint)MALLOC(lensizeof(double));

515 #endif

516 rv->k = k;

517 rv->maxwds = x;

518 }

519 FREE_DTOA_LOCK(0);

520 rv->sign = rv->wds = 0;

521 return rv;

522 }

523

524 static void

525 Bfree

526 #ifdef KR_headers

527 (v) Bigint *v;

528 #else

529 (Bigint *v)

530 #endif

531 {

532 if (v) {

533 if (v->k > Kmax)

534 free((void*)v);

535 else {

536 ACQUIRE_DTOA_LOCK(0);

537 v->next = freelist[v->k];

538 freelist[v->k] = v;

539 FREE_DTOA_LOCK(0);

540 }

541 }

542 }

543

544 #define Bcopy(x,y) memcpy((char )&x->sign, (char )&y->sign, \

545 y->wdssizeof(Long) + 2sizeof(int))

546

547 static Bigint *

548 multadd

549 #ifdef KR_headers

550 (b, m, a) Bigint *b; int m, a;

551 #else

552 (Bigint b, int m, int a) / multiply by m and add a */

553 #endif

554 {

555 int i, wds;

556 #ifdef ULLong

557 ULong *x;

558 ULLong carry, y;

559 #else

560 ULong carry, *x, y;

561 #ifdef Pack_32

562 ULong xi, z;

563 #endif

564 #endif

565 Bigint *b1;

566

567 wds = b->wds;

568 x = b->x;

569 i = 0;

570 carry = a;

571 do {

572 #ifdef ULLong

573 y = x (ULLong)m + carry;

574 carry = y >> 32;

575 *x++ = y & FFFFFFFF;

576 #else

577 #ifdef Pack_32

578 xi = *x;

579 y = (xi & 0xffff) * m + carry;

580 z = (xi >> 16) * m + (y >> 16);

581 carry = z >> 16;

582 *x++ = (z << 16) + (y & 0xffff);

583 #else

584 y = x m + carry;

585 carry = y >> 16;

586 *x++ = y & 0xffff;

587 #endif

588 #endif

589 }

590 while(++i < wds);

591 if (carry) {

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

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

594 Bcopy(b1, b);

595 Bfree(b);

596 b = b1;

597 }

598 b->x[wds++] = carry;

599 b->wds = wds;

600 }

601 return b;

602 }

603

604 static Bigint *

605 s2b

606 #ifdef KR_headers

607 (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9;

608 #else

609 (CONST char *s, int nd0, int nd, ULong y9)

610 #endif

611 {

612 Bigint *b;

613 int i, k;

614 Long x, y;

615

616 x = (nd + 8) / 9;

617 for(k = 0, y = 1; x > y; y <<= 1, k++) ;

618 #ifdef Pack_32

619 b = Balloc(k);

620 b->x[0] = y9;

621 b->wds = 1;

622 #else

623 b = Balloc(k+1);

624 b->x[0] = y9 & 0xffff;

625 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;

626 #endif

627

628 i = 9;

629 if (9 < nd0) {

630 s += 9;

631 do b = multadd(b, 10, *s++ - '0');

632 while(++i < nd0);

633 s++;

634 }

635 else

636 s += 10;

637 for(; i < nd; i++)

638 b = multadd(b, 10, *s++ - '0');

639 return b;

640 }

641

642 static int

643 hi0bits

644 #ifdef KR_headers

645 (x) register ULong x;

646 #else

647 (register ULong x)

648 #endif

649 {

650 register int k = 0;

651

652 if (!(x & 0xffff0000)) {

653 k = 16;

654 x <<= 16;

655 }

656 if (!(x & 0xff000000)) {

657 k += 8;

658 x <<= 8;

659 }

660 if (!(x & 0xf0000000)) {

661 k += 4;

662 x <<= 4;

663 }

664 if (!(x & 0xc0000000)) {

665 k += 2;

666 x <<= 2;

667 }

668 if (!(x & 0x80000000)) {

669 k++;

670 if (!(x & 0x40000000))

671 return 32;

672 }

673 return k;

674 }

675

676 static int

677 lo0bits

678 #ifdef KR_headers

679 (y) ULong *y;

680 #else

681 (ULong *y)

682 #endif

683 {

684 register int k;

685 register ULong x = *y;

686

687 if (x & 7) {

688 if (x & 1)

689 return 0;

690 if (x & 2) {

691 *y = x >> 1;

692 return 1;

693 }

694 *y = x >> 2;

695 return 2;

696 }

697 k = 0;

698 if (!(x & 0xffff)) {

699 k = 16;

700 x >>= 16;

701 }

702 if (!(x & 0xff)) {

703 k += 8;

704 x >>= 8;

705 }

706 if (!(x & 0xf)) {

707 k += 4;

708 x >>= 4;

709 }

710 if (!(x & 0x3)) {

711 k += 2;

712 x >>= 2;

713 }

714 if (!(x & 1)) {

715 k++;

716 x >>= 1;

717 if (!x)

718 return 32;

719 }

720 *y = x;

721 return k;

722 }

723

724 static Bigint *

725 i2b

726 #ifdef KR_headers

727 (i) int i;

728 #else

729 (int i)

730 #endif

731 {

732 Bigint *b;

733

734 b = Balloc(1);

735 b->x[0] = i;

736 b->wds = 1;

737 return b;

738 }

739

740 static Bigint *

741 mult

742 #ifdef KR_headers

743 (a, b) Bigint a, b;

744 #else

745 (Bigint a, Bigint b)

746 #endif

747 {

748 Bigint *c;

749 int k, wa, wb, wc;

750 ULong x, xa, xae, xb, xbe, xc, *xc0;

751 ULong y;

752 #ifdef ULLong

753 ULLong carry, z;

754 #else

755 ULong carry, z;

756 #ifdef Pack_32

757 ULong z2;

758 #endif

759 #endif

760

761 if (a->wds < b->wds) {

762 c = a;

763 a = b;

764 b = c;

765 }

766 k = a->k;

767 wa = a->wds;

768 wb = b->wds;

769 wc = wa + wb;

770 if (wc > a->maxwds)

771 k++;

772 c = Balloc(k);

773 for(x = c->x, xa = x + wc; x < xa; x++)

774 *x = 0;

775 xa = a->x;

776 xae = xa + wa;

777 xb = b->x;

778 xbe = xb + wb;

779 xc0 = c->x;

780 #ifdef ULLong

781 for(; xb < xbe; xc0++) {

782 if ((y = *xb++)) {

783 x = xa;

784 xc = xc0;

785 carry = 0;

786 do {

787 z = x++ (ULLong)y + *xc + carry;

788 carry = z >> 32;

789 *xc++ = z & FFFFFFFF;

790 }

791 while(x < xae);

792 *xc = carry;

793 }

794 }

795 #else

796 #ifdef Pack_32

797 for(; xb < xbe; xb++, xc0++) {

798 if (y = *xb & 0xffff) {

799 x = xa;

800 xc = xc0;

801 carry = 0;

802 do {

803 z = (x & 0xffff) y + (*xc & 0xffff) + carry;

804 carry = z >> 16;

805 z2 = (x++ >> 16) y + (*xc >> 16) + carry;

806 carry = z2 >> 16;

807 Storeinc(xc, z2, z);

808 }

809 while(x < xae);

810 *xc = carry;

811 }

812 if (y = *xb >> 16) {

813 x = xa;

814 xc = xc0;

815 carry = 0;

816 z2 = *xc;

817 do {

818 z = (x & 0xffff) y + (*xc >> 16) + carry;

819 carry = z >> 16;

820 Storeinc(xc, z, z2);

821 z2 = (x++ >> 16) y + (*xc & 0xffff) + carry;

822 carry = z2 >> 16;

823 }

824 while(x < xae);

825 *xc = z2;

826 }

827 }

828 #else

829 for(; xb < xbe; xc0++) {

830 if (y = *xb++) {

831 x = xa;

832 xc = xc0;

833 carry = 0;

834 do {

835 z = x++ y + *xc + carry;

836 carry = z >> 16;

837 *xc++ = z & 0xffff;

838 }

839 while(x < xae);

840 *xc = carry;

841 }

842 }

843 #endif

844 #endif

845 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;

846 c->wds = wc;

847 return c;

848 }

849

850 static Bigint *p5s;

851

852 static Bigint *

853 pow5mult

854 #ifdef KR_headers

855 (b, k) Bigint *b; int k;

856 #else

857 (Bigint *b, int k)

858 #endif

859 {

860 Bigint b1, p5, *p51;

861 int i;

862 static int p05[3] = { 5, 25, 125 };

863

864 if ((i = k & 3))

865 b = multadd(b, p05[i-1], 0);

866

867 if (!(k >>= 2))

868 return b;

869 if (!(p5 = p5s)) {

870 /* first time */

871 #ifdef MULTIPLE_THREADS

872 ACQUIRE_DTOA_LOCK(1);

873 if (!(p5 = p5s)) {

874 p5 = p5s = i2b(625);

875 p5->next = 0;

876 }

877 FREE_DTOA_LOCK(1);

878 #else

879 p5 = p5s = i2b(625);

880 p5->next = 0;

881 #endif

882 }

883 for(;;) {

884 if (k & 1) {

885 b1 = mult(b, p5);

886 Bfree(b);

887 b = b1;

888 }

889 if (!(k >>= 1))

890 break;

891 if (!(p51 = p5->next)) {

892 #ifdef MULTIPLE_THREADS

893 ACQUIRE_DTOA_LOCK(1);

894 if (!(p51 = p5->next)) {

895 p51 = p5->next = mult(p5,p5);

896 p51->next = 0;

897 }

898 FREE_DTOA_LOCK(1);

899 #else

900 p51 = p5->next = mult(p5,p5);

901 p51->next = 0;

902 #endif

903 }

904 p5 = p51;

905 }

906 return b;

907 }

908

909 static Bigint *

910 lshift

911 #ifdef KR_headers

912 (b, k) Bigint *b; int k;

913 #else

914 (Bigint *b, int k)

915 #endif

916 {

917 int i, k1, n, n1;

918 Bigint *b1;

919 ULong x, x1, *xe, z;

920

921 #ifdef Pack_32

922 n = k >> 5;

923 #else

924 n = k >> 4;

925 #endif

926 k1 = b->k;

927 n1 = n + b->wds + 1;

928 for(i = b->maxwds; n1 > i; i <<= 1)

929 k1++;

930 b1 = Balloc(k1);

931 x1 = b1->x;

932 for(i = 0; i < n; i++)

933 *x1++ = 0;

934 x = b->x;

935 xe = x + b->wds;

936 #ifdef Pack_32

937 if (k &= 0x1f) {

938 k1 = 32 - k;

939 z = 0;

940 do {

941 x1++ = x << k \| z;

942 z = *x++ >> k1;

943 }

944 while(x < xe);

945 if ((*x1 = z))

946 ++n1;

947 }

948 #else

949 if (k &= 0xf) {

950 k1 = 16 - k;

951 z = 0;

952 do {

953 x1++ = x << k & 0xffff \| z;

954 z = *x++ >> k1;

955 }

956 while(x < xe);

957 if (*x1 = z)

958 ++n1;

959 }

960 #endif

961 else do

962 x1++ = x++;

963 while(x < xe);

964 b1->wds = n1 - 1;

965 Bfree(b);

966 return b1;

967 }

968

969 static int

970 cmp

971 #ifdef KR_headers

972 (a, b) Bigint a, b;

973 #else

974 (Bigint a, Bigint b)

975 #endif

976 {

977 ULong xa, xa0, xb, xb0;

978 int i, j;

979

980 i = a->wds;

981 j = b->wds;

982 #ifdef DEBUG

983 if (i > 1 && !a->x[i-1])

984 Bug("cmp called with a->x[a->wds-1] == 0");

985 if (j > 1 && !b->x[j-1])

986 Bug("cmp called with b->x[b->wds-1] == 0");

987 #endif

988 if (i -= j)

989 return i;

990 xa0 = a->x;

991 xa = xa0 + j;

992 xb0 = b->x;

993 xb = xb0 + j;

994 for(;;) {

995 if (--xa != --xb)

996 return xa < xb ? -1 : 1;

997 if (xa <= xa0)

998 break;

999 }

1000 return 0;

1001 }

1002

1003 static Bigint *

1004 diff

1005 #ifdef KR_headers

1006 (a, b) Bigint a, b;

1007 #else

1008 (Bigint a, Bigint b)

1009 #endif

1010 {

1011 Bigint *c;

1012 int i, wa, wb;

1013 ULong xa, xae, xb, xbe, *xc;

1014 #ifdef ULLong

1015 ULLong borrow, y;

1016 #else

1017 ULong borrow, y;

1018 #ifdef Pack_32

1019 ULong z;

1020 #endif

1021 #endif

1022

1023 i = cmp(a,b);

1024 if (!i) {

1025 c = Balloc(0);

1026 c->wds = 1;

1027 c->x[0] = 0;

1028 return c;

1029 }

1030 if (i < 0) {

1031 c = a;

1032 a = b;

1033 b = c;

1034 i = 1;

1035 }

1036 else

1037 i = 0;

1038 c = Balloc(a->k);

1039 c->sign = i;

1040 wa = a->wds;

1041 xa = a->x;

1042 xae = xa + wa;

1043 wb = b->wds;

1044 xb = b->x;

1045 xbe = xb + wb;

1046 xc = c->x;

1047 borrow = 0;

1048 #ifdef ULLong

1049 do {

1050 y = (ULLong)xa++ - xb++ - borrow;

1051 borrow = y >> 32 & (ULong)1;

1052 *xc++ = y & FFFFFFFF;

1053 }

1054 while(xb < xbe);

1055 while(xa < xae) {

1056 y = *xa++ - borrow;

1057 borrow = y >> 32 & (ULong)1;

1058 *xc++ = y & FFFFFFFF;

1059 }

1060 #else

1061 #ifdef Pack_32

1062 do {

1063 y = (xa & 0xffff) - (xb & 0xffff) - borrow;

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

1065 z = (xa++ >> 16) - (xb++ >> 16) - borrow;

1066 borrow = (z & 0x10000) >> 16;

1067 Storeinc(xc, z, y);

1068 }

1069 while(xb < xbe);

1070 while(xa < xae) {

1071 y = (*xa & 0xffff) - borrow;

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

1073 z = (*xa++ >> 16) - borrow;

1074 borrow = (z & 0x10000) >> 16;

1075 Storeinc(xc, z, y);

1076 }

1077 #else

1078 do {

1079 y = xa++ - xb++ - borrow;

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

1081 *xc++ = y & 0xffff;

1082 }

1083 while(xb < xbe);

1084 while(xa < xae) {

1085 y = *xa++ - borrow;

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

1087 *xc++ = y & 0xffff;

1088 }

1089 #endif

1090 #endif

1091 while(!*--xc)

1092 wa--;

1093 c->wds = wa;

1094 return c;

1095 }

1096

1097 static double

1098 ulp

1099 #ifdef KR_headers

1100 (dx) double dx;

1101 #else

1102 (double dx)

1103 #endif

1104 {

1105 register Long L;

1106 U x, a;

1107

1108 dval(x) = dx;

1109

1110 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;

1111 #ifndef Avoid_Underflow

1112 #ifndef Sudden_Underflow

1113 if (L > 0) {

1114 #endif

1115 #endif

1116 #ifdef IBM

1117 L \|= Exp_msk1 >> 4;

1118 #endif

1119 word0(a) = L;

1120 word1(a) = 0;

1121 #ifndef Avoid_Underflow

1122 #ifndef Sudden_Underflow

1123 }

1124 else {

1125 L = -L >> Exp_shift;

1126 if (L < Exp_shift) {

1127 word0(a) = 0x80000 >> L;

1128 word1(a) = 0;

1129 }

1130 else {

1131 word0(a) = 0;

1132 L -= Exp_shift;

1133 word1(a) = L >= 31 ? 1 : 1 << 31 - L;

1134 }

1135 }

1136 #endif

1137 #endif

1138 return dval(a);

1139 }

1140

1141 static double

1142 b2d

1143 #ifdef KR_headers

1144 (a, e) Bigint a; int e;

1145 #else

1146 (Bigint a, int e)

1147 #endif

1148 {

1149 ULong xa, xa0, w, y, z;

1150 int k;

1151 U d;

1152 #ifdef VAX

1153 ULong d0, d1;

1154 #else

1155 #define d0 word0(d)

1156 #define d1 word1(d)

1157 #endif

1158

1159 xa0 = a->x;

1160 xa = xa0 + a->wds;

1161 y = *--xa;

1162 #ifdef DEBUG

1163 if (!y) Bug("zero y in b2d");

1164 #endif

1165 k = hi0bits(y);

1166 *e = 32 - k;

1167 #ifdef Pack_32

1168 if (k < Ebits) {

1169 d0 = Exp_1 \| (y >> (Ebits - k));

1170 w = xa > xa0 ? *--xa : 0;

1171 d1 = (y << ((32-Ebits) + k)) \| (w >> (Ebits - k));

1172 goto ret_d;

1173 }

1174 z = xa > xa0 ? *--xa : 0;

1175 if (k -= Ebits) {

1176 d0 = Exp_1 \| (y << k) \| (z >> (32 - k));

1177 y = xa > xa0 ? *--xa : 0;

1178 d1 = (z << k) \| (y >> (32 - k));

1179 }

1180 else {

1181 d0 = Exp_1 \| y;

1182 d1 = z;

1183 }

1184 #else

1185 if (k < Ebits + 16) {

1186 z = xa > xa0 ? *--xa : 0;

1187 d0 = Exp_1 \| (y << (k - Ebits)) \| (z >> (Ebits + 16 - k));

1188 w = xa > xa0 ? *--xa : 0;

1189 y = xa > xa0 ? *--xa : 0;

1190 d1 = (z << (k + 16 - Ebits)) \| (w << (k - Ebits)) \| (y >> (16 + Ebits - k));

1191 goto ret_d;

1192 }

1193 z = xa > xa0 ? *--xa : 0;

1194 w = xa > xa0 ? *--xa : 0;

1195 k -= Ebits + 16;

1196 d0 = Exp_1 \| y << k + 16 \| z << k \| w >> 16 - k;

1197 y = xa > xa0 ? *--xa : 0;

1198 d1 = w << k + 16 \| y << k;

1199 #endif

1200 ret_d:

1201 #ifdef VAX

1202 word0(d) = d0 >> 16 \| d0 << 16;

1203 word1(d) = d1 >> 16 \| d1 << 16;

1204 #else

1205 #undef d0

1206 #undef d1

1207 #endif

1208 return dval(d);

1209 }

1210

1211 static Bigint *

1212 d2b

1213 #ifdef KR_headers

1214 (dd, e, bits) double dd; int e, bits;

1215 #else

1216 (double dd, int e, int bits)

1217 #endif

1218 {

1219 Bigint *b;

1220 int de, k;

1221 ULong *x, y, z;

1222 #ifndef Sudden_Underflow

1223 int i;

1224 #endif

1225 #ifdef VAX

1226 ULong d0, d1;

1227 d0 = word0(d) >> 16 \| word0(d) << 16;

1228 d1 = word1(d) >> 16 \| word1(d) << 16;

1229 #else

1230 U d;

1231 dval(d) = dd;

1232 #define d0 word0(d)

1233 #define d1 word1(d)

1234 #endif

1235

1236 #ifdef Pack_32

1237 b = Balloc(1);

1238 #else

1239 b = Balloc(2);

1240 #endif

1241 x = b->x;

1242

1243 z = d0 & Frac_mask;

1244 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */

1245 #ifdef Sudden_Underflow

1246 de = (int)(d0 >> Exp_shift);

1247 #ifndef IBM

1248 z \|= Exp_msk11;

1249 #endif

1250 #else

1251 if ((de = (int)(d0 >> Exp_shift)))

1252 z \|= Exp_msk1;

1253 #endif

1254 #ifdef Pack_32

1255 if ((y = d1)) {

1256 if ((k = lo0bits(&y))) {

1257 x[0] = y \| (z << (32 - k));

1258 z >>= k;

1259 }

1260 else

1261 x[0] = y;

1262 #ifndef Sudden_Underflow

1263 i =

1264 #endif

1265 b->wds = (x[1] = z) ? 2 : 1;

1266 }

1267 else {

1268 /* This assertion fails for "1e-500" and other very

1269 * small numbers. It provides the right result (0)

1270 * though. This assert has also been removed from KJS's

1271 * version of dtoa.c.

1272 *

1273 * #ifdef DEBUG

1274 * if (!z) Bug("zero z in b2d");

1275 * #endif

1276 */

1277 k = lo0bits(&z);

1278 x[0] = z;

1279 #ifndef Sudden_Underflow

1280 i =

1281 #endif

1282 b->wds = 1;

1283 k += 32;

1284 }

1285 #else

1286 if (y = d1) {

1287 if (k = lo0bits(&y))

1288 if (k >= 16) {

1289 x[0] = y \| z << 32 - k & 0xffff;

1290 x[1] = z >> k - 16 & 0xffff;

1291 x[2] = z >> k;

1292 i = 2;

1293 }

1294 else {

1295 x[0] = y & 0xffff;

1296 x[1] = y >> 16 \| z << 16 - k & 0xffff;

1297 x[2] = z >> k & 0xffff;

1298 x[3] = z >> k+16;

1299 i = 3;

1300 }

1301 else {

1302 x[0] = y & 0xffff;

1303 x[1] = y >> 16;

1304 x[2] = z & 0xffff;

1305 x[3] = z >> 16;

1306 i = 3;

1307 }

1308 }

1309 else {

1310 #ifdef DEBUG

1311 if (!z)

1312 Bug("Zero passed to d2b");

1313 #endif

1314 k = lo0bits(&z);

1315 if (k >= 16) {

1316 x[0] = z;

1317 i = 0;

1318 }

1319 else {

1320 x[0] = z & 0xffff;

1321 x[1] = z >> 16;

1322 i = 1;

1323 }

1324 k += 32;

1325 }

1326 while(!x[i])

1327 --i;

1328 b->wds = i + 1;

1329 #endif

1330 #ifndef Sudden_Underflow

1331 if (de) {

1332 #endif

1333 #ifdef IBM

1334 *e = (de - Bias - (P-1) << 2) + k;

1335 bits = 4P + 8 - k - hi0bits(word0(d) & Frac_mask);

1336 #else

1337 *e = de - Bias - (P-1) + k;

1338 *bits = P - k;

1339 #endif

1340 #ifndef Sudden_Underflow

1341 }

1342 else {

1343 *e = de - Bias - (P-1) + 1 + k;

1344 #ifdef Pack_32

1345 bits = 32i - hi0bits(x[i-1]);

1346 #else

1347 bits = (i+2)16 - hi0bits(x[i]);

1348 #endif

1349 }

1350 #endif

1351 return b;

1352 }

1353 #undef d0

1354 #undef d1

1355

1356 static double

1357 ratio

1358 #ifdef KR_headers

1359 (a, b) Bigint a, b;

1360 #else

1361 (Bigint a, Bigint b)

1362 #endif

1363 {

1364 U da, db;

1365 int k, ka, kb;

1366

1367 dval(da) = b2d(a, &ka);

1368 dval(db) = b2d(b, &kb);

1369 #ifdef Pack_32

1370 k = ka - kb + 32*(a->wds - b->wds);

1371 #else

1372 k = ka - kb + 16*(a->wds - b->wds);

1373 #endif

1374 #ifdef IBM

1375 if (k > 0) {

1376 word0(da) += (k >> 2)*Exp_msk1;

1377 if (k &= 3)

1378 dval(da) *= 1 << k;

1379 }

1380 else {

1381 k = -k;

1382 word0(db) += (k >> 2)*Exp_msk1;

1383 if (k &= 3)

1384 dval(db) *= 1 << k;

1385 }

1386 #else

1387 if (k > 0)

1388 word0(da) += k*Exp_msk1;

1389 else {

1390 k = -k;

1391 word0(db) += k*Exp_msk1;

1392 }

1393 #endif

1394 return dval(da) / dval(db);

1395 }

1396

1397 static CONST double

1398 tens[] = {

1399 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,

1400 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,

1401 1e20, 1e21, 1e22

1402 #ifdef VAX

1403 , 1e23, 1e24

1404 #endif

1405 };

1406

1407 static CONST double

1408 #ifdef IEEE_Arith

1409 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };

1410 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,

1411 #ifdef Avoid_Underflow

1412 9007199254740992.*9007199254740992.e-256

1413 /* = 2^106 * 1e-53 */

1414 #else

1415 1e-256

1416 #endif

1417 };

1418 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */

1419 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */

1420 #define Scale_Bit 0x10

1421 #define n_bigtens 5

1422 #else

1423 #ifdef IBM

1424 bigtens[] = { 1e16, 1e32, 1e64 };

1425 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };

1426 #define n_bigtens 3

1427 #else

1428 bigtens[] = { 1e16, 1e32 };

1429 static CONST double tinytens[] = { 1e-16, 1e-32 };

1430 #define n_bigtens 2

1431 #endif

1432 #endif

1433

1434 #ifndef IEEE_Arith

1435 #undef INFNAN_CHECK

1436 #endif

1437

1438 #ifdef INFNAN_CHECK

1439

1440 #ifndef NAN_WORD0

1441 #define NAN_WORD0 0x7ff80000

1442 #endif

1443

1444 #ifndef NAN_WORD1

1445 #define NAN_WORD1 0

1446 #endif

1447

1448 static int

1449 match

1450 #ifdef KR_headers

1451 (sp, t) char *sp, t;

1452 #else

1453 (CONST char *sp, char t)

1454 #endif

1455 {

1456 int c, d;

1457 CONST char s = sp;

1458

1459 while(d = *t++) {

1460 if ((c = *++s) >= 'A' && c <= 'Z')

1461 c += 'a' - 'A';

1462 if (c != d)

1463 return 0;

1464 }

1465 *sp = s + 1;

1466 return 1;

1467 }

1468

1469 #ifndef No_Hex_NaN

1470 static void

1471 hexnan

1472 #ifdef KR_headers

1473 (rvp, sp) double rvp; CONST char *sp;

1474 #else

1475 (double rvp, CONST char *sp)

1476 #endif

1477 {

1478 ULong c, x[2];

1479 CONST char *s;

1480 int havedig, udx0, xshift;

1481

1482 x[0] = x[1] = 0;

1483 havedig = xshift = 0;

1484 udx0 = 1;

1485 s = *sp;

1486 while(c = (CONST unsigned char)++s) {

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

1488 c -= '0';

1489 else if (c >= 'a' && c <= 'f')

1490 c += 10 - 'a';

1491 else if (c >= 'A' && c <= 'F')

1492 c += 10 - 'A';

1493 else if (c <= ' ') {

1494 if (udx0 && havedig) {

1495 udx0 = 0;

1496 xshift = 1;

1497 }

1498 continue;

1499 }

1500 else if (/(/ c == ')' && havedig) {

1501 *sp = s + 1;

1502 break;

1503 }

1504 else

1505 return; /* invalid form: don't change sp /

1506 havedig = 1;

1507 if (xshift) {

1508 xshift = 0;

1509 x[0] = x[1];

1510 x[1] = 0;

1511 }

1512 if (udx0)

1513 x[0] = (x[0] << 4) \| (x[1] >> 28);

1514 x[1] = (x[1] << 4) \| c;

1515 }

1516 if ((x[0] &= 0xfffff) \|\| x[1]) {

1517 word0(*rvp) = Exp_mask \| x[0];

1518 word1(*rvp) = x[1];

1519 }

1520 }

1521 #endif /No_Hex_NaN/

1522 #endif /* INFNAN_CHECK */

1523

1524 double

1525 strtod

1526 #ifdef KR_headers

1527 (s00, se) CONST char s00; char *se;

1528 #else

1529 (CONST char s00, char *se)

1530 #endif

1531 {

1532 #ifdef Avoid_Underflow

1533 int scale;

1534 #endif

1535 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,

1536 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;

1537 CONST char s, s0, *s1;

1538 double aadj;

1539 U aadj1, adj, rv, rv0;

1540 Long L;

1541 ULong y, z;

1542 Bigint bb = NULL, bb1, bd = NULL, bd0, bs = NULL, delta = NULL;

1543 #ifdef SET_INEXACT

1544 int inexact, oldinexact;

1545 #endif

1546 #ifdef Honor_FLT_ROUNDS

1547 int rounding;

1548 #endif

1549 #ifdef USE_LOCALE

1550 CONST char *s2;

1551 #endif

1552

1553 sign = nz0 = nz = 0;

1554 dval(rv) = 0.;

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

1556 case '-':

1557 sign = 1;

1558 /* no break */

1559 case '+':

1560 if (*++s)

1561 goto break2;

1562 /* no break */

1563 case 0:

1564 goto ret0;

1565 case '\t':

1566 case '\n':

1567 case '\v':

1568 case '\f':

1569 case '\r':

1570 case ' ':

1571 continue;

1572 default:

1573 goto break2;

1574 }

1575 break2:

1576 if (*s == '0') {

1577 nz0 = 1;

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

1579 if (!*s)

1580 goto ret;

1581 }

1582 s0 = s;

1583 y = z = 0;

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

1585 if (nd < 9)

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

1587 else if (nd < 16)

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

1589 nd0 = nd;

1590 #ifdef USE_LOCALE

1591 s1 = localeconv()->decimal_point;

1592 if (c == *s1) {

1593 c = '.';

1594 if (*++s1) {

1595 s2 = s;

1596 for(;;) {

1597 if (++s2 != s1) {

1598 c = 0;

1599 break;

1600 }

1601 if (!*++s1) {

1602 s = s2;

1603 break;

1604 }

1605 }

1606 }

1607 }

1608 #endif

1609 if (c == '.') {

1610 c = *++s;

1611 if (!nd) {

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

1613 nz++;

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

1615 s0 = s;

1616 nf += nz;

1617 nz = 0;

1618 goto have_dig;

1619 }

1620 goto dig_done;

1621 }

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

1623 have_dig:

1624 nz++;

1625 if (c -= '0') {

1626 nf += nz;

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

1628 if (nd++ < 9)

1629 y *= 10;

1630 else if (nd <= DBL_DIG + 1)

1631 z *= 10;

1632 if (nd++ < 9)

1633 y = 10*y + c;

1634 else if (nd <= DBL_DIG + 1)

1635 z = 10*z + c;

1636 nz = 0;

1637 }

1638 }

1639 }

1640 dig_done:

1641 e = 0;

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

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

1644 goto ret0;

1645 }

1646 s00 = s;

1647 esign = 0;

1648 switch(c = *++s) {

1649 case '-':

1650 esign = 1;

1651 case '+':

1652 c = *++s;

1653 }

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

1655 while(c == '0')

1656 c = *++s;

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

1658 L = c - '0';

1659 s1 = s;

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

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

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

1663 /* Avoid confusion from exponents

1664 * so large that e might overflow.

1665 */

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

1667 else

1668 e = (int)L;

1669 if (esign)

1670 e = -e;

1671 }

1672 else

1673 e = 0;

1674 }

1675 else

1676 s = s00;

1677 }

1678 if (!nd) {

1679 if (!nz && !nz0) {

1680 #ifdef INFNAN_CHECK

1681 /* Check for Nan and Infinity */

1682 switch(c) {

1683 case 'i':

1684 case 'I':

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

1686 --s;

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

1688 ++s;

1689 word0(rv) = 0x7ff00000;

1690 word1(rv) = 0;

1691 goto ret;

1692 }

1693 break;

1694 case 'n':

1695 case 'N':

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

1697 word0(rv) = NAN_WORD0;

1698 word1(rv) = NAN_WORD1;

1699 #ifndef No_Hex_NaN

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

1701 hexnan(&rv, &s);

1702 #endif

1703 goto ret;

1704 }

1705 }

1706 #endif /* INFNAN_CHECK */

1707 ret0:

1708 s = s00;

1709 sign = 0;

1710 }

1711 goto ret;

1712 }

1713 e1 = e -= nf;

1714

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

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

1717 * after is the integer represented by those digits times

1718 * 10*e /

1719

1720 if (!nd0)

1721 nd0 = nd;

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

1723 dval(rv) = y;

1724 if (k > 9) {

1725 #ifdef SET_INEXACT

1726 if (k > DBL_DIG)

1727 oldinexact = get_inexact();

1728 #endif

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

1730 }

1731 bd0 = 0;

1732 if (nd <= DBL_DIG

1733 #ifndef RND_PRODQUOT

1734 #ifndef Honor_FLT_ROUNDS

1735 && Flt_Rounds == 1

1736 #endif

1737 #endif

1738 ) {

1739 if (!e)

1740 goto ret;

1741 if (e > 0) {

1742 if (e <= Ten_pmax) {

1743 #ifdef VAX

1744 goto vax_ovfl_check;

1745 #else

1746 #ifdef Honor_FLT_ROUNDS

1747 /* round correctly FLT_ROUNDS = 2 or 3 */

1748 if (sign) {

1749 rv = -rv;

1750 sign = 0;

1751 }

1752 #endif

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

1754 goto ret;

1755 #endif

1756 }

1757 i = DBL_DIG - nd;

1758 if (e <= Ten_pmax + i) {

1759 /* A fancier test would sometimes let us do

1760 * this for larger i values.

1761 */

1762 #ifdef Honor_FLT_ROUNDS

1763 /* round correctly FLT_ROUNDS = 2 or 3 */

1764 if (sign) {

1765 rv = -rv;

1766 sign = 0;

1767 }

1768 #endif

1769 e -= i;

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

1771 #ifdef VAX

1772 /* VAX exponent range is so narrow we must

1773 * worry about overflow here...

1774 */

1775 vax_ovfl_check:

1776 word0(rv) -= P*Exp_msk1;

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

1778 if ((word0(rv) & Exp_mask)

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

1780 goto ovfl;

1781 word0(rv) += P*Exp_msk1;

1782 #else

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

1784 #endif

1785 goto ret;

1786 }

1787 }

1788 #ifndef Inaccurate_Divide

1789 else if (e >= -Ten_pmax) {

1790 #ifdef Honor_FLT_ROUNDS

1791 /* round correctly FLT_ROUNDS = 2 or 3 */

1792 if (sign) {

1793 rv = -rv;

1794 sign = 0;

1795 }

1796 #endif

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

1798 goto ret;

1799 }

1800 #endif

1801 }

1802 e1 += nd - k;

1803

1804 #ifdef IEEE_Arith

1805 #ifdef SET_INEXACT

1806 inexact = 1;

1807 if (k <= DBL_DIG)

1808 oldinexact = get_inexact();

1809 #endif

1810 #ifdef Avoid_Underflow

1811 scale = 0;

1812 #endif

1813 #ifdef Honor_FLT_ROUNDS

1814 if ((rounding = Flt_Rounds) >= 2) {

1815 if (sign)

1816 rounding = rounding == 2 ? 0 : 2;

1817 else

1818 if (rounding != 2)

1819 rounding = 0;

1820 }

1821 #endif

1822 #endif /IEEE_Arith/

1823

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

1825

1826 if (e1 > 0) {

1827 if ((i = e1 & 15))

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

1829 if (e1 &= ~15) {

1830 if (e1 > DBL_MAX_10_EXP) {

1831 ovfl:

1832 #ifndef NO_ERRNO

1833 errno = ERANGE;

1834 #endif

1835 /* Can't trust HUGE_VAL */

1836 #ifdef IEEE_Arith

1837 #ifdef Honor_FLT_ROUNDS

1838 switch(rounding) {

1839 case 0: /* toward 0 */

1840 case 3: /* toward -infinity */

1841 word0(rv) = Big0;

1842 word1(rv) = Big1;

1843 break;

1844 default:

1845 word0(rv) = Exp_mask;

1846 word1(rv) = 0;

1847 }

1848 #else /Honor_FLT_ROUNDS/

1849 word0(rv) = Exp_mask;

1850 word1(rv) = 0;

1851 #endif /Honor_FLT_ROUNDS/

1852 #ifdef SET_INEXACT

1853 /* set overflow bit */

1854 dval(rv0) = 1e300;

1855 dval(rv0) *= dval(rv0);

1856 #endif

1857 #else /IEEE_Arith/

1858 word0(rv) = Big0;

1859 word1(rv) = Big1;

1860 #endif /IEEE_Arith/

1861 if (bd0)

1862 goto retfree;

1863 goto ret;

1864 }

1865 e1 >>= 4;

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

1867 if (e1 & 1)

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

1869 /* The last multiplication could overflow. */

1870 word0(rv) -= P*Exp_msk1;

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

1872 if ((z = word0(rv) & Exp_mask)

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

1874 goto ovfl;

1875 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {

1876 /* set to largest number */

1877 /* (Can't trust DBL_MAX) */

1878 word0(rv) = Big0;

1879 word1(rv) = Big1;

1880 }

1881 else

1882 word0(rv) += P*Exp_msk1;

1883 }

1884 }

1885 else if (e1 < 0) {

1886 e1 = -e1;

1887 if ((i = e1 & 15))

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

1889 if (e1 >>= 4) {

1890 if (e1 >= 1 << n_bigtens)

1891 goto undfl;

1892 #ifdef Avoid_Underflow

1893 if (e1 & Scale_Bit)

1894 scale = 2*P;

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

1896 if (e1 & 1)

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

1898 if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)

1899 >> Exp_shift)) > 0) {

1900 /* scaled rv is denormal; zap j low bits */

1901 if (j >= 32) {

1902 word1(rv) = 0;

1903 if (j >= 53)

1904 word0(rv) = (P+2)*Exp_msk1;

1905 else

1906 word0(rv) &= 0xffffffff << (j-32);

1907 }

1908 else

1909 word1(rv) &= 0xffffffff << j;

1910 }

1911 #else

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

1913 if (e1 & 1)

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

1915 /* The last multiplication could underflow. */

1916 dval(rv0) = dval(rv);

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

1918 if (!dval(rv)) {

1919 dval(rv) = 2.*dval(rv0);

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

1921 #endif

1922 if (!dval(rv)) {

1923 undfl:

1924 dval(rv) = 0.;

1925 #ifndef NO_ERRNO

1926 errno = ERANGE;

1927 #endif

1928 if (bd0)

1929 goto retfree;

1930 goto ret;

1931 }

1932 #ifndef Avoid_Underflow

1933 word0(rv) = Tiny0;

1934 word1(rv) = Tiny1;

1935 /* The refinement below will clean

1936 * this approximation up.

1937 */

1938 }

1939 #endif

1940 }

1941 }

1942

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

1944

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

1946

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

1948

1949 for(;;) {

1950 bd = Balloc(bd0->k);

1951 Bcopy(bd, bd0);

1952 bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */

1953 bs = i2b(1);

1954

1955 if (e >= 0) {

1956 bb2 = bb5 = 0;

1957 bd2 = bd5 = e;

1958 }

1959 else {

1960 bb2 = bb5 = -e;

1961 bd2 = bd5 = 0;

1962 }

1963 if (bbe >= 0)

1964 bb2 += bbe;

1965 else

1966 bd2 -= bbe;

1967 bs2 = bb2;

1968 #ifdef Honor_FLT_ROUNDS

1969 if (rounding != 1)

1970 bs2++;

1971 #endif

1972 #ifdef Avoid_Underflow

1973 j = bbe - scale;

1974 i = j + bbbits - 1; /* logb(rv) */

1975 if (i < Emin) /* denormal */

1976 j += P - Emin;

1977 else

1978 j = P + 1 - bbbits;

1979 #else /Avoid_Underflow/

1980 #ifdef Sudden_Underflow

1981 #ifdef IBM

1982 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);

1983 #else

1984 j = P + 1 - bbbits;

1985 #endif

1986 #else /Sudden_Underflow/

1987 j = bbe;

1988 i = j + bbbits - 1; /* logb(rv) */

1989 if (i < Emin) /* denormal */

1990 j += P - Emin;

1991 else

1992 j = P + 1 - bbbits;

1993 #endif /Sudden_Underflow/

1994 #endif /Avoid_Underflow/

1995 bb2 += j;

1996 bd2 += j;

1997 #ifdef Avoid_Underflow

1998 bd2 += scale;

1999 #endif

2000 i = bb2 < bd2 ? bb2 : bd2;

2001 if (i > bs2)

2002 i = bs2;

2003 if (i > 0) {

2004 bb2 -= i;

2005 bd2 -= i;

2006 bs2 -= i;

2007 }

2008 if (bb5 > 0) {

2009 bs = pow5mult(bs, bb5);

2010 bb1 = mult(bs, bb);

2011 Bfree(bb);

2012 bb = bb1;

2013 }

2014 if (bb2 > 0)

2015 bb = lshift(bb, bb2);

2016 if (bd5 > 0)

2017 bd = pow5mult(bd, bd5);

2018 if (bd2 > 0)

2019 bd = lshift(bd, bd2);

2020 if (bs2 > 0)

2021 bs = lshift(bs, bs2);

2022 delta = diff(bb, bd);

2023 dsign = delta->sign;

2024 delta->sign = 0;

2025 i = cmp(delta, bs);

2026 #ifdef Honor_FLT_ROUNDS

2027 if (rounding != 1) {

2028 if (i < 0) {

2029 /* Error is less than an ulp */

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

2031 /* exact */

2032 #ifdef SET_INEXACT

2033 inexact = 0;

2034 #endif

2035 break;

2036 }

2037 if (rounding) {

2038 if (dsign) {

2039 dval(adj) = 1.;

2040 goto apply_adj;

2041 }

2042 }

2043 else if (!dsign) {

2044 dval(adj) = -1.;

2045 if (!word1(rv)

2046 && !(word0(rv) & Frac_mask)) {

2047 y = word0(rv) & Exp_mask;

2048 #ifdef Avoid_Underflow

2049 if (!scale \|\| y > 2PExp_msk1)

2050 #else

2051 if (y)

2052 #endif

2053 {

2054 delta = lshift(delta,Log2P);

2055 if (cmp(delta, bs) <= 0)

2056 dval(adj) = -0.5;

2057 }

2058 }

2059 apply_adj:

2060 #ifdef Avoid_Underflow

2061 if (scale && (y = word0(rv) & Exp_mask)

2062 <= 2PExp_msk1)

2063 word0(adj) += (2P+1)Exp_msk1 - y;

2064 #else

2065 #ifdef Sudden_Underflow

2066 if ((word0(rv) & Exp_mask) <=

2067 P*Exp_msk1) {

2068 word0(rv) += P*Exp_msk1;

2069 dval(rv) += dval(adj)*ulp(dval(r v));

2070 word0(rv) -= P*Exp_msk1;

2071 }

2072 else

2073 #endif /Sudden_Underflow/

2074 #endif /Avoid_Underflow/

2075 dval(rv) += dval(adj)*ulp(dval(rv));

2076 }

2077 break;

2078 }

2079 dval(adj) = ratio(delta, bs);

2080 if (dval(adj) < 1.)

2081 dval(adj) = 1.;

2082 if (dval(adj) <= 0x7ffffffe) {

2083 /* adj = rounding ? ceil(adj) : floor(adj); */

2084 y = dval(adj);

2085 if (y != dval(adj)) {

2086 if (!((rounding>>1) ^ dsign))

2087 y++;

2088 dval(adj) = y;

2089 }

2090 }

2091 #ifdef Avoid_Underflow

2092 if (scale && (y = word0(rv) & Exp_mask) <= 2PExp_msk1)

2093 word0(adj) += (2P+1)Exp_msk1 - y;

2094 #else

2095 #ifdef Sudden_Underflow

2096 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {

2097 word0(rv) += P*Exp_msk1;

2098 dval(adj) *= ulp(dval(rv));

2099 if (dsign)

2100 dval(rv) += dval(adj);

2101 else

2102 dval(rv) -= dval(adj);

2103 word0(rv) -= P*Exp_msk1;

2104 goto cont;

2105 }

2106 #endif /Sudden_Underflow/

2107 #endif /Avoid_Underflow/

2108 dval(adj) *= ulp(dval(rv));

2109 if (dsign)

2110 dval(rv) += dval(adj);

2111 else

2112 dval(rv) -= dval(adj);

2113 goto cont;

2114 }

2115 #endif /Honor_FLT_ROUNDS/

2116

2117 if (i < 0) {

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

2119 * special case of mantissa a power of two.

2120 */

2121 if (dsign \|\| word1(rv) \|\| word0(rv) & Bndry_mask

2122 #ifdef IEEE_Arith

2123 #ifdef Avoid_Underflow

2124 \|\| (word0(rv) & Exp_mask) <= (2P+1)Exp_msk1

2125 #else

2126 \|\| (word0(rv) & Exp_mask) <= Exp_msk1

2127 #endif

2128 #endif

2129 ) {

2130 #ifdef SET_INEXACT

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

2132 inexact = 0;

2133 #endif

2134 break;

2135 }

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

2137 /* exact result */

2138 #ifdef SET_INEXACT

2139 inexact = 0;

2140 #endif

2141 break;

2142 }

2143 delta = lshift(delta,Log2P);

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

2145 goto drop_down;

2146 break;

2147 }

2148 if (i == 0) {

2149 /* exactly half-way between */

2150 if (dsign) {

2151 if ((word0(rv) & Bndry_mask1) == Bndry_mask1

2152 && word1(rv) == (

2153 #ifdef Avoid_Underflow

2154 (scale && (y = word0(rv) & Exp_mask) <= 2PExp_msk1)

2155 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :

2156 #endif

2157 0xffffffff)) {

2158 /boundary case -- increment exponent/

2159 word0(rv) = (word0(rv) & Exp_mask)

2160 + Exp_msk1

2161 #ifdef IBM

2162 \| Exp_msk1 >> 4

2163 #endif

2164 ;

2165 word1(rv) = 0;

2166 #ifdef Avoid_Underflow

2167 dsign = 0;

2168 #endif

2169 break;

2170 }

2171 }

2172 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {

2173 drop_down:

2174 /* boundary case -- decrement exponent */

2175 #ifdef Sudden_Underflow /{{/

2176 L = word0(rv) & Exp_mask;

2177 #ifdef IBM

2178 if (L < Exp_msk1)

2179 #else

2180 #ifdef Avoid_Underflow

2181 if (L <= (scale ? (2P+1)Exp_msk1 : Exp_msk1))

2182 #else

2183 if (L <= Exp_msk1)

2184 #endif /Avoid_Underflow/

2185 #endif /IBM/

2186 goto undfl;

2187 L -= Exp_msk1;

2188 #else /Sudden_Underflow}{/

2189 #ifdef Avoid_Underflow

2190 if (scale) {

2191 L = word0(rv) & Exp_mask;

2192 if (L <= (2P+1)Exp_msk1) {

2193 if (L > (P+2)*Exp_msk1)

2194 /* round even ==> */

2195 /* accept rv */

2196 break;

2197 /* rv = smallest denormal */

2198 goto undfl;

2199 }

2200 }

2201 #endif /Avoid_Underflow/

2202 L = (word0(rv) & Exp_mask) - Exp_msk1;

2203 #endif /Sudden_Underflow}}/

2204 word0(rv) = L \| Bndry_mask1;

2205 word1(rv) = 0xffffffff;

2206 #ifdef IBM

2207 goto cont;

2208 #else

2209 break;

2210 #endif

2211 }

2212 #ifndef ROUND_BIASED

2213 if (!(word1(rv) & LSB))

2214 break;

2215 #endif

2216 if (dsign)

2217 dval(rv) += ulp(dval(rv));

2218 #ifndef ROUND_BIASED

2219 else {

2220 dval(rv) -= ulp(dval(rv));

2221 #ifndef Sudden_Underflow

2222 if (!dval(rv))

2223 goto undfl;

2224 #endif

2225 }

2226 #ifdef Avoid_Underflow

2227 dsign = 1 - dsign;

2228 #endif

2229 #endif

2230 break;

2231 }

2232 if ((aadj = ratio(delta, bs)) <= 2.) {

2233 if (dsign)

2234 aadj = dval(aadj1) = 1.;

2235 else if (word1(rv) \|\| word0(rv) & Bndry_mask) {

2236 #ifndef Sudden_Underflow

2237 if (word1(rv) == Tiny1 && !word0(rv))

2238 goto undfl;

2239 #endif

2240 aadj = 1.;

2241 dval(aadj1) = -1.;

2242 }

2243 else {

2244 /* special case -- power of FLT_RADIX to be */

2245 /* rounded down... */

2246

2247 if (aadj < 2./FLT_RADIX)

2248 aadj = 1./FLT_RADIX;

2249 else

2250 aadj *= 0.5;

2251 dval(aadj1) = -aadj;

2252 }

2253 }

2254 else {

2255 aadj *= 0.5;

2256 dval(aadj1) = dsign ? aadj : -aadj;

2257 #ifdef Check_FLT_ROUNDS

2258 switch(Rounding) {

2259 case 2: /* towards +infinity */

2260 dval(aadj1) -= 0.5;

2261 break;

2262 case 0: /* towards 0 */

2263 case 3: /* towards -infinity */

2264 dval(aadj1) += 0.5;

2265 }

2266 #else

2267 if (Flt_Rounds == 0)

2268 dval(aadj1) += 0.5;

2269 #endif /Check_FLT_ROUNDS/

2270 }

2271 y = word0(rv) & Exp_mask;

2272

2273 /* Check for overflow */

2274

2275 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {

2276 dval(rv0) = dval(rv);

2277 word0(rv) -= P*Exp_msk1;

2278 dval(adj) = dval(aadj1) * ulp(dval(rv));

2279 dval(rv) += dval(adj);

2280 if ((word0(rv) & Exp_mask) >=

2281 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {

2282 if (word0(rv0) == Big0 && word1(rv0) == Big1)

2283 goto ovfl;

2284 word0(rv) = Big0;

2285 word1(rv) = Big1;

2286 goto cont;

2287 }

2288 else

2289 word0(rv) += P*Exp_msk1;

2290 }

2291 else {

2292 #ifdef Avoid_Underflow

2293 if (scale && y <= 2PExp_msk1) {

2294 if (aadj <= 0x7fffffff) {

2295 if ((z = aadj) <= 0)

2296 z = 1;

2297 aadj = z;

2298 dval(aadj1) = dsign ? aadj : -aadj;

2299 }

2300 word0(aadj1) += (2P+1)Exp_msk1 - y;

2301 }

2302 dval(adj) = dval(aadj1) * ulp(dval(rv));

2303 dval(rv) += dval(adj);

2304 #else

2305 #ifdef Sudden_Underflow

2306 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {

2307 dval(rv0) = dval(rv);

2308 word0(rv) += P*Exp_msk1;

2309 dval(adj) = dval(aadj1) * ulp(dval(rv));

2310 dval(rv) += dval(adj);

2311 #ifdef IBM

2312 if ((word0(rv) & Exp_mask) < P*Exp_msk1)

2313 #else

2314 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)

2315 #endif

2316 {

2317 if (word0(rv0) == Tiny0

2318 && word1(rv0) == Tiny1)

2319 goto undfl;

2320 word0(rv) = Tiny0;

2321 word1(rv) = Tiny1;

2322 goto cont;

2323 }

2324 else

2325 word0(rv) -= P*Exp_msk1;

2326 }

2327 else {

2328 dval(adj) = dval(aadj1) * ulp(dval(rv));

2329 dval(rv) += dval(adj);

2330 }

2331 #else /Sudden_Underflow/

2332 /* Compute adj so that the IEEE rounding rules will

2333 * correctly round rv + adj in some half-way cases.

2334 * If rv * ulp(rv) is denormalized (i.e.,

2335 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid

2336 * trouble from bits lost to denormalization;

2337 * example: 1.2e-307 .

2338 */

2339 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {

2340 dval(aadj1) = (double)(int)(aadj + 0.5);

2341 if (!dsign)

2342 dval(aadj1) = -dval(aadj1);

2343 }

2344 dval(adj) = dval(aadj1) * ulp(dval(rv));

2345 dval(rv) += dval(adj);

2346 #endif /Sudden_Underflow/

2347 #endif /Avoid_Underflow/

2348 }

2349 z = word0(rv) & Exp_mask;

2350 #ifndef SET_INEXACT

2351 #ifdef Avoid_Underflow

2352 if (!scale)

2353 #endif

2354 if (y == z) {

2355 /* Can we stop now? */

2356 L = (Long)aadj;

2357 aadj -= L;

2358 /* The tolerances below are conservative. */

2359 if (dsign \|\| word1(rv) \|\| word0(rv) & Bndry_mask) {

2360 if (aadj < .4999999 \|\| aadj > .5000001)

2361 break;

2362 }

2363 else if (aadj < .4999999/FLT_RADIX)

2364 break;

2365 }

2366 #endif

2367 cont:

2368 Bfree(bb);

2369 Bfree(bd);

2370 Bfree(bs);

2371 Bfree(delta);

2372 }

2373 #ifdef SET_INEXACT

2374 if (inexact) {

2375 if (!oldinexact) {

2376 word0(rv0) = Exp_1 + (70 << Exp_shift);

2377 word1(rv0) = 0;

2378 dval(rv0) += 1.;

2379 }

2380 }

2381 else if (!oldinexact)

2382 clear_inexact();

2383 #endif

2384 #ifdef Avoid_Underflow

2385 if (scale) {

2386 word0(rv0) = Exp_1 - 2PExp_msk1;

2387 word1(rv0) = 0;

2388 dval(rv) *= dval(rv0);

2389 #ifndef NO_ERRNO

2390 /* try to avoid the bug of testing an 8087 register value */

2391 if (word0(rv) == 0 && word1(rv) == 0)

2392 errno = ERANGE;

2393 #endif

2394 }

2395 #endif /* Avoid_Underflow */

2396 #ifdef SET_INEXACT

2397 if (inexact && !(word0(rv) & Exp_mask)) {

2398 /* set underflow bit */

2399 dval(rv0) = 1e-300;

2400 dval(rv0) *= dval(rv0);

2401 }

2402 #endif

2403 retfree:

2404 Bfree(bb);

2405 Bfree(bd);

2406 Bfree(bs);

2407 Bfree(bd0);

2408 Bfree(delta);

2409 ret:

2410 if (se)

2411 se = (char )s;

2412 return sign ? -dval(rv) : dval(rv);

2413 }

2414

2415 static int

2416 quorem

2417 #ifdef KR_headers

2418 (b, S) Bigint b, S;

2419 #else

2420 (Bigint b, Bigint S)

2421 #endif

2422 {

2423 int n;

2424 ULong bx, bxe, q, sx, sxe;

2425 #ifdef ULLong

2426 ULLong borrow, carry, y, ys;

2427 #else

2428 ULong borrow, carry, y, ys;

2429 #ifdef Pack_32

2430 ULong si, z, zs;

2431 #endif

2432 #endif

2433

2434 n = S->wds;

2435 #ifdef DEBUG

2436 /debug/ if (b->wds > n)

2437 /debug/ Bug("oversize b in quorem");

2438 #endif

2439 if (b->wds < n)

2440 return 0;

2441 sx = S->x;

2442 sxe = sx + --n;

2443 bx = b->x;

2444 bxe = bx + n;

2445 q = bxe / (sxe + 1); /* ensure q <= true quotient */

2446 #ifdef DEBUG

2447 /debug/ if (q > 9)

2448 /debug/ Bug("oversized quotient in quorem");

2449 #endif

2450 if (q) {

2451 borrow = 0;

2452 carry = 0;

2453 do {

2454 #ifdef ULLong

2455 ys = sx++ (ULLong)q + carry;

2456 carry = ys >> 32;

2457 y = *bx - (ys & FFFFFFFF) - borrow;

2458 borrow = y >> 32 & (ULong)1;

2459 *bx++ = y & FFFFFFFF;

2460 #else

2461 #ifdef Pack_32

2462 si = *sx++;

2463 ys = (si & 0xffff) * q + carry;

2464 zs = (si >> 16) * q + (ys >> 16);

2465 carry = zs >> 16;

2466 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;

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

2468 z = (*bx >> 16) - (zs & 0xffff) - borrow;

2469 borrow = (z & 0x10000) >> 16;

2470 Storeinc(bx, z, y);

2471 #else

2472 ys = sx++ q + carry;

2473 carry = ys >> 16;

2474 y = *bx - (ys & 0xffff) - borrow;

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

2476 *bx++ = y & 0xffff;

2477 #endif

2478 #endif

2479 }

2480 while(sx <= sxe);

2481 if (!*bxe) {

2482 bx = b->x;

2483 while(--bxe > bx && !*bxe)

2484 --n;

2485 b->wds = n;

2486 }

2487 }

2488 if (cmp(b, S) >= 0) {

2489 q++;

2490 borrow = 0;

2491 carry = 0;

2492 bx = b->x;

2493 sx = S->x;

2494 do {

2495 #ifdef ULLong

2496 ys = *sx++ + carry;

2497 carry = ys >> 32;

2498 y = *bx - (ys & FFFFFFFF) - borrow;

2499 borrow = y >> 32 & (ULong)1;

2500 *bx++ = y & FFFFFFFF;

2501 #else

2502 #ifdef Pack_32

2503 si = *sx++;

2504 ys = (si & 0xffff) + carry;

2505 zs = (si >> 16) + (ys >> 16);

2506 carry = zs >> 16;

2507 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;

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

2509 z = (*bx >> 16) - (zs & 0xffff) - borrow;

2510 borrow = (z & 0x10000) >> 16;

2511 Storeinc(bx, z, y);

2512 #else

2513 ys = *sx++ + carry;

2514 carry = ys >> 16;

2515 y = *bx - (ys & 0xffff) - borrow;

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

2517 *bx++ = y & 0xffff;

2518 #endif

2519 #endif

2520 }

2521 while(sx <= sxe);

2522 bx = b->x;

2523 bxe = bx + n;

2524 if (!*bxe) {

2525 while(--bxe > bx && !*bxe)

2526 --n;

2527 b->wds = n;

2528 }

2529 }

2530 return q;

2531 }

2532

2533 #ifndef MULTIPLE_THREADS

2534 static char *dtoa_result;

2535 #endif

2536

2537 static char *

2538 #ifdef KR_headers

2539 rv_alloc(i) int i;

2540 #else

2541 rv_alloc(int i)

2542 #endif

2543 {

2544 int j, k, *r;

2545

2546 j = sizeof(ULong);

2547 for(k = 0;

2548 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i;

2549 j <<= 1)

2550 k++;

2551 r = (int*)Balloc(k);

2552 *r = k;

2553 return

2554 #ifndef MULTIPLE_THREADS

2555 dtoa_result =

2556 #endif

2557 (char *)(r+1);

2558 }

2559

2560 static char *

2561 #ifdef KR_headers

2562 nrv_alloc(s, rve, n) char s, *rve; int n;

2563 #else

2564 nrv_alloc(const char s, char *rve, int n)

2565 #endif

2566 {

2567 char rv, t;

2568

2569 t = rv = rv_alloc(n);

2570 while ((t = s++)) t++;

2571 if (rve)

2572 *rve = t;

2573 return rv;

2574 }

2575

2576 /* freedtoa(s) must be used to free values s returned by dtoa

2577 * when MULTIPLE_THREADS is #defined. It should be used in all cases,

2578 * but for consistency with earlier versions of dtoa, it is optional

2579 * when MULTIPLE_THREADS is not defined.

2580 */

2581

2582 void

2583 #ifdef KR_headers

2584 freedtoa(s) char *s;

2585 #else

2586 freedtoa(char *s)

2587 #endif

2588 {

2589 Bigint b = (Bigint )((int *)s - 1);

2590 b->maxwds = 1 << (b->k = (int)b);

2591 Bfree(b);

2592 #ifndef MULTIPLE_THREADS

2593 if (s == dtoa_result)

2594 dtoa_result = 0;

2595 #endif

2596 }

2597

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

2599 *

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

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

2602 *

2603 * Modifications:

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

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

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

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

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

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

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

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

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

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

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

2615 * as would satisfaction of the stopping test with strict

2616 * inequality.

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

2618 * quantities.

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

2620 * we use floating-point arithmetic rather than resorting

2621 * to multiple-precision integers.

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

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

2624 * multiple-precision integer arithmetic only if we cannot

2625 * guarantee that the floating-point calculation has given

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

2627 * "uniformly" distributed input, the probability is

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

2629 * calculation.

2630 */

2631

2632 char *

2633 dtoa

2634 #ifdef KR_headers

2635 (dd, mode, ndigits, decpt, sign, rve)

2636 double dd; int mode, ndigits, decpt, sign; char **rve;

2637 #else

2638 (double dd, int mode, int ndigits, int decpt, int sign, char **rve)

2639 #endif

2640 {

2641 /* Arguments ndigits, decpt, sign are similar to those

2642 of ecvt and fcvt; trailing zeros are suppressed from

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

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

2645 then *decpt is set to 9999.

2646

2647 mode:

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

2649 and rounded to nearest.

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

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

2652 1e23 whereas mode 1 gives 9.999999999999999e22.

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

2654 return value similar to that of ecvt, except

2655 that trailing zeros are suppressed.

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

2657 gives a return value similar to that from fcvt,

2658 except that trailing zeros are suppressed, and

2659 ndigits can be negative.

2660 4,5 ==> similar to 2 and 3, respectively, but (in

2661 round-nearest mode) with the tests of mode 0 to

2662 possibly return a shorter string that rounds to d.

2663 With IEEE arithmetic and compilation with

2664 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same

2665 as modes 2 and 3 when FLT_ROUNDS != 1.

2666 6-9 ==> Debugging modes similar to mode - 4: don't try

2667 fast floating-point estimate (if applicable).

2668

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

2670

2671 Sufficient space is allocated to the return value

2672 to hold the suppressed trailing zeros.

2673 */

2674

2675 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,

2676 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,

2677 spec_case, try_quick, bias_round_up;

2678 Long L;

2679 #ifndef Sudden_Underflow

2680 int denorm;

2681 ULong x;

2682 #endif

2683 Bigint b, b1, delta, mlo, mhi, S;

2684 double ds;

2685 U d2, eps;

2686 char s, s0;

2687 #ifdef Honor_FLT_ROUNDS

2688 int rounding;

2689 #endif

2690 #ifdef SET_INEXACT

2691 int inexact, oldinexact;

2692 #endif

2693 U d;

2694 dval(d) = dd;

2695

2696 /* In mode 2 and 3 we bias rounding up when there are ties. */

2697 bias_round_up = mode == 2 \|\| mode == 3;

2698

2699 ilim = ilim1 = 0; /* to avoid Google3 compiler warnings */

2700

2701 #ifndef MULTIPLE_THREADS

2702 if (dtoa_result) {

2703 freedtoa(dtoa_result);

2704 dtoa_result = 0;

2705 }

2706 #endif

2707

2708 if (word0(d) & Sign_bit) {

2709 /* set sign for everything, including 0's and NaNs */

2710 *sign = 1;

2711 word0(d) &= ~Sign_bit; /* clear sign bit */

2712 }

2713 else

2714 *sign = 0;

2715

2716 #if defined(IEEE_Arith) + defined(VAX)

2717 #ifdef IEEE_Arith

2718 if ((word0(d) & Exp_mask) == Exp_mask)

2719 #else

2720 if (word0(d) == 0x8000)

2721 #endif

2722 {

2723 /* Infinity or NaN */

2724 *decpt = 9999;

2725 #ifdef IEEE_Arith

2726 if (!word1(d) && !(word0(d) & 0xfffff))

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

2728 #endif

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

2730 }

2731 #endif

2732 #ifdef IBM

2733 dval(d) += 0; /* normalize */

2734 #endif

2735 if (!dval(d)) {

2736 *decpt = 1;

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

2738 }

2739

2740 #ifdef SET_INEXACT

2741 try_quick = oldinexact = get_inexact();

2742 inexact = 1;

2743 #endif

2744 #ifdef Honor_FLT_ROUNDS

2745 if ((rounding = Flt_Rounds) >= 2) {

2746 if (*sign)

2747 rounding = rounding == 2 ? 0 : 2;

2748 else

2749 if (rounding != 2)

2750 rounding = 0;

2751 }

2752 #endif

2753

2754 b = d2b(dval(d), &be, &bbits);

2755 #ifdef Sudden_Underflow

2756 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));

2757 #else

2758 if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {

2759 #endif

2760 dval(d2) = dval(d);

2761 word0(d2) &= Frac_mask1;

2762 word0(d2) \|= Exp_11;

2763 #ifdef IBM

2764 if (j = 11 - hi0bits(word0(d2) & Frac_mask))

2765 dval(d2) /= 1 << j;

2766 #endif

2767

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

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

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

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

2772 *

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

2774 *

2775 * k = (i - Bias)*0.301029995663981

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

2777 *

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

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

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

2781 * to compensate for any error in the multiplication of

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

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

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

2785 * Hence we adjust the constant term to 0.1760912590558.

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

2787 * but this is probably not worthwhile.)

2788 */

2789

2790 i -= Bias;

2791 #ifdef IBM

2792 i <<= 2;

2793 i += j;

2794 #endif

2795 #ifndef Sudden_Underflow

2796 denorm = 0;

2797 }

2798 else {

2799 /* d is denormalized */

2800

2801 i = bbits + be + (Bias + (P-1) - 1);

2802 x = i > 32 ? (word0(d) << (64 - i)) \| (word1(d) >> (i - 32))

2803 : word1(d) << (32 - i);

2804 dval(d2) = x;

2805 word0(d2) -= 31Exp_msk1; / adjust exponent */

2806 i -= (Bias + (P-1) - 1) + 1;

2807 denorm = 1;

2808 }

2809 #endif

2810 ds = (dval(d2)-1.5)0.289529654602168 + 0.1760912590558 + i0.3010299956 63981;

2811 k = (int)ds;

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

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

2814 k_check = 1;

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

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

2817 k--;

2818 k_check = 0;

2819 }

2820 j = bbits - i - 1;

2821 if (j >= 0) {

2822 b2 = 0;

2823 s2 = j;

2824 }

2825 else {

2826 b2 = -j;

2827 s2 = 0;

2828 }

2829 if (k >= 0) {

2830 b5 = 0;

2831 s5 = k;

2832 s2 += k;

2833 }

2834 else {

2835 b2 -= k;

2836 b5 = -k;

2837 s5 = 0;

2838 }

2839 if (mode < 0 \|\| mode > 9)

2840 mode = 0;

2841

2842 #ifndef SET_INEXACT

2843 #ifdef Check_FLT_ROUNDS

2844 try_quick = Rounding == 1;

2845 #else

2846 try_quick = 1;

2847 #endif

2848 #endif /SET_INEXACT/

2849

2850 if (mode > 5) {

2851 mode -= 4;

2852 try_quick = 0;

2853 }

2854 leftright = 1;

2855 switch(mode) {

2856 case 0:

2857 case 1:

2858 ilim = ilim1 = -1;

2859 i = 18;

2860 ndigits = 0;

2861 break;

2862 case 2:

2863 leftright = 0;

2864 /* no break */

2865 case 4:

2866 if (ndigits <= 0)

2867 ndigits = 1;

2868 ilim = ilim1 = i = ndigits;

2869 break;

2870 case 3:

2871 leftright = 0;

2872 /* no break */

2873 case 5:

2874 i = ndigits + k + 1;

2875 ilim = i;

2876 ilim1 = i - 1;

2877 if (i <= 0)

2878 i = 1;

2879 }

2880 s = s0 = rv_alloc(i);

2881

2882 #ifdef Honor_FLT_ROUNDS

2883 if (mode > 1 && rounding != 1)

2884 leftright = 0;

2885 #endif

2886

2887 if (ilim >= 0 && ilim <= Quick_max && try_quick) {

2888

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

2890

2891 i = 0;

2892 dval(d2) = dval(d);

2893 k0 = k;

2894 ilim0 = ilim;

2895 ieps = 2; /* conservative */

2896 if (k > 0) {

2897 ds = tens[k&0xf];

2898 j = k >> 4;

2899 if (j & Bletch) {

2900 /* prevent overflows */

2901 j &= Bletch - 1;

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

2903 ieps++;

2904 }

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

2906 if (j & 1) {

2907 ieps++;

2908 ds *= bigtens[i];

2909 }

2910 dval(d) /= ds;

2911 }

2912 else if ((j1 = -k)) {

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

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

2915 if (j & 1) {

2916 ieps++;

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

2918 }

2919 }

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

2921 if (ilim1 <= 0)

2922 goto fast_failed;

2923 ilim = ilim1;

2924 k--;

2925 dval(d) *= 10.;

2926 ieps++;

2927 }

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

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

2930 if (ilim == 0) {

2931 S = mhi = 0;

2932 dval(d) -= 5.;

2933 if (dval(d) > dval(eps))

2934 goto one_digit;

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

2936 goto no_digits;

2937 goto fast_failed;

2938 }

2939 #ifndef No_leftright

2940 if (leftright) {

2941 /* Use Steele & White method of only

2942 * generating digits needed.

2943 */

2944 dval(eps) = 0.5/tens[ilim-1] - dval(eps);

2945 for(i = 0;;) {

2946 L = dval(d);

2947 dval(d) -= L;

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

2949 if (dval(d) < dval(eps))

2950 goto ret1;

2951 if (1. - dval(d) < dval(eps))

2952 goto bump_up;

2953 if (++i >= ilim)

2954 break;

2955 dval(eps) *= 10.;

2956 dval(d) *= 10.;

2957 }

2958 }

2959 else {

2960 #endif

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

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

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

2964 L = (Long)(dval(d));

2965 if (!(dval(d) -= L))

2966 ilim = i;

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

2968 if (i == ilim) {

2969 if (dval(d) > 0.5 + dval(eps))

2970 goto bump_up;

2971 else if (dval(d) < 0.5 - dval(eps)) {

2972 while(*--s == '0');

2973 s++;

2974 goto ret1;

2975 }

2976 break;

2977 }

2978 }

2979 #ifndef No_leftright

2980 }

2981 #endif

2982 fast_failed:

2983 s = s0;

2984 dval(d) = dval(d2);

2985 k = k0;

2986 ilim = ilim0;

2987 }

2988

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

2990

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

2992 /* Yes. */

2993 ds = tens[k];

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

2995 S = mhi = 0;

2996 if (ilim < 0 \|\| dval(d) < 5ds \|\| ((dval(d) == 5ds) && !bias_round_up))

2997 goto no_digits;

2998 goto one_digit;

2999 }

3000

3001 /* Limit looping by the number of digits to produce.

3002 * Firefox had a crash bug because some plugins reduce

3003 * the precision of double arithmetic. With reduced

3004 * precision "dval(d) -= L*ds" might be imprecise and

3005 * d might not become zero and the loop might not

3006 * terminate.

3007 *

3008 * See https://bugzilla.mozilla.org/show_bug.cgi?id=358569

3009 */

3010 for(i = 1; i <= k+1; i++, dval(d) *= 10.) {

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

3012 dval(d) -= L*ds;

3013 #ifdef Check_FLT_ROUNDS

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

3015 if (dval(d) < 0) {

3016 L--;

3017 dval(d) += ds;

3018 }

3019 #endif

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

3021 if (!dval(d)) {

3022 #ifdef SET_INEXACT

3023 inexact = 0;

3024 #endif

3025 break;

3026 }

3027 if (i == ilim) {

3028 #ifdef Honor_FLT_ROUNDS

3029 if (mode > 1)

3030 switch(rounding) {

3031 case 0: goto ret1;

3032 case 2: goto bump_up;

3033 }

3034 #endif

3035 dval(d) += dval(d);

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

3037 bump_up:

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

3039 if (s == s0) {

3040 k++;

3041 *s = '0';

3042 break;

3043 }

3044 ++*s++;

3045 }

3046 break;

3047 }

3048 }

3049 goto ret1;

3050 }

3051

3052 m2 = b2;

3053 m5 = b5;

3054 mhi = mlo = 0;

3055 if (leftright) {

3056 i =

3057 #ifndef Sudden_Underflow

3058 denorm ? be + (Bias + (P-1) - 1 + 1) :

3059 #endif

3060 #ifdef IBM

3061 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);

3062 #else

3063 1 + P - bbits;

3064 #endif

3065 b2 += i;

3066 s2 += i;

3067 mhi = i2b(1);

3068 }

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

3070 i = m2 < s2 ? m2 : s2;

3071 b2 -= i;

3072 m2 -= i;

3073 s2 -= i;

3074 }

3075 if (b5 > 0) {

3076 if (leftright) {

3077 if (m5 > 0) {

3078 mhi = pow5mult(mhi, m5);

3079 b1 = mult(mhi, b);

3080 Bfree(b);

3081 b = b1;

3082 }

3083 if ((j = b5 - m5))

3084 b = pow5mult(b, j);

3085 }

3086 else

3087 b = pow5mult(b, b5);

3088 }

3089 S = i2b(1);

3090 if (s5 > 0)

3091 S = pow5mult(S, s5);

3092

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

3094

3095 spec_case = 0;

3096 if ((mode < 2 \|\| leftright)

3097 #ifdef Honor_FLT_ROUNDS

3098 && rounding == 1

3099 #endif

3100 ) {

3101 if (!word1(d) && !(word0(d) & Bndry_mask)

3102 #ifndef Sudden_Underflow

3103 && word0(d) & (Exp_mask & ~Exp_msk1)

3104 #endif

3105 ) {

3106 /* The special case */

3107 b2 += Log2P;

3108 s2 += Log2P;

3109 spec_case = 1;

3110 }

3111 }

3112

3113 /* Arrange for convenient computation of quotients:

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

3115 *

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

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

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

3119 */

3120 #ifdef Pack_32

3121 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f))

3122 i = 32 - i;

3123 #else

3124 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf))

3125 i = 16 - i;

3126 #endif

3127 if (i > 4) {

3128 i -= 4;

3129 b2 += i;

3130 m2 += i;

3131 s2 += i;

3132 }

3133 else if (i < 4) {

3134 i += 28;

3135 b2 += i;

3136 m2 += i;

3137 s2 += i;

3138 }

3139 if (b2 > 0)

3140 b = lshift(b, b2);

3141 if (s2 > 0)

3142 S = lshift(S, s2);

3143 if (k_check) {

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

3145 k--;

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

3147 if (leftright)

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

3149 ilim = ilim1;

3150 }

3151 }

3152 if (ilim <= 0 && (mode == 3 \|\| mode == 5)) {

3153 S = multadd(S, 5, 0);

3154 if (ilim < 0 \|\| cmp(b, S) < 0 \|\| ((cmp(b, S) == 0) && !bias_roun d_up)) {

3155 /* no digits, fcvt style */

3156 no_digits:

3157 k = -1 - ndigits;

3158 goto ret;

3159 }

3160 one_digit:

3161 *s++ = '1';

3162 k++;

3163 goto ret;

3164 }

3165 if (leftright) {

3166 if (m2 > 0)

3167 mhi = lshift(mhi, m2);

3168

3169 /* Compute mlo -- check for special case

3170 * that d is a normalized power of 2.

3171 */

3172

3173 mlo = mhi;

3174 if (spec_case) {

3175 mhi = Balloc(mhi->k);

3176 Bcopy(mhi, mlo);

3177 mhi = lshift(mhi, Log2P);

3178 }

3179

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

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

3182 /* Do we yet have the shortest decimal string

3183 * that will round to d?

3184 */

3185 j = cmp(b, mlo);

3186 delta = diff(S, mhi);

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

3188 Bfree(delta);

3189 #ifndef ROUND_BIASED

3190 if (j1 == 0 && mode != 1 && !(word1(d) & 1)

3191 #ifdef Honor_FLT_ROUNDS

3192 && rounding >= 1

3193 #endif

3194 ) {

3195 if (dig == '9')

3196 goto round_9_up;

3197 if (j > 0)

3198 dig++;

3199 #ifdef SET_INEXACT

3200 else if (!b->x[0] && b->wds <= 1)

3201 inexact = 0;

3202 #endif

3203 *s++ = dig;

3204 goto ret;

3205 }

3206 #endif

3207 if (j < 0 \|\| (j == 0 && mode != 1

3208 #ifndef ROUND_BIASED

3209 && !(word1(d) & 1)

3210 #endif

3211 )) {

3212 if (!b->x[0] && b->wds <= 1) {

3213 #ifdef SET_INEXACT

3214 inexact = 0;

3215 #endif

3216 goto accept_dig;

3217 }

3218 #ifdef Honor_FLT_ROUNDS

3219 if (mode > 1)

3220 switch(rounding) {

3221 case 0: goto accept_dig;

3222 case 2: goto keep_dig;

3223 }

3224 #endif /Honor_FLT_ROUNDS/

3225 if (j1 > 0) {

3226 b = lshift(b, 1);

3227 j1 = cmp(b, S);

3228 if ((j1 > 0 \|\| (j1 == 0 && ((dig & 1) \|\| bias_round_up)))

3229 && dig++ == '9')

3230 goto round_9_up;

3231 }

3232 accept_dig:

3233 *s++ = dig;

3234 goto ret;

3235 }

3236 if (j1 > 0) {

3237 #ifdef Honor_FLT_ROUNDS

3238 if (!rounding)

3239 goto accept_dig;

3240 #endif

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

3242 round_9_up:

3243 *s++ = '9';

3244 goto roundoff;

3245 }

3246 *s++ = dig + 1;

3247 goto ret;

3248 }

3249 #ifdef Honor_FLT_ROUNDS

3250 keep_dig:

3251 #endif

3252 *s++ = dig;

3253 if (i == ilim)

3254 break;

3255 b = multadd(b, 10, 0);

3256 if (mlo == mhi)

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

3258 else {

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

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

3261 }

3262 }

3263 }

3264 else

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

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

3267 if (!b->x[0] && b->wds <= 1) {

3268 #ifdef SET_INEXACT

3269 inexact = 0;

3270 #endif

3271 goto ret;

3272 }

3273 if (i >= ilim)

3274 break;

3275 b = multadd(b, 10, 0);

3276 }

3277

3278 /* Round off last digit */

3279

3280 #ifdef Honor_FLT_ROUNDS

3281 switch(rounding) {

3282 case 0: goto trimzeros;

3283 case 2: goto roundoff;

3284 }

3285 #endif

3286 b = lshift(b, 1);

3287 j = cmp(b, S);

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

3289 roundoff:

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

3291 if (s == s0) {

3292 k++;

3293 *s++ = '1';

3294 goto ret;

3295 }

3296 ++*s++;

3297 }

3298 else {

3299 /* trimzeros: (never used) */

3300 while(*--s == '0');

3301 s++;

3302 }

3303 ret:

3304 Bfree(S);

3305 if (mhi) {

3306 if (mlo && mlo != mhi)

3307 Bfree(mlo);

3308 Bfree(mhi);

3309 }

3310 ret1:

3311 #ifdef SET_INEXACT

3312 if (inexact) {

3313 if (!oldinexact) {

3314 word0(d) = Exp_1 + (70 << Exp_shift);

3315 word1(d) = 0;

3316 dval(d) += 1.;

3317 }

3318 }

3319 else if (!oldinexact)

3320 clear_inexact();

3321 #endif

3322 Bfree(b);

3323 *s = 0;

3324 *decpt = k + 1;

3325 if (rve)

3326 *rve = s;

3327 return s0;

3328 }

3329 #ifdef __cplusplus

3330 }

3331 #endif

OLD	NEW

« no previous file with comments | « src/third_party/dtoa/COPYING ('k') | test/cctest/test-conversions.cc » ('j') | no next file with comments »