base/third_party/dmg_fp/dtoa.cc - Issue 1647803004: Move base to DEPS

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

Issue 1647803004: Move base to DEPS (Closed) Base URL: git@github.com:domokit/mojo.git@master

Patch Set: Created 4 years, 10 months 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. Similarly, if you

107 * want something other than the system's free() to be called to

108 * recycle memory acquired from MALLOC, #define FREE to be the

109 * name of the alternate routine. (FREE or free is only called in

110 * pathological cases, e.g., in a dtoa call after a dtoa return in

111 * mode 3 with thousands of digits requested.)

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

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

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

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

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

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

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

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

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

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

122 * pointers, PRIVATE_MEM >= 7112 appears adequate.

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

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

125 * when INFNAN_CHECK is #defined, strtod checks

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

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

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

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

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

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

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

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

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

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

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

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

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

139 * and NAN_WORD1 are used instead.

140 * #define MULTIPLE_THREADS if the system offers preemptively scheduled

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

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

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

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

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

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

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

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

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

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

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

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

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

154 * Sudden_Underflow.

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 * #define NO_HEX_FP to omit recognition of hexadecimal floating-point

173 * values by strtod.

174 * #define NO_STRTOD_BIGCOMP (on IEEE-arithmetic systems only for now)

175 * to disable logic for "fast" testing of very long input strings

176 * to strtod. This testing proceeds by initially truncating the

177 * input string, then if necessary comparing the whole string with

178 * a decimal expansion to decide close cases. This logic is only

179 * used for input more than STRTOD_DIGLIM digits long (default 40).

180 */

181

182 #define IEEE_8087

183 #define NO_HEX_FP

184

185 #ifndef Long

186 #if __LP64__

187 #define Long int

188 #else

189 #define Long long

190 #endif

191 #endif

192 #ifndef ULong

193 typedef unsigned Long ULong;

194 #endif

195

196 #ifdef DEBUG

197 #include "stdio.h"

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

199 #endif

200

201 #include "stdlib.h"

202 #include "string.h"

203

204 #ifdef USE_LOCALE

205 #include "locale.h"

206 #endif

207

208 #ifdef Honor_FLT_ROUNDS

209 #ifndef Trust_FLT_ROUNDS

210 #include <fenv.h>

211 #endif

212 #endif

213

214 #ifdef MALLOC

215 #ifdef KR_headers

216 extern char *MALLOC();

217 #else

218 extern void *MALLOC(size_t);

219 #endif

220 #else

221 #define MALLOC malloc

222 #endif

223

224 #ifndef Omit_Private_Memory

225 #ifndef PRIVATE_MEM

226 #define PRIVATE_MEM 2304

227 #endif

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

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

230 #endif

231

232 #undef IEEE_Arith

233 #undef Avoid_Underflow

234 #ifdef IEEE_MC68k

235 #define IEEE_Arith

236 #endif

237 #ifdef IEEE_8087

238 #define IEEE_Arith

239 #endif

240

241 #ifdef IEEE_Arith

242 #ifndef NO_INFNAN_CHECK

243 #undef INFNAN_CHECK

244 #define INFNAN_CHECK

245 #endif

246 #else

247 #undef INFNAN_CHECK

248 #define NO_STRTOD_BIGCOMP

249 #endif

250

251 #include "errno.h"

252

253 #ifdef Bad_float_h

254

255 #ifdef IEEE_Arith

256 #define DBL_DIG 15

257 #define DBL_MAX_10_EXP 308

258 #define DBL_MAX_EXP 1024

259 #define FLT_RADIX 2

260 #endif /IEEE_Arith/

261

262 #ifdef IBM

263 #define DBL_DIG 16

264 #define DBL_MAX_10_EXP 75

265 #define DBL_MAX_EXP 63

266 #define FLT_RADIX 16

267 #define DBL_MAX 7.2370055773322621e+75

268 #endif

269

270 #ifdef VAX

271 #define DBL_DIG 16

272 #define DBL_MAX_10_EXP 38

273 #define DBL_MAX_EXP 127

274 #define FLT_RADIX 2

275 #define DBL_MAX 1.7014118346046923e+38

276 #endif

277

278 #ifndef LONG_MAX

279 #define LONG_MAX 2147483647

280 #endif

281

282 #else /* ifndef Bad_float_h */

283 #include "float.h"

284 #endif /* Bad_float_h */

285

286 #ifndef __MATH_H__

287 #include "math.h"

288 #endif

289

290 namespace dmg_fp {

291

292 #ifndef CONST

293 #ifdef KR_headers

294 #define CONST /* blank */

295 #else

296 #define CONST const

297 #endif

298 #endif

299

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

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

302 #endif

303

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

305

306 #ifdef IEEE_8087

307 #define word0(x) (x)->L[1]

308 #define word1(x) (x)->L[0]

309 #else

310 #define word0(x) (x)->L[0]

311 #define word1(x) (x)->L[1]

312 #endif

313 #define dval(x) (x)->d

314

315 #ifndef STRTOD_DIGLIM

316 #define STRTOD_DIGLIM 40

317 #endif

318

319 #ifdef DIGLIM_DEBUG

320 extern int strtod_diglim;

321 #else

322 #define strtod_diglim STRTOD_DIGLIM

323 #endif

324

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

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

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

328 */

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

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

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

332 #else

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

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

335 #endif

336

337 /* #define P DBL_MANT_DIG */

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

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

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

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

342

343 #ifdef IEEE_Arith

344 #define Exp_shift 20

345 #define Exp_shift1 20

346 #define Exp_msk1 0x100000

347 #define Exp_msk11 0x100000

348 #define Exp_mask 0x7ff00000

349 #define P 53

350 #define Nbits 53

351 #define Bias 1023

352 #define Emax 1023

353 #define Emin (-1022)

354 #define Exp_1 0x3ff00000

355 #define Exp_11 0x3ff00000

356 #define Ebits 11

357 #define Frac_mask 0xfffff

358 #define Frac_mask1 0xfffff

359 #define Ten_pmax 22

360 #define Bletch 0x10

361 #define Bndry_mask 0xfffff

362 #define Bndry_mask1 0xfffff

363 #define LSB 1

364 #define Sign_bit 0x80000000

365 #define Log2P 1

366 #define Tiny0 0

367 #define Tiny1 1

368 #define Quick_max 14

369 #define Int_max 14

370 #ifndef NO_IEEE_Scale

371 #define Avoid_Underflow

372 #ifdef Flush_Denorm /* debugging option */

373 #undef Sudden_Underflow

374 #endif

375 #endif

376

377 #ifndef Flt_Rounds

378 #ifdef FLT_ROUNDS

379 #define Flt_Rounds FLT_ROUNDS

380 #else

381 #define Flt_Rounds 1

382 #endif

383 #endif /Flt_Rounds/

384

385 #ifdef Honor_FLT_ROUNDS

386 #undef Check_FLT_ROUNDS

387 #define Check_FLT_ROUNDS

388 #else

389 #define Rounding Flt_Rounds

390 #endif

391

392 #else /* ifndef IEEE_Arith */

393 #undef Check_FLT_ROUNDS

394 #undef Honor_FLT_ROUNDS

395 #undef SET_INEXACT

396 #undef Sudden_Underflow

397 #define Sudden_Underflow

398 #ifdef IBM

399 #undef Flt_Rounds

400 #define Flt_Rounds 0

401 #define Exp_shift 24

402 #define Exp_shift1 24

403 #define Exp_msk1 0x1000000

404 #define Exp_msk11 0x1000000

405 #define Exp_mask 0x7f000000

406 #define P 14

407 #define Nbits 56

408 #define Bias 65

409 #define Emax 248

410 #define Emin (-260)

411 #define Exp_1 0x41000000

412 #define Exp_11 0x41000000

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

414 #define Frac_mask 0xffffff

415 #define Frac_mask1 0xffffff

416 #define Bletch 4

417 #define Ten_pmax 22

418 #define Bndry_mask 0xefffff

419 #define Bndry_mask1 0xffffff

420 #define LSB 1

421 #define Sign_bit 0x80000000

422 #define Log2P 4

423 #define Tiny0 0x100000

424 #define Tiny1 0

425 #define Quick_max 14

426 #define Int_max 15

427 #else /* VAX */

428 #undef Flt_Rounds

429 #define Flt_Rounds 1

430 #define Exp_shift 23

431 #define Exp_shift1 7

432 #define Exp_msk1 0x80

433 #define Exp_msk11 0x800000

434 #define Exp_mask 0x7f80

435 #define P 56

436 #define Nbits 56

437 #define Bias 129

438 #define Emax 126

439 #define Emin (-129)

440 #define Exp_1 0x40800000

441 #define Exp_11 0x4080

442 #define Ebits 8

443 #define Frac_mask 0x7fffff

444 #define Frac_mask1 0xffff007f

445 #define Ten_pmax 24

446 #define Bletch 2

447 #define Bndry_mask 0xffff007f

448 #define Bndry_mask1 0xffff007f

449 #define LSB 0x10000

450 #define Sign_bit 0x8000

451 #define Log2P 1

452 #define Tiny0 0x80

453 #define Tiny1 0

454 #define Quick_max 15

455 #define Int_max 15

456 #endif /* IBM, VAX */

457 #endif /* IEEE_Arith */

458

459 #ifndef IEEE_Arith

460 #define ROUND_BIASED

461 #endif

462

463 #ifdef RND_PRODQUOT

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

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

466 #ifdef KR_headers

467 extern double rnd_prod(), rnd_quot();

468 #else

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

470 #endif

471 #else

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

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

474 #endif

475

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

477 #define Big1 0xffffffff

478

479 #ifndef Pack_32

480 #define Pack_32

481 #endif

482

483 typedef struct BCinfo BCinfo;

484 struct

485 BCinfo { int dp0, dp1, dplen, dsign, e0, inexact, nd, nd0, rounding, scale, uflc hk; };

486

487 #ifdef KR_headers

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

489 #else

490 #define FFFFFFFF 0xffffffffUL

491 #endif

492

493 #ifdef NO_LONG_LONG

494 #undef ULLong

495 #ifdef Just_16

496 #undef Pack_32

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

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

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

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

501 */

502 #endif

503 #else /* long long available */

504 #ifndef Llong

505 #define Llong long long

506 #endif

507 #ifndef ULLong

508 #define ULLong unsigned Llong

509 #endif

510 #endif /* NO_LONG_LONG */

511

512 #ifndef MULTIPLE_THREADS

513 #define ACQUIRE_DTOA_LOCK(n) /nothing/

514 #define FREE_DTOA_LOCK(n) /nothing/

515 #endif

516

517 #define Kmax 7

518

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

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

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

522

523 struct

524 Bigint {

525 struct Bigint *next;

526 int k, maxwds, sign, wds;

527 ULong x[1];

528 };

529

530 typedef struct Bigint Bigint;

531

532 static Bigint *freelist[Kmax+1];

533

534 static Bigint *

535 Balloc

536 #ifdef KR_headers

537 (k) int k;

538 #else

539 (int k)

540 #endif

541 {

542 int x;

543 Bigint *rv;

544 #ifndef Omit_Private_Memory

545 unsigned int len;

546 #endif

547

548 ACQUIRE_DTOA_LOCK(0);

549 /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */

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

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

552 rv = freelist[k];

553 freelist[k] = rv->next;

554 }

555 else {

556 x = 1 << k;

557 #ifdef Omit_Private_Memory

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

559 #else

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

561 /sizeof(double);

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

563 rv = (Bigint*)pmem_next;

564 pmem_next += len;

565 }

566 else

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

568 #endif

569 rv->k = k;

570 rv->maxwds = x;

571 }

572 FREE_DTOA_LOCK(0);

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

574 return rv;

575 }

576

577 static void

578 Bfree

579 #ifdef KR_headers

580 (v) Bigint *v;

581 #else

582 (Bigint *v)

583 #endif

584 {

585 if (v) {

586 if (v->k > Kmax)

587 #ifdef FREE

588 FREE((void*)v);

589 #else

590 free((void*)v);

591 #endif

592 else {

593 ACQUIRE_DTOA_LOCK(0);

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

595 freelist[v->k] = v;

596 FREE_DTOA_LOCK(0);

597 }

598 }

599 }

600

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

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

603

604 static Bigint *

605 multadd

606 #ifdef KR_headers

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

608 #else

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

610 #endif

611 {

612 int i, wds;

613 #ifdef ULLong

614 ULong *x;

615 ULLong carry, y;

616 #else

617 ULong carry, *x, y;

618 #ifdef Pack_32

619 ULong xi, z;

620 #endif

621 #endif

622 Bigint *b1;

623

624 wds = b->wds;

625 x = b->x;

626 i = 0;

627 carry = a;

628 do {

629 #ifdef ULLong

630 y = x (ULLong)m + carry;

631 carry = y >> 32;

632 *x++ = y & FFFFFFFF;

633 #else

634 #ifdef Pack_32

635 xi = *x;

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

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

638 carry = z >> 16;

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

640 #else

641 y = x m + carry;

642 carry = y >> 16;

643 *x++ = y & 0xffff;

644 #endif

645 #endif

646 }

647 while(++i < wds);

648 if (carry) {

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

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

651 Bcopy(b1, b);

652 Bfree(b);

653 b = b1;

654 }

655 b->x[wds++] = (ULong)carry;

656 b->wds = wds;

657 }

