third_party/dmg_fp/dtoa.cc - Issue 10963: Fix the mac build to add dmg_fp to the libbase.a.

Side by Side Diff: third_party/dmg_fp/dtoa.cc

Issue 10963: Fix the mac build to add dmg_fp to the libbase.a. (Closed)

Patch Set: Created 12 years, 1 month ago

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

Jump to:

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. Unless Trust_FLT_ROUNDS

76 * is also #defined, fegetround() will be queried for the rounding mode.

77 * Note that both FLT_ROUNDS and fegetround() are specified by the C99

78 * standard (and are specified to be consistent, with fesetround()

79 * affecting the value of FLT_ROUNDS), but that some (Linux) systems

80 * do not work correctly in this regard, so using fegetround() is more

81 * portable than using FLT_FOUNDS directly.

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

83 * and Honor_FLT_ROUNDS is not #defined.

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

85 * that use extended-precision instructions to compute rounded

86 * products and quotients) with IBM.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

102 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.

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

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

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

106 * directly -- and assumed always to succeed.

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

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

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

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

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

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

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

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

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

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

117 * pointers, PRIVATE_MEM >= 7112 appears adequate.

118 * #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK

119 * #defined automatically on IEEE systems. On such systems,

120 * when INFNAN_CHECK is #defined, strtod checks

121 * for Infinity and NaN (case insensitively). On some systems

122 * (e.g., some HP systems), it may be necessary to #define NAN_WORD0

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

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

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

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

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

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

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

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

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

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

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

134 * and NAN_WORD1 are used instead.

135 * #define MULTIPLE_THREADS if the system offers preemptively scheduled

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

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

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

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

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

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

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

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

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

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

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

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

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

149 * Sudden_Underflow.

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

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

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

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

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

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

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

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

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

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

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

161 * suitable wrapper, that defines two functions,

162 * int get_inexact(void);

163 * void clear_inexact(void);

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

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

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

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

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

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

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

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

172 */

173

174 #ifndef Long

175 #define Long long

176 #endif

177 #ifndef ULong

178 typedef unsigned Long ULong;

179 #endif

180

181 #ifdef DEBUG

182 #include "stdio.h"

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

184 #endif

185

186 #include "stdlib.h"

187 #include "string.h"

188

189 #ifdef USE_LOCALE

190 #include "locale.h"

191 #endif

192

193 #ifdef MALLOC

194 #ifdef KR_headers

195 extern char *MALLOC();

196 #else

197 extern void *MALLOC(size_t);

198 #endif

199 #else

200 #define MALLOC malloc

201 #endif

202

203 #ifndef Omit_Private_Memory

204 #ifndef PRIVATE_MEM

205 #define PRIVATE_MEM 2304

206 #endif

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

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

209 #endif

210

211 #undef IEEE_Arith

212 #undef Avoid_Underflow

213 #ifdef IEEE_MC68k

214 #define IEEE_Arith

215 #endif

216 #ifdef IEEE_8087

217 #define IEEE_Arith

218 #endif

219

220 #ifdef IEEE_Arith

221 #ifndef NO_INFNAN_CHECK

222 #undef INFNAN_CHECK

223 #define INFNAN_CHECK

224 #endif

225 #else

226 #undef INFNAN_CHECK

227 #endif

228

229 #include "errno.h"

230

231 #ifdef Bad_float_h

232

233 #ifdef IEEE_Arith

234 #define DBL_DIG 15

235 #define DBL_MAX_10_EXP 308

236 #define DBL_MAX_EXP 1024

237 #define FLT_RADIX 2

238 #endif /IEEE_Arith/

239

240 #ifdef IBM

241 #define DBL_DIG 16

242 #define DBL_MAX_10_EXP 75

243 #define DBL_MAX_EXP 63

244 #define FLT_RADIX 16

245 #define DBL_MAX 7.2370055773322621e+75

246 #endif

247

248 #ifdef VAX

249 #define DBL_DIG 16

250 #define DBL_MAX_10_EXP 38

251 #define DBL_MAX_EXP 127

252 #define FLT_RADIX 2

253 #define DBL_MAX 1.7014118346046923e+38

254 #endif

255

256 #ifndef LONG_MAX

257 #define LONG_MAX 2147483647

258 #endif

259

260 #else /* ifndef Bad_float_h */

261 #include "float.h"

262 #endif /* Bad_float_h */

263

264 #ifndef __MATH_H__

265 #include "math.h"

266 #endif

267

268 namespace dmg_fp {

269

270 #ifndef CONST

271 #ifdef KR_headers

272 #define CONST /* blank */

273 #else

274 #define CONST const

275 #endif

276 #endif

277

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

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

280 #endif

281

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

283

284 #ifdef YES_ALIAS

285 #define dval(x) x

286 #ifdef IEEE_8087

287 #define word0(x) ((ULong *)&x)[1]

288 #define word1(x) ((ULong *)&x)[0]

289 #else

290 #define word0(x) ((ULong *)&x)[0]

291 #define word1(x) ((ULong *)&x)[1]

292 #endif

293 #else

294 #ifdef IEEE_8087

295 #define word0(x) ((U*)&x)->L[1]

296 #define word1(x) ((U*)&x)->L[0]

297 #else

298 #define word0(x) ((U*)&x)->L[0]

299 #define word1(x) ((U*)&x)->L[1]

300 #endif

301 #define dval(x) ((U*)&x)->d

302 #endif

303

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

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

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

307 */

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

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

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

311 #else

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

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

314 #endif

315

316 /* #define P DBL_MANT_DIG */

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

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

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

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

321

322 #ifdef IEEE_Arith

323 #define Exp_shift 20

324 #define Exp_shift1 20

325 #define Exp_msk1 0x100000

326 #define Exp_msk11 0x100000

327 #define Exp_mask 0x7ff00000

328 #define P 53

329 #define Bias 1023

330 #define Emin (-1022)

331 #define Exp_1 0x3ff00000

332 #define Exp_11 0x3ff00000

333 #define Ebits 11

334 #define Frac_mask 0xfffff

335 #define Frac_mask1 0xfffff

336 #define Ten_pmax 22

337 #define Bletch 0x10

338 #define Bndry_mask 0xfffff

339 #define Bndry_mask1 0xfffff

340 #define LSB 1

341 #define Sign_bit 0x80000000

342 #define Log2P 1

343 #define Tiny0 0

344 #define Tiny1 1

345 #define Quick_max 14

346 #define Int_max 14

347 #ifndef NO_IEEE_Scale

348 #define Avoid_Underflow

349 #ifdef Flush_Denorm /* debugging option */

350 #undef Sudden_Underflow

351 #endif

352 #endif

353

354 #ifndef Flt_Rounds

355 #ifdef FLT_ROUNDS

356 #define Flt_Rounds FLT_ROUNDS

357 #else

358 #define Flt_Rounds 1

359 #endif

360 #endif /Flt_Rounds/

361

362 #ifdef Honor_FLT_ROUNDS

363 #undef Check_FLT_ROUNDS

364 #define Check_FLT_ROUNDS

365 #else

366 #define Rounding Flt_Rounds

367 #endif

368

369 #else /* ifndef IEEE_Arith */

370 #undef Check_FLT_ROUNDS

371 #undef Honor_FLT_ROUNDS

372 #undef SET_INEXACT

373 #undef Sudden_Underflow

374 #define Sudden_Underflow

375 #ifdef IBM

376 #undef Flt_Rounds

377 #define Flt_Rounds 0

378 #define Exp_shift 24

379 #define Exp_shift1 24

380 #define Exp_msk1 0x1000000

381 #define Exp_msk11 0x1000000

382 #define Exp_mask 0x7f000000

383 #define P 14

384 #define Bias 65

385 #define Exp_1 0x41000000

386 #define Exp_11 0x41000000

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

388 #define Frac_mask 0xffffff

389 #define Frac_mask1 0xffffff

390 #define Bletch 4

391 #define Ten_pmax 22

392 #define Bndry_mask 0xefffff

393 #define Bndry_mask1 0xffffff

394 #define LSB 1

395 #define Sign_bit 0x80000000

396 #define Log2P 4

397 #define Tiny0 0x100000

398 #define Tiny1 0

399 #define Quick_max 14

400 #define Int_max 15

401 #else /* VAX */

402 #undef Flt_Rounds

403 #define Flt_Rounds 1

404 #define Exp_shift 23

405 #define Exp_shift1 7

406 #define Exp_msk1 0x80

407 #define Exp_msk11 0x800000

408 #define Exp_mask 0x7f80

409 #define P 56

410 #define Bias 129

411 #define Exp_1 0x40800000

412 #define Exp_11 0x4080

413 #define Ebits 8

414 #define Frac_mask 0x7fffff

415 #define Frac_mask1 0xffff007f

416 #define Ten_pmax 24

417 #define Bletch 2

418 #define Bndry_mask 0xffff007f

419 #define Bndry_mask1 0xffff007f

420 #define LSB 0x10000

421 #define Sign_bit 0x8000

422 #define Log2P 1

423 #define Tiny0 0x80

424 #define Tiny1 0

425 #define Quick_max 15

426 #define Int_max 15

427 #endif /* IBM, VAX */

428 #endif /* IEEE_Arith */

429

430 #ifndef IEEE_Arith

431 #define ROUND_BIASED

432 #endif

433

434 #ifdef RND_PRODQUOT

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

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

437 #ifdef KR_headers

438 extern double rnd_prod(), rnd_quot();

439 #else

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

441 #endif

442 #else

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

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

445 #endif

446

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

448 #define Big1 0xffffffff

449

450 #ifndef Pack_32

451 #define Pack_32

452 #endif

453

454 #ifdef KR_headers

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

456 #else

457 #define FFFFFFFF 0xffffffffUL

458 #endif

459

460 #ifdef NO_LONG_LONG

461 #undef ULLong

462 #ifdef Just_16

463 #undef Pack_32

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

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

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

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

468 */

469 #endif

470 #else /* long long available */

471 #ifndef Llong

472 #define Llong long long

473 #endif

474 #ifndef ULLong

475 #define ULLong unsigned Llong

476 #endif

477 #endif /* NO_LONG_LONG */

478

479 #ifndef MULTIPLE_THREADS

480 #define ACQUIRE_DTOA_LOCK(n) /nothing/

481 #define FREE_DTOA_LOCK(n) /nothing/

482 #endif

483

484 #define Kmax 15

485

486 double strtod(const char s00, char *se);

487 char *dtoa(double d, int mode, int ndigits,

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

489

490 struct

491 Bigint {

492 struct Bigint *next;

493 int k, maxwds, sign, wds;

494 ULong x[1];

495 };

496

497 typedef struct Bigint Bigint;

498

499 static Bigint *freelist[Kmax+1];

500

501 static Bigint *

502 Balloc

503 #ifdef KR_headers

504 (k) int k;

505 #else

506 (int k)

507 #endif

508 {

509 int x;

510 Bigint *rv;

511 #ifndef Omit_Private_Memory

512 unsigned int len;

513 #endif

514

515 ACQUIRE_DTOA_LOCK(0);

516 if (rv = freelist[k]) {

517 freelist[k] = rv->next;

518 }

519 else {

520 x = 1 << k;

521 #ifdef Omit_Private_Memory

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

523 #else

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

525 /sizeof(double);

526 if (pmem_next - private_mem + len <= PRIVATE_mem) {

527 rv = (Bigint*)pmem_next;

528 pmem_next += len;

529 }

530 else

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

532 #endif

533 rv->k = k;

534 rv->maxwds = x;

535 }

536 FREE_DTOA_LOCK(0);

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

538 return rv;

539 }

540

541 static void

542 Bfree

543 #ifdef KR_headers

544 (v) Bigint *v;

545 #else

546 (Bigint *v)

547 #endif

548 {

549 if (v) {

550 ACQUIRE_DTOA_LOCK(0);

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

552 freelist[v->k] = v;

553 FREE_DTOA_LOCK(0);

554 }

555 }

556

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

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

559

560 static Bigint *

561 multadd

562 #ifdef KR_headers

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

564 #else

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

566 #endif

567 {

568 int i, wds;

569 #ifdef ULLong

570 ULong *x;

571 ULLong carry, y;

572 #else

573 ULong carry, *x, y;

574 #ifdef Pack_32

575 ULong xi, z;

576 #endif

577 #endif

578 Bigint *b1;

579

580 wds = b->wds;

581 x = b->x;

582 i = 0;

583 carry = a;

584 do {

585 #ifdef ULLong

586 y = x (ULLong)m + carry;

587 carry = y >> 32;

588 *x++ = y & FFFFFFFF;

589 #else

590 #ifdef Pack_32

591 xi = *x;

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

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

594 carry = z >> 16;

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

596 #else

597 y = x m + carry;

598 carry = y >> 16;

599 *x++ = y & 0xffff;

600 #endif

601 #endif

602 }

603 while(++i < wds);

604 if (carry) {

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

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

607 Bcopy(b1, b);

608 Bfree(b);

609 b = b1;

610 }

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

612 b->wds = wds;

613 }

614 return b;

615 }

616

617 static Bigint *

618 s2b

619 #ifdef KR_headers

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

621 #else

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

623 #endif

624 {

625 Bigint *b;

626 int i, k;

627 Long x, y;

628

629 x = (nd + 8) / 9;

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

631 #ifdef Pack_32

632 b = Balloc(k);

633 b->x[0] = y9;

634 b->wds = 1;

635 #else

636 b = Balloc(k+1);

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

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

639 #endif

640

641 i = 9;

642 if (9 < nd0) {

643 s += 9;

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

645 while(++i < nd0);

646 s++;

647 }

648 else

649 s += 10;

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

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

652 return b;

653 }

654

655 static int

656 hi0bits

657 #ifdef KR_headers

658 (x) register ULong x;

659 #else

660 (register ULong x)

661 #endif

662 {

663 register int k = 0;

664

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

666 k = 16;

667 x <<= 16;

668 }

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

670 k += 8;

671 x <<= 8;

672 }

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

