1 /****************************************************************
3 * The author of this software is David M. Gay.
5 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
7 * Permission to use, copy, modify, and distribute this software for any
8 * purpose without fee is hereby granted, provided that this entire notice
9 * is included in all copies of any software which is or includes a copy
10 * or modification of this software and in all copies of the supporting
11 * documentation for such software.
13 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
14 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
15 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
16 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
18 ***************************************************************/
20 /****************************************************************
21 * This is dtoa.c by David M. Gay, downloaded from
22 * http://www.netlib.org/fp/dtoa.c on April 15, 2009 and modified for
23 * inclusion into the Python core by Mark E. T. Dickinson and Eric V. Smith.
25 * Please remember to check http://www.netlib.org/fp regularly (and especially
26 * before any Python release) for bugfixes and updates.
28 * The major modifications from Gay's original code are as follows:
30 * 0. The original code has been specialized to Python's needs by removing
31 * many of the #ifdef'd sections. In particular, code to support VAX and
32 * IBM floating-point formats, hex NaNs, hex floats, locale-aware
33 * treatment of the decimal point, and setting of the inexact flag have
36 * 1. We use PyMem_Malloc and PyMem_Free in place of malloc and free.
38 * 2. The public functions strtod, dtoa and freedtoa all now have
41 * 3. Instead of assuming that PyMem_Malloc always succeeds, we thread
42 * PyMem_Malloc failures through the code. The functions
44 * Balloc, multadd, s2b, i2b, mult, pow5mult, lshift, diff, d2b
46 * of return type *Bigint all return NULL to indicate a malloc failure.
47 * Similarly, rv_alloc and nrv_alloc (return type char *) return NULL on
48 * failure. bigcomp now has return type int (it used to be void) and
49 * returns -1 on failure and 0 otherwise. _Py_dg_dtoa returns NULL
50 * on failure. _Py_dg_strtod indicates failure due to malloc failure
51 * by returning -1.0, setting errno=ENOMEM and *se to s00.
53 * 4. The static variable dtoa_result has been removed. Callers of
54 * _Py_dg_dtoa are expected to call _Py_dg_freedtoa to free
55 * the memory allocated by _Py_dg_dtoa.
57 * 5. The code has been reformatted to better fit with Python's
58 * C style guide (PEP 7).
60 * 6. A bug in the memory allocation has been fixed: to avoid FREEing memory
61 * that hasn't been MALLOC'ed, private_mem should only be used when k <=
64 * 7. _Py_dg_strtod has been modified so that it doesn't accept strings with
67 ***************************************************************/
69 /* Please send bug reports for the original dtoa.c code to David M. Gay (dmg
70 * at acm dot org, with " at " changed at "@" and " dot " changed to ".").
71 * Please report bugs for this modified version using the Python issue tracker
72 * (http://bugs.python.org). */
74 /* On a machine with IEEE extended-precision registers, it is
75 * necessary to specify double-precision (53-bit) rounding precision
76 * before invoking strtod or dtoa. If the machine uses (the equivalent
77 * of) Intel 80x87 arithmetic, the call
78 * _control87(PC_53, MCW_PC);
79 * does this with many compilers. Whether this or another call is
80 * appropriate depends on the compiler; for this to work, it may be
81 * necessary to #include "float.h" or another system-dependent header
85 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
87 * This strtod returns a nearest machine number to the input decimal
88 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
89 * broken by the IEEE round-even rule. Otherwise ties are broken by
90 * biased rounding (add half and chop).
92 * Inspired loosely by William D. Clinger's paper "How to Read Floating
93 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
97 * 1. We only require IEEE, IBM, or VAX double-precision
98 * arithmetic (not IEEE double-extended).
99 * 2. We get by with floating-point arithmetic in a case that
100 * Clinger missed -- when we're computing d * 10^n
101 * for a small integer d and the integer n is not too
102 * much larger than 22 (the maximum integer k for which
103 * we can represent 10^k exactly), we may be able to
104 * compute (d*10^k) * 10^(e-k) with just one roundoff.
105 * 3. Rather than a bit-at-a-time adjustment of the binary
106 * result in the hard case, we use floating-point
107 * arithmetic to determine the adjustment to within
108 * one bit; only in really hard cases do we need to
109 * compute a second residual.
110 * 4. Because of 3., we don't need a large table of powers of 10
111 * for ten-to-e (just some small tables, e.g. of 10^k
115 /* Linking of Python's #defines to Gay's #defines starts here. */
119 /* if PY_NO_SHORT_FLOAT_REPR is defined, then don't even try to compile
120 the following code */
121 #ifndef PY_NO_SHORT_FLOAT_REPR
125 #define MALLOC PyMem_Malloc
126 #define FREE PyMem_Free
128 /* This code should also work for ARM mixed-endian format on little-endian
129 machines, where doubles have byte order 45670123 (in increasing address
130 order, 0 being the least significant byte). */
131 #ifdef DOUBLE_IS_LITTLE_ENDIAN_IEEE754
134 #if defined(DOUBLE_IS_BIG_ENDIAN_IEEE754) || \
135 defined(DOUBLE_IS_ARM_MIXED_ENDIAN_IEEE754)
138 #if defined(IEEE_8087) + defined(IEEE_MC68k) != 1
139 #error "Exactly one of IEEE_8087 or IEEE_MC68k should be defined."
142 /* The code below assumes that the endianness of integers matches the
143 endianness of the two 32-bit words of a double. Check this. */
144 #if defined(WORDS_BIGENDIAN) && (defined(DOUBLE_IS_LITTLE_ENDIAN_IEEE754) || \
145 defined(DOUBLE_IS_ARM_MIXED_ENDIAN_IEEE754))
146 #error "doubles and ints have incompatible endianness"
149 #if !defined(WORDS_BIGENDIAN) && defined(DOUBLE_IS_BIG_ENDIAN_IEEE754)
150 #error "doubles and ints have incompatible endianness"
154 #if defined(HAVE_UINT32_T) && defined(HAVE_INT32_T)
155 typedef PY_UINT32_T ULong
;
156 typedef PY_INT32_T Long
;
158 #error "Failed to find an exact-width 32-bit integer type"
161 #if defined(HAVE_UINT64_T)
162 #define ULLong PY_UINT64_T
172 /* End Python #define linking */
175 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
179 #define PRIVATE_MEM 2304
181 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
182 static double private_mem
[PRIVATE_mem
], *pmem_next
= private_mem
;
188 typedef union { double d
; ULong L
[2]; } U
;
191 #define word0(x) (x)->L[1]
192 #define word1(x) (x)->L[0]
194 #define word0(x) (x)->L[0]
195 #define word1(x) (x)->L[1]
197 #define dval(x) (x)->d
199 #ifndef STRTOD_DIGLIM
200 #define STRTOD_DIGLIM 40
203 /* maximum permitted exponent value for strtod; exponents larger than
204 MAX_ABS_EXP in absolute value get truncated to +-MAX_ABS_EXP. MAX_ABS_EXP
205 should fit into an int. */
207 #define MAX_ABS_EXP 19999U
210 /* The following definition of Storeinc is appropriate for MIPS processors.
211 * An alternative that might be better on some machines is
212 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
214 #if defined(IEEE_8087)
215 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
216 ((unsigned short *)a)[0] = (unsigned short)c, a++)
218 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
219 ((unsigned short *)a)[1] = (unsigned short)c, a++)
222 /* #define P DBL_MANT_DIG */
223 /* Ten_pmax = floor(P*log(2)/log(5)) */
224 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
225 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
226 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
229 #define Exp_shift1 20
230 #define Exp_msk1 0x100000
231 #define Exp_msk11 0x100000
232 #define Exp_mask 0x7ff00000
238 #define Exp_1 0x3ff00000
239 #define Exp_11 0x3ff00000
241 #define Frac_mask 0xfffff
242 #define Frac_mask1 0xfffff
245 #define Bndry_mask 0xfffff
246 #define Bndry_mask1 0xfffff
248 #define Sign_bit 0x80000000
257 #define Flt_Rounds FLT_ROUNDS
261 #endif /*Flt_Rounds*/
263 #define Rounding Flt_Rounds
265 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
266 #define Big1 0xffffffff
268 /* struct BCinfo is used to pass information from _Py_dg_strtod to bigcomp */
270 typedef struct BCinfo BCinfo
;
273 int dp0
, dp1
, dplen
, dsign
, e0
, nd
, nd0
, scale
;
276 #define FFFFFFFF 0xffffffffUL
280 /* struct Bigint is used to represent arbitrary-precision integers. These
281 integers are stored in sign-magnitude format, with the magnitude stored as
282 an array of base 2**32 digits. Bigints are always normalized: if x is a
283 Bigint then x->wds >= 1, and either x->wds == 1 or x[wds-1] is nonzero.
285 The Bigint fields are as follows:
287 - next is a header used by Balloc and Bfree to keep track of lists
288 of freed Bigints; it's also used for the linked list of
289 powers of 5 of the form 5**2**i used by pow5mult.
290 - k indicates which pool this Bigint was allocated from
291 - maxwds is the maximum number of words space was allocated for
292 (usually maxwds == 2**k)
293 - sign is 1 for negative Bigints, 0 for positive. The sign is unused
294 (ignored on inputs, set to 0 on outputs) in almost all operations
295 involving Bigints: a notable exception is the diff function, which
296 ignores signs on inputs but sets the sign of the output correctly.
297 - wds is the actual number of significant words
298 - x contains the vector of words (digits) for this Bigint, from least
299 significant (x[0]) to most significant (x[wds-1]).
305 int k
, maxwds
, sign
, wds
;
309 typedef struct Bigint Bigint
;
311 /* Memory management: memory is allocated from, and returned to, Kmax+1 pools
312 of memory, where pool k (0 <= k <= Kmax) is for Bigints b with b->maxwds ==
313 1 << k. These pools are maintained as linked lists, with freelist[k]
314 pointing to the head of the list for pool k.
316 On allocation, if there's no free slot in the appropriate pool, MALLOC is
317 called to get more memory. This memory is not returned to the system until
318 Python quits. There's also a private memory pool that's allocated from
319 in preference to using MALLOC.
321 For Bigints with more than (1 << Kmax) digits (which implies at least 1233
322 decimal digits), memory is directly allocated using MALLOC, and freed using
325 XXX: it would be easy to bypass this memory-management system and
326 translate each call to Balloc into a call to PyMem_Malloc, and each
327 Bfree to PyMem_Free. Investigate whether this has any significant
328 performance on impact. */
330 static Bigint
*freelist
[Kmax
+1];
332 /* Allocate space for a Bigint with up to 1<<k digits */
341 if (k
<= Kmax
&& (rv
= freelist
[k
]))
342 freelist
[k
] = rv
->next
;
345 len
= (sizeof(Bigint
) + (x
-1)*sizeof(ULong
) + sizeof(double) - 1)
347 if (k
<= Kmax
&& pmem_next
- private_mem
+ len
<= PRIVATE_mem
) {
348 rv
= (Bigint
*)pmem_next
;
352 rv
= (Bigint
*)MALLOC(len
*sizeof(double));
359 rv
->sign
= rv
->wds
= 0;
363 /* Free a Bigint allocated with Balloc */
372 v
->next
= freelist
[v
->k
];
378 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
379 y->wds*sizeof(Long) + 2*sizeof(int))
381 /* Multiply a Bigint b by m and add a. Either modifies b in place and returns
382 a pointer to the modified b, or Bfrees b and returns a pointer to a copy.
383 On failure, return NULL. In this case, b will have been already freed. */
386 multadd(Bigint
*b
, int m
, int a
) /* multiply by m and add a */
404 y
= *x
* (ULLong
)m
+ carry
;
406 *x
++ = (ULong
)(y
& FFFFFFFF
);
409 y
= (xi
& 0xffff) * m
+ carry
;
410 z
= (xi
>> 16) * m
+ (y
>> 16);
412 *x
++ = (z
<< 16) + (y
& 0xffff);
417 if (wds
>= b
->maxwds
) {
427 b
->x
[wds
++] = (ULong
)carry
;
433 /* convert a string s containing nd decimal digits (possibly containing a
434 decimal separator at position nd0, which is ignored) to a Bigint. This
435 function carries on where the parsing code in _Py_dg_strtod leaves off: on
436 entry, y9 contains the result of converting the first 9 digits. Returns
440 s2b(const char *s
, int nd0
, int nd
, ULong y9
, int dplen
)
447 for(k
= 0, y
= 1; x
> y
; y
<<= 1, k
++) ;
458 b
= multadd(b
, 10, *s
++ - '0');
467 b
= multadd(b
, 10, *s
++ - '0');
474 /* count leading 0 bits in the 32-bit integer x. */
481 if (!(x
& 0xffff0000)) {
485 if (!(x
& 0xff000000)) {
489 if (!(x
& 0xf0000000)) {
493 if (!(x
& 0xc0000000)) {
497 if (!(x
& 0x80000000)) {
499 if (!(x
& 0x40000000))
505 /* count trailing 0 bits in the 32-bit integer y, and shift y right by that
551 /* convert a small nonnegative integer to a Bigint */
566 /* multiply two Bigints. Returns a new Bigint, or NULL on failure. Ignores
567 the signs of a and b. */
570 mult(Bigint
*a
, Bigint
*b
)
574 ULong
*x
, *xa
, *xae
, *xb
, *xbe
, *xc
, *xc0
;
583 if (a
->wds
< b
->wds
) {
597 for(x
= c
->x
, xa
= x
+ wc
; x
< xa
; x
++)
605 for(; xb
< xbe
; xc0
++) {
611 z
= *x
++ * (ULLong
)y
+ *xc
+ carry
;
613 *xc
++ = (ULong
)(z
& FFFFFFFF
);
620 for(; xb
< xbe
; xb
++, xc0
++) {
621 if (y
= *xb
& 0xffff) {
626 z
= (*x
& 0xffff) * y
+ (*xc
& 0xffff) + carry
;
628 z2
= (*x
++ >> 16) * y
+ (*xc
>> 16) + carry
;
641 z
= (*x
& 0xffff) * y
+ (*xc
>> 16) + carry
;
644 z2
= (*x
++ >> 16) * y
+ (*xc
& 0xffff) + carry
;
652 for(xc0
= c
->x
, xc
= xc0
+ wc
; wc
> 0 && !*--xc
; --wc
) ;
657 /* p5s is a linked list of powers of 5 of the form 5**(2**i), i >= 2 */
661 /* multiply the Bigint b by 5**k. Returns a pointer to the result, or NULL on
662 failure; if the returned pointer is distinct from b then the original
663 Bigint b will have been Bfree'd. Ignores the sign of b. */
666 pow5mult(Bigint
*b
, int k
)
668 Bigint
*b1
, *p5
, *p51
;
670 static int p05
[3] = { 5, 25, 125 };
673 b
= multadd(b
, p05
[i
-1], 0);
716 /* shift a Bigint b left by k bits. Return a pointer to the shifted result,
717 or NULL on failure. If the returned pointer is distinct from b then the
718 original b will have been Bfree'd. Ignores the sign of b. */
721 lshift(Bigint
*b
, int k
)
725 ULong
*x
, *x1
, *xe
, z
;
730 for(i
= b
->maxwds
; n1
> i
; i
<<= 1)
738 for(i
= 0; i
< n
; i
++)
761 /* Do a three-way compare of a and b, returning -1 if a < b, 0 if a == b and
762 1 if a > b. Ignores signs of a and b. */
765 cmp(Bigint
*a
, Bigint
*b
)
767 ULong
*xa
, *xa0
, *xb
, *xb0
;
773 if (i
> 1 && !a
->x
[i
-1])
774 Bug("cmp called with a->x[a->wds-1] == 0");
775 if (j
> 1 && !b
->x
[j
-1])
776 Bug("cmp called with b->x[b->wds-1] == 0");
786 return *xa
< *xb
? -1 : 1;
793 /* Take the difference of Bigints a and b, returning a new Bigint. Returns
794 NULL on failure. The signs of a and b are ignored, but the sign of the
795 result is set appropriately. */
798 diff(Bigint
*a
, Bigint
*b
)
802 ULong
*xa
, *xae
, *xb
, *xbe
, *xc
;
841 y
= (ULLong
)*xa
++ - *xb
++ - borrow
;
842 borrow
= y
>> 32 & (ULong
)1;
843 *xc
++ = (ULong
)(y
& FFFFFFFF
);
848 borrow
= y
>> 32 & (ULong
)1;
849 *xc
++ = (ULong
)(y
& FFFFFFFF
);
853 y
= (*xa
& 0xffff) - (*xb
& 0xffff) - borrow
;
854 borrow
= (y
& 0x10000) >> 16;
855 z
= (*xa
++ >> 16) - (*xb
++ >> 16) - borrow
;
856 borrow
= (z
& 0x10000) >> 16;
861 y
= (*xa
& 0xffff) - borrow
;
862 borrow
= (y
& 0x10000) >> 16;
863 z
= (*xa
++ >> 16) - borrow
;
864 borrow
= (z
& 0x10000) >> 16;
874 /* Given a positive normal double x, return the difference between x and the next
875 double up. Doesn't give correct results for subnormals. */
883 L
= (word0(x
) & Exp_mask
) - (P
-1)*Exp_msk1
;
889 /* Convert a Bigint to a double plus an exponent */
892 b2d(Bigint
*a
, int *e
)
894 ULong
*xa
, *xa0
, w
, y
, z
;
902 if (!y
) Bug("zero y in b2d");
907 word0(&d
) = Exp_1
| y
>> (Ebits
- k
);
908 w
= xa
> xa0
? *--xa
: 0;
909 word1(&d
) = y
<< ((32-Ebits
) + k
) | w
>> (Ebits
- k
);
912 z
= xa
> xa0
? *--xa
: 0;
914 word0(&d
) = Exp_1
| y
<< k
| z
>> (32 - k
);
915 y
= xa
> xa0
? *--xa
: 0;
916 word1(&d
) = z
<< k
| y
>> (32 - k
);
919 word0(&d
) = Exp_1
| y
;
926 /* Convert a double to a Bigint plus an exponent. Return NULL on failure.
928 Given a finite nonzero double d, return an odd Bigint b and exponent *e
929 such that fabs(d) = b * 2**e. On return, *bbits gives the number of
930 significant bits of b; that is, 2**(*bbits-1) <= b < 2**(*bbits).
932 If d is zero, then b == 0, *e == -1010, *bbits = 0.
937 d2b(U
*d
, int *e
, int *bits
)
949 z
= word0(d
) & Frac_mask
;
950 word0(d
) &= 0x7fffffff; /* clear sign bit, which we ignore */
951 if ((de
= (int)(word0(d
) >> Exp_shift
)))
953 if ((y
= word1(d
))) {
954 if ((k
= lo0bits(&y
))) {
955 x
[0] = y
| z
<< (32 - k
);
961 b
->wds
= (x
[1] = z
) ? 2 : 1;
971 *e
= de
- Bias
- (P
-1) + k
;
975 *e
= de
- Bias
- (P
-1) + 1 + k
;
976 *bits
= 32*i
- hi0bits(x
[i
-1]);
981 /* Compute the ratio of two Bigints, as a double. The result may have an
982 error of up to 2.5 ulps. */
985 ratio(Bigint
*a
, Bigint
*b
)
990 dval(&da
) = b2d(a
, &ka
);
991 dval(&db
) = b2d(b
, &kb
);
992 k
= ka
- kb
+ 32*(a
->wds
- b
->wds
);
994 word0(&da
) += k
*Exp_msk1
;
997 word0(&db
) += k
*Exp_msk1
;
999 return dval(&da
) / dval(&db
);
1004 1e0
, 1e1
, 1e2
, 1e3
, 1e4
, 1e5
, 1e6
, 1e7
, 1e8
, 1e9
,
1005 1e10
, 1e11
, 1e12
, 1e13
, 1e14
, 1e15
, 1e16
, 1e17
, 1e18
, 1e19
,
1010 bigtens
[] = { 1e16
, 1e32
, 1e64
, 1e128
, 1e256
};
1011 static const double tinytens
[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1012 9007199254740992.*9007199254740992.e
-256
1013 /* = 2^106 * 1e-256 */
1015 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1016 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1017 #define Scale_Bit 0x10
1026 dshift(Bigint
*b
, int p2
)
1028 int rv
= hi0bits(b
->x
[b
->wds
-1]) - 4;
1034 /* special case of Bigint division. The quotient is always in the range 0 <=
1035 quotient < 10, and on entry the divisor S is normalized so that its top 4
1036 bits (28--31) are zero and bit 27 is set. */
1039 quorem(Bigint
*b
, Bigint
*S
)
1042 ULong
*bx
, *bxe
, q
, *sx
, *sxe
;
1044 ULLong borrow
, carry
, y
, ys
;
1046 ULong borrow
, carry
, y
, ys
;
1052 /*debug*/ if (b
->wds
> n
)
1053 /*debug*/ Bug("oversize b in quorem");
1061 q
= *bxe
/ (*sxe
+ 1); /* ensure q <= true quotient */
1063 /*debug*/ if (q
> 9)
1064 /*debug*/ Bug("oversized quotient in quorem");
1071 ys
= *sx
++ * (ULLong
)q
+ carry
;
1073 y
= *bx
- (ys
& FFFFFFFF
) - borrow
;
1074 borrow
= y
>> 32 & (ULong
)1;
1075 *bx
++ = (ULong
)(y
& FFFFFFFF
);
1078 ys
= (si
& 0xffff) * q
+ carry
;
1079 zs
= (si
>> 16) * q
+ (ys
>> 16);
1081 y
= (*bx
& 0xffff) - (ys
& 0xffff) - borrow
;
1082 borrow
= (y
& 0x10000) >> 16;
1083 z
= (*bx
>> 16) - (zs
& 0xffff) - borrow
;
1084 borrow
= (z
& 0x10000) >> 16;
1091 while(--bxe
> bx
&& !*bxe
)
1096 if (cmp(b
, S
) >= 0) {
1106 y
= *bx
- (ys
& FFFFFFFF
) - borrow
;
1107 borrow
= y
>> 32 & (ULong
)1;
1108 *bx
++ = (ULong
)(y
& FFFFFFFF
);
1111 ys
= (si
& 0xffff) + carry
;
1112 zs
= (si
>> 16) + (ys
>> 16);
1114 y
= (*bx
& 0xffff) - (ys
& 0xffff) - borrow
;
1115 borrow
= (y
& 0x10000) >> 16;
1116 z
= (*bx
>> 16) - (zs
& 0xffff) - borrow
;
1117 borrow
= (z
& 0x10000) >> 16;
1125 while(--bxe
> bx
&& !*bxe
)
1134 /* return 0 on success, -1 on failure */
1137 bigcomp(U
*rv
, const char *s0
, BCinfo
*bc
)
1140 int b2
, bbits
, d2
, dd
, dig
, dsign
, i
, j
, nd
, nd0
, p2
, p5
, speccase
;
1145 p5
= nd
+ bc
->e0
- 1;
1147 if (rv
->d
== 0.) { /* special case: value near underflow-to-zero */
1148 /* threshold was rounded to zero */
1154 word0(rv
) = (P
+2) << Exp_shift
;
1165 b
= d2b(rv
, &p2
, &bbits
);
1170 /* floor(log2(rv)) == bbits - 1 + p2 */
1171 /* Check for denormal case. */
1173 if (i
> (j
= P
- Emin
- 1 + p2
)) {
1189 /* Arrange for convenient computation of quotients:
1190 * shift left if necessary so divisor has 4 leading 0 bits.
1193 d
= pow5mult(d
, p5
);
1200 b
= pow5mult(b
, -p5
);
1215 if ((b2
+= i
) > 0) {
1222 if ((d2
+= i
) > 0) {
1230 /* Now b/d = exactly half-way between the two floating-point values */
1231 /* on either side of the input string. Compute first digit of b/d. */
1233 if (!(dig
= quorem(b
,d
))) {
1234 b
= multadd(b
, 10, 0); /* very unlikely */
1242 /* Compare b/d with s0 */
1245 dd
= 9999; /* silence gcc compiler warning */
1246 for(i
= 0; i
< nd0
; ) {
1247 if ((dd
= s0
[i
++] - '0' - dig
))
1249 if (!b
->x
[0] && b
->wds
== 1) {
1254 b
= multadd(b
, 10, 0);
1261 for(j
= bc
->dp1
; i
++ < nd
;) {
1262 if ((dd
= s0
[j
++] - '0' - dig
))
1264 if (!b
->x
[0] && b
->wds
== 1) {
1269 b
= multadd(b
, 10, 0);
1276 if (b
->x
[0] || b
->wds
> 1)
1286 if (!dsign
) /* does not happen for round-near */
1288 dval(rv
) -= ulp(rv
);
1293 dval(rv
) += ulp(rv
);
1297 /* Exact half-way case: apply round-even rule. */
1298 if (word1(rv
) & 1) {
1309 _Py_dg_strtod(const char *s00
, char **se
)
1311 int bb2
, bb5
, bbe
, bd2
, bd5
, bbbits
, bs2
, c
, e
, e1
, error
;
1312 int esign
, i
, j
, k
, nd
, nd0
, nf
, nz
, nz0
, sign
;
1313 const char *s
, *s0
, *s1
;
1315 U aadj2
, adj
, rv
, rv0
;
1318 Bigint
*bb
, *bb1
, *bd
, *bd0
, *bs
, *delta
;
1320 sign
= nz0
= nz
= bc
.dplen
= 0;
1322 for(s
= s00
;;s
++) switch(*s
) {
1332 /* modify original dtoa.c so that it doesn't accept leading whitespace
1347 while(*++s
== '0') ;
1353 for(nd
= nf
= 0; (c
= *s
) >= '0' && c
<= '9'; nd
++, s
++)
1359 bc
.dp0
= bc
.dp1
= s
- s0
;
1363 bc
.dplen
= bc
.dp1
- bc
.dp0
;
1365 for(; c
== '0'; c
= *++s
)
1367 if (c
> '0' && c
<= '9') {
1375 for(; c
>= '0' && c
<= '9'; c
= *++s
) {
1380 for(i
= 1; i
< nz
; i
++)
1383 else if (nd
<= DBL_DIG
+ 1)
1387 else if (nd
<= DBL_DIG
+ 1)
1395 if (c
== 'e' || c
== 'E') {
1396 if (!nd
&& !nz
&& !nz0
) {
1407 if (c
>= '0' && c
<= '9') {
1410 if (c
> '0' && c
<= '9') {
1413 while((c
= *++s
) >= '0' && c
<= '9')
1415 if (s
- s1
> 8 || L
> MAX_ABS_EXP
)
1416 /* Avoid confusion from exponents
1417 * so large that e might overflow.
1419 e
= (int)MAX_ABS_EXP
; /* safe for 16 bit ints */
1439 bc
.e0
= e1
= e
-= nf
;
1441 /* Now we have nd0 digits, starting at s0, followed by a
1442 * decimal point, followed by nd-nd0 digits. The number we're
1443 * after is the integer represented by those digits times
1448 k
= nd
< DBL_DIG
+ 1 ? nd
: DBL_DIG
+ 1;
1451 dval(&rv
) = tens
[k
- 9] * dval(&rv
) + z
;
1460 if (e
<= Ten_pmax
) {
1461 dval(&rv
) *= tens
[e
];
1465 if (e
<= Ten_pmax
+ i
) {
1466 /* A fancier test would sometimes let us do
1467 * this for larger i values.
1470 dval(&rv
) *= tens
[i
];
1471 dval(&rv
) *= tens
[e
];
1475 else if (e
>= -Ten_pmax
) {
1476 dval(&rv
) /= tens
[-e
];
1484 /* Get starting approximation = rv * 10**e1 */
1488 dval(&rv
) *= tens
[i
];
1490 if (e1
> DBL_MAX_10_EXP
) {
1493 /* Can't trust HUGE_VAL */
1494 word0(&rv
) = Exp_mask
;
1499 for(j
= 0; e1
> 1; j
++, e1
>>= 1)
1501 dval(&rv
) *= bigtens
[j
];
1502 /* The last multiplication could overflow. */
1503 word0(&rv
) -= P
*Exp_msk1
;
1504 dval(&rv
) *= bigtens
[j
];
1505 if ((z
= word0(&rv
) & Exp_mask
)
1506 > Exp_msk1
*(DBL_MAX_EXP
+Bias
-P
))
1508 if (z
> Exp_msk1
*(DBL_MAX_EXP
+Bias
-1-P
)) {
1509 /* set to largest number */
1510 /* (Can't trust DBL_MAX) */
1515 word0(&rv
) += P
*Exp_msk1
;
1521 dval(&rv
) /= tens
[i
];
1523 if (e1
>= 1 << n_bigtens
)
1527 for(j
= 0; e1
> 0; j
++, e1
>>= 1)
1529 dval(&rv
) *= tinytens
[j
];
1530 if (bc
.scale
&& (j
= 2*P
+ 1 - ((word0(&rv
) & Exp_mask
)
1531 >> Exp_shift
)) > 0) {
1532 /* scaled rv is denormal; clear j low bits */
1536 word0(&rv
) = (P
+2)*Exp_msk1
;
1538 word0(&rv
) &= 0xffffffff << (j
-32);
1541 word1(&rv
) &= 0xffffffff << j
;
1552 /* Now the hard part -- adjusting rv to the correct value.*/
1554 /* Put digits into bd: true value = bd * 10^e */
1557 bc
.nd0
= nd0
; /* Only needed if nd > STRTOD_DIGLIM, but done here */
1558 /* to silence an erroneous warning about bc.nd0 */
1559 /* possibly not being initialized. */
1560 if (nd
> STRTOD_DIGLIM
) {
1561 /* ASSERT(STRTOD_DIGLIM >= 18); 18 == one more than the */
1562 /* minimum number of decimal digits to distinguish double values */
1563 /* in IEEE arithmetic. */
1568 if (--j
<= bc
.dp1
&& j
>= bc
.dp0
)
1578 if (nd
< 9) { /* must recompute y */
1580 for(i
= 0; i
< nd0
; ++i
)
1581 y
= 10*y
+ s0
[i
] - '0';
1582 for(j
= bc
.dp1
; i
< nd
; ++i
)
1583 y
= 10*y
+ s0
[j
++] - '0';
1586 bd0
= s2b(s0
, nd0
, nd
, y
, bc
.dplen
);
1591 bd
= Balloc(bd0
->k
);
1597 bb
= d2b(&rv
, &bbe
, &bbbits
); /* rv = bb * 2^bbe */
1625 i
= j
+ bbbits
- 1; /* logb(rv) */
1626 if (i
< Emin
) /* denormal */
1633 i
= bb2
< bd2
? bb2
: bd2
;
1642 bs
= pow5mult(bs
, bb5
);
1660 bb
= lshift(bb
, bb2
);
1669 bd
= pow5mult(bd
, bd5
);
1678 bd
= lshift(bd
, bd2
);
1687 bs
= lshift(bs
, bs2
);
1695 delta
= diff(bb
, bd
);
1696 if (delta
== NULL
) {
1703 bc
.dsign
= delta
->sign
;
1706 if (bc
.nd
> nd
&& i
<= 0) {
1708 break; /* Must use bigcomp(). */
1711 i
= -1; /* Discarded digits make delta smaller. */
1716 /* Error is less than half an ulp -- check for
1717 * special case of mantissa a power of two.
1719 if (bc
.dsign
|| word1(&rv
) || word0(&rv
) & Bndry_mask
1720 || (word0(&rv
) & Exp_mask
) <= (2*P
+1)*Exp_msk1
1724 if (!delta
->x
[0] && delta
->wds
<= 1) {
1728 delta
= lshift(delta
,Log2P
);
1729 if (delta
== NULL
) {
1736 if (cmp(delta
, bs
) > 0)
1741 /* exactly half-way between */
1743 if ((word0(&rv
) & Bndry_mask1
) == Bndry_mask1
1746 (y
= word0(&rv
) & Exp_mask
) <= 2*P
*Exp_msk1
) ?
1747 (0xffffffff & (0xffffffff << (2*P
+1-(y
>>Exp_shift
)))) :
1749 /*boundary case -- increment exponent*/
1750 word0(&rv
) = (word0(&rv
) & Exp_mask
)
1758 else if (!(word0(&rv
) & Bndry_mask
) && !word1(&rv
)) {
1760 /* boundary case -- decrement exponent */
1762 L
= word0(&rv
) & Exp_mask
;
1763 if (L
<= (2*P
+1)*Exp_msk1
) {
1764 if (L
> (P
+2)*Exp_msk1
)
1765 /* round even ==> */
1768 /* rv = smallest denormal */
1774 L
= (word0(&rv
) & Exp_mask
) - Exp_msk1
;
1775 word0(&rv
) = L
| Bndry_mask1
;
1776 word1(&rv
) = 0xffffffff;
1779 if (!(word1(&rv
) & LSB
))
1782 dval(&rv
) += ulp(&rv
);
1784 dval(&rv
) -= ulp(&rv
);
1791 bc
.dsign
= 1 - bc
.dsign
;
1794 if ((aadj
= ratio(delta
, bs
)) <= 2.) {
1797 else if (word1(&rv
) || word0(&rv
) & Bndry_mask
) {
1798 if (word1(&rv
) == Tiny1
&& !word0(&rv
)) {
1807 /* special case -- power of FLT_RADIX to be */
1808 /* rounded down... */
1810 if (aadj
< 2./FLT_RADIX
)
1811 aadj
= 1./FLT_RADIX
;
1819 aadj1
= bc
.dsign
? aadj
: -aadj
;
1820 if (Flt_Rounds
== 0)
1823 y
= word0(&rv
) & Exp_mask
;
1825 /* Check for overflow */
1827 if (y
== Exp_msk1
*(DBL_MAX_EXP
+Bias
-1)) {
1828 dval(&rv0
) = dval(&rv
);
1829 word0(&rv
) -= P
*Exp_msk1
;
1830 adj
.d
= aadj1
* ulp(&rv
);
1832 if ((word0(&rv
) & Exp_mask
) >=
1833 Exp_msk1
*(DBL_MAX_EXP
+Bias
-P
)) {
1834 if (word0(&rv0
) == Big0
&& word1(&rv0
) == Big1
)
1841 word0(&rv
) += P
*Exp_msk1
;
1844 if (bc
.scale
&& y
<= 2*P
*Exp_msk1
) {
1845 if (aadj
<= 0x7fffffff) {
1846 if ((z
= (ULong
)aadj
) <= 0)
1849 aadj1
= bc
.dsign
? aadj
: -aadj
;
1851 dval(&aadj2
) = aadj1
;
1852 word0(&aadj2
) += (2*P
+1)*Exp_msk1
- y
;
1853 aadj1
= dval(&aadj2
);
1855 adj
.d
= aadj1
* ulp(&rv
);
1858 z
= word0(&rv
) & Exp_mask
;
1862 /* Can we stop now? */
1865 /* The tolerances below are conservative. */
1866 if (bc
.dsign
|| word1(&rv
) || word0(&rv
) & Bndry_mask
) {
1867 if (aadj
< .4999999 || aadj
> .5000001)
1870 else if (aadj
< .4999999/FLT_RADIX
)
1886 error
= bigcomp(&rv
, s0
, &bc
);
1892 word0(&rv0
) = Exp_1
- 2*P
*Exp_msk1
;
1894 dval(&rv
) *= dval(&rv0
);
1895 /* try to avoid the bug of testing an 8087 register value */
1896 if (!(word0(&rv
) & Exp_mask
))
1902 return sign
? -dval(&rv
) : dval(&rv
);
1918 sizeof(Bigint
) - sizeof(ULong
) - sizeof(int) + j
<= (unsigned)i
;
1921 r
= (int*)Balloc(k
);
1925 return (char *)(r
+1);
1929 nrv_alloc(char *s
, char **rve
, int n
)
1937 while((*t
= *s
++)) t
++;
1943 /* freedtoa(s) must be used to free values s returned by dtoa
1944 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
1945 * but for consistency with earlier versions of dtoa, it is optional
1946 * when MULTIPLE_THREADS is not defined.
1950 _Py_dg_freedtoa(char *s
)
1952 Bigint
*b
= (Bigint
*)((int *)s
- 1);
1953 b
->maxwds
= 1 << (b
->k
= *(int*)b
);
1957 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1959 * Inspired by "How to Print Floating-Point Numbers Accurately" by
1960 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
1963 * 1. Rather than iterating, we use a simple numeric overestimate
1964 * to determine k = floor(log10(d)). We scale relevant
1965 * quantities using O(log2(k)) rather than O(k) multiplications.
1966 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1967 * try to generate digits strictly left to right. Instead, we
1968 * compute with fewer bits and propagate the carry if necessary
1969 * when rounding the final digit up. This is often faster.
1970 * 3. Under the assumption that input will be rounded nearest,
1971 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1972 * That is, we allow equality in stopping tests when the
1973 * round-nearest rule will give the same floating-point value
1974 * as would satisfaction of the stopping test with strict
1976 * 4. We remove common factors of powers of 2 from relevant
1978 * 5. When converting floating-point integers less than 1e16,
1979 * we use floating-point arithmetic rather than resorting
1980 * to multiple-precision integers.
1981 * 6. When asked to produce fewer than 15 digits, we first try
1982 * to get by with floating-point arithmetic; we resort to
1983 * multiple-precision integer arithmetic only if we cannot
1984 * guarantee that the floating-point calculation has given
1985 * the correctly rounded result. For k requested digits and
1986 * "uniformly" distributed input, the probability is
1987 * something like 10^(k-15) that we must resort to the Long
1991 /* Additional notes (METD): (1) returns NULL on failure. (2) to avoid memory
1992 leakage, a successful call to _Py_dg_dtoa should always be matched by a
1993 call to _Py_dg_freedtoa. */
1996 _Py_dg_dtoa(double dd
, int mode
, int ndigits
,
1997 int *decpt
, int *sign
, char **rve
)
1999 /* Arguments ndigits, decpt, sign are similar to those
2000 of ecvt and fcvt; trailing zeros are suppressed from
2001 the returned string. If not null, *rve is set to point
2002 to the end of the return value. If d is +-Infinity or NaN,
2003 then *decpt is set to 9999.
2006 0 ==> shortest string that yields d when read in
2007 and rounded to nearest.
2008 1 ==> like 0, but with Steele & White stopping rule;
2009 e.g. with IEEE P754 arithmetic , mode 0 gives
2010 1e23 whereas mode 1 gives 9.999999999999999e22.
2011 2 ==> max(1,ndigits) significant digits. This gives a
2012 return value similar to that of ecvt, except
2013 that trailing zeros are suppressed.
2014 3 ==> through ndigits past the decimal point. This
2015 gives a return value similar to that from fcvt,
2016 except that trailing zeros are suppressed, and
2017 ndigits can be negative.
2018 4,5 ==> similar to 2 and 3, respectively, but (in
2019 round-nearest mode) with the tests of mode 0 to
2020 possibly return a shorter string that rounds to d.
2021 With IEEE arithmetic and compilation with
2022 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
2023 as modes 2 and 3 when FLT_ROUNDS != 1.
2024 6-9 ==> Debugging modes similar to mode - 4: don't try
2025 fast floating-point estimate (if applicable).
2027 Values of mode other than 0-9 are treated as mode 0.
2029 Sufficient space is allocated to the return value
2030 to hold the suppressed trailing zeros.
2033 int bbits
, b2
, b5
, be
, dig
, i
, ieps
, ilim
, ilim0
, ilim1
,
2034 j
, j1
, k
, k0
, k_check
, leftright
, m2
, m5
, s2
, s5
,
2035 spec_case
, try_quick
;
2039 Bigint
*b
, *b1
, *delta
, *mlo
, *mhi
, *S
;
2044 /* set pointers to NULL, to silence gcc compiler warnings and make
2045 cleanup easier on error */
2046 mlo
= mhi
= b
= S
= 0;
2050 if (word0(&u
) & Sign_bit
) {
2051 /* set sign for everything, including 0's and NaNs */
2053 word0(&u
) &= ~Sign_bit
; /* clear sign bit */
2058 /* quick return for Infinities, NaNs and zeros */
2059 if ((word0(&u
) & Exp_mask
) == Exp_mask
)
2061 /* Infinity or NaN */
2063 if (!word1(&u
) && !(word0(&u
) & 0xfffff))
2064 return nrv_alloc("Infinity", rve
, 8);
2065 return nrv_alloc("NaN", rve
, 3);
2069 return nrv_alloc("0", rve
, 1);
2072 /* compute k = floor(log10(d)). The computation may leave k
2073 one too large, but should never leave k too small. */
2074 b
= d2b(&u
, &be
, &bbits
);
2077 if ((i
= (int)(word0(&u
) >> Exp_shift1
& (Exp_mask
>>Exp_shift1
)))) {
2078 dval(&d2
) = dval(&u
);
2079 word0(&d2
) &= Frac_mask1
;
2080 word0(&d2
) |= Exp_11
;
2082 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2083 * log10(x) = log(x) / log(10)
2084 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2085 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2087 * This suggests computing an approximation k to log10(d) by
2089 * k = (i - Bias)*0.301029995663981
2090 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2092 * We want k to be too large rather than too small.
2093 * The error in the first-order Taylor series approximation
2094 * is in our favor, so we just round up the constant enough
2095 * to compensate for any error in the multiplication of
2096 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2097 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2098 * adding 1e-13 to the constant term more than suffices.
2099 * Hence we adjust the constant term to 0.1760912590558.
2100 * (We could get a more accurate k by invoking log10,
2101 * but this is probably not worthwhile.)
2108 /* d is denormalized */
2110 i
= bbits
+ be
+ (Bias
+ (P
-1) - 1);
2111 x
= i
> 32 ? word0(&u
) << (64 - i
) | word1(&u
) >> (i
- 32)
2112 : word1(&u
) << (32 - i
);
2114 word0(&d2
) -= 31*Exp_msk1
; /* adjust exponent */
2115 i
-= (Bias
+ (P
-1) - 1) + 1;
2118 ds
= (dval(&d2
)-1.5)*0.289529654602168 + 0.1760912590558 +
2119 i
*0.301029995663981;
2121 if (ds
< 0. && ds
!= k
)
2122 k
--; /* want k = floor(ds) */
2124 if (k
>= 0 && k
<= Ten_pmax
) {
2125 if (dval(&u
) < tens
[k
])
2148 if (mode
< 0 || mode
> 9)
2158 ilim
= ilim1
= -1; /* Values for cases 0 and 1; done here to */
2159 /* silence erroneous "gcc -Wall" warning. */
2172 ilim
= ilim1
= i
= ndigits
;
2178 i
= ndigits
+ k
+ 1;
2190 if (ilim
>= 0 && ilim
<= Quick_max
&& try_quick
) {
2192 /* Try to get by with floating-point arithmetic. */
2195 dval(&d2
) = dval(&u
);
2198 ieps
= 2; /* conservative */
2203 /* prevent overflows */
2205 dval(&u
) /= bigtens
[n_bigtens
-1];
2208 for(; j
; j
>>= 1, i
++)
2215 else if ((j1
= -k
)) {
2216 dval(&u
) *= tens
[j1
& 0xf];
2217 for(j
= j1
>> 4; j
; j
>>= 1, i
++)
2220 dval(&u
) *= bigtens
[i
];
2223 if (k_check
&& dval(&u
) < 1. && ilim
> 0) {
2231 dval(&eps
) = ieps
*dval(&u
) + 7.;
2232 word0(&eps
) -= (P
-1)*Exp_msk1
;
2236 if (dval(&u
) > dval(&eps
))
2238 if (dval(&u
) < -dval(&eps
))
2243 /* Use Steele & White method of only
2244 * generating digits needed.
2246 dval(&eps
) = 0.5/tens
[ilim
-1] - dval(&eps
);
2250 *s
++ = '0' + (int)L
;
2251 if (dval(&u
) < dval(&eps
))
2253 if (1. - dval(&u
) < dval(&eps
))
2262 /* Generate ilim digits, then fix them up. */
2263 dval(&eps
) *= tens
[ilim
-1];
2264 for(i
= 1;; i
++, dval(&u
) *= 10.) {
2265 L
= (Long
)(dval(&u
));
2266 if (!(dval(&u
) -= L
))
2268 *s
++ = '0' + (int)L
;
2270 if (dval(&u
) > 0.5 + dval(&eps
))
2272 else if (dval(&u
) < 0.5 - dval(&eps
)) {
2283 dval(&u
) = dval(&d2
);
2288 /* Do we have a "small" integer? */
2290 if (be
>= 0 && k
<= Int_max
) {
2293 if (ndigits
< 0 && ilim
<= 0) {
2295 if (ilim
< 0 || dval(&u
) <= 5*ds
)
2299 for(i
= 1;; i
++, dval(&u
) *= 10.) {
2300 L
= (Long
)(dval(&u
) / ds
);
2302 *s
++ = '0' + (int)L
;
2307 dval(&u
) += dval(&u
);
2308 if (dval(&u
) > ds
|| (dval(&u
) == ds
&& L
& 1)) {
2328 denorm
? be
+ (Bias
+ (P
-1) - 1 + 1) :
2336 if (m2
> 0 && s2
> 0) {
2337 i
= m2
< s2
? m2
: s2
;
2345 mhi
= pow5mult(mhi
, m5
);
2354 if ((j
= b5
- m5
)) {
2361 b
= pow5mult(b
, b5
);
2370 S
= pow5mult(S
, s5
);
2375 /* Check for special case that d is a normalized power of 2. */
2378 if ((mode
< 2 || leftright
)
2380 if (!word1(&u
) && !(word0(&u
) & Bndry_mask
)
2381 && word0(&u
) & (Exp_mask
& ~Exp_msk1
)
2383 /* The special case */
2390 /* Arrange for convenient computation of quotients:
2391 * shift left if necessary so divisor has 4 leading 0 bits.
2393 * Perhaps we should just compute leading 28 bits of S once
2394 * and for all and pass them and a shift to quorem, so it
2395 * can do shifts and ors to compute the numerator for q.
2397 if ((i
= ((s5
? 32 - hi0bits(S
->x
[S
->wds
-1]) : 1) + s2
) & 0x1f))
2417 b
= multadd(b
, 10, 0); /* we botched the k estimate */
2421 mhi
= multadd(mhi
, 10, 0);
2428 if (ilim
<= 0 && (mode
== 3 || mode
== 5)) {
2430 /* no digits, fcvt style */
2436 S
= multadd(S
, 5, 0);
2449 mhi
= lshift(mhi
, m2
);
2454 /* Compute mlo -- check for special case
2455 * that d is a normalized power of 2.
2460 mhi
= Balloc(mhi
->k
);
2464 mhi
= lshift(mhi
, Log2P
);
2470 dig
= quorem(b
,S
) + '0';
2471 /* Do we yet have the shortest decimal string
2472 * that will round to d?
2475 delta
= diff(S
, mhi
);
2478 j1
= delta
->sign
? 1 : cmp(b
, delta
);
2480 if (j1
== 0 && mode
!= 1 && !(word1(&u
) & 1)
2489 if (j
< 0 || (j
== 0 && mode
!= 1
2492 if (!b
->x
[0] && b
->wds
<= 1) {
2500 if ((j1
> 0 || (j1
== 0 && dig
& 1))
2509 if (dig
== '9') { /* possible if i == 1 */
2520 b
= multadd(b
, 10, 0);
2524 mlo
= mhi
= multadd(mhi
, 10, 0);
2529 mlo
= multadd(mlo
, 10, 0);
2532 mhi
= multadd(mhi
, 10, 0);
2540 *s
++ = dig
= quorem(b
,S
) + '0';
2541 if (!b
->x
[0] && b
->wds
<= 1) {
2546 b
= multadd(b
, 10, 0);
2551 /* Round off last digit */
2557 if (j
> 0 || (j
== 0 && dig
& 1)) {
2574 if (mlo
&& mlo
!= mhi
)
2588 if (mlo
&& mlo
!= mhi
)
2595 _Py_dg_freedtoa(s0
);
2602 #endif /* PY_NO_SHORT_FLOAT_REPR */