658 return b;

659 }

660

661 static Bigint *

662 s2b

663 #ifdef KR_headers

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

665 #else

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

667 #endif

668 {

669 Bigint *b;

670 int i, k;

671 Long x, y;

672

673 x = (nd + 8) / 9;

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

675 #ifdef Pack_32

676 b = Balloc(k);

677 b->x[0] = y9;

678 b->wds = 1;

679 #else

680 b = Balloc(k+1);

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

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

683 #endif

684

685 i = 9;

686 if (9 < nd0) {

687 s += 9;

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

689 while(++i < nd0);

690 s += dplen;

691 }

692 else

693 s += dplen + 9;

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

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

696 return b;

697 }

698

699 static int

700 hi0bits

701 #ifdef KR_headers

702 (x) ULong x;

703 #else

704 (ULong x)

705 #endif

706 {

707 int k = 0;

708

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

710 k = 16;

711 x <<= 16;

712 }

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

714 k += 8;

715 x <<= 8;

716 }

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

718 k += 4;

719 x <<= 4;

720 }

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

722 k += 2;

723 x <<= 2;

724 }

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

726 k++;

727 if (!(x & 0x40000000))

728 return 32;

729 }

730 return k;

731 }

732

733 static int

734 lo0bits

735 #ifdef KR_headers

736 (y) ULong *y;

737 #else

738 (ULong *y)

739 #endif

740 {

741 int k;

742 ULong x = *y;

743

744 if (x & 7) {

745 if (x & 1)

746 return 0;

747 if (x & 2) {

748 *y = x >> 1;

749 return 1;

750 }

751 *y = x >> 2;

752 return 2;

753 }

754 k = 0;

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

756 k = 16;

757 x >>= 16;

758 }

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

760 k += 8;

761 x >>= 8;

762 }

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

764 k += 4;

765 x >>= 4;

766 }

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

768 k += 2;

769 x >>= 2;

770 }

771 if (!(x & 1)) {

772 k++;

773 x >>= 1;

774 if (!x)

775 return 32;

776 }

777 *y = x;

778 return k;

779 }

780

781 static Bigint *

782 i2b

783 #ifdef KR_headers

784 (i) int i;

785 #else

786 (int i)

787 #endif

788 {

789 Bigint *b;

790

791 b = Balloc(1);

792 b->x[0] = i;

793 b->wds = 1;

794 return b;

795 }

796

797 static Bigint *

798 mult

799 #ifdef KR_headers

800 (a, b) Bigint a, b;

801 #else

802 (Bigint a, Bigint b)

803 #endif

804 {

805 Bigint *c;

806 int k, wa, wb, wc;

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

808 ULong y;

809 #ifdef ULLong

810 ULLong carry, z;

811 #else

812 ULong carry, z;

813 #ifdef Pack_32

814 ULong z2;

815 #endif

816 #endif

817

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

819 c = a;

820 a = b;

821 b = c;

822 }

823 k = a->k;

824 wa = a->wds;

825 wb = b->wds;

826 wc = wa + wb;

827 if (wc > a->maxwds)

828 k++;

829 c = Balloc(k);

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

831 *x = 0;

832 xa = a->x;

833 xae = xa + wa;

834 xb = b->x;

835 xbe = xb + wb;

836 xc0 = c->x;

837 #ifdef ULLong

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

839 y = *xb++;

840 if (y) {

841 x = xa;

842 xc = xc0;

843 carry = 0;

844 do {

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

846 carry = z >> 32;

847 *xc++ = z & FFFFFFFF;

848 }

849 while(x < xae);

850 *xc = (ULong)carry;

851 }

852 }

853 #else

854 #ifdef Pack_32

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

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

857 x = xa;

858 xc = xc0;

859 carry = 0;

860 do {

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

862 carry = z >> 16;

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

864 carry = z2 >> 16;

865 Storeinc(xc, z2, z);

866 }

867 while(x < xae);

868 *xc = carry;

869 }

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

871 x = xa;

872 xc = xc0;

873 carry = 0;

874 z2 = *xc;

875 do {

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

877 carry = z >> 16;

878 Storeinc(xc, z, z2);

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

880 carry = z2 >> 16;

881 }

882 while(x < xae);

883 *xc = z2;

884 }

885 }

886 #else

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

888 if (y = *xb++) {

889 x = xa;

890 xc = xc0;

891 carry = 0;

892 do {

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

894 carry = z >> 16;

895 *xc++ = z & 0xffff;

896 }

897 while(x < xae);

898 *xc = carry;

899 }

900 }

901 #endif

902 #endif

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

904 c->wds = wc;

905 return c;

906 }

907

908 static Bigint *p5s;

909

910 static Bigint *

911 pow5mult

912 #ifdef KR_headers

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

914 #else

915 (Bigint *b, int k)

916 #endif

917 {

918 Bigint b1, p5, *p51;

919 int i;

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

921

922 i = k & 3;

923 if (i)

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

925

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

927 return b;

928 p5 = p5s;

929 if (!p5) {

930 /* first time */

931 #ifdef MULTIPLE_THREADS

932 ACQUIRE_DTOA_LOCK(1);

933 p5 = p5s;

934 if (!p5) {

935 p5 = p5s = i2b(625);

936 p5->next = 0;

937 }

938 FREE_DTOA_LOCK(1);

939 #else

940 p5 = p5s = i2b(625);

941 p5->next = 0;

942 #endif

943 }

944 for(;;) {

945 if (k & 1) {

946 b1 = mult(b, p5);

947 Bfree(b);

948 b = b1;

949 }

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

951 break;

952 p51 = p5->next;

953 if (!p51) {

954 #ifdef MULTIPLE_THREADS

955 ACQUIRE_DTOA_LOCK(1);

956 p51 = p5->next;

957 if (!p51) {

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

959 p51->next = 0;

960 }

961 FREE_DTOA_LOCK(1);

962 #else

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

964 p51->next = 0;

965 #endif

966 }

967 p5 = p51;

968 }

969 return b;

970 }

971

972 static Bigint *

973 lshift

974 #ifdef KR_headers

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

976 #else

977 (Bigint *b, int k)

978 #endif

979 {

980 int i, k1, n, n1;

981 Bigint *b1;

982 ULong x, x1, *xe, z;

983

984 #ifdef Pack_32

985 n = k >> 5;

986 #else

987 n = k >> 4;

988 #endif

989 k1 = b->k;

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

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

992 k1++;

993 b1 = Balloc(k1);

994 x1 = b1->x;

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

996 *x1++ = 0;

997 x = b->x;

998 xe = x + b->wds;

999 #ifdef Pack_32

1000 if (k &= 0x1f) {

1001 k1 = 32 - k;

1002 z = 0;

1003 do {

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

1005 z = *x++ >> k1;

1006 }

1007 while(x < xe);

1008 *x1 = z;

1009 if (*x1)

1010 ++n1;

1011 }

1012 #else

1013 if (k &= 0xf) {

1014 k1 = 16 - k;

1015 z = 0;

1016 do {

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

1018 z = *x++ >> k1;

1019 }

1020 while(x < xe);

1021 if (*x1 = z)

1022 ++n1;

1023 }

1024 #endif

1025 else do

1026 x1++ = x++;

1027 while(x < xe);

1028 b1->wds = n1 - 1;

1029 Bfree(b);

1030 return b1;

1031 }

1032

1033 static int

1034 cmp

1035 #ifdef KR_headers

1036 (a, b) Bigint a, b;

1037 #else

1038 (Bigint a, Bigint b)

1039 #endif

1040 {

1041 ULong xa, xa0, xb, xb0;

1042 int i, j;

1043

1044 i = a->wds;

1045 j = b->wds;

1046 #ifdef DEBUG

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

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

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

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

1051 #endif

1052 if (i -= j)

1053 return i;

1054 xa0 = a->x;

1055 xa = xa0 + j;

1056 xb0 = b->x;

1057 xb = xb0 + j;

1058 for(;;) {

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

1060 return xa < xb ? -1 : 1;

1061 if (xa <= xa0)

1062 break;

1063 }

1064 return 0;

1065 }

1066

1067 static Bigint *

1068 diff

1069 #ifdef KR_headers

1070 (a, b) Bigint a, b;

1071 #else

1072 (Bigint a, Bigint b)

1073 #endif

1074 {

1075 Bigint *c;

1076 int i, wa, wb;

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

1078 #ifdef ULLong

1079 ULLong borrow, y;

1080 #else

1081 ULong borrow, y;

1082 #ifdef Pack_32

1083 ULong z;

1084 #endif

1085 #endif

1086

1087 i = cmp(a,b);

1088 if (!i) {

1089 c = Balloc(0);

1090 c->wds = 1;

1091 c->x[0] = 0;

1092 return c;

1093 }

1094 if (i < 0) {

1095 c = a;

1096 a = b;

1097 b = c;

1098 i = 1;

1099 }

1100 else

1101 i = 0;

1102 c = Balloc(a->k);

1103 c->sign = i;

1104 wa = a->wds;

1105 xa = a->x;

1106 xae = xa + wa;

1107 wb = b->wds;

1108 xb = b->x;

1109 xbe = xb + wb;

1110 xc = c->x;

1111 borrow = 0;

1112 #ifdef ULLong

1113 do {

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

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

1116 *xc++ = y & FFFFFFFF;

1117 }

1118 while(xb < xbe);

1119 while(xa < xae) {

1120 y = *xa++ - borrow;

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

1122 *xc++ = y & FFFFFFFF;

1123 }

1124 #else

1125 #ifdef Pack_32

1126 do {

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

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

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

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

1131 Storeinc(xc, z, y);

1132 }

1133 while(xb < xbe);

1134 while(xa < xae) {

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

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

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

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

1139 Storeinc(xc, z, y);

1140 }

1141 #else

1142 do {

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

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

1145 *xc++ = y & 0xffff;

1146 }

1147 while(xb < xbe);

1148 while(xa < xae) {

1149 y = *xa++ - borrow;

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

1151 *xc++ = y & 0xffff;

1152 }

1153 #endif

1154 #endif

1155 while(!*--xc)

1156 wa--;

1157 c->wds = wa;

1158 return c;

1159 }

1160

1161 static double

1162 ulp

1163 #ifdef KR_headers

1164 (x) U *x;

1165 #else

1166 (U *x)

1167 #endif

1168 {

1169 Long L;

1170 U u;

1171

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

1173 #ifndef Avoid_Underflow

1174 #ifndef Sudden_Underflow

1175 if (L > 0) {

1176 #endif

1177 #endif

1178 #ifdef IBM

1179 L \|= Exp_msk1 >> 4;

1180 #endif

1181 word0(&u) = L;

1182 word1(&u) = 0;

1183 #ifndef Avoid_Underflow

1184 #ifndef Sudden_Underflow

1185 }

1186 else {

1187 L = -L >> Exp_shift;

1188 if (L < Exp_shift) {

1189 word0(&u) = 0x80000 >> L;

1190 word1(&u) = 0;

1191 }

1192 else {

1193 word0(&u) = 0;

1194 L -= Exp_shift;

1195 word1(&u) = L >= 31 ? 1 : 1 << 31 - L;

1196 }

1197 }

1198 #endif

1199 #endif

1200 return dval(&u);

1201 }

1202

1203 static double

1204 b2d

1205 #ifdef KR_headers

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

1207 #else

1208 (Bigint a, int e)

1209 #endif

1210 {

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

1212 int k;

1213 U d;

1214 #ifdef VAX

1215 ULong d0, d1;

1216 #else

1217 #define d0 word0(&d)

1218 #define d1 word1(&d)

1219 #endif

1220

1221 xa0 = a->x;

1222 xa = xa0 + a->wds;

1223 y = *--xa;

1224 #ifdef DEBUG

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

1226 #endif

1227 k = hi0bits(y);

1228 *e = 32 - k;

1229 #ifdef Pack_32

1230 if (k < Ebits) {

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

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

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

1234 goto ret_d;

1235 }

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

1237 if (k -= Ebits) {

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

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

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

1241 }

1242 else {

1243 d0 = Exp_1 \| y;

1244 d1 = z;

1245 }

1246 #else

1247 if (k < Ebits + 16) {

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

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

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

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

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

1253 goto ret_d;

1254 }

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

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

1257 k -= Ebits + 16;

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

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

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

1261 #endif

1262 ret_d:

1263 #ifdef VAX

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

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

1266 #else

1267 #undef d0

1268 #undef d1

1269 #endif

1270 return dval(&d);

1271 }

1272

1273 static Bigint *

1274 d2b

1275 #ifdef KR_headers

1276 (d, e, bits) U d; int e, *bits;

1277 #else

1278 (U d, int e, int *bits)

1279 #endif