674 k += 4;

675 x <<= 4;

676 }

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

678 k += 2;

679 x <<= 2;

680 }

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

682 k++;

683 if (!(x & 0x40000000))

684 return 32;

685 }

686 return k;

687 }

688

689 static int

690 lo0bits

691 #ifdef KR_headers

692 (y) ULong *y;

693 #else

694 (ULong *y)

695 #endif

696 {

697 register int k;

698 register ULong x = *y;

699

700 if (x & 7) {

701 if (x & 1)

702 return 0;

703 if (x & 2) {

704 *y = x >> 1;

705 return 1;

706 }

707 *y = x >> 2;

708 return 2;

709 }

710 k = 0;

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

712 k = 16;

713 x >>= 16;

714 }

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

716 k += 8;

717 x >>= 8;

718 }

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

720 k += 4;

721 x >>= 4;

722 }

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

724 k += 2;

725 x >>= 2;

726 }

727 if (!(x & 1)) {

728 k++;

729 x >>= 1;

730 if (!x)

731 return 32;

732 }

733 *y = x;

734 return k;

735 }

736

737 static Bigint *

738 i2b

739 #ifdef KR_headers

740 (i) int i;

741 #else

742 (int i)

743 #endif

744 {

745 Bigint *b;

746

747 b = Balloc(1);

748 b->x[0] = i;

749 b->wds = 1;

750 return b;

751 }

752

753 static Bigint *

754 mult

755 #ifdef KR_headers

756 (a, b) Bigint a, b;

757 #else

758 (Bigint a, Bigint b)

759 #endif

760 {

761 Bigint *c;

762 int k, wa, wb, wc;

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

764 ULong y;

765 #ifdef ULLong

766 ULLong carry, z;

767 #else

768 ULong carry, z;

769 #ifdef Pack_32

770 ULong z2;

771 #endif

772 #endif

773

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

775 c = a;

776 a = b;

777 b = c;

778 }

779 k = a->k;

780 wa = a->wds;

781 wb = b->wds;

782 wc = wa + wb;

783 if (wc > a->maxwds)

784 k++;

785 c = Balloc(k);

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

787 *x = 0;

788 xa = a->x;

789 xae = xa + wa;

790 xb = b->x;

791 xbe = xb + wb;

792 xc0 = c->x;

793 #ifdef ULLong

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

795 if (y = *xb++) {

796 x = xa;

797 xc = xc0;

798 carry = 0;

799 do {

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

801 carry = z >> 32;

802 *xc++ = z & FFFFFFFF;

803 }

804 while(x < xae);

805 *xc = carry;

806 }

807 }

808 #else

809 #ifdef Pack_32

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

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

812 x = xa;

813 xc = xc0;

814 carry = 0;

815 do {

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

817 carry = z >> 16;

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

819 carry = z2 >> 16;

820 Storeinc(xc, z2, z);

821 }

822 while(x < xae);

823 *xc = carry;

824 }

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

826 x = xa;

827 xc = xc0;

828 carry = 0;

829 z2 = *xc;

830 do {

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

832 carry = z >> 16;

833 Storeinc(xc, z, z2);

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

835 carry = z2 >> 16;

836 }

837 while(x < xae);

838 *xc = z2;

839 }

840 }

841 #else

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

843 if (y = *xb++) {

844 x = xa;

845 xc = xc0;

846 carry = 0;

847 do {

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

849 carry = z >> 16;

850 *xc++ = z & 0xffff;

851 }

852 while(x < xae);

853 *xc = carry;

854 }

855 }

856 #endif

857 #endif

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

859 c->wds = wc;

860 return c;

861 }

862

863 static Bigint *p5s;

864

865 static Bigint *

866 pow5mult

867 #ifdef KR_headers

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

869 #else

870 (Bigint *b, int k)

871 #endif

872 {

873 Bigint b1, p5, *p51;

874 int i;

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

876

877 if (i = k & 3)

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

879

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

881 return b;

882 if (!(p5 = p5s)) {

883 /* first time */

884 #ifdef MULTIPLE_THREADS

885 ACQUIRE_DTOA_LOCK(1);

886 if (!(p5 = p5s)) {

887 p5 = p5s = i2b(625);

888 p5->next = 0;

889 }

890 FREE_DTOA_LOCK(1);

891 #else

892 p5 = p5s = i2b(625);

893 p5->next = 0;

894 #endif

895 }

896 for(;;) {

897 if (k & 1) {

898 b1 = mult(b, p5);

899 Bfree(b);

900 b = b1;

901 }

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

903 break;

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

905 #ifdef MULTIPLE_THREADS

906 ACQUIRE_DTOA_LOCK(1);

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

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

909 p51->next = 0;

910 }

911 FREE_DTOA_LOCK(1);

912 #else

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

914 p51->next = 0;

915 #endif

916 }

917 p5 = p51;

918 }

919 return b;

