1 /****************************************************************
2 *
3 * The author of this software is David M. Gay.
4 *
5 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
6 *
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.
12 *
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.
17 *
18 ***************************************************************/
19
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.
24 *
25 * Please remember to check http://www.netlib.org/fp regularly (and especially
26 * before any Python release) for bugfixes and updates.
27 *
28 * The major modifications from Gay's original code are as follows:
29 *
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
34 * been removed.
35 *
36 * 1. We use PyMem_Malloc and PyMem_Free in place of malloc and free.
37 *
38 * 2. The public functions strtod, dtoa and freedtoa all now have
39 * a _Py_dg_ prefix.
40 *
41 * 3. Instead of assuming that PyMem_Malloc always succeeds, we thread
42 * PyMem_Malloc failures through the code. The functions
43 *
44 * Balloc, multadd, s2b, i2b, mult, pow5mult, lshift, diff, d2b
45 *
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.
52 *
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.
56 *
57 * 5. The code has been reformatted to better fit with Python's
58 * C style guide (PEP 7).
59 *
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 <=
62 * Kmax.
63 *
64 * 7. _Py_dg_strtod has been modified so that it doesn't accept strings with
65 * leading whitespace.
66 *
67 * 8. A corner case where _Py_dg_dtoa didn't strip trailing zeros has been
68 * fixed. (bugs.python.org/issue40780)
69 *
70 ***************************************************************/
71
72 /* Please send bug reports for the original dtoa.c code to David M. Gay (dmg
73 * at acm dot org, with " at " changed at "@" and " dot " changed to ".").
74 * Please report bugs for this modified version using the Python issue tracker
75 * as detailed at (https://devguide.python.org/triage/issue-tracker/). */
76
77 /* On a machine with IEEE extended-precision registers, it is
78 * necessary to specify double-precision (53-bit) rounding precision
79 * before invoking strtod or dtoa. If the machine uses (the equivalent
80 * of) Intel 80x87 arithmetic, the call
81 * _control87(PC_53, MCW_PC);
82 * does this with many compilers. Whether this or another call is
83 * appropriate depends on the compiler; for this to work, it may be
84 * necessary to #include "float.h" or another system-dependent header
85 * file.
86 */
87
88 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
89 *
90 * This strtod returns a nearest machine number to the input decimal
91 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
92 * broken by the IEEE round-even rule. Otherwise ties are broken by
93 * biased rounding (add half and chop).
94 *
95 * Inspired loosely by William D. Clinger's paper "How to Read Floating
96 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
97 *
98 * Modifications:
99 *
100 * 1. We only require IEEE, IBM, or VAX double-precision
101 * arithmetic (not IEEE double-extended).
102 * 2. We get by with floating-point arithmetic in a case that
103 * Clinger missed -- when we're computing d * 10^n
104 * for a small integer d and the integer n is not too
105 * much larger than 22 (the maximum integer k for which
106 * we can represent 10^k exactly), we may be able to
107 * compute (d*10^k) * 10^(e-k) with just one roundoff.
108 * 3. Rather than a bit-at-a-time adjustment of the binary
109 * result in the hard case, we use floating-point
110 * arithmetic to determine the adjustment to within
111 * one bit; only in really hard cases do we need to
112 * compute a second residual.
113 * 4. Because of 3., we don't need a large table of powers of 10
114 * for ten-to-e (just some small tables, e.g. of 10^k
115 * for 0 <= k <= 22).
116 */
117
118 /* Linking of Python's #defines to Gay's #defines starts here. */
119
120 #include "Python.h"
121 #include "pycore_dtoa.h" // _PY_SHORT_FLOAT_REPR
122 #include "pycore_pystate.h" // _PyInterpreterState_GET()
123 #include <stdlib.h> // exit()
124
125 /* if _PY_SHORT_FLOAT_REPR == 0, then don't even try to compile
126 the following code */
127 #if _PY_SHORT_FLOAT_REPR == 1
128
129 #include "float.h"
130
131 #define MALLOC PyMem_Malloc
132 #define FREE PyMem_Free
133
134 /* This code should also work for ARM mixed-endian format on little-endian
135 machines, where doubles have byte order 45670123 (in increasing address
136 order, 0 being the least significant byte). */
137 #ifdef DOUBLE_IS_LITTLE_ENDIAN_IEEE754
138 # define IEEE_8087
139 #endif
140 #if defined(DOUBLE_IS_BIG_ENDIAN_IEEE754) || \
141 defined(DOUBLE_IS_ARM_MIXED_ENDIAN_IEEE754)
142 # define IEEE_MC68k
143 #endif
144 #if defined(IEEE_8087) + defined(IEEE_MC68k) != 1
145 #error "Exactly one of IEEE_8087 or IEEE_MC68k should be defined."
146 #endif
147
148 /* The code below assumes that the endianness of integers matches the
149 endianness of the two 32-bit words of a double. Check this. */
150 #if defined(WORDS_BIGENDIAN) && (defined(DOUBLE_IS_LITTLE_ENDIAN_IEEE754) || \
151 defined(DOUBLE_IS_ARM_MIXED_ENDIAN_IEEE754))
152 #error "doubles and ints have incompatible endianness"
153 #endif
154
155 #if !defined(WORDS_BIGENDIAN) && defined(DOUBLE_IS_BIG_ENDIAN_IEEE754)
156 #error "doubles and ints have incompatible endianness"
157 #endif
158
159
160 // ULong is defined in pycore_dtoa.h.
161 typedef int32_t Long;
162 typedef uint64_t ULLong;
163
164 #undef DEBUG
165 #ifdef Py_DEBUG
166 #define DEBUG
167 #endif
168
169 /* End Python #define linking */
170
171 #ifdef DEBUG
172 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
173 #endif
174
175 typedef union { double d; ULong L[2]; } U;
176
177 #ifdef IEEE_8087
178 #define word0(x) (x)->L[1]
179 #define word1(x) (x)->L[0]
180 #else
181 #define word0(x) (x)->L[0]
182 #define word1(x) (x)->L[1]
183 #endif
184 #define dval(x) (x)->d
185
186 #ifndef STRTOD_DIGLIM
187 #define STRTOD_DIGLIM 40
188 #endif
189
190 /* maximum permitted exponent value for strtod; exponents larger than
191 MAX_ABS_EXP in absolute value get truncated to +-MAX_ABS_EXP. MAX_ABS_EXP
192 should fit into an int. */
193 #ifndef MAX_ABS_EXP
194 #define MAX_ABS_EXP 1100000000U
195 #endif
196 /* Bound on length of pieces of input strings in _Py_dg_strtod; specifically,
197 this is used to bound the total number of digits ignoring leading zeros and
198 the number of digits that follow the decimal point. Ideally, MAX_DIGITS
199 should satisfy MAX_DIGITS + 400 < MAX_ABS_EXP; that ensures that the
200 exponent clipping in _Py_dg_strtod can't affect the value of the output. */
201 #ifndef MAX_DIGITS
202 #define MAX_DIGITS 1000000000U
203 #endif
204
205 /* Guard against trying to use the above values on unusual platforms with ints
206 * of width less than 32 bits. */
207 #if MAX_ABS_EXP > INT_MAX
208 #error "MAX_ABS_EXP should fit in an int"
209 #endif
210 #if MAX_DIGITS > INT_MAX
211 #error "MAX_DIGITS should fit in an int"
212 #endif
213
214 /* The following definition of Storeinc is appropriate for MIPS processors.
215 * An alternative that might be better on some machines is
216 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
217 */
218 #if defined(IEEE_8087)
219 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
220 ((unsigned short *)a)[0] = (unsigned short)c, a++)
221 #else
222 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
223 ((unsigned short *)a)[1] = (unsigned short)c, a++)
224 #endif
225
226 /* #define P DBL_MANT_DIG */
227 /* Ten_pmax = floor(P*log(2)/log(5)) */
228 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
229 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
230 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
231
232 #define Exp_shift 20
233 #define Exp_shift1 20
234 #define Exp_msk1 0x100000
235 #define Exp_msk11 0x100000
236 #define Exp_mask 0x7ff00000
237 #define P 53
238 #define Nbits 53
239 #define Bias 1023
240 #define Emax 1023
241 #define Emin (-1022)
242 #define Etiny (-1074) /* smallest denormal is 2**Etiny */
243 #define Exp_1 0x3ff00000
244 #define Exp_11 0x3ff00000
245 #define Ebits 11
246 #define Frac_mask 0xfffff
247 #define Frac_mask1 0xfffff
248 #define Ten_pmax 22
249 #define Bletch 0x10
250 #define Bndry_mask 0xfffff
251 #define Bndry_mask1 0xfffff
252 #define Sign_bit 0x80000000
253 #define Log2P 1
254 #define Tiny0 0
255 #define Tiny1 1
256 #define Quick_max 14
257 #define Int_max 14
258
259 #ifndef Flt_Rounds
260 #ifdef FLT_ROUNDS
261 #define Flt_Rounds FLT_ROUNDS
262 #else
263 #define Flt_Rounds 1
264 #endif
265 #endif /*Flt_Rounds*/
266
267 #define Rounding Flt_Rounds
268
269 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
270 #define Big1 0xffffffff
271
272 /* Bits of the representation of positive infinity. */
273
274 #define POSINF_WORD0 0x7ff00000
275 #define POSINF_WORD1 0
276
277 /* struct BCinfo is used to pass information from _Py_dg_strtod to bigcomp */
278
279 typedef struct BCinfo BCinfo;
280 struct
281 BCinfo {
282 int e0, nd, nd0, scale;
283 };
284
285 #define FFFFFFFF 0xffffffffUL
286
287 /* struct Bigint is used to represent arbitrary-precision integers. These
288 integers are stored in sign-magnitude format, with the magnitude stored as
289 an array of base 2**32 digits. Bigints are always normalized: if x is a
290 Bigint then x->wds >= 1, and either x->wds == 1 or x[wds-1] is nonzero.
291
292 The Bigint fields are as follows:
293
294 - next is a header used by Balloc and Bfree to keep track of lists
295 of freed Bigints; it's also used for the linked list of
296 powers of 5 of the form 5**2**i used by pow5mult.
297 - k indicates which pool this Bigint was allocated from
298 - maxwds is the maximum number of words space was allocated for
299 (usually maxwds == 2**k)
300 - sign is 1 for negative Bigints, 0 for positive. The sign is unused
301 (ignored on inputs, set to 0 on outputs) in almost all operations
302 involving Bigints: a notable exception is the diff function, which
303 ignores signs on inputs but sets the sign of the output correctly.
304 - wds is the actual number of significant words
305 - x contains the vector of words (digits) for this Bigint, from least
306 significant (x[0]) to most significant (x[wds-1]).
307 */
308
309 // struct Bigint is defined in pycore_dtoa.h.
310 typedef struct Bigint Bigint;
311
312 #if !defined(Py_GIL_DISABLED) && !defined(Py_USING_MEMORY_DEBUGGER)
313
314 /* Memory management: memory is allocated from, and returned to, Kmax+1 pools
315 of memory, where pool k (0 <= k <= Kmax) is for Bigints b with b->maxwds ==
316 1 << k. These pools are maintained as linked lists, with freelist[k]
317 pointing to the head of the list for pool k.
318
319 On allocation, if there's no free slot in the appropriate pool, MALLOC is
320 called to get more memory. This memory is not returned to the system until
321 Python quits. There's also a private memory pool that's allocated from
322 in preference to using MALLOC.
323
324 For Bigints with more than (1 << Kmax) digits (which implies at least 1233
325 decimal digits), memory is directly allocated using MALLOC, and freed using
326 FREE.
327
328 XXX: it would be easy to bypass this memory-management system and
329 translate each call to Balloc into a call to PyMem_Malloc, and each
330 Bfree to PyMem_Free. Investigate whether this has any significant
331 performance on impact. */
332
333 #define freelist interp->dtoa.freelist
334 #define private_mem interp->dtoa.preallocated
335 #define pmem_next interp->dtoa.preallocated_next
336
337 /* Allocate space for a Bigint with up to 1<<k digits */
338
339 static Bigint *
Balloc(int k)340 Balloc(int k)
341 {
342 int x;
343 Bigint *rv;
344 unsigned int len;
345 PyInterpreterState *interp = _PyInterpreterState_GET();
346
347 if (k <= Bigint_Kmax && (rv = freelist[k]))
348 freelist[k] = rv->next;
349 else {
350 x = 1 << k;
351 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
352 /sizeof(double);
353 if (k <= Bigint_Kmax &&
354 pmem_next - private_mem + len <= (Py_ssize_t)Bigint_PREALLOC_SIZE
355 ) {
356 rv = (Bigint*)pmem_next;
357 pmem_next += len;
358 }
359 else {
360 rv = (Bigint*)MALLOC(len*sizeof(double));
361 if (rv == NULL)
362 return NULL;
363 }
364 rv->k = k;
365 rv->maxwds = x;
366 }
367 rv->sign = rv->wds = 0;
368 return rv;
369 }
370
371 /* Free a Bigint allocated with Balloc */
372
373 static void
Bfree(Bigint * v)374 Bfree(Bigint *v)
375 {
376 if (v) {
377 if (v->k > Bigint_Kmax)
378 FREE((void*)v);
379 else {
380 PyInterpreterState *interp = _PyInterpreterState_GET();
381 v->next = freelist[v->k];
382 freelist[v->k] = v;
383 }
384 }
385 }
386
387 #undef pmem_next
388 #undef private_mem
389 #undef freelist
390
391 #else
392
393 /* Alternative versions of Balloc and Bfree that use PyMem_Malloc and
394 PyMem_Free directly in place of the custom memory allocation scheme above.
395 These are provided for the benefit of memory debugging tools like
396 Valgrind. */
397
398 /* Allocate space for a Bigint with up to 1<<k digits */
399
400 static Bigint *
Balloc(int k)401 Balloc(int k)
402 {
403 int x;
404 Bigint *rv;
405 unsigned int len;
406
407 x = 1 << k;
408 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
409 /sizeof(double);
410
411 rv = (Bigint*)MALLOC(len*sizeof(double));
412 if (rv == NULL)
413 return NULL;
414
415 rv->k = k;
416 rv->maxwds = x;
417 rv->sign = rv->wds = 0;
418 return rv;
419 }
420
421 /* Free a Bigint allocated with Balloc */
422
423 static void
Bfree(Bigint * v)424 Bfree(Bigint *v)
425 {
426 if (v) {
427 FREE((void*)v);
428 }
429 }
430
431 #endif /* !defined(Py_GIL_DISABLED) && !defined(Py_USING_MEMORY_DEBUGGER) */
432
433 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
434 y->wds*sizeof(Long) + 2*sizeof(int))
435
436 /* Multiply a Bigint b by m and add a. Either modifies b in place and returns
437 a pointer to the modified b, or Bfrees b and returns a pointer to a copy.
438 On failure, return NULL. In this case, b will have been already freed. */
439
440 static Bigint *
multadd(Bigint * b,int m,int a)441 multadd(Bigint *b, int m, int a) /* multiply by m and add a */
442 {
443 int i, wds;
444 ULong *x;
445 ULLong carry, y;
446 Bigint *b1;
447
448 wds = b->wds;
449 x = b->x;
450 i = 0;
451 carry = a;
452 do {
453 y = *x * (ULLong)m + carry;
454 carry = y >> 32;
455 *x++ = (ULong)(y & FFFFFFFF);
456 }
457 while(++i < wds);
458 if (carry) {
459 if (wds >= b->maxwds) {
460 b1 = Balloc(b->k+1);
461 if (b1 == NULL){
462 Bfree(b);
463 return NULL;
464 }
465 Bcopy(b1, b);
466 Bfree(b);
467 b = b1;
468 }
469 b->x[wds++] = (ULong)carry;
470 b->wds = wds;
471 }
472 return b;
473 }
474
475 /* convert a string s containing nd decimal digits (possibly containing a
476 decimal separator at position nd0, which is ignored) to a Bigint. This
477 function carries on where the parsing code in _Py_dg_strtod leaves off: on
478 entry, y9 contains the result of converting the first 9 digits. Returns
479 NULL on failure. */
480
481 static Bigint *
s2b(const char * s,int nd0,int nd,ULong y9)482 s2b(const char *s, int nd0, int nd, ULong y9)
483 {
484 Bigint *b;
485 int i, k;
486 Long x, y;
487
488 x = (nd + 8) / 9;
489 for(k = 0, y = 1; x > y; y <<= 1, k++) ;
490 b = Balloc(k);
491 if (b == NULL)
492 return NULL;
493 b->x[0] = y9;
494 b->wds = 1;
495
496 if (nd <= 9)
497 return b;
498
499 s += 9;
500 for (i = 9; i < nd0; i++) {
501 b = multadd(b, 10, *s++ - '0');
502 if (b == NULL)
503 return NULL;
504 }
505 s++;
506 for(; i < nd; i++) {
507 b = multadd(b, 10, *s++ - '0');
508 if (b == NULL)
509 return NULL;
510 }
511 return b;
512 }
513
514 /* count leading 0 bits in the 32-bit integer x. */
515
516 static int
hi0bits(ULong x)517 hi0bits(ULong x)
518 {
519 int k = 0;
520
521 if (!(x & 0xffff0000)) {
522 k = 16;
523 x <<= 16;
524 }
525 if (!(x & 0xff000000)) {
526 k += 8;
527 x <<= 8;
528 }
529 if (!(x & 0xf0000000)) {
530 k += 4;
531 x <<= 4;
532 }
533 if (!(x & 0xc0000000)) {
534 k += 2;
535 x <<= 2;
536 }
537 if (!(x & 0x80000000)) {
538 k++;
539 if (!(x & 0x40000000))
540 return 32;
541 }
542 return k;
543 }
544
545 /* count trailing 0 bits in the 32-bit integer y, and shift y right by that
546 number of bits. */
547
548 static int
lo0bits(ULong * y)549 lo0bits(ULong *y)
550 {
551 int k;
552 ULong x = *y;
553
554 if (x & 7) {
555 if (x & 1)
556 return 0;
557 if (x & 2) {
558 *y = x >> 1;
559 return 1;
560 }
561 *y = x >> 2;
562 return 2;
563 }
564 k = 0;
565 if (!(x & 0xffff)) {
566 k = 16;
567 x >>= 16;
568 }
569 if (!(x & 0xff)) {
570 k += 8;
571 x >>= 8;
572 }
573 if (!(x & 0xf)) {
574 k += 4;
575 x >>= 4;
576 }
577 if (!(x & 0x3)) {
578 k += 2;
579 x >>= 2;
580 }
581 if (!(x & 1)) {
582 k++;
583 x >>= 1;
584 if (!x)
585 return 32;
586 }
587 *y = x;
588 return k;
589 }
590
591 /* convert a small nonnegative integer to a Bigint */
592
593 static Bigint *
i2b(int i)594 i2b(int i)
595 {
596 Bigint *b;
597
598 b = Balloc(1);
599 if (b == NULL)
600 return NULL;
601 b->x[0] = i;
602 b->wds = 1;
603 return b;
604 }
605
606 /* multiply two Bigints. Returns a new Bigint, or NULL on failure. Ignores
607 the signs of a and b. */
608
609 static Bigint *
mult(Bigint * a,Bigint * b)610 mult(Bigint *a, Bigint *b)
611 {
612 Bigint *c;
613 int k, wa, wb, wc;
614 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
615 ULong y;
616 ULLong carry, z;
617
618 if ((!a->x[0] && a->wds == 1) || (!b->x[0] && b->wds == 1)) {
619 c = Balloc(0);
620 if (c == NULL)
621 return NULL;
622 c->wds = 1;
623 c->x[0] = 0;
624 return c;
625 }
626
627 if (a->wds < b->wds) {
628 c = a;
629 a = b;
630 b = c;
631 }
632 k = a->k;
633 wa = a->wds;
634 wb = b->wds;
635 wc = wa + wb;
636 if (wc > a->maxwds)
637 k++;
638 c = Balloc(k);
639 if (c == NULL)
640 return NULL;
641 for(x = c->x, xa = x + wc; x < xa; x++)
642 *x = 0;
643 xa = a->x;
644 xae = xa + wa;
645 xb = b->x;
646 xbe = xb + wb;
647 xc0 = c->x;
648 for(; xb < xbe; xc0++) {
649 if ((y = *xb++)) {
650 x = xa;
651 xc = xc0;
652 carry = 0;
653 do {
654 z = *x++ * (ULLong)y + *xc + carry;
655 carry = z >> 32;
656 *xc++ = (ULong)(z & FFFFFFFF);
657 }
658 while(x < xae);
659 *xc = (ULong)carry;
660 }
661 }
662 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
663 c->wds = wc;
664 return c;
665 }
666
667 #ifndef Py_USING_MEMORY_DEBUGGER
668
669 /* multiply the Bigint b by 5**k. Returns a pointer to the result, or NULL on
670 failure; if the returned pointer is distinct from b then the original
671 Bigint b will have been Bfree'd. Ignores the sign of b. */
672
673 static Bigint *
pow5mult(Bigint * b,int k)674 pow5mult(Bigint *b, int k)
675 {
676 Bigint *b1, *p5, **p5s;
677 int i;
678 static const int p05[3] = { 5, 25, 125 };
679
680 // For double-to-string conversion, the maximum value of k is limited by
681 // DBL_MAX_10_EXP (308), the maximum decimal base-10 exponent for binary64.
682 // For string-to-double conversion, the extreme case is constrained by our
683 // hardcoded exponent limit before we underflow of -512, adjusted by
684 // STRTOD_DIGLIM-DBL_DIG-1, giving a maximum of k=535.
685 assert(0 <= k && k < 1024);
686
687 if ((i = k & 3)) {
688 b = multadd(b, p05[i-1], 0);
689 if (b == NULL)
690 return NULL;
691 }
692
693 if (!(k >>= 2))
694 return b;
695 PyInterpreterState *interp = _PyInterpreterState_GET();
696 p5s = interp->dtoa.p5s;
697 for(;;) {
698 assert(p5s != interp->dtoa.p5s + Bigint_Pow5size);
699 p5 = *p5s;
700 p5s++;
701 if (k & 1) {
702 b1 = mult(b, p5);
703 Bfree(b);
704 b = b1;
705 if (b == NULL)
706 return NULL;
707 }
708 if (!(k >>= 1))
709 break;
710 }
711 return b;
712 }
713
714 #else
715
716 /* Version of pow5mult that doesn't cache powers of 5. Provided for
717 the benefit of memory debugging tools like Valgrind. */
718
719 static Bigint *
pow5mult(Bigint * b,int k)720 pow5mult(Bigint *b, int k)
721 {
722 Bigint *b1, *p5, *p51;
723 int i;
724 static const int p05[3] = { 5, 25, 125 };
725
726 if ((i = k & 3)) {
727 b = multadd(b, p05[i-1], 0);
728 if (b == NULL)
729 return NULL;
730 }
731
732 if (!(k >>= 2))
733 return b;
734 p5 = i2b(625);
735 if (p5 == NULL) {
736 Bfree(b);
737 return NULL;
738 }
739
740 for(;;) {
741 if (k & 1) {
742 b1 = mult(b, p5);
743 Bfree(b);
744 b = b1;
745 if (b == NULL) {
746 Bfree(p5);
747 return NULL;
748 }
749 }
750 if (!(k >>= 1))
751 break;
752 p51 = mult(p5, p5);
753 Bfree(p5);
754 p5 = p51;
755 if (p5 == NULL) {
756 Bfree(b);
757 return NULL;
758 }
759 }
760 Bfree(p5);
761 return b;
762 }
763
764 #endif /* Py_USING_MEMORY_DEBUGGER */
765
766 /* shift a Bigint b left by k bits. Return a pointer to the shifted result,
767 or NULL on failure. If the returned pointer is distinct from b then the
768 original b will have been Bfree'd. Ignores the sign of b. */
769
770 static Bigint *
lshift(Bigint * b,int k)771 lshift(Bigint *b, int k)
772 {
773 int i, k1, n, n1;
774 Bigint *b1;
775 ULong *x, *x1, *xe, z;
776
777 if (!k || (!b->x[0] && b->wds == 1))
778 return b;
779
780 n = k >> 5;
781 k1 = b->k;
782 n1 = n + b->wds + 1;
783 for(i = b->maxwds; n1 > i; i <<= 1)
784 k1++;
785 b1 = Balloc(k1);
786 if (b1 == NULL) {
787 Bfree(b);
788 return NULL;
789 }
790 x1 = b1->x;
791 for(i = 0; i < n; i++)
792 *x1++ = 0;
793 x = b->x;
794 xe = x + b->wds;
795 if (k &= 0x1f) {
796 k1 = 32 - k;
797 z = 0;
798 do {
799 *x1++ = *x << k | z;
800 z = *x++ >> k1;
801 }
802 while(x < xe);
803 if ((*x1 = z))
804 ++n1;
805 }
806 else do
807 *x1++ = *x++;
808 while(x < xe);
809 b1->wds = n1 - 1;
810 Bfree(b);
811 return b1;
812 }
813
814 /* Do a three-way compare of a and b, returning -1 if a < b, 0 if a == b and
815 1 if a > b. Ignores signs of a and b. */
816
817 static int
cmp(Bigint * a,Bigint * b)818 cmp(Bigint *a, Bigint *b)
819 {
820 ULong *xa, *xa0, *xb, *xb0;
821 int i, j;
822
823 i = a->wds;
824 j = b->wds;
825 #ifdef DEBUG
826 if (i > 1 && !a->x[i-1])
827 Bug("cmp called with a->x[a->wds-1] == 0");
828 if (j > 1 && !b->x[j-1])
829 Bug("cmp called with b->x[b->wds-1] == 0");
830 #endif
831 if (i -= j)
832 return i;
833 xa0 = a->x;
834 xa = xa0 + j;
835 xb0 = b->x;
836 xb = xb0 + j;
837 for(;;) {
838 if (*--xa != *--xb)
839 return *xa < *xb ? -1 : 1;
840 if (xa <= xa0)
841 break;
842 }
843 return 0;
844 }
845
846 /* Take the difference of Bigints a and b, returning a new Bigint. Returns
847 NULL on failure. The signs of a and b are ignored, but the sign of the
848 result is set appropriately. */
849
850 static Bigint *
diff(Bigint * a,Bigint * b)851 diff(Bigint *a, Bigint *b)
852 {
853 Bigint *c;
854 int i, wa, wb;
855 ULong *xa, *xae, *xb, *xbe, *xc;
856 ULLong borrow, y;
857
858 i = cmp(a,b);
859 if (!i) {
860 c = Balloc(0);
861 if (c == NULL)
862 return NULL;
863 c->wds = 1;
864 c->x[0] = 0;
865 return c;
866 }
867 if (i < 0) {
868 c = a;
869 a = b;
870 b = c;
871 i = 1;
872 }
873 else
874 i = 0;
875 c = Balloc(a->k);
876 if (c == NULL)
877 return NULL;
878 c->sign = i;
879 wa = a->wds;
880 xa = a->x;
881 xae = xa + wa;
882 wb = b->wds;
883 xb = b->x;
884 xbe = xb + wb;
885 xc = c->x;
886 borrow = 0;
887 do {
888 y = (ULLong)*xa++ - *xb++ - borrow;
889 borrow = y >> 32 & (ULong)1;
890 *xc++ = (ULong)(y & FFFFFFFF);
891 }
892 while(xb < xbe);
893 while(xa < xae) {
894 y = *xa++ - borrow;
895 borrow = y >> 32 & (ULong)1;
896 *xc++ = (ULong)(y & FFFFFFFF);
897 }
898 while(!*--xc)
899 wa--;
900 c->wds = wa;
901 return c;
902 }
903
904 /* Given a positive normal double x, return the difference between x and the
905 next double up. Doesn't give correct results for subnormals. */
906
907 static double
ulp(U * x)908 ulp(U *x)
909 {
910 Long L;
911 U u;
912
913 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
914 word0(&u) = L;
915 word1(&u) = 0;
916 return dval(&u);
917 }
918
919 /* Convert a Bigint to a double plus an exponent */
920
921 static double
b2d(Bigint * a,int * e)922 b2d(Bigint *a, int *e)
923 {
924 ULong *xa, *xa0, w, y, z;
925 int k;
926 U d;
927
928 xa0 = a->x;
929 xa = xa0 + a->wds;
930 y = *--xa;
931 #ifdef DEBUG
932 if (!y) Bug("zero y in b2d");
933 #endif
934 k = hi0bits(y);
935 *e = 32 - k;
936 if (k < Ebits) {
937 word0(&d) = Exp_1 | y >> (Ebits - k);
938 w = xa > xa0 ? *--xa : 0;
939 word1(&d) = y << ((32-Ebits) + k) | w >> (Ebits - k);
940 goto ret_d;
941 }
942 z = xa > xa0 ? *--xa : 0;
943 if (k -= Ebits) {
944 word0(&d) = Exp_1 | y << k | z >> (32 - k);
945 y = xa > xa0 ? *--xa : 0;
946 word1(&d) = z << k | y >> (32 - k);
947 }
948 else {
949 word0(&d) = Exp_1 | y;
950 word1(&d) = z;
951 }
952 ret_d:
953 return dval(&d);
954 }
955
956 /* Convert a scaled double to a Bigint plus an exponent. Similar to d2b,
957 except that it accepts the scale parameter used in _Py_dg_strtod (which
958 should be either 0 or 2*P), and the normalization for the return value is
959 different (see below). On input, d should be finite and nonnegative, and d
960 / 2**scale should be exactly representable as an IEEE 754 double.
961
962 Returns a Bigint b and an integer e such that
963
964 dval(d) / 2**scale = b * 2**e.
965
966 Unlike d2b, b is not necessarily odd: b and e are normalized so
967 that either 2**(P-1) <= b < 2**P and e >= Etiny, or b < 2**P
968 and e == Etiny. This applies equally to an input of 0.0: in that
969 case the return values are b = 0 and e = Etiny.
970
971 The above normalization ensures that for all possible inputs d,
972 2**e gives ulp(d/2**scale).
973
974 Returns NULL on failure.
975 */
976
977 static Bigint *
sd2b(U * d,int scale,int * e)978 sd2b(U *d, int scale, int *e)
979 {
980 Bigint *b;
981
982 b = Balloc(1);
983 if (b == NULL)
984 return NULL;
985
986 /* First construct b and e assuming that scale == 0. */
987 b->wds = 2;
988 b->x[0] = word1(d);
989 b->x[1] = word0(d) & Frac_mask;
990 *e = Etiny - 1 + (int)((word0(d) & Exp_mask) >> Exp_shift);
991 if (*e < Etiny)
992 *e = Etiny;
993 else
994 b->x[1] |= Exp_msk1;
995
996 /* Now adjust for scale, provided that b != 0. */
997 if (scale && (b->x[0] || b->x[1])) {
998 *e -= scale;
999 if (*e < Etiny) {
1000 scale = Etiny - *e;
1001 *e = Etiny;
1002 /* We can't shift more than P-1 bits without shifting out a 1. */
1003 assert(0 < scale && scale <= P - 1);
1004 if (scale >= 32) {
1005 /* The bits shifted out should all be zero. */
1006 assert(b->x[0] == 0);
1007 b->x[0] = b->x[1];
1008 b->x[1] = 0;
1009 scale -= 32;
1010 }
1011 if (scale) {
1012 /* The bits shifted out should all be zero. */
1013 assert(b->x[0] << (32 - scale) == 0);
1014 b->x[0] = (b->x[0] >> scale) | (b->x[1] << (32 - scale));
1015 b->x[1] >>= scale;
1016 }
1017 }
1018 }
1019 /* Ensure b is normalized. */
1020 if (!b->x[1])
1021 b->wds = 1;
1022
1023 return b;
1024 }
1025
1026 /* Convert a double to a Bigint plus an exponent. Return NULL on failure.
1027
1028 Given a finite nonzero double d, return an odd Bigint b and exponent *e
1029 such that fabs(d) = b * 2**e. On return, *bbits gives the number of
1030 significant bits of b; that is, 2**(*bbits-1) <= b < 2**(*bbits).
1031
1032 If d is zero, then b == 0, *e == -1010, *bbits = 0.
1033 */
1034
1035 static Bigint *
d2b(U * d,int * e,int * bits)1036 d2b(U *d, int *e, int *bits)
1037 {
1038 Bigint *b;
1039 int de, k;
1040 ULong *x, y, z;
1041 int i;
1042
1043 b = Balloc(1);
1044 if (b == NULL)
1045 return NULL;
1046 x = b->x;
1047
1048 z = word0(d) & Frac_mask;
1049 word0(d) &= 0x7fffffff; /* clear sign bit, which we ignore */
1050 if ((de = (int)(word0(d) >> Exp_shift)))
1051 z |= Exp_msk1;
1052 if ((y = word1(d))) {
1053 if ((k = lo0bits(&y))) {
1054 x[0] = y | z << (32 - k);
1055 z >>= k;
1056 }
1057 else
1058 x[0] = y;
1059 i =
1060 b->wds = (x[1] = z) ? 2 : 1;
1061 }
1062 else {
1063 k = lo0bits(&z);
1064 x[0] = z;
1065 i =
1066 b->wds = 1;
1067 k += 32;
1068 }
1069 if (de) {
1070 *e = de - Bias - (P-1) + k;
1071 *bits = P - k;
1072 }
1073 else {
1074 *e = de - Bias - (P-1) + 1 + k;
1075 *bits = 32*i - hi0bits(x[i-1]);
1076 }
1077 return b;
1078 }
1079
1080 /* Compute the ratio of two Bigints, as a double. The result may have an
1081 error of up to 2.5 ulps. */
1082
1083 static double
ratio(Bigint * a,Bigint * b)1084 ratio(Bigint *a, Bigint *b)
1085 {
1086 U da, db;
1087 int k, ka, kb;
1088
1089 dval(&da) = b2d(a, &ka);
1090 dval(&db) = b2d(b, &kb);
1091 k = ka - kb + 32*(a->wds - b->wds);
1092 if (k > 0)
1093 word0(&da) += k*Exp_msk1;
1094 else {
1095 k = -k;
1096 word0(&db) += k*Exp_msk1;
1097 }
1098 return dval(&da) / dval(&db);
1099 }
1100
1101 static const double
1102 tens[] = {
1103 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1104 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1105 1e20, 1e21, 1e22
1106 };
1107
1108 static const double
1109 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1110 static const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1111 9007199254740992.*9007199254740992.e-256
1112 /* = 2^106 * 1e-256 */
1113 };
1114 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1115 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1116 #define Scale_Bit 0x10
1117 #define n_bigtens 5
1118
1119 #define ULbits 32
1120 #define kshift 5
1121 #define kmask 31
1122
1123
1124 static int
dshift(Bigint * b,int p2)1125 dshift(Bigint *b, int p2)
1126 {
1127 int rv = hi0bits(b->x[b->wds-1]) - 4;
1128 if (p2 > 0)
1129 rv -= p2;
1130 return rv & kmask;
1131 }
1132
1133 /* special case of Bigint division. The quotient is always in the range 0 <=
1134 quotient < 10, and on entry the divisor S is normalized so that its top 4
1135 bits (28--31) are zero and bit 27 is set. */
1136
1137 static int
quorem(Bigint * b,Bigint * S)1138 quorem(Bigint *b, Bigint *S)
1139 {
1140 int n;
1141 ULong *bx, *bxe, q, *sx, *sxe;
1142 ULLong borrow, carry, y, ys;
1143
1144 n = S->wds;
1145 #ifdef DEBUG
1146 /*debug*/ if (b->wds > n)
1147 /*debug*/ Bug("oversize b in quorem");
1148 #endif
1149 if (b->wds < n)
1150 return 0;
1151 sx = S->x;
1152 sxe = sx + --n;
1153 bx = b->x;
1154 bxe = bx + n;
1155 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1156 #ifdef DEBUG
1157 /*debug*/ if (q > 9)
1158 /*debug*/ Bug("oversized quotient in quorem");
1159 #endif
1160 if (q) {
1161 borrow = 0;
1162 carry = 0;
1163 do {
1164 ys = *sx++ * (ULLong)q + carry;
1165 carry = ys >> 32;
1166 y = *bx - (ys & FFFFFFFF) - borrow;
1167 borrow = y >> 32 & (ULong)1;
1168 *bx++ = (ULong)(y & FFFFFFFF);
1169 }
1170 while(sx <= sxe);
1171 if (!*bxe) {
1172 bx = b->x;
1173 while(--bxe > bx && !*bxe)
1174 --n;
1175 b->wds = n;
1176 }
1177 }
1178 if (cmp(b, S) >= 0) {
1179 q++;
1180 borrow = 0;
1181 carry = 0;
1182 bx = b->x;
1183 sx = S->x;
1184 do {
1185 ys = *sx++ + carry;
1186 carry = ys >> 32;
1187 y = *bx - (ys & FFFFFFFF) - borrow;
1188 borrow = y >> 32 & (ULong)1;
1189 *bx++ = (ULong)(y & FFFFFFFF);
1190 }
1191 while(sx <= sxe);
1192 bx = b->x;
1193 bxe = bx + n;
1194 if (!*bxe) {
1195 while(--bxe > bx && !*bxe)
1196 --n;
1197 b->wds = n;
1198 }
1199 }
1200 return q;
1201 }
1202
1203 /* sulp(x) is a version of ulp(x) that takes bc.scale into account.
1204
1205 Assuming that x is finite and nonnegative (positive zero is fine
1206 here) and x / 2^bc.scale is exactly representable as a double,
1207 sulp(x) is equivalent to 2^bc.scale * ulp(x / 2^bc.scale). */
1208
1209 static double
sulp(U * x,BCinfo * bc)1210 sulp(U *x, BCinfo *bc)
1211 {
1212 U u;
1213
1214 if (bc->scale && 2*P + 1 > (int)((word0(x) & Exp_mask) >> Exp_shift)) {
1215 /* rv/2^bc->scale is subnormal */
1216 word0(&u) = (P+2)*Exp_msk1;
1217 word1(&u) = 0;
1218 return u.d;
1219 }
1220 else {
1221 assert(word0(x) || word1(x)); /* x != 0.0 */
1222 return ulp(x);
1223 }
1224 }
1225
1226 /* The bigcomp function handles some hard cases for strtod, for inputs
1227 with more than STRTOD_DIGLIM digits. It's called once an initial
1228 estimate for the double corresponding to the input string has
1229 already been obtained by the code in _Py_dg_strtod.
1230
1231 The bigcomp function is only called after _Py_dg_strtod has found a
1232 double value rv such that either rv or rv + 1ulp represents the
1233 correctly rounded value corresponding to the original string. It
1234 determines which of these two values is the correct one by
1235 computing the decimal digits of rv + 0.5ulp and comparing them with
1236 the corresponding digits of s0.
1237
1238 In the following, write dv for the absolute value of the number represented
1239 by the input string.
1240
1241 Inputs:
1242
1243 s0 points to the first significant digit of the input string.
1244
1245 rv is a (possibly scaled) estimate for the closest double value to the
1246 value represented by the original input to _Py_dg_strtod. If
1247 bc->scale is nonzero, then rv/2^(bc->scale) is the approximation to
1248 the input value.
1249
1250 bc is a struct containing information gathered during the parsing and
1251 estimation steps of _Py_dg_strtod. Description of fields follows:
1252
1253 bc->e0 gives the exponent of the input value, such that dv = (integer
1254 given by the bd->nd digits of s0) * 10**e0
1255
1256 bc->nd gives the total number of significant digits of s0. It will
1257 be at least 1.
1258
1259 bc->nd0 gives the number of significant digits of s0 before the
1260 decimal separator. If there's no decimal separator, bc->nd0 ==
1261 bc->nd.
1262
1263 bc->scale is the value used to scale rv to avoid doing arithmetic with
1264 subnormal values. It's either 0 or 2*P (=106).
1265
1266 Outputs:
1267
1268 On successful exit, rv/2^(bc->scale) is the closest double to dv.
1269
1270 Returns 0 on success, -1 on failure (e.g., due to a failed malloc call). */
1271
1272 static int
bigcomp(U * rv,const char * s0,BCinfo * bc)1273 bigcomp(U *rv, const char *s0, BCinfo *bc)
1274 {
1275 Bigint *b, *d;
1276 int b2, d2, dd, i, nd, nd0, odd, p2, p5;
1277
1278 nd = bc->nd;
1279 nd0 = bc->nd0;
1280 p5 = nd + bc->e0;
1281 b = sd2b(rv, bc->scale, &p2);
1282 if (b == NULL)
1283 return -1;
1284
1285 /* record whether the lsb of rv/2^(bc->scale) is odd: in the exact halfway
1286 case, this is used for round to even. */
1287 odd = b->x[0] & 1;
1288
1289 /* left shift b by 1 bit and or a 1 into the least significant bit;
1290 this gives us b * 2**p2 = rv/2^(bc->scale) + 0.5 ulp. */
1291 b = lshift(b, 1);
1292 if (b == NULL)
1293 return -1;
1294 b->x[0] |= 1;
1295 p2--;
1296
1297 p2 -= p5;
1298 d = i2b(1);
1299 if (d == NULL) {
1300 Bfree(b);
1301 return -1;
1302 }
1303 /* Arrange for convenient computation of quotients:
1304 * shift left if necessary so divisor has 4 leading 0 bits.
1305 */
1306 if (p5 > 0) {
1307 d = pow5mult(d, p5);
1308 if (d == NULL) {
1309 Bfree(b);
1310 return -1;
1311 }
1312 }
1313 else if (p5 < 0) {
1314 b = pow5mult(b, -p5);
1315 if (b == NULL) {
1316 Bfree(d);
1317 return -1;
1318 }
1319 }
1320 if (p2 > 0) {
1321 b2 = p2;
1322 d2 = 0;
1323 }
1324 else {
1325 b2 = 0;
1326 d2 = -p2;
1327 }
1328 i = dshift(d, d2);
1329 if ((b2 += i) > 0) {
1330 b = lshift(b, b2);
1331 if (b == NULL) {
1332 Bfree(d);
1333 return -1;
1334 }
1335 }
1336 if ((d2 += i) > 0) {
1337 d = lshift(d, d2);
1338 if (d == NULL) {
1339 Bfree(b);
1340 return -1;
1341 }
1342 }
1343
1344 /* Compare s0 with b/d: set dd to -1, 0, or 1 according as s0 < b/d, s0 ==
1345 * b/d, or s0 > b/d. Here the digits of s0 are thought of as representing
1346 * a number in the range [0.1, 1). */
1347 if (cmp(b, d) >= 0)
1348 /* b/d >= 1 */
1349 dd = -1;
1350 else {
1351 i = 0;
1352 for(;;) {
1353 b = multadd(b, 10, 0);
1354 if (b == NULL) {
1355 Bfree(d);
1356 return -1;
1357 }
1358 dd = s0[i < nd0 ? i : i+1] - '0' - quorem(b, d);
1359 i++;
1360
1361 if (dd)
1362 break;
1363 if (!b->x[0] && b->wds == 1) {
1364 /* b/d == 0 */
1365 dd = i < nd;
1366 break;
1367 }
1368 if (!(i < nd)) {
1369 /* b/d != 0, but digits of s0 exhausted */
1370 dd = -1;
1371 break;
1372 }
1373 }
1374 }
1375 Bfree(b);
1376 Bfree(d);
1377 if (dd > 0 || (dd == 0 && odd))
1378 dval(rv) += sulp(rv, bc);
1379 return 0;
1380 }
1381
1382
1383 double
_Py_dg_strtod(const char * s00,char ** se)1384 _Py_dg_strtod(const char *s00, char **se)
1385 {
1386 int bb2, bb5, bbe, bd2, bd5, bs2, c, dsign, e, e1, error;
1387 int esign, i, j, k, lz, nd, nd0, odd, sign;
1388 const char *s, *s0, *s1;
1389 double aadj, aadj1;
1390 U aadj2, adj, rv, rv0;
1391 ULong y, z, abs_exp;
1392 Long L;
1393 BCinfo bc;
1394 Bigint *bb = NULL, *bd = NULL, *bd0 = NULL, *bs = NULL, *delta = NULL;
1395 size_t ndigits, fraclen;
1396 double result;
1397
1398 dval(&rv) = 0.;
1399
1400 /* Start parsing. */
1401 c = *(s = s00);
1402
1403 /* Parse optional sign, if present. */
1404 sign = 0;
1405 switch (c) {
1406 case '-':
1407 sign = 1;
1408 /* fall through */
1409 case '+':
1410 c = *++s;
1411 }
1412
1413 /* Skip leading zeros: lz is true iff there were leading zeros. */
1414 s1 = s;
1415 while (c == '0')
1416 c = *++s;
1417 lz = s != s1;
1418
1419 /* Point s0 at the first nonzero digit (if any). fraclen will be the
1420 number of digits between the decimal point and the end of the
1421 digit string. ndigits will be the total number of digits ignoring
1422 leading zeros. */
1423 s0 = s1 = s;
1424 while ('0' <= c && c <= '9')
1425 c = *++s;
1426 ndigits = s - s1;
1427 fraclen = 0;
1428
1429 /* Parse decimal point and following digits. */
1430 if (c == '.') {
1431 c = *++s;
1432 if (!ndigits) {
1433 s1 = s;
1434 while (c == '0')
1435 c = *++s;
1436 lz = lz || s != s1;
1437 fraclen += (s - s1);
1438 s0 = s;
1439 }
1440 s1 = s;
1441 while ('0' <= c && c <= '9')
1442 c = *++s;
1443 ndigits += s - s1;
1444 fraclen += s - s1;
1445 }
1446
1447 /* Now lz is true if and only if there were leading zero digits, and
1448 ndigits gives the total number of digits ignoring leading zeros. A
1449 valid input must have at least one digit. */
1450 if (!ndigits && !lz) {
1451 if (se)
1452 *se = (char *)s00;
1453 goto parse_error;
1454 }
1455
1456 /* Range check ndigits and fraclen to make sure that they, and values
1457 computed with them, can safely fit in an int. */
1458 if (ndigits > MAX_DIGITS || fraclen > MAX_DIGITS) {
1459 if (se)
1460 *se = (char *)s00;
1461 goto parse_error;
1462 }
1463 nd = (int)ndigits;
1464 nd0 = (int)ndigits - (int)fraclen;
1465
1466 /* Parse exponent. */
1467 e = 0;
1468 if (c == 'e' || c == 'E') {
1469 s00 = s;
1470 c = *++s;
1471
1472 /* Exponent sign. */
1473 esign = 0;
1474 switch (c) {
1475 case '-':
1476 esign = 1;
1477 /* fall through */
1478 case '+':
1479 c = *++s;
1480 }
1481
1482 /* Skip zeros. lz is true iff there are leading zeros. */
1483 s1 = s;
1484 while (c == '0')
1485 c = *++s;
1486 lz = s != s1;
1487
1488 /* Get absolute value of the exponent. */
1489 s1 = s;
1490 abs_exp = 0;
1491 while ('0' <= c && c <= '9') {
1492 abs_exp = 10*abs_exp + (c - '0');
1493 c = *++s;
1494 }
1495
1496 /* abs_exp will be correct modulo 2**32. But 10**9 < 2**32, so if
1497 there are at most 9 significant exponent digits then overflow is
1498 impossible. */
1499 if (s - s1 > 9 || abs_exp > MAX_ABS_EXP)
1500 e = (int)MAX_ABS_EXP;
1501 else
1502 e = (int)abs_exp;
1503 if (esign)
1504 e = -e;
1505
1506 /* A valid exponent must have at least one digit. */
1507 if (s == s1 && !lz)
1508 s = s00;
1509 }
1510
1511 /* Adjust exponent to take into account position of the point. */
1512 e -= nd - nd0;
1513 if (nd0 <= 0)
1514 nd0 = nd;
1515
1516 /* Finished parsing. Set se to indicate how far we parsed */
1517 if (se)
1518 *se = (char *)s;
1519
1520 /* If all digits were zero, exit with return value +-0.0. Otherwise,
1521 strip trailing zeros: scan back until we hit a nonzero digit. */
1522 if (!nd)
1523 goto ret;
1524 for (i = nd; i > 0; ) {
1525 --i;
1526 if (s0[i < nd0 ? i : i+1] != '0') {
1527 ++i;
1528 break;
1529 }
1530 }
1531 e += nd - i;
1532 nd = i;
1533 if (nd0 > nd)
1534 nd0 = nd;
1535
1536 /* Summary of parsing results. After parsing, and dealing with zero
1537 * inputs, we have values s0, nd0, nd, e, sign, where:
1538 *
1539 * - s0 points to the first significant digit of the input string
1540 *
1541 * - nd is the total number of significant digits (here, and
1542 * below, 'significant digits' means the set of digits of the
1543 * significand of the input that remain after ignoring leading
1544 * and trailing zeros).
1545 *
1546 * - nd0 indicates the position of the decimal point, if present; it
1547 * satisfies 1 <= nd0 <= nd. The nd significant digits are in
1548 * s0[0:nd0] and s0[nd0+1:nd+1] using the usual Python half-open slice
1549 * notation. (If nd0 < nd, then s0[nd0] contains a '.' character; if
1550 * nd0 == nd, then s0[nd0] could be any non-digit character.)
1551 *
1552 * - e is the adjusted exponent: the absolute value of the number
1553 * represented by the original input string is n * 10**e, where
1554 * n is the integer represented by the concatenation of
1555 * s0[0:nd0] and s0[nd0+1:nd+1]
1556 *
1557 * - sign gives the sign of the input: 1 for negative, 0 for positive
1558 *
1559 * - the first and last significant digits are nonzero
1560 */
1561
1562 /* put first DBL_DIG+1 digits into integer y and z.
1563 *
1564 * - y contains the value represented by the first min(9, nd)
1565 * significant digits
1566 *
1567 * - if nd > 9, z contains the value represented by significant digits
1568 * with indices in [9, min(16, nd)). So y * 10**(min(16, nd) - 9) + z
1569 * gives the value represented by the first min(16, nd) sig. digits.
1570 */
1571
1572 bc.e0 = e1 = e;
1573 y = z = 0;
1574 for (i = 0; i < nd; i++) {
1575 if (i < 9)
1576 y = 10*y + s0[i < nd0 ? i : i+1] - '0';
1577 else if (i < DBL_DIG+1)
1578 z = 10*z + s0[i < nd0 ? i : i+1] - '0';
1579 else
1580 break;
1581 }
1582
1583 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1584 dval(&rv) = y;
1585 if (k > 9) {
1586 dval(&rv) = tens[k - 9] * dval(&rv) + z;
1587 }
1588 if (nd <= DBL_DIG
1589 && Flt_Rounds == 1
1590 ) {
1591 if (!e)
1592 goto ret;
1593 if (e > 0) {
1594 if (e <= Ten_pmax) {
1595 dval(&rv) *= tens[e];
1596 goto ret;
1597 }
1598 i = DBL_DIG - nd;
1599 if (e <= Ten_pmax + i) {
1600 /* A fancier test would sometimes let us do
1601 * this for larger i values.
1602 */
1603 e -= i;
1604 dval(&rv) *= tens[i];
1605 dval(&rv) *= tens[e];
1606 goto ret;
1607 }
1608 }
1609 else if (e >= -Ten_pmax) {
1610 dval(&rv) /= tens[-e];
1611 goto ret;
1612 }
1613 }
1614 e1 += nd - k;
1615
1616 bc.scale = 0;
1617
1618 /* Get starting approximation = rv * 10**e1 */
1619
1620 if (e1 > 0) {
1621 if ((i = e1 & 15))
1622 dval(&rv) *= tens[i];
1623 if (e1 &= ~15) {
1624 if (e1 > DBL_MAX_10_EXP)
1625 goto ovfl;
1626 e1 >>= 4;
1627 for(j = 0; e1 > 1; j++, e1 >>= 1)
1628 if (e1 & 1)
1629 dval(&rv) *= bigtens[j];
1630 /* The last multiplication could overflow. */
1631 word0(&rv) -= P*Exp_msk1;
1632 dval(&rv) *= bigtens[j];
1633 if ((z = word0(&rv) & Exp_mask)
1634 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1635 goto ovfl;
1636 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1637 /* set to largest number */
1638 /* (Can't trust DBL_MAX) */
1639 word0(&rv) = Big0;
1640 word1(&rv) = Big1;
1641 }
1642 else
1643 word0(&rv) += P*Exp_msk1;
1644 }
1645 }
1646 else if (e1 < 0) {
1647 /* The input decimal value lies in [10**e1, 10**(e1+16)).
1648
1649 If e1 <= -512, underflow immediately.
1650 If e1 <= -256, set bc.scale to 2*P.
1651
1652 So for input value < 1e-256, bc.scale is always set;
1653 for input value >= 1e-240, bc.scale is never set.
1654 For input values in [1e-256, 1e-240), bc.scale may or may
1655 not be set. */
1656
1657 e1 = -e1;
1658 if ((i = e1 & 15))
1659 dval(&rv) /= tens[i];
1660 if (e1 >>= 4) {
1661 if (e1 >= 1 << n_bigtens)
1662 goto undfl;
1663 if (e1 & Scale_Bit)
1664 bc.scale = 2*P;
1665 for(j = 0; e1 > 0; j++, e1 >>= 1)
1666 if (e1 & 1)
1667 dval(&rv) *= tinytens[j];
1668 if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
1669 >> Exp_shift)) > 0) {
1670 /* scaled rv is denormal; clear j low bits */
1671 if (j >= 32) {
1672 word1(&rv) = 0;
1673 if (j >= 53)
1674 word0(&rv) = (P+2)*Exp_msk1;
1675 else
1676 word0(&rv) &= 0xffffffff << (j-32);
1677 }
1678 else
1679 word1(&rv) &= 0xffffffff << j;
1680 }
1681 if (!dval(&rv))
1682 goto undfl;
1683 }
1684 }
1685
1686 /* Now the hard part -- adjusting rv to the correct value.*/
1687
1688 /* Put digits into bd: true value = bd * 10^e */
1689
1690 bc.nd = nd;
1691 bc.nd0 = nd0; /* Only needed if nd > STRTOD_DIGLIM, but done here */
1692 /* to silence an erroneous warning about bc.nd0 */
1693 /* possibly not being initialized. */
1694 if (nd > STRTOD_DIGLIM) {
1695 /* ASSERT(STRTOD_DIGLIM >= 18); 18 == one more than the */
1696 /* minimum number of decimal digits to distinguish double values */
1697 /* in IEEE arithmetic. */
1698
1699 /* Truncate input to 18 significant digits, then discard any trailing
1700 zeros on the result by updating nd, nd0, e and y suitably. (There's
1701 no need to update z; it's not reused beyond this point.) */
1702 for (i = 18; i > 0; ) {
1703 /* scan back until we hit a nonzero digit. significant digit 'i'
1704 is s0[i] if i < nd0, s0[i+1] if i >= nd0. */
1705 --i;
1706 if (s0[i < nd0 ? i : i+1] != '0') {
1707 ++i;
1708 break;
1709 }
1710 }
1711 e += nd - i;
1712 nd = i;
1713 if (nd0 > nd)
1714 nd0 = nd;
1715 if (nd < 9) { /* must recompute y */
1716 y = 0;
1717 for(i = 0; i < nd0; ++i)
1718 y = 10*y + s0[i] - '0';
1719 for(; i < nd; ++i)
1720 y = 10*y + s0[i+1] - '0';
1721 }
1722 }
1723 bd0 = s2b(s0, nd0, nd, y);
1724 if (bd0 == NULL)
1725 goto failed_malloc;
1726
1727 /* Notation for the comments below. Write:
1728
1729 - dv for the absolute value of the number represented by the original
1730 decimal input string.
1731
1732 - if we've truncated dv, write tdv for the truncated value.
1733 Otherwise, set tdv == dv.
1734
1735 - srv for the quantity rv/2^bc.scale; so srv is the current binary
1736 approximation to tdv (and dv). It should be exactly representable
1737 in an IEEE 754 double.
1738 */
1739
1740 for(;;) {
1741
1742 /* This is the main correction loop for _Py_dg_strtod.
1743
1744 We've got a decimal value tdv, and a floating-point approximation
1745 srv=rv/2^bc.scale to tdv. The aim is to determine whether srv is
1746 close enough (i.e., within 0.5 ulps) to tdv, and to compute a new
1747 approximation if not.
1748
1749 To determine whether srv is close enough to tdv, compute integers
1750 bd, bb and bs proportional to tdv, srv and 0.5 ulp(srv)
1751 respectively, and then use integer arithmetic to determine whether
1752 |tdv - srv| is less than, equal to, or greater than 0.5 ulp(srv).
1753 */
1754
1755 bd = Balloc(bd0->k);
1756 if (bd == NULL) {
1757 goto failed_malloc;
1758 }
1759 Bcopy(bd, bd0);
1760 bb = sd2b(&rv, bc.scale, &bbe); /* srv = bb * 2^bbe */
1761 if (bb == NULL) {
1762 goto failed_malloc;
1763 }
1764 /* Record whether lsb of bb is odd, in case we need this
1765 for the round-to-even step later. */
1766 odd = bb->x[0] & 1;
1767
1768 /* tdv = bd * 10**e; srv = bb * 2**bbe */
1769 bs = i2b(1);
1770 if (bs == NULL) {
1771 goto failed_malloc;
1772 }
1773
1774 if (e >= 0) {
1775 bb2 = bb5 = 0;
1776 bd2 = bd5 = e;
1777 }
1778 else {
1779 bb2 = bb5 = -e;
1780 bd2 = bd5 = 0;
1781 }
1782 if (bbe >= 0)
1783 bb2 += bbe;
1784 else
1785 bd2 -= bbe;
1786 bs2 = bb2;
1787 bb2++;
1788 bd2++;
1789
1790 /* At this stage bd5 - bb5 == e == bd2 - bb2 + bbe, bb2 - bs2 == 1,
1791 and bs == 1, so:
1792
1793 tdv == bd * 10**e = bd * 2**(bbe - bb2 + bd2) * 5**(bd5 - bb5)
1794 srv == bb * 2**bbe = bb * 2**(bbe - bb2 + bb2)
1795 0.5 ulp(srv) == 2**(bbe-1) = bs * 2**(bbe - bb2 + bs2)
1796
1797 It follows that:
1798
1799 M * tdv = bd * 2**bd2 * 5**bd5
1800 M * srv = bb * 2**bb2 * 5**bb5
1801 M * 0.5 ulp(srv) = bs * 2**bs2 * 5**bb5
1802
1803 for some constant M. (Actually, M == 2**(bb2 - bbe) * 5**bb5, but
1804 this fact is not needed below.)
1805 */
1806
1807 /* Remove factor of 2**i, where i = min(bb2, bd2, bs2). */
1808 i = bb2 < bd2 ? bb2 : bd2;
1809 if (i > bs2)
1810 i = bs2;
1811 if (i > 0) {
1812 bb2 -= i;
1813 bd2 -= i;
1814 bs2 -= i;
1815 }
1816
1817 /* Scale bb, bd, bs by the appropriate powers of 2 and 5. */
1818 if (bb5 > 0) {
1819 bs = pow5mult(bs, bb5);
1820 if (bs == NULL) {
1821 goto failed_malloc;
1822 }
1823 Bigint *bb1 = mult(bs, bb);
1824 Bfree(bb);
1825 bb = bb1;
1826 if (bb == NULL) {
1827 goto failed_malloc;
1828 }
1829 }
1830 if (bb2 > 0) {
1831 bb = lshift(bb, bb2);
1832 if (bb == NULL) {
1833 goto failed_malloc;
1834 }
1835 }
1836 if (bd5 > 0) {
1837 bd = pow5mult(bd, bd5);
1838 if (bd == NULL) {
1839 goto failed_malloc;
1840 }
1841 }
1842 if (bd2 > 0) {
1843 bd = lshift(bd, bd2);
1844 if (bd == NULL) {
1845 goto failed_malloc;
1846 }
1847 }
1848 if (bs2 > 0) {
1849 bs = lshift(bs, bs2);
1850 if (bs == NULL) {
1851 goto failed_malloc;
1852 }
1853 }
1854
1855 /* Now bd, bb and bs are scaled versions of tdv, srv and 0.5 ulp(srv),
1856 respectively. Compute the difference |tdv - srv|, and compare
1857 with 0.5 ulp(srv). */
1858
1859 delta = diff(bb, bd);
1860 if (delta == NULL) {
1861 goto failed_malloc;
1862 }
1863 dsign = delta->sign;
1864 delta->sign = 0;
1865 i = cmp(delta, bs);
1866 if (bc.nd > nd && i <= 0) {
1867 if (dsign)
1868 break; /* Must use bigcomp(). */
1869
1870 /* Here rv overestimates the truncated decimal value by at most
1871 0.5 ulp(rv). Hence rv either overestimates the true decimal
1872 value by <= 0.5 ulp(rv), or underestimates it by some small
1873 amount (< 0.1 ulp(rv)); either way, rv is within 0.5 ulps of
1874 the true decimal value, so it's possible to exit.
1875
1876 Exception: if scaled rv is a normal exact power of 2, but not
1877 DBL_MIN, then rv - 0.5 ulp(rv) takes us all the way down to the
1878 next double, so the correctly rounded result is either rv - 0.5
1879 ulp(rv) or rv; in this case, use bigcomp to distinguish. */
1880
1881 if (!word1(&rv) && !(word0(&rv) & Bndry_mask)) {
1882 /* rv can't be 0, since it's an overestimate for some
1883 nonzero value. So rv is a normal power of 2. */
1884 j = (int)(word0(&rv) & Exp_mask) >> Exp_shift;
1885 /* rv / 2^bc.scale = 2^(j - 1023 - bc.scale); use bigcomp if
1886 rv / 2^bc.scale >= 2^-1021. */
1887 if (j - bc.scale >= 2) {
1888 dval(&rv) -= 0.5 * sulp(&rv, &bc);
1889 break; /* Use bigcomp. */
1890 }
1891 }
1892
1893 {
1894 bc.nd = nd;
1895 i = -1; /* Discarded digits make delta smaller. */
1896 }
1897 }
1898
1899 if (i < 0) {
1900 /* Error is less than half an ulp -- check for
1901 * special case of mantissa a power of two.
1902 */
1903 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask
1904 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
1905 ) {
1906 break;
1907 }
1908 if (!delta->x[0] && delta->wds <= 1) {
1909 /* exact result */
1910 break;
1911 }
1912 delta = lshift(delta,Log2P);
1913 if (delta == NULL) {
1914 goto failed_malloc;
1915 }
1916 if (cmp(delta, bs) > 0)
1917 goto drop_down;
1918 break;
1919 }
1920 if (i == 0) {
1921 /* exactly half-way between */
1922 if (dsign) {
1923 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
1924 && word1(&rv) == (
1925 (bc.scale &&
1926 (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) ?
1927 (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
1928 0xffffffff)) {
1929 /*boundary case -- increment exponent*/
1930 word0(&rv) = (word0(&rv) & Exp_mask)
1931 + Exp_msk1
1932 ;
1933 word1(&rv) = 0;
1934 /* dsign = 0; */
1935 break;
1936 }
1937 }
1938 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
1939 drop_down:
1940 /* boundary case -- decrement exponent */
1941 if (bc.scale) {
1942 L = word0(&rv) & Exp_mask;
1943 if (L <= (2*P+1)*Exp_msk1) {
1944 if (L > (P+2)*Exp_msk1)
1945 /* round even ==> */
1946 /* accept rv */
1947 break;
1948 /* rv = smallest denormal */
1949 if (bc.nd > nd)
1950 break;
1951 goto undfl;
1952 }
1953 }
1954 L = (word0(&rv) & Exp_mask) - Exp_msk1;
1955 word0(&rv) = L | Bndry_mask1;
1956 word1(&rv) = 0xffffffff;
1957 break;
1958 }
1959 if (!odd)
1960 break;
1961 if (dsign)
1962 dval(&rv) += sulp(&rv, &bc);
1963 else {
1964 dval(&rv) -= sulp(&rv, &bc);
1965 if (!dval(&rv)) {
1966 if (bc.nd >nd)
1967 break;
1968 goto undfl;
1969 }
1970 }
1971 /* dsign = 1 - dsign; */
1972 break;
1973 }
1974 if ((aadj = ratio(delta, bs)) <= 2.) {
1975 if (dsign)
1976 aadj = aadj1 = 1.;
1977 else if (word1(&rv) || word0(&rv) & Bndry_mask) {
1978 if (word1(&rv) == Tiny1 && !word0(&rv)) {
1979 if (bc.nd >nd)
1980 break;
1981 goto undfl;
1982 }
1983 aadj = 1.;
1984 aadj1 = -1.;
1985 }
1986 else {
1987 /* special case -- power of FLT_RADIX to be */
1988 /* rounded down... */
1989
1990 if (aadj < 2./FLT_RADIX)
1991 aadj = 1./FLT_RADIX;
1992 else
1993 aadj *= 0.5;
1994 aadj1 = -aadj;
1995 }
1996 }
1997 else {
1998 aadj *= 0.5;
1999 aadj1 = dsign ? aadj : -aadj;
2000 if (Flt_Rounds == 0)
2001 aadj1 += 0.5;
2002 }
2003 y = word0(&rv) & Exp_mask;
2004
2005 /* Check for overflow */
2006
2007 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
2008 dval(&rv0) = dval(&rv);
2009 word0(&rv) -= P*Exp_msk1;
2010 adj.d = aadj1 * ulp(&rv);
2011 dval(&rv) += adj.d;
2012 if ((word0(&rv) & Exp_mask) >=
2013 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
2014 if (word0(&rv0) == Big0 && word1(&rv0) == Big1) {
2015 goto ovfl;
2016 }
2017 word0(&rv) = Big0;
2018 word1(&rv) = Big1;
2019 goto cont;
2020 }
2021 else
2022 word0(&rv) += P*Exp_msk1;
2023 }
2024 else {
2025 if (bc.scale && y <= 2*P*Exp_msk1) {
2026 if (aadj <= 0x7fffffff) {
2027 if ((z = (ULong)aadj) <= 0)
2028 z = 1;
2029 aadj = z;
2030 aadj1 = dsign ? aadj : -aadj;
2031 }
2032 dval(&aadj2) = aadj1;
2033 word0(&aadj2) += (2*P+1)*Exp_msk1 - y;
2034 aadj1 = dval(&aadj2);
2035 }
2036 adj.d = aadj1 * ulp(&rv);
2037 dval(&rv) += adj.d;
2038 }
2039 z = word0(&rv) & Exp_mask;
2040 if (bc.nd == nd) {
2041 if (!bc.scale)
2042 if (y == z) {
2043 /* Can we stop now? */
2044 L = (Long)aadj;
2045 aadj -= L;
2046 /* The tolerances below are conservative. */
2047 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
2048 if (aadj < .4999999 || aadj > .5000001)
2049 break;
2050 }
2051 else if (aadj < .4999999/FLT_RADIX)
2052 break;
2053 }
2054 }
2055 cont:
2056 Bfree(bb); bb = NULL;
2057 Bfree(bd); bd = NULL;
2058 Bfree(bs); bs = NULL;
2059 Bfree(delta); delta = NULL;
2060 }
2061 if (bc.nd > nd) {
2062 error = bigcomp(&rv, s0, &bc);
2063 if (error)
2064 goto failed_malloc;
2065 }
2066
2067 if (bc.scale) {
2068 word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
2069 word1(&rv0) = 0;
2070 dval(&rv) *= dval(&rv0);
2071 }
2072
2073 ret:
2074 result = sign ? -dval(&rv) : dval(&rv);
2075 goto done;
2076
2077 parse_error:
2078 result = 0.0;
2079 goto done;
2080
2081 failed_malloc:
2082 errno = ENOMEM;
2083 result = -1.0;
2084 goto done;
2085
2086 undfl:
2087 result = sign ? -0.0 : 0.0;
2088 goto done;
2089
2090 ovfl:
2091 errno = ERANGE;
2092 /* Can't trust HUGE_VAL */
2093 word0(&rv) = Exp_mask;
2094 word1(&rv) = 0;
2095 result = sign ? -dval(&rv) : dval(&rv);
2096 goto done;
2097
2098 done:
2099 Bfree(bb);
2100 Bfree(bd);
2101 Bfree(bs);
2102 Bfree(bd0);
2103 Bfree(delta);
2104 return result;
2105
2106 }
2107
2108 static char *
rv_alloc(int i)2109 rv_alloc(int i)
2110 {
2111 int j, k, *r;
2112
2113 j = sizeof(ULong);
2114 for(k = 0;
2115 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (unsigned)i;
2116 j <<= 1)
2117 k++;
2118 r = (int*)Balloc(k);
2119 if (r == NULL)
2120 return NULL;
2121 *r = k;
2122 return (char *)(r+1);
2123 }
2124
2125 static char *
nrv_alloc(const char * s,char ** rve,int n)2126 nrv_alloc(const char *s, char **rve, int n)
2127 {
2128 char *rv, *t;
2129
2130 rv = rv_alloc(n);
2131 if (rv == NULL)
2132 return NULL;
2133 t = rv;
2134 while((*t = *s++)) t++;
2135 if (rve)
2136 *rve = t;
2137 return rv;
2138 }
2139
2140 /* freedtoa(s) must be used to free values s returned by dtoa
2141 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
2142 * but for consistency with earlier versions of dtoa, it is optional
2143 * when MULTIPLE_THREADS is not defined.
2144 */
2145
2146 void
_Py_dg_freedtoa(char * s)2147 _Py_dg_freedtoa(char *s)
2148 {
2149 Bigint *b = (Bigint *)((int *)s - 1);
2150 b->maxwds = 1 << (b->k = *(int*)b);
2151 Bfree(b);
2152 }
2153
2154 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
2155 *
2156 * Inspired by "How to Print Floating-Point Numbers Accurately" by
2157 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
2158 *
2159 * Modifications:
2160 * 1. Rather than iterating, we use a simple numeric overestimate
2161 * to determine k = floor(log10(d)). We scale relevant
2162 * quantities using O(log2(k)) rather than O(k) multiplications.
2163 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
2164 * try to generate digits strictly left to right. Instead, we
2165 * compute with fewer bits and propagate the carry if necessary
2166 * when rounding the final digit up. This is often faster.
2167 * 3. Under the assumption that input will be rounded nearest,
2168 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
2169 * That is, we allow equality in stopping tests when the
2170 * round-nearest rule will give the same floating-point value
2171 * as would satisfaction of the stopping test with strict
2172 * inequality.
2173 * 4. We remove common factors of powers of 2 from relevant
2174 * quantities.
2175 * 5. When converting floating-point integers less than 1e16,
2176 * we use floating-point arithmetic rather than resorting
2177 * to multiple-precision integers.
2178 * 6. When asked to produce fewer than 15 digits, we first try
2179 * to get by with floating-point arithmetic; we resort to
2180 * multiple-precision integer arithmetic only if we cannot
2181 * guarantee that the floating-point calculation has given
2182 * the correctly rounded result. For k requested digits and
2183 * "uniformly" distributed input, the probability is
2184 * something like 10^(k-15) that we must resort to the Long
2185 * calculation.
2186 */
2187
2188 /* Additional notes (METD): (1) returns NULL on failure. (2) to avoid memory
2189 leakage, a successful call to _Py_dg_dtoa should always be matched by a
2190 call to _Py_dg_freedtoa. */
2191
2192 char *
_Py_dg_dtoa(double dd,int mode,int ndigits,int * decpt,int * sign,char ** rve)2193 _Py_dg_dtoa(double dd, int mode, int ndigits,
2194 int *decpt, int *sign, char **rve)
2195 {
2196 /* Arguments ndigits, decpt, sign are similar to those
2197 of ecvt and fcvt; trailing zeros are suppressed from
2198 the returned string. If not null, *rve is set to point
2199 to the end of the return value. If d is +-Infinity or NaN,
2200 then *decpt is set to 9999.
2201
2202 mode:
2203 0 ==> shortest string that yields d when read in
2204 and rounded to nearest.
2205 1 ==> like 0, but with Steele & White stopping rule;
2206 e.g. with IEEE P754 arithmetic , mode 0 gives
2207 1e23 whereas mode 1 gives 9.999999999999999e22.
2208 2 ==> max(1,ndigits) significant digits. This gives a
2209 return value similar to that of ecvt, except
2210 that trailing zeros are suppressed.
2211 3 ==> through ndigits past the decimal point. This
2212 gives a return value similar to that from fcvt,
2213 except that trailing zeros are suppressed, and
2214 ndigits can be negative.
2215 4,5 ==> similar to 2 and 3, respectively, but (in
2216 round-nearest mode) with the tests of mode 0 to
2217 possibly return a shorter string that rounds to d.
2218 With IEEE arithmetic and compilation with
2219 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
2220 as modes 2 and 3 when FLT_ROUNDS != 1.
2221 6-9 ==> Debugging modes similar to mode - 4: don't try
2222 fast floating-point estimate (if applicable).
2223
2224 Values of mode other than 0-9 are treated as mode 0.
2225
2226 Sufficient space is allocated to the return value
2227 to hold the suppressed trailing zeros.
2228 */
2229
2230 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
2231 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
2232 spec_case, try_quick;
2233 Long L;
2234 int denorm;
2235 ULong x;
2236 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
2237 U d2, eps, u;
2238 double ds;
2239 char *s, *s0;
2240
2241 /* set pointers to NULL, to silence gcc compiler warnings and make
2242 cleanup easier on error */
2243 mlo = mhi = S = 0;
2244 s0 = 0;
2245
2246 u.d = dd;
2247 if (word0(&u) & Sign_bit) {
2248 /* set sign for everything, including 0's and NaNs */
2249 *sign = 1;
2250 word0(&u) &= ~Sign_bit; /* clear sign bit */
2251 }
2252 else
2253 *sign = 0;
2254
2255 /* quick return for Infinities, NaNs and zeros */
2256 if ((word0(&u) & Exp_mask) == Exp_mask)
2257 {
2258 /* Infinity or NaN */
2259 *decpt = 9999;
2260 if (!word1(&u) && !(word0(&u) & 0xfffff))
2261 return nrv_alloc("Infinity", rve, 8);
2262 return nrv_alloc("NaN", rve, 3);
2263 }
2264 if (!dval(&u)) {
2265 *decpt = 1;
2266 return nrv_alloc("0", rve, 1);
2267 }
2268
2269 /* compute k = floor(log10(d)). The computation may leave k
2270 one too large, but should never leave k too small. */
2271 b = d2b(&u, &be, &bbits);
2272 if (b == NULL)
2273 goto failed_malloc;
2274 if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
2275 dval(&d2) = dval(&u);
2276 word0(&d2) &= Frac_mask1;
2277 word0(&d2) |= Exp_11;
2278
2279 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2280 * log10(x) = log(x) / log(10)
2281 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2282 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2283 *
2284 * This suggests computing an approximation k to log10(d) by
2285 *
2286 * k = (i - Bias)*0.301029995663981
2287 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2288 *
2289 * We want k to be too large rather than too small.
2290 * The error in the first-order Taylor series approximation
2291 * is in our favor, so we just round up the constant enough
2292 * to compensate for any error in the multiplication of
2293 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2294 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2295 * adding 1e-13 to the constant term more than suffices.
2296 * Hence we adjust the constant term to 0.1760912590558.
2297 * (We could get a more accurate k by invoking log10,
2298 * but this is probably not worthwhile.)
2299 */
2300
2301 i -= Bias;
2302 denorm = 0;
2303 }
2304 else {
2305 /* d is denormalized */
2306
2307 i = bbits + be + (Bias + (P-1) - 1);
2308 x = i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32)
2309 : word1(&u) << (32 - i);
2310 dval(&d2) = x;
2311 word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
2312 i -= (Bias + (P-1) - 1) + 1;
2313 denorm = 1;
2314 }
2315 ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 +
2316 i*0.301029995663981;
2317 k = (int)ds;
2318 if (ds < 0. && ds != k)
2319 k--; /* want k = floor(ds) */
2320 k_check = 1;
2321 if (k >= 0 && k <= Ten_pmax) {
2322 if (dval(&u) < tens[k])
2323 k--;
2324 k_check = 0;
2325 }
2326 j = bbits - i - 1;
2327 if (j >= 0) {
2328 b2 = 0;
2329 s2 = j;
2330 }
2331 else {
2332 b2 = -j;
2333 s2 = 0;
2334 }
2335 if (k >= 0) {
2336 b5 = 0;
2337 s5 = k;
2338 s2 += k;
2339 }
2340 else {
2341 b2 -= k;
2342 b5 = -k;
2343 s5 = 0;
2344 }
2345 if (mode < 0 || mode > 9)
2346 mode = 0;
2347
2348 try_quick = 1;
2349
2350 if (mode > 5) {
2351 mode -= 4;
2352 try_quick = 0;
2353 }
2354 leftright = 1;
2355 ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
2356 /* silence erroneous "gcc -Wall" warning. */
2357 switch(mode) {
2358 case 0:
2359 case 1:
2360 i = 18;
2361 ndigits = 0;
2362 break;
2363 case 2:
2364 leftright = 0;
2365 /* fall through */
2366 case 4:
2367 if (ndigits <= 0)
2368 ndigits = 1;
2369 ilim = ilim1 = i = ndigits;
2370 break;
2371 case 3:
2372 leftright = 0;
2373 /* fall through */
2374 case 5:
2375 i = ndigits + k + 1;
2376 ilim = i;
2377 ilim1 = i - 1;
2378 if (i <= 0)
2379 i = 1;
2380 }
2381 s0 = rv_alloc(i);
2382 if (s0 == NULL)
2383 goto failed_malloc;
2384 s = s0;
2385
2386
2387 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2388
2389 /* Try to get by with floating-point arithmetic. */
2390
2391 i = 0;
2392 dval(&d2) = dval(&u);
2393 k0 = k;
2394 ilim0 = ilim;
2395 ieps = 2; /* conservative */
2396 if (k > 0) {
2397 ds = tens[k&0xf];
2398 j = k >> 4;
2399 if (j & Bletch) {
2400 /* prevent overflows */
2401 j &= Bletch - 1;
2402 dval(&u) /= bigtens[n_bigtens-1];
2403 ieps++;
2404 }
2405 for(; j; j >>= 1, i++)
2406 if (j & 1) {
2407 ieps++;
2408 ds *= bigtens[i];
2409 }
2410 dval(&u) /= ds;
2411 }
2412 else if ((j1 = -k)) {
2413 dval(&u) *= tens[j1 & 0xf];
2414 for(j = j1 >> 4; j; j >>= 1, i++)
2415 if (j & 1) {
2416 ieps++;
2417 dval(&u) *= bigtens[i];
2418 }
2419 }
2420 if (k_check && dval(&u) < 1. && ilim > 0) {
2421 if (ilim1 <= 0)
2422 goto fast_failed;
2423 ilim = ilim1;
2424 k--;
2425 dval(&u) *= 10.;
2426 ieps++;
2427 }
2428 dval(&eps) = ieps*dval(&u) + 7.;
2429 word0(&eps) -= (P-1)*Exp_msk1;
2430 if (ilim == 0) {
2431 S = mhi = 0;
2432 dval(&u) -= 5.;
2433 if (dval(&u) > dval(&eps))
2434 goto one_digit;
2435 if (dval(&u) < -dval(&eps))
2436 goto no_digits;
2437 goto fast_failed;
2438 }
2439 if (leftright) {
2440 /* Use Steele & White method of only
2441 * generating digits needed.
2442 */
2443 dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
2444 for(i = 0;;) {
2445 L = (Long)dval(&u);
2446 dval(&u) -= L;
2447 *s++ = '0' + (int)L;
2448 if (dval(&u) < dval(&eps))
2449 goto ret1;
2450 if (1. - dval(&u) < dval(&eps))
2451 goto bump_up;
2452 if (++i >= ilim)
2453 break;
2454 dval(&eps) *= 10.;
2455 dval(&u) *= 10.;
2456 }
2457 }
2458 else {
2459 /* Generate ilim digits, then fix them up. */
2460 dval(&eps) *= tens[ilim-1];
2461 for(i = 1;; i++, dval(&u) *= 10.) {
2462 L = (Long)(dval(&u));
2463 if (!(dval(&u) -= L))
2464 ilim = i;
2465 *s++ = '0' + (int)L;
2466 if (i == ilim) {
2467 if (dval(&u) > 0.5 + dval(&eps))
2468 goto bump_up;
2469 else if (dval(&u) < 0.5 - dval(&eps)) {
2470 while(*--s == '0');
2471 s++;
2472 goto ret1;
2473 }
2474 break;
2475 }
2476 }
2477 }
2478 fast_failed:
2479 s = s0;
2480 dval(&u) = dval(&d2);
2481 k = k0;
2482 ilim = ilim0;
2483 }
2484
2485 /* Do we have a "small" integer? */
2486
2487 if (be >= 0 && k <= Int_max) {
2488 /* Yes. */
2489 ds = tens[k];
2490 if (ndigits < 0 && ilim <= 0) {
2491 S = mhi = 0;
2492 if (ilim < 0 || dval(&u) <= 5*ds)
2493 goto no_digits;
2494 goto one_digit;
2495 }
2496 for(i = 1;; i++, dval(&u) *= 10.) {
2497 L = (Long)(dval(&u) / ds);
2498 dval(&u) -= L*ds;
2499 *s++ = '0' + (int)L;
2500 if (!dval(&u)) {
2501 break;
2502 }
2503 if (i == ilim) {
2504 dval(&u) += dval(&u);
2505 if (dval(&u) > ds || (dval(&u) == ds && L & 1)) {
2506 bump_up:
2507 while(*--s == '9')
2508 if (s == s0) {
2509 k++;
2510 *s = '0';
2511 break;
2512 }
2513 ++*s++;
2514 }
2515 else {
2516 /* Strip trailing zeros. This branch was missing from the
2517 original dtoa.c, leading to surplus trailing zeros in
2518 some cases. See bugs.python.org/issue40780. */
2519 while (s > s0 && s[-1] == '0') {
2520 --s;
2521 }
2522 }
2523 break;
2524 }
2525 }
2526 goto ret1;
2527 }
2528
2529 m2 = b2;
2530 m5 = b5;
2531 if (leftright) {
2532 i =
2533 denorm ? be + (Bias + (P-1) - 1 + 1) :
2534 1 + P - bbits;
2535 b2 += i;
2536 s2 += i;
2537 mhi = i2b(1);
2538 if (mhi == NULL)
2539 goto failed_malloc;
2540 }
2541 if (m2 > 0 && s2 > 0) {
2542 i = m2 < s2 ? m2 : s2;
2543 b2 -= i;
2544 m2 -= i;
2545 s2 -= i;
2546 }
2547 if (b5 > 0) {
2548 if (leftright) {
2549 if (m5 > 0) {
2550 mhi = pow5mult(mhi, m5);
2551 if (mhi == NULL)
2552 goto failed_malloc;
2553 b1 = mult(mhi, b);
2554 Bfree(b);
2555 b = b1;
2556 if (b == NULL)
2557 goto failed_malloc;
2558 }
2559 if ((j = b5 - m5)) {
2560 b = pow5mult(b, j);
2561 if (b == NULL)
2562 goto failed_malloc;
2563 }
2564 }
2565 else {
2566 b = pow5mult(b, b5);
2567 if (b == NULL)
2568 goto failed_malloc;
2569 }
2570 }
2571 S = i2b(1);
2572 if (S == NULL)
2573 goto failed_malloc;
2574 if (s5 > 0) {
2575 S = pow5mult(S, s5);
2576 if (S == NULL)
2577 goto failed_malloc;
2578 }
2579
2580 /* Check for special case that d is a normalized power of 2. */
2581
2582 spec_case = 0;
2583 if ((mode < 2 || leftright)
2584 ) {
2585 if (!word1(&u) && !(word0(&u) & Bndry_mask)
2586 && word0(&u) & (Exp_mask & ~Exp_msk1)
2587 ) {
2588 /* The special case */
2589 b2 += Log2P;
2590 s2 += Log2P;
2591 spec_case = 1;
2592 }
2593 }
2594
2595 /* Arrange for convenient computation of quotients:
2596 * shift left if necessary so divisor has 4 leading 0 bits.
2597 *
2598 * Perhaps we should just compute leading 28 bits of S once
2599 * and for all and pass them and a shift to quorem, so it
2600 * can do shifts and ors to compute the numerator for q.
2601 */
2602 #define iInc 28
2603 i = dshift(S, s2);
2604 b2 += i;
2605 m2 += i;
2606 s2 += i;
2607 if (b2 > 0) {
2608 b = lshift(b, b2);
2609 if (b == NULL)
2610 goto failed_malloc;
2611 }
2612 if (s2 > 0) {
2613 S = lshift(S, s2);
2614 if (S == NULL)
2615 goto failed_malloc;
2616 }
2617 if (k_check) {
2618 if (cmp(b,S) < 0) {
2619 k--;
2620 b = multadd(b, 10, 0); /* we botched the k estimate */
2621 if (b == NULL)
2622 goto failed_malloc;
2623 if (leftright) {
2624 mhi = multadd(mhi, 10, 0);
2625 if (mhi == NULL)
2626 goto failed_malloc;
2627 }
2628 ilim = ilim1;
2629 }
2630 }
2631 if (ilim <= 0 && (mode == 3 || mode == 5)) {
2632 if (ilim < 0) {
2633 /* no digits, fcvt style */
2634 no_digits:
2635 k = -1 - ndigits;
2636 goto ret;
2637 }
2638 else {
2639 S = multadd(S, 5, 0);
2640 if (S == NULL)
2641 goto failed_malloc;
2642 if (cmp(b, S) <= 0)
2643 goto no_digits;
2644 }
2645 one_digit:
2646 *s++ = '1';
2647 k++;
2648 goto ret;
2649 }
2650 if (leftright) {
2651 if (m2 > 0) {
2652 mhi = lshift(mhi, m2);
2653 if (mhi == NULL)
2654 goto failed_malloc;
2655 }
2656
2657 /* Compute mlo -- check for special case
2658 * that d is a normalized power of 2.
2659 */
2660
2661 mlo = mhi;
2662 if (spec_case) {
2663 mhi = Balloc(mhi->k);
2664 if (mhi == NULL)
2665 goto failed_malloc;
2666 Bcopy(mhi, mlo);
2667 mhi = lshift(mhi, Log2P);
2668 if (mhi == NULL)
2669 goto failed_malloc;
2670 }
2671
2672 for(i = 1;;i++) {
2673 dig = quorem(b,S) + '0';
2674 /* Do we yet have the shortest decimal string
2675 * that will round to d?
2676 */
2677 j = cmp(b, mlo);
2678 delta = diff(S, mhi);
2679 if (delta == NULL)
2680 goto failed_malloc;
2681 j1 = delta->sign ? 1 : cmp(b, delta);
2682 Bfree(delta);
2683 if (j1 == 0 && mode != 1 && !(word1(&u) & 1)
2684 ) {
2685 if (dig == '9')
2686 goto round_9_up;
2687 if (j > 0)
2688 dig++;
2689 *s++ = dig;
2690 goto ret;
2691 }
2692 if (j < 0 || (j == 0 && mode != 1
2693 && !(word1(&u) & 1)
2694 )) {
2695 if (!b->x[0] && b->wds <= 1) {
2696 goto accept_dig;
2697 }
2698 if (j1 > 0) {
2699 b = lshift(b, 1);
2700 if (b == NULL)
2701 goto failed_malloc;
2702 j1 = cmp(b, S);
2703 if ((j1 > 0 || (j1 == 0 && dig & 1))
2704 && dig++ == '9')
2705 goto round_9_up;
2706 }
2707 accept_dig:
2708 *s++ = dig;
2709 goto ret;
2710 }
2711 if (j1 > 0) {
2712 if (dig == '9') { /* possible if i == 1 */
2713 round_9_up:
2714 *s++ = '9';
2715 goto roundoff;
2716 }
2717 *s++ = dig + 1;
2718 goto ret;
2719 }
2720 *s++ = dig;
2721 if (i == ilim)
2722 break;
2723 b = multadd(b, 10, 0);
2724 if (b == NULL)
2725 goto failed_malloc;
2726 if (mlo == mhi) {
2727 mlo = mhi = multadd(mhi, 10, 0);
2728 if (mlo == NULL)
2729 goto failed_malloc;
2730 }
2731 else {
2732 mlo = multadd(mlo, 10, 0);
2733 if (mlo == NULL)
2734 goto failed_malloc;
2735 mhi = multadd(mhi, 10, 0);
2736 if (mhi == NULL)
2737 goto failed_malloc;
2738 }
2739 }
2740 }
2741 else
2742 for(i = 1;; i++) {
2743 *s++ = dig = quorem(b,S) + '0';
2744 if (!b->x[0] && b->wds <= 1) {
2745 goto ret;
2746 }
2747 if (i >= ilim)
2748 break;
2749 b = multadd(b, 10, 0);
2750 if (b == NULL)
2751 goto failed_malloc;
2752 }
2753
2754 /* Round off last digit */
2755
2756 b = lshift(b, 1);
2757 if (b == NULL)
2758 goto failed_malloc;
2759 j = cmp(b, S);
2760 if (j > 0 || (j == 0 && dig & 1)) {
2761 roundoff:
2762 while(*--s == '9')
2763 if (s == s0) {
2764 k++;
2765 *s++ = '1';
2766 goto ret;
2767 }
2768 ++*s++;
2769 }
2770 else {
2771 while(*--s == '0');
2772 s++;
2773 }
2774 ret:
2775 Bfree(S);
2776 if (mhi) {
2777 if (mlo && mlo != mhi)
2778 Bfree(mlo);
2779 Bfree(mhi);
2780 }
2781 ret1:
2782 Bfree(b);
2783 *s = 0;
2784 *decpt = k + 1;
2785 if (rve)
2786 *rve = s;
2787 return s0;
2788 failed_malloc:
2789 if (S)
2790 Bfree(S);
2791 if (mlo && mlo != mhi)
2792 Bfree(mlo);
2793 if (mhi)
2794 Bfree(mhi);
2795 if (b)
2796 Bfree(b);
2797 if (s0)
2798 _Py_dg_freedtoa(s0);
2799 return NULL;
2800 }
2801
2802 #endif // _PY_SHORT_FLOAT_REPR == 1
2803
2804 PyStatus
_PyDtoa_Init(PyInterpreterState * interp)2805 _PyDtoa_Init(PyInterpreterState *interp)
2806 {
2807 #if _PY_SHORT_FLOAT_REPR == 1 && !defined(Py_USING_MEMORY_DEBUGGER)
2808 Bigint **p5s = interp->dtoa.p5s;
2809
2810 // 5**4 = 625
2811 Bigint *p5 = i2b(625);
2812 if (p5 == NULL) {
2813 return PyStatus_NoMemory();
2814 }
2815 p5s[0] = p5;
2816
2817 // compute 5**8, 5**16, 5**32, ..., 5**512
2818 for (Py_ssize_t i = 1; i < Bigint_Pow5size; i++) {
2819 p5 = mult(p5, p5);
2820 if (p5 == NULL) {
2821 return PyStatus_NoMemory();
2822 }
2823 p5s[i] = p5;
2824 }
2825
2826 #endif
2827 return PyStatus_Ok();
2828 }
2829
2830 void
_PyDtoa_Fini(PyInterpreterState * interp)2831 _PyDtoa_Fini(PyInterpreterState *interp)
2832 {
2833 #if _PY_SHORT_FLOAT_REPR == 1 && !defined(Py_USING_MEMORY_DEBUGGER)
2834 Bigint **p5s = interp->dtoa.p5s;
2835 for (Py_ssize_t i = 0; i < Bigint_Pow5size; i++) {
2836 Bigint *p5 = p5s[i];
2837 p5s[i] = NULL;
2838 Bfree(p5);
2839 }
2840 #endif
2841 }
2842