/rpython/translator/c/src/dtoa.c
C | 2976 lines | 2123 code | 236 blank | 617 comment | 610 complexity | 2a066e10378be37507d34f354cd9073d MD5 | raw file
Large files files are truncated, but you can click here to view the full file
- /****************************************************************
- *
- * The author of this software is David M. Gay.
- *
- * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
- *
- * Permission to use, copy, modify, and distribute this software for any
- * purpose without fee is hereby granted, provided that this entire notice
- * is included in all copies of any software which is or includes a copy
- * or modification of this software and in all copies of the supporting
- * documentation for such software.
- *
- * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
- * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
- * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
- * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
- *
- ***************************************************************/
- /****************************************************************
- * This is dtoa.c by David M. Gay, downloaded from
- * http://www.netlib.org/fp/dtoa.c on April 15, 2009 and modified for
- * inclusion into the Python core by Mark E. T. Dickinson and Eric V. Smith.
- *
- * Please remember to check http://www.netlib.org/fp regularly (and especially
- * before any Python release) for bugfixes and updates.
- *
- * The major modifications from Gay's original code are as follows:
- *
- * 0. The original code has been specialized to Python's needs by removing
- * many of the #ifdef'd sections. In particular, code to support VAX and
- * IBM floating-point formats, hex NaNs, hex floats, locale-aware
- * treatment of the decimal point, and setting of the inexact flag have
- * been removed.
- *
- * 1. We use PyMem_Malloc and PyMem_Free in place of malloc and free.
- *
- * 2. The public functions strtod, dtoa and freedtoa all now have
- * a _Py_dg_ prefix.
- *
- * 3. Instead of assuming that PyMem_Malloc always succeeds, we thread
- * PyMem_Malloc failures through the code. The functions
- *
- * Balloc, multadd, s2b, i2b, mult, pow5mult, lshift, diff, d2b
- *
- * of return type *Bigint all return NULL to indicate a malloc failure.
- * Similarly, rv_alloc and nrv_alloc (return type char *) return NULL on
- * failure. bigcomp now has return type int (it used to be void) and
- * returns -1 on failure and 0 otherwise. __Py_dg_dtoa returns NULL
- * on failure. __Py_dg_strtod indicates failure due to malloc failure
- * by returning -1.0, setting errno=ENOMEM and *se to s00.
- *
- * 4. The static variable dtoa_result has been removed. Callers of
- * __Py_dg_dtoa are expected to call __Py_dg_freedtoa to free
- * the memory allocated by __Py_dg_dtoa.
- *
- * 5. The code has been reformatted to better fit with Python's
- * C style guide (PEP 7).
- *
- * 6. A bug in the memory allocation has been fixed: to avoid FREEing memory
- * that hasn't been MALLOC'ed, private_mem should only be used when k <=
- * Kmax.
- *
- * 7. __Py_dg_strtod has been modified so that it doesn't accept strings with
- * leading whitespace.
- *
- ***************************************************************/
- /* Please send bug reports for the original dtoa.c code to David M. Gay (dmg
- * at acm dot org, with " at " changed at "@" and " dot " changed to ".").
- * Please report bugs for this modified version using the Python issue tracker
- * (http://bugs.python.org). */
- /* On a machine with IEEE extended-precision registers, it is
- * necessary to specify double-precision (53-bit) rounding precision
- * before invoking strtod or dtoa. If the machine uses (the equivalent
- * of) Intel 80x87 arithmetic, the call
- * _control87(PC_53, MCW_PC);
- * does this with many compilers. Whether this or another call is
- * appropriate depends on the compiler; for this to work, it may be
- * necessary to #include "float.h" or another system-dependent header
- * file.
- */
- /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
- *
- * This strtod returns a nearest machine number to the input decimal
- * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
- * broken by the IEEE round-even rule. Otherwise ties are broken by
- * biased rounding (add half and chop).
- *
- * Inspired loosely by William D. Clinger's paper "How to Read Floating
- * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
- *
- * Modifications:
- *
- * 1. We only require IEEE, IBM, or VAX double-precision
- * arithmetic (not IEEE double-extended).
- * 2. We get by with floating-point arithmetic in a case that
- * Clinger missed -- when we're computing d * 10^n
- * for a small integer d and the integer n is not too
- * much larger than 22 (the maximum integer k for which
- * we can represent 10^k exactly), we may be able to
- * compute (d*10^k) * 10^(e-k) with just one roundoff.
- * 3. Rather than a bit-at-a-time adjustment of the binary
- * result in the hard case, we use floating-point
- * arithmetic to determine the adjustment to within
- * one bit; only in really hard cases do we need to
- * compute a second residual.
- * 4. Because of 3., we don't need a large table of powers of 10
- * for ten-to-e (just some small tables, e.g. of 10^k
- * for 0 <= k <= 22).
- */
- /* Linking of Python's #defines to Gay's #defines starts here. */
- /* Begin PYPY hacks */
- /* #include "Python.h" */
- #define HAVE_UINT32_T
- #define HAVE_INT32_T
- #define HAVE_UINT64_T
- #define PY_UINT32_T unsigned int
- #define PY_INT32_T int
- #define PY_UINT64_T unsigned long long
- #include <stdlib.h>
- #include <errno.h>
- #include <assert.h>
- #include <stdio.h>
- #include <string.h>
- #include "src/asm.h"
- #define PyMem_Malloc malloc
- #define PyMem_Free free
- /* End PYPY hacks */
- /* if PY_NO_SHORT_FLOAT_REPR is defined, then don't even try to compile
- the following code */
- #ifndef PY_NO_SHORT_FLOAT_REPR
- #include <float.h>
- #define MALLOC PyMem_Malloc
- #define FREE PyMem_Free
- /* This code should also work for ARM mixed-endian format on little-endian
- machines, where doubles have byte order 45670123 (in increasing address
- order, 0 being the least significant byte). */
- #ifdef DOUBLE_IS_LITTLE_ENDIAN_IEEE754
- # define IEEE_8087
- #endif
- #if defined(DOUBLE_IS_BIG_ENDIAN_IEEE754) || \
- defined(DOUBLE_IS_ARM_MIXED_ENDIAN_IEEE754)
- # define IEEE_MC68k
- #endif
- #if defined(IEEE_8087) + defined(IEEE_MC68k) != 1
- #error "Exactly one of IEEE_8087 or IEEE_MC68k should be defined."
- #endif
- /* The code below assumes that the endianness of integers matches the
- endianness of the two 32-bit words of a double. Check this. */
- #if defined(WORDS_BIGENDIAN) && (defined(DOUBLE_IS_LITTLE_ENDIAN_IEEE754) || \
- defined(DOUBLE_IS_ARM_MIXED_ENDIAN_IEEE754))
- #error "doubles and ints have incompatible endianness"
- #endif
- #if !defined(WORDS_BIGENDIAN) && defined(DOUBLE_IS_BIG_ENDIAN_IEEE754)
- #error "doubles and ints have incompatible endianness"
- #endif
- #if defined(HAVE_UINT32_T) && defined(HAVE_INT32_T)
- typedef PY_UINT32_T ULong;
- typedef PY_INT32_T Long;
- #else
- #error "Failed to find an exact-width 32-bit integer type"
- #endif
- #if defined(HAVE_UINT64_T)
- #define ULLong PY_UINT64_T
- #else
- #undef ULLong
- #endif
- #undef DEBUG
- #ifdef Py_DEBUG
- #define DEBUG
- #endif
- /* End Python #define linking */
- #ifdef DEBUG
- #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
- #endif
- #ifndef PRIVATE_MEM
- #define PRIVATE_MEM 2304
- #endif
- #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
- static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
- #ifdef __cplusplus
- extern "C" {
- #endif
- typedef union { double d; ULong L[2]; } U;
- #ifdef IEEE_8087
- #define word0(x) (x)->L[1]
- #define word1(x) (x)->L[0]
- #else
- #define word0(x) (x)->L[0]
- #define word1(x) (x)->L[1]
- #endif
- #define dval(x) (x)->d
- #ifndef STRTOD_DIGLIM
- #define STRTOD_DIGLIM 40
- #endif
- /* maximum permitted exponent value for strtod; exponents larger than
- MAX_ABS_EXP in absolute value get truncated to +-MAX_ABS_EXP. MAX_ABS_EXP
- should fit into an int. */
- #ifndef MAX_ABS_EXP
- #define MAX_ABS_EXP 19999U
- #endif
- /* The following definition of Storeinc is appropriate for MIPS processors.
- * An alternative that might be better on some machines is
- * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
- */
- #if defined(IEEE_8087)
- #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
- ((unsigned short *)a)[0] = (unsigned short)c, a++)
- #else
- #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
- ((unsigned short *)a)[1] = (unsigned short)c, a++)
- #endif
- /* #define P DBL_MANT_DIG */
- /* Ten_pmax = floor(P*log(2)/log(5)) */
- /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
- /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
- /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
- #define Exp_shift 20
- #define Exp_shift1 20
- #define Exp_msk1 0x100000
- #define Exp_msk11 0x100000
- #define Exp_mask 0x7ff00000
- #define P 53
- #define Nbits 53
- #define Bias 1023
- #define Emax 1023
- #define Emin (-1022)
- #define Etiny (-1074) /* smallest denormal is 2**Etiny */
- #define Exp_1 0x3ff00000
- #define Exp_11 0x3ff00000
- #define Ebits 11
- #define Frac_mask 0xfffff
- #define Frac_mask1 0xfffff
- #define Ten_pmax 22
- #define Bletch 0x10
- #define Bndry_mask 0xfffff
- #define Bndry_mask1 0xfffff
- #define Sign_bit 0x80000000
- #define Log2P 1
- #define Tiny0 0
- #define Tiny1 1
- #define Quick_max 14
- #define Int_max 14
- #ifndef Flt_Rounds
- #ifdef FLT_ROUNDS
- #define Flt_Rounds FLT_ROUNDS
- #else
- #define Flt_Rounds 1
- #endif
- #endif /*Flt_Rounds*/
- #define Rounding Flt_Rounds
- #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
- #define Big1 0xffffffff
- /* struct BCinfo is used to pass information from __Py_dg_strtod to bigcomp */
- typedef struct BCinfo BCinfo;
- struct
- BCinfo {
- int e0, nd, nd0, scale;
- };
- #define FFFFFFFF 0xffffffffUL
- #define Kmax 7
- /* struct Bigint is used to represent arbitrary-precision integers. These
- integers are stored in sign-magnitude format, with the magnitude stored as
- an array of base 2**32 digits. Bigints are always normalized: if x is a
- Bigint then x->wds >= 1, and either x->wds == 1 or x[wds-1] is nonzero.
- The Bigint fields are as follows:
- - next is a header used by Balloc and Bfree to keep track of lists
- of freed Bigints; it's also used for the linked list of
- powers of 5 of the form 5**2**i used by pow5mult.
- - k indicates which pool this Bigint was allocated from
- - maxwds is the maximum number of words space was allocated for
- (usually maxwds == 2**k)
- - sign is 1 for negative Bigints, 0 for positive. The sign is unused
- (ignored on inputs, set to 0 on outputs) in almost all operations
- involving Bigints: a notable exception is the diff function, which
- ignores signs on inputs but sets the sign of the output correctly.
- - wds is the actual number of significant words
- - x contains the vector of words (digits) for this Bigint, from least
- significant (x[0]) to most significant (x[wds-1]).
- */
- struct
- Bigint {
- struct Bigint *next;
- int k, maxwds, sign, wds;
- ULong x[1];
- };
- typedef struct Bigint Bigint;
- #ifndef Py_USING_MEMORY_DEBUGGER
- /* Memory management: memory is allocated from, and returned to, Kmax+1 pools
- of memory, where pool k (0 <= k <= Kmax) is for Bigints b with b->maxwds ==
- 1 << k. These pools are maintained as linked lists, with freelist[k]
- pointing to the head of the list for pool k.
- On allocation, if there's no free slot in the appropriate pool, MALLOC is
- called to get more memory. This memory is not returned to the system until
- Python quits. There's also a private memory pool that's allocated from
- in preference to using MALLOC.
- For Bigints with more than (1 << Kmax) digits (which implies at least 1233
- decimal digits), memory is directly allocated using MALLOC, and freed using
- FREE.
- XXX: it would be easy to bypass this memory-management system and
- translate each call to Balloc into a call to PyMem_Malloc, and each
- Bfree to PyMem_Free. Investigate whether this has any significant
- performance on impact. */
- static Bigint *freelist[Kmax+1];
- /* Allocate space for a Bigint with up to 1<<k digits */
- static Bigint *
- Balloc(int k)
- {
- int x;
- Bigint *rv;
- unsigned int len;
- if (k <= Kmax && (rv = freelist[k]))
- freelist[k] = rv->next;
- else {
- x = 1 << k;
- len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
- /sizeof(double);
- if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
- rv = (Bigint*)pmem_next;
- pmem_next += len;
- }
- else {
- rv = (Bigint*)MALLOC(len*sizeof(double));
- if (rv == NULL)
- return NULL;
- }
- rv->k = k;
- rv->maxwds = x;
- }
- rv->sign = rv->wds = 0;
- return rv;
- }
- /* Free a Bigint allocated with Balloc */
- static void
- Bfree(Bigint *v)
- {
- if (v) {
- if (v->k > Kmax)
- FREE((void*)v);
- else {
- v->next = freelist[v->k];
- freelist[v->k] = v;
- }
- }
- }
- #else
- /* Alternative versions of Balloc and Bfree that use PyMem_Malloc and
- PyMem_Free directly in place of the custom memory allocation scheme above.
- These are provided for the benefit of memory debugging tools like
- Valgrind. */
- /* Allocate space for a Bigint with up to 1<<k digits */
- static Bigint *
- Balloc(int k)
- {
- int x;
- Bigint *rv;
- unsigned int len;
- x = 1 << k;
- len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
- /sizeof(double);
- rv = (Bigint*)MALLOC(len*sizeof(double));
- if (rv == NULL)
- return NULL;
- rv->k = k;
- rv->maxwds = x;
- rv->sign = rv->wds = 0;
- return rv;
- }
- /* Free a Bigint allocated with Balloc */
- static void
- Bfree(Bigint *v)
- {
- if (v) {
- FREE((void*)v);
- }
- }
- #endif /* Py_USING_MEMORY_DEBUGGER */
- #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
- y->wds*sizeof(Long) + 2*sizeof(int))
- /* Multiply a Bigint b by m and add a. Either modifies b in place and returns
- a pointer to the modified b, or Bfrees b and returns a pointer to a copy.
- On failure, return NULL. In this case, b will have been already freed. */
- static Bigint *
- multadd(Bigint *b, int m, int a) /* multiply by m and add a */
- {
- int i, wds;
- #ifdef ULLong
- ULong *x;
- ULLong carry, y;
- #else
- ULong carry, *x, y;
- ULong xi, z;
- #endif
- Bigint *b1;
- wds = b->wds;
- x = b->x;
- i = 0;
- carry = a;
- do {
- #ifdef ULLong
- y = *x * (ULLong)m + carry;
- carry = y >> 32;
- *x++ = (ULong)(y & FFFFFFFF);
- #else
- xi = *x;
- y = (xi & 0xffff) * m + carry;
- z = (xi >> 16) * m + (y >> 16);
- carry = z >> 16;
- *x++ = (z << 16) + (y & 0xffff);
- #endif
- }
- while(++i < wds);
- if (carry) {
- if (wds >= b->maxwds) {
- b1 = Balloc(b->k+1);
- if (b1 == NULL){
- Bfree(b);
- return NULL;
- }
- Bcopy(b1, b);
- Bfree(b);
- b = b1;
- }
- b->x[wds++] = (ULong)carry;
- b->wds = wds;
- }
- return b;
- }
- /* convert a string s containing nd decimal digits (possibly containing a
- decimal separator at position nd0, which is ignored) to a Bigint. This
- function carries on where the parsing code in __Py_dg_strtod leaves off: on
- entry, y9 contains the result of converting the first 9 digits. Returns
- NULL on failure. */
- static Bigint *
- s2b(const char *s, int nd0, int nd, ULong y9)
- {
- Bigint *b;
- int i, k;
- Long x, y;
- x = (nd + 8) / 9;
- for(k = 0, y = 1; x > y; y <<= 1, k++) ;
- b = Balloc(k);
- if (b == NULL)
- return NULL;
- b->x[0] = y9;
- b->wds = 1;
- if (nd <= 9)
- return b;
- s += 9;
- for (i = 9; i < nd0; i++) {
- b = multadd(b, 10, *s++ - '0');
- if (b == NULL)
- return NULL;
- }
- s++;
- for(; i < nd; i++) {
- b = multadd(b, 10, *s++ - '0');
- if (b == NULL)
- return NULL;
- }
- return b;
- }
- /* count leading 0 bits in the 32-bit integer x. */
- static int
- hi0bits(ULong x)
- {
- int k = 0;
- if (!(x & 0xffff0000)) {
- k = 16;
- x <<= 16;
- }
- if (!(x & 0xff000000)) {
- k += 8;
- x <<= 8;
- }
- if (!(x & 0xf0000000)) {
- k += 4;
- x <<= 4;
- }
- if (!(x & 0xc0000000)) {
- k += 2;
- x <<= 2;
- }
- if (!(x & 0x80000000)) {
- k++;
- if (!(x & 0x40000000))
- return 32;
- }
- return k;
- }
- /* count trailing 0 bits in the 32-bit integer y, and shift y right by that
- number of bits. */
- static int
- lo0bits(ULong *y)
- {
- int k;
- ULong x = *y;
- if (x & 7) {
- if (x & 1)
- return 0;
- if (x & 2) {
- *y = x >> 1;
- return 1;
- }
- *y = x >> 2;
- return 2;
- }
- k = 0;
- if (!(x & 0xffff)) {
- k = 16;
- x >>= 16;
- }
- if (!(x & 0xff)) {
- k += 8;
- x >>= 8;
- }
- if (!(x & 0xf)) {
- k += 4;
- x >>= 4;
- }
- if (!(x & 0x3)) {
- k += 2;
- x >>= 2;
- }
- if (!(x & 1)) {
- k++;
- x >>= 1;
- if (!x)
- return 32;
- }
- *y = x;
- return k;
- }
- /* convert a small nonnegative integer to a Bigint */
- static Bigint *
- i2b(int i)
- {
- Bigint *b;
- b = Balloc(1);
- if (b == NULL)
- return NULL;
- b->x[0] = i;
- b->wds = 1;
- return b;
- }
- /* multiply two Bigints. Returns a new Bigint, or NULL on failure. Ignores
- the signs of a and b. */
- static Bigint *
- mult(Bigint *a, Bigint *b)
- {
- Bigint *c;
- int k, wa, wb, wc;
- ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
- ULong y;
- #ifdef ULLong
- ULLong carry, z;
- #else
- ULong carry, z;
- ULong z2;
- #endif
- if ((!a->x[0] && a->wds == 1) || (!b->x[0] && b->wds == 1)) {
- c = Balloc(0);
- if (c == NULL)
- return NULL;
- c->wds = 1;
- c->x[0] = 0;
- return c;
- }
- if (a->wds < b->wds) {
- c = a;
- a = b;
- b = c;
- }
- k = a->k;
- wa = a->wds;
- wb = b->wds;
- wc = wa + wb;
- if (wc > a->maxwds)
- k++;
- c = Balloc(k);
- if (c == NULL)
- return NULL;
- for(x = c->x, xa = x + wc; x < xa; x++)
- *x = 0;
- xa = a->x;
- xae = xa + wa;
- xb = b->x;
- xbe = xb + wb;
- xc0 = c->x;
- #ifdef ULLong
- for(; xb < xbe; xc0++) {
- if ((y = *xb++)) {
- x = xa;
- xc = xc0;
- carry = 0;
- do {
- z = *x++ * (ULLong)y + *xc + carry;
- carry = z >> 32;
- *xc++ = (ULong)(z & FFFFFFFF);
- }
- while(x < xae);
- *xc = (ULong)carry;
- }
- }
- #else
- for(; xb < xbe; xb++, xc0++) {
- if (y = *xb & 0xffff) {
- x = xa;
- xc = xc0;
- carry = 0;
- do {
- z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
- carry = z >> 16;
- z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
- carry = z2 >> 16;
- Storeinc(xc, z2, z);
- }
- while(x < xae);
- *xc = carry;
- }
- if (y = *xb >> 16) {
- x = xa;
- xc = xc0;
- carry = 0;
- z2 = *xc;
- do {
- z = (*x & 0xffff) * y + (*xc >> 16) + carry;
- carry = z >> 16;
- Storeinc(xc, z, z2);
- z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
- carry = z2 >> 16;
- }
- while(x < xae);
- *xc = z2;
- }
- }
- #endif
- for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
- c->wds = wc;
- return c;
- }
- #ifndef Py_USING_MEMORY_DEBUGGER
- /* p5s is a linked list of powers of 5 of the form 5**(2**i), i >= 2 */
- static Bigint *p5s;
- /* multiply the Bigint b by 5**k. Returns a pointer to the result, or NULL on
- failure; if the returned pointer is distinct from b then the original
- Bigint b will have been Bfree'd. Ignores the sign of b. */
- static Bigint *
- pow5mult(Bigint *b, int k)
- {
- Bigint *b1, *p5, *p51;
- int i;
- static int p05[3] = { 5, 25, 125 };
- if ((i = k & 3)) {
- b = multadd(b, p05[i-1], 0);
- if (b == NULL)
- return NULL;
- }
- if (!(k >>= 2))
- return b;
- p5 = p5s;
- if (!p5) {
- /* first time */
- p5 = i2b(625);
- if (p5 == NULL) {
- Bfree(b);
- return NULL;
- }
- p5s = p5;
- p5->next = 0;
- }
- for(;;) {
- if (k & 1) {
- b1 = mult(b, p5);
- Bfree(b);
- b = b1;
- if (b == NULL)
- return NULL;
- }
- if (!(k >>= 1))
- break;
- p51 = p5->next;
- if (!p51) {
- p51 = mult(p5,p5);
- if (p51 == NULL) {
- Bfree(b);
- return NULL;
- }
- p51->next = 0;
- p5->next = p51;
- }
- p5 = p51;
- }
- return b;
- }
- #else
- /* Version of pow5mult that doesn't cache powers of 5. Provided for
- the benefit of memory debugging tools like Valgrind. */
- static Bigint *
- pow5mult(Bigint *b, int k)
- {
- Bigint *b1, *p5, *p51;
- int i;
- static int p05[3] = { 5, 25, 125 };
- if ((i = k & 3)) {
- b = multadd(b, p05[i-1], 0);
- if (b == NULL)
- return NULL;
- }
- if (!(k >>= 2))
- return b;
- p5 = i2b(625);
- if (p5 == NULL) {
- Bfree(b);
- return NULL;
- }
- for(;;) {
- if (k & 1) {
- b1 = mult(b, p5);
- Bfree(b);
- b = b1;
- if (b == NULL) {
- Bfree(p5);
- return NULL;
- }
- }
- if (!(k >>= 1))
- break;
- p51 = mult(p5, p5);
- Bfree(p5);
- p5 = p51;
- if (p5 == NULL) {
- Bfree(b);
- return NULL;
- }
- }
- Bfree(p5);
- return b;
- }
- #endif /* Py_USING_MEMORY_DEBUGGER */
- /* shift a Bigint b left by k bits. Return a pointer to the shifted result,
- or NULL on failure. If the returned pointer is distinct from b then the
- original b will have been Bfree'd. Ignores the sign of b. */
- static Bigint *
- lshift(Bigint *b, int k)
- {
- int i, k1, n, n1;
- Bigint *b1;
- ULong *x, *x1, *xe, z;
- if (!k || (!b->x[0] && b->wds == 1))
- return b;
- n = k >> 5;
- k1 = b->k;
- n1 = n + b->wds + 1;
- for(i = b->maxwds; n1 > i; i <<= 1)
- k1++;
- b1 = Balloc(k1);
- if (b1 == NULL) {
- Bfree(b);
- return NULL;
- }
- x1 = b1->x;
- for(i = 0; i < n; i++)
- *x1++ = 0;
- x = b->x;
- xe = x + b->wds;
- if (k &= 0x1f) {
- k1 = 32 - k;
- z = 0;
- do {
- *x1++ = *x << k | z;
- z = *x++ >> k1;
- }
- while(x < xe);
- if ((*x1 = z))
- ++n1;
- }
- else do
- *x1++ = *x++;
- while(x < xe);
- b1->wds = n1 - 1;
- Bfree(b);
- return b1;
- }
- /* Do a three-way compare of a and b, returning -1 if a < b, 0 if a == b and
- 1 if a > b. Ignores signs of a and b. */
- static int
- cmp(Bigint *a, Bigint *b)
- {
- ULong *xa, *xa0, *xb, *xb0;
- int i, j;
- i = a->wds;
- j = b->wds;
- #ifdef DEBUG
- if (i > 1 && !a->x[i-1])
- Bug("cmp called with a->x[a->wds-1] == 0");
- if (j > 1 && !b->x[j-1])
- Bug("cmp called with b->x[b->wds-1] == 0");
- #endif
- if (i -= j)
- return i;
- xa0 = a->x;
- xa = xa0 + j;
- xb0 = b->x;
- xb = xb0 + j;
- for(;;) {
- if (*--xa != *--xb)
- return *xa < *xb ? -1 : 1;
- if (xa <= xa0)
- break;
- }
- return 0;
- }
- /* Take the difference of Bigints a and b, returning a new Bigint. Returns
- NULL on failure. The signs of a and b are ignored, but the sign of the
- result is set appropriately. */
- static Bigint *
- diff(Bigint *a, Bigint *b)
- {
- Bigint *c;
- int i, wa, wb;
- ULong *xa, *xae, *xb, *xbe, *xc;
- #ifdef ULLong
- ULLong borrow, y;
- #else
- ULong borrow, y;
- ULong z;
- #endif
- i = cmp(a,b);
- if (!i) {
- c = Balloc(0);
- if (c == NULL)
- return NULL;
- c->wds = 1;
- c->x[0] = 0;
- return c;
- }
- if (i < 0) {
- c = a;
- a = b;
- b = c;
- i = 1;
- }
- else
- i = 0;
- c = Balloc(a->k);
- if (c == NULL)
- return NULL;
- c->sign = i;
- wa = a->wds;
- xa = a->x;
- xae = xa + wa;
- wb = b->wds;
- xb = b->x;
- xbe = xb + wb;
- xc = c->x;
- borrow = 0;
- #ifdef ULLong
- do {
- y = (ULLong)*xa++ - *xb++ - borrow;
- borrow = y >> 32 & (ULong)1;
- *xc++ = (ULong)(y & FFFFFFFF);
- }
- while(xb < xbe);
- while(xa < xae) {
- y = *xa++ - borrow;
- borrow = y >> 32 & (ULong)1;
- *xc++ = (ULong)(y & FFFFFFFF);
- }
- #else
- do {
- y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
- borrow = (y & 0x10000) >> 16;
- z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
- borrow = (z & 0x10000) >> 16;
- Storeinc(xc, z, y);
- }
- while(xb < xbe);
- while(xa < xae) {
- y = (*xa & 0xffff) - borrow;
- borrow = (y & 0x10000) >> 16;
- z = (*xa++ >> 16) - borrow;
- borrow = (z & 0x10000) >> 16;
- Storeinc(xc, z, y);
- }
- #endif
- while(!*--xc)
- wa--;
- c->wds = wa;
- return c;
- }
- /* Given a positive normal double x, return the difference between x and the
- next double up. Doesn't give correct results for subnormals. */
- static double
- ulp(U *x)
- {
- Long L;
- U u;
- L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
- word0(&u) = L;
- word1(&u) = 0;
- return dval(&u);
- }
- /* Convert a Bigint to a double plus an exponent */
- static double
- b2d(Bigint *a, int *e)
- {
- ULong *xa, *xa0, w, y, z;
- int k;
- U d;
- xa0 = a->x;
- xa = xa0 + a->wds;
- y = *--xa;
- #ifdef DEBUG
- if (!y) Bug("zero y in b2d");
- #endif
- k = hi0bits(y);
- *e = 32 - k;
- if (k < Ebits) {
- word0(&d) = Exp_1 | y >> (Ebits - k);
- w = xa > xa0 ? *--xa : 0;
- word1(&d) = y << ((32-Ebits) + k) | w >> (Ebits - k);
- goto ret_d;
- }
- z = xa > xa0 ? *--xa : 0;
- if (k -= Ebits) {
- word0(&d) = Exp_1 | y << k | z >> (32 - k);
- y = xa > xa0 ? *--xa : 0;
- word1(&d) = z << k | y >> (32 - k);
- }
- else {
- word0(&d) = Exp_1 | y;
- word1(&d) = z;
- }
- ret_d:
- return dval(&d);
- }
- /* Convert a scaled double to a Bigint plus an exponent. Similar to d2b,
- except that it accepts the scale parameter used in __Py_dg_strtod (which
- should be either 0 or 2*P), and the normalization for the return value is
- different (see below). On input, d should be finite and nonnegative, and d
- / 2**scale should be exactly representable as an IEEE 754 double.
- Returns a Bigint b and an integer e such that
- dval(d) / 2**scale = b * 2**e.
- Unlike d2b, b is not necessarily odd: b and e are normalized so
- that either 2**(P-1) <= b < 2**P and e >= Etiny, or b < 2**P
- and e == Etiny. This applies equally to an input of 0.0: in that
- case the return values are b = 0 and e = Etiny.
- The above normalization ensures that for all possible inputs d,
- 2**e gives ulp(d/2**scale).
- Returns NULL on failure.
- */
- static Bigint *
- sd2b(U *d, int scale, int *e)
- {
- Bigint *b;
- b = Balloc(1);
- if (b == NULL)
- return NULL;
-
- /* First construct b and e assuming that scale == 0. */
- b->wds = 2;
- b->x[0] = word1(d);
- b->x[1] = word0(d) & Frac_mask;
- *e = Etiny - 1 + (int)((word0(d) & Exp_mask) >> Exp_shift);
- if (*e < Etiny)
- *e = Etiny;
- else
- b->x[1] |= Exp_msk1;
- /* Now adjust for scale, provided that b != 0. */
- if (scale && (b->x[0] || b->x[1])) {
- *e -= scale;
- if (*e < Etiny) {
- scale = Etiny - *e;
- *e = Etiny;
- /* We can't shift more than P-1 bits without shifting out a 1. */
- assert(0 < scale && scale <= P - 1);
- if (scale >= 32) {
- /* The bits shifted out should all be zero. */
- assert(b->x[0] == 0);
- b->x[0] = b->x[1];
- b->x[1] = 0;
- scale -= 32;
- }
- if (scale) {
- /* The bits shifted out should all be zero. */
- assert(b->x[0] << (32 - scale) == 0);
- b->x[0] = (b->x[0] >> scale) | (b->x[1] << (32 - scale));
- b->x[1] >>= scale;
- }
- }
- }
- /* Ensure b is normalized. */
- if (!b->x[1])
- b->wds = 1;
- return b;
- }
- /* Convert a double to a Bigint plus an exponent. Return NULL on failure.
- Given a finite nonzero double d, return an odd Bigint b and exponent *e
- such that fabs(d) = b * 2**e. On return, *bbits gives the number of
- significant bits of b; that is, 2**(*bbits-1) <= b < 2**(*bbits).
- If d is zero, then b == 0, *e == -1010, *bbits = 0.
- */
- static Bigint *
- d2b(U *d, int *e, int *bits)
- {
- Bigint *b;
- int de, k;
- ULong *x, y, z;
- int i;
- b = Balloc(1);
- if (b == NULL)
- return NULL;
- x = b->x;
- z = word0(d) & Frac_mask;
- word0(d) &= 0x7fffffff; /* clear sign bit, which we ignore */
- if ((de = (int)(word0(d) >> Exp_shift)))
- z |= Exp_msk1;
- if ((y = word1(d))) {
- if ((k = lo0bits(&y))) {
- x[0] = y | z << (32 - k);
- z >>= k;
- }
- else
- x[0] = y;
- i =
- b->wds = (x[1] = z) ? 2 : 1;
- }
- else {
- k = lo0bits(&z);
- x[0] = z;
- i =
- b->wds = 1;
- k += 32;
- }
- if (de) {
- *e = de - Bias - (P-1) + k;
- *bits = P - k;
- }
- else {
- *e = de - Bias - (P-1) + 1 + k;
- *bits = 32*i - hi0bits(x[i-1]);
- }
- return b;
- }
- /* Compute the ratio of two Bigints, as a double. The result may have an
- error of up to 2.5 ulps. */
- static double
- ratio(Bigint *a, Bigint *b)
- {
- U da, db;
- int k, ka, kb;
- dval(&da) = b2d(a, &ka);
- dval(&db) = b2d(b, &kb);
- k = ka - kb + 32*(a->wds - b->wds);
- if (k > 0)
- word0(&da) += k*Exp_msk1;
- else {
- k = -k;
- word0(&db) += k*Exp_msk1;
- }
- return dval(&da) / dval(&db);
- }
- static const double
- tens[] = {
- 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
- 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
- 1e20, 1e21, 1e22
- };
- static const double
- bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
- static const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
- 9007199254740992.*9007199254740992.e-256
- /* = 2^106 * 1e-256 */
- };
- /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
- /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
- #define Scale_Bit 0x10
- #define n_bigtens 5
- #define ULbits 32
- #define kshift 5
- #define kmask 31
- static int
- dshift(Bigint *b, int p2)
- {
- int rv = hi0bits(b->x[b->wds-1]) - 4;
- if (p2 > 0)
- rv -= p2;
- return rv & kmask;
- }
- /* special case of Bigint division. The quotient is always in the range 0 <=
- quotient < 10, and on entry the divisor S is normalized so that its top 4
- bits (28--31) are zero and bit 27 is set. */
- static int
- quorem(Bigint *b, Bigint *S)
- {
- int n;
- ULong *bx, *bxe, q, *sx, *sxe;
- #ifdef ULLong
- ULLong borrow, carry, y, ys;
- #else
- ULong borrow, carry, y, ys;
- ULong si, z, zs;
- #endif
- n = S->wds;
- #ifdef DEBUG
- /*debug*/ if (b->wds > n)
- /*debug*/ Bug("oversize b in quorem");
- #endif
- if (b->wds < n)
- return 0;
- sx = S->x;
- sxe = sx + --n;
- bx = b->x;
- bxe = bx + n;
- q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
- #ifdef DEBUG
- /*debug*/ if (q > 9)
- /*debug*/ Bug("oversized quotient in quorem");
- #endif
- if (q) {
- borrow = 0;
- carry = 0;
- do {
- #ifdef ULLong
- ys = *sx++ * (ULLong)q + carry;
- carry = ys >> 32;
- y = *bx - (ys & FFFFFFFF) - borrow;
- borrow = y >> 32 & (ULong)1;
- *bx++ = (ULong)(y & FFFFFFFF);
- #else
- si = *sx++;
- ys = (si & 0xffff) * q + carry;
- zs = (si >> 16) * q + (ys >> 16);
- carry = zs >> 16;
- y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
- borrow = (y & 0x10000) >> 16;
- z = (*bx >> 16) - (zs & 0xffff) - borrow;
- borrow = (z & 0x10000) >> 16;
- Storeinc(bx, z, y);
- #endif
- }
- while(sx <= sxe);
- if (!*bxe) {
- bx = b->x;
- while(--bxe > bx && !*bxe)
- --n;
- b->wds = n;
- }
- }
- if (cmp(b, S) >= 0) {
- q++;
- borrow = 0;
- carry = 0;
- bx = b->x;
- sx = S->x;
- do {
- #ifdef ULLong
- ys = *sx++ + carry;
- carry = ys >> 32;
- y = *bx - (ys & FFFFFFFF) - borrow;
- borrow = y >> 32 & (ULong)1;
- *bx++ = (ULong)(y & FFFFFFFF);
- #else
- si = *sx++;
- ys = (si & 0xffff) + carry;
- zs = (si >> 16) + (ys >> 16);
- carry = zs >> 16;
- y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
- borrow = (y & 0x10000) >> 16;
- z = (*bx >> 16) - (zs & 0xffff) - borrow;
- borrow = (z & 0x10000) >> 16;
- Storeinc(bx, z, y);
- #endif
- }
- while(sx <= sxe);
- bx = b->x;
- bxe = bx + n;
- if (!*bxe) {
- while(--bxe > bx && !*bxe)
- --n;
- b->wds = n;
- }
- }
- return q;
- }
- /* sulp(x) is a version of ulp(x) that takes bc.scale into account.
- Assuming that x is finite and nonnegative (positive zero is fine
- here) and x / 2^bc.scale is exactly representable as a double,
- sulp(x) is equivalent to 2^bc.scale * ulp(x / 2^bc.scale). */
- static double
- sulp(U *x, BCinfo *bc)
- {
- U u;
- if (bc->scale && 2*P + 1 > (int)((word0(x) & Exp_mask) >> Exp_shift)) {
- /* rv/2^bc->scale is subnormal */
- word0(&u) = (P+2)*Exp_msk1;
- word1(&u) = 0;
- return u.d;
- }
- else {
- assert(word0(x) || word1(x)); /* x != 0.0 */
- return ulp(x);
- }
- }
- /* The bigcomp function handles some hard cases for strtod, for inputs
- with more than STRTOD_DIGLIM digits. It's called once an initial
- estimate for the double corresponding to the input string has
- already been obtained by the code in __Py_dg_strtod.
- The bigcomp function is only called after __Py_dg_strtod has found a
- double value rv such that either rv or rv + 1ulp represents the
- correctly rounded value corresponding to the original string. It
- determines which of these two values is the correct one by
- computing the decimal digits of rv + 0.5ulp and comparing them with
- the corresponding digits of s0.
- In the following, write dv for the absolute value of the number represented
- by the input string.
- Inputs:
- s0 points to the first significant digit of the input string.
- rv is a (possibly scaled) estimate for the closest double value to the
- value represented by the original input to __Py_dg_strtod. If
- bc->scale is nonzero, then rv/2^(bc->scale) is the approximation to
- the input value.
- bc is a struct containing information gathered during the parsing and
- estimation steps of __Py_dg_strtod. Description of fields follows:
- bc->e0 gives the exponent of the input value, such that dv = (integer
- given by the bd->nd digits of s0) * 10**e0
- bc->nd gives the total number of significant digits of s0. It will
- be at least 1.
- bc->nd0 gives the number of significant digits of s0 before the
- decimal separator. If there's no decimal separator, bc->nd0 ==
- bc->nd.
- bc->scale is the value used to scale rv to avoid doing arithmetic with
- subnormal values. It's either 0 or 2*P (=106).
- Outputs:
- On successful exit, rv/2^(bc->scale) is the closest double to dv.
- Returns 0 on success, -1 on failure (e.g., due to a failed malloc call). */
- static int
- bigcomp(U *rv, const char *s0, BCinfo *bc)
- {
- Bigint *b, *d;
- int b2, d2, dd, i, nd, nd0, odd, p2, p5;
- nd = bc->nd;
- nd0 = bc->nd0;
- p5 = nd + bc->e0;
- b = sd2b(rv, bc->scale, &p2);
- if (b == NULL)
- return -1;
- /* record whether the lsb of rv/2^(bc->scale) is odd: in the exact halfway
- case, this is used for round to even. */
- odd = b->x[0] & 1;
- /* left shift b by 1 bit and or a 1 into the least significant bit;
- this gives us b * 2**p2 = rv/2^(bc->scale) + 0.5 ulp. */
- b = lshift(b, 1);
- if (b == NULL)
- return -1;
- b->x[0] |= 1;
- p2--;
- p2 -= p5;
- d = i2b(1);
- if (d == NULL) {
- Bfree(b);
- return -1;
- }
- /* Arrange for convenient computation of quotients:
- * shift left if necessary so divisor has 4 leading 0 bits.
- */
- if (p5 > 0) {
- d = pow5mult(d, p5);
- if (d == NULL) {
- Bfree(b);
- return -1;
- }
- }
- else if (p5 < 0) {
- b = pow5mult(b, -p5);
- if (b == NULL) {
- Bfree(d);
- return -1;
- }
- }
- if (p2 > 0) {
- b2 = p2;
- d2 = 0;
- }
- else {
- b2 = 0;
- d2 = -p2;
- }
- i = dshift(d, d2);
- if ((b2 += i) > 0) {
- b = lshift(b, b2);
- if (b == NULL) {
- Bfree(d);
- return -1;
- }
- }
- if ((d2 += i) > 0) {
- d = lshift(d, d2);
- if (d == NULL) {
- Bfree(b);
- return -1;
- }
- }
- /* Compare s0 with b/d: set dd to -1, 0, or 1 according as s0 < b/d, s0 ==
- * b/d, or s0 > b/d. Here the digits of s0 are thought of as representing
- * a number in the range [0.1, 1). */
- if (cmp(b, d) >= 0)
- /* b/d >= 1 */
- dd = -1;
- else {
- i = 0;
- for(;;) {
- b = multadd(b, 10, 0);
- if (b == NULL) {
- Bfree(d);
- return -1;
- }
- dd = s0[i < nd0 ? i : i+1] - '0' - quorem(b, d);
- i++;
- if (dd)
- break;
- if (!b->x[0] && b->wds == 1) {
- /* b/d == 0 */
- dd = i < nd;
- break;
- }
- if (!(i < nd)) {
- /* b/d != 0, but digits of s0 exhausted */
- dd = -1;
- break;
- }
- }
- }
- Bfree(b);
- Bfree(d);
- if (dd > 0 || (dd == 0 && odd))
- dval(rv) += sulp(rv, bc);
- return 0;
- }
- static double
- __Py_dg_strtod(const char *s00, char **se)
- {
- int bb2, bb5, bbe, bd2, bd5, bs2, c, dsign, e, e1, error;
- int esign, i, j, k, lz, nd, nd0, odd, sign;
- const char *s, *s0, *s1;
- double aadj, aadj1;
- U aadj2, adj, rv, rv0;
- ULong y, z, abs_exp;
- Long L;
- BCinfo bc;
- Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
- dval(&rv) = 0.;
- /* Start parsing. */
- c = *(s = s00);
- /* Parse optional sign, if present. */
- sign = 0;
- switch (c) {
- case '-':
- sign = 1;
- /* no break */
- case '+':
- c = *++s;
- }
- /* Skip leading zeros: lz is true iff there were leading zeros. */
- s1 = s;
- while (c == '0')
- c = *++s;
- lz = s != s1;
- /* Point s0 at the first nonzero digit (if any). nd0 will be the position
- of the point relative to s0. nd will be the total number of digits
- ignoring leading zeros. */
- s0 = s1 = s;
- while ('0' <= c && c <= '9')
- c = *++s;
- nd0 = nd = s - s1;
- /* Parse decimal point and following digits. */
- if (c == '.') {
- c = *++s;
- if (!nd) {
- s1 = s;
- while (c == '0')
- c = *++s;
- lz = lz || s != s1;
- nd0 -= s - s1;
- s0 = s;
- }
- s1 = s;
- while ('0' <= c && c <= '9')
- c = *++s;
- nd += s - s1;
- }
- /* Now lz is true if and only if there were leading zero digits, and nd
- gives the total number of digits ignoring leading zeros. A valid input
- must have at least one digit. */
- if (!nd && !lz) {
- if (se)
- *se = (char *)s00;
- goto parse_error;
- }
- /* Parse exponent. */
- e = 0;
- if (c == 'e' || c == 'E') {
- s00 = s;
- c = *++s;
- /* Exponent sign. */
- esign = 0;
- switch (c) {
- case '-':
- esign = 1;
- /* no break */
- case '+':
- c = *++s;
- }
- /* Skip zeros. lz is true iff there are leading zeros. */
- s1 = s;
- while (c == '0')
- c = *++s;
- lz = s != s1;
- /* Get absolute value of the exponent. */
- s1 = s;
- abs_exp = 0;
- while ('0' <= c && c <= '9') {
- abs_exp = 10*abs_exp + (c - '0');
- c = *++s;
- }
- /* abs_exp will be correct modulo 2**32. But 10**9 < 2**32, so if
- there are at most 9 significant exponent digits then overflow is
- impossible. */
- if (s - s1 > 9 || abs_exp > MAX_ABS_EXP)
- e = (int)MAX_ABS_EXP;
- else
- e = (int)abs_exp;
- if (esign)
- e = -e;
- /* A valid exponent must have at least one digit. */
- if (s == s1 && !lz)
- s = s00;
- }
- /* Adjust exponent to take into account position of the point. */
- e -= nd - nd0;
- if (nd0 <= 0)
- nd0 = nd;
- /* Finished parsing. Set se to indicate how far we parsed */
- if (se)
- *se = (char *)s;
- /* If all digits were zero, exit with return value +-0.0. Otherwise,
- strip trailing zeros: scan back until we hit a nonzero digit. */
- if (!nd)
- goto ret;
- for (i = nd; i > 0; ) {
- --i;
- if (s0[i < nd0 ? i : i+1] != '0') {
- ++i;
- break;
- }
- }
- e += nd - i;
- nd = i;
- if (nd0 > nd)
- nd0 = nd;
- /* Summary of parsing results. After parsing, and dealing with zero
- * inputs, we have values s0, nd0, nd, e, sign, where:
- *
- * - s0 points to the first significant digit of the input string
- *
- * - nd is the total number of significant digits (here, and
- * below, 'significant digits' means the set of digits of the
- * significand of the input that remain after ignoring leading
- * and trailing zeros).
- *
- * - nd0 indicates the position of the decimal point, if present; it
- * satisfies 1 <= nd0 <= nd. The nd significant digits are in
- * s0[0:nd0] and s0[nd0+1:nd+1] using the usual Python half-open slice
- * notation. (If nd0 < nd, then s0[nd0] contains a '.' character; if
- * nd0 == nd, then s0[nd0] could be any non-digit character.)
- *
- * - e is the adjusted exponent: the absolute value of the number
- * represented by the original input string is n * 10**e, where
- * n is the integer represented by the concatenation of
- * s0[0:nd0] and s0[nd0+1:nd+1]
- *
- * - sign gives the sign of the input: 1 for negative, 0 for positive
- *
- * - the first and last significant digits are nonzero
- */
- /* put first DBL_DIG+1 digits into integer y and z.
- *
- * - y contains the value represented by the first min(9, nd)
- * significant digits
- *
- * - if nd > 9, z contains the value represented by significant digits
- * with indices in [9, min(16, nd)). So y * 10**(min(16, nd) - 9) + z
- * gives the value represented by the first min(16, nd) sig. digits.
- */
- bc.e0 = e1 = e;
- y = z = 0;
- for (i = 0; i < nd; i++) {
- if (i < 9)
- y = 10*y + s0[i < nd0 ? i : i+1] - '0';
- else if (i < DBL_DIG+1)
- z = 10*z + s0[i < nd0 ? i : i+1] - '0';
- else
- break;
- }
- k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
- dval(&rv) = y;
- if (k > 9) {
- dval(&rv) = tens[k - 9] * dval(&rv) + z;
- }
- bd0 = 0;
- if (nd <= DBL_DIG
- && Flt_Rounds == 1
- ) {
- if (!e)
- goto ret;
- if (e > 0) {
- if (e <= Ten_pmax) {
- dval(&rv) *= tens[e];
- goto ret;
- }
- i = DBL_DIG - nd;
- if (e <= Ten_pmax + i) {
- /* A fancier test would sometimes let us do
- * this for larger i values.
- */
- e -= i;
- dval(&rv) *= tens[i];
- dval(&rv) *= tens[e];
- goto ret;
- }
- }
- else if (e >= -Ten_pmax) {
- dval(&rv) /= tens[-e];
- goto ret;
- }
- }
- e1 += nd - k;
- bc.scale = 0;
- /* Get starting approximation = rv * 10**e1 */
- if (e1 > 0) {
- if ((i = e1 & 15))
- dval(&rv) *= tens[i];
- if (e1 &= ~15) {
- if (e1 > DBL_MAX_10_EXP)
- goto ovfl;
- e1 >>= 4;
- for(j = 0; e1 > 1; j++, e1 >>= 1)
- if (e1 & 1)
- dval(&rv) *= bigtens[j];
- /* The last multiplication could overflow. */
- word0(&rv) -= P*Exp_msk1;
- dval(&rv) *= bigtens[j];
- if ((z = word0(&rv) & Exp_mask)
- > Exp_msk1*(DBL_MAX_EXP+Bias-P))
- goto ovfl;
- if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
- /* set to largest number */
- /* (Can't trust DBL_MAX) */
- word0(&rv) = Big0;
- word1(&rv) = Big1;
- }
- else
- word0(&rv) += P*Exp_msk1;
- }
- }
- else if (e1 < 0) {
- /* The input decimal value lies in [10**e1, 10**(e1+16)).
- If e1 <= -512, underflow immediately.
- If e1 <= -256, set bc.scale to 2*P.
- So for input value < 1e-256, bc.scale is always set;
- for input value >= 1e-240, bc.scale is never set.
- For input values in [1e-256, 1e-240), bc.scale may or may
- not be set. */
- e1 = -e1;
- if ((i = e1 & 15))
- dval(&rv) /= tens[i];
- if (e1 >>= 4) {
- if (e1 >= 1 << n_bigtens)
- goto undfl;
- if (e1 & Scale_Bit)
- bc.scale = 2*P;
- for(j = 0; e1 > 0; j++, e1 >>= 1)
- if (e1 & 1)
- dval(&rv) *= tinytens[j];
- if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
- >> Exp_shift)) > 0) {
- /* scaled rv is denormal; clear j low bits */
- if (j >= 32) {
- word1(&rv) = 0;
- if (j >= 53)
- word0(&rv) = (P+2)*Exp_msk1;
- else
- word0(&rv) &= 0xffffffff << (j-32);
- }
- else
- word1(&rv) &= 0xffffffff << j;
- }
- if (!dval(&rv))
- goto undfl;
- }
- }
- /* Now the hard part -- adjusting rv to the correct value.*/
- /* Put digits into bd: true value = bd * 10^e */
- bc.nd = nd;
- bc.nd0 = nd0; /* Only needed if nd > STRTOD_DIGLIM, but done here */
- /* to silence an erroneous warning about bc.nd0 */
- /* possibly not being initialized. */
- if (nd > STRTOD_DIGLIM) {
- /* ASSERT(STRTOD_DIGLIM >= 18); 18 == one more than the */
- /* minimum number of decimal digits to distinguish double values */
- /* in IEEE arithmetic. */
- /* Truncate input to 18 significant digits, then discard any trailing
- zeros on the result by updating nd, nd0, e and y suitably. (There's
- no need to update z; it's not reused beyond this point.) */
- for (i = 18; i > 0; ) {
- /* scan back until we hit a nonzero digit. significant digit 'i'
- is s0[i] if i < nd0, s0[i+1] if i >= nd0. */
- --i;
- if (s0[i < nd0 ? i : i+1] != '0') {
- ++i;
- break;
- }
- }
- e += nd - i;
- nd = i;
- if (nd0 > nd)
- nd0 = nd;
- if (nd < 9) { /* must recompute y */
- y = 0;
- for(i = 0; i < nd0; ++i)
- y = 10*y + s0[i] - '0';
- for(; i < nd; ++i)
- y = 10*y + s0[i+1] - '0';
- }
- }
- bd0 = s2b(s0, nd0, nd, y);
- if (bd0 == NULL)
- goto failed_malloc;
- /* Notation for the comments below. Write:
- - dv for the absolute value of the number represented by the original
- decimal input string.
- - if we've truncated dv, write tdv for the truncated value.
- Otherwise, set tdv == dv.
- - srv for the quantity rv/2^bc.scale; so srv is the current binary
- approximation to tdv (and dv). It should be exactly representable
- in an IEEE 754 double.
- */
- for(;;) {
- /* This is the main correction loop for __Py_dg_strtod.
- We've got a decimal value tdv, and a floating-point approximation
- srv=rv/2^bc.scale to tdv. The aim is to determine whether srv is
- close enough (i.e., within 0.5 ulps) to tdv, and to compute a new
- approximation if not.
- To determine whether srv is close enough to tdv, compute integers
- bd, bb and b…
Large files files are truncated, but you can click here to view the full file