+++ /dev/null
-/****************************************************************\r
- *\r
- * The author of this software is David M. Gay.\r
- *\r
- * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.\r
- *\r
- * Permission to use, copy, modify, and distribute this software for any\r
- * purpose without fee is hereby granted, provided that this entire notice\r
- * is included in all copies of any software which is or includes a copy\r
- * or modification of this software and in all copies of the supporting\r
- * documentation for such software.\r
- *\r
- * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED\r
- * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY\r
- * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY\r
- * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.\r
- *\r
- ***************************************************************/\r
-\r
-/****************************************************************\r
- * This is dtoa.c by David M. Gay, downloaded from\r
- * http://www.netlib.org/fp/dtoa.c on April 15, 2009 and modified for\r
- * inclusion into the Python core by Mark E. T. Dickinson and Eric V. Smith.\r
- *\r
- * Please remember to check http://www.netlib.org/fp regularly (and especially\r
- * before any Python release) for bugfixes and updates.\r
- *\r
- * The major modifications from Gay's original code are as follows:\r
- *\r
- * 0. The original code has been specialized to Python's needs by removing\r
- * many of the #ifdef'd sections. In particular, code to support VAX and\r
- * IBM floating-point formats, hex NaNs, hex floats, locale-aware\r
- * treatment of the decimal point, and setting of the inexact flag have\r
- * been removed.\r
- *\r
- * 1. We use PyMem_Malloc and PyMem_Free in place of malloc and free.\r
- *\r
- * 2. The public functions strtod, dtoa and freedtoa all now have\r
- * a _Py_dg_ prefix.\r
- *\r
- * 3. Instead of assuming that PyMem_Malloc always succeeds, we thread\r
- * PyMem_Malloc failures through the code. The functions\r
- *\r
- * Balloc, multadd, s2b, i2b, mult, pow5mult, lshift, diff, d2b\r
- *\r
- * of return type *Bigint all return NULL to indicate a malloc failure.\r
- * Similarly, rv_alloc and nrv_alloc (return type char *) return NULL on\r
- * failure. bigcomp now has return type int (it used to be void) and\r
- * returns -1 on failure and 0 otherwise. _Py_dg_dtoa returns NULL\r
- * on failure. _Py_dg_strtod indicates failure due to malloc failure\r
- * by returning -1.0, setting errno=ENOMEM and *se to s00.\r
- *\r
- * 4. The static variable dtoa_result has been removed. Callers of\r
- * _Py_dg_dtoa are expected to call _Py_dg_freedtoa to free\r
- * the memory allocated by _Py_dg_dtoa.\r
- *\r
- * 5. The code has been reformatted to better fit with Python's\r
- * C style guide (PEP 7).\r
- *\r
- * 6. A bug in the memory allocation has been fixed: to avoid FREEing memory\r
- * that hasn't been MALLOC'ed, private_mem should only be used when k <=\r
- * Kmax.\r
- *\r
- * 7. _Py_dg_strtod has been modified so that it doesn't accept strings with\r
- * leading whitespace.\r
- *\r
- ***************************************************************/\r
-\r
-/* Please send bug reports for the original dtoa.c code to David M. Gay (dmg\r
- * at acm dot org, with " at " changed at "@" and " dot " changed to ".").\r
- * Please report bugs for this modified version using the Python issue tracker\r
- * (http://bugs.python.org). */\r
-\r
-/* On a machine with IEEE extended-precision registers, it is\r
- * necessary to specify double-precision (53-bit) rounding precision\r
- * before invoking strtod or dtoa. If the machine uses (the equivalent\r
- * of) Intel 80x87 arithmetic, the call\r
- * _control87(PC_53, MCW_PC);\r
- * does this with many compilers. Whether this or another call is\r
- * appropriate depends on the compiler; for this to work, it may be\r
- * necessary to #include "float.h" or another system-dependent header\r
- * file.\r
- */\r
-\r
-/* strtod for IEEE-, VAX-, and IBM-arithmetic machines.\r
- *\r
- * This strtod returns a nearest machine number to the input decimal\r
- * string (or sets errno to ERANGE). With IEEE arithmetic, ties are\r
- * broken by the IEEE round-even rule. Otherwise ties are broken by\r
- * biased rounding (add half and chop).\r
- *\r
- * Inspired loosely by William D. Clinger's paper "How to Read Floating\r
- * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].\r
- *\r
- * Modifications:\r
- *\r
- * 1. We only require IEEE, IBM, or VAX double-precision\r
- * arithmetic (not IEEE double-extended).\r
- * 2. We get by with floating-point arithmetic in a case that\r
- * Clinger missed -- when we're computing d * 10^n\r
- * for a small integer d and the integer n is not too\r
- * much larger than 22 (the maximum integer k for which\r
- * we can represent 10^k exactly), we may be able to\r
- * compute (d*10^k) * 10^(e-k) with just one roundoff.\r
- * 3. Rather than a bit-at-a-time adjustment of the binary\r
- * result in the hard case, we use floating-point\r
- * arithmetic to determine the adjustment to within\r
- * one bit; only in really hard cases do we need to\r
- * compute a second residual.\r
- * 4. Because of 3., we don't need a large table of powers of 10\r
- * for ten-to-e (just some small tables, e.g. of 10^k\r
- * for 0 <= k <= 22).\r
- */\r
-\r
-/* Linking of Python's #defines to Gay's #defines starts here. */\r
-\r
-#include "Python.h"\r
-\r
-/* if PY_NO_SHORT_FLOAT_REPR is defined, then don't even try to compile\r
- the following code */\r
-#ifndef PY_NO_SHORT_FLOAT_REPR\r
-\r
-#include "float.h"\r
-\r
-#define MALLOC PyMem_Malloc\r
-#define FREE PyMem_Free\r
-\r
-/* This code should also work for ARM mixed-endian format on little-endian\r
- machines, where doubles have byte order 45670123 (in increasing address\r
- order, 0 being the least significant byte). */\r
-#ifdef DOUBLE_IS_LITTLE_ENDIAN_IEEE754\r
-# define IEEE_8087\r
-#endif\r
-#if defined(DOUBLE_IS_BIG_ENDIAN_IEEE754) || \\r
- defined(DOUBLE_IS_ARM_MIXED_ENDIAN_IEEE754)\r
-# define IEEE_MC68k\r
-#endif\r
-#if defined(IEEE_8087) + defined(IEEE_MC68k) != 1\r
-#error "Exactly one of IEEE_8087 or IEEE_MC68k should be defined."\r
-#endif\r
-\r
-/* The code below assumes that the endianness of integers matches the\r
- endianness of the two 32-bit words of a double. Check this. */\r
-#if defined(WORDS_BIGENDIAN) && (defined(DOUBLE_IS_LITTLE_ENDIAN_IEEE754) || \\r
- defined(DOUBLE_IS_ARM_MIXED_ENDIAN_IEEE754))\r
-#error "doubles and ints have incompatible endianness"\r
-#endif\r
-\r
-#if !defined(WORDS_BIGENDIAN) && defined(DOUBLE_IS_BIG_ENDIAN_IEEE754)\r
-#error "doubles and ints have incompatible endianness"\r
-#endif\r
-\r
-\r
-#if defined(HAVE_UINT32_T) && defined(HAVE_INT32_T)\r
-typedef PY_UINT32_T ULong;\r
-typedef PY_INT32_T Long;\r
-#else\r
-#error "Failed to find an exact-width 32-bit integer type"\r
-#endif\r
-\r
-#if defined(HAVE_UINT64_T)\r
-#define ULLong PY_UINT64_T\r
-#else\r
-#undef ULLong\r
-#endif\r
-\r
-#undef DEBUG\r
-#ifdef Py_DEBUG\r
-#define DEBUG\r
-#endif\r
-\r
-/* End Python #define linking */\r
-\r
-#ifdef DEBUG\r
-#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}\r
-#endif\r
-\r
-#ifndef PRIVATE_MEM\r
-#define PRIVATE_MEM 2304\r
-#endif\r
-#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))\r
-static double private_mem[PRIVATE_mem], *pmem_next = private_mem;\r
-\r
-#ifdef __cplusplus\r
-extern "C" {\r
-#endif\r
-\r
-typedef union { double d; ULong L[2]; } U;\r
-\r
-#ifdef IEEE_8087\r
-#define word0(x) (x)->L[1]\r
-#define word1(x) (x)->L[0]\r
-#else\r
-#define word0(x) (x)->L[0]\r
-#define word1(x) (x)->L[1]\r
-#endif\r
-#define dval(x) (x)->d\r
-\r
-#ifndef STRTOD_DIGLIM\r
-#define STRTOD_DIGLIM 40\r
-#endif\r
-\r
-/* maximum permitted exponent value for strtod; exponents larger than\r
- MAX_ABS_EXP in absolute value get truncated to +-MAX_ABS_EXP. MAX_ABS_EXP\r
- should fit into an int. */\r
-#ifndef MAX_ABS_EXP\r
-#define MAX_ABS_EXP 1100000000U\r
-#endif\r
-/* Bound on length of pieces of input strings in _Py_dg_strtod; specifically,\r
- this is used to bound the total number of digits ignoring leading zeros and\r
- the number of digits that follow the decimal point. Ideally, MAX_DIGITS\r
- should satisfy MAX_DIGITS + 400 < MAX_ABS_EXP; that ensures that the\r
- exponent clipping in _Py_dg_strtod can't affect the value of the output. */\r
-#ifndef MAX_DIGITS\r
-#define MAX_DIGITS 1000000000U\r
-#endif\r
-\r
-/* Guard against trying to use the above values on unusual platforms with ints\r
- * of width less than 32 bits. */\r
-#if MAX_ABS_EXP > INT_MAX\r
-#error "MAX_ABS_EXP should fit in an int"\r
-#endif\r
-#if MAX_DIGITS > INT_MAX\r
-#error "MAX_DIGITS should fit in an int"\r
-#endif\r
-\r
-/* The following definition of Storeinc is appropriate for MIPS processors.\r
- * An alternative that might be better on some machines is\r
- * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)\r
- */\r
-#if defined(IEEE_8087)\r
-#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \\r
- ((unsigned short *)a)[0] = (unsigned short)c, a++)\r
-#else\r
-#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \\r
- ((unsigned short *)a)[1] = (unsigned short)c, a++)\r
-#endif\r
-\r
-/* #define P DBL_MANT_DIG */\r
-/* Ten_pmax = floor(P*log(2)/log(5)) */\r
-/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */\r
-/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */\r
-/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */\r
-\r
-#define Exp_shift 20\r
-#define Exp_shift1 20\r
-#define Exp_msk1 0x100000\r
-#define Exp_msk11 0x100000\r
-#define Exp_mask 0x7ff00000\r
-#define P 53\r
-#define Nbits 53\r
-#define Bias 1023\r
-#define Emax 1023\r
-#define Emin (-1022)\r
-#define Etiny (-1074) /* smallest denormal is 2**Etiny */\r
-#define Exp_1 0x3ff00000\r
-#define Exp_11 0x3ff00000\r
-#define Ebits 11\r
-#define Frac_mask 0xfffff\r
-#define Frac_mask1 0xfffff\r
-#define Ten_pmax 22\r
-#define Bletch 0x10\r
-#define Bndry_mask 0xfffff\r
-#define Bndry_mask1 0xfffff\r
-#define Sign_bit 0x80000000\r
-#define Log2P 1\r
-#define Tiny0 0\r
-#define Tiny1 1\r
-#define Quick_max 14\r
-#define Int_max 14\r
-\r
-#ifndef Flt_Rounds\r
-#ifdef FLT_ROUNDS\r
-#define Flt_Rounds FLT_ROUNDS\r
-#else\r
-#define Flt_Rounds 1\r
-#endif\r
-#endif /*Flt_Rounds*/\r
-\r
-#define Rounding Flt_Rounds\r
-\r
-#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))\r
-#define Big1 0xffffffff\r
-\r
-/* struct BCinfo is used to pass information from _Py_dg_strtod to bigcomp */\r
-\r
-typedef struct BCinfo BCinfo;\r
-struct\r
-BCinfo {\r
- int e0, nd, nd0, scale;\r
-};\r
-\r
-#define FFFFFFFF 0xffffffffUL\r
-\r
-#define Kmax 7\r
-\r
-/* struct Bigint is used to represent arbitrary-precision integers. These\r
- integers are stored in sign-magnitude format, with the magnitude stored as\r
- an array of base 2**32 digits. Bigints are always normalized: if x is a\r
- Bigint then x->wds >= 1, and either x->wds == 1 or x[wds-1] is nonzero.\r
-\r
- The Bigint fields are as follows:\r
-\r
- - next is a header used by Balloc and Bfree to keep track of lists\r
- of freed Bigints; it's also used for the linked list of\r
- powers of 5 of the form 5**2**i used by pow5mult.\r
- - k indicates which pool this Bigint was allocated from\r
- - maxwds is the maximum number of words space was allocated for\r
- (usually maxwds == 2**k)\r
- - sign is 1 for negative Bigints, 0 for positive. The sign is unused\r
- (ignored on inputs, set to 0 on outputs) in almost all operations\r
- involving Bigints: a notable exception is the diff function, which\r
- ignores signs on inputs but sets the sign of the output correctly.\r
- - wds is the actual number of significant words\r
- - x contains the vector of words (digits) for this Bigint, from least\r
- significant (x[0]) to most significant (x[wds-1]).\r
-*/\r
-\r
-struct\r
-Bigint {\r
- struct Bigint *next;\r
- int k, maxwds, sign, wds;\r
- ULong x[1];\r
-};\r
-\r
-typedef struct Bigint Bigint;\r
-\r
-#ifndef Py_USING_MEMORY_DEBUGGER\r
-\r
-/* Memory management: memory is allocated from, and returned to, Kmax+1 pools\r
- of memory, where pool k (0 <= k <= Kmax) is for Bigints b with b->maxwds ==\r
- 1 << k. These pools are maintained as linked lists, with freelist[k]\r
- pointing to the head of the list for pool k.\r
-\r
- On allocation, if there's no free slot in the appropriate pool, MALLOC is\r
- called to get more memory. This memory is not returned to the system until\r
- Python quits. There's also a private memory pool that's allocated from\r
- in preference to using MALLOC.\r
-\r
- For Bigints with more than (1 << Kmax) digits (which implies at least 1233\r
- decimal digits), memory is directly allocated using MALLOC, and freed using\r
- FREE.\r
-\r
- XXX: it would be easy to bypass this memory-management system and\r
- translate each call to Balloc into a call to PyMem_Malloc, and each\r
- Bfree to PyMem_Free. Investigate whether this has any significant\r
- performance on impact. */\r
-\r
-static Bigint *freelist[Kmax+1];\r
-\r
-/* Allocate space for a Bigint with up to 1<<k digits */\r
-\r
-static Bigint *\r
-Balloc(int k)\r
-{\r
- int x;\r
- Bigint *rv;\r
- unsigned int len;\r
-\r
- if (k <= Kmax && (rv = freelist[k]))\r
- freelist[k] = rv->next;\r
- else {\r
- x = 1 << k;\r
- len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)\r
- /sizeof(double);\r
- if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {\r
- rv = (Bigint*)pmem_next;\r
- pmem_next += len;\r
- }\r
- else {\r
- rv = (Bigint*)MALLOC(len*sizeof(double));\r
- if (rv == NULL)\r
- return NULL;\r
- }\r
- rv->k = k;\r
- rv->maxwds = x;\r
- }\r
- rv->sign = rv->wds = 0;\r
- return rv;\r
-}\r
-\r
-/* Free a Bigint allocated with Balloc */\r
-\r
-static void\r
-Bfree(Bigint *v)\r
-{\r
- if (v) {\r
- if (v->k > Kmax)\r
- FREE((void*)v);\r
- else {\r
- v->next = freelist[v->k];\r
- freelist[v->k] = v;\r
- }\r
- }\r
-}\r
-\r
-#else\r
-\r
-/* Alternative versions of Balloc and Bfree that use PyMem_Malloc and\r
- PyMem_Free directly in place of the custom memory allocation scheme above.\r
- These are provided for the benefit of memory debugging tools like\r
- Valgrind. */\r
-\r
-/* Allocate space for a Bigint with up to 1<<k digits */\r
-\r
-static Bigint *\r
-Balloc(int k)\r
-{\r
- int x;\r
- Bigint *rv;\r
- unsigned int len;\r
-\r
- x = 1 << k;\r
- len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)\r
- /sizeof(double);\r
-\r
- rv = (Bigint*)MALLOC(len*sizeof(double));\r
- if (rv == NULL)\r
- return NULL;\r
-\r
- rv->k = k;\r
- rv->maxwds = x;\r
- rv->sign = rv->wds = 0;\r
- return rv;\r
-}\r
-\r
-/* Free a Bigint allocated with Balloc */\r
-\r
-static void\r
-Bfree(Bigint *v)\r
-{\r
- if (v) {\r
- FREE((void*)v);\r
- }\r
-}\r
-\r
-#endif /* Py_USING_MEMORY_DEBUGGER */\r
-\r
-#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \\r
- y->wds*sizeof(Long) + 2*sizeof(int))\r
-\r
-/* Multiply a Bigint b by m and add a. Either modifies b in place and returns\r
- a pointer to the modified b, or Bfrees b and returns a pointer to a copy.\r
- On failure, return NULL. In this case, b will have been already freed. */\r
-\r
-static Bigint *\r
-multadd(Bigint *b, int m, int a) /* multiply by m and add a */\r
-{\r
- int i, wds;\r
-#ifdef ULLong\r
- ULong *x;\r
- ULLong carry, y;\r
-#else\r
- ULong carry, *x, y;\r
- ULong xi, z;\r
-#endif\r
- Bigint *b1;\r
-\r
- wds = b->wds;\r
- x = b->x;\r
- i = 0;\r
- carry = a;\r
- do {\r
-#ifdef ULLong\r
- y = *x * (ULLong)m + carry;\r
- carry = y >> 32;\r
- *x++ = (ULong)(y & FFFFFFFF);\r
-#else\r
- xi = *x;\r
- y = (xi & 0xffff) * m + carry;\r
- z = (xi >> 16) * m + (y >> 16);\r
- carry = z >> 16;\r
- *x++ = (z << 16) + (y & 0xffff);\r
-#endif\r
- }\r
- while(++i < wds);\r
- if (carry) {\r
- if (wds >= b->maxwds) {\r
- b1 = Balloc(b->k+1);\r
- if (b1 == NULL){\r
- Bfree(b);\r
- return NULL;\r
- }\r
- Bcopy(b1, b);\r
- Bfree(b);\r
- b = b1;\r
- }\r
- b->x[wds++] = (ULong)carry;\r
- b->wds = wds;\r
- }\r
- return b;\r
-}\r
-\r
-/* convert a string s containing nd decimal digits (possibly containing a\r
- decimal separator at position nd0, which is ignored) to a Bigint. This\r
- function carries on where the parsing code in _Py_dg_strtod leaves off: on\r
- entry, y9 contains the result of converting the first 9 digits. Returns\r
- NULL on failure. */\r
-\r
-static Bigint *\r
-s2b(const char *s, int nd0, int nd, ULong y9)\r
-{\r
- Bigint *b;\r
- int i, k;\r
- Long x, y;\r
-\r
- x = (nd + 8) / 9;\r
- for(k = 0, y = 1; x > y; y <<= 1, k++) ;\r
- b = Balloc(k);\r
- if (b == NULL)\r
- return NULL;\r
- b->x[0] = y9;\r
- b->wds = 1;\r
-\r
- if (nd <= 9)\r
- return b;\r
-\r
- s += 9;\r
- for (i = 9; i < nd0; i++) {\r
- b = multadd(b, 10, *s++ - '0');\r
- if (b == NULL)\r
- return NULL;\r
- }\r
- s++;\r
- for(; i < nd; i++) {\r
- b = multadd(b, 10, *s++ - '0');\r
- if (b == NULL)\r
- return NULL;\r
- }\r
- return b;\r
-}\r
-\r
-/* count leading 0 bits in the 32-bit integer x. */\r
-\r
-static int\r
-hi0bits(ULong x)\r
-{\r
- int k = 0;\r
-\r
- if (!(x & 0xffff0000)) {\r
- k = 16;\r
- x <<= 16;\r
- }\r
- if (!(x & 0xff000000)) {\r
- k += 8;\r
- x <<= 8;\r
- }\r
- if (!(x & 0xf0000000)) {\r
- k += 4;\r
- x <<= 4;\r
- }\r
- if (!(x & 0xc0000000)) {\r
- k += 2;\r
- x <<= 2;\r
- }\r
- if (!(x & 0x80000000)) {\r
- k++;\r
- if (!(x & 0x40000000))\r
- return 32;\r
- }\r
- return k;\r
-}\r
-\r
-/* count trailing 0 bits in the 32-bit integer y, and shift y right by that\r
- number of bits. */\r
-\r
-static int\r
-lo0bits(ULong *y)\r
-{\r
- int k;\r
- ULong x = *y;\r
-\r
- if (x & 7) {\r
- if (x & 1)\r
- return 0;\r
- if (x & 2) {\r
- *y = x >> 1;\r
- return 1;\r
- }\r
- *y = x >> 2;\r
- return 2;\r
- }\r
- k = 0;\r
- if (!(x & 0xffff)) {\r
- k = 16;\r
- x >>= 16;\r
- }\r
- if (!(x & 0xff)) {\r
- k += 8;\r
- x >>= 8;\r
- }\r
- if (!(x & 0xf)) {\r
- k += 4;\r
- x >>= 4;\r
- }\r
- if (!(x & 0x3)) {\r
- k += 2;\r
- x >>= 2;\r
- }\r
- if (!(x & 1)) {\r
- k++;\r
- x >>= 1;\r
- if (!x)\r
- return 32;\r
- }\r
- *y = x;\r
- return k;\r
-}\r
-\r
-/* convert a small nonnegative integer to a Bigint */\r
-\r
-static Bigint *\r
-i2b(int i)\r
-{\r
- Bigint *b;\r
-\r
- b = Balloc(1);\r
- if (b == NULL)\r
- return NULL;\r
- b->x[0] = i;\r
- b->wds = 1;\r
- return b;\r
-}\r
-\r
-/* multiply two Bigints. Returns a new Bigint, or NULL on failure. Ignores\r
- the signs of a and b. */\r
-\r
-static Bigint *\r
-mult(Bigint *a, Bigint *b)\r
-{\r
- Bigint *c;\r
- int k, wa, wb, wc;\r
- ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;\r
- ULong y;\r
-#ifdef ULLong\r
- ULLong carry, z;\r
-#else\r
- ULong carry, z;\r
- ULong z2;\r
-#endif\r
-\r
- if ((!a->x[0] && a->wds == 1) || (!b->x[0] && b->wds == 1)) {\r
- c = Balloc(0);\r
- if (c == NULL)\r
- return NULL;\r
- c->wds = 1;\r
- c->x[0] = 0;\r
- return c;\r
- }\r
-\r
- if (a->wds < b->wds) {\r
- c = a;\r
- a = b;\r
- b = c;\r
- }\r
- k = a->k;\r
- wa = a->wds;\r
- wb = b->wds;\r
- wc = wa + wb;\r
- if (wc > a->maxwds)\r
- k++;\r
- c = Balloc(k);\r
- if (c == NULL)\r
- return NULL;\r
- for(x = c->x, xa = x + wc; x < xa; x++)\r
- *x = 0;\r
- xa = a->x;\r
- xae = xa + wa;\r
- xb = b->x;\r
- xbe = xb + wb;\r
- xc0 = c->x;\r
-#ifdef ULLong\r
- for(; xb < xbe; xc0++) {\r
- if ((y = *xb++)) {\r
- x = xa;\r
- xc = xc0;\r
- carry = 0;\r
- do {\r
- z = *x++ * (ULLong)y + *xc + carry;\r
- carry = z >> 32;\r
- *xc++ = (ULong)(z & FFFFFFFF);\r
- }\r
- while(x < xae);\r
- *xc = (ULong)carry;\r
- }\r
- }\r
-#else\r
- for(; xb < xbe; xb++, xc0++) {\r
- if (y = *xb & 0xffff) {\r
- x = xa;\r
- xc = xc0;\r
- carry = 0;\r
- do {\r
- z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;\r
- carry = z >> 16;\r
- z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;\r
- carry = z2 >> 16;\r
- Storeinc(xc, z2, z);\r
- }\r
- while(x < xae);\r
- *xc = carry;\r
- }\r
- if (y = *xb >> 16) {\r
- x = xa;\r
- xc = xc0;\r
- carry = 0;\r
- z2 = *xc;\r
- do {\r
- z = (*x & 0xffff) * y + (*xc >> 16) + carry;\r
- carry = z >> 16;\r
- Storeinc(xc, z, z2);\r
- z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;\r
- carry = z2 >> 16;\r
- }\r
- while(x < xae);\r
- *xc = z2;\r
- }\r
- }\r
-#endif\r
- for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;\r
- c->wds = wc;\r
- return c;\r
-}\r
-\r
-#ifndef Py_USING_MEMORY_DEBUGGER\r
-\r
-/* p5s is a linked list of powers of 5 of the form 5**(2**i), i >= 2 */\r
-\r
-static Bigint *p5s;\r
-\r
-/* multiply the Bigint b by 5**k. Returns a pointer to the result, or NULL on\r
- failure; if the returned pointer is distinct from b then the original\r
- Bigint b will have been Bfree'd. Ignores the sign of b. */\r
-\r
-static Bigint *\r
-pow5mult(Bigint *b, int k)\r
-{\r
- Bigint *b1, *p5, *p51;\r
- int i;\r
- static int p05[3] = { 5, 25, 125 };\r
-\r
- if ((i = k & 3)) {\r
- b = multadd(b, p05[i-1], 0);\r
- if (b == NULL)\r
- return NULL;\r
- }\r
-\r
- if (!(k >>= 2))\r
- return b;\r
- p5 = p5s;\r
- if (!p5) {\r
- /* first time */\r
- p5 = i2b(625);\r
- if (p5 == NULL) {\r
- Bfree(b);\r
- return NULL;\r
- }\r
- p5s = p5;\r
- p5->next = 0;\r
- }\r
- for(;;) {\r
- if (k & 1) {\r
- b1 = mult(b, p5);\r
- Bfree(b);\r
- b = b1;\r
- if (b == NULL)\r
- return NULL;\r
- }\r
- if (!(k >>= 1))\r
- break;\r
- p51 = p5->next;\r
- if (!p51) {\r
- p51 = mult(p5,p5);\r
- if (p51 == NULL) {\r
- Bfree(b);\r
- return NULL;\r
- }\r
- p51->next = 0;\r
- p5->next = p51;\r
- }\r
- p5 = p51;\r
- }\r
- return b;\r
-}\r
-\r
-#else\r
-\r
-/* Version of pow5mult that doesn't cache powers of 5. Provided for\r
- the benefit of memory debugging tools like Valgrind. */\r
-\r
-static Bigint *\r
-pow5mult(Bigint *b, int k)\r
-{\r
- Bigint *b1, *p5, *p51;\r
- int i;\r
- static int p05[3] = { 5, 25, 125 };\r
-\r
- if ((i = k & 3)) {\r
- b = multadd(b, p05[i-1], 0);\r
- if (b == NULL)\r
- return NULL;\r
- }\r
-\r
- if (!(k >>= 2))\r
- return b;\r
- p5 = i2b(625);\r
- if (p5 == NULL) {\r
- Bfree(b);\r
- return NULL;\r
- }\r
-\r
- for(;;) {\r
- if (k & 1) {\r
- b1 = mult(b, p5);\r
- Bfree(b);\r
- b = b1;\r
- if (b == NULL) {\r
- Bfree(p5);\r
- return NULL;\r
- }\r
- }\r
- if (!(k >>= 1))\r
- break;\r
- p51 = mult(p5, p5);\r
- Bfree(p5);\r
- p5 = p51;\r
- if (p5 == NULL) {\r
- Bfree(b);\r
- return NULL;\r
- }\r
- }\r
- Bfree(p5);\r
- return b;\r
-}\r
-\r
-#endif /* Py_USING_MEMORY_DEBUGGER */\r
-\r
-/* shift a Bigint b left by k bits. Return a pointer to the shifted result,\r
- or NULL on failure. If the returned pointer is distinct from b then the\r
- original b will have been Bfree'd. Ignores the sign of b. */\r
-\r
-static Bigint *\r
-lshift(Bigint *b, int k)\r
-{\r
- int i, k1, n, n1;\r
- Bigint *b1;\r
- ULong *x, *x1, *xe, z;\r
-\r
- if (!k || (!b->x[0] && b->wds == 1))\r
- return b;\r
-\r
- n = k >> 5;\r
- k1 = b->k;\r
- n1 = n + b->wds + 1;\r
- for(i = b->maxwds; n1 > i; i <<= 1)\r
- k1++;\r
- b1 = Balloc(k1);\r
- if (b1 == NULL) {\r
- Bfree(b);\r
- return NULL;\r
- }\r
- x1 = b1->x;\r
- for(i = 0; i < n; i++)\r
- *x1++ = 0;\r
- x = b->x;\r
- xe = x + b->wds;\r
- if (k &= 0x1f) {\r
- k1 = 32 - k;\r
- z = 0;\r
- do {\r
- *x1++ = *x << k | z;\r
- z = *x++ >> k1;\r
- }\r
- while(x < xe);\r
- if ((*x1 = z))\r
- ++n1;\r
- }\r
- else do\r
- *x1++ = *x++;\r
- while(x < xe);\r
- b1->wds = n1 - 1;\r
- Bfree(b);\r
- return b1;\r
-}\r
-\r
-/* Do a three-way compare of a and b, returning -1 if a < b, 0 if a == b and\r
- 1 if a > b. Ignores signs of a and b. */\r
-\r
-static int\r
-cmp(Bigint *a, Bigint *b)\r
-{\r
- ULong *xa, *xa0, *xb, *xb0;\r
- int i, j;\r
-\r
- i = a->wds;\r
- j = b->wds;\r
-#ifdef DEBUG\r
- if (i > 1 && !a->x[i-1])\r
- Bug("cmp called with a->x[a->wds-1] == 0");\r
- if (j > 1 && !b->x[j-1])\r
- Bug("cmp called with b->x[b->wds-1] == 0");\r
-#endif\r
- if (i -= j)\r
- return i;\r
- xa0 = a->x;\r
- xa = xa0 + j;\r
- xb0 = b->x;\r
- xb = xb0 + j;\r
- for(;;) {\r
- if (*--xa != *--xb)\r
- return *xa < *xb ? -1 : 1;\r
- if (xa <= xa0)\r
- break;\r
- }\r
- return 0;\r
-}\r
-\r
-/* Take the difference of Bigints a and b, returning a new Bigint. Returns\r
- NULL on failure. The signs of a and b are ignored, but the sign of the\r
- result is set appropriately. */\r
-\r
-static Bigint *\r
-diff(Bigint *a, Bigint *b)\r
-{\r
- Bigint *c;\r
- int i, wa, wb;\r
- ULong *xa, *xae, *xb, *xbe, *xc;\r
-#ifdef ULLong\r
- ULLong borrow, y;\r
-#else\r
- ULong borrow, y;\r
- ULong z;\r
-#endif\r
-\r
- i = cmp(a,b);\r
- if (!i) {\r
- c = Balloc(0);\r
- if (c == NULL)\r
- return NULL;\r
- c->wds = 1;\r
- c->x[0] = 0;\r
- return c;\r
- }\r
- if (i < 0) {\r
- c = a;\r
- a = b;\r
- b = c;\r
- i = 1;\r
- }\r
- else\r
- i = 0;\r
- c = Balloc(a->k);\r
- if (c == NULL)\r
- return NULL;\r
- c->sign = i;\r
- wa = a->wds;\r
- xa = a->x;\r
- xae = xa + wa;\r
- wb = b->wds;\r
- xb = b->x;\r
- xbe = xb + wb;\r
- xc = c->x;\r
- borrow = 0;\r
-#ifdef ULLong\r
- do {\r
- y = (ULLong)*xa++ - *xb++ - borrow;\r
- borrow = y >> 32 & (ULong)1;\r
- *xc++ = (ULong)(y & FFFFFFFF);\r
- }\r
- while(xb < xbe);\r
- while(xa < xae) {\r
- y = *xa++ - borrow;\r
- borrow = y >> 32 & (ULong)1;\r
- *xc++ = (ULong)(y & FFFFFFFF);\r
- }\r
-#else\r
- do {\r
- y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;\r
- borrow = (y & 0x10000) >> 16;\r
- z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;\r
- borrow = (z & 0x10000) >> 16;\r
- Storeinc(xc, z, y);\r
- }\r
- while(xb < xbe);\r
- while(xa < xae) {\r
- y = (*xa & 0xffff) - borrow;\r
- borrow = (y & 0x10000) >> 16;\r
- z = (*xa++ >> 16) - borrow;\r
- borrow = (z & 0x10000) >> 16;\r
- Storeinc(xc, z, y);\r
- }\r
-#endif\r
- while(!*--xc)\r
- wa--;\r
- c->wds = wa;\r
- return c;\r
-}\r
-\r
-/* Given a positive normal double x, return the difference between x and the\r
- next double up. Doesn't give correct results for subnormals. */\r
-\r
-static double\r
-ulp(U *x)\r
-{\r
- Long L;\r
- U u;\r
-\r
- L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;\r
- word0(&u) = L;\r
- word1(&u) = 0;\r
- return dval(&u);\r
-}\r
-\r
-/* Convert a Bigint to a double plus an exponent */\r
-\r
-static double\r
-b2d(Bigint *a, int *e)\r
-{\r
- ULong *xa, *xa0, w, y, z;\r
- int k;\r
- U d;\r
-\r
- xa0 = a->x;\r
- xa = xa0 + a->wds;\r
- y = *--xa;\r
-#ifdef DEBUG\r
- if (!y) Bug("zero y in b2d");\r
-#endif\r
- k = hi0bits(y);\r
- *e = 32 - k;\r
- if (k < Ebits) {\r
- word0(&d) = Exp_1 | y >> (Ebits - k);\r
- w = xa > xa0 ? *--xa : 0;\r
- word1(&d) = y << ((32-Ebits) + k) | w >> (Ebits - k);\r
- goto ret_d;\r
- }\r
- z = xa > xa0 ? *--xa : 0;\r
- if (k -= Ebits) {\r
- word0(&d) = Exp_1 | y << k | z >> (32 - k);\r
- y = xa > xa0 ? *--xa : 0;\r
- word1(&d) = z << k | y >> (32 - k);\r
- }\r
- else {\r
- word0(&d) = Exp_1 | y;\r
- word1(&d) = z;\r
- }\r
- ret_d:\r
- return dval(&d);\r
-}\r
-\r
-/* Convert a scaled double to a Bigint plus an exponent. Similar to d2b,\r
- except that it accepts the scale parameter used in _Py_dg_strtod (which\r
- should be either 0 or 2*P), and the normalization for the return value is\r
- different (see below). On input, d should be finite and nonnegative, and d\r
- / 2**scale should be exactly representable as an IEEE 754 double.\r
-\r
- Returns a Bigint b and an integer e such that\r
-\r
- dval(d) / 2**scale = b * 2**e.\r
-\r
- Unlike d2b, b is not necessarily odd: b and e are normalized so\r
- that either 2**(P-1) <= b < 2**P and e >= Etiny, or b < 2**P\r
- and e == Etiny. This applies equally to an input of 0.0: in that\r
- case the return values are b = 0 and e = Etiny.\r
-\r
- The above normalization ensures that for all possible inputs d,\r
- 2**e gives ulp(d/2**scale).\r
-\r
- Returns NULL on failure.\r
-*/\r
-\r
-static Bigint *\r
-sd2b(U *d, int scale, int *e)\r
-{\r
- Bigint *b;\r
-\r
- b = Balloc(1);\r
- if (b == NULL)\r
- return NULL;\r
- \r
- /* First construct b and e assuming that scale == 0. */\r
- b->wds = 2;\r
- b->x[0] = word1(d);\r
- b->x[1] = word0(d) & Frac_mask;\r
- *e = Etiny - 1 + (int)((word0(d) & Exp_mask) >> Exp_shift);\r
- if (*e < Etiny)\r
- *e = Etiny;\r
- else\r
- b->x[1] |= Exp_msk1;\r
-\r
- /* Now adjust for scale, provided that b != 0. */\r
- if (scale && (b->x[0] || b->x[1])) {\r
- *e -= scale;\r
- if (*e < Etiny) {\r
- scale = Etiny - *e;\r
- *e = Etiny;\r
- /* We can't shift more than P-1 bits without shifting out a 1. */\r
- assert(0 < scale && scale <= P - 1);\r
- if (scale >= 32) {\r
- /* The bits shifted out should all be zero. */\r
- assert(b->x[0] == 0);\r
- b->x[0] = b->x[1];\r
- b->x[1] = 0;\r
- scale -= 32;\r
- }\r
- if (scale) {\r
- /* The bits shifted out should all be zero. */\r
- assert(b->x[0] << (32 - scale) == 0);\r
- b->x[0] = (b->x[0] >> scale) | (b->x[1] << (32 - scale));\r
- b->x[1] >>= scale;\r
- }\r
- }\r
- }\r
- /* Ensure b is normalized. */\r
- if (!b->x[1])\r
- b->wds = 1;\r
-\r
- return b;\r
-}\r
-\r
-/* Convert a double to a Bigint plus an exponent. Return NULL on failure.\r
-\r
- Given a finite nonzero double d, return an odd Bigint b and exponent *e\r
- such that fabs(d) = b * 2**e. On return, *bbits gives the number of\r
- significant bits of b; that is, 2**(*bbits-1) <= b < 2**(*bbits).\r
-\r
- If d is zero, then b == 0, *e == -1010, *bbits = 0.\r
- */\r
-\r
-static Bigint *\r
-d2b(U *d, int *e, int *bits)\r
-{\r
- Bigint *b;\r
- int de, k;\r
- ULong *x, y, z;\r
- int i;\r
-\r
- b = Balloc(1);\r
- if (b == NULL)\r
- return NULL;\r
- x = b->x;\r
-\r
- z = word0(d) & Frac_mask;\r
- word0(d) &= 0x7fffffff; /* clear sign bit, which we ignore */\r
- if ((de = (int)(word0(d) >> Exp_shift)))\r
- z |= Exp_msk1;\r
- if ((y = word1(d))) {\r
- if ((k = lo0bits(&y))) {\r
- x[0] = y | z << (32 - k);\r
- z >>= k;\r
- }\r
- else\r
- x[0] = y;\r
- i =\r
- b->wds = (x[1] = z) ? 2 : 1;\r
- }\r
- else {\r
- k = lo0bits(&z);\r
- x[0] = z;\r
- i =\r
- b->wds = 1;\r
- k += 32;\r
- }\r
- if (de) {\r
- *e = de - Bias - (P-1) + k;\r
- *bits = P - k;\r
- }\r
- else {\r
- *e = de - Bias - (P-1) + 1 + k;\r
- *bits = 32*i - hi0bits(x[i-1]);\r
- }\r
- return b;\r
-}\r
-\r
-/* Compute the ratio of two Bigints, as a double. The result may have an\r
- error of up to 2.5 ulps. */\r
-\r
-static double\r
-ratio(Bigint *a, Bigint *b)\r
-{\r
- U da, db;\r
- int k, ka, kb;\r
-\r
- dval(&da) = b2d(a, &ka);\r
- dval(&db) = b2d(b, &kb);\r
- k = ka - kb + 32*(a->wds - b->wds);\r
- if (k > 0)\r
- word0(&da) += k*Exp_msk1;\r
- else {\r
- k = -k;\r
- word0(&db) += k*Exp_msk1;\r
- }\r
- return dval(&da) / dval(&db);\r
-}\r
-\r
-static const double\r
-tens[] = {\r
- 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,\r
- 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,\r
- 1e20, 1e21, 1e22\r
-};\r
-\r
-static const double\r
-bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };\r
-static const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,\r
- 9007199254740992.*9007199254740992.e-256\r
- /* = 2^106 * 1e-256 */\r
-};\r
-/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */\r
-/* flag unnecessarily. It leads to a song and dance at the end of strtod. */\r
-#define Scale_Bit 0x10\r
-#define n_bigtens 5\r
-\r
-#define ULbits 32\r
-#define kshift 5\r
-#define kmask 31\r
-\r
-\r
-static int\r
-dshift(Bigint *b, int p2)\r
-{\r
- int rv = hi0bits(b->x[b->wds-1]) - 4;\r
- if (p2 > 0)\r
- rv -= p2;\r
- return rv & kmask;\r
-}\r
-\r
-/* special case of Bigint division. The quotient is always in the range 0 <=\r
- quotient < 10, and on entry the divisor S is normalized so that its top 4\r
- bits (28--31) are zero and bit 27 is set. */\r
-\r
-static int\r
-quorem(Bigint *b, Bigint *S)\r
-{\r
- int n;\r
- ULong *bx, *bxe, q, *sx, *sxe;\r
-#ifdef ULLong\r
- ULLong borrow, carry, y, ys;\r
-#else\r
- ULong borrow, carry, y, ys;\r
- ULong si, z, zs;\r
-#endif\r
-\r
- n = S->wds;\r
-#ifdef DEBUG\r
- /*debug*/ if (b->wds > n)\r
- /*debug*/ Bug("oversize b in quorem");\r
-#endif\r
- if (b->wds < n)\r
- return 0;\r
- sx = S->x;\r
- sxe = sx + --n;\r
- bx = b->x;\r
- bxe = bx + n;\r
- q = *bxe / (*sxe + 1); /* ensure q <= true quotient */\r
-#ifdef DEBUG\r
- /*debug*/ if (q > 9)\r
- /*debug*/ Bug("oversized quotient in quorem");\r
-#endif\r
- if (q) {\r
- borrow = 0;\r
- carry = 0;\r
- do {\r
-#ifdef ULLong\r
- ys = *sx++ * (ULLong)q + carry;\r
- carry = ys >> 32;\r
- y = *bx - (ys & FFFFFFFF) - borrow;\r
- borrow = y >> 32 & (ULong)1;\r
- *bx++ = (ULong)(y & FFFFFFFF);\r
-#else\r
- si = *sx++;\r
- ys = (si & 0xffff) * q + carry;\r
- zs = (si >> 16) * q + (ys >> 16);\r
- carry = zs >> 16;\r
- y = (*bx & 0xffff) - (ys & 0xffff) - borrow;\r
- borrow = (y & 0x10000) >> 16;\r
- z = (*bx >> 16) - (zs & 0xffff) - borrow;\r
- borrow = (z & 0x10000) >> 16;\r
- Storeinc(bx, z, y);\r
-#endif\r
- }\r
- while(sx <= sxe);\r
- if (!*bxe) {\r
- bx = b->x;\r
- while(--bxe > bx && !*bxe)\r
- --n;\r
- b->wds = n;\r
- }\r
- }\r
- if (cmp(b, S) >= 0) {\r
- q++;\r
- borrow = 0;\r
- carry = 0;\r
- bx = b->x;\r
- sx = S->x;\r
- do {\r
-#ifdef ULLong\r
- ys = *sx++ + carry;\r
- carry = ys >> 32;\r
- y = *bx - (ys & FFFFFFFF) - borrow;\r
- borrow = y >> 32 & (ULong)1;\r
- *bx++ = (ULong)(y & FFFFFFFF);\r
-#else\r
- si = *sx++;\r
- ys = (si & 0xffff) + carry;\r
- zs = (si >> 16) + (ys >> 16);\r
- carry = zs >> 16;\r
- y = (*bx & 0xffff) - (ys & 0xffff) - borrow;\r
- borrow = (y & 0x10000) >> 16;\r
- z = (*bx >> 16) - (zs & 0xffff) - borrow;\r
- borrow = (z & 0x10000) >> 16;\r
- Storeinc(bx, z, y);\r
-#endif\r
- }\r
- while(sx <= sxe);\r
- bx = b->x;\r
- bxe = bx + n;\r
- if (!*bxe) {\r
- while(--bxe > bx && !*bxe)\r
- --n;\r
- b->wds = n;\r
- }\r
- }\r
- return q;\r
-}\r
-\r
-/* sulp(x) is a version of ulp(x) that takes bc.scale into account.\r
-\r
- Assuming that x is finite and nonnegative (positive zero is fine\r
- here) and x / 2^bc.scale is exactly representable as a double,\r
- sulp(x) is equivalent to 2^bc.scale * ulp(x / 2^bc.scale). */\r
-\r
-static double\r
-sulp(U *x, BCinfo *bc)\r
-{\r
- U u;\r
-\r
- if (bc->scale && 2*P + 1 > (int)((word0(x) & Exp_mask) >> Exp_shift)) {\r
- /* rv/2^bc->scale is subnormal */\r
- word0(&u) = (P+2)*Exp_msk1;\r
- word1(&u) = 0;\r
- return u.d;\r
- }\r
- else {\r
- assert(word0(x) || word1(x)); /* x != 0.0 */\r
- return ulp(x);\r
- }\r
-}\r
-\r
-/* The bigcomp function handles some hard cases for strtod, for inputs\r
- with more than STRTOD_DIGLIM digits. It's called once an initial\r
- estimate for the double corresponding to the input string has\r
- already been obtained by the code in _Py_dg_strtod.\r
-\r
- The bigcomp function is only called after _Py_dg_strtod has found a\r
- double value rv such that either rv or rv + 1ulp represents the\r
- correctly rounded value corresponding to the original string. It\r
- determines which of these two values is the correct one by\r
- computing the decimal digits of rv + 0.5ulp and comparing them with\r
- the corresponding digits of s0.\r
-\r
- In the following, write dv for the absolute value of the number represented\r
- by the input string.\r
-\r
- Inputs:\r
-\r
- s0 points to the first significant digit of the input string.\r
-\r
- rv is a (possibly scaled) estimate for the closest double value to the\r
- value represented by the original input to _Py_dg_strtod. If\r
- bc->scale is nonzero, then rv/2^(bc->scale) is the approximation to\r
- the input value.\r
-\r
- bc is a struct containing information gathered during the parsing and\r
- estimation steps of _Py_dg_strtod. Description of fields follows:\r
-\r
- bc->e0 gives the exponent of the input value, such that dv = (integer\r
- given by the bd->nd digits of s0) * 10**e0\r
-\r
- bc->nd gives the total number of significant digits of s0. It will\r
- be at least 1.\r
-\r
- bc->nd0 gives the number of significant digits of s0 before the\r
- decimal separator. If there's no decimal separator, bc->nd0 ==\r
- bc->nd.\r
-\r
- bc->scale is the value used to scale rv to avoid doing arithmetic with\r
- subnormal values. It's either 0 or 2*P (=106).\r
-\r
- Outputs:\r
-\r
- On successful exit, rv/2^(bc->scale) is the closest double to dv.\r
-\r
- Returns 0 on success, -1 on failure (e.g., due to a failed malloc call). */\r
-\r
-static int\r
-bigcomp(U *rv, const char *s0, BCinfo *bc)\r
-{\r
- Bigint *b, *d;\r
- int b2, d2, dd, i, nd, nd0, odd, p2, p5;\r
-\r
- nd = bc->nd;\r
- nd0 = bc->nd0;\r
- p5 = nd + bc->e0;\r
- b = sd2b(rv, bc->scale, &p2);\r
- if (b == NULL)\r
- return -1;\r
-\r
- /* record whether the lsb of rv/2^(bc->scale) is odd: in the exact halfway\r
- case, this is used for round to even. */\r
- odd = b->x[0] & 1;\r
-\r
- /* left shift b by 1 bit and or a 1 into the least significant bit;\r
- this gives us b * 2**p2 = rv/2^(bc->scale) + 0.5 ulp. */\r
- b = lshift(b, 1);\r
- if (b == NULL)\r
- return -1;\r
- b->x[0] |= 1;\r
- p2--;\r
-\r
- p2 -= p5;\r
- d = i2b(1);\r
- if (d == NULL) {\r
- Bfree(b);\r
- return -1;\r
- }\r
- /* Arrange for convenient computation of quotients:\r
- * shift left if necessary so divisor has 4 leading 0 bits.\r
- */\r
- if (p5 > 0) {\r
- d = pow5mult(d, p5);\r
- if (d == NULL) {\r
- Bfree(b);\r
- return -1;\r
- }\r
- }\r
- else if (p5 < 0) {\r
- b = pow5mult(b, -p5);\r
- if (b == NULL) {\r
- Bfree(d);\r
- return -1;\r
- }\r
- }\r
- if (p2 > 0) {\r
- b2 = p2;\r
- d2 = 0;\r
- }\r
- else {\r
- b2 = 0;\r
- d2 = -p2;\r
- }\r
- i = dshift(d, d2);\r
- if ((b2 += i) > 0) {\r
- b = lshift(b, b2);\r
- if (b == NULL) {\r
- Bfree(d);\r
- return -1;\r
- }\r
- }\r
- if ((d2 += i) > 0) {\r
- d = lshift(d, d2);\r
- if (d == NULL) {\r
- Bfree(b);\r
- return -1;\r
- }\r
- }\r
-\r
- /* Compare s0 with b/d: set dd to -1, 0, or 1 according as s0 < b/d, s0 ==\r
- * b/d, or s0 > b/d. Here the digits of s0 are thought of as representing\r
- * a number in the range [0.1, 1). */\r
- if (cmp(b, d) >= 0)\r
- /* b/d >= 1 */\r
- dd = -1;\r
- else {\r
- i = 0;\r
- for(;;) {\r
- b = multadd(b, 10, 0);\r
- if (b == NULL) {\r
- Bfree(d);\r
- return -1;\r
- }\r
- dd = s0[i < nd0 ? i : i+1] - '0' - quorem(b, d);\r
- i++;\r
-\r
- if (dd)\r
- break;\r
- if (!b->x[0] && b->wds == 1) {\r
- /* b/d == 0 */\r
- dd = i < nd;\r
- break;\r
- }\r
- if (!(i < nd)) {\r
- /* b/d != 0, but digits of s0 exhausted */\r
- dd = -1;\r
- break;\r
- }\r
- }\r
- }\r
- Bfree(b);\r
- Bfree(d);\r
- if (dd > 0 || (dd == 0 && odd))\r
- dval(rv) += sulp(rv, bc);\r
- return 0;\r
-}\r
-\r
-double\r
-_Py_dg_strtod(const char *s00, char **se)\r
-{\r
- int bb2, bb5, bbe, bd2, bd5, bs2, c, dsign, e, e1, error;\r
- int esign, i, j, k, lz, nd, nd0, odd, sign;\r
- const char *s, *s0, *s1;\r
- double aadj, aadj1;\r
- U aadj2, adj, rv, rv0;\r
- ULong y, z, abs_exp;\r
- Long L;\r
- BCinfo bc;\r
- Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;\r
- size_t ndigits, fraclen;\r
-\r
- dval(&rv) = 0.;\r
-\r
- /* Start parsing. */\r
- c = *(s = s00);\r
-\r
- /* Parse optional sign, if present. */\r
- sign = 0;\r
- switch (c) {\r
- case '-':\r
- sign = 1;\r
- /* no break */\r
- case '+':\r
- c = *++s;\r
- }\r
-\r
- /* Skip leading zeros: lz is true iff there were leading zeros. */\r
- s1 = s;\r
- while (c == '0')\r
- c = *++s;\r
- lz = s != s1;\r
-\r
- /* Point s0 at the first nonzero digit (if any). fraclen will be the\r
- number of digits between the decimal point and the end of the\r
- digit string. ndigits will be the total number of digits ignoring\r
- leading zeros. */\r
- s0 = s1 = s;\r
- while ('0' <= c && c <= '9')\r
- c = *++s;\r
- ndigits = s - s1;\r
- fraclen = 0;\r
-\r
- /* Parse decimal point and following digits. */\r
- if (c == '.') {\r
- c = *++s;\r
- if (!ndigits) {\r
- s1 = s;\r
- while (c == '0')\r
- c = *++s;\r
- lz = lz || s != s1;\r
- fraclen += (s - s1);\r
- s0 = s;\r
- }\r
- s1 = s;\r
- while ('0' <= c && c <= '9')\r
- c = *++s;\r
- ndigits += s - s1;\r
- fraclen += s - s1;\r
- }\r
-\r
- /* Now lz is true if and only if there were leading zero digits, and\r
- ndigits gives the total number of digits ignoring leading zeros. A\r
- valid input must have at least one digit. */\r
- if (!ndigits && !lz) {\r
- if (se)\r
- *se = (char *)s00;\r
- goto parse_error;\r
- }\r
-\r
- /* Range check ndigits and fraclen to make sure that they, and values\r
- computed with them, can safely fit in an int. */\r
- if (ndigits > MAX_DIGITS || fraclen > MAX_DIGITS) {\r
- if (se)\r
- *se = (char *)s00;\r
- goto parse_error;\r
- }\r
- nd = (int)ndigits;\r
- nd0 = (int)ndigits - (int)fraclen;\r
-\r
- /* Parse exponent. */\r
- e = 0;\r
- if (c == 'e' || c == 'E') {\r
- s00 = s;\r
- c = *++s;\r
-\r
- /* Exponent sign. */\r
- esign = 0;\r
- switch (c) {\r
- case '-':\r
- esign = 1;\r
- /* no break */\r
- case '+':\r
- c = *++s;\r
- }\r
-\r
- /* Skip zeros. lz is true iff there are leading zeros. */\r
- s1 = s;\r
- while (c == '0')\r
- c = *++s;\r
- lz = s != s1;\r
-\r
- /* Get absolute value of the exponent. */\r
- s1 = s;\r
- abs_exp = 0;\r
- while ('0' <= c && c <= '9') {\r
- abs_exp = 10*abs_exp + (c - '0');\r
- c = *++s;\r
- }\r
-\r
- /* abs_exp will be correct modulo 2**32. But 10**9 < 2**32, so if\r
- there are at most 9 significant exponent digits then overflow is\r
- impossible. */\r
- if (s - s1 > 9 || abs_exp > MAX_ABS_EXP)\r
- e = (int)MAX_ABS_EXP;\r
- else\r
- e = (int)abs_exp;\r
- if (esign)\r
- e = -e;\r
-\r
- /* A valid exponent must have at least one digit. */\r
- if (s == s1 && !lz)\r
- s = s00;\r
- }\r
-\r
- /* Adjust exponent to take into account position of the point. */\r
- e -= nd - nd0;\r
- if (nd0 <= 0)\r
- nd0 = nd;\r
-\r
- /* Finished parsing. Set se to indicate how far we parsed */\r
- if (se)\r
- *se = (char *)s;\r
-\r
- /* If all digits were zero, exit with return value +-0.0. Otherwise,\r
- strip trailing zeros: scan back until we hit a nonzero digit. */\r
- if (!nd)\r
- goto ret;\r
- for (i = nd; i > 0; ) {\r
- --i;\r
- if (s0[i < nd0 ? i : i+1] != '0') {\r
- ++i;\r
- break;\r
- }\r
- }\r
- e += nd - i;\r
- nd = i;\r
- if (nd0 > nd)\r
- nd0 = nd;\r
-\r
- /* Summary of parsing results. After parsing, and dealing with zero\r
- * inputs, we have values s0, nd0, nd, e, sign, where:\r
- *\r
- * - s0 points to the first significant digit of the input string\r
- *\r
- * - nd is the total number of significant digits (here, and\r
- * below, 'significant digits' means the set of digits of the\r
- * significand of the input that remain after ignoring leading\r
- * and trailing zeros).\r
- *\r
- * - nd0 indicates the position of the decimal point, if present; it\r
- * satisfies 1 <= nd0 <= nd. The nd significant digits are in\r
- * s0[0:nd0] and s0[nd0+1:nd+1] using the usual Python half-open slice\r
- * notation. (If nd0 < nd, then s0[nd0] contains a '.' character; if\r
- * nd0 == nd, then s0[nd0] could be any non-digit character.)\r
- *\r
- * - e is the adjusted exponent: the absolute value of the number\r
- * represented by the original input string is n * 10**e, where\r
- * n is the integer represented by the concatenation of\r
- * s0[0:nd0] and s0[nd0+1:nd+1]\r
- *\r
- * - sign gives the sign of the input: 1 for negative, 0 for positive\r
- *\r
- * - the first and last significant digits are nonzero\r
- */\r
-\r
- /* put first DBL_DIG+1 digits into integer y and z.\r
- *\r
- * - y contains the value represented by the first min(9, nd)\r
- * significant digits\r
- *\r
- * - if nd > 9, z contains the value represented by significant digits\r
- * with indices in [9, min(16, nd)). So y * 10**(min(16, nd) - 9) + z\r
- * gives the value represented by the first min(16, nd) sig. digits.\r
- */\r
-\r
- bc.e0 = e1 = e;\r
- y = z = 0;\r
- for (i = 0; i < nd; i++) {\r
- if (i < 9)\r
- y = 10*y + s0[i < nd0 ? i : i+1] - '0';\r
- else if (i < DBL_DIG+1)\r
- z = 10*z + s0[i < nd0 ? i : i+1] - '0';\r
- else\r
- break;\r
- }\r
-\r
- k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;\r
- dval(&rv) = y;\r
- if (k > 9) {\r
- dval(&rv) = tens[k - 9] * dval(&rv) + z;\r
- }\r
- bd0 = 0;\r
- if (nd <= DBL_DIG\r
- && Flt_Rounds == 1\r
- ) {\r
- if (!e)\r
- goto ret;\r
- if (e > 0) {\r
- if (e <= Ten_pmax) {\r
- dval(&rv) *= tens[e];\r
- goto ret;\r
- }\r
- i = DBL_DIG - nd;\r
- if (e <= Ten_pmax + i) {\r
- /* A fancier test would sometimes let us do\r
- * this for larger i values.\r
- */\r
- e -= i;\r
- dval(&rv) *= tens[i];\r
- dval(&rv) *= tens[e];\r
- goto ret;\r
- }\r
- }\r
- else if (e >= -Ten_pmax) {\r
- dval(&rv) /= tens[-e];\r
- goto ret;\r
- }\r
- }\r
- e1 += nd - k;\r
-\r
- bc.scale = 0;\r
-\r
- /* Get starting approximation = rv * 10**e1 */\r
-\r
- if (e1 > 0) {\r
- if ((i = e1 & 15))\r
- dval(&rv) *= tens[i];\r
- if (e1 &= ~15) {\r
- if (e1 > DBL_MAX_10_EXP)\r
- goto ovfl;\r
- e1 >>= 4;\r
- for(j = 0; e1 > 1; j++, e1 >>= 1)\r
- if (e1 & 1)\r
- dval(&rv) *= bigtens[j];\r
- /* The last multiplication could overflow. */\r
- word0(&rv) -= P*Exp_msk1;\r
- dval(&rv) *= bigtens[j];\r
- if ((z = word0(&rv) & Exp_mask)\r
- > Exp_msk1*(DBL_MAX_EXP+Bias-P))\r
- goto ovfl;\r
- if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {\r
- /* set to largest number */\r
- /* (Can't trust DBL_MAX) */\r
- word0(&rv) = Big0;\r
- word1(&rv) = Big1;\r
- }\r
- else\r
- word0(&rv) += P*Exp_msk1;\r
- }\r
- }\r
- else if (e1 < 0) {\r
- /* The input decimal value lies in [10**e1, 10**(e1+16)).\r
-\r
- If e1 <= -512, underflow immediately.\r
- If e1 <= -256, set bc.scale to 2*P.\r
-\r
- So for input value < 1e-256, bc.scale is always set;\r
- for input value >= 1e-240, bc.scale is never set.\r
- For input values in [1e-256, 1e-240), bc.scale may or may\r
- not be set. */\r
-\r
- e1 = -e1;\r
- if ((i = e1 & 15))\r
- dval(&rv) /= tens[i];\r
- if (e1 >>= 4) {\r
- if (e1 >= 1 << n_bigtens)\r
- goto undfl;\r
- if (e1 & Scale_Bit)\r
- bc.scale = 2*P;\r
- for(j = 0; e1 > 0; j++, e1 >>= 1)\r
- if (e1 & 1)\r
- dval(&rv) *= tinytens[j];\r
- if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)\r
- >> Exp_shift)) > 0) {\r
- /* scaled rv is denormal; clear j low bits */\r
- if (j >= 32) {\r
- word1(&rv) = 0;\r
- if (j >= 53)\r
- word0(&rv) = (P+2)*Exp_msk1;\r
- else\r
- word0(&rv) &= 0xffffffff << (j-32);\r
- }\r
- else\r
- word1(&rv) &= 0xffffffff << j;\r
- }\r
- if (!dval(&rv))\r
- goto undfl;\r
- }\r
- }\r
-\r
- /* Now the hard part -- adjusting rv to the correct value.*/\r
-\r
- /* Put digits into bd: true value = bd * 10^e */\r
-\r
- bc.nd = nd;\r
- bc.nd0 = nd0; /* Only needed if nd > STRTOD_DIGLIM, but done here */\r
- /* to silence an erroneous warning about bc.nd0 */\r
- /* possibly not being initialized. */\r
- if (nd > STRTOD_DIGLIM) {\r
- /* ASSERT(STRTOD_DIGLIM >= 18); 18 == one more than the */\r
- /* minimum number of decimal digits to distinguish double values */\r
- /* in IEEE arithmetic. */\r
-\r
- /* Truncate input to 18 significant digits, then discard any trailing\r
- zeros on the result by updating nd, nd0, e and y suitably. (There's\r
- no need to update z; it's not reused beyond this point.) */\r
- for (i = 18; i > 0; ) {\r
- /* scan back until we hit a nonzero digit. significant digit 'i'\r
- is s0[i] if i < nd0, s0[i+1] if i >= nd0. */\r
- --i;\r
- if (s0[i < nd0 ? i : i+1] != '0') {\r
- ++i;\r
- break;\r
- }\r
- }\r
- e += nd - i;\r
- nd = i;\r
- if (nd0 > nd)\r
- nd0 = nd;\r
- if (nd < 9) { /* must recompute y */\r
- y = 0;\r
- for(i = 0; i < nd0; ++i)\r
- y = 10*y + s0[i] - '0';\r
- for(; i < nd; ++i)\r
- y = 10*y + s0[i+1] - '0';\r
- }\r
- }\r
- bd0 = s2b(s0, nd0, nd, y);\r
- if (bd0 == NULL)\r
- goto failed_malloc;\r
-\r
- /* Notation for the comments below. Write:\r
-\r
- - dv for the absolute value of the number represented by the original\r
- decimal input string.\r
-\r
- - if we've truncated dv, write tdv for the truncated value.\r
- Otherwise, set tdv == dv.\r
-\r
- - srv for the quantity rv/2^bc.scale; so srv is the current binary\r
- approximation to tdv (and dv). It should be exactly representable\r
- in an IEEE 754 double.\r
- */\r
-\r
- for(;;) {\r
-\r
- /* This is the main correction loop for _Py_dg_strtod.\r
-\r
- We've got a decimal value tdv, and a floating-point approximation\r
- srv=rv/2^bc.scale to tdv. The aim is to determine whether srv is\r
- close enough (i.e., within 0.5 ulps) to tdv, and to compute a new\r
- approximation if not.\r
-\r
- To determine whether srv is close enough to tdv, compute integers\r
- bd, bb and bs proportional to tdv, srv and 0.5 ulp(srv)\r
- respectively, and then use integer arithmetic to determine whether\r
- |tdv - srv| is less than, equal to, or greater than 0.5 ulp(srv).\r
- */\r
-\r
- bd = Balloc(bd0->k);\r
- if (bd == NULL) {\r
- Bfree(bd0);\r
- goto failed_malloc;\r
- }\r
- Bcopy(bd, bd0);\r
- bb = sd2b(&rv, bc.scale, &bbe); /* srv = bb * 2^bbe */\r
- if (bb == NULL) {\r
- Bfree(bd);\r
- Bfree(bd0);\r
- goto failed_malloc;\r
- }\r
- /* Record whether lsb of bb is odd, in case we need this\r
- for the round-to-even step later. */\r
- odd = bb->x[0] & 1;\r
-\r
- /* tdv = bd * 10**e; srv = bb * 2**bbe */\r
- bs = i2b(1);\r
- if (bs == NULL) {\r
- Bfree(bb);\r
- Bfree(bd);\r
- Bfree(bd0);\r
- goto failed_malloc;\r
- }\r
-\r
- if (e >= 0) {\r
- bb2 = bb5 = 0;\r
- bd2 = bd5 = e;\r
- }\r
- else {\r
- bb2 = bb5 = -e;\r
- bd2 = bd5 = 0;\r
- }\r
- if (bbe >= 0)\r
- bb2 += bbe;\r
- else\r
- bd2 -= bbe;\r
- bs2 = bb2;\r
- bb2++;\r
- bd2++;\r
-\r
- /* At this stage bd5 - bb5 == e == bd2 - bb2 + bbe, bb2 - bs2 == 1,\r
- and bs == 1, so:\r
-\r
- tdv == bd * 10**e = bd * 2**(bbe - bb2 + bd2) * 5**(bd5 - bb5)\r
- srv == bb * 2**bbe = bb * 2**(bbe - bb2 + bb2)\r
- 0.5 ulp(srv) == 2**(bbe-1) = bs * 2**(bbe - bb2 + bs2)\r
-\r
- It follows that:\r
-\r
- M * tdv = bd * 2**bd2 * 5**bd5\r
- M * srv = bb * 2**bb2 * 5**bb5\r
- M * 0.5 ulp(srv) = bs * 2**bs2 * 5**bb5\r
-\r
- for some constant M. (Actually, M == 2**(bb2 - bbe) * 5**bb5, but\r
- this fact is not needed below.)\r
- */\r
-\r
- /* Remove factor of 2**i, where i = min(bb2, bd2, bs2). */\r
- i = bb2 < bd2 ? bb2 : bd2;\r
- if (i > bs2)\r
- i = bs2;\r
- if (i > 0) {\r
- bb2 -= i;\r
- bd2 -= i;\r
- bs2 -= i;\r
- }\r
-\r
- /* Scale bb, bd, bs by the appropriate powers of 2 and 5. */\r
- if (bb5 > 0) {\r
- bs = pow5mult(bs, bb5);\r
- if (bs == NULL) {\r
- Bfree(bb);\r
- Bfree(bd);\r
- Bfree(bd0);\r
- goto failed_malloc;\r
- }\r
- bb1 = mult(bs, bb);\r
- Bfree(bb);\r
- bb = bb1;\r
- if (bb == NULL) {\r
- Bfree(bs);\r
- Bfree(bd);\r
- Bfree(bd0);\r
- goto failed_malloc;\r
- }\r
- }\r
- if (bb2 > 0) {\r
- bb = lshift(bb, bb2);\r
- if (bb == NULL) {\r
- Bfree(bs);\r
- Bfree(bd);\r
- Bfree(bd0);\r
- goto failed_malloc;\r
- }\r
- }\r
- if (bd5 > 0) {\r
- bd = pow5mult(bd, bd5);\r
- if (bd == NULL) {\r
- Bfree(bb);\r
- Bfree(bs);\r
- Bfree(bd0);\r
- goto failed_malloc;\r
- }\r
- }\r
- if (bd2 > 0) {\r
- bd = lshift(bd, bd2);\r
- if (bd == NULL) {\r
- Bfree(bb);\r
- Bfree(bs);\r
- Bfree(bd0);\r
- goto failed_malloc;\r
- }\r
- }\r
- if (bs2 > 0) {\r
- bs = lshift(bs, bs2);\r
- if (bs == NULL) {\r
- Bfree(bb);\r
- Bfree(bd);\r
- Bfree(bd0);\r
- goto failed_malloc;\r
- }\r
- }\r
-\r
- /* Now bd, bb and bs are scaled versions of tdv, srv and 0.5 ulp(srv),\r
- respectively. Compute the difference |tdv - srv|, and compare\r
- with 0.5 ulp(srv). */\r
-\r
- delta = diff(bb, bd);\r
- if (delta == NULL) {\r
- Bfree(bb);\r
- Bfree(bs);\r
- Bfree(bd);\r
- Bfree(bd0);\r
- goto failed_malloc;\r
- }\r
- dsign = delta->sign;\r
- delta->sign = 0;\r
- i = cmp(delta, bs);\r
- if (bc.nd > nd && i <= 0) {\r
- if (dsign)\r
- break; /* Must use bigcomp(). */\r
-\r
- /* Here rv overestimates the truncated decimal value by at most\r
- 0.5 ulp(rv). Hence rv either overestimates the true decimal\r
- value by <= 0.5 ulp(rv), or underestimates it by some small\r
- amount (< 0.1 ulp(rv)); either way, rv is within 0.5 ulps of\r
- the true decimal value, so it's possible to exit.\r
-\r
- Exception: if scaled rv is a normal exact power of 2, but not\r
- DBL_MIN, then rv - 0.5 ulp(rv) takes us all the way down to the\r
- next double, so the correctly rounded result is either rv - 0.5\r
- ulp(rv) or rv; in this case, use bigcomp to distinguish. */\r
-\r
- if (!word1(&rv) && !(word0(&rv) & Bndry_mask)) {\r
- /* rv can't be 0, since it's an overestimate for some\r
- nonzero value. So rv is a normal power of 2. */\r
- j = (int)(word0(&rv) & Exp_mask) >> Exp_shift;\r
- /* rv / 2^bc.scale = 2^(j - 1023 - bc.scale); use bigcomp if\r
- rv / 2^bc.scale >= 2^-1021. */\r
- if (j - bc.scale >= 2) {\r
- dval(&rv) -= 0.5 * sulp(&rv, &bc);\r
- break; /* Use bigcomp. */\r
- }\r
- }\r
-\r
- {\r
- bc.nd = nd;\r
- i = -1; /* Discarded digits make delta smaller. */\r
- }\r
- }\r
-\r
- if (i < 0) {\r
- /* Error is less than half an ulp -- check for\r
- * special case of mantissa a power of two.\r
- */\r
- if (dsign || word1(&rv) || word0(&rv) & Bndry_mask\r
- || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1\r
- ) {\r
- break;\r
- }\r
- if (!delta->x[0] && delta->wds <= 1) {\r
- /* exact result */\r
- break;\r
- }\r
- delta = lshift(delta,Log2P);\r
- if (delta == NULL) {\r
- Bfree(bb);\r
- Bfree(bs);\r
- Bfree(bd);\r
- Bfree(bd0);\r
- goto failed_malloc;\r
- }\r
- if (cmp(delta, bs) > 0)\r
- goto drop_down;\r
- break;\r
- }\r
- if (i == 0) {\r
- /* exactly half-way between */\r
- if (dsign) {\r
- if ((word0(&rv) & Bndry_mask1) == Bndry_mask1\r
- && word1(&rv) == (\r
- (bc.scale &&\r
- (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) ?\r
- (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :\r
- 0xffffffff)) {\r
- /*boundary case -- increment exponent*/\r
- word0(&rv) = (word0(&rv) & Exp_mask)\r
- + Exp_msk1\r
- ;\r
- word1(&rv) = 0;\r
- dsign = 0;\r
- break;\r
- }\r
- }\r
- else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {\r
- drop_down:\r
- /* boundary case -- decrement exponent */\r
- if (bc.scale) {\r
- L = word0(&rv) & Exp_mask;\r
- if (L <= (2*P+1)*Exp_msk1) {\r
- if (L > (P+2)*Exp_msk1)\r
- /* round even ==> */\r
- /* accept rv */\r
- break;\r
- /* rv = smallest denormal */\r
- if (bc.nd > nd)\r
- break;\r
- goto undfl;\r
- }\r
- }\r
- L = (word0(&rv) & Exp_mask) - Exp_msk1;\r
- word0(&rv) = L | Bndry_mask1;\r
- word1(&rv) = 0xffffffff;\r
- break;\r
- }\r
- if (!odd)\r
- break;\r
- if (dsign)\r
- dval(&rv) += sulp(&rv, &bc);\r
- else {\r
- dval(&rv) -= sulp(&rv, &bc);\r
- if (!dval(&rv)) {\r
- if (bc.nd >nd)\r
- break;\r
- goto undfl;\r
- }\r
- }\r
- dsign = 1 - dsign;\r
- break;\r
- }\r
- if ((aadj = ratio(delta, bs)) <= 2.) {\r
- if (dsign)\r
- aadj = aadj1 = 1.;\r
- else if (word1(&rv) || word0(&rv) & Bndry_mask) {\r
- if (word1(&rv) == Tiny1 && !word0(&rv)) {\r
- if (bc.nd >nd)\r
- break;\r
- goto undfl;\r
- }\r
- aadj = 1.;\r
- aadj1 = -1.;\r
- }\r
- else {\r
- /* special case -- power of FLT_RADIX to be */\r
- /* rounded down... */\r
-\r
- if (aadj < 2./FLT_RADIX)\r
- aadj = 1./FLT_RADIX;\r
- else\r
- aadj *= 0.5;\r
- aadj1 = -aadj;\r
- }\r
- }\r
- else {\r
- aadj *= 0.5;\r
- aadj1 = dsign ? aadj : -aadj;\r
- if (Flt_Rounds == 0)\r
- aadj1 += 0.5;\r
- }\r
- y = word0(&rv) & Exp_mask;\r
-\r
- /* Check for overflow */\r
-\r
- if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {\r
- dval(&rv0) = dval(&rv);\r
- word0(&rv) -= P*Exp_msk1;\r
- adj.d = aadj1 * ulp(&rv);\r
- dval(&rv) += adj.d;\r
- if ((word0(&rv) & Exp_mask) >=\r
- Exp_msk1*(DBL_MAX_EXP+Bias-P)) {\r
- if (word0(&rv0) == Big0 && word1(&rv0) == Big1) {\r
- Bfree(bb);\r
- Bfree(bd);\r
- Bfree(bs);\r
- Bfree(bd0);\r
- Bfree(delta);\r
- goto ovfl;\r
- }\r
- word0(&rv) = Big0;\r
- word1(&rv) = Big1;\r
- goto cont;\r
- }\r
- else\r
- word0(&rv) += P*Exp_msk1;\r
- }\r
- else {\r
- if (bc.scale && y <= 2*P*Exp_msk1) {\r
- if (aadj <= 0x7fffffff) {\r
- if ((z = (ULong)aadj) <= 0)\r
- z = 1;\r
- aadj = z;\r
- aadj1 = dsign ? aadj : -aadj;\r
- }\r
- dval(&aadj2) = aadj1;\r
- word0(&aadj2) += (2*P+1)*Exp_msk1 - y;\r
- aadj1 = dval(&aadj2);\r
- }\r
- adj.d = aadj1 * ulp(&rv);\r
- dval(&rv) += adj.d;\r
- }\r
- z = word0(&rv) & Exp_mask;\r
- if (bc.nd == nd) {\r
- if (!bc.scale)\r
- if (y == z) {\r
- /* Can we stop now? */\r
- L = (Long)aadj;\r
- aadj -= L;\r
- /* The tolerances below are conservative. */\r
- if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) {\r
- if (aadj < .4999999 || aadj > .5000001)\r
- break;\r
- }\r
- else if (aadj < .4999999/FLT_RADIX)\r
- break;\r
- }\r
- }\r
- cont:\r
- Bfree(bb);\r
- Bfree(bd);\r
- Bfree(bs);\r
- Bfree(delta);\r
- }\r
- Bfree(bb);\r
- Bfree(bd);\r
- Bfree(bs);\r
- Bfree(bd0);\r
- Bfree(delta);\r
- if (bc.nd > nd) {\r
- error = bigcomp(&rv, s0, &bc);\r
- if (error)\r
- goto failed_malloc;\r
- }\r
-\r
- if (bc.scale) {\r
- word0(&rv0) = Exp_1 - 2*P*Exp_msk1;\r
- word1(&rv0) = 0;\r
- dval(&rv) *= dval(&rv0);\r
- }\r
-\r
- ret:\r
- return sign ? -dval(&rv) : dval(&rv);\r
-\r
- parse_error:\r
- return 0.0;\r
-\r
- failed_malloc:\r
- errno = ENOMEM;\r
- return -1.0;\r
-\r
- undfl:\r
- return sign ? -0.0 : 0.0;\r
-\r
- ovfl:\r
- errno = ERANGE;\r
- /* Can't trust HUGE_VAL */\r
- word0(&rv) = Exp_mask;\r
- word1(&rv) = 0;\r
- return sign ? -dval(&rv) : dval(&rv);\r
-\r
-}\r
-\r
-static char *\r
-rv_alloc(int i)\r
-{\r
- int j, k, *r;\r
-\r
- j = sizeof(ULong);\r
- for(k = 0;\r
- sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (unsigned)i;\r
- j <<= 1)\r
- k++;\r
- r = (int*)Balloc(k);\r
- if (r == NULL)\r
- return NULL;\r
- *r = k;\r
- return (char *)(r+1);\r
-}\r
-\r
-static char *\r
-nrv_alloc(char *s, char **rve, int n)\r
-{\r
- char *rv, *t;\r
-\r
- rv = rv_alloc(n);\r
- if (rv == NULL)\r
- return NULL;\r
- t = rv;\r
- while((*t = *s++)) t++;\r
- if (rve)\r
- *rve = t;\r
- return rv;\r
-}\r
-\r
-/* freedtoa(s) must be used to free values s returned by dtoa\r
- * when MULTIPLE_THREADS is #defined. It should be used in all cases,\r
- * but for consistency with earlier versions of dtoa, it is optional\r
- * when MULTIPLE_THREADS is not defined.\r
- */\r
-\r
-void\r
-_Py_dg_freedtoa(char *s)\r
-{\r
- Bigint *b = (Bigint *)((int *)s - 1);\r
- b->maxwds = 1 << (b->k = *(int*)b);\r
- Bfree(b);\r
-}\r
-\r
-/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.\r
- *\r
- * Inspired by "How to Print Floating-Point Numbers Accurately" by\r
- * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].\r
- *\r
- * Modifications:\r
- * 1. Rather than iterating, we use a simple numeric overestimate\r
- * to determine k = floor(log10(d)). We scale relevant\r
- * quantities using O(log2(k)) rather than O(k) multiplications.\r
- * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't\r
- * try to generate digits strictly left to right. Instead, we\r
- * compute with fewer bits and propagate the carry if necessary\r
- * when rounding the final digit up. This is often faster.\r
- * 3. Under the assumption that input will be rounded nearest,\r
- * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.\r
- * That is, we allow equality in stopping tests when the\r
- * round-nearest rule will give the same floating-point value\r
- * as would satisfaction of the stopping test with strict\r
- * inequality.\r
- * 4. We remove common factors of powers of 2 from relevant\r
- * quantities.\r
- * 5. When converting floating-point integers less than 1e16,\r
- * we use floating-point arithmetic rather than resorting\r
- * to multiple-precision integers.\r
- * 6. When asked to produce fewer than 15 digits, we first try\r
- * to get by with floating-point arithmetic; we resort to\r
- * multiple-precision integer arithmetic only if we cannot\r
- * guarantee that the floating-point calculation has given\r
- * the correctly rounded result. For k requested digits and\r
- * "uniformly" distributed input, the probability is\r
- * something like 10^(k-15) that we must resort to the Long\r
- * calculation.\r
- */\r
-\r
-/* Additional notes (METD): (1) returns NULL on failure. (2) to avoid memory\r
- leakage, a successful call to _Py_dg_dtoa should always be matched by a\r
- call to _Py_dg_freedtoa. */\r
-\r
-char *\r
-_Py_dg_dtoa(double dd, int mode, int ndigits,\r
- int *decpt, int *sign, char **rve)\r
-{\r
- /* Arguments ndigits, decpt, sign are similar to those\r
- of ecvt and fcvt; trailing zeros are suppressed from\r
- the returned string. If not null, *rve is set to point\r
- to the end of the return value. If d is +-Infinity or NaN,\r
- then *decpt is set to 9999.\r
-\r
- mode:\r
- 0 ==> shortest string that yields d when read in\r
- and rounded to nearest.\r
- 1 ==> like 0, but with Steele & White stopping rule;\r
- e.g. with IEEE P754 arithmetic , mode 0 gives\r
- 1e23 whereas mode 1 gives 9.999999999999999e22.\r
- 2 ==> max(1,ndigits) significant digits. This gives a\r
- return value similar to that of ecvt, except\r
- that trailing zeros are suppressed.\r
- 3 ==> through ndigits past the decimal point. This\r
- gives a return value similar to that from fcvt,\r
- except that trailing zeros are suppressed, and\r
- ndigits can be negative.\r
- 4,5 ==> similar to 2 and 3, respectively, but (in\r
- round-nearest mode) with the tests of mode 0 to\r
- possibly return a shorter string that rounds to d.\r
- With IEEE arithmetic and compilation with\r
- -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same\r
- as modes 2 and 3 when FLT_ROUNDS != 1.\r
- 6-9 ==> Debugging modes similar to mode - 4: don't try\r
- fast floating-point estimate (if applicable).\r
-\r
- Values of mode other than 0-9 are treated as mode 0.\r
-\r
- Sufficient space is allocated to the return value\r
- to hold the suppressed trailing zeros.\r
- */\r
-\r
- int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,\r
- j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,\r
- spec_case, try_quick;\r
- Long L;\r
- int denorm;\r
- ULong x;\r
- Bigint *b, *b1, *delta, *mlo, *mhi, *S;\r
- U d2, eps, u;\r
- double ds;\r
- char *s, *s0;\r
-\r
- /* set pointers to NULL, to silence gcc compiler warnings and make\r
- cleanup easier on error */\r
- mlo = mhi = S = 0;\r
- s0 = 0;\r
-\r
- u.d = dd;\r
- if (word0(&u) & Sign_bit) {\r
- /* set sign for everything, including 0's and NaNs */\r
- *sign = 1;\r
- word0(&u) &= ~Sign_bit; /* clear sign bit */\r
- }\r
- else\r
- *sign = 0;\r
-\r
- /* quick return for Infinities, NaNs and zeros */\r
- if ((word0(&u) & Exp_mask) == Exp_mask)\r
- {\r
- /* Infinity or NaN */\r
- *decpt = 9999;\r
- if (!word1(&u) && !(word0(&u) & 0xfffff))\r
- return nrv_alloc("Infinity", rve, 8);\r
- return nrv_alloc("NaN", rve, 3);\r
- }\r
- if (!dval(&u)) {\r
- *decpt = 1;\r
- return nrv_alloc("0", rve, 1);\r
- }\r
-\r
- /* compute k = floor(log10(d)). The computation may leave k\r
- one too large, but should never leave k too small. */\r
- b = d2b(&u, &be, &bbits);\r
- if (b == NULL)\r
- goto failed_malloc;\r
- if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {\r
- dval(&d2) = dval(&u);\r
- word0(&d2) &= Frac_mask1;\r
- word0(&d2) |= Exp_11;\r
-\r
- /* log(x) ~=~ log(1.5) + (x-1.5)/1.5\r
- * log10(x) = log(x) / log(10)\r
- * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))\r
- * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)\r
- *\r
- * This suggests computing an approximation k to log10(d) by\r
- *\r
- * k = (i - Bias)*0.301029995663981\r
- * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );\r
- *\r
- * We want k to be too large rather than too small.\r
- * The error in the first-order Taylor series approximation\r
- * is in our favor, so we just round up the constant enough\r
- * to compensate for any error in the multiplication of\r
- * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,\r
- * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,\r
- * adding 1e-13 to the constant term more than suffices.\r
- * Hence we adjust the constant term to 0.1760912590558.\r
- * (We could get a more accurate k by invoking log10,\r
- * but this is probably not worthwhile.)\r
- */\r
-\r
- i -= Bias;\r
- denorm = 0;\r
- }\r
- else {\r
- /* d is denormalized */\r
-\r
- i = bbits + be + (Bias + (P-1) - 1);\r
- x = i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32)\r
- : word1(&u) << (32 - i);\r
- dval(&d2) = x;\r
- word0(&d2) -= 31*Exp_msk1; /* adjust exponent */\r
- i -= (Bias + (P-1) - 1) + 1;\r
- denorm = 1;\r
- }\r
- ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 +\r
- i*0.301029995663981;\r
- k = (int)ds;\r
- if (ds < 0. && ds != k)\r
- k--; /* want k = floor(ds) */\r
- k_check = 1;\r
- if (k >= 0 && k <= Ten_pmax) {\r
- if (dval(&u) < tens[k])\r
- k--;\r
- k_check = 0;\r
- }\r
- j = bbits - i - 1;\r
- if (j >= 0) {\r
- b2 = 0;\r
- s2 = j;\r
- }\r
- else {\r
- b2 = -j;\r
- s2 = 0;\r
- }\r
- if (k >= 0) {\r
- b5 = 0;\r
- s5 = k;\r
- s2 += k;\r
- }\r
- else {\r
- b2 -= k;\r
- b5 = -k;\r
- s5 = 0;\r
- }\r
- if (mode < 0 || mode > 9)\r
- mode = 0;\r
-\r
- try_quick = 1;\r
-\r
- if (mode > 5) {\r
- mode -= 4;\r
- try_quick = 0;\r
- }\r
- leftright = 1;\r
- ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */\r
- /* silence erroneous "gcc -Wall" warning. */\r
- switch(mode) {\r
- case 0:\r
- case 1:\r
- i = 18;\r
- ndigits = 0;\r
- break;\r
- case 2:\r
- leftright = 0;\r
- /* no break */\r
- case 4:\r
- if (ndigits <= 0)\r
- ndigits = 1;\r
- ilim = ilim1 = i = ndigits;\r
- break;\r
- case 3:\r
- leftright = 0;\r
- /* no break */\r
- case 5:\r
- i = ndigits + k + 1;\r
- ilim = i;\r
- ilim1 = i - 1;\r
- if (i <= 0)\r
- i = 1;\r
- }\r
- s0 = rv_alloc(i);\r
- if (s0 == NULL)\r
- goto failed_malloc;\r
- s = s0;\r
-\r
-\r
- if (ilim >= 0 && ilim <= Quick_max && try_quick) {\r
-\r
- /* Try to get by with floating-point arithmetic. */\r
-\r
- i = 0;\r
- dval(&d2) = dval(&u);\r
- k0 = k;\r
- ilim0 = ilim;\r
- ieps = 2; /* conservative */\r
- if (k > 0) {\r
- ds = tens[k&0xf];\r
- j = k >> 4;\r
- if (j & Bletch) {\r
- /* prevent overflows */\r
- j &= Bletch - 1;\r
- dval(&u) /= bigtens[n_bigtens-1];\r
- ieps++;\r
- }\r
- for(; j; j >>= 1, i++)\r
- if (j & 1) {\r
- ieps++;\r
- ds *= bigtens[i];\r
- }\r
- dval(&u) /= ds;\r
- }\r
- else if ((j1 = -k)) {\r
- dval(&u) *= tens[j1 & 0xf];\r
- for(j = j1 >> 4; j; j >>= 1, i++)\r
- if (j & 1) {\r
- ieps++;\r
- dval(&u) *= bigtens[i];\r
- }\r
- }\r
- if (k_check && dval(&u) < 1. && ilim > 0) {\r
- if (ilim1 <= 0)\r
- goto fast_failed;\r
- ilim = ilim1;\r
- k--;\r
- dval(&u) *= 10.;\r
- ieps++;\r
- }\r
- dval(&eps) = ieps*dval(&u) + 7.;\r
- word0(&eps) -= (P-1)*Exp_msk1;\r
- if (ilim == 0) {\r
- S = mhi = 0;\r
- dval(&u) -= 5.;\r
- if (dval(&u) > dval(&eps))\r
- goto one_digit;\r
- if (dval(&u) < -dval(&eps))\r
- goto no_digits;\r
- goto fast_failed;\r
- }\r
- if (leftright) {\r
- /* Use Steele & White method of only\r
- * generating digits needed.\r
- */\r
- dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);\r
- for(i = 0;;) {\r
- L = (Long)dval(&u);\r
- dval(&u) -= L;\r
- *s++ = '0' + (int)L;\r
- if (dval(&u) < dval(&eps))\r
- goto ret1;\r
- if (1. - dval(&u) < dval(&eps))\r
- goto bump_up;\r
- if (++i >= ilim)\r
- break;\r
- dval(&eps) *= 10.;\r
- dval(&u) *= 10.;\r
- }\r
- }\r
- else {\r
- /* Generate ilim digits, then fix them up. */\r
- dval(&eps) *= tens[ilim-1];\r
- for(i = 1;; i++, dval(&u) *= 10.) {\r
- L = (Long)(dval(&u));\r
- if (!(dval(&u) -= L))\r
- ilim = i;\r
- *s++ = '0' + (int)L;\r
- if (i == ilim) {\r
- if (dval(&u) > 0.5 + dval(&eps))\r
- goto bump_up;\r
- else if (dval(&u) < 0.5 - dval(&eps)) {\r
- while(*--s == '0');\r
- s++;\r
- goto ret1;\r
- }\r
- break;\r
- }\r
- }\r
- }\r
- fast_failed:\r
- s = s0;\r
- dval(&u) = dval(&d2);\r
- k = k0;\r
- ilim = ilim0;\r
- }\r
-\r
- /* Do we have a "small" integer? */\r
-\r
- if (be >= 0 && k <= Int_max) {\r
- /* Yes. */\r
- ds = tens[k];\r
- if (ndigits < 0 && ilim <= 0) {\r
- S = mhi = 0;\r
- if (ilim < 0 || dval(&u) <= 5*ds)\r
- goto no_digits;\r
- goto one_digit;\r
- }\r
- for(i = 1;; i++, dval(&u) *= 10.) {\r
- L = (Long)(dval(&u) / ds);\r
- dval(&u) -= L*ds;\r
- *s++ = '0' + (int)L;\r
- if (!dval(&u)) {\r
- break;\r
- }\r
- if (i == ilim) {\r
- dval(&u) += dval(&u);\r
- if (dval(&u) > ds || (dval(&u) == ds && L & 1)) {\r
- bump_up:\r
- while(*--s == '9')\r
- if (s == s0) {\r
- k++;\r
- *s = '0';\r
- break;\r
- }\r
- ++*s++;\r
- }\r
- break;\r
- }\r
- }\r
- goto ret1;\r
- }\r
-\r
- m2 = b2;\r
- m5 = b5;\r
- if (leftright) {\r
- i =\r
- denorm ? be + (Bias + (P-1) - 1 + 1) :\r
- 1 + P - bbits;\r
- b2 += i;\r
- s2 += i;\r
- mhi = i2b(1);\r
- if (mhi == NULL)\r
- goto failed_malloc;\r
- }\r
- if (m2 > 0 && s2 > 0) {\r
- i = m2 < s2 ? m2 : s2;\r
- b2 -= i;\r
- m2 -= i;\r
- s2 -= i;\r
- }\r
- if (b5 > 0) {\r
- if (leftright) {\r
- if (m5 > 0) {\r
- mhi = pow5mult(mhi, m5);\r
- if (mhi == NULL)\r
- goto failed_malloc;\r
- b1 = mult(mhi, b);\r
- Bfree(b);\r
- b = b1;\r
- if (b == NULL)\r
- goto failed_malloc;\r
- }\r
- if ((j = b5 - m5)) {\r
- b = pow5mult(b, j);\r
- if (b == NULL)\r
- goto failed_malloc;\r
- }\r
- }\r
- else {\r
- b = pow5mult(b, b5);\r
- if (b == NULL)\r
- goto failed_malloc;\r
- }\r
- }\r
- S = i2b(1);\r
- if (S == NULL)\r
- goto failed_malloc;\r
- if (s5 > 0) {\r
- S = pow5mult(S, s5);\r
- if (S == NULL)\r
- goto failed_malloc;\r
- }\r
-\r
- /* Check for special case that d is a normalized power of 2. */\r
-\r
- spec_case = 0;\r
- if ((mode < 2 || leftright)\r
- ) {\r
- if (!word1(&u) && !(word0(&u) & Bndry_mask)\r
- && word0(&u) & (Exp_mask & ~Exp_msk1)\r
- ) {\r
- /* The special case */\r
- b2 += Log2P;\r
- s2 += Log2P;\r
- spec_case = 1;\r
- }\r
- }\r
-\r
- /* Arrange for convenient computation of quotients:\r
- * shift left if necessary so divisor has 4 leading 0 bits.\r
- *\r
- * Perhaps we should just compute leading 28 bits of S once\r
- * and for all and pass them and a shift to quorem, so it\r
- * can do shifts and ors to compute the numerator for q.\r
- */\r
-#define iInc 28\r
- i = dshift(S, s2);\r
- b2 += i;\r
- m2 += i;\r
- s2 += i;\r
- if (b2 > 0) {\r
- b = lshift(b, b2);\r
- if (b == NULL)\r
- goto failed_malloc;\r
- }\r
- if (s2 > 0) {\r
- S = lshift(S, s2);\r
- if (S == NULL)\r
- goto failed_malloc;\r
- }\r
- if (k_check) {\r
- if (cmp(b,S) < 0) {\r
- k--;\r
- b = multadd(b, 10, 0); /* we botched the k estimate */\r
- if (b == NULL)\r
- goto failed_malloc;\r
- if (leftright) {\r
- mhi = multadd(mhi, 10, 0);\r
- if (mhi == NULL)\r
- goto failed_malloc;\r
- }\r
- ilim = ilim1;\r
- }\r
- }\r
- if (ilim <= 0 && (mode == 3 || mode == 5)) {\r
- if (ilim < 0) {\r
- /* no digits, fcvt style */\r
- no_digits:\r
- k = -1 - ndigits;\r
- goto ret;\r
- }\r
- else {\r
- S = multadd(S, 5, 0);\r
- if (S == NULL)\r
- goto failed_malloc;\r
- if (cmp(b, S) <= 0)\r
- goto no_digits;\r
- }\r
- one_digit:\r
- *s++ = '1';\r
- k++;\r
- goto ret;\r
- }\r
- if (leftright) {\r
- if (m2 > 0) {\r
- mhi = lshift(mhi, m2);\r
- if (mhi == NULL)\r
- goto failed_malloc;\r
- }\r
-\r
- /* Compute mlo -- check for special case\r
- * that d is a normalized power of 2.\r
- */\r
-\r
- mlo = mhi;\r
- if (spec_case) {\r
- mhi = Balloc(mhi->k);\r
- if (mhi == NULL)\r
- goto failed_malloc;\r
- Bcopy(mhi, mlo);\r
- mhi = lshift(mhi, Log2P);\r
- if (mhi == NULL)\r
- goto failed_malloc;\r
- }\r
-\r
- for(i = 1;;i++) {\r
- dig = quorem(b,S) + '0';\r
- /* Do we yet have the shortest decimal string\r
- * that will round to d?\r
- */\r
- j = cmp(b, mlo);\r
- delta = diff(S, mhi);\r
- if (delta == NULL)\r
- goto failed_malloc;\r
- j1 = delta->sign ? 1 : cmp(b, delta);\r
- Bfree(delta);\r
- if (j1 == 0 && mode != 1 && !(word1(&u) & 1)\r
- ) {\r
- if (dig == '9')\r
- goto round_9_up;\r
- if (j > 0)\r
- dig++;\r
- *s++ = dig;\r
- goto ret;\r
- }\r
- if (j < 0 || (j == 0 && mode != 1\r
- && !(word1(&u) & 1)\r
- )) {\r
- if (!b->x[0] && b->wds <= 1) {\r
- goto accept_dig;\r
- }\r
- if (j1 > 0) {\r
- b = lshift(b, 1);\r
- if (b == NULL)\r
- goto failed_malloc;\r
- j1 = cmp(b, S);\r
- if ((j1 > 0 || (j1 == 0 && dig & 1))\r
- && dig++ == '9')\r
- goto round_9_up;\r
- }\r
- accept_dig:\r
- *s++ = dig;\r
- goto ret;\r
- }\r
- if (j1 > 0) {\r
- if (dig == '9') { /* possible if i == 1 */\r
- round_9_up:\r
- *s++ = '9';\r
- goto roundoff;\r
- }\r
- *s++ = dig + 1;\r
- goto ret;\r
- }\r
- *s++ = dig;\r
- if (i == ilim)\r
- break;\r
- b = multadd(b, 10, 0);\r
- if (b == NULL)\r
- goto failed_malloc;\r
- if (mlo == mhi) {\r
- mlo = mhi = multadd(mhi, 10, 0);\r
- if (mlo == NULL)\r
- goto failed_malloc;\r
- }\r
- else {\r
- mlo = multadd(mlo, 10, 0);\r
- if (mlo == NULL)\r
- goto failed_malloc;\r
- mhi = multadd(mhi, 10, 0);\r
- if (mhi == NULL)\r
- goto failed_malloc;\r
- }\r
- }\r
- }\r
- else\r
- for(i = 1;; i++) {\r
- *s++ = dig = quorem(b,S) + '0';\r
- if (!b->x[0] && b->wds <= 1) {\r
- goto ret;\r
- }\r
- if (i >= ilim)\r
- break;\r
- b = multadd(b, 10, 0);\r
- if (b == NULL)\r
- goto failed_malloc;\r
- }\r
-\r
- /* Round off last digit */\r
-\r
- b = lshift(b, 1);\r
- if (b == NULL)\r
- goto failed_malloc;\r
- j = cmp(b, S);\r
- if (j > 0 || (j == 0 && dig & 1)) {\r
- roundoff:\r
- while(*--s == '9')\r
- if (s == s0) {\r
- k++;\r
- *s++ = '1';\r
- goto ret;\r
- }\r
- ++*s++;\r
- }\r
- else {\r
- while(*--s == '0');\r
- s++;\r
- }\r
- ret:\r
- Bfree(S);\r
- if (mhi) {\r
- if (mlo && mlo != mhi)\r
- Bfree(mlo);\r
- Bfree(mhi);\r
- }\r
- ret1:\r
- Bfree(b);\r
- *s = 0;\r
- *decpt = k + 1;\r
- if (rve)\r
- *rve = s;\r
- return s0;\r
- failed_malloc:\r
- if (S)\r
- Bfree(S);\r
- if (mlo && mlo != mhi)\r
- Bfree(mlo);\r
- if (mhi)\r
- Bfree(mhi);\r
- if (b)\r
- Bfree(b);\r
- if (s0)\r
- _Py_dg_freedtoa(s0);\r
- return NULL;\r
-}\r
-#ifdef __cplusplus\r
-}\r
-#endif\r
-\r
-#endif /* PY_NO_SHORT_FLOAT_REPR */\r
projects / mirror_edk2.git / blobdiff