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