920 }

921

922 static Bigint *

923 lshift

924 #ifdef KR_headers

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

926 #else

927 (Bigint *b, int k)

928 #endif

929 {

930 int i, k1, n, n1;

931 Bigint *b1;

932 ULong x, x1, *xe, z;

933

934 #ifdef Pack_32

935 n = k >> 5;

936 #else

937 n = k >> 4;

938 #endif

939 k1 = b->k;

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

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

942 k1++;

943 b1 = Balloc(k1);

944 x1 = b1->x;

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

946 *x1++ = 0;

947 x = b->x;

948 xe = x + b->wds;

949 #ifdef Pack_32

950 if (k &= 0x1f) {

951 k1 = 32 - k;

952 z = 0;

953 do {

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

955 z = *x++ >> k1;

956 }

957 while(x < xe);

958 if (*x1 = z)

959 ++n1;

960 }

961 #else

962 if (k &= 0xf) {

963 k1 = 16 - k;

964 z = 0;

965 do {

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

967 z = *x++ >> k1;

968 }

969 while(x < xe);

970 if (*x1 = z)

971 ++n1;

972 }

973 #endif

974 else do

975 x1++ = x++;

976 while(x < xe);

977 b1->wds = n1 - 1;

978 Bfree(b);

979 return b1;

980 }

981

982 static int

983 cmp

984 #ifdef KR_headers

985 (a, b) Bigint a, b;

986 #else

987 (Bigint a, Bigint b)

988 #endif

989 {

990 ULong xa, xa0, xb, xb0;

991 int i, j;

992

993 i = a->wds;

994 j = b->wds;

995 #ifdef DEBUG

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

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

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

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

1000 #endif

1001 if (i -= j)

1002 return i;

1003 xa0 = a->x;

1004 xa = xa0 + j;

1005 xb0 = b->x;

1006 xb = xb0 + j;

1007 for(;;) {

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

1009 return xa < xb ? -1 : 1;

1010 if (xa <= xa0)

1011 break;

1012 }

1013 return 0;

1014 }

1015

1016 static Bigint *

1017 diff

1018 #ifdef KR_headers

1019 (a, b) Bigint a, b;

1020 #else

1021 (Bigint a, Bigint b)

1022 #endif

1023 {

1024 Bigint *c;

1025 int i, wa, wb;

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

1027 #ifdef ULLong

1028 ULLong borrow, y;

1029 #else

1030 ULong borrow, y;

1031 #ifdef Pack_32

1032 ULong z;

1033 #endif

1034 #endif

1035

1036 i = cmp(a,b);

1037 if (!i) {

1038 c = Balloc(0);

1039 c->wds = 1;

1040 c->x[0] = 0;

1041 return c;

1042 }

1043 if (i < 0) {

1044 c = a;

1045 a = b;

1046 b = c;

1047 i = 1;

1048 }

1049 else

1050 i = 0;

1051 c = Balloc(a->k);

1052 c->sign = i;

1053 wa = a->wds;

1054 xa = a->x;

1055 xae = xa + wa;

1056 wb = b->wds;

1057 xb = b->x;

1058 xbe = xb + wb;

1059 xc = c->x;

1060 borrow = 0;

1061 #ifdef ULLong

1062 do {

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

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

1065 *xc++ = y & FFFFFFFF;

1066 }

1067 while(xb < xbe);

1068 while(xa < xae) {

1069 y = *xa++ - borrow;

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

1071 *xc++ = y & FFFFFFFF;

1072 }

1073 #else

1074 #ifdef Pack_32

1075 do {

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

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

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

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

1080 Storeinc(xc, z, y);

1081 }

1082 while(xb < xbe);

1083 while(xa < xae) {

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

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

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

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

1088 Storeinc(xc, z, y);

1089 }

1090 #else

1091 do {

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

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

1094 *xc++ = y & 0xffff;

1095 }

1096 while(xb < xbe);

1097 while(xa < xae) {

1098 y = *xa++ - borrow;

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

1100 *xc++ = y & 0xffff;

1101 }

1102 #endif

1103 #endif

1104 while(!*--xc)

1105 wa--;

1106 c->wds = wa;

1107 return c;

1108 }

1109

1110 static double

1111 ulp

1112 #ifdef KR_headers

1113 (x) double x;

1114 #else

1115 (double x)

1116 #endif

1117 {

1118 register Long L;

1119 double a;

1120

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

1122 #ifndef Avoid_Underflow

1123 #ifndef Sudden_Underflow

1124 if (L > 0) {

1125 #endif

1126 #endif

1127 #ifdef IBM

1128 L \|= Exp_msk1 >> 4;

1129 #endif

1130 word0(a) = L;

1131 word1(a) = 0;

1132 #ifndef Avoid_Underflow

1133 #ifndef Sudden_Underflow

1134 }

1135 else {

1136 L = -L >> Exp_shift;

1137 if (L < Exp_shift) {

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

1139 word1(a) = 0;

1140 }

1141 else {

1142 word0(a) = 0;

1143 L -= Exp_shift;

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

1145 }

1146 }

1147 #endif

1148 #endif

1149 return dval(a);

1150 }

1151

1152 static double

1153 b2d

1154 #ifdef KR_headers

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

1156 #else

1157 (Bigint a, int e)

1158 #endif

1159 {

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

1161 int k;

1162 double d;

1163 #ifdef VAX

1164 ULong d0, d1;

1165 #else

1166 #define d0 word0(d)

1167 #define d1 word1(d)

1168 #endif

1169

1170 xa0 = a->x;

1171 xa = xa0 + a->wds;

1172 y = *--xa;

1173 #ifdef DEBUG

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

1175 #endif

1176 k = hi0bits(y);

1177 *e = 32 - k;

1178 #ifdef Pack_32

1179 if (k < Ebits) {

1180 d0 = Exp_1 \| y >> Ebits - k;

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

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

1183 goto ret_d;

1184 }

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

1186 if (k -= Ebits) {

1187 d0 = Exp_1 \| y << k \| z >> 32 - k;

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

1189 d1 = z << k \| y >> 32 - k;

1190 }

1191 else {

1192 d0 = Exp_1 \| y;

1193 d1 = z;

1194 }

1195 #else

1196 if (k < Ebits + 16) {

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

1198 d0 = Exp_1 \| y << k - Ebits \| z >> Ebits + 16 - k;

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

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

1201 d1 = z << k + 16 - Ebits \| w << k - Ebits \| y >> 16 + Ebits - k;

1202 goto ret_d;

1203 }

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

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

1206 k -= Ebits + 16;

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

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

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

1210 #endif

1211 ret_d:

1212 #ifdef VAX

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

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

1215 #else

1216 #undef d0

1217 #undef d1

1218 #endif

1219 return dval(d);

1220 }

1221

1222 static Bigint *

1223 d2b

1224 #ifdef KR_headers

1225 (d, e, bits) double d; int e, bits;

1226 #else

1227 (double d, int e, int bits)

1228 #endif

1229 {

1230 Bigint *b;

1231 int de, k;

1232 ULong *x, y, z;

1233 #ifndef Sudden_Underflow

1234 int i;

1235 #endif

1236 #ifdef VAX

1237 ULong d0, d1;

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

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

1240 #else

1241 #define d0 word0(d)

1242 #define d1 word1(d)

1243 #endif

1244

1245 #ifdef Pack_32

1246 b = Balloc(1);

1247 #else

1248 b = Balloc(2);

1249 #endif

1250 x = b->x;

1251

1252 z = d0 & Frac_mask;

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

1254 #ifdef Sudden_Underflow

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

1256 #ifndef IBM

1257 z \|= Exp_msk11;

1258 #endif

1259 #else

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

1261 z \|= Exp_msk1;

1262 #endif

1263 #ifdef Pack_32

1264 if (y = d1) {

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

1266 x[0] = y \| z << 32 - k;

1267 z >>= k;

1268 }

1269 else

1270 x[0] = y;

1271 #ifndef Sudden_Underflow

1272 i =

1273 #endif

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

1275 }

1276 else {

1277 #ifdef DEBUG

1278 if (!z)

1279 Bug("Zero passed to d2b");

1280 #endif

1281 k = lo0bits(&z);

1282 x[0] = z;

1283 #ifndef Sudden_Underflow

1284 i =

1285 #endif

1286 b->wds = 1;

1287 k += 32;

1288 }

1289 #else

1290 if (y = d1) {

1291 if (k = lo0bits(&y))

1292 if (k >= 16) {

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

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

1295 x[2] = z >> k;

1296 i = 2;

1297 }

1298 else {

1299 x[0] = y & 0xffff;

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

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

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

1303 i = 3;

1304 }

1305 else {

1306 x[0] = y & 0xffff;

1307 x[1] = y >> 16;

1308 x[2] = z & 0xffff;

1309 x[3] = z >> 16;

1310 i = 3;

1311 }

1312 }

1313 else {

1314 #ifdef DEBUG

1315 if (!z)

1316 Bug("Zero passed to d2b");

1317 #endif

1318 k = lo0bits(&z);

1319 if (k >= 16) {

1320 x[0] = z;

1321 i = 0;

1322 }

1323 else {

1324 x[0] = z & 0xffff;

1325 x[1] = z >> 16;

1326 i = 1;

1327 }

1328 k += 32;

1329 }

1330 while(!x[i])

1331 --i;

1332 b->wds = i + 1;

1333 #endif

1334 #ifndef Sudden_Underflow

