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