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