1335 if (de) {

1336 #endif

1337 #ifdef IBM

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

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

1340 #else

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

1342 *bits = P - k;

1343 #endif

1344 #ifndef Sudden_Underflow

1345 }

1346 else {

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

1348 #ifdef Pack_32

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

1350 #else

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

1352 #endif

1353 }

1354 #endif

1355 return b;

1356 }

1357 #undef d0

1358 #undef d1

1359

1360 static double

1361 ratio

1362 #ifdef KR_headers

1363 (a, b) Bigint a, b;

1364 #else

1365 (Bigint a, Bigint b)

1366 #endif

1367 {

1368 double da, db;

1369 int k, ka, kb;

1370

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

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

1373 #ifdef Pack_32

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

1375 #else

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

1377 #endif

1378 #ifdef IBM

1379 if (k > 0) {

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

1381 if (k &= 3)

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

1383 }

1384 else {

1385 k = -k;

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

1387 if (k &= 3)

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

1389 }

1390 #else

1391 if (k > 0)

1392 word0(da) += k*Exp_msk1;

1393 else {

1394 k = -k;

1395 word0(db) += k*Exp_msk1;

1396 }

1397 #endif

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

1399 }

1400

1401 static CONST double

1402 tens[] = {

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

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

1405 1e20, 1e21, 1e22

1406 #ifdef VAX

1407 , 1e23, 1e24

1408 #endif

1409 };

1410

1411 static CONST double

1412 #ifdef IEEE_Arith

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

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

1415 #ifdef Avoid_Underflow

1416 9007199254740992.*9007199254740992.e-256

1417 /* = 2^106 * 1e-256 */

1418 #else

1419 1e-256

1420 #endif

1421 };

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

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

1424 #define Scale_Bit 0x10

1425 #define n_bigtens 5

1426 #else

1427 #ifdef IBM

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

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

1430 #define n_bigtens 3

1431 #else

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

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

1434 #define n_bigtens 2

1435 #endif

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 /* allow optional initial 0x or 0X */

1487 while((c = (CONST unsigned char)(s+1)) && c <= ' ')

1488 ++s;

1489 if (s[1] == '0' && (s[2] == 'x' \|\| s[2] == 'X'))

1490 s += 2;

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

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

1493 c -= '0';

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

1495 c += 10 - 'a';

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

1497 c += 10 - 'A';

1498 else if (c <= ' ') {

1499 if (udx0 && havedig) {

1500 udx0 = 0;

1501 xshift = 1;

1502 }

1503 continue;

1504 }

1505 #ifdef GDTOA_NON_PEDANTIC_NANCHECK

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

1507 *sp = s + 1;

1508 break;

1509 }

1510 else

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

1512 #else

1513 else {

1514 do {

1515 if (/(/ c == ')') {

1516 *sp = s + 1;

1517 break;

1518 }

1519 } while(c = *++s);

1520 break;

1521 }

1522 #endif

1523 havedig = 1;

1524 if (xshift) {

1525 xshift = 0;

1526 x[0] = x[1];

1527 x[1] = 0;

1528 }

1529 if (udx0)

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

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

1532 }

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

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

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

1536 }

1537 }

1538 #endif /No_Hex_NaN/

1539 #endif /* INFNAN_CHECK */

1540

1541 double

1542 strtod

1543 #ifdef KR_headers

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

1545 #else

1546 (CONST char s00, char *se)

1547 #endif

1548 {

1549 #ifdef Avoid_Underflow

1550 int scale;

1551 #endif

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

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

1554 CONST char s, s0, *s1;

1555 double aadj, aadj1, adj, rv, rv0;

1556 Long L;

1557 ULong y, z;

1558 Bigint bb, bb1, bd, bd0, bs, delta;

1559 #ifdef SET_INEXACT

1560 int inexact, oldinexact;

1561 #endif

1562 #ifdef Honor_FLT_ROUNDS /{/

1563 int Rounding;

1564 #ifdef Trust_FLT_ROUNDS /{{ only define this if FLT_ROUNDS really works! /

1565 Rounding = Flt_Rounds;

1566 #else /}{/

1567 Rounding = 1;

1568 switch(fegetround()) {

1569 case FE_TOWARDZERO: Rounding = 0; break;

1570 case FE_UPWARD: Rounding = 2; break;

1571 case FE_DOWNWARD: Rounding = 3;

1572 }

1573 #endif /}}/

1574 #endif /}/

1575 #ifdef USE_LOCALE

1576 CONST char *s2;

1577 #endif

1578

1579 sign = nz0 = nz = 0;

1580 dval(rv) = 0.;

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

1582 case '-':

1583 sign = 1;

1584 /* no break */

1585 case '+':

1586 if (*++s)

1587 goto break2;

1588 /* no break */

1589 case 0:

1590 goto ret0;

1591 case '\t':

1592 case '\n':

1593 case '\v':

1594 case '\f':

1595 case '\r':

1596 case ' ':

1597 continue;

1598 default:

1599 goto break2;

1600 }

1601 break2:

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

1603 nz0 = 1;

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

1605 if (!*s)

1606 goto ret;

1607 }

1608 s0 = s;

1609 y = z = 0;

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

1611 if (nd < 9)

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

1613 else if (nd < 16)

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

1615 nd0 = nd;

1616 #ifdef USE_LOCALE

1617 s1 = localeconv()->decimal_point;

1618 if (c == *s1) {

1619 c = '.';

1620 if (*++s1) {

1621 s2 = s;

1622 for(;;) {

1623 if (++s2 != s1) {

1624 c = 0;

1625 break;

1626 }

1627 if (!*++s1) {

1628 s = s2;

1629 break;

1630 }

1631 }

1632 }

1633 }

1634 #endif

1635 if (c == '.') {

1636 c = *++s;

1637 if (!nd) {

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

1639 nz++;

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

1641 s0 = s;

1642 nf += nz;

1643 nz = 0;

1644 goto have_dig;

1645 }

1646 goto dig_done;

1647 }

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

1649 have_dig:

1650 nz++;

1651 if (c -= '0') {

1652 nf += nz;

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

1654 if (nd++ < 9)

1655 y *= 10;

1656 else if (nd <= DBL_DIG + 1)

1657 z *= 10;

1658 if (nd++ < 9)

1659 y = 10*y + c;

1660 else if (nd <= DBL_DIG + 1)

1661 z = 10*z + c;

1662 nz = 0;

1663 }

1664 }

1665 }

1666 dig_done:

1667 e = 0;

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

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

1670 goto ret0;

1671 }

1672 s00 = s;

1673 esign = 0;

1674 switch(c = *++s) {

1675 case '-':

1676 esign = 1;

1677 case '+':

1678 c = *++s;

1679 }

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

1681 while(c == '0')

1682 c = *++s;

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

1684 L = c - '0';

1685 s1 = s;

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

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

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

1689 /* Avoid confusion from exponents

1690 * so large that e might overflow.

1691 */

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

1693 else

1694 e = (int)L;

1695 if (esign)

1696 e = -e;

1697 }

1698 else

1699 e = 0;

1700 }

1701 else

1702 s = s00;

1703 }

1704 if (!nd) {

1705 if (!nz && !nz0) {

1706 #ifdef INFNAN_CHECK

1707 /* Check for Nan and Infinity */

1708 switch(c) {

1709 case 'i':

1710 case 'I':

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

1712 --s;

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

1714 ++s;

1715 word0(rv) = 0x7ff00000;

1716 word1(rv) = 0;

1717 goto ret;

1718 }

1719 break;

1720 case 'n':

1721 case 'N':

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

1723 word0(rv) = NAN_WORD0;

1724 word1(rv) = NAN_WORD1;

1725 #ifndef No_Hex_NaN

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

1727 hexnan(&rv, &s);

1728 #endif

1729 goto ret;

1730 }

1731 }

1732 #endif /* INFNAN_CHECK */

1733 ret0:

1734 s = s00;

1735 sign = 0;

1736 }

1737 goto ret;

1738 }

1739 e1 = e -= nf;

1740

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

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

1743 * after is the integer represented by those digits times

1744 * 10*e /

1745

1746 if (!nd0)

1747 nd0 = nd;

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

1749 dval(rv) = y;

1750 if (k > 9) {

1751 #ifdef SET_INEXACT

1752 if (k > DBL_DIG)

1753 oldinexact = get_inexact();

1754 #endif

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

1756 }

1757 bd0 = 0;

