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 (d*10^k) * 10^(e-k) with just one roundoff. 54 * 3. Rather than a bit-at-a-time adjustment of the binary 55 * result in the hard case, we use floating-point 56 * arithmetic to determine the adjustment to within 57 * one bit; only in really hard cases do we need to 58 * compute a second residual. 59 * 4. Because of 3., we don't need a large table of powers of 10 60 * for ten-to-e (just some small tables, e.g. of 10^k 61 * for 0 <= k <= 22). 62 */ 63 64 /* 65 * #define IEEE_8087 for IEEE-arithmetic machines where the least 66 * significant byte has the lowest address. 67 * #define IEEE_MC68k for IEEE-arithmetic machines where the most 68 * significant byte has the lowest address. 69 * #define Long int on machines with 32-bit ints and 64-bit longs. 70 * #define IBM for IBM mainframe-style floating-point arithmetic. 71 * #define VAX for VAX-style floating-point arithmetic (D_floating). 72 * #define No_leftright to omit left-right logic in fast floating-point 73 * computation of dtoa. 74 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 75 * and strtod and dtoa should round accordingly. 76 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 77 * and Honor_FLT_ROUNDS is not #defined. 78 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines 79 * that use extended-precision instructions to compute rounded 80 * products and quotients) with IBM. 81 * #define ROUND_BIASED for IEEE-format with biased rounding. 82 * #define Inaccurate_Divide for IEEE-format with correctly rounded 83 * products but inaccurate quotients, e.g., for Intel i860. 84 * #define NO_LONG_LONG on machines that do not have a "long long" 85 * integer type (of >= 64 bits). On such machines, you can 86 * #define Just_16 to store 16 bits per 32-bit Long when doing 87 * high-precision integer arithmetic. Whether this speeds things 88 * up or slows things down depends on the machine and the number 89 * being converted. If long long is available and the name is 90 * something other than "long long", #define Llong to be the name, 91 * and if "unsigned Llong" does not work as an unsigned version of 92 * Llong, #define #ULLong to be the corresponding unsigned type. 93 * #define KR_headers for old-style C function headers. 94 * #define Bad_float_h if your system lacks a float.h or if it does not 95 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP, 96 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX. 97 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n) 98 * if memory is available and otherwise does something you deem 99 * appropriate. If MALLOC is undefined, malloc will be invoked 100 * directly -- and assumed always to succeed. 101 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making 102 * memory allocations from a private pool of memory when possible. 103 * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes, 104 * unless #defined to be a different length. This default length 105 * suffices to get rid of MALLOC calls except for unusual cases, 106 * such as decimal-to-binary conversion of a very long string of 107 * digits. The longest string dtoa can return is about 751 bytes 108 * long. For conversions by strtod of strings of 800 digits and 109 * all dtoa conversions in single-threaded executions with 8-byte 110 * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte 111 * pointers, PRIVATE_MEM >= 7112 appears adequate. 112 * #define INFNAN_CHECK on IEEE systems to cause strtod to check for 113 * Infinity and NaN (case insensitively). On some systems (e.g., 114 * some HP systems), it may be necessary to #define NAN_WORD0 115 * appropriately -- to the most significant word of a quiet NaN. 116 * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.) 117 * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined, 118 * strtod also accepts (case insensitively) strings of the form 119 * NaN(x), where x is a string of hexadecimal digits and spaces; 120 * if there is only one string of hexadecimal digits, it is taken 121 * for the 52 fraction bits of the resulting NaN; if there are two 122 * or more strings of hex digits, the first is for the high 20 bits, 123 * the second and subsequent for the low 32 bits, with intervening 124 * white space ignored; but if this results in none of the 52 125 * fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0 126 * and NAN_WORD1 are used instead. 127 * #define MULTIPLE_THREADS if the system offers preemptively scheduled 128 * multiple threads. In this case, you must provide (or suitably 129 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed 130 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed 131 * in pow5mult, ensures lazy evaluation of only one copy of high 132 * powers of 5; omitting this lock would introduce a small 133 * probability of wasting memory, but would otherwise be harmless.) 134 * You must also invoke freedtoa(s) to free the value s returned by 135 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined. 136 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that 137 * avoids underflows on inputs whose result does not underflow. 138 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format 139 * floating-point numbers and flushes underflows to zero rather 140 * than implementing gradual underflow, then you must also #define 141 * Sudden_Underflow. 142 * #define YES_ALIAS to permit aliasing certain double values with 143 * arrays of ULongs. This leads to slightly better code with 144 * some compilers and was always used prior to 19990916, but it 145 * is not strictly legal and can cause trouble with aggressively 146 * optimizing compilers (e.g., gcc 2.95.1 under -O2). 147 * #define USE_LOCALE to use the current locale's decimal_point value. 148 * #define SET_INEXACT if IEEE arithmetic is being used and extra 149 * computation should be done to set the inexact flag when the 150 * result is inexact and avoid setting inexact when the result 151 * is exact. In this case, dtoa.c must be compiled in 152 * an environment, perhaps provided by #include "dtoa.c" in a 153 * suitable wrapper, that defines two functions, 154 * int get_inexact(void); 155 * void clear_inexact(void); 156 * such that get_inexact() returns a nonzero value if the 157 * inexact bit is already set, and clear_inexact() sets the 158 * inexact bit to 0. When SET_INEXACT is #defined, strtod 159 * also does extra computations to set the underflow and overflow 160 * flags when appropriate (i.e., when the result is tiny and 161 * inexact or when it is a numeric value rounded to +-infinity). 162 * #define NO_ERRNO if strtod should not assign errno = ERANGE when 163 * the result overflows to +-Infinity or underflows to 0. 164 */ 165 166 #ifndef Long 167 #define Long long 168 #endif 169 #ifndef ULong 170 typedef unsigned Long ULong; 171 #endif 172 173 #ifdef DEBUG 174 #include "stdio.h" 175 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} 176 #endif 177 178 #include "stdlib.h" 179 #include "string.h" 180 181 #ifdef USE_LOCALE 182 #include "locale.h" 183 #endif 184 185 #ifdef MALLOC 186 #ifdef KR_headers 187 extern char *MALLOC(); 188 #else 189 extern void *MALLOC(size_t); 190 #endif 191 #else 192 #define MALLOC malloc 193 #endif 194 195 #ifndef Omit_Private_Memory 196 #ifndef PRIVATE_MEM 197 #define PRIVATE_MEM 2304 198 #endif 199 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double)) 200 static double private_mem[PRIVATE_mem], *pmem_next = private_mem; 201 #endif 202 203 #undef IEEE_Arith 204 #undef Avoid_Underflow 205 #ifdef IEEE_MC68k 206 #define IEEE_Arith 207 #endif 208 #ifdef IEEE_8087 209 #define IEEE_Arith 210 #endif 211 212 #include "errno.h" 213 214 #ifdef Bad_float_h 215 216 #ifdef IEEE_Arith 217 #define DBL_DIG 15 218 #define DBL_MAX_10_EXP 308 219 #define DBL_MAX_EXP 1024 220 #define FLT_RADIX 2 221 #endif /*IEEE_Arith*/ 222 223 #ifdef IBM 224 #define DBL_DIG 16 225 #define DBL_MAX_10_EXP 75 226 #define DBL_MAX_EXP 63 227 #define FLT_RADIX 16 228 #define DBL_MAX 7.2370055773322621e+75 229 #endif 230 231 #ifdef VAX 232 #define DBL_DIG 16 233 #define DBL_MAX_10_EXP 38 234 #define DBL_MAX_EXP 127 235 #define FLT_RADIX 2 236 #define DBL_MAX 1.7014118346046923e+38 237 #endif 238 239 #ifndef LONG_MAX 240 #define LONG_MAX 2147483647 241 #endif 242 243 #else /* ifndef Bad_float_h */ 244 #include "float.h" 245 #endif /* Bad_float_h */ 246 247 #ifndef __MATH_H__ 248 #include "math.h" 249 #endif 250 251 #ifdef __cplusplus 252 extern "C" { 253 #endif 254 255 #ifndef CONST 256 #ifdef KR_headers 257 #define CONST /* blank */ 258 #else 259 #define CONST const 260 #endif 261 #endif 262 263 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1 264 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined. 265 #endif 266 267 typedef union { double d; ULong L[2]; } U; 268 269 #ifdef YES_ALIAS 270 #define dval(x) x 271 #ifdef IEEE_8087 272 #define word0(x) ((ULong *)&x)[1] 273 #define word1(x) ((ULong *)&x)[0] 274 #else 275 #define word0(x) ((ULong *)&x)[0] 276 #define word1(x) ((ULong *)&x)[1] 277 #endif 278 #else 279 #ifdef IEEE_8087 280 #define word0(x) ((U*)&x)->L[1] 281 #define word1(x) ((U*)&x)->L[0] 282 #else 283 #define word0(x) ((U*)&x)->L[0] 284 #define word1(x) ((U*)&x)->L[1] 285 #endif 286 #define dval(x) ((U*)&x)->d 287 #endif 288 289 /* The following definition of Storeinc is appropriate for MIPS processors. 290 * An alternative that might be better on some machines is 291 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) 292 */ 293 #if defined(IEEE_8087) + defined(VAX) 294 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \ 295 ((unsigned short *)a)[0] = (unsigned short)c, a++) 296 #else 297 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \ 298 ((unsigned short *)a)[1] = (unsigned short)c, a++) 299 #endif 300 301 /* #define P DBL_MANT_DIG */ 302 /* Ten_pmax = floor(P*log(2)/log(5)) */ 303 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ 304 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ 305 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ 306 307 #ifdef IEEE_Arith 308 #define Exp_shift 20 309 #define Exp_shift1 20 310 #define Exp_msk1 0x100000 311 #define Exp_msk11 0x100000 312 #define Exp_mask 0x7ff00000 313 #define P 53 314 #define Bias 1023 315 #define Emin (-1022) 316 #define Exp_1 0x3ff00000 317 #define Exp_11 0x3ff00000 318 #define Ebits 11 319 #define Frac_mask 0xfffff 320 #define Frac_mask1 0xfffff 321 #define Ten_pmax 22 322 #define Bletch 0x10 323 #define Bndry_mask 0xfffff 324 #define Bndry_mask1 0xfffff 325 #define LSB 1 326 #define Sign_bit 0x80000000 327 #define Log2P 1 328 #define Tiny0 0 329 #define Tiny1 1 330 #define Quick_max 14 331 #define Int_max 14 332 #ifndef NO_IEEE_Scale 333 #define Avoid_Underflow 334 #ifdef Flush_Denorm /* debugging option */ 335 #undef Sudden_Underflow 336 #endif 337 #endif 338 339 #ifndef Flt_Rounds 340 #ifdef FLT_ROUNDS 341 #define Flt_Rounds FLT_ROUNDS 342 #else 343 #define Flt_Rounds 1 344 #endif 345 #endif /*Flt_Rounds*/ 346 347 #ifdef Honor_FLT_ROUNDS 348 #define Rounding rounding 349 #undef Check_FLT_ROUNDS 350 #define Check_FLT_ROUNDS 351 #else 352 #define Rounding Flt_Rounds 353 #endif 354 355 #else /* ifndef IEEE_Arith */ 356 #undef Check_FLT_ROUNDS 357 #undef Honor_FLT_ROUNDS 358 #undef SET_INEXACT 359 #undef Sudden_Underflow 360 #define Sudden_Underflow 361 #ifdef IBM 362 #undef Flt_Rounds 363 #define Flt_Rounds 0 364 #define Exp_shift 24 365 #define Exp_shift1 24 366 #define Exp_msk1 0x1000000 367 #define Exp_msk11 0x1000000 368 #define Exp_mask 0x7f000000 369 #define P 14 370 #define Bias 65 371 #define Exp_1 0x41000000 372 #define Exp_11 0x41000000 373 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */ 374 #define Frac_mask 0xffffff 375 #define Frac_mask1 0xffffff 376 #define Bletch 4 377 #define Ten_pmax 22 378 #define Bndry_mask 0xefffff 379 #define Bndry_mask1 0xffffff 380 #define LSB 1 381 #define Sign_bit 0x80000000 382 #define Log2P 4 383 #define Tiny0 0x100000 384 #define Tiny1 0 385 #define Quick_max 14 386 #define Int_max 15 387 #else /* VAX */ 388 #undef Flt_Rounds 389 #define Flt_Rounds 1 390 #define Exp_shift 23 391 #define Exp_shift1 7 392 #define Exp_msk1 0x80 393 #define Exp_msk11 0x800000 394 #define Exp_mask 0x7f80 395 #define P 56 396 #define Bias 129 397 #define Exp_1 0x40800000 398 #define Exp_11 0x4080 399 #define Ebits 8 400 #define Frac_mask 0x7fffff 401 #define Frac_mask1 0xffff007f 402 #define Ten_pmax 24 403 #define Bletch 2 404 #define Bndry_mask 0xffff007f 405 #define Bndry_mask1 0xffff007f 406 #define LSB 0x10000 407 #define Sign_bit 0x8000 408 #define Log2P 1 409 #define Tiny0 0x80 410 #define Tiny1 0 411 #define Quick_max 15 412 #define Int_max 15 413 #endif /* IBM, VAX */ 414 #endif /* IEEE_Arith */ 415 416 #ifndef IEEE_Arith 417 #define ROUND_BIASED 418 #endif 419 420 #ifdef RND_PRODQUOT 421 #define rounded_product(a,b) a = rnd_prod(a, b) 422 #define rounded_quotient(a,b) a = rnd_quot(a, b) 423 #ifdef KR_headers 424 extern double rnd_prod(), rnd_quot(); 425 #else 426 extern double rnd_prod(double, double), rnd_quot(double, double); 427 #endif 428 #else 429 #define rounded_product(a,b) a *= b 430 #define rounded_quotient(a,b) a /= b 431 #endif 432 433 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) 434 #define Big1 0xffffffff 435 436 #ifndef Pack_32 437 #define Pack_32 438 #endif 439 440 #ifdef KR_headers 441 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff) 442 #else 443 #define FFFFFFFF 0xffffffffUL 444 #endif 445 446 #ifdef NO_LONG_LONG 447 #undef ULLong 448 #ifdef Just_16 449 #undef Pack_32 450 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long. 451 * This makes some inner loops simpler and sometimes saves work 452 * during multiplications, but it often seems to make things slightly 453 * slower. Hence the default is now to store 32 bits per Long. 454 */ 455 #endif 456 #else /* long long available */ 457 #ifndef Llong 458 #define Llong long long 459 #endif 460 #ifndef ULLong 461 #define ULLong unsigned Llong 462 #endif 463 #endif /* NO_LONG_LONG */ 464 465 #ifndef MULTIPLE_THREADS 466 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/ 467 #define FREE_DTOA_LOCK(n) /*nothing*/ 468 #endif 469 470 #define Kmax 15 471 472 #ifdef __cplusplus 473 extern "C" double strtod(const char *s00, char **se); 474 extern "C" char *dtoa(double d, int mode, int ndigits, 475 int *decpt, int *sign, char **rve); 476 #endif 477 478 struct 479 Bigint { 480 struct Bigint *next; 481 int k, maxwds, sign, wds; 482 ULong x[1]; 483 }; 484 485 typedef struct Bigint Bigint; 486 487 static Bigint *freelist[Kmax+1]; 488 489 static Bigint * 490 Balloc 491 #ifdef KR_headers 492 (k) int k; 493 #else 494 (int k) 495 #endif 496 { 497 int x; 498 Bigint *rv; 499 #ifndef Omit_Private_Memory 500 unsigned int len; 501 #endif 502 503 ACQUIRE_DTOA_LOCK(0); 504 /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0). */ 505 /* but this case seems very unlikely. */ 506 if (k <= Kmax && (rv = freelist[k])) { 507 freelist[k] = rv->next; 508 } 509 else { 510 x = 1 << k; 511 #ifdef Omit_Private_Memory 512 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong)); 513 #else 514 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1) 515 /sizeof(double); 516 if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) { 517 rv = (Bigint*)pmem_next; 518 pmem_next += len; 519 } 520 else 521 rv = (Bigint*)MALLOC(len*sizeof(double)); 522 #endif 523 rv->k = k; 524 rv->maxwds = x; 525 } 526 FREE_DTOA_LOCK(0); 527 rv->sign = rv->wds = 0; 528 return rv; 529 } 530 531 static void 532 Bfree 533 #ifdef KR_headers 534 (v) Bigint *v; 535 #else 536 (Bigint *v) 537 #endif 538 { 539 if (v) { 540 if (v->k > Kmax) 541 free((void*)v); 542 else { 543 ACQUIRE_DTOA_LOCK(0); 544 v->next = freelist[v->k]; 545 freelist[v->k] = v; 546 FREE_DTOA_LOCK(0); 547 } 548 } 549 } 550 551 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \ 552 y->wds*sizeof(Long) + 2*sizeof(int)) 553 554 static Bigint * 555 multadd 556 #ifdef KR_headers 557 (b, m, a) Bigint *b; int m, a; 558 #else 559 (Bigint *b, int m, int a) /* multiply by m and add a */ 560 #endif 561 { 562 int i, wds; 563 #ifdef ULLong 564 ULong *x; 565 ULLong carry, y; 566 #else 567 ULong carry, *x, y; 568 #ifdef Pack_32 569 ULong xi, z; 570 #endif 571 #endif 572 Bigint *b1; 573 574 wds = b->wds; 575 x = b->x; 576 i = 0; 577 carry = a; 578 do { 579 #ifdef ULLong 580 y = *x * (ULLong)m + carry; 581 carry = y >> 32; 582 *x++ = y & FFFFFFFF; 583 #else 584 #ifdef Pack_32 585 xi = *x; 586 y = (xi & 0xffff) * m + carry; 587 z = (xi >> 16) * m + (y >> 16); 588 carry = z >> 16; 589 *x++ = (z << 16) + (y & 0xffff); 590 #else 591 y = *x * m + carry; 592 carry = y >> 16; 593 *x++ = y & 0xffff; 594 #endif 595 #endif 596 } 597 while(++i < wds); 598 if (carry) { 599 if (wds >= b->maxwds) { 600 b1 = Balloc(b->k+1); 601 Bcopy(b1, b); 602 Bfree(b); 603 b = b1; 604 } 605 b->x[wds++] = carry; 606 b->wds = wds; 607 } 608 return b; 609 } 610 611 static Bigint * 612 s2b 613 #ifdef KR_headers 614 (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9; 615 #else 616 (CONST char *s, int nd0, int nd, ULong y9) 617 #endif 618 { 619 Bigint *b; 620 int i, k; 621 Long x, y; 622 623 x = (nd + 8) / 9; 624 for(k = 0, y = 1; x > y; y <<= 1, k++) ; 625 #ifdef Pack_32 626 b = Balloc(k); 627 b->x[0] = y9; 628 b->wds = 1; 629 #else 630 b = Balloc(k+1); 631 b->x[0] = y9 & 0xffff; 632 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1; 633 #endif 634 635 i = 9; 636 if (9 < nd0) { 637 s += 9; 638 do b = multadd(b, 10, *s++ - '0'); 639 while(++i < nd0); 640 s++; 641 } 642 else 643 s += 10; 644 for(; i < nd; i++) 645 b = multadd(b, 10, *s++ - '0'); 646 return b; 647 } 648 649 static int 650 hi0bits 651 #ifdef KR_headers 652 (x) register ULong x; 653 #else 654 (register ULong x) 655 #endif 656 { 657 register int k = 0; 658 659 if (!(x & 0xffff0000)) { 660 k = 16; 661 x <<= 16; 662 } 663 if (!(x & 0xff000000)) { 664 k += 8; 665 x <<= 8; 666 } 667 if (!(x & 0xf0000000)) { 668 k += 4; 669 x <<= 4; 670 } 671 if (!(x & 0xc0000000)) { 672 k += 2; 673 x <<= 2; 674 } 675 if (!(x & 0x80000000)) { 676 k++; 677 if (!(x & 0x40000000)) 678 return 32; 679 } 680 return k; 681 } 682 683 static int 684 lo0bits 685 #ifdef KR_headers 686 (y) ULong *y; 687 #else 688 (ULong *y) 689 #endif 690 { 691 register int k; 692 register ULong x = *y; 693 694 if (x & 7) { 695 if (x & 1) 696 return 0; 697 if (x & 2) { 698 *y = x >> 1; 699 return 1; 700 } 701 *y = x >> 2; 702 return 2; 703 } 704 k = 0; 705 if (!(x & 0xffff)) { 706 k = 16; 707 x >>= 16; 708 } 709 if (!(x & 0xff)) { 710 k += 8; 711 x >>= 8; 712 } 713 if (!(x & 0xf)) { 714 k += 4; 715 x >>= 4; 716 } 717 if (!(x & 0x3)) { 718 k += 2; 719 x >>= 2; 720 } 721 if (!(x & 1)) { 722 k++; 723 x >>= 1; 724 if (!x) 725 return 32; 726 } 727 *y = x; 728 return k; 729 } 730 731 static Bigint * 732 i2b 733 #ifdef KR_headers 734 (i) int i; 735 #else 736 (int i) 737 #endif 738 { 739 Bigint *b; 740 741 b = Balloc(1); 742 b->x[0] = i; 743 b->wds = 1; 744 return b; 745 } 746 747 static Bigint * 748 mult 749 #ifdef KR_headers 750 (a, b) Bigint *a, *b; 751 #else 752 (Bigint *a, Bigint *b) 753 #endif 754 { 755 Bigint *c; 756 int k, wa, wb, wc; 757 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; 758 ULong y; 759 #ifdef ULLong 760 ULLong carry, z; 761 #else 762 ULong carry, z; 763 #ifdef Pack_32 764 ULong z2; 765 #endif 766 #endif 767 768 if (a->wds < b->wds) { 769 c = a; 770 a = b; 771 b = c; 772 } 773 k = a->k; 774 wa = a->wds; 775 wb = b->wds; 776 wc = wa + wb; 777 if (wc > a->maxwds) 778 k++; 779 c = Balloc(k); 780 for(x = c->x, xa = x + wc; x < xa; x++) 781 *x = 0; 782 xa = a->x; 783 xae = xa + wa; 784 xb = b->x; 785 xbe = xb + wb; 786 xc0 = c->x; 787 #ifdef ULLong 788 for(; xb < xbe; xc0++) { 789 if ((y = *xb++)) { 790 x = xa; 791 xc = xc0; 792 carry = 0; 793 do { 794 z = *x++ * (ULLong)y + *xc + carry; 795 carry = z >> 32; 796 *xc++ = z & FFFFFFFF; 797 } 798 while(x < xae); 799 *xc = carry; 800 } 801 } 802 #else 803 #ifdef Pack_32 804 for(; xb < xbe; xb++, xc0++) { 805 if (y = *xb & 0xffff) { 806 x = xa; 807 xc = xc0; 808 carry = 0; 809 do { 810 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 811 carry = z >> 16; 812 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 813 carry = z2 >> 16; 814 Storeinc(xc, z2, z); 815 } 816 while(x < xae); 817 *xc = carry; 818 } 819 if (y = *xb >> 16) { 820 x = xa; 821 xc = xc0; 822 carry = 0; 823 z2 = *xc; 824 do { 825 z = (*x & 0xffff) * y + (*xc >> 16) + carry; 826 carry = z >> 16; 827 Storeinc(xc, z, z2); 828 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 829 carry = z2 >> 16; 830 } 831 while(x < xae); 832 *xc = z2; 833 } 834 } 835 #else 836 for(; xb < xbe; xc0++) { 837 if (y = *xb++) { 838 x = xa; 839 xc = xc0; 840 carry = 0; 841 do { 842 z = *x++ * y + *xc + carry; 843 carry = z >> 16; 844 *xc++ = z & 0xffff; 845 } 846 while(x < xae); 847 *xc = carry; 848 } 849 } 850 #endif 851 #endif 852 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ; 853 c->wds = wc; 854 return c; 855 } 856 857 static Bigint *p5s; 858 859 static Bigint * 860 pow5mult 861 #ifdef KR_headers 862 (b, k) Bigint *b; int k; 863 #else 864 (Bigint *b, int k) 865 #endif 866 { 867 Bigint *b1, *p5, *p51; 868 int i; 869 static int p05[3] = { 5, 25, 125 }; 870 871 if ((i = k & 3)) 872 b = multadd(b, p05[i-1], 0); 873 874 if (!(k >>= 2)) 875 return b; 876 if (!(p5 = p5s)) { 877 /* first time */ 878 #ifdef MULTIPLE_THREADS 879 ACQUIRE_DTOA_LOCK(1); 880 if (!(p5 = p5s)) { 881 p5 = p5s = i2b(625); 882 p5->next = 0; 883 } 884 FREE_DTOA_LOCK(1); 885 #else 886 p5 = p5s = i2b(625); 887 p5->next = 0; 888 #endif 889 } 890 for(;;) { 891 if (k & 1) { 892 b1 = mult(b, p5); 893 Bfree(b); 894 b = b1; 895 } 896 if (!(k >>= 1)) 897 break; 898 if (!(p51 = p5->next)) { 899 #ifdef MULTIPLE_THREADS 900 ACQUIRE_DTOA_LOCK(1); 901 if (!(p51 = p5->next)) { 902 p51 = p5->next = mult(p5,p5); 903 p51->next = 0; 904 } 905 FREE_DTOA_LOCK(1); 906 #else 907 p51 = p5->next = mult(p5,p5); 908 p51->next = 0; 909 #endif 910 } 911 p5 = p51; 912 } 913 return b; 914 } 915 916 static Bigint * 917 lshift 918 #ifdef KR_headers 919 (b, k) Bigint *b; int k; 920 #else 921 (Bigint *b, int k) 922 #endif 923 { 924 int i, k1, n, n1; 925 Bigint *b1; 926 ULong *x, *x1, *xe, z; 927 928 #ifdef Pack_32 929 n = k >> 5; 930 #else 931 n = k >> 4; 932 #endif 933 k1 = b->k; 934 n1 = n + b->wds + 1; 935 for(i = b->maxwds; n1 > i; i <<= 1) 936 k1++; 937 b1 = Balloc(k1); 938 x1 = b1->x; 939 for(i = 0; i < n; i++) 940 *x1++ = 0; 941 x = b->x; 942 xe = x + b->wds; 943 #ifdef Pack_32 944 if (k &= 0x1f) { 945 k1 = 32 - k; 946 z = 0; 947 do { 948 *x1++ = *x << k | z; 949 z = *x++ >> k1; 950 } 951 while(x < xe); 952 if ((*x1 = z)) 953 ++n1; 954 } 955 #else 956 if (k &= 0xf) { 957 k1 = 16 - k; 958 z = 0; 959 do { 960 *x1++ = *x << k & 0xffff | z; 961 z = *x++ >> k1; 962 } 963 while(x < xe); 964 if (*x1 = z) 965 ++n1; 966 } 967 #endif 968 else do 969 *x1++ = *x++; 970 while(x < xe); 971 b1->wds = n1 - 1; 972 Bfree(b); 973 return b1; 974 } 975 976 static int 977 cmp 978 #ifdef KR_headers 979 (a, b) Bigint *a, *b; 980 #else 981 (Bigint *a, Bigint *b) 982 #endif 983 { 984 ULong *xa, *xa0, *xb, *xb0; 985 int i, j; 986 987 i = a->wds; 988 j = b->wds; 989 #ifdef DEBUG 990 if (i > 1 && !a->x[i-1]) 991 Bug("cmp called with a->x[a->wds-1] == 0"); 992 if (j > 1 && !b->x[j-1]) 993 Bug("cmp called with b->x[b->wds-1] == 0"); 994 #endif 995 if (i -= j) 996 return i; 997 xa0 = a->x; 998 xa = xa0 + j; 999 xb0 = b->x; 1000 xb = xb0 + j; 1001 for(;;) { 1002 if (*--xa != *--xb) 1003 return *xa < *xb ? -1 : 1; 1004 if (xa <= xa0) 1005 break; 1006 } 1007 return 0; 1008 } 1009 1010 static Bigint * 1011 diff 1012 #ifdef KR_headers 1013 (a, b) Bigint *a, *b; 1014 #else 1015 (Bigint *a, Bigint *b) 1016 #endif 1017 { 1018 Bigint *c; 1019 int i, wa, wb; 1020 ULong *xa, *xae, *xb, *xbe, *xc; 1021 #ifdef ULLong 1022 ULLong borrow, y; 1023 #else 1024 ULong borrow, y; 1025 #ifdef Pack_32 1026 ULong z; 1027 #endif 1028 #endif 1029 1030 i = cmp(a,b); 1031 if (!i) { 1032 c = Balloc(0); 1033 c->wds = 1; 1034 c->x[0] = 0; 1035 return c; 1036 } 1037 if (i < 0) { 1038 c = a; 1039 a = b; 1040 b = c; 1041 i = 1; 1042 } 1043 else 1044 i = 0; 1045 c = Balloc(a->k); 1046 c->sign = i; 1047 wa = a->wds; 1048 xa = a->x; 1049 xae = xa + wa; 1050 wb = b->wds; 1051 xb = b->x; 1052 xbe = xb + wb; 1053 xc = c->x; 1054 borrow = 0; 1055 #ifdef ULLong 1056 do { 1057 y = (ULLong)*xa++ - *xb++ - borrow; 1058 borrow = y >> 32 & (ULong)1; 1059 *xc++ = y & FFFFFFFF; 1060 } 1061 while(xb < xbe); 1062 while(xa < xae) { 1063 y = *xa++ - borrow; 1064 borrow = y >> 32 & (ULong)1; 1065 *xc++ = y & FFFFFFFF; 1066 } 1067 #else 1068 #ifdef Pack_32 1069 do { 1070 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow; 1071 borrow = (y & 0x10000) >> 16; 1072 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow; 1073 borrow = (z & 0x10000) >> 16; 1074 Storeinc(xc, z, y); 1075 } 1076 while(xb < xbe); 1077 while(xa < xae) { 1078 y = (*xa & 0xffff) - borrow; 1079 borrow = (y & 0x10000) >> 16; 1080 z = (*xa++ >> 16) - borrow; 1081 borrow = (z & 0x10000) >> 16; 1082 Storeinc(xc, z, y); 1083 } 1084 #else 1085 do { 1086 y = *xa++ - *xb++ - borrow; 1087 borrow = (y & 0x10000) >> 16; 1088 *xc++ = y & 0xffff; 1089 } 1090 while(xb < xbe); 1091 while(xa < xae) { 1092 y = *xa++ - borrow; 1093 borrow = (y & 0x10000) >> 16; 1094 *xc++ = y & 0xffff; 1095 } 1096 #endif 1097 #endif 1098 while(!*--xc) 1099 wa--; 1100 c->wds = wa; 1101 return c; 1102 } 1103 1104 static double 1105 ulp 1106 #ifdef KR_headers 1107 (x) double x; 1108 #else 1109 (double x) 1110 #endif 1111 { 1112 register Long L; 1113 double a; 1114 1115 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; 1116 #ifndef Avoid_Underflow 1117 #ifndef Sudden_Underflow 1118 if (L > 0) { 1119 #endif 1120 #endif 1121 #ifdef IBM 1122 L |= Exp_msk1 >> 4; 1123 #endif 1124 word0(a) = L; 1125 word1(a) = 0; 1126 #ifndef Avoid_Underflow 1127 #ifndef Sudden_Underflow 1128 } 1129 else { 1130 L = -L >> Exp_shift; 1131 if (L < Exp_shift) { 1132 word0(a) = 0x80000 >> L; 1133 word1(a) = 0; 1134 } 1135 else { 1136 word0(a) = 0; 1137 L -= Exp_shift; 1138 word1(a) = L >= 31 ? 1 : 1 << 31 - L; 1139 } 1140 } 1141 #endif 1142 #endif 1143 return dval(a); 1144 } 1145 1146 static double 1147 b2d 1148 #ifdef KR_headers 1149 (a, e) Bigint *a; int *e; 1150 #else 1151 (Bigint *a, int *e) 1152 #endif 1153 { 1154 ULong *xa, *xa0, w, y, z; 1155 int k; 1156 double d; 1157 #ifdef VAX 1158 ULong d0, d1; 1159 #else 1160 #define d0 word0(d) 1161 #define d1 word1(d) 1162 #endif 1163 1164 xa0 = a->x; 1165 xa = xa0 + a->wds; 1166 y = *--xa; 1167 #ifdef DEBUG 1168 if (!y) Bug("zero y in b2d"); 1169 #endif 1170 k = hi0bits(y); 1171 *e = 32 - k; 1172 #ifdef Pack_32 1173 if (k < Ebits) { 1174 d0 = Exp_1 | (y >> (Ebits - k)); 1175 w = xa > xa0 ? *--xa : 0; 1176 d1 = (y << ((32-Ebits) + k)) | (w >> (Ebits - k)); 1177 goto ret_d; 1178 } 1179 z = xa > xa0 ? *--xa : 0; 1180 if (k -= Ebits) { 1181 d0 = Exp_1 | (y << k) | (z >> (32 - k)); 1182 y = xa > xa0 ? *--xa : 0; 1183 d1 = (z << k) | (y >> (32 - k)); 1184 } 1185 else { 1186 d0 = Exp_1 | y; 1187 d1 = z; 1188 } 1189 #else 1190 if (k < Ebits + 16) { 1191 z = xa > xa0 ? *--xa : 0; 1192 d0 = Exp_1 | (y << (k - Ebits)) | (z >> (Ebits + 16 - k)); 1193 w = xa > xa0 ? *--xa : 0; 1194 y = xa > xa0 ? *--xa : 0; 1195 d1 = (z << (k + 16 - Ebits)) | (w << (k - Ebits)) | (y >> (16 + Ebits - k)); 1196 goto ret_d; 1197 } 1198 z = xa > xa0 ? *--xa : 0; 1199 w = xa > xa0 ? *--xa : 0; 1200 k -= Ebits + 16; 1201 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; 1202 y = xa > xa0 ? *--xa : 0; 1203 d1 = w << k + 16 | y << k; 1204 #endif 1205 ret_d: 1206 #ifdef VAX 1207 word0(d) = d0 >> 16 | d0 << 16; 1208 word1(d) = d1 >> 16 | d1 << 16; 1209 #else 1210 #undef d0 1211 #undef d1 1212 #endif 1213 return dval(d); 1214 } 1215 1216 static Bigint * 1217 d2b 1218 #ifdef KR_headers 1219 (d, e, bits) double d; int *e, *bits; 1220 #else 1221 (double d, int *e, int *bits) 1222 #endif 1223 { 1224 Bigint *b; 1225 int de, k; 1226 ULong *x, y, z; 1227 #ifndef Sudden_Underflow 1228 int i; 1229 #endif 1230 #ifdef VAX 1231 ULong d0, d1; 1232 d0 = word0(d) >> 16 | word0(d) << 16; 1233 d1 = word1(d) >> 16 | word1(d) << 16; 1234 #else 1235 #define d0 word0(d) 1236 #define d1 word1(d) 1237 #endif 1238 1239 #ifdef Pack_32 1240 b = Balloc(1); 1241 #else 1242 b = Balloc(2); 1243 #endif 1244 x = b->x; 1245 1246 z = d0 & Frac_mask; 1247 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 1248 #ifdef Sudden_Underflow 1249 de = (int)(d0 >> Exp_shift); 1250 #ifndef IBM 1251 z |= Exp_msk11; 1252 #endif 1253 #else 1254 if ((de = (int)(d0 >> Exp_shift))) 1255 z |= Exp_msk1; 1256 #endif 1257 #ifdef Pack_32 1258 if ((y = d1)) { 1259 if ((k = lo0bits(&y))) { 1260 x[0] = y | (z << (32 - k)); 1261 z >>= k; 1262 } 1263 else 1264 x[0] = y; 1265 #ifndef Sudden_Underflow 1266 i = 1267 #endif 1268 b->wds = (x[1] = z) ? 2 : 1; 1269 } 1270 else { 1271 /* This assertion fails for "1e-500" and other very 1272 * small numbers. It provides the right result (0) 1273 * though. This assert has also been removed from KJS's 1274 * version of dtoa.c. 1275 * 1276 * #ifdef DEBUG 1277 * if (!z) Bug("zero z in b2d"); 1278 * #endif 1279 */ 1280 k = lo0bits(&z); 1281 x[0] = z; 1282 #ifndef Sudden_Underflow 1283 i = 1284 #endif 1285 b->wds = 1; 1286 k += 32; 1287 } 1288 #else 1289 if (y = d1) { 1290 if (k = lo0bits(&y)) 1291 if (k >= 16) { 1292 x[0] = y | z << 32 - k & 0xffff; 1293 x[1] = z >> k - 16 & 0xffff; 1294 x[2] = z >> k; 1295 i = 2; 1296 } 1297 else { 1298 x[0] = y & 0xffff; 1299 x[1] = y >> 16 | z << 16 - k & 0xffff; 1300 x[2] = z >> k & 0xffff; 1301 x[3] = z >> k+16; 1302 i = 3; 1303 } 1304 else { 1305 x[0] = y & 0xffff; 1306 x[1] = y >> 16; 1307 x[2] = z & 0xffff; 1308 x[3] = z >> 16; 1309 i = 3; 1310 } 1311 } 1312 else { 1313 #ifdef DEBUG 1314 if (!z) 1315 Bug("Zero passed to d2b"); 1316 #endif 1317 k = lo0bits(&z); 1318 if (k >= 16) { 1319 x[0] = z; 1320 i = 0; 1321 } 1322 else { 1323 x[0] = z & 0xffff; 1324 x[1] = z >> 16; 1325 i = 1; 1326 } 1327 k += 32; 1328 } 1329 while(!x[i]) 1330 --i; 1331 b->wds = i + 1; 1332 #endif 1333 #ifndef Sudden_Underflow 1334 if (de) { 1335 #endif 1336 #ifdef IBM 1337 *e = (de - Bias - (P-1) << 2) + k; 1338 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask); 1339 #else 1340 *e = de - Bias - (P-1) + k; 1341 *bits = P - k; 1342 #endif 1343 #ifndef Sudden_Underflow 1344 } 1345 else { 1346 *e = de - Bias - (P-1) + 1 + k; 1347 #ifdef Pack_32 1348 *bits = 32*i - hi0bits(x[i-1]); 1349 #else 1350 *bits = (i+2)*16 - hi0bits(x[i]); 1351 #endif 1352 } 1353 #endif 1354 return b; 1355 } 1356 #undef d0 1357 #undef d1 1358 1359 static double 1360 ratio 1361 #ifdef KR_headers 1362 (a, b) Bigint *a, *b; 1363 #else 1364 (Bigint *a, Bigint *b) 1365 #endif 1366 { 1367 double da, db; 1368 int k, ka, kb; 1369 1370 dval(da) = b2d(a, &ka); 1371 dval(db) = b2d(b, &kb); 1372 #ifdef Pack_32 1373 k = ka - kb + 32*(a->wds - b->wds); 1374 #else 1375 k = ka - kb + 16*(a->wds - b->wds); 1376 #endif 1377 #ifdef IBM 1378 if (k > 0) { 1379 word0(da) += (k >> 2)*Exp_msk1; 1380 if (k &= 3) 1381 dval(da) *= 1 << k; 1382 } 1383 else { 1384 k = -k; 1385 word0(db) += (k >> 2)*Exp_msk1; 1386 if (k &= 3) 1387 dval(db) *= 1 << k; 1388 } 1389 #else 1390 if (k > 0) 1391 word0(da) += k*Exp_msk1; 1392 else { 1393 k = -k; 1394 word0(db) += k*Exp_msk1; 1395 } 1396 #endif 1397 return dval(da) / dval(db); 1398 } 1399 1400 static CONST double 1401 tens[] = { 1402 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1403 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1404 1e20, 1e21, 1e22 1405 #ifdef VAX 1406 , 1e23, 1e24 1407 #endif 1408 }; 1409 1410 static CONST double 1411 #ifdef IEEE_Arith 1412 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; 1413 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1414 #ifdef Avoid_Underflow 1415 9007199254740992.*9007199254740992.e-256 1416 /* = 2^106 * 1e-53 */ 1417 #else 1418 1e-256 1419 #endif 1420 }; 1421 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */ 1422 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */ 1423 #define Scale_Bit 0x10 1424 #define n_bigtens 5 1425 #else 1426 #ifdef IBM 1427 bigtens[] = { 1e16, 1e32, 1e64 }; 1428 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 }; 1429 #define n_bigtens 3 1430 #else 1431 bigtens[] = { 1e16, 1e32 }; 1432 static CONST double tinytens[] = { 1e-16, 1e-32 }; 1433 #define n_bigtens 2 1434 #endif 1435 #endif 1436 1437 #ifndef IEEE_Arith 1438 #undef INFNAN_CHECK 1439 #endif 1440 1441 #ifdef INFNAN_CHECK 1442 1443 #ifndef NAN_WORD0 1444 #define NAN_WORD0 0x7ff80000 1445 #endif 1446 1447 #ifndef NAN_WORD1 1448 #define NAN_WORD1 0 1449 #endif 1450 1451 static int 1452 match 1453 #ifdef KR_headers 1454 (sp, t) char **sp, *t; 1455 #else 1456 (CONST char **sp, char *t) 1457 #endif 1458 { 1459 int c, d; 1460 CONST char *s = *sp; 1461 1462 while(d = *t++) { 1463 if ((c = *++s) >= 'A' && c <= 'Z') 1464 c += 'a' - 'A'; 1465 if (c != d) 1466 return 0; 1467 } 1468 *sp = s + 1; 1469 return 1; 1470 } 1471 1472 #ifndef No_Hex_NaN 1473 static void 1474 hexnan 1475 #ifdef KR_headers 1476 (rvp, sp) double *rvp; CONST char **sp; 1477 #else 1478 (double *rvp, CONST char **sp) 1479 #endif 1480 { 1481 ULong c, x[2]; 1482 CONST char *s; 1483 int havedig, udx0, xshift; 1484 1485 x[0] = x[1] = 0; 1486 havedig = xshift = 0; 1487 udx0 = 1; 1488 s = *sp; 1489 while(c = *(CONST unsigned char*)++s) { 1490 if (c >= '0' && c <= '9') 1491 c -= '0'; 1492 else if (c >= 'a' && c <= 'f') 1493 c += 10 - 'a'; 1494 else if (c >= 'A' && c <= 'F') 1495 c += 10 - 'A'; 1496 else if (c <= ' ') { 1497 if (udx0 && havedig) { 1498 udx0 = 0; 1499 xshift = 1; 1500 } 1501 continue; 1502 } 1503 else if (/*(*/ c == ')' && havedig) { 1504 *sp = s + 1; 1505 break; 1506 } 1507 else 1508 return; /* invalid form: don't change *sp */ 1509 havedig = 1; 1510 if (xshift) { 1511 xshift = 0; 1512 x[0] = x[1]; 1513 x[1] = 0; 1514 } 1515 if (udx0) 1516 x[0] = (x[0] << 4) | (x[1] >> 28); 1517 x[1] = (x[1] << 4) | c; 1518 } 1519 if ((x[0] &= 0xfffff) || x[1]) { 1520 word0(*rvp) = Exp_mask | x[0]; 1521 word1(*rvp) = x[1]; 1522 } 1523 } 1524 #endif /*No_Hex_NaN*/ 1525 #endif /* INFNAN_CHECK */ 1526 1527 double 1528 strtod 1529 #ifdef KR_headers 1530 (s00, se) CONST char *s00; char **se; 1531 #else 1532 (CONST char *s00, char **se) 1533 #endif 1534 { 1535 #ifdef Avoid_Underflow 1536 int scale; 1537 #endif 1538 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, 1539 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; 1540 CONST char *s, *s0, *s1; 1541 double aadj, aadj1, adj, rv, rv0; 1542 Long L; 1543 ULong y, z; 1544 Bigint *bb = NULL, *bb1, *bd = NULL, *bd0, *bs = NULL, *delta = NULL; 1545 #ifdef SET_INEXACT 1546 int inexact, oldinexact; 1547 #endif 1548 #ifdef Honor_FLT_ROUNDS 1549 int rounding; 1550 #endif 1551 #ifdef USE_LOCALE 1552 CONST char *s2; 1553 #endif 1554 1555 sign = nz0 = nz = 0; 1556 dval(rv) = 0.; 1557 for(s = s00;;s++) switch(*s) { 1558 case '-': 1559 sign = 1; 1560 /* no break */ 1561 case '+': 1562 if (*++s) 1563 goto break2; 1564 /* no break */ 1565 case 0: 1566 goto ret0; 1567 case '\t': 1568 case '\n': 1569 case '\v': 1570 case '\f': 1571 case '\r': 1572 case ' ': 1573 continue; 1574 default: 1575 goto break2; 1576 } 1577 break2: 1578 if (*s == '0') { 1579 nz0 = 1; 1580 while(*++s == '0') ; 1581 if (!*s) 1582 goto ret; 1583 } 1584 s0 = s; 1585 y = z = 0; 1586 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 1587 if (nd < 9) 1588 y = 10*y + c - '0'; 1589 else if (nd < 16) 1590 z = 10*z + c - '0'; 1591 nd0 = nd; 1592 #ifdef USE_LOCALE 1593 s1 = localeconv()->decimal_point; 1594 if (c == *s1) { 1595 c = '.'; 1596 if (*++s1) { 1597 s2 = s; 1598 for(;;) { 1599 if (*++s2 != *s1) { 1600 c = 0; 1601 break; 1602 } 1603 if (!*++s1) { 1604 s = s2; 1605 break; 1606 } 1607 } 1608 } 1609 } 1610 #endif 1611 if (c == '.') { 1612 c = *++s; 1613 if (!nd) { 1614 for(; c == '0'; c = *++s) 1615 nz++; 1616 if (c > '0' && c <= '9') { 1617 s0 = s; 1618 nf += nz; 1619 nz = 0; 1620 goto have_dig; 1621 } 1622 goto dig_done; 1623 } 1624 for(; c >= '0' && c <= '9'; c = *++s) { 1625 have_dig: 1626 nz++; 1627 if (c -= '0') { 1628 nf += nz; 1629 for(i = 1; i < nz; i++) 1630 if (nd++ < 9) 1631 y *= 10; 1632 else if (nd <= DBL_DIG + 1) 1633 z *= 10; 1634 if (nd++ < 9) 1635 y = 10*y + c; 1636 else if (nd <= DBL_DIG + 1) 1637 z = 10*z + c; 1638 nz = 0; 1639 } 1640 } 1641 } 1642 dig_done: 1643 e = 0; 1644 if (c == 'e' || c == 'E') { 1645 if (!nd && !nz && !nz0) { 1646 goto ret0; 1647 } 1648 s00 = s; 1649 esign = 0; 1650 switch(c = *++s) { 1651 case '-': 1652 esign = 1; 1653 case '+': 1654 c = *++s; 1655 } 1656 if (c >= '0' && c <= '9') { 1657 while(c == '0') 1658 c = *++s; 1659 if (c > '0' && c <= '9') { 1660 L = c - '0'; 1661 s1 = s; 1662 while((c = *++s) >= '0' && c <= '9') 1663 L = 10*L + c - '0'; 1664 if (s - s1 > 8 || L > 19999) 1665 /* Avoid confusion from exponents 1666 * so large that e might overflow. 1667 */ 1668 e = 19999; /* safe for 16 bit ints */ 1669 else 1670 e = (int)L; 1671 if (esign) 1672 e = -e; 1673 } 1674 else 1675 e = 0; 1676 } 1677 else 1678 s = s00; 1679 } 1680 if (!nd) { 1681 if (!nz && !nz0) { 1682 #ifdef INFNAN_CHECK 1683 /* Check for Nan and Infinity */ 1684 switch(c) { 1685 case 'i': 1686 case 'I': 1687 if (match(&s,"nf")) { 1688 --s; 1689 if (!match(&s,"inity")) 1690 ++s; 1691 word0(rv) = 0x7ff00000; 1692 word1(rv) = 0; 1693 goto ret; 1694 } 1695 break; 1696 case 'n': 1697 case 'N': 1698 if (match(&s, "an")) { 1699 word0(rv) = NAN_WORD0; 1700 word1(rv) = NAN_WORD1; 1701 #ifndef No_Hex_NaN 1702 if (*s == '(') /*)*/ 1703 hexnan(&rv, &s); 1704 #endif 1705 goto ret; 1706 } 1707 } 1708 #endif /* INFNAN_CHECK */ 1709 ret0: 1710 s = s00; 1711 sign = 0; 1712 } 1713 goto ret; 1714 } 1715 e1 = e -= nf; 1716 1717 /* Now we have nd0 digits, starting at s0, followed by a 1718 * decimal point, followed by nd-nd0 digits. The number we're 1719 * after is the integer represented by those digits times 1720 * 10**e */ 1721 1722 if (!nd0) 1723 nd0 = nd; 1724 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 1725 dval(rv) = y; 1726 if (k > 9) { 1727 #ifdef SET_INEXACT 1728 if (k > DBL_DIG) 1729 oldinexact = get_inexact(); 1730 #endif 1731 dval(rv) = tens[k - 9] * dval(rv) + z; 1732 } 1733 bd0 = 0; 1734 if (nd <= DBL_DIG 1735 #ifndef RND_PRODQUOT 1736 #ifndef Honor_FLT_ROUNDS 1737 && Flt_Rounds == 1 1738 #endif 1739 #endif 1740 ) { 1741 if (!e) 1742 goto ret; 1743 if (e > 0) { 1744 if (e <= Ten_pmax) { 1745 #ifdef VAX 1746 goto vax_ovfl_check; 1747 #else 1748 #ifdef Honor_FLT_ROUNDS 1749 /* round correctly FLT_ROUNDS = 2 or 3 */ 1750 if (sign) { 1751 rv = -rv; 1752 sign = 0; 1753 } 1754 #endif 1755 /* rv = */ rounded_product(dval(rv), tens[e]); 1756 goto ret; 1757 #endif 1758 } 1759 i = DBL_DIG - nd; 1760 if (e <= Ten_pmax + i) { 1761 /* A fancier test would sometimes let us do 1762 * this for larger i values. 1763 */ 1764 #ifdef Honor_FLT_ROUNDS 1765 /* round correctly FLT_ROUNDS = 2 or 3 */ 1766 if (sign) { 1767 rv = -rv; 1768 sign = 0; 1769 } 1770 #endif 1771 e -= i; 1772 dval(rv) *= tens[i]; 1773 #ifdef VAX 1774 /* VAX exponent range is so narrow we must 1775 * worry about overflow here... 1776 */ 1777 vax_ovfl_check: 1778 word0(rv) -= P*Exp_msk1; 1779 /* rv = */ rounded_product(dval(rv), tens[e]); 1780 if ((word0(rv) & Exp_mask) 1781 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 1782 goto ovfl; 1783 word0(rv) += P*Exp_msk1; 1784 #else 1785 /* rv = */ rounded_product(dval(rv), tens[e]); 1786 #endif 1787 goto ret; 1788 } 1789 } 1790 #ifndef Inaccurate_Divide 1791 else if (e >= -Ten_pmax) { 1792 #ifdef Honor_FLT_ROUNDS 1793 /* round correctly FLT_ROUNDS = 2 or 3 */ 1794 if (sign) { 1795 rv = -rv; 1796 sign = 0; 1797 } 1798 #endif 1799 /* rv = */ rounded_quotient(dval(rv), tens[-e]); 1800 goto ret; 1801 } 1802 #endif 1803 } 1804 e1 += nd - k; 1805 1806 #ifdef IEEE_Arith 1807 #ifdef SET_INEXACT 1808 inexact = 1; 1809 if (k <= DBL_DIG) 1810 oldinexact = get_inexact(); 1811 #endif 1812 #ifdef Avoid_Underflow 1813 scale = 0; 1814 #endif 1815 #ifdef Honor_FLT_ROUNDS 1816 if ((rounding = Flt_Rounds) >= 2) { 1817 if (sign) 1818 rounding = rounding == 2 ? 0 : 2; 1819 else 1820 if (rounding != 2) 1821 rounding = 0; 1822 } 1823 #endif 1824 #endif /*IEEE_Arith*/ 1825 1826 /* Get starting approximation = rv * 10**e1 */ 1827 1828 if (e1 > 0) { 1829 if ((i = e1 & 15)) 1830 dval(rv) *= tens[i]; 1831 if (e1 &= ~15) { 1832 if (e1 > DBL_MAX_10_EXP) { 1833 ovfl: 1834 #ifndef NO_ERRNO 1835 errno = ERANGE; 1836 #endif 1837 /* Can't trust HUGE_VAL */ 1838 #ifdef IEEE_Arith 1839 #ifdef Honor_FLT_ROUNDS 1840 switch(rounding) { 1841 case 0: /* toward 0 */ 1842 case 3: /* toward -infinity */ 1843 word0(rv) = Big0; 1844 word1(rv) = Big1; 1845 break; 1846 default: 1847 word0(rv) = Exp_mask; 1848 word1(rv) = 0; 1849 } 1850 #else /*Honor_FLT_ROUNDS*/ 1851 word0(rv) = Exp_mask; 1852 word1(rv) = 0; 1853 #endif /*Honor_FLT_ROUNDS*/ 1854 #ifdef SET_INEXACT 1855 /* set overflow bit */ 1856 dval(rv0) = 1e300; 1857 dval(rv0) *= dval(rv0); 1858 #endif 1859 #else /*IEEE_Arith*/ 1860 word0(rv) = Big0; 1861 word1(rv) = Big1; 1862 #endif /*IEEE_Arith*/ 1863 if (bd0) 1864 goto retfree; 1865 goto ret; 1866 } 1867 e1 >>= 4; 1868 for(j = 0; e1 > 1; j++, e1 >>= 1) 1869 if (e1 & 1) 1870 dval(rv) *= bigtens[j]; 1871 /* The last multiplication could overflow. */ 1872 word0(rv) -= P*Exp_msk1; 1873 dval(rv) *= bigtens[j]; 1874 if ((z = word0(rv) & Exp_mask) 1875 > Exp_msk1*(DBL_MAX_EXP+Bias-P)) 1876 goto ovfl; 1877 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { 1878 /* set to largest number */ 1879 /* (Can't trust DBL_MAX) */ 1880 word0(rv) = Big0; 1881 word1(rv) = Big1; 1882 } 1883 else 1884 word0(rv) += P*Exp_msk1; 1885 } 1886 } 1887 else if (e1 < 0) { 1888 e1 = -e1; 1889 if ((i = e1 & 15)) 1890 dval(rv) /= tens[i]; 1891 if (e1 >>= 4) { 1892 if (e1 >= 1 << n_bigtens) 1893 goto undfl; 1894 #ifdef Avoid_Underflow 1895 if (e1 & Scale_Bit) 1896 scale = 2*P; 1897 for(j = 0; e1 > 0; j++, e1 >>= 1) 1898 if (e1 & 1) 1899 dval(rv) *= tinytens[j]; 1900 if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask) 1901 >> Exp_shift)) > 0) { 1902 /* scaled rv is denormal; zap j low bits */ 1903 if (j >= 32) { 1904 word1(rv) = 0; 1905 if (j >= 53) 1906 word0(rv) = (P+2)*Exp_msk1; 1907 else 1908 word0(rv) &= 0xffffffff << (j-32); 1909 } 1910 else 1911 word1(rv) &= 0xffffffff << j; 1912 } 1913 #else 1914 for(j = 0; e1 > 1; j++, e1 >>= 1) 1915 if (e1 & 1) 1916 dval(rv) *= tinytens[j]; 1917 /* The last multiplication could underflow. */ 1918 dval(rv0) = dval(rv); 1919 dval(rv) *= tinytens[j]; 1920 if (!dval(rv)) { 1921 dval(rv) = 2.*dval(rv0); 1922 dval(rv) *= tinytens[j]; 1923 #endif 1924 if (!dval(rv)) { 1925 undfl: 1926 dval(rv) = 0.; 1927 #ifndef NO_ERRNO 1928 errno = ERANGE; 1929 #endif 1930 if (bd0) 1931 goto retfree; 1932 goto ret; 1933 } 1934 #ifndef Avoid_Underflow 1935 word0(rv) = Tiny0; 1936 word1(rv) = Tiny1; 1937 /* The refinement below will clean 1938 * this approximation up. 1939 */ 1940 } 1941 #endif 1942 } 1943 } 1944 1945 /* Now the hard part -- adjusting rv to the correct value.*/ 1946 1947 /* Put digits into bd: true value = bd * 10^e */ 1948 1949 bd0 = s2b(s0, nd0, nd, y); 1950 1951 for(;;) { 1952 bd = Balloc(bd0->k); 1953 Bcopy(bd, bd0); 1954 bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */ 1955 bs = i2b(1); 1956 1957 if (e >= 0) { 1958 bb2 = bb5 = 0; 1959 bd2 = bd5 = e; 1960 } 1961 else { 1962 bb2 = bb5 = -e; 1963 bd2 = bd5 = 0; 1964 } 1965 if (bbe >= 0) 1966 bb2 += bbe; 1967 else 1968 bd2 -= bbe; 1969 bs2 = bb2; 1970 #ifdef Honor_FLT_ROUNDS 1971 if (rounding != 1) 1972 bs2++; 1973 #endif 1974 #ifdef Avoid_Underflow 1975 j = bbe - scale; 1976 i = j + bbbits - 1; /* logb(rv) */ 1977 if (i < Emin) /* denormal */ 1978 j += P - Emin; 1979 else 1980 j = P + 1 - bbbits; 1981 #else /*Avoid_Underflow*/ 1982 #ifdef Sudden_Underflow 1983 #ifdef IBM 1984 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 1985 #else 1986 j = P + 1 - bbbits; 1987 #endif 1988 #else /*Sudden_Underflow*/ 1989 j = bbe; 1990 i = j + bbbits - 1; /* logb(rv) */ 1991 if (i < Emin) /* denormal */ 1992 j += P - Emin; 1993 else 1994 j = P + 1 - bbbits; 1995 #endif /*Sudden_Underflow*/ 1996 #endif /*Avoid_Underflow*/ 1997 bb2 += j; 1998 bd2 += j; 1999 #ifdef Avoid_Underflow 2000 bd2 += scale; 2001 #endif 2002 i = bb2 < bd2 ? bb2 : bd2; 2003 if (i > bs2) 2004 i = bs2; 2005 if (i > 0) { 2006 bb2 -= i; 2007 bd2 -= i; 2008 bs2 -= i; 2009 } 2010 if (bb5 > 0) { 2011 bs = pow5mult(bs, bb5); 2012 bb1 = mult(bs, bb); 2013 Bfree(bb); 2014 bb = bb1; 2015 } 2016 if (bb2 > 0) 2017 bb = lshift(bb, bb2); 2018 if (bd5 > 0) 2019 bd = pow5mult(bd, bd5); 2020 if (bd2 > 0) 2021 bd = lshift(bd, bd2); 2022 if (bs2 > 0) 2023 bs = lshift(bs, bs2); 2024 delta = diff(bb, bd); 2025 dsign = delta->sign; 2026 delta->sign = 0; 2027 i = cmp(delta, bs); 2028 #ifdef Honor_FLT_ROUNDS 2029 if (rounding != 1) { 2030 if (i < 0) { 2031 /* Error is less than an ulp */ 2032 if (!delta->x[0] && delta->wds <= 1) { 2033 /* exact */ 2034 #ifdef SET_INEXACT 2035 inexact = 0; 2036 #endif 2037 break; 2038 } 2039 if (rounding) { 2040 if (dsign) { 2041 adj = 1.; 2042 goto apply_adj; 2043 } 2044 } 2045 else if (!dsign) { 2046 adj = -1.; 2047 if (!word1(rv) 2048 && !(word0(rv) & Frac_mask)) { 2049 y = word0(rv) & Exp_mask; 2050 #ifdef Avoid_Underflow 2051 if (!scale || y > 2*P*Exp_msk1) 2052 #else 2053 if (y) 2054 #endif 2055 { 2056 delta = lshift(delta,Log2P); 2057 if (cmp(delta, bs) <= 0) 2058 adj = -0.5; 2059 } 2060 } 2061 apply_adj: 2062 #ifdef Avoid_Underflow 2063 if (scale && (y = word0(rv) & Exp_mask) 2064 <= 2*P*Exp_msk1) 2065 word0(adj) += (2*P+1)*Exp_msk1 - y; 2066 #else 2067 #ifdef Sudden_Underflow 2068 if ((word0(rv) & Exp_mask) <= 2069 P*Exp_msk1) { 2070 word0(rv) += P*Exp_msk1; 2071 dval(rv) += adj*ulp(dval(rv)); 2072 word0(rv) -= P*Exp_msk1; 2073 } 2074 else 2075 #endif /*Sudden_Underflow*/ 2076 #endif /*Avoid_Underflow*/ 2077 dval(rv) += adj*ulp(dval(rv)); 2078 } 2079 break; 2080 } 2081 adj = ratio(delta, bs); 2082 if (adj < 1.) 2083 adj = 1.; 2084 if (adj <= 0x7ffffffe) { 2085 /* adj = rounding ? ceil(adj) : floor(adj); */ 2086 y = adj; 2087 if (y != adj) { 2088 if (!((rounding>>1) ^ dsign)) 2089 y++; 2090 adj = y; 2091 } 2092 } 2093 #ifdef Avoid_Underflow 2094 if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1) 2095 word0(adj) += (2*P+1)*Exp_msk1 - y; 2096 #else 2097 #ifdef Sudden_Underflow 2098 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 2099 word0(rv) += P*Exp_msk1; 2100 adj *= ulp(dval(rv)); 2101 if (dsign) 2102 dval(rv) += adj; 2103 else 2104 dval(rv) -= adj; 2105 word0(rv) -= P*Exp_msk1; 2106 goto cont; 2107 } 2108 #endif /*Sudden_Underflow*/ 2109 #endif /*Avoid_Underflow*/ 2110 adj *= ulp(dval(rv)); 2111 if (dsign) 2112 dval(rv) += adj; 2113 else 2114 dval(rv) -= adj; 2115 goto cont; 2116 } 2117 #endif /*Honor_FLT_ROUNDS*/ 2118 2119 if (i < 0) { 2120 /* Error is less than half an ulp -- check for 2121 * special case of mantissa a power of two. 2122 */ 2123 if (dsign || word1(rv) || word0(rv) & Bndry_mask 2124 #ifdef IEEE_Arith 2125 #ifdef Avoid_Underflow 2126 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1 2127 #else 2128 || (word0(rv) & Exp_mask) <= Exp_msk1 2129 #endif 2130 #endif 2131 ) { 2132 #ifdef SET_INEXACT 2133 if (!delta->x[0] && delta->wds <= 1) 2134 inexact = 0; 2135 #endif 2136 break; 2137 } 2138 if (!delta->x[0] && delta->wds <= 1) { 2139 /* exact result */ 2140 #ifdef SET_INEXACT 2141 inexact = 0; 2142 #endif 2143 break; 2144 } 2145 delta = lshift(delta,Log2P); 2146 if (cmp(delta, bs) > 0) 2147 goto drop_down; 2148 break; 2149 } 2150 if (i == 0) { 2151 /* exactly half-way between */ 2152 if (dsign) { 2153 if ((word0(rv) & Bndry_mask1) == Bndry_mask1 2154 && word1(rv) == ( 2155 #ifdef Avoid_Underflow 2156 (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1) 2157 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : 2158 #endif 2159 0xffffffff)) { 2160 /*boundary case -- increment exponent*/ 2161 word0(rv) = (word0(rv) & Exp_mask) 2162 + Exp_msk1 2163 #ifdef IBM 2164 | Exp_msk1 >> 4 2165 #endif 2166 ; 2167 word1(rv) = 0; 2168 #ifdef Avoid_Underflow 2169 dsign = 0; 2170 #endif 2171 break; 2172 } 2173 } 2174 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { 2175 drop_down: 2176 /* boundary case -- decrement exponent */ 2177 #ifdef Sudden_Underflow /*{{*/ 2178 L = word0(rv) & Exp_mask; 2179 #ifdef IBM 2180 if (L < Exp_msk1) 2181 #else 2182 #ifdef Avoid_Underflow 2183 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1)) 2184 #else 2185 if (L <= Exp_msk1) 2186 #endif /*Avoid_Underflow*/ 2187 #endif /*IBM*/ 2188 goto undfl; 2189 L -= Exp_msk1; 2190 #else /*Sudden_Underflow}{*/ 2191 #ifdef Avoid_Underflow 2192 if (scale) { 2193 L = word0(rv) & Exp_mask; 2194 if (L <= (2*P+1)*Exp_msk1) { 2195 if (L > (P+2)*Exp_msk1) 2196 /* round even ==> */ 2197 /* accept rv */ 2198 break; 2199 /* rv = smallest denormal */ 2200 goto undfl; 2201 } 2202 } 2203 #endif /*Avoid_Underflow*/ 2204 L = (word0(rv) & Exp_mask) - Exp_msk1; 2205 #endif /*Sudden_Underflow}}*/ 2206 word0(rv) = L | Bndry_mask1; 2207 word1(rv) = 0xffffffff; 2208 #ifdef IBM 2209 goto cont; 2210 #else 2211 break; 2212 #endif 2213 } 2214 #ifndef ROUND_BIASED 2215 if (!(word1(rv) & LSB)) 2216 break; 2217 #endif 2218 if (dsign) 2219 dval(rv) += ulp(dval(rv)); 2220 #ifndef ROUND_BIASED 2221 else { 2222 dval(rv) -= ulp(dval(rv)); 2223 #ifndef Sudden_Underflow 2224 if (!dval(rv)) 2225 goto undfl; 2226 #endif 2227 } 2228 #ifdef Avoid_Underflow 2229 dsign = 1 - dsign; 2230 #endif 2231 #endif 2232 break; 2233 } 2234 if ((aadj = ratio(delta, bs)) <= 2.) { 2235 if (dsign) 2236 aadj = aadj1 = 1.; 2237 else if (word1(rv) || word0(rv) & Bndry_mask) { 2238 #ifndef Sudden_Underflow 2239 if (word1(rv) == Tiny1 && !word0(rv)) 2240 goto undfl; 2241 #endif 2242 aadj = 1.; 2243 aadj1 = -1.; 2244 } 2245 else { 2246 /* special case -- power of FLT_RADIX to be */ 2247 /* rounded down... */ 2248 2249 if (aadj < 2./FLT_RADIX) 2250 aadj = 1./FLT_RADIX; 2251 else 2252 aadj *= 0.5; 2253 aadj1 = -aadj; 2254 } 2255 } 2256 else { 2257 aadj *= 0.5; 2258 aadj1 = dsign ? aadj : -aadj; 2259 #ifdef Check_FLT_ROUNDS 2260 switch(Rounding) { 2261 case 2: /* towards +infinity */ 2262 aadj1 -= 0.5; 2263 break; 2264 case 0: /* towards 0 */ 2265 case 3: /* towards -infinity */ 2266 aadj1 += 0.5; 2267 } 2268 #else 2269 if (Flt_Rounds == 0) 2270 aadj1 += 0.5; 2271 #endif /*Check_FLT_ROUNDS*/ 2272 } 2273 y = word0(rv) & Exp_mask; 2274 2275 /* Check for overflow */ 2276 2277 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 2278 dval(rv0) = dval(rv); 2279 word0(rv) -= P*Exp_msk1; 2280 adj = aadj1 * ulp(dval(rv)); 2281 dval(rv) += adj; 2282 if ((word0(rv) & Exp_mask) >= 2283 Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 2284 if (word0(rv0) == Big0 && word1(rv0) == Big1) 2285 goto ovfl; 2286 word0(rv) = Big0; 2287 word1(rv) = Big1; 2288 goto cont; 2289 } 2290 else 2291 word0(rv) += P*Exp_msk1; 2292 } 2293 else { 2294 #ifdef Avoid_Underflow 2295 if (scale && y <= 2*P*Exp_msk1) { 2296 if (aadj <= 0x7fffffff) { 2297 if ((z = aadj) <= 0) 2298 z = 1; 2299 aadj = z; 2300 aadj1 = dsign ? aadj : -aadj; 2301 } 2302 word0(aadj1) += (2*P+1)*Exp_msk1 - y; 2303 } 2304 adj = aadj1 * ulp(dval(rv)); 2305 dval(rv) += adj; 2306 #else 2307 #ifdef Sudden_Underflow 2308 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 2309 dval(rv0) = dval(rv); 2310 word0(rv) += P*Exp_msk1; 2311 adj = aadj1 * ulp(dval(rv)); 2312 dval(rv) += adj; 2313 #ifdef IBM 2314 if ((word0(rv) & Exp_mask) < P*Exp_msk1) 2315 #else 2316 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) 2317 #endif 2318 { 2319 if (word0(rv0) == Tiny0 2320 && word1(rv0) == Tiny1) 2321 goto undfl; 2322 word0(rv) = Tiny0; 2323 word1(rv) = Tiny1; 2324 goto cont; 2325 } 2326 else 2327 word0(rv) -= P*Exp_msk1; 2328 } 2329 else { 2330 adj = aadj1 * ulp(dval(rv)); 2331 dval(rv) += adj; 2332 } 2333 #else /*Sudden_Underflow*/ 2334 /* Compute adj so that the IEEE rounding rules will 2335 * correctly round rv + adj in some half-way cases. 2336 * If rv * ulp(rv) is denormalized (i.e., 2337 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 2338 * trouble from bits lost to denormalization; 2339 * example: 1.2e-307 . 2340 */ 2341 if (y <= (P-1)*Exp_msk1 && aadj > 1.) { 2342 aadj1 = (double)(int)(aadj + 0.5); 2343 if (!dsign) 2344 aadj1 = -aadj1; 2345 } 2346 adj = aadj1 * ulp(dval(rv)); 2347 dval(rv) += adj; 2348 #endif /*Sudden_Underflow*/ 2349 #endif /*Avoid_Underflow*/ 2350 } 2351 z = word0(rv) & Exp_mask; 2352 #ifndef SET_INEXACT 2353 #ifdef Avoid_Underflow 2354 if (!scale) 2355 #endif 2356 if (y == z) { 2357 /* Can we stop now? */ 2358 L = (Long)aadj; 2359 aadj -= L; 2360 /* The tolerances below are conservative. */ 2361 if (dsign || word1(rv) || word0(rv) & Bndry_mask) { 2362 if (aadj < .4999999 || aadj > .5000001) 2363 break; 2364 } 2365 else if (aadj < .4999999/FLT_RADIX) 2366 break; 2367 } 2368 #endif 2369 cont: 2370 Bfree(bb); 2371 Bfree(bd); 2372 Bfree(bs); 2373 Bfree(delta); 2374 } 2375 #ifdef SET_INEXACT 2376 if (inexact) { 2377 if (!oldinexact) { 2378 word0(rv0) = Exp_1 + (70 << Exp_shift); 2379 word1(rv0) = 0; 2380 dval(rv0) += 1.; 2381 } 2382 } 2383 else if (!oldinexact) 2384 clear_inexact(); 2385 #endif 2386 #ifdef Avoid_Underflow 2387 if (scale) { 2388 word0(rv0) = Exp_1 - 2*P*Exp_msk1; 2389 word1(rv0) = 0; 2390 dval(rv) *= dval(rv0); 2391 #ifndef NO_ERRNO 2392 /* try to avoid the bug of testing an 8087 register value */ 2393 if (word0(rv) == 0 && word1(rv) == 0) 2394 errno = ERANGE; 2395 #endif 2396 } 2397 #endif /* Avoid_Underflow */ 2398 #ifdef SET_INEXACT 2399 if (inexact && !(word0(rv) & Exp_mask)) { 2400 /* set underflow bit */ 2401 dval(rv0) = 1e-300; 2402 dval(rv0) *= dval(rv0); 2403 } 2404 #endif 2405 retfree: 2406 Bfree(bb); 2407 Bfree(bd); 2408 Bfree(bs); 2409 Bfree(bd0); 2410 Bfree(delta); 2411 ret: 2412 if (se) 2413 *se = (char *)s; 2414 return sign ? -dval(rv) : dval(rv); 2415 } 2416 2417 static int 2418 quorem 2419 #ifdef KR_headers 2420 (b, S) Bigint *b, *S; 2421 #else 2422 (Bigint *b, Bigint *S) 2423 #endif 2424 { 2425 int n; 2426 ULong *bx, *bxe, q, *sx, *sxe; 2427 #ifdef ULLong 2428 ULLong borrow, carry, y, ys; 2429 #else 2430 ULong borrow, carry, y, ys; 2431 #ifdef Pack_32 2432 ULong si, z, zs; 2433 #endif 2434 #endif 2435 2436 n = S->wds; 2437 #ifdef DEBUG 2438 /*debug*/ if (b->wds > n) 2439 /*debug*/ Bug("oversize b in quorem"); 2440 #endif 2441 if (b->wds < n) 2442 return 0; 2443 sx = S->x; 2444 sxe = sx + --n; 2445 bx = b->x; 2446 bxe = bx + n; 2447 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ 2448 #ifdef DEBUG 2449 /*debug*/ if (q > 9) 2450 /*debug*/ Bug("oversized quotient in quorem"); 2451 #endif 2452 if (q) { 2453 borrow = 0; 2454 carry = 0; 2455 do { 2456 #ifdef ULLong 2457 ys = *sx++ * (ULLong)q + carry; 2458 carry = ys >> 32; 2459 y = *bx - (ys & FFFFFFFF) - borrow; 2460 borrow = y >> 32 & (ULong)1; 2461 *bx++ = y & FFFFFFFF; 2462 #else 2463 #ifdef Pack_32 2464 si = *sx++; 2465 ys = (si & 0xffff) * q + carry; 2466 zs = (si >> 16) * q + (ys >> 16); 2467 carry = zs >> 16; 2468 y = (*bx & 0xffff) - (ys & 0xffff) - borrow; 2469 borrow = (y & 0x10000) >> 16; 2470 z = (*bx >> 16) - (zs & 0xffff) - borrow; 2471 borrow = (z & 0x10000) >> 16; 2472 Storeinc(bx, z, y); 2473 #else 2474 ys = *sx++ * q + carry; 2475 carry = ys >> 16; 2476 y = *bx - (ys & 0xffff) - borrow; 2477 borrow = (y & 0x10000) >> 16; 2478 *bx++ = y & 0xffff; 2479 #endif 2480 #endif 2481 } 2482 while(sx <= sxe); 2483 if (!*bxe) { 2484 bx = b->x; 2485 while(--bxe > bx && !*bxe) 2486 --n; 2487 b->wds = n; 2488 } 2489 } 2490 if (cmp(b, S) >= 0) { 2491 q++; 2492 borrow = 0; 2493 carry = 0; 2494 bx = b->x; 2495 sx = S->x; 2496 do { 2497 #ifdef ULLong 2498 ys = *sx++ + carry; 2499 carry = ys >> 32; 2500 y = *bx - (ys & FFFFFFFF) - borrow; 2501 borrow = y >> 32 & (ULong)1; 2502 *bx++ = y & FFFFFFFF; 2503 #else 2504 #ifdef Pack_32 2505 si = *sx++; 2506 ys = (si & 0xffff) + carry; 2507 zs = (si >> 16) + (ys >> 16); 2508 carry = zs >> 16; 2509 y = (*bx & 0xffff) - (ys & 0xffff) - borrow; 2510 borrow = (y & 0x10000) >> 16; 2511 z = (*bx >> 16) - (zs & 0xffff) - borrow; 2512 borrow = (z & 0x10000) >> 16; 2513 Storeinc(bx, z, y); 2514 #else 2515 ys = *sx++ + carry; 2516 carry = ys >> 16; 2517 y = *bx - (ys & 0xffff) - borrow; 2518 borrow = (y & 0x10000) >> 16; 2519 *bx++ = y & 0xffff; 2520 #endif 2521 #endif 2522 } 2523 while(sx <= sxe); 2524 bx = b->x; 2525 bxe = bx + n; 2526 if (!*bxe) { 2527 while(--bxe > bx && !*bxe) 2528 --n; 2529 b->wds = n; 2530 } 2531 } 2532 return q; 2533 } 2534 2535 #ifndef MULTIPLE_THREADS 2536 static char *dtoa_result; 2537 #endif 2538 2539 static char * 2540 #ifdef KR_headers 2541 rv_alloc(i) int i; 2542 #else 2543 rv_alloc(int i) 2544 #endif 2545 { 2546 int j, k, *r; 2547 2548 j = sizeof(ULong); 2549 for(k = 0; 2550 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i; 2551 j <<= 1) 2552 k++; 2553 r = (int*)Balloc(k); 2554 *r = k; 2555 return 2556 #ifndef MULTIPLE_THREADS 2557 dtoa_result = 2558 #endif 2559 (char *)(r+1); 2560 } 2561 2562 static char * 2563 #ifdef KR_headers 2564 nrv_alloc(s, rve, n) char *s, **rve; int n; 2565 #else 2566 nrv_alloc(const char *s, char **rve, int n) 2567 #endif 2568 { 2569 char *rv, *t; 2570 2571 t = rv = rv_alloc(n); 2572 while ((*t = *s++)) t++; 2573 if (rve) 2574 *rve = t; 2575 return rv; 2576 } 2577 2578 /* freedtoa(s) must be used to free values s returned by dtoa 2579 * when MULTIPLE_THREADS is #defined. It should be used in all cases, 2580 * but for consistency with earlier versions of dtoa, it is optional 2581 * when MULTIPLE_THREADS is not defined. 2582 */ 2583 2584 void 2585 #ifdef KR_headers 2586 freedtoa(s) char *s; 2587 #else 2588 freedtoa(char *s) 2589 #endif 2590 { 2591 Bigint *b = (Bigint *)((int *)s - 1); 2592 b->maxwds = 1 << (b->k = *(int*)b); 2593 Bfree(b); 2594 #ifndef MULTIPLE_THREADS 2595 if (s == dtoa_result) 2596 dtoa_result = 0; 2597 #endif 2598 } 2599 2600 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. 2601 * 2602 * Inspired by "How to Print Floating-Point Numbers Accurately" by 2603 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126]. 2604 * 2605 * Modifications: 2606 * 1. Rather than iterating, we use a simple numeric overestimate 2607 * to determine k = floor(log10(d)). We scale relevant 2608 * quantities using O(log2(k)) rather than O(k) multiplications. 2609 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't 2610 * try to generate digits strictly left to right. Instead, we 2611 * compute with fewer bits and propagate the carry if necessary 2612 * when rounding the final digit up. This is often faster. 2613 * 3. Under the assumption that input will be rounded nearest, 2614 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. 2615 * That is, we allow equality in stopping tests when the 2616 * round-nearest rule will give the same floating-point value 2617 * as would satisfaction of the stopping test with strict 2618 * inequality. 2619 * 4. We remove common factors of powers of 2 from relevant 2620 * quantities. 2621 * 5. When converting floating-point integers less than 1e16, 2622 * we use floating-point arithmetic rather than resorting 2623 * to multiple-precision integers. 2624 * 6. When asked to produce fewer than 15 digits, we first try 2625 * to get by with floating-point arithmetic; we resort to 2626 * multiple-precision integer arithmetic only if we cannot 2627 * guarantee that the floating-point calculation has given 2628 * the correctly rounded result. For k requested digits and 2629 * "uniformly" distributed input, the probability is 2630 * something like 10^(k-15) that we must resort to the Long 2631 * calculation. 2632 */ 2633 2634 char * 2635 dtoa 2636 #ifdef KR_headers 2637 (d, mode, ndigits, decpt, sign, rve) 2638 double d; int mode, ndigits, *decpt, *sign; char **rve; 2639 #else 2640 (double d, int mode, int ndigits, int *decpt, int *sign, char **rve) 2641 #endif 2642 { 2643 /* Arguments ndigits, decpt, sign are similar to those 2644 of ecvt and fcvt; trailing zeros are suppressed from 2645 the returned string. If not null, *rve is set to point 2646 to the end of the return value. If d is +-Infinity or NaN, 2647 then *decpt is set to 9999. 2648 2649 mode: 2650 0 ==> shortest string that yields d when read in 2651 and rounded to nearest. 2652 1 ==> like 0, but with Steele & White stopping rule; 2653 e.g. with IEEE P754 arithmetic , mode 0 gives 2654 1e23 whereas mode 1 gives 9.999999999999999e22. 2655 2 ==> max(1,ndigits) significant digits. This gives a 2656 return value similar to that of ecvt, except 2657 that trailing zeros are suppressed. 2658 3 ==> through ndigits past the decimal point. This 2659 gives a return value similar to that from fcvt, 2660 except that trailing zeros are suppressed, and 2661 ndigits can be negative. 2662 4,5 ==> similar to 2 and 3, respectively, but (in 2663 round-nearest mode) with the tests of mode 0 to 2664 possibly return a shorter string that rounds to d. 2665 With IEEE arithmetic and compilation with 2666 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same 2667 as modes 2 and 3 when FLT_ROUNDS != 1. 2668 6-9 ==> Debugging modes similar to mode - 4: don't try 2669 fast floating-point estimate (if applicable). 2670 2671 Values of mode other than 0-9 are treated as mode 0. 2672 2673 Sufficient space is allocated to the return value 2674 to hold the suppressed trailing zeros. 2675 */ 2676 2677 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1, 2678 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5, 2679 spec_case, try_quick, bias_round_up; 2680 Long L; 2681 #ifndef Sudden_Underflow 2682 int denorm; 2683 ULong x; 2684 #endif 2685 Bigint *b, *b1, *delta, *mlo, *mhi, *S; 2686 double d2, ds, eps; 2687 char *s, *s0; 2688 #ifdef Honor_FLT_ROUNDS 2689 int rounding; 2690 #endif 2691 #ifdef SET_INEXACT 2692 int inexact, oldinexact; 2693 #endif 2694 2695 /* In mode 2 and 3 we bias rounding up when there are ties. */ 2696 bias_round_up = mode == 2 || mode == 3; 2697 2698 ilim = ilim1 = 0; /* to avoid Google3 compiler warnings */ 2699 2700 #ifndef MULTIPLE_THREADS 2701 if (dtoa_result) { 2702 freedtoa(dtoa_result); 2703 dtoa_result = 0; 2704 } 2705 #endif 2706 2707 if (word0(d) & Sign_bit) { 2708 /* set sign for everything, including 0's and NaNs */ 2709 *sign = 1; 2710 word0(d) &= ~Sign_bit; /* clear sign bit */ 2711 } 2712 else 2713 *sign = 0; 2714 2715 #if defined(IEEE_Arith) + defined(VAX) 2716 #ifdef IEEE_Arith 2717 if ((word0(d) & Exp_mask) == Exp_mask) 2718 #else 2719 if (word0(d) == 0x8000) 2720 #endif 2721 { 2722 /* Infinity or NaN */ 2723 *decpt = 9999; 2724 #ifdef IEEE_Arith 2725 if (!word1(d) && !(word0(d) & 0xfffff)) 2726 return nrv_alloc("Infinity", rve, 8); 2727 #endif 2728 return nrv_alloc("NaN", rve, 3); 2729 } 2730 #endif 2731 #ifdef IBM 2732 dval(d) += 0; /* normalize */ 2733 #endif 2734 if (!dval(d)) { 2735 *decpt = 1; 2736 return nrv_alloc("0", rve, 1); 2737 } 2738 2739 #ifdef SET_INEXACT 2740 try_quick = oldinexact = get_inexact(); 2741 inexact = 1; 2742 #endif 2743 #ifdef Honor_FLT_ROUNDS 2744 if ((rounding = Flt_Rounds) >= 2) { 2745 if (*sign) 2746 rounding = rounding == 2 ? 0 : 2; 2747 else 2748 if (rounding != 2) 2749 rounding = 0; 2750 } 2751 #endif 2752 2753 b = d2b(dval(d), &be, &bbits); 2754 #ifdef Sudden_Underflow 2755 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)); 2756 #else 2757 if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) { 2758 #endif 2759 dval(d2) = dval(d); 2760 word0(d2) &= Frac_mask1; 2761 word0(d2) |= Exp_11; 2762 #ifdef IBM 2763 if (j = 11 - hi0bits(word0(d2) & Frac_mask)) 2764 dval(d2) /= 1 << j; 2765 #endif 2766 2767 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 2768 * log10(x) = log(x) / log(10) 2769 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) 2770 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) 2771 * 2772 * This suggests computing an approximation k to log10(d) by 2773 * 2774 * k = (i - Bias)*0.301029995663981 2775 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); 2776 * 2777 * We want k to be too large rather than too small. 2778 * The error in the first-order Taylor series approximation 2779 * is in our favor, so we just round up the constant enough 2780 * to compensate for any error in the multiplication of 2781 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, 2782 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, 2783 * adding 1e-13 to the constant term more than suffices. 2784 * Hence we adjust the constant term to 0.1760912590558. 2785 * (We could get a more accurate k by invoking log10, 2786 * but this is probably not worthwhile.) 2787 */ 2788 2789 i -= Bias; 2790 #ifdef IBM 2791 i <<= 2; 2792 i += j; 2793 #endif 2794 #ifndef Sudden_Underflow 2795 denorm = 0; 2796 } 2797 else { 2798 /* d is denormalized */ 2799 2800 i = bbits + be + (Bias + (P-1) - 1); 2801 x = i > 32 ? (word0(d) << (64 - i)) | (word1(d) >> (i - 32)) 2802 : word1(d) << (32 - i); 2803 dval(d2) = x; 2804 word0(d2) -= 31*Exp_msk1; /* adjust exponent */ 2805 i -= (Bias + (P-1) - 1) + 1; 2806 denorm = 1; 2807 } 2808 #endif 2809 ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981; 2810 k = (int)ds; 2811 if (ds < 0. && ds != k) 2812 k--; /* want k = floor(ds) */ 2813 k_check = 1; 2814 if (k >= 0 && k <= Ten_pmax) { 2815 if (dval(d) < tens[k]) 2816 k--; 2817 k_check = 0; 2818 } 2819 j = bbits - i - 1; 2820 if (j >= 0) { 2821 b2 = 0; 2822 s2 = j; 2823 } 2824 else { 2825 b2 = -j; 2826 s2 = 0; 2827 } 2828 if (k >= 0) { 2829 b5 = 0; 2830 s5 = k; 2831 s2 += k; 2832 } 2833 else { 2834 b2 -= k; 2835 b5 = -k; 2836 s5 = 0; 2837 } 2838 if (mode < 0 || mode > 9) 2839 mode = 0; 2840 2841 #ifndef SET_INEXACT 2842 #ifdef Check_FLT_ROUNDS 2843 try_quick = Rounding == 1; 2844 #else 2845 try_quick = 1; 2846 #endif 2847 #endif /*SET_INEXACT*/ 2848 2849 if (mode > 5) { 2850 mode -= 4; 2851 try_quick = 0; 2852 } 2853 leftright = 1; 2854 switch(mode) { 2855 case 0: 2856 case 1: 2857 ilim = ilim1 = -1; 2858 i = 18; 2859 ndigits = 0; 2860 break; 2861 case 2: 2862 leftright = 0; 2863 /* no break */ 2864 case 4: 2865 if (ndigits <= 0) 2866 ndigits = 1; 2867 ilim = ilim1 = i = ndigits; 2868 break; 2869 case 3: 2870 leftright = 0; 2871 /* no break */ 2872 case 5: 2873 i = ndigits + k + 1; 2874 ilim = i; 2875 ilim1 = i - 1; 2876 if (i <= 0) 2877 i = 1; 2878 } 2879 s = s0 = rv_alloc(i); 2880 2881 #ifdef Honor_FLT_ROUNDS 2882 if (mode > 1 && rounding != 1) 2883 leftright = 0; 2884 #endif 2885 2886 if (ilim >= 0 && ilim <= Quick_max && try_quick) { 2887 2888 /* Try to get by with floating-point arithmetic. */ 2889 2890 i = 0; 2891 dval(d2) = dval(d); 2892 k0 = k; 2893 ilim0 = ilim; 2894 ieps = 2; /* conservative */ 2895 if (k > 0) { 2896 ds = tens[k&0xf]; 2897 j = k >> 4; 2898 if (j & Bletch) { 2899 /* prevent overflows */ 2900 j &= Bletch - 1; 2901 dval(d) /= bigtens[n_bigtens-1]; 2902 ieps++; 2903 } 2904 for(; j; j >>= 1, i++) 2905 if (j & 1) { 2906 ieps++; 2907 ds *= bigtens[i]; 2908 } 2909 dval(d) /= ds; 2910 } 2911 else if ((j1 = -k)) { 2912 dval(d) *= tens[j1 & 0xf]; 2913 for(j = j1 >> 4; j; j >>= 1, i++) 2914 if (j & 1) { 2915 ieps++; 2916 dval(d) *= bigtens[i]; 2917 } 2918 } 2919 if (k_check && dval(d) < 1. && ilim > 0) { 2920 if (ilim1 <= 0) 2921 goto fast_failed; 2922 ilim = ilim1; 2923 k--; 2924 dval(d) *= 10.; 2925 ieps++; 2926 } 2927 dval(eps) = ieps*dval(d) + 7.; 2928 word0(eps) -= (P-1)*Exp_msk1; 2929 if (ilim == 0) { 2930 S = mhi = 0; 2931 dval(d) -= 5.; 2932 if (dval(d) > dval(eps)) 2933 goto one_digit; 2934 if (dval(d) < -dval(eps)) 2935 goto no_digits; 2936 goto fast_failed; 2937 } 2938 #ifndef No_leftright 2939 if (leftright) { 2940 /* Use Steele & White method of only 2941 * generating digits needed. 2942 */ 2943 dval(eps) = 0.5/tens[ilim-1] - dval(eps); 2944 for(i = 0;;) { 2945 L = dval(d); 2946 dval(d) -= L; 2947 *s++ = '0' + (int)L; 2948 if (dval(d) < dval(eps)) 2949 goto ret1; 2950 if (1. - dval(d) < dval(eps)) 2951 goto bump_up; 2952 if (++i >= ilim) 2953 break; 2954 dval(eps) *= 10.; 2955 dval(d) *= 10.; 2956 } 2957 } 2958 else { 2959 #endif 2960 /* Generate ilim digits, then fix them up. */ 2961 dval(eps) *= tens[ilim-1]; 2962 for(i = 1;; i++, dval(d) *= 10.) { 2963 L = (Long)(dval(d)); 2964 if (!(dval(d) -= L)) 2965 ilim = i; 2966 *s++ = '0' + (int)L; 2967 if (i == ilim) { 2968 if (dval(d) > 0.5 + dval(eps)) 2969 goto bump_up; 2970 else if (dval(d) < 0.5 - dval(eps)) { 2971 while(*--s == '0'); 2972 s++; 2973 goto ret1; 2974 } 2975 break; 2976 } 2977 } 2978 #ifndef No_leftright 2979 } 2980 #endif 2981 fast_failed: 2982 s = s0; 2983 dval(d) = dval(d2); 2984 k = k0; 2985 ilim = ilim0; 2986 } 2987 2988 /* Do we have a "small" integer? */ 2989 2990 if (be >= 0 && k <= Int_max) { 2991 /* Yes. */ 2992 ds = tens[k]; 2993 if (ndigits < 0 && ilim <= 0) { 2994 S = mhi = 0; 2995 if (ilim < 0 || dval(d) < 5*ds || ((dval(d) == 5*ds) && !bias_round_up)) 2996 goto no_digits; 2997 goto one_digit; 2998 } 2999 3000 /* Limit looping by the number of digits to produce. 3001 * Firefox had a crash bug because some plugins reduce 3002 * the precision of double arithmetic. With reduced 3003 * precision "dval(d) -= L*ds" might be imprecise and 3004 * d might not become zero and the loop might not 3005 * terminate. 3006 * 3007 * See https://bugzilla.mozilla.org/show_bug.cgi?id=358569 3008 */ 3009 for(i = 1; i <= k+1; i++, dval(d) *= 10.) { 3010 L = (Long)(dval(d) / ds); 3011 dval(d) -= L*ds; 3012 #ifdef Check_FLT_ROUNDS 3013 /* If FLT_ROUNDS == 2, L will usually be high by 1 */ 3014 if (dval(d) < 0) { 3015 L--; 3016 dval(d) += ds; 3017 } 3018 #endif 3019 *s++ = '0' + (int)L; 3020 if (!dval(d)) { 3021 #ifdef SET_INEXACT 3022 inexact = 0; 3023 #endif 3024 break; 3025 } 3026 if (i == ilim) { 3027 #ifdef Honor_FLT_ROUNDS 3028 if (mode > 1) 3029 switch(rounding) { 3030 case 0: goto ret1; 3031 case 2: goto bump_up; 3032 } 3033 #endif 3034 dval(d) += dval(d); 3035 if (dval(d) > ds || (dval(d) == ds && ((L & 1) || bias_round_up))) { 3036 bump_up: 3037 while(*--s == '9') 3038 if (s == s0) { 3039 k++; 3040 *s = '0'; 3041 break; 3042 } 3043 ++*s++; 3044 } 3045 break; 3046 } 3047 } 3048 goto ret1; 3049 } 3050 3051 m2 = b2; 3052 m5 = b5; 3053 mhi = mlo = 0; 3054 if (leftright) { 3055 i = 3056 #ifndef Sudden_Underflow 3057 denorm ? be + (Bias + (P-1) - 1 + 1) : 3058 #endif 3059 #ifdef IBM 3060 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3); 3061 #else 3062 1 + P - bbits; 3063 #endif 3064 b2 += i; 3065 s2 += i; 3066 mhi = i2b(1); 3067 } 3068 if (m2 > 0 && s2 > 0) { 3069 i = m2 < s2 ? m2 : s2; 3070 b2 -= i; 3071 m2 -= i; 3072 s2 -= i; 3073 } 3074 if (b5 > 0) { 3075 if (leftright) { 3076 if (m5 > 0) { 3077 mhi = pow5mult(mhi, m5); 3078 b1 = mult(mhi, b); 3079 Bfree(b); 3080 b = b1; 3081 } 3082 if ((j = b5 - m5)) 3083 b = pow5mult(b, j); 3084 } 3085 else 3086 b = pow5mult(b, b5); 3087 } 3088 S = i2b(1); 3089 if (s5 > 0) 3090 S = pow5mult(S, s5); 3091 3092 /* Check for special case that d is a normalized power of 2. */ 3093 3094 spec_case = 0; 3095 if ((mode < 2 || leftright) 3096 #ifdef Honor_FLT_ROUNDS 3097 && rounding == 1 3098 #endif 3099 ) { 3100 if (!word1(d) && !(word0(d) & Bndry_mask) 3101 #ifndef Sudden_Underflow 3102 && word0(d) & (Exp_mask & ~Exp_msk1) 3103 #endif 3104 ) { 3105 /* The special case */ 3106 b2 += Log2P; 3107 s2 += Log2P; 3108 spec_case = 1; 3109 } 3110 } 3111 3112 /* Arrange for convenient computation of quotients: 3113 * shift left if necessary so divisor has 4 leading 0 bits. 3114 * 3115 * Perhaps we should just compute leading 28 bits of S once 3116 * and for all and pass them and a shift to quorem, so it 3117 * can do shifts and ors to compute the numerator for q. 3118 */ 3119 #ifdef Pack_32 3120 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)) 3121 i = 32 - i; 3122 #else 3123 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)) 3124 i = 16 - i; 3125 #endif 3126 if (i > 4) { 3127 i -= 4; 3128 b2 += i; 3129 m2 += i; 3130 s2 += i; 3131 } 3132 else if (i < 4) { 3133 i += 28; 3134 b2 += i; 3135 m2 += i; 3136 s2 += i; 3137 } 3138 if (b2 > 0) 3139 b = lshift(b, b2); 3140 if (s2 > 0) 3141 S = lshift(S, s2); 3142 if (k_check) { 3143 if (cmp(b,S) < 0) { 3144 k--; 3145 b = multadd(b, 10, 0); /* we botched the k estimate */ 3146 if (leftright) 3147 mhi = multadd(mhi, 10, 0); 3148 ilim = ilim1; 3149 } 3150 } 3151 if (ilim <= 0 && (mode == 3 || mode == 5)) { 3152 S = multadd(S, 5, 0); 3153 if (ilim < 0 || cmp(b, S) < 0 || ((cmp(b, S) == 0) && !bias_round_up)) { 3154 /* no digits, fcvt style */ 3155 no_digits: 3156 k = -1 - ndigits; 3157 goto ret; 3158 } 3159 one_digit: 3160 *s++ = '1'; 3161 k++; 3162 goto ret; 3163 } 3164 if (leftright) { 3165 if (m2 > 0) 3166 mhi = lshift(mhi, m2); 3167 3168 /* Compute mlo -- check for special case 3169 * that d is a normalized power of 2. 3170 */ 3171 3172 mlo = mhi; 3173 if (spec_case) { 3174 mhi = Balloc(mhi->k); 3175 Bcopy(mhi, mlo); 3176 mhi = lshift(mhi, Log2P); 3177 } 3178 3179 for(i = 1;;i++) { 3180 dig = quorem(b,S) + '0'; 3181 /* Do we yet have the shortest decimal string 3182 * that will round to d? 3183 */ 3184 j = cmp(b, mlo); 3185 delta = diff(S, mhi); 3186 j1 = delta->sign ? 1 : cmp(b, delta); 3187 Bfree(delta); 3188 #ifndef ROUND_BIASED 3189 if (j1 == 0 && mode != 1 && !(word1(d) & 1) 3190 #ifdef Honor_FLT_ROUNDS 3191 && rounding >= 1 3192 #endif 3193 ) { 3194 if (dig == '9') 3195 goto round_9_up; 3196 if (j > 0) 3197 dig++; 3198 #ifdef SET_INEXACT 3199 else if (!b->x[0] && b->wds <= 1) 3200 inexact = 0; 3201 #endif 3202 *s++ = dig; 3203 goto ret; 3204 } 3205 #endif 3206 if (j < 0 || (j == 0 && mode != 1 3207 #ifndef ROUND_BIASED 3208 && !(word1(d) & 1) 3209 #endif 3210 )) { 3211 if (!b->x[0] && b->wds <= 1) { 3212 #ifdef SET_INEXACT 3213 inexact = 0; 3214 #endif 3215 goto accept_dig; 3216 } 3217 #ifdef Honor_FLT_ROUNDS 3218 if (mode > 1) 3219 switch(rounding) { 3220 case 0: goto accept_dig; 3221 case 2: goto keep_dig; 3222 } 3223 #endif /*Honor_FLT_ROUNDS*/ 3224 if (j1 > 0) { 3225 b = lshift(b, 1); 3226 j1 = cmp(b, S); 3227 if ((j1 > 0 || (j1 == 0 && ((dig & 1) || bias_round_up))) 3228 && dig++ == '9') 3229 goto round_9_up; 3230 } 3231 accept_dig: 3232 *s++ = dig; 3233 goto ret; 3234 } 3235 if (j1 > 0) { 3236 #ifdef Honor_FLT_ROUNDS 3237 if (!rounding) 3238 goto accept_dig; 3239 #endif 3240 if (dig == '9') { /* possible if i == 1 */ 3241 round_9_up: 3242 *s++ = '9'; 3243 goto roundoff; 3244 } 3245 *s++ = dig + 1; 3246 goto ret; 3247 } 3248 #ifdef Honor_FLT_ROUNDS 3249 keep_dig: 3250 #endif 3251 *s++ = dig; 3252 if (i == ilim) 3253 break; 3254 b = multadd(b, 10, 0); 3255 if (mlo == mhi) 3256 mlo = mhi = multadd(mhi, 10, 0); 3257 else { 3258 mlo = multadd(mlo, 10, 0); 3259 mhi = multadd(mhi, 10, 0); 3260 } 3261 } 3262 } 3263 else 3264 for(i = 1;; i++) { 3265 *s++ = dig = quorem(b,S) + '0'; 3266 if (!b->x[0] && b->wds <= 1) { 3267 #ifdef SET_INEXACT 3268 inexact = 0; 3269 #endif 3270 goto ret; 3271 } 3272 if (i >= ilim) 3273 break; 3274 b = multadd(b, 10, 0); 3275 } 3276 3277 /* Round off last digit */ 3278 3279 #ifdef Honor_FLT_ROUNDS 3280 switch(rounding) { 3281 case 0: goto trimzeros; 3282 case 2: goto roundoff; 3283 } 3284 #endif 3285 b = lshift(b, 1); 3286 j = cmp(b, S); 3287 if (j > 0 || (j == 0 && ((dig & 1) || bias_round_up))) { 3288 roundoff: 3289 while(*--s == '9') 3290 if (s == s0) { 3291 k++; 3292 *s++ = '1'; 3293 goto ret; 3294 } 3295 ++*s++; 3296 } 3297 else { 3298 /* trimzeros: (never used) */ 3299 while(*--s == '0'); 3300 s++; 3301 } 3302 ret: 3303 Bfree(S); 3304 if (mhi) { 3305 if (mlo && mlo != mhi) 3306 Bfree(mlo); 3307 Bfree(mhi); 3308 } 3309 ret1: 3310 #ifdef SET_INEXACT 3311 if (inexact) { 3312 if (!oldinexact) { 3313 word0(d) = Exp_1 + (70 << Exp_shift); 3314 word1(d) = 0; 3315 dval(d) += 1.; 3316 } 3317 } 3318 else if (!oldinexact) 3319 clear_inexact(); 3320 #endif 3321 Bfree(b); 3322 *s = 0; 3323 *decpt = k + 1; 3324 if (rve) 3325 *rve = s; 3326 return s0; 3327 } 3328 #ifdef __cplusplus 3329 } 3330 #endif 3331
