ยปCore Development>Code coverage>Python/dtoa.c

Python code coverage for Python/dtoa.c

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