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 Etiny (-1074) /* smallest denormal is 2**Etiny */
239 #define Exp_1 0x3ff00000
240 #define Exp_11 0x3ff00000
242 #define Frac_mask 0xfffff
243 #define Frac_mask1 0xfffff
246 #define Bndry_mask 0xfffff
247 #define Bndry_mask1 0xfffff
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 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 #ifndef Py_USING_MEMORY_DEBUGGER
313 /* Memory management: memory is allocated from, and returned to, Kmax+1 pools
314 of memory, where pool k (0 <= k <= Kmax) is for Bigints b with b->maxwds ==
315 1 << k. These pools are maintained as linked lists, with freelist[k]
316 pointing to the head of the list for pool k.
318 On allocation, if there's no free slot in the appropriate pool, MALLOC is
319 called to get more memory. This memory is not returned to the system until
320 Python quits. There's also a private memory pool that's allocated from
321 in preference to using MALLOC.
323 For Bigints with more than (1 << Kmax) digits (which implies at least 1233
324 decimal digits), memory is directly allocated using MALLOC, and freed using
327 XXX: it would be easy to bypass this memory-management system and
328 translate each call to Balloc into a call to PyMem_Malloc, and each
329 Bfree to PyMem_Free. Investigate whether this has any significant
330 performance on impact. */
332 static Bigint
*freelist
[Kmax
+1];
334 /* Allocate space for a Bigint with up to 1<<k digits */
343 if (k
<= Kmax
&& (rv
= freelist
[k
]))
344 freelist
[k
] = rv
->next
;
347 len
= (sizeof(Bigint
) + (x
-1)*sizeof(ULong
) + sizeof(double) - 1)
349 if (k
<= Kmax
&& pmem_next
- private_mem
+ len
<= PRIVATE_mem
) {
350 rv
= (Bigint
*)pmem_next
;
354 rv
= (Bigint
*)MALLOC(len
*sizeof(double));
361 rv
->sign
= rv
->wds
= 0;
365 /* Free a Bigint allocated with Balloc */
374 v
->next
= freelist
[v
->k
];
382 /* Alternative versions of Balloc and Bfree that use PyMem_Malloc and
383 PyMem_Free directly in place of the custom memory allocation scheme above.
384 These are provided for the benefit of memory debugging tools like
387 /* Allocate space for a Bigint with up to 1<<k digits */
397 len
= (sizeof(Bigint
) + (x
-1)*sizeof(ULong
) + sizeof(double) - 1)
400 rv
= (Bigint
*)MALLOC(len
*sizeof(double));
406 rv
->sign
= rv
->wds
= 0;
410 /* Free a Bigint allocated with Balloc */
420 #endif /* Py_USING_MEMORY_DEBUGGER */
422 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
423 y->wds*sizeof(Long) + 2*sizeof(int))
425 /* Multiply a Bigint b by m and add a. Either modifies b in place and returns
426 a pointer to the modified b, or Bfrees b and returns a pointer to a copy.
427 On failure, return NULL. In this case, b will have been already freed. */
430 multadd(Bigint
*b
, int m
, int a
) /* multiply by m and add a */
448 y
= *x
* (ULLong
)m
+ carry
;
450 *x
++ = (ULong
)(y
& FFFFFFFF
);
453 y
= (xi
& 0xffff) * m
+ carry
;
454 z
= (xi
>> 16) * m
+ (y
>> 16);
456 *x
++ = (z
<< 16) + (y
& 0xffff);
461 if (wds
>= b
->maxwds
) {
471 b
->x
[wds
++] = (ULong
)carry
;
477 /* convert a string s containing nd decimal digits (possibly containing a
478 decimal separator at position nd0, which is ignored) to a Bigint. This
479 function carries on where the parsing code in _Py_dg_strtod leaves off: on
480 entry, y9 contains the result of converting the first 9 digits. Returns
484 s2b(const char *s
, int nd0
, int nd
, ULong y9
)
491 for(k
= 0, y
= 1; x
> y
; y
<<= 1, k
++) ;
502 for (i
= 9; i
< nd0
; i
++) {
503 b
= multadd(b
, 10, *s
++ - '0');
509 b
= multadd(b
, 10, *s
++ - '0');
516 /* count leading 0 bits in the 32-bit integer x. */
523 if (!(x
& 0xffff0000)) {
527 if (!(x
& 0xff000000)) {
531 if (!(x
& 0xf0000000)) {
535 if (!(x
& 0xc0000000)) {
539 if (!(x
& 0x80000000)) {
541 if (!(x
& 0x40000000))
547 /* count trailing 0 bits in the 32-bit integer y, and shift y right by that
593 /* convert a small nonnegative integer to a Bigint */
608 /* multiply two Bigints. Returns a new Bigint, or NULL on failure. Ignores
609 the signs of a and b. */
612 mult(Bigint
*a
, Bigint
*b
)
616 ULong
*x
, *xa
, *xae
, *xb
, *xbe
, *xc
, *xc0
;
625 if ((!a
->x
[0] && a
->wds
== 1) || (!b
->x
[0] && b
->wds
== 1)) {
634 if (a
->wds
< b
->wds
) {
648 for(x
= c
->x
, xa
= x
+ wc
; x
< xa
; x
++)
656 for(; xb
< xbe
; xc0
++) {
662 z
= *x
++ * (ULLong
)y
+ *xc
+ carry
;
664 *xc
++ = (ULong
)(z
& FFFFFFFF
);
671 for(; xb
< xbe
; xb
++, xc0
++) {
672 if (y
= *xb
& 0xffff) {
677 z
= (*x
& 0xffff) * y
+ (*xc
& 0xffff) + carry
;
679 z2
= (*x
++ >> 16) * y
+ (*xc
>> 16) + carry
;
692 z
= (*x
& 0xffff) * y
+ (*xc
>> 16) + carry
;
695 z2
= (*x
++ >> 16) * y
+ (*xc
& 0xffff) + carry
;
703 for(xc0
= c
->x
, xc
= xc0
+ wc
; wc
> 0 && !*--xc
; --wc
) ;
708 #ifndef Py_USING_MEMORY_DEBUGGER
710 /* p5s is a linked list of powers of 5 of the form 5**(2**i), i >= 2 */
714 /* multiply the Bigint b by 5**k. Returns a pointer to the result, or NULL on
715 failure; if the returned pointer is distinct from b then the original
716 Bigint b will have been Bfree'd. Ignores the sign of b. */
719 pow5mult(Bigint
*b
, int k
)
721 Bigint
*b1
, *p5
, *p51
;
723 static int p05
[3] = { 5, 25, 125 };
726 b
= multadd(b
, p05
[i
-1], 0);
771 /* Version of pow5mult that doesn't cache powers of 5. Provided for
772 the benefit of memory debugging tools like Valgrind. */
775 pow5mult(Bigint
*b
, int k
)
777 Bigint
*b1
, *p5
, *p51
;
779 static int p05
[3] = { 5, 25, 125 };
782 b
= multadd(b
, p05
[i
-1], 0);
819 #endif /* Py_USING_MEMORY_DEBUGGER */
821 /* shift a Bigint b left by k bits. Return a pointer to the shifted result,
822 or NULL on failure. If the returned pointer is distinct from b then the
823 original b will have been Bfree'd. Ignores the sign of b. */
826 lshift(Bigint
*b
, int k
)
830 ULong
*x
, *x1
, *xe
, z
;
832 if (!k
|| (!b
->x
[0] && b
->wds
== 1))
838 for(i
= b
->maxwds
; n1
> i
; i
<<= 1)
846 for(i
= 0; i
< n
; i
++)
869 /* Do a three-way compare of a and b, returning -1 if a < b, 0 if a == b and
870 1 if a > b. Ignores signs of a and b. */
873 cmp(Bigint
*a
, Bigint
*b
)
875 ULong
*xa
, *xa0
, *xb
, *xb0
;
881 if (i
> 1 && !a
->x
[i
-1])
882 Bug("cmp called with a->x[a->wds-1] == 0");
883 if (j
> 1 && !b
->x
[j
-1])
884 Bug("cmp called with b->x[b->wds-1] == 0");
894 return *xa
< *xb
? -1 : 1;
901 /* Take the difference of Bigints a and b, returning a new Bigint. Returns
902 NULL on failure. The signs of a and b are ignored, but the sign of the
903 result is set appropriately. */
906 diff(Bigint
*a
, Bigint
*b
)
910 ULong
*xa
, *xae
, *xb
, *xbe
, *xc
;
949 y
= (ULLong
)*xa
++ - *xb
++ - borrow
;
950 borrow
= y
>> 32 & (ULong
)1;
951 *xc
++ = (ULong
)(y
& FFFFFFFF
);
956 borrow
= y
>> 32 & (ULong
)1;
957 *xc
++ = (ULong
)(y
& FFFFFFFF
);
961 y
= (*xa
& 0xffff) - (*xb
& 0xffff) - borrow
;
962 borrow
= (y
& 0x10000) >> 16;
963 z
= (*xa
++ >> 16) - (*xb
++ >> 16) - borrow
;
964 borrow
= (z
& 0x10000) >> 16;
969 y
= (*xa
& 0xffff) - borrow
;
970 borrow
= (y
& 0x10000) >> 16;
971 z
= (*xa
++ >> 16) - borrow
;
972 borrow
= (z
& 0x10000) >> 16;
982 /* Given a positive normal double x, return the difference between x and the
983 next double up. Doesn't give correct results for subnormals. */
991 L
= (word0(x
) & Exp_mask
) - (P
-1)*Exp_msk1
;
997 /* Convert a Bigint to a double plus an exponent */
1000 b2d(Bigint
*a
, int *e
)
1002 ULong
*xa
, *xa0
, w
, y
, z
;
1010 if (!y
) Bug("zero y in b2d");
1015 word0(&d
) = Exp_1
| y
>> (Ebits
- k
);
1016 w
= xa
> xa0
? *--xa
: 0;
1017 word1(&d
) = y
<< ((32-Ebits
) + k
) | w
>> (Ebits
- k
);
1020 z
= xa
> xa0
? *--xa
: 0;
1022 word0(&d
) = Exp_1
| y
<< k
| z
>> (32 - k
);
1023 y
= xa
> xa0
? *--xa
: 0;
1024 word1(&d
) = z
<< k
| y
>> (32 - k
);
1027 word0(&d
) = Exp_1
| y
;
1034 /* Convert a scaled double to a Bigint plus an exponent. Similar to d2b,
1035 except that it accepts the scale parameter used in _Py_dg_strtod (which
1036 should be either 0 or 2*P), and the normalization for the return value is
1037 different (see below). On input, d should be finite and nonnegative, and d
1038 / 2**scale should be exactly representable as an IEEE 754 double.
1040 Returns a Bigint b and an integer e such that
1042 dval(d) / 2**scale = b * 2**e.
1044 Unlike d2b, b is not necessarily odd: b and e are normalized so
1045 that either 2**(P-1) <= b < 2**P and e >= Etiny, or b < 2**P
1046 and e == Etiny. This applies equally to an input of 0.0: in that
1047 case the return values are b = 0 and e = Etiny.
1049 The above normalization ensures that for all possible inputs d,
1050 2**e gives ulp(d/2**scale).
1052 Returns NULL on failure.
1056 sd2b(U
*d
, int scale
, int *e
)
1064 /* First construct b and e assuming that scale == 0. */
1067 b
->x
[1] = word0(d
) & Frac_mask
;
1068 *e
= Etiny
- 1 + (int)((word0(d
) & Exp_mask
) >> Exp_shift
);
1072 b
->x
[1] |= Exp_msk1
;
1074 /* Now adjust for scale, provided that b != 0. */
1075 if (scale
&& (b
->x
[0] || b
->x
[1])) {
1080 /* We can't shift more than P-1 bits without shifting out a 1. */
1081 assert(0 < scale
&& scale
<= P
- 1);
1083 /* The bits shifted out should all be zero. */
1084 assert(b
->x
[0] == 0);
1090 /* The bits shifted out should all be zero. */
1091 assert(b
->x
[0] << (32 - scale
) == 0);
1092 b
->x
[0] = (b
->x
[0] >> scale
) | (b
->x
[1] << (32 - scale
));
1097 /* Ensure b is normalized. */
1104 /* Convert a double to a Bigint plus an exponent. Return NULL on failure.
1106 Given a finite nonzero double d, return an odd Bigint b and exponent *e
1107 such that fabs(d) = b * 2**e. On return, *bbits gives the number of
1108 significant bits of b; that is, 2**(*bbits-1) <= b < 2**(*bbits).
1110 If d is zero, then b == 0, *e == -1010, *bbits = 0.
1114 d2b(U
*d
, int *e
, int *bits
)
1126 z
= word0(d
) & Frac_mask
;
1127 word0(d
) &= 0x7fffffff; /* clear sign bit, which we ignore */
1128 if ((de
= (int)(word0(d
) >> Exp_shift
)))
1130 if ((y
= word1(d
))) {
1131 if ((k
= lo0bits(&y
))) {
1132 x
[0] = y
| z
<< (32 - k
);
1138 b
->wds
= (x
[1] = z
) ? 2 : 1;
1148 *e
= de
- Bias
- (P
-1) + k
;
1152 *e
= de
- Bias
- (P
-1) + 1 + k
;
1153 *bits
= 32*i
- hi0bits(x
[i
-1]);
1158 /* Compute the ratio of two Bigints, as a double. The result may have an
1159 error of up to 2.5 ulps. */
1162 ratio(Bigint
*a
, Bigint
*b
)
1167 dval(&da
) = b2d(a
, &ka
);
1168 dval(&db
) = b2d(b
, &kb
);
1169 k
= ka
- kb
+ 32*(a
->wds
- b
->wds
);
1171 word0(&da
) += k
*Exp_msk1
;
1174 word0(&db
) += k
*Exp_msk1
;
1176 return dval(&da
) / dval(&db
);
1181 1e0
, 1e1
, 1e2
, 1e3
, 1e4
, 1e5
, 1e6
, 1e7
, 1e8
, 1e9
,
1182 1e10
, 1e11
, 1e12
, 1e13
, 1e14
, 1e15
, 1e16
, 1e17
, 1e18
, 1e19
,
1187 bigtens
[] = { 1e16
, 1e32
, 1e64
, 1e128
, 1e256
};
1188 static const double tinytens
[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1189 9007199254740992.*9007199254740992.e
-256
1190 /* = 2^106 * 1e-256 */
1192 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1193 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1194 #define Scale_Bit 0x10
1203 dshift(Bigint
*b
, int p2
)
1205 int rv
= hi0bits(b
->x
[b
->wds
-1]) - 4;
1211 /* special case of Bigint division. The quotient is always in the range 0 <=
1212 quotient < 10, and on entry the divisor S is normalized so that its top 4
1213 bits (28--31) are zero and bit 27 is set. */
1216 quorem(Bigint
*b
, Bigint
*S
)
1219 ULong
*bx
, *bxe
, q
, *sx
, *sxe
;
1221 ULLong borrow
, carry
, y
, ys
;
1223 ULong borrow
, carry
, y
, ys
;
1229 /*debug*/ if (b
->wds
> n
)
1230 /*debug*/ Bug("oversize b in quorem");
1238 q
= *bxe
/ (*sxe
+ 1); /* ensure q <= true quotient */
1240 /*debug*/ if (q
> 9)
1241 /*debug*/ Bug("oversized quotient in quorem");
1248 ys
= *sx
++ * (ULLong
)q
+ carry
;
1250 y
= *bx
- (ys
& FFFFFFFF
) - borrow
;
1251 borrow
= y
>> 32 & (ULong
)1;
1252 *bx
++ = (ULong
)(y
& FFFFFFFF
);
1255 ys
= (si
& 0xffff) * q
+ carry
;
1256 zs
= (si
>> 16) * q
+ (ys
>> 16);
1258 y
= (*bx
& 0xffff) - (ys
& 0xffff) - borrow
;
1259 borrow
= (y
& 0x10000) >> 16;
1260 z
= (*bx
>> 16) - (zs
& 0xffff) - borrow
;
1261 borrow
= (z
& 0x10000) >> 16;
1268 while(--bxe
> bx
&& !*bxe
)
1273 if (cmp(b
, S
) >= 0) {
1283 y
= *bx
- (ys
& FFFFFFFF
) - borrow
;
1284 borrow
= y
>> 32 & (ULong
)1;
1285 *bx
++ = (ULong
)(y
& FFFFFFFF
);
1288 ys
= (si
& 0xffff) + carry
;
1289 zs
= (si
>> 16) + (ys
>> 16);
1291 y
= (*bx
& 0xffff) - (ys
& 0xffff) - borrow
;
1292 borrow
= (y
& 0x10000) >> 16;
1293 z
= (*bx
>> 16) - (zs
& 0xffff) - borrow
;
1294 borrow
= (z
& 0x10000) >> 16;
1302 while(--bxe
> bx
&& !*bxe
)
1310 /* sulp(x) is a version of ulp(x) that takes bc.scale into account.
1312 Assuming that x is finite and nonnegative (positive zero is fine
1313 here) and x / 2^bc.scale is exactly representable as a double,
1314 sulp(x) is equivalent to 2^bc.scale * ulp(x / 2^bc.scale). */
1317 sulp(U
*x
, BCinfo
*bc
)
1321 if (bc
->scale
&& 2*P
+ 1 > (int)((word0(x
) & Exp_mask
) >> Exp_shift
)) {
1322 /* rv/2^bc->scale is subnormal */
1323 word0(&u
) = (P
+2)*Exp_msk1
;
1328 assert(word0(x
) || word1(x
)); /* x != 0.0 */
1333 /* The bigcomp function handles some hard cases for strtod, for inputs
1334 with more than STRTOD_DIGLIM digits. It's called once an initial
1335 estimate for the double corresponding to the input string has
1336 already been obtained by the code in _Py_dg_strtod.
1338 The bigcomp function is only called after _Py_dg_strtod has found a
1339 double value rv such that either rv or rv + 1ulp represents the
1340 correctly rounded value corresponding to the original string. It
1341 determines which of these two values is the correct one by
1342 computing the decimal digits of rv + 0.5ulp and comparing them with
1343 the corresponding digits of s0.
1345 In the following, write dv for the absolute value of the number represented
1346 by the input string.
1350 s0 points to the first significant digit of the input string.
1352 rv is a (possibly scaled) estimate for the closest double value to the
1353 value represented by the original input to _Py_dg_strtod. If
1354 bc->scale is nonzero, then rv/2^(bc->scale) is the approximation to
1357 bc is a struct containing information gathered during the parsing and
1358 estimation steps of _Py_dg_strtod. Description of fields follows:
1360 bc->e0 gives the exponent of the input value, such that dv = (integer
1361 given by the bd->nd digits of s0) * 10**e0
1363 bc->nd gives the total number of significant digits of s0. It will
1366 bc->nd0 gives the number of significant digits of s0 before the
1367 decimal separator. If there's no decimal separator, bc->nd0 ==
1370 bc->scale is the value used to scale rv to avoid doing arithmetic with
1371 subnormal values. It's either 0 or 2*P (=106).
1375 On successful exit, rv/2^(bc->scale) is the closest double to dv.
1377 Returns 0 on success, -1 on failure (e.g., due to a failed malloc call). */
1380 bigcomp(U
*rv
, const char *s0
, BCinfo
*bc
)
1383 int b2
, d2
, dd
, i
, nd
, nd0
, odd
, p2
, p5
;
1385 dd
= 0; /* silence compiler warning about possibly unused variable */
1389 b
= sd2b(rv
, bc
->scale
, &p2
);
1393 /* record whether the lsb of rv/2^(bc->scale) is odd: in the exact halfway
1394 case, this is used for round to even. */
1397 /* left shift b by 1 bit and or a 1 into the least significant bit;
1398 this gives us b * 2**p2 = rv/2^(bc->scale) + 0.5 ulp. */
1411 /* Arrange for convenient computation of quotients:
1412 * shift left if necessary so divisor has 4 leading 0 bits.
1415 d
= pow5mult(d
, p5
);
1422 b
= pow5mult(b
, -p5
);
1437 if ((b2
+= i
) > 0) {
1444 if ((d2
+= i
) > 0) {
1452 /* Compare s0 with b/d: set dd to -1, 0, or 1 according as s0 < b/d, s0 ==
1453 * b/d, or s0 > b/d. Here the digits of s0 are thought of as representing
1454 * a number in the range [0.1, 1). */
1461 b
= multadd(b
, 10, 0);
1466 dd
= s0
[i
< nd0
? i
: i
+1] - '0' - quorem(b
, d
);
1471 if (!b
->x
[0] && b
->wds
== 1) {
1477 /* b/d != 0, but digits of s0 exhausted */
1485 if (dd
> 0 || (dd
== 0 && odd
))
1486 dval(rv
) += sulp(rv
, bc
);
1491 _Py_dg_strtod(const char *s00
, char **se
)
1493 int bb2
, bb5
, bbe
, bd2
, bd5
, bs2
, c
, dsign
, e
, e1
, error
;
1494 int esign
, i
, j
, k
, lz
, nd
, nd0
, odd
, sign
;
1495 const char *s
, *s0
, *s1
;
1497 U aadj2
, adj
, rv
, rv0
;
1498 ULong y
, z
, abs_exp
;
1501 Bigint
*bb
, *bb1
, *bd
, *bd0
, *bs
, *delta
;
1505 /* Start parsing. */
1508 /* Parse optional sign, if present. */
1518 /* Skip leading zeros: lz is true iff there were leading zeros. */
1524 /* Point s0 at the first nonzero digit (if any). nd0 will be the position
1525 of the point relative to s0. nd will be the total number of digits
1526 ignoring leading zeros. */
1528 while ('0' <= c
&& c
<= '9')
1532 /* Parse decimal point and following digits. */
1544 while ('0' <= c
&& c
<= '9')
1549 /* Now lz is true if and only if there were leading zero digits, and nd
1550 gives the total number of digits ignoring leading zeros. A valid input
1551 must have at least one digit. */
1558 /* Parse exponent. */
1560 if (c
== 'e' || c
== 'E') {
1564 /* Exponent sign. */
1574 /* Skip zeros. lz is true iff there are leading zeros. */
1580 /* Get absolute value of the exponent. */
1583 while ('0' <= c
&& c
<= '9') {
1584 abs_exp
= 10*abs_exp
+ (c
- '0');
1588 /* abs_exp will be correct modulo 2**32. But 10**9 < 2**32, so if
1589 there are at most 9 significant exponent digits then overflow is
1591 if (s
- s1
> 9 || abs_exp
> MAX_ABS_EXP
)
1592 e
= (int)MAX_ABS_EXP
;
1598 /* A valid exponent must have at least one digit. */
1603 /* Adjust exponent to take into account position of the point. */
1608 /* Finished parsing. Set se to indicate how far we parsed */
1612 /* If all digits were zero, exit with return value +-0.0. Otherwise,
1613 strip trailing zeros: scan back until we hit a nonzero digit. */
1616 for (i
= nd
; i
> 0; ) {
1618 if (s0
[i
< nd0
? i
: i
+1] != '0') {
1628 /* Summary of parsing results. After parsing, and dealing with zero
1629 * inputs, we have values s0, nd0, nd, e, sign, where:
1631 * - s0 points to the first significant digit of the input string
1633 * - nd is the total number of significant digits (here, and
1634 * below, 'significant digits' means the set of digits of the
1635 * significand of the input that remain after ignoring leading
1636 * and trailing zeros).
1638 * - nd0 indicates the position of the decimal point, if present; it
1639 * satisfies 1 <= nd0 <= nd. The nd significant digits are in
1640 * s0[0:nd0] and s0[nd0+1:nd+1] using the usual Python half-open slice
1641 * notation. (If nd0 < nd, then s0[nd0] contains a '.' character; if
1642 * nd0 == nd, then s0[nd0] could be any non-digit character.)
1644 * - e is the adjusted exponent: the absolute value of the number
1645 * represented by the original input string is n * 10**e, where
1646 * n is the integer represented by the concatenation of
1647 * s0[0:nd0] and s0[nd0+1:nd+1]
1649 * - sign gives the sign of the input: 1 for negative, 0 for positive
1651 * - the first and last significant digits are nonzero
1654 /* put first DBL_DIG+1 digits into integer y and z.
1656 * - y contains the value represented by the first min(9, nd)
1657 * significant digits
1659 * - if nd > 9, z contains the value represented by significant digits
1660 * with indices in [9, min(16, nd)). So y * 10**(min(16, nd) - 9) + z
1661 * gives the value represented by the first min(16, nd) sig. digits.
1666 for (i
= 0; i
< nd
; i
++) {
1668 y
= 10*y
+ s0
[i
< nd0
? i
: i
+1] - '0';
1669 else if (i
< DBL_DIG
+1)
1670 z
= 10*z
+ s0
[i
< nd0
? i
: i
+1] - '0';
1675 k
= nd
< DBL_DIG
+ 1 ? nd
: DBL_DIG
+ 1;
1678 dval(&rv
) = tens
[k
- 9] * dval(&rv
) + z
;
1687 if (e
<= Ten_pmax
) {
1688 dval(&rv
) *= tens
[e
];
1692 if (e
<= Ten_pmax
+ i
) {
1693 /* A fancier test would sometimes let us do
1694 * this for larger i values.
1697 dval(&rv
) *= tens
[i
];
1698 dval(&rv
) *= tens
[e
];
1702 else if (e
>= -Ten_pmax
) {
1703 dval(&rv
) /= tens
[-e
];
1711 /* Get starting approximation = rv * 10**e1 */
1715 dval(&rv
) *= tens
[i
];
1717 if (e1
> DBL_MAX_10_EXP
)
1720 for(j
= 0; e1
> 1; j
++, e1
>>= 1)
1722 dval(&rv
) *= bigtens
[j
];
1723 /* The last multiplication could overflow. */
1724 word0(&rv
) -= P
*Exp_msk1
;
1725 dval(&rv
) *= bigtens
[j
];
1726 if ((z
= word0(&rv
) & Exp_mask
)
1727 > Exp_msk1
*(DBL_MAX_EXP
+Bias
-P
))
1729 if (z
> Exp_msk1
*(DBL_MAX_EXP
+Bias
-1-P
)) {
1730 /* set to largest number */
1731 /* (Can't trust DBL_MAX) */
1736 word0(&rv
) += P
*Exp_msk1
;
1740 /* The input decimal value lies in [10**e1, 10**(e1+16)).
1742 If e1 <= -512, underflow immediately.
1743 If e1 <= -256, set bc.scale to 2*P.
1745 So for input value < 1e-256, bc.scale is always set;
1746 for input value >= 1e-240, bc.scale is never set.
1747 For input values in [1e-256, 1e-240), bc.scale may or may
1752 dval(&rv
) /= tens
[i
];
1754 if (e1
>= 1 << n_bigtens
)
1758 for(j
= 0; e1
> 0; j
++, e1
>>= 1)
1760 dval(&rv
) *= tinytens
[j
];
1761 if (bc
.scale
&& (j
= 2*P
+ 1 - ((word0(&rv
) & Exp_mask
)
1762 >> Exp_shift
)) > 0) {
1763 /* scaled rv is denormal; clear j low bits */
1767 word0(&rv
) = (P
+2)*Exp_msk1
;
1769 word0(&rv
) &= 0xffffffff << (j
-32);
1772 word1(&rv
) &= 0xffffffff << j
;
1779 /* Now the hard part -- adjusting rv to the correct value.*/
1781 /* Put digits into bd: true value = bd * 10^e */
1784 bc
.nd0
= nd0
; /* Only needed if nd > STRTOD_DIGLIM, but done here */
1785 /* to silence an erroneous warning about bc.nd0 */
1786 /* possibly not being initialized. */
1787 if (nd
> STRTOD_DIGLIM
) {
1788 /* ASSERT(STRTOD_DIGLIM >= 18); 18 == one more than the */
1789 /* minimum number of decimal digits to distinguish double values */
1790 /* in IEEE arithmetic. */
1792 /* Truncate input to 18 significant digits, then discard any trailing
1793 zeros on the result by updating nd, nd0, e and y suitably. (There's
1794 no need to update z; it's not reused beyond this point.) */
1795 for (i
= 18; i
> 0; ) {
1796 /* scan back until we hit a nonzero digit. significant digit 'i'
1797 is s0[i] if i < nd0, s0[i+1] if i >= nd0. */
1799 if (s0
[i
< nd0
? i
: i
+1] != '0') {
1808 if (nd
< 9) { /* must recompute y */
1810 for(i
= 0; i
< nd0
; ++i
)
1811 y
= 10*y
+ s0
[i
] - '0';
1813 y
= 10*y
+ s0
[i
+1] - '0';
1816 bd0
= s2b(s0
, nd0
, nd
, y
);
1820 /* Notation for the comments below. Write:
1822 - dv for the absolute value of the number represented by the original
1823 decimal input string.
1825 - if we've truncated dv, write tdv for the truncated value.
1826 Otherwise, set tdv == dv.
1828 - srv for the quantity rv/2^bc.scale; so srv is the current binary
1829 approximation to tdv (and dv). It should be exactly representable
1830 in an IEEE 754 double.
1835 /* This is the main correction loop for _Py_dg_strtod.
1837 We've got a decimal value tdv, and a floating-point approximation
1838 srv=rv/2^bc.scale to tdv. The aim is to determine whether srv is
1839 close enough (i.e., within 0.5 ulps) to tdv, and to compute a new
1840 approximation if not.
1842 To determine whether srv is close enough to tdv, compute integers
1843 bd, bb and bs proportional to tdv, srv and 0.5 ulp(srv)
1844 respectively, and then use integer arithmetic to determine whether
1845 |tdv - srv| is less than, equal to, or greater than 0.5 ulp(srv).
1848 bd
= Balloc(bd0
->k
);
1854 bb
= sd2b(&rv
, bc
.scale
, &bbe
); /* srv = bb * 2^bbe */
1860 /* Record whether lsb of bb is odd, in case we need this
1861 for the round-to-even step later. */
1864 /* tdv = bd * 10**e; srv = bb * 2**bbe */
1889 /* At this stage bd5 - bb5 == e == bd2 - bb2 + bbe, bb2 - bs2 == 1,
1892 tdv == bd * 10**e = bd * 2**(bbe - bb2 + bd2) * 5**(bd5 - bb5)
1893 srv == bb * 2**bbe = bb * 2**(bbe - bb2 + bb2)
1894 0.5 ulp(srv) == 2**(bbe-1) = bs * 2**(bbe - bb2 + bs2)
1898 M * tdv = bd * 2**bd2 * 5**bd5
1899 M * srv = bb * 2**bb2 * 5**bb5
1900 M * 0.5 ulp(srv) = bs * 2**bs2 * 5**bb5
1902 for some constant M. (Actually, M == 2**(bb2 - bbe) * 5**bb5, but
1903 this fact is not needed below.)
1906 /* Remove factor of 2**i, where i = min(bb2, bd2, bs2). */
1907 i
= bb2
< bd2
? bb2
: bd2
;
1916 /* Scale bb, bd, bs by the appropriate powers of 2 and 5. */
1918 bs
= pow5mult(bs
, bb5
);
1936 bb
= lshift(bb
, bb2
);
1945 bd
= pow5mult(bd
, bd5
);
1954 bd
= lshift(bd
, bd2
);
1963 bs
= lshift(bs
, bs2
);
1972 /* Now bd, bb and bs are scaled versions of tdv, srv and 0.5 ulp(srv),
1973 respectively. Compute the difference |tdv - srv|, and compare
1974 with 0.5 ulp(srv). */
1976 delta
= diff(bb
, bd
);
1977 if (delta
== NULL
) {
1984 dsign
= delta
->sign
;
1987 if (bc
.nd
> nd
&& i
<= 0) {
1989 break; /* Must use bigcomp(). */
1991 /* Here rv overestimates the truncated decimal value by at most
1992 0.5 ulp(rv). Hence rv either overestimates the true decimal
1993 value by <= 0.5 ulp(rv), or underestimates it by some small
1994 amount (< 0.1 ulp(rv)); either way, rv is within 0.5 ulps of
1995 the true decimal value, so it's possible to exit.
1997 Exception: if scaled rv is a normal exact power of 2, but not
1998 DBL_MIN, then rv - 0.5 ulp(rv) takes us all the way down to the
1999 next double, so the correctly rounded result is either rv - 0.5
2000 ulp(rv) or rv; in this case, use bigcomp to distinguish. */
2002 if (!word1(&rv
) && !(word0(&rv
) & Bndry_mask
)) {
2003 /* rv can't be 0, since it's an overestimate for some
2004 nonzero value. So rv is a normal power of 2. */
2005 j
= (int)(word0(&rv
) & Exp_mask
) >> Exp_shift
;
2006 /* rv / 2^bc.scale = 2^(j - 1023 - bc.scale); use bigcomp if
2007 rv / 2^bc.scale >= 2^-1021. */
2008 if (j
- bc
.scale
>= 2) {
2009 dval(&rv
) -= 0.5 * sulp(&rv
, &bc
);
2010 break; /* Use bigcomp. */
2016 i
= -1; /* Discarded digits make delta smaller. */
2021 /* Error is less than half an ulp -- check for
2022 * special case of mantissa a power of two.
2024 if (dsign
|| word1(&rv
) || word0(&rv
) & Bndry_mask
2025 || (word0(&rv
) & Exp_mask
) <= (2*P
+1)*Exp_msk1
2029 if (!delta
->x
[0] && delta
->wds
<= 1) {
2033 delta
= lshift(delta
,Log2P
);
2034 if (delta
== NULL
) {
2041 if (cmp(delta
, bs
) > 0)
2046 /* exactly half-way between */
2048 if ((word0(&rv
) & Bndry_mask1
) == Bndry_mask1
2051 (y
= word0(&rv
) & Exp_mask
) <= 2*P
*Exp_msk1
) ?
2052 (0xffffffff & (0xffffffff << (2*P
+1-(y
>>Exp_shift
)))) :
2054 /*boundary case -- increment exponent*/
2055 word0(&rv
) = (word0(&rv
) & Exp_mask
)
2063 else if (!(word0(&rv
) & Bndry_mask
) && !word1(&rv
)) {
2065 /* boundary case -- decrement exponent */
2067 L
= word0(&rv
) & Exp_mask
;
2068 if (L
<= (2*P
+1)*Exp_msk1
) {
2069 if (L
> (P
+2)*Exp_msk1
)
2070 /* round even ==> */
2073 /* rv = smallest denormal */
2079 L
= (word0(&rv
) & Exp_mask
) - Exp_msk1
;
2080 word0(&rv
) = L
| Bndry_mask1
;
2081 word1(&rv
) = 0xffffffff;
2087 dval(&rv
) += sulp(&rv
, &bc
);
2089 dval(&rv
) -= sulp(&rv
, &bc
);
2099 if ((aadj
= ratio(delta
, bs
)) <= 2.) {
2102 else if (word1(&rv
) || word0(&rv
) & Bndry_mask
) {
2103 if (word1(&rv
) == Tiny1
&& !word0(&rv
)) {
2112 /* special case -- power of FLT_RADIX to be */
2113 /* rounded down... */
2115 if (aadj
< 2./FLT_RADIX
)
2116 aadj
= 1./FLT_RADIX
;
2124 aadj1
= dsign
? aadj
: -aadj
;
2125 if (Flt_Rounds
== 0)
2128 y
= word0(&rv
) & Exp_mask
;
2130 /* Check for overflow */
2132 if (y
== Exp_msk1
*(DBL_MAX_EXP
+Bias
-1)) {
2133 dval(&rv0
) = dval(&rv
);
2134 word0(&rv
) -= P
*Exp_msk1
;
2135 adj
.d
= aadj1
* ulp(&rv
);
2137 if ((word0(&rv
) & Exp_mask
) >=
2138 Exp_msk1
*(DBL_MAX_EXP
+Bias
-P
)) {
2139 if (word0(&rv0
) == Big0
&& word1(&rv0
) == Big1
) {
2152 word0(&rv
) += P
*Exp_msk1
;
2155 if (bc
.scale
&& y
<= 2*P
*Exp_msk1
) {
2156 if (aadj
<= 0x7fffffff) {
2157 if ((z
= (ULong
)aadj
) <= 0)
2160 aadj1
= dsign
? aadj
: -aadj
;
2162 dval(&aadj2
) = aadj1
;
2163 word0(&aadj2
) += (2*P
+1)*Exp_msk1
- y
;
2164 aadj1
= dval(&aadj2
);
2166 adj
.d
= aadj1
* ulp(&rv
);
2169 z
= word0(&rv
) & Exp_mask
;
2173 /* Can we stop now? */
2176 /* The tolerances below are conservative. */
2177 if (dsign
|| word1(&rv
) || word0(&rv
) & Bndry_mask
) {
2178 if (aadj
< .4999999 || aadj
> .5000001)
2181 else if (aadj
< .4999999/FLT_RADIX
)
2197 error
= bigcomp(&rv
, s0
, &bc
);
2203 word0(&rv0
) = Exp_1
- 2*P
*Exp_msk1
;
2205 dval(&rv
) *= dval(&rv0
);
2209 return sign
? -dval(&rv
) : dval(&rv
);
2219 return sign
? -0.0 : 0.0;
2223 /* Can't trust HUGE_VAL */
2224 word0(&rv
) = Exp_mask
;
2226 return sign
? -dval(&rv
) : dval(&rv
);
2237 sizeof(Bigint
) - sizeof(ULong
) - sizeof(int) + j
<= (unsigned)i
;
2240 r
= (int*)Balloc(k
);
2244 return (char *)(r
+1);
2248 nrv_alloc(char *s
, char **rve
, int n
)
2256 while((*t
= *s
++)) t
++;
2262 /* freedtoa(s) must be used to free values s returned by dtoa
2263 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
2264 * but for consistency with earlier versions of dtoa, it is optional
2265 * when MULTIPLE_THREADS is not defined.
2269 _Py_dg_freedtoa(char *s
)
2271 Bigint
*b
= (Bigint
*)((int *)s
- 1);
2272 b
->maxwds
= 1 << (b
->k
= *(int*)b
);
2276 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
2278 * Inspired by "How to Print Floating-Point Numbers Accurately" by
2279 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
2282 * 1. Rather than iterating, we use a simple numeric overestimate
2283 * to determine k = floor(log10(d)). We scale relevant
2284 * quantities using O(log2(k)) rather than O(k) multiplications.
2285 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
2286 * try to generate digits strictly left to right. Instead, we
2287 * compute with fewer bits and propagate the carry if necessary
2288 * when rounding the final digit up. This is often faster.
2289 * 3. Under the assumption that input will be rounded nearest,
2290 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
2291 * That is, we allow equality in stopping tests when the
2292 * round-nearest rule will give the same floating-point value
2293 * as would satisfaction of the stopping test with strict
2295 * 4. We remove common factors of powers of 2 from relevant
2297 * 5. When converting floating-point integers less than 1e16,
2298 * we use floating-point arithmetic rather than resorting
2299 * to multiple-precision integers.
2300 * 6. When asked to produce fewer than 15 digits, we first try
2301 * to get by with floating-point arithmetic; we resort to
2302 * multiple-precision integer arithmetic only if we cannot
2303 * guarantee that the floating-point calculation has given
2304 * the correctly rounded result. For k requested digits and
2305 * "uniformly" distributed input, the probability is
2306 * something like 10^(k-15) that we must resort to the Long
2310 /* Additional notes (METD): (1) returns NULL on failure. (2) to avoid memory
2311 leakage, a successful call to _Py_dg_dtoa should always be matched by a
2312 call to _Py_dg_freedtoa. */
2315 _Py_dg_dtoa(double dd
, int mode
, int ndigits
,
2316 int *decpt
, int *sign
, char **rve
)
2318 /* Arguments ndigits, decpt, sign are similar to those
2319 of ecvt and fcvt; trailing zeros are suppressed from
2320 the returned string. If not null, *rve is set to point
2321 to the end of the return value. If d is +-Infinity or NaN,
2322 then *decpt is set to 9999.
2325 0 ==> shortest string that yields d when read in
2326 and rounded to nearest.
2327 1 ==> like 0, but with Steele & White stopping rule;
2328 e.g. with IEEE P754 arithmetic , mode 0 gives
2329 1e23 whereas mode 1 gives 9.999999999999999e22.
2330 2 ==> max(1,ndigits) significant digits. This gives a
2331 return value similar to that of ecvt, except
2332 that trailing zeros are suppressed.
2333 3 ==> through ndigits past the decimal point. This
2334 gives a return value similar to that from fcvt,
2335 except that trailing zeros are suppressed, and
2336 ndigits can be negative.
2337 4,5 ==> similar to 2 and 3, respectively, but (in
2338 round-nearest mode) with the tests of mode 0 to
2339 possibly return a shorter string that rounds to d.
2340 With IEEE arithmetic and compilation with
2341 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
2342 as modes 2 and 3 when FLT_ROUNDS != 1.
2343 6-9 ==> Debugging modes similar to mode - 4: don't try
2344 fast floating-point estimate (if applicable).
2346 Values of mode other than 0-9 are treated as mode 0.
2348 Sufficient space is allocated to the return value
2349 to hold the suppressed trailing zeros.
2352 int bbits
, b2
, b5
, be
, dig
, i
, ieps
, ilim
, ilim0
, ilim1
,
2353 j
, j1
, k
, k0
, k_check
, leftright
, m2
, m5
, s2
, s5
,
2354 spec_case
, try_quick
;
2358 Bigint
*b
, *b1
, *delta
, *mlo
, *mhi
, *S
;
2363 /* set pointers to NULL, to silence gcc compiler warnings and make
2364 cleanup easier on error */
2365 mlo
= mhi
= b
= S
= 0;
2369 if (word0(&u
) & Sign_bit
) {
2370 /* set sign for everything, including 0's and NaNs */
2372 word0(&u
) &= ~Sign_bit
; /* clear sign bit */
2377 /* quick return for Infinities, NaNs and zeros */
2378 if ((word0(&u
) & Exp_mask
) == Exp_mask
)
2380 /* Infinity or NaN */
2382 if (!word1(&u
) && !(word0(&u
) & 0xfffff))
2383 return nrv_alloc("Infinity", rve
, 8);
2384 return nrv_alloc("NaN", rve
, 3);
2388 return nrv_alloc("0", rve
, 1);
2391 /* compute k = floor(log10(d)). The computation may leave k
2392 one too large, but should never leave k too small. */
2393 b
= d2b(&u
, &be
, &bbits
);
2396 if ((i
= (int)(word0(&u
) >> Exp_shift1
& (Exp_mask
>>Exp_shift1
)))) {
2397 dval(&d2
) = dval(&u
);
2398 word0(&d2
) &= Frac_mask1
;
2399 word0(&d2
) |= Exp_11
;
2401 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2402 * log10(x) = log(x) / log(10)
2403 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2404 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2406 * This suggests computing an approximation k to log10(d) by
2408 * k = (i - Bias)*0.301029995663981
2409 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2411 * We want k to be too large rather than too small.
2412 * The error in the first-order Taylor series approximation
2413 * is in our favor, so we just round up the constant enough
2414 * to compensate for any error in the multiplication of
2415 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2416 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2417 * adding 1e-13 to the constant term more than suffices.
2418 * Hence we adjust the constant term to 0.1760912590558.
2419 * (We could get a more accurate k by invoking log10,
2420 * but this is probably not worthwhile.)
2427 /* d is denormalized */
2429 i
= bbits
+ be
+ (Bias
+ (P
-1) - 1);
2430 x
= i
> 32 ? word0(&u
) << (64 - i
) | word1(&u
) >> (i
- 32)
2431 : word1(&u
) << (32 - i
);
2433 word0(&d2
) -= 31*Exp_msk1
; /* adjust exponent */
2434 i
-= (Bias
+ (P
-1) - 1) + 1;
2437 ds
= (dval(&d2
)-1.5)*0.289529654602168 + 0.1760912590558 +
2438 i
*0.301029995663981;
2440 if (ds
< 0. && ds
!= k
)
2441 k
--; /* want k = floor(ds) */
2443 if (k
>= 0 && k
<= Ten_pmax
) {
2444 if (dval(&u
) < tens
[k
])
2467 if (mode
< 0 || mode
> 9)
2477 ilim
= ilim1
= -1; /* Values for cases 0 and 1; done here to */
2478 /* silence erroneous "gcc -Wall" warning. */
2491 ilim
= ilim1
= i
= ndigits
;
2497 i
= ndigits
+ k
+ 1;
2509 if (ilim
>= 0 && ilim
<= Quick_max
&& try_quick
) {
2511 /* Try to get by with floating-point arithmetic. */
2514 dval(&d2
) = dval(&u
);
2517 ieps
= 2; /* conservative */
2522 /* prevent overflows */
2524 dval(&u
) /= bigtens
[n_bigtens
-1];
2527 for(; j
; j
>>= 1, i
++)
2534 else if ((j1
= -k
)) {
2535 dval(&u
) *= tens
[j1
& 0xf];
2536 for(j
= j1
>> 4; j
; j
>>= 1, i
++)
2539 dval(&u
) *= bigtens
[i
];
2542 if (k_check
&& dval(&u
) < 1. && ilim
> 0) {
2550 dval(&eps
) = ieps
*dval(&u
) + 7.;
2551 word0(&eps
) -= (P
-1)*Exp_msk1
;
2555 if (dval(&u
) > dval(&eps
))
2557 if (dval(&u
) < -dval(&eps
))
2562 /* Use Steele & White method of only
2563 * generating digits needed.
2565 dval(&eps
) = 0.5/tens
[ilim
-1] - dval(&eps
);
2569 *s
++ = '0' + (int)L
;
2570 if (dval(&u
) < dval(&eps
))
2572 if (1. - dval(&u
) < dval(&eps
))
2581 /* Generate ilim digits, then fix them up. */
2582 dval(&eps
) *= tens
[ilim
-1];
2583 for(i
= 1;; i
++, dval(&u
) *= 10.) {
2584 L
= (Long
)(dval(&u
));
2585 if (!(dval(&u
) -= L
))
2587 *s
++ = '0' + (int)L
;
2589 if (dval(&u
) > 0.5 + dval(&eps
))
2591 else if (dval(&u
) < 0.5 - dval(&eps
)) {
2602 dval(&u
) = dval(&d2
);
2607 /* Do we have a "small" integer? */
2609 if (be
>= 0 && k
<= Int_max
) {
2612 if (ndigits
< 0 && ilim
<= 0) {
2614 if (ilim
< 0 || dval(&u
) <= 5*ds
)
2618 for(i
= 1;; i
++, dval(&u
) *= 10.) {
2619 L
= (Long
)(dval(&u
) / ds
);
2621 *s
++ = '0' + (int)L
;
2626 dval(&u
) += dval(&u
);
2627 if (dval(&u
) > ds
|| (dval(&u
) == ds
&& L
& 1)) {
2647 denorm
? be
+ (Bias
+ (P
-1) - 1 + 1) :
2655 if (m2
> 0 && s2
> 0) {
2656 i
= m2
< s2
? m2
: s2
;
2664 mhi
= pow5mult(mhi
, m5
);
2673 if ((j
= b5
- m5
)) {
2680 b
= pow5mult(b
, b5
);
2689 S
= pow5mult(S
, s5
);
2694 /* Check for special case that d is a normalized power of 2. */
2697 if ((mode
< 2 || leftright
)
2699 if (!word1(&u
) && !(word0(&u
) & Bndry_mask
)
2700 && word0(&u
) & (Exp_mask
& ~Exp_msk1
)
2702 /* The special case */
2709 /* Arrange for convenient computation of quotients:
2710 * shift left if necessary so divisor has 4 leading 0 bits.
2712 * Perhaps we should just compute leading 28 bits of S once
2713 * and for all and pass them and a shift to quorem, so it
2714 * can do shifts and ors to compute the numerator for q.
2716 if ((i
= ((s5
? 32 - hi0bits(S
->x
[S
->wds
-1]) : 1) + s2
) & 0x1f))
2736 b
= multadd(b
, 10, 0); /* we botched the k estimate */
2740 mhi
= multadd(mhi
, 10, 0);
2747 if (ilim
<= 0 && (mode
== 3 || mode
== 5)) {
2749 /* no digits, fcvt style */
2755 S
= multadd(S
, 5, 0);
2768 mhi
= lshift(mhi
, m2
);
2773 /* Compute mlo -- check for special case
2774 * that d is a normalized power of 2.
2779 mhi
= Balloc(mhi
->k
);
2783 mhi
= lshift(mhi
, Log2P
);
2789 dig
= quorem(b
,S
) + '0';
2790 /* Do we yet have the shortest decimal string
2791 * that will round to d?
2794 delta
= diff(S
, mhi
);
2797 j1
= delta
->sign
? 1 : cmp(b
, delta
);
2799 if (j1
== 0 && mode
!= 1 && !(word1(&u
) & 1)
2808 if (j
< 0 || (j
== 0 && mode
!= 1
2811 if (!b
->x
[0] && b
->wds
<= 1) {
2819 if ((j1
> 0 || (j1
== 0 && dig
& 1))
2828 if (dig
== '9') { /* possible if i == 1 */
2839 b
= multadd(b
, 10, 0);
2843 mlo
= mhi
= multadd(mhi
, 10, 0);
2848 mlo
= multadd(mlo
, 10, 0);
2851 mhi
= multadd(mhi
, 10, 0);
2859 *s
++ = dig
= quorem(b
,S
) + '0';
2860 if (!b
->x
[0] && b
->wds
<= 1) {
2865 b
= multadd(b
, 10, 0);
2870 /* Round off last digit */
2876 if (j
> 0 || (j
== 0 && dig
& 1)) {
2893 if (mlo
&& mlo
!= mhi
)
2907 if (mlo
&& mlo
!= mhi
)
2914 _Py_dg_freedtoa(s0
);
2921 #endif /* PY_NO_SHORT_FLOAT_REPR */