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. Unless Trust_FLT_ROUNDS 76 * is also #defined, fegetround() will be queried for the rounding mode. 77 * Note that both FLT_ROUNDS and fegetround() are specified by the C99 78 * standard (and are specified to be consistent, with fesetround() 79 * affecting the value of FLT_ROUNDS), but that some (Linux) systems 80 * do not work correctly in this regard, so using fegetround() is more 81 * portable than using FLT_FOUNDS directly. 82 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 83 * and Honor_FLT_ROUNDS is not #defined. 84 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines 85 * that use extended-precision instructions to compute rounded 86 * products and quotients) with IBM. 87 * #define ROUND_BIASED for IEEE-format with biased rounding. 88 * #define Inaccurate_Divide for IEEE-format with correctly rounded 89 * products but inaccurate quotients, e.g., for Intel i860. 90 * #define NO_LONG_LONG on machines that do not have a "long long" 91 * integer type (of >= 64 bits). On such machines, you can 92 * #define Just_16 to store 16 bits per 32-bit Long when doing 93 * high-precision integer arithmetic. Whether this speeds things 94 * up or slows things down depends on the machine and the number 95 * being converted. If long long is available and the name is 96 * something other than "long long", #define Llong to be the name, 97 * and if "unsigned Llong" does not work as an unsigned version of 98 * Llong, #define #ULLong to be the corresponding unsigned type. 99 * #define KR_headers for old-style C function headers. 100 * #define Bad_float_h if your system lacks a float.h or if it does not 101 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP, 102 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX. 103 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n) 104 * if memory is available and otherwise does something you deem 105 * appropriate. If MALLOC is undefined, malloc will be invoked 106 * directly -- and assumed always to succeed. Similarly, if you 107 * want something other than the system's free() to be called to 108 * recycle memory acquired from MALLOC, #define FREE to be the 109 * name of the alternate routine. (FREE or free is only called in 110 * pathological cases, e.g., in a dtoa call after a dtoa return in 111 * mode 3 with thousands of digits requested.) 112 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making 113 * memory allocations from a private pool of memory when possible. 114 * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes, 115 * unless #defined to be a different length. This default length 116 * suffices to get rid of MALLOC calls except for unusual cases, 117 * such as decimal-to-binary conversion of a very long string of 118 * digits. The longest string dtoa can return is about 751 bytes 119 * long. For conversions by strtod of strings of 800 digits and 120 * all dtoa conversions in single-threaded executions with 8-byte 121 * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte 122 * pointers, PRIVATE_MEM >= 7112 appears adequate. 123 * #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK 124 * #defined automatically on IEEE systems. On such systems, 125 * when INFNAN_CHECK is #defined, strtod checks 126 * for Infinity and NaN (case insensitively). On some systems 127 * (e.g., some HP systems), it may be necessary to #define NAN_WORD0 128 * appropriately -- to the most significant word of a quiet NaN. 129 * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.) 130 * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined, 131 * strtod also accepts (case insensitively) strings of the form 132 * NaN(x), where x is a string of hexadecimal digits and spaces; 133 * if there is only one string of hexadecimal digits, it is taken 134 * for the 52 fraction bits of the resulting NaN; if there are two 135 * or more strings of hex digits, the first is for the high 20 bits, 136 * the second and subsequent for the low 32 bits, with intervening 137 * white space ignored; but if this results in none of the 52 138 * fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0 139 * and NAN_WORD1 are used instead. 140 * #define MULTIPLE_THREADS if the system offers preemptively scheduled 141 * multiple threads. In this case, you must provide (or suitably 142 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed 143 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed 144 * in pow5mult, ensures lazy evaluation of only one copy of high 145 * powers of 5; omitting this lock would introduce a small 146 * probability of wasting memory, but would otherwise be harmless.) 147 * You must also invoke freedtoa(s) to free the value s returned by 148 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined. 149 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that 150 * avoids underflows on inputs whose result does not underflow. 151 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format 152 * floating-point numbers and flushes underflows to zero rather 153 * than implementing gradual underflow, then you must also #define 154 * Sudden_Underflow. 155 * #define USE_LOCALE to use the current locale's decimal_point value. 156 * #define SET_INEXACT if IEEE arithmetic is being used and extra 157 * computation should be done to set the inexact flag when the 158 * result is inexact and avoid setting inexact when the result 159 * is exact. In this case, dtoa.c must be compiled in 160 * an environment, perhaps provided by #include "dtoa.c" in a 161 * suitable wrapper, that defines two functions, 162 * int get_inexact(void); 163 * void clear_inexact(void); 164 * such that get_inexact() returns a nonzero value if the 165 * inexact bit is already set, and clear_inexact() sets the 166 * inexact bit to 0. When SET_INEXACT is #defined, strtod 167 * also does extra computations to set the underflow and overflow 168 * flags when appropriate (i.e., when the result is tiny and 169 * inexact or when it is a numeric value rounded to +-infinity). 170 * #define NO_ERRNO if strtod should not assign errno = ERANGE when 171 * the result overflows to +-Infinity or underflows to 0. 172 * #define NO_HEX_FP to omit recognition of hexadecimal floating-point 173 * values by strtod. 174 * #define NO_STRTOD_BIGCOMP (on IEEE-arithmetic systems only for now) 175 * to disable logic for "fast" testing of very long input strings 176 * to strtod. This testing proceeds by initially truncating the 177 * input string, then if necessary comparing the whole string with 178 * a decimal expansion to decide close cases. This logic is only 179 * used for input more than STRTOD_DIGLIM digits long (default 40). 180 */ 181 182 #define IEEE_8087 183 #define NO_HEX_FP 184 185 #ifndef Long 186 #if __LP64__ 187 #define Long int 188 #else 189 #define Long long 190 #endif 191 #endif 192 #ifndef ULong 193 typedef unsigned Long ULong; 194 #endif 195 196 #ifdef DEBUG 197 #include "stdio.h" 198 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} 199 #endif 200 201 #include "stdlib.h" 202 #include "string.h" 203 204 #ifdef USE_LOCALE 205 #include "locale.h" 206 #endif 207 208 #ifdef Honor_FLT_ROUNDS 209 #ifndef Trust_FLT_ROUNDS 210 #include <fenv.h> 211 #endif 212 #endif 213 214 #ifdef MALLOC 215 #ifdef KR_headers 216 extern char *MALLOC(); 217 #else 218 extern void *MALLOC(size_t); 219 #endif 220 #else 221 #define MALLOC malloc 222 #endif 223 224 #ifndef Omit_Private_Memory 225 #ifndef PRIVATE_MEM 226 #define PRIVATE_MEM 2304 227 #endif 228 #define PRIVATE_mem ((unsigned)((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))) 229 static double private_mem[PRIVATE_mem], *pmem_next = private_mem; 230 #endif 231 232 #undef IEEE_Arith 233 #undef Avoid_Underflow 234 #ifdef IEEE_MC68k 235 #define IEEE_Arith 236 #endif 237 #ifdef IEEE_8087 238 #define IEEE_Arith 239 #endif 240 241 #ifdef IEEE_Arith 242 #ifndef NO_INFNAN_CHECK 243 #undef INFNAN_CHECK 244 #define INFNAN_CHECK 245 #endif 246 #else 247 #undef INFNAN_CHECK 248 #define NO_STRTOD_BIGCOMP 249 #endif 250 251 #include "errno.h" 252 253 #ifdef Bad_float_h 254 255 #ifdef IEEE_Arith 256 #define DBL_DIG 15 257 #define DBL_MAX_10_EXP 308 258 #define DBL_MAX_EXP 1024 259 #define FLT_RADIX 2 260 #endif /*IEEE_Arith*/ 261 262 #ifdef IBM 263 #define DBL_DIG 16 264 #define DBL_MAX_10_EXP 75 265 #define DBL_MAX_EXP 63 266 #define FLT_RADIX 16 267 #define DBL_MAX 7.2370055773322621e+75 268 #endif 269 270 #ifdef VAX 271 #define DBL_DIG 16 272 #define DBL_MAX_10_EXP 38 273 #define DBL_MAX_EXP 127 274 #define FLT_RADIX 2 275 #define DBL_MAX 1.7014118346046923e+38 276 #endif 277 278 #ifndef LONG_MAX 279 #define LONG_MAX 2147483647 280 #endif 281 282 #else /* ifndef Bad_float_h */ 283 #include "float.h" 284 #endif /* Bad_float_h */ 285 286 #ifndef __MATH_H__ 287 #include "math.h" 288 #endif 289 290 namespace dmg_fp { 291 292 #ifndef CONST 293 #ifdef KR_headers 294 #define CONST /* blank */ 295 #else 296 #define CONST const 297 #endif 298 #endif 299 300 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1 301 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined. 302 #endif 303 304 typedef union { double d; ULong L[2]; } U; 305 306 #ifdef IEEE_8087 307 #define word0(x) (x)->L[1] 308 #define word1(x) (x)->L[0] 309 #else 310 #define word0(x) (x)->L[0] 311 #define word1(x) (x)->L[1] 312 #endif 313 #define dval(x) (x)->d 314 315 #ifndef STRTOD_DIGLIM 316 #define STRTOD_DIGLIM 40 317 #endif 318 319 #ifdef DIGLIM_DEBUG 320 extern int strtod_diglim; 321 #else 322 #define strtod_diglim STRTOD_DIGLIM 323 #endif 324 325 /* The following definition of Storeinc is appropriate for MIPS processors. 326 * An alternative that might be better on some machines is 327 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) 328 */ 329 #if defined(IEEE_8087) + defined(VAX) 330 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \ 331 ((unsigned short *)a)[0] = (unsigned short)c, a++) 332 #else 333 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \ 334 ((unsigned short *)a)[1] = (unsigned short)c, a++) 335 #endif 336 337 /* #define P DBL_MANT_DIG */ 338 /* Ten_pmax = floor(P*log(2)/log(5)) */ 339 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ 340 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ 341 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ 342 343 #ifdef IEEE_Arith 344 #define Exp_shift 20 345 #define Exp_shift1 20 346 #define Exp_msk1 0x100000 347 #define Exp_msk11 0x100000 348 #define Exp_mask 0x7ff00000 349 #define P 53 350 #define Nbits 53 351 #define Bias 1023 352 #define Emax 1023 353 #define Emin (-1022) 354 #define Exp_1 0x3ff00000 355 #define Exp_11 0x3ff00000 356 #define Ebits 11 357 #define Frac_mask 0xfffff 358 #define Frac_mask1 0xfffff 359 #define Ten_pmax 22 360 #define Bletch 0x10 361 #define Bndry_mask 0xfffff 362 #define Bndry_mask1 0xfffff 363 #define LSB 1 364 #define Sign_bit 0x80000000 365 #define Log2P 1 366 #define Tiny0 0 367 #define Tiny1 1 368 #define Quick_max 14 369 #define Int_max 14 370 #ifndef NO_IEEE_Scale 371 #define Avoid_Underflow 372 #ifdef Flush_Denorm /* debugging option */ 373 #undef Sudden_Underflow 374 #endif 375 #endif 376 377 #ifndef Flt_Rounds 378 #ifdef FLT_ROUNDS 379 #define Flt_Rounds FLT_ROUNDS 380 #else 381 #define Flt_Rounds 1 382 #endif 383 #endif /*Flt_Rounds*/ 384 385 #ifdef Honor_FLT_ROUNDS 386 #undef Check_FLT_ROUNDS 387 #define Check_FLT_ROUNDS 388 #else 389 #define Rounding Flt_Rounds 390 #endif 391 392 #else /* ifndef IEEE_Arith */ 393 #undef Check_FLT_ROUNDS 394 #undef Honor_FLT_ROUNDS 395 #undef SET_INEXACT 396 #undef Sudden_Underflow 397 #define Sudden_Underflow 398 #ifdef IBM 399 #undef Flt_Rounds 400 #define Flt_Rounds 0 401 #define Exp_shift 24 402 #define Exp_shift1 24 403 #define Exp_msk1 0x1000000 404 #define Exp_msk11 0x1000000 405 #define Exp_mask 0x7f000000 406 #define P 14 407 #define Nbits 56 408 #define Bias 65 409 #define Emax 248 410 #define Emin (-260) 411 #define Exp_1 0x41000000 412 #define Exp_11 0x41000000 413 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */ 414 #define Frac_mask 0xffffff 415 #define Frac_mask1 0xffffff 416 #define Bletch 4 417 #define Ten_pmax 22 418 #define Bndry_mask 0xefffff 419 #define Bndry_mask1 0xffffff 420 #define LSB 1 421 #define Sign_bit 0x80000000 422 #define Log2P 4 423 #define Tiny0 0x100000 424 #define Tiny1 0 425 #define Quick_max 14 426 #define Int_max 15 427 #else /* VAX */ 428 #undef Flt_Rounds 429 #define Flt_Rounds 1 430 #define Exp_shift 23 431 #define Exp_shift1 7 432 #define Exp_msk1 0x80 433 #define Exp_msk11 0x800000 434 #define Exp_mask 0x7f80 435 #define P 56 436 #define Nbits 56 437 #define Bias 129 438 #define Emax 126 439 #define Emin (-129) 440 #define Exp_1 0x40800000 441 #define Exp_11 0x4080 442 #define Ebits 8 443 #define Frac_mask 0x7fffff 444 #define Frac_mask1 0xffff007f 445 #define Ten_pmax 24 446 #define Bletch 2 447 #define Bndry_mask 0xffff007f 448 #define Bndry_mask1 0xffff007f 449 #define LSB 0x10000 450 #define Sign_bit 0x8000 451 #define Log2P 1 452 #define Tiny0 0x80 453 #define Tiny1 0 454 #define Quick_max 15 455 #define Int_max 15 456 #endif /* IBM, VAX */ 457 #endif /* IEEE_Arith */ 458 459 #ifndef IEEE_Arith 460 #define ROUND_BIASED 461 #endif 462 463 #ifdef RND_PRODQUOT 464 #define rounded_product(a,b) a = rnd_prod(a, b) 465 #define rounded_quotient(a,b) a = rnd_quot(a, b) 466 #ifdef KR_headers 467 extern double rnd_prod(), rnd_quot(); 468 #else 469 extern double rnd_prod(double, double), rnd_quot(double, double); 470 #endif 471 #else 472 #define rounded_product(a,b) a *= b 473 #define rounded_quotient(a,b) a /= b 474 #endif 475 476 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) 477 #define Big1 0xffffffff 478 479 #ifndef Pack_32 480 #define Pack_32 481 #endif 482 483 typedef struct BCinfo BCinfo; 484 struct 485 BCinfo { int dp0, dp1, dplen, dsign, e0, inexact, nd, nd0, rounding, scale, uflchk; }; 486 487 #ifdef KR_headers 488 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff) 489 #else 490 #define FFFFFFFF 0xffffffffUL 491 #endif 492 493 #ifdef NO_LONG_LONG 494 #undef ULLong 495 #ifdef Just_16 496 #undef Pack_32 497 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long. 498 * This makes some inner loops simpler and sometimes saves work 499 * during multiplications, but it often seems to make things slightly 500 * slower. Hence the default is now to store 32 bits per Long. 501 */ 502 #endif 503 #else /* long long available */ 504 #ifndef Llong 505 #define Llong long long 506 #endif 507 #ifndef ULLong 508 #define ULLong unsigned Llong 509 #endif 510 #endif /* NO_LONG_LONG */ 511 512 #ifndef MULTIPLE_THREADS 513 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/ 514 #define FREE_DTOA_LOCK(n) /*nothing*/ 515 #endif 516 517 #define Kmax 7 518 519 double strtod(const char *s00, char **se); 520 char *dtoa(double d, int mode, int ndigits, 521 int *decpt, int *sign, char **rve); 522 523 struct 524 Bigint { 525 struct Bigint *next; 526 int k, maxwds, sign, wds; 527 ULong x[1]; 528 }; 529 530 typedef struct Bigint Bigint; 531 532 static Bigint *freelist[Kmax+1]; 533 534 static Bigint * 535 Balloc 536 #ifdef KR_headers 537 (k) int k; 538 #else 539 (int k) 540 #endif 541 { 542 int x; 543 Bigint *rv; 544 #ifndef Omit_Private_Memory 545 unsigned int len; 546 #endif 547 548 ACQUIRE_DTOA_LOCK(0); 549 /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */ 550 /* but this case seems very unlikely. */ 551 if (k <= Kmax && freelist[k]) { 552 rv = freelist[k]; 553 freelist[k] = rv->next; 554 } 555 else { 556 x = 1 << k; 557 #ifdef Omit_Private_Memory 558 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong)); 559 #else 560 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1) 561 /sizeof(double); 562 if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) { 563 rv = (Bigint*)pmem_next; 564 pmem_next += len; 565 } 566 else 567 rv = (Bigint*)MALLOC(len*sizeof(double)); 568 #endif 569 rv->k = k; 570 rv->maxwds = x; 571 } 572 FREE_DTOA_LOCK(0); 573 rv->sign = rv->wds = 0; 574 return rv; 575 } 576 577 static void 578 Bfree 579 #ifdef KR_headers 580 (v) Bigint *v; 581 #else 582 (Bigint *v) 583 #endif 584 { 585 if (v) { 586 if (v->k > Kmax) 587 #ifdef FREE 588 FREE((void*)v); 589 #else 590 free((void*)v); 591 #endif 592 else { 593 ACQUIRE_DTOA_LOCK(0); 594 v->next = freelist[v->k]; 595 freelist[v->k] = v; 596 FREE_DTOA_LOCK(0); 597 } 598 } 599 } 600 601 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \ 602 y->wds*sizeof(Long) + 2*sizeof(int)) 603 604 static Bigint * 605 multadd 606 #ifdef KR_headers 607 (b, m, a) Bigint *b; int m, a; 608 #else 609 (Bigint *b, int m, int a) /* multiply by m and add a */ 610 #endif 611 { 612 int i, wds; 613 #ifdef ULLong 614 ULong *x; 615 ULLong carry, y; 616 #else 617 ULong carry, *x, y; 618 #ifdef Pack_32 619 ULong xi, z; 620 #endif 621 #endif 622 Bigint *b1; 623 624 wds = b->wds; 625 x = b->x; 626 i = 0; 627 carry = a; 628 do { 629 #ifdef ULLong 630 y = *x * (ULLong)m + carry; 631 carry = y >> 32; 632 *x++ = y & FFFFFFFF; 633 #else 634 #ifdef Pack_32 635 xi = *x; 636 y = (xi & 0xffff) * m + carry; 637 z = (xi >> 16) * m + (y >> 16); 638 carry = z >> 16; 639 *x++ = (z << 16) + (y & 0xffff); 640 #else 641 y = *x * m + carry; 642 carry = y >> 16; 643 *x++ = y & 0xffff; 644 #endif 645 #endif 646 } 647 while(++i < wds); 648 if (carry) { 649 if (wds >= b->maxwds) { 650 b1 = Balloc(b->k+1); 651 Bcopy(b1, b); 652 Bfree(b); 653 b = b1; 654 } 655 b->x[wds++] = (ULong)carry; 656 b->wds = wds; 657 } 658 return b; 659 } 660 661 static Bigint * 662 s2b 663 #ifdef KR_headers 664 (s, nd0, nd, y9, dplen) CONST char *s; int nd0, nd, dplen; ULong y9; 665 #else 666 (CONST char *s, int nd0, int nd, ULong y9, int dplen) 667 #endif 668 { 669 Bigint *b; 670 int i, k; 671 Long x, y; 672 673 x = (nd + 8) / 9; 674 for(k = 0, y = 1; x > y; y <<= 1, k++) ; 675 #ifdef Pack_32 676 b = Balloc(k); 677 b->x[0] = y9; 678 b->wds = 1; 679 #else 680 b = Balloc(k+1); 681 b->x[0] = y9 & 0xffff; 682 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1; 683 #endif 684 685 i = 9; 686 if (9 < nd0) { 687 s += 9; 688 do b = multadd(b, 10, *s++ - '0'); 689 while(++i < nd0); 690 s += dplen; 691 } 692 else 693 s += dplen + 9; 694 for(; i < nd; i++) 695 b = multadd(b, 10, *s++ - '0'); 696 return b; 697 } 698 699 static int 700 hi0bits 701 #ifdef KR_headers 702 (x) ULong x; 703 #else 704 (ULong x) 705 #endif 706 { 707 int k = 0; 708 709 if (!(x & 0xffff0000)) { 710 k = 16; 711 x <<= 16; 712 } 713 if (!(x & 0xff000000)) { 714 k += 8; 715 x <<= 8; 716 } 717 if (!(x & 0xf0000000)) { 718 k += 4; 719 x <<= 4; 720 } 721 if (!(x & 0xc0000000)) { 722 k += 2; 723 x <<= 2; 724 } 725 if (!(x & 0x80000000)) { 726 k++; 727 if (!(x & 0x40000000)) 728 return 32; 729 } 730 return k; 731 } 732 733 static int 734 lo0bits 735 #ifdef KR_headers 736 (y) ULong *y; 737 #else 738 (ULong *y) 739 #endif 740 { 741 int k; 742 ULong x = *y; 743 744 if (x & 7) { 745 if (x & 1) 746 return 0; 747 if (x & 2) { 748 *y = x >> 1; 749 return 1; 750 } 751 *y = x >> 2; 752 return 2; 753 } 754 k = 0; 755 if (!(x & 0xffff)) { 756 k = 16; 757 x >>= 16; 758 } 759 if (!(x & 0xff)) { 760 k += 8; 761 x >>= 8; 762 } 763 if (!(x & 0xf)) { 764 k += 4; 765 x >>= 4; 766 } 767 if (!(x & 0x3)) { 768 k += 2; 769 x >>= 2; 770 } 771 if (!(x & 1)) { 772 k++; 773 x >>= 1; 774 if (!x) 775 return 32; 776 } 777 *y = x; 778 return k; 779 } 780 781 static Bigint * 782 i2b 783 #ifdef KR_headers 784 (i) int i; 785 #else 786 (int i) 787 #endif 788 { 789 Bigint *b; 790 791 b = Balloc(1); 792 b->x[0] = i; 793 b->wds = 1; 794 return b; 795 } 796 797 static Bigint * 798 mult 799 #ifdef KR_headers 800 (a, b) Bigint *a, *b; 801 #else 802 (Bigint *a, Bigint *b) 803 #endif 804 { 805 Bigint *c; 806 int k, wa, wb, wc; 807 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; 808 ULong y; 809 #ifdef ULLong 810 ULLong carry, z; 811 #else 812 ULong carry, z; 813 #ifdef Pack_32 814 ULong z2; 815 #endif 816 #endif 817 818 if (a->wds < b->wds) { 819 c = a; 820 a = b; 821 b = c; 822 } 823 k = a->k; 824 wa = a->wds; 825 wb = b->wds; 826 wc = wa + wb; 827 if (wc > a->maxwds) 828 k++; 829 c = Balloc(k); 830 for(x = c->x, xa = x + wc; x < xa; x++) 831 *x = 0; 832 xa = a->x; 833 xae = xa + wa; 834 xb = b->x; 835 xbe = xb + wb; 836 xc0 = c->x; 837 #ifdef ULLong 838 for(; xb < xbe; xc0++) { 839 y = *xb++; 840 if (y) { 841 x = xa; 842 xc = xc0; 843 carry = 0; 844 do { 845 z = *x++ * (ULLong)y + *xc + carry; 846 carry = z >> 32; 847 *xc++ = z & FFFFFFFF; 848 } 849 while(x < xae); 850 *xc = (ULong)carry; 851 } 852 } 853 #else 854 #ifdef Pack_32 855 for(; xb < xbe; xb++, xc0++) { 856 if (y = *xb & 0xffff) { 857 x = xa; 858 xc = xc0; 859 carry = 0; 860 do { 861 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 862 carry = z >> 16; 863 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 864 carry = z2 >> 16; 865 Storeinc(xc, z2, z); 866 } 867 while(x < xae); 868 *xc = carry; 869 } 870 if (y = *xb >> 16) { 871 x = xa; 872 xc = xc0; 873 carry = 0; 874 z2 = *xc; 875 do { 876 z = (*x & 0xffff) * y + (*xc >> 16) + carry; 877 carry = z >> 16; 878 Storeinc(xc, z, z2); 879 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 880 carry = z2 >> 16; 881 } 882 while(x < xae); 883 *xc = z2; 884 } 885 } 886 #else 887 for(; xb < xbe; xc0++) { 888 if (y = *xb++) { 889 x = xa; 890 xc = xc0; 891 carry = 0; 892 do { 893 z = *x++ * y + *xc + carry; 894 carry = z >> 16; 895 *xc++ = z & 0xffff; 896 } 897 while(x < xae); 898 *xc = carry; 899 } 900 } 901 #endif 902 #endif 903 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ; 904 c->wds = wc; 905 return c; 906 } 907 908 static Bigint *p5s; 909 910 static Bigint * 911 pow5mult 912 #ifdef KR_headers 913 (b, k) Bigint *b; int k; 914 #else 915 (Bigint *b, int k) 916 #endif 917 { 918 Bigint *b1, *p5, *p51; 919 int i; 920 static int p05[3] = { 5, 25, 125 }; 921 922 i = k & 3; 923 if (i) 924 b = multadd(b, p05[i-1], 0); 925 926 if (!(k >>= 2)) 927 return b; 928 p5 = p5s; 929 if (!p5) { 930 /* first time */ 931 #ifdef MULTIPLE_THREADS 932 ACQUIRE_DTOA_LOCK(1); 933 p5 = p5s; 934 if (!p5) { 935 p5 = p5s = i2b(625); 936 p5->next = 0; 937 } 938 FREE_DTOA_LOCK(1); 939 #else 940 p5 = p5s = i2b(625); 941 p5->next = 0; 942 #endif 943 } 944 for(;;) { 945 if (k & 1) { 946 b1 = mult(b, p5); 947 Bfree(b); 948 b = b1; 949 } 950 if (!(k >>= 1)) 951 break; 952 p51 = p5->next; 953 if (!p51) { 954 #ifdef MULTIPLE_THREADS 955 ACQUIRE_DTOA_LOCK(1); 956 p51 = p5->next; 957 if (!p51) { 958 p51 = p5->next = mult(p5,p5); 959 p51->next = 0; 960 } 961 FREE_DTOA_LOCK(1); 962 #else 963 p51 = p5->next = mult(p5,p5); 964 p51->next = 0; 965 #endif 966 } 967 p5 = p51; 968 } 969 return b; 970 } 971 972 static Bigint * 973 lshift 974 #ifdef KR_headers 975 (b, k) Bigint *b; int k; 976 #else 977 (Bigint *b, int k) 978 #endif 979 { 980 int i, k1, n, n1; 981 Bigint *b1; 982 ULong *x, *x1, *xe, z; 983 984 #ifdef Pack_32 985 n = k >> 5; 986 #else 987 n = k >> 4; 988 #endif 989 k1 = b->k; 990 n1 = n + b->wds + 1; 991 for(i = b->maxwds; n1 > i; i <<= 1) 992 k1++; 993 b1 = Balloc(k1); 994 x1 = b1->x; 995 for(i = 0; i < n; i++) 996 *x1++ = 0; 997 x = b->x; 998 xe = x + b->wds; 999 #ifdef Pack_32 1000 if (k &= 0x1f) { 1001 k1 = 32 - k; 1002 z = 0; 1003 do { 1004 *x1++ = *x << k | z; 1005 z = *x++ >> k1; 1006 } 1007 while(x < xe); 1008 *x1 = z; 1009 if (*x1) 1010 ++n1; 1011 } 1012 #else 1013 if (k &= 0xf) { 1014 k1 = 16 - k; 1015 z = 0; 1016 do { 1017 *x1++ = *x << k & 0xffff | z; 1018 z = *x++ >> k1; 1019 } 1020 while(x < xe); 1021 if (*x1 = z) 1022 ++n1; 1023 } 1024 #endif 1025 else do 1026 *x1++ = *x++; 1027 while(x < xe); 1028 b1->wds = n1 - 1; 1029 Bfree(b); 1030 return b1; 1031 } 1032 1033 static int 1034 cmp 1035 #ifdef KR_headers 1036 (a, b) Bigint *a, *b; 1037 #else 1038 (Bigint *a, Bigint *b) 1039 #endif 1040 { 1041 ULong *xa, *xa0, *xb, *xb0; 1042 int i, j; 1043 1044 i = a->wds; 1045 j = b->wds; 1046 #ifdef DEBUG 1047 if (i > 1 && !a->x[i-1]) 1048 Bug("cmp called with a->x[a->wds-1] == 0"); 1049 if (j > 1 && !b->x[j-1]) 1050 Bug("cmp called with b->x[b->wds-1] == 0"); 1051 #endif 1052 if (i -= j) 1053 return i; 1054 xa0 = a->x; 1055 xa = xa0 + j; 1056 xb0 = b->x; 1057 xb = xb0 + j; 1058 for(;;) { 1059 if (*--xa != *--xb) 1060 return *xa < *xb ? -1 : 1; 1061 if (xa <= xa0) 1062 break; 1063 } 1064 return 0; 1065 } 1066 1067 static Bigint * 1068 diff 1069 #ifdef KR_headers 1070 (a, b) Bigint *a, *b; 1071 #else 1072 (Bigint *a, Bigint *b) 1073 #endif 1074 { 1075 Bigint *c; 1076 int i, wa, wb; 1077 ULong *xa, *xae, *xb, *xbe, *xc; 1078 #ifdef ULLong 1079 ULLong borrow, y; 1080 #else 1081 ULong borrow, y; 1082 #ifdef Pack_32 1083 ULong z; 1084 #endif 1085 #endif 1086 1087 i = cmp(a,b); 1088 if (!i) { 1089 c = Balloc(0); 1090 c->wds = 1; 1091 c->x[0] = 0; 1092 return c; 1093 } 1094 if (i < 0) { 1095 c = a; 1096 a = b; 1097 b = c; 1098 i = 1; 1099 } 1100 else 1101 i = 0; 1102 c = Balloc(a->k); 1103 c->sign = i; 1104 wa = a->wds; 1105 xa = a->x; 1106 xae = xa + wa; 1107 wb = b->wds; 1108 xb = b->x; 1109 xbe = xb + wb; 1110 xc = c->x; 1111 borrow = 0; 1112 #ifdef ULLong 1113 do { 1114 y = (ULLong)*xa++ - *xb++ - borrow; 1115 borrow = y >> 32 & (ULong)1; 1116 *xc++ = y & FFFFFFFF; 1117 } 1118 while(xb < xbe); 1119 while(xa < xae) { 1120 y = *xa++ - borrow; 1121 borrow = y >> 32 & (ULong)1; 1122 *xc++ = y & FFFFFFFF; 1123 } 1124 #else 1125 #ifdef Pack_32 1126 do { 1127 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow; 1128 borrow = (y & 0x10000) >> 16; 1129 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow; 1130 borrow = (z & 0x10000) >> 16; 1131 Storeinc(xc, z, y); 1132 } 1133 while(xb < xbe); 1134 while(xa < xae) { 1135 y = (*xa & 0xffff) - borrow; 1136 borrow = (y & 0x10000) >> 16; 1137 z = (*xa++ >> 16) - borrow; 1138 borrow = (z & 0x10000) >> 16; 1139 Storeinc(xc, z, y); 1140 } 1141 #else 1142 do { 1143 y = *xa++ - *xb++ - borrow; 1144 borrow = (y & 0x10000) >> 16; 1145 *xc++ = y & 0xffff; 1146 } 1147 while(xb < xbe); 1148 while(xa < xae) { 1149 y = *xa++ - borrow; 1150 borrow = (y & 0x10000) >> 16; 1151 *xc++ = y & 0xffff; 1152 } 1153 #endif 1154 #endif 1155 while(!*--xc) 1156 wa--; 1157 c->wds = wa; 1158 return c; 1159 } 1160 1161 static double 1162 ulp 1163 #ifdef KR_headers 1164 (x) U *x; 1165 #else 1166 (U *x) 1167 #endif 1168 { 1169 Long L; 1170 U u; 1171 1172 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; 1173 #ifndef Avoid_Underflow 1174 #ifndef Sudden_Underflow 1175 if (L > 0) { 1176 #endif 1177 #endif 1178 #ifdef IBM 1179 L |= Exp_msk1 >> 4; 1180 #endif 1181 word0(&u) = L; 1182 word1(&u) = 0; 1183 #ifndef Avoid_Underflow 1184 #ifndef Sudden_Underflow 1185 } 1186 else { 1187 L = -L >> Exp_shift; 1188 if (L < Exp_shift) { 1189 word0(&u) = 0x80000 >> L; 1190 word1(&u) = 0; 1191 } 1192 else { 1193 word0(&u) = 0; 1194 L -= Exp_shift; 1195 word1(&u) = L >= 31 ? 1 : 1 << 31 - L; 1196 } 1197 } 1198 #endif 1199 #endif 1200 return dval(&u); 1201 } 1202 1203 static double 1204 b2d 1205 #ifdef KR_headers 1206 (a, e) Bigint *a; int *e; 1207 #else 1208 (Bigint *a, int *e) 1209 #endif 1210 { 1211 ULong *xa, *xa0, w, y, z; 1212 int k; 1213 U d; 1214 #ifdef VAX 1215 ULong d0, d1; 1216 #else 1217 #define d0 word0(&d) 1218 #define d1 word1(&d) 1219 #endif 1220 1221 xa0 = a->x; 1222 xa = xa0 + a->wds; 1223 y = *--xa; 1224 #ifdef DEBUG 1225 if (!y) Bug("zero y in b2d"); 1226 #endif 1227 k = hi0bits(y); 1228 *e = 32 - k; 1229 #ifdef Pack_32 1230 if (k < Ebits) { 1231 d0 = Exp_1 | y >> (Ebits - k); 1232 w = xa > xa0 ? *--xa : 0; 1233 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k); 1234 goto ret_d; 1235 } 1236 z = xa > xa0 ? *--xa : 0; 1237 if (k -= Ebits) { 1238 d0 = Exp_1 | y << k | z >> (32 - k); 1239 y = xa > xa0 ? *--xa : 0; 1240 d1 = z << k | y >> (32 - k); 1241 } 1242 else { 1243 d0 = Exp_1 | y; 1244 d1 = z; 1245 } 1246 #else 1247 if (k < Ebits + 16) { 1248 z = xa > xa0 ? *--xa : 0; 1249 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; 1250 w = xa > xa0 ? *--xa : 0; 1251 y = xa > xa0 ? *--xa : 0; 1252 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; 1253 goto ret_d; 1254 } 1255 z = xa > xa0 ? *--xa : 0; 1256 w = xa > xa0 ? *--xa : 0; 1257 k -= Ebits + 16; 1258 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; 1259 y = xa > xa0 ? *--xa : 0; 1260 d1 = w << k + 16 | y << k; 1261 #endif 1262 ret_d: 1263 #ifdef VAX 1264 word0(&d) = d0 >> 16 | d0 << 16; 1265 word1(&d) = d1 >> 16 | d1 << 16; 1266 #else 1267 #undef d0 1268 #undef d1 1269 #endif 1270 return dval(&d); 1271 } 1272 1273 static Bigint * 1274 d2b 1275 #ifdef KR_headers 1276 (d, e, bits) U *d; int *e, *bits; 1277 #else 1278 (U *d, int *e, int *bits) 1279 #endif 1280 { 1281 Bigint *b; 1282 int de, k; 1283 ULong *x, y, z; 1284 #ifndef Sudden_Underflow 1285 int i; 1286 #endif 1287 #ifdef VAX 1288 ULong d0, d1; 1289 d0 = word0(d) >> 16 | word0(d) << 16; 1290 d1 = word1(d) >> 16 | word1(d) << 16; 1291 #else 1292 #define d0 word0(d) 1293 #define d1 word1(d) 1294 #endif 1295 1296 #ifdef Pack_32 1297 b = Balloc(1); 1298 #else 1299 b = Balloc(2); 1300 #endif 1301 x = b->x; 1302 1303 z = d0 & Frac_mask; 1304 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 1305 #ifdef Sudden_Underflow 1306 de = (int)(d0 >> Exp_shift); 1307 #ifndef IBM 1308 z |= Exp_msk11; 1309 #endif 1310 #else 1311 de = (int)(d0 >> Exp_shift); 1312 if (de) 1313 z |= Exp_msk1; 1314 #endif 1315 #ifdef Pack_32 1316 y = d1; 1317 if (y) { 1318 k = lo0bits(&y); 1319 if (k) { 1320 x[0] = y | z << (32 - k); 1321 z >>= k; 1322 } 1323 else 1324 x[0] = y; 1325 x[1] = z; 1326 b->wds = x[1] ? 2 : 1; 1327 #ifndef Sudden_Underflow 1328 i = b->wds; 1329 #endif 1330 } 1331 else { 1332 k = lo0bits(&z); 1333 x[0] = z; 1334 #ifndef Sudden_Underflow 1335 i = 1336 #endif 1337 b->wds = 1; 1338 k += 32; 1339 } 1340 #else 1341 if (y = d1) { 1342 if (k = lo0bits(&y)) 1343 if (k >= 16) { 1344 x[0] = y | z << 32 - k & 0xffff; 1345 x[1] = z >> k - 16 & 0xffff; 1346 x[2] = z >> k; 1347 i = 2; 1348 } 1349 else { 1350 x[0] = y & 0xffff; 1351 x[1] = y >> 16 | z << 16 - k & 0xffff; 1352 x[2] = z >> k & 0xffff; 1353 x[3] = z >> k+16; 1354 i = 3; 1355 } 1356 else { 1357 x[0] = y & 0xffff; 1358 x[1] = y >> 16; 1359 x[2] = z & 0xffff; 1360 x[3] = z >> 16; 1361 i = 3; 1362 } 1363 } 1364 else { 1365 #ifdef DEBUG 1366 if (!z) 1367 Bug("Zero passed to d2b"); 1368 #endif 1369 k = lo0bits(&z); 1370 if (k >= 16) { 1371 x[0] = z; 1372 i = 0; 1373 } 1374 else { 1375 x[0] = z & 0xffff; 1376 x[1] = z >> 16; 1377 i = 1; 1378 } 1379 k += 32; 1380 } 1381 while(!x[i]) 1382 --i; 1383 b->wds = i + 1; 1384 #endif 1385 #ifndef Sudden_Underflow 1386 if (de) { 1387 #endif 1388 #ifdef IBM 1389 *e = (de - Bias - (P-1) << 2) + k; 1390 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask); 1391 #else 1392 *e = de - Bias - (P-1) + k; 1393 *bits = P - k; 1394 #endif 1395 #ifndef Sudden_Underflow 1396 } 1397 else { 1398 *e = de - Bias - (P-1) + 1 + k; 1399 #ifdef Pack_32 1400 *bits = 32*i - hi0bits(x[i-1]); 1401 #else 1402 *bits = (i+2)*16 - hi0bits(x[i]); 1403 #endif 1404 } 1405 #endif 1406 return b; 1407 } 1408 #undef d0 1409 #undef d1 1410 1411 static double 1412 ratio 1413 #ifdef KR_headers 1414 (a, b) Bigint *a, *b; 1415 #else 1416 (Bigint *a, Bigint *b) 1417 #endif 1418 { 1419 U da, db; 1420 int k, ka, kb; 1421 1422 dval(&da) = b2d(a, &ka); 1423 dval(&db) = b2d(b, &kb); 1424 #ifdef Pack_32 1425 k = ka - kb + 32*(a->wds - b->wds); 1426 #else 1427 k = ka - kb + 16*(a->wds - b->wds); 1428 #endif 1429 #ifdef IBM 1430 if (k > 0) { 1431 word0(&da) += (k >> 2)*Exp_msk1; 1432 if (k &= 3) 1433 dval(&da) *= 1 << k; 1434 } 1435 else { 1436 k = -k; 1437 word0(&db) += (k >> 2)*Exp_msk1; 1438 if (k &= 3) 1439 dval(&db) *= 1 << k; 1440 } 1441 #else 1442 if (k > 0) 1443 word0(&da) += k*Exp_msk1; 1444 else { 1445 k = -k; 1446 word0(&db) += k*Exp_msk1; 1447 } 1448 #endif 1449 return dval(&da) / dval(&db); 1450 } 1451 1452 static CONST double 1453 tens[] = { 1454 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1455 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1456 1e20, 1e21, 1e22 1457 #ifdef VAX 1458 , 1e23, 1e24 1459 #endif 1460 }; 1461 1462 static CONST double 1463 #ifdef IEEE_Arith 1464 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; 1465 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1466 #ifdef Avoid_Underflow 1467 9007199254740992.*9007199254740992.e-256 1468 /* = 2^106 * 1e-256 */ 1469 #else 1470 1e-256 1471 #endif 1472 }; 1473 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */ 1474 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */ 1475 #define Scale_Bit 0x10 1476 #define n_bigtens 5 1477 #else 1478 #ifdef IBM 1479 bigtens[] = { 1e16, 1e32, 1e64 }; 1480 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 }; 1481 #define n_bigtens 3 1482 #else 1483 bigtens[] = { 1e16, 1e32 }; 1484 static CONST double tinytens[] = { 1e-16, 1e-32 }; 1485 #define n_bigtens 2 1486 #endif 1487 #endif 1488 1489 #undef Need_Hexdig 1490 #ifdef INFNAN_CHECK 1491 #ifndef No_Hex_NaN 1492 #define Need_Hexdig 1493 #endif 1494 #endif 1495 1496 #ifndef Need_Hexdig 1497 #ifndef NO_HEX_FP 1498 #define Need_Hexdig 1499 #endif 1500 #endif 1501 1502 #ifdef Need_Hexdig /*{*/ 1503 static unsigned char hexdig[256]; 1504 1505 static void 1506 #ifdef KR_headers 1507 htinit(h, s, inc) unsigned char *h; unsigned char *s; int inc; 1508 #else 1509 htinit(unsigned char *h, unsigned char *s, int inc) 1510 #endif 1511 { 1512 int i, j; 1513 for(i = 0; (j = s[i]) !=0; i++) 1514 h[j] = (unsigned char)(i + inc); 1515 } 1516 1517 static void 1518 #ifdef KR_headers 1519 hexdig_init() 1520 #else 1521 hexdig_init(void) 1522 #endif 1523 { 1524 #define USC (unsigned char *) 1525 htinit(hexdig, USC "0123456789", 0x10); 1526 htinit(hexdig, USC "abcdef", 0x10 + 10); 1527 htinit(hexdig, USC "ABCDEF", 0x10 + 10); 1528 } 1529 #endif /* } Need_Hexdig */ 1530 1531 #ifdef INFNAN_CHECK 1532 1533 #ifndef NAN_WORD0 1534 #define NAN_WORD0 0x7ff80000 1535 #endif 1536 1537 #ifndef NAN_WORD1 1538 #define NAN_WORD1 0 1539 #endif 1540 1541 static int 1542 match 1543 #ifdef KR_headers 1544 (sp, t) char **sp, *t; 1545 #else 1546 (CONST char **sp, CONST char *t) 1547 #endif 1548 { 1549 int c, d; 1550 CONST char *s = *sp; 1551 1552 for(d = *t++; d; d = *t++) { 1553 if ((c = *++s) >= 'A' && c <= 'Z') 1554 c += 'a' - 'A'; 1555 if (c != d) 1556 return 0; 1557 } 1558 *sp = s + 1; 1559 return 1; 1560 } 1561 1562 #ifndef No_Hex_NaN 1563 static void 1564 hexnan 1565 #ifdef KR_headers 1566 (rvp, sp) U *rvp; CONST char **sp; 1567 #else 1568 (U *rvp, CONST char **sp) 1569 #endif 1570 { 1571 ULong c, x[2]; 1572 CONST char *s; 1573 int c1, havedig, udx0, xshift; 1574 1575 if (!hexdig['0']) 1576 hexdig_init(); 1577 x[0] = x[1] = 0; 1578 havedig = xshift = 0; 1579 udx0 = 1; 1580 s = *sp; 1581 /* allow optional initial 0x or 0X */ 1582 for(c = *(CONST unsigned char*)(s+1); c && c <= ' '; c = *(CONST unsigned char*)(s+1)) 1583 ++s; 1584 if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')) 1585 s += 2; 1586 for(c = *(CONST unsigned char*)++s; c; c = *(CONST unsigned char*)++s) { 1587 c1 = hexdig[c]; 1588 if (c1) 1589 c = c1 & 0xf; 1590 else if (c <= ' ') { 1591 if (udx0 && havedig) { 1592 udx0 = 0; 1593 xshift = 1; 1594 } 1595 continue; 1596 } 1597 #ifdef GDTOA_NON_PEDANTIC_NANCHECK 1598 else if (/*(*/ c == ')' && havedig) { 1599 *sp = s + 1; 1600 break; 1601 } 1602 else 1603 return; /* invalid form: don't change *sp */ 1604 #else 1605 else { 1606 do { 1607 if (/*(*/ c == ')') { 1608 *sp = s + 1; 1609 break; 1610 } 1611 c = *++s; 1612 } while(c); 1613 break; 1614 } 1615 #endif 1616 havedig = 1; 1617 if (xshift) { 1618 xshift = 0; 1619 x[0] = x[1]; 1620 x[1] = 0; 1621 } 1622 if (udx0) 1623 x[0] = (x[0] << 4) | (x[1] >> 28); 1624 x[1] = (x[1] << 4) | c; 1625 } 1626 if ((x[0] &= 0xfffff) || x[1]) { 1627 word0(rvp) = Exp_mask | x[0]; 1628 word1(rvp) = x[1]; 1629 } 1630 } 1631 #endif /*No_Hex_NaN*/ 1632 #endif /* INFNAN_CHECK */ 1633 1634 #ifdef Pack_32 1635 #define ULbits 32 1636 #define kshift 5 1637 #define kmask 31 1638 #else 1639 #define ULbits 16 1640 #define kshift 4 1641 #define kmask 15 1642 #endif 1643 #ifndef NO_HEX_FP /*{*/ 1644 1645 static void 1646 #ifdef KR_headers 1647 rshift(b, k) Bigint *b; int k; 1648 #else 1649 rshift(Bigint *b, int k) 1650 #endif 1651 { 1652 ULong *x, *x1, *xe, y; 1653 int n; 1654 1655 x = x1 = b->x; 1656 n = k >> kshift; 1657 if (n < b->wds) { 1658 xe = x + b->wds; 1659 x += n; 1660 if (k &= kmask) { 1661 n = 32 - k; 1662 y = *x++ >> k; 1663 while(x < xe) { 1664 *x1++ = (y | (*x << n)) & 0xffffffff; 1665 y = *x++ >> k; 1666 } 1667 if ((*x1 = y) !=0) 1668 x1++; 1669 } 1670 else 1671 while(x < xe) 1672 *x1++ = *x++; 1673 } 1674 if ((b->wds = x1 - b->x) == 0) 1675 b->x[0] = 0; 1676 } 1677 1678 static ULong 1679 #ifdef KR_headers 1680 any_on(b, k) Bigint *b; int k; 1681 #else 1682 any_on(Bigint *b, int k) 1683 #endif 1684 { 1685 int n, nwds; 1686 ULong *x, *x0, x1, x2; 1687 1688 x = b->x; 1689 nwds = b->wds; 1690 n = k >> kshift; 1691 if (n > nwds) 1692 n = nwds; 1693 else if (n < nwds && (k &= kmask)) { 1694 x1 = x2 = x[n]; 1695 x1 >>= k; 1696 x1 <<= k; 1697 if (x1 != x2) 1698 return 1; 1699 } 1700 x0 = x; 1701 x += n; 1702 while(x > x0) 1703 if (*--x) 1704 return 1; 1705 return 0; 1706 } 1707 1708 enum { /* rounding values: same as FLT_ROUNDS */ 1709 Round_zero = 0, 1710 Round_near = 1, 1711 Round_up = 2, 1712 Round_down = 3 1713 }; 1714 1715 static Bigint * 1716 #ifdef KR_headers 1717 increment(b) Bigint *b; 1718 #else 1719 increment(Bigint *b) 1720 #endif 1721 { 1722 ULong *x, *xe; 1723 Bigint *b1; 1724 1725 x = b->x; 1726 xe = x + b->wds; 1727 do { 1728 if (*x < (ULong)0xffffffffL) { 1729 ++*x; 1730 return b; 1731 } 1732 *x++ = 0; 1733 } while(x < xe); 1734 { 1735 if (b->wds >= b->maxwds) { 1736 b1 = Balloc(b->k+1); 1737 Bcopy(b1,b); 1738 Bfree(b); 1739 b = b1; 1740 } 1741 b->x[b->wds++] = 1; 1742 } 1743 return b; 1744 } 1745 1746 void 1747 #ifdef KR_headers 1748 gethex(sp, rvp, rounding, sign) 1749 CONST char **sp; U *rvp; int rounding, sign; 1750 #else 1751 gethex( CONST char **sp, U *rvp, int rounding, int sign) 1752 #endif 1753 { 1754 Bigint *b; 1755 CONST unsigned char *decpt, *s0, *s, *s1; 1756 Long e, e1; 1757 ULong L, lostbits, *x; 1758 int big, denorm, esign, havedig, k, n, nbits, up, zret; 1759 #ifdef IBM 1760 int j; 1761 #endif 1762 enum { 1763 #ifdef IEEE_Arith /*{{*/ 1764 emax = 0x7fe - Bias - P + 1, 1765 emin = Emin - P + 1 1766 #else /*}{*/ 1767 emin = Emin - P, 1768 #ifdef VAX 1769 emax = 0x7ff - Bias - P + 1 1770 #endif 1771 #ifdef IBM 1772 emax = 0x7f - Bias - P 1773 #endif 1774 #endif /*}}*/ 1775 }; 1776 #ifdef USE_LOCALE 1777 int i; 1778 #ifdef NO_LOCALE_CACHE 1779 const unsigned char *decimalpoint = (unsigned char*) 1780 localeconv()->decimal_point; 1781 #else 1782 const unsigned char *decimalpoint; 1783 static unsigned char *decimalpoint_cache; 1784 if (!(s0 = decimalpoint_cache)) { 1785 s0 = (unsigned char*)localeconv()->decimal_point; 1786 if ((decimalpoint_cache = (unsigned char*) 1787 MALLOC(strlen((CONST char*)s0) + 1))) { 1788 strcpy((char*)decimalpoint_cache, (CONST char*)s0); 1789 s0 = decimalpoint_cache; 1790 } 1791 } 1792 decimalpoint = s0; 1793 #endif 1794 #endif 1795 1796 if (!hexdig['0']) 1797 hexdig_init(); 1798 havedig = 0; 1799 s0 = *(CONST unsigned char **)sp + 2; 1800 while(s0[havedig] == '0') 1801 havedig++; 1802 s0 += havedig; 1803 s = s0; 1804 decpt = 0; 1805 zret = 0; 1806 e = 0; 1807 if (hexdig[*s]) 1808 havedig++; 1809 else { 1810 zret = 1; 1811 #ifdef USE_LOCALE 1812 for(i = 0; decimalpoint[i]; ++i) { 1813 if (s[i] != decimalpoint[i]) 1814 goto pcheck; 1815 } 1816 decpt = s += i; 1817 #else 1818 if (*s != '.') 1819 goto pcheck; 1820 decpt = ++s; 1821 #endif 1822 if (!hexdig[*s]) 1823 goto pcheck; 1824 while(*s == '0') 1825 s++; 1826 if (hexdig[*s]) 1827 zret = 0; 1828 havedig = 1; 1829 s0 = s; 1830 } 1831 while(hexdig[*s]) 1832 s++; 1833 #ifdef USE_LOCALE 1834 if (*s == *decimalpoint && !decpt) { 1835 for(i = 1; decimalpoint[i]; ++i) { 1836 if (s[i] != decimalpoint[i]) 1837 goto pcheck; 1838 } 1839 decpt = s += i; 1840 #else 1841 if (*s == '.' && !decpt) { 1842 decpt = ++s; 1843 #endif 1844 while(hexdig[*s]) 1845 s++; 1846 }/*}*/ 1847 if (decpt) 1848 e = -(((Long)(s-decpt)) << 2); 1849 pcheck: 1850 s1 = s; 1851 big = esign = 0; 1852 switch(*s) { 1853 case 'p': 1854 case 'P': 1855 switch(*++s) { 1856 case '-': 1857 esign = 1; 1858 /* no break */ 1859 case '+': 1860 s++; 1861 } 1862 if ((n = hexdig[*s]) == 0 || n > 0x19) { 1863 s = s1; 1864 break; 1865 } 1866 e1 = n - 0x10; 1867 while((n = hexdig[*++s]) !=0 && n <= 0x19) { 1868 if (e1 & 0xf8000000) 1869 big = 1; 1870 e1 = 10*e1 + n - 0x10; 1871 } 1872 if (esign) 1873 e1 = -e1; 1874 e += e1; 1875 } 1876 *sp = (char*)s; 1877 if (!havedig) 1878 *sp = (char*)s0 - 1; 1879 if (zret) 1880 goto retz1; 1881 if (big) { 1882 if (esign) { 1883 #ifdef IEEE_Arith 1884 switch(rounding) { 1885 case Round_up: 1886 if (sign) 1887 break; 1888 goto ret_tiny; 1889 case Round_down: 1890 if (!sign) 1891 break; 1892 goto ret_tiny; 1893 } 1894 #endif 1895 goto retz; 1896 #ifdef IEEE_Arith 1897 ret_tiny: 1898 #ifndef NO_ERRNO 1899 errno = ERANGE; 1900 #endif 1901 word0(rvp) = 0; 1902 word1(rvp) = 1; 1903 return; 1904 #endif /* IEEE_Arith */ 1905 } 1906 switch(rounding) { 1907 case Round_near: 1908 goto ovfl1; 1909 case Round_up: 1910 if (!sign) 1911 goto ovfl1; 1912 goto ret_big; 1913 case Round_down: 1914 if (sign) 1915 goto ovfl1; 1916 goto ret_big; 1917 } 1918 ret_big: 1919 word0(rvp) = Big0; 1920 word1(rvp) = Big1; 1921 return; 1922 } 1923 n = s1 - s0 - 1; 1924 for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1) 1925 k++; 1926 b = Balloc(k); 1927 x = b->x; 1928 n = 0; 1929 L = 0; 1930 #ifdef USE_LOCALE 1931 for(i = 0; decimalpoint[i+1]; ++i); 1932 #endif 1933 while(s1 > s0) { 1934 #ifdef USE_LOCALE 1935 if (*--s1 == decimalpoint[i]) { 1936 s1 -= i; 1937 continue; 1938 } 1939 #else 1940 if (*--s1 == '.') 1941 continue; 1942 #endif 1943 if (n == ULbits) { 1944 *x++ = L; 1945 L = 0; 1946 n = 0; 1947 } 1948 L |= (hexdig[*s1] & 0x0f) << n; 1949 n += 4; 1950 } 1951 *x++ = L; 1952 b->wds = n = x - b->x; 1953 n = ULbits*n - hi0bits(L); 1954 nbits = Nbits; 1955 lostbits = 0; 1956 x = b->x; 1957 if (n > nbits) { 1958 n -= nbits; 1959 if (any_on(b,n)) { 1960 lostbits = 1; 1961 k = n - 1; 1962 if (x[k>>kshift] & 1 << (k & kmask)) { 1963 lostbits = 2; 1964 if (k > 0 && any_on(b,k)) 1965 lostbits = 3; 1966 } 1967 } 1968 rshift(b, n); 1969 e += n; 1970 } 1971 else if (n < nbits) { 1972 n = nbits - n; 1973 b = lshift(b, n); 1974 e -= n; 1975 x = b->x; 1976 } 1977 if (e > Emax) { 1978 ovfl: 1979 Bfree(b); 1980 ovfl1: 1981 #ifndef NO_ERRNO 1982 errno = ERANGE; 1983 #endif 1984 word0(rvp) = Exp_mask; 1985 word1(rvp) = 0; 1986 return; 1987 } 1988 denorm = 0; 1989 if (e < emin) { 1990 denorm = 1; 1991 n = emin - e; 1992 if (n >= nbits) { 1993 #ifdef IEEE_Arith /*{*/ 1994 switch (rounding) { 1995 case Round_near: 1996 if (n == nbits && (n < 2 || any_on(b,n-1))) 1997 goto ret_tiny; 1998 break; 1999 case Round_up: 2000 if (!sign) 2001 goto ret_tiny; 2002 break; 2003 case Round_down: 2004 if (sign) 2005 goto ret_tiny; 2006 } 2007 #endif /* } IEEE_Arith */ 2008 Bfree(b); 2009 retz: 2010 #ifndef NO_ERRNO 2011 errno = ERANGE; 2012 #endif 2013 retz1: 2014 rvp->d = 0.; 2015 return; 2016 } 2017 k = n - 1; 2018 if (lostbits) 2019 lostbits = 1; 2020 else if (k > 0) 2021 lostbits = any_on(b,k); 2022 if (x[k>>kshift] & 1 << (k & kmask)) 2023 lostbits |= 2; 2024 nbits -= n; 2025 rshift(b,n); 2026 e = emin; 2027 } 2028 if (lostbits) { 2029 up = 0; 2030 switch(rounding) { 2031 case Round_zero: 2032 break; 2033 case Round_near: 2034 if (lostbits & 2 2035 && (lostbits & 1) | (x[0] & 1)) 2036 up = 1; 2037 break; 2038 case Round_up: 2039 up = 1 - sign; 2040 break; 2041 case Round_down: 2042 up = sign; 2043 } 2044 if (up) { 2045 k = b->wds; 2046 b = increment(b); 2047 x = b->x; 2048 if (denorm) { 2049 #if 0 2050 if (nbits == Nbits - 1 2051 && x[nbits >> kshift] & 1 << (nbits & kmask)) 2052 denorm = 0; /* not currently used */ 2053 #endif 2054 } 2055 else if (b->wds > k 2056 || ((n = nbits & kmask) !=0 2057 && hi0bits(x[k-1]) < 32-n)) { 2058 rshift(b,1); 2059 if (++e > Emax) 2060 goto ovfl; 2061 } 2062 } 2063 } 2064 #ifdef IEEE_Arith 2065 if (denorm) 2066 word0(rvp) = b->wds > 1 ? b->x[1] & ~0x100000 : 0; 2067 else 2068 word0(rvp) = (b->x[1] & ~0x100000) | ((e + 0x3ff + 52) << 20); 2069 word1(rvp) = b->x[0]; 2070 #endif 2071 #ifdef IBM 2072 if ((j = e & 3)) { 2073 k = b->x[0] & ((1 << j) - 1); 2074 rshift(b,j); 2075 if (k) { 2076 switch(rounding) { 2077 case Round_up: 2078 if (!sign) 2079 increment(b); 2080 break; 2081 case Round_down: 2082 if (sign) 2083 increment(b); 2084 break; 2085 case Round_near: 2086 j = 1 << (j-1); 2087 if (k & j && ((k & (j-1)) | lostbits)) 2088 increment(b); 2089 } 2090 } 2091 } 2092 e >>= 2; 2093 word0(rvp) = b->x[1] | ((e + 65 + 13) << 24); 2094 word1(rvp) = b->x[0]; 2095 #endif 2096 #ifdef VAX 2097 /* The next two lines ignore swap of low- and high-order 2 bytes. */ 2098 /* word0(rvp) = (b->x[1] & ~0x800000) | ((e + 129 + 55) << 23); */ 2099 /* word1(rvp) = b->x[0]; */ 2100 word0(rvp) = ((b->x[1] & ~0x800000) >> 16) | ((e + 129 + 55) << 7) | (b->x[1] << 16); 2101 word1(rvp) = (b->x[0] >> 16) | (b->x[0] << 16); 2102 #endif 2103 Bfree(b); 2104 } 2105 #endif /*}!NO_HEX_FP*/ 2106 2107 static int 2108 #ifdef KR_headers 2109 dshift(b, p2) Bigint *b; int p2; 2110 #else 2111 dshift(Bigint *b, int p2) 2112 #endif 2113 { 2114 int rv = hi0bits(b->x[b->wds-1]) - 4; 2115 if (p2 > 0) 2116 rv -= p2; 2117 return rv & kmask; 2118 } 2119 2120 static int 2121 quorem 2122 #ifdef KR_headers 2123 (b, S) Bigint *b, *S; 2124 #else 2125 (Bigint *b, Bigint *S) 2126 #endif 2127 { 2128 int n; 2129 ULong *bx, *bxe, q, *sx, *sxe; 2130 #ifdef ULLong 2131 ULLong borrow, carry, y, ys; 2132 #else 2133 ULong borrow, carry, y, ys; 2134 #ifdef Pack_32 2135 ULong si, z, zs; 2136 #endif 2137 #endif 2138 2139 n = S->wds; 2140 #ifdef DEBUG 2141 /*debug*/ if (b->wds > n) 2142 /*debug*/ Bug("oversize b in quorem"); 2143 #endif 2144 if (b->wds < n) 2145 return 0; 2146 sx = S->x; 2147 sxe = sx + --n; 2148 bx = b->x; 2149 bxe = bx + n; 2150 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ 2151 #ifdef DEBUG 2152 /*debug*/ if (q > 9) 2153 /*debug*/ Bug("oversized quotient in quorem"); 2154 #endif 2155 if (q) { 2156 borrow = 0; 2157 carry = 0; 2158 do { 2159 #ifdef ULLong 2160 ys = *sx++ * (ULLong)q + carry; 2161 carry = ys >> 32; 2162 y = *bx - (ys & FFFFFFFF) - borrow; 2163 borrow = y >> 32 & (ULong)1; 2164 *bx++ = y & FFFFFFFF; 2165 #else 2166 #ifdef Pack_32 2167 si = *sx++; 2168 ys = (si & 0xffff) * q + carry; 2169 zs = (si >> 16) * q + (ys >> 16); 2170 carry = zs >> 16; 2171 y = (*bx & 0xffff) - (ys & 0xffff) - borrow; 2172 borrow = (y & 0x10000) >> 16; 2173 z = (*bx >> 16) - (zs & 0xffff) - borrow; 2174 borrow = (z & 0x10000) >> 16; 2175 Storeinc(bx, z, y); 2176 #else 2177 ys = *sx++ * q + carry; 2178 carry = ys >> 16; 2179 y = *bx - (ys & 0xffff) - borrow; 2180 borrow = (y & 0x10000) >> 16; 2181 *bx++ = y & 0xffff; 2182 #endif 2183 #endif 2184 } 2185 while(sx <= sxe); 2186 if (!*bxe) { 2187 bx = b->x; 2188 while(--bxe > bx && !*bxe) 2189 --n; 2190 b->wds = n; 2191 } 2192 } 2193 if (cmp(b, S) >= 0) { 2194 q++; 2195 borrow = 0; 2196 carry = 0; 2197 bx = b->x; 2198 sx = S->x; 2199 do { 2200 #ifdef ULLong 2201 ys = *sx++ + carry; 2202 carry = ys >> 32; 2203 y = *bx - (ys & FFFFFFFF) - borrow; 2204 borrow = y >> 32 & (ULong)1; 2205 *bx++ = y & FFFFFFFF; 2206 #else 2207 #ifdef Pack_32 2208 si = *sx++; 2209 ys = (si & 0xffff) + carry; 2210 zs = (si >> 16) + (ys >> 16); 2211 carry = zs >> 16; 2212 y = (*bx & 0xffff) - (ys & 0xffff) - borrow; 2213 borrow = (y & 0x10000) >> 16; 2214 z = (*bx >> 16) - (zs & 0xffff) - borrow; 2215 borrow = (z & 0x10000) >> 16; 2216 Storeinc(bx, z, y); 2217 #else 2218 ys = *sx++ + carry; 2219 carry = ys >> 16; 2220 y = *bx - (ys & 0xffff) - borrow; 2221 borrow = (y & 0x10000) >> 16; 2222 *bx++ = y & 0xffff; 2223 #endif 2224 #endif 2225 } 2226 while(sx <= sxe); 2227 bx = b->x; 2228 bxe = bx + n; 2229 if (!*bxe) { 2230 while(--bxe > bx && !*bxe) 2231 --n; 2232 b->wds = n; 2233 } 2234 } 2235 return q; 2236 } 2237 2238 #ifndef NO_STRTOD_BIGCOMP 2239 2240 static void 2241 bigcomp 2242 #ifdef KR_headers 2243 (rv, s0, bc) 2244 U *rv; CONST char *s0; BCinfo *bc; 2245 #else 2246 (U *rv, CONST char *s0, BCinfo *bc) 2247 #endif 2248 { 2249 Bigint *b, *d; 2250 int b2, bbits, d2, dd, dig, dsign, i, j, nd, nd0, p2, p5, speccase; 2251 2252 dsign = bc->dsign; 2253 nd = bc->nd; 2254 nd0 = bc->nd0; 2255 p5 = nd + bc->e0 - 1; 2256 dd = speccase = 0; 2257 #ifndef Sudden_Underflow 2258 if (rv->d == 0.) { /* special case: value near underflow-to-zero */ 2259 /* threshold was rounded to zero */ 2260 b = i2b(1); 2261 p2 = Emin - P + 1; 2262 bbits = 1; 2263 #ifdef Avoid_Underflow 2264 word0(rv) = (P+2) << Exp_shift; 2265 #else 2266 word1(rv) = 1; 2267 #endif 2268 i = 0; 2269 #ifdef Honor_FLT_ROUNDS 2270 if (bc->rounding == 1) 2271 #endif 2272 { 2273 speccase = 1; 2274 --p2; 2275 dsign = 0; 2276 goto have_i; 2277 } 2278 } 2279 else 2280 #endif 2281 b = d2b(rv, &p2, &bbits); 2282 #ifdef Avoid_Underflow 2283 p2 -= bc->scale; 2284 #endif 2285 /* floor(log2(rv)) == bbits - 1 + p2 */ 2286 /* Check for denormal case. */ 2287 i = P - bbits; 2288 if (i > (j = P - Emin - 1 + p2)) { 2289 #ifdef Sudden_Underflow 2290 Bfree(b); 2291 b = i2b(1); 2292 p2 = Emin; 2293 i = P - 1; 2294 #ifdef Avoid_Underflow 2295 word0(rv) = (1 + bc->scale) << Exp_shift; 2296 #else 2297 word0(rv) = Exp_msk1; 2298 #endif 2299 word1(rv) = 0; 2300 #else 2301 i = j; 2302 #endif 2303 } 2304 #ifdef Honor_FLT_ROUNDS 2305 if (bc->rounding != 1) { 2306 if (i > 0) 2307 b = lshift(b, i); 2308 if (dsign) 2309 b = increment(b); 2310 } 2311 else 2312 #endif 2313 { 2314 b = lshift(b, ++i); 2315 b->x[0] |= 1; 2316 } 2317 #ifndef Sudden_Underflow 2318 have_i: 2319 #endif 2320 p2 -= p5 + i; 2321 d = i2b(1); 2322 /* Arrange for convenient computation of quotients: 2323 * shift left if necessary so divisor has 4 leading 0 bits. 2324 */ 2325 if (p5 > 0) 2326 d = pow5mult(d, p5); 2327 else if (p5 < 0) 2328 b = pow5mult(b, -p5); 2329 if (p2 > 0) { 2330 b2 = p2; 2331 d2 = 0; 2332 } 2333 else { 2334 b2 = 0; 2335 d2 = -p2; 2336 } 2337 i = dshift(d, d2); 2338 if ((b2 += i) > 0) 2339 b = lshift(b, b2); 2340 if ((d2 += i) > 0) 2341 d = lshift(d, d2); 2342 2343 /* Now b/d = exactly half-way between the two floating-point values */ 2344 /* on either side of the input string. Compute first digit of b/d. */ 2345 2346 dig = quorem(b,d); 2347 if (!dig) { 2348 b = multadd(b, 10, 0); /* very unlikely */ 2349 dig = quorem(b,d); 2350 } 2351 2352 /* Compare b/d with s0 */ 2353 2354 for(i = 0; i < nd0; ) { 2355 dd = s0[i++] - '0' - dig; 2356 if (dd) 2357 goto ret; 2358 if (!b->x[0] && b->wds == 1) { 2359 if (i < nd) 2360 dd = 1; 2361 goto ret; 2362 } 2363 b = multadd(b, 10, 0); 2364 dig = quorem(b,d); 2365 } 2366 for(j = bc->dp1; i++ < nd;) { 2367 dd = s0[j++] - '0' - dig; 2368 if (dd) 2369 goto ret; 2370 if (!b->x[0] && b->wds == 1) { 2371 if (i < nd) 2372 dd = 1; 2373 goto ret; 2374 } 2375 b = multadd(b, 10, 0); 2376 dig = quorem(b,d); 2377 } 2378 if (b->x[0] || b->wds > 1) 2379 dd = -1; 2380 ret: 2381 Bfree(b); 2382 Bfree(d); 2383 #ifdef Honor_FLT_ROUNDS 2384 if (bc->rounding != 1) { 2385 if (dd < 0) { 2386 if (bc->rounding == 0) { 2387 if (!dsign) 2388 goto retlow1; 2389 } 2390 else if (dsign) 2391 goto rethi1; 2392 } 2393 else if (dd > 0) { 2394 if (bc->rounding == 0) { 2395 if (dsign) 2396 goto rethi1; 2397 goto ret1; 2398 } 2399 if (!dsign) 2400 goto rethi1; 2401 dval(rv) += 2.*ulp(rv); 2402 } 2403 else { 2404 bc->inexact = 0; 2405 if (dsign) 2406 goto rethi1; 2407 } 2408 } 2409 else 2410 #endif 2411 if (speccase) { 2412 if (dd <= 0) 2413 rv->d = 0.; 2414 } 2415 else if (dd < 0) { 2416 if (!dsign) /* does not happen for round-near */ 2417 retlow1: 2418 dval(rv) -= ulp(rv); 2419 } 2420 else if (dd > 0) { 2421 if (dsign) { 2422 rethi1: 2423 dval(rv) += ulp(rv); 2424 } 2425 } 2426 else { 2427 /* Exact half-way case: apply round-even rule. */ 2428 if (word1(rv) & 1) { 2429 if (dsign) 2430 goto rethi1; 2431 goto retlow1; 2432 } 2433 } 2434 2435 #ifdef Honor_FLT_ROUNDS 2436 ret1: 2437 #endif 2438 return; 2439 } 2440 #endif /* NO_STRTOD_BIGCOMP */ 2441 2442 double 2443 strtod 2444 #ifdef KR_headers 2445 (s00, se) CONST char *s00; char **se; 2446 #else 2447 (CONST char *s00, char **se) 2448 #endif 2449 { 2450 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, e, e1; 2451 int esign, i, j, k, nd, nd0, nf, nz, nz0, sign; 2452 CONST char *s, *s0, *s1; 2453 double aadj, aadj1; 2454 Long L; 2455 U aadj2, adj, rv, rv0; 2456 ULong y, z; 2457 BCinfo bc; 2458 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; 2459 #ifdef SET_INEXACT 2460 int oldinexact; 2461 #endif 2462 #ifdef Honor_FLT_ROUNDS /*{*/ 2463 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */ 2464 bc.rounding = Flt_Rounds; 2465 #else /*}{*/ 2466 bc.rounding = 1; 2467 switch(fegetround()) { 2468 case FE_TOWARDZERO: bc.rounding = 0; break; 2469 case FE_UPWARD: bc.rounding = 2; break; 2470 case FE_DOWNWARD: bc.rounding = 3; 2471 } 2472 #endif /*}}*/ 2473 #endif /*}*/ 2474 #ifdef USE_LOCALE 2475 CONST char *s2; 2476 #endif 2477 2478 sign = nz0 = nz = bc.dplen = bc.uflchk = 0; 2479 dval(&rv) = 0.; 2480 for(s = s00;;s++) switch(*s) { 2481 case '-': 2482 sign = 1; 2483 /* no break */ 2484 case '+': 2485 if (*++s) 2486 goto break2; 2487 /* no break */ 2488 case 0: 2489 goto ret0; 2490 case '\t': 2491 case '\n': 2492 case '\v': 2493 case '\f': 2494 case '\r': 2495 case ' ': 2496 continue; 2497 default: 2498 goto break2; 2499 } 2500 break2: 2501 if (*s == '0') { 2502 #ifndef NO_HEX_FP /*{*/ 2503 switch(s[1]) { 2504 case 'x': 2505 case 'X': 2506 #ifdef Honor_FLT_ROUNDS 2507 gethex(&s, &rv, bc.rounding, sign); 2508 #else 2509 gethex(&s, &rv, 1, sign); 2510 #endif 2511 goto ret; 2512 } 2513 #endif /*}*/ 2514 nz0 = 1; 2515 while(*++s == '0') ; 2516 if (!*s) 2517 goto ret; 2518 } 2519 s0 = s; 2520 y = z = 0; 2521 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 2522 if (nd < 9) 2523 y = 10*y + c - '0'; 2524 else if (nd < 16) 2525 z = 10*z + c - '0'; 2526 nd0 = nd; 2527 bc.dp0 = bc.dp1 = s - s0; 2528 #ifdef USE_LOCALE 2529 s1 = localeconv()->decimal_point; 2530 if (c == *s1) { 2531 c = '.'; 2532 if (*++s1) { 2533 s2 = s; 2534 for(;;) { 2535 if (*++s2 != *s1) { 2536 c = 0; 2537 break; 2538 } 2539 if (!*++s1) { 2540 s = s2; 2541 break; 2542 } 2543 } 2544 } 2545 } 2546 #endif 2547 if (c == '.') { 2548 c = *++s; 2549 bc.dp1 = s - s0; 2550 bc.dplen = bc.dp1 - bc.dp0; 2551 if (!nd) { 2552 for(; c == '0'; c = *++s) 2553 nz++; 2554 if (c > '0' && c <= '9') { 2555 s0 = s; 2556 nf += nz; 2557 nz = 0; 2558 goto have_dig; 2559 } 2560 goto dig_done; 2561 } 2562 for(; c >= '0' && c <= '9'; c = *++s) { 2563 have_dig: 2564 nz++; 2565 if (c -= '0') { 2566 nf += nz; 2567 for(i = 1; i < nz; i++) 2568 if (nd++ < 9) 2569 y *= 10; 2570 else if (nd <= DBL_DIG + 1) 2571 z *= 10; 2572 if (nd++ < 9) 2573 y = 10*y + c; 2574 else if (nd <= DBL_DIG + 1) 2575 z = 10*z + c; 2576 nz = 0; 2577 } 2578 } 2579 } 2580 dig_done: 2581 e = 0; 2582 if (c == 'e' || c == 'E') { 2583 if (!nd && !nz && !nz0) { 2584 goto ret0; 2585 } 2586 s00 = s; 2587 esign = 0; 2588 switch(c = *++s) { 2589 case '-': 2590 esign = 1; 2591 case '+': 2592 c = *++s; 2593 } 2594 if (c >= '0' && c <= '9') { 2595 while(c == '0') 2596 c = *++s; 2597 if (c > '0' && c <= '9') { 2598 L = c - '0'; 2599 s1 = s; 2600 while((c = *++s) >= '0' && c <= '9') 2601 L = 10*L + c - '0'; 2602 if (s - s1 > 8 || L > 19999) 2603 /* Avoid confusion from exponents 2604 * so large that e might overflow. 2605 */ 2606 e = 19999; /* safe for 16 bit ints */ 2607 else 2608 e = (int)L; 2609 if (esign) 2610 e = -e; 2611 } 2612 else 2613 e = 0; 2614 } 2615 else 2616 s = s00; 2617 } 2618 if (!nd) { 2619 if (!nz && !nz0) { 2620 #ifdef INFNAN_CHECK 2621 /* Check for Nan and Infinity */ 2622 if (!bc.dplen) 2623 switch(c) { 2624 case 'i': 2625 case 'I': 2626 if (match(&s,"nf")) { 2627 --s; 2628 if (!match(&s,"inity")) 2629 ++s; 2630 word0(&rv) = 0x7ff00000; 2631 word1(&rv) = 0; 2632 goto ret; 2633 } 2634 break; 2635 case 'n': 2636 case 'N': 2637 if (match(&s, "an")) { 2638 word0(&rv) = NAN_WORD0; 2639 word1(&rv) = NAN_WORD1; 2640 #ifndef No_Hex_NaN 2641 if (*s == '(') /*)*/ 2642 hexnan(&rv, &s); 2643 #endif 2644 goto ret; 2645 } 2646 } 2647 #endif /* INFNAN_CHECK */ 2648 ret0: 2649 s = s00; 2650 sign = 0; 2651 } 2652 goto ret; 2653 } 2654 bc.e0 = e1 = e -= nf; 2655 2656 /* Now we have nd0 digits, starting at s0, followed by a 2657 * decimal point, followed by nd-nd0 digits. The number we're 2658 * after is the integer represented by those digits times 2659 * 10**e */ 2660 2661 if (!nd0) 2662 nd0 = nd; 2663 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 2664 dval(&rv) = y; 2665 if (k > 9) { 2666 #ifdef SET_INEXACT 2667 if (k > DBL_DIG) 2668 oldinexact = get_inexact(); 2669 #endif 2670 dval(&rv) = tens[k - 9] * dval(&rv) + z; 2671 } 2672 bd0 = 0; 2673 if (nd <= DBL_DIG 2674 #ifndef RND_PRODQUOT 2675 #ifndef Honor_FLT_ROUNDS 2676 && Flt_Rounds == 1 2677 #endif 2678 #endif 2679 ) { 2680 if (!e) 2681 goto ret; 2682 if (e > 0) { 2683 if (e <= Ten_pmax) { 2684 #ifdef VAX 2685 goto vax_ovfl_check; 2686 #else 2687 #ifdef Honor_FLT_ROUNDS 2688 /* round correctly FLT_ROUNDS = 2 or 3 */ 2689 if (sign) { 2690 rv.d = -rv.d; 2691 sign = 0; 2692 } 2693 #endif 2694 /* rv = */ rounded_product(dval(&rv), tens[e]); 2695 goto ret; 2696 #endif 2697 } 2698 i = DBL_DIG - nd; 2699 if (e <= Ten_pmax + i) { 2700 /* A fancier test would sometimes let us do 2701 * this for larger i values. 2702 */ 2703 #ifdef Honor_FLT_ROUNDS 2704 /* round correctly FLT_ROUNDS = 2 or 3 */ 2705 if (sign) { 2706 rv.d = -rv.d; 2707 sign = 0; 2708 } 2709 #endif 2710 e -= i; 2711 dval(&rv) *= tens[i]; 2712 #ifdef VAX 2713 /* VAX exponent range is so narrow we must 2714 * worry about overflow here... 2715 */ 2716 vax_ovfl_check: 2717 word0(&rv) -= P*Exp_msk1; 2718 /* rv = */ rounded_product(dval(&rv), tens[e]); 2719 if ((word0(&rv) & Exp_mask) 2720 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 2721 goto ovfl; 2722 word0(&rv) += P*Exp_msk1; 2723 #else 2724 /* rv = */ rounded_product(dval(&rv), tens[e]); 2725 #endif 2726 goto ret; 2727 } 2728 } 2729 #ifndef Inaccurate_Divide 2730 else if (e >= -Ten_pmax) { 2731 #ifdef Honor_FLT_ROUNDS 2732 /* round correctly FLT_ROUNDS = 2 or 3 */ 2733 if (sign) { 2734 rv.d = -rv.d; 2735 sign = 0; 2736 } 2737 #endif 2738 /* rv = */ rounded_quotient(dval(&rv), tens[-e]); 2739 goto ret; 2740 } 2741 #endif 2742 } 2743 e1 += nd - k; 2744 2745 #ifdef IEEE_Arith 2746 #ifdef SET_INEXACT 2747 bc.inexact = 1; 2748 if (k <= DBL_DIG) 2749 oldinexact = get_inexact(); 2750 #endif 2751 #ifdef Avoid_Underflow 2752 bc.scale = 0; 2753 #endif 2754 #ifdef Honor_FLT_ROUNDS 2755 if (bc.rounding >= 2) { 2756 if (sign) 2757 bc.rounding = bc.rounding == 2 ? 0 : 2; 2758 else 2759 if (bc.rounding != 2) 2760 bc.rounding = 0; 2761 } 2762 #endif 2763 #endif /*IEEE_Arith*/ 2764 2765 /* Get starting approximation = rv * 10**e1 */ 2766 2767 if (e1 > 0) { 2768 i = e1 & 15; 2769 if (i) 2770 dval(&rv) *= tens[i]; 2771 if (e1 &= ~15) { 2772 if (e1 > DBL_MAX_10_EXP) { 2773 ovfl: 2774 #ifndef NO_ERRNO 2775 errno = ERANGE; 2776 #endif 2777 /* Can't trust HUGE_VAL */ 2778 #ifdef IEEE_Arith 2779 #ifdef Honor_FLT_ROUNDS 2780 switch(bc.rounding) { 2781 case 0: /* toward 0 */ 2782 case 3: /* toward -infinity */ 2783 word0(&rv) = Big0; 2784 word1(&rv) = Big1; 2785 break; 2786 default: 2787 word0(&rv) = Exp_mask; 2788 word1(&rv) = 0; 2789 } 2790 #else /*Honor_FLT_ROUNDS*/ 2791 word0(&rv) = Exp_mask; 2792 word1(&rv) = 0; 2793 #endif /*Honor_FLT_ROUNDS*/ 2794 #ifdef SET_INEXACT 2795 /* set overflow bit */ 2796 dval(&rv0) = 1e300; 2797 dval(&rv0) *= dval(&rv0); 2798 #endif 2799 #else /*IEEE_Arith*/ 2800 word0(&rv) = Big0; 2801 word1(&rv) = Big1; 2802 #endif /*IEEE_Arith*/ 2803 goto ret; 2804 } 2805 e1 >>= 4; 2806 for(j = 0; e1 > 1; j++, e1 >>= 1) 2807 if (e1 & 1) 2808 dval(&rv) *= bigtens[j]; 2809 /* The last multiplication could overflow. */ 2810 word0(&rv) -= P*Exp_msk1; 2811 dval(&rv) *= bigtens[j]; 2812 if ((z = word0(&rv) & Exp_mask) 2813 > Exp_msk1*(DBL_MAX_EXP+Bias-P)) 2814 goto ovfl; 2815 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { 2816 /* set to largest number */ 2817 /* (Can't trust DBL_MAX) */ 2818 word0(&rv) = Big0; 2819 word1(&rv) = Big1; 2820 } 2821 else 2822 word0(&rv) += P*Exp_msk1; 2823 } 2824 } 2825 else if (e1 < 0) { 2826 e1 = -e1; 2827 i = e1 & 15; 2828 if (i) 2829 dval(&rv) /= tens[i]; 2830 if (e1 >>= 4) { 2831 if (e1 >= 1 << n_bigtens) 2832 goto undfl; 2833 #ifdef Avoid_Underflow 2834 if (e1 & Scale_Bit) 2835 bc.scale = 2*P; 2836 for(j = 0; e1 > 0; j++, e1 >>= 1) 2837 if (e1 & 1) 2838 dval(&rv) *= tinytens[j]; 2839 if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask) 2840 >> Exp_shift)) > 0) { 2841 /* scaled rv is denormal; clear j low bits */ 2842 if (j >= 32) { 2843 word1(&rv) = 0; 2844 if (j >= 53) 2845 word0(&rv) = (P+2)*Exp_msk1; 2846 else 2847 word0(&rv) &= 0xffffffff << (j-32); 2848 } 2849 else 2850 word1(&rv) &= 0xffffffff << j; 2851 } 2852 #else 2853 for(j = 0; e1 > 1; j++, e1 >>= 1) 2854 if (e1 & 1) 2855 dval(&rv) *= tinytens[j]; 2856 /* The last multiplication could underflow. */ 2857 dval(&rv0) = dval(&rv); 2858 dval(&rv) *= tinytens[j]; 2859 if (!dval(&rv)) { 2860 dval(&rv) = 2.*dval(&rv0); 2861 dval(&rv) *= tinytens[j]; 2862 #endif 2863 if (!dval(&rv)) { 2864 undfl: 2865 dval(&rv) = 0.; 2866 #ifndef NO_ERRNO 2867 errno = ERANGE; 2868 #endif 2869 goto ret; 2870 } 2871 #ifndef Avoid_Underflow 2872 word0(&rv) = Tiny0; 2873 word1(&rv) = Tiny1; 2874 /* The refinement below will clean 2875 * this approximation up. 2876 */ 2877 } 2878 #endif 2879 } 2880 } 2881 2882 /* Now the hard part -- adjusting rv to the correct value.*/ 2883 2884 /* Put digits into bd: true value = bd * 10^e */ 2885 2886 bc.nd = nd; 2887 #ifndef NO_STRTOD_BIGCOMP 2888 bc.nd0 = nd0; /* Only needed if nd > strtod_diglim, but done here */ 2889 /* to silence an erroneous warning about bc.nd0 */ 2890 /* possibly not being initialized. */ 2891 if (nd > strtod_diglim) { 2892 /* ASSERT(strtod_diglim >= 18); 18 == one more than the */ 2893 /* minimum number of decimal digits to distinguish double values */ 2894 /* in IEEE arithmetic. */ 2895 i = j = 18; 2896 if (i > nd0) 2897 j += bc.dplen; 2898 for(;;) { 2899 if (--j <= bc.dp1 && j >= bc.dp0) 2900 j = bc.dp0 - 1; 2901 if (s0[j] != '0') 2902 break; 2903 --i; 2904 } 2905 e += nd - i; 2906 nd = i; 2907 if (nd0 > nd) 2908 nd0 = nd; 2909 if (nd < 9) { /* must recompute y */ 2910 y = 0; 2911 for(i = 0; i < nd0; ++i) 2912 y = 10*y + s0[i] - '0'; 2913 for(j = bc.dp1; i < nd; ++i) 2914 y = 10*y + s0[j++] - '0'; 2915 } 2916 } 2917 #endif 2918 bd0 = s2b(s0, nd0, nd, y, bc.dplen); 2919 2920 for(;;) { 2921 bd = Balloc(bd0->k); 2922 Bcopy(bd, bd0); 2923 bb = d2b(&rv, &bbe, &bbbits); /* rv = bb * 2^bbe */ 2924 bs = i2b(1); 2925 2926 if (e >= 0) { 2927 bb2 = bb5 = 0; 2928 bd2 = bd5 = e; 2929 } 2930 else { 2931 bb2 = bb5 = -e; 2932 bd2 = bd5 = 0; 2933 } 2934 if (bbe >= 0) 2935 bb2 += bbe; 2936 else 2937 bd2 -= bbe; 2938 bs2 = bb2; 2939 #ifdef Honor_FLT_ROUNDS 2940 if (bc.rounding != 1) 2941 bs2++; 2942 #endif 2943 #ifdef Avoid_Underflow 2944 j = bbe - bc.scale; 2945 i = j + bbbits - 1; /* logb(rv) */ 2946 if (i < Emin) /* denormal */ 2947 j += P - Emin; 2948 else 2949 j = P + 1 - bbbits; 2950 #else /*Avoid_Underflow*/ 2951 #ifdef Sudden_Underflow 2952 #ifdef IBM 2953 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 2954 #else 2955 j = P + 1 - bbbits; 2956 #endif 2957 #else /*Sudden_Underflow*/ 2958 j = bbe; 2959 i = j + bbbits - 1; /* logb(rv) */ 2960 if (i < Emin) /* denormal */ 2961 j += P - Emin; 2962 else 2963 j = P + 1 - bbbits; 2964 #endif /*Sudden_Underflow*/ 2965 #endif /*Avoid_Underflow*/ 2966 bb2 += j; 2967 bd2 += j; 2968 #ifdef Avoid_Underflow 2969 bd2 += bc.scale; 2970 #endif 2971 i = bb2 < bd2 ? bb2 : bd2; 2972 if (i > bs2) 2973 i = bs2; 2974 if (i > 0) { 2975 bb2 -= i; 2976 bd2 -= i; 2977 bs2 -= i; 2978 } 2979 if (bb5 > 0) { 2980 bs = pow5mult(bs, bb5); 2981 bb1 = mult(bs, bb); 2982 Bfree(bb); 2983 bb = bb1; 2984 } 2985 if (bb2 > 0) 2986 bb = lshift(bb, bb2); 2987 if (bd5 > 0) 2988 bd = pow5mult(bd, bd5); 2989 if (bd2 > 0) 2990 bd = lshift(bd, bd2); 2991 if (bs2 > 0) 2992 bs = lshift(bs, bs2); 2993 delta = diff(bb, bd); 2994 bc.dsign = delta->sign; 2995 delta->sign = 0; 2996 i = cmp(delta, bs); 2997 #ifndef NO_STRTOD_BIGCOMP 2998 if (bc.nd > nd && i <= 0) { 2999 if (bc.dsign) 3000 break; /* Must use bigcomp(). */ 3001 #ifdef Honor_FLT_ROUNDS 3002 if (bc.rounding != 1) { 3003 if (i < 0) 3004 break; 3005 } 3006 else 3007 #endif 3008 { 3009 bc.nd = nd; 3010 i = -1; /* Discarded digits make delta smaller. */ 3011 } 3012 } 3013 #endif 3014 #ifdef Honor_FLT_ROUNDS 3015 if (bc.rounding != 1) { 3016 if (i < 0) { 3017 /* Error is less than an ulp */ 3018 if (!delta->x[0] && delta->wds <= 1) { 3019 /* exact */ 3020 #ifdef SET_INEXACT 3021 bc.inexact = 0; 3022 #endif 3023 break; 3024 } 3025 if (bc.rounding) { 3026 if (bc.dsign) { 3027 adj.d = 1.; 3028 goto apply_adj; 3029 } 3030 } 3031 else if (!bc.dsign) { 3032 adj.d = -1.; 3033 if (!word1(&rv) 3034 && !(word0(&rv) & Frac_mask)) { 3035 y = word0(&rv) & Exp_mask; 3036 #ifdef Avoid_Underflow 3037 if (!bc.scale || y > 2*P*Exp_msk1) 3038 #else 3039 if (y) 3040 #endif 3041 { 3042 delta = lshift(delta,Log2P); 3043 if (cmp(delta, bs) <= 0) 3044 adj.d = -0.5; 3045 } 3046 } 3047 apply_adj: 3048 #ifdef Avoid_Underflow 3049 if (bc.scale && (y = word0(&rv) & Exp_mask) 3050 <= 2*P*Exp_msk1) 3051 word0(&adj) += (2*P+1)*Exp_msk1 - y; 3052 #else 3053 #ifdef Sudden_Underflow 3054 if ((word0(&rv) & Exp_mask) <= 3055 P*Exp_msk1) { 3056 word0(&rv) += P*Exp_msk1; 3057 dval(&rv) += adj.d*ulp(dval(&rv)); 3058 word0(&rv) -= P*Exp_msk1; 3059 } 3060 else 3061 #endif /*Sudden_Underflow*/ 3062 #endif /*Avoid_Underflow*/ 3063 dval(&rv) += adj.d*ulp(&rv); 3064 } 3065 break; 3066 } 3067 adj.d = ratio(delta, bs); 3068 if (adj.d < 1.) 3069 adj.d = 1.; 3070 if (adj.d <= 0x7ffffffe) { 3071 /* adj = rounding ? ceil(adj) : floor(adj); */ 3072 y = adj.d; 3073 if (y != adj.d) { 3074 if (!((bc.rounding>>1) ^ bc.dsign)) 3075 y++; 3076 adj.d = y; 3077 } 3078 } 3079 #ifdef Avoid_Underflow 3080 if (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) 3081 word0(&adj) += (2*P+1)*Exp_msk1 - y; 3082 #else 3083 #ifdef Sudden_Underflow 3084 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) { 3085 word0(&rv) += P*Exp_msk1; 3086 adj.d *= ulp(dval(&rv)); 3087 if (bc.dsign) 3088 dval(&rv) += adj.d; 3089 else 3090 dval(&rv) -= adj.d; 3091 word0(&rv) -= P*Exp_msk1; 3092 goto cont; 3093 } 3094 #endif /*Sudden_Underflow*/ 3095 #endif /*Avoid_Underflow*/ 3096 adj.d *= ulp(&rv); 3097 if (bc.dsign) { 3098 if (word0(&rv) == Big0 && word1(&rv) == Big1) 3099 goto ovfl; 3100 dval(&rv) += adj.d; 3101 } 3102 else 3103 dval(&rv) -= adj.d; 3104 goto cont; 3105 } 3106 #endif /*Honor_FLT_ROUNDS*/ 3107 3108 if (i < 0) { 3109 /* Error is less than half an ulp -- check for 3110 * special case of mantissa a power of two. 3111 */ 3112 if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask 3113 #ifdef IEEE_Arith 3114 #ifdef Avoid_Underflow 3115 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1 3116 #else 3117 || (word0(&rv) & Exp_mask) <= Exp_msk1 3118 #endif 3119 #endif 3120 ) { 3121 #ifdef SET_INEXACT 3122 if (!delta->x[0] && delta->wds <= 1) 3123 bc.inexact = 0; 3124 #endif 3125 break; 3126 } 3127 if (!delta->x[0] && delta->wds <= 1) { 3128 /* exact result */ 3129 #ifdef SET_INEXACT 3130 bc.inexact = 0; 3131 #endif 3132 break; 3133 } 3134 delta = lshift(delta,Log2P); 3135 if (cmp(delta, bs) > 0) 3136 goto drop_down; 3137 break; 3138 } 3139 if (i == 0) { 3140 /* exactly half-way between */ 3141 if (bc.dsign) { 3142 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1 3143 && word1(&rv) == ( 3144 #ifdef Avoid_Underflow 3145 (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) 3146 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : 3147 #endif 3148 0xffffffff)) { 3149 /*boundary case -- increment exponent*/ 3150 word0(&rv) = (word0(&rv) & Exp_mask) 3151 + Exp_msk1 3152 #ifdef IBM 3153 | Exp_msk1 >> 4 3154 #endif 3155 ; 3156 word1(&rv) = 0; 3157 #ifdef Avoid_Underflow 3158 bc.dsign = 0; 3159 #endif 3160 break; 3161 } 3162 } 3163 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) { 3164 drop_down: 3165 /* boundary case -- decrement exponent */ 3166 #ifdef Sudden_Underflow /*{{*/ 3167 L = word0(&rv) & Exp_mask; 3168 #ifdef IBM 3169 if (L < Exp_msk1) 3170 #else 3171 #ifdef Avoid_Underflow 3172 if (L <= (bc.scale ? (2*P+1)*Exp_msk1 : Exp_msk1)) 3173 #else 3174 if (L <= Exp_msk1) 3175 #endif /*Avoid_Underflow*/ 3176 #endif /*IBM*/ 3177 { 3178 if (bc.nd >nd) { 3179 bc.uflchk = 1; 3180 break; 3181 } 3182 goto undfl; 3183 } 3184 L -= Exp_msk1; 3185 #else /*Sudden_Underflow}{*/ 3186 #ifdef Avoid_Underflow 3187 if (bc.scale) { 3188 L = word0(&rv) & Exp_mask; 3189 if (L <= (2*P+1)*Exp_msk1) { 3190 if (L > (P+2)*Exp_msk1) 3191 /* round even ==> */ 3192 /* accept rv */ 3193 break; 3194 /* rv = smallest denormal */ 3195 if (bc.nd >nd) { 3196 bc.uflchk = 1; 3197 break; 3198 } 3199 goto undfl; 3200 } 3201 } 3202 #endif /*Avoid_Underflow*/ 3203 L = (word0(&rv) & Exp_mask) - Exp_msk1; 3204 #endif /*Sudden_Underflow}}*/ 3205 word0(&rv) = L | Bndry_mask1; 3206 word1(&rv) = 0xffffffff; 3207 #ifdef IBM 3208 goto cont; 3209 #else 3210 break; 3211 #endif 3212 } 3213 #ifndef ROUND_BIASED 3214 if (!(word1(&rv) & LSB)) 3215 break; 3216 #endif 3217 if (bc.dsign) 3218 dval(&rv) += ulp(&rv); 3219 #ifndef ROUND_BIASED 3220 else { 3221 dval(&rv) -= ulp(&rv); 3222 #ifndef Sudden_Underflow 3223 if (!dval(&rv)) { 3224 if (bc.nd >nd) { 3225 bc.uflchk = 1; 3226 break; 3227 } 3228 goto undfl; 3229 } 3230 #endif 3231 } 3232 #ifdef Avoid_Underflow 3233 bc.dsign = 1 - bc.dsign; 3234 #endif 3235 #endif 3236 break; 3237 } 3238 if ((aadj = ratio(delta, bs)) <= 2.) { 3239 if (bc.dsign) 3240 aadj = aadj1 = 1.; 3241 else if (word1(&rv) || word0(&rv) & Bndry_mask) { 3242 #ifndef Sudden_Underflow 3243 if (word1(&rv) == Tiny1 && !word0(&rv)) { 3244 if (bc.nd >nd) { 3245 bc.uflchk = 1; 3246 break; 3247 } 3248 goto undfl; 3249 } 3250 #endif 3251 aadj = 1.; 3252 aadj1 = -1.; 3253 } 3254 else { 3255 /* special case -- power of FLT_RADIX to be */ 3256 /* rounded down... */ 3257 3258 if (aadj < 2./FLT_RADIX) 3259 aadj = 1./FLT_RADIX; 3260 else 3261 aadj *= 0.5; 3262 aadj1 = -aadj; 3263 } 3264 } 3265 else { 3266 aadj *= 0.5; 3267 aadj1 = bc.dsign ? aadj : -aadj; 3268 #ifdef Check_FLT_ROUNDS 3269 switch(bc.rounding) { 3270 case 2: /* towards +infinity */ 3271 aadj1 -= 0.5; 3272 break; 3273 case 0: /* towards 0 */ 3274 case 3: /* towards -infinity */ 3275 aadj1 += 0.5; 3276 } 3277 #else 3278 if (Flt_Rounds == 0) 3279 aadj1 += 0.5; 3280 #endif /*Check_FLT_ROUNDS*/ 3281 } 3282 y = word0(&rv) & Exp_mask; 3283 3284 /* Check for overflow */ 3285 3286 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 3287 dval(&rv0) = dval(&rv); 3288 word0(&rv) -= P*Exp_msk1; 3289 adj.d = aadj1 * ulp(&rv); 3290 dval(&rv) += adj.d; 3291 if ((word0(&rv) & Exp_mask) >= 3292 Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 3293 if (word0(&rv0) == Big0 && word1(&rv0) == Big1) 3294 goto ovfl; 3295 word0(&rv) = Big0; 3296 word1(&rv) = Big1; 3297 goto cont; 3298 } 3299 else 3300 word0(&rv) += P*Exp_msk1; 3301 } 3302 else { 3303 #ifdef Avoid_Underflow 3304 if (bc.scale && y <= 2*P*Exp_msk1) { 3305 if (aadj <= 0x7fffffff) { 3306 if ((z = (ULong)aadj) <= 0) 3307 z = 1; 3308 aadj = z; 3309 aadj1 = bc.dsign ? aadj : -aadj; 3310 } 3311 dval(&aadj2) = aadj1; 3312 word0(&aadj2) += (2*P+1)*Exp_msk1 - y; 3313 aadj1 = dval(&aadj2); 3314 } 3315 adj.d = aadj1 * ulp(&rv); 3316 dval(&rv) += adj.d; 3317 #else 3318 #ifdef Sudden_Underflow 3319 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) { 3320 dval(&rv0) = dval(&rv); 3321 word0(&rv) += P*Exp_msk1; 3322 adj.d = aadj1 * ulp(&rv); 3323 dval(&rv) += adj.d; 3324 #ifdef IBM 3325 if ((word0(&rv) & Exp_mask) < P*Exp_msk1) 3326 #else 3327 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) 3328 #endif 3329 { 3330 if (word0(&rv0) == Tiny0 3331 && word1(&rv0) == Tiny1) { 3332 if (bc.nd >nd) { 3333 bc.uflchk = 1; 3334 break; 3335 } 3336 goto undfl; 3337 } 3338 word0(&rv) = Tiny0; 3339 word1(&rv) = Tiny1; 3340 goto cont; 3341 } 3342 else 3343 word0(&rv) -= P*Exp_msk1; 3344 } 3345 else { 3346 adj.d = aadj1 * ulp(&rv); 3347 dval(&rv) += adj.d; 3348 } 3349 #else /*Sudden_Underflow*/ 3350 /* Compute adj so that the IEEE rounding rules will 3351 * correctly round rv + adj in some half-way cases. 3352 * If rv * ulp(rv) is denormalized (i.e., 3353 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 3354 * trouble from bits lost to denormalization; 3355 * example: 1.2e-307 . 3356 */ 3357 if (y <= (P-1)*Exp_msk1 && aadj > 1.) { 3358 aadj1 = (double)(int)(aadj + 0.5); 3359 if (!bc.dsign) 3360 aadj1 = -aadj1; 3361 } 3362 adj.d = aadj1 * ulp(&rv); 3363 dval(&rv) += adj.d; 3364 #endif /*Sudden_Underflow*/ 3365 #endif /*Avoid_Underflow*/ 3366 } 3367 z = word0(&rv) & Exp_mask; 3368 #ifndef SET_INEXACT 3369 if (bc.nd == nd) { 3370 #ifdef Avoid_Underflow 3371 if (!bc.scale) 3372 #endif 3373 if (y == z) { 3374 /* Can we stop now? */ 3375 L = (Long)aadj; 3376 aadj -= L; 3377 /* The tolerances below are conservative. */ 3378 if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask) { 3379 if (aadj < .4999999 || aadj > .5000001) 3380 break; 3381 } 3382 else if (aadj < .4999999/FLT_RADIX) 3383 break; 3384 } 3385 } 3386 #endif 3387 cont: 3388 Bfree(bb); 3389 Bfree(bd); 3390 Bfree(bs); 3391 Bfree(delta); 3392 } 3393 Bfree(bb); 3394 Bfree(bd); 3395 Bfree(bs); 3396 Bfree(bd0); 3397 Bfree(delta); 3398 #ifndef NO_STRTOD_BIGCOMP 3399 if (bc.nd > nd) 3400 bigcomp(&rv, s0, &bc); 3401 #endif 3402 #ifdef SET_INEXACT 3403 if (bc.inexact) { 3404 if (!oldinexact) { 3405 word0(&rv0) = Exp_1 + (70 << Exp_shift); 3406 word1(&rv0) = 0; 3407 dval(&rv0) += 1.; 3408 } 3409 } 3410 else if (!oldinexact) 3411 clear_inexact(); 3412 #endif 3413 #ifdef Avoid_Underflow 3414 if (bc.scale) { 3415 word0(&rv0) = Exp_1 - 2*P*Exp_msk1; 3416 word1(&rv0) = 0; 3417 dval(&rv) *= dval(&rv0); 3418 #ifndef NO_ERRNO 3419 /* try to avoid the bug of testing an 8087 register value */ 3420 #ifdef IEEE_Arith 3421 if (!(word0(&rv) & Exp_mask)) 3422 #else 3423 if (word0(&rv) == 0 && word1(&rv) == 0) 3424 #endif 3425 errno = ERANGE; 3426 #endif 3427 } 3428 #endif /* Avoid_Underflow */ 3429 #ifdef SET_INEXACT 3430 if (bc.inexact && !(word0(&rv) & Exp_mask)) { 3431 /* set underflow bit */ 3432 dval(&rv0) = 1e-300; 3433 dval(&rv0) *= dval(&rv0); 3434 } 3435 #endif 3436 ret: 3437 if (se) 3438 *se = (char *)s; 3439 return sign ? -dval(&rv) : dval(&rv); 3440 } 3441 3442 #ifndef MULTIPLE_THREADS 3443 static char *dtoa_result; 3444 #endif 3445 3446 static char * 3447 #ifdef KR_headers 3448 rv_alloc(i) int i; 3449 #else 3450 rv_alloc(int i) 3451 #endif 3452 { 3453 int j, k, *r; 3454 3455 j = sizeof(ULong); 3456 for(k = 0; 3457 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (size_t)i; 3458 j <<= 1) 3459 k++; 3460 r = (int*)Balloc(k); 3461 *r = k; 3462 return 3463 #ifndef MULTIPLE_THREADS 3464 dtoa_result = 3465 #endif 3466 (char *)(r+1); 3467 } 3468 3469 static char * 3470 #ifdef KR_headers 3471 nrv_alloc(s, rve, n) char *s, **rve; int n; 3472 #else 3473 nrv_alloc(CONST char *s, char **rve, int n) 3474 #endif 3475 { 3476 char *rv, *t; 3477 3478 t = rv = rv_alloc(n); 3479 for(*t = *s++; *t; *t = *s++) t++; 3480 if (rve) 3481 *rve = t; 3482 return rv; 3483 } 3484 3485 /* freedtoa(s) must be used to free values s returned by dtoa 3486 * when MULTIPLE_THREADS is #defined. It should be used in all cases, 3487 * but for consistency with earlier versions of dtoa, it is optional 3488 * when MULTIPLE_THREADS is not defined. 3489 */ 3490 3491 void 3492 #ifdef KR_headers 3493 freedtoa(s) char *s; 3494 #else 3495 freedtoa(char *s) 3496 #endif 3497 { 3498 Bigint *b = (Bigint *)((int *)s - 1); 3499 b->maxwds = 1 << (b->k = *(int*)b); 3500 Bfree(b); 3501 #ifndef MULTIPLE_THREADS 3502 if (s == dtoa_result) 3503 dtoa_result = 0; 3504 #endif 3505 } 3506 3507 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. 3508 * 3509 * Inspired by "How to Print Floating-Point Numbers Accurately" by 3510 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126]. 3511 * 3512 * Modifications: 3513 * 1. Rather than iterating, we use a simple numeric overestimate 3514 * to determine k = floor(log10(d)). We scale relevant 3515 * quantities using O(log2(k)) rather than O(k) multiplications. 3516 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't 3517 * try to generate digits strictly left to right. Instead, we 3518 * compute with fewer bits and propagate the carry if necessary 3519 * when rounding the final digit up. This is often faster. 3520 * 3. Under the assumption that input will be rounded nearest, 3521 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. 3522 * That is, we allow equality in stopping tests when the 3523 * round-nearest rule will give the same floating-point value 3524 * as would satisfaction of the stopping test with strict 3525 * inequality. 3526 * 4. We remove common factors of powers of 2 from relevant 3527 * quantities. 3528 * 5. When converting floating-point integers less than 1e16, 3529 * we use floating-point arithmetic rather than resorting 3530 * to multiple-precision integers. 3531 * 6. When asked to produce fewer than 15 digits, we first try 3532 * to get by with floating-point arithmetic; we resort to 3533 * multiple-precision integer arithmetic only if we cannot 3534 * guarantee that the floating-point calculation has given 3535 * the correctly rounded result. For k requested digits and 3536 * "uniformly" distributed input, the probability is 3537 * something like 10^(k-15) that we must resort to the Long 3538 * calculation. 3539 */ 3540 3541 char * 3542 dtoa 3543 #ifdef KR_headers 3544 (dd, mode, ndigits, decpt, sign, rve) 3545 double dd; int mode, ndigits, *decpt, *sign; char **rve; 3546 #else 3547 (double dd, int mode, int ndigits, int *decpt, int *sign, char **rve) 3548 #endif 3549 { 3550 /* Arguments ndigits, decpt, sign are similar to those 3551 of ecvt and fcvt; trailing zeros are suppressed from 3552 the returned string. If not null, *rve is set to point 3553 to the end of the return value. If d is +-Infinity or NaN, 3554 then *decpt is set to 9999. 3555 3556 mode: 3557 0 ==> shortest string that yields d when read in 3558 and rounded to nearest. 3559 1 ==> like 0, but with Steele & White stopping rule; 3560 e.g. with IEEE P754 arithmetic , mode 0 gives 3561 1e23 whereas mode 1 gives 9.999999999999999e22. 3562 2 ==> max(1,ndigits) significant digits. This gives a 3563 return value similar to that of ecvt, except 3564 that trailing zeros are suppressed. 3565 3 ==> through ndigits past the decimal point. This 3566 gives a return value similar to that from fcvt, 3567 except that trailing zeros are suppressed, and 3568 ndigits can be negative. 3569 4,5 ==> similar to 2 and 3, respectively, but (in 3570 round-nearest mode) with the tests of mode 0 to 3571 possibly return a shorter string that rounds to d. 3572 With IEEE arithmetic and compilation with 3573 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same 3574 as modes 2 and 3 when FLT_ROUNDS != 1. 3575 6-9 ==> Debugging modes similar to mode - 4: don't try 3576 fast floating-point estimate (if applicable). 3577 3578 Values of mode other than 0-9 are treated as mode 0. 3579 3580 Sufficient space is allocated to the return value 3581 to hold the suppressed trailing zeros. 3582 */ 3583 3584 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1, 3585 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5, 3586 spec_case, try_quick; 3587 Long L; 3588 #ifndef Sudden_Underflow 3589 int denorm; 3590 ULong x; 3591 #endif 3592 Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S; 3593 U d2, eps, u; 3594 double ds; 3595 char *s, *s0; 3596 #ifdef SET_INEXACT 3597 int inexact, oldinexact; 3598 #endif 3599 #ifdef Honor_FLT_ROUNDS /*{*/ 3600 int Rounding; 3601 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */ 3602 Rounding = Flt_Rounds; 3603 #else /*}{*/ 3604 Rounding = 1; 3605 switch(fegetround()) { 3606 case FE_TOWARDZERO: Rounding = 0; break; 3607 case FE_UPWARD: Rounding = 2; break; 3608 case FE_DOWNWARD: Rounding = 3; 3609 } 3610 #endif /*}}*/ 3611 #endif /*}*/ 3612 3613 #ifndef MULTIPLE_THREADS 3614 if (dtoa_result) { 3615 freedtoa(dtoa_result); 3616 dtoa_result = 0; 3617 } 3618 #endif 3619 3620 u.d = dd; 3621 if (word0(&u) & Sign_bit) { 3622 /* set sign for everything, including 0's and NaNs */ 3623 *sign = 1; 3624 word0(&u) &= ~Sign_bit; /* clear sign bit */ 3625 } 3626 else 3627 *sign = 0; 3628 3629 #if defined(IEEE_Arith) + defined(VAX) 3630 #ifdef IEEE_Arith 3631 if ((word0(&u) & Exp_mask) == Exp_mask) 3632 #else 3633 if (word0(&u) == 0x8000) 3634 #endif 3635 { 3636 /* Infinity or NaN */ 3637 *decpt = 9999; 3638 #ifdef IEEE_Arith 3639 if (!word1(&u) && !(word0(&u) & 0xfffff)) 3640 return nrv_alloc("Infinity", rve, 8); 3641 #endif 3642 return nrv_alloc("NaN", rve, 3); 3643 } 3644 #endif 3645 #ifdef IBM 3646 dval(&u) += 0; /* normalize */ 3647 #endif 3648 if (!dval(&u)) { 3649 *decpt = 1; 3650 return nrv_alloc("0", rve, 1); 3651 } 3652 3653 #ifdef SET_INEXACT 3654 try_quick = oldinexact = get_inexact(); 3655 inexact = 1; 3656 #endif 3657 #ifdef Honor_FLT_ROUNDS 3658 if (Rounding >= 2) { 3659 if (*sign) 3660 Rounding = Rounding == 2 ? 0 : 2; 3661 else 3662 if (Rounding != 2) 3663 Rounding = 0; 3664 } 3665 #endif 3666 3667 b = d2b(&u, &be, &bbits); 3668 i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)); 3669 #ifndef Sudden_Underflow 3670 if (i) { 3671 #endif 3672 dval(&d2) = dval(&u); 3673 word0(&d2) &= Frac_mask1; 3674 word0(&d2) |= Exp_11; 3675 #ifdef IBM 3676 if (j = 11 - hi0bits(word0(&d2) & Frac_mask)) 3677 dval(&d2) /= 1 << j; 3678 #endif 3679 3680 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 3681 * log10(x) = log(x) / log(10) 3682 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) 3683 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) 3684 * 3685 * This suggests computing an approximation k to log10(d) by 3686 * 3687 * k = (i - Bias)*0.301029995663981 3688 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); 3689 * 3690 * We want k to be too large rather than too small. 3691 * The error in the first-order Taylor series approximation 3692 * is in our favor, so we just round up the constant enough 3693 * to compensate for any error in the multiplication of 3694 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, 3695 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, 3696 * adding 1e-13 to the constant term more than suffices. 3697 * Hence we adjust the constant term to 0.1760912590558. 3698 * (We could get a more accurate k by invoking log10, 3699 * but this is probably not worthwhile.) 3700 */ 3701 3702 i -= Bias; 3703 #ifdef IBM 3704 i <<= 2; 3705 i += j; 3706 #endif 3707 #ifndef Sudden_Underflow 3708 denorm = 0; 3709 } 3710 else { 3711 /* d is denormalized */ 3712 3713 i = bbits + be + (Bias + (P-1) - 1); 3714 x = i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32) 3715 : word1(&u) << (32 - i); 3716 dval(&d2) = x; 3717 word0(&d2) -= 31*Exp_msk1; /* adjust exponent */ 3718 i -= (Bias + (P-1) - 1) + 1; 3719 denorm = 1; 3720 } 3721 #endif 3722 ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981; 3723 k = (int)ds; 3724 if (ds < 0. && ds != k) 3725 k--; /* want k = floor(ds) */ 3726 k_check = 1; 3727 if (k >= 0 && k <= Ten_pmax) { 3728 if (dval(&u) < tens[k]) 3729 k--; 3730 k_check = 0; 3731 } 3732 j = bbits - i - 1; 3733 if (j >= 0) { 3734 b2 = 0; 3735 s2 = j; 3736 } 3737 else { 3738 b2 = -j; 3739 s2 = 0; 3740 } 3741 if (k >= 0) { 3742 b5 = 0; 3743 s5 = k; 3744 s2 += k; 3745 } 3746 else { 3747 b2 -= k; 3748 b5 = -k; 3749 s5 = 0; 3750 } 3751 if (mode < 0 || mode > 9) 3752 mode = 0; 3753 3754 #ifndef SET_INEXACT 3755 #ifdef Check_FLT_ROUNDS 3756 try_quick = Rounding == 1; 3757 #else 3758 try_quick = 1; 3759 #endif 3760 #endif /*SET_INEXACT*/ 3761 3762 if (mode > 5) { 3763 mode -= 4; 3764 try_quick = 0; 3765 } 3766 leftright = 1; 3767 ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */ 3768 /* silence erroneous "gcc -Wall" warning. */ 3769 switch(mode) { 3770 case 0: 3771 case 1: 3772 i = 18; 3773 ndigits = 0; 3774 break; 3775 case 2: 3776 leftright = 0; 3777 /* no break */ 3778 case 4: 3779 if (ndigits <= 0) 3780 ndigits = 1; 3781 ilim = ilim1 = i = ndigits; 3782 break; 3783 case 3: 3784 leftright = 0; 3785 /* no break */ 3786 case 5: 3787 i = ndigits + k + 1; 3788 ilim = i; 3789 ilim1 = i - 1; 3790 if (i <= 0) 3791 i = 1; 3792 } 3793 s = s0 = rv_alloc(i); 3794 3795 #ifdef Honor_FLT_ROUNDS 3796 if (mode > 1 && Rounding != 1) 3797 leftright = 0; 3798 #endif 3799 3800 if (ilim >= 0 && ilim <= Quick_max && try_quick) { 3801 3802 /* Try to get by with floating-point arithmetic. */ 3803 3804 i = 0; 3805 dval(&d2) = dval(&u); 3806 k0 = k; 3807 ilim0 = ilim; 3808 ieps = 2; /* conservative */ 3809 if (k > 0) { 3810 ds = tens[k&0xf]; 3811 j = k >> 4; 3812 if (j & Bletch) { 3813 /* prevent overflows */ 3814 j &= Bletch - 1; 3815 dval(&u) /= bigtens[n_bigtens-1]; 3816 ieps++; 3817 } 3818 for(; j; j >>= 1, i++) 3819 if (j & 1) { 3820 ieps++; 3821 ds *= bigtens[i]; 3822 } 3823 dval(&u) /= ds; 3824 } 3825 else { 3826 j1 = -k; 3827 if (j1) { 3828 dval(&u) *= tens[j1 & 0xf]; 3829 for(j = j1 >> 4; j; j >>= 1, i++) 3830 if (j & 1) { 3831 ieps++; 3832 dval(&u) *= bigtens[i]; 3833 } 3834 } 3835 } 3836 if (k_check && dval(&u) < 1. && ilim > 0) { 3837 if (ilim1 <= 0) 3838 goto fast_failed; 3839 ilim = ilim1; 3840 k--; 3841 dval(&u) *= 10.; 3842 ieps++; 3843 } 3844 dval(&eps) = ieps*dval(&u) + 7.; 3845 word0(&eps) -= (P-1)*Exp_msk1; 3846 if (ilim == 0) { 3847 S = mhi = 0; 3848 dval(&u) -= 5.; 3849 if (dval(&u) > dval(&eps)) 3850 goto one_digit; 3851 if (dval(&u) < -dval(&eps)) 3852 goto no_digits; 3853 goto fast_failed; 3854 } 3855 #ifndef No_leftright 3856 if (leftright) { 3857 /* Use Steele & White method of only 3858 * generating digits needed. 3859 */ 3860 dval(&eps) = 0.5/tens[ilim-1] - dval(&eps); 3861 for(i = 0;;) { 3862 L = (long)dval(&u); 3863 dval(&u) -= L; 3864 *s++ = '0' + (char)L; 3865 if (dval(&u) < dval(&eps)) 3866 goto ret1; 3867 if (1. - dval(&u) < dval(&eps)) 3868 goto bump_up; 3869 if (++i >= ilim) 3870 break; 3871 dval(&eps) *= 10.; 3872 dval(&u) *= 10.; 3873 } 3874 } 3875 else { 3876 #endif 3877 /* Generate ilim digits, then fix them up. */ 3878 dval(&eps) *= tens[ilim-1]; 3879 for(i = 1;; i++, dval(&u) *= 10.) { 3880 L = (Long)(dval(&u)); 3881 if (!(dval(&u) -= L)) 3882 ilim = i; 3883 *s++ = '0' + (char)L; 3884 if (i == ilim) { 3885 if (dval(&u) > 0.5 + dval(&eps)) 3886 goto bump_up; 3887 else if (dval(&u) < 0.5 - dval(&eps)) { 3888 while(*--s == '0') {} 3889 s++; 3890 goto ret1; 3891 } 3892 break; 3893 } 3894 } 3895 #ifndef No_leftright 3896 } 3897 #endif 3898 fast_failed: 3899 s = s0; 3900 dval(&u) = dval(&d2); 3901 k = k0; 3902 ilim = ilim0; 3903 } 3904 3905 /* Do we have a "small" integer? */ 3906 3907 if (be >= 0 && k <= Int_max) { 3908 /* Yes. */ 3909 ds = tens[k]; 3910 if (ndigits < 0 && ilim <= 0) { 3911 S = mhi = 0; 3912 if (ilim < 0 || dval(&u) <= 5*ds) 3913 goto no_digits; 3914 goto one_digit; 3915 } 3916 for(i = 1; i <= k + 1; i++, dval(&u) *= 10.) { 3917 L = (Long)(dval(&u) / ds); 3918 dval(&u) -= L*ds; 3919 #ifdef Check_FLT_ROUNDS 3920 /* If FLT_ROUNDS == 2, L will usually be high by 1 */ 3921 if (dval(&u) < 0) { 3922 L--; 3923 dval(&u) += ds; 3924 } 3925 #endif 3926 *s++ = '0' + (char)L; 3927 if (!dval(&u)) { 3928 #ifdef SET_INEXACT 3929 inexact = 0; 3930 #endif 3931 break; 3932 } 3933 if (i == ilim) { 3934 #ifdef Honor_FLT_ROUNDS 3935 if (mode > 1) 3936 switch(Rounding) { 3937 case 0: goto ret1; 3938 case 2: goto bump_up; 3939 } 3940 #endif 3941 dval(&u) += dval(&u); 3942 if (dval(&u) > ds || (dval(&u) == ds && L & 1)) { 3943 bump_up: 3944 while(*--s == '9') 3945 if (s == s0) { 3946 k++; 3947 *s = '0'; 3948 break; 3949 } 3950 ++*s++; 3951 } 3952 break; 3953 } 3954 } 3955 goto ret1; 3956 } 3957 3958 m2 = b2; 3959 m5 = b5; 3960 mhi = mlo = 0; 3961 if (leftright) { 3962 i = 3963 #ifndef Sudden_Underflow 3964 denorm ? be + (Bias + (P-1) - 1 + 1) : 3965 #endif 3966 #ifdef IBM 3967 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3); 3968 #else 3969 1 + P - bbits; 3970 #endif 3971 b2 += i; 3972 s2 += i; 3973 mhi = i2b(1); 3974 } 3975 if (m2 > 0 && s2 > 0) { 3976 i = m2 < s2 ? m2 : s2; 3977 b2 -= i; 3978 m2 -= i; 3979 s2 -= i; 3980 } 3981 if (b5 > 0) { 3982 if (leftright) { 3983 if (m5 > 0) { 3984 mhi = pow5mult(mhi, m5); 3985 b1 = mult(mhi, b); 3986 Bfree(b); 3987 b = b1; 3988 } 3989 j = b5 - m5; 3990 if (j) 3991 b = pow5mult(b, j); 3992 } 3993 else 3994 b = pow5mult(b, b5); 3995 } 3996 S = i2b(1); 3997 if (s5 > 0) 3998 S = pow5mult(S, s5); 3999 4000 /* Check for special case that d is a normalized power of 2. */ 4001 4002 spec_case = 0; 4003 if ((mode < 2 || leftright) 4004 #ifdef Honor_FLT_ROUNDS 4005 && Rounding == 1 4006 #endif 4007 ) { 4008 if (!word1(&u) && !(word0(&u) & Bndry_mask) 4009 #ifndef Sudden_Underflow 4010 && word0(&u) & (Exp_mask & ~Exp_msk1) 4011 #endif 4012 ) { 4013 /* The special case */ 4014 b2 += Log2P; 4015 s2 += Log2P; 4016 spec_case = 1; 4017 } 4018 } 4019 4020 /* Arrange for convenient computation of quotients: 4021 * shift left if necessary so divisor has 4 leading 0 bits. 4022 * 4023 * Perhaps we should just compute leading 28 bits of S once 4024 * and for all and pass them and a shift to quorem, so it 4025 * can do shifts and ors to compute the numerator for q. 4026 */ 4027 #ifdef Pack_32 4028 i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f; 4029 if (i) 4030 i = 32 - i; 4031 #define iInc 28 4032 #else 4033 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) 4034 i = 16 - i; 4035 #define iInc 12 4036 #endif 4037 i = dshift(S, s2); 4038 b2 += i; 4039 m2 += i; 4040 s2 += i; 4041 if (b2 > 0) 4042 b = lshift(b, b2); 4043 if (s2 > 0) 4044 S = lshift(S, s2); 4045 if (k_check) { 4046 if (cmp(b,S) < 0) { 4047 k--; 4048 b = multadd(b, 10, 0); /* we botched the k estimate */ 4049 if (leftright) 4050 mhi = multadd(mhi, 10, 0); 4051 ilim = ilim1; 4052 } 4053 } 4054 if (ilim <= 0 && (mode == 3 || mode == 5)) { 4055 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) { 4056 /* no digits, fcvt style */ 4057 no_digits: 4058 k = -1 - ndigits; 4059 goto ret; 4060 } 4061 one_digit: 4062 *s++ = '1'; 4063 k++; 4064 goto ret; 4065 } 4066 if (leftright) { 4067 if (m2 > 0) 4068 mhi = lshift(mhi, m2); 4069 4070 /* Compute mlo -- check for special case 4071 * that d is a normalized power of 2. 4072 */ 4073 4074 mlo = mhi; 4075 if (spec_case) { 4076 mhi = Balloc(mhi->k); 4077 Bcopy(mhi, mlo); 4078 mhi = lshift(mhi, Log2P); 4079 } 4080 4081 for(i = 1;;i++) { 4082 dig = quorem(b,S) + '0'; 4083 /* Do we yet have the shortest decimal string 4084 * that will round to d? 4085 */ 4086 j = cmp(b, mlo); 4087 delta = diff(S, mhi); 4088 j1 = delta->sign ? 1 : cmp(b, delta); 4089 Bfree(delta); 4090 #ifndef ROUND_BIASED 4091 if (j1 == 0 && mode != 1 && !(word1(&u) & 1) 4092 #ifdef Honor_FLT_ROUNDS 4093 && Rounding >= 1 4094 #endif 4095 ) { 4096 if (dig == '9') 4097 goto round_9_up; 4098 if (j > 0) 4099 dig++; 4100 #ifdef SET_INEXACT 4101 else if (!b->x[0] && b->wds <= 1) 4102 inexact = 0; 4103 #endif 4104 *s++ = (char)dig; 4105 goto ret; 4106 } 4107 #endif 4108 if (j < 0 || (j == 0 && mode != 1 4109 #ifndef ROUND_BIASED 4110 && !(word1(&u) & 1) 4111 #endif 4112 )) { 4113 if (!b->x[0] && b->wds <= 1) { 4114 #ifdef SET_INEXACT 4115 inexact = 0; 4116 #endif 4117 goto accept_dig; 4118 } 4119 #ifdef Honor_FLT_ROUNDS 4120 if (mode > 1) 4121 switch(Rounding) { 4122 case 0: goto accept_dig; 4123 case 2: goto keep_dig; 4124 } 4125 #endif /*Honor_FLT_ROUNDS*/ 4126 if (j1 > 0) { 4127 b = lshift(b, 1); 4128 j1 = cmp(b, S); 4129 if ((j1 > 0 || (j1 == 0 && dig & 1)) 4130 && dig++ == '9') 4131 goto round_9_up; 4132 } 4133 accept_dig: 4134 *s++ = (char)dig; 4135 goto ret; 4136 } 4137 if (j1 > 0) { 4138 #ifdef Honor_FLT_ROUNDS 4139 if (!Rounding) 4140 goto accept_dig; 4141 #endif 4142 if (dig == '9') { /* possible if i == 1 */ 4143 round_9_up: 4144 *s++ = '9'; 4145 goto roundoff; 4146 } 4147 *s++ = (char)dig + 1; 4148 goto ret; 4149 } 4150 #ifdef Honor_FLT_ROUNDS 4151 keep_dig: 4152 #endif 4153 *s++ = (char)dig; 4154 if (i == ilim) 4155 break; 4156 b = multadd(b, 10, 0); 4157 if (mlo == mhi) 4158 mlo = mhi = multadd(mhi, 10, 0); 4159 else { 4160 mlo = multadd(mlo, 10, 0); 4161 mhi = multadd(mhi, 10, 0); 4162 } 4163 } 4164 } 4165 else 4166 for(i = 1;; i++) { 4167 dig = quorem(b,S) + '0'; 4168 *s++ = (char)dig; 4169 if (!b->x[0] && b->wds <= 1) { 4170 #ifdef SET_INEXACT 4171 inexact = 0; 4172 #endif 4173 goto ret; 4174 } 4175 if (i >= ilim) 4176 break; 4177 b = multadd(b, 10, 0); 4178 } 4179 4180 /* Round off last digit */ 4181 4182 #ifdef Honor_FLT_ROUNDS 4183 switch(Rounding) { 4184 case 0: goto trimzeros; 4185 case 2: goto roundoff; 4186 } 4187 #endif 4188 b = lshift(b, 1); 4189 j = cmp(b, S); 4190 if (j > 0 || (j == 0 && dig & 1)) { 4191 roundoff: 4192 while(*--s == '9') 4193 if (s == s0) { 4194 k++; 4195 *s++ = '1'; 4196 goto ret; 4197 } 4198 ++*s++; 4199 } 4200 else { 4201 #ifdef Honor_FLT_ROUNDS 4202 trimzeros: 4203 #endif 4204 while(*--s == '0') {} 4205 s++; 4206 } 4207 ret: 4208 Bfree(S); 4209 if (mhi) { 4210 if (mlo && mlo != mhi) 4211 Bfree(mlo); 4212 Bfree(mhi); 4213 } 4214 ret1: 4215 #ifdef SET_INEXACT 4216 if (inexact) { 4217 if (!oldinexact) { 4218 word0(&u) = Exp_1 + (70 << Exp_shift); 4219 word1(&u) = 0; 4220 dval(&u) += 1.; 4221 } 4222 } 4223 else if (!oldinexact) 4224 clear_inexact(); 4225 #endif 4226 Bfree(b); 4227 *s = 0; 4228 *decpt = k + 1; 4229 if (rve) 4230 *rve = s; 4231 return s0; 4232 } 4233 4234 } // namespace dmg_fp 4235