1758 if (nd <= DBL_DIG

1759 #ifndef RND_PRODQUOT

1760 #ifndef Honor_FLT_ROUNDS

1761 && Flt_Rounds == 1

1762 #endif

1763 #endif

1764 ) {

1765 if (!e)

1766 goto ret;

1767 if (e > 0) {

1768 if (e <= Ten_pmax) {

1769 #ifdef VAX

1770 goto vax_ovfl_check;

1771 #else

1772 #ifdef Honor_FLT_ROUNDS

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

1774 if (sign) {

1775 rv = -rv;

1776 sign = 0;

1777 }

1778 #endif

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

1780 goto ret;

1781 #endif

1782 }

1783 i = DBL_DIG - nd;

1784 if (e <= Ten_pmax + i) {

1785 /* A fancier test would sometimes let us do

1786 * this for larger i values.

1787 */

1788 #ifdef Honor_FLT_ROUNDS

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

1790 if (sign) {

1791 rv = -rv;

1792 sign = 0;

1793 }

1794 #endif

1795 e -= i;

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

1797 #ifdef VAX

1798 /* VAX exponent range is so narrow we must

1799 * worry about overflow here...

1800 */

1801 vax_ovfl_check:

1802 word0(rv) -= P*Exp_msk1;

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

1804 if ((word0(rv) & Exp_mask)

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

1806 goto ovfl;

1807 word0(rv) += P*Exp_msk1;

1808 #else

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

1810 #endif

1811 goto ret;

1812 }

1813 }

1814 #ifndef Inaccurate_Divide

1815 else if (e >= -Ten_pmax) {

1816 #ifdef Honor_FLT_ROUNDS

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

1818 if (sign) {

1819 rv = -rv;

1820 sign = 0;

1821 }

1822 #endif

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

1824 goto ret;

1825 }

1826 #endif

1827 }

1828 e1 += nd - k;

1829

1830 #ifdef IEEE_Arith

1831 #ifdef SET_INEXACT

1832 inexact = 1;

1833 if (k <= DBL_DIG)

1834 oldinexact = get_inexact();

1835 #endif

1836 #ifdef Avoid_Underflow

1837 scale = 0;

1838 #endif

1839 #ifdef Honor_FLT_ROUNDS

1840 if (Rounding >= 2) {

1841 if (sign)

1842 Rounding = Rounding == 2 ? 0 : 2;

1843 else

1844 if (Rounding != 2)

1845 Rounding = 0;

1846 }

1847 #endif

1848 #endif /IEEE_Arith/

1849

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

1851

1852 if (e1 > 0) {

1853 if (i = e1 & 15)

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

1855 if (e1 &= ~15) {

1856 if (e1 > DBL_MAX_10_EXP) {

1857 ovfl:

1858 #ifndef NO_ERRNO

1859 errno = ERANGE;

1860 #endif

1861 /* Can't trust HUGE_VAL */

1862 #ifdef IEEE_Arith

1863 #ifdef Honor_FLT_ROUNDS

1864 switch(Rounding) {

1865 case 0: /* toward 0 */

1866 case 3: /* toward -infinity */

1867 word0(rv) = Big0;

1868 word1(rv) = Big1;

1869 break;

1870 default:

1871 word0(rv) = Exp_mask;

1872 word1(rv) = 0;

1873 }

1874 #else /Honor_FLT_ROUNDS/

1875 word0(rv) = Exp_mask;

1876 word1(rv) = 0;

1877 #endif /Honor_FLT_ROUNDS/

1878 #ifdef SET_INEXACT

1879 /* set overflow bit */

1880 dval(rv0) = 1e300;

1881 dval(rv0) *= dval(rv0);

1882 #endif

1883 #else /IEEE_Arith/

1884 word0(rv) = Big0;

1885 word1(rv) = Big1;

1886 #endif /IEEE_Arith/

1887 if (bd0)

1888 goto retfree;

1889 goto ret;

1890 }

1891 e1 >>= 4;

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

1893 if (e1 & 1)

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

1895 /* The last multiplication could overflow. */

1896 word0(rv) -= P*Exp_msk1;

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

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

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

1900 goto ovfl;

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

1902 /* set to largest number */

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

1904 word0(rv) = Big0;

1905 word1(rv) = Big1;

1906 }

1907 else

1908 word0(rv) += P*Exp_msk1;

1909 }

1910 }

1911 else if (e1 < 0) {

1912 e1 = -e1;

1913 if (i = e1 & 15)

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

1915 if (e1 >>= 4) {

1916 if (e1 >= 1 << n_bigtens)

1917 goto undfl;

1918 #ifdef Avoid_Underflow

1919 if (e1 & Scale_Bit)

1920 scale = 2*P;

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

1922 if (e1 & 1)

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

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

1925 >> Exp_shift)) > 0) {

1926 /* scaled rv is denormal; clear j low bits */

1927 if (j >= 32) {

1928 word1(rv) = 0;

1929 if (j >= 53)

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

1931 else

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

1933 }

1934 else

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

1936 }

1937 #else

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

1939 if (e1 & 1)

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

1941 /* The last multiplication could underflow. */

1942 dval(rv0) = dval(rv);

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

1944 if (!dval(rv)) {

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

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

1947 #endif

1948 if (!dval(rv)) {

1949 undfl:

1950 dval(rv) = 0.;

1951 #ifndef NO_ERRNO

1952 errno = ERANGE;

1953 #endif

1954 if (bd0)

1955 goto retfree;

1956 goto ret;

1957 }

1958 #ifndef Avoid_Underflow

1959 word0(rv) = Tiny0;

1960 word1(rv) = Tiny1;

1961 /* The refinement below will clean

1962 * this approximation up.

1963 */

1964 }

1965 #endif

1966 }

1967 }

1968

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

1970

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

1972

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

1974

1975 for(;;) {

1976 bd = Balloc(bd0->k);

1977 Bcopy(bd, bd0);

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

1979 bs = i2b(1);

1980

1981 if (e >= 0) {

1982 bb2 = bb5 = 0;

1983 bd2 = bd5 = e;

1984 }

1985 else {

1986 bb2 = bb5 = -e;

1987 bd2 = bd5 = 0;

1988 }

1989 if (bbe >= 0)

1990 bb2 += bbe;

1991 else

1992 bd2 -= bbe;

1993 bs2 = bb2;

1994 #ifdef Honor_FLT_ROUNDS

1995 if (Rounding != 1)

1996 bs2++;

1997 #endif

1998 #ifdef Avoid_Underflow

1999 j = bbe - scale;

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

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

2002 j += P - Emin;

2003 else

2004 j = P + 1 - bbbits;

2005 #else /Avoid_Underflow/

2006 #ifdef Sudden_Underflow

2007 #ifdef IBM

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

2009 #else

2010 j = P + 1 - bbbits;

2011 #endif

2012 #else /Sudden_Underflow/

2013 j = bbe;

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

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

2016 j += P - Emin;

2017 else

2018 j = P + 1 - bbbits;

2019 #endif /Sudden_Underflow/

2020 #endif /Avoid_Underflow/

2021 bb2 += j;

2022 bd2 += j;

2023 #ifdef Avoid_Underflow

2024 bd2 += scale;

2025 #endif

2026 i = bb2 < bd2 ? bb2 : bd2;

2027 if (i > bs2)

2028 i = bs2;

2029 if (i > 0) {

2030 bb2 -= i;

2031 bd2 -= i;

2032 bs2 -= i;

2033 }

2034 if (bb5 > 0) {

2035 bs = pow5mult(bs, bb5);

2036 bb1 = mult(bs, bb);

2037 Bfree(bb);

2038 bb = bb1;

2039 }

2040 if (bb2 > 0)

2041 bb = lshift(bb, bb2);

2042 if (bd5 > 0)

2043 bd = pow5mult(bd, bd5);

2044 if (bd2 > 0)

2045 bd = lshift(bd, bd2);

2046 if (bs2 > 0)

2047 bs = lshift(bs, bs2);

2048 delta = diff(bb, bd);

2049 dsign = delta->sign;

2050 delta->sign = 0;

2051 i = cmp(delta, bs);

2052 #ifdef Honor_FLT_ROUNDS

2053 if (Rounding != 1) {

2054 if (i < 0) {

2055 /* Error is less than an ulp */

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

2057 /* exact */

2058 #ifdef SET_INEXACT

2059 inexact = 0;

2060 #endif

2061 break;

2062 }

2063 if (Rounding) {

2064 if (dsign) {

2065 adj = 1.;

2066 goto apply_adj;

2067 }

2068 }

2069 else if (!dsign) {

2070 adj = -1.;

2071 if (!word1(rv)

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

2073 y = word0(rv) & Exp_mask;

2074 #ifdef Avoid_Underflow

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

2076 #else

2077 if (y)

2078 #endif

2079 {

2080 delta = lshift(delta,Log2P);

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

2082 adj = -0.5;

2083 }

2084 }

2085 apply_adj:

2086 #ifdef Avoid_Underflow

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

2088 <= 2PExp_msk1)

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

2090 #else

2091 #ifdef Sudden_Underflow

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

2093 P*Exp_msk1) {

2094 word0(rv) += P*Exp_msk1;

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

2096 word0(rv) -= P*Exp_msk1;

2097 }

2098 else

2099 #endif /Sudden_Underflow/

2100 #endif /Avoid_Underflow/

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

2102 }

2103 break;

2104 }

2105 adj = ratio(delta, bs);

2106 if (adj < 1.)

2107 adj = 1.;

2108 if (adj <= 0x7ffffffe) {

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

2110 y = adj;

2111 if (y != adj) {

2112 if (!((Rounding>>1) ^ dsign))

2113 y++;

2114 adj = y;

2115 }

2116 }

2117 #ifdef Avoid_Underflow

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

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

2120 #else

2121 #ifdef Sudden_Underflow

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

2123 word0(rv) += P*Exp_msk1;

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

2125 if (dsign)

2126 dval(rv) += adj;

2127 else

2128 dval(rv) -= adj;

2129 word0(rv) -= P*Exp_msk1;

2130 goto cont;

2131 }

2132 #endif /Sudden_Underflow/

2133 #endif /Avoid_Underflow/

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

2135 if (dsign) {

2136 if (word0(rv) == Big0 && word1(rv) == Big1)

2137 goto ovfl;

2138 dval(rv) += adj;

2139 }

2140 else

2141 dval(rv) -= adj;

2142 goto cont;

2143 }

2144 #endif /Honor_FLT_ROUNDS/