1280 {

1281 Bigint *b;

1282 int de, k;

1283 ULong *x, y, z;

1284 #ifndef Sudden_Underflow

1285 int i;

1286 #endif

1287 #ifdef VAX

1288 ULong d0, d1;

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

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

1291 #else

1292 #define d0 word0(d)

1293 #define d1 word1(d)

1294 #endif

1295

1296 #ifdef Pack_32

1297 b = Balloc(1);

1298 #else

1299 b = Balloc(2);

1300 #endif

1301 x = b->x;

1302

1303 z = d0 & Frac_mask;

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

1305 #ifdef Sudden_Underflow

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

1307 #ifndef IBM

1308 z \|= Exp_msk11;

1309 #endif

1310 #else

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

1312 if (de)

1313 z \|= Exp_msk1;

1314 #endif

1315 #ifdef Pack_32

1316 y = d1;

1317 if (y) {

1318 k = lo0bits(&y);

1319 if (k) {

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

1321 z >>= k;

1322 }

1323 else

1324 x[0] = y;

1325 x[1] = z;

1326 b->wds = x[1] ? 2 : 1;

1327 #ifndef Sudden_Underflow

1328 i = b->wds;

1329 #endif

1330 }

1331 else {

1332 k = lo0bits(&z);

1333 x[0] = z;

1334 #ifndef Sudden_Underflow

1335 i =

1336 #endif

1337 b->wds = 1;

1338 k += 32;

1339 }

1340 #else

1341 if (y = d1) {

1342 if (k = lo0bits(&y))

1343 if (k >= 16) {

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

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

1346 x[2] = z >> k;

1347 i = 2;

1348 }

1349 else {

1350 x[0] = y & 0xffff;

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

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

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

1354 i = 3;

1355 }

1356 else {

1357 x[0] = y & 0xffff;

1358 x[1] = y >> 16;

1359 x[2] = z & 0xffff;

1360 x[3] = z >> 16;

1361 i = 3;

1362 }

1363 }

1364 else {

1365 #ifdef DEBUG

1366 if (!z)

1367 Bug("Zero passed to d2b");

1368 #endif

1369 k = lo0bits(&z);

1370 if (k >= 16) {

1371 x[0] = z;

1372 i = 0;

1373 }

1374 else {

1375 x[0] = z & 0xffff;

1376 x[1] = z >> 16;

1377 i = 1;

1378 }

1379 k += 32;

1380 }

1381 while(!x[i])

1382 --i;

1383 b->wds = i + 1;

1384 #endif

1385 #ifndef Sudden_Underflow

1386 if (de) {

1387 #endif

1388 #ifdef IBM

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

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

1391 #else

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

1393 *bits = P - k;

1394 #endif

1395 #ifndef Sudden_Underflow

1396 }

1397 else {

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

1399 #ifdef Pack_32

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

1401 #else

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

1403 #endif

1404 }

1405 #endif

1406 return b;

1407 }

1408 #undef d0

1409 #undef d1

1410

1411 static double

1412 ratio

1413 #ifdef KR_headers

1414 (a, b) Bigint a, b;

1415 #else

1416 (Bigint a, Bigint b)

1417 #endif

1418 {

1419 U da, db;

1420 int k, ka, kb;

1421

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

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

1424 #ifdef Pack_32

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

1426 #else

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

1428 #endif

1429 #ifdef IBM

1430 if (k > 0) {

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

1432 if (k &= 3)

1433 dval(&da) *= 1 << k;

1434 }

1435 else {

1436 k = -k;

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

1438 if (k &= 3)

1439 dval(&db) *= 1 << k;

1440 }

1441 #else

1442 if (k > 0)

1443 word0(&da) += k*Exp_msk1;

1444 else {

1445 k = -k;

1446 word0(&db) += k*Exp_msk1;

1447 }

1448 #endif

1449 return dval(&da) / dval(&db);

1450 }

1451

1452 static CONST double

1453 tens[] = {

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

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

1456 1e20, 1e21, 1e22

1457 #ifdef VAX

1458 , 1e23, 1e24

1459 #endif

1460 };

1461

1462 static CONST double

1463 #ifdef IEEE_Arith

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

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

1466 #ifdef Avoid_Underflow

1467 9007199254740992.*9007199254740992.e-256

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

1469 #else

1470 1e-256

1471 #endif

1472 };

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

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

1475 #define Scale_Bit 0x10

1476 #define n_bigtens 5

1477 #else

1478 #ifdef IBM

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

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

1481 #define n_bigtens 3

1482 #else

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

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

1485 #define n_bigtens 2

1486 #endif

1487 #endif

1488

1489 #undef Need_Hexdig

1490 #ifdef INFNAN_CHECK

1491 #ifndef No_Hex_NaN

1492 #define Need_Hexdig

1493 #endif

1494 #endif

1495

1496 #ifndef Need_Hexdig

1497 #ifndef NO_HEX_FP

1498 #define Need_Hexdig

1499 #endif

1500 #endif

1501

1502 #ifdef Need_Hexdig /{/

1503 static unsigned char hexdig[256];

1504

1505 static void

1506 #ifdef KR_headers

1507 htinit(h, s, inc) unsigned char h; unsigned char s; int inc;

1508 #else

1509 htinit(unsigned char h, unsigned char s, int inc)

1510 #endif

1511 {

1512 int i, j;

1513 for(i = 0; (j = s[i]) !=0; i++)

1514 h[j] = (unsigned char)(i + inc);

1515 }

1516

1517 static void

1518 #ifdef KR_headers

1519 hexdig_init()

1520 #else

1521 hexdig_init(void)

1522 #endif

1523 {

1524 #define USC (unsigned char *)

1525 htinit(hexdig, USC "0123456789", 0x10);

1526 htinit(hexdig, USC "abcdef", 0x10 + 10);

1527 htinit(hexdig, USC "ABCDEF", 0x10 + 10);

1528 }

1529 #endif /* } Need_Hexdig */

1530

1531 #ifdef INFNAN_CHECK

1532

1533 #ifndef NAN_WORD0

1534 #define NAN_WORD0 0x7ff80000

1535 #endif

1536

1537 #ifndef NAN_WORD1

1538 #define NAN_WORD1 0

1539 #endif

1540

1541 static int

1542 match

1543 #ifdef KR_headers

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

1545 #else

1546 (CONST char *sp, CONST char t)

1547 #endif

1548 {

1549 int c, d;

1550 CONST char s = sp;

1551

1552 for(d = t++; d; d = t++) {

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

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

1555 if (c != d)

1556 return 0;

1557 }

1558 *sp = s + 1;

1559 return 1;

1560 }

1561

1562 #ifndef No_Hex_NaN

1563 static void

1564 hexnan

1565 #ifdef KR_headers

1566 (rvp, sp) U rvp; CONST char *sp;

1567 #else

1568 (U rvp, CONST char *sp)

1569 #endif

1570 {

1571 ULong c, x[2];

1572 CONST char *s;

1573 int c1, havedig, udx0, xshift;

1574

1575 if (!hexdig['0'])

1576 hexdig_init();

1577 x[0] = x[1] = 0;

1578 havedig = xshift = 0;

1579 udx0 = 1;

1580 s = *sp;

1581 /* allow optional initial 0x or 0X */

1582 for(c = (CONST unsigned char)(s+1); c && c <= ' '; c = (CONST unsigne d char)(s+1))

1583 ++s;

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

1585 s += 2;

1586 for(c = (CONST unsigned char)++s; c; c = (CONST unsigned char)++s) {

1587 c1 = hexdig[c];

1588 if (c1)

1589 c = c1 & 0xf;

1590 else if (c <= ' ') {

1591 if (udx0 && havedig) {

1592 udx0 = 0;

1593 xshift = 1;

1594 }

1595 continue;

1596 }

1597 #ifdef GDTOA_NON_PEDANTIC_NANCHECK

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

1599 *sp = s + 1;

1600 break;

1601 }

1602 else

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

1604 #else

1605 else {

1606 do {

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

1608 *sp = s + 1;

1609 break;

1610 }

1611 c = *++s;

1612 } while(c);

1613 break;

1614 }

1615 #endif

1616 havedig = 1;

1617 if (xshift) {

1618 xshift = 0;

1619 x[0] = x[1];

1620 x[1] = 0;

1621 }

1622 if (udx0)

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

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

1625 }

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

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

1628 word1(rvp) = x[1];

1629 }

1630 }

1631 #endif /No_Hex_NaN/

1632 #endif /* INFNAN_CHECK */

1633

1634 #ifdef Pack_32

1635 #define ULbits 32

1636 #define kshift 5

1637 #define kmask 31

1638 #else

1639 #define ULbits 16

1640 #define kshift 4

1641 #define kmask 15

1642 #endif

1643 #ifndef NO_HEX_FP /{/

1644

1645 static void

1646 #ifdef KR_headers

1647 rshift(b, k) Bigint *b; int k;

1648 #else

1649 rshift(Bigint *b, int k)

1650 #endif

1651 {

1652 ULong x, x1, *xe, y;

1653 int n;

1654

1655 x = x1 = b->x;

1656 n = k >> kshift;

1657 if (n < b->wds) {

1658 xe = x + b->wds;

1659 x += n;

1660 if (k &= kmask) {

1661 n = 32 - k;

1662 y = *x++ >> k;

1663 while(x < xe) {

1664 x1++ = (y \| (x << n)) & 0xffffffff;

1665 y = *x++ >> k;

1666 }

1667 if ((*x1 = y) !=0)

1668 x1++;

1669 }

1670 else

1671 while(x < xe)

1672 x1++ = x++;

1673 }

1674 if ((b->wds = x1 - b->x) == 0)

1675 b->x[0] = 0;

1676 }

1677

1678 static ULong

1679 #ifdef KR_headers

1680 any_on(b, k) Bigint *b; int k;

1681 #else

1682 any_on(Bigint *b, int k)

1683 #endif

1684 {

1685 int n, nwds;

1686 ULong x, x0, x1, x2;

1687

1688 x = b->x;

1689 nwds = b->wds;

1690 n = k >> kshift;

1691 if (n > nwds)

1692 n = nwds;

1693 else if (n < nwds && (k &= kmask)) {

1694 x1 = x2 = x[n];

1695 x1 >>= k;

1696 x1 <<= k;

1697 if (x1 != x2)

1698 return 1;

1699 }

1700 x0 = x;

1701 x += n;

1702 while(x > x0)

1703 if (*--x)

1704 return 1;

1705 return 0;

1706 }

1707

1708 enum { /* rounding values: same as FLT_ROUNDS */

1709 Round_zero = 0,

1710 Round_near = 1,

1711 Round_up = 2,

1712 Round_down = 3

1713 };

1714

1715 static Bigint *

1716 #ifdef KR_headers

1717 increment(b) Bigint *b;

1718 #else

1719 increment(Bigint *b)

1720 #endif

1721 {

1722 ULong x, xe;

1723 Bigint *b1;

1724

1725 x = b->x;

1726 xe = x + b->wds;

1727 do {

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

1729 ++*x;

1730 return b;

1731 }

1732 *x++ = 0;

1733 } while(x < xe);

1734 {

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

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

1737 Bcopy(b1,b);

1738 Bfree(b);

1739 b = b1;

1740 }

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

1742 }

1743 return b;

1744 }

1745

1746 void

1747 #ifdef KR_headers

1748 gethex(sp, rvp, rounding, sign)

1749 CONST char *sp; U rvp; int rounding, sign;

1750 #else

1751 gethex( CONST char *sp, U rvp, int rounding, int sign)

1752 #endif

1753 {

1754 Bigint *b;

1755 CONST unsigned char decpt, s0, s, s1;

1756 Long e, e1;

1757 ULong L, lostbits, *x;

1758 int big, denorm, esign, havedig, k, n, nbits, up, zret;

1759 #ifdef IBM

1760 int j;

1761 #endif

1762 enum {

1763 #ifdef IEEE_Arith /{{/

1764 emax = 0x7fe - Bias - P + 1,

1765 emin = Emin - P + 1

1766 #else /}{/

1767 emin = Emin - P,

1768 #ifdef VAX

1769 emax = 0x7ff - Bias - P + 1

1770 #endif

1771 #ifdef IBM

1772 emax = 0x7f - Bias - P

1773 #endif

1774 #endif /}}/

1775 };

1776 #ifdef USE_LOCALE

1777 int i;

1778 #ifdef NO_LOCALE_CACHE

1779 const unsigned char decimalpoint = (unsigned char)

1780 localeconv()->decimal_point;

1781 #else

1782 const unsigned char *decimalpoint;

1783 static unsigned char *decimalpoint_cache;

1784 if (!(s0 = decimalpoint_cache)) {

1785 s0 = (unsigned char*)localeconv()->decimal_point;

1786 if ((decimalpoint_cache = (unsigned char*)

1787 MALLOC(strlen((CONST char*)s0) + 1))) {

1788 strcpy((char)decimalpoint_cache, (CONST char)s0);

1789 s0 = decimalpoint_cache;

1790 }

1791 }

1792 decimalpoint = s0;

1793 #endif

1794 #endif

1795

1796 if (!hexdig['0'])

1797 hexdig_init();

1798 havedig = 0;

1799 s0 = (CONST unsigned char *)sp + 2;

1800 while(s0[havedig] == '0')

1801 havedig++;

1802 s0 += havedig;

1803 s = s0;

1804 decpt = 0;

1805 zret = 0;

1806 e = 0;

1807 if (hexdig[*s])

1808 havedig++;

1809 else {

1810 zret = 1;

1811 #ifdef USE_LOCALE

1812 for(i = 0; decimalpoint[i]; ++i) {

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

1814 goto pcheck;

1815 }

1816 decpt = s += i;

1817 #else

1818 if (*s != '.')

1819 goto pcheck;

1820 decpt = ++s;

1821 #endif

1822 if (!hexdig[*s])

1823 goto pcheck;

1824 while(*s == '0')

1825 s++;

1826 if (hexdig[*s])

1827 zret = 0;

1828 havedig = 1;

1829 s0 = s;

1830 }

1831 while(hexdig[*s])

1832 s++;

1833 #ifdef USE_LOCALE

1834 if (s == decimalpoint && !decpt) {

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

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

1837 goto pcheck;

1838 }

1839 decpt = s += i;

1840 #else

1841 if (*s == '.' && !decpt) {

1842 decpt = ++s;

1843 #endif

1844 while(hexdig[*s])

1845 s++;

1846 }/}/

