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 758 static Bigint * 759 pow5mult 760 #ifdef KR_headers 761 (b, k) Bigint *b; int k; 762 #else 763 (Bigint *b, int k) 764 #endif 765 { 766 Bigint *b1, *p5, *p51; 767 int i; 768 static const int p05[3] = { 5, 25, 125 }; 769 770 if (b == BIGINT_INVALID) 771 return b; 772 773 if ((i = k & 3) != 0) 774 b = multadd(b, p05[i-1], 0); 775 776 if (!(k = (unsigned int) k >> 2)) 777 return b; 778 if (!(p5 = p5s)) { 779 /* first time */ 780 p5 = i2b(625); 781 if (p5 == BIGINT_INVALID) { 782 Bfree(b); 783 return p5; 784 } 785 p5s = p5; 786 p5->next = 0; 787 } 788 for(;;) { 789 if (k & 1) { 790 b1 = mult(b, p5); 791 Bfree(b); 792 b = b1; 793 } 794 if (!(k = (unsigned int) k >> 1)) 795 break; 796 if (!(p51 = p5->next)) { 797 p51 = mult(p5,p5); 798 if (p51 == BIGINT_INVALID) { 799 Bfree(b); 800 return p51; 801 } 802 p5->next = p51; 803 p51->next = 0; 804 } 805 p5 = p51; 806 } 807 return b; 808 } 809 810 static Bigint * 811 lshift 812 #ifdef KR_headers 813 (b, k) Bigint *b; int k; 814 #else 815 (Bigint *b, int k) 816 #endif 817 { 818 int i, k1, n, n1; 819 Bigint *b1; 820 ULong *x, *x1, *xe, z; 821 822 if (b == BIGINT_INVALID) 823 return b; 824 825 #ifdef Pack_32 826 n = (unsigned int)k >> 5; 827 #else 828 n = (unsigned int)k >> 4; 829 #endif 830 k1 = b->k; 831 n1 = n + b->wds + 1; 832 for(i = b->maxwds; n1 > i; i <<= 1) 833 k1++; 834 b1 = Balloc(k1); 835 if (b1 == BIGINT_INVALID) { 836 Bfree(b); 837 return b1; 838 } 839 x1 = b1->x; 840 for(i = 0; i < n; i++) 841 *x1++ = 0; 842 x = b->x; 843 xe = x + b->wds; 844 #ifdef Pack_32 845 if (k &= 0x1f) { 846 k1 = 32 - k; 847 z = 0; 848 do { 849 *x1++ = *x << k | z; 850 z = *x++ >> k1; 851 } 852 while(x < xe); 853 if ((*x1 = z) != 0) 854 ++n1; 855 } 856 #else 857 if (k &= 0xf) { 858 k1 = 16 - k; 859 z = 0; 860 do { 861 *x1++ = *x << k & 0xffff | z; 862 z = *x++ >> k1; 863 } 864 while(x < xe); 865 if (*x1 = z) 866 ++n1; 867 } 868 #endif 869 else do 870 *x1++ = *x++; 871 while(x < xe); 872 b1->wds = n1 - 1; 873 Bfree(b); 874 return b1; 875 } 876 877 static int 878 cmp 879 #ifdef KR_headers 880 (a, b) Bigint *a, *b; 881 #else 882 (Bigint *a, Bigint *b) 883 #endif 884 { 885 ULong *xa, *xa0, *xb, *xb0; 886 int i, j; 887 888 if (a == BIGINT_INVALID || b == BIGINT_INVALID) 889 #ifdef DEBUG 890 Bug("cmp called with a or b invalid"); 891 #else 892 return 0; /* equal - the best we can do right now */ 893 #endif 894 895 i = a->wds; 896 j = b->wds; 897 #ifdef DEBUG 898 if (i > 1 && !a->x[i-1]) 899 Bug("cmp called with a->x[a->wds-1] == 0"); 900 if (j > 1 && !b->x[j-1]) 901 Bug("cmp called with b->x[b->wds-1] == 0"); 902 #endif 903 if (i -= j) 904 return i; 905 xa0 = a->x; 906 xa = xa0 + j; 907 xb0 = b->x; 908 xb = xb0 + j; 909 for(;;) { 910 if (*--xa != *--xb) 911 return *xa < *xb ? -1 : 1; 912 if (xa <= xa0) 913 break; 914 } 915 return 0; 916 } 917 918 static Bigint * 919 diff 920 #ifdef KR_headers 921 (a, b) Bigint *a, *b; 922 #else 923 (Bigint *a, Bigint *b) 924 #endif 925 { 926 Bigint *c; 927 int i, wa, wb; 928 Long borrow, y; /* We need signed shifts here. */ 929 ULong *xa, *xae, *xb, *xbe, *xc; 930 #ifdef Pack_32 931 Long z; 932 #endif 933 934 if (a == BIGINT_INVALID || b == BIGINT_INVALID) 935 return BIGINT_INVALID; 936 937 i = cmp(a,b); 938 if (!i) { 939 c = Balloc(0); 940 if (c != BIGINT_INVALID) { 941 c->wds = 1; 942 c->x[0] = 0; 943 } 944 return c; 945 } 946 if (i < 0) { 947 c = a; 948 a = b; 949 b = c; 950 i = 1; 951 } 952 else 953 i = 0; 954 c = Balloc(a->k); 955 if (c == BIGINT_INVALID) 956 return c; 957 c->sign = i; 958 wa = a->wds; 959 xa = a->x; 960 xae = xa + wa; 961 wb = b->wds; 962 xb = b->x; 963 xbe = xb + wb; 964 xc = c->x; 965 borrow = 0; 966 #ifdef Pack_32 967 do { 968 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow; 969 borrow = (ULong)y >> 16; 970 Sign_Extend(borrow, y); 971 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow; 972 borrow = (ULong)z >> 16; 973 Sign_Extend(borrow, z); 974 Storeinc(xc, z, y); 975 } 976 while(xb < xbe); 977 while(xa < xae) { 978 y = (*xa & 0xffff) + borrow; 979 borrow = (ULong)y >> 16; 980 Sign_Extend(borrow, y); 981 z = (*xa++ >> 16) + borrow; 982 borrow = (ULong)z >> 16; 983 Sign_Extend(borrow, z); 984 Storeinc(xc, z, y); 985 } 986 #else 987 do { 988 y = *xa++ - *xb++ + borrow; 989 borrow = y >> 16; 990 Sign_Extend(borrow, y); 991 *xc++ = y & 0xffff; 992 } 993 while(xb < xbe); 994 while(xa < xae) { 995 y = *xa++ + borrow; 996 borrow = y >> 16; 997 Sign_Extend(borrow, y); 998 *xc++ = y & 0xffff; 999 } 1000 #endif 1001 while(!*--xc) 1002 wa--; 1003 c->wds = wa; 1004 return c; 1005 } 1006 1007 static double 1008 ulp 1009 #ifdef KR_headers 1010 (_x) double _x; 1011 #else 1012 (double _x) 1013 #endif 1014 { 1015 _double x; 1016 Long L; 1017 _double a; 1018 1019 value(x) = _x; 1020 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; 1021 #ifndef Sudden_Underflow 1022 if (L > 0) { 1023 #endif 1024 #ifdef IBM 1025 L |= Exp_msk1 >> 4; 1026 #endif 1027 word0(a) = L; 1028 word1(a) = 0; 1029 #ifndef Sudden_Underflow 1030 } 1031 else { 1032 L = (ULong)-L >> Exp_shift; 1033 if (L < Exp_shift) { 1034 word0(a) = 0x80000 >> L; 1035 word1(a) = 0; 1036 } 1037 else { 1038 word0(a) = 0; 1039 L -= Exp_shift; 1040 word1(a) = L >= 31 ? 1 : 1 << (31 - L); 1041 } 1042 } 1043 #endif 1044 return value(a); 1045 } 1046 1047 static double 1048 b2d 1049 #ifdef KR_headers 1050 (a, e) Bigint *a; int *e; 1051 #else 1052 (Bigint *a, int *e) 1053 #endif 1054 { 1055 ULong *xa, *xa0, w, y, z; 1056 int k; 1057 _double d; 1058 #ifdef VAX 1059 ULong d0, d1; 1060 #else 1061 #define d0 word0(d) 1062 #define d1 word1(d) 1063 #endif 1064 1065 if (a == BIGINT_INVALID) 1066 return NAN; 1067 1068 xa0 = a->x; 1069 xa = xa0 + a->wds; 1070 y = *--xa; 1071 #ifdef DEBUG 1072 if (!y) Bug("zero y in b2d"); 1073 #endif 1074 k = hi0bits(y); 1075 *e = 32 - k; 1076 #ifdef Pack_32 1077 if (k < Ebits) { 1078 d0 = Exp_1 | y >> (Ebits - k); 1079 w = xa > xa0 ? *--xa : 0; 1080 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k); 1081 goto ret_d; 1082 } 1083 z = xa > xa0 ? *--xa : 0; 1084 if (k -= Ebits) { 1085 d0 = Exp_1 | y << k | z >> (32 - k); 1086 y = xa > xa0 ? *--xa : 0; 1087 d1 = z << k | y >> (32 - k); 1088 } 1089 else { 1090 d0 = Exp_1 | y; 1091 d1 = z; 1092 } 1093 #else 1094 if (k < Ebits + 16) { 1095 z = xa > xa0 ? *--xa : 0; 1096 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; 1097 w = xa > xa0 ? *--xa : 0; 1098 y = xa > xa0 ? *--xa : 0; 1099 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; 1100 goto ret_d; 1101 } 1102 z = xa > xa0 ? *--xa : 0; 1103 w = xa > xa0 ? *--xa : 0; 1104 k -= Ebits + 16; 1105 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; 1106 y = xa > xa0 ? *--xa : 0; 1107 d1 = w << k + 16 | y << k; 1108 #endif 1109 ret_d: 1110 #ifdef VAX 1111 word0(d) = d0 >> 16 | d0 << 16; 1112 word1(d) = d1 >> 16 | d1 << 16; 1113 #else 1114 #undef d0 1115 #undef d1 1116 #endif 1117 return value(d); 1118 } 1119 1120 static Bigint * 1121 d2b 1122 #ifdef KR_headers 1123 (_d, e, bits) double d; int *e, *bits; 1124 #else 1125 (double _d, int *e, int *bits) 1126 #endif 1127 { 1128 Bigint *b; 1129 int de, i, k; 1130 ULong *x, y, z; 1131 _double d; 1132 #ifdef VAX 1133 ULong d0, d1; 1134 #endif 1135 1136 value(d) = _d; 1137 #ifdef VAX 1138 d0 = word0(d) >> 16 | word0(d) << 16; 1139 d1 = word1(d) >> 16 | word1(d) << 16; 1140 #else 1141 #define d0 word0(d) 1142 #define d1 word1(d) 1143 #endif 1144 1145 #ifdef Pack_32 1146 b = Balloc(1); 1147 #else 1148 b = Balloc(2); 1149 #endif 1150 if (b == BIGINT_INVALID) 1151 return b; 1152 x = b->x; 1153 1154 z = d0 & Frac_mask; 1155 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 1156 #ifdef Sudden_Underflow 1157 de = (int)(d0 >> Exp_shift); 1158 #ifndef IBM 1159 z |= Exp_msk11; 1160 #endif 1161 #else 1162 if ((de = (int)(d0 >> Exp_shift)) != 0) 1163 z |= Exp_msk1; 1164 #endif 1165 #ifdef Pack_32 1166 if ((y = d1) != 0) { 1167 if ((k = lo0bits(&y)) != 0) { 1168 x[0] = y | z << (32 - k); 1169 z >>= k; 1170 } 1171 else 1172 x[0] = y; 1173 i = b->wds = (x[1] = z) ? 2 : 1; 1174 } 1175 else { 1176 #ifdef DEBUG 1177 if (!z) 1178 Bug("Zero passed to d2b"); 1179 #endif 1180 k = lo0bits(&z); 1181 x[0] = z; 1182 i = b->wds = 1; 1183 k += 32; 1184 } 1185 #else 1186 if (y = d1) { 1187 if (k = lo0bits(&y)) 1188 if (k >= 16) { 1189 x[0] = y | z << 32 - k & 0xffff; 1190 x[1] = z >> k - 16 & 0xffff; 1191 x[2] = z >> k; 1192 i = 2; 1193 } 1194 else { 1195 x[0] = y & 0xffff; 1196 x[1] = y >> 16 | z << 16 - k & 0xffff; 1197 x[2] = z >> k & 0xffff; 1198 x[3] = z >> k+16; 1199 i = 3; 1200 } 1201 else { 1202 x[0] = y & 0xffff; 1203 x[1] = y >> 16; 1204 x[2] = z & 0xffff; 1205 x[3] = z >> 16; 1206 i = 3; 1207 } 1208 } 1209 else { 1210 #ifdef DEBUG 1211 if (!z) 1212 Bug("Zero passed to d2b"); 1213 #endif 1214 k = lo0bits(&z); 1215 if (k >= 16) { 1216 x[0] = z; 1217 i = 0; 1218 } 1219 else { 1220 x[0] = z & 0xffff; 1221 x[1] = z >> 16; 1222 i = 1; 1223 } 1224 k += 32; 1225 } 1226 while(!x[i]) 1227 --i; 1228 b->wds = i + 1; 1229 #endif 1230 #ifndef Sudden_Underflow 1231 if (de) { 1232 #endif 1233 #ifdef IBM 1234 *e = (de - Bias - (P-1) << 2) + k; 1235 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask); 1236 #else 1237 *e = de - Bias - (P-1) + k; 1238 *bits = P - k; 1239 #endif 1240 #ifndef Sudden_Underflow 1241 } 1242 else { 1243 *e = de - Bias - (P-1) + 1 + k; 1244 #ifdef Pack_32 1245 *bits = 32*i - hi0bits(x[i-1]); 1246 #else 1247 *bits = (i+2)*16 - hi0bits(x[i]); 1248 #endif 1249 } 1250 #endif 1251 return b; 1252 } 1253 #undef d0 1254 #undef d1 1255 1256 static double 1257 ratio 1258 #ifdef KR_headers 1259 (a, b) Bigint *a, *b; 1260 #else 1261 (Bigint *a, Bigint *b) 1262 #endif 1263 { 1264 _double da, db; 1265 int k, ka, kb; 1266 1267 if (a == BIGINT_INVALID || b == BIGINT_INVALID) 1268 return NAN; /* for lack of better value ? */ 1269 1270 value(da) = b2d(a, &ka); 1271 value(db) = b2d(b, &kb); 1272 #ifdef Pack_32 1273 k = ka - kb + 32*(a->wds - b->wds); 1274 #else 1275 k = ka - kb + 16*(a->wds - b->wds); 1276 #endif 1277 #ifdef IBM 1278 if (k > 0) { 1279 word0(da) += (k >> 2)*Exp_msk1; 1280 if (k &= 3) 1281 da *= 1 << k; 1282 } 1283 else { 1284 k = -k; 1285 word0(db) += (k >> 2)*Exp_msk1; 1286 if (k &= 3) 1287 db *= 1 << k; 1288 } 1289 #else 1290 if (k > 0) 1291 word0(da) += k*Exp_msk1; 1292 else { 1293 k = -k; 1294 word0(db) += k*Exp_msk1; 1295 } 1296 #endif 1297 return value(da) / value(db); 1298 } 1299 1300 static CONST double 1301 tens[] = { 1302 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1303 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1304 1e20, 1e21, 1e22 1305 #ifdef VAX 1306 , 1e23, 1e24 1307 #endif 1308 }; 1309 1310 #ifdef IEEE_Arith 1311 static CONST double bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; 1312 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 }; 1313 #define n_bigtens 5 1314 #else 1315 #ifdef IBM 1316 static CONST double bigtens[] = { 1e16, 1e32, 1e64 }; 1317 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 }; 1318 #define n_bigtens 3 1319 #else 1320 static CONST double bigtens[] = { 1e16, 1e32 }; 1321 static CONST double tinytens[] = { 1e-16, 1e-32 }; 1322 #define n_bigtens 2 1323 #endif 1324 #endif 1325 1326 double 1327 strtod 1328 #ifdef KR_headers 1329 (s00, se) CONST char *s00; char **se; 1330 #else 1331 (CONST char *s00, char **se) 1332 #endif 1333 { 1334 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, 1335 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; 1336 CONST char *s, *s0, *s1; 1337 double aadj, aadj1, adj; 1338 _double rv, rv0; 1339 Long L; 1340 ULong y, z; 1341 Bigint *bb1, *bd0; 1342 Bigint *bb = NULL, *bd = NULL, *bs = NULL, *delta = NULL;/* pacify gcc */ 1343 1344 #ifdef ANDROID_CHANGES 1345 CONST char decimal_point = '.'; 1346 #else /* ANDROID_CHANGES */ 1347 #ifndef KR_headers 1348 CONST char decimal_point = localeconv()->decimal_point[0]; 1349 #else 1350 CONST char decimal_point = '.'; 1351 #endif 1352 1353 #endif /* ANDROID_CHANGES */ 1354 1355 sign = nz0 = nz = 0; 1356 value(rv) = 0.; 1357 1358 1359 for(s = s00; isspace((unsigned char) *s); s++) 1360 ; 1361 1362 if (*s == '-') { 1363 sign = 1; 1364 s++; 1365 } else if (*s == '+') { 1366 s++; 1367 } 1368 1369 if (*s == '\0') { 1370 s = s00; 1371 goto ret; 1372 } 1373 1374 /* "INF" or "INFINITY" */ 1375 if (tolower((unsigned char)*s) == 'i' && strncasecmp(s, "inf", 3) == 0) { 1376 if (strncasecmp(s + 3, "inity", 5) == 0) 1377 s += 8; 1378 else 1379 s += 3; 1380 1381 value(rv) = HUGE_VAL; 1382 goto ret; 1383 } 1384 1385 #ifdef IEEE_Arith 1386 /* "NAN" or "NAN(n-char-sequence-opt)" */ 1387 if (tolower((unsigned char)*s) == 'n' && strncasecmp(s, "nan", 3) == 0) { 1388 /* Build a quiet NaN. */ 1389 word0(rv) = NAN_WORD0; 1390 word1(rv) = NAN_WORD1; 1391 s+= 3; 1392 1393 /* Don't interpret (n-char-sequence-opt), for now. */ 1394 if (*s == '(') { 1395 s0 = s; 1396 for (s++; *s != ')' && *s != '\0'; s++) 1397 ; 1398 if (*s == ')') 1399 s++; /* Skip over closing paren ... */ 1400 else 1401 s = s0; /* ... otherwise go back. */ 1402 } 1403 1404 goto ret; 1405 } 1406 #endif 1407 1408 if (*s == '0') { 1409 nz0 = 1; 1410 while(*++s == '0') ; 1411 if (!*s) 1412 goto ret; 1413 } 1414 s0 = s; 1415 y = z = 0; 1416 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 1417 if (nd < 9) 1418 y = 10*y + c - '0'; 1419 else if (nd < 16) 1420 z = 10*z + c - '0'; 1421 nd0 = nd; 1422 if (c == decimal_point) { 1423 c = *++s; 1424 if (!nd) { 1425 for(; c == '0'; c = *++s) 1426 nz++; 1427 if (c > '0' && c <= '9') { 1428 s0 = s; 1429 nf += nz; 1430 nz = 0; 1431 goto have_dig; 1432 } 1433 goto dig_done; 1434 } 1435 for(; c >= '0' && c <= '9'; c = *++s) { 1436 have_dig: 1437 nz++; 1438 if (c -= '0') { 1439 nf += nz; 1440 for(i = 1; i < nz; i++) 1441 if (nd++ < 9) 1442 y *= 10; 1443 else if (nd <= DBL_DIG + 1) 1444 z *= 10; 1445 if (nd++ < 9) 1446 y = 10*y + c; 1447 else if (nd <= DBL_DIG + 1) 1448 z = 10*z + c; 1449 nz = 0; 1450 } 1451 } 1452 } 1453 dig_done: 1454 e = 0; 1455 if (c == 'e' || c == 'E') { 1456 if (!nd && !nz && !nz0) { 1457 s = s00; 1458 goto ret; 1459 } 1460 s00 = s; 1461 esign = 0; 1462 switch(c = *++s) { 1463 case '-': 1464 esign = 1; 1465 /* FALLTHROUGH */ 1466 case '+': 1467 c = *++s; 1468 } 1469 if (c >= '0' && c <= '9') { 1470 while(c == '0') 1471 c = *++s; 1472 if (c > '0' && c <= '9') { 1473 L = c - '0'; 1474 s1 = s; 1475 while((c = *++s) >= '0' && c <= '9') 1476 L = 10*L + c - '0'; 1477 if (s - s1 > 8 || L > 19999) 1478 /* Avoid confusion from exponents 1479 * so large that e might overflow. 1480 */ 1481 e = 19999; /* safe for 16 bit ints */ 1482 else 1483 e = (int)L; 1484 if (esign) 1485 e = -e; 1486 } 1487 else 1488 e = 0; 1489 } 1490 else 1491 s = s00; 1492 } 1493 if (!nd) { 1494 if (!nz && !nz0) 1495 s = s00; 1496 goto ret; 1497 } 1498 e1 = e -= nf; 1499 1500 /* Now we have nd0 digits, starting at s0, followed by a 1501 * decimal point, followed by nd-nd0 digits. The number we're 1502 * after is the integer represented by those digits times 1503 * 10**e */ 1504 1505 if (!nd0) 1506 nd0 = nd; 1507 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 1508 value(rv) = y; 1509 if (k > 9) 1510 value(rv) = tens[k - 9] * value(rv) + z; 1511 bd0 = 0; 1512 if (nd <= DBL_DIG 1513 #ifndef RND_PRODQUOT 1514 && FLT_ROUNDS == 1 1515 #endif 1516 ) { 1517 if (!e) 1518 goto ret; 1519 if (e > 0) { 1520 if (e <= Ten_pmax) { 1521 #ifdef VAX 1522 goto vax_ovfl_check; 1523 #else 1524 /* value(rv) = */ rounded_product(value(rv), 1525 tens[e]); 1526 goto ret; 1527 #endif 1528 } 1529 i = DBL_DIG - nd; 1530 if (e <= Ten_pmax + i) { 1531 /* A fancier test would sometimes let us do 1532 * this for larger i values. 1533 */ 1534 e -= i; 1535 value(rv) *= tens[i]; 1536 #ifdef VAX 1537 /* VAX exponent range is so narrow we must 1538 * worry about overflow here... 1539 */ 1540 vax_ovfl_check: 1541 word0(rv) -= P*Exp_msk1; 1542 /* value(rv) = */ rounded_product(value(rv), 1543 tens[e]); 1544 if ((word0(rv) & Exp_mask) 1545 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 1546 goto ovfl; 1547 word0(rv) += P*Exp_msk1; 1548 #else 1549 /* value(rv) = */ rounded_product(value(rv), 1550 tens[e]); 1551 #endif 1552 goto ret; 1553 } 1554 } 1555 #ifndef Inaccurate_Divide 1556 else if (e >= -Ten_pmax) { 1557 /* value(rv) = */ rounded_quotient(value(rv), 1558 tens[-e]); 1559 goto ret; 1560 } 1561 #endif 1562 } 1563 e1 += nd - k; 1564 1565 /* Get starting approximation = rv * 10**e1 */ 1566 1567 if (e1 > 0) { 1568 if ((i = e1 & 15) != 0) 1569 value(rv) *= tens[i]; 1570 if (e1 &= ~15) { 1571 if (e1 > DBL_MAX_10_EXP) { 1572 ovfl: 1573 errno = ERANGE; 1574 value(rv) = HUGE_VAL; 1575 if (bd0) 1576 goto retfree; 1577 goto ret; 1578 } 1579 if ((e1 = (unsigned int)e1 >> 4) != 0) { 1580 for(j = 0; e1 > 1; j++, 1581 e1 = (unsigned int)e1 >> 1) 1582 if (e1 & 1) 1583 value(rv) *= bigtens[j]; 1584 /* The last multiplication could overflow. */ 1585 word0(rv) -= P*Exp_msk1; 1586 value(rv) *= bigtens[j]; 1587 if ((z = word0(rv) & Exp_mask) 1588 > Exp_msk1*(DBL_MAX_EXP+Bias-P)) 1589 goto ovfl; 1590 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { 1591 /* set to largest number */ 1592 /* (Can't trust DBL_MAX) */ 1593 word0(rv) = Big0; 1594 word1(rv) = Big1; 1595 } 1596 else 1597 word0(rv) += P*Exp_msk1; 1598 } 1599 } 1600 } 1601 else if (e1 < 0) { 1602 e1 = -e1; 1603 if ((i = e1 & 15) != 0) 1604 value(rv) /= tens[i]; 1605 if (e1 &= ~15) { 1606 e1 = (unsigned int)e1 >> 4; 1607 if (e1 >= 1 << n_bigtens) 1608 goto undfl; 1609 for(j = 0; e1 > 1; j++, 1610 e1 = (unsigned int)e1 >> 1) 1611 if (e1 & 1) 1612 value(rv) *= tinytens[j]; 1613 /* The last multiplication could underflow. */ 1614 value(rv0) = value(rv); 1615 value(rv) *= tinytens[j]; 1616 if (!value(rv)) { 1617 value(rv) = 2.*value(rv0); 1618 value(rv) *= tinytens[j]; 1619 if (!value(rv)) { 1620 undfl: 1621 value(rv) = 0.; 1622 errno = ERANGE; 1623 if (bd0) 1624 goto retfree; 1625 goto ret; 1626 } 1627 word0(rv) = Tiny0; 1628 word1(rv) = Tiny1; 1629 /* The refinement below will clean 1630 * this approximation up. 1631 */ 1632 } 1633 } 1634 } 1635 1636 /* Now the hard part -- adjusting rv to the correct value.*/ 1637 1638 /* Put digits into bd: true value = bd * 10^e */ 1639 1640 bd0 = s2b(s0, nd0, nd, y); 1641 1642 for(;;) { 1643 bd = Balloc(bd0->k); 1644 Bcopy(bd, bd0); 1645 bb = d2b(value(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */ 1646 bs = i2b(1); 1647 1648 if (e >= 0) { 1649 bb2 = bb5 = 0; 1650 bd2 = bd5 = e; 1651 } 1652 else { 1653 bb2 = bb5 = -e; 1654 bd2 = bd5 = 0; 1655 } 1656 if (bbe >= 0) 1657 bb2 += bbe; 1658 else 1659 bd2 -= bbe; 1660 bs2 = bb2; 1661 #ifdef Sudden_Underflow 1662 #ifdef IBM 1663 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 1664 #else 1665 j = P + 1 - bbbits; 1666 #endif 1667 #else 1668 i = bbe + bbbits - 1; /* logb(rv) */ 1669 if (i < Emin) /* denormal */ 1670 j = bbe + (P-Emin); 1671 else 1672 j = P + 1 - bbbits; 1673 #endif 1674 bb2 += j; 1675 bd2 += j; 1676 i = bb2 < bd2 ? bb2 : bd2; 1677 if (i > bs2) 1678 i = bs2; 1679 if (i > 0) { 1680 bb2 -= i; 1681 bd2 -= i; 1682 bs2 -= i; 1683 } 1684 if (bb5 > 0) { 1685 bs = pow5mult(bs, bb5); 1686 bb1 = mult(bs, bb); 1687 Bfree(bb); 1688 bb = bb1; 1689 } 1690 if (bb2 > 0) 1691 bb = lshift(bb, bb2); 1692 if (bd5 > 0) 1693 bd = pow5mult(bd, bd5); 1694 if (bd2 > 0) 1695 bd = lshift(bd, bd2); 1696 if (bs2 > 0) 1697 bs = lshift(bs, bs2); 1698 delta = diff(bb, bd); 1699 dsign = delta->sign; 1700 delta->sign = 0; 1701 i = cmp(delta, bs); 1702 if (i < 0) { 1703 /* Error is less than half an ulp -- check for 1704 * special case of mantissa a power of two. 1705 */ 1706 if (dsign || word1(rv) || word0(rv) & Bndry_mask) 1707 break; 1708 delta = lshift(delta,Log2P); 1709 if (cmp(delta, bs) > 0) 1710 goto drop_down; 1711 break; 1712 } 1713 if (i == 0) { 1714 /* exactly half-way between */ 1715 if (dsign) { 1716 if ((word0(rv) & Bndry_mask1) == Bndry_mask1 1717 && word1(rv) == 0xffffffff) { 1718 /*boundary case -- increment exponent*/ 1719 word0(rv) = (word0(rv) & Exp_mask) 1720 + Exp_msk1 1721 #ifdef IBM 1722 | Exp_msk1 >> 4 1723 #endif 1724 ; 1725 word1(rv) = 0; 1726 break; 1727 } 1728 } 1729 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { 1730 drop_down: 1731 /* boundary case -- decrement exponent */ 1732 #ifdef Sudden_Underflow 1733 L = word0(rv) & Exp_mask; 1734 #ifdef IBM 1735 if (L < Exp_msk1) 1736 #else 1737 if (L <= Exp_msk1) 1738 #endif 1739 goto undfl; 1740 L -= Exp_msk1; 1741 #else 1742 L = (word0(rv) & Exp_mask) - Exp_msk1; 1743 #endif 1744 word0(rv) = L | Bndry_mask1; 1745 word1(rv) = 0xffffffff; 1746 #ifdef IBM 1747 goto cont; 1748 #else 1749 break; 1750 #endif 1751 } 1752 #ifndef ROUND_BIASED 1753 if (!(word1(rv) & LSB)) 1754 break; 1755 #endif 1756 if (dsign) 1757 value(rv) += ulp(value(rv)); 1758 #ifndef ROUND_BIASED 1759 else { 1760 value(rv) -= ulp(value(rv)); 1761 #ifndef Sudden_Underflow 1762 if (!value(rv)) 1763 goto undfl; 1764 #endif 1765 } 1766 #endif 1767 break; 1768 } 1769 if ((aadj = ratio(delta, bs)) <= 2.) { 1770 if (dsign) 1771 aadj = aadj1 = 1.; 1772 else if (word1(rv) || word0(rv) & Bndry_mask) { 1773 #ifndef Sudden_Underflow 1774 if (word1(rv) == Tiny1 && !word0(rv)) 1775 goto undfl; 1776 #endif 1777 aadj = 1.; 1778 aadj1 = -1.; 1779 } 1780 else { 1781 /* special case -- power of FLT_RADIX to be */ 1782 /* rounded down... */ 1783 1784 if (aadj < 2./FLT_RADIX) 1785 aadj = 1./FLT_RADIX; 1786 else 1787 aadj *= 0.5; 1788 aadj1 = -aadj; 1789 } 1790 } 1791 else { 1792 aadj *= 0.5; 1793 aadj1 = dsign ? aadj : -aadj; 1794 #ifdef Check_FLT_ROUNDS 1795 switch(FLT_ROUNDS) { 1796 case 2: /* towards +infinity */ 1797 aadj1 -= 0.5; 1798 break; 1799 case 0: /* towards 0 */ 1800 case 3: /* towards -infinity */ 1801 aadj1 += 0.5; 1802 } 1803 #else 1804 if (FLT_ROUNDS == 0) 1805 aadj1 += 0.5; 1806 #endif 1807 } 1808 y = word0(rv) & Exp_mask; 1809 1810 /* Check for overflow */ 1811 1812 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 1813 value(rv0) = value(rv); 1814 word0(rv) -= P*Exp_msk1; 1815 adj = aadj1 * ulp(value(rv)); 1816 value(rv) += adj; 1817 if ((word0(rv) & Exp_mask) >= 1818 Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 1819 if (word0(rv0) == Big0 && word1(rv0) == Big1) 1820 goto ovfl; 1821 word0(rv) = Big0; 1822 word1(rv) = Big1; 1823 goto cont; 1824 } 1825 else 1826 word0(rv) += P*Exp_msk1; 1827 } 1828 else { 1829 #ifdef Sudden_Underflow 1830 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 1831 value(rv0) = value(rv); 1832 word0(rv) += P*Exp_msk1; 1833 adj = aadj1 * ulp(value(rv)); 1834 value(rv) += adj; 1835 #ifdef IBM 1836 if ((word0(rv) & Exp_mask) < P*Exp_msk1) 1837 #else 1838 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) 1839 #endif 1840 { 1841 if (word0(rv0) == Tiny0 1842 && word1(rv0) == Tiny1) 1843 goto undfl; 1844 word0(rv) = Tiny0; 1845 word1(rv) = Tiny1; 1846 goto cont; 1847 } 1848 else 1849 word0(rv) -= P*Exp_msk1; 1850 } 1851 else { 1852 adj = aadj1 * ulp(value(rv)); 1853 value(rv) += adj; 1854 } 1855 #else 1856 /* Compute adj so that the IEEE rounding rules will 1857 * correctly round rv + adj in some half-way cases. 1858 * If rv * ulp(rv) is denormalized (i.e., 1859 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 1860 * trouble from bits lost to denormalization; 1861 * example: 1.2e-307 . 1862 */ 1863 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) { 1864 aadj1 = (double)(int)(aadj + 0.5); 1865 if (!dsign) 1866 aadj1 = -aadj1; 1867 } 1868 adj = aadj1 * ulp(value(rv)); 1869 value(rv) += adj; 1870 #endif 1871 } 1872 z = word0(rv) & Exp_mask; 1873 if (y == z) { 1874 /* Can we stop now? */ 1875 L = aadj; 1876 aadj -= L; 1877 /* The tolerances below are conservative. */ 1878 if (dsign || word1(rv) || word0(rv) & Bndry_mask) { 1879 if (aadj < .4999999 || aadj > .5000001) 1880 break; 1881 } 1882 else if (aadj < .4999999/FLT_RADIX) 1883 break; 1884 } 1885 cont: 1886 Bfree(bb); 1887 Bfree(bd); 1888 Bfree(bs); 1889 Bfree(delta); 1890 } 1891 retfree: 1892 Bfree(bb); 1893 Bfree(bd); 1894 Bfree(bs); 1895 Bfree(bd0); 1896 Bfree(delta); 1897 ret: 1898 if (se) 1899 /* LINTED interface specification */ 1900 *se = (char *)s; 1901 return sign ? -value(rv) : value(rv); 1902 } 1903 1904 static int 1905 quorem 1906 #ifdef KR_headers 1907 (b, S) Bigint *b, *S; 1908 #else 1909 (Bigint *b, Bigint *S) 1910 #endif 1911 { 1912 int n; 1913 Long borrow, y; 1914 ULong carry, q, ys; 1915 ULong *bx, *bxe, *sx, *sxe; 1916 #ifdef Pack_32 1917 Long z; 1918 ULong si, zs; 1919 #endif 1920 1921 if (b == BIGINT_INVALID || S == BIGINT_INVALID) 1922 return 0; 1923 1924 n = S->wds; 1925 #ifdef DEBUG 1926 /*debug*/ if (b->wds > n) 1927 /*debug*/ Bug("oversize b in quorem"); 1928 #endif 1929 if (b->wds < n) 1930 return 0; 1931 sx = S->x; 1932 sxe = sx + --n; 1933 bx = b->x; 1934 bxe = bx + n; 1935 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ 1936 #ifdef DEBUG 1937 /*debug*/ if (q > 9) 1938 /*debug*/ Bug("oversized quotient in quorem"); 1939 #endif 1940 if (q) { 1941 borrow = 0; 1942 carry = 0; 1943 do { 1944 #ifdef Pack_32 1945 si = *sx++; 1946 ys = (si & 0xffff) * q + carry; 1947 zs = (si >> 16) * q + (ys >> 16); 1948 carry = zs >> 16; 1949 y = (*bx & 0xffff) - (ys & 0xffff) + borrow; 1950 borrow = (ULong)y >> 16; 1951 Sign_Extend(borrow, y); 1952 z = (*bx >> 16) - (zs & 0xffff) + borrow; 1953 borrow = (ULong)z >> 16; 1954 Sign_Extend(borrow, z); 1955 Storeinc(bx, z, y); 1956 #else 1957 ys = *sx++ * q + carry; 1958 carry = ys >> 16; 1959 y = *bx - (ys & 0xffff) + borrow; 1960 borrow = y >> 16; 1961 Sign_Extend(borrow, y); 1962 *bx++ = y & 0xffff; 1963 #endif 1964 } 1965 while(sx <= sxe); 1966 if (!*bxe) { 1967 bx = b->x; 1968 while(--bxe > bx && !*bxe) 1969 --n; 1970 b->wds = n; 1971 } 1972 } 1973 if (cmp(b, S) >= 0) { 1974 q++; 1975 borrow = 0; 1976 carry = 0; 1977 bx = b->x; 1978 sx = S->x; 1979 do { 1980 #ifdef Pack_32 1981 si = *sx++; 1982 ys = (si & 0xffff) + carry; 1983 zs = (si >> 16) + (ys >> 16); 1984 carry = zs >> 16; 1985 y = (*bx & 0xffff) - (ys & 0xffff) + borrow; 1986 borrow = (ULong)y >> 16; 1987 Sign_Extend(borrow, y); 1988 z = (*bx >> 16) - (zs & 0xffff) + borrow; 1989 borrow = (ULong)z >> 16; 1990 Sign_Extend(borrow, z); 1991 Storeinc(bx, z, y); 1992 #else 1993 ys = *sx++ + carry; 1994 carry = ys >> 16; 1995 y = *bx - (ys & 0xffff) + borrow; 1996 borrow = y >> 16; 1997 Sign_Extend(borrow, y); 1998 *bx++ = y & 0xffff; 1999 #endif 2000 } 2001 while(sx <= sxe); 2002 bx = b->x; 2003 bxe = bx + n; 2004 if (!*bxe) { 2005 while(--bxe > bx && !*bxe) 2006 --n; 2007 b->wds = n; 2008 } 2009 } 2010 return q; 2011 } 2012 2013 /* freedtoa(s) must be used to free values s returned by dtoa 2014 * when MULTIPLE_THREADS is #defined. It should be used in all cases, 2015 * but for consistency with earlier versions of dtoa, it is optional 2016 * when MULTIPLE_THREADS is not defined. 2017 */ 2018 2019 void 2020 #ifdef KR_headers 2021 freedtoa(s) char *s; 2022 #else 2023 freedtoa(char *s) 2024 #endif 2025 { 2026 free(s); 2027 } 2028 2029 2030 2031 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. 2032 * 2033 * Inspired by "How to Print Floating-Point Numbers Accurately" by 2034 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101]. 2035 * 2036 * Modifications: 2037 * 1. Rather than iterating, we use a simple numeric overestimate 2038 * to determine k = floor(log10(d)). We scale relevant 2039 * quantities using O(log2(k)) rather than O(k) multiplications. 2040 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't 2041 * try to generate digits strictly left to right. Instead, we 2042 * compute with fewer bits and propagate the carry if necessary 2043 * when rounding the final digit up. This is often faster. 2044 * 3. Under the assumption that input will be rounded nearest, 2045 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. 2046 * That is, we allow equality in stopping tests when the 2047 * round-nearest rule will give the same floating-point value 2048 * as would satisfaction of the stopping test with strict 2049 * inequality. 2050 * 4. We remove common factors of powers of 2 from relevant 2051 * quantities. 2052 * 5. When converting floating-point integers less than 1e16, 2053 * we use floating-point arithmetic rather than resorting 2054 * to multiple-precision integers. 2055 * 6. When asked to produce fewer than 15 digits, we first try 2056 * to get by with floating-point arithmetic; we resort to 2057 * multiple-precision integer arithmetic only if we cannot 2058 * guarantee that the floating-point calculation has given 2059 * the correctly rounded result. For k requested digits and 2060 * "uniformly" distributed input, the probability is 2061 * something like 10^(k-15) that we must resort to the Long 2062 * calculation. 2063 */ 2064 2065 char * 2066 __dtoa 2067 #ifdef KR_headers 2068 (_d, mode, ndigits, decpt, sign, rve) 2069 double _d; int mode, ndigits, *decpt, *sign; char **rve; 2070 #else 2071 (double _d, int mode, int ndigits, int *decpt, int *sign, char **rve) 2072 #endif 2073 { 2074 /* Arguments ndigits, decpt, sign are similar to those 2075 of ecvt and fcvt; trailing zeros are suppressed from 2076 the returned string. If not null, *rve is set to point 2077 to the end of the return value. If d is +-Infinity or NaN, 2078 then *decpt is set to 9999. 2079 2080 mode: 2081 0 ==> shortest string that yields d when read in 2082 and rounded to nearest. 2083 1 ==> like 0, but with Steele & White stopping rule; 2084 e.g. with IEEE P754 arithmetic , mode 0 gives 2085 1e23 whereas mode 1 gives 9.999999999999999e22. 2086 2 ==> max(1,ndigits) significant digits. This gives a 2087 return value similar to that of ecvt, except 2088 that trailing zeros are suppressed. 2089 3 ==> through ndigits past the decimal point. This 2090 gives a return value similar to that from fcvt, 2091 except that trailing zeros are suppressed, and 2092 ndigits can be negative. 2093 4-9 should give the same return values as 2-3, i.e., 2094 4 <= mode <= 9 ==> same return as mode 2095 2 + (mode & 1). These modes are mainly for 2096 debugging; often they run slower but sometimes 2097 faster than modes 2-3. 2098 4,5,8,9 ==> left-to-right digit generation. 2099 6-9 ==> don't try fast floating-point estimate 2100 (if applicable). 2101 2102 Values of mode other than 0-9 are treated as mode 0. 2103 2104 Sufficient space is allocated to the return value 2105 to hold the suppressed trailing zeros. 2106 */ 2107 2108 int bbits, b2, b5, be, dig, i, ieps, ilim0, 2109 j, jj1, k, k0, k_check, leftright, m2, m5, s2, s5, 2110 try_quick; 2111 int ilim = 0, ilim1 = 0, spec_case = 0; /* pacify gcc */ 2112 Long L; 2113 #ifndef Sudden_Underflow 2114 int denorm; 2115 ULong x; 2116 #endif 2117 Bigint *b, *b1, *delta, *mhi, *S; 2118 Bigint *mlo = NULL; /* pacify gcc */ 2119 double ds; 2120 char *s, *s0; 2121 Bigint *result = NULL; 2122 int result_k = 0; 2123 _double d, d2, eps; 2124 2125 value(d) = _d; 2126 2127 if (word0(d) & Sign_bit) { 2128 /* set sign for everything, including 0's and NaNs */ 2129 *sign = 1; 2130 word0(d) &= ~Sign_bit; /* clear sign bit */ 2131 } 2132 else 2133 *sign = 0; 2134 2135 #if defined(IEEE_Arith) + defined(VAX) 2136 #ifdef IEEE_Arith 2137 if ((word0(d) & Exp_mask) == Exp_mask) 2138 #else 2139 if (word0(d) == 0x8000) 2140 #endif 2141 { 2142 /* Infinity or NaN */ 2143 *decpt = 9999; 2144 s = 2145 #ifdef IEEE_Arith 2146 !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" : 2147 #endif 2148 "NaN"; 2149 result = Balloc(strlen(s)+1); 2150 if (result == BIGINT_INVALID) 2151 return NULL; 2152 s0 = (char *)(void *)result; 2153 strcpy(s0, s); 2154 if (rve) 2155 *rve = 2156 #ifdef IEEE_Arith 2157 s0[3] ? s0 + 8 : 2158 #endif 2159 s0 + 3; 2160 return s0; 2161 } 2162 #endif 2163 #ifdef IBM 2164 value(d) += 0; /* normalize */ 2165 #endif 2166 if (!value(d)) { 2167 *decpt = 1; 2168 result = Balloc(2); 2169 if (result == BIGINT_INVALID) 2170 return NULL; 2171 s0 = (char *)(void *)result; 2172 strcpy(s0, "0"); 2173 if (rve) 2174 *rve = s0 + 1; 2175 return s0; 2176 } 2177 2178 b = d2b(value(d), &be, &bbits); 2179 #ifdef Sudden_Underflow 2180 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)); 2181 #else 2182 if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) != 0) { 2183 #endif 2184 value(d2) = value(d); 2185 word0(d2) &= Frac_mask1; 2186 word0(d2) |= Exp_11; 2187 #ifdef IBM 2188 if (j = 11 - hi0bits(word0(d2) & Frac_mask)) 2189 value(d2) /= 1 << j; 2190 #endif 2191 2192 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 2193 * log10(x) = log(x) / log(10) 2194 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) 2195 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) 2196 * 2197 * This suggests computing an approximation k to log10(d) by 2198 * 2199 * k = (i - Bias)*0.301029995663981 2200 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); 2201 * 2202 * We want k to be too large rather than too small. 2203 * The error in the first-order Taylor series approximation 2204 * is in our favor, so we just round up the constant enough 2205 * to compensate for any error in the multiplication of 2206 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, 2207 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, 2208 * adding 1e-13 to the constant term more than suffices. 2209 * Hence we adjust the constant term to 0.1760912590558. 2210 * (We could get a more accurate k by invoking log10, 2211 * but this is probably not worthwhile.) 2212 */ 2213 2214 i -= Bias; 2215 #ifdef IBM 2216 i <<= 2; 2217 i += j; 2218 #endif 2219 #ifndef Sudden_Underflow 2220 denorm = 0; 2221 } 2222 else { 2223 /* d is denormalized */ 2224 2225 i = bbits + be + (Bias + (P-1) - 1); 2226 x = i > 32 ? word0(d) << (64 - i) | word1(d) >> (i - 32) 2227 : word1(d) << (32 - i); 2228 value(d2) = x; 2229 word0(d2) -= 31*Exp_msk1; /* adjust exponent */ 2230 i -= (Bias + (P-1) - 1) + 1; 2231 denorm = 1; 2232 } 2233 #endif 2234 ds = (value(d2)-1.5)*0.289529654602168 + 0.1760912590558 + 2235 i*0.301029995663981; 2236 k = (int)ds; 2237 if (ds < 0. && ds != k) 2238 k--; /* want k = floor(ds) */ 2239 k_check = 1; 2240 if (k >= 0 && k <= Ten_pmax) { 2241 if (value(d) < tens[k]) 2242 k--; 2243 k_check = 0; 2244 } 2245 j = bbits - i - 1; 2246 if (j >= 0) { 2247 b2 = 0; 2248 s2 = j; 2249 } 2250 else { 2251 b2 = -j; 2252 s2 = 0; 2253 } 2254 if (k >= 0) { 2255 b5 = 0; 2256 s5 = k; 2257 s2 += k; 2258 } 2259 else { 2260 b2 -= k; 2261 b5 = -k; 2262 s5 = 0; 2263 } 2264 if (mode < 0 || mode > 9) 2265 mode = 0; 2266 try_quick = 1; 2267 if (mode > 5) { 2268 mode -= 4; 2269 try_quick = 0; 2270 } 2271 leftright = 1; 2272 switch(mode) { 2273 case 0: 2274 case 1: 2275 ilim = ilim1 = -1; 2276 i = 18; 2277 ndigits = 0; 2278 break; 2279 case 2: 2280 leftright = 0; 2281 /* FALLTHROUGH */ 2282 case 4: 2283 if (ndigits <= 0) 2284 ndigits = 1; 2285 ilim = ilim1 = i = ndigits; 2286 break; 2287 case 3: 2288 leftright = 0; 2289 /* FALLTHROUGH */ 2290 case 5: 2291 i = ndigits + k + 1; 2292 ilim = i; 2293 ilim1 = i - 1; 2294 if (i <= 0) 2295 i = 1; 2296 } 2297 j = sizeof(ULong); 2298 for(result_k = 0; (int)(sizeof(Bigint) - sizeof(ULong)) + j <= i; 2299 j <<= 1) result_k++; 2300 // this is really a ugly hack, the code uses Balloc 2301 // instead of malloc, but casts the result into a char* 2302 // it seems the only reason to do that is due to the 2303 // complicated way the block size need to be computed 2304 // buuurk.... 2305 result = Balloc(result_k); 2306 if (result == BIGINT_INVALID) { 2307 Bfree(b); 2308 return NULL; 2309 } 2310 s = s0 = (char *)(void *)result; 2311 2312 if (ilim >= 0 && ilim <= Quick_max && try_quick) { 2313 2314 /* Try to get by with floating-point arithmetic. */ 2315 2316 i = 0; 2317 value(d2) = value(d); 2318 k0 = k; 2319 ilim0 = ilim; 2320 ieps = 2; /* conservative */ 2321 if (k > 0) { 2322 ds = tens[k&0xf]; 2323 j = (unsigned int)k >> 4; 2324 if (j & Bletch) { 2325 /* prevent overflows */ 2326 j &= Bletch - 1; 2327 value(d) /= bigtens[n_bigtens-1]; 2328 ieps++; 2329 } 2330 for(; j; j = (unsigned int)j >> 1, i++) 2331 if (j & 1) { 2332 ieps++; 2333 ds *= bigtens[i]; 2334 } 2335 value(d) /= ds; 2336 } 2337 else if ((jj1 = -k) != 0) { 2338 value(d) *= tens[jj1 & 0xf]; 2339 for(j = (unsigned int)jj1 >> 4; j; 2340 j = (unsigned int)j >> 1, i++) 2341 if (j & 1) { 2342 ieps++; 2343 value(d) *= bigtens[i]; 2344 } 2345 } 2346 if (k_check && value(d) < 1. && ilim > 0) { 2347 if (ilim1 <= 0) 2348 goto fast_failed; 2349 ilim = ilim1; 2350 k--; 2351 value(d) *= 10.; 2352 ieps++; 2353 } 2354 value(eps) = ieps*value(d) + 7.; 2355 word0(eps) -= (P-1)*Exp_msk1; 2356 if (ilim == 0) { 2357 S = mhi = 0; 2358 value(d) -= 5.; 2359 if (value(d) > value(eps)) 2360 goto one_digit; 2361 if (value(d) < -value(eps)) 2362 goto no_digits; 2363 goto fast_failed; 2364 } 2365 #ifndef No_leftright 2366 if (leftright) { 2367 /* Use Steele & White method of only 2368 * generating digits needed. 2369 */ 2370 value(eps) = 0.5/tens[ilim-1] - value(eps); 2371 for(i = 0;;) { 2372 L = value(d); 2373 value(d) -= L; 2374 *s++ = '0' + (int)L; 2375 if (value(d) < value(eps)) 2376 goto ret1; 2377 if (1. - value(d) < value(eps)) 2378 goto bump_up; 2379 if (++i >= ilim) 2380 break; 2381 value(eps) *= 10.; 2382 value(d) *= 10.; 2383 } 2384 } 2385 else { 2386 #endif 2387 /* Generate ilim digits, then fix them up. */ 2388 value(eps) *= tens[ilim-1]; 2389 for(i = 1;; i++, value(d) *= 10.) { 2390 L = value(d); 2391 value(d) -= L; 2392 *s++ = '0' + (int)L; 2393 if (i == ilim) { 2394 if (value(d) > 0.5 + value(eps)) 2395 goto bump_up; 2396 else if (value(d) < 0.5 - value(eps)) { 2397 while(*--s == '0'); 2398 s++; 2399 goto ret1; 2400 } 2401 break; 2402 } 2403 } 2404 #ifndef No_leftright 2405 } 2406 #endif 2407 fast_failed: 2408 s = s0; 2409 value(d) = value(d2); 2410 k = k0; 2411 ilim = ilim0; 2412 } 2413 2414 /* Do we have a "small" integer? */ 2415 2416 if (be >= 0 && k <= Int_max) { 2417 /* Yes. */ 2418 ds = tens[k]; 2419 if (ndigits < 0 && ilim <= 0) { 2420 S = mhi = 0; 2421 if (ilim < 0 || value(d) <= 5*ds) 2422 goto no_digits; 2423 goto one_digit; 2424 } 2425 for(i = 1;; i++) { 2426 L = value(d) / ds; 2427 value(d) -= L*ds; 2428 #ifdef Check_FLT_ROUNDS 2429 /* If FLT_ROUNDS == 2, L will usually be high by 1 */ 2430 if (value(d) < 0) { 2431 L--; 2432 value(d) += ds; 2433 } 2434 #endif 2435 *s++ = '0' + (int)L; 2436 if (i == ilim) { 2437 value(d) += value(d); 2438 if (value(d) > ds || (value(d) == ds && L & 1)) { 2439 bump_up: 2440 while(*--s == '9') 2441 if (s == s0) { 2442 k++; 2443 *s = '0'; 2444 break; 2445 } 2446 ++*s++; 2447 } 2448 break; 2449 } 2450 if (!(value(d) *= 10.)) 2451 break; 2452 } 2453 goto ret1; 2454 } 2455 2456 m2 = b2; 2457 m5 = b5; 2458 mhi = mlo = 0; 2459 if (leftright) { 2460 if (mode < 2) { 2461 i = 2462 #ifndef Sudden_Underflow 2463 denorm ? be + (Bias + (P-1) - 1 + 1) : 2464 #endif 2465 #ifdef IBM 2466 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3); 2467 #else 2468 1 + P - bbits; 2469 #endif 2470 } 2471 else { 2472 j = ilim - 1; 2473 if (m5 >= j) 2474 m5 -= j; 2475 else { 2476 s5 += j -= m5; 2477 b5 += j; 2478 m5 = 0; 2479 } 2480 if ((i = ilim) < 0) { 2481 m2 -= i; 2482 i = 0; 2483 } 2484 } 2485 b2 += i; 2486 s2 += i; 2487 mhi = i2b(1); 2488 } 2489 if (m2 > 0 && s2 > 0) { 2490 i = m2 < s2 ? m2 : s2; 2491 b2 -= i; 2492 m2 -= i; 2493 s2 -= i; 2494 } 2495 if (b5 > 0) { 2496 if (leftright) { 2497 if (m5 > 0) { 2498 mhi = pow5mult(mhi, m5); 2499 b1 = mult(mhi, b); 2500 Bfree(b); 2501 b = b1; 2502 } 2503 if ((j = b5 - m5) != 0) 2504 b = pow5mult(b, j); 2505 } 2506 else 2507 b = pow5mult(b, b5); 2508 } 2509 S = i2b(1); 2510 if (s5 > 0) 2511 S = pow5mult(S, s5); 2512 2513 /* Check for special case that d is a normalized power of 2. */ 2514 2515 if (mode < 2) { 2516 if (!word1(d) && !(word0(d) & Bndry_mask) 2517 #ifndef Sudden_Underflow 2518 && word0(d) & Exp_mask 2519 #endif 2520 ) { 2521 /* The special case */ 2522 b2 += Log2P; 2523 s2 += Log2P; 2524 spec_case = 1; 2525 } 2526 else 2527 spec_case = 0; 2528 } 2529 2530 /* Arrange for convenient computation of quotients: 2531 * shift left if necessary so divisor has 4 leading 0 bits. 2532 * 2533 * Perhaps we should just compute leading 28 bits of S once 2534 * and for all and pass them and a shift to quorem, so it 2535 * can do shifts and ors to compute the numerator for q. 2536 */ 2537 if (S == BIGINT_INVALID) { 2538 i = 0; 2539 } else { 2540 #ifdef Pack_32 2541 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) != 0) 2542 i = 32 - i; 2543 #else 2544 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) 2545 i = 16 - i; 2546 #endif 2547 } 2548 2549 if (i > 4) { 2550 i -= 4; 2551 b2 += i; 2552 m2 += i; 2553 s2 += i; 2554 } 2555 else if (i < 4) { 2556 i += 28; 2557 b2 += i; 2558 m2 += i; 2559 s2 += i; 2560 } 2561 if (b2 > 0) 2562 b = lshift(b, b2); 2563 if (s2 > 0) 2564 S = lshift(S, s2); 2565 if (k_check) { 2566 if (cmp(b,S) < 0) { 2567 k--; 2568 b = multadd(b, 10, 0); /* we botched the k estimate */ 2569 if (leftright) 2570 mhi = multadd(mhi, 10, 0); 2571 ilim = ilim1; 2572 } 2573 } 2574 if (ilim <= 0 && mode > 2) { 2575 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) { 2576 /* no digits, fcvt style */ 2577 no_digits: 2578 k = -1 - ndigits; 2579 goto ret; 2580 } 2581 one_digit: 2582 *s++ = '1'; 2583 k++; 2584 goto ret; 2585 } 2586 if (leftright) { 2587 if (m2 > 0) 2588 mhi = lshift(mhi, m2); 2589 2590 /* Compute mlo -- check for special case 2591 * that d is a normalized power of 2. 2592 */ 2593 2594 mlo = mhi; 2595 if (spec_case) { 2596 mhi = Balloc(mhi->k); 2597 Bcopy(mhi, mlo); 2598 mhi = lshift(mhi, Log2P); 2599 } 2600 2601 for(i = 1;;i++) { 2602 dig = quorem(b,S) + '0'; 2603 /* Do we yet have the shortest decimal string 2604 * that will round to d? 2605 */ 2606 j = cmp(b, mlo); 2607 delta = diff(S, mhi); 2608 jj1 = delta->sign ? 1 : cmp(b, delta); 2609 Bfree(delta); 2610 #ifndef ROUND_BIASED 2611 if (jj1 == 0 && !mode && !(word1(d) & 1)) { 2612 if (dig == '9') 2613 goto round_9_up; 2614 if (j > 0) 2615 dig++; 2616 *s++ = dig; 2617 goto ret; 2618 } 2619 #endif 2620 if (j < 0 || (j == 0 && !mode 2621 #ifndef ROUND_BIASED 2622 && !(word1(d) & 1) 2623 #endif 2624 )) { 2625 if (jj1 > 0) { 2626 b = lshift(b, 1); 2627 jj1 = cmp(b, S); 2628 if ((jj1 > 0 || (jj1 == 0 && dig & 1)) 2629 && dig++ == '9') 2630 goto round_9_up; 2631 } 2632 *s++ = dig; 2633 goto ret; 2634 } 2635 if (jj1 > 0) { 2636 if (dig == '9') { /* possible if i == 1 */ 2637 round_9_up: 2638 *s++ = '9'; 2639 goto roundoff; 2640 } 2641 *s++ = dig + 1; 2642 goto ret; 2643 } 2644 *s++ = dig; 2645 if (i == ilim) 2646 break; 2647 b = multadd(b, 10, 0); 2648 if (mlo == mhi) 2649 mlo = mhi = multadd(mhi, 10, 0); 2650 else { 2651 mlo = multadd(mlo, 10, 0); 2652 mhi = multadd(mhi, 10, 0); 2653 } 2654 } 2655 } 2656 else 2657 for(i = 1;; i++) { 2658 *s++ = dig = quorem(b,S) + '0'; 2659 if (i >= ilim) 2660 break; 2661 b = multadd(b, 10, 0); 2662 } 2663 2664 /* Round off last digit */ 2665 2666 b = lshift(b, 1); 2667 j = cmp(b, S); 2668 if (j > 0 || (j == 0 && dig & 1)) { 2669 roundoff: 2670 while(*--s == '9') 2671 if (s == s0) { 2672 k++; 2673 *s++ = '1'; 2674 goto ret; 2675 } 2676 ++*s++; 2677 } 2678 else { 2679 while(*--s == '0'); 2680 s++; 2681 } 2682 ret: 2683 Bfree(S); 2684 if (mhi) { 2685 if (mlo && mlo != mhi) 2686 Bfree(mlo); 2687 Bfree(mhi); 2688 } 2689 ret1: 2690 Bfree(b); 2691 if (s == s0) { /* don't return empty string */ 2692 *s++ = '0'; 2693 k = 0; 2694 } 2695 *s = 0; 2696 *decpt = k + 1; 2697 if (rve) 2698 *rve = s; 2699 return s0; 2700 } 2701 #ifdef __cplusplus 2702 } 2703 #endif 2704