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