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. */
121 /* if PY_NO_SHORT_FLOAT_REPR is defined, then don't even try to compile
122 the following code */
123 #ifndef PY_NO_SHORT_FLOAT_REPR
127 #define MALLOC PyMem_Malloc
128 #define FREE PyMem_Free
130 /* This code should also work for ARM mixed-endian format on little-endian
131 machines, where doubles have byte order 45670123 (in increasing address
132 order, 0 being the least significant byte). */
133 #ifdef DOUBLE_IS_LITTLE_ENDIAN_IEEE754
136 #if defined(DOUBLE_IS_BIG_ENDIAN_IEEE754) || \
137 defined(DOUBLE_IS_ARM_MIXED_ENDIAN_IEEE754)
140 #if defined(IEEE_8087) + defined(IEEE_MC68k) != 1
141 #error "Exactly one of IEEE_8087 or IEEE_MC68k should be defined."
144 /* The code below assumes that the endianness of integers matches the
145 endianness of the two 32-bit words of a double. Check this. */
146 #if defined(WORDS_BIGENDIAN) && (defined(DOUBLE_IS_LITTLE_ENDIAN_IEEE754) || \
147 defined(DOUBLE_IS_ARM_MIXED_ENDIAN_IEEE754))
148 #error "doubles and ints have incompatible endianness"
151 #if !defined(WORDS_BIGENDIAN) && defined(DOUBLE_IS_BIG_ENDIAN_IEEE754)
152 #error "doubles and ints have incompatible endianness"
156 #if defined(HAVE_UINT32_T) && defined(HAVE_INT32_T)
157 typedef PY_UINT32_T ULong
;
158 typedef PY_INT32_T Long
;
160 #error "Failed to find an exact-width 32-bit integer type"
163 #if defined(HAVE_UINT64_T)
164 #define ULLong PY_UINT64_T
174 /* End Python #define linking */
177 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
181 #define PRIVATE_MEM 2304
183 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
184 static double private_mem
[PRIVATE_mem
], *pmem_next
= private_mem
;
190 typedef union { double d
; ULong L
[2]; } U
;
193 #define word0(x) (x)->L[1]
194 #define word1(x) (x)->L[0]
196 #define word0(x) (x)->L[0]
197 #define word1(x) (x)->L[1]
199 #define dval(x) (x)->d
201 #ifndef STRTOD_DIGLIM
202 #define STRTOD_DIGLIM 40
206 extern int strtod_diglim
;
208 #define strtod_diglim STRTOD_DIGLIM
211 /* The following definition of Storeinc is appropriate for MIPS processors.
212 * An alternative that might be better on some machines is
213 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
215 #if defined(IEEE_8087)
216 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
217 ((unsigned short *)a)[0] = (unsigned short)c, a++)
219 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
220 ((unsigned short *)a)[1] = (unsigned short)c, a++)
223 /* #define P DBL_MANT_DIG */
224 /* Ten_pmax = floor(P*log(2)/log(5)) */
225 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
226 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
227 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
230 #define Exp_shift1 20
231 #define Exp_msk1 0x100000
232 #define Exp_msk11 0x100000
233 #define Exp_mask 0x7ff00000
239 #define Exp_1 0x3ff00000
240 #define Exp_11 0x3ff00000
242 #define Frac_mask 0xfffff
243 #define Frac_mask1 0xfffff
246 #define Bndry_mask 0xfffff
247 #define Bndry_mask1 0xfffff
249 #define Sign_bit 0x80000000
258 #define Flt_Rounds FLT_ROUNDS
262 #endif /*Flt_Rounds*/
264 #define Rounding Flt_Rounds
266 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
267 #define Big1 0xffffffff
269 /* struct BCinfo is used to pass information from _Py_dg_strtod to bigcomp */
271 typedef struct BCinfo BCinfo
;
274 int dp0
, dp1
, dplen
, dsign
, e0
, inexact
;
275 int nd
, nd0
, rounding
, scale
, uflchk
;
278 #define FFFFFFFF 0xffffffffUL
282 /* struct Bigint is used to represent arbitrary-precision integers. These
283 integers are stored in sign-magnitude format, with the magnitude stored as
284 an array of base 2**32 digits. Bigints are always normalized: if x is a
285 Bigint then x->wds >= 1, and either x->wds == 1 or x[wds-1] is nonzero.
287 The Bigint fields are as follows:
289 - next is a header used by Balloc and Bfree to keep track of lists
290 of freed Bigints; it's also used for the linked list of
291 powers of 5 of the form 5**2**i used by pow5mult.
292 - k indicates which pool this Bigint was allocated from
293 - maxwds is the maximum number of words space was allocated for
294 (usually maxwds == 2**k)
295 - sign is 1 for negative Bigints, 0 for positive. The sign is unused
296 (ignored on inputs, set to 0 on outputs) in almost all operations
297 involving Bigints: a notable exception is the diff function, which
298 ignores signs on inputs but sets the sign of the output correctly.
299 - wds is the actual number of significant words
300 - x contains the vector of words (digits) for this Bigint, from least
301 significant (x[0]) to most significant (x[wds-1]).
307 int k
, maxwds
, sign
, wds
;
311 typedef struct Bigint Bigint
;
313 /* Memory management: memory is allocated from, and returned to, Kmax+1 pools
314 of memory, where pool k (0 <= k <= Kmax) is for Bigints b with b->maxwds ==
315 1 << k. These pools are maintained as linked lists, with freelist[k]
316 pointing to the head of the list for pool k.
318 On allocation, if there's no free slot in the appropriate pool, MALLOC is
319 called to get more memory. This memory is not returned to the system until
320 Python quits. There's also a private memory pool that's allocated from
321 in preference to using MALLOC.
323 For Bigints with more than (1 << Kmax) digits (which implies at least 1233
324 decimal digits), memory is directly allocated using MALLOC, and freed using
327 XXX: it would be easy to bypass this memory-management system and
328 translate each call to Balloc into a call to PyMem_Malloc, and each
329 Bfree to PyMem_Free. Investigate whether this has any significant
330 performance on impact. */
332 static Bigint
*freelist
[Kmax
+1];
334 /* Allocate space for a Bigint with up to 1<<k digits */
343 if (k
<= Kmax
&& (rv
= freelist
[k
]))
344 freelist
[k
] = rv
->next
;
347 len
= (sizeof(Bigint
) + (x
-1)*sizeof(ULong
) + sizeof(double) - 1)
349 if (k
<= Kmax
&& pmem_next
- private_mem
+ len
<= PRIVATE_mem
) {
350 rv
= (Bigint
*)pmem_next
;
354 rv
= (Bigint
*)MALLOC(len
*sizeof(double));
361 rv
->sign
= rv
->wds
= 0;
365 /* Free a Bigint allocated with Balloc */
374 v
->next
= freelist
[v
->k
];
380 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
381 y->wds*sizeof(Long) + 2*sizeof(int))
383 /* Multiply a Bigint b by m and add a. Either modifies b in place and returns
384 a pointer to the modified b, or Bfrees b and returns a pointer to a copy.
385 On failure, return NULL. In this case, b will have been already freed. */
388 multadd(Bigint
*b
, int m
, int a
) /* multiply by m and add a */
406 y
= *x
* (ULLong
)m
+ carry
;
408 *x
++ = (ULong
)(y
& FFFFFFFF
);
411 y
= (xi
& 0xffff) * m
+ carry
;
412 z
= (xi
>> 16) * m
+ (y
>> 16);
414 *x
++ = (z
<< 16) + (y
& 0xffff);
419 if (wds
>= b
->maxwds
) {
429 b
->x
[wds
++] = (ULong
)carry
;
435 /* convert a string s containing nd decimal digits (possibly containing a
436 decimal separator at position nd0, which is ignored) to a Bigint. This
437 function carries on where the parsing code in _Py_dg_strtod leaves off: on
438 entry, y9 contains the result of converting the first 9 digits. Returns
442 s2b(const char *s
, int nd0
, int nd
, ULong y9
, int dplen
)
449 for(k
= 0, y
= 1; x
> y
; y
<<= 1, k
++) ;
460 b
= multadd(b
, 10, *s
++ - '0');
469 b
= multadd(b
, 10, *s
++ - '0');
476 /* count leading 0 bits in the 32-bit integer x. */
483 if (!(x
& 0xffff0000)) {
487 if (!(x
& 0xff000000)) {
491 if (!(x
& 0xf0000000)) {
495 if (!(x
& 0xc0000000)) {
499 if (!(x
& 0x80000000)) {
501 if (!(x
& 0x40000000))
507 /* count trailing 0 bits in the 32-bit integer y, and shift y right by that
553 /* convert a small nonnegative integer to a Bigint */
568 /* multiply two Bigints. Returns a new Bigint, or NULL on failure. Ignores
569 the signs of a and b. */
572 mult(Bigint
*a
, Bigint
*b
)
576 ULong
*x
, *xa
, *xae
, *xb
, *xbe
, *xc
, *xc0
;
585 if (a
->wds
< b
->wds
) {
599 for(x
= c
->x
, xa
= x
+ wc
; x
< xa
; x
++)
607 for(; xb
< xbe
; xc0
++) {
613 z
= *x
++ * (ULLong
)y
+ *xc
+ carry
;
615 *xc
++ = (ULong
)(z
& FFFFFFFF
);
622 for(; xb
< xbe
; xb
++, xc0
++) {
623 if (y
= *xb
& 0xffff) {
628 z
= (*x
& 0xffff) * y
+ (*xc
& 0xffff) + carry
;
630 z2
= (*x
++ >> 16) * y
+ (*xc
>> 16) + carry
;
643 z
= (*x
& 0xffff) * y
+ (*xc
>> 16) + carry
;
646 z2
= (*x
++ >> 16) * y
+ (*xc
& 0xffff) + carry
;
654 for(xc0
= c
->x
, xc
= xc0
+ wc
; wc
> 0 && !*--xc
; --wc
) ;
659 /* p5s is a linked list of powers of 5 of the form 5**(2**i), i >= 2 */
663 /* multiply the Bigint b by 5**k. Returns a pointer to the result, or NULL on
664 failure; if the returned pointer is distinct from b then the original
665 Bigint b will have been Bfree'd. Ignores the sign of b. */
668 pow5mult(Bigint
*b
, int k
)
670 Bigint
*b1
, *p5
, *p51
;
672 static int p05
[3] = { 5, 25, 125 };
675 b
= multadd(b
, p05
[i
-1], 0);
718 /* shift a Bigint b left by k bits. Return a pointer to the shifted result,
719 or NULL on failure. If the returned pointer is distinct from b then the
720 original b will have been Bfree'd. Ignores the sign of b. */
723 lshift(Bigint
*b
, int k
)
727 ULong
*x
, *x1
, *xe
, z
;
732 for(i
= b
->maxwds
; n1
> i
; i
<<= 1)
740 for(i
= 0; i
< n
; i
++)
763 /* Do a three-way compare of a and b, returning -1 if a < b, 0 if a == b and
764 1 if a > b. Ignores signs of a and b. */
767 cmp(Bigint
*a
, Bigint
*b
)
769 ULong
*xa
, *xa0
, *xb
, *xb0
;
775 if (i
> 1 && !a
->x
[i
-1])
776 Bug("cmp called with a->x[a->wds-1] == 0");
777 if (j
> 1 && !b
->x
[j
-1])
778 Bug("cmp called with b->x[b->wds-1] == 0");
788 return *xa
< *xb
? -1 : 1;
795 /* Take the difference of Bigints a and b, returning a new Bigint. Returns
796 NULL on failure. The signs of a and b are ignored, but the sign of the
797 result is set appropriately. */
800 diff(Bigint
*a
, Bigint
*b
)
804 ULong
*xa
, *xae
, *xb
, *xbe
, *xc
;
843 y
= (ULLong
)*xa
++ - *xb
++ - borrow
;
844 borrow
= y
>> 32 & (ULong
)1;
845 *xc
++ = (ULong
)(y
& FFFFFFFF
);
850 borrow
= y
>> 32 & (ULong
)1;
851 *xc
++ = (ULong
)(y
& FFFFFFFF
);
855 y
= (*xa
& 0xffff) - (*xb
& 0xffff) - borrow
;
856 borrow
= (y
& 0x10000) >> 16;
857 z
= (*xa
++ >> 16) - (*xb
++ >> 16) - borrow
;
858 borrow
= (z
& 0x10000) >> 16;
863 y
= (*xa
& 0xffff) - borrow
;
864 borrow
= (y
& 0x10000) >> 16;
865 z
= (*xa
++ >> 16) - borrow
;
866 borrow
= (z
& 0x10000) >> 16;
876 /* Given a positive normal double x, return the difference between x and the next
877 double up. Doesn't give correct results for subnormals. */
885 L
= (word0(x
) & Exp_mask
) - (P
-1)*Exp_msk1
;
891 /* Convert a Bigint to a double plus an exponent */
894 b2d(Bigint
*a
, int *e
)
896 ULong
*xa
, *xa0
, w
, y
, z
;
904 if (!y
) Bug("zero y in b2d");
909 word0(&d
) = Exp_1
| y
>> (Ebits
- k
);
910 w
= xa
> xa0
? *--xa
: 0;
911 word1(&d
) = y
<< ((32-Ebits
) + k
) | w
>> (Ebits
- k
);
914 z
= xa
> xa0
? *--xa
: 0;
916 word0(&d
) = Exp_1
| y
<< k
| z
>> (32 - k
);
917 y
= xa
> xa0
? *--xa
: 0;
918 word1(&d
) = z
<< k
| y
>> (32 - k
);
921 word0(&d
) = Exp_1
| y
;
928 /* Convert a double to a Bigint plus an exponent. Return NULL on failure.
930 Given a finite nonzero double d, return an odd Bigint b and exponent *e
931 such that fabs(d) = b * 2**e. On return, *bbits gives the number of
932 significant bits of e; that is, 2**(*bbits-1) <= b < 2**(*bbits).
934 If d is zero, then b == 0, *e == -1010, *bbits = 0.
939 d2b(U
*d
, int *e
, int *bits
)
951 z
= word0(d
) & Frac_mask
;
952 word0(d
) &= 0x7fffffff; /* clear sign bit, which we ignore */
953 if ((de
= (int)(word0(d
) >> Exp_shift
)))
955 if ((y
= word1(d
))) {
956 if ((k
= lo0bits(&y
))) {
957 x
[0] = y
| z
<< (32 - k
);
963 b
->wds
= (x
[1] = z
) ? 2 : 1;
973 *e
= de
- Bias
- (P
-1) + k
;
977 *e
= de
- Bias
- (P
-1) + 1 + k
;
978 *bits
= 32*i
- hi0bits(x
[i
-1]);
983 /* Compute the ratio of two Bigints, as a double. The result may have an
984 error of up to 2.5 ulps. */
987 ratio(Bigint
*a
, Bigint
*b
)
992 dval(&da
) = b2d(a
, &ka
);
993 dval(&db
) = b2d(b
, &kb
);
994 k
= ka
- kb
+ 32*(a
->wds
- b
->wds
);
996 word0(&da
) += k
*Exp_msk1
;
999 word0(&db
) += k
*Exp_msk1
;
1001 return dval(&da
) / dval(&db
);
1006 1e0
, 1e1
, 1e2
, 1e3
, 1e4
, 1e5
, 1e6
, 1e7
, 1e8
, 1e9
,
1007 1e10
, 1e11
, 1e12
, 1e13
, 1e14
, 1e15
, 1e16
, 1e17
, 1e18
, 1e19
,
1012 bigtens
[] = { 1e16
, 1e32
, 1e64
, 1e128
, 1e256
};
1013 static const double tinytens
[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1014 9007199254740992.*9007199254740992.e
-256
1015 /* = 2^106 * 1e-256 */
1017 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1018 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1019 #define Scale_Bit 0x10
1028 dshift(Bigint
*b
, int p2
)
1030 int rv
= hi0bits(b
->x
[b
->wds
-1]) - 4;
1036 /* special case of Bigint division. The quotient is always in the range 0 <=
1037 quotient < 10, and on entry the divisor S is normalized so that its top 4
1038 bits (28--31) are zero and bit 27 is set. */
1041 quorem(Bigint
*b
, Bigint
*S
)
1044 ULong
*bx
, *bxe
, q
, *sx
, *sxe
;
1046 ULLong borrow
, carry
, y
, ys
;
1048 ULong borrow
, carry
, y
, ys
;
1054 /*debug*/ if (b
->wds
> n
)
1055 /*debug*/ Bug("oversize b in quorem");
1063 q
= *bxe
/ (*sxe
+ 1); /* ensure q <= true quotient */
1065 /*debug*/ if (q
> 9)
1066 /*debug*/ Bug("oversized quotient in quorem");
1073 ys
= *sx
++ * (ULLong
)q
+ carry
;
1075 y
= *bx
- (ys
& FFFFFFFF
) - borrow
;
1076 borrow
= y
>> 32 & (ULong
)1;
1077 *bx
++ = (ULong
)(y
& FFFFFFFF
);
1080 ys
= (si
& 0xffff) * q
+ carry
;
1081 zs
= (si
>> 16) * q
+ (ys
>> 16);
1083 y
= (*bx
& 0xffff) - (ys
& 0xffff) - borrow
;
1084 borrow
= (y
& 0x10000) >> 16;
1085 z
= (*bx
>> 16) - (zs
& 0xffff) - borrow
;
1086 borrow
= (z
& 0x10000) >> 16;
1093 while(--bxe
> bx
&& !*bxe
)
1098 if (cmp(b
, S
) >= 0) {
1108 y
= *bx
- (ys
& FFFFFFFF
) - borrow
;
1109 borrow
= y
>> 32 & (ULong
)1;
1110 *bx
++ = (ULong
)(y
& FFFFFFFF
);
1113 ys
= (si
& 0xffff) + carry
;
1114 zs
= (si
>> 16) + (ys
>> 16);
1116 y
= (*bx
& 0xffff) - (ys
& 0xffff) - borrow
;
1117 borrow
= (y
& 0x10000) >> 16;
1118 z
= (*bx
>> 16) - (zs
& 0xffff) - borrow
;
1119 borrow
= (z
& 0x10000) >> 16;
1127 while(--bxe
> bx
&& !*bxe
)
1136 /* return 0 on success, -1 on failure */
1139 bigcomp(U
*rv
, const char *s0
, BCinfo
*bc
)
1142 int b2
, bbits
, d2
, dd
, dig
, dsign
, i
, j
, nd
, nd0
, p2
, p5
, speccase
;
1147 p5
= nd
+ bc
->e0
- 1;
1149 if (rv
->d
== 0.) { /* special case: value near underflow-to-zero */
1150 /* threshold was rounded to zero */
1156 word0(rv
) = (P
+2) << Exp_shift
;
1167 b
= d2b(rv
, &p2
, &bbits
);
1172 /* floor(log2(rv)) == bbits - 1 + p2 */
1173 /* Check for denormal case. */
1175 if (i
> (j
= P
- Emin
- 1 + p2
)) {
1191 /* Arrange for convenient computation of quotients:
1192 * shift left if necessary so divisor has 4 leading 0 bits.
1195 d
= pow5mult(d
, p5
);
1202 b
= pow5mult(b
, -p5
);
1217 if ((b2
+= i
) > 0) {
1224 if ((d2
+= i
) > 0) {
1232 /* Now b/d = exactly half-way between the two floating-point values */
1233 /* on either side of the input string. Compute first digit of b/d. */
1235 if (!(dig
= quorem(b
,d
))) {
1236 b
= multadd(b
, 10, 0); /* very unlikely */
1244 /* Compare b/d with s0 */
1247 dd
= 9999; /* silence gcc compiler warning */
1248 for(i
= 0; i
< nd0
; ) {
1249 if ((dd
= s0
[i
++] - '0' - dig
))
1251 if (!b
->x
[0] && b
->wds
== 1) {
1256 b
= multadd(b
, 10, 0);
1263 for(j
= bc
->dp1
; i
++ < nd
;) {
1264 if ((dd
= s0
[j
++] - '0' - dig
))
1266 if (!b
->x
[0] && b
->wds
== 1) {
1271 b
= multadd(b
, 10, 0);
1278 if (b
->x
[0] || b
->wds
> 1)
1288 if (!dsign
) /* does not happen for round-near */
1290 dval(rv
) -= ulp(rv
);
1295 dval(rv
) += ulp(rv
);
1299 /* Exact half-way case: apply round-even rule. */
1300 if (word1(rv
) & 1) {
1311 _Py_dg_strtod(const char *s00
, char **se
)
1313 int bb2
, bb5
, bbe
, bd2
, bd5
, bbbits
, bs2
, c
, e
, e1
, error
;
1314 int esign
, i
, j
, k
, nd
, nd0
, nf
, nz
, nz0
, sign
;
1315 const char *s
, *s0
, *s1
;
1318 U aadj2
, adj
, rv
, rv0
;
1321 Bigint
*bb
, *bb1
, *bd
, *bd0
, *bs
, *delta
;
1323 sign
= nz0
= nz
= bc
.dplen
= bc
.uflchk
= 0;
1325 for(s
= s00
;;s
++) switch(*s
) {
1335 /* modify original dtoa.c so that it doesn't accept leading whitespace
1350 while(*++s
== '0') ;
1356 for(nd
= nf
= 0; (c
= *s
) >= '0' && c
<= '9'; nd
++, s
++)
1362 bc
.dp0
= bc
.dp1
= s
- s0
;
1366 bc
.dplen
= bc
.dp1
- bc
.dp0
;
1368 for(; c
== '0'; c
= *++s
)
1370 if (c
> '0' && c
<= '9') {
1378 for(; c
>= '0' && c
<= '9'; c
= *++s
) {
1383 for(i
= 1; i
< nz
; i
++)
1386 else if (nd
<= DBL_DIG
+ 1)
1390 else if (nd
<= DBL_DIG
+ 1)
1398 if (c
== 'e' || c
== 'E') {
1399 if (!nd
&& !nz
&& !nz0
) {
1410 if (c
>= '0' && c
<= '9') {
1413 if (c
> '0' && c
<= '9') {
1416 while((c
= *++s
) >= '0' && c
<= '9')
1418 if (s
- s1
> 8 || L
> 19999)
1419 /* Avoid confusion from exponents
1420 * so large that e might overflow.
1422 e
= 19999; /* safe for 16 bit ints */
1442 bc
.e0
= e1
= e
-= nf
;
1444 /* Now we have nd0 digits, starting at s0, followed by a
1445 * decimal point, followed by nd-nd0 digits. The number we're
1446 * after is the integer represented by those digits times
1451 k
= nd
< DBL_DIG
+ 1 ? nd
: DBL_DIG
+ 1;
1454 dval(&rv
) = tens
[k
- 9] * dval(&rv
) + z
;
1463 if (e
<= Ten_pmax
) {
1464 dval(&rv
) *= tens
[e
];
1468 if (e
<= Ten_pmax
+ i
) {
1469 /* A fancier test would sometimes let us do
1470 * this for larger i values.
1473 dval(&rv
) *= tens
[i
];
1474 dval(&rv
) *= tens
[e
];
1478 else if (e
>= -Ten_pmax
) {
1479 dval(&rv
) /= tens
[-e
];
1487 /* Get starting approximation = rv * 10**e1 */
1491 dval(&rv
) *= tens
[i
];
1493 if (e1
> DBL_MAX_10_EXP
) {
1496 /* Can't trust HUGE_VAL */
1497 word0(&rv
) = Exp_mask
;
1502 for(j
= 0; e1
> 1; j
++, e1
>>= 1)
1504 dval(&rv
) *= bigtens
[j
];
1505 /* The last multiplication could overflow. */
1506 word0(&rv
) -= P
*Exp_msk1
;
1507 dval(&rv
) *= bigtens
[j
];
1508 if ((z
= word0(&rv
) & Exp_mask
)
1509 > Exp_msk1
*(DBL_MAX_EXP
+Bias
-P
))
1511 if (z
> Exp_msk1
*(DBL_MAX_EXP
+Bias
-1-P
)) {
1512 /* set to largest number */
1513 /* (Can't trust DBL_MAX) */
1518 word0(&rv
) += P
*Exp_msk1
;
1524 dval(&rv
) /= tens
[i
];
1526 if (e1
>= 1 << n_bigtens
)
1530 for(j
= 0; e1
> 0; j
++, e1
>>= 1)
1532 dval(&rv
) *= tinytens
[j
];
1533 if (bc
.scale
&& (j
= 2*P
+ 1 - ((word0(&rv
) & Exp_mask
)
1534 >> Exp_shift
)) > 0) {
1535 /* scaled rv is denormal; clear j low bits */
1539 word0(&rv
) = (P
+2)*Exp_msk1
;
1541 word0(&rv
) &= 0xffffffff << (j
-32);
1544 word1(&rv
) &= 0xffffffff << j
;
1555 /* Now the hard part -- adjusting rv to the correct value.*/
1557 /* Put digits into bd: true value = bd * 10^e */
1560 bc
.nd0
= nd0
; /* Only needed if nd > strtod_diglim, but done here */
1561 /* to silence an erroneous warning about bc.nd0 */
1562 /* possibly not being initialized. */
1563 if (nd
> strtod_diglim
) {
1564 /* ASSERT(strtod_diglim >= 18); 18 == one more than the */
1565 /* minimum number of decimal digits to distinguish double values */
1566 /* in IEEE arithmetic. */
1571 if (--j
<= bc
.dp1
&& j
>= bc
.dp0
)
1581 if (nd
< 9) { /* must recompute y */
1583 for(i
= 0; i
< nd0
; ++i
)
1584 y
= 10*y
+ s0
[i
] - '0';
1585 for(j
= bc
.dp1
; i
< nd
; ++i
)
1586 y
= 10*y
+ s0
[j
++] - '0';
1589 bd0
= s2b(s0
, nd0
, nd
, y
, bc
.dplen
);
1594 bd
= Balloc(bd0
->k
);
1600 bb
= d2b(&rv
, &bbe
, &bbbits
); /* rv = bb * 2^bbe */
1628 i
= j
+ bbbits
- 1; /* logb(rv) */
1629 if (i
< Emin
) /* denormal */
1636 i
= bb2
< bd2
? bb2
: bd2
;
1645 bs
= pow5mult(bs
, bb5
);
1663 bb
= lshift(bb
, bb2
);
1672 bd
= pow5mult(bd
, bd5
);
1681 bd
= lshift(bd
, bd2
);
1690 bs
= lshift(bs
, bs2
);
1698 delta
= diff(bb
, bd
);
1699 if (delta
== NULL
) {
1706 bc
.dsign
= delta
->sign
;
1709 if (bc
.nd
> nd
&& i
<= 0) {
1711 break; /* Must use bigcomp(). */
1714 i
= -1; /* Discarded digits make delta smaller. */
1719 /* Error is less than half an ulp -- check for
1720 * special case of mantissa a power of two.
1722 if (bc
.dsign
|| word1(&rv
) || word0(&rv
) & Bndry_mask
1723 || (word0(&rv
) & Exp_mask
) <= (2*P
+1)*Exp_msk1
1727 if (!delta
->x
[0] && delta
->wds
<= 1) {
1731 delta
= lshift(delta
,Log2P
);
1732 if (delta
== NULL
) {
1739 if (cmp(delta
, bs
) > 0)
1744 /* exactly half-way between */
1746 if ((word0(&rv
) & Bndry_mask1
) == Bndry_mask1
1749 (y
= word0(&rv
) & Exp_mask
) <= 2*P
*Exp_msk1
) ?
1750 (0xffffffff & (0xffffffff << (2*P
+1-(y
>>Exp_shift
)))) :
1752 /*boundary case -- increment exponent*/
1753 word0(&rv
) = (word0(&rv
) & Exp_mask
)
1761 else if (!(word0(&rv
) & Bndry_mask
) && !word1(&rv
)) {
1763 /* boundary case -- decrement exponent */
1765 L
= word0(&rv
) & Exp_mask
;
1766 if (L
<= (2*P
+1)*Exp_msk1
) {
1767 if (L
> (P
+2)*Exp_msk1
)
1768 /* round even ==> */
1771 /* rv = smallest denormal */
1779 L
= (word0(&rv
) & Exp_mask
) - Exp_msk1
;
1780 word0(&rv
) = L
| Bndry_mask1
;
1781 word1(&rv
) = 0xffffffff;
1784 if (!(word1(&rv
) & LSB
))
1787 dval(&rv
) += ulp(&rv
);
1789 dval(&rv
) -= ulp(&rv
);
1798 bc
.dsign
= 1 - bc
.dsign
;
1801 if ((aadj
= ratio(delta
, bs
)) <= 2.) {
1804 else if (word1(&rv
) || word0(&rv
) & Bndry_mask
) {
1805 if (word1(&rv
) == Tiny1
&& !word0(&rv
)) {
1816 /* special case -- power of FLT_RADIX to be */
1817 /* rounded down... */
1819 if (aadj
< 2./FLT_RADIX
)
1820 aadj
= 1./FLT_RADIX
;
1828 aadj1
= bc
.dsign
? aadj
: -aadj
;
1829 if (Flt_Rounds
== 0)
1832 y
= word0(&rv
) & Exp_mask
;
1834 /* Check for overflow */
1836 if (y
== Exp_msk1
*(DBL_MAX_EXP
+Bias
-1)) {
1837 dval(&rv0
) = dval(&rv
);
1838 word0(&rv
) -= P
*Exp_msk1
;
1839 adj
.d
= aadj1
* ulp(&rv
);
1841 if ((word0(&rv
) & Exp_mask
) >=
1842 Exp_msk1
*(DBL_MAX_EXP
+Bias
-P
)) {
1843 if (word0(&rv0
) == Big0
&& word1(&rv0
) == Big1
)
1850 word0(&rv
) += P
*Exp_msk1
;
1853 if (bc
.scale
&& y
<= 2*P
*Exp_msk1
) {
1854 if (aadj
<= 0x7fffffff) {
1855 if ((z
= (ULong
)aadj
) <= 0)
1858 aadj1
= bc
.dsign
? aadj
: -aadj
;
1860 dval(&aadj2
) = aadj1
;
1861 word0(&aadj2
) += (2*P
+1)*Exp_msk1
- y
;
1862 aadj1
= dval(&aadj2
);
1864 adj
.d
= aadj1
* ulp(&rv
);
1867 z
= word0(&rv
) & Exp_mask
;
1871 /* Can we stop now? */
1874 /* The tolerances below are conservative. */
1875 if (bc
.dsign
|| word1(&rv
) || word0(&rv
) & Bndry_mask
) {
1876 if (aadj
< .4999999 || aadj
> .5000001)
1879 else if (aadj
< .4999999/FLT_RADIX
)
1895 error
= bigcomp(&rv
, s0
, &bc
);
1901 word0(&rv0
) = Exp_1
- 2*P
*Exp_msk1
;
1903 dval(&rv
) *= dval(&rv0
);
1904 /* try to avoid the bug of testing an 8087 register value */
1905 if (!(word0(&rv
) & Exp_mask
))
1911 return sign
? -dval(&rv
) : dval(&rv
);
1927 sizeof(Bigint
) - sizeof(ULong
) - sizeof(int) + j
<= (unsigned)i
;
1930 r
= (int*)Balloc(k
);
1934 return (char *)(r
+1);
1938 nrv_alloc(char *s
, char **rve
, int n
)
1946 while((*t
= *s
++)) t
++;
1952 /* freedtoa(s) must be used to free values s returned by dtoa
1953 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
1954 * but for consistency with earlier versions of dtoa, it is optional
1955 * when MULTIPLE_THREADS is not defined.
1959 _Py_dg_freedtoa(char *s
)
1961 Bigint
*b
= (Bigint
*)((int *)s
- 1);
1962 b
->maxwds
= 1 << (b
->k
= *(int*)b
);
1966 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1968 * Inspired by "How to Print Floating-Point Numbers Accurately" by
1969 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
1972 * 1. Rather than iterating, we use a simple numeric overestimate
1973 * to determine k = floor(log10(d)). We scale relevant
1974 * quantities using O(log2(k)) rather than O(k) multiplications.
1975 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1976 * try to generate digits strictly left to right. Instead, we
1977 * compute with fewer bits and propagate the carry if necessary
1978 * when rounding the final digit up. This is often faster.
1979 * 3. Under the assumption that input will be rounded nearest,
1980 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1981 * That is, we allow equality in stopping tests when the
1982 * round-nearest rule will give the same floating-point value
1983 * as would satisfaction of the stopping test with strict
1985 * 4. We remove common factors of powers of 2 from relevant
1987 * 5. When converting floating-point integers less than 1e16,
1988 * we use floating-point arithmetic rather than resorting
1989 * to multiple-precision integers.
1990 * 6. When asked to produce fewer than 15 digits, we first try
1991 * to get by with floating-point arithmetic; we resort to
1992 * multiple-precision integer arithmetic only if we cannot
1993 * guarantee that the floating-point calculation has given
1994 * the correctly rounded result. For k requested digits and
1995 * "uniformly" distributed input, the probability is
1996 * something like 10^(k-15) that we must resort to the Long
2000 /* Additional notes (METD): (1) returns NULL on failure. (2) to avoid memory
2001 leakage, a successful call to _Py_dg_dtoa should always be matched by a
2002 call to _Py_dg_freedtoa. */
2005 _Py_dg_dtoa(double dd
, int mode
, int ndigits
,
2006 int *decpt
, int *sign
, char **rve
)
2008 /* Arguments ndigits, decpt, sign are similar to those
2009 of ecvt and fcvt; trailing zeros are suppressed from
2010 the returned string. If not null, *rve is set to point
2011 to the end of the return value. If d is +-Infinity or NaN,
2012 then *decpt is set to 9999.
2015 0 ==> shortest string that yields d when read in
2016 and rounded to nearest.
2017 1 ==> like 0, but with Steele & White stopping rule;
2018 e.g. with IEEE P754 arithmetic , mode 0 gives
2019 1e23 whereas mode 1 gives 9.999999999999999e22.
2020 2 ==> max(1,ndigits) significant digits. This gives a
2021 return value similar to that of ecvt, except
2022 that trailing zeros are suppressed.
2023 3 ==> through ndigits past the decimal point. This
2024 gives a return value similar to that from fcvt,
2025 except that trailing zeros are suppressed, and
2026 ndigits can be negative.
2027 4,5 ==> similar to 2 and 3, respectively, but (in
2028 round-nearest mode) with the tests of mode 0 to
2029 possibly return a shorter string that rounds to d.
2030 With IEEE arithmetic and compilation with
2031 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
2032 as modes 2 and 3 when FLT_ROUNDS != 1.
2033 6-9 ==> Debugging modes similar to mode - 4: don't try
2034 fast floating-point estimate (if applicable).
2036 Values of mode other than 0-9 are treated as mode 0.
2038 Sufficient space is allocated to the return value
2039 to hold the suppressed trailing zeros.
2042 int bbits
, b2
, b5
, be
, dig
, i
, ieps
, ilim
, ilim0
, ilim1
,
2043 j
, j1
, k
, k0
, k_check
, leftright
, m2
, m5
, s2
, s5
,
2044 spec_case
, try_quick
;
2048 Bigint
*b
, *b1
, *delta
, *mlo
, *mhi
, *S
;
2053 /* set pointers to NULL, to silence gcc compiler warnings and make
2054 cleanup easier on error */
2055 mlo
= mhi
= b
= S
= 0;
2059 if (word0(&u
) & Sign_bit
) {
2060 /* set sign for everything, including 0's and NaNs */
2062 word0(&u
) &= ~Sign_bit
; /* clear sign bit */
2067 /* quick return for Infinities, NaNs and zeros */
2068 if ((word0(&u
) & Exp_mask
) == Exp_mask
)
2070 /* Infinity or NaN */
2072 if (!word1(&u
) && !(word0(&u
) & 0xfffff))
2073 return nrv_alloc("Infinity", rve
, 8);
2074 return nrv_alloc("NaN", rve
, 3);
2078 return nrv_alloc("0", rve
, 1);
2081 /* compute k = floor(log10(d)). The computation may leave k
2082 one too large, but should never leave k too small. */
2083 b
= d2b(&u
, &be
, &bbits
);
2086 if ((i
= (int)(word0(&u
) >> Exp_shift1
& (Exp_mask
>>Exp_shift1
)))) {
2087 dval(&d2
) = dval(&u
);
2088 word0(&d2
) &= Frac_mask1
;
2089 word0(&d2
) |= Exp_11
;
2091 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2092 * log10(x) = log(x) / log(10)
2093 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2094 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2096 * This suggests computing an approximation k to log10(d) by
2098 * k = (i - Bias)*0.301029995663981
2099 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2101 * We want k to be too large rather than too small.
2102 * The error in the first-order Taylor series approximation
2103 * is in our favor, so we just round up the constant enough
2104 * to compensate for any error in the multiplication of
2105 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2106 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2107 * adding 1e-13 to the constant term more than suffices.
2108 * Hence we adjust the constant term to 0.1760912590558.
2109 * (We could get a more accurate k by invoking log10,
2110 * but this is probably not worthwhile.)
2117 /* d is denormalized */
2119 i
= bbits
+ be
+ (Bias
+ (P
-1) - 1);
2120 x
= i
> 32 ? word0(&u
) << (64 - i
) | word1(&u
) >> (i
- 32)
2121 : word1(&u
) << (32 - i
);
2123 word0(&d2
) -= 31*Exp_msk1
; /* adjust exponent */
2124 i
-= (Bias
+ (P
-1) - 1) + 1;
2127 ds
= (dval(&d2
)-1.5)*0.289529654602168 + 0.1760912590558 +
2128 i
*0.301029995663981;
2130 if (ds
< 0. && ds
!= k
)
2131 k
--; /* want k = floor(ds) */
2133 if (k
>= 0 && k
<= Ten_pmax
) {
2134 if (dval(&u
) < tens
[k
])
2157 if (mode
< 0 || mode
> 9)
2167 ilim
= ilim1
= -1; /* Values for cases 0 and 1; done here to */
2168 /* silence erroneous "gcc -Wall" warning. */
2181 ilim
= ilim1
= i
= ndigits
;
2187 i
= ndigits
+ k
+ 1;
2199 if (ilim
>= 0 && ilim
<= Quick_max
&& try_quick
) {
2201 /* Try to get by with floating-point arithmetic. */
2204 dval(&d2
) = dval(&u
);
2207 ieps
= 2; /* conservative */
2212 /* prevent overflows */
2214 dval(&u
) /= bigtens
[n_bigtens
-1];
2217 for(; j
; j
>>= 1, i
++)
2224 else if ((j1
= -k
)) {
2225 dval(&u
) *= tens
[j1
& 0xf];
2226 for(j
= j1
>> 4; j
; j
>>= 1, i
++)
2229 dval(&u
) *= bigtens
[i
];
2232 if (k_check
&& dval(&u
) < 1. && ilim
> 0) {
2240 dval(&eps
) = ieps
*dval(&u
) + 7.;
2241 word0(&eps
) -= (P
-1)*Exp_msk1
;
2245 if (dval(&u
) > dval(&eps
))
2247 if (dval(&u
) < -dval(&eps
))
2252 /* Use Steele & White method of only
2253 * generating digits needed.
2255 dval(&eps
) = 0.5/tens
[ilim
-1] - dval(&eps
);
2259 *s
++ = '0' + (int)L
;
2260 if (dval(&u
) < dval(&eps
))
2262 if (1. - dval(&u
) < dval(&eps
))
2271 /* Generate ilim digits, then fix them up. */
2272 dval(&eps
) *= tens
[ilim
-1];
2273 for(i
= 1;; i
++, dval(&u
) *= 10.) {
2274 L
= (Long
)(dval(&u
));
2275 if (!(dval(&u
) -= L
))
2277 *s
++ = '0' + (int)L
;
2279 if (dval(&u
) > 0.5 + dval(&eps
))
2281 else if (dval(&u
) < 0.5 - dval(&eps
)) {
2292 dval(&u
) = dval(&d2
);
2297 /* Do we have a "small" integer? */
2299 if (be
>= 0 && k
<= Int_max
) {
2302 if (ndigits
< 0 && ilim
<= 0) {
2304 if (ilim
< 0 || dval(&u
) <= 5*ds
)
2308 for(i
= 1;; i
++, dval(&u
) *= 10.) {
2309 L
= (Long
)(dval(&u
) / ds
);
2311 *s
++ = '0' + (int)L
;
2316 dval(&u
) += dval(&u
);
2317 if (dval(&u
) > ds
|| (dval(&u
) == ds
&& L
& 1)) {
2337 denorm
? be
+ (Bias
+ (P
-1) - 1 + 1) :
2345 if (m2
> 0 && s2
> 0) {
2346 i
= m2
< s2
? m2
: s2
;
2354 mhi
= pow5mult(mhi
, m5
);
2363 if ((j
= b5
- m5
)) {
2370 b
= pow5mult(b
, b5
);
2379 S
= pow5mult(S
, s5
);
2384 /* Check for special case that d is a normalized power of 2. */
2387 if ((mode
< 2 || leftright
)
2389 if (!word1(&u
) && !(word0(&u
) & Bndry_mask
)
2390 && word0(&u
) & (Exp_mask
& ~Exp_msk1
)
2392 /* The special case */
2399 /* Arrange for convenient computation of quotients:
2400 * shift left if necessary so divisor has 4 leading 0 bits.
2402 * Perhaps we should just compute leading 28 bits of S once
2403 * and for all and pass them and a shift to quorem, so it
2404 * can do shifts and ors to compute the numerator for q.
2406 if ((i
= ((s5
? 32 - hi0bits(S
->x
[S
->wds
-1]) : 1) + s2
) & 0x1f))
2426 b
= multadd(b
, 10, 0); /* we botched the k estimate */
2430 mhi
= multadd(mhi
, 10, 0);
2437 if (ilim
<= 0 && (mode
== 3 || mode
== 5)) {
2439 /* no digits, fcvt style */
2445 S
= multadd(S
, 5, 0);
2458 mhi
= lshift(mhi
, m2
);
2463 /* Compute mlo -- check for special case
2464 * that d is a normalized power of 2.
2469 mhi
= Balloc(mhi
->k
);
2473 mhi
= lshift(mhi
, Log2P
);
2479 dig
= quorem(b
,S
) + '0';
2480 /* Do we yet have the shortest decimal string
2481 * that will round to d?
2484 delta
= diff(S
, mhi
);
2487 j1
= delta
->sign
? 1 : cmp(b
, delta
);
2489 if (j1
== 0 && mode
!= 1 && !(word1(&u
) & 1)
2498 if (j
< 0 || (j
== 0 && mode
!= 1
2501 if (!b
->x
[0] && b
->wds
<= 1) {
2509 if ((j1
> 0 || (j1
== 0 && dig
& 1))
2518 if (dig
== '9') { /* possible if i == 1 */
2529 b
= multadd(b
, 10, 0);
2533 mlo
= mhi
= multadd(mhi
, 10, 0);
2538 mlo
= multadd(mlo
, 10, 0);
2541 mhi
= multadd(mhi
, 10, 0);
2549 *s
++ = dig
= quorem(b
,S
) + '0';
2550 if (!b
->x
[0] && b
->wds
<= 1) {
2555 b
= multadd(b
, 10, 0);
2560 /* Round off last digit */
2566 if (j
> 0 || (j
== 0 && dig
& 1)) {
2583 if (mlo
&& mlo
!= mhi
)
2597 if (mlo
&& mlo
!= mhi
)
2604 _Py_dg_freedtoa(s0
);
2611 #endif /* PY_NO_SHORT_FLOAT_REPR */