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
204 extern int strtod_diglim
;
206 #define strtod_diglim STRTOD_DIGLIM
209 /* The following definition of Storeinc is appropriate for MIPS processors.
210 * An alternative that might be better on some machines is
211 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
213 #if defined(IEEE_8087)
214 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
215 ((unsigned short *)a)[0] = (unsigned short)c, a++)
217 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
218 ((unsigned short *)a)[1] = (unsigned short)c, a++)
221 /* #define P DBL_MANT_DIG */
222 /* Ten_pmax = floor(P*log(2)/log(5)) */
223 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
224 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
225 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
228 #define Exp_shift1 20
229 #define Exp_msk1 0x100000
230 #define Exp_msk11 0x100000
231 #define Exp_mask 0x7ff00000
237 #define Exp_1 0x3ff00000
238 #define Exp_11 0x3ff00000
240 #define Frac_mask 0xfffff
241 #define Frac_mask1 0xfffff
244 #define Bndry_mask 0xfffff
245 #define Bndry_mask1 0xfffff
247 #define Sign_bit 0x80000000
256 #define Flt_Rounds FLT_ROUNDS
260 #endif /*Flt_Rounds*/
262 #define Rounding Flt_Rounds
264 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
265 #define Big1 0xffffffff
268 #define NAN_WORD0 0x7ff80000
276 /* struct BCinfo is used to pass information from _Py_dg_strtod to bigcomp */
278 typedef struct BCinfo BCinfo
;
281 int dp0
, dp1
, dplen
, dsign
, e0
, inexact
;
282 int nd
, nd0
, rounding
, scale
, uflchk
;
285 #define FFFFFFFF 0xffffffffUL
289 /* struct Bigint is used to represent arbitrary-precision integers. These
290 integers are stored in sign-magnitude format, with the magnitude stored as
291 an array of base 2**32 digits. Bigints are always normalized: if x is a
292 Bigint then x->wds >= 1, and either x->wds == 1 or x[wds-1] is nonzero.
294 The Bigint fields are as follows:
296 - next is a header used by Balloc and Bfree to keep track of lists
297 of freed Bigints; it's also used for the linked list of
298 powers of 5 of the form 5**2**i used by pow5mult.
299 - k indicates which pool this Bigint was allocated from
300 - maxwds is the maximum number of words space was allocated for
301 (usually maxwds == 2**k)
302 - sign is 1 for negative Bigints, 0 for positive. The sign is unused
303 (ignored on inputs, set to 0 on outputs) in almost all operations
304 involving Bigints: a notable exception is the diff function, which
305 ignores signs on inputs but sets the sign of the output correctly.
306 - wds is the actual number of significant words
307 - x contains the vector of words (digits) for this Bigint, from least
308 significant (x[0]) to most significant (x[wds-1]).
314 int k
, maxwds
, sign
, wds
;
318 typedef struct Bigint Bigint
;
320 /* Memory management: memory is allocated from, and returned to, Kmax+1 pools
321 of memory, where pool k (0 <= k <= Kmax) is for Bigints b with b->maxwds ==
322 1 << k. These pools are maintained as linked lists, with freelist[k]
323 pointing to the head of the list for pool k.
325 On allocation, if there's no free slot in the appropriate pool, MALLOC is
326 called to get more memory. This memory is not returned to the system until
327 Python quits. There's also a private memory pool that's allocated from
328 in preference to using MALLOC.
330 For Bigints with more than (1 << Kmax) digits (which implies at least 1233
331 decimal digits), memory is directly allocated using MALLOC, and freed using
334 XXX: it would be easy to bypass this memory-management system and
335 translate each call to Balloc into a call to PyMem_Malloc, and each
336 Bfree to PyMem_Free. Investigate whether this has any significant
337 performance on impact. */
339 static Bigint
*freelist
[Kmax
+1];
341 /* Allocate space for a Bigint with up to 1<<k digits */
350 if (k
<= Kmax
&& (rv
= freelist
[k
]))
351 freelist
[k
] = rv
->next
;
354 len
= (sizeof(Bigint
) + (x
-1)*sizeof(ULong
) + sizeof(double) - 1)
356 if (k
<= Kmax
&& pmem_next
- private_mem
+ len
<= PRIVATE_mem
) {
357 rv
= (Bigint
*)pmem_next
;
361 rv
= (Bigint
*)MALLOC(len
*sizeof(double));
368 rv
->sign
= rv
->wds
= 0;
372 /* Free a Bigint allocated with Balloc */
381 v
->next
= freelist
[v
->k
];
387 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
388 y->wds*sizeof(Long) + 2*sizeof(int))
390 /* Multiply a Bigint b by m and add a. Either modifies b in place and returns
391 a pointer to the modified b, or Bfrees b and returns a pointer to a copy.
392 On failure, return NULL. In this case, b will have been already freed. */
395 multadd(Bigint
*b
, int m
, int a
) /* multiply by m and add a */
413 y
= *x
* (ULLong
)m
+ carry
;
415 *x
++ = (ULong
)(y
& FFFFFFFF
);
418 y
= (xi
& 0xffff) * m
+ carry
;
419 z
= (xi
>> 16) * m
+ (y
>> 16);
421 *x
++ = (z
<< 16) + (y
& 0xffff);
426 if (wds
>= b
->maxwds
) {
436 b
->x
[wds
++] = (ULong
)carry
;
442 /* convert a string s containing nd decimal digits (possibly containing a
443 decimal separator at position nd0, which is ignored) to a Bigint. This
444 function carries on where the parsing code in _Py_dg_strtod leaves off: on
445 entry, y9 contains the result of converting the first 9 digits. Returns
449 s2b(const char *s
, int nd0
, int nd
, ULong y9
, int dplen
)
456 for(k
= 0, y
= 1; x
> y
; y
<<= 1, k
++) ;
467 b
= multadd(b
, 10, *s
++ - '0');
476 b
= multadd(b
, 10, *s
++ - '0');
483 /* count leading 0 bits in the 32-bit integer x. */
490 if (!(x
& 0xffff0000)) {
494 if (!(x
& 0xff000000)) {
498 if (!(x
& 0xf0000000)) {
502 if (!(x
& 0xc0000000)) {
506 if (!(x
& 0x80000000)) {
508 if (!(x
& 0x40000000))
514 /* count trailing 0 bits in the 32-bit integer y, and shift y right by that
560 /* convert a small nonnegative integer to a Bigint */
575 /* multiply two Bigints. Returns a new Bigint, or NULL on failure. Ignores
576 the signs of a and b. */
579 mult(Bigint
*a
, Bigint
*b
)
583 ULong
*x
, *xa
, *xae
, *xb
, *xbe
, *xc
, *xc0
;
592 if (a
->wds
< b
->wds
) {
606 for(x
= c
->x
, xa
= x
+ wc
; x
< xa
; x
++)
614 for(; xb
< xbe
; xc0
++) {
620 z
= *x
++ * (ULLong
)y
+ *xc
+ carry
;
622 *xc
++ = (ULong
)(z
& FFFFFFFF
);
629 for(; xb
< xbe
; xb
++, xc0
++) {
630 if (y
= *xb
& 0xffff) {
635 z
= (*x
& 0xffff) * y
+ (*xc
& 0xffff) + carry
;
637 z2
= (*x
++ >> 16) * y
+ (*xc
>> 16) + carry
;
650 z
= (*x
& 0xffff) * y
+ (*xc
>> 16) + carry
;
653 z2
= (*x
++ >> 16) * y
+ (*xc
& 0xffff) + carry
;
661 for(xc0
= c
->x
, xc
= xc0
+ wc
; wc
> 0 && !*--xc
; --wc
) ;
666 /* p5s is a linked list of powers of 5 of the form 5**(2**i), i >= 2 */
670 /* multiply the Bigint b by 5**k. Returns a pointer to the result, or NULL on
671 failure; if the returned pointer is distinct from b then the original
672 Bigint b will have been Bfree'd. Ignores the sign of b. */
675 pow5mult(Bigint
*b
, int k
)
677 Bigint
*b1
, *p5
, *p51
;
679 static int p05
[3] = { 5, 25, 125 };
682 b
= multadd(b
, p05
[i
-1], 0);
725 /* shift a Bigint b left by k bits. Return a pointer to the shifted result,
726 or NULL on failure. If the returned pointer is distinct from b then the
727 original b will have been Bfree'd. Ignores the sign of b. */
730 lshift(Bigint
*b
, int k
)
734 ULong
*x
, *x1
, *xe
, z
;
739 for(i
= b
->maxwds
; n1
> i
; i
<<= 1)
747 for(i
= 0; i
< n
; i
++)
770 /* Do a three-way compare of a and b, returning -1 if a < b, 0 if a == b and
771 1 if a > b. Ignores signs of a and b. */
774 cmp(Bigint
*a
, Bigint
*b
)
776 ULong
*xa
, *xa0
, *xb
, *xb0
;
782 if (i
> 1 && !a
->x
[i
-1])
783 Bug("cmp called with a->x[a->wds-1] == 0");
784 if (j
> 1 && !b
->x
[j
-1])
785 Bug("cmp called with b->x[b->wds-1] == 0");
795 return *xa
< *xb
? -1 : 1;
802 /* Take the difference of Bigints a and b, returning a new Bigint. Returns
803 NULL on failure. The signs of a and b are ignored, but the sign of the
804 result is set appropriately. */
807 diff(Bigint
*a
, Bigint
*b
)
811 ULong
*xa
, *xae
, *xb
, *xbe
, *xc
;
850 y
= (ULLong
)*xa
++ - *xb
++ - borrow
;
851 borrow
= y
>> 32 & (ULong
)1;
852 *xc
++ = (ULong
)(y
& FFFFFFFF
);
857 borrow
= y
>> 32 & (ULong
)1;
858 *xc
++ = (ULong
)(y
& FFFFFFFF
);
862 y
= (*xa
& 0xffff) - (*xb
& 0xffff) - borrow
;
863 borrow
= (y
& 0x10000) >> 16;
864 z
= (*xa
++ >> 16) - (*xb
++ >> 16) - borrow
;
865 borrow
= (z
& 0x10000) >> 16;
870 y
= (*xa
& 0xffff) - borrow
;
871 borrow
= (y
& 0x10000) >> 16;
872 z
= (*xa
++ >> 16) - borrow
;
873 borrow
= (z
& 0x10000) >> 16;
883 /* Given a positive normal double x, return the difference between x and the next
884 double up. Doesn't give correct results for subnormals. */
892 L
= (word0(x
) & Exp_mask
) - (P
-1)*Exp_msk1
;
898 /* Convert a Bigint to a double plus an exponent */
901 b2d(Bigint
*a
, int *e
)
903 ULong
*xa
, *xa0
, w
, y
, z
;
911 if (!y
) Bug("zero y in b2d");
916 word0(&d
) = Exp_1
| y
>> (Ebits
- k
);
917 w
= xa
> xa0
? *--xa
: 0;
918 word1(&d
) = y
<< ((32-Ebits
) + k
) | w
>> (Ebits
- k
);
921 z
= xa
> xa0
? *--xa
: 0;
923 word0(&d
) = Exp_1
| y
<< k
| z
>> (32 - k
);
924 y
= xa
> xa0
? *--xa
: 0;
925 word1(&d
) = z
<< k
| y
>> (32 - k
);
928 word0(&d
) = Exp_1
| y
;
935 /* Convert a double to a Bigint plus an exponent. Return NULL on failure.
937 Given a finite nonzero double d, return an odd Bigint b and exponent *e
938 such that fabs(d) = b * 2**e. On return, *bbits gives the number of
939 significant bits of e; that is, 2**(*bbits-1) <= b < 2**(*bbits).
941 If d is zero, then b == 0, *e == -1010, *bbits = 0.
946 d2b(U
*d
, int *e
, int *bits
)
958 z
= word0(d
) & Frac_mask
;
959 word0(d
) &= 0x7fffffff; /* clear sign bit, which we ignore */
960 if ((de
= (int)(word0(d
) >> Exp_shift
)))
962 if ((y
= word1(d
))) {
963 if ((k
= lo0bits(&y
))) {
964 x
[0] = y
| z
<< (32 - k
);
970 b
->wds
= (x
[1] = z
) ? 2 : 1;
980 *e
= de
- Bias
- (P
-1) + k
;
984 *e
= de
- Bias
- (P
-1) + 1 + k
;
985 *bits
= 32*i
- hi0bits(x
[i
-1]);
990 /* Compute the ratio of two Bigints, as a double. The result may have an
991 error of up to 2.5 ulps. */
994 ratio(Bigint
*a
, Bigint
*b
)
999 dval(&da
) = b2d(a
, &ka
);
1000 dval(&db
) = b2d(b
, &kb
);
1001 k
= ka
- kb
+ 32*(a
->wds
- b
->wds
);
1003 word0(&da
) += k
*Exp_msk1
;
1006 word0(&db
) += k
*Exp_msk1
;
1008 return dval(&da
) / dval(&db
);
1013 1e0
, 1e1
, 1e2
, 1e3
, 1e4
, 1e5
, 1e6
, 1e7
, 1e8
, 1e9
,
1014 1e10
, 1e11
, 1e12
, 1e13
, 1e14
, 1e15
, 1e16
, 1e17
, 1e18
, 1e19
,
1019 bigtens
[] = { 1e16
, 1e32
, 1e64
, 1e128
, 1e256
};
1020 static const double tinytens
[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1021 9007199254740992.*9007199254740992.e
-256
1022 /* = 2^106 * 1e-256 */
1024 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1025 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1026 #define Scale_Bit 0x10
1029 /* case insensitive string match, for recognising 'inf[inity]' and
1033 match(const char **sp
, char *t
)
1036 const char *s
= *sp
;
1039 if ((c
= *++s
) >= 'A' && c
<= 'Z')
1054 dshift(Bigint
*b
, int p2
)
1056 int rv
= hi0bits(b
->x
[b
->wds
-1]) - 4;
1062 /* special case of Bigint division. The quotient is always in the range 0 <=
1063 quotient < 10, and on entry the divisor S is normalized so that its top 4
1064 bits (28--31) are zero and bit 27 is set. */
1067 quorem(Bigint
*b
, Bigint
*S
)
1070 ULong
*bx
, *bxe
, q
, *sx
, *sxe
;
1072 ULLong borrow
, carry
, y
, ys
;
1074 ULong borrow
, carry
, y
, ys
;
1080 /*debug*/ if (b
->wds
> n
)
1081 /*debug*/ Bug("oversize b in quorem");
1089 q
= *bxe
/ (*sxe
+ 1); /* ensure q <= true quotient */
1091 /*debug*/ if (q
> 9)
1092 /*debug*/ Bug("oversized quotient in quorem");
1099 ys
= *sx
++ * (ULLong
)q
+ carry
;
1101 y
= *bx
- (ys
& FFFFFFFF
) - borrow
;
1102 borrow
= y
>> 32 & (ULong
)1;
1103 *bx
++ = (ULong
)(y
& FFFFFFFF
);
1106 ys
= (si
& 0xffff) * q
+ carry
;
1107 zs
= (si
>> 16) * q
+ (ys
>> 16);
1109 y
= (*bx
& 0xffff) - (ys
& 0xffff) - borrow
;
1110 borrow
= (y
& 0x10000) >> 16;
1111 z
= (*bx
>> 16) - (zs
& 0xffff) - borrow
;
1112 borrow
= (z
& 0x10000) >> 16;
1119 while(--bxe
> bx
&& !*bxe
)
1124 if (cmp(b
, S
) >= 0) {
1134 y
= *bx
- (ys
& FFFFFFFF
) - borrow
;
1135 borrow
= y
>> 32 & (ULong
)1;
1136 *bx
++ = (ULong
)(y
& FFFFFFFF
);
1139 ys
= (si
& 0xffff) + carry
;
1140 zs
= (si
>> 16) + (ys
>> 16);
1142 y
= (*bx
& 0xffff) - (ys
& 0xffff) - borrow
;
1143 borrow
= (y
& 0x10000) >> 16;
1144 z
= (*bx
>> 16) - (zs
& 0xffff) - borrow
;
1145 borrow
= (z
& 0x10000) >> 16;
1153 while(--bxe
> bx
&& !*bxe
)
1162 /* return 0 on success, -1 on failure */
1165 bigcomp(U
*rv
, const char *s0
, BCinfo
*bc
)
1168 int b2
, bbits
, d2
, dd
, dig
, dsign
, i
, j
, nd
, nd0
, p2
, p5
, speccase
;
1173 p5
= nd
+ bc
->e0
- 1;
1175 if (rv
->d
== 0.) { /* special case: value near underflow-to-zero */
1176 /* threshold was rounded to zero */
1182 word0(rv
) = (P
+2) << Exp_shift
;
1193 b
= d2b(rv
, &p2
, &bbits
);
1198 /* floor(log2(rv)) == bbits - 1 + p2 */
1199 /* Check for denormal case. */
1201 if (i
> (j
= P
- Emin
- 1 + p2
)) {
1217 /* Arrange for convenient computation of quotients:
1218 * shift left if necessary so divisor has 4 leading 0 bits.
1221 d
= pow5mult(d
, p5
);
1228 b
= pow5mult(b
, -p5
);
1243 if ((b2
+= i
) > 0) {
1250 if ((d2
+= i
) > 0) {
1258 /* Now b/d = exactly half-way between the two floating-point values */
1259 /* on either side of the input string. Compute first digit of b/d. */
1261 if (!(dig
= quorem(b
,d
))) {
1262 b
= multadd(b
, 10, 0); /* very unlikely */
1270 /* Compare b/d with s0 */
1273 dd
= 9999; /* silence gcc compiler warning */
1274 for(i
= 0; i
< nd0
; ) {
1275 if ((dd
= s0
[i
++] - '0' - dig
))
1277 if (!b
->x
[0] && b
->wds
== 1) {
1282 b
= multadd(b
, 10, 0);
1289 for(j
= bc
->dp1
; i
++ < nd
;) {
1290 if ((dd
= s0
[j
++] - '0' - dig
))
1292 if (!b
->x
[0] && b
->wds
== 1) {
1297 b
= multadd(b
, 10, 0);
1304 if (b
->x
[0] || b
->wds
> 1)
1314 if (!dsign
) /* does not happen for round-near */
1316 dval(rv
) -= ulp(rv
);
1321 dval(rv
) += ulp(rv
);
1325 /* Exact half-way case: apply round-even rule. */
1326 if (word1(rv
) & 1) {
1337 _Py_dg_strtod(const char *s00
, char **se
)
1339 int bb2
, bb5
, bbe
, bd2
, bd5
, bbbits
, bs2
, c
, e
, e1
, error
;
1340 int esign
, i
, j
, k
, nd
, nd0
, nf
, nz
, nz0
, sign
;
1341 const char *s
, *s0
, *s1
;
1344 U aadj2
, adj
, rv
, rv0
;
1347 Bigint
*bb
, *bb1
, *bd
, *bd0
, *bs
, *delta
;
1349 sign
= nz0
= nz
= bc
.dplen
= bc
.uflchk
= 0;
1351 for(s
= s00
;;s
++) switch(*s
) {
1361 /* modify original dtoa.c so that it doesn't accept leading whitespace
1376 while(*++s
== '0') ;
1382 for(nd
= nf
= 0; (c
= *s
) >= '0' && c
<= '9'; nd
++, s
++)
1388 bc
.dp0
= bc
.dp1
= s
- s0
;
1392 bc
.dplen
= bc
.dp1
- bc
.dp0
;
1394 for(; c
== '0'; c
= *++s
)
1396 if (c
> '0' && c
<= '9') {
1404 for(; c
>= '0' && c
<= '9'; c
= *++s
) {
1409 for(i
= 1; i
< nz
; i
++)
1412 else if (nd
<= DBL_DIG
+ 1)
1416 else if (nd
<= DBL_DIG
+ 1)
1424 if (c
== 'e' || c
== 'E') {
1425 if (!nd
&& !nz
&& !nz0
) {
1436 if (c
>= '0' && c
<= '9') {
1439 if (c
> '0' && c
<= '9') {
1442 while((c
= *++s
) >= '0' && c
<= '9')
1444 if (s
- s1
> 8 || L
> 19999)
1445 /* Avoid confusion from exponents
1446 * so large that e might overflow.
1448 e
= 19999; /* safe for 16 bit ints */
1462 /* Check for Nan and Infinity */
1467 if (match(&s
,"nf")) {
1469 if (!match(&s
,"inity"))
1471 word0(&rv
) = 0x7ff00000;
1478 if (match(&s
, "an")) {
1479 word0(&rv
) = NAN_WORD0
;
1480 word1(&rv
) = NAN_WORD1
;
1490 bc
.e0
= e1
= e
-= nf
;
1492 /* Now we have nd0 digits, starting at s0, followed by a
1493 * decimal point, followed by nd-nd0 digits. The number we're
1494 * after is the integer represented by those digits times
1499 k
= nd
< DBL_DIG
+ 1 ? nd
: DBL_DIG
+ 1;
1502 dval(&rv
) = tens
[k
- 9] * dval(&rv
) + z
;
1511 if (e
<= Ten_pmax
) {
1512 dval(&rv
) *= tens
[e
];
1516 if (e
<= Ten_pmax
+ i
) {
1517 /* A fancier test would sometimes let us do
1518 * this for larger i values.
1521 dval(&rv
) *= tens
[i
];
1522 dval(&rv
) *= tens
[e
];
1526 else if (e
>= -Ten_pmax
) {
1527 dval(&rv
) /= tens
[-e
];
1535 /* Get starting approximation = rv * 10**e1 */
1539 dval(&rv
) *= tens
[i
];
1541 if (e1
> DBL_MAX_10_EXP
) {
1544 /* Can't trust HUGE_VAL */
1545 word0(&rv
) = Exp_mask
;
1550 for(j
= 0; e1
> 1; j
++, e1
>>= 1)
1552 dval(&rv
) *= bigtens
[j
];
1553 /* The last multiplication could overflow. */
1554 word0(&rv
) -= P
*Exp_msk1
;
1555 dval(&rv
) *= bigtens
[j
];
1556 if ((z
= word0(&rv
) & Exp_mask
)
1557 > Exp_msk1
*(DBL_MAX_EXP
+Bias
-P
))
1559 if (z
> Exp_msk1
*(DBL_MAX_EXP
+Bias
-1-P
)) {
1560 /* set to largest number */
1561 /* (Can't trust DBL_MAX) */
1566 word0(&rv
) += P
*Exp_msk1
;
1572 dval(&rv
) /= tens
[i
];
1574 if (e1
>= 1 << n_bigtens
)
1578 for(j
= 0; e1
> 0; j
++, e1
>>= 1)
1580 dval(&rv
) *= tinytens
[j
];
1581 if (bc
.scale
&& (j
= 2*P
+ 1 - ((word0(&rv
) & Exp_mask
)
1582 >> Exp_shift
)) > 0) {
1583 /* scaled rv is denormal; clear j low bits */
1587 word0(&rv
) = (P
+2)*Exp_msk1
;
1589 word0(&rv
) &= 0xffffffff << (j
-32);
1592 word1(&rv
) &= 0xffffffff << j
;
1603 /* Now the hard part -- adjusting rv to the correct value.*/
1605 /* Put digits into bd: true value = bd * 10^e */
1608 bc
.nd0
= nd0
; /* Only needed if nd > strtod_diglim, but done here */
1609 /* to silence an erroneous warning about bc.nd0 */
1610 /* possibly not being initialized. */
1611 if (nd
> strtod_diglim
) {
1612 /* ASSERT(strtod_diglim >= 18); 18 == one more than the */
1613 /* minimum number of decimal digits to distinguish double values */
1614 /* in IEEE arithmetic. */
1619 if (--j
<= bc
.dp1
&& j
>= bc
.dp0
)
1629 if (nd
< 9) { /* must recompute y */
1631 for(i
= 0; i
< nd0
; ++i
)
1632 y
= 10*y
+ s0
[i
] - '0';
1633 for(j
= bc
.dp1
; i
< nd
; ++i
)
1634 y
= 10*y
+ s0
[j
++] - '0';
1637 bd0
= s2b(s0
, nd0
, nd
, y
, bc
.dplen
);
1642 bd
= Balloc(bd0
->k
);
1648 bb
= d2b(&rv
, &bbe
, &bbbits
); /* rv = bb * 2^bbe */
1676 i
= j
+ bbbits
- 1; /* logb(rv) */
1677 if (i
< Emin
) /* denormal */
1684 i
= bb2
< bd2
? bb2
: bd2
;
1693 bs
= pow5mult(bs
, bb5
);
1711 bb
= lshift(bb
, bb2
);
1720 bd
= pow5mult(bd
, bd5
);
1729 bd
= lshift(bd
, bd2
);
1738 bs
= lshift(bs
, bs2
);
1746 delta
= diff(bb
, bd
);
1747 if (delta
== NULL
) {
1754 bc
.dsign
= delta
->sign
;
1757 if (bc
.nd
> nd
&& i
<= 0) {
1759 break; /* Must use bigcomp(). */
1762 i
= -1; /* Discarded digits make delta smaller. */
1767 /* Error is less than half an ulp -- check for
1768 * special case of mantissa a power of two.
1770 if (bc
.dsign
|| word1(&rv
) || word0(&rv
) & Bndry_mask
1771 || (word0(&rv
) & Exp_mask
) <= (2*P
+1)*Exp_msk1
1775 if (!delta
->x
[0] && delta
->wds
<= 1) {
1779 delta
= lshift(delta
,Log2P
);
1780 if (delta
== NULL
) {
1787 if (cmp(delta
, bs
) > 0)
1792 /* exactly half-way between */
1794 if ((word0(&rv
) & Bndry_mask1
) == Bndry_mask1
1797 (y
= word0(&rv
) & Exp_mask
) <= 2*P
*Exp_msk1
) ?
1798 (0xffffffff & (0xffffffff << (2*P
+1-(y
>>Exp_shift
)))) :
1800 /*boundary case -- increment exponent*/
1801 word0(&rv
) = (word0(&rv
) & Exp_mask
)
1809 else if (!(word0(&rv
) & Bndry_mask
) && !word1(&rv
)) {
1811 /* boundary case -- decrement exponent */
1813 L
= word0(&rv
) & Exp_mask
;
1814 if (L
<= (2*P
+1)*Exp_msk1
) {
1815 if (L
> (P
+2)*Exp_msk1
)
1816 /* round even ==> */
1819 /* rv = smallest denormal */
1827 L
= (word0(&rv
) & Exp_mask
) - Exp_msk1
;
1828 word0(&rv
) = L
| Bndry_mask1
;
1829 word1(&rv
) = 0xffffffff;
1832 if (!(word1(&rv
) & LSB
))
1835 dval(&rv
) += ulp(&rv
);
1837 dval(&rv
) -= ulp(&rv
);
1846 bc
.dsign
= 1 - bc
.dsign
;
1849 if ((aadj
= ratio(delta
, bs
)) <= 2.) {
1852 else if (word1(&rv
) || word0(&rv
) & Bndry_mask
) {
1853 if (word1(&rv
) == Tiny1
&& !word0(&rv
)) {
1864 /* special case -- power of FLT_RADIX to be */
1865 /* rounded down... */
1867 if (aadj
< 2./FLT_RADIX
)
1868 aadj
= 1./FLT_RADIX
;
1876 aadj1
= bc
.dsign
? aadj
: -aadj
;
1877 if (Flt_Rounds
== 0)
1880 y
= word0(&rv
) & Exp_mask
;
1882 /* Check for overflow */
1884 if (y
== Exp_msk1
*(DBL_MAX_EXP
+Bias
-1)) {
1885 dval(&rv0
) = dval(&rv
);
1886 word0(&rv
) -= P
*Exp_msk1
;
1887 adj
.d
= aadj1
* ulp(&rv
);
1889 if ((word0(&rv
) & Exp_mask
) >=
1890 Exp_msk1
*(DBL_MAX_EXP
+Bias
-P
)) {
1891 if (word0(&rv0
) == Big0
&& word1(&rv0
) == Big1
)
1898 word0(&rv
) += P
*Exp_msk1
;
1901 if (bc
.scale
&& y
<= 2*P
*Exp_msk1
) {
1902 if (aadj
<= 0x7fffffff) {
1903 if ((z
= (ULong
)aadj
) <= 0)
1906 aadj1
= bc
.dsign
? aadj
: -aadj
;
1908 dval(&aadj2
) = aadj1
;
1909 word0(&aadj2
) += (2*P
+1)*Exp_msk1
- y
;
1910 aadj1
= dval(&aadj2
);
1912 adj
.d
= aadj1
* ulp(&rv
);
1915 z
= word0(&rv
) & Exp_mask
;
1919 /* Can we stop now? */
1922 /* The tolerances below are conservative. */
1923 if (bc
.dsign
|| word1(&rv
) || word0(&rv
) & Bndry_mask
) {
1924 if (aadj
< .4999999 || aadj
> .5000001)
1927 else if (aadj
< .4999999/FLT_RADIX
)
1943 error
= bigcomp(&rv
, s0
, &bc
);
1949 word0(&rv0
) = Exp_1
- 2*P
*Exp_msk1
;
1951 dval(&rv
) *= dval(&rv0
);
1952 /* try to avoid the bug of testing an 8087 register value */
1953 if (!(word0(&rv
) & Exp_mask
))
1959 return sign
? -dval(&rv
) : dval(&rv
);
1975 sizeof(Bigint
) - sizeof(ULong
) - sizeof(int) + j
<= (unsigned)i
;
1978 r
= (int*)Balloc(k
);
1982 return (char *)(r
+1);
1986 nrv_alloc(char *s
, char **rve
, int n
)
1994 while((*t
= *s
++)) t
++;
2000 /* freedtoa(s) must be used to free values s returned by dtoa
2001 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
2002 * but for consistency with earlier versions of dtoa, it is optional
2003 * when MULTIPLE_THREADS is not defined.
2007 _Py_dg_freedtoa(char *s
)
2009 Bigint
*b
= (Bigint
*)((int *)s
- 1);
2010 b
->maxwds
= 1 << (b
->k
= *(int*)b
);
2014 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
2016 * Inspired by "How to Print Floating-Point Numbers Accurately" by
2017 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
2020 * 1. Rather than iterating, we use a simple numeric overestimate
2021 * to determine k = floor(log10(d)). We scale relevant
2022 * quantities using O(log2(k)) rather than O(k) multiplications.
2023 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
2024 * try to generate digits strictly left to right. Instead, we
2025 * compute with fewer bits and propagate the carry if necessary
2026 * when rounding the final digit up. This is often faster.
2027 * 3. Under the assumption that input will be rounded nearest,
2028 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
2029 * That is, we allow equality in stopping tests when the
2030 * round-nearest rule will give the same floating-point value
2031 * as would satisfaction of the stopping test with strict
2033 * 4. We remove common factors of powers of 2 from relevant
2035 * 5. When converting floating-point integers less than 1e16,
2036 * we use floating-point arithmetic rather than resorting
2037 * to multiple-precision integers.
2038 * 6. When asked to produce fewer than 15 digits, we first try
2039 * to get by with floating-point arithmetic; we resort to
2040 * multiple-precision integer arithmetic only if we cannot
2041 * guarantee that the floating-point calculation has given
2042 * the correctly rounded result. For k requested digits and
2043 * "uniformly" distributed input, the probability is
2044 * something like 10^(k-15) that we must resort to the Long
2048 /* Additional notes (METD): (1) returns NULL on failure. (2) to avoid memory
2049 leakage, a successful call to _Py_dg_dtoa should always be matched by a
2050 call to _Py_dg_freedtoa. */
2053 _Py_dg_dtoa(double dd
, int mode
, int ndigits
,
2054 int *decpt
, int *sign
, char **rve
)
2056 /* Arguments ndigits, decpt, sign are similar to those
2057 of ecvt and fcvt; trailing zeros are suppressed from
2058 the returned string. If not null, *rve is set to point
2059 to the end of the return value. If d is +-Infinity or NaN,
2060 then *decpt is set to 9999.
2063 0 ==> shortest string that yields d when read in
2064 and rounded to nearest.
2065 1 ==> like 0, but with Steele & White stopping rule;
2066 e.g. with IEEE P754 arithmetic , mode 0 gives
2067 1e23 whereas mode 1 gives 9.999999999999999e22.
2068 2 ==> max(1,ndigits) significant digits. This gives a
2069 return value similar to that of ecvt, except
2070 that trailing zeros are suppressed.
2071 3 ==> through ndigits past the decimal point. This
2072 gives a return value similar to that from fcvt,
2073 except that trailing zeros are suppressed, and
2074 ndigits can be negative.
2075 4,5 ==> similar to 2 and 3, respectively, but (in
2076 round-nearest mode) with the tests of mode 0 to
2077 possibly return a shorter string that rounds to d.
2078 With IEEE arithmetic and compilation with
2079 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
2080 as modes 2 and 3 when FLT_ROUNDS != 1.
2081 6-9 ==> Debugging modes similar to mode - 4: don't try
2082 fast floating-point estimate (if applicable).
2084 Values of mode other than 0-9 are treated as mode 0.
2086 Sufficient space is allocated to the return value
2087 to hold the suppressed trailing zeros.
2090 int bbits
, b2
, b5
, be
, dig
, i
, ieps
, ilim
, ilim0
, ilim1
,
2091 j
, j1
, k
, k0
, k_check
, leftright
, m2
, m5
, s2
, s5
,
2092 spec_case
, try_quick
;
2096 Bigint
*b
, *b1
, *delta
, *mlo
, *mhi
, *S
;
2101 /* set pointers to NULL, to silence gcc compiler warnings and make
2102 cleanup easier on error */
2103 mlo
= mhi
= b
= S
= 0;
2107 if (word0(&u
) & Sign_bit
) {
2108 /* set sign for everything, including 0's and NaNs */
2110 word0(&u
) &= ~Sign_bit
; /* clear sign bit */
2115 /* quick return for Infinities, NaNs and zeros */
2116 if ((word0(&u
) & Exp_mask
) == Exp_mask
)
2118 /* Infinity or NaN */
2120 if (!word1(&u
) && !(word0(&u
) & 0xfffff))
2121 return nrv_alloc("Infinity", rve
, 8);
2122 return nrv_alloc("NaN", rve
, 3);
2126 return nrv_alloc("0", rve
, 1);
2129 /* compute k = floor(log10(d)). The computation may leave k
2130 one too large, but should never leave k too small. */
2131 b
= d2b(&u
, &be
, &bbits
);
2134 if ((i
= (int)(word0(&u
) >> Exp_shift1
& (Exp_mask
>>Exp_shift1
)))) {
2135 dval(&d2
) = dval(&u
);
2136 word0(&d2
) &= Frac_mask1
;
2137 word0(&d2
) |= Exp_11
;
2139 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2140 * log10(x) = log(x) / log(10)
2141 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2142 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2144 * This suggests computing an approximation k to log10(d) by
2146 * k = (i - Bias)*0.301029995663981
2147 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2149 * We want k to be too large rather than too small.
2150 * The error in the first-order Taylor series approximation
2151 * is in our favor, so we just round up the constant enough
2152 * to compensate for any error in the multiplication of
2153 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2154 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2155 * adding 1e-13 to the constant term more than suffices.
2156 * Hence we adjust the constant term to 0.1760912590558.
2157 * (We could get a more accurate k by invoking log10,
2158 * but this is probably not worthwhile.)
2165 /* d is denormalized */
2167 i
= bbits
+ be
+ (Bias
+ (P
-1) - 1);
2168 x
= i
> 32 ? word0(&u
) << (64 - i
) | word1(&u
) >> (i
- 32)
2169 : word1(&u
) << (32 - i
);
2171 word0(&d2
) -= 31*Exp_msk1
; /* adjust exponent */
2172 i
-= (Bias
+ (P
-1) - 1) + 1;
2175 ds
= (dval(&d2
)-1.5)*0.289529654602168 + 0.1760912590558 +
2176 i
*0.301029995663981;
2178 if (ds
< 0. && ds
!= k
)
2179 k
--; /* want k = floor(ds) */
2181 if (k
>= 0 && k
<= Ten_pmax
) {
2182 if (dval(&u
) < tens
[k
])
2205 if (mode
< 0 || mode
> 9)
2215 ilim
= ilim1
= -1; /* Values for cases 0 and 1; done here to */
2216 /* silence erroneous "gcc -Wall" warning. */
2229 ilim
= ilim1
= i
= ndigits
;
2235 i
= ndigits
+ k
+ 1;
2247 if (ilim
>= 0 && ilim
<= Quick_max
&& try_quick
) {
2249 /* Try to get by with floating-point arithmetic. */
2252 dval(&d2
) = dval(&u
);
2255 ieps
= 2; /* conservative */
2260 /* prevent overflows */
2262 dval(&u
) /= bigtens
[n_bigtens
-1];
2265 for(; j
; j
>>= 1, i
++)
2272 else if ((j1
= -k
)) {
2273 dval(&u
) *= tens
[j1
& 0xf];
2274 for(j
= j1
>> 4; j
; j
>>= 1, i
++)
2277 dval(&u
) *= bigtens
[i
];
2280 if (k_check
&& dval(&u
) < 1. && ilim
> 0) {
2288 dval(&eps
) = ieps
*dval(&u
) + 7.;
2289 word0(&eps
) -= (P
-1)*Exp_msk1
;
2293 if (dval(&u
) > dval(&eps
))
2295 if (dval(&u
) < -dval(&eps
))
2300 /* Use Steele & White method of only
2301 * generating digits needed.
2303 dval(&eps
) = 0.5/tens
[ilim
-1] - dval(&eps
);
2307 *s
++ = '0' + (int)L
;
2308 if (dval(&u
) < dval(&eps
))
2310 if (1. - dval(&u
) < dval(&eps
))
2319 /* Generate ilim digits, then fix them up. */
2320 dval(&eps
) *= tens
[ilim
-1];
2321 for(i
= 1;; i
++, dval(&u
) *= 10.) {
2322 L
= (Long
)(dval(&u
));
2323 if (!(dval(&u
) -= L
))
2325 *s
++ = '0' + (int)L
;
2327 if (dval(&u
) > 0.5 + dval(&eps
))
2329 else if (dval(&u
) < 0.5 - dval(&eps
)) {
2340 dval(&u
) = dval(&d2
);
2345 /* Do we have a "small" integer? */
2347 if (be
>= 0 && k
<= Int_max
) {
2350 if (ndigits
< 0 && ilim
<= 0) {
2352 if (ilim
< 0 || dval(&u
) <= 5*ds
)
2356 for(i
= 1;; i
++, dval(&u
) *= 10.) {
2357 L
= (Long
)(dval(&u
) / ds
);
2359 *s
++ = '0' + (int)L
;
2364 dval(&u
) += dval(&u
);
2365 if (dval(&u
) > ds
|| (dval(&u
) == ds
&& L
& 1)) {
2385 denorm
? be
+ (Bias
+ (P
-1) - 1 + 1) :
2393 if (m2
> 0 && s2
> 0) {
2394 i
= m2
< s2
? m2
: s2
;
2402 mhi
= pow5mult(mhi
, m5
);
2411 if ((j
= b5
- m5
)) {
2418 b
= pow5mult(b
, b5
);
2427 S
= pow5mult(S
, s5
);
2432 /* Check for special case that d is a normalized power of 2. */
2435 if ((mode
< 2 || leftright
)
2437 if (!word1(&u
) && !(word0(&u
) & Bndry_mask
)
2438 && word0(&u
) & (Exp_mask
& ~Exp_msk1
)
2440 /* The special case */
2447 /* Arrange for convenient computation of quotients:
2448 * shift left if necessary so divisor has 4 leading 0 bits.
2450 * Perhaps we should just compute leading 28 bits of S once
2451 * and for all and pass them and a shift to quorem, so it
2452 * can do shifts and ors to compute the numerator for q.
2454 if ((i
= ((s5
? 32 - hi0bits(S
->x
[S
->wds
-1]) : 1) + s2
) & 0x1f))
2474 b
= multadd(b
, 10, 0); /* we botched the k estimate */
2478 mhi
= multadd(mhi
, 10, 0);
2485 if (ilim
<= 0 && (mode
== 3 || mode
== 5)) {
2487 /* no digits, fcvt style */
2493 S
= multadd(S
, 5, 0);
2506 mhi
= lshift(mhi
, m2
);
2511 /* Compute mlo -- check for special case
2512 * that d is a normalized power of 2.
2517 mhi
= Balloc(mhi
->k
);
2521 mhi
= lshift(mhi
, Log2P
);
2527 dig
= quorem(b
,S
) + '0';
2528 /* Do we yet have the shortest decimal string
2529 * that will round to d?
2532 delta
= diff(S
, mhi
);
2535 j1
= delta
->sign
? 1 : cmp(b
, delta
);
2537 if (j1
== 0 && mode
!= 1 && !(word1(&u
) & 1)
2546 if (j
< 0 || (j
== 0 && mode
!= 1
2549 if (!b
->x
[0] && b
->wds
<= 1) {
2557 if ((j1
> 0 || (j1
== 0 && dig
& 1))
2566 if (dig
== '9') { /* possible if i == 1 */
2577 b
= multadd(b
, 10, 0);
2581 mlo
= mhi
= multadd(mhi
, 10, 0);
2586 mlo
= multadd(mlo
, 10, 0);
2589 mhi
= multadd(mhi
, 10, 0);
2597 *s
++ = dig
= quorem(b
,S
) + '0';
2598 if (!b
->x
[0] && b
->wds
<= 1) {
2603 b
= multadd(b
, 10, 0);
2608 /* Round off last digit */
2614 if (j
> 0 || (j
== 0 && dig
& 1)) {
2631 if (mlo
&& mlo
!= mhi
)
2645 if (mlo
&& mlo
!= mhi
)
2652 _Py_dg_freedtoa(s0
);
2659 #endif /* PY_NO_SHORT_FLOAT_REPR */