dtoa.c revision 81612e877870c52bae7c590076eec642b9803138
1/**************************************************************** 2 * 3 * The author of this software is David M. Gay. 4 * 5 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies. 6 * 7 * Permission to use, copy, modify, and distribute this software for any 8 * purpose without fee is hereby granted, provided that this entire notice 9 * is included in all copies of any software which is or includes a copy 10 * or modification of this software and in all copies of the supporting 11 * documentation for such software. 12 * 13 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED 14 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY 15 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY 16 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. 17 * 18 ***************************************************************/ 19 20/**************************************************************** 21 * This is dtoa.c by David M. Gay, downloaded from 22 * http://www.netlib.org/fp/dtoa.c on April 15, 2009 and modified for 23 * inclusion into the Python core by Mark E. T. Dickinson and Eric V. Smith. 24 * 25 * Please remember to check http://www.netlib.org/fp regularly (and especially 26 * before any Python release) for bugfixes and updates. 27 * 28 * The major modifications from Gay's original code are as follows: 29 * 30 * 0. The original code has been specialized to Python's needs by removing 31 * many of the #ifdef'd sections. In particular, code to support VAX and 32 * IBM floating-point formats, hex NaNs, hex floats, locale-aware 33 * treatment of the decimal point, and setting of the inexact flag have 34 * been removed. 35 * 36 * 1. We use PyMem_Malloc and PyMem_Free in place of malloc and free. 37 * 38 * 2. The public functions strtod, dtoa and freedtoa all now have 39 * a _Py_dg_ prefix. 40 * 41 * 3. Instead of assuming that PyMem_Malloc always succeeds, we thread 42 * PyMem_Malloc failures through the code. The functions 43 * 44 * Balloc, multadd, s2b, i2b, mult, pow5mult, lshift, diff, d2b 45 * 46 * of return type *Bigint all return NULL to indicate a malloc failure. 47 * Similarly, rv_alloc and nrv_alloc (return type char *) return NULL on 48 * failure. bigcomp now has return type int (it used to be void) and 49 * returns -1 on failure and 0 otherwise. _Py_dg_dtoa returns NULL 50 * on failure. _Py_dg_strtod indicates failure due to malloc failure 51 * by returning -1.0, setting errno=ENOMEM and *se to s00. 52 * 53 * 4. The static variable dtoa_result has been removed. Callers of 54 * _Py_dg_dtoa are expected to call _Py_dg_freedtoa to free 55 * the memory allocated by _Py_dg_dtoa. 56 * 57 * 5. The code has been reformatted to better fit with Python's 58 * C style guide (PEP 7). 59 * 60 * 6. A bug in the memory allocation has been fixed: to avoid FREEing memory 61 * that hasn't been MALLOC'ed, private_mem should only be used when k <= 62 * Kmax. 63 * 64 * 7. _Py_dg_strtod has been modified so that it doesn't accept strings with 65 * leading whitespace. 66 * 67 ***************************************************************/ 68 69/* Please send bug reports for the original dtoa.c code to David M. Gay (dmg 70 * at acm dot org, with " at " changed at "@" and " dot " changed to "."). 71 * Please report bugs for this modified version using the Python issue tracker 72 * (http://bugs.python.org). */ 73 74/* On a machine with IEEE extended-precision registers, it is 75 * necessary to specify double-precision (53-bit) rounding precision 76 * before invoking strtod or dtoa. If the machine uses (the equivalent 77 * of) Intel 80x87 arithmetic, the call 78 * _control87(PC_53, MCW_PC); 79 * does this with many compilers. Whether this or another call is 80 * appropriate depends on the compiler; for this to work, it may be 81 * necessary to #include "float.h" or another system-dependent header 82 * file. 83 */ 84 85/* strtod for IEEE-, VAX-, and IBM-arithmetic machines. 86 * 87 * This strtod returns a nearest machine number to the input decimal 88 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are 89 * broken by the IEEE round-even rule. Otherwise ties are broken by 90 * biased rounding (add half and chop). 91 * 92 * Inspired loosely by William D. Clinger's paper "How to Read Floating 93 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. 94 * 95 * Modifications: 96 * 97 * 1. We only require IEEE, IBM, or VAX double-precision 98 * arithmetic (not IEEE double-extended). 99 * 2. We get by with floating-point arithmetic in a case that 100 * Clinger missed -- when we're computing d * 10^n 101 * for a small integer d and the integer n is not too 102 * much larger than 22 (the maximum integer k for which 103 * we can represent 10^k exactly), we may be able to 104 * compute (d*10^k) * 10^(e-k) with just one roundoff. 105 * 3. Rather than a bit-at-a-time adjustment of the binary 106 * result in the hard case, we use floating-point 107 * arithmetic to determine the adjustment to within 108 * one bit; only in really hard cases do we need to 109 * compute a second residual. 110 * 4. Because of 3., we don't need a large table of powers of 10 111 * for ten-to-e (just some small tables, e.g. of 10^k 112 * for 0 <= k <= 22). 113 */ 114 115/* Linking of Python's #defines to Gay's #defines starts here. */ 116 117#include "Python.h" 118 119/* if PY_NO_SHORT_FLOAT_REPR is defined, then don't even try to compile 120 the following code */ 121#ifndef PY_NO_SHORT_FLOAT_REPR 122 123#include "float.h" 124 125#define MALLOC PyMem_Malloc 126#define FREE PyMem_Free 127 128/* This code should also work for ARM mixed-endian format on little-endian 129 machines, where doubles have byte order 45670123 (in increasing address 130 order, 0 being the least significant byte). */ 131#ifdef DOUBLE_IS_LITTLE_ENDIAN_IEEE754 132# define IEEE_8087 133#endif 134#if defined(DOUBLE_IS_BIG_ENDIAN_IEEE754) || \ 135 defined(DOUBLE_IS_ARM_MIXED_ENDIAN_IEEE754) 136# define IEEE_MC68k 137#endif 138#if defined(IEEE_8087) + defined(IEEE_MC68k) != 1 139#error "Exactly one of IEEE_8087 or IEEE_MC68k should be defined." 140#endif 141 142/* The code below assumes that the endianness of integers matches the 143 endianness of the two 32-bit words of a double. Check this. */ 144#if defined(WORDS_BIGENDIAN) && (defined(DOUBLE_IS_LITTLE_ENDIAN_IEEE754) || \ 145 defined(DOUBLE_IS_ARM_MIXED_ENDIAN_IEEE754)) 146#error "doubles and ints have incompatible endianness" 147#endif 148 149#if !defined(WORDS_BIGENDIAN) && defined(DOUBLE_IS_BIG_ENDIAN_IEEE754) 150#error "doubles and ints have incompatible endianness" 151#endif 152 153 154#if defined(HAVE_UINT32_T) && defined(HAVE_INT32_T) 155typedef PY_UINT32_T ULong; 156typedef PY_INT32_T Long; 157#else 158#error "Failed to find an exact-width 32-bit integer type" 159#endif 160 161#if defined(HAVE_UINT64_T) 162#define ULLong PY_UINT64_T 163#else 164#undef ULLong 165#endif 166 167#undef DEBUG 168#ifdef Py_DEBUG 169#define DEBUG 170#endif 171 172/* End Python #define linking */ 173 174#ifdef DEBUG 175#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} 176#endif 177 178#ifndef PRIVATE_MEM 179#define PRIVATE_MEM 2304 180#endif 181#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double)) 182static double private_mem[PRIVATE_mem], *pmem_next = private_mem; 183 184#ifdef __cplusplus 185extern "C" { 186#endif 187 188typedef union { double d; ULong L[2]; } U; 189 190#ifdef IEEE_8087 191#define word0(x) (x)->L[1] 192#define word1(x) (x)->L[0] 193#else 194#define word0(x) (x)->L[0] 195#define word1(x) (x)->L[1] 196#endif 197#define dval(x) (x)->d 198 199#ifndef STRTOD_DIGLIM 200#define STRTOD_DIGLIM 40 201#endif 202 203/* maximum permitted exponent value for strtod; exponents larger than 204 MAX_ABS_EXP in absolute value get truncated to +-MAX_ABS_EXP. MAX_ABS_EXP 205 should fit into an int. */ 206#ifndef MAX_ABS_EXP 207#define MAX_ABS_EXP 19999U 208#endif 209 210/* The following definition of Storeinc is appropriate for MIPS processors. 211 * An alternative that might be better on some machines is 212 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) 213 */ 214#if defined(IEEE_8087) 215#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \ 216 ((unsigned short *)a)[0] = (unsigned short)c, a++) 217#else 218#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \ 219 ((unsigned short *)a)[1] = (unsigned short)c, a++) 220#endif 221 222/* #define P DBL_MANT_DIG */ 223/* Ten_pmax = floor(P*log(2)/log(5)) */ 224/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ 225/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ 226/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ 227 228#define Exp_shift 20 229#define Exp_shift1 20 230#define Exp_msk1 0x100000 231#define Exp_msk11 0x100000 232#define Exp_mask 0x7ff00000 233#define P 53 234#define Nbits 53 235#define Bias 1023 236#define Emax 1023 237#define Emin (-1022) 238#define Exp_1 0x3ff00000 239#define Exp_11 0x3ff00000 240#define Ebits 11 241#define Frac_mask 0xfffff 242#define Frac_mask1 0xfffff 243#define Ten_pmax 22 244#define Bletch 0x10 245#define Bndry_mask 0xfffff 246#define Bndry_mask1 0xfffff 247#define LSB 1 248#define Sign_bit 0x80000000 249#define Log2P 1 250#define Tiny0 0 251#define Tiny1 1 252#define Quick_max 14 253#define Int_max 14 254 255#ifndef Flt_Rounds 256#ifdef FLT_ROUNDS 257#define Flt_Rounds FLT_ROUNDS 258#else 259#define Flt_Rounds 1 260#endif 261#endif /*Flt_Rounds*/ 262 263#define Rounding Flt_Rounds 264 265#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) 266#define Big1 0xffffffff 267 268/* struct BCinfo is used to pass information from _Py_dg_strtod to bigcomp */ 269 270typedef struct BCinfo BCinfo; 271struct 272BCinfo { 273 int dp0, dp1, dplen, dsign, e0, nd, nd0, scale; 274}; 275 276#define FFFFFFFF 0xffffffffUL 277 278#define Kmax 7 279 280/* struct Bigint is used to represent arbitrary-precision integers. These 281 integers are stored in sign-magnitude format, with the magnitude stored as 282 an array of base 2**32 digits. Bigints are always normalized: if x is a 283 Bigint then x->wds >= 1, and either x->wds == 1 or x[wds-1] is nonzero. 284 285 The Bigint fields are as follows: 286 287 - next is a header used by Balloc and Bfree to keep track of lists 288 of freed Bigints; it's also used for the linked list of 289 powers of 5 of the form 5**2**i used by pow5mult. 290 - k indicates which pool this Bigint was allocated from 291 - maxwds is the maximum number of words space was allocated for 292 (usually maxwds == 2**k) 293 - sign is 1 for negative Bigints, 0 for positive. The sign is unused 294 (ignored on inputs, set to 0 on outputs) in almost all operations 295 involving Bigints: a notable exception is the diff function, which 296 ignores signs on inputs but sets the sign of the output correctly. 297 - wds is the actual number of significant words 298 - x contains the vector of words (digits) for this Bigint, from least 299 significant (x[0]) to most significant (x[wds-1]). 300*/ 301 302struct 303Bigint { 304 struct Bigint *next; 305 int k, maxwds, sign, wds; 306 ULong x[1]; 307}; 308 309typedef struct Bigint Bigint; 310 311/* Memory management: memory is allocated from, and returned to, Kmax+1 pools 312 of memory, where pool k (0 <= k <= Kmax) is for Bigints b with b->maxwds == 313 1 << k. These pools are maintained as linked lists, with freelist[k] 314 pointing to the head of the list for pool k. 315 316 On allocation, if there's no free slot in the appropriate pool, MALLOC is 317 called to get more memory. This memory is not returned to the system until 318 Python quits. There's also a private memory pool that's allocated from 319 in preference to using MALLOC. 320 321 For Bigints with more than (1 << Kmax) digits (which implies at least 1233 322 decimal digits), memory is directly allocated using MALLOC, and freed using 323 FREE. 324 325 XXX: it would be easy to bypass this memory-management system and 326 translate each call to Balloc into a call to PyMem_Malloc, and each 327 Bfree to PyMem_Free. Investigate whether this has any significant 328 performance on impact. */ 329 330static Bigint *freelist[Kmax+1]; 331 332/* Allocate space for a Bigint with up to 1<<k digits */ 333 334static Bigint * 335Balloc(int k) 336{ 337 int x; 338 Bigint *rv; 339 unsigned int len; 340 341 if (k <= Kmax && (rv = freelist[k])) 342 freelist[k] = rv->next; 343 else { 344 x = 1 << k; 345 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1) 346 /sizeof(double); 347 if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) { 348 rv = (Bigint*)pmem_next; 349 pmem_next += len; 350 } 351 else { 352 rv = (Bigint*)MALLOC(len*sizeof(double)); 353 if (rv == NULL) 354 return NULL; 355 } 356 rv->k = k; 357 rv->maxwds = x; 358 } 359 rv->sign = rv->wds = 0; 360 return rv; 361} 362 363/* Free a Bigint allocated with Balloc */ 364 365static void 366Bfree(Bigint *v) 367{ 368 if (v) { 369 if (v->k > Kmax) 370 FREE((void*)v); 371 else { 372 v->next = freelist[v->k]; 373 freelist[v->k] = v; 374 } 375 } 376} 377 378#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \ 379 y->wds*sizeof(Long) + 2*sizeof(int)) 380 381/* Multiply a Bigint b by m and add a. Either modifies b in place and returns 382 a pointer to the modified b, or Bfrees b and returns a pointer to a copy. 383 On failure, return NULL. In this case, b will have been already freed. */ 384 385static Bigint * 386multadd(Bigint *b, int m, int a) /* multiply by m and add a */ 387{ 388 int i, wds; 389#ifdef ULLong 390 ULong *x; 391 ULLong carry, y; 392#else 393 ULong carry, *x, y; 394 ULong xi, z; 395#endif 396 Bigint *b1; 397 398 wds = b->wds; 399 x = b->x; 400 i = 0; 401 carry = a; 402 do { 403#ifdef ULLong 404 y = *x * (ULLong)m + carry; 405 carry = y >> 32; 406 *x++ = (ULong)(y & FFFFFFFF); 407#else 408 xi = *x; 409 y = (xi & 0xffff) * m + carry; 410 z = (xi >> 16) * m + (y >> 16); 411 carry = z >> 16; 412 *x++ = (z << 16) + (y & 0xffff); 413#endif 414 } 415 while(++i < wds); 416 if (carry) { 417 if (wds >= b->maxwds) { 418 b1 = Balloc(b->k+1); 419 if (b1 == NULL){ 420 Bfree(b); 421 return NULL; 422 } 423 Bcopy(b1, b); 424 Bfree(b); 425 b = b1; 426 } 427 b->x[wds++] = (ULong)carry; 428 b->wds = wds; 429 } 430 return b; 431} 432 433/* convert a string s containing nd decimal digits (possibly containing a 434 decimal separator at position nd0, which is ignored) to a Bigint. This 435 function carries on where the parsing code in _Py_dg_strtod leaves off: on 436 entry, y9 contains the result of converting the first 9 digits. Returns 437 NULL on failure. */ 438 439static Bigint * 440s2b(const char *s, int nd0, int nd, ULong y9, int dplen) 441{ 442 Bigint *b; 443 int i, k; 444 Long x, y; 445 446 x = (nd + 8) / 9; 447 for(k = 0, y = 1; x > y; y <<= 1, k++) ; 448 b = Balloc(k); 449 if (b == NULL) 450 return NULL; 451 b->x[0] = y9; 452 b->wds = 1; 453 454 i = 9; 455 if (9 < nd0) { 456 s += 9; 457 do { 458 b = multadd(b, 10, *s++ - '0'); 459 if (b == NULL) 460 return NULL; 461 } while(++i < nd0); 462 s += dplen; 463 } 464 else 465 s += dplen + 9; 466 for(; i < nd; i++) { 467 b = multadd(b, 10, *s++ - '0'); 468 if (b == NULL) 469 return NULL; 470 } 471 return b; 472} 473 474/* count leading 0 bits in the 32-bit integer x. */ 475 476static int 477hi0bits(ULong x) 478{ 479 int k = 0; 480 481 if (!(x & 0xffff0000)) { 482 k = 16; 483 x <<= 16; 484 } 485 if (!(x & 0xff000000)) { 486 k += 8; 487 x <<= 8; 488 } 489 if (!(x & 0xf0000000)) { 490 k += 4; 491 x <<= 4; 492 } 493 if (!(x & 0xc0000000)) { 494 k += 2; 495 x <<= 2; 496 } 497 if (!(x & 0x80000000)) { 498 k++; 499 if (!(x & 0x40000000)) 500 return 32; 501 } 502 return k; 503} 504 505/* count trailing 0 bits in the 32-bit integer y, and shift y right by that 506 number of bits. */ 507 508static int 509lo0bits(ULong *y) 510{ 511 int k; 512 ULong x = *y; 513 514 if (x & 7) { 515 if (x & 1) 516 return 0; 517 if (x & 2) { 518 *y = x >> 1; 519 return 1; 520 } 521 *y = x >> 2; 522 return 2; 523 } 524 k = 0; 525 if (!(x & 0xffff)) { 526 k = 16; 527 x >>= 16; 528 } 529 if (!(x & 0xff)) { 530 k += 8; 531 x >>= 8; 532 } 533 if (!(x & 0xf)) { 534 k += 4; 535 x >>= 4; 536 } 537 if (!(x & 0x3)) { 538 k += 2; 539 x >>= 2; 540 } 541 if (!(x & 1)) { 542 k++; 543 x >>= 1; 544 if (!x) 545 return 32; 546 } 547 *y = x; 548 return k; 549} 550 551/* convert a small nonnegative integer to a Bigint */ 552 553static Bigint * 554i2b(int i) 555{ 556 Bigint *b; 557 558 b = Balloc(1); 559 if (b == NULL) 560 return NULL; 561 b->x[0] = i; 562 b->wds = 1; 563 return b; 564} 565 566/* multiply two Bigints. Returns a new Bigint, or NULL on failure. Ignores 567 the signs of a and b. */ 568 569static Bigint * 570mult(Bigint *a, Bigint *b) 571{ 572 Bigint *c; 573 int k, wa, wb, wc; 574 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; 575 ULong y; 576#ifdef ULLong 577 ULLong carry, z; 578#else 579 ULong carry, z; 580 ULong z2; 581#endif 582 583 if (a->wds < b->wds) { 584 c = a; 585 a = b; 586 b = c; 587 } 588 k = a->k; 589 wa = a->wds; 590 wb = b->wds; 591 wc = wa + wb; 592 if (wc > a->maxwds) 593 k++; 594 c = Balloc(k); 595 if (c == NULL) 596 return NULL; 597 for(x = c->x, xa = x + wc; x < xa; x++) 598 *x = 0; 599 xa = a->x; 600 xae = xa + wa; 601 xb = b->x; 602 xbe = xb + wb; 603 xc0 = c->x; 604#ifdef ULLong 605 for(; xb < xbe; xc0++) { 606 if ((y = *xb++)) { 607 x = xa; 608 xc = xc0; 609 carry = 0; 610 do { 611 z = *x++ * (ULLong)y + *xc + carry; 612 carry = z >> 32; 613 *xc++ = (ULong)(z & FFFFFFFF); 614 } 615 while(x < xae); 616 *xc = (ULong)carry; 617 } 618 } 619#else 620 for(; xb < xbe; xb++, xc0++) { 621 if (y = *xb & 0xffff) { 622 x = xa; 623 xc = xc0; 624 carry = 0; 625 do { 626 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 627 carry = z >> 16; 628 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 629 carry = z2 >> 16; 630 Storeinc(xc, z2, z); 631 } 632 while(x < xae); 633 *xc = carry; 634 } 635 if (y = *xb >> 16) { 636 x = xa; 637 xc = xc0; 638 carry = 0; 639 z2 = *xc; 640 do { 641 z = (*x & 0xffff) * y + (*xc >> 16) + carry; 642 carry = z >> 16; 643 Storeinc(xc, z, z2); 644 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 645 carry = z2 >> 16; 646 } 647 while(x < xae); 648 *xc = z2; 649 } 650 } 651#endif 652 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ; 653 c->wds = wc; 654 return c; 655} 656 657/* p5s is a linked list of powers of 5 of the form 5**(2**i), i >= 2 */ 658 659static Bigint *p5s; 660 661/* multiply the Bigint b by 5**k. Returns a pointer to the result, or NULL on 662 failure; if the returned pointer is distinct from b then the original 663 Bigint b will have been Bfree'd. Ignores the sign of b. */ 664 665static Bigint * 666pow5mult(Bigint *b, int k) 667{ 668 Bigint *b1, *p5, *p51; 669 int i; 670 static int p05[3] = { 5, 25, 125 }; 671 672 if ((i = k & 3)) { 673 b = multadd(b, p05[i-1], 0); 674 if (b == NULL) 675 return NULL; 676 } 677 678 if (!(k >>= 2)) 679 return b; 680 p5 = p5s; 681 if (!p5) { 682 /* first time */ 683 p5 = i2b(625); 684 if (p5 == NULL) { 685 Bfree(b); 686 return NULL; 687 } 688 p5s = p5; 689 p5->next = 0; 690 } 691 for(;;) { 692 if (k & 1) { 693 b1 = mult(b, p5); 694 Bfree(b); 695 b = b1; 696 if (b == NULL) 697 return NULL; 698 } 699 if (!(k >>= 1)) 700 break; 701 p51 = p5->next; 702 if (!p51) { 703 p51 = mult(p5,p5); 704 if (p51 == NULL) { 705 Bfree(b); 706 return NULL; 707 } 708 p51->next = 0; 709 p5->next = p51; 710 } 711 p5 = p51; 712 } 713 return b; 714} 715 716/* shift a Bigint b left by k bits. Return a pointer to the shifted result, 717 or NULL on failure. If the returned pointer is distinct from b then the 718 original b will have been Bfree'd. Ignores the sign of b. */ 719 720static Bigint * 721lshift(Bigint *b, int k) 722{ 723 int i, k1, n, n1; 724 Bigint *b1; 725 ULong *x, *x1, *xe, z; 726 727 n = k >> 5; 728 k1 = b->k; 729 n1 = n + b->wds + 1; 730 for(i = b->maxwds; n1 > i; i <<= 1) 731 k1++; 732 b1 = Balloc(k1); 733 if (b1 == NULL) { 734 Bfree(b); 735 return NULL; 736 } 737 x1 = b1->x; 738 for(i = 0; i < n; i++) 739 *x1++ = 0; 740 x = b->x; 741 xe = x + b->wds; 742 if (k &= 0x1f) { 743 k1 = 32 - k; 744 z = 0; 745 do { 746 *x1++ = *x << k | z; 747 z = *x++ >> k1; 748 } 749 while(x < xe); 750 if ((*x1 = z)) 751 ++n1; 752 } 753 else do 754 *x1++ = *x++; 755 while(x < xe); 756 b1->wds = n1 - 1; 757 Bfree(b); 758 return b1; 759} 760 761/* Do a three-way compare of a and b, returning -1 if a < b, 0 if a == b and 762 1 if a > b. Ignores signs of a and b. */ 763 764static int 765cmp(Bigint *a, Bigint *b) 766{ 767 ULong *xa, *xa0, *xb, *xb0; 768 int i, j; 769 770 i = a->wds; 771 j = b->wds; 772#ifdef DEBUG 773 if (i > 1 && !a->x[i-1]) 774 Bug("cmp called with a->x[a->wds-1] == 0"); 775 if (j > 1 && !b->x[j-1]) 776 Bug("cmp called with b->x[b->wds-1] == 0"); 777#endif 778 if (i -= j) 779 return i; 780 xa0 = a->x; 781 xa = xa0 + j; 782 xb0 = b->x; 783 xb = xb0 + j; 784 for(;;) { 785 if (*--xa != *--xb) 786 return *xa < *xb ? -1 : 1; 787 if (xa <= xa0) 788 break; 789 } 790 return 0; 791} 792 793/* Take the difference of Bigints a and b, returning a new Bigint. Returns 794 NULL on failure. The signs of a and b are ignored, but the sign of the 795 result is set appropriately. */ 796 797static Bigint * 798diff(Bigint *a, Bigint *b) 799{ 800 Bigint *c; 801 int i, wa, wb; 802 ULong *xa, *xae, *xb, *xbe, *xc; 803#ifdef ULLong 804 ULLong borrow, y; 805#else 806 ULong borrow, y; 807 ULong z; 808#endif 809 810 i = cmp(a,b); 811 if (!i) { 812 c = Balloc(0); 813 if (c == NULL) 814 return NULL; 815 c->wds = 1; 816 c->x[0] = 0; 817 return c; 818 } 819 if (i < 0) { 820 c = a; 821 a = b; 822 b = c; 823 i = 1; 824 } 825 else 826 i = 0; 827 c = Balloc(a->k); 828 if (c == NULL) 829 return NULL; 830 c->sign = i; 831 wa = a->wds; 832 xa = a->x; 833 xae = xa + wa; 834 wb = b->wds; 835 xb = b->x; 836 xbe = xb + wb; 837 xc = c->x; 838 borrow = 0; 839#ifdef ULLong 840 do { 841 y = (ULLong)*xa++ - *xb++ - borrow; 842 borrow = y >> 32 & (ULong)1; 843 *xc++ = (ULong)(y & FFFFFFFF); 844 } 845 while(xb < xbe); 846 while(xa < xae) { 847 y = *xa++ - borrow; 848 borrow = y >> 32 & (ULong)1; 849 *xc++ = (ULong)(y & FFFFFFFF); 850 } 851#else 852 do { 853 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow; 854 borrow = (y & 0x10000) >> 16; 855 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow; 856 borrow = (z & 0x10000) >> 16; 857 Storeinc(xc, z, y); 858 } 859 while(xb < xbe); 860 while(xa < xae) { 861 y = (*xa & 0xffff) - borrow; 862 borrow = (y & 0x10000) >> 16; 863 z = (*xa++ >> 16) - borrow; 864 borrow = (z & 0x10000) >> 16; 865 Storeinc(xc, z, y); 866 } 867#endif 868 while(!*--xc) 869 wa--; 870 c->wds = wa; 871 return c; 872} 873 874/* Given a positive normal double x, return the difference between x and the next 875 double up. Doesn't give correct results for subnormals. */ 876 877static double 878ulp(U *x) 879{ 880 Long L; 881 U u; 882 883 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; 884 word0(&u) = L; 885 word1(&u) = 0; 886 return dval(&u); 887} 888 889/* Convert a Bigint to a double plus an exponent */ 890 891static double 892b2d(Bigint *a, int *e) 893{ 894 ULong *xa, *xa0, w, y, z; 895 int k; 896 U d; 897 898 xa0 = a->x; 899 xa = xa0 + a->wds; 900 y = *--xa; 901#ifdef DEBUG 902 if (!y) Bug("zero y in b2d"); 903#endif 904 k = hi0bits(y); 905 *e = 32 - k; 906 if (k < Ebits) { 907 word0(&d) = Exp_1 | y >> (Ebits - k); 908 w = xa > xa0 ? *--xa : 0; 909 word1(&d) = y << ((32-Ebits) + k) | w >> (Ebits - k); 910 goto ret_d; 911 } 912 z = xa > xa0 ? *--xa : 0; 913 if (k -= Ebits) { 914 word0(&d) = Exp_1 | y << k | z >> (32 - k); 915 y = xa > xa0 ? *--xa : 0; 916 word1(&d) = z << k | y >> (32 - k); 917 } 918 else { 919 word0(&d) = Exp_1 | y; 920 word1(&d) = z; 921 } 922 ret_d: 923 return dval(&d); 924} 925 926/* Convert a double to a Bigint plus an exponent. Return NULL on failure. 927 928 Given a finite nonzero double d, return an odd Bigint b and exponent *e 929 such that fabs(d) = b * 2**e. On return, *bbits gives the number of 930 significant bits of b; that is, 2**(*bbits-1) <= b < 2**(*bbits). 931 932 If d is zero, then b == 0, *e == -1010, *bbits = 0. 933 */ 934 935 936static Bigint * 937d2b(U *d, int *e, int *bits) 938{ 939 Bigint *b; 940 int de, k; 941 ULong *x, y, z; 942 int i; 943 944 b = Balloc(1); 945 if (b == NULL) 946 return NULL; 947 x = b->x; 948 949 z = word0(d) & Frac_mask; 950 word0(d) &= 0x7fffffff; /* clear sign bit, which we ignore */ 951 if ((de = (int)(word0(d) >> Exp_shift))) 952 z |= Exp_msk1; 953 if ((y = word1(d))) { 954 if ((k = lo0bits(&y))) { 955 x[0] = y | z << (32 - k); 956 z >>= k; 957 } 958 else 959 x[0] = y; 960 i = 961 b->wds = (x[1] = z) ? 2 : 1; 962 } 963 else { 964 k = lo0bits(&z); 965 x[0] = z; 966 i = 967 b->wds = 1; 968 k += 32; 969 } 970 if (de) { 971 *e = de - Bias - (P-1) + k; 972 *bits = P - k; 973 } 974 else { 975 *e = de - Bias - (P-1) + 1 + k; 976 *bits = 32*i - hi0bits(x[i-1]); 977 } 978 return b; 979} 980 981/* Compute the ratio of two Bigints, as a double. The result may have an 982 error of up to 2.5 ulps. */ 983 984static double 985ratio(Bigint *a, Bigint *b) 986{ 987 U da, db; 988 int k, ka, kb; 989 990 dval(&da) = b2d(a, &ka); 991 dval(&db) = b2d(b, &kb); 992 k = ka - kb + 32*(a->wds - b->wds); 993 if (k > 0) 994 word0(&da) += k*Exp_msk1; 995 else { 996 k = -k; 997 word0(&db) += k*Exp_msk1; 998 } 999 return dval(&da) / dval(&db); 1000} 1001 1002static const double 1003tens[] = { 1004 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1005 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1006 1e20, 1e21, 1e22 1007}; 1008 1009static const double 1010bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; 1011static const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1012 9007199254740992.*9007199254740992.e-256 1013 /* = 2^106 * 1e-256 */ 1014}; 1015/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */ 1016/* flag unnecessarily. It leads to a song and dance at the end of strtod. */ 1017#define Scale_Bit 0x10 1018#define n_bigtens 5 1019 1020#define ULbits 32 1021#define kshift 5 1022#define kmask 31 1023 1024 1025static int 1026dshift(Bigint *b, int p2) 1027{ 1028 int rv = hi0bits(b->x[b->wds-1]) - 4; 1029 if (p2 > 0) 1030 rv -= p2; 1031 return rv & kmask; 1032} 1033 1034/* special case of Bigint division. The quotient is always in the range 0 <= 1035 quotient < 10, and on entry the divisor S is normalized so that its top 4 1036 bits (28--31) are zero and bit 27 is set. */ 1037 1038static int 1039quorem(Bigint *b, Bigint *S) 1040{ 1041 int n; 1042 ULong *bx, *bxe, q, *sx, *sxe; 1043#ifdef ULLong 1044 ULLong borrow, carry, y, ys; 1045#else 1046 ULong borrow, carry, y, ys; 1047 ULong si, z, zs; 1048#endif 1049 1050 n = S->wds; 1051#ifdef DEBUG 1052 /*debug*/ if (b->wds > n) 1053 /*debug*/ Bug("oversize b in quorem"); 1054#endif 1055 if (b->wds < n) 1056 return 0; 1057 sx = S->x; 1058 sxe = sx + --n; 1059 bx = b->x; 1060 bxe = bx + n; 1061 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ 1062#ifdef DEBUG 1063 /*debug*/ if (q > 9) 1064 /*debug*/ Bug("oversized quotient in quorem"); 1065#endif 1066 if (q) { 1067 borrow = 0; 1068 carry = 0; 1069 do { 1070#ifdef ULLong 1071 ys = *sx++ * (ULLong)q + carry; 1072 carry = ys >> 32; 1073 y = *bx - (ys & FFFFFFFF) - borrow; 1074 borrow = y >> 32 & (ULong)1; 1075 *bx++ = (ULong)(y & FFFFFFFF); 1076#else 1077 si = *sx++; 1078 ys = (si & 0xffff) * q + carry; 1079 zs = (si >> 16) * q + (ys >> 16); 1080 carry = zs >> 16; 1081 y = (*bx & 0xffff) - (ys & 0xffff) - borrow; 1082 borrow = (y & 0x10000) >> 16; 1083 z = (*bx >> 16) - (zs & 0xffff) - borrow; 1084 borrow = (z & 0x10000) >> 16; 1085 Storeinc(bx, z, y); 1086#endif 1087 } 1088 while(sx <= sxe); 1089 if (!*bxe) { 1090 bx = b->x; 1091 while(--bxe > bx && !*bxe) 1092 --n; 1093 b->wds = n; 1094 } 1095 } 1096 if (cmp(b, S) >= 0) { 1097 q++; 1098 borrow = 0; 1099 carry = 0; 1100 bx = b->x; 1101 sx = S->x; 1102 do { 1103#ifdef ULLong 1104 ys = *sx++ + carry; 1105 carry = ys >> 32; 1106 y = *bx - (ys & FFFFFFFF) - borrow; 1107 borrow = y >> 32 & (ULong)1; 1108 *bx++ = (ULong)(y & FFFFFFFF); 1109#else 1110 si = *sx++; 1111 ys = (si & 0xffff) + carry; 1112 zs = (si >> 16) + (ys >> 16); 1113 carry = zs >> 16; 1114 y = (*bx & 0xffff) - (ys & 0xffff) - borrow; 1115 borrow = (y & 0x10000) >> 16; 1116 z = (*bx >> 16) - (zs & 0xffff) - borrow; 1117 borrow = (z & 0x10000) >> 16; 1118 Storeinc(bx, z, y); 1119#endif 1120 } 1121 while(sx <= sxe); 1122 bx = b->x; 1123 bxe = bx + n; 1124 if (!*bxe) { 1125 while(--bxe > bx && !*bxe) 1126 --n; 1127 b->wds = n; 1128 } 1129 } 1130 return q; 1131} 1132 1133/* version of ulp(x) that takes bc.scale into account. 1134 1135 Assuming that x is finite and nonzero, and x / 2^bc.scale is exactly 1136 representable as a double, sulp(x) is equivalent to 2^bc.scale * ulp(x / 1137 2^bc.scale). */ 1138 1139static double 1140sulp(U *x, BCinfo *bc) 1141{ 1142 U u; 1143 1144 if (bc->scale && 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift) > 0) { 1145 /* rv/2^bc->scale is subnormal */ 1146 word0(&u) = (P+2)*Exp_msk1; 1147 word1(&u) = 0; 1148 return u.d; 1149 } 1150 else 1151 return ulp(x); 1152} 1153 1154/* return 0 on success, -1 on failure */ 1155 1156static int 1157bigcomp(U *rv, const char *s0, BCinfo *bc) 1158{ 1159 Bigint *b, *d; 1160 int b2, bbits, d2, dd, dig, dsign, i, j, nd, nd0, p2, p5, speccase; 1161 1162 dsign = bc->dsign; 1163 nd = bc->nd; 1164 nd0 = bc->nd0; 1165 p5 = nd + bc->e0; 1166 speccase = 0; 1167 if (rv->d == 0.) { /* special case: value near underflow-to-zero */ 1168 /* threshold was rounded to zero */ 1169 b = i2b(1); 1170 if (b == NULL) 1171 return -1; 1172 p2 = Emin - P + 1; 1173 bbits = 1; 1174 word0(rv) = (P+2) << Exp_shift; 1175 i = 0; 1176 { 1177 speccase = 1; 1178 --p2; 1179 dsign = 0; 1180 goto have_i; 1181 } 1182 } 1183 else 1184 { 1185 b = d2b(rv, &p2, &bbits); 1186 if (b == NULL) 1187 return -1; 1188 } 1189 p2 -= bc->scale; 1190 /* floor(log2(rv)) == bbits - 1 + p2 */ 1191 /* Check for denormal case. */ 1192 i = P - bbits; 1193 if (i > (j = P - Emin - 1 + p2)) { 1194 i = j; 1195 } 1196 { 1197 b = lshift(b, ++i); 1198 if (b == NULL) 1199 return -1; 1200 b->x[0] |= 1; 1201 } 1202 have_i: 1203 p2 -= p5 + i; 1204 d = i2b(1); 1205 if (d == NULL) { 1206 Bfree(b); 1207 return -1; 1208 } 1209 /* Arrange for convenient computation of quotients: 1210 * shift left if necessary so divisor has 4 leading 0 bits. 1211 */ 1212 if (p5 > 0) { 1213 d = pow5mult(d, p5); 1214 if (d == NULL) { 1215 Bfree(b); 1216 return -1; 1217 } 1218 } 1219 else if (p5 < 0) { 1220 b = pow5mult(b, -p5); 1221 if (b == NULL) { 1222 Bfree(d); 1223 return -1; 1224 } 1225 } 1226 if (p2 > 0) { 1227 b2 = p2; 1228 d2 = 0; 1229 } 1230 else { 1231 b2 = 0; 1232 d2 = -p2; 1233 } 1234 i = dshift(d, d2); 1235 if ((b2 += i) > 0) { 1236 b = lshift(b, b2); 1237 if (b == NULL) { 1238 Bfree(d); 1239 return -1; 1240 } 1241 } 1242 if ((d2 += i) > 0) { 1243 d = lshift(d, d2); 1244 if (d == NULL) { 1245 Bfree(b); 1246 return -1; 1247 } 1248 } 1249 1250 /* Now 10*b/d = exactly half-way between the two floating-point values 1251 on either side of the input string. If b >= d, round down. */ 1252 if (cmp(b, d) >= 0) { 1253 dd = -1; 1254 goto ret; 1255 } 1256 1257 /* Compute first digit of 10*b/d. */ 1258 b = multadd(b, 10, 0); 1259 if (b == NULL) { 1260 Bfree(d); 1261 return -1; 1262 } 1263 dig = quorem(b, d); 1264 assert(dig < 10); 1265 1266 /* Compare b/d with s0 */ 1267 1268 assert(nd > 0); 1269 dd = 9999; /* silence gcc compiler warning */ 1270 for(i = 0; i < nd0; ) { 1271 if ((dd = s0[i++] - '0' - dig)) 1272 goto ret; 1273 if (!b->x[0] && b->wds == 1) { 1274 if (i < nd) 1275 dd = 1; 1276 goto ret; 1277 } 1278 b = multadd(b, 10, 0); 1279 if (b == NULL) { 1280 Bfree(d); 1281 return -1; 1282 } 1283 dig = quorem(b,d); 1284 } 1285 for(j = bc->dp1; i++ < nd;) { 1286 if ((dd = s0[j++] - '0' - dig)) 1287 goto ret; 1288 if (!b->x[0] && b->wds == 1) { 1289 if (i < nd) 1290 dd = 1; 1291 goto ret; 1292 } 1293 b = multadd(b, 10, 0); 1294 if (b == NULL) { 1295 Bfree(d); 1296 return -1; 1297 } 1298 dig = quorem(b,d); 1299 } 1300 if (b->x[0] || b->wds > 1) 1301 dd = -1; 1302 ret: 1303 Bfree(b); 1304 Bfree(d); 1305 if (speccase) { 1306 if (dd <= 0) 1307 rv->d = 0.; 1308 } 1309 else if (dd < 0) { 1310 if (!dsign) /* does not happen for round-near */ 1311 retlow1: 1312 dval(rv) -= sulp(rv, bc); 1313 } 1314 else if (dd > 0) { 1315 if (dsign) { 1316 rethi1: 1317 dval(rv) += sulp(rv, bc); 1318 } 1319 } 1320 else { 1321 /* Exact half-way case: apply round-even rule. */ 1322 if (word1(rv) & 1) { 1323 if (dsign) 1324 goto rethi1; 1325 goto retlow1; 1326 } 1327 } 1328 1329 return 0; 1330} 1331 1332double 1333_Py_dg_strtod(const char *s00, char **se) 1334{ 1335 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, e, e1, error; 1336 int esign, i, j, k, nd, nd0, nf, nz, nz0, sign; 1337 const char *s, *s0, *s1; 1338 double aadj, aadj1; 1339 U aadj2, adj, rv, rv0; 1340 ULong y, z, L; 1341 BCinfo bc; 1342 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; 1343 1344 sign = nz0 = nz = bc.dplen = 0; 1345 dval(&rv) = 0.; 1346 for(s = s00;;s++) switch(*s) { 1347 case '-': 1348 sign = 1; 1349 /* no break */ 1350 case '+': 1351 if (*++s) 1352 goto break2; 1353 /* no break */ 1354 case 0: 1355 goto ret0; 1356 /* modify original dtoa.c so that it doesn't accept leading whitespace 1357 case '\t': 1358 case '\n': 1359 case '\v': 1360 case '\f': 1361 case '\r': 1362 case ' ': 1363 continue; 1364 */ 1365 default: 1366 goto break2; 1367 } 1368 break2: 1369 if (*s == '0') { 1370 nz0 = 1; 1371 while(*++s == '0') ; 1372 if (!*s) 1373 goto ret; 1374 } 1375 s0 = s; 1376 y = z = 0; 1377 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 1378 if (nd < 9) 1379 y = 10*y + c - '0'; 1380 else if (nd < 16) 1381 z = 10*z + c - '0'; 1382 nd0 = nd; 1383 bc.dp0 = bc.dp1 = s - s0; 1384 if (c == '.') { 1385 c = *++s; 1386 bc.dp1 = s - s0; 1387 bc.dplen = bc.dp1 - bc.dp0; 1388 if (!nd) { 1389 for(; c == '0'; c = *++s) 1390 nz++; 1391 if (c > '0' && c <= '9') { 1392 s0 = s; 1393 nf += nz; 1394 nz = 0; 1395 goto have_dig; 1396 } 1397 goto dig_done; 1398 } 1399 for(; c >= '0' && c <= '9'; c = *++s) { 1400 have_dig: 1401 nz++; 1402 if (c -= '0') { 1403 nf += nz; 1404 for(i = 1; i < nz; i++) 1405 if (nd++ < 9) 1406 y *= 10; 1407 else if (nd <= DBL_DIG + 1) 1408 z *= 10; 1409 if (nd++ < 9) 1410 y = 10*y + c; 1411 else if (nd <= DBL_DIG + 1) 1412 z = 10*z + c; 1413 nz = 0; 1414 } 1415 } 1416 } 1417 dig_done: 1418 e = 0; 1419 if (c == 'e' || c == 'E') { 1420 if (!nd && !nz && !nz0) { 1421 goto ret0; 1422 } 1423 s00 = s; 1424 esign = 0; 1425 switch(c = *++s) { 1426 case '-': 1427 esign = 1; 1428 case '+': 1429 c = *++s; 1430 } 1431 if (c >= '0' && c <= '9') { 1432 while(c == '0') 1433 c = *++s; 1434 if (c > '0' && c <= '9') { 1435 L = c - '0'; 1436 s1 = s; 1437 while((c = *++s) >= '0' && c <= '9') 1438 L = 10*L + c - '0'; 1439 if (s - s1 > 8 || L > MAX_ABS_EXP) 1440 /* Avoid confusion from exponents 1441 * so large that e might overflow. 1442 */ 1443 e = (int)MAX_ABS_EXP; /* safe for 16 bit ints */ 1444 else 1445 e = (int)L; 1446 if (esign) 1447 e = -e; 1448 } 1449 else 1450 e = 0; 1451 } 1452 else 1453 s = s00; 1454 } 1455 if (!nd) { 1456 if (!nz && !nz0) { 1457 ret0: 1458 s = s00; 1459 sign = 0; 1460 } 1461 goto ret; 1462 } 1463 bc.e0 = e1 = e -= nf; 1464 1465 /* Now we have nd0 digits, starting at s0, followed by a 1466 * decimal point, followed by nd-nd0 digits. The number we're 1467 * after is the integer represented by those digits times 1468 * 10**e */ 1469 1470 if (!nd0) 1471 nd0 = nd; 1472 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 1473 dval(&rv) = y; 1474 if (k > 9) { 1475 dval(&rv) = tens[k - 9] * dval(&rv) + z; 1476 } 1477 bd0 = 0; 1478 if (nd <= DBL_DIG 1479 && Flt_Rounds == 1 1480 ) { 1481 if (!e) 1482 goto ret; 1483 if (e > 0) { 1484 if (e <= Ten_pmax) { 1485 dval(&rv) *= tens[e]; 1486 goto ret; 1487 } 1488 i = DBL_DIG - nd; 1489 if (e <= Ten_pmax + i) { 1490 /* A fancier test would sometimes let us do 1491 * this for larger i values. 1492 */ 1493 e -= i; 1494 dval(&rv) *= tens[i]; 1495 dval(&rv) *= tens[e]; 1496 goto ret; 1497 } 1498 } 1499 else if (e >= -Ten_pmax) { 1500 dval(&rv) /= tens[-e]; 1501 goto ret; 1502 } 1503 } 1504 e1 += nd - k; 1505 1506 bc.scale = 0; 1507 1508 /* Get starting approximation = rv * 10**e1 */ 1509 1510 if (e1 > 0) { 1511 if ((i = e1 & 15)) 1512 dval(&rv) *= tens[i]; 1513 if (e1 &= ~15) { 1514 if (e1 > DBL_MAX_10_EXP) { 1515 ovfl: 1516 errno = ERANGE; 1517 /* Can't trust HUGE_VAL */ 1518 word0(&rv) = Exp_mask; 1519 word1(&rv) = 0; 1520 goto ret; 1521 } 1522 e1 >>= 4; 1523 for(j = 0; e1 > 1; j++, e1 >>= 1) 1524 if (e1 & 1) 1525 dval(&rv) *= bigtens[j]; 1526 /* The last multiplication could overflow. */ 1527 word0(&rv) -= P*Exp_msk1; 1528 dval(&rv) *= bigtens[j]; 1529 if ((z = word0(&rv) & Exp_mask) 1530 > Exp_msk1*(DBL_MAX_EXP+Bias-P)) 1531 goto ovfl; 1532 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { 1533 /* set to largest number */ 1534 /* (Can't trust DBL_MAX) */ 1535 word0(&rv) = Big0; 1536 word1(&rv) = Big1; 1537 } 1538 else 1539 word0(&rv) += P*Exp_msk1; 1540 } 1541 } 1542 else if (e1 < 0) { 1543 e1 = -e1; 1544 if ((i = e1 & 15)) 1545 dval(&rv) /= tens[i]; 1546 if (e1 >>= 4) { 1547 if (e1 >= 1 << n_bigtens) 1548 goto undfl; 1549 if (e1 & Scale_Bit) 1550 bc.scale = 2*P; 1551 for(j = 0; e1 > 0; j++, e1 >>= 1) 1552 if (e1 & 1) 1553 dval(&rv) *= tinytens[j]; 1554 if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask) 1555 >> Exp_shift)) > 0) { 1556 /* scaled rv is denormal; clear j low bits */ 1557 if (j >= 32) { 1558 word1(&rv) = 0; 1559 if (j >= 53) 1560 word0(&rv) = (P+2)*Exp_msk1; 1561 else 1562 word0(&rv) &= 0xffffffff << (j-32); 1563 } 1564 else 1565 word1(&rv) &= 0xffffffff << j; 1566 } 1567 if (!dval(&rv)) { 1568 undfl: 1569 dval(&rv) = 0.; 1570 errno = ERANGE; 1571 goto ret; 1572 } 1573 } 1574 } 1575 1576 /* Now the hard part -- adjusting rv to the correct value.*/ 1577 1578 /* Put digits into bd: true value = bd * 10^e */ 1579 1580 bc.nd = nd; 1581 bc.nd0 = nd0; /* Only needed if nd > STRTOD_DIGLIM, but done here */ 1582 /* to silence an erroneous warning about bc.nd0 */ 1583 /* possibly not being initialized. */ 1584 if (nd > STRTOD_DIGLIM) { 1585 /* ASSERT(STRTOD_DIGLIM >= 18); 18 == one more than the */ 1586 /* minimum number of decimal digits to distinguish double values */ 1587 /* in IEEE arithmetic. */ 1588 i = j = 18; 1589 if (i > nd0) 1590 j += bc.dplen; 1591 for(;;) { 1592 if (--j <= bc.dp1 && j >= bc.dp0) 1593 j = bc.dp0 - 1; 1594 if (s0[j] != '0') 1595 break; 1596 --i; 1597 } 1598 e += nd - i; 1599 nd = i; 1600 if (nd0 > nd) 1601 nd0 = nd; 1602 if (nd < 9) { /* must recompute y */ 1603 y = 0; 1604 for(i = 0; i < nd0; ++i) 1605 y = 10*y + s0[i] - '0'; 1606 for(j = bc.dp1; i < nd; ++i) 1607 y = 10*y + s0[j++] - '0'; 1608 } 1609 } 1610 bd0 = s2b(s0, nd0, nd, y, bc.dplen); 1611 if (bd0 == NULL) 1612 goto failed_malloc; 1613 1614 for(;;) { 1615 bd = Balloc(bd0->k); 1616 if (bd == NULL) { 1617 Bfree(bd0); 1618 goto failed_malloc; 1619 } 1620 Bcopy(bd, bd0); 1621 bb = d2b(&rv, &bbe, &bbbits); /* rv = bb * 2^bbe */ 1622 if (bb == NULL) { 1623 Bfree(bd); 1624 Bfree(bd0); 1625 goto failed_malloc; 1626 } 1627 bs = i2b(1); 1628 if (bs == NULL) { 1629 Bfree(bb); 1630 Bfree(bd); 1631 Bfree(bd0); 1632 goto failed_malloc; 1633 } 1634 1635 if (e >= 0) { 1636 bb2 = bb5 = 0; 1637 bd2 = bd5 = e; 1638 } 1639 else { 1640 bb2 = bb5 = -e; 1641 bd2 = bd5 = 0; 1642 } 1643 if (bbe >= 0) 1644 bb2 += bbe; 1645 else 1646 bd2 -= bbe; 1647 bs2 = bb2; 1648 j = bbe - bc.scale; 1649 i = j + bbbits - 1; /* logb(rv) */ 1650 if (i < Emin) /* denormal */ 1651 j += P - Emin; 1652 else 1653 j = P + 1 - bbbits; 1654 bb2 += j; 1655 bd2 += j; 1656 bd2 += bc.scale; 1657 i = bb2 < bd2 ? bb2 : bd2; 1658 if (i > bs2) 1659 i = bs2; 1660 if (i > 0) { 1661 bb2 -= i; 1662 bd2 -= i; 1663 bs2 -= i; 1664 } 1665 if (bb5 > 0) { 1666 bs = pow5mult(bs, bb5); 1667 if (bs == NULL) { 1668 Bfree(bb); 1669 Bfree(bd); 1670 Bfree(bd0); 1671 goto failed_malloc; 1672 } 1673 bb1 = mult(bs, bb); 1674 Bfree(bb); 1675 bb = bb1; 1676 if (bb == NULL) { 1677 Bfree(bs); 1678 Bfree(bd); 1679 Bfree(bd0); 1680 goto failed_malloc; 1681 } 1682 } 1683 if (bb2 > 0) { 1684 bb = lshift(bb, bb2); 1685 if (bb == NULL) { 1686 Bfree(bs); 1687 Bfree(bd); 1688 Bfree(bd0); 1689 goto failed_malloc; 1690 } 1691 } 1692 if (bd5 > 0) { 1693 bd = pow5mult(bd, bd5); 1694 if (bd == NULL) { 1695 Bfree(bb); 1696 Bfree(bs); 1697 Bfree(bd0); 1698 goto failed_malloc; 1699 } 1700 } 1701 if (bd2 > 0) { 1702 bd = lshift(bd, bd2); 1703 if (bd == NULL) { 1704 Bfree(bb); 1705 Bfree(bs); 1706 Bfree(bd0); 1707 goto failed_malloc; 1708 } 1709 } 1710 if (bs2 > 0) { 1711 bs = lshift(bs, bs2); 1712 if (bs == NULL) { 1713 Bfree(bb); 1714 Bfree(bd); 1715 Bfree(bd0); 1716 goto failed_malloc; 1717 } 1718 } 1719 delta = diff(bb, bd); 1720 if (delta == NULL) { 1721 Bfree(bb); 1722 Bfree(bs); 1723 Bfree(bd); 1724 Bfree(bd0); 1725 goto failed_malloc; 1726 } 1727 bc.dsign = delta->sign; 1728 delta->sign = 0; 1729 i = cmp(delta, bs); 1730 if (bc.nd > nd && i <= 0) { 1731 if (bc.dsign) 1732 break; /* Must use bigcomp(). */ 1733 { 1734 bc.nd = nd; 1735 i = -1; /* Discarded digits make delta smaller. */ 1736 } 1737 } 1738 1739 if (i < 0) { 1740 /* Error is less than half an ulp -- check for 1741 * special case of mantissa a power of two. 1742 */ 1743 if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask 1744 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1 1745 ) { 1746 break; 1747 } 1748 if (!delta->x[0] && delta->wds <= 1) { 1749 /* exact result */ 1750 break; 1751 } 1752 delta = lshift(delta,Log2P); 1753 if (delta == NULL) { 1754 Bfree(bb); 1755 Bfree(bs); 1756 Bfree(bd); 1757 Bfree(bd0); 1758 goto failed_malloc; 1759 } 1760 if (cmp(delta, bs) > 0) 1761 goto drop_down; 1762 break; 1763 } 1764 if (i == 0) { 1765 /* exactly half-way between */ 1766 if (bc.dsign) { 1767 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1 1768 && word1(&rv) == ( 1769 (bc.scale && 1770 (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) ? 1771 (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : 1772 0xffffffff)) { 1773 /*boundary case -- increment exponent*/ 1774 word0(&rv) = (word0(&rv) & Exp_mask) 1775 + Exp_msk1 1776 ; 1777 word1(&rv) = 0; 1778 bc.dsign = 0; 1779 break; 1780 } 1781 } 1782 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) { 1783 drop_down: 1784 /* boundary case -- decrement exponent */ 1785 if (bc.scale) { 1786 L = word0(&rv) & Exp_mask; 1787 if (L <= (2*P+1)*Exp_msk1) { 1788 if (L > (P+2)*Exp_msk1) 1789 /* round even ==> */ 1790 /* accept rv */ 1791 break; 1792 /* rv = smallest denormal */ 1793 if (bc.nd >nd) 1794 break; 1795 goto undfl; 1796 } 1797 } 1798 L = (word0(&rv) & Exp_mask) - Exp_msk1; 1799 word0(&rv) = L | Bndry_mask1; 1800 word1(&rv) = 0xffffffff; 1801 break; 1802 } 1803 if (!(word1(&rv) & LSB)) 1804 break; 1805 if (bc.dsign) 1806 dval(&rv) += ulp(&rv); 1807 else { 1808 dval(&rv) -= ulp(&rv); 1809 if (!dval(&rv)) { 1810 if (bc.nd >nd) 1811 break; 1812 goto undfl; 1813 } 1814 } 1815 bc.dsign = 1 - bc.dsign; 1816 break; 1817 } 1818 if ((aadj = ratio(delta, bs)) <= 2.) { 1819 if (bc.dsign) 1820 aadj = aadj1 = 1.; 1821 else if (word1(&rv) || word0(&rv) & Bndry_mask) { 1822 if (word1(&rv) == Tiny1 && !word0(&rv)) { 1823 if (bc.nd >nd) 1824 break; 1825 goto undfl; 1826 } 1827 aadj = 1.; 1828 aadj1 = -1.; 1829 } 1830 else { 1831 /* special case -- power of FLT_RADIX to be */ 1832 /* rounded down... */ 1833 1834 if (aadj < 2./FLT_RADIX) 1835 aadj = 1./FLT_RADIX; 1836 else 1837 aadj *= 0.5; 1838 aadj1 = -aadj; 1839 } 1840 } 1841 else { 1842 aadj *= 0.5; 1843 aadj1 = bc.dsign ? aadj : -aadj; 1844 if (Flt_Rounds == 0) 1845 aadj1 += 0.5; 1846 } 1847 y = word0(&rv) & Exp_mask; 1848 1849 /* Check for overflow */ 1850 1851 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 1852 dval(&rv0) = dval(&rv); 1853 word0(&rv) -= P*Exp_msk1; 1854 adj.d = aadj1 * ulp(&rv); 1855 dval(&rv) += adj.d; 1856 if ((word0(&rv) & Exp_mask) >= 1857 Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 1858 if (word0(&rv0) == Big0 && word1(&rv0) == Big1) 1859 goto ovfl; 1860 word0(&rv) = Big0; 1861 word1(&rv) = Big1; 1862 goto cont; 1863 } 1864 else 1865 word0(&rv) += P*Exp_msk1; 1866 } 1867 else { 1868 if (bc.scale && y <= 2*P*Exp_msk1) { 1869 if (aadj <= 0x7fffffff) { 1870 if ((z = (ULong)aadj) <= 0) 1871 z = 1; 1872 aadj = z; 1873 aadj1 = bc.dsign ? aadj : -aadj; 1874 } 1875 dval(&aadj2) = aadj1; 1876 word0(&aadj2) += (2*P+1)*Exp_msk1 - y; 1877 aadj1 = dval(&aadj2); 1878 } 1879 adj.d = aadj1 * ulp(&rv); 1880 dval(&rv) += adj.d; 1881 } 1882 z = word0(&rv) & Exp_mask; 1883 if (bc.nd == nd) { 1884 if (!bc.scale) 1885 if (y == z) { 1886 /* Can we stop now? */ 1887 L = (Long)aadj; 1888 aadj -= L; 1889 /* The tolerances below are conservative. */ 1890 if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask) { 1891 if (aadj < .4999999 || aadj > .5000001) 1892 break; 1893 } 1894 else if (aadj < .4999999/FLT_RADIX) 1895 break; 1896 } 1897 } 1898 cont: 1899 Bfree(bb); 1900 Bfree(bd); 1901 Bfree(bs); 1902 Bfree(delta); 1903 } 1904 Bfree(bb); 1905 Bfree(bd); 1906 Bfree(bs); 1907 Bfree(bd0); 1908 Bfree(delta); 1909 if (bc.nd > nd) { 1910 error = bigcomp(&rv, s0, &bc); 1911 if (error) 1912 goto failed_malloc; 1913 } 1914 1915 if (bc.scale) { 1916 word0(&rv0) = Exp_1 - 2*P*Exp_msk1; 1917 word1(&rv0) = 0; 1918 dval(&rv) *= dval(&rv0); 1919 /* try to avoid the bug of testing an 8087 register value */ 1920 if (!(word0(&rv) & Exp_mask)) 1921 errno = ERANGE; 1922 } 1923 ret: 1924 if (se) 1925 *se = (char *)s; 1926 return sign ? -dval(&rv) : dval(&rv); 1927 1928 failed_malloc: 1929 if (se) 1930 *se = (char *)s00; 1931 errno = ENOMEM; 1932 return -1.0; 1933} 1934 1935static char * 1936rv_alloc(int i) 1937{ 1938 int j, k, *r; 1939 1940 j = sizeof(ULong); 1941 for(k = 0; 1942 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (unsigned)i; 1943 j <<= 1) 1944 k++; 1945 r = (int*)Balloc(k); 1946 if (r == NULL) 1947 return NULL; 1948 *r = k; 1949 return (char *)(r+1); 1950} 1951 1952static char * 1953nrv_alloc(char *s, char **rve, int n) 1954{ 1955 char *rv, *t; 1956 1957 rv = rv_alloc(n); 1958 if (rv == NULL) 1959 return NULL; 1960 t = rv; 1961 while((*t = *s++)) t++; 1962 if (rve) 1963 *rve = t; 1964 return rv; 1965} 1966 1967/* freedtoa(s) must be used to free values s returned by dtoa 1968 * when MULTIPLE_THREADS is #defined. It should be used in all cases, 1969 * but for consistency with earlier versions of dtoa, it is optional 1970 * when MULTIPLE_THREADS is not defined. 1971 */ 1972 1973void 1974_Py_dg_freedtoa(char *s) 1975{ 1976 Bigint *b = (Bigint *)((int *)s - 1); 1977 b->maxwds = 1 << (b->k = *(int*)b); 1978 Bfree(b); 1979} 1980 1981/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. 1982 * 1983 * Inspired by "How to Print Floating-Point Numbers Accurately" by 1984 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126]. 1985 * 1986 * Modifications: 1987 * 1. Rather than iterating, we use a simple numeric overestimate 1988 * to determine k = floor(log10(d)). We scale relevant 1989 * quantities using O(log2(k)) rather than O(k) multiplications. 1990 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't 1991 * try to generate digits strictly left to right. Instead, we 1992 * compute with fewer bits and propagate the carry if necessary 1993 * when rounding the final digit up. This is often faster. 1994 * 3. Under the assumption that input will be rounded nearest, 1995 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. 1996 * That is, we allow equality in stopping tests when the 1997 * round-nearest rule will give the same floating-point value 1998 * as would satisfaction of the stopping test with strict 1999 * inequality. 2000 * 4. We remove common factors of powers of 2 from relevant 2001 * quantities. 2002 * 5. When converting floating-point integers less than 1e16, 2003 * we use floating-point arithmetic rather than resorting 2004 * to multiple-precision integers. 2005 * 6. When asked to produce fewer than 15 digits, we first try 2006 * to get by with floating-point arithmetic; we resort to 2007 * multiple-precision integer arithmetic only if we cannot 2008 * guarantee that the floating-point calculation has given 2009 * the correctly rounded result. For k requested digits and 2010 * "uniformly" distributed input, the probability is 2011 * something like 10^(k-15) that we must resort to the Long 2012 * calculation. 2013 */ 2014 2015/* Additional notes (METD): (1) returns NULL on failure. (2) to avoid memory 2016 leakage, a successful call to _Py_dg_dtoa should always be matched by a 2017 call to _Py_dg_freedtoa. */ 2018 2019char * 2020_Py_dg_dtoa(double dd, int mode, int ndigits, 2021 int *decpt, int *sign, char **rve) 2022{ 2023 /* Arguments ndigits, decpt, sign are similar to those 2024 of ecvt and fcvt; trailing zeros are suppressed from 2025 the returned string. If not null, *rve is set to point 2026 to the end of the return value. If d is +-Infinity or NaN, 2027 then *decpt is set to 9999. 2028 2029 mode: 2030 0 ==> shortest string that yields d when read in 2031 and rounded to nearest. 2032 1 ==> like 0, but with Steele & White stopping rule; 2033 e.g. with IEEE P754 arithmetic , mode 0 gives 2034 1e23 whereas mode 1 gives 9.999999999999999e22. 2035 2 ==> max(1,ndigits) significant digits. This gives a 2036 return value similar to that of ecvt, except 2037 that trailing zeros are suppressed. 2038 3 ==> through ndigits past the decimal point. This 2039 gives a return value similar to that from fcvt, 2040 except that trailing zeros are suppressed, and 2041 ndigits can be negative. 2042 4,5 ==> similar to 2 and 3, respectively, but (in 2043 round-nearest mode) with the tests of mode 0 to 2044 possibly return a shorter string that rounds to d. 2045 With IEEE arithmetic and compilation with 2046 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same 2047 as modes 2 and 3 when FLT_ROUNDS != 1. 2048 6-9 ==> Debugging modes similar to mode - 4: don't try 2049 fast floating-point estimate (if applicable). 2050 2051 Values of mode other than 0-9 are treated as mode 0. 2052 2053 Sufficient space is allocated to the return value 2054 to hold the suppressed trailing zeros. 2055 */ 2056 2057 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1, 2058 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5, 2059 spec_case, try_quick; 2060 Long L; 2061 int denorm; 2062 ULong x; 2063 Bigint *b, *b1, *delta, *mlo, *mhi, *S; 2064 U d2, eps, u; 2065 double ds; 2066 char *s, *s0; 2067 2068 /* set pointers to NULL, to silence gcc compiler warnings and make 2069 cleanup easier on error */ 2070 mlo = mhi = b = S = 0; 2071 s0 = 0; 2072 2073 u.d = dd; 2074 if (word0(&u) & Sign_bit) { 2075 /* set sign for everything, including 0's and NaNs */ 2076 *sign = 1; 2077 word0(&u) &= ~Sign_bit; /* clear sign bit */ 2078 } 2079 else 2080 *sign = 0; 2081 2082 /* quick return for Infinities, NaNs and zeros */ 2083 if ((word0(&u) & Exp_mask) == Exp_mask) 2084 { 2085 /* Infinity or NaN */ 2086 *decpt = 9999; 2087 if (!word1(&u) && !(word0(&u) & 0xfffff)) 2088 return nrv_alloc("Infinity", rve, 8); 2089 return nrv_alloc("NaN", rve, 3); 2090 } 2091 if (!dval(&u)) { 2092 *decpt = 1; 2093 return nrv_alloc("0", rve, 1); 2094 } 2095 2096 /* compute k = floor(log10(d)). The computation may leave k 2097 one too large, but should never leave k too small. */ 2098 b = d2b(&u, &be, &bbits); 2099 if (b == NULL) 2100 goto failed_malloc; 2101 if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) { 2102 dval(&d2) = dval(&u); 2103 word0(&d2) &= Frac_mask1; 2104 word0(&d2) |= Exp_11; 2105 2106 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 2107 * log10(x) = log(x) / log(10) 2108 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) 2109 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) 2110 * 2111 * This suggests computing an approximation k to log10(d) by 2112 * 2113 * k = (i - Bias)*0.301029995663981 2114 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); 2115 * 2116 * We want k to be too large rather than too small. 2117 * The error in the first-order Taylor series approximation 2118 * is in our favor, so we just round up the constant enough 2119 * to compensate for any error in the multiplication of 2120 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, 2121 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, 2122 * adding 1e-13 to the constant term more than suffices. 2123 * Hence we adjust the constant term to 0.1760912590558. 2124 * (We could get a more accurate k by invoking log10, 2125 * but this is probably not worthwhile.) 2126 */ 2127 2128 i -= Bias; 2129 denorm = 0; 2130 } 2131 else { 2132 /* d is denormalized */ 2133 2134 i = bbits + be + (Bias + (P-1) - 1); 2135 x = i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32) 2136 : word1(&u) << (32 - i); 2137 dval(&d2) = x; 2138 word0(&d2) -= 31*Exp_msk1; /* adjust exponent */ 2139 i -= (Bias + (P-1) - 1) + 1; 2140 denorm = 1; 2141 } 2142 ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + 2143 i*0.301029995663981; 2144 k = (int)ds; 2145 if (ds < 0. && ds != k) 2146 k--; /* want k = floor(ds) */ 2147 k_check = 1; 2148 if (k >= 0 && k <= Ten_pmax) { 2149 if (dval(&u) < tens[k]) 2150 k--; 2151 k_check = 0; 2152 } 2153 j = bbits - i - 1; 2154 if (j >= 0) { 2155 b2 = 0; 2156 s2 = j; 2157 } 2158 else { 2159 b2 = -j; 2160 s2 = 0; 2161 } 2162 if (k >= 0) { 2163 b5 = 0; 2164 s5 = k; 2165 s2 += k; 2166 } 2167 else { 2168 b2 -= k; 2169 b5 = -k; 2170 s5 = 0; 2171 } 2172 if (mode < 0 || mode > 9) 2173 mode = 0; 2174 2175 try_quick = 1; 2176 2177 if (mode > 5) { 2178 mode -= 4; 2179 try_quick = 0; 2180 } 2181 leftright = 1; 2182 ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */ 2183 /* silence erroneous "gcc -Wall" warning. */ 2184 switch(mode) { 2185 case 0: 2186 case 1: 2187 i = 18; 2188 ndigits = 0; 2189 break; 2190 case 2: 2191 leftright = 0; 2192 /* no break */ 2193 case 4: 2194 if (ndigits <= 0) 2195 ndigits = 1; 2196 ilim = ilim1 = i = ndigits; 2197 break; 2198 case 3: 2199 leftright = 0; 2200 /* no break */ 2201 case 5: 2202 i = ndigits + k + 1; 2203 ilim = i; 2204 ilim1 = i - 1; 2205 if (i <= 0) 2206 i = 1; 2207 } 2208 s0 = rv_alloc(i); 2209 if (s0 == NULL) 2210 goto failed_malloc; 2211 s = s0; 2212 2213 2214 if (ilim >= 0 && ilim <= Quick_max && try_quick) { 2215 2216 /* Try to get by with floating-point arithmetic. */ 2217 2218 i = 0; 2219 dval(&d2) = dval(&u); 2220 k0 = k; 2221 ilim0 = ilim; 2222 ieps = 2; /* conservative */ 2223 if (k > 0) { 2224 ds = tens[k&0xf]; 2225 j = k >> 4; 2226 if (j & Bletch) { 2227 /* prevent overflows */ 2228 j &= Bletch - 1; 2229 dval(&u) /= bigtens[n_bigtens-1]; 2230 ieps++; 2231 } 2232 for(; j; j >>= 1, i++) 2233 if (j & 1) { 2234 ieps++; 2235 ds *= bigtens[i]; 2236 } 2237 dval(&u) /= ds; 2238 } 2239 else if ((j1 = -k)) { 2240 dval(&u) *= tens[j1 & 0xf]; 2241 for(j = j1 >> 4; j; j >>= 1, i++) 2242 if (j & 1) { 2243 ieps++; 2244 dval(&u) *= bigtens[i]; 2245 } 2246 } 2247 if (k_check && dval(&u) < 1. && ilim > 0) { 2248 if (ilim1 <= 0) 2249 goto fast_failed; 2250 ilim = ilim1; 2251 k--; 2252 dval(&u) *= 10.; 2253 ieps++; 2254 } 2255 dval(&eps) = ieps*dval(&u) + 7.; 2256 word0(&eps) -= (P-1)*Exp_msk1; 2257 if (ilim == 0) { 2258 S = mhi = 0; 2259 dval(&u) -= 5.; 2260 if (dval(&u) > dval(&eps)) 2261 goto one_digit; 2262 if (dval(&u) < -dval(&eps)) 2263 goto no_digits; 2264 goto fast_failed; 2265 } 2266 if (leftright) { 2267 /* Use Steele & White method of only 2268 * generating digits needed. 2269 */ 2270 dval(&eps) = 0.5/tens[ilim-1] - dval(&eps); 2271 for(i = 0;;) { 2272 L = (Long)dval(&u); 2273 dval(&u) -= L; 2274 *s++ = '0' + (int)L; 2275 if (dval(&u) < dval(&eps)) 2276 goto ret1; 2277 if (1. - dval(&u) < dval(&eps)) 2278 goto bump_up; 2279 if (++i >= ilim) 2280 break; 2281 dval(&eps) *= 10.; 2282 dval(&u) *= 10.; 2283 } 2284 } 2285 else { 2286 /* Generate ilim digits, then fix them up. */ 2287 dval(&eps) *= tens[ilim-1]; 2288 for(i = 1;; i++, dval(&u) *= 10.) { 2289 L = (Long)(dval(&u)); 2290 if (!(dval(&u) -= L)) 2291 ilim = i; 2292 *s++ = '0' + (int)L; 2293 if (i == ilim) { 2294 if (dval(&u) > 0.5 + dval(&eps)) 2295 goto bump_up; 2296 else if (dval(&u) < 0.5 - dval(&eps)) { 2297 while(*--s == '0'); 2298 s++; 2299 goto ret1; 2300 } 2301 break; 2302 } 2303 } 2304 } 2305 fast_failed: 2306 s = s0; 2307 dval(&u) = dval(&d2); 2308 k = k0; 2309 ilim = ilim0; 2310 } 2311 2312 /* Do we have a "small" integer? */ 2313 2314 if (be >= 0 && k <= Int_max) { 2315 /* Yes. */ 2316 ds = tens[k]; 2317 if (ndigits < 0 && ilim <= 0) { 2318 S = mhi = 0; 2319 if (ilim < 0 || dval(&u) <= 5*ds) 2320 goto no_digits; 2321 goto one_digit; 2322 } 2323 for(i = 1;; i++, dval(&u) *= 10.) { 2324 L = (Long)(dval(&u) / ds); 2325 dval(&u) -= L*ds; 2326 *s++ = '0' + (int)L; 2327 if (!dval(&u)) { 2328 break; 2329 } 2330 if (i == ilim) { 2331 dval(&u) += dval(&u); 2332 if (dval(&u) > ds || (dval(&u) == ds && L & 1)) { 2333 bump_up: 2334 while(*--s == '9') 2335 if (s == s0) { 2336 k++; 2337 *s = '0'; 2338 break; 2339 } 2340 ++*s++; 2341 } 2342 break; 2343 } 2344 } 2345 goto ret1; 2346 } 2347 2348 m2 = b2; 2349 m5 = b5; 2350 if (leftright) { 2351 i = 2352 denorm ? be + (Bias + (P-1) - 1 + 1) : 2353 1 + P - bbits; 2354 b2 += i; 2355 s2 += i; 2356 mhi = i2b(1); 2357 if (mhi == NULL) 2358 goto failed_malloc; 2359 } 2360 if (m2 > 0 && s2 > 0) { 2361 i = m2 < s2 ? m2 : s2; 2362 b2 -= i; 2363 m2 -= i; 2364 s2 -= i; 2365 } 2366 if (b5 > 0) { 2367 if (leftright) { 2368 if (m5 > 0) { 2369 mhi = pow5mult(mhi, m5); 2370 if (mhi == NULL) 2371 goto failed_malloc; 2372 b1 = mult(mhi, b); 2373 Bfree(b); 2374 b = b1; 2375 if (b == NULL) 2376 goto failed_malloc; 2377 } 2378 if ((j = b5 - m5)) { 2379 b = pow5mult(b, j); 2380 if (b == NULL) 2381 goto failed_malloc; 2382 } 2383 } 2384 else { 2385 b = pow5mult(b, b5); 2386 if (b == NULL) 2387 goto failed_malloc; 2388 } 2389 } 2390 S = i2b(1); 2391 if (S == NULL) 2392 goto failed_malloc; 2393 if (s5 > 0) { 2394 S = pow5mult(S, s5); 2395 if (S == NULL) 2396 goto failed_malloc; 2397 } 2398 2399 /* Check for special case that d is a normalized power of 2. */ 2400 2401 spec_case = 0; 2402 if ((mode < 2 || leftright) 2403 ) { 2404 if (!word1(&u) && !(word0(&u) & Bndry_mask) 2405 && word0(&u) & (Exp_mask & ~Exp_msk1) 2406 ) { 2407 /* The special case */ 2408 b2 += Log2P; 2409 s2 += Log2P; 2410 spec_case = 1; 2411 } 2412 } 2413 2414 /* Arrange for convenient computation of quotients: 2415 * shift left if necessary so divisor has 4 leading 0 bits. 2416 * 2417 * Perhaps we should just compute leading 28 bits of S once 2418 * and for all and pass them and a shift to quorem, so it 2419 * can do shifts and ors to compute the numerator for q. 2420 */ 2421 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)) 2422 i = 32 - i; 2423#define iInc 28 2424 i = dshift(S, s2); 2425 b2 += i; 2426 m2 += i; 2427 s2 += i; 2428 if (b2 > 0) { 2429 b = lshift(b, b2); 2430 if (b == NULL) 2431 goto failed_malloc; 2432 } 2433 if (s2 > 0) { 2434 S = lshift(S, s2); 2435 if (S == NULL) 2436 goto failed_malloc; 2437 } 2438 if (k_check) { 2439 if (cmp(b,S) < 0) { 2440 k--; 2441 b = multadd(b, 10, 0); /* we botched the k estimate */ 2442 if (b == NULL) 2443 goto failed_malloc; 2444 if (leftright) { 2445 mhi = multadd(mhi, 10, 0); 2446 if (mhi == NULL) 2447 goto failed_malloc; 2448 } 2449 ilim = ilim1; 2450 } 2451 } 2452 if (ilim <= 0 && (mode == 3 || mode == 5)) { 2453 if (ilim < 0) { 2454 /* no digits, fcvt style */ 2455 no_digits: 2456 k = -1 - ndigits; 2457 goto ret; 2458 } 2459 else { 2460 S = multadd(S, 5, 0); 2461 if (S == NULL) 2462 goto failed_malloc; 2463 if (cmp(b, S) <= 0) 2464 goto no_digits; 2465 } 2466 one_digit: 2467 *s++ = '1'; 2468 k++; 2469 goto ret; 2470 } 2471 if (leftright) { 2472 if (m2 > 0) { 2473 mhi = lshift(mhi, m2); 2474 if (mhi == NULL) 2475 goto failed_malloc; 2476 } 2477 2478 /* Compute mlo -- check for special case 2479 * that d is a normalized power of 2. 2480 */ 2481 2482 mlo = mhi; 2483 if (spec_case) { 2484 mhi = Balloc(mhi->k); 2485 if (mhi == NULL) 2486 goto failed_malloc; 2487 Bcopy(mhi, mlo); 2488 mhi = lshift(mhi, Log2P); 2489 if (mhi == NULL) 2490 goto failed_malloc; 2491 } 2492 2493 for(i = 1;;i++) { 2494 dig = quorem(b,S) + '0'; 2495 /* Do we yet have the shortest decimal string 2496 * that will round to d? 2497 */ 2498 j = cmp(b, mlo); 2499 delta = diff(S, mhi); 2500 if (delta == NULL) 2501 goto failed_malloc; 2502 j1 = delta->sign ? 1 : cmp(b, delta); 2503 Bfree(delta); 2504 if (j1 == 0 && mode != 1 && !(word1(&u) & 1) 2505 ) { 2506 if (dig == '9') 2507 goto round_9_up; 2508 if (j > 0) 2509 dig++; 2510 *s++ = dig; 2511 goto ret; 2512 } 2513 if (j < 0 || (j == 0 && mode != 1 2514 && !(word1(&u) & 1) 2515 )) { 2516 if (!b->x[0] && b->wds <= 1) { 2517 goto accept_dig; 2518 } 2519 if (j1 > 0) { 2520 b = lshift(b, 1); 2521 if (b == NULL) 2522 goto failed_malloc; 2523 j1 = cmp(b, S); 2524 if ((j1 > 0 || (j1 == 0 && dig & 1)) 2525 && dig++ == '9') 2526 goto round_9_up; 2527 } 2528 accept_dig: 2529 *s++ = dig; 2530 goto ret; 2531 } 2532 if (j1 > 0) { 2533 if (dig == '9') { /* possible if i == 1 */ 2534 round_9_up: 2535 *s++ = '9'; 2536 goto roundoff; 2537 } 2538 *s++ = dig + 1; 2539 goto ret; 2540 } 2541 *s++ = dig; 2542 if (i == ilim) 2543 break; 2544 b = multadd(b, 10, 0); 2545 if (b == NULL) 2546 goto failed_malloc; 2547 if (mlo == mhi) { 2548 mlo = mhi = multadd(mhi, 10, 0); 2549 if (mlo == NULL) 2550 goto failed_malloc; 2551 } 2552 else { 2553 mlo = multadd(mlo, 10, 0); 2554 if (mlo == NULL) 2555 goto failed_malloc; 2556 mhi = multadd(mhi, 10, 0); 2557 if (mhi == NULL) 2558 goto failed_malloc; 2559 } 2560 } 2561 } 2562 else 2563 for(i = 1;; i++) { 2564 *s++ = dig = quorem(b,S) + '0'; 2565 if (!b->x[0] && b->wds <= 1) { 2566 goto ret; 2567 } 2568 if (i >= ilim) 2569 break; 2570 b = multadd(b, 10, 0); 2571 if (b == NULL) 2572 goto failed_malloc; 2573 } 2574 2575 /* Round off last digit */ 2576 2577 b = lshift(b, 1); 2578 if (b == NULL) 2579 goto failed_malloc; 2580 j = cmp(b, S); 2581 if (j > 0 || (j == 0 && dig & 1)) { 2582 roundoff: 2583 while(*--s == '9') 2584 if (s == s0) { 2585 k++; 2586 *s++ = '1'; 2587 goto ret; 2588 } 2589 ++*s++; 2590 } 2591 else { 2592 while(*--s == '0'); 2593 s++; 2594 } 2595 ret: 2596 Bfree(S); 2597 if (mhi) { 2598 if (mlo && mlo != mhi) 2599 Bfree(mlo); 2600 Bfree(mhi); 2601 } 2602 ret1: 2603 Bfree(b); 2604 *s = 0; 2605 *decpt = k + 1; 2606 if (rve) 2607 *rve = s; 2608 return s0; 2609 failed_malloc: 2610 if (S) 2611 Bfree(S); 2612 if (mlo && mlo != mhi) 2613 Bfree(mlo); 2614 if (mhi) 2615 Bfree(mhi); 2616 if (b) 2617 Bfree(b); 2618 if (s0) 2619 _Py_dg_freedtoa(s0); 2620 return NULL; 2621} 2622#ifdef __cplusplus 2623} 2624#endif 2625 2626#endif /* PY_NO_SHORT_FLOAT_REPR */ 2627