]>
git.proxmox.com Git - mirror_edk2.git/blob - AppPkg/Applications/Python/Python-2.7.10/Python/dtoa.c
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 1100000000U
209 /* Bound on length of pieces of input strings in _Py_dg_strtod; specifically,
210 this is used to bound the total number of digits ignoring leading zeros and
211 the number of digits that follow the decimal point. Ideally, MAX_DIGITS
212 should satisfy MAX_DIGITS + 400 < MAX_ABS_EXP; that ensures that the
213 exponent clipping in _Py_dg_strtod can't affect the value of the output. */
215 #define MAX_DIGITS 1000000000U
218 /* Guard against trying to use the above values on unusual platforms with ints
219 * of width less than 32 bits. */
220 #if MAX_ABS_EXP > INT_MAX
221 #error "MAX_ABS_EXP should fit in an int"
223 #if MAX_DIGITS > INT_MAX
224 #error "MAX_DIGITS should fit in an int"
227 /* The following definition of Storeinc is appropriate for MIPS processors.
228 * An alternative that might be better on some machines is
229 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
231 #if defined(IEEE_8087)
232 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
233 ((unsigned short *)a)[0] = (unsigned short)c, a++)
235 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
236 ((unsigned short *)a)[1] = (unsigned short)c, a++)
239 /* #define P DBL_MANT_DIG */
240 /* Ten_pmax = floor(P*log(2)/log(5)) */
241 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
242 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
243 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
246 #define Exp_shift1 20
247 #define Exp_msk1 0x100000
248 #define Exp_msk11 0x100000
249 #define Exp_mask 0x7ff00000
255 #define Etiny (-1074) /* smallest denormal is 2**Etiny */
256 #define Exp_1 0x3ff00000
257 #define Exp_11 0x3ff00000
259 #define Frac_mask 0xfffff
260 #define Frac_mask1 0xfffff
263 #define Bndry_mask 0xfffff
264 #define Bndry_mask1 0xfffff
265 #define Sign_bit 0x80000000
274 #define Flt_Rounds FLT_ROUNDS
278 #endif /*Flt_Rounds*/
280 #define Rounding Flt_Rounds
282 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
283 #define Big1 0xffffffff
285 /* struct BCinfo is used to pass information from _Py_dg_strtod to bigcomp */
287 typedef struct BCinfo BCinfo
;
290 int e0
, nd
, nd0
, scale
;
293 #define FFFFFFFF 0xffffffffUL
297 /* struct Bigint is used to represent arbitrary-precision integers. These
298 integers are stored in sign-magnitude format, with the magnitude stored as
299 an array of base 2**32 digits. Bigints are always normalized: if x is a
300 Bigint then x->wds >= 1, and either x->wds == 1 or x[wds-1] is nonzero.
302 The Bigint fields are as follows:
304 - next is a header used by Balloc and Bfree to keep track of lists
305 of freed Bigints; it's also used for the linked list of
306 powers of 5 of the form 5**2**i used by pow5mult.
307 - k indicates which pool this Bigint was allocated from
308 - maxwds is the maximum number of words space was allocated for
309 (usually maxwds == 2**k)
310 - sign is 1 for negative Bigints, 0 for positive. The sign is unused
311 (ignored on inputs, set to 0 on outputs) in almost all operations
312 involving Bigints: a notable exception is the diff function, which
313 ignores signs on inputs but sets the sign of the output correctly.
314 - wds is the actual number of significant words
315 - x contains the vector of words (digits) for this Bigint, from least
316 significant (x[0]) to most significant (x[wds-1]).
322 int k
, maxwds
, sign
, wds
;
326 typedef struct Bigint Bigint
;
328 #ifndef Py_USING_MEMORY_DEBUGGER
330 /* Memory management: memory is allocated from, and returned to, Kmax+1 pools
331 of memory, where pool k (0 <= k <= Kmax) is for Bigints b with b->maxwds ==
332 1 << k. These pools are maintained as linked lists, with freelist[k]
333 pointing to the head of the list for pool k.
335 On allocation, if there's no free slot in the appropriate pool, MALLOC is
336 called to get more memory. This memory is not returned to the system until
337 Python quits. There's also a private memory pool that's allocated from
338 in preference to using MALLOC.
340 For Bigints with more than (1 << Kmax) digits (which implies at least 1233
341 decimal digits), memory is directly allocated using MALLOC, and freed using
344 XXX: it would be easy to bypass this memory-management system and
345 translate each call to Balloc into a call to PyMem_Malloc, and each
346 Bfree to PyMem_Free. Investigate whether this has any significant
347 performance on impact. */
349 static Bigint
*freelist
[Kmax
+1];
351 /* Allocate space for a Bigint with up to 1<<k digits */
360 if (k
<= Kmax
&& (rv
= freelist
[k
]))
361 freelist
[k
] = rv
->next
;
364 len
= (sizeof(Bigint
) + (x
-1)*sizeof(ULong
) + sizeof(double) - 1)
366 if (k
<= Kmax
&& pmem_next
- private_mem
+ len
<= PRIVATE_mem
) {
367 rv
= (Bigint
*)pmem_next
;
371 rv
= (Bigint
*)MALLOC(len
*sizeof(double));
378 rv
->sign
= rv
->wds
= 0;
382 /* Free a Bigint allocated with Balloc */
391 v
->next
= freelist
[v
->k
];
399 /* Alternative versions of Balloc and Bfree that use PyMem_Malloc and
400 PyMem_Free directly in place of the custom memory allocation scheme above.
401 These are provided for the benefit of memory debugging tools like
404 /* Allocate space for a Bigint with up to 1<<k digits */
414 len
= (sizeof(Bigint
) + (x
-1)*sizeof(ULong
) + sizeof(double) - 1)
417 rv
= (Bigint
*)MALLOC(len
*sizeof(double));
423 rv
->sign
= rv
->wds
= 0;
427 /* Free a Bigint allocated with Balloc */
437 #endif /* Py_USING_MEMORY_DEBUGGER */
439 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
440 y->wds*sizeof(Long) + 2*sizeof(int))
442 /* Multiply a Bigint b by m and add a. Either modifies b in place and returns
443 a pointer to the modified b, or Bfrees b and returns a pointer to a copy.
444 On failure, return NULL. In this case, b will have been already freed. */
447 multadd(Bigint
*b
, int m
, int a
) /* multiply by m and add a */
465 y
= *x
* (ULLong
)m
+ carry
;
467 *x
++ = (ULong
)(y
& FFFFFFFF
);
470 y
= (xi
& 0xffff) * m
+ carry
;
471 z
= (xi
>> 16) * m
+ (y
>> 16);
473 *x
++ = (z
<< 16) + (y
& 0xffff);
478 if (wds
>= b
->maxwds
) {
488 b
->x
[wds
++] = (ULong
)carry
;
494 /* convert a string s containing nd decimal digits (possibly containing a
495 decimal separator at position nd0, which is ignored) to a Bigint. This
496 function carries on where the parsing code in _Py_dg_strtod leaves off: on
497 entry, y9 contains the result of converting the first 9 digits. Returns
501 s2b(const char *s
, int nd0
, int nd
, ULong y9
)
508 for(k
= 0, y
= 1; x
> y
; y
<<= 1, k
++) ;
519 for (i
= 9; i
< nd0
; i
++) {
520 b
= multadd(b
, 10, *s
++ - '0');
526 b
= multadd(b
, 10, *s
++ - '0');
533 /* count leading 0 bits in the 32-bit integer x. */
540 if (!(x
& 0xffff0000)) {
544 if (!(x
& 0xff000000)) {
548 if (!(x
& 0xf0000000)) {
552 if (!(x
& 0xc0000000)) {
556 if (!(x
& 0x80000000)) {
558 if (!(x
& 0x40000000))
564 /* count trailing 0 bits in the 32-bit integer y, and shift y right by that
610 /* convert a small nonnegative integer to a Bigint */
625 /* multiply two Bigints. Returns a new Bigint, or NULL on failure. Ignores
626 the signs of a and b. */
629 mult(Bigint
*a
, Bigint
*b
)
633 ULong
*x
, *xa
, *xae
, *xb
, *xbe
, *xc
, *xc0
;
642 if ((!a
->x
[0] && a
->wds
== 1) || (!b
->x
[0] && b
->wds
== 1)) {
651 if (a
->wds
< b
->wds
) {
665 for(x
= c
->x
, xa
= x
+ wc
; x
< xa
; x
++)
673 for(; xb
< xbe
; xc0
++) {
679 z
= *x
++ * (ULLong
)y
+ *xc
+ carry
;
681 *xc
++ = (ULong
)(z
& FFFFFFFF
);
688 for(; xb
< xbe
; xb
++, xc0
++) {
689 if (y
= *xb
& 0xffff) {
694 z
= (*x
& 0xffff) * y
+ (*xc
& 0xffff) + carry
;
696 z2
= (*x
++ >> 16) * y
+ (*xc
>> 16) + carry
;
709 z
= (*x
& 0xffff) * y
+ (*xc
>> 16) + carry
;
712 z2
= (*x
++ >> 16) * y
+ (*xc
& 0xffff) + carry
;
720 for(xc0
= c
->x
, xc
= xc0
+ wc
; wc
> 0 && !*--xc
; --wc
) ;
725 #ifndef Py_USING_MEMORY_DEBUGGER
727 /* p5s is a linked list of powers of 5 of the form 5**(2**i), i >= 2 */
731 /* multiply the Bigint b by 5**k. Returns a pointer to the result, or NULL on
732 failure; if the returned pointer is distinct from b then the original
733 Bigint b will have been Bfree'd. Ignores the sign of b. */
736 pow5mult(Bigint
*b
, int k
)
738 Bigint
*b1
, *p5
, *p51
;
740 static int p05
[3] = { 5, 25, 125 };
743 b
= multadd(b
, p05
[i
-1], 0);
788 /* Version of pow5mult that doesn't cache powers of 5. Provided for
789 the benefit of memory debugging tools like Valgrind. */
792 pow5mult(Bigint
*b
, int k
)
794 Bigint
*b1
, *p5
, *p51
;
796 static int p05
[3] = { 5, 25, 125 };
799 b
= multadd(b
, p05
[i
-1], 0);
836 #endif /* Py_USING_MEMORY_DEBUGGER */
838 /* shift a Bigint b left by k bits. Return a pointer to the shifted result,
839 or NULL on failure. If the returned pointer is distinct from b then the
840 original b will have been Bfree'd. Ignores the sign of b. */
843 lshift(Bigint
*b
, int k
)
847 ULong
*x
, *x1
, *xe
, z
;
849 if (!k
|| (!b
->x
[0] && b
->wds
== 1))
855 for(i
= b
->maxwds
; n1
> i
; i
<<= 1)
863 for(i
= 0; i
< n
; i
++)
886 /* Do a three-way compare of a and b, returning -1 if a < b, 0 if a == b and
887 1 if a > b. Ignores signs of a and b. */
890 cmp(Bigint
*a
, Bigint
*b
)
892 ULong
*xa
, *xa0
, *xb
, *xb0
;
898 if (i
> 1 && !a
->x
[i
-1])
899 Bug("cmp called with a->x[a->wds-1] == 0");
900 if (j
> 1 && !b
->x
[j
-1])
901 Bug("cmp called with b->x[b->wds-1] == 0");
911 return *xa
< *xb
? -1 : 1;
918 /* Take the difference of Bigints a and b, returning a new Bigint. Returns
919 NULL on failure. The signs of a and b are ignored, but the sign of the
920 result is set appropriately. */
923 diff(Bigint
*a
, Bigint
*b
)
927 ULong
*xa
, *xae
, *xb
, *xbe
, *xc
;
966 y
= (ULLong
)*xa
++ - *xb
++ - borrow
;
967 borrow
= y
>> 32 & (ULong
)1;
968 *xc
++ = (ULong
)(y
& FFFFFFFF
);
973 borrow
= y
>> 32 & (ULong
)1;
974 *xc
++ = (ULong
)(y
& FFFFFFFF
);
978 y
= (*xa
& 0xffff) - (*xb
& 0xffff) - borrow
;
979 borrow
= (y
& 0x10000) >> 16;
980 z
= (*xa
++ >> 16) - (*xb
++ >> 16) - borrow
;
981 borrow
= (z
& 0x10000) >> 16;
986 y
= (*xa
& 0xffff) - borrow
;
987 borrow
= (y
& 0x10000) >> 16;
988 z
= (*xa
++ >> 16) - borrow
;
989 borrow
= (z
& 0x10000) >> 16;
999 /* Given a positive normal double x, return the difference between x and the
1000 next double up. Doesn't give correct results for subnormals. */
1008 L
= (word0(x
) & Exp_mask
) - (P
-1)*Exp_msk1
;
1014 /* Convert a Bigint to a double plus an exponent */
1017 b2d(Bigint
*a
, int *e
)
1019 ULong
*xa
, *xa0
, w
, y
, z
;
1027 if (!y
) Bug("zero y in b2d");
1032 word0(&d
) = Exp_1
| y
>> (Ebits
- k
);
1033 w
= xa
> xa0
? *--xa
: 0;
1034 word1(&d
) = y
<< ((32-Ebits
) + k
) | w
>> (Ebits
- k
);
1037 z
= xa
> xa0
? *--xa
: 0;
1039 word0(&d
) = Exp_1
| y
<< k
| z
>> (32 - k
);
1040 y
= xa
> xa0
? *--xa
: 0;
1041 word1(&d
) = z
<< k
| y
>> (32 - k
);
1044 word0(&d
) = Exp_1
| y
;
1051 /* Convert a scaled double to a Bigint plus an exponent. Similar to d2b,
1052 except that it accepts the scale parameter used in _Py_dg_strtod (which
1053 should be either 0 or 2*P), and the normalization for the return value is
1054 different (see below). On input, d should be finite and nonnegative, and d
1055 / 2**scale should be exactly representable as an IEEE 754 double.
1057 Returns a Bigint b and an integer e such that
1059 dval(d) / 2**scale = b * 2**e.
1061 Unlike d2b, b is not necessarily odd: b and e are normalized so
1062 that either 2**(P-1) <= b < 2**P and e >= Etiny, or b < 2**P
1063 and e == Etiny. This applies equally to an input of 0.0: in that
1064 case the return values are b = 0 and e = Etiny.
1066 The above normalization ensures that for all possible inputs d,
1067 2**e gives ulp(d/2**scale).
1069 Returns NULL on failure.
1073 sd2b(U
*d
, int scale
, int *e
)
1081 /* First construct b and e assuming that scale == 0. */
1084 b
->x
[1] = word0(d
) & Frac_mask
;
1085 *e
= Etiny
- 1 + (int)((word0(d
) & Exp_mask
) >> Exp_shift
);
1089 b
->x
[1] |= Exp_msk1
;
1091 /* Now adjust for scale, provided that b != 0. */
1092 if (scale
&& (b
->x
[0] || b
->x
[1])) {
1097 /* We can't shift more than P-1 bits without shifting out a 1. */
1098 assert(0 < scale
&& scale
<= P
- 1);
1100 /* The bits shifted out should all be zero. */
1101 assert(b
->x
[0] == 0);
1107 /* The bits shifted out should all be zero. */
1108 assert(b
->x
[0] << (32 - scale
) == 0);
1109 b
->x
[0] = (b
->x
[0] >> scale
) | (b
->x
[1] << (32 - scale
));
1114 /* Ensure b is normalized. */
1121 /* Convert a double to a Bigint plus an exponent. Return NULL on failure.
1123 Given a finite nonzero double d, return an odd Bigint b and exponent *e
1124 such that fabs(d) = b * 2**e. On return, *bbits gives the number of
1125 significant bits of b; that is, 2**(*bbits-1) <= b < 2**(*bbits).
1127 If d is zero, then b == 0, *e == -1010, *bbits = 0.
1131 d2b(U
*d
, int *e
, int *bits
)
1143 z
= word0(d
) & Frac_mask
;
1144 word0(d
) &= 0x7fffffff; /* clear sign bit, which we ignore */
1145 if ((de
= (int)(word0(d
) >> Exp_shift
)))
1147 if ((y
= word1(d
))) {
1148 if ((k
= lo0bits(&y
))) {
1149 x
[0] = y
| z
<< (32 - k
);
1155 b
->wds
= (x
[1] = z
) ? 2 : 1;
1165 *e
= de
- Bias
- (P
-1) + k
;
1169 *e
= de
- Bias
- (P
-1) + 1 + k
;
1170 *bits
= 32*i
- hi0bits(x
[i
-1]);
1175 /* Compute the ratio of two Bigints, as a double. The result may have an
1176 error of up to 2.5 ulps. */
1179 ratio(Bigint
*a
, Bigint
*b
)
1184 dval(&da
) = b2d(a
, &ka
);
1185 dval(&db
) = b2d(b
, &kb
);
1186 k
= ka
- kb
+ 32*(a
->wds
- b
->wds
);
1188 word0(&da
) += k
*Exp_msk1
;
1191 word0(&db
) += k
*Exp_msk1
;
1193 return dval(&da
) / dval(&db
);
1198 1e0
, 1e1
, 1e2
, 1e3
, 1e4
, 1e5
, 1e6
, 1e7
, 1e8
, 1e9
,
1199 1e10
, 1e11
, 1e12
, 1e13
, 1e14
, 1e15
, 1e16
, 1e17
, 1e18
, 1e19
,
1204 bigtens
[] = { 1e16
, 1e32
, 1e64
, 1e128
, 1e256
};
1205 static const double tinytens
[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1206 9007199254740992.*9007199254740992.e
-256
1207 /* = 2^106 * 1e-256 */
1209 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1210 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1211 #define Scale_Bit 0x10
1220 dshift(Bigint
*b
, int p2
)
1222 int rv
= hi0bits(b
->x
[b
->wds
-1]) - 4;
1228 /* special case of Bigint division. The quotient is always in the range 0 <=
1229 quotient < 10, and on entry the divisor S is normalized so that its top 4
1230 bits (28--31) are zero and bit 27 is set. */
1233 quorem(Bigint
*b
, Bigint
*S
)
1236 ULong
*bx
, *bxe
, q
, *sx
, *sxe
;
1238 ULLong borrow
, carry
, y
, ys
;
1240 ULong borrow
, carry
, y
, ys
;
1246 /*debug*/ if (b
->wds
> n
)
1247 /*debug*/ Bug("oversize b in quorem");
1255 q
= *bxe
/ (*sxe
+ 1); /* ensure q <= true quotient */
1257 /*debug*/ if (q
> 9)
1258 /*debug*/ Bug("oversized quotient in quorem");
1265 ys
= *sx
++ * (ULLong
)q
+ carry
;
1267 y
= *bx
- (ys
& FFFFFFFF
) - borrow
;
1268 borrow
= y
>> 32 & (ULong
)1;
1269 *bx
++ = (ULong
)(y
& FFFFFFFF
);
1272 ys
= (si
& 0xffff) * q
+ carry
;
1273 zs
= (si
>> 16) * q
+ (ys
>> 16);
1275 y
= (*bx
& 0xffff) - (ys
& 0xffff) - borrow
;
1276 borrow
= (y
& 0x10000) >> 16;
1277 z
= (*bx
>> 16) - (zs
& 0xffff) - borrow
;
1278 borrow
= (z
& 0x10000) >> 16;
1285 while(--bxe
> bx
&& !*bxe
)
1290 if (cmp(b
, S
) >= 0) {
1300 y
= *bx
- (ys
& FFFFFFFF
) - borrow
;
1301 borrow
= y
>> 32 & (ULong
)1;
1302 *bx
++ = (ULong
)(y
& FFFFFFFF
);
1305 ys
= (si
& 0xffff) + carry
;
1306 zs
= (si
>> 16) + (ys
>> 16);
1308 y
= (*bx
& 0xffff) - (ys
& 0xffff) - borrow
;
1309 borrow
= (y
& 0x10000) >> 16;
1310 z
= (*bx
>> 16) - (zs
& 0xffff) - borrow
;
1311 borrow
= (z
& 0x10000) >> 16;
1319 while(--bxe
> bx
&& !*bxe
)
1327 /* sulp(x) is a version of ulp(x) that takes bc.scale into account.
1329 Assuming that x is finite and nonnegative (positive zero is fine
1330 here) and x / 2^bc.scale is exactly representable as a double,
1331 sulp(x) is equivalent to 2^bc.scale * ulp(x / 2^bc.scale). */
1334 sulp(U
*x
, BCinfo
*bc
)
1338 if (bc
->scale
&& 2*P
+ 1 > (int)((word0(x
) & Exp_mask
) >> Exp_shift
)) {
1339 /* rv/2^bc->scale is subnormal */
1340 word0(&u
) = (P
+2)*Exp_msk1
;
1345 assert(word0(x
) || word1(x
)); /* x != 0.0 */
1350 /* The bigcomp function handles some hard cases for strtod, for inputs
1351 with more than STRTOD_DIGLIM digits. It's called once an initial
1352 estimate for the double corresponding to the input string has
1353 already been obtained by the code in _Py_dg_strtod.
1355 The bigcomp function is only called after _Py_dg_strtod has found a
1356 double value rv such that either rv or rv + 1ulp represents the
1357 correctly rounded value corresponding to the original string. It
1358 determines which of these two values is the correct one by
1359 computing the decimal digits of rv + 0.5ulp and comparing them with
1360 the corresponding digits of s0.
1362 In the following, write dv for the absolute value of the number represented
1363 by the input string.
1367 s0 points to the first significant digit of the input string.
1369 rv is a (possibly scaled) estimate for the closest double value to the
1370 value represented by the original input to _Py_dg_strtod. If
1371 bc->scale is nonzero, then rv/2^(bc->scale) is the approximation to
1374 bc is a struct containing information gathered during the parsing and
1375 estimation steps of _Py_dg_strtod. Description of fields follows:
1377 bc->e0 gives the exponent of the input value, such that dv = (integer
1378 given by the bd->nd digits of s0) * 10**e0
1380 bc->nd gives the total number of significant digits of s0. It will
1383 bc->nd0 gives the number of significant digits of s0 before the
1384 decimal separator. If there's no decimal separator, bc->nd0 ==
1387 bc->scale is the value used to scale rv to avoid doing arithmetic with
1388 subnormal values. It's either 0 or 2*P (=106).
1392 On successful exit, rv/2^(bc->scale) is the closest double to dv.
1394 Returns 0 on success, -1 on failure (e.g., due to a failed malloc call). */
1397 bigcomp(U
*rv
, const char *s0
, BCinfo
*bc
)
1400 int b2
, d2
, dd
, i
, nd
, nd0
, odd
, p2
, p5
;
1405 b
= sd2b(rv
, bc
->scale
, &p2
);
1409 /* record whether the lsb of rv/2^(bc->scale) is odd: in the exact halfway
1410 case, this is used for round to even. */
1413 /* left shift b by 1 bit and or a 1 into the least significant bit;
1414 this gives us b * 2**p2 = rv/2^(bc->scale) + 0.5 ulp. */
1427 /* Arrange for convenient computation of quotients:
1428 * shift left if necessary so divisor has 4 leading 0 bits.
1431 d
= pow5mult(d
, p5
);
1438 b
= pow5mult(b
, -p5
);
1453 if ((b2
+= i
) > 0) {
1460 if ((d2
+= i
) > 0) {
1468 /* Compare s0 with b/d: set dd to -1, 0, or 1 according as s0 < b/d, s0 ==
1469 * b/d, or s0 > b/d. Here the digits of s0 are thought of as representing
1470 * a number in the range [0.1, 1). */
1477 b
= multadd(b
, 10, 0);
1482 dd
= s0
[i
< nd0
? i
: i
+1] - '0' - quorem(b
, d
);
1487 if (!b
->x
[0] && b
->wds
== 1) {
1493 /* b/d != 0, but digits of s0 exhausted */
1501 if (dd
> 0 || (dd
== 0 && odd
))
1502 dval(rv
) += sulp(rv
, bc
);
1507 _Py_dg_strtod(const char *s00
, char **se
)
1509 int bb2
, bb5
, bbe
, bd2
, bd5
, bs2
, c
, dsign
, e
, e1
, error
;
1510 int esign
, i
, j
, k
, lz
, nd
, nd0
, odd
, sign
;
1511 const char *s
, *s0
, *s1
;
1513 U aadj2
, adj
, rv
, rv0
;
1514 ULong y
, z
, abs_exp
;
1517 Bigint
*bb
, *bb1
, *bd
, *bd0
, *bs
, *delta
;
1518 size_t ndigits
, fraclen
;
1522 /* Start parsing. */
1525 /* Parse optional sign, if present. */
1535 /* Skip leading zeros: lz is true iff there were leading zeros. */
1541 /* Point s0 at the first nonzero digit (if any). fraclen will be the
1542 number of digits between the decimal point and the end of the
1543 digit string. ndigits will be the total number of digits ignoring
1546 while ('0' <= c
&& c
<= '9')
1551 /* Parse decimal point and following digits. */
1559 fraclen
+= (s
- s1
);
1563 while ('0' <= c
&& c
<= '9')
1569 /* Now lz is true if and only if there were leading zero digits, and
1570 ndigits gives the total number of digits ignoring leading zeros. A
1571 valid input must have at least one digit. */
1572 if (!ndigits
&& !lz
) {
1578 /* Range check ndigits and fraclen to make sure that they, and values
1579 computed with them, can safely fit in an int. */
1580 if (ndigits
> MAX_DIGITS
|| fraclen
> MAX_DIGITS
) {
1586 nd0
= (int)ndigits
- (int)fraclen
;
1588 /* Parse exponent. */
1590 if (c
== 'e' || c
== 'E') {
1594 /* Exponent sign. */
1604 /* Skip zeros. lz is true iff there are leading zeros. */
1610 /* Get absolute value of the exponent. */
1613 while ('0' <= c
&& c
<= '9') {
1614 abs_exp
= 10*abs_exp
+ (c
- '0');
1618 /* abs_exp will be correct modulo 2**32. But 10**9 < 2**32, so if
1619 there are at most 9 significant exponent digits then overflow is
1621 if (s
- s1
> 9 || abs_exp
> MAX_ABS_EXP
)
1622 e
= (int)MAX_ABS_EXP
;
1628 /* A valid exponent must have at least one digit. */
1633 /* Adjust exponent to take into account position of the point. */
1638 /* Finished parsing. Set se to indicate how far we parsed */
1642 /* If all digits were zero, exit with return value +-0.0. Otherwise,
1643 strip trailing zeros: scan back until we hit a nonzero digit. */
1646 for (i
= nd
; i
> 0; ) {
1648 if (s0
[i
< nd0
? i
: i
+1] != '0') {
1658 /* Summary of parsing results. After parsing, and dealing with zero
1659 * inputs, we have values s0, nd0, nd, e, sign, where:
1661 * - s0 points to the first significant digit of the input string
1663 * - nd is the total number of significant digits (here, and
1664 * below, 'significant digits' means the set of digits of the
1665 * significand of the input that remain after ignoring leading
1666 * and trailing zeros).
1668 * - nd0 indicates the position of the decimal point, if present; it
1669 * satisfies 1 <= nd0 <= nd. The nd significant digits are in
1670 * s0[0:nd0] and s0[nd0+1:nd+1] using the usual Python half-open slice
1671 * notation. (If nd0 < nd, then s0[nd0] contains a '.' character; if
1672 * nd0 == nd, then s0[nd0] could be any non-digit character.)
1674 * - e is the adjusted exponent: the absolute value of the number
1675 * represented by the original input string is n * 10**e, where
1676 * n is the integer represented by the concatenation of
1677 * s0[0:nd0] and s0[nd0+1:nd+1]
1679 * - sign gives the sign of the input: 1 for negative, 0 for positive
1681 * - the first and last significant digits are nonzero
1684 /* put first DBL_DIG+1 digits into integer y and z.
1686 * - y contains the value represented by the first min(9, nd)
1687 * significant digits
1689 * - if nd > 9, z contains the value represented by significant digits
1690 * with indices in [9, min(16, nd)). So y * 10**(min(16, nd) - 9) + z
1691 * gives the value represented by the first min(16, nd) sig. digits.
1696 for (i
= 0; i
< nd
; i
++) {
1698 y
= 10*y
+ s0
[i
< nd0
? i
: i
+1] - '0';
1699 else if (i
< DBL_DIG
+1)
1700 z
= 10*z
+ s0
[i
< nd0
? i
: i
+1] - '0';
1705 k
= nd
< DBL_DIG
+ 1 ? nd
: DBL_DIG
+ 1;
1708 dval(&rv
) = tens
[k
- 9] * dval(&rv
) + z
;
1717 if (e
<= Ten_pmax
) {
1718 dval(&rv
) *= tens
[e
];
1722 if (e
<= Ten_pmax
+ i
) {
1723 /* A fancier test would sometimes let us do
1724 * this for larger i values.
1727 dval(&rv
) *= tens
[i
];
1728 dval(&rv
) *= tens
[e
];
1732 else if (e
>= -Ten_pmax
) {
1733 dval(&rv
) /= tens
[-e
];
1741 /* Get starting approximation = rv * 10**e1 */
1745 dval(&rv
) *= tens
[i
];
1747 if (e1
> DBL_MAX_10_EXP
)
1750 for(j
= 0; e1
> 1; j
++, e1
>>= 1)
1752 dval(&rv
) *= bigtens
[j
];
1753 /* The last multiplication could overflow. */
1754 word0(&rv
) -= P
*Exp_msk1
;
1755 dval(&rv
) *= bigtens
[j
];
1756 if ((z
= word0(&rv
) & Exp_mask
)
1757 > Exp_msk1
*(DBL_MAX_EXP
+Bias
-P
))
1759 if (z
> Exp_msk1
*(DBL_MAX_EXP
+Bias
-1-P
)) {
1760 /* set to largest number */
1761 /* (Can't trust DBL_MAX) */
1766 word0(&rv
) += P
*Exp_msk1
;
1770 /* The input decimal value lies in [10**e1, 10**(e1+16)).
1772 If e1 <= -512, underflow immediately.
1773 If e1 <= -256, set bc.scale to 2*P.
1775 So for input value < 1e-256, bc.scale is always set;
1776 for input value >= 1e-240, bc.scale is never set.
1777 For input values in [1e-256, 1e-240), bc.scale may or may
1782 dval(&rv
) /= tens
[i
];
1784 if (e1
>= 1 << n_bigtens
)
1788 for(j
= 0; e1
> 0; j
++, e1
>>= 1)
1790 dval(&rv
) *= tinytens
[j
];
1791 if (bc
.scale
&& (j
= 2*P
+ 1 - ((word0(&rv
) & Exp_mask
)
1792 >> Exp_shift
)) > 0) {
1793 /* scaled rv is denormal; clear j low bits */
1797 word0(&rv
) = (P
+2)*Exp_msk1
;
1799 word0(&rv
) &= 0xffffffff << (j
-32);
1802 word1(&rv
) &= 0xffffffff << j
;
1809 /* Now the hard part -- adjusting rv to the correct value.*/
1811 /* Put digits into bd: true value = bd * 10^e */
1814 bc
.nd0
= nd0
; /* Only needed if nd > STRTOD_DIGLIM, but done here */
1815 /* to silence an erroneous warning about bc.nd0 */
1816 /* possibly not being initialized. */
1817 if (nd
> STRTOD_DIGLIM
) {
1818 /* ASSERT(STRTOD_DIGLIM >= 18); 18 == one more than the */
1819 /* minimum number of decimal digits to distinguish double values */
1820 /* in IEEE arithmetic. */
1822 /* Truncate input to 18 significant digits, then discard any trailing
1823 zeros on the result by updating nd, nd0, e and y suitably. (There's
1824 no need to update z; it's not reused beyond this point.) */
1825 for (i
= 18; i
> 0; ) {
1826 /* scan back until we hit a nonzero digit. significant digit 'i'
1827 is s0[i] if i < nd0, s0[i+1] if i >= nd0. */
1829 if (s0
[i
< nd0
? i
: i
+1] != '0') {
1838 if (nd
< 9) { /* must recompute y */
1840 for(i
= 0; i
< nd0
; ++i
)
1841 y
= 10*y
+ s0
[i
] - '0';
1843 y
= 10*y
+ s0
[i
+1] - '0';
1846 bd0
= s2b(s0
, nd0
, nd
, y
);
1850 /* Notation for the comments below. Write:
1852 - dv for the absolute value of the number represented by the original
1853 decimal input string.
1855 - if we've truncated dv, write tdv for the truncated value.
1856 Otherwise, set tdv == dv.
1858 - srv for the quantity rv/2^bc.scale; so srv is the current binary
1859 approximation to tdv (and dv). It should be exactly representable
1860 in an IEEE 754 double.
1865 /* This is the main correction loop for _Py_dg_strtod.
1867 We've got a decimal value tdv, and a floating-point approximation
1868 srv=rv/2^bc.scale to tdv. The aim is to determine whether srv is
1869 close enough (i.e., within 0.5 ulps) to tdv, and to compute a new
1870 approximation if not.
1872 To determine whether srv is close enough to tdv, compute integers
1873 bd, bb and bs proportional to tdv, srv and 0.5 ulp(srv)
1874 respectively, and then use integer arithmetic to determine whether
1875 |tdv - srv| is less than, equal to, or greater than 0.5 ulp(srv).
1878 bd
= Balloc(bd0
->k
);
1884 bb
= sd2b(&rv
, bc
.scale
, &bbe
); /* srv = bb * 2^bbe */
1890 /* Record whether lsb of bb is odd, in case we need this
1891 for the round-to-even step later. */
1894 /* tdv = bd * 10**e; srv = bb * 2**bbe */
1919 /* At this stage bd5 - bb5 == e == bd2 - bb2 + bbe, bb2 - bs2 == 1,
1922 tdv == bd * 10**e = bd * 2**(bbe - bb2 + bd2) * 5**(bd5 - bb5)
1923 srv == bb * 2**bbe = bb * 2**(bbe - bb2 + bb2)
1924 0.5 ulp(srv) == 2**(bbe-1) = bs * 2**(bbe - bb2 + bs2)
1928 M * tdv = bd * 2**bd2 * 5**bd5
1929 M * srv = bb * 2**bb2 * 5**bb5
1930 M * 0.5 ulp(srv) = bs * 2**bs2 * 5**bb5
1932 for some constant M. (Actually, M == 2**(bb2 - bbe) * 5**bb5, but
1933 this fact is not needed below.)
1936 /* Remove factor of 2**i, where i = min(bb2, bd2, bs2). */
1937 i
= bb2
< bd2
? bb2
: bd2
;
1946 /* Scale bb, bd, bs by the appropriate powers of 2 and 5. */
1948 bs
= pow5mult(bs
, bb5
);
1966 bb
= lshift(bb
, bb2
);
1975 bd
= pow5mult(bd
, bd5
);
1984 bd
= lshift(bd
, bd2
);
1993 bs
= lshift(bs
, bs2
);
2002 /* Now bd, bb and bs are scaled versions of tdv, srv and 0.5 ulp(srv),
2003 respectively. Compute the difference |tdv - srv|, and compare
2004 with 0.5 ulp(srv). */
2006 delta
= diff(bb
, bd
);
2007 if (delta
== NULL
) {
2014 dsign
= delta
->sign
;
2017 if (bc
.nd
> nd
&& i
<= 0) {
2019 break; /* Must use bigcomp(). */
2021 /* Here rv overestimates the truncated decimal value by at most
2022 0.5 ulp(rv). Hence rv either overestimates the true decimal
2023 value by <= 0.5 ulp(rv), or underestimates it by some small
2024 amount (< 0.1 ulp(rv)); either way, rv is within 0.5 ulps of
2025 the true decimal value, so it's possible to exit.
2027 Exception: if scaled rv is a normal exact power of 2, but not
2028 DBL_MIN, then rv - 0.5 ulp(rv) takes us all the way down to the
2029 next double, so the correctly rounded result is either rv - 0.5
2030 ulp(rv) or rv; in this case, use bigcomp to distinguish. */
2032 if (!word1(&rv
) && !(word0(&rv
) & Bndry_mask
)) {
2033 /* rv can't be 0, since it's an overestimate for some
2034 nonzero value. So rv is a normal power of 2. */
2035 j
= (int)(word0(&rv
) & Exp_mask
) >> Exp_shift
;
2036 /* rv / 2^bc.scale = 2^(j - 1023 - bc.scale); use bigcomp if
2037 rv / 2^bc.scale >= 2^-1021. */
2038 if (j
- bc
.scale
>= 2) {
2039 dval(&rv
) -= 0.5 * sulp(&rv
, &bc
);
2040 break; /* Use bigcomp. */
2046 i
= -1; /* Discarded digits make delta smaller. */
2051 /* Error is less than half an ulp -- check for
2052 * special case of mantissa a power of two.
2054 if (dsign
|| word1(&rv
) || word0(&rv
) & Bndry_mask
2055 || (word0(&rv
) & Exp_mask
) <= (2*P
+1)*Exp_msk1
2059 if (!delta
->x
[0] && delta
->wds
<= 1) {
2063 delta
= lshift(delta
,Log2P
);
2064 if (delta
== NULL
) {
2071 if (cmp(delta
, bs
) > 0)
2076 /* exactly half-way between */
2078 if ((word0(&rv
) & Bndry_mask1
) == Bndry_mask1
2081 (y
= word0(&rv
) & Exp_mask
) <= 2*P
*Exp_msk1
) ?
2082 (0xffffffff & (0xffffffff << (2*P
+1-(y
>>Exp_shift
)))) :
2084 /*boundary case -- increment exponent*/
2085 word0(&rv
) = (word0(&rv
) & Exp_mask
)
2093 else if (!(word0(&rv
) & Bndry_mask
) && !word1(&rv
)) {
2095 /* boundary case -- decrement exponent */
2097 L
= word0(&rv
) & Exp_mask
;
2098 if (L
<= (2*P
+1)*Exp_msk1
) {
2099 if (L
> (P
+2)*Exp_msk1
)
2100 /* round even ==> */
2103 /* rv = smallest denormal */
2109 L
= (word0(&rv
) & Exp_mask
) - Exp_msk1
;
2110 word0(&rv
) = L
| Bndry_mask1
;
2111 word1(&rv
) = 0xffffffff;
2117 dval(&rv
) += sulp(&rv
, &bc
);
2119 dval(&rv
) -= sulp(&rv
, &bc
);
2129 if ((aadj
= ratio(delta
, bs
)) <= 2.) {
2132 else if (word1(&rv
) || word0(&rv
) & Bndry_mask
) {
2133 if (word1(&rv
) == Tiny1
&& !word0(&rv
)) {
2142 /* special case -- power of FLT_RADIX to be */
2143 /* rounded down... */
2145 if (aadj
< 2./FLT_RADIX
)
2146 aadj
= 1./FLT_RADIX
;
2154 aadj1
= dsign
? aadj
: -aadj
;
2155 if (Flt_Rounds
== 0)
2158 y
= word0(&rv
) & Exp_mask
;
2160 /* Check for overflow */
2162 if (y
== Exp_msk1
*(DBL_MAX_EXP
+Bias
-1)) {
2163 dval(&rv0
) = dval(&rv
);
2164 word0(&rv
) -= P
*Exp_msk1
;
2165 adj
.d
= aadj1
* ulp(&rv
);
2167 if ((word0(&rv
) & Exp_mask
) >=
2168 Exp_msk1
*(DBL_MAX_EXP
+Bias
-P
)) {
2169 if (word0(&rv0
) == Big0
&& word1(&rv0
) == Big1
) {
2182 word0(&rv
) += P
*Exp_msk1
;
2185 if (bc
.scale
&& y
<= 2*P
*Exp_msk1
) {
2186 if (aadj
<= 0x7fffffff) {
2187 if ((z
= (ULong
)aadj
) <= 0)
2190 aadj1
= dsign
? aadj
: -aadj
;
2192 dval(&aadj2
) = aadj1
;
2193 word0(&aadj2
) += (2*P
+1)*Exp_msk1
- y
;
2194 aadj1
= dval(&aadj2
);
2196 adj
.d
= aadj1
* ulp(&rv
);
2199 z
= word0(&rv
) & Exp_mask
;
2203 /* Can we stop now? */
2206 /* The tolerances below are conservative. */
2207 if (dsign
|| word1(&rv
) || word0(&rv
) & Bndry_mask
) {
2208 if (aadj
< .4999999 || aadj
> .5000001)
2211 else if (aadj
< .4999999/FLT_RADIX
)
2227 error
= bigcomp(&rv
, s0
, &bc
);
2233 word0(&rv0
) = Exp_1
- 2*P
*Exp_msk1
;
2235 dval(&rv
) *= dval(&rv0
);
2239 return sign
? -dval(&rv
) : dval(&rv
);
2249 return sign
? -0.0 : 0.0;
2253 /* Can't trust HUGE_VAL */
2254 word0(&rv
) = Exp_mask
;
2256 return sign
? -dval(&rv
) : dval(&rv
);
2267 sizeof(Bigint
) - sizeof(ULong
) - sizeof(int) + j
<= (unsigned)i
;
2270 r
= (int*)Balloc(k
);
2274 return (char *)(r
+1);
2278 nrv_alloc(char *s
, char **rve
, int n
)
2286 while((*t
= *s
++)) t
++;
2292 /* freedtoa(s) must be used to free values s returned by dtoa
2293 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
2294 * but for consistency with earlier versions of dtoa, it is optional
2295 * when MULTIPLE_THREADS is not defined.
2299 _Py_dg_freedtoa(char *s
)
2301 Bigint
*b
= (Bigint
*)((int *)s
- 1);
2302 b
->maxwds
= 1 << (b
->k
= *(int*)b
);
2306 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
2308 * Inspired by "How to Print Floating-Point Numbers Accurately" by
2309 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
2312 * 1. Rather than iterating, we use a simple numeric overestimate
2313 * to determine k = floor(log10(d)). We scale relevant
2314 * quantities using O(log2(k)) rather than O(k) multiplications.
2315 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
2316 * try to generate digits strictly left to right. Instead, we
2317 * compute with fewer bits and propagate the carry if necessary
2318 * when rounding the final digit up. This is often faster.
2319 * 3. Under the assumption that input will be rounded nearest,
2320 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
2321 * That is, we allow equality in stopping tests when the
2322 * round-nearest rule will give the same floating-point value
2323 * as would satisfaction of the stopping test with strict
2325 * 4. We remove common factors of powers of 2 from relevant
2327 * 5. When converting floating-point integers less than 1e16,
2328 * we use floating-point arithmetic rather than resorting
2329 * to multiple-precision integers.
2330 * 6. When asked to produce fewer than 15 digits, we first try
2331 * to get by with floating-point arithmetic; we resort to
2332 * multiple-precision integer arithmetic only if we cannot
2333 * guarantee that the floating-point calculation has given
2334 * the correctly rounded result. For k requested digits and
2335 * "uniformly" distributed input, the probability is
2336 * something like 10^(k-15) that we must resort to the Long
2340 /* Additional notes (METD): (1) returns NULL on failure. (2) to avoid memory
2341 leakage, a successful call to _Py_dg_dtoa should always be matched by a
2342 call to _Py_dg_freedtoa. */
2345 _Py_dg_dtoa(double dd
, int mode
, int ndigits
,
2346 int *decpt
, int *sign
, char **rve
)
2348 /* Arguments ndigits, decpt, sign are similar to those
2349 of ecvt and fcvt; trailing zeros are suppressed from
2350 the returned string. If not null, *rve is set to point
2351 to the end of the return value. If d is +-Infinity or NaN,
2352 then *decpt is set to 9999.
2355 0 ==> shortest string that yields d when read in
2356 and rounded to nearest.
2357 1 ==> like 0, but with Steele & White stopping rule;
2358 e.g. with IEEE P754 arithmetic , mode 0 gives
2359 1e23 whereas mode 1 gives 9.999999999999999e22.
2360 2 ==> max(1,ndigits) significant digits. This gives a
2361 return value similar to that of ecvt, except
2362 that trailing zeros are suppressed.
2363 3 ==> through ndigits past the decimal point. This
2364 gives a return value similar to that from fcvt,
2365 except that trailing zeros are suppressed, and
2366 ndigits can be negative.
2367 4,5 ==> similar to 2 and 3, respectively, but (in
2368 round-nearest mode) with the tests of mode 0 to
2369 possibly return a shorter string that rounds to d.
2370 With IEEE arithmetic and compilation with
2371 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
2372 as modes 2 and 3 when FLT_ROUNDS != 1.
2373 6-9 ==> Debugging modes similar to mode - 4: don't try
2374 fast floating-point estimate (if applicable).
2376 Values of mode other than 0-9 are treated as mode 0.
2378 Sufficient space is allocated to the return value
2379 to hold the suppressed trailing zeros.
2382 int bbits
, b2
, b5
, be
, dig
, i
, ieps
, ilim
, ilim0
, ilim1
,
2383 j
, j1
, k
, k0
, k_check
, leftright
, m2
, m5
, s2
, s5
,
2384 spec_case
, try_quick
;
2388 Bigint
*b
, *b1
, *delta
, *mlo
, *mhi
, *S
;
2393 /* set pointers to NULL, to silence gcc compiler warnings and make
2394 cleanup easier on error */
2399 if (word0(&u
) & Sign_bit
) {
2400 /* set sign for everything, including 0's and NaNs */
2402 word0(&u
) &= ~Sign_bit
; /* clear sign bit */
2407 /* quick return for Infinities, NaNs and zeros */
2408 if ((word0(&u
) & Exp_mask
) == Exp_mask
)
2410 /* Infinity or NaN */
2412 if (!word1(&u
) && !(word0(&u
) & 0xfffff))
2413 return nrv_alloc("Infinity", rve
, 8);
2414 return nrv_alloc("NaN", rve
, 3);
2418 return nrv_alloc("0", rve
, 1);
2421 /* compute k = floor(log10(d)). The computation may leave k
2422 one too large, but should never leave k too small. */
2423 b
= d2b(&u
, &be
, &bbits
);
2426 if ((i
= (int)(word0(&u
) >> Exp_shift1
& (Exp_mask
>>Exp_shift1
)))) {
2427 dval(&d2
) = dval(&u
);
2428 word0(&d2
) &= Frac_mask1
;
2429 word0(&d2
) |= Exp_11
;
2431 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2432 * log10(x) = log(x) / log(10)
2433 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2434 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2436 * This suggests computing an approximation k to log10(d) by
2438 * k = (i - Bias)*0.301029995663981
2439 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2441 * We want k to be too large rather than too small.
2442 * The error in the first-order Taylor series approximation
2443 * is in our favor, so we just round up the constant enough
2444 * to compensate for any error in the multiplication of
2445 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2446 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2447 * adding 1e-13 to the constant term more than suffices.
2448 * Hence we adjust the constant term to 0.1760912590558.
2449 * (We could get a more accurate k by invoking log10,
2450 * but this is probably not worthwhile.)
2457 /* d is denormalized */
2459 i
= bbits
+ be
+ (Bias
+ (P
-1) - 1);
2460 x
= i
> 32 ? word0(&u
) << (64 - i
) | word1(&u
) >> (i
- 32)
2461 : word1(&u
) << (32 - i
);
2463 word0(&d2
) -= 31*Exp_msk1
; /* adjust exponent */
2464 i
-= (Bias
+ (P
-1) - 1) + 1;
2467 ds
= (dval(&d2
)-1.5)*0.289529654602168 + 0.1760912590558 +
2468 i
*0.301029995663981;
2470 if (ds
< 0. && ds
!= k
)
2471 k
--; /* want k = floor(ds) */
2473 if (k
>= 0 && k
<= Ten_pmax
) {
2474 if (dval(&u
) < tens
[k
])
2497 if (mode
< 0 || mode
> 9)
2507 ilim
= ilim1
= -1; /* Values for cases 0 and 1; done here to */
2508 /* silence erroneous "gcc -Wall" warning. */
2521 ilim
= ilim1
= i
= ndigits
;
2527 i
= ndigits
+ k
+ 1;
2539 if (ilim
>= 0 && ilim
<= Quick_max
&& try_quick
) {
2541 /* Try to get by with floating-point arithmetic. */
2544 dval(&d2
) = dval(&u
);
2547 ieps
= 2; /* conservative */
2552 /* prevent overflows */
2554 dval(&u
) /= bigtens
[n_bigtens
-1];
2557 for(; j
; j
>>= 1, i
++)
2564 else if ((j1
= -k
)) {
2565 dval(&u
) *= tens
[j1
& 0xf];
2566 for(j
= j1
>> 4; j
; j
>>= 1, i
++)
2569 dval(&u
) *= bigtens
[i
];
2572 if (k_check
&& dval(&u
) < 1. && ilim
> 0) {
2580 dval(&eps
) = ieps
*dval(&u
) + 7.;
2581 word0(&eps
) -= (P
-1)*Exp_msk1
;
2585 if (dval(&u
) > dval(&eps
))
2587 if (dval(&u
) < -dval(&eps
))
2592 /* Use Steele & White method of only
2593 * generating digits needed.
2595 dval(&eps
) = 0.5/tens
[ilim
-1] - dval(&eps
);
2599 *s
++ = '0' + (int)L
;
2600 if (dval(&u
) < dval(&eps
))
2602 if (1. - dval(&u
) < dval(&eps
))
2611 /* Generate ilim digits, then fix them up. */
2612 dval(&eps
) *= tens
[ilim
-1];
2613 for(i
= 1;; i
++, dval(&u
) *= 10.) {
2614 L
= (Long
)(dval(&u
));
2615 if (!(dval(&u
) -= L
))
2617 *s
++ = '0' + (int)L
;
2619 if (dval(&u
) > 0.5 + dval(&eps
))
2621 else if (dval(&u
) < 0.5 - dval(&eps
)) {
2632 dval(&u
) = dval(&d2
);
2637 /* Do we have a "small" integer? */
2639 if (be
>= 0 && k
<= Int_max
) {
2642 if (ndigits
< 0 && ilim
<= 0) {
2644 if (ilim
< 0 || dval(&u
) <= 5*ds
)
2648 for(i
= 1;; i
++, dval(&u
) *= 10.) {
2649 L
= (Long
)(dval(&u
) / ds
);
2651 *s
++ = '0' + (int)L
;
2656 dval(&u
) += dval(&u
);
2657 if (dval(&u
) > ds
|| (dval(&u
) == ds
&& L
& 1)) {
2677 denorm
? be
+ (Bias
+ (P
-1) - 1 + 1) :
2685 if (m2
> 0 && s2
> 0) {
2686 i
= m2
< s2
? m2
: s2
;
2694 mhi
= pow5mult(mhi
, m5
);
2703 if ((j
= b5
- m5
)) {
2710 b
= pow5mult(b
, b5
);
2719 S
= pow5mult(S
, s5
);
2724 /* Check for special case that d is a normalized power of 2. */
2727 if ((mode
< 2 || leftright
)
2729 if (!word1(&u
) && !(word0(&u
) & Bndry_mask
)
2730 && word0(&u
) & (Exp_mask
& ~Exp_msk1
)
2732 /* The special case */
2739 /* Arrange for convenient computation of quotients:
2740 * shift left if necessary so divisor has 4 leading 0 bits.
2742 * Perhaps we should just compute leading 28 bits of S once
2743 * and for all and pass them and a shift to quorem, so it
2744 * can do shifts and ors to compute the numerator for q.
2764 b
= multadd(b
, 10, 0); /* we botched the k estimate */
2768 mhi
= multadd(mhi
, 10, 0);
2775 if (ilim
<= 0 && (mode
== 3 || mode
== 5)) {
2777 /* no digits, fcvt style */
2783 S
= multadd(S
, 5, 0);
2796 mhi
= lshift(mhi
, m2
);
2801 /* Compute mlo -- check for special case
2802 * that d is a normalized power of 2.
2807 mhi
= Balloc(mhi
->k
);
2811 mhi
= lshift(mhi
, Log2P
);
2817 dig
= quorem(b
,S
) + '0';
2818 /* Do we yet have the shortest decimal string
2819 * that will round to d?
2822 delta
= diff(S
, mhi
);
2825 j1
= delta
->sign
? 1 : cmp(b
, delta
);
2827 if (j1
== 0 && mode
!= 1 && !(word1(&u
) & 1)
2836 if (j
< 0 || (j
== 0 && mode
!= 1
2839 if (!b
->x
[0] && b
->wds
<= 1) {
2847 if ((j1
> 0 || (j1
== 0 && dig
& 1))
2856 if (dig
== '9') { /* possible if i == 1 */
2867 b
= multadd(b
, 10, 0);
2871 mlo
= mhi
= multadd(mhi
, 10, 0);
2876 mlo
= multadd(mlo
, 10, 0);
2879 mhi
= multadd(mhi
, 10, 0);
2887 *s
++ = dig
= quorem(b
,S
) + '0';
2888 if (!b
->x
[0] && b
->wds
<= 1) {
2893 b
= multadd(b
, 10, 0);
2898 /* Round off last digit */
2904 if (j
> 0 || (j
== 0 && dig
& 1)) {
2921 if (mlo
&& mlo
!= mhi
)
2935 if (mlo
&& mlo
!= mhi
)
2942 _Py_dg_freedtoa(s0
);
2949 #endif /* PY_NO_SHORT_FLOAT_REPR */