1847 if (decpt)

1848 e = -(((Long)(s-decpt)) << 2);

1849 pcheck:

1850 s1 = s;

1851 big = esign = 0;

1852 switch(*s) {

1853 case 'p':

1854 case 'P':

1855 switch(*++s) {

1856 case '-':

1857 esign = 1;

1858 /* no break */

1859 case '+':

1860 s++;

1861 }

1862 if ((n = hexdig[*s]) == 0 \|\| n > 0x19) {

1863 s = s1;

1864 break;

1865 }

1866 e1 = n - 0x10;

1867 while((n = hexdig[*++s]) !=0 && n <= 0x19) {

1868 if (e1 & 0xf8000000)

1869 big = 1;

1870 e1 = 10*e1 + n - 0x10;

1871 }

1872 if (esign)

1873 e1 = -e1;

1874 e += e1;

1875 }

1876 sp = (char)s;

1877 if (!havedig)

1878 sp = (char)s0 - 1;

1879 if (zret)

1880 goto retz1;

1881 if (big) {

1882 if (esign) {

1883 #ifdef IEEE_Arith

1884 switch(rounding) {

1885 case Round_up:

1886 if (sign)

1887 break;

1888 goto ret_tiny;

1889 case Round_down:

1890 if (!sign)

1891 break;

1892 goto ret_tiny;

1893 }

1894 #endif

1895 goto retz;

1896 #ifdef IEEE_Arith

1897 ret_tiny:

1898 #ifndef NO_ERRNO

1899 errno = ERANGE;

1900 #endif

1901 word0(rvp) = 0;

1902 word1(rvp) = 1;

1903 return;

1904 #endif /* IEEE_Arith */

1905 }

1906 switch(rounding) {

1907 case Round_near:

1908 goto ovfl1;

1909 case Round_up:

1910 if (!sign)

1911 goto ovfl1;

1912 goto ret_big;

1913 case Round_down:

1914 if (sign)

1915 goto ovfl1;

1916 goto ret_big;

1917 }

1918 ret_big:

1919 word0(rvp) = Big0;

1920 word1(rvp) = Big1;

1921 return;

1922 }

1923 n = s1 - s0 - 1;

1924 for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1)

1925 k++;

1926 b = Balloc(k);

1927 x = b->x;

1928 n = 0;

1929 L = 0;

1930 #ifdef USE_LOCALE

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

1932 #endif

1933 while(s1 > s0) {

1934 #ifdef USE_LOCALE

1935 if (*--s1 == decimalpoint[i]) {

1936 s1 -= i;

1937 continue;

1938 }

1939 #else

1940 if (*--s1 == '.')

1941 continue;

1942 #endif

1943 if (n == ULbits) {

1944 *x++ = L;

1945 L = 0;

1946 n = 0;

1947 }

1948 L \|= (hexdig[*s1] & 0x0f) << n;

1949 n += 4;

1950 }

1951 *x++ = L;

1952 b->wds = n = x - b->x;

1953 n = ULbits*n - hi0bits(L);

1954 nbits = Nbits;

1955 lostbits = 0;

1956 x = b->x;

1957 if (n > nbits) {

1958 n -= nbits;

1959 if (any_on(b,n)) {

1960 lostbits = 1;

1961 k = n - 1;

1962 if (x[k>>kshift] & 1 << (k & kmask)) {

1963 lostbits = 2;

1964 if (k > 0 && any_on(b,k))

1965 lostbits = 3;

1966 }

1967 }

1968 rshift(b, n);

1969 e += n;

1970 }

1971 else if (n < nbits) {

1972 n = nbits - n;

1973 b = lshift(b, n);

1974 e -= n;

1975 x = b->x;

1976 }

1977 if (e > Emax) {

1978 ovfl:

1979 Bfree(b);

1980 ovfl1:

1981 #ifndef NO_ERRNO

1982 errno = ERANGE;

1983 #endif

1984 word0(rvp) = Exp_mask;

1985 word1(rvp) = 0;

1986 return;

1987 }

1988 denorm = 0;

1989 if (e < emin) {

1990 denorm = 1;

1991 n = emin - e;

1992 if (n >= nbits) {

1993 #ifdef IEEE_Arith /{/

1994 switch (rounding) {

1995 case Round_near:

1996 if (n == nbits && (n < 2 \|\| any_on(b,n-1)))

1997 goto ret_tiny;

1998 break;

1999 case Round_up:

2000 if (!sign)

2001 goto ret_tiny;

2002 break;

2003 case Round_down:

2004 if (sign)

2005 goto ret_tiny;

2006 }

2007 #endif /* } IEEE_Arith */

2008 Bfree(b);

2009 retz:

2010 #ifndef NO_ERRNO

2011 errno = ERANGE;

2012 #endif

2013 retz1:

2014 rvp->d = 0.;

2015 return;

2016 }

2017 k = n - 1;

2018 if (lostbits)

2019 lostbits = 1;

2020 else if (k > 0)

2021 lostbits = any_on(b,k);

2022 if (x[k>>kshift] & 1 << (k & kmask))

2023 lostbits \|= 2;

2024 nbits -= n;

2025 rshift(b,n);

2026 e = emin;

2027 }

2028 if (lostbits) {

2029 up = 0;

2030 switch(rounding) {

2031 case Round_zero:

2032 break;

2033 case Round_near:

2034 if (lostbits & 2

2035 && (lostbits & 1) \| (x[0] & 1))

2036 up = 1;

2037 break;

2038 case Round_up:

2039 up = 1 - sign;

2040 break;

2041 case Round_down:

2042 up = sign;

2043 }

2044 if (up) {

2045 k = b->wds;

2046 b = increment(b);

2047 x = b->x;

2048 if (denorm) {

2049 #if 0

2050 if (nbits == Nbits - 1

2051 && x[nbits >> kshift] & 1 << (nbits & kmask))

2052 denorm = 0; /* not currently used */

2053 #endif

2054 }

2055 else if (b->wds > k

2056 \|\| ((n = nbits & kmask) !=0

2057 && hi0bits(x[k-1]) < 32-n)) {

2058 rshift(b,1);

2059 if (++e > Emax)

2060 goto ovfl;

2061 }

2062 }

2063 }

2064 #ifdef IEEE_Arith

2065 if (denorm)

2066 word0(rvp) = b->wds > 1 ? b->x[1] & ~0x100000 : 0;

2067 else

2068 word0(rvp) = (b->x[1] & ~0x100000) \| ((e + 0x3ff + 52) << 20);

2069 word1(rvp) = b->x[0];

2070 #endif

2071 #ifdef IBM

2072 if ((j = e & 3)) {

2073 k = b->x[0] & ((1 << j) - 1);

2074 rshift(b,j);

2075 if (k) {

2076 switch(rounding) {

2077 case Round_up:

2078 if (!sign)

2079 increment(b);

2080 break;

2081 case Round_down:

2082 if (sign)

2083 increment(b);

2084 break;

2085 case Round_near:

2086 j = 1 << (j-1);

2087 if (k & j && ((k & (j-1)) \| lostbits))

2088 increment(b);

2089 }

2090 }

2091 }

2092 e >>= 2;

2093 word0(rvp) = b->x[1] \| ((e + 65 + 13) << 24);

2094 word1(rvp) = b->x[0];

2095 #endif

2096 #ifdef VAX

2097 /* The next two lines ignore swap of low- and high-order 2 bytes. */

2098 /* word0(rvp) = (b->x[1] & ~0x800000) \| ((e + 129 + 55) << 23); */

2099 /* word1(rvp) = b->x[0]; */

2100 word0(rvp) = ((b->x[1] & ~0x800000) >> 16) \| ((e + 129 + 55) << 7) \| (b- >x[1] << 16);

2101 word1(rvp) = (b->x[0] >> 16) \| (b->x[0] << 16);

2102 #endif

2103 Bfree(b);

2104 }

2105 #endif /}!NO_HEX_FP/

2106

2107 static int

2108 #ifdef KR_headers

2109 dshift(b, p2) Bigint *b; int p2;

2110 #else

2111 dshift(Bigint *b, int p2)

2112 #endif

2113 {

2114 int rv = hi0bits(b->x[b->wds-1]) - 4;

2115 if (p2 > 0)

2116 rv -= p2;

2117 return rv & kmask;

2118 }

2119

2120 static int

2121 quorem

2122 #ifdef KR_headers

2123 (b, S) Bigint b, S;

2124 #else

2125 (Bigint b, Bigint S)

2126 #endif

2127 {

2128 int n;

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

2130 #ifdef ULLong

2131 ULLong borrow, carry, y, ys;

2132 #else

2133 ULong borrow, carry, y, ys;

2134 #ifdef Pack_32

2135 ULong si, z, zs;

2136 #endif

2137 #endif

2138

2139 n = S->wds;

2140 #ifdef DEBUG

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

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

2143 #endif

2144 if (b->wds < n)

2145 return 0;

2146 sx = S->x;

2147 sxe = sx + --n;

2148 bx = b->x;

2149 bxe = bx + n;

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

2151 #ifdef DEBUG

2152 /debug/ if (q > 9)

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

2154 #endif

2155 if (q) {

2156 borrow = 0;

2157 carry = 0;

2158 do {

2159 #ifdef ULLong

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

2161 carry = ys >> 32;

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

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

2164 *bx++ = y & FFFFFFFF;

2165 #else

2166 #ifdef Pack_32

2167 si = *sx++;

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

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

2170 carry = zs >> 16;

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

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

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

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

2175 Storeinc(bx, z, y);

2176 #else

2177 ys = sx++ q + carry;

2178 carry = ys >> 16;

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

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

2181 *bx++ = y & 0xffff;

2182 #endif

2183 #endif

2184 }

2185 while(sx <= sxe);

2186 if (!*bxe) {

2187 bx = b->x;

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

2189 --n;

2190 b->wds = n;

2191 }

2192 }

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

2194 q++;

2195 borrow = 0;

2196 carry = 0;

2197 bx = b->x;

2198 sx = S->x;

2199 do {

2200 #ifdef ULLong

2201 ys = *sx++ + carry;

2202 carry = ys >> 32;

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

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

2205 *bx++ = y & FFFFFFFF;

2206 #else

2207 #ifdef Pack_32

2208 si = *sx++;

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

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

2211 carry = zs >> 16;

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

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

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

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

2216 Storeinc(bx, z, y);

2217 #else

2218 ys = *sx++ + carry;

2219 carry = ys >> 16;

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

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

2222 *bx++ = y & 0xffff;

2223 #endif

2224 #endif

2225 }

2226 while(sx <= sxe);

2227 bx = b->x;

2228 bxe = bx + n;

2229 if (!*bxe) {

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

2231 --n;

2232 b->wds = n;

2233 }

2234 }

2235 return q;

2236 }

2237

2238 #ifndef NO_STRTOD_BIGCOMP

2239

2240 static void

2241 bigcomp

2242 #ifdef KR_headers

2243 (rv, s0, bc)

2244 U rv; CONST char s0; BCinfo *bc;

2245 #else

2246 (U rv, CONST char s0, BCinfo *bc)

2247 #endif

