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, *bb1, *bd, *bd0, *bs, *delta;
1518 size_t ndigits, fraclen;
1519
1520 dval(&rv) = 0.;
1521
1522 /* Start parsing. */
1523 c = *(s = s00);
1524
1525 /* Parse optional sign, if present. */
1526 sign = 0;
1527 switch (c) {
1528 case '-':
1529 sign = 1;
1530 /* no break */
1531 case '+':
1532 c = *++s;
1533 }
1534
1535 /* Skip leading zeros: lz is true iff there were leading zeros. */
1536 s1 = s;
1537 while (c == '0')
1538 c = *++s;
1539 lz = s != s1;
1540
1541 /* Point s0 at the first nonzero digit (if any). fraclen will be the
1542 number of digits between the decimal point and the end of the
1543 digit string. ndigits will be the total number of digits ignoring
1544 leading zeros. */
1545 s0 = s1 = s;
1546 while ('0' <= c && c <= '9')
1547 c = *++s;
1548 ndigits = s - s1;
1549 fraclen = 0;
1550
1551 /* Parse decimal point and following digits. */
1552 if (c == '.') {
1553 c = *++s;
1554 if (!ndigits) {
1555 s1 = s;
1556 while (c == '0')
1557 c = *++s;
1558 lz = lz || s != s1;
1559 fraclen += (s - s1);
1560 s0 = s;
1561 }
1562 s1 = s;
1563 while ('0' <= c && c <= '9')
1564 c = *++s;
1565 ndigits += s - s1;
1566 fraclen += s - s1;
1567 }
1568
1569 /* Now lz is true if and only if there were leading zero digits, and
1570 ndigits gives the total number of digits ignoring leading zeros. A
1571 valid input must have at least one digit. */
1572 if (!ndigits && !lz) {
1573 if (se)
1574 *se = (char *)s00;
1575 goto parse_error;
1576 }
1577
1578 /* Range check ndigits and fraclen to make sure that they, and values
1579 computed with them, can safely fit in an int. */
1580 if (ndigits > MAX_DIGITS || fraclen > MAX_DIGITS) {
1581 if (se)
1582 *se = (char *)s00;
1583 goto parse_error;
1584 }
1585 nd = (int)ndigits;
1586 nd0 = (int)ndigits - (int)fraclen;
1587
1588 /* Parse exponent. */
1589 e = 0;
1590 if (c == 'e' || c == 'E') {
1591 s00 = s;
1592 c = *++s;
1593
1594 /* Exponent sign. */
1595 esign = 0;
1596 switch (c) {
1597 case '-':
1598 esign = 1;
1599 /* no break */
1600 case '+':
1601 c = *++s;
1602 }
1603
1604 /* Skip zeros. lz is true iff there are leading zeros. */
1605 s1 = s;
1606 while (c == '0')
1607 c = *++s;
1608 lz = s != s1;
1609
1610 /* Get absolute value of the exponent. */
1611 s1 = s;
1612 abs_exp = 0;
1613 while ('0' <= c && c <= '9') {
1614 abs_exp = 10*abs_exp + (c - '0');
1615 c = *++s;
1616 }
1617
1618 /* abs_exp will be correct modulo 2**32. But 10**9 < 2**32, so if
1619 there are at most 9 significant exponent digits then overflow is
1620 impossible. */
1621 if (s - s1 > 9 || abs_exp > MAX_ABS_EXP)
1622 e = (int)MAX_ABS_EXP;
1623 else
1624 e = (int)abs_exp;
1625 if (esign)
1626 e = -e;
1627
1628 /* A valid exponent must have at least one digit. */
1629 if (s == s1 && !lz)
1630 s = s00;
1631 }
1632
1633 /* Adjust exponent to take into account position of the point. */
1634 e -= nd - nd0;
1635 if (nd0 <= 0)
1636 nd0 = nd;
1637
1638 /* Finished parsing. Set se to indicate how far we parsed */
1639 if (se)
1640 *se = (char *)s;
1641
1642 /* If all digits were zero, exit with return value +-0.0. Otherwise,
1643 strip trailing zeros: scan back until we hit a nonzero digit. */
1644 if (!nd)
1645 goto ret;
1646 for (i = nd; i > 0; ) {
1647 --i;
1648 if (s0[i < nd0 ? i : i+1] != '0') {
1649 ++i;
1650 break;
1651 }
1652 }
1653 e += nd - i;
1654 nd = i;
1655 if (nd0 > nd)
1656 nd0 = nd;
1657
1658 /* Summary of parsing results. After parsing, and dealing with zero
1659 * inputs, we have values s0, nd0, nd, e, sign, where:
1660 *
1661 * - s0 points to the first significant digit of the input string
1662 *
1663 * - nd is the total number of significant digits (here, and
1664 * below, 'significant digits' means the set of digits of the
1665 * significand of the input that remain after ignoring leading
1666 * and trailing zeros).
1667 *
1668 * - nd0 indicates the position of the decimal point, if present; it
1669 * satisfies 1 <= nd0 <= nd. The nd significant digits are in
1670 * s0[0:nd0] and s0[nd0+1:nd+1] using the usual Python half-open slice
1671 * notation. (If nd0 < nd, then s0[nd0] contains a '.' character; if
1672 * nd0 == nd, then s0[nd0] could be any non-digit character.)
1673 *
1674 * - e is the adjusted exponent: the absolute value of the number
1675 * represented by the original input string is n * 10**e, where
1676 * n is the integer represented by the concatenation of
1677 * s0[0:nd0] and s0[nd0+1:nd+1]
1678 *
1679 * - sign gives the sign of the input: 1 for negative, 0 for positive
1680 *
1681 * - the first and last significant digits are nonzero
1682 */
1683
1684 /* put first DBL_DIG+1 digits into integer y and z.
1685 *
1686 * - y contains the value represented by the first min(9, nd)
1687 * significant digits
1688 *
1689 * - if nd > 9, z contains the value represented by significant digits
1690 * with indices in [9, min(16, nd)). So y * 10**(min(16, nd) - 9) + z
1691 * gives the value represented by the first min(16, nd) sig. digits.
1692 */
1693
1694 bc.e0 = e1 = e;
1695 y = z = 0;
1696 for (i = 0; i < nd; i++) {
1697 if (i < 9)
1698 y = 10*y + s0[i < nd0 ? i : i+1] - '0';
1699 else if (i < DBL_DIG+1)
1700 z = 10*z + s0[i < nd0 ? i : i+1] - '0';
1701 else
1702 break;
1703 }
1704
1705 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1706 dval(&rv) = y;
1707 if (k > 9) {
1708 dval(&rv) = tens[k - 9] * dval(&rv) + z;
1709 }
1710 bd0 = 0;
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 Bfree(bd0);
1881 goto failed_malloc;
1882 }
1883 Bcopy(bd, bd0);
1884 bb = sd2b(&rv, bc.scale, &bbe); /* srv = bb * 2^bbe */
1885 if (bb == NULL) {
1886 Bfree(bd);
1887 Bfree(bd0);
1888 goto failed_malloc;
1889 }
1890 /* Record whether lsb of bb is odd, in case we need this
1891 for the round-to-even step later. */
1892 odd = bb->x[0] & 1;
1893
1894 /* tdv = bd * 10**e; srv = bb * 2**bbe */
1895 bs = i2b(1);
1896 if (bs == NULL) {
1897 Bfree(bb);
1898 Bfree(bd);
1899 Bfree(bd0);
1900 goto failed_malloc;
1901 }
1902
1903 if (e >= 0) {
1904 bb2 = bb5 = 0;
1905 bd2 = bd5 = e;
1906 }
1907 else {
1908 bb2 = bb5 = -e;
1909 bd2 = bd5 = 0;
1910 }
1911 if (bbe >= 0)
1912 bb2 += bbe;
1913 else
1914 bd2 -= bbe;
1915 bs2 = bb2;
1916 bb2++;
1917 bd2++;
1918
1919 /* At this stage bd5 - bb5 == e == bd2 - bb2 + bbe, bb2 - bs2 == 1,
1920 and bs == 1, so:
1921
1922 tdv == bd * 10**e = bd * 2**(bbe - bb2 + bd2) * 5**(bd5 - bb5)
1923 srv == bb * 2**bbe = bb * 2**(bbe - bb2 + bb2)
1924 0.5 ulp(srv) == 2**(bbe-1) = bs * 2**(bbe - bb2 + bs2)
1925
1926 It follows that:
1927
1928 M * tdv = bd * 2**bd2 * 5**bd5
1929 M * srv = bb * 2**bb2 * 5**bb5
1930 M * 0.5 ulp(srv) = bs * 2**bs2 * 5**bb5
1931
1932 for some constant M. (Actually, M == 2**(bb2 - bbe) * 5**bb5, but
1933 this fact is not needed below.)
1934 */
1935
1936 /* Remove factor of 2**i, where i = min(bb2, bd2, bs2). */
1937 i = bb2 < bd2 ? bb2 : bd2;
1938 if (i > bs2)
1939 i = bs2;
1940 if (i > 0) {
1941 bb2 -= i;
1942 bd2 -= i;
1943 bs2 -= i;
1944 }
1945
1946 /* Scale bb, bd, bs by the appropriate powers of 2 and 5. */
1947 if (bb5 > 0) {
1948 bs = pow5mult(bs, bb5);
1949 if (bs == NULL) {
1950 Bfree(bb);
1951 Bfree(bd);
1952 Bfree(bd0);
1953 goto failed_malloc;
1954 }
1955 bb1 = mult(bs, bb);
1956 Bfree(bb);
1957 bb = bb1;
1958 if (bb == NULL) {
1959 Bfree(bs);
1960 Bfree(bd);
1961 Bfree(bd0);
1962 goto failed_malloc;
1963 }
1964 }
1965 if (bb2 > 0) {
1966 bb = lshift(bb, bb2);
1967 if (bb == NULL) {
1968 Bfree(bs);
1969 Bfree(bd);
1970 Bfree(bd0);
1971 goto failed_malloc;
1972 }
1973 }
1974 if (bd5 > 0) {
1975 bd = pow5mult(bd, bd5);
1976 if (bd == NULL) {
1977 Bfree(bb);
1978 Bfree(bs);
1979 Bfree(bd0);
1980 goto failed_malloc;
1981 }
1982 }
1983 if (bd2 > 0) {
1984 bd = lshift(bd, bd2);
1985 if (bd == NULL) {
1986 Bfree(bb);
1987 Bfree(bs);
1988 Bfree(bd0);
1989 goto failed_malloc;
1990 }
1991 }
1992 if (bs2 > 0) {
1993 bs = lshift(bs, bs2);
1994 if (bs == NULL) {
1995 Bfree(bb);
1996 Bfree(bd);
1997 Bfree(bd0);
1998 goto failed_malloc;
1999 }
2000 }
2001
2002 /* Now bd, bb and bs are scaled versions of tdv, srv and 0.5 ulp(srv),
2003 respectively. Compute the difference |tdv - srv|, and compare
2004 with 0.5 ulp(srv). */
2005
2006 delta = diff(bb, bd);
2007 if (delta == NULL) {
2008 Bfree(bb);
2009 Bfree(bs);
2010 Bfree(bd);
2011 Bfree(bd0);
2012 goto failed_malloc;
2013 }
2014 dsign = delta->sign;
2015 delta->sign = 0;
2016 i = cmp(delta, bs);
2017 if (bc.nd > nd && i <= 0) {
2018 if (dsign)
2019 break; /* Must use bigcomp(). */
2020
2021 /* Here rv overestimates the truncated decimal value by at most
2022 0.5 ulp(rv). Hence rv either overestimates the true decimal
2023 value by <= 0.5 ulp(rv), or underestimates it by some small
2024 amount (< 0.1 ulp(rv)); either way, rv is within 0.5 ulps of
2025 the true decimal value, so it's possible to exit.
2026
2027 Exception: if scaled rv is a normal exact power of 2, but not
2028 DBL_MIN, then rv - 0.5 ulp(rv) takes us all the way down to the
2029 next double, so the correctly rounded result is either rv - 0.5
2030 ulp(rv) or rv; in this case, use bigcomp to distinguish. */
2031
2032 if (!word1(&rv) && !(word0(&rv) & Bndry_mask)) {
2033 /* rv can't be 0, since it's an overestimate for some
2034 nonzero value. So rv is a normal power of 2. */
2035 j = (int)(word0(&rv) & Exp_mask) >> Exp_shift;
2036 /* rv / 2^bc.scale = 2^(j - 1023 - bc.scale); use bigcomp if
2037 rv / 2^bc.scale >= 2^-1021. */
2038 if (j - bc.scale >= 2) {
2039 dval(&rv) -= 0.5 * sulp(&rv, &bc);
2040 break; /* Use bigcomp. */
2041 }
2042 }
2043
2044 {
2045 bc.nd = nd;
2046 i = -1; /* Discarded digits make delta smaller. */
2047 }
2048 }
2049
2050 if (i < 0) {
2051 /* Error is less than half an ulp -- check for
2052 * special case of mantissa a power of two.
2053 */
2054 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask
2055 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
2056 ) {
2057 break;
2058 }
2059 if (!delta->x[0] && delta->wds <= 1) {
2060 /* exact result */
2061 break;
2062 }
2063 delta = lshift(delta,Log2P);
2064 if (delta == NULL) {
2065 Bfree(bb);
2066 Bfree(bs);
2067 Bfree(bd);
2068 Bfree(bd0);
2069 goto failed_malloc;
2070 }
2071 if (cmp(delta, bs) > 0)
2072 goto drop_down;
2073 break;
2074 }
2075 if (i == 0) {
2076 /* exactly half-way between */
2077 if (dsign) {
2078 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
2079 && word1(&rv) == (
2080 (bc.scale &&
2081 (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) ?
2082 (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
2083 0xffffffff)) {
2084 /*boundary case -- increment exponent*/
2085 word0(&rv) = (word0(&rv) & Exp_mask)
2086 + Exp_msk1
2087 ;
2088 word1(&rv) = 0;
2089 dsign = 0;
2090 break;
2091 }
2092 }
2093 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
2094 drop_down:
2095 /* boundary case -- decrement exponent */
2096 if (bc.scale) {
2097 L = word0(&rv) & Exp_mask;
2098 if (L <= (2*P+1)*Exp_msk1) {
2099 if (L > (P+2)*Exp_msk1)
2100 /* round even ==> */
2101 /* accept rv */
2102 break;
2103 /* rv = smallest denormal */
2104 if (bc.nd > nd)
2105 break;
2106 goto undfl;
2107 }
2108 }
2109 L = (word0(&rv) & Exp_mask) - Exp_msk1;
2110 word0(&rv) = L | Bndry_mask1;
2111 word1(&rv) = 0xffffffff;
2112 break;
2113 }
2114 if (!odd)
2115 break;
2116 if (dsign)
2117 dval(&rv) += sulp(&rv, &bc);
2118 else {
2119 dval(&rv) -= sulp(&rv, &bc);
2120 if (!dval(&rv)) {
2121 if (bc.nd >nd)
2122 break;
2123 goto undfl;
2124 }
2125 }
2126 dsign = 1 - dsign;
2127 break;
2128 }
2129 if ((aadj = ratio(delta, bs)) <= 2.) {
2130 if (dsign)
2131 aadj = aadj1 = 1.;
2132 else if (word1(&rv) || word0(&rv) & Bndry_mask) {
2133 if (word1(&rv) == Tiny1 && !word0(&rv)) {
2134 if (bc.nd >nd)
2135 break;
2136 goto undfl;
2137 }
2138 aadj = 1.;
2139 aadj1 = -1.;
2140 }
2141 else {
2142 /* special case -- power of FLT_RADIX to be */
2143 /* rounded down... */
2144
2145 if (aadj < 2./FLT_RADIX)
2146 aadj = 1./FLT_RADIX;
2147 else
2148 aadj *= 0.5;
2149 aadj1 = -aadj;
2150 }
2151 }
2152 else {
2153 aadj *= 0.5;
2154 aadj1 = dsign ? aadj : -aadj;
2155 if (Flt_Rounds == 0)
2156 aadj1 += 0.5;
2157 }
2158 y = word0(&rv) & Exp_mask;
2159
2160 /* Check for overflow */
2161
2162 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
2163 dval(&rv0) = dval(&rv);
2164 word0(&rv) -= P*Exp_msk1;
2165 adj.d = aadj1 * ulp(&rv);
2166 dval(&rv) += adj.d;
2167 if ((word0(&rv) & Exp_mask) >=
2168 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
2169 if (word0(&rv0) == Big0 && word1(&rv0) == Big1) {
2170 Bfree(bb);
2171 Bfree(bd);
2172 Bfree(bs);
2173 Bfree(bd0);
2174 Bfree(delta);
2175 goto ovfl;
2176 }
2177 word0(&rv) = Big0;
2178 word1(&rv) = Big1;
2179 goto cont;
2180 }
2181 else
2182 word0(&rv) += P*Exp_msk1;
2183 }
2184 else {
2185 if (bc.scale && y <= 2*P*Exp_msk1) {
2186 if (aadj <= 0x7fffffff) {
2187 if ((z = (ULong)aadj) <= 0)
2188 z = 1;
2189 aadj = z;
2190 aadj1 = dsign ? aadj : -aadj;
2191 }
2192 dval(&aadj2) = aadj1;
2193 word0(&aadj2) += (2*P+1)*Exp_msk1 - y;
2194 aadj1 = dval(&aadj2);
2195 }
2196 adj.d = aadj1 * ulp(&rv);
2197 dval(&rv) += adj.d;
2198 }
2199 z = word0(&rv) & Exp_mask;
2200 if (bc.nd == nd) {
2201 if (!bc.scale)
2202 if (y == z) {
2203 /* Can we stop now? */
2204 L = (Long)aadj;
2205 aadj -= L;
2206 /* The tolerances below are conservative. */
2207 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
2208 if (aadj < .4999999 || aadj > .5000001)
2209 break;
2210 }
2211 else if (aadj < .4999999/FLT_RADIX)
2212 break;
2213 }
2214 }
2215 cont:
2216 Bfree(bb);
2217 Bfree(bd);
2218 Bfree(bs);
2219 Bfree(delta);
2220 }
2221 Bfree(bb);
2222 Bfree(bd);
2223 Bfree(bs);
2224 Bfree(bd0);
2225 Bfree(delta);
2226 if (bc.nd > nd) {
2227 error = bigcomp(&rv, s0, &bc);
2228 if (error)
2229 goto failed_malloc;
2230 }
2231
2232 if (bc.scale) {
2233 word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
2234 word1(&rv0) = 0;
2235 dval(&rv) *= dval(&rv0);
2236 }
2237
2238 ret:
2239 return sign ? -dval(&rv) : dval(&rv);
2240
2241 parse_error:
2242 return 0.0;
2243
2244 failed_malloc:
2245 errno = ENOMEM;
2246 return -1.0;
2247
2248 undfl:
2249 return sign ? -0.0 : 0.0;
2250
2251 ovfl:
2252 errno = ERANGE;
2253 /* Can't trust HUGE_VAL */
2254 word0(&rv) = Exp_mask;
2255 word1(&rv) = 0;
2256 return sign ? -dval(&rv) : dval(&rv);
2257
2258 }
2259
2260 static char *
rv_alloc(int i)2261 rv_alloc(int i)
2262 {
2263 int j, k, *r;
2264
2265 j = sizeof(ULong);
2266 for(k = 0;
2267 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (unsigned)i;
2268 j <<= 1)
2269 k++;
2270 r = (int*)Balloc(k);
2271 if (r == NULL)
2272 return NULL;
2273 *r = k;
2274 return (char *)(r+1);
2275 }
2276
2277 static char *
nrv_alloc(char * s,char ** rve,int n)2278 nrv_alloc(char *s, char **rve, int n)
2279 {
2280 char *rv, *t;
2281
2282 rv = rv_alloc(n);
2283 if (rv == NULL)
2284 return NULL;
2285 t = rv;
2286 while((*t = *s++)) t++;
2287 if (rve)
2288 *rve = t;
2289 return rv;
2290 }
2291
2292 /* freedtoa(s) must be used to free values s returned by dtoa
2293 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
2294 * but for consistency with earlier versions of dtoa, it is optional
2295 * when MULTIPLE_THREADS is not defined.
2296 */
2297
2298 void
_Py_dg_freedtoa(char * s)2299 _Py_dg_freedtoa(char *s)
2300 {
2301 Bigint *b = (Bigint *)((int *)s - 1);
2302 b->maxwds = 1 << (b->k = *(int*)b);
2303 Bfree(b);
2304 }
2305
2306 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
2307 *
2308 * Inspired by "How to Print Floating-Point Numbers Accurately" by
2309 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
2310 *
2311 * Modifications:
2312 * 1. Rather than iterating, we use a simple numeric overestimate
2313 * to determine k = floor(log10(d)). We scale relevant
2314 * quantities using O(log2(k)) rather than O(k) multiplications.
2315 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
2316 * try to generate digits strictly left to right. Instead, we
2317 * compute with fewer bits and propagate the carry if necessary
2318 * when rounding the final digit up. This is often faster.
2319 * 3. Under the assumption that input will be rounded nearest,
2320 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
2321 * That is, we allow equality in stopping tests when the
2322 * round-nearest rule will give the same floating-point value
2323 * as would satisfaction of the stopping test with strict
2324 * inequality.
2325 * 4. We remove common factors of powers of 2 from relevant
2326 * quantities.
2327 * 5. When converting floating-point integers less than 1e16,
2328 * we use floating-point arithmetic rather than resorting
2329 * to multiple-precision integers.
2330 * 6. When asked to produce fewer than 15 digits, we first try
2331 * to get by with floating-point arithmetic; we resort to
2332 * multiple-precision integer arithmetic only if we cannot
2333 * guarantee that the floating-point calculation has given
2334 * the correctly rounded result. For k requested digits and
2335 * "uniformly" distributed input, the probability is
2336 * something like 10^(k-15) that we must resort to the Long
2337 * calculation.
2338 */
2339
2340 /* Additional notes (METD): (1) returns NULL on failure. (2) to avoid memory
2341 leakage, a successful call to _Py_dg_dtoa should always be matched by a
2342 call to _Py_dg_freedtoa. */
2343
2344 char *
_Py_dg_dtoa(double dd,int mode,int ndigits,int * decpt,int * sign,char ** rve)2345 _Py_dg_dtoa(double dd, int mode, int ndigits,
2346 int *decpt, int *sign, char **rve)
2347 {
2348 /* Arguments ndigits, decpt, sign are similar to those
2349 of ecvt and fcvt; trailing zeros are suppressed from
2350 the returned string. If not null, *rve is set to point
2351 to the end of the return value. If d is +-Infinity or NaN,
2352 then *decpt is set to 9999.
2353
2354 mode:
2355 0 ==> shortest string that yields d when read in
2356 and rounded to nearest.
2357 1 ==> like 0, but with Steele & White stopping rule;
2358 e.g. with IEEE P754 arithmetic , mode 0 gives
2359 1e23 whereas mode 1 gives 9.999999999999999e22.
2360 2 ==> max(1,ndigits) significant digits. This gives a
2361 return value similar to that of ecvt, except
2362 that trailing zeros are suppressed.
2363 3 ==> through ndigits past the decimal point. This
2364 gives a return value similar to that from fcvt,
2365 except that trailing zeros are suppressed, and
2366 ndigits can be negative.
2367 4,5 ==> similar to 2 and 3, respectively, but (in
2368 round-nearest mode) with the tests of mode 0 to
2369 possibly return a shorter string that rounds to d.
2370 With IEEE arithmetic and compilation with
2371 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
2372 as modes 2 and 3 when FLT_ROUNDS != 1.
2373 6-9 ==> Debugging modes similar to mode - 4: don't try
2374 fast floating-point estimate (if applicable).
2375
2376 Values of mode other than 0-9 are treated as mode 0.
2377
2378 Sufficient space is allocated to the return value
2379 to hold the suppressed trailing zeros.
2380 */
2381
2382 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
2383 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
2384 spec_case, try_quick;
2385 Long L;
2386 int denorm;
2387 ULong x;
2388 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
2389 U d2, eps, u;
2390 double ds;
2391 char *s, *s0;
2392
2393 /* set pointers to NULL, to silence gcc compiler warnings and make
2394 cleanup easier on error */
2395 mlo = mhi = S = 0;
2396 s0 = 0;
2397
2398 u.d = dd;
2399 if (word0(&u) & Sign_bit) {
2400 /* set sign for everything, including 0's and NaNs */
2401 *sign = 1;
2402 word0(&u) &= ~Sign_bit; /* clear sign bit */
2403 }
2404 else
2405 *sign = 0;
2406
2407 /* quick return for Infinities, NaNs and zeros */
2408 if ((word0(&u) & Exp_mask) == Exp_mask)
2409 {
2410 /* Infinity or NaN */
2411 *decpt = 9999;
2412 if (!word1(&u) && !(word0(&u) & 0xfffff))
2413 return nrv_alloc("Infinity", rve, 8);
2414 return nrv_alloc("NaN", rve, 3);
2415 }
2416 if (!dval(&u)) {
2417 *decpt = 1;
2418 return nrv_alloc("0", rve, 1);
2419 }
2420
2421 /* compute k = floor(log10(d)). The computation may leave k
2422 one too large, but should never leave k too small. */
2423 b = d2b(&u, &be, &bbits);
2424 if (b == NULL)
2425 goto failed_malloc;
2426 if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
2427 dval(&d2) = dval(&u);
2428 word0(&d2) &= Frac_mask1;
2429 word0(&d2) |= Exp_11;
2430
2431 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2432 * log10(x) = log(x) / log(10)
2433 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2434 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2435 *
2436 * This suggests computing an approximation k to log10(d) by
2437 *
2438 * k = (i - Bias)*0.301029995663981
2439 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2440 *
2441 * We want k to be too large rather than too small.
2442 * The error in the first-order Taylor series approximation
2443 * is in our favor, so we just round up the constant enough
2444 * to compensate for any error in the multiplication of
2445 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2446 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2447 * adding 1e-13 to the constant term more than suffices.
2448 * Hence we adjust the constant term to 0.1760912590558.
2449 * (We could get a more accurate k by invoking log10,
2450 * but this is probably not worthwhile.)
2451 */
2452
2453 i -= Bias;
2454 denorm = 0;
2455 }
2456 else {
2457 /* d is denormalized */
2458
2459 i = bbits + be + (Bias + (P-1) - 1);
2460 x = i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32)
2461 : word1(&u) << (32 - i);
2462 dval(&d2) = x;
2463 word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
2464 i -= (Bias + (P-1) - 1) + 1;
2465 denorm = 1;
2466 }
2467 ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 +
2468 i*0.301029995663981;
2469 k = (int)ds;
2470 if (ds < 0. && ds != k)
2471 k--; /* want k = floor(ds) */
2472 k_check = 1;
2473 if (k >= 0 && k <= Ten_pmax) {
2474 if (dval(&u) < tens[k])
2475 k--;
2476 k_check = 0;
2477 }
2478 j = bbits - i - 1;
2479 if (j >= 0) {
2480 b2 = 0;
2481 s2 = j;
2482 }
2483 else {
2484 b2 = -j;
2485 s2 = 0;
2486 }
2487 if (k >= 0) {
2488 b5 = 0;
2489 s5 = k;
2490 s2 += k;
2491 }
2492 else {
2493 b2 -= k;
2494 b5 = -k;
2495 s5 = 0;
2496 }
2497 if (mode < 0 || mode > 9)
2498 mode = 0;
2499
2500 try_quick = 1;
2501
2502 if (mode > 5) {
2503 mode -= 4;
2504 try_quick = 0;
2505 }
2506 leftright = 1;
2507 ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
2508 /* silence erroneous "gcc -Wall" warning. */
2509 switch(mode) {
2510 case 0:
2511 case 1:
2512 i = 18;
2513 ndigits = 0;
2514 break;
2515 case 2:
2516 leftright = 0;
2517 /* no break */
2518 case 4:
2519 if (ndigits <= 0)
2520 ndigits = 1;
2521 ilim = ilim1 = i = ndigits;
2522 break;
2523 case 3:
2524 leftright = 0;
2525 /* no break */
2526 case 5:
2527 i = ndigits + k + 1;
2528 ilim = i;
2529 ilim1 = i - 1;
2530 if (i <= 0)
2531 i = 1;
2532 }
2533 s0 = rv_alloc(i);
2534 if (s0 == NULL)
2535 goto failed_malloc;
2536 s = s0;
2537
2538
2539 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2540
2541 /* Try to get by with floating-point arithmetic. */
2542
2543 i = 0;
2544 dval(&d2) = dval(&u);
2545 k0 = k;
2546 ilim0 = ilim;
2547 ieps = 2; /* conservative */
2548 if (k > 0) {
2549 ds = tens[k&0xf];
2550 j = k >> 4;
2551 if (j & Bletch) {
2552 /* prevent overflows */
2553 j &= Bletch - 1;
2554 dval(&u) /= bigtens[n_bigtens-1];
2555 ieps++;
2556 }
2557 for(; j; j >>= 1, i++)
2558 if (j & 1) {
2559 ieps++;
2560 ds *= bigtens[i];
2561 }
2562 dval(&u) /= ds;
2563 }
2564 else if ((j1 = -k)) {
2565 dval(&u) *= tens[j1 & 0xf];
2566 for(j = j1 >> 4; j; j >>= 1, i++)
2567 if (j & 1) {
2568 ieps++;
2569 dval(&u) *= bigtens[i];
2570 }
2571 }
2572 if (k_check && dval(&u) < 1. && ilim > 0) {
2573 if (ilim1 <= 0)
2574 goto fast_failed;
2575 ilim = ilim1;
2576 k--;
2577 dval(&u) *= 10.;
2578 ieps++;
2579 }
2580 dval(&eps) = ieps*dval(&u) + 7.;
2581 word0(&eps) -= (P-1)*Exp_msk1;
2582 if (ilim == 0) {
2583 S = mhi = 0;
2584 dval(&u) -= 5.;
2585 if (dval(&u) > dval(&eps))
2586 goto one_digit;
2587 if (dval(&u) < -dval(&eps))
2588 goto no_digits;
2589 goto fast_failed;
2590 }
2591 if (leftright) {
2592 /* Use Steele & White method of only
2593 * generating digits needed.
2594 */
2595 dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
2596 for(i = 0;;) {
2597 L = (Long)dval(&u);
2598 dval(&u) -= L;
2599 *s++ = '0' + (int)L;
2600 if (dval(&u) < dval(&eps))
2601 goto ret1;
2602 if (1. - dval(&u) < dval(&eps))
2603 goto bump_up;
2604 if (++i >= ilim)
2605 break;
2606 dval(&eps) *= 10.;
2607 dval(&u) *= 10.;
2608 }
2609 }
2610 else {
2611 /* Generate ilim digits, then fix them up. */
2612 dval(&eps) *= tens[ilim-1];
2613 for(i = 1;; i++, dval(&u) *= 10.) {
2614 L = (Long)(dval(&u));
2615 if (!(dval(&u) -= L))
2616 ilim = i;
2617 *s++ = '0' + (int)L;
2618 if (i == ilim) {
2619 if (dval(&u) > 0.5 + dval(&eps))
2620 goto bump_up;
2621 else if (dval(&u) < 0.5 - dval(&eps)) {
2622 while(*--s == '0');
2623 s++;
2624 goto ret1;
2625 }
2626 break;
2627 }
2628 }
2629 }
2630 fast_failed:
2631 s = s0;
2632 dval(&u) = dval(&d2);
2633 k = k0;
2634 ilim = ilim0;
2635 }
2636
2637 /* Do we have a "small" integer? */
2638
2639 if (be >= 0 && k <= Int_max) {
2640 /* Yes. */
2641 ds = tens[k];
2642 if (ndigits < 0 && ilim <= 0) {
2643 S = mhi = 0;
2644 if (ilim < 0 || dval(&u) <= 5*ds)
2645 goto no_digits;
2646 goto one_digit;
2647 }
2648 for(i = 1;; i++, dval(&u) *= 10.) {
2649 L = (Long)(dval(&u) / ds);
2650 dval(&u) -= L*ds;
2651 *s++ = '0' + (int)L;
2652 if (!dval(&u)) {
2653 break;
2654 }
2655 if (i == ilim) {
2656 dval(&u) += dval(&u);
2657 if (dval(&u) > ds || (dval(&u) == ds && L & 1)) {
2658 bump_up:
2659 while(*--s == '9')
2660 if (s == s0) {
2661 k++;
2662 *s = '0';
2663 break;
2664 }
2665 ++*s++;
2666 }
2667 break;
2668 }
2669 }
2670 goto ret1;
2671 }
2672
2673 m2 = b2;
2674 m5 = b5;
2675 if (leftright) {
2676 i =
2677 denorm ? be + (Bias + (P-1) - 1 + 1) :
2678 1 + P - bbits;
2679 b2 += i;
2680 s2 += i;
2681 mhi = i2b(1);
2682 if (mhi == NULL)
2683 goto failed_malloc;
2684 }
2685 if (m2 > 0 && s2 > 0) {
2686 i = m2 < s2 ? m2 : s2;
2687 b2 -= i;
2688 m2 -= i;
2689 s2 -= i;
2690 }
2691 if (b5 > 0) {
2692 if (leftright) {
2693 if (m5 > 0) {
2694 mhi = pow5mult(mhi, m5);
2695 if (mhi == NULL)
2696 goto failed_malloc;
2697 b1 = mult(mhi, b);
2698 Bfree(b);
2699 b = b1;
2700 if (b == NULL)
2701 goto failed_malloc;
2702 }
2703 if ((j = b5 - m5)) {
2704 b = pow5mult(b, j);
2705 if (b == NULL)
2706 goto failed_malloc;
2707 }
2708 }
2709 else {
2710 b = pow5mult(b, b5);
2711 if (b == NULL)
2712 goto failed_malloc;
2713 }
2714 }
2715 S = i2b(1);
2716 if (S == NULL)
2717 goto failed_malloc;
2718 if (s5 > 0) {
2719 S = pow5mult(S, s5);
2720 if (S == NULL)
2721 goto failed_malloc;
2722 }
2723
2724 /* Check for special case that d is a normalized power of 2. */
2725
2726 spec_case = 0;
2727 if ((mode < 2 || leftright)
2728 ) {
2729 if (!word1(&u) && !(word0(&u) & Bndry_mask)
2730 && word0(&u) & (Exp_mask & ~Exp_msk1)
2731 ) {
2732 /* The special case */
2733 b2 += Log2P;
2734 s2 += Log2P;
2735 spec_case = 1;
2736 }
2737 }
2738
2739 /* Arrange for convenient computation of quotients:
2740 * shift left if necessary so divisor has 4 leading 0 bits.
2741 *
2742 * Perhaps we should just compute leading 28 bits of S once
2743 * and for all and pass them and a shift to quorem, so it
2744 * can do shifts and ors to compute the numerator for q.
2745 */
2746 #define iInc 28
2747 i = dshift(S, s2);
2748 b2 += i;
2749 m2 += i;
2750 s2 += i;
2751 if (b2 > 0) {
2752 b = lshift(b, b2);
2753 if (b == NULL)
2754 goto failed_malloc;
2755 }
2756 if (s2 > 0) {
2757 S = lshift(S, s2);
2758 if (S == NULL)
2759 goto failed_malloc;
2760 }
2761 if (k_check) {
2762 if (cmp(b,S) < 0) {
2763 k--;
2764 b = multadd(b, 10, 0); /* we botched the k estimate */
2765 if (b == NULL)
2766 goto failed_malloc;
2767 if (leftright) {
2768 mhi = multadd(mhi, 10, 0);
2769 if (mhi == NULL)
2770 goto failed_malloc;
2771 }
2772 ilim = ilim1;
2773 }
2774 }
2775 if (ilim <= 0 && (mode == 3 || mode == 5)) {
2776 if (ilim < 0) {
2777 /* no digits, fcvt style */
2778 no_digits:
2779 k = -1 - ndigits;
2780 goto ret;
2781 }
2782 else {
2783 S = multadd(S, 5, 0);
2784 if (S == NULL)
2785 goto failed_malloc;
2786 if (cmp(b, S) <= 0)
2787 goto no_digits;
2788 }
2789 one_digit:
2790 *s++ = '1';
2791 k++;
2792 goto ret;
2793 }
2794 if (leftright) {
2795 if (m2 > 0) {
2796 mhi = lshift(mhi, m2);
2797 if (mhi == NULL)
2798 goto failed_malloc;
2799 }
2800
2801 /* Compute mlo -- check for special case
2802 * that d is a normalized power of 2.
2803 */
2804
2805 mlo = mhi;
2806 if (spec_case) {
2807 mhi = Balloc(mhi->k);
2808 if (mhi == NULL)
2809 goto failed_malloc;
2810 Bcopy(mhi, mlo);
2811 mhi = lshift(mhi, Log2P);
2812 if (mhi == NULL)
2813 goto failed_malloc;
2814 }
2815
2816 for(i = 1;;i++) {
2817 dig = quorem(b,S) + '0';
2818 /* Do we yet have the shortest decimal string
2819 * that will round to d?
2820 */
2821 j = cmp(b, mlo);
2822 delta = diff(S, mhi);
2823 if (delta == NULL)
2824 goto failed_malloc;
2825 j1 = delta->sign ? 1 : cmp(b, delta);
2826 Bfree(delta);
2827 if (j1 == 0 && mode != 1 && !(word1(&u) & 1)
2828 ) {
2829 if (dig == '9')
2830 goto round_9_up;
2831 if (j > 0)
2832 dig++;
2833 *s++ = dig;
2834 goto ret;
2835 }
2836 if (j < 0 || (j == 0 && mode != 1
2837 && !(word1(&u) & 1)
2838 )) {
2839 if (!b->x[0] && b->wds <= 1) {
2840 goto accept_dig;
2841 }
2842 if (j1 > 0) {
2843 b = lshift(b, 1);
2844 if (b == NULL)
2845 goto failed_malloc;
2846 j1 = cmp(b, S);
2847 if ((j1 > 0 || (j1 == 0 && dig & 1))
2848 && dig++ == '9')
2849 goto round_9_up;
2850 }
2851 accept_dig:
2852 *s++ = dig;
2853 goto ret;
2854 }
2855 if (j1 > 0) {
2856 if (dig == '9') { /* possible if i == 1 */
2857 round_9_up:
2858 *s++ = '9';
2859 goto roundoff;
2860 }
2861 *s++ = dig + 1;
2862 goto ret;
2863 }
2864 *s++ = dig;
2865 if (i == ilim)
2866 break;
2867 b = multadd(b, 10, 0);
2868 if (b == NULL)
2869 goto failed_malloc;
2870 if (mlo == mhi) {
2871 mlo = mhi = multadd(mhi, 10, 0);
2872 if (mlo == NULL)
2873 goto failed_malloc;
2874 }
2875 else {
2876 mlo = multadd(mlo, 10, 0);
2877 if (mlo == NULL)
2878 goto failed_malloc;
2879 mhi = multadd(mhi, 10, 0);
2880 if (mhi == NULL)
2881 goto failed_malloc;
2882 }
2883 }
2884 }
2885 else
2886 for(i = 1;; i++) {
2887 *s++ = dig = quorem(b,S) + '0';
2888 if (!b->x[0] && b->wds <= 1) {
2889 goto ret;
2890 }
2891 if (i >= ilim)
2892 break;
2893 b = multadd(b, 10, 0);
2894 if (b == NULL)
2895 goto failed_malloc;
2896 }
2897
2898 /* Round off last digit */
2899
2900 b = lshift(b, 1);
2901 if (b == NULL)
2902 goto failed_malloc;
2903 j = cmp(b, S);
2904 if (j > 0 || (j == 0 && dig & 1)) {
2905 roundoff:
2906 while(*--s == '9')
2907 if (s == s0) {
2908 k++;
2909 *s++ = '1';
2910 goto ret;
2911 }
2912 ++*s++;
2913 }
2914 else {
2915 while(*--s == '0');
2916 s++;
2917 }
2918 ret:
2919 Bfree(S);
2920 if (mhi) {
2921 if (mlo && mlo != mhi)
2922 Bfree(mlo);
2923 Bfree(mhi);
2924 }
2925 ret1:
2926 Bfree(b);
2927 *s = 0;
2928 *decpt = k + 1;
2929 if (rve)
2930 *rve = s;
2931 return s0;
2932 failed_malloc:
2933 if (S)
2934 Bfree(S);
2935 if (mlo && mlo != mhi)
2936 Bfree(mlo);
2937 if (mhi)
2938 Bfree(mhi);
2939 if (b)
2940 Bfree(b);
2941 if (s0)
2942 _Py_dg_freedtoa(s0);
2943 return NULL;
2944 }
2945 #ifdef __cplusplus
2946 }
2947 #endif
2948
2949 #endif /* PY_NO_SHORT_FLOAT_REPR */
2950