2145

2146 if (i < 0) {

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

2148 * special case of mantissa a power of two.

2149 */

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

2151 #ifdef IEEE_Arith

2152 #ifdef Avoid_Underflow

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

2154 #else

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

2156 #endif

2157 #endif

2158 ) {

2159 #ifdef SET_INEXACT

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

2161 inexact = 0;

2162 #endif

2163 break;

2164 }

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

2166 /* exact result */

2167 #ifdef SET_INEXACT

2168 inexact = 0;

2169 #endif

2170 break;

2171 }

2172 delta = lshift(delta,Log2P);

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

2174 goto drop_down;

2175 break;

2176 }

2177 if (i == 0) {

2178 /* exactly half-way between */

2179 if (dsign) {

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

2181 && word1(rv) == (

2182 #ifdef Avoid_Underflow

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

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

2185 #endif

2186 0xffffffff)) {

2187 /boundary case -- increment exponent/

2188 word0(rv) = (word0(rv) & Exp_mask)

2189 + Exp_msk1

2190 #ifdef IBM

2191 \| Exp_msk1 >> 4

2192 #endif

2193 ;

2194 word1(rv) = 0;

2195 #ifdef Avoid_Underflow

2196 dsign = 0;

2197 #endif

2198 break;

2199 }

2200 }

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

2202 drop_down:

2203 /* boundary case -- decrement exponent */

2204 #ifdef Sudden_Underflow /{{/

2205 L = word0(rv) & Exp_mask;

2206 #ifdef IBM

2207 if (L < Exp_msk1)

2208 #else

2209 #ifdef Avoid_Underflow

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

2211 #else

2212 if (L <= Exp_msk1)

2213 #endif /Avoid_Underflow/

2214 #endif /IBM/

2215 goto undfl;

2216 L -= Exp_msk1;

2217 #else /Sudden_Underflow}{/

2218 #ifdef Avoid_Underflow

2219 if (scale) {

2220 L = word0(rv) & Exp_mask;

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

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

2223 /* round even ==> */

2224 /* accept rv */

2225 break;

2226 /* rv = smallest denormal */

2227 goto undfl;

2228 }

2229 }

2230 #endif /Avoid_Underflow/

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

2232 #endif /Sudden_Underflow}}/

2233 word0(rv) = L \| Bndry_mask1;

2234 word1(rv) = 0xffffffff;

2235 #ifdef IBM

2236 goto cont;

2237 #else

2238 break;

2239 #endif

2240 }

2241 #ifndef ROUND_BIASED

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

2243 break;

2244 #endif

2245 if (dsign)

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

2247 #ifndef ROUND_BIASED

2248 else {

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

2250 #ifndef Sudden_Underflow

2251 if (!dval(rv))

2252 goto undfl;

2253 #endif

2254 }

2255 #ifdef Avoid_Underflow

2256 dsign = 1 - dsign;

2257 #endif

2258 #endif

2259 break;

2260 }

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

2262 if (dsign)

2263 aadj = aadj1 = 1.;

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

2265 #ifndef Sudden_Underflow

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

2267 goto undfl;

2268 #endif

2269 aadj = 1.;

2270 aadj1 = -1.;

2271 }

2272 else {

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

2274 /* rounded down... */

2275

2276 if (aadj < 2./FLT_RADIX)

2277 aadj = 1./FLT_RADIX;

2278 else

2279 aadj *= 0.5;

2280 aadj1 = -aadj;

2281 }

2282 }

2283 else {

2284 aadj *= 0.5;

2285 aadj1 = dsign ? aadj : -aadj;

2286 #ifdef Check_FLT_ROUNDS

2287 switch(Rounding) {

2288 case 2: /* towards +infinity */

2289 aadj1 -= 0.5;

2290 break;

2291 case 0: /* towards 0 */

2292 case 3: /* towards -infinity */

2293 aadj1 += 0.5;

2294 }

2295 #else

2296 if (Flt_Rounds == 0)

2297 aadj1 += 0.5;

2298 #endif /Check_FLT_ROUNDS/

2299 }

2300 y = word0(rv) & Exp_mask;

2301

2302 /* Check for overflow */

2303

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

2305 dval(rv0) = dval(rv);

2306 word0(rv) -= P*Exp_msk1;

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

2308 dval(rv) += adj;

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

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

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

2312 goto ovfl;

2313 word0(rv) = Big0;

2314 word1(rv) = Big1;

2315 goto cont;

2316 }

2317 else

2318 word0(rv) += P*Exp_msk1;

2319 }

2320 else {

2321 #ifdef Avoid_Underflow

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

2323 if (aadj <= 0x7fffffff) {

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

2325 z = 1;

2326 aadj = z;

2327 aadj1 = dsign ? aadj : -aadj;

2328 }

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

2330 }

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

2332 dval(rv) += adj;

2333 #else

2334 #ifdef Sudden_Underflow

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

2336 dval(rv0) = dval(rv);

2337 word0(rv) += P*Exp_msk1;

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

2339 dval(rv) += adj;

2340 #ifdef IBM

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

2342 #else

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

2344 #endif

2345 {

2346 if (word0(rv0) == Tiny0

2347 && word1(rv0) == Tiny1)

2348 goto undfl;

2349 word0(rv) = Tiny0;

2350 word1(rv) = Tiny1;

2351 goto cont;

2352 }

2353 else

2354 word0(rv) -= P*Exp_msk1;

2355 }

2356 else {

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

2358 dval(rv) += adj;

2359 }

2360 #else /Sudden_Underflow/

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

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

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

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

2365 * trouble from bits lost to denormalization;

2366 * example: 1.2e-307 .

2367 */

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

2369 aadj1 = (double)(int)(aadj + 0.5);

2370 if (!dsign)

2371 aadj1 = -aadj1;

2372 }

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

2374 dval(rv) += adj;

2375 #endif /Sudden_Underflow/

2376 #endif /Avoid_Underflow/

2377 }

2378 z = word0(rv) & Exp_mask;

2379 #ifndef SET_INEXACT

2380 #ifdef Avoid_Underflow

2381 if (!scale)

2382 #endif

2383 if (y == z) {

2384 /* Can we stop now? */

2385 L = (Long)aadj;

2386 aadj -= L;

2387 /* The tolerances below are conservative. */

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

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

2390 break;

2391 }

2392 else if (aadj < .4999999/FLT_RADIX)

2393 break;

2394 }

2395 #endif

2396 cont:

2397 Bfree(bb);

2398 Bfree(bd);

2399 Bfree(bs);

2400 Bfree(delta);

2401 }

2402 #ifdef SET_INEXACT

2403 if (inexact) {

2404 if (!oldinexact) {

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

2406 word1(rv0) = 0;

2407 dval(rv0) += 1.;

2408 }

2409 }

2410 else if (!oldinexact)

2411 clear_inexact();

2412 #endif

2413 #ifdef Avoid_Underflow

2414 if (scale) {

2415 word0(rv0) = Exp_1 - 2PExp_msk1;

2416 word1(rv0) = 0;

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

2418 #ifndef NO_ERRNO

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

2420 #ifdef IEEE_Arith

2421 if (!(word0(rv) & Exp_mask))

2422 #else

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

2424 #endif

2425 errno = ERANGE;

2426 #endif

2427 }

2428 #endif /* Avoid_Underflow */

2429 #ifdef SET_INEXACT

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

2431 /* set underflow bit */

2432 dval(rv0) = 1e-300;

2433 dval(rv0) *= dval(rv0);

2434 }

2435 #endif

2436 retfree:

2437 Bfree(bb);

2438 Bfree(bd);

2439 Bfree(bs);

2440 Bfree(bd0);

2441 Bfree(delta);

2442 ret:

2443 if (se)

2444 se = (char )s;

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

2446 }

2447

2448 static int

2449 quorem

2450 #ifdef KR_headers

2451 (b, S) Bigint b, S;

2452 #else

2453 (Bigint b, Bigint S)

2454 #endif

2455 {

2456 int n;

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

2458 #ifdef ULLong

2459 ULLong borrow, carry, y, ys;

2460 #else

2461 ULong borrow, carry, y, ys;

2462 #ifdef Pack_32

2463 ULong si, z, zs;

2464 #endif

2465 #endif

2466

2467 n = S->wds;

2468 #ifdef DEBUG

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

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

2471 #endif

2472 if (b->wds < n)

2473 return 0;

2474 sx = S->x;

2475 sxe = sx + --n;

2476 bx = b->x;

2477 bxe = bx + n;

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

2479 #ifdef DEBUG

2480 /debug/ if (q > 9)

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

2482 #endif

2483 if (q) {

2484 borrow = 0;

2485 carry = 0;

2486 do {

2487 #ifdef ULLong

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

2489 carry = ys >> 32;

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

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

2492 *bx++ = y & FFFFFFFF;

2493 #else

2494 #ifdef Pack_32

2495 si = *sx++;

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

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

2498 carry = zs >> 16;

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

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

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

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

2503 Storeinc(bx, z, y);

2504 #else

2505 ys = sx++ q + carry;

2506 carry = ys >> 16;

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

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

2509 *bx++ = y & 0xffff;

2510 #endif

2511 #endif

2512 }

2513 while(sx <= sxe);

2514 if (!*bxe) {

2515 bx = b->x;

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

2517 --n;

2518 b->wds = n;

2519 }

2520 }

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

2522 q++;

2523 borrow = 0;

2524 carry = 0;

2525 bx = b->x;

2526 sx = S->x;

2527 do {

2528 #ifdef ULLong

2529 ys = *sx++ + carry;

2530 carry = ys >> 32;

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

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

2533 *bx++ = y & FFFFFFFF;

2534 #else

2535 #ifdef Pack_32

2536 si = *sx++;

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

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

2539 carry = zs >> 16;

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

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

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

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

2544 Storeinc(bx, z, y);

2545 #else

2546 ys = *sx++ + carry;

2547 carry = ys >> 16;

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

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

2550 *bx++ = y & 0xffff;

2551 #endif

2552 #endif

2553 }

2554 while(sx <= sxe);

2555 bx = b->x;

2556 bxe = bx + n;

2557 if (!*bxe) {

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

2559 --n;

2560 b->wds = n;

2561 }

2562 }