2248 {

2249 Bigint b, d;

2250 int b2, bbits, d2, dd, dig, dsign, i, j, nd, nd0, p2, p5, speccase;

2251

2252 dsign = bc->dsign;

2253 nd = bc->nd;

2254 nd0 = bc->nd0;

2255 p5 = nd + bc->e0 - 1;

2256 dd = speccase = 0;

2257 #ifndef Sudden_Underflow

2258 if (rv->d == 0.) { /* special case: value near underflow-to-zero */

2259 /* threshold was rounded to zero */

2260 b = i2b(1);

2261 p2 = Emin - P + 1;

2262 bbits = 1;

2263 #ifdef Avoid_Underflow

2264 word0(rv) = (P+2) << Exp_shift;

2265 #else

2266 word1(rv) = 1;

2267 #endif

2268 i = 0;

2269 #ifdef Honor_FLT_ROUNDS

2270 if (bc->rounding == 1)

2271 #endif

2272 {

2273 speccase = 1;

2274 --p2;

2275 dsign = 0;

2276 goto have_i;

2277 }

2278 }

2279 else

2280 #endif

2281 b = d2b(rv, &p2, &bbits);

2282 #ifdef Avoid_Underflow

2283 p2 -= bc->scale;

2284 #endif

2285 /* floor(log2(rv)) == bbits - 1 + p2 */

2286 /* Check for denormal case. */

2287 i = P - bbits;

2288 if (i > (j = P - Emin - 1 + p2)) {

2289 #ifdef Sudden_Underflow

2290 Bfree(b);

2291 b = i2b(1);

2292 p2 = Emin;

2293 i = P - 1;

2294 #ifdef Avoid_Underflow

2295 word0(rv) = (1 + bc->scale) << Exp_shift;

2296 #else

2297 word0(rv) = Exp_msk1;

2298 #endif

2299 word1(rv) = 0;

2300 #else

2301 i = j;

2302 #endif

2303 }

2304 #ifdef Honor_FLT_ROUNDS

2305 if (bc->rounding != 1) {

2306 if (i > 0)

2307 b = lshift(b, i);

2308 if (dsign)

2309 b = increment(b);

2310 }

2311 else

2312 #endif

2313 {

2314 b = lshift(b, ++i);

2315 b->x[0] \|= 1;

2316 }

2317 #ifndef Sudden_Underflow

2318 have_i:

2319 #endif

2320 p2 -= p5 + i;

2321 d = i2b(1);

2322 /* Arrange for convenient computation of quotients:

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

2324 */

2325 if (p5 > 0)

2326 d = pow5mult(d, p5);

2327 else if (p5 < 0)

2328 b = pow5mult(b, -p5);

2329 if (p2 > 0) {

2330 b2 = p2;

2331 d2 = 0;

2332 }

2333 else {

2334 b2 = 0;

2335 d2 = -p2;

2336 }

2337 i = dshift(d, d2);

2338 if ((b2 += i) > 0)

2339 b = lshift(b, b2);

2340 if ((d2 += i) > 0)

2341 d = lshift(d, d2);

2342

2343 /* Now b/d = exactly half-way between the two floating-point values */

2344 /* on either side of the input string. Compute first digit of b/d. */

2345

2346 dig = quorem(b,d);

2347 if (!dig) {

2348 b = multadd(b, 10, 0); /* very unlikely */

2349 dig = quorem(b,d);

2350 }

2351

2352 /* Compare b/d with s0 */

2353

2354 for(i = 0; i < nd0; ) {

2355 dd = s0[i++] - '0' - dig;

2356 if (dd)

2357 goto ret;

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

2359 if (i < nd)

2360 dd = 1;

2361 goto ret;

2362 }

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

2364 dig = quorem(b,d);

2365 }

2366 for(j = bc->dp1; i++ < nd;) {

2367 dd = s0[j++] - '0' - dig;

2368 if (dd)

2369 goto ret;

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

2371 if (i < nd)

2372 dd = 1;

2373 goto ret;

2374 }

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

2376 dig = quorem(b,d);

2377 }

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

2379 dd = -1;

2380 ret:

2381 Bfree(b);

2382 Bfree(d);

2383 #ifdef Honor_FLT_ROUNDS

2384 if (bc->rounding != 1) {

2385 if (dd < 0) {

2386 if (bc->rounding == 0) {

2387 if (!dsign)

2388 goto retlow1;

2389 }

2390 else if (dsign)

2391 goto rethi1;

2392 }

2393 else if (dd > 0) {

2394 if (bc->rounding == 0) {

2395 if (dsign)

2396 goto rethi1;

2397 goto ret1;

2398 }

2399 if (!dsign)

2400 goto rethi1;

2401 dval(rv) += 2.*ulp(rv);

2402 }

2403 else {

2404 bc->inexact = 0;

2405 if (dsign)

2406 goto rethi1;

2407 }

2408 }

2409 else

2410 #endif

2411 if (speccase) {

2412 if (dd <= 0)

2413 rv->d = 0.;

2414 }

2415 else if (dd < 0) {

2416 if (!dsign) /* does not happen for round-near */

2417 retlow1:

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

2419 }

2420 else if (dd > 0) {

2421 if (dsign) {

2422 rethi1:

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

2424 }

2425 }

2426 else {

2427 /* Exact half-way case: apply round-even rule. */

2428 if (word1(rv) & 1) {

2429 if (dsign)

2430 goto rethi1;

2431 goto retlow1;

2432 }

2433 }

2434

2435 #ifdef Honor_FLT_ROUNDS

2436 ret1:

2437 #endif

2438 return;

2439 }

2440 #endif /* NO_STRTOD_BIGCOMP */

2441

2442 double

2443 strtod

2444 #ifdef KR_headers

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

2446 #else

2447 (CONST char s00, char *se)

2448 #endif

2449 {

2450 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, e, e1;

2451 int esign, i, j, k, nd, nd0, nf, nz, nz0, sign;

2452 CONST char s, s0, *s1;

2453 double aadj, aadj1;

2454 Long L;

2455 U aadj2, adj, rv, rv0;

2456 ULong y, z;

2457 BCinfo bc;

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

2459 #ifdef SET_INEXACT

2460 int oldinexact;

2461 #endif

2462 #ifdef Honor_FLT_ROUNDS /{/

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

2464 bc.rounding = Flt_Rounds;

2465 #else /}{/

2466 bc.rounding = 1;

2467 switch(fegetround()) {

2468 case FE_TOWARDZERO: bc.rounding = 0; break;

2469 case FE_UPWARD: bc.rounding = 2; break;

2470 case FE_DOWNWARD: bc.rounding = 3;

2471 }

2472 #endif /}}/

2473 #endif /}/

2474 #ifdef USE_LOCALE

2475 CONST char *s2;

2476 #endif

2477

2478 sign = nz0 = nz = bc.dplen = bc.uflchk = 0;

2479 dval(&rv) = 0.;

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

2481 case '-':

2482 sign = 1;

2483 /* no break */

2484 case '+':

2485 if (*++s)

2486 goto break2;

2487 /* no break */

2488 case 0:

2489 goto ret0;

2490 case '\t':

2491 case '\n':

2492 case '\v':

2493 case '\f':

2494 case '\r':

2495 case ' ':

2496 continue;

2497 default:

2498 goto break2;

2499 }

2500 break2:

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

2502 #ifndef NO_HEX_FP /{/

2503 switch(s[1]) {

2504 case 'x':

2505 case 'X':

2506 #ifdef Honor_FLT_ROUNDS

2507 gethex(&s, &rv, bc.rounding, sign);

2508 #else

2509 gethex(&s, &rv, 1, sign);

2510 #endif

2511 goto ret;

2512 }

2513 #endif /}/

2514 nz0 = 1;

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

2516 if (!*s)

2517 goto ret;

2518 }

2519 s0 = s;

2520 y = z = 0;

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

2522 if (nd < 9)

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

2524 else if (nd < 16)

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

2526 nd0 = nd;

2527 bc.dp0 = bc.dp1 = s - s0;

2528 #ifdef USE_LOCALE

2529 s1 = localeconv()->decimal_point;

2530 if (c == *s1) {

2531 c = '.';

2532 if (*++s1) {

2533 s2 = s;

2534 for(;;) {

2535 if (++s2 != s1) {

2536 c = 0;

2537 break;

2538 }

2539 if (!*++s1) {

2540 s = s2;

2541 break;

2542 }

2543 }

2544 }

2545 }

2546 #endif

2547 if (c == '.') {

2548 c = *++s;

2549 bc.dp1 = s - s0;

2550 bc.dplen = bc.dp1 - bc.dp0;

2551 if (!nd) {

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

2553 nz++;

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

2555 s0 = s;

2556 nf += nz;

2557 nz = 0;

2558 goto have_dig;

2559 }

2560 goto dig_done;

2561 }

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

2563 have_dig:

2564 nz++;

2565 if (c -= '0') {

2566 nf += nz;

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

2568 if (nd++ < 9)

2569 y *= 10;

2570 else if (nd <= DBL_DIG + 1)

2571 z *= 10;

2572 if (nd++ < 9)

2573 y = 10*y + c;

2574 else if (nd <= DBL_DIG + 1)

2575 z = 10*z + c;

2576 nz = 0;

2577 }

2578 }

2579 }

2580 dig_done:

2581 e = 0;

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

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

2584 goto ret0;

2585 }

2586 s00 = s;

2587 esign = 0;

2588 switch(c = *++s) {

2589 case '-':

2590 esign = 1;

2591 case '+':

2592 c = *++s;

2593 }

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

2595 while(c == '0')

2596 c = *++s;

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

2598 L = c - '0';

2599 s1 = s;

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

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

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

2603 /* Avoid confusion from exponents

2604 * so large that e might overflow.

2605 */

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

2607 else

2608 e = (int)L;

2609 if (esign)

2610 e = -e;

2611 }

2612 else

2613 e = 0;

2614 }

2615 else

2616 s = s00;

2617 }

2618 if (!nd) {

2619 if (!nz && !nz0) {

2620 #ifdef INFNAN_CHECK

2621 /* Check for Nan and Infinity */

2622 if (!bc.dplen)

2623 switch(c) {

2624 case 'i':

2625 case 'I':

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

2627 --s;

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

2629 ++s;

2630 word0(&rv) = 0x7ff00000;

2631 word1(&rv) = 0;

2632 goto ret;

2633 }

2634 break;

2635 case 'n':

2636 case 'N':

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

2638 word0(&rv) = NAN_WORD0;

2639 word1(&rv) = NAN_WORD1;

2640 #ifndef No_Hex_NaN

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

2642 hexnan(&rv, &s);

2643 #endif

2644 goto ret;

2645 }

2646 }

2647 #endif /* INFNAN_CHECK */

2648 ret0:

2649 s = s00;

2650 sign = 0;

2651 }

2652 goto ret;

2653 }

2654 bc.e0 = e1 = e -= nf;

2655

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

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

2658 * after is the integer represented by those digits times

2659 * 10*e /

2660

2661 if (!nd0)

2662 nd0 = nd;

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

2664 dval(&rv) = y;

2665 if (k > 9) {

2666 #ifdef SET_INEXACT

2667 if (k > DBL_DIG)

2668 oldinexact = get_inexact();

2669 #endif

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

2671 }

2672 bd0 = 0;

2673 if (nd <= DBL_DIG

2674 #ifndef RND_PRODQUOT

2675 #ifndef Honor_FLT_ROUNDS

2676 && Flt_Rounds == 1

2677 #endif

2678 #endif

2679 ) {

2680 if (!e)

2681 goto ret;

2682 if (e > 0) {

2683 if (e <= Ten_pmax) {

2684 #ifdef VAX

2685 goto vax_ovfl_check;

2686 #else

2687 #ifdef Honor_FLT_ROUNDS

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

2689 if (sign) {

2690 rv.d = -rv.d;

2691 sign = 0;

2692 }

2693 #endif

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

2695 goto ret;

2696 #endif

2697 }

2698 i = DBL_DIG - nd;

2699 if (e <= Ten_pmax + i) {

2700 /* A fancier test would sometimes let us do

2701 * this for larger i values.

2702 */

2703 #ifdef Honor_FLT_ROUNDS

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

2705 if (sign) {

2706 rv.d = -rv.d;

2707 sign = 0;

2708 }

2709 #endif

2710 e -= i;

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

2712 #ifdef VAX

2713 /* VAX exponent range is so narrow we must

2714 * worry about overflow here...

2715 */

2716 vax_ovfl_check:

2717 word0(&rv) -= P*Exp_msk1;

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

2719 if ((word0(&rv) & Exp_mask)

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

2721 goto ovfl;

2722 word0(&rv) += P*Exp_msk1;

2723 #else

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

2725 #endif

2726 goto ret;

2727 }

2728 }

2729 #ifndef Inaccurate_Divide

2730 else if (e >= -Ten_pmax) {

2731 #ifdef Honor_FLT_ROUNDS

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

2733 if (sign) {

2734 rv.d = -rv.d;

2735 sign = 0;

2736 }

2737 #endif

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

2739 goto ret;

2740 }

2741 #endif

2742 }

2743 e1 += nd - k;

2744

2745 #ifdef IEEE_Arith

2746 #ifdef SET_INEXACT

2747 bc.inexact = 1;

2748 if (k <= DBL_DIG)

2749 oldinexact = get_inexact();

2750 #endif

2751 #ifdef Avoid_Underflow

2752 bc.scale = 0;

2753 #endif

2754 #ifdef Honor_FLT_ROUNDS

2755 if (bc.rounding >= 2) {

2756 if (sign)

2757 bc.rounding = bc.rounding == 2 ? 0 : 2;

2758 else

2759 if (bc.rounding != 2)

2760 bc.rounding = 0;

2761 }

2762 #endif

2763 #endif /IEEE_Arith/

2764

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

2766

2767 if (e1 > 0) {

2768 i = e1 & 15;

2769 if (i)

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

2771 if (e1 &= ~15) {

2772 if (e1 > DBL_MAX_10_EXP) {

2773 ovfl:

2774 #ifndef NO_ERRNO

2775 errno = ERANGE;

2776 #endif

2777 /* Can't trust HUGE_VAL */

2778 #ifdef IEEE_Arith

2779 #ifdef Honor_FLT_ROUNDS

2780 switch(bc.rounding) {

2781 case 0: /* toward 0 */

2782 case 3: /* toward -infinity */

2783 word0(&rv) = Big0;

2784 word1(&rv) = Big1;

2785 break;

2786 default:

2787 word0(&rv) = Exp_mask;

2788 word1(&rv) = 0;

2789 }

2790 #else /Honor_FLT_ROUNDS/

2791 word0(&rv) = Exp_mask;

2792 word1(&rv) = 0;

2793 #endif /Honor_FLT_ROUNDS/

2794 #ifdef SET_INEXACT

2795 /* set overflow bit */

2796 dval(&rv0) = 1e300;

2797 dval(&rv0) *= dval(&rv0);

2798 #endif

2799 #else /IEEE_Arith/

2800 word0(&rv) = Big0;

2801 word1(&rv) = Big1;

2802 #endif /IEEE_Arith/

2803 goto ret;

2804 }

2805 e1 >>= 4;

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

2807 if (e1 & 1)

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

2809 /* The last multiplication could overflow. */

2810 word0(&rv) -= P*Exp_msk1;

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

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

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

2814 goto ovfl;

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

2816 /* set to largest number */

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

2818 word0(&rv) = Big0;

2819 word1(&rv) = Big1;

2820 }

2821 else

2822 word0(&rv) += P*Exp_msk1;

2823 }

2824 }

