1 /* $NetBSD: strtod.c,v 1.45.2.1 2005/04/19 13:35:54 tron Exp $ */ 2 3 /**************************************************************** 4 * 5 * The author of this software is David M. Gay. 6 * 7 * Copyright (c) 1991 by AT&T. 8 * 9 * Permission to use, copy, modify, and distribute this software for any 10 * purpose without fee is hereby granted, provided that this entire notice 11 * is included in all copies of any software which is or includes a copy 12 * or modification of this software and in all copies of the supporting 13 * documentation for such software. 14 * 15 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED 16 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY 17 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY 18 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. 19 * 20 ***************************************************************/ 21 22 /* Please send bug reports to 23 David M. Gay 24 AT&T Bell Laboratories, Room 2C-463 25 600 Mountain Avenue 26 Murray Hill, NJ 07974-2070 27 U.S.A. 28 dmg (at) research.att.com or research!dmg 29 */ 30 31 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines. 32 * 33 * This strtod returns a nearest machine number to the input decimal 34 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are 35 * broken by the IEEE round-even rule. Otherwise ties are broken by 36 * biased rounding (add half and chop). 37 * 38 * Inspired loosely by William D. Clinger's paper "How to Read Floating 39 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. 40 * 41 * Modifications: 42 * 43 * 1. We only require IEEE, IBM, or VAX double-precision 44 * arithmetic (not IEEE double-extended). 45 * 2. We get by with floating-point arithmetic in a case that 46 * Clinger missed -- when we're computing d * 10^n 47 * for a small integer d and the integer n is not too 48 * much larger than 22 (the maximum integer k for which 49 * we can represent 10^k exactly), we may be able to 50 * compute (d*10^k) * 10^(e-k) with just one roundoff. 51 * 3. Rather than a bit-at-a-time adjustment of the binary 52 * result in the hard case, we use floating-point 53 * arithmetic to determine the adjustment to within 54 * one bit; only in really hard cases do we need to 55 * compute a second residual. 56 * 4. Because of 3., we don't need a large table of powers of 10 57 * for ten-to-e (just some small tables, e.g. of 10^k 58 * for 0 <= k <= 22). 59 */ 60 61 /* 62 * #define IEEE_LITTLE_ENDIAN for IEEE-arithmetic machines where the least 63 * significant byte has the lowest address. 64 * #define IEEE_BIG_ENDIAN for IEEE-arithmetic machines where the most 65 * significant byte has the lowest address. 66 * #define Long int on machines with 32-bit ints and 64-bit longs. 67 * #define Sudden_Underflow for IEEE-format machines without gradual 68 * underflow (i.e., that flush to zero on underflow). 69 * #define IBM for IBM mainframe-style floating-point arithmetic. 70 * #define VAX for VAX-style floating-point arithmetic. 71 * #define Unsigned_Shifts if >> does treats its left operand as unsigned. 72 * #define No_leftright to omit left-right logic in fast floating-point 73 * computation of dtoa. 74 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3. 75 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines 76 * that use extended-precision instructions to compute rounded 77 * products and quotients) with IBM. 78 * #define ROUND_BIASED for IEEE-format with biased rounding. 79 * #define Inaccurate_Divide for IEEE-format with correctly rounded 80 * products but inaccurate quotients, e.g., for Intel i860. 81 * #define Just_16 to store 16 bits per 32-bit Long when doing high-precision 82 * integer arithmetic. Whether this speeds things up or slows things 83 * down depends on the machine and the number being converted. 84 * #define KR_headers for old-style C function headers. 85 * #define Bad_float_h if your system lacks a float.h or if it does not 86 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP, 87 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX. 88 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n) 89 * if memory is available and otherwise does something you deem 90 * appropriate. If MALLOC is undefined, malloc will be invoked 91 * directly -- and assumed always to succeed. 92 */ 93 94 #ifdef ANDROID_CHANGES 95 #include <pthread.h> 96 #define mutex_lock(x) pthread_mutex_lock(x) 97 #define mutex_unlock(x) pthread_mutex_unlock(x) 98 #endif 99 100 #include <sys/cdefs.h> 101 #if defined(LIBC_SCCS) && !defined(lint) 102 __RCSID("$NetBSD: strtod.c,v 1.45.2.1 2005/04/19 13:35:54 tron Exp $"); 103 #endif /* LIBC_SCCS and not lint */ 104 105 #define Unsigned_Shifts 106 #if defined(__m68k__) || defined(__sparc__) || defined(__i386__) || \ 107 defined(__mips__) || defined(__ns32k__) || defined(__alpha__) || \ 108 defined(__powerpc__) || defined(__sh__) || defined(__x86_64__) || \ 109 defined(__hppa__) || \ 110 (defined(__arm__) && defined(__VFP_FP__)) 111 #include <endian.h> 112 #if BYTE_ORDER == BIG_ENDIAN 113 #define IEEE_BIG_ENDIAN 114 #else 115 #define IEEE_LITTLE_ENDIAN 116 #endif 117 #endif 118 119 #if defined(__arm__) && !defined(__VFP_FP__) 120 /* 121 * Although the CPU is little endian the FP has different 122 * byte and word endianness. The byte order is still little endian 123 * but the word order is big endian. 124 */ 125 #define IEEE_BIG_ENDIAN 126 #endif 127 128 #ifdef __vax__ 129 #define VAX 130 #endif 131 132 #if defined(__hppa__) || defined(__mips__) || defined(__sh__) 133 #define NAN_WORD0 0x7ff40000 134 #else 135 #define NAN_WORD0 0x7ff80000 136 #endif 137 #define NAN_WORD1 0 138 139 #define Long int32_t 140 #define ULong u_int32_t 141 142 #ifdef DEBUG 143 #include "stdio.h" 144 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} 145 #endif 146 147 #ifdef __cplusplus 148 #include "malloc.h" 149 #include "memory.h" 150 #else 151 #ifndef KR_headers 152 #include "stdlib.h" 153 #include "string.h" 154 #ifndef ANDROID_CHANGES 155 #include "locale.h" 156 #endif /* ANDROID_CHANGES */ 157 #else 158 #include "malloc.h" 159 #include "memory.h" 160 #endif 161 #endif 162 #ifndef ANDROID_CHANGES 163 #include "extern.h" 164 #include "reentrant.h" 165 #endif /* ANDROID_CHANGES */ 166 167 #ifdef MALLOC 168 #ifdef KR_headers 169 extern char *MALLOC(); 170 #else 171 extern void *MALLOC(size_t); 172 #endif 173 #else 174 #define MALLOC malloc 175 #endif 176 177 #include "ctype.h" 178 #include "errno.h" 179 #include "float.h" 180 181 #ifndef __MATH_H__ 182 #include "math.h" 183 #endif 184 185 #ifdef __cplusplus 186 extern "C" { 187 #endif 188 189 #ifndef CONST 190 #ifdef KR_headers 191 #define CONST /* blank */ 192 #else 193 #define CONST const 194 #endif 195 #endif 196 197 #ifdef Unsigned_Shifts 198 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000; 199 #else 200 #define Sign_Extend(a,b) /*no-op*/ 201 #endif 202 203 #if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN) + defined(VAX) + \ 204 defined(IBM) != 1 205 Exactly one of IEEE_LITTLE_ENDIAN IEEE_BIG_ENDIAN, VAX, or 206 IBM should be defined. 207 #endif 208 209 typedef union { 210 double d; 211 ULong ul[2]; 212 } _double; 213 #define value(x) ((x).d) 214 #ifdef IEEE_LITTLE_ENDIAN 215 #define word0(x) ((x).ul[1]) 216 #define word1(x) ((x).ul[0]) 217 #else 218 #define word0(x) ((x).ul[0]) 219 #define word1(x) ((x).ul[1]) 220 #endif 221 222 /* The following definition of Storeinc is appropriate for MIPS processors. 223 * An alternative that might be better on some machines is 224 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) 225 */ 226 #if defined(IEEE_LITTLE_ENDIAN) + defined(VAX) + defined(__arm__) 227 #define Storeinc(a,b,c) \ 228 (((u_short *)(void *)a)[1] = \ 229 (u_short)b, ((u_short *)(void *)a)[0] = (u_short)c, a++) 230 #else 231 #define Storeinc(a,b,c) \ 232 (((u_short *)(void *)a)[0] = \ 233 (u_short)b, ((u_short *)(void *)a)[1] = (u_short)c, a++) 234 #endif 235 236 /* #define P DBL_MANT_DIG */ 237 /* Ten_pmax = floor(P*log(2)/log(5)) */ 238 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ 239 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ 240 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ 241 242 #if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN) 243 #define Exp_shift 20 244 #define Exp_shift1 20 245 #define Exp_msk1 0x100000 246 #define Exp_msk11 0x100000 247 #define Exp_mask 0x7ff00000 248 #define P 53 249 #define Bias 1023 250 #define IEEE_Arith 251 #define Emin (-1022) 252 #define Exp_1 0x3ff00000 253 #define Exp_11 0x3ff00000 254 #define Ebits 11 255 #define Frac_mask 0xfffff 256 #define Frac_mask1 0xfffff 257 #define Ten_pmax 22 258 #define Bletch 0x10 259 #define Bndry_mask 0xfffff 260 #define Bndry_mask1 0xfffff 261 #define LSB 1 262 #define Sign_bit 0x80000000 263 #define Log2P 1 264 #define Tiny0 0 265 #define Tiny1 1 266 #define Quick_max 14 267 #define Int_max 14 268 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */ 269 #else 270 #undef Sudden_Underflow 271 #define Sudden_Underflow 272 #ifdef IBM 273 #define Exp_shift 24 274 #define Exp_shift1 24 275 #define Exp_msk1 0x1000000 276 #define Exp_msk11 0x1000000 277 #define Exp_mask 0x7f000000 278 #define P 14 279 #define Bias 65 280 #define Exp_1 0x41000000 281 #define Exp_11 0x41000000 282 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */ 283 #define Frac_mask 0xffffff 284 #define Frac_mask1 0xffffff 285 #define Bletch 4 286 #define Ten_pmax 22 287 #define Bndry_mask 0xefffff 288 #define Bndry_mask1 0xffffff 289 #define LSB 1 290 #define Sign_bit 0x80000000 291 #define Log2P 4 292 #define Tiny0 0x100000 293 #define Tiny1 0 294 #define Quick_max 14 295 #define Int_max 15 296 #else /* VAX */ 297 #define Exp_shift 23 298 #define Exp_shift1 7 299 #define Exp_msk1 0x80 300 #define Exp_msk11 0x800000 301 #define Exp_mask 0x7f80 302 #define P 56 303 #define Bias 129 304 #define Exp_1 0x40800000 305 #define Exp_11 0x4080 306 #define Ebits 8 307 #define Frac_mask 0x7fffff 308 #define Frac_mask1 0xffff007f 309 #define Ten_pmax 24 310 #define Bletch 2 311 #define Bndry_mask 0xffff007f 312 #define Bndry_mask1 0xffff007f 313 #define LSB 0x10000 314 #define Sign_bit 0x8000 315 #define Log2P 1 316 #define Tiny0 0x80 317 #define Tiny1 0 318 #define Quick_max 15 319 #define Int_max 15 320 #endif 321 #endif 322 323 #ifndef IEEE_Arith 324 #define ROUND_BIASED 325 #endif 326 327 #ifdef RND_PRODQUOT 328 #define rounded_product(a,b) a = rnd_prod(a, b) 329 #define rounded_quotient(a,b) a = rnd_quot(a, b) 330 #ifdef KR_headers 331 extern double rnd_prod(), rnd_quot(); 332 #else 333 extern double rnd_prod(double, double), rnd_quot(double, double); 334 #endif 335 #else 336 #define rounded_product(a,b) a *= b 337 #define rounded_quotient(a,b) a /= b 338 #endif 339 340 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) 341 #define Big1 0xffffffff 342 343 #ifndef Just_16 344 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long. 345 * This makes some inner loops simpler and sometimes saves work 346 * during multiplications, but it often seems to make things slightly 347 * slower. Hence the default is now to store 32 bits per Long. 348 */ 349 #ifndef Pack_32 350 #define Pack_32 351 #endif 352 #endif 353 354 #define Kmax 15 355 356 #ifdef __cplusplus 357 extern "C" double strtod(const char *s00, char **se); 358 extern "C" char *__dtoa(double d, int mode, int ndigits, 359 int *decpt, int *sign, char **rve); 360 #endif 361 362 struct 363 Bigint { 364 struct Bigint *next; 365 int k, maxwds, sign, wds; 366 ULong x[1]; 367 }; 368 369 typedef struct Bigint Bigint; 370 371 static Bigint *freelist[Kmax+1]; 372 373 #ifdef ANDROID_CHANGES 374 static pthread_mutex_t freelist_mutex = PTHREAD_MUTEX_INITIALIZER; 375 #else 376 #ifdef _REENTRANT 377 static mutex_t freelist_mutex = MUTEX_INITIALIZER; 378 #endif 379 #endif 380 381 /* Special value used to indicate an invalid Bigint value, 382 * e.g. when a memory allocation fails. The idea is that we 383 * want to avoid introducing NULL checks everytime a bigint 384 * computation is performed. Also the NULL value can also be 385 * already used to indicate "value not initialized yet" and 386 * returning NULL might alter the execution code path in 387 * case of OOM. 388 */ 389 #define BIGINT_INVALID ((Bigint *)&bigint_invalid_value) 390 391 static const Bigint bigint_invalid_value; 392 393 394 /* Return BIGINT_INVALID on allocation failure. 395 * 396 * Most of the code here depends on the fact that this function 397 * never returns NULL. 398 */ 399 static Bigint * 400 Balloc 401 #ifdef KR_headers 402 (k) int k; 403 #else 404 (int k) 405 #endif 406 { 407 int x; 408 Bigint *rv; 409 410 mutex_lock(&freelist_mutex); 411 412 if ((rv = freelist[k]) != NULL) { 413 freelist[k] = rv->next; 414 } 415 else { 416 x = 1 << k; 417 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(Long)); 418 if (rv == NULL) { 419 rv = BIGINT_INVALID; 420 goto EXIT; 421 } 422 rv->k = k; 423 rv->maxwds = x; 424 } 425 rv->sign = rv->wds = 0; 426 EXIT: 427 mutex_unlock(&freelist_mutex); 428 429 return rv; 430 } 431 432 static void 433 Bfree 434 #ifdef KR_headers 435 (v) Bigint *v; 436 #else 437 (Bigint *v) 438 #endif 439 { 440 if (v && v != BIGINT_INVALID) { 441 mutex_lock(&freelist_mutex); 442 443 v->next = freelist[v->k]; 444 freelist[v->k] = v; 445 446 mutex_unlock(&freelist_mutex); 447 } 448 } 449 450 #define Bcopy_valid(x,y) memcpy(&(x)->sign, &(y)->sign, \ 451 (y)->wds*sizeof(Long) + 2*sizeof(int)) 452 453 #define Bcopy(x,y) Bcopy_ptr(&(x),(y)) 454 455 static void 456 Bcopy_ptr(Bigint **px, Bigint *y) 457 { 458 if (*px == BIGINT_INVALID) 459 return; /* no space to store copy */ 460 if (y == BIGINT_INVALID) { 461 Bfree(*px); /* invalid input */ 462 *px = BIGINT_INVALID; 463 } else { 464 Bcopy_valid(*px,y); 465 } 466 } 467 468 static Bigint * 469 multadd 470 #ifdef KR_headers 471 (b, m, a) Bigint *b; int m, a; 472 #else 473 (Bigint *b, int m, int a) /* multiply by m and add a */ 474 #endif 475 { 476 int i, wds; 477 ULong *x, y; 478 #ifdef Pack_32 479 ULong xi, z; 480 #endif 481 Bigint *b1; 482 483 if (b == BIGINT_INVALID) 484 return b; 485 486 wds = b->wds; 487 x = b->x; 488 i = 0; 489 do { 490 #ifdef Pack_32 491 xi = *x; 492 y = (xi & 0xffff) * m + a; 493 z = (xi >> 16) * m + (y >> 16); 494 a = (int)(z >> 16); 495 *x++ = (z << 16) + (y & 0xffff); 496 #else 497 y = *x * m + a; 498 a = (int)(y >> 16); 499 *x++ = y & 0xffff; 500 #endif 501 } 502 while(++i < wds); 503 if (a) { 504 if (wds >= b->maxwds) { 505 b1 = Balloc(b->k+1); 506 if (b1 == BIGINT_INVALID) { 507 Bfree(b); 508 return b1; 509 } 510 Bcopy_valid(b1, b); 511 Bfree(b); 512 b = b1; 513 } 514 b->x[wds++] = a; 515 b->wds = wds; 516 } 517 return b; 518 } 519 520 static Bigint * 521 s2b 522 #ifdef KR_headers 523 (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9; 524 #else 525 (CONST char *s, int nd0, int nd, ULong y9) 526 #endif 527 { 528 Bigint *b; 529 int i, k; 530 Long x, y; 531 532 x = (nd + 8) / 9; 533 for(k = 0, y = 1; x > y; y <<= 1, k++) ; 534 #ifdef Pack_32 535 b = Balloc(k); 536 if (b == BIGINT_INVALID) 537 return b; 538 b->x[0] = y9; 539 b->wds = 1; 540 #else 541 b = Balloc(k+1); 542 if (b == BIGINT_INVALID) 543 return b; 544 545 b->x[0] = y9 & 0xffff; 546 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1; 547 #endif 548 549 i = 9; 550 if (9 < nd0) { 551 s += 9; 552 do b = multadd(b, 10, *s++ - '0'); 553 while(++i < nd0); 554 s++; 555 } 556 else 557 s += 10; 558 for(; i < nd; i++) 559 b = multadd(b, 10, *s++ - '0'); 560 return b; 561 } 562 563 static int 564 hi0bits 565 #ifdef KR_headers 566 (x) ULong x; 567 #else 568 (ULong x) 569 #endif 570 { 571 int k = 0; 572 573 if (!(x & 0xffff0000)) { 574 k = 16; 575 x <<= 16; 576 } 577 if (!(x & 0xff000000)) { 578 k += 8; 579 x <<= 8; 580 } 581 if (!(x & 0xf0000000)) { 582 k += 4; 583 x <<= 4; 584 } 585 if (!(x & 0xc0000000)) { 586 k += 2; 587 x <<= 2; 588 } 589 if (!(x & 0x80000000)) { 590 k++; 591 if (!(x & 0x40000000)) 592 return 32; 593 } 594 return k; 595 } 596 597 static int 598 lo0bits 599 #ifdef KR_headers 600 (y) ULong *y; 601 #else 602 (ULong *y) 603 #endif 604 { 605 int k; 606 ULong x = *y; 607 608 if (x & 7) { 609 if (x & 1) 610 return 0; 611 if (x & 2) { 612 *y = x >> 1; 613 return 1; 614 } 615 *y = x >> 2; 616 return 2; 617 } 618 k = 0; 619 if (!(x & 0xffff)) { 620 k = 16; 621 x >>= 16; 622 } 623 if (!(x & 0xff)) { 624 k += 8; 625 x >>= 8; 626 } 627 if (!(x & 0xf)) { 628 k += 4; 629 x >>= 4; 630 } 631 if (!(x & 0x3)) { 632 k += 2; 633 x >>= 2; 634 } 635 if (!(x & 1)) { 636 k++; 637 x >>= 1; 638 if (!x & 1) 639 return 32; 640 } 641 *y = x; 642 return k; 643 } 644 645 static Bigint * 646 i2b 647 #ifdef KR_headers 648 (i) int i; 649 #else 650 (int i) 651 #endif 652 { 653 Bigint *b; 654 655 b = Balloc(1); 656 if (b != BIGINT_INVALID) { 657 b->x[0] = i; 658 b->wds = 1; 659 } 660 return b; 661 } 662 663 static Bigint * 664 mult 665 #ifdef KR_headers 666 (a, b) Bigint *a, *b; 667 #else 668 (Bigint *a, Bigint *b) 669 #endif 670 { 671 Bigint *c; 672 int k, wa, wb, wc; 673 ULong carry, y, z; 674 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; 675 #ifdef Pack_32 676 ULong z2; 677 #endif 678 679 if (a == BIGINT_INVALID || b == BIGINT_INVALID) 680 return BIGINT_INVALID; 681 682 if (a->wds < b->wds) { 683 c = a; 684 a = b; 685 b = c; 686 } 687 k = a->k; 688 wa = a->wds; 689 wb = b->wds; 690 wc = wa + wb; 691 if (wc > a->maxwds) 692 k++; 693 c = Balloc(k); 694 if (c == BIGINT_INVALID) 695 return c; 696 for(x = c->x, xa = x + wc; x < xa; x++) 697 *x = 0; 698 xa = a->x; 699 xae = xa + wa; 700 xb = b->x; 701 xbe = xb + wb; 702 xc0 = c->x; 703 #ifdef Pack_32 704 for(; xb < xbe; xb++, xc0++) { 705 if ((y = *xb & 0xffff) != 0) { 706 x = xa; 707 xc = xc0; 708 carry = 0; 709 do { 710 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 711 carry = z >> 16; 712 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 713 carry = z2 >> 16; 714 Storeinc(xc, z2, z); 715 } 716 while(x < xae); 717 *xc = carry; 718 } 719 if ((y = *xb >> 16) != 0) { 720 x = xa; 721 xc = xc0; 722 carry = 0; 723 z2 = *xc; 724 do { 725 z = (*x & 0xffff) * y + (*xc >> 16) + carry; 726 carry = z >> 16; 727 Storeinc(xc, z, z2); 728 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 729 carry = z2 >> 16; 730 } 731 while(x < xae); 732 *xc = z2; 733 } 734 } 735 #else 736 for(; xb < xbe; xc0++) { 737 if (y = *xb++) { 738 x = xa; 739 xc = xc0; 740 carry = 0; 741 do { 742 z = *x++ * y + *xc + carry; 743 carry = z >> 16; 744 *xc++ = z & 0xffff; 745 } 746 while(x < xae); 747 *xc = carry; 748 } 749 } 750 #endif 751 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ; 752 c->wds = wc; 753 return c; 754 } 755 756 static Bigint *p5s; 757 static pthread_mutex_t p5s_mutex = PTHREAD_MUTEX_INITIALIZER; 758 759 static Bigint * 760 pow5mult 761 #ifdef KR_headers 762 (b, k) Bigint *b; int k; 763 #else 764 (Bigint *b, int k) 765 #endif 766 { 767 Bigint *b1, *p5, *p51; 768 int i; 769 static const int p05[3] = { 5, 25, 125 }; 770 771 if (b == BIGINT_INVALID) 772 return b; 773 774 if ((i = k & 3) != 0) 775 b = multadd(b, p05[i-1], 0); 776 777 if (!(k = (unsigned int) k >> 2)) 778 return b; 779 mutex_lock(&p5s_mutex); 780 if (!(p5 = p5s)) { 781 /* first time */ 782 p5 = i2b(625); 783 if (p5 == BIGINT_INVALID) { 784 Bfree(b); 785 mutex_unlock(&p5s_mutex); 786 return p5; 787 } 788 p5s = p5; 789 p5->next = 0; 790 } 791 for(;;) { 792 if (k & 1) { 793 b1 = mult(b, p5); 794 Bfree(b); 795 b = b1; 796 } 797 if (!(k = (unsigned int) k >> 1)) 798 break; 799 if (!(p51 = p5->next)) { 800 p51 = mult(p5,p5); 801 if (p51 == BIGINT_INVALID) { 802 Bfree(b); 803 mutex_unlock(&p5s_mutex); 804 return p51; 805 } 806 p5->next = p51; 807 p51->next = 0; 808 } 809 p5 = p51; 810 } 811 mutex_unlock(&p5s_mutex); 812 return b; 813 } 814 815 static Bigint * 816 lshift 817 #ifdef KR_headers 818 (b, k) Bigint *b; int k; 819 #else 820 (Bigint *b, int k) 821 #endif 822 { 823 int i, k1, n, n1; 824 Bigint *b1; 825 ULong *x, *x1, *xe, z; 826 827 if (b == BIGINT_INVALID) 828 return b; 829 830 #ifdef Pack_32 831 n = (unsigned int)k >> 5; 832 #else 833 n = (unsigned int)k >> 4; 834 #endif 835 k1 = b->k; 836 n1 = n + b->wds + 1; 837 for(i = b->maxwds; n1 > i; i <<= 1) 838 k1++; 839 b1 = Balloc(k1); 840 if (b1 == BIGINT_INVALID) { 841 Bfree(b); 842 return b1; 843 } 844 x1 = b1->x; 845 for(i = 0; i < n; i++) 846 *x1++ = 0; 847 x = b->x; 848 xe = x + b->wds; 849 #ifdef Pack_32 850 if (k &= 0x1f) { 851 k1 = 32 - k; 852 z = 0; 853 do { 854 *x1++ = *x << k | z; 855 z = *x++ >> k1; 856 } 857 while(x < xe); 858 if ((*x1 = z) != 0) 859 ++n1; 860 } 861 #else 862 if (k &= 0xf) { 863 k1 = 16 - k; 864 z = 0; 865 do { 866 *x1++ = *x << k & 0xffff | z; 867 z = *x++ >> k1; 868 } 869 while(x < xe); 870 if (*x1 = z) 871 ++n1; 872 } 873 #endif 874 else do 875 *x1++ = *x++; 876 while(x < xe); 877 b1->wds = n1 - 1; 878 Bfree(b); 879 return b1; 880 } 881 882 static int 883 cmp 884 #ifdef KR_headers 885 (a, b) Bigint *a, *b; 886 #else 887 (Bigint *a, Bigint *b) 888 #endif 889 { 890 ULong *xa, *xa0, *xb, *xb0; 891 int i, j; 892 893 if (a == BIGINT_INVALID || b == BIGINT_INVALID) 894 #ifdef DEBUG 895 Bug("cmp called with a or b invalid"); 896 #else 897 return 0; /* equal - the best we can do right now */ 898 #endif 899 900 i = a->wds; 901 j = b->wds; 902 #ifdef DEBUG 903 if (i > 1 && !a->x[i-1]) 904 Bug("cmp called with a->x[a->wds-1] == 0"); 905 if (j > 1 && !b->x[j-1]) 906 Bug("cmp called with b->x[b->wds-1] == 0"); 907 #endif 908 if (i -= j) 909 return i; 910 xa0 = a->x; 911 xa = xa0 + j; 912 xb0 = b->x; 913 xb = xb0 + j; 914 for(;;) { 915 if (*--xa != *--xb) 916 return *xa < *xb ? -1 : 1; 917 if (xa <= xa0) 918 break; 919 } 920 return 0; 921 } 922 923 static Bigint * 924 diff 925 #ifdef KR_headers 926 (a, b) Bigint *a, *b; 927 #else 928 (Bigint *a, Bigint *b) 929 #endif 930 { 931 Bigint *c; 932 int i, wa, wb; 933 Long borrow, y; /* We need signed shifts here. */ 934 ULong *xa, *xae, *xb, *xbe, *xc; 935 #ifdef Pack_32 936 Long z; 937 #endif 938 939 if (a == BIGINT_INVALID || b == BIGINT_INVALID) 940 return BIGINT_INVALID; 941 942 i = cmp(a,b); 943 if (!i) { 944 c = Balloc(0); 945 if (c != BIGINT_INVALID) { 946 c->wds = 1; 947 c->x[0] = 0; 948 } 949 return c; 950 } 951 if (i < 0) { 952 c = a; 953 a = b; 954 b = c; 955 i = 1; 956 } 957 else 958 i = 0; 959 c = Balloc(a->k); 960 if (c == BIGINT_INVALID) 961 return c; 962 c->sign = i; 963 wa = a->wds; 964 xa = a->x; 965 xae = xa + wa; 966 wb = b->wds; 967 xb = b->x; 968 xbe = xb + wb; 969 xc = c->x; 970 borrow = 0; 971 #ifdef Pack_32 972 do { 973 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow; 974 borrow = (ULong)y >> 16; 975 Sign_Extend(borrow, y); 976 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow; 977 borrow = (ULong)z >> 16; 978 Sign_Extend(borrow, z); 979 Storeinc(xc, z, y); 980 } 981 while(xb < xbe); 982 while(xa < xae) { 983 y = (*xa & 0xffff) + borrow; 984 borrow = (ULong)y >> 16; 985 Sign_Extend(borrow, y); 986 z = (*xa++ >> 16) + borrow; 987 borrow = (ULong)z >> 16; 988 Sign_Extend(borrow, z); 989 Storeinc(xc, z, y); 990 } 991 #else 992 do { 993 y = *xa++ - *xb++ + borrow; 994 borrow = y >> 16; 995 Sign_Extend(borrow, y); 996 *xc++ = y & 0xffff; 997 } 998 while(xb < xbe); 999 while(xa < xae) { 1000 y = *xa++ + borrow; 1001 borrow = y >> 16; 1002 Sign_Extend(borrow, y); 1003 *xc++ = y & 0xffff; 1004 } 1005 #endif 1006 while(!*--xc) 1007 wa--; 1008 c->wds = wa; 1009 return c; 1010 } 1011 1012 static double 1013 ulp 1014 #ifdef KR_headers 1015 (_x) double _x; 1016 #else 1017 (double _x) 1018 #endif 1019 { 1020 _double x; 1021 Long L; 1022 _double a; 1023 1024 value(x) = _x; 1025 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; 1026 #ifndef Sudden_Underflow 1027 if (L > 0) { 1028 #endif 1029 #ifdef IBM 1030 L |= Exp_msk1 >> 4; 1031 #endif 1032 word0(a) = L; 1033 word1(a) = 0; 1034 #ifndef Sudden_Underflow 1035 } 1036 else { 1037 L = (ULong)-L >> Exp_shift; 1038 if (L < Exp_shift) { 1039 word0(a) = 0x80000 >> L; 1040 word1(a) = 0; 1041 } 1042 else { 1043 word0(a) = 0; 1044 L -= Exp_shift; 1045 word1(a) = L >= 31 ? 1 : 1 << (31 - L); 1046 } 1047 } 1048 #endif 1049 return value(a); 1050 } 1051 1052 static double 1053 b2d 1054 #ifdef KR_headers 1055 (a, e) Bigint *a; int *e; 1056 #else 1057 (Bigint *a, int *e) 1058 #endif 1059 { 1060 ULong *xa, *xa0, w, y, z; 1061 int k; 1062 _double d; 1063 #ifdef VAX 1064 ULong d0, d1; 1065 #else 1066 #define d0 word0(d) 1067 #define d1 word1(d) 1068 #endif 1069 1070 if (a == BIGINT_INVALID) 1071 return NAN; 1072 1073 xa0 = a->x; 1074 xa = xa0 + a->wds; 1075 y = *--xa; 1076 #ifdef DEBUG 1077 if (!y) Bug("zero y in b2d"); 1078 #endif 1079 k = hi0bits(y); 1080 *e = 32 - k; 1081 #ifdef Pack_32 1082 if (k < Ebits) { 1083 d0 = Exp_1 | y >> (Ebits - k); 1084 w = xa > xa0 ? *--xa : 0; 1085 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k); 1086 goto ret_d; 1087 } 1088 z = xa > xa0 ? *--xa : 0; 1089 if (k -= Ebits) { 1090 d0 = Exp_1 | y << k | z >> (32 - k); 1091 y = xa > xa0 ? *--xa : 0; 1092 d1 = z << k | y >> (32 - k); 1093 } 1094 else { 1095 d0 = Exp_1 | y; 1096 d1 = z; 1097 } 1098 #else 1099 if (k < Ebits + 16) { 1100 z = xa > xa0 ? *--xa : 0; 1101 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; 1102 w = xa > xa0 ? *--xa : 0; 1103 y = xa > xa0 ? *--xa : 0; 1104 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; 1105 goto ret_d; 1106 } 1107 z = xa > xa0 ? *--xa : 0; 1108 w = xa > xa0 ? *--xa : 0; 1109 k -= Ebits + 16; 1110 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; 1111 y = xa > xa0 ? *--xa : 0; 1112 d1 = w << k + 16 | y << k; 1113 #endif 1114 ret_d: 1115 #ifdef VAX 1116 word0(d) = d0 >> 16 | d0 << 16; 1117 word1(d) = d1 >> 16 | d1 << 16; 1118 #else 1119 #undef d0 1120 #undef d1 1121 #endif 1122 return value(d); 1123 } 1124 1125 static Bigint * 1126 d2b 1127 #ifdef KR_headers 1128 (_d, e, bits) double d; int *e, *bits; 1129 #else 1130 (double _d, int *e, int *bits) 1131 #endif 1132 { 1133 Bigint *b; 1134 int de, i, k; 1135 ULong *x, y, z; 1136 _double d; 1137 #ifdef VAX 1138 ULong d0, d1; 1139 #endif 1140 1141 value(d) = _d; 1142 #ifdef VAX 1143 d0 = word0(d) >> 16 | word0(d) << 16; 1144 d1 = word1(d) >> 16 | word1(d) << 16; 1145 #else 1146 #define d0 word0(d) 1147 #define d1 word1(d) 1148 #endif 1149 1150 #ifdef Pack_32 1151 b = Balloc(1); 1152 #else 1153 b = Balloc(2); 1154 #endif 1155 if (b == BIGINT_INVALID) 1156 return b; 1157 x = b->x; 1158 1159 z = d0 & Frac_mask; 1160 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 1161 #ifdef Sudden_Underflow 1162 de = (int)(d0 >> Exp_shift); 1163 #ifndef IBM 1164 z |= Exp_msk11; 1165 #endif 1166 #else 1167 if ((de = (int)(d0 >> Exp_shift)) != 0) 1168 z |= Exp_msk1; 1169 #endif 1170 #ifdef Pack_32 1171 if ((y = d1) != 0) { 1172 if ((k = lo0bits(&y)) != 0) { 1173 x[0] = y | z << (32 - k); 1174 z >>= k; 1175 } 1176 else 1177 x[0] = y; 1178 i = b->wds = (x[1] = z) ? 2 : 1; 1179 } 1180 else { 1181 #ifdef DEBUG 1182 if (!z) 1183 Bug("Zero passed to d2b"); 1184 #endif 1185 k = lo0bits(&z); 1186 x[0] = z; 1187 i = b->wds = 1; 1188 k += 32; 1189 } 1190 #else 1191 if (y = d1) { 1192 if (k = lo0bits(&y)) 1193 if (k >= 16) { 1194 x[0] = y | z << 32 - k & 0xffff; 1195 x[1] = z >> k - 16 & 0xffff; 1196 x[2] = z >> k; 1197 i = 2; 1198 } 1199 else { 1200 x[0] = y & 0xffff; 1201 x[1] = y >> 16 | z << 16 - k & 0xffff; 1202 x[2] = z >> k & 0xffff; 1203 x[3] = z >> k+16; 1204 i = 3; 1205 } 1206 else { 1207 x[0] = y & 0xffff; 1208 x[1] = y >> 16; 1209 x[2] = z & 0xffff; 1210 x[3] = z >> 16; 1211 i = 3; 1212 } 1213 } 1214 else { 1215 #ifdef DEBUG 1216 if (!z) 1217 Bug("Zero passed to d2b"); 1218 #endif 1219 k = lo0bits(&z); 1220 if (k >= 16) { 1221 x[0] = z; 1222 i = 0; 1223 } 1224 else { 1225 x[0] = z & 0xffff; 1226 x[1] = z >> 16; 1227 i = 1; 1228 } 1229 k += 32; 1230 } 1231 while(!x[i]) 1232 --i; 1233 b->wds = i + 1; 1234 #endif 1235 #ifndef Sudden_Underflow 1236 if (de) { 1237 #endif 1238 #ifdef IBM 1239 *e = (de - Bias - (P-1) << 2) + k; 1240 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask); 1241 #else 1242 *e = de - Bias - (P-1) + k; 1243 *bits = P - k; 1244 #endif 1245 #ifndef Sudden_Underflow 1246 } 1247 else { 1248 *e = de - Bias - (P-1) + 1 + k; 1249 #ifdef Pack_32 1250 *bits = 32*i - hi0bits(x[i-1]); 1251 #else 1252 *bits = (i+2)*16 - hi0bits(x[i]); 1253 #endif 1254 } 1255 #endif 1256 return b; 1257 } 1258 #undef d0 1259 #undef d1 1260 1261 static double 1262 ratio 1263 #ifdef KR_headers 1264 (a, b) Bigint *a, *b; 1265 #else 1266 (Bigint *a, Bigint *b) 1267 #endif 1268 { 1269 _double da, db; 1270 int k, ka, kb; 1271 1272 if (a == BIGINT_INVALID || b == BIGINT_INVALID) 1273 return NAN; /* for lack of better value ? */ 1274 1275 value(da) = b2d(a, &ka); 1276 value(db) = b2d(b, &kb); 1277 #ifdef Pack_32 1278 k = ka - kb + 32*(a->wds - b->wds); 1279 #else 1280 k = ka - kb + 16*(a->wds - b->wds); 1281 #endif 1282 #ifdef IBM 1283 if (k > 0) { 1284 word0(da) += (k >> 2)*Exp_msk1; 1285 if (k &= 3) 1286 da *= 1 << k; 1287 } 1288 else { 1289 k = -k; 1290 word0(db) += (k >> 2)*Exp_msk1; 1291 if (k &= 3) 1292 db *= 1 << k; 1293 } 1294 #else 1295 if (k > 0) 1296 word0(da) += k*Exp_msk1; 1297 else { 1298 k = -k; 1299 word0(db) += k*Exp_msk1; 1300 } 1301 #endif 1302 return value(da) / value(db); 1303 } 1304 1305 static CONST double 1306 tens[] = { 1307 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1308 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1309 1e20, 1e21, 1e22 1310 #ifdef VAX 1311 , 1e23, 1e24 1312 #endif 1313 }; 1314 1315 #ifdef IEEE_Arith 1316 static CONST double bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; 1317 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 }; 1318 #define n_bigtens 5 1319 #else 1320 #ifdef IBM 1321 static CONST double bigtens[] = { 1e16, 1e32, 1e64 }; 1322 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 }; 1323 #define n_bigtens 3 1324 #else 1325 static CONST double bigtens[] = { 1e16, 1e32 }; 1326 static CONST double tinytens[] = { 1e-16, 1e-32 }; 1327 #define n_bigtens 2 1328 #endif 1329 #endif 1330 1331 double 1332 strtod 1333 #ifdef KR_headers 1334 (s00, se) CONST char *s00; char **se; 1335 #else 1336 (CONST char *s00, char **se) 1337 #endif 1338 { 1339 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, 1340 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; 1341 CONST char *s, *s0, *s1; 1342 double aadj, aadj1, adj; 1343 _double rv, rv0; 1344 Long L; 1345 ULong y, z; 1346 Bigint *bb1, *bd0; 1347 Bigint *bb = NULL, *bd = NULL, *bs = NULL, *delta = NULL;/* pacify gcc */ 1348 1349 #ifdef ANDROID_CHANGES 1350 CONST char decimal_point = '.'; 1351 #else /* ANDROID_CHANGES */ 1352 #ifndef KR_headers 1353 CONST char decimal_point = localeconv()->decimal_point[0]; 1354 #else 1355 CONST char decimal_point = '.'; 1356 #endif 1357 1358 #endif /* ANDROID_CHANGES */ 1359 1360 sign = nz0 = nz = 0; 1361 value(rv) = 0.; 1362 1363 1364 for(s = s00; isspace((unsigned char) *s); s++) 1365 ; 1366 1367 if (*s == '-') { 1368 sign = 1; 1369 s++; 1370 } else if (*s == '+') { 1371 s++; 1372 } 1373 1374 if (*s == '\0') { 1375 s = s00; 1376 goto ret; 1377 } 1378 1379 /* "INF" or "INFINITY" */ 1380 if (tolower((unsigned char)*s) == 'i' && strncasecmp(s, "inf", 3) == 0) { 1381 if (strncasecmp(s + 3, "inity", 5) == 0) 1382 s += 8; 1383 else 1384 s += 3; 1385 1386 value(rv) = HUGE_VAL; 1387 goto ret; 1388 } 1389 1390 #ifdef IEEE_Arith 1391 /* "NAN" or "NAN(n-char-sequence-opt)" */ 1392 if (tolower((unsigned char)*s) == 'n' && strncasecmp(s, "nan", 3) == 0) { 1393 /* Build a quiet NaN. */ 1394 word0(rv) = NAN_WORD0; 1395 word1(rv) = NAN_WORD1; 1396 s+= 3; 1397 1398 /* Don't interpret (n-char-sequence-opt), for now. */ 1399 if (*s == '(') { 1400 s0 = s; 1401 for (s++; *s != ')' && *s != '\0'; s++) 1402 ; 1403 if (*s == ')') 1404 s++; /* Skip over closing paren ... */ 1405 else 1406 s = s0; /* ... otherwise go back. */ 1407 } 1408 1409 goto ret; 1410 } 1411 #endif 1412 1413 if (*s == '0') { 1414 nz0 = 1; 1415 while(*++s == '0') ; 1416 if (!*s) 1417 goto ret; 1418 } 1419 s0 = s; 1420 y = z = 0; 1421 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 1422 if (nd < 9) 1423 y = 10*y + c - '0'; 1424 else if (nd < 16) 1425 z = 10*z + c - '0'; 1426 nd0 = nd; 1427 if (c == decimal_point) { 1428 c = *++s; 1429 if (!nd) { 1430 for(; c == '0'; c = *++s) 1431 nz++; 1432 if (c > '0' && c <= '9') { 1433 s0 = s; 1434 nf += nz; 1435 nz = 0; 1436 goto have_dig; 1437 } 1438 goto dig_done; 1439 } 1440 for(; c >= '0' && c <= '9'; c = *++s) { 1441 have_dig: 1442 nz++; 1443 if (c -= '0') { 1444 nf += nz; 1445 for(i = 1; i < nz; i++) 1446 if (nd++ < 9) 1447 y *= 10; 1448 else if (nd <= DBL_DIG + 1) 1449 z *= 10; 1450 if (nd++ < 9) 1451 y = 10*y + c; 1452 else if (nd <= DBL_DIG + 1) 1453 z = 10*z + c; 1454 nz = 0; 1455 } 1456 } 1457 } 1458 dig_done: 1459 e = 0; 1460 if (c == 'e' || c == 'E') { 1461 if (!nd && !nz && !nz0) { 1462 s = s00; 1463 goto ret; 1464 } 1465 s00 = s; 1466 esign = 0; 1467 switch(c = *++s) { 1468 case '-': 1469 esign = 1; 1470 /* FALLTHROUGH */ 1471 case '+': 1472 c = *++s; 1473 } 1474 if (c >= '0' && c <= '9') { 1475 while(c == '0') 1476 c = *++s; 1477 if (c > '0' && c <= '9') { 1478 L = c - '0'; 1479 s1 = s; 1480 while((c = *++s) >= '0' && c <= '9') 1481 L = 10*L + c - '0'; 1482 if (s - s1 > 8 || L > 19999) 1483 /* Avoid confusion from exponents 1484 * so large that e might overflow. 1485 */ 1486 e = 19999; /* safe for 16 bit ints */ 1487 else 1488 e = (int)L; 1489 if (esign) 1490 e = -e; 1491 } 1492 else 1493 e = 0; 1494 } 1495 else 1496 s = s00; 1497 } 1498 if (!nd) { 1499 if (!nz && !nz0) 1500 s = s00; 1501 goto ret; 1502 } 1503 e1 = e -= nf; 1504 1505 /* Now we have nd0 digits, starting at s0, followed by a 1506 * decimal point, followed by nd-nd0 digits. The number we're 1507 * after is the integer represented by those digits times 1508 * 10**e */ 1509 1510 if (!nd0) 1511 nd0 = nd; 1512 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 1513 value(rv) = y; 1514 if (k > 9) 1515 value(rv) = tens[k - 9] * value(rv) + z; 1516 bd0 = 0; 1517 if (nd <= DBL_DIG 1518 #ifndef RND_PRODQUOT 1519 && FLT_ROUNDS == 1 1520 #endif 1521 ) { 1522 if (!e) 1523 goto ret; 1524 if (e > 0) { 1525 if (e <= Ten_pmax) { 1526 #ifdef VAX 1527 goto vax_ovfl_check; 1528 #else 1529 /* value(rv) = */ rounded_product(value(rv), 1530 tens[e]); 1531 goto ret; 1532 #endif 1533 } 1534 i = DBL_DIG - nd; 1535 if (e <= Ten_pmax + i) { 1536 /* A fancier test would sometimes let us do 1537 * this for larger i values. 1538 */ 1539 e -= i; 1540 value(rv) *= tens[i]; 1541 #ifdef VAX 1542 /* VAX exponent range is so narrow we must 1543 * worry about overflow here... 1544 */ 1545 vax_ovfl_check: 1546 word0(rv) -= P*Exp_msk1; 1547 /* value(rv) = */ rounded_product(value(rv), 1548 tens[e]); 1549 if ((word0(rv) & Exp_mask) 1550 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 1551 goto ovfl; 1552 word0(rv) += P*Exp_msk1; 1553 #else 1554 /* value(rv) = */ rounded_product(value(rv), 1555 tens[e]); 1556 #endif 1557 goto ret; 1558 } 1559 } 1560 #ifndef Inaccurate_Divide 1561 else if (e >= -Ten_pmax) { 1562 /* value(rv) = */ rounded_quotient(value(rv), 1563 tens[-e]); 1564 goto ret; 1565 } 1566 #endif 1567 } 1568 e1 += nd - k; 1569 1570 /* Get starting approximation = rv * 10**e1 */ 1571 1572 if (e1 > 0) { 1573 if ((i = e1 & 15) != 0) 1574 value(rv) *= tens[i]; 1575 if (e1 &= ~15) { 1576 if (e1 > DBL_MAX_10_EXP) { 1577 ovfl: 1578 errno = ERANGE; 1579 value(rv) = HUGE_VAL; 1580 if (bd0) 1581 goto retfree; 1582 goto ret; 1583 } 1584 if ((e1 = (unsigned int)e1 >> 4) != 0) { 1585 for(j = 0; e1 > 1; j++, 1586 e1 = (unsigned int)e1 >> 1) 1587 if (e1 & 1) 1588 value(rv) *= bigtens[j]; 1589 /* The last multiplication could overflow. */ 1590 word0(rv) -= P*Exp_msk1; 1591 value(rv) *= bigtens[j]; 1592 if ((z = word0(rv) & Exp_mask) 1593 > Exp_msk1*(DBL_MAX_EXP+Bias-P)) 1594 goto ovfl; 1595 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { 1596 /* set to largest number */ 1597 /* (Can't trust DBL_MAX) */ 1598 word0(rv) = Big0; 1599 word1(rv) = Big1; 1600 } 1601 else 1602 word0(rv) += P*Exp_msk1; 1603 } 1604 } 1605 } 1606 else if (e1 < 0) { 1607 e1 = -e1; 1608 if ((i = e1 & 15) != 0) 1609 value(rv) /= tens[i]; 1610 if (e1 &= ~15) { 1611 e1 = (unsigned int)e1 >> 4; 1612 if (e1 >= 1 << n_bigtens) 1613 goto undfl; 1614 for(j = 0; e1 > 1; j++, 1615 e1 = (unsigned int)e1 >> 1) 1616 if (e1 & 1) 1617 value(rv) *= tinytens[j]; 1618 /* The last multiplication could underflow. */ 1619 value(rv0) = value(rv); 1620 value(rv) *= tinytens[j]; 1621 if (!value(rv)) { 1622 value(rv) = 2.*value(rv0); 1623 value(rv) *= tinytens[j]; 1624 if (!value(rv)) { 1625 undfl: 1626 value(rv) = 0.; 1627 errno = ERANGE; 1628 if (bd0) 1629 goto retfree; 1630 goto ret; 1631 } 1632 word0(rv) = Tiny0; 1633 word1(rv) = Tiny1; 1634 /* The refinement below will clean 1635 * this approximation up. 1636 */ 1637 } 1638 } 1639 } 1640 1641 /* Now the hard part -- adjusting rv to the correct value.*/ 1642 1643 /* Put digits into bd: true value = bd * 10^e */ 1644 1645 bd0 = s2b(s0, nd0, nd, y); 1646 1647 for(;;) { 1648 bd = Balloc(bd0->k); 1649 Bcopy(bd, bd0); 1650 bb = d2b(value(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */ 1651 bs = i2b(1); 1652 1653 if (e >= 0) { 1654 bb2 = bb5 = 0; 1655 bd2 = bd5 = e; 1656 } 1657 else { 1658 bb2 = bb5 = -e; 1659 bd2 = bd5 = 0; 1660 } 1661 if (bbe >= 0) 1662 bb2 += bbe; 1663 else 1664 bd2 -= bbe; 1665 bs2 = bb2; 1666 #ifdef Sudden_Underflow 1667 #ifdef IBM 1668 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 1669 #else 1670 j = P + 1 - bbbits; 1671 #endif 1672 #else 1673 i = bbe + bbbits - 1; /* logb(rv) */ 1674 if (i < Emin) /* denormal */ 1675 j = bbe + (P-Emin); 1676 else 1677 j = P + 1 - bbbits; 1678 #endif 1679 bb2 += j; 1680 bd2 += j; 1681 i = bb2 < bd2 ? bb2 : bd2; 1682 if (i > bs2) 1683 i = bs2; 1684 if (i > 0) { 1685 bb2 -= i; 1686 bd2 -= i; 1687 bs2 -= i; 1688 } 1689 if (bb5 > 0) { 1690 bs = pow5mult(bs, bb5); 1691 bb1 = mult(bs, bb); 1692 Bfree(bb); 1693 bb = bb1; 1694 } 1695 if (bb2 > 0) 1696 bb = lshift(bb, bb2); 1697 if (bd5 > 0) 1698 bd = pow5mult(bd, bd5); 1699 if (bd2 > 0) 1700 bd = lshift(bd, bd2); 1701 if (bs2 > 0) 1702 bs = lshift(bs, bs2); 1703 delta = diff(bb, bd); 1704 dsign = delta->sign; 1705 delta->sign = 0; 1706 i = cmp(delta, bs); 1707 if (i < 0) { 1708 /* Error is less than half an ulp -- check for 1709 * special case of mantissa a power of two. 1710 */ 1711 if (dsign || word1(rv) || word0(rv) & Bndry_mask) 1712 break; 1713 delta = lshift(delta,Log2P); 1714 if (cmp(delta, bs) > 0) 1715 goto drop_down; 1716 break; 1717 } 1718 if (i == 0) { 1719 /* exactly half-way between */ 1720 if (dsign) { 1721 if ((word0(rv) & Bndry_mask1) == Bndry_mask1 1722 && word1(rv) == 0xffffffff) { 1723 /*boundary case -- increment exponent*/ 1724 word0(rv) = (word0(rv) & Exp_mask) 1725 + Exp_msk1 1726 #ifdef IBM 1727 | Exp_msk1 >> 4 1728 #endif 1729 ; 1730 word1(rv) = 0; 1731 break; 1732 } 1733 } 1734 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { 1735 drop_down: 1736 /* boundary case -- decrement exponent */ 1737 #ifdef Sudden_Underflow 1738 L = word0(rv) & Exp_mask; 1739 #ifdef IBM 1740 if (L < Exp_msk1) 1741 #else 1742 if (L <= Exp_msk1) 1743 #endif 1744 goto undfl; 1745 L -= Exp_msk1; 1746 #else 1747 L = (word0(rv) & Exp_mask) - Exp_msk1; 1748 #endif 1749 word0(rv) = L | Bndry_mask1; 1750 word1(rv) = 0xffffffff; 1751 #ifdef IBM 1752 goto cont; 1753 #else 1754 break; 1755 #endif 1756 } 1757 #ifndef ROUND_BIASED 1758 if (!(word1(rv) & LSB)) 1759 break; 1760 #endif 1761 if (dsign) 1762 value(rv) += ulp(value(rv)); 1763 #ifndef ROUND_BIASED 1764 else { 1765 value(rv) -= ulp(value(rv)); 1766 #ifndef Sudden_Underflow 1767 if (!value(rv)) 1768 goto undfl; 1769 #endif 1770 } 1771 #endif 1772 break; 1773 } 1774 if ((aadj = ratio(delta, bs)) <= 2.) { 1775 if (dsign) 1776 aadj = aadj1 = 1.; 1777 else if (word1(rv) || word0(rv) & Bndry_mask) { 1778 #ifndef Sudden_Underflow 1779 if (word1(rv) == Tiny1 && !word0(rv)) 1780 goto undfl; 1781 #endif 1782 aadj = 1.; 1783 aadj1 = -1.; 1784 } 1785 else { 1786 /* special case -- power of FLT_RADIX to be */ 1787 /* rounded down... */ 1788 1789 if (aadj < 2./FLT_RADIX) 1790 aadj = 1./FLT_RADIX; 1791 else 1792 aadj *= 0.5; 1793 aadj1 = -aadj; 1794 } 1795 } 1796 else { 1797 aadj *= 0.5; 1798 aadj1 = dsign ? aadj : -aadj; 1799 #ifdef Check_FLT_ROUNDS 1800 switch(FLT_ROUNDS) { 1801 case 2: /* towards +infinity */ 1802 aadj1 -= 0.5; 1803 break; 1804 case 0: /* towards 0 */ 1805 case 3: /* towards -infinity */ 1806 aadj1 += 0.5; 1807 } 1808 #else 1809 if (FLT_ROUNDS == 0) 1810 aadj1 += 0.5; 1811 #endif 1812 } 1813 y = word0(rv) & Exp_mask; 1814 1815 /* Check for overflow */ 1816 1817 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 1818 value(rv0) = value(rv); 1819 word0(rv) -= P*Exp_msk1; 1820 adj = aadj1 * ulp(value(rv)); 1821 value(rv) += adj; 1822 if ((word0(rv) & Exp_mask) >= 1823 Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 1824 if (word0(rv0) == Big0 && word1(rv0) == Big1) 1825 goto ovfl; 1826 word0(rv) = Big0; 1827 word1(rv) = Big1; 1828 goto cont; 1829 } 1830 else 1831 word0(rv) += P*Exp_msk1; 1832 } 1833 else { 1834 #ifdef Sudden_Underflow 1835 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 1836 value(rv0) = value(rv); 1837 word0(rv) += P*Exp_msk1; 1838 adj = aadj1 * ulp(value(rv)); 1839 value(rv) += adj; 1840 #ifdef IBM 1841 if ((word0(rv) & Exp_mask) < P*Exp_msk1) 1842 #else 1843 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) 1844 #endif 1845 { 1846 if (word0(rv0) == Tiny0 1847 && word1(rv0) == Tiny1) 1848 goto undfl; 1849 word0(rv) = Tiny0; 1850 word1(rv) = Tiny1; 1851 goto cont; 1852 } 1853 else 1854 word0(rv) -= P*Exp_msk1; 1855 } 1856 else { 1857 adj = aadj1 * ulp(value(rv)); 1858 value(rv) += adj; 1859 } 1860 #else 1861 /* Compute adj so that the IEEE rounding rules will 1862 * correctly round rv + adj in some half-way cases. 1863 * If rv * ulp(rv) is denormalized (i.e., 1864 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 1865 * trouble from bits lost to denormalization; 1866 * example: 1.2e-307 . 1867 */ 1868 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) { 1869 aadj1 = (double)(int)(aadj + 0.5); 1870 if (!dsign) 1871 aadj1 = -aadj1; 1872 } 1873 adj = aadj1 * ulp(value(rv)); 1874 value(rv) += adj; 1875 #endif 1876 } 1877 z = word0(rv) & Exp_mask; 1878 if (y == z) { 1879 /* Can we stop now? */ 1880 L = aadj; 1881 aadj -= L; 1882 /* The tolerances below are conservative. */ 1883 if (dsign || word1(rv) || word0(rv) & Bndry_mask) { 1884 if (aadj < .4999999 || aadj > .5000001) 1885 break; 1886 } 1887 else if (aadj < .4999999/FLT_RADIX) 1888 break; 1889 } 1890 cont: 1891 Bfree(bb); 1892 Bfree(bd); 1893 Bfree(bs); 1894 Bfree(delta); 1895 } 1896 retfree: 1897 Bfree(bb); 1898 Bfree(bd); 1899 Bfree(bs); 1900 Bfree(bd0); 1901 Bfree(delta); 1902 ret: 1903 if (se) 1904 /* LINTED interface specification */ 1905 *se = (char *)s; 1906 return sign ? -value(rv) : value(rv); 1907 } 1908 1909 static int 1910 quorem 1911 #ifdef KR_headers 1912 (b, S) Bigint *b, *S; 1913 #else 1914 (Bigint *b, Bigint *S) 1915 #endif 1916 { 1917 int n; 1918 Long borrow, y; 1919 ULong carry, q, ys; 1920 ULong *bx, *bxe, *sx, *sxe; 1921 #ifdef Pack_32 1922 Long z; 1923 ULong si, zs; 1924 #endif 1925 1926 if (b == BIGINT_INVALID || S == BIGINT_INVALID) 1927 return 0; 1928 1929 n = S->wds; 1930 #ifdef DEBUG 1931 /*debug*/ if (b->wds > n) 1932 /*debug*/ Bug("oversize b in quorem"); 1933 #endif 1934 if (b->wds < n) 1935 return 0; 1936 sx = S->x; 1937 sxe = sx + --n; 1938 bx = b->x; 1939 bxe = bx + n; 1940 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ 1941 #ifdef DEBUG 1942 /*debug*/ if (q > 9) 1943 /*debug*/ Bug("oversized quotient in quorem"); 1944 #endif 1945 if (q) { 1946 borrow = 0; 1947 carry = 0; 1948 do { 1949 #ifdef Pack_32 1950 si = *sx++; 1951 ys = (si & 0xffff) * q + carry; 1952 zs = (si >> 16) * q + (ys >> 16); 1953 carry = zs >> 16; 1954 y = (*bx & 0xffff) - (ys & 0xffff) + borrow; 1955 borrow = (ULong)y >> 16; 1956 Sign_Extend(borrow, y); 1957 z = (*bx >> 16) - (zs & 0xffff) + borrow; 1958 borrow = (ULong)z >> 16; 1959 Sign_Extend(borrow, z); 1960 Storeinc(bx, z, y); 1961 #else 1962 ys = *sx++ * q + carry; 1963 carry = ys >> 16; 1964 y = *bx - (ys & 0xffff) + borrow; 1965 borrow = y >> 16; 1966 Sign_Extend(borrow, y); 1967 *bx++ = y & 0xffff; 1968 #endif 1969 } 1970 while(sx <= sxe); 1971 if (!*bxe) { 1972 bx = b->x; 1973 while(--bxe > bx && !*bxe) 1974 --n; 1975 b->wds = n; 1976 } 1977 } 1978 if (cmp(b, S) >= 0) { 1979 q++; 1980 borrow = 0; 1981 carry = 0; 1982 bx = b->x; 1983 sx = S->x; 1984 do { 1985 #ifdef Pack_32 1986 si = *sx++; 1987 ys = (si & 0xffff) + carry; 1988 zs = (si >> 16) + (ys >> 16); 1989 carry = zs >> 16; 1990 y = (*bx & 0xffff) - (ys & 0xffff) + borrow; 1991 borrow = (ULong)y >> 16; 1992 Sign_Extend(borrow, y); 1993 z = (*bx >> 16) - (zs & 0xffff) + borrow; 1994 borrow = (ULong)z >> 16; 1995 Sign_Extend(borrow, z); 1996 Storeinc(bx, z, y); 1997 #else 1998 ys = *sx++ + carry; 1999 carry = ys >> 16; 2000 y = *bx - (ys & 0xffff) + borrow; 2001 borrow = y >> 16; 2002 Sign_Extend(borrow, y); 2003 *bx++ = y & 0xffff; 2004 #endif 2005 } 2006 while(sx <= sxe); 2007 bx = b->x; 2008 bxe = bx + n; 2009 if (!*bxe) { 2010 while(--bxe > bx && !*bxe) 2011 --n; 2012 b->wds = n; 2013 } 2014 } 2015 return q; 2016 } 2017 2018 /* freedtoa(s) must be used to free values s returned by dtoa 2019 * when MULTIPLE_THREADS is #defined. It should be used in all cases, 2020 * but for consistency with earlier versions of dtoa, it is optional 2021 * when MULTIPLE_THREADS is not defined. 2022 */ 2023 2024 void 2025 #ifdef KR_headers 2026 freedtoa(s) char *s; 2027 #else 2028 freedtoa(char *s) 2029 #endif 2030 { 2031 free(s); 2032 } 2033 2034 2035 2036 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. 2037 * 2038 * Inspired by "How to Print Floating-Point Numbers Accurately" by 2039 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101]. 2040 * 2041 * Modifications: 2042 * 1. Rather than iterating, we use a simple numeric overestimate 2043 * to determine k = floor(log10(d)). We scale relevant 2044 * quantities using O(log2(k)) rather than O(k) multiplications. 2045 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't 2046 * try to generate digits strictly left to right. Instead, we 2047 * compute with fewer bits and propagate the carry if necessary 2048 * when rounding the final digit up. This is often faster. 2049 * 3. Under the assumption that input will be rounded nearest, 2050 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. 2051 * That is, we allow equality in stopping tests when the 2052 * round-nearest rule will give the same floating-point value 2053 * as would satisfaction of the stopping test with strict 2054 * inequality. 2055 * 4. We remove common factors of powers of 2 from relevant 2056 * quantities. 2057 * 5. When converting floating-point integers less than 1e16, 2058 * we use floating-point arithmetic rather than resorting 2059 * to multiple-precision integers. 2060 * 6. When asked to produce fewer than 15 digits, we first try 2061 * to get by with floating-point arithmetic; we resort to 2062 * multiple-precision integer arithmetic only if we cannot 2063 * guarantee that the floating-point calculation has given 2064 * the correctly rounded result. For k requested digits and 2065 * "uniformly" distributed input, the probability is 2066 * something like 10^(k-15) that we must resort to the Long 2067 * calculation. 2068 */ 2069 2070 __LIBC_HIDDEN__ char * 2071 __dtoa 2072 #ifdef KR_headers 2073 (_d, mode, ndigits, decpt, sign, rve) 2074 double _d; int mode, ndigits, *decpt, *sign; char **rve; 2075 #else 2076 (double _d, int mode, int ndigits, int *decpt, int *sign, char **rve) 2077 #endif 2078 { 2079 /* Arguments ndigits, decpt, sign are similar to those 2080 of ecvt and fcvt; trailing zeros are suppressed from 2081 the returned string. If not null, *rve is set to point 2082 to the end of the return value. If d is +-Infinity or NaN, 2083 then *decpt is set to 9999. 2084 2085 mode: 2086 0 ==> shortest string that yields d when read in 2087 and rounded to nearest. 2088 1 ==> like 0, but with Steele & White stopping rule; 2089 e.g. with IEEE P754 arithmetic , mode 0 gives 2090 1e23 whereas mode 1 gives 9.999999999999999e22. 2091 2 ==> max(1,ndigits) significant digits. This gives a 2092 return value similar to that of ecvt, except 2093 that trailing zeros are suppressed. 2094 3 ==> through ndigits past the decimal point. This 2095 gives a return value similar to that from fcvt, 2096 except that trailing zeros are suppressed, and 2097 ndigits can be negative. 2098 4-9 should give the same return values as 2-3, i.e., 2099 4 <= mode <= 9 ==> same return as mode 2100 2 + (mode & 1). These modes are mainly for 2101 debugging; often they run slower but sometimes 2102 faster than modes 2-3. 2103 4,5,8,9 ==> left-to-right digit generation. 2104 6-9 ==> don't try fast floating-point estimate 2105 (if applicable). 2106 2107 Values of mode other than 0-9 are treated as mode 0. 2108 2109 Sufficient space is allocated to the return value 2110 to hold the suppressed trailing zeros. 2111 */ 2112 2113 int bbits, b2, b5, be, dig, i, ieps, ilim0, 2114 j, jj1, k, k0, k_check, leftright, m2, m5, s2, s5, 2115 try_quick; 2116 int ilim = 0, ilim1 = 0, spec_case = 0; /* pacify gcc */ 2117 Long L; 2118 #ifndef Sudden_Underflow 2119 int denorm; 2120 ULong x; 2121 #endif 2122 Bigint *b, *b1, *delta, *mhi, *S; 2123 Bigint *mlo = NULL; /* pacify gcc */ 2124 double ds; 2125 char *s, *s0; 2126 Bigint *result = NULL; 2127 int result_k = 0; 2128 _double d, d2, eps; 2129 2130 value(d) = _d; 2131 2132 if (word0(d) & Sign_bit) { 2133 /* set sign for everything, including 0's and NaNs */ 2134 *sign = 1; 2135 word0(d) &= ~Sign_bit; /* clear sign bit */ 2136 } 2137 else 2138 *sign = 0; 2139 2140 #if defined(IEEE_Arith) + defined(VAX) 2141 #ifdef IEEE_Arith 2142 if ((word0(d) & Exp_mask) == Exp_mask) 2143 #else 2144 if (word0(d) == 0x8000) 2145 #endif 2146 { 2147 /* Infinity or NaN */ 2148 *decpt = 9999; 2149 s = 2150 #ifdef IEEE_Arith 2151 !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" : 2152 #endif 2153 "NaN"; 2154 result = Balloc(strlen(s)+1); 2155 if (result == BIGINT_INVALID) 2156 return NULL; 2157 s0 = (char *)(void *)result; 2158 strcpy(s0, s); 2159 if (rve) 2160 *rve = 2161 #ifdef IEEE_Arith 2162 s0[3] ? s0 + 8 : 2163 #endif 2164 s0 + 3; 2165 return s0; 2166 } 2167 #endif 2168 #ifdef IBM 2169 value(d) += 0; /* normalize */ 2170 #endif 2171 if (!value(d)) { 2172 *decpt = 1; 2173 result = Balloc(2); 2174 if (result == BIGINT_INVALID) 2175 return NULL; 2176 s0 = (char *)(void *)result; 2177 strcpy(s0, "0"); 2178 if (rve) 2179 *rve = s0 + 1; 2180 return s0; 2181 } 2182 2183 b = d2b(value(d), &be, &bbits); 2184 #ifdef Sudden_Underflow 2185 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)); 2186 #else 2187 if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) != 0) { 2188 #endif 2189 value(d2) = value(d); 2190 word0(d2) &= Frac_mask1; 2191 word0(d2) |= Exp_11; 2192 #ifdef IBM 2193 if (j = 11 - hi0bits(word0(d2) & Frac_mask)) 2194 value(d2) /= 1 << j; 2195 #endif 2196 2197 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 2198 * log10(x) = log(x) / log(10) 2199 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) 2200 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) 2201 * 2202 * This suggests computing an approximation k to log10(d) by 2203 * 2204 * k = (i - Bias)*0.301029995663981 2205 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); 2206 * 2207 * We want k to be too large rather than too small. 2208 * The error in the first-order Taylor series approximation 2209 * is in our favor, so we just round up the constant enough 2210 * to compensate for any error in the multiplication of 2211 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, 2212 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, 2213 * adding 1e-13 to the constant term more than suffices. 2214 * Hence we adjust the constant term to 0.1760912590558. 2215 * (We could get a more accurate k by invoking log10, 2216 * but this is probably not worthwhile.) 2217 */ 2218 2219 i -= Bias; 2220 #ifdef IBM 2221 i <<= 2; 2222 i += j; 2223 #endif 2224 #ifndef Sudden_Underflow 2225 denorm = 0; 2226 } 2227 else { 2228 /* d is denormalized */ 2229 2230 i = bbits + be + (Bias + (P-1) - 1); 2231 x = i > 32 ? word0(d) << (64 - i) | word1(d) >> (i - 32) 2232 : word1(d) << (32 - i); 2233 value(d2) = x; 2234 word0(d2) -= 31*Exp_msk1; /* adjust exponent */ 2235 i -= (Bias + (P-1) - 1) + 1; 2236 denorm = 1; 2237 } 2238 #endif 2239 ds = (value(d2)-1.5)*0.289529654602168 + 0.1760912590558 + 2240 i*0.301029995663981; 2241 k = (int)ds; 2242 if (ds < 0. && ds != k) 2243 k--; /* want k = floor(ds) */ 2244 k_check = 1; 2245 if (k >= 0 && k <= Ten_pmax) { 2246 if (value(d) < tens[k]) 2247 k--; 2248 k_check = 0; 2249 } 2250 j = bbits - i - 1; 2251 if (j >= 0) { 2252 b2 = 0; 2253 s2 = j; 2254 } 2255 else { 2256 b2 = -j; 2257 s2 = 0; 2258 } 2259 if (k >= 0) { 2260 b5 = 0; 2261 s5 = k; 2262 s2 += k; 2263 } 2264 else { 2265 b2 -= k; 2266 b5 = -k; 2267 s5 = 0; 2268 } 2269 if (mode < 0 || mode > 9) 2270 mode = 0; 2271 try_quick = 1; 2272 if (mode > 5) { 2273 mode -= 4; 2274 try_quick = 0; 2275 } 2276 leftright = 1; 2277 switch(mode) { 2278 case 0: 2279 case 1: 2280 ilim = ilim1 = -1; 2281 i = 18; 2282 ndigits = 0; 2283 break; 2284 case 2: 2285 leftright = 0; 2286 /* FALLTHROUGH */ 2287 case 4: 2288 if (ndigits <= 0) 2289 ndigits = 1; 2290 ilim = ilim1 = i = ndigits; 2291 break; 2292 case 3: 2293 leftright = 0; 2294 /* FALLTHROUGH */ 2295 case 5: 2296 i = ndigits + k + 1; 2297 ilim = i; 2298 ilim1 = i - 1; 2299 if (i <= 0) 2300 i = 1; 2301 } 2302 j = sizeof(ULong); 2303 for(result_k = 0; (int)(sizeof(Bigint) - sizeof(ULong)) + j <= i; 2304 j <<= 1) result_k++; 2305 // this is really a ugly hack, the code uses Balloc 2306 // instead of malloc, but casts the result into a char* 2307 // it seems the only reason to do that is due to the 2308 // complicated way the block size need to be computed 2309 // buuurk.... 2310 result = Balloc(result_k); 2311 if (result == BIGINT_INVALID) { 2312 Bfree(b); 2313 return NULL; 2314 } 2315 s = s0 = (char *)(void *)result; 2316 2317 if (ilim >= 0 && ilim <= Quick_max && try_quick) { 2318 2319 /* Try to get by with floating-point arithmetic. */ 2320 2321 i = 0; 2322 value(d2) = value(d); 2323 k0 = k; 2324 ilim0 = ilim; 2325 ieps = 2; /* conservative */ 2326 if (k > 0) { 2327 ds = tens[k&0xf]; 2328 j = (unsigned int)k >> 4; 2329 if (j & Bletch) { 2330 /* prevent overflows */ 2331 j &= Bletch - 1; 2332 value(d) /= bigtens[n_bigtens-1]; 2333 ieps++; 2334 } 2335 for(; j; j = (unsigned int)j >> 1, i++) 2336 if (j & 1) { 2337 ieps++; 2338 ds *= bigtens[i]; 2339 } 2340 value(d) /= ds; 2341 } 2342 else if ((jj1 = -k) != 0) { 2343 value(d) *= tens[jj1 & 0xf]; 2344 for(j = (unsigned int)jj1 >> 4; j; 2345 j = (unsigned int)j >> 1, i++) 2346 if (j & 1) { 2347 ieps++; 2348 value(d) *= bigtens[i]; 2349 } 2350 } 2351 if (k_check && value(d) < 1. && ilim > 0) { 2352 if (ilim1 <= 0) 2353 goto fast_failed; 2354 ilim = ilim1; 2355 k--; 2356 value(d) *= 10.; 2357 ieps++; 2358 } 2359 value(eps) = ieps*value(d) + 7.; 2360 word0(eps) -= (P-1)*Exp_msk1; 2361 if (ilim == 0) { 2362 S = mhi = 0; 2363 value(d) -= 5.; 2364 if (value(d) > value(eps)) 2365 goto one_digit; 2366 if (value(d) < -value(eps)) 2367 goto no_digits; 2368 goto fast_failed; 2369 } 2370 #ifndef No_leftright 2371 if (leftright) { 2372 /* Use Steele & White method of only 2373 * generating digits needed. 2374 */ 2375 value(eps) = 0.5/tens[ilim-1] - value(eps); 2376 for(i = 0;;) { 2377 L = value(d); 2378 value(d) -= L; 2379 *s++ = '0' + (int)L; 2380 if (value(d) < value(eps)) 2381 goto ret1; 2382 if (1. - value(d) < value(eps)) 2383 goto bump_up; 2384 if (++i >= ilim) 2385 break; 2386 value(eps) *= 10.; 2387 value(d) *= 10.; 2388 } 2389 } 2390 else { 2391 #endif 2392 /* Generate ilim digits, then fix them up. */ 2393 value(eps) *= tens[ilim-1]; 2394 for(i = 1;; i++, value(d) *= 10.) { 2395 L = value(d); 2396 value(d) -= L; 2397 *s++ = '0' + (int)L; 2398 if (i == ilim) { 2399 if (value(d) > 0.5 + value(eps)) 2400 goto bump_up; 2401 else if (value(d) < 0.5 - value(eps)) { 2402 while(*--s == '0'); 2403 s++; 2404 goto ret1; 2405 } 2406 break; 2407 } 2408 } 2409 #ifndef No_leftright 2410 } 2411 #endif 2412 fast_failed: 2413 s = s0; 2414 value(d) = value(d2); 2415 k = k0; 2416 ilim = ilim0; 2417 } 2418 2419 /* Do we have a "small" integer? */ 2420 2421 if (be >= 0 && k <= Int_max) { 2422 /* Yes. */ 2423 ds = tens[k]; 2424 if (ndigits < 0 && ilim <= 0) { 2425 S = mhi = 0; 2426 if (ilim < 0 || value(d) <= 5*ds) 2427 goto no_digits; 2428 goto one_digit; 2429 } 2430 for(i = 1;; i++) { 2431 L = value(d) / ds; 2432 value(d) -= L*ds; 2433 #ifdef Check_FLT_ROUNDS 2434 /* If FLT_ROUNDS == 2, L will usually be high by 1 */ 2435 if (value(d) < 0) { 2436 L--; 2437 value(d) += ds; 2438 } 2439 #endif 2440 *s++ = '0' + (int)L; 2441 if (i == ilim) { 2442 value(d) += value(d); 2443 if (value(d) > ds || (value(d) == ds && L & 1)) { 2444 bump_up: 2445 while(*--s == '9') 2446 if (s == s0) { 2447 k++; 2448 *s = '0'; 2449 break; 2450 } 2451 ++*s++; 2452 } 2453 break; 2454 } 2455 if (!(value(d) *= 10.)) 2456 break; 2457 } 2458 goto ret1; 2459 } 2460 2461 m2 = b2; 2462 m5 = b5; 2463 mhi = mlo = 0; 2464 if (leftright) { 2465 if (mode < 2) { 2466 i = 2467 #ifndef Sudden_Underflow 2468 denorm ? be + (Bias + (P-1) - 1 + 1) : 2469 #endif 2470 #ifdef IBM 2471 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3); 2472 #else 2473 1 + P - bbits; 2474 #endif 2475 } 2476 else { 2477 j = ilim - 1; 2478 if (m5 >= j) 2479 m5 -= j; 2480 else { 2481 s5 += j -= m5; 2482 b5 += j; 2483 m5 = 0; 2484 } 2485 if ((i = ilim) < 0) { 2486 m2 -= i; 2487 i = 0; 2488 } 2489 } 2490 b2 += i; 2491 s2 += i; 2492 mhi = i2b(1); 2493 } 2494 if (m2 > 0 && s2 > 0) { 2495 i = m2 < s2 ? m2 : s2; 2496 b2 -= i; 2497 m2 -= i; 2498 s2 -= i; 2499 } 2500 if (b5 > 0) { 2501 if (leftright) { 2502 if (m5 > 0) { 2503 mhi = pow5mult(mhi, m5); 2504 b1 = mult(mhi, b); 2505 Bfree(b); 2506 b = b1; 2507 } 2508 if ((j = b5 - m5) != 0) 2509 b = pow5mult(b, j); 2510 } 2511 else 2512 b = pow5mult(b, b5); 2513 } 2514 S = i2b(1); 2515 if (s5 > 0) 2516 S = pow5mult(S, s5); 2517 2518 /* Check for special case that d is a normalized power of 2. */ 2519 2520 if (mode < 2) { 2521 if (!word1(d) && !(word0(d) & Bndry_mask) 2522 #ifndef Sudden_Underflow 2523 && word0(d) & Exp_mask 2524 #endif 2525 ) { 2526 /* The special case */ 2527 b2 += Log2P; 2528 s2 += Log2P; 2529 spec_case = 1; 2530 } 2531 else 2532 spec_case = 0; 2533 } 2534 2535 /* Arrange for convenient computation of quotients: 2536 * shift left if necessary so divisor has 4 leading 0 bits. 2537 * 2538 * Perhaps we should just compute leading 28 bits of S once 2539 * and for all and pass them and a shift to quorem, so it 2540 * can do shifts and ors to compute the numerator for q. 2541 */ 2542 if (S == BIGINT_INVALID) { 2543 i = 0; 2544 } else { 2545 #ifdef Pack_32 2546 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) != 0) 2547 i = 32 - i; 2548 #else 2549 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) 2550 i = 16 - i; 2551 #endif 2552 } 2553 2554 if (i > 4) { 2555 i -= 4; 2556 b2 += i; 2557 m2 += i; 2558 s2 += i; 2559 } 2560 else if (i < 4) { 2561 i += 28; 2562 b2 += i; 2563 m2 += i; 2564 s2 += i; 2565 } 2566 if (b2 > 0) 2567 b = lshift(b, b2); 2568 if (s2 > 0) 2569 S = lshift(S, s2); 2570 if (k_check) { 2571 if (cmp(b,S) < 0) { 2572 k--; 2573 b = multadd(b, 10, 0); /* we botched the k estimate */ 2574 if (leftright) 2575 mhi = multadd(mhi, 10, 0); 2576 ilim = ilim1; 2577 } 2578 } 2579 if (ilim <= 0 && mode > 2) { 2580 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) { 2581 /* no digits, fcvt style */ 2582 no_digits: 2583 k = -1 - ndigits; 2584 goto ret; 2585 } 2586 one_digit: 2587 *s++ = '1'; 2588 k++; 2589 goto ret; 2590 } 2591 if (leftright) { 2592 if (m2 > 0) 2593 mhi = lshift(mhi, m2); 2594 2595 /* Compute mlo -- check for special case 2596 * that d is a normalized power of 2. 2597 */ 2598 2599 mlo = mhi; 2600 if (spec_case) { 2601 mhi = Balloc(mhi->k); 2602 Bcopy(mhi, mlo); 2603 mhi = lshift(mhi, Log2P); 2604 } 2605 2606 for(i = 1;;i++) { 2607 dig = quorem(b,S) + '0'; 2608 /* Do we yet have the shortest decimal string 2609 * that will round to d? 2610 */ 2611 j = cmp(b, mlo); 2612 delta = diff(S, mhi); 2613 jj1 = delta->sign ? 1 : cmp(b, delta); 2614 Bfree(delta); 2615 #ifndef ROUND_BIASED 2616 if (jj1 == 0 && !mode && !(word1(d) & 1)) { 2617 if (dig == '9') 2618 goto round_9_up; 2619 if (j > 0) 2620 dig++; 2621 *s++ = dig; 2622 goto ret; 2623 } 2624 #endif 2625 if (j < 0 || (j == 0 && !mode 2626 #ifndef ROUND_BIASED 2627 && !(word1(d) & 1) 2628 #endif 2629 )) { 2630 if (jj1 > 0) { 2631 b = lshift(b, 1); 2632 jj1 = cmp(b, S); 2633 if ((jj1 > 0 || (jj1 == 0 && dig & 1)) 2634 && dig++ == '9') 2635 goto round_9_up; 2636 } 2637 *s++ = dig; 2638 goto ret; 2639 } 2640 if (jj1 > 0) { 2641 if (dig == '9') { /* possible if i == 1 */ 2642 round_9_up: 2643 *s++ = '9'; 2644 goto roundoff; 2645 } 2646 *s++ = dig + 1; 2647 goto ret; 2648 } 2649 *s++ = dig; 2650 if (i == ilim) 2651 break; 2652 b = multadd(b, 10, 0); 2653 if (mlo == mhi) 2654 mlo = mhi = multadd(mhi, 10, 0); 2655 else { 2656 mlo = multadd(mlo, 10, 0); 2657 mhi = multadd(mhi, 10, 0); 2658 } 2659 } 2660 } 2661 else 2662 for(i = 1;; i++) { 2663 *s++ = dig = quorem(b,S) + '0'; 2664 if (i >= ilim) 2665 break; 2666 b = multadd(b, 10, 0); 2667 } 2668 2669 /* Round off last digit */ 2670 2671 b = lshift(b, 1); 2672 j = cmp(b, S); 2673 if (j > 0 || (j == 0 && dig & 1)) { 2674 roundoff: 2675 while(*--s == '9') 2676 if (s == s0) { 2677 k++; 2678 *s++ = '1'; 2679 goto ret; 2680 } 2681 ++*s++; 2682 } 2683 else { 2684 while(*--s == '0'); 2685 s++; 2686 } 2687 ret: 2688 Bfree(S); 2689 if (mhi) { 2690 if (mlo && mlo != mhi) 2691 Bfree(mlo); 2692 Bfree(mhi); 2693 } 2694 ret1: 2695 Bfree(b); 2696 if (s == s0) { /* don't return empty string */ 2697 *s++ = '0'; 2698 k = 0; 2699 } 2700 *s = 0; 2701 *decpt = k + 1; 2702 if (rve) 2703 *rve = s; 2704 return s0; 2705 } 2706 #ifdef __cplusplus 2707 } 2708 #endif 2709