2563 return q;

2564 }

2565

2566 #ifndef MULTIPLE_THREADS

2567 static char *dtoa_result;

2568 #endif

2569

2570 static char *

2571 #ifdef KR_headers

2572 rv_alloc(i) int i;

2573 #else

2574 rv_alloc(int i)

2575 #endif

2576 {

2577 int j, k, *r;

2578

2579 j = sizeof(ULong);

2580 for(k = 0;

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

2582 j <<= 1)

2583 k++;

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

2585 *r = k;

2586 return

2587 #ifndef MULTIPLE_THREADS

2588 dtoa_result =

2589 #endif

2590 (char *)(r+1);

2591 }

2592

2593 static char *

2594 #ifdef KR_headers

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

2596 #else

2597 nrv_alloc(char s, char *rve, int n)

2598 #endif

2599 {

2600 char rv, t;

2601

2602 t = rv = rv_alloc(n);

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

2604 if (rve)

2605 *rve = t;

2606 return rv;

2607 }

2608

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

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

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

2612 * when MULTIPLE_THREADS is not defined.

2613 */

2614

2615 void

2616 #ifdef KR_headers

2617 freedtoa(s) char *s;

2618 #else

2619 freedtoa(char *s)

2620 #endif

2621 {

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

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

2624 Bfree(b);

2625 #ifndef MULTIPLE_THREADS

2626 if (s == dtoa_result)

2627 dtoa_result = 0;

2628 #endif

2629 }

2630

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

2632 *

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

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

2635 *

2636 * Modifications:

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

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

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

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

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

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

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

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

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

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

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

2648 * as would satisfaction of the stopping test with strict

2649 * inequality.

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

2651 * quantities.

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

2653 * we use floating-point arithmetic rather than resorting

2654 * to multiple-precision integers.

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

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

2657 * multiple-precision integer arithmetic only if we cannot

2658 * guarantee that the floating-point calculation has given

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

2660 * "uniformly" distributed input, the probability is

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

2662 * calculation.

2663 */

2664

2665 char *

2666 dtoa

2667 #ifdef KR_headers

2668 (d, mode, ndigits, decpt, sign, rve)

2669 double d; int mode, ndigits, decpt, sign; char **rve;

2670 #else

2671 (double d, int mode, int ndigits, int decpt, int sign, char **rve)

2672 #endif

2673 {

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

2675 of ecvt and fcvt; trailing zeros are suppressed from

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

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

2678 then *decpt is set to 9999.

2679

2680 mode:

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

2682 and rounded to nearest.

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

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

2685 1e23 whereas mode 1 gives 9.999999999999999e22.

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

2687 return value similar to that of ecvt, except

2688 that trailing zeros are suppressed.

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

2690 gives a return value similar to that from fcvt,

2691 except that trailing zeros are suppressed, and

2692 ndigits can be negative.

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

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

2695 possibly return a shorter string that rounds to d.

2696 With IEEE arithmetic and compilation with

2697 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same

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

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

2700 fast floating-point estimate (if applicable).

2701

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

2703

2704 Sufficient space is allocated to the return value

2705 to hold the suppressed trailing zeros.

2706 */

2707

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

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

2710 spec_case, try_quick;

2711 Long L;

2712 #ifndef Sudden_Underflow

2713 int denorm;

2714 ULong x;

2715 #endif

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

2717 double d2, ds, eps;

2718 char s, s0;

2719 #ifdef SET_INEXACT

2720 int inexact, oldinexact;

2721 #endif

2722 #ifdef Honor_FLT_ROUNDS /{/

2723 int Rounding;

2724 #ifdef Trust_FLT_ROUNDS /{{ only define this if FLT_ROUNDS really works! /

2725 Rounding = Flt_Rounds;

2726 #else /}{/

2727 Rounding = 1;

2728 switch(fegetround()) {

2729 case FE_TOWARDZERO: Rounding = 0; break;

2730 case FE_UPWARD: Rounding = 2; break;

2731 case FE_DOWNWARD: Rounding = 3;

2732 }

2733 #endif /}}/

2734 #endif /}/

2735

2736 #ifndef MULTIPLE_THREADS

2737 if (dtoa_result) {

2738 freedtoa(dtoa_result);

2739 dtoa_result = 0;

2740 }

2741 #endif

2742

2743 if (word0(d) & Sign_bit) {

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

2745 *sign = 1;

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

2747 }

2748 else

2749 *sign = 0;

2750

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

2752 #ifdef IEEE_Arith

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

2754 #else

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

2756 #endif

2757 {

2758 /* Infinity or NaN */

2759 *decpt = 9999;

2760 #ifdef IEEE_Arith

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

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

2763 #endif

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

2765 }

2766 #endif

2767 #ifdef IBM

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

2769 #endif

2770 if (!dval(d)) {

2771 *decpt = 1;

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

2773 }

2774

2775 #ifdef SET_INEXACT

2776 try_quick = oldinexact = get_inexact();

2777 inexact = 1;

2778 #endif

2779 #ifdef Honor_FLT_ROUNDS

2780 if (Rounding >= 2) {

2781 if (*sign)

2782 Rounding = Rounding == 2 ? 0 : 2;

2783 else

2784 if (Rounding != 2)

2785 Rounding = 0;

2786 }

2787 #endif

2788

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

2790 #ifdef Sudden_Underflow

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

2792 #else

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

2794 #endif

2795 dval(d2) = dval(d);

2796 word0(d2) &= Frac_mask1;

2797 word0(d2) \|= Exp_11;

2798 #ifdef IBM

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

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

2801 #endif

2802

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

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

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

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

2807 *

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

2809 *

2810 * k = (i - Bias)*0.301029995663981

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

2812 *

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

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

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

2816 * to compensate for any error in the multiplication of

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

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

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

2820 * Hence we adjust the constant term to 0.1760912590558.

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

2822 * but this is probably not worthwhile.)

2823 */

2824

2825 i -= Bias;

2826 #ifdef IBM

2827 i <<= 2;

2828 i += j;

2829 #endif

2830 #ifndef Sudden_Underflow

2831 denorm = 0;

2832 }

2833 else {

2834 /* d is denormalized */

2835

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

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

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

2839 dval(d2) = x;

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

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

2842 denorm = 1;

2843 }

2844 #endif

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

2846 k = (int)ds;

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

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

2849 k_check = 1;

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

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

2852 k--;

2853 k_check = 0;

2854 }

2855 j = bbits - i - 1;

2856 if (j >= 0) {

2857 b2 = 0;

2858 s2 = j;

2859 }

2860 else {

2861 b2 = -j;

2862 s2 = 0;

2863 }

2864 if (k >= 0) {

2865 b5 = 0;

2866 s5 = k;

2867 s2 += k;

2868 }

2869 else {

2870 b2 -= k;

2871 b5 = -k;

2872 s5 = 0;

2873 }

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

2875 mode = 0;

2876

2877 #ifndef SET_INEXACT

2878 #ifdef Check_FLT_ROUNDS

2879 try_quick = Rounding == 1;

2880 #else

2881 try_quick = 1;

2882 #endif

2883 #endif /SET_INEXACT/

2884

2885 if (mode > 5) {

2886 mode -= 4;

2887 try_quick = 0;

2888 }

2889 leftright = 1;

2890 switch(mode) {

2891 case 0:

2892 case 1:

2893 ilim = ilim1 = -1;

2894 i = 18;

2895 ndigits = 0;

2896 break;

2897 case 2:

2898 leftright = 0;

2899 /* no break */

2900 case 4:

2901 if (ndigits <= 0)

2902 ndigits = 1;

2903 ilim = ilim1 = i = ndigits;

2904 break;

2905 case 3:

2906 leftright = 0;

2907 /* no break */

2908 case 5:

2909 i = ndigits + k + 1;

2910 ilim = i;

2911 ilim1 = i - 1;

2912 if (i <= 0)

2913 i = 1;

2914 }

2915 s = s0 = rv_alloc(i);

2916

2917 #ifdef Honor_FLT_ROUNDS

2918 if (mode > 1 && Rounding != 1)

2919 leftright = 0;

2920 #endif

2921

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

2923

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

2925

2926 i = 0;

2927 dval(d2) = dval(d);

2928 k0 = k;

2929 ilim0 = ilim;

2930 ieps = 2; /* conservative */

2931 if (k > 0) {

2932 ds = tens[k&0xf];

2933 j = k >> 4;

2934 if (j & Bletch) {

2935 /* prevent overflows */

2936 j &= Bletch - 1;

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

2938 ieps++;

2939 }

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

2941 if (j & 1) {

2942 ieps++;

2943 ds *= bigtens[i];

2944 }

2945 dval(d) /= ds;

2946 }

