nspr/pr/src/misc/prdtoa.c - Issue 2078763002: Delete bundled copy of NSS and replace with README.

Side by Side Diff: nspr/pr/src/misc/prdtoa.c

Issue 2078763002: Delete bundled copy of NSS and replace with README. (Closed) Base URL: https://chromium.googlesource.com/chromium/deps/nss@master

Patch Set: Delete bundled copy of NSS and replace with README. Created 4 years, 6 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:

View unified diff | Download patch

OLD	NEW
	(Empty)
1 /* -- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -- */

2 /* This Source Code Form is subject to the terms of the Mozilla Public

3 * License, v. 2.0. If a copy of the MPL was not distributed with this

4 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */

5

6 /*

7 * This file is based on the third-party code dtoa.c. We minimize our

8 * modifications to third-party code to make it easy to merge new versions.

9 * The author of dtoa.c was not willing to add the parentheses suggested by

10 * GCC, so we suppress these warnings.

11 */

12 #if (__GNUC__ > 4) \|\| (__GNUC__ == 4 && __GNUC_MINOR__ >= 2)

13 #pragma GCC diagnostic ignored "-Wparentheses"

14 #endif

15

16 #include "primpl.h"

17 #include "prbit.h"

18

19 #define MULTIPLE_THREADS

20 #define ACQUIRE_DTOA_LOCK(n) PR_Lock(dtoa_lock[n])

21 #define FREE_DTOA_LOCK(n) PR_Unlock(dtoa_lock[n])

22

23 static PRLock *dtoa_lock[2];

24

25 void _PR_InitDtoa(void)

26 {

27 dtoa_lock[0] = PR_NewLock();

28 dtoa_lock[1] = PR_NewLock();

29 }

30

31 void _PR_CleanupDtoa(void)

32 {

33 PR_DestroyLock(dtoa_lock[0]);

34 dtoa_lock[0] = NULL;

35 PR_DestroyLock(dtoa_lock[1]);

36 dtoa_lock[1] = NULL;

37

38 /* FIXME: deal with freelist and p5s. */

39 }

40

41 #if !defined(__ARM_EABI__) \

42 && (defined(__arm) \|\| defined(__arm__) \|\| defined(__arm26__) \

43 \|\| defined(__arm32__))

44 #define IEEE_ARM

45 #elif defined(IS_LITTLE_ENDIAN)

46 #define IEEE_8087

47 #else

48 #define IEEE_MC68k

49 #endif

50

51 #define Long PRInt32

52 #define ULong PRUint32

53 #define NO_LONG_LONG

54

55 #define No_Hex_NaN

56

57 /****************************************************************

58 *

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

60 *

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

62 *

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

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

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

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

67 * documentation for such software.

68 *

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

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

71 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY

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

73 *

74 ***************************************************************/

75

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

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

78

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

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

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

82 * of) Intel 80x87 arithmetic, the call

83 * _control87(PC_53, MCW_PC);

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

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

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

87 * file.

88 */

89

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

91 *

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

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

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

95 * biased rounding (add half and chop).

96 *

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

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

99 *

100 * Modifications:

101 *

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

103 * arithmetic (not IEEE double-extended).

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

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

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

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

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

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

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

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

112 * arithmetic to determine the adjustment to within

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

114 * compute a second residual.

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

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

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

118 */

119

120 /*

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

122 * significant byte has the lowest address.

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

124 * significant byte has the lowest address.

125 * #define IEEE_ARM for IEEE-arithmetic machines where the two words

126 * in a double are stored in big endian order but the two shorts

127 * in a word are still stored in little endian order.

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

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

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

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

132 * computation of dtoa.

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

134 * and strtod and dtoa should round accordingly.

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

136 * and Honor_FLT_ROUNDS is not #defined.

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

138 * that use extended-precision instructions to compute rounded

139 * products and quotients) with IBM.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

155 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.

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

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

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

159 * directly -- and assumed always to succeed. Similarly, if you

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

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

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

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

164 * mode 3 with thousands of digits requested.)

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

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

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

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

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

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

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

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

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

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

175 * pointers, PRIVATE_MEM >= 7112 appears adequate.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

190 * and NAN_WORD1 are used instead.

191 * #define MULTIPLE_THREADS if the system offers preemptively scheduled

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

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

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

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

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

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

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

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

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

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

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

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

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

205 * Sudden_Underflow.

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

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

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

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

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

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

212 * suitable wrapper, that defines two functions,

213 * int get_inexact(void);

214 * void clear_inexact(void);

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

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

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

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

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

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

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

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

223 */

224

225 #ifndef Long

226 #define Long long

227 #endif

228 #ifndef ULong

229 typedef unsigned Long ULong;

230 #endif

231

232 #ifdef DEBUG

233 #include "stdio.h"

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

235 #endif

236

237 #include "stdlib.h"

238 #include "string.h"

239

240 #ifdef USE_LOCALE

241 #include "locale.h"

242 #endif

243

244 #ifdef MALLOC

245 #ifdef KR_headers

246 extern char *MALLOC();

247 #else

248 extern void *MALLOC(size_t);

249 #endif

250 #else

251 #define MALLOC malloc

252 #endif

253

254 #ifndef Omit_Private_Memory

255 #ifndef PRIVATE_MEM

256 #define PRIVATE_MEM 2304

257 #endif

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

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

260 #endif

261

262 #undef IEEE_Arith

263 #undef Avoid_Underflow

264 #ifdef IEEE_MC68k

265 #define IEEE_Arith

266 #endif

267 #ifdef IEEE_8087

268 #define IEEE_Arith

269 #endif

270 #ifdef IEEE_ARM

271 #define IEEE_Arith

272 #endif

273

274 #include "errno.h"

275

276 #ifdef Bad_float_h

277

278 #ifdef IEEE_Arith

279 #define DBL_DIG 15

280 #define DBL_MAX_10_EXP 308

281 #define DBL_MAX_EXP 1024

282 #define FLT_RADIX 2

283 #endif /IEEE_Arith/

284

285 #ifdef IBM

286 #define DBL_DIG 16

287 #define DBL_MAX_10_EXP 75

288 #define DBL_MAX_EXP 63

289 #define FLT_RADIX 16

290 #define DBL_MAX 7.2370055773322621e+75

291 #endif

292

293 #ifdef VAX

294 #define DBL_DIG 16

295 #define DBL_MAX_10_EXP 38

296 #define DBL_MAX_EXP 127

297 #define FLT_RADIX 2

298 #define DBL_MAX 1.7014118346046923e+38

299 #endif

300

301 #ifndef LONG_MAX

302 #define LONG_MAX 2147483647

303 #endif

304

305 #else /* ifndef Bad_float_h */

306 #include "float.h"

307 /*

308 * MacOS 10.2 defines the macro FLT_ROUNDS to an internal function

309 * which does not exist on 10.1. We can safely #define it to 1 here

310 * to allow 10.2 builds to run on 10.1, since we can't use fesetround()

311 * (which does not exist on 10.1 either).

312 */

313 #if defined(XP_MACOSX) && (!defined(MAC_OS_X_VERSION_10_2) \|\| \

314 MAC_OS_X_VERSION_MIN_REQUIRED < MAC_OS_X_VERSION_10_2)

315 #undef FLT_ROUNDS

316 #define FLT_ROUNDS 1

317 #endif /* DT < 10.2 */

318 #endif /* Bad_float_h */

319

320 #ifndef __MATH_H__

321 #include "math.h"

322 #endif

323

324 #ifdef __cplusplus