2825 else if (e1 < 0) {

2826 e1 = -e1;

2827 i = e1 & 15;

2828 if (i)

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

2830 if (e1 >>= 4) {

2831 if (e1 >= 1 << n_bigtens)

2832 goto undfl;

2833 #ifdef Avoid_Underflow

2834 if (e1 & Scale_Bit)

2835 bc.scale = 2*P;

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

2837 if (e1 & 1)

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

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

2840 >> Exp_shift)) > 0) {

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

2842 if (j >= 32) {

2843 word1(&rv) = 0;

2844 if (j >= 53)

2845 word0(&rv) = (P+2)*Exp_msk1;

2846 else

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

2848 }

2849 else

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

2851 }

2852 #else

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

2854 if (e1 & 1)

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

2856 /* The last multiplication could underflow. */

2857 dval(&rv0) = dval(&rv);

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

2859 if (!dval(&rv)) {

2860 dval(&rv) = 2.*dval(&rv0);

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

2862 #endif

2863 if (!dval(&rv)) {

2864 undfl:

2865 dval(&rv) = 0.;

2866 #ifndef NO_ERRNO

2867 errno = ERANGE;

2868 #endif

2869 goto ret;

2870 }

2871 #ifndef Avoid_Underflow

2872 word0(&rv) = Tiny0;

2873 word1(&rv) = Tiny1;

2874 /* The refinement below will clean

2875 * this approximation up.

2876 */

2877 }

2878 #endif

2879 }

2880 }

2881

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

2883

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

2885

2886 bc.nd = nd;

2887 #ifndef NO_STRTOD_BIGCOMP

2888 bc.nd0 = nd0; /* Only needed if nd > strtod_diglim, but done here */

2889 /* to silence an erroneous warning about bc.nd0 */

2890 /* possibly not being initialized. */

2891 if (nd > strtod_diglim) {

2892 /* ASSERT(strtod_diglim >= 18); 18 == one more than the */

2893 /* minimum number of decimal digits to distinguish double values */

2894 /* in IEEE arithmetic. */

2895 i = j = 18;

2896 if (i > nd0)

2897 j += bc.dplen;

2898 for(;;) {

2899 if (--j <= bc.dp1 && j >= bc.dp0)

2900 j = bc.dp0 - 1;

2901 if (s0[j] != '0')

2902 break;

2903 --i;

2904 }

2905 e += nd - i;

2906 nd = i;

2907 if (nd0 > nd)

2908 nd0 = nd;

2909 if (nd < 9) { /* must recompute y */

2910 y = 0;

2911 for(i = 0; i < nd0; ++i)

2912 y = 10*y + s0[i] - '0';

2913 for(j = bc.dp1; i < nd; ++i)

2914 y = 10*y + s0[j++] - '0';

2915 }

2916 }

2917 #endif

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

2919

2920 for(;;) {

2921 bd = Balloc(bd0->k);

2922 Bcopy(bd, bd0);

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

2924 bs = i2b(1);

2925

2926 if (e >= 0) {

2927 bb2 = bb5 = 0;

2928 bd2 = bd5 = e;

2929 }

2930 else {

2931 bb2 = bb5 = -e;

2932 bd2 = bd5 = 0;

2933 }

2934 if (bbe >= 0)

2935 bb2 += bbe;

2936 else

2937 bd2 -= bbe;

2938 bs2 = bb2;

2939 #ifdef Honor_FLT_ROUNDS

2940 if (bc.rounding != 1)

2941 bs2++;

2942 #endif

2943 #ifdef Avoid_Underflow

2944 j = bbe - bc.scale;

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

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

2947 j += P - Emin;

2948 else

2949 j = P + 1 - bbbits;

2950 #else /Avoid_Underflow/

2951 #ifdef Sudden_Underflow

2952 #ifdef IBM

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

2954 #else

2955 j = P + 1 - bbbits;

2956 #endif

2957 #else /Sudden_Underflow/

2958 j = bbe;

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

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

2961 j += P - Emin;

2962 else

2963 j = P + 1 - bbbits;

2964 #endif /Sudden_Underflow/

2965 #endif /Avoid_Underflow/

2966 bb2 += j;

2967 bd2 += j;

2968 #ifdef Avoid_Underflow

2969 bd2 += bc.scale;

2970 #endif

2971 i = bb2 < bd2 ? bb2 : bd2;

2972 if (i > bs2)

2973 i = bs2;

2974 if (i > 0) {

2975 bb2 -= i;

2976 bd2 -= i;

2977 bs2 -= i;

2978 }

2979 if (bb5 > 0) {

2980 bs = pow5mult(bs, bb5);

2981 bb1 = mult(bs, bb);

2982 Bfree(bb);

2983 bb = bb1;

2984 }

2985 if (bb2 > 0)

2986 bb = lshift(bb, bb2);

2987 if (bd5 > 0)

2988 bd = pow5mult(bd, bd5);

2989 if (bd2 > 0)

2990 bd = lshift(bd, bd2);

2991 if (bs2 > 0)

2992 bs = lshift(bs, bs2);

2993 delta = diff(bb, bd);

2994 bc.dsign = delta->sign;

2995 delta->sign = 0;

2996 i = cmp(delta, bs);

2997 #ifndef NO_STRTOD_BIGCOMP

2998 if (bc.nd > nd && i <= 0) {

2999 if (bc.dsign)

3000 break; /* Must use bigcomp(). */

3001 #ifdef Honor_FLT_ROUNDS

3002 if (bc.rounding != 1) {

3003 if (i < 0)

3004 break;

3005 }

3006 else

3007 #endif

3008 {

3009 bc.nd = nd;

3010 i = -1; /* Discarded digits make delta smaller. */

3011 }

3012 }

3013 #endif

3014 #ifdef Honor_FLT_ROUNDS

3015 if (bc.rounding != 1) {

3016 if (i < 0) {

3017 /* Error is less than an ulp */

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

3019 /* exact */

3020 #ifdef SET_INEXACT

3021 bc.inexact = 0;

3022 #endif

3023 break;

3024 }

3025 if (bc.rounding) {

3026 if (bc.dsign) {

3027 adj.d = 1.;

3028 goto apply_adj;

3029 }

3030 }

3031 else if (!bc.dsign) {

3032 adj.d = -1.;

3033 if (!word1(&rv)

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

3035 y = word0(&rv) & Exp_mask;

3036 #ifdef Avoid_Underflow

3037 if (!bc.scale \|\| y > 2PExp_msk 1)

3038 #else

3039 if (y)

3040 #endif

3041 {

3042 delta = lshift(delta,Log2P);

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

3044 adj.d = -0.5;

3045 }

3046 }

3047 apply_adj:

3048 #ifdef Avoid_Underflow

3049 if (bc.scale && (y = word0(&rv) & Exp_ma sk)

3050 <= 2PExp_msk1)

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

3052 #else

3053 #ifdef Sudden_Underflow

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

3055 P*Exp_msk1) {

3056 word0(&rv) += P*Exp_msk1;

3057 dval(&rv) += adj.d*ulp(dval(&rv) );

3058 word0(&rv) -= P*Exp_msk1;

3059 }

3060 else

3061 #endif /Sudden_Underflow/

3062 #endif /Avoid_Underflow/

3063 dval(&rv) += adj.d*ulp(&rv);

3064 }

3065 break;

3066 }

3067 adj.d = ratio(delta, bs);

3068 if (adj.d < 1.)

3069 adj.d = 1.;

3070 if (adj.d <= 0x7ffffffe) {

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

3072 y = adj.d;

3073 if (y != adj.d) {

3074 if (!((bc.rounding>>1) ^ bc.dsign))

3075 y++;

3076 adj.d = y;

3077 }

3078 }

3079 #ifdef Avoid_Underflow

3080 if (bc.scale && (y = word0(&rv) & Exp_mask) <= 2PExp_m sk1)

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

3082 #else

3083 #ifdef Sudden_Underflow

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

3085 word0(&rv) += P*Exp_msk1;

3086 adj.d *= ulp(dval(&rv));

3087 if (bc.dsign)

3088 dval(&rv) += adj.d;

3089 else

3090 dval(&rv) -= adj.d;

3091 word0(&rv) -= P*Exp_msk1;

3092 goto cont;

3093 }

3094 #endif /Sudden_Underflow/

3095 #endif /Avoid_Underflow/

3096 adj.d *= ulp(&rv);

3097 if (bc.dsign) {

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

3099 goto ovfl;

3100 dval(&rv) += adj.d;

3101 }

3102 else

3103 dval(&rv) -= adj.d;

3104 goto cont;

3105 }

3106 #endif /Honor_FLT_ROUNDS/

3107

3108 if (i < 0) {

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

3110 * special case of mantissa a power of two.

3111 */

3112 if (bc.dsign \|\| word1(&rv) \|\| word0(&rv) & Bndry_mask

3113 #ifdef IEEE_Arith

3114 #ifdef Avoid_Underflow

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

3116 #else

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

3118 #endif

3119 #endif

3120 ) {

3121 #ifdef SET_INEXACT

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

3123 bc.inexact = 0;

3124 #endif

3125 break;

3126 }

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

3128 /* exact result */

3129 #ifdef SET_INEXACT

3130 bc.inexact = 0;

3131 #endif

3132 break;

3133 }

3134 delta = lshift(delta,Log2P);

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

3136 goto drop_down;

3137 break;

3138 }

3139 if (i == 0) {

3140 /* exactly half-way between */

3141 if (bc.dsign) {

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

3143 && word1(&rv) == (

3144 #ifdef Avoid_Underflow

3145 (bc.scale && (y = word0(&rv) & Exp_mask) <= 2PExp_msk1 )

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

3147 #endif

3148 0xffffffff)) {

3149 /boundary case -- increment exponent/

3150 word0(&rv) = (word0(&rv) & Exp_mask)

3151 + Exp_msk1

3152 #ifdef IBM

3153 \| Exp_msk1 >> 4

3154 #endif

3155 ;

3156 word1(&rv) = 0;

3157 #ifdef Avoid_Underflow

3158 bc.dsign = 0;

3159 #endif

3160 break;

3161 }

3162 }

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

3164 drop_down:

3165 /* boundary case -- decrement exponent */

3166 #ifdef Sudden_Underflow /{{/

3167 L = word0(&rv) & Exp_mask;

3168 #ifdef IBM

3169 if (L < Exp_msk1)

3170 #else

3171 #ifdef Avoid_Underflow

3172 if (L <= (bc.scale ? (2P+1)Exp_msk1 : Exp_msk1 ))

3173 #else

3174 if (L <= Exp_msk1)

3175 #endif /Avoid_Underflow/

3176 #endif /IBM/

3177 {

3178 if (bc.nd >nd) {

3179 bc.uflchk = 1;

3180 break;

3181 }

3182 goto undfl;

3183 }

3184 L -= Exp_msk1;

3185 #else /Sudden_Underflow}{/

3186 #ifdef Avoid_Underflow

3187 if (bc.scale) {

3188 L = word0(&rv) & Exp_mask;

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

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

3191 /* round even ==> */

3192 /* accept rv */

3193 break;

3194 /* rv = smallest denormal */

3195 if (bc.nd >nd) {

3196 bc.uflchk = 1;

3197 break;

3198 }

3199 goto undfl;

3200 }

3201 }

3202 #endif /Avoid_Underflow/

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

3204 #endif /Sudden_Underflow}}/

3205 word0(&rv) = L \| Bndry_mask1;

3206 word1(&rv) = 0xffffffff;

3207 #ifdef IBM

3208 goto cont;

3209 #else

3210 break;

3211 #endif

3212 }

3213 #ifndef ROUND_BIASED

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

3215 break;

3216 #endif

3217 if (bc.dsign)

3218 dval(&rv) += ulp(&rv);

3219 #ifndef ROUND_BIASED

