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 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
)
447 for(k
= 0, y
= 1; x
> y
; y
<<= 1, k
++) ;
458 for (i
= 9; i
< nd0
; i
++) {
459 b
= multadd(b
, 10, *s
++ - '0');
465 b
= multadd(b
, 10, *s
++ - '0');
472 /* count leading 0 bits in the 32-bit integer x. */
479 if (!(x
& 0xffff0000)) {
483 if (!(x
& 0xff000000)) {
487 if (!(x
& 0xf0000000)) {
491 if (!(x
& 0xc0000000)) {
495 if (!(x
& 0x80000000)) {
497 if (!(x
& 0x40000000))
503 /* count trailing 0 bits in the 32-bit integer y, and shift y right by that
549 /* convert a small nonnegative integer to a Bigint */
564 /* multiply two Bigints. Returns a new Bigint, or NULL on failure. Ignores
565 the signs of a and b. */
568 mult(Bigint
*a
, Bigint
*b
)
572 ULong
*x
, *xa
, *xae
, *xb
, *xbe
, *xc
, *xc0
;
581 if (a
->wds
< b
->wds
) {
595 for(x
= c
->x
, xa
= x
+ wc
; x
< xa
; x
++)
603 for(; xb
< xbe
; xc0
++) {
609 z
= *x
++ * (ULLong
)y
+ *xc
+ carry
;
611 *xc
++ = (ULong
)(z
& FFFFFFFF
);
618 for(; xb
< xbe
; xb
++, xc0
++) {
619 if (y
= *xb
& 0xffff) {
624 z
= (*x
& 0xffff) * y
+ (*xc
& 0xffff) + carry
;
626 z2
= (*x
++ >> 16) * y
+ (*xc
>> 16) + carry
;
639 z
= (*x
& 0xffff) * y
+ (*xc
>> 16) + carry
;
642 z2
= (*x
++ >> 16) * y
+ (*xc
& 0xffff) + carry
;
650 for(xc0
= c
->x
, xc
= xc0
+ wc
; wc
> 0 && !*--xc
; --wc
) ;
655 /* p5s is a linked list of powers of 5 of the form 5**(2**i), i >= 2 */
659 /* multiply the Bigint b by 5**k. Returns a pointer to the result, or NULL on
660 failure; if the returned pointer is distinct from b then the original
661 Bigint b will have been Bfree'd. Ignores the sign of b. */
664 pow5mult(Bigint
*b
, int k
)
666 Bigint
*b1
, *p5
, *p51
;
668 static int p05
[3] = { 5, 25, 125 };
671 b
= multadd(b
, p05
[i
-1], 0);
714 /* shift a Bigint b left by k bits. Return a pointer to the shifted result,
715 or NULL on failure. If the returned pointer is distinct from b then the
716 original b will have been Bfree'd. Ignores the sign of b. */
719 lshift(Bigint
*b
, int k
)
723 ULong
*x
, *x1
, *xe
, z
;
728 for(i
= b
->maxwds
; n1
> i
; i
<<= 1)
736 for(i
= 0; i
< n
; i
++)
759 /* Do a three-way compare of a and b, returning -1 if a < b, 0 if a == b and
760 1 if a > b. Ignores signs of a and b. */
763 cmp(Bigint
*a
, Bigint
*b
)
765 ULong
*xa
, *xa0
, *xb
, *xb0
;
771 if (i
> 1 && !a
->x
[i
-1])
772 Bug("cmp called with a->x[a->wds-1] == 0");
773 if (j
> 1 && !b
->x
[j
-1])
774 Bug("cmp called with b->x[b->wds-1] == 0");
784 return *xa
< *xb
? -1 : 1;
791 /* Take the difference of Bigints a and b, returning a new Bigint. Returns
792 NULL on failure. The signs of a and b are ignored, but the sign of the
793 result is set appropriately. */
796 diff(Bigint
*a
, Bigint
*b
)
800 ULong
*xa
, *xae
, *xb
, *xbe
, *xc
;
839 y
= (ULLong
)*xa
++ - *xb
++ - borrow
;
840 borrow
= y
>> 32 & (ULong
)1;
841 *xc
++ = (ULong
)(y
& FFFFFFFF
);
846 borrow
= y
>> 32 & (ULong
)1;
847 *xc
++ = (ULong
)(y
& FFFFFFFF
);
851 y
= (*xa
& 0xffff) - (*xb
& 0xffff) - borrow
;
852 borrow
= (y
& 0x10000) >> 16;
853 z
= (*xa
++ >> 16) - (*xb
++ >> 16) - borrow
;
854 borrow
= (z
& 0x10000) >> 16;
859 y
= (*xa
& 0xffff) - borrow
;
860 borrow
= (y
& 0x10000) >> 16;
861 z
= (*xa
++ >> 16) - borrow
;
862 borrow
= (z
& 0x10000) >> 16;
872 /* Given a positive normal double x, return the difference between x and the next
873 double up. Doesn't give correct results for subnormals. */
881 L
= (word0(x
) & Exp_mask
) - (P
-1)*Exp_msk1
;
887 /* Convert a Bigint to a double plus an exponent */
890 b2d(Bigint
*a
, int *e
)
892 ULong
*xa
, *xa0
, w
, y
, z
;
900 if (!y
) Bug("zero y in b2d");
905 word0(&d
) = Exp_1
| y
>> (Ebits
- k
);
906 w
= xa
> xa0
? *--xa
: 0;
907 word1(&d
) = y
<< ((32-Ebits
) + k
) | w
>> (Ebits
- k
);
910 z
= xa
> xa0
? *--xa
: 0;
912 word0(&d
) = Exp_1
| y
<< k
| z
>> (32 - k
);
913 y
= xa
> xa0
? *--xa
: 0;
914 word1(&d
) = z
<< k
| y
>> (32 - k
);
917 word0(&d
) = Exp_1
| y
;
924 /* Convert a double to a Bigint plus an exponent. Return NULL on failure.
926 Given a finite nonzero double d, return an odd Bigint b and exponent *e
927 such that fabs(d) = b * 2**e. On return, *bbits gives the number of
928 significant bits of b; that is, 2**(*bbits-1) <= b < 2**(*bbits).
930 If d is zero, then b == 0, *e == -1010, *bbits = 0.
935 d2b(U
*d
, int *e
, int *bits
)
947 z
= word0(d
) & Frac_mask
;
948 word0(d
) &= 0x7fffffff; /* clear sign bit, which we ignore */
949 if ((de
= (int)(word0(d
) >> Exp_shift
)))
951 if ((y
= word1(d
))) {
952 if ((k
= lo0bits(&y
))) {
953 x
[0] = y
| z
<< (32 - k
);
959 b
->wds
= (x
[1] = z
) ? 2 : 1;
969 *e
= de
- Bias
- (P
-1) + k
;
973 *e
= de
- Bias
- (P
-1) + 1 + k
;
974 *bits
= 32*i
- hi0bits(x
[i
-1]);
979 /* Compute the ratio of two Bigints, as a double. The result may have an
980 error of up to 2.5 ulps. */
983 ratio(Bigint
*a
, Bigint
*b
)
988 dval(&da
) = b2d(a
, &ka
);
989 dval(&db
) = b2d(b
, &kb
);
990 k
= ka
- kb
+ 32*(a
->wds
- b
->wds
);
992 word0(&da
) += k
*Exp_msk1
;
995 word0(&db
) += k
*Exp_msk1
;
997 return dval(&da
) / dval(&db
);
1002 1e0
, 1e1
, 1e2
, 1e3
, 1e4
, 1e5
, 1e6
, 1e7
, 1e8
, 1e9
,
1003 1e10
, 1e11
, 1e12
, 1e13
, 1e14
, 1e15
, 1e16
, 1e17
, 1e18
, 1e19
,
1008 bigtens
[] = { 1e16
, 1e32
, 1e64
, 1e128
, 1e256
};
1009 static const double tinytens
[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1010 9007199254740992.*9007199254740992.e
-256
1011 /* = 2^106 * 1e-256 */
1013 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1014 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1015 #define Scale_Bit 0x10
1024 dshift(Bigint
*b
, int p2
)
1026 int rv
= hi0bits(b
->x
[b
->wds
-1]) - 4;
1032 /* special case of Bigint division. The quotient is always in the range 0 <=
1033 quotient < 10, and on entry the divisor S is normalized so that its top 4
1034 bits (28--31) are zero and bit 27 is set. */
1037 quorem(Bigint
*b
, Bigint
*S
)
1040 ULong
*bx
, *bxe
, q
, *sx
, *sxe
;
1042 ULLong borrow
, carry
, y
, ys
;
1044 ULong borrow
, carry
, y
, ys
;
1050 /*debug*/ if (b
->wds
> n
)
1051 /*debug*/ Bug("oversize b in quorem");
1059 q
= *bxe
/ (*sxe
+ 1); /* ensure q <= true quotient */
1061 /*debug*/ if (q
> 9)
1062 /*debug*/ Bug("oversized quotient in quorem");
1069 ys
= *sx
++ * (ULLong
)q
+ carry
;
1071 y
= *bx
- (ys
& FFFFFFFF
) - borrow
;
1072 borrow
= y
>> 32 & (ULong
)1;
1073 *bx
++ = (ULong
)(y
& FFFFFFFF
);
1076 ys
= (si
& 0xffff) * q
+ carry
;
1077 zs
= (si
>> 16) * q
+ (ys
>> 16);
1079 y
= (*bx
& 0xffff) - (ys
& 0xffff) - borrow
;
1080 borrow
= (y
& 0x10000) >> 16;
1081 z
= (*bx
>> 16) - (zs
& 0xffff) - borrow
;
1082 borrow
= (z
& 0x10000) >> 16;
1089 while(--bxe
> bx
&& !*bxe
)
1094 if (cmp(b
, S
) >= 0) {
1104 y
= *bx
- (ys
& FFFFFFFF
) - borrow
;
1105 borrow
= y
>> 32 & (ULong
)1;
1106 *bx
++ = (ULong
)(y
& FFFFFFFF
);
1109 ys
= (si
& 0xffff) + carry
;
1110 zs
= (si
>> 16) + (ys
>> 16);
1112 y
= (*bx
& 0xffff) - (ys
& 0xffff) - borrow
;
1113 borrow
= (y
& 0x10000) >> 16;
1114 z
= (*bx
>> 16) - (zs
& 0xffff) - borrow
;
1115 borrow
= (z
& 0x10000) >> 16;
1123 while(--bxe
> bx
&& !*bxe
)
1131 /* sulp(x) is a version of ulp(x) that takes bc.scale into account.
1133 Assuming that x is finite and nonnegative (positive zero is fine
1134 here) and x / 2^bc.scale is exactly representable as a double,
1135 sulp(x) is equivalent to 2^bc.scale * ulp(x / 2^bc.scale). */
1138 sulp(U
*x
, BCinfo
*bc
)
1142 if (bc
->scale
&& 2*P
+ 1 > (int)((word0(x
) & Exp_mask
) >> Exp_shift
)) {
1143 /* rv/2^bc->scale is subnormal */
1144 word0(&u
) = (P
+2)*Exp_msk1
;
1149 assert(word0(x
) || word1(x
)); /* x != 0.0 */
1154 /* The bigcomp function handles some hard cases for strtod, for inputs
1155 with more than STRTOD_DIGLIM digits. It's called once an initial
1156 estimate for the double corresponding to the input string has
1157 already been obtained by the code in _Py_dg_strtod.
1159 The bigcomp function is only called after _Py_dg_strtod has found a
1160 double value rv such that either rv or rv + 1ulp represents the
1161 correctly rounded value corresponding to the original string. It
1162 determines which of these two values is the correct one by
1163 computing the decimal digits of rv + 0.5ulp and comparing them with
1164 the corresponding digits of s0.
1166 In the following, write dv for the absolute value of the number represented
1167 by the input string.
1171 s0 points to the first significant digit of the input string.
1173 rv is a (possibly scaled) estimate for the closest double value to the
1174 value represented by the original input to _Py_dg_strtod. If
1175 bc->scale is nonzero, then rv/2^(bc->scale) is the approximation to
1178 bc is a struct containing information gathered during the parsing and
1179 estimation steps of _Py_dg_strtod. Description of fields follows:
1181 bc->dsign is 1 if rv < decimal value, 0 if rv >= decimal value. In
1182 normal use, it should almost always be 1 when bigcomp is entered.
1184 bc->e0 gives the exponent of the input value, such that dv = (integer
1185 given by the bd->nd digits of s0) * 10**e0
1187 bc->nd gives the total number of significant digits of s0. It will
1190 bc->nd0 gives the number of significant digits of s0 before the
1191 decimal separator. If there's no decimal separator, bc->nd0 ==
1194 bc->scale is the value used to scale rv to avoid doing arithmetic with
1195 subnormal values. It's either 0 or 2*P (=106).
1199 On successful exit, rv/2^(bc->scale) is the closest double to dv.
1201 Returns 0 on success, -1 on failure (e.g., due to a failed malloc call). */
1204 bigcomp(U
*rv
, const char *s0
, BCinfo
*bc
)
1207 int b2
, bbits
, d2
, dd
, i
, nd
, nd0
, odd
, p2
, p5
;
1209 dd
= 0; /* silence compiler warning about possibly unused variable */
1214 /* special case because d2b doesn't handle 0.0 */
1218 p2
= Emin
- P
+ 1; /* = -1074 for IEEE 754 binary64 */
1222 b
= d2b(rv
, &p2
, &bbits
);
1227 /* now rv/2^(bc->scale) = b * 2**p2, and b has bbits significant bits */
1229 /* Replace (b, p2) by (b << i, p2 - i), with i the largest integer such
1230 that b << i has at most P significant bits and p2 - i >= Emin - P +
1233 if (i
> p2
- (Emin
- P
+ 1))
1234 i
= p2
- (Emin
- P
+ 1);
1235 /* increment i so that we shift b by an extra bit; then or-ing a 1 into
1236 the lsb of b gives us rv/2^(bc->scale) + 0.5ulp. */
1240 /* record whether the lsb of rv/2^(bc->scale) is odd: in the exact halfway
1241 case, this is used for round to even. */
1251 /* Arrange for convenient computation of quotients:
1252 * shift left if necessary so divisor has 4 leading 0 bits.
1255 d
= pow5mult(d
, p5
);
1262 b
= pow5mult(b
, -p5
);
1277 if ((b2
+= i
) > 0) {
1284 if ((d2
+= i
) > 0) {
1292 /* if b >= d, round down */
1293 if (cmp(b
, d
) >= 0) {
1298 /* Compare b/d with s0 */
1299 for(i
= 0; i
< nd0
; i
++) {
1300 b
= multadd(b
, 10, 0);
1305 dd
= *s0
++ - '0' - quorem(b
, d
);
1308 if (!b
->x
[0] && b
->wds
== 1) {
1315 for(; i
< nd
; i
++) {
1316 b
= multadd(b
, 10, 0);
1321 dd
= *s0
++ - '0' - quorem(b
, d
);
1324 if (!b
->x
[0] && b
->wds
== 1) {
1330 if (b
->x
[0] || b
->wds
> 1)
1335 if (dd
> 0 || (dd
== 0 && odd
))
1336 dval(rv
) += sulp(rv
, bc
);
1341 _Py_dg_strtod(const char *s00
, char **se
)
1343 int bb2
, bb5
, bbe
, bd2
, bd5
, bbbits
, bs2
, c
, e
, e1
, error
;
1344 int esign
, i
, j
, k
, nd
, nd0
, nf
, nz
, nz0
, sign
;
1345 const char *s
, *s0
, *s1
;
1347 U aadj2
, adj
, rv
, rv0
;
1351 Bigint
*bb
, *bb1
, *bd
, *bd0
, *bs
, *delta
;
1353 sign
= nz0
= nz
= 0;
1355 for(s
= s00
;;s
++) switch(*s
) {
1365 /* modify original dtoa.c so that it doesn't accept leading whitespace
1380 while(*++s
== '0') ;
1385 for(nd
= nf
= 0; (c
= *s
) >= '0' && c
<= '9'; nd
++, s
++)
1391 for(; c
== '0'; c
= *++s
)
1393 if (c
> '0' && c
<= '9') {
1401 for(; c
>= '0' && c
<= '9'; c
= *++s
) {
1413 if (c
== 'e' || c
== 'E') {
1414 if (!nd
&& !nz
&& !nz0
) {
1425 if (c
>= '0' && c
<= '9') {
1428 if (c
> '0' && c
<= '9') {
1431 while((c
= *++s
) >= '0' && c
<= '9')
1432 abse
= 10*abse
+ c
- '0';
1433 if (s
- s1
> 8 || abse
> MAX_ABS_EXP
)
1434 /* Avoid confusion from exponents
1435 * so large that e might overflow.
1437 e
= (int)MAX_ABS_EXP
; /* safe for 16 bit ints */
1461 /* strip trailing zeros */
1462 for (i
= nd
; i
> 0; ) {
1463 /* scan back until we hit a nonzero digit. significant digit 'i'
1464 is s0[i] if i < nd0, s0[i+1] if i >= nd0. */
1466 if (s0
[i
< nd0
? i
: i
+1] != '0') {
1476 /* Now we have nd0 digits, starting at s0, followed by a
1477 * decimal point, followed by nd-nd0 digits. The number we're
1478 * after is the integer represented by those digits times
1483 /* Summary of parsing results. The parsing stage gives values
1484 * s0, nd0, nd, e, sign, where:
1486 * - s0 points to the first significant digit of the input string s00;
1488 * - nd is the total number of significant digits (here, and
1489 * below, 'significant digits' means the set of digits of the
1490 * significand of the input that remain after ignoring leading
1491 * and trailing zeros.
1493 * - nd0 indicates the position of the decimal point (if
1494 * present): so the nd significant digits are in s0[0:nd0] and
1495 * s0[nd0+1:nd+1] using the usual Python half-open slice
1496 * notation. (If nd0 < nd, then s0[nd0] necessarily contains
1497 * a '.' character; if nd0 == nd, then it could be anything.)
1499 * - e is the adjusted exponent: the absolute value of the number
1500 * represented by the original input string is n * 10**e, where
1501 * n is the integer represented by the concatenation of
1502 * s0[0:nd0] and s0[nd0+1:nd+1]
1504 * - sign gives the sign of the input: 1 for negative, 0 for positive
1506 * - the first and last significant digits are nonzero
1509 /* put first DBL_DIG+1 digits into integer y and z.
1511 * - y contains the value represented by the first min(9, nd)
1512 * significant digits
1514 * - if nd > 9, z contains the value represented by significant digits
1515 * with indices in [9, min(16, nd)). So y * 10**(min(16, nd) - 9) + z
1516 * gives the value represented by the first min(16, nd) sig. digits.
1520 for (i
= 0; i
< nd
; i
++) {
1522 y
= 10*y
+ s0
[i
< nd0
? i
: i
+1] - '0';
1523 else if (i
< DBL_DIG
+1)
1524 z
= 10*z
+ s0
[i
< nd0
? i
: i
+1] - '0';
1529 k
= nd
< DBL_DIG
+ 1 ? nd
: DBL_DIG
+ 1;
1532 dval(&rv
) = tens
[k
- 9] * dval(&rv
) + z
;
1541 if (e
<= Ten_pmax
) {
1542 dval(&rv
) *= tens
[e
];
1546 if (e
<= Ten_pmax
+ i
) {
1547 /* A fancier test would sometimes let us do
1548 * this for larger i values.
1551 dval(&rv
) *= tens
[i
];
1552 dval(&rv
) *= tens
[e
];
1556 else if (e
>= -Ten_pmax
) {
1557 dval(&rv
) /= tens
[-e
];
1565 /* Get starting approximation = rv * 10**e1 */
1569 dval(&rv
) *= tens
[i
];
1571 if (e1
> DBL_MAX_10_EXP
) {
1574 /* Can't trust HUGE_VAL */
1575 word0(&rv
) = Exp_mask
;
1580 for(j
= 0; e1
> 1; j
++, e1
>>= 1)
1582 dval(&rv
) *= bigtens
[j
];
1583 /* The last multiplication could overflow. */
1584 word0(&rv
) -= P
*Exp_msk1
;
1585 dval(&rv
) *= bigtens
[j
];
1586 if ((z
= word0(&rv
) & Exp_mask
)
1587 > Exp_msk1
*(DBL_MAX_EXP
+Bias
-P
))
1589 if (z
> Exp_msk1
*(DBL_MAX_EXP
+Bias
-1-P
)) {
1590 /* set to largest number */
1591 /* (Can't trust DBL_MAX) */
1596 word0(&rv
) += P
*Exp_msk1
;
1602 dval(&rv
) /= tens
[i
];
1604 if (e1
>= 1 << n_bigtens
)
1608 for(j
= 0; e1
> 0; j
++, e1
>>= 1)
1610 dval(&rv
) *= tinytens
[j
];
1611 if (bc
.scale
&& (j
= 2*P
+ 1 - ((word0(&rv
) & Exp_mask
)
1612 >> Exp_shift
)) > 0) {
1613 /* scaled rv is denormal; clear j low bits */
1617 word0(&rv
) = (P
+2)*Exp_msk1
;
1619 word0(&rv
) &= 0xffffffff << (j
-32);
1622 word1(&rv
) &= 0xffffffff << j
;
1633 /* Now the hard part -- adjusting rv to the correct value.*/
1635 /* Put digits into bd: true value = bd * 10^e */
1638 bc
.nd0
= nd0
; /* Only needed if nd > STRTOD_DIGLIM, but done here */
1639 /* to silence an erroneous warning about bc.nd0 */
1640 /* possibly not being initialized. */
1641 if (nd
> STRTOD_DIGLIM
) {
1642 /* ASSERT(STRTOD_DIGLIM >= 18); 18 == one more than the */
1643 /* minimum number of decimal digits to distinguish double values */
1644 /* in IEEE arithmetic. */
1646 /* Truncate input to 18 significant digits, then discard any trailing
1647 zeros on the result by updating nd, nd0, e and y suitably. (There's
1648 no need to update z; it's not reused beyond this point.) */
1649 for (i
= 18; i
> 0; ) {
1650 /* scan back until we hit a nonzero digit. significant digit 'i'
1651 is s0[i] if i < nd0, s0[i+1] if i >= nd0. */
1653 if (s0
[i
< nd0
? i
: i
+1] != '0') {
1662 if (nd
< 9) { /* must recompute y */
1664 for(i
= 0; i
< nd0
; ++i
)
1665 y
= 10*y
+ s0
[i
] - '0';
1667 y
= 10*y
+ s0
[i
+1] - '0';
1670 bd0
= s2b(s0
, nd0
, nd
, y
);
1675 bd
= Balloc(bd0
->k
);
1681 bb
= d2b(&rv
, &bbe
, &bbbits
); /* rv = bb * 2^bbe */
1709 i
= j
+ bbbits
- 1; /* logb(rv) */
1710 if (i
< Emin
) /* denormal */
1717 i
= bb2
< bd2
? bb2
: bd2
;
1726 bs
= pow5mult(bs
, bb5
);
1744 bb
= lshift(bb
, bb2
);
1753 bd
= pow5mult(bd
, bd5
);
1762 bd
= lshift(bd
, bd2
);
1771 bs
= lshift(bs
, bs2
);
1779 delta
= diff(bb
, bd
);
1780 if (delta
== NULL
) {
1787 bc
.dsign
= delta
->sign
;
1790 if (bc
.nd
> nd
&& i
<= 0) {
1792 break; /* Must use bigcomp(). */
1794 /* Here rv overestimates the truncated decimal value by at most
1795 0.5 ulp(rv). Hence rv either overestimates the true decimal
1796 value by <= 0.5 ulp(rv), or underestimates it by some small
1797 amount (< 0.1 ulp(rv)); either way, rv is within 0.5 ulps of
1798 the true decimal value, so it's possible to exit.
1800 Exception: if scaled rv is a normal exact power of 2, but not
1801 DBL_MIN, then rv - 0.5 ulp(rv) takes us all the way down to the
1802 next double, so the correctly rounded result is either rv - 0.5
1803 ulp(rv) or rv; in this case, use bigcomp to distinguish. */
1805 if (!word1(&rv
) && !(word0(&rv
) & Bndry_mask
)) {
1806 /* rv can't be 0, since it's an overestimate for some
1807 nonzero value. So rv is a normal power of 2. */
1808 j
= (int)(word0(&rv
) & Exp_mask
) >> Exp_shift
;
1809 /* rv / 2^bc.scale = 2^(j - 1023 - bc.scale); use bigcomp if
1810 rv / 2^bc.scale >= 2^-1021. */
1811 if (j
- bc
.scale
>= 2) {
1812 dval(&rv
) -= 0.5 * sulp(&rv
, &bc
);
1819 i
= -1; /* Discarded digits make delta smaller. */
1824 /* Error is less than half an ulp -- check for
1825 * special case of mantissa a power of two.
1827 if (bc
.dsign
|| word1(&rv
) || word0(&rv
) & Bndry_mask
1828 || (word0(&rv
) & Exp_mask
) <= (2*P
+1)*Exp_msk1
1832 if (!delta
->x
[0] && delta
->wds
<= 1) {
1836 delta
= lshift(delta
,Log2P
);
1837 if (delta
== NULL
) {
1844 if (cmp(delta
, bs
) > 0)
1849 /* exactly half-way between */
1851 if ((word0(&rv
) & Bndry_mask1
) == Bndry_mask1
1854 (y
= word0(&rv
) & Exp_mask
) <= 2*P
*Exp_msk1
) ?
1855 (0xffffffff & (0xffffffff << (2*P
+1-(y
>>Exp_shift
)))) :
1857 /*boundary case -- increment exponent*/
1858 word0(&rv
) = (word0(&rv
) & Exp_mask
)
1866 else if (!(word0(&rv
) & Bndry_mask
) && !word1(&rv
)) {
1868 /* boundary case -- decrement exponent */
1870 L
= word0(&rv
) & Exp_mask
;
1871 if (L
<= (2*P
+1)*Exp_msk1
) {
1872 if (L
> (P
+2)*Exp_msk1
)
1873 /* round even ==> */
1876 /* rv = smallest denormal */
1882 L
= (word0(&rv
) & Exp_mask
) - Exp_msk1
;
1883 word0(&rv
) = L
| Bndry_mask1
;
1884 word1(&rv
) = 0xffffffff;
1887 if (!(word1(&rv
) & LSB
))
1890 dval(&rv
) += ulp(&rv
);
1892 dval(&rv
) -= ulp(&rv
);
1899 bc
.dsign
= 1 - bc
.dsign
;
1902 if ((aadj
= ratio(delta
, bs
)) <= 2.) {
1905 else if (word1(&rv
) || word0(&rv
) & Bndry_mask
) {
1906 if (word1(&rv
) == Tiny1
&& !word0(&rv
)) {
1915 /* special case -- power of FLT_RADIX to be */
1916 /* rounded down... */
1918 if (aadj
< 2./FLT_RADIX
)
1919 aadj
= 1./FLT_RADIX
;
1927 aadj1
= bc
.dsign
? aadj
: -aadj
;
1928 if (Flt_Rounds
== 0)
1931 y
= word0(&rv
) & Exp_mask
;
1933 /* Check for overflow */
1935 if (y
== Exp_msk1
*(DBL_MAX_EXP
+Bias
-1)) {
1936 dval(&rv0
) = dval(&rv
);
1937 word0(&rv
) -= P
*Exp_msk1
;
1938 adj
.d
= aadj1
* ulp(&rv
);
1940 if ((word0(&rv
) & Exp_mask
) >=
1941 Exp_msk1
*(DBL_MAX_EXP
+Bias
-P
)) {
1942 if (word0(&rv0
) == Big0
&& word1(&rv0
) == Big1
)
1949 word0(&rv
) += P
*Exp_msk1
;
1952 if (bc
.scale
&& y
<= 2*P
*Exp_msk1
) {
1953 if (aadj
<= 0x7fffffff) {
1954 if ((z
= (ULong
)aadj
) <= 0)
1957 aadj1
= bc
.dsign
? aadj
: -aadj
;
1959 dval(&aadj2
) = aadj1
;
1960 word0(&aadj2
) += (2*P
+1)*Exp_msk1
- y
;
1961 aadj1
= dval(&aadj2
);
1963 adj
.d
= aadj1
* ulp(&rv
);
1966 z
= word0(&rv
) & Exp_mask
;
1970 /* Can we stop now? */
1973 /* The tolerances below are conservative. */
1974 if (bc
.dsign
|| word1(&rv
) || word0(&rv
) & Bndry_mask
) {
1975 if (aadj
< .4999999 || aadj
> .5000001)
1978 else if (aadj
< .4999999/FLT_RADIX
)
1994 error
= bigcomp(&rv
, s0
, &bc
);
2000 word0(&rv0
) = Exp_1
- 2*P
*Exp_msk1
;
2002 dval(&rv
) *= dval(&rv0
);
2003 /* try to avoid the bug of testing an 8087 register value */
2004 if (!(word0(&rv
) & Exp_mask
))
2010 return sign
? -dval(&rv
) : dval(&rv
);
2026 sizeof(Bigint
) - sizeof(ULong
) - sizeof(int) + j
<= (unsigned)i
;
2029 r
= (int*)Balloc(k
);
2033 return (char *)(r
+1);
2037 nrv_alloc(char *s
, char **rve
, int n
)
2045 while((*t
= *s
++)) t
++;
2051 /* freedtoa(s) must be used to free values s returned by dtoa
2052 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
2053 * but for consistency with earlier versions of dtoa, it is optional
2054 * when MULTIPLE_THREADS is not defined.
2058 _Py_dg_freedtoa(char *s
)
2060 Bigint
*b
= (Bigint
*)((int *)s
- 1);
2061 b
->maxwds
= 1 << (b
->k
= *(int*)b
);
2065 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
2067 * Inspired by "How to Print Floating-Point Numbers Accurately" by
2068 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
2071 * 1. Rather than iterating, we use a simple numeric overestimate
2072 * to determine k = floor(log10(d)). We scale relevant
2073 * quantities using O(log2(k)) rather than O(k) multiplications.
2074 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
2075 * try to generate digits strictly left to right. Instead, we
2076 * compute with fewer bits and propagate the carry if necessary
2077 * when rounding the final digit up. This is often faster.
2078 * 3. Under the assumption that input will be rounded nearest,
2079 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
2080 * That is, we allow equality in stopping tests when the
2081 * round-nearest rule will give the same floating-point value
2082 * as would satisfaction of the stopping test with strict
2084 * 4. We remove common factors of powers of 2 from relevant
2086 * 5. When converting floating-point integers less than 1e16,
2087 * we use floating-point arithmetic rather than resorting
2088 * to multiple-precision integers.
2089 * 6. When asked to produce fewer than 15 digits, we first try
2090 * to get by with floating-point arithmetic; we resort to
2091 * multiple-precision integer arithmetic only if we cannot
2092 * guarantee that the floating-point calculation has given
2093 * the correctly rounded result. For k requested digits and
2094 * "uniformly" distributed input, the probability is
2095 * something like 10^(k-15) that we must resort to the Long
2099 /* Additional notes (METD): (1) returns NULL on failure. (2) to avoid memory
2100 leakage, a successful call to _Py_dg_dtoa should always be matched by a
2101 call to _Py_dg_freedtoa. */
2104 _Py_dg_dtoa(double dd
, int mode
, int ndigits
,
2105 int *decpt
, int *sign
, char **rve
)
2107 /* Arguments ndigits, decpt, sign are similar to those
2108 of ecvt and fcvt; trailing zeros are suppressed from
2109 the returned string. If not null, *rve is set to point
2110 to the end of the return value. If d is +-Infinity or NaN,
2111 then *decpt is set to 9999.
2114 0 ==> shortest string that yields d when read in
2115 and rounded to nearest.
2116 1 ==> like 0, but with Steele & White stopping rule;
2117 e.g. with IEEE P754 arithmetic , mode 0 gives
2118 1e23 whereas mode 1 gives 9.999999999999999e22.
2119 2 ==> max(1,ndigits) significant digits. This gives a
2120 return value similar to that of ecvt, except
2121 that trailing zeros are suppressed.
2122 3 ==> through ndigits past the decimal point. This
2123 gives a return value similar to that from fcvt,
2124 except that trailing zeros are suppressed, and
2125 ndigits can be negative.
2126 4,5 ==> similar to 2 and 3, respectively, but (in
2127 round-nearest mode) with the tests of mode 0 to
2128 possibly return a shorter string that rounds to d.
2129 With IEEE arithmetic and compilation with
2130 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
2131 as modes 2 and 3 when FLT_ROUNDS != 1.
2132 6-9 ==> Debugging modes similar to mode - 4: don't try
2133 fast floating-point estimate (if applicable).
2135 Values of mode other than 0-9 are treated as mode 0.
2137 Sufficient space is allocated to the return value
2138 to hold the suppressed trailing zeros.
2141 int bbits
, b2
, b5
, be
, dig
, i
, ieps
, ilim
, ilim0
, ilim1
,
2142 j
, j1
, k
, k0
, k_check
, leftright
, m2
, m5
, s2
, s5
,
2143 spec_case
, try_quick
;
2147 Bigint
*b
, *b1
, *delta
, *mlo
, *mhi
, *S
;
2152 /* set pointers to NULL, to silence gcc compiler warnings and make
2153 cleanup easier on error */
2154 mlo
= mhi
= b
= S
= 0;
2158 if (word0(&u
) & Sign_bit
) {
2159 /* set sign for everything, including 0's and NaNs */
2161 word0(&u
) &= ~Sign_bit
; /* clear sign bit */
2166 /* quick return for Infinities, NaNs and zeros */
2167 if ((word0(&u
) & Exp_mask
) == Exp_mask
)
2169 /* Infinity or NaN */
2171 if (!word1(&u
) && !(word0(&u
) & 0xfffff))
2172 return nrv_alloc("Infinity", rve
, 8);
2173 return nrv_alloc("NaN", rve
, 3);
2177 return nrv_alloc("0", rve
, 1);
2180 /* compute k = floor(log10(d)). The computation may leave k
2181 one too large, but should never leave k too small. */
2182 b
= d2b(&u
, &be
, &bbits
);
2185 if ((i
= (int)(word0(&u
) >> Exp_shift1
& (Exp_mask
>>Exp_shift1
)))) {
2186 dval(&d2
) = dval(&u
);
2187 word0(&d2
) &= Frac_mask1
;
2188 word0(&d2
) |= Exp_11
;
2190 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2191 * log10(x) = log(x) / log(10)
2192 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2193 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2195 * This suggests computing an approximation k to log10(d) by
2197 * k = (i - Bias)*0.301029995663981
2198 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2200 * We want k to be too large rather than too small.
2201 * The error in the first-order Taylor series approximation
2202 * is in our favor, so we just round up the constant enough
2203 * to compensate for any error in the multiplication of
2204 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2205 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2206 * adding 1e-13 to the constant term more than suffices.
2207 * Hence we adjust the constant term to 0.1760912590558.
2208 * (We could get a more accurate k by invoking log10,
2209 * but this is probably not worthwhile.)
2216 /* d is denormalized */
2218 i
= bbits
+ be
+ (Bias
+ (P
-1) - 1);
2219 x
= i
> 32 ? word0(&u
) << (64 - i
) | word1(&u
) >> (i
- 32)
2220 : word1(&u
) << (32 - i
);
2222 word0(&d2
) -= 31*Exp_msk1
; /* adjust exponent */
2223 i
-= (Bias
+ (P
-1) - 1) + 1;
2226 ds
= (dval(&d2
)-1.5)*0.289529654602168 + 0.1760912590558 +
2227 i
*0.301029995663981;
2229 if (ds
< 0. && ds
!= k
)
2230 k
--; /* want k = floor(ds) */
2232 if (k
>= 0 && k
<= Ten_pmax
) {
2233 if (dval(&u
) < tens
[k
])
2256 if (mode
< 0 || mode
> 9)
2266 ilim
= ilim1
= -1; /* Values for cases 0 and 1; done here to */
2267 /* silence erroneous "gcc -Wall" warning. */
2280 ilim
= ilim1
= i
= ndigits
;
2286 i
= ndigits
+ k
+ 1;
2298 if (ilim
>= 0 && ilim
<= Quick_max
&& try_quick
) {
2300 /* Try to get by with floating-point arithmetic. */
2303 dval(&d2
) = dval(&u
);
2306 ieps
= 2; /* conservative */
2311 /* prevent overflows */
2313 dval(&u
) /= bigtens
[n_bigtens
-1];
2316 for(; j
; j
>>= 1, i
++)
2323 else if ((j1
= -k
)) {
2324 dval(&u
) *= tens
[j1
& 0xf];
2325 for(j
= j1
>> 4; j
; j
>>= 1, i
++)
2328 dval(&u
) *= bigtens
[i
];
2331 if (k_check
&& dval(&u
) < 1. && ilim
> 0) {
2339 dval(&eps
) = ieps
*dval(&u
) + 7.;
2340 word0(&eps
) -= (P
-1)*Exp_msk1
;
2344 if (dval(&u
) > dval(&eps
))
2346 if (dval(&u
) < -dval(&eps
))
2351 /* Use Steele & White method of only
2352 * generating digits needed.
2354 dval(&eps
) = 0.5/tens
[ilim
-1] - dval(&eps
);
2358 *s
++ = '0' + (int)L
;
2359 if (dval(&u
) < dval(&eps
))
2361 if (1. - dval(&u
) < dval(&eps
))
2370 /* Generate ilim digits, then fix them up. */
2371 dval(&eps
) *= tens
[ilim
-1];
2372 for(i
= 1;; i
++, dval(&u
) *= 10.) {
2373 L
= (Long
)(dval(&u
));
2374 if (!(dval(&u
) -= L
))
2376 *s
++ = '0' + (int)L
;
2378 if (dval(&u
) > 0.5 + dval(&eps
))
2380 else if (dval(&u
) < 0.5 - dval(&eps
)) {
2391 dval(&u
) = dval(&d2
);
2396 /* Do we have a "small" integer? */
2398 if (be
>= 0 && k
<= Int_max
) {
2401 if (ndigits
< 0 && ilim
<= 0) {
2403 if (ilim
< 0 || dval(&u
) <= 5*ds
)
2407 for(i
= 1;; i
++, dval(&u
) *= 10.) {
2408 L
= (Long
)(dval(&u
) / ds
);
2410 *s
++ = '0' + (int)L
;
2415 dval(&u
) += dval(&u
);
2416 if (dval(&u
) > ds
|| (dval(&u
) == ds
&& L
& 1)) {
2436 denorm
? be
+ (Bias
+ (P
-1) - 1 + 1) :
2444 if (m2
> 0 && s2
> 0) {
2445 i
= m2
< s2
? m2
: s2
;
2453 mhi
= pow5mult(mhi
, m5
);
2462 if ((j
= b5
- m5
)) {
2469 b
= pow5mult(b
, b5
);
2478 S
= pow5mult(S
, s5
);
2483 /* Check for special case that d is a normalized power of 2. */
2486 if ((mode
< 2 || leftright
)
2488 if (!word1(&u
) && !(word0(&u
) & Bndry_mask
)
2489 && word0(&u
) & (Exp_mask
& ~Exp_msk1
)
2491 /* The special case */
2498 /* Arrange for convenient computation of quotients:
2499 * shift left if necessary so divisor has 4 leading 0 bits.
2501 * Perhaps we should just compute leading 28 bits of S once
2502 * and for all and pass them and a shift to quorem, so it
2503 * can do shifts and ors to compute the numerator for q.
2505 if ((i
= ((s5
? 32 - hi0bits(S
->x
[S
->wds
-1]) : 1) + s2
) & 0x1f))
2525 b
= multadd(b
, 10, 0); /* we botched the k estimate */
2529 mhi
= multadd(mhi
, 10, 0);
2536 if (ilim
<= 0 && (mode
== 3 || mode
== 5)) {
2538 /* no digits, fcvt style */
2544 S
= multadd(S
, 5, 0);
2557 mhi
= lshift(mhi
, m2
);
2562 /* Compute mlo -- check for special case
2563 * that d is a normalized power of 2.
2568 mhi
= Balloc(mhi
->k
);
2572 mhi
= lshift(mhi
, Log2P
);
2578 dig
= quorem(b
,S
) + '0';
2579 /* Do we yet have the shortest decimal string
2580 * that will round to d?
2583 delta
= diff(S
, mhi
);
2586 j1
= delta
->sign
? 1 : cmp(b
, delta
);
2588 if (j1
== 0 && mode
!= 1 && !(word1(&u
) & 1)
2597 if (j
< 0 || (j
== 0 && mode
!= 1
2600 if (!b
->x
[0] && b
->wds
<= 1) {
2608 if ((j1
> 0 || (j1
== 0 && dig
& 1))
2617 if (dig
== '9') { /* possible if i == 1 */
2628 b
= multadd(b
, 10, 0);
2632 mlo
= mhi
= multadd(mhi
, 10, 0);
2637 mlo
= multadd(mlo
, 10, 0);
2640 mhi
= multadd(mhi
, 10, 0);
2648 *s
++ = dig
= quorem(b
,S
) + '0';
2649 if (!b
->x
[0] && b
->wds
<= 1) {
2654 b
= multadd(b
, 10, 0);
2659 /* Round off last digit */
2665 if (j
> 0 || (j
== 0 && dig
& 1)) {
2682 if (mlo
&& mlo
!= mhi
)
2696 if (mlo
&& mlo
!= mhi
)
2703 _Py_dg_freedtoa(s0
);
2710 #endif /* PY_NO_SHORT_FLOAT_REPR */