325 extern "C" {

326 #endif

327

328 #ifndef CONST

329 #ifdef KR_headers

330 #define CONST /* blank */

331 #else

332 #define CONST const

333 #endif

334 #endif

335

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

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

338 #endif

339

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

341

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

343 #ifdef IEEE_8087

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

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

346 #else

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

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

349 #endif

350

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

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

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

354 */

355 #if defined(IEEE_8087) + defined(IEEE_ARM) + defined(VAX)

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

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

358 #else

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

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

361 #endif

362

363 /* #define P DBL_MANT_DIG */

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

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

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

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

368

369 #ifdef IEEE_Arith

370 #define Exp_shift 20

371 #define Exp_shift1 20

372 #define Exp_msk1 0x100000

373 #define Exp_msk11 0x100000

374 #define Exp_mask 0x7ff00000

375 #define P 53

376 #define Bias 1023

377 #define Emin (-1022)

378 #define Exp_1 0x3ff00000

379 #define Exp_11 0x3ff00000

380 #define Ebits 11

381 #define Frac_mask 0xfffff

382 #define Frac_mask1 0xfffff

383 #define Ten_pmax 22

384 #define Bletch 0x10

385 #define Bndry_mask 0xfffff

386 #define Bndry_mask1 0xfffff

387 #define LSB 1

388 #define Sign_bit 0x80000000

389 #define Log2P 1

390 #define Tiny0 0

391 #define Tiny1 1

392 #define Quick_max 14

393 #define Int_max 14

394 #ifndef NO_IEEE_Scale

395 #define Avoid_Underflow

396 #ifdef Flush_Denorm /* debugging option */

397 #undef Sudden_Underflow

398 #endif

399 #endif

400

401 #ifndef Flt_Rounds

402 #ifdef FLT_ROUNDS

403 #define Flt_Rounds FLT_ROUNDS

404 #else

405 #define Flt_Rounds 1

406 #endif

407 #endif /Flt_Rounds/

408

409 #ifdef Honor_FLT_ROUNDS

410 #define Rounding rounding

411 #undef Check_FLT_ROUNDS

412 #define Check_FLT_ROUNDS

413 #else

414 #define Rounding Flt_Rounds

415 #endif

416

417 #else /* ifndef IEEE_Arith */

418 #undef Check_FLT_ROUNDS

419 #undef Honor_FLT_ROUNDS

420 #undef SET_INEXACT

421 #undef Sudden_Underflow

422 #define Sudden_Underflow

423 #ifdef IBM

424 #undef Flt_Rounds

425 #define Flt_Rounds 0

426 #define Exp_shift 24

427 #define Exp_shift1 24

428 #define Exp_msk1 0x1000000

429 #define Exp_msk11 0x1000000

430 #define Exp_mask 0x7f000000

431 #define P 14

432 #define Bias 65

433 #define Exp_1 0x41000000

434 #define Exp_11 0x41000000

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

436 #define Frac_mask 0xffffff

437 #define Frac_mask1 0xffffff

438 #define Bletch 4

439 #define Ten_pmax 22

440 #define Bndry_mask 0xefffff

441 #define Bndry_mask1 0xffffff

442 #define LSB 1

443 #define Sign_bit 0x80000000

444 #define Log2P 4

445 #define Tiny0 0x100000

446 #define Tiny1 0

447 #define Quick_max 14

448 #define Int_max 15

449 #else /* VAX */

450 #undef Flt_Rounds

451 #define Flt_Rounds 1

452 #define Exp_shift 23

453 #define Exp_shift1 7

454 #define Exp_msk1 0x80

455 #define Exp_msk11 0x800000

456 #define Exp_mask 0x7f80

457 #define P 56

458 #define Bias 129

459 #define Exp_1 0x40800000

460 #define Exp_11 0x4080

461 #define Ebits 8

462 #define Frac_mask 0x7fffff

463 #define Frac_mask1 0xffff007f

464 #define Ten_pmax 24

465 #define Bletch 2

466 #define Bndry_mask 0xffff007f

467 #define Bndry_mask1 0xffff007f

468 #define LSB 0x10000

469 #define Sign_bit 0x8000

470 #define Log2P 1

471 #define Tiny0 0x80

472 #define Tiny1 0

473 #define Quick_max 15

474 #define Int_max 15

475 #endif /* IBM, VAX */

476 #endif /* IEEE_Arith */

477

478 #ifndef IEEE_Arith

479 #define ROUND_BIASED

480 #endif

481

482 #ifdef RND_PRODQUOT

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

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

485 #ifdef KR_headers

486 extern double rnd_prod(), rnd_quot();

487 #else

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

489 #endif

490 #else

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

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

493 #endif

494

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

496 #define Big1 0xffffffff

497

498 #ifndef Pack_32

499 #define Pack_32

500 #endif

501

502 #ifdef KR_headers

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

504 #else

505 #define FFFFFFFF 0xffffffffUL

506 #endif

507

508 #ifdef NO_LONG_LONG

509 #undef ULLong

510 #ifdef Just_16

511 #undef Pack_32

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

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

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

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

516 */

517 #endif

518 #else /* long long available */

519 #ifndef Llong

520 #define Llong long long

521 #endif

522 #ifndef ULLong

523 #define ULLong unsigned Llong

524 #endif

525 #endif /* NO_LONG_LONG */

526

527 #ifndef MULTIPLE_THREADS

528 #define ACQUIRE_DTOA_LOCK(n) /nothing/

529 #define FREE_DTOA_LOCK(n) /nothing/

530 #endif

531

532 #define Kmax 7

533

534 struct

535 Bigint {

536 struct Bigint *next;

537 int k, maxwds, sign, wds;

538 ULong x[1];

539 };

540

541 typedef struct Bigint Bigint;

542

543 static Bigint *freelist[Kmax+1];

544

545 static Bigint *

546 Balloc

547 #ifdef KR_headers

548 (k) int k;

549 #else

550 (int k)

551 #endif

552 {

553 int x;

554 Bigint *rv;

555 #ifndef Omit_Private_Memory

556 unsigned int len;

557 #endif

558

559 ACQUIRE_DTOA_LOCK(0);

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

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

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

563 freelist[k] = rv->next;

564 else {

565 x = 1 << k;

566 #ifdef Omit_Private_Memory

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

568 #else

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

570 /sizeof(double);

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

572 rv = (Bigint*)pmem_next;

573 pmem_next += len;

574 }

575 else

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

577 #endif

578 rv->k = k;

579 rv->maxwds = x;

580 }

581 FREE_DTOA_LOCK(0);

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

583 return rv;

584 }

585

586 static void

587 Bfree

588 #ifdef KR_headers

589 (v) Bigint *v;

590 #else

591 (Bigint *v)

592 #endif

593 {

594 if (v) {

595 if (v->k > Kmax)

596 #ifdef FREE

597 FREE((void*)v);

598 #else

599 free((void*)v);

600 #endif

601 else {

602 ACQUIRE_DTOA_LOCK(0);

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

604 freelist[v->k] = v;

605 FREE_DTOA_LOCK(0);

606 }

607 }

608 }

609

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

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

612

613 static Bigint *

614 multadd

615 #ifdef KR_headers

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

617 #else

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

619 #endif

620 {

621 int i, wds;

622 #ifdef ULLong

623 ULong *x;

624 ULLong carry, y;

625 #else

626 ULong carry, *x, y;

627 #ifdef Pack_32

628 ULong xi, z;

629 #endif

630 #endif

631 Bigint *b1;

632

633 wds = b->wds;

634 x = b->x;

635 i = 0;

636 carry = a;

637 do {

638 #ifdef ULLong

639 y = x (ULLong)m + carry;

640 carry = y >> 32;

641 *x++ = y & FFFFFFFF;

642 #else

643 #ifdef Pack_32

644 xi = *x;

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

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

647 carry = z >> 16;

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

649 #else

650 y = x m + carry;

651 carry = y >> 16;

652 *x++ = y & 0xffff;

653 #endif

654 #endif

655 }

656 while(++i < wds);

657 if (carry) {

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

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

660 Bcopy(b1, b);

661 Bfree(b);

662 b = b1;

663 }

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

665 b->wds = wds;

666 }

667 return b;

668 }

669

670 static Bigint *

671 s2b

672 #ifdef KR_headers

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

674 #else

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

676 #endif

677 {

678 Bigint *b;

679 int i, k;

680 Long x, y;

681

682 x = (nd + 8) / 9;

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

684 #ifdef Pack_32

685 b = Balloc(k);

686 b->x[0] = y9;

687 b->wds = 1;

688 #else

689 b = Balloc(k+1);

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

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

692 #endif

693

694 i = 9;

695 if (9 < nd0) {

696 s += 9;

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

698 while(++i < nd0);

699 s++;

700 }

701 else

702 s += 10;

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

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

705 return b;

706 }

707

708 static int

709 hi0bits

710 #ifdef KR_headers

711 (x) register ULong x;

712 #else

713 (register ULong x)

714 #endif

715 {

716 #ifdef PR_HAVE_BUILTIN_BITSCAN32

717 return( (!x) ? 32 : pr_bitscan_clz32(x) );

718 #else

719 register int k = 0;

720

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

722 k = 16;

723 x <<= 16;

724 }

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

726 k += 8;

727 x <<= 8;

728 }

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

730 k += 4;

731 x <<= 4;

732 }

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

734 k += 2;

735 x <<= 2;

736 }

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

738 k++;

739 if (!(x & 0x40000000))

740 return 32;

741 }

742 return k;

743 #endif /* PR_HAVE_BUILTIN_BITSCAN32 */

744 }

745

746 static int

747 lo0bits

748 #ifdef KR_headers

749 (y) ULong *y;

750 #else

751 (ULong *y)

752 #endif

753 {

754 #ifdef PR_HAVE_BUILTIN_BITSCAN32

755 int k;

756 ULong x = *y;

757

758 if (x>1)

759 *y = ( x >> (k = pr_bitscan_ctz32(x)) );

760 else

761 k = ((x ^ 1) << 5);

762 #else

763 register int k;

764 register ULong x = *y;

765

766 if (x & 7) {

767 if (x & 1)

768 return 0;

769 if (x & 2) {

770 *y = x >> 1;

771 return 1;

772 }

773 *y = x >> 2;

774 return 2;

775 }

776 k = 0;

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

778 k = 16;

779 x >>= 16;

780 }

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

782 k += 8;

783 x >>= 8;

784 }

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

786 k += 4;

787 x >>= 4;

788 }

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

790 k += 2;

791 x >>= 2;

792 }

793 if (!(x & 1)) {

794 k++;

795 x >>= 1;

796 if (!x)

797 return 32;

798 }

799 *y = x;

800 #endif /* PR_HAVE_BUILTIN_BITSCAN32 */

801 return k;

802 }

803

804 static Bigint *

805 i2b

806 #ifdef KR_headers

807 (i) int i;

808 #else

809 (int i)

810 #endif

811 {

812 Bigint *b;

813

814 b = Balloc(1);

815 b->x[0] = i;

816 b->wds = 1;

817 return b;

818 }

819

820 static Bigint *

821 mult

822 #ifdef KR_headers

823 (a, b) Bigint a, b;

824 #else

825 (Bigint a, Bigint b)

826 #endif

827 {

828 Bigint *c;

829 int k, wa, wb, wc;

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

831 ULong y;

832 #ifdef ULLong

833 ULLong carry, z;

834 #else

835 ULong carry, z;

836 #ifdef Pack_32

837 ULong z2;

838 #endif

839 #endif

840

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

842 c = a;

843 a = b;

844 b = c;

845 }

846 k = a->k;

847 wa = a->wds;

848 wb = b->wds;

849 wc = wa + wb;

850 if (wc > a->maxwds)

851 k++;

852 c = Balloc(k);

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

854 *x = 0;

855 xa = a->x;

856 xae = xa + wa;

857 xb = b->x;

858 xbe = xb + wb;

859 xc0 = c->x;

860 #ifdef ULLong

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

862 if (y = *xb++) {

863 x = xa;

864 xc = xc0;

865 carry = 0;

866 do {

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

868 carry = z >> 32;

869 *xc++ = z & FFFFFFFF;

870 }

871 while(x < xae);

872 *xc = carry;

873 }

874 }

875 #else

876 #ifdef Pack_32

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

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

879 x = xa;

880 xc = xc0;

881 carry = 0;

882 do {

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

884 carry = z >> 16;

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

886 carry = z2 >> 16;

887 Storeinc(xc, z2, z);

888 }

889 while(x < xae);

890 *xc = carry;

891 }

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

893 x = xa;

894 xc = xc0;

895 carry = 0;

896 z2 = *xc;

897 do {

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

899 carry = z >> 16;

900 Storeinc(xc, z, z2);

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

902 carry = z2 >> 16;

903 }

904 while(x < xae);

905 *xc = z2;

906 }

907 }

908 #else

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

