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