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