3220 else {

3221 dval(&rv) -= ulp(&rv);

3222 #ifndef Sudden_Underflow

3223 if (!dval(&rv)) {

3224 if (bc.nd >nd) {

3225 bc.uflchk = 1;

3226 break;

3227 }

3228 goto undfl;

3229 }

3230 #endif

3231 }

3232 #ifdef Avoid_Underflow

3233 bc.dsign = 1 - bc.dsign;

3234 #endif

3235 #endif

3236 break;

3237 }

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

3239 if (bc.dsign)

3240 aadj = aadj1 = 1.;

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

3242 #ifndef Sudden_Underflow

3243 if (word1(&rv) == Tiny1 && !word0(&rv)) {

3244 if (bc.nd >nd) {

3245 bc.uflchk = 1;

3246 break;

3247 }

3248 goto undfl;

3249 }

3250 #endif

3251 aadj = 1.;

3252 aadj1 = -1.;

3253 }

3254 else {

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

3256 /* rounded down... */

3257

3258 if (aadj < 2./FLT_RADIX)

3259 aadj = 1./FLT_RADIX;

3260 else

3261 aadj *= 0.5;

3262 aadj1 = -aadj;

3263 }

3264 }

3265 else {

3266 aadj *= 0.5;

3267 aadj1 = bc.dsign ? aadj : -aadj;

3268 #ifdef Check_FLT_ROUNDS

3269 switch(bc.rounding) {

3270 case 2: /* towards +infinity */

3271 aadj1 -= 0.5;

3272 break;

3273 case 0: /* towards 0 */

3274 case 3: /* towards -infinity */

3275 aadj1 += 0.5;

3276 }

3277 #else

3278 if (Flt_Rounds == 0)

3279 aadj1 += 0.5;

3280 #endif /Check_FLT_ROUNDS/

3281 }

3282 y = word0(&rv) & Exp_mask;

3283

3284 /* Check for overflow */

3285

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

3287 dval(&rv0) = dval(&rv);

3288 word0(&rv) -= P*Exp_msk1;

3289 adj.d = aadj1 * ulp(&rv);

3290 dval(&rv) += adj.d;

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

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

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

3294 goto ovfl;

3295 word0(&rv) = Big0;

3296 word1(&rv) = Big1;

3297 goto cont;

3298 }

3299 else

3300 word0(&rv) += P*Exp_msk1;

3301 }

3302 else {

3303 #ifdef Avoid_Underflow

3304 if (bc.scale && y <= 2PExp_msk1) {

3305 if (aadj <= 0x7fffffff) {

3306 if ((z = (ULong)aadj) <= 0)

3307 z = 1;

3308 aadj = z;

3309 aadj1 = bc.dsign ? aadj : -aadj;

3310 }

3311 dval(&aadj2) = aadj1;

3312 word0(&aadj2) += (2P+1)Exp_msk1 - y;

3313 aadj1 = dval(&aadj2);

3314 }

3315 adj.d = aadj1 * ulp(&rv);

3316 dval(&rv) += adj.d;

3317 #else

3318 #ifdef Sudden_Underflow

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

3320 dval(&rv0) = dval(&rv);

3321 word0(&rv) += P*Exp_msk1;

3322 adj.d = aadj1 * ulp(&rv);

3323 dval(&rv) += adj.d;

3324 #ifdef IBM

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

3326 #else

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

3328 #endif

3329 {

3330 if (word0(&rv0) == Tiny0

3331 && word1(&rv0) == Tiny1) {

3332 if (bc.nd >nd) {

3333 bc.uflchk = 1;

3334 break;

3335 }

3336 goto undfl;

3337 }

3338 word0(&rv) = Tiny0;

3339 word1(&rv) = Tiny1;

3340 goto cont;

3341 }

3342 else

3343 word0(&rv) -= P*Exp_msk1;

3344 }

3345 else {

3346 adj.d = aadj1 * ulp(&rv);

3347 dval(&rv) += adj.d;

3348 }

3349 #else /Sudden_Underflow/

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

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

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

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

3354 * trouble from bits lost to denormalization;

3355 * example: 1.2e-307 .

3356 */

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

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

3359 if (!bc.dsign)

3360 aadj1 = -aadj1;

3361 }

3362 adj.d = aadj1 * ulp(&rv);

3363 dval(&rv) += adj.d;

3364 #endif /Sudden_Underflow/

3365 #endif /Avoid_Underflow/

3366 }

3367 z = word0(&rv) & Exp_mask;

3368 #ifndef SET_INEXACT

3369 if (bc.nd == nd) {

3370 #ifdef Avoid_Underflow

3371 if (!bc.scale)

3372 #endif

3373 if (y == z) {

3374 /* Can we stop now? */

3375 L = (Long)aadj;

3376 aadj -= L;

3377 /* The tolerances below are conservative. */

3378 if (bc.dsign \|\| word1(&rv) \|\| word0(&rv) & Bndry_mask) {

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

3380 break;

3381 }

3382 else if (aadj < .4999999/FLT_RADIX)

3383 break;

3384 }

3385 }

3386 #endif

3387 cont:

3388 Bfree(bb);

3389 Bfree(bd);

3390 Bfree(bs);

3391 Bfree(delta);

3392 }

3393 Bfree(bb);

3394 Bfree(bd);

3395 Bfree(bs);

3396 Bfree(bd0);

3397 Bfree(delta);

3398 #ifndef NO_STRTOD_BIGCOMP

3399 if (bc.nd > nd)

3400 bigcomp(&rv, s0, &bc);

3401 #endif

3402 #ifdef SET_INEXACT

3403 if (bc.inexact) {

3404 if (!oldinexact) {

3405 word0(&rv0) = Exp_1 + (70 << Exp_shift);

3406 word1(&rv0) = 0;

3407 dval(&rv0) += 1.;

3408 }

3409 }

3410 else if (!oldinexact)

3411 clear_inexact();

3412 #endif

3413 #ifdef Avoid_Underflow

3414 if (bc.scale) {

3415 word0(&rv0) = Exp_1 - 2PExp_msk1;

3416 word1(&rv0) = 0;

3417 dval(&rv) *= dval(&rv0);

3418 #ifndef NO_ERRNO

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

3420 #ifdef IEEE_Arith

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

3422 #else

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

3424 #endif

3425 errno = ERANGE;

3426 #endif

3427 }

3428 #endif /* Avoid_Underflow */

3429 #ifdef SET_INEXACT

3430 if (bc.inexact && !(word0(&rv) & Exp_mask)) {

3431 /* set underflow bit */

3432 dval(&rv0) = 1e-300;

3433 dval(&rv0) *= dval(&rv0);

3434 }

3435 #endif

3436 ret:

3437 if (se)

3438 se = (char )s;

3439 return sign ? -dval(&rv) : dval(&rv);

3440 }

3441

3442 #ifndef MULTIPLE_THREADS

3443 static char *dtoa_result;

3444 #endif

3445

3446 static char *

3447 #ifdef KR_headers

3448 rv_alloc(i) int i;

3449 #else

3450 rv_alloc(int i)

3451 #endif

3452 {

3453 int j, k, *r;

3454

3455 j = sizeof(ULong);

3456 for(k = 0;

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

3458 j <<= 1)

3459 k++;

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

3461 *r = k;

3462 return

3463 #ifndef MULTIPLE_THREADS

3464 dtoa_result =

3465 #endif

3466 (char *)(r+1);

3467 }

3468

3469 static char *

3470 #ifdef KR_headers

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

3472 #else

3473 nrv_alloc(CONST char s, char *rve, int n)

3474 #endif

3475 {

3476 char rv, t;

3477

3478 t = rv = rv_alloc(n);

3479 for(t = s++; t; t = *s++) t++;

3480 if (rve)

3481 *rve = t;

3482 return rv;

3483 }

3484

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

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

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

3488 * when MULTIPLE_THREADS is not defined.

3489 */

3490

3491 void

3492 #ifdef KR_headers

3493 freedtoa(s) char *s;

3494 #else

3495 freedtoa(char *s)

3496 #endif

3497 {

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

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

3500 Bfree(b);

3501 #ifndef MULTIPLE_THREADS

3502 if (s == dtoa_result)

3503 dtoa_result = 0;

3504 #endif

3505 }

3506

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

3508 *

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

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

3511 *

3512 * Modifications:

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

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

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

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

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

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

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

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

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

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

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

3524 * as would satisfaction of the stopping test with strict

3525 * inequality.

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

3527 * quantities.

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

3529 * we use floating-point arithmetic rather than resorting

3530 * to multiple-precision integers.

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

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

3533 * multiple-precision integer arithmetic only if we cannot

3534 * guarantee that the floating-point calculation has given

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

3536 * "uniformly" distributed input, the probability is

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

3538 * calculation.

3539 */

3540

3541 char *

3542 dtoa

3543 #ifdef KR_headers

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

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

3546 #else

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

3548 #endif

3549 {

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

3551 of ecvt and fcvt; trailing zeros are suppressed from

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

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

3554 then *decpt is set to 9999.

3555

3556 mode:

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

3558 and rounded to nearest.

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

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

3561 1e23 whereas mode 1 gives 9.999999999999999e22.

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

3563 return value similar to that of ecvt, except

3564 that trailing zeros are suppressed.

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

3566 gives a return value similar to that from fcvt,

3567 except that trailing zeros are suppressed, and

3568 ndigits can be negative.

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

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

3571 possibly return a shorter string that rounds to d.

3572 With IEEE arithmetic and compilation with

3573 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same

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

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

3576 fast floating-point estimate (if applicable).

3577

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

3579

3580 Sufficient space is allocated to the return value

3581 to hold the suppressed trailing zeros.

3582 */

3583

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

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

3586 spec_case, try_quick;

3587 Long L;

3588 #ifndef Sudden_Underflow

3589 int denorm;

3590 ULong x;

3591 #endif

3592 Bigint b, b1, delta, mlo = NULL, mhi, S;

3593 U d2, eps, u;

3594 double ds;

3595 char s, s0;

3596 #ifdef SET_INEXACT

3597 int inexact, oldinexact;

3598 #endif

3599 #ifdef Honor_FLT_ROUNDS /{/

3600 int Rounding;

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

3602 Rounding = Flt_Rounds;

3603 #else /}{/

3604 Rounding = 1;

3605 switch(fegetround()) {

3606 case FE_TOWARDZERO: Rounding = 0; break;

3607 case FE_UPWARD: Rounding = 2; break;

3608 case FE_DOWNWARD: Rounding = 3;

3609 }

3610 #endif /}}/

3611 #endif /}/

3612

3613 #ifndef MULTIPLE_THREADS

3614 if (dtoa_result) {

3615 freedtoa(dtoa_result);

3616 dtoa_result = 0;

3617 }

3618 #endif

3619

3620 u.d = dd;

3621 if (word0(&u) & Sign_bit) {

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

3623 *sign = 1;

3624 word0(&u) &= ~Sign_bit; /* clear sign bit */

3625 }

3626 else

3627 *sign = 0;

3628

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

3630 #ifdef IEEE_Arith

3631 if ((word0(&u) & Exp_mask) == Exp_mask)

3632 #else

3633 if (word0(&u) == 0x8000)

3634 #endif

3635 {

3636 /* Infinity or NaN */

3637 *decpt = 9999;

3638 #ifdef IEEE_Arith

3639 if (!word1(&u) && !(word0(&u) & 0xfffff))

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

3641 #endif

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

3643 }

3644 #endif

3645 #ifdef IBM

3646 dval(&u) += 0; /* normalize */

3647 #endif

3648 if (!dval(&u)) {

3649 *decpt = 1;

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

3651 }

3652

3653 #ifdef SET_INEXACT

3654 try_quick = oldinexact = get_inexact();

3655 inexact = 1;

3656 #endif

3657 #ifdef Honor_FLT_ROUNDS

3658 if (Rounding >= 2) {

3659 if (*sign)

3660 Rounding = Rounding == 2 ? 0 : 2;

3661 else

3662 if (Rounding != 2)

3663 Rounding = 0;

3664 }

3665 #endif

3666

3667 b = d2b(&u, &be, &bbits);

3668 i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1));

3669 #ifndef Sudden_Underflow

3670 if (i) {

3671 #endif

3672 dval(&d2) = dval(&u);

3673 word0(&d2) &= Frac_mask1;

3674 word0(&d2) \|= Exp_11;

3675 #ifdef IBM

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

3677 dval(&d2) /= 1 << j;

3678 #endif

3679

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

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

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

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

3684 *

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

3686 *

3687 * k = (i - Bias)*0.301029995663981

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

3689 *

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

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

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

3693 * to compensate for any error in the multiplication of

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

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

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

3697 * Hence we adjust the constant term to 0.1760912590558.

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

3699 * but this is probably not worthwhile.)

3700 */

3701

3702 i -= Bias;

3703 #ifdef IBM

3704 i <<= 2;

3705 i += j;

3706 #endif

3707 #ifndef Sudden_Underflow

3708 denorm = 0;

3709 }

3710 else {

3711 /* d is denormalized */

3712

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

3714 x = i > 32 ? word0(&u) << (64 - i) \| word1(&u) >> (i - 32)

3715 : word1(&u) << (32 - i);

3716 dval(&d2) = x;

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

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

3719 denorm = 1;

3720 }

