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
/* Standard NaN used by _Py_dg_stdnan. */
270
#define NAN_WORD0 0x7ff80000
273
/* Bits of the representation of positive infinity. */
275
#define POSINF_WORD0 0x7ff00000
276
#define POSINF_WORD1 0
278
/* struct BCinfo is used to pass information from _Py_dg_strtod to bigcomp */
280
typedef struct BCinfo BCinfo;
283
int e0, nd, nd0, scale;
286
#define FFFFFFFF 0xffffffffUL
290
/* struct Bigint is used to represent arbitrary-precision integers. These
291
integers are stored in sign-magnitude format, with the magnitude stored as
292
an array of base 2**32 digits. Bigints are always normalized: if x is a
293
Bigint then x->wds >= 1, and either x->wds == 1 or x[wds-1] is nonzero.
295
The Bigint fields are as follows:
297
- next is a header used by Balloc and Bfree to keep track of lists
298
of freed Bigints; it's also used for the linked list of
299
powers of 5 of the form 5**2**i used by pow5mult.
300
- k indicates which pool this Bigint was allocated from
301
- maxwds is the maximum number of words space was allocated for
302
(usually maxwds == 2**k)
303
- sign is 1 for negative Bigints, 0 for positive. The sign is unused
304
(ignored on inputs, set to 0 on outputs) in almost all operations
305
involving Bigints: a notable exception is the diff function, which
306
ignores signs on inputs but sets the sign of the output correctly.
307
- wds is the actual number of significant words
308
- x contains the vector of words (digits) for this Bigint, from least
309
significant (x[0]) to most significant (x[wds-1]).
315
int k, maxwds, sign, wds;
319
typedef struct Bigint Bigint;
321
#ifndef Py_USING_MEMORY_DEBUGGER
323
/* Memory management: memory is allocated from, and returned to, Kmax+1 pools
324
of memory, where pool k (0 <= k <= Kmax) is for Bigints b with b->maxwds ==
325
1 << k. These pools are maintained as linked lists, with freelist[k]
326
pointing to the head of the list for pool k.
328
On allocation, if there's no free slot in the appropriate pool, MALLOC is
329
called to get more memory. This memory is not returned to the system until
330
Python quits. There's also a private memory pool that's allocated from
331
in preference to using MALLOC.
333
For Bigints with more than (1 << Kmax) digits (which implies at least 1233
334
decimal digits), memory is directly allocated using MALLOC, and freed using
337
XXX: it would be easy to bypass this memory-management system and
338
translate each call to Balloc into a call to PyMem_Malloc, and each
339
Bfree to PyMem_Free. Investigate whether this has any significant
340
performance on impact. */
342
static Bigint *freelist[Kmax+1];
344
/* Allocate space for a Bigint with up to 1<<k digits */
353
if (k <= Kmax && (rv = freelist[k]))
354
freelist[k] = rv->next;
357
len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
359
if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
360
rv = (Bigint*)pmem_next;
364
rv = (Bigint*)MALLOC(len*sizeof(double));
371
rv->sign = rv->wds = 0;
375
/* Free a Bigint allocated with Balloc */
384
v->next = freelist[v->k];
392
/* Alternative versions of Balloc and Bfree that use PyMem_Malloc and
393
PyMem_Free directly in place of the custom memory allocation scheme above.
394
These are provided for the benefit of memory debugging tools like
397
/* Allocate space for a Bigint with up to 1<<k digits */
407
len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
410
rv = (Bigint*)MALLOC(len*sizeof(double));
416
rv->sign = rv->wds = 0;
420
/* Free a Bigint allocated with Balloc */
430
#endif /* Py_USING_MEMORY_DEBUGGER */
432
#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
433
y->wds*sizeof(Long) + 2*sizeof(int))
435
/* Multiply a Bigint b by m and add a. Either modifies b in place and returns
436
a pointer to the modified b, or Bfrees b and returns a pointer to a copy.
437
On failure, return NULL. In this case, b will have been already freed. */
440
multadd(Bigint *b, int m, int a) /* multiply by m and add a */
458
y = *x * (ULLong)m + carry;
460
*x++ = (ULong)(y & FFFFFFFF);
463
y = (xi & 0xffff) * m + carry;
464
z = (xi >> 16) * m + (y >> 16);
466
*x++ = (z << 16) + (y & 0xffff);
471
if (wds >= b->maxwds) {
481
b->x[wds++] = (ULong)carry;
487
/* convert a string s containing nd decimal digits (possibly containing a
488
decimal separator at position nd0, which is ignored) to a Bigint. This
489
function carries on where the parsing code in _Py_dg_strtod leaves off: on
490
entry, y9 contains the result of converting the first 9 digits. Returns
494
s2b(const char *s, int nd0, int nd, ULong y9)
501
for(k = 0, y = 1; x > y; y <<= 1, k++) ;
512
for (i = 9; i < nd0; i++) {
513
b = multadd(b, 10, *s++ - '0');
519
b = multadd(b, 10, *s++ - '0');
526
/* count leading 0 bits in the 32-bit integer x. */
533
if (!(x & 0xffff0000)) {
537
if (!(x & 0xff000000)) {
541
if (!(x & 0xf0000000)) {
545
if (!(x & 0xc0000000)) {
549
if (!(x & 0x80000000)) {
551
if (!(x & 0x40000000))
557
/* count trailing 0 bits in the 32-bit integer y, and shift y right by that
603
/* convert a small nonnegative integer to a Bigint */
618
/* multiply two Bigints. Returns a new Bigint, or NULL on failure. Ignores
619
the signs of a and b. */
622
mult(Bigint *a, Bigint *b)
626
ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
635
if ((!a->x[0] && a->wds == 1) || (!b->x[0] && b->wds == 1)) {
644
if (a->wds < b->wds) {
658
for(x = c->x, xa = x + wc; x < xa; x++)
666
for(; xb < xbe; xc0++) {
672
z = *x++ * (ULLong)y + *xc + carry;
674
*xc++ = (ULong)(z & FFFFFFFF);
681
for(; xb < xbe; xb++, xc0++) {
682
if (y = *xb & 0xffff) {
687
z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
689
z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
702
z = (*x & 0xffff) * y + (*xc >> 16) + carry;
705
z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
713
for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
718
#ifndef Py_USING_MEMORY_DEBUGGER
720
/* p5s is a linked list of powers of 5 of the form 5**(2**i), i >= 2 */
724
/* multiply the Bigint b by 5**k. Returns a pointer to the result, or NULL on
725
failure; if the returned pointer is distinct from b then the original
726
Bigint b will have been Bfree'd. Ignores the sign of b. */
729
pow5mult(Bigint *b, int k)
731
Bigint *b1, *p5, *p51;
733
static int p05[3] = { 5, 25, 125 };
736
b = multadd(b, p05[i-1], 0);
781
/* Version of pow5mult that doesn't cache powers of 5. Provided for
782
the benefit of memory debugging tools like Valgrind. */
785
pow5mult(Bigint *b, int k)
787
Bigint *b1, *p5, *p51;
789
static int p05[3] = { 5, 25, 125 };
792
b = multadd(b, p05[i-1], 0);
829
#endif /* Py_USING_MEMORY_DEBUGGER */
831
/* shift a Bigint b left by k bits. Return a pointer to the shifted result,
832
or NULL on failure. If the returned pointer is distinct from b then the
833
original b will have been Bfree'd. Ignores the sign of b. */
836
lshift(Bigint *b, int k)
840
ULong *x, *x1, *xe, z;
842
if (!k || (!b->x[0] && b->wds == 1))
848
for(i = b->maxwds; n1 > i; i <<= 1)
856
for(i = 0; i < n; i++)
879
/* Do a three-way compare of a and b, returning -1 if a < b, 0 if a == b and
880
1 if a > b. Ignores signs of a and b. */
883
cmp(Bigint *a, Bigint *b)
885
ULong *xa, *xa0, *xb, *xb0;
891
if (i > 1 && !a->x[i-1])
892
Bug("cmp called with a->x[a->wds-1] == 0");
893
if (j > 1 && !b->x[j-1])
894
Bug("cmp called with b->x[b->wds-1] == 0");
904
return *xa < *xb ? -1 : 1;
911
/* Take the difference of Bigints a and b, returning a new Bigint. Returns
912
NULL on failure. The signs of a and b are ignored, but the sign of the
913
result is set appropriately. */
916
diff(Bigint *a, Bigint *b)
920
ULong *xa, *xae, *xb, *xbe, *xc;
959
y = (ULLong)*xa++ - *xb++ - borrow;
960
borrow = y >> 32 & (ULong)1;
961
*xc++ = (ULong)(y & FFFFFFFF);
966
borrow = y >> 32 & (ULong)1;
967
*xc++ = (ULong)(y & FFFFFFFF);
971
y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
972
borrow = (y & 0x10000) >> 16;
973
z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
974
borrow = (z & 0x10000) >> 16;
979
y = (*xa & 0xffff) - borrow;
980
borrow = (y & 0x10000) >> 16;
981
z = (*xa++ >> 16) - borrow;
982
borrow = (z & 0x10000) >> 16;
992
/* Given a positive normal double x, return the difference between x and the
993
next double up. Doesn't give correct results for subnormals. */
1001
L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
1007
/* Convert a Bigint to a double plus an exponent */
1010
b2d(Bigint *a, int *e)
1012
ULong *xa, *xa0, w, y, z;
1020
if (!y) Bug("zero y in b2d");
1025
word0(&d) = Exp_1 | y >> (Ebits - k);
1026
w = xa > xa0 ? *--xa : 0;
1027
word1(&d) = y << ((32-Ebits) + k) | w >> (Ebits - k);
1030
z = xa > xa0 ? *--xa : 0;
1032
word0(&d) = Exp_1 | y << k | z >> (32 - k);
1033
y = xa > xa0 ? *--xa : 0;
1034
word1(&d) = z << k | y >> (32 - k);
1037
word0(&d) = Exp_1 | y;
1044
/* Convert a scaled double to a Bigint plus an exponent. Similar to d2b,
1045
except that it accepts the scale parameter used in _Py_dg_strtod (which
1046
should be either 0 or 2*P), and the normalization for the return value is
1047
different (see below). On input, d should be finite and nonnegative, and d
1048
/ 2**scale should be exactly representable as an IEEE 754 double.
1050
Returns a Bigint b and an integer e such that
1052
dval(d) / 2**scale = b * 2**e.
1054
Unlike d2b, b is not necessarily odd: b and e are normalized so
1055
that either 2**(P-1) <= b < 2**P and e >= Etiny, or b < 2**P
1056
and e == Etiny. This applies equally to an input of 0.0: in that
1057
case the return values are b = 0 and e = Etiny.
1059
The above normalization ensures that for all possible inputs d,
1060
2**e gives ulp(d/2**scale).
1062
Returns NULL on failure.
1066
sd2b(U *d, int scale, int *e)
1074
/* First construct b and e assuming that scale == 0. */
1077
b->x[1] = word0(d) & Frac_mask;
1078
*e = Etiny - 1 + (int)((word0(d) & Exp_mask) >> Exp_shift);
1082
b->x[1] |= Exp_msk1;
1084
/* Now adjust for scale, provided that b != 0. */
1085
if (scale && (b->x[0] || b->x[1])) {
1090
/* We can't shift more than P-1 bits without shifting out a 1. */
1091
assert(0 < scale && scale <= P - 1);
1093
/* The bits shifted out should all be zero. */
1094
assert(b->x[0] == 0);
1100
/* The bits shifted out should all be zero. */
1101
assert(b->x[0] << (32 - scale) == 0);
1102
b->x[0] = (b->x[0] >> scale) | (b->x[1] << (32 - scale));
1107
/* Ensure b is normalized. */
1114
/* Convert a double to a Bigint plus an exponent. Return NULL on failure.
1116
Given a finite nonzero double d, return an odd Bigint b and exponent *e
1117
such that fabs(d) = b * 2**e. On return, *bbits gives the number of
1118
significant bits of b; that is, 2**(*bbits-1) <= b < 2**(*bbits).
1120
If d is zero, then b == 0, *e == -1010, *bbits = 0.
1124
d2b(U *d, int *e, int *bits)
1136
z = word0(d) & Frac_mask;
1137
word0(d) &= 0x7fffffff; /* clear sign bit, which we ignore */
1138
if ((de = (int)(word0(d) >> Exp_shift)))
1140
if ((y = word1(d))) {
1141
if ((k = lo0bits(&y))) {
1142
x[0] = y | z << (32 - k);
1148
b->wds = (x[1] = z) ? 2 : 1;
1158
*e = de - Bias - (P-1) + k;
1162
*e = de - Bias - (P-1) + 1 + k;
1163
*bits = 32*i - hi0bits(x[i-1]);
1168
/* Compute the ratio of two Bigints, as a double. The result may have an
1169
error of up to 2.5 ulps. */
1172
ratio(Bigint *a, Bigint *b)
1177
dval(&da) = b2d(a, &ka);
1178
dval(&db) = b2d(b, &kb);
1179
k = ka - kb + 32*(a->wds - b->wds);
1181
word0(&da) += k*Exp_msk1;
1184
word0(&db) += k*Exp_msk1;
1186
return dval(&da) / dval(&db);
1191
1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1192
1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1197
bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1198
static const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1199
9007199254740992.*9007199254740992.e-256
1200
/* = 2^106 * 1e-256 */
1202
/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1203
/* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1204
#define Scale_Bit 0x10
1213
dshift(Bigint *b, int p2)
1215
int rv = hi0bits(b->x[b->wds-1]) - 4;
1221
/* special case of Bigint division. The quotient is always in the range 0 <=
1222
quotient < 10, and on entry the divisor S is normalized so that its top 4
1223
bits (28--31) are zero and bit 27 is set. */
1226
quorem(Bigint *b, Bigint *S)
1229
ULong *bx, *bxe, q, *sx, *sxe;
1231
ULLong borrow, carry, y, ys;
1233
ULong borrow, carry, y, ys;
1239
/*debug*/ if (b->wds > n)
1240
/*debug*/ Bug("oversize b in quorem");
1248
q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1250
/*debug*/ if (q > 9)
1251
/*debug*/ Bug("oversized quotient in quorem");
1258
ys = *sx++ * (ULLong)q + carry;
1260
y = *bx - (ys & FFFFFFFF) - borrow;
1261
borrow = y >> 32 & (ULong)1;
1262
*bx++ = (ULong)(y & FFFFFFFF);
1265
ys = (si & 0xffff) * q + carry;
1266
zs = (si >> 16) * q + (ys >> 16);
1268
y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
1269
borrow = (y & 0x10000) >> 16;
1270
z = (*bx >> 16) - (zs & 0xffff) - borrow;
1271
borrow = (z & 0x10000) >> 16;
1278
while(--bxe > bx && !*bxe)
1283
if (cmp(b, S) >= 0) {
1293
y = *bx - (ys & FFFFFFFF) - borrow;
1294
borrow = y >> 32 & (ULong)1;
1295
*bx++ = (ULong)(y & FFFFFFFF);
1298
ys = (si & 0xffff) + carry;
1299
zs = (si >> 16) + (ys >> 16);
1301
y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
1302
borrow = (y & 0x10000) >> 16;
1303
z = (*bx >> 16) - (zs & 0xffff) - borrow;
1304
borrow = (z & 0x10000) >> 16;
1312
while(--bxe > bx && !*bxe)
1320
/* sulp(x) is a version of ulp(x) that takes bc.scale into account.
1322
Assuming that x is finite and nonnegative (positive zero is fine
1323
here) and x / 2^bc.scale is exactly representable as a double,
1324
sulp(x) is equivalent to 2^bc.scale * ulp(x / 2^bc.scale). */
1327
sulp(U *x, BCinfo *bc)
1331
if (bc->scale && 2*P + 1 > (int)((word0(x) & Exp_mask) >> Exp_shift)) {
1332
/* rv/2^bc->scale is subnormal */
1333
word0(&u) = (P+2)*Exp_msk1;
1338
assert(word0(x) || word1(x)); /* x != 0.0 */
1343
/* The bigcomp function handles some hard cases for strtod, for inputs
1344
with more than STRTOD_DIGLIM digits. It's called once an initial
1345
estimate for the double corresponding to the input string has
1346
already been obtained by the code in _Py_dg_strtod.
1348
The bigcomp function is only called after _Py_dg_strtod has found a
1349
double value rv such that either rv or rv + 1ulp represents the
1350
correctly rounded value corresponding to the original string. It
1351
determines which of these two values is the correct one by
1352
computing the decimal digits of rv + 0.5ulp and comparing them with
1353
the corresponding digits of s0.
1355
In the following, write dv for the absolute value of the number represented
1356
by the input string.
1360
s0 points to the first significant digit of the input string.
1362
rv is a (possibly scaled) estimate for the closest double value to the
1363
value represented by the original input to _Py_dg_strtod. If
1364
bc->scale is nonzero, then rv/2^(bc->scale) is the approximation to
1367
bc is a struct containing information gathered during the parsing and
1368
estimation steps of _Py_dg_strtod. Description of fields follows:
1370
bc->e0 gives the exponent of the input value, such that dv = (integer
1371
given by the bd->nd digits of s0) * 10**e0
1373
bc->nd gives the total number of significant digits of s0. It will
1376
bc->nd0 gives the number of significant digits of s0 before the
1377
decimal separator. If there's no decimal separator, bc->nd0 ==
1380
bc->scale is the value used to scale rv to avoid doing arithmetic with
1381
subnormal values. It's either 0 or 2*P (=106).
1385
On successful exit, rv/2^(bc->scale) is the closest double to dv.
1387
Returns 0 on success, -1 on failure (e.g., due to a failed malloc call). */
1390
bigcomp(U *rv, const char *s0, BCinfo *bc)
1393
int b2, d2, dd, i, nd, nd0, odd, p2, p5;
1398
b = sd2b(rv, bc->scale, &p2);
1402
/* record whether the lsb of rv/2^(bc->scale) is odd: in the exact halfway
1403
case, this is used for round to even. */
1406
/* left shift b by 1 bit and or a 1 into the least significant bit;
1407
this gives us b * 2**p2 = rv/2^(bc->scale) + 0.5 ulp. */
1420
/* Arrange for convenient computation of quotients:
1421
* shift left if necessary so divisor has 4 leading 0 bits.
1424
d = pow5mult(d, p5);
1431
b = pow5mult(b, -p5);
1446
if ((b2 += i) > 0) {
1453
if ((d2 += i) > 0) {
1461
/* Compare s0 with b/d: set dd to -1, 0, or 1 according as s0 < b/d, s0 ==
1462
* b/d, or s0 > b/d. Here the digits of s0 are thought of as representing
1463
* a number in the range [0.1, 1). */
1470
b = multadd(b, 10, 0);
1475
dd = s0[i < nd0 ? i : i+1] - '0' - quorem(b, d);
1480
if (!b->x[0] && b->wds == 1) {
1486
/* b/d != 0, but digits of s0 exhausted */
1494
if (dd > 0 || (dd == 0 && odd))
1495
dval(rv) += sulp(rv, bc);
1499
/* Return a 'standard' NaN value.
1501
There are exactly two quiet NaNs that don't arise by 'quieting' signaling
1502
NaNs (see IEEE 754-2008, section 6.2.1). If sign == 0, return the one whose
1503
sign bit is cleared. Otherwise, return the one whose sign bit is set.
1507
_Py_dg_stdnan(int sign)
1510
word0(&rv) = NAN_WORD0;
1511
word1(&rv) = NAN_WORD1;
1513
word0(&rv) |= Sign_bit;
1517
/* Return positive or negative infinity, according to the given sign (0 for
1518
* positive infinity, 1 for negative infinity). */
1521
_Py_dg_infinity(int sign)
1524
word0(&rv) = POSINF_WORD0;
1525
word1(&rv) = POSINF_WORD1;
1526
return sign ? -dval(&rv) : dval(&rv);
1530
_Py_dg_strtod(const char *s00, char **se)
1532
int bb2, bb5, bbe, bd2, bd5, bs2, c, dsign, e, e1, error;
1533
int esign, i, j, k, lz, nd, nd0, odd, sign;
1534
const char *s, *s0, *s1;
1536
U aadj2, adj, rv, rv0;
1537
ULong y, z, abs_exp;
1540
Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1544
/* Start parsing. */
1547
/* Parse optional sign, if present. */
1557
/* Skip leading zeros: lz is true iff there were leading zeros. */
1563
/* Point s0 at the first nonzero digit (if any). nd0 will be the position
1564
of the point relative to s0. nd will be the total number of digits
1565
ignoring leading zeros. */
1567
while ('0' <= c && c <= '9')
1571
/* Parse decimal point and following digits. */
1583
while ('0' <= c && c <= '9')
1588
/* Now lz is true if and only if there were leading zero digits, and nd
1589
gives the total number of digits ignoring leading zeros. A valid input
1590
must have at least one digit. */
1597
/* Parse exponent. */
1599
if (c == 'e' || c == 'E') {
1603
/* Exponent sign. */
1613
/* Skip zeros. lz is true iff there are leading zeros. */
1619
/* Get absolute value of the exponent. */
1622
while ('0' <= c && c <= '9') {
1623
abs_exp = 10*abs_exp + (c - '0');
1627
/* abs_exp will be correct modulo 2**32. But 10**9 < 2**32, so if
1628
there are at most 9 significant exponent digits then overflow is
1630
if (s - s1 > 9 || abs_exp > MAX_ABS_EXP)
1631
e = (int)MAX_ABS_EXP;
1637
/* A valid exponent must have at least one digit. */
1642
/* Adjust exponent to take into account position of the point. */
1647
/* Finished parsing. Set se to indicate how far we parsed */
1651
/* If all digits were zero, exit with return value +-0.0. Otherwise,
1652
strip trailing zeros: scan back until we hit a nonzero digit. */
1655
for (i = nd; i > 0; ) {
1657
if (s0[i < nd0 ? i : i+1] != '0') {
1667
/* Summary of parsing results. After parsing, and dealing with zero
1668
* inputs, we have values s0, nd0, nd, e, sign, where:
1670
* - s0 points to the first significant digit of the input string
1672
* - nd is the total number of significant digits (here, and
1673
* below, 'significant digits' means the set of digits of the
1674
* significand of the input that remain after ignoring leading
1675
* and trailing zeros).
1677
* - nd0 indicates the position of the decimal point, if present; it
1678
* satisfies 1 <= nd0 <= nd. The nd significant digits are in
1679
* s0[0:nd0] and s0[nd0+1:nd+1] using the usual Python half-open slice
1680
* notation. (If nd0 < nd, then s0[nd0] contains a '.' character; if
1681
* nd0 == nd, then s0[nd0] could be any non-digit character.)
1683
* - e is the adjusted exponent: the absolute value of the number
1684
* represented by the original input string is n * 10**e, where
1685
* n is the integer represented by the concatenation of
1686
* s0[0:nd0] and s0[nd0+1:nd+1]
1688
* - sign gives the sign of the input: 1 for negative, 0 for positive
1690
* - the first and last significant digits are nonzero
1693
/* put first DBL_DIG+1 digits into integer y and z.
1695
* - y contains the value represented by the first min(9, nd)
1696
* significant digits
1698
* - if nd > 9, z contains the value represented by significant digits
1699
* with indices in [9, min(16, nd)). So y * 10**(min(16, nd) - 9) + z
1700
* gives the value represented by the first min(16, nd) sig. digits.
1705
for (i = 0; i < nd; i++) {
1707
y = 10*y + s0[i < nd0 ? i : i+1] - '0';
1708
else if (i < DBL_DIG+1)
1709
z = 10*z + s0[i < nd0 ? i : i+1] - '0';
1714
k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1717
dval(&rv) = tens[k - 9] * dval(&rv) + z;
1726
if (e <= Ten_pmax) {
1727
dval(&rv) *= tens[e];
1731
if (e <= Ten_pmax + i) {
1732
/* A fancier test would sometimes let us do
1733
* this for larger i values.
1736
dval(&rv) *= tens[i];
1737
dval(&rv) *= tens[e];
1741
else if (e >= -Ten_pmax) {
1742
dval(&rv) /= tens[-e];
1750
/* Get starting approximation = rv * 10**e1 */
1754
dval(&rv) *= tens[i];
1756
if (e1 > DBL_MAX_10_EXP)
1759
for(j = 0; e1 > 1; j++, e1 >>= 1)
1761
dval(&rv) *= bigtens[j];
1762
/* The last multiplication could overflow. */
1763
word0(&rv) -= P*Exp_msk1;
1764
dval(&rv) *= bigtens[j];
1765
if ((z = word0(&rv) & Exp_mask)
1766
> Exp_msk1*(DBL_MAX_EXP+Bias-P))
1768
if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1769
/* set to largest number */
1770
/* (Can't trust DBL_MAX) */
1775
word0(&rv) += P*Exp_msk1;
1779
/* The input decimal value lies in [10**e1, 10**(e1+16)).
1781
If e1 <= -512, underflow immediately.
1782
If e1 <= -256, set bc.scale to 2*P.
1784
So for input value < 1e-256, bc.scale is always set;
1785
for input value >= 1e-240, bc.scale is never set.
1786
For input values in [1e-256, 1e-240), bc.scale may or may
1791
dval(&rv) /= tens[i];
1793
if (e1 >= 1 << n_bigtens)
1797
for(j = 0; e1 > 0; j++, e1 >>= 1)
1799
dval(&rv) *= tinytens[j];
1800
if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
1801
>> Exp_shift)) > 0) {
1802
/* scaled rv is denormal; clear j low bits */
1806
word0(&rv) = (P+2)*Exp_msk1;
1808
word0(&rv) &= 0xffffffff << (j-32);
1811
word1(&rv) &= 0xffffffff << j;
1818
/* Now the hard part -- adjusting rv to the correct value.*/
1820
/* Put digits into bd: true value = bd * 10^e */
1823
bc.nd0 = nd0; /* Only needed if nd > STRTOD_DIGLIM, but done here */
1824
/* to silence an erroneous warning about bc.nd0 */
1825
/* possibly not being initialized. */
1826
if (nd > STRTOD_DIGLIM) {
1827
/* ASSERT(STRTOD_DIGLIM >= 18); 18 == one more than the */
1828
/* minimum number of decimal digits to distinguish double values */
1829
/* in IEEE arithmetic. */
1831
/* Truncate input to 18 significant digits, then discard any trailing
1832
zeros on the result by updating nd, nd0, e and y suitably. (There's
1833
no need to update z; it's not reused beyond this point.) */
1834
for (i = 18; i > 0; ) {
1835
/* scan back until we hit a nonzero digit. significant digit 'i'
1836
is s0[i] if i < nd0, s0[i+1] if i >= nd0. */
1838
if (s0[i < nd0 ? i : i+1] != '0') {
1847
if (nd < 9) { /* must recompute y */
1849
for(i = 0; i < nd0; ++i)
1850
y = 10*y + s0[i] - '0';
1852
y = 10*y + s0[i+1] - '0';
1855
bd0 = s2b(s0, nd0, nd, y);
1859
/* Notation for the comments below. Write:
1861
- dv for the absolute value of the number represented by the original
1862
decimal input string.
1864
- if we've truncated dv, write tdv for the truncated value.
1865
Otherwise, set tdv == dv.
1867
- srv for the quantity rv/2^bc.scale; so srv is the current binary
1868
approximation to tdv (and dv). It should be exactly representable
1869
in an IEEE 754 double.
1874
/* This is the main correction loop for _Py_dg_strtod.
1876
We've got a decimal value tdv, and a floating-point approximation
1877
srv=rv/2^bc.scale to tdv. The aim is to determine whether srv is
1878
close enough (i.e., within 0.5 ulps) to tdv, and to compute a new
1879
approximation if not.
1881
To determine whether srv is close enough to tdv, compute integers
1882
bd, bb and bs proportional to tdv, srv and 0.5 ulp(srv)
1883
respectively, and then use integer arithmetic to determine whether
1884
|tdv - srv| is less than, equal to, or greater than 0.5 ulp(srv).
1887
bd = Balloc(bd0->k);
1893
bb = sd2b(&rv, bc.scale, &bbe); /* srv = bb * 2^bbe */
1899
/* Record whether lsb of bb is odd, in case we need this
1900
for the round-to-even step later. */
1903
/* tdv = bd * 10**e; srv = bb * 2**bbe */
1928
/* At this stage bd5 - bb5 == e == bd2 - bb2 + bbe, bb2 - bs2 == 1,
1931
tdv == bd * 10**e = bd * 2**(bbe - bb2 + bd2) * 5**(bd5 - bb5)
1932
srv == bb * 2**bbe = bb * 2**(bbe - bb2 + bb2)
1933
0.5 ulp(srv) == 2**(bbe-1) = bs * 2**(bbe - bb2 + bs2)
1937
M * tdv = bd * 2**bd2 * 5**bd5
1938
M * srv = bb * 2**bb2 * 5**bb5
1939
M * 0.5 ulp(srv) = bs * 2**bs2 * 5**bb5
1941
for some constant M. (Actually, M == 2**(bb2 - bbe) * 5**bb5, but
1942
this fact is not needed below.)
1945
/* Remove factor of 2**i, where i = min(bb2, bd2, bs2). */
1946
i = bb2 < bd2 ? bb2 : bd2;
1955
/* Scale bb, bd, bs by the appropriate powers of 2 and 5. */
1957
bs = pow5mult(bs, bb5);
1975
bb = lshift(bb, bb2);
1984
bd = pow5mult(bd, bd5);
1993
bd = lshift(bd, bd2);
2002
bs = lshift(bs, bs2);
2011
/* Now bd, bb and bs are scaled versions of tdv, srv and 0.5 ulp(srv),
2012
respectively. Compute the difference |tdv - srv|, and compare
2013
with 0.5 ulp(srv). */
2015
delta = diff(bb, bd);
2016
if (delta == NULL) {
2023
dsign = delta->sign;
2026
if (bc.nd > nd && i <= 0) {
2028
break; /* Must use bigcomp(). */
2030
/* Here rv overestimates the truncated decimal value by at most
2031
0.5 ulp(rv). Hence rv either overestimates the true decimal
2032
value by <= 0.5 ulp(rv), or underestimates it by some small
2033
amount (< 0.1 ulp(rv)); either way, rv is within 0.5 ulps of
2034
the true decimal value, so it's possible to exit.
2036
Exception: if scaled rv is a normal exact power of 2, but not
2037
DBL_MIN, then rv - 0.5 ulp(rv) takes us all the way down to the
2038
next double, so the correctly rounded result is either rv - 0.5
2039
ulp(rv) or rv; in this case, use bigcomp to distinguish. */
2041
if (!word1(&rv) && !(word0(&rv) & Bndry_mask)) {
2042
/* rv can't be 0, since it's an overestimate for some
2043
nonzero value. So rv is a normal power of 2. */
2044
j = (int)(word0(&rv) & Exp_mask) >> Exp_shift;
2045
/* rv / 2^bc.scale = 2^(j - 1023 - bc.scale); use bigcomp if
2046
rv / 2^bc.scale >= 2^-1021. */
2047
if (j - bc.scale >= 2) {
2048
dval(&rv) -= 0.5 * sulp(&rv, &bc);
2049
break; /* Use bigcomp. */
2055
i = -1; /* Discarded digits make delta smaller. */
2060
/* Error is less than half an ulp -- check for
2061
* special case of mantissa a power of two.
2063
if (dsign || word1(&rv) || word0(&rv) & Bndry_mask
2064
|| (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
2068
if (!delta->x[0] && delta->wds <= 1) {
2072
delta = lshift(delta,Log2P);
2073
if (delta == NULL) {
2080
if (cmp(delta, bs) > 0)
2085
/* exactly half-way between */
2087
if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
2090
(y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) ?
2091
(0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
2093
/*boundary case -- increment exponent*/
2094
word0(&rv) = (word0(&rv) & Exp_mask)
2102
else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
2104
/* boundary case -- decrement exponent */
2106
L = word0(&rv) & Exp_mask;
2107
if (L <= (2*P+1)*Exp_msk1) {
2108
if (L > (P+2)*Exp_msk1)
2109
/* round even ==> */
2112
/* rv = smallest denormal */
2118
L = (word0(&rv) & Exp_mask) - Exp_msk1;
2119
word0(&rv) = L | Bndry_mask1;
2120
word1(&rv) = 0xffffffff;
2126
dval(&rv) += sulp(&rv, &bc);
2128
dval(&rv) -= sulp(&rv, &bc);
2135
/* dsign = 1 - dsign; */
2138
if ((aadj = ratio(delta, bs)) <= 2.) {
2141
else if (word1(&rv) || word0(&rv) & Bndry_mask) {
2142
if (word1(&rv) == Tiny1 && !word0(&rv)) {
2151
/* special case -- power of FLT_RADIX to be */
2152
/* rounded down... */
2154
if (aadj < 2./FLT_RADIX)
2155
aadj = 1./FLT_RADIX;
2163
aadj1 = dsign ? aadj : -aadj;
2164
if (Flt_Rounds == 0)
2167
y = word0(&rv) & Exp_mask;
2169
/* Check for overflow */
2171
if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
2172
dval(&rv0) = dval(&rv);
2173
word0(&rv) -= P*Exp_msk1;
2174
adj.d = aadj1 * ulp(&rv);
2176
if ((word0(&rv) & Exp_mask) >=
2177
Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
2178
if (word0(&rv0) == Big0 && word1(&rv0) == Big1) {
2191
word0(&rv) += P*Exp_msk1;
2194
if (bc.scale && y <= 2*P*Exp_msk1) {
2195
if (aadj <= 0x7fffffff) {
2196
if ((z = (ULong)aadj) <= 0)
2199
aadj1 = dsign ? aadj : -aadj;
2201
dval(&aadj2) = aadj1;
2202
word0(&aadj2) += (2*P+1)*Exp_msk1 - y;
2203
aadj1 = dval(&aadj2);
2205
adj.d = aadj1 * ulp(&rv);
2208
z = word0(&rv) & Exp_mask;
2212
/* Can we stop now? */
2215
/* The tolerances below are conservative. */
2216
if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
2217
if (aadj < .4999999 || aadj > .5000001)
2220
else if (aadj < .4999999/FLT_RADIX)
2236
error = bigcomp(&rv, s0, &bc);
2242
word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
2244
dval(&rv) *= dval(&rv0);
2248
return sign ? -dval(&rv) : dval(&rv);
2258
return sign ? -0.0 : 0.0;
2262
/* Can't trust HUGE_VAL */
2263
word0(&rv) = Exp_mask;
2265
return sign ? -dval(&rv) : dval(&rv);
2276
sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (unsigned)i;
2279
r = (int*)Balloc(k);
2283
return (char *)(r+1);
2287
nrv_alloc(char *s, char **rve, int n)
2295
while((*t = *s++)) t++;
2301
/* freedtoa(s) must be used to free values s returned by dtoa
2302
* when MULTIPLE_THREADS is #defined. It should be used in all cases,
2303
* but for consistency with earlier versions of dtoa, it is optional
2304
* when MULTIPLE_THREADS is not defined.
2308
_Py_dg_freedtoa(char *s)
2310
Bigint *b = (Bigint *)((int *)s - 1);
2311
b->maxwds = 1 << (b->k = *(int*)b);
2315
/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
2317
* Inspired by "How to Print Floating-Point Numbers Accurately" by
2318
* Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
2321
* 1. Rather than iterating, we use a simple numeric overestimate
2322
* to determine k = floor(log10(d)). We scale relevant
2323
* quantities using O(log2(k)) rather than O(k) multiplications.
2324
* 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
2325
* try to generate digits strictly left to right. Instead, we
2326
* compute with fewer bits and propagate the carry if necessary
2327
* when rounding the final digit up. This is often faster.
2328
* 3. Under the assumption that input will be rounded nearest,
2329
* mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
2330
* That is, we allow equality in stopping tests when the
2331
* round-nearest rule will give the same floating-point value
2332
* as would satisfaction of the stopping test with strict
2334
* 4. We remove common factors of powers of 2 from relevant
2336
* 5. When converting floating-point integers less than 1e16,
2337
* we use floating-point arithmetic rather than resorting
2338
* to multiple-precision integers.
2339
* 6. When asked to produce fewer than 15 digits, we first try
2340
* to get by with floating-point arithmetic; we resort to
2341
* multiple-precision integer arithmetic only if we cannot
2342
* guarantee that the floating-point calculation has given
2343
* the correctly rounded result. For k requested digits and
2344
* "uniformly" distributed input, the probability is
2345
* something like 10^(k-15) that we must resort to the Long
2349
/* Additional notes (METD): (1) returns NULL on failure. (2) to avoid memory
2350
leakage, a successful call to _Py_dg_dtoa should always be matched by a
2351
call to _Py_dg_freedtoa. */
2354
_Py_dg_dtoa(double dd, int mode, int ndigits,
2355
int *decpt, int *sign, char **rve)
2357
/* Arguments ndigits, decpt, sign are similar to those
2358
of ecvt and fcvt; trailing zeros are suppressed from
2359
the returned string. If not null, *rve is set to point
2360
to the end of the return value. If d is +-Infinity or NaN,
2361
then *decpt is set to 9999.
2364
0 ==> shortest string that yields d when read in
2365
and rounded to nearest.
2366
1 ==> like 0, but with Steele & White stopping rule;
2367
e.g. with IEEE P754 arithmetic , mode 0 gives
2368
1e23 whereas mode 1 gives 9.999999999999999e22.
2369
2 ==> max(1,ndigits) significant digits. This gives a
2370
return value similar to that of ecvt, except
2371
that trailing zeros are suppressed.
2372
3 ==> through ndigits past the decimal point. This
2373
gives a return value similar to that from fcvt,
2374
except that trailing zeros are suppressed, and
2375
ndigits can be negative.
2376
4,5 ==> similar to 2 and 3, respectively, but (in
2377
round-nearest mode) with the tests of mode 0 to
2378
possibly return a shorter string that rounds to d.
2379
With IEEE arithmetic and compilation with
2380
-DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
2381
as modes 2 and 3 when FLT_ROUNDS != 1.
2382
6-9 ==> Debugging modes similar to mode - 4: don't try
2383
fast floating-point estimate (if applicable).
2385
Values of mode other than 0-9 are treated as mode 0.
2387
Sufficient space is allocated to the return value
2388
to hold the suppressed trailing zeros.
2391
int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
2392
j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
2393
spec_case, try_quick;
2397
Bigint *b, *b1, *delta, *mlo, *mhi, *S;
2402
/* set pointers to NULL, to silence gcc compiler warnings and make
2403
cleanup easier on error */
2408
if (word0(&u) & Sign_bit) {
2409
/* set sign for everything, including 0's and NaNs */
2411
word0(&u) &= ~Sign_bit; /* clear sign bit */
2416
/* quick return for Infinities, NaNs and zeros */
2417
if ((word0(&u) & Exp_mask) == Exp_mask)
2419
/* Infinity or NaN */
2421
if (!word1(&u) && !(word0(&u) & 0xfffff))
2422
return nrv_alloc("Infinity", rve, 8);
2423
return nrv_alloc("NaN", rve, 3);
2427
return nrv_alloc("0", rve, 1);
2430
/* compute k = floor(log10(d)). The computation may leave k
2431
one too large, but should never leave k too small. */
2432
b = d2b(&u, &be, &bbits);
2435
if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
2436
dval(&d2) = dval(&u);
2437
word0(&d2) &= Frac_mask1;
2438
word0(&d2) |= Exp_11;
2440
/* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2441
* log10(x) = log(x) / log(10)
2442
* ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2443
* log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2445
* This suggests computing an approximation k to log10(d) by
2447
* k = (i - Bias)*0.301029995663981
2448
* + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2450
* We want k to be too large rather than too small.
2451
* The error in the first-order Taylor series approximation
2452
* is in our favor, so we just round up the constant enough
2453
* to compensate for any error in the multiplication of
2454
* (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2455
* and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2456
* adding 1e-13 to the constant term more than suffices.
2457
* Hence we adjust the constant term to 0.1760912590558.
2458
* (We could get a more accurate k by invoking log10,
2459
* but this is probably not worthwhile.)
2466
/* d is denormalized */
2468
i = bbits + be + (Bias + (P-1) - 1);
2469
x = i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32)
2470
: word1(&u) << (32 - i);
2472
word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
2473
i -= (Bias + (P-1) - 1) + 1;
2476
ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 +
2477
i*0.301029995663981;
2479
if (ds < 0. && ds != k)
2480
k--; /* want k = floor(ds) */
2482
if (k >= 0 && k <= Ten_pmax) {
2483
if (dval(&u) < tens[k])
2506
if (mode < 0 || mode > 9)
2516
ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
2517
/* silence erroneous "gcc -Wall" warning. */
2530
ilim = ilim1 = i = ndigits;
2536
i = ndigits + k + 1;
2548
if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2550
/* Try to get by with floating-point arithmetic. */
2553
dval(&d2) = dval(&u);
2556
ieps = 2; /* conservative */
2561
/* prevent overflows */
2563
dval(&u) /= bigtens[n_bigtens-1];
2566
for(; j; j >>= 1, i++)
2573
else if ((j1 = -k)) {
2574
dval(&u) *= tens[j1 & 0xf];
2575
for(j = j1 >> 4; j; j >>= 1, i++)
2578
dval(&u) *= bigtens[i];
2581
if (k_check && dval(&u) < 1. && ilim > 0) {
2589
dval(&eps) = ieps*dval(&u) + 7.;
2590
word0(&eps) -= (P-1)*Exp_msk1;
2594
if (dval(&u) > dval(&eps))
2596
if (dval(&u) < -dval(&eps))
2601
/* Use Steele & White method of only
2602
* generating digits needed.
2604
dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
2608
*s++ = '0' + (int)L;
2609
if (dval(&u) < dval(&eps))
2611
if (1. - dval(&u) < dval(&eps))
2620
/* Generate ilim digits, then fix them up. */
2621
dval(&eps) *= tens[ilim-1];
2622
for(i = 1;; i++, dval(&u) *= 10.) {
2623
L = (Long)(dval(&u));
2624
if (!(dval(&u) -= L))
2626
*s++ = '0' + (int)L;
2628
if (dval(&u) > 0.5 + dval(&eps))
2630
else if (dval(&u) < 0.5 - dval(&eps)) {
2641
dval(&u) = dval(&d2);
2646
/* Do we have a "small" integer? */
2648
if (be >= 0 && k <= Int_max) {
2651
if (ndigits < 0 && ilim <= 0) {
2653
if (ilim < 0 || dval(&u) <= 5*ds)
2657
for(i = 1;; i++, dval(&u) *= 10.) {
2658
L = (Long)(dval(&u) / ds);
2660
*s++ = '0' + (int)L;
2665
dval(&u) += dval(&u);
2666
if (dval(&u) > ds || (dval(&u) == ds && L & 1)) {
2686
denorm ? be + (Bias + (P-1) - 1 + 1) :
2694
if (m2 > 0 && s2 > 0) {
2695
i = m2 < s2 ? m2 : s2;
2703
mhi = pow5mult(mhi, m5);
2712
if ((j = b5 - m5)) {
2719
b = pow5mult(b, b5);
2728
S = pow5mult(S, s5);
2733
/* Check for special case that d is a normalized power of 2. */
2736
if ((mode < 2 || leftright)
2738
if (!word1(&u) && !(word0(&u) & Bndry_mask)
2739
&& word0(&u) & (Exp_mask & ~Exp_msk1)
2741
/* The special case */
2748
/* Arrange for convenient computation of quotients:
2749
* shift left if necessary so divisor has 4 leading 0 bits.
2751
* Perhaps we should just compute leading 28 bits of S once
2752
* and for all and pass them and a shift to quorem, so it
2753
* can do shifts and ors to compute the numerator for q.
2773
b = multadd(b, 10, 0); /* we botched the k estimate */
2777
mhi = multadd(mhi, 10, 0);
2784
if (ilim <= 0 && (mode == 3 || mode == 5)) {
2786
/* no digits, fcvt style */
2792
S = multadd(S, 5, 0);
2805
mhi = lshift(mhi, m2);
2810
/* Compute mlo -- check for special case
2811
* that d is a normalized power of 2.
2816
mhi = Balloc(mhi->k);
2820
mhi = lshift(mhi, Log2P);
2826
dig = quorem(b,S) + '0';
2827
/* Do we yet have the shortest decimal string
2828
* that will round to d?
2831
delta = diff(S, mhi);
2834
j1 = delta->sign ? 1 : cmp(b, delta);
2836
if (j1 == 0 && mode != 1 && !(word1(&u) & 1)
2845
if (j < 0 || (j == 0 && mode != 1
2848
if (!b->x[0] && b->wds <= 1) {
2856
if ((j1 > 0 || (j1 == 0 && dig & 1))
2865
if (dig == '9') { /* possible if i == 1 */
2876
b = multadd(b, 10, 0);
2880
mlo = mhi = multadd(mhi, 10, 0);
2885
mlo = multadd(mlo, 10, 0);
2888
mhi = multadd(mhi, 10, 0);
2896
*s++ = dig = quorem(b,S) + '0';
2897
if (!b->x[0] && b->wds <= 1) {
2902
b = multadd(b, 10, 0);
2907
/* Round off last digit */
2913
if (j > 0 || (j == 0 && dig & 1)) {
2930
if (mlo && mlo != mhi)
2944
if (mlo && mlo != mhi)
2951
_Py_dg_freedtoa(s0);
2958
#endif /* PY_NO_SHORT_FLOAT_REPR */
added
removed
Python/dtoa.c