910 if (y = *xb++) {

911 x = xa;

912 xc = xc0;

913 carry = 0;

914 do {

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

916 carry = z >> 16;

917 *xc++ = z & 0xffff;

918 }

919 while(x < xae);

920 *xc = carry;

921 }

922 }

923 #endif

924 #endif

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

926 c->wds = wc;

927 return c;

928 }

929

930 static Bigint *p5s;

931

932 static Bigint *

933 pow5mult

934 #ifdef KR_headers

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

936 #else

937 (Bigint *b, int k)

938 #endif

939 {

940 Bigint b1, p5, *p51;

941 int i;

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

943

944 if (i = k & 3)

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

946

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

948 return b;

949 if (!(p5 = p5s)) {

950 /* first time */

951 #ifdef MULTIPLE_THREADS

952 ACQUIRE_DTOA_LOCK(1);

953 if (!(p5 = p5s)) {

954 p5 = p5s = i2b(625);

955 p5->next = 0;

956 }

957 FREE_DTOA_LOCK(1);

958 #else

959 p5 = p5s = i2b(625);

960 p5->next = 0;

961 #endif

962 }

963 for(;;) {

964 if (k & 1) {

965 b1 = mult(b, p5);

966 Bfree(b);

967 b = b1;

968 }

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

970 break;

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

972 #ifdef MULTIPLE_THREADS

973 ACQUIRE_DTOA_LOCK(1);

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

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

976 p51->next = 0;

977 }

978 FREE_DTOA_LOCK(1);

979 #else

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

981 p51->next = 0;

982 #endif

983 }

984 p5 = p51;

985 }

986 return b;

987 }

988

989 static Bigint *

990 lshift

991 #ifdef KR_headers

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

993 #else

994 (Bigint *b, int k)

995 #endif

996 {

997 int i, k1, n, n1;

998 Bigint *b1;

999 ULong x, x1, *xe, z;

1000

1001 #ifdef Pack_32

1002 n = k >> 5;

1003 #else

1004 n = k >> 4;

1005 #endif

1006 k1 = b->k;

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

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

1009 k1++;

1010 b1 = Balloc(k1);

1011 x1 = b1->x;

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

1013 *x1++ = 0;

1014 x = b->x;

1015 xe = x + b->wds;

1016 #ifdef Pack_32

1017 if (k &= 0x1f) {

1018 k1 = 32 - k;

1019 z = 0;

1020 do {

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

1022 z = *x++ >> k1;

1023 }

1024 while(x < xe);

1025 if (*x1 = z)

1026 ++n1;

1027 }

1028 #else

1029 if (k &= 0xf) {

1030 k1 = 16 - k;

1031 z = 0;

1032 do {

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

1034 z = *x++ >> k1;

1035 }

1036 while(x < xe);

1037 if (*x1 = z)

1038 ++n1;

1039 }

1040 #endif

1041 else do

1042 x1++ = x++;

1043 while(x < xe);

1044 b1->wds = n1 - 1;

1045 Bfree(b);

1046 return b1;

1047 }

1048

1049 static int

1050 cmp

1051 #ifdef KR_headers

1052 (a, b) Bigint a, b;

1053 #else

1054 (Bigint a, Bigint b)

1055 #endif

1056 {

1057 ULong xa, xa0, xb, xb0;

1058 int i, j;

1059

1060 i = a->wds;

1061 j = b->wds;

1062 #ifdef DEBUG

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

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

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

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

1067 #endif

1068 if (i -= j)

1069 return i;

1070 xa0 = a->x;

1071 xa = xa0 + j;

1072 xb0 = b->x;

1073 xb = xb0 + j;

1074 for(;;) {

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

1076 return xa < xb ? -1 : 1;

1077 if (xa <= xa0)

1078 break;

1079 }

1080 return 0;

1081 }

1082

1083 static Bigint *

1084 diff

1085 #ifdef KR_headers

1086 (a, b) Bigint a, b;

1087 #else

1088 (Bigint a, Bigint b)

1089 #endif

1090 {

1091 Bigint *c;

1092 int i, wa, wb;

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

1094 #ifdef ULLong

1095 ULLong borrow, y;

1096 #else

1097 ULong borrow, y;

1098 #ifdef Pack_32

1099 ULong z;

1100 #endif

1101 #endif

1102

1103 i = cmp(a,b);

1104 if (!i) {

1105 c = Balloc(0);

1106 c->wds = 1;

1107 c->x[0] = 0;

1108 return c;

1109 }

1110 if (i < 0) {

1111 c = a;

1112 a = b;

1113 b = c;

1114 i = 1;

1115 }

1116 else

1117 i = 0;

1118 c = Balloc(a->k);

1119 c->sign = i;

1120 wa = a->wds;

1121 xa = a->x;

1122 xae = xa + wa;

1123 wb = b->wds;

1124 xb = b->x;

1125 xbe = xb + wb;

1126 xc = c->x;

1127 borrow = 0;

1128 #ifdef ULLong

1129 do {

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

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

1132 *xc++ = y & FFFFFFFF;

1133 }

1134 while(xb < xbe);

1135 while(xa < xae) {

1136 y = *xa++ - borrow;

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

1138 *xc++ = y & FFFFFFFF;

1139 }

1140 #else

1141 #ifdef Pack_32

1142 do {

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

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

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

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

1147 Storeinc(xc, z, y);

1148 }

1149 while(xb < xbe);

1150 while(xa < xae) {

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

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

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

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

1155 Storeinc(xc, z, y);

1156 }

1157 #else

1158 do {

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

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

1161 *xc++ = y & 0xffff;

1162 }

1163 while(xb < xbe);

1164 while(xa < xae) {

1165 y = *xa++ - borrow;

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

1167 *xc++ = y & 0xffff;

1168 }

1169 #endif

1170 #endif

1171 while(!*--xc)

1172 wa--;

1173 c->wds = wa;

1174 return c;

1175 }

1176

1177 static double

1178 ulp

1179 #ifdef KR_headers

1180 (dx) double dx;

1181 #else

1182 (double dx)

1183 #endif

1184 {

1185 register Long L;

1186 U x, a;

1187

1188 dval(x) = dx;

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

1190 #ifndef Avoid_Underflow

1191 #ifndef Sudden_Underflow

1192 if (L > 0) {

1193 #endif

1194 #endif

1195 #ifdef IBM

1196 L \|= Exp_msk1 >> 4;

1197 #endif

1198 word0(a) = L;

1199 word1(a) = 0;

1200 #ifndef Avoid_Underflow

1201 #ifndef Sudden_Underflow

1202 }

1203 else {

1204 L = -L >> Exp_shift;

1205 if (L < Exp_shift) {

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

1207 word1(a) = 0;

1208 }

1209 else {

1210 word0(a) = 0;

1211 L -= Exp_shift;

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

1213 }

1214 }

1215 #endif

1216 #endif

1217 return dval(a);

1218 }

1219

1220 static double

1221 b2d

1222 #ifdef KR_headers

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

1224 #else

1225 (Bigint a, int e)

1226 #endif

1227 {

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

1229 int k;

1230 U d;

1231 #ifdef VAX

1232 ULong d0, d1;

1233 #else

1234 #define d0 word0(d)

1235 #define d1 word1(d)

1236 #endif

1237

1238 xa0 = a->x;

1239 xa = xa0 + a->wds;

1240 y = *--xa;

1241 #ifdef DEBUG

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

1243 #endif

1244 k = hi0bits(y);

1245 *e = 32 - k;

1246 #ifdef Pack_32

1247 if (k < Ebits) {

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

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

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

1251 goto ret_d;

1252 }

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

1254 if (k -= Ebits) {

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

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

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

1258 }

1259 else {

1260 d0 = Exp_1 \| y;

1261 d1 = z;

1262 }

1263 #else

1264 if (k < Ebits + 16) {

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

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

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

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

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

1270 goto ret_d;

1271 }

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

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

1274 k -= Ebits + 16;

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

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

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

1278 #endif

1279 ret_d:

1280 #ifdef VAX

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

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

1283 #else

1284 #undef d0

1285 #undef d1

1286 #endif

1287 return dval(d);

1288 }

1289

1290 static Bigint *

1291 d2b

1292 #ifdef KR_headers

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

1294 #else

1295 (double dd, int e, int bits)

1296 #endif

1297 {

1298 U d;

1299 Bigint *b;

1300 int de, k;

1301 ULong *x, y, z;

1302 #ifndef Sudden_Underflow

1303 int i;

1304 #endif

1305 #ifdef VAX

1306 ULong d0, d1;

1307 #endif

1308

1309 dval(d) = dd;

1310 #ifdef VAX

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

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

1313 #else

1314 #define d0 word0(d)

1315 #define d1 word1(d)

1316 #endif

1317

1318 #ifdef Pack_32

1319 b = Balloc(1);

1320 #else

1321 b = Balloc(2);

1322 #endif

1323 x = b->x;

1324

1325 z = d0 & Frac_mask;

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

1327 #ifdef Sudden_Underflow

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

1329 #ifndef IBM

1330 z \|= Exp_msk11;

1331 #endif

1332 #else

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

1334 z \|= Exp_msk1;

1335 #endif

1336 #ifdef Pack_32

1337 if (y = d1) {

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

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

1340 z >>= k;

1341 }

1342 else

1343 x[0] = y;

1344 #ifndef Sudden_Underflow

1345 i =

1346 #endif

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

1348 }

1349 else {

1350 k = lo0bits(&z);

1351 x[0] = z;

1352 #ifndef Sudden_Underflow

1353 i =

1354 #endif

1355 b->wds = 1;

1356 k += 32;

1357 }

1358 #else

1359 if (y = d1) {

1360 if (k = lo0bits(&y))

1361 if (k >= 16) {

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

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

1364 x[2] = z >> k;

1365 i = 2;

1366 }

1367 else {

1368 x[0] = y & 0xffff;

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

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

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

1372 i = 3;

1373 }

1374 else {

1375 x[0] = y & 0xffff;

1376 x[1] = y >> 16;

1377 x[2] = z & 0xffff;

1378 x[3] = z >> 16;

1379 i = 3;

1380 }

1381 }

1382 else {

1383 #ifdef DEBUG

1384 if (!z)

1385 Bug("Zero passed to d2b");

1386 #endif

1387 k = lo0bits(&z);

1388 if (k >= 16) {

1389 x[0] = z;

1390 i = 0;

1391 }

1392 else {

1393 x[0] = z & 0xffff;

1394 x[1] = z >> 16;

1395 i = 1;

1396 }

1397 k += 32;

1398 }

