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