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