1399 while(!x[i])

1400 --i;

1401 b->wds = i + 1;

1402 #endif

1403 #ifndef Sudden_Underflow

1404 if (de) {

1405 #endif

1406 #ifdef IBM

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

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

1409 #else

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

1411 *bits = P - k;

1412 #endif

1413 #ifndef Sudden_Underflow

1414 }

1415 else {

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

1417 #ifdef Pack_32

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

1419 #else

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

1421 #endif

1422 }

1423 #endif

1424 return b;

1425 }

1426 #undef d0

1427 #undef d1

1428

1429 static double

1430 ratio

1431 #ifdef KR_headers

1432 (a, b) Bigint a, b;

1433 #else

1434 (Bigint a, Bigint b)

1435 #endif

1436 {

1437 U da, db;

1438 int k, ka, kb;

1439

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

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

1442 #ifdef Pack_32

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

1444 #else

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

1446 #endif

1447 #ifdef IBM

1448 if (k > 0) {

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

1450 if (k &= 3)

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

1452 }

1453 else {

1454 k = -k;

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

1456 if (k &= 3)

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

1458 }

1459 #else

1460 if (k > 0)

1461 word0(da) += k*Exp_msk1;

1462 else {

1463 k = -k;

1464 word0(db) += k*Exp_msk1;

1465 }

1466 #endif

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

1468 }

1469

1470 static CONST double

1471 tens[] = {

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

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

1474 1e20, 1e21, 1e22

1475 #ifdef VAX

1476 , 1e23, 1e24

1477 #endif

1478 };

1479

1480 static CONST double

1481 #ifdef IEEE_Arith

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

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

1484 #ifdef Avoid_Underflow

1485 9007199254740992.*9007199254740992.e-256

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

1487 #else

1488 1e-256

1489 #endif

1490 };

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

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

1493 #define Scale_Bit 0x10

1494 #define n_bigtens 5

1495 #else

1496 #ifdef IBM

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

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

1499 #define n_bigtens 3

1500 #else

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

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

1503 #define n_bigtens 2

1504 #endif

1505 #endif

1506

1507 #ifndef IEEE_Arith

1508 #undef INFNAN_CHECK

1509 #endif

1510

1511 #ifdef INFNAN_CHECK

1512

1513 #ifndef NAN_WORD0

1514 #define NAN_WORD0 0x7ff80000

1515 #endif

1516

1517 #ifndef NAN_WORD1

1518 #define NAN_WORD1 0

1519 #endif

1520

1521 static int

1522 match

1523 #ifdef KR_headers

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

1525 #else

1526 (CONST char *sp, char t)

1527 #endif

1528 {

1529 int c, d;

1530 CONST char s = sp;

1531

1532 while(d = *t++) {

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

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

1535 if (c != d)

1536 return 0;

1537 }

1538 *sp = s + 1;

1539 return 1;

1540 }

1541

1542 #ifndef No_Hex_NaN

1543 static void

1544 hexnan

1545 #ifdef KR_headers

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

1547 #else

1548 (double rvp, CONST char *sp)

1549 #endif

1550 {

1551 ULong c, x[2];

1552 CONST char *s;

1553 int havedig, udx0, xshift;

1554

1555 x[0] = x[1] = 0;

1556 havedig = xshift = 0;

1557 udx0 = 1;

1558 s = *sp;

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

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

1561 c -= '0';

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

1563 c += 10 - 'a';

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

1565 c += 10 - 'A';

1566 else if (c <= ' ') {

1567 if (udx0 && havedig) {

1568 udx0 = 0;

1569 xshift = 1;

1570 }

1571 continue;

1572 }

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

1574 *sp = s + 1;

1575 break;

1576 }

1577 else

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

1579 havedig = 1;

1580 if (xshift) {

1581 xshift = 0;

1582 x[0] = x[1];

1583 x[1] = 0;

1584 }

1585 if (udx0)

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

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

1588 }

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

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

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

1592 }

1593 }

1594 #endif /No_Hex_NaN/

1595 #endif /* INFNAN_CHECK */

1596

1597 PR_IMPLEMENT(double)

1598 PR_strtod

1599 #ifdef KR_headers

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

1601 #else

1602 (CONST char s00, char *se)

1603 #endif

1604 {

1605 #ifdef Avoid_Underflow

1606 int scale;

1607 #endif

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

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

1610 CONST char s, s0, *s1;

1611 double aadj, aadj1, adj;

1612 U aadj2, rv, rv0;

1613 Long L;

1614 ULong y, z;

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

1616 #ifdef SET_INEXACT

1617 int inexact, oldinexact;

1618 #endif

1619 #ifdef Honor_FLT_ROUNDS

1620 int rounding;

1621 #endif

1622 #ifdef USE_LOCALE

1623 CONST char *s2;

1624 #endif

1625

1626 if (!_pr_initialized) _PR_ImplicitInitialization();

1627

1628 sign = nz0 = nz = 0;

1629 dval(rv) = 0.;

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

1631 case '-':

1632 sign = 1;

1633 /* no break */

1634 case '+':

1635 if (*++s)

1636 goto break2;

1637 /* no break */

1638 case 0:

1639 goto ret0;

1640 case '\t':

1641 case '\n':

1642 case '\v':

1643 case '\f':

1644 case '\r':

1645 case ' ':

1646 continue;

1647 default:

1648 goto break2;

1649 }

1650 break2:

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

1652 nz0 = 1;

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

1654 if (!*s)

1655 goto ret;

1656 }

1657 s0 = s;

1658 y = z = 0;

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

1660 if (nd < 9)

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

1662 else if (nd < 16)

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

1664 nd0 = nd;

1665 #ifdef USE_LOCALE

1666 s1 = localeconv()->decimal_point;

1667 if (c == *s1) {

1668 c = '.';

1669 if (*++s1) {

1670 s2 = s;

1671 for(;;) {

1672 if (++s2 != s1) {

1673 c = 0;

1674 break;

1675 }

1676 if (!*++s1) {

1677 s = s2;

1678 break;

1679 }

1680 }

1681 }

1682 }

1683 #endif

1684 if (c == '.') {

1685 c = *++s;

1686 if (!nd) {

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

1688 nz++;

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

1690 s0 = s;

1691 nf += nz;

1692 nz = 0;

1693 goto have_dig;

1694 }

1695 goto dig_done;

1696 }

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

1698 have_dig:

1699 nz++;

1700 if (c -= '0') {

1701 nf += nz;

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

1703 if (nd++ < 9)

1704 y *= 10;

1705 else if (nd <= DBL_DIG + 1)

1706 z *= 10;

1707 if (nd++ < 9)

1708 y = 10*y + c;

1709 else if (nd <= DBL_DIG + 1)

1710 z = 10*z + c;

1711 nz = 0;

1712 }

1713 }

1714 }

1715 dig_done:

1716 if (nd > 64 * 1024)

1717 goto ret0;

1718 e = 0;

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

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

1721 goto ret0;

1722 }

1723 s00 = s;

1724 esign = 0;

1725 switch(c = *++s) {

1726 case '-':

1727 esign = 1;

1728 case '+':

1729 c = *++s;

1730 }

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

1732 while(c == '0')

1733 c = *++s;

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

1735 L = c - '0';

1736 s1 = s;

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

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

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

1740 /* Avoid confusion from exponents

1741 * so large that e might overflow.

1742 */

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

1744 else

1745 e = (int)L;

1746 if (esign)

1747 e = -e;

1748 }

1749 else

1750 e = 0;

1751 }

1752 else

1753 s = s00;

1754 }

1755 if (!nd) {

1756 if (!nz && !nz0) {

1757 #ifdef INFNAN_CHECK

1758 /* Check for Nan and Infinity */

1759 switch(c) {

1760 case 'i':

1761 case 'I':

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

1763 --s;

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

1765 ++s;

1766 word0(rv) = 0x7ff00000;

1767 word1(rv) = 0;

1768 goto ret;

1769 }

1770 break;

1771 case 'n':

1772 case 'N':

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

1774 word0(rv) = NAN_WORD0;

1775 word1(rv) = NAN_WORD1;

1776 #ifndef No_Hex_NaN

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

1778 hexnan(&rv, &s);

1779 #endif

1780 goto ret;

1781 }

1782 }

1783 #endif /* INFNAN_CHECK */

1784 ret0:

1785 s = s00;

1786 sign = 0;

1787 }

1788 goto ret;

1789 }

1790 e1 = e -= nf;

1791

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

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

1794 * after is the integer represented by those digits times

1795 * 10*e /

1796

1797 if (!nd0)

1798 nd0 = nd;

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

1800 dval(rv) = y;

1801 if (k > 9) {

1802 #ifdef SET_INEXACT

1803 if (k > DBL_DIG)

1804 oldinexact = get_inexact();

1805 #endif

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

1807 }

1808 bd0 = 0;

1809 if (nd <= DBL_DIG

1810 #ifndef RND_PRODQUOT

1811 #ifndef Honor_FLT_ROUNDS

1812 && Flt_Rounds == 1

1813 #endif

1814 #endif

1815 ) {

1816 if (!e)

1817 goto ret;

1818 if (e > 0) {

1819 if (e <= Ten_pmax) {

1820 #ifdef VAX

1821 goto vax_ovfl_check;

1822 #else

1823 #ifdef Honor_FLT_ROUNDS

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

1825 if (sign) {

1826 rv = -rv;

1827 sign = 0;

1828 }

1829 #endif

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

1831 goto ret;

1832 #endif

1833 }

1834 i = DBL_DIG - nd;

1835 if (e <= Ten_pmax + i) {

1836 /* A fancier test would sometimes let us do

1837 * this for larger i values.

1838 */

1839 #ifdef Honor_FLT_ROUNDS

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

1841 if (sign) {

1842 rv = -rv;

1843 sign = 0;

1844 }

1845 #endif

1846 e -= i;

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

1848 #ifdef VAX

1849 /* VAX exponent range is so narrow we must

1850 * worry about overflow here...

1851 */

1852 vax_ovfl_check:

1853 word0(rv) -= P*Exp_msk1;

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

1855 if ((word0(rv) & Exp_mask)

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

1857 goto ovfl;

1858 word0(rv) += P*Exp_msk1;

1859 #else

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

1861 #endif

1862 goto ret;

1863 }

1864 }

1865 #ifndef Inaccurate_Divide

1866 else if (e >= -Ten_pmax) {

1867 #ifdef Honor_FLT_ROUNDS

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

1869 if (sign) {

1870 rv = -rv;

1871 sign = 0;

1872 }

1873 #endif

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

1875 goto ret;

1876 }

1877 #endif

1878 }

1879 e1 += nd - k;

1880

1881 #ifdef IEEE_Arith

1882 #ifdef SET_INEXACT

1883 inexact = 1;

1884 if (k <= DBL_DIG)

1885 oldinexact = get_inexact();

1886 #endif

1887 #ifdef Avoid_Underflow

1888 scale = 0;

1889 #endif

1890 #ifdef Honor_FLT_ROUNDS

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

1892 if (sign)

1893 rounding = rounding == 2 ? 0 : 2;

1894 else

1895 if (rounding != 2)

1896 rounding = 0;

1897 }

