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