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