1898 #endif

1899 #endif /IEEE_Arith/

1900

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

1902

1903 if (e1 > 0) {

1904 if (i = e1 & 15)

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

1906 if (e1 &= ~15) {

1907 if (e1 > DBL_MAX_10_EXP) {

1908 ovfl:

1909 #ifndef NO_ERRNO

1910 PR_SetError(PR_RANGE_ERROR, 0);

1911 #endif

1912 /* Can't trust HUGE_VAL */

1913 #ifdef IEEE_Arith

1914 #ifdef Honor_FLT_ROUNDS

1915 switch(rounding) {

1916 case 0: /* toward 0 */

1917 case 3: /* toward -infinity */

1918 word0(rv) = Big0;

1919 word1(rv) = Big1;

1920 break;

1921 default:

1922 word0(rv) = Exp_mask;

1923 word1(rv) = 0;

1924 }

1925 #else /Honor_FLT_ROUNDS/

1926 word0(rv) = Exp_mask;

1927 word1(rv) = 0;

1928 #endif /Honor_FLT_ROUNDS/

1929 #ifdef SET_INEXACT

1930 /* set overflow bit */

1931 dval(rv0) = 1e300;

1932 dval(rv0) *= dval(rv0);

1933 #endif

1934 #else /IEEE_Arith/

1935 word0(rv) = Big0;

1936 word1(rv) = Big1;

1937 #endif /IEEE_Arith/

1938 if (bd0)

1939 goto retfree;

1940 goto ret;

1941 }

1942 e1 >>= 4;

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

1944 if (e1 & 1)

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

1946 /* The last multiplication could overflow. */

1947 word0(rv) -= P*Exp_msk1;

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

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

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

1951 goto ovfl;

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

1953 /* set to largest number */

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

1955 word0(rv) = Big0;

1956 word1(rv) = Big1;

1957 }

1958 else

1959 word0(rv) += P*Exp_msk1;

1960 }

1961 }

1962 else if (e1 < 0) {

1963 e1 = -e1;

1964 if (i = e1 & 15)

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

1966 if (e1 >>= 4) {

1967 if (e1 >= 1 << n_bigtens)

1968 goto undfl;

1969 #ifdef Avoid_Underflow

1970 if (e1 & Scale_Bit)

1971 scale = 2*P;

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

1973 if (e1 & 1)

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

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

1976 >> Exp_shift)) > 0) {

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

1978 if (j >= 32) {

1979 word1(rv) = 0;

1980 if (j >= 53)

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

1982 else

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

1984 }

1985 else

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

1987 }

1988 #else

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

1990 if (e1 & 1)

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

1992 /* The last multiplication could underflow. */

1993 dval(rv0) = dval(rv);

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

1995 if (!dval(rv)) {

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

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

1998 #endif

1999 if (!dval(rv)) {

2000 undfl:

2001 dval(rv) = 0.;

2002 #ifndef NO_ERRNO

2003 PR_SetError(PR_RANGE_ERROR, 0);

2004 #endif

2005 if (bd0)

2006 goto retfree;

2007 goto ret;

2008 }

2009 #ifndef Avoid_Underflow

2010 word0(rv) = Tiny0;

2011 word1(rv) = Tiny1;

2012 /* The refinement below will clean

2013 * this approximation up.

2014 */

2015 }

2016 #endif

2017 }

2018 }

2019

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

2021

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

2023

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

2025

2026 for(;;) {

2027 bd = Balloc(bd0->k);

2028 Bcopy(bd, bd0);

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

2030 bs = i2b(1);

2031

2032 if (e >= 0) {

2033 bb2 = bb5 = 0;

2034 bd2 = bd5 = e;

2035 }

2036 else {

2037 bb2 = bb5 = -e;

2038 bd2 = bd5 = 0;

2039 }

2040 if (bbe >= 0)

2041 bb2 += bbe;

2042 else

2043 bd2 -= bbe;

2044 bs2 = bb2;

2045 #ifdef Honor_FLT_ROUNDS

2046 if (rounding != 1)

2047 bs2++;

2048 #endif

2049 #ifdef Avoid_Underflow

2050 j = bbe - scale;

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

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

2053 j += P - Emin;

2054 else

2055 j = P + 1 - bbbits;

2056 #else /Avoid_Underflow/

2057 #ifdef Sudden_Underflow

2058 #ifdef IBM

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

2060 #else

2061 j = P + 1 - bbbits;

2062 #endif

2063 #else /Sudden_Underflow/

2064 j = bbe;

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

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

2067 j += P - Emin;

2068 else

2069 j = P + 1 - bbbits;

2070 #endif /Sudden_Underflow/

2071 #endif /Avoid_Underflow/

2072 bb2 += j;

2073 bd2 += j;

2074 #ifdef Avoid_Underflow

2075 bd2 += scale;

2076 #endif

2077 i = bb2 < bd2 ? bb2 : bd2;

2078 if (i > bs2)

2079 i = bs2;

2080 if (i > 0) {

2081 bb2 -= i;

2082 bd2 -= i;

2083 bs2 -= i;

2084 }

2085 if (bb5 > 0) {

2086 bs = pow5mult(bs, bb5);

2087 bb1 = mult(bs, bb);

2088 Bfree(bb);

2089 bb = bb1;

2090 }

2091 if (bb2 > 0)

2092 bb = lshift(bb, bb2);

2093 if (bd5 > 0)

2094 bd = pow5mult(bd, bd5);

2095 if (bd2 > 0)

2096 bd = lshift(bd, bd2);

2097 if (bs2 > 0)

2098 bs = lshift(bs, bs2);

2099 delta = diff(bb, bd);

2100 dsign = delta->sign;

2101 delta->sign = 0;

2102 i = cmp(delta, bs);

2103 #ifdef Honor_FLT_ROUNDS

2104 if (rounding != 1) {

2105 if (i < 0) {

2106 /* Error is less than an ulp */

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

2108 /* exact */

2109 #ifdef SET_INEXACT

2110 inexact = 0;

2111 #endif

2112 break;

2113 }

2114 if (rounding) {

2115 if (dsign) {

2116 adj = 1.;

2117 goto apply_adj;

2118 }

2119 }

2120 else if (!dsign) {

2121 adj = -1.;

2122 if (!word1(rv)

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

2124 y = word0(rv) & Exp_mask;

2125 #ifdef Avoid_Underflow

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

2127 #else

2128 if (y)

2129 #endif

2130 {

2131 delta = lshift(delta,Log2P);

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

2133 adj = -0.5;

2134 }

2135 }

2136 apply_adj:

2137 #ifdef Avoid_Underflow

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

2139 <= 2PExp_msk1)

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

2141 #else

2142 #ifdef Sudden_Underflow

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

2144 P*Exp_msk1) {

2145 word0(rv) += P*Exp_msk1;

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

2147 word0(rv) -= P*Exp_msk1;

2148 }

2149 else

2150 #endif /Sudden_Underflow/

2151 #endif /Avoid_Underflow/

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

2153 }

2154 break;

2155 }

2156 adj = ratio(delta, bs);

2157 if (adj < 1.)

2158 adj = 1.;

2159 if (adj <= 0x7ffffffe) {

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

2161 y = adj;

2162 if (y != adj) {

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

2164 y++;

2165 adj = y;

2166 }

2167 }

2168 #ifdef Avoid_Underflow

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

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

2171 #else

2172 #ifdef Sudden_Underflow

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

2174 word0(rv) += P*Exp_msk1;

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

2176 if (dsign)

2177 dval(rv) += adj;

2178 else

2179 dval(rv) -= adj;

2180 word0(rv) -= P*Exp_msk1;

2181 goto cont;

2182 }

2183 #endif /Sudden_Underflow/

2184 #endif /Avoid_Underflow/

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

2186 if (dsign)

2187 dval(rv) += adj;

2188 else

2189 dval(rv) -= adj;

2190 goto cont;

2191 }

2192 #endif /Honor_FLT_ROUNDS/

2193

2194 if (i < 0) {

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

2196 * special case of mantissa a power of two.

2197 */

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

2199 #ifdef IEEE_Arith

2200 #ifdef Avoid_Underflow

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

2202 #else

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

2204 #endif

2205 #endif

2206 ) {

2207 #ifdef SET_INEXACT

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

2209 inexact = 0;

2210 #endif

2211 break;

2212 }

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

2214 /* exact result */

2215 #ifdef SET_INEXACT

2216 inexact = 0;

2217 #endif

2218 break;

2219 }

2220 delta = lshift(delta,Log2P);

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

2222 goto drop_down;

2223 break;

2224 }

2225 if (i == 0) {

2226 /* exactly half-way between */

2227 if (dsign) {

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

2229 && word1(rv) == (

2230 #ifdef Avoid_Underflow

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

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

2233 #endif

2234 0xffffffff)) {

2235 /boundary case -- increment exponent/

2236 word0(rv) = (word0(rv) & Exp_mask)

2237 + Exp_msk1

2238 #ifdef IBM

2239 \| Exp_msk1 >> 4

2240 #endif

2241 ;

2242 word1(rv) = 0;

2243 #ifdef Avoid_Underflow

2244 dsign = 0;

2245 #endif

2246 break;

2247 }

2248 }

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