2947 else if (j1 = -k) {

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

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

2950 if (j & 1) {

2951 ieps++;

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

2953 }

2954 }

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

2956 if (ilim1 <= 0)

2957 goto fast_failed;

2958 ilim = ilim1;

2959 k--;

2960 dval(d) *= 10.;

2961 ieps++;

2962 }

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

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

2965 if (ilim == 0) {

2966 S = mhi = 0;

2967 dval(d) -= 5.;

2968 if (dval(d) > dval(eps))

2969 goto one_digit;

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

2971 goto no_digits;

2972 goto fast_failed;

2973 }

2974 #ifndef No_leftright

2975 if (leftright) {

2976 /* Use Steele & White method of only

2977 * generating digits needed.

2978 */

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

2980 for(i = 0;;) {

2981 L = dval(d);

2982 dval(d) -= L;

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

2984 if (dval(d) < dval(eps))

2985 goto ret1;

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

2987 goto bump_up;

2988 if (++i >= ilim)

2989 break;

2990 dval(eps) *= 10.;

2991 dval(d) *= 10.;

2992 }

2993 }

2994 else {

2995 #endif

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

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

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

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

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

3001 ilim = i;

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

3003 if (i == ilim) {

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

3005 goto bump_up;

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

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

3008 s++;

3009 goto ret1;

3010 }

3011 break;

3012 }

3013 }

3014 #ifndef No_leftright

3015 }

3016 #endif

3017 fast_failed:

3018 s = s0;

3019 dval(d) = dval(d2);

3020 k = k0;

3021 ilim = ilim0;

3022 }

3023

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

3025

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

3027 /* Yes. */

3028 ds = tens[k];

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

3030 S = mhi = 0;

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

3032 goto no_digits;

3033 goto one_digit;

3034 }

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

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

3037 dval(d) -= L*ds;

3038 #ifdef Check_FLT_ROUNDS

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

3040 if (dval(d) < 0) {

3041 L--;

3042 dval(d) += ds;

3043 }

3044 #endif

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

3046 if (!dval(d)) {

3047 #ifdef SET_INEXACT

3048 inexact = 0;

3049 #endif

3050 break;

3051 }

3052 if (i == ilim) {

3053 #ifdef Honor_FLT_ROUNDS

3054 if (mode > 1)

3055 switch(Rounding) {

3056 case 0: goto ret1;

3057 case 2: goto bump_up;

3058 }

3059 #endif

3060 dval(d) += dval(d);

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

3062 bump_up:

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

3064 if (s == s0) {

3065 k++;

3066 *s = '0';

3067 break;

3068 }

3069 ++*s++;

3070 }

3071 break;

3072 }

3073 }

3074 goto ret1;

3075 }

3076

3077 m2 = b2;

3078 m5 = b5;

3079 mhi = mlo = 0;

3080 if (leftright) {

3081 i =

3082 #ifndef Sudden_Underflow

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

3084 #endif

3085 #ifdef IBM

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

3087 #else

3088 1 + P - bbits;

3089 #endif

3090 b2 += i;

3091 s2 += i;

3092 mhi = i2b(1);

3093 }

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

3095 i = m2 < s2 ? m2 : s2;

3096 b2 -= i;

3097 m2 -= i;

3098 s2 -= i;

3099 }

3100 if (b5 > 0) {

3101 if (leftright) {

3102 if (m5 > 0) {

3103 mhi = pow5mult(mhi, m5);

3104 b1 = mult(mhi, b);

3105 Bfree(b);

3106 b = b1;

3107 }

3108 if (j = b5 - m5)

3109 b = pow5mult(b, j);

3110 }

3111 else

3112 b = pow5mult(b, b5);

3113 }

3114 S = i2b(1);

3115 if (s5 > 0)

3116 S = pow5mult(S, s5);

3117

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

3119

3120 spec_case = 0;

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

3122 #ifdef Honor_FLT_ROUNDS

3123 && Rounding == 1

3124 #endif

3125 ) {

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

3127 #ifndef Sudden_Underflow

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

3129 #endif

3130 ) {

3131 /* The special case */

3132 b2 += Log2P;

3133 s2 += Log2P;

3134 spec_case = 1;

3135 }

3136 }

3137

3138 /* Arrange for convenient computation of quotients:

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

3140 *

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

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

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

3144 */

3145 #ifdef Pack_32

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

3147 i = 32 - i;

3148 #else

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

3150 i = 16 - i;

3151 #endif

3152 if (i > 4) {

3153 i -= 4;

3154 b2 += i;

3155 m2 += i;

3156 s2 += i;

3157 }

3158 else if (i < 4) {

3159 i += 28;

3160 b2 += i;

3161 m2 += i;

3162 s2 += i;

3163 }

3164 if (b2 > 0)

3165 b = lshift(b, b2);

3166 if (s2 > 0)

3167 S = lshift(S, s2);

3168 if (k_check) {

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

3170 k--;

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

3172 if (leftright)

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

3174 ilim = ilim1;

3175 }

3176 }

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

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

3179 /* no digits, fcvt style */

3180 no_digits:

3181 k = -1 - ndigits;

3182 goto ret;

3183 }

3184 one_digit:

3185 *s++ = '1';

3186 k++;

3187 goto ret;

3188 }

3189 if (leftright) {

3190 if (m2 > 0)

3191 mhi = lshift(mhi, m2);

3192

3193 /* Compute mlo -- check for special case

3194 * that d is a normalized power of 2.

3195 */

3196

3197 mlo = mhi;

3198 if (spec_case) {

3199 mhi = Balloc(mhi->k);

3200 Bcopy(mhi, mlo);

3201 mhi = lshift(mhi, Log2P);

3202 }

3203

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

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

3206 /* Do we yet have the shortest decimal string

3207 * that will round to d?

3208 */

3209 j = cmp(b, mlo);

3210 delta = diff(S, mhi);

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

3212 Bfree(delta);

3213 #ifndef ROUND_BIASED

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

3215 #ifdef Honor_FLT_ROUNDS

3216 && Rounding >= 1

3217 #endif

3218 ) {

3219 if (dig == '9')

3220 goto round_9_up;

3221 if (j > 0)

3222 dig++;

3223 #ifdef SET_INEXACT

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

3225 inexact = 0;

3226 #endif

3227 *s++ = dig;

3228 goto ret;

3229 }

3230 #endif

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

3232 #ifndef ROUND_BIASED

3233 && !(word1(d) & 1)

3234 #endif

3235 ) {

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

3237 #ifdef SET_INEXACT

3238 inexact = 0;

3239 #endif

3240 goto accept_dig;

3241 }

3242 #ifdef Honor_FLT_ROUNDS

3243 if (mode > 1)

3244 switch(Rounding) {

3245 case 0: goto accept_dig;

3246 case 2: goto keep_dig;

3247 }

3248 #endif /Honor_FLT_ROUNDS/

3249 if (j1 > 0) {

3250 b = lshift(b, 1);

3251 j1 = cmp(b, S);

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

3253 && dig++ == '9')

3254 goto round_9_up;

3255 }

3256 accept_dig:

3257 *s++ = dig;

3258 goto ret;

3259 }

3260 if (j1 > 0) {

3261 #ifdef Honor_FLT_ROUNDS

3262 if (!Rounding)

3263 goto accept_dig;

3264 #endif

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

3266 round_9_up:

3267 *s++ = '9';

3268 goto roundoff;

3269 }

3270 *s++ = dig + 1;

3271 goto ret;

3272 }

3273 #ifdef Honor_FLT_ROUNDS

3274 keep_dig:

3275 #endif

3276 *s++ = dig;

3277 if (i == ilim)

3278 break;

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

3280 if (mlo == mhi)

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

3282 else {

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

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

3285 }

3286 }

3287 }

3288 else

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

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

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

3292 #ifdef SET_INEXACT

3293 inexact = 0;

3294 #endif

3295 goto ret;

3296 }

3297 if (i >= ilim)

3298 break;

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

3300 }

3301

3302 /* Round off last digit */

3303

3304 #ifdef Honor_FLT_ROUNDS

3305 switch(Rounding) {

3306 case 0: goto trimzeros;

3307 case 2: goto roundoff;

3308 }

3309 #endif

3310 b = lshift(b, 1);

3311 j = cmp(b, S);

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

3313 roundoff:

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

3315 if (s == s0) {

3316 k++;

3317 *s++ = '1';

3318 goto ret;

3319 }

3320 ++*s++;

3321 }

3322 else {

3323 trimzeros:

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

3325 s++;

3326 }

3327 ret:

3328 Bfree(S);

3329 if (mhi) {

3330 if (mlo && mlo != mhi)

3331 Bfree(mlo);

3332 Bfree(mhi);

3333 }

3334 ret1:

3335 #ifdef SET_INEXACT

3336 if (inexact) {

3337 if (!oldinexact) {

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

3339 word1(d) = 0;

3340 dval(d) += 1.;

3341 }

3342 }

3343 else if (!oldinexact)

3344 clear_inexact();

3345 #endif

3346 Bfree(b);

3347 *s = 0;

3348 *decpt = k + 1;

3349 if (rve)

3350 *rve = s;

3351 return s0;

3352 }

3353

3354 } // namespace dmg_fp

OLD	NEW

« no previous file with comments | « third_party/dmg_fp/dmg_fp.xcodeproj/project.pbxproj ('k') | third_party/dmg_fp/g_fmt.cc » ('j') | no next file with comments »