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