2250 drop_down:

2251 /* boundary case -- decrement exponent */

2252 #ifdef Sudden_Underflow /{{/

2253 L = word0(rv) & Exp_mask;

2254 #ifdef IBM

2255 if (L < Exp_msk1)

2256 #else

2257 #ifdef Avoid_Underflow

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

2259 #else

2260 if (L <= Exp_msk1)

2261 #endif /Avoid_Underflow/

2262 #endif /IBM/

2263 goto undfl;

2264 L -= Exp_msk1;

2265 #else /Sudden_Underflow}{/

2266 #ifdef Avoid_Underflow

2267 if (scale) {

2268 L = word0(rv) & Exp_mask;

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

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

2271 /* round even ==> */

2272 /* accept rv */

2273 break;

2274 /* rv = smallest denormal */

2275 goto undfl;

2276 }

2277 }

2278 #endif /Avoid_Underflow/

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

2280 #endif /Sudden_Underflow}}/

2281 word0(rv) = L \| Bndry_mask1;

2282 word1(rv) = 0xffffffff;

2283 #ifdef IBM

2284 goto cont;

2285 #else

2286 break;

2287 #endif

2288 }

2289 #ifndef ROUND_BIASED

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

2291 break;

2292 #endif

2293 if (dsign)

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

2295 #ifndef ROUND_BIASED

2296 else {

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

2298 #ifndef Sudden_Underflow

2299 if (!dval(rv))

2300 goto undfl;

2301 #endif

2302 }

2303 #ifdef Avoid_Underflow

2304 dsign = 1 - dsign;

2305 #endif

2306 #endif

2307 break;

2308 }

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

2310 if (dsign)

2311 aadj = aadj1 = 1.;

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

2313 #ifndef Sudden_Underflow

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

2315 goto undfl;

2316 #endif

2317 aadj = 1.;

2318 aadj1 = -1.;

2319 }

2320 else {

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

2322 /* rounded down... */

2323

2324 if (aadj < 2./FLT_RADIX)

2325 aadj = 1./FLT_RADIX;

2326 else

2327 aadj *= 0.5;

2328 aadj1 = -aadj;

2329 }

2330 }

2331 else {

2332 aadj *= 0.5;

2333 aadj1 = dsign ? aadj : -aadj;

2334 #ifdef Check_FLT_ROUNDS

2335 switch(Rounding) {

2336 case 2: /* towards +infinity */

2337 aadj1 -= 0.5;

2338 break;

2339 case 0: /* towards 0 */

2340 case 3: /* towards -infinity */

2341 aadj1 += 0.5;

2342 }

2343 #else

2344 if (Flt_Rounds == 0)

2345 aadj1 += 0.5;

2346 #endif /Check_FLT_ROUNDS/

2347 }

2348 y = word0(rv) & Exp_mask;

2349

2350 /* Check for overflow */

2351

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

2353 dval(rv0) = dval(rv);

2354 word0(rv) -= P*Exp_msk1;

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

2356 dval(rv) += adj;

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

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

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

2360 goto ovfl;

2361 word0(rv) = Big0;

2362 word1(rv) = Big1;

2363 goto cont;

2364 }

2365 else

2366 word0(rv) += P*Exp_msk1;

2367 }

2368 else {

2369 #ifdef Avoid_Underflow

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

2371 if (aadj <= 0x7fffffff) {

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

2373 z = 1;

2374 aadj = z;

2375 aadj1 = dsign ? aadj : -aadj;

2376 }

2377 dval(aadj2) = aadj1;

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

2379 aadj1 = dval(aadj2);

2380 }

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

2382 dval(rv) += adj;

2383 #else

2384 #ifdef Sudden_Underflow

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

2386 dval(rv0) = dval(rv);

2387 word0(rv) += P*Exp_msk1;

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

2389 dval(rv) += adj;

2390 #ifdef IBM

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

2392 #else

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

2394 #endif

2395 {

2396 if (word0(rv0) == Tiny0

2397 && word1(rv0) == Tiny1)

2398 goto undfl;

2399 word0(rv) = Tiny0;

2400 word1(rv) = Tiny1;

2401 goto cont;

2402 }

2403 else

2404 word0(rv) -= P*Exp_msk1;

2405 }

2406 else {

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

2408 dval(rv) += adj;

2409 }

2410 #else /Sudden_Underflow/

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

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

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

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

2415 * trouble from bits lost to denormalization;

2416 * example: 1.2e-307 .

2417 */

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

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

2420 if (!dsign)

2421 aadj1 = -aadj1;

2422 }

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

2424 dval(rv) += adj;

2425 #endif /Sudden_Underflow/

2426 #endif /Avoid_Underflow/

2427 }

2428 z = word0(rv) & Exp_mask;

2429 #ifndef SET_INEXACT

2430 #ifdef Avoid_Underflow

2431 if (!scale)

2432 #endif

2433 if (y == z) {

2434 /* Can we stop now? */

2435 L = (Long)aadj;

2436 aadj -= L;

2437 /* The tolerances below are conservative. */

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

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

2440 break;

2441 }

2442 else if (aadj < .4999999/FLT_RADIX)

2443 break;

2444 }

2445 #endif

2446 cont:

2447 Bfree(bb);

2448 Bfree(bd);

2449 Bfree(bs);

2450 Bfree(delta);

2451 }

2452 #ifdef SET_INEXACT

2453 if (inexact) {

2454 if (!oldinexact) {

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

2456 word1(rv0) = 0;

2457 dval(rv0) += 1.;

2458 }

2459 }

2460 else if (!oldinexact)

2461 clear_inexact();

2462 #endif

2463 #ifdef Avoid_Underflow

2464 if (scale) {

2465 word0(rv0) = Exp_1 - 2PExp_msk1;

2466 word1(rv0) = 0;

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

2468 #ifndef NO_ERRNO

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

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

2471 PR_SetError(PR_RANGE_ERROR, 0);

2472 #endif

2473 }

2474 #endif /* Avoid_Underflow */

2475 #ifdef SET_INEXACT

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

2477 /* set underflow bit */

2478 dval(rv0) = 1e-300;

2479 dval(rv0) *= dval(rv0);

2480 }

2481 #endif

2482 retfree:

2483 Bfree(bb);

2484 Bfree(bd);

2485 Bfree(bs);

2486 Bfree(bd0);

2487 Bfree(delta);

2488 ret:

2489 if (se)

2490 se = (char )s;

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

2492 }

2493

2494 static int

2495 quorem

2496 #ifdef KR_headers

2497 (b, S) Bigint b, S;

2498 #else

2499 (Bigint b, Bigint S)

2500 #endif

2501 {

2502 int n;

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

2504 #ifdef ULLong

2505 ULLong borrow, carry, y, ys;

2506 #else

2507 ULong borrow, carry, y, ys;

2508 #ifdef Pack_32

2509 ULong si, z, zs;

2510 #endif

2511 #endif

2512

2513 n = S->wds;

2514 #ifdef DEBUG

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

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

2517 #endif

2518 if (b->wds < n)

2519 return 0;

2520 sx = S->x;

2521 sxe = sx + --n;

2522 bx = b->x;

2523 bxe = bx + n;

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

2525 #ifdef DEBUG

2526 /debug/ if (q > 9)

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

2528 #endif

2529 if (q) {

2530 borrow = 0;

2531 carry = 0;

2532 do {

2533 #ifdef ULLong

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

2535 carry = ys >> 32;

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

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

2538 *bx++ = y & FFFFFFFF;

2539 #else

2540 #ifdef Pack_32

2541 si = *sx++;

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

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

2544 carry = zs >> 16;

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

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

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

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

2549 Storeinc(bx, z, y);

2550 #else

2551 ys = sx++ q + carry;

2552 carry = ys >> 16;

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

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

2555 *bx++ = y & 0xffff;

2556 #endif

2557 #endif

2558 }

2559 while(sx <= sxe);

2560 if (!*bxe) {

2561 bx = b->x;

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

2563 --n;

2564 b->wds = n;

2565 }

2566 }

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

2568 q++;

2569 borrow = 0;

2570 carry = 0;

2571 bx = b->x;

2572 sx = S->x;

2573 do {

2574 #ifdef ULLong

2575 ys = *sx++ + carry;

2576 carry = ys >> 32;

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

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

2579 *bx++ = y & FFFFFFFF;

2580 #else

2581 #ifdef Pack_32

2582 si = *sx++;

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

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

2585 carry = zs >> 16;

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

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

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

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

2590 Storeinc(bx, z, y);

2591 #else

2592 ys = *sx++ + carry;

2593 carry = ys >> 16;

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

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

2596 *bx++ = y & 0xffff;

2597 #endif

2598 #endif

2599 }

2600 while(sx <= sxe);

2601 bx = b->x;

2602 bxe = bx + n;

2603 if (!*bxe) {

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

2605 --n;

2606 b->wds = n;

2607 }

2608 }

2609 return q;

2610 }

2611

2612 #ifndef MULTIPLE_THREADS

2613 static char *dtoa_result;

2614 #endif

2615

2616 static char *

2617 #ifdef KR_headers

2618 rv_alloc(i) int i;

2619 #else

2620 rv_alloc(int i)

2621 #endif

2622 {

2623 int j, k, *r;

2624

2625 j = sizeof(ULong);

2626 for(k = 0;

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

2628 j <<= 1)

2629 k++;

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

2631 *r = k;

2632 return

2633 #ifndef MULTIPLE_THREADS

2634 dtoa_result =

2635 #endif

2636 (char *)(r+1);

2637 }

2638

2639 static char *

2640 #ifdef KR_headers

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

2642 #else

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

2644 #endif

2645 {

2646 char rv, t;

2647

2648 t = rv = rv_alloc(n);

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

2650 if (rve)

2651 *rve = t;

2652 return rv;

2653 }

2654

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

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

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

2658 * when MULTIPLE_THREADS is not defined.

2659 */

2660

2661 static void

2662 #ifdef KR_headers

2663 freedtoa(s) char *s;

2664 #else

2665 freedtoa(char *s)

2666 #endif

2667 {

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

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

2670 Bfree(b);

2671 #ifndef MULTIPLE_THREADS

2672 if (s == dtoa_result)

2673 dtoa_result = 0;

2674 #endif

2675 }