3721 #endif

3722 ds = (dval(&d2)-1.5)0.289529654602168 + 0.1760912590558 + i0.301029995 663981;

3723 k = (int)ds;

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

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

3726 k_check = 1;

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

3728 if (dval(&u) < tens[k])

3729 k--;

3730 k_check = 0;

3731 }

3732 j = bbits - i - 1;

3733 if (j >= 0) {

3734 b2 = 0;

3735 s2 = j;

3736 }

3737 else {

3738 b2 = -j;

3739 s2 = 0;

3740 }

3741 if (k >= 0) {

3742 b5 = 0;

3743 s5 = k;

3744 s2 += k;

3745 }

3746 else {

3747 b2 -= k;

3748 b5 = -k;

3749 s5 = 0;

3750 }

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

3752 mode = 0;

3753

3754 #ifndef SET_INEXACT

3755 #ifdef Check_FLT_ROUNDS

3756 try_quick = Rounding == 1;

3757 #else

3758 try_quick = 1;

3759 #endif

3760 #endif /SET_INEXACT/

3761

3762 if (mode > 5) {

3763 mode -= 4;

3764 try_quick = 0;

3765 }

3766 leftright = 1;

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

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

3769 switch(mode) {

3770 case 0:

3771 case 1:

3772 i = 18;

3773 ndigits = 0;

3774 break;

3775 case 2:

3776 leftright = 0;

3777 /* no break */

3778 case 4:

3779 if (ndigits <= 0)

3780 ndigits = 1;

3781 ilim = ilim1 = i = ndigits;

3782 break;

3783 case 3:

3784 leftright = 0;

3785 /* no break */

3786 case 5:

3787 i = ndigits + k + 1;

3788 ilim = i;

3789 ilim1 = i - 1;

3790 if (i <= 0)

3791 i = 1;

3792 }

3793 s = s0 = rv_alloc(i);

3794

3795 #ifdef Honor_FLT_ROUNDS

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

3797 leftright = 0;

3798 #endif

3799

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

3801

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

3803

3804 i = 0;

3805 dval(&d2) = dval(&u);

3806 k0 = k;

3807 ilim0 = ilim;

3808 ieps = 2; /* conservative */

3809 if (k > 0) {

3810 ds = tens[k&0xf];

3811 j = k >> 4;

3812 if (j & Bletch) {

3813 /* prevent overflows */

3814 j &= Bletch - 1;

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

3816 ieps++;

3817 }

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

3819 if (j & 1) {

3820 ieps++;

3821 ds *= bigtens[i];

3822 }

3823 dval(&u) /= ds;

3824 }

3825 else {

3826 j1 = -k;

3827 if (j1) {

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

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

3830 if (j & 1) {

3831 ieps++;

3832 dval(&u) *= bigtens[i];

3833 }

3834 }

3835 }

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

3837 if (ilim1 <= 0)

3838 goto fast_failed;

3839 ilim = ilim1;

3840 k--;

3841 dval(&u) *= 10.;

3842 ieps++;

3843 }

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

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

3846 if (ilim == 0) {

3847 S = mhi = 0;

3848 dval(&u) -= 5.;

3849 if (dval(&u) > dval(&eps))

3850 goto one_digit;

3851 if (dval(&u) < -dval(&eps))

3852 goto no_digits;

3853 goto fast_failed;

3854 }

3855 #ifndef No_leftright

3856 if (leftright) {

3857 /* Use Steele & White method of only

3858 * generating digits needed.

3859 */

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

3861 for(i = 0;;) {

3862 L = (long)dval(&u);

3863 dval(&u) -= L;

3864 *s++ = '0' + (char)L;

3865 if (dval(&u) < dval(&eps))

3866 goto ret1;

3867 if (1. - dval(&u) < dval(&eps))

3868 goto bump_up;

3869 if (++i >= ilim)

3870 break;

3871 dval(&eps) *= 10.;

3872 dval(&u) *= 10.;

3873 }

3874 }

3875 else {

3876 #endif

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

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

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

3880 L = (Long)(dval(&u));

3881 if (!(dval(&u) -= L))

3882 ilim = i;

3883 *s++ = '0' + (char)L;

3884 if (i == ilim) {

3885 if (dval(&u) > 0.5 + dval(&eps))

3886 goto bump_up;

3887 else if (dval(&u) < 0.5 - dval(&eps)) {

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

3889 s++;

3890 goto ret1;

3891 }

3892 break;

3893 }

3894 }

3895 #ifndef No_leftright

3896 }

3897 #endif

3898 fast_failed:

3899 s = s0;

3900 dval(&u) = dval(&d2);

3901 k = k0;

3902 ilim = ilim0;

3903 }

3904

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

3906

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

3908 /* Yes. */

3909 ds = tens[k];

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

3911 S = mhi = 0;

3912 if (ilim < 0 \|\| dval(&u) <= 5*ds)

3913 goto no_digits;

3914 goto one_digit;

3915 }

3916 for(i = 1; i <= k + 1; i++, dval(&u) *= 10.) {

3917 L = (Long)(dval(&u) / ds);

3918 dval(&u) -= L*ds;

3919 #ifdef Check_FLT_ROUNDS

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

3921 if (dval(&u) < 0) {

3922 L--;

3923 dval(&u) += ds;

3924 }

3925 #endif

3926 *s++ = '0' + (char)L;

3927 if (!dval(&u)) {

3928 #ifdef SET_INEXACT

3929 inexact = 0;

3930 #endif

3931 break;

3932 }

3933 if (i == ilim) {

3934 #ifdef Honor_FLT_ROUNDS

3935 if (mode > 1)

3936 switch(Rounding) {

3937 case 0: goto ret1;

3938 case 2: goto bump_up;

3939 }

3940 #endif

3941 dval(&u) += dval(&u);

3942 if (dval(&u) > ds \|\| (dval(&u) == ds && L & 1)) {

3943 bump_up:

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

3945 if (s == s0) {

3946 k++;

3947 *s = '0';

3948 break;

3949 }

3950 ++*s++;

3951 }

3952 break;

3953 }

3954 }

3955 goto ret1;

3956 }

3957

3958 m2 = b2;

3959 m5 = b5;

3960 mhi = mlo = 0;

3961 if (leftright) {

3962 i =

3963 #ifndef Sudden_Underflow

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

3965 #endif

3966 #ifdef IBM

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

3968 #else

3969 1 + P - bbits;

3970 #endif

3971 b2 += i;

3972 s2 += i;

3973 mhi = i2b(1);

3974 }

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

3976 i = m2 < s2 ? m2 : s2;

3977 b2 -= i;

3978 m2 -= i;

3979 s2 -= i;

3980 }

3981 if (b5 > 0) {

3982 if (leftright) {

3983 if (m5 > 0) {

3984 mhi = pow5mult(mhi, m5);

3985 b1 = mult(mhi, b);

3986 Bfree(b);

3987 b = b1;

3988 }

3989 j = b5 - m5;

3990 if (j)

3991 b = pow5mult(b, j);

3992 }

3993 else

3994 b = pow5mult(b, b5);

3995 }

3996 S = i2b(1);

3997 if (s5 > 0)

3998 S = pow5mult(S, s5);

3999

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

4001

4002 spec_case = 0;

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

4004 #ifdef Honor_FLT_ROUNDS

4005 && Rounding == 1

4006 #endif

4007 ) {

4008 if (!word1(&u) && !(word0(&u) & Bndry_mask)

4009 #ifndef Sudden_Underflow

4010 && word0(&u) & (Exp_mask & ~Exp_msk1)

4011 #endif

4012 ) {

4013 /* The special case */

4014 b2 += Log2P;

4015 s2 += Log2P;

4016 spec_case = 1;

4017 }

4018 }

4019

4020 /* Arrange for convenient computation of quotients:

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

4022 *

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

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

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

4026 */

4027 #ifdef Pack_32

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

4029 if (i)

4030 i = 32 - i;

4031 #define iInc 28

4032 #else

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

4034 i = 16 - i;

4035 #define iInc 12

4036 #endif

4037 i = dshift(S, s2);

4038 b2 += i;

4039 m2 += i;

4040 s2 += i;

4041 if (b2 > 0)

4042 b = lshift(b, b2);

4043 if (s2 > 0)

4044 S = lshift(S, s2);

4045 if (k_check) {

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

4047 k--;

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

4049 if (leftright)

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

4051 ilim = ilim1;

4052 }

4053 }

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

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

4056 /* no digits, fcvt style */

4057 no_digits:

4058 k = -1 - ndigits;

4059 goto ret;

4060 }

4061 one_digit:

4062 *s++ = '1';

4063 k++;

4064 goto ret;

4065 }

4066 if (leftright) {

4067 if (m2 > 0)

4068 mhi = lshift(mhi, m2);

4069

4070 /* Compute mlo -- check for special case

4071 * that d is a normalized power of 2.

4072 */

4073

4074 mlo = mhi;

4075 if (spec_case) {

4076 mhi = Balloc(mhi->k);

4077 Bcopy(mhi, mlo);

4078 mhi = lshift(mhi, Log2P);

4079 }

4080

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

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

4083 /* Do we yet have the shortest decimal string

4084 * that will round to d?

4085 */

4086 j = cmp(b, mlo);

4087 delta = diff(S, mhi);

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

4089 Bfree(delta);

4090 #ifndef ROUND_BIASED

4091 if (j1 == 0 && mode != 1 && !(word1(&u) & 1)

4092 #ifdef Honor_FLT_ROUNDS

4093 && Rounding >= 1

4094 #endif

4095 ) {

4096 if (dig == '9')

4097 goto round_9_up;

4098 if (j > 0)

4099 dig++;

4100 #ifdef SET_INEXACT

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

4102 inexact = 0;

4103 #endif

4104 *s++ = (char)dig;

4105 goto ret;

4106 }

4107 #endif

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

4109 #ifndef ROUND_BIASED

4110 && !(word1(&u) & 1)

4111 #endif

4112 )) {

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

4114 #ifdef SET_INEXACT

4115 inexact = 0;

4116 #endif

4117 goto accept_dig;

4118 }

4119 #ifdef Honor_FLT_ROUNDS

4120 if (mode > 1)

4121 switch(Rounding) {

4122 case 0: goto accept_dig;

4123 case 2: goto keep_dig;

4124 }

4125 #endif /Honor_FLT_ROUNDS/

4126 if (j1 > 0) {

4127 b = lshift(b, 1);

4128 j1 = cmp(b, S);

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

4130 && dig++ == '9')

4131 goto round_9_up;

4132 }

4133 accept_dig:

4134 *s++ = (char)dig;

4135 goto ret;

4136 }

4137 if (j1 > 0) {

4138 #ifdef Honor_FLT_ROUNDS

4139 if (!Rounding)

4140 goto accept_dig;

4141 #endif

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

4143 round_9_up:

4144 *s++ = '9';

4145 goto roundoff;

4146 }

4147 *s++ = (char)dig + 1;

4148 goto ret;

4149 }

4150 #ifdef Honor_FLT_ROUNDS

4151 keep_dig:

4152 #endif

4153 *s++ = (char)dig;

4154 if (i == ilim)

4155 break;

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

4157 if (mlo == mhi)

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

4159 else {

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

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

4162 }

4163 }

4164 }

4165 else

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

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

4168 *s++ = (char)dig;

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

4170 #ifdef SET_INEXACT

4171 inexact = 0;

4172 #endif

4173 goto ret;

4174 }

4175 if (i >= ilim)

4176 break;

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

4178 }

4179

4180 /* Round off last digit */

4181

4182 #ifdef Honor_FLT_ROUNDS

4183 switch(Rounding) {

4184 case 0: goto trimzeros;

4185 case 2: goto roundoff;

4186 }

4187 #endif

4188 b = lshift(b, 1);

4189 j = cmp(b, S);

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

4191 roundoff:

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

4193 if (s == s0) {

4194 k++;

4195 *s++ = '1';

4196 goto ret;

4197 }

4198 ++*s++;

4199 }

4200 else {

4201 #ifdef Honor_FLT_ROUNDS

4202 trimzeros:

4203 #endif

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

4205 s++;

4206 }

4207 ret:

4208 Bfree(S);

4209 if (mhi) {

4210 if (mlo && mlo != mhi)

4211 Bfree(mlo);

4212 Bfree(mhi);

4213 }

4214 ret1:

4215 #ifdef SET_INEXACT

4216 if (inexact) {

4217 if (!oldinexact) {

4218 word0(&u) = Exp_1 + (70 << Exp_shift);

4219 word1(&u) = 0;

4220 dval(&u) += 1.;

4221 }

4222 }

4223 else if (!oldinexact)

4224 clear_inexact();

4225 #endif

4226 Bfree(b);

4227 *s = 0;

4228 *decpt = k + 1;

4229 if (rve)

4230 *rve = s;

4231 return s0;

4232 }

4233

4234 } // namespace dmg_fp

OLD	NEW

« no previous file with comments | « base/third_party/dmg_fp/dmg_fp.h ('k') | base/third_party/dmg_fp/dtoa_wrapper.cc » ('j') | no next file with comments »