2676

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

2678 *

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

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

2681 *

2682 * Modifications:

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

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

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

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

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

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

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

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

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

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

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

2694 * as would satisfaction of the stopping test with strict

2695 * inequality.

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

2697 * quantities.

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

2699 * we use floating-point arithmetic rather than resorting

2700 * to multiple-precision integers.

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

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

2703 * multiple-precision integer arithmetic only if we cannot

2704 * guarantee that the floating-point calculation has given

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

2706 * "uniformly" distributed input, the probability is

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

2708 * calculation.

2709 */

2710

2711 static char *

2712 dtoa

2713 #ifdef KR_headers

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

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

2716 #else

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

2718 #endif

2719 {

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

2721 of ecvt and fcvt; trailing zeros are suppressed from

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

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

2724 then *decpt is set to 9999.

2725

2726 mode:

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

2728 and rounded to nearest.

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

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

2731 1e23 whereas mode 1 gives 9.999999999999999e22.

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

2733 return value similar to that of ecvt, except

2734 that trailing zeros are suppressed.

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

2736 gives a return value similar to that from fcvt,

2737 except that trailing zeros are suppressed, and

2738 ndigits can be negative.

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

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

2741 possibly return a shorter string that rounds to d.

2742 With IEEE arithmetic and compilation with

2743 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same

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

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

2746 fast floating-point estimate (if applicable).

2747

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

2749

2750 Sufficient space is allocated to the return value

2751 to hold the suppressed trailing zeros.

2752 */

2753

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

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

2756 spec_case, try_quick;

2757 Long L;

2758 #ifndef Sudden_Underflow

2759 int denorm;

2760 ULong x;

2761 #endif

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

2763 U d, d2, eps;

2764 double ds;

2765 char s, s0;

2766 #ifdef Honor_FLT_ROUNDS

2767 int rounding;

2768 #endif

2769 #ifdef SET_INEXACT

2770 int inexact, oldinexact;

2771 #endif

2772

2773 #ifndef MULTIPLE_THREADS

2774 if (dtoa_result) {

2775 freedtoa(dtoa_result);

2776 dtoa_result = 0;

2777 }

2778 #endif

2779

2780 dval(d) = dd;

2781 if (word0(d) & Sign_bit) {

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

2783 *sign = 1;

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

2785 }

2786 else

2787 *sign = 0;

2788

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

2790 #ifdef IEEE_Arith

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

2792 #else

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

2794 #endif

2795 {

2796 /* Infinity or NaN */

2797 *decpt = 9999;

2798 #ifdef IEEE_Arith

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

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

2801 #endif

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

2803 }

2804 #endif

2805 #ifdef IBM

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

2807 #endif

2808 if (!dval(d)) {

2809 *decpt = 1;

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

2811 }

2812

2813 #ifdef SET_INEXACT

2814 try_quick = oldinexact = get_inexact();

2815 inexact = 1;

2816 #endif

2817 #ifdef Honor_FLT_ROUNDS

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

2819 if (*sign)

2820 rounding = rounding == 2 ? 0 : 2;

2821 else

2822 if (rounding != 2)

2823 rounding = 0;

2824 }

2825 #endif

2826

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

2828 #ifdef Sudden_Underflow

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

2830 #else

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

2832 #endif

2833 dval(d2) = dval(d);

2834 word0(d2) &= Frac_mask1;

2835 word0(d2) \|= Exp_11;

2836 #ifdef IBM

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

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

2839 #endif

2840

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

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

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

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

2845 *

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

2847 *

2848 * k = (i - Bias)*0.301029995663981

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

2850 *

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

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

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

2854 * to compensate for any error in the multiplication of

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

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

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

2858 * Hence we adjust the constant term to 0.1760912590558.

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

2860 * but this is probably not worthwhile.)

2861 */

2862

2863 i -= Bias;

2864 #ifdef IBM

2865 i <<= 2;

2866 i += j;

2867 #endif

2868 #ifndef Sudden_Underflow

2869 denorm = 0;

2870 }

2871 else {

2872 /* d is denormalized */

2873

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

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

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

2877 dval(d2) = x;

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

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

2880 denorm = 1;

2881 }

2882 #endif

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

2884 k = (int)ds;

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

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

2887 k_check = 1;

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

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

2890 k--;

2891 k_check = 0;

2892 }

2893 j = bbits - i - 1;

2894 if (j >= 0) {

2895 b2 = 0;

2896 s2 = j;

2897 }

2898 else {

2899 b2 = -j;

2900 s2 = 0;

2901 }

2902 if (k >= 0) {

2903 b5 = 0;

2904 s5 = k;

2905 s2 += k;

2906 }

2907 else {

2908 b2 -= k;

2909 b5 = -k;

2910 s5 = 0;

2911 }

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

2913 mode = 0;

2914

2915 #ifndef SET_INEXACT

2916 #ifdef Check_FLT_ROUNDS

2917 try_quick = Rounding == 1;

2918 #else

2919 try_quick = 1;

2920 #endif

2921 #endif /SET_INEXACT/

2922

2923 if (mode > 5) {

2924 mode -= 4;

2925 try_quick = 0;

2926 }

2927 leftright = 1;

2928 switch(mode) {

2929 case 0:

2930 case 1:

2931 ilim = ilim1 = -1;

2932 i = 18;

2933 ndigits = 0;

2934 break;

2935 case 2:

2936 leftright = 0;

2937 /* no break */

2938 case 4:

2939 if (ndigits <= 0)

2940 ndigits = 1;

2941 ilim = ilim1 = i = ndigits;

2942 break;

2943 case 3:

2944 leftright = 0;

2945 /* no break */

2946 case 5:

2947 i = ndigits + k + 1;

2948 ilim = i;

2949 ilim1 = i - 1;

2950 if (i <= 0)

2951 i = 1;

2952 }

2953 s = s0 = rv_alloc(i);

2954

2955 #ifdef Honor_FLT_ROUNDS

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

2957 leftright = 0;

2958 #endif

2959

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

2961

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

2963

2964 i = 0;

2965 dval(d2) = dval(d);

2966 k0 = k;

2967 ilim0 = ilim;

2968 ieps = 2; /* conservative */

2969 if (k > 0) {

2970 ds = tens[k&0xf];

2971 j = k >> 4;

2972 if (j & Bletch) {

2973 /* prevent overflows */

2974 j &= Bletch - 1;

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

2976 ieps++;

2977 }

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

2979 if (j & 1) {

2980 ieps++;

2981 ds *= bigtens[i];

2982 }

2983 dval(d) /= ds;

2984 }

2985 else if (j1 = -k) {

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

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

2988 if (j & 1) {

2989 ieps++;

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

2991 }

2992 }

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

2994 if (ilim1 <= 0)

2995 goto fast_failed;

2996 ilim = ilim1;

2997 k--;

2998 dval(d) *= 10.;

2999 ieps++;

3000 }

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

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

3003 if (ilim == 0) {

3004 S = mhi = 0;

3005 dval(d) -= 5.;

3006 if (dval(d) > dval(eps))

3007 goto one_digit;

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

3009 goto no_digits;

3010 goto fast_failed;

3011 }

3012 #ifndef No_leftright

3013 if (leftright) {

3014 /* Use Steele & White method of only

3015 * generating digits needed.

3016 */

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

3018 for(i = 0;;) {

3019 L = dval(d);

3020 dval(d) -= L;

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

3022 if (dval(d) < dval(eps))

3023 goto ret1;

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

3025 goto bump_up;

3026 if (++i >= ilim)

3027 break;

3028 dval(eps) *= 10.;

3029 dval(d) *= 10.;

3030 }

3031 }

3032 else {

3033 #endif

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

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

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

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

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

3039 ilim = i;

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

3041 if (i == ilim) {

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

3043 goto bump_up;

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

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

3046 s++;

3047 goto ret1;

3048 }

3049 break;

3050 }

3051 }

3052 #ifndef No_leftright

3053 }

3054 #endif

3055 fast_failed:

3056 s = s0;

3057 dval(d) = dval(d2);

3058 k = k0;

3059 ilim = ilim0;

3060 }

3061

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

3063

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

3065 /* Yes. */

3066 ds = tens[k];

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

3068 S = mhi = 0;

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

3070 goto no_digits;

3071 goto one_digit;

3072 }

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

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

3075 dval(d) -= L*ds;

3076 #ifdef Check_FLT_ROUNDS

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

3078 if (dval(d) < 0) {

3079 L--;

3080 dval(d) += ds;

3081 }

3082 #endif

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

3084 if (!dval(d)) {

3085 #ifdef SET_INEXACT

3086 inexact = 0;

3087 #endif

3088 break;

3089 }

3090 if (i == ilim) {

3091 #ifdef Honor_FLT_ROUNDS

3092 if (mode > 1)

3093 switch(rounding) {

3094 case 0: goto ret1;

3095 case 2: goto bump_up;

3096 }

3097 #endif

3098 dval(d) += dval(d);

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

3100 bump_up:

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

3102 if (s == s0) {

3103 k++;

3104 *s = '0';

3105 break;

3106 }

3107 ++*s++;

3108 }

3109 break;

3110 }

3111 }

3112 goto ret1;

3113 }

3114

3115 m2 = b2;

3116 m5 = b5;

3117 mhi = mlo = 0;

3118 if (leftright) {

3119 i =

3120 #ifndef Sudden_Underflow

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

3122 #endif

3123 #ifdef IBM

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

3125 #else

3126 1 + P - bbits;

3127 #endif

3128 b2 += i;

3129 s2 += i;

3130 mhi = i2b(1);

3131 }

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

3133 i = m2 < s2 ? m2 : s2;

3134 b2 -= i;

3135 m2 -= i;

3136 s2 -= i;

3137 }

3138 if (b5 > 0) {

3139 if (leftright) {

3140 if (m5 > 0) {

3141 mhi = pow5mult(mhi, m5);

3142 b1 = mult(mhi, b);

3143 Bfree(b);

3144 b = b1;

3145 }

3146 if (j = b5 - m5)

3147 b = pow5mult(b, j);

3148 }

3149 else

3150 b = pow5mult(b, b5);

3151 }

3152 S = i2b(1);

3153 if (s5 > 0)

3154 S = pow5mult(S, s5);

3155

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

3157

3158 spec_case = 0;

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

3160 #ifdef Honor_FLT_ROUNDS

3161 && rounding == 1

3162 #endif

3163 ) {

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

3165 #ifndef Sudden_Underflow

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

3167 #endif

3168 ) {

3169 /* The special case */

3170 b2 += Log2P;

3171 s2 += Log2P;

3172 spec_case = 1;

3173 }

3174 }

3175

3176 /* Arrange for convenient computation of quotients:

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

3178 *

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

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

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

3182 */

3183 #ifdef Pack_32

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

3185 i = 32 - i;

3186 #else

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

3188 i = 16 - i;

3189 #endif

3190 if (i > 4) {

3191 i -= 4;

3192 b2 += i;

3193 m2 += i;

3194 s2 += i;

3195 }

3196 else if (i < 4) {

3197 i += 28;

3198 b2 += i;

3199 m2 += i;

3200 s2 += i;

3201 }

3202 if (b2 > 0)

3203 b = lshift(b, b2);

3204 if (s2 > 0)

3205 S = lshift(S, s2);

3206 if (k_check) {

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

3208 k--;

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

3210 if (leftright)

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

3212 ilim = ilim1;

3213 }

3214 }

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

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

3217 /* no digits, fcvt style */

3218 no_digits:

3219 k = -1 - ndigits;

3220 goto ret;

3221 }

3222 one_digit:

3223 *s++ = '1';

3224 k++;

3225 goto ret;

3226 }

3227 if (leftright) {

3228 if (m2 > 0)

3229 mhi = lshift(mhi, m2);

3230

3231 /* Compute mlo -- check for special case

3232 * that d is a normalized power of 2.

3233 */

3234

3235 mlo = mhi;

3236 if (spec_case) {

3237 mhi = Balloc(mhi->k);

3238 Bcopy(mhi, mlo);

3239 mhi = lshift(mhi, Log2P);

3240 }

3241

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

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

3244 /* Do we yet have the shortest decimal string

3245 * that will round to d?

3246 */

3247 j = cmp(b, mlo);

3248 delta = diff(S, mhi);

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

3250 Bfree(delta);

3251 #ifndef ROUND_BIASED

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

3253 #ifdef Honor_FLT_ROUNDS

3254 && rounding >= 1

3255 #endif

3256 ) {

3257 if (dig == '9')

3258 goto round_9_up;

3259 if (j > 0)

3260 dig++;

3261 #ifdef SET_INEXACT

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

3263 inexact = 0;

3264 #endif

3265 *s++ = dig;

3266 goto ret;

3267 }

3268 #endif

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

3270 #ifndef ROUND_BIASED

3271 && !(word1(d) & 1)

3272 #endif

3273 ) {

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

3275 #ifdef SET_INEXACT

3276 inexact = 0;

3277 #endif

3278 goto accept_dig;

3279 }

3280 #ifdef Honor_FLT_ROUNDS

3281 if (mode > 1)

3282 switch(rounding) {

3283 case 0: goto accept_dig;

3284 case 2: goto keep_dig;

3285 }

3286 #endif /Honor_FLT_ROUNDS/

3287 if (j1 > 0) {

3288 b = lshift(b, 1);

3289 j1 = cmp(b, S);

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

3291 && dig++ == '9')

3292 goto round_9_up;

3293 }

3294 accept_dig:

3295 *s++ = dig;

3296 goto ret;

3297 }

3298 if (j1 > 0) {

3299 #ifdef Honor_FLT_ROUNDS

3300 if (!rounding)

3301 goto accept_dig;

3302 #endif

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

3304 round_9_up:

3305 *s++ = '9';

3306 goto roundoff;

3307 }

3308 *s++ = dig + 1;

3309 goto ret;

3310 }

3311 #ifdef Honor_FLT_ROUNDS

3312 keep_dig:

3313 #endif

3314 *s++ = dig;

3315 if (i == ilim)

3316 break;

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

3318 if (mlo == mhi)

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

3320 else {

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

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

3323 }

3324 }

3325 }

3326 else

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

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

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

3330 #ifdef SET_INEXACT

3331 inexact = 0;

3332 #endif

3333 goto ret;

3334 }

3335 if (i >= ilim)

3336 break;

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

3338 }

3339

3340 /* Round off last digit */

3341

3342 #ifdef Honor_FLT_ROUNDS

3343 switch(rounding) {

3344 case 0: goto trimzeros;

3345 case 2: goto roundoff;

3346 }

3347 #endif

3348 b = lshift(b, 1);

3349 j = cmp(b, S);

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

3351 roundoff:

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

3353 if (s == s0) {

3354 k++;

3355 *s++ = '1';

3356 goto ret;

3357 }

3358 ++*s++;

3359 }

3360 else {

3361 #ifdef Honor_FLT_ROUNDS

3362 trimzeros:

3363 #endif

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

3365 s++;

3366 }

3367 ret:

3368 Bfree(S);

3369 if (mhi) {

3370 if (mlo && mlo != mhi)

3371 Bfree(mlo);

3372 Bfree(mhi);

3373 }

3374 ret1:

3375 #ifdef SET_INEXACT

3376 if (inexact) {

3377 if (!oldinexact) {

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

3379 word1(d) = 0;

3380 dval(d) += 1.;

3381 }

3382 }

3383 else if (!oldinexact)

3384 clear_inexact();

3385 #endif

3386 Bfree(b);

3387 *s = 0;

3388 *decpt = k + 1;

3389 if (rve)

3390 *rve = s;

3391 return s0;

3392 }

3393 #ifdef __cplusplus

3394 }

3395 #endif

3396

3397 PR_IMPLEMENT(PRStatus)

3398 PR_dtoa(PRFloat64 d, PRIntn mode, PRIntn ndigits,

3399 PRIntn decpt, PRIntn sign, char *rve, char buf, PRSize bufsize)

3400 {

3401 char *result;

3402 PRSize resultlen;

3403 PRStatus rv = PR_FAILURE;

3404

3405 if (!_pr_initialized) _PR_ImplicitInitialization();

3406

3407 if (mode < 0 \|\| mode > 3) {

3408 PR_SetError(PR_INVALID_ARGUMENT_ERROR, 0);

3409 return rv;

3410 }

3411 result = dtoa(d, mode, ndigits, decpt, sign, rve);

3412 if (!result) {

3413 PR_SetError(PR_OUT_OF_MEMORY_ERROR, 0);

3414 return rv;

3415 }

3416 resultlen = strlen(result)+1;

3417 if (bufsize < resultlen) {

3418 PR_SetError(PR_BUFFER_OVERFLOW_ERROR, 0);

3419 } else {

3420 memcpy(buf, result, resultlen);

3421 if (rve) {

3422 rve = buf + (rve - result);

3423 }

3424 rv = PR_SUCCESS;

3425 }

3426 freedtoa(result);

3427 return rv;

3428 }

3429

3430 /*

3431 ** conversion routines for floating point

3432 ** prcsn - number of digits of precision to generate floating

3433 ** point value.

3434 ** This should be reparameterized so that you can send in a

3435 ** prcn for the positive and negative ranges. For now,

3436 ** conform to the ECMA JavaScript spec which says numbers

3437 ** less than 1e-6 are in scientific notation.

3438 ** Also, the ECMA spec says that there should always be a

3439 ** '+' or '-' after the 'e' in scientific notation

3440 */

3441 PR_IMPLEMENT(void)

3442 PR_cnvtf(char *buf, int bufsz, int prcsn, double dfval)

3443 {

3444 PRIntn decpt, sign, numdigits;

3445 char num, nump;

3446 char *bufp = buf;

3447 char *endnum;

3448 U fval;

3449

3450 dval(fval) = dfval;

3451 /* If anything fails, we store an empty string in 'buf' */

3452 num = (char*)PR_MALLOC(bufsz);

3453 if (num == NULL) {

3454 buf[0] = '\0';

3455 return;

3456 }

3457 /* XXX Why use mode 1? */

3458 if (PR_dtoa(dval(fval),1,prcsn,&decpt,&sign,&endnum,num,bufsz)

3459 == PR_FAILURE) {

3460 buf[0] = '\0';

3461 goto done;

3462 }

3463 numdigits = endnum - num;

3464 nump = num;

3465

3466 if (sign &&

3467 !(word0(fval) == Sign_bit && word1(fval) == 0) &&

3468 !((word0(fval) & Exp_mask) == Exp_mask &&

3469 (word1(fval) \|\| (word0(fval) & 0xfffff)))) {

3470 *bufp++ = '-';

3471 }

3472

3473 if (decpt == 9999) {

3474 while ((bufp++ = nump++) != 0) {} /* nothing to execute */

3475 goto done;

3476 }

3477

3478 if (decpt > (prcsn+1) \|\| decpt < -(prcsn-1) \|\| decpt < -5) {

3479 bufp++ = nump++;

3480 if (numdigits != 1) {

3481 *bufp++ = '.';

3482 }

3483

3484 while (*nump != '\0') {

3485 bufp++ = nump++;

3486 }

3487 *bufp++ = 'e';

3488 PR_snprintf(bufp, bufsz - (bufp - buf), "%+d", decpt-1);

3489 } else if (decpt >= 0) {

3490 if (decpt == 0) {

3491 *bufp++ = '0';

3492 } else {

3493 while (decpt--) {

3494 if (*nump != '\0') {

3495 bufp++ = nump++;

3496 } else {

3497 *bufp++ = '0';

3498 }

3499 }

3500 }

3501 if (*nump != '\0') {

3502 *bufp++ = '.';

3503 while (*nump != '\0') {

3504 bufp++ = nump++;

3505 }

3506 }

3507 *bufp++ = '\0';

3508 } else if (decpt < 0) {

3509 *bufp++ = '0';

3510 *bufp++ = '.';

3511 while (decpt++) {

3512 *bufp++ = '0';

3513 }

3514

3515 while (*nump != '\0') {

3516 bufp++ = nump++;

3517 }

3518 *bufp++ = '\0';

3519 }

3520 done:

3521 PR_DELETE(num);

3522 }

OLD	NEW

« no previous file with comments | « nspr/pr/src/misc/prcountr.c ('k') | nspr/pr/src/misc/prenv.c » ('j') | no next file with comments »